Change soil albedo to depend of soil moisture (#862)

* Change soil albedo to depend of soil moisture

PAR and NIR albedo are moved from the EnergyHydrologyParameters
to the auxiliary vars via AtmosDrivenFluxBC. The parameters now
hold PAR/NIR wet and dry as a field of float. The old method
of creating the parameters with PAR/NIR albedo is marked as deprecated,
and the constructor just sets both the wet and dry values to the
passed albedo.

* Use relative saturation

* Make relative saturation in albedo calc not include ice

* Add soil albedos to diagnostics

* Make experiments use ground albedo maps

* Improve diagnostic var names

* Update docstring for PrognosticSoilConditions

* Make update_albedo! use effective_saturation

* Combine par and nir albedo diagnostics

* Refactor soil albedo calculation

update_albedo is no longer called conditionally. Now it is dispatched
to a function that either updates albedo or a function that does nothing.

update_albedo takes in model.params instead of each param individually.

update_albedo now just uses the top layer of if delta_z is larger than
the top depth.

top_depth is now a parameter instead of hard coded

cache var names are changed. Physically irrelevent vars changed to scratch vars

Delete integrated soil-canopy test

Previosuly, PAR and NIR albedo were held in a PrognosticSoilConditions
struct, which acted as an intermediary between the soil and canopy models.

Now, PAR and NIR albedo are stored in the cache, and the struct is only
used for dispatch.

The integrated soil-canopy test created a PrognosticSoilConditions
struct and filled in with a PAR and NIR albedo. Then it called
ground_albedo_PAR and ground_albedo_NIR to check the values.

Both of those functions now read from the cache. In order to test them, a
SoilCanopy model would need to be created, initialized, and then updated.
This seems to complex for a test, so this test is deleted.

Make run_fluxnet work with new SoilEnergy

Make update_albedo! comments better

Apply suggestions from code review

Co-authored-by: Gabriele Bozzola <gbozzola@caltech.edu>

Change top_depth name to alebdo_calc_top_thickness

Add convenience constructor for soil albedo maps

Fix calls to create_soil_albedo_vars

Fix spelling errors

Simplify and correct calculate_albedo!

Add citation and fix albedo variable constructor error

* Remake soil_canopy_lsm test

This test was deleted because it would not function with previous changes.
As requested, this test was recreated to work with the changes. The previous
test simply set an albedo it the top bc struct and then checked it using
the ground_albedo functions. With the changes, this would only work with an
initialized model. Because this is an integrated model test, I thought it
made sense to test the SoilCanopy model instead of just creating a
Soil model and checking if the function works. This test, like the
old version, calls the ground_albedo functions and checks their value.
This does seem very involved for a test, so maybe it should just create
a EnergyHydrology model and call the function on its Y and p

* Make suggested changes from Gabriele

change create_soil_albedo_vars to clm_soil_albedo_properties
and use map in it to reduce code duplication

made stylistic change when taking mean of PAR and NIR albedo

change <:FT to just FT

Move salb to short diagnostics

Also make other requested changes such as changing
surface_liquid_fraction to surface_eff_sat and making docstring more
specific

Fix albedo constructor and change benchmark sampling

Move salb disgnostic and reduce profiling samples

revert samples

* Merge soil param convenience constructors

Fix merge issue

fix merge issue 2

* Fix albedo params in snowy_land.jl

experiments/README.md

@@ -20,7 +20,7 @@ Many of the outputs from the experiments are saved using a non-destructive appro
 output directory of `example_outdir`. When run, if the directory does
 not exist relative to the working directory, then `example_outdir` is created. The results of the latest experiment
 run will be saved into `example_outdir/output_xxx`, where xxxx is an increasing counter. The lower indices store the outputs of previous runs. The latest results are also linked to
-in `example_outdir/output_active`, which can be assumed to always contain the most recent ouput. More details on this style of output directory handeling
+in `example_outdir/output_active`, which can be assumed to always contain the most recent output. More details on this style of output directory handling
 can be found [here](https://clima.github.io/ClimaUtilities.jl/dev/outputpathgenerator/#ActiveLinkStyle-(Non-Destructive))
 
 ### benchmarks
@@ -30,7 +30,7 @@ When a benchmark is run, its outputs are saved using `ActiveLinkStyle` into (ben
 
 ### long_runs
 
-The expirements in experiments/long_runs save its image outputs to (experiment_name)\_longrun\_(device_suffix).
+The experiments in experiments/long_runs save its image outputs to (experiment_name)\_longrun\_(device_suffix).
 All other outputs are saved using `ActiveLinkStyle` to (experiment_name)\_longrun\_(device_suffix)/global_diagnostics.
 
 ### integrated and standalone

experiments/benchmarks/land.jl

@@ -96,8 +96,10 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -110,9 +112,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
+
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -137,6 +142,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
 
     # Energy Balance model
     ac_canopy = FT(2.5e3)
+
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)

experiments/integrated/global/global_parameters.jl

@@ -13,8 +13,10 @@ spatially_varying_soil_params =
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
+    PAR_albedo_dry,
+    NIR_albedo_dry,
+    PAR_albedo_wet,
+    NIR_albedo_wet,
     f_max,
 ) = spatially_varying_soil_params
 soil_params = Soil.EnergyHydrologyParameters(
@@ -27,8 +29,10 @@ soil_params = Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo = PAR_albedo,
-    NIR_albedo = NIR_albedo,
+    PAR_albedo_dry = PAR_albedo_dry,
+    NIR_albedo_dry = NIR_albedo_dry,
+    PAR_albedo_wet = PAR_albedo_wet,
+    NIR_albedo_wet = NIR_albedo_wet,
 );
 f_over = FT(3.28) # 1/m
 R_sb = FT(1.484e-4 / 1000) # m/s

experiments/long_runs/land.jl

@@ -97,8 +97,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -111,9 +113,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
+
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;

experiments/long_runs/land_region.jl

@@ -100,8 +100,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -114,8 +116,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s

experiments/long_runs/snowy_land.jl

@@ -99,8 +99,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -113,8 +115,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s

experiments/long_runs/soil.jl

@@ -94,8 +94,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -108,8 +110,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
 
     f_over = FT(3.28) # 1/m
@@ -154,7 +158,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     # The approximation arises because the porosity, residual water fraction,
     # and van Genuchtem parameters are spatially varying but treated constant
     # in solving for equilibrium. Finally, we make a plausible but total guess
-    # for the water table depth using the TOPMODEL parameters. 
+    # for the water table depth using the TOPMODEL parameters.
     function hydrostatic_profile(
         lat::FT,
         z::FT,

ext/CreateParametersExt.jl

@@ -279,7 +279,9 @@ that certain parameters must have the same type (e.g, if a field is
 supplied for porosity, it must be supplied for all other parameters
 defined in the interior of the domain). Some parameters are defined only
 on the surface of the domain (e.g albedo), while other are defined everywhere
-(e.g. porosity). These are indicated with types `F` and `SF`.
+(e.g. porosity). These are indicated with types `F` and `SF`. If both dry/wet albedos
+and general albedos are given as keywords, the dry/wet albedos will override the general
+albedos.
 
 Please see the EnergyHydrologyParameters documentation for a complete list.
 """
@@ -319,12 +321,19 @@ function EnergyHydrologyParameters(
     K_sat::F,
     S_s::F,
     θ_r::F,
-    PAR_albedo::SF = 0.2,
-    NIR_albedo::SF = 0.4,
+    PAR_albedo_dry::SF = nothing,
+    NIR_albedo_dry::SF = nothing,
+    PAR_albedo_wet::SF = nothing,
+    NIR_albedo_wet::SF = nothing,
+    PAR_albedo::SFD = 0.2,
+    NIR_albedo::SFD = 0.4,
+    albedo_calc_top_thickness::TD = 0.07,
     kwargs...,
 ) where {
     F <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
-    SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
+    SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field, Nothing},
+    SFD <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
+    TD <: AbstractFloat,
     C,
 }
     earth_param_set = LP.LandParameters(toml_dict)
@@ -386,9 +395,19 @@ function EnergyHydrologyParameters(
     parameters = CP.get_parameter_values(toml_dict, name_map, "Land")
     PSE = typeof(earth_param_set)
     FT = CP.float_type(toml_dict)
-    EnergyHydrologyParameters{FT, F, SF, C, PSE}(;
-        PAR_albedo,
-        NIR_albedo,
+    if isnothing(PAR_albedo_dry)
+        PAR_albedo_dry = FT.(PAR_albedo)
+        NIR_albedo_dry = FT.(NIR_albedo)
+        PAR_albedo_wet = FT.(PAR_albedo)
+        NIR_albedo_wet = FT.(NIR_albedo)
+    end
+    albedo_calc_top_thickness = FT(albedo_calc_top_thickness)
+    EnergyHydrologyParameters{FT, F, typeof(PAR_albedo_dry), C, PSE}(;
+        PAR_albedo_wet,
+        NIR_albedo_wet,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        albedo_calc_top_thickness,
         ν,
         ν_ss_om,
         ν_ss_quartz,

land_longrun_gpu/global_diagnostics/output_active

@@ -0,0 +1 @@
+output_0110

lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl

@@ -21,6 +21,12 @@ function run_fluxnet(
     # Now we set up the model. For the soil model, we pick
     # a model type and model args:
     soil_domain = domain.land_domain
+    PAR_albedo =
+        FT(0.5) * params.soil.PAR_albedo_wet +
+        FT(0.5) * params.soil.PAR_albedo_dry
+    NIR_albedo =
+        FT(0.5) * params.soil.NIR_albedo_wet +
+        FT(0.5) * params.soil.NIR_albedo_dry
     soil_ps = Soil.EnergyHydrologyParameters(
         FT;
         ν = params.soil.ν,
@@ -34,8 +40,8 @@ function run_fluxnet(
         z_0m = params.soil.z_0m,
         z_0b = params.soil.z_0b,
         emissivity = params.soil.emissivity,
-        PAR_albedo = params.soil.PAR_albedo,
-        NIR_albedo = params.soil.NIR_albedo,
+        PAR_albedo = PAR_albedo,
+        NIR_albedo = NIR_albedo,
     )
 
     soil_args = (domain = soil_domain, parameters = soil_ps)

src/diagnostics/default_diagnostics.jl

@@ -177,6 +177,7 @@ function default_diagnostics(
             "lhf",
             "shf",
             "ghf",
+            "salb",
         ]
     elseif output_vars == :short
         soilcanopy_diagnostics = [

src/diagnostics/define_diagnostics.jl

@@ -605,6 +605,18 @@ function define_diagnostics!(land_model)
             compute_infiltration!(out, Y, p, t, land_model),
     )
 
+    # Soil albedo
+    add_diagnostic_variable!(
+        short_name = "salb",
+        long_name = "Soil Albedo",
+        standard_name = "surface albedo",
+        units = "",
+        comments = "The mean of PAR and NIR albedo, which are calculated as α_soil_band = α_band_dry * (1 - S_e) + α_band_wet * S_e.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_albedo!(out, Y, p, t, land_model),
+    )
+
+
     # Soil hydraulic conductivity
     add_diagnostic_variable!(
         short_name = "shc",

src/diagnostics/land_compute_methods.jl

@@ -173,6 +173,20 @@ end
 @diagnostic_compute "soil_net_radiation" Union{SoilCanopyModel, LandModel} p.soil.R_n
 @diagnostic_compute "soil_temperature" Union{SoilCanopyModel, LandModel} p.soil.T
 
+function compute_soil_albedo!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel{FT},
+) where {FT}
+    if isnothing(out)
+        return (p.soil.PAR_albedo .+ p.soil.NIR_albedo) ./ 2
+    else
+        @. out = (p.soil.PAR_albedo + p.soil.NIR_albedo) / 2
+    end
+end
+
 # Soil - Turbulent Fluxes
 @diagnostic_compute "soil_latent_heat_flux" Union{SoilCanopyModel, LandModel} p.soil.turbulent_fluxes.lhf
 @diagnostic_compute "soil_sensible_heat_flux" Union{SoilCanopyModel, LandModel} p.soil.turbulent_fluxes.shf