small changes to bring model into alignment with paper

CliMA · Nov 25, 2024 · e5599d3 · e5599d3
1 parent f69a484
commit e5599d3
Show file tree

Hide file tree

Showing 10 changed files with 48 additions and 57 deletions.
diff --git a/src/integrated/land.jl b/src/integrated/land.jl
@@ -481,15 +481,11 @@ function soil_boundary_fluxes!(
 ) where {FT}
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
-    influx = @. p.drivers.P_liq + p.snow.water_runoff * p.snow.snow_cover_fraction + p.excess_water_flux + (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
-    Soil.Runoff.update_runoff!(
-        p,
-        bc.runoff,
-        influx,
-        Y,
-        t,
-        soil,
-    )
+    influx = @. p.drivers.P_liq +
+       p.snow.water_runoff * p.snow.snow_cover_fraction +
+       p.excess_water_flux +
+       (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
+    Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
     @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =
         (1 - p.snow.snow_cover_fraction) * (

diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
@@ -370,9 +370,15 @@ function soil_boundary_fluxes!(
 ) where {FT}
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
-    Soil.Runoff.update_runoff!(p, bc.runoff, p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq, Y, t, soil)
-    @. p.soil.top_bc.water =
-        p.soil.infiltration
+    Soil.Runoff.update_runoff!(
+        p,
+        bc.runoff,
+        p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq,
+        Y,
+        t,
+        soil,
+    )
+    @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =
         -p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf
 end

diff --git a/src/integrated/soil_snow_model.jl b/src/integrated/soil_snow_model.jl
@@ -387,18 +387,12 @@ function soil_boundary_fluxes!(
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
     p.soil.R_n .= net_radiation(bc.radiation, soil, Y, p, t)
-    influx = @. p.drivers.P_liq + p.snow.water_runoff * p.snow.snow_cover_fraction +
-        p.excess_water_flux + (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
-    Soil.Runoff.update_runoff!(
-        p,
-        bc.runoff,
-        influx,
-        Y,
-        t,
-        soil,
-    )
-    @. p.soil.top_bc.water =
-        p.soil.infiltration
+    influx = @. p.drivers.P_liq +
+       p.snow.water_runoff * p.snow.snow_cover_fraction +
+       p.excess_water_flux +
+       (1 - p.snow.snow_cover_fraction) * p.soil.turbulent_fluxes.vapor_flux_liq
+    Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
+    @. p.soil.top_bc.water = p.soil.infiltration
 
     @. p.soil.top_bc.heat =
         (1 - p.snow.snow_cover_fraction) * (

diff --git a/src/standalone/Snow/snow_parameterizations.jl b/src/standalone/Snow/snow_parameterizations.jl
@@ -173,7 +173,8 @@ function snow_thermal_conductivity(
     _ρ_ice = FT(LP.ρ_cloud_ice(parameters.earth_param_set))
     κ_ice = parameters.κ_ice
     return _κ_air +
-           (FT(0.07) * (ρ_snow/_ρ_ice) + FT(0.93) * (ρ_snow/_ρ_ice)^2) * (κ_ice - _κ_air)
+           (FT(0.07) * (ρ_snow / _ρ_ice) + FT(0.93) * (ρ_snow / _ρ_ice)^2) *
+           (κ_ice - _κ_air)
 end
 
 """

diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
@@ -825,10 +825,16 @@ function soil_boundary_fluxes!(
 )
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, model, Y, p, t)
     p.soil.R_n .= net_radiation(bc.radiation, model, Y, p, t)
-    update_runoff!(p, bc.runoff, p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq, Y, t, model)
+    update_runoff!(
+        p,
+        bc.runoff,
+        p.drivers.P_liq .+ p.soil.turbulent_fluxes.vapor_flux_liq,
+        Y,
+        t,
+        model,
+    )
     # We do not model the energy flux from infiltration.
-    @. p.soil.top_bc.water =
-        p.soil.infiltration
+    @. p.soil.top_bc.water = p.soil.infiltration
     @. p.soil.top_bc.heat =
         p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf
 end

diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl
@@ -556,7 +556,7 @@ function ClimaLand.make_update_aux(model::EnergyHydrology)
         @. p.soil.κ = thermal_conductivity(
             model.parameters.κ_dry,
             kersten_number(
-                Y.soil.θ_i,
+                Y.soil.θ_i / (ν - θ_r),
                 relative_saturation(p.soil.θ_l, Y.soil.θ_i, ν),
                 α,
                 β,

diff --git a/src/standalone/Soil/soil_heat_parameterizations.jl b/src/standalone/Soil/soil_heat_parameterizations.jl
@@ -214,7 +214,7 @@ end
 
 """
     kersten_number(
-        θ_i::FT,
+        θ_i_over_ν::FT,
         S_r::FT,
         α::FT,
         β::FT,
@@ -227,24 +227,21 @@ Compute the expression for the Kersten number, using the Balland
 and Arp model.
 """
 function kersten_number(
-    θ_i::FT,
+    θ_i_over_ν::FT,
     S_r::FT,
     α::FT,
     β::FT,
     ν_ss_om::FT,
     ν_ss_quartz::FT,
     ν_ss_gravel::FT,
 ) where {FT}
-    if θ_i < eps(FT)
-        K_e =
-            S_r^((FT(1) + ν_ss_om - α * ν_ss_quartz - ν_ss_gravel) / FT(2)) *
-            (
-                (FT(1) + exp(-β * S_r))^(-FT(3)) -
-                ((FT(1) - S_r) / FT(2))^FT(3)
-            )^(FT(1) - ν_ss_om)
-    else
-        K_e = S_r^(FT(1) + ν_ss_om)
-    end
+    K_e_unfrozen =
+        S_r^((FT(1) + ν_ss_om - α * ν_ss_quartz - ν_ss_gravel) / FT(2)) *
+        (
+            (FT(1) + exp(-β * S_r))^(-FT(3)) - ((FT(1) - S_r) / FT(2))^FT(3)
+        )^(FT(1) - ν_ss_om)
+    K_e_frozen = S_r^(FT(1) + ν_ss_om)
+    K_e = (1 - θ_i_over_ν) * K_e_unfrozen + θ_i_over_ν * K_e_frozen
     return K_e
 end
 

diff --git a/test/standalone/Snow/parameterizations.jl b/test/standalone/Snow/parameterizations.jl
@@ -73,7 +73,8 @@ for FT in (Float32, Float64)
         @test specific_heat_capacity(FT(0.0), parameters) == _cp_i
         @test snow_thermal_conductivity(ρ_snow, parameters) ==
               κ_air +
-              (FT(0.07) * (ρ_snow/_ρ_i) + FT(0.93) * (ρ_snow/_ρ_i)^2) * (κ_ice - κ_air)
+              (FT(0.07) * (ρ_snow / _ρ_i) + FT(0.93) * (ρ_snow / _ρ_i)^2) *
+              (κ_ice - κ_air)
         @test runoff_timescale.(z, Ksat, FT(Δt)) ≈ max.(Δt, z ./ Ksat)
 
 

diff --git a/test/standalone/Soil/soil_bc.jl b/test/standalone/Soil/soil_bc.jl
@@ -172,7 +172,7 @@ for FT in (Float32, Float64)
             thermal_conductivity.(
                 parameters.κ_dry,
                 kersten_number.(
-                    θ_i,
+                    θ_i / (ν - θ_r),
                     relative_saturation.(ϑ_c, θ_i, ν),
                     parameters.α,
                     parameters.β,

diff --git a/test/standalone/Soil/soil_parameterizations.jl b/test/standalone/Soil/soil_parameterizations.jl
@@ -129,34 +129,24 @@ for FT in (Float32, Float64)
         @test relative_saturation(FT(0.25), FT(0.05), FT(0.4)) ==
               FT((0.25 + 0.05) / 0.4)
 
-        # ice fraction = 0
         @test kersten_number(
-            FT(0),
+            FT(0.05),
             FT(0.75),
             parameters.α,
             parameters.β,
             ν_ss_om,
             ν_ss_quartz,
             ν_ss_gravel,
-        ) ≈ FT(
+        ) ≈
+              FT(0.05 * 0.75^(FT(1) + 0.1)) +
+              FT(1 - 0.05) * FT(
             0.75^((FT(1) + 0.1 - 0.24 * 0.1 - 0.1) / FT(2)) *
             (
                 (FT(1) + exp(-18.3 * 0.75))^(-FT(3)) -
                 ((FT(1) - 0.75) / FT(2))^FT(3)
             )^(FT(1) - 0.1),
         )
 
-        # ice fraction ~= 0
-        @test kersten_number(
-            FT(0.05),
-            FT(0.75),
-            parameters.α,
-            parameters.β,
-            ν_ss_om,
-            ν_ss_quartz,
-            ν_ss_gravel,
-        ) == FT(0.75^(FT(1) + 0.1))
-
         @test thermal_conductivity(FT(1.5), FT(0.7287), FT(0.7187)) ==
               FT(0.7287 * 0.7187 + (FT(1) - 0.7287) * 1.5)