Resources Numerical accuracy Calculation policy
Builds a runoff hydrograph, a stage and storage table and an outlet rating, routes the storm through the pond, sizes an orifice, and checks the result.
64-bit arithmetic in SI units. The results depend on the time step (1 to 10 minutes), on tables that are interpolated linearly, and on empirical coefficients written to two or three digits. Feet are converted with exact international-foot factors.
12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-civilcraft-detention.mjs
The depth of runoff from a depth of rain.
S = 25400 ÷ CN − 254; Q = (P − 0.2S)² ÷ (P + 0.8S) for P > 0.2S
One closed-form expression in millimeters. 25400 and 254 are the exact metric forms of 1000 and 10 inches. The method itself is empirical.
export function runoffDepth(P, cn) {
if (!(cn > 0 && cn <= 100)) return 0;
const S = 25400 / cn - 254, Ia = 0.2 * S;
return P > Ia ? ((P - Ia) ** 2) / (P + 0.8 * S) : 0;
}Flow against time from a catchment for a storm, with its peak and volume.
Tp = Δt ÷ 2 + 0.6·Tc; qp = 0.208·A ÷ Tp; base = 2.67·Tp; Q(t) = Σ ΔQₖ·u(t − k·Δt)
Runoff increments are differences of the cumulative runoff depth, then convolved with the triangle in one sum per time step. The peak factor 0.208 and the base 2.67 are three-digit constants. Volume is the sum of flow × step, a rectangle rule, so peak and volume depend on the time step.
export function scsHydrograph({ area, cn, tc, rain, duration, dt = 5, tail = null } = {}) {
const n = Math.max(1, Math.round(duration / dt)), incr = [];
let cumP = 0, cumQ = 0;
for (let k = 0; k < n; k++) {
cumP += (rainAt(rain, k * dt) * dt) / 60;
const q = runoffDepth(cumP, cn); incr.push(q - cumQ); cumQ = q;
}
const tpH = (dt / 2 + 0.6 * tc) / 60, qp = (0.208 * (area / 1e6)) / tpH, tpMin = tpH * 60, base = 2.67 * tpMin;
const total = duration + base + (tail ?? 0), m = Math.ceil(total / dt) + 1;
const uh = (t) => (t <= 0 || t >= base ? 0 : t < tpMin ? (qp * t) / tpMin : (qp * (base - t)) / (base - tpMin)); // m³/s per mm, a pulse of one step's rain
const q = new Array(m).fill(0);
for (let i = 0; i < m; i++) { let s = 0; for (let k = 0; k < n; k++) if (incr[k] > 0) s += incr[k] * uh((i - k) * dt - 0); q[i] = s; }
// the unit hydrograph above is for one step of rain starting at k·dt; shift by half a step is inside Tp = dt/2 + lag
let vol = 0, peak = 0, tp = 0;
for (let i = 0; i < m; i++) { vol += q[i] * dt * 60; if (q[i] > peak) { peak = q[i]; tp = i * dt; } }
return { dt, q, volume: vol, depth: cumQ, peak, tPeak: tp };
}A trapezoid hydrograph from the rational method peak.
Qp = C·i·A ÷ 360 (i in mm/h, A in hectares, Q in m³/s)
360 is the exact unit factor. The reported peak is the formula's value, while the volume is summed from the sampled ordinates, so the two are consistent only to the time step.
export function rationalHydrograph({ c, area, intensity, tc, duration, dt = 5 } = {}) {
const qp = (c * intensity * (area / 1e4)) / 360, D = Math.max(duration, tc), total = D + tc, m = Math.ceil(total / dt) + 1, q = [];
// rises linearly over tc to the peak, holds until the storm ends, falls linearly over tc
for (let i = 0; i < m; i++) {
const t = i * dt;
q.push(t <= tc ? (qp * t) / tc : t <= D ? qp : t < D + tc ? (qp * (D + tc - t)) / tc : 0);
}
let vol = 0; for (const v of q) vol += v * dt * 60;
return { dt, q, volume: vol, depth: (vol / area) * 1000, peak: qp, tPeak: tc };
}Water surface area and stored volume at each level of the pond.
slopes: area = A₀ + P₀·m·h + π·m²·h²; storage = A₀·h + P₀·m·h² ÷ 2 + π·m²·h³ ÷ 3; table: V += h ÷ 3·(A₁ + A₂ + √(A₁·A₂))
The side-slope form is evaluated directly at each level, not accumulated, and its storage is the exact integral of its area. A table of areas uses the frustum formula between rows and a running sum. Levels are bottom + depth × k ÷ n.
export function stageStorage(def) {
if (def.kind === 'slopes') {
const out = [], { bottom, area0, perimeter0, side, top, step } = def, m = side || 0;
const n = Math.max(1, Math.ceil((top - bottom) / (step || (top - bottom) / 20) - 1e-9));
for (let k = 0; k <= n; k++) {
const h = Math.min(top - bottom, (k * (top - bottom)) / n);
out.push({ z: bottom + h, area: area0 + perimeter0 * m * h + Math.PI * m * m * h * h, storage: area0 * h + (perimeter0 * m * h * h) / 2 + (Math.PI * m * m * h ** 3) / 3 });
}
return out;
}
const rows = [...(def.rows || [])].sort((a, b) => a.z - b.z), out = [];
let s = 0;
for (let k = 0; k < rows.length; k++) {
if (k) { // frustum of a pyramid between two areas: h/3·(A1 + A2 + √(A1·A2)), the standard contour-area volume
const a1 = rows[k - 1].area, a2 = rows[k].area, h = rows[k].z - rows[k - 1].z;
s += (h / 3) * (a1 + a2 + Math.sqrt(a1 * a2));
}
out.push({ z: rows[k].z, area: rows[k].area, storage: s });
}
return out;
}Area and volume below each level inside a polygon on the ground surface.
storage(z) = Σ cells below z of (z − ground) × cell area
The ground is sampled at cell centers; the default cell gives about 40,000 cells over the polygon's bounding box, so area and volume are known to the cell size. Levels are a running sum of the step, compared with the top using an absolute 1e-9.
export function groundStorage(ground, polygon, { bottom, top, step, cell } = {}) {
let x0 = Infinity, y0 = Infinity, x1 = Infinity * -1, y1 = -Infinity;
for (const p of polygon) { x0 = Math.min(x0, p.x); y0 = Math.min(y0, p.y); x1 = Math.max(x1, p.x); y1 = Math.max(y1, p.y); }
const c = cell > 0 ? cell : Math.sqrt(((x1 - x0) * (y1 - y0)) / 40000), A = c * c;
const inside = (x, y) => { let k = false; for (let i = 0, j = polygon.length - 1; i < polygon.length; j = i++) { const p = polygon[i], q = polygon[j]; if ((p.y > y) !== (q.y > y) && x < ((q.x - p.x) * (y - p.y)) / (q.y - p.y) + p.x) k = !k; } return k; };
const zs = [];
for (let y = y0 + c / 2; y < y1; y += c) for (let x = x0 + c / 2; x < x1; x += c) if (inside(x, y)) { const z = ground.elevationAt(x, y); if (z !== null && z !== undefined && Number.isFinite(z)) zs.push(z); }
if (!zs.length) return [];
const lo = bottom ?? Math.min(...zs), hi = top ?? Math.max(...zs), st = step || (hi - lo) / 20, out = [];
for (let z = lo; z <= hi + 1e-9; z += st) {
let area = 0, vol = 0;
for (const v of zs) if (v < z) { area += A; vol += (z - v) * A; }
out.push({ z, area, storage: vol });
}
return out;
}Discharge of an orifice, pipe, weir, spillway, V-notch or riser at a water level.
orifice: Cd·A·√(2g·(h − D/2)); weir: Cw·(L − 0.2·n·H)·H^1.5; V-notch: (8 ÷ 15)·Cd·√(2g)·tan(θ/2)·H^2.5
g is 9.80665. Below a head of one diameter the orifice uses a 1.5-power rise that meets the submerged formula at head = D; that blend is the code's own, stated in its comment. The contracted weir length is clamped at zero. The riser takes the smaller of weir and orifice flow.
function deviceFlow(d, z) {
switch (d.kind) {
case 'orifice': case 'pipe': {
const D = d.diameter ?? Math.sqrt((4 * d.area) / Math.PI), A = d.area ?? (Math.PI * D * D) / 4, head = z - d.invert, cd = d.cd ?? 0.6;
if (head <= 0) return 0;
const full = cd * A * Math.sqrt(2 * G * (D / 2)); // flow when the water just covers the opening (head to its centre = D/2)
// below that the opening is partly covered: a weir-like rise, exact at both ends (zero at the invert, `full` at head = D)
return head < D ? full * (head / D) ** 1.5 : cd * A * Math.sqrt(2 * G * (head - D / 2));
}
case 'weir': case 'spillway': {
const H = z - d.crest; if (H <= 0) return 0;
const broad = d.kind === 'spillway' || d.broad, cw = d.cw ?? (broad ? 1.66 : 1.84);
return cw * (broad ? d.length : Math.max(0, d.length - 0.2 * (d.contractions ?? 2) * H)) * H ** 1.5;
}
case 'vnotch': {
const H = z - d.invert; if (H <= 0) return 0;
return (8 / 15) * (d.cd ?? 0.58) * Math.sqrt(2 * G) * Math.tan(((d.angle ?? 90) * Math.PI) / 360) * H ** 2.5;
}
case 'riser': { // a box or pipe riser: weir over the top rim until it submerges, then an orifice into the barrel (the smaller governs)
const H = z - d.crest; if (H <= 0) return 0;
const per = d.perimeter ?? Math.PI * (d.diameter ?? 1), qw = (d.cw ?? 1.66) * per * H ** 1.5;
const A = d.barrelArea ?? (Math.PI * (d.barrel ?? d.diameter ?? 0.5) ** 2) / 4, qo = (d.cd ?? 0.6) * A * Math.sqrt(2 * G * Math.max(1e-9, z - (d.barrelInvert ?? d.crest - 1)));
return Math.min(qw, qo);
}
default: return 0;
}
}A value read from a table between two rows, used for storage, level and outflow.
Linear between rows found by bisection on the row index. Below the first row the first value is returned; above the last the last segment is extended in a straight line. A zero row spacing divides by 1 instead.
function interp(xs, ys, x) {
const n = xs.length;
if (x <= xs[0]) return ys[0]; if (x >= xs[n - 1]) return ys[n - 1] + ((x - xs[n - 1]) * (ys[n - 1] - ys[n - 2])) / (xs[n - 1] - xs[n - 2] || 1);
let lo = 0, hi = n - 1; while (hi - lo > 1) { const m = (lo + hi) >> 1; if (xs[m] <= x) lo = m; else hi = m; }
return ys[lo] + ((ys[hi] - ys[lo]) * (x - xs[lo])) / (xs[hi] - xs[lo] || 1);
}Storage and outflow at each time step from the inflow hydrograph.
2S₂ ÷ Δt + O₂ = I₁ + I₂ + (2S₁ ÷ Δt − O₁)
Each step reads the outflow from the storage-indication curve by linear interpolation, so no iteration is needed. The right-hand side is carried from step to step, so rounding accumulates over the storm. Outflow and storage are clamped at zero.
/** Storage–outflow table routing. s[] ascending storage (m³), o[] outflow (m³/s). Returns storage and outflow per step. */
export function routeTable(inflow, dt, { s, o, s0 = 0 } = {}) {
const T = dt * 60, ind = s.map((v, k) => (2 * v) / T + o[k]); // storage-indication curve
const outS = [s0], outO = [interp(s, o, s0)];
let ind1 = (2 * s0) / T - outO[0];
for (let i = 1; i < inflow.length; i++) {
const rhs = inflow[i - 1] + inflow[i] + ind1; // 2S₂/Δt + O₂
const O2 = Math.max(0, interp(ind, o, rhs)), S2 = Math.max(0, ((rhs - O2) * T) / 2);
outS.push(S2); outO.push(O2); ind1 = rhs - 2 * O2; // 2S₂/Δt − O₂
}
return { storage: outS, q: outO };
}Peak inflow and outflow, the highest level, the greatest storage, attenuation, draw-down time and the volumes in and out.
attenuation = 1 − peak out ÷ peak in; volume = Σ (q₁ + q₂) ÷ 2 × Δt
Peaks are the largest sampled ordinates, so they are known to the time step. Levels come back from storage by linear interpolation in the stage table. Volumes use the trapezoid rule, unlike the hydrograph's own volume, which uses rectangles.
export function routePond(inflow, dt, storage, devices, { z0 = null, tailEmpty = 1e-4 } = {}) {
const zs = storage.map((r) => r.z), S = storage.map((r) => r.storage), O = zs.map((z) => outletFlow(devices, z));
const zStart = z0 ?? zs[0], s0 = interp(zs, S, zStart);
const r = routeTable(inflow, dt, { s: S, o: O, s0 });
const z = r.storage.map((v) => interp(S, zs, v));
let peakOut = 0, tPeak = 0, peakIn = 0, zMax = -Infinity, sMax = 0;
for (let i = 0; i < r.q.length; i++) {
if (r.q[i] > peakOut) { peakOut = r.q[i]; tPeak = i * dt; }
if (inflow[i] > peakIn) peakIn = inflow[i];
if (z[i] > zMax) zMax = z[i]; if (r.storage[i] > sMax) sMax = r.storage[i];
}
// draw-down: the time from the highest level until the storage is back to the starting storage (within a thousandth of the maximum)
let drawdown = null;
const iMax = r.storage.indexOf(sMax);
for (let i = iMax; i < r.storage.length; i++) if (r.storage[i] <= s0 + Math.max(1e-6, sMax * 1e-3)) { drawdown = (i - iMax) * dt; break; }
return { dt, z, q: r.q, storage: r.storage, peakOut, peakIn, tPeak, zMax, sMax, attenuation: peakIn > 0 ? 1 - peakOut / peakIn : 0, drawdown, s0,
volumeIn: inflow.reduce((a, v, i) => a + (i ? ((inflow[i - 1] + v) / 2) * dt * 60 : 0), 0), volumeOut: r.q.reduce((a, v, i) => a + (i ? ((r.q[i - 1] + v) / 2) * dt * 60 : 0), 0) };
}The orifice diameter at which the peak outflow just reaches the allowed flow.
Bisection on the diameter between 0.02 m and 2 m, 50 halvings, each one a full routing. The bracket ends near 2e-15 m wide, far below the time-step accuracy of the peak it is matching. The lower end of the bracket is returned, so the target is not exceeded.
export function sizeOrifice(inflow, dt, storage, { target, z0 = null, base = null, cd = 0.6, lo = 0.02, hi = 2, devices = [] } = {}) {
const inv = base ?? storage[0].z, run = (dia) => routePond(inflow, dt, storage, [...devices, { kind: 'orifice', invert: inv, diameter: dia, cd }], { z0 });
if (run(lo).peakOut > target) return { diameter: lo, ...pick(run(lo)), met: false, reason: 'too-small' };
if (run(hi).peakOut < target) return { diameter: hi, ...pick(run(hi)), met: true, reason: 'larger-ok' };
let a = lo, b = hi;
for (let k = 0; k < 50; k++) { const m = (a + b) / 2; if (run(m).peakOut > target) b = m; else a = m; }
return { diameter: a, ...pick(run(a)), met: true };
function pick(r) { return { peakOut: r.peakOut, zMax: r.zMax, sMax: r.sMax }; }
}Overtopping, peak above the target, no attenuation, slow draw-down and a dry pond.
Overtopping uses an absolute 1e-9 m. The target comparison allows 0.01 % over. No attenuation means less than 5 %. Draw-down is compared with 72 hours unless another limit is given.
export function checkPond(res, { top = null, target = null, drawdownMax = 72 * 60, minAttenuation = 0.05 } = {}) {
const out = [];
if (top !== null && res.zMax > top + 1e-9) out.push({ code: 'overtop', zMax: res.zMax, top, by: res.zMax - top });
if (target !== null && res.peakOut > target * 1.0001) out.push({ code: 'over-target', peakOut: res.peakOut, target });
if (res.peakIn > 0 && res.attenuation < minAttenuation) out.push({ code: 'no-attenuation', attenuation: res.attenuation });
if (res.drawdown === null) out.push({ code: 'draindown', drawdown: null, max: drawdownMax });
else if (res.drawdown > drawdownMax) out.push({ code: 'draindown', drawdown: res.drawdown, max: drawdownMax });
if (res.sMax <= 1e-9) out.push({ code: 'dry' });
return out;
}Converts a saved pond to SI, runs it, and converts the results back for display.
1 ft = 0.3048 m; 1 ft² = 0.09290304 m²; 1 acre = 4046.8564224 m²; 1 in = 25.4 mm; 1 ft³/s = 0.028316846592 m³/s
Every factor is the exact international-foot value written out in full. units.js in the same build defines the project foot as the US survey foot (1200/3937 m), so lengths here are converted with a factor that differs from it by 2 parts per million. The time step is the storm duration ÷ 60, rounded to a whole minute between 1 and 10.
/** Run a saved pond: convert to SI, build the inflow, route it, check it; results come back in SI and in the project's units. `hydro` = the project's IDF settings (idf a, b, c) for the IDF storm. */
export function runPond(pond, units = 'ft', { hydro = null } = {}) {
const p = cleanPond(pond, units); if (!p) return null;
const m = units === 'm', L = m ? 1 : 0.3048, A = m ? 1 : 0.09290304, areaM2 = m ? p.area : p.area * 4046.8564224, depthMM = (v) => (m ? v : v * 25.4), Q = m ? 1 : 0.028316846592, V = L * L * L;
const dt = Math.max(1, Math.min(10, Math.round(p.storm.duration / 60)));
let rain;
if (p.storm.mode === 'idf' && hydro?.idf) {
const f = (tc) => hydro.idf.a / Math.pow(tc + hydro.idf.b, hydro.idf.c);
rain = alternatingBlock(f, p.storm.duration, dt).map((r) => ({ t: r.t, i: depthMM(r.i) }));
} else rain = (depthMM(p.storm.depth) / p.storm.duration) * 60;
const hy = scsHydrograph({ area: areaM2, cn: p.cn, tc: p.tc, rain, duration: p.storm.duration, dt });
const s = p.storage, rows = stageStorage({ kind: 'slopes', bottom: s.bottom * L, area0: s.area0 * A, perimeter0: s.perimeter0 * L, side: s.side, top: s.top * L, step: Math.max(0.02, (s.top - s.bottom) * L / 30) });
const devs = p.devices.map((d) => (d.kind === 'orifice' ? { ...d, invert: d.invert * L, diameter: d.diameter * L } : d.kind === 'vnotch' ? { ...d, invert: d.invert * L } : { ...d, crest: d.crest * L, length: d.length * L }));
const r = routePond(hy.q, hy.dt, rows, devs, { z0: (p.start ?? s.bottom) * L });
const checks = checkPond(r, { top: s.top * L, target: p.target ? p.target * Q : null, drawdownMax: p.drawdownMax * 60 });
return { pond: p, hydro: hy, rows, devices: devs, route: r, checks, conv: { len: L, area: A, vol: V, flow: Q }, runoffVolume: hy.volume,
show: { peakIn: hy.peak / Q, peakOut: r.peakOut / Q, zMax: r.zMax / L, storageMax: r.sMax / V, capacity: rows[rows.length - 1].storage / V, runoff: hy.volume / V } };
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.