diff --git a/run.sh b/run.sh
index 7cd3d64eb8e3095d27d4df1f4ab619602fee180b..fefbafe06748fc5c19b4c5b76d7b47253c1fa35c 100755
--- a/run.sh
+++ b/run.sh
@@ -6,6 +6,6 @@ then
     rm -r temp_results
 fi
 
-./run.jl -o temp_results -s 2
+./run.jl -o temp_results
 
 src/analysis/analyse_nature.R temp_results
diff --git a/src/Persephone.jl b/src/Persephone.jl
index 1369288a797e8b53aca4d9b8f8d64cd029af729a..3d86d2dd5cb70ddcc31abb5ab3c99b4bb1704bf2 100644
--- a/src/Persephone.jl
+++ b/src/Persephone.jl
@@ -21,6 +21,7 @@ using
     LoggingExtras,
     #MacroTools, #http://fluxml.ai/MacroTools.jl/stable/utilities/
     Random,
+    StableRNGs,
     TOML
 
 ## define exported functions and variables
@@ -56,7 +57,7 @@ export
     initialise,
     stepsimulation!,
     createevent!,
-    finalise
+    finalise!
 
 ## include all module files (note that the order matters - if file
 ## b references something from file a, it must be included later)
diff --git a/src/core/input.jl b/src/core/input.jl
index 872486ce5692589eb7dc248d4ccbb9a42c6f243d..814d72d166d2ef8f8b4056559de80d988f312640 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -33,7 +33,7 @@ Precedence: commandline parameters - user config file - default values
 function getsettings(configfile::String, seed::Union{Int64,Nothing}=nothing)
     # read in and merge configurations from the commandline, the default config file
     # and a user-supplied config file
-    defaults = TOML.parsefile(configfile)
+    defaults::Dict{String, Dict{String, Any}} = TOML.parsefile(configfile)
     commandline = parsecommandline()
     if haskey(commandline, "configfile") && isfile(commandline["configfile"])
         configs = TOML.parsefile(commandline["configfile"])
diff --git a/src/core/output.jl b/src/core/output.jl
index dca2801548619f954afde8dec1e0f644206bea76..d0d76261160e72d687c824fbce6204d76b315507 100644
--- a/src/core/output.jl
+++ b/src/core/output.jl
@@ -42,7 +42,7 @@ function setupdatadir(model::AgentBasedModel)
     simulationlogger = TeeLogger(ConsoleLogger(logfile, loglevel),
                                  ConsoleLogger(stdout, loglevel))
     global_logger(simulationlogger)
-    @debug "Setting up output directory $(@param(core.outdir))"
+    @debug "Setting up output directory $(@param(core.outdir))."
     # Export a copy of the current parameter settings to the output folder.
     # This can be used to replicate this exact run in future, and also
     # records the current time and git commit.
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index 05a90d9e389b0687881331257e906209cef499f0..d47c372b05092f85c2c84d09ddafb2b9c3606b55 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -20,7 +20,6 @@ configuration file and overriding the `seed` parameter.
 function initialise(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
     @info "Simulation run started at $(Dates.now())."
     settings = getsettings(config, seed)
-    Random.seed!(settings["core"]["seed"])
     events = Vector{FarmEvent}()
     dataoutputs = Vector{DataOutput}()
     landscape = initlandscape(settings["core"]["landcovermap"], settings["core"]["farmfieldsmap"])
@@ -32,7 +31,7 @@ function initialise(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing
                                   :events=>events)
     @debug "Setting up model."
     model = AgentBasedModel(Union{Farmer,Animal,FarmPlot}, space, properties=properties,
-                            rng=Random.Xoshiro(settings["core"]["seed"]), warn=false)
+                            rng=StableRNG(settings["core"]["seed"]), warn=false)
     setupdatadir(model)
     initfarms!(model)
     initfields!(model)
@@ -61,15 +60,15 @@ function stepsimulation!(model::AgentBasedModel)
 end
 
 """
-    finalise(model)
+    finalise!(model)
 
 Wrap up the simulation. Currently doesn't do anything except print some information.
 """
-function finalise(model::AgentBasedModel)
+function finalise!(model::AgentBasedModel)
     @info "Simulated $(model.date-@param(core.startdate))."
     @info "Simulation run completed at $(Dates.now()),\nwrote output to $(@param(core.outdir))."
     #XXX is there anything to do here?
-    #genocide!(model)
+    model
 end
 
 """
@@ -82,5 +81,5 @@ function simulate(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
     model = initialise(config, seed)
     runtime = Dates.value(@param(core.enddate)-@param(core.startdate))+1
     step!(model, dummystep, stepsimulation!, runtime)
-    finalise(model)
+    finalise!(model)
 end
diff --git a/src/nature/ecologicaldata.jl b/src/nature/ecologicaldata.jl
index 564097a05e00c9c32a6c8fb121bd8247d2132a6e..aa21296c577a4476f92d7657f94daad34260c0a5 100644
--- a/src/nature/ecologicaldata.jl
+++ b/src/nature/ecologicaldata.jl
@@ -21,9 +21,10 @@ end
 """
     savepopulationdata(model)
 
-Print a comma-separated set of lines to `populations.csv`, giving the current date
-and population size for each animal species. May be called never, daily, monthly,
-yearly, or at the end of a simulation, depending on the parameter `nature.popoutfreq`.
+Return a comma-separated set of lines (to be printed to `populations.csv`), giving
+the current date and population size for each animal species. May be called never,
+daily, monthly, yearly, or at the end of a simulation, depending on the parameter
+`nature.popoutfreq`.
 """
 function savepopulationdata(model::AgentBasedModel)
     pops = Dict{String,Int}(s=>0 for s = @param(nature.targetspecies))
@@ -41,10 +42,10 @@ end
 """
     saveindividualdata(model)
 
-Print a comma-separated set of lines to `individuals.csv`, listing all properties
-of all animal individuals in the model. May be called never, daily, monthly, yearly, or
-at the end of a simulation, depending on the parameter `nature.indoutfreq`.
-WARNING: Produces very big files!
+Return a comma-separated set of lines (to be printed to `individuals.csv`), listing
+all properties of all animal individuals in the model. May be called never, daily,
+monthly, yearly, or at the end of a simulation, depending on the parameter
+`nature.indoutfreq`. WARNING: Produces very big files!
 """
 function saveindividualdata(model::AgentBasedModel)
     data = ""
diff --git a/src/nature/nature.jl b/src/nature/nature.jl
index 3606cd1edc5d6ebd469a3e8c41366cebc5bf2100..8010f135f4099dcdd6285f9dc58eb05efba5c23a 100644
--- a/src/nature/nature.jl
+++ b/src/nature/nature.jl
@@ -156,7 +156,7 @@ variables:
 
 Several utility macros can be used within the body of `@phase` as a short-hand for
 common expressions: `@trait`, `@setphase`, `@respond`, `@here`, `@kill`,
-`@reproduce`, `@neighbours`.
+`@reproduce`, `@neighbours`, `@rand`.
 
 Note that the first phase that is defined in a species definition block will be
 the phase that animals are assigned at birth, unless the variable `phase` is
@@ -197,6 +197,7 @@ end
 Switch this animal over to a different phase. This can only be used nested within `@phase`.
 """
 macro setphase(newphase)
+    #XXX make this usable in the top part of a species definition?
     :($(esc(:animal)).traits["phase"] = $(String(newphase)))
 end
 
@@ -364,4 +365,14 @@ macro countanimals(args...)
     :(countanimals($(esc(:pos)), $(esc(:model)); $(map(esc, args)...)))
 end
 
+"""
+    @rand(args...)
+
+Return a random number or element from the sample, using the model RNG.
+This is a utility wrapper that can only be used nested within `@phase` or `@habitat`.
+"""
+macro rand(args...)
+    :(rand($(esc(:model)).rng, $(map(esc, args)...)))
+end
+
 ##TODO add movement macros
diff --git a/src/nature/populations.jl b/src/nature/populations.jl
index dbe6924388356171630e3f0a8ae47e6d807ff46a..a76008898b6fc116ca017883001bc0378cecb62b 100644
--- a/src/nature/populations.jl
+++ b/src/nature/populations.jl
@@ -40,15 +40,15 @@ function initpopulation(habitatdescriptor::Function; phase::Union{String,Nothing
         (!isnothing(phase)) && (species["phase"] = phase)
         width, height = size(model.landscape)
         while n == 0 || n < popsize
-            for x in shuffle!(Vector(1:width))
-                for y in shuffle!(Vector(1:height))
+            for x in shuffle!(model.rng, Vector(1:width))
+                for y in shuffle!(model.rng, Vector(1:height))
                     if habitatdescriptor((x,y), model)
                         if pairs
                             add_agent!((x,y), Animal, model, deepcopy(species), female, 0)
                             add_agent!((x,y), Animal, model, deepcopy(species), male, 0)
                             n += 2
                         else
-                            sex = asexual ? hermaphrodite : rand([male, female])
+                            sex = asexual ? hermaphrodite : rand(model.rng, [male, female])
                             add_agent!((x,y), Animal, model, deepcopy(species), sex, 0)
                             n += 1
                         end
@@ -88,7 +88,7 @@ Produce one or more offspring for the given animal at its current location.
 """
 function reproduce!(animal::Animal, model::AgentBasedModel, n::Int64=1)
     for i in 1:n
-        sex = (animal.sex == hermaphrodite) ? hermaphrodite : rand([male, female])
+        sex = (animal.sex == hermaphrodite) ? hermaphrodite : rand(model.rng, [male, female])
         # We need to generate a fresh species dict here
         species = @eval $(Symbol(animal.traits["name"]))($model)
         add_agent!(animal.pos, Animal, model, species, sex, 0)
@@ -103,7 +103,7 @@ Kill this animal, optionally with a given percentage probability.
 Returns true if the animal dies, false if not.
 """
 function kill!(animal::Animal, model::AgentBasedModel, probability::Float64=1.0, cause::String="")
-    if rand() < probability
+    if rand(model.rng) < probability
         postfix = isempty(cause) ? "." : " from $cause."
         @debug "$(animalid(animal)) has died$(postfix)"
         kill_agent!(animal, model)
diff --git a/src/nature/species/wolpertinger.jl b/src/nature/species/wolpertinger.jl
index a4bc779578e7be0abf97682b67a238ef866aa05b..af1f1b6efc3c9ef04390ee907b5bb4122a7104f8 100644
--- a/src/nature/species/wolpertinger.jl
+++ b/src/nature/species/wolpertinger.jl
@@ -25,11 +25,12 @@ of a deer.
     """
     @phase lifephase begin
         direction = Tuple(rand([-1,1], 2))
-        for i in 1:rand(1:@trait(maxspeed))
+        for i in 1:@rand(1:@trait(maxspeed))
             walk!(animal, direction, model; ifempty=false)
         end
 
-        if rand() < @trait(fecundity)  && @countanimals(species="Wolpertinger") < @trait(crowding)
+        if @rand() < @trait(fecundity) &&
+            @countanimals(species="Wolpertinger") < @trait(crowding)
             @reproduce
         end
 
diff --git a/src/nature/species/wyvern.jl b/src/nature/species/wyvern.jl
index dc6be0b625f49784a395ba875f6fe0243d59488d..f5a8fde4e8c21b98abc814efa3a7c0c0781b6717 100644
--- a/src/nature/species/wyvern.jl
+++ b/src/nature/species/wyvern.jl
@@ -30,15 +30,15 @@ legs, but that doesn't make it any less dangerous...
             # check if a wolpertinger is in pouncing distance
             if a.traits["name"] == "Wolpertinger"
                 move_agent!(animal, a.pos, model)
-                if rand() < @trait(huntsuccess)
+                if @rand() < @trait(huntsuccess)
                     @debug "$(animalid(animal)) killed $(animalid(a))."
                     kill_agent!(a, model)
                     @goto reproduce
                 end
-            elseif a.traits["name"] == "Wyvern" && rand() < @trait(aggression)
+            elseif a.traits["name"] == "Wyvern" && @rand() < @trait(aggression)
                 # wyverns also fight against each other if they get too close
                 move_agent!(animal, a.pos, model)
-                outcome = rand()
+                outcome = @rand()
                 if outcome < 0.4
                     @debug "$(animalid(animal)) killed $(animalid(a)) in a fight."
                     kill_agent!(a, model)
@@ -49,7 +49,7 @@ legs, but that doesn't make it any less dangerous...
             end
         end
         # check if a wolpertinger is in seeing distance, or walk in a random direction
-        direction = Tuple(rand([-1,1], 2))
+        direction = Tuple(@rand([-1,1], 2))
         for a in @neighbours(@trait(vision))
             if a.traits["name"] == "Wolpertinger"
                 direction = get_direction(animal.pos, a.pos, model)
@@ -61,7 +61,7 @@ legs, but that doesn't make it any less dangerous...
         end
         # reproduce every once in a blue moon
         @label reproduce
-        (rand() < @trait(fecundity)) && @reproduce
+        (@rand() < @trait(fecundity)) && @reproduce
         # hibernate from November to March
         if monthday(model.date) == (11,1)
             @trait(phase) = "winter"
diff --git a/src/parameters.toml b/src/parameters.toml
index d0281848d3e4ef0f142e8ca34a597531937adfe0..9f8f8e0e39bd513abc2200a9e1a368552ae03757 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -13,10 +13,10 @@ farmfieldsmap = "data/fields_jena.tif" # location of the field geometry map
 outdir = "results" # location and name of the output folder
 overwrite = "ask" # overwrite the output directory? (true/false/"ask")
 loglevel = "debug" # verbosity level: "debug", "info", "quiet"
-seed = 0 # seed value for the RNG (0 -> random value)
+seed = 2 # seed value for the RNG (0 -> random value)
 # dates to start and end the simulation
 startdate = 2022-01-01
-#enddate = 2022-01-02
+#enddate = 2022-03-31
 enddate = 2022-12-31
 
 [farm]