Abstract

The purpose of this study is to investigate the saltwater intrusion phenomenon in the alluvial aquifer of Katapola, on Amorgos Island, under current and future climatic conditions and to provide groundwater management options for alleviating this problem. To this end, a groundwater flow model was developed and the sharp-interface approximation combined with the Ghyben–Herzberg equation was used. A correction factor that accounts for the hydrodynamic dispersion occurring at the brackish zone was also incorporated in the analysis. The model results show that under the current pumping strategy, the saltwater intrusion front extent is vast, posing a serious threat to the quality of groundwater used for drinking and irrigation in the area. The management goal is to find the alternative pumping scenarios for the existing well network that will prevent further spreading of saltwater intrusion. Several water management scenarios were developed, taking into account the effects of climate change, the increase in water supply demand and the expected population growth. The results indicate that controlling the propagation of seawater intrusion in Katapola necessitates the periodic deactivation of most of the pumping wells and the design of alternative plans in order to meet the increasing water demand.

INTRODUCTION

Aegean islands' limited groundwater resources set them among the most sensitive water districts of Greece. Some of the main inhibiting factors affecting groundwater resources are the islands' limited surface areas, complex geological structure, sharp hilly morphology and low storativity combined with low annual rainfall and high average temperatures (Kourmoulis et al. 1992). The planning and management of such aquifers requires special attention, in order to avoid the salinization of the limited freshwater sources. The proximity to the sea and overexploitation of coastal aquifers are conducive to the effect of seawater intrusion. Since the study of Ghyben (1888) and Herzberg (1901), extensive research has been carried out leading to the understanding of the mechanisms that govern seawater intrusion (Bear & Cheng 2010). According to the width of the transmission zone between seawater and freshwater, two approaches have been developed to describe saltwater intrusion mathematically. The first approach is based on the sharp-interface approximation. This approach, neglects mixing and assumes steady state conditions, often leading to the overestimation of the saltwater intrusion extend. It is applied with more accuracy when the width of the freshwater–saltwater transition zone is small relative to the thickness of the aquifer. The second approach is the density-dependent approach, according to which the two fluids (freshwater and seawater) are miscible, thus forming a transition zone of significant width (Pinder & Cooper 1970). The width is dictated by three phenomena: the advection of water towards the sea, the recirculation of seawater and mixed water and the hydrodynamic dispersion (Bear & Cheng 2010).

Numerous groundwater flow models have been developed over the years, taking into account the common sharpness of the fluids' interface, such as Princeton Transport Code (PTC) (Pinder et al. 1997). More recently, models solving the coupled variable-density groundwater flow and solute-transport equation have also been developed, such as FEFLOW, SEAWAT–MODFLOW and SUTRA, among others. The appropriateness of each method and corresponding model depends on the complexity of the problem at sight and the available data (Karatzas & Dokou 2015).

When sharp interface models are used together with the Dupuit approximation, the problem is simplified, facilitating model development. For this reason, sharp interface models are widely employed to study saltwater intrusion. Karatzas & Dokou (2015) used the sharp-interface approach with a particle swarm optimization algorithm and PTC to select an optimal management strategy for the saltwater intrusion problem in Malia, Crete. Papadopoulou et al. (2005) have also applied this approach in conjunction with PTC to simulate the salinity zone in karstified systems, using recharging scenarios and sensitivity analysis to suggest water management solutions.

As mentioned above, particularly in Cyclades, groundwater is the main water resource in many semi-arid coastal regions and water demand, especially during summer months, can be very high. Groundwater withdrawals for meeting this demand, often cause seawater intrusion and degradation of water quality of coastal aquifers (Kourakos & Mantoglou 2011). Extensive research has been made for Naxos Island, in an attempt to study the parameters that control the groundwater flow regime in fractured aquifers (Partsinevelou et al. 2011) but also to propose methods for facing long term water scarcity (Tsakiris & Spiliotis 2011). More recently, Giannikopoulou et al. (2015) worked towards proposing drought mitigation solutions for Syros Island, by combining water balance modelling, hazard analysis and risk and cost effectiveness analysis.

Previous studies on groundwater flow and saltwater intrusion in Greek islands have considered both of the aforementioned methods to describe the saline zone mathematically. Kourakos & Mantoglou (2011, 2015) proposed a coupling of a simulation model (SEAWAT) and an optimization algorithm (genetic algorithm) to maximize pumping rates subject to environmental constraints for Santorini Island, which exhibits the same semi-arid coastal conditions as Amorgos Island. This simulation-optimization coupling method was also applied by Evelpidou et al. (2010) to manage coastal aquifers in Paros Island. Kopsiaftis et al. (2009) used the variable-density model FEFLOW to study an alluvial unconfined aquifer with volcanic sediments in Santorini. The same mathematical modeling tool was combined with geophysical methods by Kourgialas et al. (2016) to develop a saltwater intrusion prediction tool for the prevention of freshwater in agricultural coastal regions. Mantoglou (2003) used the sharp-interface approach with the Ghyben–Herzberg approximation and MODFLOW to study the shallow, coastal aquifer of Vathi, Kalymnos. In Hersonissos, Crete, a more recent study suggested using PTC with the Algorithm of Pattern Extraction (ALOPEX) to optimize the area's pumping activities (Stratis et al. 2015).

Most of the above research, that uses the sharp interface assumption, neglects mixing and implicitly assumes that saltwater remains static. Consequently, this approximation overestimates the penetration of the saltwater front, as has been shown by Dokou & Karatzas (2012) when they compared the two main saltwater intrusion approaches. To quantify the error introduced by this method, Pool & Carrera (2011) performed numerical, three-dimensional variable-density flow simulations based on the formulation of Strack (1976). They concluded that the Ghyben–Herzberg estimate of the interface depth can be corrected by an empirically derived dispersion factor. This new term incorporates transverse dispersity and the aquifer thickness in order to estimate the critical pumping rates in the field. This correction has not been evaluated in many field applications. An exception is the recent work by Beebe et al. (2016) where an analytical model was used and evaluated against observations from Canada, the United States, and Australia to assess its utility as an initial approximation of sea water extent. Beebe et al. (2016) concluded that using the correction by Pool & Carrera (2011) did not markedly improve the results. In contrast, this work shows that applying this correction in the Katapola aquifer improves the results and provides a more realistic representation of the saltwater intrusion front than the original Ghyben–Herzberg formulation. This is confirmed by a direct comparison against chloride concentration data for a variety of wells and time periods.

The work presented here, evaluates this correction using chloride concentration field data and examines how various climatic change scenarios would affect the groundwater resources and saltwater intrusion phenomenon, using a model that has been corrected to account for dispersivity effects. Specifically, the present study focuses on the development of a saltwater intrusion model for the alluvial coastal aquifer of Katapola, in Amorgos Island, using a combination of the sharp interface approximation and the Ghyben–Herzberg relationship corrected by the dispersion factor proposed by Pool & Carrera (2011). The system of partially differential equations is solved using the PTC model, a groundwater flow simulator that combines finite element and finite difference methods. An evaluation of the correction factor is performed by comparing against chloride concentration data, concluding that the correction improves the estimation of the saltwater front location. The system's response to changing climatic conditions, pumping rates and population growth is also examined by applying various scenarios to the calibrated model.

METHODS

For the simulation of the groundwater flow in the aquifer of Katapola, the PTC was employed. In order to solve numerically a system of three-dimensional equations, PTC uses a combination of finite element and finite difference methods. These equations represent groundwater flow through a porous medium, Equation (1), and groundwater pore velocity, Equation (2).  
formula
(1)
where : specific storage and W: source/sink term.  
formula
(2)
where h is the hydraulic head [L], is the hydraulic conductivity [LT−1], v is the pore velocity vector [LT−1] and n is the effective porosity.

PTC employs a unique splitting algorithm for solving the fully three-dimensional equations, reducing the computational burden significantly (Babu et al. 1997). The algorithm involves discretizing the domain into approximately parallel horizontal layers. Within each layer a finite element discretization (Pinder & Gray 1977) is used, allowing for accurate representation of irregular domains. The layers are connected vertically by a finite difference discretization. During a given time iteration, the computations are divided into two steps. In the first step, all the horizontal finite element discretization is solved. In the second step, the vertical equations which link the layers are solved (Babu et al. 1997).

The model also accommodated two types of boundary conditions for the flow equation. The specified head boundary (Dirichlet) of zero meters was used at the coastline and a specified flow condition (Neumann), was used for pumping wells as well as for the lateral flow boundaries. To predict the saltwater intrusion front, the sharp-interface approach was used, in combination with the Ghyben–Herzberg relationship:  
formula
(3)
where is the head between the groundwater table and the sea level [L], is the head between the sea level and the seawater–freshwater interface [L], is the density of saline water (1.025 kg/m3) and is the density of freshwater (1.000 kg/m3).
The term was then multiplied by the dispersion factor , as introduced by Pool & Carrera (2011)  
formula
(4)
where is aquifer thickness [L] and is transverse dispersity [L] which can be estimated as a fraction of longitudinal dispersivity :  
formula
(5)
When no data are available, longitudinal dispersivity is typically calculated by the empirical equation of Xu & Eckstein (1995)  
formula
(6)
where L is a characteristic length; in this case it represents the average distance between the coastline and the active pumping wells [L].
Therefore, the Ghyben–Herzberg approximation to the interface depth is rewritten as:  
formula
(7)

Pool & Carrera (2011) compared the salt mass fraction distributions derived from numerical models and the results from the corrected analytical solution, to validate the applicability of this formulation. The corrected equation can be employed to compute the extent of saltwater intrusion in a more accurate way, for a broad range of conditions, and can therefore lead to more realistic and effective management of coastal aquifers.

Site description

Figure 1 depicts the study area of Katapola Bay in the Aegean island of Amorgos, a land area of just 2.2 km2. The elevation in this area varies from 0 to 100 m above sea level and rises intensely along the perimeter of the region, which consists mainly of carbonate rocks. The study area comprises three communities, Katapola, Rachidi and Ksilokeratidi that have approximately 595 permanent residents. However, during summer periods, the population may rise up to 4,000, due to the fact that it is a highly touristic area. Heterogeneous farmlands and natural vegetation covers most of the study area. The alluvial coastal aquifer, which consists mainly of sand and clay depositions, occupies an area of 0.7 km2 (Papadopoulos & Paritsis 1997). The alluvial zone is laid over a series of flysch, which consists mainly of clay schists that act as an impermeable substrate of the domain (Siaka et al. 2016).

Figure 1

The study area of Katapola bay in the Aegean island of Amorgos.

Figure 1

The study area of Katapola bay in the Aegean island of Amorgos.

Discretization and basic parameters

The vertical discretization of the model area was based on a study by Papadopoulos & Paritsis (1997). The conceptual model comprises an unconfined aquifer that is about 20 m deep and was simulated using one numerical layer, simplifying the model in a two-dimensional form. The elevation of the model's layer follows the topography of the surface (0–100 m). In order to discretize the domain horizontally, a triangular finite-element mesh consisting of 1,999 elements and 1,078 nodes was created.

The PTC model run for a period of 18 years (from September 1997 to April 2015). Each year was divided into two management periods: a dry period consisting of 5 months (May–September) and a wet period consisting of 7 months (October–April). This division was based on the average monthly precipitation values, obtained from the weather station of Naxos, which is close to Amorgos Island (Figure 2).

Figure 2

Rainfall data for the time period 1997–2015.

Figure 2

Rainfall data for the time period 1997–2015.

The rainfall values were also used to estimate groundwater recharge (infiltration) in the model. A typical value for the Cyclades Islands of 15% infiltration was provided by the Greek Ministry of Development (2003–2008). This value was initially applied and then fine-tuned during calibration of the model (Siaka et al. 2016).

The model area includes two geologic formations (Figure 3) with different values of hydraulic conductivity (K), porosity and storativity factor. Information about the ranges of these values was obtained from the River Basin Management Draft of Aegean Islands' Water District. The final values of the above mentioned parameters that were chosen during model calibration are presented in Table 1.

Table 1

Parametric values of hydraulic conductivity, porosity and storativity for the two main geologic formations of the study area

 Geologic formation
Alluvial (al)Flysch (ft)
Hydraulic conductivity 25 m/d 5 10−6 m/d 
Porosity 0.45 0.03 
Storativity 0.06 0.08 
 Geologic formation
Alluvial (al)Flysch (ft)
Hydraulic conductivity 25 m/d 5 10−6 m/d 
Porosity 0.45 0.03 
Storativity 0.06 0.08 
Figure 3

Geological map of the study area (al: alluvial depositions, ft: flysch layer, ft.k: limestone) and model domain (black polygon).

Figure 3

Geological map of the study area (al: alluvial depositions, ft: flysch layer, ft.k: limestone) and model domain (black polygon).

Initial and boundary conditions

A first-type flow boundary condition was set along the coastline of the study area. A change of elevation datum was performed using the bottom of the aquifer as a reference level, at 0 m. Then, the value of 20 m was added to all elevations and hydraulic heads, setting the sea level at 20 m, and topography from 20 to 120 m. The initial conditions for groundwater flow are defined by hydraulic head values, measured in September 1997, for 10 wells and then interpolated over the entire model domain (Table 2).

Table 2

Initial head values, measured in September 1997 at 10 wells

Initial heads [m]
W1W2W3W4W5W6W7W8W9W10
22 20.8 20.2 21.3 21.7 22.7 20.7 20.2 20.3 21 
Initial heads [m]
W1W2W3W4W5W6W7W8W9W10
22 20.8 20.2 21.3 21.7 22.7 20.7 20.2 20.3 21 

Pumping wells were also used as boundary conditions for flow (second-type). As recorded, 20 out of the 27 currently active pumping wells are private and they are used for irrigation purposes only. Most of them pump less than 20 m3/d, with the exception of well F36, that supplies water to an organic nursery crop (100 m3/d) and F17 that pumps approximately 60 m3/d. Code letter F is used for those wells that have been registered since the 2005 study of IGME (Institute of Geology & Mineral Exploration), while those, with letter code Π, have been registered since 2012. Three municipal wells (G14, F34 & W3) pump 100–400 m3 of water daily for drinking supply. G14 is the main supplier of drinking water in Katapola. W3 also supplies water to the municipality of Chora. Another well that pumps a relatively large amount of water for irrigation is W9, which pumps 100 m3/d. Pumping for irrigation is assumed to occur during the dry period, while pumping used for drinking purposes is assumed to reach its highest values during the same period, due to the rise of tourist activity.

Model calibration

Model calibration was performed manually by changing model input parameter values to produce simulated values that match the measured data as closely as possible. Hydraulic conductivity, well pumping rates, groundwater recharge (infiltration) due to precipitation and, most importantly, lateral boundary conditions (inflows) to the model, were the main input parameters used to calibrate the model. For hydraulic conductivities, well pumping rates and groundwater recharge, the initial values used have been reported above. These values were fine tuned to better match hydraulic head results to measured values. Regarding the inflow rates (m3/d), they were defined on four different parts of the model boundaries (southern, south-east and eastern boundaries), as depicted in Figure 4(b), based on hydrogeological considerations. Each boundary was given a different inflow rate coefficient that depends on the extent of precipitation on upstream areas, the upstream geological formations, as well as the recorded hydraulic heads of nearby wells. Most groundwater has been reported to inflow mainly from the southern, southeast and eastern boundary of the domain. Thus, these boundaries were given higher inflow coefficients that were fine-tuned to be 20%, 50% and 20% of the total with inflow rates estimated at 101, 469 and 313 m3/d on average, respectively, for all periods. Their maximum values were 303, 980 and 908 m3/d and their minimum values 21, 117 and 67 m3/d, respectively, depending on seasonal variability. The boundary condition at the northern part, contributes the least as groundwater inflow (10% of the total), due to the limited permeable formations in the upstream geology of this area. Thus, the average flowrate for this boundary condition is 102 m3/d, with an annual minimum of 13 m3/d and 303 m3/d maximum.

Figure 4

Colored contours of the hydraulic heads for (a) wet period 2015 and (b) dry period 2014. Boundary conditions for specified head (coastline) and specified flow (inflows demonstrated by arrows) are also shown. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

Figure 4

Colored contours of the hydraulic heads for (a) wet period 2015 and (b) dry period 2014. Boundary conditions for specified head (coastline) and specified flow (inflows demonstrated by arrows) are also shown. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

To evaluate the agreement between the model results and the field data, model evaluation statistics were calculated as suggested in the work of Moriasi et al. (2007). The calibrated model's metrics are: Coefficient of Determination (R2) = 0.963 and Root-Mean-Square Error (RMSE) = 1.016 m, Normalized Root-Mean-Square Deviation (NRMSE) = 0.0527, Nash-Sutcliffe Efficiency (NSE) = 0.916 which show a good agreement between simulated and measured hydraulic heads. The percent bias (PBIAS) index is negative (−2.237), indicating that model results are slightly overestimated but still within acceptable limits. The hydraulic head output is visualized by colored contours in PTC, creating a pattern that provides insight into the direction and magnitude of groundwater flow throughout the year. As observed from Figure 5, which demonstrates the comparison of calibrated hydraulic heads and existing data, for five different periods and four different wells, G14 and W9 values appear to have the greatest deviation from field measurements. This indicates that the measurements were either erroneous or that there are more pumping wells near G14 and W9 that have not been recorded. Alternatively, it may be assumed that there are locally complex hydrogeological characteristics that the model cannot predict. Figure 4 depicts the colored contours of the hydraulic heads for wet period 2015 (Figure 4(a)) and dry period 2014 (Figure 4(b)). The flysch layer acts as an impermeable barrier for flow, inhibiting hydraulic connection between the southern and the northern parts of the domain. The area around the schist layer also displayed low velocities, in contrast with the eastern part of the domain. Furthermore, the influence of pumping wells G14 (used for drinking purposes) and W9 (irrigation) on groundwater flow, is clearly evident. The model calibration showed that over-pumping takes place by these two wells that abstract approximately 900–1,000 m3 of water per day.

Figure 5

Comparison of calibrated hydraulic heads with the existing data, for five different periods and four different wells.

Figure 5

Comparison of calibrated hydraulic heads with the existing data, for five different periods and four different wells.

RESULTS AND DISCUSSION

Saltwater intrusion estimation

In order to examine the number of wells that pump water of impaired quality, the collected data of were compared with the 250 mg/L chloride concentration, which is the parametric value specified by the European Community directive 98/83/EC on the quality of water intended for human consumption. Then, these values were compared to the saltwater intrusion front, as estimated by the Ghyben–Herzberg relationship. In this case, the depth of the aquifer is assumed only 20 m below sea level, thus, the hydraulic head of freshwater at the interface is 0.5 m (or 20.5 m measured from the aquifer's bottom). The saltwater intrusion front is corrected by the coefficient which was estimated as a function of transverse dispersivity and the thickness of the aquifer (20 m) using Equation (4). Based on this correction, the new hydraulic head of freshwater at the interface is . This value falls within the recorded range of values found by field-scale tracer experiments conducted in an alluvial aquifer in France by Courtois et al. (2000). While it is generally admitted that the pore scale is mainly influenced by grain size, the field scale (macroscopic) longitudinal dispersivity is influenced by factors such as the variance of the log-transformed hydraulic conductivity and the correlation length of the hydraulic conductivities (Domenico & Schwartz 1998). Transverse dispersivity is an important parameter in the correction of the Ghyben–Herzberg equation; but in the absence of direct hydraulic conductivity measurements at the Katapola site, a more sophisticated approach for the estimation of this parameter was not feasible. Therefore, the estimation of transverse dispersivity had to be based only on information found in the literature. Dispersivity values can vary by several orders of magnitude, thus, a sensitivity analysis was conducted to show the effect of the selection of transverse dispersivity values to the calculation of hf. For this purpose, values corresponding to the ranges provided by Domenico & Schwartz (1998) and Courtois et al. (2000) were used. These values span at least two orders of magnitude (0.01 m, 0.1 m, 1 m and 2 m). Values greater than 2 m are not typical for transverse dispersivity, according to the same literature. The resulting hf values are 0.36 m, 0.29 m, 0.19 m and 0.16 m, respectively. The above results show that the choice of transverse dispersivity will have an effect on the estimation of the saltwater intrusion front location, as expected, with an upper bound of 0.36 m and a lower bound of 0.16 m. The 0.2 m contour used in this work represents a value close but less than the average of the above ranges. If a different hf value (in between these bounds) was used, the main conclusions of this work would still hold, but the blue contour in Figures 710 would be shifted further away (0.36 m and 0.29 m) or closer to the coast (0.19 m and 0.16 m), representing a more intense or less intense saltwater intrusion problem, respectively.

Based on the above, it is assumed that as long as every well's hydraulic head remains at least 0.2 m above sea level (20.2 m from the aquifer bottom), the well is not threatened by saltwater intrusion. In order to test this assumption, the hydraulic heads of 13 different wells were compared with the 250 mg/L chloride concentration measured at seven different time periods (1997, 2000, 2005, 2007 and 2013–2015). A correct pairing was observed between the hydraulic head values and the saltwater intrusion front (Figure 6), further validating the choice of transverse dispersivity value used in this work. It is evident from Figure 6, that all wells with chloride concentrations higher than 250 mg/L exhibit hydraulic heads lower than the 20.2 m threshold, except F17 (yellow bar in Figure 6) that has a very high chloride concentration of 808 mg/L, indicating a high level of salinization, while the hydraulic head is slightly over the threshold of 20.2 m (20.5 m measured in April 2007). In addition, it is observed that the majority of pumping wells demonstrate hydraulic head values lower than 20.2 m mainly during the summer periods, when extensive pumping occurs, for all years from 1997 to 2015. It is also quite alarming that G14's water samples seem to exceed the admissible chloride concentration in almost every dry period. The water quality of G14 is of utmost importance, since G14 is the main supplier of drinking water for domestic use in Katapola. Thus, the measured quality of this well is an indicator of the quality of groundwater in Katapola basin, in general. It is also observed that the majority of wells that pump water of impaired quality are either closer to the coastline or closer to wells that seem to over-pump.

Figure 6

Comparison between the values of chloride concentration and hydraulic heads for different time periods (Note: the parametric value for chloride concentration is 250 mg/L and the correspoding hydraulic head threshold for salinization is 20.2 m).

Figure 6

Comparison between the values of chloride concentration and hydraulic heads for different time periods (Note: the parametric value for chloride concentration is 250 mg/L and the correspoding hydraulic head threshold for salinization is 20.2 m).

A visual comparison between the two different fronts is provided in Figure 7. The purple area represents the saltwater intrusion front as estimated by the corrected Ghyben–Herzberg equation. The pink area, is the saltwater intrusion front as estimated from the original Ghyben–Herzberg equation. It is observed that during the wet periods, the saltwater intrusion front is limited near the coastline, whereas, during the dry periods, the front moves inland, affecting a large number of water supply wells (only five wells remain completely unaffected), while reaching the boundaries of the study area. This extreme behavior can be attributed to various factors, such as: the decreased hydraulic gradient during the summer months due to minimal or no groundwater infiltration and the excessive water abstraction, the aquifer's small dimensions and the lack of data for the northeastern part of the study area.

Figure 7

Change of critical heads for the dry (a) and wet (b) period of 2015 (red line: saltwater intrusion front according to Ghyben–Herzberg equation, blue line: saltwater intrusion front according to the corrected Ghyben–Herzberg equation). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

Figure 7

Change of critical heads for the dry (a) and wet (b) period of 2015 (red line: saltwater intrusion front according to Ghyben–Herzberg equation, blue line: saltwater intrusion front according to the corrected Ghyben–Herzberg equation). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

Future scenarios

The final aim of this case study is to examine the behavior of the saltwater intrusion front in the future and how the coastal aquifer of Katapola will respond to climate change (represented by reduced rainfall and subsequent groundwater recharge). The precipitation data trend required was based on the research of Tolika et al. (2007) who examined future precipitation projections from ensembles of RCMs, for various emission scenarios of the Intergovernmental Panel on Climate Change (IPCC). The future scenario Α1Β was chosen for this work, as a worst case scenario, since A1B is characterized as a pessimistic scenario in which CΟ2 concentration will increase up to 815 ppm until the end of the century and on global scale temperature will rise from 2.5 °C to 4.5 °C. In this study, predictions are made up to the year 2050.

Scenario Ι

The first case scenario aims to predict how the saltwater intrusion front will develop in the future, taking climate change into account, if pumping rates remain to current levels. The worst case scenario for precipitation projection over the Greek region suggests a 10.5% reduction during the winter periods and a 34.3% reduction during the summer periods (Tolika et al. 2010). These data were used as input in the groundwater model and the new saltwater intrusion front was estimated. The results are shown in Figure 8.

Figure 8

Depiction of the saltwater for (a) dry and (b) wet expected conditions in 2050 (red line: saltwater intrusion front according to Ghyben–Herzberg equation, blue line: saltwater intrusion front according to corrected Ghyben–Herzberg equation). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

Figure 8

Depiction of the saltwater for (a) dry and (b) wet expected conditions in 2050 (red line: saltwater intrusion front according to Ghyben–Herzberg equation, blue line: saltwater intrusion front according to corrected Ghyben–Herzberg equation). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/ws.2017.160

Comparing projected results for 2050 to the current conditions (2015), it is observed that the projected change in the precipitation rate greatly affects the hydraulic head levels both in the wet and dry periods. Particularly, the saltwater front threatens to cover the biggest part of the study area, and reaches the east and north boundaries of the domain during the dry period of 2050, proving groundwater, as a water source, to be unreliable (Figure 8(a)). It is still clear that the southeast area receives the largest amount of inflows from the upstream of the domain, prohibiting further spreading of the front in that part. The flysch layer also continues to act as an impermeable boundary, preventing hydraulic communication between the northern and the southern regions. Similarly, to a lower but nevertheless considerable extent, future wet periods that seemed to be relatively unaffected under the current conditions, will also incur the spread of the brackish zone, in the area close to the coastline, where more than three active wells are located.

Scenario ΙΙ

In the second scenario, the effect of reduced extraction rates was examined. The pumping reduction method was applied only on those wells that pump over 100 m3/for irrigation or drinking purposes. The new extension of the saltwater intrusion front was predicted for a 25%, 50% and a 90% reduction in pumping rates (Figure 9). Since the problem of groundwater quality degradation becomes more acute during the summer months, the reduced pumping rate scenarios were examined only for the dry periods until 2050. Four wells, two of which extract water for drinking purposes, are recorded to pump more than 100 m3/d (G14: 1,000 m3/d, W3: 100 m3/d, F34: 150 m3/d and W9: 900 m3/d). As a result, 25% and 50% reduction means that the total abstraction rate will not exceed 1,612.5 and 1,075 m3/d, accordingly. It is observed that only the 90% decrease sub-scenario alleviates the problem of the saltwater front's extent completely. However, the implementation of such a scenario suggests that one of the main water supply wells in the area (G14) will be abstracting less than 40 m3 of water per day and that the total water abstraction in the area will not exceed 215 m3/d. This restriction is extreme and its adoption would be unrealistic. But, even the 25% or 50% reduction of pumping rates seem to have a significant effect on the saltwater intrusion extend and could potentially be considered by the water management stakeholders in combination with alternative water resource solutions, to supply the required amount of water in the Katapola area.

Figure 9

The critical heads of the dry period 2050, with a 25%, 50% and 90% reduction of pumping rates.

Figure 9

The critical heads of the dry period 2050, with a 25%, 50% and 90% reduction of pumping rates.

Scenario ΙΙΙ

The third scenario considers the system's response to population growth and, thus, to increased pumping rates, if no other water supply method is applied. Based on the Hellenic Statistical Authority (HAS) and the studies of Papadopoulos & Paritsis (1997) and Kourmoulis (1983), local population in Katapola region, is expected to reach approximately 1,734 residents by 2050, and even exceed 8,000 residents during the summer. Accordingly, the daily water demand will be greatly affected for each period. Specifically, the total daily demand in water supply is expected to rise by 28% by 2050. Figure 10 shows how the saltwater intrusion front will develop in this case for the dry period of 2050 (Figure 10(a)) and for the wet period (Figure 10(b)). This result clearly indicates that the current pumping strategies are not sustainable and alternative water management scenarios need to be considered in the near future.

Figure 10

The critical heads for the (a) dry and (b) wet period of 2050.

Figure 10

The critical heads for the (a) dry and (b) wet period of 2050.

Scenario ΙV

The fourth scenario examines two cases of zero abstraction. In the first case (Figure 11(a)), all wells are shut down except for the main G14. In the second case (Figure 11(b)), all wells are deactivated for the dry periods, to estimate the natural effect of saltwater intrusion. Groundwater recharge and lateral inflow rates remain constant. Figure 11 clearly shows that the deactivation of G14 greatly affects the saline zone location in the southern part of the coastal area. This phenomenon can be explained considering the geological formation of slate, which acts as an impermeable barrier for flow, and inhibits hydraulic communication between the southern and the northern part of the domain. Overall, the full deactivation of pumping wells limits the front close to the coastal zone, however, every well that operates near the coastline will abstract water of low quality.

Figure 11

The critical heads for the dry period of 2050, with (a) partial and (b) full deactivation of the pumping wells.

Figure 11

The critical heads for the dry period of 2050, with (a) partial and (b) full deactivation of the pumping wells.

CONCLUSIONS

A groundwater flow simulation of the alluvial aquifer in the area of Katapola was conducted with acceptable accuracy, given that the model's hydraulic head values deviate 1 m maximum from the corresponding field measurements. The most important wells in the area are G14, which is the main water abstraction well for drinking purposes and W9, the main water abstraction well for irrigation purposes. The model showed that these over-pumping wells, greatly affect not only neighboring areas, but also the largest part of the study area. This result demonstrates how sensitive the groundwater system is to high pumping rates, mainly during the dry (summer) periods.

Regarding the saltwater intrusion, the use of the empirical correction factor led to more accurate predictions of the interface location and, thus, to a more realistic and efficient simulation of the groundwater system. It was found that, during the wet periods, the saltwater front is limited near the coastline, due to higher recharge rates and lower pumping rates. During the dry periods, the front moves further inland and, in extreme cases of drought, seems to reach the boundaries of the study area, posing a serious threat to groundwater quality.

Given the continuously growing touristic activity in the area, it is necessary to design and implement a fully controlled pumping strategy during the future dry periods. In particular, a reduction of 25%, 50% and 90% in pumping was tested for the study area. Only the 90% reduction in pumping rate scenario seemed to alleviate the problem of extended saltwater intrusion, in the southern part of Katapola. Such restriction is deemed unrealistic, thus a more moderate reduction could be applied in combination with alternative water sources. In any case, there are several wells still threatened by the saltwater front in the northern part. In addition, it is evident that the saltwater front is directly affected by rainfall (and subsequently groundwater recharge), which does not exceed 400 mm annually, as derived from the weather station of Naxos. Therefore, the effect of climate change on the front, which projects a considerable reduction in rainfall during the summer months, will exacerbate the problem, even if well abstraction is limited. Specifically, the worst case scenario suggests that, due to population growth, all wells will pump water of impaired quality by 2050, if no water management plan is implemented to mitigate the problem. On the other hand, wet periods may be seen as an opportunity to store water, as most of the wells remain unaffected by the seawater front.

For the points illustrated in the proceeding discussion, it may be assumed that an integrated groundwater management would include the periodic deactivation of most pumping wells. The replenishing water for irrigation could derive from a combination of methods such as water containment barriers, artificial recharge of the aquifer or rainwater recovery. The additional water for drinking purposes could derive from desalination, a method that has already be implemented in other Aegean islands.

REFERENCES

REFERENCES
Babu
,
D.
,
Pinder
,
G. F.
,
Niemi
,
A.
,
Ahlfeld
,
D. P.
&
Stothoff
,
S. A.
1997
Chemical Transport by Three-Dimensional Groundwater Flows
,
84-WR-3
.
Princeton University
,
Princeton, NJ, USA
.
Bear
,
J.
&
Cheng
,
A. H. D.
2010
Modeling groundwater flow and contaminant transport
.
Transport in Porous Media
23
(
1
),
850
.
Beebe
,
C. R.
,
Ferguson
,
G.
,
Gleeson
,
T.
,
Morgan
,
L. K.
&
Werner
,
A. D.
2016
Application of an analytical solution as a screening tool for sea water intrusion
.
Groundwater
54
,
709
718
.
Courtois
,
N.
,
Gerbaux-Francois
,
O.
,
Grenier
,
C.
,
Maugis
,
P.
,
Mouche
,
E.
,
de Fouquet
,
C.
,
Goblet
,
P.
&
Ledoux
,
E.
2000
Characterization of dispersion in an alluvial aquifer by tracing techniques and stochastic modeling
. In:
Calibration and Reliability in Groundwater Modelling
(Proceedings of the ModelCARE 99 Conference held at Zurich, Switzerland, 1999), pp. 84–89
.
Domenico
,
P. A.
&
Schwartz
,
F. W.
1998
Physical and Chemical Hydrogeology
, 2nd edn.
John Wiley & Sons Inc
.,
New York, NY, USA
.
Evelpidou
,
N.
,
Poulos
,
S. E.
&
Vassilopoulos
,
A.
2010
Paros Island (Cyclades, Aegean Sea) Coastal Zone: Natural Processes and Dynamics
.
Springer
,
Dordrecht, The Netherlands
, pp.
285
296
.
Ghyben
,
W. B.
1888
Nota in verband met de voorgenomen putboring nabij Amsterdam
.
Tijdschrift van het Kononklijk Instituut van Ingenieurs The Hague (Note in connection with the proposed wellbore near Amsterdam. Journal of the Institute of Engineers Kononklijk The Hague)
9
,
8
22
.
Giannikopoulou
,
A. S.
,
Gad
,
F. K.
,
Kampragou
,
E.
&
Assimacopoulos
,
D.
2015
Risk-based assessment of drought mitigation options: the case of Syros Island, Greece
.
Water Resources Management
31
(
2
),
655
669
.
Herzberg
,
B.
1901
Die WasserversQrgung einiger Nordseebader
.
J. Gasbeleuchtung und Wasserversorgung (The water supply of some North Sea resorts. Gas Lighting and Water Supply)
44
,
815
819
,
842
844
.
Kopsiaftis
,
G.
,
Mantoglou
,
Α.
&
Giannoulopoulos
,
P.
2009
Variable density coastal aquifer models with application to an aquifer on Thira Island
.
Desalination
237
(
1–3
),
65
80
.
Kourakos
,
G.
&
Mantoglou
,
A.
2011
Simulation and multi-objective management of coastal aquifers in semi-arid regions
.
Water Resources Management
25
(
4
),
1063
1074
.
Kourakos
,
G.
&
Mantoglou
,
A.
2015
An efficient simulation-optimization coupling for management of coastal aquifers
.
Hydrogeology Journal
23
(
6
),
1167
1179
.
Kourgialas
,
N.
,
Dokou
,
Z.
,
Karatzas
,
G. P.
,
Panagopoulos
,
G.
,
Soupios
,
P.
,
Vafidis
,
A.
,
Manoutsoglou
,
E.
&
Schafmeister
,
M.
2016
Saltwater intrusion in an intensively irrigated agricultural area: combining density-dependent modeling and geophysical methods
.
Environmental Earth Sciences
75
,
15
.
Kourmoulis
,
Ν.
1983
Υδρογεωλογική Έρευνα Κυκλάδων ΙΙ Αμοργός (Hydrogeological Study of Cyclades II Amorgos)
.
Υδρολογικές και υδρογεωλογικές έρευνες
35
,
1
27
.
Kourmoulis
,
Ν.
,
Papapetros
,
P.
&
Charmanidis
,
F.
1992
Η συμβολή του Ι.Γ.Μ.Ε. στην επίλυση των υδρευτικών προβλημάτων των νησιών του νοτίου Αιγαίου (IGME contribution to the solution of water problems in S. Aegean islands). http://www.igme.gr/index.php/9-igme/18-i-vivliothiki
(
accessed 3 December 2016
).
Moriasi
,
D.
,
Arnold
,
J.
,
Van Liew
,
M.
,
Bingner
,
R.
,
Harmel
,
R.
&
Veith
,
T.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Transactions of the American Society of Agricultural and Biological Engineers
St. Josef MI.
50
(
3
),
885
900
.
Papadopoulos
,
N.
&
Paritsis
,
S.
1997
A study of water containment barriers in Cyclades – Final Geological – Hydrogeological Study, Athens, Greece
.
Papadopoulou
,
M. P.
,
Karatzas
,
G. P.
,
Koukadaki
,
M.
&
Trichakis
,
Y.
2005
Modeling the saltwater intrusion phenomenon in coastal aquifers–a case study in the industrial zone of Herakleio in Crete
.
Global NEST Journal
7
(
2
),
197
203
.
Partsinevelou
,
A. S.
,
Lozios
,
S.
,
Stournaras
,
G.
2011
Fracture pattern description and analysis of the hard rock hydrogeological environment in Naxos Island, Hellas
. In:
Lambrakis
,
N.
,
Stournaras
,
G.
&
Katsanou
,
K.
(eds),
Advances in the Research of Aquatic Environment
.
Springer, Berlin & Heidelberg
,
Germany
, pp.
45
52
.
Pinder
,
G. F.
&
Cooper
,
H. H.
1970
A numerical technique for calculating the transient position of the saltwater front
.
Water Resources Research
6
(
3
),
875
882
.
Pinder
,
G. F.
&
Gray
,
W. G.
1977
Finite Element Simulation in Surface and Subsurface Hydrology
.
Academic Press
,
New York, NY, USA
, p.
295
.
Pinder
,
G. F.
,
Niemi
,
A.
,
Ahlfeld
,
D. P.
&
Stothoff
,
S. A.
1997
Chemical Transport by Three-Dimensional Groundwater Flows
.
Princeton University
,
Princeton, NJ, USA
,
84-WR-3
.
Siaka
,
M.
,
Dokou
,
Z.
&
Karatzas
,
G. P.
2016
A study of groundwater flow and saltwater intrusion at the alluvial aquifer of Katapola, at the island of Amorgos, Greece
. In:
Paper Presented to 2nd EWaS International Conference
,
Chania, Crete
,
Greece
.
Stratis
,
P. N.
,
Dokou
,
Z.
,
Karatzas
,
G. P.
,
Papadopoulou
,
E. P.
&
Saridakis
,
Y. G.
2015
Stochastic optimization and numerical simulation for pumping management of the hersonissos freshwater coastal aquifer in Crete
.
Applied Water Science
7
(
5
),
2425
2435
.
Tolika
,
K.
,
Maheras
,
P.
,
Vafiadis
,
M.
,
Flocas
,
H. A.
&
Arseni-Papadimitriou
,
A.
2007
Simulation of seasonal precipitation and raindays over Greece: a statistical downscaling technique based on artificial neural networks (ANNs)
.
International Journal of Climatology
27
,
861
881
.
Tolika
,
C. K.
,
Zanis
,
P.
&
Anagnostopoulou
,
C.
2010
Regional climate change scenarios for Greece: future temperature and precipitation projections from ensembles of RCMs
.
Global NEST Journal
14
,
407
421
.
Tsakiris
,
G.
&
Spiliotis
,
M.
2011
Planning against long term water scarcity: a fuzzy multicriteria approach
.
Water Resources Management
25
(
4
),
1103
1129
.