Flood forecasting with machine learning models in an operational framework

In plain words

Rivers flood, and the people downstream need to know early and where the water will reach. This paper describes the flood warning system Google runs for large rivers, and the machine learning parts inside it. The system predicts how high a river will get at a measuring station (a gauge) over the next day or two, using the station's recent readings, readings from stations upstream and rainfall estimates from satellites. When a predicted level crosses a warning level for that station, a second set of models estimates which streets and fields will be under water, and how deep.

The inundation part comes in two versions. One only says whether a spot floods, using past floods as examples. The other also estimates water depth, but needs a ground elevation map of the area. The paper also describes the system in use in India and Bangladesh during the 2021 monsoon season, when alerts went out by search results, phone notifications and map icons.

TL;DR

Google's operational river flood warning system uses two learned stage-forecast models (a per-gauge linear model and a shared LSTM) and two learned inundation models (Thresholding, and the new Manifold model that also gives depth). The LSTM with precipitation is the one used for most gauges. Author-reported skill is high: median NSE 0.987 and persistent-NSE 0.671 for the LSTM, and median F1 of about 69% for inundation extent in 1-year cross-validation. Code and the main stage data are not public, so the card is the only route to re-implement the models.

Problem

Flood warnings need two numbers: river stage at the target gauge over the lead time, and the area or depth that would flood §1 §2. Physics-based hydraulic models need a digital elevation model (DEM), bathymetry estimates and per-gauge calibration; the authors report that the hydraulic model was discarded because it was less accurate and less scalable than the ML models, and much costlier to run §2.3 p. 10. They also cite an earlier comparison in which the hydraulic model needed about eight times more manual work §4 p. 27. Discharge, which conceptual hydrology models often use, is scarce and uncertain for stream gauges, so the system forecasts stage directly §2 p. 5. Large rivers are the target: slow response, and stage data available in near real time §2 p. 5.

Method

The system has four subsystems: data management, stage forecasting, inundation modelling and alert distribution §2 fig. 1.

Data management §2.1 p. 6–7. Stage readings are hourly for the deployments described. Each reading is checked for typing errors (for example a misplaced decimal point), tested against record statistics and against the previous reading, and suspicious values are removed or kept by ad-hoc criteria. Short gaps are filled by linear interpolation. Precipitation comes from IMERG Early Run, which has 0.1° and 30-minute resolution and a stated latency of 12 hours. Negative and missing values are removed, an upper rain-intensity threshold of 200 mm/h is applied (what happens to values above it is not stated), and the rain is averaged over each gauge's drainage area and integrated to hourly.

Stage forecast, Linear model §2.2a fig. 2a p. 7. Inputs are past stages at the target gauge and its upstream gauges (typically 2–5 gauges) over a 72-hour lookback. A multiple linear regression is trained with MSE plus an L2 penalty on the coefficients. One model is trained per gauge and per hourly lead time up to that gauge's maximal lead time; a 24-hour lead time means 24 models for that gauge.

Stage forecast, LSTM §2.2b fig. 2b tab. 1 p. 8. The model has two LSTMs with a state handoff between them. The hindcast LSTM is sequence-to-one. It reads the past at each hour: IMERG precipitation, the near-real-time stage, and five features made by a gauge-specific linear layer from the upstream stages over an upstream lookback. Its final cell and hidden states go through a fully connected layer, and the result initialises the forecast LSTM. The forecast LSTM runs one step per hour of lead time up to the maximum. Its outputs go to a dedicated head per time step that produces a countable mixture of asymmetric Laplacians (CMAL), following Klotz et al. §2.2b. The loss is the negative log-likelihood of observed stage at each forecast step, summed over forecast steps only; the hindcast contributes only through its initial state. All LSTM weights are shared across gauges (regional calibration). Only the linear upstream combiner is gauge-specific. Training uses one maximal lead time for all gauges, the largest gauge-specific value §2.2b.

Warning threshold §2.2 p. 10. For each gauge, the maximal forecast stage between the forecast's current time and the gauge's maximal lead time is compared with a gauge-specific warning threshold. If it exceeds the threshold, an alert is issued and that maximal stage drives the inundation model. External stage forecasts from national agencies go through the same check §2.2 p. 10.

Inundation, general §2.3 p. 10. The inundation models map a forecast stage to a flood extent over an area of interest (AOI) around each gauge. Stage is taken as steady for each forecast, so earlier stages are ignored.

Thresholding model §2.3a fig. 3a p. 10–11. Each AOI pixel is treated as a binary classifier of wet (positive) against dry, with gauge stage as input. For each pixel the model learns a stage threshold by repeatedly picking the threshold with the best ratio of true wet to false wet events among events above it, then discarding those events and continuing. It stops when the ratio drops below a global minimal ratio parameter and keeps the previous threshold. A global parameter, named minimal ratio in one sentence and minimal threshold in another, is tuned for an F-score. The authors show that for each value of the minimal ratio, no other set of pixel thresholds has both higher precision and higher recall, so the set is Pareto optimal §2.3a p. 11. If the forecast stage exceeds all training stages, the extent starts from the largest historical flood and grows in every direction, by a distance linear in the stage excess §2.3a p. 11.

# Thresholding model, one AOI pixel  [§2.3a]
# Input:  stages (N,), wet (N,) bool for this pixel, minimal_ratio
# Output: threshold stage for this pixel

def fit_pixel_threshold(stages, wet, minimal_ratio):
    # Start with every training event  [§2.3a]
    remaining = all_events
    # Keep the last accepted threshold; none yet  [§2.3a]
    accepted = not stated
    while remaining is not empty:
        # Find the threshold with the best wet/false-wet ratio among
        # the remaining events above it  [§2.3a]
        candidate, ratio = best_ratio_threshold(stages, wet, remaining)
        # Stop once the ratio falls below minimal_ratio  [§2.3a]
        if ratio < minimal_ratio:
            break
        accepted = candidate
        # Discard events above the accepted threshold, then iterate  [§2.3a]
        remaining = remaining where stage < candidate
    return accepted

Manifold model §2.3b fig. 3b p. 11–12. This model needs a DEM for the AOI and produces both an extent and a depth map. It has two parts.

Flood extent to water height §2.3b.1 p. 11–12. Given a DEM and a wet/dry map, the goal is a smooth water surface with no interior minima or maxima inside flooded regions. The paper uses a heuristic: trace the boundary of each flooded region (OpenCV findContours), read DEM heights on the boundary, interpolate inside by solving the Laplace equation, and remove DEM outlier pixels that cause high Laplace tension. The result is stored on a grid where each cell covers 32×32 DEM pixels.

# Flood extent to water height  [§2.3b.1]
# Input:  dem (H, W) in metres, wet (H, W) bool
# Output: water_height on a 32x32-DEM-pixel grid, metres above sea level

def extent_to_height(dem, wet):
    # Trace the boundary of each flooded region  [§2.3b.1]
    boundaries = find_contours(wet)
    # Read DEM heights along those boundaries  [§2.3b.1]
    boundary_heights = dem[boundaries]
    # Drop DEM outlier pixels with high Laplace tension  [§2.3b.1]
    boundary_heights = remove_outliers(boundary_heights)  # rule not stated
    # Interpolate inside each region with the Laplace equation  [§2.3b.1]
    surface = laplace_interpolate(boundary_heights, wet)
    # Store on the coarse grid: one cell per 32x32 DEM pixels  [§2.3b.1]
    return downsample(surface, factor=32)

Gauge stage to depth §2.3b.2 p. 12. Precomputation: run the Thresholding model on every historical event, giving an extent map per gauge stage, then run the extent-to-height step on each to get one water-height map per training event. At run time, for a forecast stage, the model interpolates per pixel, piecewise-linearly between the water-height maps of the two nearest training stages. A pixel is dry if its water height is below the DEM and wet otherwise; for wet pixels, depth is water height minus DEM. If the stage is above all training stages, the stage difference is added to every pixel of the highest map.

# Manifold model at run time: stage to extent and depth  [§2.3b.2]
# Input:  stage (scalar), height_maps (per training event),
#         train_stages (sorted), dem (H, W)
# Output: wet (H, W) bool, depth (H, W) in metres for wet pixels

def manifold_inundation(stage, height_maps, train_stages, dem):
    # Interpolate each pixel between the two nearest training maps  [§2.3b.2]
    water_height = piecewise_linear(stage, train_stages, height_maps)
    # Above the top training stage, add the stage gap to every pixel  [§2.3b.2]
    # (applies only when stage exceeds max(train_stages))
    # Wet if water height reaches the DEM  [§2.3b.2]
    wet = water_height >= dem
    # Depth of wet pixels: water height minus DEM  [§2.3b.2]
    depth = water_height - dem
    return wet, depth

Hydraulic model, benchmark only §2.3c fig. 3c p. 13. This is a St. Venant-based model under steady-state assumptions, run offline for many upstream discharges with fixed roughness and downstream slope. Results are stored in a hash table keyed by gauge stage, and the nearest stage's map is returned. River bathymetry is estimated by fitting three equal-sided trapezoids to cross sections. It was used for benchmarking only and not deployed §3.2 p. 23.

Training and validation of inundation §2.3 p. 15. Extent maps come from Sentinel-1 SAR ground-range-detected data, one image every few days per AOI, and each pixel is classed wet or dry by a Gaussian-mixture classifier. Sentinel-2 maps calibrate it: NDWI = (B3 − B8)/(B3 + B8), wet when the value is above 0 §2.3 p. 15. DEMs are built per AOI from high-resolution optical satellite imagery by stereographic imaging, and retrained yearly where floods reshape the terrain §2.3 p. 15. The Manifold model is trained and validated on extent only, since no ground-truth depth exists §2.3 p. 15. Manifold is preferred for a gauge when a DEM exists, because it also yields depth §2.3 p. 15.

Stage model training §2.2 p. 9. One-year leave-out cross-validation: each fold holds out one year. For operation the models are retrained on the full data.

Alerts §2.4 fig. 4 p. 16–19. Alerts give the forecast stage change at the gauge, the inundation map and depth where available. They go to authorities, relief agencies and the public. Public channels are Google Search, phone push notifications to phones in the forecast flood area, and Google Maps icons. Alerts were sent in nine languages for India and Bangladesh §2.5 p. 19.

Hyperparameters

Results

Author-reported. Setup: gauges in India and Bangladesh; 1-year cross-validation over 2014–2020 for stage §3.1 p. 21; 1-year cross-validation and leave-extreme-out over 2016–2020 for inundation §3.2 p. 24. Persistent-NSE uses the last observed stage as the baseline forecast §3.1 p. 21.

Stage forecast, 167 gauges tab. 2 §3.1 p. 21–22:

Source: tab. 2. Differences marked significant by the paper's Wilcoxon signed-rank test on paired samples tab. 2 p. 22. The text adds that the LSTM beats the Linear model on Persistent-NSE for 95% of gauges, and that precipitation improved the LSTM for 89% and the Linear model for 65% of gauges §3.1 p. 22. Persistent-NSE rises with watershed area; medians for the three bins (<10,000, 10,000–100,000 and >100,000 km²) are plotted with 10–90% quantile bars, without printed values fig. 9 p. 22.

Inundation, 126 AOIs with a DEM tab. 2 §3.2 p. 24–25:

Thresholding has a higher F1 than Manifold in 70% (CV) and 57% (leave-extreme-out) of AOIs, and the paper reports the median differences as significant by Wilcoxon signed-rank §3.2 p. 25 tab. 2. The Manifold model also gives depth, which is not evaluated §3.2 p. 15.

Thresholding against the hydraulic model, 11 AOIs, trained 2016–2019 and validated on 2020 tab. 3 §3.2 p. 23–24:

Source: tab. 3. Thresholding is higher on 9 of 11 AOIs §3.2 p. 23.

Deployment, 2021 monsoon, India and Bangladesh §2.5 p. 18–19 (author-reported):

• 376 target gauges; stage models on 167 of them, the rest use external forecasts.
• Inundation maps for 233 gauges: 174 Thresholding, 59 Manifold. A further 143 gauges had alerts without maps.
• Coverage of about 287,000 km² of flood-prone area. §2.5 p. 18
• About 115 million notifications sent, reaching about 22 million people, from 8 July to 10 September 2021 fig. 6 p. 19.
• Stage benchmark: median RMSE 0.1 m at 24 h and 0.2 m at 48 h, against 0.3 m and 0.5 m for US NWS forecasts in a 2017 study. The paper warns the basins differ §4 p. 26.

Paper vs code

No official code. The paper's code-availability statement says the software was developed at Google and the code is proprietary [code availability] p. 30. Searched: the paper and abstract link no repository and no project page; GitHub search for "flood forecasting" in the google-research organisation returned google-research/flood-forecasting. That repository does not cite this paper, implements a different LSTM family (Mean-Embedding and Handoff forecast models) and has no Linear, Thresholding or Manifold model, so it is not treated as this paper's code. A GitHub search for "flood inundation manifold forecasting" returned nothing.

Pitfalls

• Persistent-NSE is the honest yardstick. Raw NSE is near 0.98–0.99 for every stage model, because large rivers change slowly and the forecasts use past stage as input §3.1 p. 22. Compare models with Persistent-NSE, and note that its lower tail (negative values for some gauges) is wide tab. 2.
• Stage LSTM loss only sees forecast steps. The hindcast LSTM only sets the forecast LSTM's initial state; the loss runs over the forecast outputs §2.2b p. 8. The text does not say what the forecast LSTM receives at each step besides the state; fig. 2b shows no per-step input. Treat that input as not stated.
• Shared maximal lead time. The LSTM trains with one maximal lead time for all gauges (the largest gauge value), while operation uses each gauge's own maximal lead time §2.2b tab. 1. How the extra steps are handled at inference is not stated.
• Which forecast value drives the alert. Alerts use the maximal forecast stage between the current time and the gauge's lead time §2.2 p. 10. The paper does not say whether that stage is the mean, median or another quantile of the CMAL distribution, nor the number of mixture components §2.2b p. 8.
• Linear model count. One model per gauge per hourly lead time, so the number of models grows with the maximal lead time §2.2a p. 7.
• Stage quality control has unstated criteria. Suspicious readings are removed or kept by "a set of ad-hoc criteria", and "short" gaps are interpolated with no length given §2.1 p. 7. Implement these with your own thresholds, and expect the results to depend on them.
• Precipitation latency. The source states a 12-hour IMERG latency but does not say how the real-time system handled it §2.1 p. 7.
• Thresholding parameters are missing. The minimal ratio, minimal threshold and the expansion coefficient for stages above history are not given §2.3a p. 11–12. The candidate-threshold grid used in the search is not stated either. The "Pareto optimal" claim is for the fitted per-pixel thresholds §2.3a p. 11.
• Manifold outlier removal is unspecified. The rule for dropping DEM outliers that cause high Laplace tension is not stated §2.3b.1 p. 12. The 32×32 grid is a coarse depth map: depth is only as fine as that grid §2.3b.1 p. 12.
• Depth is not validated. No ground-truth depth exists. The Manifold model is trained and tested only on extent; the authors say accurate extent suggests reasonable depth, and that depth accuracy needs more validation §2.3 p. 15.
• Depth above history is a heuristic. Above the highest training stage, the stage gap is added uniformly to the highest map §2.3b.2 p. 12. This is not checked against floods above the record.
• Resolution inconsistency. Maps are computed at 16×16 m² but evaluated at 64×64 m² §2.5 p. 19. The Discussion describes 16 m alerts §4 p. 25, which the metrics do not cover.
• Population figure conflicts. The abstract says the system covers more than 350 million people; Section 2.5 says "close to 245,000,000" [abstract] §2.5 p. 18. Use the body figure and mark the conflict.
• Years differ between tasks. Stage training uses 2014–2020, inundation 2016–2020, and the hydraulic comparison trains on 2016–2019 and validates on 2020 §2.5 §3.1 §3.2. Do not pool them.
• Labels are classifier output. Inundation labels come from a Gaussian-mixture classifier on Sentinel-1, calibrated against Sentinel-2 NDWI > 0 §2.3 p. 15. Scores measure agreement with these labels, not with ground-truth water.
• Leave-extreme-out has a gap. Events within 30 cm of the highest stage are excluded from training, not just the top event §3.2 p. 24. The exclusion rule is one of the things that makes the scores hard to compare.
• Reproduction is blocked. The code is proprietary and the stage data belong to national agencies and cannot be distributed. Precipitation (NASA IMERG) and Sentinel data are public; the DEM source imagery is not named [code availability] [data availability] p. 30 §2.3 p. 15.
• Table S1 is not in the source PDF. The per-gauge maximal lead times and population figures refer to it tab. 1 §2.5. Those values are not used here.
• Operational costs are hidden. The hydraulic model needed heavy manual calibration per gauge §2.3c p. 13. The Thresholding model needs only stage history and satellite-derived maps, but still needs a labelled history of past floods per AOI §2.3a p. 10.

Lineage

• Builds on: LSTM rainfall-runoff modelling by Kratzert et al. (2018) and related regional LSTM work §2.2b p. 8 §4 p. 26–27; the state-handoff structure of Gauch et al. (2021) §2.2b p. 8; the CMAL uncertainty head of Klotz et al., arXiv:2012.14295 §2.2b p. 8.
• Compared against: the physics-based hydraulic model of Ben-Haim et al. (2019) §2.3c p. 13; US NWS stage forecasts in Zalenski et al. (2017), cited only for a rough benchmark §4 p. 26.
• Companion work: Nevo et al., arXiv:2012.00671, on ML flood forecasting at scale and reach, cited for the effort comparison §4 p. 27.
• Superseded by: not known from wiki content.