Using Stochastic Dynamic Programming to Support Water Resources Management in the Ziya River Basin, China

The optimal reservoir operation rules for the multiple major reservoirs of the Ziya River basin were found by using an SDP approach (as discussed in a later section). Including all five major reservoirs (Huangbizhuang and Gangnan reservoirs on the Hutuo River northwest of Shijiazhuang and Lincheng, Zhuzhuang, and Dongwushi reservoirs south of Shijiazhuang) as separate states into the model would make the model computationally intractable owing to the well-known curse of dimensionality. Therefore, the reservoirs were aggregated, assuming that the releases from all reservoirs can be diverted to any point in the downstream basin. This assumption is realistic because the primary rivers of the downstream basin are strongly linked with smaller canals. Conceptually, this single reservoir receives the combined runoff from the basin situated upstream the reservoirs. The water users of the upstream basin were divided into three groups; agricultural maize production, industry, and domestic. The water demand of these users can be satisfied with either river runoff or groundwater (no surface water storage).

Unused upstream runoff flows to the aggregated reservoir, where it is stored for the downstream water users, which are grouped into agricultural wheat-maize production, domestic, and industry. Reservoir releases generate hydropower up to a turbine capacity, and basin outlets flow to the Baiyangdian Lake and to the Bohai Gulf. The downstream users also have access to groundwater and water from the middle route of the SNWTP. Water from the reservoir can be diverted into the canal and sent to the city of Beijing. This setup makes it possible to test the impact of the new canal on optimal management of the reservoir.

The natural runoff to the reservoirs at the western boundary of the NCP was estimated from the seven subcatchments shown in Fig. 2. A simple conceptual water balance model based on the Budyko framework was used to simulate daily runoff (Zhang et al. 2008). The Hargreaves method was used to estimate daily evapotranspiration from daily maximum and minimum temperatures. The simple model is based on Budyko’s assumption that rainfall and available energy determine the long-term annual evapotranspiration from a catchment (Budyko 1958; Zhang et al. 2008). Budyko derived a water balance model (the Budyko curve) based on this assumption. This empirical relationship describes runoff as the fraction of precipitation not used to satisfy the evaporative demand (Zhang et al. 2008; Roderick and Farquhar 2011). Zhang et al. (2008) applied the Budyko framework for shorter time steps and allowed for water storage between time steps. A set of calibration parameters, including evapotranspiration efficiency and groundwater storage, controls the resulting discharge (Zhang et al. 2008).

The Yehe River tributary of the Hutuo River in the Taihang Mountains (gray shaded in Fig. 2) was used as a calibration subcatchment, with in situ measured discharge from the Pingshan station (Number 30912428) [Ministry of Water Resources (MWR) Bureau of Hydrology 2011]. This catchment has a relatively low population density and small reservoirs, so the flow regime was assumed to be close to natural. The runoff from this catchment was autocalibrated with daily measured discharge from 1971–1980, 1983–1991, and 2006–2008 [Ministry of Water Resources (MWR) Bureau of Hydrology 2011] with the objective of maximizing the monthly Nash-Sutcliffe efficiency (NSE). The resulting calibration parameters were applied to the other six mountainous subcatchments. The calibrated rainfall-runoff model was forced with 51 years (1958–2008) of daily precipitation and maximum and minimum temperature measurements from six weather stations (China Meteorological Administration 2009). Spatial interpolation of weather station data was performed with Thiessen polygons.

The combined simulated daily runoff from all upstream subcatchments (Fig. 2) was aggregated to monthly time steps and normalized. Each monthly discharge was assigned an inflow class [dry (0–20th percentile), normal (20th–80th percentile), or wet (80th–100th percentile)], and the transition probability matrix, ${\mathbf{P}}_{k,l}$, was determined by counting the number of transitions from each inflow class in month $k$ to each inflow class in month $l$. This aggregated data set was used to represent the stochastic properties of the natural reservoir inflow, and hence, used as input to the optimization model. The Markov chain was validated to ensure second-order stationarity (Loucks and van Beek 2005). The Markov chain was used to generate a synthetic runoff time series, which was tested against the measured runoff for stationary mean and variance.

For a system with $M$ users denoted with index $m$, the total costs, $tc$, of satisfying the demand of the users in each stage can be formulated as the cost, $c$, of allocating or curtailing ($ct$) the user $x$ amounts of water

In the Ziya River basin, the water sources are surface water ($sw$), groundwater ($gw$), and water from the SNWTP middle route ($sn$); the supply costs arise from pumping, treatment, and conveyance. Alternatively, the users can be curtailed, which leads to curtailment costs (as discussed in a later section). Although surface water is often associated with the lowest direct cost to the users, the absence of surface water will force the users to switch to curtailment or alternative, more expensive, sources. With scarce surface water resources and limited storage capacity combined with large seasonal and annual precipitation variations, present releases will decrease future surface water availability. Thus, present costs should be balanced against future uncertain costs. These costs are coupled in time because the reservoirs allow storage for future use (dynamic problem). Therefore, the one step ahead allocation problem for a given stage cannot be solved independently from the other stages. Incorporation of uncertain reservoir inflow adds stochasticity to the system.

The authors used SDP to develop rational reservoir operation rules in the uncertain environment, a method known in the hydropower sector as the water value method (Stage and Larsson 1961; Stedinger et al. 1984; Pereira et al. 1998; Loucks and van Beek 2005; Wolfgang et al. 2009). The backward recursive SDP equation runs in monthly time steps (stages) and calculates the minimum of the sum of immediate and expected future costs by using the classical Bellman formulation [Eq. (2)] for all discretized reservoir states (water levels). The objective is to minimize the total costs given by the optimal value function $F_t^{*}(V_t,Q_t^k)$ while satisfying water demands [Eqs. (3)–(7)]where $V_t$ = reservoir storage at stage $t$; $n$ = sources shown in Eq. (1); $p_{kl}$ = reservoir inflow serial correlation that describes the probability of a transition from a flow class $k$ reservoir inflow, $Q_t^k$, in stage $t$ to a flow class $l$ reservoir inflow, $Q_{t+1}^l$, in stage $t+1$. Because the future cost function is convex, future costs can be constrained from below by a set of linear constraints derived from the discrete sampling of total costs, $F_{t+1}^{*}$, as described by Pereira and Pinto (1991). The optimal solution will have a future cost that falls on one of these linear constraints. The set of linear constraints on the future cost is integrated with the remainder of the optimization problem to form one single linear program for the one step ahead optimization problem [Eq. (8)]. Three flow classes were defined for every month, as described previously, owing to the length of the simulated runoff series. The benefits from hydropower production (reservoir releases $r$ times marginal production benefit, $b_{hp}$) are subtracted from the immediate costs. The optimization problem is defined by the objective function [Eq. (2)] subject to constraints on water demand fulfillment [Eq. (3)], water balance of the reservoir [Eq. (4)], water balance of reservoir releases [Eq. (5)], upstream releases [Eq. (6)], and releases to Beijing [Eq. (7)]where $s{w}_{m,t}$, $g{w}_{m,t}$, and $s{n}_{m,t}$ = surface water, groundwater, and SNWTP water allocations, respectively; $c{t}_{m,t}$ = water curtailment of user $m$ in stage $t$; and $d_{m,t}$ = given demand of user $m$. The dynamics of the system is described by a water balance equation [Eq. (4)] linking reservoir storage in $t+1$ to reservoir storage in $t$, where $\Sigma s{w}_{u,t}$ is the sum of surface water allocations in stage $t$ to the upstream user, $u$ ($U$ users). The total reservoir releases consist of release through the hydropower turbines, $r_t$, and spills, $s_t$, exceeding the turbine capacity. The combined reservoir releases $(r_t+s_t)$ must equal the combined surface water allocations to the $D$ downstream users (denoted $d$) and unused outflow, $Q_{\text{out}}$, from the basin [Eq. (5)]. The upstream water users are constrained by the upstream runoff [Eq. (6)] and the combined allocations to the Beijing user (reservoir releases from Ziya River basin, $s{w}_{\text{Beijing},t}$, and SNWTP water from Yangtze River, $s{n}_{\text{Beijing},t}$) are constrained by the capacity of the middle route of the SNWTP, $c{p}_{\text{SNWTP}}$ [Eq. (7)]. The piecewise linear future cost function, $FC$, was added as $h$ linear pieces, as shown in Eq. (8), with slopes ${\lambda }_h$ representing the shadow prices logged in $t+1$, $h$ an index of the discrete reservoir storage levels, and $F{C}_h$ the future costs logged in $t+1$. The linear problem optimizer CPLEXLP (IBM 2013) was chosen for its high efficiency.

The outcome of the optimization [Eq. (2)] is two (stages by states) matrices with minimum total costs (CNY) and shadow prices (${\mathrm{CNY}/\mathrm{m}}^3$), with CNY representing 2005 Chinese yuan. As alternative to directly logging $\lambda$, the water values, $\theta {(\mathrm{CNY}/\mathrm{m}}^3)$, can also be found as the first derivative of the total costs with respect to the discretized reservoir volumes (Stage and Larsson 1961)

To avoid the effect of the end condition ($\mathrm{future}\mathrm{cost}=0$), the model is run to equilibrium by looping the annual input data until the interannual differences in water values become insignificant. The initial year furthest from the end condition is used as the equilibrium water value table. The equilibrium water value tables, one for each inflow class, are used to guide a forward moving simulation phase similar to the optimization phase. From a given reservoir storage at time $t$, the water values are used as future costs to determine the allocation of water between current use and storage for the future. A perfect foresight dynamic program (DP) with a single future cost function (similar setup as the SDP optimizer) is used to evaluate the performance. With perfect foresight, this test will show exactly how well it is possible to manage the water resources; the test can be used to benchmark the performance of the SDP optimization.

A central input for the optimal value function [Eq. (2)] is the curtailment cost of the water users. This cost describes the marginal cost in monetary units (${\mathrm{CNY}/\mathrm{m}}^3$) for not satisfying a user’s water demand. For agricultural users, the calculations were based on literature estimates of the crop water use efficiency, WUE (Table 1), in the lower part of the basin (Mo et al. 2009). The curtailment cost was estimated aswhere $cc$ = curtailment cost of the crop and $p$ = producer’s price of the crop, scaled to 2005 prices by using the consumer price index (World DataBank 2013). The irrigation demands and schedule for the individual crops are based on field interviews from March 2013 of 22 farmers distributed within the Ziya River basin (Table 1) and fall in the same range as in the literature (Liu et al. 2001; Mo et al. 2009). The total irrigated areas were extracted as the land use classes “Dryland Cropland and Pasture” and “Irrigated Cropland and Pasture” in the USGS Euroasia Landcover Map (USGS 2013). These land use classes were evaluated as irrigated based on field observations.

Industrial water demands were estimated from Hai River statistics and scaled with the upstream and downstream areas (Berkoff 2003; Moiwo et al. 2009). The industrial curtailment costs were based on a previous study of the industrial water value in the Hai River basin (World Bank 2001) and not differentiated between upper and lower areas (Table 2).

Domestic demands were estimated from provincial per capita water consumption statistics, as listed in Table 2 [National Bureau of Statistics of China (NBSC) 2011]. The demands were scaled to the upper and lower basins with the 2007 Landscan population density map (Berkoff 2003; Bright et al. 2008). The curtailment costs were based on a previous study of the water scarcity damage costs for the domestic users in the Hai River basin (World Bank 2001). The water demand and curtailment costs of the Beijing user were based on the same study.

The ecosystem water demands could have been included as regular water users in the optimization framework, but owing to a lack of ecosystem water values, ecosystem water requirements were added as demand constraints, based on an estimate of the deficit in the water balance of the Baiyangdian Lake, as listed in Table 2 (People’s Government of Hebei Province 2012).

The supply costs of surface water and SNWTP water were set to zero. Benefits from hydropower production will favor surface water allocations over SNWTP allocations to the downstream users. A constant groundwater pumping price of ${0.4\mathrm{CNY}/\mathrm{m}}^3$ was used, based on field interviews. This price represents only the electricity costs for pumping: pumping costs from irrigation canal to field were not included.

The SNWTP water available for allocation to the Ziya River basin was estimated from the data in Table 2. Different expected water transfer rates for the middle route are reported in the literature, including ${5\mathrm{km}}^3/\mathrm{year}$ (Jia et al. 2012), 9 to ${13\mathrm{km}}^3/\mathrm{year}$ (Berkoff 2003), and ${9.5\mathrm{km}}^3/\mathrm{year}$, as presented in Table 2 (Water-Technology.Net 2013). As a boundary condition for the model, the fraction of water available for allocation has been defined as ${9.5\mathrm{km}}^3/\mathrm{year}$ minus ${7.4\mathrm{km}}^3/\mathrm{year}$ for the 100 major cities in the north (Wang and Ma 1999). The remaining ${2.1\mathrm{km}}^3/\mathrm{year}$ were distributed evenly to the arable land of the NCP and the share available to Ziya was found. The Beijing water deficit of ${1\mathrm{km}}^3/\mathrm{year}$ (Ivanova 2011) was considered part of the ${7.4\mathrm{km}}^3/\mathrm{year}$ value and added to the SNWTP water available in the model.

Monte Carlo simulations are used to assess the uncertainty of economic models. The input uncertainties were initially estimated and a set of samples were generated with Latin hypercube sampling (LHS). The water demands and curtailment costs were assumed to be normally distributed, with SDs of 20% around the estimated value. The hydropower benefits were assumed to be uniformly distributed with 80% uncertainty. Each uncertain input is examined once per sample. With 17 uncertain parameters (eight demands, eight curtailment costs, and hydropower benefits) a sample size of $n=50$ was found to be sufficient. The $n$ results were used to estimate SDs.

The rainfall-runoff model was autocalibrated to measured runoff in the calibration catchment shown in Fig. 2. Despite the absence of major reservoirs in the calibration catchment, the measured runoff included delayed peaks occurring in the dry winter months, as shown for the early years in Fig. 3. The authors expect these peaks to be a result of reservoir releases, because the timing fits the normal irrigation practices in the region. In certain later years, peak summer discharge is very low, despite the occurrence of precipitation events similar to the early period. The highest achievable monthly NSE for all of the 7,700 overlapping days is 0.47 (calibration target), which increases to 0.64 if the winter months are not used for calculation of the NSE. The water balance error, $(\mathrm{sim}-\mathrm{obs})/\mathrm{obs}$, is 10%. The calibration catchment contains multiple smaller reservoirs to serve irrigation agriculture and almost 2 million people. Therefore, the observed discharge may deviate significantly from discharge under natural conditions. Before using these presented results in actual decision making, the modeling framework should be updated with more realistic estimates of the natural water availability, preferably by using observed river discharge time series.

The assumption of stable conditions in the basin is a prerequisite for reaching steady-state water value tables. However, the simulated runoff shows a decreasing tendency over the 51 years. Plotting the accumulated mean precipitation for the three weather stations used for the Shanxi Province (Fig. 4) reveals that the precipitation has been decreasing over the period. This has previously been discussed in the literature (Chen 2010; Sun et al. 2010; Cao et al. 2013). Based on a simple manual fit, the simulation was split into two climate periods with a shift in 1980, in which each period is assumed to have stationary precipitation.

The average flow of each inflow class is shown in Fig. 5. The timing of the peak flow is similar in the period 1980–2008, but the magnitude of the flow is lower. In particular, the dry winter months from November to February are drier in the period 1980–2008.

The backward recursive SDP algorithm [Eq. (2)] was run with a looped 10-year sequence of the annual input data to reach interannual equilibrium water value tables for each climate period. The resulting water value tables for three scenario runs are shown in Fig. 6. The water values are highest at low reservoir states and in the dry periods, with water values above ${0.3\mathrm{CNY}/\mathrm{m}}^3$. In wet months and at high reservoir states, the water value drops below ${0.2\mathrm{CNY}/\mathrm{m}}^3$. Comparing different inflow classes shows that a dry flow class results in higher water values (more conservative reservoir management) than in the normal or wet flow class. The lower reservoir inflow in the period 1980–2008 results in increased water values, and therefore, more conservative reservoir management with lower reservoir releases. The large blue areas with water values around ${0.4\mathrm{CNY}/\mathrm{m}}^3$ in the two scenarios with unlimited groundwater pumping are caused by the groundwater pumping price, which is lower than any of the users’ curtailment costs. In the scenario with the partly finished middle route of the SNWTP, the Beijing user is constrained to surface and SNWTP water alone, because the groundwater in this area is already fully exploited. If no surface water is available, the user can only be curtailed. This increases the water value whenever the reservoir is close to empty, and the demand cannot be satisfied with the in-stage runoff alone. After the SNWTP middle route has been completed, the water diverted from the Yangtze River can completely satisfy the Beijing demand, so the water value at low states becomes lower. Based on the SNWTP data in Table 2, a limit of $109\mathrm{million}{\mathrm{m}}^3/\mathrm{month}$ was placed on the SNWTP water from Yangtze River. In the last scenario, an annual average sustainable groundwater pumping limit for the NCP users is also introduced. A study has modeled the annual NCP groundwater recharge rate as ${17.77\mathrm{km}}^3/\mathrm{year}$ (Liu et al. 2011), which is scaled to the share of the Ziya River basin of the NCP (${3.43\mathrm{km}}^3/\mathrm{year}$) and distributed evenly to obtain a monthly limit. This groundwater pumping limit causes the users with the lowest curtailment costs (the farmers) to be curtailed and increases the water value to ${2–2.5\mathrm{CNY}/\mathrm{m}}^3$.

The expected equilibrium water values were used to drive a 51-year reservoir operation simulation phase. From a starting volume at a given stage with known reservoir state and inflow class, the corresponding water value vector (all possible states in $t+1$) represents benefits of storing water for future use in a forward moving optimization algorithm. The resulting reservoir management and perfect foresight DP solutions are shown for the three scenarios in Fig. 7. In general, the SDP simulations show the same trends as the DP solution, but the apparent interannual cycles primarily after 1980 are not captured by the SDP simulation. The DP solution will save water for these sequences of dry years, whereas the simple Markov chain runoff serial correlation will contain the same probabilities for each year. Therefore, the SDP solution releases more water in the beginning of the dry years and ends up with lower storage than the DP solution. Introducing a groundwater pumping limit (Fig. 7) will make the reservoir management more conservative because higher water values cause more water to be stored for the highest value water uses. Additionally, there is higher storage in the SDP solution than the DP solution. In Figs. 7(a and c), the SDP solution shows higher reservoir states relative to the DP solution. This is caused by nonzero transition probabilities to a low inflow state; therefore, the SDP model saves water rather than curtailing expensive users.

The constraints can be modified to enable the evaluation of a variety of case setups and policy scenarios. Table 3 shows the total costs with SDP and DP for 12 different scenarios. From the difference in total costs between the scenarios, it is possible to calculate the average shadow price of water allocated to ecosystems (ecosystem flow constraint, as indicated in Table 2) or the SNWTP water. The SDs of the results have been calculated from the Monte Carlo simulations. The water diverted from the Yangtze River along the SNWTP middle route will lower the total costs with an average of ${4.6\mathrm{CNY}/\mathrm{m}}^3$ (${\mathrm{SD}=1.3\mathrm{CNY}/\mathrm{m}}^3$). Forcing a minimum ecosystem flow of $100{\mathrm{million}\mathrm{m}}^3$ to fill up the Baiyangdian Lake in July will, on average, cost ${0.41\mathrm{CNY}/\mathrm{m}}^3$ (${\mathrm{SD}=0.13\mathrm{CNY}/\mathrm{m}}^3$) if the current practice continues with unlimited groundwater pumping. This indicates that the water allocated to Baiyangdian Lake will be substituted with groundwater pumping, which is available at a fixed cost of ${0.4\mathrm{CNY}/\mathrm{m}}^3$. If an average groundwater pumping limit is introduced, the average shadow price of the ecosystem water will be ${2.78\mathrm{CNY}/\mathrm{m}}^3$ (${\mathrm{SD}=1.04\mathrm{CNY}/\mathrm{m}}^3$), indicating that the diversion will cause curtailment of the farmers. The differences between the objective values of SDP and perfect foresight DP are between 0.3 and 4.7% (average 1.8%) for scenarios without a groundwater limit and between 3.6 and 5.7% (average 4.3%) for scenarios with a groundwater limit.

The impact of describing the runoff serial correlation can be found by comparing with a model run based on average monthly runoff from the 51 years. For the three scenarios presented in Figs. 6 and 7, the total objective values with runoff serial correlation were 3.70 billion $\mathrm{CNY}/\mathrm{year}$ ($\mathrm{SD}=0.6\mathrm{billion}\mathrm{CNY}/\mathrm{year}$), $3.09\mathrm{billion}\mathrm{CNY}/\mathrm{year}$ ($\mathrm{SD}=0.5\mathrm{billion}\mathrm{CNY}/\mathrm{year}$), and $11.39\mathrm{billion}\mathrm{CNY}/\mathrm{year}$ ($\mathrm{SD}=3.7\mathrm{billion}\mathrm{CNY}/\mathrm{year}$), as listed in Table 3. Using average monthly flows instead, these total costs become 3.70, 3.21, and $12.02\mathrm{billion}\mathrm{CNY}/\mathrm{year}$. The total costs with perfect foresight of 3.53, 3.06, and $10.75\mathrm{billion}\mathrm{CNY}/\mathrm{year}$ for these scenarios indicate the direct gain of using the stochastic representation of the runoff. The tradeoff for the higher accuracy of the SDP model is a computation time three times longer (66 s) than the model based on average flows (23 s).