SigmaColumn.hpp File Reference

G-steps 4 and 5 — closure B: the σ-coordinate explicit unsaturated column with a moving lower boundary (the water table). More...

#include "SoilCharacteristic.hpp"

Include dependency graph for SigmaColumn.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	openswmm::twoD::sigma::ColumnStep
	What one column sweep needs and what it reports back. More...

Namespaces
namespace	openswmm
namespace	openswmm::twoD
namespace	openswmm::twoD::sigma

Functions
void	openswmm::twoD::sigma::advanceColumn (const soil::Params &p, double *theta, int m, int stride, ColumnStep &st) noexcept
	Advance one σ column by st.dt.
double	openswmm::twoD::sigma::columnDtLimit (const soil::Params &p, const double *theta, int m, int stride, double L, double Ldot, double c_col, bool capillary) noexcept
void	openswmm::twoD::sigma::seedHydrostatic (const soil::Params &p, double *theta, int m, int stride, double L) noexcept
	Seed a column to hydrostatic equilibrium above a table at depth L.
double	openswmm::twoD::sigma::columnStorage (const double *theta, int m, int stride, double L) noexcept
	Total water in the column (m of water) — Σ θ_j·L·Δσ.

Detailed Description

G-steps 4 and 5 — closure B: the σ-coordinate explicit unsaturated column with a moving lower boundary (the water table).

Coordinate and state

The column occupies z ∈ [0, L] measured DOWNWARD from the ground surface, L = z_s − h_g. Map to σ = z/L ∈ [0,1] with a fixed layer count m; layer j owns σ ∈ [j/m, (j+1)/m], Δσ = 1/m. Layers stretch and compress as the table moves — no regridding, no variable layer counts, fixed SoA stride.

The conserved quantity is water depth per layer, w_j = θ_j·L·Δσ (metres of water). Evolving w and deriving θ = w/(L·Δσ) from the NEW L is what makes the compression exact: a column whose faces move but whose fluxes vanish keeps every w_j unchanged, so total water is invariant to machine precision. That is plan gate 5, and it holds by construction here rather than by cancellation.

Fluxes

w_j ← w_j + Δt·( F_{j−1/2} − F_{j+1/2} )          [downward-positive]
 
F_{j+1/2} = K(θ_up)                                  advective, upwind
          − θ_donor·σ_{j+1/2}·L̇                      grid motion (ALE)
          − ( D̄_{j+1/2} / (L·Δσ) )·( θ_{j+1} − θ_j ) optional capillary

• Gravity is downward, so the advective donor is always the layer above.
• The grid term is upwinded on its own sign: with L̇ > 0 (table falling, column growing) the faces sweep downward and the donor is the layer below; with L̇ < 0 it is the layer above.
• D̄ is the harmonic mean of the two layers' diffusivities. Off by default (D-N3).

Boundaries and the handover

• Top (σ=0): F_top = q⁺ − q_ET, with the top layer's acceptance capping q⁺ — the excess is rejected and returned to the caller for the surface state (Hortonian / saturation rejection).
• Bottom (σ=1): F_bot = q₀_phys − θ_bot·L̇. The first term is the physical Darcy flux across the table; the second is the handover — the water swept across the moving boundary. Both together are what the saturated zone gains, which is why SubsurfaceSolver books F_bot and not q₀_phys as the recharge.

Note

Specific yield. Making the handover explicit forces the saturated update's storage coefficient to be θ_s − θ_bot, not θ_s. The plan's §2.1 writes θ_s; θ_s − θ_bot is what actually conserves once the column retains water above a falling table, and it is the textbook specific yield. This is the one place this implementation is deliberately more precise than the plan text, and it is called out again in SubsurfaceSolver::fireCell.

Author

Caleb Buahin caleb.nosp@m..bua.nosp@m.hin@g.nosp@m.mail.nosp@m..com

Copyright

Copyright (c) 2026 Caleb Buahin. All rights reserved.

License\n Apache-2.0