Groundwater numerical modelling

The focus of the numerical modelling for the Clarence-Moreton Bioregional Assessment is the coal seam gas (CSG) development in the Walloon Coal Measures in the Richmond river basin (Figure 3). The Richmond river basin overlies a sedimentary basin with an alternating sequence of aquifers and aquitards (Figure 4) in which the alluvial deposits of the Richmond River and its tributaries are incised. The Walloon Coal Measures outcrop in the west of the study area. In the proposed CSG development area close to Casino, the Walloon Coal Measures are overlain by hundreds of metres of younger sedimentary bedrock, including several hydrostratigraphic units that are considered to be aquitards (Raiber et al., 2016). In the northern part of the Richmond river basin, the Lamington Volcanics form a topographic high, and these high-elevation basalts are associated with very high groundwater recharge rates.

The Richmond River system is perennial, with the exception of some smaller tributaries in the west such as Shannon Brook. Unlike the perennial streams, these smaller tributaries in the western part of the Richmond river basin are not incised into the Lamington Volcanics, and flow in Shannon Brook, for instance, is observed 94% of the year. Downstream from Casino, the Richmond River is tidally influenced.

Figure 3

Figure 3 Three-dimensional geological model of the Clarence-Moreton bioregion

Viewed from the south-east; the vertical extent is from –2500 to +1400 m Australian Height Datum (AHD); the north–south extent is 320 km; the maximum east–west extent is 140 km; the vertical exaggeration is 10

Source: Raiber et al. (2016)

The main causal pathway group identified in companion product 2.3 for the Clarence-Moreton bioregion (Raiber et al., 2016) considered in the numerical groundwater modelling is the ‘Subsurface depressurisation and dewatering’ causal pathway group. Dewatering applies to coal mines and, given that no coal mines are being modelled in the baseline or CRDP, it will not be discussed further in relation to the Clarence-Moreton bioregion. Subsurface depressurisation of the Walloon Coal Measures near Casino, may lead to drawdowns in the Walloon Coal Measures and overlying aquifers including the alluvial aquifers. A change in groundwater pressures and levels in the alluvium may result in changes to surface water – groundwater interactions and thus a reduction of streamflow in the Richmond River or its tributaries.

Figure 4

Figure 4 Fence diagram through the Richmond river basin showing geometric and thickness relationships between alluvial, volcanic and sedimentary bedrock hydrostratigraphic units

Different model types and model codes have been adopted in bioregional assessments (BAs) depending on the specific requirements of each subregion or bioregion. The companion submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2016) stipulates that a groundwater model must deliver spatially explicit outputs that can evaluate the direct/indirect impacts on groundwater resources as well as providing input to other BA modelling tasks, including surface water modelling, uncertainty analysis and receptor impact modelling. A formal assessment of the modelling approach adopted in the Clarence-Moreton bioregion was conducted according to the five fit-for-purpose criteria listed in Table 3.

Table 3 Assessment of groundwater numerical modelling approach in the Clarence-Moreton bioregion

Fit-for-purpose assessment criteria


Criteria met

1. Prediction of hydrological response variables

Probabilistic estimates of hydrological change at mode nodes


Integration with receptor impact modelling


Linkage with surface water numerical models


2. Design and construction

Modelling objectives stated


Model confidence level


Modelling approach


3. Integration with sensitivity and uncertainty analysis workflow

Formally address uncertainty






4. Water balance components

Conceptual model agreement


5. Transparent and reproducible model outputs

Model data repository


Model code and executables


Pre- and post- processing scripts


Transparent description of conceptual model and traceable process for deriving model parameter/data

Yes Prediction of hydrological response variables

The objective of groundwater modelling that is undertaken as part of a BA is to probabilistically assess the impacts of additional coal resource development at model nodes associated with water-dependent assets. The impact is quantified by the change in some predefined hydrological response variables such as fluxes, hydraulic heads or variables depending on fluxes or hydraulic heads. Hydrological response variables need to be decided upon prior to the sensitivity analysis, for example, drawdown at model node (x, y, z) at time t. The hydrological response variables adopted in the current assessment include additional drawdown (dmax, maximum difference in drawdown between the coal resource development pathway (CRDP) and baseline, due to the additional coal resource development), year of maximum change (tmax) and leakage/baseflow changes at a number of model nodes.

In order to quantify the uncertainties associated with model predictions, groundwater modelling in BA is run probabilistically and not deterministically. Consequently, this means that model predictions cannot be presented as unique values but as probability distributions. Mismatches in spatial and temporal scale between the regional nature of the modelling and the point-scale nature of model nodes mean that the modelling is not able to capture fine-scale complexities of impacts upon model nodes. In BA, the focus is on the impact of the additional coal resource development, which is attributed to the difference in results between two possible futures – the baseline and the CRDP. Hence, numerical groundwater modelling results are primarily presented as the difference in groundwater drawdown between the baseline and CRDP. This approach reduces the uncertainty as the prediction is formulated as the difference between the output of two model results (Barnett et al., 2012).

The linkage between the groundwater model and the Australian Water Resources Assessment (AWRA) landscape model (AWRA-L) (Vaze et al., 2015) is depicted in Figure 5. AWRA-L was used to generate the temporal trends of diffuse recharge for the groundwater model. The river stage elevation time series for the groundwater model was derived from the runoff outputs of AWRA-L, which were coupled with rating curves derived from measured data sourced from the NSW streamflow dataset from NSW Office of Water. The groundwater model was used to generate the leakage/baseflow change at a number of gauges along the Richmond River due to the additional coal resource development. These flux changes can then be fed back into the AWRA-L model predictions to estimate changes of relevant hydrological response variables. Subsequently, changes of hydrological response variables can be used to assess the impacts on model nodes associated with surface water-dependent assets.

Although the groundwater model was able to deliver input for receptor impact modelling, receptor impact modelling was not considered in this current BA. Analysis of the potential impacts of the additional coal resource development have shown minimal hydrological changes, particularly at the uppermost aquifer where median additional drawdowns does not exceed 0.2 m. The median predicted change in surface water – groundwater flux at simulated model nodes, including both leakage and baseflow, does not exceed 0.01 GL/year as shown in companion product 2.6.1 for the Clarence-Moreton bioregion (Gilfedder et al., 2016). As a result of these minor hydrological impacts, ecological impacts are expected to be below the detectable limit and hence the decision was to not carry out receptor impact modelling in the Clarence-Moreton bioregion.

Figure 5

Figure 5 Model sequence for the Clarence-Moreton bioregion

AWRA-L = Australian Water Resources Assessment landscape module; baseline = baseline coal resource development; CRDP = coal resource development pathway; CRDP = baseline + additional coal resource development; green rectangles represent models; orange rectangles are input and/or output data of models Design and construction

Following the guidelines of companion submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2016), a numerical groundwater model was developed to assess the potential impacts of the additional coal resource development on the groundwater system of the Clarence-Moreton bioregion. As reported in companion product 2.3 for the Clarence-Moreton bioregion (Raiber et al., 2016), the additional coal resource development in this BA only includes one CSG development, the West Casino Gas Project. The groundwater model mainly focused on the Richmond river basin where the West Casino Gas Project is likely to occur. The conceptualisation and the likely causal pathways that underpinned the model are identified in companion product 2.3 for the Clarence-Moreton bioregion (Raiber et al., 2016).

The major assumptions and model choices related to model development are described in Section The groundwater model is considered the first regional model for the area and is based on a three-dimensional geological model that was specifically developed for the current BA (see companion product 2.3 for the Clarence-Moreton bioregion (Raiber et al., 2016)). MODFLOW-NWT (Niswonger et al., 2011) was used to run the model in order to minimise the numerical instabilities that can be caused by cell drying and rewetting in the unconfined layer of the model. The effect of the assumptions and model choices on predictions are discussed in Section

Groundwater modelling simulations were run under historical conditions for the 30-year period from 1 January 1983 to 31 December 2012. The historical climate time series was repeated three times to create a 90-year time series and modified to be consistent with a median future climate projection as described in companion submethodology M06 (as listed in Table 1) for surface water modelling (Viney, 2016). This 30-year period was repeated to ensure that the effect of droughts and floods does not confound the comparison between time periods. Monthly stress periods were used throughout the simulation to provide monthly change of baseflow required by the surface water model. This led to a total number of 1441 stress periods from 1983 to 2102, including the first steady-state stress period.

Final outputs of predicted impacts were produced with the aid of emulators as described in companion submethodology M07 (as listed in Table 1) for groundwater modelling (Crosbie et al., 2016) and M09 (as listed in Table 1) for propagation of uncertainty through models (Peeters et al., 2016). The main goal of the groundwater model was to provide realisations that can be used to train the emulators for the model nodes. Emulators are computationally efficient approximation statistical tools compared to the slower process-based groundwater model. As described in companion submethodology M09 (Peeters et al., 2016), the application of model emulation is to overcome the computational inefficiency of the original groundwater model, but still being able to capture the important groundwater flow features.

For those model nodes that are associated with groundwater-dependent assets and hence directly impacted by groundwater extraction, such as stock and domestic groundwater bores, the groundwater model directly generates the drawdown impacts. If receptor impact modelling is deemed necessary, then the groundwater model outputs will directly feed into the receptor impact model. For other model nodes that are indirectly impacted by groundwater extraction through an interaction with any of the various surface water features, a surface water model is required to simulate the impact of the additional coal resource development. Due to minimal regulation in the Richmond river basin, the rainfall-runoff model AWRA-L (Hafeez et al., 2015) was deemed adequate to model the Richmond River flows. The modelling results for these model nodes are presented in companion product 2.6.1 for the Clarence-Moreton bioregion (Gilfedder et al., 2016). Integration with sensitivity and uncertainty analysis workflow

Figure 6

Figure 6 Uncertainty analysis workflow

ABC MCMC = Approximate Bayesian Computation 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 model node
    2. objective function tailored to each hydrological response variable at each model node
  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 model node.

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 are carried out, sampling extensively from this parameter range. For each evaluation the corresponding simulated change in hydrological response variables at the model nodes 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 model node.

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 into 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 two models in the model chain for the Clarence-Moreton bioregion 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 and Section provide details of the implementation of this uncertainty propagation workflow for the MODFLOW model. The uncertainty analysis for AWRA-L is in Section and Section of companion product 2.6.1 for the Clarence-Moreton bioregion (Gilfedder et al., 2016). 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. Water balance components

The water balance was reported for a defined control volume in BA that includes all hydrologically connected changes predicted by the surface water and groundwater models. The water balance components (e.g. recharge, evapotranspiration, baseflow, licensed extractions, upward flow from deeper groundwater and change in storage) were compared with estimates described in the regional-scale conceptual model and localised groundwater models to provide confidence in model predictions, thus meeting the modelling guiding principles (Barnett et al., 2012). Companion product 2.5 for the Clarence-Moreton bioregion (Cui et al., 2016) compares model estimates of water extraction from CSG with available local and regional estimates. Parameters associated with water balance components are discussed in Section and Section Transparent and reproducible model outputs

The BA requirement that model results must be reproducible means that the models need to be run as part of a documented workflow that records the provenance of the input data, executables and outputs. This is achieved through the use of scripting in BAs. All pre-processing, model runs and post-processing is done using scripts that are made available along with the products at www.bioregionalassessments.gov.au; this ensures that all model inputs, parameters, executables and outputs are traceable, meeting the modelling guiding principles related to transparency (Barnett et al., 2012).

The current model described in this product is designed specifically for delivering a probabilistic assessment of the impacts of the additional coal resource development in the Clarence-Moreton bioregion. In its current form, the model is neither suited to address any other water management questions nor provide deterministic predictions in this bioregion. The validity of all assumptions in this product should be re-assessed for any other application of the current model.

Last updated:
23 October 2018
Thumbnail images of the Clarence-Moreton bioregion

Product Finalisation date