The model
A water balance of 52 sewage works, tested on a year it never saw.
How it works
Water in = water treated + water stored + water spilled. The works treats what it can; the storm tank takes the rest; when the tank is full, it spills.
How good it is
Works spilling each day, Oct 2025 – Oct 2026
Grey: recorded. Red: predicted.
Are its chances right?
On the dashed line = right.
Does the physics match reality?
Learned values vs published figures. Green band = within 20%.
Data
| Rain, river flow | Environment Agency |
| Treated flow, spill monitors | Thames Water |
| Evaporation, 7-day rain forecast | Open-Meteo, ECMWF |
| Capacity, population | Oxford Rivers Improvement Campaign |
| Rivers | OS Open Rivers |
Formulae
Every works is simulated hour by hour (\(t\)); results are reported per day (\(d\)). Flows in m³/h, rain and evaporation in mm. Fitted per works: \(F, A, A_s, \tau_f, D_{\max}, k\) and \(\tau_g\) or \(q_0\). Measured: \(C\). Sourced: \(V\).
1 · The sewage works, hour by hour
Sewage from homes and businesses: a ±30% daily cycle peaking at noon; \(h_t\) is the local hour.
$$Q^{\text{sewage}}_t = F\,\Big[1 + 0.3\cos\!\Big(\frac{2\pi\,(h_t - 12)}{24}\Big)\Big]$$Soil: the moisture deficit \(D\) (mm) is filled by rain \(P_t\) and emptied by evaporation \(E_t\), more slowly as the soil dries. Rain that falls on saturated soil recharges groundwater (\(R_t\)).
$$\tilde D_t = D_{t-1} + E_t\Big(1 - \frac{D_{t-1}}{D_{\max}}\Big) - P_t,\qquad D_t = \min\!\big(D_{\max},\ \max(0,\ \tilde D_t)\big),\qquad R_t = \max(0,\ -\tilde D_t)$$Rain runoff, delayed by \(\tau_f\) hours in the sewers. Hard surfaces \(A\) always drain; land \(A_s\) drains only as it saturates. Areas in hectares: 1 mm on 1 ha = 10 m³.
$$\bar P_t = a_f\,\bar P_{t-1} + (1 - a_f)\,P_t,\quad a_f = e^{-1/\tau_f},\qquad Q^{\text{runoff}}_t = 10\,\Big(A + A_s\Big(1 - \frac{D_t}{D_{\max}}\Big)\Big)\,\bar P_t$$Groundwater leaking into the pipes, in whichever of two forms fits each works better: a store \(G\) draining over \(\tau_g\) days, or driven by Thames flow at Farmoor \(Q^{\text{F}}_t\) (m³/s) above \(q_0\).
$$G_t = a_g\,G_{t-1} + R_t,\quad a_g = e^{-1/(24\,\tau_g)},\quad Q^{\text{gw}}_t = k\,G_t \qquad\text{or}\qquad Q^{\text{gw}}_t = k\,\max\!\big(0,\ Q^{\text{F}}_t - q_0\big)$$Treat, store or spill. Inflow \(I_t\) (\(m_t\) is the inflow error, 1 when fitting). The works treats up to capacity \(C\); the storm tank \(V\) takes the excess and is pumped back when inflow drops; only what a full tank can't hold spills. \(S_t\) tank level, \(T_t\) treated, \(X_t\) spill.
$$I_t = m_t\,\big(Q^{\text{sewage}}_t + Q^{\text{runoff}}_t + Q^{\text{gw}}_t\big)$$ $$S_t = \min\!\big(V,\ \max(0,\ S_{t-1} + I_t - C)\big),\qquad X_t = \max\!\big(0,\ S_{t-1} + I_t - C - V\big),\qquad T_t = \min(I_t,\ C) + \max(0,\ S_{t-1} - S_t)$$Capacity and tank. \(C\) is measured, not fitted: the 95th percentile of the works' daily treated flow \(T^{\text{obs}}_d\) on training days. \(V\) is the published tank volume where one exists, otherwise the design rule of 2 hours at capacity.
$$C = \tfrac{1}{24}\,\mathrm{P}_{95}\big(T^{\text{obs}}_d\big),\qquad V = 2\ \text{h} \times C$$2 · Fitting, and the chance of a spill
Fit to daily treated flow only, never to spill records. Errors are in logs, so a 10% miss counts the same at every works, with a robust loss (\(s = 0.1\)). On days at capacity (\(T^{\text{obs}}_d \ge 0.9 \times 24C\)) the meter only shows that inflow was at least that much, so only too little inflow is penalised.
$$e_d = \begin{cases} \ln \hat T_d - \ln T^{\text{obs}}_d & \text{below capacity} \\ \min\!\big(0,\ \ln \hat I_d - \ln T^{\text{obs}}_d\big) & \text{at capacity} \end{cases}\qquad \hat\theta = \arg\min_\theta \sum_d s^2\Big(\sqrt{1 + (e_d/s)^2} - 1\Big)$$Inflow error, from the fit's misses on days below capacity: its size \(\sigma\) (robust) and day-to-day persistence \(\rho\) (lag-1 correlation). Each run draws a fresh daily error.
$$\sigma = 1.4826\ \operatorname*{median}_d \big|e_d - \operatorname{median}(e)\big|,\qquad \varepsilon_d = \rho\,\varepsilon_{d-1} + \sqrt{1 - \rho^2}\,\sigma\,z_d,\quad z_d \sim \mathcal N(0, 1),\qquad m_t = e^{\varepsilon_{d(t)} - \sigma^2/2}$$Chance of a spill on day \(d\), expected spill-hours and spill volume, over \(N = 200\) runs. The 7-day forecast runs 4 errors on each of the 51 ECMWF rain-forecast members, \(N = 204\).
$$\hat p_d = \frac{1}{N}\sum_{n=1}^{N} \mathbf 1\Big[\textstyle\sum_{t \in d} X^{(n)}_t > 0\Big],\qquad \bar H_d = \frac{1}{N}\sum_{n=1}^{N} \#\big\{t \in d : X^{(n)}_t > 0\big\},\qquad \bar X_d = \frac{1}{N}\sum_{n=1}^{N} \sum_{t \in d} X^{(n)}_t$$What if. The day is rerun from the works' state that morning, with the day's hourly rain pattern scaled to the new total \(R'\) and your tank \(V\) or capacity \(C\).
$$P'_t = P_t\,\frac{R'}{\sum_{t \in d} P_t}$$3 · River water quality
Concentration in 100 m segment \(j\) on day \(d\): background, plus each works \(w\) upstream that spills that day, diluted by the river flow \(Q_j\) and decaying over the travel time \(t_{wj}\).
$$c_j = c_{\text{bg}} + \sum_{w\ \text{spilling}} c_0\,\frac{r_w}{Q_j + r_w}\,e^{-k\,t_{wj}}$$Travel time and river flow. Travel at \(u = 0.5\) m/s along the river, up to 40 km. Gauged rivers use their nearest EA gauge; elsewhere, the gauges' median runoff per km² \(q_d\) times the catchment area, estimated as 1.99 km² per km of river upstream \(L^{\text{up}}_j\).
$$t_{wj} = \frac{L_{wj}}{u},\qquad Q_j = \begin{cases} Q^{\text{gauge}}_d & \text{gauged river} \\ 1.99\,L^{\text{up}}_j\,q_d & \text{elsewhere} \end{cases},\qquad q_d = \operatorname*{median}_g \frac{Q_{g,d}}{\text{area}_g}$$Spill rate of each works, from the physics model over the test year (m³/s); 0.37 × capacity at works the model never spills. A day counts as a spill day if the monitor recorded one; on forecast days, and days with no usable monitor record, if \(\hat p_d \ge 0.5\).
$$r_w = \frac{\sum_d \bar X_d}{3600\,\sum_d \bar H_d}$$| In the spill \(c_0\) | Background \(c_{\text{bg}}\) | Loss \(k\) (per day) | Unsafe from | |
|---|---|---|---|---|
| E. coli (per 100 mL) | 1.5 × 10⁶ | 100 | 0.5 Oct–Mar, 1.1 Apr–Sep | 235: not safe to swim |
| Ammonia (mg N/L) | 4 | 0.04 | 0.3 | 0.5: harmful to fish |
4 · Statistical backup
One small logistic model per overflow, fitted to spill records. It covers overflows without treated-flow data and gives a second opinion. Rain is held in a store with memory \(\tau\), and the Thames flow at Farmoor measures how wet the catchment is; \(\Lambda(z) = 1/(1 + e^{-z})\).
$$F_t = a\,F_{t-1} + (1 - a)\,\tau P_t,\quad a = e^{-1/\tau},\quad \tau \in \{6, 12, 24, 48, 96\}\ \text{h}$$ $$z_t = b_0 + b_1 \ln Q^{\text{F}}_t + b_2 F_t,\qquad \Pr(\text{spilling in hour } t) = \Lambda(z_t),\qquad \Pr(\text{spill on day } d) = \Lambda\big(c_0 + c_1 \max_{t \in d} z_t\big)$$5 · Regulator watchlist
Test A, spilling below the permit: a day spilling 20+ hours (\(H^{\text{obs}}_d\)) while treating under 80% of the permitted flow (FFT). The permit asks for 92%; 80% leaves room for outlet metering and daily averaging.
$$H^{\text{obs}}_d \ge 20\ \text{h}\quad\text{and}\quad T^{\text{obs}}_d \lt 0.8 \times \text{FFT}$$Test B, a spill the weather doesn't explain: a recorded spill that both models, given the rain that fell, called under 5% likely.
$$H^{\text{obs}}_d \ge 1\ \text{h}\quad\text{and}\quad \hat p^{\text{physics}}_d \lt 0.05\quad\text{and}\quad \hat p^{\text{statistical}}_d \lt 0.05$$Lead level. Strong: 10+ all-day spill days with 30%+ of them failing Test A, or 10+ days failing Test B. Moderate: 5+ days failing Test A, or 3+ failing Test B.
6 · How accuracy is scored
On days the model never saw; \(y_i = 1\) if the monitor recorded a spill.
Brier score and skill: the mean squared error of the chances, and the improvement over always predicting each works' usual spill rate (climatology).
$$\text{BS} = \frac{1}{n}\sum_{i=1}^{n} (\hat p_i - y_i)^2,\qquad \text{Skill} = 1 - \frac{\text{BS}}{\text{BS}_{\text{clim}}}$$Accuracy: a spill is called when \(\hat p \ge 0.5\).
$$\text{Accuracy} = \frac{1}{n}\sum_{i=1}^{n} \mathbf 1\big[\,\mathbf 1[\hat p_i \ge 0.5] = y_i\,\big]$$AUC: share of (spill day, dry day) pairs where the spill day gets the higher chance; \(n_1\) spill days, \(n_0\) dry days, \(R_1\) sum of the spill days' ranks.
$$\text{AUC} = \frac{1}{n_1 n_0}\sum_{i:\,y_i=1}\ \sum_{j:\,y_j=0}\Big(\mathbf 1[\hat p_i > \hat p_j] + \tfrac12\,\mathbf 1[\hat p_i = \hat p_j]\Big) = \frac{R_1 - \frac{n_1(n_1+1)}{2}}{n_1\,n_0}$$Reliability: for days grouped into bins \(b\) by predicted chance, compare the mean prediction with the share that spilled.
$$\bar{\hat p}_b = \frac{1}{|b|}\sum_{i\in b}\hat p_i \quad\text{vs}\quad \bar y_b = \frac{1}{|b|}\sum_{i\in b} y_i$$Fitted only to treated flow; spill records used only for testing. River quality: sewage diluted by river flow and fading downstream; unsafe at E. coli 235/100 mL or ammonia 0.5 mg/L.