Bayesian inversion of 3D groundwater flow

within the Sydney-Gunnedah-Bowen Basin

@BenRMather @MullerDietmar @TheCraigONeill @adambeallgeo @wv2017 @LouisMoresi

Ben Mather, Dietmar Muller, Craig O'Neill, Adam Beall, Willem Vervoort, Louis Moresi

The driest inhabited continent

Australian population

Squeeze on surface water

Surface water reserves have become increasingly stretched. With longer droughts, more bushfire days and a growing population, we must seek alternate sources of water.

Groundwater resources

Much is unknown about the dynamics of groundwater flow in Australia, including:

How geology affects local to regional flow patterns
Timescales of groundwater flow through major aquifers
Hydraulic pressure within different layers
Whether saltwater intrusion affects coastal aquifers

Understanding these factors can help us better manage Australia's groundwater resources.

Great Artesian Basin - GA

Coupled groundwater-thermal numerical models of the
Sydney-Gunnedah-Bowen Basin

We have:

Temperature gradients measured in boreholes
Recharge estimates from Cl mass balance
Hydraulic pressure measured in wells
3D geological model of the basin

We want:

Groundwater flow velocities within major aquifers
Capacity of aquifers to store groundwater

Numerical model - coupling heat flow and fluid flow

Advection-diffusion equation:

Coupling fluid flow with heat flow:

Solve steady-state Darcy flow to find q
Substitute q into advection-diffusion equation to solve heat equation.
Repeat step 2 until convergence is reached.

Implemented using Underworld 2

hydraulic = blue

thermal = red

\nabla ( k \nabla T) + H - C \mathbf{q} \cdot \nabla T = 0

\nabla \cdot \mathbf{q} =\nabla \cdot (k_h \nabla h ) = 0

T_0

h_0

\frac{\partial h}{\partial x}=0

\frac{\partial T}{\partial x}=0

T_1

\frac{\partial h}{\partial y}=0

lower BC

upper BC

side BC

c_p \frac{\partial T}{\partial t} = \nabla ( k_T \nabla T ) - C \mathbf{q} \cdot \nabla T + H

Bayes' Theorem

Formally describes the link between observations, model, & prior information.
Where these intersect is called the posterior

P(\mathbf{m}|\mathbf{d}) \propto P(\mathbf{d}|\mathbf{m}) \cdot P(\mathbf{m})

posterior

likelihood

prior

model

data

example of an ill-posed problem

example of a well-posed problem

Data-driven

Use the data to "drive" the parameter exploration.
Infer what input parameters you need to satisfy your data and prior information

Input parameters

Model being solved

Compare data & priors

\mathbf{m} : [\kappa_1,\kappa_2,\kappa_3,\ldots, \kappa_n]

\nabla ( k \nabla T) + H - C \mathbf{q} \cdot \nabla T = 0

P(\mathbf{m}|\mathbf{d}) \propto P(\mathbf{d}|\mathbf{m}) \cdot P(\mathbf{m})

FORWARD MODEL

Prior

P(\mathbf{m})

Likelihood

P(\mathbf{d}|\mathbf{m})

Posterior

P(\mathbf{m}|\mathbf{d})

REPEAT

\nabla \cdot \mathbf{q} =\nabla \cdot (k_h \nabla h ) = 0

Sampling the posterior

The posterior must be properly sampled to avoid local minima e.g. MCMC, Newton method, etc.

\mathbf{m}_1

Misfit

S(\mathbf{m})

Global
minimum

Local
maximum

Local
minimum

P(\mathbf{m}|\mathbf{d}) = A \, \exp (-S(\mathbf{m}))

S(\mathbf{m}) = \frac{1}{2} \sum_{i} \frac{\vert g^i(\mathbf{m}) - d^i \vert^2}{(\sigma^i_d)^2} + \frac{1}{2} \sum_{j} \frac{\vert \mathrm{m}^j - \mathrm{m}^j_p \vert^2}{(\sigma^j_p)^2}

A "misfit" function is chosen to find the minimum misfit between data and priors.