Merge pull request #779 from CliMA/kd/snow_soil_option2

Integrated snow and soil modeling - Option 2

.buildkite/pipeline.yml

@@ -43,6 +43,10 @@ steps:
 
   - group: "Experiments"
     steps:
+      - label: "Snow/Soil Ozark"
+        command: "julia --color=yes --project=.buildkite experiments/integrated/fluxnet/snow_soil/simulation.jl"
+        artifact_paths: "experiments/integrated/fluxnet/snow_soil/*.png"
+
       - label: "Snow Col de Porte"
         command: "julia --color=yes --project=.buildkite experiments/standalone/Snow/snowmip_simulation.jl cdp"
         artifact_paths: "experiments/standalone/Snow/*png"

docs/list_tutorials.jl

@@ -1,9 +1,17 @@
 tutorials = [
-    "Running simulations" => [
+    "Running land simulations" => [
         "Bucket LSM" => [
             "standalone/Bucket/bucket_tutorial.jl",
             "standalone/Bucket/coupled_bucket.jl",
         ],
+        "Integrated soil+canopy modeling" => [
+            "Coupled Canopy and Soil" => "integrated/soil_canopy_tutorial.jl",
+        ],
+        "Handling interactions between model components" => [
+            "Adjusting boundary conditions for the soil" => "integrated/handling_soil_fluxes.jl",
+        ],
+    ],
+    "Running standalone component simulations" => [
         "Soil modeling" => [
             "Boundary conditions" => "standalone/Soil/boundary_conditions.jl",
             "Richards Equation" => "standalone/Soil/richards_equation.jl",
@@ -17,9 +25,6 @@ tutorials = [
         "Canopy modeling" => [
             "Standalone Canopy" => "standalone/Canopy/canopy_tutorial.jl",
         ],
-        "Integrated soil+canopy modeling" => [
-            "Coupled Canopy and Soil" => "integrated/soil_canopy_tutorial.jl",
-        ],
         "Snow Modeling" => [
             "standalone/Snow/base_tutorial.jl",
             "standalone/Snow/data_tutorial.jl",

docs/tutorials/integrated/handling_soil_fluxes.jl

@@ -0,0 +1,56 @@
+# # Background
+# When solving the partial differential equations describing the storage of water, ice, and energy in soil,
+# boundary conditions must be set for the soil liquid water and energy.
+# The tutorial on [boundary conditions for the soil model](@ref Soil/boundary_conditions.md) describes
+# the various boundary condition options we currently support. Many of these are suitable for use when
+# replicating laboratory experiments or other standalone soil settings - for example, a particular experiment
+# may fix the temperature or water content at the top of a soil column. However, when we wish to model the
+# soil interacting with the atmosphere, in standalone or integrated models, the particular boundary condition
+# we use is one that computes the sensible, latent, and radiative fluxes in addition to evaporation and sublimation;
+# this option is called the `AtmosDrivenFluxBC`. This boundary condition type includes the atmospheric and radiative
+# forcings, the runoff parameterization, and finally, a Tuple which indicates which components of the land are present.
+# The last argument is what we will describe here. For more information on how to supply the prescribed forcing
+# data, please see the tutorial linked [here](@ref shared_utilities/driver_tutorial.jl).
+
+# # Adjusting the boundary conditions for the soil model when run as part of an integrated land model
+
+# The presence of other land components (a canopy, snow) affects the boundary conditions of the soil
+# and changes them relative to what they would be if only bare soil was interacting with the atmosphere.
+# For example, the canopy and snow absorb radiation that might otherwise be absorbed by the soil, the
+# canopy intercepts precipitation, and the snow can melt and contribute to soil infiltration. Because of this,
+# we need a way to indicate to the soil model that the boundary conditions must be computed in a different way.
+# In all cases, the same `AtmosDrivenFluxBC` is used with the exception of the field `prognostic_land_components`,
+# which is a Tuple of the symbols associated with each component in the land model. For example,
+# if you are simulating the soil model in standalone mode, you would set
+
+# `prognostic_land_components = (:soil,)`
+
+# `top_soil_boundary_condition = ClimaLand.Soil.AtmosDrivenFluxBC(atmos_forcing, radiative_forcing, runoff_parameterization, prognostic_land_components)`
+
+# while, if you are simulating both a prognostic canopy and soil model, you would set
+
+# `prognostic_land_components = (:canopy, :soil)`
+
+# `top_soil_boundary_condition = ClimaLand.Soil.AtmosDrivenFluxBC(atmos_forcing, radiative_forcing, runoff_parameterization, prognostic_land_components)`
+
+# Note that the land components are *always* in alphabetical order, and the symbols associated with each
+# land model component are *always* the same as the name of that component, i.e., 
+
+# `ClimaLand.name(snow_model) = :snow`
+
+# `ClimaLand.name(canopy_model) = :canopy`
+
+# `ClimaLand.name(soilco2_model) = :soilco2`
+
+# `ClimaLand.name(soil_model) = :soil`
+
+# Currently, the only supported options are: `(:soil,), (:canopy, :soil, :soilco2), (:snow, :soil)`.
+
+# # How it works
+# When the soil model updates the boundary conditions, it calls the function
+# `soil_boundary_fluxes!(bc, ClimaLand.TopBoundary(), soil_model, Δz_top, Y, p, t)`, where `bc`
+# is the boundary condition being used at the top of the domain.
+# This function has multiple methods depending on the boundary condition type `typeof(bc)`. When that type is `AtmosDrivenFluxBC`,
+# the function then calls `soil_boundary_fluxes!(bc, Val(bc.prognostic_land_components), soil, Y, p, t)`.
+# Again multiple dispatch is used, this time to compute the fluxes correctly according to the value of the `prognostic_land_components`,
+# i.e., based on which components are present.

experiments/integrated/fluxnet/met_drivers_FLUXNET.jl

@@ -109,20 +109,44 @@ esat =
     )
 e = @. esat - drivers.VPD.values
 q = @. 0.622 * e ./ (drivers.PA.values - 0.378 * e)
-
+RH = @. e / esat
+
+"""
+    snow_precip_fraction(air_temp, hum)
+
+Estimate the fraction of precipitation that is in snow form,
+given the air temperature at the surface in K and the relative humidity
+(between 0 and 1).
+
+See Jennings, K.S., Winchell, T.S., Livneh, B. et al. 
+Spatial variation of the rain–snow temperature threshold across the 
+Northern Hemisphere. Nat Commun 9, 1148 (2018). 
+https://doi.org/10.1038/s41467-018-03629-7
+"""
+function snow_precip_fraction(air_temp, hum)
+    air_temp_C = air_temp - 273.15
+    α = -10.04
+    β = 1.41
+    γ = 0.09
+    snow_frac = (1.0 / (1.0 + exp(α + β * air_temp_C + γ * hum)))
+    return snow_frac
+end
+snow_frac = snow_precip_fraction.(drivers.TA.values[:], RH[:])
 # Create interpolators for each atmospheric driver needed for PrescribedAtmosphere and for
 # PrescribedRadiation
 seconds = FT.(0:DATA_DT:((length(UTC_DATETIME) - 1) * DATA_DT));
 
-precip = TimeVaryingInput(seconds, -FT.(drivers.P.values[:]); context) # m/s
+P_liq = -FT.(drivers.P.values[:] .* (1 .- snow_frac))
+precip = TimeVaryingInput(seconds, P_liq; context) # m/s
 atmos_q = TimeVaryingInput(seconds, FT.(q[:]); context)
 atmos_T = TimeVaryingInput(seconds, FT.(drivers.TA.values[:]); context)
 atmos_p = TimeVaryingInput(seconds, FT.(drivers.PA.values[:]); context)
 atmos_co2 = TimeVaryingInput(seconds, FT.(drivers.CO2.values[:]); context)
 atmos_u = TimeVaryingInput(seconds, FT.(drivers.WS.values[:]); context)
 LW_IN = TimeVaryingInput(seconds, FT.(drivers.LW_IN.values[:]); context)
 SW_IN = TimeVaryingInput(seconds, FT.(drivers.SW_IN.values[:]); context)
-snow_precip = TimeVaryingInput((t) -> FT(0))
+P_snow = -FT.(drivers.P.values[:] .* snow_frac)
+snow_precip = TimeVaryingInput(seconds, P_snow; context) # m/s
 
 # Construct the drivers. For the start date we will use the UTC time at the
 # start of the simulation

experiments/integrated/fluxnet/snow_soil/parameters_fluxnet.jl

@@ -0,0 +1,49 @@
+## 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.
+data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv"
+site_ID = "US-MOz"
+# Timezone (offset from UTC in hrs)
+time_offset = 7
+
+# Height of sensor on flux tower
+atmos_h = FT(32)
+
+# Site latitude and longitude
+lat = FT(38.7441) # degree
+long = FT(-92.2000) # degree
+# Soil parameters
+soil_ν = FT(0.5) # m3/m3
+soil_K_sat = FT(0.45 / 3600 / 100) # m/s,
+soil_S_s = FT(1e-3) # 1/m, guess
+soil_vg_n = FT(2.0) # unitless
+soil_vg_α = FT(2.0) # inverse meters
+θ_r = FT(0.09) # m3/m3, 
+
+# Soil heat transfer parameters; not needed for hydrology only test
+ν_ss_quartz = FT(0.2)
+ν_ss_om = FT(0.0)
+ν_ss_gravel = FT(0.4);
+κ_quartz = FT(7.7) # W/m/K
+κ_minerals = FT(2.5) # W/m/K
+κ_om = FT(0.25) # W/m/K
+κ_liq = FT(0.57) # W/m/K
+κ_ice = FT(2.29) # W/m/K
+κ_air = FT(0.025); #W/m/K
+ρp = FT(2700); # kg/m^3
+κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals)
+ρc_ds = FT((1 - soil_ν) * 4e6); # J/m^3/K
+z_0m_soil = FT(0.01)
+z_0b_soil = FT(0.001)
+soil_ϵ = FT(0.98)
+soil_α_PAR = FT(0.25)
+soil_α_NIR = FT(0.25)
+
+# 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)