Conservative variably saturated SlabLand + evaporation-front humidity#323
Conservative variably saturated SlabLand + evaporation-front humidity#323glwagner wants to merge 47 commits into
Conversation
…land Scaffolding for the new `VariablySaturatedBucketHydrology` closure that follows in a later commit: - `VanGenuchtenRetention`, `VanGenuchtenConductivity` — retention curve Π_m(𝒮) and Mualem hydraulic conductivity K(𝒮). - `NoDeepLiquidFlux`, `FreeDrainageFlux`, `DarcyDeepLiquidFlux`, `LinearReservoirDrainage` — bottom-boundary liquid flux closures, all positive upward. - `NoRunoff`, `InfiltrationCapacityRunoff` — surface- and subsurface- runoff diagnostics. Runoff is a storage/input bookkeeping export, not an atmosphere-facing flux. All closures are pure @inline functions of scalar state, callable from per-cell kernels. No allocations, type-stable, fully GPU-ready. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Conservative slab hydrology replacing BucketHydrology's clamped P − E
update with a signed-flux budget
dMˡᵃ/dt = Jˡ_b − Jˡ_s − Jᵛ − Rᴹ_lat
driven by the new flux accumulators (`vapor_flux`,
`liquid_precipitation_flux`) and the deep-flux/runoff sub-closures from
the previous commit. Storage variable is the augmented liquid fraction
`ϑˡ = θˡ + Sₛ·max(Π,0)`, so `Mˡᵃ > Mˡᵃ⁺` is admitted and corresponds to
saturated positive-pressure storage.
Adds a closure-extensible diagnostics slot on SlabLand:
* `diagnostic_variables(closure)` / `initial_diagnostic(closure, name, grid)`
on AbstractHydrology and AbstractEnergyBalance.
* `SlabLand.diagnostics::NamedTuple` of closure-owned Fields. With legacy
closures (BucketHydrology + SlabEnergy/ForceRestoreEnergy) the tuple is
empty; with VariablySaturatedBucketHydrology it carries deep_liquid_flux,
surface_liquid_flux, surface_runoff, subsurface_runoff, and
water_storage_tendency (consumed next by the new energy closure).
Verified that a single Jᵛ = 0.01 kg m⁻² s⁻¹ step gives bit-exact
dM/dt = -0.01 and consistent saturation diagnostic.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Force-restore energy closure with M-dependent heat capacity and conservative water-energy advection. Replaces ForceRestoreEnergy's ∂T/∂t = Q/C + (Tᵈ − T)/τ with dEˡᵃ/dt = Λᵈᵉᵉᵖ(Tᵈᵉᵉᵖ − Tˡᵃ) + advective deep − Jᴱ_s − cˡ(Tˡᵃ−Tᵣ) Rᴹ_lat dTˡᵃ/dt = (dEˡᵃ/dt − cˡ(Tˡᵃ − Tᵣ) dMˡᵃ/dt) / C(Mˡᵃ) with C(Mˡᵃ) = C_dry + cˡ Mˡᵃ updated every step. The conservative correction guarantees that adding or removing water at the slab temperature leaves Tˡᵃ unchanged — verified numerically (drainage at T_g over 1000 s gives ΔT = 0 bit-exact). Reads `water_storage_tendency`, `deep_liquid_flux`, `surface_liquid_flux`, and `subsurface_runoff` from `land.diagnostics` (published by the hydrology step) and `surface_energy_flux`, `liquid_precipitation_temperature` from `land.fluxes` (the latter written by the interface in a later commit). Caller supplies *exactly one* of `deep_conductance` (W m⁻² K⁻¹) or `deep_time_scale` (s); supplying both or neither raises ArgumentError. With `deep_time_scale` the conductance is C(Mˡᵃ)/τᵈᵉᵉᵖ — same convention as ForceRestoreEnergy. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
A new humidity formulation alongside SkinHumidity that solves the
atmosphere-facing specific humidity qⁱⁿ from a dry-layer vapor balance
against an unresolved evaporation front at diagnostic depth δᵛ(𝒮):
Gᵉ(qᵉ − qⁱⁿ) = -ρᵃᵗ u★ q★ (Jᵉ = Jᵃ)
with Gᵉ = ρᵃᵗ Dᵛ_eff/max(δᵛ, δᵛ_min) and source humidity qᵉ at the
*front* temperature Tᵉ = Tⁱⁿ + χ(Tˡᵃ − Tⁱⁿ), χ = clip(δᵛ/ℓᵀ, 0, 1). The
wet branch δᵛ ≤ δᵛ_min collapses to qⁱⁿ = qᵛ⁺(Tⁱⁿ), so the saturated-
skin limit reproduces the existing similarity behavior exactly.
The formulation plugs into compute_interface_state via the same Δq-
multiplied SkinHumidity solver pattern (kept finite as Δq → 0). It pairs
with SkinTemperature(DiffusiveFlux(δ = ℓᵀ, κ = κᵀ)) on the temperature
side; Λⁱⁿ = κᵀ/ℓᵀ couples Tˡᵃ to Tⁱⁿ via the existing energy fixed point.
Sub-closures added:
* StorageBasedEvaporationFrontDepth — δᵛ(𝒮) diagnostic
* DryLayerVaporPistonVelocity — Dᵛ_eff with :constant or
:millington_quirk tortuosity
* UnitWaterActivity — aᵉ ≡ 1 (Kelvin activity deferred)
Wired into atmosphere_land_fluxes.jl: EvaporationFrontHumidity pulls
saturation from the materialized hydrology state and bulk temperature
from the energy state (same hooks SkinHumidity uses).
Verified bit-exact against the analytic vapor divider for dry and wet
branches.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Covers the §8 acceptance tests from the PR plan: VariablySaturatedBucketHydrology * augmented-storage diagnostics: 𝒮 at M = 200, 400, 401 with M⁺ = 400 * no-flux conservation: M(t) = M₀ exactly over 10 steps * constant evaporation: dM/dt = -Jᵛ bit-exact * infiltration-capacity runoff: Pˡ < cap fully infiltrates, Pˡ > cap produces R_sfc = Pˡ − cap and storage gains only cap * free drainage at saturation: dM/dt = -ρˡ K_sat WaterCoupledForceRestoreEnergy * analytic force-restore relaxation T(t) = Tᵈ + (T₀−Tᵈ)exp(−Λt/C) matches forward-Euler to 1e-3 K over 100 steps with ΛΔt/C = 5e-4 * advective drainage at slab temperature leaves T invariant to 1e-10 K while M decreases (the conservative cˡ(T−Tᵣ) dM/dt correction) EvaporationFrontHumidity * wet limit (𝒮 ≥ 𝒮ᶜ): qⁱⁿ = qᵛ⁺(Tⁱⁿ) bit-exact * vapor divider: qⁱⁿ matches (Gᵉ qᵉ Δq + Jᵃ qᵃᵗ)/(Gᵉ Δq + Jᵃ) to ε_machine * Tᵉ interpolation: dry case (𝒮 = 0) gives qⁱⁿ < qᵛ⁺(Tⁱⁿ) because Tᵉ < Tⁱⁿ * Gᵉ → 0 (vanishing diffusivity): qⁱⁿ → qᵃᵗ All existing SlabLand tests (35) continue to pass — the diagnostics slot is backward-compatible with the legacy BucketHydrology/SlabEnergy configuration. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Extends `update_net_fluxes!` in slab_land.jl to populate the new signed flux fields (`vapor_flux`, `surface_energy_flux`, `liquid_precipitation_flux`) alongside the legacy positive-part fields (`evaporation`, `precipitation`, `net_energy_flux`). The kernel emits all six in one pass; downstream closures pull whichever subset they declared via `flux_variables(closure)`. Backward compatibility is preserved — closures that don't declare the new keys see zero allocator overhead and the legacy fields keep their existing semantics. Adds `test/test_breeze_evaporation_front_land.jl`: 16 assertions over a tiny coupled Breeze + variably-saturated SlabLand setup. Verifies that * the model constructs (atmosphere, land, evaporation-front interface all wired correctly); * time stepping advances without NaNs; * the new signed `vapor_flux` is populated and nonzero in evaporative conditions; * drydown over 100 short steps measurably reduces `Mˡᵃ`. Uses the default `BulkTemperature()` for the interface temperature since the existing `SkinTemperature(DiffusiveFlux)` solver expects ocean-side `reference_density`/`heat_capacity` properties — adding land-side support is deferred per the PR plan's §13. With `Tⁱⁿ = Tˡᵃ`, the χ-interpolation in the humidity formulation degrades gracefully (Tᵉ = Tˡᵃ) and the wet/dry contrast comes from variable dry-layer piston velocity `Gᵉ = ρᵃᵗ Dᵛ_eff / max(δᵛ, δᵛ_min)`. Pre-existing `test_breeze_coupling.jl` (14 assertions) continues to pass unchanged. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
* `evaporation_front_drydown.jl` — standalone SlabLand with prescribed constant Jᵛ; verifies the conservative water budget produces the analytic M(t) = M₀ − Jᵛ t straight line over 60 days. Shows the saturation diagnostic dropping below 𝒮ᶜ and the bulk T relaxing toward Tᵈᵉᵉᵖ. * `diurnal_force_restore_slab_land.jl` — standalone two-column slab (wet vs dry) under sinusoidal diurnal forcing. Verifies that WaterCoupledForceRestoreEnergy's Mˡᵃ-dependent C(Mˡᵃ) produces the expected wet/dry diurnal-amplitude ratio of 2.12 (= C_wet/C_dry = (C_dry + cˡ M_wet)/C_dry) bit-exactly. * `breeze_over_evaporation_front_slab_land.jl` — 2D Breeze atmosphere coupled to a wet-center/dry-edge SlabLand via EvaporationFrontHumidity. Demonstrates the headline wet/dry contrast: the wet patch evaporates ~112× faster than the dry edge after 60 short steps, with no β(𝒮) scaling prescribed — the contrast emerges from variable δᵛ(𝒮) and the dry-layer piston velocity wᵈ = Dᵛ_eff/max(δᵛ, δᵛ_min). Each example runs in seconds on CPU. The standalone examples are independent of Breeze and so can run as `@example` blocks in the docs; the coupled example exercises the full atmosphere ↔ interface ↔ hydrology ↔ energy stack end-to-end. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
* New docs page `docs/src/land/evaporation_front_slab_land.md` walks through the 10 topics from the PR plan's §10: state variables, Tˡᵃ-vs-Tⁱⁿ distinction, augmented ϑˡ vs physical θˡ, conservative budgets, positive-upward flux convention with the legacy net_energy_flux exception, runoff bookkeeping, the vapor-balance derivation of qⁱⁿ, parameter-choice guidance, a minimal target configuration, and the out-of-scope deferrals. * New "Variably-saturated slab land" subsection in notation.md registers the 20 new symbols (D, ν, ϑˡ, θˡ, θʳ, Sₛ, Π, h, K, κᵀ, Λᵈᵉᵉᵖ, T_r, Tⁱⁿ, qⁱⁿ, Tᵉ, qᵉ, δᵛ, χ, η, ℓᵀ, Dᵛ, wᵈ) with code bindings and units. * Three new Literate examples registered in `docs/make.jl` so they build under Documenter, plus a "Land" page section linking to the new docs page. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
The existing direct-call site in test_slab_land.jl uses the old (Q, P, E, interface_fluxes, Jʳⁿ) positional signature; the kernel now also takes the new signed Jv/Es/Pl outputs (added in 46f33f9 for new-closure consumers). Pass `nothing` at each new position to preserve the existing behavior verifying the legacy fields. All 50+ land-side tests (test_slab_land + new test_variably_saturated + test_water_coupled + test_evaporation_front + test_breeze_evaporation + test_breeze_coupling) pass. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Per review: instead of adding three new Literate examples, convert the existing `breeze_over_slab_land.jl` and `era5_forced_slab_land.jl` to exercise the new closures. * `breeze_over_slab_land.jl` — `BucketHydrology`/`SlabEnergy` replaced with `VariablySaturatedBucketHydrology` + `WaterCoupledForceRestoreEnergy`, and the `FractionalHumidity(CriticalSaturation(0.75))` humidity formulation replaced with `EvaporationFrontHumidity`. The wet-Gaussian / dry-edge contrast that drove the sea-breeze circulation now emerges from the dry-layer piston velocity `wᵈ = Dᵛ_eff / max(δᵛ, δᵛ_min)` with `δᵛ` diagnostic of saturation, rather than from a prescribed `β(𝒮)` scaling. * `era5_forced_slab_land.jl` — same closure swap, plus an explicit `EvaporationFrontHumidity` interface so the Yellowstone domain exercises the new humidity formulation against an actual atmosphere. Removes the three placeholder examples from the previous commit (`evaporation_front_drydown.jl`, `diurnal_force_restore_slab_land.jl`, `breeze_over_evaporation_front_slab_land.jl`) and their docs/make.jl registrations. Updates the land docs page to point at the two converted examples instead. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
`KelvinWaterActivity` was a confusing name — "Kelvin" referred to the Kelvin equation (William Thomson, the temperature-unit Kelvin) for vapor-pressure reduction over a curved meniscus, not the temperature scale. Rename the deferred follow-up closure to `MatricPotentialActivity`, which directly describes the physics: matric-suction-driven `aᵉ < 1` via `aᵉ = exp(g Π / (Rᵛ Tᵉ))`. The Kelvin equation attribution stays in the docstring for citation. Updates the PR plan §2 and §13, the docs page's "Out of scope" list, and the EvaporationFrontHumidity docstring. No code change — the only shipped activity model is still `UnitWaterActivity`. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Lifts the load-bearing design material from the working-tree PR plan into the user-facing docs: * §2 walks through the continuous PDE balances and the augmented- liquid-fraction depth integration that gives the slab its single- prognostic conservative storage variable (was only in the plan). * §3 splits surface vs subsurface runoff bookkeeping. * §4 derives the dry-layer vapor balance and the Δq-multiplied Picard step, including the χ interpolation for the evaporation-front source temperature `Tᵉ` (Picard map was only in the plan). * §5 collects numerical-robustness notes — wet/dry branch handling, bounded-humidity clamp, Sₛ stiffness, convergence tolerances (was only in the plan). * §6 promotes the parameter-choice table; §7 keeps the worked configuration; §8 links the two coupled examples; §9 lists the follow-up deferrals. The plan .md in the repo root is no longer load-bearing — it's now working notes that can stay untracked or be discarded. Menu title in `docs/make.jl` updated to "SlabLand tutorial" so it appears clearly as a tutorial rather than an API-reference page. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Five fixes from the code review: 1. **Explicit imports in `evaporation_front_humidity.jl`** — file used `AtmosphericThermodynamics.*` and `Oceananigans.defaults.FloatType` without an explicit `using` at the top, relying on the include order in `InterfaceComputations.jl`. Now matches the explicit-imports convention used by `interface_states.jl`. 2. **Conservative water-storage tendency** — the hydrology kernel applied a positivity floor (`Mnew = max(M + Δt·dMdt, 0)`) but wrote the *pre-floor* `dMdt` into `water_storage_tendency`. The energy step reads that diagnostic for its conservative `cˡ(T−Tᵣ)·dM/dt` correction, so the two were inconsistent when the floor activated. Now records `(Mnew − Mij)/Δt`, with a comment documenting that the floor silently destroys a sliver of mass on activation (positivity-preserving but not strictly conservative — a known limitation). 3. **Tortuosity-model singleton dispatch** — replaced `tortuosity_model::Symbol` (which triggered `if v.tortuosity_model === :millington_quirk` in the kernel call path) with two singleton types `ConstantTortuosity()` and `MillingtonQuirk()`. The kernel now dispatches at compile time, matching the AGENTS.md rule against short-circuiting `if`/`else` inside kernels. Updated call sites in examples and tests. 4. **Guard `advect_surface_liquid_energy=true`** — the surface advective term reads `land.fluxes.liquid_precipitation_temperature`, but the atmosphere–land flux assembly does not yet write that field. Until the precip temperature is plumbed from the atmospheric state, the constructor throws an `ArgumentError` with a pointer to the followup. The default remains `false`, so this only triggers when a user opts into an incomplete code path. 5. **Removed working-tree planning .md files** — the three plans that sat untracked in the repo root (their content is now in the SlabLand docs tutorial) have been deleted. All existing tests (50 land-side assertions + 14 breeze-coupling assertions + 16 new coupled-Breeze assertions) continue to pass. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
The example body was already using `VariablySaturatedBucketHydrology + WaterCoupledForceRestoreEnergy + EvaporationFrontHumidity` (per the earlier example-conversion commit), but the file's top-level docstring still advertised the old `SlabEnergy + BucketHydrology` configuration. Update the header so the tutorial intro matches the code, and link to the SlabLand docs page for the derivation. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is 📢 Thoughts on this report? Let us know! |
Contributor-facing planning page for everything after PR #323: * Central decision rule — `Metadata` is for satellite remote-sensing products only; reanalysis (ERA5), hybrid products (SoilGrids, GLEAM), and tower observations (FLUXNET) get their own abstractions. Learned maps need training provenance, not product metadata. * Tiered roadmap with effort estimates: Tier 0 (finish PR #323 — 200 LOC), Tier 1 (property providers — 400 LOC), Tier 2 (satellite adapters: MCD43 / MOD21 / GEDI / SoilGrids / etc.), Tier 3 (observation operators: LST, microwave, ET, TWS, towers), Tier 4 (training infrastructure). * Property-by-property classification of every closure parameter as "satellite Metadata?" / "implementation?" with current status. * Two refinements over the earlier strategy doc: - Canopy-height-driven roughness is a *roughness* concern (one-line formula `ℓ_m ≈ a h_c` consuming GEDI canopy height), not a vegetation closure. Belongs in `SimilarityTheoryFluxes` as a property provider. - Aerodynamic roughness, albedo, and emissivity all live outside `SlabLand` — in the flux closure and the radiation module respectively. * Architectural invariants `SlabLand` must keep ruling out: product names in kernels, hard upper clamps on `Mˡᵃ`, observations as state. * Per-tier acceptance criteria. * Out-of-scope list with rationale. Renders under "Land" → "Follow-up roadmap" alongside the SlabLand tutorial. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Closure renames: - VariablySaturatedBucketHydrology -> VariablySaturatedHydrology - WaterCoupledForceRestoreEnergy -> WaterCoupledEnergy (files, includes, exports, tests, docs; drop redundant "bucket" prose where it described the new closure, keep it only for the legacy Manabe BucketHydrology) Fix radiation->land coupling (apply_air_land_radiative_fluxes.jl): route the radiative flux to whichever accumulator the land energy closure declares -- legacy net_energy_flux (+ΣQ_rad, positive into slab) or the new surface_energy_flux (-ΣQ_rad, positive upward). Previously hard-coded net_energy_flux, which crashed any SlabLand using WaterCoupledEnergy with radiation (both converted examples). Closes roadmap Tier-0 gap #10. Examples: fix stale flux-field names in era5_forced_slab_land.jl (net_energy_flux/evaporation/precipitation -> surface_energy_flux/ vapor_flux/liquid_precipitation_flux) and maximum_water_storage -> Mˡᵃ⁺ in breeze_over_slab_land.jl. Both now run end-to-end (verified: full sims + animations). Work around #326 in breeze_over_slab_land.jl: WaterCoupledEnergy's finite deep restoring destabilizes the two-way-coupled LES (runaway near-surface cooling ~day 3). Set deep_conductance = 0 (no deep restoring); documented in follow_up_roadmap.md Tier 0. Not a timestep/heat-capacity/restoring-stiffness issue -- implicit treatment doesn't help (Δt/τ ≈ 1.4e-4); it's a coupled surface-flux feedback (prescribed-atmosphere ERA5 is stable with the same closures and restoring). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Root-caused the WaterCoupledEnergy coupled-LES instability (#326): it is the deep-restoring *target temperature*, not the conductance or numerics. The runaway starts at the dry edges (smallest heat capacity); `deep_temperature = 290` sat ~30 K below the dry surface's natural daytime equilibrium (~320 K), and holding the thin low-C surface that far below equilibrium creates a marginally-unstable cold-surface/warm-air regime that a chaotic LES gust trips into a 2Δz runaway. Restoring on with `deep_temperature = 315` (same Λᵈ) is stable, confirming it is the target. Replace the `deep_conductance = 0` workaround with `deep_temperature = 310` + `deep_time_scale = 12hours`: keeps the (physical) deep restoring rather than disabling it, and runs to completion (verified; land settles ~300-305 K). Roadmap "Known issue" entry rewritten with the root cause and the discriminating runs. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The CI docs build (Julia 1.12) failed building breeze_over_slab_land.jl: the
Float32 atmosphere Simulation's clock could not be reconciled with the
EarthSystemModel's default Float64 clock ("the simulation clock tracks time as
Float32 but the EarthSystemModel clock uses Float64"). Pass a matching
`clock = Clock{eltype(grid)}(time = 0)` to AtmosphereLandModel (forwarded to
EarthSystemModel). The Float64 era5 example is unaffected.
Surfaced only in CI: local verification ran on Julia 1.11 (different resolved
deps) and did not trigger the clock-consistency check.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Roadmap Tier-0 #10: the radiation→land coupling path was previously untested — the existing breeze coupling test passes no `radiation`, so `apply_air_land_radiative_fluxes!` was a no-op there. Add a regression test: build the coupled SlabLand (WaterCoupledEnergy ⇒ `surface_energy_flux`) with a lightweight `PrescribedRadiation` carrying known downwelling SW/LW, and assert that after `update_state!` the radiative flux is *added* to the turbulent `surface_energy_flux` rather than overwriting it (or being overwritten): Es_with_rad ≈ Es_no_rad − ΣQ_rad with ΣQ_rad reconstructed from the kernel's stored diagnostics. Guards both the flux-assembly/radiation call ordering and the accumulator/sign choice the radiation-coupling fix introduced. Verified locally: full file passes (20/20). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
With the breeze clock fix, the docs build now reaches Documenter's linkcheck stage (previously masked by the earlier example failure). Two doc links point at `blob/main/test/test_water_coupled_energy.jl` and `blob/main/test/test_evaporation_front_humidity.jl`, which 404 during this PR's CI because those files only land on `main` at merge. Add a `linkcheck_ignore` regex for the repo's own `blob/main/test/` links; they resolve post-merge. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add `workflow_dispatch` so the docs build can be re-run on demand (e.g. to re-verify after a docs-only fix). Also re-triggers the Documentation workflow, which did not fire on the previous docs-only commit. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…odel # Conflicts: # examples/breeze_over_slab_land.jl # src/Lands/slab_land.jl
The docs build's `install_veros()` step failed at `pixi install` with "failed to solve the pypi requirements" — the conda solve pulled tqdm==4.68.1 while veros `main` pins tqdm==4.67.3 (requirements.txt), which pixi can't reconcile between the conda and pip resolves. Pre-pin tqdm==4.67.3 in conda (available on conda-forge) so both agree, mirroring the existing requests==2.34.0 pin. This failure also affects `main`'s docs build; the fix lands there when this merges. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
|
A few initial comments before going through everything:
Most of what is listed here should be out-of-scope for the For the Van Genuchten retention curve, it might be nice if you can make use of FreezeCurves.jl like Terrarium does, to avoid more code duplication. |
| maximum_front_depth = 0.05, | ||
| critical_saturation = 0.5, | ||
| front_depth_exponent = 2), | ||
| vapor_exchange = DryLayerVaporPistonVelocity( |
There was a problem hiding this comment.
What does "dry-layer piston velocity" mean?
There was a problem hiding this comment.
The piston velocity is the velocity scale U that appears in a bulk formulation like Jᵛ = U Δq where Jᵛ is moisture flux and Δq is the difference between the atmosphere and surface specific humidity.
I believe this flux formulation introduces the notion of a dry surface layer / evaporation front, so that the soil specific humidity is saturated at some depth, separated from the surface by a "dry layer". In this case the surface specific humidity must be computed by balancing fluxes, similar to the way we compute skin temperature through an energy balance.
There was a problem hiding this comment.
Oh, that's cool. Would be great to have something like that in Terrarium too then. Does this surface water balance in the end also include the precipitation flux? Or is that separate? I guess conceptually in my head it's separate, but it depends a bit on what exactly this dry surface layer represents.
There was a problem hiding this comment.
I have to think about this, but I believe it should not depend on precipitation (which is liquid water, whereas this model balances vapor fluxes). That said, if it is precipitating significantly, the surface should become saturated whch would make the surface vapor pressure equal to the saturation specific humidity. This would be good to show in a simple test / example.
There was a problem hiding this comment.
"Surface vapor pressure" here means the vapor pressure of the dry layer? If this is determined in part by the soil moisture content of the uppermost soil layer, then I would agree that this should happen automatically.
| StorageBasedEvaporationFrontDepth(maximum_front_depth, critical_saturation, | ||
| front_depth_exponent) | ||
|
|
||
| Diagnostic evaporation-front depth `δᵛ` as a function of land saturation `𝒮`: |
There was a problem hiding this comment.
The notation for saturation is not included in the table above. It's also a bit confusing to use \mathcal{S} when you use S_s for specific storage. I would suggest using consistent (but distinct) symbols for all volumetric contents (upper bounded by constituent fraction) and saturation fractions (relative to constituent).
There was a problem hiding this comment.
Hmm yes \mathcal{S} should appear in the table. I don't think we should use S_s for specific storage. It conflicts with salinity S.
There was a problem hiding this comment.
L148 of notation.md in this PR:
| ``S_s`` | `specific_storage` | specific storage | Saturated pressure-storage coefficient (m⁻¹) || struct DryLayerVaporPistonVelocity{FT, T} | ||
| minimum_front_depth :: FT | ||
| molecular_diffusivity :: FT | ||
| tortuosity_model :: T |
There was a problem hiding this comment.
Since "model" already has a specific meaning in the Oceananigans/NumericalEarth ecosystem, I would suggest maybe not referring to this as such.
| MillingtonQuirk() | ||
|
|
||
| Millington–Quirk tortuosity: `Dᵛ_eff = Dᵛ₀ · θᵍ^(10/3) / ν²` where | ||
| `θᵍ = ν − θˡ` is the gas-filled pore fraction. Reduces vapor diffusivity in |
There was a problem hiding this comment.
I don't understand what the purpose of this is. Are you trying to explicitly simulate water vapor diffusion between the pore space and the atmosphere? That's a level of detail far beyond slab land model... AFAIK most spatially resolving land/hydrology models don't even do this. Can you provide a reference?
ClimaLand does not do this, as far as I am aware. A quick search of the codebase seems to indicate that tortuosity is only explicitly considered in the CO2 model.
There was a problem hiding this comment.
This is not part of SlabLand though. This is a surface parameterization for computing moisture fluxes between an atmosphere and any land model (including a multilayer land model). This interface parameterization can be used with Terrarium or ClimaLand.
I think we can simplify the model for vapor diffusion. But I don't have an answer for you re: "reference". We don't have references for many things here.
There was a problem hiding this comment.
Ah, ok I missed that. If it's just a general part of the interface flux computation, then that's fine. But I would want to see a reference (or maybe a derivation, if possible) I think to better understand this formulation and what kinds of assumptions it makes. I haven't really seen it before, which is why I am asking.
There was a problem hiding this comment.
I think we should omit "tortuosity" modeling from this PR, and use simple constant vapor diffusion coefficient for this dry layer parameterization, what do you think?
There was a problem hiding this comment.
Here is a reference on different models for the vapor diffusion coefficient in soil:
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011WR011653
There was a problem hiding this comment.
Thanks (I guess Claude...), that helps, especially the Vanderborght et al. reference. I also managed to find the relevant section of the CLM docs which further clarified it.
I guess I am just used to the simpler formulation of the surface humidity flux in terms of a single transfer coefficient, i.e.
There was a problem hiding this comment.
I'm not sure the resistance formulation is consistent with MO theory.
I prepared this comment with the help of Claude but it was not solely due to Claude. I feel weird that you give credit just to Claude here.
There was a problem hiding this comment.
It just read very much like Claude style and not your usual tone, which felt very weird to me.
There was a problem hiding this comment.
Ok, looking at this paper: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1002/2014JD022314
and specifically equation 4, I think that the typical way of formulating this in terms of soil resistances can be consistent with MO theory, provided that
The inconsistency only arises if a "resistance correction" is applied after the solve, which some models do (I believe).
So I think one can think of this DryLayerPistonVelocity as a rewriting of the typical resistance formulation. The advantage is that it reveals the role that a vapor diffusion coefficient has to play in the formulation. It also avoids introducing the "resistance" which could be an advantage or disadvantage depending on your perspective. Basically this way of expressing the concept of a prognostic surface specific humidity is more consistent with how we express the concept of skin temperatures for air-ice and air-sea fluxes.
There was a problem hiding this comment.
It always bothered me that the surface humidity fluxes were handled differently than the ground heat flux in many of these models' formulations, so I'm happy to have something more consistent.
| thermal_exchange_depth :: FT | ||
| porosity :: FT | ||
| water_activity :: A | ||
| phase :: Φ |
There was a problem hiding this comment.
Why is this here? Evaporation is by definition liquid, no? From solid phase it would be sublimation.
There was a problem hiding this comment.
I think it is important that this works for wet or frozen land. We can change the name "evaporation".
There was a problem hiding this comment.
So in principle this should also handle sublimation?
There was a problem hiding this comment.
That would make sense to me. For example, soil with a dry surface layer which is frozen at some depth. I guess in that case the vapor flux is very small, but it would be consistent to treat it in the same way as wet soil with a dry surface layer?
There was a problem hiding this comment.
When you say "dry surface layer", you mean an infinitesimally small layer of solid material at the surface, yes? And you mean that this dry surface layer would be sitting on top of a frozen layer? I think in that case the humidity flux should be calculated in the same way but with different vapor pressures coefficients depending on the skin temperature (frozen vs. unfrozen). See this function in Terrarium.
I guess it should also work the same way for snow too.
There was a problem hiding this comment.
For a snow surface there is no dry layer. The dry layer does not have to be infinitesimally thin, necessarily, but there is an assumption that the layer does not have "memory" (which is why an instantaneous flux balance can be used)
| end | ||
| end | ||
|
|
||
| @inline upwind_deep_internal_energy(Jˡb, Tᵈ, T, Tᵣ, cˡ) = |
There was a problem hiding this comment.
This is something that we're actually missing at the moment in Terrarium, so would be nice to add it there too.
|
One thought that occurred to me while reviewing this is that it might be helpful to adopt some additional standards regarding AI use in drafting issues and PRs. I found myself several times wondering which of the design choices here were consciously made by @glwagner in the prompt/plan vs. which originated from the LLM. This is not so important in cases where the specific formulations or conventions are common or obvious, but it definitely matters when they aren't. I think it would help me at least for this to be made more transparent. I am not sure though what the best way to do this is without creating too much work. Maybe one idea would be to always include your prompt at the top or bottom of the PR? What do you think @glwagner, @simone-silvestri, @maximilian-gelbrecht, et al? |
The pressure-based saturation specific humidity qₛ = εᵈᵛ⁻¹ pᵛ⁺/(p-(1-εᵈᵛ⁻¹)pᵛ⁺) turns negative once the saturation vapor pressure pᵛ⁺ exceeds p/(1-εᵈᵛ⁻¹) ≈ 2.6p (super-boiling interface temperatures). A negative interface humidity drives a runaway spurious-condensation instability in the coupled air-sea fluxes. This arises in practice when an ocean temperature is supplied in Kelvin and read as °C (→ ~566 K). Cap pᵛ⁺ below p so qₛ stays in [0,1); inert in the physical regime where pᵛ⁺ ≪ p. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The main difficulty is that this was not generated by a single prompt, but quite a few prompts over several days and a lot of iteration. What is the additional information that you need to write an effective review? |
Yes, that's what I kind of expected. That is indeed difficult. But maybe there was an initial design or implementation plan?
Basically I need to know what to pay closer attention to. I trust you @glwagner to make reasonable design decisions, check equations carefully, keep track of references, etc. But I don't trust Claude to do that. So if some equations or formulations were invented by Claude without any clear basis, then I think those parts should be subject to greater scrutiny. |
Addresses review feedback on the 𝒮-vs-S_s notation clash by rationalizing the
land height/length family onto a single base symbol `h` distinguished by
superscripts:
- slab depth `D` → `h^{la}` ("depth of prognostic land"); field stays `slab_depth`.
- specific storage `S_s` (m⁻¹) → storage height `h^{ss}` (m): store the
reciprocal length `1/S_s` instead of the inverse-length coefficient. Field
`specific_storage` → `storage_height`; values invert (1e-3 m⁻¹ → 1000 m).
`ϑˡ = θˡ + max(Π,0)/h^{ss}`, saturated `Π = (M−M⁺) h^{ss}/(ρˡ h^{la})`.
- `𝒮` (saturation) kept — the clash is resolved since `S_s` is no longer an "S".
Behavior-neutral (the reciprocal swap is exact); test_variably_saturated_hydrology
passes 10/10. Salinity `Sₛ` (interface_states.jl) and vapor diffusivity `Dᵛ`
are unrelated and untouched.
Updates the hydrology source, notation table, the SlabLand tutorial (equations,
parameter table, reciprocal-aware guidance prose), the follow-up roadmap, both
land examples, and the affected tests.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…odel # Conflicts: # ext/NumericalEarthVerosExt/veros_ocean_simulation.jl
The closure parameterizes phase-agnostic Fickian vapor transport through an
unsaturated surface layer (the "dry surface layer" of soil-evaporation theory),
so "evaporation" — which presumes the liquid phase — is too narrow once the
source at the layer's base can sublimate (frozen soil). Addresses review comment
on evaporation/sublimation, and unifies terminology with the existing
`DryLayerVaporPistonVelocity`.
- `EvaporationFrontHumidity` → `DryLayerHumidity`
- `StorageBasedEvaporationFrontDepth` → `StorageBasedDryLayerDepth`
- fields: `evaporation_front_depth` → `dry_layer_depth`;
`{maximum,minimum}_front_depth` → `{maximum,minimum}_dry_layer_depth`;
`front_depth_exponent` → `dry_layer_exponent`
- `DryLayerVaporPistonVelocity`, `δᵛ` unchanged
- file renames: evaporation_front_humidity.jl → dry_layer_humidity.jl;
test_evaporation_front_humidity.jl → test_dry_layer_humidity.jl;
test_breeze_evaporation_front_land.jl → test_breeze_dry_layer_land.jl
- tutorial page slug kept (avoids doc-URL churn); content updated to "dry layer"
Verified: coupled test passes 20/20.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Two fixes let AtmosphereLandModel/AtmosphereOceanModel couple a Breeze *compressible* terrain-following atmosphere (not just anelastic): - interpolate_state!: get near-surface density and surface pressure via the dynamics-generic dynamics_density(dynamics) / surface_pressure(dynamics) accessors instead of reaching into dynamics.reference_state, which is anelastic-only (nothing for CompressibleDynamics). - surface_layer_height: wrap the zspacing read in @allowscalar; on a terrain-following grid the vertical spacing is a GPU OffsetArray. Validated to iteration 0 on H100 with a coupled terrain LES (correct wet/dry land saturation + cold-air-outbreak surface fluxes). A separate GPU-codegen issue in Breeze's energy-flux BC (𝒬_to_Jᶿ) remains for the full run. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…as NamedTuples - `Oceananigans.prognostic_fields(::SlabLand)` returns the math-named NamedTuple `(; T, Mˡᵃ)` of the prognostics (`saturation` is diagnostic). SlabLand had no field/state accessor before. - `Oceananigans.prognostic_state` / `restore_prognostic_state!(::SlabLand)` make SlabLand checkpointable (it previously had no method, so an ESM-with-SlabLand could not be restored). Verified by an in-memory save/perturb/restore round-trip. - `_water_coupled_step!`: replace the 18 positional args (and the Rule-1-violating `*_field` hybrid names) with `(prognostic, fluxes, diagnostics, energy, …)` — passing `prognostic_fields(land)`, `land.fluxes`, `land.diagnostics` as NamedTuples and reading scalars from `energy`. Self-documenting at both the launch and the body; addresses review comment on water_coupled_energy.jl:137. Behavior-neutral: test_water_coupled_energy passes 3/3. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Locks in the new SlabLand field/state interface: asserts `prognostic_fields` returns the math-named `(; T, Mˡᵃ)` aliasing `temperature`/`water_storage`, and that `prognostic_state` → perturb → `restore_prognostic_state!` recovers the prognostics (deepcopy decouples the snapshot the way on-disk serialization does), plus the `nothing` no-op. Passes 8/8. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The previous `(; T, Mˡᵃ)` mixed a bare and a component-decorated symbol. Use bare symbols within the land model's own namespace — matching Oceananigans (`prognostic_fields(ocean)` is `(; u, v, w, T, S)`, not `Tᵒᶜ`) — and reserve the `ˡᵃ` superscript for cross-component contexts (coupling, the all-components notation table). The `hˡᵃ` slab depth is unaffected: its `ˡᵃ` disambiguates within land (vs hydraulic head `h`), not across components. Updates the kernel read (`prognostic.M`) and the checkpoint test's key assertion. Behavior-neutral; test_water_coupled_energy and test_slab_land pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…_tendency Addresses review comment on water_coupled_energy.jl:173 (separate the tendency calculation from the forward-Euler step). The two land slab-temperature closures (ForceRestoreEnergy — incl. SlabEnergy — and WaterCoupledEnergy) had near-identical "compute ∂T/∂t then T += Δt·∂T/∂t" kernels. Collapse them to: - one shared kernel `_step_land_temperature!` + `step_land_temperature!` launcher (in energy_balance.jl), doing `T += Δt · temperature_tendency(...)`; - a `temperature_tendency(i, j, grid, energy, prognostic, fluxes, diagnostics, time)` generic (kernel-form, (i,j,grid) first), dispatched on the energy closure; - each closure provides its `temperature_tendency` method (the physics) and a one-line `time_step!` delegating to the shared launcher. 2 bespoke step kernels → 1 shared kernel + 2 tendency methods. Behavior-neutral: test_water_coupled_energy 3/3, test_slab_land (incl. Force restore + checkpoint) green. SlabOcean's `_step_slab_temperature!` is the same shape and can join this dispatch in a follow-up (would move the generic to a shared module). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The `#####` banner before the humidity solver announced an "Interface-state hooks" section that no longer has any code under it (the hooks live in atmosphere_land_fluxes.jl), and its description was stale. Remove it; the remaining `#####` banners are the usual Oceananigans file-section dividers. Addresses review comment dry_layer_humidity.jl:226. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…d tutorial Justification (review comments :88/:295): in-soil vapor diffusion only dominates within a thin dry surface layer (Philip & de Vries 1957; Tang & Riley 2013), so land models — including multi-layer Richards models like GEOtop (Endrizzi et al. 2014) and ClimaLand — transport only liquid water in the column and represent the dry-layer limitation as a surface condition (Vanderborght et al. 2017). DryLayerHumidity is that surface condition: the DSL resistance δᵛ/Dᵛ_eff of Yamanaka et al. (1997) and Swenson & Lawrence (2014) (the CLM5 scheme, also in ClimaLand with τ=1 and δ_max = 15 mm), closed by the vapor-flux balance of Ye & Pielke (1993) instead of a prescribed α or β. The series solution is algebraically the same expression ClimaLand evaluates. References now woven into the file banner and docstrings; eleven entries added to NumericalEarth.bib (ye1993.pdf lives at the repo root). Minimality: drop the speculative water-activity hook (UnitWaterActivity always returned 1 from a placeholder Π = 0) and the unused (Tᵉ, pᵃᵗ) arguments of effective_vapor_diffusivity. The Kelvin-equation pore humidity hₛ = exp(gΠ/(RᵛTᵉ)) — Ye & Pielke's hₛ, the α factor of CLM5/ClimaLand — remains a documented follow-up (matters only at extreme dryness). Tutorial: new §4.3 "Where this sits in the literature" (two-stage evaporation, why Richards models neglect in-column vapor, the lineage of each piece) and §7 "The closure menu" surveying every energy / hydrology / interface-humidity closure (FractionalHumidity = β method, SkinHumidity = constant resistance, DryLayerHumidity = DSL) plus PrescribedLand. The assembled example is now a runnable @example (was a plain fence calling the nonexistent tortuosity_model = :millington_quirk), verified locally. Verified: NumericalEarth + Breeze ext precompile; tutorial example runs; test_breeze_dry_layer_land.jl 20/20. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
glwagner
force-pushed
the
glw/systematic-land-model
branch
from
c34c67e to
0f2753e
Compare
The development process is really the outcome of a proposal-review loop, so it isn't really uniformly possible to filter for things that we "can't trust because Claude did it". I think as a result, we should take as close a look as possible to all parts of the review. If possible, I would approach the review more from a practical perspective: quickly review everything that you possibly can, and for things that can't be reviewed by reading the code (and for which validation isn't provided), propose or ask for empirical validation. In terms of checking equations, I think we usually take the approach to merge first and validate over a longer period of time. The process of validation will involve a mix of experiments / simulations, plus reading docs and checking math. |
Addresses review comment :295 ("I don't recognize this formulation"):
spell out, step by step, that the coded Δq-multiplied expression is the
two-conductances-in-series solution of the surface vapor-flux balance
Jᵉ = Jᵃ (eq. 12b of Ye & Pielke 1993 with hₛ = 1), with the similarity
flux linearized as Gᵃ = Jᵃ/Δq about the previous Picard iterate and the
balance multiplied through by Δq to stay finite as Δq → 0. Also records
the Gᵉ → ∞ / Gᵉ → 0 limits and that the fixed point satisfies the true
nonlinear balance.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…odel # Conflicts: # docs/src/NumericalEarth.bib
- §1: liquid_precipitation_flux is positive DOWNWARD (per slab_land.jl); give it its own sign row instead of lumping it under "leaves slab". update_net_fluxes! takes the coupled model, not the land. - §1: "§3.3 below" pointed at a section that does not exist → §2.3. - §3: runoff(land) does not exist; point at the land.diagnostics fields. - §4: repair two sentences garbled by the evaporation-front → dry-layer rename; unify J^soil → Jᵉ with the code and the solver banner; qˢ → qⁱⁿ in the β-method strawman. - §4.2/§5: undefined symbols pᵉ/pⁱⁿ → pᵃᵗ (the code evaluates qᵛ⁺ at atmospheric pressure); step 8 used undefined wᵛ → state the converged flux as Jᵛ = −ρᵃᵗ u★ q★. - §7: define Q, C, τ used in the energy-closure table. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…utput renders An @example block displays only its final value, so the informative SlabLand tree (closure summaries, flux/diagnostic keys) was invisible — only the DryLayerHumidity summary appeared. Split at the land/interface boundary so both render. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…calEarth/NumericalEarth.jl into glw/systematic-land-model
…, not developers
Notation: the component-superscript bullet now states the rule — bare
symbols within a component's own namespace (prognostic_fields(land) =
(; T, M)), the ˡᵃ decoration only where land variables appear alongside
other components. The land table drops the inconsistent M^{la}/M^{la+}
decoration (T was already bare), and the tutorial's §1 explains why the
page itself uses Tˡᵃ/Mˡᵃ (cross-component context throughout).
Tutorial: remove development-history framing — the docs describe the
current implementation instantaneously. "signed (new)"/"legacy
(pre-2026)" flux table → conventions keyed by which closures read them;
"in this PR"/"is a follow-up"/"deferred" phrasing removed throughout;
§10 "Out of scope" (a developer PR-scope list) rewritten as
"Limitations" for users judging whether SlabLand fits their problem,
with developer items (implicit Darcy, two-node solve, unimplemented
closure names) dropped. Same sweep through the dry_layer_humidity and
hydraulic_functions docstrings/banners; the WaterCoupledEnergy error
pointer now targets the roadmap page instead of the renamed §10.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…U kernels
The shared `_step_land_temperature!` kernel receives the energy closure
whole (for `temperature_tendency` dispatch), but `ForceRestoreEnergy` can
hold Field-valued heat capacities and `WaterCoupledEnergy` a Field/
FieldTimeSeries deep temperature. Without `Adapt.adapt_structure` rules
the closure reached the CUDA kernel un-adapted and failed the isbits
check ("SlabLand property providers" on the GPU pipeline). The previous
bespoke kernels dodged this by passing each property as its own argument.
Add adapt rules for both closures (the `ThreeEquationHeatFlux` pattern)
and `using Adapt: Adapt` in Lands.
Verified on an H100 (gpu-prod, GPU_TEST=true): test_slab_land.jl all
green including property providers 5/5 on GPU{CUDABackend}; CPU run also
green (43/43).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Tutorial §4 restructured per review: §4.1 states the physical model on its own — the source humidity is explicitly the saturation specific humidity at the front temperature, qᵉ = qᵛ⁺(Tᵉ, pᵃᵗ); the dry-layer Fick flux balances the Monin–Obukhov flux Jᵃ(Tⁱⁿ, qⁱⁿ), a nonlinear equation for qⁱⁿ. §4.2 then describes the numerical solution: Picard iteration, the linearized conductance Gᵃ = Jᵃ/Δq anchored at the previous iterate qⁱⁿ⁻, and the Δq-multiplied series solution. Previously the model and solver were interleaved (qⁱⁿ⁻ appeared inside the model description). All displayed equations now use \frac instead of inline division (§2 storage/head/saturation displays included). The source banner mirrors the model/solver split. Implementation review findings: - Tutorial §5 claimed a qⁱⁿ clip and solver settings (ε_q = 1e-10, ε_T = 1e-4, 30 iterations, ω = 0.7 under-relaxation) that do not exist. Real solver: drift |Δu★|+|Δθ★|+|Δq★| < 1e-8 (default), maxiter 100, no relaxation, no clamp (same as SkinHumidity). §5 now states the implemented behavior. - compute_interface_humidity: replace the early-return `if` wet branch with a single branchless ifelse return (kernel rules); the wet branch no longer computes an unused front-saturation humidity; the wᵈ intermediate folds into Gᵉ. Behavior unchanged. - Algebra audited: sign conventions, D == 0 guard, and the Millington–Quirk θᵍ → 0 degeneracy is unreachable (requires 𝒮 ≥ 𝒮ᶜ, which takes the wet branch). Struct field alignment tidied. Verified: test_breeze_dry_layer_land 20/20, test_slab_land all green, land docs pages build with citations and KaTeX rendering checked. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The tutorial described two-stage evaporation but never demonstrated it. New §9 runs a single-column SlabLand under a constant warm, dry prescribed atmosphere (305 K, q = 0.004, 3 m/s) with prescribed radiation (SW 250 / LW 350 W m⁻², albedo 0.3, ε = 0.95), drainage and runoff off, for twelve days at Δt = 10 min — seconds of compute in the docs build. The rendered figure shows the stage-1 demand-limited plateau (≈ 3.4 mm/day, slab steady at 298 K under evaporative cooling), the 𝒮ᶜ crossing near day 8, the stage-2 evaporation collapse through the growing dry layer, and the concurrent slab warming to ≈ 304 K as latent heat repartitions to sensible — with no prescribed β(𝒮). The notch at the crossing is called out as §5's hard wet/dry branch switch. Sections renumber: Examples → §10, Limitations → §11 (internal refs updated). Verified: full @example chain runs end-to-end and the figure renders in a local Documenter build. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the fill!(parent(...)) initialization with the public set! API: set!(slab; T, M) for the land, and set!(field, t -> value) for the prescribed-atmosphere and radiation FieldTimeSeries (their spatial dimensions are Flat/Nothing, so the setter takes time alone; a constant in t is steady forcing). This drops the reach into parent(...) and is more compact. Conciseness: reuse the §8 interface_humidity instead of redefining an identical DryLayerHumidity, and split the setup into a slab block and a forcing/coupling block. Clarity: the initial-storage comment now reads M = 0.9 Mˡᵃ⁺ (was an imprecise "𝒮 ≈ 0.9"). Accuracy: the stage-1 steady-state note no longer claims a two-term radiation+sensible balance (the deep restoring toward 295 K also cools at 298 K) — it now says evaporative cooling removes most of the net surface heating. Behavior identical (same 𝒮/E/T trajectory); full @example chain runs and the figure renders in a local Documenter build. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…orial §9 Mirror the existing set!(::SlabLand; T, M) convenience for the two prescribed components so the case study can initialize each model with one top-level call: set!(atmosphere; u, v, T, q, p) set!(radiation; downwelling_shortwave, downwelling_longwave) Their fields are FieldTimeSeries, which only accept set!(fts, ::Function); a helper wraps a Number as a constant (in space and time) via a varargs lambda, forwards functions/Fields unchanged, and treats `nothing` as a no-op. Each setter refreshes the interpolated state with update_state!. Both follow the qualified-Oceananigans.set! pattern already used by SlabLand (no new imports). Tutorial §9 now uses set!(slab; T, M), set!(atmosphere; ...), and set!(radiation; ...) — and builds radiation from the grid constructor instead of allocating FieldTimeSeries by hand. The forcing setup drops from ~20 lines to ~6, off the parent(...) internal and onto the public API. Tests: set! testsets added to test_prescribed_atmosphere.jl (6/6) and test_radiations.jl (4/4), both download-free; verified locally. Tutorial @example chain runs (identical 𝒮/E/T trajectory) and the figure renders. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…field, val) Drop the constant-in-(space,time) lambda wrapper. A prescribed field is a FieldTimeSeries whose data lives in its per-time `Field` slices, and `set!(::Field, value)` already handles a Number, Field, function, or operation — exactly as the Oceananigans/Breeze model `set!` relies on. So the helper just loops the time slices and calls `set!(fts[n], value)`; no manual number broadcast. Behavior identical: tutorial §9 chain reproduces the same 𝒮/E/T trajectory, set! testsets still pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2 participants
Summary
VariablySaturatedBucketHydrology+WaterCoupledForceRestoreEnergyEvaporationFrontHumidityto replace the clamped bucket, theconstant-
Cforce-restore, and the prescribed-β(𝒮)surface humidity.All three plug into the existing
SlabLandandComponentInterfacesscaffolding via per-closure declarations — legacy
BucketHydrology + SlabEnergy + FractionalHumidityconfigurations stay unchanged.ϑˡ = θˡ + Sₛ·max(Π, 0), soMˡᵃ > Mˡᵃ⁺is admitted as positivepressure head rather than clamped. The surface humidity
qⁱⁿis solvedby a dry-layer vapor-flux balance against an unresolved evaporation
front at saturation-dependent depth
δᵛ(𝒮), collapsing exactly to asaturated skin when
𝒮 ≥ 𝒮ᶜ.breeze_over_slab_land.jlandera5_forced_slab_land.jltothe new closures; adds a tutorial-style docs page; extends
notation.mdwith 22 new symbols.See
docs/src/land/evaporation_front_slab_land.mdfor the SlabLandtutorial (continuous balances, depth-integrated reduction, the Picard
fixed-point algorithm, parameter guidance, numerical-robustness notes,
and a worked configuration).
What changed at the closure level
VariablySaturatedBucketHydrologydMˡᵃ/dt = Jˡ_b − Jˡ_s − Jᵛ − Rᴹ_latbudget; augmented liquid fraction storageWaterCoupledForceRestoreEnergyMˡᵃ-dependentC(Mˡᵃ) = C_dry + cˡ Mˡᵃ; conservativedTˡᵃ/dt(adding water at slab T leaves T invariant)VanGenuchtenRetention/VanGenuchtenConductivityΠ(𝒮)and MualemK(𝒮)NoDeepLiquidFlux,FreeDrainageFlux,DarcyDeepLiquidFlux,LinearReservoirDrainageNoRunoff,InfiltrationCapacityRunoffEvaporationFrontHumidityqⁱⁿ; sub-closuresStorageBasedEvaporationFrontDepth,DryLayerVaporPistonVelocity,ConstantTortuosity/MillingtonQuirk,UnitWaterActivitySlabLandgains a closure-extensiblediagnostics::NamedTupleslot andthe flux assembly now populates signed flux fields (
vapor_flux,surface_energy_flux,liquid_precipitation_flux) alongside the legacypositive-part fields (
evaporation,precipitation,net_energy_flux).Backward-compatible — legacy configurations get an empty diagnostics
tuple and only the legacy flux fields are written/read.
Verified analytically
conservation
T(t) = Tᵈᵉᵉᵖ + (T₀ − Tᵈᵉᵉᵖ)·exp(−Λt/C)qⁱⁿ = (Gᵉ qᵉ Δq + Jᵃ qᵃᵗ) / (Gᵉ Δq + Jᵃ)(ε_machine)qⁱⁿ = qᵛ⁺(Tⁱⁿ)JᵛOut of scope (deferred to follow-ups)
Listed in
docs/src/land/evaporation_front_slab_land.md§9:Eˡᵃas a first-class state variable (needed for phase change).MatricPotentialActivity(Kelvin-equation suction-driven vapor reduction).SkinTemperature(DiffusiveFlux)solve soTⁱⁿ ≠ Tˡᵃ.solve, implicit deep Darcy, subgrid tile blending.
Test plan
test/test_slab_land.jl— 35 existing assertionstest/test_breeze_coupling.jl— 14 existing assertionstest/test_variably_saturated_bucket_hydrology.jl— diagnostics +conservation (10)
test/test_water_coupled_force_restore_energy.jl— analyticforce-restore + advective invariance (3)
test/test_evaporation_front_humidity.jl— vapor divider, wetbranch, Tᵉ interpolation, edge cases (5)
test/test_breeze_evaporation_front_land.jl— coupled wiring +drydown (16)
🤖 Generated with Claude Code