diff --git a/CHANGELOG.md b/CHANGELOG.md
index d552f6ac5de137fbc3d2763c4408ba8616d88608..fe6816cbfe23f1b0789b3cd15123a06d4ab18eba 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -21,6 +21,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Add AquaCrop crop model
+
+- Allow multiple crop models to be used in one simulation.  The
+  settings "crop.cropmodel" and "crop.cropdirectory" are now
+  comma-separated lists of crop models and their data directories.
+
+- Simple linear crop height estimation for AquaCrop plants from
+  biomass (AquaCrop does not model plant height)
+
+- Read soil type map, controlled with the setting "world.soiltypesmap"
+
+- Landscape `Pixel`s store their soil type (`enum SoilType`)
+
+- `FarmPlot`s store the most common soil type of their landscape
+  `Pixel`s (AquaCrop needs the soil type as an input parameter).
+
 ### Changed
 
 ### Deprecated
diff --git a/Project.toml b/Project.toml
index 5189975b4a886607b0ef54092816a9b8f3937796..b50507ca1b5c7633cec893eb93bb043dbb1e0a96 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Daniel Vedder <daniel.vedder@idiv.de>"]
 version = "0.6.1"
 
 [deps]
+AquaCrop = "8f16cebd-c0b4-44a3-857f-c101f083256c"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -14,6 +15,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -28,4 +30,6 @@ TiffImages = "731e570b-9d59-4bfa-96dc-6df516fadf69"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+AquaCrop = "0.1.0"
+ImageMagick = "1.2.1"
 julia = ">= 1.10"
diff --git a/data/crops/aquacrop/crop_data.csv b/data/crops/aquacrop/crop_data.csv
new file mode 100644
index 0000000000000000000000000000000000000000..89ed9154613b335c9a68ab5b7011a1f3d8eaa61a
--- /dev/null
+++ b/data/crops/aquacrop/crop_data.csv
@@ -0,0 +1,5 @@
+persefone_name,aquacrop_name,crop_group,min_sow_date,max_sow_date
+winter rape,rapeseed,grain,20 August,25 August
+winter wheat,wheatGDD,grain,15 October,31 October
+maize,maize,grain,15 April,30 April
+winter barley,barleyGDD,grain,15 September,30 September
\ No newline at end of file
diff --git a/src/Persefone.jl b/src/Persefone.jl
index 2b519088adf0df6e06f3679f77756462df875393..fa7be8b2a923660bd10b204c53eb8fb6dadb2fbe 100644
--- a/src/Persefone.jl
+++ b/src/Persefone.jl
@@ -132,6 +132,25 @@ The supertype of all agents in the model (animal species, farmer types, farmplot
 """
 abstract type ModelAgent end
 
+"""
+    AbstractCropType
+
+The abstract supertype of all crop types in the model.  Each crop
+model has to define a type `CropType <: AbstractCropType`.
+"""
+abstract type AbstractCropType end
+
+"""
+    AbstractCropState
+
+The abstract supertype of all crop states in the model.  Each crop
+model has to define a type `CropState <: AbstractCropState`.
+"""
+abstract type AbstractCropState end
+
+function makecropstate end
+setsoiltype!(cropstate::AbstractCropState, soiltype) = nothing
+
 function stepagent! end
 
 ## include all module files (note that the order matters - if file
@@ -147,6 +166,7 @@ include("crop/farmplot.jl")
 include("crop/cropmodels.jl")
 include("crop/almass.jl")
 include("crop/simplecrop.jl")
+include("crop/aquacrop.jl")
 
 include("farm/farm.jl")
 include("farm/farmdata.jl")
diff --git a/src/core/input.jl b/src/core/input.jl
index db2bdef3750c7748cfde8b4b8615bae32e65041f..9437f0e61b136fa5d76459aea41e2aec507db33e 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -20,7 +20,7 @@ const PARAMFILE = joinpath(pkgdir(Persefone), "src/parameters.toml")
 """
 The crop models that can be used in the simulation.
 """
-const AVAILABLE_CROPMODELS = ["almass", "simple"]
+const AVAILABLE_CROPMODELS = ["almass", "simple", "aquacrop"]
 
 """
     @param(domainparam)
@@ -104,21 +104,23 @@ function preprocessparameters(settings::Dict{String,Any}, defaultoutdir::String)
             Base.error("Couldn't find map directory $(settings["world.mapdirectory"]).")
         end
     end
-    if !isdir(settings["crop.cropdirectory"])
-        if isdir(abspath(pkgdir(@__MODULE__), settings["crop.cropdirectory"]))
-            settings["crop.cropdirectory"] = abspath(pkgdir(@__MODULE__), settings["crop.cropdirectory"])
-            @debug "Using package directory to load crop data: $(settings["crop.cropdirectory"])."
-        else
-            Base.error("Couldn't find crop directory $(settings["crop.cropdirectory"]).")
+    # check cropmodel settings
+    cropmodels = string.(split(settings["crop.cropmodel"], ","))
+    for cropmodel in cropmodels
+        if !(cropmodel in AVAILABLE_CROPMODELS)
+            error("cropmodel = \"$cropmodel\", but has to be one of: $AVAILABLE_CROPMODELS")
+        end
+    end
+    cropdirs = string.(split(settings["crop.cropdirectory"], ","))
+    for cropdir in cropdirs
+        if !isdir(cropdir)
+            error("Can't access crop model data dir $cropdir")
         end
     end
     # sanity checks
     if settings["core.startdate"] > settings["core.enddate"]
         Base.error("Enddate is earlier than startdate.")
     end
-    if !(settings["crop.cropmodel"] in AVAILABLE_CROPMODELS)
-        Base.error("crop.cropmodel = \"$(settings["crop.cropmodel"])\", but has to be one of: $AVAILABLE_CROPMODELS")
-    end
     #FIXME enable parallelisation
     # if !isempty(settings["internal.scanparams"]) && (nprocs() < 2)
     #     @warn "To parallelise multiple simulation runs, use `julia -p <n>`."
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index e57f075b9b722ae78c920bd058b4a40d50413523..2aa11cedc69826e9961fea71abf2ecb649ba436c 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -12,7 +12,7 @@ This is the heart of the model - a struct that holds all data and state
 for one simulation run. It is created by [`initialise`](@ref) and passed
 as input to most model functions.
 """
-@kwdef mutable struct AgricultureModel{Tcroptype, Tcropstate} <: SimulationModel
+@kwdef mutable struct AgricultureModel <: SimulationModel
     settings = Dict{String, Any}()
     rng::AbstractRNG
     logger::AbstractLogger
@@ -20,9 +20,9 @@ as input to most model functions.
     date::Date
     landscape = Matrix{Pixel}()
     weather::Weather
-    crops = Dict{String, Tcroptype}()
+    crops = Dict{String, AbstractCropType}()
     farmers = Vector{Farmer}()
-    farmplots = Vector{FarmPlot{Tcropstate}}()
+    farmplots = Vector{FarmPlot}()
     animals = Vector{Union{Animal, Nothing}}()
     migrants = Vector{Pair{Animal, AnnualDate}}()
     events = Vector{FarmEvent}()
@@ -134,14 +134,15 @@ function initmodel(settings::Dict{String, Any})
     with_logger(logger) do
         landscape = initlandscape(settings["world.mapdirectory"],
                                   settings["world.landcovermap"],
-                                  settings["world.farmfieldsmap"])
+                                  settings["world.farmfieldsmap"],
+                                  settings["world.soiltypesmap"])
         weather = initweather(joinpath(settings["world.mapdirectory"],
                                        settings["world.weatherfile"]),
                               settings["core.startdate"],
                               settings["core.enddate"])
-        crops, Tcroptype, Tcropstate = initcropmodel(settings["crop.cropmodel"],
-                                                     settings["crop.cropdirectory"])
-        model = AgricultureModel{Tcroptype, Tcropstate}(;
+        crops = initcropmodels(settings["crop.cropmodel"],
+                               settings["crop.cropdirectory"])
+        model = AgricultureModel(;
             settings,
             rng = StableRNG(settings["core.seed"]),
             logger,
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index 59b94b9f0cf6ace3c5ebd6869dc10aa56e6b3542..de804b38ba31654d80ce9d2ba9402bab2813f6be 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -43,10 +43,13 @@ const GROWTHFILE = "crop_growth_curves.csv"
 using Persefone:
     @rand,
     @u_str,
+    AbstractCropState,
+    AbstractCropType,
     AnnualDate,
     Management,
     cm,
     SimulationModel,
+    SoilType,
     fertiliser,
     maxtemp,
     mintemp,
@@ -60,6 +63,7 @@ import Persefone:
     cropheight,
     cropcover,
     cropyield,
+    makecropstate,
     sow!,
     harvest!,
     isharvestable,
@@ -116,7 +120,7 @@ end
 The type struct for all crops. Currently follows the crop growth model as
 implemented in ALMaSS.
 """
-struct CropType
+struct CropType <: AbstractCropType
     #FIXME this needs thinking about. The sowing and harvest dates belong in the farm model,
     # not here. Also, we need to harmonise crops across the crop growth models.
     name::String
@@ -155,7 +159,7 @@ cropname(ct::CropType) = ct.name
 The state data for an ALMaSS vegetation point calculation.  Usually
 part of a `FarmPlot`.
 """
-@kwdef mutable struct CropState
+@kwdef mutable struct CropState <: AbstractCropState
     # Class in original ALMaSS code:
     #     `VegElement` from `Landscape/Elements.h`, line 601
     croptype::CropType
@@ -200,6 +204,12 @@ function isharvestable(cs::CropState)
     return cs.ddegs >= last(gdd)
 end
 
+function makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType)
+    phase = month(model.date) < 3 ? janfirst : marchfirst
+    return CropState(; croptype, phase)
+end
+
+
 # Constant in original ALMaSS code:
 #     `EL_VEG_START_LAIT` from `Landscape/Elements.cpp`, line 238-239
 const VEG_START_LAIT::Float64 = 1.08
@@ -398,7 +408,7 @@ Parse a CSV file containing the required parameter values for each crop
 (as produced from the original ALMaSS file by `convert_almass_data.py`).
 """
 function readcropparameters(cropdirectory::String)
-    @debug "Reading crop parameters"
+    @debug "Reading crop parameters for ALMaSS crop model"
     # TODO: the second-last column (sowingdensity) uses `String` as
     # type.  This is because the entry for "peas/beans" has a value of
     # "35/80" for the sowingdensity.
@@ -407,7 +417,7 @@ function readcropparameters(cropdirectory::String)
     growthdata = CSV.File(joinpath(cropdirectory, GROWTHFILE), missingstring="NA",
                           types=[Int,String,String,GrowthPhase,String,
                                  Float64,Float64,Float64,Float64])
-    croptypes = Dict{String,CropType}()
+    croptypes = Dict{String,AbstractCropType}()
     for crop in cropdata
         cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
                                               x.crop_name == crop.name)
diff --git a/src/crop/aquacrop.jl b/src/crop/aquacrop.jl
new file mode 100644
index 0000000000000000000000000000000000000000..9a6dc7233b46af0ed0cbf22c8e21ab62300bbed0
--- /dev/null
+++ b/src/crop/aquacrop.jl
@@ -0,0 +1,253 @@
+# Persefone.jl wrapper for AquaCrop.jl
+
+module AquaCropWrapper
+
+const CROPFILE = "crop_data.csv"
+
+import AquaCrop
+using AquaCrop: ton
+import CSV
+using Dates: Date
+using DataFrames: DataFrame
+using Unitful: @u_str
+
+using Persefone:
+    AbstractCropState,
+    AbstractCropType,
+    AnnualDate,
+    cm,
+    daynumber,
+    SoilType,
+    SimulationModel
+
+import Persefone:
+    stepagent!,
+    croptype,
+    cropname,
+    cropheight,
+    cropcover,
+    cropyield,
+    makecropstate,
+    sow!,
+    harvest!,
+    isharvestable,
+    setsoiltype!
+
+using Unitful: @u_str
+
+# TODO: read crop names directly from AquaCrop.jl
+# names extracted by running this command in the AquaCrop.jl src dir:
+#     grep -nir 'crop_type ==' src/ | cut -d '=' -f 3 | sed 's/$/,/' | sort | uniq
+const AQUACROP_CROPNAMES = [
+    "alfalfaGDD",
+    "barley",
+    "barleyGDD",
+    "cotton",
+    "cottonGDD",
+    "drybean",
+    "drybeanGDD",
+    "maize",
+    "maizeGDD",
+    "oat",
+    "paddyrice",
+    "paddyriceGDD",
+    "potato",
+    "potatoGDD",
+    "quinoa",
+    "rapeseed",
+    "sorghum",
+    "sorghumGDD",
+    "soybean",
+    "soybeanGDD",
+    "sugarbeet",
+    "sugarbeetGDD",
+    "sugarcane",
+    "sunflower",
+    "sunflowerGDD",
+    "tef",
+    "tomato",
+    "tomatoGDD",
+    "wheat",
+    "wheatGDD",
+]
+
+@kwdef struct CropType <: AbstractCropType
+    name::String
+    aquacrop_name::String
+    group::String
+    minsowdate::Union{Missing,AnnualDate}
+    maxsowdate::Union{Missing,AnnualDate}
+end
+
+cropname(ct::CropType) = ct.name
+
+"""
+    readcropparameters(cropdirectory)
+
+Parse a CSV file containing some parameters required to map AquaCrop
+crop names to Persefone crop names as well as additional crop data
+needed for Persefone (cropgroup, minsowdate, maxsowdate).
+"""
+function readcropparameters(cropdirectory::String)
+    @debug "Reading extra crop parameters for AquaCrop crop model"
+    croptypes = Dict{String, AbstractCropType}()
+    df = CSV.read(joinpath(cropdirectory, CROPFILE), DataFrame;
+                  missingstring="NA", types=[String, String, String, AnnualDate, AnnualDate])
+    for row in eachrow(df)
+        croptypes[row.persefone_name] =
+            CropType(name = row.persefone_name,
+                     aquacrop_name = row.aquacrop_name,
+                     group = row.crop_group,
+                     minsowdate = row.min_sow_date,
+                     maxsowdate = row.max_sow_date)
+    end
+    return croptypes
+end
+
+function make_aquacrop_cropfield(croptype::CropType,
+                                 model::SimulationModel,
+                                 soiltype::SoilType)
+    # TODO: put this constant somewhere else?
+    MIN_DATASET_LEN = 366
+    # TODO: get this mapping for soil types from a CSV file?
+    soiltype_to_aquacrop = Dict(soiltype => replace(string(soiltype), "_" => " ")
+                                for soiltype in instances(SoilType))
+    aquacrop_soiltype = soiltype_to_aquacrop[soiltype]
+
+    start_date = model.date            # TODO: maybe `today(model)`
+    end_date = model.weather.lastdate  # TODO: maybe `lastdate(model)`
+    start_daynum = daynumber(model.weather, start_date)
+    end_daynum = daynumber(model.weather, end_date)
+    dayrange = start_daynum:end_daynum
+    numdays = length(dayrange)
+    if numdays >= MIN_DATASET_LEN
+        input_Tmin = @view model.weather.mintemp[dayrange]
+        input_Tmax = @view model.weather.maxtemp[dayrange]
+        input_ETo = @view model.weather.evapotranspiration[dayrange]
+        input_Rain = @view model.weather.precipitation[dayrange]
+    else
+        # Append fake data so that there are always more than MIN_DATASET_LEN datapoints
+        input_Tmin = append!(model.weather.mintemp[dayrange], fill(10.0, MIN_DATASET_LEN - numdays))
+        input_Tmax = append!(model.weather.maxtemp[dayrange], fill(16.0, MIN_DATASET_LEN - numdays))
+        input_ETo = append!(model.weather.evapotranspiration[dayrange], fill(1.0, MIN_DATASET_LEN - numdays))
+        input_Rain = append!(model.weather.precipitation[dayrange], fill(1.0, MIN_DATASET_LEN - numdays))
+    end
+
+    # Generate the keyword object for the AquaCrop simulation
+    kwargs = (
+        ## Necessary keywords
+        # runtype
+        runtype = AquaCrop.NoFileRun(),
+        # project input
+        Simulation_DayNr1 = start_date,
+        Simulation_DayNrN = end_date,
+        Crop_Day1 = start_date,    # TODO: originally start_date + Week(1)
+        Crop_DayN = end_date,
+        # soil
+        soil_type = aquacrop_soiltype,
+        # crop
+        crop_type = croptype.aquacrop_name,
+        # Climate
+        InitialClimDate = start_date,
+        ## Optional keyworkds
+        # Climate
+        Tmin = input_Tmin,
+        Tmax = input_Tmax,
+        ETo = input_ETo,
+        Rain = input_Rain,
+        # change soil properties
+        soil_layers = Dict("Thickness" => 5.0)
+    )
+
+    aquacrop_cropfield, all_ok = AquaCrop.start_cropfield(; kwargs...)
+    if ! all_ok.logi
+        @error "Failed calling AquaCrop.start_cropfield()\nAquaCrop error: $(all_ok.msg)"
+    end
+
+    return aquacrop_cropfield
+end
+
+mutable struct CropState <: AbstractCropState
+    croptype::CropType
+    soiltype::SoilType
+    cropstate::AquaCrop.AquaCropField
+
+    CropState(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
+        new(croptype, soiltype, make_aquacrop_cropfield(croptype, model, soiltype))
+end
+
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+function cropheight(cs::CropState)
+    # AquaCrop does not explicitly model crop height, so we have to
+    # estimate crop height from biomass
+    biomass = get_aquacrop_biomass(cs)
+    # TODO: linear model, fixed params
+    a = 2cm
+    b = 6u"cm * ha" / ton
+    return a + b * biomass
+end
+cropcover(cs::CropState) = AquaCrop.canopycover(cs.cropstate)
+cropyield(cs::CropState) = AquaCrop.dryyield(cs.cropstate)  # TODO: there is also freshyield
+function isharvestable(cs::CropState)
+    stage = get_aquacrop_stage(cs)
+    return !ismissing(stage) && stage == 0
+end
+makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
+    CropState(croptype, model, soiltype)
+function setsoiltype!(cs::CropState, soiltype::SoilType)
+    # TODO: this does not affect the actual soil type of the state of
+    # the aquacrop simulation
+    cs.soiltype = soiltype
+    return nothing
+end
+
+function get_aquacrop_stage(cs::CropState)
+    stages = cs.cropstate.dayout.Stage
+    length(stages) > 0 ? last(stages) : missing
+end
+
+function get_aquacrop_biomass(cs::CropState)
+    biomasses = cs.cropstate.dayout.Biomass
+    length(biomasses) > 0 ? last(biomasses) : missing
+end
+
+"""
+    stepagent!(cropstate, model)
+
+Update a crop state by one day.
+"""
+function stepagent!(cs::CropState, model::SimulationModel)
+    AquaCrop.dailyupdate!(cs.cropstate)
+end
+
+"""
+    sow!(cropstate, model, cropname)
+
+Change the cropstate to sow the specified crop.
+"""
+function sow!(cs::CropState, model::SimulationModel, cropname::String)
+    if cropname ∉ keys(model.crops)
+        @error "cropname \"$cropname\" not found"
+    end
+    new_croptype = model.crops[cropname]
+    if typeof(new_croptype) !== CropType
+        error("Cannot change crop model of AquaCrop crop state.")
+    end
+
+    cs.croptype = new_croptype
+    # cs.soiltype stays the way it is
+    cs.cropstate = make_aquacrop_cropfield(cs.croptype, model, soiltype)
+end
+
+"""
+    harvest!(cropstate, model)
+
+Harvest the crop of this cropstate.
+"""
+function harvest!(cs::CropState, model::SimulationModel)
+    # TODO: implement this
+    return 0.0u"g/m^2"
+end
+
+end # module
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 92722a11be435ccfe5d75a1900331bf45a6142cd..5aafe247deebb95d2e67a9b04eb6c1ab9880b270 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -3,27 +3,42 @@
 ### Functions for switchable crop models.
 ###
 
+using StatsBase: mode
+
 """
-    initfields!(cropmodel, cropfile, growthfile)
+    initcropmodels(cropmodels, cropdirectory)
 
-Initialise the crop model.  Returns the crop types available in the
-simulation, as well as the types `Tcroptype` and `Tcropstate`.
+Initialise the crop models given as a comma-delimited string
+(e.g. "almass,aquacrop"). The crop model parameters are read from the
+`cropdirectories`, also a comma-delimited string. Returns the crop
+types available in the simulation.
 """
-function initcropmodel(cropmodel::AbstractString, cropdirectory::AbstractString)
-    if cropmodel == "almass"
-        Tcroptype = ALMaSS.CropType
-        Tcropstate = ALMaSS.CropState
-        crops = ALMaSS.readcropparameters(cropdirectory)
-    elseif cropmodel == "simple"
-        Tcroptype = SimpleCrop.CropType
-        Tcropstate = SimpleCrop.CropState
-        crops_almass = ALMaSS.readcropparameters(cropdirectory)
-        crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
-                     for (name, ct) in crops_almass)
-    else
-        error("initcropmodel: no implementation for crop model '$cropmodel'")
+function initcropmodels(cropmodels::AbstractString,
+                        cropdirectories::AbstractString)
+    models = string.(split(cropmodels, ","))
+    dirs = string.(split(cropdirectories, ","))
+    if length(models) != length(dirs)
+        error("Different number of cropmodels and cropdirectories.\n" *
+              "cropmodels = $cropmodels\ncropdirectories = $cropdirectories")
+    end
+
+    allcrops = Dict{String, AbstractCropType}()
+    for (i, cropmodel) in enumerate(models)
+        cropdirectory = dirs[i]
+        if cropmodel == "almass"
+            crops = ALMaSS.readcropparameters(cropdirectory)
+        elseif cropmodel == "simple"
+            crops_almass = ALMaSS.readcropparameters(cropdirectory)
+            crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
+                         for (name, ct) in crops_almass)
+        elseif cropmodel == "aquacrop"
+            crops = AquaCropWrapper.readcropparameters(cropdirectory)
+        else
+            error("initcropmodel: no implementation for crop model '$cropmodel'")
+        end
+        merge!(allcrops, crops)
     end
-    return crops, Tcroptype, Tcropstate
+    return allcrops
 end
 
 """
@@ -32,6 +47,12 @@ end
 Initialise the farm plots in the simulation model.
 """
 function initfields!(model::SimulationModel)
+    # TODO: we take the soiltype from the first pixel in the farmplot,
+    # then later we fix the soiltypes of the farmplots and cropstates
+    # to be the most common soiltype of a farmplot.  It would be
+    # better to first collect all the data for pixels belonging
+    # together, and then create the farmplots and cropstates with the
+    # correct soiltype directly
     convertid = Dict{Int64,Int64}()
     width, height = size(model.landscape)
     for x in 1:width
@@ -45,36 +66,22 @@ function initfields!(model::SimulationModel)
                 model.landscape[x,y].fieldid = objectid
                 push!(model.farmplots[objectid].pixels, (x,y))
             else
-                cropstate = makecropstate(model)
-                fp = FarmPlot(length(model.farmplots) + 1, [(x, y)], -1, cropstate)
+                soiltype = model.landscape[x,y].soiltype
+                cropstate = makecropstate(model.crops["natural grass"], model, soiltype)
+                fp = FarmPlot(length(model.farmplots) + 1, [(x, y)], -1, soiltype, cropstate)
                 push!(model.farmplots, fp)
                 model.landscape[x,y].fieldid = fp.id
                 convertid[rawid] = fp.id
             end
         end
     end
-    @info "Initialised $(length(model.farmplots)) farm plots."
-end
 
-"""
-    makecropstate(model, cropmodel)
-
-An internal utility function to initialise one instance of the configured crop growth model.
-"""
-function makecropstate(model::SimulationModel)
-    if @param(crop.cropmodel) == "almass"
-        phase = (month(model.date) < 3 ? ALMaSS.janfirst : ALMaSS.marchfirst)
-        cs = ALMaSS.CropState(
-            croptype = model.crops["natural grass"],
-            phase = phase,
-        )
-    elseif @param(crop.cropmodel) == "simple"
-        cs = SimpleCrop.CropState(
-            model.crops["natural grass"],
-            0.0m
-        )
-    else
-        Base.error("Unhandled crop model '$(@param(crop.cropmodel))' in makecropstate().")
+    # Adjust farmplot soil type to most common soil type of its
+    # pixels, and set the cropstate soil type to this as well.
+    for fp in model.farmplots
+        soiltype = mode(map(p -> model.landscape[p[1], p[2]].soiltype, fp.pixels))
+        setsoiltype!(fp, soiltype)
     end
-    return cs
+
+    @info "Initialised $(length(model.farmplots)) farm plots."
 end
diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl
index c44d6ad85d2c56f37784ae7df9c9cab89f3d361d..b39e5bde5012334da22bcb85f9359170f83b8ad0 100644
--- a/src/crop/farmplot.jl
+++ b/src/crop/farmplot.jl
@@ -8,26 +8,33 @@
 
 A struct representing a single field, on which a crop can be grown.
 """
-mutable struct FarmPlot{T} <: ModelAgent
+mutable struct FarmPlot <: ModelAgent
     const id::Int64
     pixels::Vector{Tuple{Int64, Int64}}
     farmer::Int64
-    cropstate::T
+    soiltype::SoilType
+    cropstate::AbstractCropState
 end
 
-croptype(f::FarmPlot{T}) where {T} = croptype(f.cropstate)
-cropname(f::FarmPlot{T}) where {T} = cropname(croptype(f))
-cropheight(f::FarmPlot{T}) where {T} = cropheight(f.cropstate)
-cropcover(f::FarmPlot{T}) where {T} = cropcover(f.cropstate)
-cropyield(f::FarmPlot{T}) where {T} = cropyield(f.cropstate)
-isharvestable(f::FarmPlot{T}) where {T} = isharvestable(f.cropstate)
+croptype(f::FarmPlot) = croptype(f.cropstate)
+cropname(f::FarmPlot) = cropname(croptype(f))
+cropheight(f::FarmPlot) = cropheight(f.cropstate)
+cropcover(f::FarmPlot) = cropcover(f.cropstate)
+cropyield(f::FarmPlot) = cropyield(f.cropstate)
+isharvestable(f::FarmPlot) = isharvestable(f.cropstate)
+
+function setsoiltype!(f::FarmPlot, soiltype::SoilType)
+    f.soiltype = soiltype
+    setsoiltype!(f.cropstate, soiltype)
+    return nothing
+end
 
 """
     stepagent!(farmplot, model)
 
 Update a farm plot by one day.
 """
-function stepagent!(farmplot::FarmPlot{T}, model::SimulationModel) where T
+function stepagent!(farmplot::FarmPlot, model::SimulationModel)
     stepagent!(farmplot.cropstate, model)
 end
 
@@ -37,8 +44,13 @@ end
 Sow the specified crop on the farmplot.
 """
 function sow!(farmplot::FarmPlot, model::SimulationModel, cropname::String)
+    new_croptype = model.crops[cropname]
     createevent!(model, farmplot.pixels, sowing)
-    sow!(farmplot.cropstate, model, cropname)
+    if typeof(new_croptype) === typeof(croptype(farmplot))
+        sow!(farmplot.cropstate, model, cropname)
+    else
+        farmplot.cropstate = makecropstate(new_croptype, model, farmplot.soiltype)
+    end
     @debug "Farmer $(farmplot.farmer) sowed $(cropname) on farmplot $(farmplot.id)."
 end
 
@@ -47,7 +59,7 @@ end
 
 Harvest the crop of this farmplot.
 """
-function harvest!(farmplot::FarmPlot{T}, model::SimulationModel) where T
+function harvest!(farmplot::FarmPlot, model::SimulationModel)
     createevent!(model, farmplot.pixels, harvesting)
     harvest!(farmplot.cropstate, model)  # TODO: multiply with area to return units of `g`
     @debug "Farmer $(farmplot.farmer) harvested $(cropname(farmplot)) from farmplot $(farmplot.id)."
diff --git a/src/crop/simplecrop.jl b/src/crop/simplecrop.jl
index bc4f11dc1c9cd4b52d18274ff626af7e73cddabd..2fd9c29cafb0144db39f29b043130a8c18ea8a1a 100644
--- a/src/crop/simplecrop.jl
+++ b/src/crop/simplecrop.jl
@@ -1,10 +1,12 @@
 module SimpleCrop
 
 using Persefone:
+    AbstractCropState,
+    AbstractCropType,
     AnnualDate,
-    FarmPlot,
     Length,
     cm,
+    SoilType,
     SimulationModel
 
 import Persefone:
@@ -14,6 +16,7 @@ import Persefone:
     cropheight,
     cropcover,
     cropyield,
+    makecropstate,
     sow!,
     harvest!,
     isharvestable
@@ -21,7 +24,7 @@ import Persefone:
 using Unitful: @u_str
 
 # TODO: alternatively just use ALMaSS.CropType ?
-struct CropType
+struct CropType <: AbstractCropType
     name::String
     group::String
     minsowdate::Union{Missing,AnnualDate}
@@ -30,7 +33,7 @@ end
 
 cropname(ct::CropType) = ct.name
 
-mutable struct CropState
+mutable struct CropState <: AbstractCropState
     croptype::CropType
     height::Length{Float64}
 end
@@ -41,6 +44,9 @@ cropheight(cs::CropState) = cs.height
 cropcover(cs::CropState) = 0.0
 cropyield(cs::CropState) = 0.0  # TODO: units?
 isharvestable(cs::CropState) = true
+makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
+    CropState(croptype, 0.0cm)
+
 
 """
     stepagent!(cropstate, model)
diff --git a/src/parameters.toml b/src/parameters.toml
index b743ca5882c71730f481206826029e56320d45c5..1b7d4201b3460daf34f6b674cda7bee31fb0da93 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -23,10 +23,11 @@ enddate = 2023-12-31 # last day of the simulation
 #enddate = 2022-01-02 # last day of the simulation (test value)
 
 [world]
-mapdirectory = "data/regions/jena" # the directory in which all geographic data are stored
+mapdirectory = "data/regions/jena/" # the directory in which all geographic data are stored
 mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # name of the landcover map in the map directory
 farmfieldsmap = "fields.tif" # name of the field geometry map in the map directory
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 weatherfile = "weather.csv" # name of the weather data file in the map directory
 	
 [farm]
@@ -45,5 +46,5 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
 insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting insect growth
 	
 [crop]
-cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass" # the directory storing all data files for the selected crop model
+cropmodel = "almass,aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropdirectory = "data/crops/almass/,data/crops/aquacrop/" # the directory storing all data files for the selected crop model
diff --git a/src/world/landscape.jl b/src/world/landscape.jl
index bf8757d932493f4773d22055a18755cf66944d6b..cdd6d188a03978557988f5273369d7be7f2dbe65 100644
--- a/src/world/landscape.jl
+++ b/src/world/landscape.jl
@@ -9,9 +9,55 @@ using Printf
 "The land cover classes encoded in the Mundialis Sentinel data."
 @enum LandCover nodata forest grass water builtup soil agriculture
 
+"The soil type of a Pixel or FarmPlot"
+@enum SoilType begin
+    nosoildata = 0   # 0
+    sand             # 1
+    loamy_sand       # 2
+    sandy_loam       # 3
+    loam             # 4
+    silt_loam        # 5
+    silt             # 6
+    sandy_clay_loam  # 7
+    clay_loam        # 8
+    silty_clay_loam  # 9
+    sandy_clay       # 10
+    silty_clay       # 11
+    clay             # 12
+    nosoil           # 13
+end
+
+# TODO: move this into a CSV file and read it from there
+"""Bodenatlas soil type id, corresponding Persefone soil type, and numbers to Persefone SoilType enum and the
+   original Bodenatlas description of the soil type"""
+const df_soiltypes_bodenatlas = DataFrame(
+    NamedTuple{(:bodenatlas_id, :bodenatlas_name, :persefone_soiltype), Tuple{Int, String, SoilType}}.([
+        ( 1, "Abbauflächen",      nosoil),
+        ( 2, "Gewässer",          nosoil),
+        ( 3, "Lehmsande (ls)",    loamy_sand),
+        ( 4, "Lehmschluffe (lu)", silt_loam),
+        ( 5, "Moore",             nosoil),
+        ( 6, "Normallehme (ll)",  loam),
+        ( 7, "Reinsande (ss)",    sand),
+        ( 8, "Sandlehme (sl)",    sandy_loam),
+        ( 9, "Schluffsande (us)", sandy_loam),
+        (10, "Schlufftone (ut)",  silty_clay),
+        (11, "Siedlung",          nosoil),
+        (12, "Tonlehme (tl)",     clay_loam),
+        (13, "Tonschluffe (tu)",  silty_clay_loam),
+        (14, "Watt",              nosoil),
+    ])
+)
+
+"Map a Bodenatlas soil type integer to a Persefone `SoilType` enum"
+const soiltype_bodenatlas_to_persefone =
+    Dict{Int,SoilType}(row.bodenatlas_id => row.persefone_soiltype for row in eachrow(df_soiltypes_bodenatlas))
+
+
 "The types of management event that can be simulated"
 @enum Management tillage sowing fertiliser pesticide harvesting
 
+
 """
     Pixel
 
@@ -21,21 +67,27 @@ in a single object. The model landscape consists of a matrix of pixels.
 """
 mutable struct Pixel
     landcover::LandCover          # land cover class at this position
+    soiltype::SoilType            # soil type class at this position
     fieldid::Union{Missing,Int64} # ID of the farmplot (if any) at this position
     events::Vector{Management}    # management events that have been applied to this pixel
     animals::Vector{Int64}        # IDs of animals currently at this position
     territories::Vector{String}   # IDs of animals that claim this pixel as part of their territory
 end
 
+Pixel(landcover::LandCover, soiltype::SoilType, fieldid::Union{Missing,Int64}) =
+    Pixel(landcover, soiltype, fieldid, Vector{Management}(), Vector{Int64}(), Vector{String}())
 Pixel(landcover::LandCover, fieldid::Union{Missing,Int64}) =
-    Pixel(landcover, fieldid, Vector{Management}(), Vector{Int64}(), Vector{Int64}())
+    Pixel(landcover, nosoildata, fieldid)
 Pixel(landcover::LandCover) = Pixel(landcover, missing)
 
-function Base.show(io::IO, ::MIME"text/plain", mat::T) where T <: AbstractMatrix{Pixel}
-    max_fieldid = maximum(skipmissing(map(x -> getfield(x, :fieldid), mat)); init=0)
+showlandscape(mat::T) where T <: AbstractMatrix{Pixel} = showlandscape(stdout, mat)
+
+function showlandscape(io::IO, mat::T) where T <: AbstractMatrix{Pixel}
     println(io, "Matrix{Pixel}:")
+
     println(io, "  LandCover:")
     nrow, ncol = size(mat)
+    max_fieldid = maximum(skipmissing(map(x -> getfield(x, :fieldid), mat)); init=0)
     for i in axes(mat, 1)
         print(io, "    ")
         for j in axes(mat, 2)
@@ -49,12 +101,33 @@ function Base.show(io::IO, ::MIME"text/plain", mat::T) where T <: AbstractMatrix
         end
         println(io)
     end
-
     # print legend
     legend = join(map(x -> "$(uppercase(first(string(x)))) = $x", instances(Persefone.LandCover)), ", ")
     println(io, "Legend: ", legend)
 
+    soiltype_to_str(s::SoilType) = replace(string(s), r"^soiltype_" => "")
+    soiltype_to_char(s::SoilType) = uppercase(first(soiltype_to_str(s)))
+    println(io)
+    println(io, "  Soil types:")
+    for i in axes(mat, 1)
+        print(io, "    ")
+        for j in axes(mat, 2)
+            charid = soiltype_to_char(mat[i, j].soiltype)
+            fieldid = if ismissing(mat[i, j].fieldid)
+                repeat(" ", ndigits(max_fieldid))
+            else
+                @sprintf("%*s", ndigits(max_fieldid), mat[i, j].fieldid)
+            end
+            print(io, charid, fieldid, " ")
+        end
+        println(io)
+    end
+    # print legend
+    legend = join(map(x -> "$(soiltype_to_char(x)) = $(soiltype_to_str(x))", instances(Persefone.SoilType)), ", ")
+    println(io, "Legend: ", legend)
+
     # print number of unique animals, events, territories
+    println(io)
     nanimals = length(unique(reduce(vcat, map(x -> getfield(x, :animals), mat))))
     nevents = length(unique(reduce(vcat, map(x -> getfield(x, :events), mat))))
     nterritories = length(unique(reduce(vcat, map(x -> getfield(x, :territories), mat))))
@@ -78,22 +151,26 @@ mutable struct FarmEvent
 end
 
 """
-    initlandscape(directory, landcovermap, farmfieldsmap)
+    initlandscape(directory, landcovermap, farmfieldsmap, soiltypesmap)
 
 Initialise the model landscape based on the map files specified in the
 configuration. Returns a matrix of pixels.
 """
-function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String)
+function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String, soiltypesmap::String)
     @debug "Initialising landscape"
     landcovermap = joinpath(directory, landcovermap)
     farmfieldsmap = joinpath(directory, farmfieldsmap)
-    #TODO replace errors with exception
-    !(isfile(landcovermap)) && Base.error("Landcover map $(landcovermap) doesn't exist.")
-    !(isfile(farmfieldsmap)) && Base.error("Field map $(farmfieldsmap) doesn't exist.")
+    soiltypesmap = joinpath(directory, soiltypesmap)
+    !(isfile(landcovermap)) && error("Landcover map $(landcovermap) doesn't exist.")
+    !(isfile(farmfieldsmap)) && error("Field map $(farmfieldsmap) doesn't exist.")
+    !(isfile(soiltypesmap)) && error("Soil type map $(soiltypesmap) doesn't exist.")
     # read in TIFFs, replacing missing values with 0
     landcover = coalesce(GeoArrays.read(landcovermap), 0)
     farmfields = coalesce(GeoArrays.read(farmfieldsmap), 0)
-    (size(landcover) != size(farmfields)) && Base.error("Input map sizes don't match.")
+    soiltypes = coalesce(GeoArrays.read(soiltypesmap), 0)
+    if size(landcover) != size(farmfields) || size(farmfields) != size(soiltypes)
+        error("Input map sizes don't match.")
+    end
     width, height = size(landcover)[1:2]
     landscape = Matrix{Pixel}(undef, width, height)
     for x in 1:width
@@ -101,8 +178,9 @@ function initlandscape(directory::String, landcovermap::String, farmfieldsmap::S
             # convert the numeric landcover value to LandCover, then create the Pixel object
             lcv = LandCover(Int(landcover[x,y][1]/10))
             ff = Int64(farmfields[x,y][1])
+            st = soiltype_bodenatlas_to_persefone[Int(soiltypes[x, y][1])]
             (iszero(ff)) && (ff = missing)
-            landscape[x,y] = Pixel(lcv, ff)
+            landscape[x,y] = Pixel(lcv, st, ff)
         end
     end
     return landscape
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index b412030f5b822b968c616f76524c96174d54d6d7..d7b79fbb6c0ee99cc0afc8a74ca777651eba2b46 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -3,10 +3,13 @@
 ### These are the tests for the crop growth models.
 ###
 
+import Unitful
+
 const TESTPARAM_ALMASS = joinpath(pkgdir(Persefone), "test", "test_parameters_almass.toml")
 const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameters_simplecrop.toml")
+const TESTPARAM_AQUACROP = joinpath(pkgdir(Persefone), "test", "test_parameters_aquacrop.toml")
 
-@testset for configfile in (TESTPARAM_ALMASS, TESTPARAM_SIMPLECROP)
+@testset for configfile in (TESTPARAM_ALMASS, TESTPARAM_SIMPLECROP, TESTPARAM_AQUACROP)
     @testset "Model initialisation" begin
         model = initialise(; configfile)
         @test model isa AgricultureModel
@@ -15,6 +18,7 @@ const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameter
         model = initialise(; configfile)
         stepsimulation!(model)
         @test model isa AgricultureModel
+        # TODO: test stepagent!(), sow!(), harvest!()
     end
 end
 
@@ -27,9 +31,9 @@ end
     mature = false
     events = Ps.Management[]
     force_growth = false
-    fp = FarmPlot(id, pixels, farmer, Ps.ALMaSS.CropState(croptype=ct, phase=Ps.ALMaSS.janfirst))
+    fp = FarmPlot(id, pixels, farmer, Ps.nosoildata, Ps.ALMaSS.CropState(croptype=ct, phase=Ps.ALMaSS.janfirst))
     @test fp isa FarmPlot
-    @test fp isa FarmPlot{Ps.ALMaSS.CropState}
+    @test fp.cropstate isa Ps.ALMaSS.CropState
     @test croptype(fp) isa Ps.ALMaSS.CropType
     @test cropname(fp) isa String
     @test cropheight(fp) isa Length{Float64}
@@ -42,12 +46,33 @@ end
     id = 0
     pixels = [(0, 0)]
     farmer = 0
-    fp = FarmPlot(id, pixels, farmer, Ps.SimpleCrop.CropState(ct, 0.0cm))
+    fp = FarmPlot(id, pixels, farmer, Ps.nosoildata, Ps.SimpleCrop.CropState(ct, 0.0cm))
     @test fp isa FarmPlot
-    @test fp isa FarmPlot{Ps.SimpleCrop.CropState}
+    @test fp.cropstate isa Ps.SimpleCrop.CropState
     @test croptype(fp) isa Ps.SimpleCrop.CropType
     @test cropname(fp) isa String
     @test cropheight(fp) isa Length{Float64}
     @test cropcover(fp) isa Float64
     @test cropyield(fp) isa Float64
 end
+
+@testset "Submodule AquaCrop" begin
+    model = initialise(; configfile=TESTPARAM_AQUACROP)
+
+    PsAC = Ps.AquaCropWrapper
+    ct = PsAC.CropType("olive tree", "alfalfaGDD", "no_crop_group", AnnualDate(4, 1), AnnualDate(5, 1))
+    id = 0
+    pixels = [(0, 0)]
+    farmer = 0
+    soiltype = Ps.sand
+    fp = FarmPlot(id, pixels, farmer, soiltype, PsAC.CropState(ct, model, soiltype))
+    @test fp isa FarmPlot
+    @test fp.cropstate isa PsAC.CropState
+    @test croptype(fp) isa PsAC.CropType
+    @test cropname(fp) isa String
+    # TODO: better testing for cropheight
+    ch = cropheight(fp)
+    @test ismissing(ch) || Unitful.dimension(ch) == Unitful.𝐋
+    @test cropcover(fp) isa Float64
+    @test Unitful.dimension(cropyield(fp)) == Unitful.𝐌 * Unitful.𝐋^-2
+end
diff --git a/test/landscape_tests.jl b/test/landscape_tests.jl
index 7f1fe4266a15ca6651d0057510a0cf2229336c32..3b1044633d801c92fb1de6bc2a7f5519aaad0890 100644
--- a/test/landscape_tests.jl
+++ b/test/landscape_tests.jl
@@ -22,7 +22,7 @@ end
 @testset "Utility functions" begin
     model = inittestmodel()
     iobuf = IOBuffer()
-    show(iobuf, MIME"text/plain"(), model.landscape)
+    Ps.showlandscape(iobuf, model.landscape)
     @test length(take!(iobuf)) > 0
 end
 
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index 0a786cabd2cd81eacd45fd12d9930ee561d57566..84ba55ddcfb720f3dbb438d4bfafda66d7bb796b 100644
--- a/test/nature_tests.jl
+++ b/test/nature_tests.jl
@@ -47,9 +47,9 @@ end) # end eval
 @testset "Habitat macros" begin
     # set up the testing landscape
     model = inittestmodel()
-    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [], [], [])
+    model.landscape[6,6] = Pixel(Ps.agriculture, Ps.sand, 1, [], [], [])
     fp = Ps.FarmPlot(
-        1, [(6,6)], 1,
+        1, [(6,6)], 1, Ps.nosoildata,
         Ps.ALMaSS.CropState(
             croptype = model.crops["winter wheat"],
             phase = Ps.ALMaSS.janfirst,
diff --git a/test/paramscan.toml b/test/paramscan.toml
index 565e6855f876dbc6474137cd22b0353c384485f8..fecac929e70eb2a73d99bc58aef5af8eaf690b17 100644
--- a/test/paramscan.toml
+++ b/test/paramscan.toml
@@ -29,4 +29,4 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
 	
 [crop]
 cropmodel = "simple" # crop growth model to use, "simple" or "almass"
-
+cropdirectory = "../data/crops/almass/"
diff --git a/test/runtests.jl b/test/runtests.jl
index a655c01849151bd01c0fdffea176676d1c11ac1e..3ca43643d9b4ba553548c48a840e6d7136ca35c1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -37,7 +37,8 @@ function inittestmodel(smallmap=true)
     else
         landscape = Ps.initlandscape(TESTSETTINGS["world.mapdirectory"],
                                      TESTSETTINGS["world.landcovermap"],
-                                     TESTSETTINGS["world.farmfieldsmap"])
+                                     TESTSETTINGS["world.farmfieldsmap"],
+                                     TESTSETTINGS["world.soiltypesmap"])
     end
     weather = Ps.initweather(joinpath(TESTSETTINGS["world.mapdirectory"],
                                       TESTSETTINGS["world.weatherfile"]),
@@ -45,7 +46,7 @@ function inittestmodel(smallmap=true)
                              TESTSETTINGS["core.enddate"])
     # TODO: support other crop models besides ALMaSS
     crops = Ps.ALMaSS.readcropparameters(TESTSETTINGS["crop.cropdirectory"])
-    model = AgricultureModel{Ps.ALMaSS.CropType,Ps.ALMaSS.CropState}(;
+    model = AgricultureModel(;
         settings = copy(TESTSETTINGS),
         rng = StableRNG(TESTSETTINGS["core.seed"]),
         logger = global_logger(),
diff --git a/test/test_parameters.toml b/test/test_parameters.toml
index abc3dab62641bf5518f6f65d497d7bf64cd883e2..a9d6b9b7b3261efcb93fc66d47c0c481fe02a2e0 100644
--- a/test/test_parameters.toml
+++ b/test/test_parameters.toml
@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use
@@ -41,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
 	
 [crop]
 cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_almass.toml b/test/test_parameters_almass.toml
index be08e4af053ca7368bd0418b937f624d3ccc6781..c497679a626ec5d4e9073d0015fee7aaf409b4f5 100644
--- a/test/test_parameters_almass.toml
+++ b/test/test_parameters_almass.toml
@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use
@@ -41,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
 
 [crop]
 cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_aquacrop.toml b/test/test_parameters_aquacrop.toml
new file mode 100644
index 0000000000000000000000000000000000000000..13d5935e40c10ac8894a292b1fa5547bf084d951
--- /dev/null
+++ b/test/test_parameters_aquacrop.toml
@@ -0,0 +1,45 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This configuration file is used for the test suite.
+###
+
+#XXX remember that changes here may break tests!
+
+[core]
+configfile = "test_parameters_aquacrop.toml" # location of the configuration file
+outdir = "results_testsuite" # location and name of the output folder
+logoutput = "both"
+overwrite = true # overwrite the output directory? (true/false/"ask")
+csvoutput = true # save collected data in CSV files
+visualise = true # generate result graphs
+storedata = true # keep collected data in memory
+loglevel = "warn" # verbosity level: "debug", "info", "warn"
+processors = 6 # number of processors to use on parallel runs
+seed = 1 # seed value for the RNG (0 -> random value)
+# dates to start and end the simulation
+startdate = 2022-02-01
+enddate = 2022-03-31
+
+[world]
+mapdirectory = "data/regions/jena/" # the directory in which all geographic data are stored
+mapresolution = 10 # map resolution in meters
+landcovermap = "landcover.tif" # location of the landcover map
+farmfieldsmap = "fields.tif" # location of the field geometry map
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
+weatherfile = "weather.csv" # location of the weather data file
+
+[farm]
+farmmodel = "BasicFarmer" # which version of the farm model to use
+setaside = 0.04 # proportion of farm area set aside as fallow
+fieldoutfreq = "daily" # output frequency for crop/field data, daily/monthly/yearly/end/never
+
+[nature]
+targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate - example species
+#targetspecies = ["Skylark"] # list of target species to simulate
+popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
+indoutfreq = "daily" # output frequency individual-level data, daily/monthly/yearly/end/never
+insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
+
+[crop]
+cropmodel = "almass,aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropdirectory = "../data/crops/almass/,../data/crops/aquacrop/"
diff --git a/test/test_parameters_simplecrop.toml b/test/test_parameters_simplecrop.toml
index 1bcf23b8d91bba31d8840cea91d00425e75c8773..f2fd15c830714dd97991234092a2ae8ecd8b33d1 100644
--- a/test/test_parameters_simplecrop.toml
+++ b/test/test_parameters_simplecrop.toml
@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use
@@ -41,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
 
 [crop]
 cropmodel = "simple" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model