move soil layers to problem (#113)

* move soil layers to problem

* remove field

Author: rafaqz

src/Microclimate.jl

@@ -38,7 +38,8 @@ export AbstractApparentHeatCapacity, BonacinaStep, TanhSmoothed, Gaussian, Weste
 export AbstractSoilProperties, CampbelldeVriesSoilProperties
 
 # Soil hydraulics
-export AbstractSoilHydraulicsModel, CampbellSoilHydraulics
+export AbstractSoilHydraulicsModel, CampbellSoilHydraulics, CampbellHydraulicProfile
+export example_campbell_hydraulic_profile
 
 # Soil moisture strategy
 export AbstractSoilMoistureStrategy, PrescribedSoilMoisture, DynamicSoilMoisture

src/simulation.jl

@@ -137,13 +137,13 @@ end
 
 function CommonSolve.init(mp::MicroProblem)
     (; hours, depths, heights,
-       soil_profile, soil_properties_model, soil_hydraulic_model, snow_model,
+       soil_properties_model, soil_hydraulic_model, snow_model,
        vapour_pressure_equation, boundary_layer_model,
        radiation, evaporation_model, soil_energy_model,
        config) = mp.model
     days = mp.days
     longwave_model = radiation.longwave_model
-    (; site, environment_minmax, environment_daily, environment_hourly,
+    (; site, soil_profile, environment_minmax, environment_daily, environment_hourly,
        initial_soil_temperature, initial_soil_moisture,
        initial_snow_depth, initial_snow_temperature, initial_snow_density) = mp.inputs
     n_snow = n_snow_nodes(snow_model)
@@ -349,14 +349,14 @@ end
 function solve_soil!(cache::MicroCache)
     mp = cache.problem
     (; hours, depths, heights,
-       soil_profile, soil_properties_model, soil_hydraulic_model, snow_model,
+       soil_properties_model, soil_hydraulic_model, snow_model,
        vapour_pressure_equation, boundary_layer_model,
        radiation, evaporation_model, soil_energy_model,
        config) = mp.model
     days = mp.days
     longwave_model = radiation.longwave_model
     soil_freezing_model = soil_energy_model.freezing_model
-    (; site, environment_daily, environment_hourly,
+    (; site, soil_profile, environment_daily, environment_hourly,
        initial_soil_temperature, initial_soil_moisture,
        initial_snow_depth, initial_snow_temperature, initial_snow_density) = mp.inputs
     n_snow = n_snow_nodes(snow_model)
@@ -372,7 +372,7 @@ function solve_soil!(cache::MicroCache)
     (; convergence, rainfall_schedule, max_surface_pool) = config
     time_mode = mp.time_mode
     moisture_mode = config.soil_moisture_strategy
-    (; campbell_b_parameter, air_entry_water_potential) = soil_hydraulic_model
+    (; campbell_b_parameter, air_entry_water_potential) = soil_profile.hydraulics
     (; bulk_density, mineral_density) = soil_profile
     init_soil_wetness!(moisture_mode)
 

src/soil_hydraulics/campbell.jl

@@ -1,25 +1,45 @@
 """
-    CampbellSoilHydraulics(; ...)
+    CampbellHydraulicProfile(; air_entry_water_potential,
+                              saturated_hydraulic_conductivity,
+                              campbell_b_parameter,
+                              root_density)
 
-Soil hydraulic parameters for Campbell's (1985) soil water balance model.
-Holds parameters only — the moisture-solver strategy lives on
+Per-depth hydraulic parameters for Campbell's (1985) soil water balance
+model. Lives on `SoilProfile.hydraulics` (per-location structural data),
+not on `CampbellSoilHydraulics` (which now carries only the plant-side
+scalars).
+
+Each field is a vector sized to match `MicroModel.depths`.
+"""
+@kwdef struct CampbellHydraulicProfile{AEWP,SHC,CBP,RDen}
+    air_entry_water_potential::AEWP
+    saturated_hydraulic_conductivity::SHC
+    campbell_b_parameter::CBP
+    root_density::RDen
+end
+
+"""
+    CampbellSoilHydraulics(; root_resistance, stomatal_closure_potential,
+                              leaf_resistance, stomatal_stability_parameter,
+                              root_radius)
+
+Plant-side scalar parameters for Campbell's (1985) soil water balance
+model. Holds parameters only — the moisture-solver strategy lives on
 `MicroConfig.soil_moisture_strategy` (`PrescribedSoilMoisture` /
 `DynamicSoilMoisture`), with solver tuning (`moisture_tolerance`,
 `moisture_max_iterations`, `moisture_timestep`) carried on
 `DynamicSoilMoisture` itself.
 
-`bulk_density` and `mineral_density` are not stored here — they live on
-`MicroModel.soil_profile::SoilProfile` and are passed in alongside this
-model wherever they're needed.
+Per-depth soil hydraulic parameters (`air_entry_water_potential`,
+`saturated_hydraulic_conductivity`, `campbell_b_parameter`, `root_density`)
+live on `SoilProfile.hydraulics::CampbellHydraulicProfile`, alongside
+`bulk_density` and `mineral_density` — all per-location data on one
+struct, all passed in alongside this model wherever needed.
 
 # References
 Campbell, G. S. (1985). Soil Physics with BASIC. Elsevier.
 """
-@kwdef struct CampbellSoilHydraulics{AEWP,SHC,CBP,RDen,RRes,SCP,LRes,SSP,RRad} <: AbstractSoilHydraulicsModel
-    air_entry_water_potential::AEWP
-    saturated_hydraulic_conductivity::SHC
-    campbell_b_parameter::CBP
-    root_density::RDen
+@kwdef struct CampbellSoilHydraulics{RRes,SCP,LRes,SSP,RRad} <: AbstractSoilHydraulicsModel
     root_resistance::RRes
     stomatal_closure_potential::SCP
     leaf_resistance::LRes
@@ -28,23 +48,28 @@ Campbell, G. S. (1985). Soil Physics with BASIC. Elsevier.
 end
 
 # TODO move real defaults to the struct keywords
-function example_soil_hydraulic_model(depths=DEFAULT_DEPTHS;
-    # soil hydraulic parameters
-    air_entry_water_potential = fill(0.7, length(depths))u"J/kg", #air entry potential
-    saturated_hydraulic_conductivity = fill(0.0058, length(depths))u"kg*s/m^3", #saturated conductivity
-    campbell_b_parameter = fill(1.7, length(depths)), #soil 'b' parameter
-    # plant parameters
-    root_density = [0, 0, 8.2, 8.0, 7.8, 7.4, 7.1, 6.4, 5.8, 4.8, 4.0, 1.8, 0.9, 0.6, 0.8, 0.4, 0.4, 0, 0] * 1e4u"m/m^3", # root density at each node (from Campell 1985 Soil Physics with Basic, p. 131)
-    root_resistance = 2.5e+10u"m^3/kg/s", # resistance per unit length of root
+function example_campbell_hydraulic_profile(depths=DEFAULT_DEPTHS;
+    air_entry_water_potential = fill(0.7, length(depths))u"J/kg", # air entry potential
+    saturated_hydraulic_conductivity = fill(0.0058, length(depths))u"kg*s/m^3", # saturated conductivity
+    campbell_b_parameter = fill(1.7, length(depths)), # soil 'b' parameter
+    root_density = [0, 0, 8.2, 8.0, 7.8, 7.4, 7.1, 6.4, 5.8, 4.8, 4.0, 1.8, 0.9, 0.6, 0.8, 0.4, 0.4, 0, 0] * 1e4u"m/m^3", # root density at each node (Campbell 1985, p. 131)
+)
+    CampbellHydraulicProfile(;
+        air_entry_water_potential, saturated_hydraulic_conductivity,
+        campbell_b_parameter, root_density,
+    )
+end
+
+function example_soil_hydraulic_model(;
+    root_resistance = 2.5e+10u"m^3/kg/s",     # resistance per unit length of root
     stomatal_closure_potential = -1500.0u"J/kg", # critical leaf water potential for stomatal closure
-    leaf_resistance = 2.0e6u"m^4/kg/s", # resistance per unit length of leaf
-    stomatal_stability_parameter = 10.0, # stability parameter, -
-    root_radius = 0.001u"m", # root radius, m
+    leaf_resistance = 2.0e6u"m^4/kg/s",       # resistance per unit length of leaf
+    stomatal_stability_parameter = 10.0,      # stability parameter, -
+    root_radius = 0.001u"m",                  # root radius, m
 )
     CampbellSoilHydraulics(;
-        air_entry_water_potential, saturated_hydraulic_conductivity, campbell_b_parameter,
-        root_density, root_resistance, stomatal_closure_potential, leaf_resistance, stomatal_stability_parameter,
-        root_radius,
+        root_resistance, stomatal_closure_potential, leaf_resistance,
+        stomatal_stability_parameter, root_radius,
     )
 end
 
@@ -105,11 +130,12 @@ function infiltration_step!(buffers, soil_hydraulic_model::CampbellSoilHydraulic
     relative_humidity_local = local_relative_humidity
 
     dt = moisture_timestep
-    saturated_conductivity = soil_hydraulic_model.saturated_hydraulic_conductivity
-    campbell_b = soil_hydraulic_model.campbell_b_parameter
+    hydraulic_profile = soil_profile.hydraulics
+    saturated_conductivity = hydraulic_profile.saturated_hydraulic_conductivity
+    campbell_b = hydraulic_profile.campbell_b_parameter
     bulk_density = soil_profile.bulk_density
     mineral_density = soil_profile.mineral_density
-    root_density = soil_hydraulic_model.root_density
+    root_density = hydraulic_profile.root_density
     root_resistance_param = soil_hydraulic_model.root_resistance
     stomatal_closure_potential = soil_hydraulic_model.stomatal_closure_potential
     leaf_resistance = soil_hydraulic_model.leaf_resistance
@@ -133,7 +159,7 @@ function infiltration_step!(buffers, soil_hydraulic_model::CampbellSoilHydraulic
     vapor_diffusivity = 2.4e-5u"m^2/s"
 
     # Convert to negative absolute value; fill positions 1..num_layers, extend boundary
-    air_entry_potential[1:num_layers] .= -abs.(soil_hydraulic_model.air_entry_water_potential)
+    air_entry_potential[1:num_layers] .= -abs.(hydraulic_profile.air_entry_water_potential)
     air_entry_potential[num_layers+1] = air_entry_potential[num_layers]
 
     # Saturation water content, m3/m3; extend boundary

src/types.jl

@@ -1,24 +1,29 @@
 """
-    SoilProfile(; bulk_density, mineral_density)
+    SoilProfile(; bulk_density, mineral_density, hydraulics)
 
-Per-depth soil column profile read by both the soil properties model
-(thermal) and the soil hydraulic model. Sized to match `MicroModel.depths`.
+Per-depth soil column profile read by the soil properties model (thermal)
+and the soil hydraulic model. Sized to match `MicroModel.depths`.
 
 - `bulk_density` — dry soil bulk density (kg/m³) at each depth node
 - `mineral_density` — soil mineral density (kg/m³) at each depth node
+- `hydraulics` — per-depth hydraulic parameters (e.g. `CampbellHydraulicProfile`)
+  matched to the chosen `soil_hydraulic_model`
 """
-@kwdef struct SoilProfile{BD,MD}
+@kwdef struct SoilProfile{BD,MD,H}
     bulk_density::BD
     mineral_density::MD
+    hydraulics::H
 end
 
 function example_soil_profile(depths=DEFAULT_DEPTHS;
     bulk_density = 2.56u"Mg/m^3",
     mineral_density = 2.560u"Mg/m^3",
+    hydraulics = example_campbell_hydraulic_profile(depths),
 )
     SoilProfile(;
         bulk_density = bulk_density isa AbstractVector ? bulk_density : fill(bulk_density, length(depths)),
         mineral_density = mineral_density isa AbstractVector ? mineral_density : fill(mineral_density, length(depths)),
+        hydraulics,
     )
 end
 
@@ -64,7 +69,7 @@ end
 
 """
     MicroModel(; hours, depths, heights,
-                 soil_profile, soil_properties_model, soil_hydraulic_model,
+                 soil_properties_model, soil_hydraulic_model,
                  radiation=RadiationModel(),
                  snow_model=NoSnow(),
                  vapour_pressure_equation=GoffGratch(),
@@ -77,8 +82,6 @@ Constant-across-runs scientific description of the simulation:
 
 - solver geometry (`hours`, `depths`, `heights`) # TODO generalise to sub-hourly
 - physical-process models:
-    - `soil_profile::SoilProfile` — per-depth `bulk_density` and `mineral_density`
-      profiles read by both the soil properties and soil hydraulic models
     - `soil_properties_model::AbstractSoilProperties` (e.g. `CampbelldeVriesSoilProperties(...)`)
     - `soil_hydraulic_model::AbstractSoilHydraulicsModel` (e.g. `CampbellSoilHydraulics(...)`)
     - `radiation::RadiationModel` — bundle of `solar_radiation_model`,
@@ -91,17 +94,12 @@ Constant-across-runs scientific description of the simulation:
       (carries the phase-transition `freezing_model` and ODE solver settings)
 - iteration/data-delivery strategy in `config::MicroConfig`
 
-The days of year to simulate live on `MicroProblem`, not here — they're a
-per-run choice (e.g. monthly mid-month days vs daily 1:365) that can vary
-without changing the model.
-
 Combine with a `MicroInputs` via `MicroProblem(model, inputs; days)` to run.
 """
-@kwdef struct MicroModel{H,Dep,Ht,SPR,SPM,SHM,RAD,SNM,VPE,BLM,EVM,SEM,C}
+@kwdef struct MicroModel{H,Dep,Ht,SPM,SHM,RAD,SNM,VPE,BLM,EVM,SEM,C}
     hours::H = DEFAULT_HOURS # hour of day for solar_radiation
     depths::Dep = DEFAULT_DEPTHS # soil nodes - keep spacing close near the surface
     heights::Ht = [0.01, 2]u"m" # air nodes for temperature, wind speed and humidity profile, last height is reference height for weather data
-    soil_profile::SPR
     soil_properties_model::SPM
     soil_hydraulic_model::SHM
     radiation::RAD = RadiationModel()
@@ -114,14 +112,20 @@ Combine with a `MicroInputs` via `MicroProblem(model, inputs; days)` to run.
 end
 
 """
-    MicroInputs(; site, environment_minmax, environment_daily, environment_hourly=nothing,
+    MicroInputs(; site, soil_profile, environment_minmax, environment_daily,
+                  environment_hourly=nothing,
                   initial_soil_temperature=nothing, initial_soil_moisture,
                   initial_snow_depth=0.0u"cm", initial_snow_temperature=u"K"(0.0u"°C"),
                   initial_snow_density=nothing)
 
-Per-run input data: site, environment forcings, and initial conditions.
+Per-run input data: site, soil column profile, environment forcings, and
+initial conditions.
 
 - `site::Site` — properties of the place
+- `soil_profile::SoilProfile` — per-depth `bulk_density` and `mineral_density`
+  read by both the soil properties and soil hydraulic models. Per-location
+  structural soil-column data; lives here (not on `MicroModel`) because it
+  varies between sites without changing the physics.
 - `environment_minmax`, `environment_daily`, `environment_hourly` — input forcings
 - `initial_*` — initial conditions (`initial_soil_temperature=nothing` falls back to
   the day-mean reference air temperature)
@@ -130,8 +134,9 @@ Per-run input data: site, environment forcings, and initial conditions.
 Pass a fresh `MicroInputs` to `reinit!(cache, inputs)` to re-solve the same
 model on different data without rebuilding the cache.
 """
-@kwdef struct MicroInputs{S,EMM,EH,ED,IST,ISM,ISD,ISNT,ISND}
+@kwdef struct MicroInputs{S,SP<:SoilProfile,EMM,EH,ED,IST,ISM,ISD,ISNT,ISND}
     site::S
+    soil_profile::SP
     environment_minmax::EMM
     environment_hourly::EH = nothing
     environment_daily::ED
@@ -176,7 +181,7 @@ function example_microclimate_problem(;
     site = example_site(),
     soil_profile = example_soil_profile(depths),
     soil_properties_model = example_soil_properties_model(),
-    soil_hydraulic_model = example_soil_hydraulic_model(depths),
+    soil_hydraulic_model = example_soil_hydraulic_model(),
     snow_model = NoSnow(),
     environment_minmax = example_monthly_weather(),
     environment_daily = example_daily_environment(days),
@@ -185,10 +190,11 @@ function example_microclimate_problem(;
     initial_soil_temperature = fill(u"K"(7.741667u"°C"), length(depths)),
     initial_soil_moisture = fill(0.42 * 0.25, length(depths)),
 )
-    model = MicroModel(; hours, depths, heights,
-        soil_profile, soil_properties_model, soil_hydraulic_model, snow_model, config)
+    model = MicroModel(; 
+        hours, depths, heights, soil_properties_model, soil_hydraulic_model, snow_model, config
+    )
     inputs = MicroInputs(;
-        site, environment_minmax, environment_daily, environment_hourly,
+        site, soil_profile, environment_minmax, environment_daily, environment_hourly,
         initial_soil_temperature, initial_soil_moisture,
     )
     MicroProblem(model, inputs; days)

test/example_constructors.jl

@@ -9,7 +9,9 @@ using Test
     @test example_hourly_environment() isa HourlyTimeseries
     @test example_soil_properties_model() isa CampbelldeVriesSoilProperties
     @test example_soil_hydraulic_model() isa CampbellSoilHydraulics
+    @test example_campbell_hydraulic_profile() isa CampbellHydraulicProfile
     @test example_soil_profile() isa SoilProfile
+    @test example_soil_profile().hydraulics isa CampbellHydraulicProfile
     @test example_microclimate_problem() isa MicroProblem
 end
 

test/micro_testrun_daily.jl

@@ -82,14 +82,14 @@ _runmoist = Bool(Int(microinput[:runmoist]))
 soil_profile = SoilProfile(;
     bulk_density = (DataFrame(CSV.File("$testdir/data/init_daily/BD.csv"))[:, 2] * 1.0u"Mg/m^3"),
     mineral_density = (DataFrame(CSV.File("$testdir/data/init_daily/DD.csv"))[:, 2] * 1.0u"Mg/m^3"),
+    hydraulics = CampbellHydraulicProfile(;
+        air_entry_water_potential = (DataFrame(CSV.File("$testdir/data/init_daily/PE.csv"))[:, 2] * 1.0u"J/kg"),
+        saturated_hydraulic_conductivity = (DataFrame(CSV.File("$testdir/data/init_daily/KS.csv"))[:, 2] * 1.0u"kg*s/m^3"),
+        campbell_b_parameter = (DataFrame(CSV.File("$testdir/data/init_daily/BB.csv"))[:, 2] * 1.0),
+        root_density = DataFrame(CSV.File("$testdir/data/init_daily/L.csv"))[:, 2] * u"m/m^3",
+    ),
 )
 soil_hydraulic_model = CampbellSoilHydraulics(;
-    # soil hydraulic parameters
-    air_entry_water_potential = (DataFrame(CSV.File("$testdir/data/init_daily/PE.csv"))[:, 2] * 1.0u"J/kg"),
-    saturated_hydraulic_conductivity = (DataFrame(CSV.File("$testdir/data/init_daily/KS.csv"))[:, 2] * 1.0u"kg*s/m^3"),
-    campbell_b_parameter = (DataFrame(CSV.File("$testdir/data/init_daily/BB.csv"))[:, 2] * 1.0),
-    # plant parameters
-    root_density = DataFrame(CSV.File("$testdir/data/init_daily/L.csv"))[:, 2] * u"m/m^3",
     root_resistance = microinput[:RW] * u"m^3/kg/s",
     stomatal_closure_potential = -microinput[:PC] * u"J/kg",
     leaf_resistance = microinput[:RL] * u"m^4/kg/s",
@@ -150,14 +150,14 @@ model = MicroModel(;
     depths, # soil nodes - keep spacing close near the surface
     heights, # air nodes for temperature, wind speed and humidity profile
     radiation,
-    soil_profile,
     soil_properties_model,
     soil_hydraulic_model,
     boundary_layer_model,
     config,
 )
 inputs = MicroInputs(;
     site,
+    soil_profile,
     environment_minmax,
     environment_daily,
     environment_hourly,

test/micro_testrun_monthly.jl

@@ -102,9 +102,10 @@ _runmoist = Bool(Int(microinput[:runmoist]))
 soil_profile = SoilProfile(;
     bulk_density = fill(bulk_density, length(depths)),
     mineral_density = fill(mineral_density, length(depths)),
+    hydraulics = example_campbell_hydraulic_profile(depths;
+        root_density = fill(0.0, length(depths))u"m/m^3"),
 )
-soil_hydraulic_model = example_soil_hydraulic_model(depths;
-    root_density = fill(0.0, length(depths))u"m/m^3")
+soil_hydraulic_model = example_soil_hydraulic_model()
 radiation = RadiationModel(;
     solar_radiation_model = SolarProblem(;
         diffuse_model = Bool(Int(microinput[:IUV])) ? SolarRadiation.ChandrasekharScattering() : SolarRadiation.NoScattering(),
@@ -134,14 +135,14 @@ model = MicroModel(;
     depths,
     heights, # air nodes for temperature, wind speed and humidity profile
     radiation,
-    soil_profile,
     soil_properties_model,
     soil_hydraulic_model,
     boundary_layer_model,
     config,
 )
 inputs = MicroInputs(;
     site,
+    soil_profile,
     environment_minmax,
     environment_daily,
     environment_hourly,