change albedo parameterization to enforce physicality (#1124)

Author: kmdeck

experiments/calibration/forward_model.jl

@@ -64,7 +64,7 @@ function setup_prob(
         #        α_leaf_scaler,
         #        α_soil_scaler,
         α_0,
-        α_horizon,
+        Δα,
         k,
         beta_snow,
         x0_snow,
@@ -318,10 +318,10 @@ function setup_prob(
     )
 
     # Snow model
-    # α_snow = Snow.ZenithAngleAlbedoModel(α_0, α_horizon, k)
+    # α_snow = Snow.ZenithAngleAlbedoModel(α_0, Δα, k)
     α_snow = Snow.ZenithAngleAlbedoModel(
         FT(α_0),
-        FT(α_horizon),
+        FT(Δα),
         FT(k);
         β = FT(beta_snow),
         x0 = FT(x0_snow),

experiments/calibration/priors.jl

@@ -14,7 +14,7 @@
 #    EKP.constrained_gaussian("α_leaf_scaler", 1, 0.15, 0.5, 1.5);
 
 prior_α_0 = EKP.constrained_gaussian("α_0", 0.5, 0.2, 0.2, 0.8);
-prior_α_horizon = EKP.constrained_gaussian("α_horizon", 0.85, 0.1, 0.7, 1);
+prior_Δα = EKP.constrained_gaussian("Δα", 0.2, 0.1, 0.0, 1.0);
 prior_k = EKP.constrained_gaussian("k", 10, 5, 2, 25);
 
 prior_beta_snow = EKP.constrained_gaussian("beta_snow", 0.4, 0.2, 0.1, 0.8);
@@ -35,7 +35,7 @@ prior = EKP.combine_distributions([
     #    # prior_α_snow,
     #    prior_α_soil_scale,
     prior_α_0,
-    prior_α_horizon,
+    prior_Δα,
     prior_k,
     prior_beta_snow,
     prior_x0_snow,

src/standalone/Snow/Snow.jl

@@ -108,9 +108,9 @@ end
 
 Establishes the albedo parameterization where albedo
 depends on the cosine of the zenith angle of the sun, as
-    α = f(x) * [α_0 + (α_horizon - α_0)*exp(-k*cos(θs))],
+    α = f(x) * [α_0 + Δα*exp(-k*cos(θs))],
 
-where cos θs is the cosine of the zenith angle, α_0, α_1, and k 
+where cos θs is the cosine of the zenith angle, α_0, Δα, and k 
 are free parameters. The factor out front is a function of 
 x = ρ_snow/ρ_liq, of the form f(x) = min(1 - β(x-x0), 1). The parameters
 x0 ∈ [0,1] and β ∈ [0,1] are free. Choose β = 0 to remove this dependence on snow density.
@@ -122,9 +122,9 @@ 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"
+    "Free parameter controlling the snow albedo when θs = 90∘"
+    Δα::FT
+    "Rate at which albedo drops to its minimum value with zenith angle"
     k::FT
     "Rate governing how snow albedo changes with snow density, a proxy for grain size and liquid water content, ∈ [0,1]"
     β::FT
@@ -134,14 +134,14 @@ end
 
 function ZenithAngleAlbedoModel(
     α_0::FT,
-    α_horizon::FT,
+    Δα::FT,
     k::FT;
     β = FT(0),
     x0 = FT(0.2),
 ) where {FT}
     @assert 0 ≤ x0 ≤ 1
     @assert 0 ≤ β ≤ 1
-    ZenithAngleAlbedoModel(α_0, α_horizon, k, β, x0)
+    ZenithAngleAlbedoModel(α_0, Δα, k, β, x0)
 end
 
 """

src/standalone/Snow/snow_parameterizations.jl

@@ -47,10 +47,9 @@ function update_snow_albedo!(
     FT = eltype(earth_param_set)
     _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
     @. α =
-        min(1 - m.β * (p.snow.ρ_snow / _ρ_liq - m.x0), 1) * (
-            m.α_0 +
-            (m.α_horizon - m.α_0) * exp(-m.k * max(p.drivers.cosθs, eps(FT)))
-        )
+        min(1 - m.β * (p.snow.ρ_snow / _ρ_liq - m.x0), 1) *
+        (m.α_0 + m.Δα * exp(-m.k * max(p.drivers.cosθs, eps(FT))))
+    @. α = max(min(α, 1), 0)
 end
 """
     update_snow_cover_fraction!(x::FT; z0 = FT(1e-1), β_scf = FT(2))::FT where {FT}

test/standalone/Snow/parameterizations.jl

@@ -129,7 +129,7 @@ for FT in (Float32, Float64)
                min(m.β_scf * depth / FT(0.05) / (depth / FT(0.05) + 1), 1)
         )
 
-        m = Snow.ZenithAngleAlbedoModel(FT(0.6), FT(0.7), FT(2); β = FT(0.7))
+        m = Snow.ZenithAngleAlbedoModel(FT(0.6), FT(0.1), FT(2); β = FT(0.7))
         update_snow_albedo!(α_snow, m, nothing, p, 0.0, param_set)
         @test all(
             @. α_snow ≈
@@ -144,7 +144,7 @@ for FT in (Float32, Float64)
             for b in β
                 m = Snow.ZenithAngleAlbedoModel(
                     FT(0.6),
-                    FT(0.6),
+                    FT(0.0),
                     FT(0);
                     β = b,
                     x0 = x,