add energy flux of infiltration to soil (#1164)

Author: kmdeck

docs/src/APIs/Snow.md

@@ -23,8 +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
+ClimaLand.Snow.energy_flux_falling_rain
+ClimaLand.Snow.energy_flux_falling_snow
 ```
 
 ## Computing fluxes for snow

experiments/integrated/fluxnet/fluxnet_simulation.jl

@@ -13,7 +13,7 @@ timestepper = CTS.ARS111();
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 6,
+        max_iters = 3,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );

experiments/integrated/fluxnet/snow_soil/simulation.jl

@@ -253,27 +253,28 @@ _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set))
 _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set))
 fig = Figure(size = (1600, 1200), fontsize = 26)
 ax1 = Axis(fig[2, 1], ylabel = "ΔEnergy (J/A)", xlabel = "Days")
-function compute_surface_energy_fluxes(p)
-    ρe_flux_falling_snow =
-        Snow.volumetric_energy_flux_falling_snow(atmos, p, land.snow.parameters)
-    ρe_flux_falling_rain =
-        Snow.volumetric_energy_flux_falling_rain(atmos, p, land.snow.parameters)
+function compute_atmos_energy_fluxes(p)
+    e_flux_falling_snow =
+        Snow.energy_flux_falling_snow(atmos, p, land.snow.parameters)
+    e_flux_falling_rain =
+        Snow.energy_flux_falling_rain(atmos, p, land.snow.parameters)
 
     return @. (1 - p.snow.snow_cover_fraction) * (
                   p.soil.turbulent_fluxes.lhf +
                   p.soil.turbulent_fluxes.shf +
-                  p.soil.R_n
+                  p.soil.R_n +
+                  e_flux_falling_rain
               ) +
               p.snow.snow_cover_fraction * (
                   p.snow.turbulent_fluxes.lhf +
                   p.snow.turbulent_fluxes.shf +
                   p.snow.R_n +
-                  ρe_flux_falling_rain
+                  e_flux_falling_rain
               ) +
-              ρe_flux_falling_snow
+              e_flux_falling_snow
 end
 
-function compute_surface_water_vol_fluxes(p)
+function compute_atmos_water_vol_fluxes(p)
     return @. p.drivers.P_snow +
               p.drivers.P_liq +
               (1 - p.snow.snow_cover_fraction) * (
@@ -283,11 +284,37 @@ function compute_surface_water_vol_fluxes(p)
               p.snow.snow_cover_fraction * p.snow.turbulent_fluxes.vapor_flux
 end
 
+function compute_energy_of_runoff(p)
+    liquid_influx = @. p.snow.water_runoff * p.snow.snow_cover_fraction +
+       (1 - p.snow.snow_cover_fraction) * p.drivers.P_liq
+    e_flux_falling_rain =
+        Soil.volumetric_internal_energy_liq.(p.drivers.T, earth_param_set) .*
+        p.drivers.P_liq
+    influx_energy = @. p.snow.energy_runoff * p.snow.snow_cover_fraction +
+       (1 - p.snow.snow_cover_fraction) * e_flux_falling_rain
+    runoff_fraction = @. 1 - ClimaLand.Soil.compute_infiltration_fraction(
+        p.soil.infiltration,
+        liquid_influx,
+    )
+    return runoff_fraction .* influx_energy
+end
+
+function compute_runoff(p)
+    liquid_influx = @. p.snow.water_runoff * p.snow.snow_cover_fraction +
+       (1 - p.snow.snow_cover_fraction) * p.drivers.P_liq
+    runoff_fraction = @. 1 - ClimaLand.Soil.compute_infiltration_fraction(
+        p.soil.infiltration,
+        liquid_influx,
+    )
+    return runoff_fraction .* liquid_influx
+end
+
 ΔE_expected =
     cumsum(
         -1 .* [
             parent(
-                compute_surface_energy_fluxes(sv.saveval[k]) .-
+                compute_atmos_energy_fluxes(sv.saveval[k]) #.- compute_energy_of_runoff(sv.saveval[k])
+                .-
                 sv.saveval[k].soil.bottom_bc.heat,
             )[end] for k in 1:1:(length(sv.t) - 1)
         ],
@@ -300,8 +327,9 @@ E_measured = [
     cumsum(
         -1 .* [
             parent(
-                compute_surface_water_vol_fluxes(sv.saveval[k]) .-
-                sv.saveval[k].soil.bottom_bc.water .+ sv.saveval[k].soil.R_s,
+                compute_atmos_water_vol_fluxes(sv.saveval[k]) .-
+                compute_runoff(sv.saveval[k]) .-
+                sv.saveval[k].soil.bottom_bc.water,
             )[end] for k in 1:1:(length(sv.t) - 1)
         ],
     ) * (sv.t[2] - sv.t[1])
@@ -332,7 +360,6 @@ lines!(
 lines!(ax4, daily[2:end], ΔW_expected, label = "Expected")
 axislegend(ax4, position = :rt)
 
-
 ax3 = Axis(fig[2, 2], ylabel = "ΔE/E", xlabel = "Days", yscale = log10)
 lines!(
     ax3,

src/ClimaLand.jl

@@ -391,7 +391,11 @@ using .Pond
 import .Pond: surface_runoff
 include("standalone/Soil/Soil.jl")
 using .Soil
-import .Soil: soil_boundary_fluxes!, sublimation_source
+import .Soil:
+    soil_boundary_fluxes!,
+    sublimation_source,
+    compute_liquid_influx,
+    compute_infiltration_energy_flux
 import .Soil.Biogeochemistry: soil_temperature, soil_moisture
 include("standalone/Snow/Snow.jl")
 using .Snow

src/integrated/land.jl

@@ -479,38 +479,113 @@ the presence of the canopy modifies the soil BC.
 function soil_boundary_fluxes!(
     bc::AtmosDrivenFluxBC,
     prognostic_land_components::Val{(:canopy, :snow, :soil, :soilco2)},
-    soil::EnergyHydrology{FT},
+    soil::EnergyHydrology,
     Y,
     p,
     t,
-) where {FT}
+)
     turbulent_fluxes!(p.soil.turbulent_fluxes, bc.atmos, soil, Y, p, t)
-    # influx = maximum possible rate of infiltration given precip, snowmelt, evaporation/condensation
-    # but if this exceeds infiltration capacity of the soil, runoff will
-    # be generated.
-    # Use top_bc.water as temporary storage to avoid allocation
-    influx = p.soil.top_bc.water
-    @. influx =
-        p.snow.water_runoff * p.snow.snow_cover_fraction +
+    # Liquid influx is a combination of precipitation and snowmelt in general
+    liquid_influx =
+        Soil.compute_liquid_influx(p, soil, prognostic_land_components)
+    # This partitions the influx into runoff and infiltration
+    Soil.update_infiltration_water_flux!(
+        p,
+        bc.runoff,
+        liquid_influx,
+        Y,
+        t,
+        soil,
+    )
+    # This computes the energy of the infiltrating water
+    infiltration_energy_flux = Soil.compute_infiltration_energy_flux(
+        p,
+        bc.runoff,
+        bc.atmos,
+        prognostic_land_components,
+        liquid_influx,
+        soil,
+        Y,
+        t,
+    )
+    @. p.soil.top_bc.water =
+        p.soil.infiltration +
         p.excess_water_flux +
         (1 - p.snow.snow_cover_fraction) *
-        (p.drivers.P_liq + p.soil.turbulent_fluxes.vapor_flux_liq)
-    # The update_runoff! function computes how much actually infiltrates
-    # given influx and our runoff model bc.runoff, and updates
-    # p.soil.infiltration in place
-    Soil.Runoff.update_runoff!(p, bc.runoff, influx, Y, t, soil)
-    @. p.soil.top_bc.water = p.soil.infiltration
+        p.soil.turbulent_fluxes.vapor_flux_liq
+    # The actual boundary condition is a mix of liquid water infiltration and
+    # evaporation. The infiltration already has accounted for snow cover fraction,
+    # because the influx it is computed from has accounted for that.
+    # The last term, `excess water flux`, arises when snow melts in a timestep but
+    # has a nonzero sublimation which was applied for the entire step.
     @. p.soil.top_bc.heat =
         (1 - p.snow.snow_cover_fraction) * (
             p.soil.R_n +
             p.soil.turbulent_fluxes.lhf +
             p.soil.turbulent_fluxes.shf
         ) +
         p.excess_heat_flux +
-        p.snow.snow_cover_fraction * p.ground_heat_flux
+        p.snow.snow_cover_fraction * p.ground_heat_flux +
+        infiltration_energy_flux
+
     return nothing
 end
 
+"""
+   compute_liquid_influx(p,
+                         model,
+                         prognostic_land_components::Val{(:canopy, :snow, :soil, :soilco2)},
+    ) 
+
+Returns the liquid water volume flux at the surface of the soil; uses
+the same method as the soil+snow integrated model.
+"""
+function Soil.compute_liquid_influx(
+    p,
+    model,
+    prognostic_land_components::Val{(:canopy, :snow, :soil, :soilco2)},
+)
+    Soil.compute_liquid_influx(p, model, Val((:snow, :soil)))
+end
+
+"""
+    compute_infiltration_energy_flux(
+        p,
+        runoff,
+        atmos,
+        prognostic_land_components:::Val{(:canopy, :snow, :soil, :soilco2)},
+        liquid_influx,
+        model::EnergyHydrology,
+        Y,
+        t,
+    )
+
+Computes the energy associated with infiltration of
+liquid water into the soil; uses the same method as
+the soil+snow integrated model.
+"""
+function Soil.compute_infiltration_energy_flux(
+    p,
+    runoff,
+    atmos,
+    prognostic_land_components::Val{(:canopy, :snow, :soil, :soilco2)},
+    liquid_influx,
+    model::EnergyHydrology,
+    Y,
+    t,
+)
+    Soil.compute_infiltration_energy_flux(
+        p,
+        runoff,
+        atmos,
+        Val((:snow, :soil)),
+        liquid_influx,
+        model,
+        Y,
+        t,
+    )
+end
+
 function snow_boundary_fluxes!(
     bc::Snow.AtmosDrivenSnowBC,
     prognostic_land_components::Val{(:canopy, :snow, :soil, :soilco2)},
@@ -535,19 +610,19 @@ function snow_boundary_fluxes!(
             p.snow.water_runoff
         ) * p.snow.snow_cover_fraction
 
-    ρ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)
+    e_flux_falling_snow =
+        Snow.energy_flux_falling_snow(bc.atmos, p, model.parameters)
+    e_flux_falling_rain =
+        Snow.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 =
-        ρe_flux_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 +
-            ρe_flux_falling_rain
+            e_flux_falling_rain
         ) * p.snow.snow_cover_fraction
     return nothing
 end