Open-Source Monte Carlo Tool Brings Uncertainty Into Wellhead Protection Zones

Groundwater quietly supplies roughly half of the world’s drinking water, yet it remains one of the most poorly protected resources on the planet. Once an aquifer is contaminated, remediation is notoriously difficult and expensive, which makes preventing pollution in the first place the only realistic strategy. A central instrument of that strategy is the wellhead protection area, or WHPA: a mapped zone around a pumping well within which land use is restricted so that contaminants never reach the water being extracted. The trouble is that drawing the boundary of such a zone requires knowing how water moves through the subsurface, and the subsurface is fundamentally uncertain. A new open-access study published in Discover Geoscience by Nadine Gärtner, Maryam Zamzami, and Andreas Lindhe, researchers at Chalmers University of Technology and KTH Royal Institute of Technology in Sweden, presents a practical, uncertainty-aware workflow that puts probabilistic WHPA delineation within reach of the small and medium-sized water suppliers who need it most.

The core problem the researchers tackle is that conventional WHPA delineation relies on deterministic models built from fixed parameter values. Hydraulic conductivity, aquifer thickness, and effective porosity are each estimated from sparse field data, and a single ‘best estimate’ of each is fed into a model that then produces one crisp boundary. Studies dating back to the 1990s have shown that such deterministic approaches can produce overly optimistic delineations that mask the true range of plausible capture zones. Because WHPA boundaries carry real consequences, including land-use restrictions and potential compensation for affected landowners, an artificially narrow zone can leave a drinking water source exposed while an inflated one can impose unnecessary burdens. Recent legislative shifts, notably in Sweden, now demand risk-based approaches that explicitly account for uncertainty, but the advanced numerical tools capable of doing so, such as MODFLOW coupled with groundwater modeling systems, demand data, expertise, and computational resources that many rural utilities simply do not have.

The workflow developed by the Swedish team combines two established techniques in a novel, accessible package. The first is the Analytic Element Method, or AEM, pioneered by Otto Strack, which represents groundwater flow by superimposing elementary analytical solutions for wells, rivers, lakes, and other hydrogeological features. Because AEM requires no spatial grid, it avoids grid-related numerical artifacts and lets modelers specify inputs directly in terms of real features and boundary conditions, striking a balance between the simplicity of closed-form analytical methods and the flexibility of full numerical models. The second ingredient is Monte Carlo simulation, a decades-old statistical technique that propagates uncertainty by repeatedly sampling input parameters from probability distributions and running the model anew for each sample. Coupling the two produces not a single capture zone but an entire ensemble of equally plausible ones, each reflecting a different combination of aquifer properties.

Technically, the researchers built their workflow around TimML, an open-source Python and Fortran AEM package created by Mark Bakker, whose open code made it possible to add a dedicated Monte Carlo component. The implementation is delivered as three Jupyter notebooks covering the full pipeline: a pre-processor that defines uncertain inputs as probability distributions, a sampling engine that runs TimML repeatedly with different parameter sets, and a post-processor that visualizes the results in a GIS environment. Three parameters were treated as uncertain because they are typically poorly constrained yet strongly influence travel-time-based capture zones: hydraulic conductivity, saturated aquifer thickness, and effective porosity. Hydraulic conductivity was represented as lognormally distributed, a standard choice in hydrogeology because the parameter is strictly positive and often spans orders of magnitude; the authors verified this assumption with quantile-quantile plots and a Shapiro-Wilk test on the log-transformed data from 70 Hazen-based estimates. Thickness and porosity were assigned truncated normal distributions bounded by physically plausible limits, with the bounds set at roughly the mean plus or minus three standard deviations.

The treatment of uncertainty is conceptually careful. The authors distinguish aleatory uncertainty, arising from natural spatial variability in the aquifer, from epistemic uncertainty, which stems from limited measurements and imperfect knowledge of site conditions. In the present implementation, the aquifer is represented as a single hydrogeological unit with effective properties, so the Monte Carlo ensemble primarily captures epistemic uncertainty in the effective parameterization. Where conservative estimates were needed, the team computed upper confidence limits of the mean using t-distributions, and for lognormal data applied a Cox-modified method in log space. For each Monte Carlo realization, a unique parameter set is sampled and assigned to the TimML model, reverse particle tracking generates pathlines for the chosen travel times, and the procedure repeats until an ensemble of plausible WHPA realizations accumulates. Latin Hypercube Sampling could reduce the number of runs, but ordinary runtimes proved manageable on available hardware.

The output is summarized through percentile-based envelopes rather than a single boundary. Particle locations along the simulated pathlines are pooled across all realizations, and percentiles of their distances from the pumping well define nested WHPA polygons. Lower percentiles yield larger, more precautionary zones because only points beyond the percentile distance are excluded before polygon construction, while higher percentiles yield tighter, more central delineations. The researchers computed the 50th, 75th, 95th, and 99th percentile envelopes using a convex hull algorithm, exported them as georeferenced GIS layers in the Swedish SWEREF 99 coordinate system, and thereby produced maps that water managers can overlay directly on land-use data. This percentile framework turns an abstract statistical ensemble into something a municipal planner can read, compare, and discuss with stakeholders.

The case study applied the workflow to the Varnum aquifer near Borås in southwest Sweden, an unconfined glaciofluvial delta deposit forming part of the Rångedala esker. The aquifer spans roughly three square kilometers in a valley at about 170 to 175 meters above sea level, with deposit thicknesses between 20 and 55 meters, and consists of fine sand overlying medium to coarse sand on bedrock. Ten years of historical head measurements show a stable water table, justifying a steady-state model under the Dupuit-Forchheimer approximation. Crucially, the team benchmarked their probabilistic results against an existing deterministic numerical model of the site that had been developed, calibrated, and officially approved by the Municipality of Borås. The digitized 100-day pathlines from that numerical model fell squarely within the central portion of the probabilistic ensemble, most consistent with the 50th to 75th percentile envelopes, exactly what one would expect if a deterministic parameter set represents one plausible draw near the center of the assumed distributions.

The magnitude of the uncertainty effect is striking. At the 50th percentile, the 100-day WHPA covered 32,260 square meters, or about 3.23 hectares, concentrated around the well. At the 75th percentile it ballooned to roughly 15.6 hectares, and at the 95th and 99th percentiles it reached 61.7 and 102.1 hectares respectively, a more than thirtyfold expansion across the range. As the envelope widens, it progressively intersects agricultural land, residential areas, cemeteries, and peat bogs, each bringing new potential contaminant sources into scope, from septic systems and household chemicals to leaching from burial grounds and elevated dissolved organic carbon in peatland environments. The authors emphasize that this exposes a genuine trade-off: the safest decision under uncertainty may differ from the statistically optimal one, and choosing a very conservative percentile pulls farmland, forest, and peatland into the protection zone with corresponding management implications, such as fertilizer restrictions. They also caution that final zoning in practice is often adjusted to follow roads and parcel boundaries, which can amplify the effect of the percentile choice on the ultimate protected area.

The researchers are candid about the limitations. Results depend on the probability distributions assigned to the inputs, which at data-poor sites should be read as plausible representations of uncertainty rather than definitive aquifer characterizations, and the parameters were sampled independently, a deliberate simplification that correlated sampling could later relax. The single-unit aquifer representation matches common practice for small supplies, where the realistic alternative is usually a deterministic model under the same homogeneous assumption rather than a heterogeneity-resolving model. Computational cost grows with ensemble size, and location-specific capture probabilities would require additional post-processing. Yet the value proposition is clear: the workflow is a screening-level first step for sites with limited but not absent data, requiring only basic hydrogeological and GIS skills, and integrating preprocessing, modeling, and GIS-ready output in one Python environment reduces the software fragmentation that undermines transparency and reproducibility. As the authors note, risk communication research has long recognized that getting the numbers right is only the beginning; making uncertainty visible, as these percentile envelopes do, allows decision-makers and communities to debate precaution versus land-use burden concretely. For the thousands of small and medium-sized water supplies worldwide that currently rely on a single deterministic boundary despite substantial subsurface uncertainty, this open-source workflow offers a defensible, feasible path toward risk-based groundwater protection.

Subject of Research: Probabilistic delineation of groundwater wellhead protection areas using Monte Carlo simulation and the Analytic Element Method

Article Title: Accessible probabilistic modeling of wellhead protection areas for small and medium-sized water supplies

Article References: Gärtner, N., Zamzami, M., & Lindhe, A. (2026). Accessible probabilistic modeling of wellhead protection areas for small and medium-sized water supplies. Discover Geoscience, 4(1), Article 378. https://doi.org/10.1007/s44288-026-00741-w

Image Credits: AI Generated

DOI: 10.1007/s44288-026-00741-w

Keywords: groundwater, wellhead protection areas, Monte Carlo simulation, analytic element method, hydrogeology, drinking water, uncertainty analysis, open-source software, TimML, water policy, risk assessment, GIS