# Optimizing Reservoir-Stream-Aquifer Interactions for Conjunctive Use and Hydropower Production

Authors: Hala Fayad; Richard C. Peralta; Ali Forghani
Journal: Advances in Civil Engineering
 (2012)

---

## Abstract

Conjunctive management of water resources involves coordinating use of surface water and groundwater resources. Very few simulation/optimization (S-O) models for stream-aquifer system management have included detailed interactions between groundwater, streams, and reservoir storage. This paper presents an S-O model doing that via artificial neural network simulators and genetic algorithm optimizer for multiobjective conjunctive water use problems. The model simultaneously addresses all significant flows including reservoir-stream-diversion-aquifer interactions in a more detailed manner than previous models. The model simultaneously maximizes total water provided and hydropower production. A penalty function implicitly poses constraints on state variables. The model effectively finds feasible optimal solutions and the Pareto optimum. Illustrated is application for planning water resource and minihydropower system development.


---

## Body

## 1. Introduction

Conjunctive management of water resources involves coordinating use of surface water and groundwater resources, often to address demands of competing water users. Conjunctive management can reduce temporal water deficits by using groundwater when surface water supplies are scarce. Where surface water and groundwater resources are hydraulically connected, understanding the interaction between the two resources is important.Water resources management cannot be optimized without the ability to simulate the system’s response to management and hydrologic stimuli. Simulation models predict physical system response to natural and managed stimuli. On the other hand, a simulation-optimization (S-O) model directly computes the water management strategy that best satisfies desired goals without causing unacceptable system responses. It employs user-specified objectives to drive the search for an optimal strategy.The embedding technique and response matrix approach are the two methods generally used to incorporate groundwater simulation abilities within an S-O model [1–6]. For transient problems in which relatively few cells include decision or state variables, the response matrix method is more practical and has been generally used.The response matrix method relies upon linear systems theory and superposition [7]. Confined aquifers are linear systems. Nonlinear systems, such as unconfined aquifers, can also be treated as linear if changes in saturated thickness and transmissivity in time are proportionally small with respect to initial values [8]. Many researchers adopted linear treatment to satisfactorily address nonlinear systems [5, 6, 9–11]. Others used cycling or sequential linear programming (SLP) to adapt linear systems theory to address nonlinear groundwater problems [2–4, 6, 12–14].Most reported conjunctive use of S-O models uses classical optimizers relying upon simplex and gradient search approaches. O’Mara and Duloy [15] used an agricultural production economic model to plan management of a river-aquifer system. Walker [16] coordinated surface water and groundwater systems serving independent irrigation areas. Peralta et al. [17, 18] and Da Conceicao Morais et al.[19]  used crop production functions to allocate surface water and groundwater. Ejaz and Peralta [20] optimized conjunctive use with surface water quality constraints. Belaineh et al. [11] included reservoir management rules and detailed interactions between reservoir storage and releases, surface water diversions, and groundwater pumping. Barlow et al. [5] evaluated tradeoffs between groundwater withdrawal and streamflow depletion.Because most conjunctive use problems are nonlinear, classical nonlinear programming optimizers have been used often [18, 21–25]. However, because the convexity of the objective function and feasible region cannot always be ensured, global optimality of computed solutions is not always guaranteed. Furthermore, the process of computing derivatives can consume much CPU time.Consequently, researchers employed heuristic optimizers that do not require computing derivatives. Such methods include simulated annealing [23, 26], and Genetic Algorithms (GAs) [27–38].Heuristic optimizers are coupled with groundwater and other simulation models. Both normal numerical finite simulators and surrogate simulators are used. Surrogates such as artificial neural networks (ANNs) are used to reduce computational burden during optimization. ANNs are trained using field data or a calibrated numerical model. A trained ANN simulates much more quickly than the original numerical model.Rogers et al. [30] used GA optimization with substitute simulators, artificial neural networks (ANNs). The ANNs suitably represented state variable response to management. Rao et al. [23] used SA with ANNs to perform combinatorial optimization for simple nonlinear, nonconvex, conjunctive use allocation. Karamouz et al. [33] addressed conjunctive use and water quality using GAs and ANNs. Safavi et al. [37] used an ANN to represent surface water-groundwater interaction and a GA to minimize unsatisfied demand in three irrigation systems. However, they did not consider reservoir storage.Some of the above efforts involved multiple measures of performance or objectives that could be optimized simultaneously. Multiple objective (MO) optimization permits evaluating tradeoffs between conflicting goals and aids identifying acceptable compromise solutions. When considering all objectives simultaneously in MO problems, a set of solutions exists that are superior to all other solutions in the search space. The superior solutions are known as “Pareto Optimal” solutions or “Non Dominated” solutions [39, 40].This paper presents a multiobjective simulation/optimization model that uses artificial neural networks simulation, GA optimization, and constraints for optimizing conjunctive use of a nonlinear hydraulically linked reservoir-stream-aquifer system. Trained using results of simulations from a spatially distributed finite difference simulation model, the ANNs represent groundwater head, reservoir level and stream stage, response to boundary conditions, and the water management strategy being optimized. State variables are discouraged from having undesirable values by penalty functions. Penalty functions are usually derived from experience and are subject to imprecision and uncertainty. Decision variables include spatially distributed groundwater pumping, reservoir diversion, stream diversion, and reservoir release rates. The linked ANN-GA model robustly addresses multiple objectives for a complicated nonlinear reservoir-stream-aquifer system and includes more significant flows than previously reported ANN-based models (Table1 and Figure 1).Table 1
Parameters associated with the streams, lake, and aquifer.
ParametersValuesReservoir/lakeArea2.3 × 106 m2Capacity42 × 106 m3Cutoff elevation29.5 mTurbine elevation20.0 mInitial stage29.0 mDiversion channel capacity14115 m3/dayPenstock Capacity (release to hydropower)200,000 m3/dayHydraulic Conductivity of lakebed0.04 m/dayThickness of lakebed sediment0.9 mStreamsWidth9.1 mThickness of streambed sediment0.3 mSlope of streambed segments0.0002Hydraulic conductivity of streambed sediment0.04 m/dayManning’s roughness coefficient0.02Streambed conductance180 m2/dayStream inflow to lakeTotal length6090 mInflows in stress period 1225840 m3/dayInflows in stress period 256460 m3/dayInflows in stress period 3112920 m3/dayStream outflow from lakeTotal length2435 mFirst diversion channel capacity8750 m3/daySecond diversion channel capacity14115 m3/dayAquiferHydraulic conductivity of unconfined layer61 m/DaySpecific yield of the unconfined layer0.15Transmissivity of the confined layer371 m2/dayStorage coefficient of the confined layer0.0001Ratio of horizontal to vertical hydraulic conductivity100Groundwater recharge0.8 × 10−3m/dayPrecipitation rate2.2 × 10−3m/dayEvaporation rate1.74 × 10−3m/dayIrrigated field650×103 m2Figure 1
Conceptual view of the hypothetical study area.


## 2. Study Area

The proposed model is applied to a hypothetical study area based upon that of Cheng and Anderson [41], modified to more comprehensively illustrate conjunctive use issues. The area includes an influent stream, multipurpose reservoir and discharges, effluent stream and diversions, groundwater extraction wells, and municipal, industrial, and irrigation users (Figure 1). Wells pump from a homogeneous anisotropic two-layer aquifer system represented by a finite difference model (Figure 2). The total study area size is 11.15 × 106 m2.Figure 2
Discretization and cell types of the hypothetical study area.The Figure3 conceptual cross-section illustrates the saturated hydraulic connection between the 18.2 m deep reservoir and the aquifer. The 9.1 m wide stream is in excellent hydraulic connection with the aquifer and penetrates only the unconfined layer. Wells extract only from the confined aquifer layer. From the reservoir, flows one diversion and a release that passes through turbines, producing hydropower before re-entering the stream. Table 1 data show that the reservoir can provide a maximum head of 9.5 m (29.5 m cutoff elevation minus 20 m). Table 1 also contains other distinguishing characteristics and parameters of the study area.Figure 3
Distorted cross section of study area at reservoir location.Table2 shows different industrial, municipal, and irrigation water demands in the region and identifies sources that can satisfy them. The irrigated area receives water from the first stream diversion. The reservoir diversion, the second stream diversion, and the pumping wells deliver water to four industrial areas and one small rural town. Conveyance losses are insignificant because water flows through new pipes and lined concrete canals. The four industrial and municipal users either consume water completely or discharge unconsumed water outside the study area after treatment. Assumed losses for the irrigated area total 30 percent of applied water (18 percent percolates as recharge to the unconfined aquifer and 12 percent returns to the stream as return flow).Table 2
Water demands in the study area.
Water userCharacteristicSourceDemand (m3/day)Industry 1Pumping well-11660Industry 2Reservoir diversion1383Industry 3Pumping well-24150Industry 4Stream diversion-21245Irrigated area650,000 m2Stream diversion-13873Town600 personsPumping well-2 and reservoir diversion180


## 3. Methodology

### 3.1. Simulation of Reservoir-Stream-Aquifer Interactions

#### 3.1.1. Introduction

The artificial neural networks used here are trained to produce some of the same outputs as a finite difference numerical simulation model. After training, the ANNs are embedded within the presented S-O model. The GA heuristic optimizer uses the ANNs during optimization, rather than using the finite difference model. The following explanation of the process begins by discussing the original simulator.


#### 3.1.2. Direct Simulation Using Finite Difference Model

Initially employed to simulate reservoir-stream-aquifer flows is the USGS MODFLOW finite difference groundwater flow model [42], with the LAK2 lake package [43], and STR stream routing package [44]. LAK2 simulates reservoir-groundwater interaction and reservoir water stage. STR computes flow in the stream, stream stage, and stream-aquifer seepage. Seepage is a function of stream stage and aquifer head. Together, these are termed the MODLAKE model. MODLAKE simulated the system under transient conditions using three sixty-day stress periods and Table 1 parameter values.Input values for MODLAKE simulations include pumping rates at wells 1 and 2, stream diversion rates 1 and 2, reservoir diversion, and reservoir release. For 3 stress periods, these total 18 values. For each period, seven important computed output system states descriptors are heads at both wells in both confined and unconfined strata, stream stage at two control locations, and reservoir stage.The discharge from the reservoir that can occur via the emergency spillway is not a model input value. Occurring only when reservoir stage is above the Table1 cut-off value, it flows to the stream (Figure 1). That discharge is computed via rating equation:(1)Qemergency,k=Const×[hlake,k-Cut-offelevation]Exponent,
where Const and Exponent parameter values are 1385649 and 1.66667, respectively. These values are obtained following Cheng and Anderson [41] where stream outflow is governed by Manning’s equation.


#### 3.1.3. Indirect Simulation Using Artificial Neural Networks

To develop the data base for ANN training, many MODLAKE simulations were performed for a range of input flow values. For each simulation, the above-mentioned 18 flow values were inputs, and 21 state variables were recorded outputs. Note that, during subsequent optimization, the S-O model will determine optimal values for these 18 decision variables, and will ensure that these 21 state variables have satisfactory values.One ANN was trained to represent each of the 21 state variables. Table3 shows the number of input, hidden, and output nodes for ANNs describing period 1 state variables. Because there are six input flows in period one, Table 3 shows that the input layer has six nodes. The number of hidden nodes within the ANN hidden layer was determined through trial and error. Because a separate ANN was developed for each state variable, there is only one output node per ANN.Table 3
Evaluation results for 7 ANNs in stress period 1.
State variableNumber of nodes in ANNCorrelationRRatio of test to training RMSLake stage6-6-10.9991.049Head 1 layer 16-4-10.9980.804Head 2 layer 16-4-10.9920.416Head 1 layer 26-6-10.9990.488Head 2 layer 26-8-10.9950.366Stream stage 16-7-10.9011.183Stream stage 26-3-10.9680.812Networks for the first stress period were developed using input data from the first stress period only (six input nodes corresponding to six decision variables). Networks for the second stress period were developed using input data from the first and second stress periods. Therefore, all networks in this period had 12 input nodes, corresponding to 12 decision variables. Networks for the third period required 18 input nodes.Here, the NeuralSIM commercial package was used to train and validate the 21 ANNs for the 21 state variables. NeuralSIM uses the cascade method of network construction with an adaptive gradient learning rule. The back-propagation learning method was used.To adequately train ANNs for the highly nonlinear reservoir-stream-aquifer system, many data were used. About 400 sets of decision variables were created randomly within a broad range, and the same number of MODLAKE simulations was performed. Then, additional 300 sets of decision variables were developed for a narrow range, and were simulated using MODLAKE. The combined sets of decision and state variables constituted the ANN training data.Training involves adjusting ANN parameters until ANN output is close to desired (target) values. The above MODLAKE output provided the target values. According to NeuralSIM, an indicator of an adequate ANN is that training set performance and test set performance are fairly similar. Figure4 displays ANN-predicted versus MODLAKE-simulated (actual) lake stage in stress period 1.Figure 4
ANN-predicted versus actual lake stage.Two means of comparing ANN and MODLAKE outputs were used linear correlation(R) and root mean squared error (RMS).Linear correlation values were 0.875 to 0.999 for all ANNs. This high correlation indicates that the real-world target is highly explained by the ANN-predicted values. Also, for almost all networks, the ratio of the RMS of the test data to the RMS of the training data was less than 1.45. Table3 provides the evaluation results for 7 ANNs in stress period 1. Hsu et al. [45] showed detail ANN training.


### 3.2. Optimization Problem

Water resources allocation problems can involve complicated hydrological, environmental, and economical constraints and conflicting management objectives. The posed study area has adequate water supply for existing water uses, but planners want to evaluate best to increase water use while producing hydropower. These two objectives conflict.


#### 3.2.1. Water Supply Objective

For an area relying on an aquifer and a stream for its water supply, the objective is to maximize the water provided from pumping wells and from surface diversions. Surface diversions include both diversions from streams and diversions from surface water reservoirs. Equation (2) is the objective function of maximizing water provided to users:(2)Z1=Max⁡(1T)∑k=1T[∑â=1NPCâ,kPpâ,k+∑ê=1NSCê,kSdê,kstr+∑î=1NRCî,kRdî,kres],
where Z1is the average flow rate delivered during all periods, pâ,k is the rate of water pumped from cell â during period k, dê,kstr is the rate of stream water diverted at diversion point ê during period k, and dî,kres is the rate of reservoir water diverted at diversion point î during time period k. NP, NS, and NR are total numbers of groundwater pumping, stream diversions, and reservoir diversions, respectively. Câ,kP, Cê,kD, and Cî,kR are weighting coefficients assigned to groundwater pumping, stream diversion, and reservoir diversion, respectively.Optimal solutions are sensitive to the weighting coefficients. Careful selection of coefficient values permits their use for economic optimization and emphasizing or de-emphasizing specific decision variables. A variable can be made ineffective in the objective function by setting its respective weighting coefficient equal to zero.


#### 3.2.2. Hydropower Production Objective

Water released from the lake through a turbine produces hydropower as per the following:(3)Hydropower=epγQhp,k86,400hg,k,
where Hydropower is in kilowatts (kw), ep is the overall efficiency of the power plant (the product of hydraulic and turbine efficiencies), γ is the specific weight of water (9.81 KN/m3), Qhp,k is the reservoir release flow rate (m3/d), and hg is the gross head difference between the reservoir water surface and the point where water passes through the turbine (m).Rearranging yields (4), the objective function for maximizing the average hydropower produced during the management era:(4)Z2=Max⁡epγ×∑k=1T((Qhp,k86,400)(hlake,k-hturbine)(1T)),
where, for this site, ep is 0.65, hturbine is the elevation of the point where water passes through the turbine (20 m), and hlake,k is the lake stage at end of time period k.


#### 3.2.3. Bounds and Constraints

Constraints enforce physical laws, legal rights, and management goals. The most simple constraints are lower and/or upper bounds on decision variables (Table4). Upper limits on well pumping rates assure that pumpage from a given well does not exceed the capacity of the well, pump, and pipe. Upper limits on reservoir and stream diversions prevent exceeding the capacity of the diversion canals. Lower limits on stream diversions consider legal rights. The upper bound on reservoir release (Qhp,k) equals the smaller of the turbine capacity and the penstock capacity. The lower bound reflects legal or economic reasons.Table 4
Bounds on decision variables in stress period 1.
Decision variableLower bound (L) (m3/day)Upper bound (U) (m3/day)Pumping well 13102820Pumping well 2308460Stream diversion 1758460Stream diversion 23602820Reservoir diversion282820Reservoir release57680195370More complicated constraints involving decision variables are demand constraints and constraints limiting total stream diversion. Demand constraints ensure that total water delivered to each user at least satisfies existing demand:(5)∑k=1T(pâ,k+dê,kstr+dî,kres)≥∑k=1TDemandk,
where Demandk is the total current demand for water in time period k. The intent is to possibly provide more water for future development.Upper bounds on cumulative stream diversions downstream of the reservoir ensure that stream diversions do not exceed the available water in the stream. In other words, the sum of stream diversions downstream of the reservoir must be always less than or equal the reservoir release to the stream.Each state variable that must be bounded during optimization is represented within the S-O model as an ANN. The GA represents their bounds using the penalty method. Table5 shows the imposed bounds on state variables.Table 5
Bounds on state variables.
ConstraintEquationConstraintStress period 11 (Lake stage)825.1 ≤ lake stage ≤ 29.3 m2 (Head 1 layer 1)626.2 ≤ head 1 layer 1 ≤ 30.2 m3 (Head 2 layer 1)622.9 ≤ head 2 layer 1 ≤ 29 m4 (Head 1 layer 2)620.7 ≤ head 1 layer 2 ≤ 29 m5 (Head 2 layer 2)617.1 ≤ head 2 layer 2 ≤ 27.1 m6 (Stream stage 1)720.5 ≤ stream stage 1 ≤ 20.6 m7 (Stream stage 2)720.2 ≤ stream stage 2 ≤ 20.3 mStress period 28 (Lake stage)821.9 ≤ lake stage ≤ 29.6 m9 (Head 1 layer 1)621.9 ≤ head 1 layer 1 ≤ 28 m10 (Head 2 layer 1)620.6 ≤ head 2 layer 1 ≤ 27.1 m11 (Head 1 layer 2)616.5 ≤ head 1 layer 2 ≤ 26.5 m12 (Head 1 layer 2)615.2 ≤ head 2 layer 2 ≤ 26.2 m13 (Stream stage 1)720 ≤ stream stage 1 ≤ 20.1 m14 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mStress period 315 (Lake stage)822.6 ≤ lake stage ≤ 28 m16 (Head 1 layer 1)618.9 ≤ head 1 layer 1 ≤ 27.4 m17 (Head 2 layer 1)618.9 ≤ head 2 layer 1 ≤ 25.6 m18 (Head 1 layer 2)614.3 ≤ head 1 layer 2 ≤ 25.3 m19 (Head 2 layer 2)614.6 ≤ head 2 layer 2 ≤ 24.4 m20 (Stream stage 1)719.9 ≤ stream stage 1 ≤ 20 m21 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mLimits on reservoir stage ensure that the water does not bypass the turbine and that the reservoir does not face environmental dangers:(6)hlake,kL≤hlake,k≤hlake,kU.
To prevent unacceptable drawdown, the hydraulic head in the aquifer at each well â should be within lower (L) and upper (U) bounds for each time period:(7)hâ,kL≤hâ,k≤hâ,kU.
Lower limit on stream stage at each control location û ensures adequate flow for recreational, navigational, legal, and environmental reasons:(8)stageû,k≥stageû,kL.


### 3.3. Genetic Algorithm and Modeling Methodology

Genetic algorithms (GAs) can be powerful and robust optimizers for nonlinear problems, if using appropriate optimization parameters and enough simulations. A simple genetic algorithm (SGA) is a search procedure based on the mechanics of natural selection and genetics. It uses selection and recombination operations while searching a solution space to identify the best solution. SGAs have been generally applied to single-objective optimization problems.However, many optimization problems have multiple objectives. When optimizing with GAs, multiple objectives have usually been addressed by either combining all objectives into a scalar objective function (consisting of a sum of weighted multiple objectives) or by converting some objectives into constraints having thresholds and penalty functions. Other methods for handling multiple objectives includeε-constraint method, goal programming, and Multi objective genetic algorithms (MOGAs).The complexity of conjunctive use problems has led some researchers to use MOGAs. MOGAs also facilitate linking nonlinear optimizers and groundwater quantity and quality simulators. An efficient MOGAs is the nondominated sorting genetic algorithm (NSGA) [46]. NSGA uses a ranking selection method to emphasize current nondominated (Pareto optimal) points and a sharing function method to maintain diversity in the population.For simulation, the presented model uses coupled ANN and finite difference simulation to represent more detailed flow interactions, than other reported models for multiobjective reservoir-stream-aquifer system management. The S-O model uses NSGA for optimization. It includes objective or fitness functions, constraints, 21 linked ANNs, and the optimization algorithm. Using the optimizer requires setting five parameters values properly. These are mutation probability (Pmutate), cross-over probability (Pcross), niche size (σshare), population size, and maximum number of generations. σshare is dependent on the number of peaks in the solution space. We used a mutation probability of 0.01 and a cross-over rate of 0.6. The stopping criterion was 500 generations. The initial population size was 2758 individuals. The number of considered peaks was 100. These values were derived via sizing criteria suggested by Mahfoud [47] and sensitivity analysis. Continuous decision variables were represented by binary coding. NSGA selection of strategies employed the stochastic remainder proportionate method [46].Figure5 shows the simulation/optimization modeling procedure.(1)
Generate numerous water management strategies. Each strategy is a set of pumping, diversions, and reservoir release rates. Run MODLAKE for each strategy. Store variable system responses for all stress periods.(2)
Using strategies and results of step 1, train an ANN to enable it to predict state variable response to other sets of decision variables.(3)
Apply the S-O model to compute optimal strategies. The model calls the ANN trained in step 2.(4)
Confirm that output strategies indeed satisfy all constraints. This is necessary because the ANNs are not perfectly accurate. To do this, run MODLAKE for the optimal water management strategies and check the state variable results.Figure 5
S-O model processing.


## 3.1. Simulation of Reservoir-Stream-Aquifer Interactions

### 3.1.1. Introduction

The artificial neural networks used here are trained to produce some of the same outputs as a finite difference numerical simulation model. After training, the ANNs are embedded within the presented S-O model. The GA heuristic optimizer uses the ANNs during optimization, rather than using the finite difference model. The following explanation of the process begins by discussing the original simulator.


### 3.1.2. Direct Simulation Using Finite Difference Model

Initially employed to simulate reservoir-stream-aquifer flows is the USGS MODFLOW finite difference groundwater flow model [42], with the LAK2 lake package [43], and STR stream routing package [44]. LAK2 simulates reservoir-groundwater interaction and reservoir water stage. STR computes flow in the stream, stream stage, and stream-aquifer seepage. Seepage is a function of stream stage and aquifer head. Together, these are termed the MODLAKE model. MODLAKE simulated the system under transient conditions using three sixty-day stress periods and Table 1 parameter values.Input values for MODLAKE simulations include pumping rates at wells 1 and 2, stream diversion rates 1 and 2, reservoir diversion, and reservoir release. For 3 stress periods, these total 18 values. For each period, seven important computed output system states descriptors are heads at both wells in both confined and unconfined strata, stream stage at two control locations, and reservoir stage.The discharge from the reservoir that can occur via the emergency spillway is not a model input value. Occurring only when reservoir stage is above the Table1 cut-off value, it flows to the stream (Figure 1). That discharge is computed via rating equation:(1)Qemergency,k=Const×[hlake,k-Cut-offelevation]Exponent,
where Const and Exponent parameter values are 1385649 and 1.66667, respectively. These values are obtained following Cheng and Anderson [41] where stream outflow is governed by Manning’s equation.


### 3.1.3. Indirect Simulation Using Artificial Neural Networks

To develop the data base for ANN training, many MODLAKE simulations were performed for a range of input flow values. For each simulation, the above-mentioned 18 flow values were inputs, and 21 state variables were recorded outputs. Note that, during subsequent optimization, the S-O model will determine optimal values for these 18 decision variables, and will ensure that these 21 state variables have satisfactory values.One ANN was trained to represent each of the 21 state variables. Table3 shows the number of input, hidden, and output nodes for ANNs describing period 1 state variables. Because there are six input flows in period one, Table 3 shows that the input layer has six nodes. The number of hidden nodes within the ANN hidden layer was determined through trial and error. Because a separate ANN was developed for each state variable, there is only one output node per ANN.Table 3
Evaluation results for 7 ANNs in stress period 1.
State variableNumber of nodes in ANNCorrelationRRatio of test to training RMSLake stage6-6-10.9991.049Head 1 layer 16-4-10.9980.804Head 2 layer 16-4-10.9920.416Head 1 layer 26-6-10.9990.488Head 2 layer 26-8-10.9950.366Stream stage 16-7-10.9011.183Stream stage 26-3-10.9680.812Networks for the first stress period were developed using input data from the first stress period only (six input nodes corresponding to six decision variables). Networks for the second stress period were developed using input data from the first and second stress periods. Therefore, all networks in this period had 12 input nodes, corresponding to 12 decision variables. Networks for the third period required 18 input nodes.Here, the NeuralSIM commercial package was used to train and validate the 21 ANNs for the 21 state variables. NeuralSIM uses the cascade method of network construction with an adaptive gradient learning rule. The back-propagation learning method was used.To adequately train ANNs for the highly nonlinear reservoir-stream-aquifer system, many data were used. About 400 sets of decision variables were created randomly within a broad range, and the same number of MODLAKE simulations was performed. Then, additional 300 sets of decision variables were developed for a narrow range, and were simulated using MODLAKE. The combined sets of decision and state variables constituted the ANN training data.Training involves adjusting ANN parameters until ANN output is close to desired (target) values. The above MODLAKE output provided the target values. According to NeuralSIM, an indicator of an adequate ANN is that training set performance and test set performance are fairly similar. Figure4 displays ANN-predicted versus MODLAKE-simulated (actual) lake stage in stress period 1.Figure 4
ANN-predicted versus actual lake stage.Two means of comparing ANN and MODLAKE outputs were used linear correlation(R) and root mean squared error (RMS).Linear correlation values were 0.875 to 0.999 for all ANNs. This high correlation indicates that the real-world target is highly explained by the ANN-predicted values. Also, for almost all networks, the ratio of the RMS of the test data to the RMS of the training data was less than 1.45. Table3 provides the evaluation results for 7 ANNs in stress period 1. Hsu et al. [45] showed detail ANN training.


## 3.1.1. Introduction

The artificial neural networks used here are trained to produce some of the same outputs as a finite difference numerical simulation model. After training, the ANNs are embedded within the presented S-O model. The GA heuristic optimizer uses the ANNs during optimization, rather than using the finite difference model. The following explanation of the process begins by discussing the original simulator.


## 3.1.2. Direct Simulation Using Finite Difference Model

Initially employed to simulate reservoir-stream-aquifer flows is the USGS MODFLOW finite difference groundwater flow model [42], with the LAK2 lake package [43], and STR stream routing package [44]. LAK2 simulates reservoir-groundwater interaction and reservoir water stage. STR computes flow in the stream, stream stage, and stream-aquifer seepage. Seepage is a function of stream stage and aquifer head. Together, these are termed the MODLAKE model. MODLAKE simulated the system under transient conditions using three sixty-day stress periods and Table 1 parameter values.Input values for MODLAKE simulations include pumping rates at wells 1 and 2, stream diversion rates 1 and 2, reservoir diversion, and reservoir release. For 3 stress periods, these total 18 values. For each period, seven important computed output system states descriptors are heads at both wells in both confined and unconfined strata, stream stage at two control locations, and reservoir stage.The discharge from the reservoir that can occur via the emergency spillway is not a model input value. Occurring only when reservoir stage is above the Table1 cut-off value, it flows to the stream (Figure 1). That discharge is computed via rating equation:(1)Qemergency,k=Const×[hlake,k-Cut-offelevation]Exponent,
where Const and Exponent parameter values are 1385649 and 1.66667, respectively. These values are obtained following Cheng and Anderson [41] where stream outflow is governed by Manning’s equation.


## 3.1.3. Indirect Simulation Using Artificial Neural Networks

To develop the data base for ANN training, many MODLAKE simulations were performed for a range of input flow values. For each simulation, the above-mentioned 18 flow values were inputs, and 21 state variables were recorded outputs. Note that, during subsequent optimization, the S-O model will determine optimal values for these 18 decision variables, and will ensure that these 21 state variables have satisfactory values.One ANN was trained to represent each of the 21 state variables. Table3 shows the number of input, hidden, and output nodes for ANNs describing period 1 state variables. Because there are six input flows in period one, Table 3 shows that the input layer has six nodes. The number of hidden nodes within the ANN hidden layer was determined through trial and error. Because a separate ANN was developed for each state variable, there is only one output node per ANN.Table 3
Evaluation results for 7 ANNs in stress period 1.
State variableNumber of nodes in ANNCorrelationRRatio of test to training RMSLake stage6-6-10.9991.049Head 1 layer 16-4-10.9980.804Head 2 layer 16-4-10.9920.416Head 1 layer 26-6-10.9990.488Head 2 layer 26-8-10.9950.366Stream stage 16-7-10.9011.183Stream stage 26-3-10.9680.812Networks for the first stress period were developed using input data from the first stress period only (six input nodes corresponding to six decision variables). Networks for the second stress period were developed using input data from the first and second stress periods. Therefore, all networks in this period had 12 input nodes, corresponding to 12 decision variables. Networks for the third period required 18 input nodes.Here, the NeuralSIM commercial package was used to train and validate the 21 ANNs for the 21 state variables. NeuralSIM uses the cascade method of network construction with an adaptive gradient learning rule. The back-propagation learning method was used.To adequately train ANNs for the highly nonlinear reservoir-stream-aquifer system, many data were used. About 400 sets of decision variables were created randomly within a broad range, and the same number of MODLAKE simulations was performed. Then, additional 300 sets of decision variables were developed for a narrow range, and were simulated using MODLAKE. The combined sets of decision and state variables constituted the ANN training data.Training involves adjusting ANN parameters until ANN output is close to desired (target) values. The above MODLAKE output provided the target values. According to NeuralSIM, an indicator of an adequate ANN is that training set performance and test set performance are fairly similar. Figure4 displays ANN-predicted versus MODLAKE-simulated (actual) lake stage in stress period 1.Figure 4
ANN-predicted versus actual lake stage.Two means of comparing ANN and MODLAKE outputs were used linear correlation(R) and root mean squared error (RMS).Linear correlation values were 0.875 to 0.999 for all ANNs. This high correlation indicates that the real-world target is highly explained by the ANN-predicted values. Also, for almost all networks, the ratio of the RMS of the test data to the RMS of the training data was less than 1.45. Table3 provides the evaluation results for 7 ANNs in stress period 1. Hsu et al. [45] showed detail ANN training.


## 3.2. Optimization Problem

Water resources allocation problems can involve complicated hydrological, environmental, and economical constraints and conflicting management objectives. The posed study area has adequate water supply for existing water uses, but planners want to evaluate best to increase water use while producing hydropower. These two objectives conflict.


### 3.2.1. Water Supply Objective

For an area relying on an aquifer and a stream for its water supply, the objective is to maximize the water provided from pumping wells and from surface diversions. Surface diversions include both diversions from streams and diversions from surface water reservoirs. Equation (2) is the objective function of maximizing water provided to users:(2)Z1=Max⁡(1T)∑k=1T[∑â=1NPCâ,kPpâ,k+∑ê=1NSCê,kSdê,kstr+∑î=1NRCî,kRdî,kres],
where Z1is the average flow rate delivered during all periods, pâ,k is the rate of water pumped from cell â during period k, dê,kstr is the rate of stream water diverted at diversion point ê during period k, and dî,kres is the rate of reservoir water diverted at diversion point î during time period k. NP, NS, and NR are total numbers of groundwater pumping, stream diversions, and reservoir diversions, respectively. Câ,kP, Cê,kD, and Cî,kR are weighting coefficients assigned to groundwater pumping, stream diversion, and reservoir diversion, respectively.Optimal solutions are sensitive to the weighting coefficients. Careful selection of coefficient values permits their use for economic optimization and emphasizing or de-emphasizing specific decision variables. A variable can be made ineffective in the objective function by setting its respective weighting coefficient equal to zero.


### 3.2.2. Hydropower Production Objective

Water released from the lake through a turbine produces hydropower as per the following:(3)Hydropower=epγQhp,k86,400hg,k,
where Hydropower is in kilowatts (kw), ep is the overall efficiency of the power plant (the product of hydraulic and turbine efficiencies), γ is the specific weight of water (9.81 KN/m3), Qhp,k is the reservoir release flow rate (m3/d), and hg is the gross head difference between the reservoir water surface and the point where water passes through the turbine (m).Rearranging yields (4), the objective function for maximizing the average hydropower produced during the management era:(4)Z2=Max⁡epγ×∑k=1T((Qhp,k86,400)(hlake,k-hturbine)(1T)),
where, for this site, ep is 0.65, hturbine is the elevation of the point where water passes through the turbine (20 m), and hlake,k is the lake stage at end of time period k.


### 3.2.3. Bounds and Constraints

Constraints enforce physical laws, legal rights, and management goals. The most simple constraints are lower and/or upper bounds on decision variables (Table4). Upper limits on well pumping rates assure that pumpage from a given well does not exceed the capacity of the well, pump, and pipe. Upper limits on reservoir and stream diversions prevent exceeding the capacity of the diversion canals. Lower limits on stream diversions consider legal rights. The upper bound on reservoir release (Qhp,k) equals the smaller of the turbine capacity and the penstock capacity. The lower bound reflects legal or economic reasons.Table 4
Bounds on decision variables in stress period 1.
Decision variableLower bound (L) (m3/day)Upper bound (U) (m3/day)Pumping well 13102820Pumping well 2308460Stream diversion 1758460Stream diversion 23602820Reservoir diversion282820Reservoir release57680195370More complicated constraints involving decision variables are demand constraints and constraints limiting total stream diversion. Demand constraints ensure that total water delivered to each user at least satisfies existing demand:(5)∑k=1T(pâ,k+dê,kstr+dî,kres)≥∑k=1TDemandk,
where Demandk is the total current demand for water in time period k. The intent is to possibly provide more water for future development.Upper bounds on cumulative stream diversions downstream of the reservoir ensure that stream diversions do not exceed the available water in the stream. In other words, the sum of stream diversions downstream of the reservoir must be always less than or equal the reservoir release to the stream.Each state variable that must be bounded during optimization is represented within the S-O model as an ANN. The GA represents their bounds using the penalty method. Table5 shows the imposed bounds on state variables.Table 5
Bounds on state variables.
ConstraintEquationConstraintStress period 11 (Lake stage)825.1 ≤ lake stage ≤ 29.3 m2 (Head 1 layer 1)626.2 ≤ head 1 layer 1 ≤ 30.2 m3 (Head 2 layer 1)622.9 ≤ head 2 layer 1 ≤ 29 m4 (Head 1 layer 2)620.7 ≤ head 1 layer 2 ≤ 29 m5 (Head 2 layer 2)617.1 ≤ head 2 layer 2 ≤ 27.1 m6 (Stream stage 1)720.5 ≤ stream stage 1 ≤ 20.6 m7 (Stream stage 2)720.2 ≤ stream stage 2 ≤ 20.3 mStress period 28 (Lake stage)821.9 ≤ lake stage ≤ 29.6 m9 (Head 1 layer 1)621.9 ≤ head 1 layer 1 ≤ 28 m10 (Head 2 layer 1)620.6 ≤ head 2 layer 1 ≤ 27.1 m11 (Head 1 layer 2)616.5 ≤ head 1 layer 2 ≤ 26.5 m12 (Head 1 layer 2)615.2 ≤ head 2 layer 2 ≤ 26.2 m13 (Stream stage 1)720 ≤ stream stage 1 ≤ 20.1 m14 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mStress period 315 (Lake stage)822.6 ≤ lake stage ≤ 28 m16 (Head 1 layer 1)618.9 ≤ head 1 layer 1 ≤ 27.4 m17 (Head 2 layer 1)618.9 ≤ head 2 layer 1 ≤ 25.6 m18 (Head 1 layer 2)614.3 ≤ head 1 layer 2 ≤ 25.3 m19 (Head 2 layer 2)614.6 ≤ head 2 layer 2 ≤ 24.4 m20 (Stream stage 1)719.9 ≤ stream stage 1 ≤ 20 m21 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mLimits on reservoir stage ensure that the water does not bypass the turbine and that the reservoir does not face environmental dangers:(6)hlake,kL≤hlake,k≤hlake,kU.
To prevent unacceptable drawdown, the hydraulic head in the aquifer at each well â should be within lower (L) and upper (U) bounds for each time period:(7)hâ,kL≤hâ,k≤hâ,kU.
Lower limit on stream stage at each control location û ensures adequate flow for recreational, navigational, legal, and environmental reasons:(8)stageû,k≥stageû,kL.


## 3.2.1. Water Supply Objective

For an area relying on an aquifer and a stream for its water supply, the objective is to maximize the water provided from pumping wells and from surface diversions. Surface diversions include both diversions from streams and diversions from surface water reservoirs. Equation (2) is the objective function of maximizing water provided to users:(2)Z1=Max⁡(1T)∑k=1T[∑â=1NPCâ,kPpâ,k+∑ê=1NSCê,kSdê,kstr+∑î=1NRCî,kRdî,kres],
where Z1is the average flow rate delivered during all periods, pâ,k is the rate of water pumped from cell â during period k, dê,kstr is the rate of stream water diverted at diversion point ê during period k, and dî,kres is the rate of reservoir water diverted at diversion point î during time period k. NP, NS, and NR are total numbers of groundwater pumping, stream diversions, and reservoir diversions, respectively. Câ,kP, Cê,kD, and Cî,kR are weighting coefficients assigned to groundwater pumping, stream diversion, and reservoir diversion, respectively.Optimal solutions are sensitive to the weighting coefficients. Careful selection of coefficient values permits their use for economic optimization and emphasizing or de-emphasizing specific decision variables. A variable can be made ineffective in the objective function by setting its respective weighting coefficient equal to zero.


## 3.2.2. Hydropower Production Objective

Water released from the lake through a turbine produces hydropower as per the following:(3)Hydropower=epγQhp,k86,400hg,k,
where Hydropower is in kilowatts (kw), ep is the overall efficiency of the power plant (the product of hydraulic and turbine efficiencies), γ is the specific weight of water (9.81 KN/m3), Qhp,k is the reservoir release flow rate (m3/d), and hg is the gross head difference between the reservoir water surface and the point where water passes through the turbine (m).Rearranging yields (4), the objective function for maximizing the average hydropower produced during the management era:(4)Z2=Max⁡epγ×∑k=1T((Qhp,k86,400)(hlake,k-hturbine)(1T)),
where, for this site, ep is 0.65, hturbine is the elevation of the point where water passes through the turbine (20 m), and hlake,k is the lake stage at end of time period k.


## 3.2.3. Bounds and Constraints

Constraints enforce physical laws, legal rights, and management goals. The most simple constraints are lower and/or upper bounds on decision variables (Table4). Upper limits on well pumping rates assure that pumpage from a given well does not exceed the capacity of the well, pump, and pipe. Upper limits on reservoir and stream diversions prevent exceeding the capacity of the diversion canals. Lower limits on stream diversions consider legal rights. The upper bound on reservoir release (Qhp,k) equals the smaller of the turbine capacity and the penstock capacity. The lower bound reflects legal or economic reasons.Table 4
Bounds on decision variables in stress period 1.
Decision variableLower bound (L) (m3/day)Upper bound (U) (m3/day)Pumping well 13102820Pumping well 2308460Stream diversion 1758460Stream diversion 23602820Reservoir diversion282820Reservoir release57680195370More complicated constraints involving decision variables are demand constraints and constraints limiting total stream diversion. Demand constraints ensure that total water delivered to each user at least satisfies existing demand:(5)∑k=1T(pâ,k+dê,kstr+dî,kres)≥∑k=1TDemandk,
where Demandk is the total current demand for water in time period k. The intent is to possibly provide more water for future development.Upper bounds on cumulative stream diversions downstream of the reservoir ensure that stream diversions do not exceed the available water in the stream. In other words, the sum of stream diversions downstream of the reservoir must be always less than or equal the reservoir release to the stream.Each state variable that must be bounded during optimization is represented within the S-O model as an ANN. The GA represents their bounds using the penalty method. Table5 shows the imposed bounds on state variables.Table 5
Bounds on state variables.
ConstraintEquationConstraintStress period 11 (Lake stage)825.1 ≤ lake stage ≤ 29.3 m2 (Head 1 layer 1)626.2 ≤ head 1 layer 1 ≤ 30.2 m3 (Head 2 layer 1)622.9 ≤ head 2 layer 1 ≤ 29 m4 (Head 1 layer 2)620.7 ≤ head 1 layer 2 ≤ 29 m5 (Head 2 layer 2)617.1 ≤ head 2 layer 2 ≤ 27.1 m6 (Stream stage 1)720.5 ≤ stream stage 1 ≤ 20.6 m7 (Stream stage 2)720.2 ≤ stream stage 2 ≤ 20.3 mStress period 28 (Lake stage)821.9 ≤ lake stage ≤ 29.6 m9 (Head 1 layer 1)621.9 ≤ head 1 layer 1 ≤ 28 m10 (Head 2 layer 1)620.6 ≤ head 2 layer 1 ≤ 27.1 m11 (Head 1 layer 2)616.5 ≤ head 1 layer 2 ≤ 26.5 m12 (Head 1 layer 2)615.2 ≤ head 2 layer 2 ≤ 26.2 m13 (Stream stage 1)720 ≤ stream stage 1 ≤ 20.1 m14 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mStress period 315 (Lake stage)822.6 ≤ lake stage ≤ 28 m16 (Head 1 layer 1)618.9 ≤ head 1 layer 1 ≤ 27.4 m17 (Head 2 layer 1)618.9 ≤ head 2 layer 1 ≤ 25.6 m18 (Head 1 layer 2)614.3 ≤ head 1 layer 2 ≤ 25.3 m19 (Head 2 layer 2)614.6 ≤ head 2 layer 2 ≤ 24.4 m20 (Stream stage 1)719.9 ≤ stream stage 1 ≤ 20 m21 (Stream stage 2)719.7 ≤ stream stage 2 ≤ 19.8 mLimits on reservoir stage ensure that the water does not bypass the turbine and that the reservoir does not face environmental dangers:(6)hlake,kL≤hlake,k≤hlake,kU.
To prevent unacceptable drawdown, the hydraulic head in the aquifer at each well â should be within lower (L) and upper (U) bounds for each time period:(7)hâ,kL≤hâ,k≤hâ,kU.
Lower limit on stream stage at each control location û ensures adequate flow for recreational, navigational, legal, and environmental reasons:(8)stageû,k≥stageû,kL.


## 3.3. Genetic Algorithm and Modeling Methodology

Genetic algorithms (GAs) can be powerful and robust optimizers for nonlinear problems, if using appropriate optimization parameters and enough simulations. A simple genetic algorithm (SGA) is a search procedure based on the mechanics of natural selection and genetics. It uses selection and recombination operations while searching a solution space to identify the best solution. SGAs have been generally applied to single-objective optimization problems.However, many optimization problems have multiple objectives. When optimizing with GAs, multiple objectives have usually been addressed by either combining all objectives into a scalar objective function (consisting of a sum of weighted multiple objectives) or by converting some objectives into constraints having thresholds and penalty functions. Other methods for handling multiple objectives includeε-constraint method, goal programming, and Multi objective genetic algorithms (MOGAs).The complexity of conjunctive use problems has led some researchers to use MOGAs. MOGAs also facilitate linking nonlinear optimizers and groundwater quantity and quality simulators. An efficient MOGAs is the nondominated sorting genetic algorithm (NSGA) [46]. NSGA uses a ranking selection method to emphasize current nondominated (Pareto optimal) points and a sharing function method to maintain diversity in the population.For simulation, the presented model uses coupled ANN and finite difference simulation to represent more detailed flow interactions, than other reported models for multiobjective reservoir-stream-aquifer system management. The S-O model uses NSGA for optimization. It includes objective or fitness functions, constraints, 21 linked ANNs, and the optimization algorithm. Using the optimizer requires setting five parameters values properly. These are mutation probability (Pmutate), cross-over probability (Pcross), niche size (σshare), population size, and maximum number of generations. σshare is dependent on the number of peaks in the solution space. We used a mutation probability of 0.01 and a cross-over rate of 0.6. The stopping criterion was 500 generations. The initial population size was 2758 individuals. The number of considered peaks was 100. These values were derived via sizing criteria suggested by Mahfoud [47] and sensitivity analysis. Continuous decision variables were represented by binary coding. NSGA selection of strategies employed the stochastic remainder proportionate method [46].Figure5 shows the simulation/optimization modeling procedure.(1)
Generate numerous water management strategies. Each strategy is a set of pumping, diversions, and reservoir release rates. Run MODLAKE for each strategy. Store variable system responses for all stress periods.(2)
Using strategies and results of step 1, train an ANN to enable it to predict state variable response to other sets of decision variables.(3)
Apply the S-O model to compute optimal strategies. The model calls the ANN trained in step 2.(4)
Confirm that output strategies indeed satisfy all constraints. This is necessary because the ANNs are not perfectly accurate. To do this, run MODLAKE for the optimal water management strategies and check the state variable results.Figure 5
S-O model processing.


## 4. Results and Discussion

The new S-O model is applied to the hypothetical study area discussed previously. The coupling of finite difference and neural networks to optimize all interactions of a reservoir-stream-aquifer system is unique. Here, the aim is to maximize water supply and hydropower production simultaneously. The first objective (2) maximizes average water provided from pumping wells, stream diversions, and reservoir diversion. The weighting coefficient for each component was 1.0, emphasizing all decision variables equally. The second objective function (4) maximizes hydropower production by maximizing reservoir releases thru the turbine.The two objectives are noncommensurate and conflicting. Maximizing provided water leads to maximizing water diverted from the reservoir. Therefore, the amount of water available for the turbine to produce hydropower decreases. This leads to a drop in hydropower production.Figure6 shows the trade-off (Pareto optimum) curve between the two objectives. Providing the most water produces the least average power and total energy.Figure 6
Average delivered water versus produced hydropower and energy.Figure5 suggests that, after optimization, one must confirm that the proposed model predicts system responses with acceptable accuracy. To do this, five Pareto optimal strategies were used as inputs to MODLAKE. All tested state variables satisfied the constraints listed in Table 5, meaning that optimal strategies are mathematically feasible and considered implementable. Table 6 contrasts ANN optimizer and MODLAKE outputs. Negative values in Table 6 indicate networks overestimating state variable values. The maximum absolute error for all networks and all tests is 0.5 m, corresponding to the head in the unconfined layer at well 1 during the third stress period. Since this state variable is constrained between 18.9 and 27.4 m, the error 0.5 m constitutes 5.8%. The slight lake stage inaccuracy affects hydropower production. The maximum error among all tests was 1.03 KW. For this demonstration, that is acceptable.Table 6
Residual errors (MODLAKE minus ANN-GA) for five tested strategies.
Test 1 (m)Test 2 (m)Test 3 (m)Test 4 (m)Test 5 (m)Stress period 1Lake stage−0.06−0.04−0.09−0.08−0.08Head 1 layer 1−0.13−0.1−0.13−0.11−0.14Head 2 layer 1−0.08−0.06−0.08−0.09−0.1Head 1 layer 20.20.20.20.130.22Head 2 layer 2−0.13−0.1−0.15−0.22−0.25Stream stage 10.00.00.00.00.0Stream stage 20.00.00.00.00.0Stress period 2Lake stage0.020.080.050.030.06Head 1 layer 1−0.010.08−0.040−0.04Head 2 layer 10.180.180.050.080.08Head 1 layer 2−0.080.1−0.03−0.1−0.05Head 2 layer 20.080.060.040.090.07Stream stage 10.020.030.020.020.02Stream stage 20.040.030.040.040.04Stress period 3Lake stage−0.020.050.030.020.04Head 1 layer 1−0.5−0.29−0.43−0.45−0.49Head 2 layer 1−0.150.090.00−0.030.08Head 1 layer 2−0.32−0.16−0.22−0.24−0.21Head 2 layer 2−0.42−0.35−0.36−0.35−0.35Stream stage 1−0.01−0.01−0.02−0.02−0.02Stream stage 20.010.030.030.020.03Results show the ability of the presented S-O model to find optimal solutions that will not violate constraints. This shows penalty function method effectiveness for handling constraints for multiobjective optimization.


## 5. Conclusion

This paper presents a simulation-optimization (S-O) model based on artificial neural networks and genetic algorithms for solving multiobjective nonlinear, conjunctive use problems. The model simultaneously addresses all significant flows and represents reservoir-stream-aquifer-diversion interactions in a more detailed manner than previously reported models. Using trained ANNs as state variable simulators is adequately accurate. The model robustly developed the tradeoff curve between the objectives of maximizing provided water and hydropower production. The S-O model identified optimal solutions that would not violate system constraints.



---

*Source: 910203-2012-04-29.xml*