### Ben Mather PRO

Computational Geophysicist at the Sydney Informatics Hub, University of Sydney

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

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

Australian population

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

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

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

- 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

- 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

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.

\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**

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

**3) Temperature gradients**

**1) Chloride concentration**

**2) Hydraulic Pressure Head**

**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

**1) Chloride concentration**

**3) Temperature gradients**

**2) Hydraulic Pressure Head**

**2) Hydraulic Pressure Head**

Groundwater

observation well

Cl recharge observation

layer 1

layer 2

layer 3

Misfit

**2) Hydraulic Pressure Head**

**3) Temperature gradients**

**1) Chloride concentration**

**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**

**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.

**Dr. Ben Mather**

Madsen Building, School of Geosciences,

The University of Sydney, NSW 2006

By Ben Mather

AESC 2021 conference

- 735

Computational Geophysicist at the Sydney Informatics Hub, University of Sydney