There is a paradox at work when simulating groundwater systems. Our models solve complex partial differential equations which ensure that groundwater flow and contaminant transport rigorously obey the laws of physics. However the coefficients that appear in these equations are only vaguely known. These represent the hydraulic properties of the media through which groundwater flows. Is there any other field of engineering or science where such exactness is applied to media whose properties are so poorly known?

It is not just the properties of the subsurface that are poorly known. Many of the processes that affect the movement of groundwater, and the solutes that it carries, are also poorly known. This applies particularly to recharge processes. Complex temporal and spatial variability of recharge characterise most groundwater systems.

Numerical simulation cannot capture all of these details correctly. Furthermore, the greater the detail that it attempts to represent, the less accurately can it represent it. Little is known about the details of the details.

Consider for a moment the domain of a typical groundwater model.

Many groundwater model domains span areas that are underlain by many different rock types. Flow may be governed by inter-granular porosity in some of these rocks, and by fractures in others. Fractures may be wider near the surface, but may be lined with clay, a by-product of weathering. Shallow alluvial systems may be imposed on this motley association of rocks near rivers and creeks.

Conditions vary rapidly with depth, especially at shallow depths. Weathering status and style vary. Pressure varies. Permeability is sensitive to all of these things.

Geostatistics, as a field of study, recognises that subsurface hydraulic properties are known with certainty at only a handful of locations. So it attempts to characterise them stochastically. If they can be characterised stochastically, then they can populate models stochastically. Models can then make stochastic predictions.

Conceptually, this makes sense.

A problem is that any statistical characterisation of subsurface hydraulic property heterogeneity is itself uncertain. Assumptions such as "stationarity" (that statistical properties are uniform over considerable distances) are generally treated with contempt by geological media whose genesis is an outcome of chaotic environmental processes. To make matters worse, the subsurface features that may exert most influence on groundwater flow (for example faults or old alluvial channels) are often discrete and continuous. At a particular study site, these features may or may not exist. If they exist, their dispositions and connectedness are, on the one hand, of great importance but, on the other hand, are difficult to characterise statistically.

Most groundwater models solve a discretised form of the partial differential equations which groundwater flow is assumed to obey. Discretisation is based on a grid or mesh. It generally also uses a small number of layers to represent media in which significant vertical heterogeneity may prevail.

With discretisation, partial differential equations become large matrix equations. The values of matrix elements may depend on the solutions of these equations. So these matrix equations must be solved iteratively. For large systems, they may be difficult to solve.

Use of a matrix equation in place of a partial differential equation assumes that the same laws that govern flow into and out of tiny representative elementary volumes (REVs) can be applied to model grid cells that encapsulate millions of REVs. It also assumes that we can find hydraulic properties to use at the grid scale that ensure obedience to these equations - regardless of the direction of groundwater flow, and regardless of the heterogeneity of processes and properties that prevail within a single model cell. These grid scale properties cannot be measured. However they can potentially be estimated from properties that are measured at a smaller scale.

The literature (particular the petroleum reservoir literature) documents many studies on upscaling. It concludes that upscaling can only be approximate, and can only be done properly under circumstances where the nature of sub-grid-cell hydraulic property heterogeneity is well understood. This rarely applies to groundwater modelling.

If our goal is to simulate subsurface processes with integrity, we have zero hope.

But is this the goal of decision-support modelling? For some it probably is. This is disappointing. It creates a fertile ground for arguments about whether this aspect or that aspect of a numerical model is correct or incorrect.

Nothing in a numerical groundwater model is correct. Everything is approximate. Furthermore, the asymptote of greater and greater modelling detail (often seen as the path to simulation integrity) is not reality. We know too little about the subsurface to achieve anything like this. Nor is it possible to achieve integrity of "stochastic simulation", wherein everything that we do not know is endowed with "correct", appropriately upscaled, statistical representation.

Pursuit of "integrity of simulation" is a waste of time and money, and a recipe for irresolvable and perpetual argument.

Fortunately, "integrity of simulation" is not a metric by which decision-support groundwater modelling should be judged.  Furthermore, "integrity of simulation" is not something that is necessary to satisfy the metrics that actually apply to decision-support groundwater modelling.