Example: Pipeline-Leak Community Health Risk

This chapter builds a complete prophys model step by step. Its defining idea: the hazard source is not a fixed point. A pipeline can fail anywhere along its route, so the leak location is a continuous, reparameterized random variable — not a hand-picked point. On top of that, wind speed and direction are two more random inputs. All three are marginalized out together by Monte Carlo, through a Gaussian-plume dispersion model and a dose-response curve, to give the expected health-complaint risk everywhere on the map and the tail exposure at individual residences. The full runnable source is examples/gas_dispersion/model.py.

Step 1 — The route and receptors

Everything spatial lives in one site frame. The pipeline is a Polyline — a single connected route whose vertices could themselves be a Param (optimizable siting) — and the residences are a PointList:

import jax.numpy as jnp
import prophys as prp

site = prp.Frame("site", units="m")
pipeline = prp.Polyline(
    jnp.array([[-550., -320.], [-150., -60.], [280., 40.], [650., -90.], [1000., 180.]]),
    site,
)
residences = prp.PointList(
    jnp.array([[420., 300.], [640., 120.], [-250., 380.], [520., -260.]]), site
)

Step 2 — Uncertainty enters at the inputs

The leak location is leak_frac ~ Uniform(0, 1) pushed through Polyline.point_at — a reparameterized draw of a position anywhere along the route, differentiable in both the fractional arc-length and the route geometry. Wind speed is Weibull, wind direction von Mises (the circular analog of a Gaussian — a plain Gaussian would be wrong because 0° and 360° are the same direction):

leak_frac      = prp.RandomVariable("leak_frac", prp.Uniform(0.0, 1.0))
wind_speed     = prp.RandomVariable("wind_speed", prp.Weibull(scale=6.0, concentration=2.0))
wind_direction = prp.RandomVariable("wind_direction", prp.VonMises(loc=prp.deg2rad(35.0), kappa=3.0))

leak_xy = pipeline.point_at(leak_frac)   # leak position, per sample

Every quantity downstream of these three leaves is itself random. The compiled model marginalizes them out by Monte Carlo, and — because all three sample by reparameterization — gradients flow through the marginalization back to any trainable parameter.

Step 3 — The physics: a Gaussian plume, symbolically

Ground-level concentration follows the standard Gaussian dispersion model. Each receptor is rotated into plume-aligned coordinates using the random wind direction, the plume spread grows with downwind distance via differentiable TableLookup1D Pasquill–Gifford coefficients, and the emission rate carries a trainable, bounded abatement fraction:

delta = receptor_coords - leak_xy
dx, dy = delta[..., 0], delta[..., 1]
downwind  =  dx * prp.cos(wind_direction) + dy * prp.sin(wind_direction)
crosswind = -dx * prp.sin(wind_direction) + dy * prp.cos(wind_direction)

sigma_y = prp.TableLookup1D([50, 200, 500, 1000, 2000], [8, 25, 55, 100, 180])(prp.relu(downwind) + 1e-3)
sigma_z = prp.TableLookup1D([50, 200, 500, 1000, 2000], [5, 15, 30, 50, 80])(prp.relu(downwind) + 1e-3)

abatement = prp.Param("abatement", init=0.2, bounds=(0.0, 0.9))
Q = (1.0 - abatement) * prp.Param("emission_rate", init=5.0)
conc = (Q / (2 * jnp.pi * u * sigma_y * sigma_z)) \
       * prp.exp(-0.5 * (crosswind / sigma_y) ** 2) \
       * prp.exp(-0.5 * (release_height / sigma_z) ** 2)

Every line is symbolic — jax.grad differentiates through the rotation, the table lookups, and the exponentials in one pass, and the whole thing is shape-agnostic: the same function serves a single residence (broadcasting against the Monte-Carlo sample axis) or a dense grid slice.

Step 4 — Uncertain attributes, tail risk, and a correlation

Exposure at a residence is LogNormal around the physical concentration (multiplicative monitoring noise with a learnable sigma). Because every attribute is marginalized over the joint (leak, wind) draws, the engine gives not just the expectation but the Conditional Value-at-Risk — the mean of the worst 5 % of outcomes:

exposure = prp.UncertainAttribute(
    "exposure_R1",
    prp.LogNormal(mu=prp.log(prp.clip(conc_R1, 1e-6, None)), sigma=noise_sigma),
    unit="g/m3",
)
...
compiled.expectation("exposure_R1")      # mean exposure
compiled.cvar("exposure_R1", alpha=0.95) # tail exposure

Two nearby residences share the same weather and leak, so their exposures are correlated even conditional on the mean model — captured with a trainable Correlation (a bivariate Gaussian copula) in an AttributeInteractions.

Step 5 — Risk landscape, calibration, design, export

Because every quantity is an ordinary symbolic expression, evaluating it over a dense Grid while sampling the same random inputs turns it into the marginalized risk landscape at the top of this page — the hazard smeared along the whole route, blown downwind. The compiled model then calibrates the monitoring-noise sigma against observed exposures, optimizes the abatement fraction (its [0, 0.9] box enforced structurally by the bounded Param, not a penalty), and exports to the neutral ModelPackage.

Primitives exercised

Category	Used for
Polyline	The pipeline route; point_at makes “where along it” a reparameterized draw
PointList / Grid	Residence receptors; dense risk-map evaluation
Uniform, Weibull, VonMises	Leak location, wind speed, and wind direction as RandomVariable inputs
TableLookup1D, PiecewiseLinear	Dispersion coefficients; dose-response
LogNormal, Bernoulli, Correlation	Exposure noise; complaint likelihood; shared-weather coupling
expectation / cvar, calibrate(), optimize(), ModelPackage	Expected & tail risk, training, design, export

Running it

python examples/gas_dispersion/model.py