Snow model constructor, Soil/snow integrated models (#1220)

Author: kmdeck

experiments/integrated/fluxnet/snow_soil/parameters_fluxnet.jl

@@ -1,3 +1,11 @@
+## Some site parameters have been taken from
+## Ma, S., Baldocchi, D. D., Xu, L., Hehn, T. (2007)
+## Inter-Annual Variability In Carbon Dioxide Exchange Of An
+## Oak/Grass Savanna And Open Grassland In California, Agricultural
+## And Forest Meteorology, 147(3-4), 157-171. https://doi.org/10.1016/j.agrformet.2007.07.008 
+## CLM 5.0 Tech Note: https://www2.cesm.ucar.edu/models/cesm2/land/CLM50_Tech_Note.pdf
+# Bonan, G. Climate change and terrestrial ecosystem modeling. Cambridge University Press, 2019.
+
 import SciMLBase
 import ClimaTimeSteppers as CTS
 import ClimaComms
@@ -19,6 +27,11 @@ import ClimaParams
 climaland_dir = pkgdir(ClimaLand)
 
 FT = Float64
+site_ID = "US-MOz"
+time_offset = 7 # offset from UTC in hrs
+lat = FT(38.7441) # degree
+long = FT(-92.2000) # degree
+
 earth_param_set = LP.LandParameters(FT)
 # Utility functions for reading in and filling fluxnet data
 include(joinpath(climaland_dir, "experiments/integrated/fluxnet/data_tools.jl"))
@@ -30,7 +43,7 @@ zmax = FT(0)
 dz_bottom = FT(1.0)
 dz_top = FT(0.04)
 
-land_domain = Column(;
+domain = Column(;
     zlim = (zmin, zmax),
     nelements = nelements,
     dz_tuple = (dz_bottom, dz_top),
@@ -43,13 +56,16 @@ N_days = N_days_spinup + 360
 dt = FT(900)
 tf = t0 + FT(3600 * 24 * N_days)
 
-# Read in the parameters for the site
-include(
-    joinpath(
-        climaland_dir,
-        "experiments/integrated/fluxnet/snow_soil/parameters_fluxnet.jl",
-    ),
-)
+
+data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv"
+# Height of sensor on flux tower
+atmos_h = FT(32)
+
+# Required to use the same met_drivers_FLUXNET script, even though we currently have no canopy
+h_leaf = FT(0)
+f_root_to_shoot = FT(0)
+plant_ν = FT(0)
+
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
 include(
@@ -58,57 +74,54 @@ include(
         "experiments/integrated/fluxnet/met_drivers_FLUXNET.jl",
     ),
 )
-
-
-# Now we set up the model. For the soil model, we pick
-# a model type and model args:
-soil_domain = land_domain
-soil_albedo = Soil.ConstantTwoBandSoilAlbedo{FT}(;
-    PAR_albedo = soil_α_PAR,
-    NIR_albedo = soil_α_NIR,
+forcing = (; atmos, radiation)
+prognostic_land_components = (:snow, :soil)
+α_soil = Soil.ConstantTwoBandSoilAlbedo{FT}(;
+    PAR_albedo = FT(0.25),
+    NIR_albedo = FT(0.25),
 )
-soil_ps = Soil.EnergyHydrologyParameters(
-    FT;
-    ν = soil_ν,
-    ν_ss_om,
-    ν_ss_quartz,
-    ν_ss_gravel,
-    hydrology_cm = vanGenuchten{FT}(; α = soil_vg_α, n = soil_vg_n),
-    K_sat = soil_K_sat,
-    S_s = soil_S_s,
-    θ_r,
-    z_0m = z_0m_soil,
-    z_0b = z_0b_soil,
-    emissivity = soil_ϵ,
-    albedo = soil_albedo,
+runoff = ClimaLand.Soil.SurfaceRunoff()
+retention_parameters = (;
+    ν = FT(0.5),
+    K_sat = FT(0.45 / 3600 / 100),
+    hydrology_cm = vanGenuchten{FT}(; α = FT(2), n = FT(2)),
+    θ_r = FT(0.09),
 )
-
-soil_args = (parameters = soil_ps,)
-soil_model_type = Soil.EnergyHydrology
-
-ρ = 300.0
-α = 0.8
-snow_parameters = SnowParameters{FT}(
-    dt;
-    α_snow = Snow.ConstantAlbedoModel(α),
-    density = Snow.MinimumDensityModel(ρ),
-    earth_param_set = earth_param_set,
-);
-snow_args = (; parameters = snow_parameters);
-snow_model_type = Snow.SnowModel
-land_input = (
-    atmos = atmos,
-    radiation = radiation,
-    domain = land_domain,
-    runoff = ClimaLand.Soil.SurfaceRunoff(),
+composition_parameters =
+    (; ν_ss_quartz = FT(0.2), ν_ss_om = FT(0.0), ν_ss_gravel = FT(0.4))
+S_s = FT(1e-3)
+z_0m = FT(0.01)
+z_0b = FT(0.001)
+emissivity = FT(0.98)
+soil_model = Soil.EnergyHydrology(
+    FT,
+    domain,
+    forcing,
+    earth_param_set;
+    prognostic_land_components,
+    albedo = α_soil,
+    runoff,
+    retention_parameters,
+    S_s,
+    composition_parameters,
+    z_0m,
+    z_0b,
+    emissivity,
 )
-land = ClimaLand.SoilSnowModel{FT}(;
-    land_args = land_input,
-    soil_model_type = soil_model_type,
-    soil_args = soil_args,
-    snow_args = snow_args,
-    snow_model_type = snow_model_type,
+α_snow = Snow.ConstantAlbedoModel(0.8)
+density = Snow.MinimumDensityModel(300.0)
+snow_model = Snow.SnowModel(
+    FT,
+    ClimaLand.Domains.obtain_surface_domain(domain),
+    forcing,
+    earth_param_set,
+    dt;
+    prognostic_land_components,
+    α_snow,
+    density,
 )
+
+land = ClimaLand.SoilSnowModel(; snow = snow_model, soil = soil_model)
 Y, p, cds = initialize(land)
 exp_tendency! = make_exp_tendency(land)
 imp_tendency! = make_imp_tendency(land)
@@ -131,7 +144,7 @@ T_0 =
     volumetric_heat_capacity.(
         Y.soil.ϑ_l,
         Y.soil.θ_i,
-        soil_ps.ρc_ds,
+        soil_model.parameters.ρc_ds,
         earth_param_set,
     )
 Y.soil.ρe_int =

experiments/integrated/fluxnet/snow_soil/simulation.jl

@@ -95,6 +95,7 @@ atmos = PrescribedAtmosphere(
     h_atmos,
     earth_param_set,
 )
+forcing = (; atmos, radiation)
 # Snow data
 snow_data = NCDataset(joinpath(path, obs_filename))
 

experiments/standalone/Snow/process_snowmip.jl

@@ -50,28 +50,28 @@ ndays = (tf - t0) / 3600 / 24
 
 domain = ClimaLand.Domains.Point(; z_sfc = FT(0))
 
-density_model = NeuralSnow.NeuralDepthModel(FT)
-#density_model = Snow.MinimumDensityModel(ρ)
+density = NeuralSnow.NeuralDepthModel(FT)
+#density = Snow.MinimumDensityModel(ρ)
+α_snow = Snow.ConstantAlbedoModel(α)
 
-parameters = SnowParameters{FT}(
-    Δt;
-    α_snow = Snow.ConstantAlbedoModel(α),
-    density = density_model,
-    earth_param_set = param_set,
-)
 model = ClimaLand.Snow.SnowModel(
-    parameters = parameters,
-    domain = domain,
-    boundary_conditions = ClimaLand.Snow.AtmosDrivenSnowBC(atmos, radiation),
+    FT,
+    domain,
+    forcing,
+    earth_param_set,
+    Δt;
+    density,
+    α_snow,
 )
+
 Y, p, coords = ClimaLand.initialize(model)
 
 # Set initial conditions
 Y.snow.S .= FT(SWE[1]) # first data point
 Y.snow.S_l .= 0 # this is a guess
 Y.snow.Z .= FT(depths[1]) #first depth value - comment out if using MinimumDensityModel instead of NeuralDepthModel 
 Y.snow.U .=
-    ClimaLand.Snow.energy_from_q_l_and_swe(FT(SWE[1]), FT(0), parameters) # with q_l = 0
+    ClimaLand.Snow.energy_from_q_l_and_swe(FT(SWE[1]), FT(0), model.parameters) # with q_l = 0
 
 set_initial_cache! = ClimaLand.make_set_initial_cache(model)
 set_initial_cache!(p, Y, t0)

experiments/standalone/Snow/snowmip_simulation.jl

@@ -15,7 +15,11 @@ using ClimaCore.Operators: column_integral_definite!
     end
 
 A concrete type of land model used for simulating systems with
-snow and soil (and eventually rivers).
+snow and soil.
+
+The inner constructor checks that the two models are consistent with
+respect to the forcing (atmos, radiation), the parameters, the domain,
+and the prognostic land components of the model.
 $(DocStringExtensions.FIELDS)
 """
 struct SoilSnowModel{
@@ -27,75 +31,25 @@ struct SoilSnowModel{
     snow::SnM
     "The soil model to be used"
     soil::SoM
-end
-
-"""
-    SoilSnowModel{FT}(;
-        land_args::NamedTuple = (;),
-        snow_model_type::Type{SnM},
-        snow_args::NamedTuple = (;),
-        soil_model_type::Type{SoM},
-        soil_args::NamedTuple = (;),
-        ) where {
-            FT,
-            SnM <: Snow.SnowModel{FT},
-            SoM <: Soil.EnergyHydrology{FT},
-            }
-
-A constructor for the `LandHydrology`, which takes in the concrete model
-type and required arguments for each component, constructs those models,
-and constructs the `SoilSnowModel` from them.
-
-Each component model is constructed with everything it needs to be stepped
-forward in time, including boundary conditions, source terms, and interaction
-terms.
-"""
-function SoilSnowModel{FT}(;
-    land_args::NamedTuple = (;),
-    snow_model_type::Type{SnM},
-    snow_args::NamedTuple = (;),
-    soil_model_type::Type{SoM},
-    soil_args::NamedTuple = (;),
-) where {FT, SnM <: Snow.SnowModel, SoM <: Soil.EnergyHydrology}
-    (; atmos, radiation, domain) = land_args
-    prognostic_land_components = (:snow, :soil)
-    if :runoff ∈ propertynames(land_args)
-        runoff_model = land_args.runoff
-    else
-        runoff_model = ClimaLand.Soil.Runoff.NoRunoff()
+    function SoilSnowModel(; snow, soil)
+        prognostic_land_components = (:snow, :soil)
+        top_soil_bc = soil.boundary_conditions.top
+        snow_bc = snow.boundary_conditions
+        @assert top_soil_bc.prognostic_land_components ==
+                prognostic_land_components
+        @assert snow_bc.prognostic_land_components == prognostic_land_components
+
+        @assert top_soil_bc.atmos == snow_bc.atmos
+        @assert top_soil_bc.radiation == snow_bc.radiation
+
+        @assert Domains.obtain_surface_domain(soil.domain) == snow.domain
+
+        @assert soil.parameters.earth_param_set ==
+                snow.parameters.earth_param_set
+        FT = eltype(soil.parameters.earth_param_set)
+        new{FT, typeof(snow), typeof(soil)}(snow, soil)
     end
-    top_bc = ClimaLand.AtmosDrivenFluxBC(
-        atmos,
-        radiation,
-        runoff_model,
-        prognostic_land_components,
-    )
-    sources = (Soil.PhaseChange{FT}(),)
-    zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
-    boundary_conditions = (;
-        top = top_bc,
-        bottom = Soil.WaterHeatBC(;
-            water = Soil.FreeDrainage(),
-            heat = zero_flux,
-        ),
-    )
-    soil = soil_model_type{FT}(;
-        boundary_conditions = boundary_conditions,
-        sources = sources,
-        domain = domain,
-        soil_args...,
-    )
-    snow = snow_model_type(;
-        boundary_conditions = Snow.AtmosDrivenSnowBC(
-            atmos,
-            radiation,
-            prognostic_land_components,
-        ),
-        domain = Domains.obtain_surface_domain(domain),
-        snow_args...,
-    )
 
-    return SoilSnowModel{FT, typeof(snow), typeof(soil)}(snow, soil)
 end
 
 """