Cyclic storage system (CSS) is defined as physically interconnected and operationally integrated surface water and groundwater subsystems with full direct interactions between the subsystems. Mathematical development and implementation of a CSS model is very complex and all previous works are fully case dependent with a minimum possibility of generalization. This article proposes an integrated development environment called CSSDev, which assists researchers to create and design object-oriented CSS models more easily. Using CSSDev, researchers may skip regeneration of repetitive simulation codes for common elements of a CSS. CSSDev employs NSGA-II to optimally select the design parameters of the models. Two objective functions of the optimization problem are system's total costs and total loss associated with the development alternatives. A real-world large-scale CSS has been modeled and optimized to illustrate the performance of CSSDev. The final Pareto-front is presented and two selected solutions from the set of optimal non-dominated ones are evaluated and discussed.

INTRODUCTION

Water planners often work on development and management of projects that integrate utilization of surface and groundwater, the process which is commonly referred to as conjunctive use of surface and groundwater (de Wrachien & Fasso 2002). Although the first attempts in joint operation of surface and groundwater systems goes back to the early 1960s (Buras 1963; Buras & Bear 1964; Burt 1964), the conjunctive use of surface and groundwater received growing attention from the late 1990s along with the introduction of integrated water resources management.

Based on type of groundwater simulation model, conjunctive use management models may also be classified as lumped or distributed systems. In lumped models, the groundwater system is often treated as a simple storage cell, similar to a surface reservoir. Complex hydraulic connection between surface and groundwater is disregarded and stream–aquifer interactions are not addressed. In a distributed groundwater management approach, a distributed simulation model is employed to evaluate the aquifer response to external stresses (excitations) in the domain. Distributed parameter models are normally used to increase the accuracy of predictions and to achieve a higher degree of spatial resolution. In distributed modeling approach, different simulation models with varying accuracies have been employed to address the response of the groundwater storage systems to excitations and their interactions with surface water bodies.

For a distributed modeling scheme, the simulation model is coupled with the optimization model either by ‘embedding method (EM)’ or by ‘unit response matrix (URM)’ method. In EM, finite difference or finite element approximations of the governing groundwater flow equation are directly used as a set of constraints in the management model. In embedding approach, drawdowns are calculated at many grid points where the information has no economic interest (Gupta et al. 1996). This approach is somewhat inefficient and has limited application because of dimensionality problems.

In URM, on the other hand, response of the aquifer system to unit perturbations at any selected points in the domain is evaluated using an offline approach. The URMs are generally preferred for unsteady-flow optimization because they use constraint equations to restrict system response only at user-specified locations and times (Peralta et al. 1995). Unit responses are assembled to form response matrices to be included in the management model as groundwater flow constraints. Mathematically, it is possible to embed a fully distributed groundwater simulator into the optimization module to form a complete, embedded, simulation-optimization model. However, its solution for a large-scale, real-world, non-linear, and non-convex system requires extensive computational cost. Although recent advances in computing methodologies, such as parallelization and response surface methods, have increased the promising opportunities for computationally expensive simulation-optimization frameworks (Alba & Troya 2002; Tang et al. 2007), for a long-term planning problem of cyclic storage system (CSS), the unit response method may provide reasonable results. Therefore, it was decided to replace the simulation model by the URM method. Alimohammadi et al. (2009) proposed a modified unit response method (MURM) with satisfactory performance to be used in unconfined aquifers under moderate drawdown (Jahanpour et al. 2013).

Addressing a special version of active long cycle conjunctive use system, the general concept of CSS was first introduced by Lettenmaier & Burges (1982). Alimohammadi et al. (2009) extended this definition to differentiate between conjunctive use and CSS. They defined a CSS as physically interconnected and operationally integrated surface water and groundwater subsystems with full direct interactions between the subsystems. Based on their new definition, surface and subsurface impoundment subsystems might be treated as competing and potentially interconnected parallel storage facilities that minimize most of the problems associated with large-scale surface impoundments. In a holistic view, CSS forms a network structure of all possible water transfers between system components (i.e., surface reservoir, river, aquifer, and demand area) taking into account natural interactions and operational policies imposed on the system. Considering the problem as an optimization model, the ecological, environmental, and hydrological constraints may ensure feasibility and sustainability of the resulting solutions.

The lumped modeling approach for optimum design of a large-scale CSS was carried out by Afshar et al. (2008). Alimohammadi et al. (2009) presented a distributed parameter approach for optimum planning of CSSs. Afshar et al. (2010) developed a hybrid two-stage genetic algorithm (GA)-linear programing algorithm to optimize design and operation of large-scale CSS as a single objective problem, not considering the multi-objective design and operation problem. In addition to untouched multi-objective CSS, formulation and implementation of a large-scale CSS as an optimization model is a challenge for researchers and practitioners. Researchers have independently developed models of various scales and resolutions using various programing platforms and/or development environments, often being obliged to rewrite simulation or optimization codes for similar conjunctive use problems. An integrated development environment (IDE) can assist researchers to create and design CSS models faster and more easily in a unique format. Using an IDE, researchers may skip regeneration of repetitive simulation codes for common elements of a CSS, such as surface reservoirs and wells. They may simply add an instance object and then modify its properties. Furthermore, such an IDE makes it easier to share, view, understand, and even develop CSS models already built by other researchers. To ensure integrity, which is a significant feature of sustainable development, such IDE should comprehend surface water, groundwater and the ecosystems through which they flow. Furthermore, the IDE should be constantly updated over time to include more environmental, hydrological, ecological, and economic features of CSS.

This article describes a prototype object-oriented IDE called CSSDev, which is developed to create, view, modify, share, and optimally design multi-objective CSSs. CSSDev has a simulation-optimization structure. CSSs with different scales may easily be generated, viewed, and modified by simulation module of CSSDev which is an object-oriented development environment with a graphic user interface. The object-oriented paradigm of CSSDev guarantees the transparency required to comprehend models that are already generated by others. Initially, based on the specific case study, user generates a CSS model in the simulation environment. Once the model is completely developed and all properties of its objects are assigned, the simulation model is prepared to be coupled with the optimization module. Subsequently, a multi-objective optimization module optimally designs the CSS design parameters such as reservoirs and water transfer systems capacity, providing multiple optimum alternative plans as a Pareto-front.

Object-oriented structure of the simulation model, including classes of objects, objects properties, cost terms, and loss functions are explained in the model development section. Next, the optimization module is explained. A real-world large-scale CSS is developed and optimized by CSSDev followed by related discussion. A conclusion is also presented at the end of this paper.

MODEL DEVELOPMENT

Definition of classes is one of the fundamental tenets of object-oriented paradigm. In CSSDev, six classes of objects (i.e., surface reservoir (dam), river, aquifer, well, demand area, and allocation) are defined, each modeling one of the identified physical component of a CSS (Table 1).

Table 1

Graphical presentations and abbreviation of CSSDev classes

Class of objects Notation Abbreviation 
Surface Reservoir (dam)  RES 
River  RIV 
Well  WELL 
Aquifer  AQU 
Demand area  DEM 
Allocation   
Class of objects Notation Abbreviation 
Surface Reservoir (dam)  RES 
River  RIV 
Well  WELL 
Aquifer  AQU 
Demand area  DEM 
Allocation   

To create a CSS model, several instance objects of the aforementioned classes should be added to the model manually. In this phase of model development, the user is asked to draw a layout plan of a water resources system using the graphical objects provided in CSSDev environment (Table 1). Once completed, the user may proceed to assign properties for each object of the model. Properties of model objects are explained in the following subsections. The collection of the objects as a whole forms a CSS, and objects sharing common properties are said to constitute a class. Every object inherits its modeling code from the class it belongs to, whereas certain properties of the object are modifiable and can be customized manually by the user. Classes of CSSDev are explained in the following subsections.

Allocation class

Although allocation objects should be added to a CSS model subsequent to all other objects, to better explain the development procedure of CSSDev, this class is discussed prior to others. Allocation objects determine water allocations between objects in the model during the planning time. Every allocation object is shown in the form of which stands for the water allocated from the source object (SOURCEi) to the target object (TARGETj) during period t. For example, stands for the amount of water allocated from river #2 to well #3 in the 10th period. When presented without the time indicator, is an array built of NT (model periods) members corresponding to the volumes of water transferred from the source to target objects during the entire planning horizon. There are seven predefined time-dependent water allocation systems in the context of CSSDev:

  • (1) release from reservoir to river ();

  • (2) release from reservoir to demand area ();

  • (3) transfer from reservoir to recharge well to artificially recharge the aquifer ();

  • (4) water transfer from aquifer to reservoir ();

  • (5) river diversion to recharge well to artificially recharge of the aquifer ();

  • (6) river diversion to demand area (); and

  • (7) water transfer from aquifer to demand area ().

User may manually add or remove allocation objects between source and target objects. Allocations are treated as decision variables, because their values are unknown to the user and have to be determined through an optimization procedure. Any feasible set of allocations value form a potential development scenario of the CSS. A multi-objective optimization algorithm is used to obtain development scenarios which minimize both system total present value costs (PVC) and the loss associated with system deficits. The multi-objective optimization problem can be formulated as: 
formula
1a
 
formula
1b
Objective functions defined by Equation (1a) and (1b) outline a set of conflicting objectives which may receive different priorities from the decision-makers and cases under consideration. Every optimum development scenario, which is considered as a solution of the optimization problem, determines the optimum design capacities for all components and water allocations to all elements of the system. For any development scenario, the present value of the system construction, PVCconst in billion rials (BR) and operational costs, PVCop (in BR), form the first objective function as the system PVC (in BR) 
formula
2

The second objective, Equation (1b), accounts for unsatisfied demand. This objective has been explained in the subsection ‘Demand area class’.

PVCconst is the sum of construction costs of the reservoirs and water transfer subsystems which will be used to allocate water from different sources to the sinks 
formula
3

NRES is the number of reservoirs in the system, is the construction cost of the ith reservoir (BR), NALL is the number of allocation objects, and is the construction cost for the ith allocation system in BR. Construction cost of any element is to be defined by the users as a function of its capacity. Without loss of generalities, construction cost of wells and pipelines from wells to demand areas are disregarded and only their operational costs (pumping costs) are included in the model. This study assumes that the demand area is very close to the pumping wells where the cost of the pipeline may not impose significant cost on the system.

System operational cost (PVCop) consists of present value of groundwater pumping cost (), groundwater recharge cost (), and operation, maintenance, and replacement costs of allocation systems and reservoirs. As shown in Equation (4), the groundwater pumping cost is a function of consumed energy (Basagaoglu & Marino 1999). 
formula
4
 
formula
5

is the operation cost of groundwater pumping (BR), NT is the number of operational seasons in the operating horizon, NK is the number of pumping wells, uelif is the energy required to pump a unit volume of water to a unit height (kWh), ucen is the unit cost of energy (BR per kWh), is the pump efficiency in well k, is the seasonal interest rate, is the initial groundwater level in the pumping well k (m), is the drawdown in the pumping well k at the beginning of period t (m), is the water pumped from well k during period t in million cubic meters (MCM), kqv is the conversion factor (discharge to volume), is the operation cost of groundwater recharge, is the number of recharging wells, is the unit recharging cost in the recharge well l in period t (BR per MCM), and is the total recharge to well l during period t (MCM).

Operation cost of the ith surface reservoir () is assumed as a predefined fraction of its construction cost (Equation (6)). Operation cost of the remaining five types of allocation subsystems (i.e., , , , , ) is calculated by Equation (7), in which is the unit operation cost coefficient of the allocation system (defined by the user) and is the allocation value in period t determined by the solution to the model. User assigns appropriate cost coefficients for all allocation objects in a model. 
formula
6
 
formula
7

Surface reservoir class

A surface reservoir object models a dam in a CSS. Minimum and maximum allowable storage of the ith surface reservoir in the model ( and , respectively) are user-defined constant properties. Other constant properties of reservoirs are the coefficients of surface-storage equations (linear functions). User should adjust these coefficients so that the linear relationship can properly describe the variation of reservoir's surface as a function. Time-dependent properties consist of evaporation rate per period and natural inflow to the reservoir, which should be assigned by the user for all periods. Equation (8) models storage variation of the ith surface reservoir during the planning time. 
formula
8
where and are storages of the ith surface reservoir in the periods t and t + 1, respectively, is the natural inflow to the reservoir in period t, is the volumetric evaporation loss, and is the sum of all allocations related to the reservoir. All terms in Equation (8) are in MCM. These allocations are consist of releases (to demand areas, recharge wells, rivers) and incoming (from aquifer) in period t. Evaporation loss in period t is estimated by multiplying evaporation rate of surface reservoir i in time step t (), in meters, by mean surface area of the reservoir for that period () in km2 
formula
9

River class

River flow is commonly regulated by surface reservoirs, defined by an allocation object . River flow may be diverted to demand areas and/or to wells to artificially recharge the aquifer . Outflow from the ith river reach in period t is estimated as (all terms are in MCM) 
formula
10
where ND is the number of demand areas and is the hydraulic interaction between the ith river and aquifer in period t (Equation (11)). To simplify estimation of the interaction a uniform rectangular channel is used to represent the river geometry (McDonald & Harbaugh 1988) 
formula
11
where 
formula
12
In Equation (11), (106m2/month) is the river conductance for the ith reach which is a function of its semi-pervious stream bed hydraulic conductivity (, 106m/month), its length , width and thickness (all in meters). and are the water level and elevation of semi-pervious stream bed elevation (in meters). It is assumed that a fraction () of total water delivered to a demand area will return to the river as return flow in period t (): 
formula
13
Constraint defined by Equation (14) limits the river outflow from a river () to a predefined maximum value () in all time periods. Constraint defined by Equation (15) controls the river flow to provide required environmental flow which is the quantity of water required to sustain freshwater and downstream ecosystems in all time periods. All volumes are in MCM. 
formula
14
 
formula
15

Well class

A well can be used to artificially recharge the aquifer and/or to pump water to demand areas or back to surface reservoirs. In conjunctive use systems, aquifer plays a significant role in supplying water to satisfy different demands. However, to maintain sustainability, aquifer storage should be preserved by limiting pumping and supervising groundwater level variations in wells during the entire planning horizon. In CSSDev, discharge and/or recharge from/to wells are controlled by constraints defined by Equations (16) and (17), where and are user-defined maximum allowable values (in MCM) of pumping and recharge from kth and to lth wells in period t, respectively 
formula
16
Drawdown below the minimum permissible piezometric levels, which may lead to increased pumping costs, land subsidence, infiltration of poor quality water, and drying up of springs and shallow wells (de Wrachien & Fasso 2002), is in contrast to the idea of sustainable development. Constraint (18) is imposed on the groundwater level fluctuation where and are user-defined minimum and maximum values (in meters) of water drawdown in well k in period t, respectively. 
formula
17
 
formula
18
MURM method (Alimohammadi et al. 2009) is employed to approximate groundwater level variations in the wells. Response matrices may be developed offline for all point, linear, and distributed excitations using a numerical flow model such as MODFLOW (McDonald & Harbaugh 1988) and be imported to CSSDev for each excitation element. Excitation components are categorized into point, linear, and distributed elements. Pumping and recharge wells are samples of point excitation. Interaction between aquifer and river reach may be treated as linear excitation, while rainfall on the aquifer and deep percolation from demand areas are considered as distributed excitation elements. In MURM method, the excitation (drawdown) at excited well k, in period n, is estimated as 
formula
19
where is drawdown in well k at the end of period n in meters. , , , and are modified unit response coefficients of excited well k for unit excitation from pumping, recharge, rainfall on aquifer and recharge, respectively, in the demand area in period t. In Equation (19), is the number of pumping wells, refers to the number of recharging wells , , , and are the associated excitations in period t in MCM, respectively.

Demand area class

A demand area may include domestic, agricultural, and industrial sectors. User determines time-dependent water demand values for each demand area. In the proposed CSS, demands can be met via the surface reservoir, aquifer pumping, and river diversion (Equation (20)) 
formula
20
where is total water supply to dth demand area in period t in MCM. NS, NW, and NRIV are number of surface reservoirs, wells, and river reaches which contribute to supplying water to the dth demand area, respectively, and all volumes are in MCM. For every development plan proposed by CSSDev, a loss function is calculated based on water deficit and/or surplus. Loss function of Equation (21) penalizes both insufficient and excessive water supply to demand areas while Equation (22) only considers deficits as losses and does not penalize excessive supplies. User may choose between Equations (21) and (22) to calculate the loss function. 
formula
21
 
formula
22

Stedinger (1978) believed it is unrealistic to penalize positive deviations (excessive water supply) from demand values, while Klemes (1978) considered a quadratic form of loss function in which positive deviations from demand values are penalized as well as negative deviations. Various studies have considered hydrological, economical, and even political measures in calculation of the loss function; however, the appropriate loss function is case dependent and may change from one management plan to another. In this study, the economic aspect of CSS development has been regarded as a separate objective and disregarded in the loss function.

Aquifer class

Conjunctive use relies on utilization of aquifer as a parallel storage to surface reservoirs to supply water demands in an efficient and more reliable manner. However, uncontrolled exploitation of any aquifer will most likely lead to persistent negative results such as continuous water-level drawdown, progressive water-quality deterioration, and gradual reduction of the aquifer storage capability. These irreversible damages jeopardize the ability of future generations to meet their own needs and thus, are in opposition to the concept of sustainability which is based on intergenerational equity (Loucks 1997). To avoid overexploitation of the aquifer, constraint number 23 ensures that total aquifer abstraction volume is not greater than the total recharge volume over the planning horizon 
formula
23
In addition to artificial recharge, aquifer is naturally recharged from river flow (Equation (11)), precipitation (Equation (24)), and deep percolation from irrigated area. Without loss of generality, the annual recharge resulting from precipitation is assumed to be a given fraction (Seep) of its total volume 
formula
24
where is the water recharged to the aquifer via precipitation in period t (MCM). User determines aquifer surface area ( in km2), precipitation heights ( in meters) and seepage ratio (). It is also assumed that a prespecified fraction of the total water delivered to demand areas will percolate into the aquifer ( in MCM). Through recharge wells, aquifer is artificially recharged with water allocated from surface reservoirs and rivers. Wells are also used to pump water from the aquifer to the demand areas or back to the reservoirs. As discussed, all natural and artificial recharges and extraction to/from the aquifer alter groundwater level in wells.

OPTIMIZATION

Having completed the simulation model, user may proceed to the optimization step. As discussed, allocation values over planning horizon are decision variables of the optimization problem. Both objective functions, namely the loss function and PVC, are functions of allocation values. Owing to the conflicting nature of the two objectives, there is certainly no one best plan that keeps both of them at their absolute optimum values. In other words, there will be no development plan of a CSS by which both objective functions are at their global minimal points. Dealing with problems associated with conflicting and non-commensurable objective functions, a multi-objective optimization methodology is preferred to a single objective one.

Multi-objective optimization does not yield a single optimal solution, but identifies a set of technologically efficient and acceptable trade-off plans (called Pareto-front) among conflicting goals and interests. All solutions in the Pareto-front are non-dominated. Solutions in the Pareto-front are not inferior to other solutions in both objectives, and, furthermore, these solutions are better than others in at least one objective. Depending on the type of optimization problem, different shapes of Pareto-front can be expected. In the current study, both objectives of the optimization problem will be minimized; consequently, the Pareto-front is expected to have a convex shape. CSSDev employs NSGA-II to achieve the Pareto-front of optimum design alternatives of a CSS model. NSGA-II is well known as a prevailing and robust multi-objective GA for optimization of non-linear and non-convex problems (Deb et al. 2002). Simulated binary crossover (SBX) method (Deb & Agrawal 1995) is used as the crossover operator. Evolutionary process is continued until a stopping criterion is met. Predefined stopping criteria are reaching the maximum number of generations and time limitation.

To discriminate a solution from a Pareto-front is a political decision (Loucks et al. 1981). Multi-objective analyses should assist those responsible for making these political decisions by illustrating the range of optimum decisions and the impacts of the alternative and competing plans. Following this idea, CSSDev only presents the Pareto-front of the optimum solution and leaves the selection of the final solution to users.

CASE STUDY

CSSDev was used to develop and optimize a real-world CSS model which has been previously solved by Alimohammadi et al. (2009) by a single objective optimization method. The study area is the Kineh-Vars reservoir (36 °7′18″N, 49 °4′13″E) and its irrigating area located at the downstream of the Abhar-rud watershed in Zanjan province in west-central Iran. The area consists of one surface reservoir, one river which is hydraulically connected to the aquifer, three supply and recharging wells, one demand area, and the underlying aquifer (Figure 1(a)). Figure 1(b) shows a graphic presentation of the CSSDev object-oriented model which was developed based on the features of the study area. All feasible allocation objects between model objects are assigned by default (Figure 1(b)). Having created all required model objects, user may proceed to assign their properties as discussed in the section ‘Model development’.

Figure 1

The study area.

Figure 1

The study area.

The planning horizon of this real-world, large-scale CSS was considered to be 10 years (from fall 1990 to summer 2000) containing 40 seasonal periods (NT = 40). Seasonally varying data are presented in Table 2. The model presented in this paper is deterministic; consequently, it accepts constant values for inflow to reservoir and demands. Users, though, may use historic information or other time series such as maximum or minimum flows or demand as input values; however, the model does not support stochastic calculations. The problem under consideration consists of six excitation points (three pumping and recharging wells), two surface excitations (deep percolation from the demand area and rainfall infiltration to the aquifer), and one linear excitation (river). Three pumping and recharging cells are considered as excited units. Response matrices have already been developed and verified for this case (Alimohammadi et al. 2009).

Table 2

Seasonally varying data

 Fall Winter Spring Summer Annual 
(MCM) 2.34 6.54 7.37 0.58 16.83 
2.85 5.96 33.59 1.16 43.56 
3.59 8.63 9.48 0.87 22.57 
5.15 15.7 15.32 0.94 37.11 
15.84 13.6 19.58 1.76 50.78 
4.96 11.75 49.2 1.45 67.36 
4.95 7.89 6.17 0.91 19.92 
3.96 10.41 17.85 0.93 33.15 
2.77 5.51 2.57 0.52 11.37 
1.84 5.17 5.8 0.24 13.05 
(mm) 185.5 77.1 348.5 717.1 1,328.2 
(MCM) 0.262 0.262 0.542 0.542 1.608 
(mm) 96.91 136.52 123.00 8.58 365.00 
(MCM/season) 2.886 1.456 11.024 10.634 26.000 
 Fall Winter Spring Summer Annual 
(MCM) 2.34 6.54 7.37 0.58 16.83 
2.85 5.96 33.59 1.16 43.56 
3.59 8.63 9.48 0.87 22.57 
5.15 15.7 15.32 0.94 37.11 
15.84 13.6 19.58 1.76 50.78 
4.96 11.75 49.2 1.45 67.36 
4.95 7.89 6.17 0.91 19.92 
3.96 10.41 17.85 0.93 33.15 
2.77 5.51 2.57 0.52 11.37 
1.84 5.17 5.8 0.24 13.05 
(mm) 185.5 77.1 348.5 717.1 1,328.2 
(MCM) 0.262 0.262 0.542 0.542 1.608 
(mm) 96.91 136.52 123.00 8.58 365.00 
(MCM/season) 2.886 1.456 11.024 10.634 26.000 

The demand area may receive water from the aquifer, the river and directly from the reservoir. It was assumed that 10% of the total water delivered to the demand area would percolate into the aquifer ( = 10%) and 10% would return to downstream of the river as irrigation return flow. Pumping wells are also considered as recharging wells with maximum allowed pumping and recharging rates of 3 MCM/season. Maximum and minimum allowed drawdown in wells, during any time step, are −10 and 10 m, respectfully. Initial drawdown of all wells is considered to be 10 m. Table 3 shows construction cost equations and operation cost coefficients for the case study.

Table 3

Construction cost equations and operation cost coefficients for the case study

  Construction cost equations (BR) Operation cost coefficient 
Reservoir   
Reservoir to demand area   
Reservoir to artificial recharge area   
Aquifer to reservoir   
River to demand area   
River to artificial recharge area   
  Construction cost equations (BR) Operation cost coefficient 
Reservoir   
Reservoir to demand area   
Reservoir to artificial recharge area   
Aquifer to reservoir   
River to demand area   
River to artificial recharge area   

RESULTS AND DISCUSSION

After developing the object-oriented model (Figure 1(b)) and assigning objects properties regarding data provided in Table 2, the optimization procedure triggered. Equation (21) was considered as the loss function of the optimization model. Employing CSSDev, the problem was solved for a planning horizon of 40 seasonal periods. With a population size of 200 chromosomes, NSGA-II algorithm stopped after nearly 20,000 generations which accounts for about 4 million function evaluation. Benefiting from Microsoft's PCP technology in the platform resulted in 30% reduction in computing time. Pareto-fronts of the last 800 generations and their diversities are presented in Figure 2.

Figure 2

Pareto-front of the last 800 generations.

Figure 2

Pareto-front of the last 800 generations.

Although it is for the last 800 generations, Figure 2 highlights the diversity of the solutions and the extent of the solution domain. It is shown that the non-dominated solutions form convex fronts and move gradually toward the left corner as solution proceeds. It clearly shows that the solutions with the first and second objective functions ranging (0–100) and (76–100) are included in the solutions, respectively.

Figure 3 shows the last Pareto-front with the set of optimal non-dominated solutions after 4 million function evaluation. The Pareto-front presented in Figure 3 consists of 52 non-dominated solutions. Each solution is an optimum design plan for the CSS under consideration. None of these solutions dominates others in the Pareto-front. Solutions with any degree of infeasibility are excluded from the Pareto-front; thus remaining solutions are all feasible alternatives for the CSS under consideration.

Figure 3

Final Pareto-front.

Figure 3

Final Pareto-front.

The conflicting nature of two objective functions may easily be observed in the pattern of the Pareto-front (Figure 3) where one objective function value decreases as the other one increases. Each solution in the Pareto-front is the most economic feasible design plan that can be proposed when a specific amount of loss function is accepted. For example, if the decision-makers accept the value of 17.68 for the loss function, then the most economical design is determined by the solution with total PVC of 77.78 BR.

Two solutions from the Pareto-front are selected to discuss in detail. The solutions, which are the first and last ones in the Pareto-front, are marked by circles in Figure 3. They are two extreme solutions, identified by solution number 1 and 2, with the least and most loss values, respectively. Solution number 1 is the most expensive (PVC = 84.69 BR) alternative with zero loss, whereas the solution number 2 has the lowest cost (PVC = 76.45 BR) and the highest loss value of 59.94 which corresponds to a deficit of 23% of total water demands during the planning horizon (i.e., 260 MCM). Table 4 presents resulting design capacities along with construction and operation costs of the selected solutions.

Table 4

Design capacities and costs of the objects

 Solution #1 Solution #2 
PVCconst(BR) PVCop(BR) Capacitya (MCM/Season) PVCconst(BR) PVCop(BR) Capacitya (MCM/Season) 
Reservoir 44.170 2.650 9.293 49.092 2.946 12.813 
Reservoir to demand area 0.0 0.0 0.0 17.133 0.007 4.853 
Reservoir to artificial recharge area 0.0 0.0 0.0 0.0 0.0 0.0 
Aquifer to reservoir 0.0 0.0 0.0 0.0 0.0 0.0 
River to demand area 20.484 5.172 9.357 1.05 0.337 0.402 
River to artificial recharge area 7.656 0.838 6.046 2.627 0.564 2.024 
Pumping wells – 1.629 N/A – 1.41 – 
Recharge wells – 2.095 N/A – 1.283 – 
Sum 72.310 12.384  69.902 6.547  
 84.694  76.449    
 Solution #1 Solution #2 
PVCconst(BR) PVCop(BR) Capacitya (MCM/Season) PVCconst(BR) PVCop(BR) Capacitya (MCM/Season) 
Reservoir 44.170 2.650 9.293 49.092 2.946 12.813 
Reservoir to demand area 0.0 0.0 0.0 17.133 0.007 4.853 
Reservoir to artificial recharge area 0.0 0.0 0.0 0.0 0.0 0.0 
Aquifer to reservoir 0.0 0.0 0.0 0.0 0.0 0.0 
River to demand area 20.484 5.172 9.357 1.05 0.337 0.402 
River to artificial recharge area 7.656 0.838 6.046 2.627 0.564 2.024 
Pumping wells – 1.629 N/A – 1.41 – 
Recharge wells – 2.095 N/A – 1.283 – 
Sum 72.310 12.384  69.902 6.547  
 84.694  76.449    

aReservoir capacity is in MCM.

To verify the results, law of conservation of mass was checked for results obtained from solution 1. Model subsystems (i.e., reservoir, aquifer, river, and demand area) were considered for the entire operation period of 40 seasons. Table 5 shows the cumulative values of system variables within the subsystems. Evaluation of mass balance of the transferred water within the subsystems indicated no apparent violation to the law of mass conservation.

Table 5

Mass balance verification in the subsystems (solution 1)

Reservoir Aquifer 
Initial storage 2.100 Total pumping from wells 116.143 
Final storage 0.000 Aquifer pumping to demand area 116.143 
Natural inflow to the reservoir 315.700 Aquifer pumping to reservoir 0.000 
Allocation from aquifer to reservoir 0.000 Seepage from river to aquifer 8.746 
Allocation to demand area 0.000 Total artificial recharge 66.892 
Allocation to wells 0.000 Percolation from precipitation 14.600 
Allocation to river 312.969 Percolation from demand area 26.000 
Evaporation 4.831   
Total reservoir inputs + initial storage 317.800 Total aquifer inputs 116.238 
Total reservoir outputs + final storage 317.800 Total aquifer outputs 116.143 
River Demand area 
Inflow from reservoir 312.969 Allocation from reservoir to demand area 0.000 
Allocation to demand area 143.857 River allocation to demand area 143.857 
Allocation to wells 66.892 Aquifer pumping to demand area 116.143 
Returned from demand area 26.000 Deficit 0.000 
Seepage from river to aquifer 8.746   
Outflow 119.473   
Total river inputs 338.969 Total water demand 260.000 
Total river outputs 338.969 Total water supply 260.000 
Reservoir Aquifer 
Initial storage 2.100 Total pumping from wells 116.143 
Final storage 0.000 Aquifer pumping to demand area 116.143 
Natural inflow to the reservoir 315.700 Aquifer pumping to reservoir 0.000 
Allocation from aquifer to reservoir 0.000 Seepage from river to aquifer 8.746 
Allocation to demand area 0.000 Total artificial recharge 66.892 
Allocation to wells 0.000 Percolation from precipitation 14.600 
Allocation to river 312.969 Percolation from demand area 26.000 
Evaporation 4.831   
Total reservoir inputs + initial storage 317.800 Total aquifer inputs 116.238 
Total reservoir outputs + final storage 317.800 Total aquifer outputs 116.143 
River Demand area 
Inflow from reservoir 312.969 Allocation from reservoir to demand area 0.000 
Allocation to demand area 143.857 River allocation to demand area 143.857 
Allocation to wells 66.892 Aquifer pumping to demand area 116.143 
Returned from demand area 26.000 Deficit 0.000 
Seepage from river to aquifer 8.746   
Outflow 119.473   
Total river inputs 338.969 Total water demand 260.000 
Total river outputs 338.969 Total water supply 260.000 

All values are given in MCM.

Both solutions have found it unjustified to transfer water directly from reservoir to artificial recharge area. Alternatively, that transfer has been made possible through the river to the recharging facilities. It means that the natural runoff, after being regulated in the surface reservoir, was released to the river to be diverted to artificial recharge area. Allocation from aquifer to reservoir also seems to be unjustified and deactivated in both solutions.

In addition to differences in total cost and total loss values, the selected solutions reveal some pronounced differences in the design objects and their capacities. The most obvious difference is observed in the allocation class (i.e., reservoir to the demand area). Specifically speaking, solution number 1 does not allow water allocation and transfer from the reservoir to demand area (zero capacity in Table 4), whereas solution number 2 has utilized this facility in its arrangement with a seasonal capacity of 4.853 MCM per season. In other words, selection of different solutions from the set of non-dominated optimum solutions may both affect the operational strategies and the system design characteristics. The non-linearity of the system is further illustrated by having larger reservoir capacity (12.813 MCM) for solution number 2 with larger total deficit compared to solution number 1 which has fully satisfied the demand on the entire planning horizon with smaller reservoir capacity (9.293 MCM). Regardless of these differences in the final design capacities, the results show how these small capacities have been efficiently utilized to avoid any unnecessary development in both solutions. Time variation of reservoir storage for both solutions is presented in Figure 4. As illustrated, on number of periods (mainly in spring season) both reservoirs are almost at their full capacities of 9.293 and 12.813 MCM.

Figure 4

RES1 storage variation, obtained by solutions 1 and 2.

Figure 4

RES1 storage variation, obtained by solutions 1 and 2.

Figures 5 and 6 compare the water supply to the demand area from various sources for the two selected solutions. As expected, the total water supply from the first solution with the highest total cost has fully satisfied the demand during the entire planning horizon with zero loss function (or deficit). The second solution, however, fails to satisfy the demand in some periods resulting in large loss function with almost 57% deficit for the entire horizon.

Figure 5

Supply to DEM1, obtained from solution 1.

Figure 5

Supply to DEM1, obtained from solution 1.

In solution number 1, 55% of the total demand is satisfied through the river diversion, and the remaining 45% is met by pumping water from the aquifer (Figure 5).

As shown in Figure 6, solution 2 supplies water to demand area from river (0.05%), from aquifer (42.5%), and from reservoir (42.5%).

Figure 6

Supply to DEM1, obtained from solution 2.

Figure 6

Supply to DEM1, obtained from solution 2.

Solution 2 shows a tendency to supply most of the demands from the reservoir (42.5%) and leave the river flow to recharge the aquifer or to deliver the environmental flow of the downstream.

There is a fairly equal tendency in both solutions to use aquifer to supply demands. None of the solutions was allowed to overexploit the aquifer due to the model constraints, which explicitly limit the groundwater use (Equations (16) and (23)) and fluctuation of water level in the wells (Equation (18)) to predefined values.

Figures 5 and 6 reveal that the aquifer plays a significant role in both solutions for supplying demands. This result may be used to emphasize the importance and performance of the parallel usage of aquifer, reservoir, and river in the proposed CSS. Although a cyclic system is more costly compared to a non-cyclic system, its performance and reliability in supplying the prespecified demand is significantly higher. The decision-maker is free to choose the most desirable solution considering the trade-off between the total cost and reliability of the system and/or its performance.

As illustrated in Table 5, approximately 58% of the aquifer recharge for solution number 1 is accomplished through artificial recharge while the next 42% comes from water seeping into the aquifer through precipitation (12%), deep percolation of irrigation water (22.5%), and the river–aquifer interaction (7.5%). Such a significant value discloses the importance of considering river–aquifer interactions in any conjunctive water resources management. Seepage from the river to the aquifer, as a natural phenomenon, may be considered as a free aquifer recharging plan, which decreases the need for other expensive artificial methods. This is particularly highlighted in the rivers with high hydraulic conductivity.

Solution 1 is the most expensive design alternative with zero deficit. As illustrated in Figure 5, seasonal demands for the entire planning horizon are fully satisfied using water from different sources with varying quantities. One may further appreciate the performance of the proposed CSS by pinpointing the allocations distribution from surface and groundwater sources during different seasons of the planning horizon. The third and fourth seasons are jointly satisfied from surface and groundwater sources with varying percentage. However, for the first and second seasons in the entire horizon quite contradicting strategies are observed. As an example, demand of the first season is completely satisfied with the surface water allocation whereas that of the 33rd season receives its water needs from the groundwater source. As another example, one may observe that only groundwater (surface water) is used to satisfy demand in the 14th (34th) seasons of the horizon. In addition, as demanded by Equation (15) (as a model constraint), ecological or environmental demand for the downstream river reach were met during the entire planning periods.

As an important model constraint, Equation (23) ensures that, at the end of the planning horizon, the total pumping from the aquifer should not exceed its total recharge. This is to ensure that any operation policy will not result in unwanted groundwater drawdown after the simulation period.

Seasonal variations of water allocated for artificial recharge from surface water, planned by solutions 1 and 2, are presented in Figures 7 and 8, respectively. As for solution 1, during the first 14 seasons the direct aquifer recharge remains very small with a jump starting from the 15th season where reservoir storage builds up and inflow to the reservoir increases. Although during periods 25–29, for example, the surface reservoir is not full, considerable amounts of already regulated water are transferred for artificial recharge. This is mainly for providing empty space for the relatively high flow coming into the reservoir in the next couple of seasons. It may be expected to have the highest recharge during the 23th period where the inflow is at its maximum; however, the high groundwater level in all three wells restricts further recharge which will violate the associated constraints.

Figure 7

RIV1 allocations to WELL1,2,3, obtained from solution 1.

Figure 7

RIV1 allocations to WELL1,2,3, obtained from solution 1.

Figure 8

RIV1 allocations to WELL1,2,3, obtained from solution 2.

Figure 8

RIV1 allocations to WELL1,2,3, obtained from solution 2.

As expected, seasonal allocation from different wells to the demand area follows a periodic trend. During the first 16 seasons, where the groundwater level is not high, groundwater extraction is comparatively lower than the next 20 periods (Figure 9).

Figure 9

Allocations from WELL1,2,3 to DEM1 , obtained from solution 1.

Figure 9

Allocations from WELL1,2,3 to DEM1 , obtained from solution 1.

The constraints which limit the groundwater pumping and recharge of each well to a predefined maximum value (3 MCM/season), have been observed (Figures 9 and 10 for solutions 1 and 2, respectively).

Figure 10

Allocations from WELL1,2,3 to DEM1 , obtained from solution 2.

Figure 10

Allocations from WELL1,2,3 to DEM1 , obtained from solution 2.

Also for both solutions, as shown in Figures 11 and 12, water levels fluctuate within the predefined range of [−10, 10 m], as constrained by Equation (18).

Figure 11

Water level fluctuation in WELL1,2,3, obtained from solution 1.

Figure 11

Water level fluctuation in WELL1,2,3, obtained from solution 1.

Figure 12

Water level fluctuation in WELL1,2,3, obtained from solution 2.

Figure 12

Water level fluctuation in WELL1,2,3, obtained from solution 2.

CSSDev was developed by Microsoft Visual Basic® 2010 as a standalone Windows application. It requires Microsoft® .Net Framework 4 to execute. CSSDev is equipped with Microsoft's Parallel Computing Platform (PCP) in order to take advantage of all cores of multi-core CPUs. PCP technology shows its power particularly while running complex algorithms which deal with large sets of data. As mentioned at the beginning of this section, the number of function evaluations during a run is significantly high (i.e., 4 million). This process took almost 6 hours and 45 minutes to complete without PCP enabled. A 30% decrease in time consumption was observed when activating PCP option of CSSDev and the optimization process finished in 4 hours and 45 minutes. Nevertheless, the PCP is an option to be enabled as needed. For larger and more detailed cases the PCP option might be a necessity.

CONCLUSIONS

Realizing the mathematical complexity of development and implementation of a CSS model, this article presented the prototype of CSSDev, which is an object-oriented development environment for CSSs. In the proposed platform, CSSs with different scales may easily be generated, viewed, and modified by built-in simulation classes of CSSDev with high transparency.

The prototype of CSSDev supports six built-in classes of objects which are surface storage, river, well, demand area, aquifer, and allocation. These classes can be used to develop object-oriented models in the graphical interface of CSSDev. It was shown that using CSSDev, researchers may skip regeneration of repetitive simulation codes for common elements of a CSS in both surface water and groundwater subsystems. The built-in multi-objective optimization algorithm (NSGA-II) provides the researchers and users with the final Pareto-front with set of non-dominated optimal solutions for further analysis to discriminate between the solutions and final decision-making.

Sustainability of the plans optimized by CSSDev was ensured imposing multiple constraints on ecological demands and conservation of aquifer storage in a multiple period real-world CSS. Using the built-in classes of the objects, such as aquifer, reservoir, and river classes, development of a new CSS model will be computationally efficient and technically simpler. However, implementation of the proposed object-oriented model requires some problem and site specific data and relationships. Some of them, namely coefficients of the response functions, calls for a validated groundwater model. These data and/or validated models for data generation for developing response functions' coefficients are common to all conjunctive use models.

REFERENCES

REFERENCES
Afshar
A.
Ostadrahimi
L.
Ardeshir
A.
Alimohammadi
S.
2008
Lumped approach to a multi-period–multi-reservoir cyclic storage system optimization
.
Water Resour. Manage.
22
(
12
),
1741
1760
.
Afshar
A.
Zahraei
A.
Mariño
M.
2010
Large-scale nonlinear conjunctive use optimization problem: decomposition algorithm
.
J. Water Resour. Plann. Manage.
136
(
1
),
59
71
.
Alimohammadi
S.
Afshar
A.
Marino
M. A.
2009
Cyclic storage systems optimization: semidistributed parameter approach
.
J. Am. Water Works Ass.
101
(
2
),
90
103
.
Basagaoglu
H.
Marino
M. A.
1999
Joint management of surface and ground water supplies
.
Ground Water
37
(
2
),
214
222
.
Buras
N.
1963
Conjunctive operation of dams and aquifers
.
J. Hydraul. Div.
89
(
6
),
111
131
.
Buras
N.
Bear
J.
1964
Optimal Utilization of a Coastal Aquifer
.
Sixth Congress of Agricultural Engineering
,
Lausanne
, pp.
238
248
.
Deb
K.
Agrawal
R. B.
1995
Simulated binary crossover for continuous search space
.
Complex Syst.
9
(
2
),
115
148
.
Deb
K.
Pratap
A.
Agarwal
S.
Meyarivan
T.
2002
A fast and elitist multiobjective genetic algorithm: NSGA-II
.
IEEE Trans. Evol. Comput.
6
(
2
),
182
197
.
Jahanpour
M. A.
Afshar
A.
Alimohammadi
S.
2013
Optimum management of cyclic storage systems: a simulation–optimization approach
.
J. Am. Water Works Ass.
105
(
11
),
E671
E683
.
Lettenmaier
D. P.
Burges
S. J.
1982
Cyclic storage: a preliminary assessment
.
Ground Water
20
(
3
),
278
288
.
Loucks
D. P.
1997
Quantifying trends in system sustainability
.
Hydrolog. Sci. J.
42
(
4
),
513
530
.
Loucks
D. P.
Stedinger
J. R.
Haith
D. A.
1981
Water Resource System Planning and Analysis
.
Prentice Hall
,
Upper Saddle River, NJ
.
McDonald
M. G.
Harbaugh
A.
1988
A Modular Three-Dimensional Finite-Difference Groundwater Flow Model. Techniques of Water-Resources Investigations
. US Geological Survey
,
Reston, VA
,
988
pp.
Peralta
R.
Cantiller
R.
Terry
J.
1995
Optimal large-scale conjunctive water-use planning: case study
.
J. Water Resour. Plann. Manage.
121
(
6
),
471
478
.