An ensemble-based approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty

Coulon, Cécile; White, Jeremy T.; Pryet, Alexandre; Gatel, Laura; Lemieux, Jean-Michel

doi:https://doi.org/10.5194/hess-2023-38

Preprints

https://doi.org/10.5194/hess-2023-38

Preprints

07 Mar 2023

| 07 Mar 2023

Status: a revised version of this preprint is currently under review for the journal HESS.

An ensemble-based approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty

Cécile Coulon, Jeremy T. White, Alexandre Pryet, Laura Gatel, and Jean-Michel 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 real-world coastal aquifers, and for the incorporation of climate uncertainty into OUU approaches. An ensemble-based OUU approach was developed, considering parameter, observation and climate uncertainty, and was implemented in a real-world island aquifer in the Magdalen Islands (Quebec, Canada). A sharp-interface seawater intrusion model was developed using MODFLOW-SWI2 and a prior parameter ensemble was generated, containing multiple equally plausible realizations. Ensemble-based 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. Multi-objective 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 multi-objective, ensemble-based 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 risk-averse stances. Incorporating climate uncertainty in the OUU was also critical since it reduced the maximum allowable pumping rates for users with a risk-averse stance. The workflow used tools adapted to very high-dimensional, nonlinear models and optimization problems, to facilitate its implementation in a wide range of real-world settings.

Received: 01 Feb 2023 – Discussion started: 07 Mar 2023

Download & links

Cécile Coulon et al.

Status: final response (author comments only)

RC1:
'Comment on hess-2023-38', 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 well-organized, 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 well-written, 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 well-explained 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 500-year 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/hess-2023-38-RC1
- 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. History-matching was indeed implemented in a Bayesian framework, using the assumptions of multivariate Gaussian prior and posterior distributions, as encoded within the ensemble-based 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 model-to-measurement 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 multi-Gaussian distributions and using a prior parameter covariance matrix where the diagonal elements contain the variances of the parameters, and the off-diagonal elements contain the covariances. Note that parameters varying over several orders of magnitude were all log-transformed. It was further assumed that all parameters were statistically independent (all off-diagonal elements are null), with the exception of pilot point parameters which were spatially correlated (non-null off-diagonal 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 pilot-point 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 FOSM-based optimization under uncertainty (Coulon et al 2022) (https://doi.org/10.1029/2021WR031793). In order to compare the ensemble-based 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 PESTPP-IES 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 long-lasting.
  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 sharp-interface approximation through pairing with an advective-dispersive-based variable-density model. In the context of lateral seawater intrusion, methodologies have also been developed to optimize pumping using a coupled sharp-interface/advective-dispersive approach (e.g. Christelis and Mantoglou 2018 - https://doi.org/10.1007/s11269-018-2116-0; Dey and Prakash 2022, https://doi.org/10.1007/s11269-022-03145-w) and this could be explored in the context of freshwater lenses.
  
  Citation: https://doi.org/10.5194/hess-2023-38-AC1
RC2:
'Comment on hess-2023-38', Sreekanth Janardhanan, 27 Jun 2023

Overview:
Authors of the manuscript “An ensemble-based approach for pumping optimization in an island aquifer considering parameter, observation and climate uncertainty” applied a multi-objective optimisation under uncertainty approach available through PEST-MOU tool kit to find Pareto-optimal 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 well-written and well-organised. 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 steady-state? 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 200-year 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 Pareto-optimal 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 NSGA-II to make it computationally feasible and it may have affected the pareto-optimality? 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/hess-2023-38-RC2
- 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 steady-state? Could you please explain that in this section.
  
  Climate change predictive simulations were conducted under steady-state conditions, with boundary conditions representative of 2050 sea-level and recharge conditions. All predictive simulations were conducted in steady-state 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, steady-state 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 200-year simulation period?
  
  The pumping optimization under uncertainty was conducted under steady-state 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 PESTPP-MOU that was used, each individual in intermediate generations was mapped to the nearest individual, in decision-variable space, at which constraint probability distribution functions (PDFs) were previously evaluated, in a minimum-Euclidean-distance 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 worst-case 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 worst-case realizations could still underestimate the worst-case 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 model-independent tool for evolutionary constrained multi-objective optimization under uncertainty, Environmental Modelling & Software, 149, 10.1016/j.envsoft.2022.105316, 2022.
  
  Citation: https://doi.org/10.5194/hess-2023-38-AC2

Cécile Coulon et al.

Model code and software

swi2-ensembles Cécile Coulon, Jeremy T. White, Alexandre Pryet, Laura Gatel, and Jean-Michel Lemieux https://doi.org/10.5281/zenodo.7574457

Cécile Coulon et al.

Viewed

Total article views: 705 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
539	148	18	705	6	7

HTML: 539
PDF: 148
XML: 18
Total: 705
BibTeX: 6
EndNote: 7

Views and downloads (calculated since 07 Mar 2023)

Month	HTML	PDF	XML	Total
Mar 2023	279	52	6	337
Apr 2023	85	25	2	112
May 2023	67	19	3	89
Jun 2023	30	18	3	51
Jul 2023	53	15	3	71
Aug 2023	15	14	0	29
Sep 2023	10	5	1	16

Cumulative views and downloads (calculated since 07 Mar 2023)

Month	HTML	PDF	XML	Total
Mar 2023	279	52	6	337
Apr 2023	85	25	2	112
May 2023	67	19	3	89
Jun 2023	30	18	3	51
Jul 2023	53	15	3	71
Aug 2023	15	14	0	29
Sep 2023	10	5	1	16

Viewed (geographical distribution)

Total article views: 658 (including HTML, PDF, and XML) Thereof 658 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 11 Sep 2023

Short summary

In coastal areas, groundwater managers require information on the risk of well salinization associated with various pumping scenarios. We developed a modeling approach to identify the optimal tradeoff between groundwater pumping and probability of salinization, considering model parameter and historical observation uncertainty as well as uncertainty in sea level and recharge projections. The workflow can be implemented in a wide range of coastal settings.


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