add and use new free drainage BC (#1161)

Author: kmdeck

docs/src/APIs/Soil.md

@@ -72,6 +72,7 @@ ClimaLand.Soil.HeatFluxBC
 ClimaLand.Soil.WaterFluxBC
 ClimaLand.Soil.TemperatureStateBC
 ClimaLand.Soil.FreeDrainage
+ClimaLand.Soil.EnergyWaterFreeDrainage
 ClimaLand.Soil.RichardsAtmosDrivenFluxBC
 ClimaLand.Soil.AtmosDrivenFluxBC
 ClimaLand.Soil.WaterHeatBC

experiments/long_runs/soil.jl

@@ -122,13 +122,8 @@ function setup_model(FT, start_date, domain, earth_param_set)
         (:soil,),
     )
     zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
-    boundary_conditions = (;
-        top = top_bc,
-        bottom = Soil.WaterHeatBC(;
-            water = Soil.FreeDrainage(),
-            heat = zero_flux,
-        ),
-    )
+    boundary_conditions =
+        (; top = top_bc, bottom = Soil.EnergyWaterFreeDrainage())
     soil = soil_model_type(;
         boundary_conditions = boundary_conditions,
         sources = sources,

src/integrated/land.jl

@@ -114,13 +114,8 @@ function LandModel{FT}(;
     )
 
     zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
-    boundary_conditions = (;
-        top = top_bc,
-        bottom = Soil.WaterHeatBC(;
-            water = Soil.FreeDrainage(),
-            heat = zero_flux,
-        ),
-    )
+    boundary_conditions =
+        (; top = top_bc, bottom = Soil.EnergyWaterFreeDrainage())
     soil = soil_model_type(;
         boundary_conditions = boundary_conditions,
         sources = sources,

src/integrated/soil_canopy_model.jl

@@ -89,13 +89,8 @@ function SoilCanopyModel{FT}(;
         prognostic_land_components,
     )
     zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
-    boundary_conditions = (;
-        top = top_bc,
-        bottom = Soil.WaterHeatBC(;
-            water = Soil.FreeDrainage(),
-            heat = zero_flux,
-        ),
-    )
+    boundary_conditions =
+        (; top = top_bc, bottom = Soil.EnergyWaterFreeDrainage())
     soil = soil_model_type(;
         boundary_conditions = boundary_conditions,
         sources = sources,

src/standalone/Soil/boundary_conditions.jl

@@ -10,6 +10,7 @@ using ClimaCore: Geometry
 export TemperatureStateBC,
     MoistureStateBC,
     FreeDrainage,
+    EnergyWaterFreeDrainage,
     HeatFluxBC,
     WaterFluxBC,
     AtmosDrivenFluxBC,
@@ -48,9 +49,14 @@ end
 
 """
     FreeDrainage <: AbstractWaterBC
+
 A concrete type of soil boundary condition, for use at
 the BottomBoundary only, where the flux is set to be
-`F = -K∇h = -K`.
+`F = -K∇h = -K`. 
+
+This is not tied to any boundary condition for the heat equation. 
+To account for the energy flux resulting from free drainage of liquid
+water, please see `EnergyWaterFreeDrainage`.
 """
 struct FreeDrainage <: AbstractWaterBC end
 
@@ -569,6 +575,39 @@ function WaterHeatBC(; water, heat)
     return WaterHeatBC{typeof(water), typeof(heat)}(water, heat)
 end
 
+"""
+    EnergyWaterFreeDrainage <: AbstractEnergyHydrologyBC
+
+A concrete type of soil boundary condition, for use at
+the BottomBoundary only, where the fluxes are set to be
+`F_liq = -K∇h = -K`, `F_energy = -K ρe_liq`.
+
+That is, this enforces that the free drainage boundary condition
+for liquid water is paired the the corresponding loss of energy
+that that entails.
+"""
+struct EnergyWaterFreeDrainage <: AbstractEnergyHydrologyBC end
+
+function soil_boundary_fluxes!(
+    bc::EnergyWaterFreeDrainage,
+    boundary::ClimaLand.BottomBoundary,
+    soil::EnergyHydrology,
+    Δz,
+    Y,
+    p,
+    t,
+)
+    FT = eltype(Δz)
+    K_c = Fields.level(p.soil.K, 1)
+    T_c = Fields.level(p.soil.T, 1)
+    @. p.soil.bottom_bc.water = -1 * K_c
+    @. p.soil.bottom_bc.heat =
+        -1 *
+        K_c *
+        volumetric_internal_energy_liq(T_c, soil.parameters.earth_param_set)
+    return nothing
+end
+
 """
     AtmosDrivenFluxBC{
         A <: AbstractAtmosphericDrivers,

test/standalone/Soil/soil_bc.jl

@@ -192,6 +192,83 @@ for FT in (Float32, Float64)
 
         @test abs(flux_expected - flux_int) < eps(FT)
     end
+
+    @testset "Test water+energy free drainage BC, FT = $FT" begin
+        earth_param_set = LP.LandParameters(FT)
+        ν = FT(0.495)
+        K_sat = FT(0.0443 / 3600 / 100) # m/s
+        S_s = FT(1e-3) #inverse meters
+        vg_n = FT(2.0)
+        vg_α = FT(2.6) # inverse meters
+        hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
+        θ_r = FT(0.1)
+        ν_ss_om = FT(0.0)
+        ν_ss_quartz = FT(1.0)
+        ν_ss_gravel = FT(0.0)
+
+        parameters = Soil.EnergyHydrologyParameters(
+            FT;
+            ν,
+            ν_ss_om,
+            ν_ss_quartz,
+            ν_ss_gravel,
+            hydrology_cm,
+            K_sat,
+            S_s,
+            θ_r,
+        )
+
+        zmax = FT(0)
+        zmin = FT(-1)
+        nelems = 20
+        domain = Column(; zlim = (zmin, zmax), nelements = nelems)
+
+        zero_water_flux_bc = WaterFluxBC((p, t) -> 0.0)
+        zero_heat_flux_bc = HeatFluxBC((p, t) -> 0.0)
+        boundary_conditions = (;
+            top = WaterHeatBC(;
+                water = zero_water_flux_bc,
+                heat = zero_heat_flux_bc,
+            ),
+            bottom = EnergyWaterFreeDrainage(),
+        )
+        sources = ()
+        soil = Soil.EnergyHydrology{FT}(;
+            parameters,
+            domain,
+            boundary_conditions,
+            sources,
+        )
+
+        Y, p, coords = initialize(soil)
+        Y.soil.ϑ_l .= ν / 2
+        Y.soil.θ_i .= 0
+        T = coords.subsurface.z .* 10 .+ FT(290.0)
+        ρc_s =
+            Soil.volumetric_heat_capacity.(
+                Y.soil.ϑ_l,
+                Y.soil.θ_i,
+                parameters.ρc_ds,
+                parameters.earth_param_set,
+            )
+        Y.soil.ρe_int .=
+            Soil.volumetric_internal_energy.(
+                Y.soil.θ_i,
+                ρc_s,
+                T,
+                parameters.earth_param_set,
+            )
+        set_initial_cache! = ClimaLand.make_set_initial_cache(soil)
+        set_initial_cache!(p, Y, 0.0)
+        flux_expected_water = -1 .* ClimaCore.Fields.level(p.soil.K, 1)
+        T1 = ClimaCore.Fields.level(p.soil.T, 1)
+        flux_expected_energy =
+            flux_expected_water .*
+            Soil.volumetric_internal_energy_liq.(T1, parameters.earth_param_set)
+        @test all(parent(p.soil.bottom_bc.heat) .≈ parent(flux_expected_energy))
+        @test all(parent(p.soil.bottom_bc.water) .≈ parent(flux_expected_water))
+
+    end
 end