Impacts of Water Clarity Variability on Temperature and Biogeochemistry in the Chesapeake Bay

Estuarine water clarity depends on the concentrations of aquatic constituents, such as colored dissolved organic matter, phytoplankton, inorganic suspended solids, and detritus, which are influenced by variations in riverine inputs. These constituents directly affect temperature because when water is opaque, sunlight heats a shallower layer of the water compared to when it is clear. Despite the importance of accurately predicting temperature variability, many numerical modeling studies do not adequately account for this key process. In this study, we quantify the effect of water clarity on heating by comparing two simulations of a hydrodynamic-biogeochemical model of the Chesapeake Bay for the years 2001–2005, in which (1) water clarity is constant in space and time for the computation of solar heating, compared to (2) a simulation where water clarity varies with modeled concentrations of light-attenuating materials. In the variable water clarity simulation, the water is more opaque, particularly in the northern region of the Bay. This decrease in water clarity reduces the total heat, phytoplankton, and nitrate throughout the Bay. During the spring and summer months, surface temperatures in the northern Bay are warmer by 0.1 °C and bottom temperatures are colder by 0.2 °C in the variable light attenuation simulation. Warmer surface temperatures encourage phytoplankton growth and nutrient uptake near the head of the Bay, and fewer nutrients are transported downstream. These impacts are greater during higher river flow years, when differences in temperature, nutrients, phytoplankton, and zooplankton extend further seaward compared to other years. This study demonstrates the consequences of utilizing different light calculations for estuarine heating and biogeochemistry.


Introduction
Rivers flowing through watersheds collect sediment, colored dissolved organic material (CDOM), and nutrients from the surrounding landscape. These materials eventually reach estuaries and coastal ocean environments where they affect water clarity. CDOM is a major light absorber in the energetic blue region of the visible spectrum, whereas sediment contributes to light scattering in the water column (Mobley 1994). Nutrients affect water clarity indirectly by fueling phytoplankton production of organic matter, which attenuates light near the water's surface through light absorption by pigments as well as light scattering. When water is clear, sunlight heats water at greater depths compared to when the water is turbid. Since sunlight is the main source of heat in aquatic environments, variations in water clarity influence water temperature.
Temperature range and variability have wide-ranging implications for estuarine biogeochemistry and ecology. Many microbial processes are regulated by temperature (Apple et al. 2006;Lomas et al. 2002), which in turn impacts primary productivity and nutrient cycling. The timing of temperature shifts is also important for phytoplankton phenology ) and the initiation of harmful algal blooms (Anderson et al. 2012) that can dominate the microbial community under specific environmental conditions. Seasonal transitions are critical for larger animals with complex life cycles, such as fish and shellfish that utilize estuarine waters as feeding, nursery, or wintering habitat. While species' distributions may roughly follow temperature distributions based on their optimal growth ranges and thermal tolerances, episodic exposure to extreme hot or cold temperatures can lead to significant mortalities (Moore and Jarvis J. C. 2008;Bauer and Miller 2010). Studies investigating long-term high-frequency temperature records of species' body temperature emphasize the importance of examining environmental variability at temporal and spatial scales relevant to organisms for understanding population changes (Helmuth et al. 2006).
Despite the importance of accurately predicting temperatures, many estuarine numerical models are ill-equipped to study the interaction between shortwave heating and water clarity. In many coupled hydrodynamic-biogeochemical models, the hydrodynamic component of the model simulates physical processes such as heating and circulation. The resulting temperature, salinity, and other environmental parameters can be used in calculations of biogeochemical processes, resulting in spatially and temporally varying simulations of light-attenuating materials. The influence of these materials on water clarity is unaccounted for in the model's physical calculations in this typical "one-way" coupled configuration (Mobley and Boss 2012). Instead, the thermodynamics of in-water solar heating are usually prescribed in a way that is constant in space and time by assigning a characteristic Jerlov water type, a classification system for natural waters based on the water body's optical attenuation depth (Jerlov 1968). For example, this is the default configuration for the widely used Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams J. C. 2005). Moreover, it is difficult to assess the prevalence of this constant light attenuation assumption because the details for solar heating in the hydrodynamic routine are often not reported in model documentation publications, even in studies with very complex optics for ecosystem calculations (del Barrio et al. 2014). In addition, most estuarine water quality models are run with offline coupling (Cerco et al. 2010;Testa et al. 2014) prohibiting interaction between the hydrodynamic and biogeochemical calculations.
By contrast, in a "two-way" coupled model, the spatiotemporally varying water clarity parameters from the biogeochemical calculations are used in calculations of in-water solar heating. In a numerical modeling study of the Hudson River plume, accounting for variable light attenuation by phytoplankton resulted in warmer surface temperatures and intensified temperature stratification compared to a constant light attenuation simulation (Cahill et al. 2008). These effects were also reported in a similar study of Monterey Bay, CA (Jolliff and Smith 2014), where the enhanced temperature stratification stimulated a 27% increase in modeled phytoplankton biomass for the variable light attenuation simulation. Under intensified stratification, which coincides with nutrient replete conditions, phytoplankton growth is restricted to a near-surface layer, enabling sustained exposure to light and improving nutrient utilization for growth. While it is difficult to find observational analogs for these twin modeling experiments (i.e., simulations with and without variable light attenuation), observations of concurrent changes in water clarity and temperature provide further evidence that high abundances of light-attenuating materials affect heating in coastal environments (Fournier et al. 2017).
To further investigate the relationship between water clarity and temperature, we perform a case study in the Chesapeake Bay, comparing simulations of a coupled hydrodynamic-biogeochemical model with and without variable light attenuation in the solar heating calculations. We perform this comparison to quantify the effect of twoway coupling between hydrodynamic and biogeochemical calculations. In the following two sections, we introduce  , the river forcing for simulations in this study the Chesapeake Bay estuary and the numerical model that we used in this study. In the "Results and Discussion" sections, we present how seasonal, year-to-year, and episodic variations in water clarity affect temperature and temperature-dependent biogeochemical processes. Guided by these results, we infer how long-term changes in water clarity might affect the Bay's hydrography and ecosystem.

The Chesapeake Bay
The Chesapeake Bay experiences considerable seasonal variability in physical forcing. Surface temperatures in the Bay are near zero during the coldest months of the year (January to February) and can exceed 30 • C in the warmest months (July to August). During the late winter and early spring months, melting snow and ice in the Chesapeake Bay watershed contribute large amounts of water through the landscape and into the estuary (Fig. 1). Freshwater plumes from rivers flow into the Bay's main channel carrying sediment, CDOM, and nutrients (Schubel and Pritchard 1986). The dynamics of these plumes determines the extent and timing of light attenuation impacts on heating and biogeochemistry. The spring freshet is associated with the annual onset of the heating and stratification of Bay waters, as buoyant freshwater flowing over saltier estuarine waters is warmed by sunlight (Murphy et al. 2011). This stratification persists through the summer until the fall, when the Bay is less stratified due to wind forcing and cooling air temperatures. These processes decrease water column stability, particularly during fall storms when strong winds impart the mixing energy required to destratify the water column (Goodrich et al. 1987).
The annual cycle of phytoplankton abundance in the Bay is strongly influenced by variations in nutrient inputs and light availability. Nitrate concentrations are highest in the tributaries and in the upper Bay. Nitrate is consumed by phytoplankton or transformed through nutrient cycling as it is transported from the riverine inputs to the mouth of the Bay. During the late winter to spring months, substantial riverine nutrient inputs and longer days provide the necessary ingredients for phytoplankton growth. The productive season continues until the late summer and fall months, when shorter days, less nutrient inputs, and increased grazing from zooplankton suppress phytoplankton growth. Annual flow conditions have been linked to the spatial extent and magnitude of phytoplankton biomass in the Bay, owing to the associated magnitude and distribution of nutrient concentrations. High-flow "wet" years exhibit higher phytoplankton abundance extending further down the Bay (seaward) compared to dry years (Harding et al. 2016).
Water clarity in the Chesapeake Bay has declined in recent decades. In the mesohaline region, shallower Secchi depths from 1965-2015 are correlated with increasing chlorophyll-a concentration (Harding et al. 2019). However, this pattern does not generally hold for the other regions. In the oligohaline region, water clarity trends are opposite the chlorophyll-a trends. Increases in chlorophyll-a are accompanied by improvements in water clarity (deeper Secchi depth), while declines in chlorophyll-a are coincident with shallower Secchi depths. Considering the importance of CDOM in determining water clarity in the Chesapeake Bay, especially in the oligohaline zone, these Secchi depth trends may be indicating the seasonal trends in CDOM Fig. 2 a Map of Chesapeake Bay with locations of stations referenced in this study. Station markers are colored by salinity zone, from oligohaline (low salinity) to mesohaline (moderately salty) to polyhaline (salty). b Transect view of the Chesapeake Bay, following the blue line in panel a absorption over the past several decades. In a statistical analysis of water clarity measurements, Testa et al. (2019) found CDOM to be a key driver of water clarity spatial variability in the Chesapeake Bay, while total suspended solids and chlorophyll-a explained the temporal variability in different regions of the Bay.
The Chesapeake Bay water temperature record, though highly variable, indicates Bay-wide warming in recent decades. In the 1990s, water temperature was about 1 • C warmer than in the 1960s (Najjar et al. 2010). The heating changes associated with these shifts in water clarity likely contribute to the observed variability and trend in the historic temperature record. However, the contribution to changes in heating by water clarity remains understudied. Numerous modeling efforts of the Chesapeake Bay seek to accurately predict the estuary's hydrography and biogoechemistry (see Irby et al. (2016) for a recent review), but none employs a "two-way" coupling that allows variations in water clarity to affect solar heating.
In this study, we report results primarily from the stations along the main stem of the Bay (shown in Fig. 2), which correspond to a subset of sites where the Chesapeake Bay Program conducts repeat water quality monitoring measurements. Cruise data from these and other sites were used in evaluation for the model used in this study, ChesROMS-ECB (Feng et al. 2015;Irby et al. 2016;Da et al. 2018) and are synthesized in other studies of the Chesapeake Bay. The northernmost stations are highly influenced by the Susquehanna River, which is the largest source of freshwater to the Bay. Further downstream, salinity generally increases until it reaches its highest level at the mouth, where Bay waters mix with the salty Atlantic Ocean. Below, we describe the longitudinal gradient along the Bay in terms of salinity, with fresher waters in the oligohaline stations (CB2.1, CB2.2, CB3.1), the upper (CB3.2, CB3.3C, CB4.1C, CB4.2C, CB4.3C) and lower mesohaline (CB5.1, CB5.2, CB5.3, CB5.4, CB5.5) stations in the middle salinity zones, and the saltier polyhaline stations (CB6.1, CB6.2, CB6.3, CB7.3, CB7.4) closest to the Bay mouth.

Methods
ChesROMS-ECB is a modeling system that combines hydrodynamics from the Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams J. C. 2005) with an estuarine-carbon-biogeochemistry (ECB) model of the Chesapeake Bay Feng et al. 2015). The ECB and ROMS are coupled at every time step, which is 1 min in model time. Atmospheric forcing variables (e.g., winds, solar radiation, air temperature) are prescribed based on the North American Regional Reanalysis (Mesinger et al. 2006). Daily freshwater, nitrogen, and carbon inputs are derived from the Dynamic Land Ecosystem Model (DLEM), a terrestrial biogeochemical-hydrological model previously used to study North American East Coast hydrology , nitrogen fluxes , and carbon fluxes . Sediment inputs are obtained from the Chesapeake Bay Program's watershed model (Shenk and Linker 2013). In our model configuration, terrestrial nutrients and sediment enter the Bay at the sites of the 10 largest rivers.
The incident solar radiation at the surface of the water is propagated through the water column in two separate calculations. The ROMS hydrodynamic model calculates the in-water light profile to determine the contribution of solar radiation to water heating. The ECB model calculates a different light profile for use in biogeochemical processes, such as phytoplankton growth. The main consequence of this inconsistency in light calculations is that variations in light-attenuating materials are not accounted for in the heating calculations.
In the ROMS hydrodynamic model, the vertical profile of light in water follows the formulation of Paulson and Simpson (1977), where the downwelling irradiance I (W m −2 ) at depth z (defined positive upward) is given by where I (0) (W m −2 ) is the incident downwelling irradiance at the surface less reflected solar radiation and R, ζ 1 [m], ζ 2 [m] are the best fit parameters to measure irradiance profiles for wavelengths 400 to 1000 nm. The parameter R is a fraction of the energy associated with the attenuation coefficient, (1/ζ 1 ), and the remainder of the energy, In the default configuration of the ROMS, R, ζ 1 , and ζ 2 are assigned to Jerlov Type I, representing clear ocean waters. The user is tasked with assigning the water type parameter to one of 9 Jerlov water types. The attenuation depth ζ 1 represents the attenuation length scale for nearinfrared portion of the spectrum, and ζ 2 is the attenuation depth for the visible wavelengths. The attenuation length scale is the distance where the intensity of incident irradiance is diminished to (1/e) of the initial value, or in other words, the depth where about 37% of the incident light is remaining. These parameters apply to the entire model domain and are fixed for the duration of the simulation.
Our biogeochemical model, ECB, calculates photosynthetically available radiation (PAR; including wavelengths 400 to 700 nm) for use in the biogeochemical formulations. PAR decreases approximately exponentially with depth: where I (0) (W m −2 ) is the irradiance at the water surface for the visible and near-infrared wavelengths (400 to 1000nm), and P ARf rac = 0.43 is the fraction of solar radiation that is assumed to be used in photosynthesis (400 to 700 nm) (Fasham et al. 1990). The attenuation coefficient, k bgc (m −1 ), varies as a function of the modeled particulate and dissolved materials in the water adapted from (Xu and Chao 2005;Xu and Hood 2006), with revised coefficients to match observations during our modeling study period (Feng et al. 2015): where [T SS] is the concentration of total suspended solids, and salinity, S, is a proxy for colored dissolved organic matter. The resuspension of inorganic sediment is represented implicitly through the imposition of a lower bound (0.6 m −1 ). It will be shown later that this lower bound is only activated in lower parts of the Bay. When k bgc > 0.6 m −1 , it is calculated by the parameterized expression in Eq. 3 which can be rewritten as the sum of each modeled attenuator: The attenuation by pure water is 0.04 m −1 (Fasham et al. 1990), k CDOM = 1.36 − 0.057S, and the contribution by the TSS term is separated into its component parts; i.e., inorganic suspended solids (ISS), detrital carbon (detC), phytoplankton (phyto), and zooplankton (zoo). We attribute the constant 1.4 from Eq. 3, minus the attenuation by pure water, to the CDOM term based on the interpretation by (Xu and Chao 2005) of this constant as the "attenuation due to pure water and CDOM in freshwater." The negative sign of the salinity term describes the inverse relationship between salinity and light attenuation by CDOM. CDOM in the estuarine environment is more abundant in fresher water because it mostly originates from terrestrial organic matter. This spatially and temporally varying light calculation more realistically represents underwater light in the natural environment, compared to the constant light attenuation scheme used in the modeled hydrodynamics described above.
In order to investigate the effects of variable light attenuation on solar heating, we changed the light profile in the hydrodynamic calculations to match the light profile in the biogeochemical calculations (Fig. 3, "Experimental Run: bgc-light"). We assigned the PAR from the ECB model (Eq. 2) as the solar energy for wavelengths 400-700 nm in the ROMS hydrodynamic model (second term in Eq. 1). The remaining energy from wavelengths 700-1000 nm is attenuated over a very shallow depth due to the rapid light attenuation for these wavelengths by water and particles (Jerlov 1968;Doxaran et al. 2007). For comparison, we also ran a reference simulation, "standard" run, where the light calculations for heating follow the standard approximation in the default ROMS configuration. We assigned the parameters in Eq. 1 to values for Jerlov Type 7, which is most consistent with in situ measurements of light attenuation for the Chesapeake Bay (Xu and Chao 2005). For the remainder of this manuscript, we assign the term k h to refer to the light attenuation coefficient for calculating shortwave heating in the hydrodynamic calculations and k bgc to refer to the light attenuation coefficient for PAR in the biogeochemical calculations. For the "bgc-light" run, k h = k bgc because the irradiance profile from the biogeochemical calculation is used in the solar heating calculation. For the "standard" run, k h =0.66 m −1 according to the coastal Jerlov water type used here. Fig. 3 Schematic of experimental and control runs. In the control run, shortwave (SW) solar heating in the hydrodynamic calculations affects water temperature, which is used in the biogeochemical calculations to determine the concentrations of light-attenuating particles (e.g., phytoplankton). In the experimental run "bgc-light," the irradiance calculated in the biogeochemistry model is then used in calculations for solar heating in the "bgc-light" run, establishing a two-way coupling that allows the biogeochemistry to affect the hydrodynamics at every time step The changes in heating between the two model runs affect biogeochemical processes according to temperaturedependent parameterizations, further detailed in Irby et al. (2018) and Da et al. (2018). The phytoplankton growth rate and zooplankton grazing rate parameterizations in ChesROMS-ECB are temperature-dependent, based on experimental observations of the temperature dependence of Chesapeake Bay microbial processes (Lomas et al. 2002). The maximum phytoplankton growth rate increases exponentially for temperatures above 20 • C, which suggests that small changes in temperature can correspond to large changes in phytoplankton growth under optimal nutrient and light conditions. The maximum zooplankton grazing rate is also parameterized as an exponentially increasing function of temperature.
Model dynamics were spun up for 1 year, and then simulations were run for 2001-2005 to simulate variable environmental conditions. Riverine fluxes, oceanic conditions on the continental shelf, and atmospheric forcings (e.g., winds, precipitation, atmospheric pressure, relative humidity, surface air temperature, shortwave, and longwave radiation at the sea surface) are the same in the two model runs. Based on the DLEM modeled freshwater discharge  of all Chesapeake Bay rivers from 1980-2014: 2001-2002 were dry years with total river discharge less than the 25th percentile of all years; 2003-2004 were wet years with river discharge greater than the 75th percentile; and 2005 was very close to an average year. For the remainder of this manuscript, we report phytoplankton and zooplankton in terms of carbon, after having converted our model output in nitrogen units to carbon using the Redfield ratio. Additionally, we define surface results as values averaged from the free surface to 1-m depth.

Light Attenuation
In the "bgc-light" run, k h values are generally larger than 0.66 m −1 , indicating lower water clarity than Jerlov Type 7. This results in solar heating at deeper depths in the "standard" run compared to the "bgc-light" simulation. Differences in light attenuation for heating calculations in the two runs is the difference between k h in the "bgc-light" simulation minus 0.66 m −1 , the attenuation coefficient for the "standard" simulation ( Fig. 4). For mainstem stations in the oligohaline region, the difference in surface k h ranges from 0.8 to 1.8 m −1 , which translates to a difference in the attenuation depth of −0.8m to −1.1m (k −1 h "bgc-light" minus k −1 h "standard"). Throughout the year, differences in the monthly averaged attenuation depth range from 0.4 to 0.8 m shallower in the "bgc-light" run for the upper mesohaline stations, and 0.1 to 0.5 m shallower for the lower mesohaline stations. The annual average modeled attenuation coefficient for the polyhaline region is 0.66 m −1 , resulting in small differences between the two runs.

Heating and Temperature
Here, we report changes in the heat content of the Bay (H ), which we define as the volume integration of water temperature (T [ • C]) over the Bay: where the reference density is ρ 0 = 1020 kg m −3 , specific heat is c p = 4000 J kg −1 • C −1 , and V is the surfaceto-bottom volume of the Bay, including the tributaries but excluding the continental shelf. Since water temperatures are generally warmer than 0 • C, heat content is a positive quantity everywhere. The heat content is lowest in the winter, increases to a peak in the summer months and then declines again during the fall months (Fig. 5a). The heat content of the Bay is smaller in the "bgc-light" run from roughly April through September every year. The heat contents of the two model runs are similar from October through April.
In the "bgc-light" run, monthly average (2001)(2002)(2003)(2004)(2005) surface temperatures are generally warmer late-winter to springtime, and colder from summer to autumn at the mainstem stations compared to the "standard" simulation ( Fig. 6a). Surface temperatures at stations in the upper mesohaline region station are warmer by 0.07 • C in May, and colder by 0.02 • C in September. While the mean temperature change in June in this region is 0.06 • C, it is accompanied by a 0.2 • C increase in the temperature amplitude, with an average increase in 0.18 • C in the maximum temperature and 0.02 • C decrease in the minimum temperature. Average surface temperature differences in other regions are smaller in magnitude but exhibit similar seasonality.
While daytime surface temperatures are generally warmer in the "bgc-light" run, decreased water clarity shades deeper waters from solar heating, which decreases  subsurface temperatures. The monthly average bottom temperature difference reaches −0.2 • C in July for the upper mesohaline stations, with small changes during the wintertime (Fig. 6e). The surface-to-bottom temperature difference increases by approximately 0.2 • C from May to August in this region, mostly due to the decrease in bottom temperatures (Fig. 6f). This increase in temperature stratification does not have a considerable impact on the density stratification, due to the dominant contribution of salinity in determining density in the estuary.

Biogeochemistry
The total integrated phytoplankton in the Bay is lowest in the winter, increases to a peak in the summer months, and then declines again during the fall months (Fig. 5b). In the "bgc-light" run, there is less total phytoplankton in the Bay throughout the spring to summer months. When integrated over different salinity ranges, the decreases in phytoplankton over the 5-25 psu salinity range are larger than the small increases in phytoplankton between 0 and 5 psu (Appendix, Fig. 15). Differences in nitrate concentration integrated over the 0-5 psu salinity range account for most of the change in nitrate during the summer months. There are coincident increases in total phytoplankton nitrate uptake in this salinity range during these times. For most months, there is less nitrate uptake in the 5-25 psu salinity range.
During the productive months, warmer surface temperatures in the oligohaline stations are accompanied by higher phytoplankton concentrations in the "bgc-light" run compared to the standard run. Relative differences for average surface concentrations range from + 1-2% from May to October (Fig. 7d). This small increase in phytoplankton is accompanied by −3 to −8% changes in surface nitrate concentrations for the upper mesohaline stations from May to September (Fig. 7e). In other zones, the relative change in average surface phytoplankton for 2001-2005 is small for most months (< 1%). There is less surface zooplankton in the "bgc-light" run compared to the standard run for all months of the year. The largest relative change in surface zooplankton is a 9% decline averaged over mainstem stations in August (Fig. 7f), which is one month before the modeled peak in zooplankton concentrations in September (Fig. 7c).

Impacts of Variable Light Attenuation on Interannual Variability of Summer Temperature and Biogeochemistry
The spatial extent of differences in summertime (June to August) k h varies from year to year, with larger values throughout the Bay in 2003 in the "bgc-light" run, in contrast to 2001 when the difference is mostly confined to oligo-and upper mesohaline regions (Fig. 8a-c). The spatial extent of colder summertime bottom temperatures is also larger in 2003 compared to 2001 and the average of all the other years ( Fig. 8d-f). Declines in water temperature are present within the water column (Appendix, Fig. 16), with the largest differences occurring around 5-m depth.
The year-to-year variations in phytoplankton and zooplankton differences between the two model runs follow similar patterns as the light attenuation differences. In 2001 and 2005, increases in phytoplankton are confined to the northern reaches of the Bay, with decreases in nutrients near zero by 250 km from the mouth of the Bay, near station 3.2 (Fig. 9c). For 2003, the increases in phytoplankton extend much further downstream, and the change in nutrients is found as far as nearly 150 km from the bay mouth, near

Impacts of Variable Light Attenuation During a Spring Freshet and Fall Storms
The most intense spring freshet during our study period occurred during the first week of April 2005 when a large volume of freshwater flowed into the Bay. The peak daily discharge from the Susquehanna River, which is the largest source of freshwater into the Bay, reached 1.15 × 10 4 m 3 s −1 on April 5, 2005 (Fig. 1). The water was nearly opaque in the upper Bay between April 4 and 7, with an average light attenuation length scale (k −1 h ) of 0.2 m for oligo and upper mesohaline stations in the "bgc-light" run. As the ISS settled out of surface waters, attenuation depths increased to 0.7 m by April 18-21. This event marked a warming period for surface waters, with modeled temperatures for all mainstem stations starting at 10.5 • C early on the 18th and reaching 14.6 • C by afternoon on April 20 for the "standard" run.
The largest Bay-wide hourly temperature differences between the two simulations occurred during this April 18-21 event. Differences in daily maximum temperatures at Fig. 9 Absolute difference in average June through August surface a phytoplankton c nitrate and e zooplankton concentrations along the Bay mainstem. Relative difference in surface b phytoplankton d nitrate and f zooplankton ("bgc-light" minus "standard"). Mainstem locations shown in Fig. 2 individual stations ranged from 0.1 • C to 1.2 • C, and bottom minimum temperatures were colder by up to 0.4 • C in the upper Bay (Fig. 10). During these days, the hourly surface temperatures averaged over the oligohaline stations were up to 0.4 • C warmer in the "bgc-light" run, and up to 0.6 • C warmer for the mesohaline stations.
The largest freshwater influx for our simulation period occurs in the fall of 2004, with elevated daily freshwater Susquehanna River input starting in mid-September, reaching 1.40 × 10 4 m 3 s −1 on Sept 20 (Fig. 1). From September 18 to the 21, surface temperatures cooled about 2 • C, as the average surface temperature for all mainstem stations decreased from an average of 25.5 • C to 23.2 • C for the weeks before and after this period in the "standard" run.
The average surface temperature difference between the two simulations is + 0.01 • C for all mainstem stations during the week following the peak in freshwater discharge, September 21 to 27. This small change in the average is accompanied by an increase in the diurnal temperature amplitude. This is most apparent in the upper mesohaline stations, with an increase of 0.2 • C in the daily afternoon temperature maximum and a decrease of To further examine the effect of changing water clarity on heating during these events, we compare the dynamics at station CB3.2 during the spring and fall (Fig. 11). At station 3.2 from April 18 to 21, the temperature difference between the surface to 5 m increased by 0.6 • C, with a +0.3 • C difference in temperature for the upper 1 m and −0.3 • C difference at 5 m between the two model runs. During these days, there was an averaged 5% more surface phytoplankton in the "bgc-light" run compared to the "standard" run. The largest relative difference was an +18% change in the afternoon of the 18th, as surface ISS concentrations settled out from the morning into the afternoon.
During the fall event, from September 21 to 27, there are instances of increased temperature stratification during this period (Fig. 11b), but it is not sustained as in the springtime example. In general, temperature stratification is negative during these days, as cold, freshwater flows above warm, salty water, except in the afternoons. At station 3.2, surface temperatures of the run "bgc-light" are warmer than the "standard" by up to 0.3 • C during the afternoon hours. The increase in temperatures is confined to a shallow layer near the surface, with negative temperature differences below, increasing the temperature difference between the surface to 5 m by up to 0.5 • C during those hours. Surface phytoplankton concentrations are greater by up to 15% as the ISS settles out of the water.

Discussion
Seasonal variations in water clarity in different regions of the Bay are influenced by the seasonal dynamics of the aquatic constituents. In the oligohaline region, large pulses of sediment are coincident with low water clarity events (i.e., high light attenuation) during the spring freshet. While light attenuation by ISS makes up 30% of the total k bgc in April at the oligohaline mainstem stations, the largest contributor is CDOM for the years 2001-2005 (Appendix, Fig. 12). During the summertime, the relative contribution by CDOM is larger, at an average contribution of about 70% from June to August, with phytoplankton contributing approximately 10%. Heavy rain during autumn and winter storms again increases delivery of sediments to the Bay, which increases the contribution by ISS during those months. In the meso-and polyhaline regions, the relative contribution by phytoplankton is larger while the contribution by ISS is smaller compared to the oligohaline stations. The k bgc calculated by the parameterized expression in Eq. 3 is sometimes less than 0.6 m −1 in the meso-and polyhaline regions. This occurs at times when salinity, the modeled proxy for CDOM, in the polyhaline region is high and the phytoplankton concentration is low. During times when k bgc > 0.6 m −1 , modeled contributors to light attenuation are CDOM, phytoplankton, and detrital carbon, listed in order of decreasing importance. In situ measurements of light absorption in the Bay confirm the importance of CDOM on optical properties through all regions of the Bay. The light absorption coefficient for CDOM is highest in the oligohaline regions, decreasing linearly along a conservative mixing line from low to high salinity waters (Rochelle-Newall and Fisher 2002), with evidence of seasonal variability extending into the shelf waters of the Mid Atlantic Bight (Del Vecchio and Blough 2004).
These processes contribute to the differences in water clarity between the two model runs. The spring and fall sediment inputs contribute to peaks in k h during April and September in the oligohaline region for the "bgc-light" run. A seasonal cycle for the difference in k −1 h is most apparent in the mesohaline region, with a peak in k h for the "bgclight" run in April and May and minima in the winter months (Fig. 4). During the summertime, average river flow is lowest in 2001 and highest in 2003 (Fig. 1), which affects the spatial extent and magnitude of differences in k h between the two model runs. The most widespread differences in the penetration depth of solar heating occur during the high-flow year of 2003 (Appendix, Fig. 13). Interannual changes in river flow have previously been linked to phytoplankton dynamics in the Bay (Harding et al. 2016), with increased nutrient loadings during high-flow years being associated with increased chl-a and shift in flora to larger cells. In addition to increased light attenuation by phytoplankton, the modeled contribution by CDOM is also greater when there is more freshwater, further decreasing water clarity.
The increased light attenuation of the "bgc-light" simulation warms a thin surface layer while decreasing solar heating for deeper waters. Because the warmer surface layer is a small fraction of the volume of the Bay, most of the Bay is generally colder, and the heat content is smaller from April to September (Fig. 5). Since the two model runs were forced with the same atmospheric conditions, river fluxes, and oceanic boundary conditions, differences in heat content must originate from changes in surface heat fluxes at the air-water interface or horizontal heat fluxes (at the mouth of the Bay and from the rivers). Surface heat fluxes include the evaporative or latent heat flux and sensible heat flux, the conductive heat flux from the water to the air. Monthly climatology of the heat fluxes for our simulation years shows latent heat flux is the largest heat sink in the Bay throughout the year (Appendix, Fig. 14). Latent heat flux is more negative in the "bgc-light" run, due to the higher surface temperatures of this run. Latent heat flux acts as a larger heat sink from March to June, causing the relative decrease in heat content during these months. From July to September, the latent, sensible, and horizontal fluxes are less negative in the "bgc-light" run, narrowing the gap in heat content between the two simulations during these months. In the natural environment, changes in surface water heat fluxes affect the ambient air temperature. While our model simulations are forced with prescribed atmospheric conditions, Jolliff and Smith (2014) demonstrate how changes in water clarity affect air temperature in a coupled ocean-atmosphere system. For fully coupled modeling systems, water clarity variability is important for modeling air temperature variability as well.
Summertime flow is generally a good predictor of the extent and intensity of the change in subsurface temperatures, except in 2005 which is the 2nd driest summer but exhibits similar subsurface temperature changes as the wettest summers, 2003 and 2004. The largest spring freshet occurred in 2005, which suggests springtime flow conditions and associated changes in water clarity are factoring in to the intensity and extent of summertime subsurface temperature differences between the two simulations. Previous studies have demonstrated links between winter to spring freshwater flow and summertime stratification, e.g., (Murphy et al. 2011), which supports this hypothesis. In our model runs, colder temperatures in the "bgc-light" run reach further seaward along the mainstem during years with higher flow (Appendix, Fig. 16). Due to the exponential dependence of both phytoplankton growth and zooplankton grazing on temperature, colder temperatures affect zooplankton concentrations both directly (by changing temperature) and indirectly (by decreasing zooplankton food availability).
In the "bgc-light" run, increased phytoplankton nutrient uptake in low salinity waters decreases the amount of nitrate that is transported downstream to higher salinity waters. This is supported by our Bay-integrated and surface results. Increases in phytoplankton in fresher waters are concurrent with declines in nitrate in saltier waters. Because fewer nutrients are transported downstream, there is less nitrate uptake and phytoplankton in the Bay. The colder temperatures and fewer phytoplankton of the "bgc-light" run contribute to the decline in zooplankton during the late summer to early fall months. Summertime flow conditions determine the magnitude and seaward extent of surface phytoplankton concentrations and the associated changes in nutrients. Increases in phytoplankton are confined to the northern Bay during years with the lowest June to August flow (e.g., 2001 and 2005), while the increases in phytoplankton extent are found further downstream in high flow years (e.g., 2003). These differences in the spatial distribution of the phytoplankton bloom also influence the location and extent of zooplankton concentrations. Relative and absolute differences in surface zooplankton between the "bgc-light" and the "standard" runs along the Bay mainstem are larger and found further downstream during higher flow summers (e.g., 2003-2004).
During the spring freshet and fall storms, differences in temperature and phytoplankton averaged over several days are much less than the hourly differences. This is in part due to diurnal variations in surface processes. During the daytime hours, stratification is enhanced by surface warming while vertical mixing episodically homogenizes the upper water column during the evenings or wind events. This vertical mixing is more frequent during the fall, whereas during the spring, sustained density stratification allows for the development of greater temperature gradients. In the case of the spring freshet example, differences in hourly surface phytoplankton concentrations at the mainstem stations span −33% to +60%, with a difference of only +2% averaged over those days and all the stations. Additionally, at station CB3.2 large relative increases in surface phytoplankton concentrations occurred after a period of high ISS concentrations. This suggests that as ISS settle out of the surface water, the warmer surface temperatures in the "bgc-light" run encourages growth, allowing phytoplankton to further utilize the abundant nutrients in the water.

Concluding Remarks
By comparing our two simulations, we showed how accounting for water clarity variability in the solar heating calculations can affect the estuarine system. In our "standard" run, we had a "one-way" coupled configuration, in which the physical parameters influenced biogeochemistry, but the light attenuation coefficient for solar heating, k h , was constant at all locations and times. In the "bgc-light" run, a "two-way" coupled configuration, the variations in k h resulted in greater surface temperature variability and decreased bottom temperatures.
The primary implication of this work for future water quality and climate change modeling studies is that running the models in a "one-way" coupled configuration underestimates temperature variability, both temporally and spatially. As we showed in our spring freshet example, the largest episodic differences in temperature may occur during a seasonal transition. Additionally, temperature shifts may be sustained over the summer months, as demonstrated by the colder summertime bottom temperatures in the "bgclight" run. Changes in vertical temperature gradients will have greater impacts on circulation and biogeochemistry in water bodies where temperature plays a large role in density stratification, as in Cahill et al. (2008) and Jolliff and Smith (2014), in contrast to the salinity-stratified Chesapeake Bay.
For models that assume a constant light attenuation coefficient for solar heating that is much different from the observed light attenuation in the water body, temperature differences and impacts on biogeochemistry will be larger in magnitude than reported in this study. It is difficult to assess the prevalence of this inconsistency because many modeling studies do not report the irradiance calculations used for solar heating. We recommend other modelers report their hydrodynamic light attenuation assumption, Jerlov type or otherwise, in their model documentation. For our model runs, the differences in temperature between the two simulations are much smaller than the difference between the model and observations, and neither systematically agrees better with the observations. Therefore, there was no change in the model agreement with in situ measurements of temperature from discrete bi-weekly to monthly measurements by the Chesapeake Bay Program along mainstem stations (see Appendix B).
In this study, we used a parameterization of light attenuation that was highly dependent on salinity, a proxy for CDOM. This additive model of light-attenuating materials allowed us to determine the contributions by different modeled biogeochemical constituents in the oligoand mesohaline regions of the Bay. While it likely adequately represents the spatial pattern of light absorption by CDOM (Rochelle-Newall and Fisher 2002), this parameterization predicted a light attenuation coefficient that was sometimes smaller than typical values of k d from in situ measurements in the polyhaline region. This is in part due to the lack of explicit resuspension of inorganic sediment in our implementation of ChesROMS-ECB. Here, resuspension is represented implicitly through the imposition of a lower bound on k bgc (0.6m −1 ). Future studies should investigate and incorporate the processes that affect water clarity in the lower reaches of the Bay, such as wave-induced sediment resuspension and coastal erosion.
Projected increases in precipitation and streamflow during the winter and spring due to climate change (Najjar et al. 2010;Irby et al. 2018) may bring more intense spring freshet events, depending on the suddenness of snow and ice melt. As shown in our study, decreased water clarity events may be followed by warmer surface temperatures during the spring-to-summer transition from cold to warm temperatures (Fig. 6a). These changes are in addition to the long-term temperature rise that has already been observed and is expected to continue. While warmer temperatures generally increase phytoplankton growth, it can also shift the distributions of specific phytoplankton taxa which could have implications for trophic interactions. As demonstrated in our study, increasing surface temperatures in the upper Bay can affect the distribution of nutrients by decreasing nutrient transport into the meso-and polyhaline regions of the Bay. This resulted in a decline in phytoplankton in these regions of the Bay, decreasing zooplankton populations. These indirect effects on higher trophic levels are poorly understood and warrant further study.
We demonstrated the impacts of variable water clarity on solar heating by showing the difference between the "bgc-light" minus "standard" simulations. The temperature changes associated with increasing water clarity would be the opposite of what we reported in this manuscript; i.e., more solar heating at deeper depths. While the recent recovery of submerged aquatic vegetation (SAV) in the Chesapeake Bay demonstrates the success of management efforts (Lefcheck et al. 2018), warmer subsurface and bottom temperatures in shallow areas may be an unintended consequence of water clarity improvements. Along with rising temperatures from climate change, these changes together may threaten recent gains in SAV restoration. While spectrally resolved light calculations have been incorporated into modeling studies for SAV habitat prediction (del Barrio et al. 2014), the impact of improved calculations for solar heating are unknown. On the other hand, we report how decreases in water clarity may be linked to colder bottom temperatures. Unusually cold winter temperatures can impact other organisms of importance to the Bay as well, such as mortality for blue crabs during severe winters (Bauer and Miller 2010). The interaction between water clarity and solar heating should be incorporated in future modeling studies involving species' habitability to better simulate environmental temperature variability. to values for the "standard" simulation. Green line and right hand side y-axis are the difference between the two simulations ("bgc-light" minus "standard")

Fig. 15
Daily average change in Bay-integrated a phytoplankton carbon, b nitrate, and c nitrate uptake by phytoplankton ("bgc-light" minus "standard"). Different colors correspond to different salinity ranges for integration Fig. 16 Difference in June to August water temperature at 5m depth along the mainstem ("bgc-light" minus "standard"). Mainstem locations shown in Fig. 2