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:

  1. Solve steady-state Darcy flow to find q
  2. Substitute q into advection-diffusion equation to solve heat equation.
  3. 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.

Sample the posterior - more!

\mathbf{m}_1

Misfit

S(\mathbf{m})

Global
minimum

Local
maximum

Local
minimum

More thoroughly sampling the posterior will identify local min and max

We initialise a MCMC "random walk" at multiple starting points generated from a Latin Hypercube to ensure we are not trapped by local minima.

After N samples the topology of   is understood and gradient-optimisation is employed to quickly descend to

S(\mathbf{m})
\min S(\mathbf{m})

Random Sampling

Latin Hypercube

Non-uniqueness

There may be many solutions that fit the same set of observations.

Data Assimilation

3) Temperature gradients

1) Chloride concentration

2) Hydraulic Pressure Head

Data Assimilation

1) Chloride concentration

Chloride is excluded from evapotranspiration and leaches downward into the water table to become a proxy for recharge.

R = \frac{P \cdot Cl_p}{Cl_{gw}}

Chloride concentration of rainfall

Chloride concentration of groundwater

mean annual rainfall

Crosbie et al. 2018

Data Assimilation

1) Chloride concentration

3) Temperature gradients

2) Hydraulic Pressure Head

Data Assimilation

2) Hydraulic Pressure Head

Groundwater

observation well

Cl recharge observation

layer 1

layer 2

layer 3

Misfit

Data Assimilation

2) Hydraulic Pressure Head

3) Temperature gradients

1) Chloride concentration

Data Assimilation

3) Temperature gradients

Temperature borehole

Misfit

layer 1

layer 2

layer 3

Coal seam

Thermal conductivity

Temperature field with groundwater velocity vectors

Fluid flow off

Fluid flow on

Difference

Data Assimilation

3) Temperature gradients

2) Hydraulic Pressure Head

1) Chloride concentration

High discharge rates off the continental shelf

- up to 3 m/day

 

Very long residence times inland (383 ka) compared to coastal aquifers (182 a)

Groundwater flow velocities for each

geological layer in our 3D model

Misfit with respect to model iterations.

We reach convergence at approximately 8500 iterations.

Parameter adjustment

Future steps

  • Integrate the total catchment area for specific discharge points
  • Calculate the groundwater velocity accounting for increased extraction into the future
  • Identify specific wells that are suitable for groundwater extraction.
  • Scale up workflow to include the Great Artesian Basin

Groundwater velocity (log m/s) beneath major aquitards - i.e. volcanics and coal seams

Water level 2000-2020

Thank you

Dr. Ben Mather

Madsen Building, School of Geosciences,

The University of Sydney, NSW 2006

benmather.info  |  @BenRMather