the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An ensemblebased approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty
Cécile Coulon
Jeremy T. White
Alexandre Pryet
Laura Gatel
JeanMichel Lemieux
Abstract. In coastal zones, a major objective of groundwater management is often to determine sustainable pumping rates avoiding well salinization. Understanding how model and climate uncertainties affect optimal management solutions is essential to provide groundwater managers information about salinization risk, and is facilitated by the use of optimization under uncertainty (OUU) methods. However, guidelines are missing for the widespread implementation of OUU in realworld coastal aquifers, and for the incorporation of climate uncertainty into OUU approaches. An ensemblebased OUU approach was developed, considering parameter, observation and climate uncertainty, and was implemented in a realworld island aquifer in the Magdalen Islands (Quebec, Canada). A sharpinterface seawater intrusion model was developed using MODFLOWSWI2 and a prior parameter ensemble was generated, containing multiple equally plausible realizations. Ensemblebased history matching was conducted using an iterative ensemble smoother and yielded a posterior parameter ensemble conveying both parameter and observation uncertainty. 2050 sea level and recharge ensembles were generated and incorporated to generate a predictive parameter ensemble conveying parameter, observation and climate uncertainty. Multiobjective OUU was then conducted, aiming to both maximize pumping rates and minimize probability of well salinization. As a result, the optimal tradeoff between pumping and probability of salinization was quantified, considering parameter, historical observation and future climate uncertainty simultaneously. The multiobjective, ensemblebased OUU led to optimal pumping rates that were very different from a previous deterministic OUU, and close to the current and projected water demand for riskaverse stances. Incorporating climate uncertainty in the OUU was also critical since it reduced the maximum allowable pumping rates for users with a riskaverse stance. The workflow used tools adapted to very highdimensional, nonlinear models and optimization problems, to facilitate its implementation in a wide range of realworld settings.
 Preprint
(2028 KB)  Metadata XML
 BibTeX
 EndNote
Cécile Coulon et al.
Status: final response (author comments only)

RC1: 'Comment on hess202338', Anonymous Referee #1, 03 Apr 2023
In this paper, an approach to optimizing pumping rates on Grande Entrée Island is proposed. This approach considers uncertainty in observations, model parameters, and climate inputs. The manuscript is wellorganized, with a clear objective, and the results and conclusions provide answers to the study's objectives. Part of the methodology consists of techniques that were previously proposed and tested in other studies. I think the following comments/suggestions should be addressed:
 I suggest that the authors review the manuscript for English usage mistakes. While the article is generally wellwritten, some grammar errors are confusing, and I had to reread parts of the article to fully understand the ideas.
 In the first part of the approach, you move from a prior parameter ensemble to a posterior parameter ensemble, using observations. This seems to me like Bayesian model calibration, but you never mention the word "Bayesian" in the manuscript. Why is that? Perhaps you could explain this procedure from a Bayesian perspective.
 Following up on the previous point, is it convenient to distinguish between ensembles and pdfs? Clearly, an ensemble is always different from a PDF. Also, given that the statistical moments of your ensembles sometimes vary by an order of magnitude, what effect do you think this difference will have on the conclusions of your study?
 How do you perform the sampling? Do you assume that your parameters are mutually independent?
 Why were the initial transition zone width and longitudinal dispersivity not considered in the first part of the approach? I see that you calibrated 58 parameters, so what is the reason for not including the other two parameters?
 In the sentence that starts on line 152, you write that you account for spatial correlations between pilot points. However, I think this part is not wellexplained in the document. Are you using kriging somehow? Or does it mean that the samples you draw from a parameter x_i are conditioned on the values of a parameter x_j? Please explain the method you use to account for spatial correlation.
 On line 166, you mention that 500year simulations were carried out. Why 500? Why not 100 or 1000?
 Could you please include the mean absolute error metric in the results of Figures 6a and 6b? In Fig. 6b, the vertical axis refers to the residual (Y_sim – Y_obs). Isn't this residual large? In the caption, it is written that it is a scatter plot, so which is correct?
 What does the top horizontal axis, "Probability density," mean in Figure 10?
 In your study, you consider three sources of uncertainty (observations, model parameters, and climate forcing). However, you always rely on the same numerical model, which is not very accurate (according to Figure 6). Is the model's conceptualization another source of uncertainty? How would you consider this source of uncertainty? I don't expect you to include this analysis in the current study, which is already very comprehensive, but I wanted to highlight this source of uncertainty for your future studies.
In addition, please consider the following minor comments.
 Define the acronyms TDEM and ERT in figure 1.
 Caption of figure 1 says “The BC implemented in MODFLOW are shown” but they are not shown in the figure, they are described in the caption. Please rephrase.
 In line 334 you use the words “former” and “latter”. However, given that in the same section you speak about the results of your current study and the results of your previous study, the usage of those words can be confusing. Could you consider using other words?
Citation: https://doi.org/10.5194/hess202338RC1 
AC1: 'Reply on RC1', Cécile Coulon, 06 May 2023
We thank the anonymous reviewer for their comments, which give us the opportunity to clarify certain aspects of the manuscript. We will address all the comments in the final draft. Here is our response to the main comments.
2. Historymatching was indeed implemented in a Bayesian framework, using the assumptions of multivariate Gaussian prior and posterior distributions, as encoded within the ensemblebased data assimilation approach. We defined the prior parameter probability distribution from site characterization and expert knowledge and generated an ensemble from the prior distribution using standard sampling techniques. We then conditioned the prior ensemble with the information contained in listed observations by minimizing a modeltomeasurement fit objective function (which is inversely proportional to the likelihood function), to yield an ensemble of posterior parameter realizations. We will add this description to the manuscript.
3. Thank you for bringing up this interesting point. We agree that the statistical moments of the parameter ensemble and the PDF will generally be different, especially for small ensemble sizes. We think it is useful to present both in the manuscript, to show that although ensembles provide more reliable estimations of PDFs than FOSM analysis, they remain samples of, and therefore approximations of the PDFs. The factor that will most influence the outcome of optimization under uncertainty is whether the statistical moments of the constraint ensembles converged relative to the ensemble size, since the probabilities of constraint violation are directly used in the optimization algorithm (Eq 4, 5). We tested prior parameter ensembles with sizes ranging from 50 to 1000 realizations and found that the statistical moments of the 50% seawater salinity ensembles (mean, standard deviation, 5^{th} and 95^{th} percentiles) converged after 200 realizations.
4, 6. Random parameter fields were sampled from the prior parameter PDFs, assuming they can be described by multiGaussian distributions and using a prior parameter covariance matrix where the diagonal elements contain the variances of the parameters, and the offdiagonal elements contain the covariances. Note that parameters varying over several orders of magnitude were all logtransformed. It was further assumed that all parameters were statistically independent (all offdiagonal elements are null), with the exception of pilot point parameters which were spatially correlated (nonnull offdiagonal elements). An exponential variogram with a range equal to 3 times the pilot point spacing (i.e., 500 m) was used to describe the spatial correlation between the hydraulic conductivities at pilot point locations. We also note that as part of the pilotpoint parameterization, ordinary kriging was used to interpolate grid cell values from pilot point locations to the model cells.
5. In the same study area, the previous deterministic parameter estimation in Coulon et al (2021) (https://doi.org/10.1016/j.jhydrol.2021.126509) did not consider the initial transition zone width M and longitudinal dispersivity α_{L}, and these parameters were only considered for the FOSMbased optimization under uncertainty (Coulon et al 2022) (https://doi.org/10.1029/2021WR031793). In order to compare the ensemblebased approach to the deterministic approach, we decided to conserve this strategy. Furthermore, α_{L} and M are uninformed by the head and interface observations used in history matching (but may influence the predicted simulated salinity contours).
7. We agree that this is not clear in the manuscript. We ran the simulations until heads and interface elevations were stable close to pumping wells, and initial testing showed this was achieved under 500 years.
8. The y axis of Figure 6b should also be labeled “Simulated (m)”, thank you for finding this error. The range of simulated values is large for each observation; however, this is partly explained by the cutoff of PESTPPIES after iteration 2. Over successive IES iterations, the goodness of fit increases and the ensemble diversity (and therefore the posterior parameter variance) decreases. Using a small number of iterations is recommended when using IES, but we acknowledge that the choice of the cutoff iteration is subjective. We preferred to maintain a large ensemble diversity and possibly overestimate posterior parameter variance, rather than taking the risk of underestimating posterior parameter variance and/or risking biases in the parameter estimates arising from model error phenomena. We thought being conservative was appropriate since there are no alternative drinking water sources on the island and the consequences of well salinization can be longlasting.
9. This wasn’t clear, the top horizontal axis is the probability density function associated with the 2050 water demand.
10. We agree that the model conceptualization is an additional source of uncertainty. This is discussed in more detail in Coulon et al (2022) but should be mentioned in this paper as well. Using a sharp interface approach was a simplification of mixing processes which could result in increased conceptual uncertainty. However, the posterior parameter values display physically plausible values that are coherent with prior parameter distribution and the information in the observations was assimilated in appropriate ways; these are the two indicators available to detect the potential for conceptual model uncertainty issues. The Doherty and Christensen (2011) (https://doi.org/10.1029/2011WR010763) model pairing methodology could be used to more explicitly investigate the potential for conceptual model issues surrounding the use of the sharpinterface approximation through pairing with an advectivedispersivebased variabledensity model. In the context of lateral seawater intrusion, methodologies have also been developed to optimize pumping using a coupled sharpinterface/advectivedispersive approach (e.g. Christelis and Mantoglou 2018  https://doi.org/10.1007/s1126901821160; Dey and Prakash 2022, https://doi.org/10.1007/s1126902203145w) and this could be explored in the context of freshwater lenses.
Citation: https://doi.org/10.5194/hess202338AC1

RC2: 'Comment on hess202338', Sreekanth Janardhanan, 27 Jun 2023
Overview:
Authors of the manuscript “An ensemblebased approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty” applied a multiobjective optimisation under uncertainty approach available through PESTMOU tool kit to find Paretooptimal solutions for management of sea water intrusion in a coastal island considering maximisation of pumping and reliability of meeting the salinity constraint as conflicting objectives. The manuscript is wellwritten and wellorganised. It is a contribution of interest to readers of HESS and advances the science in methodologies for optimising groundwater management. I have only a few minor comments and hence recommend accepting the manuscript after minor revisions.
Individual comments:
Section 3.3: It is not clear if the climate change predictive runs were also steadystate? Could you please explain that in this section.
Line 210: It is not very clear how values are resampled. Could you please explain.
Line 215: Ah OK. So the are uncorrelated with current recharge – does this mean were randomly added to realisations in current recharge ensemble ?
Line 245: :500 year initial simulations with zero pumping” – Is this for reaching a new equilibrium corresponding to the new recharge and other factors corresponding to future climate? I assume this steady state can be considerably different from the historical steady state with recharge and sea level changes?
Line 250: What is the significance of 200year simulation period?
Line 265: “Prediction ensemble was reevaluated every 10 generations and reused in intermediate generations” – how was objective function calculated in the intermediate generations?
Figure 7: Figure 7 shows considerable correlation between recharge and K values. Perhaps discuss the implication of this in predictive simulations with a different recharge regime (although noting that the mean doesn’t change much) that does not consider correlations with historical recharge. Is it likely to bias predictions?
Figures 10 and 11: Looks like the final convergence of the Paretooptimal front hasn’t been achieved. There are several solutions that are dominated by other solutions for both objective functions. I assume you had to optimise the number of iterations (150) and population (30) of NSGAII to make it computationally feasible and it may have affected the paretooptimality? From a practical point of view, these still give valuable solutions.
Line 375: Ah, I see the answer to the above question explained here.
line 420 – 425: Yes, stack ordering approaches may relieve the computational burden. Especially in this case, if you pick realisations from the tail end of recharges for worst climate scenarios it would do the job?Citation: https://doi.org/10.5194/hess202338RC2 
AC2: 'Reply on RC2', Cécile Coulon, 21 Jul 2023
We thank Dr Sreekanth Janardhanan for their comments, which we will address in the final draft. Here is our response to the comments:
 Section 3.3: It is not clear if the climate change predictive runs were also steadystate? Could you please explain that in this section.
Climate change predictive simulations were conducted under steadystate conditions, with boundary conditions representative of 2050 sealevel and recharge conditions. All predictive simulations were conducted in steadystate conditions because the storage parameters were unconstrained by the history matching, and therefore remained highly uncertain. In reality, since climate change effects are a transient process and this steady state may never be reached, this approach can be viewed as being conservative.
 Line 210: It is not very clear how values are resampled. Could you please explain. Line 215: Ah OK. So the are uncorrelated with current recharge – does this mean were randomly added to realisations in current recharge ensemble ?
The recharge perturbations ΔR were indeed uncorrelated with the current recharge values R_{current, MODFLOW}. The ΔR and R_{current, MODFLOW} realizations were randomly paired together to generate the R_{2050, MODFLOW ensemble}.
 Line 245: :500 year initial simulations with zero pumping” – Is this for reaching a new equilibrium corresponding to the new recharge and other factors corresponding to future climate? I assume this steady state can be considerably different from the historical steady state with recharge and sea level changes?
We will justify the length of the simulations in the manuscript, as this question was also raised by Reviewer 1. Initial, steadystate simulations with zero pumping were run until heads and interface elevations were stable close to pumping wells. Initial testing showed this was achieved within 500 years. As shown in Figure 9, there is variability both in the steady states obtained with the posterior ensemble and the 2050 predictive ensemble. For the interface elevation under pumping wells, the mean steady state value is similar in both cases, but there is twice as much variability in the initial steady states when climate change is accounted for.
 Line 250: What is the significance of 200year simulation period?
The pumping optimization under uncertainty was conducted under steadystate conditions, using long transient simulations to allow the freshwater lens to reach a new steady state under the pumping rates tested. This new steady state was found to be achieved within 200 years.
 Line 265: “Prediction ensemble was reevaluated every 10 generations and reused in intermediate generations” – how was objective function calculated in the intermediate generations?
In the version 5.1.24 of PESTPPMOU that was used, each individual in intermediate generations was mapped to the nearest individual, in decisionvariable space, at which constraint probability distribution functions (PDFs) were previously evaluated, in a minimumEuclideandistance sense. The constraint PDFs of the latter were translated to the former, by differencing the simulated constraint values between the two individuals, assuming these values represented the mean of the PDFs. This approach assumes that individuals close to each other in decision variable space have similar constraint PDFs (PEST++ Development Team, 2022; White et al., 2022).
 Figure 7: Figure 7 shows considerable correlation between recharge and K values. Perhaps discuss the implication of this in predictive simulations with a different recharge regime (although noting that the mean doesn’t change much) that does not consider correlations with historical recharge. Is it likely to bias predictions?
Parameter correlations can be identified during history matching. Recharge (R) and hydraulic conductivity (K) parameters are known to be correlated (Anderson et al., 2015), which can result in realizations where high current R is paired with high K values or low current R is paired with low K values. The current R values stayed paired with their corresponding K field and recharge perturbations ΔR were randomly added to them. However, no correlation was assumed between current and future recharge. With this assumption, high future R could be paired with low K values, and low future R with high K values; therefore, the tails of the constraint PDFs can be explored more thoroughly. While this assumption might overestimate the constraint uncertainty, it can be viewed as being conservative. Furthermore, bias in future recharge predictions could be caused by having excessive confidence in the R/K correlation learned during history matching.
 line 420 – 425: Yes, stack ordering approaches may relieve the computational burden. Especially in this case, if you pick realisations from the tail end of recharges for worst climate scenarios it would do the job?
We agree that worstcase scenarios could be explored by running optimizations on the realizations with the lowest historical recharge values associated with the most extreme recharge decrease scenarios (i.e., the most important ΔR perturbations). However, realizations remain samples of PDFs, and the ΔR realizations randomly sampled from the ΔR probability distribution function may not be representative of its extreme percentiles. Therefore, these worstcase realizations could still underestimate the worstcase scenario.
References:
Anderson, M.P., Woessner, W.W., Hunt, R.J.,: Applied groundwater modeling: simulation of flow and advective transport. Academic press. 2015
PEST++ Development Team: PEST++ Manual, Version 5.1.23, 2022.
White, J. T., Knowling, M. J., Fienen, M. N., Siade, A., Rea, O., and Martinez, G.: A modelindependent tool for evolutionary constrained multiobjective optimization under uncertainty, Environmental Modelling & Software, 149, 10.1016/j.envsoft.2022.105316, 2022.Citation: https://doi.org/10.5194/hess202338AC2

AC2: 'Reply on RC2', Cécile Coulon, 21 Jul 2023
Cécile Coulon et al.
Model code and software
swi2ensembles Cécile Coulon, Jeremy T. White, Alexandre Pryet, Laura Gatel, and JeanMichel Lemieux https://doi.org/10.5281/zenodo.7574457
Cécile Coulon et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

539  148  18  705  6  7 
 HTML: 539
 PDF: 148
 XML: 18
 Total: 705
 BibTeX: 6
 EndNote: 7
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1