Applying Genetic Algorithm and Neural Network to the Conjunctive Use of Surface and Subsurface Water

DOI 10.1007/s11269-013-0418-9

Applying Genetic Algorithm and Neural Network

to the Conjunctive Use of Surface

and Subsurface Water

Received: 4 July 2011 / Accepted: 5 August 2013 / Published online: 8 October 2013

© Springer Science+Business Media Dordrecht 2013

Abstract The conjunctive use of surface and subsurface water is one of the most effective ways to increase water supply reliability with minimal cost and environmen-tal impact. This study presents a novel stepwise optimization model for optimizing the conjunctive use of surface and subsurface water resource management. At each time step, the proposed model decomposes the nonlinear conjunctive use problem into a linear surface water allocation sub-problem and a nonlinear groundwater simulation sub-problem. Instead of using a nonlinear algorithm to solve the entire problem, this decomposition approach integrates a linear algorithm with greater computational efficiency. Specifically, this study proposes a hybrid approach con-sisting of Genetic Algorithm (GA), Artificial Neural Network (ANN), and Linear Programming (LP) to solve the decomposed two-level problem. The top level uses GA to determine the optimal pumping rates and link the lower level sub-problem, while LP determines the optimal surface water allocation, and ANN performs the groundwater simulation. Because the optimization computation requires many groundwater simulations, the ANN instead of traditional numerical simulation greatly reduces the computational burden. The high computing performance of both LP and ANN significantly increase the computational efficiency of entire model. This study examines four case studies to determine the supply efficiencies under different operation models. Unlike the high interaction between climate conditions

Y. W. Chen_{· L. C. Chang (}

_B

Department of the Civil Engineering, National Chiao Tung University, 1001 University Rd., Hsinchu 300, Hsinchu, Taiwan, Republic of China e-mail: lcchang31938@gmail.com

C. W. Huang

Department of Bioenvironmental Systems Engineering, 1, Section 4, Roosevelt Rd., Taipei 106, Republic of China H. J. Chu

Department of Geomatics, National Cheng Kung University, 1 University Rd., Tainan 701, Republic of China

and surface water resource, groundwater resources are more stable than the surface water resources for water supply. First, results indicate that adding an groundwater system whose supply productivity is just 8.67 % of the entire water requirement with a surface water supply first (SWSF) policy can significantly decrease the shortage index (SI) from 2.93 to 1.54. Second, the proposed model provides a more efficient conjunctive use policy than the SWSF policy, achieving further decrease from 1.54 to 1.13 or 0.79, depending on the groundwater rule curves. Finally, because of the usage of the hybrid framework, GA, LP, and ANN, the computational efficiency of proposed model is higher than other models with a purebred architecture or traditional groundwater numerical simulations. Therefore, the proposed model can be used to solve complicated large field problems. The proposed model is a valuable tool for conjunctive use operation planning.

Keywords Conjunctive use management· Genetic algorithm · Artificial neural network· Linear program and hybrid architecture

1 Introduction

Due to increasing water demands and global warming, water shortages are occurring more frequently. Expanding human populations and economic development have increased the water demand (Chang et al.2009; Braga et al.1985; Rosegrant and Cai 2002; Jenkins et al. 2004; Pulido-Velazquez et al. 2006). Global warming also causes extreme climate conditions, such as droughts and floods, which are occurring more often than in previous decades (Tung et al. 2006). These extreme climate conditions greatly influence water resource stability. The two conditions mentioned above increase the shortage risk of water supply in general. Because of the high interaction between climate conditions and surface water resources, groundwater is more stable than surface water. Thus, the highly efficient conjunctive use of surface and subsurface water is a good way to reduce the water shortage risk with minimal environmental impact for an area with plentiful groundwater aquifer (Mishra and Desai2006). Therefore, developing a conjunctive use management model (CUMM) to manage the surface and subsurface water efficiently and simultaneously has become an important issue (Vedula et al.2005; Pulido-Velazquez et al.2006).

Traditionally, the simulation-optimization method is one of the most widely used approaches for water resource management planning (Yeh 1985; Dhar and Datta 2008). The simulation-optimization method simulates the water resource using a numerical simulation model. This simulation model is then embedded in an optimization algorithm to determine the optimal operation policy of water supply. However, this simulation-optimization method for developing a CUMM requires the construction of numerical simulation models for both a multi-reservoir system and a groundwater system. Because the behaviors of reservoirs are linear, the implementation of the simulation model for the multi-reservoir system can be formed as a set of linear equations.

Although the model implementation for the multi-reservoir system is simple, the implementation for the groundwater system especially for unconfined aquifer is more complicated. Gorelick (1983) determined that the groundwater simulation approaches used in groundwater management models can be classified as embedded

or response matrix methods, depending on how the management model performs groundwater modeling. The embedded method uses a full function numerical model for groundwater simulation to represent the relationship between the pumping rates and groundwater head. Wang and Zheng (1998) embedded MODFLOW, one of the most commonly used models for groundwater flow simulation, in an optimal ground-water management model. Chang et al. (1992,2009) and Hsiao and Chang (2002) also embedded ISOQUAD, a finite element model for groundwater simulation proposed by Pinder and Frind (1972), in their optimal management model. This embedded method describes the groundwater system with less simplification than the response matrix method, and can be directly applied to confined or unconfined aquifers. The embedded method computes the state variables, hydraulic heads, for all numerical grids and provides detail description of the groundwater system. However, it also requires a large amount of computational resources. A management model generally requires only the information at the monitoring locations, and not all the numerical grids. From the management model perspective, the embedded method uses a great amount of computational resources to compute information that are not required. Therefore, the embedded method is not a computationally efficient approach for conjunctive use management problems.

The response matrix method simplifies the groundwater simulation as a linear matrix equation and computes only the information required for a management model. Psilovikos (2006) used the response matrix method to simulate the ground-water system for a confined aquifer. The advantage of the response matrix method is its high computational efficiency compared to the embedded method. However, due to its linear assumption, the response matrix method is not directly applicable to nonlinear systems, such as unconfined aquifers. However, unconfined aquifers are normally the top layers of groundwater system and provide much more water resources than confined aquifers.

Based on the discussion above, the embedded method and response matrix approach both have their limitations for conjunctive use management problems with unconfined aquifers. An Artificial Neural Network (ANN) is a computationally efficient alternative for addressing the problem. Rogers and Dowla (1994) integrated ANN and the Nonlinear Programming (NLP) into a groundwater remediation model. They used an ANN model to represent the groundwater system behaviors, including groundwater flow and pollutant transport. Based on previous research, Rogers et al. (1995) applied a combination of an ANN with the Genetic Algorithm (GA). The GA avoids the convergence problem of NLP induced by initial solutions. Coppola et al. (2003) proposed an ANN model for field groundwater simulation. Chang et al. (2005) applied an ANN model in a real-time reservoir operation. Chu and Chang (2009) combined the recursive ANN with Constrained Differential Dynamic Programming (CDDP) for a groundwater resource planning problem. They used the ANN model in the planning model instead of using ISOQUAD, effectively reducing the computational time by up to 94.5 %.

For the optimization algorithm of solving surface water allocation problems, LP is the most popular choice because of the linearity of reservoir systems (Labadie2004; Wei and Hsu2008). In the optimization of groundwater management or conjunctive use problems, because of the nonlinearity of groundwater systems, algorithms which have ability to deal with the nonlinearity such as NLP are required (Rogers and Dowla1994). Without the derivation of optimization formulation, GA also becomes an attractive alternative (Rogers et al.1995; Wang et al.2011).

The behaviors of multi-reservoir system and groundwater system in the CUMM are linear and nonlinear, respectively. Therefore, the CUMM optimization algorithm used must deal with mixed linear and nonlinear problems. These problems are complicated and difficult to solve using single algorithms. Thus, this study proposes a decomposition method to decompose the problem into several sub-problems based on problem characteristics such as integer, linear, nonlinear or time-variant charac-teristics. A hybrid architecture of optimal algorithms can then solve the set of sub-problems and each algorithm for solving different sub-sub-problems can be chosen based on the characteristics of sub-problem itself. For example, LP is usually the best choice for linear sub-problems because of its high searching efficiency and low computa-tional burden. Watkins and McKinney (1998) applied two decomposition methods, generalized Benders decomposition (GBD) and outer approximation (OA), to a mixed-integer nonlinear programming (MINLP) problem. They decomposed the MINLP problem into a mixed-integer linear programming (MILP) sub-problem and a nonlinear programming (NLP) sub-problem. Hsiao and Chang (2002) and Chang et al. (2009) also used decomposition methods in their groundwater management models. They proposed optimal groundwater management models to minimize the total costs, including fixed costs for well network design and operational costs for groundwater pumping. The models in both studies simultaneously determined the optimal design of well networks and the associated time-variant pumping quantities. Because the determination of a network design contains integer characteristics and the determination of pumping quantities contains time-variant characteristics, the entire problem becomes a mixed integer and time-variant problem. They decom-posed the entire problem into an integer sub-problem, which is a network design problem, and another time-variant sub-problem, which is a pumping quantity deter-mination problem. Based on the sub-problem characteristics, the algorithms for both sub-problems respectively are GA, a good choice for integer programming problems, and CDDP (Chang et al.1992), which a dynamic programming-based algorithm and also the most powerful algorithm for time-variant problems. These studies show that the decomposition method can successfully decompose the complicated water resource management problem into several sub-problems. The proposed hybrid architecture, unlike a purebred architecture, can combine the advantages of multiple algorithms for solving these sub-problems.

This study proposes a stepwise optimization CUMM for unconfined aquifer. The proposed model uses the decomposition method to decompose the entire opti-mization problem into a linear sub-problem (a surface water system), and another nonlinear sub-problem (an unconfined aquifer system). The proposed model uses LP and GA to solve these two sub-problems. This study uses a recursive ANN model to simulate the behavior of unconfined aquifer.

2 Methodology

This section describes the methodology of developing the proposed CUMM. Section 2.1. demonstrate the model architecture of CUMM. Section2.2. uses the recursive ANN model to simulate the behavior of an unconfined aquifer. Section2.3. further describes some basic operating concepts, including the rule curve operating method, the index value theory, the index balance method (IBM), and the

mathe-matical formulation for CUMM. Section2.4. presents a hybrid architecture including ANN, LP, and GA, which is then used to determine the supply strategies for CUMM. 2.1 Model Architecture

Figure1shows the model architecture of the proposed CUMM. In Fig.1, the model architecture contains three elements, include GA, LP and ANN. The ANN instead of numerical modeling for unconfined aquifers is trained to simulate an unconfined groundwater system (Eq. 12) and the construction of ANN is further described in Section2.2. The combination of GA and LP is the optimization solver of CUMM which determine the optimal pumping rates for well networks and release quan-tities for reservoirs. Because the behavior of reservoir system and the behavior of unconfined aquifers respectively are linear and nonlinear, a decomposition method is applied to decompose the CUMM into two sub-problems, linear and nonlinear parts. GA and LP are respectively used to solve the different parts. The combination of two algorithms can remain the high computational efficiency of LP for the linear part and the robust, flexibility and nonlinear ability of GA for the nonlinear part. 2.2 Groundwater ANN Model Development

The governing equation for unconfined aquifer is a nonlinear equation which describes the temporal and spatial variations of groundwater levels. In traditional groundwater numerical simulation, only the groundwater levels of the first time step are directly assigned as the initial condition. The groundwater levels of later time steps are recursively calculated based on the previous time step. This subsection pro-poses a recursive ANN model to replace the usage of traditional numerical models. Due to the recursive procedure of the time variation problem mentioned above, the recursive ANN model also simulates as the recursive procedure mentioned above.

Figure2shows a structural diagram of the recursive ANN model. The structure can be divided into two stages: a training stage (Fig. 2a), and a prediction stage (Fig.2b). The training stage has a structure similar to other ANN models. The input nodes of the structure contain state variables and decision variables of the t-th time step. The output nodes only contain state variables of the next time step. In CUMM,

Genetic Algorithm

(Decision variable of primary sub-problem: pumping rates)

Artificial Neural Network (groundwater transfer function)

Linear Programming (Decision variable of secondary

sub-problem: reservoir supply amounts)

Output layer (5 nodes) Hidden layer (20 nodes) Input layer (10 nodes) ... t=t+1 ... (a) (b) ( )t h ( )t h P( )t P( )t ( )t+1 h ( )t+1 h

Fig. 2 Recursive ANN structure diagram a Training stage; b Prediction stage

the state variables mean the groundwater levels for each observation well, h(t) and

h (t + 1), and the decision variables mean the pumping rates for each pumping well, P (t). The training stage normally collects data pairs between input nodes and output nodes from field data observations or numerical model generation. The prediction stage modifies the structure of the recursive ANN model (Fig.2b). Similar to the recursive procedure mentioned above, the results calculated from output nodes are recursively assigned to the input nodes for the next step simulation, except for the first time step.

Figure 3 illustrates the creation procedure of the recursive ANN model. This procedure can be divided into two parts: data creation for ANN model training and verification and ANN model training and verification. The first part contains two steps: data creation by groundwater numerical model and data preprocessing for ANN. The first step uses the numerical groundwater model, MODFLOW 96, to generate data for ANN training and verification. Figure4depicts the groundwater system structure. Five pumping wells, well #A to #E, and five observation wells, wells #1 to #5, and five observation wells, wells #1 to #5, were installed on a homogeneous and isotropic unconfined aquifer whose area is 17× 17 (km2_{) (shown as Table 1).} The aquifer thickness was 110 (m). The aquifer was divided into 170× 170 cells whose area was 100× 100 m2 _{each. The left side and right side of the aquifer were} assigned as the Dirichlet boundary conditions, which were respectively assigned as 100 and 80 (m). The other two sides were assigned as no flow boundary condition (Fig.4). The length of simulation time step was 10 (days) and the entire simulation period was 365 time steps (10 years total). Table1shows the hydrogeological and other simulation parameters.

The training data set was generated by repeating MODFLOW simulations with different pumping rates. A random number generator was used to create ten series

ANN training data and verification data creation Data creation by groundwater numerical model

Groundwater numerical model construction

Pumping rates randomly generation for each wells

Groundwater numerical model simulation

Data preprocessing for ANN

Data collection

Groundwater levesl for each oberservation wells Pumping rates for each

pumping wells

Data normalization

ANN groundwater simulation model training and verification

ANN model training

ANN verification fail

ANN groundwater simulation model construction completing

success ANN structure or parameter

modification

Fig. 3 Flowchart of training and verification of groundwater ANN simulation model

of time varying pumping rates with 365 time steps for the five pumping wells, and ten matrixes labeled as P365×5were created to store the pumping rates. After repeating MODFLOW simulations, the groundwater levels, including the initial time step and other 365 time steps observed from the five observation wells, were stored in another ten matrixes labeled as h366×5.

The data preprocessing step can be divided into two sub-steps, including data collection and data normalization (Fig. 3). The data collection sub-step collects training data and verification data from the matrixes, P365×5and h366×5. In the data normalization sub-step, a linear transformation transfers the variable values varied between 1 to -1 is applied. Both the aquifer thickness (110 m) and maximum pumping rate (0.5 cms) are used to normalize h366×5and P365×5, as Eq.1shows.

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ PiN,365×5= Pi 365×5× 2 0.5 (cms) − 1 hi N,366×5= hi 366×5× 2 110(m) − 1 , for i = 1, ..., 10 (1)

Fig. 4 Map of groundwater system

The ANN groundwater simulation model training and verification stage uses the MATLAB Neural Network toolbox to construct the proposed ANN model. Table2

lists the ANN parameters used in this study. The proposed model uses a three-layer backward propagation network (BPN), including an input three-layer, a hidden layer, and an output layer. The number of nodes in the hidden layer is 20 and the transfer function of hidden layer and output layer are respectively tan-sigmoid transfer function (Eq. 2), and linear transfer function (Eq.3). The assignment for

Table 1 Table of aquifer

properties and simulation parameters Parameter Value Simulation parameter Simulation area 17× 17 (km2) Cell area 100× 100 (m2₎ Number of cells 170× 170 Length of simulation time step 10(days) Entire simulation period 3, 650 (days)

= 10 years = 365 time steps Maximum pumping rate 0.26 × 106(m3/ten days)

= 0.5 (cms) Aquifer property

Aquifer thickness 110(m) Specific yield, (sy) 0.2

Porosity, (n) 0.2

Horizontal hydraulic conductivity, 0.0001 (m/s) (Kh)

Vertical hydraulic conductivity, 0.0001 (m/s) (K_v)

Table 2 Table of ANN

parameters Parameter Value

Node number for input layer 10 Node number for hidden layer 20 Node number for output layer 5

Hidden layer transfer function tan-sigmoid transfer function Output layer transfer function Linear transfer function Training algorithm BFGS

RMSE criteria for training 10−7

the number of output layer is based on the size of the observation wells. To avoid the affect of over fitting, a try-and-error approach is used to determine the proper number of hidden nodes. The number of hidden nodes increases gradually when the ability of the ANN model with fewer hidden nodes can not deal with the complexity and nonlinearity of data. Until the model errors for both the training set and the validation set can satisfy the model criteria, the current number of hidden nodes is the proper number. The BFGS method, a qusai-Newton method, is the training algorithm. The ANN convergence criterion is that the root mean-square error (RMSE) must be smaller than 10−7.

f(x) = e

x_{− e}−x

ex_{+ e}−x (2)

f(x) = x (3)

After finishing the training stage, model verification involving another set of data pairs is necessary to verify model reliability. The ANN model structure after the training stage becomes the prediction stage (Fig. 2b). Therefore, by applying the initial groundwater levels and other sets of pumping rates, the RMSE between the simulated levels and the observation levels is an index of model reliability.

Two types of time series for pumping created by an uniform random number, type A, and a step function, type B, are used to train and verify the robustness and ability. The type A was generated by the uniform random number generator which varied from 0 to 0.5 (cms). The type B was generated by the step function generator which pumps 0.3 (cms) for 18 continuous time steps and stops for the next 18 time steps. Each verification cases simulate 3,650 time steps (100 years). The verification results in Table3indicate that the relative root-mean-square error (relative RMSE) of a 100-year simulation is less than 3.0 %. These result confirm the reliability of the

Table 3 ANN model verification error

Verification of case #A Verification of case #B

RMSE (m) Related RMSE RMSE (m) Related RMSE

(%) (%) Well #1 0.44 0.56 1.47 1.84 Well #2 0.47 0.67 1.98 2.75 Well #3 1.31 1.93 1.24 1.79 Well #4 0.44 0.56 1.53 1.92 Well #5 0.46 0.65 1.93 2.68

proposed ANN model and indicate that the ANN groundwater simulation model is well trained.

On a computer with Intel(R) Xeon(R) CPU E5620 @2.40GHz, the computational burden of the groundwater simulation with 25, 000 time steps for the MODFLOW model with 170× 170 grids consumes about 131 seconds. However, the computa-tional burden of the simulation with the same time steps for the groundwater ANN model just requires two seconds. The ANN model is about 65 times faster than the MODFLOW model. The high computational efficiency of the ANN model make the proposed conjunctive-use model attractive in a operation with large number of time steps.

2.3 Formulation of Conjunctive-Use Management Model (CUMM)

This section describes three operational concepts: Rule Curve Operating Method (Section 2.3.1), Index Value Theory (Section 2.3.2) and Index Balancing Method (Section2.3.3). The Rule Curve Operating Method is traditionally used to determine the total supply quantity based on current reservoir storage. The Index Value Theory linearly demonstrates the status of reservoir storage. The Index Balancing Method is used to allocate the supply quantities between different reservoirs in a multi-reservoir system. Because this study develops a CUMM, it applies the Rule Curve Operating Method and Index Value Theory to groundwater systems and extends the Index Balancing Method to allocate the supply quantities between surface water system and groundwater system. Section 2.3.4 describes the formulation of the proposed CUMM.

2.3.1 Rule Curve Operating Method

The rule curve operating method is one of the most widely used approaches in practice because engineers can easily operate the reservoirs basing on the predefined rule curves. These rule curves vertically divide the entire reservoir volume into several layers. Each layer associates a discounting ratio, which varies from 0 to 1, and the ratio for upper layer is normally larger than the ratio for lower layer. The ratio is calculated based on the sum of reservoir storage and inflow quantity (Eq.4).

Dt₌ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ γ1_Dt _{i f ν}2_{> V}t_{+ IF}t_{≥ ν}1 γ2_Dt _{i f ν}3_{> V}t_{+ IF}t_{≥ ν}2 ... γNr_Dt _{i f} νNr+1> Vt+ IFt≥ νNr (4)

In Eq.4, Vt_{is the reservoir storage at the t-th time step, I F}t_{is the reservoir inflow} quantity at the same time step,νi _{is the volume of the reservoir storage below the} bottom of the i-th layer,γi_{is the discounting ratio associated with the i-th layer, N}

r is the number of the rule curves, Dt is the original water demand, and Dt _{is the} discounted water demand and also is the supply quantity for the current time step. The discounting ratios for upper layers are normally equal to 1.0 and those for lower layers are less than 1.0. The reason for the discounting ratio assignment is that the reservoir can fully satisfy the water demand when the reservoir storage is plentiful. Otherwise, the reservoir should save water to go through the entire draught period.

The upper layers and the lower layers can respectively be called the full supply layer (FSL) and partial supply layer (PSL).

The rule curve operating method also can be applied to a multi-reservoir system. The equation of discounting ratio determining should be modified as Eq. 5. The storage quantities, inflow quantities, and reservoir volumes of each reservoir are summed to form the equivalent reservoir values for the entire system.

Dt= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ γ1_Dt _{i f} Ns j=1ν2j > Ns j=1  Vt j+ IFtj  ≥Ns j=1ν1j γ2_Dt _{i f} Ns j=1ν3j > Ns j=1  Vtj+ IF t j  ≥Ns j=1ν2j ... γNr_Dt _{i f} Ns j=1νNjr+1> Ns j=1  Vt j+ IFtj  ≥Ns j=1νNjr (5)

2.3.2 Index Value Theory

The index value for a reservoir is linearly normalized by the volume of the reservoir storage and its associating rule curves (Eq.6). In Eq.6, It

SW, jmeans the index value at the t-th time step for the j-th reservoir, Vt

jmeans the volume of the water storage at the same time step,νi

jmeans the reservoir storage at the bottom of the i-th layer andνij+1means the reservoir storage at the bottom of the i+ 1-th layer (i.e., the top of the i-th layer). Therefore, an index value between i and i+ 1 means that the water level occupies the i-th layer.

It_{SW, j}= i +Vtj− ν i j  /νi+1 j − ν i j  , i = 1, ..., Nr (6) ItGW= i + h t_{− μ}i_{/ μ}i+1_{− μ}i_{, i = 1, ..., N} r (7)

The concept of the Index Value Theory can be further applied in a conjunctive use system. The definition of the index value of a groundwater system can be also written as Eq.7. In Eq.7, It

GWmeans the index value of groundwater system at the t-th time step, ht_{means the average groundwater level at the same time step,μ}i_means the average groundwater level at the bottom of the i-th layer andμi+1_{means the} average level at the top of the same layer. Therefore, an index value between i and

i+ 1 means that the average level occupies the i-th layer. 2.3.3 Index Balancing Method

Equations 4 and 5 calculate the discounted water demand. In a multi-reservoir system, the supply quantities for each reservoir and the pumping quantities for groundwater aquifers must be further allocated. This study uses the index balance method (IBM) to allocate the supply quantities in a multi-reservoir system or a conjunctive use system.

The IBM is an operating policy in which the index values of each reservoir should equal to the values of other reservoirs during operation. The goal of the IBM is to avoid overuse of single reservoir in contrast to the entire water resource system, and to guarantee that the storage of each reservoir has similar status. The IBM of

a conjunctive use system is similar to that of a multi-reservoir system. The index value of a surface water system and the value of a groundwater system should be the same or nearly equal.

2.3.4 CUMM Formulation

The operation of the proposed CUMM is based on the predefined operation rule curve. The discounted ratio and discounted water demand for each time step is calculated directly using the rule curve. Although the discounted demand at each time step can be calculated, the supply quantities for each reservoir and the pumping quantities for groundwater system must obey the IBM and the pre-calculated dis-counted demand (Eqs.4and5). The determination of supply quantities and pumping quantities can be formed as a stepwise optimization problem. Equation8through20

show this formulation. 

Problem #1 (Original problem for CUMM) Objective function: ZCU MMt = min  Xt_{, P}t Nd  i=1 ωSH,iSH Rti+ ωG AP Ns  m=1 Ns  n=1 G APtSW,mn+ (8) ωG APG APtSY S+ ωSL ⎛ ⎝Ns j=1 SLtSW, j+ SLtGW ⎞ ⎠ , ωSH,i> ωG AP> ωSL s.t. SH Rt i= SHti/ Nd k=1 Dt k  , ∀i ∈ Nd (9) Dti= ∀ j∈Ns j=1 λijXtj+ ∀k∈Np k=1 λikPtk+ SH t i, ∀i ∈ Nd, Xtj∈ Xt, Ptk∈ P t (10) Vtj+1= V t j+ ∀m∈NI F m=1 λjmI Fmt − ∀i∈Nd i=1 λijXtj− OF t j, ∀ j ∈ Ns (11) ht+1_{= f} AN N  ht_{, P}t (12) G APtSW,mn= |Itsw,m− Istw,n|, ∀m, n ∈ Ns (13) G APt SY S= |ItSW− IGWt | (14)

It SW, j= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1+V t j−ν1j ν2 j−ν1j i fν 2 j> V t j≥ ν1j 2+V t j−ν2j ν3 j−ν2j i fν 3 j> Vtj≥ ν2j ... Nr− 1 + Vt j−νNrj νNr +1 j −νNrj i fνNr+1 j > Vtj≥ ν Nr j , ∀ j ∈ Ns (15) ItGW= ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1+ ht_{− μ}1_{/ μ}2_{− μ}1 _{i fμ}2_{> h}t_{≥ μ}1 2+ ht− μ2_{/ μ}3_{− μ}2 _{i fμ}3_{> h}t_{≥ μ}2 ... Nr+ ht− μNr  / μNr+1− μNr _{i f}μNr+1> _ht≥ μNr (16) SLtSW, j= Vj,max− Vtj Vj,max− Vj,min, ∀ j ∈ Ns (17) SLtGW= hmax− ht hmax− hmin (18) Vj,min≤ Vtj≤ Vj,max, ∀ j ∈ Ns (19) hmin≤ ht≤ hmax (20)

where Xt_{and P}t_{respectively represent the water supply quantities and} ground-water pumping quantities at the t-th time step and also are the decision variables for [Problem #1]. The terms Nd, Ns, Np, NI F, and Nr are the number of water demands, reservoirs, pumping wells, inflow nodes and rule curve levels, respectively. The lengths of Xt _{and P}t _{are N}

s and Np, respectively. SH Rti is the ratio of the water shortage to water demand for the i-th demand node (Eq.9). G APt

SW,mnmeans the index value discrepancy between the m-th and the n-th reservoirs (Eq. 13).

G APtSY Smeans index value discrepancy between the surface water system and the groundwater system (Eq. 14). SLtSW, j means the ratio of spare volume to system capacity for the j-th reservoir (Eq.17). SLt

GWmeans ratio of spare volume to system capacity for groundwater system (Eq.18). Dt

imeans the discounted water demand for the i-th demand node and is already pre-determined based on Eq.4or5. SHt i means water shortage for the i-th demand, while Vtj means water storage of the

j-th reservoir. The terms ht and ht respectively means the groundwater levels of

different observation points and the average groundwater level, while I Ft

mand OFtj respectively mean the inflow quantity for the m-th inflow node and the overflow quantity for the j-th reservoir. It

SW, j and ItSW respectively mean the index value

of the j-th reservoir and the average index value for entire surface water system.

It

GW means the index value for groundwater water system.νmj andμmrespectively mean the volume of the m-th level of rule curve for the j-th reservoir and that for groundwater system. Vj,min and Vj,maxare the minimum and maximum volumes of the j-th reservoir, respectively, while hminand hmaxare the minimum and maximum average groundwater levels. The termλij means the active index between the i-th

node and the j-th node. If the i-th node and the j-th node are connected, the value of

λijis 1; otherwise, the value is 0.

Equation 8 is an objective function that sums shortage rates for each demand (SH Rtj), the summation of discrepancy of index value (G AP

t

SW,mn and G APtSY S) and the summation of the ratio of spare volume to system capacity (SLtSW, j and

SLt

GW). Three weight values,ωSH, j,ωG APandωSL, represent the priority between different sub-objectives. The shortage rate of each demand is further defined by Eq.9

in which the water shortage, SHt

j, is defined through Eq.10, a mass balance equation for each demand. In this study, the priority of the summation of shortage rates is higher than the other two sub-objectives. The second sub-objective, which minimizes the summation of discrepancy of index value, tries to balance the index value among reservoirs or different systems. G APt

SW,mnand G APtSY S, respectively are defined in

Eqs.13and14, are the discrepancy of index value between different reservoirs or different systems. The third term of objective function, minimizing the summation of the ratio of spare volume to system capacity, tries to keep more water resource in the supply system and also be further formed as Eqs.17and18.

Equations 11 and 12, which are the continuity equations of reservoirs and groundwater system, present the behavior of reservoirs and groundwater system. The implementation of Eq.12repuires a recursive ANN model instead of a traditional numerical model or a response matrix, as introduced in Section 2.2. Equations 19

and 20 respectively are constraint functions for reservoirs and the groundwater system. These constraint functions require that the reservoir levels or groundwater levels must vary in predefined ranges.

Since the unconfined groundwater transfer function, Eq. 12, is a nonlinear constraint, the optimal problem, [Problem #1], becomes a complicated nonlinear optimization problem. However, this problem can be reformulated as a linear sub-problem and another nonlinear sub-sub-problem by exploring the sub-problem structure. The original problem is denoted as [Problem #1], and the decomposed problem is denoted as [Problem #2].

Problem #2 (Decomposed problem for CUMM) Objective function of primary sub-problem

ZCU MMt = min Pt  JPt+ ωG APG APtSY S+ ωSLSLtGW  (21) s.t. Equations12,14,16,18and20

Objective function of secondary sub-problem

JPt= min  Xt Nd  i=1 ωSH,iSH Rti+ ωG AP Ns  m=1 Ns  n=1 G APtSW,mn+ ωSL Ns  j=1 SLtSW, j (22) s.t. Equations9,10,11,13,15,17and19

[Problem #2] contains a primary sub-problem and a secondary sub-problem. The primary sub-problem, which determines the pumping quantities ( Pt_{) on an} unconfined aquifer system, is still a nonlinear problem. However, the secondary sub-problem, which just deals with the surface water system under the pre-determined pumping quantities ( Pt_{), is a linear problem.}

2.4 Hybrid Architecture for Solving CUMM

As described previously, this study uses a hybrid architecture to solve [Problem #2]. Figure1shows the model architecture of the proposed CUMM. In Fig.1, GA and LP are respectively used to solve the primary problem and the secondary sub-problem. The ANN described in Section1is also adapted to simulate the unconfined

Copy, Crossover and Mutation Operator Penalty Function Yes Yes

Randomly create an initial population (GA for primary sub-problem)

Decode chromosomes (pumping rates, P)

Calculate groundwater level by using ANN

Is the groundwater level satisfied water level constrain?

Calculate J*(P) by using LP (Secondary sub-problem)

Calculate objective function of Primary sub-problem

Is the GA convergent?

Hybrid algorithm stop for t-th time step

Does the simulation reach the last time step?

No No No Stop Yes t=t+1

groundwater system (Eq. 12). Figure5 shows the computation steps of the hybrid architecture, as illustrated in detail below.

1. The initial population in GA is generated for the solution of the primary sub-problem. Each chromosome is encoded as a binary string.

2. Decoding each binary string, chromosome, into a set of decimal value, a set of pumping rates ( Pt_).

3. The groundwater ANN model (Eq.12) computes the average groundwater level (ht+1_{) for the next time step reacted by assigning the set of pumping rates} decoded from each chromosome. If the average level is not satisfied the ground-water level constrain (Eq.20), a penalty function will be added; otherwise, LP further solves the secondary sub-problem to obtain the objective value, J

 Pt

, which is also the first term of Eq.21. Because the second and the third term of Eq.21are functions of the average groundwater level, the objective function value of secondary sub-problem for each chromosome can be determined based on the Eqs.12,14,16and18.

4. After calculating the fitness (the value of objective function in primary sub-problem), examining the convergence criteria. If it is not convergent, apply selection, crossover, and mutation operators to generate a new population of chromosomes and return to Step 2; otherwise, stop entire process.

Using the proposed hybrid architecture has several benefits. First, solving surface water allocation problem (the secondary sub-problem) is computationally more efficient using LP than other algorithms, such as GA. Second, because the value

of J( Pt_{) in the objective function of the primary sub-problem is the optimal value}

of the secondary sub-problem, the gradient of J( Pt) is difficult to be derived. Therefore, GA is a better choice than traditional gradient-based optimal algorithms such as NLP. Therefore, the GA-LP hybrid architecture for the entire problem, the primary and the secondary sub-problems, not only achieves higher computa-tional performance for the linear part, but also has ability to deal with the non-linear part.

3 Results and Discussion

This section introduces the study cases and their associated setting in Section 3.1. Section3.3then describes the operation results of the cases and compares different cases to show the advantages of the proposed CUMM.

3.1 Cases Introduction

The length of operation time step is ten days, same as the length of simulation time step used in the recursive ANN model. Table4compares four numerical cases with the same water requirement, 15× 106_(m3_{/ten days). Case #1 only uses surface} water resources to satisfy the water requirement, and thus a pure surface water operation case (refers to SW case). The formulation of supply strategy determination for the SW case is determined using only the secondary sub-problem in Problem #2. Because the objective function of the secondary sub-problem is conditioned with the

Table 4 Numerical case list table

Operation form Groundwater rule curve Model usage Case # 1 Surface water independent operation No rule curves SWAM

(SW)

Case # 2 Decoupled conjunctive use operation No rule curves CUMM (DCU)

Case # 3 Coupled conjunctive use operation Groundwater rule curve #1 CUMM (CCU)

Case # 4 Coupled conjunctive use operation Groundwater rule curve #2 CUMM (CCU)

predefined pumping quantities, Pt_{, the value of P}t_{is directly assigned as zero in Case} #1, a SW case.

Cases #2, #3 and #4 are conjunctive use cases. Case #2 uses both surface water and groundwater resources in a prearranged supply priority called surface water supply first (SWSF) policy. This case is also called a decoupling conjunctive use (refers to DCU) case. The groundwater system supplies only during the water demand cannot be satisfied by only using the surface water system. The determination of the supply strategy for the surface water system in this case is equal to that in Case #1. After determining the supply quantities of the surface water system, the pumping quantities can be assigned as either the gap between the surface water supply quantities and the water demand or the maximum pumping productivity constraint. The DCU model has been widely applied in Taiwan because of its ease of implementation.

Cases #3 and #4 use both surface water and groundwater resources with the pro-posed CUMM. Contrast to the Case #2, the supply priority which is not predefined before model operating is determined by the CUMM. Different operation rule curves of the groundwater resource which define the priority of resource usage are applied in Cases #3 and #4. Both of them are coupled conjunctive use (refers to CCU) cases.

Fig. 6 System diagram of

surface water system (SW case) Required demand 1.5 X 107_(m3_/ten-day) 1 SEA max storage 7.0X107(m3) 1.0 2.0 3.0 4.0 5.0 1 1.0 2.0 3.0 4.0 5.0 B r i o v r e s e R A r i o v r e s e R max storage 5.0X107(m3) min storage 0(m3) min storage 0(m3)

• Case #1:

The surface water system contains two reservoirs, reservoir A and reservoir B, with maximum storages of 70× 106 _{and 50× 10}6 _(m3_{), respectively. These two} reservoirs supply water a water demand node with a parallel connection. Figure6

is a system diagram. The Fig. 7 shows the variations of inflow quantities for two reservoirs. The average inflow quantities for the two reservoirs are 8.4 × 106 _and

4.95 × 106_(m3_{/ten days), respectively. The average inflow quantity of entire surface}

system, 13.35 × 106 _(m3_{/ten days) is a little smaller than the water requirement} quantity, 15× 106 _(m3_{/ten days). Figures 8 and 9 show the operation rule curve} diagrams of reservoirs A and B, respectively. In this case, the discounting ratios for each layers are assigned as 1.0 and each layers are FSLs.

• Case #2, #3, #4:

The other three cases are all conjunctive use operation cases. The two-reservoir surface water system mentioned above and an additional five-well groundwater system construct the entire supply system (Fig.10). This study uses a recursive ANN model as the simulation model for the five-well groundwater system, as described in Section 2.2. The reservoir assignment of the surface water system, including maximum capacities, inflow quantities, and operation rule curves, is the same in all four cases, from Cases #1 to #4. In the groundwater system, the maximum pumping rate of each well is 0.26 × 106 _(m3_{/ten days). Therefore, the maximum} supply productivity of the entire groundwater system is 1.30 × 106_(m3_{/ten days).}

Because Cases #3 and #4 are both CCU cases, their supply priorities are de-termined by comparing the gap of SI between surface and groundwater system, which is related to the assigned operation rule curves for both systems. Groundwater

0 5 10 15 20 0 50 100 150 200 250 300 350 400 Inflow Quantity (X10 6 m 3/ten-days) Reservoir A 0 5 10 15 20 0 50 100 150 200 250 300 350 400 Inflow Quantity (X10 6 m 3/ten-days)

Time Step No.

Reservoir B

0 10,000,000 20,000,000 30,000,000 40,000,000 50,000,000 60,000,000 70,000,000 80,000,000 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 Stora g e (m 3)

time step (10 days) Reservoir A rule curves

Index=4.0

Index=3.0

Index=2.0

Index=1.0

Fig. 8 The diagram of operation rule curve for reservoir A

0 10,000,000 20,000,000 30,000,000 40,000,000 50,000,000 60,000,000 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 Sto rag e (m 3)

time step (10 days) Reservoir B rule curves

Index=4.0

Index=3.0

Index=2.0

Index=1.0

Fig. 9 The diagram of operation rule curve for reservoir B

Fig. 10 System diagram of

conjunctive use system(CU system) Required demand 1.5 X 107_(m3_/ten-day) 1 SEA max storage 7.0X107_(m3₎ 1.0 2.0 3.0 4.0 5.0 1 1.0 2.0 3.0 4.0 5.0 B r i o v r e s e R A r i o v r e s e R max storage 5.0X107_(m3₎ min storage 0(m3₎ min storage 0(m3₎ Groundwater system max supply 1.3X106_(m3_/ten-day)

70.00 75.00 80.00 85.00 90.00 95.00 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 A v erag e g roundw ater lev el (m )

time step (10 days) Operation rule curve I

Index=4.0 Index=3.0 Index=2.0

Index=1.0 Index=5.0

Fig. 11 Operation rule cuve #1 of groundwater system

operation rule curves #1 and #2 (Figs.11and12, respectively) are applied to Cases #3 and #4. The third layer, whose index value varies from 3.0 to 4.0, in Fig. 11is thinner than the layer in Fig.12. On the contrary, the first layer, whose index value varies from 1.0 to 2.0, in Fig.11is thicker than the layer in Fig.12. Figure13plots the relationship between the average groundwater level and the groundwater index value for both rule curves #1 and #2. For example, if the groundwater level is 85 (m), the index values based on different rule curves are 2.0 and 3.67, and if the groundwater level is 75 (m), the index values are 1.33 and 3. Therefore, the index value based on the rule curve #1 is always smaller than the value based on the rule curve #2 for the same groundwater level. The lower index values mean the lower priority of groundwater usage in rule curve #1.

70.00 75.00 80.00 85.00 90.00 95.00 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 A v erag e g roundw ater lev el (m )

time step (10 days) Operation rule curve II

Index=4.0

Index=3.0 Index=2.0 Index=1.0 Index=5.0

0 1 2 3 4 5 70 75 80 85 90 Groundwater level (m) Inde x va lu e Rule curve #1 Rule curve #2

Fig. 13 The relationship diagram between groundwater level and index value for rule curve #1 and #2

3.2 Model Parameters for Genetic Algorithm

Table5lists the model parameters for GA. The GA population contains 50 chro-mosomes. The mutation rate and crossover rate respectively are 0.01 and 0.8. The proposed model based on two stopping criteria. The first one is the model will be stopped while the model generations is exceeded the maximum generation. The second one is the model will be stopped while the optimal fitness of each generation is continuously fixed without improvement in a specific length of generations. The maximum generation of the proposed model is 200 and the number of non-improved generation is 20.

3.3 Case Results

Table6lists the simulation results of four cases. • Case #1:

In Case #1 (SW), because the average supply amount of surface water system,

13.55 × 106 _(m3_{/ten days), is 90.33 % of the water requirement quantity, 15 × 10}6

(m3_{/ten days), the average supply amount is almost equal to the average inflow}

quantity mentioned above. The shortage index (SI) of Case #1 is 2.93 and the number of supply shortage time steps (labels as SSTSN) is 135.

Table 5 Table of GA

parameters Parameter Value

Size of GA population 50

Mutation rate 0.01

Crossover rate 0.8

Maximum generation 200

Number of non-improved generations 20 for stopping criteria

Table 6 Table of simulation results

Case type Case #1 Case #2

SW DCU

Shortage index (SI) 2.93 1.54

Number of supply shortage time step (SSTSN) 135 125 Average requirement quantity (m3_{/ten days) 15.00 15.00}

Average supply amount of surface water system 13.55 (90.33 %) 13.55 (90.33 %) (m3_{/ten days)}

Maximum supply productivity of groundwater system 0 (0 %) 1.30 (8.67 %) (m3/ten days)

Average supply amount of groundwater system 0 (0 %) 0.48 (3.20 %) (m3_{/ten days)}

Average production efficiency of groundwater system 0 % 36.87 % Average supply amount of the conjunctive use system 13.55 (90.33 %) 14.03 (93.53 %)

(m3/ten days)

Case type Case #3 Case #4

CCU CCU

Shortage index (SI) 1.13 0.79

Number of supply shortage time step (SSTSN) 111 89 Average requirement quantity (m3_{/ten days) 15.00 15.00}

Average supply amount of surface water system 13.48 (89.87 %) 13.53 (90.20 %) (m3_{/ten days)}

Maximum supply productivity of groundwater system 1.30 (8.67 %) 1.30 (8.67 %) (m3/ten days)

Average supply amount of groundwater system 0.68 (4.53 %) 0.89 (5.93 %) (m3_{/ten days)}

Average production efficiency of groundwater system 52.31 % 68.46 % Average supply amount of the conjunctive use system 14.16 (94.40 %) 14.420 (96.13 %)

(m3/ten days)

• Case #2:

In Case #2 (DCU), the average supply amount of surface water system is also

13.55 × 106 _(m3_{/ten days) and the operation behavior is the same as that of Case}

#1. The groundwater system is worked within the water demand cannot be satisfied by using only the surface water system. The average supply amount of groundwater system is 0.48 × 106 _(m3_{/ten days), which is just 3.20 % of the entire requirement.} Because of the additional groundwater resource, the SI of Case #2 decreases from 2.93 to 1.54 and the SSTSN decreases from 135 to 125. Although the supply ratio of groundwater system is just 3.20 %, the additional groundwater resource still can significantly decrease the SI in a DCU case.

A comparison between Cases #1 and #2 indicates that a small but stable water resource also can significantly relax the deficit status of water supply shortage. Com-paring the average pumping quantity, 0.48 × 106 _(m3_{/ten days), and the maximum} productivity constrain, 1.30 × 106 _(m3_{/ten days), the production efficiency of the} groundwater system in Case #2 is just 36.87 %. Because the production efficiency caused by the SWSF policy is still low, the production efficiency of this groundwater system remains great potential to improve with various conjunctive use policies.

• Case #3 and #4:

In Cases #3 and #4, both CCU cases, the operation priority between two systems is determined by the assignment of rule curves for surface water system and groundwa-ter system. In both cases, the average supply amounts of surface wagroundwa-ter respectively are 13.48 × 106 _(m3_{/ten days), which is 89.87 % of the entire water requirement,} and 13.53 × 106 _(m3_{/ten days), 90.20 %. The amounts of groundwater respectively} are 0.68 × 106 _(m3_{/ten days), 4.53 %, and 0.89 × 10}6_(m3_{/ten days), 5.93 %. The SI} values respectively are 1.13 and 0.79, while the SSTSNs are 111 and 89.

A comparison of the DCU case, Case #2A comparison of the DCU case, Cases #3 and #4, shows that the average supply ratios of groundwater system increase from 3.20 % to 4.53 % and to 5.93 %. The production efficiencies also increase from 36.87 % to 52.31 % and to 68.46 %. The SIs further decrease from 1.54 to 1.13 and to 0.79. These results indicate that the proposed conjunctive use policy can increase the production efficiency of the groundwater system and also decrease the system SI to relax the shortage status of water supply.

The improvement of productivity efficiency for the groundwater systems between Cases #2 and Case #4 is more significant than that between Case #2 and Case #3. The reason of the difference on improvement is that the different height of operation rule curves of groundwater system (Figs.11 and12). Figure14compares SW and GW index values for the operation results of Case #3 and #4, while Fig.15 shows the associated groundwater level. Because of the effect of the operation rule curves (Fig.13), the slope of GW level and index value for rule curve #1 is steeper than the slope for rule curve #2 when the index values vary from 3.0 to 4.0 (the third layer of the rule curve). Besides, because the IBM is the supply principle to allocate quantities between groundwater system and surface water system, the groundwater usage in Case #3 has lower priority and shorter pumping period than in Case #4. Counting the number of time steps whose GW index value is higher than the associated SW index value from Fig.14, the number of Case #3, which is 187 within 365 time steps, is less than that of Case #4, which is 247.

3 4 5

Index value

SW Index value(Case #3) GW Index value(Case #3)

SW Index value(Case #4) GW Index value(Case #4)

0 1 2

1 51 101 151 201 251 301 351

Index value

Time step No.

80 85 90

Groundwater Level (m)

Average water level of Case #3 Average water level of Case #4

70 75 80

1 51 101 151 201 251 301 351

Groundwater Level (m)

Time step No.

Fig. 15 The variation diagram of groundwater level for Case #3 and #4

Figure 16 demonstrates the convergency of GA of Case #3 for the 15-th time step. In the figure, the red line demonstrates the variation of the optimal fitnesses in each generations and the blue points demonstrates the fitness values for each chromosomes. The model converges at the 59-th generation while the optimal fitnesses remain 1.1228 from the 40-th to the 59-th generation. For the convergency of GA, the optimal fitness is stable indicate that the proposed model has had high possibility to find the optimal solution.

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0 10 20 30 40 50 60 Fitness Generation Time Step = 15 Optimal Fitness Chromosome Fitness

4 Conclusion

The CUMM, the proposed stepwise optimal model for conjunctive use management, has been developed in this study and a hybrid architecture, the combination of LP, GA and a recursive ANN, is also proposed in CUMM. In the hybrid architecture, LP is applied to determine the optimal amounts of supply quantities for each reservoir, GA is applied to connect a surface water system and a groundwater system and is also used to determine the optimal amounts of pumping quantities for each wells, and the recursive ANN is trained to represent the behaviors of an unconfined aquifer. The hybrid architecture used in the CUMM has several advantages. First, for groundwater simulation, the usage of the recursive ANN can remain both the high computational efficiency of the response matrix method and the nonlinear simulation ability of the embedded method for unconfined aquifers. In simulations of unconfined aquifers with large number of time steps, the recurse ANN is 65 times faster than the MODFLOW model and, besides, the error of the ANN related to the MODFLOW model remains small. Therefore, the high computational efficiency and the nonlinear simulation ability make the proposed model have the potential to apply on a real case planning for unconfined aquifer. Second, for supply policy determination, the optimal algorithm combination in the hybrid architecture, LP and GA, can remain the high computational efficiency of LP for solving the linear sub-problem and the robust, flexibility and nonlinear ability of GA for solving the nonlinear sub-problem. Based on not only the high computational efficiency of the recursive ANN groundwater simulation but also the advantages of the hybrid optimal algorithms, the proposed hybrid architecture can be further extending to a field case. The usage of the operation rule curves for multi-reservoir systems is the most widely-used approach to determine the priority between different reservoirs and allocate reservoir supply quantities. In this study, the operation rule curves were further applied to determine the priority between multi-reservoir and groundwater systems. Compared with the DCU case, the proposed CCU model, CUMM, could improve the SI from 1.54 to 1.13 or 0.79. The different improvement efficiencies depend on the definition of the groundwater operation rule curves. The thicker the upper layer of groundwater system, the lower the SI for the CCU model. The comparison of different rule curves shows that the different thickness of the upper layer defined by rule curves significantly affected water usage priority between reservoir system and groundwater system. Comparing Case #3 and #4 shows that the average production efficiency for groundwater system significantly increased from 52.31 % to 68.46 %.

The operation rule curves for traditional reservoir operation typically vary with seasons. Because water resources are more plentiful during flood seasons, the heights of the operation rule curves are lower, increasing the usage priority of reservoir storage. Because the heights of the operation rule curves for groundwater system remain static over seasons, the variation of the groundwater rule curves allows the CUMM to further improve the system SI in future study.

Finally, this study uses the operation rule curve theory as a new type of framework for conjunctive use management. Because the operation rule curve theory is widely used for practical reservoir management, the framework is easy to understand for not only academic researchers, but also engineers, who are responsible for practical water management.

Acknowledgements The authors would like to thank the National Science Council of Republic of China (Taiwan), which supported this research under contract NSC 99-2918-I-009-010-.

References

Braga BPF, Conejo JGL, Becker L, Yeh WWG (1985) Capacity expansion of Sao-Paulo water-supply. J Water Resour Plan Manag-ASCE 111:238–252

Chang LC, Chen YW, Yeh MS (2009) Optimizing system capacity expansion schedules for ground-water supply. Water Resour Res 45:W07407

Chang LC, Shoemaker CA, Liu PLF (1992) Optimal time-varying pumping rates for groundwa-ter remediation—application of a constrained optimal-control algorithm. Wagroundwa-ter Resour Res 28:3157–3173

Chang Y-T, Chang L-C, Chang F-J (2005) Intelligent control for modeling of real-time reservoir operation, part II: artificial neural network with operating rule curves. Hydrol Process 19:1431– 1444

Chu HJ, Chang LC (2009) Optimal control algorithm and neural network for dynamic groundwater management. Hydrol Process 23:2765–2773

Coppola E, Szidarovszky F, Poulton M, Charles E (2003) Artificial neural network approach for predicting transient water levels in a multilayered groundwater system under variable state, pumping, and climate conditions. J Hydrol Eng 8:348–360

Dhar A, Datta B (2008) Optimal operation of reservoirs for downstream water quality control using linked simulation optimization. Hydrol Process 22:842–853. doi:10.1002/hyp.6651

Gorelick SM (1983) A review of distributed parameter groundwater-management modeling meth-ods. Water Resour Res 19:305–319

Hsiao CT, Chang LC (2002) Dynamic optimal groundwater management with inclusion of fixed costs. J Water Resour Plan Manag 128:57

Jenkins MW, Lund JR, Howitt RE, Draper AJ, Msangi SM, Tanaka SK, Ritzema RS, Marques GF (2004) Optimization of California’s water supply system: Results and insights. J Water Resour Plan Manag-ASCE 130:271–280

Labadie JW (2004) Optimal operation of multireservoir systems: state-of-the-art review. J Water Resour Plan Manag 130(2):93–111

Mishra AK, Desai VR (2006) Drought forecasting using feed-forward recursive neural network. Ecol Model 198:127–138

Pinder GF, Frind EO (1972) Application of galerkin’s procedure to aquifer analysis. Water Resour Res 8:108–120

Psilovikos A (2006) Response matrix minimization used in groundwater management with mathe-matical programming: a case study in a transboundary aquifer in northern greece. Water Resour Manag 20:277–290

Pulido-Velazquez M, Andreu J, Sahuquillo A (2006) Economic optimization of conjunctive use of surface water and groundwater at the basin scale. J Water Resour Plan Manag-ASCE 132: 454–467

Rogers LL, Dowla FU (1994) Optimization of groundwater remediation using artificial neural networks with parallel solute transport modeling. Water Resour Res 30:457–481

Rogers LL, Dowla FU, Johnson VM (1995) Optimal field-scale groundwater remediation using neural networks and the genetic algorithm. Environ Sci Technol 29:1145–1155

Rosegrant MW, Cai XM (2002) Global water demand and supply projections part—2. Results and prospects to 2025. Water Int 27:170–182

Tung CP, Lee TY, Yang YC (2006) Modelling climate-change impacts on stream temperature of formosan landlocked salmon habitat. Hydrol Process 20:1629–1649

Vedula S, Mujumdar PP, Sekhar GC (2005) Conjunctive use modeling for multicrop irrigation. Agric Water Manag 73:193–221

Wang K-W, Chang L-C, Chang F-J (2011) Multi-tier interactive genetic algorithms for the optimiza-tion of long-term reservoir operaoptimiza-tion. Adv Water Resour 34:1343–1351

Wang M, Zheng C (1998) Ground water management optimization using genetic algorithms and simulated annealing: formulation and comparison. J Am Water Resour Assoc 34:519–530

Watkins DW, McKinney DC (1998) Decomposition methods for water resources optimization mod-els with fixed costs. Adv Water Resour 21:283–295

Wei C-C, Hsu N-S (2008) Multireservoir real-time operations for flood control using balanced water level index method. J Environ Manag 88:1624–1639

Yeh WWG (1985) Reservoir management and operations models—a state-of-the-art review. Water Resour Res 21:1797–1818