2.6.2.1.2 Groundwater numerical modelling


Figure 3

Figure 3 Conceptual block diagram of the Gloucester subregion

The conceptual understanding of the Gloucester subregion is summarised in Figure 3 and companion product 2.3 for the Gloucester subregion (Dawes et al., 2018). The Gloucester Basin is considered to be a geologically closed basin with three main hydrogeological units:

  • surface alluvium up to 15 m thick, a semi-confined to unconfined aquifer
  • shallow weathered and fractured rocks up to 150 m thick, a confined to semi-confined aquifer
  • interburden units alternating with coal seams to a maximum depth of about 2500 m, only considered to be water-bearing strata.

The shallow weathered and fractured rock layer (SRL) underlies the alluvium entirely, and outcrops extensively across the rest of the surface of the Gloucester subregion. Both the Avon and Karuah rivers are unregulated streams (see Section 2.1.6 in companion product 2.1-2.2 for the Gloucester subregion (Frery et al., 2018)), connected with local groundwater (McVicar et al., 2014, Section 1.1.6). The river system is mostly gaining, with baseflow estimated to be about one-tenth of total streamflow (Dawes et al., 2018, Section 2.3.2.4). The alluvial aquifer only receives water from the river system during high flow and flood events.

The main causal pathways inferred in Dawes et al. (2018), Section 2.3.5.3, for coal seam gas operations are aquifer depressurisation and inter-aquifer connectivity, while for open-cut mines the disruption of natural surface water drainage and inter-aquifer connectivity, including pumping to dewater mines by lowering the watertable, are listed as main causal pathways. Any water extraction in either the coal seams, interburden or the shallow weathered and fractured rocks has the potential to affect groundwater levels in the SRL and alluvium. A change in groundwater levels in the alluvium may affect the surface water – groundwater exchange flux and thus streamflow. Open-cut coal mines have a more direct impact on surface water flow as all rainfall within the mine footprint area is contained on site and no longer contributes to runoff. The role of faults and fractures to increase or decrease inter-aquifer connectivity is highlighted as an important knowledge gap.

As outlined in companion submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2016), different model types and model codes are chosen in a bioregional assessment (BA), depending on the specific requirements of each subregion. The main goal of each groundwater model in BA remains, however, to deliver spatially explicit model outputs that are used as inputs to other BA models, including surface water modelling, uncertainty analysis and receptor impact modelling, and to directly evaluate change to water resources. Table 3 lists the criteria a groundwater model in BA needs to satisfy to be considered fit for purpose for BA. Beneath the table, these fit-for-purpose criteria are discussed briefly for the numerical modelling approach taken in the Gloucester subregion. The remainder of this product describes in greater detail the numerical modelling, the underlying assumptions and their effect on predictions.

Table 3 Assessment of groundwater numerical modelling approach in the Gloucester subregion


Fit-for-purpose assessment criteria

Components

Prediction of hydrological response variables

Probabilistic estimates of hydrological change at receptors

Integration with receptor impact modelling

Integration with surface water numerical models

Design and construction

Modelling objectives stated

Model confidence level

Modelling approach

Integration with sensitivity and uncertainty analysis workflow

Formally address uncertainty

Parameterisation

Convergence

Water balance components

Conceptual model agreement

Transparent and reproducible model outputs

Model data repository

Model code and executables

Pre- and post-processing scripts

2.6.2.1.2.1 Prediction of hydrological response variables

The objective of the numerical modelling undertaken as part of a BA is to probabilistically assess hydrological changes arising from coal resource development at water-dependent assets and receptors (see companion submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2015)). The groundwater and surface water modelling predicts changes in hydrological response variables, the hydrological characteristics of the system or landscape class that potentially change due to coal resource development. These hydrological response variables are the input for receptor impact models that will evaluate how the change in hydrology and hydrogeology results in a change in the economic value or ecology of assets.

In order to probabilistically assess a change through the uncertainty analysis, the receptor locations and hydrological response variables need to be defined explicitly, both in their spatial and temporal support as in what state variable of the numerical model they relate to. The receptors directly linked to groundwater in the receptor register (see companion product 1.3 (McVicar et al., 2015) for the Gloucester subregion) are point locations associated with bores, basic water rights and groundwater (stock and domestic bores), water access rights and supply and monitoring infrastructure, were provided by the NSW Office of Water (2013). Groundwater-dependent ecosystems for which spatial coordinates are available and, at a minimum which aquifer or hydrogeological unit they are associated with, also were identified by the NSW Office of Water (Kuginis et al., In prep.). Table 4 summarises these groundwater receptors and their locations are shown in Figure 4.

Table 4 Summary of groundwater receptors in the Gloucester subregion


Zone

Basic water right (stock and domestic)

Water access right

Supply and monitoring infrastructure

Groundwater-dependent ecosystem

Total

Surface weathered and fractured rock layer

20

16

86

0

122

Alluvium

1

1

28

32

62

Outside model domain

1

1

14

39 (*)

55

Within mine lease

0

1

23

0

24

Total

22

19

151

71

263

(*) outside alluvium

Data: Bioregional Assessment Programme (Dataset 1)

The four types of receptors (basic water right, water access right, supply and monitoring, and groundwater-dependent ecosystem) are assigned to the zone they are located in. Bores are assigned to the surface weathered and fractured rock layer if they are outside the alluvium or, if they are in the alluvium, have a reported bore depth in excess of 15 m. Bores situated in the alluvium with no recorded bore depth or depth less than 15 m are assigned to the alluvium. 90 out of the 122 bores assigned to the surface weathered and fractured rock layer have recorded bore depths. Only 15 of these have depths in excess of 150 m, which are limited to supply and monitoring infrastructure bores. While it is very likely that these are bores associated with the exploration and monitoring of coal seam gas, the receptor database does however not record the owner of each bore. 24 groundwater receptors, 1 water access right and 23 supply and monitoring bores, appear to be located within a mine lease. These are excluded from the numerical modelling as they are beyond the resolution of the regional model and out of scope for the BAs.

Figure 4

Figure 4 Groundwater receptor locations in the Gloucester subregion

ACRD = additional coal resource development

Data: Bioregional Assessment Programme (Dataset 1, Dataset 2, Dataset 3)

Groundwater-dependent ecosystems are only included in the list of predictions of the numerical modelling if they are situated in the alluvium. The regional watertable outside the alluvium is hosted in the surface weathered and fractured rock layer and is mostly found at depths greater than 5 to 10 m (Parsons Brinckerhoff, 2015, Fig. 7.9). Groundwater-dependent ecosystems outside the alluvium are therefore likely to depend on local groundwater flow conditions beyond the resolution of the regional-scale groundwater modelling. For this reason, groundwater-dependent ecosystems outside the alluvium are not formally included as receptors in the numerical modelling. The regional-scale hydrological change resulting from the numerical modelling however, can be combined with local hydrogeological information to assess the change at the groundwater-dependent ecosystem location.

The components of the water balance are not specified as a hydrological response variable of relevance for a receptor impact model and are therefore not reported in this product. Companion product 2.5 for the Gloucester subregion (Herron et al., 2018) does summarise the water balance for the subregion and the change in these balances for selected periods in the future.

The hydrological response variables for groundwater are maximum difference in drawdown (dmax) for one realisation within an ensemble of groundwater modelling runs, obtained by choosing the maximum of the time series of differences between two futures, and year of maximum change (tmax) at receptor locations, where drawdown is defined as the difference in groundwater level between the baseline and coal resource development. The difference in drawdown between CRDP and baseline is due to additional coal resource development.

For surface water, nine hydrological response variables are defined in submethodology M06 (as listed in Table 1) for surface water modelling (Viney, 2016) at 30 nodes along the stream network (companion product 2.6.1 for the Gloucester subregion (Zhang et al., 2018)).

Simulating the change in hydrological response variables at the various receptor locations necessitates the development of an integrated surface water – groundwater model. Groundwater and surface water, however, operate at very different spatial and temporal scales. The surface water obviously is bound to the river channel and floodplain. Streamflow is very responsive to individual rainfall events, requiring at least a daily temporal resolution. The groundwater in the alluvium is bound to the alluvial sedimentary deposits, which form a strip along the rivers of about 15 m thick. Groundwater levels in the alluvial respond to changes in rainfall and river stage, albeit at a longer timescale than surface water (see Section 2.1.5.2 in Frery et al. (2018)). Capturing this dynamic in a numerical model necessitates at minimum a monthly resolution. The deeper hydrogeological units, the SRL, interburden and coal seams are much more spatially extensive, both horizontally and vertically. The groundwater dynamics are very slow, with limited and delayed response to recharge events or flood events in the shallow weathered and fractured rock layer (Parsons Brinckerhoff, 2015, pp. 44–45). To simulate this part of the groundwater system, a high temporal resolution is not required.

While fully coupled surface water – groundwater model codes are available (e.g. HydroGeoSphere, Brunner and Simmons, 2012), their use is not justified within BA due to the high data requirements for parameterisation and due to operational constraints. The latter relates mainly to the general numerical instability of such models and long runtimes which would severely limit a probabilistic uncertainty analysis that requires the models to be evaluated 100s of times with vastly different parameter sets.

For this Assessment, a pragmatic coupling of three models was developed, consisting of a regional groundwater model and an alluvial groundwater model to simulate the impact on the groundwater systems, and a rainfall-runoff model to simulate the impact on the surface water systems of the subregion (Figure 5). The individual models have different spatial and temporal resolution which requires a set of customised processing steps to up or downscale model data to allow the models to be linked.

Figure 5

Figure 5 Model sequence for the Gloucester subregion

GW AEM: regional analytic element groundwater model; GW ALV: alluvial MODFLOW groundwater model; AWRA-L: rainfall-runoff model; SRL: surface weathered and fractured rock layer; dmax: maximum difference in drawdown for one realisation within an ensemble of groundwater modelling runs, obtained by choosing the maximum of the time series of differences between two futures; tmax: year of maximum change; Δh: change in groundwater level; ΔQb: change in surface water – groundwater interaction flux; Qt: total streamflow; ΔHRV: change in hydrological response variable; CRDP: coal resource development pathway; SW: surface water

The regional groundwater model is an analytic element model (referred to as GW AEM), designed to simulate the change in drawdown at the receptors associated with the groundwater bores in the Gloucester geological basin weathered zone, and to provide the change in groundwater level underneath the Avon and Karuah alluvium. The latter provides the lower boundary condition for the alluvial groundwater models. For both alluvial systems a MODFLOW model was developed (referred to as GW ALV) to simulate the change in drawdown on receptors associated with the alluvium and the change in surface water – groundwater flux. This flux is taken into account in the Australian Water Resources Assessment landscape (AWRA-L) surface water model generated streamflow. The change in a number of hydrological response variables is modelled at surface water receptor locations. The modelling of river management or routing of streamflow through the river network with a river model is not necessary as the salient features of streamflow can be simulated solely with a rainfall-runoff model (see companion submethodology M06 (as listed in Table 1) for surface water modelling (Viney, 2016)).

Figure 5 shows in more detail the sequencing of the different models. In the GW AE baseline coal resource development (baseline) model, the impact of historical coal mines and coal mines commercially producing coal as of December 2012 are simulated. The GW AE CRDP model simulates the impact of the CRDP, which is the impact of the baseline coal resource developments as well as those that are expected to begin commercial production after 2012. The difference in simulated drawdown of those two runs will be the simulated impact of the CRDP on the economic receptors in the shallow weathered and fractured rock layer of the Gloucester geological basin.

The GW AE baseline model and GW AE CRDP model simulated impacts underneath the alluvium feed into the alluvial groundwater models for the Avon and Karuah rivers. The difference in simulated drawdown of those two model runs is the simulated impact of the CRDP on the economic and ecological receptors in the Avon and Karuah alluvium. The GW ALV models for the Avon and Karuah rivers also simulate the time series of the change in surface water – groundwater exchange flux, increment Q b left parenthesis t right parenthesis, for the surface water catchments associated with receptor nodes in the AWRA-L model as:

increment Q b left parenthesis t right parenthesis equals Q b subscript baseline end subscript left parenthesis t right parenthesis minus Q b subscript CRDP end subscript left parenthesis t right parenthesis

(1)

The AWRA-L baseline run simulates streamflow at surface water receptors incorporating the effect of existing and approved open-cut coal mines. The AWRA-L CRDP run simulates streamflow at the surface water receptors incorporating the effect of existing and approved open-cut coal mines plus the additional coal resource development. The total streamflow difference, increment Q t left parenthesis t right parenthesis is obtained as:

increment Q t left parenthesis t right parenthesis equals Q t subscript baseline end subscript left parenthesis t right parenthesis minus Q t subscript CRDP end subscript left parenthesis t right parenthesis minus increment Q b left parenthesis t right parenthesis

(2)

The time series of increment Q t left parenthesis t right parenthesis are summarised in the nine hydrological response variables to highlight different aspects of the hydrograph. These hydrological response variables will inform the receptor impact models for the surface water receptors.

2.6.2.1.2.2 Design and construction

According to the Australian groundwater modelling guidelines (Barnett et al., 2012), it is essential to design and construct the groundwater model in function of clearly stated objectives and to provide a model confidence level. The objective of the modelling is explicitly stated in the previous section. The model confidence level is an a priori categorisation of a groundwater model to reflect its predictive capability in function of the model complexity, prediction timeframe and data availability. As clarified in submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2016), the groundwater models in the BAs are all classified as level 1, the lowest level, as they are required to make predictions of unprecedented stresses at time frames longer than periods with data available to constrain the model.

The objectives of the numerical modelling are not to simulate the state of groundwater in the future under baseline and coal resource development conditions, but to quantify the difference between those two futures. This is a very important nuance to the modelling objectives as it allows to make a number of simplifying assumptions based on the principle of superposition (Reilly et al., 1987). The principle of superposition means that for linear systems, the solution to a problem involving multiple inputs (or stresses) is equal to the sum of the solutions to a set of simpler individual problems that form the composite problem. To simulate the effect of change in stress, such as depressurisation and dewatering for coal resource development, it is therefore sufficient to only know the change in stress. It is not necessary to know the initial conditions in the aquifer or the other fluxes and stresses, provided these do not change (Barlow and Leake, 2012). The principle of superposition underpins most of the pumping test interpretations (Kruseman and de Ridder, 1994); aquifer parameters are inferred from the change in stress (pumping rate) and change in groundwater level (drawdown).

The principle of superposition is only valid for linear systems, that is, systems where the response to a change in stress is proportional to the change in the stress. In other words, where a doubling of stress will result in a doubling of the response. In groundwater flow dynamics this condition is satisfied for confined aquifers. Unconfined aquifers are not strictly linear, as the transmissivity depends on the saturated thickness. Reilly et al. (1987) and Rassam et al. (2004) do show, however, that the concepts are still valid for mild violations of the linearity conditions.

As such, this concept can be implemented in any groundwater modelling code. The choice for analytic element modelling for this regional model is driven by the capability of this code to simulate flow through and the effect of linear features. This is essential to represent flow through faults and fractures. An example of using an analytic element model to represent flow through faults is presented in Department of the Environment (2014, p. 25, Figure 2.4). Simulating this behaviour in finite difference groundwater models, such as MODFLOW-2005, is more challenging. The current analytic element software is very flexible and makes stochastically varying crucial aspects of the conceptualisation, such as the number, position and nature of faults or the number and position of coal seams, relatively straightforward. This means that in the uncertainty analysis it becomes possible to explore more of the conceptual model uncertainty than would be practically feasible with finite difference model codes.

As the surface weathered and fractured rock layer in the Gloucester Basin can be considered as a largely confined flow system (McVicar et al., 2014, p. 56 in Section 1.1.4), the principle of superposition can be applied to directly estimate the change in groundwater levels and fluxes due to coal resource development with the analytic element model. This model approach is very similar to the approach used in Leake et al. (2008). The alluvial system is unconfined and the linearity assumption is likely to be violated. This warranted the design and construction of local MODFLOW models for the Avon and Karuah systems, with the change in boundary conditions due to coal resource development provided by the analytic element model. The MODFLOW models can be considered local refinements or child models of the analytic element model, similar to the approach taken in Abrams et al. (2016).

Further technical detail of the conceptualisation, parameterisation and implementation are documented in Section 2.6.1 of companion product for the Gloucester subregion (Zhang et al., 2018) for the AWRA-L model, and Section 2.6.2.3 for the GW AEM and alluvial MODFLOW models.

2.6.2.1.2.3 Integration with sensitivity and uncertainty analysis workflow

Figure 6

Figure 6 Uncertainty analysis workflow

ABC MCMC= Approximate Bayesian Computating Markov chain Monte Carlo; HRV = hydrological response variable

Submethodology M09 (as listed in Table 1) (Peeters et al., 2016) discusses in detail the propagation of uncertainty through numerical models in the BAs. Figure 6 summarises the uncertainty propagation workflow which consists of four major steps:

  1. Design of experiment: large number of model chain evaluations with a wide range of parameter values
  2. Train emulators for:
    1. each hydrological response variable at each receptor location
    2. objective function tailored to each hydrological response variable at each receptor location
  3. Create posterior parameter probability distribution through Approximate Bayesian Computing Markov chain Monte Carlo
  4. Sample the posterior parameter probability distribution to generate the posterior probability distribution for each hydrological response variable at each receptor location.

The first step is to identify the parameters of the model chain to include in the uncertainty analysis and to define a wide range that represents the plausible range of the parameters. A large number of model chain evaluations is carried out, sampling extensively from this parameter range. For each evaluation the corresponding simulated change in hydrological response variables at the receptor locations is stored, together with the simulated equivalents to the observations. The latter are summarised into objective functions, tailored to each hydrological response variable.

This information forms the basis for the subsequent uncertainty analysis. In the uncertainty analysis, the prior parameter distributions, the most likely range of the parameter values based on data and expert knowledge, is constrained with the available relevant data using the Approximate Bayesian Computation methodology. This results in a posterior parameter distribution, tailored to a specific hydrological response variable, which subsequently can be sampled to generate a probability distribution at each receptor location.

This type of uncertainty analysis requires a very large number of model evaluations, which is practically not feasible. This is the main reason that the original model chain in the uncertainty analysis is replaced by emulators, statistical functions that closely mimic the effect of parameter values on predictions. These emulators take little time to evaluate and are straightforward to integrate in the uncertainty analysis workflow.

In order for the model chain to be amenable for incorporation into this uncertainty analysis it needs to be scripted so that parameter values can be changed in an automated fashion, be able to be evaluated from a command line on high performance computers and, most importantly, be numerically stable so that the model converges for a wide range of parameter values.

The three models in the model chain for the Gloucester subregion have text files as input files and can be executed from the command line. The robustness of each model is tested through a stress-test in which a selection of extreme parameter combinations is evaluated. While this does not guarantee that all model evaluations will converge, it provides confidence that the majority of parameter combinations will.

Section 2.6.2.7 and Section 2.6.2.8 provides the details of the implementation of this uncertainty propagation workflow for the GW AEM and alluvial MODFLOW models. The uncertainty analysis for the AWRA-L model is in sections 2.6.1.5 and 2.6.1.6 of companion product 2.6.1 for the Gloucester subregion (Zhang et al., 2018). These sections also have a qualitative uncertainty analysis that provides a structured discussion of the assumptions and model choices not included in the numerical uncertainty analysis and the perceived effect on the predictions.

2.6.2.1.2.4 Water balance components

A secondary objective of the numerical models is to inform the water balance assessment (companion product 2.5 for the Gloucester subregion (Herron et al., 2018)). The MODFLOW models and AWRA model produce estimates of the water balance under baseline and coal resource development futures and can therefore be used in that assessment. The analytic element model, however, only simulates the change in stress due to coal resource development. Its model output therefore has no information on other components of the regional water balance such as recharge or lateral exchange fluxes.

Notwithstanding this limitation, the simulated coal seam gas water production rates will be compared to the estimates of water production by the main coal seam gas proponent to formally constrain the parameter values of the analytic element model. Similarly, estimates of historical groundwater fluxes in the Avon and Karuah alluvial systems will be used to constrain the MODFLOW models.

2.6.2.1.2.5 Transparent and reproducible model outputs

An over-arching requirement of the BAs is for all model outputs to be transparent and reproducible.

Input data, model files (including the pre- and post-processing scripts and executables) and results are available at www.bioregionalassessments.gov.au.

As the evaluation of the model chain is a highly automated and scripted process, it is possible to reproduce the results reported in this product using these scripts and executables, provided the computational resources are available.

Last updated:
7 January 2019
Thumbnail of the Gloucester subregion

Product Finalisation date

2018
PRODUCT CONTENTS