Stochastic Modeling — Heterogeneity, Uncertainty, and the Recursive-Statistics Architecture

19.1 What Stochastic Modeling Is — Single Realizations and Ensembles

Stochastic modeling in groundwater practice covers a broader territory than Monte Carlo simulation alone. Both single-realization analysis and ensemble simulation belong here; IGW-NET supports both.

19.1.1 The two forms of stochastic analysis

Two kinds of stochastic work address different questions:

These are both stochastic analyses — both treat spatial inputs as random variables. Single-realization shows the effect of randomness (the qualitative structure of heterogeneous flow); ensemble quantifies uncertainty under that randomness (the quantitative spread of possible outcomes). Many projects benefit from both: start with a single realization to understand heterogeneity effects; escalate to ensemble when uncertainty quantification is needed.

19.1.2 Why single-realization is valuable on its own

A common misconception: "stochastic modeling requires many realizations." Not always. A single well-chosen realization can be a valuable analysis by itself:

• It shows realistic flow physics. Head contours distort around low-K zones; flow vectors bend and channel through high-K corridors; particle paths twist and diverge rather than following smooth uniform-K trajectories. These are visual features of real aquifers that uniform-K models systematically miss.
• It teaches the role of heterogeneity. For students, clients, or collaborators unfamiliar with groundwater heterogeneity, a single realization is the most direct demonstration. "This is what your aquifer actually looks like, roughly, in terms of its flow structure. This is why smooth-K predictions can be misleading." Synthetic mode (Ch. 1 §1.3.2) is especially effective here — enter it via Utilities → Go to Synthetic Case Area, and you can toggle heterogeneity on and off with a single click, watching the dramatic effect of K structure on transport immediately, without real-site complications obscuring the lesson. The emergence of macrodispersion from small-scale heterogeneity — the nonlinear scale effect — is vivid right in front of the class. Work through Tutorial 15 (Single-Realization Stochastic Flow) and Tutorial 17 (MC Transport) for the full recipe.
• It is one hypothesis — useful in its own right. The real K field in the ground is unknown; a realization is one physically-plausible arrangement consistent with the available statistics. For planning, design, or "what if the aquifer is like this?" questions, one realistic hypothesis often does what a full ensemble would.
• It can be run quickly. A single realization costs essentially what a single deterministic simulation costs — no ensemble overhead. For exploratory work, this is the right level of effort.

Move to ensemble Monte Carlo when the question specifically requires uncertainty quantification — probabilistic capture zones, breakthrough curve distributions, confidence bounds for risk-based decisions. Don't run an ensemble when a single realization would answer your question.

19.1.3 What the rest of this chapter covers

Section 19.2 establishes spatially correlated random fields as the conceptual framework uniting SGS (continuous) and T-PROGS (categorical). Section 19.3 covers SGS setup in detail: which variables can be randomized, what parameters govern the random field, how per-zone configuration works. Section 19.4 walks through single-realization analysis. Section 19.5 covers the full ensemble Monte Carlo workflow. Section 19.6 is the distinctive IGW-NET contribution: the recursive ensemble-statistics architecture that makes arbitrary realization counts feasible at full model resolution. Section 19.7 covers result interpretation — mean fields, variance maps, probability distributions, macrodispersion, probabilistic capture zones. Section 19.8 connects back to calibration (Ch. 18) for parameter uncertainty as stochastic input. Section 19.9 collects common pitfalls, including the critical one that stochastic analysis does not quantify structural uncertainty.

19.2 Spatially Correlated Random Fields — The Unifying Framework

IGW-NET's two main stochastic methods — Sequential Gaussian Simulation (this chapter) and T-PROGS transition probability (Ch. 17) — are different operational techniques within a common conceptual framework: both produce spatially correlated random fields. Understanding the framework helps clarify when to use which method.

19.2.1 What "spatially correlated random field" means

A spatially correlated random field is a collection of values, one at each cell in the domain, with two key properties:

• Random. Values are drawn from probability distributions, not prescribed deterministically. Running the generation with different random seeds produces different-looking fields.
• Spatially correlated. Neighboring cell values are not independent — they're related through a specified correlation structure. Cells within the correlation scale tend to have similar values; cells beyond the correlation scale are effectively independent.

The correlation structure is what makes the field "realistic." Completely independent random values at every cell would produce pure noise — a salt-and-pepper pattern with no geological meaning. Real aquifers don't look like that. They have regions of high-K sand grading to regions of low-K clay, with characteristic length scales (the correlation scale) determined by depositional processes. The correlation structure in the random field captures that reality.

19.2.2 Continuous vs. categorical random fields

Two fundamental types of spatially correlated random fields:

19.2.3 The common workflow pattern

Both methods share the same high-level operational pattern:

From §19.3 onward, the focus is on SGS specifically. T-PROGS is covered in its own chapter (Ch. 17). Where appropriate, ensemble-statistics sections apply to both methods identically.

19.3 Sequential Gaussian Simulation — Continuous Random Fields

Sequential Gaussian Simulation (SGS) is the standard geostatistical method for generating continuous random fields. IGW-NET supports SGS for any continuous spatial input variable — K is the most common use, but porosity, recharge, retardation, and elevations are equally valid.

19.3.1 Which variables can be random fields

In IGW-NET's zone-level configuration, any of these variables can be flagged as random:

Each zone can independently configure which of these are random. A common setup: K random everywhere; porosity random in the AQ zone only (where it matters most for transport); elevations deterministic from the Data Center. Flexibility matches the geological and hydrologic reality of the site.

19.3.2 The random field parameters

Each random field is specified by a small set of statistical parameters, configured per-zone-per-variable:

19.3.3 The correlation-scale choice matters

Of all the parameters, the correlation scale (Lambda) is the most impactful on model behavior and the most often misconfigured. Small LambdaX (say 10 m in a 10-km model) gives near-independent noise — salt-and-pepper heterogeneity that doesn't look like geology. Very large LambdaX (say 10 km in a 10-km model) gives fields that are essentially constant across the domain — not much heterogeneity at all.

19.3.4 Log-normal vs. normal K

Hydraulic conductivity is conventionally modeled as lognormally distributed — log K is Gaussian, K itself is lognormal. This matches empirical observation: K in natural aquifers spans many orders of magnitude, with values clustered in the middle of log space and tails extending to very high and very low values. SGS in IGW-NET operates on log K by default for hydraulic conductivity, converting to K via exponentiation for the flow simulation. The Mean and Variance parameters, for K, refer to the log-K mean and variance respectively.

Other variables (porosity, recharge) may be closer to normally distributed depending on their underlying physics. IGW-NET's default treatment is generally sensible; advanced users can override the distribution type where needed.

19.4 Single-Realization Analysis — What Heterogeneity Does

Before escalating to ensemble Monte Carlo, spend time with single-realization simulations. They show what heterogeneity does to flow and transport in ways that uniform-K modeling cannot, and they're cheap — essentially the cost of one deterministic simulation. Single-realization analysis is the gateway to stochastic thinking.

19.4.1 What to expect from a single realization

Compare the single-realization result against a uniform-K baseline run of the same model. You'll see, qualitatively:

• Head contours irregular — they distort around low-K bodies where flow is impeded and concentrate in high-K channels where flow runs freely
• Flow vectors bend and channel — rather than following smooth parallel paths from source to sink, vectors deflect around low-K zones and concentrate in high-K pathways, producing spatially-variable magnitudes
• Particle paths tortuous — paths twist as they navigate the random structure, with some particles moving quickly through high-K corridors and others taking meandering slow routes
• Travel times spread — particles that started together arrive at outlets at very different times, reflecting the diversity of path conductivity experienced

19.4.2 What a single realization is and is not

Critical framing:

• Is: one physically-plausible arrangement of the aquifer consistent with the specified statistics
• Is: a visualization of heterogeneity effects that uniform-K modeling cannot show
• Is: a useful hypothesis for planning, design, or exploratory analysis
• Is NOT: "the truth" about the actual aquifer in the ground
• Is NOT: a prediction with quantified uncertainty — without an ensemble, you don't know how this specific realization compares to other possible realizations
• Is NOT: sufficient for regulatory or risk-based decisions that require uncertainty bounds

19.4.3 When single-realization work is enough

Single realization suffices when the question is:

• "What does flow through a realistic heterogeneous aquifer look like?" (teaching, client communication)
• "Can low-K zones plausibly deflect a contaminant plume in unexpected directions?" (exploratory)
• "Does the heterogeneity structure produce preferential flow that could route contamination to a specific well?" (screening)
• "Is the flow field qualitatively sensitive to heterogeneity?" (testing whether ensemble analysis is warranted)

Escalate to ensemble when the question requires probability or risk language — confidence intervals, exceedance probabilities, percentile bounds, risk-weighted decisions. Those require the ensemble (§19.5).

19.5 Ensemble Monte Carlo Simulation

Ensemble Monte Carlo runs many realizations — tens to thousands — and aggregates them into ensemble statistics. This is where uncertainty quantification happens. The setup extends the single-realization setup with a realization count and ensemble output options.

19.5.1 The Monte Carlo setup workflow

19.5.2 Watching the ensemble form

19.5.3 How many realizations

Common question without a rigid answer. Guidelines:

• For mean fields (ensemble mean head, ensemble mean concentration): 50-100 realizations often converge adequately. Watch the mean field visualization; when it stops changing noticeably as new realizations add, you've converged.
• For variance fields (head variance map): 100-500 realizations typical. Variance converges more slowly than mean.
• For probability tails (exceedance probabilities, 5th/95th percentile): 500-1000+ realizations. The rarer the event you're trying to quantify, the more realizations you need to see it reliably.
• For probabilistic capture zones: 100-500 is typical; the fuzzy boundary structure stabilizes at this range.

Because IGW-NET's recursive-statistics architecture (§19.6) has no memory penalty for larger realization counts, there's no strong reason to under-run. When in doubt, run more. Watch convergence in the real-time visualization and stop when statistics stabilize.

19.6 The Recursive-Statistics Architecture — Why This Works at Scale

Ensemble Monte Carlo at fine grid resolution is traditionally prohibitive because storing the ensemble is prohibitive. IGW-NET's recursive-statistics architecture eliminates that constraint. This section explains the mechanism and why it matters.

19.6.1 The storage problem

Consider a realistic 3D groundwater model: 10⁶ cells (say 1000 × 1000 × 1 layer for a 2D model, or 200 × 200 × 25 layers for 3D). You run 1000 realizations of Monte Carlo simulation. The naive approach stores every realization's head field, then post-processes the ensemble to compute statistics:

• Storage per realization per variable: 10⁶ × 8 bytes = 8 MB
• Storage for 1000 realizations of heads alone: 8 GB
• Plus concentrations, fluxes, other outputs: quickly into tens of gigabytes
• Post-processing to compute statistics: read the whole ensemble back; compute mean, variance, percentiles; write statistics maps

At these sizes, the workflow starts choking. Disk I/O becomes the bottleneck. Memory pressure causes swapping. The user can no longer interactively explore the results — they're locked into a batch-processing pattern. And it only gets worse with finer grids (the grid cell count scales with the cube of linear resolution in 3D) or more realizations.

19.6.2 The recursive solution

IGW-NET never stores the ensemble. Instead, each realization updates running statistics in place, then is discarded. The running statistics have the same memory footprint as a single realization's output, regardless of how many realizations have been processed.

For the ensemble mean at cell (i,j,k) after N realizations:

mu_N = mu_{N-1} + (h_N − mu_{N-1}) / N

Each realization's head h_N is incorporated into the running mean mu and then discarded — no need to ever store h_N again.

For the ensemble variance, IGW-NET uses Welford's online algorithm, which computes running variance numerically stably without accumulating round-off errors (naive cumulative-sum-of-squares algorithms become unstable for large N):

M2_N = M2_{N-1} + (h_N − mu_{N-1}) × (h_N − mu_N)

variance_N = M2_N / N (or M2_N / (N−1) for sample variance)

Probability-of-exceedance fields work similarly: maintain a count of realizations exceeding a specified threshold at each cell, divided by N, updated per realization. Percentiles require slightly more work but are also computable via online algorithms or streaming quantile estimators.

19.6.3 What you can do because of this

Practical consequences of the recursive-statistics architecture:

• Run as many realizations as you want. 10, 100, 1000, 10,000 — the memory footprint is the same. Run whatever you need for statistical convergence; stop when visualizations stabilize.
• High-resolution 3D ensembles are tractable. Where traditional tools force you to coarsen the grid to fit the ensemble in memory, IGW-NET runs at full model resolution regardless of realization count.
• Real-time convergence monitoring. Watch the mean field smooth and the variance map stabilize as realizations accumulate. Decide when to stop based on what you see, not on a pre-specified count that might be too small or too big.
• Iterative exploration. Try a setup with 20 realizations for quick feedback; if you like what you see, continue with 200 more to refine statistics. The workflow is responsive at every step.
• Consistent with interactive modeling philosophy. IGW-NET's whole design philosophy (Ch. 18 §18.1.3) is about making model iteration fast and interactive. Stochastic analysis fits into this philosophy rather than forcing a batch-processing exception.

19.6.4 What trades off

The recursive approach has minor limitations worth noting:

• Specific percentile information is limited. Running statistics give mean, variance, and probability-of-exceedance at pre-specified thresholds naturally. Arbitrary percentiles (e.g. "what's the 27th percentile at this cell?") require either storing more information or accepting approximations.
• Retroactive analysis is limited. Once the ensemble is done, you cannot go back and ask "what do the realizations with low K look like in aggregate?" — those realizations aren't stored. For questions like this, selectively record subsets of realizations at specific points of interest (IGW-NET's Monitoring Probability checkbox does this for wells) or save a smaller ensemble sample.
• Debugging a specific unusual realization is harder since you can't revisit it. Use a deterministic seed for reproducibility if you want to re-run a specific case.

In practice these trade-offs are small compared to the massive gain in ensemble capacity.

19.7 Result Interpretation — From Ensemble to Insight

Running the ensemble is the mechanical step. Interpreting the ensemble statistics is where the value is extracted. This section covers the main output types and their interpretation.

19.7.1 Mean head field

The ensemble mean head at each cell is the expected head value averaging over the uncertainty in K. It's smoother than any single realization because the ensemble averages out realization-specific heterogeneity patterns. For most purposes (regional water budget, typical head pattern, average flux estimates), the mean is what you want.

19.7.2 Head variance map

Typical interpretation patterns:

• Low variance near boundaries: BCs pin the head regardless of K, so variance is minimal
• Higher variance in the interior: K has more room to affect heads when no BC is nearby
• Peaks near pumping: drawdown depends sensitively on local K, so regions near heavy pumping often have high variance
• Peaks at SW-GW boundaries: leakance and K both affect exchange, compounding uncertainty

19.7.3 Probability distributions at monitoring wells

Where Monitoring Probability was enabled, IGW-NET records head (and concentration, if active) values from every realization at that specific well, producing a probability distribution. The chart shows a histogram or density curve — "at this well, the head is most likely between X and Y, with some probability of being lower or higher."

Typical uses:

• "What's the probability that head at this well exceeds a specific threshold?" — read off the exceedance probability from the cumulative distribution
• "What's the 5th–95th percentile range for head at this well?" — confidence interval for engineering design
• "Is the well likely to go dry under a specified scenario?" — compute P(head below well bottom) from the distribution

19.7.4 Macrodispersion in transport ensembles

Macrodispersion is the excess spreading of the ensemble mean plume beyond any single realization's plume. Different realizations route contaminants along different preferential pathways; averaging them produces a smoother, wider mean plume. Macrodispersion is a fundamental property of heterogeneous aquifers — it cannot be captured by a uniform-K model with a tuned dispersivity parameter; it's a real physical phenomenon emerging from K heterogeneity.

Concentration variance at the plume fringe is diagnostic of contamination risk. Low variance in the plume interior (always contaminated); low variance far ahead of the plume (always clean); high variance in the transition zone (maybe contaminated, maybe not). The variance map highlights locations where small K-field differences produce large concentration differences — where the contamination risk is genuinely uncertain.

19.7.5 Probabilistic capture zones

A probabilistic capture zone replaces the deterministic "this polygon is inside the capture zone" with "each location has probability P of being captured." Useful applications:

• Wellhead protection area designation: protect all locations with capture probability above a chosen threshold (e.g., 10% for conservative, 50% for less conservative)
• Contaminant source tracking: given contamination at a specific well, what's the probability each potential source reached it?
• Source-control prioritization: among known contamination sources, rank by capture probability for the affected well — highest probability = most likely contributor

19.8 Connection to Calibration — Parameter Uncertainty as Stochastic Input

Stochastic modeling and calibration (Ch. 18) naturally connect. UCODE parameter estimation produces not just optimal parameter values but parameter uncertainty — standard errors, correlation matrices, confidence intervals on parameters. Those uncertainty estimates are directly usable as inputs to stochastic analysis.

19.8.1 UCODE output as stochastic input

After running Automatic Calibration with UCODE (Ch. 18 §18.7), you have:

• Optimized parameter values (e.g., K multiplier = 0.3)
• Parameter standard errors (e.g., K multiplier = 0.3 ± 0.05)
• Parameter correlation matrix (e.g., K multiplier and Recharge multiplier correlated at 0.8)

These together define a multivariate distribution of plausible parameter values consistent with the calibration data. You can use this distribution to drive stochastic analysis:

• Sample parameter values from the UCODE-derived distributions (accounting for correlations)
• For each sampled parameter set, run the model
• Aggregate outputs across samples into ensemble statistics

The result is uncertainty quantification that propagates calibration uncertainty into simulation outputs. A prediction (say, head at a new well site) now has a confidence interval derived from the calibration's own uncertainty — not an arbitrary guess.

19.8.2 Combined parameter and heterogeneity uncertainty

Advanced stochastic analyses combine multiple uncertainty sources:

• Parameter uncertainty (from UCODE): uncertainty in calibrated multipliers
• Heterogeneity uncertainty (from SGS/T-PROGS): uncertainty in the spatial K distribution
• Boundary condition uncertainty (where applicable): uncertainty in prescribed heads or fluxes

Each realization draws from all three sources jointly: one K multiplier sampled from UCODE distribution, one SGS K field generated with a different random seed, one BC value sampled if that's in play. The resulting ensemble captures all three sources of uncertainty together.

IGW-NET supports this compound stochastic analysis natively — the isRandom flags work at multiple levels (multiplier-level for UCODE-derived parameter uncertainty; zone-level for heterogeneity; boundary-level for BC uncertainty). The recursive-statistics architecture works identically regardless of how many uncertainty sources are active.

19.8.3 Parameter uncertainty vs. structural uncertainty

Before stochastic analysis, use IGW-NET's real-time modeling and visualization to verify the structure is defensible. Stochastic analysis then quantifies the residual parameter and heterogeneity uncertainty within that sound structure.

19.9 Common Pitfalls

Stochastic modeling has its own failure modes, distinct from deterministic modeling pitfalls. This section collects the common ones.

19.9.1 Single realization mistaken for the answer

A single realization is one hypothesis — one possible arrangement consistent with the statistics, not "the truth." Users who generate a single realization and interpret its specific patterns as real are over-reading. The patterns in one realization reflect that specific random seed; a different seed would produce different patterns.

Fix: either frame the single realization explicitly as "one plausible hypothesis" for exploratory or communication purposes, or escalate to ensemble analysis. Don't report single-realization results as "the stochastic answer."

19.9.2 Correlation scale misconfiguration

Setting LambdaX too small produces near-independent noise; too large produces near-constant fields. Both break the realism of the random field. Users who leave defaults without thinking about the geological setting often get unrealistic heterogeneity patterns.

Fix: match LambdaX, Y, Z to the geological feature lengths. Horizontal 100-2000 m for most sedimentary aquifers; vertical 1-10 m reflecting layer thicknesses. §19.3.3's tip callout covers this.

19.9.3 Too few realizations

Running 10 realizations and declaring the statistics converged is almost always wrong. Mean fields need dozens; variance maps hundreds; probability tails thousands. Premature stopping produces ensemble statistics that reflect the specific realizations chosen more than the underlying distribution.

Fix: watch convergence in IGW-NET's real-time visualization (§19.6 enables this). Continue until statistics visibly stop changing. With the recursive-statistics architecture, running more realizations has no memory penalty.

19.9.4 Treating realizations as independent when they shouldn't be

If parameter uncertainty and heterogeneity uncertainty are both active, sampling them independently may miss important correlations. For example, if UCODE finds K and Recharge strongly correlated (the non-uniqueness from Ch. 18 §18.8), sampling them independently produces (K, R) pairs that contradict the calibration. Honoring the correlation structure matters.

Fix: sample multi-variable parameter distributions with their covariance structure, not as independent univariate distributions. UCODE-derived distributions have correlation matrices — use them.

19.9.5 Structural uncertainty ignored

Covered in §19.8.3 and in the key callout. Monte Carlo quantifies parameter uncertainty within a fixed structure; structural uncertainty is bigger and not captured. Modelers who report tight confidence intervals from a Monte Carlo ensemble without acknowledging that the structure could be wrong are overselling.

Fix: always frame Monte Carlo results as "parameter uncertainty within this structure" and separately address structural uncertainty through sensitivity analyses (try alternative structures, compare results) or through explicit discussion of the structural assumptions.

19.9.6 Realistic appearance without observational backing

SGS produces visually impressive K fields that look realistic even when the specified statistics aren't based on anything real. Users who set variance, correlation, and mean to plausible-looking numbers without empirical justification get plausible-looking results that may have no defensible statistical basis.

Fix: ground the random field parameters in data where possible — borehole spacing analysis for correlation scales, pump test variability for K variance, published regional studies for typical values. When data is limited, explicitly acknowledge the parameters are assumptions.

19.9.7 Using ensemble mean as "the answer"

The ensemble mean is one useful statistic but not "the answer." Planning decisions based on the mean alone ignore the spread — decisions robust to the mean may fail under realistic realizations. Conversely, worst-case planning based on the tail alone ignores typicality — overly conservative decisions waste resources.

Fix: use the full statistical output. Plan for the mean behavior; stress-test against percentile bounds; make decisions explicitly risk-informed by the probability structure. The whole point of the ensemble is to enable this richer decision-making.