FLUXOS

2D Hydrodynamic & Solute Transport Modelling Framework

A high-performance, open-source platform for physics-based flood simulation,
solute transport, and soil infiltration modelling

Diogo Costa

University of Évora, Portugal

What is FLUXOS?

FLUXOS is an open-source, physics-based 2D hydrodynamic model for event-based flood simulation, coupled with solute transport and soil infiltration.

                            Key Innovation:

                            Supports both regular Cartesian and unstructured triangular meshes with full GPU acceleration, enabling fast, high-resolution simulations on complex terrain.

Physics-based — 2D Shallow Water Equations (dynamic wave)
Multi-physics — ADE solute transport + soil infiltration
Scalable — OpenMP, MPI, and CUDA parallelization
Open — GNU GPL v3+ license

┌─────────────────────────────────┐
│        INPUT DATA                │
│  DEM + Boundary Conditions      │
│  + Soil Map + Meteo Forcing     │
└────────────────┬────────────────┘
                │
                ▼
┌─────────────────────────────────┐
│  2D SHALLOW WATER EQUATIONS    │
│  Roe flux + wetting/drying     │
└────────────────┬────────────────┘
                │
        ┌───────┴───────┐
        │               │
        ▼               ▼
┌──────────────┐ ┌────────────────┐
│ ADE SOLVER   │ │ SOIL           │
│ Solute       │ │ INFILTRATION   │
│ Transport    │ │ Horton model   │
└──────────────┘ └────────────────┘

Core Capabilities

🌊

2D Hydrodynamics

Full dynamic wave solution of the 2D shallow water equations. Roe's approximate Riemann solver with robust wetting/drying and back-flooding capture.

🧪

Solute Transport

Mass-conservative ADE solver with first-order upwind advection, Fickian diffusion, and monotonicity clamp. Two-pass fully parallel formulation.

🌱

Soil Infiltration

Horton decay model with USDA 12-class soil texture classification (Rawls et al. 1983). Per-cell infiltration tracking.

⚡

GPU Acceleration

CUDA kernels for both regular and triangular meshes. Typical 10–50x speedup over single-core CPU. Supports Pascal GPUs and newer.

🌐

WebGL 3D Viewer

Interactive browser-based visualization with particle tracking, satellite imagery, concentration overlays, and real-time animation controls.

📊

Multi-Format Output

ASCII grids, VTK/ParaView (.vtu + .pvd), KML/KMZ for Google Earth, and WebGL data export for web hosting.

Coupled Host Models

FLUXOS can operate standalone or coupled with other modelling frameworks:

Integration	Domain	Description
Standalone	Overland Flow	Event-based 2D flood simulation with precipitation or inflow forcing
WINTRA	Nutrient Release	Runoff-soil interaction algorithm for snowmelt-driven nutrient transport

                    Originally developed for Canadian Prairie hydrology (snowmelt flooding, fill-and-spill dynamics), FLUXOS has evolved into a general-purpose 2D overland flow platform applicable to any terrain and climate.
                

01

Mesh & Numerical Methods

Regular Cartesian Mesh

📈 Structured Grid

Uniform cell size (configurable dxy)
Input: ESRI ASCII DEM format (.asc)
Armadillo matrix storage (arma::Mat)
Efficient row-column indexing
Output: ASCII text grids per timestep

                            Best for: Rectangular domains with uniform resolution requirements. Simple DEM-based setups.
                        

🔧 Solver Details

Flux scheme	Roe's approximate Riemann
Time stepping	CFL-adaptive
Wetting/drying	Velocity desingularization
Friction	Manning's equation
Parallelization	OpenMP, MPI, CUDA

🕑 Performance

Rosa_2m test domain (859×618, 2m cells):
~16 minutes on single CPU core

Unstructured Triangular Mesh

🔶 Edge-Based Finite Volume

Arbitrary triangular elements
Input: Gmsh (.msh) or Triangle (.node/.ele)
Flat vector storage (std::vector)
Output: VTK XML (.vtu) + ParaView (.pvd)

🕑 Performance

Same Rosa_2m domain:
~3 seconds (vs 16 min on regular mesh!)

🧮 Advanced Numerics

Reconstruction	MUSCL with least-squares gradients
Limiter	Barth-Jespersen (monotonicity)
Well-balanced	Hydrostatic reconstruction
Dry fronts	Ritter dry-front solution
CFL	Inradius-based adaptive

                            Best for: Complex geometries, channel networks, slope-adaptive refinement, and when VTK/ParaView output is desired.
                        

SWE: Conservative Vector Form

2D Shallow Water Equations — Hyperbolic Conservation Law

∂U/∂t + ∂F(U)/∂x + ∂G(U)/∂y = S

Saint-Venant (1871) · Manning (1891) · Elder (1959) after Boussinesq (1877)

State & Flux Vectors

U = [h, hu, hv]^T

F(U) = [hu, hu² + gh²/2, huv]^T

G(U) = [hv, huv, hv² + gh²/2]^T

Source Terms

S = [R − I, −ghS_fx + τ_xx, −ghS_fy + τ_yy]^T

R − I: rainfall minus infiltration (mass source)
S_f: Manning friction slope = n²|u|u / h^4/3
τ: turbulent stress via eddy viscosity ν_t = c_v u_* h + ν_mol

                    h = water depth  |  u, v = depth-averaged velocities  |  qx = hu, qy = hv (unit discharge)  |  g = 9.81 m/s²  |  n = Manning coefficient
                

Why conservation form?

Ensures mass and momentum are preserved exactly across shocks, wetting fronts, and hydraulic jumps — essential for flood simulation where discontinuities arise naturally.

Finite Volume Discretization

Cell-Centered Godunov Scheme

Godunov (1959); LeVeque (2002)

Domain divided into non-overlapping control volumes (cells)
Fluxes computed at cell faces via Riemann solver
Explicit forward-Euler time integration

Discrete Update Rule

Uⁿ⁺¹_i = Uⁿ_i − (Δt / A_i) ∑_faces F_face · L_face

A_i = cell area, L_face = face length
F_face = numerical flux from Roe solver

CFL Adaptive Time-Stepping

Courant, Friedrichs & Lewy (1928)

Time step is recomputed every iteration as the global minimum over all wet cells:

Regular mesh

Δt = CFL · Δx / (|u| + √(gh))

Triangular mesh

Δt = CFL · r_in / (|u| + c)

r_in = triangle inradius (inscribed circle)
c = √(gh) = wave celerity

Why Godunov-type FVM?

Naturally handles shocks and wetting/drying without artificial viscosity. CFL-adaptive stepping (CFL < 1; typical 0.5) guarantees stability while maximizing computational efficiency.

Roe Approximate Riemann Solver

Roe Flux Structure

Roe (1981); Toro (2001)

F = ½(F_L + F_R) − ½ R|Λ|R⁻¹ · (U_R − U_L)

3 wave speeds: λ₁ = û − ĉ, λ₂ = û, λ₃ = û + ĉ
Roe averages (depth-weighted):

û = (u_L√h_L + u_R√h_R) / (√h_L + √h_R)

ĉ = √(g · (h_L + h_R) / 2)

R, Λ = right eigenvectors & eigenvalues of Roe matrix
Upwind dissipation captures shocks & rarefactions

Hydrostatic Reconstruction

Audusse et al. (2004) — preserves lake-at-rest (C-property)

z_b,face = max(z_b,L, z_b,R)

h_recon = max(0, η − z_b,face)

Reconstructed states fed into Roe solver
Prevents spurious flow over dry bed steps
Pressure source correction at each cell:

S_p = ½ g (h²_orig − h²_recon)

Balances flux and topography gradients exactly

Why Roe + hydrostatic reconstruction?

Roe’s solver resolves all 3 wave families (shocks, rarefactions, shear) with minimal numerical diffusion. Hydrostatic reconstruction preserves the lake-at-rest state exactly, preventing spurious flows over irregular bathymetry.

Wetting/Drying & Bed Friction

Ritter Dry-Front Solution

Ritter (1892)

At wet-dry interfaces, analytical dam-break solution replaces Roe flux:

q_front = 0.296 · Δz · √(g · Δz)

h_front = 0.444 · Δz

Δz = water surface − dry bed elevation
Applied when one cell is wet, neighbor is dry
Volume limiter: flux capped at available water

Velocity Desingularization

Kurganov & Petrova (2007)

u = 2hq / (h² + max(h², ε²))

Avoids division by zero as h → 0. Smoothly transitions to u = 0 in dry areas.

Semi-Implicit Manning Friction

Implicit treatment avoids stiff ODE instability at small depths:

qⁿ⁺¹ = qⁿ / (1 + c_f |u| Δt / h)

c_f = g n² / h^1/3

Applied after the explicit SWE update step
Unconditionally stable (denominator ≥ 1)
Naturally decelerates flow as h → 0
No additional time-step restriction

Why Ritter + semi-implicit friction?

Ritter’s analytical dam-break replaces the Riemann solver at wet-dry fronts, avoiding zero-depth singularities. Semi-implicit Manning friction is unconditionally stable — no stiff time-step restriction at shallow depths. Mass balance limiter guarantees h ≥ 0.

MUSCL Reconstruction & Limiter

Unstructured triangular mesh only — achieves 2nd-order spatial accuracy while preserving monotonicity.

MUSCL Reconstruction

van Leer (1979)

φ_face = φ_i + Φ_i · (∇φ_i · (x_face − x_i))

Gradient ∇φ via weighted least-squares over cell neighbours
Distance-weighted (1/d²) for numerical consistency
LSQ inverse matrix precomputed at mesh setup
Φ_i ∈ [0, 1] = slope limiter (see right)

Edge-Normal Rotation

Rotate (q_x, q_y) into face-normal frame
Solve 1D Riemann problem (Roe/HLL)
Rotate flux back to Cartesian frame

Barth-Jespersen Limiter

Barth & Jespersen (1989)

Φ_i = min_edges min(1, α_e)

α_e = (φ_max − φ_i) / δ_e if δ_e > 0

α_e = (φ_min − φ_i) / δ_e if δ_e < 0

δ_e = reconstructed − cell-center value
φ_max, φ_min = extrema over cell & neighbors
Ensures no new extrema are created (monotonicity)
Φ = 1 → full 2nd order; Φ = 0 → 1st order (near discontinuity)

Why MUSCL + Barth-Jespersen?

MUSCL achieves 2nd-order spatial accuracy on unstructured meshes, critical for resolving flow features on complex terrain. Barth-Jespersen prevents spurious oscillations near discontinuities while preserving maximum accuracy in smooth regions.

ADE Solver: Two-Pass Formulation

Advection-Diffusion Equation ∂(hC)/∂t + ∇·(uhC) = ∇·(Dh∇C) solved with a two-pass explicit scheme:

Fick (1855); Leonard (1991)

Pass 1: Face Fluxes

F_face = q_face · C_upwind · L − D · h_face · L · ΔC/Δx

First-order upwind advection (unconditionally TVD)
Fickian diffusion: central-difference gradient
Wet-dry edges: flux only if both cells wet
Read-only access to old C field — fully parallel

Pass 2: Mass Update & Monotonicity

M_new = max(0, M_old + net_flux · Δt / A)

C_new = min(M_new / h_new, max(C_neighbors))

M_old = C_old · h_old · A (stored mass)
Positivity enforced: M ≥ 0
Monotonicity clamp: C cannot exceed neighbor maximum
Eliminates artificial extrema at wetting fronts

Why two-pass upwind ADE?

Two-pass design eliminates read-write race conditions, enabling fully parallel GPU execution. First-order upwind is unconditionally TVD, preventing unphysical concentration overshoots. Mass-conservative by construction on both mesh types.

Soil Infiltration: Horton Decay Model

Horton Infiltration Rate

f(t) = K_s + (f₀ − K_s) · exp(−k · t_wet)

Horton (1933); Rawls, Brakensiek & Miller (1983)

Per-Timestep Algorithm

Read cell soil type → look up f₀, K_s, k from USDA table
If h > 0: increment cumulative wet time t_wet += Δt
Compute rate: f(t) = K_s + (f₀ − K_s) exp(−k t_wet)
Potential infiltration: I_pot = f(t) · Δt
Limit by available water: I = min(I_pot, h)
Update depth: h_new = h − I
Store diagnostic rate: soil_infil_rate = I / Δt

USDA Soil Texture Classes

Class	K_s (mm/h)	f₀ (mm/h)
Sand	210.0	630.0
Loamy Sand	61.0	183.0
Sandy Loam	26.0	78.0
Loam	13.0	39.0
Silt Loam	6.8	20.4
Sandy Clay Loam	4.3	12.9
Clay Loam	2.3	6.9
Silty Clay	0.9	2.7
Clay	0.6	1.8

Why Horton decay?

Computationally cheap (one exponential per cell per timestep), widely validated for overland flow, requires only standard USDA soil parameters. Per-cell independent — trivially parallelizable. Backward compatible: no config section = module disabled.

Numerical Solver Summary

Both mesh types share the same physics — the triangular mesh adds 2nd-order spatial accuracy via MUSCL reconstruction.

Feature	Regular Cartesian Mesh	Unstructured Triangular Mesh
Flux solver	Roe approximate Riemann	Rotated Roe / HLL Riemann
Spatial order	1st order (piecewise constant)	2nd order (MUSCL + Barth-Jespersen)
Reconstruction	None (cell-averaged values)	Least-squares gradient + limiter
Wetting/drying	Hydrostatic reconstruction + Ritter	Hydrostatic reconstruction + Ritter
Bed friction	Semi-implicit Manning	Semi-implicit Manning
Time stepping	CFL · Δx / (\|u\| + c)	CFL · r_in / (\|u\| + c)
Mass limiter	Scale outflow per cell	Each edge ≤ 1/3 cell volume
Turbulence	Eddy viscosity (ν_t = c_vu_*h)	Eddy viscosity (ν_t = c_vu_*h)
ADE transport	Two-pass (upwind + diffusion)	Two-pass (upwind + diffusion)
Soil infiltration	Horton decay (per cell)	Horton decay (per cell)

                    Same conservation laws, same physical models — choose regular mesh for speed, triangular mesh for complex geometries and higher accuracy.
                

Key References: Saint-Venant (1871) • Godunov (1959) • Roe (1981) • van Leer (1979) • Barth & Jespersen (1989) • Audusse et al. (2004) • Horton (1933) • Rawls et al. (1983) • Toro (2001) • LeVeque (2002)

02

Parallelization & Performance

Parallelization Options

Technology	Memory Model	Use Case	CMake Flag
OpenMP	Shared memory	Multi-core workstations (default)	Always enabled
MPI	Distributed memory	HPC clusters, 100+ cores	`-DUSE_MPI=ON`
CUDA	GPU device memory	NVIDIA GPU acceleration	`-DUSE_CUDA=ON`
Hybrid	All combined	GPU clusters with multi-node	All flags combined

Regular Mesh Parallelization

OpenMP: loop collapse over rows/cols
MPI: 2D Cartesian domain decomposition with ghost cell exchange
CUDA: 2D thread grid (1 thread/cell)

Triangular Mesh Parallelization

OpenMP: edge & cell parallel loops
MPI: graph-based partitioning (METIS)
CUDA: 7 specialized kernels

CUDA GPU Kernels

7 specialized kernels for unstructured triangular mesh (edge-based finite volume):

Kernel Pipeline

1. Wet/Dry Classification

▼

2. Gradient Computation

▼

3. Barth-Jespersen Limiting

▼

4. Edge-Based Roe Flux

▼

5. Flux Accumulation (atomicAdd)

▼

6. State Update

▼

7. CFL Reduction

⚡ Performance

Speedup	10–50x over single CPU core
Min GPU	Compute Capability 6.0+
Architecture	Pascal, Volta, Turing, Ampere
CUDA version	11.0+

                            Key design: Race-free flux accumulation via atomicAdd for shared edges between triangular cells.
                        

03

Applications & Tools

Application Areas

🌊 Flood Risk Assessment

Flood extent & depth mapping
Event timing prediction
Back-flooding dynamics
Infrastructure vulnerability

❄️ Snowmelt Flooding

Spring freshet events
Prairie fill-and-spill
Frozen soil effects
Nutrient transport (WINTRA)

🧪 Solute Transport

Contaminant plume tracking
Point source dispersion
Tracer experiments
Water quality coupling

🌱 Soil & Infiltration

Infiltration excess runoff
Soil moisture dynamics
Land-use change impacts
Variable contributing areas

🌎 Wetland Dynamics

Depressional storage
Connectivity analysis
Flow pathway mapping
Fill-spill-merge cascades

🌡️ Climate Impact Studies

Extreme precipitation events
Scenario modelling
Land-use change
Hydrological response

Technical Stack

Languages & Libraries

C++17	Core simulation engine
CUDA C++	GPU kernels
Python	Pre/post-processing tools
Armadillo	Linear algebra (regular mesh)
nlohmann/json	Configuration parsing
WebGL/JS	3D visualization

Build & Deployment

CMake 3.10+	Build system
OpenMP	Shared-memory parallelism
OpenMPI	Distributed parallelism
CUDA 11.0+	GPU acceleration
Gmsh	Mesh generation (trimesh)

                            All compile options controlled via CMake flags:

                            -DUSE_TRIMESH=ON -DUSE_CUDA=ON -DUSE_MPI=ON

WebGL 3D Viewer

Interactive browser-based 3D visualization — no plugins or backend required:

🎨 Terrain Rendering

3D heightmap from DEM
Hillshade lighting
Satellite imagery overlay
Adjustable height exaggeration

🌊 Flow Visualization

GPU particle tracking (65k+)
Water depth colormaps
Concentration overlays
Velocity-driven animation

🎮 Interactive Controls

Time slider & playback speed
Orbit camera (mouse/touch)
Toggle layers on/off
Fullscreen mode

                    Fully self-contained static files — host on any web server, GitHub Pages, or open directly in a browser. No backend needed.
                

Summary

⚡ Fast

GPU-accelerated with CUDA. OpenMP + MPI for HPC. 10–50x speedup on modern GPUs.

🔧 Flexible

Regular or triangular meshes. JSON-driven configuration with optional ADE transport and soil infiltration modules.

🌐 Visual

WebGL 3D viewer with particle tracking, satellite imagery, and real-time animation.

                            FLUXOS bridges physics-based hydrodynamics with modern GPU computing

                            — enabling fast, high-resolution flood and solute transport simulations

Thank You

Questions & Resources

GitHub: github.com/ue-hydro/FLUXOS_cpp

Documentation: fluxos-cpp.readthedocs.io

Contact: diogo.costa@uevora.pt

FLUXOS

Contents

What is FLUXOS?

Core Capabilities

2D Hydrodynamics

Solute Transport

Soil Infiltration

GPU Acceleration

WebGL 3D Viewer

Multi-Format Output

Coupled Host Models

Mesh & Numerical Methods

Regular Cartesian Mesh

📈 Structured Grid

🔧 Solver Details

🕑 Performance

Unstructured Triangular Mesh

🔶 Edge-Based Finite Volume

🕑 Performance

🧮 Advanced Numerics

SWE: Conservative Vector Form

State & Flux Vectors

Source Terms

Finite Volume Discretization

Cell-Centered Godunov Scheme

CFL Adaptive Time-Stepping

Roe Approximate Riemann Solver

Roe Flux Structure

Hydrostatic Reconstruction

Wetting/Drying & Bed Friction

Ritter Dry-Front Solution

Velocity Desingularization

Semi-Implicit Manning Friction

MUSCL Reconstruction & Limiter

MUSCL Reconstruction

Edge-Normal Rotation

Barth-Jespersen Limiter

ADE Solver: Two-Pass Formulation

Pass 1: Face Fluxes

Pass 2: Mass Update & Monotonicity

Soil Infiltration: Horton Decay Model

Per-Timestep Algorithm

USDA Soil Texture Classes

Numerical Solver Summary

Parallelization & Performance

Parallelization Options

Regular Mesh Parallelization

Triangular Mesh Parallelization

CUDA GPU Kernels

Kernel Pipeline

⚡ Performance

Applications & Tools

Application Areas

🌊 Flood Risk Assessment

❄️ Snowmelt Flooding

🧪 Solute Transport

🌱 Soil & Infiltration

🌎 Wetland Dynamics

🌡️ Climate Impact Studies

Technical Stack

Languages & Libraries

Build & Deployment

WebGL 3D Viewer

🎨 Terrain Rendering

🌊 Flow Visualization

🎮 Interactive Controls

Summary

⚡ Fast

🔧 Flexible

🌐 Visual

Thank You