Ben Mather PRO
Computational Geophysicist at the Sydney Informatics Hub, University of Sydney
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.
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
Bayesian inversion of 3D groundwater flow
By Ben Mather
Bayesian inversion of 3D groundwater flow
AESC 2021 conference
