add energy from rain to snow (#1176)

Author: kmdeck

.buildkite/Manifest-v1.11.toml

@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.5"
+julia_version = "1.11.4"
 manifest_format = "2.0"
 project_hash = "c057205b5dfa45d173ab89f45dff9c1d2ab4a915"
 
@@ -472,7 +472,7 @@ uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 version = "0.2.14"
 
 [[deps.ClimaLand]]
-deps = ["ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaTimeSteppers", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
+deps = ["ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaTimeSteppers", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LazyBroadcast", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.16.2"
@@ -2245,7 +2245,7 @@ version = "2.5.4+0"
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.5+0"
+version = "0.8.1+4"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]

Project.toml

@@ -14,6 +14,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+LazyBroadcast = "9dccce8e-a116-406d-9fcc-a88ed4f510c8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -23,8 +24,8 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [weakdeps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
-ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
@@ -51,6 +52,7 @@ HTTP = "1.10"
 Insolation = "0.9.2"
 Interpolations = "0.15.1, 0.16"
 LazyArtifacts = "1"
+LazyBroadcast = "1.0.0"
 LinearAlgebra = "1"
 NCDatasets = "0.14"
 SciMLBase = "2.53"

docs/src/APIs/Snow.md

@@ -23,6 +23,8 @@ ClimaLand.Snow.compute_water_runoff
 ClimaLand.Snow.energy_from_q_l_and_swe
 ClimaLand.Snow.energy_from_T_and_swe
 ClimaLand.Snow.snow_cover_fraction
+ClimaLand.Snow.volumetric_energy_flux_falling_rain
+ClimaLand.Snow.volumetric_energy_flux_falling_snow
 ```
 
 ## Computing fluxes for snow

src/ClimaLand.jl

@@ -2,6 +2,7 @@ module ClimaLand
 using DocStringExtensions
 
 using ClimaCore
+using LazyBroadcast: @lazy
 import ClimaCore: Fields, Spaces
 
 include("shared_utilities/Parameters.jl")

src/integrated/land.jl

@@ -570,20 +570,19 @@ function snow_boundary_fluxes!(
             p.snow.water_runoff
         ) * p.snow.snow_cover_fraction
 
-    # I think we want dU/dt to include energy of falling snow.
-    # otherwise snow can fall but energy wont change
-    # We are assuming that the sensible heat portion of snow is negligible.
-    _LH_f0 = FT(LP.LH_f0(model.parameters.earth_param_set))
-    _ρ_liq = FT(LP.ρ_cloud_liq(model.parameters.earth_param_set))
-    ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
+    ρe_flux_falling_snow =
+        Snow.volumetric_energy_flux_falling_snow(bc.atmos, p, model.parameters)
+    ρe_flux_falling_rain =
+        Snow.volumetric_energy_flux_falling_rain(bc.atmos, p, model.parameters)
 
     # positive fluxes are TOWARDS atmos, but R_n positive if snow absorbs energy
     @. p.snow.total_energy_flux =
-        P_snow * ρe_falling_snow +
+        ρe_flux_falling_snow +
         (
             p.snow.turbulent_fluxes.lhf +
             p.snow.turbulent_fluxes.shf +
-            p.snow.R_n - p.snow.energy_runoff - p.ground_heat_flux
+            p.snow.R_n - p.snow.energy_runoff - p.ground_heat_flux +
+            ρe_flux_falling_rain
         ) * p.snow.snow_cover_fraction
     return nothing
 end

src/integrated/soil_snow_model.jl

@@ -312,20 +312,19 @@ function snow_boundary_fluxes!(
             p.snow.water_runoff
         ) * p.snow.snow_cover_fraction
 
-    # I think we want dU/dt to include energy of falling snow.
-    # otherwise snow can fall but energy wont change
-    # We are assuming that the sensible heat portion of snow is negligible.
-    _LH_f0 = FT(LP.LH_f0(model.parameters.earth_param_set))
-    _ρ_liq = FT(LP.ρ_cloud_liq(model.parameters.earth_param_set))
-    ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
+    ρe_flux_falling_snow =
+        Snow.volumetric_energy_flux_falling_snow(bc.atmos, p, model.parameters)
+    ρe_flux_falling_rain =
+        Snow.volumetric_energy_flux_falling_rain(bc.atmos, p, model.parameters)
 
     # positive fluxes are TOWARDS atmos
     @. p.snow.total_energy_flux =
-        P_snow * ρe_falling_snow +
+        ρe_flux_falling_snow +
         (
             p.snow.turbulent_fluxes.lhf +
             p.snow.turbulent_fluxes.shf +
-            p.snow.R_n - p.snow.energy_runoff - p.ground_heat_flux
+            p.snow.R_n - p.snow.energy_runoff - p.ground_heat_flux +
+            ρe_flux_falling_rain
         ) * p.snow.snow_cover_fraction
 end
 

src/standalone/Snow/Snow.jl

@@ -3,6 +3,7 @@ module Snow
 using DocStringExtensions
 import ...Parameters as LP
 using ClimaCore
+using LazyBroadcast: @lazy
 using Thermodynamics
 using ClimaLand
 using ClimaLand:

src/standalone/Snow/boundary_fluxes.jl

@@ -102,10 +102,6 @@ function snow_boundary_fluxes!(
     p,
     t,
 ) where {FT}
-    parameters = model.parameters
-    _LH_f0 = FT(LP.LH_f0(parameters.earth_param_set))
-    _ρ_liq = FT(LP.ρ_cloud_liq(parameters.earth_param_set))
-
     bc = model.boundary_conditions
 
     turbulent_fluxes!(p.snow.turbulent_fluxes, bc.atmos, model, Y, p, t)
@@ -125,17 +121,18 @@ function snow_boundary_fluxes!(
             p.snow.water_runoff
         ) * p.snow.snow_cover_fraction
 
-    # We are assuming that the sensible heat portion of snow is negligible.
-    ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
+    ρe_flux_falling_snow =
+        volumetric_energy_flux_falling_snow(bc.atmos, p, model.parameters)
+    ρe_flux_falling_rain =
+        volumetric_energy_flux_falling_rain(bc.atmos, p, model.parameters)
 
     # positive fluxes are TOWARDS atmos
-    # This does not include the energy in rain yet
     @. p.snow.total_energy_flux =
-        P_snow * ρe_falling_snow +
+        ρe_flux_falling_snow +
         (
             p.snow.turbulent_fluxes.lhf +
             p.snow.turbulent_fluxes.shf +
-            p.snow.R_n - p.snow.energy_runoff
+            p.snow.R_n - p.snow.energy_runoff + ρe_flux_falling_rain
         ) * p.snow.snow_cover_fraction
     return nothing
 

src/standalone/Snow/snow_parameterizations.jl

@@ -11,7 +11,9 @@ export snow_surface_temperature,
     update_snow_cover_fraction!,
     snow_bulk_density,
     phase_change_flux,
-    update_snow_albedo!
+    update_snow_albedo!,
+    volumetric_energy_flux_falling_snow,
+    volumetric_energy_flux_falling_rain
 
 """
     update_snow_albedo!(α, m::ConstantAlbedoModel, Y, p, t, earth_param_set)
@@ -477,6 +479,45 @@ function energy_from_T_and_swe(S::FT, T::FT, parameters) where {FT}
 
 end
 
+"""
+    volumetric_energy_flux_falling_snow(atmos, p, parameters)
+
+Returns the volumetric energy flux of falling snow for a PrescribedAtmosphere,
+approximated as ρe_snow * P_snow, where ρe_snow = -LH_f0 * _ρ_liq. 
+This is a negative internal energy, due the to negative contribution of
+the latent heat of melting to the energy of the snow,
+ and it neglects the sensible heat portion of the snow. The overall flux
+ is per unit volume of liquid water, and P_snow is expressed as 
+the volume flux of liquid water resulting from the snow.
+
+This method can be extended to coupled simulations, where atmos is of type
+CoupledAtmosphere, and the energy flux of the falling snow is passed in the
+cache `p`. In that case, this should specify `atmos::PrescribedAtmosphere`.
+"""
+function volumetric_energy_flux_falling_snow(atmos, p, parameters)
+    _LH_f0 = LP.LH_f0(parameters.earth_param_set)
+    _ρ_liq = LP.ρ_cloud_liq(parameters.earth_param_set)
+    ρe_snow = -_LH_f0 * _ρ_liq
+    return @lazy @. ρe_snow * p.drivers.P_snow # per unit vol of liquid water
+end
+
+"""
+    volumetric_energy_flux_falling_rain(atmos, p, parameters)
+
+Returns the volumetric energy flux of falling rain for a PrescribedAtmosphere,
+approximated as ρ_l e_l(T_atmos) * P_liq. This is per unit volume of liquid water, 
+and P_liq is expressed as the volume flux of liquid water resulting from the rain.
+
+This method can be extended to coupled simulations, where atmos is of type
+CoupledAtmosphere, and the energy flux of the falling rain is passed in the
+cache `p`.  In that case, this should specify `atmos::PrescribedAtmosphere`.
+"""
+function volumetric_energy_flux_falling_rain(atmos, p, parameters)
+    return @lazy @. volumetric_internal_energy_liq(p.drivers.T, parameters) *
+                    p.drivers.P_liq
+end
+
+
 """
     update_density_and_depth!(ρ_snow, z_snow, density::MinimumDensityModel, Y, p, params::SnowParameters)