Resources Numerical accuracy Calculation policy

Activated sludge plant steady state (ASM1)

Works out the steady concentrations in a chain of treatment tanks with recycles and a clarifier, the effluent and waste sludge, the sludge age and the mass balances.

How the precision is kept

64-bit arithmetic. The steady state is reached by an implicit iteration with a finite-difference Jacobian and a dense linear solve, stopped at a relative rate tolerance. COD and nitrogen balances are then computed separately from the final state and reported with their residuals.

Open item. Handle the null that solveLinear returns for a singular matrix; the catch around the solve only handles a thrown error.

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.

process/env/water.js CadCraft 0.8.2 64-bit 8 calculations Reviewed Nov 5, 2026

Tests tools/process-test/test-env.mjs

Stoichiometric matrix

64-bit

The 8 by 13 table of how much of each component each process makes or uses.

Each coefficient is a short 64-bit expression of the parameters, rebuilt on every call. The oxygen equivalents 2.86 and 4.57 and the alkalinity divisor 14 are written as decimal constants.

process/env/water.js
lines 37 to 51
CadCraft 0.8.2

/** The stoichiometric matrix as 8 rows (processes) of 13 coefficients, for the parameter set `par`. */
export function STOICH(par = PAR) {
  const { YA, YH, fP, iXB, iXP } = par, z = () => Object.fromEntries(COMP.map((k) => [k, 0]));
  const rows = [];
  let r;
  r = z(); Object.assign(r, { SS: -1 / YH, XBH: 1, SO: -(1 - YH) / YH, SNH: -iXB, SALK: -iXB / 14 }); rows.push(r);                                     // 1 aerobic growth, heterotrophs
  r = z(); Object.assign(r, { SS: -1 / YH, XBH: 1, SNO: -(1 - YH) / (2.86 * YH), SNH: -iXB, SALK: (1 - YH) / (14 * 2.86 * YH) - iXB / 14 }); rows.push(r);   // 2 anoxic growth
  r = z(); Object.assign(r, { XBA: 1, SNO: 1 / YA, SO: -(4.57 - YA) / YA, SNH: -iXB - 1 / YA, SALK: -iXB / 14 - 1 / (7 * YA) }); rows.push(r);          // 3 aerobic growth, autotrophs
  r = z(); Object.assign(r, { XP: fP, XS: 1 - fP, XBH: -1, XND: iXB - fP * iXP }); rows.push(r);                                                     // 4 decay, heterotrophs
  r = z(); Object.assign(r, { XP: fP, XS: 1 - fP, XBA: -1, XND: iXB - fP * iXP }); rows.push(r);                                                     // 5 decay, autotrophs
  r = z(); Object.assign(r, { SNH: 1, SND: -1, SALK: 1 / 14 }); rows.push(r);                                                                         // 6 ammonification
  r = z(); Object.assign(r, { SS: 1, XS: -1 }); rows.push(r);                                                                                         // 7 hydrolysis, organics
  r = z(); Object.assign(r, { SND: 1, XND: -1 }); rows.push(r);                                                                                       // 8 hydrolysis, organic N
  return rows;
}

Process rates

64-bit

The rate of each of the 8 biological processes at given concentrations.

r₁ = µH × SS ÷ (KS + SS) × SO ÷ (KOH + SO) × XBH (and seven more Monod products)

Negative concentrations are treated as 0 before the Monod terms. The ratio XS ÷ XBH is taken as 0 when XBH is at or below 1e-12 and the last rate as 0 when XS is at or below 1e-12; both thresholds are absolute, in g/m³.

process/env/water.js
lines 53 to 61
CadCraft 0.8.2

/** The 8 process rates (g/m3/d) at concentrations c (object by COMP name). */
export function rates(c, par = PAR) {
  const p = par, SS = pos(c.SS), SO = pos(c.SO), SNO = pos(c.SNO), SNH = pos(c.SNH), SND = pos(c.SND), XS = pos(c.XS), XBH = pos(c.XBH), XBA = pos(c.XBA), XND = pos(c.XND);
  const fS = SS / (p.KS + SS), fO = SO / (p.KOH + SO), fAn = p.KOH / (p.KOH + SO), fNO = SNO / (p.KNO + SNO);
  const r1 = p.muH * fS * fO * XBH, r2 = p.muH * fS * fAn * fNO * p.etag * XBH, r3 = p.muA * SNH / (p.KNH + SNH) * SO / (p.KOA + SO) * XBA;
  const r4 = p.bH * XBH, r5 = p.bA * XBA, r6 = p.ka * SND * XBH;
  const xr = XBH > 1e-12 ? XS / XBH : 0, r7 = p.kh * (xr / (p.KX + xr)) * (fO + p.etah * fAn * fNO) * XBH, r8 = XS > 1e-12 ? r7 * XND / XS : 0;
  return [r1, r2, r3, r4, r5, r6, r7, r8];
}

Reaction rate of each component

64-bit

The net production of each of the 13 components from the 8 process rates.

dcₖ/dt = Σⱼ νⱼₖ × ρⱼ

Eight products added in a fixed order per component. Terms of opposite sign cancel near steady state, which is the expected behavior of a balance.

process/env/water.js
lines 62 to 65
CadCraft 0.8.2

/** Reaction rate of every component (object), g/m3/d. */
export function derivs(c, par = PAR) {
  const rho = rates(c, par), M = STOICH(par), out = {}; for (const k of COMP) { let s = 0; for (let j = 0; j < 8; j++) s += M[j][k] * rho[j]; out[k] = s; } return out;
}

Clarifier split

64-bit

The effluent and underflow concentrations from the clarifier feed and the capture fraction a person gives.

eff = (1 − capture) × c; und = (Qc × c − Qe × eff) ÷ Qu

A mass balance on each particulate component; soluble components pass unchanged. The underflow is a difference of two flows of mass divided by the underflow, with no rounding.

process/env/water.js
lines 84 to 87
CadCraft 0.8.2

const clarify = (cf) => {                                   // feed concentration vector → { eff, und }
  const eff = cf.map((v, i) => (isP[i] ? (1 - plant.capture) * v : v)), und = cf.map((v, i) => (isP[i] && Qu > 0 ? (Qc * v - Qe * eff[i]) / Qu : v));
  return { eff, und };
};

Tank mass balances

64-bit

The rate of change of every component in every tank: flow in and out, reaction and aeration.

dc/dt = (Σ Q_in × c_in − Q × c) ÷ V + Σⱼ νⱼ × ρⱼ + KLa × (SO_sat − SO)

Transport, reaction and aeration are added as three 64-bit terms per component. Oxygen saturation is the constant 8 g/m³; a tank with a fixed oxygen value has its oxygen rate set to 0.

process/env/water.js
lines 88 to 103
CadCraft 0.8.2

const f = (_t, y) => {
  const c = Array.from({ length: n }, (_, i) => y.slice(i * 13, i * 13 + 13)), { und } = clarify(c[n - 1]);
  const dy = new Array(y.length).fill(0);
  for (let i = 0; i < n; i++) {
    const o = obj(c[i]), rho = rates(o, par), t = T[i], q = Qt / t.V;
    for (let k = 0; k < 13; k++) {
      const key = COMP[k]; let tr;
      if (i === 0) tr = (Qin * infl[k] + Qr * und[k] + Qint * c[n - 1][k] - Qt * c[0][k]) / t.V; else tr = q * (c[i - 1][k] - c[i][k]);
      let rx = 0; for (let j = 0; j < 8; j++) rx += M[j][key] * rho[j];
      let aer = 0; if (key === 'SO' && t.aerated) aer = (t.KLa || 0) * (SO_SAT - c[i][k]);
      dy[i * 13 + k] = tr + rx + aer;
      if (key === 'SO' && t.SOfixed != null) dy[i * 13 + k] = 0;
    }
  }
  return dy;
};

Steady-state iteration

64-bit

Advances the plant state with implicit steps until nothing changes any more.

(I ÷ h − J) × Δy = F(y); stop when max |Fᵢ| ÷ (|yᵢ| + 1) < tol

The Jacobian is a forward difference with step 1e-6 × max(1, |y|), so the steps are approximate but the stop test is on the true rates: 1e-7 by default, with a cap of 4000 iterations and a warning when the cap is reached. The catch around the linear solve handles a thrown error, while solveLinear in process/numeric.js returns null for a singular matrix. Concentrations other than alkalinity are clamped at 0 after each step.

process/env/water.js
lines 107 to 119
CadCraft 0.8.2

const N = y.length, scale = (v) => Math.abs(v) + 1, resid = (yy) => { const d = f(0, yy); let m = 0; for (let i = 0; i < N; i++) m = Math.max(m, Math.abs(d[i]) / scale(yy[i])); return m; };
const maxIter = 4000, tol = plant.tol ?? 1e-7; let day = 0, conv = false, h = 1e-5, r0 = resid(y), it = 0;
while (it < maxIter) {
  const F = f(0, y), Jm = Array.from({ length: N }, () => new Array(N).fill(0));
  for (let j = 0; j < N; j++) { const e = 1e-6 * Math.max(1, Math.abs(y[j])), yp = y.slice(); yp[j] += e; const Fp = f(0, yp); for (let i = 0; i < N; i++) Jm[i][j] = (Fp[i] - F[i]) / e; }
  const A = Jm.map((row, i) => row.map((v, j) => (i === j ? 1 / h : 0) - v));
  let dy; try { dy = solveLinear(A, F); } catch { h /= 4; it++; continue; }
  let ch = 0; for (let i = 0; i < N; i++) ch = Math.max(ch, Math.abs(dy[i]) / scale(y[i]));
  if (ch > 0.2 && h > 1e-9) { h /= 2; it++; continue; }
  y = y.map((v, i) => (COMP[i % 13] === 'SALK' ? v + dy[i] : Math.max(0, v + dy[i]))); day += h; if (ch < 0.05) h = Math.min(h * 2, 1e9);
  r0 = resid(y); it++;
  if (r0 < tol) { conv = true; break; }
}

COD, nitrogen and oxygen balances, sludge age

64-bit

Checks that what enters the plant equals what leaves plus what reacts, and gives the sludge retention time.

residual = in − out + reaction; rel = residual ÷ (|in| + 10⁻¹²); SRT = Σ V × solids ÷ (Qw × Xw + Qe × Xe)

Computed from the final state and the rate functions, not from the iteration's own numbers. The residual is a difference of large nearly equal totals, which is what it is meant to measure; it is reported unrounded and raises an error above 1e-3 on a converged run.

process/env/water.js
lines 124 to 137
CadCraft 0.8.2

let reactCOD = 0, oxygenUsed = 0, n2 = 0, supplied = 0;
T.forEach((t, i) => {
  const rho = rates(c[i], par);
  for (let j = 0; j < 8; j++) { let s = 0; for (const k of COD_SPECIES) s += M[j][k]; reactCOD += t.V * rho[j] * s; oxygenUsed += t.V * rho[j] * -M[j].SO; }
  n2 += t.V * rho[1] * (1 - par.YH) / (2.86 * par.YH);
  if (t.aerated) supplied += t.V * (t.KLa || 0) * (SO_SAT - c[i].SO);
});
const codOf = (o) => COD_SPECIES.reduce((s, k) => s + (o[k] ?? 0), 0), nOf = (o) => (o.SNO ?? 0) + (o.SNH ?? 0) + (o.SND ?? 0) + (o.XND ?? 0) + par.iXB * ((o.XBH ?? 0) + (o.XBA ?? 0)) + par.iXP * (o.XP ?? 0);
const codIn = Qin * codOf(plant.influent || {}), codOut = Qe * codOf(effluent) + Qw * codOf(waste);
const nIn = Qin * nOf(plant.influent || {}), nOut = Qe * nOf(effluent) + Qw * nOf(waste);
const codRes = codIn - codOut + reactCOD, nRes = nIn - nOut - n2;
const balances = { cod: { in: codIn, out: codOut, reaction: reactCOD, residual: codRes, rel: codRes / (Math.abs(codIn) + 1e-12) }, n: { in: nIn, out: nOut, n2Lost: n2, residual: nRes, rel: nRes / (Math.abs(nIn) + 1e-12) }, oxygen: { used: oxygenUsed, supplied, residual: supplied - oxygenUsed } };
const solids = (o) => SOLIDS.reduce((s, k) => s + o[k], 0), inv = c.reduce((s, ci, i) => s + T[i].V * solids(ci), 0), leaving = Qw * solids(waste) + Qe * solids(effluent);
const srt = leaving > 0 ? inv / leaving : Infinity;

Oxygen used and supplied

64-bit

The oxygen the reactions use next to what the aerators supply, in kg per day.

kg/d = g/d ÷ 1000

Two divisions by 1000; the returned numbers are unrounded. toFixed(1) is used only in the warning text, and the warning fires when supply is below 95 % of use.

process/env/water.js
lines 160 to 165
CadCraft 0.8.2

/** Oxygen: requirement (kg/d) next to what the aerators supply, as Problems when they disagree. */
export function oxygenReport(res) {
  const used = res.balances.oxygen.used / 1000, supplied = res.balances.oxygen.supplied / 1000, ps = [];
  if (res.converged && supplied < 0.95 * used) ps.push(problem('oxygen', 'warning', `Aeration supplies ${supplied.toFixed(1)} kg O2/d but the reactions use ${used.toFixed(1)} kg/d.`, [{ kind: 'plant', id: 'plant' }]));
  return { usedKgPerDay: used, suppliedKgPerDay: supplied, problems: ps };
}

In detail

Method: IAWQ Activated Sludge Model No. 1 (13 components, 8 processes) as listed in the file's SOURCES; pseudo-transient continuation, (I/h − J)·Δy = F.
Number types: 64-bit numbers in plain arrays and objects.
Tolerances: Stop when the largest |rate| ÷ (|value| + 1) is below 1e-7 (plant.tol) or after 4000 iterations. Jacobian step 1e-6 × max(1, |y|). Step halved when the largest relative change exceeds 0.2. Balance error flagged above 1e-3 relative. Ratios guarded at 1e-12 g/m³.
Limits: The file's REVIEW constant states the matrix and parameter set are not confirmed against the ASM1 report. No settler model, no temperature correction.

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.

Every policy that takes this approach