Groundwater withdrawal in randomly heterogeneous coastal aquifers

We analyze the combined effects of aquifer heterogeneity and pumping operations on seawater intrusion (SWI), a phenomenon which is threatening coastal aquifers worldwide. Our investigation is set within a probabilistic framework and relies on a numerical Monte Carlo approach targeting transient variable-density flow and solute transport in a three-dimensional randomly heterogeneous porous domain. The geological setting is patterned after the Argentona river basin, in the Maresme region of Catalonia (Spain). Our numerical study is concerned with exploring the effects of (a) random heterogeneity of the domain on SWI in combination with (b) a variety of groundwater withdrawal schemes. The latter have been designed by varying the screen location along the vertical direction and the distance of the wellbore from the coastline and from the location of the freshwater–saltwater mixing zone which is in place prior to pumping. For each random realization of the aquifer permeability field and for each pumping scheme, a quantitative depiction of SWI phenomena is inferred from an original set of metrics characterizing (a) the inland penetration of the saltwater wedge and (b) the width of the mixing zone across the whole three-dimensional system. Our results indicate that the stochastic nature of the system heterogeneity significantly affects the statistical description of the main features of the seawater wedge either in the presence or in the absence of pumping, yielding a general reduction of toe penetration and an increase of the width of the mixing zone. Simultaneous extraction of fresh and saltwater from two screens along the same wellbore located, prior to pumping, within the freshwater–saltwater mixing zone is effective in limiting SWI in the context of groundwater resources exploitation.

1 Introduction

Groundwater resources in coastal aquifers are seriously threatened by seawater intrusion (SWI), which can deteriorate the quality of freshwater aquifers, thus limiting their potential use. This situation is particularly exacerbated within areas with intense anthropogenic activities, which are associated with competitive uses of groundwater in connection, e.g., with agricultural processes, industrial processes and/or high urban water supply demand. Critical SWI scenarios are attained when seawater reaches extraction wells designed for urban freshwater supply, with severe environmental, social and economic implications (e.g., Custodio, 2010; Mas-Pla et al., 2014; Mazi et al., 2014).

The development of effective strategies for sustainable use of groundwater resources in coastal regions should be based on a comprehensive understanding of SWI phenomena. This challenging problem has been originally studied by assuming a static equilibrium between freshwater (FW) and seawater (SW) and a sharp FW–SW interface, where FW and SW are considered as immiscible fluids. Under these hypotheses, the vertical position of the FW–SW interface below the sea level, z, is given by the Ghyben–Herzberg solution, $z=h_{\mathrm{F}}{\mathit{\rho }}_{\mathrm{F}}/\mathrm{\Delta }\mathit{\rho }$, where h_F is the FW head above sea level and $\mathrm{\Delta }\mathit{\rho }={\mathit{\rho }}_{\mathrm{S}}-{\mathit{\rho }}_{\mathrm{F}}$ is the density contrast, ρ_S and ρ_F, respectively, being SW and FW density. Starting from these works, several analytical and semi-analytical expressions have been developed to describe SWI under diverse flow configurations (e.g., Strack, 1976; Dagan and Zeitoun, 1998; Bruggeman, 1999; Cheng et al., 2000; Bakker, 2006; Nordbotten and Celia, 2006; Park et al., 2009). In this broad context, Strack (1976) evaluated the maximum (or critical) pumping rate to avoid encroachment of SW in FW pumping wells. All of these sharp-interface-based solutions neglect a key aspect of SWI phenomena, i.e., the formation of a transition zone where mixing between FW and SW takes place and fluid density varies with salt concentration. As a consequence, these models may significantly overestimate the actual penetration length of the SW wedge, leading to an excessively conservative evaluation of the critical pumping rate (Gingerich and Voss, 2005; Pool and Carrera, 2011; Llopis-Albert and Pulido-Velazquez, 2013).

A more realistic approach relies on the formulation and solution of a variable-density problem, in which SW and FW are considered as miscible fluids and groundwater density depends on salt concentration. The complexity of the problem typically prevents finding a solution via analytical or semi-analytical methods, with a few notable exceptions (Henry, 1964; Dentz et al., 2006; Bolster et al., 2007; Zidane et al., 2012; Fahs et al., 2014). Henry (1964) developed a semi-analytical solution for a variable-density diffusion problem in a (vertical) two-dimensional homogeneous and isotropic domain. Dentz et al. (2006) and Bolster et al. (2007) applied perturbation techniques to solve analytically the Henry problem for a range of (small and intermediate) values of the Péclet number, which characterizes the relative strength of convective and dispersive transport mechanisms. Zidane et al. (2012) solved the Henry problem for realistic (small) values of the diffusion coefficient. Finally, Fahs et al. (2014) presented a semi-analytical solution for a square porous cavity system subject to diverse salt concentrations at its vertical walls. Practical applicability of these solutions is quite limited, due to their markedly simplified characteristics. Various numerical codes have been proposed to solve variable-density flow and transport equations (e.g., Voss and Provost, 2002; Ackerer et al., 2004; Soto Meca et al., 2007; Ackerer and Younes, 2008; Albets-Chico and Kassinos, 2013). Numerical simulations can provide valuable insights into the effects of dispersion on SWI, a feature that is typically neglected in analytical and semi-analytical solutions. Abarca et al. (2007) introduced a modified Henry problem to account for dispersive solute transport and anisotropy in hydraulic conductivity. Kerrou and Renard (2010) analyzed the dispersive Henry problem within two- and three-dimensional randomly heterogeneous aquifer systems. These authors relied on computational analyses performed on a single three-dimensional realization, invoking ergodicity assumptions. Lu et al. (2013) performed a set of laboratory experiments and numerical simulations to investigate the effect of geological stratification on SW–FW mixing. Riva et al. (2015) considered the same setting as in Abarca et al. (2007) and studied the way quantification of uncertainty associated with SWI features is influenced by lack of knowledge of four key dimensionless parameters controlling the process, i.e., gravity number, permeability anisotropy ratio and transverse and longitudinal Péclet numbers. Enhancement of mixing in the presence of tidal fluctuations and/or FW table oscillations has been analyzed by Ataie-Ashtiani et al. (1999), Lu et al. (2009, 2015) and Pool et al. (2014) in homogeneous aquifers and by Pool et al. (2015) in randomly heterogeneous three-dimensional systems (under ergodic conditions). Recent reviews on the topic are offered by Werner et al. (2013) and Ketabchi et al. (2016).

Several numerical modeling studies have been performed with the aim of identifying the most effective strategy for the exploitation of groundwater resources in domains mimicking the behavior of specific sites. Dausman and Langevin (2005) examined the influence of hydrologic stresses and water management practices on SWI in a superficial aquifer (Broward County, USA) by developing a variable-density model formed by two homogeneous hydrogeological units. Werner and Gallagher (2006) studied SWI in the coastal aquifer of the Pioneer Valley (Australia), illustrating the advantages of combining hydrogeological and hydrochemical analyses to understand salinization processes. Misut and Voss (2007) analyzed the impact of aquifer storage and recovery practices on the transition zone associated with the salt water wedge in the New York City aquifer, which was modeled as a perfectly stratified system. Cobaner et al. (2012) studied the effect of transient pumping rates from multiple wells on SWI in the Gosku deltaic plain (Turkey) by means of a three-dimensional heterogeneous model, calibrated using head and salinity data.

All the aforementioned field-scale contributions are framed within a deterministic approach, where the system attributes (e.g., permeability) are known (or determined via an inverse modeling procedure), so that the impact of uncertainty of hydrogeological properties on target environmental (or engineering) performance metrics is not considered. Exclusive reliance on a deterministic approach is in stark contrast with the widely documented and recognized observation that a complete knowledge of aquifer properties is unfeasible. This is due to a number of reasons, including observation uncertainties and data availability, i.e., available data are most often too scarce or too sparse to yield an accurate depiction of the subsurface system in all of its relevant details. Stochastic approaches enable us not only to provide predictions (in terms of best estimates) of quantities of interest, but also to quantify the uncertainty associated with such predictions. The latter can then be transferred, for example, into probabilistic risk assessment, management and protection protocols for environmental systems and water resources.

Only a few contributions studying SWI within a stochastic framework have been published to date. The vast majority of these works considers idealized synthetic showcases and/or simplified systems. To the best of our knowledge, only two studies (Lecca and Cau, 2009; Kerrou et al., 2013) have analyzed the transient behavior of a realistic costal aquifer within a probabilistic framework. Lecca and Cau (2009) evaluated SWI in the Oristano (Italy) aquifer by considering a stratified system where the aquitard is characterized by random heterogeneity. Kerrou et al. (2013) analyzed the effects of uncertainty in permeability and distribution of pumping rates on SWI in the Korba aquifer (Tunisia). Both works characterize the uncertainty of SWI phenomena in terms of the planar extent of the system characterized by a target probability of exceedance of a given threshold concentration.

Our stochastic numerical study has been designed to mimic the general behavior of the Argentona aquifer, in the Maresme region of Catalonia (Spain). This area, as well as other Mediterranean deltaic sites, is particularly vulnerable to SWI (Custodio, 2010). We note that the objective of this study is not the quantification of the SWI dynamics of a particular field site. Our emphasis is on the analysis of the impact of random aquifer heterogeneity and withdrawals on SWI patterns in a realistic scenario. Key elements of novelty in our work include the introduction and the study of a new set of metrics, aimed at investigating quantitatively the effects of random heterogeneity on the three-dimensional extent of the SW wedge penetration and of the FW–SW mixing zone. The impact of random heterogeneity is assessed in combination with three diverse withdrawal scenarios, designed by varying (i) the distance of the wellbore from the coastline and from the FW–SW mixing zone and (ii) the screen location along the vertical direction. We frame our analysis within a numerical Monte Carlo (MC) approach.

Figure 1 (a) Location of the Argentona river basin. (b) Two-dimensional, constant-density model (modified from Rodriguez Fernandez, 2015). Pumping wells (brown dots) and three-dimensional model location (blue shaded area) are also depicted. (c) Planar view of the three-dimensional variable-density flow and transport model. The unstructured grid and diverse recharge areas are also shown.

Section 2 provides a general description of the field site, the mathematical model adopted to simulate flow and transport phenomena in three-dimensional heterogeneous systems and the numerical settings. Section 3 illustrates the key results of our work and Sect. 4 contains our concluding remarks.

2 Materials and methods

2.1 Site description

To consider a realistic scenario that is relevant to SWI problems in highly exploited aquifers, we cast our analysis in a setting inspired from the Argentona river basin, located in the region of Maresme, in Catalonia, Spain (see Fig. 1a). This aquifer is a typical deltaic site, characterized by shallow sedimentary units and a flat topography. As such, it has a strategic value for anthropogenic (including agricultural, industrial and touristic) activities. The geological formation hosting the groundwater resource is mainly composed of a granitic Permian unit. A secondary unit of quaternary sediments is concentrated along the Argentona river.

Rodriguez Fernandez (2015) developed a conceptual and numerical model to simulate transient two-dimensional (horizontal) constant-density flow in the Argentona river basin across the area of about 35 km² depicted in Fig. 1b. The aquifer is heavily exploited, as shown by the large number of wells included in Fig. 1b, mainly located along the Argentona river. The latter is a torrential ephemeral stream, in which water flows only after heavy rain events. On these bases, it is assumed that the only water intake to the river comes from surface runoff of both granitic and quaternary units. On the basis of transient hydraulic head measurements available at 21 observation wells (not shown in Fig. 1b) in the period 2006–2013, Rodriguez Fernandez (2015) characterized the site through a uniform permeability, whose value was estimated as k_B = 1.77 × 10⁻¹¹ m², and provided estimates of temporally and spatially variable recharge rates, according to land use or cover (see Fig. 1c).

Our numerical analysis focuses on the coastal portion of the Argentona basin (Fig. 1c). This region extends for about 2.5 km along the coast (i.e., along the full width of the basin) and up to 750 m inland from the coast. The offshore extension of the aquifer is not considered. The size of the study area has been selected on the basis of the results of preliminary simulations on larger domains, which highlighted that salt concentration values were appreciable only within a narrow (less than 400 m wide) region close to the coast. The vertical thickness of the domain ranges from 50 m along the coast up to 60 m at the inland boundary, the underlying clay sequence being considered to represent the impermeable bottom of the aquifer. SWI is simulated by means of a three-dimensional variable-density flow and solute transport model based on the well-established finite element USGS SUTRA code (Voss and Provost, 2002) over the 8-year time window 2006–2013. Details of the mathematical and numerical model are discussed in the following sections.

2.2 Flow and transport equations

Fluid flow is governed by mass conservation and Darcy's law:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial \left(\mathit{\varphi }\mathit{\rho }\right)}{\partial t}}+{\mathrm{\nabla }}^{\mathrm{T}}\left(\mathit{\rho }\mathit{q}\right)=\mathrm{0}{\displaystyle }& {\displaystyle \frac{\partial \left(\mathit{\varphi }\mathit{\rho }\right)}{\partial t}}={\displaystyle \frac{\mathrm{1}}{g}}{S}_{\mathrm{S}}{\displaystyle \frac{\partial p}{\partial t}}+\mathit{\varphi }{\displaystyle \frac{\partial \mathit{\rho }}{\partial C}}{\displaystyle \frac{\partial C}{\partial t}},\\ \text{(2)}& {\displaystyle }& {\displaystyle }\mathit{q}=-{\displaystyle \frac{\mathit{k}}{\mathit{\mu }}}\left(\mathrm{\nabla }p+\mathit{\rho }g\mathrm{\nabla }z\right),\end{array}$$

where q (L T⁻¹) is the specific discharge vector, with components q_x, q_y and q_z, respectively, along x, y (see Fig. 1c) and z axes (z denoting the vertical direction); C (–) is solute concentration, or solute mass fraction, expressed as mass of solute per unit mass of fluid; k (L²) is the diagonal permeability tensor, with components k₁₁ = k_x, k₂₂ = k_y and k₃₃ = k_z, respectively, along directions x, y and z; ϕ (–) is aquifer porosity; ρ (M L⁻³) and μ (M L⁻¹ T⁻¹), respectively, are fluid density and dynamic viscosity; p (M L⁻¹ T⁻²) is fluid pressure; g (L T⁻²) is gravity; S_S (L⁻¹) is the specific storage coefficient, and the superscript T denotes transpose. Equations (1) and (2) must be solved jointly with the advection–dispersion equation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \left(\mathit{\varphi }\mathit{\rho }C\right)}{\partial t}}+{\mathrm{\nabla }}^{\mathrm{T}}\left(\mathit{\rho }C\mathit{q}\right)-{\mathrm{\nabla }}^{\mathrm{T}}\left[\mathit{\rho }\mathbf{D}\mathrm{\nabla }C\right]=\mathrm{0}.\end{array}$$

D (L² T⁻¹) is the dispersion tensor defined as

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathbf{D}=\left(\mathit{\varphi }{D}_{\mathrm{m}}+{\mathit{\alpha }}_{\mathrm{T}}\left|\mathit{q}\right|\right)\mathbf{I}+\left({\mathit{\alpha }}_{\mathrm{L}}-{\mathit{\alpha }}_{\mathrm{T}}\right){\displaystyle \frac{\mathit{q}{\mathit{q}}^{\mathrm{T}}}{\left|\mathit{q}\right|}},\end{array}$$

where D_m (L² T⁻¹) is the molecular diffusion coefficient, and α_L and α_T (L), respectively, are the longitudinal and transverse dispersivity coefficients. Closure of the system of Eqs. (1)–(4) requires a relationship between fluid properties (ρ and μ) and solute concentration C. Fluid viscosity can be assumed constant in typical SWI settings and the following model has been shown to be accurate at describing the evolution of ρ with C (Kolditz et al., 1998):

$$\begin{array}{}\text{(5)}& \mathit{\rho }={\mathit{\rho }}_{\mathrm{F}}+\left({\mathit{\rho }}_{\mathrm{S}}-{\mathit{\rho }}_{\mathrm{F}}\right){\displaystyle \frac{C}{C_{\mathrm{S}}}},\end{array}$$

where C_S is SW concentration. Table 1 lists aquifer and fluid parameters adopted.

2.3 Numerical model

The domain depicted in Fig. 1c is discretized through an unstructured three-dimensional grid formed by 101 632 hexahedral elements. The resolution of the mesh increases towards the sea, where our analysis requires the highest spatial detail. The element size along both horizontal directions x and y ranges from 10 m (in a 200 m wide region along the coast) to 60 m (close to the inland boundary). Element size along the vertical direction is 2 and 4 m, respectively, within the first 10 m from the top surface and across the remaining 40 m. These choices are consistent with constraints for numerical stability, namely Λ_L≤4α_L, Λ_L being the distance between element sides measured along a flow line and α_L=5 m. In the absence of transverse dispersivity estimates and in analogy with main findings of previous works (e.g., Cobaner et al., 2012), we set ${\mathit{\alpha }}_{\mathrm{T}}={\mathit{\alpha }}_{\mathrm{L}}/\mathrm{10}$.

Table 1 Parameters adopted in the numerical model.

A sequential Gaussian simulation algorithm (Deutsch and Journel, 1998) is employed to generate unconditional log permeability fields, Y(x)=ln k(x), characterized by a given mean 〈Y〉=ln k_B and variogram structure. Since no Y data are available, for the purpose of our simulations we assume that the spatial structure of Y is described by a spherical variogram, with moderate variance, i.e., ${\mathit{\sigma }}_Y^{\mathrm{2}}=\mathrm{1.0}$, and isotropic correlation scale λ=100 m. This value is consistent with the observation that the integral scale of log conductivity and transmissivity values inferred worldwide using traditional (such as exponential and spherical) variograms tends to increase with the length scale of the sampling window at a rate of about 1∕10 (Neuman et al., 2008). The MC approach allows us to provide a probabilistic description of the SWI processes associated with space–time distributions of flux and solute concentration obtained in each random realization of the permeability field by solving the transient flow and transport Eqs. (1)–(5). Note also that λ>5Δ in our model, Δ being the largest cell size within a 200 m wide region along the coast. A model satisfying this condition has the benefit of (a) yielding a proper reconstruction of the spatial correlation structure of Y and (b) limiting the occurrence of excessive variations between values of aquifer properties across neighboring cells (Kerrou and Renard, 2010) in the region where SWI takes place.

Equations (1)–(5) are solved jointly, adopting the following boundary conditions. The lateral boundaries perpendicular to the coastline and the base of the aquifer are impervious; along the inland boundary we set a time-dependent prescribed hydraulic head, h_F (with $h_{\mathrm{F}}=z+p/({\mathit{\rho }}_{\mathrm{F}}g$) ranging between 1.2 and 3.5 m), and C=C_F; along the coastal boundary we impose a prescribed head, h_S=0 (with $h_{\mathrm{S}}=z+p/\left({\mathit{\rho }}_{\mathrm{S}}g\right)$), and

$$\begin{array}{}\text{(6)}& {\displaystyle }\left(\mathit{q}C-\mathbf{D}\mathrm{\nabla }C\right)\cdot \mathit{n}=\left\{\begin{array}{ll}\mathit{q}\cdot \mathbf{n}C& \text{if}\mathit{q}\cdot \mathit{n}<\mathrm{0}\\ \mathit{q}\cdot \mathit{n}{C}_{\mathrm{S}}& \text{if}\mathit{q}\cdot \mathit{n}>\mathrm{0}\end{array}\right.,\end{array}$$

n being the normal vector pointing inward along the boundary; i.e., the concentration of the water entering and leaving the system across the coastal boundary is, respectively, equal to C_S and C. The area of interest is characterized by five recharge zones, identified on the basis of the land use (inferred from the SIGPAC, 2005, dataset). The total recharge varies slightly in time (on a monthly basis), with mean value equal to 7.6 L s⁻¹ (see Fig. 1c). Initial conditions are set as h=0 and C=0, as commonly assumed in the literature (e.g., Jakovovic et al., 2016) and flow and transport are simulated for an 8-year period (2006–2013) with a uniform time step Δt=1 day. Varying the initial conditions (e.g., setting, h=2.4 m, corresponding to the mean value of h set along the inland boundary), did not lead to significantly different results at the end of the simulated 8-year period, t=t_w (details not shown).

A pumping well is activated at time t_w and flow and transport are simulated for 8 additional months. We consider three diverse pumping schemes, reflecting realistic engineered operational settings. The borehole is located along the vertical cross-section B–B^′, as shown in Fig. 2a. Scheme 1 (S1) and Scheme 2 (S2) consider the pumping well to be placed 180 m away from the shoreline, i.e., at y^′ = y∕λ = 1.8, outside of the mixing zone which is in place prior to pumping in all MC realizations. In S1 (Fig. 2b) the well is screened in the upper part of the aquifer, the screen starting from 2 m below the top and extending across a total thickness of 15 m. The double-negative hydraulic barrier is a technology employed to manage (or possibly prevent) SWI. It relies on the joint use of an FW (production well) and an SW (scavenger well) pumping well. We employ this setting in S2 (Fig. 2c), which is designed by adding an additional screen in the lower part of the aquifer to S1, starting from 34 m below the ground surface and extending across a total thickness of 12 m. In S2, production and scavenger wells are located along the same borehole, complying with technical and economic requirements typically associated with field applications (Pool and Carrera, 2009; Saravanan et al., 2014; Mas-Pla et al., 2014). This scheme is also particularly appealing when applied to renewable energy resources, such as the pressure retarded osmosis, which allows conversion of the chemical energy of two fluids (FW and SW) into mechanical and electrical energy (Panyor, 2006). Scheme 3 (S3, Fig. 2d) shares the same operational design of S2, but the borehole is moved seawards, at y^′ = 0.8, i.e., within the mixing zone in place before pumping. In all pumping schemes, the withdrawal rate is constant in time and uniformly distributed along the well screens, where pumping is implemented by setting a flux-type condition. We set a total withdrawal rate Q = 5 L s⁻¹ at the upper well screen in S1, S2 and S3; an additional extraction at the same rate Q is imposed along the lower screen in S2 and S3. During the pumping period, transport is simulated using a time step Δt = 30 min for the first month and with increasing Δt progressively up to a maximum value Δt = 120 min for the remaining 7 months, as the system showed progressively smoother variations while approaching steady state. In Sect. 3 all schemes are compared against a benchmark case, Scheme 0 (S0), where no pumping wells are active.

Figure 2 (a) Planar view of the location of the boreholes. Vertical cross section B–B^′ and well-screen location for scheme (b) S1, (c) S2 and (d) S3.

Figure 3 Graphical representation of the SWI metrics defined in Sect. 2.4.

2.4 Local and global SWI metrics

To provide a comprehensive characterization of SWI phenomena across the whole three-dimensional domain, we quantify the extent of the SW wedge and of the associated mixing zone on the basis of seven metrics, as illustrated in Fig. 3.

For each MC realization and along each vertical cross-section perpendicular to the coast, we evaluate (i) toe penetration, $L_{\mathrm{T}}^{\prime }$, measured as the distance from the coast of the iso-concentration curve C∕C_S = 0.5 at the bottom of the aquifer; (ii) solute spreading at the toe, $L_{\mathrm{S}}^{\prime }$, evaluated as the separation distance along the direction perpendicular to the coast (y axis in Fig. 1c) between the iso-concentration curves C∕C_S = 0.25 and C∕C_S = 0.75 at the bottom of the aquifer; (iii) mean width of the mixing zone, $W_{\mathrm{MZ}}^{\prime }$, i.e., the spatially averaged vertical distance between iso-concentration curves C∕C_S = 0.25 and C∕C_S = 0.75 within the region $\mathrm{0.2}{L}_{\mathrm{T}}^{\prime }$ ≤ y^′ ≤ $\mathrm{0.8}{L}_{\mathrm{T}}^{\prime }$. All quantities $L_{\mathrm{T}}^{\prime }$, $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$ are dimensionless, corresponding to distances rescaled by λ, the correlation scale of Y. We also analyze the (dimensionless) areal extent of SW penetration and solute spreading at the bottom of the aquifer by evaluating (iv) $A_{\mathrm{T}}^{\prime }$ and (v) $A_{\mathrm{S}}^{\prime }$, as obtained by integrating $L_{\mathrm{T}}^{\prime }$ and $L_{\mathrm{S}}^{\prime }$ along the coastline. Finally, we quantify SWI across the whole thickness of the aquifer by evaluating the dimensionless volumes enclosed between (vi) the sea boundary and the iso-concentration surface C∕C_S = 0.5, here termed $V_{\mathrm{T}}^{\prime }$, and (vii) the isosurfaces C∕C_S = 0.25 and C∕C_S = 0.75, here termed $V_{\mathrm{S}}^{\prime }$. First- and second-order statistical moments of each of these metrics are evaluated within the adopted MC framework. With the settings employed here, a flow and transport simulation on a single system realization is associated with a computational cost of about 12 h (on an Intel Core^™ i7-6900K CPU@ 3.20 GHz processor). Results illustrated in the following sections are based on a suite of 400 MC simulations (i.e., 100 MC simulation for each scenario, S0–S3). Details concerning the analysis of convergence of the first and second statistical moments of the metrics here introduced are provided in Appendix A.

Figure 4 Permeability distribution (color contour plots), streamlines (black dashed curves) and iso-concentration curves C∕C_S = 0.25, 0.5 and 0.75 (solid red curves) along the cross-section B–B^′ evaluated at t=t_w within (b) the equivalent homogeneous aquifer and (c) one random realization of k. Vertical exaggeration = 5.

3 Results and discussion

3.1 Effects of three-dimensional heterogeneity on SWI

The effects of heterogeneity on SWI are inferred by comparing the results of our MC simulations against those obtained for an equivalent homogeneous aquifer, characterized by an effective permeability, k_ef. The latter is here evaluated as $k_{\mathrm{ef}}=e^{〈Y〉+{\mathit{\sigma }}_Y^{\mathrm{2}}/\mathrm{6}}=\mathrm{2.09}\times {\mathrm{10}}^{-\mathrm{11}}$ m² (Ababou, 1996). A first clear effect of heterogeneity is noted on the structure of the three-dimensional flow field. This is elucidated in Fig. 4, where we depict the permeability color map and streamlines (dashed curves) obtained at the end of the 8-year simulation period along the vertical cross-section B–B^′ perpendicular to the coastline. Note that henceforth, results are depicted in terms of dimensionless spatial coordinates, y^′ = y∕λ and z^′ = z∕λ. Flow within the equivalent homogeneous aquifer (Fig. 4b) is essentially horizontal at locations far from the zone where SWI phenomena occur. Otherwise, the vertical flux component, q_z, is non-negligible throughout the whole domain for the heterogeneous system (Fig. 4c) and streamlines tend to focus towards regions characterized by large permeability values. Figure 4 also depicts iso-concentration curves C∕C_S = 0.25, 0.5 and 0.75 within the transition zone (red curves). It can be noted that iso-concentration profiles tend to be sub-parallel to streamlines directed towards the seaside boundary. In the homogeneous domain the slope of these curves varies mildly and in a gradual manner from the top to the bottom of the aquifer. Their slope in the heterogeneous domain is irregular and markedly influenced by the spatial arrangement of permeability. In other words, the way streamlines are refracted at the boundary between two blocks of contrasting permeability drives the local pattern of concentration contour lines. As a consequence, iso-concentration curves tend to become sub-vertical when solute is transitioning from regions characterized by high k values to zones associated with moderate to small k, a sub-horizontal pattern being observed when transitioning from low to high k values.

Figure 5 depicts isolines C∕C_S = 0.5 at the bottom of the aquifer (Fig. 5a) and along three vertical cross-sections, selected to exemplify the general pattern observed in the system (Fig. 5b–d), evaluated for (i) the 100 heterogeneous realizations analyzed (dotted blue curves), (ii) the equivalent homogeneous system (denoted as Hom; solid red curve) and (iii) the configuration obtained by averaging the concentration fields across all MC heterogeneous realizations (denoted as Ens; solid blue curve). These results suggest that iso-concentration curves exhibit considerably large spatial variations within a single realization. Comparison of the results obtained for Ens and Hom reveals that the mean wedge penetration at the bottom layer is slightly overestimated by the solution computed within the equivalent homogeneous system. Otherwise, along the coastal vertical boundary, the extent of the area with (mean) relative concentration larger than 0.5 is generally underestimated by Hom. These outcomes (hereafter called rotation effects) are consistent with previous literature findings associated with (a) deterministic models (Abarca et al., 2007) or (b) stochastic analyses performed under ergodicity assumptions (Kerrou and Renard, 2010, and Pool et al., 2015). Abarca et al. (2007) showed that a similar effect can be observed in a homogenous domain when considering increasing values of the dispersion coefficient. Kerrou and Renard (2010) and Pool et al. (2015) noted that the strength of the rotation effect increases with the variance of the log permeability field.

Figure 5 MC-based (dotted blue curves) iso-concentration lines C∕C_S = 0.5 at t=t_w along (a) the bottom of the aquifer and (b–d) three vertical cross-sections perpendicular to the coast. The ensemble-average concentration curves (solid blue curve) and the results obtained within the equivalent homogeneous system (solid red curve) are also reported. Vertical exaggeration = 5.

Figure 6 Metrics ${\mathit{\xi }}^{\prime \prime }={\mathit{\xi }}^{\prime }/{\mathit{\xi }}^{\prime \mathrm{Hom}}$ evaluated at t=t_w within each MC heterogeneous realization (dots) and on the basis of the ensemble-average concentration field (blue line). Ensemble mean metrics (solid black line) and confidence intervals of width equal to ± 1 standard deviation of ξ^′′ about their mean (dashed black lines) are also depicted. Metrics in (a–c) are computed along the vertical cross-section B–B^′.

In the following, we analyze values of ${\mathit{\xi }}^{\prime \prime }={\mathit{\xi }}^{\prime }/{\mathit{\xi }}^{\prime \mathrm{Hom}}$, ξ^′ representing a given SWI metric (as listed in Sect. 2.4) and ξ^′Hom being the corresponding metric obtained on the equivalent homogeneous system. Figure 6 depicts the range of variability of ξ^′′ across the set of MC realizations (symbols). Each plot is complemented by the depiction of (i) the (ensemble) average of ξ^′′, as evaluated through all MC realizations, 〈ξ^′′〉 (black solid line); (ii) the confidence intervals, $〈{\mathit{\xi }}^{\prime \prime }〉\pm {\mathit{\sigma }}_{{\mathit{\xi }}^{\prime \prime }}$ (black dashed lines), ${\mathit{\sigma }}_{{\mathit{\xi }}^{\prime \prime }}$ being the standard deviation of ξ^′′; (iii) ξ^′′Ens, evaluated on the basis of the ensemble-average concentration field (blue line). Inspection of Fig. 6a–c complements the qualitative analysis of Fig. 5 and suggests that heterogeneity causes (on average) a slight reduction of the toe penetration, $〈L_{\mathrm{T}}^{\prime \prime }〉$ and $L_{\mathrm{T}}^{\prime \prime \mathrm{Ens}}$ being slightly less than 1, and an enlargement of the mixing zone, $〈L_{\mathrm{S}}^{\prime \prime }〉$ and $〈W_{\mathrm{MZ}}^{\prime \prime }〉$ being larger than 1. The results of Fig. 6a–c also emphasize that, while $L_{\mathrm{T}}^{\prime \prime \mathrm{Ens}}$ calculated from the ensemble-average concentration distributions is virtually indistinguishable from $〈L_{\mathrm{T}}^{\prime \prime }〉$, $L_{\mathrm{S}}^{\prime \prime \mathrm{Ens}}$ and $W_{\mathrm{MZ}}^{\prime \prime \mathrm{Ens}}$ markedly overestimate $〈L_{\mathrm{S}}^{\prime \prime }〉$ and $〈W_{\mathrm{MZ}}^{\prime \prime }〉$, respectively, as they visibly lie outside of the corresponding confidence intervals of width ${\mathit{\sigma }}_{{\mathit{\xi }}^{\prime \prime }}$. Note that even as Fig. 6a–c have been computed along cross-section B–B^′, qualitatively similar results have been obtained along all vertical cross-sections, as also suggested by the behavior of the dimensionless areal extent of toe penetration, $A_{\mathrm{T}}^{\prime \prime }$ (Fig. 6d), and solute spreading, $A_{\mathrm{S}}^{\prime \prime }$ (Fig. 6e), as well as of the dimensionless volumes $V_{\mathrm{T}}^{\prime \prime }$ and $V_{\mathrm{S}}^{\prime \prime }$ (Fig. 6f–g). Overall, the results in Fig. 6 clearly indicate that the ensemble-average concentration field can provide accurate estimates of the average wedge penetration while rendering biased estimates of quantities characterizing mixing. Our findings are consistent with previous studies (e.g., Dentz and Carrera, 2005; Pool et al., 2015) in showing that an analysis relying on the ensemble concentration field tends to overestimate significantly the degree of mixing and spreading of the solute because it combines the uncertainty associated with sample-to-sample variations of (a) the solute center of mass and (b) the actual spreading. Our results also highlight that the extent of this overestimation decreases whenever integral rather than local quantities are considered: $V_{\mathrm{S}}^{\prime \prime \mathrm{Ens}}$ in Fig. 6g is indeed about half of $L_{\mathrm{S}}^{\prime \prime \mathrm{Ens}}$ and $A_{\mathrm{S}}^{\prime \prime \mathrm{Ens}}$.

3.2 Effects of pumping on SWI in three-dimensional heterogeneous media

Here, we investigate SWI phenomena in heterogeneous systems when a pumping scheme is activated as described in Sect. 2.3. Figure 7 collects contour maps of relative concentration C∕C_S obtained for all schemes along cross-section B–B^′ at the end of the 8-month pumping period, with reference to (i) Hom (the equivalent homogeneous system, left column) and (ii) Ens (the ensemble-average concentration field, right column). Contour lines C∕C_S = 0.25, 0.50 and 0.75 are also highlighted. Extracting FW according to the engineered solution S1 (Fig. 7c–d) results in a slight landward displacement of the SWI wedge and in an enlargement of the transition zone, as compared to S0 (Fig. 7a–b), where no pumping is activated. This behavior can be observed in the homogeneous as well as (on average) in the heterogeneous settings, and is related to the general decrease of the piezometric head within the inland side caused by pumping which, in turn, favors SWI. One can also note that the partially penetrating well induces non-horizontal flow (in the proximity of the well) thus enhancing mixing along the vertical direction. Wedge penetration and solute spreading are further enlarged in S2 (see Fig. 7e–f), where the total extracted volume is increased with respect to S1 through an additional pumping rate at the bottom of the aquifer. However, when the pumping well is operating within the transition zone (S3, Fig. 7g–h) the SW wedge tends to recede, approaching the pumping well location. In this case, the lower screen acts as a barrier limiting the extent of the SWI at the bottom of the aquifer.

Figure 7 Concentration distribution (color contour plots) and iso-concentration lines C∕C_S = 0.25, 0.5 and 0.75 (solid red curves) along the cross-section B–B^′ shown in Fig. 2a evaluated at the end of the pumping period within the equivalent homogeneous aquifer (left column) and the ensemble-average concentration field (right column). Black vertical lines represent the location of the well screens. Vertical exaggeration = 5.

Figure 8 Temporal evolution of ${\mathit{\xi }}^{\prime }=L_{\mathrm{T}}^{\prime }$, $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$ evaluated along the vertical cross-section B–B^′ shown in Fig. 2a during the pumping period, for all pumping schemes within the equivalent homogeneous domain (red lines) and the ensemble-average concentration field (blue lines). MC-based mean values of ξ^′ (black lines) and confidence intervals of width equal to ± 1 standard deviation of ξ^′ about their mean are also depicted.

The combined effects of groundwater withdrawal and stochastic heterogeneity on SWI are investigated quantitatively through the analysis of the temporal evolution of the seven metrics introduced in Sect. 2.4.

Figure 8 illustrates mean values of ξ^′ (with ${\mathit{\xi }}^{\prime }=L_{\mathrm{T}}^{\prime }$, $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$) computed for all pumping schemes considered across the collection of all heterogeneous realizations (black line) along the vertical cross-section B–B^′ versus dimensionless time $t^{\prime }=\left(t-t_w\right)Q/{\mathit{\lambda }}^{\mathrm{3}}$. Confidence intervals $〈{\mathit{\xi }}^{\prime }〉\pm {\mathit{\sigma }}_{{\mathit{\xi }}^{\prime }}$ are also depicted (dashed lines). As additional terms of comparison, Fig. 8 also includes ξ^′Hom (red line) and ξ^′Ens (blue line), respectively, evaluated in the equivalent homogeneous domain and in the ensemble-average concentration field. Figure 8 shows that the toe penetration, $L_{\mathrm{T}}^{\prime }$, and the solute spreading, as quantified by $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$, do not vary significantly with time in the absence of pumping (S0) because stationary boundary conditions are imposed for $t^{\prime }>\mathrm{0}$. Pumping schemes S1 and S2 cause the progressive inland displacement of the toe, together with an overall increase of spreading at the bottom of the aquifer. This phenomenon is more severe in S2, where SW and FW are simultaneously extracted. On the other hand, the vertical width of the mixing zone, e.g., $W_{\mathrm{MZ}}^{\prime }$, is not significantly affected by the pumping within schemes S1–S2. Configuration S3, in which the pumping well is located within the mixing zone, leads to the most pronounced changes in the shape and position of the SW wedge. The toe penetration first decreases rapidly in time and then stabilizes around the well location. Quantities $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$ show an early time increase, suggesting that a rapid displacement of the wedge enhances mixing and spreading. As the toe stabilizes over time, both $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$ decrease, reaching values equal to (or slightly less than) those detected in the absence of pumping. Consistent with the observations of Sect. 3.1, Fig. 8 supports the findings that (i) an equivalent homogeneous domain is typically characterized by the largest toe penetration and the smallest vertical width of the mixing zone, the only exception being observed in S3 at late times, when $L_{\mathrm{T}}^{\prime \mathrm{Hom}}\cong L_{\mathrm{T}}^{\prime \mathrm{Ens}}<〈L_{\mathrm{T}}^{\prime }〉$; (ii) basing the characterization of the mixing zone on ensemble-average concentrations enhances the actual effects of heterogeneity and yields overestimated mixing zone widths. Figure 8 also highlights that the extent of the confidence interval associated with toe penetration is approximately constant in time and does not depend significantly on the analyzed flow configuration. Otherwise, the confidence intervals associated with our mixing zone metrics depend on the pumping scenario and tend to increase with time if FW and SW are simultaneously extracted, especially in S2.

Figure 9 Isolines C∕C_S = 0.5 at the end of the pumping period along the bottom of the aquifer for the four schemes S0–S3 evaluated (a) within the equivalent homogeneous domain and (b) from the ensemble-average concentration field.

The fully three-dimensional nature of the analyzed problem is exemplified in Fig. 9, where we depict isolines C∕C_S = 0.5 along the bottom of the aquifer at the end of the pumping period for the equivalent homogeneous system (Fig. 9a) and as a result of ensemble averaging across the collection of heterogeneous fields (Fig. 9b). The spatial pattern of iso-concentration contours associated with pumping schemes S1–S3 departs from the corresponding results associated with S0 within a range extending for about 10λ around the well, along the direction parallel to the coastline.

Figure 10 Temporal evolution of ${\mathit{\xi }}^{\ast }=\left({\mathit{\xi }}^{\prime }-{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }\right)/{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }$, where ${\mathit{\xi }}^{\prime }=A_{\mathrm{T}}^{\prime }$, $A_{\mathrm{S}}^{\prime }$, evaluated for all pumping schemes within the equivalent homogeneous domain (red lines) and the ensemble-average concentration field (blue lines). MC-based mean values of ξ^∗ (black lines) and confidence intervals of width equal to ± 1 standard deviation of ξ^∗ about their mean are also depicted.

Figure 11 Temporal evolution of ${\mathit{\xi }}^{\ast }=\left({\mathit{\xi }}^{\prime }-{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }\right)/{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }$, where ${\mathit{\xi }}^{\prime }=V_{\mathrm{T}}^{\prime }$, $V_{\mathrm{S}}^{\prime }$, evaluated for all pumping schemes within the equivalent homogeneous domain (red lines) and the ensemble-average concentration field (blue lines). MC-based mean values of ξ^∗ (black lines) and confidence intervals of width equal to ± 1 standard deviation of ξ^∗ about their mean are also depicted.

We quantify the global effect of pumping on the three-dimensional SW wedge by evaluating, for each pumping scenario and each MC realization, the temporal evolution of ${\mathit{\xi }}^{\ast }=\left({\mathit{\xi }}^{\prime }-{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }\right)/\phantom{\left({\mathit{\xi }}^{\prime }-{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }\right){\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }}{\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }$, i.e., the relative percentage variation of areal extent, ${\mathit{\xi }}^{\prime }=A_{\mathrm{T}}^{\prime },A_{\mathrm{S}}^{\prime }$, and volumetric extent, ${\mathit{\xi }}^{\prime }=V_{\mathrm{T}}^{\prime },V_{\mathrm{S}}^{\prime }$, of wedge penetration and mixing zone with respect to their counterparts computed in the absence of pumping ${\mathit{\xi }}_{\mathrm{S}\mathrm{0}}^{\prime }$. Figure 10 shows that the penetration area, $A_{\mathrm{T}}^{\ast }$, and the solute spreading, as quantified through $A_{\mathrm{S}}^{\ast }$, at the bottom of the aquifer tend to increase with pumping time. Corresponding results for $V_{\mathrm{T}}^{\ast }$ and $V_{\mathrm{S}}^{\ast }$ are depicted in Fig. 11. The scenario where these quantities are most affected by pumping is S2, where the FW and SW pumping wells are located outside the transition zone (in place prior to pumping). Placing both wells within the transition zone prior to pumping yields a significant decrease of the effect of pumping on both areal and volumetric metrics. We further note that the mean penetration area and volume can be accurately determined by the solution obtained within the equivalent homogeneous domain and is associated with a relatively small uncertainty, as quantified by the confidence intervals. Heterogeneity effects are clearly visible on spreading along the bottom of the aquifer. Figure 10d–f highlight that $A_{\mathrm{S}}^{\ast \mathrm{Hom}}$ and $A_{\mathrm{S}}^{\ast \mathrm{Ens}}$ are not accurate approximations of $〈A_{\mathrm{S}}^{\ast }〉$. Spreading uncertainty, as quantified by ${\mathit{\sigma }}_{A_{\mathrm{S}}^{\ast }}$, is significantly larger than ${\mathit{\sigma }}_{A_{\mathrm{T}}^{\ast }}$, especially in pumping scenario S2. On the other hand, the spreading of the SWI volume, as quantified by $V_{\mathrm{S}}^{\ast }$ (Fig. 11d–f), can be accurately estimated by the homogeneous equivalent system, as well as by the ensemble mean concentration field, and it is characterized by a relatively low uncertainty.

Our analyses document that for a stochastically heterogeneous aquifer an operational scheme of the kind engineered in S3 (i) is particularly efficient for the reduction of SWI maximum penetration (localized at the bottom of the aquifer) and (ii) is advantageous in controlling the extent of the volume of the SW wedge, as compared to the double-negative barrier implemented in S2. However, this withdrawal system may lead to the salinization of the FW extracted from the upper screen due to upconing effects. This aspect is further analyzed in Fig. 12 where the temporal evolution of C_T∕C_S is depicted, C_T being the salt concentration associated with the total mass of fluid extracted by the upper screen (production well) in S3. Results for the equivalent homogeneous system (red curve), each of the heterogeneous MC realizations (grey dashed curves) and the ensuing ensemble-average concentration 〈C_T〉∕C_S (blue curve) are depicted. Vertical bars represent the 95 % confidence interval around the ensemble mean, evaluated as $〈C_{\mathrm{T}}〉/C_{\mathrm{S}}\pm \mathrm{1.96}{\mathit{\sigma }}_{C_{\mathrm{T}}/C_{\mathrm{S}}}/\sqrt{n}$. Here ${\mathit{\sigma }}_{C_{\mathrm{T}}/C_{\mathrm{S}}}$ is the standard deviation of C_T∕C_S and n=100 is the number of MC realizations forming the collection of samples. The width of these confidence intervals can serve as a metric to quantify the order of magnitude of the uncertainty associated with the estimated mean. The equivalent homogeneous system significantly underestimates the ratio 〈C_T〉∕C_S. Moreover, Fig. 12 highlights the marked variability of the results across the MC space. As such, it reinforces the conclusion that the mean value 〈C_T〉 is an intrinsically weak indicator of the actual salt concentration at the producing well.

Figure 12 Temporal evolution of C_T∕C_S for pumping scheme S3 within the equivalent homogeneous system (red curve) and for each heterogeneous realization (gray dashed curves). The ensemble-average (blue curve) dimensionless concentration and its 95 % confidence interval are also depicted.

4 Conclusions

We investigate quantitatively the role of stochastic heterogeneity and groundwater withdrawal on seawater intrusion (SWI) in a randomly heterogeneous coastal aquifer through a suite of numerical Monte Carlo (MC) simulations of transient variable-density flow and transport in a three-dimensional domain. Our work attempts to include the effects of random heterogeneity of permeability and groundwater withdrawal within a realistic and relevant scenario. For this purpose, the numerical model has been tailored to the general hydrogeological setting of a coastal aquifer, i.e., the Argentona river basin (Spain), a region which is plagued by SWI. To account for the inherent uncertainty associated with aquifer hydrogeological properties, we conceptualize our target system as a heterogeneous medium whose permeability is a random function of space. The SWI phenomenon is studied through the analysis of (a) the general pattern of iso-concentration curves and (b) a set of seven dimensionless metrics describing the toe penetration and the extent of the mixing or transition zone. We compare results obtained across a collection of n=100 MC realizations and for an equivalent homogeneous system. Our work leads to the following major conclusions.

Heterogeneity of the system affects the SW wedge along all directions both in the presence and absence of pumping. On average, our heterogeneous system is characterized by toe penetration and extent of the mixing zone that are, respectively, smaller and larger than their counterparts computed in the equivalent homogeneous system.

Ensemble (i.e., across the MC realizations) mean values of linear, $L_{\mathrm{T}}^{\prime }$, areal, $A_{\mathrm{T}}^{\prime }$, and volumetric, $V_{\mathrm{T}}^{\prime }$, metrics representing 50 % of SWI penetration virtually coincide with their counterparts evaluated in the ensemble-average concentration field.

Ensemble mean values of linear, $L_{\mathrm{S}}^{\prime }$ and $W_{\mathrm{MZ}}^{\prime }$, and areal, $A_{\mathrm{S}}^{\prime }$, metrics representing the mixing and spreading of SWI penetration are markedly overestimated by their counterparts evaluated on the basis of the ensemble-average concentration field. Therefore, average concentration fields, typically estimated through interpolation of available concentration data, cannot be employed to provide reliable estimates of solute spreading and mixing.

All of the tested pumping schemes lead to an increased SW wedge volume compared to the scenario where pumping is absent. The key aspect controlling the effects of groundwater withdrawal on SWI is the position of the wellbore with respect to the location of the saltwater wedge in place prior to pumping. The toe penetration decreases or increases depending on whether the well is initially (i.e., before pumping) located within or outside the seawater intruded region, respectively. The water withdrawal scheme that is most efficient for the reduction of the maximum inland penetration of the seawater toe is the one according to which freshwater and saltwater are, respectively, extracted from the top and the bottom of the same borehole, initially located within the SW wedge. This result suggests the potential effectiveness of the so-called “negative barriers” in limiting intrusion, even when considering the uncertainty effects stemming from our incomplete knowledge of permeability spatial distributions.

Salt concentration, C_T, of water pumped from the producing well is strongly affected by permeability heterogeneity. Our MC simulations document that C_T can vary by more than 2 orders of magnitude amongst individual realizations, even in the moderately heterogeneous aquifer considered in this study. As such, relying solely on values of C_T obtained through an effective homogeneous system or on the ensemble-average estimates 〈C_T〉 does not yield reliable quantification of the actual salt concentration at the producing well.

Future development of our work includes the analysis of the influence of the degree of heterogeneity and of the functional format of the covariance structure of the permeability field on SWI, also in the presence of multiple pumping wells.

Appendix A

Here, we assess the stability of the MC-based first- and second-order (statistical) moments associated with the metrics introduced in Sect. 2.4. Figures A1 and A2, respectively, depict the sample mean, 〈ξ^′′〉, and standard deviation, σ(ξ^′′), of all metrics evaluated at the end of the 8-year period versus the number of MC simulations, n. The estimated 95 % confidence intervals, computed according to Eqs. (3) and (8)

Figure A1 Sample mean of the seven dimensionless metrics at t=t_w and associated 95 % confidence intervals versus ensemble size, n.

Figure A2 Sample standard deviation of the seven dimensionless metrics at t=t_w and associated 95 % confidence intervals versus ensemble size, n.

of Ballio and Guadagnini (2004), are also depicted. Figures A1 and A2 show that the oscillations displayed by the quantities of interest are in general limited and do not hamper the strength of the main message of our work. As expected, moments of integral quantities (i.e., $A_{\mathrm{T}}^{\prime \prime }$, $A_{\mathrm{S}}^{\prime \prime }$, $V_{\mathrm{T}}^{\prime \prime }$, $V_{\mathrm{S}}^{\prime \prime }$, $W_{\mathrm{MZ}}^{\prime \prime }$) tend to stabilize faster than their counterparts evaluated for local quantities ($L_{\mathrm{T}}^{\prime \prime }$ and $L_{\mathrm{S}}^{\prime \prime }$).