Works out the headwater above a culvert, the water surface above a bridge, and the flows and water levels where river reaches join, for the steady water surface profile.
How the precision is kept
All arithmetic is in 64-bit numbers. The results are set by the method, not by the number type: weir and orifice formulas with typed coefficients, a blend between them, and root finding. Culvert headwater is found by bisection with a fixed number of halvings; bridge water surfaces are found with the Brent root finder of river-solve.js at its default tolerance of 1e-6 length units. Gravity and the Manning constant come from the rounded table in river.js.
Open item. Pass the Brent root finder a tolerance scaled to the flow depth for the bridge water surfaces instead of its fixed default of 1e-6 length units.
Calculations under this policy
8 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.
The flow area, water surface width and full-barrel hydraulic radius of one culvert barrel at a given depth.
circle: θ = 2·acos((r − y) ÷ r); A = r²(θ − sin θ) ÷ 2; T = 2√(y(D − y)); R = D ÷ 4. box: A = B·y; R = B·D ÷ (2(B + D))
64-bit closed forms. The acos argument stays inside (−1, 1) because depths at or below zero and at or above the rise return before it is evaluated. For a very shallow depth θ − sin θ subtracts two nearly equal numbers and loses digits; the square root for the width is guarded with max(0, …).
river-structures.js
lines 54 to 67
CivilCraft 0.8.1
/** Flow area, top width and the full-barrel hydraulic radius of one barrel at depth y. */
function barrel(c) {
const D = c.rise;
if (c.shape === 'circle') {
const r = D / 2;
return {
D, B: D, full: Math.PI * r * r, R: D / 4,
area: (y) => { if (y <= 0) return 0; if (y >= D) return Math.PI * r * r; const th = 2 * Math.acos((r - y) / r); return (r * r * (th - Math.sin(th))) / 2; },
top: (y) => { if (y <= 0 || y >= D) return 0; return 2 * Math.sqrt(Math.max(0, y * (D - y))); },
};
}
const B = c.span;
return { D, B, full: B * D, R: (B * D) / (2 * (B + D)), area: (y) => B * clamp(y, 0, D), top: (y) => (y > 0 && y < D ? B : 0) };
}
Critical depth in a barrel
64-bit
The depth at which the flow in one barrel is critical, capped at the rise.
Q²·T(y) − g·A(y)³ = 0
Solved with the Brent root finder between 1e-9 × rise and 0.999 × rise, with tolerance 1e-9 × rise, so the tolerance scales with the barrel. If the function is still positive at 0.999 × rise, or the root finder does not converge, the rise is returned.
river-structures.js
lines 69 to 76
CivilCraft 0.8.1
/** Critical depth in one barrel for flow qb (Q²T = gA³; the depth is capped at the rise). */
function barrelCritical(br, qb, g) {
const f = (y) => qb * qb * br.top(y) - g * Math.pow(br.area(y), 3);
const hi = br.D * 0.999;
if (f(hi) > 0) return br.D; // still supercritical at the crown: the barrel runs full
const r = brent(f, 1e-9 * br.D, hi, 1e-9 * br.D);
return r.ok ? r.x : br.D;
}
Headwater for one barrel's flow
64-bit
The inlet-control and outlet-control headwater depths for a flow in one barrel, and which one governs.
64-bit expressions with Math.pow for the 2/3 and 4/3 powers. The divisions by length, span, area and coefficients are not guarded here; the model reader in river.js only accepts positive values for them. g and k come from the rounded unit table (32.174 and 1.486 in feet).
The headwater at which the flow through the barrels plus the flow over the road equals the total flow.
n·Qbarrel(HW) + Cr·Lr·(HW − crest)^1.5 = Q
Two nested bisections with fixed counts and no tolerance test: 60 halvings of the headwater bracket (which starts at the rise and doubles until it holds the flow), and inside each step 70 halvings to invert the barrel headwater. Sixty halvings shrink the bracket by 2⁶⁰, below the spacing of 64-bit numbers at that size. The barrel count is made a whole number with a bitwise OR; it is a count, not a measured value.
river-structures.js
lines 95 to 116
CivilCraft 0.8.1
export function culvertHeadwater(c, q, tw, units = 'm') {
const g = G(units), br = barrel(c), nb = Math.max(1, c.barrels | 0), flags = [];
const roadQ = (hw) => { const over = c.invUp + hw - c.roadCrest; return Number.isFinite(c.roadCrest) && over > 0 ? c.roadCw * c.roadLength * Math.pow(over, 1.5) : 0; };
// barrel flow for a headwater: invert hwOfBarrelFlow by bisection (it grows with the flow)
const qbOf = (hw) => {
let lo = 0, hi = Math.max(1e-6, q / nb);
for (let i = 0; i < 80 && hwOfBarrelFlow(c, br, hi, tw, units).hw < hw; i++) hi *= 2;
for (let i = 0; i < 70; i++) { const m = (lo + hi) / 2; if (hwOfBarrelFlow(c, br, m, tw, units).hw < hw) lo = m; else hi = m; }
return (lo + hi) / 2;
};
const total = (hw) => nb * qbOf(hw) + roadQ(hw);
let lo = 0, hi = br.D;
for (let i = 0; i < 60 && total(hi) < q; i++) hi *= 2;
for (let i = 0; i < 60; i++) { const m = (lo + hi) / 2; if (total(m) < q) lo = m; else hi = m; }
const hw = (lo + hi) / 2, qRoad = roadQ(hw), qCulvert = Math.max(0, q - qRoad);
const at = hwOfBarrelFlow(c, br, qCulvert / nb, tw, units);
if (qRoad > 1e-9) flags.push('overtop');
if (hw > 1.5 * br.D) flags.push('submerged');
if (at.v > (units === 'ft' ? 15 : 4.5)) flags.push('fast'); // a typical scour limit: a hint, not a rule
void g;
return { hw, control: at.control, qCulvert, qRoad, hwIn: at.hwIn, hwOut: at.hwOut, dc: at.dc, v: at.v, flags };
}
Bridge opening section
64-bit
The downstream cross-section clipped to the abutments, with each pier as a strip that carries no flow.
z(o) = z₀ + (o − o₀) ÷ (o₁ − o₀) × (z₁ − z₀)
Ground levels at the two abutments are linear interpolations in 64-bit. Section points are kept when they lie more than 1e-9 inside the abutments; both that and the 1e-12 used in the lookup are absolute offsets, not scaled to the section width.
river-structures.js
lines 121 to 139
CivilCraft 0.8.1
/** The downstream section clipped to the abutments, piers as areas that carry no flow. */
export function openingSection(secD, b) {
const lo = Math.max(secD.pts[0].o, b.abutL), hi = Math.min(secD.pts[secD.pts.length - 1].o, b.abutR);
const zAt = (o) => {
const p = secD.pts;
for (let i = 1; i < p.length; i++) {
if (o <= p[i].o + 1e-12 && p[i].o > p[i - 1].o) { const f = (o - p[i - 1].o) / (p[i].o - p[i - 1].o); return p[i - 1].z + f * (p[i].z - p[i - 1].z); }
}
return p[p.length - 1].z;
};
const pts = [{ o: lo, z: zAt(lo) }];
for (const p of secD.pts) if (p.o > lo + 1e-9 && p.o < hi - 1e-9) pts.push({ o: p.o, z: p.z });
pts.push({ o: hi, z: zAt(hi) });
const bl = Math.max(lo, Math.min(hi, secD.banks[0])), br = Math.max(lo, Math.min(hi, secD.banks[1]));
const banks = br > bl ? [bl, br] : [lo, hi];
const ineff = (b.piers || []).map((p) => ({ o0: Math.max(lo, p.o - p.w / 2), o1: Math.min(hi, p.o + p.w / 2), trigger: null })).filter((r) => r.o1 > r.o0);
const n = nZonesFor(secD).map((z) => ({ o: Math.max(lo, z.o0), n: z.n })).filter((z) => z.o < hi);
return { ...secD, id: secD.id + '-bridge', pts, banks, n: n.length ? n : [{ o: lo, n: 0.035 }], ineff, levees: null, alpha: null };
}
Water surface above a bridge
64-bit
The water surface in the opening and at the section above the bridge, in low flow, pressure flow or flow over the deck.
Each water surface is a Brent root of a 64-bit energy balance, at the root finder's default tolerance of 1e-6 in model length units (absolute). Brackets are widened by doubling up to 60 times; a failed solve falls back to critical depth and adds the flag noconv. Reach lengths are floored at 1e-6, and the velocity head is infinite when the area is zero.
river-structures.js
lines 146 to 194
CivilCraft 0.8.1
export function bridgeUpstream(b, secD, secU, q, wsD, units = 'm') {
const g = G(units), flags = [];
const geoD = geometry(secD), geoU = geometry(secU), opening = openingSection(secD, b), geoB = geometry(opening);
const dn = stateAt(geoD, wsD, q, units), eD = dn.ws + dn.vh;
const critB = criticalWs(geoB, q, units), critU = criticalWs(geoU, q, units);
const LB = Math.max(1e-6, b.rsB - b.rsD), LU = Math.max(1e-6, b.rsU - b.rsB);
// opening, from the tailwater (expansion on the way out)
const fB = (w) => { const s = stateAt(geoB, w, q, units); const sf = Math.pow(q / (0.5 * (s.h.K + dn.h.K)), 2); return s.ws + s.vh - (eD + LB * sf + b.ce * Math.max(0, s.vh - dn.vh)); };
let wsB, choked = false;
if (!(fB(critB) < 0)) { wsB = critB; choked = true; flags.push('choked'); } else {
let hi = Math.max(critB, eD) + Math.max(1e-3, eD - geoB.minZ);
for (let k = 0; k < 60 && !(fB(hi) > 0); k++) hi = critB + (hi - critB) * 2;
const r = brent(fB, critB, hi);
wsB = r.ok ? r.x : critB;
if (!r.ok) flags.push('noconv');
}
const sB = stateAt(geoB, wsB, q, units), eB = sB.ws + sB.vh;
// upstream section (contraction on the way in)
const fU = (w) => { const s = stateAt(geoU, w, q, units); const sf = Math.pow(q / (0.5 * (s.h.K + sB.h.K)), 2); return s.ws + s.vh - (eB + LU * sf + b.cc * Math.max(0, sB.vh - s.vh)); };
let hi = Math.max(critU, eB) + Math.max(1e-3, eB - geoU.minZ);
for (let k = 0; k < 60 && !(fU(hi) > 0); k++) hi = critU + (hi - critU) * 2;
const rU = fU(critU) < 0 ? brent(fU, critU, hi) : { ok: true, x: critU };
let wsU = rU.ok ? rU.x : critU, mode = 'low', qPress = q, qWeir = 0;
const uS = (w) => stateAt(geoU, w, q, units);
if (wsU > b.lowChord) {
// pressure (and weir) flow
const An = hydraulics(geoB, b.lowChord, units).A, zb = geoB.minZ, Z = b.lowChord - zb;
const Lw = Number.isFinite(b.weirLength) ? b.weirLength : secU.pts[secU.pts.length - 1].o - secU.pts[0].o;
const flowAt2 = (w) => {
const s = uS(w), egl = s.ws + s.vh;
const head = wsD >= b.lowChord ? egl - wsD : egl - (zb + Z / 2);
const qp = head > 0 ? b.cd * An * Math.sqrt(2 * g * head) : 0;
const qw = egl > b.deckTop ? b.cw * Lw * Math.pow(egl - b.deckTop, 1.5) : 0;
return { qp, qw };
};
const f = (w) => { const r = flowAt2(w); return r.qp + r.qw - q; };
let top = Math.max(wsU, b.lowChord) + 1;
for (let k = 0; k < 60 && !(f(top) > 0); k++) top = b.lowChord + (top - b.lowChord) * 2;
const r = f(b.lowChord) < 0 ? brent(f, b.lowChord, top) : { ok: true, x: b.lowChord };
if (r.ok && r.x > wsU) wsU = r.x;
const fl = flowAt2(wsU);
qPress = fl.qp; qWeir = fl.qw;
mode = qWeir > 1e-9 ? 'weir' : 'pressure';
flags.push(mode === 'weir' ? 'overtop' : 'pressure');
if (!r.ok) flags.push('noconv');
}
const sU = uS(wsU);
return { wsU, wsB, egl: sU.ws + sU.vh, mode, qPress, qWeir, choked, flags, openingArea: hydraulics(geoB, Math.min(wsB, b.lowChord), units).A };
}
Stepping a profile across a structure
64-bit
Takes the water surface below a structure, works out the level above it and uses that as the start of the next reach upstream.
The tailwater depth is one 64-bit subtraction, floored at zero. The upstream level is passed on unrounded as the downstream boundary of the next reach.
river-structures.js
lines 234 to 247
CivilCraft 0.8.1
const q = flowAt(profile, s.rs), wsD = dRow.wse;
let wsU, info;
if (s.kind === 'culvert') {
const tw = Math.max(0, wsD - s.invDown), r = culvertHeadwater(s, q, tw, units);
wsU = s.invUp + r.hw;
info = { id: s.id, kind: 'culvert', rs: s.rs, q, tw: wsD, hwWse: wsU, control: r.control, qCulvert: r.qCulvert, qRoad: r.qRoad, v: r.v, dc: r.dc, flags: r.flags };
} else {
const r = bridgeUpstream({ ...s, rsD: secD.rs, rsU: secU.rs, rsB: s.rs }, secD, secU, q, wsD, units);
wsU = r.wsU;
info = { id: s.id, kind: 'bridge', rs: s.rs, q, tw: wsD, hwWse: wsU, control: r.mode, qCulvert: r.qPress, qRoad: r.qWeir, v: 0, dc: r.wsB, flags: r.flags };
}
out.push(info);
for (const f of info.flags) flagsTotal[f] = (flagsTotal[f] || 0) + 1;
down = { type: 'ws', wse: wsU };
Flows and water levels at a junction
64-bit
Adds the flows of the reaches that join and carries the water surface at the head of the downstream reach up to them.
Q = Σ Q of the joining reaches; WS of each joining reach's downstream end = WS at the head of the reach below
Flows are added in list order in 64-bit; the water surface is copied, not recomputed. No angle or momentum term is applied.
river-structures.js
lines 278 to 295
CivilCraft 0.8.1
const qOf = (id, pi) => { // the flow of profile index pi at the head of reach id (sum of the ones joining it)
const r = reaches.get(id), k = kids.get(id) || [];
return k.length ? k.reduce((s, c) => s + qOf(c, pi), 0) : r.model.profiles[pi]?.q ?? 0;
};
const run = (id, downFor) => { // downFor(pi) → boundary of this reach's profile pi
const r = reaches.get(id);
const m = { ...r.model, profiles: r.model.profiles.map((p, pi) => ({ ...p, q: qOf(id, pi), changes: (kids.get(id) || []).length ? [] : p.changes, down: downFor(pi) })) };
const hasStruct = (m.structures || []).length > 0;
const res = hasStruct ? solveStructures(m) : { units: m.units, issues: [], profiles: m.profiles.map((p) => solveProfile({ ...m, regime: m.regime || 'sub' }, p)) };
out[id] = res; order.push(id);
const k = kids.get(id) || [];
if (k.length) {
const heads = res.profiles.map((p) => p.rows[p.rows.length - 1]);
const j = { into: id, ups: k, profiles: heads.map((h, pi) => ({ name: m.profiles[pi].name, q: h.q, wse: h.wse, parts: k.map((c) => qOf(c, pi)) })) };
junctions.push(j);
for (const c of k) run(c, (pi) => ({ type: 'ws', wse: heads[pi].wse }));
}
};
In detail
Method
Culvert: inlet control (weir Q = Cw·B·HW^1.5, orifice Q = Cd·A·√(2g(HW − D/2)), blended between HW/D = 1 and 1.5) against outlet control (energy with entrance and friction losses); the larger governs. Bridge: energy balance through the clipped opening, then pressure and weir flow. Junction: flows add and the downstream water surface is carried up. The file header says the published design-chart polynomials are not used.
Number types
64-bit numbers throughout.
Tolerances
Brent tolerance 1e-9 × barrel rise for the critical depth in a barrel; Brent default 1e-6 (absolute, in model length units) for bridge water surfaces; bisection with 60 outer and 70 inner halvings and no tolerance test for culvert headwater; 1e-12 and 1e-9 (absolute, offset units) when clipping a section to the abutments; reach lengths floored at 1e-6; flows below 1e-9 over the road or deck count as none.
Limits
Subcritical profiles only. Default coefficients are typed round numbers (1.7 or 3.1 for weirs, 0.6 and 0.8 for orifices). Junctions use equal water surfaces with no angle or momentum term. Gravity in feet is 32.174 and the Manning constant is 1.486, both rounded.
Approach: 64-bit, stated limit
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.