From 049fc42b4154136a80dc69002944211f22b41f09 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 28 Sep 2026 23:28:33 +0200 Subject: [PATCH 1/2] Fit ET not being applied to soil water tendency MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The coupled SoilEnergyWaterCarbon compute_tendencies! was not passing the evapotranspiration scheme to the Richards hydrology, so ET was never removed from the soil water in LandModel. Add an optional surface_hydrology argument, mirroring closure!, and pass it from LandModel. Also add a regression test checking that the top-layer saturation tendency drops by E / (Δz * porosity) when ground evaporation is applied. Fixes #203 --- src/models/coupled/land_model.jl | 3 ++- src/processes/soil/soil_coupled.jl | 10 ++++++++-- test/coupled_models/land_model_tests.jl | 21 +++++++++++++++++++++ 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/src/models/coupled/land_model.jl b/src/models/coupled/land_model.jl index 9d61c5a7c..c71d22d5e 100644 --- a/src/models/coupled/land_model.jl +++ b/src/models/coupled/land_model.jl @@ -140,7 +140,8 @@ end function compute_tendencies!(state, model::LandModel) grid = get_grid(model) compute_tendencies!(state, grid, model.surface_hydrology) - compute_tendencies!(state, grid, model.soil, model.constants) + # Pass the surface hydrology so that evapotranspiration is applied as a sink in the soil water tendency + compute_tendencies!(state, grid, model.soil, model.constants, model.surface_hydrology) # Snow tendencies after surface hydrology; no-op without snow (`nothing`) compute_tendencies!(state, grid, model.snow, model.constants, model.atmosphere) compute_tendencies!(state, grid, model.vegetation, model.constants, model.atmosphere) diff --git a/src/processes/soil/soil_coupled.jl b/src/processes/soil/soil_coupled.jl index 64a998558..6f97a07ad 100644 --- a/src/processes/soil/soil_coupled.jl +++ b/src/processes/soil/soil_coupled.jl @@ -87,14 +87,20 @@ end Compute tendencies for soil energy, water, and carbon state variables on `grid` based on the given values in `constants`. + +An optional `surface_hydrology` process may be supplied (by the coupled `LandModel`) so that its +evapotranspiration scheme (`get_evapotranspiration`) is applied as a sink term in the soil +water tendency. Without it (standalone soil), no evapotranspiration is removed from the soil. """ function compute_tendencies!( state, grid, soil::SoilEnergyWaterCarbon, - constants::PhysicalConstants + constants::PhysicalConstants, + surface_hydrology::Optional{AbstractSurfaceHydrology} = nothing ) + evapotranspiration = isnothing(surface_hydrology) ? nothing : get_evapotranspiration(surface_hydrology) # TODO: consider implementing fused kernel here? - compute_tendencies!(state, grid, soil.hydrology, soil, constants) + compute_tendencies!(state, grid, soil.hydrology, soil, constants, evapotranspiration) compute_tendencies!(state, grid, soil.biogeochem, soil, constants) compute_tendencies!(state, grid, soil.energy, soil, constants) return nothing diff --git a/test/coupled_models/land_model_tests.jl b/test/coupled_models/land_model_tests.jl index 548fb0771..7346c4daf 100644 --- a/test/coupled_models/land_model_tests.jl +++ b/test/coupled_models/land_model_tests.jl @@ -29,6 +29,27 @@ using Oceananigans.BoundaryConditions: BoundaryCondition, Flux energy_top_bc = integrator.state.internal_energy.boundary_conditions.top @test isa(energy_top_bc, BoundaryCondition{<:Flux}) @test energy_top_bc.condition == integrator.state.ground_heat_flux + # Check that ground evaporation is applied as a sink in the soil water tendency (issue #203). + # Compute the tendencies once with zero and once with nonzero ground evaporation; the difference + # in the top-layer saturation tendency must equal E / (Δz * porosity). The evaporation flux is set + # after the boundary conditions since the SEB solve recomputes it from the skin temperature. + state = integrator.state + function top_saturation_tendency(E) + Terrarium.reset_tendencies!(state) + compute_boundary_conditions!(state, land) + set!(state.evaporation_ground, E) + compute_tendencies!(state, land) + return Array(interior(state.tendencies.saturation_water_ice))[1, 1, end] + end + E = 2.0e-8 + dsat_no_ET = top_saturation_tendency(0.0) + dsat_ET = top_saturation_tendency(E) + Δz_top = Terrarium.Δzᵃᵃᶜ(1, 1, grid.Nz, grid) + strat = Terrarium.get_stratigraphy(soil) + bgc = Terrarium.get_biogeochemistry(soil) + por = Terrarium.porosity(1, 1, grid.Nz, grid, get_fields(state, strat, bgc), strat, bgc) + @test dsat_ET < dsat_no_ET + @test dsat_no_ET - dsat_ET ≈ E / (Δz_top * por) # Advance one timestep timestep!(integrator, 60.0) @test all(isfinite.(integrator.state.saturation_water_ice)) From a7c7e6d7fd39381c94bd1949c5515c520f5f40fd Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 29 Sep 2026 10:32:26 +0200 Subject: [PATCH 2/2] Add units in test comment --- test/coupled_models/land_model_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/coupled_models/land_model_tests.jl b/test/coupled_models/land_model_tests.jl index 7346c4daf..659488ab9 100644 --- a/test/coupled_models/land_model_tests.jl +++ b/test/coupled_models/land_model_tests.jl @@ -41,7 +41,7 @@ using Oceananigans.BoundaryConditions: BoundaryCondition, Flux compute_tendencies!(state, land) return Array(interior(state.tendencies.saturation_water_ice))[1, 1, end] end - E = 2.0e-8 + E = 2.0e-8 # m/s dsat_no_ET = top_saturation_tendency(0.0) dsat_ET = top_saturation_tendency(E) Δz_top = Terrarium.Δzᵃᵃᶜ(1, 1, grid.Nz, grid)