Application of the Iterative Ensemble Smoother Method and Cloud Computing: A Groundwater Modeling Case Study

Abstract

Numerical groundwater modelling to support mining decisions is often challenging and time consuming. Simulation of open pit mining for model calibration or prediction requires models that include unsaturated flow, large magnitude hydraulic gradients and often require transient simulations with time varying material properties and boundary conditions. This combination of factors typically results in models with long simulation times and/or some level of numerical instability. In modelling practice, long run times and instability can result in reduced effort for predictive uncertainty analysis, and ultimately decrease the value of the decision-support modelling. This study presents an early application of the Iterative Ensemble Smoother (IES) method of calibration-constrained uncertainty analysis to a mining groundwater flow model. The challenges of mining models and uncertainty quantification were addressed using the IES method and facilitated by highly parallelized cloud computing. The project was an open pit mine in South Australia that required predictions of pit water levels and inflow rates to guide the design of a proposed pumped hydro energy storage system. The IES calibration successfully produced 150 model parameter realizations that acceptably reproduced groundwater observations. The flexibility of the IES method allowed for the inclusion of 1493 adjustable parameters and geostatistical realizations of hydraulic conductivity fields to be included in the analysis. Through the geostatistical realizations and IES analysis, alternative conceptual models of fractured rock aquifer orientation and connections could be conditioned to observation data and used for predictive uncertainty analysis. Importantly, the IES method out-performed finite difference methods when model simulations contained small magnitude numerical instabilities.

1. Introduction

Model predictive uncertainty analysis is a fundamental component of informative modeling for decision-making [1]. Although this is recognized by the groundwater modelling community [2], in practice, uncertainty analysis is often lacking due to practical limitations of model complexity, client expectations and consulting time/budget constraints [3,4]. Ideally an uncertainty analysis approach follows a Bayesian method of considering the prior uncertainty in model parameters as statistical distributions that are refined through the collection of observation data [5,6]. However, the application of rigorous Monte Carlo sampling methods [5,6,7] to do this analysis with computationally expensive numerical simulations of groundwater flow is not currently feasible. Commonly applied methods for Monte Carlo analysis with efficient workarounds based on subspace techniques [8] require a computational effort that scales with the dimensionality of the parameters, limiting the feasible parameterization of computationally expensive models, and potentially inducing predictive error through calibration with a simplified parameterization scheme [9,10]. The use of paired simple and complex models has been advanced to address some of these limitations [10,11]. However, recent advances in predictive uncertainty techniques in the field of petroleum reservoir modeling [12,13] provide Iterative Ensemble Smother (IES) methods for rigorous Monte Carlo sampling-based uncertainty analysis. IES methods are computationally efficient and not limited by the dimensionality of the parameterization, so a paired simple complex approach is not needed. Recently, the IES method has been incorporated into a freely available model independent software package [14] that has been used in this study.

Even with the efficiencies of the IES method, the Monte Carlo analysis of a computationally expensive numerical model is demanding and required a high degree of parallelization. High throughput computing employing cloud computing resources was used for this study and has been noted by previous studies to provide practitioners with access to high performance computing capability on an as needed basis [15,16,17] and make computationally demanding methods feasible.

This paper presents a case study of an early application of the IES method for calibration-constrained uncertainty analysis with a computationally expensive real-world model built to support mining decisions. Calibration of an ensemble of model parameter sets, representing alternative conceptual models of fractured rock aquifer orientations, was conducted using a simulation of open pit mine development conditioned on a relatively extensive calibration dataset.

The experience of this project demonstrates that the IES methods enabled with cloud computing resources can enable rigorous model uncertainty analysis despite complex, challenging models and finite time constraints. It was found in this study that the IES method was also more amenable for conducting uncertainty analysis with models containing numerical oscillations that hinder the application of methods-based finite difference model sensitivities.

2. Materials and Methods

The Kanmantoo Copper Mine in South Australia, located approximately 50 km south east of the city of Adelaide (Figure 1), operated from 1970 to 1976. In 2011, mining was recommenced and conducted from 2011 to 2019; an open pit, 600 m wide and 400 m deep, was excavated. The purpose of the groundwater modeling project described in this paper was to aid the design of a proposed pumped hydro energy storage (PHES) scheme, for the decommissioned open pit at the Kanmantoo Copper Mine. The PHES project is due to start in 2021. Under the PHES scheme, a pond in the base of the open pit, and an upper storage pond, will be used to store energy to support development of renewable energy projects in the region. The water level recovery in the pit at the end of mining and rate of groundwater inflow into the open pit during the operation of the PHES project, have large implications for the project design. So, the main predictions of interest are recovery pit water levels two years after mine dewatering has ended, and water levels in the open pit and surrounding aquifers during and after the proposed fifty-year PHES system operation.

For the purposes of conceptual model development and calibration, groundwater levels recorded in the five years prior to the recent mining activity (pre-2011) were assumed to be representative of steady state.

The following subsections summarize the development of the predictive groundwater model, model calibration and the application of uncertainty analysis. Observation data available at the site included pre-2011 static water-level measurements and timeseries of water-level fluctuations during the recent mining activity (2011–2019). A numerical groundwater flow model was constructed based on a conceptual model of the local hydrogeology and hydrology. The model was developed to simulate predictions of interest; pre 2011 static water levels; and water-level changes during the recent mining activity. The latter two simulations were used for model calibration purposes. Initial model calibration attempts highlighted the challenges posed by small-scale numerical instabilities in the model and was additional motivation for application of the IES method. Fine scale spatially variable parameterization for IES analysis allowed for heterogeneity and geostatistical simulations of alternative starting models, which enabled better representation of the fractured rock aquifer conceptual model.

2.1. Groundwater Conceptual Model

The geology at the Kanmantoo Mine is highly complex with extensive fracturing and metamorphoses [18]. Fresh bedrock is overlain by highly weathered bedrock. Groundwater mainly occurs in discrete fracture zones within the fresh bedrock, or in highly heterogeneous, strip aquifers. The key aspects and interpretations of the hydro-stratigraphy that were considered during development of the numerical predictive model were:

• Groundwater flow is controlled by discrete fracture zones within relatively low hydraulic conductivity geologic units.
• Shallow weathered rock has reduced hydraulic conductivity due to weathering filling fractures, or greater hydraulic conductivity due to lower loading and larger fracture apertures [19].
• Fracture aperture and resulting hydraulic conductivity is likely to be inversely proportional to depth. Deeper rocks carry greater loads that can compress fractures [19].
• The orientation of large-scale faulting in the area is approximately north–south [18]. It was interpreted that the local-scale fractures and subsequent strip aquifers, also trend in a north–south direction.
• An old tailings storage facility used during historical mining, is located north of the open pit. It is interpreted that this region has higher hydraulic conductivity than the surrounding weathered rock. This feature may have an impact on some of the upper storage pond designs in the predictive simulations.

Watershed boundaries and potential surface water locations were identified based on flow accumulation modeling, and regional groundwater flow patterns were characterized by contouring publicly available groundwater-level data [20] (Figure 1). Regionally, groundwater levels predominantly follow the topography and trend towards the southeast. The watershed and sub watershed boundaries combined with regional water-level data and were used to select an approximately 10 km by 10 km model domain and set model boundary conditions.

2.2. Observation Data

Pre-2011 groundwater-level observations were available from 32 industrial monitoring bores near the mine site and 71 locations across the model domain from public databases. Transient groundwater levels measured during the recent mining activity (2011–2018) were collated from 16 locations near the mine site.

Transient observations of groundwater changes due to mining advancement, provides information about the local groundwater flow system around the open pit, and are likely to reduce uncertainty in predictions of mine inflow. The observation data was used to create both steady state and transient targets for model calibration. Pre-2011 static water levels were simulated as steady state for calibration. Progressive open pit mining was simulated in a transient model for calibration to observed water-level changes during recent mining activity (2011–2018).

2.3. Numerical Model

A three-dimensional numerical model of groundwater flow was constructed using the FEFLOW [21] finite element software. The regional study area (RSA) comprising the entire model domain and was selected so that where possible, boundaries aligned with surface water divides for the larger watershed or subwatersheds and no flow boundary conditions could be applied. Cauchy type boundary conditions were applied to boundaries that did not align with surface water divides. Reference water levels for the Cauchy boundary conditions were based on the contoured regional water-level observations. A local study area (LSA) with dimension of 3.5 km by 3 km was defined around the open pit and proposed upper storage ponds. Model mesh refinement and fine scale model parameterization occurred in the LSA.

2.3.1. Mesh and Layers

The FEFLOW three-dimensional (3D) model construction process occurred in two stages. Initially, a plan view, two-dimensional (2D) triangular mesh was constructed using the triangle [22] mesh generator and, was then extended to three-dimensions through layer definition. The initial 2D mesh was based on linear and point features that included: observed surface water; outlines of three proposed upper storage pond designs; 30 m contours of five different stages of mine pit development; and two water source wells. The 2D model was extended to a thirteen-layer, 3D model with a combination of flat and undulating layers mirroring topography to maintain a 10m minimum layer thickness.

During model development, the mesh was iteratively refined and altered to achieve simulation convergence. Nodes in the final 2D mesh were spaced approximately 400 m apart in the RSA and refined to approximately 15 m spacing near the open pit. The mesh had 11,000 nodes and was a combination of triangles and tetrahedra (Figure 2). Model layer thickness varied from 10 m to 30 m over the depth of the open pit.

Surface elevation in the RSA ranged from 500 m above the Australian Height Datum (AHD) in the north-west, to 75 m AHD in the south and, was approximately 200 m AHD near the open pit. The base of pit elevation varied during the simulated period from a pre-2011 level of 100 m, to 30 m at the end of the transient calibration period (2018), to the final pit design depth of −160 m expected to be reached later in 2019. This was a change in depth of 260 m over eight years.

2.3.2. Boundary Conditions

Boundary conditions are shown on Figure 2 and included natural no-flow conditions at hydrologic divides, and Cauchy boundaries elsewhere to capture the regional inflow/outflow. The model base consisted of a no-flow boundary to replicate impermeable bedrock. Recharge was applied at the model surface, and discharge into creeks and the pit lake was represented by seepage face boundary conditions. Groundwater flow was conceptualized as recharging at surface and, discharging into the mine pit, creeks, rivers, and the Cauchy boundary on the south side of the model domain.

2.3.3. Representation of Mining

The excavation of 260 m from the open pit from 2011 to 2019 was represented in the model in eight stages. Mine development and dewatering were represented using time-varying material properties and boundary conditions.

Dewatering of the mine was simulated using time-varying, constrained, specified head boundary conditions that were set at nodes on the base of the pit for each stage. At the start of each stage, the nodes along the base were activated and assigned with conditions that linearly changed from the elevation of the previous stage to the elevation of the current stage over the elapsed time span. Boundary conditions were constrained to only release water from the system and inflows were precluded. In this way, the eight mining stages were simulated as specified head boundaries lowering the water table to the mining stage surface level.

To simulate excavation of sediment from the model domain as pit depths increased, material property values of the excavated sediment were changed from background values to values of high hydraulic conductivity, low specific storage, and low unsaturated porosity to represent open pit conditions.

The property values were linearly transitioned from initial values to values representing open pit conditions over the span of each stage.

The transient simulation required two hours of computation time. No model convergence issued were identified in initial testing.

2.4. Initial PEST Calibration

An initial model calibration using PEST software [23] was conducted with a small number of zone-based parameters to identify reasonable, average parameter values, to use as starting points for further calibration and uncertainty analysis. The PEST calibration was conducted in two parts: steady state and transient.

2.4.1. Calibration Data Processing

Changes in measured values over space and time can yield more information to calibrations than the use of absolute measurement values alone [1]. Groundwater observation data was processed into calibration targets with this concept in mind. The pre-2011 groundwater-level measurements were processed into spatial groundwater-level gradient calibration targets in addition to the absolute water levels. The difference between every unique combination of groundwater-level measurement was used. The transient observation data was processed into timeseries of changes from the first absolute measurement (drawdown) and, timeseries of changes in groundwater level between successive measurements. The total set of calibration targets of each type is listed in

Table 1. Processed calibration targets for the groundwater model.

Observation Group	Number of Observations
Public Data Groundwater Levels	70
Public Data Depth to Water	70
Public Data Level Gradients	2211
Monitoring Data Groundwater Levels	32
Monitoring Data Level Gradients	351
Monitoring Data Drawdown	482
Monitoring Data Temporal Changes	466
Total Calibration Targets	3682

.

2.4.2. Initial PEST Calibration Parameterization

Four unique zones were defined in the model for hydraulic parameter assignment:

• Weathered rock from surface to 30 m depth.
• Fresh rock below 30 m depth.
• Compressed rock below 30 m elevation (~170 m depth near the open pit).
• Disturbed rock at the old tailing’s storage facility.

Adjustable model parameters for the initial PEST calibration included values of horizontal and vertical hydraulic conductivity, specific storage and, specific yield for each of the four zones. Uniform recharge and conductance values at four Cauchy boundaries were also included for a total of 21 parameters.

2.4.3. Initial PEST Calibration Results

The PEST steady state calibration to the pre-2011 static head observations and gradients, reduced observed-to-simulated misfit by altering the recharge parameter. However, other parameters remained largely unchanged. The transient calibration, including simulation of mining progression, resulted in no reduction in observed-to-simulated misfit. A PEST utility JACTEST [24] was used to examine the finite difference derivatives. The horizontal hydraulic conductivity parameter in the fresh rock was selected to test derivatives because the outputs of both steady state and transient simulations were sensitive to hydraulic conductivity. The derivative test was achieved by running the model seven times with small parameter perturbations. If the relationship between a model parameter and the simulated observations is well represented by the finite difference approach, then running the model several times with small perturbations will produce simulation results that have a linear trend that is reasonably well approximated by the slope derived from two or three of the points.

Figure 3 shows the simulation results for a transient drawdown observation and a static hydraulic head measurement. Seven model runs with one percent perturbations to the fresh rock horizontal hydraulic conductivity parameter were used for this test. The consistent downward trend in simulated head values (Figure 3b) indicates that a finite difference approximation to the model captured local model behavior for the steady state simulation. The absence of a trend in the drawdown simulations (Figure 3a) suggests that a finite difference approximation of model behavior using only two or three points may be a poor approximation of local model behavior. The poor finite difference approximations explain the inability of the initial PEST-based transient calibration to improve observed-to-simulated misfit. The lack of linear trend in the simulated drawdowns is likely due to numerical instabilities in the model. This drawdown prediction and associated observation have absolute values exceeding 100 m and are in relatively close proximity to the active open pit mine. Consequently, in the numerical simulation this prediction is based on a point in space in time that has been influenced by the transition of many model elements in at least five layers from saturated to unsaturated conditions and time varying material property and boundary conditions representing the simulation of mine development. For these reasons, the relatively small-scale oscillations in predicted drawdown value between simulations with slightly different hydraulic conductivity values (Figure 2a) are not unexpected with this nonlinear model.

To proceed with model calibration and uncertainty analysis, three options were considered:

(1)

Iterative refinement of model layering and mesh until a more stable simulation of mine development could be achieved.

(2)

Discarding the transient data from the model calibration and uncertainty analysis.

(3)

Proceeding with an IES model calibration and uncertainty analysis technique, not based on the finite difference approximation, with the existing model, and determining whether the observed numerical oscillations would hinder the IES process through trial.

Given that the run time of the transient model was two hours, a model with sufficient layers to provide the vertical refinement necessary to produce a stable solution for finite difference-based methods, was anticipated to have a prohibitively long runtime. The transient observation data was interpreted to contain important information for constraining the predictive uncertainty, so discarding it was unappealing.

It was decided to proceed with the existing numerical model and use the IES method for calibration and parameter uncertainty analysis described in the following sections. Use of the IES method with the current model was deemed appropriate because the model simulations with individual parameter sets produced reasonable results and convergence, and the magnitude of the numerical oscillations observed in the finite difference testing was likely to be much smaller than the combined effects of measurement and model parameter uncertainty.

The only aspect of the initial PEST calibration that was carried forward into the IES analysis was a revised initial value of mean recharge. The initial calibration was done as a check on the model simulation rather than any necessary part of the following IES analysis.

2.5. IES Methodology

The IES method of generating alternative, calibration-constrained, parameter realizations begins with an initial ensemble of alternative model parameter sets that are considered reasonable based on background system knowledge, or samples from the prior in Bayesian terminology. The calibration simulation is run for each of the model parameter sets in the ensemble. Based on the parameter and simulated observation variability about mean values for each parameter set, an approximation to a Jacobian matrix of parameter sensitivities can be calculated; see Equation 4 in [12]. The approximated Jacobian matrix is used in a regularized Levenberg Marquart algorithm to calculate update vectors for the entire ensemble of parameter sets that minimize the sum square error between observations and model simulations similar to traditional calibration with PEST [1]. Regularization is applied to penalize parameter values departing from the initial ensemble values unless necessary to reduce measurement misfit. This process is iteratively repeated until a sufficient match to the observation data has been achieved with the ensemble of parameter sets. The result of the process of the IES analysis is a collection of alternative parameter sets that are consistent with background system knowledge and observations, or samples of the posterior distribution in Bayesian terminology. The approximation of model parameter sensitivities based on a small number of simulations in the ensemble has been observed to cause a collapse of the ensemble into unrealistically similar sets [12], however the application of regularization is designed to prevent this. Furthermore, localization methods [13] are intended to deal with spurious correlations that are introduced through estimating sensitivities with a small ensemble. The IES method implemented in PEST++ IES [14] takes advantage of several numerical simplifications described in [12] and follows Algorithm 1 given in the appendix of [12].

2.6. Parameterization for IES

The IES method is not limited by number of parameters. Therefore, a geostatistical approach was adopted, where alternate hydraulic conductivity fields in the LSA were simulated. This approach captured the heterogeneous, fractured rock, strip aquifers in the conceptual model. Pilot point parameterization on a regular 100 m grid was used to define the hydraulic conductivity in the fresh rock over the LSA and enabled geostatistical simulations (Figure 4B). The fresh rock hydraulic conductivity was used for the primary geostatistical simulations because it appears at the same depth at which mining occurred during the transient calibration period.

In the RSA, fresh rock horizontal hydraulic conductivity values were defined at pilot point intervals of approximately 1000 m. An additional set of pilot points was used to define the other parameters and had a spacing of 500 m in the LSA and 1000 m in the RSA. (Figure 4A). The second pilot point set defined:

(1)

Multiplier fields to decrease the hydraulic conductivity field in the weathered rock zone.

(2)

Multiplier fields to decrease the hydraulic conductivity field in the deep rock zone.

(3)

Specific storage parameter field.

(4)

Drainable porosity parameter field.

(5)

Recharge parameter field.

Vertical to horizontal hydraulic conductivity anisotropy was defined as a uniform zone-based parameter for the fresh, weathered and deep rock zones and, four independent conductance parameters were defined for the Cauchy boundary. The IES analysis had a total of 1493 adjustable parameters (

Table 2. Iterative Ensemble Smoother (IES) model parameters.

Parameter Group	Number of Parameters
Horizontal K Pilot points	811
Weathered K Multiplier Pilot points	135
Deep K Multiplier Pilot points	135
Specific Storage Pilot points	135
Drainable Porosity Pilot points	135
Recharge Pilot points	135
Boundary Condition Conductance’s	4
Vertical K anisotropy multiplier	3
Total Calibration Parameters	1493

).

2.7. Initial Parameter Generation and IES Runs

The sequential Gaussian simulation method [25], as implemented in the PEST utility FIELDGEN [26], was used to generate horizontal hydraulic conductivity parameter values for the dense grid of pilot points in the LSA (Figure 4B). The geostatistical structure was defined by two variograms to produce realizations of hydraulic conductivity fields consistent with the conceptual model of strip aquifers with approximately north–south orientation (

Table 3. Geostatistical structure used for hydraulic conductivity field generation.

Variogram	Type	Range(m)	Anisotropy	Bearing	Weight
1	Exponential	1000	10	−17° to 17°	0.4
2	Spherical	2500	5	−17° to 17°	0.2

).

All other parameters were randomly drawn from normal or uniform distributions (

Table 4. Initial parameter values and bounds.

Parameter	Unit	Initial Value	Upper Bound	Lower Bound
Fresh Rock Kh	m/s	1.2 × 10−8	1 × 10−4	1 × 10−9
Weathered Rock Kh multiplier	None	8.3	15	1
Deep Rock Kh multiplier	None	0.5	1	1 × 10−3
Recharge	Mm/year	12	35	1
Specific storage	1/m	2 × 10−7	2 × 10−6	2 × 10−8
Drainable Porosity	None	0.3	0.35	0.1
Weathered Rock Vertical Anisotropy	None	1	1	0.01
Fresh Rock Vertical Anisotropy	None	1	1	0.01
Deep Rock Vertical Anisotropy	None	1	1	0.01
North Boundary Conductance	1/s	1 × 10−5	1 × 10−3	1 × 10−9
South Boundary Conductance	1/s	1 × 10−5	1 × 10−3	1 × 10−9
West Boundary Conductance	1/s	1 × 10−5	1 × 10−3	1 × 10−9

).

Three hundred initial parameter realizations were produced by using the PEST utility RANDPAR to generate the RSA hydraulic conductivity, specific storage, drainable porosity, recharge, vertical anisotropy, multiplier field and Cauchy boundary condition conductance parameters. These parameter realizations were then merged with the LSA fresh rock hydraulic conductivity pilot point values produced by the geostatistical field generation. The first IES run iterated once with the full three hundred realizations and the best one hundred and fifty realizations with the lowest measurement misfit were selected for two additional iterations.

2.8. Cloud Computing Implementation

The Amazon EC2 cloud computing service was used to provide computational resources for the IES analysis. Distribution of model resource files and initiation of PEST++ slave processes was automated using PPAPI software [27]. Based on previous experience [16], fewer high CPU instances were selected for the analysis, and the computation was distributed over a network of up to nine cloud computing instances with thirty-six CPUs each for a total maximum of 324 CPUs. This model setup allowed up to three hundred model parameter sets to be run in parallel.

3. Results

The IES method, as implemented in PEST++ IES, includes a base parameter realization, where parameter values are set at initial specified values. The purpose of the base parameters is to take the place of the “calibrated model” in the analysis by starting at the user defined, most likely parameter values, rather than at random realization values. The base parameters for each parameter and zone were set as homogenous values based on the results of the steady state PEST calibration. Measurement objective function is defined as the weighted sum squared difference between simulated and observed data. Data weighting was done so that each of the measurement groups listed in Table 1 contributed equally to the measurement objective function when the initial base parameter set was used. The first iteration of IES reduced the mean measurement objective function value for all realizations by 78%. The second iteration reduced the mean measurement objective function by 91% relative to the starting value and, the third iteration caused a negligible objective function reduction.

3.1. Parameter Values

Figure 5 shows the parameter field progression during each iteration, from initial geostatistical fields, to calibrated models with similar features. The base and three hydraulic conductivity parameter realizations are shown (Figure 5).

All the parameter fields, including the horizontal hydraulic conductivity shown in Figure 5, became more similar with increasing calibration iterations. This similarity shows the large degree of heterogeneity in the parameter fields necessary to reproduce the calibration dataset.

3.2. Model Predictions

The water level in the pit lake at the proposed start of the PHES scheme, two years after mining cessation, was a prediction of interest.

For this predictive simulation, the simulation of mine development was run from 2011 to the 2019 final pit design. The material properties in the zones of the excavated pit had transitioned to a drainable porosity of 1.0 to simulate the open hole storage. Pit level recovery was simulated by removing seepage face boundary conditions that simulated open pit dewatering and, through applying boundary conditions to mimic evaporation from the surface of the rising pit lake. The simulation was run inside a python script. Additionally, at each timestep, the hydraulic head in the open pit was compared to mine pit geometry to estimate pit lake and saturated pit wall surface area. The surface area and an evaporation flux parameter were used to calculate a total evaporative loss rate from the system which was then applied using a well boundary condition at the base of the lake. Fifty parameter realizations with the lowest measurement misfit from the third IES iteration were used in the predictive modeling for the pit-level simulation. The use of fifty alternative parameter realizations for predictive modeling allowed for the predictive uncertainty to be expressed statistically. Figure 6 shows the predicted depth of the mine pit lake.

The mean predicted recovery depth in the pit from fifty alternative parameter realizations was 30.5 m with a standard deviation of 14.8 m. Predicted water depth varied from 1.6 m to 68.6 m. Therefore, it can be inferred that a groundwater-fed pit lake is likely to form in the decommissioned open pit within two years after mining ends. However, scenarios with low pit lake levels such as the minimum predicted depth of 1.6 m, indicate low groundwater flow to the pit, and cannot be excluded based on the current observation data.

4. Discussion

This study presents an application of the Iterative Ensemble Smother (IES) method of calibration-constrained, parameter uncertainty analysis, to a complex, groundwater model with numerical instabilities and a long run time. The results demonstrate that the IES method can be successfully applied to models with instabilities that hinder finite difference-based calibration and uncertainty analysis. High numbers of model parameters can be used with IES methods to allow for geostatistical simulations of parameter fields, and the number of model runs required per iteration is much lower than the number of parameters. The parameter realizations produced in this study are visually similar (Figure 5, bottom row) and intuitively may not be representing a full range of calibration-constrained variability. Regularization and localization for the IES method is an area of active research [13,14], and further developments may provide guidance on ensuring that maximum variability between parameter sets is maintained while the parameters are altered to reduce observation misfit. In this study, regularization was used to penalize the departure of parameter values from the initial generated value unless necessary to fit data. Localization methods to remove spurious correlations were not used. The necessity of subjective localization when using the IES methods decreases with an increasing number of model simulations, or ensemble size, per iteration [13]. The minimum number of recommended ensembles is recommended to be equivalent to the number of uniquely identifiable pieces of information in the calibration dataset [12]. In this study, a larger ensemble, of two to three times the estimated recommended minimum, was used to try to manage any spurious correlations [13]. With a model requiring more than two hours to execute, a highly parallelized approach to the IES analysis was necessary and, as this was an industrial consulting project without access to government or industry computing clusters, cloud computing was essential to the project. The initial effort in setting up the cloud computing was far outweighed by the benefit of parallelization for the IES analysis and predictive runs. Extensions to the implementation of the IES method to take advantage of dynamic cloud computing environments with low costs [17] is a potential direction for future development.

Based on the resulting similarity in the parameter fields, localization methods may have produced a better result. However, the variability in the predictions was sufficient to identify most likely outcomes for post mining pit water-level and low pit water-level outcomes that cannot be discounted on the basis of observation data and need to be considered as plausible. Cloud computing resources enabled using a larger number of realizations for the IES method and predictions in a reasonable timeframe, which can improve IES results [14] and allow a greater assessment of predictive uncertainty. In this case, a thorough predictive uncertainty analysis could be conducted because the IES method provided a one-step calibration and parameter uncertainty analysis and, the small-scale model instabilities did not prevent the IES method from improving the fit to observation data. IES analysis has been facilitated with model-independent software [14], so the methods described in this study are easily applicable to a wide range of models with variable spatial and temporal scales. The ensemble size necessary for appropriate approximation of sensitivities increases with the amount of information in the calibration dataset, so analysis on rich datasets may require large ensembles, though the computational burden can be managed using cloud computing, as was done in this study. The findings of this study do indicate that parameter localization [13] needs to be considered and applied to preserve variability even with larger than minimally necessary ensembles. This project was conducted in a consulting environment with a finite budget and timeframe. If more time was required to iteratively refine the calibration model so that finite difference methods could be applied, model run times would likely increase and prevent much effort being directed towards predictive uncertainty analysis. This would mean that decision makers would not be equipped with an understanding of the range in plausible values for important predictions for risk-based decisions.

5. Conclusions

This project has demonstrated the successful application of the IES method for assessing alternative conceptual models of connected fractured rock strip aquifers in a real-world groundwater modeling study to support mining decisions. The IES method is not limited by the number of parameters and allowed for geostatistical simulations of alternative parameter sets on a fine grid of pilot points. The outcome of this project, particularly given the numerically difficult model, supports further use of the IES method combined with geostatistical simulations of alternative conceptual models wherever possible. Where alternative conceptual models cannot be represented with grid- or cell-based geostatistical simulations, such as more zone-based geological models, then the application of the IES method could be applied in parallel to several alternative geological conceptual models to assess the impact of both the location of geological zone boundaries and the heterogeneous distribution of model properties within each geological zone on predictions of interest.