Change soil albedo to depend of soil moisture

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/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

src/integrated/land.jl

@@ -17,7 +17,7 @@ export LandModel
         snow::SnM
     end
 
-A concrete type of land model used for simulating systems with 
+A concrete type of land model used for simulating systems with
 soil, canopy, snow, soilco2.
 $(DocStringExtensions.FIELDS)
 """
@@ -129,10 +129,8 @@ function LandModel{FT}(;
 
     transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}()
     ground_conditions =
-        PrognosticGroundConditions{FT, typeof(soil.parameters.PAR_albedo)}(
+        PrognosticGroundConditions{FT, typeof(snow.parameters.α_snow)}(
             snow.parameters.α_snow,
-            soil.parameters.PAR_albedo,
-            soil.parameters.NIR_albedo,
         )
     if :energy in propertynames(canopy_component_args)
         energy_model = canopy_component_types.energy(
@@ -260,7 +258,7 @@ lsm_aux_domain_names(m::LandModel) = (
 A method which makes a function; the returned function
 updates the additional auxiliary variables for the integrated model,
 as well as updates the boundary auxiliary variables for all component
-models. 
+models.
 
 This function is called each ode function evaluation, prior to the tendency function
 evaluation.
@@ -293,7 +291,7 @@ function make_update_boundary_fluxes(
             t,
         )
 
-        # Effective (radiative) land properties 
+        # Effective (radiative) land properties
         set_eff_land_radiation_properties!(
             p,
             land.soil.parameters.earth_param_set,
@@ -396,8 +394,8 @@ function lsm_radiant_energy_fluxes!(
     T_canopy =
         ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p, t)
 
-    α_soil_PAR = land.soil.parameters.PAR_albedo
-    α_soil_NIR = land.soil.parameters.NIR_albedo
+    α_soil_PAR = p.soil.PAR_albedo
+    α_soil_NIR = p.soil.NIR_albedo
     ϵ_soil = land.soil.parameters.emissivity
     T_soil = ClimaLand.Domains.top_center_to_surface(p.soil.T)
     α_snow_NIR = land.snow.parameters.α_snow
@@ -579,8 +577,6 @@ struct PrognosticGroundConditions{
     F <: Union{FT, ClimaCore.Fields.Field},
 } <: Canopy.AbstractGroundConditions
     α_snow::FT
-    α_soil_PAR::F
-    α_soil_NIR::F
 end
 
 """
@@ -601,7 +597,7 @@ function Canopy.ground_albedo_PAR(
     p,
     t,
 )
-    return @. (1 - p.snow.snow_cover_fraction) * ground.α_soil_PAR +
+    return @. (1 - p.snow.snow_cover_fraction) * p.soil.PAR_albedo +
               p.snow.snow_cover_fraction * ground.α_snow
 end
 
@@ -623,7 +619,7 @@ function Canopy.ground_albedo_NIR(
     p,
     t,
 )
-    return @. (1 - p.snow.snow_cover_fraction) * ground.α_soil_NIR +
+    return @. (1 - p.snow.snow_cover_fraction) * p.soil.NIR_albedo +
               p.snow.snow_cover_fraction * ground.α_snow
 end
 
@@ -632,7 +628,7 @@ end
 
 Returns the "drivers", or forcing variables, for the LandModel.
 
-These consist of atmospheric and radiative forcing, as well as 
+These consist of atmospheric and radiative forcing, as well as
 soil organic carbon.
 """
 function ClimaLand.get_drivers(model::LandModel)