Extent and Causes of Chesapeake Bay Warming

Coastal environments such as the Chesapeake Bay have long been impacted by eutrophication stressors resulting from human activities, and these impacts are now being compounded by global warming trends. However, there are few studies documenting long‐term estuarine temperature change and the relative contributions of rivers, the atmosphere, and the ocean. In this study, Chesapeake Bay warming, since 1985, is quantified using a combination of cruise observations and model outputs, and the relative contributions to that warming are estimated via numerical sensitivity experiments with a watershed–estuarine modeling system. Throughout the Bay’s main stem, similar warming rates are found at the surface and bottom between the late 1980s and late 2010s (0.02 ± 0.02°C/year, mean ± 1 standard error), with elevated summer rates (0.04 ± 0.01°C/year) and lower rates of winter warming (0.01 ± 0.01°C/year). Most (~85%) of this estuarine warming is driven by atmospheric effects. The secondary influence of ocean warming increases with proximity to the Bay mouth, where it accounts for more than half of summer warming in bottom waters. Sea level rise has slightly reduced summer warming, and the influence of riverine warming has been limited to the heads of tidal tributaries. Future rates of warming in Chesapeake Bay will depend not only on global atmospheric trends, but also on regional circulation patterns in mid‐Atlantic waters, which are currently warming faster than the atmosphere.


INTRODUCTION
Increasing temperatures in the marine environment are a globally recognized phenomenon, with the world's oceans accounting for approximately 93% of the total increase in global heat content due to the effects of anthropogenic warming in the atmosphere since 1970 (Levitus et al. 2012;Rhein et al. 2013).
Compared to the open ocean, warming patterns in the coastal ocean are much less uniform spatially and temporally (Belkin 2009;Lima and Wethey 2012;Liao et al. 2015;Alexander et al. 2018), which can make them more difficult to detect (Henson et al. 2016). Warming in coastal habitats is known to affect rates of biogeochemical cycling that can drive negative impacts like hypoxia (Breitburg et al. 2018;Fennel and Testa 2019). Rising temperatures over the past few decades have been linked to worsened water quality conditions in well-studied coastal environments, such as the Baltic Sea (Carstensen et al. 2014;Lennartz et al. 2014) and the Chesapeake Bay (Du et al. 2018). As global temperatures continue to increase, projections of regional impacts on coastal environments have sought to assess how future warming may continue to negatively impact water quality by compounding effects of eutrophication in locations like the Chesapeake Bay (Irby et al. 2018;Ni et al. 2019), the Baltic Sea (Saraiva et al. 2019;W ahlstr€ om et al. 2020), and the Gulf of Mexico ), among many others (Breitburg et al. 2018). The Chesapeake Bay is a partially mixed, highly productive estuary that is the largest in the continental United States (U.S.) (Kemp et al. 2005;Cloern et al. 2014). Over the past several decades, substantial changes to Chesapeake Bay waters have been documented, including the extraction of natural resources, increased sediment and nutrient inputs, loss of historic habitat, and substantial increases in annual hypoxia and anoxia (Hagy et al. 2004;Kemp et al. 2005;Bever et al. 2013;Pan et al. 2021). Intensive management efforts to reduce excessive nutrient and sediment runoff entering the Bay (Ator et al. 2020) are outlined by a Total Maximum Daily Load (TMDL; USEPA 2010) and the estuary is home to a comprehensive long-term observational network that tracks physical, chemical, and biological parameters (Olson 2012;Chesapeake Bay Program DataHub 2020). Determining the role of warming in this highly managed estuary is important when assessing projected changes in the distribution of numerous biologically, commercially, and recreationally important species. Quantifying the degree of long-term warming in the Bay, as well as the associated mechanisms, will also help improve projections of future estuarine conditions. This is critical information for managers who need to assess future water quality goal attainment. It is highly likely that additional nutrient reductions will be necessary in the future to mitigate the impacts of temperature increases throughout the Bay (Irby and Friedrichs 2019).
For context, historical rates of warming over past centuries in the Chesapeake Bay are markedly lower than more recent estimates. Long-term warming events in the paleotemperature record have been documented, such as a 10°C increase in Bay temperatures over a period of 2,000 years following changes in the Gulf Stream position (Balsam and Heusser 1976); however, temperature records over the past~2,000 years have shown that the 20th Century Chesapeake Bay experienced more rapid rates of warming than at any time in the previous two millennia (Cronin et al. 2003). Important features of Holocene climate in the Bay region included multidecadal oscillatory cycles (substantially less influenced by anthropogenic activity) as well as changes in circulation along the continental shelf of the U.S. East Coast (Cronin et al. 2003(Cronin et al. , 2010. Warming in the Mid-Atlantic Bight is also a contemporary phenomenon that may be a result of shifts in the Gulf Stream position (Shearman and Lentz 2010;Andres 2016;Alexander et al. 2020) and could have implications for current shelf warming (Forsyth et al. 2015;Wallace et al. 2018;Chen et al. 2020). Studies of recent Chesapeake Bay temperatures have provided various estimates of warming using a variety of methods. Projected surface warming based on records of early 20th Century Baltimore Harbor data from a handful of observing stations was posited (and discounted) by Brady (1976), who suggested that, by 2012, average Bay temperatures could be 2.2°C higher than those observed in the first half of the 20th Century. Increasing temperatures (0.8°C-1.1°C) were found Bay-wide at the surface during winter and spring and at the bottom during winter, spring, and summer in the latter half of the 20th Century (Preston 2004). Longer and more continuous records at two surface-water observing stations also showed longterm warming over this time period (Najjar et al. 2010). Using satellite data from 1984 to 2010, Ding and Elmore (2015) showed that annual warming patterns in northern Bay surface waters were generally higher in the main stem than the tributaries, suggesting a possible oceanic influence.
Despite the relative breadth of research on the warming of the Chesapeake Bay, there are pressing issues that remain unresolved and have significant consequences for long-term Bay management. To date, it remains unclear whether significant temperature trends can be estimated from the available monthly in situ data, or whether the observational network may have implicit biases that could limit our understanding of past warming. Such biases may include interannual variability in weather conditions or differences in the timing of sampling, especially with sampling frequently being postponed until after major wind or storm events, when the water column is well mixed. Additionally, spatial discrepancies in the rates of increasing air and water temperatures (Ding and Elmore 2015) suggest that there is a need to determine a mechanistic connection between drivers of temperature change in the Chesapeake Bay and the relative responses over time. Identifying causal mechanisms that govern observed warming will provide greater insights into how future climate projections may evolve.
Mechanistic models can help resolve these issues by complementing analyses of temperature trends derived from in situ data and identifying the relative drivers of these trends. Such models have been applied previously to other coastal systems to gauge the importance of different temperature stressors. Kniebusch et al. (2019) used a numerical modeling framework of the Baltic Sea to attribute and quantify the relative sources of observed warming since the 1850s. Sensitivity tests of water temperature influences in the San Francisco Bay delta area have been evaluated over recent years to help inform projections of future changes driven by global climate trends (Vroom et al. 2017). To date, the authors are unaware of published research quantifying the impacts of causal mechanisms on temperature change within the Chesapeake Bay.
This study utilizes the extensive observational network of in situ data along with a watershed-estuarine modeling system forced by realistic atmospheric and oceanic inputs to quantify and better understand the causes of warming in the Chesapeake Bay. Using this approach, a more robust estimate of the recent observed temperature trends and the causality of said trends can be more precisely determined. This study provides a comprehensive estimate of temperature changes that can be applied to the 2025 TMDL (USEPA 2010) management timeline, and offers a retrospection of anthropogenically influenced warming that can inform future assumptions of temperature patterns and associated biological impacts.

Chesapeake Bay Program Data
The Chesapeake Bay is home to an extensive observational network that, for more than three decades, has collected monthly and semi-monthly physical and biogeochemical in situ data. This CBP Water Quality Database (Olson 2012; Chesapeake Bay Program DataHub 2020), available through an online web portal http://data.chesapeakebay.net/WaterQuality), acts as a repository for temperature and salinity in situ data collected at numerous main stem stations since 1984. Temperatures and salinities are typically measured from the surface to the bottom at one to 2-m intervals using YSI 6600 sensors that have manufacturer-reported accuracies of AE0.15°C and AE1%, respectively (Olson 2012). The temporal frequency of sampling typically includes twice-monthly observations from April to September and monthly observations from October to March; however, in reality, the observational sampling frequency can vary based on weather activity within the Bay. This study uses observational data from the 20 main stem stations ( Figure 1) that are most consistently sampled and encompass the full length of the main stem Bay's deep channel.

Linked Watershed-Estuarine-Hydrodynamic Model
A three-dimensional hydrodynamic model is used to determine the causal mechanisms of observed Bay warming. The modeling system is an implementation of the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams 2005) developed for the Chesapeake Bay (ChesROMS; Xu et al. 2012). Ches-ROMS utilizes a curvilinear grid with an average horizontal grid cell resolution of approximately 1.8 km in the main stem Bay and has 20 topographyfollowing vertical levels ( Figure 1). The version of the model used here has been described and evaluated in Oligohaline N Mesohaline S Mesohaline Polyhaline previous research (Feng et al. 2015;Irby et al. 2016;Da et al. 2018;St-Laurent et al. 2020). Minor updates have been made, including a modification of the Jerlov water type parameter (Jerlov 1976;Paulson and Simpson 1977), which determines the attenuation of solar radiation as a function of depth. This parameter was increased from Jerlov type 3 to 5 in order to more accurately represent estuarine temperatures in the main stem Bay (Kim et al. 2020). ChesROMS receives daily riverine inputs of freshwater discharge and temperature from the CBP's Watershed Model for the years 1985-2014 (Shenk and Linker, 2013). The latest iteration of the Watershed Model includes a refined representation of observed changes in land use and management practices that have occurred over this time period (Easton et al. 2017). When Watershed Model results are unavailable (after 2014), an annual climatology of daily riverine temperatures (2010-2014) was used. In addition, daily U.S. Geological Survey (USGS) freshwater discharge estimates (USGS 2020) were scaled by river-specific factors derived from the Watershed Model, in order to correct for discharge occurring downstream of the USGS estimates.
Using the linked Watershed-Estuarine-Hydrodynamic model described above, a reference run was conducted for the years 1985-2019 using realistic forcing (Model Atmospheric and Oceanic Forcing). This reference run was extensively evaluated with Chesapeake Bay data (Model Evaluation) and used as a baseline for the sensitivity experiments designed to isolate the mechanisms driving temperature trends in the Bay (Model Sensitivity Experiments).

Model Atmospheric and Oceanic Forcing
In addition to forcing at river mouths, the ROMS reference run received forcing from the atmosphere and open ocean. Atmospheric forcing includes surface air temperature, downwelling longwave radiation, net shortwave radiation, precipitation, relative humidity, winds, and air pressure obtained from ERA5, a fifthgeneration atmospheric reanalysis dataset with a resolution of 0.25°(C3S 2017). Open ocean forcing includes temperature, salinity, and sea surface height. Open ocean temperature and salinity for the reference run were derived from the World Ocean Database (Boyer et al. 2018). Data from 2008 to 2018 were used to produce monthly climatologies at each grid cell along the model's oceanic boundary, and were assumed to be representative of a central year (2013). The 2013 baseline was then projected backward to 1985 and forward through 2019 using longterm linear trends derived from the same database but utilizing a longer period of time. Specifically, ocean boundary temperature changes were calculated by performing seasonal linear regressions of depthaveraged (upper 35 m) World Ocean Database values at the model's open boundary, spanning the years 1985-2018 and a latitudinal interval from 36°to 37.8°N. Trends thus vary seasonally but not spatially across the model boundary ( Figure 2c) and agree well with previous estimates of continental shelf warming (Wallace et al. 2018). As in previous ChesROMS implementations St-Laurent et al. 2020), sea surface height forcing at the open boundary was derived using hourly observations from coastal stations at Duck, North Carolina and Lewes, Delaware and tidal harmonics from the Advanced Circulation model (Luettich et al. 1992).

Model Evaluation
Modeled water temperatures and salinities were statistically evaluated at the surface and bottom against long-term in situ data collected at the main stem stations (Figure 1) grouped by salinity regime: oligohaline (salinity <13), northern mesohaline (13 ≤ salinity < 18), southern mesohaline (18 ≤ salinity < 22), and polyhaline (salinity ≥22), as well as for the Bay as a whole. Hourly model outputs from the reference run were selected that were nearest in time to each individual recorded observation and spatially interpolated to the station's location using the four neighboring grid cells. Statistical metrics of model skill are summarized using target diagrams, which illustrate bias, unbiased root-mean squared difference (RMSD) and Total RMSD (note that unbiased RMSD is defined as the Total RMSD computed after removing the mean, and mathematically Total RMSD is equal to the square root of the sum of the bias squared and unbiased RMSD squared, Hofmann et al. 2008;Jolliff et al. 2009). Here these statistics are normalized by the standard deviation of the observations, rendering them unitless and enabling the representation of multiple variables on a single diagram (Jolliff et al. 2009). In target diagrams, the distance along the x-and y-axes denote, respectively, the normalized unbiased RMSD and normalized bias. The distance of a given symbol from the origin mathematically represents total normalized RMSD; thus, a perfect model fit would be represented by a symbol at the origin (0,0).

Observed and Modeled Estimates of Long-Term Change in Estuarine Temperature
The overall rate of observed estuarine bottom and surface temperature change over the study period was computed for the full Bay as well as for individual salinity regimes. Temperature trends were calculated using a first-order autoregressive (AR1) model that accounted for interannual autocorrelation using the Cochrane-Orcutt methodology (Cochrane and Orcutt 1949). This AR1 model was applied to monthly averaged data from the 20 stations used in this analysis (Figure 1), to produce temperature change estimates from observations (referred to as Obs(stn)). For months that recorded a single observation, that individual point was used as the monthly value; for months that included two observations, the average of the two was used as the representative value for that month at a given station. Trends calculated from these spatially grouped monthly values were averaged to produce seasonal (May-October and November-April) and annual estimates of change.
These temperature trends were also estimated in two ways from the reference run. In order to better understand the limitations of assuming that a single temperature observation (or average of two observations) represents the monthly average temperature, the reference run was used to estimate the temperature trend using model results that matched the time and location of the available observations (i.e., Model-Ref(stn)) as well as using the monthly and spatial average of the main stem grid cells (ModelRef(all)) whose average depths are ≥5 m. Comparison of Mod-elRef(stn) and ModelRef(all) allows an assessment of (1) aliasing in the data resulting from the fact that data are collected at most twice per month and sampling is frequently influenced by weather conditions (i.e., cruises are postponed until after storms have passed) and (2) the limitation of a finite number of stations to represent salinity regimes and the bay as a whole. The difference in results between ModelRef (stn) and ModelRef(all) provides a further uncertainty estimate on how well the observationally derived trends (Obs(stn)) represent actual Bay warming. Standard errors were computed for all trend estimates calculated from observations and the reference run (Obs(stn), ModelRef(stn), and ModelRef(all)).

Model Sensitivity Experiments
The first five years of the reference run (1985)(1986)(1987)(1988)(1989), BASE) served as the baseline for five temperature sensitivity experiments. These experiments were designed to isolate the relative magnitude of the mechanisms responsible for the long-term changes in Seasonal variability in long-term change (1985-2019) computed using linear regression for (a) air temperature, (b) downwelling longwave radiation, and (c) ocean boundary temperatures. Boxplot whiskers in (a,b) represent the range of spatial variability in atmospheric trends based on the entirety of grid cells used for ChesROMS atmospheric forcings (n = 238). Seasonal ocean temperature trends were interpolated by month and represent the average change along the ChesROMS model open boundary, and have no spatial variability. Boxplots (a,b) show the median (red line), 25th and 75th percentiles (top and bottom edges of blue box), and whiskers that cover the full range of calculated trends. Errorbars (c) represent seasonal standard errors averaged across the model's ocean boundary.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
estuarine temperature, as well as their combined effects (Table 1). Each scenario had identical forcings to the reference run (i.e., discharge, precipitation, winds, humidity, etc.), except for the modification of a particular forcing field that was altered to represent conditions of a time period 30 years later (2015-2019). These experiments altered: (1) air temperature and downwelling longwave radiation (AtmTemp), (2) ocean temperature (OceanTemp), (3) sea level (SeaLevel), and (4) river temperature (RiverTemp). A final experiment combined all four changes simultaneously (Combined). Both air temperature and downwelling long-wave radiation were changed in the same experiment, because these changes are not independent and are highly correlated in space and time. Prior to each experiment, there was a two-year simulation period that allowed the model to adjust to the temperature forcing perturbation.
The perturbations used in the five sensitivity experiments were obtained by multiplying the seasonally varying long-term trends ( Figure 2) by a period of 30 years. Temporal spline fits were applied to the monthly changes in air temperature, longwave radiation, and ocean temperature in order to apply smooth daily changes in these forcing fields. Two of these atmospheric fields were associated with generally positive long-term trends over the time period 1985-2019: air temperature (0.02 AE 0.01°C/year) and downwelling longwave radiation (0.21 AE 0.06 W/m 2 / year); however, these trends were also characterized by considerable seasonal variability (Figure 2a, 2b). The air temperature trends agree well with previous estimates derived from annual weather station data (Ding and Elmore 2015) and data from the U.S. Historical Climate Network (USHCN; Menne et al. 2009; see Figure S1). This agreement also extends to the cooling trends present in March and November (Figure 2a), which are clearly visible in the observations ( Figure S1). Over this 30-year period, a linear regression of sea surface heights averaged across the open boundary showed an increase of 0.15 AE 0.02 m (4.9 AE 0.5 mm/year), consistent with estimates from Boon and Mitchell (2015) and Ezer and Atkinson (2015). Changes in river temperature used previous estimates of average long-term USGS stream temperature trends (0.28°C per decade) throughout the Chesapeake Bay watershed (Rice and Jastram 2015); over the 30-year period these amounted to +0.84°C (Table 1). These changes in sea surface height and river temperatures were both added as constant values, as there was no clear spatial variability in these trends. In each experiment, the 30-year perturbations were applied daily to the 1985-1989 baseline simulation.
Finally, for comparison with the temperature trends described above (Obs(stn), ModelRef(stn), and ModelRef(all)), temperature trends were also computed from the Combined sensitivity experiment. This estimate relies on the same spatial extent and monthly averages as in the ModelRef(all) trends. Using the same atmospheric and riverine forcing altered by a delta computed from the long-term observed trend estimates, this final estimate of estuarine temperature change removes the impact of shorter term changes in rainfall and temperature (e.g., the impact of a recent year being anomalously wet/dry or hot/cold.)

Model Skill Assessment
The reference simulation (1985-2019) demonstrates considerable skill in capturing temporal variability in temperature and salinity along the Chesapeake Bay's main stem (Table 2). Throughout the simulation period ChesROMS slightly overestimates temperatures by 0.6°C at both the surface and bottom, and this bias is lower in the southern Bay than the northern regions. The RMSD of surface temperatures (0.7°C) is slightly lower than that for bottom temperatures (0.8°C). Both represent a relatively small model bias compared to the range of observed Bay temperatures over the course of an average year (standard deviation = 8°C). Model skill in simulating temperatures varies both temporally (Figure 3a, 3b) and spatially throughout the Bay (Figure 3c, 3d). Temporal variability in temperature is well reproduced throughout the Bay at both the surface and bottom, with normalized RMSD and bias typically <0.25. There is relatively little temporal variability in temperature throughout the Bay compared to spatial variability; therefore, it is not surprising that the normalized temporal RMSD (Figure 3a, 3b) is lower than normalized spatial RMSD (Figure 3c, 3d). Model results are also biased higher in the spring and summer months at the surface, and in the winter and spring at the bottom.
Modeled salinities agree well with observations ( Table 2), suggesting that ChesROMS is simulating hydrodynamic processes with reasonably high skill. Average surface and bottom salinities are very near the average value of observations, with biases of approximately À0.03 and 0.25 units, respectively (Table 2). The model generally overestimates surface salinity in the oligohaline and underestimates it in the polyhaline. At the bottom, the model underestimates salinity in the northern Bay regions and overestimates salinity in the southern Bay regions. Although there are no distinct regional patterns in terms of overall skill, Total RMSD is lowest in the northern mesohaline at the surface and bottom.
The model also has variable model skill in terms of its ability to reproduce the observed temporal and spatial patterns of salinity ( Figure 3). Seasonal variability in salinity is greater than that of temperature, and can be more difficult for the model to reproduce, particularly at the bottom in the southern mesohaline region (Figure 3a, 3b). Modeled salinities have a lower temporal bias at the surface than at the bottom although there is a greater range of RMSD in surface waters. In terms of spatial variability, salinity skill is similar at the surface and bottom and is greater than temperature skill in all months of the year with normalized RMSD and normalized bias both <0.2 (Figure 3c, 3d).

Observed and Modeled Bay Temperature Trends
In general, all three estimates of temperature change -Obs(stn), ModelRef(stn), and ModelRef(all) suggest that both surface and bottom warming are strongest in the months between mid-spring and early fall, whereas weaker levels of long-term warming, or even potential cooling, occur in the cooler months of the year (Table 3; Figure S1 and Figure 4). The observed and modeled trend error bars overlap for all months at both the surface and bottom, with the exception of the month of June at the bottom when the observed trend is considerably larger (Obs (stn) = 0.06°C/year) than the modeled trend (Model- Figure 4). The most substantial differences between modeled trend estimates occur in October at both the surface and bottom, when ModelRef(stn) is significantly less than ModelRef(all) (Figure 4). Over the 35-year study time period, the total change in bottom temperature calculated from observations (Obs(stn)) and model results (ModelRef(all)) was similar: 0.77 AE 0.32°C and 0.60 AE 0.28°C, respectively (Table 3). Over this same time period, estimates of surface temperature increases were estimated to be similar to those at the TABLE 2. Summary of model skill in simulating observed monthly regional and full Bay main stem temperatures (°C) and salinities (unitless). Number of data points = n 9 12 months 9 35 years, where n = 4 for the oligohaline and polyhaline regions, n = 6 for the north mesohaline and south mesohaline regions, and n = 20 for the whole bay. and ModelRef(all), respectively. All three trend calculations also indicate that both the mean and seasonal cycle of long-term temperature change is fairly consistent across all Bay regions ( Figure S2).

Results of Sensitivity Experiments
The results of the sensitivity experiments assessing the total and individual impacts of the four mechanisms responsible for the increased estuarine temperatures are discussed below. For all sensitivity results, estimates of temperature change are accompanied by their respective standard deviations calculated in the same way as the Combined sensitivity experiment (Model Sensitivity Experiments). the surface and the bottom, with an overall change of 0.7 AE 0.5°C over the full 30-year period studied, where the numbers following the AE sign refer to the standard deviation of the temperature change across all the main stem grid points (Tables 4 and 5; Figure 5). At both the surface ( Figure 6; Table 4) and the bottom ( Figure 6; Regionally, surface-to-bottom differences in warming are most apparent in the polyhaline region of the lower Bay main stem, where the average increase in bottom water temperatures is~0.1°C greater than that at the surface ( Figure 5). Overall changes in average annual temperatures range from 0.65 AE 0.38°C in oligohaline bottom waters to 0.84 AE 0.55°C in polyhaline bottom waters and at the surface from 0.72 AE 0.43°C in oligohaline waters to 0.74 AE 0.50°C in polyhaline waters (Tables 4 and 5; Figure 5). Greater differences between surface  Red circles denote temperature change calculated solely using observational data (see Figure 1 for station locations). Blue circles represent temperature change based on model outputs extracted from the same time and location as observations. Black diamonds represent temperature change calculated using daily model results from all grid cells in the main stem Bay. Error bars represent AE1 standard error for the spatially averaged full Bay region (temperature changes broken apart by region are shown in Figure S2). and bottom warming are more apparent in certain months than others ( Figure 6). In cooler months when there is relatively little bottom warming, there is a smaller difference between increases in surface and bottom temperatures (Tables 4 and 5).
There is a greater temperature increase in surface waters than bottom waters in the spring season (April-May, Figure 6), when there are also substantial trends in atmospheric temperatures and downwelling longwave radiation (Figure 2a, 2b). In addition, the greatest amount of bottom warming occurs in summer months, when it can be twice as high as surface warming. On average during the period from May to October, a time of particular interest since it encompasses the bottom hypoxia season, there is more warming in the shallower, southernmost extent of the main stem than in the rest of the Bay (Figure 6).
Experiments: AtmTemp, OceanTemp, SeaLevel, RiverTemp. Together, increasing air temperatures and downwelling longwave radiation substantially contribute to the total increase in Bay water temperatures across all regions of the Bay, at  ure 7; AtmTemp = red line). At the surface, atmospheric changes are responsible for the majority of total temperature change (Combined) in all regions of the Bay. This is not the case at the bottom, where the contribution of atmospheric changes to total estuarine warming is reduced in the southern half of the Bay, particularly during summer months (Figure 7g, 7h). The influence of increasing ocean temperatures in the main stem is substantial in certain seasons and locations (Figure 7; OceanTemp = orange line). The greatest impact of ocean warming is from July to September (Figure 7). At the surface, this impact is never greater than the impact of atmospheric change (AtmTemp), but at the bottom, the impact of ocean temperature change exceeds that of atmospheric change in both the southern mesohaline and polyhaline regions during summer months (0.33 AE 0.26°C and 0.69 AE 0.49°C, respectively; Table 5). Increasing ocean temperatures exert a greater impact on surface and bottom warming for a longer period of time with increasing proximity to the Bay mouth, with the peak of its influence typically occurring in August for all regions (Figure 7). In the southern mesohaline region, ocean warming accounts for over half of the bottom warming during July and August, and for nearly all of the combined changes in the polyhaline region during these same months (Figure 7g, 7h).
Sea level rise also contributes to long-term estuarine temperature change. At the surface, the effect of sea level rise has little to no impact on water temperatures in any region of the Bay apart from the oligohaline region, where it contributes a slight cooling effect in summer months and warming effect in winter months (Table 4; Figure 7a; SeaLevel = light blue line). The summer cooling effect of sea level rise is more pronounced at the bottom where its impacts are evident from the spring to early fall throughout the Bay (Figure 7e-7h; Table 5); however, the average Bay-wide contribution to summer bottom water temperature change is on the order of À0.1°C (Table 5). In all regions of the Bay at both the surface and the bottom, the timing of the maximum cooling effect of sea level rise occurs approximately one to two months earlier than the maximum warming effect due to oceanic change (Figure 7).
River warming plays only a small role in increasing main stem surface and bottom temperatures (Figure 7; RiverTemp = green line); however, a small degree of warming can be seen in the upper parts of the largest tributaries (Susquehanna, Potomac and James; Figure 5). In the oligohaline region, increasing river temperatures impact the main stem warming mostly at the surface when discharge is greatest Surface Bottom (a) (b) 30-year change in temperature, ºC  Table 4; Figure 7). When the discharge is lower during summer months, the effect of increasing river temperatures in the oligohaline region is slightly reduced (Figure 7; Tables 4 and 5). The contribution of individual warming mechanisms can be summarized by their percent contribution to total long-term change in Bay temperatures in warmer months (Figure 8a-8d) and cooler months (Figure 8e-8h). This percentage represents the relative change in temperature of an individual experiment (AtmTemp, OceanTemp, SeaLevel, RiverTemp) divided by the total temperature change (Combined). At the surface, long-term temperature change is driven almost completely by atmospheric changes, except in the summer in the southern Bay (Figure 7d) where ocean temperature change plays a substantial role. Bottom waters are primarily warmed by changes in atmospheric conditions (78%-87%), with a smaller contribution from increasing ocean temperatures in cooler months (9%; Figure 8f) that becomes more pronounced during the summer (26%), particularly in mesohaline and polyhaline Bay regions (Figure 8b). A slight cooling effect from sea level rise exists during the warmer months throughout the Bay that produces a net negative contribution to warming bottom waters (À6%; Figure 8c). The contributions of river-driven warming are largely confined to locations near the heads of major tributaries that do not greatly impact the main stem (<2%; Figure 8d, 8h).
Altogether, the summation of temperature changes in each individual sensitivity experiment is approximately equal to the combination of all warming mechanisms together, suggesting that the applied changes are independent. When comparing the temperature change over the entire Bay in the Combined experiment and the summation of temperature changes in the four sensitivity experiments, the differences are less than or equal to AE0.03°C (Tables 4 and 5, final column). When looking at individual regions, depths, and seasons, these differences are greatest in southern Bay bottom waters during warmer months, but are still less than or equal to AE0.05°C (Table 5).

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
within the Chesapeake Bay is critical for improving our understanding of future estuarine warming and its impacts on hypoxic conditions throughout the 21st century. Here this is obtained by combining longterm temperature observations with output from a linked watershed-estuarine mechanistic model. The Bay is typically sampled monthly or semi-monthly along its main stem; in recent years two samples per month have only been taken during summer months, but twice-monthly sampling occurred from March to October in early years of the dataset. Using these data to calculate temperature trends for a specific month can lead to aliased and biased estimates of long-term change, since the data are often limited by observational constraints including weather and scheduling. Model outputs offer the possibility of

May -October
November -April  filling in these spatio-temporal gaps and assessing the level of uncertainty due to these constraints. Here trends calculated using model output only at the times and locations for which data are available were compared with those associated with daily output at all grid cells in the mainstem (ModelRef(stn) and ModelRef(all), Figure 4; Figure S2). The standard errors of both estimates in some cases overlap each other, suggesting that the observational sampling resolution may be sufficient to capture long-term estuarine temperature change on a month-by-month basis. However, when aggregated by warmer and cooler periods of the year, the dissimilarities between these trends (ModelRef(stn) and ModelRef(all)) become more apparent, particularly in the warmer months from May to October (Table 3). This stretch of warmer months encompasses the seasonal shifts between states of cooler and warmer weather in May and October and is also the most active period for thunderstorms and the possibility of tropical storms (Lucy et al. 1979;Dolan et al. 1988). Water column temperatures can change rapidly after a significant summer wind-mixing event, especially during months when temperature stratification is likely. The availability of more frequent data would help reduce biases that may, for example, suggest that the water column is not warming as fast as it is in reality, perhaps due to postponing sampling until after strong wind events. In contrast, long-term changes in Bay temperatures computed from the model via ModelRef(all), are not biased by the potential sampling frequency issues described above, but have other sources of uncertainty ( Figure 3). For example, trends computed from Model-Ref(all) may be considerably different if the most recent years are particularly warm and wet than if they are relatively cool and dry. Temperature change computed from the sensitivity experiment assuming all four mechanisms (Combined) alleviates this issue by looking at long-term averaged trends that are less influenced by specific weather events. Overall, given the uncertainties in observational and modeled estimates of long-term estuarine temperature change, the fact that these estimates are all roughly equal to 0.02°C/year highlights the robustness of this result.
The estimates of long-term change in Bayaveraged temperature derived from this study are also in alignment with previous literature estimates. Preston (2004) calculated annual temperature trends of 0.016-0.021°C/year over a slightly earlier time period , whereas Ding and Elmore (2015) found surface trends that ranged from 0.005 to 0.175°C/year in the upper Chesapeake Bay over the period 1984-2010. The estimates from this study are slightly higher than some of these earlier estimates, likely because these are based on information encapsulating more recent years (1985-2019) that include recent periods of accelerated warming in the Mid-Atlantic region (Liao et al. 2015;Wallace et al. 2018). Our Chesapeake Bay warming estimates are also considerably greater than recent estimates of global sea surface temperature trends, which range between 0.007 and 0.012°C/year based on several large-scale reconstructed datasets (Hausfather et al. 2017).
Our work also demonstrates that although longterm warming does not differ substantially between the surface and bottom main stem, the long-term rate of warming waters throughout the Bay does vary both spatially and seasonally (Figures 5 and 6). Overall warming within the main stem Bay over the study period exhibits a clear seasonal cycle that is characterized by substantially greater-than-average warming from May to October and comparatively less warming or none at all from November to April (Table 3). This finding is in contrast with elevated rates of winter and spring warming reported by Preston (2004), but in agreement with those of Testa et al. (2018), who found significantly increasing surface temperature trends in the main stem Bay in June and July and insignificant trends during most of the remainder of the year.
More differentiated warming patterns are present in the upper tributaries and shoals adjacent to the main stem Bay ( Figure 5), but the majority of these areas follow the seasonal patterns present in air temperatures more closely, along with that of river warming (Figure 8). In agreement with Ding and Elmore (2015), this study finds comparatively greater increases in surface temperatures in the mainstem Bay than in tidal tributaries ( Figure 5): The greatest increase in the surface and bottom temperatures is typically found in the polyhaline region, whereas some of the smallest increases in mainstem temperatures are located in the northern mesohaline Bay ( Figure 5; Tables 4 and 5). This information may prove useful in identifying ideal locations for the deployment of continuous monitoring instruments designed to detect estuarine temperature change, which would also improve estimates of total Bay hypoxia (Bever et al. 2018).

Impacts of Atmospheric Change on Chesapeake Bay Temperatures
Changes in atmospheric temperatures and downwelling longwave radiation are highly correlated, occur in unison, and are fairly spatially homogeneous across the Bay (Figure 2a, 2b). Increasing air temperatures are associated with a positive feedback of absolute humidity given a constant relative humidity JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION (Soden and Held 2006), which can promote increased impacts of longwave radiation that are already amplified by rising concentrations of greenhouse gases. Air temperatures in the Chesapeake Bay region have increased faster (~0.025°C/year) than corresponding global averages (0.015-0.016°C/year over the period 1979Hartmann et al. 2013).
Atmospheric changes vary seasonally (Figure 2a; Figure S1), which can result in large differences in Bay temperature trends between months (Figure 7). On average, changes in air temperature and downwelling longwave radiation account for more than 80% of long-term estuarine bottom water warming; however, this percent varies by region (Figure 8). Atmospheric warming accounts for nearly all change in a given year and for the majority of the mainstem in the oligohaline and surface mesohaline regions of the Bay, but nearer the Bay mouth, this dominance decreases. Intriguingly, the lesser contribution of atmospheric warming to higher Bay temperatures is not solely due to the greater impact of Mid-Atlantic Bight warming in the summer months. In fact, the atmospheric forcing contributes a smaller absolute change in total Bay warming during July and August (Figure 7). This is partially due to the fact that the trends in atmospheric forcings above the Chesapeake Bay are weaker in July and August compared to April-June and September-October ( Figure S1). In addition, the gradient between water and atmospheric temperatures is lessened by the influence of ocean warming, which acts to decrease the ability of the Bay to act as a heat sink for warmer atmospheric conditions. This dynamic is strongest in the southern Bay, but is apparent even in the northern mesohaline region (Figure 7b, 7f).

Impacts of Sea Level Rise and Mid-Atlantic Bight Warming on Chesapeake Bay Temperatures
Chesapeake Bay temperatures are increasing at a rate (0.24 AE 0.15°C per decade; Table 3) that is more than twice as fast as the globally averaged temperature of the upper 75 m of the world's oceans (0.11 AE 0.02°C per decade; Rhein et al. 2013), but is in alignment with previous warming estimates along the continental shelf in the Mid-Atlantic Bight (0.26°C per decade, Forsyth et al. 2015;0.1-0.4°C per decade, Kavanaugh et al. 2017;0.43°C per decade, Wallace et al. 2018). Increases in offshore ocean temperatures used in the model experiments in this study (Figure 2c) are approximately twice the increases in atmospheric temperatures, and nearly quadruple the increases in globally averaged surface ocean temperatures. Atmospheric conditions driven by the position of the jet stream can produce short-term anomalies in shelf temperatures along the Mid-Atlantic Bight, but are likely not responsible for long-term warming patterns (Chen et al. 2014). An investigation of broader Northwest Atlantic Shelf warming over the time period of this study found that approximately two-thirds of observed warming may be attributable to multidecadal cycles, and that the remainder is driven by external forcings resulting from anthropogenic warming (Chen et al. 2020). Additional research has suggested that changes along the shelf in the Mid-Atlantic Bight are more similar to observed trends over Labrador, from which circulation is driven southwards toward the shelf adjacent to the Chesapeake Bay mouth (Shearman and Lentz 2010;Forsyth et al 2015;Wallace et al. 2018), and may be influenced by the greater intrusion of warm-core rings from the Gulf Stream (Gangopadhyay et al. 2019). These changes in offshore temperatures are central to the large regional variation in Bay warming, particularly with respect to the seasonality of summer bottom mesohaline and polyhaline anomalies. Seasonal changes in offshore warming patterns (i.e., increased rates of warming outside of summer months), could also result in greater impacts to bottom Bay waters during the winter and spring.
The mechanism driving this Mid-Atlantic Bight warming is inextricably linked with regionally elevated rates of sea level rise (Boon 2012;Sallenger et al. 2012). Numerous previous studies (Boon 2012;Ezer et al. 2013;Domingues et al. 2018;Ezer 2019) have documented changes with respect to the shifting position of the Gulf Stream that help to make the Mid-Atlantic Bight a global hotspot of sea level rise. Within the Bay, sea level rise causes the Bay to be cooler in the warmer months and warmer in the cooler months (Tables 4 and 5). Since sea level rise increases the total volume of the Chesapeake Bay, a possible explanation for this result is an increased thermal inertia, where the Bay would take longer to warm in the spring, and longer to cool in the fall (Tables 4 and  5; St-Laurent et al. 2019). However, this effect is very small compared to the atmospheric and oceanic warming mechanisms (Figure 8). Future sea level rise could amplify the observed summer cooling effect, but it will likely still only slightly offset the additional impacts of shelf warming. Future Bay bottom warming may be largely dependent upon additional changes in offshore currents, which are dependent upon both the rate of additional heat uptake by the oceans and changes to heat transport and large-scale circulation patterns along the Mid-Atlantic Bight (Saba et al. 2016;Thomas et al. 2017;Alexander et al. 2020).

Implications of Chesapeake Bay Warming
Previous studies have identified decreasing oxygen solubility as an important mechanism for future anticipated climate change effects (Li et al 2015;Irby et al. 2018;Ni et al. 2019). Summer warming in mesohaline bottom waters (1.01°C-1.16°C; Table 5) during the height of annually recurring hypoxia and anoxia has the effect of decreasing dissolved oxygen levels on the order of 0.1-0.2 mg/L based on the temperature dependence of oxygen solubility. Established management goals to improve levels of bottomdissolved oxygen (USEPA 2010) allow for an exceedance of water quality minimum standards (Irby and Friedrichs 2019) that such a change in solubility may produce. However, relatively spatially homogeneous warming in designated deep-water habitats will essentially lower the entire distribution of dissolved oxygen levels toward a threshold of water quality noncompliance. Other significant factors relating Bay oxygen levels to increased respiration rates and reduced nutrient loads entering the Bay must also be studied in order to determine the comprehensive cost of warming on Bay water quality and management efforts. Better quantification of this climate change component of temperature variability could provide insights into how the success of nutrient reductions Zhang et al. 2018) may have been limited, and the consequent economic costs of climate change on Bay restoration efforts.
Future estuarine warming has the potential to affect other temperature thresholds critical to endemic flora and fauna in the Chesapeake Bay. Shifts in the range of native seagrasses have been documented in the polyhaline region of the Bay and temperature change has been identified as the driver behind changes in distributional patterns (Richardson et al. 2018). Increasing temperatures can produce cascades of varying effects on aquatic organisms, including those that impact copepod mortality (Pierson et al. 2016), fish recruitment and timing based on those same zooplankton populations (Millette et al. 2020), metabolic tolerance of increased temperatures for demersal species (Marcek et al. 2019), and Bay habitat suitability for juvenile sandbar sharks (Crear et al. 2020). These examples only scratch the surface of additional identified impacts of temperature on marine organisms not listed here; more are summarized in Najjar et al. (2010).
Although the role of anthropogenic emissions in global warming over the course of the last century is well understood (Bindoff et al. 2013), it is not certain to what extent the observed warming signal in the Chesapeake Bay over the past 30 years is controlled by global and regional processes, natural variability arising from multidecadal oscillations, or other unknown factors. However, a comparison of global climate model simulations with observations in the northeastern U.S. indicates that anthropogenic factors are dominant (Kunkel et al. 2013).
Understanding the role that natural variability has had on observed warming over the past 30 years can also help determine how the future of projected warming may be exacerbated by multidecadal processes. The length of time necessary to discern a clear climate change signal in such a highly variable coastal environment is also a critical factor when differentiating global signals from interannual fluctuations.
Based on the results presented here, the likely drivers of future temperature change in the Bay can be better identified and total future temperature change can be more accurately projected. The future trajectory of warming in the Chesapeake Bay may be strongly influenced by the changes that occur offshore. Changes to both river inputs and Bay waters are likely to be primarily controlled by changes in atmospheric forcing that are predicted to continue increasing, but the response of ocean currents and their exchange with the Bay can be an important influence as well.

CONCLUSIONS
This study analyzed the mechanisms causing changes in Chesapeake Bay temperatures over the period 1985-2019, including atmospheric forcings, increasing ocean temperatures, sea level rise, and riverine warming. Despite inherent uncertainties in estuarine temperature trend estimates computed from monthly/semi-monthly data, regionally aggregated observed and modeled temperatures produced consistent seasonally varying patterns of Bay warming. Sensitivity experiments simulating impacts of observed warming mechanisms produced changes in temperature (0.02°C/year) that have high explanatory power of observed trends, with elevated rates in warmer months of the year (0.04°C/year) and lower rates in cooler months (0.01°C/year). Changes in atmospheric conditions, including both increasing air temperatures and downwelling longwave radiation, have driven much of the monthly variability of Bay temperature trends, except in the southern Bay where additional summer warming has resulted from warmer waters advecting into the Bay from the continental shelf. Sea level rise has effected a slight cooling trend during spring and summer months that is greater than the inconsequential impact of river warming; however, this cooling effect is insufficient to counteract the impacts of atmospheric and oceanic warming.
Based on these conclusions, future warming patterns in the Chesapeake Bay are likely to be determined by the evolution of combined increases in atmospheric warming and longwave radiation as well as continental shelf circulation patterns. Studies of future warming in the Chesapeake Bay must thus consider not only impacts of atmospheric change, but also change in the adjacent coastal ocean. If future climate projections express high variability in continental shelf circulation and temperature patterns, it may be difficult to constrain future increases in estuarine temperatures. As summer cooling induced by sea level rise accelerates in the coming decades, this impact on estuarine temperature may become more important and should not be overlooked. Uncertainties in climate projections related to these primary mechanisms underlying historical Bay warming should inform the manner in which future studies frame our confidence in future physical, biogeochemical, and ecological outcomes.

SUPPORTING INFORMATION
Additional supporting information may be found online under the Supporting Information tab for this article: Includes additional regional details for water temperature trends and monthly atmospheric forcing estimates. ACKNOWLEDGMENTS This paper is the result of research funded by the National Oceanic and Atmospheric Administration's National Centers for Coastal Ocean Science under award NA16NOS4780207 to the Virginia Institute of Marine Science. The authors acknowledge William & Mary Research Computing for providing computational resources and/or technical support that have contributed to the results reported within this paper (https://www.wm.edu/it/rc). CBP Watershed Model output was provided by Gopal Bhatt and others at the Environmental Protection Agency's Chesapeake Bay Program Office. Thank you to Dr. Grace Chiu for statistical consultation provided through the VIMS ESTDatS lab. We also thank the two anonymous reviewers whose comments and suggestions significantly improved the paper. The model results used in the manuscript are permanently archived at the W&M ScholarWorks data repository associated with this article and are available for free download (https://doi.org/10.25773/c774-a366