Numerical solution effort in dependence on the optimization horizon
1. Introduction
The sustainable management of the water resources and a safe supply of drinking water will play a key role for the development of the human prosperity in the following decades. The fast growth of many cities puts a large pressure on the local water resources, especially in regions with arid or semiarid climate. A research project at Fraunhofer IOSB and AST has aimed to investigate ways for economic and sustainable use of the available water resources in the region of the capital of China, Beijing [1].
A main issue of the project is to develop components for a modelbased decision support system (DSS), which will assist the local water authority in management, maintenance and extension of the water supply system at hand [2]. This paper deals with the derivation of suitable management strategies for a mid (till long) term horizon based on assumptions for future environmental and socioeconomic conditions, which are provided by other modules of the DSS. The general structure of the proposed optimal control DSS is shown in Fig. 1.
An overview about several DSS concepts and implementations is provided in [3, 4]. A challenge of the given problem is the large area of the water supply system. The water management has to consider the total water resources of five river basins with an area of 16,800 km² as well as large groundwater storage in the plain with an area of 6,300 km². The main portion of the annual precipitation (85%) in this semiarid region is falling from June to September leading to a highly uneven distribution throughout the year. The formerly abundant groundwater resources have been overexploited over the last decades resulting in a strong decline of the groundwater head (up to 40 m). Five reservoirs are important for the management of the surface water in the considered area, where the two largest account for about 90% of the total storage capacity of roughly 9x10^{9} m^{³}. The water is distributed to the customers using rivers and artificial transport ways (channels, pipes) of a total length of about 400 km.
A common approach for policy generation in this field of application is to use mathematical programming techniques based on a dynamic model of the essential elements of the water allocation and distribution system [3]. The great impact of the groundwater storage for the supply system at hand requires a more detailed description of the groundwater flow dynamics compared to other known DSS implementations. Therefore a 3D FiniteElement model of the plain region has been developed. However, a direct integration of this 3D model into an optimal control framework is not possible due to its computational costs. A trajectory based model reduction scheme is proposed, which guarantees a very fast response of the DSS in combination with a specially tailored nonlinear programming algorithm.
The chapter is organized as follows: In section 2 the essential models (surface and groundwater models) and the model reduction approach are described. While the formulation of the optimal control problem is subject of section 3, the numerical solution of the large scale structured nonlinear programming problem is described in section 4. First results of the optimal water management approach are presented in section 5.
2. Water allocation model
The water allocation model can be divided into the surface water model and the groundwater model. The parts of the water allocation model are described in the sequel.
2.1. Surface water model
The surface water model has to comprise all important elements for the allocation, storage and distribution of water within the considered region. The intended field of application of the decision support system under development embraces management and upgrading strategies for the mid and long term range. This implies a significant simplification for the process model, because the retention time along the different transport elements (river reaches, channel or pipelines) is less than the desired minimum time step for decision of one day. Therefore, a simple static approach for flow processes is sufficient and the use of sophisticated models for the dynamics of wave propagation (like e. g. SaintVenantEquations) with respect to control decisions is avoided. In this case, the flow characteristics are represented by simple lag elements of the first order combined with deadtime elements. Given
where T_{1} is the only parameter and indicates the time lag; t indicates the time. The following relationship represents the deadtime element:
In this case the parameter is the deadtime
The surface water system is described as a directed graph. The edges characterize the transport elements and introduce only a variable for the discharge. The nodes represent reservoirs, lakes, points of water supply or extraction and simple junction points. Every node constitutes a balance equation involving the edges linked with and possibly the storage volume. The sole nonlinearity results from modeling the evaporation from the water surface of the storages (volumeareacurve), which is described by a piecewise polynomial approach. At time step
where
Channels with a very low slope are modeled as water storage. The level dependent upper bound for the channel outflow is derived from a steady state levelflow relation like e.g. ChezyManning friction formula and is directly added as constraint to the optimization problem.
The Structure of the Beijing water supply system is shown in Fig. 2. Firstly, there are the four main reservoirs Miyun, Huairou, Baihebao and Guanting. Further sources are groundwater storages and water transfers. Secondly, there are the water transportation systems such as channels and rivers. Miyun reservoir and Huairou reservoir are connected to Beijing by the MiyunBeijing water diversion. In the simulation model the arrows describe hydraulic behavior of water flow. Baihebao and Guanting reservoir are connected by tunnel and river Guishui. From Guanting water runs into the Yongding river water diversion system to Beijing. Existing retention areas for flood control are also considered in the simulation model.
The surface water from channels and rivers is delivered to the customers in different ways, directly or through the surface waterworks. Groundwater is distributed to the customers through ground waterworks as well as motorpumped wells. Therefore, the waterworks build the third part of the Beijing water supply system. The last category is made up of all the customer groups (agriculture, industry, households and environment). To complete up the cycle, catchments area models are integrated in the system to take into account precipitation and evapotranspiration.
The surface water simulation model has been implemented in Matlab/Simulink using the toolbox “WaterLib” [5] and contains the most important elements of the drinking water supply system of Beijing. The simulation model was developed to reach a sufficient accuracy as well as a high simulation speed. A oneyear simulation is carried out in a simulation time of several seconds. This is a very important condition for using the model in the decision support system.
2.2. Groundwater water model
The most important water resource in the considered area is groundwater that is modeled by a dynamic spatially distributed finite element groundwater model. The governing equation for groundwater flow is Darcy's law [6] describing slow streams through unconfined aquifers. Combining Darcy's law with mass conservation yields the partial differential equation (4) which is a diffusion equation.
In (4) denotes
The partial differential equation (4) is an initialboundary value problem which has to be solved numerically for
The main task with respect to the groundwater model is the parameterization of the largescaled model covering an area of 6,300 km². On the one hand the time independent soil parameters k_{f}, S_{0} have to be estimated and generalized for the whole domain Ω by a (small) set of measured values. On the other hand the source / sink terms Q_{expl}, Q_{rech} have to be calculated time dependent. For these calculations time dependent maps of precipitation and water demand are needed. The water demand is splitted into the three user groups households, industry and agriculture (see [8] for details). This parameterization issue is supported by powerful geographical information systems (GIS).
2.2.1. Hydrogeological conditions and derivation of hydrogeological parameters
The groundwater model area is located at the northern part of the North China Plain (NCP), which is the largest alluvial plain of eastern Asia. The NCP is a basin with quaternary aged surficial deposits (loess, sand, gravel and boulder, silt and clay). According to the hydrogeological profiles of Beijing the quaternary system in this region is fairly complicated. A great variety of different sedimentary facies exists with different thicknesses ranging from several tens of meters around the piedmont area to 150  350 meters in the northern central part of the NCP [9]. Groundwater is exploited in the layers of quaternary deposits, i.e. in the loose stratum/porous aquifers with high to very high water storage capacities. From the Taihang Mountains in the west to east there are two main geomorphological units in the model area: the piedmont plain below the mountain escarpments and the flood plain. In the piedmont plain the aquifers structure is coarse and becomes finer from west to east. In the flood plain the structure of aquifer is fine with silt sand, clay and silt interlay and in areas of ancient rivers and paleochannels the aquifer is composed mainly of gravels and coarse sands with good permeability. Therefore the distribution of groundwater in the Beijing region is inhomogeneous. Regions of high abundance and high yielding porous groundwater aquifers are the piedmont plains and the northeastern districts of Miyun, Huairou and Shunyi whereas less yielding aquifers are found in the Yangqing and Tong districts. In the transition zone from the Taihang Mountains to the NCP the quarternary sediments with low thickness of e. g. some tens of meters are lying on the older rock formations of the regions.
In the mountainous districts unstable groundwater distributions were assumed in dependence on the form of the rocks with geological discontinuities (fractures, joints, dissolution features) and the groundwater flow. In the transition area from the Taihang and Yanshan Mountains to the NCP stratigraphic sequences of various ages ranging from archaean metamorphic rocks to quaternary are documented in the geological and hydrogeological maps. A detailed description of the geological and hydrogeological conditions can be found in [10].
On the base of a conceptual geological model and a structured horizontal (2D) groundwater model, a horizontal and vertical structured 3D  groundwater model was developed describing the saturated zone till approx. 200 m depth below ground surface (bgs.) in the area of the quaternary sediments of the NCP. In addition the borehole data from approx. 125 drillings situated in the model area were used in the groundwater model. Although a quite homogeneous distribution of the boreholes was given, one measurement point represents an area of about 50 km^{2} which is only a rare database for modelling subsurface conditions.
There can be found strong variations in structure and thickness of the loose stratum sediments in the model area. The evaluation of all data (borehole data, geological and hydrogeological maps and profiles, ground water levels from observation wells, literature etc.) shows that the large number of the water bearing layers can be summarised in up to three essential ground water aquifers according to present knowledge on regional level. These aquifer systems are from top to down:
Aquifer I: Shallow aquifer in approx. 5 m to 30 m depth bgs.
Aquifer II: Primary aquifer, till approx. 120 m depth bgs.
Aquifer III in the depth area of approx. 120/140 m to 200/260 m bgs
The aquifers are separated by less permeable layers or aquitards, above all fine sands, silts and clays. It can be assumed that the three essential aquifers are not completely independent from each other, i.e. a groundwater exchange takes place between them in a certain range. Where low permeable layers or aquitards are absent or have a low thickness two aquifers can form a hydraulic unity as in the area of piedmont plains. Thus in the piedmont plains only one porous aquifer between the unsaturated loess top set layers and the bedrock was assumed. In regions, in which the separating layers have bigger thickness and larger extension, local confined aquifers can appear. Because of morphology and evolution processes perched aquifers can appear within the loess deposits. All these local effects are summarised in the above mentioned three essential aquifers.
The piedmont areas of the Taihang Mountains and the Yanshan Mountains along the western boundary and the northern/northeastern boundary of the area are the areas where groundwater inflow into the plains contributes to the groundwater recharge of the confined and unconfined aquifers. Because of the multilayered geological structure of the loose stratum in the model area, consisting of loess, alluvial loess, sandgravelcobbleboulder sediments with more or less mighty clay and silt inclusions, the details of hydrogeological conditions are complicated. Therefore the used spatial distribution of the hydraulic conductivity and specific yield (both important for good model results) are adequate phenomenological descriptions of mean values and not a detailed representation of local conditions in reality.
The aquifer characteristics were determined mainly on the base of the interpretation and evaluation of the above mentioned borehole data. The borehole data represent the geological layers (boring logs) at single points. They show high variations from one point to another. In particular the hydrogeological parameters k_{f} and S_{0} were deduced from these borehole data due to the following procedure:
In a first step the local single point data have to be transformed in values representative for an associated area (mesoscale values) [11]. For this step the information and data from thematic maps (e.g. Beijing hydrogeological maps lithology map, water abundance map, etc.) were included. The validity range of these mesoscale values could be specified with the informations from a water abundance map representing the water yield and water storage capability. These data resulted primarily on measurements in water exploitation wells. With this approach mesoscale values (for k_{f} and S_{0}) and their spatial distribution could be determined.
In a second step the discrete mesoscale values were interpolated within the model area
In a third step the absolute values of the smoothed spatial distribution of k_{f} and S_{0} were adapted to fulfil the water budget of the model area.
On the basis of the derived values a fine tuning was realized in order to get minimal differences between calculated and measured groundwater surface map.
The k_{f} values of the aquifers range from 2.0x10^{3}m/s in the region of the piedmont plains and the alluvial fan plains to 0.1x10^{4}m/s in the flood plains. The loose stratum depositions can be classified according to German Institute for Standardization Guideline 18130 as permeable to very permeable.
In [12] a mean storativity in the range of 0.08 and 0.18 has been estimated. Due to hydrogeological investigations of borehole data a regionalisation of the hydrogeological parameters could be performed, yielding a mean storativity of 0.13. The values of S_{0} changes from 0.19 (piedmont plains) to 0.03 (flood plains). This spatial distribution of k_{f}  and S_{0}values was also a basis for the regionalisation of the inflow boundary conditions.
The inflow from the north and west into the model area is governed by the transition zone between the mountain terrain and the plain where a high conductivity can be assumed. Here the inflow consists of the surface water run off from the mountains that depends on the precipitation rate. But due to investigations even after a number of dry years the total inflow did not disappear, but decreased from a long term mean value of about 0.7x10^{9}m^{3}/a to 0.4x10^{9} m^{3}/a could be observed by the Chinese partners. This amount of water could not only result from rainfall runoff from the mountainous domain. Therefore it was assumed a split of the horizontal groundwater recharge from groundwater inflow into the model area (see Fig. 4). One part is a raindependent contribution which in 'normal years', with a mean precipitation rate of about 590 mm/a is in the range of 0.4x10^{9} m^{3}/a. This value is scaled in dependency of the mean precipitation rate of the current year. In a dry year with a precipitation rate of 380 mm/a, for instance, the rain dependent inflow to 0.26x10^{9} m^{3}/a. For this contribution the imagination is that a part of the precipitation of the mountain slopes infiltrate on the surface and percolate down to the bedrock basis. On the relatively impermeable bedrock the water flows as shallow groundwater aquifer in the loose stratum with low thickness into the model area. The loose stratum depositions in the piedmont plains consist of boulder, gravels and sand with inclusions of local loess and loessloam lenses. These loose stratum depositions are well permeable and this inflow contribution from the mountains is asssumed to be directly dependent on the precipitation.
The second contribution to the horizontal groundwater recharge from the mountainous area is of about 0.3x10^{9} m^{3}/a and it was implemented deeper. Here the underlying idea is that this component corresponds to the part of the precipitation which infiltrates in the mountainous area through clefts into the deeper rock formations. The groundwater flow system in the bedrock consists of macroscopic structures and cavities linked with each other, as for example fracture networks, faults, layer joints, dissolution features and conduits etc. The groundwater inflow depths were assumed according to the distribution of the water bearing carbonate rocks, sandstones and crystalline rocks in the hydrogeological map of Beijing. These socalled deep inflows are not dependent directly on the precipitation and change only in terms decades.
The total horizontal groundwater recharge (inflow) ranges from 0.55 to 0.75x10^{9} m^{3}/a. The quantitative split is to some extend arbitrary and based only on plausibility considerations since none of the components can be measured directly.
2.2.2. Derivation of timedependent groundwater recharge and exploitation data
In order to determine the timedependent groundwater recharge and exploitation we follow the subsequent procedure:
It starts with stating a long term mean value budget for the considered area for getting an idea which fluxes are in which order.
In a second step the budget data are regionalized by means of land use maps.
Finally the temporal distribution is taken into account by implementing the crop water need (agricultural water demand) and precipitation as dynamic input data such that in the end due to balancing all data, e.g. irrigation, evapotranspiration, etc, become time dependent quantities.
The first task in setting up models covering the water resources of an area of this size is to construct a water budget. The water budget is a theoretical device that supports structuring the water resource system and identifying the most important water fluxes. Here fluxes into and out of the system has to be collected as well as the water fluxes within the model area. The intension must be to realize the relation of water fluxes to each other, to quantify them, separating the more important from negligible water fluxes and to estimate the error that happens due to neglecting them. Since most of the quantities in the water budget are not independent from each other, the quantification of the water budget must be an iterative process. Fig. 5 illustrates the water budget of the region of Beijing. All the before mentioned fluxes are entried. The width of the arrows corresponds to the quantity of the fluxes.
Since the groundwater model is spatially distributed the input data for the groundwater recharge and the exploitation have to be spatially distributed as well. The above mentioned water budget can be regarded as a lumped parameter model for the whole region. The next step is to relate these data to locations by adapting the water balance with respect to specific regions.
An approach which is followed quite often is to regionalize by means of land use maps. The land use of the considered region is depicted in Fig. 6 showing eleven land use classes. These classes can be summed up to the following four classes:urban and paved areas,
agriculturally used and irrigated areas,
water areas and
noncultivated areas
Regionalization of the water budget means to adapt and evaluate the balance equation
for each of these classes. It denotes
2.2.3. Software realisation of the parameter calculation within the DSS
The DSS works on the base of socalled
The consistency of data can only be assured by the user. So if we use precipitation rates for groundwater model and for the surface water model for instance in general the information is provided in different data set and formats. In this case the user is responsible, that the data for the groundwater model match the data for the surface water model, i.e. that at same time periods in same regions the same precipitation rates are considered.
Nevertheless, the user can be supported by a parameterization software, as it was realized within the Scenario Wizard of the Beijing DSS to ensure consistency as good as possible. With respect to the groundwater model the scenario is complete and consistent only if the before mentioned data are available for the complete simulation horizon:
Initial conditions
Boundary conditions (well fields and inflow)
Spatially distributed groundwater recharge
Spatially distributed exploitation
In order to obtain a complete and consistent set of input data for the groundwater model the user has to pass through the groundwater panel of the DSS Scenario Wizard and a socalled groundwater project is created. A groundwater project is a part of the scenario containing all groundwater relevant data.
Within the GW panel the user has to execute five subpanels such that in the end at least a complete set of input data for the groundwater model is generated.
The first subpanel calculates the surface water recharge. For this a weighting map is required that determines how much precipitation becomes surface water in direct or indirect manner. On the other hand a reliable precipitation map for the entire model area is necessary. Since the precipitation is not constant over the year we also need temporal weights that define the temporal distribution of the precipitation rate. The Scenario wizard provides a number of precipitation maps of the past which can be also used for future scenarios. The surface water recharge is not a direct input data for the groundwater model but it is required to determine the distribution of exploitation and groundwater recharge in time and space.
In the second step the agricultural irrigation from groundwater is and the total spatially distributed exploitation from groundwater is derived. For these calculations a map of the agricultural used areas (irrigated areas) is needed as well as corresponding information about the crop water need/ agricultural water demand. In addition some information about the irrigation from surface water and waste water are requested. Since the agricultural water demand is not constant during the year a temporal distribution is needed as well.
In the third subpanel the minimal surface water inflow and the diffuse losses in urban areas is computed and therefore maps of water areas and urban areas are requested. The resulting information is incorporated into the calculation of the groundwater recharge from surface water areas and from urban areas.
In this subpanel all input data with respect to the well fields and the inflow parameters have to be defined and entered in to table which asks for yearly data of the exploitation rates from the different well fields and the upper and deeper inflow into the model area.
In the last subpanel the generation of the spatial and temporal distribution of the groundwater recharge data is performed. Here the user again has to provide a weighting map that determines how much of the supplied water becomes groundwater recharge and how much becomes evapotranspiration. The rest of the required data was already determined in the before described subpanels. This step finishes the creation of the groundwater project and all the generated data are documented in a specific groundwater scenario report which is saved within the corresponding project directory. After the simulation of a scenario with a groundwater simulation run a report about simulation results is saved in the project directory as well.
All computed data are saved as ASCII grid maps which can be read by any GIS system. Figure 7 shows the graphical user interface generating groundwater maps.
In order to gain the project data accessible to the FEFLOW simulation model the spatially distributed data have to be projected on the finite element mesh. But in this step no new information is gained therefore it is not a procedure which is important for the scenario generation but it is only a technical requirement.
2.3. Derivation of reduced groundwater model
For the optimization of the water allocation system the full 3D FEM model (with > 150,000 nodes) is not very suited due to the mentioned large computational time. As for the optimization task a prediction of the hydraulic head (groundwater level) at a set of representative and fixed points is sufficient, an inputoutput model (e.g. a linear state space model) with considerably smaller order n (e.g. n < 50) than the original FEM model has to be derived.
Methods for model reduction of those large scale systems have gained increasing importance in the last few years [13]. Two main classes of methods for model reduction can be identified, namely methods based on singular value decomposition (SVD) and Krylov based methods. SVD based methods are suited for linear systems and nonlinear systems of an order n < 500 (e. g. balanced truncation for linear systems, proper orthogonal decomposition (POD) for nonlinear systems). Most of these methods have favourable properties like global error bounds and preservation of stability [13]. Krylov based methods are numerically very efficient as only matrix multiplications and no matrix factorization or inversion are needed. Hence they are also suited for largescale systems. Unfortunately, global error bounds and preservation of stability cannot be guaranteed. Hence actual research is focused on the development of concepts which combine elements of SVD and Krylov based methods [13].
All of these approaches have in common that they aim to approximate the state vector
The basic idea of the proposed model reduction method is sketched in Fig. 8. We assume the existence of a reference scenario which means that the time dependent input parameters uref(t) of the FEM groundwater model (especially groundwater exploitation Qexpl and recharge Qrech) are determined for the whole optimization horizon. In practical cases these reference scenarios are mostly available or can be generated by plausible assumptions. Hence the task consists in the derivation of a model which approximates the behavior of the full FEM model in the case that the input parameters u differ from uref(t). This model is gained by identification techniques: Test signals (e.g. steps) are added to the reference input uref(t) (dimension p) and the corresponding deviations from the reference output yref(t) (dimension q) are identified. Doing this separately for every component of the input/output vectors u and y, we finally merge the (p x q) single inputsingle output (SISO) models to a multi inputmulti output (MIMO) model. For the groundwater model, the input parameters are e. g. cumulated (e.g. spatially integrated) exploitation of certain regions or cumulated exploitations of large well fields. The output parameters of the groundwater model are the hydraulic head at representative points (“observation wells”). In our application 13 input and 13 output parameters have been defined by the users: The inputs consist by 9 counties and 4 wellfields, the 13 output parameters are 12 observation wells and the mean hydraulic head of the whole area of the water supply system. As the slow stream groundwater flow can be interpreted as diffusion process (cf. equation (4)) only nearby located input and output parameters (e.g. regions/wellfields and the corresponding observation wells) have some correlation and a SISO model with these input/output combinations can be gained. Due to this physical reason the number of relevant SISO models is relatively small and hence the resulting MIMO model of relatively low dimension (n < 50) which is appropriate for the optimization problem. This proposed approach can be called trajectory and identification based model reduction. A similar approach (for a nonlinear large scale system) is found in [14].
3. Water resources management as optimal control problem
The water resources allocation problem is formulated as a discretetime optimal control problem:
subject to
The state variables
The process equations (8) consist of the balance equations of the storage nodes and the reduced groundwater model. The balance equations of the nonstorage nodes are formulated as general equality constraints (9). The objective function (6) contains the goals of the water management, which are primarily the fulfillment of the customer demand, the compliance with targets for the reservoir and groundwater storage volume and the delivery of water with respect to environmental purposes. Therefore quadratic terms are formulated, which penalize the deviations from desired values, like e.g. for the demand deficit of the demand node j:
where
The inequality constraints (10) follow from the technical capabilities of the water distribution system and rules to guarantee safe operation, which are simple bounds for the control variables:
as well as constraints for the reservoir volume
and the hydraulic head
With respect to the practical applicability selected parts of the inequality constraints (10) can be relaxed in order to avoid infeasible optimization problem with respect to unrealistic management demands.
4. Numerical solution of the optimal control problem
The optimal control problem is numerically solved as large scale structured nonlinear programming problem (NLP):
The optimization variables are the state and control variables of the several stages in time:
The state equations of the process model are incorporated as equality constraints:
The advantage of this problem formulation with
One alternative approach consists of eliminating the state variables from the optimization problem. The reduced dimension of the according nonlinear programming problem with (
For the numerical solution of large scale nonlinear programming problems interior point (IP) solver have become popular during the last years because of their superior behavior for NLPs with many inequality constraints. In this approach the objective function is expanded by adding barrier terms for the inequality constraints:
The solution of the original NLP (15) results from the subsequent solution of (17) with a decaying sequence of
A typical water management problem (horizon of five years, discretization of one month) has about 8000 optimization variables, 5500 equality constraints and 7200 inequality constraints. The numerical solution takes approximately 60 iterations and a calculation time of 10 seconds on an Intel Core 2 Duo CPU (2.5 GHz). Table 1 shows the linear dependency of the numerical solution effort from the number of time steps within the optimization horizon.








30days  7898  56  52  9.6 
10 days  23880  66  145  35.1 
5 days  47853  56  286  61.7 
5. First results of the optimal water management approach
The proposed concept for optimal water management is applied to the Beijing region. The region has precipitation, which varies geographically, seasonally and yearly. Eightyfive percent of rainfall falls between July and September. Groundwater is the most important source of water for the Beijing region, covers about 5070%. Beijing has suffered from over exploitation of this source over the years. Surface water supply in the Beijing region depends mainly on upstream inflows of the major river systems Chaobai, North Grand Canal and Yongding. Aside from problems such as excessive withdrawal and water quality deterioration of surface waters, the lack of regional coordination leads to issues such as uncoordinated withdrawals. Besides these problems, the water supply system is subject to other common problems, such as rapid population growth and urbanization, decentralized reservoir/groundwater management, changing attitude towards sustainability and attribution to greater attention of environmental issues.
The optimal water management system is evaluated for three scenarios based on the same set of input data, which are a combination of historical rainfall measurements and customer demands reflecting the predicted development of population size and economic growth. The objective function contains four quadratic terms in order to penalize deviations from the desired final water level of largest reservoir in the system (Miyun reservoir), from the desired final average hydraulic head of the groundwater storage as well as the deficit of the customer demand separated into two groups for household/industry and agriculture. The deficit of the delivered water must be less than 5 % for domestic/industrial clients and less than 25 % for agricultural clients. For the third scenario it is assumed that water can be transferred to the considered region up to an annual amount of 300 Mio. m³ starting from the third year within the optimization horizon.
The overall water demand exceeds noticeable the natural sources. The reduced groundwater model has been derived under the assumption that this overconsumption was covered by the groundwater storage. This leads to a strong reduction of the average hydraulic head of the groundwater storage of about 9 meters during 5 years (scenario 1, see Fig. 10). Using the initial state of the Miyun reservoir and the final value of the reference trajectory for the average groundwater head as target in the objective function, for the base scenario (scenario 1) there are only small deviations from these values in combination with a minor demand deficit observable. In the second scenario the increase of the target value for the final groundwater hydraulic head at 5 m and a corresponding shift of penalty coefficients in order to keep this value lead to a better spreading of the overdraft over the different storages as well as to the customers.
Fig. 9 shows the course of the water level for two reservoirs. Because of its capacity of 4.4x10^{9} m³ the Miyun reservoir plays an important role for the long term management of the overall system. As can be seen, the years with aboveaverage precipitation produce only a medium rise of the water level. The second scenario, which attempts to reduce the decline of the average hydraulic head of the groundwater storage, results in an increased release of this reservoir compared to the base scenario, which corresponds to a change of the final water level from 142 m above sea level to 137 m. In the third scenario a part of this release is replaced by water from outside of the considered region, which reduces the water level decrease at about 2 m. The second reservoir is situated at a channel from the Miyunreservoir to the city of Beijing and serves only as intermediate storage. The admissible range for water management is completely utilized by the optimal control approach. The shift in the management target causes a different operating strategy because water from the connected channel is taken to replace groundwater abstractions in the nearby regions. The change of the management target for scenario 2 and scenario 3 induce also an increase of the overall demand deficit from 0.1 % to up to 3.5 % (scenario 2), which is nearly 6.8x10^{8} m³ over the full horizon.
Fig. 10 shows the time plot of the mean groundwater level of the optimization scenarios. It is obvious that in scenarios 2 and 3 the aimed increase of the target value for the final groundwater hydraulic head at 5 m is nearly achieved. In Fig. 11 the impact of the different strategies to the exemplary input 1 (exploitation in a certain region) and exemplary output 5 (groundwater level at a defined observation point) can be studied. The exploitation is clearly decreased in scenarios 2 and 3 which correspond to an increase of the groundwater level at the observation point. Finally, from Fig. 12 it can be seen that the performance of the drastically reduced groundwater model is good, reflecting the fact that the original FEM model with more than 100,000 nodes has been reduced to a state space model with 36 states.
6. Conclusions
In this paper an optimal control approach as a component of a decision support system (DSS) for the management of the total water resources (surface water, groundwater, external water resources) in a fast developing region under a critical water shortage has been presented. It has been proven that in spite of the large area which has to be managed and the corresponding complex surface and groundwater models the optimization problem could be solved in an appropriate computation time (~ minutes). This could be achieved by a drastic reduction of the complex groundwater model to a state space model of relatively low dimension (n < 50). The user (i.e. water allocation decision maker) is enabled to select from a number of predefined performance criteria as well as to assign constraints to the elements of the water allocation system in order to specify the management targets according to his/her needs. The performance of the proposed concept is demonstrated by close to reality optimization scenarios, whereby the benefit of a new strategic channel has been investigated with a planning horizon of 5 years. Actually the developed DSS component is used in a first version by the decision makers. Future work will focus on the application and adaptation of the developed concept and software for the water resources management of further regions with critical water shortage.
Acknowledgements
This work was supported by the German Federal Ministry of Education and Research (BMBF).
References
 1.
Bernard T., Krol O., Linke H., Rauschenbach T. Optimal Management of Regional Water Supply Systems Using a Reduced FiniteElement. at Automatisierungstechnik, 2009; (57)12: 593600.  2.
Rauschenbach T., Steusloff H., Birkle M. Modeling of Beijing municipality’s water resources as basis for an optimal water allocation system. In: Proceedings of Beijing Conference on Sustainable Water Management. Beijing, China; 2007.  3.
Mays L. W. Water resources system management tools. New York: McGrawHill Companies; 2005.  4.
McKinney D. C. International Survey of Decision Support Systems for Integrated Water Management. Technical Report, IRG Project No. 1673000, Bucharest, 2004.  5.
Rauschenbach T.WaterLib 2.0  Simulation and Management of Surface WaterSystems. In: 2008 Promotion Conference on International Advanced Water Technology.Beijing, China, 2008.  6.
Bear J. Dynamic of Fluids in Porous Media, New York, Dover Publications; 1988.  7.
Software package FEFLOW: http://www.feflow.info/ (accessed 5 June 2012)  8.
Bernard T., Linke H., Krol O. A Concept for the long term Optimization of regional Water Supply Systems using a reduced Finite Element Groundwater Model, VDIBerichte 1980, Düsseldorf: VDIVerlag, 2007; p751762  9.
ChenJ.Y., TangC.Y., ShenY.J., SakuraY., KondonA.,ShimadaJ. Use of water balance calculationand tritium to examine the dropdown of groundwatertable in the piedmont of the north china plain (ncp). Environmental Geology, 44:564571, 2003.  10.
Z. Han. Groundwater resources protection andaquifer recovery in china. Environmental Geology,44:106111, 2003.  11.
Zippel M., Hannappel S. Ermittlung des Grundwasserdargebotes der Berliner Wasserwerke mittels regionaler numerischer Grundwasserströmungsmodelle. Grundwasser, 4:195207, 2008.  12.
FosterS., GardunoH., EvansR., OlsonD., TianY.,ZhangW., HanZ. Quarternary aquifer of thenorth china plain. Hydrogeology Journal, 12:8193,2004.  13.
Antoulas A. C. Approximation of LargeScale Dynamical Systems. Philadelphia: SIAM Press, 2005.  14.
Wolfrum P., Vargas A., Gallivan M., Allgöwer F.: Complexity reduction of a thin film deposition model using a trajectory based nonlinear model reduction technique. American Control Conference, Portland, USA, 2005.  15.
Wächter A., Biegler L. T. On the implementation of a primaldual interior point filter line search algorithm for largescale nonlinear programming. Mathematical Programming, 2006: 106(1):25–57.  16.
Franke R., Arnold E. The solver Omuses/HQP for structured largescale constrained optimization: algorithm, implementation and example application, Sixth SIAM Conference on Optimization, Atlanta, USA, 1999.  17.
Griewank A., Juedes D., Utke J. ADOLC: A package for the automatic differentiation of algorithms written in C/C++. ACM Trans. Math. Software, 1996: 22(2): 131167.