The competing impacts of climate change and nutrient reductions on dissolved oxygen in Chesapeake Bay

The Chesapeake Bay region is projected to experience changes in temperature, sea level, and precipitation as a result of climate change. This research uses an estuarine-watershed hydrodynamic–biogeochemical modeling system along with projected mid-21st-century changes in temperature, freshwater flow, and sea level rise to explore the impact climate change may have on future Chesapeake Bay dissolved-oxygen (DO) concentrations and the potential success of nutrient reductions in attaining mandated estuarine water quality improvements. Results indicate that warming bay waters will decrease oxygen solubility year-round, while also increasing oxygen utilization via respiration and remineralization, primarily impacting bottom oxygen in the spring. Rising sea level will increase estuarine circulation, reducing residence time in bottom waters and increasing stratification. As a result, oxygen concentrations in bottom waters are projected to increase, while oxygen concentrations at mid-depths (3 < DO < 5 mg L⁻¹) will typically decrease. Changes in precipitation are projected to deliver higher winter and spring freshwater flow and nutrient loads, fueling increased primary production. Together, these multiple climate impacts will lower DO throughout the Chesapeake Bay and negatively impact progress towards meeting water quality standards associated with the Chesapeake Bay Total Maximum Daily Load. However, this research also shows that the potential impacts of climate change will be significantly smaller than improvements in DO expected in response to the required nutrient reductions, especially at the anoxic and hypoxic levels. Overall, increased temperature exhibits the strongest control on the change in future DO concentrations, primarily due to decreased solubility, while sea level rise is expected to exert a small positive impact and increased winter river flow is anticipated to exert a small negative impact.

1 Introduction

Global climate change is projected to alter the world's marine environments with coastal and estuarine systems bearing exacerbated impacts. Rising temperatures and sea levels, along with changes in precipitation patterns, have the potential to dramatically alter water quality conditions in these highly productive and increasingly human-influenced systems (Najjar et al., 2010; Altieri and Gedan, 2015). While there are multiple metrics with which to evaluate water quality, dissolved-oxygen (DO) concentrations are widely used to identify systems in distress. Large volumes of hypoxic water (generally considered to be waters with DO < 2 mg L⁻¹), commonly referred to as dead zones, can be found in many coastal oceans and estuaries around the world (Diaz and Rosenberg, 2008). As the climate continues to change, it is important to evaluate the impact these changes will have on DO concentrations in critical coastal environments like the Chesapeake Bay.

Climate change is generally predicted to have a net negative effect on DO in coastal waters through changes in temperature, sea level, and precipitation (Boesch et al., 2007; Meier et al., 2011; Altieri and Gedan, 2015). Higher temperatures impact both the timing and rates of biological functions, while also potentially driving changes in oxygen production and consumption (Winder and Sommer, 2012). Although increased temperature is not anticipated to have a major effect on estuarine stratification, which is primarily controlled by salinity in systems such as the Chesapeake Bay (Murphy et al., 2011), the increased temperature will act to reduce the amount of oxygen a given volume of water can hold via decreased solubility. Sea level rise (SLR) can act to increase estuarine circulation, water column stratification, and water body volume (Chua and Xu, 2014). These impacts are counteractive, as increasing volume and circulation can bring in high-oxygen water from the coastal ocean, while increased stratification inhibits downward mixing of the high-oxygen water from the surface waters. Stronger estuarine circulation generally also leads to shorter residence times that typically increase oxygen concentrations (Hong and Shen, 2012; Du and Shen, 2015). In addition, over much of the mid-Atlantic region, annual precipitation, and thus river discharge, has been increasing (Tian et al., 2015; Yang et al., 2015a, b). In the future, precipitation is most likely to increase most during the winter/spring and in the northern part of the region (Najjar et al., 2009; IPCC Annex I, 2013), delivering higher river flows and nutrient loads that fuel spring productivity and produce more organic matter available for summer decomposition (Najjar et al., 2010). Changes in nutrient loading and hydrologic conditions can also alter the bay's phytoplankton composition, changing the biomass available for eventual decomposition (Harding et al., 2016).

Compounding the complicated process of projecting future water quality conditions are nutrient management efforts such as the Chesapeake Bay 2010 Total Maximum Daily Load (TMDL; USEPA, 2010) that was developed to improve water quality conditions in the bay by decreasing nutrient and sediment loads. These nutrient management efforts should be fully implemented by 2025 with the ultimate goal of reducing summer hypoxia (Keisman and Shenk, 2013). Examining the potential impact of climate change in light of these mandated nutrient reductions is important because the multiple impacts of climate change have the potential to render current nutrient reduction goals inadequate (Justic et al., 2007; Meier et al., 2013; Lennartz et al., 2014; Altieri and Gedan, 2015). Furthermore, assessing the science behind climate change impacts is critical for policies like the Chesapeake Bay TMDL that are prone to legal challenges (McCormick et al., 2017).

While much of the discussion around water quality regulations focuses on hypoxia (DO < 2 mg L⁻¹), studying low-DO water that encompasses concentrations greater than hypoxic levels (DO concentrations up to 5 mg L⁻¹) is also critical due to the impact of increases in temperature on economically important fisheries. For example, not only do temperature increases impact DO concentrations, they also increase metabolic rates in fish. This increase causes fish to experience adverse health impacts at higher and higher DO concentrations (Portner and Knust, 2007; Vaquer-Sunyer and Duarte, 2011; Buchheister et al., 2013). Furthermore, the TMDL mandates multiple levels of minimum DO concentrations at various times and locations throughout the Chesapeake Bay (USEPA, 2010; Tango and Batiuk, 2013). While much of the regulation targets traditional hypoxia, the TMDL mandates a monthly mean DO ≥ 3 mg L⁻¹ in the deep water of the bay to protect the survival and recruitment of bay anchovy eggs and larvae and a monthly mean of DO ≥ 5 mg L⁻¹ above the pycnocline to protect the growth of larval, juvenile, and adult fish and shellfish (Tango and Batiuk, 2013).

This study examines the impact of climate change on oxygen concentrations in the Chesapeake Bay by utilizing a coupled hydrodynamic–biogeochemical model that has previously been compared to other Chesapeake Bay models (Irby et al., 2016). As the Chesapeake Bay TMDL stipulates a time horizon of 2025 for full nutrient reduction implementation, this research assumes that the required nutrient management strategies are in place and limiting nutrient delivery to their full potential. Future estimates of mid-21st-century temperature, SLR, and watershed nutrient loading are applied to the model in order to examine the sensitivity of the individual and combined impacts of these environmental changes on future anoxic (< 0.2 mg L⁻¹), hypoxic (< 2 mg L⁻¹), and low-DO (2–5 mg L⁻¹) water in the Chesapeake Bay.

2 Methods

2.1 ChesROMS-ECB

The estuarine model is based on the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005) and uses the Chesapeake Bay curvilinear horizontal grid (ChesROMS) of Xu et al. (2012) with an average wet cell resolution inside the bay of 1.7 km. As in Feng et al. (2015), the model is configured to use the recursive MPDATA 3-D advection scheme for tracers, third-order upstream advection scheme for horizontal momentum, and fourth-order centered difference for momentum in the vertical, with a 20-layer vertically stretched sigma grid. The Estuarine-Carbon-Biogeochemistry (ECB) component of the model (Feng et al., 2015) was developed originally from a continental shelf application (Hofmann et al., 2011), using dissolved organic matter cycling similar to that described in Druon et al. (2010). With only single phytoplankton and zooplankton classes and only one limiting nutrient (nitrogen), the ECB model is simpler than that employed by the Chesapeake Bay Program (Cerco et al., 2010) but is more complex than simple dissolved-oxygen models that utilize a constant oxygen consumption rate (e.g., Scully, 2010; Bever et al., 2013). ChesROMS-ECB has been previously shown to adequately resolve the spatial and temporal variability of key physical and biological variables such as temperature, salinity, nitrogen, and DO (Feng et al., 2015; Irby et al., 2016).

Before using ChesROMS-ECB to determine the impact of changes in temperature on water quality parameters, the temperature dependence of the biogeochemical formulations within the model required a careful evaluation. Several biogeochemical formulations within ChesROMS-ECB did not previously include a dependence on temperature, and temperature dependence was added as part of this study (a complete list of model changes is provided in Appendix A). For example, temperature dependence was introduced to the rates for maximum phytoplankton growth, zooplankton grazing/growth, nitrification, detrital solubilization, and detrital remineralization. All modifications introduce an exponential relationship between temperature and maximum rate, except for maximum phytoplankton growth. The function for phytoplankton growth is based on Lomas et al. (2002) and employs a constant growth rate below 20 ^∘C of 2.15 d⁻¹, with an exponential maximum growth curve for temperatures above 20 ^∘C. Remineralization of the dissolved organic constituents previously included temperature dependence, but to ensure consistency, these rates were modified to match the Chesapeake-specific community respiration Q₁₀ values from Lomas et al. (2002).

An additional two changes were made to improve the light attenuation parameterization in ChesROMS-ECB. First, a minimum value of 0.6 m⁻¹ was applied to the diffuse attenuation coefficient, based on model–data comparisons (Wang et al., 2009; Son and Wang, 2015). Second, the organic portion of the total suspended solids term in the light attenuation formulation of Feng et al. (2015) was multiplied by 2, since carbon is generally considered to be roughly half of the total weight of organic matter.

Table 1 Definitions of sensitivity experiments.

To assess the relative skill of the revised model, the skill in reproducing water quality observations at 23 stations along the bay (Fig. 1, Table A1 in the Appendix) was compared to the skill of the earlier version of the model used in Feng et al. (2015) and Irby et al. (2016). The 23 stations were assigned to four regions that are functionally delineated by salinity characteristics, with Region A representing the oligohaline, Regions B and C representing the upper and lower mesohaline (and generally the lowest DO concentrations), and Region D representing the polyhaline. The updated model retained its gross skill in terms of total root mean squared difference (Table A2) compared to the version of the model evaluated in Irby et al. (2016). Specifically, the updated model improved bottom DO skill in Regions C and D, primarily due to the light attenuation modifications mentioned above (see Appendix A for details).

Figure 1 Map of the Chesapeake Bay with water quality monitoring stations (Table A1) identified by region, based primarily on salinity. A: oligohaline; B and C: upper and lower mesohaline; D: polyhaline.

Figure 2 (a) Total (bay-wide) riverine nitrogen loading into the bay, (b) time series (7-day running mean) of bottom DO with and without nutrient reductions (TMDL vs. Base) and with and without climate change (allCC vs. noCC) for the average of the stations in the upper mesohaline Region B, and (c) as in (b) but for the lower mesohaline Region C.

2.2 Chesapeake Bay Program Watershed Model

This study utilizes freshwater discharge and riverine nitrogen and sediment concentrations from the Chesapeake Bay Program's Watershed Model (version 5.3.2) that was used in the development of the 2010 TMDL (Shenk and Linker, 2013). (As in Feng et al., 2015, riverine carbon concentrations that are required as inputs to ChesROMS-ECB were obtained from the Dynamic Land Ecosystem Model; Tian et al., 2015.) This research generally assumes that the management practices required to meet the 2010 TMDL nutrient reductions in the absence of climate change (Shenk and Linker, 2013) are fully realized; however, a brief examination of the potential impact of climate change without nutrient reduction is also explored. Because the TMDL is based on a reference time period of 1993–1995 (USEPA, 2010), these are the years used in this study. Fortuitously, this period includes both relatively wet years (1993, 1994) and a dry year (1995), allowing the investigation of how future climate change impacts are affected by natural interannual variability. Simulations using the TMDL reduction in nutrient concentrations are hereafter referred to as the TMDL scenarios while the base 1993 to 1995 simulations will hereafter be referred to as the Base run (Table 1). Compared to the Base run, the TMDL scenarios include a bay-wide reduction in riverine nutrient loading of 45, 44, and 38 % for the three years (1993 to 1995), respectively (Fig. 2a).

2.3 Climate change sensitivity experiments

In this study, the sensitivity of Chesapeake Bay hypoxia to projected regional impacts for three aspects of climate change (temperature, SLR, and precipitation/rivers) is examined. A time horizon of the mid-21st century is chosen for these changes because it is far enough in the future to allow for the assumption that the TMDL nutrient reductions have been fully implemented (including nutrient transport lag effects), while also being soon enough for relatively constrained estimates of future climate change impacts.

2.3.1 Temperature

By 2050, the Chesapeake Bay region is expected to experience air temperature increases greater than the global average. Specifically, the IPCC projection of a median annual average atmospheric temperature increase for 2046-2065 relative to 1986–2005 for the Chesapeake Bay region is about 2 ^∘C (∼ 0.036 ^∘C yr⁻¹; IPCC Annex I, 2013), whereas the analogous global increase is projected to be 1.4 ^∘C (∼ 0.025 ^∘C yr⁻¹; IPCC Summary, 2013). Further research from the IPCC establishes that ocean warming tends to be 20 to 40 % lower than the rate of atmospheric warming (Collins et al., 2013). As the Chesapeake Bay is a relatively shallow, well-mixed estuary and there has recently been an observed increase in the rate of Chesapeake Bay warming (Ding and Elmore, 2015), this research utilizes a ratio between atmospheric and ocean warming that is slightly lower than the open ocean range. The 1.75 ^∘C (∼ 0.032 ^∘C yr⁻¹) increase in bay water temperature for the mid-21st century relative to the mid-1990s used in this study (Table 1) is higher than the ∼ 0.02 ^∘C yr⁻¹ observed Chesapeake Bay warming between 1949 and 2002 (Preston, 2004). However, Preston (2004) found evidence of increased warming in the late 1990s. The rate of warming used in this analysis is also consistent with projected increases from downscaled global climate models for the bay (Muhling et al., 2017). It is also slightly lower than the warming estimated using a high-resolution climate model (CM2.6) for the location of the ChesROMS open boundary (2.6 ^∘C; Saba et al., 2016) and less than the average satellite-derived rate of bay surface water warming of 0.005–0.175 ^∘C yr⁻¹ from 1984 to 2007 (Ding and Elmore, 2015).

The 1.75 ^∘C water temperature increase was applied uniformly across time and space to biogeochemical processes and oxygen solubility throughout the bay, but the temperature increase was not applied to other physical properties or processes, such as water density gradients or meteorological forcing. Thus, these increased temperatures do not impact stratification or other physical dynamics of the bay within the model. This approach implicitly assumes that the bay is shallow enough that climatic warming will occur uniformly over time. Supporting this assumption, Preston (2004) found that the surface and subsurface waters of the bay warmed at relatively similar rates, even finding that, on average, the subsurface waters warmed slightly faster than surface waters. In addition, recent trends in the intensification of early summer stratification have been found not to be due to water column temperature changes but rather are primarily due to changes in salinity as a result of SLR and altered freshwater inputs (Murphy et al., 2011). Changes in salinity along the ChesROMS open boundary on the continental shelf between the 1990s and the mid-20th century have been computed by Saba et al. (2016) to be very minor (∼ 0.2 psu) and are thus not considered here. This is consistent with our goal of examining the first-order impacts of temperature change on hypoxia (through solubility and growth/grazing/remineralization changes); the effect of warming on Chesapeake Bay hydrodynamics is being analyzed in a separate follow-on study. The temperature increase sensitivity experiment will hereafter be referred to as the TMDL + tempCC scenario since the increase in temperature is applied to the TMDL nutrient scenario (Table 1).

2.3.2 Sea level rise (SLR)

The Chesapeake Bay is expected to incur a greater increase in sea level than the global average, and the bay has experienced a recent acceleration in SLR, as has most of the mid-Atlantic coast (Sallenger et al., 2012). Boon and Mitchell (2015) found a roughly 0.1 m increase in sea level in Norfolk, Virginia between 1993 and 2014. Assuming a linear extrapolation of that rate (∼ 5 mm yr⁻¹), by the mid-21st century, Norfolk would expect an SLR of 0.3 m relative to the mid-1990s. However, the linear extrapolation ignores the projected, and recently observed, acceleration in SLR. Incorporating anticipated acceleration, Boon and Mitchell (2015) estimate an average increase in SLR by 2050 of ∼ 0.5 m in the Chesapeake Bay relative to the relative mean sea level between 1969 and 2014. This estimate is similar to that of Sweet et al. (2017), who, using downscaled global models, estimated a similar SLR in the mid-Atlantic for 2050 under an intermediate emissions scenario. Based on these recent regional estimates, this research assumes a mid-21st-century SLR of 0.5m (∼ 9 mm yr⁻¹) relative to the mid-1990s.

Table 2 Monthly freshwater discharge fractional change factor used for the TMDL + riverCC, TMDL + allCC, and Base + allCC scenarios, calculated as the ratio between the freshwater inputs in 2050 divided by the freshwater inputs in the Base run.

Model implementation of SLR follows that of Hong and Shen (2012). The 0.5 m increase was added to the free water surface layer at the outer boundary of the model grid, along the continental shelf. The vertical grid stretching parameters were not altered and the simulation required less than 6 months for the bay to equilibrate to the increased sea level. The SLR sensitivity experiment will hereafter be referred to as the TMDL + slrCC scenario since the 0.5m increase is applied to the TMDL scenario (Table 1).

2.3.3 River flow

The Chesapeake Bay watershed is expected to undergo a range of precipitation changes over the next century, with the southern portion of the watershed expected to experience a lower intensity change than the northern portion of the watershed, complicating estimates of future precipitation change and, as a result, river flow (Najjar et al., 2009). The future river loading estimates used here (Fig. 2a) were provided by the Chesapeake Bay Program and were derived from an implementation of the Chesapeake Bay Watershed Model (Sect. 2.2) that incorporated downscaled (1∕8^∘ resolution) precipitation and temperature estimates for the mid-21st century from multiple global climate model realizations (Hargreaves and Samani, 1982; Groisman et al., 2004; Reclamation, 2013). Overall, increases in precipitation over the Chesapeake Bay watershed resulted in generally greater runoff to the Chesapeake Bay, especially in the winter and spring months, even though the warmer temperatures throughout the year mitigated some of these increases via increased rates of evapotranspiration.

From these Watershed Model results, the ratio of monthly freshwater delivery to the bay for the climate change scenario relative to the Base run was calculated for the Susquehanna River (Table 2) and was applied to all rivers in ChesROMS-ECB. This is a reasonable approach given that the Susquehanna watershed accounts for > 80 % of the bay watershed area that drains directly to the main stem and is the primary source of the nutrients that drive the summer hypoxic region of the bay between the Patapsco River in the north and the Rappahannock River in the south (Hagy et al., 2004). Overall, the resulting increase in river flow applied to the model (Table 2) causes both an increase in freshwater discharge and an increase in nutrient delivery (Fig. 2a). The combined impact of increased freshwater flow and nutrient loads will hereafter be referred to as the TMDL + riverCC scenario (Table 1).

2.3.4 Combined climate change sensitivity experiment

A final sensitivity experiment that combines all three of the climate change impacts was run for both nutrient cases, i.e., the TMDL scenario (reduced nutrients) and the Base run (realistic nutrients). These experiments will hereafter be referred to as the TMDL + allCC and Base + allCC scenarios, respectively, since the combined impact of all climate change variables (temperature, SLR, and rivers) was applied (Table 1).

2.4 Dissolved-oxygen analysis

Hypoxic volume is a commonly used metric to quantify the amount of water that experiences a given level of DO concentration over a specific time (e.g., Murphy et al., 2011; Bever et al., 2013). In this study, two metrics related to hypoxic volume are computed: cumulative hypoxic volume (CHV) and hypoxic duration (HD). CHV is calculated as the sum of each day's hypoxic volume over a year (Bever et al., 2013), and HD is computed as the number of days that have a hypoxic volume greater than 1 km³. While traditional DO concentration levels of hypoxia (< 2 mg L⁻¹) and anoxia (< 0.2 mg L⁻¹) are examined, this research also considers impacts of low-DO, defined here as DO < 5 mg L⁻¹. This level is consistent with the highest DO concentrations stipulated in the Chesapeake Bay water quality standards (USEPA, 2010) and is a conservative upper bound on DO concentration found to initiate stress on marine fish (Vaquer-Sunyer and Duarte, 2008; Buchheister et al., 2013).

3 Results

The impact of nutrient reduction on bottom DO concentrations is greater than that of climate change (Fig. 2b, c). The reduction in nutrients (between the Base run and TMDL scenario) causes a general increase in bottom DO concentrations. This impact is largest during the drawdown of bottom oxygen in the spring (April–June), dampens during the course of the summer, and is lowest in winter (December–February). In Region B, the region of the bay where oxygen concentrations are lowest and most persistent, this impact is strongest in the driest year (1995), during which the increase in bottom DO exceeds 2.5 mg L⁻¹. In 1993 and 1994 the bottom DO increase is only around 1.5 mg L⁻¹ (Fig. 2). In contrast Region C, encompassing the southern extent of the hypoxic zone, experiences a greater increase in spring bottom DO than Region B in the wet years (> 2 mg L⁻¹ in 1993 and 1994) and a smaller increase in the dry year (∼ 1.5 mg L⁻¹ in 1995).

Climate change has a smaller effect on bottom DO concentrations than the TMDL nutrient reductions. Climate change has almost no impact on bottom DO during the peak of summer when bottom DO concentrations are the lowest (near zero). In the Base run (realistic nutrient inputs), the effect of climate change on spring bottom DO is a decrease of ∼ 0.6 and ∼ 0.8 mg L⁻¹ in Regions B and C, respectively. Climate change impacts bottom DO similarly in the TMDL scenario, with reductions in spring bottom DO of ∼ 0.5 mg L⁻¹ in both Regions B and C (Fig. 2). In both regions, these reductions in bottom DO are similar in all three years.

Table 3 Average change in bottom DO (mg L⁻¹) relative to the TMDL + noCC scenario for each experiment and region.

Table 4 Average (1993–1995) change in cumulative hypoxic volume (km³ days) for various oxygen concentration ranges, relative to the TMDL + noCC experiment.

Figure 3 Time series (7-day running mean) of the change in bottom DO between the TMDL climate change and no climate change scenarios for the average of the (a) upper mesohaline Region B and (b) lower mesohaline Region C.

Of the three climate factors considered (temperature, SLR, and river flow), temperature had the largest impact on bottom DO. As a result, the TMDL + allCC scenario is most similar to the TMDL + tempCC scenario (Fig. 3). In Region B, the TMDL + slrCC and the TMDL + riverCC scenarios have a smaller impact on bottom DO during the wet years of 1993 and 1994 than during the dry year of 1995. The opposite occurs in Region C, with the TMDL + slrCC and the TMDL + riverCC scenarios having a larger impact on bottom DO during the wet years of 1993 and 1994 than during the dry year of 1995. In both regions, the impact of SLR generally increases bottom DO during the spring and summer, while changes in the rivers (increased seasonality and nutrient load) suppress DO. These two essentially equal and opposite effects largely cancel each other out (Fig. 3).

Although temperature had the largest impact on bottom DO in each of the four regions considered, the magnitude of the individual impacts of climate change differed by region (Table 3) and by oxygen concentration (Table 4). Specifically, in the TMDL + allCC scenario, bottom DO decreased compared to the TDML + noCC run in all four regions, with Region A exhibiting the highest total average change (−0.58 mg L⁻¹) and the other three regions all exhibiting an average change of roughly −0.4 mg L⁻¹ (Table 3). This is primarily due to the large decreases in bottom DO in the TMDL + slrCC scenario in Region A (−0.21 mg L⁻¹), relative to the small (mostly positive) impacts due to sea level rise in the other regions. Overall, the impact of all three of the climate change factors is nearly linearly additive (Table 3).

Figure 4 Cumulative hypoxic volume for six ranges of DO concentrations for each of the study years and each of the scenarios (colored circles).

The relative impact of the three climate change factors also varied with oxygen concentration, particularly for temperature, which exerted the greatest impact on cumulative hypoxic volume at 3 < DO < 5 mg L⁻¹ (Table 4). The TMDL + slrCC scenario increased oxygen at low concentrations (DO < 3 mg L⁻¹) and decreased oxygen at higher concentrations (3 < DO < 5 mg L⁻¹). In contrast, the impact of changes in riverine nutrient loading was relatively similar for all oxygen concentrations, typically generating a small decrease in DO (Table 4).

The CHV for all of the TMDL scenarios (both with and without climate change) is less than the CHV from the Base run without climate change (Fig. 4). This pattern holds true for all six DO levels examined (< 0.2 to < 5 mg L⁻¹). At each DO level, the CHV for the dry year (1995) is much less than for the wet years (1993 and 1994) for each TMDL scenario. Furthermore, the CHV for the TMDL scenarios in the wet years is generally higher than the CHV from the Base run for the dry year (Fig. 4). The CHV in the TMDL + slrCC and TMDL + riverCC scenarios tends to track closely to the TMDL + noCC scenario, while the TMDL + tempCC scenario is most similar to the TMDL + allCC scenario (Fig. 4), as was also the case for bottom DO (Fig. 3).

Figure 5 Percent change due to climate change, relative to the improvement in CHV between the TMDL + noCC scenario and Base + noCC run. Percent change in CHV gain is defined as (TMDL + xx – TMDL + noCC)/(TMDL + noCC – Base run + noCC).

Table 5 Percent^* of 3-year average bottom DO change as a result of the temperature experiment due to solubility for each region at the surface and bottom of the water column.

^* Percent calculated as the expected change in bottom DO, as predicted by solubility, divided by the modeled change in bottom DO.

The percent change in CHV relative to the progress, or gains, made in CHV by applying the TMDL nutrient reductions varies across DO level and by scenario (Fig. 5). In general, the TMDL + slrCC scenario resulted in a ∼ 0–10 % increase in the improvement made by the TMDL scenario (here, an increase in gains is actually a decrease in CHV) across all DO levels and all years. In contrast, the TMDL + riverCC and TMDL + tempCC scenarios resulted in a degradation of the system, compared to the TMDL + noCC scenario. The TMDL + riverCC scenario consistently causes a loss of ∼ 0–5 % of the gains, with slightly larger losses in 1994 and 1995 at higher DO levels. The TMDL + riverCC scenario combines two separate, but linked, climate change impacts: increased freshwater flow (particularly in the winter) and increased nutrient loads (as a result of increased freshwater flow). While not shown, separate experiments isolating the impacts of flow and load demonstrated that the increase in nutrient load, rather than the increase in freshwater flow, caused the increase in CHV in the TMDL + riverCC scenario. The TMDL + tempCC scenario was the strongest function of DO level, with a relatively small loss of ∼ 5 % at the < 0.2 mg L⁻¹ level and a large ∼ 40 % loss at the < 5 mg L⁻¹ level (Fig. 5). The combined effect of climate change (TMDL + allCC) was a net increase in CHV of more than 50 % over the TMDL + noCC scenario in the wet years of 1993 and 1994 for DO < 5mg L⁻¹ and a corresponding 40 % increase in CHV for the dry year of 1995 (Fig. 5).

Figure 6 DO differences due to climate change (between the TMDL + noCC and TMDL + tempCC scenarios) averaged for (a) the stations in Region B and (b) the stations in Region C. The black lines are the average change expected if only solubility was impacted by an increase in temperature. The red lines are the modeled change in DO as a result of the increase in temperature affecting both solubility and biological oxygen production/demand.

As shown above, increased temperature generally maintains the greatest control on the TMDL + allCC scenario (Fig. 4). The impact of temperature on DO in this analysis is due to two factors: chemical solubility and biological oxygen demand. (The impact of temperature on DO is not due to a change in stratification, since the experiment was explicitly designed to focus on the impacts of chemical solubility and biological oxygen demand and neglect any change in stratification, which previous studies (Preston, 2004) have suggested is small.) To further isolate these two impacts, the differences in modeled DO computed with and without warming are computed considering only solubility effects and considering both solubility and biological oxygen demand (Fig. 6; Table 5). Since oxygen saturation is more sensitive to changes in temperature at low temperatures, there is a larger change in DO as a result of changes in solubility during the winter than during the summer, even though the change in temperature is constant in time. Deviations from the change in DO due to solubility can be attributed to changes in biological oxygen demand and can be estimated by comparing the simulation assuming only solubility impacts (red line in Fig. 6) with the simulation assuming temperature changes affect both solubility and biological oxygen demand (black line in Fig. 6). Overall, 65–85 % of the change in DO between the TMDL + tempCC scenario compared to the TMDL + noCC scenario is a result of temperature's impact on solubility, with solubility exerting a larger impact at the surface than at depth (Table 5). The impact of biological oxygen demand is consistently negative at depth during the spring and early summer, enhancing the initiation of hypoxic conditions (Fig. 6b).

Figure 7 Bars showing duration of hypoxic volume (> 1 km³) at each DO level for the Base + noCC run and the TMDL + noCC and TMDL + allCC nutrient scenarios.

In terms of the number of days that the bay experiences hypoxic and low-oxygen conditions each year, climate change reduces the positive impact of the nutrient reduction (Fig. 7). While there is a moderate decrease in hypoxic duration resulting from the nutrient reduction, the TMDL + allCC scenario demonstrates that when climate change is included all levels of low-DO and hypoxia initiate an average of ∼ 7 days earlier. This trend is not evident in the cessation of hypoxia and low DO, i.e., climate change does not necessarily cause hypoxia to last later in the year. While all three years exhibit a similar pattern and timeline of cessation of low-DO with < 0.2 mg L⁻¹ ceasing 3–4 months before < 5 mg L⁻¹, each year is different in terms of initiation timing. In 1993 for the Base + noCC run, all levels of DO initiate within 2 weeks of each other. This timing holds true for the TMDL scenarios as well but with anoxia lagging behind. In 1994 in the Base + noCC run, there is a steady progression from low DO to anoxia over ∼ 6 weeks. In the TMDL scenarios, that is extended to ∼ 3 months. In 1995, the TMDL nutrient reduction results in no DO < 1 mg L⁻¹ and significantly delays the onset of low-DO by up to ∼ 3 months compared to the Base run.

Figure 8 Latitudinal along-bay DO transects for the Base + noCC run and TMDL + noCC scenario for 1 July 1993 (a, c) and 1 July 1995 (b, d) and average summer (May–September) for 1993 (e, g) and 1995 (f, h); the difference in average summer DO between the TMDL + noCC and TMDL + allCC scenarios (i, j) and the percent difference in average summer DO between the TMDL + noCC and TMDL + allCC scenarios (k, l).

Nutrient reduction primarily reduces the horizontal extent of the hypoxic zone (Fig. 8). Examining a south–north transect along the main stem of the bay for 1 July 1993 (Fig. 8a, c) and 1995 (Fig. 8b, d) reveals that nutrient reduction acts to compress the southern extent of the hypoxic zone more than the northern extent. One similarity between all four subplots (a–d) is the vertical extent of the low-oxygen waters, which are capped by the pycnocline at ∼ 5 m depth. As expected, the extent and severity of anoxia and hypoxia on 1 July is much greater than the summer (May–September) average for both the Base + noCC run and TMDL + noCC scenario (Fig. 8e–h). In general, the impact of climate change is greater in the dry year (1995; Fig. 8j, l) than in the wet year (1993; Fig. 8i, k). The location of the greatest magnitude change is near the pycnocline depth (Fig. 8i, j), but the location of greatest percent change is below the pycnocline (Fig. 8k, l).

Figure 9 Percent of the entire bay that experiences a given DO level during 1993, 1994, and 1995.

The climate change sensitivity experiments cause a larger volume of the bay to experience low-DO concentrations in both wet and dry years and under both the Base + allCC and TMDL + allCC scenarios (Fig. 9). While climate change does not greatly exacerbate the volume of the bay that experiences anoxic and hypoxic conditions, climate change increases the percent of the bay experiencing conditions of DO < 5 mg L⁻¹ by up to ∼ 6 %, regardless of whether or not the TMDL nutrient reductions have occurred. Similarly, regardless of whether or not climate change occurs, the volume of the bay experiencing low-DO under nutrient reduction is ∼ 10 % lower than that in the 1993–1995 Base run nutrient conditions. Overall, the dry year (1995) results in ∼ 30–50 % as much of the bay experiencing low DO and hypoxic waters as compared to the wet years (1993, 1994).

4 Discussion

4.1 How will Chesapeake Bay DO concentrations change in the future as a result of climate change?

By the mid-21st century, low-DO conditions can be expected to begin about 1 week earlier due to climate change, with increases in volume and extent being largest at the margins and at the southern extent of the hypoxic zone. Significant impacts will be felt on water with DO concentrations in the range of 2–5 mg L⁻¹ and not only on hypoxic waters (DO < 2 mg L⁻¹).

The most consistent impact across all levels of low-DO waters due to climate change is an earlier onset of hypoxic and low-DO conditions by an average of ∼ 7 days. While an earlier onset was identified, there was no trend in the cessation of hypoxic and low-DO conditions, with climate change sometimes causing an earlier and sometimes a later cessation. Furthermore, an earlier onset of conditions is projected to occur under both nutrient-reduced and nutrient-replete futures. The pattern of earlier onset is primarily due to the additive impacts of an increase in spring biological oxygen utilization at depth and decreased solubility, both the result of the increase in temperature (Fig. 6). An analysis of climate change impacts on DO of an estuarine tributary of the Chesapeake Bay similarly found that hypoxic duration is likely to be extended in the future (Lake and Brush, 2015).

In terms of a change in the volume of low-DO waters, the relative impact of climate change increases with DO concentration (Figs. 4, 5). The improvements due to the nutrient reductions are reduced by climate change, ranging from ∼ 5 % for DO < 0.2 to ∼ 45 % for DO < 5 mg L⁻¹. The difference between impacts at anoxic levels versus waters with DO of 3–5 mg L⁻¹ is accentuated during the dry year of 1995 due to the fact that the nutrient reductions result in no modeled DO < 1 mg L⁻¹ during this year (Fig. 7), regardless of whether or not climate change is occurring. Even assuming realistic 1995 nutrient inputs, the volume and duration of anoxia under climate change in 1995 is very small.

Throughout the water column, the greatest change in DO will be at the edges of the low-DO and hypoxic zones, particularly at the southern and vertical extents (Fig. 8). Conversely, the smallest changes will occur in the anoxic waters where DO cannot be decreased further (Fig. 8). As hypoxia is capped by the pycnocline (Irby et al., 2016), the magnitude of DO change (∼ 0.5 mg L⁻¹) is not great enough to extend low-DO conditions to the DO-replete surface waters. Laterally, the largest changes in bottom DO will be in the southern extent of hypoxia and the degree of east–west compression along the main stem of the bay. Such a change would be likely to detrimentally impact demersal fish and shellfish communities along the shallow flanks of the bay and its tributaries.

4.2 How will the individual impacts of climate change (increased temperatures, SLR, river flow) affect DO concentrations in the Chesapeake Bay?

The combined impacts of climate change will cause reduced DO concentrations in the mid-21st century, with increased water temperatures being the strongest driver of this change.

In examining the individual impacts of future (mid-21st-century) temperature, SLR, and river flow on Chesapeake Bay DO concentrations, temperature exhibits the largest overall impact (Figs. 4, 5; Tables 3, 4). The decrease in DO associated with increases in temperature is also consistent with other modeling research focused on the York River estuary, a tributary of the Chesapeake Bay (Lake and Brush, 2015). The present research demonstrates the importance of temperature on solubility, as the annual average impact of temperature on oxygen saturation outpaced the impact of temperature on biological functions on average by roughly 2 : 1 in the region of the bay that experiences hypoxia (Fig. 6; Table 5). This ratio is decreased to roughly 1 : 1 during the spring/early summer drawdown of bottom DO in the main stem channel (Fig. 6). Murphy et al. (2011) similarly found that increased respiration due to increased temperature potentially plays a smaller role in changes in hypoxia than the physical and chemical changes. However, it is possible that as temperature continues to increase, the ratio of impact between solubility and biological oxygen demand may shift toward a greater influence by biological oxygen demand. This is because the additional impact of further changes in solubility will decrease as temperatures continue to rise, while biological respiration at depth and production at the surface may continue to steadily rise with increasingly warmer temperatures.

Increasing sea levels can impact future Chesapeake Bay oxygen concentrations in multiple ways. By increasing estuarine circulation and decreasing residence time (Hong and Shen, 2012; Du and Shen, 2015), rising sea levels can actually increase bottom oxygen concentrations in the most anoxic portions of the deep main stem trench. At the same time, increasing estuarine circulation increases stratification (Chua and Xu, 2014), which serves to further decrease oxygen concentrations (Lennartz et al., 2014). In the Chesapeake Bay simulations presented here, the former process dominates at depth where SLR results in higher oxygen concentrations (Fig. 3; Table 4). By contrast, the latter process (increased stratification) dominates in the mid-water column, where SLR results in lower oxygen concentrations (Table 4). In the driest year (1995) the overall impact of SLR is a decrease in the total volume of water with DO < 5 mg L⁻¹ (Fig. 5), whereas in the wettest year (1994) when stratification is already relatively high, the overall impact of SLR is dominated by the stratification effect and results in a net increase in the total volume of water with DO < 5 mg L⁻¹. The larger impact of SLR during dry years is consistent with a study from the Delaware Bay showing that high flow dampens the salinity intrusion that results from SLR (Ross et al., 2015) and with a study in San Francisco Bay finding that the impact of SLR is limited under high-flow conditions (Chua and Xu, 2014).

Future climate change will also modify freshwater and nutrient loading from the watershed to the bay, causing the largest increases in the winter months (December–February; Table 2). This increase in river nutrient loading to the bay (Fig. 2a) reduces oxygen concentrations both at the bottom (Fig. 3) as well as throughout the water column, resulting in greater cumulative hypoxic volumes (Figs. 4, 5; Table 4). However, this impact is much smaller than that of temperature (Table 4). The increase in winter precipitation, delivering both increased freshwater flow and increased nutrient loads, accounts for a larger percentage of the overall change in DO during the dry year (1995) because the low-flow conditions cause the bay to be more sensitive to changes in freshwater flow and nutrient loading.

4.3 How might climate change impact the success of the 2010 TMDL nutrient reductions?

Climate change may cause the 2010 TMDL nutrient reductions to be insufficient to meet the required water quality improvements in the Chesapeake Bay.

This research demonstrates that the improvements in Chesapeake Bay water quality due to the TMDL nutrient reductions are much greater than the deleterious impacts of ∼ 2050 climate change; however, results also indicate that by the mid-21st century, climate change will likely decrease oxygen levels and increase both hypoxic volume and hypoxic duration. Because some locations in the bay barely pass TMDL standards and others require special allowances to meet the standards (Irby and Friedrichs, 2018), even these small increases in anoxic and hypoxic conditions can cause locations that previously passed the water quality standards to fail under a changing climate. The DO minima in the TMDL regulations are based on both space and time criteria. Although the spatial dimension may not be greatly impacted at the anoxic and hypoxic levels, this research suggests that the temporal dimension will be. This could cause locations in the bay that are currently projected to pass the minimum standards to fail them in light of climate change, simply due to an extension of the hypoxic season without an expansion of hypoxic volume.

With increased temperature being the primary cause of the impact of climate change on DO concentrations, it is important to consider other potential impacts increased temperature may have on the ecosystem in the context of the success of the TMDL nutrient reductions. Temperature increases in the Chesapeake Bay are anticipated to produce temperatures outside of previously observed extremes (Muhling et al., 2017), lending increased pertinence to understanding the impact of temperature changes on meeting water quality goals. In light of this, the impact on the TMDL of a decrease in oxygen concentrations due to climate change should be viewed in conjunction with the impact increased temperature is likely to have on the species upon which the DO levels in the TMDL nutrient reductions were predicated. Multiple studies have established that increasing water temperature increases metabolic rates in fish that cause them to experience negative health impacts at higher DO concentrations than they do at lower temperatures (Breitburg, 2002; Portner and Lanning, 2009; Lapointe et al., 2014). Due to those compounding impacts and the large role temperature is expected to play in regulating future DO, it may be prudent for the TMDL to elevate the mandated minimum DO levels in an effort to protect target species. If this occurred, the impacts of climate change would likely cause an even larger failure rate of TMDL standards.

4.4 How will climate change impact DO if the TMDL nutrient reductions are not met?

Although the relative impact of climate change is similar on a reduced-nutrient future and a high-nutrient future, the degree of interannual variability in hypoxia may change in a reduced- versus high-nutrient future due to differences in the responses of oxygen to fluctuations between dry and wet years.

The relative impact of climate change on a reduced-nutrient versus a high-nutrient future is similar in terms of hypoxic volume and duration. In both a low- and high-nutrient future, the percent of the bay that experiences a given DO level is increased with climate change (Fig. 9). Furthermore, in both cases, the impact of climate change at low-DO concentrations (< 5 mg L⁻¹) is greater than that at hypoxic levels (< 2 mg L⁻¹). In terms of relative change in DO along the main stem of the bay, a high-nutrient future is expected to experience a higher (∼ 9–15 %) change in DO concentration than a low-nutrient future (∼ 6–9 %), with the largest changes in both cases occurring at the southern end of the hypoxic zone (Fig. 8).

The largest potential ecological difference between the two futures is in the dry year of 1995. In this year TMDL scenarios exhibited no anoxia in the bay, regardless of whether or not climate change was occurring. This suggests that during dry years, when the nutrient reduction may be sufficient to alleviate anoxia, climate change impacts may not be large enough to overcome the hysteric or threshold level of DO initiation similar to what has been observed with hypoxic responses to nutrient loading (Kemp et al., 2009). It may seem counterintuitive, but this suggests that the interannual variability of anoxic conditions may be exacerbated in a future with nutrient reduction because the interannual percent change in anoxic conditions will be relative to ∼ 0 % in the very dry years. Because of this, when climate change is added to the TMDL nutrient reductions, there is likely to be greater disparity in terms of anoxic volume between wet and dry years. Further intensifying the difference between wet and dry years is the potential impact of nutrient storage in the watershed during dry years that is delivered to the bay in a successive wet year, amplifying hypoxia and anoxia (Lee et al., 2016).

4.5 Methodological limitations

This research is a first-order assessment of the potential impacts that changes in climate may have on the efficacy of nutrient reduction efforts in the Chesapeake Bay; however, more robust examinations of the problem are needed in order to adequately aid in the regulatory decision making process going forward. As the present research has identified increased temperature as the largest contributor to changes in DO, future efforts should work to incorporate the impact of increased air temperature and changes in meteorological forcing on the air–sea interface and bay hydrodynamics. In addition, increased stream temperatures will likely need to be accounted for, as there is evidence that the current rates of bay warming cannot be fully explained by the observed increase in regional air temperatures (Ding and Elmore, 2015). Estimates of future precipitation indicate changes in storm intensification and extreme events that could have dramatic effects on nutrient delivery to the bay (Sinha et al., 2017), and thus these should be considered in future work as well. Finally, the atmospheric wind field will likely change in the future. Although there is significant uncertainty associated with future wind projections, the strong impact of wind on hypoxia in the Chesapeake Bay (Scully, 2010) makes this an important issue to better understand.

Due to the uncertainty in projected changes in temperature, river flow, and SLR, assessing the sensitivity of DO to multiple estimates of climate change will be important. This research establishes that the increase in temperature has the strongest control on DO, but that does not mean that DO concentrations are most sensitive to the bounds of potential 2050 temperature changes. While the high computational expense of running multiple sensitivity tests through complex coupled hydrodynamic–biogeochemical models can be prohibitive, establishing a range of uncertainty is critical to informed adaptive-management decision making.

Additional limitations are related to timing. For example, the present research assumes a discontinuity between the reduction in nutrients and the changes in climate. This is an unrealistic assumption because the nutrient reductions and climate change will continue to occur contemporaneously. These changes are also not immediate but manifest over time in a continuously evolving environment. In addition, the current approach simply identifies the potential ramifications of climate change on nutrient reduction efforts but does not establish a timeline for the water quality changes as a result of nutrient reductions to occur. This means that climate change has the potential to further limit the effectiveness of nutrient reduction efforts because the impacts of climate change may be more immediate than the impacts of nutrient reduction. To address these limitations, an effort to conduct a continuous simulation from 2015 to 2050 including both gradual implementation of the nutrient reductions and climate change impacts is currently underway.

5 Conclusions

Overall, the most striking result of this research is that the potential impact of climate change by the mid-21st century is much smaller than the impact of the 2010 TMDL nutrient reductions, particularly at anoxic and hypoxic levels. However, the decrease in DO concentrations resulting from the combined impacts of climate change may cause portions of the bay that currently meet mandated water quality standards to fail them in the future. At the most stringent DO standards, this is primarily due to an increase in hypoxic duration rather than hypoxic volume, as under climate change, the onset of hypoxic conditions is projected to initiate ∼ 7 days earlier on average across all DO concentrations 0.2–5 mg L⁻¹.

Changes in DO as a result of the increase in temperature dominate the combined climate change impact. While the influence of solubility on DO concentrations is the primary control on decreased DO throughout the year, the impact of increased biological oxygen demand is most prevalent at the bottom in the spring to early summer, contributing to the earlier initiation of hypoxic conditions. The impact of temperature is likely to affect low-oxygen tolerance of higher trophic levels as well by increasing metabolic rates, making species less tolerant at higher DO levels. This may result in the DO minimums mandated in the water quality standards being insufficient to protect key species even if the current goals are met.

Both sea level rise and changes in river flow exert a greater influence on change in DO during dry, low-streamflow years. Changes in river flow are likely to deliver higher freshwater flows during the winter and spring that will both deliver higher nutrient loads and increase estuarine circulation. These two effects impact DO concentrations oppositely, with higher loads resulting in more organic matter being available for decomposition and increased estuarine circulation delivering more oxygen-rich ocean water; however, the impact of increased loads outcompetes the greater circulation. Sea level rise exerts the only net positive impact of climate change on DO concentrations, increasing the effectiveness of the TMDL nutrient reductions by ∼ 5 % in the mesohaline. However, this positive impact is undermined by the large negative impact of temperature.

The relative effects of climate change are similar whether the TMDL nutrient reductions are achieved or not. In both cases, there is a slight increase in anoxic conditions, and the relative impact of climate change intensifies at higher DO concentrations (3–5 mg L⁻¹). The impact of the nutrient reductions on dry years is accentuated compared to the “business as usual” dry years due to the greater moderating influence sea level rise exerts during low-flow conditions. This results in anoxic and hypoxic conditions being depressed with nutrient reduction plus climate change in the dry year of 1995, but not when climate change is combined with no nutrient reduction.

Overall, this study demonstrates that climate change has the potential to limit the effectiveness of future management actions aimed at reducing nutrient inputs to the Chesapeake Bay. However, the negative impacts of climate change are smaller than the positive impacts resulting from the mandated nutrient reductions. Given that this analysis only considers a mid-21st-century time horizon and climate impacts are expected to intensify with time, it is critical to continue to examine how the bay may evolve in the future.

Appendix A:  

Before being used for climate change sensitivity experiments, the ChesROMS-ECB temperature parameterizations were reexamined and modified as necessary based on information from the literature and extensive skill assessment using data from 23 Chesapeake Bay Program Water Quality Monitoring stations (Table A1). (Data are available at http://www.chesapeakebay.net/data/downloads/cbpwaterqualitydatabase1984present, last access: September 2017.) Modifications of biological functions from the model version published in Feng et al. (2015) and used in the model comparisons published in Irby et al. (2016) are documented in Table A2. Specifically, temperature dependence was added to the zooplankton maximum growth rate, the remineralization rates of large and small detritus, and the phytoplankton growth rate at temperatures above 20 ^∘C. The maximum rate of nitrification, the temperature dependency on remineralization of semi-labile DON, and the remineralization rate of DOC at 0 ^∘C were also modified to fit with current understanding (Lomas et al., 2002).

The skill of the modified model was assessed via total root mean squared difference (RMSD; Table A3), normalized target diagrams (Joliff et al., 2009), and time series analysis (Irby, 2017). For the total RMSD calculations, the model results were compared to monthly/bimonthly observations at the stations and regions shown in Fig. 1. Results from the modified model were also compared to an earlier iteration of the model evaluated in Irby et al. (2016).

Table A1 Characteristics of observation stations.

Table A2 Modifications to formulations and parameter values from Feng et al. (2015).

^a Community respiration and zooplankton grazing temperature-dependent functions are based on a Q₁₀ of 2.1 (Lomas et al., 2002).
^b Phytoplankton growth rate at low temperatures (T < 20 ^∘C) is constant with higher temperatures following a rate based on Lomas et al. (2002) with an average Q₁₀ of between 20 ^∘C and 40 ^∘C of ∼ 2.4.

Table A3 Total RMSD (and observational mean) of surface and bottom temperature (T), salinity (S), dissolved oxygen (DO), and nitrate (NO₃) of the present model and the earlier model version used in Feng et al. (2015) and Irby et al. (2016) for the four regions (A, B, C, D) defined in Fig. 1.