Optimal Design and Prediction-Independent Verification of Groundwater Monitoring Network

Abstract

In this study, we developed a workflow that applies a complex groundwater model for purpose-driven groundwater monitoring network design and uses linear uncertainty analysis to explore the predictive dependencies and provide insights into the veracity of the monitoring design. A numerical groundwater flow model was used in a probabilistic modelling framework for obtaining the spatial distribution of predicted drawdown for a wide range of plausible combination of uncertain parameters pertaining to the deep sedimentary basin and groundwater flow processes. Reduced rank spatial prediction was used to characterize dominant trends in these spatial drawdown patterns using empirical orthogonal functions (EOF). A differential evolution algorithm was used to identify optimal locations for multi-level piezometers for collecting groundwater pressure data to minimize predictive uncertainty in groundwater drawdown. Data-worth analysis helps to explore the veracity of the design by using only the sensitivities of the observations to predictions independent of the absolute values of predictions. A 10-bore monitoring network that collects drawdown data from multiple depths at each location was designed. The data-worth analysis revealed that the design honours sensitivities of the predictions of interest to parameters. The designed network provided relatively high data-worth for minimizing uncertainty in the drawdown prediction at locations of interest.

1. Introduction

Onshore gas industry dominated by coal seam gas has been growing in Australia over the last decade. Unlike conventional gas, production of coal seam gas often involves extraction of large volumes of water from the coal beds. Where freshwater aquifers are connected to the coal beds, it is important to assess and monitor whether depressurization of coal seams cause pressure drawdown in the connected freshwater sources over the life and beyond the gas development project. Monitoring of groundwater using dedicated observation networks is important to monitoring the status quo of the resource subjected to stressors like mining and gas, irrigation, and or climate change.

Prediction and monitoring of groundwater impacts caused by the large scale onshore gas resource development activities is challenged by the fact that such activities focus on coal seams or reservoirs in deeper parts of sedimentary basins where conventional groundwater monitoring datasets are often sparse. Consequently, there are large-scale uncertainties about the hydrogeological characteristics of the aquifers and intervening aquitards involved. In data-poor areas, an iterative procedure of assessment of groundwater impacts from resource development should be adopted to inform adaptive management. The essential steps in such a management framework comprises three steps. These are: (a) predictive analysis of impacts using models informed by existing datasets; (b) identifying the types and locations for new monitoring data that can inform the predictions of interest; and (c) refining the predictions and management actions, including make-good arrangements where needed. Such an iterative approach can help to ensure that groundwater impacts are quantified in advance, whilst acknowledging the uncertainties. The validity of predictions and the extent of any potential impacts can be verified by collecting compliance monitoring data and managed adaptively with an appropriate level of certainty. This can also help regulators to ensure that resource development projects achieve minimal impact to the environment and water resources. In this study, we developed and applied a workflow of predictive modelling and monitoring network design that focus on identifying the data gaps and optimizing monitoring strategies. These two foci can progressively help reduce uncertainty and improve confidence in predicted groundwater impacts from anthropogenic activities like mining and gas development.

Numerous research studies in the past have employed optimization methods in conjunction with numerical simulation models to inform optimal monitoring strategies [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]. Many of these studies have been developed in the context of water quality monitoring and/or detection of contaminant migration from point sources [17,18,19,20,21,22,23,24,25]. The objectives of optimization considered in these studies include minimization of variance of the statistical estimator, minimization of uncertainty over the predicted spatial domain and strength of contaminant plumes, minimizing the cost of monitoring and maximization of the worth of data collected while satisfying budgetary constraints.

Such monitoring network optimization studies have previously not been developed in the context of monitoring groundwater drawdown impacts from CSG development. Previous practical approaches to design groundwater monitoring networks in regions of CSG development [26] have focused on the footprint of the planned CSG development, with the objective of achieving specified monitoring bore density within long-term affected areas and providing background monitoring and information required for supporting regional flow models. Sreekanth et al. [27] developed a theoretical basis for designing monitoring networks in the case of contaminant transport monitoring, using spatial basis functions and stochastic optimization. The proof-of-concept of the method was tested for designing a monitoring network for monitoring contamination risks in proposed well field for managed aquifer recharge. However, the use of simulation and optimization for monitoring network design when impacts can potentially arise over large spatial extents horizontally and vertically in the case of deep sedimentary basins is challenged by several factors. One of the biggest challenges of using predicted variables (groundwater heads or concentrations) for monitoring network design is that the veracity of the design is dependent on the model predictions.

Linear uncertainty analysis provides a powerful tool to identify the influence of model parameters on predictions and uncertainties [28,29]. Data-worth identification using linear uncertainty analysis have been applied in the recent past in groundwater modeling studies [30,31]. Middlemis and Peeters [32] discussed the value of relative composite sensitivity analysis as used in the above studies for informing the relative value of measurements in informing predictions. Recently, Sundell et al. [33] presented a value of information analysis (VOIA) to determine the need for additional information when assessing the effect of alternative designs in hydrogeological systems. Our study and the proposed workflow, in addition to optimization of the monitoring design, enables us to explore the predictive dependencies and veracity of the design independent of the predicted variables by employing these methods. The workflow employs linear analysis techniques that uses only the model sensitivities in computing the data-worth instead of absolute values of predictions at the designed monitoring locations.

Another major challenge is the scalability of conventional simulation-optimization methods to regional scale modelling problems. Historically, optimization of monitoring networks has largely focussed on informing monitoring decisions in groundwater remediation contexts. In the context of groundwater drawdown induced by onshore gas extraction, groundwater models need to cover spatial domains extending hundreds of kilometres horizontally and many hydrostratigraphic layers (aquifers and aquitards) vertically. Informing the design of groundwater monitoring networks using groundwater models in such contexts requires numerous simulations from regional scale models. Often, these models take several hours to complete a single simulation run [26,34].

Another major challenge is the difficulty in accounting for predictive uncertainty of these large-scale models in monitoring network design. A significant source of uncertainty stems from a paucity of data relating to hydraulic characteristics of deep aquifers, aquitards and coal seams. This is especially true in the green field sites where minimal onshore gas development has occurred and where the only available information about hydraulic characteristics is from the coal exploration bores. Observations of water levels from deep bores and continuous record of groundwater pressures are almost always absent in such regions, resulting in substantial uncertainty in the model predictions of pressures and pressure changes over large spatial scales.

Often, highly parameterized models that have hundreds or thousands of parameters are required to adequately capture the uncertainty in spatially varying hydraulic characteristics to appropriately estimate the model prediction uncertainty. In such cases, deterministic simulations using groundwater models cannot provide reliable information for major investment decisions like drilling of monitoring bores. Probabilistic approaches that help to estimate the predictive uncertainty and their spatial characteristics are better suited to identify the type of data and the locations of its measurement that will enable to reduce the uncertainty in the predictions of interest.

The main aim of this study is to address the practical challenges of design and verification of groundwater monitoring network using an integrated methodology and apply it to a complex groundwater flow model for groundwater monitoring network design and explore prediction-independent verification of the design. The remaining sections of this paper are organized in the following order. Section 2 describes the development of the proposed methodology and workflow. Section 3 describes various steps in the application of the methodology for the selected case study. Section 4 describes and discusses the results, and Section 5 presents the conclusions.

2. Materials and Methods

The workflow presented in this study for optimal design of groundwater monitoring network has three essential steps: (1) Probabilistic predictive analysis of the variables of interest (e.g., drawdown) and their spatial variability using a numerical groundwater model; (2) monitoring network design using the predictive analysis in conjunction with spatial basis function and global optimization algorithm; and (3) independent verification of the relative data-worth of observations collected at the designed optimal locations using linear analysis. The workflow is diagrammatically represented in Figure 1.

2.1. Probabilistic Predictive Analysis

A numerical groundwater model developed for simulating the variables of predictive interest forms the basis of this task. The model should ideally be developed in such a way that uncertainties in the input parameters are variables can be propagated through the model to quantify the uncertainties in the predicted variables. Calibration of the model results in identification of parameter sets that honour the observation data sets by satisfying a threshold objective function. Predictive analysis is subsequently undertaken by running Monte Carlo simulations of the model scenario for all combinations of parameters that satisfy the threshold objective function. These predictive simulations then form the essential data set for doing the optimal monitoring network design in step 2. The details of the model development, calibration and Monte Carlo simulations used in our case study is described in detail in Section 3.

2.2. Optimal Monitoring Network Design

We extended and applied the optimization method proposed in Sreekanth et al. [27]. Sreekanth et al. [27] tested the applicability of spatial basis functions in reproducing predicted patterns in the transport of contaminant plumes in an aquifer and used that information in designing a monitoring network to minimize the error in predicting the contaminant plume over the spatial domain using measurements at limited number of monitoring points in the aquifer. In this component of our present study, we applied the same method and tested the scalability and adaptability of the method to a larger dimensional problem in terms of the number of variables and patterns for which spatial basis functions were to be developed. In this study we used spatial basis functions to reproduce D-max and T-max patterns in five different model layers corresponding to different aquifers and aquitard formations. In practical terms, this provides the additional advantage of optimizing the placement of multi-level piezometers to monitoring propagation of depressurization towards risk receptors that are separated from the coal seam gas target formations by means of hydraulically distinct formations. Such sentinel monitoring networks are important to provide early warnings and identify mitigation measures before important risk receptors are impacted.

Spatial Basis Functions and Monitoring Network Optimization

The methodology for development of spatial basis functions in this study used orthogonal basis function expansion. As noted in [35], there are several choices for spatial basis functions, including empirical orthogonal functions (EOFs), Fourier basis functions, splines, and wavelets, among others. EOFs have been used in several successful model simulations in many disciplines (see [35]). Following the successful application in groundwater modelling by Sreekanth et al. [27] we use EOFs to project the groundwater simulations of predicted groundwater head drawdown in multiple model layers into a lower dimensional space. Critically, these EOFs account for the spatial dependence as well as the interlayer and model output dependencies. The theoretical basis of the development of spatial basis function using this approach is described in detail in Sreekanth et al. [34].

In this study, we used the differential evolution optimization algorithm [36] in conjunction with the spatial basis functions to optimize monitoring network. In this case, implementation of the DE algorithm in the reduced-rank space minimizes the uncertainty of the predicted drawdown pattern for the project region, regardless of what the true values of D-max and T-max across space are.

We made use of n = 500 runs of the probabilistic groundwater model to determine the optimal locations of monitoring wells. Initially, we undertook an investigation into whether the number of model runs undertaken, sufficiently covered the range of outcomes possible from the model simulation. This was undertaken by cross-validation where we sampled without replacement from our 500 simulations to create smaller sets of simulations. Increasingly larger in-sample sets of simulations were used to explore how well the constructed basis functions could predict the remaining out-of-sample model runs.

The probabilistic simulations of the impacts of CSG-induced drawdown provided the best available information regarding the range of potential impacts. Each model prediction from this set is an equally probable sample from the distribution of potential impact corresponding to a plausible realization of aquifer connectivity in the sedimentary basin. The objective of the monitoring network design was to find a set of m well locations, denoted $\mathit{w}={\left({\mathit{w}}_{\mathbf{1}},{\mathit{w}}_{\mathbf{2}},\dots ,{\mathit{w}}_{\mathit{m}}\right)}^{\mathit{T}}$ amongst all possible well locations given by

$$W\equiv \left\{w_i\in {\left\{1,2,\cdots ,p\right\}}^m:w_i\ne w_j \forall i,j\in \left\{1,\cdots ,m\right\}\right\}$$

(1)

Given m model runs each consisting of s spatial fields each containing p pixels, our aim is to solve the optimization problem [27]

$${\mathit{w}}^{*}\equiv {\mathrm{argmin}}_{w\in W}f\left(\mathit{w}|\mathit{X}\right)$$

(2)

The objective function $\mathit{f}(\mathit{w}|\mathit{X})$ is a function of the well locations and its value is conditional on the matrix, X, containing a representative sample of m possible model realizations. Herein, we define the objective function as

$$f\left(\mathit{w}|\mathit{X}\right)={\displaystyle \sum _{j=1}^m}{\displaystyle \sum _{k=1}^{sp}}{\left|{({\tilde{\mathit{x}}}_j)}_k-{\left({\mathit{x}}_j\right)}_k \right|}^2$$

(3)

where ${({\tilde{\mathit{x}}}_j)}_k$ denotes the kth element of the predicted vector where our observations are taken from model run j and ${({\tilde{\mathit{x}}}_j)}_k$ denotes the kth element of the vector from model run j.

Differential evolution [36] was used to solve the global optimization problem. The stochastic optimization procedure implemented using the DEoptim package [37] for the R statistical programming language [38] was used for this purpose.

This method was employed to identify locations for 10 monitoring bores within areas of predicted drawdown impact. The choice of the maximum number of bores was arbitrary in this case reflecting practical limitation that the total investment is often decided a priory based on available budget. The 10 locations were identified for multi-level piezometers that would collect groundwater head data from 5 different depths corresponding to 2 coal formations, 2 inter-burden formations, and the potentially impacted aquifer.

2.3. Data-Worth Analysis

Monitoring network design described in the previous section used the non-linear predictive analysis-based Monte-Carlo simulations of CSG-induced groundwater drawdown changes in multiple model layers. The data-worth analysis described in the current section employs the linear uncertainty analysis technique based on first order second moment-based methods [39] to investigate the utility and worth of the data collected at the potential monitoring bores. While the monitoring network design was based on trends in absolute values of predicted drawdown, the data-worth analysis is based on sensitivities of the prediction and future observation (monitored at the designed locations) to model parameters. Thus, data-worth analysis does not require the actual values of parameters or observations. Instead, only the sensitivities of predictions and observations with respect to parameters are used. This implies that the relative worth of future observations in informing the predictions of interest can also be calculated using these approaches. In this study, we calculated the relative worth of 100 drawdown data points measured at the 10 monitoring bore locations in informing maximum drawdown predictions.

The theoretical bases of such linear analysis approaches have been discussed in [30,31,40,41]. We applied a data-worth analysis method that is based on Schur’s complement. Schur’s complement for linear uncertainty analysis can be viewed as a form of Bayes equation under the assumption of a linear model and a multivariate Gaussian distribution to describe the distribution of parameters, forecasts, and observation noise [31,41,42,43]. The Schur’s complement approach, implemented in the software utility pyEMU [41], uses the following equation to estimate the posterior parameter covariance matrix ${\overline{\mathbf{\Sigma }}}_{\mathit{\theta }}$ as

$${\overline{\Sigma }}_{\theta }={\Sigma }_{\theta }- {\Sigma }_{\theta }{\mathit{J}}^T{\left(\mathit{J}{\Sigma }_{\theta }{\mathit{J}}^T+{\Sigma }_{ϵ}\right)}^{-1}\mathit{J}{\Sigma }_{\theta }$$

(4)

where ${\Sigma }_{\mathit{\theta }}$ is the prior parameter covariance matrix, ${\Sigma }_{\mathit{ϵ}}$ is the epistemic observation noise covariance matrix, and J is the Jacobian matrix of partial first derivatives of observations with respect to parameters. This equation highlights the behaviour of the inversion process [41].

The first term (i.e., ${\Sigma }_{\mathit{\theta }}$) represents the parameter uncertainty prior to inversion, and the second term $ie {\Sigma }_{\mathit{\theta }}{\mathit{J}}^{\mathit{T}}{\left(\mathit{J}{\Sigma }_{\mathit{\theta }}{\mathit{J}}^{\mathit{T}}+{\Sigma }_{\mathit{ϵ}}\right)}^{-\mathbf{1}}\mathit{J}{\Sigma }_{\mathit{\theta }})$) encapsulates the inversion process, through the Jacobian matrix and both parameter and observation covariance, as mapping of information from observations to parameters. The estimation of prior uncertainty for the prediction is identical across these two approaches. We used this formulation to evaluate prediction uncertainty in CSG-induced drawdown impacts at chosen risk receptors in this study. For this purpose, the ${\Sigma }_{\mathit{\theta }}$ matrix was obtained from the observed parameter covariance structure underpinned by the variogram. Similarly, the ${\Sigma }_{\mathit{ϵ}}$ matrix representing the observation error variance was obtained from the limited calibration analysis.

Prior and posterior uncertainty estimates for forecasts, defined as ${\mathit{\sigma }}^2{}_{\mathit{s}}$ and ${\overline{\mathit{\sigma }}}_{\mathit{s}}^2$ respectively, can be easily calculated by projecting the requisite parameter covariance matrix to the forecast output space using a forecast sensitivity vector

$${\sigma }^2{}_s={\mathit{y}}^T{\Sigma }_{\theta }\mathit{y}$$

(5)

and,

$${\overline{\sigma }}^2{}_s={\mathit{y}}^T{\overline{\Sigma }}_{\theta }\mathit{y}$$

(6)

where y is the vector of prediction sensitivity to each parameter. The details of the implementation of this method in pyEMU is described in [41].

3. Case Study

The Gunnedah Basin (northern NSW, Australia) has significant CSG resources and the state Government is currently considering the proposal for gas development from coal seams in the Maules Creek and Hoskissons Coal formations. The Gunnedah Basin underlies parts of the geological Surat Basin and hydrogeological Great Artesian Basin in this region which is a major source of fresh groundwater. The aquifers of the Namoi alluvium that overlies the Surat Basin formations are also important sources of fresh groundwater that are extensively used for irrigation and other uses. CSG development from the Gunnedah Basin requires extraction of large volumes of water from the coal seams. One of the potential environmental impacts of CSG development is the propagation of depressurization into the aquifers that overlie the Gunnedah Basin formations resulting in groundwater drawdown in the aquifers, and potentially affecting water availability for irrigation and other beneficial groundwater users and groundwater dependent ecosystems [36]. The extent of groundwater drawdown that may result in the aquifers largely depends on the volume of water extracted for gas production, the hydraulic characteristics of aquitards and other inter-burden formations in between the coal seams and the aquifers and water balance component of the aquifers. The plan and cross section view of the study area is shown in Figure 2.

In deep sedimentary basins with no history of gas production, field data to constrain these characteristics of the groundwater system are largely uncertain and model-based predictive analysis of potential impacts is challenged by data-scarcity resulting in substantial uncertainties over large spatial domains. Ideally, predictive analysis of impacts, in such circumstances, should be undertaken in a probabilistic fashion and the knowledge derived from uncertainty analysis should inform strategies for data collection to minimize uncertainty in future prediction of impacts in an iterative manner. In this study, we make use of probabilistic prediction of impacts and predictive uncertainty analysis to identify groundwater bores and other risk receptors that are potentially impacted by gas development. Additionally, predictive uncertainty analysis is used to identify optimal locations of future monitoring bores that can provide groundwater head measurements that will help constrain prediction of the propagation of depressurization. Furthermore, the proposed workflow also permitted us to do an independent data-worth analysis to validate the usefulness of measurements at these locations to inform the predictions of our interest. The implementation of the method for this case study is described in the following.

3.1. Numerical Groundwater Model

A 15-layer numerical groundwater model that was recently developed for the study area [36,44,45] to assess groundwater impacts from coal seam gas development in the Gunnedah Basin. We adopted this model with a modified and improved parameterization scheme in this study. MODFLOW-USG [46] was adopted as the preferred modelling platform considering the complexity of the geology comprising two distinct sedimentary basins (Surat and Gunnedah). This was especially useful to realistically represent pinching out layers and allowed model layers to follow the geometry of the geological formations and boundaries.

The extraction of CSG from the Hoskissons Coal and Maules Creek Formation (refer Figure 2) targets in the Gunnedah Basin has been proposed. Most of the important assets including stock and domestic and irrigation bores and groundwater-dependent ecosystems in the region rely on water from important aquifers in the region including alluvial and other deeper formations, such as the Pilliga Sandston. These formations were represented by independent numerical layers in the model. The formations that are present between the alluvium and the Pilliga Sandstone were conceptualized as inter-burden layers with variable hydraulic characteristics. Similarly, the formations between the Pilliga Sandstone and the two coal formations were represented in the numerical groundwater model by means of inter-burden layers with distinct effective hydraulic characteristics. The basement rock under the Maules Creek Formation is represented by means of another layer in the numerical model layer. This resulted in a numerical model architecture with nine layers to represent the hydrostratigraphy. Each inter-burden layer in the hydrostratigraphy was further divided into three layers in the numerical model grid to characterize realistic vertical gradients within aquitards. This resulted in a numerical model grid with 15 layers corresponding to the hydrostratigraphic units shown in Figure 2.

An earlier version of the model used a simplified parameterization approach which used depth-dependent hydraulic characteristics in the model layers [44]. In this study, we updated the model with a highly parameterized approach for undertaking predictive analysis of CSG-induced drawdown changes. The highly parameterized approach enabled a better representation of the potential spatial heterogeneity in the hydraulic properties of the model layers and was underpinned by measured data available from core samples [47], downhole porosity logs and geostatistical models that conform to this dataset. This approach was also more appropriate for the methods used in this study for explicit analysis of the model sensitivity and data-worth of observations and parameters to many spatially explicit predictions.

A total of 1672 parameters were used for parameterizing the groundwater model. Most of these parameters were used for characterizing the heterogeneity and spatial variability in the hydraulic properties in multiple layers of the model. Pilot points [48] were used as the parameterization device for this purpose. The remaining parameters were used to represent flood and irrigation recharge and river conductance (all of which are spatially variable). The spatial variability adopted in this parameterization was based on the data analysis and upscaling of hydraulic properties undertaken by another study that looked at aquitard properties in this gas development region [47].

Turnadge et al. [47] collected core samples of key aquitard sequences within the Gunnedah and Surat basin formations in the region and correlated measured permeability with available geophysical borehole data. Porosity–permeability relationships were derived based on this and were applied to downhole porosity logs obtained from 97 exploration wells located across the Gunnedah Basin [47]. These were then used to estimate vertical distributions of permeability at various locations across the Gunnedah Basin. They applied upscaling approaches to use these core scale aquitard vertical hydraulic conductivity (Kv) values in a regional-scale numerical groundwater model [49]. Furthermore, they used this upscaled data to fit geostatistical variogram models to characterize the spatial co-variance of hydraulic properties in the aquitard sequences. Two parametric models were used to characterize the experimental variograms: a spherical and an exponential model. The models were fitted to experimental variograms using least squares curve fitting. We adopted a spherical variogram with a Sill of 0.764, Nugget of 0.327 and range of 129 km to define the co-variance structure of hydraulic properties horizontally in different model layers based on this previous study. The highly parameterized modelling approach used in the current study is underpinned by this spatial covariance structure for characterizing prior uncertainty in hydraulic properties.

We undertook probabilistic simulation analysis to predict the impacts of a proposed CSG development project in the Gunnedah Basin that considered drilling 850 CSG wells near Narrabri in NSW. The drawdown was predicted as the difference between the baseline and CSG development cases. This was accomplished by taking the difference between model predictions of two possible states of GAB groundwater resource—one corresponding to no depressurization and the other corresponding to depressurization based on a generic case of CSG asset development. This approach assumes that all other uses of groundwater remain unchanged over this period and quantifies only the changes in flux and water balance resulting from depressurization only.

3.2. Model Calibration and Predictive Analysis

Calibration of a groundwater model using field observations (e.g., groundwater level/pressure) is most valuable when the information content of the dataset is sufficient to constrain the model. While groundwater observation from many bores were available for the Namoi alluvial and water table aquifers, only a small proportion of the bores screen deeper formations in the Surat and Gunnedah basins. Intuitively, these observations were not very useful for constraining the model when the predictions of interest were the propagation of CSG-induced drawdown in GAB and other formations. This is quantitatively demonstrated in the data-worth analysis results presented in the following section. Nonetheless, the model was used together with the parameter estimation software suite (PEST) to enable linear uncertainty and data-worth analyses, based on first order second moment (FOSM) sensitivities. However, rather than constraining the predictions to a narrow range based on calibrated parameters that are less relevant for the predictions concerned, we chose to do the predictive analysis using a Monte Carlo simulation framework which sampled the model parameters from their full prior-distribution. The probabilistic modelling analysis was used to simulate 500 plausible CSG-induced groundwater drawdown patterns in multiple model layers. Specifically, the two variables of predictive interest in this study for monitoring network design are: (i) maximum CSG-induced drawdown (D-max); and (ii) the time elapsed before we observe the maximum drawdown (T-max) following the extraction of gas. These patterns were used in the monitoring network design as described in Section 2

4. Results and Discussion

The calibration error analysis undertaken using PEST provided useful insights concerning errors in the model structure and parameters. The model calibration scatter plot and distribution of calibration error is shown in Figure 3. The distribution of errors is also used to inform the observation error variance that is required for the linear uncertainty analysis. The error distribution showed that nearly 75% of the errors are in the range −10 m to 10 m. The errors in this range are more likely to be occurring because of uncertainties in the model parameters and boundary conditions. A small proportion of large errors indicate the possibility of some model structural errors and errors in aquifer assignments.

Five hundred model simulations resulted in 500 predictions of plausible CSG-induced groundwater drawdown in different model layers corresponding to plausible realization of uncertain model parameters. The 5th, 50th, and 95th percentile of predicted drawdown in the Pilliga Sandstone is shown in Figure 4. The 50th percentile drawdown prediction indicated that CSG depressurization would not result in a drawdown over 0.2 m in the GAB aquifer. The 95th percentile, indicated drawdown up to 1.3 m in the GAB aquifer. The drawdown was calculated as the difference between the baseline and CSG development cases. The predictive analysis resulted in simulation of CSG water extraction in the range 4.4 GL to 107.1 GL over the life of the gas project. The wide range reflects the large uncertainty in estimation of CSG water production. This range encompasses the base and high case water production of 35.5 GL and 87.1 GL respectively estimated by the proponent.

As described in Section 2, the variance in the prediction arises not only from the simulations of wide range of CSG water extraction, but also because of plausible wide range of values for the hydraulic properties of aquitards that separate the coal seams from the shallower aquifers.

Given the nature of water and gas extraction from deeper sedimentary basin it is intuitive that early detection of the propagation of depressurization would require sentinel monitoring bores, ideally multi-level piezometers placed in different formations at different depths from the surface. Such a monitoring network is designed in this study using the proposed method. Subsequently, the data-worth of monitoring in existing monitoring bores and the designed network is assessed independent of the predicted variables.

4.1. Optimal Monitoring Network

Simulated patterns of CSG-induced maximum drawdown (D-max) obtained from the Monte Carlo simulation were used to design the optimal monitoring network. The contour plots of D-max for model layers 8 and 12 are shown in Figure 5. Such drawdown patterns (D-max and T-max) simulated for six model layers (layers 6, 7, 8, 10, 12, and 14) were used in the monitoring network design.

As discussed in the methods section, we used EOFs to project these model-simulated groundwater drawdown patterns into a lower dimensional space. These EOFs account for the spatial dependence as well as the interlayer and model output dependencies. Ten EOFs were used to explain most of the variability in the groundwater model simulations and were calculated using singular value decomposition of the model output. Ten EOFs were used in this study because 10 monitoring wells are designed. This allows to estimate 10 coefficients for each of the EOFs when fitting the surface to the observations at the bores.

Importantly, the first EOF accounts for the most variability of all the EOFs, followed by the second one and so on. Figure 6 shows that beyond 200 model runs, the mean squared error in reconstructing the out-of-sample model runs decreased little and indicates that over 200 model runs has adequately explored most of the patterns of variability in the model output.

Examples of spatial basis functions generated for D-max and T-max in model layer 6 are shown respectively in Figure 7 and Figure 8. Such spatial basis functions developed for the D-max and T-max patterns in the six model layers were used in conjunction with the differential evolution algorithm to identify 10 optimal locations where multi-level piezometers could be installed within the CSG development region. The choice of the locations were restricted to areas where all six model layers of interest were present. The 10 optimal locations identified using this approach are shown in Figure 9.

4.2. Data-Worth Analysis

Before applying data-worth analysis to investigate the utility of the designed we used it to explore the predictive dependencies between parameters, observations and predictions. Data-worth analysis is undertaken for specific predictions of interest. In this study we evaluated the data-worth in predicting the CSG-induced drawdown at model locations where risk receptors were identified. A total of 57 model nodes were selected for data-worth analysis. These nodes corresponded to locations, within which, risk receptors were identified and fell within 95th percentile predicted drawdown in the near surface aquifer (from all model layers). These nodes are shown in Figure 10. The data-worth analysis provided a number different insights that relates the predictions of interest to the groundwater system characteristics and the data that can inform these characteristics. Specifically, the following four insights are reported in the following: (a) relative parameter group contributions to drawdown prediction uncertainty at the receptors; (b) relative parameter contribution of hydraulic properties of intervening layers for the predictive uncertainty of vertical propagation of depressurization; (c) parameter estimatability using existing and designed future monitoring network; and (d) data-worth of the designed optimal monitoring network in comparison to the existing observation network.

Relative contribution of five different parameter groups to drawdown prediction uncertainty at the 57 locations is shown in Figure 11.

Important information that is readily gleaned from these plots is that the prediction uncertainty of drawdown has distinct contributions from different parameter groups. The specific storage groups is the predominant parameter group for many receptors. Relatively higher contributions from the specific storage (ss), river (z), and horizontal hydraulic conductivity (kh) groups are observed compared to the vertical hydraulic conductivity (kv) group. Uncertainty contributions of the z group are enhanced when the receptor is closer to the river network.

For example, for the receptor nodes 67927 in outcrop areas of layer 6 and have a small 95th percentile predicted drawdown as shown in Figure 10. It may be observed from Figure 11 that the parameter contribution of the river (z) group is highest for these receptors. This implies that the drawdown incurred at these bores depends on how fast/slow the CSG-induced flux losses from this area is replenished by additional inflows from (or reduction in base flows to) the river.

This also implies that for CSG-induced drawdown at this scale (noting that maximum drawdown is close to zero for the 50th percentile), drawdown incurred is more sensitive to how fast/slow the small flux losses from the aquifer are compensated by release of water from storage, river or recharge, than to the propagation of declining pressures in deeper formations to these aquifers.

This is because the aquitard formations that are present between the coal formations and aquifers dampen the propagation of depressurization, and pressure changes are time lagged before they reach the GAB and water table aquifers. It is important to note that hydraulic characteristics of deeper formations are important parameters that govern the propagation of CSG-induced drawdown to the upper layers. This was quantitatively evaluated by analyzing the parameter contributions to prediction uncertainty of drawdown propagation through model layers between the water table aquifer and the coal bearing formations. The hydraulic property parameter groups’ contribution to drawdown prediction at 20 locations within model layers 12 is shown in Figure 12. It shows that the propagation of CSG-induced drawdown across the intervening layers depends significantly on the hydraulic characteristics of these and adjacent layers.

The third set of analyses investigated how currently existing and future dedicated monitoring networks can inform these hydraulic characteristics and parameters of aquifers and aquitards and other processes. The top 10 parameters for which prior to posterior uncertainty reduction can be achieved by the currently existing monitoring data set and future monitoring using the designed dedicated monitoring network are compared in Figure 13.

Subsequently, we investigated the data-worth of piezometers in multiple model layers obtained from the monitoring network design in comparison to the existing observation data

The top three parameters for which maximum uncertainty reduction can be achieved using the observations of groundwater head in the existing monitoring bores are parameters affecting modelled recharge and horizontal hydraulic conductivity (rfl, rdr, and kh6491) as shown in Figure 13a Other key parameters for which uncertainty reduction can be achieved by monitoring bores in the shallow aquifers included the horizontal hydraulic conductivity (kh) and specific storage (ss). These findings reiterate the earlier result that the groundwater drawdown experienced in the shallow water table aquifer will be largely dependent on how quickly the recharge and river processes, and release of water from the storage, compensates for the water lost due to CSG development. On the other hand, the parameters for which maximum uncertainty reduction is achieved by monitoring deeper zones, especially in the intervening formations, included hydraulic characteristics around the coal seams (pilot points ss318753, kv318753, and kh318753 correspond to the coal seam model layer 14) as indicated in Figure 13b. It is also noteworthy that majority of the top 10 parameters that achieve an uncertainty reduction, correspond to vertical hydraulic conductivity. All the above findings from the data worth analysis confirms that monitoring the vertical and horizontal propagation of groundwater drawdown across multiple layers influence the impact at the risk receptors and hence having multi-level piezometers within these drawdown extents is most important for monitoring.

Subsequently, we also used the data-worth analysis to identify the relative worth of groundwater head monitoring in multilevel piezometers at the designed locations. Relative data-worth of monitoring data from different depths in informing drawdown predictions was calculated as a relative increase in prediction uncertainty owing to the removal of datasets from the calibration data set. Relative data worth was calculated for five different observation groups (i.e., data from potential multi-level piezometers in model layers 7, 8, 12, and 14 names respectively as monit7, monit8, monit12, and monit14 and the data from existing monitoring bores named the ‘obsgroup’). Figure 14 illustrates the relative worth of monitoring observations in model layers 7, 8, 12, and 14 and existing observations from shallower wells (obsgroup) in informing drawdown prediction at 20 locations within model layer 6 corresponding to the Pilliga Sandstone. It may be observed that the greatest data-worth for this purpose arise from observations in model layer 7, that lies immediately below model layer 6. This implies that measuring groundwater head in the Purlawaugh formation immediately below the Pilliga Sandstone (and evaluating it in comparison with pre-CSG heads) would be of relatively highest worth in informing the propagation of drawdown into the Pilliga Sandstone formation. It is also noteworthy that the design of the multi-level piezometer network provides relatively high data-worth for minimizing uncertainty in the drawdown prediction at most of the 20 locations.

5. Conclusions

A workflow comprising probabilistic groundwater modelling, reduced-rank spatial prediction, and linear analysis was developed and applied for designing optimal monitoring network that provided maximum data-worth in informing prediction of propagation of coal seam gas-induced groundwater drawdown changes into a shallower sandstone aquifer. The workflow uses ensemble predictions of drawdown corresponding to wide range of model parameters to design the network. The data-worth analysis of observations obtained from the network is assessed independently using linear uncertainty analysis. The linear analysis uses information from the sensitivity of the observations to model parameters and is independent of the absolute value of the predicted variables. Thus, the linear analysis enables to explore the predictive dependencies to provide insights into the veracity of the designed network. The workflow is generic and could be applied for the design of purpose-built monitoring networks for monitoring groundwater impacts from anthropogenic and other stresses.

In this study, the probabilistic groundwater flow simulations using a numerical model considered wide ranges of variability in the CSG water production rates and hydraulic characteristics of several intervening model layers to account for substantial uncertainty in the groundwater system and gas development processes at large spatial scales. This allowed us to probabilistically characterize potential groundwater drawdown patterns in the region. Reduced rank spatial prediction was used to characterize dominant trends in these spatial patterns using empirical orthogonal functions. These EOFs were then used in conjunction with the differential evolution algorithm to identify optimal locations for multi-level piezometers: helping to minimize predictive uncertainty in groundwater drawdown. The value of the monitoring bores as a network for minimizing prediction uncertainty across the region is demonstrated by designing a 10-bore monitoring network that collects drawdown data from multiple depths at each location. Data-worth analysis validates the suitability of these optimal monitoring locations. The study demonstrates that both prediction variables (e.g., drawdown distribution) and independent sensitivity information (data-worth) can provide valuable information for designing monitoring networks.

In this study, we used the data-worth based on linear uncertainty analysis as a tool for independent verification of the designed monitoring network. It is possible to include the data-worth directly in the design to ensure that monitoring wells are placed at locations and times where data-worth is maximized. There is opportunity for future work in this domain and we are currently pursuing that.