A New Method to Infer Advancement of Saline Front in Coastal Groundwater Systems by 3D: The Case of Bari (Southern Italy) Fractured Aquifer

Abstract

A new method to study 3D saline front advancement in coastal fractured aquifers has been presented. Field groundwater salinity was measured in boreholes of the Bari (Southern Italy) coastal aquifer with depth below water table. Then, the Ghyben-Herzberg freshwater/saltwater (50%) sharp interface and saline front position were determined by model simulations of the freshwater flow in groundwater. Afterward, the best-fit procedure between groundwater salinity measurements, at assigned water depth of 1.0 m in boreholes, and distances of each borehole from the modelled freshwater/saltwater saline front was used to convert each position (x, y) in groundwater to the water salinity concentration at depth of 1.0 m. Moreover, a second best-fit procedure was applied to the salinity measurements in boreholes with depth z. These results provided a grid file (x, y, z, salinity) suitable for plotting the actual Bari aquifer salinity by 3D maps. Subsequently, in order to assess effects of pumping on the saltwater-freshwater transition zone in the coastal aquifer, the Navier-Stokes (N-S) equations were applied to study transient density-driven flow and salt mass transport into freshwater of a single fracture. The rate of seawater/freshwater interface advancement given by the N-S solution was used to define the progression of saline front in Bari groundwater, starting from the actual salinity 3D map. The impact of pumping of 335 L·s⁻¹ during the transition period of 112.8 days was easily highlighted on 3D salinity maps of Bari aquifer.

1. Introduction

Groundwater over-abstractions typically can lower water table level and reduce freshwater fluxes, leading to severe saltwater intrusion problems in several coastal and metropolitan areas [1,2,3]. However, water supply is necessary for regional economic development and even for energy production. For instance, renewable sources of energy associated with the salinity gradient can be recovered in coastal areas where freshwater is mixed with saltwater. The techniques to produce electric energy from salinity gradient use the pressure-retarded reverse osmosis process [4] and associated conversion technologies. In addition, also geothermal energy production at low enthalpy should also increase at a rate of 3%–4% in Italy until 2020 [5]. Both geothermal low-enthalpy technology and the pressure-retarded reverse osmosis method require freshwater supplies (i.e., pumping) from coastal aquifers.

Subsequently, the subtraction of appreciable freshwater volumes from coastal aquifers cannot be always avoided and its real impact on groundwater should be adequately investigated. The salinity maps are useful tools for coastal groundwater management and they show how inland advancements of seawater may affect coastal aquifer quality. Chongo et al. [6] carried out an experimental work in Africa to define groundwater salinity variation map within the sedimentary formations of the Barotse sub-basin in the Western Province of Zambia. 39 boreholes were used to construct a database for the geological model development in this study area. The GEOSCENE 3D (I-GIS, 2016, I-GIS, Risskov, Denmark) [7] software was employed for visualizing, interpreting, editing and publishing geological data. Data derived from field measurements were analyzed using the SiTEM SEMDI (Hydrogeophysics-Group-1, 2001, University of Aarhus, Aarhus, Denmark) [8] and Res2Dinv (Geomoto, Penang, Malaysia) [9] software. In the study area the distribution of most boreholes was somewhat unidirectional and it resulted in data gaps that make geo-modeling quite difficult. For this reason, authors proposed the use of pseudo (i.e., imaginary) boreholes in order to reduce the data gaps. Other researchers [10] suggest the employment of empirical models such as the Archie’s law (1942) [11] to improve the best fit of experimental data (i.e., water electrical conductance) collected from wells. Lesch et al. [12] demonstrated the efficiency of the regression-based statistical method for predicting, at field scale, soil spatial salinity conditions, starting from rapidly acquired electromagnetic induction data. This regression model incorporates multiple measurements and their trends in order to increase the prediction accuracy of soil salinity maps.

In the present work a new method is presented for a rapid 3D visualization of saltwater-freshwater transition in coastal zone under new pumping stress conditions. The method requires field salinity measurements in boreholes and simulation results from two different flow models. The first model studies the steady groundwater flow in fractures by providing the Ghyben-Herzberg seawater/freshwater sharp interface toe position along the coast. Fracture transmissivity were derived from pumping and tracer (chlorophyllin) injection tests. The second flow model studies the transient density-driven flow of seawater into freshwater of a single fracture by means of the N-S equations. The result of the latter model was used to predict the transient advancement of saline front in fractures subjected to new pumping conditions.

2. Methodology

Groundwater flow modeling was addressed in a 3D set of horizontal and parallel fractures characterized by impermeable rock walls (Figure 1). Each horizontal fracture of the 3D set has a variable aperture in the x-y plane [13] of 7240 × 6300 m². Pumping and/or injection tests carried out into Bari wells, allowed the identification of the average aperture of all fractures of the set at each well position in order to determine the experimental variogram of spatial aperture covariance. Experiments on model calibrations and validations have been extensively reported by Masciopinto [14], Masciopinto et al. [13,15], and Masciopinto and Palmiotta [16]. The flow solution given by the flow model yields the freshwater discharge along the top border of domain. This freshwater discharge was used to determine the Ghyben-Herzberg sharp (50%) freshwater/saltwater interface [14]. These results defined also the saline front position along the coast.

At the second step (Figure 2), the groundwater salinity into a generic borehole of the Bari aquifer was fitted as a function of the distance of the same borehole from the saline front. For this interpolation, salinity measurements in 25 boreholes of the Bari coastal aquifer (Southern Italy) at the same depth of 1.0 m below the water table were used. The resulting best-fit equation was used to convert the generic grid node position of the computational domain with respect to the saline front into the groundwater salinity concentration (at the depth of 1.0 m). Moreover, a further best-fit procedure was applied to the salinity profiles collected into 17 boreholes in order to estimate the groundwater salinity concentration with depth z below water table. These groundwater concentrations (i.e., discontinuous values) were used in RockWorks15 (RockWare Inc., Golden, CO, USA) to obtain a solid model (i.e., continuous values) useful to visualize groundwater salinity 3D map at Bari coastal aquifer.

Figure 1. (a) Parallel set of horizontal (x-y) fractures with permeability equivalent to the real fractured medium; (b) modeled fracture apertures in each single fracture of the set: the aperture variation in the x-y plane is obtained by means of the experimental variogram derived from results of well pumping (or injection) tests.

Figure 1. (a) Parallel set of horizontal (x-y) fractures with permeability equivalent to the real fractured medium; (b) modeled fracture apertures in each single fracture of the set: the aperture variation in the x-y plane is obtained by means of the experimental variogram derived from results of well pumping (or injection) tests.

Figure 2. Flowchart of the applied methodology to produce transient 3D maps of the groundwater salinity.

Figure 2. Flowchart of the applied methodology to produce transient 3D maps of the groundwater salinity.

Finally, the N-S solutions into a single fracture were used to define the inland progression of saline front into groundwater during saltwater-freshwater transition period due to a simulated new pumping, starting from the actual salinity data. The results produce groundwater salinity values (and 3D map), at each selected simulation time.

An alternative method to obtain the same groundwater salt concentrations on the freshwater-saltwater transition is application of specific type of software, such as HydroGeoSphere (Aquanty Inc., Waterloo, ON, Canada), or similar [17], which allows 3D simulations of density-driven flow system under transient conditions. Anyway the application of this type of models to discretely fractured media is very challenging [18] due to complexity of input data required to describe preferential flow pathways [19] and due to the unknown positions (and apertures) of actual fractures and their intersections, especially at a regional scale (>1000 m).

3. The Case Study and Field Tests

The Bari’s coastal land portion (Puglia region, southern Italy) (Figure 3) has a catchment surface of around 10,000 × 7000 m² due to streams that flow into the Adriatic Sea. In this coastal area groundwater is confined within the limestone (Cretaceous) formation, and water flows in horizontal fractures in SW-NE direction, i.e., towards the coast. Geologically, the Bari aquifer is located in the Murgia region. In the tested area (coastal area of the Murgia region), rock permeability is low and discontinuous both horizontally and vertically. The geological sequence of layers observed during borehole drilling, from top downward, is Pleistocene sandstone (3–5 m thick), Cretaceous limestone (26 m thick), and Jurassic dolomite (>20 m thick). Murgia shows neo-tectonic sub-vertical fractures of the Mesozoic rocks that are relatively frequent and barely open or are sealed by calcspar and terra rossa (bulk density 1.26 ± 0.1 g/cm³). The Bari aquifer is not highly karstified. This means that calcite dissolution had taken place inside fractures of the Bari/IRSA aquifer, and an increase in hydraulic conductivity to 0.42 cm/s was found with respect to the quasi non-karstified limestone with 0.01–0.04 cm/s of hydraulic conductivity close (<1 km) to the sea coast.

Figure 3. Spatial distribution of wells (red circles and tringles) and flow simulation results with distance d of each well from position of Ghyben-Herzberg saline front in groundwater. The yellow lines (and polygon) refer to actual groundwater flow simulation. The blue lines (and polygon) consider the flow modifications due to new pumping of 335 L·s^-1. The blue dashed line is the contour head at 1 m.

Figure 3. Spatial distribution of wells (red circles and tringles) and flow simulation results with distance d of each well from position of Ghyben-Herzberg saline front in groundwater. The yellow lines (and polygon) refer to actual groundwater flow simulation. The blue lines (and polygon) consider the flow modifications due to new pumping of 335 L·s^-1. The blue dashed line is the contour head at 1 m.

3.1. Field Set Up

Groundwater macroscopic parameters such as transmissivity T [l²·t⁻¹] (l stand for length; t stands for time) and hydraulic conductivity K [l·t⁻¹] were determined by inverting the semi-analytical solution of Thiem’s equation. The results of 58 pumping tests and 10 tracer tests (Figure 4) carried out by IRSA on the Bari aquifer boreholes, were considered in this work. In particular, pumping (or injection) and tracer tests under undisturbed (or natural) gradient were carried out in order to estimate local values of both aquifer hydraulic transmissivity and groundwater specific discharge (see Table 1). The analytical solution of tracer (chlorophyllin) dispersion into the water was applied [13] to the breakthrough curves in Figure 4 in order to determine the horizontal water velocity in each borehole.

Figure 4. Breakthrough curves of relative tracer (chlorophyllin) concentrations during injection tests in ten boreholes of Bari aquifer.

Figure 4. Breakthrough curves of relative tracer (chlorophyllin) concentrations during injection tests in ten boreholes of Bari aquifer.

Table 1. Network of monitored wells at the Bari fractured aquifer: Pumping and tracer test results.

Well ID	Coordinates (UTM)		Tracer Test	Internal Diameter (m)	Piezometric Head (m above Sea Level) January–February 2012		Aquifer Hydraulic Transmissivity (m2/s)	Ground Water Specific Discharge (m/day)
	E	N
P10	652,892.8	4,553,271.3		0.300		4.480	0.09
L1-S	652,994.3	4,552,693.5	X	0.180	4.268	5.278	0.6 × 10−3 ± 1.8 × 10−4	1.3
P11	654,726.9	4,552,308.3		>1.0	1.251	1.611	0.02
P19	654,580.1	4,552,382.9		>1.0	0.77	1.110	0.02
P14	654,093.3	4,554,883.8		0.012	0.311	0.321	0.045
L2-S	653,252.5	4,555,151.7	X	>1.0	0.582	0.722	0.047 ± 0.01	0.5
P4	653,053.8	4,554,841.5		0.012	0.266	0.386	0.06
L3-S	652,431.0	4,554,429.8	X	>1.0	1.903	2.163	0.043 ± 0.02	2.6
P3	651,569.2	4,553,208.7		0.300	4.321	5.441	0.07
P16	652,360.0	4,553,741.7		0.014	2.658	3.538	0.09
L4-S	652,850.9	4,553,352.7	X	0.300	2.535	3.450	0.033 ± 0.03	2.2
P18	652,442.4	4,552,454.1		0.200		6.926	0.07
L5-S	647,930.7	4,551,813.2	X	0.325		33.649	3.0 × 10−3 ± 1.8 × 10−3	0.2
L8-S	652,094.5	4,551,194.2	X	0.080		8.532	0.02 ± 0.01	0.4
L7-S	652,237.3	4,550,862.3	X	0.080		7.880	0.5 × 10−2 ± 0.4 × 10−2	0.6
L6-S	651,974.4	4,550,961.1	X	0.080		8.892	0.11 ± 0.01	1.4
P13	651,750.4	4,554,936.6		0.300		0.807	0.02
L9	654,682.2	4,555,123.7		0.100		0.705	0.1
L10	654,572.7	4,555,144.1		0.100		0.167	0.01
L11	654,588.5	4,555,126.1		0.100		0.317	0.01
L12-S	654,679.5	4,555,109.1	X	0.300		0.360	0.01 ± 0.7 × 10−3	0.3
L13	654,777.0	4,555,188.8		0.100		0.418	1.2 × 10−5
L14-S	649,860.6	4,552,197.6	X	0.080		34.370	1.1 × 10−4 ± 5.3 × 10−5	0.6
L15	649,623.8	4,552,059.7		0.080		35.260	1.2 × 10−5
P2	651,595.4	4,551,934.3		0.300		6.760	0.025

X means that tracer test was also performed on the borehole.

Table 1. Network of monitored wells at the Bari fractured aquifer: Pumping and tracer test results.

Well ID	Coordinates (UTM)		Tracer Test	Internal Diameter (m)	Piezometric Head (m above Sea Level) January–February 2012		Aquifer Hydraulic Transmissivity (m2/s)	Ground Water Specific Discharge (m/day)
	E	N
P10	652,892.8	4,553,271.3		0.300		4.480	0.09
L1-S	652,994.3	4,552,693.5	X	0.180	4.268	5.278	0.6 × 10−3 ± 1.8 × 10−4	1.3
P11	654,726.9	4,552,308.3		>1.0	1.251	1.611	0.02
P19	654,580.1	4,552,382.9		>1.0	0.77	1.110	0.02
P14	654,093.3	4,554,883.8		0.012	0.311	0.321	0.045
L2-S	653,252.5	4,555,151.7	X	>1.0	0.582	0.722	0.047 ± 0.01	0.5
P4	653,053.8	4,554,841.5		0.012	0.266	0.386	0.06
L3-S	652,431.0	4,554,429.8	X	>1.0	1.903	2.163	0.043 ± 0.02	2.6
P3	651,569.2	4,553,208.7		0.300	4.321	5.441	0.07
P16	652,360.0	4,553,741.7		0.014	2.658	3.538	0.09
L4-S	652,850.9	4,553,352.7	X	0.300	2.535	3.450	0.033 ± 0.03	2.2
P18	652,442.4	4,552,454.1		0.200		6.926	0.07
L5-S	647,930.7	4,551,813.2	X	0.325		33.649	3.0 × 10−3 ± 1.8 × 10−3	0.2
L8-S	652,094.5	4,551,194.2	X	0.080		8.532	0.02 ± 0.01	0.4
L7-S	652,237.3	4,550,862.3	X	0.080		7.880	0.5 × 10−2 ± 0.4 × 10−2	0.6
L6-S	651,974.4	4,550,961.1	X	0.080		8.892	0.11 ± 0.01	1.4
P13	651,750.4	4,554,936.6		0.300		0.807	0.02
L9	654,682.2	4,555,123.7		0.100		0.705	0.1
L10	654,572.7	4,555,144.1		0.100		0.167	0.01
L11	654,588.5	4,555,126.1		0.100		0.317	0.01
L12-S	654,679.5	4,555,109.1	X	0.300		0.360	0.01 ± 0.7 × 10−3	0.3
L13	654,777.0	4,555,188.8		0.100		0.418	1.2 × 10−5
L14-S	649,860.6	4,552,197.6	X	0.080		34.370	1.1 × 10−4 ± 5.3 × 10−5	0.6
L15	649,623.8	4,552,059.7		0.080		35.260	1.2 × 10−5
P2	651,595.4	4,551,934.3		0.300		6.760	0.025

X means that tracer test was also performed on the borehole.

3.2. Experiment Methodology

During pumping (or injection) tests the flow rate was changed from 1 to 3 L·s⁻¹ and from 3 to 6 L·s⁻¹, on average. In addition water table depth was simultaneously monitored (and recorded) by means of the pressure probe (mini-log data logger, SIM Instrument SNC, Milan, Italy). The flow rate of each pumping (or injection) was kept constant for about 2 h; at the end of the test, pumping was stopped and the rise (or decline, after an injection) in water depth was recorded during the recovery period.

For each tracer injection, 1 m³ of water was marked with 500 g of chlorophyll powder commercially sold as E141 hydro soluble or sodium copper chlorophyllin, which is usually used for foodstuffs such as pickles, ice creams, candies and fruit juices, and therefore it is not dangerous for human health. After mixing water and chlorophyllin powder in a tank of 1 m³ volume, the traced mixture was injected with a constant flow rate of around 2 L·s⁻¹ into a borehole, using a submersible pump and a pipeline (3 inch internal diameter). The traced water was injected at an assigned water depth (3 m) for about 8 min. The injection pipeline was then removed and groundwater was sampled at regular time intervals via a sampling pump (Grundfos BTI/MP1, Downers Grove, IL, USA) placed at the assigned water depth of 3 m into the well. Each water sample was stored in a plastic (Polyvinyl chloride or PVC) black bottle and transported in the laboratory where samples were monitored for water absorbance. Groundwater was sampled as long as the background tracer value of 10⁻³ absorbance units (AU) [20] was again achieved in the well. During injection water depth was monitored and recorded using a pressure probe and the data logger (SIM Instrument, Milano, Italy). The measurement of tracer concentration in the borehole referred to a calibration curve previously determined in the laboratory using a spectrophotometer (HACH DR/2000, Hach Lange Srl, Lainate, Italy) at a frequency of 405 nm. Test results suggest an aquifer discretization of 80 smoothed (and parallel) fractures. In fact, as n [–] defines the effective aquifer porosity, it is the uniform ratio of the void-space per unit volume of aquifer and, in each cross-section, it will be:

$$n=\frac{{\displaystyle \sum _{i=1}^{N_f}2b_i}}{B}$$

(1)

where ∑2b_i [l] is the sum of all horizontal fracture apertures of the aquifer column with unitary horizontal area (1 × 1 m²) and thickness B [l], while N_f [–] is the total number of the fractures of the parallel set (see Figure 1). Assuming that all fractures of the set have the same aperture of 1.3 mm, i.e., ∑2b_i = N_f × 2b = n × B and for n = 0.35% and B = 30 m (see test results Table 1), Equation (1) suggests N_f = 80.

4. Groundwater Flow Simulation and Ghyben-Herzberg Interface Toe Position

At the regional scale (7200 × 6300 m²) every fracture belonging to the 3D parallel set (Figure 1) was discretized using a grid step size of Δx = Δy = 150 m (i.e., 49 × 43 grid nodes). The steady and non-uniform groundwater flow was addressed in a series of parallel and horizontal fractures, with each single fracture having a spatially variable aperture and impermeable rock matrix. The model implemented the approximated analytical radial flow solution to a well in order to determine the mean conductivity of aquifer fractures at each well position by using K = nb²/3 γ_w/µ, where γ_w [M·l⁻²·t⁻²] (M stands for mass) is the water specific weight and µ [M·l⁻¹·t⁻¹] is the dynamic water viscosity. The experimental variogram regarding the fracture apertures derived from well pumping tests was employed to perform the stochastic generation of apertures in each horizontal fracture belonging to the set. The discharge/head relationship in each fracture is defined as [15]

$$\left({\phi }_i-{\phi }_j\right)={Q_{ij}}^2\left[\frac{f}{2g\mathsf{\Delta }y}\frac{\mathsf{\Delta }x}{\mathsf{\Delta }y}\left(\frac{1}{{\left(2b_i\right)}^3}+\frac{1}{{\left(2b_j\right)}^3}\right)\right]$$

(2)

where f [–] is the friction factor derived from the Reynolds number; Q_ij [l³·t⁻¹] is the flow rate of each fracture between two generic grid nodes i and j; Δx [l] and Δy [l] are the grid steps; and φ = p/γ_w [l] is the water head, whilst g [l·t⁻²] is the gravity acceleration. It should be noted that Equation (2) supports flow calculations at any Reynolds number. By imposing the continuity equation, i.e., ∑Q_ij = 0, in every grid node, a system of equations was defined. The over-relaxation method was applied to solve the system of equations after that the boundary conditions were assigned. The flow simulation results enabled calculations of the seawater/freshwater 50% sharp interface positions with respect to the coastline, by applying the Ghyben-Herzberg equation. Indeed, to predict the interface toe position L [l] with respect to the coastline, the resulting groundwater outflow was managed in order to calculate the length of intrusion for every position along the coast defines as [14]

$$L-L_d=K\frac{B^2-H_s^2}{2{\delta }_{\gamma }\times Q_0}-L_d=n\frac{b_i^2}{3}\frac{{\gamma }_w}{\mu }\frac{{({\delta }_{\gamma }\times {\phi }_0)}^2-H_s^2}{2{\delta }_{\gamma }\times Q_0^i}-L_d$$

(3)

where L_d [l] is the distance of the contour head φ₀ (for instance 1 m) from the coastline given by flow simulation result; H_s [l] is the sharp interface depth at the outflow (usually set = 0); $Q_0^i$ [l²·t⁻¹] is the groundwater discharge along the coast predicted by model at grid node i; and ${\delta }_{\gamma }={\gamma }_w/\left({\gamma }_s-{\gamma }_w\right)$ [–] is the ratio of the specific weights.

Equation (3) shows that the extension $L-L_d$ of the sea intrusion is a function of both groundwater outflow and maximum freshwater thickness. Figure 3 saline fronts in groundwater before (in yellow) and after (in blue) a simulated new pumping of 335 L·s⁻¹ (assuming that reinjections of pumped water are not allowed). Figure 3 shows how the new pumping changes the saline front position with respect to the coast. Simulation results provided a total decreases of 12% of freshwater $Q_0$ discharged into the sea, i.e., from 10.8 to 9.5 m³·day⁻¹·m⁻¹ in each fracture of the parallel set, due to new (or apparent) pumping wells. The maximum distance d of pumping wells from interface (see Figure 3) decreased of about 1000–1500 m, on average, due to new pumping.

Advancement of Saline Front in a Fracture Using N-S

The salt water advancement due to transient density-driven flow in a single fracture was investigated via the solution of the N-S equations. These solutions allowed the estimation of the time period required by saltwater to reach new Ghyben-Herzberg equilibrium due to simulated pumping. The governing N-S equations can be written by considering the salt mass and momentum conservation. In the Lagrangian framework can be written:

$$\frac{D\rho }{Dt}+\frac{\partial \rho u_{\alpha }}{\partial x_{\alpha }}=0$$

(4)

and

$$\frac{D\rho u_{\alpha }}{Dt}=-\frac{\partial {\sigma }_{\alpha \beta }}{\partial x_{\beta }}+\rho g_{\alpha }\text{ and }\alpha ,\beta =x,y,z$$

(5)

where indexes α and β denoted the generic spatial component following Einstein notations, ρ [M·l⁻³] is water density, and u_α is the velocity [l·t⁻¹] of the fluid particle; x_α [l] and x_β [l] are two spatial position coordinates of the fluid particle (Einstein notation); g_α [l·t⁻²] is the component of the gravity vector. The tensor component of the Newtonian stress component σ_αβ [M·t⁻²·l⁻¹], which was applied to each particle subject to the pressure pδ_αβ [M·t⁻²·l⁻¹] can be defined:

$${\sigma }_{\alpha \beta }=p{\delta }_{\alpha \beta }-{\tau }_{\alpha \beta }$$

(6)

where laminar and turbulent stress component τ [M·t⁻²·l⁻¹] of the viscosity can be defined using [21,22]:

$${\tau }_{\alpha \beta }=\left\{\mu \times \left[\left(\frac{\partial u_{\beta }}{\partial x_{\alpha }}+\frac{\partial u_{\alpha }}{\partial x_{\beta }}\right)-\frac{2}{3}{\delta }_{\alpha \beta }\nabla \cdot u\right]\right\}+{\overline{\tau }}_{\alpha \beta }$$

(7)

where δ_αβ [–] is the Kroneker delta; ${\overline{\tau }}_{\alpha \beta }$ is the average sub-grid stress due to the fluctuating variation of velocities around the averaged values during turbulent flows. Similarly, Reynolds stress can be determined by assuming the eddy viscosity theory (Boussinesq’s hypothesis), and by the Smagorinsky constant [23].

Some authors introduce also the Tait’s equation [24], i.e., the Equation of State, in order to complete the unknowns/equations balance.

However, when the first N-S equation is applied to a single computational cell of discretized water flow domain, water flow density remains constant by changing pressure due to water incompressibility and Equation (4) reduces to:

$$\nabla ·u=0$$

(8)

and cannot be used to estimate velocity due to its undetermined form. By using the projection method the velocity u* given by N-S conservation momentum Equation (5) is included in the following (Poisson) [25] equation:

$$\nabla \cdot u^{*}=\mathsf{\Delta }t\nabla \cdot \left(\nabla p\right)$$

(9)

in order to estimate the water pressure. In Equation (9) u* is the intermediate velocity whose value is based on the viscous stress. Equation (8) is then used in a finite difference model (FDM) [16] in order to estimate the pressure of water via a separate numerical calculation. In particular, the conjugate gradient numerical method was applied to solve Equation (8) in order to calculate the water pressure by forcing the flow divergence to zero. After that the water pressure was determined, the correct velocity is obtained by Equation (9). Thus, at the next time step, the N-S momentum conservation Equation (4) is solved again to calculate the new intermediate water velocity. Moreover, in order to account for flow density variations due to the saline front advancement in the FDM code, the velocity u* derived from Equation (5) was updated at each time step according to the new distribution of water densities into the cells of the fracture domain. This water density distribution due to mass salt flow advancement in the fracture was determined at every instant of the flow simulation by solving the salt advection and dispersion equation into the freshwater of the fracture, using a hydrodynamic dispersion coefficient of 10⁻⁹ m²·s⁻¹. In this way in every cell of the domain the divergence of the salt mass flux in the x and z directions was predicted at each simulation time.

The N-S flow equations were solved in a horizontal fracture of Bari aquifer with 5 m of extent and 2b = 3 mm of aperture, which is the maximum fracture aperture derived from pumping-tests carried out on Bari wells. The computational domain x-z was discretized in cells of 0.3 mm × 50 mm of size, i.e., 10 × 100 cells. The transitory density-driven flow of seawater mixed with freshwater was performed assuming that the freshwater fracture, at the inlet section, is subjected to a constant salty water inflow at 35 g·L⁻¹ of total concentration. Moreover, the reduction of 1.3 m³·day⁻¹·m⁻¹ of freshwater discharge Q₀ in each fracture was also applied at the inflow section as the initial conditions during N-S simulations. Flow through the fracture cross-section (2b × 1 m²) was driven by a momentum (per unit volume) equal to 12,375 kg·m⁻²·day⁻¹. This was defined by the product of the reduction of freshwater flowrate along the coast (from 102.9 to 90.5 m·day⁻¹) (due to new pumping of 335 L·s⁻¹) and the freshwater density (12.4 m·day⁻¹ × 1000 kg·m⁻³). Furthermore, at t = t₀ = 0, the salt density distribution in all cells of fracture was also imposed as initial condition.

The N-S code provided the time-dependent distribution of both water pressure and velocity in a fracture, caused by changes in water density due to salt water advancement due to the new pumping.

The FDM code results (see Figure 5 and Figure 6) produced an almost regular flow with a maximum central value of velocity (and Reynolds number) and, subsequently, in both pressure gradient and water density. This occurs because the diffusive/convective salt mass flux divergence decreases from the centerline to the border of the fracture. By using a PC with a single Intel Pentium(R) 32 bits, the FDM code required 1 h and 12 min run-time (CPU) to simulate a period of 792 min. N-S solutions provided information regarding the freshwater/seawater sharp interface positions x_s [l] in the studied fracture at specific simulation time. Using twelve simulation run results the following linear equation was (0.99 of correlation) close-fitting in TableCurve2D:

$$x_s(t)=C_s+B_s\times t$$

(10)

where the position of salt water advancement in the fracture is dependent on time and two constants: C_s = 1.08 × 10⁻¹ m and B_s = 9.36 × 10⁻³ m·min⁻¹ (or 13.3 m·day⁻¹). The results of the N-S code simulations confirmed the strong influence of fracture aperture size on water velocity estimation. For a fracture with 3 cm aperture, the mean water velocity (and Reynolds number) may be 100 times higher than the value determined for a fracture with 3 mm aperture. The quasi linear trend of the saline front advancement given by Equation (10) well matches the temporal trend of the progression of toe positions given by [26] during simulations. The rate of seawater advancement given by the N-S code (i.e., 13.3 m·day⁻¹) is close to the expected value (12.4 m·day⁻¹), which is driven by the prescribed momentum. Using these N-S solutions, maximum elapsed time from the start of pumping to reach the new position of freshwater/saltwater sharp interface in the Bari groundwater is equal to 112.8 day (i.e., 1500 m/(13.3 m·day⁻¹)). Material concerning the advancement of sea movement intrusion is very useful for groundwater management practices of coastal aquifers [26,27,28].

Figure 5. N-S solutions in a fracture (of aperture 2b = 3 mm, 2b stands for the fracture aperture in the manuscript.) of the Bari aquifer. The rate of saline front advancement was used to determine the effects of new pumping of 335 L·s⁻¹ on the freshwater/saltwater transition zone in coastal aquifers.

Figure 5. N-S solutions in a fracture (of aperture 2b = 3 mm, 2b stands for the fracture aperture in the manuscript.) of the Bari aquifer. The rate of saline front advancement was used to determine the effects of new pumping of 335 L·s⁻¹ on the freshwater/saltwater transition zone in coastal aquifers.

Figure 6. N-S solution at time t = 243.03 min from the beginning of the new pumping of 335 L·s⁻¹: horizontal flow velocity magnitude map due to freshwater/saltwater mixing in the fracture, and Reynolds number distribution.

Figure 6. N-S solution at time t = 243.03 min from the beginning of the new pumping of 335 L·s⁻¹: horizontal flow velocity magnitude map due to freshwater/saltwater mixing in the fracture, and Reynolds number distribution.

5. Conversion of Distances from Interface into Groundwater Salinity Data

The conversion of each grid node distance from saline front into groundwater salinity concentration required two data sets, namely: (i) field salinity monitored in wells; and (ii) data from groundwater flow simulations.

5.1. Field Monitored Data

Seventeen wells of the Bari aquifer were monitored for water depth and specific water conductance. These measurements were carried out using multiparameter mini-probes (OTT Hydrolab, Kempetn, Germany) Hydrolab-DS5 and OCEAN SEVEN 315 (IDRONAUT Srl, Brugherio, Italy). The latter one also provided dissolved the water oxygen concentration, which varied from 3.0 mg·L⁻¹ (at the surface) to 0.3 mg·L⁻¹ or less (at 37 m depth) in Bari’s wells, and the water temperature. Each probe was previously calibrated in the laboratory in order to directly provide proper values of water salinity, temperature, pressure and dissolved oxygen vs. water depth. The relative monitored water specific conductance in 14 boreholes is displayed in Figure 7, where C/C_min is the ratio of each measurement to the minimum recorded value at each specific borehole location. Figure 7 presents anomalous salinity trends in the polluted sites that are characterized by high water temperatures.

Figure 7. Relative specific conductance vs. water depth in 14 boreholes of the Bari aquifer positioned at different distance from Ghyben-Herzberg freshwater/sweater (50%) sharp interface position (Winter 2012) and best-fit equation (dashed line).

Figure 7. Relative specific conductance vs. water depth in 14 boreholes of the Bari aquifer positioned at different distance from Ghyben-Herzberg freshwater/sweater (50%) sharp interface position (Winter 2012) and best-fit equation (dashed line).

5.2. Ghyben-Herzberg Data

At first step, an interpolation in TableCurve2D provided the best equation (correlation coefficient of 0.92–0.9) which fits groundwater salt concentrations measured in 25 boreholes at depth of 1.0 m below water table versus the modeled distance d of each borehole from the saline front into groundwater provided by Ghyben-Herzberg theory (see yellow polygon in Figure 3).

$$C_{salt}=C_{s0}+A_s\left[\exp \left(-\frac{d}{D_s}\right)\right]$$

(11)

where, for d ≤ 1500 m, the best fit constants are C_s0 = 1.54 g∙L⁻¹, A_s = 12.02 g∙L⁻¹ and D_s = 592.65 m. At distances d > 1500 m the following best-fit equation:

$$C_{salt}=\left[\exp \left(G_s-I_s\sqrt{d}\right)\right]$$

(11a)

it was preferred to Equation (11), where G_s = 2.54 log(g·L⁻¹) and I_s = 0.04177 log(g·L⁻¹∙m^−0.5) are the best fit constants.

The significance of the above interpolation is that Equations (11) and (11a) enable the calculation of the salinity at depth z = 1.0 below water table at every generic position (x, y) of the computational domain, on the basis of field monitored data and distance d given by groundwater flow model solution.

Thus, a second interpolation of groundwater salinity concentrations measured in wells as a function of the water depth was performed. In particular, the use of TableCurve2D enabled the acquisition of groundwater salt concentrations with depth z at each generic location (x, y), using the best-fit curve displayed (dashed line) in Figure 7 (with correlation coefficient of 0.93).

$$C_{salt}=C_0\times \left[\frac{E_s(z)}{F_s(z)}\right]$$

(12)

where E_s and F_s are two dimensionless interpolating functions with polynomial form:

$$E_s(z)=0.99+2.26\cdot {10}^{-3}z-5.32\cdot {10}^{-4}{z}^2+4.38\cdot {10}^{-6}{z}^3$$

(12a)

$$F_s(z)=1-1.80\cdot {10}^{-2}z-1.07\cdot {10}^{-5}{z}^2+1.03\cdot {10}^{-6}{z}^3$$

(12b)

where z [l] is water depth into borehole and C₀ is minimum salinity value at the specific borehole location. Equation (12) is valid up to a water depth of 140 m. Moreover, the following best-fit function (correlation coefficient of 0.93) could also replace Equation (12):

$$C_{salt}=C_0\times \left\{1.55\cdot \left[1+erf\left(\frac{z-5.6}{\sqrt{5800}}\right)\right]-0.5\right\}$$

(12c)

The significance of this interpolation is that Equation (12) or (12c) enables the calculation of the salinity with depth below water table at every generic position of the computational domain.

Once Equations (11) and (12) were determined, a FORTRAN90 (Visual FORTRAN Professional Edition 5.0.A, Microsoft Developer Studio 97, Microsoft Italia, Peschiera Borromeo (MI), Italy) code was used to perform groundwater salinity estimations. All grid-node positions of computational domain with respect to the freshwater/seawater interface were then converted to groundwater salinity data. A 3D file of 21,070 salt groundwater concentrations was obtained by using 10 grid steps of 10 m in z direction. This file was elaborated using RockWorks15 in order to determine the solid model of actual groundwater salinity, suitable for 3D map visualization. The reliability of this conversion method was quantitatively tested (Table 2) by comparing measured and predicted values at three different depths. The average uncertainty of these estimations is about 20% with respect to the mean salinity (=1.53 g·L⁻¹) (see Table 2) at depth of 1 m below water table and it may be reduced by increasing the number of boreholes for salinity measurements.

Table 2. Comparison between modelled (Mod) values of groundwater salinity (i.e., Rockwork15 input data at t = 0) and measurements (Mea) at three water depths into boreholes of the Bari aquifer. SD stands for the standard deviation between measured and modelled values.

ID Borehole	Actual (Winter 2012) Groundwater Salinity (g·L−1) at Different Water Depths
	1 m			11 m			31 m
	Mod	Mea	±SD	Mod	Mea	±SD	Mod	Mea	±SD
L1-S	0.9	1.54	0.46	1.06	1.65	0.41	1.31
P11	2.76	2.99	0.16	3.14			3.9
P19	3.75	3.71	0.03	4.27			5.29
L2-S	3.12	2.63	0.34	3.55			4.39
L3-S	2.11	3.40	0.91	2.4			2.97
L4	0.98	1.05	0.05	1.1	1.72	0.44	1.38	2.77	0.98
L8-S	0.5	0.71	0.15	0.52			0.65
L7-S	0.5	0.58	0.06	0.6			0.74
L6-S	0.5	0.67	0.12	0.52	0.78	0.19	0.65
L9	2.02	1.54	0.34	2.3			2.85
L14-S	0.93	0.59	0.24	1.06	0.56	0.35	1.31	1.03	0.20
L15-S	0.93	0.86	0.05	1.06	0.95	0.08	1.31	1.75	0.31
L5-S	0.5	0.63	0.09	0.52	0.68	0.12	0.65	1.50	0.60
L12	2.02	0.56	1.03	2.3	2.26	0.03	2.85
Mean		1.53	±0.32		1.23	±0.23		1.76	±0.52

Table 2. Comparison between modelled (Mod) values of groundwater salinity (i.e., Rockwork15 input data at t = 0) and measurements (Mea) at three water depths into boreholes of the Bari aquifer. SD stands for the standard deviation between measured and modelled values.

ID Borehole	Actual (Winter 2012) Groundwater Salinity (g·L−1) at Different Water Depths
	1 m			11 m			31 m
	Mod	Mea	±SD	Mod	Mea	±SD	Mod	Mea	±SD
L1-S	0.9	1.54	0.46	1.06	1.65	0.41	1.31
P11	2.76	2.99	0.16	3.14			3.9
P19	3.75	3.71	0.03	4.27			5.29
L2-S	3.12	2.63	0.34	3.55			4.39
L3-S	2.11	3.40	0.91	2.4			2.97
L4	0.98	1.05	0.05	1.1	1.72	0.44	1.38	2.77	0.98
L8-S	0.5	0.71	0.15	0.52			0.65
L7-S	0.5	0.58	0.06	0.6			0.74
L6-S	0.5	0.67	0.12	0.52	0.78	0.19	0.65
L9	2.02	1.54	0.34	2.3			2.85
L14-S	0.93	0.59	0.24	1.06	0.56	0.35	1.31	1.03	0.20
L15-S	0.93	0.86	0.05	1.06	0.95	0.08	1.31	1.75	0.31
L5-S	0.5	0.63	0.09	0.52	0.68	0.12	0.65	1.50	0.60
L12	2.02	0.56	1.03	2.3	2.26	0.03	2.85
Mean		1.53	±0.32		1.23	±0.23		1.76	±0.52

5.3. Data from Solutions of N-S Equations

At third step, the same FORTRAN90 code previously applied was implemented applying Equation (10) to update the salinity concentrations in each grid node of computational domain at three different simulation times of 30, 75 and 112.8 days. The saline front displacement during time is given by N-S solutions, by determining the new grid node positions with respect to the salinity front. In particular, as input for this calculation, 49 distances of the top border of domain from the salinity front were provided at each simulation time. A control on the maximum distance calculated was imposed into the FORTRAN90 code in order to avoid that the advancement of the saline front may exceed the distance provided by flow model (i.e., blue polygon) plotted in Figure 3. The significance of this calculation is that it enables the estimation of groundwater salinity concentrations according to the advancement of salinity front in groundwater fractures during a new pumping obtained by solutions of the Navier-Stokes equations. The results provide 3D solid models in RockWorks15. 3D maps of Figure 8 show actual salinity of Bari coastal aquifer and its progression as the consequence of the constant new simulated pumping of 335 L·s⁻¹.

Figure 8. (a) 3D map of groundwater salinity during winter 2012 (t = 0 day) and its comparison with groundwater salinity maps after (b) 30; (c) 75 and (d) 112.8 days of a simulated continuous pumping of 335 L∙s⁻¹.

Figure 8. (a) 3D map of groundwater salinity during winter 2012 (t = 0 day) and its comparison with groundwater salinity maps after (b) 30; (c) 75 and (d) 112.8 days of a simulated continuous pumping of 335 L∙s⁻¹.

6. Discussion and Conclusions

A new method to infer saline front advancement in a fractured coastal aquifer was applied to Bari groundwater. Results visualize the 3D maps of groundwater salinity concentrations, at different simulation times, due to the impact of an apparent pumping of 335 L·s⁻¹. The 3D maps are useful tools for coastal groundwater management and show how over-abstractions may affect groundwater salinity changes. The same method can be applied to others aquifers with discrete fractures and it requires two data sets: (i) field salinity data in boreholes (i.e., logs); and (ii) results of groundwater flow model simulation to represent actual 3D maps of aquifer salinity. Moreover, solutions of the N-S equations in a fracture can provide the rate (on average) of saline front advancement in fractures during time. This result allows updates of actual salt concentrations, by defining 3D salinity maps at specific simulation times during a new pumping. Although mapping based on the sharp interface approach is not representative of real conditions of the salt mass transport in fractures, the salinity front positions at the initial and final stage of the simulated pumping stress condition in this method, take into account for the real flow conditions in the fractured aquifer. Moreover, by including approximations due to the interpolation stages (correlation coefficients >0.92–0.93) authors have estimated the 20% of uncertainty of results at the Bari aquifer. This uncertainty established an acceptable distance of results from reality.

In the present work groundwater flow modeling was based on aquifer transmissivity values derived from pumping and tracer (chlorophyllin) injection tests carried out during winter 2012. At specific sites, a minor number (<30) of monitoring wells may increase limitations of results, and the uncertainty of the presented method.

The presented method is an alternative to the application of the 3D unsteady density-driven flow models at a regional scale, which may present computational limitations in discrete fractured aquifers due to complexity of input data required to describe preferential flow pathways and the unknown position (and aperture) of the fractures and their intersections, especially at a regional scale. Moreover, conventional meshed models based on finite elements, finite differences or boundary elements, are inappropriate for flow and transport simulations in fractured aquifers, i.e., not continuous porous media. This is due to impossibility to define a representative elementary volume (REV) of appropriate size in a fractured (or very heterogeneous) aquifer.

Data regarding the seawater front advancement in the Bari aquifer are highly suitable for groundwater management of coastal aquifers. Results suggest an elapsed time of 112.8 days to achieve the new Ghyben-Herzberg interface equilibrium position, starting from the beginning of the new pumping of 335 L·s⁻¹.