Next Article in Journal
Application of NSGA-II and Improved Risk Decision Method for Integrated Water Resources Management of Malian River Basin
Next Article in Special Issue
Bayesian Model Weighting: The Many Faces of Model Averaging
Previous Article in Journal
Evaluation of Temperature Profiling and Seepage Meter Methods for Quantifying Submarine Groundwater Discharge to Coastal Lagoons: Impacts of Saltwater Intrusion and the Associated Thermal Regime
Previous Article in Special Issue
Making Steppingstones out of Stumbling Blocks: A Bayesian Model Evidence Estimator with Application to Groundwater Transport Model Selection
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Groundwater Solutions Pty. Ltd., Melbourne 3031, Australia
2
Mining One Pty. Ltd., Melbourne 3000, Australia
3
Hillgrove Resources Pty. Ltd., Kanmantoo 5252, Australia
*
Author to whom correspondence should be addressed.
Water 2019, 11(8), 1649; https://doi.org/10.3390/w11081649
Submission received: 27 June 2019 / Revised: 6 August 2019 / Accepted: 7 August 2019 / Published: 9 August 2019

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.

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).

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).
All other parameters were randomly drawn from normal or uniform distributions (Table 4).
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.

Author Contributions

Conceptualization, A.V. and K.H.; methodology, K.H.; software, J.S. and K.H.; validation, K.H. and A.V.; formal analysis, K.H.; investigation, B.H.; resources, K.H.; data curation, K.H., A.V. and B.H.; writing—Original draft preparation, K.H. and E.W.; writing—Review and editing, E.W.; visualization, K.H.; supervision, K.H.; project administration, K.H.; funding acquisition, B.H.

Funding

This research received no external funding.

Acknowledgments

We would like to thank Hillgrove Resources, AGL Energy, and Mining One Consulting for the opportunity to contribute to this project for energy sustainability, to pursue emerging methods for quantitative groundwater modeling, and to share such methods in this forum.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Doherty, J. Calibration and Uncertainty Analysis for Complex Environmental Models; Watermark Numerical Computing: Brisbane, Australia, 2015; ISBN 978-0-9943786-0-6. [Google Scholar]
  2. Hunt, R.J.; Zheng, C. The Current State of Modeling. Ground Water 2012, 50, 330–333. [Google Scholar] [CrossRef] [PubMed]
  3. Delottier, H.; Pryet, A.; Dupuy, A. Why Should Practitioners be Concerned about Predictive Uncertainty of Groundwater Management Models? Water Resour. Manag. 2017, 31, 61–73. [Google Scholar] [CrossRef]
  4. Sanchez-Vila, X.; Fernàndez-Garcia, D. Debates—Stochastic subsurface hydrology from theory to practice: Why stochastic modeling has not yet permeated into practitioners? Water Resour. Res. 2016, 52, 9246–9258. [Google Scholar] [CrossRef]
  5. Tarantola, A. Inverse Problem Theory and Methods for Model Parameter Estimation; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2005. [Google Scholar]
  6. Tarantola, A. Popper, Bayes and the inverse problem. Nat. Phys. 2006, 2, 492–494. [Google Scholar] [CrossRef]
  7. Mosegaard, K.; Tarantola, A. Monte Carlo sampling of solutions to inverse problems. J. Geophys. Res. Solid Earth 1995, 100, 12431–12447. [Google Scholar] [CrossRef]
  8. Tonkin, M.; Doherty, J. Calibration-constrained Monte Carlo analysis of highly parameterized models using subspace techniques. Water Resour. Res. 2009, 45. [Google Scholar] [CrossRef] [Green Version]
  9. Knowling, M.J.; White, J.T.; Moore, C.R. Role of model parameterization in risk-based decision support: An empirical exploration. Adv. Water Resour. 2019, 128, 59–73. [Google Scholar] [CrossRef]
  10. Doherty, J.; Christensen, S. Use of paired simple and complex models to reduce predictive bias and quantify uncertainty. Water Resour. Res. 2011, 47. [Google Scholar] [CrossRef] [Green Version]
  11. Razavi, S.; Tolson, B.A.; Burn, D.H. Review of surrogate modeling in water resources. Water Resour. Res. 2012, 48. [Google Scholar] [CrossRef]
  12. Chen, Y.; Oliver, D.S. Levenberg–Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Comput. Geosci. 2013, 17, 689–703. [Google Scholar] [CrossRef]
  13. Chen, Y.; Oliver, D.S. Localization and regularization for iterative ensemble smoothers. Comput. Geosci. 2017, 21, 13–30. [Google Scholar] [CrossRef]
  14. White, J.T. A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensions. Environ. Model. Softw. 2018, 109, 191–201. [Google Scholar] [CrossRef]
  15. Hunt, R.J.; Luchette, J.; Schreuder, W.A.; Rumbaugh, J.O.; Doherty, J.; Tonkin, M.J.; Rumbaugh, D.B. Using a Cloud to Replenish Parched Groundwater Modeling Efforts. Ground Water 2010, 48, 360–365. [Google Scholar] [CrossRef] [PubMed]
  16. Hayley, K.; Schumacher, J.; MacMillan, G.J.; Boutin, L.C. Highly parameterized model calibration with cloud computing: An example of regional flow model calibration in northeast Alberta, Canada. Hydrogeol. J. 2014, 22, 729–737. [Google Scholar] [CrossRef]
  17. Hayley, K. The Present State and Future Application of Cloud Computing for Numerical Groundwater Modeling. Groundwater 2017, 55, 678–682. [Google Scholar] [CrossRef] [PubMed]
  18. Belperio, A.P.; Preiss, W.V.; Fairclough, M.C.; Gatehouse, C.G.; Gum, J.; Hough, J.; Burtt, A. Tectonic and metallogenic framework of the Cambrian Stansbury Basin—Kanmantoo Trough, South Australia. AGSO J. Aust. Geol. Geophys. 1998, 17, 183–200. [Google Scholar]
  19. Lavrov, A. Fracture permeability under normal stress: A fully computational approach. J. Pet. Explor. Prod. Technol. 2017, 7, 181–194. [Google Scholar] [CrossRef]
  20. Australian Groundwater Explorer: Water Information: Bureau of Meteorology. Available online: http://www.bom.gov.au/water/groundwater/explorer/map.shtml (accessed on 24 July 2019).
  21. Diersch, H.-J.G. FEFLOW Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer: Berlin, Germany, 2014; ISBN 978-3-642-38738-8. [Google Scholar]
  22. Shewchuk, J.R. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator. In Applied Computational Geometry Towards Geometric Engineering; Lin, M.C., Manocha, D., Eds.; Springer: Berlin, Germany, 1996; pp. 203–222. [Google Scholar]
  23. Doherty, J. PEST Model-Independent Parameter Estimation User Manual Part I: PEST, SENSAN and Global Optimisers 2016. Pesthomepage. Available online: http://www.pesthomepage.org/Downloads.php (accessed on 8 August 2019).
  24. Doherty, J. PEST Model-Independent Parameter Estimation User Manual Part II: PEST Utility Support Software 2016. Pesthomepage. Available online: http://www.pesthomepage.org/Downloads.php (accessed on 8 August 2019).
  25. Deutsch, C.; Journel, A. GSLIB: Geostatistical Software Library and User’s Guide, 2nd ed.; Oxford University Press: England, UK, 1997; ISBN 978-0-19-510015-0. [Google Scholar]
  26. Doherty, J. Groundwater Data Utilities Part B: Program Descriptions 2015. Pesthomepage. Available online: http://www.pesthomepage.org/Downloads.php (accessed on 8 August 2019).
  27. Schumacher, J.; Hayley, K.; Boutin, L.-C.; White, E. PPAPI: A Program for Groundwater Modeling Tasks in Distributed Parallel Computing Environments. Groundwater 2018, 56, 248–250. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Location of the Kanmantoo mine site and data used to develop a conceptual model of groundwater flow.
Figure 1. Location of the Kanmantoo mine site and data used to develop a conceptual model of groundwater flow.
Water 11 01649 g001
Figure 2. Numerical Model Mesh, Boundary Conditions, and Observation Data Locations. Model domain corresponds to Pink shaded area shown in Figure 1.
Figure 2. Numerical Model Mesh, Boundary Conditions, and Observation Data Locations. Model domain corresponds to Pink shaded area shown in Figure 1.
Water 11 01649 g002
Figure 3. Results of small perturbations to the bulk fresh rock horizontal hydraulic conductivity parameter for simulated (a) Drawdown and (b) Static head observation. The decreasing trend observed in static head simulations (b) indicates that the finite difference approximation captures local model behavior for the steady state simulation. The lack of a linear trend in the drawdown simulation (a) indicates that numerical oscillations likely prevent any finite difference approximation based on two or three points from approximating the local model behavior.
Figure 3. Results of small perturbations to the bulk fresh rock horizontal hydraulic conductivity parameter for simulated (a) Drawdown and (b) Static head observation. The decreasing trend observed in static head simulations (b) indicates that the finite difference approximation captures local model behavior for the steady state simulation. The lack of a linear trend in the drawdown simulation (a) indicates that numerical oscillations likely prevent any finite difference approximation based on two or three points from approximating the local model behavior.
Water 11 01649 g003
Figure 4. Pilot point locations for the IES calibration. (A) Regional study area (RSA) and local study area (LSA) pilot point locations for recharge, storage and hydraulic conductivity multiplier field parameters; (B) LSA pilot point locations for geostatistical simulations of fresh rock horizontal hydraulic conductivity fields.
Figure 4. Pilot point locations for the IES calibration. (A) Regional study area (RSA) and local study area (LSA) pilot point locations for recharge, storage and hydraulic conductivity multiplier field parameters; (B) LSA pilot point locations for geostatistical simulations of fresh rock horizontal hydraulic conductivity fields.
Water 11 01649 g004
Figure 5. Horizontal hydraulic conductivity fields on a log scale over the local study area (LSA) shown in Figure 4 for the base and three selected realizations through the three calibration iterations. This shows the progression of parameter fields from the starting model and initial geostatistical fields to calibrated models displaying similar features through iterations one and three.
Figure 5. Horizontal hydraulic conductivity fields on a log scale over the local study area (LSA) shown in Figure 4 for the base and three selected realizations through the three calibration iterations. This shows the progression of parameter fields from the starting model and initial geostatistical fields to calibrated models displaying similar features through iterations one and three.
Water 11 01649 g005
Figure 6. Predicted mine pit water depth two years after the cessation of mining using fifty parameter realizations from the IES analysis.
Figure 6. Predicted mine pit water depth two years after the cessation of mining using fifty parameter realizations from the IES analysis.
Water 11 01649 g006
Table 1. Processed calibration targets for the groundwater model.
Table 1. Processed calibration targets for the groundwater model.
Observation GroupNumber of Observations
Public Data Groundwater Levels70
Public Data Depth to Water70
Public Data Level Gradients2211
Monitoring Data Groundwater Levels32
Monitoring Data Level Gradients351
Monitoring Data Drawdown482
Monitoring Data Temporal Changes466
Total Calibration Targets3682
Table 2. Iterative Ensemble Smoother (IES) model parameters.
Table 2. Iterative Ensemble Smoother (IES) model parameters.
Parameter GroupNumber of Parameters
Horizontal K Pilot points811
Weathered K Multiplier Pilot points135
Deep K Multiplier Pilot points135
Specific Storage Pilot points135
Drainable Porosity Pilot points135
Recharge Pilot points135
Boundary Condition Conductance’s4
Vertical K anisotropy multiplier3
Total Calibration Parameters1493
Table 3. Geostatistical structure used for hydraulic conductivity field generation.
Table 3. Geostatistical structure used for hydraulic conductivity field generation.
VariogramTypeRange(m)AnisotropyBearingWeight
1Exponential100010−17° to 17°0.4
2Spherical25005−17° to 17°0.2
Table 4. Initial parameter values and bounds.
Table 4. Initial parameter values and bounds.
ParameterUnitInitial ValueUpper BoundLower Bound
Fresh Rock Khm/s1.2 × 10−81 × 10−41 × 10−9
Weathered Rock Kh multiplierNone8.3151
Deep Rock Kh multiplierNone0.511 × 10−3
RechargeMm/year12351
Specific storage1/m2 × 10−72 × 10−62 × 10−8
Drainable PorosityNone0.30.350.1
Weathered Rock Vertical AnisotropyNone110.01
Fresh Rock Vertical AnisotropyNone110.01
Deep Rock Vertical AnisotropyNone110.01
North Boundary Conductance1/s1 × 10−51 × 10−31 × 10−9
South Boundary Conductance1/s1 × 10−51 × 10−31 × 10−9
West Boundary Conductance1/s1 × 10−51 × 10−31 × 10−9

Share and Cite

MDPI and ACS Style

Hayley, K.; Valenza, A.; White, E.; Hutchison, B.; Schumacher, J. Application of the Iterative Ensemble Smoother Method and Cloud Computing: A Groundwater Modeling Case Study. Water 2019, 11, 1649. https://doi.org/10.3390/w11081649

AMA Style

Hayley K, Valenza A, White E, Hutchison B, Schumacher J. Application of the Iterative Ensemble Smoother Method and Cloud Computing: A Groundwater Modeling Case Study. Water. 2019; 11(8):1649. https://doi.org/10.3390/w11081649

Chicago/Turabian Style

Hayley, Kevin, Alexis Valenza, Emma White, Bruce Hutchison, and Jens Schumacher. 2019. "Application of the Iterative Ensemble Smoother Method and Cloud Computing: A Groundwater Modeling Case Study" Water 11, no. 8: 1649. https://doi.org/10.3390/w11081649

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop