Mass-Consistent Wind Fields

Problem

Weather APIs return a single wind vector for an entire location. That forecast is spatially uniform — the same speed and direction everywhere. In flat terrain this is a reasonable approximation. Over hills, valleys, and ridgelines it is not.

Wind accelerates as it is compressed over a ridge and decelerates as it spreads into a valley. A uniform field ignores this entirely. Particles advected through such a field all move in lockstep, regardless of the terrain beneath them. The result looks wrong and conveys no useful local information.

Model

The adjustment is based on conservation of mass. Air is incompressible at these scales, so the volume flux through each column of atmosphere must balance. Where terrain rises and the column is shallow, wind must speed up. Where terrain drops and the column is deep, wind slows down.

Formally, we solve for a Lagrange multiplier field λ that enforces this constraint. The governing equation is an elliptic PDE:

The right-hand side measures how much the initial uniform wind diverges due to terrain slope. The diffusion coefficient D = (H − h)² is the square of the column depth — deeper columns diffuse the correction more readily. Boundary conditions are λ = 0, meaning no correction at the domain edges where the forecast is trusted.

Once λ is known, the corrected wind is recovered by adding the gradient of the multiplier field, scaled by the local diffusion:

This is a Sasaki-type variational method. It produces the minimum adjustment to the initial wind that satisfies mass consistency — the corrected field stays as close to the forecast as possible while respecting the terrain.

Implementation

The solver runs server-side in a Next.js API route. Two functions handle the work.

Elevation grid

A 20×20 grid of elevation values is fetched from the Open-Meteo elevation API, centred on the queried location with a 5 km radius. The grid is cached in memory so repeated queries to nearby coordinates avoid redundant API calls. Grid spacing is converted from degrees to metres using the local latitude for the east-west direction.

Gauss-Seidel solve

The elliptic PDE is discretised on the elevation grid using central differences. Diffusion coefficients are averaged at cell interfaces (harmonic-mean stencil). The system is solved iteratively with Gauss-Seidel relaxation — up to 300 iterations, converging when the maximum update falls below 10⁻⁶. In practice, convergence takes 40 to 80 iterations for typical terrain.

Client interpolation

The corrected wind grid is returned to the client as two flat arrays (u, v) with grid bounds. The wind map component uses bilinear interpolation to sample wind at each particle position, producing smooth advection across the entire domain. Particles that drift outside the grid bounds are recycled to random positions.

Graceful fallback

If the elevation fetch or solve fails, the API returns the raw forecast wind without a grid. The client detects the missing grid and falls back to uniform advection — visually simpler, but still functional. No error is shown to the user.