2.6.2.7.3 Sensitivity analysis


In the design of experiment phase of the uncertainty analysis workflow (Figure 6, Section 2.6.2.1) a set of parameter combinations and corresponding model outputs is generated by evaluating the GW AEM and the alluvial MODFLOW models a large number of times, to train the emulators that replace the model chain in the uncertainty analysis. This training dataset is also used for the formal sensitivity analysis. Ten thousand different combinations are generated of the parameters of the GW AEM (listed in Table 5 of Section 2.6.2.6.1), the alluvial MODFLOW models (listed in Table 6 of Section 2.6.2.6.2) and the AWRA-L model (listed in Table 7 of Section 2.6.1.5.1 in companion product 2.6.1 for the Gloucester subregion (Zhang et al., 2018)).

The number of samples is considered sufficient if, emulators trained with that number of samples are able to accurately predict the model response for unsampled parameter combinations (Sacks et al., 1989). Due to the large number of predictions that are required from the model chain, it is not possible to guarantee that the sampling scheme is sufficient for each and every individual prediction. The accuracy of the emulators is, however, individually evaluated and emulators that do not satisfy the accuracy criteria are omitted from the uncertainty analysis (see Section 2.6.2.7.4).

A pragmatic choice is made to generate 10,000 parameter combinations of the entire model chain in the design of experiment. The size of the training set, 10,000 samples, is a trade-off between the sampling density of parameter space and computational resources available (both in run-time and storage space). The bioregional assessment (BA) modelling team used the Pearcey and Bragg clusters, which are part of the High Performance Computing facilities of CSIRO (CSIRO, 2016) to evaluate the model runs for the design of experiment and the subsequent training and running of the emulators for the uncertainty analysis.

The design of experiment set of parameter combinations is generated using a maximin Latin Hypercube design (see Santner et al., 2003, p. 138). Latin Hypercube sampling is a stratified random sampling of a multivariate parameter space designed to provide an even coverage of parameter space (Iman, 2008). The maximin Latin Hypercube design is generated like a standard Latin Hypercube design, one design point at a time, but with each new point selected to maximise the minimum Euclidean distance between design points in the parameter space. Points in the design span the full range of parameter values in each dimension of the parameter space, but also avoid redundancy among points by ensuring that no two samples are too close together (since nearby points are likely to have similar model output).

Figure 21 shows histograms of the 10,000 parameter combinations. The parameters are sampled uniformly on a log10 scale (K_IB_intercept, K_CS _intercept, S_IB_intercept, S_CS_intercept, KvKh, Kfh, Kfv, ne) or on a natural scale (K_IB_slope, K_CS_slope, S_IB_slope, S_CS_slope). In the design of experiment, no correlation between parameters is specified.

2.6.2.7.3.1 Analytic element model

Figure 21

Figure 21 Histograms of the Latin Hypercube sampling of the parameter space of the regional analytic element groundwater model for the sensitivity analysis for the Gloucester subregion

Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 1)

Figure 22 shows the actual hydraulic conductivities and storage assigned to the model based on this Latin Hypercube sampling for two nominal depths: 250 m and 1000 m. As mentioned in Section 2.6.2.6, the numerical model implements lower threshold values to ensure numerical stability. These are also incorporated in Figure 22.

Figure 22

Figure 22 Range of hydraulic conductivity (K) and storage (S) for interburden and coal seams at 250 m and 1000 m for the Gloucester subregion

1000 m corr incorporates the lower threshold of 1 x 10–7 m/day for hydraulic conductivity and 1 x 10–6 for storage

Data: Bioregional Assessment Programme (Dataset 1)

Ten thousand parameter combinations are evaluated, 7819 of these simulations produced dmax and tmax values; 2181 returned a ‘Not a Number’ value for at least one receptor. These failed runs represent parameter combinations for which the GW AEM could not arrive at a converging solution. Figure 23 shows box and whisker plots of the parameter ranges that gave rise to successful and failed model runs. In these plots, the box represents the interquartile range (25th to 75th percentile) and the whiskers extend from the 5th to the 95th percentile. The median is indicated by a horizontal line through the interquartile box. Outliers, points below the 5th percentile or above the 95th percentile, are depicted as black points.

Failed model runs appear to be limited for parameter combinations that lead to low vertical conductivity of the interburden (low K_IB_intercept, high K_IB_slope and low KvKh) and high interburden storage (high S_IB_intercept) in combination with elevated hydraulic conductivity/low storage of the coal seams (high K_CS_intercept, low K_CS_slope, low S_CS_intercept).

These parameter combinations result in a fast and large drawdown in the coal seam, which cannot be replenished quickly enough through the interburden because of its low Kv. The pumping rate required to maintain the prescribed head in the CSG wells becomes increasingly small, which gives rise to numerical precision errors that in turn result in ‘Not a Number’ values for the CSG pumping rate. CSG pumping is implemented as a two-stage process (see Section 2.6.2.5). Pumping rate is computed in the first stage from a prescribed head boundary condition that will be implemented in a second stage as a prescribed flux.

In the failed runs, the prescribed flux is not available and hence the CRDP run fails. A case can be made to assign a zero-prescribed CSG pumping rate to the failed runs and run these again to at least simulate the cumulative impact of coal mining in the weathered zone for these parameter combinations. Another option is to limit the range of the parameters in the sensitivity analysis run, to ensure failed runs do not occur. As both options represent a manipulation of the model runs, these options are not implemented and in the further analysis, the failed model runs are omitted. The information from the failed runs will, however, be incorporated in specifying the prior distributions of the parameters and their covariance for the uncertainty analysis (Section 2.6.2.8).

Figure 23

Figure 23 Boxplots of the parameter combinations of the regional analytic element groundwater model that lead to successful and failed model runs for the Gloucester subregion

A failed run is a run which returned a ‘Not a Number’ value for dmax and tmax for at least one receptor. Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 1)

Figure 24

Figure 24 Scatterplots of parameter values versus maximum coal seam gas water production rate (Q_csg) for the Gloucester subregion

Each dot represents a single simulation. The red lines show the median value of dmax over the parameter range spanned by the line. The green line is the maximum water production rate threshold of 1.1 ML/day based on AGL Energy Limited (2015) predicted water extraction rates. Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 1)

The values of dmax and tmax are discussed in depth in Section 2.6.2.8, when the final posterior distributions of the dmax and tmax for each receptor are available. To gain a better understanding of model behaviour, however, the effect of the parameters on the maximum CSG production rate (Figure 24) and on the dmax values for receptor GLO_156, 2 km north from the Stratford Mining Complex and 2 km south of the Rocky Hill Coal Project (Figure 25), are examined. Rather than plotting the dmax values directly, they are transformed by taking the cubed root. This transformation increases the resolution for both positive and negative values.

From Figure 24 it is apparent that the maximum CSG production rate is controlled by the hydraulic conductivity of the coal seams (K_CS_intercept and K_CS_slope). The log of the maximum water production rate increases linearly with increasing hydraulic conductivity. An increase in the intercept of the exponential decrease function leads to higher conductivity values at depth. A decrease in the slope of equation 6 in Section 2.6.2.6 will have the same effect, higher conductivity at depth.

Almost half of the parameter combinations of the design of experiment result in maximum CSG water production rates in excess of 1.1 ML/day, the threshold value above which maximum water production rates are considered unrealistically high (see Section 2.6.2.7.1). This indicates that a sizeable proportion of the wide parameter range sampled in the design of experiment (Figure 25 and Figure 22) gives rise to unrealistically high CSG water production rates. Section 2.6.2.8 describes how the parameter range is constrained in the uncertainty analysis to avoid these unrealistic parameter combinations in the posterior parameter ensemble that is used to generate the posterior prediction ensemble.

Figure 25 shows the variation of the maximum difference in drawdown, dmax, in function of the parameter values of the design of experiment. One of the most striking features in Figure 25 is the presence of negative dmax values. These values appear to be limited to parameter combinations with high K_CS_Intercept and low S_CS_intercept values. These extreme and unlikely parameter combinations cause numerical instabilities which result in failed runs or negative drawdown due to additional coal resource development. The numerical instability occurs because the large fluxes generated by high hydraulic conductivity result in extremely high hydraulic gradients due to the low storage.

The red lines are the median of dmax in the parameter range spanned by the line. For this particular receptor, GLO_156, it is apparent that the two most influential parameters are S_CS_intercept and K_CS_intercept.

The intercept of the depth–S relationship for the coal seam, S_CS_intercept, shows a marked response in which drawdown increases with decreasing storage. This is likely an effect limited to the weathered zone, where, per definition of storage, for a volume of water extracted from an aquifer, small storage values result in large drawdowns and large storage values result in small drawdowns.

Decreasing the intercept of the depth–K relationship for the coal seam, K_CS_intercept, decreases both the K of the weathered zone and the K of the coal seams. A decrease of hydraulic conductivity in the weathered zone has a result that the cone of depression becomes deeper, as the resistance is higher for water to flow laterally to replenish the water extracted by the mines.

A detailed discussion for each receptor of these relationships is beyond the scope of the report. Figure 26 does show the sensitivity index of the cube root of dmax for each parameter–receptor combination, calculated with the density-based sensitivity index introduced by Plischke et al. (2013). Larger values of the sensitivity index indicate higher sensitivity of the prediction to the parameter. The cube root transform is chosen to compensate for the skewed distribution of most of the dmax values. It is apparent that the K and S intercept of the coal seams (and thus the weathered zone) are the dominant factors for most receptors.

Figure 25

Figure 25 Scatterplots of parameter values versus cubed root of maximum difference in drawdown (dmax), due to the additional coal resource development (ACRD) for receptor GLO_156 for the Gloucester subregion

Each dot represents a single simulation. The red lines show the median value of dmax over the parameter range spanned by the line. Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 1)

Figure 26

Figure 26 Sensitivity indices of the cube root of maximum difference in drawdown (dmax) for all parameter-receptor combinations, ordered from north (top) to south (bottom)

Larger values indicate higher sensitivity of a receptor to a parameter. Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 1)

2.6.2.7.3.2 Alluvial MODFLOW models

In the design of experiment of the sensitivity analysis of the Gloucester model chain, 10,000 parameter combinations of the MODFLOW Avon and Karuah models are generated and evaluated for both models. Figure 27 shows histograms of the parameter values in the design of experiment. The hydraulic conductivity of the weathered zone (Khw) and the alluvium (Kha) and the drainbed conductance are uniformly sampled in log10 space while the other parameters, the specific yield (Sy), the recharge multiplier (Rmult) and the constant head offset (dh) are uniformly sampled in native space.

Figure 27

Figure 27 Histograms of the Latin Hypercube sampling of the parameter space of the alluvial Avon and Karuah models for the sensitivity analysis for the Gloucester subregion

(a) Saturated hydraulic conductivity of top alluvial layer (Kha (m/day)); (b) Saturated hydraulic conductivity of lower weathered layer (Khw (m/day)); (c) Specific yield of the top alluvial layer (Sy (-)); (d) Hydraulic conductance of lower boundary of drain bed (Dc (m2/day); (e) Multiplier for monthly recharge (Rmult (-)); and (f) Depth to water in the lower weathered layer (dh (m)). Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 5)

From the 10,000 runs, 723 failed to converge. Figure 28 shows that failed model runs are associated with high hydraulic conductivity values for the alluvium (Kha) and a low specific yield (Sy). This physically unrealistic combination of parameters allows for the alluvial layer to drain too quickly and become dry. This causes numerical instability preventing the model to converge.

Figure 28

Figure 28 Boxplots of the parameter combinations of the Avon and Karuah models that led to successful and failed model runs for the Gloucester subregion

(a) Saturated hydraulic conductivity of top alluvial layer (Kha (m/day)); (b) Saturated hydraulic conductivity of lower weathered layer (Khw (m/day)); (c) Specific yield of the top alluvial layer (Sy (-)); (d) Hydraulic conductance of lower boundary of drain bed (Dc (m2/day); (e) Multiplier for monthly recharge (Rmult (-)); and (f) Depth to water in the lower weathered layer (dh (m)). Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 5)

Figure 29

Figure 29 Histograms of drainage and upward fluxes simulated in the design of experiment runs for the Avon and Karuah models, compared to the target distribution, for the Gloucester subregion

Qdr = drainage flux; Qup = alluvium – weathered zone exchange flux

Data: Bioregional Assessment Programme (Dataset 5)

Figure 29 shows histograms of the simulated drainage fluxes and upward fluxes into the alluvium for the Avon and Karuah models, compared to their target range from Table 7. The simulated drainage fluxes are close to the target distributions, while the current range of parameters evaluated in the design of experiment appears to underestimate the upward flux, especially for the Karuah model. During the uncertainty analysis, a Markov chain Monte Carlo simulation is carried out to select parameter combinations that are in agreement with the target distributions.

Figure 30 shows the response of the drainage and upward fluxes to variations in parameter values for the Avon and Karuah models.

The drainage flux in both Avon and Karuah models is controlled by the hydraulic conductivity of the alluvium and the recharge multiplier, where increases in these parameters lead to increases in drainage flux. The other parameters do affect the drainage flux, but to a much lesser extent.

The magnitude of the upward flux is completely dominated by the hydraulic conductivity of the weathered zone. The direction of flow appears to be controlled primarily by the offset of the specified head boundary condition. Secondary controlling factors are the hydraulic conductivity of the alluvium and the recharge multiplier. The upward flux is a function of the hydraulic gradient between the alluvium and the weathered zone. Large dh values result in fluxes towards the weathered zone. Elevated groundwater levels in the alluvium, through low Kha or high Rmult, will also result in fluxes towards the weathered zone.

Figure 31 illustrates the response of drawdown due to additional coal resource development at the receptor locations to variations in parameter values for two representative receptors: GLO_147 in the Avon model and GLO_248 in the Karuah model. These plots include the two parameters that control drawdown in the weathered zone from the regional model, namely K_CS_intercept and S_CS_intercept.

For both receptors it is clear that the most important parameter is the hydraulic conductivity in the weathered zone (Khw), followed by the weathered zone parameters of the regional model. The hydraulic conductivity and storage in the alluvium are of lesser importance, as is the recharge multiplier and drainbed conductance and constant head offset. This indicates that the drawdown due to additional coal resource development is dominated by the magnitude of the change in flux at the bottom of the alluvium.

Figure 30

Figure 30 Scatterplots of the simulated drainage (1st and 3rd column of figures) and upward fluxes (2nd and 4th column of figures) versus parameter values for the Avon (1st and 2nd column of figures) and Karuah (3rd and 4th column of figures) MODFLOW models for the Gloucester subregion

The red lines indicate the median y-axis value over the range spanned by the width of the line. The orange band indicates the target range of the y-axis. Qdr = drainage flux; Qup = alluvium – weathered zone exchange flux. Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 5)

Figure 31

Figure 31 Scatterplots of maximum difference in drawdown (dmax) versus parameter values for receptors GLO_147 ((a) through (h)) and GLO_248 ((a) through (h)) for the Gloucester subregion

Refer to Table 5 in Section 2.6.2.6.1 for definitions of terms.

Data: Bioregional Assessment Programme (Dataset 5)

In Figure 32 this sensitivity is quantified for all dmax parameter–receptor combinations with the density-based sensitivity index introduced by Plischke et al. (2013). Rows in the figure without colour coding indicate receptors with an impact too small to analyse. The sensitivity indices confirm that, in general terms, S_CS_intercept, K_CS_intercept and Khw dominate the magnitude of dmax for most receptors. Note that the water balance targets are not able to greatly constrain these parameters.

Last updated:
5 November 2018
Thumbnail of the Gloucester subregion

Product Finalisation date

2018
PRODUCT CONTENTS