diff --git a/Manifest.toml b/Manifest.toml
index b9a4d2bab4cf15fd9461580ea22700b113281c24..1429f065147f128c8c90b28a46535aea66d5d98f 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.9.0-alpha1"
 manifest_format = "2.0"
-project_hash = "f0da7d5ff50bf02c651315826a4a291410c6b984"
+project_hash = "6a6afac02132d4ea401777f7459614dc2d5cfa37"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra"]
diff --git a/Project.toml b/Project.toml
index 9add9d69ea9cc595f6d292dc10d0940ed6d55998..60bac47b869a5ebf91b4b88001e7d526bae5ee09 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.0.1"
 Agents = "46ada45e-f475-11e8-01d0-f70cc89e6671"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
@@ -16,5 +17,5 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 
 [compat]
+Agents = ">= 5.6"
 julia = ">= 1.9"
-Agents = ">= 5.6"
\ No newline at end of file
diff --git a/src/Persephone.jl b/src/Persephone.jl
index 4c1341bd4817547d49e30a568240c8df1f4196b1..571eb1effd25c5db4aefd378e65dc7d2fd96dc6f 100644
--- a/src/Persephone.jl
+++ b/src/Persephone.jl
@@ -16,10 +16,11 @@ using
     Agents,
     ArgParse,
     Dates,
+    Distributed,
     GeoArrays, #XXX this is a big dependency - can we get rid of it?
     Logging,
     LoggingExtras,
-    #MacroTools, #http://fluxml.ai/MacroTools.jl/stable/utilities/
+    #MacroTools, #may be useful: http://fluxml.ai/MacroTools.jl/stable/utilities/
     Random,
     StableRNGs,
     TOML
diff --git a/src/core/input.jl b/src/core/input.jl
index a049e71f9859f419279355330adca5830c6bb494..d493d2109534251c14d435f5905c30f03c8daf19 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -31,24 +31,6 @@ macro param(domainparam)
     :($(esc(:model)).settings[$(domain*"."*paramname)])
 end
 
-"""
-    flattenTOML(dict)
-
-An internal utility function to convert the two-dimensional dict returned
-by `TOML.parsefile()` into a one-dimensional dict, so that instead of
-writing `settings["domain"]["param"]` one can use `settings["domain.param"]`.
-Can be reversed with `expandTOML()`.
-"""
-function flattenTOML(tomldict)
-    flatdict = Dict{String, Any}()
-    for domain in keys(tomldict)
-        for param in keys(tomldict[domain])
-            flatdict[domain*"."*param] = tomldict[domain][param]
-        end
-    end
-    flatdict
-end
-
 """
     getsettings(configfile, seed=nothing)
 
@@ -95,12 +77,31 @@ function getsettings(configfile::String, seed::Union{Int64,Nothing}=nothing)
         settings["core.outdir"] = outdir
     end
     if settings["core.startdate"] > settings["core.enddate"]
-        Base.error("Enddate is earlier than startdate.")
+        Base.error("Enddate is earlier than startdate.") #TODO replace with exception
     end
+    !isempty(scanparams) && addprocs(settings["core.processors"])
     settings["internal.scanparams"] = scanparams
     settings
 end
 
+"""
+    flattenTOML(dict)
+
+An internal utility function to convert the two-dimensional dict returned
+by `TOML.parsefile()` into a one-dimensional dict, so that instead of
+writing `settings["domain"]["param"]` one can use `settings["domain.param"]`.
+Can be reversed with `expandTOML()`.
+"""
+function flattenTOML(tomldict)
+    flatdict = Dict{String, Any}()
+    for domain in keys(tomldict)
+        for param in keys(tomldict[domain])
+            flatdict[domain*"."*param] = tomldict[domain][param]
+        end
+    end
+    flatdict
+end
+
 """
     parsecommandline()
 
diff --git a/src/core/landscape.jl b/src/core/landscape.jl
index a7f1085eae0b11412a0f49fbedb8a067454d3e02..2d2cbafd598f9b62696e47f78717bf37f03296a0 100644
--- a/src/core/landscape.jl
+++ b/src/core/landscape.jl
@@ -45,6 +45,7 @@ function initlandscape(landcovermap::String, farmfieldsmap::String)
     @debug "Initialising landscape"
     landcover = GeoArrays.read(landcovermap)
     farmfields = GeoArrays.read(farmfieldsmap)
+    #TODO replace error with exception
     (size(landcover) != size(farmfields)) && Base.error("Input map sizes don't match.")
     width, height = size(landcover)[1:2]
     landscape = Matrix{Pixel}(undef, width, height)
diff --git a/src/core/output.jl b/src/core/output.jl
index 09f621e568f26ff356b8ed493b4f52f99294c5f9..6711ff5ef871e6f5213b41980abc3cfd8cd4838c 100644
--- a/src/core/output.jl
+++ b/src/core/output.jl
@@ -23,7 +23,8 @@ function createdatadir(outdir::String, overwrite::Union{Bool,String})
             (answer == "yes") && (overwrite = true)
         end
         if !overwrite
-            Base.error("Output directory exists, will not overwrite. Aborting.")
+            #TODO replace with exception
+            Base.error("Output directory exists, will not overwrite. Aborting.") 
         else
             @warn "Overwriting existing output directory $(outdir)."
         end
@@ -76,6 +77,7 @@ function saveinputfiles(model::AgentBasedModel)
     # Copy the map files to the output folder
     lcmap = @param(core.landcovermap)
     ffmap = @param(core.farmfieldsmap)
+    #TODO replace errors with exceptions
     !(isfile(lcmap)) && Base.error("The map file $(lcmap) doesn't exist.")
     !(isfile(ffmap)) && Base.error("The map file $(ffmap) doesn't exist.")
     cp(lcmap, joinpath(@param(core.outdir), basename(lcmap)), force = true)
@@ -127,7 +129,7 @@ that want to have their output functions called regularly.
 function newdataoutput!(model::AgentBasedModel, filename::String, header::String,
                         outputfunction::Function, frequency::String)
     if !(frequency in ("daily", "monthly", "yearly", "end", "never"))
-        Base.error("Invalid frequency '$frequency' for $filename.")
+        Base.error("Invalid frequency '$frequency' for $filename.") #TODO replace with exception
     end
     ndo = DataOutput(filename, header, outputfunction, frequency)
     append!(model.dataoutputs, [ndo])
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index 0195a9714ffb724e3162a4f469196130819b93af..492d3e7433518a7567bc5acb7d1dba3ef4d7c0d4 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -3,15 +3,59 @@
 ### This file includes the core functions for initialising and running simulations.
 ###
 
-#XXX With the parameter scanning, code execution has become rather difficult to follow.
-# Can I refactor this into two clear, separate paths - one for the normal case (single
-# run) and one for parameter scanning?
+#XXX How can I make the model output during a parallel run clearer?
+
+"""
+    simulate(config=PARAMFILE, seed=nothing)
+
+Initialise one or more model objects and carry out a full simulation experiment,
+optionally specifying a configuration file and a seed for the RNG.
+
+This is the default way to run a Persephone simulation.
+"""
+function simulate(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
+    models = initialise(config, seed)
+    isa(models, Vector) ? 
+        pmap(simulate!, models) :
+        simulate!(models)
+end
+
+"""
+    simulate!(model)
+
+Carry out a complete simulation run using a pre-initialised model object.
+"""
+function simulate!(model::AgentBasedModel)
+    runtime = Dates.value(@param(core.enddate)-@param(core.startdate))+1
+    step!(model, dummystep, stepsimulation!, runtime)
+    finalise!(model)
+end
+
+"""
+    initialise(config=PARAMFILE, seed=nothing)
+
+Initialise the model: read in parameters, create the output data directory,
+and instantiate the AgentBasedModel object(s). Optionally allows specifying the
+configuration file and overriding the `seed` parameter. This returns a single
+model object, unless the config file contains multiple values for one or more
+parameters, in which case it creates a full-factorial simulation experiment
+and returns a vector of model objects.
+"""
+function initialise(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
+    @info "Simulation run started at $(Dates.now())."
+    settings = getsettings(config, seed)
+    scanparams = settings["internal.scanparams"]
+    delete!(settings, "internal.scanparams")
+    isempty(scanparams) ?
+        initmodel(settings) :
+        pmap(initmodel, paramscan(settings, scanparams))
+end
 
 """
     initmodel(settings)
 
 Initialise a model object using a ready-made settings dict. This is
-a helper function for `initialise()` and `initmodelsparallel()`.
+a helper function for `initialise()`.
 """
 function initmodel(settings::Dict{String, Any})
     @debug "Initialising model object."
@@ -38,69 +82,32 @@ function initmodel(settings::Dict{String, Any})
     end
 end
 
-"""
-    initmodelsparallel(settings)
-
-Initialise multiple model objects using ready-made settings dicts. This is
-a helper function for `initialise()`.
-"""
-function initmodelsparallel(settingsdicts::Vector{Dict{String, Any}})
-    #TODO parallelise model initialisation
-    @debug "Beginning to initialise model objects."
-    models = Vector{AgentBasedModel}()
-    for settings in settingsdicts
-        push!(models, initmodel(settings))
-    end
-    models
-end
-
 """
     paramscan(settings)
 
-Initialise a list of model objects, covering all possible parameter combinations
-given by the settings (i.e. a full-factorial experiment). This is a helper function
-for `initialise()`.
+Create a list of settings dicts, covering all possible parameter combinations
+given by the input settings (i.e. a full-factorial experiment). This is a helper
+function for `initialise()`.
 """
 function paramscan(settings::Dict{String,Any}, scanparams::Vector{String})
+    isempty(scanparams) && return [settings]
+    param = pop!(scanparams)
+    combinations = Vector{Dict{String,Any}}()
     # recursively generate a set of settings dicts covering all combinations
-    function generatecombinations(params::Vector{String})
-        (length(params) == 0) && return [settings]
-        param = pop!(params)
-        combinations = Vector{Dict{String,Any}}()
-        for comb in generatecombinations(params)
-            for value in settings[param]
-                newcombination = deepcopy(comb)
-                newcombination[param] = value
-                if comb["core.outdir"] == settings["core.outdir"]
-                    outdir = joinpath(comb["core.outdir"], "$(split(param, ".")[2])_$(value)")
-                else
-                    outdir = "$(comb["core.outdir"])_$(split(param, ".")[2])_$(value)"
-                end
-                newcombination["core.outdir"] = outdir
-                push!(combinations, newcombination)
+    for comb in paramscan(settings, scanparams)
+        for value in settings[param]
+            newcombination = deepcopy(comb)
+            newcombination[param] = value
+            if comb["core.outdir"] == settings["core.outdir"]
+                outdir = joinpath(comb["core.outdir"], "$(split(param, ".")[2])_$(value)")
+            else
+                outdir = "$(comb["core.outdir"])_$(split(param, ".")[2])_$(value)"
             end
+            newcombination["core.outdir"] = outdir
+            push!(combinations, newcombination)
         end
-        combinations
     end
-    generatecombinations(scanparams)
-end
-
-"""
-    initialise(config=PARAMFILE, seed=nothing)
-
-Initialise the model: read in parameters, create the output data directory,
-and instantiate the AgentBasedModel object(s). Optionally allows specifying the
-configuration file and overriding the `seed` parameter. This returns a single
-model object unless the config file contains multiple values for one or more
-parameters, in which case it creates a full-factorial simulation experiment
-and returns a vector of model objects.
-"""
-function initialise(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
-    @info "Simulation run started at $(Dates.now())."
-    settings = getsettings(config, seed)
-    scanparams = settings["internal.scanparams"]
-    delete!(settings, "internal.scanparams")
-    isempty(scanparams) ? initmodel(settings) : initmodelsparallel(settings, scanparameters)
+    combinations
 end
 
 """
@@ -139,34 +146,3 @@ function finalise!(model::AgentBasedModel)
         model
     end
 end
-
-"""
-    simulate!(model)
-
-Carry out a complete simulation run using a pre-initialised model object.
-"""
-function simulate!(model::AgentBasedModel)
-    runtime = Dates.value(@param(core.enddate)-@param(core.startdate))+1
-    step!(model, dummystep, stepsimulation!, runtime)
-    finalise!(model)
-end
-
-"""
-    simulate(config=PARAMFILE, seed=nothing)
-
-Initialise one or more model objects and carry out a full simulation experiment,
-optionally specifying a configuration file and a seed for the RNG.
-
-This is the default way to run a Persephone simulation.
-"""
-function simulate(config::String=PARAMFILE, seed::Union{Int64,Nothing}=nothing)
-    models = initialise(config, seed)
-    if isa(models, Vector)
-        for m in models
-            @info "Executing run $(m.settings["core.outdir"])"
-            simulate!(m) #TODO parallelise
-        end
-    else
-        simulate!(models)
-    end
-end
diff --git a/src/nature/nature.jl b/src/nature/nature.jl
index ffee8c28fce880b27d26a535c431d4573e9813e7..8f4eba6c4e243b058cabaa38466fbf1e6f14ae0f 100644
--- a/src/nature/nature.jl
+++ b/src/nature/nature.jl
@@ -25,6 +25,9 @@ by trait dictionaries passed by them during initialisation.
     age::Int32
 end
 
+#TODO If I write a `getproperty` method for `Animal`, I could get around the ugly
+# `animal.traits[property]` syntax (like Agents.jl does in `model.jl`).
+
 """
     animalid(animal)
 
@@ -155,8 +158,8 @@ variables:
     information).
 
 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`, `@rand`.
+common expressions: `@trait`, `@setphase`, `@respond`, `@kill`, `@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
@@ -215,17 +218,6 @@ macro respond(eventname, body)
     end
 end
 
-"""
-    @here(property)
-
-A utility macro to quickly access a property of the animal's current position
-(i.e. `landcover`, `fieldid`, or  `events` - see the `Pixel` struct).
-This can only be used nested within `@phase`.
-"""
-macro here(property)
-    :($(esc(:model)).landscape[$(esc(:animal)).pos...].$(property))
-end
-
 """
     @kill
 
@@ -369,7 +361,8 @@ 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`.
+This is a utility wrapper that can only be used nested within `@phase` or `@habitat`
+(or in other contexts where the `model` object is available).
 """
 macro rand(args...)
     :($(esc(:rand))($(esc(:model)).rng, $(map(esc, args)...)))
diff --git a/src/parameters.toml b/src/parameters.toml
index e3359e3d1cc5f2523a056efa8adf2945b4cd9928..ea33de601073255b42103bae3fcd2da85ecded13 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -13,6 +13,7 @@ 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", "warn"
+processors = 2 # number of processors to use on parallel runs
 seed = 2 # seed value for the RNG (0 -> random value)
 # dates to start and end the simulation
 startdate = 2022-01-01