Environ Model Assess DOI 10.1007/s10666-016-9529-z

Utilizing Successive Linearization Optimization to Control the Saltwater Intrusion Phenomenon in Unconfined Coastal Aquifers in Crete, Greece Zoi Dokou 1 & Maria Dettoraki 1 & George P. Karatzas 1 & Emmanouil A. Varouchakis 1 & Athina Pappa 1

Received: 18 November 2015 / Accepted: 28 July 2016 # Springer International Publishing Switzerland 2016

Abstract In the present work, a simulation-optimization method is employed in order to manage saltwater intrusion in two unconfined coastal aquifers in Crete, Greece. The optimization formulation seeks to maximize groundwater withdrawal rates while maintaining the saltwater intrusion front at the current location or inhibiting it at locations closer to the coast. A combination of a groundwater flow model (MODFLOW) with the Ghyben-Herzberg saltwater front approximation and a sequential implementation of the Simplex algorithm (GWM) are employed. The results show that under the current pumping strategies, the saltwater intrusion front will continue to move inland, posing a serious threat to the groundwater quality of these regions. Optimal groundwater withdrawal scenarios that take into consideration the water needs of the local communities and environmental concerns are presented and discussed. In both case studies, significant reductions in pumping are required in order for the saltwater intrusion front to retract closer to the shoreline.

* Zoi Dokou [email protected] Maria Dettoraki [email protected] George P. Karatzas [email protected] Emmanouil A. Varouchakis [email protected] Athina Pappa [email protected] 1

Technical University of Crete, School of Environmental Engineering, Polytexneioupolis, 73100 Chania, Greece

Keywords Saltwater intrusion . Groundwater management . Coastal aquifer . Optimization . Greece

1 Introduction The sustainable management of water resources has become an increasingly important issue in the last decades. Many recent research studies have focused on the management of coastal aquifers, due to their vulnerability to contamination by saline water. This phenomenon, termed as saltwater intrusion, poses a risk to groundwater quality if not managed in time. The extent and shape of the saltwater zone depends on various parameters such as the dispersive mixing, tidal effects, density effects including unstable convection, and geological characteristics, as well as aquifer hydraulic and transport properties [1]. More importantly though, the saltwater intrusion phenomenon is caused by excessive groundwater pumping activity; therefore, the development of sustainable pumping strategies is of paramount importance to protect this resource. Mathematical models can offer a great tool for achieving this goal [2]. In order to describe the evolution of the saltwater front in a coastal aquifer mathematically, there are two main approaches: the density dependent and the sharp interface approach. The density-dependent flow approach assumes that a mixing zone of significant width between freshwater and saltwater exists [3]. Researchers have successfully applied this method on real field sites [4–8, among others]. However, density-dependent models of saltwater intrusion require a substantial amount of data for model parameterization and may lead to complex and computationally demanding applications. Long simulation times complicate and in many instances even prohibit the

Dokou Z. et al.

application of simulation-optimization problems due to the large number of simulation calls to the optimizer. The sharp-interface approach on the other hand, which is usually combined with the Ghyben-Herzberg relationship, is the most simplified, widely used approach that requires far less parameters and can be very useful for initial analysis purposes. It assumes that freshwater and saltwater (two immiscible fluids of different constant densities) form a mixing zone of very small finite width [9]. This interface can, under certain circumstances (such as steady-state conditions, short time scales, and narrow transition zones relative to the field scale of the problem) be neglected. Under transient conditions though, the freshwater and saltwater dynamics control the behavior of the interface [10]. This approach has been used by many researchers [2, 10–16, among others]. This alternative approach is often used in order to overcome the computational issues arising in simulation-optimization approaches that employ density-dependent numerical models. For the development of a sustainable management strategy, selecting the appropriate optimization model that will be coupled with the simulation model is very important. Various objectives can be considered, depending on the goals of the stakeholders, such as economic, environmental, and societal aspects [2]. Three main types of optimization models have been used by researchers for managing saltwater intrusion, depending on the complexity of each problem, namely linearization methods [14, 17, 18, among others], nonlinear programming [19–23], and more recently heuristic algorithms [2, 14, 24–28]. Their main disadvantage is that they require a large number of simulation model calls, rendering the process computationally expensive. The goal of this work is to assess and manage the extent of saltwater intrusion in two unconfined coastal aquifers, located in Crete, Greece. The approach employed here is to combine a groundwater simulation model (MODFLOW) with a simple and computationally inexpensive optimization method, the sequential implementation of the Simplex algorithm of the groundwater management model (GWM). The management formulation objective is to maximize groundwater pumping while hindering saltwater intrusion in the coastal aquifer. Optimal pumping schemes that take into consideration both the water needs of the local communities and environmental concerns are presented and discussed. The contribution of this work to the field is based on the fact that the GWM model has not been widely used; thus, this work can serve as a verification of the algorithm using not just example problems but actual field case studies in semi-arid Mediterranean regions. Secondly, the verification is performed in two types of aquifers (a karstic and a porous aquifer). A comparison between the two is also provided that could be of interest to researchers.

2 Methodology 2.1 Groundwater Simulation Model The three-dimensional flow simulator Visual-MODFLOW was used for the development of the three-dimensional numerical model of each coastal aquifer of interest, based on data obtained from the literature and field measurements. VisualMODFLOW incorporates subroutines for defining the boundary conditions and the hydraulic and topographic properties of the aquifer. The governing equation for the subsurface regional flow is solved over the entire domain using a finite difference approximation. The domain of interest is subdivided into cells created by a grid of mutually perpendicular lines that may be variably spaced. A finite difference approximation for the flow equation is written horizontally for each cell and another vertically, connecting the different model layers. The resulting system of finite difference equations is solved, and hydraulic heads are calculated at the center of the cells. The thickness and properties of each layer varies depending on the geological stratification of the area of interest [16, 29–30]. 2.2 Sharp Interface Approach For the estimation of the saltwater intrusion front location, the sharp interface approach was employed in this work. According to this approach, the mixing zone between the two immiscible fluids (freshwater and saltwater) has a small finite width. The aquifer is modeled using only the standard groundwater flow equation (no transport equation is needed), and the seawater intrusion front is approximated hydraulically using the Ghyben-Herzberg relationship: ξ¼

ρf h f ≈40h f ρs −ρ f

ð1Þ

where ξ = interface depth below the sea level; hf = hydraulic head of the freshwater above the sea level; ρf = density of freshwater (1000 kg/m3); and ρ s = density of saltwater (1025 kg/m3). According to the above equation, the location of the sharp interface can be determined by the difference between the hydraulic heads of the saltwater and freshwater and their densities. For both aquifers considered in this work, given that their depth was estimated about 100 m based on boring log information, hf is calculated as 2.5 m. This means that contour of 2.5 m is the hydraulic head isolevel limit below which a zone is considered intruded by saltwater. The use of the sharp interface approximation is reasonable in regional-scale problems when the transition zone is narrow relative to the scale of the problem [31]. In addition, it is more appropriate for modeling short-term responses in aquifers where seawater is able to enter or exit the aquifer easily [10].

Saltwater Intrusion Control in Cretan Coastal Aquifers

2.3 Optimization Problem Formulation The objective of the proposed optimization algorithm is to maximize the total water withdrawn from the active pumping wells in the area while controlling the location of the saltwater intrusion front using hydraulic head constraints. The constraints are imposed at preselected control points where the calculated hydraulic head at the end of the management period should be greater than a minimum value (hmin) that is calculated by the Ghyben-Herzberg relationship and represents the desired saltwater intrusion front location based on the sharp interface approach. The maximum allowable pumping rates for all the production wells were set equal to the pumping rates that are being currently used. Based on the above, the mathematical formulation of the problem is described by as follows: n X Qi ; i ¼ 1; :::n ð2Þ max i¼1

where H is the hydraulic head at the desired location, R is the distance of this location from the pumping well, Q is the pumping rate, K is the hydraulic conductivity of the aquifer, and r is the distance from the pumping well at which the hydraulic head, h, is calculated. GWM solves the nonlinear problem by repeated linearization of the nonlinear features of the management problem (sequential linearization problem (SLP)), where response coefficients are recalculated for each algorithm iteration [32]. It is based on repeated linearization of the nonlinear features in the management problem and is implemented by recalculating the response matrix for each sequential linear program. The first-order Taylor series expansion for head is assumed to be accurate for each sequential linear program, but in contrast to the linear case, the vector of base flow rates changes at each iteration. The head at any location is estimated by hi; j;k;t ðQwÞ ¼ hv i; j;k;t ðQwv Þ

s.t. h j ≥ hmin ; j ¼ 1; :::m

ð3Þ

Qi ≤ Qimax

ð4Þ

Qi ≥ 0

ð5Þ

where n is the number of production wells, m is the number of control points where a head constrain is imposed, Qi is the pumping rate of well i, Qimax maximum allowable extraction rate for each well, hj is the hydraulic head at a control point j, and hmin is the minimum hydraulic head needed to avoid saltwater intrusion in the aquifer calculated by the GhybenHerzberg equation.

2.4 GWM Model GWM is a groundwater management process that is combined with the groundwater flow simulator MODFLOW. GWM uses a response-matrix approach to solve several types of linear, nonlinear, and mixed-binary linear groundwater management formulations. MODFLOWis used in order to calculate the change in head ateach constraintlocation thatresults from a perturbation ofa flow rate variable; these changes are used to construct the response matrix and calculate the response coefficients [32]. The response coefficients are a critical link between the physics of the groundwater flow system and the results of the management model. When optimizing pumping rates, the resulting formulation is slightly nonlinear, a fact attributed to the nonlinear relation of the groundwater heads to the pumping stress under unconfined conditions [14]. To interpret the aforementioned nonlinearity physically, one has to consider the Dupuit-Forchheimer equation: Q r ln h2 ¼ H 2 þ ð6Þ πK R

þ

XN

  ∂hv i; j;k;t ðQwv Þ Qwn −Qwvn n¼1 ∂Qwv n

ð7Þ

where the superscript v represents an iteration level, so that hvi , j , k , t is the head obtained when the set of withdrawal and injection rates Qwv is applied. For values of Qwv that are close to Qwv, the error in this approximation is small, and the values of hi , j , k , t predicted by this equation are relatively accurate. At each iteration of the SLP algorithm, a linear program is constructed on the basis of the first-order Taylor series approximation and the associated response coefficients and is solved using the simplex algorithm. Due to the fact that the head and streamflow responses may be nonlinear, the response coefficients may no longer be constant; thus, they must be recalculated at each iteration. This calculation uses a new vector of base-condition flow rates that is derived from the optimal flow rates obtained from the linear program solution from the prior iteration [32]. The sequential process is continued until two convergence criteria are met: (a) the change in flow rate values from the prior iteration to the current iteration becomes less than a fraction of the magnitude of the flow rate at the current iteration and (b) the change in objective function value becomes less than a specified fraction of the magnitude of the objective function value [32]. Because of the nonlinearity of the system, care should be taken to choose small perturbation values, to ensure adequate precision in the response matrix. In this work, several runs were performed with different perturbation values ranging from −1 to 1, both for the initial and minimum perturbation coefficients required by the SLP problem. The solution remained the same, ensuring the accurate calculation of the response matrix.

Dokou Z. et al.

3 Description of Field Sites Two coastal areas in Crete, Greece, were used as case studies in this work. In both cases, the saltwater intrusion problem is intense because the natural phenomenon has been has augmented by the overexploitation of the aquifers. The selected sites are located in Tympaki and Hersonissos. 3.1 Coastal Aquifer of Tympaki, Crete The Tympaki aquifer is located southwest of Heraklion, Crete, and is a part of the Messara valley, which constitutes the most important agricultural region of Crete. The main source of income in the region is agriculture, an activity that demands large amounts of water, especially during the summer. Greenhouse crops and olive trees constitute the irrigated area of 40,000 m2. The Tympaki basin covers an area of about 55 km2 with one river, Geropotamos, and five smaller streams, flowing through it [33]. Geologically, the area is covered by alluvial formations of broad and fine-grained materials with different hydrogeological and morphological characteristics. 3.2 Coastal Aquifer of Hersonissos, Crete The Hersonissos aquifer is located on the north coast of Heraklion, 25 km from the city of Heraklion in Crete, Greece. The Hersonissos basin covers an area of about 18 km2 and stretches for 3.8 km in W–E direction and almost 4.7 km in N–S direction. Tourism and agriculture are the main occupations in the region. Especially during the summer period water demand is high due to the above activities, leading to high pumping rates that cause a massive drop in the hydraulic head. This intensifies the existing problem of seawater intrusion. The basin is covered mainly with karstified limestones of variable hydraulic conductivity and marls, whereas along the coastal line, alluvial deposits with high permeability can be found [29]. The exact location of the two study areas on the island of Crete is presented in Fig. 1.

4 Groundwater Flow Models Development 4.1 Tympaki Aquifer In order to model and optimize the existing saltwater intrusion problem, two simulations of the groundwater flow model were implemented, one for the wet and one for the dry period, in order to examine the problem throughout a hydrological year. The first simulation covers the period between October and April, while dry period starts in May and ends in September. Results derived from the wet period were used to calibrate the model in the dry period, which was the period of major concern. The hydrological year modeled, using the

available reference data, was the year 2000–2001. In order to construct the groundwater flow model, the unconfined aquifer was discretized in a 50 × 50 mesh grid. The area was divided into three vertical layers and the elevation of the region was defined via a GIS application. Initial conditions, boundary conditions, hydraulic conductivity values, and precipitation data were also used in order to calibrate the model and simulate the groundwater flow in the region. The interpolation method used both for elevation and conductivity, is inverse distance, while steady-state flow conditions were used. In order to simulate the groundwater flow of this region, the density of the water is considered constant and there is no pollutant transport during the simulation period. For initial head, the value of 4 m was applied as an average hydraulic head based on information from the observation wells in the area. Pumping and observation wells must also be taken into consideration, with their main period of operation being during spring and summer in order to cover the irrigation needs of the region. There are 27 pumping wells, the majority of which are located at the centre of the basin, with pumping rates for the dry and wet periods presented in Fig. 2. Rainfall input data were also considered, knowing that the average annual rainfall in the area is 474 mm/year. For the wet period, 70 % of the total rainfall (331.8 mm) was considered as infiltration, while in the dry period, only 30 % (142.2 mm). In order to determine the rainfall that corresponds to each cell, both dry and wet periods’ rainfall values were divided by the total number of cells in the simulated area. The area of interest consists of alluvial formations of variable hydraulic conductivity. The main formations are the alluvial of broad material, with hydraulic conductivity of 5 × 10−4 m/s and the alluvial of fine-grained material with a value of 5 × 10−7 m/s. In the mainland and between the main two formations, there are alluvial formations with hydraulic conductivity of 8 × 10−6 m/s. In addition, there is also a smaller formation in the southeast of the area of interest, with hydraulic conductivity of 1 × 10−5 m/s. Boundary conditions of first type (controlled hydraulic head) were defined along the study area borders that vary significantly during the wet and the dry periods. Along the shoreline, the hydraulic head remains constant and equal to zero, so a first-type flow boundary condition was set. The pumping wells were considered as second-type boundary conditions with pumping rates defined by previous studies. 4.2 Hersonissos Aquifer A 6-month groundwater flow simulation (May to September) was implemented for the Hersonissos aquifer, in order to examine the saltwater intrusion behavior. The unconfined aquifer of interest was considered as one vertical layer with a mesh-grid consisting of 105 columns and 106 rows. The geology of the area consists mainly of limestone formations and

Saltwater Intrusion Control in Cretan Coastal Aquifers Fig. 1 Map of Crete, Greece, with the aquifers of interest designated (Tympaki and Hersonissos aquifers)

sand. The hydraulic conductivity values used for each geologic formation are 12.96 m/day for limestones and dolomitic limestones, 5.2 m/day for bioclastic limestones, 0.15 m/day for marl formations, 0.6 m/day for clay and 430 for alluvial deposits located near the coastline [14]. Regarding the model boundary conditions, a fixed head of 100 m was applied along the coastline to simulate the sea boundary and during the calibration process, the fixed head value of 108 m was chosen along the south side of the basin, in order to create the observed hydraulic gradient and match the measured hydraulic head values.

Fig. 2 Pumping rates for the wet and dry periods for the Tympaki aquifer (values from Paritsis [33])

The region of interest has five pumping wells that are active only during the dry period, in order to meet the irrigation needs and the demand of the population that is increased during the summer period due to tourism. The contour of 2.5 m is the hydraulic head isolevel that represents the saltwater intrusion front for this aquifer as well. However, the reference level was changed from the sea level to the bottom of the aquifer. Thus, if the aquifer bottom elevation is set to zero, then the sea level will be at 100 m. This means that the contour of 102.5 m indicates the location of the saltwater front in the Hersonissos study area.

Dokou Z. et al.

It is important to note here that in the two study areas, the hydrological conditions as well as external stresses to the aquifer have not been modified significantly over the past 15 years. Specifically, the agricultural activities have remained the same (similar crops have been planted in the regions over the years) and a ban on the construction of new wells is in effect; thus, the pumping activity in the area has remained more or less the same. The above evidence suggests that the conditions of the 2000–2001 period have not changed significantly; thus, the results of the optimization algorithm are representative of the current conditions as well.

5 Optimization Results 5.1 Tympaki Aquifer After the definition of all necessary data, the calibration of the flow field was performed, based on the hydraulic head data measured in the observation wells as well as the results of a study conducted by Paritsis [33] in the Tympaki basin. Due to the fact that hydraulic head measurements were only available for the wet season, the model calibration was performed for the wet season, and the results were used as initial conditions for the dry season, in order to estimate the saltwater intrusion front location during this season. The hydraulic heads at the end of each season are presented in Fig. 3. The dry season hydraulic head results are the ones utilized for the construction of the response matrix needed for the management scenarios presented next. As described above, the contour of 2.5 m is the toe of the saltwater intrusion wedge in this study according to the Ghyben-Herzberg approach. For the wet period, the saltwater intrusion is located about 2.3 km from the shoreline, thereby affecting two pumping wells. Respectively, for the dry period, the problem becomes very intense, with the saltwater intrusion reaching 5.8 km from the coastline. In this case, 22 pumping wells of the total of 27 of the region of interest fall inside the saltwater affected area. The purpose of this study is to investigate the current extent of saltwater intrusion, as well as to find the optimal pumping scheme, for which the saltwater intrusion will retreat, while meeting at least part of the irrigation needs of the region. According to the results for the current conditions in the Tympaki aquifer, the location of the saltwater intrusion front exhibits a large difference during the hydrological year, from wet to dry period. During the wet period, the saltwater intrusion front retracts because the pumping activity is greatly reduced and there is increased rainfall; during the dry period, the saltwater intrusion phenomenon intensifies and the front moves farther inland because all the wells are in operation with high pumping flow rates. Even in the case where all the pumping wells are shut down (best-case scenario in terms of

retreat), the saltwater intrusion front retracts but not very close to the coastline, due to the reduced rainfall activity. Four scenarios were considered in order to investigate the possibility of retreat of the saltwater intrusion front toward the shoreline. These scenarios seek to maximize pumping activity, while allowing the saltwater intrusion front to retreat toward the shoreline in prespecified locations. The locations of control points represent the worst case, namely the current situation, and are gradually moved toward the best-case scenario in terms of retreat, which represents the result of subsurface flow with zero pumping activity during the summer. This scenario is used only for comparison reasons since it is would not be desirable in practice due to the grave implications this would have on tourism and farming activities in the area. Due to excessive pumping, the distance of the control points between the worst-case (the current situation) and the best-case scenario in terms of retreat is not very large. This indicates that the saltwater intrusion phenomenon is a slow process that is not fully reversible, at least not within a short period of time. For this reason, it is better to focus on prevention of the phenomenon (to the extent possible), rather than its mitigation. In the MODFLOW environment, 14 control points were considered where the hydraulic head needs to be maintained higher than 2.5 m amsl for each retreat scenario. The corresponding points for each scenario are shown in Fig. 4. For each scenario, the optimized well pumping rates are presented in Fig. 5 and are compared to the current conditions. Comparing the results from the first scenario to the current conditions, there is an almost 60 % reduction of the total pumping rate (from 27,432 to 10,935 m3/day). The number of wells with zero pumping activity is 17, while there is a noticeable reduction in well 24, which reaches 70 %. Most of the wells are located near the 2.5 m contour; thus, they are highly affected during optimization. While the saltwater intrusion front location has not retracted too far from the initial state, the changes in pumping activity are large. This indicates the intensity of the saltwater intrusion problem in the area and the overpumping that occurs in order to meet the needs for drinking water supply and irrigation. The second scenario has a more drastic effect on the reduction of the pumping rates with 19 pumping wells closing down completely. There is a 64 % reduction in the total pumping rate (reduced to 9874 m3/day), while the reduction in well 21 is almost 6 %. As expected, as the saltwater intrusion is retracted toward the shoreline, the total pumping rate decreases, and the number of closed wells increases. The third scenario indicates a 68 % reduction in the total extraction rates (reduced to 8288 m3/day), while the number of wells with no pumping activity is 16. Although the saltwater intrusion front continues to retract toward the shoreline and the total pumping rate decreases, the number of pumping wells with no pumping activity decreases. Specifically, wells W6 and W8 have a small pumping rate of 74.3 and 23 m3/day,

Saltwater Intrusion Control in Cretan Coastal Aquifers Fig. 3 Hydraulic head contours during a the wet and b the dry period for the Tympaki aquifer

respectively, whereas in the first two scenarios, they had zero pumping activity. Well W26 reopens with a 700 m3/day, a rate equal to the upper bound. However, well W21, which is a well with high pumping activity (initial value equal to 3000 m3/ day) needs to shut down completely according to this scenario. For this reason, the total pumping rate decreases although the number of active wells increases. Comparing the results of the fourth and last scenario with the current conditions, there is a 71 % reduction in total pumping rate (reduced to 8055 m3/day). The number of wells with no pumping activity is 17, while there is a noticeable change in the pumping activity in only two wells. Well W6 has a rate of 240 m3/day, a 40 % reduction from the initial, while W24 is almost active with a rate of 6.1 m3/day. As for the distribution of the wells in the region of interest, all of the wells that are located in the central part of the study area have

no pumping activity. In contrast, the wells with no change are away from the mainland, with six of them being in the west side and close to the shoreline and the remaining four in the east side of the area. 5.2 Hersonissos Aquifer The hydraulic head model results after calibration are shown in Fig. 6. The results of the flow simulation of the Hersonissos basin lead to the conclusion that the phenomenon of saltwater intrusion is prominent in the region, with the saltwater intrusion front having advanced 1 to 2.5 km inland. Between the worst- and the best-case scenarios in terms of retreat (red and green lines in Fig. 6), three possible saltwater intrusion retreat locations were considered, each corresponding to a different scenario (scenarios 1, 2, and 3). Eight control points were

Dokou Z. et al.

Fig. 4 Control points for all retreat scenarios for the Tympaki aquifer

placed along the desired location of the saltwater intrusion front, in order to use as constraint points in the optimization formulation and observe the impact on the pumping activity. The constraint locations are shown in Fig. 7 for each scenario and the optimization results in Fig. 8. The results from the first scenario indicate a reduction in the total pumping activity of about 23 % (from 7562 to 5817.8 m3/day). The control points where the front is retracted are not located very far from the initial saltwater front location, and therefore, the difference between the two conditions (current and scenario 1) is not great. Well W2, is the only well that is affected and the reduction of its pumping rate reaches 92.5 %. The key factor that affects the rate of the reduction of the pumping rate of a well is its position in relation to the initial position of the saltwater intrusion. In addition, the reduction of the pumping activity is related to the rate of pumping. Well W2 is located very close to the saltwater intrusion front and is also one of the wells with the largest pumping rate. Both the above reasons contribute to the large reduction of its pumping rate for the first retreat scenario. Comparing the results from the second scenario with the initial conditions, there is a 54 % reduction of the total pumping rate (reduced to 3494.9 m3/day). Well W2 is affected

even more than the first scenario, and its pumping rate is further reduced. This reduction almost turns the pumping well inactive with the optimized rate being only 1.7 % of the initial rate. In addition, there is a 61 % reduction in the pumping rate of well W4, which is a well with high pumping activity (initial value equal to 2520 m3/day). The third and last scenario indicates a total reduction in the pumping activity, which reaches a very high level of 87 % (total pumping rate reduced to 1007.2 m3/day). There are two wells that become inactive (wells W2 and W3) while well W4 exhibits a 79 % reduction in its pumping rate. This scenario has a more drastic effect on the wells and their pumping activity, because of the further relocation of the control points toward the shoreline. Wells W1 and W5 are not affected by the further retreat of the saltwater intrusion front because they are both located away from the 2.5 m isoline (saltwater intrusion toe). 5.3 Sensitivity Analysis A sensitivity analysis on the solutions obtained for each scenario in each aquifer of interested was performed. The binding and nonbinding constraints were identified and their shadow

Saltwater Intrusion Control in Cretan Coastal Aquifers Fig. 5 Optimal pumping rate results for each retreat scenario for the Tympaki aquifer

prices calculated in each case, in order to assess the sensitivity of the objective function to changes in the constraints. The binding constraints are the ones that restrict the value of the objective function, because they prevent the decision variables

from taking values that further improve the optimal solution. Range analysis was used in order to determine the range of values of the right-hand side of the constraints over which the solution does not change.

Dokou Z. et al. Fig. 6 Hydraulic head (m) distribution for the dry (pumping) period for the Hersonissos aquifer

5.3.1 Tympaki Aquifer

5.3.2 Hersonissos Aquifer

For the Tympaki aquifer, the number of near-binding constraints increases as the target saltwater intrusion front moves toward the coast. For scenarios 1 and 2, only 1 out of the 14 constraints is near-binding (scenario 1: constraint 12 and scenario 2: constraint 11). The range of righthand side (RHS) values over which the solution does not change for these constraints is 2.49–2.51 m (scenario 1) and 2.49–2.55 m (scenario 3). For scenarios 3 and 4, three constraints (scenario 3: constraints 6, 11, and 13 and scenario 4: constraints 2, 3, and 10) are near-binding. The range of right-hand side values over which the solution does not change for these constraints is very limited in all cases (scenario 3: 2.49–2.5 m for all constraints and scenario 4: 2.49–2.5 m for constrains 2 and 3 and 2.49–2.53 m for constraint 10), indicating the high sensitivity of the solution to the binding constraints, where a slight change on the right-hand side will alter the optimal solution. This is also evident, if the distances to the right-hand side of each constraint are calculated. The distances for the nearbinding constraint are almost zero, while the ones for the nonbinding constraint vary as shown in Table 1. These distances become smaller as the target saltwater intrusion front moves toward the shore, because the hydraulic head values are smaller toward the coast, as expected.

For all scenarios, only one of the eight constraints (constraint 1 for scenarios 1 and 2 and constraint 6 for scenario 3) is nearbinding while the rest are satisfied by the optimal solution. The range of right-hand side values of the near-binding constraint over which the solution does not change becomes increasingly small as the target saltwater intrusion front moves toward the shore (scenario 1: 102.25–102.66 m, scenario 2: 102.4–102.56 m, and scenario 3: 102.5–102.51 m), indicating the high sensitivity of the solution to the binding constraint, especially for scenario 3, where a slight change on the righthand side will alter the optimal solution. This is also evident, from the distances to the right-hand side of each constraint, shown in Table 2.

6 Discussion and Conclusions In this work, two coastal aquifers in the island of Crete, Greece (Tympaki and Hersonissos aquifers), affected by saltwater intrusion were investigated in order to provide sustainable management strategies of this complex phenomenon. For this reason, the current extend of the problem was studied and then various retreat scenarios of the saltwater intrusion front were investigated in order to

Saltwater Intrusion Control in Cretan Coastal Aquifers

Fig. 7 Control points for the three retreat scenarios for the Hersonissos aquifer

Fig. 8 Current conditions and optimization results for the three retreat scenarios for the Hersonissos aquifer

Dokou Z. et al. Table 1 Constraint

Sensitivity analysis on optimization constraints for the Tympaki aquifer Scenario 1

Scenario 2

Scenario 3

Scenario 4

Status

Distance to RHS (m)

Status

Distance to RHS (m)

Status

Distance to RHS (m)

Status

Distance to RHS (m)

1

Satisfied

0.78

Satisfied

0.47

Satisfied

0.17

Satisfied

0.18

2

Satisfied

0.86

Satisfied

0.52

Satisfied

0.30

Near-binding

6.4 × 10−5

3 4

Satisfied Satisfied

0.91 0.66

Satisfied Satisfied

0.69 0.39

Satisfied Satisfied

0.26 0.10

Near-binding Satisfied

1.9 × 10−3 0.07

5 6 7 8

Satisfied Satisfied Satisfied Satisfied

0.62 1.75 3.08 2.73

Satisfied Satisfied Satisfied Satisfied

0.41 0.68 1.79 1.46

Satisfied Near-binding Satisfied Satisfied

0.18 6.8 × 10−5 0.60 0.58

Satisfied Satisfied

0.02 0.14

Satisfied Satisfied

0.16 0.06

9 10 11

Satisfied Satisfied Satisfied

2.22 0.85 0.20

12 13 14

Near-binding Satisfied Satisfied

1.4 × 10−6 0.14 0.66

Satisfied Satisfied Near-binding Satisfied Satisfied Satisfied

1.26 0.15 2.8 × 10−5 0.29 0.43 0.41

Satisfied Satisfied Near-binding Satisfied Near-binding Satisfied

0.03 0.001 6.4 × 10−5 0.11 1.3 × 10−4 0.39

Satisfied Near-binding Satisfied Satisfied Satisfied Satisfied

0.17 1.2 × 10−4 0.03 0.26 0.36 0.18

determine the corresponding sustainable pumping rates in each case. For this purpose, sequential linearization methods were employed in conjunction with the Ghyben-Herzberg approximation of the saltwater front. For both aquifers, the study of both the current situation and the saltwater front retreat scenarios indicate that the human interference in the phenomenon of saltwater intrusion is great due to the fact that the wells are used for irrigation needs of the areas. The location of the saltwater intrusion front exhibits a large difference during the hydrological year, from wet to dry period. During the wet period, the saltwater intrusion front retracts because the pumping activity is greatly reduced and there is increased rainfall; during the dry period, the saltwater intrusion phenomenon intensifies and the front moves farther inland because all the wells are in operation with high pumping rates. For the Tympaki aquifer, four retreat scenarios of increasing intensity were considered in order to investigate the Table 2 Constraint

possibility of the saltwater intrusion front retreat towards the shoreline. These scenarios seek to maximize pumping activity, while allowing the saltwater intrusion front to retreat toward the shoreline in prespecified locations. In all cases, a large reduction of the total pumping rate is required ranging from 60 to 71 % of the current rate. The number of wells that become inactive ranges from 16 to 19 with sometime noticeable reductions in other wells as well. In all scenarios, although the retreat of the saltwater intrusion zone is not very large, the changes in pumping activity are quite large. This indicates the intensity of the saltwater intrusion problem in the area and the overpumping that occurs in order to meet the needs for drinking water supply and irrigation. For the Hersonissos aquifer, similar results have been produced for the three retreat scenarios that were considered. A very wide range of total reduction in pumping activity is observed ranging from 23 to 87 % of the current rate. The

Sensitivity analysis on optimization constraints for the Hersonissos aquifer Scenario 1

Scenario 2

Scenario 3

Status

Distance to RHS (m)

Status

Distance to RHS (m)

Status

Distance to RHS (m)

1

Near-binding

3 × 10−7

2 3 4 5 6 7 8

Satisfied Satisfied Satisfied Satisfied Satisfied Satisfied Satisfied

0.42 0.50 0.65 0.66 0.42 0.36 0.13

Near-binding Satisfied Satisfied Satisfied Satisfied Satisfied Satisfied Satisfied

3.7 × 10−10 0.32 0.41 0.55 0.65 0.53 0.39 0.29

Satisfied Satisfied Satisfied Satisfied Satisfied Near-binding Satisfied Satisfied

0.11 0.09 0.04 0.01 0.13 1.9 × 10−6 0.16 0.35

Saltwater Intrusion Control in Cretan Coastal Aquifers

number of wells that become inactive ranges from 1 to 3 out of the total of five wells located in the study area, with noticeable reductions in other wells in most of the scenarios. In this work, it is assumed that as the saltwater front recedes during the wet season due to the seasonal increase of freshwater heads, the saline water is immediately flushed upon retreat. However, in many cases, a delay in the flushing process is expected, which cannot be captured using the sharp interface assumption. If this delay is very long, then it would lead to a larger saltwater intruded zone than the one estimated using the sharp interface method. The above studies indicate that the saltwater intrusion phenomenon is very difficult to reverse, at least not within a short period of time. In order for the saltwater intrusion front to retract closer to the shoreline large reductions in pumping activity are required in most cases. Thus, alternative water sources need to be ensured if such strategies are to be implemented, in order to cover the needs of the local communities in drinking and irrigation water. The main alternative water supplies that could be considered in the region include seawater and brackish groundwater that will be treated in desalination plants, the construction of hydraulic structures such as reservoirs that would capture the precipitation water and the use of reclaimed water through infiltration ponds, which will also alleviate the saltwater intrusion problem in the aquifer. Although these nontraditional water supply sources are frequently more expensive to develop and operate than traditional sources, the technology associated with them is advancing rapidly, and they could become viable alternatives in cases where the freshwater reserves are not sufficient. From the above, it becomes evident that it is much more cost-effective to focus on the prevention of the phenomenon to the extent possible, rather than its mitigation.

6.

7.

8.

9.

10.

11.

12.

13.

14.

15.

16.

17. 18.

References 19. 1.

2.

3.

4.

5.

Werner, A. D., Bakker, M., Post, V. E. A., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., Simmons, C. T., & Barry, A. (2013). Seawater intrusion processes, investigation, and management: recent advances and future challenges. Advances in Water Resources, 51, 3–26. Karatzas, G. P., & Dokou, Z. (2015). Optimal management of saltwater intrusion in the coastal aquifer of Malia, Crete (Greece). Hydrogeology Journal, 23(6), 1181–1194. Pinder, G. F., & Cooper, H. H. (1970). A numerical technique for calculating transient position of saltwater front. Water Resources Research, 6(3), 875–882. Milnes, E., & Renard, P. (2004). The problem of salt recycling and seawater intrusion in coastal irrigated plains: an example from the Kiti aquifer (southern Cyprus). Journal of Hydrology, 288(3–4), 327–343. Giambastiani, B. M. S., Antonellini, M., Essink, G. H. P. O., & Stuurman, R. J. (2007). Saltwater intrusion in the unconfined

20.

21.

22.

23.

24.

coastal aquifer of Ravenna (Italy): a numerical model. Journal of Hydrology, 340(1–2), 91–104. Kopsiaftis, G., Mantoglou, A., & 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. (2009). Pumping optimization of coastal aquifers based on evolutionary algorithms and surrogate modular neural network models. Advances in Water Resources, 32(4), 507–521. Dokou, Z., & 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. Reilly, T. E., & Goodman, A. S. (1985). Quantitative-analysis of saltwater fresh-water relationships in groundwater systems-a historical-perspective. Journal of Hydrology, 80(1–2), 125–160. Mantoglou, A., Papantoniou, M., & Giannoulopoulos, P. (2004). Management of coastal aquifers based on nonlinear optimization and evolutionary algorithms. Journal of Hydrology, 297(1–4), 209– 228. Mantoglou, A. (2003). Pumping management of coastal aquifers using analytical models of saltwater intrusion. Water Resources Research. doi:10.1029/2002WR001891. Guvanasen, V., Wade, S. C., & Barcelo, M. D. (2000). Simulation of regional ground water flow and salt water intrusion in Hernando County, Florida. Ground Water, 38(5), 772–783. Ababou, R., & Al-Bitar, A. (2004). Salt water intrusion with heterogeneity and uncertainty: mathematical modelling and analyses. Developments in Water Science, 55, 1559–1571. Karterakis, S. M., Karatzas, G. P., Nikolos, I. K., & Papadopoulou, M. P. (2007). Application of linear programming and differential evolutionary optimization methodologies for the solution of coastal subsurface water management problems subject to environmental criteria. Journal of Hydrology, 342(3–4), 270–282. Koukadaki, M. A., Karatzas, G. P., Papadopoulou, M. P., & Vafidis, A. (2007). Identification of the saline zone in a coastal aquifer using electrical tomography data and simulation. Water Resources Management, 21(11), 1881–1898. Papadopoulou, M. P., Varouchakis, E. A., & Karatzas, G. P. (2010). Terrain discontinuity effects in the regional flow of a complex karstified aquifer. Environmental Modeling and Assessment, 15(5), 319–328. Shamir, U., Bear, J., & Gamliel, A. (1984). Optimal annual operation of a coastal aquifer. Water Resources Research, 20, 435–444. Uddameri, V., & Kuchanur, Z. M. (2007). Simulation-optimization approach to assess groundwater availability in Refugio County, TX. Environment Geology, 51, 921–929. Gorelick, S. M., Voss, C. I., Gill, P. E., Murray, W., Saunders, M. A., & Wright, M. H. (1984). Aquifer reclamation design: the use of contaminant transport simulation combined with non-linear programming. Water Resources Research, 20(4), 415–427. Willis, R., & Finney, B. A. (1988). Planning model for optimal control of saltwater intrusion. Journal of Water Resources Planning and Management, 114(2), 163–178. Finney, B. A., Samsuhadi, & Willis, R. (1992). Quasi-threedimensional optimization model of Jakarta Basin. Journal of Water Resources Planning and Management, 118(1), 18–31. Gharbi, A., & Peralta, R. C. (1994). Integrated embedding optimization applied to salt Lake valley aquifers. Water Resources Research, 30(3), 817–832. Reinelt, P. (2005). Seawater intrusion policy analysis with a numerical spatially heterogeneous dynamic optimization model. Water Resources Research, 41(W05006). Cheng, A. H.-D., Halhal, D., Naji, A., & Ouazar, D. (2000). Pumping optimization in saltwater-intruded coastal aquifers. Water Resources Research, 36(8), 2155–2165.

Dokou Z. et al. 25.

Cai, X., McKinney, D. C., & Lasdon, L. S. (2001). Solving nonlinear water management models using a combined genetic algorithm and linear programming approach. Advances in Water Resources, 24, 667–676. 26. Dhar, A., & Datta, B. (2009). Saltwater intrusion management of coastal aquifers. I: linked simulation-optimization. Journal of Hydrologic Engineering, 14(12), 1263–1272. 27. Mantoglou, A., & Papantoniou, M. (2008). Optimal design of pumping networks in coastal aquifers using sharp interface models. Journal of Hydrology, 361, 52–63. 28. Kourakos, G., & Mantoglou, A. (2011). Simulation and multiobjective management of coastal aquifers in semi-arid regions. Water Resources Management, 25, 1063–1074. 29. Papadopoulou, M. P., Varouchakis, E. A., & Karatzas, G. P. (2009). Simulation of complex aquifer behavior using numerical and geostatistical methodologies. Desalination, 237, 42–53.

30.

Varouchakis, E. A., Karatzas, G. P., & Giannopoulos, G. P. (2014). Impact of irrigation scenarios and precipitation projections on the groundwater resources of Viannos basin in the island of Crete. Greece. Environmental Earth Sciences. doi:10.1007/s12665-0143913-2. 31. Llopis-Albert, C., & Pulido-Velazquez, D. (2014). Discussion about the validity of sharp-interface models to deal with seawater intrusion in coastal aquifers. Hydrological Processes, 28(10), 3642–3654. 32. Ahlfeld, D.P., Barlow P.M., & Mulligan, A.E. (2005). GWM—a groundwater-management process for the U.S. Geological Survey modular ground-water model (MODFLOW-2000). USA: U.S. Geological Survey Open-File Report 2005–1072. 33. Paritsis, S.N. (2005). Simulation of seawater intrusion into the Tymbaki aquifer, South Central Crete, Greece. Heraklion, Crete, Greece: Department of Management of Water Resources of the Region of Crete.