Pre-computed solar heating pipeline for paragliding thermal analysis · Terrain: SRTM 30 m · Land cover: ESA WorldCover 2021 · Calibration: ERA5 / Open-Meteo
The Thermal simulation is a diagnostic pre-computation pipeline — it does not solve fluid dynamics equations in real time. Instead, it walks through a sequence of physical steps (solar geometry → surface heating → thermal inertia → orographic flow → trigger zones) and produces a set of 2D maps, one per time step, exported as JSON to the Three.js viewer. The pipeline runs in seconds to minutes on the server.
Elevation data comes from SRTM 30 m (physics) and AWS Terrarium ~7 m tiles (visual mesh). Land cover is from ESA WorldCover 2021 (10 m), resampled to the terrain grid. The peak ground–air temperature difference is calibrated from ERA5 reanalysis (Open-Meteo API).
tf and time constant τ.
Exact solar azimuth and elevation are computed with pysolar (USNO astronomical algorithm) for each time step. The direct beam incidence on each terrain cell accounts for the slope and aspect:
A south-facing 30° slope receives significantly more direct beam than a flat or north-facing surface at the same insolation. The incidence field is Gaussian-smoothed (σ = slope_smooth, default 1.5 cells) to suppress artefacts from noisy DEM normals on steep terrain.
For each time step, a ray is cast from every terrain cell toward the sun. If any terrain point along that ray rises above the sun ray, the cell is in shadow and receives no direct beam — only diffuse sky radiation. This produces realistic cold patches on north-facing slopes and valley floors during morning and evening.
The solar beam weakens as it passes through the atmosphere, proportional to the airmass (path length). The Kasten-Young airmass formula and a Bird clear-sky transmittance model are used:
A low sun delivers much less energy than geometry alone would suggest.
Default: T_base = 0.70, T_exp = 0.678.
Cells in terrain shadow still receive scattered sky radiation. Modelled with an isotropic sky assumption and a sky view factor:
More scattering at low sun angles (higher diffuse fraction at
dawn/dusk). Default: diffuse_frac = 0.30.
The instantaneous ground temperature anomaly is the product of total insolation, the thermal factor, and the peak heating amplitude:
max_heat is the peak ground–air temperature difference
under full insolation on grassland (tf = 1). It is auto-calibrated from
ERA5 data or set manually.
The ground does not instantly reach its equilibrium temperature. Heat accumulates using a per-cell low-pass filter with time constant τ (land-cover dependent):
This causes thermals to lag behind the sun — peaking in the early afternoon rather than at solar noon. Bare rock (τ ≈ 0.5 h) heats and cools rapidly; forest (τ ≈ 3.75 h) responds slowly; water (τ ≈ 22 h) barely changes within a day.
Each grid cell has a thermal heating factor tf (scaling the
temperature anomaly) and a time constant τ (controlling how
quickly it heats and cools), derived from ESA WorldCover 2021:
| Land cover | tf (heating factor) | τ (h) | Effect |
|---|---|---|---|
| Bare rock / scree | 1.80 | 0.50 | Strongest heating, fastest response — spikes at noon |
| Built-up | 1.30 | 0.75 | Concrete/tarmac heat island |
| Grassland | 1.00 | 1.50 | Reference land cover |
| Cropland | 0.85 | 1.30 | Moderate inertia |
| Shrubland | 0.70 | 2.00 | Partial canopy — slower |
| Tree cover | 0.45 | 3.75 | Forest: weak and slow — canopy evaporation damps heating |
| Permanent water | 0.03 | 22.5 | Very high heat capacity — nearly constant temperature |
The ceiling at which a thermal parcel loses buoyancy is derived from the surface temperature anomaly and the difference between the Dry Adiabatic Lapse Rate and the Environmental Lapse Rate:
Higher ELR → less stable atmosphere → thermals reach greater altitude. At ELR = DALR the atmosphere is neutral and thermals rise indefinitely (convective instability).
Cells with a high trigger factor (ridges, land-cover boundaries) lose excess heat each time step once the temperature anomaly exceeds a release threshold, simulating periodic bubble detachment and surface cooling:
After the release step, the temperature anomaly field is advected using first-order upwind finite differences along a velocity field with two components:
The advection speed is computed dynamically each step — nearly idle at dawn (weak heating), peaking around solar noon. The time step is subdivided adaptively to keep the CFL number ≤ CFL_safety (default 0.9).
Two static diagnostic maps are derived from the simulation:
Clicking a cell in the viewer shows an estimated vertical speed profile derived from simple parcel buoyancy theory (no entrainment, no turbulence):
This is a rough order-of-magnitude estimate only. Real thermals are slower because this formula ignores entrainment of surrounding air (which dilutes buoyancy), turbulent drag, and the continuous-plume structure. Typical paragliding thermal speeds are 1–4 m/s; the formula may overestimate by a factor of 2–3 in strong conditions. Use it as a relative indicator, not an absolute number.
The scale factor max_heat is derived automatically from
ERA5 reanalysis data (Open-Meteo archive API) for the specific
simulation date and site. Clear-sky midday values of
surface_temperature − temperature_2m from the past 5 years are
averaged. This implicitly captures seasonal effects (soil moisture,
vegetation state) without explicit modelling.
For each terrain cell (i, j) and each time step, a ray is
marched from the cell toward the sun in the horizontal plane. The ray is
stepped in increments of one cell width; at each step the expected ray
height (based on the solar elevation angle and horizontal distance) is
compared with the DEM elevation. If any DEM point is higher than the
ray, the cell is flagged as shadowed. To avoid self-shadowing artefacts,
the ray starts one cell away from the source.
The ray direction is bilinearly interpolated on the DEM, so fractional-cell positions are handled correctly. Complexity is O(NX · NY · ray_length) per time step.
Surface normals are computed from the DEM using central finite differences:
One-sided differences are used at the domain boundary. The resulting cos(incidence) field is then Gaussian-smoothed with kernel radius σ = slope_smooth (default 1.5 cells) to suppress high-frequency artefacts from noisy DEM gradients on steep terrain.
Several fields (slope incidence, land-cover factor, trigger maps, advection flow) are spatially smoothed with a 2D Gaussian kernel of standard deviation σ (in grid cells). The kernel is applied separably as two 1D convolutions (x then y) for efficiency:
Boundary pixels are handled by reflection padding to avoid edge darkening. The kernel is re-normalised after truncation so the total weight sums to 1.
The ground temperature response to solar forcing is modelled as a first-order Infinite Impulse Response (low-pass) filter applied per cell at each time step:
This is the exact discrete solution to the ODE
d(dT)/dt = (dT_instant − dT) / τ, so it is exact regardless
of step size and unconditionally stable for any τ and dt.
Terrain convexity (ridges and peaks have negative curvature) is computed from the DEM Laplacian using a 5-point finite-difference stencil, after Gaussian pre-smoothing (σ = convexity_smooth, default 4 cells):
The land-cover boundary gradient is computed as the magnitude of the
gradient of the thermal factor field:
|∇tf| via central differences. Both maps are individually
Gaussian-blurred (σ = trigger_blur), normalised to [0, 1], then blended:
trigger = w·convexity + (1−w)·|∇tf|.
The temperature anomaly field is advected along the flow field (u, v) using first-order upwind finite differences — the spatial difference is taken from the upwind side to ensure stability:
Upwind differencing is monotone (no spurious oscillations) but
first-order accurate, introducing numerical diffusion proportional to
u·dx/2. This artificial diffusion slightly spreads thermal
plumes spatially, which is acceptable at typical grid spacings of
100–300 m.
The upwind scheme is stable only when the CFL condition is satisfied per
cell. The advection step is subdivided into substeps of size
dt_sub such that the maximum CFL across the domain is below
the safety threshold:
At low wind speeds or coarse grids a single substep suffices. Near solar
noon with strong heating and steep slopes (fast anabatic flow), 3–5
substeps may be used. The total advection displacement per full time
step is therefore always V_advect · dt regardless of the
substep count.
The horizontal divergence of the flow field is computed at each scenario step using central differences:
The weighting by dT means that warm-but-divergent cells and
cold-but-convergent cells both score low — only cells that are
simultaneously warm and in a convergent zone accumulate a high collector
score.
| Parameter | Default | Effect |
|---|---|---|
| max_heat | 15 K | Peak ground–air temperature anomaly under full insolation on grassland (tf = 1). Directly scales all surface heat fluxes. Auto-calibrated from ERA5. Primary knob for overall thermal strength. |
| diffuse_frac | 0.30 | Fraction of solar energy arriving as scattered sky radiation. Higher → more uniform heating across shadowed and north-facing slopes. |
| ELR (K/km) | 6.5 | Environmental Lapse Rate. Controls how high thermals rise: h_top = dT / (DALR − ELR). Higher ELR → thermals reach greater altitude. At 9.8 K/km (= DALR) the atmosphere is neutral. |
| Parameter | Default | Effect |
|---|---|---|
| T_atm base | 0.70 | Zenith atmospheric transmittance: T_atm = T_base^(AM^T_exp). Higher → more solar energy reaches the surface. Use 0.80–0.85 for high-altitude clear-sky days; lower for haze or humidity. |
| T_atm exponent | 0.678 | Controls how sharply transmittance drops at low sun angles (high airmass). Higher → stronger penalty near dawn/dusk; lower → more uniform insolation through the day. |
| σ slope (incidence smoothing) | 1.5 cells | Gaussian blur on the cos(incidence) field. Suppresses artefacts on cliff edges and noisy DEM normals. Set to 0 to disable. |
| Clear-sky threshold | 600 W/m² | Minimum direct radiation for an ERA5 hour to count as "clear-sky" in the max_heat calibration. Lower includes more marginal hours; higher keeps only sunniest midday observations. |
| Parameter | Default | Effect |
|---|---|---|
| τ hours (reference) | 1.5 h | Reference thermal time constant for grassland (tf = 1). Other land covers scale from this: rock × 0.35, forest × 2.5, water × 15. Smaller τ → faster, more peaked heating. |
| σ land cover (lc smoothing) | 4 cells | Gaussian blur on the land-cover thermal factor map. Blends abrupt forest/grassland boundaries. Set to 0 for sharp LC boundaries. |
| Parameter | Default | Effect |
|---|---|---|
| Convexity smooth σ | 4 cells | Gaussian blur on the DEM before computing curvature (Laplacian). Larger σ → peaks and ridges blend; smaller → sharper but noisier curvature map. |
| Trigger blur σ | 3 cells | Post-processing blur on convexity and LC boundary maps before blending. Controls spatial extent of trigger zones around ridges and forest edges. |
| Convexity weight | 0.60 | Blend: trigger = w·convexity + (1−w)·LC_boundary. Higher w favours terrain peaks/ridges; lower favours forest/grassland transitions. |
| Parameter | Default | Effect |
|---|---|---|
| V_advect_max | 1000 m/h | Ceiling for anabatic flow speed. Scales the entire advection velocity field. Higher → heat moves faster upslope. |
| Ref slope | 0.30 | Characteristic slope (tan) normalising the terrain gradient in the advection scaling. Lower → advection activates on gentler terrain. |
| k_collect | 1.0 | Strength of collector-suction term that pulls heat down the ΔT gradient toward trigger zones. 0 = pure orographic flow. |
| Flow smooth σ | 5 cells | Gaussian blur on the advection vector field before applying it. Reduces checkerboard artefacts from the upwind scheme. |
| CFL safety | 0.9 | Target CFL for adaptive substep splitting. Lower → more substeps, more stable. Values above 1.0 are unstable. |
| Parameter | Default | Effect |
|---|---|---|
| K_release | 0.3 h⁻¹ | Rate at which excess heat is removed from trigger-zone cells once the threshold is exceeded. Higher → thermals release more aggressively each step. |
| Release threshold | 1.5 K | Minimum temperature anomaly above which release cooling activates. Acts as a noise floor: only cells that have genuinely accumulated heat release. |
| Parameter | Default | Effect |
|---|---|---|
| Diurnal amplitude | 5 K | Half the daily temperature swing for gap-filling missing ERA5 hourly values: T(h) = T_mean + A·cos(2π(h − h_max)/24). Used only when ERA5 data has gaps. |
| Peak hour | 14 h (local) | Hour of the daily temperature maximum in the cosine gap-fill model. |
| Field | Symbol | Description |
|---|---|---|
| Thermal heat | dT(x,y,t) | Ground temperature anomaly at each time step (filtered by display threshold). Positive = warm surface. Primary indicator of thermal strength and location. |
| Thermal top | h_top(x,y,t) | Estimated ceiling altitude (m AGL) where a parcel from this cell loses buoyancy. Derived from dT and atmospheric stability. |
| Daily heat | max_t dT(x,y) | Maximum temperature anomaly over all time steps — the single hottest moment of the day per cell. |
| Trigger zones | trigger(x,y) | Static terrain convexity + land-cover boundary score. High values (bright) identify where thermals preferentially detach from the surface. |
| Collecting zones | conv(x,y) | Accumulated heat-weighted horizontal convergence. High values mark areas that funnel warm air into trigger points. |
| Flow arrows | u, v (x,y,t) | Advection velocity field (orographic + collector suction). Shows the direction thermal heat moves across the terrain at each time step. |
| Dataset | Used for | Source | Resolution |
|---|---|---|---|
| Elevation — simulation | Shadow casting, slope normals, thermal tops, orographic flow | SRTM via AWS elevation tiles (Mapzen/Tilezen) | ~30 m |
| Elevation — visual mesh | 3D terrain surface rendered in the browser | AWS Terrain Tiles, Terrarium RGB encoding (IGN RGE Alti for France) | ~7 m |
| Land cover | Thermal factor and time constant per cell | ESA WorldCover 2021 | 10 m |
| Sun position | Azimuth + elevation at each time step | pysolar (USNO astronomical algorithm) | — |
| Calibration temperature | max_heat scale factor | ERA5 reanalysis via Open-Meteo archive API | ~28 km |
The SRTM 30 m grid is used for all physics computations: slope normals, shadow ray-casting, thermal tops, and orographic flow. A coarser grid keeps the computation and exported JSON compact.
The AWS Terrarium ~7 m tiles are fetched once and used only to build the
visual terrain mesh in Three.js. At zoom level 14 the underlying source
for France is IGN RGE Alti (LiDAR-derived), giving ~4–5× higher surface
detail than SRTM. The RGB elevation encoding is:
elevation = R×256 + G + B/256 − 32768, delivered in Web
Mercator (EPSG:3857) and reprojected to WGS84 before use.
Physics overlays (heatmap, flow arrows, etc.) are rendered on the
coarser 30 m grid but remain visually correct because they cover the
same geographic extent and use depthWrite: false, floating
just above the terrain surface.