Modeling Coastal Water Clarity Using Landsat-8 and Sentinel-2

Abstract

Understanding and attributing changes to water quality is essential to the study and management of coastal ecosystems and the ecological functions they sustain (e.g., primary productivity, predation, and submerged aquatic vegetation growth). However, describing patterns of water clarity—a key aspect of water quality—over meaningful scales in space and time is challenged by high spatial and temporal variability due to natural and anthropogenic processes. Regionally tuned satellite algorithms can provide a more complete understanding of coastal water clarity changes and drivers. In this study, we used open-access satellite data and low-cost in situ methods to improve estimates of water clarity in an optically complex coastal water body. Specifically, we created a remote sensing water clarity product by compiling Landsat-8 and Sentinel-2 reflectance data with long-term Secchi depth measurements at 12 sites over 8 years in a shallow turbid coastal lagoon system in Virginia, USA. Our satellite-based model explained ∼33% of the variation in in situ water clarity. Our approach increases the spatiotemporal coverage of in situ water clarity data and improves estimates from bio-optical algorithms that overpredicted water clarity. This could lead to a better understanding of water clarity changes and drivers to better predict how water quality will change in the future.

Introduction

Measuring and predicting water quality variability is valuable for understanding and managing coastal marine and estuarine ecosystems (Álvarez-Romero et al., 2014). For example, understanding changes in water quality over space and time can aid the restoration and conservation of submerged aquatic vegetation, which is often limited by light availability (Carr et al., 2010). The health of coastal oceans are also tied to human well-being, with about 40% of the global population living within 100 km of the coast (Maul & Duedall, 2019). Coastal water quality impacts primary productivity and marine life, with associated impacts on fishing and industry (Cloern et al., 2014). Eutrophication, caused by an excess of nutrients and measured by changes in water quality metrics, can lead to hypoxia, harmful algal blooms, and fish kills (Sullivan et al., 2021).

Water clarity—a key component of water quality—has historically been measured using Secchi disks: black and white painted disks that are lowered into the water to estimate the depth at which they are no longer visible from the surface (Secchi depth [Z_SD]; Preisendorfer, 1986). Secchi disks serve as a reliable, repeatable, and low-cost tool; thus, they have contributed to extensive long-term records of overall water quality (Lee et al., 2015). For example, they are used by citizen scientists (Luis et al., 2019) to greatly expand sampling coverage, and in national assessments to measure water clarity as an ecological indicator for eutrophication (Sullivan et al., 2021). However, due to logistical constraints, Z_SD and other in situ point measurements cannot fully capture broad-scale, long-term spatiotemporal variability in water clarity. Hence for most regions and time periods, there are large gaps in records of water clarity (Lee et al., 2016; Luis et al., 2019). Coupling remotely-sensed satellite measurements with in situ measurements could provide a more complete understanding of water quality changes and drivers in coastal ecosystems, as well as a way to separate directional trends from natural variability.

Z_SD and other biogeochemical parameters can be estimated by satellites measuring water-leaving radiances in the visible spectrum, referred to as “ocean color” (Lee et al., 2016; Werdell & McClain, 2019). Ocean color data can be helpful in studying coastal waters on a synoptic scale (>100 km²), but coastal environments pose unique challenges to ocean color remote sensing; thus, most studies to link in situ and remotely sensed water clarity have been based in open ocean environments (Concha & Schott, 2016; Morel & Bélanger, 2006). High concentrations of particulate organic matter and suspended sediments, proximity to land, backscattering from shallow waters, and bubbles from breaking waves affect bio-optical measurements (Loisel et al., 2013). Validation with in situ measurements is crucial before implementing ocean color algorithms in any water body, but the unique challenges posed by the coastal zone make it especially important to develop and validate such algorithms.

A second challenge posed by satellite estimation of coastal water clarity is that the spatial resolution of traditional ocean color satellite sensors (≤1 km) is too coarse to avoid interference by landforms while capturing spatial variation across coastal waters features (e.g., river mouths, ocean inlets, coastal lagoons). Fortunately, the Landsat-8's Operational Land Imager (OLI) and the Sentinel-2's MultiSpectral Instrument (MSI) have moderate-to-high spatial resolution (10–60 m) that enable observations in coastal environments (Franz et al., 2015; Pahlevan et al., 2019). The benefit of combining these observations has been recently demonstrated. For example, NASA's Harmonized Landsat Sentinel (HLS)-2 project shows that using the two satellites in unison increases data accuracy and temporal resolution (Masek et al., 2018).

However, studies using Landsat 8 and Sentinel-2 to measure water clarity have largely focused on lakes and rivers, highlighting the need to evaluate product consistency in coastal seas (Chen et al., 2020; Manzo et al., 2015; Page et al., 2019). The Lee et al. (2015) semi-analytical Z_SD model has been applied to various water bodies, but cross comparisons of OLI and MSI have been limited in coastal seas and require further investigation for effective usage of satellite imagery for water quality monitoring (Chen et al., 2019; Lee et al., 2016; Liu et al., 2019; Luis et al., 2019).

Thus, this work evaluates the accuracy and consistency of Landsat-8 and Sentinel-2 water clarity estimates in a coastal ocean lagoon system in Virginia, USA. Satellite estimations of Z_SD estimates were compared with corresponding in situ Z_SD from long-term in situ sampling efforts and these comparisons informed the development of regionally tuned Z_SD models for Landsat 8 and Sentinel 2. These results reveal that regionally tuned satellite models can improve the estimation of overall water quality in shallow, turbid coastal areas.

Methods Study System and In Situ Data Collection

We focused our investigation on a coastal lagoon system studied by the Virginia Coast Reserve (VCR) Long Term Ecological Research Project (VCR-LTER) located on the Eastern Shore of Virginia, USA, near the southern tip of the Delmarva Peninsula (Figure 1). The VCR is the largest expanse of undeveloped coastline along the U.S. Atlantic seaboard. Low nitrogen inputs and frequent exchange with the Atlantic Ocean via inlets between barrier islands (Figure 1) results in relatively good water quality in comparison with most temperate coastal bays and estuaries in the United States and worldwide (McGlathery et al., 2001). Due to low human impacts, the VCR serves as a model system for studying the long-term impacts of climate on temperate coastal lagoons. The VCR also serves as the site of a successful large-scale eelgrass restoration project, after a massive die-off in the early 1930s due to storms and disease (Reynolds et al., 2016). Water quality monitoring in the bay has informed restoration efforts and helped quantify positive environmental effects of restoration (Orth et al., 2020). Since 1992, VCR-LTER researchers have collected Z_SD and other water quality parameters at 17 sites that include tidal flats (0–2 m depth), deep flats (2–4 m depth), and deeper oceanic inlets and channels (>4 m depth) (McGlathery & Christian, 2020; Safak et al., 2015). Sampling was carried out monthly from 1998 to 2008, and quarterly from 2008 to 2021 (McGlathery & Christian, 2020). Temperature, salinity, and dissolved oxygen (converted to apparent oxygen utilization) were measured with a YSI Datasonde (6600 V2). Total suspended solids were estimated via dry weight difference after filtering water samples with Whatman GF/F filters (0.7 μm nominal pore size), chlorophyll-a and phaeopigments were determined spectrophotometrically, and 500 mL water samples were analyzed for nutrient concentrations on a Lachat Quick-Chem 8500 flow-through analyzer.

Figure 1. Map of the 17 tested sites with in situ ZSD data. Numbers identify each site. Green markers indicate that satellite data could be collected at the site using the appropriate quality flags, while red markers indicate sites without valid satellite data.

Satellite Data Processing and Algorithm Evaluation

To retrieve remote sensing reflectances (R_rs) for the calculation of satellite-derived Z_SD, Landsat-8 and Sentinel-2 Level-1 images were collected. The Sentinel-2A and -2B satellites have a combined 5-day repeat orbit, while Landsat-8 has a 16-day repeat orbit. The OLI on Landsat-8 collects data at a 30 m spatial resolution in four spectral bands centered on wavelengths 443, 482, 561, and 655 nm. The MSI on Sentinel-2 collects data at 60 m spatial resolution with a band centered on 443 nm and at 10 m resolution for bands centered on 490, 560, and 665 nm. USGS Earth Explorer was used (U.S. Geological Survey, 2022; https://earthexplorer.usgs.gov/, accessed June 2019 through July 2022) to collect Level-1 Collection-1 Landsat-8 images (https://doi.org/10.5066/F71835S6) and the Copernicus Open Access Hub (European Space Agency, 2022; https://scihub.copernicus.eu/, accessed June 2020 through July 2022) to collect Level-1 Collection-1 Sentinel-2 images. The Landsat-8 and Sentinel-2 satellite images used in this study are listed as Tables S1 and S2 in Supporting Information S1, respectively.

R_rs were generated with two atmospheric correction processors: NASA's standard Level-2 generator atmospheric correction algorithm (l2gen) in NASA SeaDAS 8.2 (Baith et al., 2001) and ACOLITE from the Royal Belgian Institute of Natural Sciences (Version 20220222.0) (Vanhellemont, 2019, 2020; Vanhellemont & Ruddick, 2018). In l2gen, the NASA standard NIR-SWIR algorithm with bands 5 and 7 (865 and 2,201 nm) were selected for atmospheric correction of coastal waters (Wei et al., 2018). The following default Level-2 quality flags were used: CLDICE (probable cloud or ice contamination), HILT (very high or saturated radiance), and STRAYLIGHT (straylight contamination). LAND (pixel is over land) was deselected due to the default land mask being too coarse for this area and the following reflectance thresholds were used to flag and remove land pixels: R_rs(655) < 0.01, R_rs (561) < 0.012, R_rs (482) < 0.001, and R_rs (443) < 0.001. The bidirectional reflectance distribution function of Morel et al. (2002) was implemented. In ACOLITE, the default Dark Spectrum Fitting algorithm was implemented (Vanhellemont, 2019, 2020; Vanhellemont & Ruddick, 2018) with default masks for land, negative reflectances, cirrus clouds, and high reflectance. Valid R_rs from l2gen and ACOLITE were recovered at 12 of 17 in situ sampling sites (five sites were too close to land): six ocean inlet sites, two lagoon sites, and three mainland tidal creek sites (Figure 1).

Z_SD was computed by deriving inherent and apparent optical properties from R_rs. The Quasi-Analytical Algorithm (QAA, version 6) was used to derive total absorption (a) and backscattering (b_b) coefficients from R_rs with the QAA (Lee et al., 2002), and these inherent optical property (IOP) coefficients were used to calculate the diffuse attenuation coefficient (K_d) with Lee et al., 2013. K_d(530) was determined empirically by the methods of Lee et al. (2016) to fill the large spectral gap between 482 and 561 nm. The minimum K_d value, the spectral band with the lowest attenuation, and the R_rs value at the corresponding wavelength ($${R}_{\mathrm{r}\mathrm{s}}^{\mathrm{t}\mathrm{r}}$$), were then used to compute the satellite derived Secchi depth (Z_SD,sat, m) (Equation 1) (Lee et al., 2015, 2016). [Image Omitted. See PDF]

Satellite imagery within ±0–1 day of in situ sampling were used (2 OLI images, 6 MSI images) and Z_SD,sat was found by averaging pixels in a 3 × 3 box (90 × 90 m for Landsat-8, 30 × 30 m for Sentinel-2) centered around the corresponding in situ site coordinate (Figure S1 in Supporting Information S1). The impact of anomalous data points on the spatial averaging within the 3 × 3 box was evaluated by also calculating the median of the 3 × 3 box. The mean and medians were not statistically different (Landsat-8: P = 0.5; Sentinel-2: P = 0.9). We did not filter for large variance due to limited match-ups, as is done in some studies using a coefficient of variation threshold (Dogliotti et al., 2011; Mélin, 2022). L2gen yielded 11 matchups from OLI and 41 matchups from MSI. ACOLITE yielded 3 OLI and 44 MSI matchups. The number of observations varied by site and date due to varying cloud cover and choice of quality flags. We removed matchups with recorded water depths within 0.5 m of the recorded in situ Secchi depths. The workflow is synthesized in Figure 2.

Figure 2. Workflow for obtaining ZSD with l2gen in NASA SeaDAS, bio-optical algorithms, and empirical adjustments. l2gen ultimately chosen over ACOLITE for processing (Section 3.1).

Cross Sensor Comparison

Landsat-8 and Sentinel-2 Secchi depths from l2gen and the Lee et al. (2016) algorithm (Z_SD,sat) (Equation 1) from the same day at the same location were compared. This approach yielded 5 clear-sky same-day images from 2021, from which we sampled 150 random coordinates using QGIS 3.14 (QGIS.org, 2020). We compared Landsat-8 Z_SD,sat and Sentinel-2 Z_SD,sat using a Type II ordinary least squares linear regression.

Modeling for Algorithm Adjustment

We used Type II ordinary least squares linear regression to predict in situ Secchi depths (Z_SD,insitu) from satellite Secchi depth estimates (Z_SD,sat) from l2gen and ACOLITE to determine which atmospheric correction software was most suitable. We also fit two individual linear models for Landsat-8 and Sentinel-2 data to determine whether Landsat-8 and Sentinel-2 differed in their predictions of in situ Z_SD and to determine the proportion of variance of in situ Z_SD explained by each satellite separately.

Models were fit in R 4.2.1 (R Core Team, 2022). We assessed the significance of model terms using F tests to determine if the model better fit the data than the Z_SD algorithm. We checked for homogeneity of variance by plotting normalized model residuals against model predictions and individual predictors. We ensured normality of residuals using histograms and quantile-quantile plots. We tested for temporal autocorrelation using autocorrelation function analysis; no significant autocorrelation was detected.

Model Assessment

We evaluated model skill by comparing modeled Secchi depths (Z_SD,model) to in situ Secchi depths (Z_SD,insitu) using the root mean square error (RMSE) and the mean absolute percent difference (MAPD) in R using the package Metrics 0.1.4 (Hamner & Frasco, 2018). [Image Omitted. See PDF] [Image Omitted. See PDF]

Due to the lack of a sufficient number of observations for an out-of-sample validation, we assessed the uncertainty in the model coefficients by calculating the standard errors and confidence intervals with bootstrapping (n = 10,000) (R package boot 1.3–28; Canty & Ripley, 2021; Davison & Hinkley, 1997).

Spatiotemporal Analysis

We compared spatial trends of water clarity maps generated from clear sky Landsat-8 and Sentinel-2 images taken on the same day (30 January 2021). We also generated monthly water clarity maps from clear sky Landsat-8 images to demonstrate changing spatial patterns with time. To determine the benefit of increased temporal coverage, we visually compared temporal patterns in time series of combined satellite model estimates (n = 396 observations; 2013–2021) and all in situ data (n = 306 observations; 1992–2021) at three water quality sites in a lagoon (site 2), mainland creek (site 6), and ocean inlet (site 16). Specifically, generalized additive models (GAMs; Hastie & Tibshirani, 1986) were used to estimate nonlinear trends as a function of date (interannual variation), year-day (seasonal variation), and their interaction (tensor product) in R using the package mgcv 1.8–35 (Wood, 2017). By checking k-indices and P values, the following number of basis functions were selected: (a) 15 for date, (b) 15 for year day, and (c) 10 for the interaction/tensor product. Thin plate regression splines, computationally-efficient splines used to estimate smooth functions of multiple predictors, were used to model interannual variation and the interaction. We used cyclic cubic splines to model seasonal variation to avoid a discontinuity between the end and beginning of the year (Wood, 2017). To investigate seasonal trends in in situ water quality data, we used locally weighted scatterplot smoothing fits.

Results Algorithm Evaluation

Satellite estimates overpredicted Z_SD relative to their corresponding in situ values by an average factor of about 2.4 for l2gen and 1.8 for ACOLITE. l2gen and ACOLITE Z_SD both predicted in situ Z_SD (l2gen: R²_adj = 32.5%, F_1,49 = 25.05, P < 0.001; ACOLITE: R²_adj = 22.0%, F_1,43 = 13.38, P < 0.001; Figure 3) and NASA l2gen yielded larger Z_SD from the Lee et al. (2016) algorithm (mean, $$\overline{x}$$ = 1.75 m, standard deviation, s = 0.46 m) than ACOLITE ($$\overline{x}$$ = 1.30 m, s = 0.24 m). However, ACOLITE still overpredicted Z_SD relative to the corresponding in situ values ($$\overline{x}$$ = 0.77 m, s = 0.26 m).

Figure 3. ZSD predictions from l2gen (a) were generally greater than those from ACOLITE (b), although both were positively correlated with in situ measurements. Dotted line shows 1:1 relationship.

Cross-Sensor Comparisons

Z_SD estimates from the three satellites were highly correlated with a linear relationship (R²_adj = 85.3%, F_1,545 = 3,170, P < 0.001). Sentinel-2 generally yielded higher Z_SD,sat than Landsat-8 (note that most points are below the 1:1 line in Figure 4).

Figure 4. Landsat-8 ZSD,sat plotted against Sentinel-2 ZSD,sat. Sentinel-2 generally yielded higher ZSD,sat than Landsat-8. Dashed 1:1 line as in Figure 3. Linear fit in red.

Model Statistics and Assessment

Spectral shapes were similar for both atmospheric correction processors, but ACOLITE systematically yielded larger R_rs than l2gen (Figure S2 in Supporting Information S1). Although ACOLITE yielded Z_SD,sat that were closer to the 1:1 line (Z_SD,model = 0.56 Z_SD,sat + 0.04; Figure 3a) than l2gen Z_SD,sat (Z_SD,model = 0.33 Z_SD,sat + 0.19; Figure 3b), R²_adj, model skill metrics, and standard errors were better for the l2gen-based model (Table 1).

Table 1 Comparison of In Situ and Satellite Derived Z_SD From l2gen (SeaDAS) and ACOLITE

	l2gen (SeaDAS)	ACOLITE
RMSE	0.21 m	0.24 m
MAPD	25%	27%
Mean, standard error, and 95% CIs of the slope	0.33 ± 0.07 [0.22, 0.49]	0.56 ± 0.15 [0.29, 0.91]
Mean, standard error, and 95% CIs of the intercept	0.19 ± 0.11 [−0.08, 0.38]	0.04 ± 0.19 [−0.39, 0.37]
R2adj	32.5%	22.0%

Note. For each atmospheric correction method, the RMSE, MAPD, Mean, standard error, and 95% CIs of the intercept, and R²_adj are reported.

Due to differences in Landsat-8 and Sentinel-2 Z_SD,sat estimates (Figure 4), the final algorithm adjustment involved separate satellite models (Equations 4 and 5). We also used l2gen atmospheric correction due to improved results over ACOLITE (Table 1). [Image Omitted. See PDF] [Image Omitted. See PDF]

The final model (Equations 4 and 5; Figure 5) improved estimates (Landsat-8: RMSE = 0.16 m, MAPD = 19%; Sentinel-2: RMSE = 0.22 m, MAPD = 26%; Table 2) relative to the Lee et al. (2016) model (Landsat-8: RMSE = 1.04 m, MAPD = 121%; Sentinel-2: RMSE = 1.05 m, MAPD = 148%).

Figure 5. ZSD,sat from Landsat-8 (a) and Sentinel-2 (b) plotted against ZSD,insitu. The single-satellite regression models (Equations 4 and 5) are plotted as solid lines with 95% confidence intervals plotted in gray.

Table 2 Model Skill Metrics, Model Uncertainties, and Fixed-Effects Results for Satellite Models

	Landsat-8 OLI (l2gen)	Sentinel-2 MSI (l2gen)
RMSE	0.16 m	0.22 m
MAPD	19%	26%
Mean, standard error, and 95% CIs of the slope	0.25 ± 0.092 [0.10, 0.48]	0.37 ± 0.094 [0.21, 0.58]
F, P value of the slope	F = 7.30, P = 0.024	F = 19.22, P < 0.001
Mean, standard error, and 95% CIs of the intercept	0.36 ± 0.16 [0.005, 0.64]	0.10 ± 0.16 [−0.25, 0.38]
F, P value of the intercept	F = 4.74, P = 0.06	F = 0.44, P = 0.51
Mean and standard error of R2adj	38.6% ± 22%	31.8% ± 14.5%

Model skill was similar among lagoon, mainland creek, and ocean inlet sites (Figure 6). There were no clear trends between water depth and model performance (Figure S3 and Table S3 in Supporting Information S1). Ocean inlet sites are most well-represented in the model, although the MSI model tends to underestimate high Z_SD at ocean inlet sites and overestimate low Z_SD (orange in Figure 6c).

Figure 6. ZSD,sat from Landsat-8 (green) and Sentinel-2 (orange) plotted against ZSD,insitu for each site type. Sites no. 8 and 16 are lagoons, sites no. 2, 13, and 14 are mainland creeks, and sites no. 3, 5, 6, 12, 15, and 17 are ocean inlets. Landsat-8 (green) and Sentinel-2 (orange) prediction lines are plotted. Mean absolute percent difference (MAPD) and root mean square error (RMSE) are reported for each region. Sample size is not equal between site types.

Spatiotemporal Analysis

Spatial trends within the system were similar in Landsat-8 and Sentinel-2 water clarity maps taken on the same day, providing confidence that the remote sensing data is capturing robust geophysical and ecological patterns. Flows from the inlets through the deeper channels have higher Z_SD (ca. >0.8 m) than the surrounding water (ca. 0.4–0.7 m), and plumes off-shore of inlets have shallower Z_SD (0.6–0.8 m) than further offshore (>1 m). Sentinel-2 MSI had a larger Z_SD range (ca. 0.3–1 m) than Landsat-8 OLI (ca. 0.4–0.9 m) and captured finer scale spatial variability (Figure 7). Satellite images from individual dates highlight that Secchi depths in this system are highly variable in both time and space (Figure 8).

Figure 7. Maps of model-corrected Landsat-8 derived ZSD (a) and Sentinel-2 derived ZSD (b) from 30 January 2021.

Figure 8. Maps of model-corrected Landsat-8 derived Secchi depths from monthly clear sky images (2019–2021).

Season cycles were significant at all three sites (Site 2: P < 0.001, Site 6: P = 0.002, Site 16: P < 0.001), as well as the tensor product explaining the interaction between interannual and seasonal variability (Site 2: P = 0.006, Site 6: P = 0.009, Site 16: P = 0.01; Figures 9a, 9c, and 9e). Sites 2 and 16 experienced the strongest seasonality and water clarity dips around July, corresponding to a dip in in situ Z_SD and peaks in in situ particulate inorganic matter, total suspended solids, apparent oxygen utilization, chlorophyll-a, phaeopigments, PO₄, and NH₄ (Figure 10). Interannual trends were not significant at any of the 3 sites (Site 2: P = 0.20, Site 6: P = 0.66, Site 16: P = 0.26; Figures 9b, 9d, and 9f).

Figure 9. ZSD,insitu values (blue) and ZSD,model values (green: Operational Land Imager (OLI) and orange: MultiSpectral Instrument (MSI)) plotted against time (top; a, c, and e) and day of the year (bottom; b, d, and f) in a mainland creek (site 2; a–b), lagoon (site 6; c–d), and ocean inlet (site 16; e–f). Black lines show thin-plate regression splines and cyclic cubic splines for interannual (top) and seasonal (bottom) patterns from the generalized additive model, respectively. The splines were calculated from in situ, OLI, and MSI data. The 95% confidence intervals are plotted in gray.

Figure 10. Seasonal trends in water quality parameters at 17 sites as shown by locally weighted scatterplot smoothing (LOWESS) fits. All parameters except Secchi depth and water temperature have been standardized. Panel A shows a peak in temperature (blue) corresponding to a dip in Secchi depth (black) mid-summer. Water temperature data are shown in light blue and Secchi depth data are shown in gray. Panels B–D show the seasonal variations of water quality parameters as modeled by LOWESS. Particulate organic matter denoted by POM, particulate inorganic matter denoted by PIM, total suspended solids denoted by TSS, and total dissolved nitrogen denoted by TDN. Secchi depths are plotted as well, as shown by the gray data points and black smooth. Data from 1992 to 2020 water quality sampling (McGlathery & Christian, 2020).

Discussion

Landsat-8 and Sentinel-2 models were developed to increase and improve water clarity observations in the VCR-LTER. Our results demonstrate that high-resolution satellite observations can enable the estimation of water clarity in estuaries and coastal seas across a range of spatiotemporal scales, but require sensor-specific calibration and validation with in situ measurements. We anticipate that our approach can be adapted to coastal waters broadly where environmental monitoring organizations are limited to in situ data and potentially biased satellite estimates. The potential for high frequency water clarity estimation is further enhanced because we found that Landsat-8 and Sentinel-2 reflectances and Secchi depths are highly comparable.

Satellite Overestimation of Z_SD

Biases introduced during atmospheric correction could have led to satellite Z_SD overestimation. The atmospheric correction methods implemented by l2gen explained more variation in in situ Z_SD than ACOLITE (Figure 3, Table 1). This result complements previous work showing that l2gen performs better than ACOLITE for turbid waters and complex, optically shallow coastal environments (Ilori et al., 2019; Vanhellemont, 2019, 2020; Wei et al., 2018). We also found that ACOLITE consistently yielded higher reflectance values than l2gen at all four wavelength bands (Figure S2 in Supporting Information S1). These results are consistent with the findings of Ilori et al. (2019) who compared different atmospheric correction methods over optically complex coastal waters using in situ radiometric measurements from the Aerosol Robotic Network—Ocean Color (AERONET-OC). Future work could investigate other atmospheric correction algorithms such as Case 2 Regional CoastColour (C2RCC) and POLYnomial-based algorithm applied to Medium Resolution Imaging Spectrometer (POLYMER) in this water body as was done in the Chesapeake Bay (Windle et al., 2022).

Validation of remote sensing products with in situ radiometric observations would be useful in regionally adjusting the algorithm. However, using an existing long-term data set for regional adjustment, as demonstrated in this paper, has its own advantages compared to this approach. Long term datasets have repeated seasonal sampling, so statistical models are less biased by the timing of in situ sampling. Capturing a wide range of variability in in situ sampling is especially important in systems that experience dynamic change across seasons. Using existing datasets is also cost effective compared to radiometric calibration. Another source of bias could have been introduced during the derivation of IOPs. It has been found that the QAA can lead to underestimation of IOPs in turbid waters (Yang et al., 2014), which ultimately leads to an overestimation of Z_SD. Alternative R_rs to IOP algorithms, such as the “QAA turbid” (Yang et al., 2014), or alternative Z_SD algorithms could possibly yield more accurate Z_SD values.

Using Sentinel-2 Data in Conjunction With Landsat-8 Data

The strong relationship between Landsat-8 and Sentinel-2 measurements (Figure 4) and similarities in spatial water clarity patterns (Figure 7) demonstrate the compatibility of their data products for estimating water clarity. We found that Sentinel-2 consistently yields larger Z_SD values (corresponding to lower R_rs values) than Landsat-8, suggesting that differences are most likely due to inherent satellite product differences rather than environmental factors (e.g., tidal differences occurring in the temporal window between satellite overpasses).

We decided to fit two separate models for each satellite due to differences in satellite estimates and to preserve the highest spatial resolution of Sentinel-2 MSI. However, care must be taken when comparing water clarity maps with different spatial resolutions. Down-sampling Sentinel-2 MSI to the spatial resolution of Landsat-8 OLI may be necessary for certain spatial analyses. Although down-sampling did not affect match-ups used in modeling (Figure S1 in Supporting Information S1), differences are significant when models are applied to entire satellite images (Figure 7).

Another consideration regarding merging Landsat and Sentinel data is parametrizing NASA's HLS products, surface reflectances with 30 m spatial resolutions and 5 days temporal resolutions (Masek et al., 2018). Although developed primarily for terrestrial applications, HLS data have been used in water quality modeling to increase temporal resolution while minimizing errors associated with satellite product differences (Peterson et al., 2020). The workflow for harmonizing products includes common atmospheric correction, spatial co-registration, normalization of the solar and view angles, and adjustment for differences in wavelength bands (Claverie et al., 2018). This procedure corrects for Landsat-8/Sentinel-2 differences that could have contributed to the mismatch between Landsat-8 and Sentinel-2 reflectances and Secchi depths in our study.

Data Collection and Model Limitations

Our work highlights the need for in situ bio-optical and water quality measurements to fully leverage satellite imagery as an effective and reliable water quality monitoring tool, especially in dynamic waters. Increasing the number of in situ match-ups, especially to cover a large dynamic range, will strengthen the predictive power of statistical water clarity models. In this study, in situ sites were limited to the area shown in Figure 1. Application of the satellite model (adjustment to standard water clarity algorithm) outside of the VCR-LTER region and conditions introduces potential uncertainties. This could be addressed in future work by expanding the spread of in situ evaluation sites to quantify the accuracy of the model outside this region.

Although both satellites capture similar spatial patterns (Figure 7), the Sentinel-2 model captures higher Secchi depth values than the Landsat-8 model. We had limited in situ data for Landsat-8 adjustment due to its infrequent overpass relative to Sentinel-2. The discrepancy could be due to the inclusion of high (>1 m) in situ Secchi depths in the Sentinel-2 model. Although we captured a sufficient range of Secchi depths for a statistical fit, the Landsat-8 model would be strengthened by adding higher in situ Secchi depth values (>1 m). Coordinating field surveys with satellite overpasses would strengthen future remote sensing studies and provide more matchups for empirical algorithm adjustments.

Additionally, further investigation into the diurnal variability introduced by tidal forcing is important for interpreting water clarity (Shi et al., 2013). Future work could expand on our model to include tidal effects, which are different at each site. The residence times in this system range from weeks near the mainland to hours near the inlets (Safak et al., 2015), so tides more significantly affect inlet match-ups. We recommend that the current model be used to map spatiotemporal patterns and changes in water clarity across Landsat and Sentinel sensors.

Coupled In Situ/Satellite Spatiotemporal Analysis

Spatial variation in water clarity as captured by MSI and OLI (Figures 7 and 8) is affected by aquatic vegetation, salt marshes, tides, winds, bathymetry, and other factors. Submerged aquatic vegetation decreases turbidity by enhancing sediment deposition, and salt marshes (masked in white) slow down exchange with the ocean (Nardin et al., 2018). Strong winds can lead to increased mixing in the lagoon and decreased lagoon/ocean exchange (Safak et al., 2015). Patterns in water clarity can be paired with other data (e.g., bathymetry, weather data, hydrology) to better understand drivers of water clarity.

Satellite Z_SD maps (Figure 8) and in situ water quality data (Figure 10) demonstrate strong temporal variability in Secchi depth. As a highly seasonal system, we expect water clarity to change over the course of the year, but trends shown by in situ or satellite data alone may be affected by unique biases introduced by each approach. Field measurements are rarely taken during stormy winter weather when the water is more mixed and turbid, so in situ winter values may be biased toward higher measurements. In situ Z_SD depth measurements can also be biased by reduced visibility from waves, cloud cover, and sun position, as well as observer error (Pitarch, 2020). Satellite measurements are also imperfect, being affected by adjacency effects, proximity to land, seafloor backscatter, cloud cover, and aliasing due to whitecaps. There are also spatial limitations in satellite data retrievals, as shown by the inability to obtain satellite data too close to land (ca. <100 m away) or in areas affected by sunglint.

Coupling in situ Z_SD with satellite Z_SD can help alleviate these biases. For example, satellite data in our temporal analysis of three representative water quality sites may have alleviated the wintertime bias in in situ measurements. Additionally, the timing of field sampling within a season changed among years, which may have led to seasonal trends being confounded with interannual variability. Satellite data from 2013 to 2021 helped fill these gaps and create a more even spread of data across the year (Figures 9b, 9d, and 9f). An area for future work would be to conduct a rigorous analysis of Secchi depth seasonality using paired in situ and satellite observations. Interannual in situ Z_SD were also scarce at all three sites from 2013 to 2021 before satellite data was added. The combined time series showed no evidence of interannual changes, although relative constancy in water clarity may soon be changing (Parker & Crichton, 2011). This provides opportunities to further evaluate the model's predictive power with future in situ matchups. Rapid, accelerating sea level rise (Sallenger et al., 2012) and storm intensification (Fenster & Hayden, 2007; Hayden et al., 1991) may increase coastal erosion and increase turbidity (Turner et al., 2021). Likewise, increasing frequency, duration, and intensity of marine heatwaves threaten recently restored seagrass meadows, which have stabilized sediments and reduced local turbidity and chlorophyll concentrations (Aoki et al., 2020; Orth et al., 2012; Ummenhofer & Meehl, 2017). Rising water temperatures could also affect numerous parameters relevant to water clarity such as phytoplankton concentrations (Pesce et al., 2018; Trombetta et al., 2019). It is vital to couple in situ observations with satellite observations over the long term to understand and derive changes in water clarity and separate directional trends from natural variability.

Conclusion

We developed a Landsat-8/Sentinel-2 remote sensing model to estimate water clarity in an optically complex coastal water body. The application of this model increases the spatiotemporal resolution of water clarity estimation, addresses algorithm overestimates of water clarity for specific localities, and decreases errors associated with Landsat-8/Sentinel-2. We believe our approach can be implemented in dynamic coastal water bodies with limited in situ measurements; for example, as part of routine water quality monitoring. Coupling accurate satellite estimates with in situ observations over the long term is crucial to understanding coastal water clarity variability and its underlying physical and biological drivers. This understanding could help improve water clarity predictions and lead to the better management of coastal ecosystems.

Acknowledgments

We thank the numerous scientists, technicians, and students involved over the years in the water quality sampling and analysis component of the Virginia Coast Reserve LTER program. Early elements of this study were conducted by Jacob Bushey as part of undergraduate research at the University of Virginia. This work was funded by The Virginia Space Grant Consortium (award GG12317) with additional support from the National Science Foundation (NSF award 1832221).

Data Availability Statement

Datasets and software used for the analysis of in situ and satellite data were archived with the Environmental Data Portal (EDI) via https://doi.org/10.6073/pasta/fe66683665a0133b2d831e552ecbae10 under a CC-BY Attribution license and are available in this in-text citation reference: Lang et al. (2022). Satellite data was processed by NASA SeaDAS 8.2 (Baith et al., 2001; National Aeronatics and Space Administration and Ocean Biology Processing Group, 2022), available for download at https://seadas.gsfc.nasa.gov/downloads/ under the GNU General Public License (GPL), and ACOLITE from the Royal Belgian Institute of Natural Sciences (Version 20220222.0) (Royal Belgian Institute of Natural Sciences and The Remote Sensing and Ecosystem Modelling Team, 2022; Vanhellemont, 2019, 2020; Vanhellemont & Ruddick, 2018), available for download at https://odnature.naturalsciences.be/remsem/software-and-data/acolite under the GNU General Public License v2. Level-1 Collection 1 Landsat-8 were downloaded from USGS Earth Explorer, available here: https://earthexplorer.usgs.gov (U.S. Geological Survey, 2022) in compliance with U.S. Public Domain. Level-1 images used in this study can be downloaded using the bounding box [37.6501, −75.7864], [37.6077, −75.4582], [37.0596, −75.7809], [37.1570, −76.0776], date range 1 January 2013 to 30 July 2022, and Landsat Level-1 Collection-1 under Data Sets. Level-1 Collection-1 Sentinel-2 images are available to download from the Copernicus Open Access Hub (European Space Agency, 2022), available here: https://scihub.copernicus.eu/ and subject to the Legal Notice on the use of Copernicus Sentinel Data and Service Information: https://sentinels.copernicus.eu/documents/247904/690755/Sentinel_Data_Legal_Notice. Level-1 data used in this study can be downloaded by selecting the appropriate bounding box around the VCR peninsula (see above coordinates), selecting Sentinel-2 mission, and selecting S2MSI1C as the Product Type. Processed Level-2 imagery available at Lang et al., 2022 data set described above. In situ water quality data is available at this in-text citation reference: McGlathery and Christian (2020) and at http://www.vcrlter.virginia.edu/cgi-bin/showDataset.cgi?docid=knb-lter-vcr.247 under a CC-BY Attribution license. Random sampling of coordinate points was done in QGIS 3.14 (QGIS.org, 2020) at https://www.qgis.org/en/site/forusers/download.html and under the GNU General Public License (GPL). Figures and statistics were done in R 4.2.1 (R Core Team, 2022) at https://www.qgis.org/en/site/forusers/download.html and under the GNU General Public License v2. Satellite algorithms were applied in MATLAB (The MathWorks, Inc., 2022) at https://www.mathworks.com/products/matlab.html.