Simulating and understanding effects of water level fluctuations on thermal regimes in Miyun Reservoir

ABSTRACT

Water temperature dynamics in a reservoir are affected by its bathymetry, climatic conditions and hydrological processes. Miyun Reservoir in China is a large and deep reservoir that experienced a large water level decline in 1999–2004 due to low rainfall and relatively high water supply to Beijing. To study changes of stratification characteristics in Miyun Reservoir from 1998 to 2011, the one-dimensional year-round lake model MINLAKE2010 was modified by adding a new selective withdraw module and a reservoir hydrological model. Simulation results under three scenarios demonstrated that the new MINLAKE2012 model accurately predicted daily water levels and temperature dynamics during the water level fluctuation period. The water level decline led to 7.6 and 3.8°C increases in the maximum and mean bottom temperatures and about 29 days reduction in the stratification days. These simulation results provide an insight into the thermal evolution of Miyun Reservoir during the planned future water filling process.

Editor D. Koutsoyiannis Associate editor M. Acreman

KEYWORDS:

• temperature simulation
• Miyun Reservoir
• water level fluctuation
• one-dimensional model
• reservoir hydrological model
• stratification characteristics

1 Introduction

Thermal stratification is a ubiquitous phenomenon in temperate lakes (Hutchinson 1957) and reservoirs that can influence water quality and the functions of aquatic ecosystems (Owens et al. 1986). Physical (e.g., wind mixing, growth or decay of ice cover), chemical (e.g. chemical oxygen demand and nitrification) and biological (e.g., phytoplankton composition, zooplankton abundance and fish habitat) processes in an aquatic system are sensitive to water temperature distribution and thermal structure (Shimoda et al. 2011). Stratification during the summer always acts as a barrier restraining mixing of the water column and cycling of nutrients (Elçi 2008). In the period of stratification, a hypolimnetic oxygen deficit caused by thermal stratification probably reduces the metabolic rates of various organisms and thereby leads to deterioration of water quality (Wang et al. 2012).

An increase of atmospheric carbon dioxide and other greenhouse gases is projected to cause climate warming (NRC 1983, IPCC 2007). Global warming and climate change are projected to alter limnological characteristics of lakes and reservoirs, e.g., decreased ice cover period, increased surface water temperature, decreased dissolved oxygen (DO) in the hypolimnion, decreased mixing rate and increased oxythermal pressure for cold fish (Stefan and Fang 1997, Ficke et al. 2007, Fang and Stefan 2009, Jiang et al. 2012). Many studies have focused on investigating impacts of climate changes on these physical and ecological processes in individual lakes (Branco et al. 2009, Naithani et al. 2011, Shimoda et al. 2011) and lakes at regional scale (Ficke et al. 2007). However, man-made reservoirs have received less attention from water quality modellers in terms of studying the impacts of climate change on aquatic systems because complex hydrological and hydrodynamics processes in one reservoir may be different from ones in another reservoir and projections on water quality in one reservoir under future climate scenarios may not apply to any other reservoirs. Climate change projected by General Circulation Models (IPCC 2007) includes not only warmer air temperatures but also more extreme hydrological events (droughts and severe rainfall events). Climate conditions (lack of rainfall and higher air temperatures) for Miyun Reservoir in the Beijing area are possibly related to global climate warming as evidenced in some areas after 1980s (Xia et al. 2006, NRC 1983, Ma et al. 2010, Bao et al. 2012). Therefore, the study is indirectly to examine the impact of climate warming on thermal regimes in Miyun Reservoir.

In contrast to natural lakes, reservoirs typically have more complex hydrological processes related to inflow and outflow and commonly experience relatively large temporal variations in water storage and water surface elevation (Straškraba et al. 1993, Owens et al. 1998). Reservoirs may have major outlets below the surface by design for outflow whereas surface water typically flows out of natural lakes, which cause different hydrodynamics and longitudinal gradients (Straškraba et al. 1993). Therefore, thermal regimes and stratification dynamics in reservoirs are determined by not only the natural climatic conditions and bathymetry but also hydrological manipulations associated with watershed managements that may affect inflows and water demands for water supply, irrigation, etc.

Typically, in the heat balance calculation of a water body such as a reservoir, three indispensable processes, which include surface heat exchange, bottom (lakebed) heat exchange and heat flux through inflow and outflow, should be under consideration. Surface heat exchange between reservoir surface and the atmosphere (meteorological conditions) includes five heat budget components: net shortwave solar radiations H_sn, longwave radiation from atmosphere H_an, heat flux from back radiation of water body H_br, evaporation (latent heat) H_e and convection (sensible heat) H_c. All these five components are sensitive to climate conditions and have been described by Edinger et al. (1968) and Hondzo and Stefan (1993a), among others. Heat flux through inflow and outflow is essential for reservoirs although it is sometimes overlooked for thermal simulations in natural lakes by some lake models, e.g. regional lake model MINLAKE96 (Fang and Stefan 1996) and MINLAKE2010 (Fang et al. 2012a), FLAKE (Mironov et al. 2010) and LakeoneD (Joehnk and Umlauf 2001, Jöhnk et al. 2008). Climate change can modify the thermal regimes of reservoirs by not only altering surface-atmosphere interactions but also altering water budgets through precipitation and evaporation variations.

Recently, various studies have demonstrated that water level fluctuations play an important role in temperature dynamics and evolution of reservoirs, including Cannonsville Reservoir (Owens 1998b), Sooke Lake Reservoir (Nowlin et al. 2004), Lake Arancio (Naselli-Flores and Barone 2005), Lake Mead (Li et al. 2010) and Kranji Reservoir (Xing et al. 2012), among others. Water level fluctuations can result from natural hydrological change of inflows, seasonal floods and droughts and water resource management on outflows or a combination of all of these processes.

In this study, the one-dimensional (vertical) thermal lake model MINLAKE2010 was first improved by adding a new selective withdraw module and a reservoir hydrological model and then was applied to simulate long-term thermal stratification characteristics in Miyun Reservoir in China under changing water level conditions. Three modelling scenarios were developed and used to examine the ability of the new MINLAKE model for performing a continuous, long-term thermal simulation under large water level fluctuations. Long-term trends of meteorological and hydrological conditions as well as long-term characteristics of thermal regimes in Miyun Reservoir inferred from the measurements and from model simulation results were determined and analyzed. How the simulation results from the study can benefit management decisions on planned future water diversion from the Middle Routine Project (MRP) of the South-to-North Water Diversion of China was discussed. The numerical thermal simulations and predictions in Miyun Reservoir provide a good tool and information to support future adaptive management.

2 Data and methods

2.1 Study area

Miyun Reservoir is the biggest reservoir in North China and the only surface water source for the city of Beijing (Chen et al. 1998). It provides not only drinking, irrigation and industrial water for approximately 20 million people but also recreational opportunities and habitat for various aquatic organisms. Miyun Reservoir (40°25ʹN–40°38ʹN, 116°47ʹE–117°10ʹE) is located to the north of the city of Beijing and has a maximum surface area of 188 km² at the maximum pool level of 154.0 m above mean sea level (a.m.s.l.). The lowest reservoir bottom elevation is 110.5 m a.m.s.l. and the maximum water depth (H_MAX) near Baihe Dam (Fig. 1) is 43.5 m. The reservoir was created in 1960 by the interception of the Chao River and Bai River (Fig. 1). The reservoir consists of three main basins: the Bai River basin, the Chao River basin and the Inner Lake basin.

Figure 1. Map of Miyun Reservoir. The circles show the locations of the eight monitoring sites: M1 stands for Baihe station, M2 for Kuxi, M3 for ShuiJiu, M4 for Chaohe, M5 for Taoli, M6 for Henghe, M7 for Kudong and M8 for Jingou station.

Depending on inflows from upstream watersheds and water supply requirements, some reservoirs may have significant inter-annual variations in the water surface elevation. Figure 2(a) shows the time series of daily water levels in Miyun Reservoir from 1992 to 2012. Water level fluctuations in Miyun Reservoir can be divided into three periods. The first period is from 1992 to 1998, in which water level fluctuated around 150 m (H_MAX = 39.5 m). The second period is from 1999 to 2004, due to the small inflows, low precipitation and comparatively high demand on water supply, the water level decreased fast. The lowest water level was reached in July 2004, approaching 130 m a.m.s.l. The drawdown was more than 20 m, while the maximum water depth was about half of the first period. The third period is from 2005 to 2011, in which water level fluctuated around 135 m (H_MAX = 24.5 m) after a small increase in water level in 2005, but it was still much less than the water level in the first period.

Figure 2. (a) Water surface elevation from 1992 to 2011 (MSL stands for mean sea level), and (b) bathymetry curves for Miyun Reservoir (cumulative volume and horizontal area vs water level).

Figure 2(b) shows bathymetry curves (horizontal area in km² and cumulative volume in ×10⁸ m³ vs water level in m a.m.s.l.) of Miyun Reservoir and the data were obtained from the Haihe River Basin Water Information System. Based on the bathymetry, the overall stratification regime of Miyun Reservoir was first evaluated using the geometry ratio that is defined as A_s^0.25/H_MAX, where the surface area A_s is in m² and H_MAX in m. The geometry ratio is strongly correlated to the likelihood of a strong or weak stratification in a lake, which was originally proposed by Gorham and Boyce (1989) and extensively used to study stratification characteristics in regional lakes under past and future climate scenarios (Hondzo and Stefan 1993b, Fang and Stefan 2009). Polymictic lakes have high geometry ratio numbers, while strongly stratified lakes occur at the lowest numbers. The transition occurs between 3 and 5 (Fang and Stefan 1998). The geometry ratio of Miyun Reservoir increases from 2.81 m^-0.5 to 4.52 m^-0.5, when H_MAX decreases from 40 m to 20 m with corresponding changes in surface area. This means that Miyun Reservoir could have been in transition from being a stratified lake to a weakly stratified or even polymictic lake when the water level decreased from 150.0 m (1992–1998) to 135.0 m (2005–2011) a.m.s.l. This study used an improved one-dimensional (1D) model and quantified stratification characteristics of Miyun Reservoir under the variation of the water level from 1998 to 2011.

2.2 Year-round reservoir water temperature model

A 1D, year-round lake water quality model, MINLAKE2012, was developed from MINLAKE2010 that has been applied to several individual lakes (Fang et al. 2012a, 2012b) and regional lakes (Jiang et al. 2012, Fang and Stefan 2012). The MINLAKE2012 model contains a Chinese weather database including 725 weather stations in China and a user-friendly spreadsheet interface with a parameter optimization module. MINLAKE2012 has been applied to Lake Kivu in Africa (Thiery et al. 2014) and Lake Jingpo in China (Li et al. 2013).

The numerical simulation model for water temperature solves the 1D (vertical along depth) unsteady heat transport equation:

(1)

where T(z, t) is the water temperature as a function of time t and depth z, A(z) is the lake horizontal area as a function of depth z (Fig. 2 (b)), K_z is the vertical turbulent heat diffusion coefficient, H_w is the heat source term resulted from solar radiation and other factors, ρ is the density of water and c_p is the specific heat capacity of water. Description of heat source terms and the development and applications of the year-round lake model MINLAKE2010 have been presented elsewhere (Hondzo and Stefan 1993a, Fang et al. 1996, 2012a).

In order to simulate the large water level fluctuation in Miyun Reservoir, a new selective withdraw module for outflows was added into MINLAKE2012. The module simulates the physical process of the orifice outflow. It first uses total outflow, structure elevation, structure dimensions and horizontal withdrawal angle to compute upper and lower interfered limits by two criteria, i.e. Froude number (Bohan and Grace 1973) and buoyancy frequency number (Smith et al. 1987). After vertical withdrawal zone limits are calculated, velocities are determined using a quadratic shape function. Finally, outflow distribution for each layer is calculated. The detailed theoretical description of selective withdrawal can be found in other model literatures (Environmental Laboratory 1995, Cole and Wells 2004).

The MINLAKE2012 model was run with a variable horizontal layer system and a fixed time step (one day). It was assumed that water temperature is well mixed in each horizontal layer but may change from one layer to another along the depth. The variable horizontal layer system has smaller layer thicknesses for the first eight surface layers (with the layer centre depths at 0.025, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8 and 1.0 m) in order to more accurately predict snow and ice cover during the winter period (Fang et al. 1996) and thicknesses of other, deeper layers are about 1.0 m. The average Secchi disk depth used for light and solar radiation extinction coefficient computation was calculated from available measured data and equal to 3.5 m (ranging from 1.6 m to 8.3 m at the Baihe station, M1 in Fig. 1), which means that Miyun Reservoir was a mesotrophic lake classified using Secchi disk depth (USEPA 1976, Hondzo and Stefan 1993b). The extinction coefficient was then calculated using equation (2) proposed for Minnesota lakes (Hondzo and Stefan 1993a):

(2)

where z_sd is the average Secchi disk depth in m.

Water temperature modelling in Miyun Reservoir using MINLAKE2012 requires the following model input files (data): bathymetry (horizontal area versus depth, Fig. 2), model parameters and coefficients, inflow and outflow and daily meteorological data. Weather data and hydrological data for the study are presented later in a separate section. There were ice formations on Miyun Reservoir during the winter periods. Because of a lack of measured ice thickness and snow thickness data, the standard model parameters (Fang et al. 1996) for the winter module of MINLAKE2012 were used. The model calibration was conducted by adjusting model parameters/coefficients for computing surface heat flux and the vertical turbulent diffusion coefficients using the guidance of model application and calibration for MINLAKE2010 (Fang et al. 2010, 2012a). Measured water temperatures from January 2008 to December 2008 were used to calibrate the model because of the relatively smaller water level fluctuation and relatively greater amount of dates with temperature measurements in 2008. The root mean square error (RMSE) and the Nash-Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe 1970) computed between simulated and mean measured water temperatures are two model performance criteria for model calibration. After the model calibration, the calibrated model parameters were unchanged for the whole simulation period from 1998 to 2011.

2.3 Reservoir hydrological model

A reservoir hydrological model was developed and added to MINLAKE2012 by solving the water budget of the reservoir based on the principle of mass conservation, as follows:

(3)

where V is the reservoir volume (m³) as a function of time t, Q_IG is the gauged surface inflow (m³/d) from upstream rivers, Q_IU is the ungauged surface inflow (runoff, m³/d), Q_O is the total outflow (m³/d) for water supply and other uses, A_S is the surface area (km²) of the reservoir as a function of the water level (Fig. 2(b)), P is the direct precipitation (mm/d), E is the evaporation (mm/d) and G is the net groundwater source or sink (m³/d) and F is a unit conversion factor equal to 1000 for P and E in mm/d and A_s in km².

For Miyun Reservoir, G is assumed to be small and neglected, because there is no substantial groundwater into or out of Miyun Reservoir (Ma et al. 2010). However, the direct precipitation and evaporation cannot be eliminated. The evaluations over the study period (1998–2011) indicate that annual average evaporation (1218 mm) was much larger than the annual average precipitation (562 mm) (Gao 2009, Guo et al. 2011). The evaporation was calculated formulaically in the water temperature model as a function of wind speed and relative humidity (Riley and Stefan 1987).

A simple model for Q_IU was proposed:

(4)

where A_sm is the ungauged watershed area of Miyun Reservoir, C is a dimensionless runoff coefficient (a model calibration parameter) and F is the unit conversion factor used in equation (3).

Substituting equation (4) into (3), we obtain:

(5)

The time step used was one day for equation (5) and P, E and Q_IG are daily input data or calculated values. This formula was used in the MINLAKE2012 model to improve the calculation accuracy of the water budget.

2.4 Meteorological and hydrological data

Meteorological data were obtained from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/) for two weather stations (WS). Air temperature, relative humidity, wind speed, wind direction and precipitation (rainfall and snowfall) were obtained from the Surface Meteorological Dataset of China at the Miyun Weather Station (Miyun Weather Station; 40°23ʹN, 116°52ʹE, 71.8 m a.m.s.l.) and solar radiation was obtained from Beijing Weather Station (Beijing WS; 39°48ʹN, 116°28ʹE, 31.3 m a.m.s.l.) in the China radiation daily dataset. Miyun Weather Station (WS) is situated approximately 15 km south of Miyun Reservoir and it daily records air temperature, relative humidity, air pressure, wind velocity and direction at 10 m height, rainfall and snowfall. The Beijing WS is situated approximately 87 km southwest of Miyun Reservoir and it provides incoming shortwave radiation and cloud cover at a daily time interval for the model study. All necessary weather data were combined and then sewed by a Matlab program to fill the gaps of the missing data to create daily continuous meteorological data series.

The inflows and outflows of Miyun Reservoir were examined over the study period. We considered major inflows from two rivers (Bai River and Chao River, Fig. 1) and water losses due to lake surface evaporation, outflows (water demand) for a power plant (Baihe outflow) and a drinking water intake (Shuijiu outflow). We did not estimate subsurface inflows and outflows, as they are difficult to obtain and are likely to be only small contributions to the water balance. However, we estimated not-channelized surface runoff using proposed equation (4) and the calibrated runoff coefficient C was found to be 0.4.

2.5 Measured water temperatures

The measured water temperature data were obtained from the Miyun Reservoir Management Department. There are eight monitoring sites in Miyun Reservoir, i.e. Baihe (M1 in Fig. 1), Kuxi (M2), Chaohe (M3), Shuijiu (M4), Taoli (M5), Henghe (M6), Kudong (M7) and Jingou (M8). The estimated reservoir bottom elevations for these sites (based on a hydrographic contour map) are, respectively, 109.0, 112.0, 113.0, 111.5, 125.0, 117.5, 133.5 and 135.0 m a.m.s.l. Water temperatures at each monitoring site were typically measured once a month at three points along the depth (i.e. top, middle and bottom) according to the Chinese technical specifications required for monitoring surface water and waste water (Chinese Ministry of Environmental Protection 2003), but the M8 site is shallow and has only surface temperature measured. The top measuring depth was 0.5 m below the water surface, the bottom depth was 0.5 m above the reservoir bottom and the middle depth point was at half of the local water depth of each monitoring site. The water temperature measurements near the surface were made in all months throughout the year. However, the measurements for the middle and bottom depths were only made for the ice-free periods (i.e. April–October for Miyun Reservoir). The measured water temperature data are available from all monitoring sites from 1998–2004 and 2006–2011 (a total of 13 years with 204 monitoring dates) and from site M1 (Baihe station) for 2005 only.

Figure 3 shows the average measured water temperatures at the surface, middle and bottom depths from 1998 to 2011 and their standard deviations (SD) (except 2005 when measurements were only made at the M1 site). The average measured surface water temperatures for all monitoring dates were calculated from the data collected at all eight sites (M1–M8 in Fig. 1) and ranged from 0.00 to 28.66°C (SD ranged from 0.00 to 5.23°C, with a mean of 0.87°C). The average measured water temperatures at the middle depths during the ice-free periods (April–October) were calculated from the data collected at five sites (M1–M4 and M6 with similar local water depths) and ranged from 5.00 to 26.40°C (SD: 0.00–7.78°C, with a mean of 1.80°C). The average measured water temperatures at the bottom depths during the ice-free periods were calculated from the data collected at five sites (M1–M4 and M6 with similar local water depths) and ranged from 5.00 to 19.90°C (SD: 0.00–7.07°C, with a mean of 1.31°C). The measured temperatures showed good horizontal homogeneity at the surface layer, but not so good at the middle and bottom depths. However, it is still acceptable to use the 1D model to study overall stratification characteristics in Miyun Reservoir.

Figure 3. Average and standard deviation (StdDev) series of measured water temperatures at (a) the surface depth (0.5 m), (b) the middle depths and (c) the bottom depths calculated from available measurements at eight different sites M1–M8 (Fig. 1) in Miyun Reservoir for the periods 1998–2005 and 2006–2011.

3 Results and discussions

3.1 Long-term trends of meteorological and hydrological conditions

The meteorological data during 1951–2012 were investigated first. Annual mean air temperatures at the Beijing WS varied from 10.5°C to 14.1°C, with an average of 12.2°C (Fig. 4(a)). There was an obvious warming trend in Beijing from 1951 to 2011 due to the expansion of urbanization and the effects of global climate warming. The mean rate of increase was about 0.4°C/decade (R² = 0.62 for linear regression), much larger than the rate 0.074°C/decade at the global scale during the last century reported by IPCC (2007). The average air temperature in the last 10 years (2003–2012) was larger than that in the first 10 years (1951–1960) in Beijing. Figure 4(a) also shows annual air temperatures at the Miyun WS from 1989 to 2011, which are on average 1.7°C lower than those at the Beijing WS for the same time. Annual mean air temperatures at the Miyun WS varied from 10.7 to 12.4°C, with an average of 11.5°C from 1998 to 2011. The mean rate of increase was 0.089°C/decade at the Miyun WS, comparable to the IPCC’s global rate of 0.074^oC/decade. Therefore, Miyun weather was possibly affected by global climate warming and urban development in the Beijing area.

Figure 4. Inter-annual variation of (a) air temperature at the Miyun and Beijing Weather Stations, (c) total river inflow from 1990 to 2011, and (e) measured and back-calculated outflow for water supply from 1998 to 2011. Frequency distributions (histograms) of these three parameters are shown in (b), (d) and (f), respectively.

The inflow and outflow recorded during 1990–2011 were also investigated. Annual mean inflows varied from 3.7 to 35.2 m³/s, with an average of 14.4 m³/s (Fig. 4(c)), and showed a sharp decline in 1999. The average inflow of the first 9 years (1990–1998) was 22.9 m³/s, which was more than twice the average (8.5 m³/s) of the last 13 years (1999–2011). The main reason for inflow decrease was probably the much lower rainfall in the Beijing area during those 13 years. Average annual total rainfall decreased from 731 mm (1990–1998) to 555 mm (1999–2011).

Annual mean outflows varied from 8.5 to 24.7 m³/s, with an average of 13.7 m³/s (Fig. 4(e)). The average outflow of the first 6 years (1998–2003) was 19.2 m³/s, which was more than twice the average (9.5 m³/s) in the last 8 years (2004–2011). The severe decrease in outflow from Miyun Reservoir means not only less water supply for Beijing but also a great increase in the retention time and associated thermal regime modification in Miyun Reservoir. The theoretical retention time is defined as follows:

(6)

where V is the volume (m³) of water in Miyun Reservoir and Q is the average annual outflow (m³/d).

The average retention time was 1141 days (about 3 years) for 1992–1998, but 4967 days (about 14 years) for 2005–2011, which means much slower water movement in the latter period.

3.2 Evaluation of simulation results

The water temperature simulation in Miyun Reservoir was started on 1 January 1998, with an assumed initial uniform water temperature profile of 4^oC, which was a good estimation of the initial temperature profile as shown and indicated later by model results (Figs 6–8). The RMSE between simulated and mean measured surface water temperatures was 2.26°C in the calibration year of 2008. The RRMSE was 0.15 when RMSE was normalized by mean observed surface temperature. The NSE coefficient for simulated and measured water level was 0.94, which indicates an accurate water balance computation by the proposed model (equation 3). The calibrated wind sheltering coefficient was 1.0, which is the same value as recommended for large surface lakes (Hondzo and Stefan 1993a), and the multiplier for calibrating turbulent diffusion coefficient in the metalimnion was determined to be 1.0, which is larger than 0.5 recommended for deep natural lakes (Fang et al. 2012a). This is because vertical mixing in the reservoir with inflows and outflows is stronger than the mixing in deep natural lakes without inflows and outflows.

Figure 5. (a) Time series of measured and simulated water levels under the scenarios S1, S2 and S3; (b) annual water inputs (direct precipitation, ungauged runoff and river inflow) and losses (outflow and evaporation) in Miyun Reservoir from 1998 to 2011.

Figure 6. Simulated water temperatures at the surface depth (0.5 m) under (a) scenario S1, (b) scenario S2 and (c) scenario S3 compared against mean measured temperatures (including average values ±1 SD) in Miyun Reservoir for 1998–2011.

Figure 7. Simulated water temperatures at the middle depths under (a) scenario S1, (b) scenario S2 and (c) scenario S3 compared against mean measured temperatures (including average values ±1 SD) in Miyun Reservoir for 1998–2011.

Figure 8. Simulated water temperatures at the bottom depths under (a) scenario S1, (b) scenario S2 and (c) scenario S3 compared against mean measured temperatures (including average values ±1 SD) in Miyun Reservoir for 1998–2011.

Three simulation scenarios (Table 1) were established to examine the effects of water level fluctuation on water balance and thermal simulation results. The first scenario (S1) was used to represent the baseline condition, in which the reservoir hydrological model (equation 3) was not included. The inflow and outflow processes were included for scenario S1 but daily outflows were back-calculated to keep the water level unchanged during the simulation period by considering inputs from measured inflows and direct precipitation and water lost through evaporation. The initial water level on 1 January 1998 was 149.63 m and stayed the same throughout the simulation from 1998 to 2011 for scenario S1. Scenario S2 includes the inflow and outflow processes (measured inflows from major rivers and selective withdrawal for outflow), but without considering the ungauged surface inflow. Scenario S3 includes the ungauged surface inflow (runoff) estimated for ungauged watershed areas using precipitation and runoff coefficient. All three scenarios S1, S2 and S3 were run using the model MINLAKE2012 that was improved and enhanced from MINLAKE2010 for the Miyun Reservoir study.

Table 1. Description of different simulation scenarios.

Simulation scenarios	Description	With inflow and outflow	Reservoir hydrologic model
S1	Baseline (outflows calculated to keep water level unchanged)	Yes	No
S2	Include inflow and outflow process	Yes	No
S3	Include reservoir hydrologic model	Yes	Yes

Simulated water levels in Miyun Reservoir from 1998 to 2011 under the three scenarios are shown in Fig. 5(a). The mean absolute error between the measured and assumed water level from scenario S1 is 11.54 m and the standard deviation from the mean is 5.08 m. Scenario S1 was not designed to accurately simulate the water level fluctuations in Miyun Reservoir but was instead used to examine the impact of the water level decline on the simulation of water temperature and thermal stratification characteristics when results from scenario S3 were compared with the ones from scenario S1. The mean absolute error between the measured and simulated water level from scenario S2 is 1.80 m and the standard deviation from the mean is 1.03 m. The mean absolute error between the measured and simulated water level from scenario S3 is 1.15 m, which is less than 3% of the initial maximum water depth (39.13 m). The NSE coefficients between the measured and simulated daily water levels from scenarios S2 and S3 are 0.89 and 0.93. This means that the reservoir hydrological model improved the model accuracy of water level simulation in Miyun Reservoir using MINLAKE2012. Figure 5(b) shows magnitudes of the five components for water budget computation in each year during the simulation period. Inflows from major rivers and the outflow for water supply are major contributors for the water level fluctuations in Miyun Reservoir. There were nine years that had a negative water balance (more water loss through the outflow of water supply and evaporation than water input to the reservoir); and there were continuous negative water balances from 1999 to 2003 (over five years), which resulted in the lowest water level in July 2004. This was because the demand of water supply for those five years was significantly larger than the water input to Miyun Reservoir. Over the 14 simulation years, river inflows, direct precipitation and ungauged runoff are 82.7%, 13.2% and 4.1% of the total water input; and evaporation is 16.4% of the total water loss. The total water loss over the 14 years was 17.2 × 10⁸ m³ and there was 20.0 × 10⁸ m³ water loss from 1999 to 2003. The back-calculated annual outflows for scenario S1 are shown in Fig. 4(e) for comparison with measured outflows for water supply. In order to keep the water level in 1998 unchanged, available water for water supply to Beijing from 1999 to 2003 would have to be significantly lower than measured outflow. The demand for water supply from 2004 to 2011 was reasonable in comparison to water input to the reservoir.

Figures 6–8 show the time series of simulated water temperatures under the three scenarios (the top panel for S1, the middle panel for S2 and the bottom panel for S3) at three different water depths in Miyun Reservoir: surface, middle and bottom d, respectively. Mean measured water temperatures calculated from the monitoring sites ±1 SD are also given in Figs 6–8 for comparison of simulated temperatures. In general, simulated water temperatures from scenario S3 have a closer agreement with observed temperatures. In addition to computing and reporting the ±1 SD RMSE between simulated and mean measured water temperatures, the percentages of simulated temperatures that are within plus and minus one, two and three standard deviations from the mean observed temperatures were also determined and are presented in Table 2.

Table 2. Error parameters computed between measured and simulated temperatures under scenarios S1, S2 and S3.

Simulation scenarios	RMSE (°C)^1			±1×SD (%)^2			±2×SD (%)^2			±3×SD (%)^2
	S1	S2	S3	S1	S2	S3	S1	S2	S3	S1	S2	S3
Surface	1.80	1.83	1.80	40	40	41	59	62	61	78	72	70
Middle	6.37	2.81	2.93	15	45	43	35	70	70	48	78	77
Bottom	4.60	4.02	3.95	19	30	32	41	55	51	55	70	73

¹RMSE was calculated between simulated temperatures and mean measured temperatures derived from available data at the monitoring sites.

²Percentages of simulated temperatures that are within plus and minus 1, 2 and 3 standard deviations (SD) from the mean measured temperatures.

For simulating water temperatures in the surface depth (0.5 m) over 14 years, the differences among the three scenarios are very small; it is hard to see any difference of simulated surface temperatures between two scenarios from Fig. 6. The RMSEs simulated from scenarios S1, S2 and S3 are 1.80°C, 1.83°C and 1.80°C, respectively. The absolute differences between simulated daily surface temperatures between two scenarios were calculated; 75% of the absolute differences of surface temperatures between scenarios S2 and S3 are less than 0.27°C. Between scenarios S1 and S2 or S3, 75% of the absolute differences of surface temperatures are less than 0.67°C with mean absolute differences of 0.58°C. This means that simulated surface water temperatures from scenario S1 with assumed constant water level (using back-calculated daily outflows) from 1998 to 2011 are not exactly the same but similar or close to the ones from scenarios S2 and S3 that have a much better match with observed water levels. This is reasonable, because the thermal balance in the surface layer is mainly controlled by weather conditions and almost not influenced by the total water depth for relatively deep reservoirs.

For simulated water temperatures in the middle depths over 14 years, the RMSEs simulated from scenarios S1, S2 and S3 are 6.37°C, 2.81°C and 2.93°C, respectively. For scenario S1 (assumed constant water level), except the first three years (1998–2000), simulated middle-depth temperatures are significantly different from mean observed ones (Fig. 7(a)), which explains the largest RMSE. For scenarios S2 and S3, most of the simulated middle-depth temperatures (78% and 77%, respectively) are within ±3 SD of the corresponding mean measured temperatures, but only 48% for scenario S1 (Table 2). The mean absolute difference of middle-depth temperatures between scenarios S2 and S3 is only 0.27°C; and between scenarios S1 and S2 or S3 is about 2.9°C (25% and 50% of the absolute differences are greater than 5.2 and 1.8°C, respectively). Therefore, assumed constant water level for scenario S1 did not produce accurate temperature predictions at the middle depths in Miyun Reservoir after the year 2000. Water temperatures in the middle depths in a relatively deep reservoir such as Miyun Reservoir are strongly influenced by thermal stratification that is controlled by the total water depth. The maximum water depths in 1998 and 1999 were about 40.0 m and the middle depths were about 20.0 m at the monitoring sites M1–M4 and M6. Both measured and simulated water temperatures between 1998 and 1999 are a few degrees lower than those in other years. This is because the maximum water depths from 2006 to 2011 were about 26.0 m and the middle depths were about 13.0 m from the water surface at the monitoring sites M1–M4 and M6.

For simulated water temperatures in the bottom depths over 14 years, the RMSEs simulated from scenarios S1, S2 and S3 are 4.60°C, 4.02°C and 3.95°C, respectively. Table 2 shows that 55%, 70% and 73% of simulated bottom temperatures from scenarios S1, S2 and S3, respectively are within ±3 SD of the corresponding mean measured temperatures. The mean absolute difference of bottom-depth temperatures between scenarios S2 and S3 is 0.31°C. Simulated middle-depth and bottom water temperatures (Figs 7 and 8) show steady increases in the spring and summer months because of solar heating and decreases in the autumn due to cooling. Mean measured middle-depth temperatures from different monitoring sites (Fig. 7) in all the years except 2000 clearly display a similar magnitude of temperature increases due to solar heating. Figure 8 clearly shows smaller measured temperature increases in 1998–2000 when the reservoir was relatively deep (about 40 m) but much larger temperature increases after 2000 when the maximum water depth decreased to about 20 m. Mean measured bottom temperatures from different monitoring sites (Fig. 8) in 1999, 2000, 2002 and 2005 did not show much temperature increase; it is possible that these measurements were less accurate than other measurements because it might have been difficult to locate the bottom depth, or the water quality probe got into the mud at the bottom during data collection. Figures 7 and 8 show a gradual increase in the middle-depth and bottom water temperatures from 2000 to 2004 (maximum temperatures) that correspond to a decrease in the water level from about 150 to 130 m a.m.s.l (Figs 2(a) and 5(a)). Model accuracy on predicting temperatures at middle and bottom depths for scenario S3 was not high; this is because the influences of horizontal circulation patterns and inflow/outflow effects cannot be coped with in a 1D model such as MINLAKE2012. These horizontal processes are only possibly coped with using two-or three-dimensional models (Joehnk and Umlauf 2001).

3.3 Long-term characteristics of thermal regimes

Nine water temperature characteristics (maximum and mean surface water temperatures, maximum and mean bottom water temperatures, maximum and mean water temperature differences, total stratification days, maximum temperature DOY and temperature increase slope) defined in Table 3 were extracted from the daily simulation results. The simulated values of these characteristics for the period 1998–2011 under scenario S3 are listed in Table 4. Figures 9 and 10 illustrate the response of thermal characteristics to water level fluctuations. Simulated annual values of long-term characteristics under the S1 scenario are also presented in Figs 9 and 10; these simulation results show what long-term characteristics could be if there were no water level changes from 1998 to 2011. Results under the S1 scenario only reflect climatic variations from 1998 to 2011 (see Fig. 4(a) for air temperature variations in the Miyun WS), while results under the S3 scenario reflect effects from both water level changes and climatic variations.

Table 3. Parameters used to define long-term water temperature characteristics in Miyun Reservoir.

Parameter	Description
Maximum surface temperature (Max. Surf.)	Annual maximum daily water temperature at 0.5 m below the water surface
Mean surface temperature (Ave. Surf.)	Annual mean daily water temperature at 0.5 m below the water surface
Maximum bottom temperature (Max. Bot.)	Annual maximum daily water temperature at 0.5 m above the bottom
Mean bottom temperature (Ave. Bot.)	Annual mean daily water temperature at 0.5 m above the bottom
Maximum temperature difference (Max. Dif.)	Annual maximum daily water temperature difference between surface and bottom
Mean temperature difference (Ave. Dif.)	Annual mean daily water temperature difference between surface and bottom
Total stratification days (Str. Days)	Total number of days during the ice-free period when difference between surface and bottom water temperature is larger than 1°C
Maximum temperature day of year (Max. DOY)	Date (day of the year) of the maximum water temperature at the middle or bottom depths
Temperature increase slope (Sti)	Line regression slope of water temperature with respect to DOY during the period from April 1 to the date of maximum water temperature at the middle or bottom depths

Table 4. Simulated annual values of seven parameters (Table 3) used to define long-term water temperature characteristics from 1998 to 2011 under scenario S3.

Year	1998	1999	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
Max. Surf.	29.34	31.64	31.14	28.89	30.35	28.90	29.20	29.96	29.05	29.15	30.31	29.75	31.71	30.01
Ave. Surf.	15.00	14.54	14.35	14.42	14.21	14.36	14.78	14.69	14.73	15.39	14.73	15.04	14.31	15.22
Max. Bot.	13.21	12.07	15.22	16.37	18.10	19.62	19.65	19.33	19.56	20.39	19.69	19.92	19.37	19.67
Ave. Bot.	7.68	6.62	7.95	8.88	9.23	10.08	10.43	10.31	10.67	10.46	10.34	10.05	10.00	10.46
Max. Dif.	19.61	24.17	21.85	18.67	17.04	15.4	15	17.34	14.84	15.93	14.97	15.71	17.62	16.03
Ave. Dif.	7.32	7.93	6.41	5.54	4.98	4.28	4.35	4.38	4.06	4.93	4.39	4.99	4.31	4.76
Str. Days	240	221	208	206	207	190	206	199	213	205	204	204	197	191

Figure 9. Simulated annual maximum and mean (a) surface water temperatures, (b) bottom water temperatures and (c) differences of temperatures between the surface and bottom depths; and (d) simulated stratification days (Str. Days) under scenarios S1 and S3 in Miyun Reservoir (1998–2011).

Figure 10. (a) Day of the year (DOY) having the maximum water temperatures in the middle depths and bottom depths, (b) slopes of temperature increases during the period from April 1 to the day of maximum water temperatures at the middle depths and bottom depths simulated under scenarios S1 and S3.

Water level fluctuations had different effects on water temperatures at the surface (0.5 m below) and bottom (0.5 m above) depths (Fig. 9). For simulated water temperatures at 0.5 m below the water surface, the annual mean values for 1998–2011 under scenarios S1 and S3 ranged from 13.91 to 15.39°C. The differences in annual maximum and mean surface temperatures between scenarios S1 and S3 are less than 0.36 and 0.31°C, respectively, which suggests that the water level decline (1999–2004) has little influence on the surface temperatures (Fig. 9(a)). The annual cycle of the surface temperature was basically the same regardless of the large water level fluctuation (Fig. 6(c)) and was controlled and followed by the change (cycle) of air temperature and solar radiation.

The decline in water level (1999–2004) leads to an obvious rise of simulated water temperatures at the bottom depths (Fig. 9(b)). The maximum and mean bottom temperatures under scenario S1 were in the ranges 11.46–13.94°C and 6.46–7.71°C, respectively. These relatively small ranges (2.48 and 1.25°C) indicate that the inter-annual variations of weather factors would cause a small change of water temperature at the bottom depths in the deep Miyun Reservoir if the water level stays unchanged. The maximum and mean bottom temperatures under scenario S3 were in the ranges 12.07–20.39°C and 6.62–10.67°C, respectively. The much larger ranges (8.32 and 4.05°C) indicate that the water level variations had an obvious influence on water temperature at the bottom depths. For the water level decline period (1999–2004), the lower the water level was, the higher the bottom temperature was. The increases in the maximum and mean bottom temperatures from 1999 to 2004 reached 7.58°C and 3.81°C, respectively. Both increases exceed 50% of the corresponding mean bottom temperatures in 1999.

The maximum and mean temperature differences between surface and bottom under scenario S1 were in the ranges 19.91–24.32°C and 7.24–8.50°C from 1998 to 2011, with average values of 22.03 and 7.73°C, respectively when water level is assumed unchanged. The maximum and mean temperature differences under scenario S3 were in the ranges 14.84–24.17°C and 4.06–7.93°C, with average values of 17.44 and 5.19°C, respectively. The smaller temperature differences occurred when the water depths in Miyun Reservoir were smaller after 2003. The average values of maximum and mean temperature differences from 1998 to 2011 had 21% and 33% decreases, while the water depth had an average 28% decrease by comparing the corresponding results between scenarios S1 and S3. This means that the strength of the stratification weakened significantly because of the water level decline in Miyun Reservoir.

The stratification period (days) was determined to sum the days when the temperature differences between surface and bottom were larger than 1°C (Stefan et al. 1996). The duration of the stratification shortened as well when the water levels reduced from 2000 to 2011 (Fig. 9). The stratification days in 1998 and 1999 under the S3 scenario are almost the same as the days under the S1 scenario because the assumed constant water level was similar to the measured water levels in these two years (Fig. 5(a)). The stratification days under scenario S1 increased from 226 to 240 days, but under scenario S3, they decreased from 221 to 206 days (Fig. 9(d)) from 1999 to 2004 (the period of sharp decline in water level). The stratification days under scenario S1 increased by 14 days, which were probably due to climatic warming (Fig. 4(a)) because the water level was constant. The stratification days under scenario S3 decreased by 15 days, but considering the increase due to climatic conditions, the decrease due to the water level decline would only be about 29 days (just slightly less than one month). The decrease of the stratification periods was more dependent on the timing of the fall turnover, rather than the onset of stratification in spring. The standard deviations of the DOY (day of the year) of the end time and the onset of stratification are 14 days and 7 days, respectively. This result indicates that the water level fluctuation has more significant influence on the destratification in autumn than the stratification formation in spring. This conclusion is consistent with the study conducted in Cannonsville Reservoir (Owens 1998a). This is because a lake or reservoir has weak stratification in the summer and early fall and is more easily mixed by wind when its water depth becomes smaller due to the water level decline.

Figure 10(a) shows the maximum temperature DOYs, which is the day of the year when simulated maximum water temperature occurs at either the middle or bottom depths under scenarios S1 and S3. The maximum temperature DOY at the middle depths under scenarios S1 and S3 were in the ranges 283–297 and 244–285, with average values of 291 and 260 (SD of 4 and 12 days), respectively. Climatic variations (S1 scenario) only caused a small variation (14 days) on the day of maximum temperatures at the middle depths but water level fluctuations resulted much larger variation (41 days). The temperature at the middle depths reached the maximum value more than a month early (average 31 days, from late October to late September) when the reservoir had a lower water depth. The maximum temperature DOY at the bottom depths under scenarios S1 and S3 were in the ranges 305–322 and 270–318, with average values of 310 and 283 (standard deviations SD of 7 and 14), respectively. The temperature at the bottom depths also reached the maximum value about a month early (average 27 days, from early November to early October) when the water depth in Miyun Reservoir was decreased.

Figures 7 and 8 show clearly certain linear temperature increases occurring at the middle and bottom depths during the summer periods. The period of water temperature increase was assumed to start on April 1 and end on the day of maximum water temperature at the middle or bottom depths. A linear relation between water temperature and corresponding DOY was developed for each simulation year. All correlation coefficients of linear regressions between simulated water temperature and corresponding DOY at the middle and bottom depths under scenarios S1 and S3 were larger than 0.995 (SD < 0.004). Figure 10(b) shows temperature increase slopes during the summer periods calculated at the middle and bottom depths under scenarios S1 and S3. Temperature increase slopes at the middle depths under scenarios S1 and S3 were in the ranges 0.052–0.082°C/d and 0.068–0.115°C/d from 1998 to 2011, with average values of 0.062 and 0.100°C/d (SD of 0.007 and 0.014°C/d), respectively. Based on the average maximum temperature DOY (260) under the S3 scenario, the period of the temperature increase is 170 days (from 1 April), therefore, the average temperature increase at the middle depths is about 17°C (170 × 0.1°C/d). Temperature increase slopes at the bottom depths under scenarios S1 and S3 were in the ranges 0.037–0.044°C/d and 0.040–0.092°C/d, with average values of 0.040 and 0.076°C/d (SD of 0.002 and 0.018°C/d), respectively. Similarly, the average temperature increase at the bottom depths is about 14.7°C [(283–90) × 0.076°C/d]. Water temperature on 1 April in a temperate lake is about 4–5°C and then the maximum water temperature at the bottom depths is about 19–20°C, which matches the magnitude of the maximum bottom temperatures shown in Fig. 9(b) for the S3 scenario.

Temperature increase slopes (S_ti) from scenario S3 show good linear relationships with water depth (H_MAX). At the middle depths, the linear regression equation is S_ti = −0.0022H_MAX + 0.1608, with a correlation coefficient R² = 0.80 and at the bottom depths, S_ti = −0.0032H_MAX + 0.1632, with a correlation coefficient R² = 0.95. These regressions show that the slopes of temperature increase become larger when the maximum depth (H_MAX) is small. It should be mentioned that in general the slope of temperature increase is dependent not only on the maximum depth but also on the available solar radiation (geographic location of a reservoir) and the radiation attenuation coefficient that is affected by water transparency or water quality. In general, the results suggest that water temperature in the hypolimnion increased and total stratification days decreased with the water level decline in Miyun Reservoir from 1998 to 2011.

3.4 Thermal regime changes due to future water level increase

The outflow from Miyun Reservoir for water supply to the capital city of Beijing was kept at a very low level since 2004 (Fig. 4(e)). Because of the steady population increase and urban/industrial development in Beijing, the water demand actually did not decrease instead of increase. In order to solve the water shortage in Beijing and other cities in North China, the Middle Routine Project (MRP) for the South-to-North Water Diversion was planned and designed for many years and now is under construction. When the MPR is completed in 2015, the water diverted from Danjiangkou Reservoir on Hanjiang River (a tributary of Changjiang River) will be sent to Beijing and other cities in North China. At that time, Miyun Reservoir will be used to store excess water (beyond the water supply need) temporarily in the first 5 years, i.e. from 2015 to 2019.

The water from the MRP will be pumped into the Miyun Reservoir through the diversion tunnel of the power plant; its entrance or discharge point is near the Baihe Dam (Fig. 1) and at 125.5 m a.m.s.l. According to the design and reservoir bathymetry (Fig. 2(b)), the water from the MRP will be pumped into the Miyun Reservoir in a deeply submerged way at a flow rate of 10 m³/s. The water level in Miyun Reservoir will rise from about 135 m (the average water level in the last eight years, shown in Fig. 2(a)) to 150 m a.m.s.l., and this water filling process will be just an inverse process of water level decline from 1999 to 2004.

The model simulation results from 1998 to 2011 in Miyun Reservoir have indicated that the stratification characteristics were affected by the water level decline, which was caused by low rainfall, small river inflow and reservoir operations to fulfil the water supply need. The water level decline had a greater impact on the evolution of thermal regimes. The results presented in this paper provide useful information and insights to water managers and policy makers in China about what will/may happen to Miyun Reservoir after the completion of the MRP water diversion. The 1D numerical hydrothermal model can provide basic information to support analysis and prediction of water quality and the ecosystem in Miyun Reservoir and provide a good tool to support adaptive management. To understand the impact of the deeply submerged discharge from water diversion on nearby thermal regimes especially in the hypolimnion and sediment oxygen demand in Miyun Reservoir, two- or three-dimensional modelling will be useful and necessary.

4 Conclusions

In this paper the one-dimensional year-round lake water quality model MINLAKE2010 was improved by adding a new selective withdrawal module and a reservoir hydrological model in order to simulate water temperature dynamics under the large water level fluctuation in a lake or reservoir. The enhanced model MINLAKE2012 model was then used to investigate thermal evolution in Miyun Reservoir during a 14-year study period (1998–2011) when the water level declined from 150 m to about 135 m a.m.s.l. The simulation results show that the improved model MINLAKE2012 accurately predicted water level (NSE = 0.93) and temperature dynamics in Miyun Reservoir under variable water levels. The RMSE between simulated and mean measured surface water temperatures is 1.79°C. Model simulations were performed under three scenarios (Table 1); results under scenario S1 with assumed constant water level were compared with results under scenario S3 and this helps us to understand the direct impact of water level declines on thermal regimes in Miyun Reservoir.

Long-term characteristics of thermal stratification and regimes in Miyun Reservoir were extracted from daily simulation results. It was found that the water level decline had different effects on water temperature characteristics in the surface depth (0.5 m below the surface), the middle depths and the bottom depths (0.5 m above the bottom). The water level decline had little influence on the temperature in the epilimnion and the annual cycle of the surface temperature because they are primarily controlled by climatic conditions. The decline of water level led to steady rises of water temperature in the middle and bottom depths. The lower the water level was, the higher the temperature in the hypolimnion was. The increases of the maximum and mean bottom temperatures from 1999 to 2004 reached 7.59°C and 3.73°C, respectively. Both the strength and duration of the stratification weakened with water level decline. After considering the increase of the stratification days due to climatic conditions (from the S1 scenario), the decrease in the stratification days due to water level decline would only be about 29 days (just slightly less than one month) in Miyun Reservoir. The maximum temperature differences in Miyun Reservoir decreased from 23.91 to 14.88°C due to the water level decline. These changes in thermal stratification characteristics would affect the water quality and the ecosystem in the reservoir which will be studied at the next step of the project.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (No.51039002 and No.51279078) and Tsinghua University Initiative Scientific Research Program (No. 20121088082).