Bounding the Null Space: Interval-Based Uncertainty Quantification for Non-Identifiable Groundwater Models

This paper proposes an Optimization-based Bound Tightening (OBBT) framework that uses interval arithmetic and McCormick relaxations to provide guaranteed, sampling-free uncertainty bounds for non-identifiable groundwater models, while addressing challenges like non-physical rotational flow through specific sign and irrotationality constraints.

The Gist

The Big Problem: The "Missing Puzzle Pieces"

Imagine you are trying to figure out how water moves underground. You have a few clues: you know where water is being pumped out, where rain falls, and you have measurements of water levels in a couple of wells.

But the ground is huge, and you only have a few measurements. This creates a problem called non-identifiability. It's like trying to guess the exact shape of a hidden object by only touching three of its corners. There are millions of different shapes (combinations of rock types, flow speeds, and water levels) that could fit those three corners perfectly.

The Old Way (Sampling):
Most scientists try to solve this by guessing. They run thousands of computer simulations, each with a slightly different guess about the underground conditions. They look at the results and say, "Okay, the water level is probably between 5 and 10 meters."

• The Flaw: If you only guess 1,000 times, you might miss the real extreme cases. You might think the water level is safe (5–10m), but the actual reality could be 2–15m. You underestimated the danger because you didn't guess enough times.

The New Way: "Bounding the Box" (OBBT)

The authors propose a completely different approach called Optimization-based Bound Tightening (OBBT). Instead of guessing random scenarios, they treat the problem like a math puzzle with strict rules.

The Analogy: The Shrink-Wrap Box
Imagine the possible answers are floating inside a giant, transparent cardboard box.

1. The Initial Box: At first, the box is huge because we don't know much. The water level could be anywhere from 0 to 100 meters.
2. Adding Rules: We then start adding "rules" (constraints) based on physics (water flows downhill) and our actual data (we measured 7 meters here).
3. Shrinking the Box: Every time we add a rule, we can cut away the parts of the box that are impossible. We keep shrinking the box until it fits as tightly as possible around the only answers that are physically possible.
4. The Result: We don't get a list of guesses; we get a guaranteed safe range. We know for a fact the water level cannot be outside this final, tight box.

The Hurdle: The "Broken Compass"

To make this math work on a computer, the authors had to simplify the complex laws of groundwater flow. They used a mathematical trick called McCormick relaxations.

The Analogy:
Think of groundwater flow like a car driving on a road. The car (water) must always drive in the direction the road slopes (downhill).

• The Problem: When the authors simplified the math to make it faster, their "compass" got broken. The math allowed the car to drive uphill if it was a very specific, weird combination of speed and slope.
• The Consequence: Because the math allowed these "impossible" uphill drives, the computer couldn't shrink the box effectively. The box stayed huge because the computer thought, "Well, maybe the water flows uphill here, so I can't rule anything out."

The Fix: Forcing the Rules

The authors realized they had to manually tell the computer, "No, water cannot flow uphill." They added two specific fixes:

1. Flow Signs: They forced the computer to decide early on: "Is the water flowing North or South?" Once that direction is locked in, the "uphill" nonsense disappears.
2. No Swirls: They added a rule that water cannot spin in circles (like a whirlpool) without a reason. This helps the math understand the true shape of the flow.

With these fixes, the "box" finally shrinks tight, giving a reliable answer.

What They Tested

The team tested this method on three different scenarios:

1. A 1D Strip: A simple line of cells. It worked perfectly and was much better than the old "guessing" methods.
2. A 2D Grid: A flat map. This showed that without the "No Swirls" or "Flow Sign" fixes, the method fails. With the fixes, it worked well.
3. A Time-Traveling Grid: A 2D map that changes over time (like a video). They showed the method could handle water levels changing day by day, shrinking the uncertainty as time went on.

The Trade-Off

The Good News: This method gives you guaranteed safety. You don't have to worry about missing a rare, dangerous scenario because you didn't guess enough times. It finds the absolute outer limits of what is possible.

The Bad News: It is computationally expensive. It takes a long time to solve these math puzzles compared to just running a few thousand guesses. It's like using a laser cutter to trim a piece of paper instead of just using scissors. It's slower, but the result is mathematically perfect.

Summary

The paper presents a new way to handle uncertainty in groundwater models. Instead of guessing thousands of times and hoping you catch the worst-case scenario, they use strict math rules to carve away all the impossible answers, leaving behind a guaranteed "safe zone" for the water levels. They found that to make this work, they had to add extra rules to stop the math from imagining impossible water flows, but once they did, the method provided a much more reliable safety net than traditional guessing methods.

Technical Summary

Technical Summary: Bounding the Null Space via Optimization-Based Bound Tightening

1. Problem Statement

Groundwater models are frequently non-identifiable due to sparse subsurface observations. This equifinality means that many distinct combinations of parameters (e.g., transmissivity, specific yield), states (hydraulic heads), and boundary conditions can produce identical observations.

Traditional Uncertainty Quantification (UQ) methods address this by exploring a finite set of model realizations (e.g., Monte Carlo, ensemble data assimilation, or regularized inversion). The authors argue that this paradigm suffers from a structural dilemma: as system complexity increases, the number of unknowns grows, yet the computational cost of simulation limits the number of realizations that can be explored. Consequently, incomplete exploration leads to a systematic underestimation of the true range of admissible solutions. This is particularly dangerous in groundwater management, where risks often arise from rare extremes (e.g., droughts or floods) that finite samples may miss.

The paper posits that for conservative UQ, the goal should not be to sample the feasible region but to represent the extent of the feasible region itself through guaranteed outer bounds on uncertain variables.

2. Methodology: Optimization-Based Bound Tightening (OBBT)

The proposed framework replaces sampling with Optimization-Based Bound Tightening (OBBT). Instead of enumerating realizations, the method treats uncertainty as intervals and tightens them by extremizing variables over a constraint system that encodes physical laws and observations.

2.1 Mathematical Formulation

The approach formulates the groundwater flow problem as a system of constraints:

• Variables: Hydraulic heads ($h$), fluxes ($q$), recharge ($R$), transmissivity ($T$), and specific yield ($S$).
• Constraints:
  • Mass Balance: Discretized using a finite-volume scheme (steady-state and transient).
  • Darcy's Law: $q = -T \nabla h$.
  • Observations: Fixed values or bounds at specific locations.
  • Prior Bounds: Initial intervals for all variables.

2.2 Handling Nonlinearity: McCormick Relaxations

The governing equations involve bilinear terms (e.g., the product of uncertain transmissivity $T$ and head gradient $\nabla h$). Solving the resulting non-convex optimization problems globally is computationally prohibitive for large grids.

To maintain computational tractability while ensuring conservative bounds, the authors employ McCormick relaxations. These replace the bilinear equality $w = xy$ with four linear inequalities that form a convex hull (envelope) around the true nonlinear feasible region.

• Benefit: Allows the use of efficient Linear Programming (LP) solvers.
• Drawback: The relaxation admits non-physical solutions, specifically breaking the sign coupling between fluxes and head gradients. This permits "rotational flow" (cycles where water circulates against the gradient), which violates the physics of Darcy flow and prevents effective bound tightening.

2.3 The Algorithm

The OBBT algorithm (Algorithm 1) operates iteratively:

1. Initialization: Define variables on a spatio-temporal grid graph and set initial bounds.
2. Relaxation Construction: Build McCormick constraints based on current variable bounds.
3. Extremization: For each variable, solve two LPs (minimization and maximization) subject to the full constraint system.
4. Tightening: Update variable bounds with the new extrema.
5. Post-Processing: Use interval arithmetic to enforce sign consistency between fluxes and gradients.
6. Iteration: Repeat until convergence (bounds stop tightening) or a maximum iteration count is reached.

The process is embarrassingly parallel regarding the extremization of individual variables within each iteration.

3. Key Contributions and Remedies

The paper identifies that naive application of OBBT with McCormick relaxations fails to provide informative bounds because the relaxation allows non-physical rotational flow. The authors propose and evaluate two specific remedies to restore physical consistency:

1. Irrotationality Constraints: Explicitly enforcing that the sum of fluxes around any fundamental cycle in the grid graph is zero.

   • Effectiveness: Highly effective in homogeneous media, restoring informative bounds.
   • Limitation: Not generally applicable in heterogeneous or anisotropic media where zero net flux does not strictly follow from irrotationality.

2. Flow Sign Prescription: Pre-specifying the sign (direction) of flow across each interface based on a reference simulation or prior knowledge.

   • Effectiveness: Removes the non-physical quadrants of the McCormick envelope, significantly tightening bounds.
   • Generality: More robust than irrotationality constraints in heterogeneous settings, though it introduces a prior assumption about flow direction.

4. Results

The framework was tested on three numerical examples:

4.1 1D Steady-State Model

• Setup: 10 cells with source/sink wells and two observation points.
• Comparison: OBBT bounds were compared against an exact Monte Carlo sampler (drawing from the true null manifold) and MCMC.
• Findings:
  • OBBT successfully enclosed the null manifold, providing guaranteed outer bounds.
  • MCMC and finite-sample methods severely underestimated the extent of the feasible region (e.g., covering only ~45% of the marginal range even with 100,000 samples).
  • OBBT converged in 2 iterations, demonstrating that sampling is unnecessary for capturing the bounds of this non-identifiable system.

4.2 2D Steady-State Model

• Setup: 5x5 grid with a source, sink, and central observation.
• Findings:
  • Basic OBBT (no extra constraints): Failed to tighten bounds for heads or transmissivity due to the sign ambiguity allowing rotational flow.
  • OBBT with Irrotationality: Successfully resolved flow directions and tightened bounds, closely matching a pure-LP reference solution.
  • OBBT with Flow Signs: Also tightened bounds, though slightly wider than the irrotational case in this homogeneous setup. It demonstrated the trade-off between assumption strength and bound tightness.

4.3 2D Transient Model

• Setup: Hexagonal grid, 5 timesteps, recharge zone, extraction well, and fixed-head boundary.
• Findings:
  • OBBT successfully handled time-varying systems, capturing the development of a drawdown cone and the propagation of constraints backwards in time.
  • Hydraulic head bounds tightened significantly, while transmissivity bounds tightened only slightly (as expected in under-constrained parameter spaces).
  • Computational Cost: The run required ~71.7 hours on 16 cores. The authors note this is substantial but argue it is the "price of conservative coverage" compared to the risk of underestimation in sampling methods.

5. Significance and Claims

The paper claims that OBBT offers a conservative, deterministic, and physically grounded alternative to ensemble-based UQ. Its primary significance lies in its ability to provide guaranteed outer bounds on all uncertain variables without sampling, thereby avoiding the "under-exploration leads to underestimation" problem inherent in Monte Carlo methods.

The authors emphasize that while OBBT is computationally more expensive than current sampling methods, the trade-off is justified for decision-support scenarios where regulators must prove that a proposed abstraction rate is safe under all admissible parameter combinations. The method connects naturally to null space theory, characterizing the "null manifold" (the set of equivalent solutions) through interval bounds rather than point estimates.

The paper concludes that while the current implementation relies on specific assumptions (like flow sign prescription) to handle the non-convexity of Darcy flow, the framework represents a promising avenue for future research, particularly in developing more efficient mixed-integer formulations and integrating with data assimilation workflows.