Saltwater Intrusion in Coastal Aquifers: a Primary Case Study along the Adriatic Coast Investigated within a Probabilistic Framework

Environmentally sensitive areas along coastlines may be adversely affected by saltwater intrusion (SI), a condition which can be worsened by extensive groundwater extraction. Given the uncertainty of problem parameters, the risk of contamination of the vegetation capture zone needs to be cast in a probabilistic framework. In order to exemplify real situations existing along the Adriatic coast of Emilia-Romagna, a case study involving a pinewood strip and a well field drawing freshwater from an unconfined coastal aquifer was examined. On the basis of a widely adopted sharp interface formulation, key hydrogeological problem parameters were modeled as random variables, and a global sensitivity analysis was carried out to determine their influence on the position of the interface. This analysis utilized an efficient model reduction technique based on Polynomial Chaos Expansion. The risk that saltwater intrusion affects coastal vegetation was then evaluated via a two-step procedure by computing the probability that (i) the leading edge of the saltwater wedge reaches the sensitive area in the horizontal plane, and (ii) the freshwater/saltwater interface reaches the capture zone. The influence of the design parameters of the well field on the overall probability of contamination was investigated, revealing the primary role of the pumping discharge in the examined configuration.


Introduction
Contamination of freshwater bodies caused by saltwater intrusion (SI) is a global issue, affecting water quality, vegetation, and soil conditions along coastal lines.Deterioration of this freshwater resource threatens the sustainability of the water supply of coastal communities and their economic development [1,2].Freshwater stored in coastal aquifers is particularly vulnerable to degradation due to (i) its close proximity to seawater, and (ii) the significant water demand associated with coastal areas whereby groundwater is often the main source of drinking water [3].Continual degradation of groundwater quality in coastal aquifers has become an increasing concern in heavily urbanized regions across the world.In many countries, including Cyprus, Mexico, Oman, and Israel, hundreds of wells along the coastline have had to be shut down because of SI [4]; in California, saltwater from the Pacific Ocean has seeped into the Los Angeles Basin coastal aquifer, replacing freshwater [5].Extensive research is currently being carried out worldwide in order to understand the mechanisms of SI and develop strategies to control the process.For example, costly technological solutions, such as desalination, barriers or injection wells have been implemented to deal with the reduction in freshwater availability [6].
Over the last several decades SI has become a crucial issue in Italy as denoted by many studies conducted on different coastal areas [7,8].SI represents an even more serious threat along the coast of Emilia-Romagna (Figure 1), as it is combined with a strong anthropogenic subsidence, only partially limited in recent years.Though the subsidence in this area has a natural component, with an average of 0.5 mm/year [9], it is exacerbated by gas extraction [10], drainage to reclaim farmland, and groundwater over-exploitation, limited only since the late 1980s [11,12].The presence of vegetation and complex ecosystems in such areas (i.e., lagoons, marshes, and pine forests), also increases potential adverse impacts from an ecological point of view: understanding the mechanisms whereby these sensitive environments/ecosystems can be affected by SI represents a critical research challenge, as evidenced by recent literature findings [13].During the past twenty years, the coastal pine forest of Ravenna, along the Adriatic coast, has been affected by a progressive lowering in vegetation density, uprooting of tall plants, drying of pine plants, and changes in the richness of plant species.These impacts have been attributed to groundwater salinity, shallow water tables, drought, and air pollution [14,15].The trend in groundwater salinization underneath the coastal pine forest is a process that is progressing in time, as recent studies have demonstrated [16].
In recent years, numerical approaches have been increasingly adopted to represent the detailed and complex processes involved in SI [3].The characterization of the mixing zone as well as the evolution of intrusion temporally represent some applications in which the adoption of numerical models is relevant [17,18].Similarly, analytical formulations can provide a prompt insight into relevant physical processes involved and, for this reason, a significant number of studies have been conducted in this realm [19][20][21][22].In particular, the theory developed by Strack [23], based on the sharp interface approximation, has been widely adopted to predict the extension of SI in steady-state flow under the Dupuit-Forchheimer assumption [3,20,24,25].Several models describing situations of potential engineering interest have been developed and serve as a suitable reference to be re-examined from a risk analysis (RA) point of view [26,27].
A rigorous application of RA in environmental problems has been promoted by several environmental agencies and institutions [28,29].In the context of climate change, taking into account the uncertainty affecting hydrological parameters is even more important to make accurate predictions.In general, uncertainties are linked to (i) the incomplete knowledge of system dynamics (epistemic uncertainty), and (ii) the randomness which is inherent to natural phenomena (aleatory uncertainty) [30,31].This uncertainty, propagating towards the state variables of interest through a selected analytical or numerical model, has to be properly considered in a probabilistic framework.
In this work we expand upon this approach by providing a general methodology applicable to different contexts and model complexities.In particular, using examples and characteristics of the coastal aquifers found in the area of Ravenna in Emilia-Romagna, we propose a basic scheme representing the key features of steady-state SI in a probabilistic framework.First, an analytical model based on the theory of Strack [25] whereby groundwater withdrawal from a well within a phreatic aquifer is considered.Subsequently, the effects of well pumping on SI and coastal vegetation are investigated where the problem parameters exhibit associated uncertainty.This problem is addressed utilizing model reduction techniques [32,33], which provide an alternative to onerous and time consuming traditional simulation methods (e.g., Monte Carlo simulations).This methodology may be easily applied to complex case studies, assisting in the definition of proper management solutions for coastal aquifers.Regarding this point, in the exemplificative application presented here, we derive the dependence of the risk for the vegetation on the parameters characterizing the well field, quantifying the primary role of the pumping discharge.

Characterization of the Area of Interest
The Adriatic Sea is a shallow basin whose depth rises southwards and is characterized by a modest tidal excursion.The Adriatic coastline south of the Po River was originally part of its delta and alluvial plain.The province of Ravenna, bounded to the north by the Comacchio Valleys and to the south by the city of Cervia, is characterized by a high population, about 350,000 inhabitants whose density more than triples during the summer.In this area, many unpermitted and poorly-constructed wells have been drilled by bathing establishments and farms, resulting in environmental damage [35].Anthropogenic activities have altered the original ecosystem and now the once sandy coastline is characterized by river and canal mouths, lagoons, pine forests (e.g., San Vitale and Classe), heavily urbanized areas (e.g., the cities of Marina Romea and Cervia), industrial facilities, and reclaimed agricultural lands [36].The entire coast is threatened by SI in the phreatic aquifer and seawater encroachment inland along the rivers (e.g., Reno, Lamone, Fiumi Uniti, and Bevano rivers), natural and anthropogenic land subsidence, direct contamination from water bodies open to the sea, destruction of coastal dunes with consequent reduction of their barrier effect, insufficient aquifer recharge and sea level rise.In particular, over the last few decades the land surface has dropped to 4 m below mean sea level in some areas [36].This in turn has modified the rivers and regular groundwater flow regimes.These inland areas require protection from river flooding and marine intrusion with dikes and levees [36].A partially artificial drainage network was built in the coastal area to lower the phreatic level and keep the agricultural land dry.Unfortunately, this management process is complex and the presence of these modified drainage networks has not been able to mitigate SI and keep the pine forests healthy [36].
The coastal phreatic aquifer of Ravenna, the focus of this study, is primarily located within the littoral sands, and locally in the shallow marine wedge deposits.The groundwater basin is unconfined in the east, but 3-4 km from the coast it becomes overlain and confined by the most recent alluvial fine-grained continental deposits.The basin thickness varies from a maximum of 30 m in the area of San Vitale Forest to a minimum of about 15 m in the vicinity of Cervia [36].From east to west, the major sedimentary units consist of: (i) a wedge of fine-grained (fine sand to silty clays) sediments representative of a shallow marine environment; (ii) littoral sands representative of the foreshore, deep-shore, and the sand dune environments adjacent and inland of the beach.The westernmost portion of the study area is characterized by the presence of clay back-shore lagoon deposits [37,38].Some continental-type alluvial deposits (clay and silt) overlay the littoral sands and lagoon deposits in the west.Based on geologic data from slug tests and grain size analyses [36,39], the hydraulic conductivity of the aquifer is relatively high near the San Vitale Forest (1.1-24.8m/day) and Cervia (11-21 m/day), and is generally lower in the Classe Forest area (0.2-7.7 m/day).Along the entire coast, the hydraulic head measured in wells was almost always above mean sea level, with values ranging from 0.5 to 4 m in winter and from −0.5 to 3.5 m in summer [36].It was observed that water salinity (groundwater at the water table) varied seasonally.For instance, salinity values were higher (22-44 g/L) during the summer periods (e.g., near Comacchio Valleys) compared to winter months where values generally ranged between 3 and 22 g/L.Moreover, encroachment of sea water along rivers and canals, combined with subsidence and drainage from agricultural lands, contributes to enhance the salinization process [36].
Along the coast of Ravenna, pine forests are situated mostly on sandy dunes and were planted in the last century by farmers to protect their land from sea-salt spray.Plant ecologists have reported that the distribution of plant species in the coastal area varies depending on elevation, soil salinity, water salinity, and the duration of various stress periods by which the vegetation is affected [40].Based on soil and water salinity, plants along the coast can be divided into categories including very tolerant to salt (halophyte plants), tolerant/fairly, and tolerant/not tolerant to water/soil/sea-spray salinity [41].Increases in salinity within a coastal habitat will have direct and indirect effects on vegetation [15].In the study area of interest, new sand is added to the beaches every few years from dredging the offshore Adriatic Sea in order to counteract the continuous erosion along the coast.The pine forests (i.e., 50-150 m from the coast) and the adjacent farmlands are drained in winter to keep the groundwater level low.However, this increases the extent of SI in the phreatic aquifer, which is further enhanced by the transpiration of the pine trees [14].Natural pine seedling germination does not occur and new small pine trees are planted each year to replace dead plants.The most dramatic decrease in the diversity of plant species occurs where water salinity exceeds the 3 g/L threshold, above which most trees species become stressed, with the exception of a few, i.e., Pinus pinaster (William Aiton, 1789) and Pinus pinea (Carl von Linné, 1753) [15].

Theoretical Framework
In this study, the theory developed by Strack [23] was used to predict the extension of SI into a coastal phreatic aquifer for the case of steady state flow under the Dupuit-Forchheimer assumption [20,22,24,25].The sharp interface approximation and the homogeneity assumption allow for the derivation of a closed form solution that describes the main processes involved in SI.Consider a semi-infinite domain bounded by a straight coastline (Figure 2a), and divided into two zones, with freshwater in zone 1 and freshwater and saltwater in zone 2. The aquifer thickness b [L] is: where d [L] is the mean sea level (MSL) with respect to the bottom of the aquifer, ξ [L] the interface depth with respect to the MSL, h f [L] the hydraulic head and q [L 2 T −1 ] the uniform discharge to the sea per unit width of aquifer [42].In zone 2, the Ghyben-Herzberg assumption ξ is imposed to express the relationship between hydraulic head and interface depth (with respect to MSL).
Following the approach of Strack for an unconfined aquifer [23], the potential φ [L 2 ] is defined as: these functions and their first derivatives are continuous across the multiple-zone aquifer and satisfy the Laplace equation, , in the x-z plane.At the interface, substituting in Equations ( 3) or (4) the condition ξ = d, hence h f = (γ s /γ f ) ξ, the potential is obtained as: Now consider then the presence of a pumping/recharging well with discharge Q w [L 3 T −1 ], located at a distance x w from the coast (Figure 2b,c).The solution of the problem, in terms of potential, can be found by resorting to the superposition principle and the method of images as [23]: where K [LT −1 ] represents the hydraulic conductivity of the aquifer.The location of the leading edge, or "toe", of the saltwater wedge can be computed from Equations ( 5) and ( 6) as: Finally, combining expressions ( 3) and ( 4) with ( 6) allows for the determination of the elevation of the hydraulic head above the MSL: ) utilizing the Ghyben-Herzberg assumption, the interface depth is ξ = (γ f /(γ s − γ f ))(h f -d) for zone 2, and ξ = d for zone 1, whereby no saltwater is present.For multiple pumping wells, Equation ( 6) can be extended to: where (x i , y i ) are the coordinates and Q i is the pumping rate of well i, which can be taken as positive or negative to represent pumping and recharging mode, respectively.With these assumptions, Equations ( 8) and ( 9) become for zone 1 and 2, respectively: Figure 2 represents the effects on the freshwater-saltwater interface due to the presence of one or more wells.While under natural undisturbed conditions the interface maintains a state of equilibrium (Figure 2a), water extraction lowers the water table with respect to its initial position (Figure 2b, blue solid line) allowing the interface to move inland.As the pumping rate increases, the interface advances towards the well until an unstable critical situation is reached (Figure 2c, red solid line).Any further increase in the pumping rate causes a jump of the toe position landward of the well, which thereby becomes contaminated by saltwater.

Problem Formulation
In order to evaluate the effects of SI on coastal vegetation due to both natural and anthropogenic causes, a simplified conceptualization representing specific areas of the Emilia-Romagna coast was configured (Figure 3).An 'environmentally sensitive area' representing a pinewood forest (an infinitely long strip) was located at a minimum distance x* = 50 m from the coastline.The lower limit of the plants capture zone (the volume of soil affected by water uptake by roots) is defined such that the sensitive zone surrounding the plant is not impacted by the intruding saltwater and any adverse effects on the plants are prevented.The extension of this zone depends on both the length of roots and the width of the transpiration zone; its minimum elevation from the bottom of the aquifer is set to z lim = 12 m in the present case.A well field established for agricultural purposes was located inland of the pinewood, as depicted in Figure 3.The well field was comprised of three wells, each located in the middle of an area of dimensions b•e, distant a from the coast; between every pumping area there are non-overlapping zones of width equal to c, to reduce mutual interactions among the wells.The total discharge Q was equally divided among the wells.Values assumed for the design parameters of the well field are summarized in Table 1.
Following a probabilistic approach, we model the hydrogeological parameters K, q and d, as independent random variables.Based on previous studies [42], all parameters were assumed to follow a normal distribution, characterized by the average values and the coefficients of variation reported in Table 2.Note that the coefficient of variation of d was selected according to the prediction of annual sea level rise [12] and subsidence for the area of interest.
Once the uncertain parameters are defined, the extension x toe of SI in the x-y plane is evaluated through Equations ( 5) and (10).The computed values of x toe are affected by the uncertainty propagated from the input parameters through the selected model equations.Since the domain considered is symmetric, the maximum value of x toe occurs along the middle, at y = y crit .The extent of intrusion can be determined in the x-z plane corresponding to y crit , by computing the elevation of the interface z i with respect to the bottom of the aquifer via Equation (12).This is also a random variable.The risk for the vegetation, here considered as the probability that the saltwater interface reaches the sensitive zone below the plants roots, is then derived following two steps.First, the probability that plant roots occur in the zone impacted by SI (i.e., zone 2 in Figure 4) is computed as: is the cumulative distribution function of the random variables x toe .According to the data results P 1 ≈ 1, i.e., the vegetation falls into zone 2 almost certainly within the selected case study.Secondly, the probability that the saltwater interface reaches the sensitive zone below the plants roots is evaluated in correspondence of the maximum intrusion, i.e., y = y crit .This probability is conditional on the previous one as follows: ( ) ( ) (14) In this selected case study, as P 1 ≈ 1, Equation ( 14) reduces to , where i z F is the cumulative distribution function of the random variable z i .Consequently, in the following section the simplified formulation for the computation of the risk, i.e., R = P 2 , was considered.
The proposed methodology is suitable to predict the risk to which coastal vegetation is subjected, considering the uncertainty of the parameters involved.For this purpose, the evaluation of the mean position of the toe (considering the schematic scenario described above) and its variance, applying global sensitivity analysis [32], must be conducted in order to evaluate the influence (sensitivity) of each uncertain parameter.Subsequently, the risk affecting the vegetation is computed for different scenarios defined through the parameters that characterize the well field and modifying their values with respect to those reported in Table 1.In particular, variations in (i) the overall pumping discharge, Q; (ii) the distance of wells from the coast, a; and (iii) the width of non-overlapping zones, c, were considered.Table 1.Design parameters of the well field represented in Figure 3.

Model Reduction Approach for Risk Analysis
In order to perform the stochastic analysis described above, time intensive simulation methods based on the Monte Carlo (MC) approach are typically applied.For this analysis, a different, more efficient strategy was implemented that relied on model reduction via Polynomial Chaos Expansion (PCE) [43].The basic idea of model reduction methods consists in modifying the form of the selected model without introducing any simplifying assumptions in the description of the physical processes involved.For the PCE approach, a surrogate model was derived in a simple polynomial form depending only on the uncertain input parameters.The computation of the PCE approximation requires a very limited number of MC simulations performed on the original model.
Let the model response of interest, y, belong to the Hilbert space of second-order random variables.Under this assumption, y can be approximated through the PCE technique [43,44] as follows: The resulting formulation constitutes a metamodel, y ~, of y.This metamodel is a simple polynomial function which is expressed in terms of a set of independent random variables, collected in vector ζ .
Here, ( ) ( ) , is the number of terms of the expansion, and q is the maximum degree considered; a j represent unknown deterministic coefficients, while Ψ j denote the suitable multivariate polynomial basis in the Hilbert space containing the response.
A non-intrusive regression-based approach can be employed to calculate the coefficients a j [32].As mentioned above, this method requires a number N of simulations performed on the original model proportional to P.
Once the PCE is available, simple analytical post-processing allows derivation of (i) the first two moments of the approximated response, and (ii) the global sensitivity indices of Sobol, that represent dimensionless ratios among the partial variance of the response, due to a generic subset of uncertain model parameters, and the total variance of the response [32].Specifically the mean of the model response coincides with the coefficient of the zero-order term, a 0 in the Equation (15), and the total variance of the response and the generic Sobol index are given, respectively, by: for the computation of [ ] 2 α Ψ E see e.g., [45].
Furthermore, the PCE metamodel represents a suitable basis to perform Monte Carlo simulations, allowing for a more efficient means to perform risk analysis or to solve inverse problems with reduced computational time/cost.This methodology has demonstrated excellent promise for use in many other applications related to flow and transport processes in porous media (e.g., [46,47]).

Results and Discussion
Applying the methodology outlined in Section 3.2, the extension of SI was first evaluated by using the values (parameters) collected in Table 1 to represent the well field.The uncertainty in the hydrogeological parameters was modeled using probability distributions as reported in Table 2.Under this framework, the first two moments associated with x toe and the global sensitivity indices were evaluated as described in Section 3.3.
The mean value of x toe as a function of y is depicted in Figure 5a, showing the intervals of width ±σ and ±2σ surrounding it.The maximum value for the mean of x toe (i.e., the leading edge of the intrusion) is reached at y crit = 3b/2 + c, since the domain is symmetric.The maximum intrusion extent was equal to 303.6 m for the selected case study; this value is consistent with those typically observed for the area of Ravenna, and with results reported in numerical studies of the effects of subsidence on SI [35].The uncertainty in the hydrogeological parameters resulted in a standard deviation of x toe with a maximum of ~50 m (i.e., more than the 15% of the correspondent mean value).
The standard deviation σ of x toe is illustrated in Figure 5b as a function of y, showing the individual contributions due to each random parameter value ( K , q and d) in order to investigate how these different values influence SI.The results indicate that the primary influence on SI was due to the variability in q, whereby, in contrast, the mean sea level, d, was relatively inconsequential.However, the impact of d was shown to produce significant effects on SI in a long-term period of time due to climatic changes and the combined increase in natural/anthropogenic subsidence.As expected, the effect of domain heterogeneity, exemplified through the variability of K, was also significant.These results highlight the importance of proper site characterization in order to obtain accurate predictions for SI.Other studies (e.g., [22]) have relied on the theory presented by Strack [23], whereby a sensitivity-based approach was adopted in order to define vulnerability indicators for SI.However, these global sensitivity analysis models do not account for the simultaneous variability of all the parameters of interest when evaluating effects on SI.This study also evaluated he influence on SI on the risk for vegetation health due to variations of the deterministic parameters Q (overall pumping discharge), a (distance of the well field from the coast), and c (width of non-overlapping zones).
The mean and the standard deviation of x toe versus y for different values of total pumping rate Q are depicted in Figures 6a and 6b respectively.In particular, it is shown that the mean of x toe reaches its maximum at y = y crit , ranging between 280 and 365 m.This revealed a relatively high degree of variability for different values of Q i ; specifically, the largest extent of SI was reached for the maximum pumping rate considered (i.e., Q i = Q/3 = 300 m 3 /day) (Figure 6a).The standard deviation σ of x toe followed the same trend (Figure 6b), whereby values ranged between 45.7 and 74.8 m, primarily as a function of the variability of domain properties, q and K, as described above.
A similar analysis is illustrated in Figure 7a,b, with the exception that different values of a were considered to evaluate effects on SI.As previously noted, the results of this study revealed that the distance of the well field from the coast, within the range considered, seemed to be inconsequential, as the two moments of x toe did not exhibit significant changes in the extent of SI.In particular, the maximum value of the mean of x toe ranged between 294 and 308 m (i.e., trend increasing with a) (Figure 7a), and the associated standard deviation varied between 49.7 and 54 m (Figure 7b).
Finally, the influence of c (the width of the non-overlapping zones among the wells) on the first two moments of x toe is depicted in Figure 8a,b.By varying the width of the non-overlapping zones, the size of the entire domain is modified.For example, maximum values of mean x toe and standard deviations were obtained for smaller values of c , whereby y crit = 3b/2 + c tended to decrease.For c = 50 m, y crit = 500 m, with a mean value of x toe equal to 307 m and an associated standard deviation of 53.6 m; for a value of c = 200 m, y crit becomes 660 m, whereby the mean value of x toe was reduced to 297 m with a standard deviation of 50.2 m.
For this analysis the design parameters of the well field demonstrated that SI was extremely sensitive to the pumping discharge Q, in contrast to the other parameters considered (i.e., a and c) which seem to have only a minimal influence on the optimization of groundwater extraction and effect on SI.  domain heterogeneity (i.e., permeability distribution) significantly influenced the extent of SI, consistent with findings reported by other studies [20,48].

Conclusions
This study presents a computationally efficient technique to evaluate risk to ecosystems (i.e., vegetation health) along coastal regions that are impacted by saltwater intrusion processes.The methods used were applied to the coastline of Emilia-Romagna, Italy, a region that has been and continues to be impacted by saltwater intrusion.A stochastic approach to the problem was adopted to represent the random nature of the hydrogeological parameters.Computational efficiency was further increased by incorporating model reduction techniques via Polynomial Chaos Expansion theory for modeling saltwater intrusion effects and associated risk analyses.
A simplified analytical formulation was adopted to allow for the evaluation of the probability that the vegetation capture zone overlaps with the intruding saltwater wedge.Anthropogenic influence on the risk to vegetation due to saltwater intrusion was modeled through the presence of a well field, specifically accounted for through adjustment of (i) the pumping discharge; (ii) the wells distance from the coast; and (iii) the mutual distance between wells.For the simplified scenario representative of the Ravenna coast, the comprehensive results of the modeling analyses revealed the parameters most responsible for influencing the extent of saltwater intrusion.For example, increases in the overall pumping discharge of the well field produced the most significant increase in the probability that saltwater intrusion reaches the vegetation capture zone.
The computational methodology used herein can be applied to more complex systems whereby a closed-form solution is not available and would require more time-consuming numerical approaches (i.e., to resolve the shape and size of the intruding saltwater wedge).In such cases, the model reduction technique adopted for these analyses allowed for both excellent representation of the physical system and an increase of computational efficiency.Moreover, the current approach can assist in the development of appropriate management strategies for coastal aquifers.These techniques allow the parameters describing the anthropogenic influences to be considered as random variables, which provide an improved means of discriminating the parameters most responsible for saltwater intrusion and for solving optimization problems that require a specific set of environmental constraints.

Figure 2 .
Figure 2. (a) Reference unconfined aquifer with saltwater wedge; (b) Effect of a pumping well on saltwater wedge; (c) SI fronts near the pumping well (blue and red lines), in the x-y plane, for two different values of pumping rate.

Figure 3 .
Figure 3. Well field and location of pine forest in x-y plane.

Figure 4 .
Figure 4. Simplified conceptualization of the problem, in x-z plane, for y = y crit .

Figure 6 .F
Figure 6.(a) Mean value of x toe ; (b) Standard deviation σ of x toe ; (c) Cumulative distribution function i z F of the surface elevation i z computed in lim z (Figure 4); and (d) Computed risk [see Equation (14)] resulting from different values of single pumping discharge Q i .

Figure 7 .
Figure 7. (a-d) As Figure 6 but for different values of distance from the coast a.