From c4a70159fa2ef45069afa125fdabc32ba9458675 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 7 Apr 2026 17:09:04 +0200 Subject: [PATCH 01/27] Update speedy dry coupling example --- examples/Project.toml | 1 - examples/simulations/speedy_dry_land.jl | 58 ++++++++++++++----------- 2 files changed, 33 insertions(+), 26 deletions(-) diff --git a/examples/Project.toml b/examples/Project.toml index ee395dac8..0506e53c2 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -15,7 +15,6 @@ Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689" RingGrids = "d1845624-ad4f-453b-8ff4-a8db365bf3a7" SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9" -SpeedyWeatherInternals = "34489162-d270-4603-8b96-37b04f830d73" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Terrarium = "80418d68-07fa-499d-ae2b-c12e531f5cd8" diff --git a/examples/simulations/speedy_dry_land.jl b/examples/simulations/speedy_dry_land.jl index 871353059..8797a47b6 100644 --- a/examples/simulations/speedy_dry_land.jl +++ b/examples/simulations/speedy_dry_land.jl @@ -14,38 +14,52 @@ import SpeedyWeather as Speedy Speedy.SpeedyTransforms.FFTW.set_num_threads(1) """ -Naive implementation of a SpeedyWeather "dry" land model that couples to a Terrarium model. +Naive implementation of a SpeedyWeather "dry" land model based on Terrarium. """ -struct TerrariumDryLand{NF, TMI <: Terrarium.ModelIntegrator{NF}} <: Speedy.AbstractDryLand +struct TerrariumDryLand{ + NF, + LG <: Speedy.LandGeometry, + TM <: Terrarium.ModelIntegrator{NF}, + } <: Speedy.AbstractDryLand "Speedy spectral grid" spectral_grid::Speedy.SpectralGrid + "Speedy land model geometry" + geometry::LG + "Initialized Terrarium model integrator" - integrator::TMI + integrator::TM - function TerrariumDryLand(itnegrator::Terrarium.ModelIntegrator{NF, Arch, Grid}; spectral_grid_kwargs...) where {NF, Arch, Grid <: ColumnRingGrid} - spectral_grid = Speedy.SpectralGrid(integrator.model.grid.rings; NF, nlayers_soil = 1, spectral_grid_kwargs...) - return new{eltype(integrator), typeof(integrator)}(spectral_grid, integrator) + function TerrariumDryLand(integrator::Terrarium.ModelIntegrator{NF, Arch, Grid}; spectral_grid_kwargs...) where {NF, Arch, Grid <: ColumnRingGrid} + spectral_grid = Speedy.SpectralGrid(integrator.model.grid.rings; NF, spectral_grid_kwargs...) + land_grid = integrator.grid + Δz = on_architecture(CPU(), land_grid.z.Δᵃᵃᶜ) + geometry = Speedy.LandGeometry(1, Δz[end]) + return new{eltype(integrator), typeof(geometry), typeof(integrator)}(spectral_grid, geometry, integrator) end end +Speedy.variables(land::TerrariumDryLand) = ( + Speedy.PrognosticVariable(name = :soil_temperature, dims = Speedy.Grid3D(), namespace = :land), +) + function Speedy.initialize!( - progn_land::Speedy.PrognosticVariablesLand, # for dispatch - progn::Speedy.PrognosticVariables{NF}, - diagn::Speedy.DiagnosticVariables{NF}, - land::TerrariumDryLand, + ::Any, + progn::Speedy.PrognosticVariables, + diagn::Speedy.DiagnosticVariables, + land::TerrariumDryLand{NF}, model::Speedy.PrimitiveEquation, ) where {NF} - Tsoil = interior(land.integrator.state.temperature)[:, 1, end] - progn_land.soil_temperature .= Tsoil .+ NF(273.15) + Tsoil = interior(land.integrator.state.temperature)[:, 1, end] .+ 273.15 + progn.land.soil_temperature .= Tsoil Terrarium.initialize!(land.integrator) return nothing end function Speedy.timestep!( - progn::Speedy.PrognosticVariables{NF}, - diagn::Speedy.DiagnosticVariables{NF}, - land::TerrariumDryLand, + progn::Speedy.PrognosticVariables, + diagn::Speedy.DiagnosticVariables, + land::TerrariumDryLand{NF}, model::Speedy.PrimitiveEquation, ) where {NF} # get speedy state variables @@ -71,10 +85,8 @@ ring_grid = RingGrids.FullGaussianGrid(24) Nz = 30 Δz_min = 0.05 # currently the coupling is only stable with a large surface layer grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(; N = Nz, Δz_min), ring_grid) -# grid = ColumnGrid(CPU(), Float32, ExponentialSpacing(N=30)) # Initial conditions soil_initializer = SoilInitializer(eltype(grid)) - # Soil model with prescribed surface temperautre BC model = SoilModel(grid, initializer = soil_initializer) air_temperature = Field(grid, XY()) @@ -86,16 +98,15 @@ integrator = initialize(model, ForwardEuler(eltype(grid)), Tair_input, boundary_ land = TerrariumDryLand(integrator) # Set up coupled model land_sea_mask = Speedy.RockyPlanetMask(land.spectral_grid) -output = Speedy.NetCDFOutput(land.spectral_grid, Speedy.PrimitiveDryModel, path = "outputs/") +output = Speedy.NetCDFOutput(land.spectral_grid, Speedy.PrimitiveDryModel, land.geometry, path = "outputs/") time_stepping = Speedy.Leapfrog(land.spectral_grid, Δt_at_T31 = Minute(15)) primitive_dry_coupled = Speedy.PrimitiveDryModel(land.spectral_grid; land, land_sea_mask, time_stepping, output) # add soil temperature as output variable for Speedy simulation Speedy.add!(primitive_dry_coupled.output, Speedy.SoilTemperatureOutput()) # initialize coupled simulation sim_coupled = Speedy.initialize!(primitive_dry_coupled) -# Speedy.timestep!(sim) # run it -Speedy.run!(sim_coupled, period = Month(1), output = true) +Speedy.run!(sim_coupled, period = Day(2), output = true) # Soil temperature in the 5th layer (~0.54 m) Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 4], grid), title = "", size = (800, 400)) @@ -103,13 +114,10 @@ Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, Tair_fig = heatmap(sim_coupled.diagnostic_variables.grid.temp_grid[:, 8] .- 273.15, title = "Air temperature", size = (800, 400)) pres_fig = heatmap(exp.(sim_coupled.diagnostic_variables.grid.pres_grid), title = "Surface pressure", size = (800, 400)) srad_fig = heatmap(exp.(sim_coupled.diagnostic_variables.physics.surface_shortwave_down), title = "Surface shortwave down", size = (800, 400)) -# Tskin_fig = heatmap(RingGrids.Field(interior(integrator.state.skin_temperature)[:,1,end], grid), title="", size=(800,400)) -# save("plots/speedy_primitive_dry_coupled_tair_R48.png", Tair_fig, px_per_unit=1) -# save("plots/speedy_primitive_dry_coupled_tsoil_R48.png", Tsoil_fig, px_per_unit=1) # animate surface air and soil temperatures -Speedy.animate(sim, variable = "temp", coastlines = false, level = spectral_grid.nlayers, output_file = "plots/speedy_terrarium_dry_air_temperature.mp4") -Speedy.animate(sim, variable = "st", coastlines = false, level = 1, output_file = "plots/speedy_terrarium_dry_soil_temperature.mp4") +Speedy.animate(sim_coupled, variable = "temp", coastlines = false, level = spectral_grid.nlayers, output_file = "plots/speedy_terrarium_dry_air_temperature.mp4") +# Speedy.animate(sim_coupled, variable = "st", coastlines = false, level = 1, output_file = "plots/speedy_terrarium_dry_soil_temperature.mp4") # TODO: this is broken now for some reason? # pick a point somewhere in the mid-lattitudes T = interior(integrator.state.temperature)[2000, 1, :] From fbb22978cc0dba95dea046f316fdeb4fb115243d Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 21 Apr 2026 15:14:43 +0200 Subject: [PATCH 02/27] Add initial version of speedy wet coupling --- examples/simulations/speedy_dry_land.jl | 2 +- examples/simulations/speedy_wet_land.jl | 169 ++++++++++++++++++++++++ 2 files changed, 170 insertions(+), 1 deletion(-) create mode 100644 examples/simulations/speedy_wet_land.jl diff --git a/examples/simulations/speedy_dry_land.jl b/examples/simulations/speedy_dry_land.jl index 8797a47b6..a33831ef6 100644 --- a/examples/simulations/speedy_dry_land.jl +++ b/examples/simulations/speedy_dry_land.jl @@ -87,7 +87,7 @@ Nz = 30 grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(; N = Nz, Δz_min), ring_grid) # Initial conditions soil_initializer = SoilInitializer(eltype(grid)) -# Soil model with prescribed surface temperautre BC +# Soil model with prescribed surface temperature BC model = SoilModel(grid, initializer = soil_initializer) air_temperature = Field(grid, XY()) Tair_input = InputSource(grid, air_temperature; name = :air_temperature) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl new file mode 100644 index 000000000..0da2a46a3 --- /dev/null +++ b/examples/simulations/speedy_wet_land.jl @@ -0,0 +1,169 @@ +using Terrarium + +using CUDA +using Dates +using Rasters, NCDatasets +using Statistics + +using CairoMakie, GeoMakie + +import RingGrids +import SpeedyWeather as Speedy + +""" +Naive implementation of a SpeedyWeather "wet" land model based on Terrarium. +""" +struct TerrariumWetLand{ + NF, + LG <: Speedy.LandGeometry, + TM <: Terrarium.ModelIntegrator{NF}, + } <: Speedy.AbstractWetLand + "Speedy spectral grid" + spectral_grid::Speedy.SpectralGrid + + "Speedy land model geometry" + geometry::LG + + "Initialized Terrarium model integrator" + integrator::TM + + function TerrariumWetLand(integrator::Terrarium.ModelIntegrator{NF, Arch, Grid}; spectral_grid_kwargs...) where {NF, Arch, Grid <: ColumnRingGrid} + spectral_grid = Speedy.SpectralGrid(integrator.model.grid.rings; NF, spectral_grid_kwargs...) + land_grid = integrator.grid + Δz = on_architecture(CPU(), land_grid.z.Δᵃᵃᶜ) + geometry = Speedy.LandGeometry(1, Δz[end]) + return new{eltype(integrator), typeof(geometry), typeof(integrator)}(spectral_grid, geometry, integrator) + end +end + +Speedy.variables(land::TerrariumWetLand) = ( + Speedy.PrognosticVariable(name = :soil_temperature, dims = Speedy.Grid3D(), namespace = :land), +) + +function Speedy.initialize!( + ::Any, + progn::Speedy.PrognosticVariables, + diagn::Speedy.DiagnosticVariables, + land::TerrariumWetLand{NF}, + model::Speedy.PrimitiveEquation, + ) where {NF} + Terrarium.initialize!(land.integrator) + Tsoil = interior(land.integrator.state.temperature)[:, 1, end] .+ 273.15 + progn.land.soil_temperature .= Tsoil + return nothing +end + +function Speedy.timestep!( + progn::Speedy.PrognosticVariables, + diagn::Speedy.DiagnosticVariables, + land::TerrariumWetLand{NF}, + model::Speedy.PrimitiveEquation, + ) where {NF} + speedy_timestep!(progn, diagn, land) + return nothing +end + +function speedy_timestep!( + progn::Speedy.PrognosticVariables, + diagn::Speedy.DiagnosticVariables, + land::TerrariumWetLand{NF}, + ) where {NF} + # land constants + consts = land.integrator.model.constants + # get speedy state variables + Tair = @view diagn.grid.temp_grid[:, end] + humid = @view diagn.grid.humid_grid[:, end] + pres = diagn.grid.pres_grid + wind = diagn.physics.surface_wind_speed + rain = diagn.physics.rain_rate + snow = diagn.physics.snow_rate + SwIn = diagn.physics.surface_shortwave_down + LwIn = diagn.physics.surface_longwave_down + # update Terrarium input fields + state = land.integrator.state + set!(state.inputs.air_temperature, Tair) # first set directly to avoid allocating new fields + set!(state.inputs.air_temperature, state.inputs.air_temperature - 273.15) # then convert to celsius + set!(state.inputs.air_pressure, pres) # similar with pressure + set!(state.inputs.air_pressure, exp(state.inputs.air_pressure)) # take exp to get pressure in Pa + set!(state.inputs.specific_humidity, humid) + set!(state.inputs.rainfall, rain) + set!(state.inputs.snowfall, snow) + set!(state.inputs.windspeed, wind) + set!(state.inputs.surface_shortwave_down, SwIn) + set!(state.inputs.surface_longwave_down, LwIn) + # run land forward over speedy timestep interval; + # we use a smaller actual timestep to ensure stability + Terrarium.run!(land.integrator, period = progn.clock.Δt, Δt = 60.0) + # Update speedy variables + progn.land.soil_temperature .= state.skin_temperature .+ NF(273.15) + progn.land.sensible_heat_flux .= state.sensible_heat_flux + progn.land.surface_humidity_flux .= state.latent_heat_flux ./ (consts.Lsl * consts.ρₐ) + diagn.physics.surface_longwave_up .= state.surface_longwave_up + diagn.physics.surface_shortwave_up .= state.surface_shortwave_up + return nothing +end + +# quick test of default Speedy PrimitiveWetModel +ring_grid = RingGrids.FullGaussianGrid(24) +spectral_grid = Speedy.SpectralGrid(ring_grid) +primitive_wet = Speedy.PrimitiveWetModel(spectral_grid) +sim = Speedy.initialize!(primitive_wet) +Speedy.run!(sim, period = Day(1)) + +Nz = 30 +Δz_min = 0.05 # currently the coupling is only stable with a large surface layer +grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(; N = Nz, Δz_min), ring_grid) +# Initial conditions +soil_initializer = SoilInitializer(eltype(grid)) +# Land model with "prescribed" atmosphere (from the perspective of the land model at least...) +# vegetation = PrescribedVegetationCarbon(eltype(grid)) +model = LandModel(grid; initializer = soil_initializer, vegetation = nothing) +initializers = (;) +integrator = initialize(model, ForwardEuler(eltype(grid)); initializers) +# check if land model works standalone (with default atmospheric state) +timestep!(integrator, 60.0) # one step +run!(integrator, period = Hour(1)) # one hour +Terrarium.initialize!(integrator) + +# Initialize Terrarium-Speedy land model +land = TerrariumWetLand(integrator) +# Set up coupled model +land_sea_mask = Speedy.RockyPlanetMask(land.spectral_grid) +surface_heat_flux = Speedy.SurfaceHeatFlux(land.spectral_grid, land = Speedy.PrescribedLandHeatFlux()) +surface_humidity_flux = Speedy.SurfaceHumidityFlux(land.spectral_grid, land = Speedy.PrescribedLandHumidityFlux()) +output = Speedy.NetCDFOutput(land.spectral_grid, Speedy.PrimitiveDryModel, land.geometry, path = "outputs/") +time_stepping = Speedy.Leapfrog(land.spectral_grid, Δt_at_T31 = Minute(15)) +primitive_wet_coupled = Speedy.PrimitiveWetModel( + land.spectral_grid; + land, + surface_heat_flux, + surface_humidity_flux, + land_sea_mask, + time_stepping, + output +) +# add soil temperature as output variable for Speedy simulation +Speedy.add!(primitive_wet_coupled.output, Speedy.SoilTemperatureOutput()) +# initialize coupled simulation +sim_coupled = Speedy.initialize!(primitive_wet_coupled) +# run it +@run Speedy.run!(sim_coupled, steps = 1, output = true) + +# Soil temperature in the 5th layer (~0.54 m) +Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 4], grid), title = "", size = (800, 400)) +# Atmosphere variables +Tair_fig = heatmap(sim_coupled.diagnostic_variables.grid.temp_grid[:, 8] .- 273.15, title = "Air temperature", size = (800, 400)) +pres_fig = heatmap(exp.(sim_coupled.diagnostic_variables.grid.pres_grid), title = "Surface pressure", size = (800, 400)) +srad_fig = heatmap(exp.(sim_coupled.diagnostic_variables.physics.surface_shortwave_down), title = "Surface shortwave down", size = (800, 400)) + +# animate surface air and soil temperatures +Speedy.animate(sim_coupled, variable = "temp", coastlines = false, level = spectral_grid.nlayers, output_file = "plots/speedy_terrarium_dry_air_temperature.mp4") +# Speedy.animate(sim_coupled, variable = "st", coastlines = false, level = 1, output_file = "plots/speedy_terrarium_dry_soil_temperature.mp4") # TODO: this is broken now for some reason? + +# pick a point somewhere in the mid-lattitudes +T = interior(integrator.state.temperature)[2000, 1, :] +f = interior(integrator.state.liquid_water_fraction)[2000, 1, :] +zs = znodes(integrator.state.temperature) +# Plot temperature and liquid fraction profiles in upper 15 layers +Makie.scatterlines(T[(end - 15):end], zs[(end - 15):end]) +Makie.scatterlines(f, zs) From 196592fb55beeb9564c97c03147073996d2cd2a3 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 21 Apr 2026 17:25:47 +0200 Subject: [PATCH 03/27] Fix sign errors in ground heat flux --- .../surface_energy/skin_temperature.md | 2 +- .../surface_energy/surface_energy_balance.md | 4 +-- .../surface_energy/skin_temperature.jl | 2 +- test/surface_energy/skin_temperature.jl | 32 +++++++++++++++++-- 4 files changed, 34 insertions(+), 6 deletions(-) diff --git a/docs/src/processes/surface_energy/skin_temperature.md b/docs/src/processes/surface_energy/skin_temperature.md index 06536182e..8667299f1 100644 --- a/docs/src/processes/surface_energy/skin_temperature.md +++ b/docs/src/processes/surface_energy/skin_temperature.md @@ -21,7 +21,7 @@ ImplicitSkinTemperature The implicit approach requires solving a nonlinear equation at each grid point and time step. In order to avoid iteration, the current approach implementation of [`ImplicitSkinTemperature`](@ref) in Terrarium approximately solves for the skin temperature using a fixed point iteration where the ground heat flux is computed as the residual energy flux given the current prognostic state of the `skin_temperature` $T_s$, ```math -G^\star = R_{\text{net}}(T_s) - H_s(T_s) - H_l(T_s) +G^\star = R_{\text{net}}(T_s) + H_s(T_s) + H_l(T_s) ``` and the new skin temperature $T_s^\star$ is determined by setting this heat flux equal to gradient between the skin and the ground surface temperature (uppermost soil layer) $T_g$, ```math diff --git a/docs/src/processes/surface_energy/surface_energy_balance.md b/docs/src/processes/surface_energy/surface_energy_balance.md index 8ae2f8976..bb6fa35d9 100644 --- a/docs/src/processes/surface_energy/surface_energy_balance.md +++ b/docs/src/processes/surface_energy/surface_energy_balance.md @@ -16,10 +16,10 @@ The surface energy balance (SEB) describes how solar radiation, thermal radiatio ```math \begin{equation} -R_{\text{net}} = H_s + H_l + G\,. +R_{\text{net}} + H_s + H_l - G = 0\,. \end{equation} ``` -Following the standard convention of Terrarium and Oceananigans, all surface energy fluxes are defined **positive upward** (away from surface). +Following the standard convention of Terrarium and Oceananigans, all surface energy fluxes are defined **positive upward**. The negative sign in front of $G$ reflects that heat flows *towards* the surface from the uppermost ground layer. The [`SurfaceEnergyBalance`](@ref) process is responsible for computing all of the above flux terms and thus closing the energy balance between the atmosphere and land surface. Implementations of [`AbstractSurfaceEnergyBalance`](@ref) should generally include, at minimum, representations of each of the four SEB components: - An implementation of [`AbstractSkinTemperature`](@ref) that defines and updates both the skin temperature $T_s$ and the ground heat flux $G$. The skin temperature $T_s$ is the effective radiative temperature of the land surface. For an implicit approach, $T_s$ self-consistently satisfies the energy balance at each time step. For a prescribed approach, $T_s$ is given as input. The ground heat flux at the surface is derived either directly or as a residual from the energy balance. See [Skin temperature and ground heat flux](@ref) for further details. diff --git a/src/processes/surface_energy/skin_temperature.jl b/src/processes/surface_energy/skin_temperature.jl index 69e38f01b..b9db87a0b 100644 --- a/src/processes/surface_energy/skin_temperature.jl +++ b/src/processes/surface_energy/skin_temperature.jl @@ -75,7 +75,7 @@ sensible `H_s` and latent `H_l` heat flux. """ @inline function compute_ground_heat_flux(::AbstractSkinTemperature, R_net, H_s, H_l) # Compute ground heat flux as the residual of R_net and turbulent fluxes - G = R_net - H_s - H_l + G = R_net + H_s + H_l return G end diff --git a/test/surface_energy/skin_temperature.jl b/test/surface_energy/skin_temperature.jl index 79fbe1f20..515c97809 100644 --- a/test/surface_energy/skin_temperature.jl +++ b/test/surface_energy/skin_temperature.jl @@ -34,14 +34,42 @@ end # check that skin temperature converges for a large number of iterations Tskin_old = deepcopy(state.skin_temperature) resid = nothing - for i in 1:5 + balance = nothing + for i in 1:10 # compute fluxes Terrarium.compute_surface_energy_fluxes!(state, grid, seb, model.constants, model.atmosphere) # diagnose skin temperature Terrarium.update_skin_temperature!(state, grid, seb.skin_temperature) resid = maximum(abs.(state.skin_temperature - Tskin_old)) + balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux + state.ground_heat_flux Tskin_old = deepcopy(state.skin_temperature) - println("skin temperature residual at iteration $i: $resid") + println("skin temperature at iteration $i: $(state.skin_temperature[1, 1]) residual: $resid energy balance: $(balance[1, 1, 1])") end @test all(resid .< sqrt(eps())) + @test all(abs.(balance) .< sqrt(eps())) + # Cloudy and wet + set!(state.surface_shortwave_down, 150.0) + set!(state.surface_longwave_down, 50.0) + set!(state.specific_humidity, 0.01) + set!(state.air_pressure, 101_325) # standard pressure + set!(state.air_temperature, 10.0) # 10 °C + set!(state.ground_temperature, 2.0) # 2 °C + set!(state.windspeed, 1.0) # 10 m/s + set!(state.skin_temperature, 2.0) # initial skin temperature = ground temperature + compute_auxiliary!(state, grid, skin_temperature, seb) + @test all(isfinite.(state.skin_temperature)) + # check that skin temperature converges for a large number of iterations + Tskin_old = deepcopy(state.skin_temperature) + for i in 1:10 + # compute fluxes + Terrarium.compute_surface_energy_fluxes!(state, grid, seb, model.constants, model.atmosphere) + # diagnose skin temperature + Terrarium.update_skin_temperature!(state, grid, seb.skin_temperature) + resid = maximum(abs.(state.skin_temperature - Tskin_old)) + balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux + state.ground_heat_flux + Tskin_old = deepcopy(state.skin_temperature) + println("skin temperature at iteration $i: $(state.skin_temperature[1, 1]) residual: $resid energy balance: $(balance[1, 1, 1])") + end + @test all(resid .< sqrt(eps())) + @test all(abs.(balance) .< sqrt(eps())) end From dc6e33ebcb99f97170cbe6c45c11c1fcc768eda0 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 01:15:21 +0200 Subject: [PATCH 04/27] Fix bug in saturation initializer --- src/models/soil/soil_model_init.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/soil/soil_model_init.jl b/src/models/soil/soil_model_init.jl index 63e496413..e2318bb6f 100644 --- a/src/models/soil/soil_model_init.jl +++ b/src/models/soil/soil_model_init.jl @@ -147,6 +147,6 @@ end SaturationWaterTable(::Type{NF}; kwargs...) where {NF} = SaturationWaterTable{NF}(; kwargs...) function initialize!(state, ::AbstractModel, init::SaturationWaterTable{NF}) where {NF} - set!(state.saturation_water_ice, (x, z) -> z <= init.water_table_depth ? one(NF) : init.vadose_zone_saturation) + set!(state.saturation_water_ice, (x, z) -> z <= -init.water_table_depth ? one(NF) : init.vadose_zone_saturation) return nothing end From 5217bf0b67483cb11ef8166f502b3df2604f8b2a Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 01:30:41 +0200 Subject: [PATCH 05/27] Fix bug in soil moisture resistance factor Also use as default scheme for bare ground evaporation --- src/models/soil/soil_model_init.jl | 2 +- src/processes/soil/stratigraphy/soil_volume.jl | 2 ++ .../evapotranspiration/bare_ground_evaporation.jl | 2 +- .../evapotranspiration/evapotranspiration_base.jl | 4 ++-- .../evapotranspiration/ground_resistance_factor.jl | 3 ++- 5 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/models/soil/soil_model_init.jl b/src/models/soil/soil_model_init.jl index e2318bb6f..880edbb0a 100644 --- a/src/models/soil/soil_model_init.jl +++ b/src/models/soil/soil_model_init.jl @@ -140,7 +140,7 @@ Properties: $TYPEDFIELDS """ @kwdef struct SaturationWaterTable{NF} <: AbstractInitializer{NF} - vadose_zone_saturation::NF = 0.5 + vadose_zone_saturation::NF = 0.75 water_table_depth::NF = 5.0 end diff --git a/src/processes/soil/stratigraphy/soil_volume.jl b/src/processes/soil/stratigraphy/soil_volume.jl index 03f7c7e47..9237c014b 100644 --- a/src/processes/soil/stratigraphy/soil_volume.jl +++ b/src/processes/soil/stratigraphy/soil_volume.jl @@ -92,6 +92,8 @@ Alias for `SoilVolume{T, MineralOrganic{T}}` """ const MineralOrganicSoil{NF} = SoilVolume{NF, MineralOrganic{NF}} +@inline mineral_texture(solid::MineralOrganic) = solid.texture + @inline organic_fraction(solid::MineralOrganic) = solid.organic """ diff --git a/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl b/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl index da95ef0af..7835f9976 100644 --- a/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl +++ b/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl @@ -16,7 +16,7 @@ end BareGroundEvaporation( ::Type{NF}; - ground_resistance::GR = ConstantEvaporationResistanceFactor(one(NF)) + ground_resistance::GR = SoilMoistureResistanceFactor(NF) ) where {NF, GR} = BareGroundEvaporation{NF, GR}(ground_resistance) @propagate_inbounds surface_humidity_flux(i, j, grid, fields, evaporation::BareGroundEvaporation, args...) = fields.evaporation_ground[i, j] diff --git a/src/processes/surface_hydrology/evapotranspiration/evapotranspiration_base.jl b/src/processes/surface_hydrology/evapotranspiration/evapotranspiration_base.jl index 1683b1ab7..ac5da9037 100644 --- a/src/processes/surface_hydrology/evapotranspiration/evapotranspiration_base.jl +++ b/src/processes/surface_hydrology/evapotranspiration/evapotranspiration_base.jl @@ -7,8 +7,8 @@ Compute and return the evapotranspiration forcing for soil moisture at the given The ET forcing is just the `surface_humidity_flux` rescaled by the thickness of layer `k`. """ @inline function forcing(i, j, k, grid, clock, fields, evapotranspiration::AbstractEvapotranspiration, ::AbstractSoilHydrology) - let Δz = Δzᵃᵃᶜ(i, j, k, grid), - Qh = surface_humidity_flux(i, j, grid, fields, evapotranspiration) + let Δz = Δzᵃᵃᶜ(i, j, k, grid) + Qh = surface_humidity_flux(i, j, grid, fields, evapotranspiration) ∂θ∂t = -Qh / Δz # rescale by layer thickness to get water content flux return ∂θ∂t * (k == grid.Nz) end diff --git a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl index 8659c6c1f..9ae7da09b 100644 --- a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl +++ b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl @@ -46,7 +46,8 @@ ground_evaporation_resistance_factor(i, j, grid, fields, ::SoilMoistureResistanc hydrology = get_hydrology(soil) bgc = get_biogeochemistry(soil) soil = soil_volume(i, j, fgrid.Nz, grid, fields, strat, hydrology, bgc) - fc = field_capacity(get_hydraulic_properties(hydrology), soil) + texture = mineral_texture(soil) + fc = field_capacity(get_hydraulic_properties(hydrology), texture) fracs = volumetric_fractions(soil) if fracs.water < fc β = (1 - cos(π * fracs.water / fc))^2 / 4 From 5ed0e1effaf8f851f54371eb842ad5dc2d48bac5 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 01:31:07 +0200 Subject: [PATCH 06/27] Fix outdated method calls in surface energy model --- src/models/surface/surface_energy_model.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/models/surface/surface_energy_model.jl b/src/models/surface/surface_energy_model.jl index bdebdc24f..12974851f 100644 --- a/src/models/surface/surface_energy_model.jl +++ b/src/models/surface/surface_energy_model.jl @@ -41,12 +41,11 @@ function SurfaceEnergyModel( end function compute_auxiliary!(state, model::SurfaceEnergyModel) - compute_auxiliary!(state, model, model.atmosphere) - compute_auxiliary!(state, model, model.surface_energy_balance) + compute_auxiliary!(state, model.grid, model.atmosphere) + compute_auxiliary!(state, model.grid, model.surface_energy_balance, model.constants, model.atmosphere) return nothing end function compute_tendencies!(state, ::SurfaceEnergyModel) - compute_tendencies!(state, model, model.surface_energy_balance) return nothing end From f0b0310bb23c550cd3f19b62bb703770d8d1c968 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 15:22:58 +0200 Subject: [PATCH 07/27] Add sensible heat flux checks in skin temperature tests --- test/surface_energy/skin_temperature.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/test/surface_energy/skin_temperature.jl b/test/surface_energy/skin_temperature.jl index 515c97809..4822273b5 100644 --- a/test/surface_energy/skin_temperature.jl +++ b/test/surface_energy/skin_temperature.jl @@ -27,10 +27,11 @@ end set!(state.specific_humidity, 0.002) # dry conditions set!(state.air_pressure, 101_325) # standard pressure set!(state.air_temperature, 10.0) # 10 °C - set!(state.ground_temperature, 2.0) # 2 °C + set!(state.ground_temperature, 5.0) # 5 °C set!(state.windspeed, 1.0) # 1 m/s - compute_auxiliary!(state, grid, skin_temperature, seb) + compute_auxiliary!(state, model) @test all(isfinite.(state.skin_temperature)) + @test all(state.sensible_heat_flux .< 0) # check that skin temperature converges for a large number of iterations Tskin_old = deepcopy(state.skin_temperature) resid = nothing @@ -41,7 +42,7 @@ end # diagnose skin temperature Terrarium.update_skin_temperature!(state, grid, seb.skin_temperature) resid = maximum(abs.(state.skin_temperature - Tskin_old)) - balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux + state.ground_heat_flux + balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux - state.ground_heat_flux Tskin_old = deepcopy(state.skin_temperature) println("skin temperature at iteration $i: $(state.skin_temperature[1, 1]) residual: $resid energy balance: $(balance[1, 1, 1])") end @@ -52,12 +53,13 @@ end set!(state.surface_longwave_down, 50.0) set!(state.specific_humidity, 0.01) set!(state.air_pressure, 101_325) # standard pressure - set!(state.air_temperature, 10.0) # 10 °C - set!(state.ground_temperature, 2.0) # 2 °C - set!(state.windspeed, 1.0) # 10 m/s - set!(state.skin_temperature, 2.0) # initial skin temperature = ground temperature - compute_auxiliary!(state, grid, skin_temperature, seb) + set!(state.air_temperature, 5.0) # 5 °C + set!(state.ground_temperature, 10.0) # 10 °C + set!(state.windspeed, 10.0) # 10 m/s + set!(state.skin_temperature, 10.0) # initial skin temperature = ground temperature + compute_auxiliary!(state, model) @test all(isfinite.(state.skin_temperature)) + @test all(state.sensible_heat_flux .> 0) # check that skin temperature converges for a large number of iterations Tskin_old = deepcopy(state.skin_temperature) for i in 1:10 @@ -66,7 +68,7 @@ end # diagnose skin temperature Terrarium.update_skin_temperature!(state, grid, seb.skin_temperature) resid = maximum(abs.(state.skin_temperature - Tskin_old)) - balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux + state.ground_heat_flux + balance = state.surface_net_radiation + state.latent_heat_flux + state.sensible_heat_flux - state.ground_heat_flux Tskin_old = deepcopy(state.skin_temperature) println("skin temperature at iteration $i: $(state.skin_temperature[1, 1]) residual: $resid energy balance: $(balance[1, 1, 1])") end From a71ed8eff50891443bd5fee6c2d89659da07ebbc Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 15:24:12 +0200 Subject: [PATCH 08/27] Minor formatting tweaks --- .../surface_energy/turbulent_fluxes.jl | 52 ++++++------------- .../bare_ground_evaporation.jl | 3 +- .../ground_resistance_factor.jl | 7 +-- 3 files changed, 21 insertions(+), 41 deletions(-) diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl index 45c1411fa..e49a0cba2 100644 --- a/src/processes/surface_energy/turbulent_fluxes.jl +++ b/src/processes/surface_energy/turbulent_fluxes.jl @@ -89,12 +89,12 @@ Compute the sensible heat flux at `i, j` based on the current skin temperature a constants::PhysicalConstants, atmos::AbstractAtmosphere ) - let ρₐ = constants.ρₐ, # density of air - cₐ = constants.cₐ, # specific heat capacity of air - rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance - Tₛ = skin_temperature(i, j, grid, fields, skinT), # skin temperature - Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature - Q_T = (Tₛ - Tₐ) / rₐ # bulk aerodynamic temperature-gradient + let ρₐ = constants.ρₐ # density of air + cₐ = constants.cₐ # specific heat capacity of air + rₐ = aerodynamic_resistance(i, j, grid, fields, atmos) # aerodynamic resistance + Tₛ = skin_temperature(i, j, grid, fields, skinT) # skin temperature + Tₐ = air_temperature(i, j, grid, fields, atmos) # air temperature + Q_T = (Tₛ - Tₐ) / rₐ # bulk aerodynamic temperature-gradient # Calculate sensible heat flux (positive upwards) Hₛ = compute_sensible_heat_flux(tur, Q_T, ρₐ, cₐ) return Hₛ @@ -115,12 +115,12 @@ to the latent heat flux. constants::PhysicalConstants, atmos::AbstractAtmosphere ) - let L = constants.Llg, # specific latent heat of vaporization of water - ρₐ = constants.ρₐ, # density of air - Tₛ = skin_temperature(i, j, grid, fields, skinT), - rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance - Δq = compute_humidity_vpd(i, j, grid, fields, atmos, constants, Tₛ), - Q_h = Δq / rₐ # humidity flux + let L = constants.Llg # specific latent heat of vaporization of water + ρₐ = constants.ρₐ # density of air + Tₛ = skin_temperature(i, j, grid, fields, skinT) + rₐ = aerodynamic_resistance(i, j, grid, fields, atmos) # aerodynamic resistance + Δq = compute_humidity_vpd(i, j, grid, fields, atmos, constants, Tₛ) + Q_h = Δq / rₐ # humidity flux # Calculate latent heat flux (positive upwards) Hₗ = compute_latent_heat_flux(tur, Q_h, ρₐ, L) return Hₗ @@ -151,32 +151,10 @@ end # Kernels -@kernel function compute_auxiliary_kernel!( - out, grid, fields, - tur::DiagnosedTurbulentFluxes, - skinT::AbstractSkinTemperature, - constants::PhysicalConstants, - atmos::AbstractAtmosphere - ) - i, j = @index(Global, NTuple) - # compute sensible heat flux - out.sensible_heat_flux[i, j, 1] = compute_sensible_heat_flux(i, j, grid, fields, tur, skinT, constants, atmos) - # compute latent heat flux - out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, tur, skinT, constants, atmos) -end - -# TODO: Can these dispatches be standardized to reduce redundancy? -@kernel function compute_auxiliary_kernel!( - out, grid, fields, - tur::DiagnosedTurbulentFluxes, - skinT::AbstractSkinTemperature, - constants::PhysicalConstants, - atmos::AbstractAtmosphere, - evtr::AbstractEvapotranspiration - ) +@kernel function compute_auxiliary_kernel!(out, grid, fields, tur::DiagnosedTurbulentFluxes, args...) i, j = @index(Global, NTuple) # compute sensible heat flux - out.sensible_heat_flux[i, j, 1] = compute_sensible_heat_flux(i, j, grid, fields, tur, skinT, constants, atmos) + out.sensible_heat_flux[i, j, 1] = compute_sensible_heat_flux(i, j, grid, fields, tur, args...) # compute latent heat flux - out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, tur, evtr, constants) + out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, tur, args...) end diff --git a/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl b/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl index 7835f9976..91fdaaaac 100644 --- a/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl +++ b/src/processes/surface_hydrology/evapotranspiration/bare_ground_evaporation.jl @@ -58,5 +58,6 @@ end β = ground_evaporation_resistance_factor(i, j, grid, fields, evaporation.ground_resistance, soil) Δq = compute_humidity_vpd(i, j, grid, fields, atmos, constants, Ts) # Calculate water evaporation flux in m/s (positive upwards) - return out.evaporation_ground[i, j, 1] = β * Δq / rₐ + E_g = β * Δq / rₐ + return out.evaporation_ground[i, j, 1] = E_g end diff --git a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl index 9ae7da09b..eaaedf4ed 100644 --- a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl +++ b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl @@ -47,10 +47,11 @@ ground_evaporation_resistance_factor(i, j, grid, fields, ::SoilMoistureResistanc bgc = get_biogeochemistry(soil) soil = soil_volume(i, j, fgrid.Nz, grid, fields, strat, hydrology, bgc) texture = mineral_texture(soil) - fc = field_capacity(get_hydraulic_properties(hydrology), texture) fracs = volumetric_fractions(soil) - if fracs.water < fc - β = (1 - cos(π * fracs.water / fc))^2 / 4 + θfc = field_capacity(get_hydraulic_properties(hydrology), texture) + θw = fracs.water + if fracs.water < θfc + β = (1 - cos(π * θw / θfc))^2 / 4 else β = NF(1) end From 857cb7ac71a1102cc1dd4a711c76fb9c84b59b0f Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 15:24:40 +0200 Subject: [PATCH 09/27] Fix incorrect flux in speedy wet coupling --- examples/simulations/speedy_wet_land.jl | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl index 0da2a46a3..a306da194 100644 --- a/examples/simulations/speedy_wet_land.jl +++ b/examples/simulations/speedy_wet_land.jl @@ -97,7 +97,7 @@ function speedy_timestep!( # Update speedy variables progn.land.soil_temperature .= state.skin_temperature .+ NF(273.15) progn.land.sensible_heat_flux .= state.sensible_heat_flux - progn.land.surface_humidity_flux .= state.latent_heat_flux ./ (consts.Lsl * consts.ρₐ) + progn.land.surface_humidity_flux .= state.latent_heat_flux ./ consts.Lsl diagn.physics.surface_longwave_up .= state.surface_longwave_up diagn.physics.surface_shortwave_up .= state.surface_shortwave_up return nothing @@ -115,15 +115,16 @@ Nz = 30 grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(; N = Nz, Δz_min), ring_grid) # Initial conditions soil_initializer = SoilInitializer(eltype(grid)) +soil = SoilEnergyWaterCarbon(eltype(grid), hydrology = SoilHydrology(eltype(grid), RichardsEq())) # Land model with "prescribed" atmosphere (from the perspective of the land model at least...) # vegetation = PrescribedVegetationCarbon(eltype(grid)) -model = LandModel(grid; initializer = soil_initializer, vegetation = nothing) +model = LandModel(grid; initializer = soil_initializer, vegetation = nothing, soil) initializers = (;) integrator = initialize(model, ForwardEuler(eltype(grid)); initializers) # check if land model works standalone (with default atmospheric state) timestep!(integrator, 60.0) # one step -run!(integrator, period = Hour(1)) # one hour -Terrarium.initialize!(integrator) +run!(integrator, period = Hour(1), Δt = 120.0) # one hour +Terrarium.initialize!(integrator) # reinitialize before setting up atmosphere # Initialize Terrarium-Speedy land model land = TerrariumWetLand(integrator) @@ -147,10 +148,13 @@ Speedy.add!(primitive_wet_coupled.output, Speedy.SoilTemperatureOutput()) # initialize coupled simulation sim_coupled = Speedy.initialize!(primitive_wet_coupled) # run it -@run Speedy.run!(sim_coupled, steps = 1, output = true) +Speedy.run!(sim_coupled, period = Hour(1)) -# Soil temperature in the 5th layer (~0.54 m) -Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 4], grid), title = "", size = (800, 400)) +# Soil temperature in the 3rd layer +Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 2], grid), title = "", size = (800, 400)) +Tsurf_fig = heatmap(RingGrids.Field(interior(integrator.state.skin_temperature)[:, 1], grid), title = "", size = (800, 400)) +sat_fig = heatmap(RingGrids.Field(interior(integrator.state.saturation_water_ice)[:, 1, end], grid), title = "", size = (800, 400)) +E_fig = heatmap(RingGrids.Field(interior(integrator.state.evaporation_ground)[:, 1], grid), title = "", size = (800, 400)) # Atmosphere variables Tair_fig = heatmap(sim_coupled.diagnostic_variables.grid.temp_grid[:, 8] .- 273.15, title = "Air temperature", size = (800, 400)) pres_fig = heatmap(exp.(sim_coupled.diagnostic_variables.grid.pres_grid), title = "Surface pressure", size = (800, 400)) @@ -162,8 +166,10 @@ Speedy.animate(sim_coupled, variable = "temp", coastlines = false, level = spect # pick a point somewhere in the mid-lattitudes T = interior(integrator.state.temperature)[2000, 1, :] +sat = interior(integrator.state.saturation_water_ice)[2000, 1, :] f = interior(integrator.state.liquid_water_fraction)[2000, 1, :] zs = znodes(integrator.state.temperature) # Plot temperature and liquid fraction profiles in upper 15 layers Makie.scatterlines(T[(end - 15):end], zs[(end - 15):end]) +Makie.scatterlines(sat[(end - 15):end], zs[(end - 15):end]) Makie.scatterlines(f, zs) From d131c8b9a4cf7b6137904f9e6a9aae1c34f8f266 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 19:25:20 +0200 Subject: [PATCH 10/27] Enforce minimum saturation level in soil hydrology --- .../soil/hydrology/soil_hydraulic_closures.jl | 8 +++- .../hydrology/soil_hydraulic_properties.jl | 45 ++++++++++++------- .../soil/hydrology/soil_hydrology.jl | 12 ++--- .../soil/stratigraphy/soil_volume.jl | 14 +++--- 4 files changed, 51 insertions(+), 28 deletions(-) diff --git a/src/processes/soil/hydrology/soil_hydraulic_closures.jl b/src/processes/soil/hydrology/soil_hydraulic_closures.jl index 9ec5b19ad..58220fa1c 100644 --- a/src/processes/soil/hydrology/soil_hydraulic_closures.jl +++ b/src/processes/soil/hydrology/soil_hydraulic_closures.jl @@ -94,7 +94,9 @@ end ψm = ψ - ψh - ψz swrc = get_swrc(hydrology) por = porosity(i, j, k, grid, fields, strat, bgc) - vol_water_ice_content = swrc(ψm; θsat = por) + sat_res = residual_saturation(get_hydraulic_properties(hydrology)) + θres = sat_res * por + vol_water_ice_content = swrc(ψm; θsat = por, θres) saturation_water_ice[i, j, k] = vol_water_ice_content / por return nothing end @@ -114,8 +116,10 @@ end # get inverse of SWRC inv_swrc = inv(get_swrc(hydrology)) por = porosity(i, j, k, grid, fields, strat, bgc) + sat_res = residual_saturation(get_hydraulic_properties(hydrology)) + θres = sat_res * por # compute matric pressure head - ψm = inv_swrc(sat * por; θsat = por) + ψm = inv_swrc(sat * por; θsat = por, θres) # compute elevation pressure head ψz = z - z_ref # compute hydrostatic pressure head assuming impermeable lower boundary diff --git a/src/processes/soil/hydrology/soil_hydraulic_properties.jl b/src/processes/soil/hydrology/soil_hydraulic_properties.jl index b47cfc78f..09e54ded1 100644 --- a/src/processes/soil/hydrology/soil_hydraulic_properties.jl +++ b/src/processes/soil/hydrology/soil_hydraulic_properties.jl @@ -51,6 +51,13 @@ Compute the empirical field capacity of the soil. """ function field_capacity end +""" + $SIGNATURES + +Compute the (numerical) residual saturation level of the soil. +""" +function residual_saturation end + # Hydraulics parameterizations """ @@ -73,11 +80,14 @@ $TYPEDFIELDS "Hydraulic conductivity at saturation [m/s]" sat_hydraulic_cond::NF = 1.0e-5 - "Prescribed field capacity [-]" + "Constant field capacity [-]" field_capacity::NF = 0.25 - "Prescribed wilting point [-]" + "Constant wilting point [-]" wilting_point::NF = 0.05 + + "Residual (minimum) saturation level [-]" + residual::NF = 1.0e-3 end function ConstantSoilHydraulics( @@ -96,6 +106,8 @@ end @inline field_capacity(hydraulics::ConstantSoilHydraulics, args...) = hydraulics.field_capacity +@inline residual_saturation(hydraulics::ConstantSoilHydraulics, args...) = hydraulics.residual + """ $TYPEDEF @@ -127,6 +139,9 @@ $TYPEDFIELDS "Exponent of field capacity adjustment due to clay content [-]" field_capacity_exp::NF = 0.35 + + "Residual (minimum) saturation level [-]" + residual::NF = 1.0e-3 end function SoilHydraulicsSURFEX( @@ -142,6 +157,8 @@ end # TODO: this is not quite correct, SURFEX uses a hydraulic conductivity function that decreases exponentially with depth @inline saturated_hydraulic_conductivity(hydraulics::SoilHydraulicsSURFEX, args...) = hydraulics.sat_hydraulic_cond +@inline residual_saturation(hydraulics::SoilHydraulicsSURFEX, args...) = hydraulics.residual + @inline function wilting_point(hydraulics::SoilHydraulicsSURFEX, texture::SoilTexture) β_w = hydraulics.wilting_point_coef wp = β_w * sqrt(texture.clay * 100) @@ -171,10 +188,10 @@ function hydraulic_conductivity( hydraulics::AbstractSoilHydraulics{NF, RC, UnsatKLinear{NF}}, soil::SoilVolume ) where {NF, RC} - let fracs = volumetric_fractions(soil), - θw = fracs.water, # unfrozen water content - θsat = fracs.water + fracs.ice + fracs.air, # water + ice content at saturation (porosity) - K_sat = saturated_hydraulic_conductivity(hydraulics, soil.solid.texture) + let fracs = volumetric_fractions(soil) + θw = fracs.water # unfrozen water content + θsat = fracs.water + fracs.ice + fracs.air # water + ice content at saturation (porosity) + K_sat = saturated_hydraulic_conductivity(hydraulics, soil.solid.texture) K = K_sat * θw / θsat return K end @@ -205,15 +222,13 @@ function hydraulic_conductivity( soil::SoilVolume ) where {NF} # TODO: The SWRC parameters will need to also be spatially varying at some point - let n = hydraulics.swrc.n, # van Genuchten parameter `n` - fracs = volumetric_fractions(soil), - θw = fracs.water, # unfrozen water content - θwi = fracs.water + fracs.ice, # total water + ice content - θsat = porosity(soil), # porosity - f = liquid_fraction(soil), - Ω = hydraulics.unsat_hydraulic_cond.impedance, # scaling parameter for ice impedance - I_ice = 10^(-Ω * (1 - f)), # ice impedance factor - K_sat = saturated_hydraulic_conductivity(hydraulics, soil.solid.texture) + let n = hydraulics.swrc.n # van Genuchten parameter `n` + θw = water(soil) # volumetric content of unfrozen water + θsat = porosity(soil) # porosity + f = liquid_fraction(soil) # fraction of pore water that is unfrozen (equiv. θw / θwi) + Ω = hydraulics.unsat_hydraulic_cond.impedance # scaling parameter for ice impedance + I_ice = 10^(-Ω * (1 - f)) # ice impedance factor + K_sat = saturated_hydraulic_conductivity(hydraulics, soil.solid.texture) # We use `complex` types here to permit illegal state values which may occur when using adaptive time steppers. K = abs(K_sat * I_ice * sqrt(complex(θw / θsat)) * (1 - complex(1 - complex(θw / θsat)^(n / (n + 1)))^((n - 1) / n))^2) return K diff --git a/src/processes/soil/hydrology/soil_hydrology.jl b/src/processes/soil/hydrology/soil_hydrology.jl index b4709add7..883500f82 100644 --- a/src/processes/soil/hydrology/soil_hydrology.jl +++ b/src/processes/soil/hydrology/soil_hydrology.jl @@ -182,9 +182,11 @@ arising due to numerical error. This implementation scans over the saturation pr grid cell and redistributes excess water upward layer-by-layer until reaching the topmost layer, where any remaining excess water is added to the `surface_excess_water` pool. """ -@propagate_inbounds function adjust_saturation_profile!(out, i, j, grid, ::SoilHydrology{NF}) where {NF} +@propagate_inbounds function adjust_saturation_profile!(out, i, j, grid, hydrology::SoilHydrology{NF}) where {NF} sat = out.saturation_water_ice surface_excess_water = out.surface_excess_water + props = get_hydraulic_properties(hydrology) + sat_min = residual_saturation(props) field_grid = get_field_grid(grid) N = field_grid.Nz @@ -200,8 +202,8 @@ any remaining excess water is added to the `surface_excess_water` pool. # then from top to bottom for k in N:-1:2 - # calculate saturation deficit - deficit_sat = max(-sat[i, j, k], zero(NF)) + # calculate saturation deficit from residual saturation level + deficit_sat = max(-sat[i, j, k] + sat_min, zero(NF)) # add back saturation deficit and subtract from layer below sat[i, j, k] += deficit_sat sat[i, j, k - 1] -= deficit_sat * Δzᵃᵃᶜ(i, j, k, field_grid) / Δzᵃᵃᶜ(i, j, k - 1, field_grid) @@ -212,9 +214,9 @@ any remaining excess water is added to the `surface_excess_water` pool. sat[i, j, N] -= excess_sat surface_excess_water[i, j, 1] += excess_sat * Δzᵃᵃᶜ(i, j, N, field_grid) - # If the lowermost (bottom) layer has a deficit, just set to zero. + # If the lowermost (bottom) layer has a deficit, just set to the residual saturation level. # This constitutes a mass balance violation but should not happen under realistic conditions. - sat[i, j, 1] = max(sat[i, j, 1], zero(NF)) + sat[i, j, 1] = max(sat[i, j, 1], sat_min) return nothing end diff --git a/src/processes/soil/stratigraphy/soil_volume.jl b/src/processes/soil/stratigraphy/soil_volume.jl index 9237c014b..0590fd8a6 100644 --- a/src/processes/soil/stratigraphy/soil_volume.jl +++ b/src/processes/soil/stratigraphy/soil_volume.jl @@ -30,17 +30,19 @@ $FIELDS end end -porosity(soil::SoilVolume) = soil.porosity +@inline porosity(soil::SoilVolume) = soil.porosity -saturation(soil::SoilVolume) = soil.saturation +@inline saturation(soil::SoilVolume) = soil.saturation -liquid_fraction(soil::SoilVolume) = soil.liquid +@inline liquid_fraction(soil::SoilVolume) = soil.liquid -water_ice(soil::SoilVolume) = soil.porosity * soil.saturation +@inline water_ice(soil::SoilVolume) = soil.porosity * soil.saturation -organic_fraction(soil::SoilVolume) = organic_fraction(soil.solid) +@inline water(soil::SoilVolume) = soil.liquid * soil.porosity * soil.saturation -mineral_texture(soil::SoilVolume) = mineral_texture(soil.solid) +@inline organic_fraction(soil::SoilVolume) = organic_fraction(soil.solid) + +@inline mineral_texture(soil::SoilVolume) = mineral_texture(soil.solid) """ $TYPEDSIGNATURES From 80486048ff5eae68cc67a02e4011b7d8c62c0708 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 22 Apr 2026 22:39:19 +0200 Subject: [PATCH 11/27] Fix incorrect latent heat constant in coupling script --- examples/simulations/speedy_wet_land.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl index a306da194..3019ea99d 100644 --- a/examples/simulations/speedy_wet_land.jl +++ b/examples/simulations/speedy_wet_land.jl @@ -97,7 +97,7 @@ function speedy_timestep!( # Update speedy variables progn.land.soil_temperature .= state.skin_temperature .+ NF(273.15) progn.land.sensible_heat_flux .= state.sensible_heat_flux - progn.land.surface_humidity_flux .= state.latent_heat_flux ./ consts.Lsl + progn.land.surface_humidity_flux .= state.latent_heat_flux ./ consts.Llg diagn.physics.surface_longwave_up .= state.surface_longwave_up diagn.physics.surface_shortwave_up .= state.surface_shortwave_up return nothing @@ -148,7 +148,7 @@ Speedy.add!(primitive_wet_coupled.output, Speedy.SoilTemperatureOutput()) # initialize coupled simulation sim_coupled = Speedy.initialize!(primitive_wet_coupled) # run it -Speedy.run!(sim_coupled, period = Hour(1)) +Speedy.run!(sim_coupled, period = Hour(4)) # Soil temperature in the 3rd layer Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 2], grid), title = "", size = (800, 400)) From a25950f5c45a0f955a44980aac46ede991bc9b27 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Thu, 23 Apr 2026 00:14:42 +0200 Subject: [PATCH 12/27] Add more debug hooks and rename nancheck! to checkfinite! --- src/diagnostics/debugging.jl | 16 ++++++++-------- src/grids/grid_utils.jl | 3 ++- src/processes/processes.jl | 4 ++-- src/timesteppers/abstract_timestepper.jl | 7 +++++++ 4 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/diagnostics/debugging.jl b/src/diagnostics/debugging.jl index 62d3b9bc1..86e21e64e 100644 --- a/src/diagnostics/debugging.jl +++ b/src/diagnostics/debugging.jl @@ -14,25 +14,25 @@ end """ $SIGNATURES -Check whether the given `field` has any `NaN` values using `Diagnostics.hasnan` and raise an error if `NaN`s are detected. +Check whether the given `field` has any `NaN` or `Inf` values and raise an error if `NaN`s are detected. """ -nancheck!(field::AbstractField, name = nothing) = Diagnostics.hasnan(field) && error("Found NaNs in Field $name: $field") -function nancheck!(nt::NamedTuple) +checkfinite!(field::AbstractField, name = nothing) = any(!isfinite, parent(field)) && error("Found NaN/Inf values in Field $name: $field") +function checkfinite!(nt::NamedTuple) for key in keys(nt) - nancheck!(nt[key], key) + checkfinite!(nt[key], key) end - return + return nothing end """ $SIGNATURES Provides a "hook" for handling debug calls from relevant callsites. Default implementations for -`Field` and `NamedTuple` (assumed to be of `Field`s) simply forward to [`nancheck!`](@ref). +`Field` and `NamedTuple` (assumed to be of `Field`s) simply forward to [`checkfinite!`](@ref). """ @inline debughook!(args...) = nothing -@inline debughook!(field::AbstractField) = nancheck!(field) -@inline debughook!(nt::NamedTuple) = nancheck!(nt) +@inline debughook!(field::AbstractField) = checkfinite!(field) +@inline debughook!(nt::NamedTuple) = checkfinite!(nt) """ $SIGNATURES diff --git a/src/grids/grid_utils.jl b/src/grids/grid_utils.jl index ab6fd8c96..50f3ea535 100644 --- a/src/grids/grid_utils.jl +++ b/src/grids/grid_utils.jl @@ -2,7 +2,8 @@ function Oceananigans.launch!(grid::AbstractLandGrid, workspec, kernel!::Function, first_arg, other_args...; kwargs...) fgrid = get_field_grid(grid) launch!(fgrid.architecture, fgrid, get_workspec(workspec), kernel!, first_arg, grid, other_args...; kwargs...) - return debugsite!(kernel!, first_arg, other_args...) + debugsite!(kernel!, first_arg, grid, other_args...) + return nothing end """ diff --git a/src/processes/processes.jl b/src/processes/processes.jl index 21fb4770a..88111abb5 100644 --- a/src/processes/processes.jl +++ b/src/processes/processes.jl @@ -124,5 +124,5 @@ export SurfaceHydrology include("surface_hydrology/surface_hydrology.jl") # Default debug hooks -@inline debughook!(::typeof(compute_auxiliary_kernel!), out, args...) = nancheck!(out) -@inline debughook!(::typeof(compute_tendencies_kernel!), out, args...) = nancheck!(out) +@inline debughook!(::typeof(compute_auxiliary_kernel!), out, args...) = checkfinite!(out) +@inline debughook!(::typeof(compute_tendencies_kernel!), out, args...) = checkfinite!(out) diff --git a/src/timesteppers/abstract_timestepper.jl b/src/timesteppers/abstract_timestepper.jl index 6613bacea..c1d11058f 100644 --- a/src/timesteppers/abstract_timestepper.jl +++ b/src/timesteppers/abstract_timestepper.jl @@ -67,8 +67,12 @@ function explicit_step!(state, grid::AbstractLandGrid, timestepper::AbstractTime fastiterate(keys(state.prognostic)) do name # apply flux BCs, if present compute_z_bcs!(state.tendencies[name], state.prognostic[name], grid, state) + # debug site post-BC + debugsite!(explicit_step!, state.tendencies[name], name) # update prognostic state variable explicit_step!(state.prognostic[name], state.tendencies[name], grid, timestepper, Δt) + # debug site post-step + debugsite!(explicit_step!, state.prognostic[name], name) end fastiterate(state.namespaces) do ns explicit_step!(ns, grid, timestepper, Δt) @@ -139,3 +143,6 @@ end u[i, j, 1] += ∂u∂t[i, j] * Δt end end + +# Default debug hooks +@inline debughook!(::typeof(explicit_step!), field, name) = checkfinite!(field, name) From 5e1d9d3116033e585ab4fa4b2e47f0d88678c067 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Thu, 23 Apr 2026 00:15:11 +0200 Subject: [PATCH 13/27] Fix issues with infinite matric potentials in dry soil --- .../soil/hydrology/soil_hydraulic_closures.jl | 12 +++++++----- .../soil/hydrology/soil_hydraulic_properties.jl | 4 ++-- src/processes/soil/hydrology/soil_hydrology.jl | 14 ++++++++++---- .../evapotranspiration/ground_resistance_factor.jl | 8 ++++++-- 4 files changed, 25 insertions(+), 13 deletions(-) diff --git a/src/processes/soil/hydrology/soil_hydraulic_closures.jl b/src/processes/soil/hydrology/soil_hydraulic_closures.jl index 58220fa1c..88ee53c86 100644 --- a/src/processes/soil/hydrology/soil_hydraulic_closures.jl +++ b/src/processes/soil/hydrology/soil_hydraulic_closures.jl @@ -94,9 +94,7 @@ end ψm = ψ - ψh - ψz swrc = get_swrc(hydrology) por = porosity(i, j, k, grid, fields, strat, bgc) - sat_res = residual_saturation(get_hydraulic_properties(hydrology)) - θres = sat_res * por - vol_water_ice_content = swrc(ψm; θsat = por, θres) + vol_water_ice_content = swrc(ψm; θsat = por) saturation_water_ice[i, j, k] = vol_water_ice_content / por return nothing end @@ -117,9 +115,9 @@ end inv_swrc = inv(get_swrc(hydrology)) por = porosity(i, j, k, grid, fields, strat, bgc) sat_res = residual_saturation(get_hydraulic_properties(hydrology)) - θres = sat_res * por # compute matric pressure head - ψm = inv_swrc(sat * por; θsat = por, θres) + ψm = inv_swrc(max(sat * por, sat_res); θsat = por) + @assert isfinite(ψm) "NaN/infinite matric potential ψm = $ψm with sat=$sat por=$por resid=$sat_res" # compute elevation pressure head ψz = z - z_ref # compute hydrostatic pressure head assuming impermeable lower boundary @@ -151,3 +149,7 @@ end i, j, k = @index(Global, NTuple) saturation_to_pressure!(out.pressure_head, i, j, k, grid, fields, closure, args...) end + +# Default debug hooks +@inline debughook!(::typeof(pressure_to_saturation_kernel!), out, args...) = checkfinite!(out) +@inline debughook!(::typeof(saturation_to_pressure_kernel!), out, args...) = checkfinite!(out) diff --git a/src/processes/soil/hydrology/soil_hydraulic_properties.jl b/src/processes/soil/hydrology/soil_hydraulic_properties.jl index 09e54ded1..e3ddce5ed 100644 --- a/src/processes/soil/hydrology/soil_hydraulic_properties.jl +++ b/src/processes/soil/hydrology/soil_hydraulic_properties.jl @@ -87,7 +87,7 @@ $TYPEDFIELDS wilting_point::NF = 0.05 "Residual (minimum) saturation level [-]" - residual::NF = 1.0e-3 + residual::NF = 0.01 end function ConstantSoilHydraulics( @@ -141,7 +141,7 @@ $TYPEDFIELDS field_capacity_exp::NF = 0.35 "Residual (minimum) saturation level [-]" - residual::NF = 1.0e-3 + residual::NF = 0.01 end function SoilHydraulicsSURFEX( diff --git a/src/processes/soil/hydrology/soil_hydrology.jl b/src/processes/soil/hydrology/soil_hydrology.jl index 883500f82..45ceace18 100644 --- a/src/processes/soil/hydrology/soil_hydrology.jl +++ b/src/processes/soil/hydrology/soil_hydrology.jl @@ -183,14 +183,15 @@ grid cell and redistributes excess water upward layer-by-layer until reaching th any remaining excess water is added to the `surface_excess_water` pool. """ @propagate_inbounds function adjust_saturation_profile!(out, i, j, grid, hydrology::SoilHydrology{NF}) where {NF} - sat = out.saturation_water_ice - surface_excess_water = out.surface_excess_water props = get_hydraulic_properties(hydrology) sat_min = residual_saturation(props) field_grid = get_field_grid(grid) N = field_grid.Nz + sat = out.saturation_water_ice + surface_excess_water = out.surface_excess_water - # First iterate over soil layers from bottom to top + # First iterate over soil layers from bottom to top, transferring water from + # overfilled layers to the layer above for k in 1:(N - 1) # calculate excess saturation excess_sat = max(sat[i, j, k] - one(NF), zero(NF)) @@ -200,7 +201,7 @@ any remaining excess water is added to the `surface_excess_water` pool. sat[i, j, k + 1] += excess_sat * Δzᵃᵃᶜ(i, j, k, field_grid) / Δzᵃᵃᶜ(i, j, k + 1, field_grid) end - # then from top to bottom + # then from top to bottom, extracting water for underfilled cells from layers below for k in N:-1:2 # calculate saturation deficit from residual saturation level deficit_sat = max(-sat[i, j, k] + sat_min, zero(NF)) @@ -300,3 +301,8 @@ end i, j, k = @index(Global, NTuple) compute_hydraulics!(out, i, j, k, grid, fields, hydrology, args...) end + +# Default debug hooks +@inline debughook!(::typeof(compute_water_table_kernel!), out, args...) = checkfinite!(out) +@inline debughook!(::typeof(adjust_saturation_profile_kernel!), out, args...) = checkfinite!(out) +@inline debughook!(::typeof(compute_hydraulics_kernel!), out, args...) = checkfinite!(out) diff --git a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl index eaaedf4ed..e8172f3fd 100644 --- a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl +++ b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl @@ -38,7 +38,7 @@ ground_evaporation_resistance_factor(i, j, grid, fields, ::SoilMoistureResistanc @inline function ground_evaporation_resistance_factor( i, j, grid, fields, - ::SoilMoistureResistanceFactor{NF}, + res::SoilMoistureResistanceFactor{NF}, soil::AbstractSoil ) where {NF} fgrid = get_field_grid(grid) @@ -50,7 +50,11 @@ ground_evaporation_resistance_factor(i, j, grid, fields, ::SoilMoistureResistanc fracs = volumetric_fractions(soil) θfc = field_capacity(get_hydraulic_properties(hydrology), texture) θw = fracs.water - if fracs.water < θfc + return ground_evaporation_resistance_factor(res, θw, θfc) +end + +@inline function ground_evaporation_resistance_factor(::SoilMoistureResistanceFactor{NF}, θw, θfc) where {NF} + if θw < θfc β = (1 - cos(π * θw / θfc))^2 / 4 else β = NF(1) From 2b0f8749e57ad1c82f0808b9a38236f00736770f Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Thu, 23 Apr 2026 13:13:15 +0200 Subject: [PATCH 14/27] Account for residual water content in resistance factor --- .../evapotranspiration/ground_resistance_factor.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl index e8172f3fd..22f4538d5 100644 --- a/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl +++ b/src/processes/surface_hydrology/evapotranspiration/ground_resistance_factor.jl @@ -45,17 +45,20 @@ ground_evaporation_resistance_factor(i, j, grid, fields, ::SoilMoistureResistanc strat = get_stratigraphy(soil) hydrology = get_hydrology(soil) bgc = get_biogeochemistry(soil) + props = get_hydraulic_properties(hydrology) soil = soil_volume(i, j, fgrid.Nz, grid, fields, strat, hydrology, bgc) texture = mineral_texture(soil) fracs = volumetric_fractions(soil) + # Get field capacity, water content, and residual water content θfc = field_capacity(get_hydraulic_properties(hydrology), texture) θw = fracs.water - return ground_evaporation_resistance_factor(res, θw, θfc) + θres = props.residual * soil.saturation + return ground_evaporation_resistance_factor(res, θw, θfc, θres) end -@inline function ground_evaporation_resistance_factor(::SoilMoistureResistanceFactor{NF}, θw, θfc) where {NF} +@inline function ground_evaporation_resistance_factor(::SoilMoistureResistanceFactor{NF}, θw, θfc, θres) where {NF} if θw < θfc - β = (1 - cos(π * θw / θfc))^2 / 4 + β = (1 - cos(π * (θw - θres) / (θfc - θres)))^2 / 4 else β = NF(1) end From 264acb8da79b9a45b33c8133358d76ecd0cdae03 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Thu, 23 Apr 2026 18:32:58 +0200 Subject: [PATCH 15/27] Formatting changes in radiative_fluxes.jl --- .../surface_energy/radiative_fluxes.jl | 56 +++++++++---------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/processes/surface_energy/radiative_fluxes.jl b/src/processes/surface_energy/radiative_fluxes.jl index 348db8b8b..6f5f44df9 100644 --- a/src/processes/surface_energy/radiative_fluxes.jl +++ b/src/processes/surface_energy/radiative_fluxes.jl @@ -78,27 +78,46 @@ struct DiagnosedRadiativeFluxes{NF} <: AbstractRadiativeFluxes{NF} end DiagnosedRadiativeFluxes(::Type{NF}) where {NF} = DiagnosedRadiativeFluxes{NF}() """ - compute_shortwave_up(::DiagnosedRadiativeFluxes, surface_shortwave_down, α) + compute_shortwave_up(::DiagnosedRadiativeFluxes, S_down, α) -Compute outgoing shortwave radiation from the incoming `surface_shortwave_down` and albedo `α`. +Compute outgoing shortwave radiation from the incoming shortwave radiation `S_down` and albedo `α`. """ -@inline function compute_shortwave_up(::DiagnosedRadiativeFluxes, surface_shortwave_down, α) - surface_shortwave_up = α * surface_shortwave_down - return surface_shortwave_up +@inline function compute_shortwave_up(::DiagnosedRadiativeFluxes, S_down, α) + S_up = α * S_down + return S_up end """ - compute_longwave_up(::DiagnosedRadiativeFluxes, constants::PhysicalConstants, surface_longwave_down, Ts, ϵ) + compute_longwave_up(::DiagnosedRadiativeFluxes, constants::PhysicalConstants, L_down, Ts, ϵ) -Compute outgoing longwave radiation from incoming `surface_longwave_down`, surface temperature `Ts`, and emissivity `ϵ`. +Compute outgoing longwave radiation from incoming longwave radiation `L_down`, surface temperature `Ts`, and emissivity `ϵ`. """ -@inline function compute_longwave_up(::DiagnosedRadiativeFluxes, constants::PhysicalConstants, surface_longwave_down, Ts, ϵ) +@inline function compute_longwave_up(::DiagnosedRadiativeFluxes, constants::PhysicalConstants, L_down, Ts, ϵ) T = celsius_to_kelvin(constants, Ts) + L_emit = stefan_boltzmann(constants, T, ϵ) # outgoing LW radiation is the sum of the radiant emittance and the residual incoming radiation - surface_longwave_up = stefan_boltzmann(constants, T, ϵ) + (1 - ϵ) * surface_longwave_down + surface_longwave_up = L_emit + (1 - ϵ) * L_down return surface_longwave_up end +""" + $TYPEDEF + +Compute the net radiation budget given incoming and outgoing shortwave and longwave radiation. +""" +@inline function compute_surface_net_radiation( + ::AbstractRadiativeFluxes, + SW_up, + SW_down, + LW_up, + LW_down + ) + # Sum up radiation fluxes; note that by convention fluxes are positive upward + rad_net = SW_up - SW_down + LW_up - LW_down + return rad_net +end + + ## Top-level interface methods variables(::DiagnosedRadiativeFluxes) = ( @@ -189,25 +208,6 @@ end return surface_net_radiation end -# Common implementations - -""" - $TYPEDEF - -Compute the net radiation budget given incoming and outgoing shortwave and longwave radiation. -""" -@inline function compute_surface_net_radiation( - ::AbstractRadiativeFluxes, - SW_up, - SW_down, - LW_up, - LW_down - ) - # Sum up radiation fluxes; note that by convention fluxes are positive upward - rad_net = SW_up - SW_down + LW_up - LW_down - return rad_net -end - ## Kernels @kernel inbounds = true function compute_auxiliary_kernel!(out, grid, fields, rad::AbstractRadiativeFluxes, args...) From 71e5cec3cd3eb07ea6cfd4f3a08f9566ec158b92 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 24 Apr 2026 15:13:59 +0200 Subject: [PATCH 16/27] Fix typo in compute_latent_heat_flux args --- src/processes/surface_energy/turbulent_fluxes.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl index e49a0cba2..158beaafb 100644 --- a/src/processes/surface_energy/turbulent_fluxes.jl +++ b/src/processes/surface_energy/turbulent_fluxes.jl @@ -42,10 +42,10 @@ end $TYPEDSIGNATURES Compute the latent heat flux as a function of the humidity flux `Q_h` [m/s], the density `ρₐ` [kg/m³] of air, -and the specific latent heat of fusion `Lsl` [J/kg]. +and the specific latent heat of vaporization or sublimation `L` [J/kg]. """ -function compute_latent_heat_flux(::DiagnosedTurbulentFluxes, Q_h, ρₐ, Lsl) - Hₗ = Lsl * ρₐ * Q_h +function compute_latent_heat_flux(::DiagnosedTurbulentFluxes, Q_h, ρₐ, L) + Hₗ = L * ρₐ * Q_h return Hₗ end @@ -140,9 +140,9 @@ defined by `evtr` which is assumed to be already computed. evtr::AbstractEvapotranspiration, constants::PhysicalConstants ) - let L = constants.Llg, # specific latent heat of vaporization of water - ρₐ = constants.ρₐ, # density of air - Q_h = surface_humidity_flux(i, j, grid, fields, evtr) # humidity flux + let L = constants.Llg # specific latent heat of vaporization of water + ρₐ = constants.ρₐ # density of air + Q_h = surface_humidity_flux(i, j, grid, fields, evtr) # humidity flux # Calculate latent heat flux (positive upwards) Hₗ = compute_latent_heat_flux(tur, Q_h, ρₐ, L) return Hₗ From 8e8021de156618d0f60a7db1a0e9e51e126ce1bc Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 24 Apr 2026 15:14:42 +0200 Subject: [PATCH 17/27] Ensure that LAI/SAI are nonnegative in r_e calculation --- .../evapotranspiration/canopy_evapotranspiration.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl b/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl index 54a403a2b..3a431100a 100644 --- a/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl +++ b/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl @@ -138,7 +138,7 @@ for the given scheme `evapotranspiration` and process dependencies. ) # Get inputs Ts = fields.skin_temperature[i, j] # skin temperature (top of canopy) - Tg = fields.ground_temperature[i, j] # ground temeprature (top snow/soil layer) + Tg = fields.ground_temperature[i, j] # ground temperature (top snow/soil layer) gw_can = fields.canopy_water_conductance[i, j] # stomatal conductance (assumed to be defined by vegetation) # Compute VPD and resistance terms @@ -164,11 +164,11 @@ Compute the aerodynamic resistance between the ground and canopy as a function o @inline function aerodynamic_resistance( i, j, grid, fields, atmos::AbstractAtmosphere, - evapotranspiration::PALADYNCanopyEvapotranspiration, + evapotranspiration::PALADYNCanopyEvapotranspiration{NF}, ::AbstractVegetation # included just to make explicit the dependence on vegetation fields - ) - @inbounds let LAI = fields.leaf_area_index[i, j], - SAI = fields.SAI[i, j], + ) where {NF} + @inbounds let LAI = max(fields.leaf_area_index[i, j], zero(NF)), + SAI = max(fields.SAI[i, j], zero(NF)), Vₐ = windspeed(i, j, grid, fields, atmos), C = evapotranspiration.C_can # drag coefficient for the canopy rₙ = (1 - exp(-LAI - SAI)) / (C * Vₐ) From 68ba75ca61296f49658b0d7b196b1a6dfeb61d01 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sun, 26 Apr 2026 18:50:17 +0200 Subject: [PATCH 18/27] Fix formatting in prescribed_atmosphere.jl --- src/processes/atmosphere/prescribed_atmosphere.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl index 240864b5a..7093579bc 100644 --- a/src/processes/atmosphere/prescribed_atmosphere.jl +++ b/src/processes/atmosphere/prescribed_atmosphere.jl @@ -161,8 +161,8 @@ Retrieve or compute the specific_humidity at the current time step. Computes the specific humidity (vapor pressure) deficit over a surface at temperature `Ts` from the current atmospheric fields. """ @propagate_inbounds function compute_humidity_vpd(i, j, grid, fields, atmos::AbstractAtmosphere, c::PhysicalConstants, Ts = nothing) - let Δe = compute_vpd(i, j, grid, fields, atmos, c, Ts), - p = air_pressure(i, j, grid, fields, atmos) + let Δe = compute_vpd(i, j, grid, fields, atmos, c, Ts) + p = air_pressure(i, j, grid, fields, atmos) Δq = vapor_pressure_to_specific_humidity(Δe, p, c.ε) return Δq end From 6c51fc6ebb4503db3c6b87dac2ffbb2aed50b329 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sun, 26 Apr 2026 18:50:30 +0200 Subject: [PATCH 19/27] Add temporary compat bounds in examples project --- examples/Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/Project.toml b/examples/Project.toml index 0506e53c2..d5cb56b51 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -19,5 +19,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Terrarium = "80418d68-07fa-499d-ae2b-c12e531f5cd8" +[compat] +SpeedyWeather = "0.18.1" + [sources.Terrarium] path = ".." From 2eb9629800873fc287ca8d67600de437ccee0402 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sun, 26 Apr 2026 18:51:46 +0200 Subject: [PATCH 20/27] More tweaks to speedy primitive wet script --- examples/simulations/speedy_wet_land.jl | 30 +++++++++++++------------ 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl index 3019ea99d..3a4800787 100644 --- a/examples/simulations/speedy_wet_land.jl +++ b/examples/simulations/speedy_wet_land.jl @@ -93,7 +93,7 @@ function speedy_timestep!( set!(state.inputs.surface_longwave_down, LwIn) # run land forward over speedy timestep interval; # we use a smaller actual timestep to ensure stability - Terrarium.run!(land.integrator, period = progn.clock.Δt, Δt = 60.0) + Terrarium.run!(land.integrator, period = progn.clock.Δt, Δt = 300.0) # Update speedy variables progn.land.soil_temperature .= state.skin_temperature .+ NF(273.15) progn.land.sensible_heat_flux .= state.sensible_heat_flux @@ -111,11 +111,11 @@ sim = Speedy.initialize!(primitive_wet) Speedy.run!(sim, period = Day(1)) Nz = 30 -Δz_min = 0.05 # currently the coupling is only stable with a large surface layer +Δz_min = 0.05 grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(; N = Nz, Δz_min), ring_grid) # Initial conditions soil_initializer = SoilInitializer(eltype(grid)) -soil = SoilEnergyWaterCarbon(eltype(grid), hydrology = SoilHydrology(eltype(grid), RichardsEq())) +soil = SoilEnergyWaterCarbon(eltype(grid), hydrology = SoilHydrology(eltype(grid))) # Land model with "prescribed" atmosphere (from the perspective of the land model at least...) # vegetation = PrescribedVegetationCarbon(eltype(grid)) model = LandModel(grid; initializer = soil_initializer, vegetation = nothing, soil) @@ -123,7 +123,7 @@ initializers = (;) integrator = initialize(model, ForwardEuler(eltype(grid)); initializers) # check if land model works standalone (with default atmospheric state) timestep!(integrator, 60.0) # one step -run!(integrator, period = Hour(1), Δt = 120.0) # one hour +run!(integrator, period = Hour(1), Δt = 300.0) # one hour Terrarium.initialize!(integrator) # reinitialize before setting up atmosphere # Initialize Terrarium-Speedy land model @@ -148,28 +148,30 @@ Speedy.add!(primitive_wet_coupled.output, Speedy.SoilTemperatureOutput()) # initialize coupled simulation sim_coupled = Speedy.initialize!(primitive_wet_coupled) # run it -Speedy.run!(sim_coupled, period = Hour(4)) +period = Day(1) +@info "Running simulation for $period" +@time Speedy.run!(sim_coupled, period = period) +Terrarium.checkfinite!(integrator.state) -# Soil temperature in the 3rd layer +# Land variables Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 2], grid), title = "", size = (800, 400)) Tsurf_fig = heatmap(RingGrids.Field(interior(integrator.state.skin_temperature)[:, 1], grid), title = "", size = (800, 400)) -sat_fig = heatmap(RingGrids.Field(interior(integrator.state.saturation_water_ice)[:, 1, end], grid), title = "", size = (800, 400)) +Hs_fig = heatmap(RingGrids.Field(interior(integrator.state.sensible_heat_flux)[:, 1], grid), title = "", size = (800, 400)) +Hl_fig = heatmap(RingGrids.Field(interior(integrator.state.latent_heat_flux)[:, 1], grid), title = "", size = (800, 400)) E_fig = heatmap(RingGrids.Field(interior(integrator.state.evaporation_ground)[:, 1], grid), title = "", size = (800, 400)) +sat_fig = heatmap(RingGrids.Field(interior(integrator.state.saturation_water_ice)[:, 1, end], grid), title = "", size = (800, 400)) # Atmosphere variables Tair_fig = heatmap(sim_coupled.diagnostic_variables.grid.temp_grid[:, 8] .- 273.15, title = "Air temperature", size = (800, 400)) pres_fig = heatmap(exp.(sim_coupled.diagnostic_variables.grid.pres_grid), title = "Surface pressure", size = (800, 400)) -srad_fig = heatmap(exp.(sim_coupled.diagnostic_variables.physics.surface_shortwave_down), title = "Surface shortwave down", size = (800, 400)) - -# animate surface air and soil temperatures -Speedy.animate(sim_coupled, variable = "temp", coastlines = false, level = spectral_grid.nlayers, output_file = "plots/speedy_terrarium_dry_air_temperature.mp4") -# Speedy.animate(sim_coupled, variable = "st", coastlines = false, level = 1, output_file = "plots/speedy_terrarium_dry_soil_temperature.mp4") # TODO: this is broken now for some reason? +srad_fig = heatmap(sim.diagnostic_variables.physics.surface_shortwave_down, title = "Surface shortwave down", size = (800, 400)) # pick a point somewhere in the mid-lattitudes T = interior(integrator.state.temperature)[2000, 1, :] sat = interior(integrator.state.saturation_water_ice)[2000, 1, :] f = interior(integrator.state.liquid_water_fraction)[2000, 1, :] zs = znodes(integrator.state.temperature) + # Plot temperature and liquid fraction profiles in upper 15 layers -Makie.scatterlines(T[(end - 15):end], zs[(end - 15):end]) -Makie.scatterlines(sat[(end - 15):end], zs[(end - 15):end]) +Makie.scatterlines(T, zs) +Makie.scatterlines(sat, zs) Makie.scatterlines(f, zs) From bea0aebc19b50e1fbf8041c6414f0fe3671e666a Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 27 Apr 2026 15:20:06 +0200 Subject: [PATCH 21/27] Fix failing soil hydrology test --- src/processes/soil/hydrology/soil_hydrology.jl | 7 +++---- test/soil/soil_hydrology_tests.jl | 8 ++++---- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/processes/soil/hydrology/soil_hydrology.jl b/src/processes/soil/hydrology/soil_hydrology.jl index 45ceace18..1b13c87f1 100644 --- a/src/processes/soil/hydrology/soil_hydrology.jl +++ b/src/processes/soil/hydrology/soil_hydrology.jl @@ -252,10 +252,9 @@ the porosity and is thus not the same as the saturation tendency. evtr::Optional{AbstractEvapotranspiration} ) where {NF} # Compute divergence of water fluxes due to forcings only - ∂θ∂t = ( - + forcing(i, j, k, grid, clock, fields, evtr, hydrology, constants) # ET forcing - + forcing(i, j, k, grid, clock, fields, hydrology.vwc_forcing, hydrology) # generic user-defined forcing - ) + ET_loss = forcing(i, j, k, grid, clock, fields, evtr, hydrology, constants) # ET forcing + F = forcing(i, j, k, grid, clock, fields, hydrology.vwc_forcing, hydrology) # generic user-defined forcing + ∂θ∂t = ET_loss + F return ∂θ∂t end diff --git a/test/soil/soil_hydrology_tests.jl b/test/soil/soil_hydrology_tests.jl index 0455ad0c0..a11b49f69 100644 --- a/test/soil/soil_hydrology_tests.jl +++ b/test/soil/soil_hydrology_tests.jl @@ -6,7 +6,7 @@ using FreezeCurves using Oceananigans @testset "Hydraulic properties (constant)" begin - # For prescribed hyraulic properties, just check that the returned values + # For prescribed hydraulic properties, just check that the returned values # match what was set. hydraulic_props = ConstantSoilHydraulics( Float64; @@ -113,13 +113,13 @@ end compute!(∫sat_deficit) Terrarium.adjust_saturation_profile!(state, grid, hydrology) ∫sat₁ = compute!(∫sat)[1, 1, 1] - @test all(state.saturation_water_ice .>= 0) - @test all(∫sat₁ - ∫sat₀ .≈ 0) + @test all(state.saturation_water_ice .>= hydraulic_properties.residual) + @test ∫sat₁ ≈ ∫sat₀ # Case 3: Completely dry with negative saturation near surface set!(state.saturation_water_ice, (x, z) -> min(-0.1 - z, 0.0)) Terrarium.adjust_saturation_profile!(state, grid, hydrology) - @test all(state.saturation_water_ice .≈ 0) + @test all(state.saturation_water_ice .≈ hydraulic_properties.residual) end @testset "SoilHydrology: Richardson-Richards' equation" begin From 9a403d7dc1c4c099ff5e9e07811fbb133b7e54b3 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 27 Apr 2026 15:21:07 +0200 Subject: [PATCH 22/27] Change TracerGas type and default incoming radiation values --- examples/simulations/land_column.jl | 7 ++++- .../atmosphere/prescribed_atmosphere.jl | 27 ++++++++++--------- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/examples/simulations/land_column.jl b/examples/simulations/land_column.jl index 3a8cc2b40..e8df98cb9 100644 --- a/examples/simulations/land_column.jl +++ b/examples/simulations/land_column.jl @@ -14,9 +14,14 @@ vegetation = VegetationCarbon(eltype(grid)) land = LandModel(grid; soil, vegetation) # Variably saturated with water table at roughly 5 m depth initializers = ( + temperature = 15.0, saturation_water_ice = (x, z) -> min(1, 0.5 - 0.1 * z), - carbon_vegetation = 0.1, + carbon_vegetation = 0.5, ) integrator = @time initialize(land, ForwardEuler(); initializers); +# manually set atmospheric inputs to different values +set!(integrator.state.windspeed, 1.0) # 1 m/s +set!(integrator.state.specific_humidity, 1.0e-4) # kg/kg Δt = 60.0 @time timestep!(integrator, Δt) +@show integrator.state.latent_heat_flux[1, 1, 1] diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl index 7093579bc..fc45e1725 100644 --- a/src/processes/atmosphere/prescribed_atmosphere.jl +++ b/src/processes/atmosphere/prescribed_atmosphere.jl @@ -1,23 +1,24 @@ """ Generic type representing the concentration of a particular tracer gas in the atmosphere. """ -@kwdef struct TracerGas{name} - TracerGas(name::Symbol) = new{name}() +@kwdef struct TracerGas{NF, name} + "Default concentration of the gas" + concentration::NF + + TracerGas(name::Symbol, concentration::NF) where {NF} = new{NF, name}(concentration) end -Base.nameof(::TracerGas{name}) where {name} = name +Base.nameof(::TracerGas{NF, name}) where {NF, name} = name -variables(::TracerGas{name}) where {name} = ( - input(name, XY(), default = default_tracer_conc(Val(name)), units = u"ppm", desc = "Ambient atmospheric $(name) concentration in ppm"), +variables(gas::TracerGas{NF, name}) where {NF, name} = ( + input(name, XY(), default = gas.concentration, units = u"ppm", desc = "Ambient atmospheric $(name) concentration in ppm"), ) -default_tracer_conc(::Val{:CO2}) = 380 # in ppm - """ Creates a `TracerGas` for ambient CO2 with concentration prescribed by an input variable with the given name. """ -AmbientCO2(name = :CO2) = TracerGas(name) +AmbientCO2(::Type{NF}, name = :CO2) where {NF} = TracerGas(name, NF(380)) """ Creates a `NamedTuple` from the given tracer gas types. @@ -49,7 +50,7 @@ struct PrescribedAtmosphere{ IncomingRad <: AbstractIncomingRadiation, Humidity <: AbstractHumidity, Aerodynamics <: AbstractAerodynamics, - Tracers <: NamedTuple{tracernames, <:Tuple{Vararg{TracerGas}}}, + Gases <: Tuple{Vararg{TracerGas{NF}}}, } <: AbstractAtmosphere{NF, Precip, IncomingRad, Humidity, Aerodynamics} "Surface-relative altitude in meters at which the atmospheric forcings are assumed to be applied" altitude::NF @@ -70,7 +71,7 @@ struct PrescribedAtmosphere{ aerodynamics::Aerodynamics "Atmospheric tracer gases" - tracers::Tracers + tracers::NamedTuple{tracernames, Gases} end function PrescribedAtmosphere( @@ -81,7 +82,7 @@ function PrescribedAtmosphere( radiation::AbstractIncomingRadiation = LongShortWaveRadiation(), humidity::AbstractHumidity = SpecificHumidity(), aerodynamics::AbstractAerodynamics = ConstantAerodynamics(NF), - tracers::NamedTuple = TracerGases(AmbientCO2()), + tracers::NamedTuple = TracerGases(AmbientCO2(NF)), ) where {NF} return PrescribedAtmosphere(altitude, min_windspeed, precip, radiation, humidity, aerodynamics, tracers) end @@ -218,8 +219,8 @@ length [hr]. struct LongShortWaveRadiation <: AbstractIncomingRadiation end variables(::LongShortWaveRadiation) = ( - input(:surface_shortwave_down, XY(), default = 300, units = u"W/m^2", desc = "Incoming (downwelling) shortwave solar radiation"), - input(:surface_longwave_down, XY(), default = 50, units = u"W/m^2", desc = "Incoming (downwelling) longwave thermal radiation"), + input(:surface_shortwave_down, XY(), default = 341, units = u"W/m^2", desc = "Incoming (downwelling) shortwave solar radiation"), + input(:surface_longwave_down, XY(), default = 333, units = u"W/m^2", desc = "Incoming (downwelling) longwave thermal radiation"), input(:daytime_length, XY(), default = 12, units = u"hr", desc = "Number of daytime hours varying with the season and orbital parameters"), ) From 61cb33efcd0d7fef30a1d96a3c1277636b3aef4e Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 27 Apr 2026 15:30:21 +0200 Subject: [PATCH 23/27] Rename SAI to stem_area_index --- .../canopy_interception/canopy_interception.jl | 4 ++-- .../evapotranspiration/canopy_evapotranspiration.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/processes/surface_hydrology/canopy_interception/canopy_interception.jl b/src/processes/surface_hydrology/canopy_interception/canopy_interception.jl index e5d274094..c6e9643ed 100644 --- a/src/processes/surface_hydrology/canopy_interception/canopy_interception.jl +++ b/src/processes/surface_hydrology/canopy_interception/canopy_interception.jl @@ -61,7 +61,7 @@ variables(::PALADYNCanopyInterception) = ( auxiliary(:saturation_canopy_water, XY(); desc = "Fraction of the canopy saturated with water"), auxiliary(:rainfall_ground, XY(); desc = "Rainfall rate reaching the ground", units = u"m/s"), input(:leaf_area_index, XY(); desc = "Leaf Area Index", units = u"m^2/m^2"), - input(:SAI, XY(); desc = "Stem Area Index", units = u"m^2/m^2"), + input(:stem_area_index, XY(); desc = "Stem Area Index", units = u"m^2/m^2"), ) @propagate_inbounds canopy_water(i, j, grid, fields, ::PALADYNCanopyInterception) = fields.canopy_water[i, j] @@ -181,7 +181,7 @@ end # Get inputs rain = rainfall(i, j, grid, fields, atmos) LAI = fields.leaf_area_index[i, j] - SAI = fields.SAI[i, j] + SAI = fields.stem_area_index[i, j] w_can = fields.canopy_water[i, j] # Compute canopy saturation faction diff --git a/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl b/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl index 3a431100a..271ad6d9e 100644 --- a/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl +++ b/src/processes/surface_hydrology/evapotranspiration/canopy_evapotranspiration.jl @@ -168,7 +168,7 @@ Compute the aerodynamic resistance between the ground and canopy as a function o ::AbstractVegetation # included just to make explicit the dependence on vegetation fields ) where {NF} @inbounds let LAI = max(fields.leaf_area_index[i, j], zero(NF)), - SAI = max(fields.SAI[i, j], zero(NF)), + SAI = max(fields.stem_area_index[i, j], zero(NF)), Vₐ = windspeed(i, j, grid, fields, atmos), C = evapotranspiration.C_can # drag coefficient for the canopy rₙ = (1 - exp(-LAI - SAI)) / (C * Vₐ) From 2c62db90b996c18a23b093ddacc75c54796ba64b Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 27 Apr 2026 17:22:22 +0200 Subject: [PATCH 24/27] Minor refactoring of soil hydrology code --- .../soil/hydrology/soil_hydrology.jl | 22 ++-------------- .../soil/hydrology/soil_hydrology_rre.jl | 25 +++++++++++++++++++ 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/src/processes/soil/hydrology/soil_hydrology.jl b/src/processes/soil/hydrology/soil_hydrology.jl index 1b13c87f1..74f1f2172 100644 --- a/src/processes/soil/hydrology/soil_hydrology.jl +++ b/src/processes/soil/hydrology/soil_hydrology.jl @@ -140,27 +140,9 @@ end """ $TYPEDSIGNATURES -Kernel function that computes soil hydraulics and unsaturated hydraulic conductivity. +Kernel function that computes dynamic soil hydraulic properties. """ -@propagate_inbounds function compute_hydraulics!( - out, i, j, k, grid, fields, - hydrology::SoilHydrology, - strat::AbstractStratigraphy, - bgc::AbstractSoilBiogeochemistry - ) - # Get underlying grid - fgrid = get_field_grid(grid) - # compute hydraulic conductivity - @inbounds if k <= 1 - out.hydraulic_conductivity[i, j, k] = hydraulic_conductivity(i, j, 1, fgrid, fields, hydrology, strat, bgc) - elseif k >= fgrid.Nz - out.hydraulic_conductivity[i, j, k] = hydraulic_conductivity(i, j, fgrid.Nz, fgrid, fields, hydrology, strat, bgc) - out.hydraulic_conductivity[i, j, k + 1] = out.hydraulic_conductivity[i, j, k] - else - out.hydraulic_conductivity[i, j, k] = min_zᵃᵃᶠ(i, j, k, fgrid, hydraulic_conductivity, fields, hydrology, strat, bgc) - end - return nothing -end +compute_hydraulics!(out, i, j, k, grid, hydrology::SoilHydrology, args...) = nothing """ $TYPEDSIGNATURES diff --git a/src/processes/soil/hydrology/soil_hydrology_rre.jl b/src/processes/soil/hydrology/soil_hydrology_rre.jl index 021f31250..4f530d238 100644 --- a/src/processes/soil/hydrology/soil_hydrology_rre.jl +++ b/src/processes/soil/hydrology/soil_hydrology_rre.jl @@ -145,6 +145,31 @@ Compute the hydraulic conductivity at the center of the grid cell `i, j, k`. return hydraulic_conductivity(hydrology.hydraulic_properties, soil) end +""" + $TYPEDSIGNATURES + +Computes the unsaturated hydraulic conductivity for `RichardsEq` configurations of `SoilHydrology`. +""" +@propagate_inbounds function compute_hydraulics!( + out, i, j, k, grid, fields, + hydrology::SoilHydrology{NF, RichardsEq}, + strat::AbstractStratigraphy, + bgc::AbstractSoilBiogeochemistry + ) + # Get underlying grid + fgrid = get_field_grid(grid) + # compute hydraulic conductivity + @inbounds if k <= 1 + out.hydraulic_conductivity[i, j, k] = hydraulic_conductivity(i, j, 1, fgrid, fields, hydrology, strat, bgc) + elseif k >= fgrid.Nz + out.hydraulic_conductivity[i, j, k] = hydraulic_conductivity(i, j, fgrid.Nz, fgrid, fields, hydrology, strat, bgc) + out.hydraulic_conductivity[i, j, k + 1] = out.hydraulic_conductivity[i, j, k] + else + out.hydraulic_conductivity[i, j, k] = min_zᵃᵃᶠ(i, j, k, fgrid, hydraulic_conductivity, fields, hydrology, strat, bgc) + end + return nothing +end + # Kernels @kernel inbounds = true function compute_tendencies_kernel!( From 05cb25c6c4207a06f23eb3bdb75a526daabf374e Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 29 Apr 2026 14:11:04 +0200 Subject: [PATCH 25/27] Add soil moisture in wet land --- examples/simulations/speedy_wet_land.jl | 7 +++++++ src/processes/soil/hydrology/soil_hydrology.jl | 2 +- src/processes/soil/hydrology/soil_hydrology_rre.jl | 2 +- 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl index 3a4800787..1229f4600 100644 --- a/examples/simulations/speedy_wet_land.jl +++ b/examples/simulations/speedy_wet_land.jl @@ -10,6 +10,9 @@ using CairoMakie, GeoMakie import RingGrids import SpeedyWeather as Speedy +# Choose architecture based on available hardware +arch = CUDA.functional() ? GPU() : CPU() + """ Naive implementation of a SpeedyWeather "wet" land model based on Terrarium. """ @@ -38,6 +41,7 @@ end Speedy.variables(land::TerrariumWetLand) = ( Speedy.PrognosticVariable(name = :soil_temperature, dims = Speedy.Grid3D(), namespace = :land), + Speedy.PrognosticVariable(name = :soil_moisture, dims = Speedy.Grid3D(), namespace = :land), ) function Speedy.initialize!( @@ -49,7 +53,9 @@ function Speedy.initialize!( ) where {NF} Terrarium.initialize!(land.integrator) Tsoil = interior(land.integrator.state.temperature)[:, 1, end] .+ 273.15 + sat = interior(land.integrator.state.saturation_water_ice)[:, 1, end] progn.land.soil_temperature .= Tsoil + progn.land.soil_moisture .= sat return nothing end @@ -96,6 +102,7 @@ function speedy_timestep!( Terrarium.run!(land.integrator, period = progn.clock.Δt, Δt = 300.0) # Update speedy variables progn.land.soil_temperature .= state.skin_temperature .+ NF(273.15) + progn.land.soil_moisture .= interior(state.saturation_water_ice)[:, 1, end] progn.land.sensible_heat_flux .= state.sensible_heat_flux progn.land.surface_humidity_flux .= state.latent_heat_flux ./ consts.Llg diagn.physics.surface_longwave_up .= state.surface_longwave_up diff --git a/src/processes/soil/hydrology/soil_hydrology.jl b/src/processes/soil/hydrology/soil_hydrology.jl index 74f1f2172..82046c27a 100644 --- a/src/processes/soil/hydrology/soil_hydrology.jl +++ b/src/processes/soil/hydrology/soil_hydrology.jl @@ -142,7 +142,7 @@ end Kernel function that computes dynamic soil hydraulic properties. """ -compute_hydraulics!(out, i, j, k, grid, hydrology::SoilHydrology, args...) = nothing +compute_hydraulics!(out, i, j, k, grid, fields, hydrology::SoilHydrology, args...) = nothing """ $TYPEDSIGNATURES diff --git a/src/processes/soil/hydrology/soil_hydrology_rre.jl b/src/processes/soil/hydrology/soil_hydrology_rre.jl index 4f530d238..0cd3fe48f 100644 --- a/src/processes/soil/hydrology/soil_hydrology_rre.jl +++ b/src/processes/soil/hydrology/soil_hydrology_rre.jl @@ -155,7 +155,7 @@ Computes the unsaturated hydraulic conductivity for `RichardsEq` configurations hydrology::SoilHydrology{NF, RichardsEq}, strat::AbstractStratigraphy, bgc::AbstractSoilBiogeochemistry - ) + ) where {NF} # Get underlying grid fgrid = get_field_grid(grid) # compute hydraulic conductivity From 57610cb4144a897c7f5246f6623e629f8a09a5c0 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Wed, 29 Apr 2026 16:58:34 +0200 Subject: [PATCH 26/27] Fix checkfinite! call --- examples/simulations/speedy_wet_land.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/simulations/speedy_wet_land.jl b/examples/simulations/speedy_wet_land.jl index 1229f4600..5da4491a8 100644 --- a/examples/simulations/speedy_wet_land.jl +++ b/examples/simulations/speedy_wet_land.jl @@ -158,7 +158,7 @@ sim_coupled = Speedy.initialize!(primitive_wet_coupled) period = Day(1) @info "Running simulation for $period" @time Speedy.run!(sim_coupled, period = period) -Terrarium.checkfinite!(integrator.state) +Terrarium.checkfinite!(integrator.state.prognostic) # Land variables Tsoil_fig = heatmap(RingGrids.Field(interior(integrator.state.temperature)[:, 1, end - 2], grid), title = "", size = (800, 400)) From 9ab8cb427572abd1ac8fa3cb19e279e12749a60d Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sun, 3 May 2026 23:24:15 +0200 Subject: [PATCH 27/27] Increase number of iterations in SEB skinT tests --- test/surface_energy/skin_temperature.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/surface_energy/skin_temperature.jl b/test/surface_energy/skin_temperature.jl index 4822273b5..c67d5e14c 100644 --- a/test/surface_energy/skin_temperature.jl +++ b/test/surface_energy/skin_temperature.jl @@ -36,7 +36,7 @@ end Tskin_old = deepcopy(state.skin_temperature) resid = nothing balance = nothing - for i in 1:10 + for i in 1:20 # compute fluxes Terrarium.compute_surface_energy_fluxes!(state, grid, seb, model.constants, model.atmosphere) # diagnose skin temperature @@ -62,7 +62,7 @@ end @test all(state.sensible_heat_flux .> 0) # check that skin temperature converges for a large number of iterations Tskin_old = deepcopy(state.skin_temperature) - for i in 1:10 + for i in 1:20 # compute fluxes Terrarium.compute_surface_energy_fluxes!(state, grid, seb, model.constants, model.atmosphere) # diagnose skin temperature