Long-term trends, current status, and transitions of water quality in Chesapeake Bay

Coincident climatic and human effects strongly influence water-quality properties in estuarine-coastal ecosystems around the world. Time-series data for a number of ecosystems reveal high spatio-temporal variability superimposed on secular trends traceable to nutrient over-enrichment. In this paper, we present new analyses of long-term data for Chesapeake Bay directed at several goals: (1) to distinguish trends from spatio-temporal variability imposed by climatic effects; (2) to assess long-term trends of water-quality properties reflecting degradation and recovery; (3) to propose numerical water-quality criteria as targets for restoration; (4) to assess progress toward attainment of these targets. The bay has experienced multiple impairments associated with nutrient over-enrichment since World War II, e.g., low dissolved oxygen (DO), decreased water clarity, and harmful algal blooms (HAB). Anthropogenic eutrophication has been expressed as increased chlorophyll-a (chl-a) driven by accelerated nutrient loading from 1945 to 1980. Management intervention led to decreased loading thereafter, but deleterious symptoms of excess nutrients persist. Climatic effects exemplified by irregular “dry” and “wet” periods in the last 30+ years largely explain high inter-annual variability of water-quality properties, requiring adjustments to resolve long-term trends. Here, we extend these analyses at a finer temporal scale to six decades of chl-a, Secchi depth, and nitrite plus nitrate (NO2 + NO3) data to support trend analyses and the development of numerical water-quality criteria. The proposed criteria build on a conceptual model emphasizing the need to distinguish climatic and human effects in gauging progress to reverse eutrophication in estuarine-coastal ecosystems.

Modeling climatic effects. Model predictions of log 10 chl-a, Secchi depth, and NO 2 + NO 3 in low-flow, mean-flow, and high-flow conditions were developed as time series by applying GAM to long-term data to adjust for climatic effects (Fig. 4a-i). These climatic effects were computed by varying input terms for SRF and salinity for each salinity zone. Model predictions in mean-flow conditions were obtained by setting flow and salinity terms at their respective means, "dry" conditions with flow and salinity terms set at 10 th and 90 th percentiles, and "wet" conditions with flow and salinity terms set at 90 th and 10 th percentiles. Model inputs also included terms for annual TN and NO 2 + NO 3 loading to account for inter-annual variability of N-limitation 17,18 . Temporal lags were tested using auto-regression terms (AR) in generalized additive mixed models (GAMM), revealing no significant differences from GAM. Additional details on models are provided in Methods.  18 , modified to cover the time frame for data analyzed in this paper.
Model predictions of mean, monthly chl-a for the OH salinity zone in mean-flow conditions showed higher chl-a in the 1960s and 1970s than in recent years, with lower chl-a after the mid-to late-1970s that continued from 1985 to 2015 (Fig. 4a). Climatic effects on chl-a for the OH salinity zone were expressed as model predictions in low-flow, "dry" conditions that exceeded those in mean-flow or high-flow, "wet" conditions. This stimulatory effect on chl-a for the OH salinity zone reflected an alleviation of light-limitation as inputs of bio-optically active materials (e.g., suspended particulate matter, chromophoric dissolved organic material) were reduced in low-flow, "dry" conditions. Model predictions of mean, monthly chl-a for the OH salinity zone in mean-flow conditions showed a significant downward trend after the mid-1960s. Previous analyses based on long-term data aggregated at an annual scale also showed this step-down of chl-a for the OH salinity zone, ascribed to regulation of  www.nature.com/scientificreports www.nature.com/scientificreports/ orthophosphate (PO 4 3− ), increased P-limitation in the upper bay, and reduction of N uptake in the OH salinity zone leading to increased throughput of TN and NO 2 + NO 3 to MH and PH salinity zones 17,18 .
Model predictions of mean, monthly chl-a for the MH salinity zone in mean-flow conditions were highly variable, characterized by a shallow, increasing trend that continued throughout the time series (Fig. 4b). Flow-adjusted model predictions for the MH salinity zone showed higher chl-a in high-flow, "wet" conditions than in mean-flow or low-flow, "dry" conditions. This climatic effect for the MH salinity zone was opposite that for the OH salinity zone, a pattern consistent with increased TN and NO 2 + NO 3 loading in "wet" years, an alleviation of N-limitation, and stimulation of chl-a and NPP described previously 17,18 .
Fewer observations were available to support model predictions of mean, monthly chl-a for the PH salinity zone than for OH and MH salinity zones, although temporal variability of chl-a was comparable among the three salinity zones (Fig. 4c). Climatic effects consisted of higher chl-a in high-flow, "wet" conditions than in mean-flow or low-flow, "dry" conditions, similar to the MH salinity zone. Model predictions of mean, monthly chl-a for the PH salinity zone in mean-flow conditions did not indicate a secular trend for the period of record (Fig. 4c).  Table 3. Statistics for generalized additive models (GAM) of water-quality properties in Chesapeake Bay using the set of predictor variables compiled in Table 1. a GCV = generalized cross-validation score [57][58][59] . www.nature.com/scientificreports www.nature.com/scientificreports/ Climatic effects on Secchi depth were expressed as higher model predictions (=increased water clarity) in low-flow, "dry" conditions than in mean-flow or high-flow, "wet" conditions for OH, MH, and PH salinity zones ( Fig. 4d-f). Model predictions of Secchi depth in mean-flow conditions showed decreasing trends (=decreased water clarity) for OH and MH salinity zones ( Fig. 4d-e), but no trend for the PH salinity zone (Fig. 4f). Model predictions of Secchi depth in high-flow, "wet" conditions were generally lower (=decreased water clarity) than predictions in mean-flow or low-flow, "dry" conditions, but overall, Secchi depth was less sensitive to climatic effects than chl-a.
Model predictions of NO 2 + NO 3 showed strong climatic effects for OH, MH, and PH salinity zones ( Fig. 4gi). NO 2 + NO 3 was consistently higher in high-flow, "wet" conditions than in mean-flow or low-flow, "dry" conditions for all salinity zones. Model predictions of NO 2 + NO 3 in mean-flow conditions showed mixed trends, consisting of a shallow increase for the OH salinity zone, and decreases for MH and PH salinity zones. These trends were consistent with decreased N loading since the early 1980s 17,18 .
Climatic effects on chl-a, Secchi depth, and NO 2 + NO 3 derived as model predictions in low-flow, "dry" conditions, mean-flow conditions, and high-flow, "wet" conditions were aggregated at a monthly scale (Fig. 5a-i), revealing patterns consistent with predictions for the complete time series (Fig. 4a-i). Mean, monthly chl-a for the OH salinity zone was higher in low-flow, "dry" conditions than in mean-flow or high-flow, "wet" conditions ( Fig. 5a), and higher for MH and PH salinity zones in high-flow, "wet" conditions ( Fig. 5b,c). Climatic effects on other water-quality properties consisted of higher Secchi depth in low-flow, "dry" conditions ( Fig. 5d-f), and higher NO 2 + NO 3 in high-flow, "wet" conditions for OH, MH, and PH salinity zones ( Fig. 5g-i).
Numerical water-quality criteria. Model predictions aggregated at a monthly scale were used to develop numerical water-quality criteria adjusted for climatic effects. We averaged predictions for spring (March-May), summer (July-Aug), and fall (Sep-Nov) to obtain seasonal values, with predictions in mean-flow conditions shown at the top of each panel (Fig. 5a-i). Numerical criteria for chl-a, Secchi depth, and NO 2 + NO 3 based on predictions for specific time periods were tabulated as entries for each water-quality property (Tables 4-6  www.nature.com/scientificreports www.nature.com/scientificreports/ criteria, revealing significant differences among properties with respect to attaining putative targets (Tables 4-6). Spring chl-a for the OH salinity zone in recent years exceeded proposed criteria based on all data, the 1960s (bold type), or previous analyses, with percent differences from 19.6 to 86.6% (Table 4). Conversely, summer and fall chl-a for the OH salinity zone in recent years were lower than proposed criteria on several time bases, and thereby    www.nature.com/scientificreports www.nature.com/scientificreports/ in compliance, with percent differences from −4.97 to −26.7% (Table 4). Seasonal chl-a for the MH salinity zone in recent years exceeded proposed criteria based on all data, the 1960s (except fall), or previous analyses, especially in spring with percent differences from 32.7 to 346% (Table 4). Lastly, seasonal chl-a concentrations for the PH salinity zone in recent years were similar to proposed criteria based on all data, but higher than criteria based on the 1960s (except fall) or previous analyses (Table 4). These findings suggest further reductions of chl-a in spring will be required to meet proposed criteria for the OH salinity zone, and in all seasons for MH and PH salinity zones, depending on the time bases used to set criteria.
Model predictions of Secchi depth for OH, MH, and PH salinity zones in mean-flow conditions were generally lower (=decreased water clarity) in recent years (2011 to 2015) than proposed criteria based on all data, the 1960s (bold type), or 1985 to 1989 (Table 5). An exception was Secchi depth for the OH salinity zone in spring that was slightly higher in recent years (=increased water clarity), ostensibly associated with a downward trajectory of chl-a in the upper bay following a ban of PO 4 3− in detergents. Model predictions for the MH salinity zone in mean-flow conditions showed lower Secchi depth (=decreased water clarity) in recent years than proposed criteria based on all data, the 1960s (bold type), or 1985 to 1989, with percent differences from 5.25 to 30.5%. Consistent with the MH salinity zone, model predictions of Secchi depth for the PH salinity zone in mean-flow conditions were generally lower (=decreased water clarity) in recent years than proposed criteria based on all data, the 1960s (bold type), or 1985 to 1989, with percent differences from 3.77 to 22.0%.  Table 6. Numerical criteria based on flow-adjusted model predictions of NO 2 + NO 3 (units -μM) in Chesapeake Bay aggregated at a monthly scale to account for climatic effects on water-quality properties.    www.nature.com/scientificreports www.nature.com/scientificreports/ Lastly, model predictions of NO 2 + NO 3 for the OH salinity zone in mean-flow conditions in recent years (2011 to 2015) consistently exceeded proposed criteria based on all data, the 1960s, or low-flow conditions in recent years (bold type) (Table 6). Conversely, model predictions of NO 2 + NO 3 for MH and PH salinity zones in low-flow conditions in recent years met proposed criteria based on all data or the 1960s. These findings were consistent with decreased TN and NO 2 + NO 3 loading after 1980 [15][16][17][18] , and increased nutrient consumption accompanying a historical increase of chl-a in the bay 9,17,18 . Model predictions of NO 2 + NO 3 for MH and PH salinity zones in low-flow conditions in recent years were consistently lower than those in mean-flow conditions, capturing climatic effects on TN and NO 2 + NO 3 loading, shown as positive percent differences (Table 6). These model predictions guided numerical criteria as they accounted for the downward trend of annual TN and NO 2 + NO 3 loading since 1980, and for climatic effects on TN and NO 2 + NO 3 loading. Basing proposed criteria for NO 2 + NO 3 on model predictions in low-flow conditions in recent years proved effective for the OH salinity zone, but seasonal depletion of NO 2 + NO 3 in summer and fall limited the usefulness of this approach for MH and PH salinity zones. A continuing increase of model predictions of chl-a for the MH salinity zone in mean-flow conditions may explain the decrease of NO 2 + NO 3 , i.e., increased consumption by phytoplankton, suggesting that progress should not be defined by a single water-quality property.

Water-quality property
Long-term trends of water-quality properties. Model predictions of water-quality properties in mean-flow conditions were used to quantify trends of chl-a, Secchi depth, and NO 2 + NO 3 as percent changes (Fig. 6a-c). Long-term trends of chl-a differed by salinity zone and season. Spring chl-a showed increasing trends for OH and MH salinity zones, but was essentially constant for the PH salinity zone; summer chl-a showed an increasing trend for the MH salinity zone, but decreasing trends for OH and PH salinity zones; fall chl-a showed consistent, decreasing trends for OH, MH, and PH salinity zones (Fig. 6a). Percent changes from 1985 to 2015 placed at the top of individual bars showed that recent trends of chl-a sometimes differed from the complete time series (compare Fig. 4-c). Examples include chl-a for the OH salinity zone in summer and fall with decreasing trends from 1965 to 2015, but increasing trends from 1985 to 2015; the MH salinity zone in fall with a decreasing trend from 1965 to 2015, but an increasing trend from 1985 to 2015; the PH salinity zone in spring with a decreasing trend from 1985 to 2015, but nearly constant chl-a for the complete time series. Secchi depth showed consistent, decreasing trends for OH, MH, and PH salinity zones from 1967 to 2015 as percent changes, with the exception of the OH salinity zone that showed an increasing trend in spring (Fig. 6b). Percent changes of Secchi depth from 1985 to 2015 placed at the top of individual bars showed decreasing trends from −13.3 to −30.1% for OH, MH, and PH salinity zones in all seasons. NO 2 + NO 3 for the OH salinity zone showed increasing trends in all seasons from 1966 to 2015 as percent changes (Fig. 6c). A reversal of sign for percent changes from 1985 to 2015 occurred for the OH salinity zone shown at the top of individual bars. Percent changes for MH and PH salinity zones from 1966 to 2015 showed consistent, decreasing trends of NO 2 + NO 3 that continued throughout the time series.
trajectories of chl-a vs tN loading. Long-term data on chl-a were paired with annual TN loading to depict trajectories by salinity zone and season (Fig. 7a-f). Data for the OH salinity zone from 1964 to 2015 provided strong evidence of shifted trajectories of chl-a vs TN loading in spring and summer, consisting of higher chl-a in years with mid-range TN loading (1970,1971,1973,1974) than in recent years (2003,2004,2011) characterized by the highest annual TN loading in the time series (Fig. 7a,d). Much higher spring chl-a from 17-25 mg m −3 and summer chl-a from 22-35 mg m −3 occurred at TN loading of 60-80 (×10 6 ) kg yr −1 in the early 1970s, than spring chl-a from 10-13 mg m −3 and summer chl-a from 14-18 mg m −3 at higher TN loading >90 (×10 6 ) kg yr −1 in the 2000s (2003,2004,2011). Spring and summer chl-a in recent years with TN loading <50 (×10 6 ) kg yr −1 (2000)(2001)(2002)2013) were similar to chl-a at much higher TN loading in the 2000s (2003,2004,2011).
The trajectory of chl-a vs TN loading in spring for the MH salinity zone showed a weaker response of chl-a to TN loading than for the OH salinity zone (compare Fig. 7a,b). Two prominent spring chl-a maxima from 20-31 mg m −3 for the MH salinity zone occurred at moderate TN loading of 45-65 (×10 6 ) kg yr −1 in 1974 and 2013, contrasted with much lower chl-a from 2-7 mg m −3 at similar TN loading from 1964 to 1971 (Fig. 7b). The trajectory of chl-a vs TN loading in summer for the MH salinity zone showed chl-a was less sensitive to changes of TN loading than in spring, with chl-a from 6-15 mg m −3 at low to moderate TN loading of 25-60 (×10 6 ) kg yr −1 between 1964 and 1974 (compare Fig. 7b,e). Summer chl-a in this range was similar to chl-a from 9-13 mg m −3 at higher TN loading >90 (×10 6 ) kg yr −1 in the 2000s (2003,2004,2011) (Fig. 7e).
Water-quality criteria for major tributaries. Models for major tributaries in the bay were based on means, 10 th , and 90 th percentiles of salinity at each station as surrogates for flow to capture seasonal to inter-annual variability. Time series of salinity-adjusted model predictions of chl-a at nine tributary stations from 1985 to 2015 resembled analogous predictions for the main-stem bay (Fig. 8a-i). Model predictions in mean-salinity conditions showed increasing chl-a from 1985 to 2015 at the mouths of the Patuxent River (Fig. 8d) and James River (2019) 9:6709 | https://doi.org/10.1038/s41598-019-43036-6 www.nature.com/scientificreports www.nature.com/scientificreports/ (Fig. 8i), and at mid-estuary stations in the Choptank River (Fig. 8b), Potomac River (Fig. 8e), Rappahannock River (Fig. 8g), and James River (Fig. 8h). These long-term trends of chl-a were similar to those for MH and PH salinity zones in the main-stem bay (compare Fig. 4b,c). Observed vs model-fitted values of chl-a, Secchi depth, and NO 2 + NO 3 for nine tributary stations, and time series of salinity-adjusted model predictions of Secchi depth and NO 2 + NO 3 complementing analogous predictions of chl-a (Fig. 8a-i) are presented in Supplementary Material.
Numerical criteria for chl-a, Secchi depth, and NO 2 + NO 3 for nine tributary stations were based on salinity-adjusted model predictions (Fig. 8a-i; Supplementary Material, Figs S3, S5) aggregated at a monthly scale (Fig. 9a-c). Model predictions of chl-a in mean-salinity conditions ranged from 9.45-14.8 mg m −3 (Fig. 9a), compared to 5.36-14.7 mg m −3 for OH, MH, and PH salinity zones in the main-stem bay (Fig. 5a-c). Model predictions of Secchi depth in mean-salinity conditions ranged from 1.04-1.39 m at tributary stations (Fig. 9b), compared to 0.87-1.77 m for salinity zones in the main-stem bay (Fig. 5d-f). Lastly, model predictions of NO 2 + NO 3 in mean-salinity conditions ranged from 3.74-25.5 μM at tributary stations (Fig. 9c), compared to 0.690-57.4 μM for the main-stem bay (Fig. 5g-i). These proposed criteria for tributaries were similar to those for the main-stem bay. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Climatic effects on phytoplankton. An essential starting point in developing numerical water-quality criteria in the studies described here was to adjust for climatic effects. This allowed us to take account of spatio-temporal variability imposed by climatic effects, to distinguish long-term trends reflecting anthropogenic eutrophication, and to set realistic targets for restoration. Climatic effects on phytoplankton dynamics have been well described in estuarine-coastal ecosystems using long-term data from shipboard, aircraft, and satellite measurements [26][27][28][29][30][31][32][33][34][35] . Previous studies contributed to this understanding of climatic effects on nutrient loading, chl-a, floral composition, and NPP in Chesapeake Bay [17][18][19]24,33 and the Neuse and New River estuaries 34 .
A logical sequence emerged from these studies, accentuating hydrological regulation of TN and NO 2 + NO 3 loading with predictable consequences for water-quality properties and phytoplankton dynamics in the bay. Summarizing, low-flow, "dry" conditions lead to a landward shift of N-limitation toward OH and MH salinity zones, resulting in lower chl-a, lower NPP, and a decreased proportion of diatoms in the flora; conversely, high-flow, "wet" conditions extend the area of N sufficiency seaward to MH and PH salinity zones, resulting in higher chl-a, higher NPP, and an increased proportion of diatoms in the flora [17][18][19]24,33 . Dissolved and suspended materials affecting Secchi depth are similarly sensitive to climatic effects, with higher inputs of bio-optically active constituents in high-flow, "wet" conditions than in low-flow, "dry" conditions 17,18 .
We applied this logic to develop numerical water-quality criteria for the bay, using predictions conditioned on specific model inputs of flow and salinity to distinguish long-term trends from spatio-temporal variability imposed by climatic effects 17,18 . SRF and frequencies of predominant weather patterns identified "dry" and "wet" conditions ( Fig. 2) [19][20][21] ; salinity served as an explanatory variable in all statistical models, and as a proxy for Numerical water-quality criteria. Statistical models used to generate flow-adjusted predictions of chl-a, Secchi depth, and NO 2 + NO 3 allowed us: (1) to derive numerical water-quality criteria by season and salinity zone; (2) to adjust for climatic effects in establishing these criteria; (3) to propose criteria corresponding to climatic conditions reflecting decreased TN and NO 2 + NO 3 loading (MH and PH salinity zones) or increased light limitation (OH salinity zone); (4) to evaluate attainment by comparing recent values of water-quality properties to criteria based on selected time periods.
Multiple lines of scientific evidence, consisting of a historical increase of chl-a, low dissolved oxygen (DO), decreased water clarity, and harmful algal blooms (HAB), were used previously to develop numerical chl-a criteria for Chesapeake Bay 25 . We were guided by a "protective" approach to avoid impairments associated with nutrient over-enrichment and high chl-a. Sutula et al. 36,37 recently developed chl-a criteria for the San Francisco Bay estuary (SFB) using a similar approach designed to lessen the probability of ecosystem impairments such as HAB and low DO.
Numerical chl-a criteria developed in previous studies on Chesapeake Bay used a different statistical approach than we used here. We designated the 1960s as a "reference period" and computed geometric means and 90th percentiles for a period when symptoms of anthropogenic eutrophication were less evident 25 . Those studies proposed chl-a criteria for spring and summer with means from 1.4-15 mg m −3 as goals, and 90 th percentiles from 4.3-45 mg m −3 as thresholds. We later concluded that deriving numerical chl-a criteria with the 1960s as a reference period hinged on the validity of an assumption that water quality was better during that time. Although an absence of data for pristine conditions limited our options for a reference period, we believe that using the 1960s was too simple as: (1) mean, annual chl-a was higher for the OH salinity zone in the 1960s than in the mid-to late-1970s; (2) mean, annual chl-a was lower for MH and PH salinity zones in the 1960s, coinciding with persistent low-flow, "dry" conditions; (3) maxima of mean, annual chl-a for MH and PH salinity zones occurred by the mid-1980s following a decade of high flow, "wet" conditions; (4) spatio-temporal variability of mean, annual chl-a from 1995 to 2015 was driven by irregular "dry" and "wet" conditions. These several observations accentuated the need to adjust for climatic effects 17,18 rather than to base criteria on the 1960s as a reference period. www.nature.com/scientificreports www.nature.com/scientificreports/ To avoid limitations of using the 1960s as a reference period, our new analyses combined historical and monitoring data (1960s to 2015) to derive numerical criteria, focusing on climatic effects and long-term trends of water-quality properties. We based these analyses on data aggregated at monthly to seasonal scales, expanded the set of response variables to include chl-a, Secchi depth as a measure of water clarity, and NO 2 + NO 3 concentrations as a measure of nutrient over-enrichment, and considered reference periods specific to each water-quality property (Tables 4-6). Comparative data for nine tributary stations also supported criteria for chl-a, Secchi depth, and NO 2 + NO 3 revealing similar ranges for the main-stem bay and tributaries (Figs 5a-i and 9a-c, and Supplementary Material).
Proposed numerical criteria for chl-a, Secchi depth, and NO 2 + NO 3 presented in Tables 4-6 consist of specific values and their underlying bases. A conceptual diagram summarizes these criteria for the main-stem bay, providing a simple view based on salinity zone, season, and water-quality property (Fig. 10). These criteria improve www.nature.com/scientificreports www.nature.com/scientificreports/ upon earlier analyses by accounting for climatic effects and long-term trends of water-quality properties on seasonal and spatial bases.

Long-term trends, current status, and transitions.
Numerical water-quality criteria have practical applications to assess long-term trends, current status, and transitions. We previously reported a significant, decreasing trend of chl-a for the OH salinity zone after the mid-1960s using model predictions in mean-flow conditions for data aggregated at an annual scale. This trend was ascribed to a ban on PO 4 3− that enhanced P-limitation of phytoplankton in the upper bay 23 , leading to increased throughput of TN and NO 2 + NO 3 to MH and PH salinity zones 17,18 . New analyses of long-term trends suggest large reductions of TN and NO 2 + NO 3 loading in spring will be required to attain proposed chl-a criteria for OH, MH, and PH salinity zones (Table 4). A continuing, upward trend of chl-a for the MH salinity zone is especially problematic in the context of summertime hypoxia in the bay 3,16 .
Bio-optical conditions as Secchi depth have deteriorated throughout the bay, as shown by generally decreasing trends of model predictions for OH, MH, and PH salinity zones in mean-flow conditions from 1967 to 2015 (Fig. 6b). The largest percent changes of Secchi depth occurred for the MH salinity zone from 1967 to 2015, matched by trends of similar magnitude from 1985 to 2015. Some of the long-term trend of Secchi depth for the MH salinity zone can be explained by increased chl-a (Fig. 6a), but increased phytoplankton biomass does not account for decreasing trends of Secchi depth in the PH salinity zone (Fig. 6b) where chl-a has remained essentially constant from 1985 to 2015 (Fig. 6a). www.nature.com/scientificreports www.nature.com/scientificreports/ Our analyses suggest further reductions of TN and NO 2 + NO 3 loading will be required to attain chl-a criteria that constrain the magnitude and extent of the spring diatom bloom, particularly in mean-flow and high-flow years (Tables 4, 6). An important review of phytoplankton dynamics in the bay described NO 2 + NO 3 loading and density stratification associated with the spring freshet of the Susquehanna River as essential triggers of the spring bloom 35 . Consecutive features of the annual phytoplankton cycle include seasonal exhaustion of nutrients by the spring bloom, landward migration of the chl-a maximum, subsequent deposition of diatom biomass, and persistent thermal stratification in summer, culminating in deep-water DO depletion, and coinciding with a summer NPP maximum driven by N regeneration 18 . Decreasing trends of NO 2 + NO 3 for MH and PH salinity zones expressed as negative percent changes from 1985 to 2015 represent encouraging, albeit modest progress that is consistent with decreased TN and NO 2 + NO 3 loading since the early 1980s. Unfortunately, continued increases of model predictions of chl-a for the MH salinity zone in mean-flow conditions suggest recent progress to reduce TN and NO 2 + NO 3 loading has not been sufficient to lower phytoplankton biomass in the mid-bay (Fig. 4b), reversing eutrophication of the past half century 17,18 . trajectories of chl-a vs tN loading. Pronounced decadal differences in trajectories of chl-a vs TN loading expressed climatic effects and anthropogenic eutrophication in the bay 17 , exemplified by low chl-a per unit TN in the 1960s, high chl-a per unit TN in the 1970s, and relatively stable chl-a regardless of TN loading in the 2000s (Fig. 7a-f). These differences coincided with contrasting climatic conditions in the 1960s ("dry") and 1970s ("wet"), a doubling of flow-adjusted TN and NO 2 + NO 3 loading from 1945 to 1980 17 , decreased ratios of chl-a : TN from 1984 to 1992, and increased ratios in mean-flow conditions after 1994 17,18 . Notable features of trajectories included chl-a ~20 mg m −3 for the MH salinity zone that coincided with TN loading ~45 (×10 6 ) kg yr −1 in 2013, compared to chl-a < 6 mg m −3 at TN loading of 25-60 (×10 6 ) kg yr −1 in the mid-1960s to early 1970s (Fig. 7b). Similarly, chl-a ~15 mg m −3 for the PH salinity zone occurred at TN loading ~60 (×10 6 ) kg yr −1 in 1987, compared to chl-a < 4 mg m −3 at slightly higher TN loading in 1970 (Fig. 7c).
Trajectories of chl-a vs TN captured near-term climatic effects for the 2000s (Fig. 7a-f, amber shading). These data showed similar chl-a in "wet" years (e.g. 2003,2004,2011) with highest TN loading in the time series, and "dry" years (e.g. [2000][2001][2002]2013) with low TN loading similar to the 1960s. Thus, chl-a in recent years with reduced TN loading did not decrease to earlier, lower concentrations. These findings are consistent with the concept of "Return to Neverland" put forth by Duarte et al. 38 , postulating that a simple retracing of past trajectories is unlikely following a reversal of anthropogenic eutrophication. "Neverland" 38 was based on three-year means of chl-a vs TN loading, while trajectories of chl-a vs TN loading presented here were based on flow-adjusted annual means (Fig. 7a-f). Nonetheless, the non-linear responses of chl-a to reduced TN loading we observed in recent years support the applicability of this concept to the bay. Non-linear trajectories have been explained previously by long-term changes of ecosystem structure, including regeneration of legacy nutrients, changes of the cell-size distribution of phytoplankton, and altered floral composition 38,39 . Evidence for causes of non-linear responses of chl-a to TN loading in the bay consists of: (1) bio-optical properties with a change in coupling constants of Secchi depth and K D(PAR) that suggest the particle-size distribution may have shifted over time 40,41 ; (2)  Case studies: comparisons with other ecosystems. Anthropogenic eutrophication of estuarine-coastal ecosystems around the world has been extensively studied from the perspective of remediation. Analyses centered on two approaches to reverse effects of nutrient over-enrichment, "corrective" and "protective". The following paragraphs give examples of how these approaches have been implemented in several at-risk ecosystems with comparisons to Chesapeake Bay.
Novel statistical methods were used by Sutula et al. to develop chl-a criteria for SFB directed at lessening the likelihood of impairments 36,37 . SFB is characterized by high nutrient concentrations but has yet to experience deleterious impacts commonly associated with anthropogenic eutrophication. SFB has a shorter residence time than Chesapeake Bay, possibly explaining why potential impairments have yet to appear. A significant ecosystem shift occurred in SFB following introduction of the Asian clam, Potamocorbula amurensis, leading to increased benthic grazing that significantly reduced phytoplankton biomass in the northern estuary 12 . Some regions have not been affected by this invasive species, and resistance to nutrient over-enrichment elsewhere in SFB appears to be weakening, supported by observations of: (1) a three-fold increase of chl-a in the south bay during summer and fall since 1999 42 ; (2) common occurrences of HAB taxa throughout the estuary 43-45 ; (3) increased incidences of hypoxia with DO < 5 mg L −1 in the southernmost bay 46 .
These worrisome changes in SFB stimulated the development of chl-a thresholds intended to be "protective" from impairments, distinct from Chesapeake Bay where "corrective" measures are directed at reversing existing impairments. Significant relationships of chl-a to HAB abundance and low DO derived from quantile regressions indicated increased chl-a leads to an increased risk of impairments. Conditional probability analysis identified a chl-a threshold of 13 mg m −3 , below which probabilities of exceeding alert levels for HAB abundance and toxins decreased. A similar chl-a threshold of 13 to 16 mg m −3 was linked to a mandated water-quality criterion of 80% saturation for DO. Higher chl-a thresholds from 25-40 mg m −3 corresponded to 0.5 probability of exceeding alert levels for HAB abundance and DO < 5.0 mg L −1 in southerly regions of SFB. While these predictive relationships between chl-a, HAB, and low DO were sensitive to climatic effects and highly variable, "protective" chl-a thresholds provided starting points for SFB based on potential impairments 36,37 . Similar to SFB, our previous analyses to develop numerical chl-a criteria for Chesapeake Bay were based on multiple lines of scientific evidence, focusing on long-term trends and ecosystem impairments. It is noteworthy that comparable thresholds emerged for SFB 36,37 and Chesapeake Bay 25 , pointing to the merits of a "protective" approach to avoid negative consequences of anthropogenic eutrophication.
www.nature.com/scientificreports www.nature.com/scientificreports/ The Potomac River (PR) is the largest tributary of Chesapeake Bay. PR has been studied extensively since the 1960s following dense algal blooms that were stimulated by nutrient over-enrichment. PR is representative of mid-Atlantic estuaries in the United States that have experienced moderate to high levels of degradation. Riverine parts of PR have undergone extensive hydrologic modifications in the Washington, D.C. urban area. The large sub-estuary drains an extensive watershed covering several states and is highly responsive to climatic effects, such as droughts and floods, that influence the entire region. Long-term studies of PR by Jaworski et al. 47 and a synthesis by Bricker et al. 48 documented significant impacts of anthropogenic eutrophication, despite decreasing trends of N loading from the upper river basin and declining concentrations of NO 3 in surface waters. To this point, Jaworski et al. 47 estimated a 50% reduction of 1985 base-year TN loading, including 54-65% reductions of non-point sources and continued reductions of TN in wastewater effluent, would be required to meet water-quality criteria for PR. Reductions of TN loading of similar magnitude have been recommended for Chesapeake Bay, but progress to attain this goal has been modest. Other approaches in PR such as shellfish aquaculture have been proposed to complement land-based measures to decrease TN loading, but oyster restoration alone is deemed unlikely to reverse symptoms of anthropogenic eutrophication in either the sub-estuary 48 or the main-stem bay 49 .
In northern Europe, extensive management efforts have been directed at Danish coastal waters to reverse deleterious symptoms of anthropogenic eutrophication. Widespread hypoxia in the Danish straits stimulated the 1985 NPo Action Plan, and long-term studies chronicled successful reductions of nutrient loading 13,38,39,50,51 . Riemann et al. 13 described 25 years of water-quality responses in a "corrective" approach that has produced ~50% decreases of N and P loading since 1990. These dramatic decreases of N and P loading led to decreased nutrient concentrations in receiving waters, a modest decrease of chl-a, and restoration of macro-algae to deep waters. Notably, the decrease of chl-a has failed to match nutrient reductions stoichiometrically 13 , consistent with observations for other estuarine-coastal ecosystems detailed by Duarte et al. 38 and Carstensen et al. 39 Nonetheless, successful corrective actions in Danish coastal waters give hope that significant reductions of N and P inputs may reverse deleterious symptoms of anthropogenic eutrophication in other estuarine-coastal ecosystems. A key difference between Danish coastal waters and Chesapeake Bay is that decreasing trends of nutrients in the N. European ecosystem entailed changes in the relationships of nutrient loading to freshwater flow 50,51 , while relationships of TN and NO 2 + NO 3 loading to SRF in the bay have not returned to previous conditions 17,18 . expectations, restoration, possibilities. How might we use climatic effects on water-quality properties to develop numerical criteria that inform us about targets for restoration? We submit that flow-adjusted model predictions constitute guidance to identify attainable values for water-quality properties by incorporating climatic effects. In this way, numerical criteria for chl-a, Secchi depth, and NO 2 + NO 3 bracket realistic, pragmatic goals that can be compared to current conditions to gauge progress. Too often, bay-health assessments have ignored spatio-temporal variability associated with climatic effects. Our analysis of long-term data accentuates the sensitivity of water-quality properties to climatic effects, showing that assessments based on one or several years can be misleading. To this point, claims of progress toward improved water quality following low-flow, "dry" conditions are generally premature or erroneous. An apt analogy would be skepticism about "global warming" based on periodic cold winters experienced in a historical context of rising temperatures.
A shallow trajectory toward reduced TN and NO 2 + NO 3 loading in the bay 17,18 has been accompanied by troubling trajectories of chl-a vs TN loading (Fig. 7a-f), confirming that a simple return to previous conditions may prove elusive 38 . An increasing trend of model predictions of chl-a for the MH salinity zone in mean-flow conditions (Fig. 4b) provides worrisome evidence that organic matter derived from phytoplankton has not responded to decreased TN and NO 2 + NO 3 loading in the region where annual hypoxia/anoxia occurs 3,[16][17][18]24,25 . Successful efforts in Danish coastal waters produced significant, 50% reductions of nutrient inputs and changed relationships between nutrient concentrations and flow 50,51 . In contrast, modest progress to reduce nutrient inputs in the bay has consisted of <20% reductions of TN and NO 2 + NO 3 loading, and relationships between loading and SRF have not returned to previous conditions 17 . Absent such a change, climatic effects will continue to dominate seasonal to inter-annual variability of TN and NO 2 + NO 3 loading, and a "corrective" approach alone is unlikely to yield significant improvements of water quality. We suggest that specifying criteria based on sustained adherence to flow-adjusted model predictions represents an approach that takes advantage of time-series data on water-quality properties, while incorporating climatic effects that strongly influence contemporary conditions.

Conclusions
Numerical water-quality criteria were developed for Chesapeake Bay using statistical models to adjust for climatic effects. These criteria were directed at a "corrective" approach, setting proposed criteria in a domain of the reasonable, based on long-term observations. Summarizing: • Spatio-temporal variability of water-quality properties exemplified by chl-a, Secchi depth, and NO 2 + NO 3 is driven primarily by climatic effects, superimposed on long-term trends associated with anthropogenic eutrophication; • Flow-adjusted TN and NO 2 + NO 3 loading to the bay captures the course of anthropogenic eutrophication since World War II, providing a rationale to adjust time series of water-quality properties for climatic effects; • Statistical models applied to time series of chl-a, Secchi depth, and NO 2 + NO 3 supported numerical water-quality criteria for the main-stem bay and major tributaries; • Flow-adjusted model predictions were used to compute long-term trends of water-quality properties, to specify numerical criteria constituting realistic goals, and to assess progress toward attainment using comparisons with conditions in recent years; www.nature.com/scientificreports www.nature.com/scientificreports/ • This 'corrective' approach based on numerical criteria extends work on other estuarine-coastal ecosystems to incorporate climatic effects, thereby addressing spatio-temporal variability, resolving long-term trends, and quantifying improvements.
Methods study site. The focus of these studies was Chesapeake Bay in the mid-Atlantic region of the United States.
The bay is a shallow, partially mixed, temperate estuary of the Susquehanna River, with a main-stem surface area ~8,000 km 2 , receiving inputs of freshwater, sediment, and solutes from an extensive 165,000 km 2 watershed. North-south gradients of salinity, nutrients, and light penetration characterize the ecosystem, with a number of significant tributaries also contributing freshwater and solutes. These tributaries include the Patapsco, Patuxent, Potomac, Rappahannock, York, and James Rivers on the western shore, and the Choptank, Pocomoke, and Nanticoke Rivers on the eastern shore. A map showing major geographic features, salinity zones, and station locations ( Fig. 1) was produced with Surfer v. 8 (Golden Software) and customized with Adobe Photoshop v. CS6.
Data sources. Long-term data on water-quality properties for the bay and tributaries from 1960 to 2015