An object-oriented development environment to optimally design cyclic storage systems

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.

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).

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.

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 (SOURCE_i) to the target object (TARGET_j) 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 (⁠⁠).

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

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.

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).

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.

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.

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’.

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).

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.

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.

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.

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.

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.

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.

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.

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%).

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.

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).

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).

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).

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.

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.