Refactor ground heat flux into separate subprocess of SEB and add `PrescribedGroundHeatFlux` by bgroenks96 · Pull Request #200 · NumericalEarth/Terrarium.jl

bgroenks96 · 2026-09-25T12:25:14Z

Extracts the ground heat flux out of AbstractSkinTemperature and into its own
surface-energy-balance sub-process, so it gets the same Prescribed/Diagnosed
pair that the radiative and turbulent fluxes already have.

Implements docs/dev/2026-09/2026-09-01_PLAN_prescribed_surface_energy_balance.md
(Rev 2, with three deviations logged as Rev 3).

Why

G had no process of its own: the accessor, the residual closure, and roughly
eight compute_ground_heat_flux* methods all dispatched on AbstractSkinTemperature,
and the field was declared twice. It could not be prescribed without also replacing
the skin temperature scheme, which is exactly what a coupler needs to do. In the
coupled EarthSystemModel, a value written into ground_heat_flux by
InterfaceComputations did not survive, because Terrarium recomputed it as the
residual before the soil tendency read it.

What changed

New AbstractGroundHeatFlux with DiagnosedGroundHeatFlux and
PrescribedGroundHeatFlux, in src/processes/surface/ground_heat_flux.jl.
The compute_ground_heat_flux* family moves there from skin_temperature.jl.
ground_heat_flux is now declared exactly once, by the sub-process: auxiliary
when diagnosed, input when prescribed. Without snow it is the soil-top
boundary condition, so a prescribed value routes straight to the soil column.
SurfaceEnergyBalance gains a fifth sub-process field, defaulting to
DiagnosedGroundHeatFlux. Its type parameters are reordered to match the field
order, fixing the pre-existing radiative/turbulent swap.
PrescribedSurfaceEnergyBalance is a type alias for the all-prescribed
configuration, plus a convenience constructor.
New doc page for the ground heat flux, including the supported-combination table.
skin_temperature.md is retitled "Skin temperature" and its inbound @refs updated.

Note for reviewers

The atmosphere-side demand dispatches on AbstractGroundHeatFlux, not on
DiagnosedGroundHeatFlux alone, and is threaded into compute_skin_temperature
as a new ghf argument. This is what makes ImplicitSkinTemperature +
PrescribedGroundHeatFlux well posed. Without it, that pairing would invert
against R_net + H_s + H_l while storing an unrelated G. It is the main
signature change outside the new file, and it also touches
diagnose_skin_temperature_residual.

DiagnosedGroundHeatFlux still dispatches on the skin temperature scheme: the
plan described it as always storing R_net + H_s + H_l, but that is only the
PrescribedSkinTemperature case, while ImplicitSkinTemperature stores the
explicit conductive flux. Both methods moved across unchanged.

Verification

Default behavior is unchanged bit for bit. A 30-step deterministic LandModel
run dumping ground_heat_flux, skin_temperature, surface_net_radiation, both
turbulent fluxes, and the final internal_energy profile gives a maximum absolute
difference of exactly 0.0 against the base commit, over 180 values.
Full Pkg.test() suite passes; new tests in test/surface/ground_heat_flux.jl
cover variable classification, the alias and its type parameter ordering, the
no-overwrite guarantee, and energy conservation (ΔU = −G·Δt·N).
New Enzyme adjoint of a timestep! with respect to the prescribed
ground_heat_flux in test/differentiability/prescribed_seb_diff.jl.
All four supported combinations are type-stable and allocate identically
(3504 B in solve_surface_energy_balance!, zero in compute_auxiliary!).
Draft doc build is clean of docs_block errors.

Out of scope

Trimming NumericalEarthTerrariumExt's update_net_fluxes! and the coupled
end-to-end ERA5 run both live in NumericalEarth.jl and follow separately.
Moisture is not prescribed here; the vapor seam is the immediate follow-up.

🤖 Generated with Claude Code

bgroenks96 · 2026-09-25T15:18:27Z

@olivierbonte I haven't gone through this myself yet so I will also do a review here.

bgroenks96 · 2026-09-25T21:52:55Z

+land model configuration: an implicitly solved skin temperature with all surface fluxes diagnosed
+internally.
+
+Not every combination of `skin_temperature` and `ground_heat_flux` is well posed; the supported


I don't think "well posed" is the right term here.

replace by "valid"? I also don't like the term well posed in this context, I think it should be reserverd for the mathematical context where this terms orginiates from

bgroenks96 · 2026-09-25T21:56:06Z

+"""
+@propagate_inbounds function compute_ground_heat_flux_demand(i, j, grid, fields, ghf::DiagnosedGroundHeatFlux)
+    # Get individual flux terms
+    R_net = fields.surface_net_radiation[i, j, end]


This works but violates our convention of only accessing fields defined by the function arguments; i.e. those which would be returned by get_fields(state, args...). We could just add seb to the signature.

olivierbonte

What is conceptually not that clear to me: when we are prescribing G, the flux into the ground, or the compound G*, which can include the snow flux? If the first is the case, and skin temperature is also prescribed, I think you have the risk of violating the energy balance. If the compound flux is prescribed, then this needs to be redistributed between G and S in the model code.

In general, the G = R_n + H_s + H_l that is used in a lot of docstrings is a bit confusing if you factor in snow. It would maybe be good to have separate symbols for 1) the demand from the atmosphere 2) the supply by the ground system with potential snow 3) the always soil only flux

olivierbonte · 2026-09-28T12:03:53Z

+- With [`PrescribedSkinTemperature`](@ref), there is no separate conduction target, so the stored flux *is* the demand, $G = G^\star = R_{\text{net}} + H_s + H_l$.
+- With [`ImplicitSkinTemperature`](@ref), the stored flux is the explicit bare-ground conductive flux evaluated at the current skin temperature, $G = 2\kappa_s (T_g - T_s) / \Delta z_1$. This coincides with the demand only at convergence of the skin temperature solve (see [Skin temperature](@ref "Skin temperature")).
+
+```@example ghf


why is this ghf needed?

olivierbonte · 2026-09-28T12:08:03Z

+DiagnosedGroundHeatFlux
+```
+
+[`DiagnosedGroundHeatFlux`](@ref) closes the surface energy balance internally. What is stored depends on the accompanying skin temperature scheme:


reading the docs first (before looking at the code) it is unclear what this stored means. Maybe this is logical from a code persepctive, but I can't follow what this means from a process point of view

olivierbonte · 2026-09-28T12:10:31Z

+PrescribedGroundHeatFlux
+```
+
+[`PrescribedGroundHeatFlux`](@ref) declares `ground_heat_flux` as an *input* variable and computes nothing. Whatever is written into the field, typically by an external coupler that owns the atmosphere-land interface, survives untouched through `compute_auxiliary!` and reaches the ground boundary condition. This is the configuration that makes [`PrescribedSurfaceEnergyBalance`](@ref) possible; see [Surface energy balance](@ref surface_energy_balance_docs) for the supported combinations.


Suggested change

[`PrescribedGroundHeatFlux`](@ref) declares `ground_heat_flux` as an *input* variable and computes nothing. Whatever is written into the field, typically by an external coupler that owns the atmosphere-land interface, survives untouched through `compute_auxiliary!` and reaches the ground boundary condition. This is the configuration that makes [`PrescribedSurfaceEnergyBalance`](@ref) possible; see [Surface energy balance](@ref surface_energy_balance_docs) for the supported combinations.

[`PrescribedGroundHeatFlux`](@ref) declares `ground_heat_flux` as an *input* variable and computes nothing. Whatever is written into the field (e.g. by a land-atmosphere coupler), is not modified by `compute_auxiliary!`. This is the configuration that makes [`PrescribedSurfaceEnergyBalance`](@ref) possible; see [Surface energy balance](@ref surface_energy_balance_docs) for the supported combinations.

olivierbonte · 2026-09-28T12:14:03Z

+```math
+G^\star = R_{\text{net}} + H_s + H_l
+```
+where $R_{\text{net}}$ is the net radiation budget, $H_s$ is the sensible heat flux, and $H_l$ is the latent heat flux. How that demand is realized, and whether it is imposed at all, is the responsibility of the [`AbstractGroundHeatFlux`](@ref) sub-process.


Suggested change

where $R_{\text{net}}$ is the net radiation budget, $H_s$ is the sensible heat flux, and $H_l$ is the latent heat flux. How that demand is realized, and whether it is imposed at all, is the responsibility of the [`AbstractGroundHeatFlux`](@ref) sub-process.

where $G^\star$ is the ground heat flux *demanded* by the SEB, $R_{\text{net}}$ is the net radiation budget, $H_s$ is the sensible heat flux, and $H_l$ is the latent heat flux. How that demand is realized, and whether it is imposed at all, is the responsibility of the [`AbstractGroundHeatFlux`](@ref) sub-process.

olivierbonte · 2026-09-28T12:16:07Z

+```
+
+!!! warning "Over-determination"
+    Nothing enforces $G = R_{\text{net}} + H_s + H_l$ when the ground heat flux is prescribed alongside prescribed radiative and turbulent fluxes. It is the caller's responsibility to supply a mutually consistent set; Terrarium will happily integrate an inconsistent one.


Suggested change

Nothing enforces $G = R_{\text{net}} + H_s + H_l$ when the ground heat flux is prescribed alongside prescribed radiative and turbulent fluxes. It is the caller's responsibility to supply a mutually consistent set; Terrarium will happily integrate an inconsistent one.

Nothing enforces $G = R_{\text{net}} + H_s + H_l$ when the ground heat flux is prescribed alongside prescribed radiative and turbulent fluxes. It is the user's responsibility to supply a mutually consistent set.

olivierbonte · 2026-09-28T14:08:36Z

+
+Ground heat flux `G` supplied externally as an input variable, for example assembled by a coupler from
+the atmosphere-side surface energy budget. Nothing is computed: `compute_ground_heat_flux!` is a no-op,
+so whatever was written into the `ground_heat_flux` field survives to the ground (soil) top boundary


Suggested change

so whatever was written into the `ground_heat_flux` field survives to the ground (soil) top boundary

so whatever was written into the `ground_heat_flux` is passed on directly the ground (soil) top boundary

olivierbonte · 2026-09-28T14:09:23Z

+condition. The same positive-upward convention applies as for [`DiagnosedGroundHeatFlux`](@ref).
+
+Note that with all four surface fluxes prescribed, the residual identity `G = R_net + H_s + H_l` is
+*not* enforced by construction; it is the caller's responsibility to supply a consistent set.


Suggested change

*not* enforced by construction; it is the caller's responsibility to supply a consistent set.

*not* enforced by construction; it is the user's responsibility to supply a consistent set.

olivierbonte · 2026-09-28T14:17:53Z

-atmosphere-side demanded flux `G` (`= R_net + H_s + H_l`): `Ts = Tg − G/(2κg/Δzg)`. This is the no-snow
-special case (`f_snow = 0`) of the snow-aware method below; it is a separate method (rather than a default
-`snow = nothing`) purely so it can skip the unused `snow_thermal_interface`/`snow_cover_fraction` calls.
+atmosphere-side demanded flux `G` (see [`compute_ground_heat_flux_demand`](@ref)): `Ts = Tg − G/(2κg/Δzg)`.


Suggested change

atmosphere-side demanded flux `G` (see [`compute_ground_heat_flux_demand`](@ref)): `Ts = Tg − G/(2κg/Δzg)`.

atmosphere-side demanded flux `G*` (see [`compute_ground_heat_flux_demand`](@ref)): `Ts = Tg − G*/(2κg/Δzg)`.

olivierbonte · 2026-09-28T14:24:40Z

+
+[`DiagnosedGroundHeatFlux`](@ref) closes the surface energy balance internally. What is stored depends on the accompanying skin temperature scheme:
+
+- With [`PrescribedSkinTemperature`](@ref), there is no separate conduction target, so the stored flux *is* the demand, $G = G^\star = R_{\text{net}} + H_s + H_l$.


But what happens with the snow in this case? Because we are violating the energy balance if this written G is fixed and not altered by the S flux

olivierbonte · 2026-09-28T14:36:40Z

+land model configuration: an implicitly solved skin temperature with all surface fluxes diagnosed
+internally.
+
+Not every combination of `skin_temperature` and `ground_heat_flux` is well posed; the supported


replace by "valid"? I also don't like the term well posed in this context, I think it should be reserverd for the mathematical context where this terms orginiates from

Extract the ground heat flux from `AbstractSkinTemperature` into its own surface-energy-balance sub-process, `AbstractGroundHeatFlux`, with `Diagnosed` and `Prescribed` implementations mirroring the existing pairs for the radiative and turbulent fluxes. `PrescribedSurfaceEnergyBalance` becomes a type alias for a `SurfaceEnergyBalance` whose four flux sub-processes are all prescribed, plus a convenience constructor.

Co-authored-by: Olivier Bonte <65309133+olivierbonte@users.noreply.github.com>

bgroenks96 force-pushed the bg/var-domains branch from 5b8bda6 to 2d832d4 Compare September 25, 2026 13:11

bgroenks96 force-pushed the bg/prescribed-seb branch from c7b6156 to 49d9ce0 Compare September 25, 2026 13:50

bgroenks96 marked this pull request as ready for review September 25, 2026 15:17

bgroenks96 requested a review from olivierbonte September 25, 2026 15:17

bgroenks96 changed the title ~~Refactor ground heat flux into seprate subprocess of SEB and add PrescribedGroundHeatFlux~~ Refactor ground heat flux into separate subprocess of SEB and add PrescribedGroundHeatFlux Sep 25, 2026

bgroenks96 commented Sep 25, 2026

View reviewed changes

Comment thread src/processes/surface/skin_temperature.jl

bgroenks96 commented Sep 25, 2026

View reviewed changes

olivierbonte reviewed Sep 28, 2026

View reviewed changes

Comment thread docs/src/processes/surface_energy/ground_heat_flux.md Outdated

olivierbonte reviewed Sep 28, 2026

View reviewed changes

olivierbonte mentioned this pull request Sep 29, 2026

Add plan (documentation only) force-restore 3-source hydrology model #201

Open

bgroenks96 and others added 5 commits October 1, 2026 18:57

Initial version of plan for prescribed SEB

ae59800

Add worktrees subdir to gitignore

0b12e24

Update implementation plan

9781a06

Update docs/src/processes/surface_energy/ground_heat_flux.md

5909761

Co-authored-by: Olivier Bonte <65309133+olivierbonte@users.noreply.github.com>

bgroenks96 force-pushed the bg/prescribed-seb branch from 70aad11 to 5909761 Compare October 1, 2026 16:58

	where $R_{\text{net}}$ is the net radiation budget, $H_s$ is the sensible heat flux, and $H_l$ is the latent heat flux. How that demand is realized, and whether it is imposed at all, is the responsibility of the [`AbstractGroundHeatFlux`](@ref) sub-process.
	where $G^\star$ is the ground heat flux demanded by the SEB, $R_{\text{net}}$ is the net radiation budget, $H_s$ is the sensible heat flux, and $H_l$ is the latent heat flux. How that demand is realized, and whether it is imposed at all, is the responsibility of the [`AbstractGroundHeatFlux`](@ref) sub-process.

	Nothing enforces $G = R_{\text{net}} + H_s + H_l$ when the ground heat flux is prescribed alongside prescribed radiative and turbulent fluxes. It is the caller's responsibility to supply a mutually consistent set; Terrarium will happily integrate an inconsistent one.
	Nothing enforces $G = R_{\text{net}} + H_s + H_l$ when the ground heat flux is prescribed alongside prescribed radiative and turbulent fluxes. It is the user's responsibility to supply a mutually consistent set.

	so whatever was written into the `ground_heat_flux` field survives to the ground (soil) top boundary
	so whatever was written into the `ground_heat_flux` is passed on directly the ground (soil) top boundary

	not enforced by construction; it is the caller's responsibility to supply a consistent set.
	not enforced by construction; it is the user's responsibility to supply a consistent set.

	atmosphere-side demanded flux `G` (see [`compute_ground_heat_flux_demand`](@ref)): `Ts = Tg − G/(2κg/Δzg)`.
	atmosphere-side demanded flux `G` (see [`compute_ground_heat_flux_demand`](@ref)): `Ts = Tg − G/(2κg/Δzg)`.


		[`DiagnosedGroundHeatFlux`](@ref) closes the surface energy balance internally. What is stored depends on the accompanying skin temperature scheme:

		- With [`PrescribedSkinTemperature`](@ref), there is no separate conduction target, so the stored flux is the demand, $G = G^\star = R_{\text{net}} + H_s + H_l$.

Conversation

bgroenks96 commented Sep 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Why

What changed

Note for reviewers

Verification

Out of scope

Uh oh!

bgroenks96 commented Sep 25, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

olivierbonte left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

bgroenks96 commented Sep 25, 2026 •

edited

Loading