Adds a snow albedo which varies with zenith angle (#1100)

* add albedo option to snow

* address ordering issue

Author: kmdeck

.dev/Project.toml

@@ -2,4 +2,5 @@
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 
 [compat]
-JuliaFormatter = "1"
+JuliaFormatter = "1"
+

docs/src/APIs/Snow.md

@@ -31,4 +31,11 @@ ClimaLand.Snow.snow_cover_fraction
 ClimaLand.Snow.snow_boundary_fluxes!
 ClimaLand.Snow.phase_change_flux
 ClimaLand.Snow.AtmosDrivenSnowBC
-```
+```
+
+## Snow parameterizations
+
+```@docs
+ClimaLand.Snow.ZenithAngleAlbedoModel
+ClimaLand.Snow.ConstantAlbedoModel
+```

experiments/integrated/fluxnet/snow_soil/simulation.jl

@@ -87,7 +87,7 @@ soil_model_type = Soil.EnergyHydrology
 α = 0.8
 snow_parameters = SnowParameters{FT}(
     dt;
-    α_snow = α,
+    α_snow = Snow.ConstantAlbedoModel(α),
     density = Snow.MinimumDensityModel(ρ),
     earth_param_set = earth_param_set,
 );

experiments/long_runs/snowy_land.jl

@@ -312,7 +312,13 @@ function setup_prob(
     )
 
     # Snow model
-    snow_parameters = SnowParameters{FT}(Δt; earth_param_set = earth_param_set)
+    #    α_snow = Snow.ConstantAlbedoModel(FT(0.7))
+    α_snow = Snow.ZenithAngleAlbedoModel(FT(0.5), FT(0.9), FT(1))
+    snow_parameters = SnowParameters{FT}(
+        Δt;
+        earth_param_set = earth_param_set,
+        α_snow = α_snow,
+    )
     snow_args = (;
         parameters = snow_parameters,
         domain = ClimaLand.obtain_surface_domain(domain),

experiments/standalone/Snow/snowmip_simulation.jl

@@ -55,7 +55,7 @@ density_model = NeuralSnow.NeuralDepthModel(FT)
 
 parameters = SnowParameters{FT}(
     Δt;
-    α_snow = α,
+    α_snow = Snow.ConstantAlbedoModel(α),
     density = density_model,
     earth_param_set = param_set,
 )

src/integrated/land.jl

@@ -128,10 +128,7 @@ function LandModel{FT}(;
     )
 
     transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}()
-    ground_conditions =
-        PrognosticGroundConditions{FT, typeof(snow.parameters.α_snow)}(
-            snow.parameters.α_snow,
-        )
+    ground_conditions = PrognosticGroundConditions()
     if :energy in propertynames(canopy_component_args)
         energy_model = canopy_component_types.energy(
             canopy_component_args.energy.parameters,
@@ -179,6 +176,20 @@ function LandModel{FT}(;
     return LandModel{FT, typeof.(args)...}(args...)
 end
 
+"""
+   ClimaLand.land_components(land::LandModel)
+
+Returns the components of the `LandModel`. 
+
+Currently, this method is required in order to preserve an ordering in how
+we update the component models' auxiliary states. The canopy update_aux! step
+depends on snow and soil albedo, but those are only updated in the snow and soil
+update_aux! steps. So those must occur first (as controlled by the order of the components
+returned by `land_components!`.
+
+This needs to be fixed.
+"""
+ClimaLand.land_components(land::LandModel) = (:soil, :snow, :soilco2, :canopy)
 """
     lsm_aux_vars(m::LandModel)
 
@@ -416,8 +427,8 @@ function lsm_radiant_energy_fluxes!(
     ϵ_soil = land.soil.parameters.emissivity
     T_soil = ClimaLand.Domains.top_center_to_surface(p.soil.T)
 
-    α_snow_NIR = land.snow.parameters.α_snow
-    α_snow_PAR = land.snow.parameters.α_snow
+    α_snow_NIR = p.snow.α_snow
+    α_snow_PAR = p.snow.α_snow
     ϵ_snow = land.snow.parameters.ϵ_snow
     T_snow = p.snow.T_sfc
 
@@ -579,25 +590,16 @@ function ClimaLand.Soil.sublimation_source(
 end
 
 """
-     PrognosticGroundConditions{FT <: AbstractFloat, F <: Union{FT, ClimaCore.Fields.Field}} <: Canopy.AbstractGroundConditions
+     PrognosticGroundConditions <: Canopy.AbstractGroundConditions
 
 A type of Canopy.AbstractGroundConditions to use when the soil model is prognostic and
 of type `EnergyHydrology`, and the snow model is prognostic and included.
 
-The canopy model needs albedo of the ground
-in order to compute its update_aux! function, and that function must only depend on the canopy model.
-Because of this, α_snow must be stored in this struct until it is stored in the cache.
-
-Note that this struct is linked with the EnergyHydrology model. If we ever had a different
+Note that this struct is linked with the EnergyHydrology/SnowModel model. If we ever had a different
 soil model, we might need to construct a different `PrognosticGroundConditions` because
 the fields may be stored in different places.
 """
-struct PrognosticGroundConditions{
-    FT <: AbstractFloat,
-    F <: Union{FT, ClimaCore.Fields.Field},
-} <: Canopy.AbstractGroundConditions
-    α_snow::FT
-end
+struct PrognosticGroundConditions <: Canopy.AbstractGroundConditions end
 
 """
     Canopy.ground_albedo_PAR(
@@ -620,7 +622,7 @@ function Canopy.ground_albedo_PAR(
 )
     @. p.α_ground.PAR =
         (1 - p.snow.snow_cover_fraction) * p.soil.PAR_albedo +
-        p.snow.snow_cover_fraction * ground.α_snow
+        p.snow.snow_cover_fraction * p.snow.α_snow
     return p.α_ground.PAR
 end
 
@@ -645,7 +647,7 @@ function Canopy.ground_albedo_NIR(
 )
     @. p.α_ground.NIR =
         (1 - p.snow.snow_cover_fraction) * p.soil.NIR_albedo +
-        p.snow.snow_cover_fraction * ground.α_snow
+        p.snow.snow_cover_fraction * p.snow.α_snow
     return p.α_ground.NIR
 end
 

src/shared_utilities/models.jl

@@ -484,12 +484,16 @@ This can be converted to a mass using the density of liquid water.
 This includes the water in multiple phases. For example, if ice is present, the water
 volume is computed using ratio of the density of ice to the density of liquid water.
 """
-function total_liq_water_vol_per_area!(cache, model::AbstractModel, Y, p, t) end
+function total_liq_water_vol_per_area!(cache, model::AbstractModel, Y, p, t)
+    cache .= 0
+end
 
 """
     total_energy_per_area!(cache, model::AbstractModel, Y, p, t)
 
 A function which updates `cache` in place with the total energy
 per unit ground area for the `model`, computed from `Y`, `p`, and `t`.
 """
-function total_energy_per_area!(cache, model::AbstractModel, Y, p, t) end
+function total_energy_per_area!(cache, model::AbstractModel, Y, p, t)
+    cache .= 0
+end

src/standalone/Snow/Snow.jl

@@ -33,7 +33,12 @@ import ClimaLand:
     get_drivers,
     total_energy_per_area!,
     total_liq_water_vol_per_area!
-export SnowParameters, SnowModel, AtmosDrivenSnowBC, snow_boundary_fluxes!
+export SnowParameters,
+    SnowModel,
+    AtmosDrivenSnowBC,
+    snow_boundary_fluxes!,
+    ConstantAlbedoModel,
+    ZenithAngleAlbedoModel
 
 """
     AbstractSnowModel{FT} <: ClimaLand.AbstractExpModel{FT}
@@ -73,6 +78,50 @@ struct MinimumDensityModel{FT} <: AbstractDensityModel{FT}
 end
 
 
+"""
+    AbstractAlbedoModel{FT}
+
+Defines the model type for albedo parameterization
+for use within an `AbstractSnowModel` type. 
+
+These parameterizations are stored in parameters.α_snow, and 
+are used to update the value of p.snow.α_snow (the broadband
+albedo of the snow at a point).
+stored 
+"""
+abstract type AbstractAlbedoModel{FT <: AbstractFloat} end
+
+"""
+    ConstantAlbedoModel{FT <: AbstractFloat} <: AbstractAlbedoModel{FT}
+
+Establishes the albedo parameterization where albedo is treated as a
+constant spatially and temporally.
+"""
+struct ConstantAlbedoModel{FT} <: AbstractAlbedoModel{FT}
+    "Albedo of snow (unitless)"
+    α::FT
+end
+
+"""
+    ZenithAngleAlbedoModel{FT <: AbstractFloat} <: AbstractAlbedoModel{FT}
+
+Establishes the albedo parameterization where albedo only
+depends on the cosine of the zenith angle of the sun, as
+    α = α_0 + (α_horizon - α_0)*exp(-k*cos(θs))
+
+Note: If this choice is used, the field cosθs must appear in the cache
+p.drivers. This is available through the PrescribedRadiativeFluxes object.
+"""
+struct ZenithAngleAlbedoModel{FT} <: AbstractAlbedoModel{FT}
+    "Free parameter controlling the minimum snow albedo"
+    α_0::FT
+    "The albedo of snow when the sun is at the horizon"
+    α_horizon::FT
+    "Rate at which albedo drops from α_horizon to its minimum value"
+    k::FT
+end
+
+
 """
     SnowParameters{FT <: AbstractFloat, PSE}
 
@@ -90,6 +139,7 @@ $(DocStringExtensions.FIELDS)
 Base.@kwdef struct SnowParameters{
     FT <: AbstractFloat,
     DM <: AbstractDensityModel,
+    AM <: AbstractAlbedoModel,
     PSE,
 }
     "Choice of parameterization for snow density"
@@ -98,8 +148,8 @@ Base.@kwdef struct SnowParameters{
     z_0m::FT
     "Roughness length over snow for scalars (m)"
     z_0b::FT
-    "Albedo of snow (unitless)"
-    α_snow::FT
+    "Albedo parameterization for snow"
+    α_snow::AM
     "Emissivity of snow (unitless)"
     ϵ_snow::FT
     "Volumetric holding capacity of water in snow (unitless)"
@@ -121,7 +171,7 @@ end
                       density = MinimumDensityModel(200),
                       z_0m = FT(0.0024),
                       z_0b = FT(0.00024),
-                      α_snow = FT(0.8),
+                      α_snow = ConstantAlbedoModel(FT(0.8)),
                       ϵ_snow = FT(0.99),
                       θ_r = FT(0.08),
                       Ksat = FT(1e-3),
@@ -137,15 +187,20 @@ function SnowParameters{FT}(
     density::DM = MinimumDensityModel(FT(200)),
     z_0m = FT(0.0024),
     z_0b = FT(0.00024),
-    α_snow = FT(0.8),
+    α_snow::AM = ConstantAlbedoModel(FT(0.8)),
     ϵ_snow = FT(0.99),
     θ_r = FT(0.08),
     Ksat = FT(1e-3),
     κ_ice = FT(2.21),
     ΔS = FT(0.1),
     earth_param_set::PSE,
-) where {FT <: AbstractFloat, DM <: AbstractDensityModel, PSE}
-    return SnowParameters{FT, DM, PSE}(
+) where {
+    FT <: AbstractFloat,
+    DM <: AbstractDensityModel,
+    AM <: AbstractAlbedoModel,
+    PSE,
+}
+    return SnowParameters{FT, DM, AM, PSE}(
         density,
         z_0m,
         z_0b,
@@ -273,6 +328,7 @@ auxiliary_vars(snow::SnowModel) = (
     :T,
     :T_sfc,
     :z_snow,
+    :α_snow,
     :ρ_snow,
     :R_n,
     :phase_change_flux,
@@ -305,6 +361,7 @@ auxiliary_types(snow::SnowModel{FT}) where {FT} = (
     FT,
     FT,
     FT,
+    FT,
     boundary_var_types(
         snow,
         snow.boundary_conditions,
@@ -330,6 +387,7 @@ auxiliary_domain_names(snow::SnowModel) = (
     :surface,
     :surface,
     :surface,
+    :surface,
     boundary_var_domain_names(
         snow.boundary_conditions,
         ClimaLand.TopBoundary(),
@@ -342,6 +400,7 @@ ClimaLand.name(::SnowModel) = :snow
 function ClimaLand.make_update_aux(model::SnowModel{FT}) where {FT}
     function update_aux!(p, Y, t)
         parameters = model.parameters
+        update_snow_albedo!(p.snow.α_snow, parameters.α_snow, Y, p, t)
 
         # This has to happen first, since other quantities below depend
         # on it.

src/standalone/Snow/snow_parameterizations.jl

@@ -10,8 +10,31 @@ export snow_surface_temperature,
     energy_from_T_and_swe,
     snow_cover_fraction,
     snow_bulk_density,
-    phase_change_flux
+    phase_change_flux,
+    update_snow_albedo!
 
+"""
+    update_snow_albedo!(α, m::ConstantAlbedoModel, Y, p, t)
+
+Updates the snow albedo `α` in place with the current albedo,
+according to the ConstantAlbedoModel. 
+"""
+function update_snow_albedo!(α, m::ConstantAlbedoModel, Y, p, t)
+    @. α = m.α
+end
+
+"""
+    update_snow_albedo!(α, m::ZenithAngleAlbedoModel, Y, p, t)
+
+Updates the snow albedo `α` in place with the current albedo,
+according to the ZenithAngleAlbedoModel.
+"""
+function update_snow_albedo!(α, m::ZenithAngleAlbedoModel, Y, p, t)
+    FT = FTfromY(Y)
+    @. α =
+        m.α_0 +
+        (m.α_horizon - m.α_0) * exp(-m.k * max(p.drivers.cosθs, eps(FT)))
+end
 """
     snow_cover_fraction(x::FT; z0 = FT(1e-1), α = FT(2))::FT where {FT}
 
@@ -54,8 +77,8 @@ end
 
 A helper function which computes and returns the snow albedo.
 """
-function ClimaLand.surface_albedo(model::SnowModel, _...)
-    return model.parameters.α_snow
+function ClimaLand.surface_albedo(model::SnowModel, Y, p)
+    return p.snow.α_snow
 end
 
 """

test/integrated/full_land.jl

@@ -19,7 +19,7 @@ f_over = FT(3.28) # 1/m
 R_sb = FT(1.484e-4 / 1000) # m/s
 scalar_soil_params = (; f_over, R_sb)
 
-α_snow = FT(0.67)
+α_snow = Snow.ConstantAlbedoModel(FT(0.67))
 scalar_snow_params = (; α_snow, Δt)
 
 # Energy Balance model
@@ -82,6 +82,7 @@ land = global_land_model(
     LAI = LAI,
 )
 
+@test ClimaLand.land_components(land) == (:soil, :snow, :soilco2, :canopy)
 Y, p, cds = initialize(land)
 
 # Soil IC