Spring protection and sustainable management of groundwater resources in a spring field

Journal of Hydrology 582 (2020) 124498

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Spring protection and sustainable management of groundwater resources in a spring field

T

⁎

Qiankun Luoa, , Yun Yangb,c, Jiazhong Qiana, Xiuxuan Wanga, Xing Changa, Lei Maa, Fulin Lid, Jianfeng Wub a

School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China Department of Hydrosciences, School of Earth Sciences and Engineering, Nanjing University, Nanjing 210046, China c Huaihe River Commission of the Ministry of Water Resources, Bengbu 233001, China d Water Conservancy Research Institute of Shandong Province, Jinan 250013, China b

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by G. Syme, Editor-in-Chief, with the assistance of Li He, Associate Editor

Springs are the natural discharge points of groundwater. They are of considerable value for drinking water supply and some springs are of historical and tourism value. However, increasing industrial and urban expansion has led to serious problems of overdraft of groundwater resources, which in turn had caused the disappearance of springs globally. Thus, the protection of springs while at the same time allowing for the use of local groundwater resources is an important task of water resources management. In this study, we propose a new multi-objective simulation and optimization (S/O) model to find an optimal extraction strategy which balances the competitive relationship between spring outflow and groundwater extraction. In the newly developed model, the extraction rates of different groundwater wells are taken as the decision variables, and maximization of groundwater extraction and maximization of spring outflow are taken as the two objective functions. Baotu Spring which is located in northern China has been famous since the Shang Dynasty (1600 BCE). However, Baotu Spring has been drying up since the 1990s due to over extraction of the local groundwater. Thus, we take Baotu Spring as a case study to demonstrate the applicability of the newly developed model. Simultaneously, a multiobjective evolutionary algorithm (MOEA), the multi-objective fast harmony search algorithm (MOFHS) which is coupled with the commonly used groundwater flow code MODFLOW, is adopted to search the Pareto optimal solutions (the optimal extraction strategies). The optimization results of Baotu Spring field during the management period from July 2013 to June 2014 show that more groundwater resources could be extracted without threatening the sustainable outflow of Baotu Spring and Black Tiger Spring. The optimization results for different hydrological years show that climate variability (mainly in precipitation) is an important factor when choosing the optimal extraction strategy. The optimization results of Baotu Spring in northern China show that the newly developed model is a promising tool to find the optimal groundwater extraction strategies which can protect springs while maximizing the extraction of groundwater resources in a spring field.

Keywords: Groundwater management Spring protection Simulation and optimization Multi-objective optimization MOFHS Baotu Spring

1. Introduction Groundwater is one of the major resources for drinking water supply, irrigation and industry (Richey et al., 2015; Khosraxi et al., 2018). A spring is a natural discharge point of groundwater and plays a vital role in groundwater circulation systems (Elhatip and Günay, 1998; Oraseanu and Mather, 2000; Labat et al, 2002; Birk et al., 2004; Obolewski et al., 2016; Liu et al., 2018). Some springs have played a significant role in human history and cultural development (e.g., Qian et al., 2006; Sato et al., 2015). Generally, a spring field refers to the independent hydrogeological unit where springs are located (e.g., Wang ⁎

et al., 2009; Sun et al., 2014; Gao, 2016). In this hydrogeological unit, the outflow of springs is directly affected by other discharges such as groundwater extraction by wells. Therefore, if the groundwater is overdrafted in a spring field, the outflow of the spring will decrease or even stop completely. With the rapid increase of population and the development of industry, water demand is continuously growing. Especially, in dry areas where groundwater is the main water resource, over extraction of groundwater is common (e.g., Wu and Xu, 2005; Zhang et al., 2017). The Baotu Spring in Jinan, Shandong Province is a very famous tourist attraction in China. Historically Baotu Spring can be dated back

Corresponding author. E-mail address: [email protected] (Q. Luo).

https://doi.org/10.1016/j.jhydrol.2019.124498 Received 30 October 2019; Received in revised form 8 December 2019; Accepted 19 December 2019 Available online 23 December 2019 0022-1694/ © 2019 Published by Elsevier B.V.

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

management of groundwater resources in a spring field is given. Then the numerical simulation model and the multi-objective optimization model of Baotu Spring field are presented in Sections 3 and 4 respectively. An analysis of the optimization results is given in Section 5. At the end, the conclusion of our research is summarized and the future research directions are given.

to the Shang Dynasty, some 3500 years ago (Anonymous, Spring and Autumn Annals). Since then, Baotu Spring has repeatedly been mentioned in the traditional Chinese literature. However, because of the overdraft of groundwater resources for human consumption, Baotu Spring has suffered from frequent drying up since 1972. Between 1999 and 2002, Baotu Spring seized flowing for approximately 926 days (Liu et al., 2018). The local government has taken a series of actions to protect Baotu Spring, including such as importing water from the Yellow and Yangtze rivers and the artificial recharge. The water imported from other watersheds has already become a main source of water supply in Jinan city, covering about 44% of the total water use in 2017 (JWRB and JURWAB, 2017). Especially, in recent years, the groundwater abstraction at Mount La and Dayang village in the Baotu Spring field has been stopped. The question remains whether the exploitation of groundwater in the Baotu Spring field should be stopped completely or could still be allowed within a reasonable rate adapted to varying hydrological conditions. Some researchers have focused on how to ensure the outflow of springs while maximizing the sustainable groundwater extraction in the spring field (e.g., Qian et al., 2006; Kang et al. 2011; Zhang et al., 2017; Liu et al., 2018). However, most of these previous studies only analyzed the relationship between groundwater abstraction and outflow of springs (Zhang et al., 2017; Liu et al., 2018) or considered only one single objective and determined the optimal extraction strategy simply based on comparing a few cases (Wang and Zhang, 1999; Qian et al., 2006; Kang et al. 2011). In order to get a comprehensive series of optimal groundwater extraction strategies while keeping a sustainable outflow of springs under changing environmental conditions, the multi-objective simulation and optimization (S/O) method is adopted in this study. In recent decades, many multi-objective S/O methods have been successfully applied to solve groundwater management problems in many different fields (Reed et al., 2013; Maier et al., 2014; Maier et al., 2019), including optimal design of groundwater remediation system (Erickson et al., 2002; Luo et al., 2012, 2014; Yang et al., 2017; Ouyang et al., 2017), optimal design of a long-term groundwater monitoring network (Reed and Minsker, 2004; Reed et al., 2007; Kollat et al. 2008; Reed et al., 2013), controlling of seawater intrusion in coastal aquifers (Javadi et al., 2015), and optimal management of groundwater resources (Sadeghi-Tabas et al., 2017; Rajabi and Ketabchi, 2017; Alizadeh et al., 2017). The main idea of the multi-objective S/O method is to find satisfactory groundwater extraction strategies by coupling the multi-objective optimization model with the groundwater simulation model. Therefore, there are two major parts in a multi-objective S/O model: the groundwater simulation model and the multi-objective optimization model. Generally, the groundwater simulation model describes the state and character of the groundwater flow, while the multi-objective optimization model with constraints describes the optimization problem. We present a new multi-objective S/O model of groundwater management in a spring field. The main purpose of the newly developed model is to balance the competitive relationship between spring protection and groundwater extraction. Furthermore, we take Baotu Spring field in northern China as a case study to demonstrate the applicability of the newly proposed model. Comparing with the GA-based optimization algorithm, the Multi-objective Fast Harmony Search algorithm (MOFHS) (Luo et al., 2012) does not need the encoding process of the decision variables, which will save a lot of calculation time. Thus the MOFHS is more suitable for solving optimization problems with continuous decision variables. Based on the advantages in finding the Pareto optimal solutions of continuous variable optimization problems (e.g., extraction rates of groundwater wells in this study), the MOFHS is used to find the Pareto optimal solutions (the optimal extraction strategies) of the newly developed multi-objective S/O model. A brief introduction of MOFHS is given in Section 2.3 below. In the following Section 2, an overview over the newly developed multi-objective S/O model of spring protection and sustainable

2. Overview over the multi-objective S/O model of spring protection and sustainable management of groundwater resources in a spring field 2.1. Groundwater flow simulation model Generally, a groundwater flow simulation model is used to calculate the state variables (the groundwater head values in this study) for the multi-objective S/O model of spring protection and sustainable management of groundwater resources in a spring field. In this study, the three-dimensional finite-difference groundwater flow simulator MODFLOW (Harbaugh and McDonald, 1996) is adopted to solve the groundwater flow equation. In order to be conveniently called by the main optimization program, the MODFLOW simulator is changed into subroutine form. 2.2. Multi-objective optimization model of spring protection and sustainable management of groundwater resources in a spring field A multi-objective optimization problem (Rao, 1991) can be mathematically described as:

maximum y = F (x ) = (f1 (x ), f2 (x ), …, fk (x ))

(1)

subject to

gi (x ) ⩽ 0,

i = 1, 2, …, M

(2)

hi (x ) ⩽ 0,

i = 1, 2, …, P

(3)

li ⩽ x i ⩽ ui ,

i = 1, 2, …, N

(4)

where y = (y1 , y2 , …, yk ) ∈ Y , yi = fi(x) represents the ith objective function among k objectives, Y represents the objective function space; x = (x1, x2, …, x n ) ∈ X , x represents a n-dimensional decision variable vector representing a solution, X is the set of all feasible solutions which are constrained by M inequality and P equality constraints; li and ui represent the lower and upper bounds of the ith decision variable (xi), respectively. The objective functions and constraints are usually determined by managers. For instance, minimization of remediation cost and remaining mass of pollutants in aquifers are often taken as the objective functions of multi-objective optimal design of groundwater remediation systems (Singh and Minsker, 2008; Luo et al., 2012; Yang et al., 2013). The constraint conditions include the total number of wells, the minimal values of hydraulic head at special locations, and the maximal solute concentration in certain areas (Singh and Minsker, 2008; Luo et al., 2012; Yang et al., 2013). In this study, maximizing the extraction of groundwater from different groundwater wells while maximizing the outflow of springs are taken as the objective functions. The warning groundwater level at the outlet point of a spring, the maximum extraction ability of groundwater wells, and the basic water demands constitute the constraint conditions. Accordingly, the two objective functions can be mathematically expressed as: u

Maximize F1 =

v

∑ ∑ pi,k tk i=1 k=1

(5)

w

Maximize F2 =

∑ qj j=1

(6)

where F1, F2 are the objective functions for the total amount of 2

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

solutions (the optimal groundwater extraction strategies in this study).

extraction and spring outflow respectively, u is the total number of the groundwater wells, v is the total number of the extraction periods (management periods), w is the total number of springs, pi,k is the extraction amount of the ith groundwater well in the kth period, tk is the length of the kth period, and qj is the total amount of outflow from the jth spring. The constraints can be mathematically formulated as:

3. Numerical simulation of groundwater flow in Baotu Spring field 3.1. Background Baotu Spring field is located in Shandong province, China (see Fig. 2). Based on the previous studies and hydrogeological survey works, the contact zone between the older metamorphic rock complex and younger magmatic intrusion, mainly along the Yellow River and the Xiaoqing River, is considered as the northern boundary of the spring field. The Dongwu fault and the Mashan fault (see Fig. 2), which are considered as water-barriers (Qian et al., 2006; JURWAB, 2008), are taken as the eastern and western boundaries of the spring field. The watershed in the southern mountains is taken as the southern boundary of the spring field. The Baotu Spring field thus covers an area of approximately 1486 km2 (JURWAB, 2008; Wang et al., 2009). From the southern part to the northern part of Baotu Spring field, there are mountains, hilly land, piedmont inclined plain, as well as the Yellow River alluvial plain (Qian et al., 2006; JURWAB, 2008; Wang et al., 2009). The terrain of the study field is high in the south and low in the north, which ranges from 26 m above sea level (m a.s.l.) to 700 m a.s.l (Liu et al., 2018). The climatic condition of Baotu Spring field is a typical warm temperature continental monsoon climate. The annual mean precipitation is approximately 648 mm, and the precipitation is seasonally varying (JWRB and JURWAB, 2017). The wet season lasts from June to September and accounts for 70% of the annual precipitation (JURWAB, 2008). The annual mean potential evaporation of Baotu Spring field is about 1530 mm, which is larger than the annual precipitation (JURWAB, 2008). The relative differences between the precipitation and potential evaporation increase from the southern part to the northern part of the study area (JURWAB, 2008). The Yellow River, Xiaoqing River, Yufu River and Beidasha River are the four major rivers flowing in the study area (see Fig. 2).

hj > hj, L j = 1, ...w pi,min ⩽ pi, k ⩽ pi,max i = 1, ...u, k = 1, ...v u

∑ pi, k ⩾ Qb

k = 1, ...v

i=1

(7)

where hj,L is the warning groundwater level for the jth spring, w is the number of springs, pi,min is the minimum amount of water extracted from the ith groundwater well representing the basic groundwater demand, pi,max is the maximum extraction capacity of the ith groundwater well (determined by its maximum pumping capacity), and Qb is the basic water demand for daily life in the study area. 2.3. Solution by optimization algorithms The multiple objectives in the S/O model are usually conflicting with each other (Dev, 1995). Compared with the best solution of a single objective optimization problem, the Pareto optimal solutions (non-dominated solutions) are usually used to describe the optimal solutions of the multi-objective optimization model (Dev, 1995). A Pareto-optimal solution can be mathematically defined as: if x ∗ ∈ X , and if-and-only-if there is no x ∈ X satisfying F (x ) < F (x ∗) , then x* is a Pareto optimal solution in the decision variable space, X. Generally, in order to find the Pareto optimal solutions, the original optimization model should be transformed into an unconstrained form by adding penalties to the objective functions if they violate the constraint conditions (Reed et al., 2000). In this study, the MOFHS is used, which has demonstrated its superiority in solving multi-objective groundwater management problems (Luo et al., 2012). The main process of the MOFHS-based multi-objective S/O model of sustainable management of groundwater resources in a spring field is shown in Fig. 1 and described as follows: Step 1. Initialization of the MOFHS and the multi-objective optimization problem. The initial setting of the MOFHS algorithm includes setting the control parameters of the MOFHS, the objective functions, the constraint conditions and the stopping criterion of the MOFHS. The initial population of the MOFHS is also set in the first step. The initial setting of the optimization problem includes establishment of the multiobjective optimization model and the groundwater simulation model of the study area (Luo et al., 2012). Step 2. Computation of the state variables (groundwater levels in this study) of the study area by calling the MODFLOW simulator. In this study, the calculated groundwater levels are then used to judge whether the present extraction strategy can satisfy the constraints. Step 3. Calculation of the objective function values (the total amount of groundwater extraction and the total amount of spring flow in the study area) of every individual in the present population. Step 4. Ranking of the objective functions by Pareto ranking and the selection of the best solutions by Elitist selection technologies (Deb et al., 2002; Luo et al., 2012). Step 5. Generalization of the new population according to the best individual based on the harmony memory consideration and pitch adjustment rules (Luo et al., 2012). The number of individuals in the new population is the same as that in the initial population. Step 6. Evaluation of the stopping criterion. Usually the maximally allowed number of the evolutionary generations is taken as the stopping criterion of evolutionary algorithms (Luo et al., 2011; Luo et al., 2012). If the iteration process reaches the termination criterion, the iteration process of MOFHS will stop and lead to the output of the Pareto optimal

3.2. Geological and hydrological conditions of Baotu Spring field The tectonic unit of Baotu Spring field is part of the northern edge of the doming structure in Tai Mountain, the northern wing of the Lu Zhongnan anticline (JURWAB, 2008; Wang et al., 2009). The strata in the field include Archaean, Cambrian, Ordovician and Quaternary sedimentary rocks, and igneous rocks. There are six active faults in the study field including Wenzu, Dongwu, Qianfoshan, Chaomidian, Mashan, and Ganggou faults, which are approximately parallel and aligned with the north-south direction (JURWAB, 2008). The hydrogeological conditions of the Baotu Spring field are complex. The main aquifers in the field are Quaternary porous aquifers and Cambrian-Ordovician fractured karst aquifers. The Quaternary porous aquifers are mainly distributed in the Yufu River and the Beidasha River alluvial fans, and the thicknesses of the aquifers are approximately 10–30 m (Qian et al., 2006; JURWAB, 2008; Kang et al., 2011). The Cambrian-Ordovician fractured karst water is abundant in areas where the fractures and the karst are well-connected. The most important recharge of Baotu Spring field is precipitation. In the southern mountain areas where the Cambrian and Ordovician rocks dominate, the infiltration of precipitation constitutes the main source of groundwater recharge (Qian et al., 2006; JURWAB, 2008; Kang et al., 2011). Other sources of recharge in Baotu Spring field include the leakage of rivers and reservoirs and lateral inflows. The deep fractured karst aquifer is further recharged by the leakage from the upper Quaternary aquifers. The direction of groundwater flow is mainly from southeast to northwest and controlled by terrain and faults. When the groundwater flow arrives at the northern part of the field, the flow is upward due to the blockage by the igneous rocks in the north. Finally, the upward flow reaches the surface and generates the famous Jinan spring series such as 3

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

Fig. 1. Flowchart of the MOFHS-based multi-objective S/O model of spring protection and sustainable management of groundwater in spring field.

et al. (2011). Previous studies revealed that the fractured karst in Baotu Spring field is well developed and has a uniform water table (Qian et al., 2006; Kang et al., 2011). Therefore, in this study, the equivalent porous medium model is used as the hydrogeological conceptual model of Baotu Spring field and a three dimensional transient flow model is used to simulate the groundwater flow in the study area. The mathematical model of groundwater flow in the study area can be stated as:

the Baotu Spring and the Black Tiger Spring (the direction of groundwater flow and the location of Baotu Spring are shown in Fig. 2). Under natural conditions, spring flow is the main discharge of the local groundwater system. However, because of the rapid urbanization of Jinan, extraction of groundwater is nowadays the main discharge process, which has led to a dramatic decline and even drying up of the spring flow. Phreatic evaporation in regions of shallow groundwater tables is another way of discharge for the groundwater system.

∂H

⎧∇ (K ∇H ) + W = μs ∂t x , y, z ∈ Ω , t ⩾ 0 ⎪ H (x , y, z , 0)| = H0 x , y, z ∈ Ω ⎪

3.3. Hydrogeological conceptual model and mathematical model of the groundwater flow in Baotu Spring field

⎨ ⎪ ⎪ ⎩

The domain of the hydrogeological conceptual model of groundwater flow in Baotu Spring field is congruent with the study area. Therefore, the eastern and western boundaries of the conceptual model are Dongwu and Mashan faults respectively, the southern boundary is the lowest part of the middle Cambrian Zhangxia Formation, and the northern boundary is on the top of the limestone bed located 600 m below the ground (JURWAB, 2008). Based on the hydrogeological conditions of Baotu Spring field, the conceptual model includes three layers. The first layer is an unconfined aquifer consisting of the Holocene and Late Pleistocene loose sediments. The second layer is an aquitard consisting of Quaternary clay. The distribution of the second layer is congruent with the first layer. The third layer consists of the confined-unconfined fractured karst aquifers of Cambrian and Ordovician ages (Qian et al., 2006; JURWAB, 2008). The boundary conditions for the three layers are set as shown in Fig. 3 (JURWAB, 2008). The zonation of the hydrogeological parameters is also shown in Fig. 3 (JURWAB, 2008; Kang et al., 2011). The recharge infiltration ratios are set according to the JURWAB (2008) and Kang

∂H

Kn ∂n hr − H σ

=

Γ1

= 0 x , y, z ∈ Γ1, t ⩾ 0

∂H Kn ∂n Γ2

x , y, z ∈ Γ2, t ⩾ 0

(8)

where K represents the hydraulic conductivity, t represents time, H (x,y,z,t) represents the hydraulic head at time t, H0(x,y,z) represents the initial head, μs represents the specific storage coefficient, and W represents the source and sink items. The value of W is the combination of precipitation, recharge by irrigation, recharge by seepage from reservoirs, evapotranspiration from groundwater, and the well discharge from the aquifer in this study. Ω represents the scale of the study area, Γ1 and Γ3 represent impervious and third-type boundaries, respectively, Kn represents the hydraulic conductivity in the direction normal to the boundary surface, σ represents the thickness of the third-type boundary, and hr represents the specified head value of the third-type boundary. 3.4. Model identification and verification In this study, the mathematical model of groundwater flow in Baotu 4

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

Fig. 2. (a) The location and hydrogeological map of the Baotu Spring field in Shandong Province, China; (b) Geological profile of the Jinan Spring field. 1—quaternary (Q), 2—Igneous rocks (γ), 3—Ordovician (O), 4—Cambrian (Є), 5—Archean (Ar), 6—Karst groundwater recharge area (Є-O), 7—Karst groundwater abundant and discharge area (Є-O), 8—Karst groundwater buried below impermeable sandstone and shale (C-P), 9—Impermeable igneous rocks, 10—Fault, 11—Spring, 12—Watershed, 13—Water source, 14—Flow direction, 15—Groundwater observation point, 16—River. (Modified from Qian et al., 2006; Kang et al., 2011).

Fig. 3. Boundary condition and zonation of the hydrogeological parameters (after Qian et al., 2006; JURWAB, 2008; Kang et al., 2011). (a) Hydraulic conductivity and the specific yield (or specific storage) of the first and second layers; (b) Hydraulic conductivity and the specific yield (or specific storage) of the third layer. 5

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

Spring field was numerically solved by MODFLOW and the hydrogeological parameters (the hydraulic conductivity, the specific yield or specific storage of aquifers) of the simulation model were identified by the fast harmony search algorithm (FHS) as in Luo et al. (2011). The identification period is from 1 October 2012 to 30 September 2013, and the verification period is from 1 October 2013 to 30 September 2014. The observation values of groundwater head were provided by the Water Conservancy Research Institute of Shandong Province (unpublished data). Different from other parts of the Baotu Spring field, the temporal variation of groundwater head values in the southern mountain area is significantly larger than that in the northern plain (the maximum variation of observed groundwater levels in Baotu Spring field was about 95 m in 2014). Therefore, the traditional method of calibration, which uses the root mean square error (RMSE) as the objective function to identify hydrogeological parameters, would lead to unsatisfactory results. After several attempts, minimization of the sum of square of relative errors (SSRE) was taken as the objective function during the model identification stage. The mathematical form of the model identification objective in this study is as follows:

obj = min

n m i, j i, j hcal − hobs ⎛ ⎛ ⎞⎞ ⎟ i i ⎜∑ ∑ ⎜ max(hobs ⎟ ) − min( h ) obs ⎠⎠ ⎝ i=1 j=1 ⎝

Fig. 4 shows the observed and computed values of groundwater levels at three observation points (Baotu Spring, Black Tiger Spring and J65; see Fig. 2). The computed values fit the observed water levels well (Fig. 4). In the identification period, the maximum absolute errors between the computed and observed values are 0.34 m, 0.36 m, and 0.31 m at Baotu Spring, Black Tiger Spring and J65 respectively. In the verification period, the maximum absolute errors between the computed and observed values are 0.36 m, 0.37 m, and 0.39 m at Baotu Spring, Black Tiger Spring and J65 respectively. Thus, the identified hydrogeological parameters seem reasonable and can be used as the input parameters of the numerical simulation model of groundwater flow in Baotu Spring field. 4. Multi-objective optimization model of spring protection and sustainable management of groundwater resources in Baotu Spring field The multi-objective optimization model of spring protection and sustainable management of groundwater resources in Baotu Spring field is the same as described in Eqs. (5)–(7). The optimization model aims at maximizing both the groundwater extraction and spring outflow in Baotu Spring field. Therefore, the first objective of the optimization model is to maximize the total amount of groundwater extracted by the groundwater wells including Qiaozili, Cold village, Ancient city, Mount Emei, Dayang village, and Mount La during the management period. The location of the six groundwater wells is shown in Fig. 2. The second objective of the optimization model is to maximize the total amount of the spring outflow including Baotu and Black Tiger springs. The constraints include the warning water level (27.5 m a.s.l. at the outlet point of Baotu Spring), the maximum extraction capacity of every groundwater well, and the basic water demand Qb which is set to 16 × 104m3/d in our optimization process. The decision variables are the amounts of groundwater extracted from different groundwater wells. Generally, the hydrological year in Jinan is from July to the following June, with the wet season from July to September, and the dry season from March to June (JURWAB, 2008). Thus, the optimization model in this study includes three management periods: the first management period is from July to October, the second management period is from November to the following February, and the third management period is from March to June. The period from July 2013 to June 2014 was taken as the total management period in this study. Thus there are 18 decision variables (pumping rates of six groundwater wells in each of three different management periods) which correspond to pi,k in Eq. (7). The maximum and minimum values of pi,k are presented in Table 1.

(9)

where, i denotes the ith observation point, j denotes the jth observation period, n is the total number of observation points, m is the total i, j is the simulated head value at the number of observation periods, hcal i, j ith point in the jth period, hobs is the observed head value at the ith point i, j i, j ) and min(hobs ) are the maximum and the in the jth period, and max(hobs minimum head values at the ith point over the whole observation period.

5. Optimization results and the sustainable groundwater management strategies Based on the multi-objective S/O framework presented in Fig. 1, the established multi-objective optimization model, coupled with the numerical simulation model of groundwater flow, was used to find the optimal extracting strategies in Baotu Spring field. The optimization results for the period between July 2013 and June 2014 were used to demonstrate the effectiveness of the new method. The optimization results for different hydrological years can provide the water managers with the best optimal exploitation strategies. 5.1. Optimization results of the management period from July 2013 to June 2014 Fig. 5 shows the Pareto optimal solutions based on the MOFHS during the management period from July 2013 to June 2014. The maximum and minimum values of the total groundwater extraction by wells are 1.00 × 108m3/y and 0.72 × 108m3/y respectively. The maximum and minimum values of spring outflow are 0.51 × 108m3/y

Fig. 4. Comparison of the calculated and observed groundwater head values during the identification and calibration periods. (a) Baotu Spring; (b) Black tiger Spring; (c) J65. 6

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

Table 1 The maximum and minimum extractions from different water wells Unit: ×104m3/d. Water resources

Qiaozili

Cold village

Ancient city

Mount Emei

Dayang village

Mount La

Maximum value Minimum value

8 0

4 0

8 0

5 0

10 0

7 0

Fig. 6. Statistics of extraction from different groundwater wells for all of the Pareto optimal solutions, arranged in ascending order by the total extraction.

optimal solutions (green points shown in Fig. 5) are given in Fig. 6, rearranged in ascending order by the total extraction amounts. It is clear that the extraction of Qiaozili, Cold village, Ancient city, and Mount Emei almost reach the maximum pumping capacity for all of the Pareto optimal solutions, while the situations of Dayang village and Mount La are different. The extractions of Dayang village and Mount La increase parallel to the increasing total extraction, which indicates that the differences among different Pareto optimal solutions are mainly due to the extractions of Dayang village and Mount La. That is to say, Dayang village and Mount La are the most sensitive groundwater wells for the Pareto optimal solutions and these two wells have a decisive influence on the outflow of the springs. In order to further study the relationship of extraction during different management periods, the distribution of the total extraction during different management periods of the Pareto optimal solutions (green points shown in Fig. 5) is given in Fig. 7. It shows that the total extraction in the first management period is small. This is because that the first management period is a groundwater level recovery period with low initial groundwater levels in the study area. Although the precipitation is abundant, the risk of springs drying up would increase if the extraction was increased in this management period. Therefore, the total extraction in the first management period is the lowest. The total extraction in the second period lies between the values of the first and the third management periods. The reason is that there is a certain amount of groundwater storage after the first management period. Thus, the total extraction in the second management period can be larger than during the first management period. Considering that the groundwater level in the third management period should satisfy the constraint conditions (above the warning level), the total extraction is lower than during the third management period. Comparatively, the total extraction in the third management period is mainly located in the top region of Fig. 7. This is because that there is not only extra recharge (e.g., artificial recharge) for the groundwater system during the third

Fig. 5. Pareto optimal solutions based on the MOFHS algorithm between July 2013 and June 2014.

and 0.45 × 108m3/y respectively. Every green point in Fig. 5 represents a Pareto optimal solution (an optimal groundwater extraction strategy) which satisfies all the constraint conditions. The water managers can choose one of them as the preferred strategy. For instance, the Pareto solution shown as a red star can be selected as the extraction strategy between July 2013 and June 2014. Table 2 shows the extraction distribution of the six groundwater wells during the three management periods for the selected solution. Compared with the actual extraction in the study field (about 0.29 × 108m3 between July 2013 and June 2014), the total exploitation of the selected solution (about 0.90 × 108m3) is much bigger. That is to say, more of the groundwater resource can be exploited without the danger of the springs drying up during the three management periods between July 2013 and June 2014. From Table 2 it can also be summarized that the extraction is mainly concentrated in Qiaozili, Cold village, Ancient city, and Mount Emei, whereas Dayang village only provides a little groundwater during the third management period. The exploitation of Mount La is equal to zero during all three management periods. The reason for the small extraction of Dayang village and zero extraction of Mount La is because of their close hydraulic connection with the springs, which is consistent with previous studies (JURWAB, 2008). In this situation, the imported water can be delivered to Dayang village and Mount La by the public water pipe network to satisfy the basic demand of daily life. In order to further analyze the relationship between groundwater extraction of different groundwater wells and outflow of springs, the statistics of extraction from different groundwater wells for all Pareto

Table 2 The distribution of the extraction from different groundwater wells of the preselected Pareto optimal solution Unit: ×104m3/d. Management period

Qiaozili

Cold village

Ancient city

Mount Emei

Dayang village

Mount La

Total amount

The first management period The second management period The third management period

8.00 8.00 8.00

3.94 4.00 4.00

8.00 8.00 8.00

0.78 5 5

0.00 0.00 4.24

0.00 0.00 0.00

20.72 25.00 29.24

7

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

to high flow year. That is because precipitation is the major recharge of the aquifer in the study field. Fig. 8 also indicates that the ranges of the objective functions are increasing from low flow years to high flow years. That is to say, there is a wider range of extracting strategies, from which water managers can select in high flow years compared with normal and low flow years. The maximum and minimum values of the objective functions under different hydrological years can be found in Table 4. Generally, the extraction rates of different groundwater wells are varying in different hydrogeological years. Therefore, the outflow of springs may not always be the maximum under a selected optimal extraction strategy. Taking the total extraction equaling to 0.76 × 108m3/ y as an example (the black line in Fig. 8), the detailed exploitation of different groundwater wells in different hydrological years is shown in Table 5. It can be seen that the extraction in the normal flow year is almost equal to the one in the high flow year. However, the situation in the low flow year is different, and the exploitation at Dayang village and Mount La is nearly 0. In this case, the exploitation of Mount Emei must be increased to ensure the maximization of the total exploitation. In normal and high flow years, Dayang village can provide a certain amount of groundwater during the third management period. Based on the above analysis, the hydrological situation is an important factor in selecting the optimal extraction strategy. For the purpose of spring protection, water managers can select the best extraction strategy from the overlapping part (the pink part) of Fig. 8 for different hydrological years. In the overlapping area, every Pareto optimal solution satisfies the constraint conditions under different hydrological years. To further demonstrate the reliability of the Pareto optimal solutions in the pink area of Fig. 8, the extraction strategies for a low flow year in the overlapping area are recalculated for a high and normal flow year respectively. The red and blue triangles in Fig. 8 represent the calculation results. The results show that with the increase of extraction, the gaps between the recalculated results and the high and normal flow year of spring outflow are increasing. This means the spring outflow may not always be the maximum for different hydrological years. Thus, to be safe, the extraction objective equaling to 0.71 × 108m3/y (the red star in Fig. 8, which has the smallest gap among the outflow of springs) can be selected as the best extraction strategy in Baotu Spring field. The detailed distribution of extraction of different groundwater wells is shown in Table 6. However, it needs to be emphasized that the initial conditions for the optimization process play a vital role in the optimal extraction strategies. The optimization results under different hydrological years

Fig. 7. Distribution of the total extraction in different management periods, arranged in ascending order by the total extraction.

Table 3 Precipitation in different typical hydrological years in Jinan Unit: mm. Frequency

Low flow year

Normal flow year

High flow year

Precipitation

522.8

622.8

764.4

period but also more precipitation than in the second management period. Therefore, the extraction rate in the third management period is the largest among three management periods. The optimization results of the management period from July 2013 to June 2014 demonstrate that the proposed methodology can find the optimal groundwater extraction strategies in Baotu Spring field while guaranteeing the sustainable outflow of the springs. Furthermore, comparing with the actual groundwater exploitation during the optimization period, there is more groundwater available than is extracted under the preselected extraction strategy. 5.2. Optimization results under different hydrological years In order to find the Pareto optimal solutions, which satisfy the conditions for different hydrological years, the proposed multi-objective S/O model is rerun for different hydrological years. Table 3 shows the precipitation of different typical hydrological years based on the precipitation data from 1960 to 2014 in Baotu Spring field. The optimization results are shown in Fig. 8. It can be seen that at the same extraction level, the outflow of springs is increasing from low flow year

Fig. 8. Pareto optimal solutions under different hydrological years. The red star represents the final selected Pareto optimal solution. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 8

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

Table 4 The range of the objective functions under different hydrological years Unit: ×108m3/y. Objective functions The total exploitation of groundwater

The The The The

The total of spring outflow

maximum value minimum value maximum value minimum value

Low flow year

Normal flow year

High flow year

0.79 0.62 0.45 0.42

0.96 0.71 0.47 0.42

1.10 0.71 0.50 0.43

Table 5 The distribution of extraction under different hydrological years when the total exploitation of groundwater is equal to 0.76 × 108m3/y Unit: ×104m3/d. Management period

Hydrological year

Qiaozili

Cold village

Ancient city

Mount Emei

Dayang village

Mount La

The first management period

Low flow year Normal flow year High flow year Low flow year Normal flow year High flow year Low flow year Normal flow year High flow year

7.71 8.00 8.00 8.00 8.00 8.00 8.00 8.00 8.00

4.00 4.00 3.96 4.00 3.98 4.00 3.90 3.78 3.80

6.58 6.32 6.18 6.88 6.34 6.29 7.42 7.30 7.23

1.37 0.82 0.90 0.92 0.73 0.77 4.45 4.85 4.98

0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.27 1.20

0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

The second management period

The third management period

Table 6 The distribution of extraction for the final selected Pareto optimal solution Unit: ×104m3/d. Management period

Qiaozili

Cold village

Ancient city

Mount Emei

Dayang village

Mount La

Total amount

The first management period The second management period The third management period

7.72 7.99 7.99

4.00 4.00 4.00

6.12 6.04 6.67

0.55 0.96 3.14

0.00 0.00 0.00

0.00 0.00 0.00

18.39 18.99 21.80

the study areas. Besides, we did not consider the strong spatial variability of the fractured karst aquifer system in Baotu Spring field in our study. In order to get more accurate predictions of the groundwater system, the spatial variability of the fractured karst aquifers should be considered in the future studies. Accordingly, the computational burden of the multi-objective S/O model for the field application will increase dramatically when considering the parameter uncertainty related to aquifer heterogeneity by stochastic methods. New efficient calculation methods such as parallel computing and surrogate model methods could be introduced. In our study, we only optimized the compromise relationship between groundwater extraction and the outflow of springs in a spring field. Efficient reallocation of the extracted groundwater resources combined with the imported surface water should be studied in the future.

in our study are based on the hydrological conditions in June 2013. If we want to predict the optimal extraction strategy in the future, the initial hydrological conditions must be determined first. Then, the Pareto optimal solution at the lower left corner of the overlapping part of the Pareto front under different hydrological years (i.e., the red star in Fig. 8) can be selected as the final extraction strategy with the lowest risk. 6. Conclusion and future research In this study, we focused on the solution of the problem of drying-up of springs with the help of optimal design of the local groundwater extraction strategy. Based on the framework of a multi-objective S/O method, we established a new multi-objective S/O model in a spring field to find the optimal groundwater extraction strategies. The main purpose is to maximize the groundwater extraction while protecting the springs. The MOFHS algorithm coupled with the commonly used groundwater flow simulator MODFLOW was used to find the Pareto optimal solutions. The optimization results for the three management periods between July 2013 and June 2014 demonstrate that the newly developed model can provide decision makers with a series of Pareto optimal solutions, which satisfy all constraints. The optimization results for different hydrologic years indicate that climatic variability (precipitation) is an important influencing factor when choosing the best extraction strategy. Water managers can select the best solution from the lower left corner of the overlapping part of the Pareto front under different hydrological years for Baotu Spring field. It is important to note that the optimization results in this study are based on the limited data collected by our research team. However, the proposed multi-objective S/O model of groundwater management in a spring field is an open system. If there are new data in the future, the model can be rerun to find more precise Pareto optimal solutions. In addition, if the proposed model is used to solve the spring protection problems in other spring fields, only the groundwater flow simulation model and the constraint conditions need to be changed according to

CRediT authorship contribution statement Qiankun Luo: Conceptualization, Methodology, Data curation, Writing - original draft. Yun Yang: Formal analysis, Software, Writing review & editing. Jiazhong Qian: Funding acquisition, Project administration. Xiuxuan Wang: Investigation, Validation. Xing Chang: Investigation. Lei Ma: Investigation. Fulin Li: Resources. Jianfeng Wu: Supervision. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This study is financially supported by the National Natural Science Foundation of China (Nos. 41502226, 41831289, 41602256, 41572242), the Natural Science Foundation of Anhui Province (No. 9

Journal of Hydrology 582 (2020) 124498

Q. Luo, et al.

1708085QD82), the Fundamental Research Funds for the Central Universities of China (No. JZ2018HGTB0251), the Shandong Province Science and Technology Support Pilot Program for Achieving Water Ecological Civilization under Contract (No. SSTWMZCJHSD04). We thank the Water Conservancy Research Institute of Shandong Province for the assistance in collecting hydrogeological data. We also would like to thank Wolfgang Kinzelbach (Zurich) and Xiuyu Liang (Shenzhen) for their helpful suggestions on the manuscript. The authors are grateful to editor Geoff Syme and three anonymous reviewers for their invaluable comments on the manuscript.

algorithm under uncertainty. J. Hydrol. 519, 3305–3315. Maier, H.R., Kapelan, Z., Kasprzyk, J., Kollat, J., Matott, L.S., Cunha, M.C., Dancy, G.C., Gibbs, M.S., Keedwell, E., Marchi, A., Ostfeld, A., Savic, D., Solomatine, D.P., Vrugt, J.A., Zecchin, A.C., Minsker, B.S., Barbour, E.J., Kuczera, G., Pasha, F., Castelletti, A., Giuliani, M., Reed, P.M., 2014. Evolutionary algorithm and other metaheuristics in water resources: current status, research challenges and future directions. Environ. Model. Softw. 62, 271–299. Maier, H.R., Razavi, S., Kapelan, Z., Matott, L.S., Kasprzyk, J., Tolson, B.A., 2019. Introductory overview: optimization using evolutionary algorithms and other metaheuristics. Environ. Model. Softw. 114, 195–213. Obolewski, K., Strzelczak, A., Glińska-Lewczuk, K., Osadowski, Z., Astel, A., Timofte, C.M., 2016. Ecohydrological relationship between benthic communities and environmental conditions in the spring area. Environ. Eng. Manag. J. 15 (6), 1281–1291. Oraseanu, I., Mather, J., 2000. Karst hydrogeology and origin of thermal waters in the Codru Moma Mountains, Romainia. Hydrog. J. 8 (4), 379–389. Ouyang, Q., Lu, W.X., Hou, Z.Y., Zhang, Y., Li, S., Luo, J.N., 2017. Chance-constrained multi-objective optimization of groundwater remediation design at DNAPLs-contaminated sites using a multi-algorithm genetically adaptive method. J. Contam. Hydrol. 200, 15–23. Qian, J.Z., Zhan, H.B., Wu, Y.F., Li, F.L., Wang, J.Q., 2006. Fractured-karst spring-flow protections: a case study in Jinan. China. Hydroge. J. 14 (7), 1192–1205. Rajabi, M.M., Ketabchi, H., 2017. Uncertainty-based simulation-optimization using Gaussian process emulation: application to coastal groundwater management. J. Hydro. 555, 518–534. Rao, S.S., 1991. Optimization Theory and Application. Wiley Eastern Limited, New Delhi. Reed, P.M., Minsker, B.S., 2004. Striking the balance: long-term groundwater monitoring design for conflicting objectives. J. Water Resour. Plann. Manage. 130 (2), 140–149. Reed, P.M., Minsker, B.S., Valocchi, A.J., 2000. Cost-effective long-term groundwater monitoring design using a genetic algorithm and global mass interpolation. Water Resour. Res. 36 (12), 3731–3741. Reed, P.M., Kollat, J.B., Devireddy, V.K., 2007. Using interactive archives in evolutionary multi-objective optimization: a case study for long-term groundwater monitoring design. Enviro. Model. Softw. 22, 683–692. Reed, P.M., Hadaka, D., Herman, J.D., Kasprzyk, J.R., Kollat, J.B., 2013. Evolutionary multi-objective optimization in water resources: the past, present, and future. Adv. Water Resour. 51, 438–456. Richey, A.S., Thomas, B.F., Lo, M.H., Reager, J.T., Famiglietti, J., Voss, K., Swenson, S., Rodell, M., 2015. Quantifying renewable groundwater stress with GRACE. Water Resour. Res. 51, 5217–5238. Sadeghi-Tabas, S., Samadi, S.Z., Akbarpour, A., Pourreza-Bilonbi, M., 2017. Sustainable groundwater modeling using single- and multi-objective optimization algorithms. J. Hydroinform. 19 (1), 97–114. Sato, M., Ohara, K., Ono, K., Mori, Y., Osato, K., Okabe, T., Nakata, H., Yanagisawa, N., Muraoka, H., 2015. Development of Micro Grid Kalina Cycle System – The First Demonstration Plant in Hot Spring Area in Japan. Proceedings World Geothermal Congress 2015, Melbourne, Australia, 19–25 April. Singh, A., Minsker, B., 2008. Uncertainty-based multi-objective optimization of groundwater remediation design. Water Resour. Res. 44, W02404. Sun, B., Peng, Y.M., 2014. Boundary condition, water cycle and water environment changes in the Jinan spring region. Carsol. Ainica 33 (3), 271–279 (in Chinese with English abstract). Wang, G.L., Zhang, W.G., 1999. Research on management model for Jinan fissure-karst water system. J. Dalian Univ. Technol. 39 (5), 689–693 (in Chinese with English abstract). Wang, Q.B., Duan, X.M., Gao, Z.D., Fred, K.L., Zhou, Y.X., Li, W.P., 2009. Groundwater flow modelling in the Jinan karst spring area. Hydrog. Eng. Geol. 5, 53–60 (in Chinese with English abstract). Wu, Q., Xu, H., 2005. A three-dimensional model and its potential application to spring protection. Environ. Geol. 48, 551–558. Yang, Y., Wu, J.F., Luo, Q.K., Zhang, T., Wu, J.C., Wang, J.G., 2017. Effects of stochastic simulations on multiobjective optimization of groundwater remediation design under uncertainty. J. Hydrol. Eng. 22 (8), 1084–10699. Yang, Y., Wu, J.F., Sun, X.M., Wu, J.C., Zheng, C.M., 2013. A niched Pareto tabu search for multi-objective optimal design of groundwater remediation systems. J. Hydro. 490, 56–73. Zhang, Z.X., Liu, Y., Zhang, F.X., Zhang, L.X., 2017. Study on dynamic relationship of spring water in Jinan spring area based on gray relational analysis. IOP Conf. Series: Earth Environ. Sci. 128, 012068.

Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.jhydrol.2019.124498. References Alizadeh, M.R., Nikoo, M.R., Rakhshandehroo, G.R., 2017. Hydro-environmental management of groundwater resources: a fuzzy-based multi-objective compromise approach. J. Hydro. 551, 540–554. Birk, S., Liedl, R., Sauter, M., 2004. Identification of localised recharge and conduit flow by combined analysis of hydraulic and physico-chemical spring responses (Urenbrunnen, SW-Germany). J. Hydro. 286, 179–193. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multi-objective genetic algorithm: NSGA-II. IEEE Trans. Evolut. Comput. 6 (2), 182–197. Dev, K., 1995. Optimization for Engineering Design: Algorithms and Examples. PrenticeHall, New Delhi. Erickson, M., Mayer, A., Horn, J., 2002. Multi-objective optimal design of groundwater remediation systems: application of the niched Pareto genetic algorithm (NPGA). Adv. Water Resour. 25 (1), 51–65. Elhatip, H., Günay, G., 1998. Karst hydrogeology of the Kas-Kalkan springs along the Mediterranean coast of Turkey. Environ. Geol. 36 (1–2), 150–158. Gao, X.B., 2016. Geochemical Evolution of Niangziguan Karst Water System in Shanxi Province, Northern China. China University of Geosciences Press, Beijing. Harbaugh, A.W., McDonald, M.G, 1996. Programmer's documentation for MODFLOW-96, an update to the U.S. Geological Survey modular finite-difference ground-water flow model: U.S. Geological Survey Open-File Report, pp. 96–486. Javadi, A.A., Hussain, M.S., Sherif, M., Farmani, R., 2015. Multi-objective optimization of different management scenarios to control seawater intrusion in coastal aquifers. Water Resour. Manage. 29, 1843–1857. JURWAB (Jinan Urban and Rural Water Affairs Bureau), 2008. Investigation and Evaluation of Groundwater Resources in Jinan City. JWRB (Jinan Water Resources Bureau) and JURWAB (Jinan Urban and Rural Water Affairs Bureau), 2017. Water Resources Bulletin of Jinan. Kang, F.X., Jin, M.G., Qin, P.R., 2011. Sustainable yield of a karst aquifer system: a case study of Jinan springs in northern China. Hydrogeo. J. 19 (4), 851–863. Khosraxi, K., Panahi, M., Bui, T., 2018. Spatial prediction of groundwater spring potential mapping based on an adaptive neuro-fuzzy inference system and metaheuristic optimization. Hydrol. Earth Syst. Sci. 22, 4771–4792. Kollat, J.B., Reed, P.M., Kasprzyk, J.R., 2008. A new epsilon-dominance hierarchical bayesian optimization algorithm for large multiobjective monitoring network design problems. Adv. Water Resour. 31, 828–845. Labat, D., Mangin, A., Ababou, R., 2002. Rainfall-runoff relations for karstic springs: multifractal analyses. J. Hydro. 256, 176–195. Liu, X.M., Hu, L.T., Sun, K.N., 2018. Analysis of spring flow change in the Jinan city under influences of recent human activities. Proc. IAHS 379, 263–268. Luo, Q.K., Wang, P., Zhu, G.R., 2011. Fast harmony search algorithm and its application to hydrogeological parameters identification (in Chinese with English abstract). Hydrog. Eng. Geol. 38 (4), 15–19. Luo, Q.K., Wu, J.F., Sun, X.M., Yang, Y., Wu, J.C., 2012. Optimal design of a groundwater remediation system using a multi-objective fast harmony search algorithm. Hydrog. J. 20 (8), 1497–1510. Luo, Q.K., Wu, J.F., Yang, Y., Qian, J.Z., Wu, J.C., 2014. Optimal design of groundwater remediation system using a probabilistic multi-objective fast harmony search

10