Every drop, every three metres

Skybrudd answers one question: what happens to this place if a cloudburst hits it? You pick a spot on the map, choose how much rain, and a minute later you can watch the water run down the streets, pool against walls and fill the underpass. The water is computed on your own GPU, in the page, with no server doing the work.

This piece is about how that water is calculated. It covers what goes in, the equations, the numbers I had to get right, and the checks a result has to pass before anyone else gets to see it.

The ground, to the metre#

Water goes where the ground sends it, so everything starts with the terrain. Kartverket publishes a laser-scanned terrain model of most of Norway at 1 m resolution. An area in Skybrudd is a 16:9 rectangle. The result file always holds 256 × 144 cells, and the solver works on cells half that size: 512 × 288, 3 m each for the default 1.5 km area.

On top of the terrain go the things that change where water flows:

• Buildings from OpenStreetMap raise the bed by their height. A building with four storeys and no height tag is 13 m tall: 3 m per storey plus 1 m of roof. Rain on a roof runs off its edges.
• Roads are graded. Real roads are cut into slopes and built up over dips, so they are level across and smooth along. Averaging laser points smears that out and tilts the road with the hill, which lets water cross it sideways. Each drivable road's centreline profile is smoothed over 20 m and burned back in, with a cut-or-fill shoulder. No change is larger than 2.5 m, so a smoothed road never dams a valley.
• Land cover sets the roughness and how much water the ground can take.
• Drains come from NVDB, the road authority's database: gullies take up to 20 l/s, grates 30 l/s, and culverts are pipes under the embankment.

Going from 6 m to 3 m cells uncovered one problem. A footprint that rasterises with a missing cell in the middle leaves a 3 m pit walled in by roofs. Every drop that lands on those roofs drains into the pit, which fills to the eaves: 5 m of water in a hole the size of a parking space. Any patch of ground smaller than 50 m² that buildings wall in on every side now counts as roof.

The rain, shaped like a storm#

There are two scenarios. A cloudburst puts 0.5 to 6 hours of rain on the area in 10-minute steps, then runs 6 dry hours so the water can drain away. The year scenario takes the place's total precipitation for last year and drops it in 24 hours.

Rain never falls evenly, so the intensity follows a Chicago design storm around a peak at $t_p$:

Each step is the shape averaged over 20 sub-steps, then scaled so the total comes out exact. A climate factor multiplies the total; 1.4 is the usual allowance for cloudbursts towards 2100. Here is 49 mm in an hour, times 1.4:

The soil fills up#

The first version let every pervious cell soak up water at a fixed rate forever. For a short storm that is roughly right. For a year's rain in a day it is badly wrong, because a forest at 10 mm/h would swallow 240 mm and still ask for more.

Now infiltration follows Horton's curve. Dry ground takes water fast, and the rate decays towards a steady $f_c$ (5 mm/h by default) until the soil has taken its storage $S$. After that, only slow percolation continues:

The trick that makes this cheap on a GPU is that $f$ depends only on the time since the rain started, so every cell shares one curve and no cell needs a soil state of its own. A cell takes its land-cover factor times the curve, and only the water actually standing on it. The time the soil fills, $t_\text{full}$, is found once by bisection, and the integral is exact between any two times, so the step size never changes how much water soaks in. On the day of the run, seNorge's ground state scales $f_c$ and $S$ by the unsaturated fraction, and by another 0.1 when the ground is frozen.

Water from beyond the edge#

This was the biggest gap. A town in a valley also gets the rain that fell on the hills above it, and the first version simply did not. Its edges only let water out.

Now the model searches 20 km around the area on Terrarium elevation tiles (zoom 11, about 60 m), routes every pixel downhill with a priority flood that fills pits so even lakes drain somewhere, and adds up the area of every pixel whose path enters the box. That gives a list of catchments, each with an area, a place where it crosses the edge, and a time of concentration $T_c$: the longest travel time to the edge, from a velocity at every pixel.

Each catchment's rain, less the same Horton losses, becomes a Clark unit hydrograph. The HEC time–area curve spreads the excess over $T_c$, and a linear reservoir with $K = 0.6\,T_c$ smooths it:

A cloudburst is local, so over catchments larger than 1 km² the rain is reduced by the Flood Studies Report's areal reduction factor. The largest catchment above Balderheia is 18.0 km² with $T_c$ = 4.3 h, which reduces the storm above by 11 %. Its runoff enters through two edge cells:

This changes the answer a lot. In most of the eleven demo areas, more water now comes in across the edge than falls on the area itself:

Rivers get no base flow. The laser terrain already shows them at their usual level, so only the storm's rise is added.

Shallow water, local-inertial#

The solver is a local-inertial scheme from the LISFLOOD-FP family (Bates et al., 2010). It solves the 2D shallow-water equations without the advection term. That holds for floods over land, where friction dominates, and it means each face of the grid can be updated on its own, which suits a GPU.

Depth $h$ lives in cell centres, and unit discharge $q$ (m²/s) on the faces between them. With the water level $\eta = z + h$, the depth that can cross a face is

so only the water above the top of a kerb or wall gets across it. The momentum update blends the old discharge with its two neighbours along the same axis (de Almeida et al., 2012), treats gravity explicitly and Manning friction semi-implicitly:

Then come three rules that keep it honest. The velocity is capped so a waterfall off a roof stays at a Courant number of $\alpha$ or below. Outflow is limited so no cell gives away more water than it holds. Continuity adds up the rest:

where $r$ is the rain, $f$ the infiltration and $d_i$ what the drains take from the cell.

The time step follows the fastest gravity wave in the domain, and is never longer than a third of a second per metre of cell:

A 2 m deep pond in 3 m cells gives $\Delta t = 0.5 \cdot 3 / \sqrt{9.81 \cdot 2} \approx 0.34$ s, so a seven-hour storm run takes about 40 000 steps. Water leaves at the edges at normal depth, $q = h^{5/3}\sqrt{S}/n$. Rain, inflow, storage, outflow, infiltration and drains are tracked separately and balance to about $10^{-14}$ of the total.

Three bugs that looked like floods#

Explicit schemes fail quietly. The water still balances and the map still looks wet, but the numbers are wrong. When I went back to check whether the water was calculated correctly, I found three faults.

The Courant number was too high. At $\alpha = 0.7$, deep pools sloshed by metres between neighbouring cells, and the fastest velocity in a storm run reached 275 m/s. LISFLOOD-FP recommends 0.5–0.7, and 0.5 is stable in every test I have, deep pools included.

The time-step cap was fixed. A cap of 5 s regardless of cell size let thin sheet flow grow a checkerboard. On a plane, the peak depth of a test storm came out at 16 cm where the kinematic wave gives 9 to 10. The cap is now $\Delta x/3$, and a test compares the plane against the kinematic wave at 3, 6 and 12 m cells.

Deep, smooth water oscillated. The plain local-inertial scheme grows oscillations on deep water with little friction. De Almeida's blending of each face with its neighbours, $\theta = 0.7$, damps them while keeping the mass exact. A lake-at-rest test, still water over rough ground with a building standing in it, must give a fastest velocity below $10^{-9}$ m/s.

One step, nine kernels#

The same Rust code runs on the server, in the desktop app and in the browser. Compiled to WebAssembly, it runs in a worker: it fetches the terrain, buildings, drains and catchments, builds the grids, and hands them to a WGSL port of the solver. Every time step is nine compute passes:

// one time step, in dispatch order
ctrl      // one thread: next Δt from the CFL, h_max and the time left in the frame
old       // copy every face's discharge, so the blend reads the previous step
fluxes    // the momentum update of every face, with the velocity cap
limit     // per cell: scale outgoing fluxes so nothing goes below zero
depths    // continuity, rain, Horton, max-depth and arrival trackers
inflows   // pour each catchment's hydrograph into its edge cells
sinks     // gullies and grates take up to their capacity
culverts  // move water through the pipes under the roads
network   // share the municipal pipe capacity among the wet road cells

The hard part was precision. GPUs compute in 32-bit floats, which carry about seven significant digits. At 400 m above sea level, that leaves the bed with about 0.03 mm of resolution, which is enough to lose a film of rain. So the bed is uploaded relative to the lowest point in the area, and depths and the clock are carried as pairs of floats, a high part and a low part, that together hold more than 40 bits.

The GPU and the CPU solver give the same answer. On the same storm, mean max depth differs by 0.1 mm and every volume total agrees within 0.01 %.

Solve fine, store coarse#

Ditches, kerbs and the gaps between houses decide where water goes, so the solver runs at twice the file's resolution. Every 2 × 2 block of solver cells then becomes one file cell, by rules that keep what matters from each:

Halving the cell size costs about eight times the work: four times as many cells and twice as many steps, since $\Delta t$ scales with $\Delta x$. On the GPU that is still under five seconds for a storm. The difference in the results is large. Narrow channels between buildings now carry water that 6 m cells averaged away. Here are the eleven demos before and after the whole review: finer cells, Horton soil, inflow and the stability fixes.

The result file is plain binary: a JSON header, then terrain, buildings, land cover, the three trackers, and one frame of depth and velocity every ten minutes.

"FLOM" · u32 version · u32 json_len · JSON metadata · pad to 4
f32 terrain[N] · f32 building_h[N] · u8 landcover[N] · pad
f32 max_depth[N] · f32 arrival_h[N] · f32 max_q[N]
per frame:  u16 depth_mm[N] · i16 vx_cms[N] · i16 vy_cms[N]      N = 256 × 144

Trust, but run it again#

Anyone can share a run to the public map, and the run was computed on their machine, so the server does not take the file on trust. Before it keeps one, it checks:

• the grid matches the centre and width, in Norway, and the rain series matches the scenario;
• frames, totals and trackers agree, and the water balances within 0.1 %, inflow included;
• the soil storage and every catchment's inflow volume are what the model gives;
• the terrain matches Kartverket's (median at most 1.5 m off), and the catchments, worked out again from the terrain around the area, match the ones in the file;
• the last rain step and a random one are simulated again from the frame before, and must land within 5 % of the volume, 10 % of the flooded cells and 3 cm of mean depth.

For a refined run the server rebuilds the 3 m grids from Kartverket's terrain and the buildings and roads listed in the file, and spreads each file cell's water over its four solver cells by level-pooling, before it steps the run forward.

What it does not model#

The underground pipe network is the biggest gap. Gullies drain at full capacity for the whole run, which is optimistic once the pipes are full, and the municipal network is an assumed rate per area because its layout is not open data. Buildings have no doors, so water piles up against walls and flows around them. A very deep pool usually means a culvert or underpass that is missing from NVDB. Erosion, debris and blocked inlets, which do much of the damage in real cloudbursts, are not modelled at all.

The flooded-area figure currently counts any cell with more than 10 cm of water, rivers and lakes included. Now that upstream water fills the river channels, that figure flatters the flood in valleys with a big river. The map's risk layer already ignores open water. The summary number will follow.

You can try your own street at skybrudd.dagthom.as.