Long-Term Trends in Chesapeake Bay Remote Sensing Long-Term Trends in Chesapeake Bay Remote Sensing Reflectance: Implications for Water Clarity Reflectance: Implications for Water Clarity

,


Introduction
Studying long-term change in water clarity in estuaries is an integral part of assessing improvements from historically polluted conditions.Light availability in estuaries is a critical driver of primary production and ecosystem health, shaping important nursery habitats such as seagrass meadows, coral reefs, and oyster reefs.Low water clarity is often concurrent with pathogens and harmful algal blooms, impacting fisheries and human health.Some of the world's estuaries have experienced widespread eutrophication and degraded water clarity, such as the Chesapeake Bay, the Baltic Sea, and the Wadden Sea (Cloern, 2001;Dupont & Aksnes, 2013;Kemp et al., 2005;van Beusekom, 2005).Other estuaries have recovered from past degradation and experienced improved water clarity conditions in recent decades, such as Tampa Bay, the San Francisco Bay, Danish estuaries and coastal waters, and the Black Sea (Boesch, 2019;Cloern & Jassby, 2012;Greening et al., 2014;Sherwood et al., 2016).In recovered areas, the positive effects of watershed management actions mandated to reduce sediment and nutrient inputs and increase light availability for aquatic plants are now being seen.Even so, despite nutrient reductions and related Abstract While ecosystem health is improving in many estuaries worldwide following nutrient reductions, inconsistent trends in water clarity often remain.The Chesapeake Bay, a eutrophic estuary with a highly populated watershed, is a crucial testbed for these concerns.Improved efforts are needed to understand why some measurements of downstream estuarine water clarity appear to be uncorrelated with watershed management actions, and multiple metrics of clarity are needed to address this issue.To complement in situ measurements, satellite remote sensing provides an additional tool with which to assess long-term change in water clarity.In this study, remote sensing reflectance (Rrs) from the Moderate Resolution Imaging Spectroradiometer on satellite Aqua was evaluated from 2003 to 2020 at multiple wavelengths for surface waters of the Chesapeake Bay.Trends show an overall long-term decrease in Rrs in the upper estuary for all wavelengths, yet an increase in Rrs in the lower estuary for green wavelengths.Trends in band ratios show longterm decreasing red-to-green and red-to-blue ratios, yet long-term increasing green-to-blue ratios.Seasonally, trends in band ratios were relatively consistent throughout the year and along-estuary, whereas single band reflectance trends varied seasonally and along-estuary.In the lower Bay, Septembers showed the strongest decreasing trends in red reflectance, while early spring and summer had the most pronounced increasing trends in green reflectance.These trends suggest that the system has experienced a long-term reduction in suspended solids concentration and light attenuation without a systematic reduction in chlorophyll-a concentration.

Plain Language Summary
In many estuaries around the world, waters have become less clear over time because of nutrient and sediment pollution.Watershed management efforts on land have reduced the delivery of nutrient and sediment from specific regions in recent decades.Through these efforts, some estuaries are showing improvements such as more seagrass and higher oxygen levels; however, the water does not always appear clearer.The Chesapeake Bay is one such place where not all measurements of water clarity show improvement, even though the ecosystem is recovering.Satellites provide snapshots of surface water color, covering large areas with frequent repeats in time.From 2003 to 2020, we analyzed long-term trends in surface water reflectance, or the brightness of the light coming from water at specific colors of visible light.We also analyzed the change over time in ratios between reflectance at different colors, or reflectance ratios.Reflectance ratios are often used to estimate common water quality variables, such as light attenuation, chlorophyll pigment concentration, and particle concentration.Our results show that there has been a decrease in red-to-green and red-to-blue reflectance ratios, suggesting improving water clarity over time in most of the Bay.However, increasing green-to-blue reflectance ratios suggest increasing chlorophyll concentration.TURNER ET AL.
The Chesapeake Bay serves as a prime case study for past estuarine eutrophication, recent improvements, and inconsistent water clarity response requiring further analysis.The watershed of the Chesapeake Bay includes rocky uplands and low-lying coastal plains with very different underlying geologies and drainage patterns that lead to varied river chemistry (Najjar et al., 2020) and likely to diverse types of dissolved and particulate matter entering the Bay from the different tributary rivers.In this estuary, water clarity is used in regional watershed management alongside chlorophyll-a (Chl-a) and oxygen to assess the health of the estuary (Tango & Batiuk, 2013), and watershed sediment and nutrient inputs are both actively managed for reduction (Shenk & Linker, 2013).The results of in situ clarity change measurements following cleanup have been inconsistent: the Chesapeake Bay has experienced multiple, sometimes divergent, historical trends in water clarity.Despite extensive water quality management efforts and recent documentation of reductions in riverine nutrient inputs (Lefcheck et al., 2018;Zhang, Hu, et al., 2018;Zhang, Tango, et al., 2018), Secchi disk depth declined in the Chesapeake Bay from the 1980s to the 2000s (Gallegos et al., 2011;Harding et al., 2016;Testa et al., 2019;Williams et al., 2010;Zhang et al., 2015).Current knowledge of water clarity from in situ observations fails to explain these incongruities.It is crucial that we understand the causes of discrepancies between management actions and water clarity results in the Chesapeake Bay.To this end, a more thorough understanding of the spatial and temporal patterns in water clarity metrics -including remote sensing reflectance -is needed.While in situ data are limited by spatial coverage, cruise frequency, and methods changes (Bever et al., 2013;Williams et al., 2010), remote sensors have been measuring surface reflectance for multiple decades and thus can substantially complement in situ programs.
The Chesapeake Bay has served as a testbed for remote sensing research for the last 40 years (Shelley, 1976;Son et al., 2014;Son & Wang, 2012, 2015;Stumpf, 1988;Stumpf & Pennock, 1989;Tzortziou et al., 2007;M. Wang et al., 2009;Werdell et al., 2009;Williamson & Grabau, 1973).Applications of water quality remote sensing have included: seasonal variability of dissolved organic carbon (DOC) retrievals (Cao et al., 2018;Mannino et al., 2008Mannino et al., , 2014)), DOC export to the coastal ocean (Signorini et al., 2019), the relationship between tidal energy and water clarity at the Bay mouth (Shi et al., 2013), timing of high Chl-a preceding low oxygen in the mainstem Bay (Zheng & DiGiacomo, 2020), and multi-sensor analysis of sediment plume trajectories in the upper Bay estuarine turbidity maximum (Zheng et al., 2015).Past studies have demonstrated effective use of the Moderate-resolution Imaging Spectroradiometer on NASA Earth Observation satellite Aqua (MODIS-Aqua) to estimate water quality variables in the mainstem Chesapeake Bay, including light attenuation (K d ), total suspended solids (TSS), Chl-a, absorption by colored dissolved organic matter (a CDOM ), and DOC, albeit with a range of uncertainties (Table S1 in Supporting Information S1).MODIS studies in the relatively bright, turbid Chesapeake Bay have included land bands, wavelengths at which image collection is optimized for bright land surfaces, in addition to traditional ocean bands, which are optimized for dark water surfaces (Feng et al., 2018;Franz et al., 2006).
In general, remote sensing reflectance (Rrs) patterns across the visible wavelengths (400-700 nm), i.e., Rrs spectra, illustrate the apparent color and brightness of water bodies depending on concentrations of CDOM, phytoplankton, and other particles, such as organic detritus and mineral sediments (Babin et al., 2003;Bidigare et al., 1990;Hawes, 1992;Kirk, 1994;Morel & Prieur, 1977;Vaillancourt et al., 2004).In open ocean waters, peak Rrs values occur in the blue portion of the visible spectrum.However, in estuaries, light absorption by CDOM and chlorophyll combined with absorption and scattering by particles shift the peak Rrs values into the green and red portions of the visible spectrum compared to the open ocean (Roesler & Perry, 1995).Surface water Rrs spectra relate to whether waters appear green, yellow, or brown in color.For example, in the somewhat darker central channel of the Chesapeake Bay, surface waters experience peak Rrs values in the green wavelengths with a relatively low maximum Rrs (Tzortziou et al., 2007), while in the more turbid, bright sediment-rich tributary rivers, surface waters experience peak Rrs values in the yellow-to-red wavelengths with a higher maximum Rrs.
Of the past approaches for estimating Chesapeake Bay water clarity from remote sensing, the most widely used method involves the MODIS high-resolution land band at 250 m resolution for the 645 nm wavelength (Aurin et al., 2013;Crooke et al., 2017;DeLuca et al., 2018;Hasan & Benninger, 2017;Ondrusek et al., 2012;Yunus et al., 2021).In the present study, we consider more wavelengths and band ratios in addition to this widely used MODIS red-band approach.While most past studies focused on derived variables, uncertainties were high, as the use of derived variables may mask important particulate, planktonic, and dissolved contributions to long-term 10.1029/2021JC017959 3 of 20 trends (Zheng & DiGiacomo, 2017).Rrs values themselves have been used in other estuaries to study general patterns (Tao & Hill, 2019), and ratios are the underlying metric of a major amount of ocean color algorithms (Dierssen, 2010).Thus, the Rrs approach is used in the present study for Chesapeake Bay to further discern change over time.
To date, there have been few studies of long-term change in Chesapeake Bay water clarity using satellite remote sensing reflectance data.Satellite estimates have been established for many water quality variables using multiple algorithms for the Bay, and retrieval of those variables has been well-documented; however, these data are not often used to answer science questions about temporal trends in the mainstem Bay.Change over time has been heavily studied using in situ cruise observations, yet analysis of long-term change in satellite-derived clarity indices has been limited to red-to-green Chl-a (Le et al., 2013), briefly noted for Rrs(645)-based clarity estimates (Crooke et al., 2017), and analyzed for Rrs(555) in a specific mid-estuary region of the Bay (Gallegos et al., 2011).Past studies of long-term trends in satellite-derived water quality, though useful, only extended to the early 2010s, missing the recent years characterized by continued improvements in watershed management.
The objective of this work is to examine water clarity change in the Chesapeake Bay mainstem over the past two decades using Rrs and band ratios from MODIS-Aqua based observations.Here the term water clarity is used with regard to two types of water constituents, with clearer water corresponding to lower amounts of the constituents: (a) suspended particulate matter, closely related to overall turbidity (often related to red-band reflectance, red-to-green ratios, and red-to-blue ratios) and (b) phytoplankton, closely associated with chlorophyll concentration (often related to green-to-blue reflectance ratios).Since the goal of the proposed work is primarily to highlight the general magnitude and direction of any impactful long-term temporal trends in Rrs, the absolute accuracy and precision of water clarity variable estimates potentially calibrated from Rrs are less important here than overall Rrs patterns.That is, trends and patterns in the remote sensing data are useful in answering our research questions despite some potential biases in the values themselves.We aim to answer the questions: How have Chesapeake Bay Rrs(λ) and band ratios evolved over the last two decades?What do common and/or distinct trends found among groups of bands and ratios suggest in the context of long-term change with regards to water clarity?Answering these questions will allow us to better understand how major contributions to Chesapeake Bay water clarity, such as band ratios commonly associated with TSS, K d , and Chl-a, have changed over the last two decades.

In Situ Rrs Data
Field observations in the mainstem Chesapeake Bay were obtained from the SeaWiFS Bio-optical Archive and Storage System (SeaBASS) (Werdell & Bailey, 2002;Werdell et al., 2003).In situ, Rrs at the water surface was quantified from concurrent values of upwelling radiance and downwelling irradiance, using two connected radiometers pointing downward into the water column and upward at the sky, respectively.Thus, Rrs is the surface ratio of upwelling radiance emerging from water to downwelling radiative flux in air (Mobley et al., 2005;O'Reilly et al., 1998).Because upwelling radiance is measured at a specific angle and downwelling irradiance is measured on a flat plane receiving light from a half-hemispherical area above the water surface, Rrs has units of per steradian (sr −1 ).In situ Rrs measurements were reported from a wide range of times and locations throughout the mainstem Bay, including multiple seasons and multiple salinity regimes, spanning the years 2005-2014 (Figure 1, Table 1).

Satellite Rrs Data
MODIS-Aqua satellite data from January 2003 to December 2020 were used to study long-term trends in Rrs (Table 2).Aqua is a polar-orbiting NASA Earth Observation satellite, local ascending at the equator, with a 1:30pm local standard time overpass.MODIS ocean bands (bands 8-14) have high gain settings in order to sense the characteristics of dark water bodies.In contrast, the land bands used here (bands 1, 3, 4) have lower gain settings to detect brighter features on land, but can be used in brighter inland and coastal water bodies (Feng et al., 2018;Franz et al., 2006).Rrs at the water surface must be calculated from the sensor measurement of total radiance exiting the top of Earth's atmosphere through a process of atmospheric correction to remove the contributions of aerosols (Ahmad et al., 2010;Gordon & Wang, 1994).After additional corrections to remove sun glint, whitecaps, and other artifacts, less than 10% of top-of-atmosphere reflectance is contributed by ocean water (Kirk, 1994;Martin, 2014).
Data were processed from instrument data (Level-1; see Appendix A) through atmospherically corrected and spatially binned data (Level-3) using a custom merging method following Aurin et al. (2013).Traditional atmospheric correction in highly turbid waters may cause data loss.In this paper the terms "turbid" and "turbidity" are used as qualitative descriptors of sediment-laden water, not strictly with regard to the quantitative measure by a turbidity sensor.In estuaries, lakes, and rivers, atmospheric correction often aliases bright sediment-laden waters as atmospheric haze or clouds due to high emission in the infrared.In these environments, standard near-infrared atmospheric correction removes the more turbid pixels, biasing aquatic retrievals toward low-turbidity conditions and underrepresents the available data.In the Chesapeake Bay, standard processing would thus "miss" many images of conditions following storm events where high winds and/or strong precipitation caused especially turbid, bright conditions in surface waters (Figure S1 in Supporting Information S1).Therefore, in the present study a custom atmospheric correction merging method was performed, using both shortwave infrared and near-infrared, to ensure that the high-turbidity data were included in the analysis.
Following the methods of Aurin et al. (2013), we used two atmospheric correction methods, shortwave infrared for high-turbidity and near-infrared for low-turbidity pixels, to process each scene into Rrs values (analysis-ready but not spatially binned; i.e., Level-2), using each separate method and then merging those scenes.The present study diverged from Aurin et al. (2013) in that the median negative offset was not added back to Level-2 low-turbidity scenes, as it caused misalignment of Rrs values for the same band at different spatial resolutions (e.g., 555 nm).Spatially, high-turbidity and low-turbidity Rrs were processed to Level-2 at the nominal spatial resolution of each band.Rrs(645) was processed to all spatial resolutions (250 m, 500 m, and 1 km) in order to facilitate the merging method.At all spatial resolutions, scenes were merged along spatial guidelines set by a Rrs(645) threshold value of 0.01 sr −1 .Where Rrs(645) < 0.01 sr −1 , values from the low-turbidity processing method were used.Where Rrs(645) > 0.01 sr −1 , values from the high-turbidity processing method replaced values from the low-turbidity method to create the merged (Level-2B) scene.Using these merging methods, up to twice as many scenes per month could be included in monthly averages (Figure S1 in Supporting Information S1).Custom-processed Level-2B Rrs(λ) scenes were binned to Level-3 monthly composites at the respective spatial resolutions for each band to facilitate trend analysis over a consistent spatial grid.
Data quantity was somewhat limited by clouds and proximity to land.Although many scenes were excluded from the data set due to clouds and other artifacts, no marked seasonal bias in cloud cover was found.Approximately 4-8 scenes were included in each monthly composite (Figure S2 in Supporting Information S1).The number of points in monthly composites showed only very small variation between spatial resolutions, quantified by comparing Rrs(645) level-3 mapped images for 250 m, 500 m, and 1 km spatial resolutions.The slight decrease in points per month with coarsening spatial resolution was most relevant to coast-adjacent areas (Figure S3 in Supporting Information S1), for which long-term trends were not analyzed in the current study.
Spatially, analysis focused on the mainstem Bay and lower reaches of the large tributary rivers (Figure 1).Satellite data points <750 m from shore and in smaller tributaries were excluded; some points farther from shore were also excluded due to additional data quality control, such as partial to full   1).Analysis focuses mainly on upper versus lower Bay.When the term mid-Bay is used to help describe spatial patterns in results, the term mid-Bay refers to the mainstem part of the estuary between the southern end of the Chester River mouth and the southern end of the Potomac River mouth (Figure 1).Additionally, analysis quantifies the spatial percent of the Bay (percent of water surface area) experiencing significant trends over time.In climatologies, upper Bay and lower Bay include only mainstem waters and thus exclude tributary rivers.
Validation of MODIS-Aqua Rrs for the Chesapeake Bay region was performed using SeaBASS data spanning years from 2005 to 2013 (Table 1; Table 3).In situ Rrs validation data points were all collected within 6 hr of a MODIS-Aqua overpass (Figure S4 in Supporting Information S1).Satellite Rrs at station locations were extracted via spatial matchup windows depending on the spatial resolution of the relevant wavelengths, using an aerial coverage of ∼1.6, 2.25, and 1 km 2 for the three respective spatial resolutions 250 m, 500 m, and 1 km.Rrs values were compared by individual bands at each wavelength's nominal spatial resolution (Table 3, Figure S5 in Supporting Information S1) and at 1 km (Figure 2).Metrics for satellite skill assessment included the mean ratio, bias, mean absolute error, root mean squared error, mean absolute percent difference, and correlation coefficient (Text S1 in Supporting Information S1).Comparison with in situ Rrs revealed generally close matches between satellite Rrs and in situ data.Overall, satellite Rrs slightly underestimated in situ Rrs, with overestimation in the blue wavelengths (<488 nm) and underestimation in the green through red wavelengths (>488 nm) (Table 3, Figure 2).Skill was high for specific single bands.For example, the closest match to in situ Rrs was measured at 488 and 645 nm (Figure 2).There was a small effect of spatial resolution on skill of satellite Rrs retrieval.Satellite Rrs(645) was consistently slightly lower than in situ Rrs(645).Skill of satellite Rrs was minimally different among different spatial resolutions (Table S4 in Supporting Information S1, Figures S6, and S7 in Supporting Information S1).

Calculation of Long-Term Trends
Monthly composite images were used for trend analysis to maximize spatial coverage and to maintain consistent sampling intervals for each year.Trends in Rrs at each wavelength will be referenced heretofore as single band trends.Analysis of trends in single bands used the nominal spatial resolution of each band.Analysis of trends in band ratios were calculated at 1 km spatial resolution.For each band ratio, for each combination Rrs(λ 1 )/Rrs(λ 2 ), where λ 1 > λ 2 (e.g., Rrs(555)/Rrs(469) rather than vice versa).Prior to long-term trend analysis, monthly climatologies (Figure S8 in Supporting Information S1) were analyzed to explore seasonal variability in each single band and each band ratio (see Section 3.1).During the monthly climatology stage   of analysis, Rrs(667) and Rrs(678) were found to have high noise over space and time and were thus excluded from long-term trend analysis due to their patchiness and extreme variability along with associated band ratios.
Additionally, Rrs(412) was removed from long-term trend analysis due to its large deviation from observed Rrs(412) compared to other bands (37% mean APD compared to other bands at 1%-15%; Table 3), likely caused by known atmospheric correction issues in coastal waters.To determine which parts of the seasonal cycle in Rrs were most closely associated with the overall long-term trends, trends for each month from 2003 to 2020 were assessed.Seasonal adjustment had very little effect on the results of long-term trend analysis (Figures S9 and S10 in Supporting Information S1); therefore, seasonal cycles are presented here but not subtracted during trend analysis.Rather, seasons are described in terms of which months are the most aligned with the overall long-term trends for single bands and band ratios (see Section 3.4).
For each single band and each band ratio, spatially explicit trends were calculated for spatially binned and mapped pixel locations (Level-3 data) that contained data for >80% of the months in the time series from 2003 to 2020, i.e., >173 of 216 monthly composites.Trends were calculated as linear regressions using the slope of the least squares fit (i.e., Greene et al., 2019).Significant trends were assigned using 90% confidence (p < 0.1).While conventional definitions of significance frequently use 95% confidence, a more moderate confidence level is appropriate here considering the strong seasonal and interannual variability.This definition of significance is used in the present study to describe the general strength of trends over time, recognizing that some trends may still be meaningful in the broader context even if they do not meet this definition.Trends in all band ratios were examined.For quality control, unrealistic ratios (i.e., Rrs(λ 1 )/Rrs(λ 2 ) less than zero or greater than ten) were removed prior to trend analysis.Certain band ratios were examined more closely if the trends were both significant over a substantial area of the Bay (greater than 15% of the area analyzed showed trends where p < 0.1) and large in magnitude (median trend within that area greater than or equal to +/-0.003 yr −1 ).Additionally, trends were calculated in terms of relative change (% yr −1 ) using the trend per year normalized to the long-term mean value of each pixel for each single band and each band ratio: Relative change (yr −1 ) = trendRrs()(sr −1 yr −1 ) meanRrs()(sr −1 ) (1) Relative change (yr −1 ) = trend Rrs(1) Rrs(2) MODIS is well suited for time series analysis despite the potential sensor drift and/or increased noise associated with its extended time in orbit.MODIS has been used for time series analysis in estuaries, lakes, and coastal waters such as, to name just a few: the Columbia River Estuary (Hudson et al., 2017), the coastal waters of the northern Gulf of Mexico (Reisinger et al., 2017), the Patos Lagoon in Brazil (Tavora et al., 2019), fjords in Svalbard (Payne & Roesler, 2019), coastal waters in the Arctic (Chaves et al., 2015), the Bohai Sea (Shang et al., 2016), and Minnesota lakes (Knight & Voth, 2012).To correct for potential drift, MODIS-Aqua is calibrated onboard regularly using scheduled solar diffuser measurements and lunar observations (Wu et al., 2013).MODIS-Aqua is collecting data in orbit past its expected performance lifetime, yet studies show that the sensor still has sufficient signal to noise ratio and the solar bands are still well-calibrated for all visible wavelengths (Lee et al., 2019;Xiong et al., 2020).

Seasonal Cycles in Rrs Single Bands and Band Ratios
Single band Rrs showed spatially heterogeneous seasonal cycles that varied in magnitude and seasonal pattern by wavelength (Figure 3).Green band Rrs (531-555 nm) showed a smaller magnitude seasonal cycle in the upper Bay in terms of percent change, with maximum in January through April, yet a relatively dynamic seasonal cycle in lower Bay, with larger values in January and October and lower values through the summer months May to August (Figures 3d-3f ).Rrs (645) had the most distinct one-per-year seasonal cycle in both the lower and upper Bay, with highest values in winter and lowest in summer, with upper Bay April and September months showing large standard deviation among years (Figure 3g).In general, the highest standard deviation over the time series was seen for the month of September, especially in the green and red bands (531-645 nm) (Figures 3d-3g).
Band ratios showed seasonal cycles with varying magnitudes and patterns (Figure 4).Red-to-green ratios and red-to-blue ratios were slightly higher in winter months (November-March) than summer months (May-September) (Figures 4a-4d).Green-to-blue ratios were highest in both April and fall months (September to December) (Figures 4e-4h).For ratios using the red band Rrs(645), standard deviation was markedly wider in the month of September (Figures 4a-4d).Spatially, seasonal cycles were not markedly different in shape, i.e., general seasonal pattern, between upper and lower Bay for any specific ratio.However, all band ratios were higher in value in the upper Bay than in the lower Bay year-round.

Change Over Time in Single Band Rrs
Generally, Rrs at most wavelengths decreased over time in the upper Bay and increased over time in the lower Bay, although the spatial extent, magnitude, and significance of those trends varied among wavelengths (Table 4, Figure 5).Blue wavelength Rrs(443 and 469) decreased over time in the upper Bay, with strong decreasing trends covering 11%, and 20% of the surface area of waters analyzed (Table 4, Figures 5a and 5b).Rrs(488) upper Bay decreases were not as spatially extensive as other blue bands.In terms of increasing trends, increases over time in Rrs(443 and 469) were found near the Bay mouth (Figures 5a and 5b).Increases in Rrs(488) were found for most of the lower Bay (26% of water surface area), extending spatially up-estuary to Tangier Island and reaching eastto-west across most of the mainstem Bay (Figure 5c).In the green wavelengths, upper Bay decreases and lower Bay increases in Rrs(531, 547, and 555) (Figures 5d-5f) closely resembled spatial patterns in Rrs(488), except all three green wavelengths additionally showed a region of significant long-term increase in Rrs in the lower James River.Green Rrs values showed spatially extensive increasing trends in the lower Bay, with 24%-27% of the area analyzed showing a strong increase over time (p < 0.1, Figures 5d-5f, Table 4).In the red, decreases in Rrs(645) were found in 14% of the area of the Bay (Table 4), including the upper Bay extending down-estuary to just above the Choptank River mouth, near the mouth of the Patuxent River, and in the lower James River (Figure 5g).

Change Over Time in Band Ratios
For band ratios, red-to-green ratios and red-to-blue ratios generally showed consistent, widespread (greater than or equal to 15% of area analyzed) decreases, while green-to-blue ratios generally showed long-term increases (Figure 6).Strong temporal trends (median trend greater than or equal to +/-0.003 yr −1 ) were found throughout the mainstem Bay for three red-to-green band ratios, one red-to-blue ratio, and four green-to-blue ratios (Table 5).
The largest magnitude and most spatially widespread decreases of the red-to-green ratios was Rrs(645)/Rrs(531),  Median trend 2003 to 2020 of area exhibiting significant long-term increase or decrease (p < 0.1).b Median relative trend of area exhibiting significant long-term increase or decrease, relative to the long-term mean.c Percent of water area analyzed exhibiting trend (p < 0.1).

Seasonality of Trends Over Time
On a month-to-month basis, Rrs trends from 2003 to 2020 generally showed decreases over time in the upper Bay and increases over time in the lower Bay.However, in both regions of the Bay, seasonality of trends varied by wavelength (Figure 7).Upper Bay decreases in blue band reflectances (Rrs(443) through Rrs( 488)) were most pronounced (p < 0.1 for greater than 25% of the region) in June (Figures 7a-7c) and April (Figure 7b).Upper Bay decreases in green band reflectances Rrs(547) and Rrs(555) were strongest in May (Figures 7e and 7f).Upper Bay decreases in Rrs(645) were seasonally more variable than for blue and green Rrs, showing strong decreases in May and December, yet weak increases (p < 0.1 for less than 25% of the region) in February and March (Figure 7g).Lower Bay increases in blue bands were strongest in January (Figures 7a-7c), June (Figures 7a and 7b), and for Rrs(488) (more similar to green band patterns) in March, June, and July (Figure 7c).Like Rrs(488), green bands Rrs(531) and Rrs(547) showed strong increasing trends in the lower Bay in January, March, June, and July (Figures 7d and 7e).For Rrs(645) in the lower Bay, the months March and June had strong increasing trends while the month of September had a strong decreasing trend (Figure 7g).Because of variability in sign for lower Bay Rrs(645) over the year, the general long-term trend for most of the lower Bay is negligible (Figure 5g).For all  469) and (f) Rrs( 555) at 500 m, and all other bands at 1 km.Small black dots highlight significant trends (p < 0.1).Trends are expressed as relative trends normalized to the long-term mean Rrs at each location (% yr −1 ).Absolute values of median increasing and decreasing trends over time were 0.8%-1.0%yr −1 and 1.3%-2.9%yr −1 respectively (see Table 4).
single band Rrs in the lower Bay, trends over time in April and May months from 2003 to 2020 were negligible (Figure 7).
For band ratios on a month-to-month basis, long-term trends were more consistent in sign both spatially (upper Bay vs. lower Bay) and seasonally than single band trends (Figure 8).For the upper Bay, red-to-green and red-to-blue ratios show strong decreases in May and December (Figures 8a-8d) with Rrs(645)/Rrs(555) showing an additional significant decrease in April (Figure 8a).Upper Bay green-to-blue ratios showed strong increases particularly in February, June, August, and September (Figures 8e-8g), as well as increases for additional months in the case of Rrs(531)/Rrs(469) and Rrs(547)/ Rrs(469) (Figures 8f and 8g).For lower Bay red-to-green ratios, the months of April, May, September, and December showed the strongest (p < 0.1 for greater than 25% of the region) decreasing long-term trends (Figures 8a-8c).
Red-to-blue ratio Rrs(645)/Rrs(488) lower Bay trends were strongest in January, May, September, and December (Figure 8d).Lower Bay green-to-blue ratios showed the strongest increases in August, September (Figures 8e-8g), and November (Figure 8f).The ratio Rrs(488)/Rrs(469) showed the most seasonally consistent strong trends of all the ratios, increasing in the both the ), including (a-c) decreasing red-to-green ratios, (d) a decreasing red-to-blue ratios, and (e-h) increasing green-to-blue ratios.Selected trends shown were > +/− 0.003 yr −1 in magnitude and significant (p < 0.1) for >15% of the Bay.Trends are expressed as relative trends normalized to the long-term mean value at each location (% yr −1 ).Absolute values of median trends over time ranged from 0.5% to 0.9% yr −1 (see Table 5).Trends > +/− 0.003 yr −1 in magnitude that were significant (p < 0.1) for >15% of the Bay.b Median trend of area analyzed exhibiting significant longterm trend (p < 0.1).c Percent of water area analyzed exhibiting significant long-term trend (p < 0.1).

Table 5
Strong Trends a in Band Ratios Rrs(λ 1 )/Rrs(λ 2 ) upper and lower Bay throughout the year (Figure 8h).For green-to-blue ratios, though ratios were nearly always increasing over time, the magnitudes of the increasing trends were slightly more variable month-to-month in the upper Bay and more seasonally consistent in the lower Bay (Figures 8e-8h).Upper Bay spatial variability in the absolute magnitude of trends was generally higher than lower Bay spatial variability for both single bands (Figure 7) and band ratios (Figure 8).Seasonally grouped trends, i.e., winter (December-February), show similar results to month-to-month trends (Figures S11 and S12 in Supporting Information S1).

Relevance of Single Band Rrs Trends
The red-band Rrs approach has historically been the most widely used satellite-derived index of water clarity for estuaries having similar turbidity conditions as the Chesapeake Bay (Crooke et al., 2017;Constantin et al., 2016;DeLuca et al., 2018;Hudson et al., 2017;Ondrusek et al., 2012;Stumpf & Pennock, 1989;Yunus et al., 2021).In moderately turbid estuaries like the Chesapeake Bay, red wavelengths (∼640-660 nm) are most commonly used to estimate suspended sediments or TSS concentrations due to the relatively lower influence of phytoplankton and CDOM in this red portion of the visible spectrum (Tzortziou et al., 2007).Our findings using this red-band approach (Figure 5g) suggest that clarity is improving more substantially in the upper Bay than the lower Bay.
In part, the larger Rrs decrease in the upper Bay may be due to higher initial signal in the early 2000s which contained relatively more hydrologically wet years (Figures S13 and S14 in Supporting Information S1, see Section 4.2).Furthermore, larger Rrs decreases in the upper Bay may be related to the proximity of this region to sediment inputs; the upper Chesapeake Bay is sometimes considered an estuarine turbidity maximum for the entire estuary due its response to inputs from the largest river, the Susquehanna (Najjar et al., 2020).Thus, close to the source of sediments, riverine sediment reductions may be more closely related to decreasing red-band Rrs over the last two decades.
Increased green Rrs is often associated with increased green particulate backscattering.Spatially, the increasing trends observed in the present study are seen mostly in the lower Bay, down-estuary from the analysis by Gallegos et al. (2011).Possibly, increased Rrs(green) may be associated with different shaped particles or different morphologies of phytoplankton cells (i.e., Organelli et al., 2018), or with increasing colloidal organic matter or bacteria (Balch et al., 2000;Morel, 1991;Stramski & Kiefer, 1991;Stramski & Woźniak, 2005).The strong increasing trends in Rrs(green) for 2003 to 2020 in the lower Bay over large spatial coverage, especially in summer (Figure 7), suggest that an increase in green backscattering and an increase in organic matter associated with phytoplankton (e.g., Gallegos et al., 2011) may be occurring in summer more strongly than the rest of the year.
Remote sensing thus provides further evidence that Chesapeake Bay summers are getting greener, aligning with the conclusions of Testa et al. (2019) that trends in long-term shallowing Secchi depth are strongest in summer months.

Relevance of Band Ratio Trends
Red-to-green and red-to-blue ratios have been used in past estuarine studies to estimate water clarity metrics such as TSS, K d , and turbidity (Table S3 in Supporting Information S1).Red-to-green ratio algorithms are often helpful in optically complex coastal waters because they minimize the interference by non-algal particle absorption and CDOM absorption in the blue bands (Abbas et al., 2019;Ioannou et al., 2014;Le et al., 2013;Tzortziou et al., 2007).For example, the red-to-green ratio Rrs(645)/Rrs(555) has been employed to estimate TSS (Reisinger et al., 2017) and turbidity (J.Wang et al., 2021) in the coastal Gulf of Mexico and Pearl River Estuary, respectively.The corresponding red-to-green ratio in our study showed spatially widespread decreases over time for the mainstem Chesapeake Bay (Figure 6a), suggesting a potential improvement in TSS or turbidity over time.
Red-to-blue ratios have previously been used to estimate K d and TSS in coastal waters (Tables S1 and S3 in Supporting Information S1; Siswanto et al., 2011;Son & Wang, 2012;Tomlinson et al., 2019;M. Wang et al., 2009).
In the present study, results showed only one consistently decreasing red-to-blue ratio (Figure 6d).In short, our analysis of red-to-blue ratios suggests that water clarity may have improved since the early 2000s, especially in the lower Bay, yet conclusions are more strongly supported by red-to-green ratio findings.
The increases in green-to-blue ratios in the Chesapeake Bay since the early 2000s (Figures 6e-6h) suggest that more blue light is being absorbed relative to green, an indicator of increased Chl-a or increased a CDOM .Since both phytoplankton and CDOM preferentially absorb more strongly in blue wavelengths than in green, here focus is placed on a potential increase in blue light absorption over time rather than a potential increase in green light absorption.Green-to-blue band ratios algorithms have been used to retrieve both a CDOM and Chl-a (Table S3 in Supporting Information S1).Historically, green-to-blue ratios have been widely used to estimate satellite-derived Chl-a even in the optically complex waters of the Chesapeake Bay (e.g., Werdell et al., 2009).For example, Ogashawara et al. (2014) describe the OC2 and OC3 algorithms used to estimate Chl-a from high-resolution (500 m) MODIS imagery with the band ratio Rrs(555)/Rrs(469).The results of the present study point to a longterm increase in Chesapeake Bay for the band ratio relevant to this Chl-a algorithm (Figure 6e).These results suggest a long-term increase in Chl-a as measured by green colored pigment concentration throughout most of the central and lower Bay.

Implications of Rrs Long-Term Trends
Long term trends in reflectances and reflectance ratios suggest different implications for water clarity in terms of the different constituents that can be estimated from remote sensing data.While red-band reflectance, red-togreen ratios, and red-to-blue ratios suggest improving water clarity over time in terms of suspended particulate matter and overall clarity, green-to-blue ratios suggest an increase in chlorophyll concentration or organic matter.
The results of the present study show that Rrs(645), the most commonly used for water clarity variables TSS and turbidity (Table S3 in Supporting Information S1), has decreased from 2003 to 2020 in the upper Bay but not in the lower mainstem Bay (Figure 5g).These results suggest that decreasing watershed sediment inputs have the greatest impact in the turbid sections of the Bay closest to river inputs.In contrast, our decreasing red-to-green and red-to-blue ratios (Figure 6) result from trends in green and blue dominating changes in red in the lower Bay (Figure 5).Decreased sediment inputs can increase light availability, in turn increasing down-estuary organic matter production from phytoplankton (Turner et al., 2021).The results of our work support this interpretation, since Rrs in the green bands (Figures 5d-5f) and green-to-blue ratios (Figures 6e-6h) were found to be increasing over time in the lower Bay.Decreasing red-to-green and red-to-blue ratios suggest that sediment inputs may be a weaker contributor to water clarity patterns in the lower Bay in recent years compared to the early 2000s.Increasing Rrs(green) and green-to-blue ratios suggest that phytoplankton and associated organic matter may have had greater impacts on water clarity patterns in recent years in the lower Bay (CBP, 2021;Testa et al., 2019).
The relative influence of rivers varied over the past two decades, as there were several particularly wet years in the early 2000s (Figure S13 in Supporting Information S1, Harding et al., 2019).September's widespread decreases in Rrs(645) (Figure 7g), red-to-green ratios, and red-to-blue ratios (Figures 8a-8d) suggest that earlier years (2003,2004) had more turbid September months than later years in the time series (Figure S14 in Supporting Information S1).The month of September also has the largest variability in streamflow over the time series (Figure S13 in Supporting Information S1).This seasonally specific long-term decrease in single-band and band ratios associated with turbidity and TSS (Table S3 in Supporting Information S1) may thus be associated with August and September hurricanes during 2003-2005.
Single band increases in some blue wavelengths and all green wavelengths in the lower Bay may also be related to the reduced influence of CDOM from rivers absorbing light in the 2010s.CDOM is strongly associated with river discharge and proximity to river inputs, yielding a strong negative correlation between CDOM and salinity, and CDOM generally absorbs light strongly in the blue wavelengths and somewhat into the green wavelengths.Results show that all single bands from Rrs(443) to Rrs(555) consistently increased in notable regions of the lower Bay (Figure 5).The fact that the trends span the Bay mouth and up the eastern side of the lower Bay suggest that this pattern may be related to water color trends associated with ocean-sourced water, since the path for higher salinity ocean-sourced water from the ocean moves up along the eastern side of the lower Bay due to Coriolis (Norcross et al., 1962;Pritchard, 1952).Relatively lower river CDOM influence due to more dry years in the 2010s may have allowed ocean-sourced waters bright in the blue-green part of the spectrum to contribute more strongly to lower Bay reflectance patterns in recent years compared to the early 2000s.
At the head of the estuary, more sediment is now bypassing the Conowingo Dam than in the early 2000s due to reservoir infilling (Cerco & Noel, 2016;Palinkas et al., 2019;Zhang et al., 2016).However, since satellite data analysis in this study does not account for flow normalization, the wetter years in the early 2000s versus the drier years in the mid-2010s may contribute to observed decreases in red-to-blue and red-to-green ratios over the period of 2003-2020 despite the infilling of the Conowingo Dam.Generally, drier years with lower relative river influence in the latter half of the time series may be a strong driver of the long-term water clarity change in the lower Bay.

Directions for Future Work
For future remote sensing studies of Chesapeake Bay water clarity, we suggest increasing spatial coverage of in situ data, exploring time series from other satellites, looking to higher spatial and spectral resolutions in the future, and incorporating neural network approaches in addition to empirical and analytical algorithms.For in situ data (Figure 1), more spatial coverage is needed for Rrs and other measurements relevant to satellite studies of water quality.Currently, there are gaps in Rrs data spatially, including estuarine turbidity maxima of major tributaries, nearshore areas, and marsh-adjacent waters.Increasing the variety of environments in the Bay where in situ ocean optics data are collected would improve future remote sensing analysis.For time series, Landsat heritage missions may be useful for coastal waters with high enough signal-to-noise ratio (Barnes et al., 2014;Pahlevan et al., 2018), but have not yet been applied in Chesapeake Bay.For certain applications, high-spatial-resolution sensors like Landsat 8, Sentinel-2 Multi-spectral Imager and Sentinel-3 OLCI can be used with great success (e.g., Cao et al., 2018).With newer sensors using more bands in the red edge portion of the spectrum, satellite-derived chlorophyll-a fluorescence can be used to detect blooms (Gilerson et al., 2015), providing a path forward for using ocean color satellite data operationally in coastal waters (Wilson, 2011).Future work should also consider modeling more complex processes, such as particle settling velocities, using satellite remote sensing (Nasiha et al., 2019).Recently, neural network approaches have been incorporated into methods aimed at merging future sensors with past and current sensors.For example, neural network approaches have been used in modeling hyperspectral water quality retrievals (Ibrahim et al., 2016), estimating TSS using a thorough multi-algorithm switching process (Balasubramanian et al., 2020;Wei et al., 2021), and employing similar multi-algorithm switching and neural networks to derive multiple other in situ variables (Fan et al., 2021).When future applications require very fine spatial resolution, unoccupied aircraft systems are advantageous for mapping water clarity or bloom events (Windle & Silsbe, 2021).
Upcoming satellite missions have the potential to greatly improve our understanding of water clarity in optically complex estuarine waters.Future missions, including Plankton, Atmosphere, Clouds, and ocean Ecosystems (PACE), Surface Biology and Geology (SBG), and the Geosynchronous Littoral Imaging and Monitoring Radiometer (GLIMR) will bring hyperspectral and polarimetric measurements into the suite of earth observation satellites (Cawse-Nicholson et al., 2021;Cetinić et al., 2018;Petro et al., 2020;Salisbury & Mannino, 2020).Importantly, societal applications of future missions weigh heavily in the planning process, including monitoring and predicting the impacts of climate change on marine resources and food security (Uz et al., 2019;Wolny et al., 2020).These future missions will allow for improved characterization of phytoplankton community composition and better estimates of particle size distributions.In particular, polarization may allow for improved partitioning of organic matter versus inorganic matter using the real part of the refractive index and the degree of linear polarization of water leaving radiances (Gao et al., 2021;Loisel et al., 2008;Stamnes et al., 2018;Tonizzo et al., 2009).These planned developments have the potential to further our understanding of water clarity in the Chesapeake Bay using satellite remote sensing.

Conclusions
This research represents a holistic analysis of MODIS-Aqua Rrs and band ratios to quantify change over time in the Chesapeake Bay estuarine complex.Results suggest water clarity improvement, potentially associated with reduced TSS and K d , as evidenced by decreasing red-to-green and red-to-blue ratios throughout the mainstem Bay.However, an increase in green band Rrs in the lower Bay suggests a possible local increase in Chl-a.In contrast, a decrease in Rrs across all bands in the upper Bay implies a decrease near the estuarine turbidity maximum in Chl-a as well as TSS.Seasonally, trends showed summer decreases in upper Bay blue-band Rrs, early spring and summer increases in lower Bay green-band Rrs, and variable trends month-to-month in red-band Rrs, albeit a strong September decrease in the lower Bay.Throughout the year, red-to-green and red-to-blue ratios generally decreased 2003 to 2020 while green-to-blue ratios generally increased.Seasonally and spatially along-estuary, more variability was seen in trends for single band Rrs values than for trends in reflectance ratios.
Our work highlights that the optical complexity and time scale of a given research question is critical for careful algorithm selection in optically complex estuarine waters.In contrast to studying long-term trends in the Chesapeake Bay, short-term applications of remote sensing data with stable calibrations to specific water column constituents or features can provide powerful insights.However, analysis of Rrs or Rrs ratios as variables instead of other derived variables (Chl-a, TSS, K d ) should be considered when targeting temporally evolving, optically complex regions like the Chesapeake Bay.The results of this study further support the idea that at long time scales, remote sensing of water clarity can complement, but not replace, in situ monitoring and help improve and direct sustainable, long term coastal system economic and environmental management efforts.Land band Wavelength or spectral range at which a satellite sensor has been optimized to image bright land surfaces, for example, the 645 nm band on MODIS sensors.
Level-1 Satellite data at the instrument data processing level, with radiometric calibration applied only (Level-1A).
Level-2 Satellite data at the analysis-ready processing level, with atmospheric correction performed and geophysical variables calculated, but not yet spatially binned or mapped.
Level-2B Satellite data after custom merging atmospheric correction methods in the present study, not yet spatially binned or mapped.
Level-3 Satellite data that is analysis-ready and also has been spatially binned and mapped to a consistent map projection.
Ocean band Wavelength or spectral range at which a satellite sensor has been optimized to image dark water surfaces, for example, the 443 nm band on MODIS sensors.

Figure 1 .
Figure 1.Map of MODIS-Aqua satellite data extent and in situ validation stations.Long-term mean Rrs(645) at 250 m spatial resolution is mapped in color.White points indicate SeaBASS in situ data locations from 2005 to 2014.

Figure 2 .
Figure 2. Validation at 1 km spatial resolution of MODIS-Aqua Rrs(λ) with in situ Rrs(λ) observations, showing mean (points) and standard deviations (error bars) of Rrs(λ) validation data points derived from MODIS-Aqua (black dashed line) and measured in situ (green solid line) (n = 63), spanning multiple seasons and years 2005-2014.

Figure 4 .
Figure 4. Climatologies of mean monthly band ratios (Rrs(λ 1 )/Rrs(λ 2 )) over the 18-year time series.Error bars represent standard deviation over the years 2003-2020.By ratio type, climatologies include (a-c) red-to-green ratios, (d) red-to-blue ratio, and (e-h) green-to-blue ratios.Black lines indicate upper Bay mainstem climatologies and blue lines indicate lower Bay mainstem climatologies, defined as mainstem area north and south of 38°N latitude in Figure 1.

Figure 5 .
Figure 5. Trends over time in Rrs at single bands, from 2003 to 2020 for bands (a) 443 nm through (g) 645 nm according to linear least squares fits over all surface water areas with >80% of monthly images at the nominal spatial resolution of each band, i.e., (g) Rrs(645) at 250 m, (b) Rrs(469) and (f) Rrs(555) at 500 m, and all other bands at 1 km.Small black dots highlight significant trends (p < 0.1).Trends are expressed as relative trends normalized to the long-term mean Rrs at each location (% yr −1 ).Absolute values of median increasing and decreasing trends over time were 0.8%-1.0%yr −1 and 1.3%-2.9%yr −1 respectively (see Table4).

Figure 7 .
Figure 7. Climatology of trends over time 2003 to 2020 in single band Rrs.Trends are expressed as relative trends normalized to the long-term mean value (% yr −1 ).Black circles indicate mean trends for the Upper Bay and blue squares indicate mean trends for the lower Bay mainstem, as defined by mainstem waters separated by the latitude line 38°N in Figure 1.Filled points indicate strong trends (p < 0.1 for >25% of the area of the region).Shading indicates the standard deviation spatially for the trends in each region.Points above the zero line (outlined in red) indicate that Rrs for that month showed a long-term increase, while points below the zero line indicate a long-term decrease.

Figure 8 .
Figure 8. Climatology of trends over time 2003-2020 in selected band ratios (Rrs(λ 1 )/Rrs(λ 2 )).Trends are expressed as relative trends normalized to the long-term mean value (% yr −1 ).Black circles indicate mean trends for the upper Bay mainstem and blue squares indicate mean trends for the lower Bay mainstem, as defined by mainstem waters separated by the latitude line 38°N in Figure 1.Filled points indicate strong trends (p < 0.1 for >25% of the area of the region).Shading indicates the standard deviation spatially for the trends in each region.Points above the zero line (outlined in red) indicate that Rrs for that month showed a long-term increase, while points below the zero line indicate a long-term decrease.
Diffuse light attenuation coefficient of photosynthetically active radiation (PAR), a.k.a.light attenuation (m −1 ) optical Archive and Storage System (data repository)TSSTotal suspended solids, a.k.a.total suspended matter (TSM) or suspended particulate matter (SPM)

Table 2 MODIS
-Aqua Bands Analyzed in This Study 10.1029/2021JC017959 6 of 20 cloud cover, especially at 1 km spatial resolution.The regions termed upper Bay and lower Bay are defined as mainstem surface waters north and south of the latitude 38°N, respectively (Figure Paired data points using in situ observed Rrs(λ) versus corresponding daily (<6 hr) MODIS-Aqua pixels or pixel window averages.b MAE = Mean absolute error, RMSE = Root mean squared error, Mean APD = Mean absolute percent difference, and R = correlation coefficient (See Text S1 in Supporting Information S1 for calculations). a

Table 3
Validation of Satellite Rrs