Wind NS Simulation — Technical Documentation

Full 3D incompressible Navier-Stokes solver for wind flow over terrain with surface roughness

Overview

The Wind NS simulation solves the full 3D incompressible Navier-Stokes equations to model how a horizontal wind flows over real terrain. Unlike the 2D Shallow Water approach, this solver resolves the complete three-dimensional velocity field — including how wind accelerates over ridges, how it decelerates and recirculates in the lee, and how the logarithmic boundary layer develops above the surface. Surface roughness varies per land cover type, derived from the ESA WorldCover dataset, so forests slow the wind more than bare rock or water.

The simulation uses a Chorin projection method (fractional step) with semi-Lagrangian advection (stable for any CFL number) and a Jacobi pressure Poisson solver. It runs in the browser via a Web Worker on a 3D rectangular grid. Terrain cells are flagged as solid (no-flow); the wind profile is initialised with a logarithmic law anchored to each cell's land-cover roughness length.

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.

Fetch land cover — ESA WorldCover 2021 COG GeoTIFF is fetched for the bounding box and resampled to the grid. Each cell is assigned an aerodynamic roughness length z₀ from the WorldCover class. A slope-based heuristic is used as fallback when the file is unavailable.

Build 3D grid and solid mask — the NX × NY × NZ Cartesian grid is constructed above the terrain base. Each cell is flagged solid if its centre lies below the terrain surface; fluid otherwise. Per-column roughness lengths are stored for the log-profile boundary condition.

Initialise velocity with log profile — every fluid cell is given an initial velocity following the logarithmic profile anchored to the inflow speed at z_ref above the local terrain surface. Pressure is initialised to zero everywhere.

Time-march NS — the solver steps forward: adaptive time step (advective + diffusive CFL) → semi-Lagrangian advection → explicit diffusion → apply boundary conditions → compute divergence → Jacobi pressure Poisson iterations → velocity projection → re-apply boundary conditions. Repeated until the flow reaches a quasi-steady state.

Compute surface fields — at each saved frame: surface wind speed √(U²+V²), vertical wind W, kinematic pressure P, and horizontal vorticity ∂V/∂x − ∂U/∂y are extracted from the lowest fluid cell above each terrain column.

Export to IndexedDB — snapshots of the 2D surface fields and the full 3D velocity and pressure arrays are stored in the browser's IndexedDB, from where the viewer reads them.

Physics Modelled

1. Incompressible Navier-Stokes equations

The flow is governed by the 3D incompressible Navier-Stokes equations. Density is assumed constant (no buoyancy, no temperature). The coordinate system is x = East, y = South, z = Up.

Momentum: ∂U/∂t + (u·∇)U = −∂p/∂x + ν·∇²U ∂V/∂t + (u·∇)V = −∂p/∂y + ν·∇²V ∂W/∂t + (u·∇)W = −∂p/∂z + ν·∇²W Continuity (incompressibility): ∂U/∂x + ∂V/∂y + ∂W/∂z = 0 where u = (U, V, W) is the 3D velocity p is kinematic pressure (Pa / ρ = m²/s²) ν is kinematic (turbulent) viscosity (m²/s)

There is no gravity term in the momentum equations (the flow is horizontal and density is constant, so buoyancy is absent). The only vertical forcing is through pressure gradients created by the terrain geometry.

2. Logarithmic wind profile

Near a rough surface, the mean wind speed follows a logarithmic profile with height above ground. This is the standard atmospheric surface-layer law:

U(z) = U_ref · ln(z / z₀) / ln(z_ref / z₀) for z > z₀ U(z) = 0 for z ≤ z₀ where z = height above local terrain surface z₀ = aerodynamic roughness length (land-cover dependent) z_ref = reference height where speed = U_ref (typically 100 m) U_ref = reference wind speed (user parameter)

This profile is applied at simulation initialisation and at inflow boundary cells every time step. It ensures the wind is slower close to rough surfaces (forest, built-up areas) and faster over smooth surfaces (water, snow).

3. Surface roughness from land cover

The roughness length z₀ is assigned per grid column from the ESA WorldCover 10 m land cover map. When ESA data is unavailable, a terrain-slope heuristic provides a fallback estimate.

Land cover class	ESA code	z₀ (m)	Physical meaning
Permanent water	80	0.0002	Charnock relation (wave-dependent); almost frictionless
Snow / ice	70	0.001	Smooth, hard surface
Bare / sparse veg.	60	0.01	Desert, rocky alpine
Grassland	30	0.03	Short grass — reference surface
Herb. wetland	90	0.05	Reed beds and marshes
Moss / lichen	100	0.05	Alpine heath
Cropland	40	0.10	Cereal crops, row crops
Shrubland	20	0.30	Garrigue, maquis
Mangrove	95	0.50	Dense coastal vegetation
Tree cover	10	1.00	Closed-canopy forest
Built-up	50	1.50	Urban: buildings impose large drag

Fallback heuristic when ESA WorldCover is unavailable:

if z > 2400 m → Snow/ice (z₀ = 0.001 m) else if slope > 0.5 → Bare/sparse (z₀ = 0.01 m) else if slope > 0.25 → Shrubland (z₀ = 0.30 m) else if z < 500 m and slope < 0.05 → Cropland (z₀=0.10 m) else → Grassland (z₀=0.03 m)

4. Terrain blocking — solid mask

The 3D grid is divided into fluid and solid cells. A cell at height z_c above the lowest terrain point is solid if its centre lies below the local terrain surface:

z_c = (k + 0.5) · dz [cell centre height above base] SOLID[k,j,i] = 1 if z_c ≤ Z_rel[j,i] (inside terrain) SOLID[k,j,i] = 0 if z_c > Z_rel[j,i] (air)

Solid cells are enforced with a no-slip condition: U = V = W = 0. The terrain therefore acts as a physical obstacle that diverts and accelerates the flow.

Numerical Method (Chorin Projection)

Each time step applies the following sequence to advance the velocity and pressure fields:

Semi-Lagrangian advection — for each fluid cell, back-trace the Lagrangian departure point by following the velocity backward in time, then interpolate the field value there using trilinear interpolation:
x_d = x − U·dt, y_d = y + V·dt, z_d = z − W·dt
(the +V sign compensates for j increasing southward). This scheme is unconditionally stable for any Courant number, allowing time steps larger than explicit schemes.

Explicit diffusion — apply the viscous term with forward Euler:
u* = u_adv + ν · ∇²u · dt
The Laplacian uses a standard 7-point 3D stencil (central differences). The viscosity ν is clamped by a diffusion-CFL condition to prevent instability:
ν ≤ diff_safety · min(dx², dy², dz²) / dt

Boundary conditions on u* — apply inflow (log profile), no-slip on solid cells, free-slip on the top, and Neumann on outflow edges before the pressure solve.

Divergence of u* — compute:
div = ∂U*/∂x − ∂V*/∂y + ∂W*/∂z
(the −∂V*/∂y accounts for the y-axis convention where j increases southward). Non-zero divergence indicates mass imbalance that the pressure step must correct.

Pressure Poisson equation — solve ∇²p = div/dt using Jacobi iterations:

p_new[i,j,k] = ( (p_e + p_w)/dx² + (p_n + p_s)/dy² + (p_u + p_d)/dz² − div/dt ) / (2/dx² + 2/dy² + 2/dz²)

Outflow boundaries: p = 0 (Dirichlet). Solid and top: ∂p/∂n = 0 (Neumann, i.e. copy from fluid neighbour).

Velocity projection — subtract the pressure gradient to enforce incompressibility:
U = U* − dt · ∂p/∂x
V = V* + dt · ∂p/∂y (+sign: y-axis convention)
W = W* − dt · ∂p/∂z
After projection, ∇·u ≈ 0 (to the residual accuracy of the Jacobi solver).

Re-apply boundary conditions — enforce no-slip and inflow conditions after projection to prevent the pressure step from violating them.

Time step calculation

dt_adv = CFL · min(dx, dy) / |u|_max [semi-Lagrangian advective limit] dt_diff = diff_safety · min(dx², dy², dz²) / ν [diffusion stability limit] dt = min(dt_adv, dt_diff, dt_max) [smallest wins]

Because semi-Lagrangian advection is unconditionally stable, the time step is often limited by the diffusion condition or by dt_max. The CFL factor can safely exceed 1.0 (typical: 1.5–2.5) without instability from advection.

Boundary Conditions

Boundary	Velocity	Pressure
Solid terrain (any k)	No-slip: U = V = W = 0	Neumann: ∂p/∂n = 0
Inflow lateral edges	Log profile: U = U_ref·f(z,z₀), V = V_ref·f(z,z₀), W = 0	Neumann
Outflow lateral edges	Neumann: copy from interior neighbour	Dirichlet: p = 0
Top (k = NZ−1)	Free-slip: U_top = U_{NZ-2}, V_top = V_{NZ-2}, W_top = 0	Neumann
Bottom (k = 0, not solid)	Neumann: U_0 = U_1, V_0 = V_1, W_0 = 0	Neumann

The inflow log factor is computed as:

f(z, z₀) = min(1.0, ln(z / z₀) / ln(z_ref / z₀) ) for z > z₀ f(z, z₀) = 0 for z ≤ z₀

where z is the height above the local terrain surface (not above sea level). The factor is clamped to 1 so the log profile does not overshoot the reference speed above z_ref.

Parameters and Their Effects

Grid and domain

Parameter	Typical value	Effect
NX, NY	32–80	Horizontal grid resolution. Finer grids resolve narrower terrain features and tighter boundary layers but increase memory (NX·NY·NZ cells) and computation time.
NZ	20–40	Vertical layers. More layers improve the representation of the log profile near the surface and resolve vertical flow structure (separation bubbles, lee waves). Decreasing dz tightens the diffusion time-step limit.
h_domain	1000–3000 m	Domain height above the highest terrain point. Should be at least 3–5× the maximum ridge height to avoid the top boundary interfering with flow over ridges. Too small truncates the wind profile aloft.
dt_max	5–30 s	Hard cap on the time step. Prevents the adaptive time step from becoming too large in regions with very slow flow. Reduce if the solution diverges.

Wind input

Parameter	Typical value	Effect
Wind direction	0–360°	Compass bearing the wind comes from. Determines the inflow edge(s) and the direction of the log-profile inflow. Controls which slopes are windward (acceleration, uplift) and which are lee (deceleration, recirculation).
U_ref (wind speed)	2–25 m/s	Reference wind speed at z_ref above the surface. Scales the entire inflow log profile and therefore all velocities in the domain. Stronger wind amplifies all flow features proportionally (the NS equations are linear in U for steady laminar flow — but non-linear advection creates speed-dependent patterns).
z_ref	50–200 m	Height at which U_ref is specified. The log profile is anchored here. Increasing z_ref while keeping U_ref constant effectively shifts the profile upward, producing faster flow near the surface for the same reference speed.

Roughness and boundary layer

Parameter	Typical value	Effect
z₀ (log_z0)	0.001–1.5 m	Aerodynamic roughness length. Used as a scalar fallback when ESA WorldCover is not loaded. Larger z₀ → slower near-surface wind, deeper boundary layer, more drag on the flow. Small z₀ (water, snow) → fast surface wind, thin boundary layer. Overridden per-cell by the land cover map when available.
ν (visc)	1–500 m²/s	Kinematic (turbulent eddy) viscosity. Diffuses momentum laterally and vertically. Higher ν → smoother velocity field, lee recirculations are broadened and weakened, sharp shear layers are smeared. Lower ν → sharper gradients, stronger separation zones, risk of numerical oscillations if too low for the grid resolution.
diff_safety	0.1–0.25	Safety factor for the diffusion time-step limit: `ν ≤ diff_safety · Δ² / dt`. Lower values make the scheme more conservative (smaller effective ν per step). Increase only if the diffusion step is artificially limiting dt on fine grids.

Pressure solver

Parameter	Typical value	Effect
p_iter	20–80	Number of Jacobi iterations per time step. More iterations → pressure field converges closer to the true solution of ∇²p = div/dt → velocity field is more nearly divergence-free. Too few leaves residual divergence that accumulates over time and can destabilise the solution. Too many wastes computation; for most grids 30–60 is sufficient.
CFL	1.0–2.5	Courant number for the semi-Lagrangian advection step. Since semi-Lagrangian is unconditionally stable, values above 1 are safe. Higher CFL → larger time steps → faster wall-clock time per simulated second. Very large CFL (>3) may degrade accuracy of the back-traced departure point (trilinear interpolation error grows).

Output Fields

Field	Symbol	Description
Surface wind speed	√(U²+V²) at surface	Horizontal wind speed at the lowest fluid cell above the terrain. Blue→green→yellow. Acceleration over ridges, deceleration in wakes. Primary indicator for mechanical conditions at flying height.
Vertical wind	W at surface level	Vertical velocity at the first fluid cell. Red = rising (orographic uplift, lee convergence). Blue = sinking (lee side, valleys, divergence). Most directly relevant to ridge soaring and rotor detection.
Pressure	p at surface level	Kinematic pressure anomaly (Pa/ρ). High pressure (red) upwind of obstacles — stagnation zones. Low pressure (blue) in wakes and over crests (Bernoulli acceleration). Useful for identifying flow separation points.
Vorticity	∂V/∂x − ∂U/∂y at surface	Horizontal vorticity at the surface level. Red = cyclonic (CCW); blue = anticyclonic (CW). Marks shear layers, separated wakes, and the edges of recirculation zones. High vorticity → turbulence risk.

Simplifications

Incompressibility: density is constant everywhere — no sound waves, no compressibility effects. Valid for wind speeds well below ~100 m/s. Eliminates the need for an energy equation.
Constant turbulent viscosity ν: a single scalar viscosity is used rather than a spatially varying turbulence model (e.g. k-ε or Smagorinsky). The Smagorinsky model is used in the Thermal NS simulation but not here — ν must be set manually.
Neutral atmosphere: no stratification, no buoyancy. The log profile assumes a neutral (adiabatic) atmospheric surface layer. Stable or unstable stratification would modify the profile (Monin-Obukhov similarity theory) but is not implemented.
Flat top boundary: the domain top is a horizontal plane with free-slip. This reflects vertical momentum back into the domain rather than allowing it to radiate upward — an artefact for strong updrafts reaching the top.
Jacobi pressure solver: converges slowly for large grids; many iterations may be needed to reach acceptable divergence residuals. SOR or a multigrid solver would converge faster but is not implemented.
Steady inflow: wind speed and direction are fixed in time. No gusts, no diurnal variation, no mesoscale dynamics.
No terrain-following coordinates: the grid is Cartesian; terrain is represented as blocked cells. Cells that straddle the exact terrain surface are either fully solid or fully fluid — no partial-cell treatment. This slightly over-estimates the terrain height by up to one dz.

Limitations

No thermal effects — no solar heating, no buoyancy, no convective thermals. For thermally driven vertical motion, use the Thermal NS simulation.
No atmospheric stratification — gravity waves, lee waves, and wave breaking (which cause severe turbulence at altitude) require a stratified atmosphere and are absent here.
No Coriolis force — negligible for domains < 50 km but formally absent.
No moisture — clouds, precipitation, and latent heat are not represented.
No synoptic forcing — large-scale pressure gradients, fronts, and mesoscale weather systems are not modelled. The inflow is a simple uniform log profile, not a realistic atmospheric wind profile with jets or wind shear.
Slow pressure convergence — the Jacobi solver requires many more iterations than SOR or multigrid to achieve low divergence residuals, especially on fine 3D grids. The flow may not be fully divergence-free within practical iteration counts.
No wall model — the log profile is enforced at inflow boundaries and at initialisation, but there is no dynamic wall model that adjusts the near-wall velocity at internal solid faces each step. Flow separation from steep terrain faces may be under-represented.

Data Sources

Dataset	Used for	Source	Resolution
Elevation	Terrain mesh, solid mask, height above ground for log profile	AWS Terrain Tiles — Terrarium RGB encoding (IGN RGE Alti for France)	~7 m (zoom 14)
Land cover	Aerodynamic roughness length z₀ per cell	ESA WorldCover 2021 (Cloud-Optimised GeoTIFF)	10 m
Simulation output	3D velocity/pressure fields and 2D surface maps 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. ESA WorldCover is fetched as a Cloud-Optimised GeoTIFF (COG) — only the portion covering the bounding box is downloaded via HTTP range requests.

Tips for Useful Results

Allow at least 100–300 time steps for the flow to develop from the initial log profile to a quasi-steady state. The simulation is not initialised with zero velocity, so convergence is faster than a cold start.
Increase p_iter to 50–80 if the vertical wind field looks noisy or shows checkerboard patterns — these are symptoms of insufficient pressure convergence (residual divergence).
Set ν to 5–50 m²/s for mountain-scale flows (dx ≈ 200 m). At coarser resolution, increase ν to damp sub-grid instabilities. At finer resolution, ν can be reduced for sharper features.
The vertical wind (W) map is most relevant for paragliding: red patches on windward ridges indicate dynamic (orographic) lift zones; blue patches indicate sink on the lee side.
The pressure map identifies stagnation zones (high p, blue upwind) and suction peaks (low p over crests). The crest is where orographic lift transitions to lee-side sink.
Load ESA WorldCover land cover for per-cell roughness. Without it, a single scalar z₀ is used for the whole domain, which misses the contrast between forest (slow wind aloft, strong drag) and bare rock (fast near-surface wind).
Use a higher CFL (1.5–2.0) to speed up the simulation without sacrificing stability. Only reduce CFL if the back-traced departure points land outside the domain frequently (very strong winds with coarse grids).