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.

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.