25. Site flood forecasts

A design storm starts from a model's initial conditions: dry soil, an empty pond, the water table where [GROUNDWATER] puts it. A forecast for a real site cannot. The same two inches of rain do very different things to soil that has had a wet week, a water table a foot from the surface, and a pond already half full. aquaswmm.forecast is a Python subpackage that keeps a SWMM model running continuously so it carries that antecedent state from one day to the next, forces it with observed radar rain and forecast rain ensembles, corrects it with live observations, and reports the probability that site thresholds will be exceeded.

It is a first version. It runs on the unmodified EPA engine through the same aquaswmm.run the rest of the Python package uses (§18), and it is not yet in the desktop app. §25.11 lists what it does not do.

25.1 Install

pip install "aquaswmm[forecast]"

The forecast extra adds numpy (the filter and the grid extraction), ECMWF's eccodes (GRIB2 decoding; the Windows wheel carries the C library, elsewhere it comes from the eccodeslib wheel) and, before Python 3.11, tomli for the site file. import aquaswmm and import aquaswmm.forecast need none of them: hotstart files, continuous runs and forecast products work with the standard library, and the parts that need an extra say which one when they are used without it. An EPA SWMM 5.2 engine must be installed (§9).

The package's compiled module is now aquaswmm._native; aquaswmm re-exports everything in it, so existing scripts are unchanged.

25.2 Hotstart files

State passes from run to run through SWMM's own hotstart files: a run with SAVE HOTSTART "end.hsf" in [FILES] writes its final state, and a run with USE HOTSTART "end.hsf" starts from it. aquaswmm.forecast.hotstart reads, edits and writes them.

from aquaswmm.forecast import read_hotstart

h = read_hotstart("end.hsf", "site.inp")      # the model that wrote it
h.variables("groundwater", "S1")    # {'theta': 0.33, 'water_table': 4.22, 'flow': ..., 'max_infil_vol': ...}
h.get("node", "Pond", "depth", model_units=True)
h.set("node", "Pond", "depth", 3.0, model_units=True)
h.storage_volume("Pond")            # ft3 (m3 for SI) from the [STORAGE] shape
h.write("edited.hsf")

The layout is EPA's src/solver/hotstart.c (format 4, written by SWMM 5.2; format 3 from 5.1 is also read):

Kinds for get/set are subcatchment, infiltration, groundwater, snow, quality, node and link.

Nothing in the file names an object or says which subcatchments have groundwater or which nodes are storage units, so a file is always read against its model. The object order is the order objects first appear in the .inp, which is how SWMM indexes them; the tests check it against the names in the engine's .out. The reader refuses a file that does not fit:

• a stamp that is not SWMM5-HOTSTART, or formats 1 and 2;
• counts or flow units that differ from the model (the message lists each);
• a byte length that differs from what the model's layout needs, which is what happens when the counts match but a subcatchment gained or lost groundwater, a node became a storage unit, or the infiltration method changed.

Values in the file are SWMM's internal units, feet and cubic feet per second, whatever the model's flow units; model_units=True converts with the engine's own factors. Water table elevations are elevations, not depths.

An edit patches the bytes of that one value. A file read and written without an edit is the engine's file byte for byte, and rebuild() re-encodes every value from scratch and must equal the file, which the tests check on files EPA SWMM 5.2.4 wrote from the Groundwater, Detention Pond and Site Drainage samples. An edited file is accepted by the engine: raising the Detention Pond's empty storage unit SU1 to 6.00 ft and running on gives 5.916 ft at the first report step, draining from there, and a stored volume of 115,848 ft³, which is what the pyramidal shape gives for that depth.

One EPA quirk is reported rather than hidden: for each land use, 5.2's writer saves pollutants × pollutants buildup values but its reader reads pollutants, so the engine misreads hotstart files from models with more than one pollutant and any land use. The reader follows the writer, and Hotstart.warnings says so for such a model.

The 2D surface's groundwater. A 2D run with groundwater on (§20.9f) keeps its own state beside the hotstart: <model>.2d.gwstate, a small text file with each store's upper-zone moisture, water-table elevation and flows (20b §5.3). The CLI carries it the way SWMM carries a hotstart:

aquaswmm-swmm twod run site.inp --gw-state-in prev.2d.gwstate --gw-state-out next.2d.gwstate

and aquaswmm.forecast.gwstate reads and writes it. Moisture theta and the water-table elevation are the same quantities as a hotstart's groundwater record for the same subcatchment, so a chain can move them either way:

from aquaswmm.forecast import gwstate, read_hotstart

h = read_hotstart("end.hsf", "site.inp")
s = gwstate.from_hotstart(h, time_s=0)    # SWMM's aquifers -> the 2D stores
gwstate.write("start.2d.gwstate", s)

s = gwstate.read("site.2d.gwstate")       # a 2D run's end state ...
gwstate.apply_to_hotstart(s, h)           # ... into the next SWMM hotstart
h.write("next.hsf")

The flows do not carry across (SWMM stores its lateral flow per unit area, the 2D state a total); a state made from a hotstart starts them at zero, which affects only the first groundwater step.

25.3 Continuous runs

import datetime as dt
from aquaswmm.forecast import spinup, advance

s = spinup("site.inp", dt.datetime(2026, 6, 1), dt.datetime(2026, 9, 1), rain, runs_dir="runs")
s = advance(s, dt.timedelta(hours=6), rain, runs_dir="runs")
s.hotstart, s.time, s.record["continuity"]

spinup runs from the model's own initial conditions (or from warm= a state); advance continues a state for a period or to an end time. Each run sets START_DATE/END_DATE, the hotstart lines and the forced gages in a copy of the model and runs it in its own folder:

Chaining two runs gives the same state as one run over both periods: the groundwater moisture and water table after 48 h + 48 h equal those after one 96 h run to the last digit in the tests.

A run cannot be chained from if it failed, its results are not valid, it wrote no hotstart, or a quantity continuity error (runoff, groundwater, flow routing) is above continuity_limit (2 % by default) and that percentage is at least 0.01 acre-feet (hectare-metres for SI) of the largest volume in its table. The second condition matters: a dry six-hour run that drains a few cubic feet out of a pond reports 3.4 % flow routing continuity error on 0.002 acre-feet. Quality continuity is not checked.

Running the same inputs into the same folder reuses the finished run instead of running it again. A hotstart carries no clock, so a warm run must start exactly when its state ends; anything else is refused.

25.4 Rain

Each SWMM rain gage named in the site file is forced with hourly rain in a user-format rain file (FC_<gage> YYYY MM DD HH mm value), as a VOLUME gage with a 1:00 interval, in inches for US flow units and millimetres for SI. The FC_ prefix keeps a numeric gage name from looking like an NWS station to SWMM's file sniffing. Gages not in the site file keep their own data.

Observed rain is NOAA MRMS one-hour quantitative precipitation estimates on the 0.01° CONUS grid, millimetres, one gzipped GRIB2 file per hour stamped at the end of the hour it covers:

https://noaa-mrms-pds.s3.amazonaws.com/CONUS/{product}/{YYYYMMDD}/MRMS_{product}_{YYYYMMDD}-{HH}0000.grib2.gz

product is MultiSensor_QPE_01H_Pass2_00.00 (gauge-corrected, about an hour behind; the default) or RadarOnly_QPE_01H_00.00. The AWS archive starts 2020-10-14; earlier hours come from the Iowa Environmental Mesonet, https://mtarchive.geol.iastate.edu/{YYYY}/{MM}/{DD}/mrms/ncep/{name}/{name}_00.00_{YYYYMMDD}-{HH}0000.grib2.gz (GaugeCorr_QPE_01H before October 2020). Cells below zero are no coverage. An hour no archive has is missing, not zero: the chain stops before it and picks up when it appears, unless allow_missing_observed = true.

Forecast rain is:

• HRRR, 3 km, hourly cycles to 18 h (48 h at 00, 06, 12 and 18 UTC). Only the one-hour precipitation record (APCP:surface:{F-1}-{F} hour acc fcst) is fetched, by byte range from the .idx next to each file: https://noaa-hrrr-bdp-pds.s3.amazonaws.com/hrrr.{YYYYMMDD}/conus/hrrr.t{HH}z.wrfsfcf{FF}.grib2. It is one deterministic member.
• GEFS 0.25°, the control c00 and 30 perturbed members, 3-hourly from cycles at 00, 06, 12, 18 UTC: https://noaa-gefs-pds.s3.amazonaws.com/gefs.{YYYYMMDD}/{HH}/atmos/pgrb2sp25/ge{member}.t{HH}z.pgrb2s.0p25.f{FFF} (that layout from 2020-09-23). Precipitation accumulates 0-3 h then 0-6 h from each 6-hour boundary; the 3-6 h amount is the difference, and each 3-hour amount is spread evenly over its three hours.

The pipeline takes the newest cycle whose last needed file exists (HRRR up to 6 h old, GEFS up to 24 h), and the forecast starts at the analysis time, so a member's first hours are the part of the cycle still to come. An HRRR member covers what is left of its cycle, which may be shorter than forecast_hours; its run is shorter and the thresholds are judged on what it covers.

At the gages. A gage reads a point (the nearest cell), a point with radius_km (the mean over that disc) or a polygon (the area mean). A polygon is sampled on a lattice of points finer than the grid, each takes its nearest cell, and the mean over samples is the cell values weighted by how much of the polygon each covers. This works the same on the latitude-longitude MRMS and GEFS grids and HRRR's Lambert grid, using the grid's own point coordinates. [[areas]] gives a subcatchment its own polygon: the run adds a gage AREA_<subcatchment> and points the subcatchment at it.

Prefer areas to points. On 2026-09-30 the MRMS multi-sensor 24-hour total had one cell at 39.675 N 86.345 W reading 205 mm with every neighbour at zero, and the same cell held the largest 72-hour total in the region two days earlier. A point gage on it fed the model 8 inches of rain that did not fall, and the hourly values summed exactly to the 24-hour product, so no consistency check between MRMS products would have caught it. Nothing in this version removes such cells.

Times are UTC. The model runs on its own clock, UTC plus utc_offset_hours; keep that fixed (standard time or UTC), because a clock that jumps for daylight saving cannot carry state across the jump.

25.5 Observations and the ensemble Kalman filter

Observations are USGS instantaneous values, https://waterservices.usgs.gov/nwis/iv/?format=json&sites=...&parameterCd=...&startDT=...&endDT=... (public, no key): gage height 00065, discharge 00060, depth to water below land surface 72019, groundwater level 62610/62611. Values with USGS's no-data value or flagged equipment, ice, backwater, discontinued or maintenance are dropped; provisional data are used. Each observation maps to the hotstart value it measures: model = scale × reading + offset, in model units, taking the reading nearest the analysis time within tolerance_min. A stage gage over a pond is offset = gage datum − pond invert; a well reporting depth below land surface is scale = -1, offset = land surface elevation. USGS has said the legacy waterservices.usgs.gov endpoints will be retired for api.waterdata.usgs.gov; the base URL is a parameter of the reader.

The filter is a stochastic ensemble Kalman filter (Burgers et al., 1998) over an ensemble of hotstart files:

• State vector x: the hotstart values listed in [[assimilation.states]], in model units, e.g. groundwater:S1:water_table, groundwater:S1:theta, node:Pond:depth.
• Ensemble: after spin-up the hotstart is copied members times and each state gets Gaussian noise of its perturb_sd. In every cycle each member is run with the observed rain times its own lognormal factor (rain_sd_log, median 1), so the spread reflects rain error.
• Observation operator H: each observation names one hotstart value, which is what the instrument would read at the end of the cycle. It need not be a state: an observed link flow can correct the states it covaries with and is not written back.
• Update, with A and Y the state and predicted-observation anomalies and R the diagonal of error_sd²: K = A Yᵀ/(N−1) · (Y Yᵀ/(N−1) + R)⁻¹, and each member xᵢ ← xᵢ + K (y + eᵢ − H xᵢ), eᵢ ~ N(0, R).
• Bounds after the update: depths 0 … max depth, moisture wilting point … porosity, water table aquifer bottom … ground surface; or min/max per state.

No localisation, no inflation, no parameter estimation: a site model has a handful of states and the ensembles are small enough to read member by member. Each analysis writes analysis/<time>/analysis.json with the observations used, the innovations and the prior and posterior mean and spread.

Twin experiment. aquaswmm.forecast.ensemble.twin_experiment checks the filter on a model against a synthetic truth. On the test model (EPA's Groundwater_Model with a pond added, python/tests/data/groundwater_pond.inp): the truth starts with the water table at 4.0 ft and is forced by 179 mm of synthetic storms over 10 days; 16 members start at 2.5 ft (±0.5) and moisture 0.30 (±0.03), forced by the truth's rain × 0.8 × a lognormal factor (σ 0.4) per member and 6-hour cycle; the well (σ 0.05 ft) and the pond stage (σ 0.03 ft) are observed every 6 hours. Root-mean-square error of the ensemble mean against the truth, before each update, over 40 cycles:

The first cycle runs before any observation, so the last column leaves it out. The water table is pulled to the truth at the first update and stays. The pond stage improves less: between updates it is driven by the rain, and the members' rain is still biased and noisy. A shorter version (72 h, 8 members) runs in the test suite and requires the filter's water-table error to be under half the open loop's.

25.6 Forecast products

Thresholds compare a quantity from each member's .out with a value in model units:

A member exceeds when the quantity goes strictly above the value at any reporting step in the forecast window. The probability of exceedance is the share of members that do, all weighted equally (HRRR counts as one member beside the GEFS members). For each threshold the products give the probability, the earliest and median first exceedance, and per member whether it exceeded, when first, and its peak:

• forecast.json: everything, with the run metadata;
• forecast.csv: one row per threshold and member;
• forecast.html: a plain table that says it is model output, not a warning.

They are written to forecasts/<analysis time>/ and copied to latest/. A member whose run fails, or breaks the continuity limit, is left out and named in members_left_out.

25.7 Hindcasts

hindcast replays past issue times from [hindcast] start to end every every_hours, in its own hindcast/ folder: it chains observed MRMS rain to each issue time, runs the forecast members with the HRRR and GEFS forecasts issued then (from the same archives), and runs the same state with the observed rain over the window. Each threshold is scored against that observed-rain run, and, when the threshold names an observation, against the USGS record mapped to model units: Brier score and hit, miss, false alarm and correct negative counts with "yes" at probability ≥ 0.5. The results are hindcast.json and hindcast.csv.

The score against the observed-rain run isolates the rain forecast: the model and its state are the same on both sides. The score against observations includes the model's own error. Hindcasts run without the filter.

25.8 The site file and the command line

python -m aquaswmm.forecast example > site.toml   # a complete, commented file
python -m aquaswmm.forecast spinup site.toml
python -m aquaswmm.forecast run site.toml          # advance + forecast
python -m aquaswmm.forecast run site.toml --no-forecast
python -m aquaswmm.forecast hindcast site.toml

--offline uses cached data only. The exit status is 0 for done (or nothing new), 1 when a run or the data failed, with the reason on stderr, and 2 for a bad configuration.

Paths are relative to the site file; times are UTC. The configuration is checked against the model when it loads: unknown gages, subcatchments, nodes or groundwater are reported by name.

Scheduling. run is meant to be called every hour (or every day) by cron or Task Scheduler. It advances the analysis from state.json to the latest hour MRMS should have (observed_latency_hours behind now) in cycle_hours steps, then forecasts. Each step writes its results before the state moves, and state.json is replaced atomically, so an interrupted run is finished by the next one; a run with nothing new reuses the finished runs and cached rain and fetches only availability checks. Observed rain is cached as values at the gages (rain/mrms/<product>/<hour>.json, a few hundred bytes an hour), not as GRIB files; keep_raw = true keeps the files too. A spin-up of a month downloads about 720 MRMS files (roughly 0.6 GB) once. Change the gages and the cached values are fetched again; change the model and the state must be spun up again (delete state.json), which this version does not detect.

25.9 Worked example

The test model python/tests/data/groundwater_pond.inp is EPA's Groundwater_Model (docs/datasets/Groundwater_Model.inp: one 5-acre pervious subcatchment over a two-zone aquifer, ground surface 6 ft, water table 4 ft) with a pond added: subcatchment 1 and its groundwater drain to storage unit Pond (pyramidal, 60 × 40 ft base, 3:1 sides, 6 ft deep), which empties through a 4-inch orifice and, above 2.5 ft, over a 10 ft spillway weir. The location is made up: the model is placed at 33.88 N 97.38 W in north Texas only because rain fell there the week this was written.

[site]
name = "Groundwater pond, made-up location in north Texas"
model = "groundwater_pond.inp"
work_dir = "work"

[[gages]]
name = "1"
lat = 33.88
lon = -97.38
radius_km = 2.0

[[areas]]
subcatchment = "1"
polygon = [[33.87, -97.39], [33.89, -97.39], [33.89, -97.37], [33.87, -97.37]]

[rain]
forecast = ["hrrr", "gefs"]
forecast_hours = 24

[spinup]
start = "2026-09-24T00:00"

[assimilation]
members = 8
[[assimilation.states]]
key = "groundwater:1:water_table"
perturb_sd = 0.3
[[assimilation.states]]
key = "node:Pond:depth"
perturb_sd = 0.05

[[thresholds]]
name = "Pond over 1 ft"
kind = "node_depth"
object = "Pond"
value = 1.0

[[thresholds]]
name = "Spillway flowing"
kind = "link_flow"
object = "Spill"
value = 0.0

[[thresholds]]
name = "Water table above 4.5 ft"
kind = "groundwater"
object = "1"
value = 4.5

Run on 2026-10-01 at 16:18 UTC against the live archives:

• spinup chained 8 runs from a cold start at 2026-09-24 00 UTC to 2026-10-01 14 UTC (the latest hour two hours behind the clock) on MRMS Pass 2 rain: 0.05 in on 2026-09-27, 2.17 in on 09-28 and 6.11 in in the first 14 hours of 10-01, 182 hourly files. The first attempt stopped on a network timeout after four days; the next call resumed from the last finished chunk and finished in 1 min 40 s. It then seeded 8 members.
• run found the analysis already at 14 UTC, took HRRR's 14 UTC cycle (18 hours left of it; 5.4 mm at the site) and GEFS's 12 UTC cycle (31 members), and ran 32 members in about two minutes: Pond over 1 ft 100 % (32 of 32; the pond was already above it), Spillway flowing 100 %, Water table above 4.5 ft 12 % (4 of 32). Called again 40 minutes later, with no new MRMS hour yet, it reused everything and took 2 seconds.
• hindcast for issue times 2026-09-29 12 UTC, 09-30 00 UTC and 09-30 12 UTC (24-hour windows) spun up its own chain and replayed the forecasts issued then. For the first two, HRRR and GEFS had 0-20 mm and the model with observed rain stayed below every threshold: all correct negatives. For 09-30 12 UTC, HRRR had 70.5 mm and GEFS 27-103 mm (median 52.5) while MRMS measured 136 mm; the forecast gave Pond over 1 ft 38 % and Spillway flowing 13 %, and with observed rain both were exceeded, so both count as misses at the 0.5 cut. Brier scores over the three issues: 0.130, 0.255 and 0.000 (water table). Three issues at one place say nothing about skill; they show the replay and the scoring working on real archives.

There is no observation in this example: a USGS gage cannot measure a pond that does not exist, and assimilating one would only show that numbers move. The filter is demonstrated by the twin experiment (§25.5) and the pipeline test, which serves a synthetic well reading through the same code.

25.10 From Python

from aquaswmm.forecast import config, pipeline

cfg = config.load("site.toml")
pl = pipeline.Pipeline(cfg)
pl.spinup()
result = pl.forecast()              # advance + forecast
result["results"][0].probability

The pieces work on their own too: hotstart for any SWMM 5.2 model, chain.spinup/advance with your own Rain, sources.MRMS, HRRRSource, GEFSSource with a list of grids.Target, observations.USGS, enkf.EnKF over any list of hotstarts, and products.evaluate over any .out files.

25.11 Not yet

• Not in the desktop app; Python and the command line only.
• No quality control of MRMS beyond dropping no-coverage cells: single-cell artifacts (§25.4) go into the model unless an area target dilutes them.
• No bias correction or downscaling of HRRR or GEFS; GEFS 3-hour amounts are spread evenly over their hours. Members are weighted equally.
• HRRR is one member and covers only what is left of its cycle; no time-lagged HRRR ensemble.
• The filter has no localisation, inflation or parameter estimation, and does not update infiltration state, link flows or snowpack. Rain error is one lognormal factor per member per cycle.
• Hindcasts run without the filter, and need GEFS from 2020-09-23 on (the atmos/pgrb2sp25 layout).
• A changed model or gage set is not detected against an existing state.json.
• Snowmelt forcing (temperature) is not fetched; a model with snowpacks keeps its own climate data.
• The 2D surface (§20) is not run by the chain itself; its groundwater state can be carried beside the hotstart files (section 25.2), but the pipeline does not yet launch 2D runs or filter 2D states.
• USGS only; no other observation networks.