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:

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:

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

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

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 1Parameters adopted in the numerical model.

Download Print Version | Download XLSX

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

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.

Download

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

Download

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 4Permeability 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.

Download

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 5MC-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.

Download

Figure 6Metrics ${\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^′.

Download

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

Download

Figure 8Temporal 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.

Download

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 9Isolines 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.

Download

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 10Temporal 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.

Download

Figure 11Temporal 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.

Download

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 12Temporal 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.

Download