The cost of ending groundwater overdraft on the North China Plain

. Overexploitation of groundwater reserves is a major environmental problem around the world. In many river basins, groundwater and surface water are used conjunctively and joint optimization strategies are required. A hydroeconomic modeling approach is used to ﬁnd cost-optimal sustainable surface water and groundwater allocation strategies for a river basin, given an arbitrary initial groundwater level in the aquifer. A simpliﬁed management problem with conjunctive use of scarce surface water and groundwater under inﬂow and recharge uncertainty is presented. Because of head-dependent groundwater pumping costs the optimization problem is nonlinear and non-convex, and a genetic algorithm is used to solve the one-step-ahead subproblems with the objective of minimizing the sum of immediate and expected future costs. A real-world application in the water-scarce Ziya River basin in northern China is used to demonstrate the model capabilities. Persistent overdraft from the groundwater aquifers on the North China Plain has caused declining groundwater levels. The model maps the marginal cost of water in different scenarios, and the minimum cost of ending groundwater overdraft in the basin is estimated to be CNY 5.58 billion yr − 1 . The study shows that it is cost-effective to slowly recover the groundwater aquifer to a level close to the surface, while gradually lowering the groundwater value to the equilibrium at CNY 2.15 m − 3 . The model can be used to guide decision-makers to economic efﬁcient long-term sustainable management of groundwater and surface water resources.

tation of the groundwater system using spatially distributed groundwater models (Andreu et al., 1996;Harou and Lund, 2008;Marques et al., 2006;Pulido-Velázquez et al., 2006). Stochasticity is commonly represented in scenarios where a regression analysis is used to formulate operation rules; see e.g., the implicit stochastic optimization approaches reviewed by Labadie (2004). Singh (2014) reviewed the use of simulation-optimization (SO) modeling for conjunctive groundwater and surface water use. In SO-based studies, efficient groundwater simulation models are used to answer "what if" questions, while an optimization model is wrapped around the simulation model to find "what is best". Groundwater aquifers have been represented as simple deterministic box or "bathtub" models (e.g., Cai et al., 2001;Riegels et al., 2013) and as spatially distributed models (e.g., Maddock, 1972;Siegfried et al., 2009) with stochasticity (Reichard, 1995;Siegfried and Kinzelbach, 2006). While the results obtained from these methods are rich in detail, they yield only a single solution to the optimization problem.
Methods based on dynamic programming (DP, Bellman, 1957) have been used extensively to demonstrate the dynamics of conjunctive groundwater-surface water use for both deterministic (e.g., Buras, 1963;Provencher and Burt, 1994;Yang et al., 2008) and stochastic DP (SDP, e.g., Burt, 1964;Philbrick and Kitanidis, 1998;Provencher and Burt, 1994;Tsur and Graham-Tomasi, 1991) optimization problems. In DP-based methods, the original optimization problem is decomposed into subproblems, which are solved sequentially over time. The entire decision space is thereby mapped, enabling use of the results as dynamic decision rules. However, the number of subproblems grows exponentially with the number of state variables and this curse of dimensionality has frequently limited the use of DP and SDP (Labadie, 2004;Provencher and Burt, 1994;Saad and Turgeon, 1988). Although it causes loss of detail and inability to disaggregate the results, reservoir aggregation has been suggested as a possible solution strategy (Saad and Turgeon, 1988).
This study aims to answer the following two macroscale decision support questions for conjunctive groundwater and surface water management for the Ziya River basin in North China. (1) What are the minimum costs of ending groundwater overdraft? (2) What is the cost-efficient recovery strategy of the overpumped aquifer? A hydroeconomic modeling approach is used to identify the least-cost strategy to achieve sustainable groundwater abstraction, defined as the long-term average abstraction that does not exceed the longterm average recharge. To overcome the management problem similar to Harou and Lund (2008) with increased complexity caused by uncertain surface water runoff and groundwater recharge, the surface water reservoirs are aggregated. This is adequate at a macroscale (Davidsen et al., 2015) and allows the use of DP-based approaches. The cost minimization problem is solved with the water value method, a variant of SDP (Stage and Larsson, 1961;Stedinger et al., 1984) which produces dynamic tables of marginal costs linked to states, stages, and water source. Head-and rate-dependent pumping costs introduce nonlinearity in the discrete subproblems. This nonlinearity is handled with a hybrid genetic algorithm (GA) and linear programming (LP) method similar to that used by Cai et al. (2001), here applied in a coupled groundwater-surface water management problem within an SDP framework.

Study area
Northern China and particularly the North China Plain (NCP) have experienced increasing water scarcity problems over the past 50 years due to population growth, economic development, and reduced precipitation (Liu and Xia, 2004). The deficit in the water balance has historically been covered by overexploitation of the groundwater aquifer, causing a regional lowering of the groundwater table by up to 1 m yr −1 (Zheng et al., 2010).
The Ziya River basin, a part of the Hai River basin, was selected as a case study area (see Fig. 1). The upper basin is located in the Shanxi province, while the lower basin is located in the Hebei province on the NCP. The 52 300 km 2 basin has approximately 25 million inhabitants (data from 2007, Bright et al., 2008), and severe water scarcity is causing multiple conflicts. Five major reservoirs with a combined storage capacity of 3.5 km 3 are located in the basin. While reservoir rule curves and flood control volumes can easily be accommodated for in the optimization approach, current reservoir management rule curves in the case area are unavailable to the public. Instead it is assumed that the full storage capacity can be managed flexibly without consideration of storage reserved for flood protection or existing management rules. Incorporating flood storage volumes will reduce the available storage and increase water scarcity in the long dry season. In the present model setup, we therefore find the lower limit on water scarcity costs, assuming that the entire storage capacity is available for storing water. Reservoir spills will cause an economic loss, and the model tends to avoid spills by entering the rainy season with a low reservoir storage level.
A previous hydroeconomic study of the Ziya River basin was a traditional implementation of SDP on a singlereservoir system (surface water reservoir) and showed optimal water management, while disregarding dynamic groundwater storage and head-dependent groundwater pumping costs (Davidsen et al., 2015). Instead, the groundwater resource was included as a simple monthly upper allocation constraint.
In the present study, the groundwater resource is included as a dynamic aquifer box model with a storage capacity of 275 km 3 . The river basin has two aquifers (upstream and downstream) which are only connected by the river. Ideally, each aquifer should be modeled as a box model, but this extra state variable would be computationally challenging within the SDP framework. We therefore set up a box model for the downstream and most important aquifer only, and abstraction from the upstream aquifer is only bounded by an upper pumping limit corresponding to the average monthly recharge. The box model for the downstream aquifer is formulated as Infiltration + Storage = Pumping + Overflow. The groundwater overflow is only used in extreme cases, where the total pumping and available storage is less than the infiltration. The spills will go to the spill variable and leave the system as baseflow to the rivers (unavailable for allocation). The aquifer is so heavily overexploited that no significant baseflow is being created or will be created in the foreseeable future. The box model allows for more flexible management with large abstractions in dry years and increased recharge in wet years. The groundwater aquifer can thereby be used to bridge longer drought periods. Except for the groundwater box model, the conceptual model is identical to the one used by Davidsen et al. (2015). A conceptual sketch of the management problem is shown in Fig. 2. The water users are divided into groups of economic activities; irrigation agriculture, industrial users, and domestic water users. Ideally, each water user group should be characterized by flexible demand curves, but due to poor data availability a constant water demand (m 3 ) and a constant curtailment cost of not meeting the demand were used for each group (see Table 1). The water demands are assumed to be deterministic and decoupled from the stochastic runoff. This is a reasonable assumption because the rainfall on the NCP normally occurs in the summer months, while irrigation water demands are concentrated in the dry spring. The irrigation schedule is centrally planned and typically unchanged from year to year. The upstream (u) users have access to runoff and are restricted to an upper pumping limit X gw cor- Figure 2. Conceptual sketch of the simplified water management problem. Water users are located upstream (u) and downstream (d) of a surface water reservoir. Allocation decision variables for surface water (blue), SNWTP water (green), and groundwater (orange) are indicated. The conceptual sketch of the downstream dynamic aquifer shows how the total lift ( h) is composed of the thickness of the top layer, the regional groundwater drawdown, and the local Thiem steady-state groundwater drawdown.
responding to the average monthly upstream recharge, while the downstream users (d) have access to reservoir releases, water delivered through the South-to-North Water Transfer Project (SNWTP) and groundwater from the downstream aquifer.

Optimization model formulation
An SDP formulation is used to find the expected value of storing an incremental amount of surface water or groundwater, given the month of the year, the available storage in surface and groundwater reservoirs, and the inflow scenarios. The backward recursive equation calculates the sum of immediate and expected future costs for all combinations of discrete reservoir storage levels (states) and monthly time steps (stages). The immediate management cost (IC) arises from water supply and water curtailment, whereas the expected future cost (EFC) is the optimal value function in t +1 weighed by the corresponding transition probabilities. In the present setup, we decided to weigh the IC and EFC equally, but inclusion of discount rates other than zero is possible. Because of the head-and rate-dependent groundwater pumping costs, which will be described in detail later, the immediate cost depends nonlinearly on the decision variables. The objective is to minimize the total costs over the planning period, given by the optimal value function F * t V gw, t , V sw, t , Q k sw, t based on the classical Bellman formulation: with IC being subject to x sw, m, t + x gw, m, t + x SNWTP, m, t + x ct, m, t = dm m, t U u=1 x gw, u, t ≤ X gw, t See Table 2 for nomenclature. Equation (3) is the water demand fulfillment constraint; i.e., the sum of water allocation and water curtailments equals the water demand of each user. Equation (4) is the water balance of the combined surface water reservoir, while Eq. (5) is the water balance of the reservoir releases. A similar water balance for the dynamic groundwater aquifer follows in Eq. (6). The upstream surface water allocations are constrained by the upstream runoff as shown in Eq. (7), while the upstream groundwater allocations are constrained to a fixed sustainable monthly average as shown in Eq. (8). In Eq. (9), the upper and lower hard constraints on the decision variables are shown. Lastly, Eq. (10) is the marginal groundwater pumping cost, which depends on the combined downstream groundwater allocations as described later.
A rainfall-runoff model based on the Budyko Framework (Budyko, 1958;Zhang et al., 2008) has been used in a previous study to estimate the near-natural daily surface water runoff into reservoirs (Davidsen et al., 2015). The resulting 51 years  of simulated daily runoff was aggregated to monthly runoff and normalized. A Markov chain, which describes the runoff serial correlation between three flow classes defined as dry (0-20th percentile), normal (20-80th percentile), and wet (80-100th percentile), was established and validated to ensure second-order stationarity (Davidsen et al., 2015;Loucks and van Beek, 2005). The groundwater recharge is estimated from the precipitation data also used in the rainfall-runoff model. The average monthly precipitation (mm month −1 ) for each runoff class is calculated, and a simple groundwater recharge coefficient of 17.5 % of the precipitation (Wang et al., 2008) is used.
The SDP loop is initiated with EFC set to zero and will propagate backward in time through all the discrete system states as described in the objective function. For each discrete combination of states, a cost minimization subproblem will be solved. A subproblem will have the discrete reservoir storage levels (V gw, t and V sw, t ) as initial conditions and reservoir inflow given by the present inflow class in the Markov chain. The optimization algorithm will search for the optimal solution, given the costs of the immediate management (water allocations and water curtailments, including reservoir releases and groundwater pumping) which have to be balanced against the expected future costs. As the SDP algorithm propagates backward in time, the future costs will be equal to the minimum total costs from t+1, weighted by maximum monthly groundwater pumping in the upstream basin (m 3 month −1 ) q E, t unused surface water available to ecosystems, decision variable (m 3 month −1 ) Q E minimum in-stream ecosystem flow constraint (m 3 month −1 ) Bei Beijing user Q SNWTP maximum capacity of the SNWTP canal (m 3 month −1 ) the Markov chain transition probabilities. The algorithm will continue backward in time until equilibrium is reached, i.e., until the shadow prices (marginal value of storing water for future use) in two successive years remain constant. The SDP model is developed in MATLAB (MathWorks Inc., 2013) and uses the fast cplexlp (IBM, 2013) to solve the linear subproblems.
The sets of equilibrium shadow prices, referred to as the water value tables, can subsequently be used to guide optimal water resources management forward in time with unknown future runoff. In this study, the available historic runoff time series is used to demonstrate how the derived water value tables should be used in real-time operation. The simulation will be initiated from different initial groundwater aquifer storage levels, thereby demonstrating which pricing policy should be used to bring the NCP back into a sustainable state.

Dynamic groundwater aquifer
The groundwater aquifer is represented as a simple box model (see Fig. 2) with recharge and groundwater pumping determining the change in the stored volume of the aquifer (Eq. 6). The pumping is associated with a pumping cost determined by the energy needed to lift the water from the groundwater table to the land surface (Eq. 10): where P is the specific pump energy (J m −3 ), ρ is the density of water (kg m −3 ), g is the gravitational acceleration (m s −2 ), h is the head difference between groundwater table and land surface (m), and ε is the pump efficiency (−). The marginal pumping cost c gw (CNY m −3 ) is found from the average electricity price c el (CNY/Ws) in northern China: Hence this cost will vary with the stored volume in the groundwater aquifer. The present electricity price structure in China is quite complex, with the users typically paying between CNY 0.4 and 1 kWh −1 depending on power source, province, and consumer type (Li, 2012;Yu, 2011). In this study a fixed electricity price of CNY 1 kWh −1 is used. The immediate costs of supplying groundwater to a single user follow where h is found as the mean depth from the land surface to the groundwater table (see Fig. 2) between t and t + 1: where h top is the distance from the land surface to the top of the aquifer at full storage (m), S y is the specific yield (−) of the aquifer, and A is the area of the aquifer (m 2 ).
Here V gw, t+1 is a decision variable, and once substituted into Eq. (13) it is clear that the problem becomes nonlinear. In Eq. (14) the drawdown is assumed uniform over the entire aquifer. This simplification might be problematic as the local cone of depression around each well could contribute significantly to the pumping cost and thereby the optimal policy. Therefore, the steady-state Thiem drawdown (Thiem, 1906) solution is used to estimate local drawdown at the pumping wells. Local drawdown is then added to Eq. (15) to estimate total required lift: where Q w is the pumping rate of each well (m 3 month −1 ), T is the transmissivity (m 2 month −1 ), r in is the radius of influence (m), and r w is the distance from origin to the point of interest (m), which here is the radius of the well. The transmissivity is based on a hydraulic conductivity of 1.3 × 10 −6 m s −1 for silty loam (Qin et al., 2013). The hydraulic conductivity is lower than the expected average for the NCP to provide a conservative estimate of the effect of drawdown. Field interviews revealed that the wells typically reach no deeper than 200 m below surface, which results in a specific yield of 5 %. The groundwater pumping Q w is defined as the total allocated groundwater within the stage (m 3 month −1 ) and assumed evenly distributed among the number of wells in the catchment: where n w is the number of wells in the downstream basin. The even pumping distribution is a fair assumption, as field investigations showed that (1) the majority of the groundwater wells are for irrigation; (2) the timing of irrigation, crop types, and climate is homogeneous; and (3) the groundwater wells have comparable capacities. Erlendsson (2014) estimated the well density in the Ziya River basin from Google Earth to be 16 wells km −2 . Assuming that the wells are distributed evenly on a regular grid and that the radius of influence r in is 500 m, overlapping cones of depression from eight surrounding wells are included in the calculation of the local drawdown. This additional drawdown is included using the principle of superposition as also applied by Erlendsson (2014).

Solving nonlinear and non-convex subproblems
With two reservoir state variables and a climate state variable, the number of discrete states is quickly limited by the curse of dimensionality. A very fine discretization of the groundwater aquifer to allow discrete storage levels and decisions is computationally infeasible. A low number of discrete states increases the discretization error, particularly if both the initial and the end storages V gw, t+1 and V sw, t+1 are kept discrete. The discrete volumes of the large aquifer become much larger than the combined monthly demands, and storing all recharge will therefore not be sufficient to recharge to a higher discrete storage level. Similarly, the demands will be smaller than the discrete volumes, and pumping the remaining water to reach a lower discrete level would also be infeasible. Allowing free end storage in each subproblem will allow the model to pick, e.g., the optimal groundwater recharge and pumping without a requirement of meeting an exact discrete end state. With free surface water and groundwater end storages, the future cost function has three dimensions (surface water storage, groundwater storage, and expected future costs). Pereira and Pinto (1991) used Benders' decomposition approach, which employs piecewise linear approximations and requires convexity. With head-and rate-dependent pumping costs and increasing electricity price, we observed that the future cost function changes from strictly convex (very low electricity price) to strictly concave (very high electricity price). At realistic electricity prices, we observed a mix of concave and convex shapes. An alternative is to use linear interpolation with defined upper and lower bounds. However, with two state variables, interpolation between the future cost points will yield a hyperplane in three dimensions, which complicates the establishment of boundary conditions for each plane. Nonlinear optimization problems can be solved with evolutionary search methods, a subdivision of global optimizers. A widely used group of evolutionary search methods are genetic algorithms (GAs), which are found to be efficient tools for getting approximate solutions to complex nonlinear optimization problems (see, e.g., Goldberg, 1989;Reeves, 1997). GAs use a random search approach inspired by natural evolution and have been applied to the field of water resources management by, e.g., Cai et al. (2001), McKinney and Lin (1994), and Nicklow et al. (2010). Cai et al. (2001) used a combined GA and LP approach to solve a highly nonlinear surface water management problem. By fixing some of the complicating decision variables, the remaining objective function became linear and thereby solvable with LP. The GA was used to test combinations of the fixed parameters, while looking for the optimal solution. The combination yielded faster computation time than if the GA was used to estimate all the parameters.
A GA implemented in MATLAB is used to solve the cost minimization subproblems. This GA function will initially generate a set of candidate solutions known as the popula-Hydrol. Earth Syst. Sci., 20, 771-785, 2016 www.hydrol-earth-syst-sci.net/20/771/2016/ tion. Each of the candidate solutions contains a set of decision variables (sampled within the decision space), which will yield a feasible solution to the optimization problem. In MATLAB, a set of options specifies the population size, the stopping criteria (fitness limit, stall limit, function tolerance, and others), the crossover fraction, the elite count (number of top parents to be guaranteed survival), and the generation function (how the initial population is generated). The options were adjusted to achieve maximum efficiency of the GA for the present optimization problem. The computation time for one single subproblem is orders of magnitude larger than solving a simple LP. As the optimization problem became computationally heavier with increasing number of decision variables, a hybrid version of GA and LP, similar to the method used by Cai et al. (2001), was developed (see Fig. 3). Decision variables that cause nonlinearity are identified and chosen by the GA. Once these complicating decision variables are chosen, the remaining objective function becomes linear and thereby solvable with LP. In the optimization problem presented in Eq. (1), the nonlinearity is caused by the head-dependent pumping costs as explained in Eq. (13)-(14). Both the regional lowering of the groundwater table and the Thiem local drawdown cones depend on the decision variable for the stored volume in t + 1,V gw, t+1 . If V gw, t+1 is preselected, the regional drawdown is given, and the resulting groundwater pumping rate Q w can be calculated from the water balance. The groundwater pumping price is thereby also given, and the remaining optimization problem becomes linear.
For a given combination of stages, discrete states and flow classes, the objective of the GA is to minimize the total cost, TC, with the free states V gw, t+1 , V sw, t+1 being the decisions: with EFC being the expected future cost. Given initial states and once the GA has chosen the end states, the immediate cost minimization problem becomes linear and hence solvable with LP (see Fig. 3). The expected future costs are found by cubic interpolation of the discrete neighboring future cost grid points in each dimension of the matrix. The GA approaches the global optimum until fitness limit criteria are met. The total costs are stored and the algorithm continues to the next state. To reduce the computation time, the outer loop through the groundwater states is parallelized.
The performance of the GA-SDP model is compared to a fully deterministic DP, which finds the optimal solution given perfect knowledge about future inflows and groundwater recharge. The DP model uses the same algorithm as the SDP model and one-dimensional state transition matrices with p = 1 between the deterministic monthly runoff data. For low storage capacity and long timescales, the effect of the end storage volume becomes negligible.  model, the DP model was looped until the present management was no longer affected by the end of period condition.

Results
Without any regulation or consideration of the expected future costs arising from overexploitation of the groundwater aquifer, the water users will continue maximizing immediate profits (producers) or utility (consumers). Because there are only electricity costs for groundwater, the users will continue pumping groundwater until the marginal groundwater cost exceeds the curtailment cost. At CNY 1 kWh −1 the marginal cost of lifting groundwater 200 m (typical depth of wells observed in the study area) can be found with Eqs. (13)- (14) to be CNY 0.8 m −3 and thereby less than the lowest curtailment cost at CNY 2.3 m −3 . Only once the electricity price is higher than CNY 2.8 kWh −1 does the user with the lowest curtailment cost stop pumping, if the groundwater level is 200 m below the surface. The backward recursive SDP algorithm was run with a looped annual data set until equilibrium water values, i.e., no interannual changes, were obtained. The water values increase fastest during the first years, and after approximately 100 years, the annual increases become small. Due to the large storage capacity of the groundwater aquifer, equilibrium is however not achieved until after 150-180 years. sample of the resulting equilibrium water value tables is presented in Fig. 4. This figure shows the temporal variations of water values as a function of one state variable, keeping the other state variables at a fixed value. The state variables are fixed at empty, half full, and full storage respectively. During the rainy season from June to August, high precipitation rates reduce water scarcity, resulting in lower surface water values. Because the groundwater storage capacity is much larger, increased recharge can easily be stored for later use, and groundwater values are therefore not affected. Addition of stream-aquifer interactions to the model is expected to affect this behavior, but since the flow in rivers/canals in the case study area is small most of the year, and since most areas are far from a river, it is reasonable to ignore these dynamics. The water values after 1980 are clearly higher than in the period before 1980 due to increased water scarcity caused by a reduction in the regional precipitation. In contrast, the groundwater value tables are uniform, with variation only in groundwater storage. The detailed water value tables are included in the Supplement. We simulate management using the equilibrium water value tables as a pricing policy and force the system with 51 years of simulated historical runoff. Time series of the simulated groundwater storage can be seen in Fig. 5 for different initial storage scenarios. The groundwater aquifer approaches an equilibrium storage level around 260 km 3 (95 % full). If the storage in the aquifer is below this level, the average recharge will exceed average pumping until the equilibrium storage is reached. If the storage level is above equilibrium, average pumping will exceed average recharge until equilibrium is reached. In Fig. 6, the surface water and Results for the business-as-usual scenario (EFC is disregarded) and a benchmark scenario with perfect foresight (DP) are also shown.
groundwater storages are shown for a situation with equilibrium groundwater storage. In most years, the surface water storage falls below 1 km 3 , leaving space in the reservoir for the rainy season. The potential high scarcity cost of facing a dry scenario with an almost empty reservoir is avoided by pumping more groundwater. These additional pumping costs seem to be exceeded by the benefits of minimizing spills in the rainy season. To demonstrate the business-as-usual solution, the simulation model is run for a 20-year period with the present water demands and curtailment costs and with a discount rate set to infinity (i.e., zero future costs). The resulting groundwater table is continuously decreasing, as shown in Fig. 5. In the simulated management runs, water will be allocated to the users up to a point where reductions in immediate cost are compensated by increases in expected future costs. The User's price for groundwater and surface water through for a 51-year simulation based on simulated historical runoff for two initial groundwater storages. P refers to the user's price. M1 denotes results for a single surface water reservoir with a constant groundwater cost (Davidsen et al., 2015). M2 denotes results from the presented model framework with an additional dynamic groundwater aquifer. The user's price for groundwater in M2 is the immediate pumping cost added to the marginal cost from the water value tables.
user's price, which can be applied in a marginal cost pricing (MCP) scheme, is the marginal value of the last unit of water allocated to the users. The user's price is the sum of the actual pumping cost (electricity used) and the additional marginal cost given by the equilibrium water value tables. In Fig. 7, the user's prices for groundwater and surface water are shown for the 51-year simulation at and below the long-term sustainable groundwater storage level. When the groundwater storage level is close to equilibrium, the user's prices of groundwater and surface water are equal during periods with water scarcity. In wet months with reduced water scarcity, the model switches to surface water allocation only, and the groundwater user's price is undefined (gaps in the time series in Fig. 7). If the groundwater storage level is below equilibrium, the groundwater user's price will be higher, causing an increase in water curtailments and an increase in storage level, as shown in Fig. 5. Under these circumstances the surface water user's price increases up to a point where the two prices meet. With an initial aquifer storage at one-third of the aquifer capacity (100 km 3 ), the groundwater value is The results are shown for a simple drawdown model with uniform regional lowering of the groundwater table, and a more realistic drawdown model, which includes the stationary Thiem local drawdown cones.
CNY 3 m −3 (see Fig. 7). As the aquifer slowly recovers, the groundwater price decreases gradually. At the equilibrium groundwater storage level, the user's prices for groundwater is stable at around CNY 2.15 m −3 , as shown in Fig. 7. This indicates curtailment of wheat agriculture in the downstream Hebei province, which has a willingness to pay CNY 2.12 m −3 (see Table 1). The allocation pattern to this user is shown in Fig. 8; the model switches between high curtailment and high allocations, depending on water availability and storage in the reservoirs. Groundwater allocations fluctuate between satisfying 0 and 80 % of the demand. Inclusion of the steady-state Thiem drawdown cones in the optimization model increases the marginal groundwater pumping cost with increased pumping rates. Groundwater allocations are distributed more evenly over the months, which results in less local drawdown. The total curtailments remain constant, while 1 % of the total water abstraction is shifted from groundwater to surface water, if the stationary Thiem drawdown is included. Inclusion of well drawdown significantly changed the simulated management but resulted in an only slightly increased computation time.
The average total costs of the 51-year simulation for different scenarios can be seen in Table 3. The average reduction in the total costs, associated with the introduction of the SNWTP canal, can be used to estimate the expected marginal economic impact of the SNWTP water. The minimum total costs after the SNWTP is put in operation are compared to the scenario without the SNWTP (pre-2008) and divided by the allocated SNWTP water. The resulting marginal value of the SNWTP water delivered from Shijiazhuang to Beijing (2008-2014 scenario) is CNY 3.2 m −3 , while the SNWTP water from Yangtze River (post-2014 scenario) reduces the total costs with CNY 4.9 m −3 . Similarly, a comparison of the total costs for the post-2014 scenarios shows a marginal increase of CNY 0.91 m −3 as a consequence of introducing a minimum in-stream flow constraint. Table 3. Average minimum total cost (TC) and hydropower benefits (HP) over the 51-year planning period for different scenario runs. SNWTP scenarios: pre-2008 refers to the period before the central route of the SNWTP was built; 2008-2014 refers to a situation with a partly completed SNWTP central route from Shijiazhuang to Beijing; post-2014 refers to the fully operational SNWTP central route from the Yangtze River to Beijing. Scenarios: LGW is initial groundwater storage at 100 km 3 (all other scenarios are initiated at equilibrium groundwater storage); dm is 10 % higher water demands; ct is 10 % higher curtailment costs ; T is 10 % higher transmissivity; TD is Thiem steady-state drawdown; E is minimum ecosystem flow constraint; "+" is active and "−" is inactive. A local sensitivity analysis focused on the water demands and curtailment costs used directly in the objective function (Eq. 1) and the transmissivity used to estimate the local drawdown (Eq. 14). The uncertain input parameters were increased by 10 %, and the sensitivity was evaluated based on the simulation results. The resulting total costs can be seen in Table 3. The benchmark DP was run for the post-2014 scenario with Thiem drawdown and minimum ecosystem flow constraint. The minimum total cost of this run is Hydrol. Earth Syst. Sci., 20, 771-785, 2016 www.hydrol-earth-syst-sci.net/20/771/2016/ CNY 8.46 × 10 9 yr −1 . This is 1.3 % lower than the equivalent SDP run (CNY 8.56 × 10 9 yr −1 ). The minimum total costs were lowered from CNY 10.50 × 10 9 yr −1 (Davidsen et al., 2015) to CNY 8.56 × 10 9 yr −1 (18 % reduction) by allowing the groundwater aquifer to be utilized as buffer instead of a fixed monthly volume. This difference highlights the problem of defining realistic boundaries of optimization problems and shows that simple hard constraints, here fixed groundwater pumping limits, can highly limit the optimal decision space. With a dynamic groundwater aquifer, the model can mitigate dry periods and stabilize the user's price of surface water as shown in Fig. 7. Finally, policies like minimum in-stream ecosystem flow constraints can be satisfied with less impact on users with high curtailment costs. The total costs without restrictions on the groundwater pumping have been estimated to be CNY 2.98 × 10 9 ,yr −1 (Davidsen et al., 2015). To end the groundwater overdraft in the basin, the present study thus estimates a cost increase of CNY 5.58 × 10 9 yr −1 , once the groundwater aquifer is at equilibrium storage. The cost of recharging the aquifer from the present storage level below the equilibrium is significantly higher. In Table 3, the LGW scenario shows that the average cost of sustainable management from an initial storage at 100 km 3 (one-third full) is CNY 13.32 × 10 9 yr −1 .

SNWTP
From any initial groundwater reservoir storage level, the model brings the groundwater table to an equilibrium storage level at approximately 95 % of the aquifer storage capacity. Only small variations in the aquifer storage level are observed after the storage level reaches equilibrium as shown in Fig. 6. While addition of the Thiem stationary drawdown has only a small effect on total costs and total allocated water, it is clear from Fig. 8 that the additional Thiem drawdown highly impacts the allocation pattern for some of the water users. High groundwater pumping rates result in larger local drawdown and thus in higher pumping costs. This mechanism leads to a more uniform groundwater pumping strategy, which is clearly seen in Fig. 8, and results in a much more realistic management policy.

Discussion
This study presents a hydroeconomic optimization approach that provides a macroscale economic pricing policy in terms of water values for conjunctive surface water-groundwater management. The method was used to demonstrate how the water resources in the Ziya River basin should be priced over time, to reach a sustainable situation at minimum cost. We believe that the presented modeling framework has great potential use as a robust decision support tool in real-time water management. However, a number of limitations and simplifications need to be discussed.
A first limitation of the approach is the high level of simplification needed. There are two main reasons for the high level of simplification: limited data availability and the lim-itations of the SDP method. The curse of dimensionality limits the approach to 2-3 interlinked storage facilities and higher dimensional management problems will not be computationally feasible with SDP today. This limit on the number of surface water reservoirs and groundwater aquifers requires a strongly simplified representation of the real-world situation in the optimization model. The simulation phase following the optimization is not limited to the same extent, since only a single subproblem is solved at each stage. The water values determined by the SDP scheme can thus be used to simulate management using a much more spatially resolved model with a high number of users; this was not demonstrated in this study. The advantage of SDP is that it provides a complete set of pricing policies that can be applied in adaptive management, provided that the system can be simplified to a computationally feasible level. An alternative approach known as stochastic dual dynamic programming (SDDP, Pereira and Pinto, 1991;Pereira et al., 1998) has shown great potential for multi-reservoir river basin water management problems. Instead of sampling the entire decision space with the same accuracy level, SDDP samples with a variable accuracy not predefined in a grid, focusing the highest accuracy around the optimal solution. This variable accuracy makes SDDP less suitable for adaptive management. Despite the highly simplified system representation, we believe that the modeling framework provides interesting and nontrivial insights, which are extremely valuable for water resources management on the NCP.
Computation time was a limitation in this study. Three factors increased the computational load of the optimization model. (1) Inclusion of the groundwater state variable resulted in an exponential growth of the number of subproblems; (2) the non-convexity handled by the slower GA-LP formulation caused an increase in the computation time of 10-100 times a single LP; and (3) the SDP algorithm needed to run through more than 200 years to reach a steady state. A single scenario run required 4000 CPU hours and was solved in 2 weeks, using 12 cores at the high-performance computing facilities at the Technical University of Denmark. This is 50 000 times more CPU hours than a single-reservoir SDP model (Davidsen et al., 2015). Since the water value tables can be used offline in the decision making, this long computation time can be accepted.
The long computation time made the use of, e.g., Monte Carlo-based uncertainty analysis infeasible. The local sensitivity analysis showed that a 10 % increase in the curtailment costs is returned as a 6.0 % increase in the total costs, while a similar increase of the demands generates a 2.1 % increase in costs. The transmissivity can vary over many orders of magnitude because it is a log-normally distributed variable. The sensitivity of log(T ) is high: a 1.3 % change of log(T ) from the baseline value results in a 1.5 % change in the cost. At the same time, the simple system representation needed in SDP required assumptions of inflow and storage discretization, aggregation of the surface water reservoirs, generalized estimates of pumping cost, and a lumped groundwater model, which all contribute to the uncertainty. Further, poor data availability for the case study area required some rough estimates of the natural water availability, single-point demand curves, and perfect correlation between rainfall and groundwater recharge. The method-driven assumptions generally limit the decision support to the basin scale, while the simple estimates caused by poor data availability contribute to raising the general uncertainty of the model results. Given the computational challenges and the diverse and significant uncertainties, the model results should be seen as a demonstration of the model's capabilities rather than precise cost estimates. Better estimates will require access to a more comprehensive case data set and involve a complete sensitivity analysis.
Intuitively, one would expect the equilibrium groundwater storage level to be as close as possible to full capacity, while still ensuring that any incoming groundwater recharge can be stored. Finding the exact equilibrium groundwater storage level would require a very fine storage discretization, which, given the size of the groundwater storage, is computationally infeasible. Therefore the equilibrium groundwater storage level is subject to significant discretization errors. The long time steps (monthly) make the stationarity required for using the Thiem stationary drawdown method a realistic assumption.
The difference between total cost with SDP and with DP (perfect foresight) is small (1.3 %). Apart from Beijing, which has access to the SNWTP water, the remaining downstream users have unlimited access to groundwater. The large downstream groundwater aquifer serves as a buffer to the system and eliminates the economic consequences of a wrong decision. The model almost empties the reservoir every year as shown in Fig. 6, and wrong decisions are not punished with curtailment of expensive users as observed by Davidsen et al. (2015). The groundwater aquifer reduces the effect of wrong decisions by allowing the model to minimize spills from the reservoir without significant economic impact of facing a dry period with an empty reservoir. A dynamic groundwater aquifer thereby makes the decision support more robust, since it is the timing and not the amount of curtailment that is being affected.
The derived equilibrium groundwater value tables in Fig. 4 (and detailed water value tables in the Supplement) show that the groundwater values vary with groundwater storage alone and are independent of time of the year, the inflow and recharge scenario, and the storage in the surface water reservoir. This finding is important for future work, as a substitution of the groundwater values with a simpler cost function could greatly reduce the number of states and thereby the computation time. The equilibrium groundwater price, i.e., the groundwater values around the long-term equilibrium groundwater storage, can possibly be estimated from the total renewable water and the water demands ahead of the optimization, but further work is required to test this.
Further work should also address the effect of discounting of the future costs on the equilibrium water value tables and the long-term steady-state groundwater table. In the present model setup, the large groundwater aquifer storage capacity forces the backward-moving SDP algorithm to run through 200-250 model years, until the water values converge to the long-term equilibrium. Another great improvement, given the availability of the required data, would be to replace the constant water demands with elastic demand curves in the highly flexible GA-LP setup.
A significant impact of including groundwater as a dynamic aquifer is the more stable user's prices shown in Fig. 7. The user's price of groundwater consists of two parts: the immediate groundwater pumping costs (electricity costs) and the expected future costs represented by the groundwater value for the last allocated unit of water. As the model is run to equilibrium, the user's prices converge towards the long-term equilibrium at approximately CNY 2.2 m −3 . The electricity price can be used as a policy tool to internalize the user's prices of groundwater shown in Fig. 7. Stable water user's prices will ease the implementation of, e.g., an MCP scheme, which is one of the available policy options to enforce long-term sustainability of groundwater management.

Conclusions
This study describes development and application of a hydroeconomic approach to optimally manage conjunctive use of groundwater and surface water. The model determines the water allocation, reservoir operation, and groundwater pumping that minimizes the long-term sum of head-and rate-dependent groundwater pumping costs and water curtailment costs. The model is used to quantify potential savings of joint water management of the Ziya River basin in northern China, but the model can be applied to other basins as well. Estimates of natural runoff, groundwater recharge, water demands, and marginal user curtailment costs are cast into an SDP-based optimization framework. Regional and Thiem stationary drawdown is used to estimate rate-and head-dependent marginal groundwater pumping costs. The resulting optimization subproblems become nonlinear and non-convex and are solved with a hybrid GA-LP setup. A central outcome from the SDP framework is tables of shadow prices of surface and groundwater for any combination of time, inflow class, and reservoir storage. These tables represent a complete set of pricing policies for any combination of system states and can be used to guide real-time water management. Despite a significant computational demand to extract the water value tables, the method provides a suitable approach for basin-scale decision support for conjunctive groundwater and surface water management. The model provides useful insight to basin-scale scarcity-driven tradeoffs. The model outputs time series of optimal reservoir storage, groundwater pumping, water allocation, and the marginal economic value of the water resources at each time step. The Hydrol. Earth Syst. Sci., 20, 771-785, 2016 www.hydrol-earth-syst-sci.net/20/771/2016/ model is used to derive a pricing policy to bring the overexploited groundwater aquifer back to a long-term sustainable state. The economic efficient recovery policy is found by trading off the immediate costs of water scarcity with the long-term additional costs of a large groundwater head. From an initial storage at one-third of the aquifer capacity, the average costs of ending groundwater overdraft are estimated to be CNY 13.32 × 10 9 yr −1 . The long-term cost-effective reservoir policy is to slowly recover the groundwater aquifer to a level close to the surface by gradually lowering the groundwater value from an initial level of CNY 3 m −3 . Once at this sustainable state, the groundwater values are almost constant at CNY 2.15 m −3 , which suggests that wheat agriculture should generally be curtailed under periods of water scarcity. The dynamic groundwater aquifer serves as a buffer to the system and is used to bridge the water resources to multiple years. The average annual total costs are reduced by 18 % to CNY 8.56 × 10 9 yr −1 compared to a simpler formulation with fixed monthly pumping limits. The stable user's prices are suitable to guide a policy scheme based on water prices, and the method has great potential as a basin-scale decision support tool in the context of the China No. 1 Policy Document.
The Supplement related to this article is available online at doi:10.5194/hess-20-771-2016-supplement.