River water surface profiles by the standard step method
Works out the water surface, energy level, critical depth, velocities, Froude number and flow split at each river cross section for steady flow.
How the precision is kept
Section properties are closed form in 64-bit. Each water surface is found with Brent's method to 1e-6 of a length unit, and critical depth by a 200-step scan refined by golden-section search to 1e-7. In feet the Manning factor is 1.486 and gravity 32.174, both rounded; in meters gravity is 9.80665. The report and CSV round to 2, 3 or 6 decimals.
Open item. Use 1 ÷ 0.3048^(1/3) for the Manning factor and 9.80665 ÷ 0.3048 for gravity in feet (UNIT_CONST in river.js) instead of the rounded 1.486 and 32.174.
Calculations under this policy
15 calculations in 2 source files. 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.
Elevations at inserted offsets are linear interpolations. Offsets are compared with an absolute 1e-12, and a vertical wall is assigned to its side with a 1e-9 nudge.
river-solve.js
lines 53 to 101
CivilCraft 0.8.1
export function geometry(sec) {
const zones = nZonesFor(sec);
const [bL, bR] = sec.banks;
// element boundaries: zone edges + banks
const cuts = new Set([bL, bR]);
for (const z of zones) { cuts.add(z.o0); cuts.add(z.o1); }
const ineff = ineffectiveOf(sec); // 0.8 SU.957: ineffective areas and levees
for (const r of ineff) { cuts.add(r.o0); cuts.add(r.o1); }
const lo = sec.pts[0].o, hi = sec.pts[sec.pts.length - 1].o;
const edges = [...cuts].filter((o) => o >= lo && o <= hi).sort((a, b) => a - b);
// points with the element boundaries inserted
const pts = [];
const zAtO = (o) => {
for (let i = 1; i < sec.pts.length; i++) {
const p = sec.pts[i - 1], q = sec.pts[i];
if (o >= p.o - 1e-12 && o <= q.o + 1e-12) return q.o > p.o ? p.z + ((q.z - p.z) * (o - p.o)) / (q.o - p.o) : q.z;
}
return null;
};
let e = 0;
for (let i = 0; i < sec.pts.length; i++) {
const p = sec.pts[i];
while (e < edges.length && edges[e] < p.o - 1e-12) {
if (i > 0 && edges[e] > sec.pts[i - 1].o + 1e-12) { const z = zAtO(edges[e]); if (z !== null) pts.push({ o: edges[e], z }); }
e++;
}
while (e < edges.length && Math.abs(edges[e] - p.o) <= 1e-12) e++;
pts.push({ o: p.o, z: p.z });
}
// elements: [o0, o1] with region and n
const elems = [];
for (let i = 1; i < edges.length; i++) {
const o0 = edges[i - 1], o1 = edges[i], mid = (o0 + o1) / 2;
if (!(o1 > o0)) continue;
const zone = zones.find((z) => mid >= z.o0 && mid <= z.o1) || zones[zones.length - 1];
elems.push({ o0, o1, n: zone.n, reg: mid < bL ? 'lob' : mid > bR ? 'rob' : 'ch' });
}
const segs = []; // each segment with its element index
for (let i = 1; i < pts.length; i++) {
const a = pts[i - 1], b = pts[i];
// a vertical wall belongs to the piece whose bed it bounds: going down → the piece to its right, going up → to its left
const mid = b.o > a.o ? (a.o + b.o) / 2 : b.z < a.z ? a.o + 1e-9 : a.o - 1e-9;
let k = elems.findIndex((el) => mid >= el.o0 - 1e-12 && mid <= el.o1 + 1e-12);
if (k < 0) k = mid < bL ? 0 : elems.length - 1;
segs.push({ a, b, k });
}
const zs = pts.map((p) => p.z);
return { sec, pts, elems, segs, ineff, minZ: Math.min(...zs), maxZ: Math.max(...zs), zL: pts[0].z, zR: pts[pts.length - 1].z, alphaFix: sec.alpha ?? null };
}
Area, wetted perimeter, conveyance and velocity coefficient
64-bit
The hydraulic properties of a section at a water level.
A = Σ (d_a + d_b) ÷ 2 × Δo; P = Σ hypot(Δo, Δz); K = (k ÷ n) × A × (A ÷ P)^(2/3); n_c = (Σ P_i·n_i^1.5 ÷ P)^(2/3); α = max(1, A²·Σ(K_i³ ÷ A_i²) ÷ K³)
A partly wet segment is cut where the depth is zero, f = d⁺ ÷ (d⁺ − d⁻); the two depths have opposite signs, so the divisor does not cancel. Sums run over segments in order. α is floored at 1. Water above the section's ends adds vertical wall to the wetted perimeter.
river-solve.js
lines 103 to 146
CivilCraft 0.8.1
export function hydraulics(geo, ws, units = 'm') {
const k = UNIT_CONST[units === 'ft' ? 'ft' : 'm'].k;
const el = geo.elems.map(() => ({ A: 0, P: 0, T: 0 }));
let M1 = 0;
const cut = (o) => geo.ineff.some((r) => o > r.o0 - 1e-12 && o < r.o1 + 1e-12 && (r.trigger === null || ws < r.trigger)); // inside an active ineffective area
for (const { a, b, k: i } of geo.segs) {
const da = ws - a.z, db = ws - b.z, dx = b.o - a.o, e = el[i];
if (da <= 0 && db <= 0) continue;
if (cut(b.o > a.o ? (a.o + b.o) / 2 : b.z < a.z ? a.o + 1e-9 : a.o - 1e-9)) continue;
if (da >= 0 && db >= 0) {
e.A += ((da + db) / 2) * dx; e.P += Math.hypot(dx, b.z - a.z); e.T += dx;
M1 += (dx * (da * da + da * db + db * db)) / 6;
} else {
const dp = Math.max(da, db), dn = Math.min(da, db), f = dp / (dp - dn), w = f * dx;
e.A += (dp * w) / 2; e.P += Math.hypot(w, dp); e.T += w;
M1 += (dp * dp * w) / 6;
}
}
let walls = false;
if (ws > geo.zL && el.length && !cut(geo.pts[0].o + 1e-9)) { el[0].P += ws - geo.zL; walls = true; }
if (ws > geo.zR && el.length && !cut(geo.pts[geo.pts.length - 1].o - 1e-9)) { el[el.length - 1].P += ws - geo.zR; walls = true; }
const reg = { lob: { A: 0, P: 0, T: 0, K: 0 }, ch: { A: 0, P: 0, T: 0, K: 0 }, rob: { A: 0, P: 0, T: 0, K: 0 } };
let chPn = 0;
for (let i = 0; i < el.length; i++) {
const r = reg[geo.elems[i].reg], e = el[i], n = geo.elems[i].n;
r.A += e.A; r.P += e.P; r.T += e.T;
if (geo.elems[i].reg === 'ch') chPn += e.P * Math.pow(n, 1.5);
else if (e.A > 0 && e.P > 0) r.K += (k / n) * e.A * Math.pow(e.A / e.P, 2 / 3);
}
if (reg.ch.A > 0 && reg.ch.P > 0) {
const nc = Math.pow(chPn / reg.ch.P, 2 / 3);
reg.ch.K = (k / nc) * reg.ch.A * Math.pow(reg.ch.A / reg.ch.P, 2 / 3);
reg.ch.n = nc;
}
const A = reg.lob.A + reg.ch.A + reg.rob.A, K = reg.lob.K + reg.ch.K + reg.rob.K;
let alpha = 1;
if (geo.alphaFix) alpha = geo.alphaFix;
else if (K > 0) {
let s = 0;
for (const r of [reg.lob, reg.ch, reg.rob]) if (r.A > 0) s += (r.K * r.K * r.K) / (r.A * r.A);
alpha = Math.max(1, (A * A * s) / (K * K * K));
}
return { A, P: reg.lob.P + reg.ch.P + reg.rob.P, T: reg.lob.T + reg.ch.T + reg.rob.T, K, alpha, M1, walls, reg };
}
Root finding
64-bit
Solves the energy balance for the water level between two bracketing levels.
stop when |bracket ÷ 2| ≤ 2ε|x| + tol ÷ 2, tol = 1e-6
Brent's method with an absolute tolerance of 1e-6 length units plus a term scaled to the size of the level, and a cap of 200 iterations. The water level is therefore good to about 1e-6 units; a run that reaches the cap is reported as not converged.
river-solve.js
lines 151 to 178
CivilCraft 0.8.1
/** Brent's method on [a, b] with f(a)·f(b) ≤ 0. Returns { x, ok, iter }. */
export function brent(f, a, b, tol = TOL, maxIter = 200) {
let fa = f(a), fb = f(b);
if (fa === 0) return { x: a, ok: true, iter: 0 };
if (fb === 0) return { x: b, ok: true, iter: 0 };
if (fa * fb > 0) return { x: NaN, ok: false, iter: 0 };
let c = a, fc = fa, d = b - a, e = d;
for (let it = 1; it <= maxIter; it++) {
if (fb * fc > 0) { c = a; fc = fa; d = e = b - a; }
if (Math.abs(fc) < Math.abs(fb)) { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; }
const tol1 = 2 * Number.EPSILON * Math.abs(b) + tol / 2, xm = (c - b) / 2;
if (Math.abs(xm) <= tol1 || fb === 0) return { x: b, ok: true, iter: it };
if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
const s = fb / fa;
let p, q;
if (a === c) { p = 2 * xm * s; q = 1 - s; } else {
const qq = fa / fc, r = fb / fc;
p = s * (2 * xm * qq * (qq - r) - (b - a) * (r - 1)); q = (qq - 1) * (r - 1) * (s - 1);
}
if (p > 0) q = -q; else p = -p;
if (2 * p < Math.min(3 * xm * q - Math.abs(tol1 * q), Math.abs(e * q))) { e = d; d = p / q; } else { d = xm; e = d; }
} else { d = xm; e = d; }
a = b; fa = fb;
b += Math.abs(d) > tol1 ? d : xm > 0 ? tol1 : -tol1;
fb = f(b);
}
return { x: b, ok: false, iter: maxIter };
}
Minimum search
64-bit
Refines the level of least specific energy.
Golden-section search until the bracket is under the tolerance (1e-7 when called for critical depth) or 200 steps, returning the middle of the bracket. Near a minimum the energy is flat, so the level is less sharply fixed than the energy.
river-solve.js
lines 180 to 187
CivilCraft 0.8.1
function golden(f, a, b, tol = TOL) {
const r = (Math.sqrt(5) - 1) / 2;
let c = b - r * (b - a), d = a + r * (b - a), fc = f(c), fd = f(d);
for (let i = 0; i < 200 && Math.abs(b - a) > tol; i++) {
if (fc < fd) { b = d; d = c; fd = fc; c = b - r * (b - a); fc = f(c); } else { a = c; c = d; fc = fd; d = a + r * (b - a); fd = f(d); }
}
return (a + b) / 2;
}
Specific energy
64-bit
The energy level at a section for a water level and a flow.
E = WS + α·Q² ÷ (2g·A²)
One 64-bit expression; a section with no flow area returns infinity instead of dividing by zero. g is 9.80665 in meters and 32.174 in feet (river.js).
river-solve.js
lines 191 to 195
CivilCraft 0.8.1
/** Specific energy above the datum at water surface ws. */
function energy(geo, ws, q, units) {
const h = hydraulics(geo, ws, units);
return h.A > 0 ? ws + (h.alpha * q * q) / (2 * G(units) * h.A * h.A) : Infinity;
}
Critical water surface
64-bit
The water level of least specific energy at a section.
The range from the lowest point to the top of the section is scanned in 200 steps, then the best step is refined by golden-section search to 1e-7. If the minimum lies at the top the range is doubled, up to 40 times. A second, lower minimum between scan points could be missed.
river-solve.js
lines 197 to 209
CivilCraft 0.8.1
export function criticalWs(geo, q, units = 'm') {
const E = (ws) => energy(geo, ws, q, units);
const lo = geo.minZ;
let hi = Math.max(geo.maxZ, lo + 1e-3);
for (let grow = 0; grow < 40; grow++) {
const N = 200, step = (hi - lo) / N;
let bi = 1, be = Infinity;
for (let i = 1; i <= N; i++) { const v = E(lo + i * step); if (v < be) { be = v; bi = i; } }
if (bi < N) return golden(E, lo + Math.max(bi - 1, 0) * step + 1e-12, lo + (bi + 1) * step, TOL / 10);
hi = lo + (hi - lo) * 2; // the minimum is above the section: walls carry it
}
return hi;
}
Normal depth boundary
64-bit
The water level at which the section carries the flow at a given energy slope.
K(WS) = Q ÷ √S
Solved with Brent's method to 1e-6 after doubling the upper level until the conveyance is enough (up to 60 times).
river-solve.js
lines 211 to 219
CivilCraft 0.8.1
export function normalWs(geo, q, slope, units = 'm') {
if (!(slope > 0)) return null;
const target = q / Math.sqrt(slope);
const f = (ws) => hydraulics(geo, ws, units).K - target;
let hi = Math.max(geo.maxZ, geo.minZ + 1e-3);
for (let i = 0; i < 60 && f(hi) < 0; i++) hi = geo.minZ + (hi - geo.minZ) * 2;
const r = brent(f, geo.minZ + 1e-12, hi);
return r.ok ? r.x : null;
}
Velocity head and flow split
64-bit
The velocity head, energy level and the flow in the left overbank, channel and right overbank.
h_v = α·Q² ÷ (2g·A²); Q_part = Q × K_part ÷ K
One expression each; the split is zero when the section has no conveyance.
The reach length is weighted by the flows averaged over the two sections. The velocity head difference is a subtraction of two like quantities and loses digits when they are nearly equal; it is then multiplied by a coefficient below 1.
river-solve.js
lines 231 to 241
CivilCraft 0.8.1
/** Loss from up (u, the section that owns the reach) to down (d). */
function loss(up, u, d) {
const r = up.reach || {};
const ql = (u.qLob + d.qLob) / 2, qc = (u.qCh + d.qCh) / 2, qr = (u.qRob + d.qRob) / 2, qt = ql + qc + qr;
const L = qt > 0 ? ((r.lob ?? r.ch ?? 0) * ql + (r.ch ?? 0) * qc + (r.rob ?? r.ch ?? 0) * qr) / qt : r.ch ?? 0;
const Ks = u.h.K + d.h.K;
const sf = Ks > 0 ? Math.pow((u.q + d.q) / Ks, 2) : Infinity;
const C = d.hv > u.hv ? up.cc : up.ce;
const hf = L * sf, hce = C * Math.abs(u.hv - d.hv);
return { L, sf, hf, hce, he: hf + hce };
}
Results at a section
64-bit
Depth, velocities, Froude numbers, energy slope and specific force for the results table.
Fr = V ÷ √(g·A ÷ T); S_f = (Q ÷ K)²; specific force = Q² ÷ (g·A) + A·ȳ
Each value is one 64-bit expression from the section properties; divisions are guarded by tests for a positive area and top width.
Steps upstream from the downstream boundary, solving each section's water level.
The bracket runs from critical depth up to a level doubled until the balance changes sign (up to 60 times); with no bracket the section is set to critical depth and flagged. A subcritical answer with a channel Froude number above 1 + 1e-6 is flagged.
river-solve.js
lines 264 to 297
CivilCraft 0.8.1
function runSub(secs, geos, profile, units) {
const rows = [];
const n = secs.length;
let prev = null;
for (let i = 0; i < n; i++) {
const sec = secs[i], geo = geos[i], q = flowAt(profile, sec.rs);
const crit = criticalWs(geo, q, units);
let ws, flags = [], lossInfo = null;
if (i === 0) {
ws = boundaryWs(profile.down, geo, q, units);
if (ws === null || ws < crit) { flags.push('critical'); ws = crit; }
} else {
const f = (w) => { const u = state(geo, w, q, units); return u.H - (prev.H + loss(sec, u, prev).he); };
const fc = f(crit);
if (!(fc < 0)) { ws = crit; flags.push('critical'); } else {
let hi = Math.max(crit, prev.H) + Math.max(1e-3, prev.H - geo.minZ);
for (let k = 0; k < 60 && !(f(hi) > 0); k++) hi = crit + (hi - crit) * 2;
const r = brent(f, crit, hi);
ws = r.ok ? r.x : crit;
if (!r.ok) flags.push('noconv');
}
}
const st = state(geo, ws, q, units);
if (prev) lossInfo = loss(sec, st, prev);
const row = makeRow(sec, geo, st, crit, units, lossInfo ? { length: lossInfo.L, hf: lossInfo.hf, hce: lossInfo.hce } : {});
row.flags.push(...flags);
row.regime = flags.includes('critical') ? 'crit' : 'sub';
if (row.regime === 'sub' && row.frCh > 1 + 1e-6) row.flags.push('froude');
if (prev && prev.h.K > 0) { const kr = st.h.K / prev.h.K; if (kr < 0.7 || kr > 1.4) row.flags.push('kratio'); }
rows.push(row);
prev = st;
}
return rows;
}
Supercritical run
64-bit
Steps downstream from the upstream boundary, solving each section's water level below critical depth.
The lower bracket starts at a thousandth of the critical depth above the bed (at least 1e-9) and is quartered toward the bed until the balance changes sign.
river-solve.js
lines 299 to 335
CivilCraft 0.8.1
function runSuper(secs, geos, profile, units) {
const rows = new Array(secs.length);
const n = secs.length;
let prev = null; // upstream state
for (let i = n - 1; i >= 0; i--) {
const sec = secs[i], geo = geos[i], q = flowAt(profile, sec.rs);
const crit = criticalWs(geo, q, units);
let ws, flags = [], lossInfo = null;
if (i === n - 1) {
ws = boundaryWs(profile.up, geo, q, units);
if (ws === null || ws > crit) { flags.push('critical'); ws = crit; }
} else {
const up = secs[i + 1];
const f = (w) => { const d = state(geo, w, q, units); return d.H + loss(up, prev, d).he - prev.H; };
const fc = f(crit);
if (!(fc < 0)) { ws = crit; flags.push('critical'); } else {
let lo = geo.minZ + Math.max(1e-9, (crit - geo.minZ) * 1e-3);
for (let k = 0; k < 60 && !(f(lo) > 0); k++) lo = geo.minZ + (lo - geo.minZ) / 4;
const r = brent(f, lo, crit);
ws = r.ok ? r.x : crit;
if (!r.ok) flags.push('noconv');
}
}
const st = state(geo, ws, q, units);
rows[i] = makeRow(sec, geo, st, crit, units);
rows[i].flags.push(...flags);
rows[i].regime = flags.includes('critical') ? 'crit' : 'super';
if (rows[i].regime === 'super' && rows[i].frCh < 1 - 1e-6) rows[i].flags.push('froude');
if (prev) {
lossInfo = loss(secs[i + 1], prev, st);
Object.assign(rows[i + 1], { length: lossInfo.L, hf: lossInfo.hf, hce: lossInfo.hce });
if (prev.h.K > 0) { const kr = prev.h.K / st.h.K; if (kr < 0.7 || kr > 1.4) rows[i + 1].flags.push('kratio'); }
}
prev = st;
}
return rows;
}
Mixed regime
64-bit
Keeps, at each section, the subcritical or supercritical answer with the larger specific force, and marks hydraulic jumps.
The choice is a direct comparison of two 64-bit specific forces with no tolerance.
river-solve.js
lines 337 to 346
CivilCraft 0.8.1
function runMixed(secs, geos, profile, units) {
const sub = runSub(secs, geos, profile, units), sup = runSuper(secs, geos, profile, units);
const rows = sub.map((r, i) => {
const s = sup[i];
const pickSuper = s.regime === 'super' && (r.regime === 'crit' || s.specificForce > r.specificForce);
return { ...(pickSuper ? s : r), flags: [...(pickSuper ? s : r).flags] };
});
for (let i = 1; i < rows.length; i++) if (rows[i].regime === 'super' && rows[i - 1].regime !== 'super') rows[i].flags.push('jump');
return rows;
}
The number of decimals each result is shown with in the table, the text report and the CSV.
Levels, depth, velocities and Froude number are written with 3 decimals, river station, flows, area and top width with 2, and the energy slope with 6. The CSV carries these rounded values, not the 64-bit results.
Formats each section's results for the table and the CSV.
text = v.toFixed(d)
Each value is rounded with toFixed to its column's decimals; a value that is not finite is written as empty.
river-report.js
lines 47 to 49
CivilCraft 0.8.1
export function tableRows(profile, { t = (k) => k } = {}) {
return profile.rows.map((r) => [...COLUMNS.map((c) => f(r[c.key], c.d)), r.flags.map((x) => t('sx.su957.flag.' + x + '.short')).join('; ')]);
}
In detail
Method
Standard step energy balance (the file cites Chow 1959 ch. 10 and a hydraulic reference manual by Brunner, ch. 2): average-conveyance friction slope, contraction and expansion losses, conveyance by subdivision, composite channel roughness, velocity coefficient α.
Number types
64-bit numbers.
Tolerances
Brent: 1e-6 absolute plus 2ε|x|, at most 200 iterations. Golden section: 1e-7, at most 200 iterations. Critical scan: 200 steps. Bracket searches double up to 60 times. Froude flag at 1 ± 1e-6. Conveyance ratio flag outside 0.7 to 1.4. Geometry comparisons 1e-12; vertical walls nudged 1e-9.
Limits
Water levels are solved to 1e-6 units, not to 64-bit accuracy. The location of a hydraulic jump between two sections is not found. The constants k and g come from UNIT_CONST in river.js: 1.486 and 32.174 for feet.
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.