Thermal and Hydrodynamic Changes under a Warmer Climate in a Variably Stratified Hypereutrophic Reservoir

Abstract

We quantified effects of future climate warming on temperature and stability in a variably stratified, hypereutrophic reservoir with large fluctuations in water level by calibrating a 2-D model (CE-QUAL-W2, version 3.7.1, Portland State University, Portland, USA) of reservoir hydrodynamics using a time series (1992 to 2011) of inflow and air and water temperature. The model was then forced with increased air temperature projected by an ensemble of climate models that accounted for complex local topography and seasonality, with greater warming in summer. Warming increased annual evaporation rates by 2.6 to 7.9%. Water temperature increased by 0.44 (whole-reservoir; p < 0.05), 0.47 (epilimnion; p < 0.01), and 0.30 °C (hypolimnion; p < 0.05) per 1 °C increase in air temperature. Thickness of the epilimnion and hypolimnion diminished, with expansion of the metalimnion. Schmidt stability correlated with mean water depth over a wide range of depths (3.9 to 8.1 m; Adj. R² = 0.91 to 0.93; p < 0.001). Increased air temperature increased mean annual stability by 6.1 to 23.6 J m⁻² when depth was large and the reservoir stratified, but when depth was low (due to combined low inflow and, in preceding years, high withdrawals), inhibiting stratification, then water temperatures increased evenly (and more) throughout the vertical profile so change in mean annual stability was near zero (−0.1 to 1.1 J m⁻²). Combined effects of reservoir management (volume, timing, and elevation of water withdrawal) and climate warming (temperature of air and benthic sediment) can impact the hydrodynamic regime differently under variably stratified conditions with implications for release of phosphorus from sediment and vertical transport of phosphorus to the euphotic zone.

1. Introduction

Algal production, anoxia, and subsequent release of nutrients from benthic sediment are strongly controlled by timing, strength, and duration of thermal stratification. When a reservoir is stratified in summer, influx of dissolved oxygen (DO) from the atmosphere to the warmer epilimnion (with lower DO saturation) or from plant photosynthesis to the hypolimnion is limited and anoxia can occur when DO demand exceeds supply. Under anoxic conditions, ferric ions of iron reduce to the ferrous form, which leads to the dissolution of the phosphorus and iron bond and causes sediment to release soluble reactive phosphorus, iron, ammonia, and manganese to the hypolimnion [1]. Liberated phosphorus can be entrained to the euphotic zone through sustained vertical eddy diffusion, occasional entrainment, or annual mixis, resulting in increased algal growth. Decay and settling of algae into the hypolimnion and sediment feed microbial bacteria that have high DO demand, reinforcing a positive feedback.

Stratification, which is controlled in part by climate, has been altered in many systems by climate warming between 1970 and 2010 [2], and may be further impacted by future warming. Mean surface air temperature is projected to increase in the future [3], which may have opposing effects on stratification. Increased air temperature (and associated longwave radiation that enters a reservoir) heats a reservoir disproportionately at the surface, intensifying the vertical temperature gradient, increasing stability, extending the duration of stratification, and resulting in stratification at lower water levels [2,4,5,6,7,8], all of which increase the potential for anoxia and release of phosphorus from sediment [9]. Alternatively, increased evaporation rates and/or decreased inflow can lead to a long-term decrease in reservoir water depth, which can decrease stability and inhibit stratification if the reservoir is sufficiently shallow [6,10,11]. Most studies in the literature have investigated effects of climate warming on deeper reservoirs and lakes, which are monomictic or dimictic, but not variably stratified water bodies that are regularly or occasionally polymictic and have large intra-annual fluctuations in water level [2,4,5,6,7,12,13,14,15,16,17,18].

The objective of this study was to quantify the effects of a projected future warmer climate on thermal and hydrodynamic regimes in a reservoir that was intensively managed and became relatively shallow so that it ranged from being strongly to weakly stratified. Here, dam operators made decisions for volume, timing, and elevation of water withdrawal, which contributed to intra- and inter-annual variability in water depth and circulation. In turn, this complicated the determination of the parameters of the air temperature-stability relationship and the estimation of the impact of future warming on stability and duration of stratification. We calibrated a 2-D simulation model of reservoir hydrodynamics over a study period from 1992 to 2011, then forced the calibrated model with an ensemble of climate projections. We accounted for the effect of warming on other terms and feedbacks to the model in vertical (air temperature, precipitation, inflow, and benthic sediment) and lateral (inflow) directions.

The research questions were: What are the impacts of increased air temperature on water temperature, evaporation, and stability? How does warming impact the depth-stability relationship? How do impacts of warming differ when the reservoir is monomictic compared to polymictic? Are there seasonal differences in impacts of warming? We expected increased temperatures in all strata, with greater warming at the surface, which would increase stability for warmer climates. However, due to increased evaporation in a shallow reservoir with large surface area and concomitant lower water depth due to water withdrawal and management, we expected to see inhibition of stratification in some years of already low water level. We also expected greater magnitudes of change in water temperature and stability in summer compared to winter, given the extended period of high air temperature, high evaporation rates, and low precipitation rates typical during Mediterranean summers.

The study is of general interest as it tests a modeling framework for understanding how intensively managed and variably stratified reservoirs may respond to climate change. Intensively managed systems that have large seasonal changes in water level due to variations in both inflow and withdrawals may differ in their temperature-stability relationships compared with deep and natural lake systems. Management of intensively managed reservoirs is critical for municipal water supply, and understanding the impacts of climate change on their hydrodynamics and water quality is important for understanding potential climate-related risks to urban water resources.

2. Materials and Methods

2.1. Study Site

Hodges Reservoir (hereafter Hodges; Figure 1) supplies water for municipal use and is 50 km north of the city of San Diego, United States, and 20 km inland from the coast at an altitude of 70 m a.s.l. The reservoir basin is relatively shallow (8.5 m mean depth at capacity) compared to reservoirs in the region [11] and there is a relatively large interface between the hypolimnion and sediment, making it vulnerable to hypolimnetic anoxia and subsequent release of phosphorus from sediment [11,19]. At capacity, surface area is 6.0 km² and maximum depth is 29.1 m. In most years, it is warm monomictic, but can be continuous warm polymictic when decreased inflow and/or increased drawdown inhibit stratification. Water level fluctuations were large over the study period, both intra-annually (mean amplitude = 3.5 m and maximum amplitude = 13.7 m) and inter-annually (maximum amplitude = 17.2 m), due to reservoir management, so stratification was controlled primarily by water depth and secondarily by air temperature [11].

Hodges is classified as hypereutrophic, as it had a high mean Anoxic Factor (AF; 157 day year⁻¹), the number of days in a calendar year for which a sediment area equal to the surface area of a reservoir was overlaid with anoxic water (<1 mg DO L⁻¹) [19], and the lowest Secchi disk transparency depth (1.0 m) relative to three other major water supply reservoirs in San Diego County [11]. Barrett, Morena, and Sutherland reservoirs had mean annual values of 161, 60, and 149 day year⁻¹, respectively, for AF and 1.8, 1.7, and 1.2 m, respectively, for Secchi disk transparency depth over the study period. At Hodges, the mean hypolimnetic DO concentration was 0.1 mg L⁻¹ and the reservoir became almost entirely anoxic due to algal production, decay, and related metabolic processes during mixis in autumn. The AF increased and transparency decreased while the frequency of winter events with near complete anoxia increased over the study period, which may be due in part to increased expansion of urban development in a watershed with significant agricultural land use [11].

Of all water entering Hodges between 1929 to 2010, local runoff accounted for 89%, imports for 5%, other inflows for 3%, and direct precipitation for 3%. There are four inflows to Hodges, all within close proximity of each other. The San Dieguito River contributes 92% of local runoff. It drains most of the Hodges watershed, which consists of mixed land use, including intensive and extensive agricultural, low-density residential, and open space. The balance (8%) of the local runoff to Hodges is conveyed by three streams (Felicita, Green Valley, and Kit Carson Creeks) draining smaller urban watersheds.

2.2. Study Design

We used CE-QUAL-W2 (hereafter W2; version 3.7.1, Portland State University, Portland, USA [20]), to simulate reservoir evaporation, water level, temperature, stability, and stratification from 1992 to 2011. These variables are impacted by simulations of the energy balance of the reservoir done by W2 using incoming solar radiation and longwave radiation estimated from air temperature, cloud cover, humidity, and wind data. We tested the sensitivity of the model to warming using the delta method, which adds a given temperature increase (∆T) to the observed historical time series of air temperature and uses the newly adjusted air temperature series as input to the model [21]. Our results include a reference run with no change, and runs for four warming scenarios.

Higher air temperatures increase conduction of heat from the atmosphere, and increase the temperature of precipitation, inflow, and sediment, all of which affect the energy balance of the reservoir. We accounted for higher temperatures of air, precipitation, inflow, and sediment in input data. We assumed relative humidity remains constant under a warmer climate; temperature of precipitation was set equal to daily mean air temperature; temperature of inflow was set equal to mean monthly air temperature over the study period [22,23]; and temperature of sediment was set equal to mean historical air temperature over the study period [20,24,25]. Although air temperature may not be a good estimate for temperature of precipitation due to a differential in elevation between where air temperature is measured and where precipitation forms, we assumed this assumption was justified due to the small proportion (3%) of precipitation relative to the total inflow that entered the reservoir. We set the temperature of sediment equal to the historical mean rather than to the mean of projected increased air temperature to allow for a temporal lag between air temperature and sediment temperature changes.

We did not consider explicitly how volume of precipitation or inflow may change due to climate warming because of a lack of consensus in the direction and magnitude of change in the literature. However, ranges of precipitation depth and inflow volume over the study period (1992 to 2011) were wide enough to be representative of the overall long-term historical record (1929 to 2013; Appendix A.1; Figure A1). Furthermore, almost all large precipitation events resulted in large runoff volumes, most of which (>90%) was spilled over the dam. Observed reservoir mean depth ranged from 3.6 to 8.5 m over the study period, including periods of overflow and an extended drought (2002 to 2004) when water levels were at record lows. Therefore, we simulated the response of the reservoir to both forecast climate (air temperature) scenarios and to a wide range of precipitation, runoff, and operation conditions.

2.3. Model Input and Calibration Data

The W2 model requires (1) daily time series of the volume and temperature of reservoir inflow and outflow (Appendix A.1); (2) bathymetric data to create depth-area-volume relationships (DAVR) for each model segment (Appendix A.2); and (3) daily time series of meteorological data (Appendix A.3). Water temperature profiles used for model calibration were taken at the deepest point near the dam at approximately 1 m depth intervals throughout the water column. Water temperature and Secchi disk transparency depth were measured monthly from January 1992 to November 2011 by City of San Diego.

2.4. Model Parameterization and Calibration

In low water years (2002 to 2004), when mean annual maximum depth, Z, <14.3 m and mean annual mean depth, $\overline{Z}$, <4.5 m, observed stability was low and stratification did not develop. When water depths are low, sediment temperature becomes more tightly coupled with air temperature [24], most markedly during the spring and summer [26]. Our preliminary model runs showed that setting the sediment temperature equal to the mean of the historical daily mean air temperature artificially decreased the hypolimnetic temperature and artificially increased stability for warmer scenarios in low water years. Therefore, we divided the study period into smaller periods and modeled the summer (May to October) and the rest of the year (November to April) separately in low water years (Appendix A.4). For these low water years, we set the sediment temperature equal to the historical (no change run) or projected (warming scenario runs) mean of the daily maximum air temperature for the respective periods. Results for setting sediment temperature equal to the mean of the daily mean air temperature in low water years are reported to test for the impact of sediment temperature on the temperature profile of the water column.

Dam operations changed the elevation of water withdrawal from the middle of the water column in normal years to the bottom of the water column in low water years to ensure only less turbid hypolimnetic water was withdrawn for municipal water use [27]. The dam outlet structure is comprised of four downspouts (51 cm diameter cast iron pipes controlled by 51 cm gate valves) on the face of the dam for selective level draft control. Outlets are at the right side of the dam (looking downstream), not at the center line, and the overflow area is much shallower than other normal discharges. However, W2 is a two-dimensional model and does not fully reflect the three-dimensional flow dynamics. Preliminary model runs showed poor fits in the temperature profiles, particularly when the reservoir was observed to be isothermal in low water years. These poor fits led to simulated stratification when none was observed. Therefore, we raised the upper withdrawal elevation by 1.7 m to 88.3 m a.s.l and lowered the lower withdrawal elevation by 4.8 m to 72.6 m a.s.l, which improved model fit substantially, particularly in low water years. We determined that outlet elevation was the most impactful and sensitive parameter that could explain why the shapes of the temperature profile were modeled to be either isothermal (which matched observations) or stratified (which was erroneous) in low water years. The shapes of the profile could not be corrected by adjusting any other parameters (e.g., coefficients of sediment heat transfer, wind sheltering, or shading). There may be some discrepancies in the model bathymetry and other bathymetric data used to build the model due to sedimentation (Appendix A.2).

The model was further calibrated by adjusting radiation and wind coefficients until the model matched well with observed water levels and profiles of water temperature. Short-wave solar radiation was reduced by a static shade coefficient, which was calibrated to 0.94 [20], meaning that 94% of the incoming short-wave solar radiation was allowed to reach each surface segment. A shade coefficient was ideal for modeling a long, narrow reservoir, such as the study reservoir, which is in a canyon and surrounded by rolling hills. In low water years, shading increased and the value of the shading coefficient was calibrated to 80%. Radiation penetrates the surface of the water, heating the water column, and decays exponentially with depth, according to Beer’s Law:

H_s(Z) = (1 − β_S) H_s × e^−ηZ

(1)

where H_s(Z) is short-wave radiation (W m⁻²) at depth Z (m), β_S is the fraction absorbed at the water surface, H_s is short-wave radiation reaching the surface (W m⁻²), e is the base of the natural logarithm, and η is the radiation extinction coefficient in the water column (m⁻¹). The extinction coefficient, η, was estimated for each observation of Secchi disk transparency depth, Z_SD:

η = K/Z_SD

(2)

where K is a constant, 1.7 [28,29]. Secchi disk transparency depth was chronically low (mean = 1.0 ± 0.05 m) over the study period so the mean extinction coefficient value was 2.3 (± 0.1) m⁻¹, which was typical for a eutrophic system, in which values can range from 0.5 to 4.0 m⁻¹ [30]. Monthly values for the extinction coefficient were linearly interpolated to daily values for the study period. The coefficient calculated in this approach empirically accounted for total radiation extinction due to the combination of, for example, pure water, algae, suspended solids, zooplankton, and macrophytes. We did not model any biological processes using W2 to isolate the abiotic effects of warming on the study reservoir and to reduce the number of parameters for which we did not have data.

The evaporative wind speed function, f(W), in W m⁻² mm Hg⁻¹, is given as:

f(W) = a + bW^c

(3)

where a, b, and c are default empirical coefficients (9.2, 0.46, and 2, respectively) [31] and W is wind speed measured at 2 m above the ground (m s⁻¹) multiplied by the wind sheltering coefficient. Water levels and temperature profiles were insensitive to the wind sheltering coefficient and so the wind coefficient was set to 1.00, meaning that wind speed was not changed from observed values.

Model fit was determined visually and statistically. After the general shape (e.g., straight versus curved) in a temperature profile was inspected and determined to be correctly simulated, then the more precise differences between the modeled and observed temperature profiles were quantified using the absolute mean error (AME) and root mean square error (RMSE), given as:

$$\mathrm{AME}=\frac{\sum |X_{i,\mathrm{observed}}-X_{i,\mathrm{modeled}}|}{n}$$

(4)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^n{{(X}_{i,\mathrm{observed}}-X_{i,\mathrm{modeled}})}^2}{n}}$$

(5)

where X_i is the ith temperature (°C) and n is the number of observations in a profile. The AME is the average absolute deviation between modeled and observed temperature in degrees (°C) and we considered an AME of 1 °C or below to be a good fit [20]. The AME and RMSE have been used to assess the accuracy of W2 models in other studies [20].

2.5. Climate Warming Scenarios

The calibrated model was run with no change (NC) for reference, then with ∆T taken from four projections of mean air temperature for 2035 to 2064. Projections were made by the Parallel Climate Model (PCM) from the National Center for Atmospheric Research and Department of Energy, and the Geophysical Fluid Dynamics Laboratory CM2.1 model (GFDL) from the National Oceanic and Atmospheric Administration. The models were selected because they produced a more realistic representation of local climate, accounting for diverse topography and Mediterranean seasonality [32]. Projections were made for low (B1) and medium-high (A2) CO₂ emissions scenarios [33], resulting in four climate warming scenarios (PCM.B1, PCM.A2, GFDL.B1, GFDL.A2;

Table 1. Projected changes in mean air temperature in southern California for 2035 to 2064. All values are °C.

Time Period	PCM		GFDL
	B1	A2	B1	A2
	Scenario 1	Scenario 2	Scenario 3	Scenario 4
Spring/Autumn	0.8	1.2	2.1	2.3
Summer	0.8	1.3	2.3	3.1
Winter	0.6	1.0	1.6	1.7
Annual (weighted average)	0.8	1.2	2.0	2.4

Notes: All changes are positive and relative to a reference period, 1961 to 1990. B1 (low emissions) and A2 (medium-high emissions) scenarios were simulated using GFDL and PCM global climate models. Spring/Autumn refer to March to May and September to November, respectively; Summer here refers to June to August; Winter refers to December to February. Values were taken from [32].

). These warming scenarios are referred to hereafter in increasing order of warming as Scenario 1 (PCM.B1), Scenario 2 (PCM.A2), Scenario 3 (GFDL.B1), and Scenario 4 (GFDL.A2). The B1 emission scenarios assumed that CO₂ emissions will reach a maximum of ~12 GtC year⁻¹ in 2050 before dropping below ca. 2000 levels to 5.2 GtC year⁻¹ by 2100. The A2 emission scenarios assumed that CO₂ emissions will increase throughout the century to ~29 GtC year⁻¹ by 2100. Cumulative CO₂ emitted from 1990 to 2100 varied from 989 (B1) to 1773 (A2) GtC.

2.6. Data Analysis

2.6.1. Reservoir Temperature and Stability

Volume-weighted water temperature was calculated for the whole reservoir (T_R), epilimnion (T_E), and hypolimnion (T_H) for each observed and modeled water temperature profile. The thermocline was defined as the horizontal plane of the maximum temperature gradient, with a minimum of 1 °C m⁻¹, in the vertical water column [34]. The lower boundary of the epilimnion was defined as the deepest depth above the thermocline where the temperature gradient between adjacent strata was less than 0.5 °C m⁻¹ and the upper boundary of the hypolimnion was defined as the shallowest depth below the thermocline where the temperature gradient between adjacent strata was less than 0.5 °C m⁻¹. The metalimnion was defined as the water body in between the epilimnion and hypolimnion.

Stability was quantified using the Schmidt Stability Index, S (J m⁻²), which is the minimum amount of work that must be done by wind to return a reservoir to an isothermal state, calculated from differences in water density in the water column [35]:

$$\mathrm{S}=\frac{g}{A_0}\underset{0}{\overset{Z_{\mathrm{m}}}{{\displaystyle \int }}}(Z-Z^{\ast })({\rho }_Z-{\rho }^{\ast })A_Z\mathrm{dZ}$$

(6)

where g is acceleration due to gravity (m s⁻²), A₀ is surface area of the reservoir (m²), A_Z is reservoir area (m²) at depth Z (m), ρ_Z is water density (kg m⁻³) at depth Z, ρ* is mean water density (kg m⁻³) of the reservoir, and Z* is the depth (m) where the mean water density occurs. The summation is taken over all depths at approximately 1 m intervals, dZ, from the surface to the maximum depth, Z_m. Water density, ρ (kg m⁻³), was calculated as a function of temperature, T (°C), using the Thiessen-Scheel-Diesselhorst Equation [36].

Hodges was defined as thermally stratified when there was development of a thermocline and when S exceeded 100 J m⁻². Linear interpolation of S between observation dates was used to determine dates of onset of stratification and turnover. Results were aggregated to means for the calendar year (January to December) and non-stratified winter period (November to March).

2.6.2. Multiple Regression Models of Controls on Temperature and Stability

Multiple regression was used to identify the controls on changes in water temperature and stability that resulted from climate warming. Minimal adequate regression models of stability and water temperature were determined through forward and backward stepwise variable selection. Backward selection was an iterative process of deleting the least significant term in a maximal model until remaining terms were significant (α = 0.05). An appropriate model was selected based on explanatory power (R²), plausibility given physical processes, and parsimony. The relative contribution to R² from each term, alone and in combination with others, in a model was based on sequentially determined R² values, which depended on the order in which terms were entered into a regression. To control for the sequence of addition, the unique contribution to R² for an explanatory variable was given by the Lindeman, Merenda, and Gold statistic (LMG) [37,38]. The statistic can be interpreted as the average squared semi-partial correlation coefficient for a regressor, which is determined as the average of all possible permutational orderings of that regressor within the set of all regressors.

3. Results and Discussion

3.1. Model Calibration

Observed and modeled water depths (Z) under historical climate and hydrologic conditions were generally in good agreement (Figure 2). W2 correctly simulated water depth during most reservoir filling events and extended periods of drawdown and evaporation. However, an underestimation for the fill event in 1998 may have led to an underestimation in depth between 1998 and 2003.

Observed and modeled temperature profiles under historical conditions with no change were in good agreement (Figure 3, Figure 4 and Figure 5). The profile’s AME was typically below 1 °C in most months (Figure 3a). The AME was 0.84 ± 0.03 °C when averaged over the entire study period, and similar (0.78 ± 0.10 °C) for low water years (Figure 3a), confirming that W2 adequately modeled temperature for highly variable water levels. Observed and modeled epilimnetic and hypolimnetic temperatures were significantly correlated (p < 0.001), with Adj. R²’s of 0.83 and 0.42, respectively (Figure 3b,c). Monthly temperature profiles for a typical year (1995) with normal water levels are presented in Figure 4 to show the model fit to observed profiles and forecast warming of water temperature due to four projected scenarios (Table 1). When the sediment temperature was set equal to the mean of daily mean air temperature in the low water years, the shapes of the temperature profiles were erroneous, particularly in the (observed) isothermal hypolimnion in September and October, and mean AME was higher (0.95 ± 0.14 °C), suggesting that the mean of daily maximum (instead of mean) air temperature was a better assumption for sediment temperature in low water years. Temperature profiles in the month of October (1992 to 2011) are presented in Figure 5 to show the consistency of the model in simulating stratification over the study period and also to show the robustness of the model in simulating the contrasting temperature profiles in low water years, when the reservoir did not stratify and the sediment temperature was set to the mean of daily maximum air temperature.

3.2. Higher Rates of Evaporation under Warming Scenarios

Observed and modeled evaporation (E) under historical conditions were significantly correlated (Adj. R² = 0.71; p < 0.001; Figure 6a) and mean annual averages for the study period were similar (E_OBS = 154.0 cm year⁻¹; E_MOD = 159.0 cm year⁻¹). Evaporation increased by 2.6 to 7.9% for all warming scenarios relative to the reference run (E_Scenario1 = 163.1; E_Scenario2 = 165.5; E_Scenario3 = 170.6; E_Scenario4 = 171.6; all units cm year⁻¹), which was comparable to the average increase (3.6%) projected across 47 reservoirs in California under similar climate forcing [39]. Though the impact of these increases in evaporation is small for deeper water bodies, it can be substantial in shallower reservoirs, including the study reservoir, during years with low water levels ($\overline{Z}$ < 4.5 m).

Increased evaporation due to warming reduced modeled $\overline{Z}$ by 0.04 to 0.09 m (

Table 2. Mean annual (January to December) or winter (November to March; denoted by a subscript “W”) reservoir characteristics simulated over the study period (1992 to 2011), excluding low water years (2002 to 2004). Means calculated for only low water years (2002 to 2004) are denoted by a subscript “NS”, indicating that there was no stratification in those years. Means include reservoir maximum (Z and Z_NS) and mean ($\overline{Z}$ and $\overline{Z_{\mathrm{NS}}}$) depth; evaporation (E and E_NS); thicknesses of the epilimnion (L_E) and hypolimnion (L_H); depth of the thermocline (Z_T); temperatures of the whole reservoir (T_R, T_R-W, T_R-NS), epilimnion (T_E), and hypolimnion (T_H); Schmidt stability (S, S_W, and S_NS); and onset (STR_ON) and duration (STR_DUR) of stratification. Values for warming scenarios (Scenarios 1 to 4) are differences relative to the reference run with no change.

	Observed	No change	Scenario 1	Scenario 2	Scenario 3	Scenario 4
$\overline{Z}$ (m)	7.5	7.5	−0.04	−0.05	−0.09	−0.09
$\overline{Z_{\mathrm{NS}}}$ (m)	4.5	4.4	−0.08	−0.20	−0.25	−0.24
E (cm)	151.7	157.2	3.5	5.2	9.0	10.2
ENS (cm)	166.9	169.3	7.7	13.6	26.1	25.9
LE (m)	4.2	3.1	0.0	0.0	−0.1	−0.1
LH (m)	14.5	13.5	−0.1	−0.2	−0.4	−0.4
ZT (m)	6.5	6.0	−0.1	−0.1	−0.1	−0.1
TR (°C)	17.6	17.9	0.5	0.8	1.3	1.5
TR-W (°C)	14.4	14.6	0.5	0.8	1.3	1.5
TR-NS (°C)	18.7	18.4	0.6	0.9	1.6	1.8
TE (°C)	23.1	24.4	0.4	0.7	1.2	1.5
TH (°C)	14.1	14.2	0.4	0.6	1.0	1.1
S (J m−2)	173.7	168.4	6.1	9.3	17.4	23.6
SW (J m−2)	36.7	28.3	2.6	3.8	7.2	8.0
SNS (J m−2)	13.3	14.2	0.1	−0.1	−0.1	1.1
STRON (day of year)	85	91	−4	−5	−7	−9
STRDUR (day)	220	209	6	8	12	14

), averaged over the study period, after removing data from low water years. Modeled reductions in depth under a warmer climate were ~3 times larger during low water years due to higher evaporation rates, continued withdrawals for urban use, and the rapid rate of change in depth for a given change in volume at low water depths. Higher rates of E during low water years were plausible because water temperature, net radiation, and the vapor deficit were all higher in those years. The temperature of the whole reservoir was observed and modeled to be ~1 °C higher in low water years than in normal years (Table 2). The bathymetry of the reservoir is such that depth changed slowly for a given volume change at high water depths but rapidly at low water depths, so withdrawals for urban use had a larger impact on depth when Z was low (Figure 2).

3.3. A Warmer Thermal Regime

There was a strong coupling of air and water temperatures throughout the profile. Annual weighted-average warming of air temperature (T_A) by 0.8 to 2.4 °C increased water temperature by 0.5 to 1.5 °C, on average, in the whole reservoir (Table 2) and epilimnetic and hypolimnetic temperatures increased significantly with increasing air temperature across scenarios (Figure 7). Increases in water temperatures were larger for warmer scenarios. Regressions showed that a 1 °C increase in summer air temperature corresponded with an increase of 0.44 °C in the whole reservoir (T_R = 0.18 + 0.44 × T_A; Adj. R² = 0.94; p < 0.05; Figure 7), 0.47 °C in the epilimnion (T_E = 0.09 + 0.47 × T_A; Adj. R² = 0.99; p < 0.01), and 0.30 °C in the hypolimnion (T_H = 0.19 + 0.30 × T_A; Adj. R² = 0.89; p < 0.05). Both whole-reservoir temperature averaged over the year (T_R) and over winter (T_R-W) increased more in warmer scenarios (Table 2). Contrary to initial expectations, increases in T_R-W were similar to increases in T_R in all warming scenarios.

Changes in T_E explained the majority of variance in the changes in T_H in a warmer climate (LMG = 0.66) and changes in hypolimnetic thickness (L_H) explained little variance (LMG = 0.03;

Table 3. Regression results for mean annual (January to December) temperature of the hypolimnion (T_H) and Schmidt stability (S). Independent variables include temperatures of the epilimnion (T_E), hypolimnion (T_H), and air (T_A), thickness of the hypolimnion (L_H), and mean reservoir depth ($\overline{Z}$). All values are differences between a warming scenario relative to the reference run with no change. The β is the regression parameter value and LMG is the relative contribution of each variable to the final R² [37,38]. Data from low water years (2002 to 2004) were excluded given that the study reservoir did not stratify then.

Response Variable	β			LMG			Adj. R2	RMSE
	TE	LH		TE	LH
TH	0.75 ***	0.27 ***		0.66	0.03		0.68	0.26
	TA	$\overline{Z}$	TH	TA	$\overline{Z}$	TH
S	17.45 ***	65.20 ***	−11.22 ***	0.67	0.09	0.11	0.86	3.20

* p < 0.05; *** p < 0.001.

), which showed the importance of the epilimnion in mediating warming of the hypolimnion. In a relatively shallow reservoir, heating of the hypolimnion is not directly related to meteorological forcing by itself [5] but is rather a complex interaction between energy exchanges between the epilimnion and hypolimnion, which is mediated by mixing and stratification. Furthermore, heating of the hypolimnion is due primarily to warming from the atmosphere and downward transfer of heat through convection from the epilimnion [40], rather than to secondary sources of heating, including precipitation, inflow, and benthic sediment. The total input of heat to the reservoir increased under warming, but heat at the surface increased resistance to mixing, reducing the downward movement of heat across the thermocline [8].

3.4. Stability and Stratification

Mean annual stability (S) increased by 6.1 to 23.6 (J m⁻²) in normal years, by less in the winter (S_W; 2.6 to 8.0 J m⁻²), and by near zero in low water years (S_NS; −0.1 to 1.1 J m⁻²), all with generally larger increases for warmer scenarios (Table 2). In normal years, the duration of stratification increased for all warming scenarios by 6 to 14 days, with the onset of stratification moving earlier in the year by 4 to 9 days on average over the study period (Table 2). W2 correctly simulated the low stability and observed inhibition of stratification in low water years in the reference run. Stratification did not occur during these years in any of the warming scenarios, either.

Increases in S were dampened both in the winter and in low water years because the water column was isothermal. Stability increased less in low water years than in the winter of normal years because of a warmer and shallower water column. In low water years, water temperature was initially high due to preferential removal of colder hypolimnetic water through a lower elevation outlet, rather than due to heating from the air. Also, there was likely a larger influence from warmer benthic sediment on warming of the hypolimnion when the water level was low. This suggested that when the water level was low, air temperature, which we assumed impacted sediment temperature, had opposing impacts on S depending on the relative changes at the top and bottom of the temperature profile.

To understand the impact of climate warming on stability in normal years, we did separate correlation analyses between relevant reservoir variables. In Figure 8, circles are correlations between metrics under the changed climate in a management-dominated system where water level fluctuations were the primary control on hydrodynamics. In the same figure, squares are correlations between changes in metrics after climate warming relative to the reference run (with no change). The squares in Figure 8 were calculated as the correlations of the differences in the values of a given variable between the given warming scenario and the reference run. This analysis held water level fluctuation and management constant and showed only the impact of increased air temperature (and related variables) on reservoir metrics.

In the management-dominated system, S consistently correlated with $\overline{Z}$ and L_H (Figure 8 circles), but in the climate-modified system, S did not consistently correlate with any single variable (Figure 8 squares). Multiple regression results showed that changes in T_A, $\overline{Z}$, and T_H together explained the change in S (Adj. R² = 0.86; Table 3). This suggested that increased T_A (LMG = 0.67) and E, which also directly impacted $\overline{Z}$ (LMG = 0.09), were important and significant predictors of changes in S. Air temperature increased the slope of the depth-stability relationship, and mean annual stability was higher for all depths in warmer scenarios, with increased stability for warmer scenarios (Figure 9).

Climate warming changed the distribution of heat in the water column. Changes in depth and heat distribution in the profile (epilimnetic thickness L_E, L_H, and thermocline depth Z_T) were highly correlated with each other (Figure 8 squares). The metalimnion expanded as L_E, L_H, and Z_T decreased on average in all scenarios over the study period, with larger decreases in L_E, L_H, and Z_T in warmer scenarios (Table 2). Thermocline depth has increased in many lakes in response to climate warming [2]. In contrast, in the study reservoir Z_T decreased likely due to factors including lower water depth due to increased E, and high turbidity, which limited penetration of light and heat, as observed in other reservoirs [8,41]. Future research could model the impact of allochthonous nutrient supply and climate warming on algal production and the resulting direction of change in Z_T.

3.5. Implications of Changes to the Hydrodynamic Regime

Higher air temperatures increased water temperature throughout the profile and throughout the year (Figure 4). The change in temperature has implications for other water quality parameters. Increased water temperature decreases DO saturation. At water temperatures around 14 °C (mean observed T_H; Table 2), the DO saturation decreases by about 0.20 mg L⁻¹ per one degree of increase in water temperature. Hence, the DO saturation would decrease by approximately 0.06 mg L⁻¹ for the expected 0.30 °C increase in hypolimnetic water temperature, per one degree in the air temperature (Figure 7). Warmer water has been shown to correlate with increased sediment microbial metabolism [42], hypolimnetic anoxia [5], and release of phosphorus from sediment [43].

Warmer climate had variable impacts on the hydrodynamic regime of the reservoir, depending on its depth. Warmer climate increased stability and extended the duration of stratification in normal years. The combination of increased water temperatures and stability can favor nuisance cyanobacteria, increasing their ambient biomass levels and the frequency and persistence of blooms [41]. However, hydrodynamic conditions were different in low water years, due primarily to management (i.e., volume, timing, and elevation of water withdrawal) removing colder hypolimnetic water. Under low depth conditions, the reservoir became isothermal, essentially behaving as a well-mixed epilimnetic layer, and the magnitude of increase in stability from climate warming was near zero. The reservoir became sufficiently shallow that solar radiation penetrated deeper through the water body than when water level was higher, which led to higher water temperatures and evaporation rates compared with normal years. Water temperature and evaporation rates may be more impacted by climate warming due to additional processes of sediment warming in reservoirs that are actively managed and have shallower depths.

Sediment can heat or cool the hypolimnion seasonally [40] so the timing of equilibration of sediment temperature to air temperature has important effects on stability and mixing of nutrients. Equilibration appeared to be rapid in low water years, but, in early model runs when assuming instant equilibration in normal years (sediment temperature = air temperature), average warming of the hypolimnion was higher than of the epilimnion (causing smaller increases of stability), which suggested that sediment was a net heat source in the model. This seemed physically implausible because sediment does not produce heat and sediment heat exchange with water is small compared to surface heat exchange [20]. Therefore, we reported results from when sediment temperature was set equal to the historical air temperature record, recognizing that we likely underestimated sediment temperature in warming scenarios, which could result in an underestimate of warming of the hypolimnion and overestimate of stability. These findings demonstrate the importance of sediment-water energy exchange, particularly in shallower reservoirs, and future work in such systems should simulate dynamic sediment temperature and its interactions with the water column [24].

In low water years with low stability, it was possible that release of phosphorus from sediment may still have occurred, despite the lack of long-term stratification, which can result in enhanced mixing and increased delivery of nutrients to the euphotic zone. It is known that in a polymictic reservoir, short periods of stratification, on the order of hours or days, may occur [44] and cause sediment release of phosphorus if there are conditions of low DO and high rates of oxygen reduction at the sediment-water interface [45]. The study reservoir showed perennially low hypolimnetic DO concentrations during the study period [11], and previous studies report low DO concentrations after mixis in December [46] and heterograde DO profiles when the reservoir was isothermal in January [47], indicating that the hypolimnion was rich in nutrients and could support production through winter. When depth and stability are low, algal production can then be supported by vertical transport of hypolimnetic and sedimentary nutrient reserves released from the sediment [48] and exposed littoral zone [49,50]. Sediment can be in constant contact with the euphotic layer and phosphorus concentrations can remain high because of the lack of water to dilute those concentrations in that layer.

This study provides an example of a managed system in which air and water temperatures were tightly coupled, so that abiotic climate warming significantly exacerbated changes in water level and, in turn, stability. Inhibition of stratification due to forecast warmer climate is projected to turn lakes from a dimictic to monomictic regime in other settings [51], though these previous results are not for managed water bodies with high water level fluctuations. Projected future increases in water demand could severely impact reservoir storage levels in San Diego [52] and management decisions in response will strongly control depth and stability in such managed reservoirs [2,11]. Management decisions should consider the effects of warmer climate on reservoir dynamics and potential for mixing of water and phosphorus as it relates to water quality.

4. Summary and Conclusions

Warm monomictic systems have global importance given that they predominate in the tropics and lower temperate latitudes. Our model showed that climate warming could increase water temperature, stability, and duration of stratification, while reducing water depth and thicknesses of the epilimnion and hypolimnion, in years with “normal” water depths (mean depth > 4.5 m). In contrast, in years of low water depths (mean depth < 4.5 m), the reservoir was projected to respond differently to warming, due in part to the thermal dynamics of sediment, so that there were relatively larger increases in warming throughout the water column, and a relatively smaller increase in stability. The combined effects of climate warming and reservoir management on changes in water levels can significantly impact the hydrodynamic regime, with implications for nutrient cycling and algal production.