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 <vector>

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::ColumnTransfers
struct	openswmm::twoD::sigma::ColumnStep
	What one column sweep needs and what it reports back. More...
struct	openswmm::twoD::sigma::ColumnPreparation

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.
void	openswmm::twoD::sigma::prepareColumn (const soil::Params &p, const double *theta, int m, int stride, const ColumnStep &st, ColumnPreparation &prepared) noexcept
	Prepare reusable trial invariants. Same input preconditions as advanceColumn; m > 0.
void	openswmm::twoD::sigma::advanceColumn (const soil::Params &p, double *theta, int m, int stride, ColumnStep &st, const ColumnPreparation &prepared) noexcept
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
	The explicit stability bound of a column, without advancing it (§2.3): Δt_u ≤ C_col·min_j[ L·Δσ / ( c(θ_j) + |σ_j·L̇| + 2·D̄_j/(L·Δσ) ) ] with celerity c(θ) = dK/dθ.
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. SubsurfaceSolver settles the table using the actual bounded F_bot, compression surplus and residual support from the same sweep. Its physical recharge diagnostic remains q₀_phys; boundary handover is internal and is excluded from external continuity.

Note

Specific yield. θ_s − θ_bot gives the table predictor. The final SIGMA bracket pairs saturated storage θ_s·h_g with the actual column transfer rather than correcting the table after a column has collapsed to a zero-yield saturated profile.

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

Apache-2.0