From e5599d3a2245d437b5626f9064f5bf1d9453eeaf Mon Sep 17 00:00:00 2001 From: kmdeck Date: Mon, 25 Nov 2024 11:56:02 -0800 Subject: [PATCH] small changes to bring model into alignment with paper --- src/integrated/land.jl | 14 +++++-------- src/integrated/soil_canopy_model.jl | 12 ++++++++--- src/integrated/soil_snow_model.jl | 18 ++++++---------- src/standalone/Snow/snow_parameterizations.jl | 3 ++- src/standalone/Soil/boundary_conditions.jl | 12 ++++++++--- src/standalone/Soil/energy_hydrology.jl | 2 +- .../Soil/soil_heat_parameterizations.jl | 21 ++++++++----------- test/standalone/Snow/parameterizations.jl | 3 ++- test/standalone/Soil/soil_bc.jl | 2 +- .../standalone/Soil/soil_parameterizations.jl | 18 ++++------------ 10 files changed, 48 insertions(+), 57 deletions(-) diff --git a/src/integrated/land.jl b/src/integrated/land.jl index 36fec289e2..4cf155ed1c 100644 --- 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 index 20e26a45aa..1ac47fefbf 100644 --- 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 index c9ebdaa327..7c257cbafe 100644 --- 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 index 9670a1d4c9..7532339fee 100644 --- 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 index 5b58c2c876..ea21ff8577 100644 --- 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 index 929ddb3b57..206b985bf5 100644 --- 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 index dc6fb35427..274829ac96 100644 --- 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,7 +227,7 @@ 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, @@ -235,16 +235,13 @@ function kersten_number( ν_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 index f1355f56d7..8ab6aafb09 100644 --- 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 index 47c5f27f07..8bd0cdd134 100644 --- 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 index 096b8ee156..32cda80394 100644 --- a/test/standalone/Soil/soil_parameterizations.jl +++ b/test/standalone/Soil/soil_parameterizations.jl @@ -129,16 +129,17 @@ 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)) - @@ -146,17 +147,6 @@ for FT in (Float32, Float64) )^(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)