Developing a Conjunctive Use Optimization Model for Allocating Surface and Subsurface Water in an Off-Stream Artificial Lake System

This work develops a rule curve-based conjunctive use management model for optimizing the operating rules for a lake–groundwater system with off-stream storage lakes. The proposed procedure is a simulation-optimization approach that embeds an Artificial Neural Network (ANN) instead of a groundwater numerical model into a genetic algorithm (GA). The direct physical exchange between lake water with groundwater is simulated using the ANN model, which is a reduced version of a full numerical model, MODFLOW with an LAK3 module. By applying the ANN model, the proposed procedure can reduce the computational burden that is induced by the nonlinear exchange. An operating rule-based optimal conjunctive use management model for the Gaopin Artificial lakes system in Taiwan was thus developed using the proposed framework. A set of optimal solutions involves rule curves and a discount ratio. Simulation results demonstrate that the embedded ANN model can accurately simulate the nonlinear exchange of a lake with groundwater. The embedded ANN model is less computationally complex than the numerical model. This work demonstrates a methodology for reducing the computational burden of the optimal conjunctive use management model that is associated with an internal nonlinear system by using the ANN reduced model. Specifically, the concept of, and results obtained using the developed operating rule-based model incorporating five artificial lakes and considering the nonlinear exchange of those lakes with the groundwater system provides a valuable practical reference for solving related conjunctive use problems.


Introduction
Rapid population growth and economic development increase the difficulty of managing water resources [1][2][3].The operation of conjunctive use of surface and subsurface water is one of the effective strategies for mitigating water scarcity [1,2,4].By coordinating surface and ground water resource systems, the operation of conjunctive use significantly enhances the reliability of the water supply and reduces the risk of water shortage [5].
Artificial lakes constructed on an unconfined aquifer with high permeability are replacing pumping wells for conjunctive use.In Taiwan, a plan to construct an off-stream artificial lake, called Gaopin Artificial Lake, has been proposed to deliver stable water resources over the last decade [6].Since the high permeability of the aquifer makes the quantity of water exchanged between the lake and groundwater system significant, the operation of Gaopin Artificial Lake is a unique form of conjunctive use.The advantage of the off-stream design of Gaopin Artificial Lake is that it prevents the problem of sedimentation associated with in-stream surface water storage [7].The lakes can be operated with less environmental impact than pumping wells and they provide water of higher quality [8] because sedimentation is less [7].Therefore, determining how to effectively integrate the off-channel system into a conjunctive use scheme is worth determining.
Conjunctive use management models (CUMM) have been extensively used in water resource planning to achieve optimal operating strategies [1,[9][10][11][12].CUMM is developed based on a simulation-optimization approach that allocates surface water and groundwater based on quantity and/or quality to satisfy one or more objectives while satisfying certain constraints [13].To operate a conjunctive use system accurately, a CUMM must capture the behavior of physically independent and/or integrated surface water (SW) and groundwater (GW) systems, such as a reservoir (or a lake) and its connected aquifer.In cases of conjunctive use that involve pumping wells, the difficulty associated with a CUMM is the simulation of groundwater behavior.Traditionally, two methods-the response matrix method and the embedded method-are used to simulate the relationship between pumping quantity and groundwater drawdown [1,14].
The quantity and direction of exchange between the Gaopin Artificial Lake and the groundwater system vary dynamically based on the lake stages and the groundwater level.Many numerical models have therefore been developed to simulate seepage fluxes, SW/GW water table fluctuation, groundwater recharge, and other variables, for managerial purposes [15][16][17][18].For real cases, these numerical models might be complex and consume high computational burden.In the framework of the simulation-optimization approach, a simulation model that is embedded into an optimization process is called iteratively.Therefore, if the simulation model is complex and causes a high computational burden, the computational expense associated with the embedded simulation model makes it impractical [19].
The artificial neural network (ANN) has been used as a surrogate for the computationally expensive embedded method and the linear assumption of the response matrix method for GW-related simulation [1,14,20,21].Two steps-a training mode and a predicting mode-must be implemented to construct an ANN model for simulating groundwater.First, in the training step (which is associated with a large computational burden), the ANN model can be trained individually before embedding in the simulation-optimization framework.Second, in the predicting step (which can be run with high efficiency), the well-trained ANN model is embedded into a simulation-optimization modeling framework for conjunctive use management planning.Instead of full function numerical modeling, the well-trained ANN model not only captures GW-related non-linearity accurately but also reduces the computational cost significantly [1,14].
This work develops a CUMM to explore the optimal management of surface and subsurface water in an off-stream lake system.To efficiently simulate lake-groundwater exchange, the proposed CUMM integrates a recursive ANN model into a genetic algorithm (GA)-based optimization process.The recursive ANN model has been used as a surrogate model instead of numerical modeling [22] to forecast fluctuations of the water levels of lakes and surrounding aquifers following the withdrawal of water from a lake system.The proposed CUMM was used to optimize the operation of Gaopin Artificial Lakes.The construction of Gaopin Artificial Lakes, an off-channel reservoir providing an alternative water supply, was intended to decrease the severity of groundwater overuse in southern Taiwan.The lake system was designed to be hydraulically linked to the enclosing aquifer through highly-permeable geological continuities [6].The water budget of the lake system is governed by the quantity of the intake water, the water supply, and the surface-subsurface exchange of the lake system.
Since the water supply device that is considered herein is an artificial lake system, the operating rule (which is extensively used in the field of reservoir operation) is included in the proposed CUMM herein.The GA-based CUMM can be used to optimize the operating strategy for the off-stream lake system.Therefore, the proposed CUMM extends the traditional water resource management model to be applicable to an integrated lake-groundwater system.

Study Area
The Gaopin Artificial Lakes in the Ping-Tung plain were constructed for a conjunctive use project that seeks to increase the quantity and stability of the water supply in southern Taiwan (Figure 1).Th project includes a system with five artificial lakes, a water intake weir, and a waterway that connects the water intake to the lake system.The area of the Ping-Tung plain is approximately 508 ha.Based on core drilling, the possible alluvial layer is more than 100 m thick.The groundwater stored in this area is plentiful, and the groundwater levels are 10 m below ground.The soil in the study area comprises mostly sand and a clay layer from 33 to 38.5 m below ground.Based on pumping tests in the area, the horizontal hydraulic conductivities of the sand layers vary from 1 ˆ10 ´4 to 5 ˆ10 ´3 m/s, and the specific yields are from 0.05 to 0.25.The designed depth of each lake in the lake system is 12 m.The elevations of each lakebed are lower than the mean groundwater level, but the normal water levels are higher.The lakebed material is sufficiently permeable to hydraulically link the lakes to the groundwater system.Therefore, the enclosing groundwater system supplies water to the surface water system when the water level of the lake system is lower than groundwater level.Therefore, the study area has potential for the establishment of an off-stream lake-groundwater system.The Gaopin Artificial Lakes in the Ping-Tung plain were constructed for a conjunctive use project that seeks to increase the quantity and stability of the water supply in southern Taiwan (Figure 1).The project includes a system with five artificial lakes, a water intake weir, and a waterway that connects the water intake to the lake system.The area of the Ping-Tung plain is approximately 508 ha.Based on core drilling, the possible alluvial layer is more than 100 m thick.The groundwater stored in this area is plentiful, and the groundwater levels are 10 m below ground.The soil in the study area comprises mostly sand and a clay layer from 33 to 38.5 m below ground.Based on pumping tests in the area, the horizontal hydraulic conductivities of the sand layers vary from 1 × 10 −4 to 5 × 10 −3 m/s, and the specific yields are from 0.05 to 0.25.The designed depth of each lake in the lake system is 12 m.The elevations of each lakebed are lower than the mean groundwater level, but the normal water levels are higher.The lakebed material is sufficiently permeable to hydraulically link the lakes to the groundwater system.Therefore, the enclosing groundwater system supplies water to the surface water system when the water level of the lake system is lower than groundwater level.Therefore, the study area has potential for the establishment of an off-stream lake-groundwater system.

Lake-Groundwater CUMM Model Formulation
This section presents the formulation of the conjunctive use management model (CUMM) for the operation of Gaopin Artificial Lake and the procedure for computing the optimal solution to obtain the optimal water supply strategy.The formulation is composed of an objective function and the constraints that are required.The procedure for solving the optimization problem is a simple GA framework, which comprises a surrogate ANN model to simulate the exchange quantity between lake and groundwater in the lake-groundwater system.

Model Formulation
The mathematical formulation of the CUMM comprises an objective function and a set of constraints that are given by Equations ( 1)-( 8).

Lake-Groundwater CUMM Model Formulation
This section presents the formulation of the conjunctive use management model (CUMM) for the operation of Gaopin Artificial Lake and the procedure for computing the optimal solution to obtain the optimal water supply strategy.The formulation is composed of an objective function and the constraints that are required.The procedure for solving the optimization problem is a simple GA framework, which comprises a surrogate ANN model to simulate the exchange quantity between lake and groundwater in the lake-groundwater system.

Model Formulation
The mathematical formulation of the CUMM comprises an objective function and a set of constraints that are given by Equations ( 1)- (8).
Equation ( 1) is the objective function whose goal is to minimize a ten-day shortage index (SI 10 ).SI 10 is a modified version of the shortage index (SI) that was proposed by the Hydrological Engineering Canter, U.S. Army Corps of Engineers.The SI 10 and the original SI have the same formulation, but differ in duration.In Equation (1), SD t is the water shortage in period t; D t denotes the target demand in period t, and N denotes the number of periods.For the SI 10 , the duration of each period is ten days and the number of periods in one year is 36.For the original SI, the counting period is one year.The counting period of the proposed SI 10 is ten days, which is much shorter than one year.The shorter period of the shortage index makes the index more sensitive to reflect the water shortage situation.Equation ( 1) is subject to the following constraints in the form of Equation (2) to Equation (8).
S t " g pL t q or L t " g ´1 pS t q where H t and H t+1 denote the hydraulic head of the groundwater system at time steps t and t + 1, respectively; R t represents the exchange flow between the lake system and the groundwater system at time step t; OV t represents the overflow quantity of the lakes at time step t; L t represents the water level in the lakes at time step t; I t represents the inflow quantity to the lakes at time step t; O t represents the outflow from the lakes to meet water demand at time step t; B t represents the groundwater head located at system boundary at time step t; S t and S t`1 respectively denote the lake storage at the time step t-st and the following time step; d is the discount ratio which can decrease the quantity of water supply during drought period and to minimize the peak quantity of water shortage; S rc τ is the storage volume under the rule curve at time step t; HL b and HL s are the elevations of the bottom and top of the lake, respectively; HL 1 , HL 2 , T 1 , T 2 , T 3 , and T 4 are the geometric parameters of the rule curve.
The Gaopin Artificial Lake system includes a lake system and a groundwater system (Figure 2).Equation ( 2) is a simulation model of the lake-groundwater system that is used to determine the groundwater level in the next time step (H t`1 ) and the exchange flow (R t ) from the predefined current groundwater level (H t ), the groundwater level located at system boundary (B t ), inflow volume (I t ), supply quantity (O t ) and lake level (L t ).Since the groundwater is shallow, the direction of exchange flow changes dynamically.When lake levels are above the groundwater level, the exchange mechanism of the artificial lake system is that of recharge; when lake levels are lower than the groundwater level, the exchange mechanism is that of discharge.Equation ( 3) is the mass balance equation of the lake storage.Figure 2 presents a diagram of the water balance of the Gaopin Artificial Lake System.Four major elements of flow quantities-inflow, overflow, supply, and exchange flow between the groundwater system and the artificial lake system-influence the water balance.
Equation ( 4) defines the relationship between the volume of water that is stored in the lake system ( ) and the lake level ( ).Equation ( 4), which is conceptually similar to the idea of a heightarea-volume curve in reservoir operations, is the transformation equation between lake volume and lake stage.
Equation ( 5) defines the quantity of water supplied from the lake system ( ).Equation ( 6) specifies that the shortage equals the gap between the target demand ( ) and the supply quantity ( ).
Equations ( 7) and ( 8) define the rule curve that is displayed in Figure 3 (which separates the lake into two operating zones).The six geometric parameters, , , , , , and are the decision variables in the CUMM.Since the time step is ten days and the number of time steps in one hydrological year is 36, the proposed rule curve (whose length is also 36) helps to determine the supplied quantity of water at each time step.and define the elevations of the rule curve in the flooding period and the drought period, respectively.The period between and is the flooding period, and that from to is the drought period.The volume of artificial lake above the rule curve is the fully supplied zone; otherwise, the volume below the rule curve is the partial supply zone, and , is the critical volume while lake stage is located at the elevation of the rule curve.If the total volume of current storage and inflow, , exceeds , defined by the rule curve, then the supply quantity equals the demand quantity in Equation ( 5); otherwise, the supply quantity equals the demand quantity multiplied by a discount ratio (d).If the hydrological condition is one of extreme drought, then a minimize operator is applied when the total volume is too small to provide the discounted quantity of water demand.Equation ( 3) is the mass balance equation of the lake storage.Figure 2 presents a diagram of the water balance of the Gaopin Artificial Lake System.Four major elements of flow quantities-inflow, overflow, supply, and exchange flow between the groundwater system and the artificial lake system-influence the water balance.
Equation ( 4) defines the relationship between the volume of water that is stored in the lake system (S t ) and the lake level (L t ).Equation ( 4), which is conceptually similar to the idea of a height-area-volume curve in reservoir operations, is the transformation equation between lake volume and lake stage.
Equation ( 5) defines the quantity of water supplied from the lake system (O t ).Equation ( 6) specifies that the shortage equals the gap between the target demand (D t ) and the supply quantity (O t ).
Equations ( 7) and ( 8) define the rule curve that is displayed in Figure 3 (which separates the lake into two operating zones).The six geometric parameters, HL 1 , HL 2 , T 1 , T 2 , T 3 , and T 4 are the decision variables in the CUMM.Since the time step is ten days and the number of time steps in one hydrological year is 36, the proposed rule curve (whose length is also 36) helps to determine the supplied quantity of water at each time step.HL 1 and HL 2 define the elevations of the rule curve in the flooding period and the drought period, respectively.The period between T 2 and T 3 is the flooding period, and that from T 4 to T 1 is the drought period.The volume of artificial lake above the rule curve is the fully supplied zone; otherwise, the volume below the rule curve is the partial supply zone, and S rc, t is the critical volume while lake stage is located at the elevation of the rule curve.Equation ( 3) is the mass balance equation of the lake storage.Figure 2 presents a diagram of the water balance of the Gaopin Artificial Lake System.Four major elements of flow quantities-inflow, overflow, supply, and exchange flow between the groundwater system and the artificial lake system-influence the water balance.
Equation ( 4) defines the relationship between the volume of water that is stored in the lake system ( ) and the lake level ( ).Equation ( 4), which is conceptually similar to the idea of a heightarea-volume curve in reservoir operations, is the transformation equation between lake volume and lake stage.
Equation ( 5) defines the quantity of water supplied from the lake system ( ).Equation ( 6) specifies that the shortage equals the gap between the target demand ( ) and the supply quantity ( ).Equations ( 7) and ( 8) define the rule curve that is displayed in Figure 3 (which separates the lake into two operating zones).The six geometric parameters, , , , , , and are the decision variables in the CUMM.Since the time step is ten days and the number of time steps in one hydrological year is 36, the proposed rule curve (whose length is also 36) helps to determine the supplied quantity of water at each time step.and define the elevations of the rule curve in the flooding period and the drought period, respectively.The period between and is the flooding period, and that from to is the drought period.The volume of artificial lake above the rule curve is the fully supplied zone; otherwise, the volume below the rule curve is the partial supply zone, and , is the critical volume while lake stage is located at the elevation of the rule curve.If the total volume of current storage and inflow, , exceeds , defined by the rule curve, then the supply quantity equals the demand quantity in Equation ( 5); otherwise, the supply quantity equals the demand quantity multiplied by a discount ratio (d).If the hydrological condition is one of extreme drought, then a minimize operator is applied when the total volume is too small to provide the discounted quantity of water demand.If the total volume of current storage and inflow, S t `It , exceeds S rc, t defined by the rule curve, then the supply quantity equals the demand quantity in Equation ( 5); otherwise, the supply quantity equals the demand quantity multiplied by a discount ratio (d).If the hydrological condition is one of extreme drought, then a minimize operator is applied when the total volume is too small to provide the discounted quantity of water demand.

Optimal Solution
A genetic algorithm, as presented in Figure 4, was proposed to compute the optimal solution to the CUMM, which is defined by Equations ( 1)- (8).The proposed algorithm is a simple GA, which has been extensively used in engineering and management studies.This section will only outline the key procedures that are related to Equation (1) to Equation (8).The decision variables, HL 1 , HL 2 , T 1 , T 2 , T 3 , and T 4 are encoded into each chromosome, and Equation (1) to Equation ( 8) are computed in the fitness evaluation procedure.Constraints on the decision variables, given by Equations ( 7) and ( 8), were directly considered in the encoding of the chromosomes.Equation ( 2) is the simulation model of the lake-groundwater system and is the most complex and computationally-intensive equation in the CUMM.To reduce the computational loading of the algorithm, an ANN model is proposed as a reduced model surrogate for the numerical model, MOFLOW with the LAK3 module.The following section will discuss the ANN model.After Equation ( 2) is solved, other constraints and the objective function can be solved.Chromosome reproduction, crossover, and mutation were implemented using a simple GA.

Optimal Solution
A genetic algorithm, as presented in Figure 4, was proposed to compute the optimal solution to the CUMM, which is defined by Equations ( 1)- (8).The proposed algorithm is a simple GA, which has been extensively used in engineering and management studies.This section will only outline the key procedures that are related to Equation (1) to Equation (8).The decision variables, , , , , , and are encoded into each chromosome, and Equation (1) to Equation ( 8) are computed in the fitness evaluation procedure.Constraints on the decision variables, given by Equations ( 7) and ( 8), were directly considered in the encoding of the chromosomes.Equation ( 2) is the simulation model of the lake-groundwater system and is the most complex and computationallyintensive equation in the CUMM.To reduce the computational loading of the algorithm, an ANN model is proposed as a reduced model surrogate for the numerical model, MOFLOW with the LAK3 module.The following section will discuss the ANN model.After Equation ( 2) is solved, other constraints and the objective function can be solved.Chromosome reproduction, crossover, and mutation were implemented using a simple GA.

Development of Lake-Groundwater ANN Model
Rather than directly embedding a numerical model into the optimization model to simulate the nonlinear exchange between lake and groundwater, a recursive ANN model was proposed as a surrogate for the full numerical model-MODFLOW with the LAK3 package.Since the ANN is a data-mining algorithm, data concerning the dynamic and nonlinear exchange between the lake and the aquifer are required.Restated, MOFLOW with LAK3 is used to generate very many simulation results, which are the data used to train the ANN model.Figure 5 shows the flowchart of the construction of the ANN model.
Three steps are conducted to construct a well-trained ANN model.First, a numerical model was developed using MODFLOW with the LAK3 package; it was then used to simulate the dynamic and nonlinear exchange under various hydrological conditions and supply quantities.These simulation

Development of Lake-Groundwater ANN Model
Rather than directly embedding a numerical model into the optimization model to simulate the nonlinear exchange between lake and groundwater, a recursive ANN model was proposed as a surrogate for the full numerical model-MODFLOW with the LAK3 package.Since the ANN is a data-mining algorithm, data concerning the dynamic and nonlinear exchange between the lake and the aquifer are required.Restated, MOFLOW with LAK3 is used to generate very many simulation results, which are the data used to train the ANN model.Figure 5 shows the flowchart of the construction of the ANN model.
were used to train the ANN model.To establish the reliability of the ANN model, another set of verification data was used.Third, the ANN model-instead of the traditional numerical model-was embedded into the optimization model.

Simulation of Lake-Groundwater System Using LAK3 Package
MODFLOW with the LAK3 package (Merritt and Konikow 2000) was used to capture the dynamic exchange between lake water and the surrounding groundwater system.The governing equations in the LAK3 package combine Equations ( 9) and (10).
where is the exchange quantity between the lake water and the aquifer; and are the stage of the lake and the aquifer head, respectively; is the surface area of the lake; and are the hydraulic conductivity and the thickness of the lakebed, respectively; , , , and are the rates of precipitation on the lake, evaporation from the lake, runoff to the lake and withdrawal of water from the lake, respectively, at time step t; , and , are the rates of inflow and outflow, respectively, from the streams.
Equation ( 9) is based on Darcy's Law and specifies the exchange quantity between the lake and the aquifer.Equation ( 8) specifies the water budget of the lake.
The gridded area in the simulation is 8000 × 8000 m and each grid cell is 100 × 100 m (Figure 6).Each sub-lake has a depth of 12 m; therefore, the parameters and are 0 and 12 m, respectively.The lake is initially full.After the quantities of inflow and lake water usage are set, the LAK3 package can compute the water budget of the lake system with a permeable lakebed, including groundwater levels ( ), lake stages ( ), and the exchange quantity between the lake and the aquifer ( ).
Based on the previous geological investigation, the driller's logs indicated that the mostly material is coarse sand and a clay layer from 33 to 38.5 m below surface ground.Based on pumping tests, the horizontal hydraulic conductivity of the unconfined aquifer is 5 10 m/s.
In MODFLOW, the hydrogeological parameters of the aquifer and the aquitard are assumed to be homogeneous.For the aquifer, the value of the horizontal hydraulic conductivity is assigned based on the result of pumping test and the ratio between vertical and horizontal hydraulic conductivity is 1 10 .For the aquitard, the horizontal hydraulic conductivity is assumed to be 5 10 m/s.The specific yield of each layer is 0.2.The west boundary is the Chishan River, which was simulated using Three steps are conducted to construct a well-trained ANN model.First, a numerical model was developed using MODFLOW with the LAK3 package; it was then used to simulate the dynamic and nonlinear exchange under various hydrological conditions and supply quantities.These simulation results are the training data and verification data constructing the ANN.Second, the training data were used to train the ANN model.To establish the reliability of the ANN model, another set of verification data was used.Third, the ANN model-instead of the traditional numerical model-was embedded into the optimization model.

Simulation of Lake-Groundwater System Using LAK3 Package
MODFLOW with the LAK3 package (Merritt and Konikow 2000) was used to capture the dynamic exchange between lake water and the surrounding groundwater system.The governing equations in the LAK3 package combine Equations ( 9) and (10).
h l,t`1 " h l,t `∆t ˆ´P t ´Et `Rn f , t ´Wt ´Qsp,t `Qsi, t ´Qso, t Ās where Q sp is the exchange quantity between the lake water and the aquifer; h l and h a are the stage of the lake and the aquifer head, respectively; A s is the surface area of the lake; K b and b are the hydraulic conductivity and the thickness of the lakebed, respectively; P t , E t , R n f , t and W t are the rates of precipitation on the lake, evaporation from the lake, runoff to the lake and withdrawal of water from the lake, respectively, at time step t; Q si, t and Q so, t are the rates of inflow and outflow, respectively, from the streams.Equation ( 9) is based on Darcy's Law and specifies the exchange quantity between the lake and the aquifer.Equation ( 8) specifies the water budget of the lake.
The gridded area in the simulation is 8000 ˆ8000 m 2 and each grid cell is 100 ˆ100 m 2 (Figure 6).Each sub-lake has a depth of 12 m; therefore, the parameters HL b and HL s are 0 and 12 m, respectively.The lake is initially full.After the quantities of inflow and lake water usage are set, the LAK3 package can compute the water budget of the lake system with a permeable lakebed, including groundwater levels (H t ), lake stages (L t ), and the exchange quantity between the lake and the aquifer (R t ).
the river package (RIV).The boundary conditions at the east, north and south boundaries are the general head boundary (GHB), and the associated head values are interpolated from observations from an observation network close to the study area.

Lake-Groundwater ANN Model
The input variables of the numerical simulation at each time step are the groundwater levels and the lake stages at time step t ( and ), the net inflow ( ), and the groundwater level located at the system boundary ( ).Lake-groundwater ANN models were trained and verified for the simulation of lakegroundwater exchanges in each sub-lake.The lake-groundwater ANN models comprise an input layer, one or two hidden layers, and an output layer (Figure 7).A hyperbolic tangent sigmoid transfer function from −1 to 1 was used to reflect positive and negative fluxes between a lake and a groundwater system.Furthermore, the number of hidden layers and corresponding neural sectors were obtained by trial-and-error.Figure 7 presents the output variables of ANN, which are and .To represent the recursive behavior of a numerical model (the ANN model), the water balance equation and the transformation equation between lake volume and lake stage, Equations ( 2)-( 4), are applied recursively.The state variables , , and are determined by a previous run of Equation (2) to Equation ( 4) and the variables , , and are defined previously.Based on the mass balance equation (Equation ( 3)) and the transformation equation between lake volume and lake stage (Equation ( 4)), the lake volume and lake stage at the next time step ( and ) can be determined.Based on the previous geological investigation, the driller's logs indicated that the mostly material is coarse sand and a clay layer from 33 to 38.5 m below surface ground.Based on pumping tests, the horizontal hydraulic conductivity of the unconfined aquifer is 5 ˆ10 ´4 m{s.
In MODFLOW, the hydrogeological parameters of the aquifer and the aquitard are assumed to be homogeneous.For the aquifer, the value of the horizontal hydraulic conductivity is assigned based on the result of pumping test and the ratio between vertical and horizontal hydraulic conductivity is 1  10 .For the aquitard, the horizontal hydraulic conductivity is assumed to be 5 ˆ10 ´9 m{s.The specific yield of each layer is 0.2.The west boundary is the Chishan River, which was simulated using the river package (RIV).The boundary conditions at the east, north and south boundaries are the general head boundary (GHB), and the associated head values are interpolated from observations from an observation network close to the study area.

Lake-Groundwater ANN Model
The input variables of the numerical simulation at each time step are the groundwater levels and the lake stages at time step t (H t and L t ), the net inflow (I t ´Ot ), and the groundwater level located at the system boundary (B t ).The output variables of the simulation model are the groundwater level at the next time step (H t`1 ), lake stages at the next time step (L t`1 ), and the exchange flux between the lake and the aquifer (R t ).MODFLOW with the LAK3 package determines the output variables H t`1 , L t`1 , and R t recursively.To obtain training data and verification data, the combination of MODFLOW and LAK3 runs under various hydrological conditions and supply strategies.The variables I t are evaluated using the flow records from 1984 to 2003, and the variables O t are generated randomly within a predefined range.Ultimately, 8640 pairs of training data and 2880 pairs of verification data are obtained.Lake-groundwater ANN models were trained and verified for the simulation of lake-groundwater exchanges in each sub-lake.The lake-groundwater ANN models comprise an input layer, one or two hidden layers, and an output layer (Figure 7).A hyperbolic tangent sigmoid transfer function from ´1 to 1 was used to reflect positive and negative fluxes between a lake and a groundwater system.Furthermore, the number of hidden layers and corresponding neural sectors were obtained by trial-and-error.Figure 7 presents the output variables of ANN, which are R t and H t`1 .

Results and Discussion
Four cases (Table 1) were considered to evaluate the efficiency of water supply strategies that were determined using the proposed lake-groundwater CUMM.The four cases were designed to compare water supply strategies with a focus on two main considerations-the different of operation efficiencies between the lake system with impermeable lakebeds and permeable lakebeds, and the optimization strategies that implement a single rule curve for all sub-lakes or various rule curves for each sub-lake.In Case 1, the lake system is impermeable, so the water supply depends only on surface water resources.In Cases 2-4, the lake system is permeable, so groundwater is an additional resource, which reduces the risk of water shortage.In Case 2, the supply system runs without using rule curves, so the supply system fully meets water demand until the system storage is empty.In Cases 3 and 4, the lake supply system uses optimal operating rules, which are obtained using the proposed CUMM.In Case 3, a set of rule curves is optimized for all sub-lakes.In Case 4, five sets of rule curves are optimized individually.
The total simulation horizon is ten years, and the time step interval is ten days.The water demand is 8.5 (approximately 0.73 million cubic meter per day), which are set in the model.Table 2 presents the setups in the four cases, and Table 3   To represent the recursive behavior of a numerical model (the ANN model), the water balance equation and the transformation equation between lake volume and lake stage, Equations ( 2)-( 4), are applied recursively.The state variables H t , L t , and R t are determined by a previous run of Equation (2) to Equation (4) and the variables I t , O t , and OV t are defined previously.Based on the mass balance equation (Equation (3)) and the transformation equation between lake volume and lake stage (Equation ( 4)), the lake volume and lake stage at the next time step (S t`1 and L t`1 ) can be determined.The calculated variables H t`1 , S t`1 , and L t`1 are then used as input variables to the ANN model at the next time step.

Results and Discussion
Four cases (Table 1) were considered to evaluate the efficiency of water supply strategies that were determined using the proposed lake-groundwater CUMM.The four cases were designed to compare water supply strategies with a focus on two main considerations-the different of operation efficiencies between the lake system with impermeable lakebeds and permeable lakebeds, and the optimization strategies that implement a single rule curve for all sub-lakes or various rule curves for each sub-lake.In Case 1, the lake system is impermeable, so the water supply depends only on surface water resources.In Cases 2-4, the lake system is permeable, so groundwater is an additional resource, which reduces the risk of water shortage.In Case 2, the supply system runs without using rule curves, so the supply system fully meets water demand until the system storage is empty.In Cases 3 and 4, the lake supply system uses optimal operating rules, which are obtained using the proposed CUMM.In Case 3, a set of rule curves is optimized for all sub-lakes.In Case 4, five sets of rule curves are optimized individually.
The total simulation horizon is ten years, and the time step interval is ten days.The water demand is 8.5 m 3 s (approximately 0.73 million cubic meter per day), which are set in the model.

Effectiveness of Lake-Groundwater CUMM
The above results reveal that the conjunctive use of surface and subsurface water from the groundwater-lake system with permeable lakebeds (Cases 2-4) outperforms the strategy in which surface water is only used from the lakes with impermeable lakebeds (Case 1).The values of SI 10 are much lower in cases with permeable lakebeds, and the total shortage times and shortage volumes are smaller.The comparison reveals that the exchange quantity can significantly increase the stability of the water supply.
Comparing cases with and without rule curves, the results of operations in Cases 3 and 4 are better than that of Case 2, except with respect to total shortage time.The shortage time in Cases 3 and 4 are 1780 and 1810 days, respectively, which greatly exceed the shortage time in Case 2. The shortage time is longer because the purpose of partial supply is to save water during periods of slight drought and to minimize the peak shortage during periods of extreme drought, while the strategy of partial supply extends the shortage time to minimize SI 10 .The evidence demonstrates that the proposed CUMM can allocate water resources in a more cost-effective manner than the operation without using a rule curve.
In Case 2, the recharge volume from the lake system to the aquifer is 168,671 m 3 .During the flood seasons, the surface water is plentiful, so the lake stages are often above the groundwater, making the exchange mechanism one of recharge; during the drought seasons, the surface water is scarce, so the lake stages are often lower than the groundwater, making the exchange mechanism one of discharge.Over ten years of operation, the accumulated quantity of the exchange is positive, indicating recharge.Hence, the permeable lakebed not only increases the stability of the entire supply system but also increases the recharge quantity.Since the permeable lakebed connects the lake body to the aquifer system, the fact that the aquifer system increases the entire storage volume of the lake system enables a win-win, decreases the value of SI 10 , and recharges the groundwater system.
In Case 3, one operation rule and one discount ratio were optimized for five sub-lakes.The optimal T 1 , T 2 , T 3 , and T 4 were the third, the tenth, the 21st, and the 23rd time steps in a year, so the curve from April to July is at HL 1 and from August to next January is at HL 2 .The optimal HL 1 and HL 2 are 2.32 and 11.48 m high, respectively.The optimal discount ratio is approximately 82%.Since the top of the lake is 12 m high, the full supply zone is from 12 to 11.48 m high from August to January, and the partial supply zone is from 11.48 to 0 m high.
In Case 4, five operation rule curves and five discount ratios are optimized for each sub-lake.Comparing the optimal HL 2 , values of HL 2 are 11.40, 11.80, 6.07, 6.46, and 5.76 m.For lakes A and B, the elevations HL 2 are similar to HL 2 in Case 3; however, the optimal HL 2 values for lakes C, D, and E are much lower.A lower elevation of the rule curve corresponds to a larger amount of supplied water.Based on the flow direction of groundwater, in the optimal solution, lakes A and B-which are located on the upper stream-supply less water than the other lakes during drought seasons.
A comparison among the four cases revealed that Case 4, with an SI 10 of 2.37, was the most efficient in the simulation of water resource management.In this case, the recharge volume from the lakes to the aquifer is negative, meaning that the net effect of the exchange is the consumption of water from the groundwater system.
The SI 10 -based objective function is sensitive for reflecting the severity of a water shortage when optimizing the rule curve and the associated demand-discount ratio [23].Based on the objective function, a low SI 10 is expected because the SI 10 -based CUMM not only moderates the peak water shortage but also spreads it out among time steps.As in previous studies [1,14,23,24], the GA-based CUMM can be used to solve the water resource management problem, accounting for the non-linear groundwater system.The hydrological condition directly influences the optimized rule curve and the demand-discount ratio.The shape of the rule curve might be different in an area with different hydrological conditions.

Model Calibration and Validation of Lake-Groundwater ANN Model
This work successfully developed a group of ANNs as surrogate for MODFLOW with LAK3 to elucidate the lake-groundwater exchanges of each sub-lake that are obtained from lake operations (Table 4).The ANN model is trained using output data of the MODFLOW with LAK3 model with a random series of water supplies from the sub-lakes.The lake water supply (O t ) of each sub-lake was set as a time series that was generated using a random variable ranging from 2 to 10 m 3 s .The inflow (I t ) to the lake system is the available flow of the Laonung River from 1984 to 2003.The net inflow (IO t ) to the lake system was calculated as the volume of inflow (I t ) minus the volume of water supplies (O t ), which was also given as a random time series.Lake water levels in all sub-lakes, as well as upstream and downstream groundwater levels near each sub-lake were verified.The verification of lake water levels in all sub-lakes revealed the accuracy of the mass balance that was obtained as a result of the integrated simulation using the ANN and lake storage routing.The maximum RMSE of the lake water level is 0.41 m, which was that at sub-lake E downstream of the lake-groundwater system (Table 4).Therefore, our lake-groundwater ANNs can capture the volume of flows between the surface and the subsurface system.The ANNs can also accurately simulate the fluctuation of surrounding groundwater levels.The maximum RMSE is obtained for the downstream groundwater level of sub-lake E and is 0.55 m (Table 4).Accordingly, the ANN can represent the variation of lake stages and groundwater levels that are caused by lake operations.The well-trained ANN can be a surrogate for MODFLOW.The concept of using ANN as a surrogate for MODFLOW has been successfully applied to CUMMs to reduce computational cost [1,14].Unlike the ANN developed by Chen et al. [1] or Peralta et al. [14], the ANN proposed herein captures not only the fluctuation of groundwater level but also the dynamics of groundwater exchange between surface and subsurface systems.Our ANN further considers the fluxes between target sub-lakes and their surrounding groundwater system in an ANN simulation.The accurate prediction of fluxes is important for the design of a flexible operation rule that captures the potential effects of the environment on groundwater.

Conclusions
This work develops an operating rule-based conjunctive use management model that reveals the optimal operation of surface and subsurface water in an off-stream storage system.The proposed lake-groundwater CUMM represents an improvement on the previous CUMM by considering non-linearity of the exchange between a surface and a sub-surface system in a computationally cost-effective manner.The CUMM was successfully applied to a water supply project that was supported by an off-stream lake system in Taiwan, called the Gaoping Artificial Lake.The GA-based CUMM reveals the ability of the GA algorithm to solve a water resource allocation problem in a way that takes a nonlinear lake-groundwater exchange into consideration.Moreover, the CUMM integrates a lake-groundwater ANN model as a surrogate simulator for MODFLOW into the simulation-optimization framework.The ANN can capture nonlinear exchange between surface water and subsurface water in an off-stream lake system with permeable lakebeds.Therefore, the inflow and outflow quantities of each sub-lake can be accurately simulated by combining the ANN sub-model with water balance equation.
The cases considered herein highlight the key effectiveness, such as computational efficiency, flexible operation rule, and systematic view of the environmental impact of the conjunctive operation on the groundwater system of the proposed CUMM.The ANN-based model can simulate lake-groundwater exchange faster than a finite-element based numerical model can.Based on the computational efficiency of the ANN model, the CUMM can optimize the operating rule curve, which provides flexibility in improving the efficiency of the conjunctive use operation of an off-stream lake system.The CUMM can provide a rule curve and an associated demand-discount ratio.Furthermore, the proposed CUMM can provide a systematic overview of the inflow, outflow, and exchange between surface water and groundwater in the lake system.In addition to providing details of the mass balance that identify recharge/seepage between the lake and the groundwater system, the CUMM can also be used to monitor the drawdown of surrounding groundwater to further assess its environmental effects that influence the operations.
The study cases herein confirmed the ability of the CUMM to integrate an off-channel system into a conjunctive use scheme.The CUMM helps to reveal how to effectively improve the water supply performance of the system in a manner that reduces environmental impacts below those of an impermeable reservoir.Hence, this work provides a useful tool for operating an off-stream artificial lake system and a valuable reference for future conjunctive use studies.

Figure 1 .
Figure 1.The plan view and observation of the lake system.

Figure 1 .
Figure 1.The plan view and observation of the lake system.

Figure 2 .
Figure 2. The diagram of water balance for the Gaopin Artificial Lake system.

Figure 3 .
Figure 3. Decision variables which define the buffer zone of the operating rule.

Figure 2 .
Figure 2. The diagram of water balance for the Gaopin Artificial Lake system.

Figure 2 .
Figure 2. The diagram of water balance for the Gaopin Artificial Lake system.

Figure 3 .
Figure 3. Decision variables which define the buffer zone of the operating rule.

Figure 3 .
Figure 3. Decision variables which define the buffer zone of the operating rule.

Figure 4 .
Figure 4. Genetic algorithm (GA)-based optimization algorithm of the optimal conjunctive use planning model.

Figure 4 .
Figure 4. Genetic algorithm (GA)-based optimization algorithm of the optimal conjunctive use planning model.

Figure 5 .
Figure 5.The flowchart of artificial neural network (ANN) model construction.

Figure 5 .
Figure 5.The flowchart of artificial neural network (ANN) model construction.

Figure 6 .
Figure 6.Finite-difference grid and boundary of MODFLOW.
The output variables of the simulation model are the groundwater level at the next time step ( ), lake stages at the next time step ( ), and the exchange flux between the lake and the aquifer ( ).MODFLOW with the LAK3 package determines the output variables , , and recursively.To obtain training data and verification data, the combination of MODFLOW and LAK3 runs under various hydrological conditions and supply strategies.The variables are evaluated using the flow records from 1984 to 2003, and the variables are generated randomly within a predefined range.Ultimately, 8640 pairs of training data and 2880 pairs of verification data are obtained.

Figure 6 .
Figure 6.Finite-difference grid and boundary of MODFLOW.

Figure 7 .
Figure 7. Integration of the ANN model and lake storage routing in the forward simulation.

Figure 7 .
Figure 7. Integration of the ANN model and lake storage routing in the forward simulation.

Table 1 .
presents optimal solutions in each case.The values in the four cases are 21.10, 5.44, 4.26, and 2.37, and the total shortage times are 2080, 930, 1780, and 1810, respectively.List and descriptions of study cases.

Table 2 .
The results of the four case studies.
Table 2 presents the setups in the four cases, and Table 3 presents optimal solutions in each case.The SI 10 values in the four cases are 21.10, 5.44, 4.26, and 2.37, and the total shortage times are 2080, 930, 1780, and 1810, respectively.

Table 1 .
List and descriptions of study cases.

Table 2 .
The results of the four case studies.

Table 3 .
List of optimal solutions.

Table 4 .
ANN structure and the verification performance for each sub-lake.