Merge pull request #796 from CliMA/kd/soil_snow_heat_flux

Add a ground heat flux between snow and soil

NEWS.md

@@ -3,6 +3,8 @@ ClimaLand.jl Release Notes
 
 main
 --------
+- Add ground heat flux to snow-soil model PR[#796](https://github.com/CliMA/ClimaLand.jl/pull/796)
+- Add snow-soil model PR [#779](https://github.com/CliMA/ClimaLand.jl/pull/779)
 
 v0.15.0
 --------

docs/list_tutorials.jl

@@ -9,6 +9,7 @@ tutorials = [
         ],
         "Handling interactions between model components" => [
             "Adjusting boundary conditions for the soil" => "integrated/handling_soil_fluxes.jl",
+            "Adjusting boundary conditions for the snow" => "integrated/handling_snow_fluxes.jl",
         ],
     ],
     "Running standalone component simulations" => [

docs/src/APIs/Snow.md

@@ -23,4 +23,12 @@ ClimaLand.Snow.runoff_timescale
 ClimaLand.Snow.compute_water_runoff	
 ClimaLand.Snow.energy_from_q_l_and_swe
 ClimaLand.Snow.energy_from_T_and_swe
+ClimaLand.Snow.snow_cover_fraction
+```
+
+## Computing fluxes for snow
 
+```@docs
+ClimaLand.Snow.snow_boundary_fluxes!
+ClimaLand.Snow.AtmosDrivenSnowBC
+```

docs/tutorials/integrated/handling_snow_fluxes.jl

@@ -0,0 +1,64 @@
+# # Background
+# When solving the system of equations describing the storage of water and energy in snow,
+# "boundary conditions" must be set. These are only true boundary conditions
+# if one considers the equations for water and energy in snow to be discretized PDEs,
+# but they do represent the fluxes at upper boundary (the snow/atmosphere interface) and at
+# the lower boundary (the snow/soil interface), so we refer to them as boundary conditions
+# in order to use a consistent notation with the soil model.
+
+# At present, the only type of boundary condition the snow model supports is
+# `AtmosDrivenSnowBC`, which is an attempt at a compact way of indicating that with this
+# choice, the snow model exchanges water and energy with the atmosphere via sensible, latent, and radiative heat fluxes, and by
+# precipitation and sublimation/evaporation. With this choice, the snow model also has exchange fluxes with the soil (melt water,
+# and a conductive ground heat flux). This boundary condition type includes the atmospheric and radiative
+# forcings, and, importantly, 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). For an explanation of the
+# same design implementation for the Soil model, please see [here](@ref integrated/handling_soil_fluxes.jl).
+
+# # Adjusting the boundary conditions for the snow model when run as part of an integrated land model
+
+# The presence of other land components (a canopy, the soil) affects the boundary conditions of the snow
+# and changes them relative to what they would be if only bare snow was interacting with the atmosphere.
+# For example, the canopy absorbs radiation that might otherwise be absorbed by the snow, the
+# canopy intercepts precipitation, and the soil and snow interact thermally at the interface between them.
+# Because of this, we need a way to indicate to the snow model that the boundary conditions must be computed in a different way
+# depending on the components. We do this via the field in the boundary condition of type `AtmosDrivenSnowBC` via
+# a 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 snow model in standalone mode, you would set
+
+# `prognostic_land_components = (:snow,)`
+
+# `boundary_condition = ClimaLand.Snow.AtmosDrivenSnowBC(atmos_forcing, radiative_forcing, prognostic_land_components)`
+
+# while, if you are simulating both a prognostic snow and soil model, you would set
+
+# `prognostic_land_components = (:snow, :soil)`
+
+# `boundary_condition = ClimaLand.Snow.AtmosDrivenSnowBC(atmos_forcing, radiative_forcing, 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: `(:snow,)`, and `(:snow, :soil)`.
+
+# # How it works
+# When the snow model computes its boundary fluxes, it calls the function
+# `snow_boundary_fluxes!(bc, snow_model, Y, p, t)`, where `bc`
+# is the boundary condition of type `AtmosDrivenSnowBC`.
+# This function then calls `snow_boundary_fluxes!(bc, Val(bc.prognostic_land_components), snow, Y, p, t)`.
+# Here multiple dispatch is used to compute the fluxes correctly according to the value of the `prognostic_land_components`,
+# i.e., based on which components are present.
+
+# # Some important notes
+# When the snow model is run in standalone mode (`prognostic_land_components = (:snow,)`), the ground heat flux is approximate
+# as zero. When the soil and snow model are run together, this flux is computed and affects both the snow and soil.

experiments/integrated/fluxnet/snow_soil/simulation.jl

@@ -40,7 +40,7 @@ land_domain = Column(;
 t0 = FT(1800)
 N_days_spinup = 0
 N_days = N_days_spinup + 360
-dt = FT(180)
+dt = FT(900)
 tf = t0 + FT(3600 * 24 * N_days)
 
 # Read in the parameters for the site
@@ -97,7 +97,7 @@ land_input = (
     atmos = atmos,
     radiation = radiation,
     domain = land_domain,
-    runoff = ClimaLand.Soil.NoRunoff(),#SiteLevelSurfaceRunoff(),
+    runoff = ClimaLand.Soil.SurfaceRunoff(),
 )
 land = ClimaLand.LandHydrologyModel{FT}(;
     land_args = land_input,
@@ -184,14 +184,17 @@ ax1 = Axis(fig[2, 2], ylabel = "SWC", xlabel = "Days")
 lines!(
     ax1,
     daily,
-    [parent(sol.u[k].soil.ϑ_l)[end] for k in 1:1:length(sol.t)],
-    label = "2cm",
+    [parent(sol.u[k].soil.ϑ_l)[end - 2] for k in 1:1:length(sol.t)],
+    label = "10cm",
 )
 lines!(
     ax1,
     daily,
-    [parent(sol.u[k].soil.θ_i)[end] for k in 1:1:length(sol.t)],
-    label = "2cm, ice",
+    [
+        parent(sol.u[k].soil.θ_i .+ sol.u[k].soil.ϑ_l)[end - 2] for
+        k in 1:1:length(sol.t)
+    ],
+    label = "10cm, liq+ice",
 )
 
 lines!(
@@ -214,11 +217,12 @@ lines!(ax3, daily, [parent(sol.u[k].snow.S)[1] for k in 1:1:length(sol.t)])
 # Temp
 ax4 = Axis(fig[1, 1], ylabel = "T (K)")
 hidexdecorations!(ax4, ticks = false)
+lines!(ax4, seconds ./ 3600 ./ 24, drivers.TA.values[:], label = "Data, Air")
 lines!(
     ax4,
     sv.t ./ 24 ./ 3600,
-    [parent(sv.saveval[k].soil.T)[end] for k in 1:1:length(sv.t)],
-    label = "Soil, 2cm",
+    [parent(sv.saveval[k].soil.T)[end - 2] for k in 1:1:length(sv.t)],
+    label = "Model 10 cm",
 )
 
 lines!(
@@ -227,8 +231,12 @@ lines!(
     [parent(sv.saveval[k].snow.T)[1] for k in 1:1:length(sv.t)],
     label = "Snow",
 )
-lines!(ax4, seconds ./ 3600 ./ 24, drivers.TS.values[:], label = "Data, ?cm")
-lines!(ax4, seconds ./ 3600 ./ 24, drivers.TA.values[:], label = "Data, Air")
+lines!(
+    ax4,
+    seconds ./ 3600 ./ 24,
+    drivers.TS.values[:],
+    label = "Data, Unknown depth",
+)
 axislegend(ax4, position = :rt)
 CairoMakie.save(joinpath(savedir, "results.png"), fig)
 
@@ -255,7 +263,7 @@ E_measured = [
         -1 .* [
             parent(
                 sv.saveval[k].atmos_water_flux .-
-                sv.saveval[k].soil.bottom_bc.water,
+                sv.saveval[k].soil.bottom_bc.water .+ sv.saveval[k].soil.R_s,
             )[end] for k in 1:1:(length(sv.t) - 1)
         ],
     ) * (sv.t[2] - sv.t[1])

experiments/standalone/Snow/snowmip_simulation.jl

@@ -50,8 +50,7 @@ parameters =
 model = ClimaLand.Snow.SnowModel(
     parameters = parameters,
     domain = domain,
-    atmos = atmos,
-    radiation = radiation,
+    boundary_conditions = ClimaLand.Snow.AtmosDrivenSnowBC(atmos, radiation),
 )
 Y, p, coords = ClimaLand.initialize(model)
 

src/ClimaLand.jl

@@ -303,6 +303,7 @@ import .Soil: soil_boundary_fluxes!, sublimation_source
 import .Soil.Biogeochemistry: soil_temperature, soil_moisture
 include("standalone/Snow/Snow.jl")
 using .Snow
+import .Snow: snow_boundary_fluxes!
 include("standalone/Vegetation/Canopy.jl")
 using .Canopy
 using .Canopy.PlantHydraulics