Model Based Sustainable Management of Regional Water Supply Systems

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 semi-arid 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 model-based 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 socio-economic 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 semi-arid 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⁹ m^³. The water is distributed to the customers using rivers and artificial transport ways (channels, pipes) of a total length of about 400 km.

Figure 1.

Structure of the proposed decision support system (DSS) for the region of Beijing

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 Finite-Element 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 non-linear 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 non-linear 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. Saint-Venant-Equations) 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 dead-time elements. Given y and u as an output and input, respectively. The following mathematical relationship indicates a lag element of the first order.

$T_1\cdot \frac{dy}{dt}+y=u$E1

where T₁ is the only parameter and indicates the time lag; t indicates the time. The following relationship represents the dead-time element:

$y=u(t-T_t)$E2

In this case the parameter is the dead-time T_t, which is thus the measure of the time taken for water to flow in a conduit over a known distance.

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 (volume-area-curve), which is described by a piecewise polynomial approach. At time step k, the volume of a reservoir node evolves as follows:

$$V_j^{k+1}=V_j^k+\Delta t^k\left({\displaystyle \sum _{i\in {\rm E}(j)}{Q}_i^k+A_{O,j}^k{q}_{evpot,j}^k-q_{seep,j}^k+A_{O\max }{q}_{prec,j}^k}\right),$$E3

where $$A_{O,j}^k=f\left(V_j^k\right)$$, the storage volume is denoted by V, the discharge into and from the storage is denoted by Q and discrete time step is denoted by Δt. The total evaporation from the reservoir depends on the water surface A₀ and the potential evaporation q_evpot. q_seep denotes the seepage from the reservoir to the groundwater and q_prec specifies the precipitation.

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 level-flow relation like e.g. Chezy-Manning 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 Miyun-Beijing 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 motor-pumped 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.

Figure 2.

Structure of the Beijing water supply system

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 one-year 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.

$S_0\dot{h}-\nabla \cdot \left(k_f\nabla h\right)=Q_{rech}-Q_{\exp l}$E4

In (4) denotes h the hydraulic head (which corresponds to the groundwater level) and k_f the hydraulic conductivity that governs the hydrogeological properties of the soil. S₀ denotes the specific storage coefficient. The terms on the right hand side of (4) summarize all sources and sinks that coincide with the time dependent groundwater exploitation due to industry, households and agriculture (Q_expl) and recharge e. g. due to precipitation and irrigation (Q_rech) in Ω.

The partial differential equation (4) is an initial-boundary value problem which has to be solved numerically for h in the 3 dimensional model domain Ω. The groundwater model has been implemented using FEFLOW, which is a Finite Element (FEM) software specialized on subsurface flow [7]. The initial condition is h (Ω,t₀) (groundwater surface) at the initial time t0. The inflow/outflow is described by Dirichlet boundary conditions, i.e. h (∂Ω) at the boundary ∂Ω and by well boundary conditions, that define a particular volume rate into or out of Ω. The advantage of the latter one is that they are scalable. The 3D FEM model consists of more than 150,000 nodes, distributed on 25 layers (cf. Fig. 3). Huge computational costs result from this high resolution. The simulation of 5 years needs ~15 Minutes on an Intel Core 2 Duo CPU (2.5 GHz). Hence it is very time consuming to calculate optimal water allocation strategies with the 3D FEM groundwater model. This is the motivation for model reduction (see subsection 2.3).

The main task with respect to the groundwater model is the parameterization of the large-scaled model covering an area of 6,300 km². On the one hand the time independent soil parameters k_f, S₀ 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).

Figure 3.

Mesh of the 3D Finite Element groundwater model of the region of Beijing

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² 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 multi-layered geological structure of the loose stratum in the model area, consisting of loess, alluvial loess, sand-gravel-cobble-boulder 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₀ were deduced from these borehole data due to the following procedure:

1. In a first step the local single point data have to be transformed in values representative for an associated area (meso-scale 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 meso-scale 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 meso-scale values (for k_f -and S₀) and their spatial distribution could be determined.

2. In a second step the discrete meso-scale values were interpolated within the model area

3. In a third step the absolute values of the smoothed spatial distribution of k_f and S₀ 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^-3m/s in the region of the piedmont plains and the alluvial fan plains to 0.1x10^-4m/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₀ changes from 0.19 (piedmont plains) to 0.03 (flood plains). This spatial distribution of k_f - and S₀-values was also a basis for the regionalisation of the inflow boundary conditions.

Inflow 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⁹m³/a to 0.4x10⁹ m³/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 rain-dependent contribution which in 'normal years', with a mean precipitation rate of about 590 mm/a is in the range of 0.4x10⁹ m³/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⁹ m³/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.

Figure 4.

The groundwater inflow regime into the NCP

The second contribution to the horizontal groundwater recharge from the mountainous area is of about 0.3x10⁹ m³/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 so-called 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⁹ m³/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 time-dependent groundwater recharge and exploitation we follow the subsequent procedure:

1. It starts with stating a long term mean value budget for the considered area for getting an idea which fluxes are in which order.

2. In a second step the budget data are regionalized by means of land use maps.

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

• non-cultivated areas

Figure 5.

Water budget of the region of Beijing

Regionalization of the water budget means to adapt and evaluate the balance equation

$P=GWR+ET+SWR$E5

for each of these classes. It denotes P the precipitation rate, GWR is the groundwater recharge, ET the evapotranspiration SWR corresponds to the surface water runoff. There are some quantities which are relevant for every class, like infiltration and there are some which are specific for a certain group. And for other classes additional quantities have to be added, i.e. irrigation rates (IRR) for agriculturally used areas.

Figure 6.

Landuse map of the region of Beijing

2.2.3. Software realisation of the parameter calculation within the DSS

The DSS works on the base of so-called scenarios whereby a scenario can be understood as set of input data that represents either a real historical situation or an imaginable situation in future with respect to precipitation, exploitation, water use etc. This complete set of input data is applied to the hydrological model to obtain an answer to the question ‘What are the effects to the water allocation system’. The historical scenarios usually are used for the validation of the models whereas future scenarios should deliver guidelines for the best way to act in future under certain circumstances. In order to obtain reasonable results from the simulation the scenario has to be complete and consistent. Otherwise the models will gain unreasonable results. A scenario is complete if all data are available the hydrological model expects. In order to keep the system stable often default values are inserted if no data are assigned. But in this case the question is if the data are consistent.

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 so-called 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.

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

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

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

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

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

Figure 7.

Graphical user interface for the generation of input data like groundwater recharge maps

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 input-output 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 large-scale 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 x with respect to a performance criterion, e. g. minimize the deviation between original system and reduced system for a given test input. As a black box input-output model would be sufficient for our purposes there is no need to approximate the whole state vector x. Furthermore, the dimension n of a reduced model which approximates the whole state space vector x would be in most cases n > 100. With this dimension, for the given optimization problem the solution time would be unacceptably high (~ hours). Last but not least the use of the commercial software FEFLOW also prevents the application of e. g. a Krylov based method as no model representation (e. g. state space model) is provided by the software.

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 input-single output (SISO) models to a multi input-multi 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].

Figure 8.

Reduced groundwater model as a linear state space model in combination with a pre-simulated reference scenario.

3. Water resources management as optimal control problem

The water resources allocation problem is formulated as a discrete-time optimal control problem:

$$\underset{u^k\text{, k}=\text{1,}\dots \text{,K}}{\min }\left\{F\left(x^K\right)+{\displaystyle \sum _{k=0}^{K-1}{f}_0^k\left(x^k,u^k,z^k\right)}\right\}$$E6

subject to

$$x^0=x\left(t_0\right)$$E7

$$x^{k+1}=f^k\left(x^k,u^k,z^k\right)$$E8

$$h^k\left(x^k,u^k,z^k\right)=0$$E9

$$g^k\left(x^k,u^k,z^k\right)\le 0$$E10

The state variables x are the volume content of the reservoirs and channels with small slope and the states of the reduced groundwater model. Control variables u are the discharge of the transport elements as well as the water demand of the customers. The uncontrollable inputs z are the direct precipitation and the potential evaporation for the reservoirs and the flow of rivers entering the considered region, which is derived by means of rainfall-runoff-models. The number of time steps within the optimization horizon is denoted by K.

The process equations (8) consist of the balance equations of the storage nodes and the reduced groundwater model. The balance equations of the non-storage 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:

$$f_0^k\left(u_j^k\right)={\rho }_j^k\Delta t^k\frac{{\left(u_j^k-q_{dem,ref,j}^k\right)}^2}{q_{dem,ref,j}^k{}^2}$$E11

where q_dem,ref,j is the demand of and u_j is the discharge delivered to the customer. While this term applies for every time step within the optimization horizon, other terms are formulated only for the final point of the horizon, like e.g. for the desired volume content of the reservoirs.

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:

$$u_{\min }^k\le u^k\le u_{\max }^k$$E12

as well as constraints for the reservoir volume x_v:

$$x_{v,\min }^k\le x_v^k\le x_{v,\max }^k$$E13

and the hydraulic head h_hydr of the observation wells:

$$h_{hydr,\min }^k\le g^k\left(x_{gw}^k\right)\le h_{hydr,\max }^k$$E14

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 non-linear programming problem (NLP):

$$\underset{y}{\min }\left\{J\left(y\right)|h\left(y\right)=0;g\left(y\right)\le 0\right\}$$E15

The optimization variables are the state and control variables of the several stages in time:

$$y={\left[{\left(x^0\right)}^{\text{T}}{\left(u^0\right)}^{\text{T}}\dots {\left(x^{K-1}\right)}^{\text{T}}{\left(u^{K-1}\right)}^{\text{T}}{\left(x^K\right)}^T\right]}^T$$E16

The state equations of the process model are incorporated as equality constraints:

$$h\left(y\right)=\left[\begin{array}{c}{x}^0-x\left(t_0\right)\\ f^0\left(x^0,u^0,z^0\right)-x^1\\ ⋮\\ f^{K-1}\left(x^{K-1},u^{K-1},z^{K-1}\right)-x^K\end{array}\right]$$E17

The advantage of this problem formulation with (K(n+m)+n) optimization variables is the special sparsity structure with a block-diagonal Hessian-matrix and block-banded Jacobian matrices (k: number of time steps, n: number of state variables, m: number of control variables), which follows from the fact, that the process equations as coupling element between adjusting stages are linear in the state variables x^k+1. The numerical solution of this problem requires about (K(n+m)³) basic arithmetic operations.

One alternative approach consists of eliminating the state variables from the optimization problem. The reduced dimension of the according non-linear programming problem with (K m) optimization variables comes along with a loss of structure. The Hessian and Jacobian matrices are full. Because of the solution effort of order (K m)³ and the large amount of control variables (number of edges in the network description) this approach is not promising for this special field of application.

For the numerical solution of large scale non-linear 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:

$$\underset{y}{\min }\left\{J\left(y\right)+\mu {\displaystyle \sum _{\text{j}=\text{1}}^{{\text{n}}_{\text{g}}}\ln \left(-g_j\left(y\right)\right)}|h\left(y\right)=0\right\}$$E18

The solution of the original NLP (15) results from the subsequent solution of (17) with a decaying sequence of µ → 0. The identification of right active set with its combinatorial complexity is avoided. The state of the art non-linear interior point solver IPOPT is used for the application at hand [15]. The interface for multistage optimal control problems of the optimization solver HQP [16], which provides an efficient way for problem formulation along with routines for the derivation of $$\nabla J,\nabla h,\nabla g$$ by means of automatic differentiation [17], is used and coupled to IPOPT.

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.

optimization horizon	Numerical solution effort
optimization variables	Iterations	main storage [MB]	calculation time [s]
30days	7898	56	52	9.6
10 days	23880	66	145	35.1
5 days	47853	56	286	61.7

Table 1.

Numerical solution effort in dependence on the optimization horizon

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. Eighty-five percent of rainfall falls between July and September. Groundwater is the most important source of water for the Beijing region, covers about 50-70%. 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 over-consumption 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⁹ m³ the Miyun reservoir plays an important role for the long term management of the overall system. As can be seen, the years with above-average 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 Miyun-reservoir 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⁸ 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.

Figure 9.

Water level of the largest reservoir in the water distribution system (Miyun-reservoir) as well as of a reservoir with seasonal storage capacity (Huairou-reservoir).

Figure 10.

Time plot of the mean groundwater level in the optimization scenarios 1 – 3.

Figure 11.

Time plot of input 1 (exploitation in a certain region) and output 5 (groundwater level at a defined point) in the optimization scenarios 1 – 3.

Figure 12.

Time plot of output 5 of the full FEM model compared with the reduced model (optimization scenario 1).

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