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

: 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 ﬂow 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 e ﬀ ect on saltwater intrusion in a karstic aquifer. We also show results of the three-dimensional numerical simulations on groundwater salinity for di ﬀ erent scenarios with the variabilities of the important parameters and compare results with electric conductivity proﬁles measured in a well.


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.

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: 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 dss values that satisfy the following equation: where dss ij is the dimensionless scaled sensitivity of the i-th simulation with respect to the j-th parameter, and ω ii is the weight of the i-th simulation, given by 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 dss values of different simulations with respect to each parameter are accumulated as the composite scaled sensitivity (css) values and provide the total amount of information provided by the simulation for the estimation of one parameter. The css of the j-th parameter is given by where ND 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 css values identify parameters for which more information is provided by all observations for each parameter. Correlation coefficients, cor(i, j), are defined as: where cov(i, j) is the covariance between parameter i and j; var( 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.

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

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 and a different hydraulic conductivity with respect to the fractured medium. 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.

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.

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/m 3 ), respectively.

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.

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 km 2 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. 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].

Initial and Boundary Conditions
Constant-head and constant-salinity boundaries were used to represent an ocean body.

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

Initial and Boundary Conditions
Constant-head and constant-salinity boundaries were used to represent an ocean body. 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.

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

Local Sensitivity Analysis of Numerical Simulations in the Fractured Medium
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.   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. 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. 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. 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 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. 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.

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/m 3 (a) to 10 kg/m 3 (e)). A value of 10 kg/m 3 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.

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/m 3 (Figure 12a) to 10 kg/m 3 (Figure 12e). A value of 10 kg/m 3 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.
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/m 3 (a) to 10 kg/m 3 (e)). A value of 10 kg/m 3 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/m 3 (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. As can be seen from Figure 12, the salt concentration of 37 kg/m 3 (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. 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.
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 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.
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 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.
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.   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 . 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. 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 twosteps shape. The cross-section in Figure 16 is situated at a distance of about 2375 m from the inland.

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

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.