Saltwater intrusion estimation in a karstified coastal system using density-dependent modelling and comparison with the sharp-interface approach

Abstract

Saltwater intrusion is a naturally occurring phenomenon that is exacerbated significantly by excessive groundwater exploitation in coastal aquifers. In order to determine the extent of saltwater intrusion in a karstified coastal aquifer in Crete, Greece, a three-dimensional, density-dependent groundwater flow and transport model was developed and compared to the more traditional sharp-interface approach. The karstified medium was modelled using a combination of the equivalent porous medium approach (for lower-order fractures) and a discrete fracture approach (for the main fractures/faults). The model takes into consideration the geomorphologic characteristics of the karstic system, such as the depth and orientation of the fault network, and the diffusion phenomena associated with the variable densities of freshwater and saltwater—parameters that create a complex system, inducing uncertainty in the model. The model results showed that the orientation of the fractures, the pumping activity and the fluid density effects drive the seawater intrusion front asymmetrically inland.

Editor Z.W. Kundzewicz

Citation Dokou, Z. and Karatzas, G.P., 2012. Saltwater intrusion estimation in a karstified coastal system using density-dependent modelling and comparison with the sharp-interface approach. Hydrological Sciences Journal, 57 (5), 985–999.

Résumé

L'intrusion d'eau salée est un phénomène naturel, qui peut être accru de manière significative par l'exploitation excessive des eaux souterraines dans les aquifères côtiers. Afin de dèterminer l'ampleur de l'intrusion d'eau salée dans un aquifère karstique côtier en Crète (Grèce), un modèle en trois dimensions d'écoulement et de transport souterrain dépendant de la densité a été mis au point et comparé à l'approche plus traditionnelle d'interface abrupte. Le milieu karstique a été modélisé en utilisant une combinaison de l'approche de milieu poreux équivalent (pour les fractures d'ordre inférieur) et une approche discréte de fracture (pour les fractures et failles principales). Le modéle prend en compte les caractéristiques géomorphologiques du systéme karstique, telles que la profondeur et l'orientation du réseau de failles, et le phénomène de diffusion associé à la densité variable entre eau douce et eau salée—paramétres qui créent un systéme complexe, induisant de l'incertitude dans le modèle. Les résultats du modèle ont montré que l'orientation des fractures, l'activité de pompage et les effets de la densité du fluide conduisent à un front asymétrique d'intrusion d'eau de mer à l'intérieur des terres.

Keywords:

• groundwater flow and transport modelling
• saltwater intrusion
• density-dependent flow
• sharp-interface approach
• Ghyben-Herzberg relationship
• karstic systems

Mots clefs:

• modélisation de l'écoulement souterrain et du transport
• intrusion d'eau salée
• écoulement dépendant de la densité
• approche d'interface abrupte
• relation de Ghyben-Herzberg
• systèmes karstiques

INTRODUCTION

The intrusion of saltwater into coastal aquifers is a natural process that is often augmented by pumping activity disturbing the equilibrium of the saltwater–freshwater system. Continuous and unrestrained pumping activity results in the lowering of the water table in these zones and also reduces the natural subsurface flow of freshwater from the mainland towards the shoreline. The sudden water table level decline affects not only water quality but also the hydraulic connection between the surface and subsurface water resources. The consequences of the saltwater intrusion phenomenon for the environment, ecology and economy of a coastal region are very serious (Datta et al. 2009).

Coastal karstified aquifers are at greater risk of contamination from saltwater intrusion processes. The existence of irregular fissures and solution openings enables saltwater to enter the aquifer using different flowpaths to those in homogeneous aquifers (Bear et al. 1999). Thus, the development of accurate groundwater and contaminant transport models to simulate and predict the location and width of the saltwater–freshwater interface in karstified aquifers, using equations that represent the physical processes governing the phenomenon as closely as possible, is required.

There are two main approaches to describe the evolution of the saltwater front in a coastal aquifer mathematically. The most simplified, but nevertheless widely-used approach is based on the sharp-interface approximation and the Ghyben-Herzberg relationship developed independently by Ghyben (1888) and Herzberg (1901) at about the same time. According to this methodology, the physical system consists of two completely immiscible fluids with different constant densities (fresh and saline water) that form a mixing zone that is limited to an interface with a small finite width (Reilly and Goodman 1985). This approach has been applied and extended by many researchers such as: Mahesha and Nagaraja (1996), Dagan and Zeitoun (1998), Naji et al. (1998), Guvanasen et al. (2000), Brunke and Schelkes (2001), Mantoglou (2003), Mantoglou et al. (2004), Ababou and Al-Bitar (2004), Karterakis et al. (2007), Koukadaki et al. (2007), Papadopoulou et al. (2010a) and Papadopoulou et al. (2010b). The interested reader is directed to Reilly and Goodman (1985) and Diersch and Kolditz (2002) for an extensive review of works utilizing the sharp-interface approach.

The second approach, termed the density-dependent flow approach, assumes that the freshwater–saltwater system consists of a miscible fluid transporting a solute that affects the density and viscosity of the fluid. In this case the width of the mixing zone between the fresh and saline water is significant. In the early 1990s, Schincariol and Schwartz (1990) developed experiments in order to investigate the density-dependent flow in heterogeneous porous media, proving the flow sensitivity to the aquifer's hydraulic conductivity. The first attempt to model density-dependent flow in saltwater intrusion problems was made by Henry (1964). Later, various formulations and numerical methods were applied to the solution of the density-driven saltwater intrusion problem (Pinder and Copper 1970, Huyakorn and Taylor 1976, Frind 1982, Huyakorn et al. 1987, Diersch 1988, Putti and Paniconi 1995, Benson et al. 1998, Gambolati et al. 1999, Abarca et al. 2007, Servan-Camas and Tsai 2009, among others), whose validity was tested using benchmark problems.

The application of density-dependent models in real world applications has recently increased, especially regarding three-dimensional (3D) models, due to the advances in computational speed enabling the use of finer grids and complex model geometries (Abarca et al. 2007). Such field applications have been implemented at coastal sites all around the world. The earliest attempts were 2D models applied in the Biscayne aquifer near Miami, Florida, USA (Segol and Pinder 1976), in Oahu, Hawai (Voss and Souza 1987) and in the Nile Delta (Sherif et al. 1988). More recently, 3D models were developed and applied in Florida (Panday et al. 1993), in Huangheying, Longkou, China (Xue et al. 1995), in La Vinuela reservoir system in Malaga, Spain (Garcia-Arostegui 1998), in the Korba coastal plain, Tunisia (Paniconi et al. 2001, Kerrou et al. 2010), in The Netherlands (Essink 2001), in Morocco (Aharnouch and Larabi 2004), in the Kiti aquifer, southern Cyprus (Milnes and Renard 2004), in Ravena, Italy (Giambastiani et al. 2007) and in Santorini, Greece (Kopsiaftis et al. 2009, Kourakos and Mantoglou 2009).

Some researchers have investigated the tidal influence on the flow and transport processes in coastal areas (Li et al. 1997, 1999, 2009). Indeed in coastal aquifers where a significant rise and fall of the water table occurs due to tidal effects, both the inland hydraulic gradient and the salt transport into the aquifer can be affected (Li et al. 1999).

The high uncertainty associated with the presence and the geomorphology of fractures and faults within a karstified medium has tremendous effects on the modelling of subsurface flow and contaminant transport, as it results in complex flowpath geometry (Papadopoulou et al. 2010c). The fundamental problems that modellers face when trying to simulate a karstified physical system are the system's hydraulic characteristics (extreme values of aquifer hydraulic conductivity and porosity) and its distribution (network of permeable channels and conduits embedded in low permeability fractured rocks). The geometry of this network may affect the hydraulic behaviour and the contaminant mass transport, rendering the region particularly vulnerable to pollution. Park and Zhan (2003), in their analysis of groundwater hydraulic head in the vicinity of a horizontal well in fractured and porous aquifers, showed that the water flow from the matrix to a fracture cannot be ignored even if the hydraulic conductivity of the matrix is small compared to that of the fracture.

Three main approaches have been widely used in the past for modelling karstified aquifers. The first views the fractured system as an equivalent porous medium. In this approach, the spatial variations in the hydrogeological properties of the rock mass are averaged over a representative elementary volume (REV) in order to define their values at the macroscopic level. Many researchers have employed the equivalent porous medium approach over the past decades (Long et al. 1982, Berkowitz et al. 1988, Schwartz and Smith 1988, Angelini and Dragoni 1997, Larocque et al. 1999, Scanlon et al. 2003, Zhang et al. 2010, Dokou and Pinder 2011) for simulating saturated groundwater flow and transport in fractured materials.

The second approach, termed the dual porosity model, was introduced by Barenblatt et al. (1960) and represents the fractured medium using two interacting continua: the primary porosity blocks that exhibit low permeability and high storage capacity, and the secondary porosity fractures that exhibit high permeability and negligible storage capacity. An exchange of fluid occurs between the two continua represented by a leakage term (Therrien and Sudicky 1996). Bibby (1981), Huyakorn et al. (1993), Sudicky (1990), Gerke and van Genuchten (1993), Samardzioska and Popov (2005) and Zyvolovski et al. (2008), among others, have studied the applicability of this approach. A similar model, called the coupled continuum pipe-flow model, has been used by researchers in order to describe the flow and solute transport in fractured aquifers (Liedl et al. 2003, Faulkner et al. 2009, Hu 2010). In this model the flow exchange between the matrix and the fractures is controlled by an exchange rate coefficient that depends on the hydraulic head and conductivity differences in the two systems and the fracture geometry (Hu 2010).

The third approach is a discrete fracture approach that utilizes one- and two-dimensional elements to represent faults, fractures or karst networks, while neglecting the hydraulic properties of the matrix blocks. It is typically applied to fractured systems with low permeability and low connectivity between fractures (Jaquet et al. 2004). In this approach, the fracture characteristics such as their location, orientation, length and aperture, are usually created via statistical methods.

The applicability of these approaches depends on a variety of factors, such as the objective of the study, the degree of karstification, the scale of the model and the availability of different kinds of data (Scanlon et al. 2003). Some researchers have developed numerical models of coastal karstified aquifers combining two of the aforementioned methods: the equivalent porous medium and the discrete fracture approach (Clemens et al. 1996, Kiraly 1998, Graf and Therrien 2005, Quinn et al. 2006, Papadopoulou et al. 2010b).

The work presented in this paper focuses on the development of a saltwater intrusion model of a karstified coastal area located in the northern part of Crete, Greece, using both the density-dependent approach and a combination of the sharp-interface approximation and the Ghyben-Herzberg relationship, and comparing their results. For both methods, the partial differential equations are solved using the same finite element numerical model, FEFLOW, configured each time to solve a different problem. The karstified medium is modelled using a combination of the equivalent porous medium approach, to represent the lower-order fractures or channels, and a discrete fracture approach, to represent the main fractures or karst networks.

METHODOLOGY

For the purpose of modelling the karstified coastal aquifer of interest, the finite element subsurface flow and transport simulation system, FEFLOW, was employed. Two different modelling approaches to the saltwater intrusion phenomenon were used and compared: the density-dependent and the sharp-interface approaches.

Density-dependent approach

This approach assumes that the two fluids, freshwater and saline water, are miscible and that the width of the mixing zone between them is significant. The main assumptions are: (a) the porous medium is considered rigid and saturated, (b) the incompressibility condition and Boussinesq approximation hold, (c) fluid viscosity depends on concentration, (d) fluid motion can be adequately described by Darcy's law, and (e) dispersion is a property of the porous medium alone that follows the Fickian dispersion law. Furthermore, the Einstein summation convention for repeated indices is used (Diersch 1988). Under these assumptions, the saltwater intrusion phenomenon can be described by the governing density-dependent flow and transport equations:

Continuity equation:

(1)

(1)

Density difference ratio:

(2)

Fluid density equation:

(3)

Transport equation:

(4)

(4)

Darcy's law:

(5)

(5)

where q_i is Darcy's velocity; S_s is specific storage; h is the freshwater hydraulic head; Q _ρ is the fluid source/sink; Q _EB is a term related to the Boussinesq approximation; α is the density difference ratio; ρ_s is saltwater density; ρ_f is freshwater density; C is salt concentration; C_s is maximum salt concentration; φ is porosity; D_ij is the hydrodynamic dispersion tensor; Q_c is the contaminant mass source/sink; K_ij is the hydraulic conductivity tensor; f _μ is a viscosity constitutive relationship; and η _j is a gravitational unit vector. For a more detailed description of the density-dependent flow and transport equations, please see Diersch et al. (1988).

The model requires freshwater hydraulic heads as input; thus, an equivalent hydraulic head needs to be calculated when applying boundary conditions (BCs) (constant head boundary conditions along the sea boundary) and when comparing model results to field data. The transformation is performed using the following equation, assuming that saltwater density varies linearly with depth:

(6)

where h _eq is equivalent freshwater head; h is measured head; α is the density difference ratio; and z is elevation.

In order to calculate the density difference ratio for the field data, the density of the freshwater–saltwater mixture in each well was needed. If a chloride concentration measurement was taken in that same well, this measurement was used in order to calculate the density and then the equivalent freshwater head at the particular well. If not, concentration measurements from nearby wells were used instead.

To accommodate free surface conditions, the model uses a general moving mesh strategy (BASD—Best Adaptation to Stratigraphic Data), which is capable of computing movable 3D finite element meshes even under general stratigraphic heterogeneous conditions (Diersch 2002a).

In typical saltwater intrusion problems, the freshwater–saltwater boundary conditions play an important role in the phenomenon and should be handled properly. The concentrations on the upper sea boundary change dynamically, given that, when water is entering the domain, it should have a freshwater concentration, while, when it is leaving the domain, it has an unknown concentration that needs to be calculated by the model. This is handled by using a boundary condition of constant concentration, C = C^R and, at the same time, imposing an additional constraint of zero mass flux, , on the same boundary (Diersch and Kolditz 2002). This constraint results in a decrease in salt concentration in the upper parts of the aquifer due to the zero dispersive flux along the out-flowing part of the sea boundary (Kopsiaftis et al. 2009).

Sometimes, high chloride concentrations near coastal zones may not be the result of saltwater intrusion, but, instead, originate from anthropogenic activities. In order to use the correct mass boundary conditions (sources), it is important to be able to distinguish between the two. For this purpose, in this work, the Mg/Ca ratio was used. This ratio is an indicator of the origin of chloride contamination based on the fact that, in saltwater, Mg²⁺ is found in excess of Ca²⁺, while Ca²⁺ is predominant in freshwater when it is not polluted by other anthropogenic activities (Bear et al. 1999).

Sharp-interface approach

According to this approach, the mixing zone between the two immiscible fluids with different constant densities (freshwater and saline water) is limited to an interface with a small finite width. The sharp-interface approach can be used in problems when the transition zone is relatively narrow, compared to the scale of the problem (Mantoglou et al. 2004). The location of the sharp interface between the two fluids is determined by the difference between the hydraulic heads of the saline and fresh water and the volume of freshwater flowing towards the shoreline from inland. Based on the sharp-interface theory and under steady-state conditions, the physical system is modelled using the standard groundwater flow equation:

(7)

and the seawater intrusion front can be approximated hydraulically using the Ghyben-Herzberg relationship:

(8)

where ξ is the interface depth below the sea level; h _f is the hydraulic head of the freshwater above the sea level; ρ_f is the density of freshwater; and ρ_s is density of saline water.

Assuming that freshwater density is 1000 kg/m³ and saline water density is 1025 kg/m³, equation (8)(8) becomes:

(9)

Based on equation (9)(9), for each metre of freshwater hydraulic head (h _f) above the mean sea level, the interface between the saline water and freshwater is balanced at approximately 40 m below the mean sea level.

The boundary conditions used in such a model are a constant head boundary of zero metres at the coast (first type BC), a zero flux at the no flow boundaries (second type BC) and a constant or time-dependent influx or outflux at the flow boundaries (second type BC). The pumping or injection activities at the wells are also treated as a second type BC with the withdrawal or injection rates equal to the outflux or influx BC value.

Under transient conditions, the behaviour of the interface is controlled by both the freshwater and saltwater dynamics, because significant amounts of saltwater must enter or exit the aquifer following the interface movement (Mantoglou et al. 2004).

Discrete features

The FEFLOW system offers an integrated approach for modelling karstified aquifers that combines a 3D equivalent porous medium with discrete features. This approach provides large flexibility in modelling complex systems (Diersch 2002b). The equivalent porous medium is used to represent the matrix structure of the aquifer system, and the discrete features are utilized in order to represent either channels, rivers or wells (1D features) and runoff processes, fractures or faults (2D features), as shown schematically in Fig. 1.

Fig. 1 Schematic representation of the combined equivalent porosity and discrete feature approaches for modelling groundwater flow and transport.

In the work presented herein, the main fractures and faults found on the field site are represented by 2D discrete feature elements that follow the Hagen-Poiseuille law:

(10)

(11)

and the discharge of a single fracture is expressed as:

(12)

where v is average velocity in the fracture; b is fracture aperture; Q is discharge; and the other parameters are as defined previously.

CASE STUDY

Site description

The area of interest is an industrial zone (BIPE) located on the coastal aquifer east of the city of Heraklion, Crete. The coastline forms the northern boundary of the area, which extends approximately 4500 m to the south. Saltwater intrusion phenomena, exacerbated by intensive pumping due to the extensive water needs of local industry, in addition to the dumping of industrial waste directly into the aquifer, have deteriorated the quality of groundwater in the area. A topographic map of the area of interest is shown in Fig. 2.

Fig. 2 Topographic map of the area of interest.

A preliminary work by Koukadaki et al. (2007), focusing on the estimation of hydraulic conductivity obtained from field measurements of electrical resistivity, provided the geological mapping of the area, where fractures, karstified and non-karstified geological formations occur (Fig. 3). The main formations present in the area are:

Fig. 3 Geology of the study area, location of wells and industrial zone.

a.	bioclastic limestones of Neogene age (St Barbara formation): mainly marly limestones with some interlayered gypsum beds, which have an estimated hydraulic conductivity of 5 m/d;
b.	carbonate formations of the Tripolis zone: this unit includes limestones, dolomitic limestones and dolomites occupying the tectonic base of the area with an estimated hydraulic conductivity of 10.5 m/d; and
c.	alluvial deposits and sea sands of Quaternary age: present along the eastern and western coastline with an estimated hydraulic conductivity of 200 m/d (Koukadaki et al. 2007).

Model development

The horizontal discretization of the unconfined aquifer was implemented using a triangular finite element mesh consisting of 101 790 nodes and 180 018 elements. The vertical discretization of the model area was based on vertical cross-sections created using information provided by boring logs. One of them (A–A′) is shown in Fig. 4. The conceptual model comprises an unconfined aquifer that is about 80 m deep and was discretized in nine numerical layers. The elevation of the model's top layer was initially variable, following the topography of the area, and then it was adopted depending on the model results using a moving mesh strategy that accommodates free surface conditions. The elevation of the second layer was set to 1 m below sea level, while the rest of the model layers have a thickness of 10 m (Fig. 5). The industrial zone (enclosed by the circle), the location of 16 private and municipal wells situated inside the model area, the faults encountered in the area (represented by the dotted lines) and the location of the geological cross-section A–A′ are shown in Fig. 3.

Fig. 4 Geological cross-section.

Fig. 5 Vertical discretization of the conceptual model.

Boundary conditions for flow

A first-type flow boundary condition was set along the coastline. For the first layer, the hydraulic head was set to zero, and, for the remaining layers, it was set as equivalent to freshwater, assuming a saltwater density of 1025 kg/m³; thus the difference ratio, α, was equal to 0.025 for this application.

Pumping or injection wells in the area of interest were defined as second-type boundary conditions and their screens were set at various depths, from –1 to –10 m, and –10 to –20 m (second and third layer, respectively, depending on the information on their boring logs). Some of the pumping or injection rates were held constant during the simulation and some varied with time, depending on field conditions. In addition, a time-dependent influx of 0.052 m/d during the winter and 0.048 m/d during the summer months was defined on the western model boundary. The above values were a result of the calibration process discussed in the next section. An average monthly infiltration rate (30% of the average monthly rainfall rate) that varied from a maximum value of 0.872 mm/d in January to a minimum value of 0.0038 mm/d in August, based on meteorological data of the area, was also included as an inflow parameter on the top layer of the model.

Boundary conditions for mass

A first-type boundary condition for mass transport was imposed along the coastline with a constant value of 18 980 mg/L of chloride ion concentration. This value is based on field measurements on saltwater taken directly from the sea in the northern part of the area. This boundary condition was coupled with an additional constraint as discussed in the previous section. A constant mass boundary condition (with slightly lower salt concentration of 17 500 mg/L) and a corresponding constraint were also imposed along the eastern fault. The significance of this fault is great because it originates from the coast and runs along the border of the industrial area (where extremely high chloride concentrations were observed) serving as a “contaminant highway” that transports saltwater rapidly over a long distance.

Discrete features

During the initial development of the model a total of seven main fractures were included, all located in a depth that corresponds to the second model layer. The fractures were represented by triangular 2D discrete feature elements of 0.5 m thickness and 0.5 m aperture, values based on site observations conducted by an expert geologist. Their characteristics (except their length and orientation) were kept the same because of lack of data. For the fractures, diffusivity was set to 1 m²/s, the longitudinal and transverse dispersivities were 50 and 5 m, respectively, and the density ratio (where applicable) was set to 0.025.

The main fracture of the area is located along the eastern boundary of the model domain and is considered to be the main avenue of rapid transport of saltwater contamination from the coast to the industrial area. An equivalent hydraulic head boundary condition was set along this fracture that matches the BC of the coastal boundary at the second layer (0.025 m) that was calculated using equation (6)(6).

Additional parameters

The model was run in transient mode for a period of 30 years (since saltwater intrusion is a slow process), accounting for seasonal variations in pumping activities, infiltration and lateral influx rates. An automatic time step control was used in order to overcome numerical difficulties in solving the nonlinear equations. The porosity was estimated to be 0.2 for all the karstic rocks (Tripolis and marly limestone), and 0.25 for the alluvial deposits. The longitudinal and transverse dispersivity were set to 50 and 5 m, respectively, based on model calibration. The spatial distribution of the hydraulic conductivity obtained from Koukadaki et al. (2007) and used in the model is presented in Fig. 6. The hydraulic conductivity values were not held constant throughout the formations, but an inverse distance interpolation technique was used in order to obtain a smoother transition between them.

Fig. 6 Spatial distribution of hydraulic conductivity in the top model layer.

RESULTS

Model calibration

The calibration of the flow field was based on hydraulic head data measured in nine wells in November 2004. The parameters that were varied during the calibration procedure included some of the flow rates in pumping and injection wells located in the area, the lateral influx rate on the western model boundary and the constant head boundary condition of the equivalent freshwater head along part of the eastern boundary. In the field, there is rarely reliable information regarding the actual pumping or injection rates and thus changing their values in order to match the observation points is a well-justified technique. Nevertheless, the pumping or injection rate information obtained from reports was left unchanged wherever possible. The final pumping/injection rates used in the model are presented in Table 1. In addition, the time-dependent influx defined on the western model boundary was also a result of the calibration process. Finally, the constant head boundary condition of the equivalent freshwater head along the eastern boundary was slightly increased along a small segment adjacent to the industrial area (from 0.025 to 0.06 in the second layer, and from 0.25 to 0.6 in the third) to better match field observations.

Table 1  Calibrated pumping/injection rates for the wells in the area of interest

Well	Rate (m^3/d)	Model	layer	Type	Duration
G2	28.57	Field info	L3	Pumping	Yearly
G3	900	Calibrated	L4	Injection	Yearly
G5	320	Field info	L2	Pumping	Yearly
G7	21.5	Field info	L3	Pumping	Winter
G12	28.57	Field info	L3	Pumping	Yearly
G13	1000	Calibrated	L3	Injection	Yearly
G14	350	Calibrated	L3	Pumping	Summer
G23	125	Calibrated	L2	Pumping	Summer
E1	225	Calibrated	L5	Pumping	Yearly
P2	315	Calibrated	L2	Pumping	Yearly

The flow calibration results show a good agreement between measured and modelled values in most cases (Table 2, Fig. 7) with the exception of wells G3 and G13 that exhibit extremely high hydraulic head values indicating that the measurements were either erroneous or were taken under injection conditions. Alternatively, they could be a result of the complex hydrogeological characteristics of the study area, thus they cannot serve as reliable data. Regardless, an effort was made to match those high measurements up to a point, resulting in significantly increased and unrealistic injection rates for those two wells. The R ² value for the comparison between data and model results was 0.875, indicating a good fit. The flow field contours of the model area in terms of hydraulic head are presented in Fig. 8.

Table 2  Flow calibration results

Well	Field data in saltwater head (m)	Field data converted to freshwater head (m)	Model results in freshwater head (m)
G1	2.13	2.15	1.44
G3	19.35	20.0	7.76
G4	0.5	1.13	1.16
G8	0.5	1.01	1.00
G13	35.1	35.68	11.10
G23	−1.84	−1.83	−1.90
E1	−0.1	0.5	0.54
E2	3.1	3.14	3.10
P2	−6.18	−6.05	−6.1

Fig. 7 Comparison of measured and simulated hydraulic head values using the density-dependent approach.

Fig. 8 Hydraulic head distribution for model layer 3 using the density-dependent approach and location of all fractures used in the model.

The contaminant mass model calibration was based on the following concept: The main eastern fault serves as a means of rapid transport of saltwater contamination to the area close to the industrial site. Then, through three smaller fractures with an east to west orientation the contamination continues towards the industrial area. These new fractures were included in the model in addition to the already-known original fractures shown in Fig. 3, and are represented by the dotted lines in Fig. 8. The orientation of these additional fractures was towards the wells exhibiting the highest concentrations of chlorides. This was the only plausible way to explain the extremely high chloride concentrations (reaching as high as 17 300 mg/L in well G3). Their characteristics were the same as the original fractures. In order to match the field observations, the mass source rates released from these additional fractures were set different for each one and varied from 250 to 480 mg/L per day. The mass transport model calibration also included a trial and error procedure for the longitudinal and transverse dispersivities, having as a starting point the values of 20 and 2, respectively, and increasing them in order to achieve a better fit between modelled and measured chloride concentration values at the observation wells.

During the chemical analysis of the field data, a measurement of the Mg/Ca ratio was made for all samples. This ratio is an indicator of the origin of chloride contamination. The Mg/Ca ratio of the samples collected in the area of interest indicated that well G14 is solely contaminated by anthropogenic activities (probably by dumping of waste from the nearby industries). This is a very reasonable conclusion taking into consideration the fact the well G14 is located far from the eastern fault that acts as a source of saltwater pollution. In wells G1, G2 and G12, the origin of contamination is probably a mixture of anthropogenic and saltwater intrusion processes. For the rest of the wells the Mg/Ca ratio suggests that the contamination is mostly attributed to saltwater intrusion phenomena. Based on the above, a contaminant source modelled as a second type boundary condition (flux) was applied in well G14 to simulate the injection of industrial wastes.

The chloride concentration data used for the calibration process were measured in nine wells in November 2004 and in four wells in April 2005 (Table 3, Fig. 9). The model values for the wells measured in November correspond to results from the end of the dry season, while those measured in April correspond to model results obtained five months later (at the end of the wet season). The R ² value for the comparison between data and model results was 0.966, indicating a very good fit. Indeed, a good match between observed and measured values was achieved in most cases, except in wells G1, G2 and G12. These wells exhibit higher model values than measured because they are located very close to the extreme value of 16 480 mg/L (well G3) and thus they are affected by it. The saltwater intrusion extent in terms of chloride mass concentration is shown in Fig. 10. It is observed that on the northern border the saltwater intrusion phenomenon is more prevalent on the western side, where fractures exist, and the hydraulic conductivity is higher due to the alluvial deposits, reaching 545 m inland, compared for example with the middle part where the intrusion extends 365 m inland, when considering the 750 mg/L as a limit (the extent is greater if we consider a lower isochlor, as was done for the comparison of the two modelling methods in the next section). In the industrial area, the three smaller fractures added during calibration transport the saltwater for distances as large as 1.15 km from the eastern fault.

Fig. 9 Comparison of measured and simulated chloride concentration values.

Fig. 10 Saltwater intrusion extent based on chloride concentration (mg/L) for layer 3.

Table 3  Chloride concentration field measurements and model results

Well	Date collected	Cl^- Field data	Cl^- Model results
		(mg/L)	(mg/L)
G1	April 2005	670	3 596
G2	April 2005	1070	1459
G3	April 2005	17300	16480
G4	November 2004	15763	15766
G5	November 2004	9360	9430
G7	November 2004	15000	14903
G8	November 2004	13950	13727
G12	November 2004	1110	4641
G13	November 2004	7800	7905
G14	April 2005	3880	3732
G21	November 2004	282	2840
G22	November 2004	313	362
G23	November 2004	532	547

Comparison with sharp-interface approach

As mentioned before, in the sharp-interface approach, under steady-state conditions, the seawater intrusion front is approximated hydraulically using the Ghyben-Herzberg relationship (equation (9)(9)). In this case study, the depth of the modelled aquifer is 80 m below the sea level (ξ = 80 m) thus, according to the above equation, the hydraulic head of the freshwater at the interface should be h _f = 2 m. This means that, as long as the equivalent hydraulic head of freshwater remains 2 m above sea level, the movement of saltwater into the coastal aquifer can be avoided.

The sharp-interface approach was also applied using the same numerical model turning the density dependence option off and considering only the groundwater flow equations. A fine tuning of the flow boundary conditions was needed to achieve a better agreement between modelled hydraulic head values and measured equivalent head data. Specifically, a constant head boundary condition of 0 m was set along the coastline for all layers. The pumping/injection schedules of the wells in the area and also the infiltration rates were kept the same as in the density-dependent method. On the western boundary, the time-dependent influx was slightly increased reaching 0.06 m/d during the winter and 0.055 m/d during the summer months. The sharp-interface model calibration achieved an R ² value of 0.8754, indicating a good fit (Fig. 11).

Fig. 11 Comparison of measured and simulated hydraulic head values using the sharp-interface approach.

In order to examine the closeness of the extent of seawater intrusion for the two approaches, the 2 m contour line was compared with the 250 mg/L chloride concentration contour (1.3% isochlor) which is the parametric value specified by the European Community Directive 98/83/EC on the quality of water intended for human consumption.

A visual comparison between the density-dependent and sharp-interface approaches is provided in Fig. 12. Assuming that the 2 m isoline and the 250 mg/L isochlor should coincide, it is evident that the results from the two approaches differ significantly. The extent of saltwater intrusion in many regions is vastly over-estimated when utilizing the sharp-interface approach. Especially on the northern part of the model, the extent of saltwater intrusion according to the density-dependent approach is 450 m, while for the sharp-interface method it reaches 1450 m inland. Regarding the industrial area (which is the area of utmost importance in this application because of the very high chloride concentrations detected in many wells), the two approaches show very similar results with the density-dependent extent estimation being a little larger (about 150 m).

Fig. 12 Chloride concentration contour (isochlor) of 250 mg/L (density-dependent approach) and 2 m head contour (saltwater intrusion extent according to the sharp-interface approach).

Since this method is based on a hydraulic balance approach, its accuracy depends on a good estimate of the aquifer's depth. If the depth is over- or under-estimated, then the saltwater intrusion extent will also be over- or under-estimated. In the density-dependent approach inaccuracy regarding the aquifer's depth will have a less severe effect on the results because chloride concentration measurements instead of hydraulic heads are used to model the saltwater intrusion phenomenon. This is considered a drawback of the sharp-interface approach: a direct calibration using chloride measurements is infeasible since only groundwater flow simulation is performed. Consequently, this method provides information only regarding the saltwater intrusion extent and no estimation of water salinity in different points of the aquifer is performed. In contrast, the density-dependent approach allows for chloride mass calibration and provides a concentration distribution in all aquifer points.

DISCUSSION AND CONCLUSIONS

Saltwater intrusion is an important factor that can intensify water quality degradation of coastal aquifers. The physics of the saltwater intrusion system can be viewed in two ways: by the sharp-interface approach, which assumes that the density is constant over the aquifer and a very narrow transition zone exists between freshwater and saltwater; and the density-dependent approach, which assumes that the two fluids are miscible, the density changes according to salt concentration and a wider transition zone exists. The appropriateness of each approach depends on the given problem. In systems that exhibit thin transition zones, the sharp-interface formulation can yield useful results. However, there exist very few real aquifers where density does not vary over a significant portion of the aquifer; thus the applicability of the sharp-interface approach in real-world applications is limited. In areas where a wide transition zone exists, the more complex density-dependent method is the only option that yields meaningful results. However, the density-dependent method requires the knowledge of parameters, in particular dispersivity, that are difficult to estimate and are usually determined empirically or using bibliographic values (Reilly and Goodman 1987). Consequently, a sharp-interface or a density-dependent simulation can be employed depending upon the complexity of the system.

The aquifer investigated in this study is highly karstified and therefore an equivalent porosity model alone could not describe the system accurately. To this extent, a combination of the equivalent porosity method with the discrete feature method, that utilizes 2D elements to describe fractures, was employed in this work. The advantage of this approach is that fractures can be placed in any location of the aquifer, with different orientation and variable characteristics. In this work, a total of ten fractures was included, represented by 2D triangular discrete feature elements and located at the depth that corresponds to the second model layer. Seven of the fractures were the original ones mapped by site experts and shown on the geological map. In order to explain the extremely high chloride concentrations observed in wells around the industrial area, three extra smaller fractures were added with an east to west orientation towards these high concentration wells.

The flow calibration results for the density-dependent approach show a good agreement between measured and model values in most cases (R ² = 0.875) with the exception of wells G3 and G13, which exhibit extremely high hydraulic head values indicating that the measurements were either erroneous or were taken under injection conditions; thus they cannot serve as reliable data. The chloride contaminant transport results show a very good fit (R ² = 0.966) except for wells G1, G2 and G12, which exhibit higher chloride concentration values because they are located very close to the extreme value of 16 480 mg/L (well G3) and thus are affected by it. For the sharp-interface approach flow results, an R ² value of 0.8754 was achieved also indicating a good fit.

In the problem examined in this paper, the calculated sharp interface is located further inland than expected, an observation that is consistent with the findings of Hill (1988). This result may be related to the appropriateness of the 250 mg/L isochlor for the comparison with the 2 m isoline of the sharp-interface approach. The selection of a lower isochlor would yield a better match between the two approaches. The over-estimation of the extent of saltwater intrusion by the sharp-interface approach may also be related to the accuracy of the estimate of the aquifer's depth, which is of paramount importance when using the Ghyben-Hertzberg relationship. If the depth is over- or under-estimated, then the saltwater intrusion extent will consequently be over- or under-estimated too. This is not the case when applying the density-dependent approach, which utilizes chloride concentration measurements instead of hydraulic head values. Based on the above, the over-estimation of the saltwater intrusion front by the sharp-interface method is dependent upon various assumptions and conditions that do not fully apply to all groundwater systems rendering the comparison between the two approaches a difficult task.

As a final remark, in the sharp-interface approach, a direct calibration with chloride concentration measurements is not possible, since no contaminant transport simulation is performed. The only information provided by the method is in terms of saltwater intrusion extent and no detail is obtained regarding water salinity. In contrast, the density-dependent approach is more informative, allowing for chloride mass calibration and providing a concentration distribution for the entire aquifer.