Introduction

Hydra is a water infrastructure simulation platform written in Rust. It is built as a suite of domain engines sharing one toolchain: a desktop GUI, a hydra CLI, and a Rust SDK.

See Engines for what each engine covers and what “planned” means in practice.

Water Distribution Engine

The engine that ships today. It performs extended-period simulation (EPS) of hydraulic behaviour and water quality dynamics across pressurised pipe networks, computing the full time history of flows, pressures, and constituent concentrations at every node and link.

Unless a page says otherwise, the rest of this documentation describes this engine.

Hydraulics

• Head-loss formulas: Hazen-Williams, Darcy-Weisbach, Chezy-Manning (with minor losses)
• Demand models: Demand-Driven Analysis (DDA) and Pressure-Dependent Analysis (PDA)
• Emitters: pressure-dependent outflow at junctions
• Leakage: FAVAD (Fixed and Variable Area Discharge) model
• Pumps: head-curve (1/3-point and custom), constant-power, variable-speed patterns
• Valves: PRV, PSV, FCV, TCV, GPV, PBV, PCV
• Tanks: cylindrical and volume-curve geometries, overflow mode
• Controls: simple time/level/pressure controls, rule-based controls with priorities
• Solver: Global Gradient Algorithm (GGA) with sparse Cholesky factorisation

Water Quality

• Modes: chemical constituent, water age, source trace
• Transport: Lagrangian segment-based advection
• Reactions: first-order and zero-order bulk and wall decay, limiting potential, roughness correlation
• Sources: concentration, mass booster, flow-paced booster, setpoint booster
• Tank mixing: complete (CSTR), two-compartment, FIFO (plug flow), LIFO

I/O

• Input: EPANET .inp format, any 2.x release (local files, HTTP URLs)
• Output: EPANET-compatible .out binary format, .rpt text report, .json report
• Unit systems: all 11 EPANET flow unit variants (CFS, GPM, MGD, IMGD, AFD, LPS, LPM, MLD, CMH, CMD, CMS)

Relationship to EPANET

The water distribution engine’s hydraulic and quality solvers were derived by studying EPANET’s mathematical foundations. Hydra is not an EPANET clone or compatibility layer; it is a distinct solver that models the same physics. Where the two diverge, Hydra’s result is authoritative.

The same principle will apply to each engine Hydra adds: it reads the established source-model format for its domain, and models the physics independently rather than reimplementing the reference tool.

For migration guidance, see Migrating from EPANET.

See INP Format Support for current EPANET input coverage.

Engines

Hydra is a platform, not a single solver. Each modelling domain gets its own engine — its own data model, source-model format, and numerical methods — and they all sit behind one shared toolchain: the desktop GUI, the hydra CLI, and the hydra-sdk Rust library.

Available vs. planned

Every engine — shipping or not — is registered in Hydra’s engine registry with an immutable descriptor: a stable key, a display label, a badge, an accent colour, a one-line summary, and the file formats it imports. The registry is the single source of truth; applications resolve a project’s stored engine key against it rather than hardcoding names or file filters.

The descriptor also carries a status:

• Available — implemented in this distribution and usable.
• Planned — registered and reserved, with no implementation behind it. Applications present it so the full modelling scope is visible, but refuse to create projects, import models, or run simulations for it.

“Planned” and “unknown” are deliberately distinct states. A planned key resolves successfully; only a key this build has never heard of — for example a project created by a newer Hydra — is an error, and it is always surfaced as an explicit unsupported state, never silently substituted with a default engine.

The crates for the planned engines (hydra-engine-uds, hydra-engine-och) are already published as empty scaffolds, so their names and version history track the workspace from the start rather than being introduced mid-life.

Why the split matters

Two engines can claim the same file extension with wholly incompatible contents — wds and uds both read .inp, one an EPANET model and one a SWMM model. An extension is therefore a file-picker filter, never a validity test: only the owning engine’s parser can decide whether a file really is one of its models. Hydra rejects a mismatched model explicitly rather than loading a stormwater network as a plausible-looking pressurised one — see Foreign .inp dialects for what that rejection looks like in the CLI, the GUI, and the SDK.

The same separation runs through the codebase. hydra-common holds only what every engine shares — engine identity and the reportable-output contract — and carries no engine vocabulary. Engines emit neutral content fragments; hydra-report renders them without knowing which engine produced them. See Crate Layout for how the workspace is arranged.

Working with the registry from code

#![allow(unused)]
fn main() {
use hydra_sdk::common::{engine_by_key, ENGINES};

for engine in ENGINES {
    println!("{} ({}) — {:?}", engine.label, engine.key, engine.status);
}

let wds = engine_by_key("wds")?;
assert!(wds.is_available());
}

See the SDK Overview for the full engine-identity API.

Installation

GUI — Desktop Application

Download the installer for your platform from the releases page:

After installing, see Troubleshooting if macOS blocks the app from opening.

CLI — Command Line

For most users, Cargo install is the recommended path.

Option 1 — Install with Cargo (recommended)

cargo install hydra-cli

Requires Rust ≥ 1.95 (install via rustup.rs).

After installing, verify with:

hydra -V

Option 2 — Pre-built binary (no Rust required)

Download the hydra binary for your platform from the releases page and place it somewhere on your PATH.

macOS — Pre-built CLI binaries are currently not notarised. If Gatekeeper blocks the binary, remove the quarantine flag:

xattr -d com.apple.quarantine hydra

Building from Source

If you want to build Hydra yourself (e.g. to contribute or run the test suite):

Prerequisites

• Rust stable ≥ 1.95 — rustup.rs
• just — cargo install just or brew install just
• GUI only: Node.js 24, pnpm 11, Tauri CLI (cargo install tauri-cli), and the Tauri system prerequisites for your platform

git clone https://github.com/neeraip/hydra
cd hydra
just setup          # optional: install Cargo deps, frontend deps, and CLI tools (needs pnpm)
just build          # debug build
just release        # optimised release build (fat LTO, embeds the GUI frontend — needs Node/pnpm)
just test           # run the full test suite

Key Concepts

This page explains the terminology of the water distribution engine — the vocabulary used throughout Hydra’s documentation and in EPANET .inp files. If you are familiar with EPANET, most of these will be review.

Each engine brings its own vocabulary; see Engines for the full roster.

Network Elements

A water distribution network in Hydra is made up of nodes (points) and links (connections between nodes).

Nodes

Links

Hydraulic Concepts

Time and Patterns

Water Quality

File Formats

GUI

Hydra’s desktop application lets you load, run, and explore simulations without using the command line. It is the front end to every Hydra engine — today that means the water distribution engine; planned engines appear in the app but cannot yet be used.

Download and Install

Download the installer for your platform from the releases page.

macOS — “Hydra is damaged and can’t be opened”

Hydra GUI macOS releases are notarised. If Gatekeeper still shows this warning after dragging Hydra to Applications, run this once in Terminal:

xattr -cr /Applications/Hydra.app

Then open the app normally.

Basic Workflow

Hydra organises work into projects. Each project holds a network model and one or more scenarios — independent parameter sets you can run and compare.

1. Create a project — on the Projects screen, click New Project. The wizard asks for the engine first, then project details, then a review. Choosing the engine up front is deliberate: .inp belongs to both EPANET and SWMM, so the file extension cannot decide the modelling domain on your behalf. Engines marked as planned are shown but cannot be selected. With Water Distribution chosen, either import an existing EPANET .inp file or start from a blank network.
2. Configure and run — press ⌘R (macOS) or Ctrl+R (Windows/Linux), or click the Simulate button in the scenario strip at the bottom of the screen. Select which scenarios to run and confirm.
3. Explore results — after the simulation completes, the network map updates with colour-coded results. Click any node or link to inspect its time-series values (pressure, head, flow, velocity, water age, etc.). Use the timeline scrubber to step through reporting periods.

Press ⌘K (macOS) or Ctrl+K (Windows/Linux) at any time to open the command palette, which lists every action — running simulations, switching views, imports and exports — filtered as you type by substring match. Type # in the palette to find any node or link by ID and zoom to it.

Keyboard Shortcuts

Beyond ⌘R (simulate), ⌘K (command palette), and ⌘Z / ⇧⌘Z (undo/redo), the app has shortcuts for navigation and the canvas. Press ? at any time to open the in-app cheatsheet, or see the full Keyboard Shortcuts reference. Common ones (⌘ on macOS, Ctrl elsewhere):

On the canvas, single-key tools select elements and add or measure geometry (select, edit, add node, add link, measure), and annotations can be placed on the map.

Editing the Network

The Network tab provides editable tables for junctions, reservoirs, tanks, pipes, pumps, and valves, including each pipe’s initial status. Links are drawn as polylines, and their intermediate vertices are preserved and rendered on the canvas; dragging a node moves the link’s endpoint while the intermediate vertices stay fixed. Committed edits can be undone and redone with ⌘Z / ⇧⌘Z (Ctrl+Z / Ctrl+Shift+Z).

The Issues panel collects network validation findings (structural problems detected before a run) and warnings produced by the last simulation run, with links to the affected elements.

Dedicated editors are available for curves, patterns, and controls.

If a model’s coordinates use a projected coordinate system, the CRS picker on the canvas can scan the network’s coordinates and suggest matching coordinate reference systems so the network lines up with the basemap. You can also define and save custom CRS entries.

Scenarios and Comparison

Scenarios let you keep independent parameter sets side by side within one project. On the canvas, a comparison overlay can display the delta between the active scenario’s results and a baseline (the base model or another scenario).

The Analysis tab includes a system summary (key metric chips), result histograms, pipe criticality, pump energy, audit panels, and tank level charts.

Units

Choose between SI (metric) and US customary display units in Settings. This affects how values are shown and entered throughout the app; files and exports (INP, CSV, GeoJSON) always remain in the model’s native units. Settings also offers a light / dark / system theme.

Performance on Large Networks

Hydra GUI is tuned to stay responsive on larger models.

• Opening a project navigates immediately while network data finishes loading.
• Network Inspector node/link lists use virtualized rendering to avoid large DOM slowdowns.
• Basemap switching keeps network overlays attached so features remain visible while the style reloads.

Exporting and Output Files

Hydra saves simulation results inside the project folder on disk. To open the folder for a scenario, go to the Scenarios panel and click the Open in Finder icon next to the scenario name; it reveals the folder in your platform’s file manager (Finder, Explorer, or the Linux equivalent). The folder contains results.out — EPANET-compatible binary output, readable by post-processing tools.

Other formats are available from the command palette (⌘K / Ctrl+K):

• Export INP… — save the current network as an EPANET .inp file
• Export results as CSV… — save node and link result time series as CSV files (shown once results exist)
• Export results to GeoJSON — save nodes/links with attributes, including result values when available

For a plain-text .rpt report, run the exported .inp through the CLI.

Supported Networks

Any EPANET .inp file works directly — any 2.x release, no conversion needed. See INP Format Support for the full coverage list.

Troubleshooting

See the Troubleshooting page for common issues including the macOS Gatekeeper error and Windows Defender prompts.

Keyboard Shortcuts

A complete reference for the Hydra desktop GUI. Press ? at any time to open the in-app cheatsheet.

Throughout this page, the primary modifier is ⌘ (Command) on macOS and Ctrl on Windows/Linux; ⇧ is Shift. Single-key shortcuts (like the canvas tools) are ignored while you are typing in a text field.

Global

Project navigation

These require a project to be open.

On the Projects screen, ⌘F / Ctrl+F focuses the projects search box.

Simulation

⌘R opens a confirmation dialog where you select scenarios; the run itself starts on ⌘Enter (or the confirm button).

Canvas and map

These work on the Canvas view. The zoom/layout shortcuts also switch you to the Canvas view first.

Canvas tools (single keys)

Active on the Canvas view; the edit/add/measure tools apply in geographic (map) mode.

Canvas mouse

Timeline playback (single keys)

The same keys work when the timeline scrubber has focus.

Command palette

While the palette is open (⌘K / Ctrl+K):

Editors and dialogs

CLI

The hydra binary runs models through Hydra’s engines and builds report documents from the results.

hydra run <MODEL> [--engine KEY] [--results PATH] [--summary PATH]
hydra report --model <PATH> --results <PATH> [-o PATH]
hydra engines

The engine is detected from the model, never from its extension — .inp belongs to both EPANET and SWMM, so the filename cannot decide. There is no default engine: if the model does not identify one, Hydra stops and asks you to name it with --engine rather than guessing. See Engine selection.

Only the water distribution engine is implemented today; hydra engines lists what this build provides.

Upgrading from 2.x — hydra <model> <report> <output> is gone. Use hydra run <model> --summary <report> --results <output>. Running the old form prints a hint naming the replacement. See Migrating from EPANET for the full mapping.

Install

For most users, Cargo install is the recommended path on macOS, Linux, and Windows.

Option 1 — Pre-built binary (no Rust required)

Download the hydra binary for your platform from the releases page and place it on your PATH.

macOS — Pre-built CLI binaries are currently not notarised. If Gatekeeper blocks the binary, remove the quarantine flag:

xattr -d com.apple.quarantine hydra

Option 2 — Cargo (recommended)

cargo install hydra-cli

Verify the installation:

hydra -V

Basic Usage

# Run a simulation — summary goes to stdout
hydra run network.inp

# Save the summary to a file
hydra run network.inp --summary report.rpt

# Save the summary and the binary time-series results
hydra run network.inp --summary report.rpt --results output.out

hydra run takes exactly one positional argument: the model. Everything the run writes is named by a flag, so nothing depends on argument order. hydra -V prints the Hydra engine version and the CLI version on separate lines.

Output Formats

The --summary path controls the summary format:

# Plain-text summary (EPANET-style .rpt)
hydra run network.inp --summary report.rpt

# JSON summary (useful for scripts and data pipelines)
hydra run network.inp --summary report.json

# Binary results (.out) — EPANET-compatible, readable by post-processing tools
hydra run network.inp --summary report.rpt --results output.out

For configurable report documents (txt, csv, html, pdf) built from a saved .out, see Generating a report.

Running from a URL

The model may be fetched over HTTP or HTTPS:

hydra run https://example.com/network.inp
hydra run https://example.com/network.inp --summary report.rpt --results output.out

Both http:// and https:// are accepted. The fetch follows up to 10 redirects, uses a 10-second connect timeout and a 300-second overall timeout, and accepts response bodies up to 1 GiB. An HTTP 4xx response is treated as an input error (exit 1); a 5xx or network failure is an I/O error (exit 3).

Flags

hydra run

Global

Global flags may appear before or after the subcommand.

Engine selection

A model’s engine is decided by its contents, not its filename. Each engine is asked whether the model is one of its own, and exactly one positive identification is required to proceed.

There is deliberately no fallback. Handing a stormwater model to a pressurised-pipe solver would produce a confident, wrong answer rather than a failure, so Hydra stops instead of guessing.

--engine <KEY> names the engine explicitly. That is more information than detection has, so it is also the escape hatch for a sparse model that carries nothing identifying — the named engine parses it under its normal rules.

hydra run net.inp --engine wds     # skip detection
hydra engines                      # what this build provides

Generating a report

hydra report builds a report document from a completed run’s results. It is a separate step from the simulation: you point it at the model and the .out file the run produced.

hydra report --model network.inp --results output.out -o report.html

The content comes from the engine’s report blocks — named, self-contained sections such as run summary, result extremes, pump energy, service compliance, and the distribution charts. A template selects and orders them; without one you get everything that applies to the run. See Post-Simulation Analytics for what the blocks cover.

--no-timestamp exists for diffing and for reproducible builds: with it, the same inputs produce byte-identical output.

Exit Codes

Breaking change — internal errors previously exited with code 2 (the solver-error code). They now exit with the dedicated code 4; codes 0–3 are unchanged.

Reading the Report

Both report formats are summary-level. Per-node and per-link time series are written only to the binary .out file (--results).

The text report (.rpt) contains:

• Header — a Hydra version banner and the network title
• Input summary — element counts, head-loss formula, demand model, timesteps, and simulation duration
• Warnings — non-fatal diagnostics raised during the run: unbalanced hydraulics, negative pressures, and pump-head warnings
• Analysis timestamps — “Analysis begun” / “Analysis ended” markers

It does not contain per-node/link result tables, a network-status section, or an energy-usage section. Use the .out file for full results, or the JSON report’s energy block for the energy summary.

The JSON report contains the same summary-level data plus energy, flow-balance, and mass-balance blocks, in a structured format:

{
  "input": {
    "junctions": 92, "reservoirs": 1, "tanks": 2,
    "pipes": 117, "pumps": 2, "valves": 0,
    "headloss_formula": "Hazen-Williams", "demand_model": "DDA",
    "hydraulic_timestep_s": 3600.0, "quality_timestep_s": 360.0,
    "duration_s": 86400.0, "report_timestep_s": 3600.0
  },
  "warnings": [...],
  "energy": { "pumps": [...], "peak_demand_kw": 12.3 },
  "flow_balance": { ... },
  "mass_balance": { ... },
  "analysis": { "begun_epoch": "1615687166", "ended_epoch": "1615687167" }
}

The begun_epoch / ended_epoch values are strings holding raw seconds since the Unix epoch (or null if unavailable), not formatted datetimes.

For full time-series data across all nodes and links, use the binary .out format.

Diagnostics on stderr

Independently of the report, the simulation command emits warnings and errors to stderr as one JSON object per line, suitable for machine parsing in scripts and pipelines:

{"level":"warning","code":"warning/negative_pressure","message":"...","object_id":"J1","time_step":3600.0}
{"level":"error","code":"solver/hydraulic","message":"...","object_id":null,"time_step":null}

Each line has level (warning or error), a code, a human-readable message, and nullable object_id and time_step fields. The human-readable progress bar and banner (also on stderr) are suppressed by -q/--quiet and when stderr is not a terminal; the JSON diagnostics are not.

hydra report does not emit these. It shares the same exit codes but reports failures as plain error: … lines on stderr.

Troubleshooting

GUI

macOS — “Hydra is damaged and can’t be opened”

Hydra GUI macOS releases are code-signed and notarised. If Gatekeeper still shows this message, the app bundle usually has stale quarantine metadata from the download/copy step.

To open Hydra after installing it to /Applications, run this once in Terminal:

xattr -cr /Applications/Hydra.app

Then open the app normally from Finder or Spotlight.

macOS — App opens but immediately quits

This can happen if the app was launched directly from the .dmg or if macOS retained stale quarantine metadata. Move Hydra.app to /Applications first, then run the xattr -cr command above only if the app still fails to open.

Windows — “Windows protected your PC” (SmartScreen)

Click More info, then Run anyway. SmartScreen warns on unsigned executables. This will be resolved once Hydra’s Windows builds are code-signed.

Linux — AppImage does not open

Make the AppImage executable before running it:

chmod +x Hydra-*.AppImage
./Hydra-*.AppImage

If you see a FUSE-related error, install the required library:

# Ubuntu / Debian
sudo apt install libfuse2

# Fedora
sudo dnf install fuse-libs

“This is a SWMM model, not a water-distribution one”

The .inp extension is shared by EPANET and SWMM, so a file picker cannot tell them apart — only the parser can. Hydra detected a section the EPANET data model has no concept of and stopped rather than misreading the file.

Nothing is wrong with the file. It belongs to the urban drainage engine, which is planned but not yet implemented, so Hydra cannot open it today. See Foreign .inp dialects.

Canvas features disappear after changing basemap

This should not happen in current releases, but if the map style reload fails on a specific GPU/driver stack, try:

• Switch basemap once more (for example, to No basemap, then back).
• Change to another project tab and return to the canvas.
• Restart Hydra to reset the map renderer.

If the issue keeps reproducing, open an issue and include your OS version, GPU model, and whether it happens on all projects or only specific large networks.

CLI

macOS — “hydra cannot be opened because the developer cannot be verified”

Pre-built macOS CLI binaries are currently not notarised. If this warning appears for a downloaded binary, clear the quarantine attribute and try again:

xattr -d com.apple.quarantine hydra

Then move it to your PATH and run normally.

hydra: command not found

The hydra binary is not on your PATH.

• If you installed with cargo install hydra-cli, ensure ~/.cargo/bin is on your PATH:

  export PATH="$HOME/.cargo/bin:$PATH"
  

  Add this line to your shell profile (.bashrc, .zshrc, etc.) to make it permanent.

• If you downloaded a pre-built binary, move it to a directory that is already on your PATH (e.g. /usr/local/bin on macOS/Linux).

Exit code 1 — Input error

Hydra could not read or parse the network file. Common causes:

• The file path is wrong or the file does not exist.
• The .inp file contains a syntax error. Check the report for the specific line.
• A URL was provided but the server returned 4xx. Verify the URL is accessible.
• The file is a SWMM .inp, not an EPANET one — the diagnostic code is input/engine rather than input/parse. The file is fine; it belongs to an engine Hydra does not yet ship. See Foreign .inp dialects.

Exit code 2 — Solver did not converge

The hydraulic solver could not find a balanced solution for one or more time steps. This usually means the network model itself has an issue:

• Check for isolated nodes or disconnected sub-networks.
• Verify pump curves and valve settings are physically reasonable.
• Try setting UNBALANCED CONTINUE 10 in the [OPTIONS] section to let the simulation proceed past the failing step and produce a partial report for diagnosis.

Exit code 3 — I/O error

Hydra could not write output. Check that the output directory exists and that you have write permission.

Getting Help

If the steps above do not resolve your issue, open a GitHub issue with:

• The Hydra version (hydra -V)
• Your operating system and version
• A minimal .inp file that reproduces the problem (if applicable)
• The full error message or report output

→ Open an issue

INP Format Support

Hydra parses the EPANET .inp file format. This page documents which sections and keywords are supported, which are silently ignored, and where Hydra’s behaviour differs from or extends the standard.

Supported EPANET versions

Any EPANET 2.x file. The tables below are written against EPANET 2.3, which is the newest dialect Hydra understands — not a requirement your file has to meet. Everything 2.3 added over earlier releases is optional: the [LEAKAGE] section and the DISABLED suffix on a control line. A 2.0 or 2.2 file that uses neither loads unchanged and runs, with leakage simply zero on every pipe.

Older constructs are skipped rather than rejected. A legacy [ROUGHNESS] section — superseded by the roughness column in [PIPES] — is accepted as a no-op, as is any section or [OPTIONS] keyword Hydra does not recognise. Rejection is reserved for a file that is not an EPANET model at all; see Foreign .inp dialects.

Sections

Fully Supported

All data in these sections is parsed and applied to the simulation, with one exception noted below ([REPORT]).

Silently Ignored

These sections are recognised and accepted without error but produce no simulation effect. Files containing them parse cleanly.

Unknown sections (not listed in either table) are also silently ignored for forward compatibility.

An [END] marker, if present, terminates parsing: any content after the first [END] line is ignored.

Foreign .inp dialects

The .inp extension is not exclusive to EPANET — SWMM uses it too, for a wholly different data model. Because the water distribution engine ignores sections it does not recognise (above), a SWMM file would otherwise parse “successfully” into a network of junctions carrying each node’s maximum depth as its demand, joined by no links at all: a wrong answer wearing the costume of a right one.

So the parser rejects a foreign dialect up front, before any network is built. The test is positive only — it fires on the presence of a section EPANET has no concept of, never on the absence of one EPANET expects, because a valid EPANET model is not required to contain any particular section. The markers are a fixed list of SWMM-only section names ([SUBCATCHMENTS], [CONDUITS], [OUTFALLS], [RAINGAGES], [INFILTRATION], [POLLUTANTS], [XSECTIONS], and roughly two dozen more), matched on the upper-cased name.

This is reported as an engine mismatch, not a bad file — the same bytes may be a flawless model in the tool that owns them:

Once the urban drainage engine lands, such a file becomes openable rather than merely diagnosable. Until then it is a dead end: Hydra can tell you exactly what the file is, but has no engine to open it with.

OPTIONS Keywords

The [OPTIONS] keywords listed below are parsed and applied. A few EPANET keywords are not parsed — notably PRESSURE (pressure display units) and MAP — and any unknown keyword is silently ignored.

Pump Curves

A single-point pump curve (Q₁, H₁) is automatically expanded to a three-point power-function curve (0, 1.33334·H₁), (Q₁, H₁), (2·Q₁, 0), matching EPANET’s internal behaviour.

LEAKAGE Section

[LEAKAGE] was added in OWA-EPANET 2.3 and is not present in legacy EPANET 2.2 files. Each row specifies per-pipe FAVAD (Fixed and Variable Area Discharge) leakage coefficients:

[LEAKAGE]
;PipeID   C1       C2
P1        0.0002   0.5
P2        0.00015  0.6

Where C1 is the fixed-area discharge coefficient and C2 is the variable-area discharge coefficient. Standard EPANET files (without a [LEAKAGE] section) parse cleanly; leakage is simply zero for all pipes.

Differences from EPANET 2.3

Controls & Rules

Hydra supports two mechanisms for changing the network during a simulation: simple controls ([CONTROLS]) and rule-based controls ([RULES]). Both change a link’s status or setting in response to conditions; rules are more expressive and are evaluated more often.

Simple controls ([CONTROLS])

Each simple control is a single line that acts on one link. Every control begins with LINK. Keywords are case-insensitive.

Action is one of:

• OPEN or CLOSED (CLOSE is accepted as a synonym for CLOSED), or
• a numeric setting — a bare number sets a pump’s relative speed or a valve’s setting. For example, LINK PU1 1.0000 IF NODE T1 BELOW 4.0 runs pump PU1 at speed 1.0 while tank T1 is below 4.0.

When a numeric setting is given without a status, the status is inferred from the link type: a pump or pipe opens for a value greater than zero and closes at zero; a valve becomes active.

Time values accept H:MM or H:MM:SS, or a plain number. A bare number is interpreted as hours unless followed by a unit (SECONDS, MINUTES, HOURS, or DAYS). AT CLOCKTIME additionally accepts a trailing AM/PM.

Examples:

[CONTROLS]
LINK 1A OPEN IF NODE A BELOW 2.5275
LINK 1A CLOSED IF NODE A ABOVE 3.2689
LINK LINK-7491 OPEN AT TIME 18:00:00

Evaluation. Simple controls are checked once per hydraulic timestep. If several controls act on the same link in one step, the last one (in file order) wins.

A control line that does not begin with LINK or has too few tokens is silently ignored; an unknown node or link ID is a parse error.

Rule-based controls ([RULES])

A rule has the form:

RULE <label>
IF <premise>
AND|OR <premise>
...
THEN <action>
AND <action>
...
ELSE <action>
...
PRIORITY <value>

The THEN block runs when the premises are satisfied; the optional ELSE block runs otherwise. A rule must have at least one premise.

Premises

Each premise tests one object’s attribute against a value:

<object> <id> <attribute> <operator> <value>      # node or link
SYSTEM <attribute> <operator> <value>             # simulation clock

Operators: = (also ==, IS, EQUALS); <> (also !=, NOT); < (also BELOW); > (also ABOVE); <=; >=. Each operator is a single token — EPANET’s two-word IS NOT is not supported.

Combining premises

Premises are joined with AND and OR, and AND binds more tightly than OR — the condition is a disjunction of AND-clauses. For example:

IF  TANK T1 LEVEL BELOW 2
AND SYSTEM CLOCKTIME >= 22:00
OR  TANK T1 LEVEL BELOW 1

is read as (T1 level < 2 AND clocktime ≥ 22:00) OR (T1 level < 1).

Write rule clock times in 24-hour form (H:MM, e.g. 20:00). A trailing AM/PM is honoured in simple-control AT CLOCKTIME lines but is not applied inside rule premises.

Actions

THEN/ELSE actions set a link’s status or setting (rule actions cannot target nodes):

LINK <id> STATUS OPEN|CLOSED
LINK <id> SETTING <value>          # SPEED is a synonym for SETTING

The object keyword before the link ID is cosmetic — the ID is always resolved as a link. An optional IS and = are allowed for readability (THEN PUMP HSP#1 STATUS IS OPEN).

Priority and conflicts

PRIORITY <value> sets a rule’s priority (default 0). When rules that fire in the same step give conflicting instructions for one link, the higher priority wins; equal priorities are broken in favour of the earliest-defined rule. A STATUS action and a SETTING action on the same link do not conflict.

Evaluation

Rules are evaluated at a rule timestep that subdivides each hydraulic step (defaulting to one-tenth of the hydraulic step). When a rule changes the network mid-step, the hydraulics are re-solved for the remainder of the step, so a rule can react to a tank crossing a level partway through a step.

Example:

[RULES]
RULE 1
IF SYSTEM CLOCKTIME >= 6:00
AND SYSTEM CLOCKTIME < 20:00
AND TANK Tank LEVEL BELOW 97
THEN PUMP HSP#1 STATUS IS OPEN
AND PUMP HSP#2 STATUS IS OPEN
PRIORITY 1

See INP Format Support for how the [CONTROLS] and [RULES] sections fit into the wider file format, and Diagnostics & Errors for the validation errors these sections can produce.

Output Files

Hydra produces three output files. Only one of them holds full time-series results; the other two are summaries of the same run.

Only the .out file contains time series. The .rpt and .json reports are two serialisations of the same summary information — neither includes per-node or per-link result tables.

How each is produced:

• CLI — the report path (.rpt or .json) selects the text or JSON report; --output writes the .out file. See CLI.
• GUI — every run writes results.out into the scenario folder; CSV and GeoJSON exports are available from the command palette. See GUI.
• SDK — io::out_writer::write_binary_output, io::rpt_writer::build_text_report, and io::rpt_writer::build_json_report. See SDK Examples.

.out — Binary results

The .out file is the EPANET 2.3 binary output layout (format version 20012), so tools that read EPANET output files read Hydra’s. Hydra records its own metadata — such as the topology digest that detects a since-edited model — in a run.json beside the results rather than inside them, so the results file stays a format EPANET defines. Values are stored as 32-bit floats (REAL4) and 32-bit integers (INT4), little-endian; IDs and strings are fixed-width, zero-padded.

The file is written in five sections:

1. Prolog — a 15-integer header (magic number, format version, element counts, quality mode, trace node, flow-unit code, pressure-unit code, report start/step, duration) followed by the network’s static data: title lines, input/report filenames, chemical name and units, node and link IDs, link end-node indices and type codes, tank node indices and cross-section areas, and node elevations, link lengths, and link diameters.
2. Energy — per-pump summary: percent online, average efficiency, average energy per unit flow, average and peak power, and average cost, plus the network demand charge.
3. Dynamic results — one record per reporting period (not per hydraulic step), stored column-major. Each record holds, for that period:
   • Node quantities (4): demand, head, pressure, quality
   • Link quantities (8): flow, velocity, head loss, quality, status, setting, reaction rate, friction factor
4. Network reactions — average bulk, wall, tank, and source reaction rates over the run.
5. Epilog — the number of reporting periods, a warning flag, a network content digest, and a closing magic number.

Units

All numeric values are written in the unit system chosen at write time (output_units). The CLI and GUI default to the model’s declared flow units; the SDK writer takes an explicit FlowUnits argument (pass sim.net().options.flow_units for the model’s units). Pressures are written in metres (SI) or PSI (US customary).

Two categories are not unit-converted: tank cross-section areas (always internal m²) and all quality/reaction-rate values (written in their native quality units).

.rpt — Text report

A human-readable summary in EPANET report style. It contains, in order:

• A date stamp and the Hydra banner with the version number
• The network title lines
• An input summary: counts of junctions, reservoirs, tanks, pipes, pumps, and valves; head-loss formula; demand model (DDA/PDA); hydraulic timestep; hydraulic accuracy; maximum trials; the quality-analysis mode (with constituent name, quality timestep, and tolerance where applicable); specific gravity; demand multiplier; total duration; and report timestep
• An “Analysis begun” timestamp
• Warnings raised during the run, grouped by simulation time (unbalanced hydraulics, negative pressures, pump-exceeds-maximum-head)
• An “Analysis ended” timestamp

It does not contain per-node or per-link result tables. For full results, read the .out file. For the catalogue of warning types, see Diagnostics & Errors.

.json — JSON report

The same summary data as the .rpt report, in a structured form suited to scripts and pipelines. Top-level keys:

{
  "input": {
    "junctions": 92, "reservoirs": 1, "tanks": 2,
    "pipes": 117, "pumps": 2, "valves": 0,
    "headloss_formula": "Hazen-Williams", "demand_model": "DDA",
    "hydraulic_timestep_s": 3600.0, "quality_timestep_s": 360.0,
    "duration_s": 86400.0, "report_timestep_s": 3600.0
  },
  "warnings": [
    { "time": 3600.0, "code": "warning/negative_pressure", "message": "...", "object_id": "J1" }
  ],
  "energy": {
    "pumps": [
      { "pump_id": "PU1", "kwh": 120.0, "total_cost": 14.4,
        "avg_efficiency": 0.75, "max_kw": 8.1, "time_online_s": 86400.0 }
    ],
    "peak_demand_kw": 12.3
  },
  "flow_balance": {
    "total_inflow": 0.0, "total_outflow": 0.0,
    "tank_change": 0.0, "unaccounted": 0.0, "ratio": 1.0
  },
  "mass_balance": {
    "initial": 0.0, "added": 0.0, "demand": 0.0,
    "reacted": 0.0, "reacted_bulk": 0.0, "reacted_wall": 0.0,
    "reacted_tank": 0.0, "source": 0.0, "final_mass": 0.0, "ratio": 1.0
  },
  "analysis": { "begun_epoch": "1615687166", "ended_epoch": "1615687167" }
}

Notes:

• warnings[].time is the simulation time in seconds; object_id is the affected node/link ID or null.
• avg_efficiency is a fraction in [0, 1].
• flow_balance is null if unavailable, and mass_balance is null when no quality analysis was run. Balance volumes are in m³; mass values are in mg.
• begun_epoch / ended_epoch are strings holding raw seconds since the Unix epoch (or null) — not formatted datetimes.

Post-Simulation Analytics

After a run, Hydra can derive higher-level metrics from the saved results. There are three surfaces:

• Programmatic analytics (SDK) — two on-demand modules, demand reliability and service compliance, that read a saved .out file and return a structured report.
• Report blocks (SDK) — a catalog of named, self-contained content blocks that render into txt/csv/html/PDF documents. Several of them are built on the two modules above, which is how those metrics reach a generated report. See SDK Overview.
• The GUI Analysis tab — an interactive dashboard computed separately from the same results.

The CLI exposes report blocks through its report subcommand (CLI) but not the two modules directly. The GUI uses report blocks for its generated reports, while its Analysis tab computes its own dashboard rather than calling the two modules.

Report blocks

The water distribution engine publishes thirteen blocks. Each is self-contained: it carries its own heading and renders identically in every output format.

A block that does not apply to a run is not dropped: it renders as a placeholder section under its normal heading, carrying the reason. Asking for wds.pump-energy on a network with no pumps yields a Pump Energy section reading [not available: the network has no pumps]. This is deliberate — a requested section that silently vanished would be indistinguishable from one that was never requested.

Some blocks accept options — the *-thresholds pair takes its band edges, and several take a worstCount for their worst-performing tables. report_block_options(id, network) describes what a given block accepts, including labels, defaults, and bounds, so an interface can build an editor for it without hardcoding the list.

Templates

A template is the saved answer to “what goes in my report”: a document title plus an ordered list of block references. The GUI’s template builder and the CLI’s --template flag read the same JSON.

{
  "version": 1,
  "title": "Quarterly hydraulic report",
  "blocks": [
    { "id": "wds.run-summary" },
    { "id": "wds.pump-energy", "title": "Pumping cost" },
    { "id": "wds.pressure-thresholds", "options": { "edges": [20, 40, 60] } }
  ]
}

Unknown fields are ignored on read, so a template written by a newer Hydra still loads. An id that is not in the catalog is not an error: it renders as a placeholder section headed with the id itself, so a mistyped id is visible in the output rather than silently dropped.

Omit --template entirely to get every block that applies to the run.

Formats

The same document renders four ways, and the differences are presentational rather than editorial — every format carries the same sections in the same order.

Only the PDF is typeset, so only it paginates: a section heading is never left as the last thing on a page, and a table continuing onto another page repeats its column header. A long table still breaks where it falls — sections are not forced onto fresh pages.

Demand reliability (SDK)

Measures how well the network delivered the demand that was asked of it, per junction and network-wide.

Entry points (see SDK Examples):

#![allow(unused)]
fn main() {
compute_demand_reliability_from_out(out_path, &network) -> DemandReliabilityReport
compute_demand_reliability_from_out_with_options(out_path, &network, options)
}

It needs both the .out file and the loaded Network (the network supplies the demand categories and patterns used to recompute required demand). Only junction nodes are analysed.

What it computes, per reporting period and junction:

• Required demand (from the model) and delivered demand (from the .out file), each accumulated into a volume.
• Unmet volume (required − delivered, clamped at zero) and surplus volume (delivered − required, non-zero only under PDA when pressure exceeds the required-pressure threshold).
• Deficit periods — the number of periods where the shortfall exceeded the deficit tolerance — plus the longest consecutive deficit streak and the maximum instantaneous deficit rate.

Reliability ratio = served ÷ required volume, in [0, 1] (reported per node and for the whole network).

Options — DemandReliabilityOptions { deficit_tolerance }: shortfalls below this rate (m³/s) are not counted as deficit periods, though their volume still accumulates. Default 1e-9.

Units: volumes in m³, rates in m³/s, times in seconds (internal SI). The report also records the model’s demand model (DDA or PDA) — under DDA delivered ≈ required, so reliability is mainly meaningful under PDA or when checking a deficit scenario.

Service compliance (SDK)

Measures how often junction pressures stayed within acceptable bounds.

Entry point (see SDK Examples):

#![allow(unused)]
fn main() {
compute_service_compliance_from_out(out_path, thresholds) -> ServiceComplianceReport
}

It reads only the .out file (junction membership comes from the file’s node table). Only junction nodes are analysed — reservoirs and tanks are excluded so they don’t register as permanent violations.

Thresholds — ServiceComplianceThresholds { min_pressure, max_pressure }: min_pressure is required; max_pressure is optional (None disables the upper-bound check). Use ServiceComplianceThresholds::min_only(p) for a lower bound only. When set, max_pressure must be strictly greater than min_pressure.

What it computes, per period and junction:

• Sample counts: within limits, below minimum, above maximum.
• Deficit / excess integrals — pressure shortfall or excess accumulated over time (m·s).
• Worst observed deficit and excess (m), and the longest consecutive violation streak.

The summary aggregates these across all junctions and periods and reports a compliance ratio (fraction of in-limit samples). Pressures use the units stored in the .out file (metres of head for Hydra output).

The GUI Analysis tab

The Analysis tab computes its own dashboard from the scenario’s results (a histogram/summary pass, independent of the two SDK modules above). It has six panels:

Availability

Diagnostics & Errors

Hydra reports problems in three ways: exit codes (CLI process status), structured diagnostics (machine-readable JSON lines on stderr), and warnings (non-fatal issues recorded during a run). This page catalogues each.

Scope. The structured diagnostics below are emitted by the simulation command (hydra <input>). The hydra report subcommand shares the same exit codes but writes plain error: … text to stderr rather than JSON lines, so do not parse its stderr as diagnostics.

Exit codes (CLI)

Codes 0–4 are a stable contract. Note that network validation failures are input errors (exit 1), not solver errors.

Structured diagnostics (stderr)

Independently of the report, the CLI writes each warning and error to stderr as one compact JSON object per line. This stream is always emitted — it is not suppressed by -q/--quiet (which only silences the human-readable progress output).

Each line has five keys:

Example:

{"level":"warning","code":"warning/negative_pressure","message":"negative pressure at node 'J1'","object_id":"J1","time_step":3600.0}
{"level":"error","code":"solver/hydraulic","message":"...","object_id":null,"time_step":null}

Diagnostic codes

Warnings (non-fatal; the run continues):

Errors (fatal; the run stops):

Warnings

Warnings are recorded during a run and surfaced in three places: the stderr diagnostics above, the warnings array of the JSON report, and the Warnings section of the text report. There are three kinds:

• Unbalanced hydraulics — a step did not converge. Often points to an over-constrained or disconnected model; see Troubleshooting.
• Negative pressure — a node’s pressure went below zero, indicating demand that the network could not supply at that point (consider PDA).
• Pump exceeds maximum head — a pump operated past the top of its curve; check the curve and the operating point.

Validation errors

Before a simulation runs, Hydra validates the network and reports all structural problems it finds (not just the first). Each is emitted under the validation/network diagnostic code. The checks cover:

The Issues panel in the GUI surfaces these findings with links to the affected elements.

Migrating from EPANET

This page is for engineers and developers switching from EPANET to Hydra. It covers what works out of the box, what to expect numerically, and where behaviour intentionally differs.

Your .inp Files Work

Hydra parses the EPANET .inp format directly — any 2.x release. No conversion is needed. Pass your existing .inp file to the CLI or the library and Hydra will run it.

The command line is Hydra’s own, not EPANET’s. Hydra deliberately does not mimic epanet input.inp report.rpt output.out: that argument order encodes one engine and one pair of artifacts, which stops being true as Hydra adds engines. File-format compatibility and command-line compatibility are separate promises, and Hydra keeps the first.

See INP Format Support for the full section-by-section reference.

Output Formats

Expect Small Numerical Differences

Hydra and EPANET solve the same physics using the same Global Gradient Algorithm, but they follow independent numerical paths. On most networks you will see differences of less than 0.1% in head and flow values. These are not bugs; they are the expected consequence of floating-point arithmetic being non-associative.

The practical impact depends on network topology. Simple, stable networks with few controls and no quality agree to well within a rounding error. Differences grow with the number of demand nodes and control switches, because a step that lands either side of a control threshold changes what happens next.

Quality results are the most sensitive, and for a structural reason rather than a numerical one: quality integrates the hydraulic solution. A flow difference too small to notice in heads is carried into transport, where it compounds across periods. If your workflow depends on sub-percent quality agreement with EPANET output, treat both results as independent estimates of the same physical system; neither is more “correct” than the other in an absolute sense.

Hydra’s result is authoritative. If you observe a difference and suspect a Hydra bug, open a GitHub issue with a minimal reproducer.

Behavioural Differences

Unbalanced-stop mode

EPANET halts the simulation when a hydraulic step does not converge within the configured iteration limit (UNBALANCED STOP). Hydra honours this setting: when a hydraulic step is genuinely unbalanced (fails to converge), Hydra also halts and records an UnbalancedHydraulics warning. The UNBALANCED CONTINUE N option is also supported.

Because the two engines follow independent numerical paths, the step at which non-convergence first occurs can differ — so the same model may halt at different periods, or converge throughout in one engine and stop partway in the other, even though both apply the same rule.

There is a second, harder stop in both engines: if the linear system becomes singular and no control valve can be demoted to recover it, the run aborts with no result saved for that step. That differs from the unbalanced stop, which saves the failing step before ending.

Quality timestep handling

EPANET’s quality timestep can reach 0 seconds via integer truncation when hydraulic timesteps are very short. Hydra keeps the quality timestep as a real number and, when it is 0 or unset, defaults it to one-tenth of the hydraulic timestep, so it never truncates to zero. An explicitly set sub-second step is used as given.

This only matters for networks with very short hydraulic timesteps (well under 60 seconds), which is unusual in practice.

Newer EPANET Features Worth Knowing

Both are fully supported by Hydra, and both are optional — a file that uses neither still loads and runs.

FAVAD Leakage — OWA-EPANET 2.3

Per-pipe background leakage is modelled using the FAVAD (Fixed and Variable Area Discharge) model, configured via a [LEAKAGE] section in the .inp file. This section is the one genuine 2.3 addition. Older files (without [LEAKAGE]) parse cleanly; leakage is simply zero for all pipes.

Pressure-Dependent Analysis — EPA EPANET 2.2

PDA is configured exactly as in EPANET (DEMAND MODEL PDA in [OPTIONS], with MINIMUM PRESSURE, REQUIRED PRESSURE, and PRESSURE EXPONENT). No changes needed.

EPANET API Mapping

If you are migrating code that uses the EPANET Toolkit C API, the equivalent Hydra library workflow is:

See the SDK overview for complete library usage examples.

Performance

Hydra is fast enough to run large extended-period simulations interactively. The figures below are indicative — they show the order of magnitude you can expect, not a controlled lab result. Your own hardware, build flags, and network characteristics will change the numbers.

End-to-end timings

Each figure is the best of 5 runs (after a warm-up) of the full command-line workflow on one network: parse the .inp, run the complete extended-period simulation, and produce the summary report. Networks are the ones bundled in tests/benchmarks/.

Measured on an Apple M5 Pro (macOS 26.5), release build (cargo build --release, fat LTO, codegen-units = 1).

Steps is the number of reporting periods. Single-period (Steps: 1) rows are steady-state snapshots — their times are dominated by process start-up and parsing rather than the solve, so treat any sub-10 ms figure as “effectively instant.” The larger the network, the more of each hydraulic step is spent in the sparse linear solve; networks with heavy control/rule logic or water-quality transport (Micropolis, for example) cost more per step than their node count alone suggests.

Reproducing these numbers

The table is generated by a committed harness, so it can be regenerated rather than going stale:

just bench-report

This builds the release CLI and runs scripts/benchmark.py, which times each network and prints the Markdown table above. Pass --runs N to change the sample count, or --hydra PATH to benchmark a specific binary.

Building for maximum speed

The release profile already enables fat LTO and a single codegen unit. For the best local performance, build with native CPU features:

just release-native

This tunes the binary for the machine it is built on (-C target-cpu=native); such binaries are not portable to older CPUs.

Solver micro-benchmarks

For work on the solver itself, the criterion suite times the hydraulic solve step (warm and cold) in isolation:

just bench

SDK Overview

hydra-sdk is the umbrella crate for Hydra’s public API. Add it to your Cargo.toml:

[dependencies]
hydra-sdk = "6"

It re-exports every type needed to parse networks, run simulations, query results, run post-simulation analytics, and generate reports — with all internal dependency versions pre-pinned.

Modules and Key Types

Session API

The primary entry point. Import Simulation to parse, run, and query a network.

Analytics

Post-simulation analysis functions that operate on a saved .out file.

Data Model

The full network data model, mirroring the EPANET .inp structure.

I/O

#![allow(unused)]
fn main() {
use hydra_sdk::io;
}

Engine Identity

#![allow(unused)]
fn main() {
use hydra_sdk::common;
}

Every Hydra engine publishes an immutable descriptor. Applications resolve a project’s stored engine key against the registry rather than hardcoding names, colours, or file filters. See Engines for the current roster and what Planned means in practice.

import is a file-picker filter, never a validity test — wds and uds both claim .inp with incompatible contents, so only the owning engine’s parser can decide whether a file really is its model.

Reports

#![allow(unused)]
fn main() {
use hydra_sdk::{report, report_catalog, produce_report_block};
}

Report generation is split in two: the engine produces neutral content fragments, and report turns them into documents. The report layer knows nothing about engines.

Also re-exported

Beyond the tables above, hydra-sdk re-exports several supporting items:

• Version constants — HYDRA_VERSION and the per-subsystem HYDRA_*_VERSION strings.
• Runtime estimation — estimate_simulation_runtime, estimate_simulation_runtime_from_summary, and RuntimeEstimate. The millisecond-level forms estimate_simulation_runtime_millis_from_summary and classify_simulation_runtime_millis are also available when you want the raw prediction or the bucketing separately.
• Threshold binning — threshold_bands(values, edges) counts values into the bands defined by ascending edges, with the outer two unbounded so nothing is dropped. It is the same binning the wds.*-thresholds report blocks use, so an interface presenting that view counts identically.
• Threshold binning — threshold_bands, the shared band-counting used by the *-thresholds report blocks, so an interface presenting the same view counts identically.

SDK Examples

Parse an INP file and run a full simulation

use hydra_sdk::{io, Simulation, NodeQuantity, LinkQuantity};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bytes = std::fs::read("network.inp")?;
    let network = io::parse(&bytes)?;

    let mut sim = Simulation::create();
    sim.load(network)?;
    sim.run()?;

    for t in sim.snapshot_times() {
        let head = sim.get_node_result("J1", NodeQuantity::Head, t)?;
        let pressure = sim.get_node_result("J1", NodeQuantity::GaugePressure, t)?;
        let flow = sim.get_link_result("P1", LinkQuantity::Flow, t)?;
        println!("t={t:.0}s  head={head:.3}  pressure={pressure:.3}  flow={flow:.6}");
    }

    for w in sim.warnings() {
        println!("[t={:.0}s] {:?}", w.t, w.kind);
    }

    Ok(())
}

Shorthand constructor

use hydra_sdk::{io, Simulation};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bytes = std::fs::read("network.inp")?;
    let network = io::parse(&bytes)?;

    // Convenience: shorthand for Simulation::create() + sim.load(network).
    let mut sim = Simulation::from_network(network)?;
    sim.run()?;

    Ok(())
}

Step through hydraulics manually

use hydra_sdk::{io, Simulation};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bytes = std::fs::read("network.inp")?;
    let network = io::parse(&bytes)?;

    let mut sim = Simulation::create();
    sim.load(network)?;

    loop {
        let dt = sim.step_hydraulics()?;
        if dt == 0.0 { break; }
        // inspect or modify state between steps
    }

    Ok(())
}

Write output files

use hydra_sdk::{io, Simulation, WritableSimulation};
use std::fs::File;
use std::io::BufWriter;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bytes = std::fs::read("network.inp")?;
    let network = io::parse(&bytes)?;
    let mut sim = Simulation::from_network(network)?;
    sim.run()?;

    // Plain-text .rpt report
    let rpt = io::rpt_writer::build_text_report(&sim)?;
    std::fs::write("report.rpt", rpt)?;

    // JSON report
    let json = io::rpt_writer::build_json_report(&sim)?;
    std::fs::write("report.json", json)?;

    // EPANET-compatible binary .out file. The last three arguments are the
    // input/report paths recorded in the file prolog and the output unit system
    // (use the model's declared units via `sim.net().options.flow_units`).
    let mut out_file = BufWriter::new(File::create("output.out")?);
    io::out_writer::write_binary_output(
        &mut out_file,
        &sim,
        "network.inp",
        "report.rpt",
        sim.net().options.flow_units,
    )?;

    Ok(())
}

Demand reliability analysis

Post-simulation analytics operate on a saved .out file and the original Network.

use hydra_sdk::{io, compute_demand_reliability_from_out};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let bytes = std::fs::read("network.inp")?;
    let network = io::parse(&bytes)?;

    let report = compute_demand_reliability_from_out(
        std::path::Path::new("output.out"),
        &network,
    )?;

    println!("Network reliability: {:.1}%",
        report.summary.reliability_ratio() * 100.0);

    for node in &report.nodes {
        if node.reliability_ratio() < 0.99 {
            println!(
                "  {} — {:.1}% reliable, {} deficit period(s)",
                node.node_id,
                node.reliability_ratio() * 100.0,
                node.deficit_periods,
            );
        }
    }

    Ok(())
}

Pressure compliance analysis

use hydra_sdk::{compute_service_compliance_from_out, ServiceComplianceThresholds};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // Check that all nodes stay above 10 m and below 80 m.
    let thresholds = ServiceComplianceThresholds {
        min_pressure: 10.0,
        max_pressure: Some(80.0),
    };

    let report = compute_service_compliance_from_out(
        std::path::Path::new("output.out"),
        thresholds,
    )?;

    let compliant = report.nodes.iter()
        .filter(|n| n.below_min_count == 0 && n.above_max_count == 0)
        .count();
    println!("Compliant nodes: {}/{}", compliant, report.nodes.len());

    for node in &report.nodes {
        if node.below_min_count > 0 {
            println!(
                "  node {} — {} period(s) below {}m (worst deficit: {:.2}m)",
                node.node_index,
                node.below_min_count,
                thresholds.min_pressure,
                node.worst_below_min,
            );
        }
    }

    Ok(())
}

Crate Layout

Hydra is a multi-crate Rust workspace:

The two empty engine scaffolds exist so their crate names and versions track the workspace from the start, rather than being introduced mid-life — see Engines for what each engine covers. The split between hydra-common, the engines, and hydra-report is what lets a report be assembled from any engine’s output: engines emit neutral fragments, and the report layer renders them without knowing which engine produced them.

hydra-cli and hydra-gui are downstream consumers of Hydra in exactly the same way a third-party integrator would be: they depend on the umbrella crate and never import from hydra-engine-wds directly. Anyone who wants a different interface (HTTP, gRPC, Python bindings, etc.) follows the same pattern.

Specifications

Every engine is specified subsystem by subsystem, alongside the cross-cutting specs for the shared foundation and the report layer. These documents are the authoritative definitions of Hydra’s behaviour — where a spec and the implementation disagree, the spec wins.

Shared

Water Distribution (wds)

Urban Drainage (uds) and Open Channel (och)

Not yet written. Both engines are registered but planned — no behaviour has been specified, so there is nothing yet for an implementation to conform to. Their specs land before any implementation code does, as the spec-first workflow requires.

EPANET: A Conceptual and Mathematical Analysis

This document analyses EPANET, not Hydra. It is a pinned snapshot: it describes EPANET 2.3.5 at tag v2.3.5 faithfully, defects included. Where it records a bug or an inconsistency, that is a finding about EPANET — not a Hydra issue to file.

Hydra’s water distribution engine was derived from this analysis but is a distinct solver, not a reimplementation. Its own authoritative behaviour is defined by the specifications.

Introduction

OWA-EPANET is a computational engine for simulating the hydraulic and water quality behaviour of pressurised water distribution networks over time. It represents a network as a directed graph of nodes connected by links and advances a time-stepped extended-period simulation, solving at each step for pressures, flows, and constituent concentrations throughout the system. The solver combines rigorous physical models — empirical head-loss formulas, Newton–Raphson linearisation, sparse direct linear algebra, and Lagrangian advection-reaction transport — with flexible engineering constructs such as demand patterns, operational controls, and tank mixing models.

This document provides a self-contained, mathematical and conceptual description of every major subsystem: how the network is represented, how hydraulic equilibrium is computed at each time step, how demands are handled under both fixed and pressure-dependent conditions, how leakage and emitter flows enter the system, how the simulation advances through time, how control logic operates, how water quality is transported and reacted, and how energy and mass balance are tracked. The goal is to give the reader a complete algorithmic and mathematical understanding of the system; implementation-specific details such as memory layout and data-structure internals are omitted, but input/output behaviour and the public API are described at a conceptual level. The analysis describes EPANET 2.3.5 — tag v2.3.5 of the OWA-EPANET repository — and was derived from its C source, which is authoritative throughout.

Table of Contents

• EPANET: A Conceptual and Mathematical Analysis
  • Introduction
  • Table of Contents
  • 1. Network Representation
    • Node Types
    • Link Types
    • Curves and Patterns
    • The State-Vector View
  • 2. Hydraulic Simulation
    • 2.1 Head Loss in Pipes
      • Hazen–Williams Formula
      • Darcy–Weisbach Formula
      • Chezy–Manning Formula
      • Minor Losses and Total Head Loss
    • 2.2 Pump Head Gain
    • 2.3 Valve Behaviour
    • 2.4 The Global Gradient Algorithm
      • Governing Equations
      • Linearisation
      • Assembly of the Linear System
      • Flow Update
      • Convergence Criterion
    • 2.5 Sparse Linear Algebra
  • 3. Demand Models
    • 3.1 Demand-Driven Analysis
    • 3.2 Pressure-Driven Analysis
  • 4. Emitters
  • 5. Pipe Leakage — the FAVAD Model
  • 6. Time-Stepping and Tank Dynamics
    • 6.1 Extended-Period Simulation
    • 6.2 Adaptive Time Step
    • 6.3 Tank Level Update
  • 7. Control Systems
    • 7.1 Simple Controls
    • 7.2 Rule-Based Controls
  • 8. Water Quality Simulation
    • 8.1 Overview and Simulation Modes
    • 8.2 Lagrangian Segment Transport
    • 8.3 Source Terms
    • 8.4 Chemical Reactions
      • Bulk Reactions
      • Wall Reactions
      • Combined Reaction in a Segment
      • Roughness–Reaction Correlation
  • 9. Tank Mixing Models
    • Complete Mix (CSTR)
    • Two-Compartment Mix
    • FIFO Plug Flow
    • LIFO (Stacked Layers)
  • 10. Mass Balance
  • 11. Energy Tracking
  • 12. Flow Balance
  • 13. Units and Physical Constants
  • 14. Input and Output
    • Input
    • Output
    • Binary Output File Format
      • Prolog
      • Energy
      • Dynamic Results
      • Network Reactions
      • Epilog
    • API
  • 15. Cross-Cutting Engine Contracts

1. Network Representation

The physical infrastructure is represented as a directed graph. Nodes correspond to points in the network — junctions, reservoirs, and tanks — and links correspond to the conduits and devices connecting them — pipes, pumps, and valves. The orientation of a link defines a positive flow direction; negative flows simply indicate flow in the reverse direction.

Node Types

Junctions are the ordinary connection points of the network. Each junction has a fixed elevation and one or more demand categories, each associating a base demand rate with a time-varying multiplier pattern. Demand lines in the [DEMANDS] section replace the demand entered on the junction’s own line — the first entry overwrites, later entries for the same node append further categories — and demand entries naming tanks or reservoirs are silently ignored. Junctions may also carry an emitter (representing orifice or sprinkler outflow) and a water quality source. The hydraulic head at a junction is an unknown to be solved at every time step.

Reservoirs are fixed-grade nodes whose hydraulic head is always known and equal to the water surface elevation. A reservoir represents an infinitely large storage body that maintains a constant pressure boundary condition. Because its head is known, it does not appear as an unknown in the linear system; instead, it contributes boundary terms to the equations of its neighbouring junctions.

Tanks are storage nodes with a variable water level that evolves as water flows in and out. Each tank is characterised by a minimum level, a maximum level, an initial level, and a geometry: either a constant cross-sectional area or a user-defined volume-versus-elevation curve (an optional explicit minimum volume overrides the cylindrical estimate $A \times h_{\min}$; the tank’s levels must lie within a volume curve’s elevation range). Its hydraulic head at any instant equals the elevation of its water surface, which is updated after each hydraulic time step based on the net flow. A tank also carries a bulk reaction coefficient governing water quality transformations in its stored volume.

Link Types

Pipes are the primary conduits. Each pipe is characterised by its length, internal diameter, a roughness coefficient (whose interpretation depends on the chosen head-loss formula), a minor loss coefficient representing localised losses at fittings, and bulk and wall reaction coefficients for quality simulation. A pipe may also be designated as a check valve, in which case flow is permitted in only one direction; the link is treated as closed whenever the computed flow or pressure gradient would drive flow in the reverse direction.

Pumps add hydraulic head to the flow passing through them. The head added is described by a head-versus-flow curve, one of the most important user-supplied relationships in the model. Alternatively, a pump may be defined by a constant power output. A variable-speed setting scales both the head and flow axes of the pump curve according to the affinity laws.

Valves regulate flow or pressure and come in seven varieties:

• Pressure Reducing Valve (PRV): limits the hydraulic head on its downstream side to a specified setpoint when the upstream head exceeds it.
• Pressure Sustaining Valve (PSV): maintains the hydraulic head on its upstream side above a specified setpoint when the downstream head would otherwise pull it below.
• Flow Control Valve (FCV): restricts the volumetric flow through it to a specified setpoint.
• Throttle Control Valve (TCV): applies a specified head loss coefficient; it behaves as a pipe with adjustable resistance.
• General Purpose Valve (GPV): head loss as a function of flow is entirely described by a user-supplied piece-wise linear curve.
• Positional Control Valve (PCV): The loss coefficient varies with a percent-open setting, optionally governed by a user-supplied curve relating the valve opening (percent open) to the ratio of its flow coefficient to its fully-open flow coefficient, $K_v/K_{v0}$. The internal minor-loss coefficient is then derived from this ratio as $K_{m0}/(K_v/K_{v0})^2$ (see §2.3); absent a curve, the ratio is simply linear in percent open.
• Pressure Breaker Valve (PBV): imposes a fixed head-loss setpoint. When the setting exceeds the natural minor-loss head drop at the current flow, the solver forces the exact head loss; otherwise the PBV falls back to ordinary pipe resistance. PBVs have no control states — they always contribute a resistance.

Valve placement is validated at parse time: a PRV, PSV, or FCV may not connect to a tank or reservoir, and specific adjacent combinations are prohibited — two PRVs sharing a downstream node or in series, two PSVs sharing an upstream node or in series, and certain PRV/PSV/FCV adjacencies.

Curves and Patterns

User-defined curves are piece-wise linear relationships used throughout the model: pump head versus flow, pump efficiency versus flow, tank volume versus elevation, general purpose valve head loss versus flow, and positional valve opening versus flow-coefficient ratio ($K_v/K_{v0}$). Intermediate values are obtained by linear interpolation between the two bracketing data points. Behaviour outside a curve’s x-range depends on the evaluation path: the solver’s segment evaluation for custom pump head curves and general purpose valve curves extends the first or last segment as a straight line — extrapolating linearly beyond the end points — whereas the general lookup used for tank volume, pump efficiency, and similar curves clamps to the end-point values and never extrapolates. A curve’s type may be declared with an explicit tag or is otherwise inferred from its first use.

Patterns are repeating sequences of dimensionless multipliers indexed by time. They modulate base demands at junctions, pump speed settings, and constituent source concentrations over the course of the simulation. At each hydraulic time step the pattern multiplier is determined by the pattern period index: given a global pattern start offset $t_{\text{start}}$, pattern time step $\Delta t_p$, and current simulation time $t$, the number of elapsed periods is $p = \lfloor (t + t_{\text{start}}) / \Delta t_p \rfloor$. For a demand category assigned to pattern $j$ of length $L_j$, the applicable multiplier is $F_j[p \bmod L_j]$, giving each pattern an independently repeating cycle. A default pattern is available and is applied to any demand category that has no explicit pattern assigned; if no default pattern exists either, a multiplier of 1.0 is used.

Reservoir head patterns: while reservoirs are nominally fixed-grade nodes, each reservoir may optionally be assigned a time pattern. When assigned, the head at that reservoir at each hydraulic time step equals its base elevation multiplied by the current pattern multiplier. This allows time-varying source heads to represent, for example, tidal fluctuations or varying water tower levels.

Pump utilisation patterns: each pump may have a separate utilisation pattern (distinct from any energy cost pattern) that controls the pump’s speed setting at each hydraulic time step. The current pattern multiplier is applied directly as the normalised speed $\omega$; a multiplier of zero closes the pump, while a multiplier of 1.0 sets it to its rated speed. This allows pump schedules to be encoded as a time series without requiring explicit control rules.

The State-Vector View

Because each hydraulic time step is solved as a steady state (§6.1), the quantities that genuinely evolve from one step to the next form a small state vector. On the hydraulic side it comprises: the tank water levels — equivalently the stored volumes — integrated forward after each solve (§6.3), the only continuously-evolving physical state; every link’s discrete status and numeric setting, as modified by controls and rules (§7) or by pattern-driven pump speed settings (see above); and the link flows themselves, which carry over as the estimates around which the next solve’s linearisation begins — §2.2 gives the cold-start seeding used before the first solve, and on re-initialisation existing link flows are preserved as a warm start unless re-initialisation is explicitly requested or the flow is near zero (§6.1). Everything else the hydraulic engine reports — junction heads, pressures, equilibrium link flows, velocities — is derived: recomputed at every step by the Global Gradient Algorithm (§2.4) from the boundary conditions this state defines, namely the known heads of reservoirs (possibly pattern-scaled) and tanks at their current levels, together with the pattern-scaled demands.

The water-quality engine carries its own state between quality sub-steps: the ordered segment lists of each pipe, each segment holding a volume and a uniform concentration (§8.2); the mixing state of each tank as prescribed by its mixing model — a single uniform concentration for a complete-mix tank, the volumes and concentrations of the mixing and stagnant zones for a two-compartment tank, and the segment stacks of FIFO and LIFO tanks (§9); and the running mass-balance ledger that accounts for source mass added, mass removed as outflow, and mass reacted (§10).

The seam between the two engines is explicit: at each hydraulic step the solver writes nodal demands, heads, link flows, statuses, and settings to a hydraulics binary file, and the quality simulation replays this file — velocities are reconstructed from the saved flows and pipe geometry, and the flow field is held constant across the quality sub-steps within each hydraulic period (§8.1, §14). The saved record is therefore the complete hydraulic input the quality engine ever sees, which is why a saved hydraulics file can be reused across runs without recomputing the hydraulics (§14).

2. Hydraulic Simulation

2.1 Head Loss in Pipes

Hydraulic head at any node is the mechanical energy per unit weight of water:

$$H = \frac{P}{\rho g} + z$$

where $P$ is the gauge pressure, $\rho$ is the water density, $g$ is gravitational acceleration, and $z$ is the elevation above datum. Flow from node $i$ to node $j$ is driven by the head difference $H_i - H_j$.

Three empirical head-loss formulas are available, and one is selected uniformly for the entire network.

Hazen–Williams Formula

$$h_f = \frac{4.727 , L}{C^{1.852} , D^{4.871}} , Q^{1.852}$$

Here $L$ is the pipe length, $D$ is the internal diameter, $C$ is the Hazen–Williams roughness coefficient (higher values indicate smoother pipes), and $Q$ is the volumetric flow rate. The flow exponent is $n = 1.852$. This formula is empirical and strictly valid only for turbulent flow of water at ordinary temperatures.

Darcy–Weisbach Formula

$$h_f = f \cdot \frac{L}{D} \cdot \frac{V^2}{2g} = f \cdot \frac{8 L}{\pi^2 g D^5} , Q^2$$

where $V = Q / (\pi D^2 / 4)$ is the mean flow velocity and $f$ is the dimensionless Darcy friction factor, which depends on the Reynolds number $Re = VD/\nu$ and the relative roughness $\varepsilon/D$.

For laminar flow ($Re \leq 2000$) the Hagen–Poiseuille result applies:

$$f = \frac{64}{Re}$$

yielding a head loss proportional to $Q$ (linear regime).

For turbulent flow ($Re \geq 4000$) the friction factor is computed from the Swamee–Jain approximation to the Colebrook-White implicit equation:

$$f = \left[ -2 \log!\left( \frac{\varepsilon}{3.7 D} + \frac{5.74}{Re^{0.9}} \right) \right]^{-2}$$

where $\varepsilon$ is the absolute roughness. The quantity $f$ and its derivative with respect to $Q$ are evaluated simultaneously at each Newton iteration so that the linearisation of the solver (§2.4) remains consistent.

For transitional flow ($2000 < Re < 4000$), a cubic polynomial ensures continuity in $f$ and $df/dQ$ across the transition. The polynomial is anchored at both ends: $f = 64/Re$ at $Re = 2000$ (exact laminar value), and the Swamee–Jain value and its derivative at $Re = 4000$ (turbulent end). The laminar Hagen–Poiseuille branch handles flows below a pipe-geometry-dependent low-flow threshold independently and does not call this cubic.

Chezy–Manning Formula

$$h_f = \left( \frac{4 n_M}{1.49 , \pi , D^2} \right)^2 \left( \frac{D}{4} \right)^{-1.333} L , Q^2$$

where $n_M$ is the Manning roughness coefficient. (The source hardcodes the empirical constant as $1.49$, a two-figure rounding of the textbook Manning-Strickler value $1.486$, and the exponent as $-1.333$, a truncation of $-4/3$ — both roundings preserved here as the code’s actual behaviour.) The exponent on $Q$ is 2 in this formulation (as used for full circular pipes). This formula is less common for pressurised systems but is supported for completeness.

Minor Losses and Total Head Loss

Minor (local) losses due to fittings, bends, and contractions are modelled as:

$$h_{\text{minor}} = K_m , Q , |Q|$$

where $K_m$ is the minor loss coefficient (head loss per unit of $Q^2$). The sign convention ensures the loss opposes flow in either direction.

The total head loss across a pipe combining friction and minor losses is:

$$h = R , Q^n \cdot \mathrm{sign}(Q) + K_m , Q , |Q|$$

where $R$ is the friction resistance coefficient derived from whichever formula is in use and $n$ is the corresponding flow exponent (1.852 for Hazen–Williams, 2 for Darcy–Weisbach and Chezy–Manning in their simplified forms). The sign convention ensures that the expression is an odd function of $Q$: head loss is always in the direction opposing flow.

To keep the Jacobian non-singular near zero flow, whenever an element’s head-loss gradient falls below the RQTOL option (default $10^{-7}$), the gradient is floored and the element’s head loss becomes a linear function of flow — applied to Hazen–Williams and Chezy–Manning pipes, nonlinear pump curves, valve minor losses (floored at $RQtol/2$), and emitters; Darcy–Weisbach pipes rely on their laminar branch instead.

2.2 Pump Head Gain

A pump adds head to the flow. Three types of pump curves are supported.

Power-function curve: the head gain follows

$$\Delta H = h_0 - r , Q^N$$

where $h_0$ is the shutoff head (head at zero flow), $r$ is a resistance-like coefficient, and $N$ is the curve exponent. This three-parameter form fits most centrifugal pump characteristics well.

Constant-power pump: the head gain is determined by maintaining a fixed power output regardless of flow:

$$\Delta H = \frac{\text{Power}}{\gamma , Q}$$

where $\gamma = \rho g$ is the specific weight of water. As flow decreases toward zero, the head gain grows without bound. The solver handles this by monitoring the head-loss gradient $|\partial h / \partial Q| = r / Q^2$: when this gradient exceeds $C_\infty \approx 10^8$ (near-zero flow), the pump is treated as a closed link ($P_k = 1/C_\infty$, $Y_k = Q_k$); when the gradient falls below $10^{-6}$ (extremely high flow), the pump is treated as a fully open link with minimal resistance. Between these extremes, the standard linearisation $h = r/Q$ applies. Independently, a constant-power pump whose flow falls below $10^{-6}$ ft³/s is set to TEMPCLOSED status (see pump status below), because the power formula is undefined at zero flow. A constant-power pump’s nominal design flow is fixed at 1 ft³/s and is restored whenever the pump is re-opened by a status change, setting change, or control, so the $r/Q$ linearisation restarts from a sane point; and during the flow update, a Newton correction that would drive its flow strictly negative is replaced by halving the current flow instead (a correction landing exactly at zero flow is allowed through).

Custom curve: a user-defined piece-wise linear head-versus-flow curve. At any operating point the solver identifies the two adjacent data points that bracket the current flow and interpolates linearly to obtain the head gain and its derivative. For initial flow conditions before the first solve, the design flow $Q_0$ for a custom-curve pump is taken as the midpoint between the first and last flow data points on the curve (rather than a named design point).

Speed scaling via affinity laws: when a pump operates at a relative speed $\omega$ (with $\omega = 1$ being its rated speed), the affinity laws relate the scaled curve to the rated curve:

$$\Delta H(\omega, Q) = \omega^2 \cdot \Delta H_1!\left(\frac{Q}{\omega}\right)$$

where $\Delta H_1$ is the head gain at rated speed. Equivalently, the shutoff head scales as $\omega^2$ and the flow axis scales as $\omega$. In the Newton–Raphson solver, a pump is treated as a link with a negative head-loss value (a gain), and its linearised resistance coefficient $P_k$ and offset $Y_k$ are derived from the pump curve in the same algebraic framework as pipe head losses.

Three-point pump curve fitting: when a pump is specified by three operating points — the shutoff head $h_0$ (head at zero flow), the design point $(q_1, h_1)$, and the maximum-flow point $(q_2, h_2)$ — the power-function parameters are determined analytically:

$$c = \frac{\ln!\left(\dfrac{h_0 - h_2}{h_0 - h_1}\right)}{\ln!\left(\dfrac{q_2}{q_1}\right)}, \qquad b = \frac{h_0 - h_1}{q_1^{,c}}, \qquad a = h_0$$

yielding the curve $\Delta H = a - b , Q^c$. The curve is validated: it must be strictly decreasing in head ($h_0 > h_1 > h_2$), the exponent $c$ must be positive, and additionally $c \leq 20$ (an upper-bound sanity check enforced by the validator). A one-point curve $(q_1, h_1)$ is expanded to this form using $h_0 = 1.33334,h_1$, $q_2 = 2q_1$, $h_2 = 0$; the power-function fit applies only to curves of exactly one or exactly three points (the latter with their first point at zero flow) — any other curve, whether by point count or by a nonzero first flow, is treated as a custom piece-wise curve instead. Validation also derives each pump’s limits: power curves take $Q_{\max}$ as the flow at which the fitted head reaches zero and $H_{\max} = h_0$, while custom curves take their last flow and first head points — this $Q_{\max}$ is the one referenced by the XFLOW warning (§2.4).

Pump status — XHEAD: a pump in the OPEN state transitions to the XHEAD (excess head) state when the head gain required to maintain the computed flow exceeds the speed-adjusted head limit $\omega^2 H_{\max}$ — the shutoff head $h_0$ for power-function curves, the first head point for custom curves (see the curve-limit derivation above) — by more than the head tolerance $\text{Htol}$ (§2.3). In this state the pump is treated as a closed link for that iteration. The status reverts to OPEN at the start of each periodic status check, and is re-tested against the new computed operating point. For constant-power pumps, XHEAD cannot occur; instead, the pump is flagged TEMPCLOSED when the flow falls below $10^{-6}$ ft³/s (a strictly positive cutoff, since the power formula is undefined at zero flow).

Initial flow conditions: before the first Newton–Raphson solve, link flows are initialised as follows — closed links receive a negligible flow $Q_0 \approx 10^{-6}$ ft³/s; pumps receive the product of their speed setting and their design flow $Q_{\text{design}}$; all other links (pipes and valves) receive the flow corresponding to a nominal velocity of 1 ft/s through the full pipe cross-section: $Q = \pi D^2 / 4$. These initial values need not be physically consistent; the Newton–Raphson iteration converges from them to the true solution.

2.3 Valve Behaviour

The three control valves — PRV, PSV, and FCV — can inhabit one of three principal discrete states at any iteration — active, open, or closed — plus the exception states XFCV and XPRESSURE described below.

• Active: the valve enforces its design constraint. A PRV fixes the downstream head equal to its setpoint; a PSV fixes the upstream head equal to its setpoint; an FCV fixes the flow through it equal to its setpoint. When active, these valves introduce a head constraint rather than a resistance relationship, and the corresponding row of the linear system is modified accordingly.
• Open: the valve is behaving as a short section of pipe with negligible resistance; no constraint is enforced.
• Closed: the valve passes no flow.

After each Newton–Raphson iteration the hydraulic state of each control valve is examined:

• A PRV transitions from ACTIVE to OPEN if the upstream head has fallen to or below the setpoint (no pressure reduction needed), or to CLOSED if the flow through it reverses (maintaining the setpoint would require reverse flow).
• A PSV transitions from ACTIVE to OPEN if the downstream head has risen to or above the setpoint, or to CLOSED if the flow through it reverses.
• An FCV transitions to the XFCV state (it cannot enforce its setpoint) if the available head difference is insufficient to sustain the target flow, or if the flow through it turns negative.

TCV, GPV, and PCV valves do not have control states; they always contribute a resistance (head loss as a function of flow) determined by their current setting.

Precise valve status transitions: PRV and PSV states are re-evaluated after each Newton–Raphson iteration (governed by the DampLimit parameter; see §2.4); FCV transitions instead belong to the periodic linkstatus schedule (CheckFreq/MaxCheck, plus the convergence-time pass). The transition rules are:

• PRV: in the ACTIVE and OPEN states the reverse-flow test is evaluated first and takes precedence — ACTIVE → CLOSED and OPEN → CLOSED if $Q < -\varepsilon_Q$ (reverse flow). Only otherwise: ACTIVE → OPEN if $H_1 - K_m Q^2 < H_\text{set} - \varepsilon_H$ (upstream pressure insufficient to need reduction); OPEN → ACTIVE if $H_2 \geq H_\text{set} + \varepsilon_H$ (downstream pressure reaches setpoint). CLOSED → ACTIVE if $H_1 \geq H_\text{set} + \varepsilon_H$ and $H_2 < H_\text{set} - \varepsilon_H$; CLOSED → OPEN if $H_1 < H_\text{set} - \varepsilon_H$ and $H_1 > H_2 + \varepsilon_H$. The special XPRESSURE state (entered only when an active valve renders the solution matrix ill-conditioned; reported as open but unable to deliver pressure) transitions to CLOSED on reverse flow.
• PSV: symmetric to PRV — the reverse-flow test again runs first: ACTIVE → CLOSED and OPEN → CLOSED if $Q < -\varepsilon_Q$; only otherwise ACTIVE → OPEN if $H_2 + K_m Q^2 > H_\text{set} + \varepsilon_H$, and OPEN → ACTIVE if $H_1 < H_\text{set} - \varepsilon_H$. From CLOSED, the OPEN test runs first: CLOSED → OPEN if $H_2 > H_\text{set} + \varepsilon_H$ and $H_1 > H_2 + \varepsilon_H$ — when both this and the ACTIVE condition hold, the valve opens; only otherwise CLOSED → ACTIVE if $H_1 \geq H_\text{set} + \varepsilon_H$ and $H_1 > H_2 + \varepsilon_H$. The XPRESSURE state transitions to CLOSED on reverse flow.
• FCV: transitions to XFCV (cannot enforce set point) if the head difference across the valve is negative or flow is negative. A third ACTIVE→XFCV condition also exists: when the valve is active but the pressure drop across it implies a head-loss coefficient smaller than its fully-open minor-loss coefficient (i.e., the network cannot maintain even a friction-free connection without violating the setpoint), the valve also reverts to XFCV. Transitions back to ACTIVE from XFCV once the flow meets or exceeds the setting.

A PRV/PSV/FCV whose status has been fixed OPEN or CLOSED by a control (its setting voided, §7.1) bypasses all of this logic: it is coefficiented as a plain open or closed link and undergoes no status transitions until a new setting is assigned.

Here $H_\text{set}$ is the absolute head setpoint (elevation of the controlled node plus the setting in pressure-head units), $K_m Q^2$ is the minor-loss head drop at the current flow, and $\varepsilon_H$, $\varepsilon_Q$ are user-configured head and flow tolerances. Important distinction: $\varepsilon_H = \text{Htol}$ (default 0.0005 ft) and $\varepsilon_Q = \text{Qtol}$ (default 0.0001 ft³/s) are tolerances used in link status transition tests; $\varepsilon_Q$ appears nowhere else, while $\varepsilon_H$ additionally supplies the dead-band for pressure-based simple controls (§7.1). They are entirely separate from the convergence tolerance $\text{Hacc}$ (default 0.001), which governs solver termination via the relative flow-change criterion (§2.4). Confusing $\text{Qtol}$ with $\text{Hacc}$ leads to incorrect valve and check-valve status behaviour.

Linearisation coefficients for the resistance-type valve types:

Throttle Control Valve (TCV): the minor-loss coefficient is computed from the valve setting $s$ (a dimensionless loss coefficient value) and pipe diameter $D$ (in feet):

$$K_m = \frac{0.02517 , s}{D^4}$$

This factor converts from the user-supplied dimensionless loss coefficient into the internal US customary unit system (flow in ft³/s, head in ft). The coefficient enters the standard minor-loss formula $h = K_m Q |Q|$.

Pressure Breaker Valve (PBV): a PBV imposes a fixed head loss equal to its setting $h_\text{set}$ (in feet) when the setting exceeds the current minor-loss head drop $K_m Q^2$. In this active regime the solver enforces the exact head drop by assigning very large linearisation coefficients: $P_k = C_\infty$ and $Y_k = h_\text{set} \cdot C_\infty$, where $C_\infty$ is a large constant ($\approx 10^8$). Because the GGA flow update is $\Delta Q_k = P_k (H_i - H_j) - Y_k$, this drives $H_i - H_j \to Y_k / P_k = h_\text{set}$ extremely strongly on every iteration. If the current minor-loss head drop already exceeds the setting (the valve is overmatched), the PBV instead falls back to ordinary pipe treatment.

General Purpose Valve (GPV): the solver evaluates the user-supplied head-loss-vs-flow curve at the current absolute flow $|Q|$ to extract the local slope $r$ (ft per ft³/s) and zero-intercept $h_0$ (ft) of the bracketing linear segment. The linearisation coefficients are:

$$P_k = \frac{1}{r}, \qquad Y_k = \left(\frac{h_0}{r} + |Q|\right) \mathrm{sign}(Q)$$

Positional Control Valve (PCV): the percent-open setting $s$ is mapped through a user-supplied opening-to-flow-coefficient ratio curve — below the curve’s first point the ratio interpolates from the origin; above its last point $(x_n, y_n)$ it interpolates along the straight line joining that point to the literal coordinate pair $(1, 1)$ — an anchor evidently intended as $(100%, 100%)$ but written as unity on the curve’s percent-scaled axes, so for $x_n > 1$ the slope is negative and the computed ratio decreases beyond the last point (the actual arithmetic, preserved here as the code’s behaviour); the ratio is clamped to $[10^{-6}, 1]$ and the resulting $K_m$ capped at $10^8$ — to obtain the dimensionless ratio $k_{vr} = K_v / K_{v0}$, where $K_{v0}$ is the flow coefficient at full open. The curve’s $x$-axis is percent open and its $y$-axis is $K_v / K_{v0}$ as a percentage. The effective minor-loss coefficient is then:

$$K_m = \frac{K_{m0}}{k_{vr}^2}$$

where $K_{m0}$ is the fully-open minor-loss coefficient. This reflects the relationship that for a given flow, halving the effective orifice area quadruples the head loss.

Active-state matrix modifications for PRV, PSV, and FCV:

When these valves are in the ACTIVE state, $P_k$ is set to zero and the link is excluded from the standard off-diagonal assembly (which skips links with $P_k = 0$). Instead, the linear system is augmented directly to enforce the valve’s constraint:

PRV active: the downstream node head is pinned to the absolute setpoint $H_\text{set} = z_{n_2} + s_k$ by injecting a large conductance into the diagonal and a correspondingly large forcing term into the RHS:

$$A_{jj} \mathrel{+}= C_\infty, \qquad F_j \mathrel{+}= H_\text{set} \cdot C_\infty$$

The link’s $Y_k$ is set to the current flow plus the downstream node’s flow excess, maintaining approximate flow balance. When the downstream node shows a flow deficit (negative excess), that deficit is added to the upstream node’s RHS to preserve global mass conservation; a positive excess is not transferred.

PSV active: identical treatment applied to the upstream node $i$ with $H_\text{set} = z_{n_1} + s_k$. A small residual conductance ($1/C_\infty$) is also added through the off-diagonal entry to preserve matrix connectivity and avoid numerical singularity.

FCV active: the link is rendered nearly disconnected by setting $P_k = 1/C_\infty \approx 0$. The setpoint flow $Q_\text{set}$ is injected as an external demand at the upstream node and as an external supply at the downstream node, both in the flow-excess array and in the RHS vector:

$$F_i \mathrel{-}= Q_\text{set}, \qquad F_j \mathrel{+}= Q_\text{set}$$

The two sides of the FCV are effectively decoupled; the network is solved as if $Q_\text{set}$ flows through the valve as a prescribed boundary condition, and the resulting head difference across the valve is whatever the network produces with that imposed flow.

Check valve (CVPIPE) status transitions: a check valve is treated as a pipe that is permitted to carry flow only in its positive direction. After each Newton–Raphson iteration, its status is re-evaluated. Let $\Delta h = H_i - H_j$ be the head difference and $Q$ the current flow:

• If $|\Delta h| > H_{\text{tol}}$: CLOSED if $\Delta h < -H_{\text{tol}}$ (reverse head gradient); CLOSED if $Q < -Q_{\text{tol}}$ (reverse flow); otherwise OPEN.
• If $|\Delta h| \leq H_{\text{tol}}$: CLOSED if $Q < -Q_{\text{tol}}$; otherwise the current status is preserved.

This hysteresis prevents rapid cycling near the zero-flow condition.

2.4 The Global Gradient Algorithm

The hydraulic solver at each time step employs the Todini–Pilati Global Gradient Algorithm (GGA), a variant of Newton–Raphson that solves simultaneously for all unknown junction heads and then derives all link flows in a single update.

Governing Equations

Two sets of equations must be satisfied simultaneously.

Flow conservation at each junction $i$:

$$\sum_{k \in \text{in}(i)} Q_k ;-; \sum_{k \in \text{out}(i)} Q_k ;=; D_i$$

where the sums are over all links $k$ whose flow enters or leaves junction $i$, and $D_i$ is the total demand withdrawn at node $i$. The demand includes consumer base demand (scaled by patterns), emitter outflow, leakage, and — in pressure-driven mode — the pressure-dependent portion of consumer demand.

Head-loss equation for each link $k$ connecting nodes $i$ and $j$:

$$H_i - H_j = h_k(Q_k)$$

where $h_k$ is a nonlinear function of $Q_k$ (pipe friction, pump curve, or valve characteristic). For a pump, $h_k < 0$ (head gain).

The network has $n_j$ unknown junction heads and $n_l$ unknown link flows, giving $n_j + n_l$ unknowns and the same number of equations. The GGA exploits the specific algebraic structure to reduce this to a system of size $n_j$ for the heads alone.

Linearisation

At iteration $m$, the head-loss function of link $k$ is linearised around the current flow estimate $Q_k^{(m)}$ by a first-order Taylor expansion:

$$h_k(Q_k) ;\approx; h_k(Q_k^{(m)}) ;+; \frac{\partial h_k}{\partial Q_k}\bigg|_{Q^{(m)}} \left(Q_k - Q_k^{(m)}\right)$$

Two derived quantities characterise every link:

$$P_k = \frac{1}{\displaystyle\frac{\partial h_k}{\partial Q_k}\bigg|_{Q^{(m)}}}, \qquad Y_k = P_k \cdot h_k(Q_k^{(m)})$$

$P_k$ is the inverse of the head-loss gradient and has the dimensions of flow per unit head; it plays the role of a hydraulic conductance. $Y_k$ is the normalised head-loss term, representing the flow contribution from the current head-loss value.

Assembly of the Linear System

Substituting the linearised head-loss relationships into the flow-conservation equations yields a symmetric sparse linear system for the unknown heads:

$$\mathbf{A} , \mathbf{H} = \mathbf{F}$$

The coefficient matrix $\mathbf{A}$ has the structure of a weighted graph Laplacian:

$$A_{ii} = \sum_{k \ni i} P_k \qquad \text{(diagonal: sum over all links incident to node } i\text{)}$$

$$A_{ij} = -P_k \qquad \text{(off-diagonal: } k \text{ is the link connecting } i \text{ and } j\text{)}$$

The right-hand side vector $\mathbf{F}$ at junction $i$ accumulates:

$$F_i = \left(\sum_{k \in \text{out}(i)} Y_k - \sum_{k \in \text{in}(i)} Y_k\right) + \Delta_i$$

where $\Delta_i$ includes the fixed-head boundary contributions from any reservoir or tank directly connected to junction $i$ (their known heads multiply the corresponding $P_k$ and are moved to the right-hand side), plus the demand imbalance at the current iteration. Emitters, leakage, and pressure-dependent demands each add their own linearised conductance to the diagonal of $\mathbf{A}$ and their corresponding $Y$-terms to $\mathbf{F}$.

Flow Update

Once the linear system is solved for the new heads $\mathbf{H}^{(m+1)}$, the flow in each link is updated as:

$$Q_k^{(m+1)} = Q_k^{(m)} - Y_k + P_k \left( H_i^{(m+1)} - H_j^{(m+1)} \right)$$

This is exactly the Newton correction obtained by inverting the linearised head-loss equation $H_i - H_j = h_k(Q_k)$ for $Q_k$: the current-flow term $Q_k^{(m)}$ and the offset $-Y_k = -P_k h_k(Q_k^{(m)})$ carry the old operating point, while $P_k(H_i - H_j)$ applies the head-driven correction. Here $i$ is the upstream node and $j$ the downstream node of link $k$ according to the assumed positive direction. For pumps, $H_j - H_i = \Delta H_k > 0$, so the sign convention is consistent.

Emitter flows, leakage flows, and pressure-dependent demand flows are similarly updated using their own linearised head-flow relationships.

Convergence Criterion

The primary convergence criterion is flow accuracy: the sum of absolute flow changes relative to total absolute flow — both sums accumulated over all links and all nodal elements (emitter, pressure-dependent-demand, and leakage flows) — must fall below the user-specified tolerance, and no link’s hydraulic status may have changed during the most recent iteration:

$$\epsilon = \frac{\displaystyle\sum_k \left| Q_k^{(m+1)} - Q_k^{(m)} \right|}{\displaystyle\sum_k \left| Q_k^{(m+1)} \right|} \leq \epsilon_{\text{tol}}$$

When the total absolute flow $\sum_k |Q_k^{(m+1)}|$ falls at or below $\epsilon_{\text{tol}}$ (near-stagnant network), the relative formula cannot be used; in that case the solver returns the absolute flow change $\sum_k |\Delta Q_k|$ directly rather than the ratio.

Two kinds of status check run during iteration, with different scheduling:

• Valve status checks (valvestatus, governing PRV and PSV transitions): when DampLimit = 0 (the default), these run after every flow update. When DampLimit > 0, they are deferred until the relative flow error at or below DampLimit; at that point damping (relaxation factor 0.6) is also activated — the check occurs after the current iteration’s flow update, so the damping takes effect on the next iteration’s update.
• Link status checks (linkstatus, governing pumps, check valves, FCVs, and pipes adjacent to tanks): these run periodically. The first check occurs at iteration CheckFreq and repeats every CheckFreq iterations thereafter, but stops once MaxCheck iterations have been reached. This staging prevents premature status oscillation during early iterations when flows are far from convergence.

When hasconverged returns true, linkstatus and pswitch are called unconditionally (regardless of the CheckFreq schedule), alongside the valvestatus result already computed for that iteration under the DampLimit schedule. If any of them changes a link’s status the periodic-check schedule (nextcheck) is advanced and the solve continues — the iteration counter itself is never reset; it only increments toward the iteration cap. Only a convergence pass that produces no status changes terminates the loop.

The ExtraIter parameter handles networks that fail to converge because of status cycling — a cycle in which links repeatedly toggle between open and closed states without settling to a consistent configuration. When convergence is not achieved within MaxIter iterations and ExtraIter > 0, an additional ExtraIter iterations are performed during which a convergence pass exits the loop immediately, skipping the convergence-time linkstatus and pswitch calls — pumps, check valves, FCVs, and tank-adjacent pipes no longer change state at convergence. The periodic linkstatus checks are also silent during the extra iterations, but only because MaxCheck (default 10) lies far below MaxIter (default 200) — an emergent consequence of the defaults rather than an explicit suspension. valvestatus (PRV/PSV) continues to run every iteration. This allows the linear system to converge to a solution consistent with the current link configuration, even if that configuration is not the true steady state. A warning is issued that the system is unbalanced. If ExtraIter = −1, the time-step routine wrapping the solver sets a halt flag (Haltflag) as soon as the unbalanced solve returns — before the step’s results are saved; the results are then saved as usual, and the simulation terminates at the start of the next step rather than stopping mid-step.

Two supplementary convergence criteria may also be applied; rather than terminating the loop early, they make convergence stricter — each is tested only after the flow-accuracy criterion has already passed, and can withhold convergence for another iteration. FlowChangeLimit requires the maximum absolute flow change over all links and nodal elements (emitters, pressure-dependent demands, leakage) during the most recent iteration to be at or below its threshold. HeadErrorLimit requires the largest discrepancy, over all open links, between the head difference across the link and the head loss implied by its current flow to be at or below its threshold — a per-link head-loss residual, not a nodal quantity. Both default to zero, which disables them. Two further gates are always active when their models are in use: convergence additionally requires the leakage flows and — in pressure-driven mode — the pressure-dependent demands to have themselves converged. In the default configuration (fixed demands, no leakage) the sole termination criterion is $\epsilon \leq \epsilon_{\text{tol}}$ with no status change.

Damping (RelaxFactor): the Newton flow update $\Delta Q_k$ may be scaled by a relaxation factor to improve convergence stability. By default RelaxFactor = 1.0 (full Newton step). When DampLimit > 0 and the relative flow error falls at or below DampLimit, RelaxFactor is set to 0.6; because the factor is set after the current iteration’s flow update, it damps the following iteration’s updates. This under-relaxation is applied uniformly to all link flow updates, emitter flow updates, pressure-dependent demand flow updates, and leakage flow updates:

$$Q_k^{(m+1)} = Q_k^{(m)} - \text{RelaxFactor} \cdot \Delta Q_k$$

The purpose is to stabilise convergence in networks with highly nonlinear elements (e.g., active control valves) by reducing the step size when the solver is close to convergence but oscillating.

Matrix recovery (badvalve): if the Cholesky factorisation of $\mathbf{A}$ fails at a diagonal entry corresponding to node $n$, the solver checks whether an active PRV, PSV, or FCV has that node as one of its endpoints. If found, the valve’s status is forced to XPRESSURE (for PRV/PSV) or XFCV (for FCV), breaking the singularity that the active-state matrix modification introduced. The solver then retries the factorisation and solve. This recovery mechanism ensures that ill-conditioned valve configurations do not crash the solver. The recovery examines only the first valve found at the offending node; if that one is not active, the solve fails with an ill-conditioning error even when another valve there is active.

The XFLOW status (“pump exceeds maximum flow”) is not produced by the solver’s iterative status logic — the pump status evaluation checks only whether the head gain exceeds the speed-adjusted shutoff head (yielding XHEAD), and does not compare flow against the maximum-flow point of the pump curve. It is instead raised diagnostically after a step has been solved: the warning routine (writehydwarn) flags a pump as XFLOW when its flow exceeds $\text{setting} \cdot Q_{\max}$ (emitting warning WARN04), and the same status is returned through the public link-status query. XFLOW therefore never alters the solve; it is a reported warning state only.

2.5 Sparse Linear Algebra

The matrix $\mathbf{A}$ is symmetric and positive semi-definite (it is a Laplacian, plus small positive diagonal contributions from emitters and demands). Its sparsity pattern corresponds to the adjacency structure of the junction subgraph. Efficient solution is essential because this system must be solved at every Newton iteration of every hydraulic time step.

Three phases are performed:

Phase 1 — Node reordering (performed once before simulation begins): the Multiple Minimum Degree (MMD) algorithm reorders the junction indices to minimise the fill-in that occurs during Cholesky factorisation. Fill-in arises when a non-zero appears in the factor $\mathbf{L}$ at a position that was zero in $\mathbf{A}$; reordering the rows and columns can dramatically reduce the number of such positions. Parallel links (multiple links connecting the same pair of nodes, regardless of their orientation) are condensed into a single equivalent link before reordering to avoid redundancy. (The reordering routine is invoked with its multiple-elimination tolerance set to −1, so in practice one minimum-external-degree node is eliminated per update cycle — plain minimum external degree rather than true multiple elimination.)

Phase 2 — Symbolic factorisation (performed once before simulation begins): using the reordered sparsity pattern, the algorithm predetermines the exact set of non-zero positions that will appear in the lower Cholesky factor $\mathbf{L}$ (satisfying $\mathbf{A} = \mathbf{L} \mathbf{L}^\top$). These positions are stored in a compressed sparse form. From this point forward, only numerical values need to change; the structure is fixed.

Phase 3 — Numerical factorisation and solution (performed at every Newton iteration): the current values of $P_k$ are inserted into the pre-allocated arrays, the Cholesky factorisation is carried out in the stored non-zero positions, and forward and backward substitution yields the updated head vector $\mathbf{H}$.

Because the sparsity structure does not change between iterations (only the values change), Phases 1 and 2 are not repeated, making per-iteration cost proportional only to the number of non-zeros in $\mathbf{L}$.

3. Demand Models

3.1 Demand-Driven Analysis

In the default Demand-Driven Analysis (DDA) mode, all demands are treated as fixed withdrawals regardless of the pressure at the node. Each junction has one or more demand categories; within each category a base demand rate is multiplied by the current value of its associated time pattern to give the instantaneous withdrawal rate. Multiple categories are summed. A global demand multiplier ($D_\text{mult}$) is also applied uniformly to all base demands at all junctions, allowing the overall demand level to be scaled up or down without modifying individual data. If the net demand at a node is negative, the node acts as an inflow point (external source).

This model is simple and numerically robust, but it can produce physically unrealistic results for heavily stressed systems: a node with insufficient pressure will still show its full demand satisfied, possibly at a negative computed pressure.

3.2 Pressure-Driven Analysis

In Pressure-Driven Analysis (PDA) mode, the demand actually delivered depends on the available pressure at each node. The governing relationship is:

$$D(P) = \begin{cases} 0 & P \leq P_{\min} \ D_{\text{full}} \left( \dfrac{P - P_{\min}}{P_{\text{req}} - P_{\min}} \right)^{n_P} & P_{\min} < P < P_{\text{req}} \ D_{\text{full}} & P \geq P_{\text{req}} \end{cases}$$

where $P$ is the gauge pressure at the node, $P_{\min}$ is the pressure below which no demand is delivered, $P_{\text{req}}$ is the pressure at which full demand is delivered, $D_{\text{full}}$ is the requested demand, and $n_P$ is the pressure exponent. The default value of $n_P = 0.5$ corresponds to the Wagner formula, which models demand as proportional to the square root of the available pressure head above the minimum threshold. The three parameters are global — one set for the whole network — defaulting to $P_{\min} = 0$, $P_{\text{req}} = 0.1$, $n_P = 0.5$, with $P_{\text{req}} - P_{\min} \geq 0.1$ psi (or m) enforced. Pressure-dependent treatment applies only to junctions whose pattern-adjusted total demand is positive: zero- and negative-demand junctions (external inflow points) keep their fixed demand-driven values even in PDA mode.

The PDA model is incorporated into the GGA by treating the pressure-dependent component of demand as a pressure-dependent emitter at each junction (cf. §4). The demand-pressure function is inverted to express pressure as a function of demand:

$$P = P_{\min} + (P_{\text{req}} - P_{\min}) \left( \frac{D}{D_{\text{full}}} \right)^{1/n_P}$$

This is linearised and added to the diagonal of $\mathbf{A}$ and the right-hand side $\mathbf{F}$. Barrier terms prevent the numerical demand from drifting below zero or above $D_{\text{full}}$, maintaining the physical bounds throughout the iteration. The barrier is implemented as a smooth differentiable approximation (not a hard constraint) to avoid discontinuities that would break the Newton–Raphson iteration. Specifically, the signed head-loss and gradient increments from a lower barrier at $Q = 0$ take the form:

$$\Delta h = \frac{a - \sqrt{a^2 + 10^{-6}}}{2}, \qquad \Delta(\partial h/\partial Q) = \frac{10^9}{2}\left(1 - \frac{a}{\sqrt{a^2 + 10^{-6}}}\right), \qquad a = 10^9 , Q$$

which approaches a large one-sided penalty as $Q \to 0^-$ while remaining smooth and differentiable throughout. An analogous upper barrier is applied at $Q = D_{\text{full}}$. Each iteration’s change in a node’s pressure-dependent demand is additionally capped at 40% of its full demand, preventing overshoot of the bounded demand function. PDA convergence is tested by re-evaluating the demand function at the newly computed pressures and requiring agreement with the solved demand flows within $10^{-4}$ ft³/s at every node; the same pass counts pressure-deficient junctions and the total percent demand shortfall reported to the user.

4. Emitters

An emitter represents a device — a sprinkler head, orifice, or nozzle — that discharges water from a junction at a rate governed by the local pressure. The emitter flow is:

$$Q_e = K_e , P^{n_e}$$

where $K_e$ is the per-junction emitter discharge coefficient, $P$ is the gauge pressure at the junction, and $n_e$ is the pressure exponent — a single network-wide option (default 0.5, the orifice value); only the coefficient varies per junction. In US units the coefficient is defined per psi$^{n_e}$ (adjusted for specific gravity), in SI per metre of head. Emitters can also represent aggregate leakage if high spatial resolution is not required.

Within the GGA, an emitter is treated as an additional element at its junction. The head-flow relationship is inverted:

$$H - z = C_e , Q_e^{1/n_e}$$

where $H - z$ is the pressure head and $C_e = K_e^{-1/n_e}$. This is linearised at the current flow estimate and added to the matrix: the linearised conductance $P_e = 1 / (\partial h_e / \partial Q_e)$ is added to the diagonal of $\mathbf{A}$ at the relevant junction, and the corresponding $Y_e$ term is added to the right-hand side. When the head-loss gradient falls below a resistance tolerance ($10^{-7}$), the emitter relation is replaced by a linear one to keep the Jacobian bounded near zero flow.

By default emitters can admit reverse flow (suction) if the junction pressure falls below the emitter’s reference elevation. An emitter backflow option can disable this: when backflow is forbidden, a lower barrier function is applied to the head-loss gradient whenever the emitter flow would go negative. The barrier takes the same smooth differentiable form used for PDA demand bounds (§3.2):

$$\Delta h = \frac{a - \sqrt{a^2 + 10^{-6}}}{2}, \qquad \Delta(\partial h/\partial Q) = \frac{10^9}{2}\left(1 - \frac{a}{\sqrt{a^2 + 10^{-6}}}\right), \qquad a = 10^9 , Q_e$$

This adds a large one-sided penalty as $Q_e \to 0^-$ while remaining smooth and differentiable, strongly driving $Q_e \geq 0$ without creating a hard discontinuity that would break convergence.

5. Pipe Leakage — the FAVAD Model

Background leakage from deteriorated pipes — through corroded joints, stress cracks, and micro-fractures — is modelled using the FAVAD (Fixed And Variable Area Discharge) framework. Unlike a simple orifice, pipe cracks may dilate under pressure, making the effective discharge area itself pressure-dependent. The FAVAD model captures this through:

$$Q_{\text{leak}} = C_o , \frac{L}{100} \left( A_o + m H \right) \sqrt{H}$$

where $H$ is the pressure head driving leakage from the pipe, $L$ is the pipe length, $A_o$ is the fixed (zero-pressure) crack area per 100 units of pipe length, $m$ is the rate of increase of crack area with pressure head (same per-100-length basis), and $C_o = 0.6\sqrt{2g} \approx 4.815$ in pure ft units — leakage scales linearly with pipe length. The crack parameters are supplied in mm² (and mm² per metre of head): internally the coefficient absorbs the mm²→m² factor $10^{-6}$ and a further division by $0.3048^2$ (fixed-area term) or $0.3048$ (variable-area term), so with areas in their input units the effective fixed-area coefficient is $\approx 5.18 \times 10^{-5}$. Expanding this:

$$Q_{\text{leak}} = C_o \tfrac{L}{100} A_o H^{1/2} + C_o \tfrac{L}{100} m H^{3/2}$$

The two terms have different pressure exponents: the fixed-area term behaves like a standard orifice (exponent $1/2$), while the variable-area term has exponent $3/2$.

These are decomposed into two equivalent emitters at each node, one for each component:

$$H = C_{\text{fa}} , Q_{\text{fa}}^{2} \qquad \text{(fixed-area component, orifice-type, exponent } 1/2 \text{ on } H\text{)}$$

$$H = C_{\text{va}} , Q_{\text{va}}^{2/3} \qquad \text{(variable-area component, exponent } 3/2 \text{ on } H\text{)}$$

The resistance coefficients $C_{\text{fa}}$ and $C_{\text{va}}$ are determined from the FAVAD parameters. For each pipe whose both end nodes are junctions, the pipe’s leakage contribution is split equally: half is attributed to each end node (the pipe is split conceptually at its midpoint, with each half’s leakage then driven by that end node’s own pressure head rather than a single midpoint value). When one end of a pipe is a fixed-grade node (reservoir or tank), that fixed-grade end cannot accumulate leakage in the nodal model; the junction at the other end therefore receives the full pipe-length contribution rather than half. A pipe whose both end nodes are fixed-grade contributes no leakage at all — leakage is realised as a nodal demand, which fixed-grade nodes carry none of. The contributions of all pipes meeting at a given junction are aggregated: the total fixed-area conductance and variable-area conductance at the node are the sums over all incident (half- or full-) pipe contributions. The resulting nodal coefficients are then inverted to form $C_{\text{fa}}$ and $C_{\text{va}}$.

Derivation of $C_{\text{fa}}$ and $C_{\text{va}}$: let $\text{LeakCoeff1}_p = C_o A_o L_p / 100$ and $\text{LeakCoeff2}_p = C_o, m, L_p / 100$ be the full-pipe FAVAD discharge coefficients for pipe $p$ (unit conversions absorbed as above). The per-end contribution for a junction endpoint $v$ of pipe $p$ is:

$$k_{1,p,v} = \begin{cases} \tfrac{1}{2},\text{LeakCoeff1}_p & \text{both end nodes of pipe } p \text{ are junctions} \ \text{LeakCoeff1}_p & \text{exactly one end node of pipe } p \text{ is a fixed-grade node} \end{cases}$$

(with the same rule applied to $\text{LeakCoeff2}p$ to give $k{2,p,v}$). For junction $i$, the total discharge conductances are $K_{\text{fa},i} = \sum_{p \ni i} k_{1,p,i}$ and $K_{\text{va},i} = \sum_{p \ni i} k_{2,p,i}$. The resistance coefficients follow by inverting the discharge relations $Q = K H^{1/2}$ (fixed-area) and $Q = K H^{3/2}$ (variable-area):

$$C_{\text{fa},i} = 1/K_{\text{fa},i}^{2}, \qquad C_{\text{va},i} = 1/K_{\text{va},i}^{2/3}$$

(with the respective term omitted when $K = 0$).

These two emitter-like terms are linearised and incorporated into the GGA matrix assembly in exactly the same way as ordinary emitters (§4): a conductance term is added to the diagonal of $\mathbf{A}$ and a flow offset term is added to the right-hand side $\mathbf{F}$.

Leakage is non-negative by construction: the smooth lower barrier of §3.2 is applied unconditionally to both leakage components — unlike emitters, where it is conditional on the backflow option — so a node at sub-atmospheric pressure simply leaks nothing; the post-solution per-pipe leakage report likewise clamps negative pressure heads to zero. Leakage carries its own convergence gate: each node’s solved leakage must match a direct FAVAD evaluation at the current pressure, $\sqrt{H/C_{\text{fa}}} + (H/C_{\text{va}})^{3/2}$, within $10^{-4}$ ft³/s; leakage flow corrections are damped by the same relaxation factor as the flow update and enter the global relative-flow-change measure (§2.4). Before iteration begins, the nodal elements are seeded with nonzero flows — 1.0 ft³/s per emitter, 0.001 ft³/s per leakage component, and the full demand for pressure-dependent demands — and at the converged solution the reported junction demand aggregates all three outflow elements: consumer demand plus emitter flow plus leakage flow.

6. Time-Stepping and Tank Dynamics

6.1 Extended-Period Simulation

The extended-period simulation (EPS) advances the network state through a sequence of discrete hydraulic time steps of duration $\Delta t$. Within each hydraulic step the network configuration (demands, pump settings, valve statuses) is assumed constant, and a steady-state hydraulic solution is computed. The procedure at each step is:

1. Apply patterns (demands): evaluate all time patterns at the current clock time and apply the resulting multipliers to demands at junctions, reservoir heads, and pump speed settings.
2. Simple controls (controls): evaluate and apply all simple controls that trigger at this time (§7.1). May change link status or settings.
3. Solve (hydsolve): solve for hydraulic equilibrium via the Global Gradient Algorithm (§2.4).
4. Compute pump power (getallpumpsenergy): determine the current power draw and efficiency at each pump using the just-solved flow field.
5. Determine time step (timestep): compute the next time-step duration as the minimum of the nominal step, reporting interval, pattern change, tank fill/drain time, and rule evaluation (§6.2). Rule-based controls (§7.2) are evaluated within this step at intermediate rule-step intervals; if a rule fires, the time step is shortened to the firing time. Tank levels are updated within this computation (not after it).
6. Accumulate energy (addenergy): accumulate pump energy consumption (kWh) and cost over the time step. Track peak demand.
7. Update flow balance (updateflowbalance): accumulate volumetric flow balance ledger entries for the step.
8. Advance clock: advance the simulation time by the computed time-step duration.

Steps 1–8 repeat until the specified simulation duration is reached.

Initial hydraulic state: PRV/PSV/FCV valves carrying a numeric setting begin the simulation ACTIVE; any non-GPV valve initialised to a fixed OPEN/CLOSED status has its setting wiped (no setpoint is enforced until one is assigned); junctions with emitters start from an emitter-flow guess of 1.0 ft³/s; tank pseudo-links start TEMPCLOSED; and on re-initialisation existing link flows are preserved as a warm start unless re-initialisation is explicitly requested or the flow is near zero.

Important: rule-based controls are evaluated after the hydraulic solve, within the time-step computation (step 5), not before or during the solve. If a rule fires, the hydraulic step is shortened so that the next solve will reflect the new configuration — the current step’s solution is not re-computed. Tank levels are updated by ruletimestep() at each rule sub-step interval; when no rules exist, tanklevels() is called directly during timestep().

6.2 Adaptive Time Step

The hydraulic time step is not fixed; it adapts so that no physical limit is overshot. The actual step duration used is the minimum of the following quantities:

• The user-specified nominal hydraulic time step.
• The time remaining until the next reporting interval (so that results are recorded at exactly the right moments).
• For each tank, the time at which the tank would reach its minimum level (if the net outflow continues at the current rate) or its maximum level (if the net inflow continues): $\Delta t_{\text{tank}} = \Delta V_{\text{available}} / |Q_{\text{net}}|$.
• The time until the next scheduled change in any pattern (so that demand or pump multipliers change at exactly the right instant).
• The time until the next simple control activates (the earliest instant at which a timer control fires or a tank-level control crosses its threshold at the current fill/drain rate) — considering only controls whose action would actually change the link’s current status or setting; level controls on junction pressures never participate, and timer controls fire only on exact clock coincidence, so a timer the clock steps past is skipped entirely.
• When rule-based controls are present, the time at which the next rule fires — determined by sub-stepping the rule evaluation across the hydraulic period (see §6.1 and §7.2).

The end of the simulation is not one of these minimised quantities; it is enforced separately, by the main loop simply ceasing to issue steps once the clock reaches the specified duration. This adaptive strategy avoids the need for post-hoc correction of tank levels and ensures pattern changes and control actions are applied at their intended times.

Default time step derivation: if the user does not specify either the quality time step or the rule evaluation time step, both default to $\Delta t_h / 10$, capped at $\Delta t_h$. The rule time step is additionally aligned so that evaluations fall on even multiples of the rule step within each hydraulic period — the first evaluation within a period may therefore be shorter than one full rule step to achieve this alignment. The quality time step is further constrained so it never exceeds the hydraulic time step. The hydraulic time step itself is clamped at the minimum of the user-specified nominal step, the pattern time step, and the reporting time step.

6.3 Tank Level Update

After each hydraulic solution, the change in stored volume during the time step is:

$$\Delta V = Q_{\text{net}} \cdot \Delta t$$

where $Q_{\text{net}} = Q_{\text{in}} - Q_{\text{out}}$ is the net volumetric flow rate into the tank. For a constant cross-section tank with cross-sectional area $A$, the level change is:

$$\Delta h = \frac{\Delta V}{A}$$

For a tank described by a volume-elevation curve, the new volume $V_{\text{new}} = V_{\text{old}} + \Delta V$ is looked up in the curve to find the corresponding new water surface elevation.

The post-step limit handling is a predictive one-second snap, compensating for the fill/drain time being rounded to the nearest whole second: a tank within one second’s net inflow of its maximum volume is set to exactly $V_{\max}$ (the level need never actually cross the limit), and one that has fallen a second’s net outflow past its minimum is set to $V_{\min}$; the whole update, head recomputation included, is skipped when $|Q_{\text{net}}| \leq 10^{-6}$ ft³/s. A tank at its minimum is treated as a fixed-grade node (like a small reservoir at its minimum head) for the next time step. The behavioural difference between the two overflow modes is enforced during the hydraulic solve rather than by this snap. When overflow is allowed, an over-full tank keeps accepting inflow and the surplus exits freely as overflow while the level is held at the maximum; when overflow is not permitted, the inflow links are held closed (see TEMPCLOSED below) so the tank cannot exceed its maximum, and it is treated as fixed-grade at that level.

TEMPCLOSED for links at tank limits: independently of the level-clamping applied after the hydraulic step, whenever the solver performs a link-status check (at the periodic CheckFreq interval while the iteration count remains within MaxCheck, and whenever a trial solution first satisfies the convergence criteria) the status of every link adjacent to a tank is examined. If a link is carrying flow into a tank whose current head equals or exceeds the maximum level and overflow is not permitted, that link is set to TEMPCLOSED — it retains only the vanishingly small conductance assigned to closed links, and the coefficient matrix is assembled as if the link were closed. Similarly, a link carrying flow out of a tank at its minimum level is TEMPCLOSED. At the start of each such status check, all TEMPCLOSED and XHEAD links are first re-opened before the new operating point is evaluated, so the status is re-tested rather than being locked in.

7. Control Systems

7.1 Simple Controls

A simple control consists of a single condition and a single action. The condition is either a level control — a node’s pressure or hydraulic grade exceeds or falls below a specified threshold — or a timer control — the simulation clock reaches a specified time or the current time of day reaches a specified hour.

When a simple control fires, its action is applied immediately: it may open or close a link, change a pump’s speed setting, or change a valve’s setting. Status and setting actions on valves are mutually destructive: opening or closing a valve discards its setting (an “open” action converts an active PRV/PSV/FCV into a plain open link), assigning a setting to an FCV forces it ACTIVE, and assigning a setting to a closed valve whose setting was wiped re-opens it; a pump opened by status takes speed 1.0 and closed takes 0.0, and a PCV’s loss coefficient is recomputed the moment its setting changes. Tank-level and time-based controls are evaluated once at the start of each hydraulic time step; if such a control fires and changes the network configuration, the subsequent hydraulic solution reflects the new state for the entire duration of that time step. Simple controls conditioned on a junction pressure or grade are handled differently: because the controlling pressure is itself an output of the solve, they are re-checked after the hydraulic system has converged, by the same pswitch routine that runs inside the solver loop. If the converged pressure has crossed the control threshold, the link status is switched and the system is re-solved, iterating until no such control changes state.

For level controls on tanks, a small hysteresis margin is applied to prevent chattering when the tank is exactly at the control threshold under non-zero flow. The trigger condition is checked against the tank volume corresponding to the control’s grade level, with a margin equal to the current absolute value of the tank’s net demand flow rate (|NodeDemand|, in internal flow units). A low-level control fires when the current tank volume falls at or below the threshold volume plus the margin; a high-level control fires when the current tank volume reaches or exceeds the threshold volume minus the margin.

Actions are normalised at parse time: pump OPEN/CLOSED map to speeds 1/0, a numeric pump or pipe setting maps to open/closed status, GPVs accept only OPEN/CLOSED, check-valve pipes cannot be controlled at all, and clock-time triggers wrap modulo 24 hours. Both simple controls and rules may carry a DISABLED tag removing them from evaluation — though the post-convergence junction-pressure pass does not honour the flag, so a disabled pressure-conditioned control still fires there.

7.2 Rule-Based Controls

Rule-based controls support arbitrarily complex conditional logic and are well suited to representing operational strategies such as “turn on pump A if tank X level falls below 5 m AND time of day is between 22:00 and 06:00.”

A rule has the structure:

IF (premise₁) AND/OR (premise₂) … THEN (action₁, action₂, …) ELSE (action₃, action₄, …)

Premises may test:

• Node pressure, hydraulic grade, level (head above the node’s elevation), or demand
• Total system demand
• Link flow rate, status, or setting (a POWER keyword exists in the grammar’s vocabulary, but no premise using it can actually be created)
• Tank fill time (time for the tank to reach its maximum at the current inflow rate) or drain time — each false unless the tank is actually filling (respectively draining)
• Simulation time or time of day

Premises evaluate left-to-right: an OR clause is consulted only when the accumulated result is false, and subsequent AND clauses must still hold — approximating AND binding tighter than OR, except that a false result reaching an AND clause rejects the rule outright, so any later OR alternatives are never consulted; no parentheses are available. Numeric comparisons use a fixed tolerance of 0.001 in user units (with ≤/≥ shifted by it); flow premises test the absolute value of flow; status premises map internal XHEAD/TEMPCLOSED states to CLOSED and admit only IS/NOT; and a setting premise on a valve whose setting has been voided is always false. Clock-time equality premises are satisfied if the stated instant fell anywhere within the elapsed rule interval — including across midnight — while inequalities compare only the interval’s end.

Rules are evaluated not just at the start of a hydraulic time step but at intermediate rule time steps that subdivide each hydraulic step. When a rule fires and changes the network state, the current hydraulic step is shortened to end at the firing time; consistent with §6.1, the shortened step’s solution is not recomputed — the next hydraulic step is solved afresh from the new configuration. This allows the simulation to capture, for example, a pump switching on mid-step in response to a falling tank level.

When multiple rules fire simultaneously and their THEN actions conflict (e.g., two rules disagree on the status of the same pump), priority levels resolve the conflict: the rule with the strictly higher priority value wins — on a tie, the rule evaluated first keeps its action. Rule actions are normalised like simple-control actions (check-valve pipes cannot be acted on, GPVs take no numeric setting, a numeric pipe setting converts to open/closed status), and a setting action is a no-op unless it changes the setting by more than 0.001 — which also determines whether it counts as a state change that shortens the hydraulic step.

8. Water Quality Simulation

8.1 Overview and Simulation Modes

Water quality simulation is layered on top of the hydraulic solution. The flows and velocities computed during the hydraulic phase are stored and replayed during quality simulation, which advances in quality time steps that are no longer than hydraulic time steps, and are typically shorter, in order to resolve the advection of concentration fronts through pipes. The quality time step $\delta t$ is a user-specified parameter (defaulting to one-tenth of the hydraulic step), typically a small fraction of the hydraulic time step to satisfy the Courant stability condition for the advection scheme. Within each hydraulic period of duration $\Delta t_h$, the quality simulation executes multiple transport sub-steps of size $\min(\delta t, \text{remaining time})$, iterating until the full hydraulic period is consumed; within each sub-step reactions are applied first and advection second — an operator split whose ordering matters at coarse $\delta t$ — with the reaction pass skipped entirely when no pipe carries a nonzero bulk or wall coefficient and no tank a nonzero bulk coefficient. The hydraulic flow field (velocities, directions) is held constant throughout the sub-cycles. Any change of a link’s direction class (positive, negative, or negligible — flow below Q_STAGNANT classes as zero) between consecutive hydraulic periods triggers a re-sort of the nodes into topological order, and a link whose flow reverses sign additionally has its segment list reversed in place so the leading end stays downstream.

Three simulation modes are available:

• Chemical concentration: the transport and reaction of a dissolved constituent (e.g., residual chlorine) through the network.
• Water age: the average time water has been resident in the distribution system, measured from the sources. No external source or reaction is needed; the age field is initialised from each node’s initial quality $C_0$ (zero by default — nonzero initial ages are legal) and increases at a uniform rate as time passes.
• Source tracing: the percentage of water at any point in the network that originated from a designated source node. The tracer is injected at 100% at the source; all other inflowing water carries 0%.

8.2 Lagrangian Segment Transport

Advective transport of water quality through pipes is modelled with a Lagrangian moving-segment scheme. Each pipe is represented as an ordered sequence of segments. Each segment has a volume and a uniform concentration. Segments are created when new water of a different concentration enters a pipe and destroyed when they are fully flushed out the other end. At startup each pipe holds one full-volume segment at its downstream node’s initial quality; pumps and valves start with no segments; each tank starts as one segment at $(C_0, V_0)$ — split for two-compartment tanks so the mixing zone fills first.

Advection step: over a quality time step of duration $\delta t$, a volume $\mathcal{V}_k = Q_k , \delta t$ of water is swept through pipe $k$. Consumption proceeds from the downstream (leading) end of the pipe — the segments nearest the receiving node — inward, one segment at a time. The mass and volume of each consumed fraction are tracked; when the cumulative volume consumed equals $\mathcal{V}_k$, the remaining portion of the last partially-consumed segment is retained as the new leading (downstream-end) segment. The total mass and volume that exited the pipe’s downstream end are accumulated at the downstream node.

Nodal mixing: transport is organised as a single sweep over the nodes in topological order, not as separate all-pipes phases. When a node’s turn comes, the advection step above is applied to each pipe carrying flow into it, and its outflow concentration is computed as:

$$c_{\text{out},i} = \frac{\displaystyle\sum_{k \in \text{in}(i)} m_k^{\text{out}}}{\displaystyle\sum_{k \in \text{in}(i)} \mathcal{V}_k^{\text{out}}}$$

where $m_k^{\text{out}}$ and $\mathcal{V}k^{\text{out}}$ are the mass and volume that flowed out of (i.e., into node $i$ from) pipe $k$ during the quality step. This is the complete instantaneous mixing assumption: water entering a junction from all sources is instantly and uniformly blended. A junction with net negative demand (an external inflow point) adds that inflow volume at zero concentration to the denominator, diluting the blend unless a concentration source supplies its mass. The resulting concentration $c{\text{out},i}$ — after tank mixing and any source contribution — is then pushed immediately, before the next node is processed, as a new upstream (trailing) segment into each pipe carrying flow away from node $i$. Because each pipe’s upstream node is swept before its downstream node, this push occurs before that pipe’s own consumption: when $Q_k,\delta t$ exceeds the pipe’s stored volume, the downstream node draws water that entered the pipe during the same sub-step.

Topological ordering: for the advection step to be computed correctly — that is, for each node’s outflow concentration to reflect all its upstream contributions — the nodes must be processed in topological order (upstream nodes before downstream nodes). The network is sorted into such an order prior to the quality simulation. For looped networks where no pure topological sort exists, nodes that share mutual dependencies are processed with a fallback ordering strategy.

Stagnant junctions: when a junction receives no inflow at all during a quality time step — neither from links nor as external inflow (a dead-end or temporarily stagnant condition) — and the constituent is reactive, the junction concentration is updated to the arithmetic mean of the concentrations in the segment nearest to the junction on each incident link. For a pipe whose flow runs into the node (the node is the pipe’s downstream end), this is the pipe’s downstream-end segment (FirstSeg); for a pipe whose flow runs away from the node (the node is the pipe’s upstream end), this is the pipe’s upstream-end segment (LastSeg). In both cases the chosen segment is the one physically adjacent to the stagnant junction, regardless of flow direction. This heuristic prevents the unphysical drift in concentration that would otherwise occur at stagnant junctions driven purely by reaction kinetics without advective renewal.

Segment merging: after the outflow concentration of a node is computed, a new segment is pushed into each pipe carrying flow away from the node. Before a new segment is created, the node’s outflow concentration $c_\text{node}$ is compared with the concentration of the segment already sitting at the upstream end of that outflow pipe (the pipe end adjacent to the node). If the difference is less than the threshold $C_\text{tol}$ (default 0.01 — mg/L for chemicals, hours for age), the existing segment’s volume and mass are merged with the new contribution rather than creating a separate segment; the same tolerance governs inflow-segment coalescing in FIFO and LIFO tanks. Without this merging step, each transport sub-step could create a new segment boundary, causing the total segment count to grow without bound in steady flow. Merging keeps the count manageable at the cost of negligible concentration smoothing.

Segment memory management: pipe segments are allocated from an in-memory pool that grows on demand in large blocks, backed by a free-list of recycled records. When a segment is fully flushed out of a pipe it is returned to the free list rather than released to the operating system. New segments are drawn from the free list first; fresh pool entries are used only when the free list is empty. If the pool cannot be extended, an out-of-memory flag is raised and quality simulation is terminated gracefully with an error.

8.3 Source Terms

Water quality sources inject constituent into the network at designated nodes. Four injection types are available:

• Concentration source: the concentration of all water leaving the node is fixed at the specified value. At reservoirs the source replaces the outflow quality; at tanks it is added on top of the tank’s mixed outflow quality, leaving the stored quality unmodified. At junctions this type is effective only when the node has a net negative demand (i.e., the junction is itself a local inflow point, NodeDemand < 0); if the junction demand is non-negative the source contributes nothing. Inflows from other links are not overridden at junctions.
• Mass inflow booster: a fixed mass rate is added to the node continuously, regardless of flow conditions. The resulting concentration increment depends on the total flow through the node.
• Setpoint booster: if the naturally mixed concentration at the node falls below the specified setpoint, it is raised to the setpoint. If it already exceeds the setpoint, no adjustment is made.
• Flow-paced booster: a fixed concentration increment is added to the natural concentration of all water leaving the node, in proportion to the flow.

Source concentrations may vary over time via a multiplier pattern. Mass-source strengths are interpreted as mass per minute (converted internally to per-second); the other source types are concentrations in user units. A source whose baseline strength is zero is inert regardless of its pattern.

For all source types, a stagnation guard suppresses injection when the total volumetric outflow rate from the node during a quality sub-step falls at or below a small threshold Q_STAGNANT $= 1.114 \times 10^{-5}$ ft³/s (0.005 gpm). This avoids division by zero — and unstable large increments — when computing concentration increments at a near-stagnant node. The same Q_STAGNANT threshold also governs whether a link’s flow is treated as negligible for the topological-sort and transport steps (§8.1); it is one constant serving both purposes, not two. (It is distinct from the hydraulic constant QZERO $= 10^{-6}$ ft³/s, which seeds the flow of closed links and serves as the hydraulic engine’s own negligible-flow threshold — in the tank level update of §6.3 and the tank and control time-step checks of §6.2 — but plays no role in the quality engine.)

Water age requires no source. The “concentration” is initialised from each node’s initial quality $C_0$ (zero by default, as noted in §8.1) and incremented by $\delta t$ (in hours) at every quality time step, representing the elapsed time since the water entered the system from a source (reservoir or tank).

Source tracing assigns a “concentration” of 100 (representing 100%) to all water leaving the designated trace node. All water entering the network from other reservoirs carries a concentration of zero permanently; a tank’s reported node quality starts at zero, though its stored volume segments are seeded from its initial quality $C_0$ (per §8.2), and tanks store and re-release whatever tracer reaches them through their mixing models. The trace value at any pipe segment or junction then represents the fraction of that water — expressed as a percentage — that originated from the traced source.

8.4 Chemical Reactions

Chemical decay or growth of the constituent occurs simultaneously with transport. Reactions are applied independently in each pipe segment and in each tank (as a whole, subject to its mixing model).

Bulk Reactions

Bulk reactions occur within the water volume and are governed by:

$$r_{\text{bulk}} = k_b \cdot f(c)$$

where $k_b$ is the bulk reaction rate coefficient and $f(c)$ is the concentration potential, which depends on the reaction order:

• Zero order ($n = 0$): $f(c) = 1$ — the reaction proceeds at a constant rate independent of concentration.
• First order ($n = 1$): $f(c) = c$ — the reaction rate is proportional to concentration; this covers simple first-order decay (e.g., chlorine demand).
• Second order ($n = 2$): $f(c) = c^2$.
• $n$-th order with limiting concentration $C_L$: $f(c) = c^{n-1} \cdot c_{\text{potential}}$, where the potential accounts for the approach to a limiting residual (decay toward a non-zero floor, or growth toward a ceiling).
• Michaelis-Menten kinetics (negative order): $f(c) = c / (C_L + \mathrm{sign}(k_b),c)$ — the denominator is $C_L + c$ for growth ($k_b > 0$) and $C_L - c$ for decay ($k_b < 0$), so the sign is governed by the sign of the rate coefficient. This models saturation kinetics — the reaction rate is approximately first-order at low concentrations and approximately zero-order at high concentrations relative to the half-saturation constant $C_L$.

Wall Reactions

Wall reactions represent the interaction of the constituent with the pipe wall material (e.g., disinfectant demand from biofilm or iron corrosion products). The mass transfer process has two stages in series: molecular diffusion from the bulk water to the wall, and the chemical reaction at the wall surface. Both stages must be overcome for mass to be transferred.

The analysis proceeds as follows:

1. Compute the Reynolds number $Re = VD/\nu$ and the Schmidt number $Sc = \nu / \mathcal{D}$, where $\nu$ is the kinematic viscosity of water and $\mathcal{D}$ is the molecular diffusivity of the constituent.

2. Compute the Sherwood number $Sh$, which characterises the ratio of convective to diffusive mass transfer:

   • Stagnant ($Re < 1$): $Sh = 2$ (pure diffusion limit)
   • Laminar ($1 \leq Re < 2300$): Graetz–Lévêque solution for developing concentration profiles in a tube:

   $$Sh = 3.65 + \frac{0.0668 ,(D/L),Re,Sc}{1 + 0.04,[(D/L),Re,Sc]^{2/3}}$$

   • Turbulent ($Re \geq 2300$): Notter–Sleicher correlation:

   $$Sh = 0.0149 , Re^{0.88} , Sc^{1/3}$$

   (The fractional exponents $2/3$ and $1/3$ are evaluated as the truncated constants 0.667 and 0.333.)

3. The mass transfer coefficient is:

   $$k_f = \frac{Sh \cdot \mathcal{D}}{D}$$

4. For first-order wall reactions ($n_w = 1$), the wall reaction rate $k_w$ (with units of velocity $[\text{m/s}]$) and the mass-transfer coefficient $k_f$ combine in series to give an effective first-order wall decay coefficient:

   $$k_{\text{eff}} = \frac{4}{D} \cdot \frac{k_w , k_f}{k_f + |k_w|}$$

   Here $4/D$ converts from a surface-area basis to a volume basis for a circular pipe. The series combination ensures that if either the diffusion step ($k_f$) or the wall reaction step ($k_w$) is slow, it dominates the overall rate.

5. For zero-order wall reactions ($n_w = 0$), the wall demand rate $k_w$ (converted to internal units as $k_w \cdot f_u^2$ where $f_u$ is the elevation unit conversion factor) and the concentration-dependent diffusive supply rate $c \cdot k_f$ are compared independently. The effective volumetric wall rate is $\mathrm{sgn}(k_w) \cdot \min(|k_w \cdot f_u^2|,; c \cdot k_f) \cdot 4/D$. When the diffusion boundary layer cannot supply mass as fast as the wall consumes it ($c \cdot k_f < |k_w|$), the reaction becomes mass-transfer-limited and concentration-dependent despite the nominally zero-order kinetics.

Setting the molecular diffusivity to zero disables the mass-transfer stage entirely: the Sherwood analysis is skipped, first-order wall kinetics act at their intrinsic rate ($k_{\text{eff}} = 4k_w/D$ with no $k_f$ limitation), and zero-order wall demand is never transfer-limited.

Combined Reaction in a Segment

The net concentration change in a pipe segment over a quality time step $\delta t$ is:

$$\Delta c = \left( r_{\text{bulk}} + r_{\text{wall}} \right) \delta t$$

This forward-Euler update is applied uniformly for all reaction orders. After each step the updated concentration is floored at zero ($c \leftarrow \max(0,, c + \Delta c)$) so that a decay reaction can never drive a constituent negative. The quality time step is kept short relative to the reaction time scale to keep truncation error acceptably small.

Roughness–Reaction Correlation

As an alternative to specifying a wall reaction coefficient $k_w$ for each pipe individually, EPANET supports a global roughness–reaction correlation factor $R_f$. When $R_f \neq 0$, the wall coefficient for each pipe is derived automatically from its roughness parameter. The correlation formula depends on the head-loss formula in use:

• Hazen–Williams ($C$ is the HW roughness coefficient, smoother pipes have higher $C$): $$k_w = \frac{R_f}{C}$$

• Darcy–Weisbach ($\varepsilon$ is the absolute roughness, $D$ the diameter): $$k_w = \frac{R_f}{|\ln(\varepsilon/D)|}$$

• Chezy–Manning ($n_M$ is Manning’s roughness, rougher pipes have higher $n_M$): $$k_w = R_f \cdot n_M$$

In all three cases, $R_f$ has units of $[k_w \cdot \text{roughness parameter}]$, chosen so that the resulting $k_w$ is in the wall-rate units expected by the wall reaction formulas — velocity (m/s or ft/s) for first-order wall kinetics, mass per unit area per unit time for zero-order. The physical motivation is that rougher pipe surfaces tend to harbour more biofilm or corrosion products and hence exhibit higher wall demand. Any pipe whose $k_w$ is set explicitly in the input takes precedence over the correlation.

9. Tank Mixing Models

Tanks use one of four models to govern how incoming water mixes with water already in storage. The choice of mixing model can significantly affect the predicted concentration of dissolved constituents leaving the tank.

Complete Mix (CSTR)

The tank is modelled as a Continuously Stirred Tank Reactor (CSTR). All water entering the tank is assumed to mix instantly and uniformly with the existing contents. The concentration at any instant is therefore uniform throughout the tank volume. The mixing update at each quality sub-step is:

$$c_{\text{new}} = \frac{c \cdot V + c_{\text{in}} \cdot V_{\text{in}}}{V + V_{\text{in}}}$$

where $V$ is the current stored volume, $c$ is the (uniform) tank concentration, $c_{\text{in}}$ is the volume-weighted inflow concentration, and $V_{\text{in}}$ is the inflow volume during the sub-step. Bulk reactions are applied separately (not during the mixing step) in the same phase as pipe reactions. If the tank is at capacity, the stored volume is clamped at $V_{\max}$ and the excess mass is booked as overflow outflow in the mass balance.

Two-Compartment Mix

The tank volume is divided into two compartments represented as two segments: an inlet mixing zone with a maximum capacity of $V_{\text{mz}} = f \cdot V_{\max}$ (where $f$ is a user-specified fraction, typically 0.1–0.3, and $V_{\max}$ is the maximum tank volume) and a stagnant zone comprising the remaining capacity $V_{\text{sz}} = V_{\max} - V_{\text{mz}}$. All inflow enters the mixing zone; all outflow also exits from the mixing zone. Transfers between zones are directional and discrete, not bidirectional or continuous:

• Filling ($v_{\text{net}} > 0$): inflow mass mixes into the mixing zone (weighted average). If the mixing zone volume would exceed $V_{\text{mz}}$, the excess volume $v_t = \max(0,; V_{\text{mix}} + v_{\text{net}} - V_{\text{mz}})$ is transferred from the mixing zone to the stagnant zone, carrying the mixing zone’s post-mix concentration. The stagnant zone concentration is updated as a volume-weighted average of its current contents and the transferred mass. If the stagnant zone’s volume would exceed $V_{\text{sz}}$, the surplus exits as overflow (counted as outflow in the mass balance) and the stagnant zone is clamped to $V_{\text{sz}}$.

• Emptying ($v_{\text{net}} < 0$): water is drawn back from the stagnant zone into the mixing zone to compensate for the net deficit: $v_t = \min(V_{\text{stag}},; |v_{\text{net}}|)$. The mixing zone concentration is updated as a volume-weighted average of its current contents, the inflow mass, and the transferred stagnant zone water.

• No net flow ($v_{\text{net}} = 0$): no volume transfer occurs between the zones, and — because there is no net volume change over the sub-step — the mixing-zone concentration is left unchanged; any inflow mass is not blended in during this case.

The outflow concentration is always the mixing zone concentration. This model captures the behaviour of elongated tanks where short-circuiting occurs — inflow water can exit before it fully mixes with the bulk stored water.

FIFO Plug Flow

The tank is treated as a perfectly ordered pipe with no axial mixing. Water enters from one end and exits from the other in strict first-in, first-out order. The segment representation used for pipes (§8.2) is applied directly to the tank. New inflow creates a new segment at the inlet end; outflow consumes segments from the outlet end. Reactions occur within each segment. This model is appropriate for narrow, tall standpipes or tanks with well-separated inlet and outlet ports. When the tank is full, the full inflow volume is withdrawn from the outlet end and the net inflow’s mass is booked as overflow.

LIFO (Stacked Layers)

Water enters and exits from the same end of the tank, as in a stratified system. Unlike FIFO, only the net flow moves segments — simultaneous inflow and outflow are netted against each other: net inflow forms a new segment at the top (or inlet side), while net outflow removes segments from that same end in last-in, first-out order. This model approximates thermal stratification in tanks where buoyancy prevents vertical mixing, so that recently added water leaves first. When the tank is full, the net inflow volume is removed from the opposite (first) end and that mass is booked as overflow.

10. Mass Balance

The simulator maintains a running mass balance for the quality constituent throughout the simulation. At each quality time step, the following quantities are accumulated:

• Initial mass stored: the total constituent mass in all pipes and tanks at the start of the simulation.
• Mass added from sources: constituent injected at network sources — including, implicitly, every reservoir’s outflow at its own fixed quality, source or not; conversely, all mass flowing into a reservoir counts as network outflow.
• Mass removed as outflow: constituent carried out of the network — consumer withdrawals, water leaving through reservoirs, and tank overflow are all accumulated in this term (not consumer demand alone).
• Mass reacted: constituent lost (or gained) through bulk and wall reactions; computed as the integral of reaction rates over all pipe segments and tanks.
• Final mass stored: the total mass remaining in the network at the end of the simulation.

The overall mass balance ratio is computed as follows. Let $m_\text{reacted}$ be the signed total mass change due to reactions (negative for growth, positive for decay). If $m_\text{reacted} > 0$ (net decay), it is added to the output side of the ledger. If $m_\text{reacted} < 0$ (net growth), its absolute value is added to the input side:

$$\text{ratio} = \frac{\text{mass outflow} + \max(m_\text{reacted}, 0) + \text{final mass stored}}{\text{initial mass stored} + \text{mass added by sources} + \max(-m_\text{reacted}, 0)}$$

A value close to 1.0 confirms that constituent mass is being conserved to within numerical precision. A significant deviation from 1.0 indicates either a numerical error or an inconsistency in the reaction parameterisation. This diagnostic is reported at the end of the simulation.

11. Energy Tracking

For each pump in the network, the hydraulic power consumed during each time step is:

$$P_{\text{hydraulic}} = \rho g Q , \Delta H$$

where $Q$ is the flow through the pump and $\Delta H$ is the head added. The actual power drawn from the electrical supply is:

$$P_{\text{electrical}} = \frac{\rho g Q , \Delta H}{\eta}$$

where $\eta$ is the pump efficiency. If an efficiency curve (efficiency versus flow) is provided, $\eta$ is read from the curve at the current operating point; otherwise a default efficiency is assumed. When a pump operates at a speed setting $\omega \neq 1.0$ and an efficiency curve is supplied, the efficiency is further adjusted using the Sarbu–Borza speed-correction formula:

$$\eta_{\omega} = 100 - \frac{100 - \eta_1}{\omega^{0.1}}$$

where $\eta_1$ is the efficiency read from the curve at the speed-adjusted operating point and $\eta_{\omega}$ is the corrected efficiency. For $\omega < 1$ the correction enlarges the efficiency gap $(100 - \eta_1)$, so the corrected efficiency is lower than the curve value (and slightly higher for $\omega > 1$) — the formula models the penalty of off-nominal operation. The result is clamped to the range 1–100%; the global default efficiency, used when no curve applies, is 75%.

The following energy statistics are accumulated over the simulation period for each pump:

• Kilowatt-hours consumed: the time integral of $P_{\text{electrical}}$ over the simulation.
• Time-weighted average efficiency: the average of $\eta$ weighted by the fraction of time spent at each operating point.
• Maximum demand: the peak value of $P_{\text{electrical}}$ observed at any time step.
• Cost: the product of energy consumed and a unit energy price, which may itself vary over time via a cost pattern.

All of these accumulate only while a pump is open — a closed pump contributes zero power and is skipped — so the averages are weighted over online time, not total duration. The raw accumulators are normalised at report time: time online becomes percent utilisation of the duration, the efficiency sum becomes its average over online hours, kilowatt-hours divided by online hours become average kW, and total cost is expressed per day (a single-period, zero-duration analysis accounts energy as one hour of operation). These energy diagnostics are essential for assessing the operating cost of different pumping schedules or for optimising pump dispatch.

Energy cost model: the unit energy cost $c$ (cost per kWh) used to accumulate TotalCost is determined at each time step as follows. A global base cost $c_0$ is multiplied by the current value of a global energy price pattern (if one is assigned), giving a time-varying rate $c_0 f(t)$. Each pump may also carry its own cost override $c_p$ and/or its own cost pattern, which are resolved independently of the global values: if a pump’s Ecost is positive it replaces $c_0$; otherwise the global $c_0$ is used. If a pump’s Epat is assigned it replaces the global pattern multiplier; otherwise the global pattern multiplier is applied even when the pump has its own cost override. The energy cost accumulated for pump $j$ over time step $\Delta t$ (hours) is therefore:

$$\text{Cost}_j \mathrel{+}= c_j(t) \cdot P_j \cdot \Delta t$$

where $c_j(t)$ is the applicable unit rate at the current time. Cost patterns advance on the same pattern start/step clock as demand patterns.

In addition to energy cost, a global peak demand charge parameter $D_c$ (cost per peak kW) is supported. A running maximum of the simultaneous power draw across all pumps $P_{\text{max}} = \max_t \sum_j P_j(t)$ is tracked throughout the simulation. At report time the total peak demand cost $D_c \cdot P_{\text{max}}$ is added to the energy cost summary. The KwHrsPerFlow statistic is accumulated at each time step as the time integral $\sum_i (P_i / Q_i) \cdot \Delta t_i$ of the instantaneous energy intensity $P/Q$, with the flow floored at $10^{-6}$ cfs; at report time it is divided by the pump’s online hours and unit-converted — kWh per million gallons (US) or kWh per m³ (SI) — as the measure of pumping energy per unit throughput.

12. Flow Balance

In addition to the local hydraulic solution at each time step, the simulation accumulates a global volumetric flow balance over the entire simulation period. Each quantity is computed by time-integrating the corresponding flow rate and dividing by total elapsed time at the end of the run — so every component, storage change included, is reported as a period-average flow rate:

• Total inflow: water entering the network from reservoirs (fixed-grade sources).
• Consumer demand delivered: water withdrawn at junctions as consumer demand.
• Emitter outflow: water discharged through emitters.
• Leakage outflow: water lost through pipe leakage modelled by the FAVAD equations.
• Demand deficit: the volume of consumer demand that was not delivered owing to insufficient pressure (relevant only in PDA mode; zero in DDA mode).
• Storage change: the net change in volume stored in all tanks over the simulation period (positive if tanks filled overall, negative if they drained).

The flow balance ratio is defined as:

$$\text{balance ratio} = \frac{q_\text{out}}{q_\text{in}}$$

where the ledger is built from the time-averaged flows as follows. When the net tank storage flow $q_\text{stor}$ is positive (tanks filling on average), it is added to the output side: $q_\text{out} = \text{total outflow} + q_\text{stor}$. When $q_\text{stor}$ is negative (tanks draining on average), its absolute value is added to the input side: $q_\text{in} = \text{total inflow} + |q_\text{stor}|$. Here total outflow comprises consumer demand, emitter flows, leakage, and water leaving the network into reservoirs (a reservoir with non-negative net demand accrues to the outflow side); total inflow is the supply from reservoirs plus any junction node with net negative demand. The demand deficit (undelivered demand in PDA mode) is tracked as a separate component and reported alongside the ratio but is not incorporated into the ratio calculation itself; it accumulates only at junctions with positive full demand and only when positive — over-delivery elsewhere never offsets it. A ratio close to 1.0 indicates that the simulation is globally mass-conserving; an exactly balanced (including empty) system reports exactly 1.0, and zero inflow with nonzero outflow reports 0.0. The instantaneous leakage fraction — leakage as a percentage of the current period’s supply, meaning reservoir and junction inflow plus net tank drainage, and zero when that supply is zero — is also tracked and can be reported at each reporting period to identify periods of high loss.

13. Units and Physical Constants

All hydraulic computations are performed internally in a fixed set of US customary units regardless of what units the user specifies in the input file: lengths in feet (ft), diameters in feet, flows in ft³/s (cfs), heads in feet, and power in horsepower (hp). Unit conversion factors are applied once during input parsing to translate user-supplied values into internal units, and again at output time to translate results back into user-facing units.

Unit system selection: the unit system (US or SI) is inferred automatically from the chosen flow unit:

Pressure units may be overridden independently to psi, kPa, metres, bar, or feet, regardless of the primary unit system. Reported pressures in psi, kPa, and bar scale with the specific gravity (the conversion factor embeds it); metres and feet do not.

Physical constants: the default values used in the absence of user overrides are:

For kinematic viscosity, the user may supply either a multiplier (value $> 10^{-3}$, interpreted as a scale factor on the default) or an actual value in ft²/s. For molecular diffusivity, the multiplier threshold is $10^{-4}$ rather than $10^{-3}$: values greater than $10^{-4}$ are treated as scale factors on the default diffusivity, while smaller values are taken as the actual diffusivity. When the SI unit system is active, supplied actual values for both quantities are converted from m²/s to ft²/s before storage. (The API’s viscosity/diffusivity options always interpret the supplied value as a multiplier of the default — the threshold rule applies to the input file only.)

Hydraulic solver defaults: the default values for convergence parameters that apply in the absence of user specification are summarised below.

The flow accuracy (Hacc) is bounded differently depending on how it is supplied, and the two paths disagree in both range and severity. Read from [OPTIONS] ACCURACY it is silently clamped to $[10^{-5}, 0.1]$, so an out-of-range value in a file loads and runs at the nearest bound. Set through the API it is rejected outside $[10^{-8}, 0.1]$, raising the invalid-option-value error. The shared ceiling of $0.1$ is common to both; the floor is three orders of magnitude tighter on the API path, so the engine’s own bounds place the limit of what its solver is expected to reach at $10^{-8}$ while the file path stops well short of it.

Time parameter defaults: the time settings default and interlock as follows, all resolved once at the end of input processing.

The hydraulic step is then lowered to the pattern step and to the reporting step if it exceeds either, and the quality step lowered to the hydraulic step, as §6 describes. One further adjustment is a silent substitution rather than a clamp, and its direction is the opposite of what one would expect: a report start time later than the simulation duration is reset to zero, not to the duration — so a run configured to begin reporting after it ends reports everything instead of nothing.

14. Input and Output

Input

The network is described in a structured plain-text input file organised into labelled sections. Each section corresponds to a class of network object or a simulation parameter group. The parser makes two passes through the file. A first pass does two things: it counts all objects of each type so that memory can be allocated in one contiguous block, and it also captures the UNITS and HEADLOSS options from the [OPTIONS] section — so that second-pass interpretation is independent of where [OPTIONS] appears in the file; the object counting itself depends on neither option. A second pass reads and interprets all remaining data. The sections handled include: junctions, reservoirs, tanks, pipes, pumps, valves, demand categories, time patterns, head-loss and efficiency curves, simple controls, rule-based controls, water quality sources, emitter coefficients, leakage parameters, options (head-loss formula selection, flow units, demand model, tolerances), energy pricing, reaction coefficients, tank mixing model assignments, reporting options, initial status overrides, and simulation time parameters — plus initial water quality ([QUALITY]), the map/meta sections ([COORDINATES], [VERTICES], [LABELS], [BACKDROP], [TAGS]) that are round-tripped but never simulated, and a legacy [ROUGHNESS] section that is accepted and silently ignored.

Keyword matching is by prefix, and this is a compatibility fact rather than a parsing detail. One routine serves every keyword lookup in the file — section headers, option names, and option values alike — walking a keyword table in order and returning the first entry that is a prefix of the supplied token, compared without regard to case after leading blanks are skipped. Three consequences follow. Trailing characters are ignored, so [JUNCTIONS], [JUNC] and [JUNCTIONSXYZ] all select the same section; truncations are not, since the comparison runs to the end of the keyword. The keyword table exploits this deliberately, storing stems rather than words — MESS, VERI, Junc, Reser, Tank — which only resolve because the file’s full spellings extend them. And table order becomes load-bearing wherever one keyword prefixes another: the ordering is part of the format’s meaning, not an incidental property of the table.

The engine hit that last hazard once and patched around it by hand. In the report-field table Head precedes Headloss and is a prefix of it, so a [REPORT] HEADLOSS line would select head reporting instead. The parser therefore tests the token against the literal Headloss with a whole-string comparison before consulting the prefix matcher at all — a special case existing solely to defeat the general rule. It is the only such guard in the parser.

Initial link status overrides ([STATUS] section): the [STATUS] section allows the user to set the initial operational status of any link before simulation begins. For pipes, valid entries are OPEN or CLOSED. For pumps, an entry may be OPEN, CLOSED, or a numeric relative-speed setting (a value of 0 closes the pump). For valves, an entry may be OPEN, CLOSED, or a numeric setting that overrides the value from the [VALVES] section. Because [STATUS] is parsed after the link-definition sections, its entries take precedence over inline status values. During simulation, controls and rules may subsequently change these statuses.

Before simulation begins, the project undergoes a validation pass that checks: each tank satisfies $H_{\min} \leq H_{\text{init}} \leq H_{\max}$, with levels inside any volume curve’s elevation range; all patterns are non-empty; all curves are non-empty with strictly increasing $x$-values; pump curves follow the §2.2 admissible forms (one point, three points starting at zero flow, or a strictly-decreasing custom curve, with $0 < c \leq 20$); and the network holds at least two nodes and at least one tank or reservoir. Validation runs when the hydraulic solver opens, not at file load. (Valve topology, by contrast, is enforced against the §1 placement rules as each valve is defined — at input parse time and again on any API edit that changes a valve’s type or end nodes — not by this validation pass.) The unconnected-node check covers junctions only — an unlinked tank is harmless and silently permitted — and reports up to ten offenders before failing. If any check fails, the simulation does not start.

One step between parsing and validation mutates the network rather than checking it. Node indices are assigned junctions first, then tanks and reservoirs, with the junction block sized from the count taken in the first pass. When fewer junctions are created than [JUNCTIONS] lines were counted — which happens when a line’s identifier fails to register, most obviously a duplicate — the junction block is left with a gap, and every tank is shifted down to close it. That shift rewrites all references: link end nodes, simple-control node references, rule references, the quality trace node, and the identifier hash table. Because node index order is the order in which the output file’s identifier table and every per-period result vector are written (§14), this shift propagates into the output file’s layout. The path is reachable only through input that already raised a non-fatal error, so for a well-formed file the resulting order is simply junctions in file order followed by tanks and reservoirs in file order.

Alternatively, networks may be constructed entirely through the project API without reference to an input file, by calling the object-creation and property-setting operations in programmatic sequence; the in-memory network can also be regenerated as an input file, with the [LABELS] and [BACKDROP] sections copied verbatim from the original while [COORDINATES], [VERTICES], and [TAGS] are rewritten from the in-memory data. Deleting a pattern or curve renumbers all higher references — references to the deleted pattern are cleared, which for demand patterns means falling back to the default pattern at run time (source, energy, and tank-head references simply become pattern-less) — and a PCV losing its curve has its loss coefficient zeroed.

Output

Hydraulic binary file: at each hydraulic time step the solver writes nodal demands, nodal heads, link flows, link status flags, and link settings to a temporary binary file. This file is then replayed during the water quality simulation, supplying the flow field from which velocities are reconstructed for advection, without requiring the hydraulics to be recomputed. (Velocities themselves are not stored; they are derived from the saved flows and pipe geometry.) The file opens with a 2+6 INT4 header — magic number 516114521, engine version, then node/link/tank/pump/valve counts and duration — validated field-for-field when a saved file is reused, a mismatch being fatal; each snapshot holds an INT4 time, REAL4 demands, heads, flows, status codes and settings, and the next-step interval. The API can save and reuse this file across runs, and can populate the results file from saved hydraulics without a water-quality run.

Results binary file: at each reporting time step (which may be less frequent than the hydraulic time step), all computed quantities for every node and link are saved to a separate binary output file. Node quantities include hydraulic head, pressure, demand, and constituent concentration. Link quantities include flow rate, velocity, unit head loss, friction factor, and quality. This file may subsequently be post-processed by external programs (a bundled reader library exposes a standalone API over the same format). When the report-statistic option is anything other than time-series, results accumulate in a temporary file and the dynamic section holds exactly one pseudo-period of averaged/minimum/maximum/range values — flows as magnitudes, statuses collapsed to open/closed — with the period count written as 1. The text report is generated by replaying this file, so results must be saved before reporting.

Binary Output File Format

The binary output file (.out) is written in native byte order (little-endian on x86) using float (4-byte IEEE 754 single-precision, hereafter REAL4) for all floating-point values and int (4-byte signed, hereafter INT4) for integers. String fields are fixed-width arrays: IDs are 32 bytes (MAXID+1 = 32, null-terminated), title lines are 80 bytes (TITLELEN+1 = 80), and filenames are 260 bytes (MAXFNAME+1 = 260). No padding or alignment bytes exist between sections.

The file has five sections written sequentially:

Prolog

15 × INT4 header (60 bytes), then strings and arrays:

Following the ID strings:

Energy

Written immediately after the prolog, once per simulation:

Per pump ($N_p$ records of 28 bytes each):

Followed by one trailing REAL4: demand charge (peak demand × demand cost rate).

Dynamic Results

Written once per reporting period. Each period contains the following arrays, all column-major (one variable across all objects, then the next variable):

Node variables — 4 arrays of $N_n$ × REAL4, in output units:

Link variables — 8 arrays of $N_l$ × REAL4:

Bytes per period: $(4 N_n + 8 N_l) \times 4$.

Network Reactions

4 × REAL4 (16 bytes):

Epilog

3 × INT4 (12 bytes):

The total file size is: $$884 + 36 N_n + 52 N_l + 8 N_t + (28 N_p + 4) + N_{\text{periods}} \cdot 4(4 N_n + 8 N_l) + 16 + 12$$

Text status report: an optional text-format status report records, at user-specified verbosity, the convergence history of every hydraulic time step (number of iterations, peak head error, flow accuracy achieved), any link status changes during the simulation, the energy consumption and cost summary for all pumps, the final mass balance ratio for water quality, the final flow balance statistics, and — if requested — tabular node and link results at every reporting period.

API

The system exposes a complete project-handle–based API. The workflow is as follows:

1. Create a project object, which encapsulates all state for a single simulation instance.
2. Open a network description (from file or by programmatic construction).
3. Optionally, open a pre-existing hydraulics results file; if none exists, run the full hydraulic simulation first.
4. Run the hydraulic simulation either in full (computing all time steps internally) or step-by-step; in step-by-step mode the caller advances the clock one hydraulic time step at a time and may modify scalar properties and settings between steps — structural changes (adding or deleting nodes and links; deleting curves or patterns; changing a link’s type or end nodes) are refused while either solver is open — though adding a pattern or curve is permitted — and the head-loss formula cannot change while the hydraulic solver is open, nor the quality model while the quality solver is open.
5. Run the water quality simulation, either in full or step-by-step, in a similar fashion.
6. Retrieve any computed result (pressure, flow, concentration, energy, etc.) at any time step.
7. Set any scalar network property (demand, pipe roughness, pump speed, valve setting, control threshold, reaction coefficient, etc.) and re-run as desired; structural edits require the solvers to be closed first.
8. Delete the project to release all resources.

Beyond this workflow, the API surface also covers: full network editing (creating, deleting, and renaming every object class, including demands, patterns, curves, controls, and rules, plus vertex/coordinate and comment/tag metadata); per-control and per-rule enable/disable; solver introspection (iteration statistics, result indices, time to the next event); report control (generation, custom lines, copying and resetting); one-shot run drivers; and regeneration of the input file from memory. A legacy function family wraps a single global default project — the multi-instance concurrency described below applies only to the handle-based API.

Multiple project instances may coexist in the same process, enabling Monte Carlo analysis, parallel scenario evaluation, or re-entrant simulation from multiple threads (provided each thread operates on a distinct project handle and any shared file system resources are managed appropriately).

15. Cross-Cutting Engine Contracts

The preceding sections follow EPANET’s physical subsystems; this one collects the conceptual conventions that span them, each referenced back to the sections carrying its details.

The internal-units contract. All hydraulic computation runs in a fixed US customary system — lengths and heads in feet, flows in ft³/s, power in horsepower — with conversion applied at the boundaries: once at input parsing and again at output time (§13). One interior exception exists and it is easy to miss because it converts back: user-defined curves are stored untransformed, so extracting a linearised segment from a custom pump curve (§2.2) or a GPV head-loss curve (§2.3) converts the current flow into the curve’s own units, fits the bracketing segment there, and converts the resulting slope and intercept back to internal units. The arithmetic is the same either way, but a re-implementation that normalises curve points at parse time and then reuses this segment-fitting logic unchanged will double-convert. Interior constants are therefore expressed in these units: the TCV loss-coefficient factor $0.02517,s/D^4$ (§2.3), the FAVAD orifice coefficient $C_o$ (§5), the constant-power pump’s nominal design flow of 1 ft³/s (§2.2), and the quality thresholds Q_STAGNANT and QZERO (§8.3). Specific gravity couples into the pressure boundary alone: reported pressures in psi, kPa, and bar embed it in their conversion factors, while metres and feet do not (§13), and in US units the emitter coefficient is defined per psi$^{n_e}$ adjusted for specific gravity (§4).

The status-machine convention. Every link carries a discrete status re-evaluated during the solve: check-valve pipes toggle OPEN/CLOSED (§2.3), pumps pass through OPEN, XHEAD, and TEMPCLOSED (§2.2), PRVs and PSVs inhabit ACTIVE/OPEN/CLOSED plus XPRESSURE and FCVs ACTIVE/XFCV (§2.3), and tank-adjacent links are forced TEMPCLOSED at tank limits (§6.3). All transition tests share the two tolerances Htol and Qtol, which serve status logic only and are entirely separate from the convergence tolerance Hacc governing solver termination (§2.3, §13). The checks follow a re-open-then-re-test pattern — TEMPCLOSED and XHEAD links are first re-opened at the start of each status check and the new operating point re-evaluated, so no status is locked in (§2.2, §6.3) — and convergence itself requires a pass that produces no status change (§2.4).

The smooth-barrier idiom. One differentiable one-sided penalty, $\Delta h = (a - \sqrt{a^2 + 10^{-6}})/2$ with $a = 10^9 Q$ and its matching gradient increment, is the engine’s device for enforcing inequality bounds without breaking the Newton–Raphson iteration: it holds pressure-driven demands within $[0, D_{\text{full}}]$ (§3.2, the only site using both the lower and the upper barrier), drives emitter flow non-negative when backflow is forbidden (§4), and holds leakage flow non-negative with no equivalent opt-out (§5) — the same closed form at all three, approaching a hard constraint as the bound is violated while remaining smooth throughout. The upper barrier is the lower one’s mirror, $\Delta h = (a + \sqrt{a^2 + 10^{-6}})/2$ with gradient $\tfrac{1}{2}10^9(1 + a/\sqrt{a^2+10^{-6}})$. The leakage site reaches it through a private duplicate of the routine — same body, different name — rather than a shared call, so a reader tracing callers of the shared helper alone will find two sites and miss the third.

The $(P, Y)$ linearisation contract. Every head–flow element reduces at each iteration to the same pair — a conductance $P_k$ (inverse head-loss gradient) and an offset $Y_k$ — consumed identically by the GGA assembly: $P$ into the Laplacian diagonal and off-diagonals, $Y$ into the right-hand side (§2.4). Pipes under all three friction formulas with minor losses (§2.1), pumps as negative head loss (§2.2), the resistance-type valves TCV, PBV, GPV, and PCV (§2.3), emitters (§4), pressure-dependent demands (§3.2), and the two FAVAD leakage components as equivalent emitters (§5) all enter through this one algebra. Its edge cases are handled inside the same coefficients: the RQtol gradient floor keeps $P$ bounded near zero flow (§2.1), and the large constant $C_\infty \approx 10^8$ turns the pair into a near-exact constraint — the PBV’s fixed head drop, the active PRV/PSV head pinning, and the constant-power pump’s near-zero-flow closure (§2.2, §2.3). Only the ACTIVE control valves step outside it, zeroing $P_k$ and augmenting the matrix directly (§2.3).

Prefix keyword matching. Every keyword in the input file — section header, option name, option value — is resolved by one routine that returns the first table entry forming a prefix of the supplied token, case-insensitively (§14). The keyword table is written to suit it, holding stems rather than words, and the rule’s one harmful collision (Head shadowing Headloss) is defeated by a hand-written whole-string test rather than by reordering. Matching is therefore strictly looser than equality, and table order is part of the file format’s meaning.

The curve-clamping rule. User-defined curves are piece-wise linear: intermediate values interpolate linearly between the two bracketing data points, and values outside a curve’s x-range are clamped to its end-point values — curves never extrapolate (§1). The PCV opening curve is the one documented departure, interpolating from the origin below its first point and toward the $(100%, 100%)$ anchor above its last, with the resulting ratio clamped to $[10^{-6}, 1]$ (§2.3). Validation enforces that every curve is non-empty with strictly increasing $x$-values and that pump curves take one of the admissible forms of §2.2 (§14).

The pattern-clock convention. All patterns advance on one clock: the period index $p = \lfloor (t + t_{\text{start}}) / \Delta t_p \rfloor$ selects multiplier $F_j[p \bmod L_j]$, giving each pattern an independently repeating cycle (§1). The same convention drives junction demands, reservoir heads, and pump speed settings (§1), source-concentration multipliers (§8.3), and energy cost patterns (§11). The fallback for an unassigned demand category is not a special case in the arithmetic: pattern index 0 is a synthetic single-entry pattern holding the multiplier 1.0, installed before parsing begins, so a category with no pattern resolves to the default pattern and a project with no default resolves to index 0 — the unit multiplier arrives through the same lookup as every other (§1).

“One clock” is true of the result but not quite of the mechanism. Four of the five consumers read the hydraulic clock; source-concentration multipliers read the quality clock, which advances on its own finer step within each hydraulic period (§8.1). The two agree on the period index only because the hydraulic step is trimmed so that no pattern boundary falls strictly inside it — the step limit below is therefore not merely a convenience for accuracy but the thing that keeps the two clocks reporting the same pattern period. The adaptive time step conspires with this clock — the next pattern change is one of the minimised step limits, so multipliers change exactly at step boundaries, never mid-step (§6.2).

SWMM: A Conceptual and Mathematical Analysis

Hydra implements none of this yet. This document analyses SWMM itself. It is groundwork for the urban drainage engine, which is registered but planned — no implementation exists behind it, and nothing below describes behaviour you can run in Hydra today.

It is also a pinned snapshot: it describes SWMM 5.2.4 at tag OWA_v5.2.4 faithfully, defects included. Where it records a bug or an inconsistency, that is a finding about SWMM — not a Hydra issue to file.

Introduction

SWMM (Storm Water Management Model) is a computational engine for simulating the quantity and quality of runoff from urban catchments and its conveyance through drainage systems — storm sewers, sanitary and combined sewers, open channels, storage units, and flow regulators — over single events or continuous multi-year periods. Where a water distribution engine solves for a pressurised network in hydraulic equilibrium at each time step, SWMM is a rainfall-runoff-routing model: precipitation falls on subcatchments, becomes runoff after losses to infiltration, evaporation, and depression storage, and the resulting hydrographs and pollutographs are routed through a free-surface conveyance network by solving forms of the Saint-Venant equations. Flows may be driven by gravity or pumped, conduits may transition between open-channel and pressurised (surcharged) states, and the network may include backwater, flow reversal, ponding, and tidal boundary conditions.

This document provides a self-contained, mathematical and conceptual description of every major subsystem: how the physical system is represented as modelling objects, how runoff is generated on subcatchments, how infiltration and other hydrologic losses are computed, how flows are routed through the conveyance network under the steady, kinematic wave, and dynamic wave options, how special structures (pumps, orifices, weirs, outlets, culverts, force mains) behave, how control rules operate, how pollutants build up, wash off, and are transported and treated, and how continuity is tracked. The goal is to give the reader a complete algorithmic and mathematical understanding of the system; implementation-specific details such as memory layout are omitted, but input/output behaviour is described precisely. The analysis describes SWMM 5.2.4 — tag OWA_v5.2.4 of the community-maintained repository — with the EPA SWMM Reference Manuals (Volume I, Hydrology, EPA/600/R-15/162A; Volume II, Hydraulics, EPA/600/R-17/111; Volume III, Water Quality, EPA/600/R-16/093) as secondary references. Wherever the manuals and the source code disagree, the source is authoritative, and each such discrepancy is noted in place.

Table of Contents

• SWMM: A Conceptual and Mathematical Analysis
  • Introduction
  • Table of Contents
  • 1. System Representation
    • 1.1 Environmental Compartments
    • 1.2 Hydrology Objects
    • 1.3 Conveyance Nodes
    • 1.4 Conveyance Links
    • 1.5 Water Quality Objects
    • 1.6 Data Objects
    • 1.7 The State-Vector View
  • 2. Simulation Architecture
  • 3. Hydrology
    • 3.1 Meteorology
    • 3.2 Surface Runoff
    • 3.3 Infiltration
      • Horton
      • Modified Horton
      • Green–Ampt
      • Modified Green–Ampt
      • Curve Number
    • 3.4 Groundwater
      • Zone Structure and Fluxes
      • Constitutive Relations
      • Lateral Discharge to the Drainage System
      • The ODE Pair in Manual and Source
    • 3.5 Snowmelt
    • 3.6 Rainfall-Dependent Inflow/Infiltration (RDII)
  • 4. Flow Routing Theory
  • 5. Dynamic Wave Analysis
    • 5.1 Conduit Flow Update
    • 5.2 Node Head Update
    • 5.3 Picard Iteration
    • 5.4 Surcharge
    • 5.5 Flooding and Ponding
    • 5.6 Time-Step Control
    • 5.7 Initial Conditions
  • 6. Cross-Section Geometry
    • 6.1 Shape Families
      • Analytic Shapes
      • Circular Sections
      • Composite Shapes
      • Tabulated Shapes
      • Custom Shapes
      • The Section-Factor Maximum
    • 6.2 Transects
    • 6.3 Storage Geometry and Characteristic Depths
  • 7. Pumps and Flow Regulators
    • 7.1 Pumps
    • 7.2 Orifices
    • 7.3 Weirs
    • 7.4 Outlets
  • 8. Advanced Hydraulics
    • 8.1 Conduit Evaporation and Seepage
    • 8.2 Minor Losses and Force Mains
      • Minor Losses
      • Force Mains
    • 8.3 Culverts and Roadway Weirs
      • Culverts
      • Roadway Weirs
    • 8.4 Streets and Inlets
      • On-Grade Capture
      • On-Sag Capture
      • Capture Transfer and Statistics
  • 9. Water Quality
    • 9.1 Pollutants and Sources
    • 9.2 Buildup and Street Sweeping
    • 9.3 Washoff
    • 9.4 Transport and Treatment
      • Conveyance Transport
      • Treatment
  • 10. LID Controls
    • 10.1 The Generic Layered Unit
    • 10.2 Constitutive Fluxes
      • Per-Type Flux Balances
      • The Limiter Cascade
    • 10.3 Unit-Type Variants
    • 10.4 Integration
    • 10.5 Deployment and Routing of Outflows
    • 10.6 Water Quality in LIDs
  • 11. Control Rules
    • 11.1 Premises
    • 11.2 Actions and Conflict Resolution
    • 11.3 PID Controllers
    • 11.4 The Expression Language
  • 12. Continuity Accounting
    • 12.1 The Five Balances
    • 12.2 Reported Statistics
  • 13. Units and Physical Constants
  • 14. Input and Output
    • Input
      • File Grammar
      • Prefix Matching
      • Sections
      • Reporting Directives
      • Deprecated and Inert Options
      • Process Switches and Time-Step Interlocks
      • Semantic Validation
    • Interface Files
    • Output
      • Report File Structure
      • Binary Output File Layout
  • 15. The Engine as a Library
    • Run-Loop Lifecycle
    • Object and Property Access
    • Query and Mutation
    • Rainfall Injection and Checkpointing
  • 16. Cross-Cutting Engine Contracts

1. System Representation

1.1 Environmental Compartments

SWMM conceptualises an urban drainage system as water and material flows between four environmental compartments:

• The Atmosphere compartment generates precipitation and deposits pollutants onto the land surface. It is represented by rain gage objects.
• The Land Surface compartment receives precipitation as rain or snow and loses water through evaporation back to the atmosphere, infiltration into the sub-surface, and surface runoff (with its pollutant load) into the conveyance system. It is represented by subcatchment objects.
• The Sub-Surface compartment receives infiltration from the land surface and transfers a portion of it to the conveyance system as groundwater interflow. It is represented by aquifer objects.
• The Conveyance compartment is the network of channels, pipes, pumps, regulators, and storage units that carries water to outfalls or treatment. It is represented as a directed graph of nodes and links. Inflows to this compartment can come from surface runoff, groundwater interflow, rainfall-dependent infiltration/inflow, sanitary dry-weather flow, or user-defined time series.

Not every compartment need be present in a model: a pure conveyance model may be driven entirely by user-supplied inflow hydrographs, and a pure hydrology model may end at subcatchment outlets. This compartmental decomposition is the fundamental structural difference from a pressurised-network model such as EPANET, in which the network is the entire model: in SWMM the node-link conveyance graph is only one of four coupled subsystems, and areal objects (subcatchments, aquifers, snow packs) participate in the simulation without being part of the graph at all.

1.2 Hydrology Objects

Rain gages supply precipitation to one or more subcatchments. A gage’s data come from a user-supplied time series or an external rainfall file, expressed as intensity, volume, or cumulative volume over a fixed recording interval.

Subcatchments are parcels of land that receive precipitation from exactly one rain gage and generate runoff and pollutant loads. A subcatchment discharges either to a conveyance node or to another subcatchment, allowing overland flow to cascade across parcels. Each subcatchment is idealised as a rectangular plane of a given area, characteristic width, and uniform slope, partitioned into three sub-areas: an impervious fraction with depression storage, an impervious fraction without, and a pervious fraction (with depression storage). The five geometric parameters (area, imperviousness, width, slope, curb length) are rejected only for being negative; an imperviousness above 100% is not rejected but silently capped at 100%. Only the pervious fraction infiltrates. Runoff from each sub-area may optionally be re-routed onto another sub-area rather than directly to the outlet, modelling e.g. rooftops draining onto lawns. Overland flow is generated by treating each sub-area as a nonlinear reservoir (§3). A subcatchment may additionally host: a snow pack object governing snow accumulation and melt on its plowable, impervious, and pervious fractions; a groundwater connection to an aquifer; LID controls occupying a portion of its area (layered as surface, pavement, soil, storage, underdrain, and — on green roofs — drainage mat); and per-land-use pollutant buildup state.

Aquifers are two-zone (unsaturated/saturated) sub-surface reservoirs placed beneath subcatchments. They receive percolation from the surface, lose water to deep percolation and evapotranspiration, and exchange flow with a designated conveyance node through a user-parameterised groundwater flow equation — the mechanism by which baseflow and groundwater infiltration enter sewers.

Unit hydrographs (organised in groups of up to three per month, the “RTK” parameterisation) describe rainfall-dependent infiltration/inflow (RDII) — the delayed entry of stormwater into sanitary sewers through defects and illicit connections. Each unit hydrograph converts a unit of instantaneous rainfall into a triangular response defined by a fraction of rainfall volume ($R$), a time to peak ($T$), and a recession ratio ($K$).

LID controls are depth-explicit representations of low-impact-development practices — bio-retention cells, rain gardens, green roofs, infiltration trenches, permeable pavement, rain barrels/cisterns, rooftop disconnection, and vegetative swales — composed of layered storage elements (surface, soil, storage, drain) with their own governing equations (§10). They are defined once and deployed in multiple subcatchments at specified sizes.

1.3 Conveyance Nodes

Nodes are the points of the conveyance graph. Every node has an invert elevation and may receive external inflows, in addition to the runoff and groundwater delivered by hydrology objects. A direct external inflow (of flow or of any constituent) is the composite $cf,(sf \cdot TS(t) + \text{baseline} \cdot P(t))$ — an optional time series with scale factor $sf$, plus a constant baseline modulated by its own monthly/daily/hourly/weekend pattern, times a units factor $cf$ (mass-type pollutant inflows carry their own conversion factor); dry-weather sanitary inflow multiplies an average value by up to four patterns of distinct types, the weekend-hourly pattern replacing (not multiplying) the hourly one on weekends; and RDII arrives per §3.6. There are four node types:

Junctions are ordinary connection points — manholes, pipe fittings, or channel confluences — with negligible storage volume. A junction has a maximum (ground/rim) depth; when the hydraulic grade line reaches it, excess water is either lost from the system or, if ponding is enabled, stored atop the node in a user-specified ponded area and returned as capacity recovers.

Outfalls are terminal boundary nodes where water leaves the system to a receiving body. An outfall’s boundary stage may be: free (the smaller of critical and normal depth at the connecting conduit), normal (normal-depth), fixed (constant stage), tidal (a repeating 24-hour stage curve, indexed by elapsed routing time from the curve’s first hour — coinciding with clock time only for midnight starts), or a time series. For the staged variants, the stage governs only when it exceeds the critical-depth elevation; below that the node sits at critical depth, except that a conduit on a positive offset whose invert lies above the stage leaves the node at the stage height and the conduit in free fall. Outfalls may carry their own flap gate blocking reverse flow. An outfall connects to exactly one link (of any type — pumps and regulators may discharge directly to it; under steady and kinematic wave it may have no outlet links while its inflow count goes unchecked, whereas dynamic wave enforces a single connecting link that may even be an outlet link) and may optionally route its discharge back onto a subcatchment. Under steady and kinematic wave routing, any non-outfall, non-storage terminal node with no outlet links behaves identically — its inflow leaves the system without being counted as flooding; under dynamic wave such a node is an ordinary interior node whose overflow counts as flooding.

Storage units are nodes with significant free-surface storage volume — ponds, wet wells, detention basins, chambers. Their geometry is described either by a functional relation (surface area as a power function of depth), by a tabulated area-versus-depth storage curve, or by one of four analytical shapes new in 5.2 (elliptical cylinder, elliptical cone, elliptical paraboloid, rectangular pyramid). Storage units may lose water through evaporation and through seepage into the soil.

Dividers split inflow (including any node overflow) between two outflow conduits according to a prescribed rule: a cutoff divider diverts all inflow above a threshold; an overflow divider diverts whatever the non-diverted conduit declines to accept (the non-diverted link is deliberately routed first); a tabular divider uses a diverted-flow-versus-inflow curve; and a weir divider applies the weir equation at head $f,d_{max}$ — discharging $C_W (f,d_{max})^{1.5} = q_{max}f^{3/2}$ — with fraction $f = (Q_{in} - q_{min})/(q_{max} - q_{min})$, $q_{max} = C_W d_{max}^{1.5}$, switching to an orifice-like $q_{max}\sqrt{f}$ when surcharged ($f > 1$). Diverted flow is clamped to the inflow — except on the overflow path, which returns its split before the clamp is reached. The weir and tabular rules are evaluated in the user’s flow and length units (§16). Dividers are only meaningful under steady and kinematic wave routing; under dynamic wave analysis they behave as ordinary junctions, since the full momentum treatment determines the flow split naturally.

1.4 Conveyance Links

Links connect nodes and carry flow; a link’s orientation defines positive flow direction, and negative flows denote reversal. There are five link types:

Conduits are the pipes and channels of the network — the only link type with hydraulic length and the primary object of flow routing. Each conduit has a cross-sectional shape drawn from one of four descriptions:

• a library of more than twenty standard closed and open geometries (circular, rectangular, trapezoidal, egg, horseshoe, arch, elliptical, and others);
• an irregular transect (a surveyed station-elevation profile with left-bank/main-channel/right-bank roughness, in the manner of HEC-2/HEC-RAS river sections);
• a street cross-section (a curb-and-gutter roadway profile for dual-drainage street modelling); or
• a user-supplied custom shape curve of width versus depth.

Conduits carry:

• a Manning roughness coefficient;
• upstream and downstream invert offsets above their end nodes (expressible as heights or, via a global option, absolute elevations);
• an optional count of identical parallel barrels (routing solves one barrel at $Q/N$ and scales volumes and losses back by $N$);
• an optional maximum-flow limit;
• optional entrance/exit minor loss coefficients;
• an optional flap gate preventing reverse flow; and
• optional seepage and evaporation losses.

Conduit slope is drop over horizontal distance, $S_0 = \Delta z/\sqrt{L^2 - \Delta z^2}$, floored at a 0.001-ft elevation drop (or a user minimum slope) and — in the degenerate case where the drop equals or exceeds the length, which would make the horizontal distance imaginary — falling back to $\Delta z/L$ with a warning. The user minimum slope is applied last and, under steady or kinematic routing, returns immediately — before the step that gives an adverse conduit its negative sign — so a conduit shallow enough to be clamped is also silently made downhill. A negative offset at either end is silently zeroed, also with a warning. Under dynamic wave an adverse-slope conduit is silently reversed internally, with all reported flows carrying a direction multiplier so output keeps the user’s orientation. A conduit may also be designated a culvert (activating inlet-control capacity limits per FHWA HDS-5) or a force main (using Hazen–Williams or Darcy–Weisbach friction while pressurised, §8).

Pumps raise water between nodes according to a pump curve of five types:

• Type 1 (flow varies stepwise with wet-well volume);
• Type 2 (flow varies stepwise with inlet depth);
• Type 3 (flow varies continuously with delivered head);
• Type 4 (flow varies continuously with inlet depth);
• Type 5 (a variable-speed Type 3);

plus an ideal transfer pump whose outflow equals its inflow. Pumps have on/off depth setpoints and may have their speed modulated by control rules.

Orifices are openings in the side (side orifice) or bottom (bottom orifice) of a node’s wall or floor, closable to a variable degree by control rules, discharging according to the orifice equation with distinct free, submerged, and partially-open regimes. Orifice geometry is circular or rectangular; an optional flap gate prevents reverse flow.

Weirs are overflow structures of five types, each with its characteristic head-discharge exponent and coefficient, with corrections for end contractions, submergence, and surcharge:

• transverse;
• side-flow;
• V-notch (triangular);
• trapezoidal; and
• roadway (an FHWA HDS-5 embankment-overtopping weir for culvert/roadway systems).

Outlets are general-purpose head-discharge devices: their outflow is an arbitrary user-defined function of head or depth, given either as a power function or a rating curve. They model devices with bespoke ratings — vortex valves, flow-duration-control devices — that fit none of the standard structures.

(A sixth object family new in 5.2 — streets and inlets — pairs street cross-sections with FHWA HEC-22 inlet capacity calculations to model dual drainage: flow on the street surface, captured by inlets, entering the below-ground sewer. It is treated in §8.)

1.5 Water Quality Objects

Pollutants are user-defined constituents (any number) carried by runoff and routed through the conveyance system. Each has a concentration unit (mg/L, µg/L, or counts/L), optional rainfall/groundwater/RDII/dry-weather-flow background concentrations, an initial network concentration, a snow-only buildup flag, a first-order decay coefficient, and optionally a co-pollutant relationship (its concentration set as a fixed fraction of another pollutant’s).

Land uses partition a subcatchment’s area into categories (residential, industrial, …) that govern pollutant buildup during dry weather and washoff during runoff, each by a choice of functional forms (§9), plus street-sweeping parameters that periodically remove accumulated buildup.

1.6 Data Objects

Curves are tabulated x-y relations, linearly interpolated (except Type 1/2 pump curves, which are read stepwise), serving typed roles: storage (area vs. depth), diversion (diverted flow vs. inflow), tidal (stage vs. hour), pump (per the five pump types), rating (outlet discharge vs. head), control (setting vs. controller variable), shape (conduit width vs. depth), and weir coefficient curves.

Time series are timestamped value sequences used for rainfall, outfall stage, external inflows, and evaporation. Their reader is cursored on date and carries an extend flag that splits the consumers in two: outfall stage, temperature, and control-rule TIMESERIES actions hold the first or last value outside the series’ range, while external inflows and external buildup return zero — so an inflow series that ends before the simulation does silently falls to nothing rather than holding. Evaporation bypasses the reader entirely and is stepped (§3.1).

Time patterns are repeating multiplier sets — monthly, daily, hourly, and weekend-hourly — that modulate dry-weather sanitary inflows and external-inflow baselines. A dry-weather inflow holds up to four patterns, and the two halves of that arrangement disagree about what fixes a pattern’s role. Assignment is positional: the four columns of a [DWF] line fill a monthly, daily, hourly and weekend slot in order, with no check that the pattern named in a column declares the matching type. Evaluation, by contrast, dispatches on the pattern’s own declared type, ignoring the slot it occupies. The evident intent is that the two agree — the columns exist to be filled by patterns of the corresponding type — but nothing enforces it, so a mismatched pattern is neither rejected nor coerced: it simply contributes its declared type’s multiplier from whichever slot it sits in. The one case that silently contributes a multiplier of 1 is a weekend pattern placed in a weekday-active slot, since weekend dispatch yields a factor only on days 0 and 6.

Control rules are IF-THEN-ELSE statements over the simulation state (node depths, link flows, timing, …) with priorities, which switch pumps and adjust regulator settings; actions may be immediate or modulated through PID controllers (§11).

Beyond these, a model references shared parameter sets that are neither network elements nor data tables: transects and street sections (geometry shared by many conduits), aquifers (shared by many subcatchments), snow pack parameter sets, unit-hydrograph groups, and LID designs — each defined once and instantiated by reference.

1.7 The State-Vector View

SWMM is a distributed discrete-time simulation: at each time step the entire system state advances as $X_t = f(X_{t-1}, I_t, P)$ and outputs are computed as $Y_t = g(X_t, P)$, where $I_t$ are external inputs (precipitation, temperature, boundary stages, control settings) and $P$ the constant parameters. The state vector is remarkably small relative to the model’s scope. Per subcatchment: a ponded depth $d$ for each of the three sub-areas independently (§3.2), the infiltration state of the chosen method (e.g. cumulative infiltration, Horton-curve position, or moisture deficit), groundwater moisture content and saturated-zone depth, and snow pack depth/free-water/temperature/cold-content. Per conveyance node: water depth $y$. Per conduit: flow rate $q$ and flow area $a$. Per pollutant: surface buildup mass and ponded mass per subcatchment (with the last-swept time), and concentration per node and link. Everything else reported by the engine — velocities, volumes, flooding, loads — is derived from these states, the inputs, and the parameters.

The state vector closely tracks what a hotstart file checkpoints — though not exactly: the file also persists node lateral inflows, storage residence times, and regulator settings; stores link depth rather than flow area; and omits LID layer states entirely, so LID antecedent conditions are lost across a hotstart. The state vector likewise defines what initial conditions the user must supply, and marks the seam between the hydrologic and hydraulic halves of the engine, which advance on different clocks (§2).

2. Simulation Architecture

SWMM advances three clocks with independent step sizes: a runoff clock ($\Delta t_{roff}$, itself split into a wet step used while any precipitation, snow cover, surface runoff, or non-dry LID unit exists anywhere — a draining LID holds the short step long after rain ceases — and a much longer dry step used otherwise), a routing clock ($\Delta t_{rout}$, typically far shorter — seconds to a minute under dynamic wave), and a reporting clock ($\Delta t_{rpt}$). The main loop per routing step $[T,\ T + \Delta t_{rout}]$:

1. While the runoff clock lags the end of the routing step, compute hydrology for a full runoff step — precipitation, snowmelt, infiltration, evaporation, groundwater, overland flow, buildup/washoff — and advance it.
2. Route flow and quality through the conveyance network over the routing step, using the runoff results (computed on the coarser runoff grid) linearly interpolated to routing times as lateral node inflows.
3. If the reporting clock has been reached, interpolate results to the report time and write them to the output file — at most one report time is serviced per routing step.

Interpolation admits exceptions: infiltration and evaporation rates are held piecewise-constant within a runoff step (a report time inside the step receives the step-start value), groundwater elevation and soil moisture are reported at their end-of-step values, and climate-file temperatures are interpolated sinusoidally (§3.1).

Precipitation is the exception that does not merely skip interpolation — the rainfall written to the output is computed by a different rule from the rainfall that drove the runoff. The runoff side uses whatever intensity the gage holds for the current runoff step. The reporting side re-derives a value at each report time, advancing the report date by one second and then asking where it falls relative to the gage’s own interval boundaries:

$$i_{rpt} = \begin{cases} i_{current}, & t_{rpt} + 1\text{s} < t_{end},\ 0, & t_{end} \le t_{rpt} + 1\text{s} < t_{next},\ i_{next}, & t_{rpt} + 1\text{s} \ge t_{next}, \end{cases}$$

where $t_{end}$ closes the gage’s current recording interval and $t_{next}$ opens the following one. A report time landing in a gap between intervals therefore reports zero rainfall, and one landing at or past the next interval’s start reports the next interval’s intensity — a value the runoff computation for that step has not yet seen. A gage deferring to a co-gage copies its reported value, and an API rainfall override supersedes the whole rule. Co-gage status is decided at validation and short-circuits it: the first used gage already reading a given time series claims it, and any later gage naming the same series becomes its co-gage and returns at once — checked only for declaring the same rain type (fatal if not), and never subjected to the interval rules below. For a gage that does own its series, three checks apply, asymmetrically: the series may not itself be a reference to an external file (fatal); a declared recording interval longer than the series’ own smallest interval is fatal, since the gage would claim resolution the data cannot supply; while a declared interval shorter than it is merely a warning, the declared value being kept. This is why a reported rainfall series read back from the binary file can be offset by an interval from the one that produced the reported runoff. Both wet and dry runoff steps are truncated to end exactly at the next rainfall-interval boundary of any gage or the next evaporation-change date, so all forcing is constant within a step; a wet step longer than a used time-series-fed gage’s recording interval is permanently reduced to that interval with a warning (file-fed gages impose no such reduction). Guidance is wet step ≲ subcatchment time of concentration, dry step of hours to a day. Disabling a half of the engine re-times the loop rather than merely emptying it: with routing suppressed the step becomes the smaller of the wet runoff step and the report step instead of the routing step, and with runoff suppressed the climate state is still advanced each step so that evaporation continues to be applied.

Three option-driven behaviours modify the routing loop itself:

• An [EVENT] list restricts routing to date windows: between events the routing step stretches to the next runoff or report time, no lateral inflows are applied and no flow or quality routing occurs (hydrology continues; state freezes); overlapping events are clipped to the next event’s start.
• A steady-state skip option bypasses flow routing for a step when no control action fired, the previous step’s system flow error is within a tolerance (default 5%), and no node’s lateral inflow changed by more than a relative tolerance (also 5% by default; a zero↔nonzero change counts as 100%) — quality routing and outflow accounting still run.
• A rule step option evaluates control rules only at fixed intervals (the routing step is trimmed to land on them), though pump startup/shutoff depth targets still apply every step.

Lateral inflows are assembled at the start of each routing step — runoff, groundwater, and LID drains interpolated to the old routing time; external, dry-weather, and RDII inflows evaluated at the step-start date — with near-zero inflows truncated, and a negative external inflow legal and booked as an outflow removing mass at the node’s concentration.

The dual-clock design is an economy: hydrology on a 15-minute wet / 1-day dry grid costs a small fraction of the routing effort, and continuous multi-year simulations remain tractable while routing runs at whatever step stability demands (§5). The state vector of §1.7 is, up to the small deviations noted there, what is saved to and restored from a hotstart file, allowing a simulation to resume from a prior ending condition (e.g. to establish non-zero antecedent conditions before a design storm). Internally, all computation is in feet and seconds regardless of the user’s unit system, and dates are 8-byte doubles counting decimal days since 30 December 1899 (the Delphi epoch).

3. Hydrology

3.1 Meteorology

Precipitation enters through rain gages, from user time series or external files (NCDC formats, Environment Canada formats, and a standard user-prepared station format), recorded as intensity, volume, or cumulative volume on a fixed interval. User-supplied data are interpreted as start-of-interval values — end-of-interval records must be shifted back one interval — while NCDC and Canadian files carry end-of-interval stamps that SWMM converts automatically. External files are pre-collated into a binary rainfall interface file read during simulation. Radar or gridded rainfall is accommodated by one gage per grid cell or by area-weighting rainfall onto subcatchments.

Temperature (needed only for snowmelt and Hargreaves evaporation) comes from a time series, linearly interpolated, or from a daily climate file of max/min values which SWMM converts to instantaneous temperatures by sinusoidal interpolation. The minimum is placed at sunrise and the maximum three hours before sunset, both fixed by the solar declination on day $d$ of the year,

$$\delta = 0.40928\cos!\big(0.017202,(172 - d)\big),\qquad \omega_h = 3.8197,\arccos!\left(-\tan\delta,\tan\varphi\right),$$

$$h_{sr} = 12 - \omega_h + \Delta\lambda,\qquad h_{ss} = 12 + \omega_h + \Delta\lambda - 3,$$

where $\varphi$ is the site latitude and $\Delta\lambda$ a user longitude correction entered in minutes (conventionally four minutes per degree from the standard meridian). Half-sine arcs are then fitted between successive extremes, in three branches over the day:

$$T_a(h) = \begin{cases} T_{min} + \dfrac{T_{r}’}{2}\sin!\left(\dfrac{\pi,(h_{sr} - h)}{24 + h_{sr} - h_{ss}}\right), & h < h_{sr},\[10pt] T_{ave} + T_r\sin!\left(\dfrac{\pi,(h_{day} - h)}{h_{sr} - h_{ss}}\right), & h_{sr} \le h \le h_{ss},\[10pt] T_{max} - T_r\sin!\left(\dfrac{\pi,(h - h_{ss})}{24 + h_{sr} - h_{ss}}\right), & h > h_{ss}, \end{cases}$$

with $T_{ave} = (T_{min} + T_{max})/2$, $T_r = (T_{max} - T_{min})/2$, $h_{day} = (h_{sr} + h_{ss})/2$, and $T_r’$ the span from the previous day’s maximum to today’s minimum, so the overnight limb joins the two days continuously. Saturation vapour pressure, needed by the rain-on-snow melt equation, follows from the same temperature as

$$e_a = 8.1175\times10^{6},\exp!\left(\frac{-7701.544}{T_a + 405.0265}\right).$$

Evaporation applies to ponded water, groundwater, channels, storage units, and LIDs, from one of five sources: a constant, monthly averages, a time series (honouring each entry’s exact timestamp — values may vary within a day), climate-file daily values, or the Hargreaves formula. The time-series source is the one that is not interpolated: the engine holds the most recent entry’s rate until the clock passes the next timestamp, making the series a step function, where temperature and the other series interpolate. Hargreaves is computed from 7-day running averages of climate-file temperatures,

$$E = 0.0023,\frac{R_a}{\lambda},T_r^{1/2},(T_a + 17.8) \quad \text{mm/day},$$

where $T_a$ is the average temperature and $T_r$ the daily temperature range (both °C), $\lambda = 2.50 - 0.002361,T_a$ the latent heat of vaporisation, and $R_a$ the extraterrestrial radiation computed from latitude and day of year. Climate-file (pan) values are scaled by user monthly pan coefficients (≈0.7), and a DRY_ONLY switch suppresses evaporation from the land surface — subcatchment ponding and LID units — during rainfall, leaving conduit, storage-unit, and subsurface evaporation untouched. Wind speed — monthly averages (default zero) or climate-file daily values — enters only the rain-on-snow melt equation. Days missing from a climate file inherit the most recent recorded value of each variable.

Beyond the per-method infiltration patterns (§3.3), optional monthly [ADJUSTMENTS] modify the forcing itself, in four calendar-month vectors: additive offsets to temperature (rescaled by $9/5$ under SI, being a temperature difference) and to potential evaporation, and multiplicative factors on all gage rainfall (also applied during RDII preprocessing) and on saturated hydraulic conductivity (§3.3). Only the conductivity factor is range-guarded, and silently: a value $\le 0$ is replaced by $1$, so a month entered as zero means “no adjustment” rather than “no infiltration”. Per-subcatchment monthly patterns may further scale the pervious sub-area’s depression storage and Manning roughness; impervious sub-areas are never adjusted.

3.2 Surface Runoff

Each subcatchment sub-area is a nonlinear reservoir: an idealised rectangular plane of area $A$, characteristic width $W$, slope $S$, and Manning roughness $n$, holding a ponded depth $d$ with depression storage $d_s$. Mass balance and a Manning wide-channel rating for the outflow give the governing ODE

$$\frac{\partial d}{\partial t} = i - e - f - \alpha,(d - d_s)^{5/3}, \qquad \alpha = \frac{1.49,W\sqrt{S}}{A,n},$$

with $i$ the rainfall/snowmelt input, $e$ surface evaporation, $f$ infiltration (pervious sub-areas only), and outflow zero while $d \le d_s$. The wide-channel assumption sets hydraulic radius equal to $d - d_s$, whence the 5/3 exponent. Each of the three sub-areas (§1.2) integrates its own copy of this ODE — pervious and impervious sub-areas differ in roughness, depression storage, and infiltration, but share $W$ and $S$, with the impervious $\alpha$ prorated so both impervious sub-areas use $W/(A_2{+}A_3)$ — using the same adaptive fifth-order Runge–Kutta integrator as the groundwater module. The filling phase up to $d_s$ is handled analytically before the integrator engages. Subcatchment runoff is the area-weighted sum $\sum q_j A_j$. Two area caveats: run-on from upstream subcatchments and outfalls is spread over the non-LID area only, and snow packs likewise cover only the non-LID area (LID units receive raw precipitation directly). When snowmelt is simulated, the rain/snow split and catch factor apply gage-wide, so a subcatchment without a snow-pack object receives SCF-scaled snowfall as immediate liquid input.

All geometry is lumped into $\alpha$: the model has no internal spatial variation, so the width parameter is the primary shape calibration handle ($W \approx$ area / average maximum overland-flow length — the flow-path length to the drainage divide — with a skew correction for off-centre drainage; increasing $W$ sharpens and advances the hydrograph). This spatial uniformity is also what makes re-routing trivial: a fraction of impervious runoff may be directed onto the pervious sub-area, or (mutually exclusively) a fraction of pervious runoff onto the impervious sub-area with depression storage, and a subcatchment’s outflow onto another subcatchment — in each case applied like additional rainfall on the receiver, delayed by one time step per hop. Setting $n = 0$ bypasses the nonlinear routing — ponded water above depression storage converts instantly to runoff each step, though depression storage, evaporation, and infiltration still apply — which together with parameter choices lets the method emulate simple runoff-coefficient or SCS-volume models.

3.3 Infiltration

Infiltration is computed on the pervious sub-area of each subcatchment by one of five methods, selected per subcatchment. All methods share two conventions. First, the actual infiltration rate is the smaller of the potential (capacity) rate and the water available, $f = \min(f_p,\ i_a)$, where the available rate includes ponded water: $i_a = i + d/\Delta t$ with $d$ the current ponded depth (the Curve Number method excepted — since 5.2 it folds run-on into ponded depth only, so run-on infiltrates at the held rate but never advances the event’s cumulative-rainfall curve). Second, every method carries a recovery model that regenerates capacity during dry weather, so that continuous simulation across many storms is meaningful. Two distinct optional adjustments scale the constants: a monthly conductivity pattern (global or per-subcatchment) scales $f_0$, $f_\infty$, and $K_s$ — and the Green–Ampt upper-zone depth $L_u$ by its square root — while a separate monthly soil-recovery pattern scales every recovery/regeneration coefficient (and divides the Green–Ampt inter-event timer). The constants below are exact only absent both. Internally all quantities are in feet and seconds.

Horton

The classic exponential decay of capacity from an initial rate $f_0$ to an equilibrium rate $f_\infty$ ($\approx K_s$):

$$f_p = f_\infty + (f_0 - f_\infty)e^{-k_d t}$$

with decay coefficient $k_d$ (s⁻¹). SWMM does not evaluate this at wall-clock time: because actual infiltration can be rainfall-limited ($f < f_p$), it tracks an equivalent time $t_p$ on the Horton curve such that the curve’s cumulative infiltration matches what actually infiltrated. Cumulative capacity is

$$F(t_p) = f_\infty t_p + \frac{f_0 - f_\infty}{k_d}\left(1 - e^{-k_d t_p}\right)$$

The capacity actually applied over a step is not the point rate above but the step average of the cumulative curve, $f_p = \big(F(t_p + \Delta t) - F(t_p)\big)/\Delta t$, floored at $f_\infty$ — a distinction that matters whenever the step is long relative to $1/k_d$, since the point rate at the step’s start overstates what the curve delivers across it. Each wet step then either advances $t_p$ by $\Delta t$ (when infiltration proceeded at capacity, or the curve has flattened — SWMM treats $t_p \ge 16/k_d$ as flat) or solves $F(t_p^{old}) + f,\Delta t = F(t_p^{new})$ for $t_p^{new}$ by Newton–Raphson (when rainfall-limited). An optional cap $F_{max}$ makes the surface impermeable beyond a total infiltrated volume. During dry steps the state recovers along an exponential drying curve with coefficient $k_r$; the wetting- and drying-curve mapping collapses to the closed form

$$t_p \leftarrow -\frac{1}{k_d}\ln!\left[1 - e^{-k_r\Delta t}\left(1 - e^{-k_d t_p}\right)\right].$$

$k_r$ derives from a user drying time $T_{dry}$ (days) via $k_r = 3.912/T_{dry}$ (98% recovery definition). Recovery is purely empirical, independent of evaporation. When $F_{max}$ is active, the spent volume is tracked and wound back along the recovery curve during dry weather, so capacity under the cap regenerates.

Both Horton variants short-circuit on degenerate parameters, and one of the two outcomes is a trap. If $f_0 = f_\infty$ or $k_d = 0$ the method degrades to a constant capacity $f_0$, which is the sensible reading. But if $f_0 < f_\infty$ — the two rates entered the wrong way round, an easy slip — the method returns zero infiltration for the whole simulation, silently and without a validation error. Nothing distinguishes that case in the output from a genuinely impermeable surface.

Modified Horton

Akan’s reformulation replaces elapsed time as the state with cumulative excess infiltration $F_e$ — the volume infiltrated above the equilibrium rate — on the argument that only water accumulating near the surface reduces capacity, while $f_\infty$ percolates away harmlessly:

$$f_p = \max!\left(f_0 - k_d F_e,\ f_\infty\right), \qquad F_e = \sum_i \max(f_i - f_\infty,\ 0),\Delta t_i .$$

This behaves better under low-intensity rainfall (plain Horton decays capacity even when little water has actually entered the soil). The scheme is fully explicit — no Newton solve, and no step-averaging, the point rate being used directly — and dry-period recovery is a simple exponential decay of the state, $F_e \leftarrow F_e,e^{-k_r \Delta t}$. The degenerate-parameter behaviour above applies here unchanged, zero-infiltration trap included. Its $F_{max}$ handling is idiosyncratic: once active, each wet step sets $F_e \leftarrow \max(F_e, F_{max})$, so infiltration shuts off after a single wet step and only dry-weather decay restores it.

Green–Ampt

The Mein–Larson two-stage form of the sharp-wetting-front model. The soil is parameterised by saturated conductivity $K_s$, wetting-front suction head $\psi_s$, and an initial moisture deficit $\theta_d$; the engine adds the current ponded depth $d$ to the suction head throughout. Before the surface saturates, all rainfall infiltrates ($f = i_a$); saturation occurs once cumulative infiltration reaches

$$F_s = \frac{K_s,(\psi_s + d),\theta_d}{i_a - K_s}$$

(defined only while $i_a > K_s$). Thereafter capacity follows

$$f_p = K_s\left(1 + \frac{(\psi_s + d),\theta_d}{F}\right),$$

which SWMM integrates over the step in cumulative form — $F_2 = C + (\psi_s + d),\theta_d\ln!\big(F_2 + (\psi_s + d),\theta_d\big)$, solved for $F_2$ by Newton–Raphson — avoiding overshoot on long steps; sub-10-second steps with $F > 0.01,(\psi_s + d),\theta_d$ use the explicit point rate instead. Whichever branch runs, the result is floored at $F + K_s\Delta t$ and capped at the water actually available over the saturated part of the step, and a zero moisture deficit bypasses the solve for infiltration at exactly $K_s$. Recovery tracks the moisture deficit $\theta_{du}$ of an upper soil zone of fixed thickness $L_u = 4\sqrt{K_s}$ (inches, $K_s$ in in/hr): wet steps deplete the deficit by $f\Delta t/L_u$, dry steps regenerate it at rate $k_r,\theta_{dmax}$ with $k_r = \sqrt{K_s}/75$ hr⁻¹. After a dry spell longer than $T_r = 0.06/k_r$ a new event begins with $\theta_d = \theta_{du}$ and $F = 0$.

Modified Green–Ampt

A fifth selectable method (5.1.010) differing from Green–Ampt in exactly one respect: during low-intensity periods ($i_a \le K_s$) it does not reset the event state when the inter-event timer expires, so cumulative infiltration $F$ keeps building through light rain and surface saturation arrives sooner. This is also the variant invoked internally by LID surface layers and storage-node seepage.

Curve Number

An incremental adaptation of the SCS/NRCS relation $Q = P^2/(P + S_{max})$ with $S_{max} = 1000/CN - 10$ (inches). SWMM omits the usual initial-abstraction term (depression storage plays that role) and differences the cumulative form: each wet step updates event totals $P$ and $F = P - P^2/(P+S_e)$ and takes $f_p$ as $\Delta F/\Delta t$, where $S_e$ is the storage capacity remaining at the start of the current event. During rainless gaps within an event the previous rate is held so ponded water can continue to infiltrate. Remaining capacity $S$ depletes with infiltrated volume and recovers during dry weather at $k_r S_{max}$ per hour with $k_r = 1/(24,T_{dry})$; a dry spell longer than $0.06/k_r$ hours starts a new event with $S_e = S$. The curve number itself is clamped to $[10, 99]$. Because tabulated urban curve numbers already lump impervious cover, the method is documented as applying to a subcatchment treated as fully pervious.

3.4 Groundwater

Zone Structure and Fluxes

Each subcatchment may sit on an independent two-zone aquifer: an unsaturated upper zone of uniform moisture content $\theta$ and depth $d_U$, above a saturated lower zone of depth $d_L$ (the water-table height over the aquifer bottom), with $d_U = E_G - E_B - d_L$ for ground and bottom elevations $E_G, E_B$. The unknowns are $\theta$ and $d_L$. Six volumetric fluxes (per unit area) connect the zones to the surface and the conveyance system: surface infiltration $f_I$ (the §3.3 result scaled by pervious fraction, capped by upper-zone storability), upper-zone evapotranspiration $f_{EU}$, percolation $f_U$, lower-zone ET $f_{EL}$, deep percolation $f_L$, and lateral groundwater discharge $f_G$ to a designated conveyance node.

Moisture accounting reduces to a coupled ODE pair in $(\theta, d_L)$, driven by the zone flux sums $f_{UZ} = f_I - f_{EU} - f_U$ and $f_{LZ} = f_U - f_{EL} - f_L - f_G$, which SWMM integrates over each runoff time step with adaptive fifth-order Runge–Kutta, clamping $\theta \in [\theta_{WP}, \phi)$ and $d_L \in [0, E_G - E_B)$ — and jumping the water table to the surface whenever $\theta$ reaches porosity. A subcatchment’s [GROUNDWATER] line may override the shared aquifer’s bottom elevation, initial water table, initial moisture, and node threshold elevation.

Constitutive Relations

Percolation uses an exponential unsaturated-conductivity model with a finite suction-gradient factor,

$$f_U = K_s,e^{-(\phi-\theta)HCO}\left(1 + \frac{2,\psi_{TS},(\theta - \theta_{FC})}{d_U}\right)$$

where $\psi_{TS}$ is the aquifer’s tension-slope parameter, $HCO$ its conductivity slope, and $\phi$ its porosity; the flux is zero below field capacity $\theta_{FC}$ and capped at $d_U(\theta - \theta_{FC})/\Delta t$. The manual’s simpler $f_U = K_s e^{-(\phi-\theta)HCO}$ omits the gradient factor the code applies.

Evapotranspiration is drawn in priority order surface → upper zone → lower zone, the upper-zone share as a user fraction $UEF$ of potential ET (first prorated by the pervious area fraction — the only surface through which subsurface ET is exerted — and optionally rescaled by a per-aquifer monthly pattern) and the lower-zone share the complementary fraction $(1 - UEF)$ of that potential — complementary to the pattern-adjusted upper fraction, so the monthly pattern propagates into the lower zone inversely, raising the upper share depressing the lower one — itself scaled by how far the saturated zone reaches into the cutoff depth $DEL$ (declining linearly to zero once the water table falls past it) and then capped by whatever ET remains after the surface and upper-zone draws — with no subsurface ET at all during steps with surface infiltration, and no upper-zone ET at or below the wilting point. The surface draw netted off here is the runoff module’s whole evaporation total for the step — both impervious sub-areas and any LID units included — not the pervious share alone.

Deep percolation is a linear reservoir $f_L = DP, d_L/(E_G - E_B)$, capped at $d_L/\Delta t$.

Lateral Discharge to the Drainage System

The lateral discharge — the term that generates baseflow and groundwater infiltration into sewers — is the user-configurable power function

$$f_G = A1,(d_L - h^)^{B1} ;-; A2,(h_{SW} - h^)^{B2} ;+; A3,d_L,h_{SW}$$

where $h_{SW}$ is the surface-water stage at the receiving node (a fixed value or the live routed stage) and $h^$ a threshold height defaulting to the node invert. When $d_L \le h^$ the whole function returns zero — all three terms, so surface water cannot recharge a depleted aquifer through the $A2$ or $A3$ terms; the $A2$ term likewise vanishes on its own whenever $h_{SW} \le h^*$; a zero exponent degrades its term to the bare coefficient. Notably, the terms are evaluated in user length units and the result converted from user groundwater-flow units, making $A1/A2/A3$ unit-system-dependent whenever the exponents differ from 1 (custom expressions likewise run in user units). Choices of the five coefficients reproduce standard conceptualisations — a linear reservoir ($B1{=}1$, $A2{=}A3{=}0$), Dupuit–Forchheimer seepage, or Hooghoudt tile drainage — and negative $f_G$ (bank storage from channel to aquifer) is admitted when the interaction term is unused. The flux is bounded each step by what the aquifer stores, what the unsaturated zone can accept, and what the node can supply. User-defined expressions may customise both sinks — but asymmetrically: a deep-percolation expression replaces $f_L$, while a lateral-flow expression is added to the power-function $f_G$ (a pure replacement only if $A1 = A2 = A3 = 0$). The aquifer carries no water quality: infiltrate arrives at the node clean unless a constant concentration is assigned.

The ODE Pair in Manual and Source

The manual states the ODE pair inconsistently — the §5.2 derivation carries $(\phi - \theta)$ denominators where the §5.4 integration sidebar uses $\phi$, and the infiltration-cap formula flips the sign of $f_U$ between statements — and the pinned source implements a third form matching neither in full: $\partial d_L/\partial t = f_{LZ}/(\phi - \theta)$ (the §5.2 denominator), but $\partial\theta/\partial t = f_{UZ}/(E_G - E_B - d_L)$ with no $\theta f_{LZ}$ coupling term, and an infiltration cap $(E_G - E_B - d_L)(\phi - \theta)/(F_{perv},\Delta t)$ containing no $f_U$ term of either sign. The code’s forms, not either manual variant, are SWMM’s actual behaviour.

3.5 Snowmelt

Snow state is kept as depth of water equivalent and simulated per subcatchment on a three-way split that differs from the runoff sub-areas: pervious (SA1), plowable impervious (a user fraction $SNN$ of impervious area — streets and lots subject to snow removal, always fully snow-covered), and remaining impervious (SA3, rooftops). Precipitation falls as snow when air temperature $T_a \le SNOTMP$, with gage snowfall scaled by a snow catch factor $SCF$ to correct wind under-catch.

Melt is computed per surface by two regimes. During rain ($i > 0.02$ in/hr), Anderson’s energy-budget equation for saturated, radiation-free conditions applies:

$$SMELT = \left(0.001167 + 7.5\gamma U_A + 0.007,i\right)(T_a - 32) + 8.5,U_A(e_a - 0.18)$$

(in/hr; $U_A = 0.006,u$ for wind speed $u$ in mph; $\gamma = 0.000359,p_a$ with atmospheric pressure $p_a$ computed from the site’s average elevation; $e_a$ saturation vapour pressure at $T_a$). Otherwise, when $T_a \ge T_{base}$, a degree-day law $SMELT = DHM,(T_a - T_{base})$ applies, with the melt coefficient varying sinusoidally through the year between a December 21 minimum $DHMIN$ and June 21 maximum $DHMAX$; below $T_{base}$ no melt occurs and the step instead updates the cold-content account. Street de-icing is represented by lowering a surface’s $T_{base}$ rather than by explicit chemistry.

Two mechanisms delay and shape melt. The degree-day coefficient itself is not constant: it sweeps sinusoidally through the year between the user’s December and June limits,

$$DHM(d) = \tfrac{1}{2}\Big[DHM_{max}(1 + s) + DHM_{min}(1 - s)\Big],\qquad s = \sin!\big(0.0172615,(d - 81)\big),$$

so $s = +1$ near day 172 (21 June) and $-1$ near day 355 (21 December).

A cold content account $cc$ (heat deficit, in water-equivalent depth) must then be paid off before any liquid melt leaves. Its antecedent temperature index snaps to the air temperature during snowfall heavier than 0.02 in/hr and otherwise relaxes toward it,

$$ATI \leftarrow ATI + \gamma,(T_a - ATI),\qquad \gamma = 1 - (1 - TIPM)^{\Delta t/6,\text{hr}},$$

which rescales the user’s $TIPM$ weight from its nominal 6-hour basis to the actual step; $ATI$ is then capped at the surface’s base melt temperature. Below that temperature the deficit accumulates rather than melting,

$$cc \leftarrow \min\Big(cc + RNM\cdot DHM,(ATI - T_a),\Delta t,ASC,\ \ 0.007,\tfrac{1}{12},W_{snow},(T_{base} - ATI)\Big),$$

the cap expressing an assumed snow specific heat of 0.007 in w.e. per °F per foot of pack (equivalently a dimensionless $0.007/12$ per °F when deficit and pack are measured in the same unit), and the negative-melt ratio $RNM$ scaling the exchange rate. When melt does occur it is debited against the deficit first: a potential melt $SMELT$ over a step with conversion factor $\Delta t,RNM,ASC$ either clears the deficit and passes the remainder on, or is absorbed entirely and reduces $cc$ by $SMELT$ times that factor. A free-water reservoir additionally requires the pack’s liquid-holding capacity ($FWFRAC \times$ pack depth) to fill before runoff releases; the initial free-water depth is silently clamped to that same product of the entered free-water fraction and initial pack depth, so an over-specified initial value is reduced rather than refused. Partial snow cover on SA1 and SA3 is handled by areal depletion curves (fraction of area snow-covered vs. relative pack depth $WSNOW/SI$), two watershed-wide curves with Anderson’s temporary-linear-curve adjustment after fresh snowfall on partial cover — 100% cover is assumed until 25% of the new snow melts ($SBWS = AWE + 0.75,SNO/SI$); melt and cold-content exchange scale by the covered fraction $ASC$. Once the plowable surface’s depth reaches a trigger $WEPLOW$ — defaulted to $10^6$, so ploughing is off unless the user supplies a depth — its entire current depth is redistributed by five constant fractions — to the other sub-areas, to another subcatchment, out of the system, or to immediate melt.

The net result per surface, $RI = ASC \cdot SMELT + (1-ASC),i$ plus any immediate melt, replaces gage rainfall as the input to infiltration and overland flow — though the melt term is not a direct pass-through: rain falling on the snow-covered fraction is added to the free-water reservoir alongside the melt and released only once the holding capacity is exceeded, so rain-on-snow is itself delayed by the pack; the two impervious sub-area results are area-averaged onto the runoff model’s impervious sub-areas, while the pervious result carries over directly. When a pack thins below 0.001 in it is flushed as immediate melt. Snow is assumed not to alter infiltration or surface roughness.

3.6 Rainfall-Dependent Inflow/Infiltration (RDII)

RDII — stormwater entering sanitary and combined sewers through defects and illicit connections — is modelled independently of the runoff/groundwater machinery, as a rainfall-convolved inflow at designated nodes. The kernel is the RTK triangular unit hydrograph, defined by $R$ the fraction of rainfall volume entering the sewer, $T$ the time to peak, and $K$ the recession-to-peak ratio. Its base is $T_b = T(1 + K)$ and, normalised to unit rainfall depth over unit area, its ordinate at time $t$ is the triangle

$$u(t) = \frac{2}{T_b}\cdot\begin{cases} t/T, & 0 \le t \le T,\[2pt] 1 - \dfrac{t - T}{T_b - T}, & T < t < T_b,\[4pt] 0, & t \ge T_b, \end{cases}$$

so that $\int_0^{T_b} u,dt = 1$ and the peak ordinate is $2/T_b$. The engine carries the ordinate per hour while holding $T$ and $T_b$ in seconds, so its stored peak is $2\times3600/T_b$ and the convolution below yields a depth per hour, converted to a flow by the sewershed area. The RDII flow at any instant is the discrete convolution of this kernel with the depleted past rainfall,

$$q_{rdii}(t) = \sum_{k=1}^{3}\ \sum_{i} R_{m(i),k}\ u_k!\left(\left(p_i - \tfrac{1}{2}\right)\Delta t_r\right) v_i,$$

where $v_i$ is the net rain depth in past interval $i$, $p_i$ its age in processing intervals, $\Delta t_r$ the interval, and $m(i)$ the calendar month that interval fell in — so both the ordinate and the $R$ ratio are taken from the month of the rainfall, not of the response. The kernel is evaluated at interval mid-points.

Each triangle depletes its own initial abstraction store before contributing: with $IA_{max}$ the capacity and $IA_{used}$ the amount already spent, an interval of gross depth $P_i$ yields

$$v_i = \max!\left(P_i - \left(IA_{max} - IA_{used}\right),\ 0\right),$$

and the store advances by whatever was absorbed, recovering at rate $IA_r$ through dry intervals. Because observed RDII responses are multi-modal, each unit-hydrograph group sums up to three triangles of increasing duration — rapid inflow, mixed, slow infiltration — and each group may vary by calendar month. Each triangle carries an initial abstraction account ($IA_{max}$, initial depletion $IA_0$, recovery rate $IA_r$) that absorbs rainfall before convolution and regenerates in dry weather.

RDII flows are computed for the whole simulation before routing begins and written to an interface file (unless an IGNORE_RDII option suppresses the subsystem): per node, gage rainfall — with any monthly rainfall adjustment applied — is sampled onto a processing grid set to the minimum of the wet runoff step and the shortest rising or falling limb across all months and all three unit hydrographs, depleted by initial abstraction, and convolved; results are emitted at the wet runoff step and held piecewise-constant during routing, with flows below 0.0001 cfs zeroed. Monthly parameters are selected by the month each rainfall increment fell in, not the month of the response. The per-area result is scaled by a user sewershed area (which need not correspond to any subcatchment — RDII-only models are common). Each month’s three $R$ values must individually be non-negative and sum to at most 1 — enforced with a 1% slack, so a month summing to 1.01 is accepted. The R-T-K parameters have no meaningful defaults; they are calibrated against flow-monitor records with dry-weather flow subtracted.

4. Flow Routing Theory

Conveyance routing solves the one-dimensional Saint-Venant equations — continuity and momentum for gradually-varied unsteady free-surface flow —

$$\frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = 0, \qquad \frac{\partial Q}{\partial t} + \frac{\partial (Q^2/A)}{\partial x} + gA\frac{\partial H}{\partial x} + gA,S_f = 0,$$

with $A$ flow area, $Q$ flow, $H = Z + Y$ hydraulic head, and friction slope from Manning: $S_f = (n/1.486)^2,Q|U|,/,(A R^{4/3})$ (the $|U|$ making friction oppose the flow direction). SWMM offers three levels of approximation:

• Steady flow routing simply translates each conduit’s inflow hydrograph to its outlet within the step — no storage, delay, or attenuation, though evaporation/seepage losses are first subtracted — with flow area back-computed from the Manning rating and flow capped at conduit capacity. It shares kinematic wave’s topology restrictions and serves for screening and preliminary sizing.
• Kinematic wave keeps continuity but reduces momentum to $S_0 = S_f$: flow is always at Manning normal depth, $Q = \beta,\Psi(A)$ with $\beta = 1.486\sqrt{S_0}/n$ and section factor $\Psi = A R^{2/3}$. Hydrographs translate and attenuate through conduit storage, but backwater, reversal, pressurisation, and entrance/exit losses are unrepresentable. A conduit’s accepted inflow is capped at its full-flow capacity, the rejected excess remaining at the upstream node as flooding or ponding, and any node with storage limits its outflow to inflow plus stored volume per step. The network must be a directed acyclic graph with junctions limited to one outlet link of any type (storage nodes exempt — they may have several), no adverse-slope conduits (a user minimum slope rectifies them here, with a warning), and regulators permitted only as outlets of storage nodes — which also always discharge freely, taking the upstream node’s own invert as tailwater, never submerged; dividers function only here, splitting flow by their §1.3 rules. Non-storage node depths are reconstructed as the maximum over connecting conduits of end depth plus offset, capped at full depth.
• Dynamic wave solves the full pair over the general network graph — loops, multiple outfalls, backwater, reverse flow, surcharge — and is the production method (§5). Its own topology demands: at least one outfall must exist, and a dummy conduit or ideal pump must be the sole link leaving its upstream node (with no dummy link leaving a node fed only by dummy links or ideal pumps, and no storage node a dummy outflow).

At validation SWMM silently normalises geometry that would otherwise be inconsistent: a node’s maximum depth is raised (with a warning) to the crown of its highest connecting link — pumps and bottom orifices exempted entirely, the downstream node raised only by conduits, and storage nodes skipped unless they carry a surcharge depth — and a regulator whose crest sits below its downstream node’s invert has the crest raised to that invert under dynamic wave.

Kinematic wave’s numerical scheme is a weighted implicit (Wendroff) four-point difference of continuity over each conduit, with both space and time weights fixed at 0.6 — unconditionally stable for weights above 0.5, so no Courant restriction applies. Conduits are processed in topological order; at each, the known upstream flow yields a scalar nonlinear equation in downstream area, $\beta\Psi(A_2) + C_1 A_2 + C_2 = 0$, solved by bracketed Newton–Raphson to 0.1% of full area. Storage nodes (in both steady and kinematic modes) iterate a trapezoidal mass balance against their head-dependent outflow rating with under-relaxation 0.55 and tolerance 0.005 ft, executing at most 9 balance passes (a 10-cap loop counting from 1). Junction flooding sheds any net inflow surplus as overflow (optionally banked as a ponded volume re-injected as capacity recovers).

5. Dynamic Wave Analysis

The dynamic-wave engine is a staggered node-link scheme: conduits carry the momentum equation for flow, nodes carry continuity for head, and the two are advanced together by fixed-point (Picard) iteration within each time step.

5.1 Conduit Flow Update

Substituting continuity into momentum and discretising over a conduit of length $L$ (implicit backward Euler in time, end-difference in space, overbars denoting conduit-average values) gives the update SWMM actually computes:

$$Q^{t+\Delta t} = \frac{Q^t + \Delta Q_{inertia} + \Delta Q_{pressure} + \Delta Q_{loss}}{1 + \Delta Q_{friction} + \Delta Q_{losses}}$$

$$\Delta Q_{inertia} = \sigma\left[2\bar U(\bar A^{t+\Delta t} - \bar A^{t}) + \bar U^2\frac{(A_2 - A_1)\Delta t}{L}\right],\quad \Delta Q_{pressure} = -g A_w\frac{(H_2 - H_1)\Delta t}{L},\quad \Delta Q_{friction} = \frac{g,(n/1.486)^2,|\bar U|,\Delta t}{R_w^{4/3}}$$

where $\bar U$ and $\bar A$ are the mid-section velocity and area, $\bar A^{t}$ the mid-section area at the previous time step (not the previous iterate), $A_1, A_2, H_1, H_2$ the end areas and heads, and $A_w, R_w$ the upstream-weighted area and hydraulic radius defined below — the pressure and friction terms alone use the weighted values, every other term the mid-section ones. $\Delta Q_{loss}$ is the evaporation/seepage momentum term $2.5,\bar U q_L \Delta t/L$ of §8, which uses the conduit’s true length where the other terms use its lengthened one.

Friction (and entrance/exit/average local losses, treated likewise) sit in the denominator — an implicit linearisation that keeps the update stable as flows approach zero. The inertial damping factor $\sigma$ scales the inertial terms by the Froude number ($\sigma = 1$ for $Fr \le 0.5$, tapering linearly to $0$ at $Fr = 1$), suppressing the terms that destabilise trans- and supercritical flow; user options force $\sigma = 1$ (keep all inertia) or $\sigma = 0$ (the local-inertial formulation — distinct from the diffusion wave, which also drops $\partial Q/\partial t$), and closed conduits flowing full always use $\sigma = 0$. The upstream weighting of the pressure/friction areas is this same Froude-based $\sigma$ (computed before any user damping override): none at $Fr \le 0.5$, fully upstream at $Fr \ge 1$, applied only in positive, non-full, downstream-sloping flow.

Several limits then constrain the updated flow:

• The velocity used in forming the momentum terms (not the resulting flow) is capped at 50 ft/s.
• A flow that reverses sign between successive iterates is clamped to 0.001 cfs in the new direction.
• Flow out of an essentially dry node is clamped to $\pm10^{-4}$ cfs rather than zeroed.
• A user conduit flow limit, when given, caps $|Q|$ every iteration.
• A positive computed flow is limited to Manning normal flow when the water-surface slope is less than the bed slope (upstream depth below downstream) or the upstream Froude number is at least 1 — user-selectable criteria (slope, Froude, both, or neither, disabling the limit entirely), except that conduits adjoining an outfall always apply the slope test and never the Froude test; the check is skipped for full upstream ends, critical/dry flow classes, and culvert-coded conduits.

Flow classification precedes all of this, and a conduit whose both end depths reach full depth is assigned the subcritical class outright without being classified — the branch that leads to the surcharge treatment of §5.4. Special flow classes at nearly-dry or critical-depth ends substitute critical/normal depth for the nodal head on the affected end — with a linear fasnh ramp of the downstream area contribution across the band between critical and normal depth — and a conduit classed dry at both ends (or at either end alone) carries exactly zero flow for the trial while retaining a nominal $\partial Q/\partial H$. Multi-barrel conduits solve one barrel and scale back (§1.4).

5.2 Node Head Update

Where §5.1 carries momentum on the links, the nodes carry continuity. Each node integrates

$$\frac{\partial H}{\partial t} = \frac{\sum Q}{A_S}$$

discretised over the step with trapezoidal averaging of the net inflow,

$$H^{t+\Delta t} = H^{t} + \frac{\Delta V}{A_S},\qquad \Delta V = \tfrac{1}{2}\left[\left(\sum Q\right)^{t} + \left(\sum Q\right)^{t+\Delta t}\right]\Delta t,$$

where $\sum Q$ is the node’s net inflow — link flows signed by direction, plus lateral inflow, less evaporation and seepage — and $A_S$ is an assembled surface area rather than a property of the node.

That assembly is the substance of the scheme. $A_S$ sums the node’s own storage area (zero for a junction, and replaced by the user’s ponded area once a ponding-enabled node rises above its full depth) with a contribution from each connecting link. A conduit in ordinary subcritical flow gives each end node the width-weighted trapezoid of its adjacent half,

$$A_{S,1} = \frac{\big(W(y_1) + W(\bar y)\big)L}{4},\qquad A_{S,2} = \frac{\big(W(\bar y) + W(y_2)\big)L}{4},\lambda,$$

with $\bar y$ the mid-point depth and $\lambda$ the fasnh ramp of §5.1 — note this is a width-weighted average over the half-length, not half the conduit’s plan area. Flow class then reapportions:

• a critical (free-fall) end contributes nothing, and the far node instead takes $\big(W(\bar y) + W(y)\big)L/2$ — the average width over the full length;
• a dry end contributes only when the conduit has no offset there, the wet end taking its ordinary half;
• a conduit dry at both ends contributes a nominal $10^{-4}L/2$ to each.

Closed-conduit top widths are frozen at 96% of full depth (98.53% under the slot) so the contribution never collapses to zero at the crown. Weirs and outlets contribute no surface area at all; orifices contribute half of theirs to each end — the equivalent-pipe water surface for a side orifice, the bare opening area for a bottom one — dropped at storage-node and critical-class ends. Multi-barrel conduits scale their contribution by the barrel count. The total is finally floored at a user-adjustable minimum defaulting to 12.566 ft², the plan area of a 4-ft manhole, which is what keeps $\partial H/\partial t$ finite at a node with no storage and no wet links.

5.3 Picard Iteration

Within a time step, each trial runs in two phases: all true conduits solve from the last-iteration heads (an OpenMP-parallel loop under the THREADS option — nodal accumulation stays serial, so results are thread-count-invariant), then dummy conduits, pumps, and regulators solve serially in link-definition order, each immediately updating its node flows — so a pump’s available-volume clamp sees the accumulation so far, making pump/regulator results sensitive to link order. Both flows and heads are then under-relaxed against the previous iterate,

$$Q^{(m)} \leftarrow (1-\omega),Q^{(m-1)} + \omega,\tilde{Q}^{(m)},\qquad H^{(m)} \leftarrow (1-\omega),H^{(m-1)} + \omega,\tilde{H}^{(m)},\qquad \omega = 0.5,$$

where $\tilde{\cdot}$ denotes the raw result of the flow or head update and $m$ the trial index. Pumps are exempt from the flow relaxation and surcharged nodes from the head relaxation, both being solved by mechanisms that already carry their own damping. The relaxation is applied first and the result then tested: if the relaxed flow has the opposite sign to the previous iterate it is discarded and replaced by $\pm10^{-3}$ cfs in the new direction, so a reversal must pass through approximately zero on some iteration rather than jumping across it. Iteration continues until every non-outfall node satisfies $|H^{(m)} - H^{(m-1)}| \le$ the head tolerance (default 0.005 ft) — outfalls are excluded, their depths reset from the boundary condition each pass — with a hard minimum of 2 trials and a user-adjustable maximum defaulting to 8 (links between converged nodes are bypassed). Non-convergence is tallied and reported but does not halt the simulation.

5.4 Surcharge

A non-storage node is surcharged when its head exceeds the crown of its highest connecting link — orifice and weir opening tops participate, not conduits alone; closed storage nodes surcharge only when a supplementary surcharge depth is specified and the full depth exceeded, and ponded nodes never do. Two treatments exist:

• The classic EXTRAN point iteration. With no free surface there is no $A_S$ to divide by: continuity degenerates to the algebraic constraint $\sum Q = 0$, and the head is corrected by a Newton step on it,

  $$\Delta H = \frac{\kappa \sum Q}{D},\qquad D = \sum_{\text{links}}\frac{\partial Q}{\partial H},$$

  where $\sum Q$ is the net inflow as in §5.2 and $\partial Q/\partial H$ for each link falls out of the flow-update denominator of §5.1 as $g A_w\Delta t,b/(L,\text{denom})$ for $b$ barrels, positive by construction and accumulated as a positive magnitude at both end nodes. Regulators and pumps supply their own derivative (§7) on the same positive-magnitude convention, with one exception: a Type 4 pump, whose flow is a function of its inlet depth alone rather than of a head difference, contributes its derivative to the upstream node only. The stepwise pump types contribute nothing, having no derivative to give. The sign convention repays attention: the Newton step for $f(H) = \sum Q = 0$ is $-f/f’$, and $\partial(\sum Q)/\partial H$ is negative because raising a node’s head drives flow out of it — so the two negatives cancel and the correction as written is added to the head, net inflow raising it. The two conventions are individually consistent — signed derivatives with the leading minus, or positive magnitudes without it — and SWMM uses the second throughout; it is only their mixture that inverts the correction. The damping factor $\kappa$ is 0.6 at terminal upstream nodes — those with outflow links but no inflow links — and 1 elsewhere. Rather than switch discontinuously between this and the free-surface update, SWMM 5.2 blends the two over a transition band reaching 25% above the crown: writing $f_H = (H - H_{crown})/H_{crown}$ for the relative surcharge depth, the denominator becomes

  $$D = \sum\frac{\partial Q}{\partial H} + \left(\frac{A_S^{,*}}{\Delta t} - \sum\frac{\partial Q}{\partial H}\right)e^{-15 f_H}, \qquad f_H < 0.25,$$

  with $A_S^{,}$ the surface area saved from the node’s last non-surcharged state. At $f_H = 0$ the weighting is unity and $D$ reduces to $A_S^{,}/\Delta t$ — exactly the free-surface update — while by $f_H = 0.25$ the exponential has decayed to $e^{-3.75} \approx 0.024$ and the pure Newton form governs. The corrected head is floored at the crown.

• The Preissmann slot method (an optional SURCHARGE_METHOD introduced in 5.1.013 — EXTRAN remains the default): closed conduits acquire a narrow hypothetical slot above the crown — width $0.5423,e^{-(y/y_{full})^{2.4}}$ of the maximum width (Sjöberg’s formula), switching to a flat 1% above $y = 1.78,y_{full}$ where the exponential has itself decayed to almost exactly that value — so depth may exceed the crown and the ordinary free-surface equations remain valid everywhere; hydraulic radius freezes at its full-pipe value, and the special surcharge branch is never taken.

5.5 Flooding and Ponding

A non-ponded node whose head would exceed its ground (plus optional surcharge depth) is pinned there and the surplus inflow is lost as reported flooding; with ponding enabled, the surplus accumulates in a user-specified ponded area atop the node — a virtual storage whose head may rise above ground and which drains back as the system recovers.

5.6 Time-Step Control

With a variable step enabled (the 5.2.4 default — Courant factor 0.75, with PARTIAL inertial damping the companion default), the step is the minimum over conduits of the Courant time $\frac{L}{|U| + \sqrt{gA/W}}$ (computed as $\tfrac{L}{|U|}\cdot\tfrac{Fr}{1+Fr}$, which is identically that expression since $Fr = |U|/\sqrt{gA/W}$, and scaled by the user’s Courant factor, with $L$ the lengthened length of the conduit-lengthening transform below rather than the true length; conduits with $Fr \le 0.01$ or negligible flow or area exempt — and because the Froude number is defined as zero for a full closed conduit, surcharged pipes are exempt as well, under either surcharge method, so the Courant condition never constrains the step on a pressurised network) and over nodes of the time to change head by a quarter of the crown height at the recent rate (outfalls, near-dry, and surcharged nodes exempt), bounded above by the user’s fixed routing step, floored at a minimum step, quantised down to a whole millisecond, and starting the run at the minimum step. Because the scheme is iterative and semi-implicit, Courant factors above 1 are usable. The optional conduit lengthening transform trades short conduits for stability: $L’ = \max[L,\ \Delta t(\sqrt{g y_{full}} + U_{full})]$ with slope rescaled by $L/L’$ and roughness by $\sqrt{L/L’}$ — preserving the conveyance factor $\beta$ exactly (the manual’s own $\sqrt{}$-rescale of slope contradicts the code). The rule of thumb for stability is $\Delta t \approx L/\sqrt{g,y_{full}}$; the standard diagnostics are the continuity error, a per-link flow-instability index, and a capacity-limited flag raised when a conduit’s upstream end is full with HGL slope exceeding the bed slope (§12).

5.7 Initial Conditions

Default zero depths and flows; user-supplied initial conduit flows imply Manning normal depth; a node without a user depth is seeded with the average, over every link touching it, of that link’s depth plus the link’s upstream offset — the same offset applied at both ends — and only at non-outfall, non-storage nodes; conduits without an initial flow then take the mean of their two end-node depths. A hotstart file bypasses these depth-seeding heuristics (derived volumes and areas are still computed from its state).

6. Cross-Section Geometry

6.1 Shape Families

Every conduit shape must supply a consistent family of geometric functions — area $A(Y)$, top width $W(Y)$, hydraulic radius $R(Y)$, the inverses $Y(A)$ and $A(\Psi)$, and the section factor $\Psi(A) = A,R(A)^{2/3}$ with its derivative — because the routing methods consume geometry only through these. Everything downstream reduces to this interface: kinematic wave solves $\beta\Psi(A) + C_1A + C_2 = 0$ for area (§4), dynamic wave reads $A$, $R$ and $W$ each iteration (§5), and normal depth inverts $\Psi$ (§6.3). Three implementation families cover the shape library.

Analytic Shapes

The wetted perimeter is the quantity that distinguishes these; $R = A/P$ throughout, and the section-factor derivative is

$$\Psi’(A) = \left(\frac{5}{3} - \frac{2}{3},\frac{dP}{dA},R\right)R^{2/3}$$

wherever $dP/dA$ is available in closed form, and a central difference over $\pm 0.001,A_{full}$ otherwise.

For a closed rectangle of width $W$, $A = WY$, $W(Y) = W$ and $P = W + 2A/W$, with $dP/dA = 2/W$. Above the section-factor peak at $\alpha_{max} = 0.97$ (where $\alpha = A/A_{full}$) the crown is progressively enrolled into the perimeter,

$$P = W + \frac{2A}{W} + \frac{\alpha - \alpha_{max}}{1 - \alpha_{max}},W,$$

which is what makes $\Psi$ fall back to $\Psi_{full}$ as the pipe fills. An open rectangle instead subtracts the walls the user declares frictionless: with $n_{ig} \in {0,1,2}$ walls ignored, $P = W + (2 - n_{ig})Y$ and $dP/dA = (2 - n_{ig})/W$.

A trapezoid of base width $b$ and side slopes $z_1, z_2$ (horizontal per unit vertical) is described through $s = (z_1 + z_2)/2$ and $r = \sqrt{1 + z_1^2} + \sqrt{1 + z_2^2}$:

$$A(Y) = (b + sY),Y,\qquad W(Y) = b + 2sY,\qquad P(Y) = b + rY,$$

$$Y(A) = \frac{\sqrt{b^2 + 4sA} - b}{2s}\ \ (s > 0),\qquad \frac{dP}{dA} = \frac{r}{\sqrt{b^2 + 4sA}},$$

degenerating to $Y = A/b$ for a rectangle-like $s = 0$. A triangle of maximum width $W_m$ is the same algebra with $b = 0$, written through the side slope $s = W_m/(2Y_{full})$ and $r = \sqrt{1 + s^2}$:

$$A(Y) = sY^2,\quad W(Y) = 2sY,\quad R(Y) = \frac{sY}{2r},\quad Y(A) = \sqrt{A/s},\quad \frac{dP}{dA} = \frac{r}{\sqrt{As}}.$$

A parabola whose half-width is $x = c\sqrt{Y}$, fixed by $c = W_m/(2\sqrt{Y_{full}})$, has

$$A(Y) = \tfrac{4}{3},c,Y^{3/2},\qquad W(Y) = 2c\sqrt{Y},\qquad Y(A) = \left(\frac{3A}{4c}\right)^{2/3},$$

and the only shape in the library with an analytic arc-length perimeter,

$$P(Y) = \frac{c^2}{2}\left[,x\sqrt{1 + x^2} + \ln!\left(x + \sqrt{1 + x^2}\right)\right],\qquad x = \frac{2\sqrt{Y}}{c}.$$

A power-law section takes a user exponent $n$, stored inverted as $m = 1/n$, with $c = W_m/[(m+1)Y_{full}^{,m}]$:

$$A(Y) = c,Y^{m+1},\qquad W(Y) = (m+1),c,Y^{m},\qquad Y(A) = (A/c)^{1/(m+1)}.$$

Its perimeter has no closed form and is summed as polyline arc length in depth increments of $0.02,Y_{full}$ — the one shape whose $P$ is quadrature rather than formula, and whose $\Psi’$ therefore falls to the central difference.

Circular Sections

Circles are parameterised by the subtended angle $\theta$ rather than by depth, giving the exact relations

$$\frac{A}{A_{full}} = \frac{\theta - \sin\theta}{2\pi},\qquad \frac{Y}{Y_{full}} = \frac{1 - \cos(\theta/2)}{2},\qquad P = \frac{\theta,Y_{full}}{2},$$

$$\frac{\Psi}{\Psi_{full}} = \frac{(\theta - \sin\theta)^{5/3}}{2\pi,\theta^{2/3}},$$

with $A_{full} = \pi Y_{full}^2/4$ and $R_{full} = Y_{full}/4$. Recovering $\theta$ from a known area means solving $\theta - \sin\theta = 2\pi\alpha$, which SWMM does by Newton iteration — at most 40 passes to a $10^{-4}$ tolerance, with a correction limiter to stop the near-full end diverging — seeded from

$$\theta_0 = \begin{cases} 0.031715 - 12.79384,\alpha + 8.28479\sqrt{\alpha}, & \alpha \le 0.04,\[2pt] 1.2 + 5.08,(\alpha - 0.04)/0.96, & \alpha > 0.04. \end{cases}$$

That limiter is worth stating in full, because what it does and what it was written to do are not the same. The correction is truncated only when it exceeds $+1$, and the truncated value is formed by a sign-transfer helper returning unit magnitude carrying the correction’s own sign. Under a guard that has already established the correction is positive, that helper can only ever return $+1$: the sign transfer is inert as written, and the same expression would be produced by the literal. A sign-transfer idiom is what one writes for a symmetric magnitude clamp, $|\Delta\theta| \le 1$, and it is used to genuine effect elsewhere in the library’s own root-finder. The evident intent is therefore to bound the correction’s magnitude; the behaviour bounds it in one direction only, so an arbitrarily large negative correction passes through untouched. Whether that asymmetry ever changes a converged answer depends on the seed, which for $\alpha \le 0.04$ is a fitted polynomial that can overshoot; the guard’s own comment attributes it to convergence at large $\theta$.

Below $\alpha = 10^{-5}$ the iteration is skipped for the small-angle asymptote $\theta = (37.6911,\alpha)^{1/3}$, whence $Y/Y_{full} = \theta^2/16$ and $\Psi/\Psi_{full} = \theta^{13/3}/124.4797$. The inverse direction — area from section factor — mirrors this: for $\psi = \Psi/\Psi_{full} \le 0.015$, $\theta = (124.4797,\psi)^{3/13}$ and $\alpha = \theta^3/37.6911$; above it a second Newton solve runs on the section-factor relation, seeded piecewise in four bands,

$$\theta_0 = \begin{cases} 0.12103 - 55.5075,\psi + 15.62254\sqrt{\psi}, & \psi \le 0.015,\ 1.2 + 1.94,(\psi - 0.015)/0.485, & 0.015 < \psi \le 0.5,\ 3.14 + 1.03,(\psi - 0.5)/0.4, & 0.5 < \psi \le 0.90,\ 4.17 + 1.12,(\psi - 0.90)/0.176, & \psi > 0.90, \end{cases}$$

and iterated to the same $10^{-4}$ tolerance over at most 40 passes — but stabilised differently, with no correction clamp at all, instead taking the absolute value of the iterate at the head of each pass so a negative excursion reflects back rather than being truncated. Both solvers share a non-convergence policy that is deliberate rather than accidental: each saves its seed before iterating and returns that seed if all 40 passes pass without meeting the tolerance. A non-converged result is therefore the piecewise polynomial estimate above, not the last iterate — the wandering iterate is judged less trustworthy than the fit it started from. Only beyond $\alpha = 0.04$ does the circle revert to table lookup, so the near-empty regime every dynamic-wave step visits is analytic.

Two circular variants reuse this machinery. A filled circular section models sediment of depth $Y_b$ by evaluating the full circle and subtracting the buried segment, then repairing the perimeter — the buried arc leaves and the sediment surface chord $w_b$ arrives, so $R_{full} = A_{full}/(\pi Y_{full} - P_b + w_b)$. Depth is thereafter measured from the sediment surface rather than the invert, and a sediment depth at or above the diameter is fatal. Its section-factor peak stays at the circle’s, $1.08,\Psi_{full}$, computed on the revised full-flow values. A force main is geometrically a circle but carries the Hazen–Williams section factor $\Psi = A R^{0.63}$ and its own peak ratio (§8.2).

Composite Shapes

Four shapes glue two primitives at a junction depth $Y_b$. Rectangular-triangular sets a triangular invert of height $Y_b$ under a rectangle of width $W$, with $s = W/(2Y_b)$, $r = \sqrt{1+s^2}$ and $A_b = WY_b/2$:

$$A(Y) = \begin{cases} sY^2, & Y \le Y_b\ A_b + (Y - Y_b),W, & Y > Y_b\end{cases} \qquad P(Y) = \begin{cases} 2rY, & Y \le Y_b\ 2rY_b + 2(Y - Y_b), & Y > Y_b.\end{cases}$$

Rectangular-round and modified basket-handle both pair a rectangle with a circular segment of radius $r_b$ subtending

$$\theta = 2\arcsin!\left(\frac{W}{2r_b}\right),\qquad A_b = \frac{r_b^2}{2},(\theta - \sin\theta),\qquad Y_b = r_b!\left(1 - \cos\frac{\theta}{2}\right),$$

differing only in which end the segment occupies — invert for rectangular-round, crown for modified basket-handle — so the latter’s widest point sits at $Y_{full} - Y_b$ rather than at the crown. A radius smaller than half the width is silently raised to it in both, since the arcsine would otherwise be undefined — but a radius large enough to make the segment taller than the conduit is rejected rather than adjusted, the only parameter combination in either shape that is fatal. Each composite carries its own crown cutoff: 0.98 for rectangular-triangular and rectangular-round, 0.96 for modified basket-handle.

Tabulated Shapes

The legacy masonry sewer profiles — basket-handle, catenary, egg, gothic, horseshoe, semi-circular, semi-elliptical — together with the ellipses and arches have no useful closed form, so they are carried as normalised property tables and rescaled by a fixed set of full-flow constants. Every one of these shapes is fully specified by the row below plus its tables: depth is the only user input, and all other dimensions follow.

The ellipses and arches sit outside that table because they take two dimensions rather than one. Each ships a catalogue of standard US sizes stored in inches — 23 ellipse codes, 102 arch codes — whose full-flow area and hydraulic radius are tabulated directly and selected by code. Given arbitrary user axes instead, they fall back to fixed proportionality constants:

The two ellipses share their constants but apply them to different axes — the horizontal to the minor axis it takes as its full depth, the vertical to the major axis it takes as its width — so a vertical ellipse’s area depends on width alone and its depth enters only through the normalised tables. The arch is the one shape whose full area is a product of both dimensions rather than a multiple of a squared one. All three set $\Psi_{max} = \Psi_{full}$, declining the section-factor peak that §6.1’s closing subsection describes for other closed shapes, even though their area cutoffs remain below unity (0.96 for the ellipses, 0.92 for arches). Table provision varies by shape, and sorts the library into four groups rather than the two its documentation suggests. All five circular tables are 51-entry. Egg, horseshoe and basket-handle carry 26-entry area, radius and width tables for the forward direction and 51-entry depth and section-factor tables for the inverse. Gothic, catenary, semi-elliptical and semi-circular carry only 21-entry width tables alongside those same 51-entry depth and section-factor tables, recovering area by inverse lookup on the depth table and hydraulic radius from the section factor as $R = (\Psi/A)^{3/2}$. The ellipses and arch are the exception in the other direction: 26-entry area, radius and width tables and no depth or section-factor table at all, so depth comes from an inverse lookup on the area table and the section factor is computed from area and radius rather than tabulated — consistent with their declining the section-factor peak. Lookups are linear except over the two lowest depth segments, where quadratic interpolation applies (with a linear fallback if it returns a non-positive value) — the near-empty geometry of every tabulated shape depends on it.

Custom Shapes

A custom shape integrates a user width-versus-depth curve into the same 51-point tables. The curve describes a unit-height section scaled by the conduit’s full depth: it is anchored at $(0,0)$, truncated above unit height, and extended at its last width if it stops short. Sweeping upward in 50 equal depth increments accumulates area by the trapezoid rule and perimeter as polyline arc length,

$$\Delta A = \frac{w_{i-1} + w_i}{2},\Delta y, \qquad \Delta P = 2\sqrt{\Delta y^2 + \left(\frac{w_i - w_{i-1}}{2}\right)^2},$$

after which the three tables are normalised by their full-depth values. The result is a closed section: the bottom width seeds the perimeter and the top width is added on the final step, so both count as wetted perimeter, and $A_{max} = 0.96,A_{full}$. Transects (§6.2) reach the same tables from surveyed geometry.

The Section-Factor Maximum

Closed shapes embed a subtlety every routing method must respect: the section factor peaks below full depth. The table above gives both halves of that fact per shape — how far the peak rises above the full-flow value ($\Psi_{max}/\Psi_{full}$) and the relative area at which it occurs ($\alpha_{max} = A_{max}/A_{full}$). For a circle the peak is 8% above full-pipe value and sits at 97.56% of full area, which means Manning flow at about 94% depth ($0.938,Y_{full}$) exceeds full-pipe flow. The analytic and composite shapes carry their own cutoffs — 0.97 for a closed rectangle, 0.98 for rectangular-triangular and rectangular-round, 0.96 for modified basket-handle and for custom shapes, 0.92 for arches — while every open shape has $\alpha_{max} = 1$ and no peak at all. SWMM stores $\Psi_{max}$ and the area at which it occurs, interpolates the non-monotone tail linearly between them,

$$\Psi(A) = \Psi_{max} + (\Psi_{full} - \Psi_{max}),\frac{\alpha - \alpha_{max}}{1 - \alpha_{max}}, \qquad \alpha > \alpha_{max},$$

and gives the derivative the constant value $(\Psi_{full} - \Psi_{max})/[(1 - \alpha_{max})A_{full}]$ over that band. The inverse $A(\Psi)$ then faces a two-branch ambiguity, resolved by bracketing: a section factor between $\Psi_{full}$ and $\Psi_{max}$ is bracketed on $[A_{full}, A_{max}]$ and anything else on $[0, A_{max}]$, then solved by Newton–bisection to $10^{-4}A_{full}$.

6.2 Transects

Transects represent natural channels by station-elevation pairs (up to 1,500 stations; an X1-line multiplier and offset can rescale the survey) with distinct left-overbank, main-channel, and right-overbank Manning coefficients — an omitted overbank $n$ defaults to the channel’s, and NC-line values persist as defaults into subsequent transects. That persistence is worth stating precisely, because the three roughness values and the survey buffer are section-scope state rather than per-transect fields, and the reader is driven by the NC line rather than by the transect record: an NC line first finishes the previous transect — building its geometry tables — and only then overwrites whichever of the three values it supplies as positive, a zero meaning “inherit”. Two consequences follow that the HEC-2 format description does not lead one to expect.

First, since only an NC line (or the section’s end) finishes a transect, every transect but the last must be followed by an NC line; the HEC-2 economy of repeating NC only when roughness changes leaves the intervening transects untabulated, their geometry tables all zero. This is not silently wrong — a conduit referencing such a transect fails validation on zero full area or zero roughness — but the error is reported against the link, as an invalid cross-section, with nothing pointing at the transect that was skipped.

Second, the meander adjustment below is applied to the shared channel roughness in place and is never restored. That the adjustment is meant to belong to one transect only is visible in the code itself: the pre-adjustment value is saved off and written to the transect’s own roughness field, so the record keeps the unadjusted number while the shared state keeps the adjusted one. A subsequent transect that inherits the channel $n$ therefore inherits the adjusted value, and the $\sqrt{L_{factor}}$ inflation compounds down the section — the saved-and-restored discipline applied to the record but not to the variable the next transect actually reads. Preprocessing appends vertical end walls at both ends (contributing wetted perimeter) and builds the 51-point tables by sweeping depth: each depth accumulates area, width, and wetted perimeter segment-by-segment, with composite roughness handled through conveyance summation — a new conveyance segment starts at each bank-roughness change and wherever the ground re-emerges above the water line (multi-thread sections sum correctly), each contributing $K_i = (1.486/n_i)A_i R_i^{2/3}$ — and the table’s hydraulic-radius entry back-computes an effective $R$ from total conveyance as $R = (n_C K / 1.49 A)^{3/2}$. The two directions are inverses of one equation and so must share its unit constant; the library names that constant once and the forward conveyance refers to it by name, with a comment pointing at its definition. The back-computation instead hard-codes a differently-rounded literal, $1.49$ against the named $1.486$. The intent that both carry the same Manning constant is therefore explicit in the code’s own structure, and the inversion does not round-trip: the tabulated hydraulic radius comes out $(1.486/1.49)^{3/2}$ of the consistent value, about 0.40% low, for every transect and street section in the model. A meander modifier substitutes the shorter overbank (valley) length for the meandering main-channel length as the conduit’s effective length, inflating main-channel roughness by the modifier’s square root to preserve friction loss. Street cross-sections (§8) compile to transects through the same machinery.

6.3 Storage Geometry and Characteristic Depths

Storage geometry integrates the surface-area description into volume: functional curves $A = c_0 + c_1 Y^{c_2}$ integrate analytically; the four analytical shapes of §1.3 are compiled at parse time into the common quadratic $A = a_0 + a_1 Y + a_2 Y^2$ and integrate as a cubic; tabular curves trapezoidal-integrate, with depth-from-volume solved analytically per segment (or by Newton–bisection for the functional, conical, and pyramidal forms). All three relations are evaluated in the user’s length and volume units, so the functional coefficients are unit-system-dependent (§16). Below a tabular curve’s first point, area is assumed to grow linearly from zero ($V = \tfrac{a_1}{2 y_1}y^2$); above its last, area extrapolates along the final segment’s slope — the regimes governing shallow and overfull storage. Curve lookups made through the extrapolating table reader — storage surface-area, custom-inlet capture, and exfiltration bottom-area curves — behave similarly: below the first point they extrapolate proportionally through the origin, above the last along the final slope — floored at zero, so a curve whose last segment narrows extrapolates flat rather than downward; outlet rating, pump, and weir-coefficient curves instead clamp to their end values. Critical depth — needed at free outfalls and free-fall discontinuities — uses exact formulas where they exist (rectangular, triangular, parabolic, power-law) and otherwise a root search on $A^3/W = Q^2/g$ seeded by a circular-pipe approximation — interval enumeration (25 fixed depth intervals with linear interpolation) when the section’s full area lies within a factor of two of the equivalent circle’s, Ridder’s method to 0.001 ft when it does not. Normal depth inverts the section factor: $Y_N = Y(A(\Psi = Q,n/1.486\sqrt{S_0}))$.

7. Pumps and Flow Regulators

7.1 Pumps

Pumps are links whose flow comes from a user curve, in five types plus one degenerate. Writing $\omega$ for the speed setting, $V_1$ and $y_1$ for the inlet node’s volume and depth, and $H_1, H_2$ for the end-node heads, the six rules are

$$Q = \omega\cdot\begin{cases} \hat{q}(V_1) & \text{Type 1, stepwise on wet-well volume,}\ \hat{q}(y_1) & \text{Type 2, stepwise on inlet depth,}\ q(H_2 - H_1) & \text{Type 3, the centrifugal characteristic,}\ q(y_1) & \text{Type 4, an in-line depth profile,}\ q!\left(\dfrac{H_2 - H_1}{\omega^2}\right) & \text{Type 5, a variable-speed Type 3,}\ Q_{in} & \text{ideal, with }\omega\text{ applied likewise,} \end{cases}$$

where $\hat{q}$ denotes stepwise lookup — the curve’s $y$-value at the first point whose $x$ exceeds the argument, so Types 1 and 2 are step functions while 3, 4 and 5 interpolate linearly. Type 5’s division of head by $\omega^2$ before lookup, with the resulting flow then multiplied by $\omega$, is the affinity-law scaling of the rated curve. The head argument is floored at zero, and an ideal pump must be its node’s only outlet. Types 3 and 5 report $\partial Q/\partial H$ as the negated curve slope divided by $\omega$ (flow falls as head rises), and Type 4 as a forward difference over 0.001 ft; the stepwise types report none. Startup and shutoff wet-well depths latch the pump on and off around the curve; at storage inlet nodes (and the virtual wet well a Type 1 pump receives at a non-storage node) flow is clamped so the node cannot be drawn below empty ($Q \le Q_{in} + V_N/\Delta t$), while Type 2–4 pumps at non-storage nodes fall back to $Q = Q_{in}$ when the projected end-of-step depth would go negative (Type 5 is omitted from this check); pumps contribute no surface area to their nodes, and reverse flow is never allowed. Energy is tallied as $0.7457,\Delta H,Q,\Delta t/3600/8.814$ kWh (no efficiency factor).

7.2 Orifices

Orifices (side or bottom, circular or rectangular, coefficient $C_d$, optional flap gate) discharge by Torricelli:

$$Q = C_d A_O \sqrt{2gH_e}$$

where $C_d$ is the discharge coefficient, $A_O$ the opening area, and $H_e$ the effective head, which switches between free-discharge (head above opening centre/invert) and differential (submerged tailwater) regimes. An unsubmerged inlet degrades smoothly to weir behaviour: below a threshold head the flow follows

$$Q = C_W L (H_1 - Z_O)^{1.5}$$

where $H_1$ is the upstream head and $Z_O$ the opening elevation. The weir coefficient is not a user input; it is derived by requiring the two regimes to agree at the changeover, and the derivation differs by orientation.

For a bottom orifice the changeover head is set at $h_c = (C_d/0.414)(A_O/P_O)$, with $P_O$ the opening’s perimeter and 0.414 the sharp-crested weir coefficient divided by $\sqrt{2g}$. Equating the two expressions there and substituting $h_c$ collapses the free constant:

$$C_d A_O\sqrt{2g,h_c} = C_W P_O h_c^{3/2} ;;\Longrightarrow;; C_W = 0.414\sqrt{2g} \approx 3.32,$$

so a bottom orifice below its changeover behaves as a sharp-crested weir of crest length equal to the opening perimeter, whatever $C_d$ the user supplied. For a side orifice the changeover is the opening height itself and the matching is done against the centre-line head $h/2$, giving $C_W L = C_d w\sqrt{g}$ for a rectangular opening of width $w$ — that is, the same $C_d$ carried across, rescaled by $\sqrt{g}$ rather than replaced. Submergence then applies a Villemonte factor $[1 - ((H_2 - Z_O)/(H_1 - Z_O))^{1.5}]^{0.385}$ on the heads above the crest. A partially-open setting $\omega$ (sluice-gate fraction, optionally slewing at a user open/close rate) re-computes the opening area from the §6 geometry. Flap gates charge the Armco head loss $\Delta H = (4U^2/g),e^{-1.15 U/\sqrt{H_e}}$, subtracted and re-solved. Under dynamic wave an orifice masquerades as an equivalent short pipe of length $\max!\big(200\ \text{ft},\ 2,\Delta t_{rout}\sqrt{g,y_{full}}\big)$, contributing surface area to its end nodes and analytic $\partial Q/\partial H$ ($0.5,Q/H_e$ submerged, $1.5,Q/(H_1 - Z_O)$ as a weir) to the surcharge update.

7.3 Weirs

Weirs come as transverse rectangular, V-notch, trapezoidal, and side-flow (Engels, reverting to the transverse form under reverse flow). The trapezoidal type is the sum of a rectangular centre and triangular ends evaluated with two independent discharge coefficients — the second applying to the end sections alone — so it is not a single weir equation with a composite crest. Each type admits exactly one cross-section shape — open rectangular for transverse, side-flow, and roadway; triangular for V-notch; trapezoidal for trapezoidal — any other being a fatal error. The head-discharge relations are:

$$Q = C_W L_e H_e^{3/2} \qquad \text{(transverse rectangular)}$$

$$Q = C_W \tan(\theta/2)H_e^{5/2} \qquad \text{(V-notch)}$$

$$Q = C_W L_e^{0.83} H_e^{1.67} \qquad \text{(side-flow, Engels)}$$

where $C_W$ is the weir discharge coefficient, $L_e$ the effective crest length, $H_e$ the effective head, and $\theta$ the notch angle. Effective crest length subtracts end contractions ($L_e = L - 0.1,n_c H_e$, floored at zero, so a high enough head on a heavily contracted weir stops flow entirely rather than reversing it); a partially-raised crest ($\omega < 1$) turns a V-notch into a trapezoid; submergence applies Villemonte with the type’s own head exponent, save that a trapezoidal weir’s triangular end sections always take the V-notch exponent 2.5. Like the orifice, a weir acquires an equivalent length $\max!\big(200\ \text{ft},\ 2,\Delta t_{rout}\sqrt{g,y_{full}}\big)$ for its nodal surface area, but unlike the orifice its discharge equations are evaluated in the user’s length and flow units (§16), making $C_W$ unit-system-dependent — the roadway weir being the exception, computing in internal feet and instead rescaling a user-supplied $C_D$ by $1/0.552$ under SI. Weirs default to surchargeable (roadway weirs excepted): above the opening they switch to an equivalent-orifice form $Q = C_O\sqrt{H_e}$, and the matching is a two-part convention worth stating exactly, since either half alone is misleading. The coefficient is fixed by evaluating the weir equation at a head equal to the full opening height and then dividing that flow by the square root of half that height, $C_O = Q_{weir}(h)/\sqrt{h/2}$; at run time the orifice form is driven by the head measured to the opening’s mid-height — differential once the tailwater rises above it. It is the centre-line head convention, borrowed from the side orifice, that makes the two regimes agree at the changeover. A weir with surcharging disabled instead simply caps the head at its opening height and continues weir-equation flow. Under steady/kinematic routing, all regulators discharge freely — head is always computed against the upstream node’s invert as tailwater, never submerged. $\partial Q/\partial H$ is the analytic exponent-scaled ratio per type.

7.4 Outlets

Outlets are the catch-all: flow from a power function $Q = aH_e^b$ or tabulated rating curve of either upstream depth or head difference, scaled by the setting, with flap-gate reversal blocking — the vehicle for vortex valves and other bespoke devices. Both forms are evaluated in the user’s length and flow units (§16), so $a$ is unit-system-dependent whenever $b \ne 1$.

8. Advanced Hydraulics

8.1 Conduit Evaporation and Seepage

Conduit evaporation and seepage are uniformly-distributed lateral losses: $q_E = e_t W(\bar Y)$ (open channels) and $q_S = s f_c W(\bar Y)$ (seepage, with monthly adjustment $f_c$ and the width capped at the depth of maximum width, since seepage is vertical), together bounded by the conduit volume per step (dynamic wave) or the flow magnitude (steady/kinematic). The momentum equation itself gains Strelkoff’s lateral-outflow term $-\bar U q_L/2$, which after substituting continuity becomes a $+2.5,\bar U q_L,\Delta t / L$ term in the dynamic-wave flow-update numerator (with $L$ the conduit’s true length, not the lengthened one); the lost volume debits the appropriate node; kinematic wave adds $q_L L/\phi$ into its $C_2$ constant. Storage units evaporate at the potential rate times a user-supplied realisation fraction $f_E$ (1 normally, 0 for roofed units — and 0 by default, so an unconfigured storage unit does not evaporate at all), applied to the start-of-step surface area, and seep through bottom and sloped-side areas separately. The seepage law depends on how the unit was parameterised: supplying only a saturated conductivity gives a constant rate $K_s f_c$ on both surfaces, while a full suction-head/conductivity/deficit triple invokes modified Green–Ampt (§3.3). The driving head differs between the two surfaces — the bottom sees the ponded depth, the banks half the depth above the elevation at which the storage curve starts widening (and a shifted expression once it stops), with bank area taken as the surface area capped at its widest tabulated value less the bottom area. One shape is unhandled. The exfiltration initialiser branches on storage geometry and sets bottom area and bank limits for the tabular, functional, cylindrical, conical and pyramidal cases; the elliptical paraboloid — added to the storage geometry of §1.3 in 5.2.0, after exfiltration — has no case and the switch has no fallback. Since the object is allocated without zeroing, the evident intent of covering every shape is not met and a paraboloid storage unit with seepage reads uninitialised bottom and bank geometry.

8.2 Minor Losses and Force Mains

Minor Losses

Minor losses (entrance, exit, average, with velocities evaluated at the respective locations) enter the flow-update denominator as $\frac{\Delta t}{2L}\sum K_{m,i}|U_i|$ — dynamic wave only.

Force Mains

Force mains (circular, dynamic wave) swap the Manning friction term for Hazen–Williams ($\Delta Q_{friction} = 0.6g|\bar U|^{0.852}\Delta t / C_{HW}^{1.852}R_{full}^{1.1667}$ — note the 7/6 hydraulic-radius exponent) or Darcy–Weisbach ($f|\bar U|\Delta t/8R_{full}$, with Swamee–Jain $f$, laminar $64/Re$ below 2000 with $Re$ floored at 10, a linear blend from $f = 0.032$ between 2000 and 4000, and the fully-rough form above $Re = 10^{10}$) only while pressurised; partly-full flow uses an equivalent Manning $n$ — slope-dependent for Hazen–Williams, $n = 1.067,C^{-1}(D/S_0)^{0.04}$, and $n = \sqrt{f/185};D^{1/6}$ for Darcy–Weisbach with $f$ evaluated at $Re = 10^{12}$. Force mains also carry their own conduit-lengthening compensation (dividing friction by the length factor rather than rescaling $n$), and the force-main cross-section’s section factor uses the Hazen–Williams exponent $A R^{0.63}$ instead of Manning’s $R^{2/3}$ — altering its normal-flow limit.

8.3 Culverts and Roadway Weirs

Culverts

Culverts designated by an FHWA HDS-5 code get an inlet-control capacity check layered on the ordinary dynamic-wave (outlet-control) solution: unsubmerged flow from either the form-1 critical-energy equation (Ridder’s method on critical depth) or the form-2 power law, submerged flow from the quadratic HDS-5 relation, a linear transition between (SWMM places the unsubmerged limit at $H_1 < Z_1 + 0.95,Y_{full}$), and the smaller of the two flows governs. Each relation carries a slope-correction term, and its sign convention needs stating carefully because the source and the manual write it oppositely. HDS-5 puts the term on the discharge side of $HW/Y_{full}$, as $-0.5,S_O$ for ordinary inlets and $+0.7,S_O$ for mitered ones. SWMM instead carries a quantity it adds to $HW/Y_{full}$ (and subtracts in the submergence bound below), so the equivalent literals are $+0.5,S_O$ and $-0.7,S_O$ — and the ordinary case is indeed coded as $+0.5,S_O$. The mitered case, however, is coded as $-7.0,S_O$: correct in sign for this convention, but ten times the published magnitude, so a mitered culvert on any appreciable slope carries a slope correction an order of magnitude too large. The constants come from a compiled table of 57 inlet configurations storing the HDS-5 Table H-2 values (form, $K$, $M$, $c$, $Y$) — 58 rows, the first a zero placeholder since culvert codes are one-based — of which codes 5, 37, and 46 are the mitered ones; submergence begins at $y = Y_{full}(16c + Y - S_{cf})$ for the correction $S_{cf}$ above, fixing the transition band’s upper bound, while the unsubmerged bound is what the source itself calls an arbitrary limit.

Roadway Weirs

Roadway weirs apply the FHWA head-dependent coefficient only when both a road width and surface type are given (otherwise the user’s constant $C_D$); the “charts” are small digitised piecewise-linear tables — low-head coefficients looked up against absolute head in feet below $h/W_{road} = 0.15$ and against the ratio above, with submergence factors bottoming at 0.40 (paved) / 0.24 (gravel); they are typically paired in parallel with a culvert to model embankment overtopping.

8.4 Streets and Inlets

Streets and inlets (new in 5.2) implement HEC-22 dual drainage. A street section — crown width, curb height, cross slope, optional depressed gutter and backing, and a one- or two-sided flag (default two; approach flow halves and capture doubles per side; the cross and backing slopes are entered as percentages) — compiles into a §6 transect, so street conduits route like any channel. Inlet designs comprise: grates (seven standard types with open-area ratios and splash-over velocity fits, plus a generic type with user-supplied values), curb openings with three throat geometries, slotted drains, drop grates/curbs, custom-curve inlets, and the implicit combination inlet formed when one design defines both a grate and a curb. Placement rules are shape-checked: street inlets (grate/curb/combo/slotted) belong only in street cross-sections, drop inlets only in rectangular or trapezoidal open channels, custom inlets anywhere with a diversion or rating curve; an invalid placement is removed with a warning, and a conduit holds at most one inlet usage (a second definition overwrites). Each usage line adds modifiers: a replicate count (on-grade replicates evaluate sequentially, each seeing the previous one’s bypass; on-sag they multiply), a clogging percentage scaling both capture and the open area used for backflow apportioning, a per-inlet flow cap, and a local gutter depression added to the street’s continuous one. AUTOMATIC placement resolves to on-grade when the bypass node has an outgoing link, else on-sag; drop-curb inlets always compute in depth-driven (on-sag) mode capped by approach flow, their opening length quadrupled to stand for the full perimeter of the drop; and custom inlets ignore placement entirely (diversion curve = flow-driven, rating curve = depth-driven).

On-Grade Capture

On-grade capture starts from the gutter-spread relation

$$Q = \frac{0.56}{n}\sqrt{S_L},S_x^{1.67},T^{2.67}$$

where $n$ is the gutter’s Manning roughness, $S_L$ the longitudinal slope, $S_x$ the cross slope, and $T$ the spread. Inverted, this gives the spread from a known flow as $T = (Q/f)^{0.375}$ with $f = (0.56/n)\sqrt{S_L},S_x^{1.67}$, and that inverse is what the engine actually evaluates.

A depressed gutter of depression $a$ over width $W$ breaks the single-section formula, because the flow then spans two cross slopes. HEC-22 handles it through the frontal-flow ratio $E_o$ — the fraction of total gutter flow carried within the depressed width — given by

$$E_o = \left[1 + \cfrac{S_r}{\left(1 + \cfrac{S_r}{T_s/W}\right)^{2.67} - 1}\right]^{-1},\qquad S_r = \frac{S_x + a/W}{S_x},$$

where $T_s = T - W$ is the spread beyond the depressed width and $S_r$ the ratio of depressed to normal cross slope. Because $E_o$ depends on the spread and the spread depends on $E_o$, the pair is solved by fixed-point iteration: an initial $T_s$ is refined by alternately computing $E_o$, taking the side flow $Q_s = (1 - E_o)Q$, and re-inverting the spread relation on $Q_s$ — at most ten passes to a 0.01 ft tolerance — after which $T = T_s + W$. A first check short-circuits this whenever the flow fits entirely within the depressed width. The resulting spread is capped at the distance from curb to crown, and for an undepressed gutter $E_o$ collapses to the closed form $1 - (1 - W/T)^{2.67}$. Grates apply a frontal efficiency (above splash-over)

$$R_f = 1 - 0.09(V - V_o)$$

and a side efficiency

$$R_s = [1 + 0.15V^{1.8}/S_xL^{2.3}]^{-1}$$

where $V$ is the gutter velocity, $V_o$ the splash-over velocity, and $L$ the grate length. Curb openings use the equivalent slope $S_e = S_x + (a/W)E_o$, the full-capture length

$$L_T = 0.6,Q^{0.42}S_L^{0.3}(nS_e)^{-0.6}$$

and the efficiency

$$E = 1 - (1 - L/L_T)^{1.8}$$

where $L$ here is the curb-opening length. On-grade slotted drains are treated as curb openings of equal length; combination inlets capture through the curb “sweeper” (curb length beyond the grate) first, then the grate on the remainder at recomputed spread.

On-Sag Capture

On-sag capture is weir flow at shallow depth, orifice flow at depth — curb openings linearly interpolating across the transition band, grates and slotted drains switching outright at their equal-flow depths (no discontinuity either way): grate weir $3.0,P,d^{1.5}$ ($P = L_g + 2W_g$; full perimeter for drop grates) switching at $d = 1.79,A_o/P$ to orifice $0.67,A_o\sqrt{2gd}$ — with the regime test taken on the raw water level but both flows evaluated at the depression-corrected mean depth over the grate; curb weir $3.0,L,d^{1.5}$ (or $2.3(L + 1.8W)d^{1.5}$ with crest at $h + a$ when depressed and no longer than 12 ft) to orifice $0.67,hL\sqrt{2g,d_{eff}}$ above $d = 1.4h$, with throat-angle head corrections; slotted weir $2.48,L,d^{1.5}$ to orifice $0.8,Lw\sqrt{2gd}$ at $d = 2.587w$; a combination adds curb-orifice flow over the grate length once the grate is in orifice mode. (The inlet equations use HEC-22’s $g = 32.16$ ft/s², not the engine’s 32.2.)

Capture Transfer and Statistics

Captured flow transfers from the street conduit’s downstream (bypass) node to the sewer capture node each routing step, carrying pollutant mass at the bypass node’s previous-step concentration; sewer surcharge returns as backflow at the capture node’s concentration, apportioned by open-area ratio among standard inlets (by count among custom ones) sharing the node — flooding that stays inside the model, with continuity accounting corrected accordingly. Under steady/kinematic routing, on-sag capture is additionally limited to the inlet’s share of bypass-node inflow plus stored volume per step. One structural caveat: the gutter-spread factor is cached at validation from the conduit’s bed slope under a normal-flow assumption, so on-grade capture is insensitive to dynamic-wave backwater — though the reported maximum street spread is depth-based and does reflect it. Per-inlet statistics (flow/capture/backflow period counts after report start, capture efficiency at peak approach flow, average efficiency, bypass/backflow frequencies, peak flows) feed the street-flow summary, which lists every street conduit with or without an inlet.

9. Water Quality

9.1 Pollutants and Sources

A pollutant is any constituent expressible as an additive concentration (mass or organism counts per volume) — which deliberately excludes pH, conductivity, turbidity, and colour. Each carries optional background concentrations in rainfall, groundwater, RDII, and dry-weather flow; a first-order decay coefficient (days⁻¹) active in the conveyance system — which may be negative to model growth, though growth is inert on the steady-flow path; a snow-only buildup flag (de-icing chemicals); and an optional co-pollutant relation $C_{total,i} = C_i + f_{ij}C_j$ — the HSPF-style potency factor, applying to buildup/washoff loads only, with $f_{ij}$ free to exceed 1 since it bridges the two constituents’ units. Mass enters the conveyance system at a node through exactly eight paths, assembled in this order at the start of each routing step:

Two further paths move mass within the system rather than into it: an inflowing link delivers its previous-step concentration to its downstream node (§9.4), and street inlets transfer captured mass from a bypass node to a sewer capture node — with sewer surcharge returning it as backflow — at the donating node’s previous-step concentration (§8.4).

9.2 Buildup and Street Sweeping

Land uses partition each subcatchment purely for quality: each (pollutant, land use) pair owns one buildup and one washoff function. Buildup $b$ (mass per unit area or per unit curb length) grows with dry time $t$ by one of three forms:

$$b_{pow}(t) = \min!\left(B_{max},\ K_B t^{N_B}\right),\qquad b_{exp}(t) = B_{max}!\left(1 - e^{-K_B t}\right),\qquad b_{sat}(t) = \frac{B_{max},t}{K_B + t},$$

where $B_{max}$ is the maximum attainable buildup, $K_B$ a rate constant (a half-saturation time in the saturation form), and $N_B$ a power exponent. Two facts about how these reach the three forms matter to a reader of the file. The three coefficients occupy fixed columns, but the saturation form takes its half-saturation time from the third column while the power and exponential forms take their rate constant from the second — so a saturation buildup ignores its second coefficient entirely. And $N_B$ is validated to $[0.01, 10]$ only when it is positive; since all three coefficients are separately rejected for being negative, the admissible set is ${0} \cup [0.01, 10]$.

The true state, however, is the accumulated mass, not the elapsed time. Each dry step therefore inverts the chosen form to recover the equivalent time for the mass on hand,

$$t_{pow} = \left(\frac{b}{K_B}\right)^{1/N_B},\qquad t_{exp} = -\frac{1}{K_B}\ln!\left(1 - \frac{b}{B_{max}}\right),\qquad t_{sat} = \frac{b,K_B}{B_{max} - b},$$

advances it by the step, and re-evaluates the forward form. Washoff and sweeping consequently rewind the clock rather than resetting it, and buildup resumes along the same curve from wherever the remaining mass places it. (The pinned source adds a fourth, external option — a scaled user time-series loading rate capped at a maximum — that bypasses the inversion mechanism entirely.)

Each form also precomputes a time-to-maximum beyond which buildup is pinned exactly at $B_{max}$:

$$T_{pow} = \left(\frac{B_{max}}{K_B}\right)^{1/N_B},\qquad T_{exp} = -\frac{\ln 0.001}{K_B},\qquad T_{sat} = 1000,K_B,$$

with the power form swapped for a flat 3650 days when $\log_{10}(B_{max})/N_B > 3.5$ — a blow-up guard rather than a cap, since the test ignores $K_B$ and for small enough $K_B$ the computed value exceeds 3650 anyway. The admission of a zero exponent above has a consequence here: a power form with $K_B = 0$ or $N_B = 0$ takes a time-to-maximum of zero, and since buildup is pinned at $B_{max}$ once the elapsed time reaches that figure, such a land use jumps to maximum buildup after its first dry step. A POWER line written with zeroed coefficients to mean “no buildup” therefore produces the opposite. Initial buildup comes from a user areal loading or, absent one, from evaluating the buildup function over the antecedent dry days; each land use also carries an initial days-since-last-swept that offsets its first sweeping. Buildup pauses during wet steps (runoff > 0.001 in/hr), and snow-only pollutants accumulate only while snow depth is at least 0.001 in. Street sweeping runs on a per-land-use interval within a seasonal window; each pass removes the fraction (availability × efficiency) of current buildup, and is suppressed when rainfall exceeds 0.001 in/hr, when more than 0.05 in of snow lies on the plowable impervious area, or when the interval is zero.

9.3 Washoff

Three per-(pollutant, land-use) washoff models, all cut off below 0.001 in/hr of runoff:

$$w_{exp} = K_W,q^{N_W} m_B \ \ \text{(mass/hr)},\qquad w_{rat} = K_W,(f Q_{sub})^{N_W} \ \ \text{(mass/s)},\qquad w_{emc} = C_{emc},f Q_{sub},$$

where $q$ is the runoff intensity over the whole subcatchment in in/hr (or mm/hr), $m_B$ the remaining buildup mass on that land use, $Q_{sub}$ the subcatchment’s runoff rate, $f$ the land use’s area fraction, and $C_{emc}$ a constant event-mean concentration. The engine reaches all three through a concentration — it forms $c = w/(f Q_{sub})$ and re-multiplies by the exported volume — which is why the three carry different natural units.

• Exponential: first-order in remaining buildup $m_B$, driven by runoff intensity over the whole subcatchment; each step depletes buildup by $\min(w,\Delta t,\ m_B)$ before BMP removal is applied. Source-limited by construction, producing the classic first-flush hysteresis. The exponent $N_W$ generalises the original linear $k = K_W q$, whose $q$-cancellation forced concentration to decrease monotonically; the classic $K_W = 4.6$ in⁻¹ is Burdoin’s separate calibration (“half an inch of runoff in an hour removes 90% of the load”).
• Rating curve: evaluated on the land-use share of flow rather than prorated after the fact, so the rate is exactly $K_W(f,Q_{sub})^{N_W}$ and not the linear proration $f,K_W Q_{sub}^{N_W}$ — the two differ whenever $N_W \ne 1$. No inherent source limit unless buildup is also modelled as a cap, in which case exhaustion drops the load abruptly to zero. When no buildup function is paired with the washoff, there is no mass to draw down and the surface-loading balance of §12 would not close; SWMM keeps it closed by booking each washoff load as an equal, simultaneous buildup input to the ledger — a phantom source that exists only for the accounting, and which is why an EMC-only model still balances.
• EMC: the rating curve with $N_W = 1$, its coefficient absorbing the 28.3 L/ft³ conversion.

Rain and run-on loads are not simply added: they mix through the ponded water atop the subcatchment (a completely-mixed store consistent with the nonlinear reservoir), which introduces one extra state per pollutant per subcatchment — the ponded mass $M_p$. Over a step the store takes wet deposition and run-on, then loses mass to infiltration and to runoff in proportion to their volumes:

$$c_p = \frac{M_p^{,t} + C_{rain}V_{rain} + M_{runon}}{V_{in}},\qquad M_{infil} = c_p V_{infil},\qquad M_{out} = c_p V_{out},$$

$$M_p^{,t+\Delta t} = c_p,d_{pond},A_{nonLID},$$

where $V_{in}$ is the step’s total inflow volume to the store — start-of-step ponded volume plus run-on plus net precipitation — and $d_{pond}$ the new ponded depth. The two loss masses are each additionally clamped to the mass still on hand, applied in order so that the outflow share is drawn from a store already reduced by infiltration; the clamps are inert whenever the loss volumes sum to no more than $V_{in}$, which is the normal case. Writing the residual mass this way rather than as $M_p^{,t} + \text{gains} - \text{losses}$ has a consequence worth stating: the volume lost to evaporation is absent from $d_{pond}$, so its share of mass silently leaves the ledger — nothing concentrates, and the evaporated mass is neither retained nor booked as a loss; and a step with no inflow at all writes any residual ponded mass off to final storage — it does not persist for resuspension in the next storm. Per-land-use BMP removal fractions discount the washoff stream, and their area-weighted average discounts the ponded stream; the outflow concentration is total load over outflow (the pre-re-routing runoff rate drives the washoff rate, but both load streams are exported on the post-re-routing outflow volume). Loads to another subcatchment become its run-on at the next step.

9.4 Transport and Treatment

Conveyance Transport

Conveyance transport treats every conduit and storage node as a completely-mixed reactor (the WASP/QUASAR box-model lineage) rather than solving advection–dispersion. At each routing step: node inflow loads are accumulated (subcatchments, DWF, external, groundwater, RDII, plus each inflowing link at its previous concentration); non-storage nodes holding negligible volume take the flow-weighted mixture (a junction actually holding water — surcharged or ponded — updates as a mixed reactor instead); storage nodes and conduits update by the deliberately robust mixing formula

$$c(t{+}\Delta t) = \frac{c(t),V(t),e^{-K_1\Delta t} + C_{in}Q_{in}\Delta t}{V(t) + Q_{in}\Delta t}$$

chosen over the analytical CSTR solution (exact only under constant-inflow, averaged-volume assumptions) because it stays stable as volumes vanish and never overshoots a step input. Two code-level deviations from this manual form: the pinned source evaluates the decay factor as the linear truncation $(1 - K_1\Delta t)$ floored at zero — the exponential survives only on the steady-flow path — and clamps the mixed result to at most the larger of the reactor and inflow concentrations. Under dynamic wave (which yields one flow per conduit) the mixing inflow is volume-adjusted, $Q_{in} \leftarrow \max(0,\ Q_{in} + (V_2 + V_{losses} - V_1)/\Delta t)$. The dry thresholds are concrete: below 1 litre of volume or 1 mm of depth an element’s remaining mass is flushed to final storage and its concentration zeroed — unconditionally for conduits, but only in the absence of inflow for the mixed-reactor nodes (initial concentrations seed only elements wet at start; a wet no-inflow junction keeps its previous concentration). Volume-less links (pumps, regulators, dummy conduits) pass their upstream node concentration through; evaporation concentrates by $1 + V_{evap}/V$; steady-flow routing replaces conduit contents with the upstream node concentration decayed by $e^{-K_1\Delta t}$ and scaled by the evaporation factor.

Treatment

Treatment attaches a user expression to any (node, pollutant): either c = … (resulting concentration) or r = … (fractional removal applied to the inflow concentration), written over:

• pollutant symbols — a pollutant’s bare name (TSS) for its concentration and R_<name> for its fractional removal;
• hydraulic variables:
  • FLOW, in user flow units;
  • DEPTH and AREA, as old/new-step averages in user length units;
  • DT, in seconds;
  • for storage nodes, HRT — the residence time in hours, updated as $\theta \leftarrow (\theta + \Delta t),V/(V + Q_{in}\Delta t)$ — zero elsewhere;
• and the expression language’s 19 functions (sin cos tan cot asin acos atan acot sinh cosh tanh coth abs sgn sqrt log log10 exp step, case-insensitive, with + - * / ^ and scientific literals).

Domain violations do not error: square roots and logarithms of non-positive arguments, powers of non-positive bases, and NaN results all silently evaluate to zero. A subtle semantic: a referenced pollutant symbol denotes the combined-influent concentration when that pollutant’s equation at the node is removal-type, and the node’s pre-treatment concentration otherwise — equivalent only at nodes holding no volume. A pollutant with no equation at the node resolves to the influent concentration as well, and not by an explicit rule: the treatment records are zero-initialised and the removal type is the zero-valued one, so “no equation” is indistinguishable from “removal-type” to the variable lookup. This small expression language expresses constant EMCs, co-removal, concentration-switched removal, $n$-th-order kinetics, the k-C* wetland model, and quiescent gravity settling. Guardrails: treated concentration bounded by [0, untreated]; removals ≤ 1; removal-form yields zero without inflow; a treatment expression at a node overrides the pollutant’s global decay there; and co-pollutants receive no automatic co-treatment.

10. LID Controls

10.1 The Generic Layered Unit

LID units are depth-explicit layered moisture-accounting models embedded in subcatchments — a deliberate middle path between curve-number credits (no dynamics) and Richards-equation soil physics (too costly for hundreds of units). The generic unit (a bio-retention cell) stacks a surface layer (ponding depth $d_1$, void fraction $\phi_1$), a soil layer (moisture $\theta_2$ across thickness $D_2$), and a storage layer (depth $d_3$, void fraction $\phi_3$) with optional underdrain; each layer’s state advances by a flux balance of the form

$$\phi_1\frac{\partial d_1}{\partial t} = i + q_0 - e_1 - f_1 - q_1,\qquad D_2\frac{\partial \theta_2}{\partial t} = f_1 - e_2 - f_2,\qquad \phi_3\frac{\partial d_3}{\partial t} = f_2 - e_3 - f_3 - q_3 .$$

10.2 Constitutive Fluxes

The constitutive fluxes echo the engine’s own hydrology, with LID-specific parameters: surface-to-soil infiltration is modified Green–Ampt in the amended-media parameters (except permeable pavement, whose surface intake is inflow-plus-ponding capped by the clog-reduced pavement permeability — no Green–Ampt; and vegetative swales, which use the parent subcatchment’s live native infiltration — but only when that subcatchment’s model is Green–Ampt or modified Green–Ampt, a Horton or curve-number parent leaving the swale’s soil conductivity at zero); soil percolation follows the exponential form $K_{2S}e^{-k_{slope}(\phi_2 - \theta_2)}$, zero below field capacity, where $k_{slope}$ is the LID soil layer’s own conductivity-slope input, not the aquifer’s $HCO$; exfiltration to native soil is its saturated conductivity, further capped by the groundwater module’s available upper-zone storage when an aquifer is modelled; the “native” infiltration rate the swale and exfiltration limits read is not the pervious sub-area’s rate but the subcatchment’s infiltrated volume for the step spread over its whole non-LID area — impervious fraction included, so the rate is diluted — falling back to a direct infiltration evaluation only when no non-LID pervious area exists; ET cascades top-down through the layers from the same potential-ET series — but all sub-surface ET is suppressed while surface infiltration is active, the storage layer’s evaporation is zeroed while the soil or pavement layer above is saturated (in trenches, while the surface is ponded), and rain barrels evaporate nothing at all, covered or not. The underdrain is a power law $q_3 = C_{3D}h_3^{\eta_{3D}}$ ($\eta = 0.5$ recovers the orifice equation) with head cases spanning storage, saturated-soil, and ponded regimes — extended in 5.1.013 with open/close head thresholds (hysteretic on prior drain flow) and an optional multiplier-vs-head curve; notably the drain equation is evaluated in user units (head in in/mm, flow in in/hr or mm/hr), an exception to the engine’s internal ft–s convention that makes drain coefficients unit-system-dependent. Surface excess over the berm overflows — by Manning routing $\alpha(d - D_1)^{5/3}W_1/A_1$, itself capped at $(d - D_1)/\Delta t$ so a long step cannot drain more than the excess actually present, when roughness, slope, and width are all non-zero, and instantaneously otherwise.

Per-Type Flux Balances

The generic triple of §10.1 is the bio-retention cell. The rain garden is the same cell without a storage layer, but dropping that layer also changes how the fluxes are limited: instead of the cascade below it applies an unconditional equal-flux rule between percolation and exfiltration — $f_2 = f_3 = \min(f_2, f_3)$ whatever the moisture state — followed only by the surface-infiltration cap, and it zeroes storage evaporation. The other six types each carry their own balance:

Green roof. The storage layer becomes a drainage mat, sealed against the roof deck so exfiltration vanishes and the drain is Manning flow along the mat rather than a power-law orifice:

$$f_3 = 0, \qquad q_3 = \alpha_{mat},d_3^{5/3},\frac{W_1}{A_1},\phi_{mat},$$

with $\alpha_{mat} = 1.49\sqrt{S_1}/n_{mat}$ taken on the surface layer’s slope. A mat given no roughness has no $\alpha_{mat}$, and the drain then falls back not to zero but to passing the whole soil percolation rate straight through — the mat becomes transparent rather than impermeable, mirroring the instantaneous case for surface overflow.

Infiltration trench. No soil layer, so the surface discharges straight into storage and $\partial\theta_2/\partial t \equiv 0$:

$$\phi_1\frac{\partial d_1}{\partial t} = i + q_0 - e_1 - f_{13} - q_1, \qquad \phi_3\frac{\partial d_3}{\partial t} = f_{13} - e_3 - f_3 - q_3,$$

where $f_{13}$ is a single surface-to-storage flux limited by both ends.

Permeable pavement. A fourth layer sits between surface and soil, and its intake is not Green–Ampt: the surface delivers all it has, capped by the clog-reduced pavement permeability prorated by the pervious paver fraction,

$$f_1 = \min!\left(i + q_0 + \frac{d_1\phi_1}{\Delta t},\ \ K_{pave},(1 - f_{imp})\right),\qquad K_{pave} = K_{pave,0}!\left(1 - \min!\left(\frac{V_{treated}}{C_{clog}},,1\right)\right).$$

A soil layer is optional; without one the pavement percolates directly to storage.

Rain barrel. Pure storage with a sealed bottom, no evaporation and no infiltration, so the balance reduces to a two-term chain:

$$\frac{\partial d_1}{\partial t} = i + q_0 - f_{13}, \qquad \frac{\partial d_3}{\partial t} = f_{13} - q_3,$$

with $f_{13}$ limited by the barrel’s remaining freeboard plus whatever the drain is removing over the same step, so a barrel draining while it fills accepts more than its headroom alone would allow.

Rooftop disconnection. A lone surface layer whose drain is a gutter capacity that pre-empts the overflow rather than adding to it:

$$q_3 = \min!\left(C_{3D},\ q_1^{}\right), \qquad q_1 = q_1^{} - q_3,$$

where $q_1^{*}$ is the Manning surface outflow computed before the split — so a gutter wide enough to take everything leaves no overflow at all.

Vegetative swale. The one type whose geometry varies with depth. A trapezoidal channel of top width $W_1$, side slope $z$ (run per rise) and berm height $D_1$ has bottom width $b = W_1 - 2zD_1$, and at depth $d$

$$w(d) = b + 2zd,\qquad A_{flow}(d) = \phi_1 d,(b + zd),\qquad L = \frac{A_{unit}}{W_1},$$

so the wetted surface area is $L,w(d)$ and the stored volume $L,A_{flow}(d)$. Evaporation and infiltration act on the surface area rather than the unit footprint, and the state equation is written on volume,

$$\frac{\partial V}{\partial t} = Q_{in} - E,L,w(d) - f,L,w(d) - Q_{out},$$

then converted back to a depth rate by dividing by $L,w(d)$. Both widths are floored at 0.5 ft; if the declared top width and berm height imply a bottom narrower than that, the side slope is recomputed as $z = (W_1 - 0.5)/(2D_1)$ to keep the section consistent.

The Limiter Cascade

Unconstrained, the constitutive rates above would move more water than the layers hold. Every unit type therefore applies the same ordered set of min-limits, and the order is part of the model rather than an implementation detail — each flux is clipped to what its source can supply before the layer beneath is asked what it can accept, so the constraint propagates downward and then rebounds:

1. Soil percolation is capped at the drainable water above field capacity net of soil ET, $\big[(\theta_2 - \theta_{FC})D_2/\Delta t\big] - e_2$, and floored at zero.
2. Storage exfiltration is capped at what percolation delivers plus the storage already present, net of storage ET.
3. Underdrain flow is capped at the volume standing above the drain offset, plus percolation once storage is full, less exfiltration and storage ET.
4. Percolation is then re-capped at the storage layer’s remaining freeboard plus everything leaving storage.
5. Surface infiltration is finally capped at the soil layer’s remaining void volume plus everything leaving the soil.

When soil and storage are both saturated an equal-flux rule replaces the cascade: the whole chain collapses to the smallest of percolation and total storage outflow, and if percolation is the binding constraint the underdrain absorbs the remainder after exfiltration takes its share — $q_3 = f_2 - f_3$ when $f_2 > f_3$, else $f_3 = f_2$ and $q_3 = 0$. Surface infiltration is then capped at that same limiting rate, so the constraint reaches the top of the unit in one step rather than through the cascade. This is what stops a saturated bio-retention cell from draining faster than its media can pass water.

10.3 Unit-Type Variants

The other unit types are configurations of this template: rain gardens drop the storage layer; green roofs replace it with a drainage mat drained by Manning flow along the roof; infiltration trenches drop the soil layer; permeable pavement inserts a pavement layer (with block-paver area fraction, permeability, and optional sand filter); rain barrels are pure storage (void fraction forced to 1, sealed bottom) with a delayable drain valve and an optional cover flag — uncovered barrels receive direct rainfall, covered ones none; rooftop disconnection is a lone surface layer with a gutter-capacity-limited drain; vegetative swales are a lone surface layer with trapezoidal, depth-varying geometry and Manning outflow. Any gravel storage layer — and the pavement layer — may clog: conductivity declines linearly with cumulative void-volumes of inflow treated, controlled by a single clogging factor per layer — the user’s factor being a multiple of the layer’s own void volume, and the accumulator the storage layer’s total inflow but the pavement layer’s total volume treated. Pavement permeability may optionally regenerate on a fixed-day cycle (5.1.013+): regeneration does not restore the layer to new but discounts its accumulated treated volume by a user degree, so a degree of 1 resets the clock and anything less only winds it partway back.

10.4 Integration

Numerically, vegetative swales integrate their layer-state vector by the iterated trapezoidal method ($\Omega = 0.5$, 1 mm tolerance, at most 20 passes) — the LID module’s one implicit solve; every other unit type advances by a single explicit Euler step, which testing showed sufficient.

10.5 Deployment and Routing of Outflows

Deployment is per-unit-area: each unit captures a specified percentage of the subcatchment’s non-LID impervious-area runoff (percentages across units validated to sum ≤ 100%; a “capture ratio” of areas is only a sizing heuristic) and, since 5.1.013, a percentage of pervious-area runoff as well — both reduced by any internal sub-area re-routing first. A second, independent check requires the units’ combined footprint not to exceed the subcatchment’s area; when it comes within 0.1% of it, validation snaps the two to be equal. Direct rainfall always lands on the unit, but run-on from upstream subcatchments reaches LID units only when the LID footprint occupies the entire subcatchment — an exact-equality test, which is why the 0.1% snap matters; otherwise run-on bypasses all LIDs. Each unit takes an initial saturation percentage that pre-fills its soil and storage layers (and correspondingly shrinks the soil Green–Ampt deficit). Surface overflow joins subcatchment runoff, exfiltration joins infiltration, and underdrain flow is tracked separately, routable to its own subcatchment or node — defaulting to the parent subcatchment’s outlet; drain flow to a node is interpolated between runoff-step values at each routing step, while drain flow to a subcatchment arrives one runoff step delayed. Independently, a unit’s entire outflow can be returned onto the pervious area (surface flow always; drain flow only when its destination is the subcatchment’s own outlet) — silently disabled at validation when the parent subcatchment is at least 99.9% impervious, there being no pervious area to receive it. An optional per-unit detailed report file logs eight flux rates and four storage levels each runoff step, compressing dry spells to their boundary records.

10.6 Water Quality in LIDs

Water quality in LIDs is volume-based: outflows carry the subcatchment’s computed washoff concentration unchanged (with mixing corrections for direct-rainfall loads), so load reduction is proportional to runoff reduction — full capture is 100% removal. No media treatment chemistry is represented, though 5.1.013+ accepts per-pollutant percent removals applied to underdrain loads only (an empirical credit, not process chemistry).

11. Control Rules

11.1 Premises

Rules are RULE name / IF premise / {AND|OR premise}* / THEN action / {AND action}* / [ELSE action*] / PRIORITY p blocks, parsed by a strict state machine. Premises take the form object id attribute relop value-or-reference: objects are gages, nodes, links (conduit/pump/orifice/weir/outlet), or the simulation itself; attributes include node depth/max-depth/head/volume/inflow — the last being the node’s lateral inflow, not its total — link flow/depth/velocity/status/setting (conduits adding full-flow/full-depth/length/slope, new in 5.2), a link’s time-open/time-closed, a gage’s current intensity or its past-n-hours rainfall (up to 48 h, summed over completed hourly buckets so the partial current hour is excluded — and always zero for a gage no subcatchment references), and simulation time/date/clock-time/day/month/day-of-year. Values compare in user units; every time-valued comparison (elapsed time, clock time, time-open/-closed) carries a half-step tolerance window for both = and <>. Attribute applicability is enforced at evaluation, not parse: velocity and the conduit attributes return “missing” for non-conduits, status only for conduits and pumps, and setting only for pumps/orifices/weirs (an OUTLET … SETTING premise parses but is silently always false) — a missing operand makes the premise false, never an error. One consistency check does run at parse time, and only warns: a premise comparing two object attributes of different kinds — a node depth against a link flow, say — is accepted with a warning rather than refused, and the check is skipped entirely when the left-hand side is a named expression. Boolean evaluation is sequential with short-circuiting: an OR premise is evaluated only when the running result is false, so it disjoins with the immediately preceding premise — A AND B OR C evaluates as A AND (B OR C), not conventional precedence. SWMM 5.2 adds named VARIABLE and EXPRESSION declarations usable as premise left-hand sides.

11.2 Actions and Conflict Resolution

Actions set link controls only: conduit open/closed, pump on/off or speed setting, orifice/weir/outlet setting in $[0,1]$ — as a constant, a CURVE lookup (evaluated at the last-compared premise’s left-hand value — the same premise that supplies the PID set-point, though the controller reads its right-hand side), a TIMESERIES lookup at the current time, or a PID controller. Conflicts resolve through a per-link pending-action slot where a strictly higher rule priority replaces, ties keeping the earlier rule; an action “fires” only if it actually changes the link’s target setting (the changes count feeds steady-state detection), and modulated actions are excluded from the report’s action log. Rules are normally evaluated every routing step; a RULE_STEP option restricts evaluation to a fixed clock (§2).

11.3 PID Controllers

The PID form interprets its three parameters as gain $K_p$, integral time $K_i$ (minutes), and derivative time $K_d$ (minutes), applying a velocity-form update on the normalised error $e = (x_{sp} - x)/x_{sp}$ (normalised by the controlled value instead when the set-point is zero; integral term dropped when $K_i = 0$):

$$\Delta u = K_p\left[(e_0 - e_1) + \frac{e_0,\Delta t}{K_i} + K_d,\frac{e_0 - 2e_1 + e_2}{\Delta t}\right]$$

added to the link’s target setting each step, floored at 0 for all links and capped at 1 for non-pumps — with small-error dead-banding, and a “stuck-controller” reset zeroing the error history when successive errors differ by less than $10^{-4}$. The set-point is taken from the rule’s last-compared premise — so the premise both triggers the rule and defines what the controller regulates toward.

11.4 The Expression Language

VARIABLE and EXPRESSION declarations, treatment equations (§9.4), and custom groundwater flow and deep-percolation relations (§3.4) are all evaluated by one tokenizing recursive-descent parser, which compiles each expression once into a binary tree and thereafter walks it against a 1024-deep evaluation stack. Its grammar is the conventional three-level precedence

$$E \rightarrow E,{+,|,-},T,\qquad T \rightarrow T,{\times,|,\div},F,\qquad F \rightarrow F\ \verb|^|\ F \mid \text{fn}(E) \mid (E) \mid \text{literal} \mid \text{variable},$$

with unary minus admitted only where a previous token is absent or an operator, and exponentiation binding tighter than multiplication. Literals accept scientific notation. Names resolve through a caller-supplied binding table, which is what lets the same evaluator serve three unrelated variable vocabularies.

Nineteen functions are recognised, case-insensitively: sin cos tan cot asin acos atan acot sinh cosh tanh coth abs sgn sqrt log log10 exp step. Evaluation is total — no domain error can propagate to the caller. Square roots and logarithms of non-positive arguments return zero, as does a power with a non-positive base, division by zero, and cot/coth of zero; and any result that is not equal to itself (that is, NaN) is replaced by zero at the top of the walk. step returns 0 for a negative argument and 1 otherwise. The practical consequence is that an ill-posed treatment or groundwater expression silently evaluates to zero rather than announcing itself, so a mistyped relation reads as “no flux” rather than as an error.

12. Continuity Accounting

12.1 The Five Balances

Five balances are tallied over the whole run. Each reduces to the same error statistic — writing $\mathcal{I}$ for the accumulated inflow side and $\mathcal{O}$ for the outflow side,

$$\varepsilon = \begin{cases} \approx 0, & |\mathcal{I} - \mathcal{O}| < \tau,\[2pt] 100\left(1 - \dfrac{\mathcal{O}}{\mathcal{I}}\right), & \mathcal{I} > 0,\[6pt] 100\left(\dfrac{\mathcal{I}}{\mathcal{O}} - 1\right), & \mathcal{I} \le 0 < \mathcal{O}, \end{cases}$$

with the agreement threshold $\tau = 1$ ft³ for the volumetric balances and 0.001 mass units for the two quality balances. The third branch is the sign-preserving mirror used when a balance has outflow but no inflow.

Runoff, over the subcatchment surfaces and their snow packs:

$$\underbrace{V_{rain} + V_{runon} + V_{pond,0} + V_{snow,0}}{\mathcal{I}} ;;\text{vs.};; \underbrace{V{evap} + V_{infil} + V_{runoff} + V_{drain} + V_{plow} + V_{pond,f} + V_{snow,f}}_{\mathcal{O}}$$

where $V_{drain}$ is LID underdrain discharge and $V_{plow}$ the snow ploughed out of the system.

Groundwater, over the aquifers:

$$V_{infil} + V_{gw,0} ;;\text{vs.};; V_{ET,U} + V_{ET,L} + V_{perc,deep} + V_{gw,lat} + V_{gw,f}.$$

Flow routing is the one balance with signed terms. Wet-weather and RDII inflows and the initial stored volume always sit on the inflow side; final storage, flooding, evaporation and seepage always on the outflow side; but dry-weather, groundwater and external inflows, and the system outflow, are each placed by their sign — a net-negative external inflow crosses to the outflow side and a net-negative outflow crosses to the inflow side:

$$V_{sys,0} + V_{ww} + V_{rdii} + \sum_{k}\max(V_k, 0) - \min(V_{out},0) ;;\text{vs.};; V_{sys,f} + V_{flood} + V_{evap} + V_{seep} + V_{react} - \sum_{k}\min(V_k, 0) + \max(V_{out},0)$$

for $k \in {\text{dry-weather},\ \text{groundwater},\ \text{external}}$. Two notes on this one: the $V_{react}$ slot exists but is never accumulated, so reaction losses appear only in the quality balance; and under steady-flow routing the final storage omits link volumes while the initial storage includes them, which shows up as an apparent loss.

Quality, per pollutant, with the worst $|\varepsilon|$ reported as the run’s figure:

$$M_{0} + M_{dw} + M_{ww} + M_{gw} + M_{rdii} + M_{ex} ;;\text{vs.};; M_{flood} + M_{out} + M_{react} + M_{seep} + M_{f}.$$

Unlike the flow balance this one is not sign-split, because the signed case is handled at accumulation time: a negative outflow mass is re-credited into $M_{ex}$ as it occurs. Count-unit pollutants report as $\log_{10}$, and mass totals convert to user units only at reporting.

Surface loading, the buildup/washoff ledger that sits upstream of the quality balance:

$$B_{0} + B_{buildup} + B_{deposition} ;;\text{vs.};; B_{sweep} + B_{infil} + B_{bmp} + B_{washoff} + B_{f}.$$

A per-step flow error (rates in vs. rates out) feeds the steady-state skip decision (§2), not the time-step diagnostics. A continuity table prints when its error exceeds 10% — a signed comparison, so a large negative error does not by itself trigger it — or when the CONTINUITY report option is on, which it is by default; the 10% test therefore only decides anything for a run that has explicitly switched CONTINUITY off. Per-node cumulative inflow/outflow volumes are also tracked (initial volume seeding inflow; outfalls and terminal nodes counting inflow as outflow; final volume added to outflow), driving the node-inflow summary’s flow-balance column and the “highest continuity errors” top-five.

12.2 Reported Statistics

Alongside the balances, the engine accumulates summary statistics on every routing step (“hours” quantities are step-time integrals; maxima carry occurrence dates). The per-object statistics — everything through the pump entry below — are gated on the report start date and see nothing before it; the numerical-performance statistics of the last two entries are not gated and span the whole run:

• per-subcatchment water-balance totals, peak runoff, and runoff coefficient;
• groundwater flux totals and time-weighted average moisture/water-table;
• per-pollutant washoff loads;
• node average/maximum depths (plus a separate maximum sampled only at reporting intervals), flooding (hours, volume, peak overflow, peak ponding — a node “floods” when over full volume or overflowing), and surcharge (dynamic wave only; hours and clearances above crown/below rim);
• storage average/max volumes and losses;
• outfall flow-frequency, average/max flow, and total loads (plus the system-wide maximum simultaneous outfall flow);
• link maxima (|flow|, velocity, depth, capacity ratios);
• conduit time-in-flow-class across the seven dynamic-wave classes, hours normal-flow-limited, under inlet control, at full flow, capacity-limited, and full at either end;
• pump utilisation, startup count, min/avg/max flow, volume, energy (kWh), and time off each end of its curve — a pair of columns the engine evidently intends to fill for every pump type, since the summary prints both unconditionally, but which only a Type 4 pump distinguishes. Type 4 sets the flow class to the dedicated upstream-dry and downstream-dry markers according to which end of its curve was exceeded; the other four types set a generic off-curve flag whose numeric value collides with the upstream-dry marker. The statistics reader tests those markers, so for Types 1, 2, 3 and 5 all off-curve time is booked to the high end and the low-end column is structurally always zero;
• routing time-step min/avg/max with a log-binned frequency table, average iterations, percent non-converging, and percent of time in steady state; and
• top-five “highest” lists — node continuity errors, Courant-critical elements (counted as occurrences of being the step-limiting element), flow-instability indices, and most-frequently non-converging nodes.

The report file additionally carries a rainfall-file summary, an RDII sewershed summary, the control-actions log, an options echo, and (on request) per-object time-series tables — a text channel separate from the binary output.

13. Units and Physical Constants

Unit system selection: the user’s flow unit selects the entire unit system, for every quantity:

Internally computation runs in feet, square/cubic feet, cfs, and °F, with time in seconds (and dates as decimal days since 1899-12-30); conversion happens at input parsing and output writing, through a fixed factor table, and — in the handful of places catalogued in §16 — around an interior expression that is deliberately evaluated in the user’s units instead. Key examples:

Physical constants, together with the localised exceptions to the internal-units convention (the full list is catalogued in §16):

The various structure coefficients of §7 complete the constant set. Concentrations are mg/L, µg/L, or counts/L regardless of system.

Analysis-option defaults. Roughly forty scalar options are initialised before any input is read, so an INP file that omits a keyword still gets a definite value. The ones that change results:

A model that specifies no dates starts at 1 January 2004.

14. Input and Output

Input

File Grammar

Before any section means anything, the reader imposes a uniform lexical layer. A line is at most 1024 characters; a longer one is an error unless the overflow lies entirely past a semicolon, since the length is re-measured up to the first ; before the check. Everything from a ; to end of line is a comment and is cut before tokenising, so a comment may follow data on the same line and a line whose first token begins with ; is skipped whole. Tokens are separated by spaces, tabs, carriage returns and newlines, at most 40 per line; a token beginning with a double quote runs to the next quote or newline, which is the only way to carry a separator inside a value. Section headers are the tokens beginning [, matched case-insensitively and by prefix. An unrecognised bracketed token raises a keyword error and leaves the reader sectionless, silently discarding subsequent lines until the next recognised header. Reporting stops after 100 errors.

Prefix Matching

Prefix matching is not confined to section headers: one routine performs every keyword lookup in the file — section names, option names, and the enumerated values of options alike — by walking a keyword table in order and returning the index of the first entry that is a prefix of the token. Three consequences follow, and all three are part of the file contract rather than of any one section’s syntax:

• Trailing characters are ignored; truncations are not accepted. DYNWAVE is matched by the token DYNWAVEXYZ, which is accepted as dynamic-wave routing; the truncation DYN is rejected, because the comparison runs to the end of the keyword and fails on the token’s terminator. Matching is therefore strictly looser than equality: the token space each keyword accepts is unbounded to the right.
• Table order is load-bearing wherever one keyword prefixes another. The library contains exactly two such pairs, and they are ordered oppositely. Among section names [INLET prefixes [INLET_USAGE, and the table lists [INLET_USAGE first, so the longer name wins; the ordering is therefore part of the format’s meaning rather than an incidental property of the table. Among [REPORT] keywords NODE prefixes NODESTATS and is listed first, so the shorter name wins; the consequence is worked through under Reporting Directives below.
• Prefix and exact matching coexist within a single routine. Value tokens compared against ALL and NONE in [REPORT] go through a full-string comparison, not the prefix matcher, so ALLNODES is a node name there while DYNWAVEXYZ is a routing model in [OPTIONS].

The file is read twice. The first pass recognises sections and registers the identifier on each line, counting objects by type and rejecting duplicates; the second rewinds and parses the data. Two consequences follow, and both are properties of the file format rather than of the reader: forward references are legal — a conduit may name nodes defined later, because every identifier exists by the time parsing begins — and object counts are fixed before any parameter is read, so a section’s identity is established by its first token alone. [TITLE] accumulates up to three lines.

Sections

The 57 recognised section keywords, what each configures, and where this document treats its semantics:

The remaining nine — [MAP], [COORDINATES], [VERTICES], [POLYGONS], [SYMBOLS], [LABELS], [BACKDROP], [TAGS], [PROFILES] — are display metadata for the desktop interface. The engine recognises them and discards their contents: they are listed in the keyword table, so they neither raise an error nor swallow the sections that follow, but the parser has no case for them and reads nothing. They are therefore preserved verbatim through a load-and-save cycle without ever being interpreted.

Reporting Directives

[REPORT] mixes two grammars. Ten keywords are recognised; seven (DISABLED, INPUT, CONTINUITY, FLOWSTATS, CONTROLS, AVERAGES, and the dead NODESTATS) take a yes/no argument, while SUBCATCHMENTS, NODES and LINKS take either ALL, NONE, or a list of identifiers, each of which sets that object’s individual report flag and the object type’s flag to “some”. The type flags govern which objects reach the binary output file, so a [REPORT] line changes results the file carries, not merely their presentation.

NODESTATS is the casualty of the prefix rule. The reader has a branch for it that returns success without acting — the keyword is deprecated and meant to be tolerated and ignored — but NODE precedes it in the table and is a prefix of it, so the branch is unreachable. A NODESTATS line is parsed as a NODES line, and its argument as a list of node identifiers: NODESTATS YES raises an undefined-node error on the token YES, aborting the run, while NODESTATS ALL silently switches per-node reporting on for the entire network and enlarges the binary output file. Neither outcome is the intended tolerate-and-ignore: the branch expressing that intent is present in the source but cannot be reached.

Deprecated and Inert Options

Several [OPTIONS] keywords parse successfully and then do nothing. They remain part of the accepted grammar, so a file carrying them runs exactly as one without them.

FLOW_ROUTING also accepts a legacy alias table — NONE, NF, KW, EKW, DW — consulted only when the token matches no modern name. The two tables are matched positionally rather than by meaning, so NF (SWMM 4’s normal-flow routing) selects steady flow, KW kinematic wave, EKW the reassigned-to-kinematic extended model, and DW dynamic wave.

FLOW_ROUTING NONE is not a fourth routing model: it sets the same ignore-routing flag the IGNORE_ROUTING option sets and leaves the routing model at its default of dynamic wave, so dynamic-wave validation still runs over a model that will never be routed.

Process Switches and Time-Step Interlocks

Global process switches (IGNORE_RAINFALL, IGNORE_SNOWMELT, IGNORE_GROUNDWATER, IGNORE_RDII, IGNORE_ROUTING, IGNORE_QUALITY, plus the equivalent FLOW_ROUTING NONE) disable whole subsystems — quality ignoring also strips all pollutant variables from the binary output — and subsystems with no objects are ignored automatically. Time steps interlock at validation: the report step must be ≥ the routing step (fatal otherwise), the dry step is raised to the wet step, and the routing step is clamped to the wet step. Date/time conventions: INP dates are M/D/Y with -// separators (3-letter month names accepted), times decimal-hours or h:m:s; decoded times round to the nearest second; every conversion of elapsed time to a calendar date adds +1 ms (so each reporting timestamp and date-driven lookup sits 1 ms past nominal); and elapsed-time labels measure from report start.

Semantic Validation

Beyond the topology rules of §4 and the per-object parameter checks described alongside each object, validation enforces a set of cross-object consistency rules whose violation is fatal. They are worth stating together because each marks a model the engine refuses outright rather than silently accepting:

Errors occupy two disjoint numbering spaces. The core’s own catalogue runs 101–509 across some two dozen themed bands, one per subsystem or file type (runtime, subcatchment/aquifer, conduit/pump, topology, node, RDII, rain gage, treatment, curve/series, snowmelt, LID, date/time, input parser, then one band per interface-file kind, and finally API), and the OWA fork adds a parallel toolkit catalogue at 2000–2013 (§16). Both are internally consistent: every declared code carries exactly one message, with no duplicates and no orphaned text, in either space.

Eight core codes are nonetheless unreachable, for three distinct reasons. One is a demotion: the code for a conduit whose elevation drop exceeds its length survives, though that case is now handled as the warning and fallback of §1.4 instead. Six are displacement — the core’s own API errors for an invalid object type, index or name, an invalid property type or value, and an invalid time period are all superseded by the toolkit set, which the fork raises in their place while leaving the originals declared. The last is an error nothing can detect: the code for a failed read of the binary results file is unreachable because none of the four readers checks its read at all — each seeks to a computed offset and reads, ignoring the returned count — so a truncated or corrupt file yields stale buffer contents to swmm_getSavedValue (§15) rather than an error. All eight remain part of the published catalogue, numbered and messaged, even though nothing emits them.

Interface Files

Ancillary interface files carry data between runs. They are declared in [FILES], one line of mode, type, name, and the declaration grammar is less uniform than the four-state pattern of §16 suggests. Four modes are spellable — NO, SCRATCH, USE, SAVE — though the format is conventionally described as USE/SAVE only; on a SCRATCH line the supplied name is discarded in favour of a generated temporary one that is deleted when the run closes. Six types are recognised, and their mode constraints differ:

Three asymmetries follow. Every type but HOTSTART keeps a single slot, so a second line naming the same type overwrites the first rather than adding to it. The INFLOWS/OUTFLOWS mode restrictions are enforced, but reported as a wrong-number-of-items error rather than a bad keyword, which points a reader at the wrong token. And file names are made absolute against the input file’s directory for every type except RDII, whose name alone is left relative to the process’s working directory — so the identical [FILES] line resolves to different paths for different types depending on where the program was launched. A line carrying only a mode and a type, with no name, is accepted and silently does nothing.

The rainfall interface file (binary) collates external files into per-station records of (date, depth) pairs for non-zero periods only. Its layout is a 10-byte SWMM5-RAIN stamp, a gage count, then a per-gage index of station ID (80 bytes), recording interval in seconds, and the first and one-past-last byte offsets of that gage’s data; the data themselves are, per gage and per non-zero period, an 8-byte date followed by a 4-byte rain depth in inches whatever the model’s unit system. Gages are matched by station ID (shared stations share data; one station in two files is fatal), file-fed gages become volume-type at the file’s interval (NWS/Canadian formats override the declared interval and shift end-of-interval stamps), NWS accumulation codes split totals evenly across their span, and decreasing cumulative readings reset the accumulator. The routing interface file is plain text, and its header is positional rather than keyed:

SWMM5 Interface File
<title line>
<reporting time step in sec>  - reporting time step in sec
<n constituents>              - number of constituents as listed below:
FLOW <flow units word>
<pollutant id> <concentration units word>          × (n − 1)
<n nodes>                     - number of nodes as listed below:
<node id>                                          × n
Node             Year Mon Day Hr  Min Sec FLOW       <pollutant ids…>

after which each reporting period contributes one line per node — identifier, then year, month, day, hour, minute, second as integers, then flow in the file’s own units, then one concentration per pollutant. The leading integer on each of the two count lines is what the reader consumes; the trailing prose is decoration. On reading, values are linearly interpolated between bracketing periods, nodes and pollutants are matched by name with unmatched pollutants taken as zero, and flows are converted from the file’s declared units rather than the current model’s. Outflows are written only for outlet nodes, and one file cannot serve as both inflow and outflow in the same run. The hotstart file (SWMM5-HOTSTART4 stamp; versions 1–4 readable, with older versions carrying progressively less state) checkpoints approximately the §1.7 state vector (see the caveats there) — runoff state as doubles, routing state as floats, link settings re-applied through the control machinery — and its compatibility check covers object counts and flow units only: a reordered model silently loads the wrong state. The buildup block is written and read asymmetrically. Reader and writer evidently intend the same layout — one value per (land use, pollutant) pair followed by the last-swept date — and the reader implements it; the writer instead emits, for each pollutant, a record whose length is the pollutant count, drawing from a six-element scratch buffer of which only the first element was set. The two agree only when a model carries exactly one pollutant. Beyond that the file is longer than the reader expects and its buildup values are meaningless, and beyond six pollutants the writer also reads past the end of its buffer. A hotstart file therefore never round-trips surface buildup for a multi-pollutant model.

Output

Output has two faces: a text report file written for a reader, and a binary .out file written for a program.

Report File Structure

The report is not assembled at the end — it is emitted in phases as the run proceeds, so its section order is fixed by the lifecycle of §15 rather than by any layout choice. On open: the logo and version banner, then, once the input has been read, the title lines and — only if the INPUT report option is set — an echo of the parsed network, object by object. On start: the options summary. During initialisation: the rainfall-file summary, once per station, and the RDII sewershed summary. During the run: the control-actions log, one line per action that actually changed a link’s target setting, excluding curve-, time-series- and PID-modulated actions (§11).

On end the bulk arrives, in this order:

1. the five continuity balances of §12 — runoff, surface loading, groundwater, flow routing, quality — each printed when its error exceeds 10% or, as is the default, whenever the CONTINUITY option is on;
2. the numerical-performance block, printed when the FLOWSTATS option is on (also the default): the four top-five lists (node continuity errors, Courant-critical elements, flow-instability indices, most-frequently non-converging nodes) followed by the time-step summary and its log-binned histogram;
3. the per-object summary tables, in the fixed order below.

Each table draws on the statistics of §12.2 and reports one row per object. Suppression is by group rather than by table, and the grouping is not what the contents suggest. The first four tables are gated on the model having subcatchments and on rainfall not being ignored — unless snowmelt or aquifer objects exist and their own subsystems are active — so IGNORE_RAINFALL alone silences the groundwater and LID summaries too. All the remaining tables, node tables included, sit inside a single gate on the model having links and routing not being ignored: a model of nodes with no links between them prints no node depth, inflow, flooding, storage or outfall table at all. Two tables are further restricted to dynamic-wave routing — Node Surcharge and Flow Classification — the latter because the flow classes it tabulates are only assigned by that solver.

The “adjusted/actual length” column of the flow-classification table is the conduit-lengthening ratio of §5.6, reported so that a user can see which conduits the stability transform altered.

Each group is suppressed when its objects are absent or its subsystem ignored, and the whole report can be switched off. Finally, and only on an explicit request after the run has ended, the per-object time-series tables are read back out of the binary file and printed for whichever subcatchments, nodes and links were flagged in [REPORT].

Binary Output File Layout

The .out file is written as the following record sequence:

Node and link values are period-interpolated, or period-averaged on request. The average is an unweighted mean of the end-of-step values sampled over the period, not a time-weighted integral — so under a variable routing step the short steps carry the same weight as the long ones. Regulator and pump settings are exempted from it (the accumulator is overwritten with the latest value scaled by the step count, so the period reports the final setting rather than a mean), and pump flows are not interpolated across on/off transitions; subcatchment values are always interpolated and system values are current-step totals — all already in user units. Two reader caveats: per-object results appear only for objects flagged in [REPORT] (all off by default — an unconfigured run’s binary file holds only the 15 system series), and when the report start postdates the simulation start, the stored start-date field is deliberately backdated one period before the first record.

15. The Engine as a Library

SWMM is an embeddable shared library; the command-line tool is a thin progress-callback wrapper over the same public API. The surface comes in two strata that a consumer sees as one: a core of 20 exported functions carried in the EPA lineage, and the community fork’s toolkit of 64 more (§16) — plus two unexported state predicates the toolkit uses internally to enforce the phase rules below. The core is deliberately minimal, a lifecycle plus a property-code accessor pair and error reporting; the toolkit provides the typed, per-object surface most integrators actually use.

Run-Loop Lifecycle

swmm_open (parse and validate) → swmm_start(saveResults) (initialise state; collate the rainfall interface file; pre-compute RDII; read any hotstart file; saveResults = FALSE skips the binary output) → repeated swmm_step, each advancing exactly one routing step and returning elapsed decimal days, 0.0 signalling completion — or swmm_stride(seconds), which advances a fixed span of simulation time by temporarily capping the routing step — → swmm_end → optional swmm_report (write the text report) → swmm_close. The window between open and start is where pre-run modification is legal; the run-total mass-balance errors are queryable only between end and re-start. swmm_run collapses the whole sequence into one call for the batch case, and every entry point is guarded by open/started state checks that return an error code rather than faulting.

The lifecycle divides the surface into four timing classes, and this is the part an embedder must get right, because the classification is enforced per property rather than per function:

Attempting a write outside its phase returns an error rather than silently taking effect, so the boundary is discoverable at run time. The asymmetry is the point: geometry is frozen once the run starts, forcing is not.

Object and Property Access

The core exposes a single generic accessor pair keyed by object type and property code, plus swmm_getCount, swmm_getIndex and swmm_getName for enumeration and identifier resolution. The toolkit layers a typed surface over the same state, organised by object rather than by code:

• enumeration and identity — object counts by type, identifier from index and index from identifier, link connectivity and direction, subcatchment outlet connection, node and link type codes;
• parameters — get/set per object class for subcatchments, nodes, links and inlets, each under the timing rules above;
• results — current-time values per object, including per-pollutant concentrations for subcatchments, nodes and links;
• statistics — the accumulated node, storage, outfall, link and pump summaries of §12 read out mid-run, and the system-wide runoff and routing totals; these allocate, so the toolkit pairs them with explicit free functions;
• simulation control — analysis settings, unit system, simulation parameters, and the simulation date-time, both readable and settable;
• LID — process- and unit-level parameters, options, per-unit results, and the eight flux rates of §10’s detailed report.

Error reporting is uniform: every toolkit entry point returns an integer code drawn from the fork’s own 2000–2013 catalogue (§14) — fourteen codes covering out-of-bounds arguments, a project not open, a simulation running or not running, wrong object types, and bad object, pollutant, time-series, pattern, LID-unit and inlet indices — with swmm_getAPIError turning it into text, distinct from the core’s swmm_getError and swmm_getWarnings which report on the run rather than on the call.

Query and Mutation

swmm_getSavedValue re-reads an object’s binary-file results by reporting period after swmm_end — the one accessor that reads the output file rather than live state, and a link’s setting is served from the file’s capacity slot since the two share a column (§14). It reaches only objects the [REPORT] section selected, and returns 0 rather than an error for any other, because the lookup goes through the object’s report flag: that field is a boolean while the input is parsed and validated, then rewritten at initialisation to the object’s one-based position in the binary output file, zero meaning absent. Object order in the file is therefore project order filtered by the report flags, and a [REPORT] line silently changes both which objects this accessor can see and what the file contains. Mid-simulation setters change boundary forcing and controls while the model runs: gage rainfall override (taking precedence over every data source — except that a gage deferring to a shared co-gage adopts the co-gage’s value before its own override is consulted), node lateral inflow, outfall stage (converting the outfall to fixed-stage), and link target settings (conduits excluded, with no report logging; the OWA toolkit’s separate swmm_setLinkSetting is the converse — it logs each application to the report as a virtual “ToolkitAPI” rule but does not exclude conduits). Setting the routing step mid-run silently zeroes the Courant factor — disabling variable time-stepping for the remainder of the run, a side effect no error code announces.

Rainfall Injection and Checkpointing

Two distinct rainfall mechanisms exist: the property-based intensity override, and swmm_setGagePrecip, which converts a gage to the RAIN_API data source fed each step by the caller. swmm_hotstart both loads a state file before start and saves an on-demand checkpoint mid-run — the programmatic complement of the [FILES] hotstart directives.

16. Cross-Cutting Engine Contracts

The preceding sections follow SWMM’s physical subsystems; this one collects the engine-wide contracts that span them — behaviours visible only in the code’s architecture, each referenced from the sections carrying its fragments.

The unit-boundary contract. All internal computation is US customary (ft, ft², cfs, °F, seconds); conversion factors (UCF, the flow-unit table) are applied only at boundaries — input parsing, report writing, binary output, and the API get/set surface (§13). The exceptions — places where an interior computation is carried out in the user’s units and its coefficients therefore change meaning with the flow-unit selection — form a longer list than the manual acknowledges:

• the LID underdrain equation (§10);
• the groundwater lateral-flow power function and any user-written groundwater expression (§3.4);
• the weir discharge equations, including the head-dependent coefficient curve (§7.3), the roadway weir excepted;
• the outlet power function and rating curve (§7.4);
• the storage-unit area, volume, and depth relations (§6.3);
• the weir and tabular divider rules (§1.3);
• the treatment expression language, whose FLOW, DEPTH, and AREA variables are all served in user units (§9.4).

The HEC-22 inlet equations additionally use their own $g = 32.16$ ft/s² (§8.4), and the roadway weir keeps internal feet but rescales a user-supplied discharge coefficient by $1/0.552$ under SI (§8.3).

The old/new state discipline. Nearly every dynamic quantity is stored as an old/new pair, rolled at each step’s start; results read out at intermediate times by weighted interpolation between the pair. This is the mechanism that reconciles the three clocks of §2 — runoff results interpolate onto routing times, and routing results onto reporting times. Four object families carry an explicit “roll the old state” operation — subcatchments, links, nodes, and LID groups — one per level at which interpolation is read out. The continuously-integrated subsystems (infiltration, groundwater, snow) keep no old/new pair at all: their previous state is implicit in the integrator’s starting condition, which is why they are reported at end-of-step values rather than interpolated (§2).

The state-serialization schema. Three modules expose paired get-state/set-state vectors — snowpack, groundwater, and infiltration, the last dispatching to whichever of the five methods a subcatchment uses. These vectors are the persistence schema and, once defined, a compatibility contract. Their sole consumer is the hotstart file, whether written at the end of a run or as a mid-run checkpoint through the API (§14, §15); everything else in a hotstart — sub-area ponded depths, buildup, node and link state — is read and written field by field rather than through a vector. The runoff interface file is not a client of this schema despite sounding like one: it caches the same eight reporting variables per subcatchment that the binary output carries, so it replays results, not state, and a run resumed from it starts with cold antecedent conditions.

Validation as mutation. The validation pass does not merely check the model — it rewrites it: node maximum depths raised to link crowns, regulator crests raised to downstream inverts under dynamic wave only — the same model merely earning a warning under steady or kinematic routing, so the routing option silently selects between two different networks — conduit slopes floored or adverse-slope conduits reversed, offsets converted between conventions, elevations snapped, infeasible shape radii enlarged, street sections compiled into transects, and equivalent lengths/surface areas computed for orifices (§1, §4, §6–§8). A behaviourally-faithful reading of an INP file is the post-validation model, not the literal text.

Alongside the rewrites runs a thinner but equally binding class: out-of-range values that are clamped rather than refused. A subcatchment’s imperviousness is capped at 100% (§1.2); a curve number is confined to $[10, 99]$ and a monthly conductivity adjustment of zero or less becomes 1 (§3.1, §3.3); an initial snow free-water depth is capped at the pack’s holding capacity (§3.5); site elevations at or below sea level bypass the atmospheric-pressure fit (§13); orifice and weir equivalent lengths are floored at 200 ft (§7); an LID unit’s return-to-pervious routing is switched off for a ≥99.9% impervious parent (§10); and the report start is advanced to the simulation start, the lengthening step floored at zero, and the thread count capped at the machine’s maximum (§2, §13). The distinction is consequential: a rejected value produces an error the modeller sees, whereas a clamped one produces a different model that runs to completion.

Option plumbing. Forty-five scalar analysis options are parsed in one place and consumed deep inside distant modules (§2, §5, §14). The IGNORE_* family silently amputates entire subsystems, and subsystems with no objects are ignored automatically — so the effective process set is a joint function of options and object counts.

Vestigial state. A small, bounded residue of abandoned features survives in the parser and the object model, and it is bounded tightly enough to enumerate exhaustively: two options parse into globals nothing reads (SLOPE_WEIGHTING and COMPATIBILITY), one routing keyword is recognised and immediately reassigned (XKINWAVE → kinematic wave), one [REPORT] keyword is deprecated but unreachable behind a prefix collision (NODESTATS, §14), and exactly one field of the entire object model is written and never read — a conduit’s supercritical-flow flag, computed at initialisation for every conduit from a comparison of normal and critical velocity scaled by a factor of $0.3$ that the source itself marks as belonging to a modified kinematic-wave routing no longer present. None of this affects results; all of it nonetheless parses, so its presence in a file is invisible in the results.

The interface-file lifecycle. Rainfall, runoff, RDII, hotstart, and routing files all follow one four-state mode pattern — none / scratch / use / save — decoupling expensive stages so they can be precomputed once and replayed (§14). This is a single architectural idea, not five unrelated formats. The idea is uniform; its realisation is not, each type accepting a different subset of the four states and one resolving its path differently from the rest (§14).

The numerical toolkit. Four solution devices recur everywhere. Bracketed Newton–Raphson appears through a shared root-finder in the kinematic-wave continuity solve (§4), storage depth-from-volume (§6.3), and the section-factor inverse (§6.1), and open-coded inside the Horton and Green–Ampt cumulative-infiltration solves (§3.3). Ridder’s method takes the two places where the function is too awkward to differentiate — critical depth (§6.3) and the culvert form-1 energy equation (§8.3). Adaptive Cash–Karp RK5 integrates exactly two systems, the sub-area ponded depths (§3.2) and the groundwater pair (§3.4), through one shared integrator. Its step controller is the standard embedded-pair form: the fifth- and fourth-order solutions differ by $\varepsilon_{max}$, scaled by the caller’s tolerance, and the step is rescaled by $0.9,\varepsilon_{max}^{-1/4}$ on rejection and $0.9,\varepsilon_{max}^{-1/5}$ on acceptance, clamped so a single adjustment never shrinks the step below a tenth or grows it beyond five times its current value, and never advances past a 10,000-step ceiling. Tolerances are $10^{-4}$ ft on ponded depth and, for groundwater, on the moisture/depth pair. Picard successive approximation with under-relaxation carries the dynamic-wave solve (§5), the steady/kinematic storage-node balance (§4), and the LID swale (§10). A fifth device, stateful bracketed table lookup, underlies all tabulated geometry (§6). Curves and time series share one table representation with two semantic modes — sorted x-lookup versus date-cursored streams — and the cursor-stateful lookups are not thread-safe.

The threading model. OpenMP parallelism exists in exactly one place — the dynamic-wave per-conduit and per-node loops (§5) — and is off unless asked for, the thread count defaulting to 1; nodal flow accumulation stays serial, so results are thread-count-invariant. Everything else is strictly sequential and order-dependent, most consequentially the topological routing order (§4) and the link-definition-order sensitivity of pump/regulator solves (§5).

In-loop instrumentation. Continuity accounting and statistics are not post-processing: mass-balance and statistics updates are woven through the inner loops of every physical module, each with its own definition of what counts as inflow, outflow, or loss at that point (§12). A numerically-matching re-computation of SWMM’s balances must replicate these call sites, not just the ledger formulas.

One expression language, three bindings. The tokenized math-expression evaluator (19 functions, zero-on-domain-violation semantics, §9.4) serves three unrelated features — control-rule variables/expressions (§11), treatment equations (§9.4), and custom groundwater flow equations (§3.4) — each with its own variable-binding table.

The template/instance idiom. Shared parameter sets instantiated per-consumer recur throughout: snowmelt parameter sets → per-subcatchment snowpacks, LID process designs → deployed units, unit-hydrograph groups → per-node RDII, aquifers → per-subcatchment groundwater, transects/streets → per-conduit geometry (§1.6). The instance may override selected template values (aquifers being the fullest example, §3.4).

The OWA stratum. The community fork overlays the EPA core in a marked, separable layer: the toolkit API (§15), callback-driven runs, API rainfall/inflow/quality injection paths, extra state, and structs exposed for interoperability. The extra state is worth distinguishing carefully, because one of the two obvious candidates is not an addition at all: the mixed-reactor concentrations held on nodes and links are an OWA field, but the storage node’s hydraulic residence time is EPA’s — it drives the documented HRT treatment variable of §9.4 and is persisted by the hotstart file. What OWA adds alongside it is a second, separately-named residence-time field on the node structure. The physics of §2–§12 belongs to the EPA core; the pinned tag’s identity is EPA 5.2.4 physics plus this instrumentation stratum.

Error discipline. One global error code short-circuits every phase; API entry points are guarded by open/started state checks; warnings accumulate without halting (several of §4’s silent mutations announce themselves only as warnings). Keyword tables, enum orders, and report strings are maintained in positional correspondence across three files — parsing correctness is positional — a fragile invariant of the INP grammar’s implementation, and one the prefix matcher sharpens into a hazard, since table order also decides which of two nested keywords wins (§14).

Contributing

Contributions are welcome. Please read CONTRIBUTING.md before opening a pull request, in particular the Spec First workflow, which requires spec changes to land before implementation changes for any solver, model, or analytics work.

Hydra uses conventional commit messages, expects one logical change per pull request, and requires all CI checks to pass before merge.

Testing Strategy

Hydra’s integration strategy is Hydra-native and fixture-driven:

• Purpose-built fixture networks validate control, timestep, quality, reaction, and tank behaviour through physics/behaviour invariants.
• Large real-world networks are exercised with deterministic regression checks by comparing repeated Hydra runs section-by-section.

Correctness is established by physics/behaviour invariants and deterministic Hydra-vs-Hydra regression checks, not by agreement with any external tool’s output.

License

Hydra is licensed under AGPL v3. Commercial products built on Hydra must either release their source under AGPL v3 or obtain a separate commercial license.

For details, see the project LICENSE and COMMERCIAL_LICENSE.md.