Model setup of soil models (#1170)

Author: kmdeck

experiments/benchmarks/richards.jl

@@ -80,24 +80,6 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     subsurface_space = domain.space.subsurface
 
     start_date = DateTime(2008)
-    (; ν, hydrology_cm, K_sat, θ_r) =
-        ClimaLand.Soil.soil_vangenuchten_parameters(subsurface_space, FT)
-    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
-    f_over = FT(3.28) # 1/m
-    R_sb = FT(1.484e-4 / 1000) # m/s
-    runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
-        f_over = f_over,
-        f_max = ClimaLand.Soil.topmodel_fmax(surface_space, FT),
-        R_sb = R_sb,
-    )
-    soil_params = ClimaLand.Soil.RichardsParameters(;
-        hydrology_cm = hydrology_cm,
-        ν = ν,
-        K_sat = K_sat,
-        S_s = S_s,
-        θ_r = θ_r,
-    )
-
     era5_ncdata_path =
         ClimaLand.Artifacts.era5_land_forcing_data2008_path(; context)
 
@@ -113,20 +95,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 1000),
     )
-    atmos = ClimaLand.PrescribedPrecipitation{FT, typeof(precip)}(precip)
-    bottom_bc = ClimaLand.Soil.WaterFluxBC((p, t) -> 0.0)
-    bc = (;
-        top = ClimaLand.Soil.RichardsAtmosDrivenFluxBC(atmos, runoff_model),
-        bottom = bottom_bc,
-    )
-    model = ClimaLand.Soil.RichardsModel{FT}(;
-        parameters = soil_params,
-        domain = domain,
-        boundary_conditions = bc,
-        sources = (),
-        lateral_flow = false,
+    forcing = (;
+        atmos = ClimaLand.PrescribedPrecipitation{FT, typeof(precip)}(precip),
     )
 
+    model = ClimaLand.Soil.RichardsModel(FT, domain, forcing)
+
     Y, p, cds = initialize(model)
     z = ClimaCore.Fields.coordinate_field(cds.subsurface).z
     lat = ClimaCore.Fields.coordinate_field(domain.space.subsurface).lat
@@ -138,13 +112,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         α::FT,
         n::FT,
         S_s::FT,
-        fmax,
     ) where {FT}
         m = 1 - 1 / n
         zmin = FT(-50.0)
         zmax = FT(0.0)
 
-        z_∇ = FT(zmin / 5.0 + (zmax - zmin) / 2.5 * (fmax - 0.35) / 0.7)
+        z_∇ = FT(zmin / 5.0 + (zmax - zmin) / 2.5)
         if z > z_∇
             S = FT((FT(1) + (α * (z - z_∇))^n)^(-m))
             ϑ_l = S * (ν - θ_r) + θ_r
@@ -154,10 +127,14 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         return FT(ϑ_l)
     end
 
+    hydrology_cm = model.parameters.hydrology_cm
+    ν = model.parameters.ν
+    θ_r = model.parameters.S_s
+    S_s = model.parameters.θ_r
     # Set initial state values
     vg_α = hydrology_cm.α
     vg_n = hydrology_cm.n
-    Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max)
+    Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s)
     # Create model update functions
     set_initial_cache! = make_set_initial_cache(model)
     exp_tendency! = make_exp_tendency(model)

experiments/long_runs/land_region.jl

@@ -29,7 +29,7 @@ import ClimaUtilities.TimeVaryingInputs:
 import ClimaUtilities.TimeManager: ITime, date
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
-
+using ClimaCore
 using ClimaLand
 using ClimaLand.Snow
 using ClimaLand.Soil

experiments/long_runs/snowy_land.jl

@@ -30,7 +30,7 @@ import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaParams as CP
-
+using ClimaCore
 using ClimaLand
 using ClimaLand.Snow
 using ClimaLand.Soil

experiments/long_runs/soil.jl

@@ -6,13 +6,13 @@
 # land/sea mask are not yet supported by ClimaCore.
 
 # Simulation Setup
-# Number of spatial elements: 101 1in horizontal, 15 in vertical
+# Number of spatial elements: 101 in horizontal, 15 in vertical
 # Soil depth: 50 m
-# Simulation duration: 365 d
+# Simulation duration: 2-20 years, based on LONGER_RUN setting
 # Timestep: 450 s
 # Timestepper: ARS111
-# Fixed number of iterations: 1
-# Jacobian update: every new timestep
+# Fixed number of iterations: 3
+# Jacobian update: every new Newton iteration
 # Atmos forcing update: every 3 hours
 import ClimaComms
 ClimaComms.@import_required_backends
@@ -28,7 +28,7 @@ import ClimaUtilities.TimeVaryingInputs: LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import ClimaUtilities.TimeManager: ITime
 import ClimaParams as CP
-
+using ClimaCore
 using ClimaLand
 using ClimaLand.Soil
 import ClimaLand
@@ -57,13 +57,25 @@ root_path = "soil_longrun_$(device_suffix)"
 diagnostics_outdir = joinpath(root_path, "global_diagnostics")
 outdir =
     ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir)
+time_interpolation_method =
+    LONGER_RUN ? LinearInterpolation() : LinearInterpolation(PeriodicCalendar())
 
-function setup_model(FT, start_date, stop_date, domain, earth_param_set)
-    time_interpolation_method =
-        LONGER_RUN ? LinearInterpolation() :
-        LinearInterpolation(PeriodicCalendar())
-    surface_space = domain.space.surface
-    subsurface_space = domain.space.subsurface
+function setup_simulation(; greet = false)
+    # If not LONGER_RUN, run for 2 years; note that the forcing from 2008 is repeated.
+    # If LONGER run, run for 20 years, with the correct forcing each year.
+    start_date = LONGER_RUN ? DateTime(2000) : DateTime(2008)
+    stop_date = LONGER_RUN ? DateTime(2020) : DateTime(2010)
+    Δt = 450.0
+    nelements = (101, 15)
+    if greet
+        @info "Run: Global Soil Model"
+        @info "Resolution: $nelements"
+        @info "Timestep: $Δt s"
+        @info "Start Date: $start_date"
+        @info "Stop Date: $stop_date"
+    end
+    domain = ClimaLand.Domains.global_domain(FT; context, nelements)
+    earth_param_set = LP.LandParameters(FT)
     # Forcing data
     if LONGER_RUN
         era5_ncdata_path = ClimaLand.Artifacts.find_era5_year_paths(
@@ -75,81 +87,16 @@ function setup_model(FT, start_date, stop_date, domain, earth_param_set)
         era5_ncdata_path =
             ClimaLand.Artifacts.era5_land_forcing_data2008_path(; context)
     end
-    atmos, radiation = ClimaLand.prescribed_forcing_era5(
+    forcing = ClimaLand.prescribed_forcing_era5(
         era5_ncdata_path,
-        surface_space,
+        domain.space.surface,
         start_date,
         earth_param_set,
         FT;
         max_wind_speed = 25.0,
         time_interpolation_method,
     )
-
-    (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) =
-        ClimaLand.Soil.soil_composition_parameters(subsurface_space, FT)
-    (; ν, hydrology_cm, K_sat, θ_r) =
-        ClimaLand.Soil.soil_vangenuchten_parameters(subsurface_space, FT)
-    soil_albedo = Soil.CLMTwoBandSoilAlbedo{FT}(;
-        ClimaLand.Soil.clm_soil_albedo_parameters(surface_space)...,
-    )
-    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
-    soil_params = Soil.EnergyHydrologyParameters(
-        FT;
-        ν,
-        ν_ss_om,
-        ν_ss_quartz,
-        ν_ss_gravel,
-        hydrology_cm,
-        K_sat,
-        S_s,
-        θ_r,
-        albedo = soil_albedo,
-    )
-    f_over = FT(3.28) # 1/m
-    R_sb = FT(1.484e-4 / 1000) # m/s
-    runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
-        f_over = f_over,
-        f_max = ClimaLand.Soil.topmodel_fmax(surface_space, FT),
-        R_sb = R_sb,
-    )
-
-    soil_args = (domain = domain, parameters = soil_params)
-    soil_model_type = Soil.EnergyHydrology{FT}
-    sources = (Soil.PhaseChange{FT}(),)# sublimation and subsurface runoff are added automatically
-    top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(
-        atmos,
-        radiation,
-        runoff_model,
-        (:soil,),
-    )
-    zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
-    boundary_conditions =
-        (; top = top_bc, bottom = Soil.EnergyWaterFreeDrainage())
-    soil = soil_model_type(;
-        boundary_conditions = boundary_conditions,
-        sources = sources,
-        soil_args...,
-    )
-    return soil
-end
-
-function setup_simulation(; greet = false)
-    # If not LONGER_RUN, run for 2 years; note that the forcing from 2008 is repeated.
-    # If LONGER run, run for 20 years, with the correct forcing each year.
-    start_date = LONGER_RUN ? DateTime(2000) : DateTime(2008)
-    stop_date = LONGER_RUN ? DateTime(2020) : DateTime(2010)
-    Δt = 450.0
-    nelements = (101, 15)
-    if greet
-        @info "Run: Global Soil Model"
-        @info "Resolution: $nelements"
-        @info "Timestep: $Δt s"
-        @info "Start Date: $start_date"
-        @info "Stop Date: $stop_date"
-    end
-    domain = ClimaLand.Domains.global_domain(FT; context, nelements)
-    params = LP.LandParameters(FT)
-    model = setup_model(FT, start_date, stop_date, domain, params)
+    model = ClimaLand.Soil.EnergyHydrology(FT, domain, forcing, earth_param_set)
     diagnostics = ClimaLand.default_diagnostics(
         model,
         start_date;

src/shared_utilities/drivers.jl

@@ -1506,7 +1506,7 @@ function prescribed_forcing_era5(
         earth_param_set = earth_param_set,
         frac_diff = frac_diff,
     )
-    return atmos, radiation
+    return (; atmos, radiation)
 end
 
 

src/simulations/Simulations.jl

@@ -84,7 +84,7 @@ function LandSimulation(
             Dates.Month(1),
             start_date,
             ITime(Δt, epoch = start_date),
-            mask = ClimaLand.landsea_mask(ClimaLand.get_domain(model)),
+            mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)),
         ),
         ClimaLand.ReportCallback(1000),
     ),

src/standalone/Soil/Soil.jl

@@ -188,4 +188,134 @@ include("./soil_hydrology_parameterizations.jl")
 include("./spatially_varying_parameters.jl")
 include("Biogeochemistry/Biogeochemistry.jl")
 using .Biogeochemistry
+
+
+# Soil model constructor useful for working with simulations forced by
+# the atmosphere
+"""
+    EnergyHydrology(FT, domain, forcing, earth_param_set;
+                         prognostic_land_components = (:soil),
+                         albedo = Soil.CLMTwoBandSoilAlbedo{FT}(; clm_soil_albedo_parameters(domain.space.surface)...),
+                         runoff =  ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(f_over = FT(3.28),
+                                                                            R_sb = FT(1.484e-4 / 1000),
+                                                                            f_max = topmodel_fmax(domain.space.surface,FT),
+                                                                            ),
+                         retention_parameters = soil_vangenuchten_parameters(domain.space.subsurface, FT),
+                         composition_parameters = soil_composition_parameters(domain.space.subsurface, FT),
+                         S_s = ClimaCore.Fields.zeros(domain.space.subsurface) .+ 1e-3,
+                         additional_sources = (),
+                         )
+
+Creates a EnergyHydrology model with the given float type FT, domain, earth_param_set, forcing, and prognostic land components.
+
+When running the soil model in standalone mode, `prognostic_land_components = (:soil,)`, while for running integrated land models,
+this should be a list of the component models. This value of this argument must be the same across all components in the land model.
+
+Default spatially varying parameters (for retention curve parameters, composition, and specific storativity) are provided but can be
+changed with keyword arguments. 
+
+The runoff and albedo parameterizations are also provided and can be changed via keyword argument; 
+additional sources may be required in your model if the soil model will be composed with other 
+component models.
+
+TODO: Move scalar parameters to ClimaParams and obtain from earth_param_set, possibly use types in retention and composition arguments.
+"""
+function EnergyHydrology(
+    FT,
+    domain,
+    forcing,
+    earth_param_set;
+    prognostic_land_components = (:soil,),
+    albedo::AbstractSoilAlbedoParameterization = CLMTwoBandSoilAlbedo{FT}(;
+        clm_soil_albedo_parameters(domain.space.surface)...,
+    ),
+    runoff::Runoff.AbstractRunoffModel = Runoff.TOPMODELRunoff{FT}(
+        f_over = FT(3.28), # extract from EPS
+        R_sb = FT(1.484e-4 / 1000),# extract from EPS
+        f_max = topmodel_fmax(domain.space.surface, FT),
+    ),
+    retention_parameters = soil_vangenuchten_parameters(
+        domain.space.subsurface,
+        FT,
+    ),
+    composition_parameters = soil_composition_parameters(
+        domain.space.subsurface,
+        FT,
+    ),
+    S_s = ClimaCore.Fields.zeros(domain.space.subsurface) .+ 1e-3,
+    additional_sources = (),
+)
+    top_bc = AtmosDrivenFluxBC(
+        forcing.atmos,
+        forcing.radiation,
+        runoff;
+        prognostic_land_components,
+    )
+    bottom_bc = EnergyWaterFreeDrainage()
+    boundary_conditions = (; top = top_bc, bottom = bottom_bc)
+    # sublimation and subsurface runoff are added automatically
+    sources = (additional_sources..., PhaseChange{FT}())
+    parameters = EnergyHydrologyParameters(
+        FT;
+        retention_parameters...,
+        composition_parameters...,
+        albedo,
+        S_s,
+    )
+    return EnergyHydrology{FT}(;
+        parameters,
+        domain,
+        boundary_conditions,
+        sources,
+    )
+end
+
+
+"""
+    RichardsModel(FT, domain, forcing;
+                         runoff =  ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(f_over = FT(3.28),
+                                                                            R_sb = FT(1.484e-4 / 1000),
+                                                                            f_max = topmodel_fmax(domain.space.surface,FT),
+                                                                            ),
+                         retention_parameters = soil_vangenuchten_parameters(domain.space.subsurface, FT),
+                         S_s = ClimaCore.Fields.zeros(domain.space.subsurface) .+ 1e-3,
+                         )
+
+Creates a RichardsModel model with the given float type FT, domain, earth_param_set and forcing.
+
+Default spatially varying parameters (for retention curve parameters and specific storativity) are provided but can be
+changed with keyword arguments. 
+
+The runoff parameterization is also provided and can be changed via keyword argument.
+
+TODO: Move scalar parameters to ClimaParams and obtain from earth_param_set, possibly use types in retention argument.
+"""
+function RichardsModel(
+    FT,
+    domain,
+    forcing;
+    runoff::Runoff.AbstractRunoffModel = Runoff.TOPMODELRunoff{FT}(
+        f_over = FT(3.28), # extract from EPS
+        R_sb = FT(1.484e-4 / 1000),# extract from EPS
+        f_max = topmodel_fmax(domain.space.surface, FT),
+    ),
+    retention_parameters = soil_vangenuchten_parameters(
+        domain.space.subsurface,
+        FT,
+    ),
+    S_s = ClimaCore.Fields.zeros(domain.space.subsurface) .+ 1e-3,
+)
+    top_bc = RichardsAtmosDrivenFluxBC(forcing.atmos, runoff)
+    bottom_bc = WaterFluxBC((p, t) -> 0.0)
+    boundary_conditions = (; top = top_bc, bottom = bottom_bc)
+    parameters = RichardsParameters(retention_parameters..., S_s)
+    return RichardsModel{FT}(;
+        parameters,
+        domain,
+        boundary_conditions,
+        sources = (),
+        lateral_flow = false,
+    )
+end
+
 end