Wind Flow Simulation — Technical Documentation

2D Shallow Water Equations solver for orographic wind analysis over real terrain

Overview

The Wind Flow simulation solves the 2D Shallow Water Equations (SWE) to model how a horizontal wind interacts with terrain. The atmosphere is approximated as a single layer of fluid with a free surface, flowing over an orographic (terrain) floor. Where the terrain rises, the fluid layer thins; where the terrain is lower, the layer is deeper. This reduced-gravity hydraulic analogy captures terrain-blocking, orographic lifting, lee-side acceleration, and hydraulic jumps — the main mechanisms that shape wind patterns relevant to paragliding site analysis.

The simulation runs entirely in the browser using a Web Worker on a 2D grid aligned with the terrain bounding box. It is fast (seconds to a few minutes) and does not require a 3D pressure solve. Elevation data is fetched from Terrarium tiles; no land cover is used in this model.

Pipeline Overview

Fetch terrain — AWS Terrarium RGB-encoded tiles are downloaded for the bounding box, decoded to metres, and bilinearly resampled to the NX × NY simulation grid.

Build 2D grid with elevation — the terrain height Z(i,j) is stored per cell. The undisturbed layer depth at each cell is initialised to h = H₀ − Z. Cells where H₀ − Z ≤ h_min are immediately treated as dry (terrain-blocked).

Initialise flow field — the velocity is set uniformly to the inflow values (U_in, V_in) everywhere. The free-surface height η = h + Z is computed from the initial h.

Time-march SWE — the solver steps forward in time: adaptive CFL time step → continuity (conservative flux divergence) → upwind momentum advection → free-surface pressure gradient → bottom drag → turbulent diffusion → wet/dry masking → boundary conditions. Repeated until the flow reaches a quasi-steady state.

Compute derived fields — at each saved frame: wind speed |u|, Froude number Fr = |u| / √(g'·h), vertical wind proxy w_vert = U·∂Z/∂x + V·∂Z/∂y − h·div(u,v), and vorticity ζ = ∂V/∂x − ∂U/∂y are computed from the current (h, U, V) state.

Export to IndexedDB — snapshots of all fields are stored in the browser's IndexedDB as the simulation runs, enabling the viewer to display them without waiting for the full run.

Physics Modelled

1. Shallow Water Equations (SWE)

The atmosphere is treated as a single, vertically integrated layer of depth h flowing over terrain of elevation Z. The free-surface height is η = h + Z. The governing equations in conservative form are:

Continuity: ∂h/∂t = −∂(h·U)/∂x − ∂(h·V)/∂y Momentum (U): ∂U/∂t = −U·∂U/∂x − V·∂U/∂y − g'·∂η/∂x − C_d·|u|·U/h + ν·∇²U Momentum (V): ∂V/∂t = −U·∂V/∂x − V·∂V/∂y − g'·∂η/∂y − C_d·|u|·V/h + ν·∇²V where η = h + Z (free-surface elevation) |u| = √(U² + V²) (wind speed) ∇² = ∂²/∂x² + ∂²/∂y² (2D Laplacian)

The four physical terms in each momentum equation are:

Advection — U·∂U/∂x + V·∂U/∂y: the fluid carries its own momentum. Causes flow to accelerate through constrictions (mountain passes).
Free-surface pressure gradient — g'·∂η/∂x: drives flow from high to low free-surface. The terrain slope ∂Z/∂x is included via η, so flow is deflected around or over ridges.
Bottom drag — C_d·|u|·U/h: quadratic friction from the terrain surface, opposing the flow. Stronger where the layer is thin or flow is fast.
Turbulent viscosity — ν·∇²U: diffuses momentum laterally, smoothing velocity gradients and preventing sharp shocks from going numerically unstable.

2. Reduced gravity

The parameter g' (reduced gravity) replaces the true gravitational acceleration g = 9.81 m/s². In the full atmosphere, a density contrast Δρ/ρ between the flowing layer and the air above it gives g' = g · Δρ/ρ. In practice, g' controls the speed of gravity waves in the layer:

c = √(g' · h) [gravity wave phase speed]

A small g' (e.g. 0.2 m/s²) mimics a weakly stratified interface — waves travel slowly and the flow responds gradually to terrain. A large g' (e.g. 2 m/s²) stiffens the free surface, causing the layer to behave more like a rigid-lid flow with faster pressure adjustment.

3. Orographic lifting and vertical wind

Although the SWE is 2D, a proxy for the vertical wind velocity is derived from terrain-forced lifting and flow convergence:

w_vert = U·∂Z/∂x + V·∂Z/∂y − h·(∂U/∂x + ∂V/∂y) ↑ orographic lifting ↑ convergence-driven ascent

The first term is the kinematic uplift: wind hitting a rising slope must go upward. The second term is the convergence: where horizontal flow converges (∂U/∂x + ∂V/∂y < 0), mass is forced upward. Positive w_vert (red in the viewer) indicates rising air — the primary indicator for orographic thermals and ridge lift.

4. Froude number

The Froude number distinguishes two flow regimes analogous to sub-sonic / super-sonic in gas dynamics:

Fr = |u| / √(g' · h)

Fr < 1 (sub-critical): disturbances can propagate upstream. Flow smoothly adjusts to terrain; typical upwind of a ridge.
Fr > 1 (super-critical): disturbances cannot travel upstream. Flow accelerates through passes and over crests, producing hydraulic jumps on the lee side — violent, turbulent transitions where super-critical flow suddenly decelerates. In the real atmosphere this corresponds to mountain wave / rotor conditions.

5. Vorticity

The 2D vertical vorticity measures local rotation in the horizontal wind field:

ζ = ∂V/∂x − ∂U/∂y

Positive ζ (red) is cyclonic (counter-clockwise); negative ζ (blue) is anticyclonic. Strong vorticity patches appear at the edges of mountain wakes, behind ridges, and in the shear layers surrounding accelerated jets. High vorticity is associated with turbulence.

Numerical Method

The SWE are integrated with an explicit finite-difference scheme (forward Euler in time, second-order central differences in space), stabilised by a CFL time-step condition.

Time step (CFL condition) — chosen adaptively each step to keep the scheme stable:
dt = CFL · min(dx, dy) / max(c_est, 0.1)
where c_est = |u|_max + √(g'·H₀) is the fastest signal speed in the domain (advection + gravity wave). CFL is typically 0.3–0.5.

Continuity — conservative flux divergence — the depth update uses the conservative form to guarantee exact mass conservation:
h_new = h − dt · [∂(hU)/∂x + ∂(hV)/∂y]
Fluxes are computed with central differences on the interior.

Momentum — upwind advection — the advection terms use first-order upwind differences (one-sided based on the local flow direction) to prevent non-physical oscillations:
If U > 0: ∂U/∂x ≈ (U[i] − U[i−1]) / dx (backward difference)
If U < 0: ∂U/∂x ≈ (U[i+1] − U[i]) / dx (forward difference)

Pressure gradient — the free-surface gradient uses central differences:
−g' · (η[i+1,j] − η[i−1,j]) / (2·dx)
where η = h + Z is evaluated at the current time level.

Drag and diffusion — applied to the intermediate velocity after advection and pressure. Drag is quadratic in speed; diffusion uses the standard 5-point Laplacian stencil.

Wet/dry masking — cells where h < h_min are treated as dry (terrain-blocked). Velocities are zeroed and fluxes are suppressed to avoid division by near-zero depth and to represent complete terrain blocking.

Boundary conditions — inflow edges (where wind enters the domain) are held at h = H₀ − Z, U = U_in, V = V_in. Outflow edges use Neumann conditions (copy the interior value), allowing fluid to leave freely.

Boundary Conditions

The wind direction determines which edges are inflow and which are outflow. For a wind with eastward component U_in > 0, the left (west) edge is inflow and the right (east) edge is outflow. Both conditions are applied simultaneously for oblique winds.

Edge	Condition	Applied when
Left (West)	Inflow: h = H₀−Z, U = U_in, V = V_in	U_in > 0
Right (East)	Inflow: h = H₀−Z, U = U_in, V = V_in	U_in < 0
Top (North)	Inflow: h = H₀−Z, U = U_in, V = V_in	V_in < 0
Bottom (South)	Inflow: h = H₀−Z, U = U_in, V = V_in	V_in > 0
Any outflow edge	Neumann: copy from interior neighbour	Opposite sign of above

Parameters and Their Effects

Wind input

Parameter	Typical range	Effect
Wind direction	0–360°	Compass bearing the wind comes from. Determines which domain edges are inflow. Changes which slopes are windward (orographic lift) and which are lee (sink/rotors).
Wind speed (U_in, V_in)	0–30 m/s	Inflow velocity magnitude. Stronger wind → faster flow, larger Froude number, stronger orographic lifting and more pronounced hydraulic jumps. Too strong with low ν may cause instability.

Layer physics

Parameter	Typical range	Effect
H₀ (layer depth)	500–3000 m	Undisturbed atmospheric layer depth. Mountains taller than H₀ completely block the flow (h → 0, terrain-blocked). A shallow H₀ makes more of the terrain protruding and emphasises blocking effects. A deep H₀ allows flow to pass over most terrain.
g' (reduced gravity)	0.1–2.0 m/s²	Controls the gravity wave speed `c = √(g'·h)` and thus the Froude number regime. Small g' → slow wave propagation, flow is more easily blocked and deflected around obstacles. Large g' → stiff free surface, flow adjusts rapidly, fewer blocking effects. Also sets the transition wind speed for hydraulic jumps.
C_d (bottom drag)	0–0.05	Quadratic friction coefficient. Decelerates the flow proportional to speed². Higher drag → flow slows faster downstream of ridges, wake zones are wider, jet acceleration through passes is damped. Zero drag allows unphysical indefinite acceleration.
ν (viscosity)	0–2000 m²/s	Turbulent (eddy) viscosity. Diffuses lateral velocity gradients. Acts as a numerical stabiliser and physical proxy for sub-grid turbulence. Too low: sharp shocks may become unstable. Too high: flow features are over-smoothed and fine-scale patterns are lost.

Numerical settings

Parameter	Typical range	Effect
CFL	0.2–0.5	Courant number controlling the time-step size. Smaller → more stable but slower (more steps per second of simulated time). Values above 0.5 risk numerical instability with explicit upwind schemes.
h_min	0.1–5 m	Wet/dry threshold. Cells with h below this are treated as completely blocked terrain. Too large: shallow-but-valid flow is incorrectly suppressed. Too small: near-zero h causes division instability in the drag term.
NX, NY (grid)	32–256	Horizontal resolution. Finer grids resolve narrower terrain features (passes, gullies) and produce sharper gradients. Runtime scales as NX·NY per step.

Output Fields

Field	Symbol	Description
Wind speed	\|u\| = √(U²+V²)	Horizontal wind speed at the layer. Blue→cyan→yellow gradient. Acceleration through mountain passes and deceleration in wakes are clearly visible.
Vertical wind	w_vert	Proxy for vertical velocity (orographic uplift + convergence). Red = rising air (windward slopes, convergence zones). Blue = sinking air (lee slopes, divergence). Most directly relevant to paragliding ridge and thermal soaring.
Froude number	Fr = \|u\|/√(g'·h)	Blue (Fr < 0.5, sub-critical) → yellow (Fr ≈ 1, critical) → red (Fr > 1, super-critical). Red zones indicate locations prone to hydraulic jumps and strong turbulence on the lee side.
Flow depth	h	Layer thickness. Thin (dark) where terrain is high; deep (bright) over valleys. Zero/near-zero indicates terrain-blocked cells.
Vorticity	ζ = ∂V/∂x − ∂U/∂y	Red = cyclonic rotation; blue = anticyclonic. Marks shear layers at the edges of mountain wakes and along the flanks of jets.
Elevation	Z(x, y)	Terrain height for spatial reference, overlaid with a hillshade.

Simplifications

Single layer (depth-integrated): the vertical structure of the atmosphere is collapsed into one layer depth h. There is no wind shear with height, no vertical profile, and no interaction between atmospheric layers.
Hydrostatic pressure: pressure is purely hydrostatic (g'·η). Non-hydrostatic effects (important in steep terrain with fast flow) are absent.
Reduced gravity: g' is a tunable proxy for the density contrast between the flowing layer and the air above — it does not correspond to a specific atmospheric layer without careful calibration.
No vertical wind profile: the inflow is uniform with height; logarithmic boundary layer effects are not modelled.
Upwind advection (1st order): numerically diffusive; sharp fronts (shocks) are smeared over 2–3 cells.
Explicit time stepping: the CFL condition forces small time steps at high wind speeds or fine resolution, making the simulation slower than semi-Lagrangian schemes.
Uniform bottom drag: C_d is the same everywhere regardless of land cover or terrain roughness.

Limitations

No thermal effects — the simulation is purely dynamical. No solar heating, no buoyancy, no convective thermals. Temperature is absent from the model.
No 3D turbulence — turbulence is parameterised by a single constant lateral viscosity ν. The complex 3D turbulent structure of mountain wakes (rotors, eddies) is not resolved.
No atmospheric stratification — there is no stable stratification or inversion layer. The reduced gravity g' is a crude proxy but does not capture the full physics of a stratified atmosphere (gravity waves, wave breaking, critical layers).
No Coriolis effect — negligible for domains < 50 km but formally absent.
No moisture — clouds, precipitation, and latent heat release are not modelled.
Steady inflow — the inflow wind speed and direction are constant in time. Gusts, diurnal variation, and mesoscale variability are not represented.
No valley channelling detail — fine-scale channelling in narrow valleys may be under-resolved at typical grid spacings (200–500 m per cell).

Data Sources

Dataset	Used for	Source	Resolution
Elevation	Terrain height Z(i,j), layer depth initialisation, orographic lifting	AWS Terrain Tiles — Terrarium RGB encoding (IGN RGE Alti for France)	~7 m (zoom 14)
Simulation output	Snapshots of h, U, V and derived fields read by the viewer	Browser IndexedDB (written by the Web Worker, read by the viewer)	—

Elevation tiles use the Terrarium RGB encoding: elevation = R×256 + G + B/256 − 32768 metres, delivered in Web Mercator (EPSG:3857) and reprojected to WGS84 before use. No land cover data is used — surface roughness is represented by the uniform drag coefficient C_d.

Tips for Useful Results

Set H₀ to the height of the atmospheric boundary layer or the ridge height you want to block — typically 800–2000 m.
Start with g' ≈ 0.5 m/s² and adjust until the Froude number map shows sub-critical flow upwind of ridges (Fr < 1, blue), transitioning to super-critical on the crests (Fr ≈ 1, yellow). This is the most physically meaningful regime.
Increase ν if the simulation becomes visually noisy or velocity spikes appear. Decrease it for sharper, more realistic gradients on coarser grids.
The vertical wind map is the most directly actionable output for paragliding: red patches on windward slopes indicate orographic lift zones.
The Froude number map identifies hazardous lee-side zones (Fr > 1) where hydraulic jumps and mechanical turbulence are likely.
Allow at least 200–500 time steps before reading results — the flow needs time to reach a quasi-steady state after the uniform initial condition.