4D Midstream Digital Twin — engineering report

Version 0.1.0 · 2026-09-20 · midstream/ in the octopusgames workspace · dev server port 5177 (npm run dev:midstream)

Product boundary. This is a digital-twin demonstration over a fictional network with synthetic telemetry. It is not a certified pipeline simulator, leak-detection system, process-safety system, SIL-rated protection, emergency-shutdown controller, operational control system or engineering design tool. Nothing in it presents a simulated output as an operational command: model proposes, evidence is shown, the human operator decides. Every displayed value carries one of SIMULATED · OBSERVED · ESTIMATED · DERIVED · ASSUMED · PREDICTED.

1. Architecture

src/sim/runtime/     shared twin runtime kernel (no domain semantics)
  rng.ts             hash32 / randN — stateless seeded noise (seed, id, stepIndex)
  schedule.ts        piecewise input schedule with [t, until) entries and origin (scenario | operator | branch)
  trajectory.ts      fixed-step Trajectory<S, I>: stepper, snapshots every N steps, branch(i, schedule)
  invariants.ts      checkInvariants(list, state, t) → report; NaN helpers
  divergence.ts      Variable<S> registry + firstDivergences(): first persistent thresholded difference

src/sim/network/     network schema v1 (lines, segments, stations, pumps, compressors, valves, tanks, sensors)
src/data/network.ts  Synthetic Midstream Network #001 (fictional): 255 km main + 30 km branch (liquid), 40 km gas demo line
src/sim/hydraulic/   DemoHydraulicModel (inputs.ts = operator inputs + hidden faults + model params)
src/sim/sensors/     SyntheticObservationAdapter → canonical Observation records + data-quality rules
src/sim/estimation/  DemoStateEstimator (model + quality-weighted, distance-weighted residual blend)
src/sim/risk/        risk engine: evidence channels → RiskEvidence per segment, hypotheses, localisation, equipment trends
src/sim/causal/      causal chain / time-to-effect / incident timeline / causal-graph activation; ASK THE TWIN
src/sim/scenarios/   scenario catalogue (G01–G15) and counterfactual branch templates
src/twin/            MidstreamTwin orchestrator (world + expected trajectories, frames, prediction, branching, exports) and contracts
src/world/layout.ts  km → 3-D mapping, terrain, station pads, camera presets
src/scene/           React-Three-Fiber renderer (consumes twin frames only)
src/ui/              industrial UI (consumes twin frames only)
src/tests/           vitest: model, twin (acceptance 1–11), golden G01–G12
scripts/             verify.mjs (browser acceptance), capture.mjs (visual set), profile.mjs (frame pacing)
tools/               check-assets.mjs + asset-policy.json (transmission/volume/dispersion regression)

Single source of truth. MidstreamTwin owns two Trajectory objects — the synthetic WORLD (faults on) and the twin's EXPECTED model (same operator inputs, faults off) — and one AnalysisFrame per 15 s snapshot: observations, estimated pressure field, estimated segment inventory, risk evidence, hypotheses, localisation, equipment health, network summary. The renderer and every panel read frames at the viewed time; scrubbing the timeline changes only the view cursor, never history.

The twin loop as implemented

WORLD (synthetic)  ──observeAll──▶  OBSERVATIONS (+ quality)  ──estimate──▶  ESTIMATE  ──analyseRisk──▶  RiskEvidence / hypotheses / localisation
EXPECTED (twin)    ──────────────▶  expected per sensor  ─────┘                 │
                                                                                 ├──▶ 3-D overlays, panels (RENDER)
                                                                                 ├──▶ causal chain, incident timeline, ASK (ANALYSE)
                                                                                 └──▶ predict(): estimated state + estimated loss → BASELINE / A / B / C (PREDICT, COMPARE)

2. Simulation model (DemoHydraulicModel 1.0.0)

Staggered 1-D grid per line: pressure P_i and temperature T_i per cell (volume A·dx), mass flow ṁ_j per face.

dP_i/dt = (c² / V_i) · (ṁ_i − ṁ_{i+1} − ṁ_leak,i − ṁ_tee,i)        c² = K/ρ (liquid, K = ρ a²), Z R T / M (gas)
dṁ_j/dt = (A / dx) · (P_{j−1} − P_j − ΔP_friction − ρ g Δz + ΔP_pump − ΔP_valve)

• Friction Darcy–Weisbach f · dx/D · ṁ|ṁ| / (2 ρ A²), linearised implicitly (ṁ⁺ = (ṁ + coef·drive) / (1 + coef·R·|ṁ| − coef·b′)), with a viscosity proxy f·(1 + 0.0035·(25 − T)).
• Pumps affinity-law head H = s²·h0·(1 − 0.15·deg) − k·Q² with the curve slope b′ treated implicitly (explicit treatment was unstable: the slope is ~25× the friction slope at the face). PS-A = two parallel units on one face, PS-B booster, check valve when stopped.
• Compressor ratio curve r = 1 + 0.55·s²·(1 − 0.12·deg) − COMP_K·ṁ², power from the isentropic head with γ = 1.3, aftercooler to 45 °C.
• Valves K(x) = kOpen + 400·((1 − x)/x)²; stroke at a fixed rate; below 0.4 % open the face is closed. Receipt valves are sized from the design profile so design flow meets design boundary pressure.
• Leak orifice sink Cd·A·ramp(t)·√(2ρ(P − P_atm)) at a cell; the ramp represents a developing failure.
• Temperature upwind advection + relaxation to ground (τ = 6 h liquid, 3 h gas).
• Tanks level ODE per tank (source → T-101/T-102 → PS-A suction head; receipts → T-401/402; loading out of T-401–403).
• Slack line a cell cannot fall below 0.6 atm; at the floor its outflow is capped at its inflow and any mass created by the clamp is tracked in cum.slack so the audit balance stays exact.
• Numerics Δt = 1.25 s, cells 2.5 km (liquid) / 2 km (gas), wave speed 1000 m/s → CFL 0.5. 24 h of world + expected costs ≈ 0.3 s of stepping. Analysis frames every 15 s cost ≈ 1.2 ms each.

Simplifications (explicit). Isothermal-ish; no column separation, batches, DRA, viscosity-pressure coupling, valve dynamics beyond a stroke rate, pump start/stop hydraulics beyond a 25 s speed ramp, or elevation-driven slack-line vapour cavities. The expected model uses the same equations as the world plus a ±0.3 % ASSUMED friction-mismatch term.

3. Network model

34 main segments (S01–S34, 7.5 km each, 3 cells), 4 branch segments (B01–B04), 5 gas segments (G01–G05): 134 state samples. Stations: Inlet Terminal, Tank Farm (4 tanks), PS-A (P-101/P-102), BV-01, MS-01 (3 meter runs + prover), BV-02, PS-B (P-301, isolation valve V-301), BV-03, TEE-01 → branch → BD-01, BV-04, Delivery Terminal (3 tanks, loading gantry), CS-01 (C-601) → 40 km gas line → GR-01, Control Centre. 47 sensors: 20 PT, 8 FT, 4 TT, 4 VS, 4 TB, 7 LT. Corridor geography is compressed 1 km = 22 m in the 3-D world; station footprints are true scale.

4. State variables and provenance

The asset-inspection panel shows provenance, source sensor, quality, expected value, model version and residual for every number (see InfoPanel.tsx).

5. Sensor model and data quality

observeAll() produces the canonical Observation { timestamp, sensorId, assetId, variable, value, unit, quality, provenance, expected, residual, sigma }. Noise is σ_sensor · randN(seed, id, stepIndex); faults are scheduled in the same input schedule as everything else (drift, failed, frozen, spike, noisy, stale, comms loss). Quality rules: STALE (no update > 45 s), FAILED (out of physical range, or frozen ≥ 3 min), SUSPECT (implausible step, or a persistent isolated disagreement with neighbours while the mass balance closes — the cross-check flag set by the risk layer). sigma = √(σ_sensor² + σ_model²) with ASSUMED model error per kind.

6. State estimation

estimate() returns per-cell pressure = expected + Σ w_k·residual_k / Σ w_k over pressure sensors of the line, w = quality·exp(−|Δkm|/25) (GOOD 1, SUSPECT 0.3, FAILED/STALE 0), a leave-one-out cross-check per sensor, and segment inventories from the estimated field. It is honest about what it is: a weighted blend, not a Kalman filter; the interface (Estimate) is what a rigorous estimator would implement.

7. Anomaly logic (risk engine)

Per segment, seven channels are computed every frame from the last minutes of frames: pressure residual (interpolated between the bracketing GOOD/SUSPECT sensors, 2-min means), flow imbalance across the flow-meter bracket (Q_up − Q_down − offtakes − expected − d(unexplained linepack)/dt, 3-min mean), accumulated unaccounted volume (forgetting integral), persistence, spatial corroboration (GOOD sensors within 45 km agreeing in sign), rate of change, sensor quality. The score is

score = 100 · sat(0.30·pr·c·loc + 0.28·fi + 0.15·inv + 0.15·pers·loc + 0.12·spatial·loc) · (0.5 + 0.5·qualityFactor)
c   = 0.35 + 0.65·max(spatial, fi)          a lone sensor cannot carry the score
loc = 0.2 + 0.8·fi (or 0.6 without a usable bracket)   pressure evidence is weighted fully only where mass is also missing

States NORMAL < 18 ≤ WATCH < 40 ≤ ELEVATED < 65 ≤ HIGH with 5-point hysteresis; INSUFFICIENT when fewer than a third of the contributing sensors are GOOD; fail-closed to WATCH while history is shorter than 2 min. Known operational changes (operator or branch schedule entries) down-weight the mass-balance channels to 40 % for 8 min because the linepack correction is unreliable during fast transients. Only the unexplained linepack change (5-min EMA of the residual field) is subtracted — the expected flow difference already contains the expected linepack change (subtracting it twice produced false alarms during shut-ins; fixed).

Hypotheses (process anomaly, sensor drift, sensor failure, equipment degradation, demand transient, data quality) are scored from support/contradict facts and shown with WEAK/MODERATE/STRONG plus the WHY lists. Localisation combines residual arrival order (inverse-delay weighting) with the flow-meter bracket → core ± 3.5 km and broad region, labelled ESTIMATED with a confidence.

8. Risk, forecast, early warning

RiskEvidence is the renderer/UI contract (segmentId, timestamp, anomalyScore, riskState, evidence[], evidenceStrength, modelConfidence, uncertainty, persistenceMin). The demonstration risk matrix separates a likelihood indicator (score) from a consequence indicator (estimated loss over 60 min). Forecast: risk state at +5/+15/+30/+60 from the baseline prediction, time-to-next-state, leading indicators from the active channels. Early warning requires two channels ≥ 0.3, a sustained score ≥ 10 and a rising trend while the state is still NORMAL/WATCH.

9. Prediction and counterfactuals

predict(branch, fromT) builds the initial state from the EXPECTED snapshot with pressures replaced by the ESTIMATE, injects the twin's estimated loss (from the flow-imbalance channel, located at the localisation core) and steps the nominal model with the branch schedule appended. Three members (loss × 0.5 / 1 / 1.5) form the SIMULATED UNCERTAINTY BAND. Predicted risk is obtained by running the same observation/estimation/risk pipeline on noise-free predicted values against a nominal branch (every 4th frame). Branches: BASELINE, A = SIMULATE THROUGHPUT REDUCTION −20 %, B = SIMULATE SEGMENT ISOLATION km 120–127.5 (V-102 + V-301 closed, pumps stopped), C = SIMULATE PUMP SETPOINT CHANGE (stop PS-B). All start from the exact same frame; history is immutable (fingerprint() is unchanged by predicting). branchDivergences() derives what changed, when, where from the two trajectories.

A finding worth stating: a shut-in traps linepack in the enclosed segment, so the loss rate briefly rises before the segment depressurises (~25–35 min at this wall/fluid stiffness). The isolation branch therefore wins on end-of-horizon loss rate and preserved upstream inventory, not necessarily on cumulative volume within 60 min.

10. Causal analysis

causalChain() derives, from frames after the event, the first persistent occurrence of: local pressure deviation (PT ≥ 2 σ 2-min mean), flow residual (bracket ≥ 2 σ), neighbour corroboration (≥ 2 PTs), inventory imbalance (≥ 5 m³), WATCH, ELEVATED, HIGH. deriveIncidentEvents() builds the structured incident timeline (sensor first residuals, quality changes, escalations, equipment trends, operator interventions). The causal graph lights only evidence-supported nodes (graphActivation). No timestamps are hardcoded anywhere in the UI.

11. 3-D architecture

R3F scene reading runtime.viewT. Ground (heightfield with scanned grass + macro, x-ray sectioning shader around the camera focus), Pipeline (tubes that surface at pads; ribbon overlays coloured by per-cell field textures filled from frames; risk envelope; instanced sensor markers coloured by quality with replayable first-residual pulses; marker posts; optional ground-truth leak marker), Stations (procedural true-scale geometry merged per material bucket, floating roofs / gauges / couplings / valve indicators driven by the viewed world state), ControlRoomScreens (canvas textures), Vegetation (CC0 scans instanced < 260 m, baked billboard impostors beyond), Agents (instanced workers on inspection loops, road vehicles), Lighting (sun/sky/fog/weather, pooled night spot lights), adaptive DoF only while inspecting and still, AUTO quality drops dpr after two slow seconds.

12. Performance

Measured 2026-09-20 with node scripts/profile.mjs (headless Chrome, Metal, vsync off, 8 s samples, scenario G07 at 10×). Frame pacing from the app's own ring buffer; draw calls count the whole frame including post passes.

Simulation CPU is separate from rendering: one 15 s analysis frame (12 hydraulic steps for world + expected, observation, estimation, risk) costs ≈ 1.2 ms and runs at most 33× per real second at 500×; at 10× it is ≈ 0.02 ms per rendered frame. The initial 1 h 45 of history (421 frames) costs ≈ 0.3–1.1 s at load. A 60-min prediction branch (3 ensemble members + risk every 4th frame) costs 0.4–0.5 s and is queued one per animation frame. Triangle counts are dominated by near-camera CC0 tree scans (up to 400 instances per model); beyond 260 m trees are single-quad impostors. The worst 0.1 % frames at dpr 2 stay under 26 ms; the AUTO quality tier drops to dpr 1 after two seconds below 42 fps. Raw table: .build/profile/report.md. frame pacing (mean, 1 % low, 0.1 % low, worst, frames > 33/50 ms), draw calls (whole frame incl. post passes), triangles, simulation CPU and agent CPU at dpr 1 and dpr 2 for static / moving / facility / follow. Draw calls are counted by resetting gl.info in the ticker (three's own counter only reflects the last post pass). Asset regression: npm run test:assets fails on unapproved KHR_materials_transmission / volume / dispersion.

13. Tests

npm test → 43 vitest tests: model (steady operation, continuity, tank conservation, pump/valve/demand response, leak, compressor, degradation, determinism), twin (observation/quality/comms loss, normal operation, drift vs loss, propagation ordering, false-positive resistance, known operator changes, demand surge, pump degradation trend, multiple events, branching + recovery, provenance/manifest/ask/exports) and golden G01–G12 (fingerprint determinism + qualitative outcomes + event ordering). npm run verify → 42 browser checks (see .build/verify/results.json); npm run capture → 24 acceptance captures + contact sheet in .build/acceptance/.

14. Known limitations

See the in-app MODEL LIMITATIONS panel (LIMITATIONS in src/ui/AskReport.tsx); highlights: fictional network, synthetic telemetry, simplified hydraulics, same-equation expected model with a small mismatch term, weighted-blend estimator, heuristic score (not a probability), demonstration risk states, estimated (not confirmed) localisation, illustrative integrity values, simulation-only counterfactuals, mass-balance channels down-weighted after known transients, compressed corridor geography. Sensor spacing (≈ 20 km) limits pressure-based localisation to roughly the segment scale; the flow-meter bracket does the heavy lifting, and all segments inside the affected bracket (up to 7) share the flow evidence, so the elevated region is a bracket, not a single segment. Visual realism is a first pass: terrain relief and the right-of-way read weakly from the drone view, and night lighting is deliberately dark outside the stations.

15. Future real telemetry integration

The renderer and analysis consume AnalysisFrames built from Observation[] + an Estimate; the synthetic generator is one observeAll() call. Replacing it:

SCADA / HISTORIAN ──▶ ObservationAdapter (OPC UA · MQTT · historian/TSDB · CSV / Parquet replay · REST)
                       └─ produces Observation { timestamp, sensorId, assetId, variable, value, unit, quality, provenance: 'OBSERVED' }
                    ──▶ quality / validation (same rules or the site's)  ──▶ estimate(expected, obs)  ──▶ analyseRisk(...)  ──▶ frames

Adapter interface (to be implemented, not connected here): interface ObservationAdapter { poll(t): Observation[]; sensors(): SensorDef[] } — OPC UA (tag → sensorId map, quality codes → GOOD/SUSPECT/FAILED/STALE), MQTT (topic per tag), historian/TSDB (batched time-range reads for replay), CSV/Parquet replay (the exported sensors.csv schema is the same record), REST. The expected model then needs measured boundary conditions instead of the synthetic operator schedule (inlet flow / outlet pressure as Inputs channels), and DynamicsModel (Stepper<PhysState, Inputs>) can be swapped for a validated simulator without touching trajectories, branching, timeline, evidence UI or provenance. DemoStateEstimator and the risk rules are likewise behind their result contracts (Estimate, RiskEvidence).

16. Shared runtime validation (third domain)

Plant (STATE(t), integer days) and Factory (STATE(t), fractional hours) were the first two domains. Midstream tests STATE(x, t) with topology, observations and estimation. Result:

What did not generalise: spatial state samples, sensors/quality, estimation, risk evidence and localisation are all domain objects (as the architecture study predicted). Two trajectories per twin (world vs expected) is a new pattern that the runtime supported without change. The runtime stays ~250 lines with no three.js/React dependency; nothing was extracted into a shared package yet — the same five files could be lifted verbatim.