make albedo parameterization modular and make it match CLM (#1184)

Author: kmdeck

NEWS.md

@@ -3,7 +3,9 @@ ClimaLand.jl Release Notes
 
 main
 -------
-- Use new pun up initial conditions from 19 year run PR[#1196](https://github.com/CliMA/ClimaLand.jl/pull/1196)
+- ![breaking change][badge-💥breaking] Make soil albedo parameterization modular
+PR[#1184](https://github.com/CliMA/ClimaLand.jl/pull/1184)
+- Use new spun up initial conditions from 19 year run PR[#1196](https://github.com/CliMA/ClimaLand.jl/pull/1196)
 
 v0.16.3
 -------

docs/list_of_apis.jl

@@ -10,6 +10,7 @@ apis = [
         "Soil" => [
             "Soil Energy and Hydrology" => "APIs/Soil.md",
             "Soil Biogeochemistry" => "APIs/SoilBiogeochemistry.md",
+            "Soil Albedo" => "APIs/SoilAlbedo.md",
         ],
         "Canopy" => [
             "Canopy Models" => "APIs/canopy/Canopy.md",

docs/src/APIs/Soil.md

@@ -64,6 +64,14 @@ ClimaLand.Soil.subsurface_runoff_source
 ClimaLand.Soil.update_runoff!
 ```
 
+## Soil Albedo Types and Methods
+
+```@docs
+ClimaLand.Soil.CLMTwoBandSoilAlbedo
+ClimaLand.Soil.ConstantTwoBandSoilAlbedo
+ClimaLand.Soil.update_albedo!
+```
+
 ## Soil BC Methods and Types
 
 ```@docs

docs/src/APIs/SoilAlbedo.md

@@ -0,0 +1,75 @@
+## Parameterizing Soil Albedo
+
+# We currently support two parameterizations for the soil albedo. These
+# options are defined using Julia types, using simplified code, as follows.
+# We first introduce the abstract type `AbstractSoilAlbedoParameterization`,
+# and then the two albedo parameterizations are concrete examples of that
+# abstract type:
+
+```julia
+ConstantTwoBandSoilAlbedo <: AbstractSoilAlbedoParameterization
+CLMTwoBandSoilAlbedo <: AbstractSoilAlbedoParameterization
+```
+
+# Each of these types define a method for `update_albedo!`, which is
+# called in the soils tendency functions, and each stores the
+# parameters required to compute the albedo for that parameterizations.
+# For example, the parameterization based off of CLM's approach requires
+# the wet and dry soil albedos in the PAR and NIR wavelength bands. For
+# each band, the actual albedo is a linear combination of the wet and dry
+# values, with a weight depending on the soil moisture. The method
+# of update_albedo! for this parameterization computes this linear
+# combination using a helper function `albedo_from_moisture` and ultimately
+# sets
+
+```julia
+@. p.soil.PAR_albedo =
+    albedo_from_moisture(S_sfc, PAR_albedo_dry, PAR_albedo_wet)
+@. p.soil.NIR_albedo =
+    albedo_from_moisture(S_sfc, NIR_albedo_dry, NIR_albedo_wet)
+````
+where `S_sfc` is the effective saturation at the surface.
+
+### Creating a new albedo parameterization
+
+# Suppose you want to define a new parameterization which models albedo
+# as a linear combination of the albedos of quartz, organic matter, minerals
+# and water, treating NIR and PAR albedos the same.
+
+# First, create the type:
+
+```julia
+struct SoilAlbedoFromComposition{FT <: AbstractFloat} <: AbstractSoilAlbedoParameterization
+    α_quartz::FT
+    α_minerals::FT
+    α_om::FT
+    α_water::FT
+end
+```
+
+# And then create the method. For now, don't worry about the other arguments
+# and their types:
+
+```julia
+function update_albedo!(
+    bc::AtmosDrivenFluxBC,
+    albedo::SoilAlbedoFromComposition,
+    p,
+    soil_domain,
+    model_parameters,
+    )
+    # unpack parameters of the albedo model
+    (; α_quartz, α_minerals, α_om, α_water) = albedo
+    # unpack composition parameters from the soil parameters
+    (; ν_ss_om, ν_ss_quartz) = model_parameters
+    # Extract the values at the top level
+    ν_ss_om_sfc = ClimaLand.Domains.top_center_to_surface(ν_ss_om)
+    ν_ss_quartz_sfc = ClimaLand.Domains.top_center_to_surface(ν_ss_quartz)
+    # Extract the top layer's volumetric water content
+    θ_l_sfc = ClimaLand.Domains.top_center_to_surface(p.soil.θ_l)
+
+    # Compute the linear combination
+    @. p.soil.PAR_albedo = α_water * θ_l_sfc + (1-θ_l_sfc)* (α_quartz * ν_ss_quartz_sfc + α_om * ν_ss_om_sfc + α_minerals * (1-ν_ss_om - ν_ss_quartz))
+    p.soil.NIR_albedo .= p.soil.PAR_albedo 
+end
+```

docs/tutorials/integrated/soil_canopy_tutorial.jl

@@ -147,6 +147,10 @@ soil_α_PAR = FT(0.2)
 soil_α_NIR = FT(0.4)
 
 soil_domain = land_domain
+soil_albedo = ClimaLand.Soil.ConstantTwoBandSoilAlbedo{FT}(;
+    PAR_albedo = soil_α_PAR,
+    NIR_albedo = soil_α_NIR,
+)
 soil_ps = Soil.EnergyHydrologyParameters(
     FT;
     ν = soil_ν,
@@ -161,8 +165,7 @@ soil_ps = Soil.EnergyHydrologyParameters(
     z_0m = z_0m_soil,
     z_0b = z_0b_soil,
     emissivity = soil_ϵ,
-    PAR_albedo = soil_α_PAR,
-    NIR_albedo = soil_α_NIR,
+    albedo = soil_albedo,
 );
 
 soil_args = (domain = soil_domain, parameters = soil_ps)

docs/tutorials/standalone/Soil/evaporation.jl

@@ -102,8 +102,6 @@ S_s = FT(1e-3)
 ν_ss_quartz = FT(1.0)
 ν_ss_gravel = FT(0.0)
 emissivity = FT(1.0)
-PAR_albedo = FT(0.2)
-NIR_albedo = FT(0.4)
 z_0m = FT(1e-2)
 z_0b = FT(1e-2)
 d_ds = FT(0.01)
@@ -117,8 +115,6 @@ params = ClimaLand.Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
     emissivity,
     z_0m,
     z_0b,

docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl

@@ -47,8 +47,6 @@ S_s = FT(1e-3)
 ν_ss_quartz = FT(0.3)
 ν_ss_gravel = FT(0.0)
 emissivity = FT(1.0)
-PAR_albedo = FT(0.2)
-NIR_albedo = FT(0.4)
 z_0m = FT(1e-3)
 z_0b = FT(1e-4)
 d_ds = FT(0.01)# 10mm
@@ -62,8 +60,6 @@ params = ClimaLand.Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
     emissivity,
     z_0m,
     z_0b,

docs/tutorials/standalone/Soil/sublimation.jl

@@ -88,8 +88,6 @@ S_s = FT(1e-3)
 ν_ss_quartz = FT(1.0)
 ν_ss_gravel = FT(0.0)
 emissivity = FT(1.0)
-PAR_albedo = FT(0.2)
-NIR_albedo = FT(0.4)
 z_0m = FT(1e-3)
 z_0b = FT(1e-4)
 d_ds = FT(0.01)
@@ -103,8 +101,6 @@ params = ClimaLand.Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
     emissivity,
     z_0m,
     z_0b,

experiments/benchmarks/snowy_land.jl

@@ -116,6 +116,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
+    albedo = Soil.CLMTwoBandSoilAlbedo{FT}(;
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
+    )
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -126,10 +132,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo_dry = PAR_albedo_dry,
-        NIR_albedo_dry = NIR_albedo_dry,
-        PAR_albedo_wet = PAR_albedo_wet,
-        NIR_albedo_wet = NIR_albedo_wet,
+        albedo,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s

experiments/calibration/forward_model_land.jl

@@ -120,7 +120,12 @@ function setup_prob(
     #        min.(NIR_albedo_dry .* α_soil_dry_scaler .* α_soil_scaler, FT(1))
     #    PAR_albedo_wet .= min.(PAR_albedo_wet .* α_soil_scaler, FT(1))
     #    NIR_albedo_wet .= min.(NIR_albedo_wet .* α_soil_scaler, FT(1))
-
+    albedo = Soil.CLMTwoBandSoilAlbedo{FT}(;
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
+    )
     soil_params = Soil.EnergyHydrologyParameters(
         FT;
         ν,
@@ -131,10 +136,7 @@ function setup_prob(
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo_dry = PAR_albedo_dry,
-        NIR_albedo_dry = NIR_albedo_dry,
-        PAR_albedo_wet = PAR_albedo_wet,
-        NIR_albedo_wet = NIR_albedo_wet,
+        albedo,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s