NS Thermal Simulation — Technical Documentation

Full 3D Navier-Stokes atmospheric simulation for paragliding thermal analysis

Overview

The NS (Navier-Stokes) simulation solves the incompressible Boussinesq equations for atmospheric convection over real terrain. It models how solar heating of the ground creates temperature anomalies in the boundary layer, which generate buoyant thermals that rise through a stably stratified atmosphere.

The simulation runs entirely in the browser using a Web Worker, operating on a 3D rectangular grid that spans the terrain horizontally and from the lowest terrain point up to a configurable domain height vertically. Terrain cells are flagged as solid (no-flow) and heat is injected at the first fluid cell above the surface.

Pipeline Overview

Fetch terrain — AWS Terrarium RGB-encoded tiles are downloaded at the chosen zoom level, 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, resampled to the grid, and converted to per-cell thermal factor tf and time constant τ. A slope-based heuristic is used as fallback when the file is unavailable.

ERA5 calibration — Open-Meteo archive API is queried for historical 2 m air temperature and surface temperature at the site centre. The median clear-sky midday surface–air difference is used to set max_heat, or the user can override it manually.

Build 3D grid and solid mask — the NX × NY × NZ Cartesian grid is constructed. Each cell is flagged solid if its centre lies below the terrain surface; fluid otherwise. The first fluid cell above each terrain column is identified for surface heat-flux injection.

Compute solar heating maps — for each scenario time step, solar geometry, terrain shadow ray-casting, land-cover factors, and atmospheric transmittance are combined to produce the per-cell surface heat flux Q_Ks(x,y) that drives the simulation at that scenario time.

Spin-up — the solver is advanced for the configured number of spin-up steps before the first snapshot, allowing convective motions to develop from the initially quiescent (zero-velocity) field.

Scenario loop — for each output time step the corresponding heat flux is applied, then the NS solver runs the configured number of steps: Smagorinsky viscosity → semi-Lagrangian advection → explicit diffusion → buoyancy → stratification → surface flux → Red-Black SOR pressure → projection. 2D surface fields and 3D cross-section slices are snapshotted.

Export to IndexedDB — all scenario snapshots plus static fields (elevation, land-cover map) are packed and stored in the browser's IndexedDB, from where the Three.js viewer reads them.

Physics Modelled

1. Incompressible Boussinesq Navier-Stokes equations

The flow is governed by the incompressible Navier-Stokes equations under the Boussinesq approximation: density variations are neglected everywhere except in the buoyancy term, where temperature anomalies drive vertical acceleration.

∂u/∂t + (u·∇)u = −∇p + ∇·(νₜ ∇u) [momentum]
∇·u = 0 [incompressibility]
∂θ'/∂t + (u·∇)θ' = κₜ ∇²θ' − w·(DALR − ELR) + Q_surface [temperature]

Here u = (u, v, w) is the 3D velocity field (East, South, Up), p is the kinematic pressure (pressure divided by reference density), θ' is the potential temperature anomaly relative to the background profile, νₜ is the turbulent viscosity, and κₜ = νₜ / Prₜ is the turbulent thermal diffusivity.

2. Buoyancy

Warm air rises. The buoyancy term couples the temperature anomaly to the vertical momentum equation:

w* += (g / T₀) · θ' · dt

where g = 9.81 m/s², T₀ is the reference temperature in Kelvin (estimated from ERA5 reanalysis or set to 293 K), and dt is the timestep. A cell with θ' = +2 K at T₀ = 293 K produces a buoyancy acceleration of about 0.067 m/s² — enough to drive 1–3 m/s vertical velocities over several minutes.

3. Background atmospheric stability (stratification)

The atmosphere has a background temperature profile that resists vertical motion. When a parcel rises, it cools at the Dry Adiabatic Lapse Rate (DALR = 9.76 K/km). If the environment cools more slowly (ELR < DALR), the parcel becomes colder than its surroundings and is pushed back down — this is atmospheric stability.

θ'* -= w · (DALR − ELR) · dt

The stability parameter strat = DALR − ELR is positive when the atmosphere is stable (typical daytime: strat ≈ 3 K/km). When ELR approaches DALR the atmosphere is neutral and thermals penetrate much higher. If ELR > DALR the atmosphere is absolutely unstable (rare in practice).

4. Smagorinsky Large Eddy Simulation (LES) turbulence closure

Turbulent motions smaller than a grid cell are parameterised using the Smagorinsky-Lilly model. The turbulent viscosity is computed from the local strain-rate tensor:

νₜ = (Cs · Δ)² · √(2 SᵢⱼSᵢⱼ) Sᵢⱼ = ½ (∂uᵢ/∂xⱼ + ∂uⱼ/∂xᵢ) [strain-rate tensor] Δ = (dx · dy · dz)^(1/3) [filter width = geometric mean grid spacing]

The turbulent Prandtl number Prₜ = 0.7 gives the thermal diffusivity κₜ = νₜ / 0.7. The viscosity is clamped for numerical stability: νₜ ≤ 0.4 · min(dx,dy,dz)² / dt.

5. Solar heating and surface heat flux

At each scenario time step, the ground temperature anomaly is computed from sun geometry, terrain slope and shadowing, land cover, and atmospheric transmittance.

dT(x,y) = maxHeat · (I_direct + I_diffuse) · tf(x,y)

I_direct = cos(θ_inc) · T_atm · (1 − shadow)
I_diffuse = diffuseFrac · (1 − T_atm) · SVF

T_atm = T_base^(AM^T_exp) [Kasten-Czeplak atmospheric transmittance]
AM = 1 / [sin(el) + 0.50572 · (el_deg + 6.07995)^(−1.6364)] [airmass]
SVF = (1 + cos(slope)) / 2 [sky view factor for diffuse radiation]

θ_inc is the angle between the solar ray and the terrain surface normal (computed from the terrain gradient). The ground temperature anomaly is then converted to a continuous heat flux injected into the simulation:

Q_Ks(x,y) = dT(x,y) / (τ(x,y) · 3600) [K/s]

where τ is the ground thermal time constant (in hours, land-cover dependent). This flux is added each timestep to the potential temperature of the bottom fluid cell above the terrain surface.

6. Land cover

Each grid cell has a thermal factor tf (heating efficiency) and time constant τ derived from the ESA WorldCover land cover map:

Land cover	tf (heating)	τ (h)	Effect
Built-up	1.30	0.50	Strong, fast heating (heat island)
Bare/sparse	1.80	0.35	Strongest heating, fastest response — rocky terrain is the best thermal trigger
Cropland	0.85	0.85	Moderate heating
Grassland	1.00	1.00	Reference land cover
Shrubland	0.70	1.30	Moderate, slower
Moss/lichen	0.60	1.00	Moderate
Herb. wetland	0.35	1.50	Weak heating
Mangrove	0.50	2.00	Weak, slow
Tree cover	0.45	2.50	Forest: weak and slow — poor thermal source
Snow/ice	0.08	1.50	Almost no heating (high albedo)
Perm. water	0.03	15.0	Negligible heating, very high thermal mass

7. Terrain shadowing

Before the simulation, a ray-marching shadow pass determines which terrain cells are in shadow for each scenario time step. Shadowed cells receive only diffuse radiation, suppressing their surface heat flux. This creates realistic cold patches on north-facing slopes or in valleys in the morning and evening.

8. ERA5 calibration

The maximum heating parameter maxHeat can be calibrated automatically from ERA5 reanalysis data. The model fetches historical hourly 2m-temperature and surface temperature from the Open-Meteo API for the target site, computes the typical midday surface–air temperature difference, and fits maxHeat to match observations. This grounds the simulation in real meteorological data for the chosen date and location.

Numerical Method (Chorin Projection)

Each timestep applies the following sequence (Chorin fractional-step / projection method):

Smagorinsky turbulent viscosity — compute νₜ at every fluid cell from the current strain-rate tensor (9 velocity gradients via central differences).

Semi-Lagrangian advection — for each cell, backtrack the departure point x_d = x − u·dt and interpolate the field value there using trilinear interpolation. This scheme is unconditionally stable for any CFL number, allowing large timesteps.

Explicit diffusion — apply u* += νₜ ∇²u · dt (and same for θ') using second-order central differences. The viscosity is clamped to satisfy the diffusive CFL condition.

Buoyancy — add the buoyancy acceleration to the intermediate vertical velocity: w* += (g/T₀) · θ' · dt.

Stratification — damp the temperature anomaly proportional to vertical motion and stability: θ'* -= w* · (DALR − ELR) · dt.

Surface heat flux — inject Q_Ks · dt of temperature into the bottom fluid cell of each column above the terrain surface.

Boundary conditions — enforce no-slip (u=v=w=0) inside solid terrain cells and at the domain top. Horizontal boundaries are periodic.

Pressure Poisson equation — solve ∇²p = (∇·u*)/dt using Red-Black Successive Over-Relaxation (SOR) with ω ≈ 1.75. Boundary conditions are Neumann (∂p/∂n = 0) at terrain and top, periodic horizontally. The mean pressure is removed every N iterations (tunable via SOR mean-removal freq) to handle the null space of the Neumann problem.

Helmholtz projection — subtract the pressure gradient to enforce incompressibility: u = u* − ∇p · dt. This guarantees ∇·u = 0 (i.e. divergence-free velocity) in the fluid domain. Note that incompressibility is enforced ndirectly via pressure projection.

Grid and Boundary Conditions

The simulation domain is a rectangular box aligned with the terrain bounding box. The grid spacing is uniform: dx = domain_width / NX, dy = domain_height / NY, dz = domH / NZ. Terrain is represented as a solid mask: cells below (or at) the terrain surface are flagged solid, and all velocities inside them are forced to zero (no-slip). Heat is applied only to the first fluid cell above each terrain column.

Boundary	Velocity	Pressure	Temperature
Terrain surface	No-slip (u=v=w=0)	Neumann (∂p/∂n=0)	Heat flux Q_Ks applied
Domain top	w=0 (free-slip)	Neumann	Neumann (∂θ'/∂n=0)
X/Y sides	Periodic	Periodic	Periodic

Simplifications

Incompressibility (Boussinesq): density is constant except in the buoyancy term. Valid for vertical motions < ~2 km and speeds well below the speed of sound.
The ELR is assumed uniform (the same lapse rate throughout the domain), ignoring inversions and local stability variations.
The ground heat flux is spatially smooth (Gaussian-filtered land cover) and instantaneously applied, ignoring sub-surface heat conduction.
The SOR pressure solver converges to finite precision; residual divergence remains in the velocity field after the projection step.
Periodic lateral boundaries mean that thermals that leave one side re-enter from the opposite side — an artefact for small domains.

Limitations

Background wind — there is no mean horizontal flow. Orographic lift, rotor zones and wind shear are absent. All flow is purely thermally driven.
Moisture and clouds — no latent heat, no condensation, no cloud base or cumulus development. The simulation is dry throughout.
Coriolis effect — negligible at the spatial scales considered (< 10 km) but formally absent.
Radiation divergence — longwave cooling of the air column and radiative flux divergence within the atmosphere are ignored. Only surface heating is modelled.
Entrainment — the erosion of the thermal top by mixing with environmental air is handled implicitly by turbulent diffusion but no explicit entrainment parameterisation is used.
Synoptic forcing — large-scale pressure gradients, fronts, and mesoscale circulation are not represented.
Non-flat terrain top — the domain top is a flat horizontal plane, not a terrain-following coordinate. Thermals that reach the top are absorbed.

Parameters and Their Effects

Domain and grid

Parameter	Default	Effect
NX, NY	50, 50	Horizontal grid resolution. More cells → finer spatial structure, smaller thermals can be resolved, but runtime grows as NX·NY·NZ. Aspect ratio should match the terrain.
NZ	30	Vertical layers. Increasing NZ improves thermal height resolution but reduces dz, which tightens the diffusive stability limit and may require a smaller dt.
Domain height (domH)	4000 m	Vertical extent of the simulation box above the lowest terrain point. Should be at least as tall as the expected thermal top. Too small truncates rising thermals; too large spreads NZ layers thinly and coarsens vertical resolution.
Timestep (dt)	90 s	Integration step. Semi-Lagrangian advection is stable for any CFL, but buoyancy and diffusion impose practical limits. Reduce if the simulation diverges (velocities grow without bound). Larger dt speeds up computation but reduces accuracy.

Atmospheric physics

Parameter	Default	Effect
ELR	6.5 K/km	Environmental Lapse Rate. Controls atmospheric stability: `strat = DALR − ELR`. Higher ELR → less stable → thermals rise higher and faster. ELR = 9.76 K/km (= DALR) is neutral; ELR > DALR is absolutely unstable. Typical summer value: 6–8 K/km. Most impactful parameter for thermal height.
DALR	9.76 K/km	Dry Adiabatic Lapse Rate — the rate at which an unsaturated rising air parcel cools. This is a physical constant and should normally be left at 9.76 K/km. Reducing it artificially weakens the buoyancy restoring force, making the atmosphere behave as if moister (pseudo-adiabatic conditions). Only change for sensitivity studies.
Prt (Prandtl)	0.7	Turbulent Prandtl number: ratio of momentum to thermal eddy diffusivity, `κₜ = νₜ / Prt`. Higher Prt → less thermal mixing → temperature anomalies stay sharper and thermals are more intense. Lower Prt → more diffuse temperature field. Physically, the atmosphere has Prt ≈ 0.7; values of 1.0–1.5 reduce numerical diffusion of thermals in coarse grids.
Cs (Smagorinsky)	0.18	Scales subgrid turbulent viscosity and diffusivity. Higher Cs → more dissipation → thermals are weaker and more diffuse. Lower Cs → less dissipation → sharper, stronger plumes, but risk of numerical instability. Typical range: 0.10–0.25.
CFL factor	0.4	Caps the turbulent viscosity νₜ for diffusive stability: `νₜ ≤ CFL_factor · Δ² / dt`. Reducing this makes the diffusion scheme more conservative (less mixing per step); increasing it allows stronger diffusion but risks numerical instability. Values above ~0.5 may cause blow-up on coarse grids with large νₜ.
ω (SOR relaxation)	1.75	Over-relaxation factor for the pressure solver. Higher ω → faster convergence per iteration. Values above ~1.9 risk divergence. Reduce if pressure residuals are large.
Max iterations	60	Maximum SOR iterations per timestep. More iterations → better satisfaction of incompressibility (∇·u ≈ 0) but slower per step. Increase for large or fine grids.
SOR mean-removal freq	20	Remove the mean pressure from the pressure field every N SOR iterations. This handles the null space of the all-Neumann Poisson problem (pressure is defined only up to a constant). More frequent removal (smaller N) is safer but adds small overhead. Values of 10–30 are all equivalent in practice.

Solar heating

Parameter	Default	Effect
max_heat	15 K	Peak ground–air temperature difference under full insolation on grassland (tf=1). Directly scales all surface heat fluxes. Primary knob for overall thermal strength. Real values: 10–25 K on a warm summer day. Can be auto-calibrated from ERA5 data.
diffuse_frac	0.30	Fraction of solar energy arriving as scattered sky radiation (independent of direct beam and shadow). Higher diffuse fraction reduces the contrast between sunny and shadowed slopes, producing more uniform heating.
tau_hours	1.5 h	Reference ground thermal time constant, scaled per land cover type. Controls how quickly the ground temperature anomaly drives a heat flux into the atmosphere: `Q_Ks = dT / (τ · 3600)`. Smaller τ → faster, stronger flux for the same surface anomaly. Rocky terrain is already set to τ·0.35.
T_atm base	0.70	Base transmittance in the Kasten-Czeplak formula: `T_atm = T_base^(AM^T_exp)`. Represents clear-sky atmospheric transparency at the zenith (airmass = 1). Higher values (e.g. 0.80) correspond to a cleaner, drier atmosphere or high-altitude sites — more solar energy reaches the ground, driving stronger thermals. Lower values model haze, humidity, or pollution.
T_atm exponent	0.678	Exponent on the airmass in the transmittance formula. Controls how quickly transmittance drops as the sun approaches the horizon (high airmass). A higher exponent steepens the airmass penalty; a lower exponent gives a more uniform insolation throughout the day. The Kasten-Czeplak empirical value is 0.678.
σ slope (slope smoothing)	1.5 cells	Gaussian smoothing radius applied to the cos(incidence) field before computing direct solar irradiance. Reduces high-frequency artifacts from steep or noisy terrain normals. Larger σ → smoother heating patterns, less cliff/ridge edge artefacts. Set to 0 to disable smoothing.
σ land cover (lc smoothing)	4 cells	Gaussian smoothing radius applied to the land-cover thermal factor map. Blends abrupt forest/grassland boundaries, which are coarser than the simulation grid. Larger σ produces gentler spatial gradients in surface heating. Set to 0 for sharp land-cover boundaries.

Simulation time

Parameter	Default	Effect
Start / end time	10:00–16:00	UTC times at which scenarios are computed. Solar elevation and azimuth change significantly; the most intense thermals typically occur 11:00–14:00 local solar time.
Steps before snapshot	80	Number of NS timesteps integrated between each saved scenario. More steps allow the flow to develop further before sampling. Too few steps and the flow is under-developed; too many and transient features are lost.
Spin-up steps	40	Extra steps run before the first snapshot to let the flow develop from rest. Without spin-up, the first scenario shows an unrealistically quiet flow field.

Output Fields

Field	Symbol	Description
θ' surface	θ'(x,y)	Potential temperature anomaly at the first fluid cell above terrain. Positive = warm air rising; negative = cold sinking. Proxy for thermal trigger strength.
Vertical w	w(x,y)	Vertical velocity at the surface level. Positive (red) = rising air; negative (blue) = sinking. More directly related to what a paraglider feels than θ'.
Convergence	−(∂u/∂x + ∂v/∂y)	Negative horizontal divergence. Where surface winds converge, air is forced upward, triggering or feeding thermals. High convergence with high w indicates a strong thermal source.
θ'(x,y,z)	theta_3d	Full 3D temperature anomaly field, viewable on EW/NS/Z cross-section slices. Shows thermal plume structure with altitude.
w(x,y,z)	w_3d	Full 3D vertical velocity field. The updraft core visible in cross-sections corresponds to the exploitable thermal column.

Data Sources

Dataset	Used for	Source	Resolution
Elevation	Terrain mesh, solid mask, surface-flux injection point	AWS Terrain Tiles — Terrarium RGB encoding (IGN RGE Alti for France)	~7 m (zoom 14)
Land cover	Thermal factor tf, time constant τ per cell	ESA WorldCover 2021 (Cloud-Optimised GeoTIFF)	10 m
Calibration temperature	max_heat scale factor	ERA5 reanalysis via Open-Meteo archive API	~28 km
Sun position	Solar azimuth and elevation at each scenario time	Computed from lat/lon/date/time (astronomical formula)	—
Simulation output	Scenario snapshots 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, avoiding full-file downloads.

Tips for Realistic Results

Use ERA5 calibration to auto-set max_heat for the specific date and site.
Increase ELR toward 8–9 K/km on active convective days (cumulus days).
For a 10 km × 10 km domain, NX=NY=50 gives dx≈200 m — resolving thermals wider than ~400 m. Large thermals (mountain scale) are well captured; narrow convergence lines may be missed.
Run at least 80 steps before the first snapshot. At dt=90 s this is 2 hours of simulated time — enough for a well-developed convective boundary layer to form.
The simulation conserves neither energy nor temperature globally (periodic BCs + Neumann top allow net heat accumulation). Very long runs (>500 steps) may show unrealistic amplitudes.
Cross-section slices (EW/NS) are the most physically informative output: they reveal whether thermals detach from the ground, how high they reach, and whether there is organised convective structure (columns vs. sheets).