Next Article in Journal
Streambed Flux Measurement Informed by Distributed Temperature Sensing Leads to a Significantly Different Characterization of Groundwater Discharge
Next Article in Special Issue
Mapping of Groundwater Spring Potential in Karst Aquifer System Using Novel Ensemble Bivariate and Multivariate Models
Previous Article in Journal
Analysis of 15N-NO3 Via Anoxic Slurries Coupled to MIMS Analysis: An Application to Estimate Nitrification by Burrowing Macrofauna
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of the Saltwater Wedge in a Coastal Karst Aquifer with a Double Conduit Network, Numerical Simulations and Sensitivity Analysis

1
Department of Chemistry, Life Sciences and Environmental Sustainability, Parma University, Parco Area delle Scienze, 157/A, 43124 Parma, Italy
2
Department of Engineering and Architecture, Parma University, Parco Area delle Scienze, 181/A, 43124 Parma, Italy
3
Department of Earth Sciences, Pinar del Rio University “Hermanos Saiz Montes de Oca”, A4, Pinar del Rio 20100, Cuba
*
Author to whom correspondence should be addressed.
Water 2019, 11(11), 2311; https://doi.org/10.3390/w11112311
Submission received: 26 August 2019 / Revised: 28 October 2019 / Accepted: 30 October 2019 / Published: 5 November 2019
(This article belongs to the Special Issue Groundwater Modelling in Karst Areas)

Abstract

:
We investigate the long-distance salinity in a dual permeability coastal karst aquifer with a double conduit network using a three-dimensional variable-density groundwater flow and multispecies transport SEAWAT model. Sensitivity analyses were used to evaluate the impact of the parameters and boundary conditions on the modeling saltwater wedge in a karstic aquifer situated in the Cuban land territory, including hydraulic conductivity, vertical anisotropy and salinity concentration; both in the conduits network and the fractured medium. These analyses indicated that hydraulic conductivity of the fractured medium and salt concentration were the ones that have a stronger effect on saltwater intrusion in a karstic aquifer. We also show results of the three-dimensional numerical simulations on groundwater salinity for different scenarios with the variabilities of the important parameters and compare results with electric conductivity profiles measured in a well.

1. Introduction

Coastal karst aquifers have an important role as water resources. These aquifers have hydraulic links with the sea resulting in dominant or important conduits flow conditions, submarine freshwater springs and/or natural seawater intrusion into the aquifer through karst conduits. The management of saltwater intrusion into a coastal aquifer is one of the most challenging problems due to its very complex network of conduits. The contamination of fresh groundwater by saline seawater may be intensified by different natural or anthropogenic factors. The main factors are related to the climatic regime, the variation of sea-levels rise, over-pumping for water supply and agricultural activities. The flow in the conduits open to the sea depends on the hydraulic head gradient between the aquifer and the sea and it is, therefore, a function of the water density and head losses in the aquifer. Coastal groundwater resources are vulnerable to frequent seawater-freshwater exchanges in a karst aquifer, which consists of high permeability conduits and low permeability fractured medium with sinkholes and karst windows that are usually connected by well-developed subsurface conduit networks [1,2,3,4,5].
Modeling saltwater flow in freshwater systems requires numerical codes that can solve a coupled variable-density flow and transport equation in the transition zone between freshwater and saltwater. Several variable-density numerical methods have been developed and used to study seawater intrusion, including SUTRA [6] and FEFLOW [7]. SEAWAT is a widely used variable-density numerical code that solves the coupled flow and transport equations, using a finite difference method with an implicit procedure. SEAWAT is a coupled version of MODFLOW [8] and MT3DMS [9,10,11,12]. In this contest, since groundwater flow in a karst conduit is often non-laminar [13,14,15] several numerical codes, such as MODFLOW-CFPM1 [14] and CFPv2 [16,17,18] have been developed to simultaneously solve Darcy’s flow in the fractured medium, the non-laminar flow inside the karst conduits and the exchanges between both systems. However, these constant-density karst models have limitations in simulating the density-dependent seawater intrusion processes in a coastal aquifer [19]. The VDFST-CFP, developed in Reference [20], is based on a density-dependent approach to study seawater intrusion in a coastal karst aquifer with conduits. However, have some computational constraints associated with the aquifer geometry. Therefore, the variable-density SEAWAT model can still be applied, where Darcy’s equation is used to compute flow not only in the fractured medium but also in the conduit with large values of the hydraulic conductivity and effective porosity [19].
In this paper, we are interested to investigate the capabilities of the three-dimensional SEAWAT code, as in Reference [19], using a numerical model of the aquifer, based on a simple conceptualization of a study case, recently published [21]. Since fewer observational data are possible, a sensitivity analysis is important to provide insight on which are the parameters that might result in a variation of the observed quantities and which are not important [22,23,24].
Fewer studies have addressed the problem of the parameter’s sensitivity for seawater intrusion numerical simulations in a coastal karst aquifer. Shoemaker shows a sensitivity analysis of the SEAWAT code for seawater intrusion in a homogeneous porous aquifer and concluded that dispersivity is the most important parameter in the head, salinity and groundwater flow numerical simulations and observations in the wedge zone [25]. Also, Xu et al. addressed this issue using a two-dimensional model and concluded that salinity and head simulations in the karst features, such as the conduit system are critical for understanding seawater intrusion in a coastal karst aquifer [26]. They also evaluated the hydraulic conductivity sensitivity and they found that it may be biased since the conduit flow velocity is not accurately determined by Darcy’s equation as a function of the hydraulic gradient. In Reference [26], they use an improved variable-density flow and solute transport-conduit flow process model inside the conduits where Darcy’s law is not anymore satisfied.
In this paper, we performed three-dimensional numerical simulations of seawater migration into a karst aquifer with two parallels conduits. This is a simplified model of a real case corresponding to a coastal carbonate aquifer with karst conduits, recently published [21], and investigate the saltwater intrusion as a function of different parameters and boundary conditions. We performed a local sensitivity analysis and then we investigated the correlation between the different parameters and show its nonlinear dependence. We performed three-dimensional simulations of saltwater concentration using two karst conduits with high permeability and porosity and a fractured medium with a relatively low hydraulic conductivity and porosity and investigate the salt “contamination” inside the aquifer. Salt intrusion in the presence of karst conduits is a three-dimensional problem, especially in the zone close to the conduits. A similar study has been performed in Reference [19] using a two-dimensional SEAWAT model of seawater intrusion in a dual-permeability coastal karst aquifer with a single conduit and in [26] with a more complex coastal karst aquifer. Here we extend this study using a three-dimensional parallelepiped grid that contains two parallel conduits along with the entire site. The objective of the numerical model is to describe the variation of the salinity concentration observed in a monitoring well as a function of the depth. The model is based on the parameter values of a fractured medium aquifer with two karst conduits with large values of both, the hydraulic conductivity and the effective porosity.
The organization of the paper is the following. After a brief Introduction, in Section 2, we introduced the numerical setup, hydrological conditions, model discretization, and boundary conditions. We also review the method used for the sensitivity analysis. In Section 3, we reported the results of the sensitivity analysis using four parameters. The scenarios of the saltwater wedge numerical simulations with different salinity concentration and comparison with experimental data. Different boundary conditions and aquifer’s depth are also investigated. The summary and conclusions are given in Section 4.

2. Methods

To study the saltwater intrusion, density-dependent groundwater flow dynamics are needed to simulate flow in the transition zone between freshwater and saltwater. In this paper, we used a local sensitivity analysis to evaluate the derivatives of the model simulations with respect to the parameters at specified values [22,23,24]. The forward difference approximation of sensitivity is calculated as the derivative of the i-th simulation with respect to the j-th model parameters:
y i x j y i ( x + Δ x ) y i ( x ) Δ x j
where y is the value of the i-th simulation; x is the j-th estimated parameter, x is a vector of zeros except that the j-th parameter equals Δ x j and y i / x j is the derivative or sensitivity of the simulated value [22]. The sensitivities indicate the slope of a plot of a simulated value y i relative to one parameter or, approximately, how much a simulated value would change if a parameter value were changed, divided by the change in the parameter value. The parameters are considered individually. Sensitivities do not account for changes in multiple parameters [22,23,24] and can be used to indicate the importance of the observation variable (in this case, salt concentration and/or head) to the estimation of the parameter values.
Since parameters can have different units, a scaling method is used to calculate the dimensionless scaled sensitivities d s s values that satisfy the following equation:
d s s i j = ( y i x j ) | x | x j | ω i i 1 / 2
where d s s i j is the dimensionless scaled sensitivity of the i-th simulation with respect to the j-th parameter, and ω i i is the weight of the i-th simulation, given by
ω i i = 1 σ 2
where σ is the error standard-deviation. As in the Reference [25], the measurement error is normally distributed with a standard deviation of the order of 0.003 m for the head, while for the salinity measurement, the error is of the order of 0.1 PSU (Practical Salinity Unit). They are based on standard error estimates for water levels measured in a well.
The d s s values of different simulations with respect to each parameter are accumulated as the composite scaled sensitivity ( c s s ) values and provide the total amount of information provided by the simulation for the estimation of one parameter. The c s s of the j-th parameter is given by
c s s i j = ( i = 1 N D ( d s s i j ) 2 | x / N D ) 1 / 2
where N D is the number of simulated quantities, in terms of the head and salinity simulations.
Composite scaled sensitivities are used to determine the relative importance of various flow and transport parameters to reproduce observed values and as a measure of the amount of information provided by the set of observations for estimating a parameter value. Larger c s s values identify parameters for which more information is provided by all observations for each parameter. Correlation coefficients, c o r ( i , j ) , are defined as:
c o r ( i , j ) = c o v ( i , j ) v a r ( i ) 1 / 2 v a r ( j ) 1 / 2
where c o v ( i , j ) is the covariance between parameter i and j ; v a r ( j ) is the variance of the parameter j . Correlation coefficients are used to identify parameters that are extremely correlated given the observations used in the numerical simulations. Parameters correlation coefficients with values of +1.00 or −1.00 indicate parameter values that are extremely correlated and generally cannot be estimated uniquely with the observations involved; values <     0.95 indicate that unique estimates can likely be obtained [25]. In this paper, we used the local sensitivity analysis to evaluate the parameter sensitivities at one specified value for each parameter instead of a range. Moreover, the local sensitivity indices are based on the first-order derivative where it is assumed a linear relationship of simulated quantities with respect to the parameters.

2.1. Study Site

The numerical model developed in this paper is based on the conceptual model and the parameter values of a fractured and karstified carbonate aquifer of the Cuban land western territory described in a recent work of Hernàndez-Diaz et al. [21]. Carbonate aquifers are the primary source of freshwater but are often contaminated by seawater intrusion that might be due to rain regime, sea tides that may affect the hydraulic gradient (which in principle prevents the intrusion), or excessive aquifer exploitation by wells for agriculture or human consumption [27,28,29,30]. The conceptual model of the study area is simplified in the sketch of Figure 1, which represents a schematic cross-section picture of a coastal karst aquifer with two conduits networks opening to the sea. The direction of the flow is toward the sea with a hydraulic gradient of 0.001 m. At about 5 km from the sea there is a monitoring well (well 18) in which are available several electrical conductivity profiles.
A three-dimensional SEAWAT model is set up to simulate seawater intrusion and investigate the three-dimensional flow and transport in the fractured medium. We are also able to investigate all parameters variation in several locations of the aquifer and highlight the differences between measured calculated both, inside the two conduits and determine whether the depth of the conduits affects the results. Furthermore, we investigated the importance of the fractured medium, both near or away from the conduits. Thanks to the three-dimensional numerical implementation of this study, we are able to give a numerical result at any point in the grid.
To simplify the three-dimensional numerical simulations, we assume that the aquifer’s grid geometry is a parallelepiped of 10 km long, starting from the seashore where the constant head is set to 0.0 m, up to a constant head equals to 10.0 m in the opposite side of the parallelepiped, situated at the inland direction; 1.0 km wide and −45.0 m depth. Then two parallel conduits situated along the entire parallelepiped at a depth of −14 m and −24 m, respectively, in concomitance with the observed step-like shape of the Electrical Conductivity (EC) profile (Figure 2). Each of them corresponds to a small parallelepiped of area 1   m 2 and a different hydraulic conductivity with respect to the fractured medium.

2.2. Hydrological Parameter

In Table 1, we show the range of hydrological parameter values used in the numerical simulations and sensitivity analysis. These parameters were calibrated in the regional-scale groundwater flow steady-state and solute transport transient models by [21]. In any case, a calibration analysis is not performed in this paper since the head and salinity observational field data are not enough, particularly inside the conduits. Some of the parameter values of Table 1 (hydraulic conductivity both, in the conduits and in the fractured medium, salinity concentration, and vertical anisotropy) are evaluated in the local sensitivity analysis and then applied in the different scenarios later on. The values of the hydrological parameters, hydraulic conductivity, effective porosity, longitudinal dispersivity (we do not investigate specific storage and rainfall recharge) in the conduit are greater than those corresponding to the fractured medium. For simplicity, we chose two values of the hydraulic conductivity: one inside both conduits and another value for the fractured medium, hydraulic conductivity in the fractured medium (HCf) and the hydraulic conductivity inside the conduits (HCc), respectively. See Table 1 for details.
The hydraulic conductivity of the fractured medium is assigned to be 1045   m / day = 1.21 × 10 2   m / s , considering a mean velocity of 10 m/day for a typically fractured media (see Reference [31]). This is due to the fact that numerous small fractures and relatively large pores existing in the karst aquifer associated with the dissolution of carbonate rocks. For the conduits system instead, we assigned a very high value of the hydraulic conductivity as 2.4 × 10 6   m / day = 27.8   m / s based on the value of 2371 m/d as a mean velocity for a typical karst conduit.
Some of the parameters of Table 1, i.e., the porosity for both the fractured and the conduits and the longitudinal dispersivity for both, the fractured and the conduits, were not investigated in the local sensitivity analyses, but they were fixed to specific values. The effective porosity was not investigated in the local sensitivity analysis and set to its maximum value of 1.0 in both conduits, while it is set to 0.1 in the surrounding fractured medium. The longitudinal dispersivity was estimated to be 10.0 m in the fractured medium and very small (0.3 m) along both conduits since advection is dominating compare to dispersion, which is negligible in the transport equation inside the conduit (in agreement with previous simulations in karst aquifers [19]).
We also investigated the value of the vertical anisotropy defined as the ratio between the horizontal hydraulic conductivity, Kh, and the vertical one, Kv, (VA = Kh/Kv) and was set to 1.0 in the local sensitivity analysis. Constant-head and constant-salinity boundaries conditions were used to represent the ocean body. The constant-head and constant-salinity boundaries were assigned values equal to the sea level (0.0 m) and the salinity of the seawater (37 kg/m3), respectively.

3. Results and Discussions

A numerical simulation analysis of the saltwater wedge in a dual permeability coastal karst aquifer with two conduit networks has been performed using a three-dimensional variable-density groundwater flow and multi-species transport SEAWAT model. The model successfully describes the variation of the Electrical Conductivity (EC) (as indirect measure of the salinity concentration) observed in “well 18”, as a function of the depth (see Figure 2) superimposed with a different x-scale axis. It is based on the parameter values of a fractured aquifer with two karst conduits with large values of both, the hydraulic conductivity and the effective porosity.

3.1. Spatial and Temporal Discretization

The grid discretization and boundary conditions of the three-dimensional SEAWAT numerical model are shown in Figure 3. We used the finite-difference grid frame of the three-dimensional model that covers a total area of 10 km2 and 45 m depth is set to a rectangular parallelepiped (see Figure 3) of 10 km long (with 400 column), 1 km width (with 33 rows) and a thickness of 26 layers in the cross-section, for a total number of 343,200 cells. In general, a fine resolution vertical grid is required for accurately modeling the density-dependent flow and solute transport.
The horizontal discretization for each cell is set uniformly to 2 m except for columns 8, 9, 19, 11 and 16, 17 where a fine resolution of 1 m each is implemented in correspondence with both conduits. A first conduit is situated at layer 10 (− 14   m a.s.l.) and row 16, where the first step is observed in Figure 2 (and also a refinement is implemented). A second conduit is located below the first one and parallel to it at layer 16 (− 24   m a.s.l.) and row 16, where the second step is situated in Figure 2. Both conduits across the whole horizontal area of 10 km long. For simplicity, the size of the horizontal conduits is assumed to be constant all over the parallelepiped. The outlet of both karst conduits of 1   m 2   is in contact with the sea boundary on the left side where the wedge will take place.
See Figure 3, for details on the grid used for the numerical simulations: A three-dimensional view (top), a cross-section view where the conduits are situated (bottom left) and front view of the grid (bottom right). One may observe both conduits (in red) along the whole grid up to 10   km and a grid refinement near the conduits both, in the y and z directions. While in the bottom there is a three-dimensional view of the grid represented by blue lines and both conduits corresponding to the red lines at −14 m (layer n. 10) and −24 m (layer n. 16), respectively. Notice the refinement in the proximity of both conduits. The head is set constant to 10.0 m in the inland region ( x   =   10 4   m ) and constant head equal to 0.0   m in contact with the saltwater. For this reason, the first layer is not fully saturated. Figure 3 was generated using the MODFLOW’s output (in hdf5 format) on FloPy [32].

3.2. Initial and Boundary Conditions

Constant-head and constant-salinity boundaries were used to represent an ocean body. The constant-head was set to 0.0 m at the sea boundary on the left side (x = 0 m), see Figure 2, bottom left. A constant initial condition of the salt concentration was fixed to 37   kg / m 3 on the left side and uniformly distributed through the entire left face (x = 0 m). The constant head and concentration inland boundary condition on the right face (x = 10,000 m) is 10.0 m as the elevation of the spring and 0.0   kg / m 3 as uncontaminated freshwater. Figure 3 shows the groundwater flow gradient that goes from 0.0 m to 10.0 m as a result of the MODFLOW steady-state groundwater simulation (see the colored legend). The bottom (z = 0 m) and the lateral surfaces (y = 0 m and y = 1000 m) were set as no flow boundaries.
The values of the density fluid in the SEAWAT VDF packages go from 10 3   kg / m 3 for freshwater up to 1026   kg / m 3 , which corresponds to a typical seawater density. For the three-dimensional numerical simulations run of transient 20-day stress in the SEAWAT model is evaluated and then used as a starting point of a new simulation. This procedure is repeated about 20 times for a longer simulation which corresponds to 400 days.

3.3. Sensitivity Analysis

The sensitivity analysis evaluates the uncertainties of salinity and head simulations to a small variation (not bigger than 5%) of the different parameters of Table 1. It is therefore local in the parameter space. This analysis may help to understand the effects of variations and interactions of the aquifer parameters and boundary conditions on the numerical simulations. We used four different parameters in the numerical model and performed 20 numerical simulations for the local sensitivity analysis, where three of them correspond to the groundwater flow model, including hydraulic conductivities (HCf, HCc), vertical anisotropy (VA) and one corresponds to the solute transport model, the salt concentration at the boundary condition (SC). For the sensitivity analysis we used the values in agreement with Reference [21] for the effective porosity (POp, POc), and the longitudinal dispersivity (LDf, LDc).
Furthermore, in the local sensitivity analysis, the values of the composite scaled sensitivity (Css) of the parameters for head and salinity simulations were computed at several locations along both conduits network and several locations in the fractured medium. The Css were computed for the parameter values in the maximum seawater intrusion. Parameter sensitivities were calculated at several locations from column n. 2 to column n. 13 (see Figure 3, bottom left) along both conduits (layers 10 and 16, respectively and row 16), where Δ r j = 25   m   long in the horizontal axes. Column n.1 is close to the shoreline and has a constant concentration of 37   kg / m 3 (column n.400 instead, has a salt concentration of 0   kg / m 3 ). See Figure 4. The parameter sensitivities of the numerical simulations in the fractured medium were evaluated in two different locations. The first set at layer 6, 13 and 21, row 16 and along the columns from 2 to 13, just above, between the conduits and below them. The second set is similar to the previous one but in row 25, instead of row 16, that is, away from both conduits. See Figure 5.

3.3.1. Local Sensitivity Analysis of Numerical Simulations in the Conduits

Figure 4a shows the salinity Css values of all parameters for the numerical simulations along both conduits (layer 10, blue rectangle and layer 16, red one) in the local sensitivity analysis. In general, we noticed that the largest Css value corresponds to the salt concentration, SC, in the second conduit. Indeed, the effect of a variation of SC is maximum in the second conduit. This is the most important parameter for the salinity profile. Also, the hydraulic conductivity in the fractured medium, HCf, is a very important parameter for the salinity profile in the second conduit. While the HCc and the vertical anisotropy (VA) have intermediate values. Concerning the upper conduit, the largest value of Css corresponds to SC while both, HCf and HCc have similar intermediate values. Our results show that the second conduit is more sensitive than the upper conduit for the salinity profile.
For simplicity, we have considered a system with a dual-permeability aquifer and thus, we considered the same value in both conduits. Results on the Css for the head simulations in Figure 4b are similar to those obtained for the salinity in Figure 4a, but the values of Css are much smaller. In general, salinity observations are more effective in the sensitivity analysis, than hydraulic head observations, as already noticed in Reference [25].

3.3.2. Local Sensitivity Analysis of Numerical Simulations in the Fractured Medium

Figure 5 shows the values of Css for all parameters for salinity (a) and head (b) simulations calculated in the evaluated locations in the fractured medium (layer n. 6, 13, 21) and column j = 16 (above, between and below the two conduits). This Figure shows that the most important parameter is the salt concentration SC, for both, the salinity (a) and the head (b) in all the fractured medium, but also the hydraulic conductivity HCf is an important parameter. The largest value of Css indicates that SC below the conduits is the most sensitive parameter in both salinity and head numerical simulations (and also for HCf). Below both conduits, VA has intermediate Css values.
Figure 6 is similar to Figure 5 but the Css values were computed in the fractured medium, away from both conduits, more specifically at row 25 (j = 25) instead of j = 16. Their behaviour is similar in both salinity and head Css and indicates that the saltwater concentration and the hydraulic conductivity in the fractured medium are the most important parameters. When comparing with Figure 5, the vertical anisotropy is a more important parameter near the conduits. The Css values are bigger in the region below both conduits where the salt wedge is more pronounced. That means that the conduit systems have a significant impact on the numerical simulations results.
Figure 7 shows the Css composite scaled sensitivity values of selected parameters at different locations along with the first (a) and second (b) conduits as a function of the column number along the conduits. The largest Css value is found at i = 2 for the first conduit (a) and i = 6 for the second conduit (b). These locations are found at the mixing zone more than anywhere else because the salinity and head profile change significantly near the mixing zone and remains almost constant in other locations.
In Figure 8, instead, the largest Css values of selected parameters at different locations in the fractured medium are located in i = 2, 3, 10, respectively.
The computed correlation parameters coefficients and covariance matrix of all parameters are shown in Figure 9. In general, we observed that hydrological parameters of the fractured medium are positively correlated with the other parameters of the fractured medium but negatively correlated with the conduit parameters. For example, the first square on the top left side (first conduit) indicates that for example, HCf is positively correlated with itself and negatively correlated with HCc. As soon as one goes in-depth this pattern changes. Below conduits (last square first row) the correlation between HCf and HCc becomes positive. The same happens away conduits. For example, below (away) conduits, the correlation coefficient between HCf and HCc increases positively to the maximum.
In Figure 10 we show the nonlinear relationship between salinity and head for the parameters SC (salt concentration), HCc (hydraulic conductivity in the conduits, HCf (hydraulic conductivity in the fractured medium) using the locations where the Css is maximum. This nonlinearity of the derivative of the simulations with respect to the selected parameters clearly shows that the local sensitivity results are not representative of the entire parameter range (of Table 1). One reason could be that the Darcy equation which describes a laminar flow is also used to calculate the flow velocity inside both conduits where the flow is non-laminar. Similar behavior is shown in Figure 11 where this nonlinear relationship is also observed far away from both conduits.

3.4. Numerical Simulations Results of Seawater Intrusion Scenarios

After having determined the important parameters of the model, i.e., SC and HCf inside the second conduit, and using the values of the porosity and longitudinal dispersion in both, conduits and fractured medium from [21], where a series of transient simulations were performed in order to calibrate the model and determine the best choice of parameters, we investigated various scenarios of saltwater wedge shape and examined which were the parameters that most influence the extension of the wedge.
These results are presented in Figure 12 where we showed the cross-section salinity profile (at row 16) as a function of the depth, at different distances from the sea (with a constant head of 0.0 m) after a simulation of approximately 140 days. The bold red line corresponding to the salinity profile at a distance of 150 m from the sea, (i = 5), highlights how the salinity profile changes its shape as a function of the salt concentration. The variation is computed for different values of SC at the boundary, which goes from 37 kg/m3 (Figure 12a) to 10 kg/m3 (Figure 12e). A value of 10 kg/m3 is not present in nature and we used it to understand how the shape changes from the maximum value of salinity to almost freshwater. It is also interesting to notice how the two-steps like pattern disappear when the salt concentration decreases. In this Figure, we also depicted the position of the two karst conduits, which are represented as two parallel blue lines of 1 m thick each.
As can be seen from Figure 12, the salt concentration of 37 kg/m3 (a), three different step-like shapes may be identified in the salinity profile. One at the left sea boundary of the figure (x = 0) and the other two in concomitance with both conduits, similar to those appearing in Figure 2. Apparently, the presence of the two karst conduits prevent the formation of the typical salt wedge intrusion (as in Figure 13 where the salinity profile is computed away from both conduits, in row 25) and, instead, shows a step-like shape in which the freshwater coming from the aquifer conduit push away the seawater that in any case may enter into the karst conduits more or less, depending on the calibration of the different parameters such as, the hydraulic conductivity or the initial value of the heads.
We also investigated the effects on the depth and the value of the constant head in the inland. Figure 14a shows the differences between the salinity profile and wedge when the constant-head vary from 10.0 m (red line) to 8.0 m (blue line) at the inland as a function of the depth. Here we can notice that the salt intrusion is more pronounced for the red line case rather than the blue one, at the same position from the sea.
Figure 14b shows a similar step-like shape behaviour when the grid is extended from −45 m to −220 m. Here we increased the number of layers from 26 to 35 in the cross-section. In this case, the salt intrusion is overall more pronounced than the case of Figure 14a and, in particular, the blue line is ahead of the red one.
Figure 15 compares the data of Figure 14 with the Electrical conductivity (EC) data (as an indirect measure of the salinity) observed in well 18, for both scenarios with different depths. The EC data is superimposed with a different X-axis scale to the modeled salt concentration. The model corresponding to the numerical simulation of the case (Figure 15a) (red line) can reproduce the two-steps shape of the observed EC. These results support the hypothesis of two conduits at about −14 m and −24 m asl taken during the development of the conceptual model. The case (Figure 15b) does not fit very well the EC data.
We analyze the data resulting from a numerical simulation of salt intrusion similar to the case of Figure 15b when the whole grid is filled up with a salt concentration of 37 kg / m 3 . A comparison with the EC profile is presented in Figure 16a and a zoom-in (b). Here we show a cross-section at row 16 and column i = 305.
Since the salt concentration fills up the whole grid and because the groundwater flows in the direction of the sea, the salt is pushed away from the inland toward the sea and we observed a two-steps shape. The cross-section in Figure 16 is situated at a distance of about 2375 m from the inland.

4. Conclusions

In this paper, we performed a series of three-dimensional numerical simulations using the SEAWAT code to study seawater intrusion in a dual-permeability coastal karst aquifer with a double conduit network according to other recent studies coupling density-dependent flow and transport models [19]. The numerical model is based on the conceptual model of Reference [21], recently published and the parameter values of a fractured and karstified carbonate aquifer of the Cuban land western territory. Salt intrusion in the presence of karst conduits is a three-dimensional problem, especially in the zone close to the conduits.
We performed a sensitivity analysis to evaluate the impact of the parameters and boundary conditions on the modeling saltwater wedge in a karstic aquifer, including hydraulic conductivity, both, in the conduits and the fractured medium, vertical anisotropy, salinity concentration at the boundary conditions. Our results indicate that salt concentration and hydraulic conductivity of the fractured medium in the second conduit are the most important parameters, but also the hydraulic conductivity in the conduits plays an important role.
The model successfully describes the variation of the Electrical conductivity EC (as an indirect measure of the salinity concentration) observed in a well as a function of the depth and it is based on the parameter values of a fractured aquifer with two karst conduits with large value of both, the hydraulic conductivity and the effective porosity. The simple numerical model reproduces the two-steps shape of the observed EC profile. These results support the hypothesis of two conduits at about −14 m a.s.l. and −24 m a.s.l. taken during the development of the conceptual model.
We also showed results of three-dimensional numerical simulations on groundwater salinity for different scenarios where the boundary conditions and depth of the aquifer are changed, using the important parameters and compare results with the Electrical Conductivity (EC) profiles measured in a well as a function of the depth. The numerical model further demonstrated that karst conduits play an important role in influencing the distribution of water salinity.

Author Contributions

Conceptualization was developed by A.F., A.Z., E.P., R.H.-D. and F.C.; methodology was carried out by A.F., A.Z. and F.C.; software development, modelling and validation were conducted by A.F.; data collection and investigation process were carried out by R.H.D. and F.C.; data curation was managed by R.H.D. and F.C.; writing-original draft preparation, review and editing were realized by A.F., A.Z., E.P., R.H.-D. and F.C.; supervision was fulfilled by F.C.; project administration and funding acquisition were achieved by F.C.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Acknowledgments

The authors would like to thank the two anonymous reviewers for their fruitful comments.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fleury, P.; Bakalowicz, M.; de Marsily, G. Submarine springs and coastal karst aquifer: A review. J. Hydrol. 2007, 339, 79–92. [Google Scholar] [CrossRef]
  2. Bear, J. Seawater Intrusion in Coastal Aquifers; Springer: Berlin/Heidelberg, Germany, 1999. [Google Scholar]
  3. Werner, C.L.; Bakker, M.; Post, V.E.; Vandenbohede, A.; Lu, C.; Ataie-Ashtiani, B.; Simmons, C.T.; Barry, D.A. Seawater intrusion processes, investigations and management: A recent advances and future challenges. Adv. Water Resour. 2013, 51, 3–26. [Google Scholar] [CrossRef]
  4. Voss, C.I.; Souza, W.R. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone. Water Resour. Res. 1987, 23, 1851–1866. [Google Scholar] [CrossRef]
  5. Moore, W.S.; Wilson, A.M. Advective flow through the upper continental shelf driven bt storms, buoyancy, and submarine groundwater discharge. Earth Planet Sci. Lett. 2005, 235, 564–576. [Google Scholar] [CrossRef]
  6. Voss, C.I.; Provost, A.M. SUTRA, US Geological Survey Water Resources Investigation Reports; 84-4369; USGS: Reston, VA, USA, 1984.
  7. Diersch, H.J.G. FEFLOW Reference Manual; Institute for Water Resources Planning and Systems Research Ltd.: Berlin, Germany, 2002; p. 278. [Google Scholar]
  8. Harbaugh, A.W.; Banta, E.R.; Hill, M.C.; McDonald, M.G. MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water-User Guide to Modularization Concepts and the Ground-Water Flow Process: U.S. Geological Survey Open-File Report 00-92; U.S. Geological Survey: Reston, VA, USA, 2000; p. 121.
  9. Guo, W.; Langevin, C. Users’ Guide to SEAWAT: A Computer Program for Simulation of Three-Dimensional Variable-Density Ground-Water Flow; Water Resources Investigations Report; United States Geological Survey: Reston, VA, USA, 2002.
  10. Langevin, C.; Shoemaker, W.B.; Guo, W. MODFLOW-2000, the US Geological Survey Modular Ground-Water Model Documentation of the SEAWAT-2000 Version with the Variable-Density Flow Process (VDF) and the Integrated MT3DMS Transport Process (IMT); US Department of the Interior, US Geological Survey: Reston, VA, USA, 2003.
  11. Zheng, C.; Wang, P.P. MT3DMS—A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion and Chemical Reactions of Contaminants in Ground-Water Systems; Documentation and User’s Guide: U.S. Army Corps of Engineers Contract Report SERDP-99-1; University of Alabama: Tuscaloosa, AL, USA, 1999. [Google Scholar]
  12. Zheng, C. MT3DMS v5.2 Supplemental User’s Guide: TECHNICAL Report to the U.S. Army Engineer Research and Development Center, Department of Geological Sciences; University of Alabama: Tuscaloosa, AL, USA, 2006; p. 24. [Google Scholar]
  13. Davis, J.H. Hydraulic Investigation and Simulation of Groundwater Flow in the Upper Floridan Aquifer of North Central Florida and Southwestern Georgia and Delineation of Contributing Areas for Selected City of Tallahassee, Florida; 95-4296; USGS: Reston, VA, USA, 1966; p. 55.
  14. Shoemaker, W.B.; Kuniansky, E.L.; Birk, S.; Bauer, S.; Swain, E.D. Documentation of a conduit flow process (CFP) for MODFLOW-2005. U.S. Geological Survey. Tech. Methods 2008, 6, 50. [Google Scholar]
  15. Gallegos, J.J.; Hu, B.X.; Davis, H. Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP. Hydrogeol. J. 2013, 21, 1749–1760. [Google Scholar] [CrossRef]
  16. Reimann, T.; Giese, M.; Geyer, T.; Liedl, R.; Maréchal, J.C.; Shoemaker, W.B. Representation of water abstraction from a larst conduit with numerical discrete-continuum models. Hydrol. Earth Syst. Sci. 2014, 18, 227–241. [Google Scholar] [CrossRef]
  17. Xu, Z.; Hu, B.X.; Davis, H.; Cao, J. Simulating long term nitrate-N contamination processes in the Wookville Karst Plain using CFPv2 with UMT3D. J. Hydrol. 2015, 524, 72–88. [Google Scholar] [CrossRef]
  18. Xu, Z.; Hu, B.X.; Davis, H.; Kish, S. Numerical study of groundwater flow cycling controlled by seawater/freshwater interaction in a coastal karst aquifer through conduit network using CFPv2. J. Contam. Hydrol. 2015, 182, 131–145. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Xu, Z.; Hu, B.; Ye, M. Numerical modeling and sensitivity analysis of seawater intrusion in a dual permeability coastal karst aquifer with conduit networks. Hydrol. Earth Syst. Sci. 2018, 22, 221–239. [Google Scholar] [CrossRef]
  20. Xu, Z.; Hu, B.X. Development of a discrete-continuum VDFST-CFP numerical model for simulating seawater intrusion to a coastal karst aquifer with a conduit system. Water Resour. Res. 2017, 53, 688–711. [Google Scholar] [CrossRef]
  21. Hernandez-Diaz, R.; Petrella, E.; Bucci, A.; Naclerio, G.; Feo, A.; Sferra, G.; Chelli, A.; Zanini, A.; Gonzales-Hernandez, P.; Celico, F. Integrating hydrogeological and microbiological data and modelling to characterize the hydraulic features and behaviour of coastal carbonate aquifers: A case in Western Cuba. Water 2019, 11, 1989. [Google Scholar] [CrossRef]
  22. Hill, M.; Tiedemann, C. Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertaninty; Wiley: Hoboken, NJ, USA, 2006; p. 480. [Google Scholar]
  23. Hill, M.C. Methods and Guidelines for Effective Model Calibration; Report 98-4005; U.S. Geological Survey Water-Resources Investigations: Reston, VA, USA, 1998.
  24. Foglia, L.; Hill, M.C.; Mehl, S.W.; Burlando, P. Sensitivity Analysis calibration, and testing of a distributed hydrological model using error-based weighting and one objective function. Water Resour. Res. 2000, 45, W06427. [Google Scholar] [CrossRef]
  25. Shoemaker, W.B. Important observations and parameters for a salt water intrusion model. Ground Water 2004, 42, 829–840. [Google Scholar] [PubMed]
  26. Xu, Z.; Hu, B.; Xu, Z.; Wu, X. Simulating seawater intrusion in a complex coastal karst aquifer using an improved variable-density flow and solute transport-conduit flow process model. Hydrogeol. J. 2019, 27, 1277–1289. [Google Scholar] [CrossRef]
  27. Fagundo Castillo, J.R.; Gonzales Hernandez, P. Agricultural use and water quality at karstic Cuban western plain. Int. J. Speleol. 1999, 28, 175–185. [Google Scholar] [CrossRef] [Green Version]
  28. Molerio León, L.; Parise, M. Managing environmental problems in Cuban karstic aquifers. Environ. Geol. 2009, 58, 275–283. [Google Scholar] [CrossRef]
  29. Hernández, R.; Ramírez, R.; López-Portilla, M.; González, P.; Antigüedad, I.; Díaz, S. Seawater Intrusion in the Coastal Aquifer of Guanahacabibes, Pinar del Río, Cuba. In Management of Water Resources in Protected Areas; Farfán González, H., Corvea Porras, J.L., de Bustamante Gutiérrez, I., LaMoreaux, J.W., Eds.; Springer: Berlin/Heidelberg, Germany, 2013; pp. 301–308. [Google Scholar]
  30. Bucci, A.; Petrella, E.; Celico, F.; Naclerio, G. Use of molecular approaches in hydrogeological studies: The case of carbonate aquifers in southern Italy. Hydrogeol. J. 2017, 25, 1017–1031. [Google Scholar] [CrossRef]
  31. Petrella, E.; Falasca, A.; Celico, F. Natural-gradient tracer experiment in epikarst: A test study in the Acqua dei Faggi experimental site, southern Italy. Geofluids 2008, 8, 159–166. [Google Scholar] [CrossRef]
  32. Bakker, M.; Post, V.; Langevin, C.D.; Hughes, J.D.; White, J.T.; Starn, J.J.; Fienen, M.N. Scripting MODFLOW Model Development Using Python and FloPy. Groundwater 2016, 54, 733–739. [Google Scholar] [CrossRef]
Figure 1. Sketch view of the karst system.
Figure 1. Sketch view of the karst system.
Water 11 02311 g001
Figure 2. Electrical Conductivity (EC) variation with depth in “well 18” [21].
Figure 2. Electrical Conductivity (EC) variation with depth in “well 18” [21].
Water 11 02311 g002
Figure 3. Three-dimensional grid used for the numerical simulations. The two red parallel lines correspond to the karst conduits. This figure was generated using FloPy [32].
Figure 3. Three-dimensional grid used for the numerical simulations. The two red parallel lines correspond to the karst conduits. This figure was generated using FloPy [32].
Water 11 02311 g003
Figure 4. (a) Calculated composite scaled sensitivity (Css) values of all parameters (HCf, HCc, SC and VA) with respect to salinity numerical simulations in the conduit (layer n. 10, 16) and j = 16, in the local sensitivity analysis; (b) Calculated Css values of all parameters with respect to head numerical simulations in the conduit (same as (a)). HCf, HCc, SC and VA represent the hydraulic conductivity (fractured medium), the hydraulic conductivity (conduits), salt concentration and vertical anisotropy, respectively.
Figure 4. (a) Calculated composite scaled sensitivity (Css) values of all parameters (HCf, HCc, SC and VA) with respect to salinity numerical simulations in the conduit (layer n. 10, 16) and j = 16, in the local sensitivity analysis; (b) Calculated Css values of all parameters with respect to head numerical simulations in the conduit (same as (a)). HCf, HCc, SC and VA represent the hydraulic conductivity (fractured medium), the hydraulic conductivity (conduits), salt concentration and vertical anisotropy, respectively.
Water 11 02311 g004
Figure 5. Calculated Css values of all parameters with respect to: (a) Salinity simulations in the fractured medium (layer n. 6, 13, 21) and row 16; (b) Head simulations in the fractured medium (layer n. 6, 13, 21) and row 16.
Figure 5. Calculated Css values of all parameters with respect to: (a) Salinity simulations in the fractured medium (layer n. 6, 13, 21) and row 16; (b) Head simulations in the fractured medium (layer n. 6, 13, 21) and row 16.
Water 11 02311 g005
Figure 6. Calculated Css values of all parameters with respect to: (a) Salinity simulations in the fractured medium, away (j = 25) from de two conduits (layer n. 6, 13, 21); (b) Head simulations in the fractured medium, away (j = 25) from the two conduits (layer n. 6, 13, 21).
Figure 6. Calculated Css values of all parameters with respect to: (a) Salinity simulations in the fractured medium, away (j = 25) from de two conduits (layer n. 6, 13, 21); (b) Head simulations in the fractured medium, away (j = 25) from the two conduits (layer n. 6, 13, 21).
Water 11 02311 g006
Figure 7. Calculated Css values of selected parameters at different locations along the: (a) First conduit; (b) Second conduit.
Figure 7. Calculated Css values of selected parameters at different locations along the: (a) First conduit; (b) Second conduit.
Water 11 02311 g007
Figure 8. Calculated Css values of selected parameters at different locations in the fractured medium: (a) above the first conduit; (b) between the two conduits; (c) below conduits, and row 25.
Figure 8. Calculated Css values of selected parameters at different locations in the fractured medium: (a) above the first conduit; (b) between the two conduits; (c) below conduits, and row 25.
Water 11 02311 g008
Figure 9. Correlation coefficient matrix for all four parameters in different locations of the aquifer: conduits, fractured medium near and away conduits.
Figure 9. Correlation coefficient matrix for all four parameters in different locations of the aquifer: conduits, fractured medium near and away conduits.
Water 11 02311 g009
Figure 10. The nonlinear relationship between: (a) salinity and (b) head simulations with respect to the parameters SC (Salt Concentration) for both conduits.
Figure 10. The nonlinear relationship between: (a) salinity and (b) head simulations with respect to the parameters SC (Salt Concentration) for both conduits.
Water 11 02311 g010
Figure 11. Nonlinear relationship between: (a) Salinity and salt concentration; (b) Salinity and HDf both in the fractured medium away from the conduits using the locations where the Css are maximum.
Figure 11. Nonlinear relationship between: (a) Salinity and salt concentration; (b) Salinity and HDf both in the fractured medium away from the conduits using the locations where the Css are maximum.
Water 11 02311 g011
Figure 12. Salinity profile at row j = 16 as a function of the depth at different distances from the sea (0.0m). The blue line, the red bold line and the orange line correspond to a distance of 50 m, 150 m and 325 m, respectively, from the sea boundary. And a salt concentration at the boundary of: (a) 37 kg / m 3 ; (b) 35 kg / m 3 ; (c) 30 kg / m 3 ; (d) 20 kg / m 3 ; (e) 10 kg / m 3 .
Figure 12. Salinity profile at row j = 16 as a function of the depth at different distances from the sea (0.0m). The blue line, the red bold line and the orange line correspond to a distance of 50 m, 150 m and 325 m, respectively, from the sea boundary. And a salt concentration at the boundary of: (a) 37 kg / m 3 ; (b) 35 kg / m 3 ; (c) 30 kg / m 3 ; (d) 20 kg / m 3 ; (e) 10 kg / m 3 .
Water 11 02311 g012
Figure 13. Salinity profile away from the conduits (j = 25), to be compared with the previous Figure 10a. In this case the scenario is similar to a salt intrusion without the two-steps pattern, that is, the presence of both conduits does not affect the salinity profile when it is measured away from them.
Figure 13. Salinity profile away from the conduits (j = 25), to be compared with the previous Figure 10a. In this case the scenario is similar to a salt intrusion without the two-steps pattern, that is, the presence of both conduits does not affect the salinity profile when it is measured away from them.
Water 11 02311 g013
Figure 14. Differences between the salinity profile and wedge intrusion as a function of the depth with a constant-head equal to 10.0 m (red line) and 8.0 m (blue line) for an aquifer of: (a) −45 m a.s.l.; (b) −220 m a.s.l. Notice the presence of the two steps-shape in both scenarios.
Figure 14. Differences between the salinity profile and wedge intrusion as a function of the depth with a constant-head equal to 10.0 m (red line) and 8.0 m (blue line) for an aquifer of: (a) −45 m a.s.l.; (b) −220 m a.s.l. Notice the presence of the two steps-shape in both scenarios.
Water 11 02311 g014
Figure 15. Comparison of the data of Figure 14 with the EC data in well 18 of Reference [21] using two different values of the constant head in the inland, for and aquifer of (a) −45 m a.s.l.; (b) −220 m a.s.l.
Figure 15. Comparison of the data of Figure 14 with the EC data in well 18 of Reference [21] using two different values of the constant head in the inland, for and aquifer of (a) −45 m a.s.l.; (b) −220 m a.s.l.
Water 11 02311 g015
Figure 16. (a) Comparison of the EC profile with a numerical simulation where the whole three-dimensional grid is filled up with a salt concentration of 37 kg / m 3 . The cross-section is at row 16 and the column 305; (b) a zoom in the region of interest. Notice that due to the hydraulic conductivity the salt is pushed away in the direction toward the sea. It is observed the two-steps shape.
Figure 16. (a) Comparison of the EC profile with a numerical simulation where the whole three-dimensional grid is filled up with a salt concentration of 37 kg / m 3 . The cross-section is at row 16 and the column 305; (b) a zoom in the region of interest. Notice that due to the hydraulic conductivity the salt is pushed away in the direction toward the sea. It is observed the two-steps shape.
Water 11 02311 g016
Table 1. Definitions of the parameters used in this paper, the specific evaluated values in the local sensitivity analysis and evaluation ranges (the lower and upper ones) of each parameter in the global sensitivity analysis.
Table 1. Definitions of the parameters used in this paper, the specific evaluated values in the local sensitivity analysis and evaluation ranges (the lower and upper ones) of each parameter in the global sensitivity analysis.
ParameterDefinitionLowerUpperEvaluated ValueUnits
HCfHydraulic conductivity (fractured medium)10.03000.01045.0 m / d ay
HCcHydraulic conductivity (conduits) 1.0 × 10 6 3.0 × 10 6 2.4 × 10 6 m / day
VAVertical anisotropy1.010.01.0dimensionless
POfPorosity (fractured medium) 0.1dimensionless
POcPorosity (conduits) 1.0dimensionless
LDfLongitudinal dispersivity (fractured medium) 10.0 m
LDcLongitudinal dispersivity (conduits) 0.3 m
SCSalinity concentration10.037.037.0 kg / m 3
SHSpecified head boundary conditions8.010.010.0 m

Share and Cite

MDPI and ACS Style

Feo, A.; Zanini, A.; Petrella, E.; Hernàndez-Diaz, R.; Celico, F. Analysis of the Saltwater Wedge in a Coastal Karst Aquifer with a Double Conduit Network, Numerical Simulations and Sensitivity Analysis. Water 2019, 11, 2311. https://doi.org/10.3390/w11112311

AMA Style

Feo A, Zanini A, Petrella E, Hernàndez-Diaz R, Celico F. Analysis of the Saltwater Wedge in a Coastal Karst Aquifer with a Double Conduit Network, Numerical Simulations and Sensitivity Analysis. Water. 2019; 11(11):2311. https://doi.org/10.3390/w11112311

Chicago/Turabian Style

Feo, Alessandra, Andrea Zanini, Emma Petrella, Rebeca Hernàndez-Diaz, and Fulvio Celico. 2019. "Analysis of the Saltwater Wedge in a Coastal Karst Aquifer with a Double Conduit Network, Numerical Simulations and Sensitivity Analysis" Water 11, no. 11: 2311. https://doi.org/10.3390/w11112311

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop