An estuarine-tuned quasi-analytical algorithm (QAA-V): assessment and application to satellite estimates of SPM in Galveston Bay following Hurricane Harvey

The standard quasi-analytical algorithm (Lee et al., 2002) was tuned as QAA-V using a suite of synthetic data and in situ measurements to improve its performance in optically complex and shallow estuarine waters. Two modifications were applied to the standard QAA: (1) the semi-analytical relationship for obtaining remote sensing reflectance just below the water surface as a function of absorption and backscattering coefficients was updated using Hydrolight^® simulations, and (2) an empirical model of the total non-water absorption coefficient was proposed using a ratio of green to red bands of an ocean color sensor, which is known to work well in various inland and estuarine environments. The QAA-V-derived total absorption and backscattering coefficients, which were evaluated in a variety of waters ranging from highly absorbing and turbid to relatively clear shelf waters, showed satisfactory performance on a Hydrolight-simulated synthetic dataset (R² > 0.87, MRE < 17 %), an in situ estuarine and nearshore dataset (R² > 0.70, MRE < 35 %), and the NOMAD (R² > 0.90, MRE < 30 %). When compared to the standard QAA (QAA-v6), the QAA-V showed an obvious improvement with ∼ 30–40 % reduction in absolute mean relative error for the Hydrolight-simulated synthetic and in situ estuarine and nearshore datasets, respectively. The methodology of tuning QAA was applied to the VIIRS ocean color sensor and validation results suggest that the proposed methodology can also be applied to other ocean color and land-observing sensors. The QAA-V was also assessed on VIIRS imagery using a regional relationship between suspended particulate matter (SPM) and particulate backscattering coefficient at 532 nm (b_btnw532; R² = 0.89, N = 33). As a case study, the QAA-V processing chain and VIIRS imagery were used to generate a sequence of SPM maps of Galveston Bay, Texas following the unprecedented flooding of Houston and the surrounding regions due to Hurricane Harvey in August 2017. The record discharge of floodwaters through two major rivers into the bay resulted in very high SPM concentrations over several days throughout the bay, with wind forcing additionally influencing its distribution into the coastal waters of the northern Gulf of Mexico. The promising results of this study suggest that the application of QAA-V to various ocean color and land-observing satellite imagery could be used to assess the bio-optical state and water quality dynamics in a variety of coastal systems around the world.

1 Introduction

Urbanization and the associated anthropogenic stressors are of major concern for the ecosystem health and water quality of estuarine environments, cumulatively affecting the coastal and marine ecosystems through estuarine–shelf exchange processes (Haynes et al., 2007; Bricker et al., 2008; Jutterström et al., 2014). Inherent optical properties (IOPs) such as absorption and backscattering coefficients have immense potential to capture changes in the bio-optical state of an aquatic system and hence provide crucial information about regional alterations in water quality associated with terrestrial pollution (Zielinski et al., 2009), harmful algal blooms (Hu et al., 2008), floods (Álvarez-Romero et al., 2013), hurricanes (Lohrenz et al., 2008; Chen et al., 2009), seasonal cycles (D'Sa and Miller, 2003; D'Sa et al., 2006; Singh et al., 2010; Joshi and D'Sa, 2015), and even human-induced catastrophes such as oil spills (Ramsey III et al., 2011; D'Sa et al., 2016). In addition to water itself, there are three major water constituents that contribute to water-leaving radiance (L_w), namely colored dissolved organic matter (CDOM; also called “gelbstoff” or “gilvin”), suspended sediments (detritus and minerals), and phytoplankton (Naik et al., 2011). CDOM and suspended sediments are strongly associated with light absorption in the blue, and therefore high concentrations may reduce light quality for photosynthetic organisms (e.g., phytoplankton and submerged vegetation) in estuarine waters (Keith et al., 2002; Ralph et al., 2007; Pedersen et al., 2012). The effects of water turbidity caused by dissolved and particulate components on physical and behavioral changes in aquatic species have been well reported in the literature (Wang et al., 2008; Kjelland et al., 2015). Collectively, these water constituents attenuate incoming light, while a fraction of it is backscattered out of water by the water itself and particles. Therefore, deciphering L_w (or remote sensing reflectance; R_rs) to separate the individual contributions of optically active components may provide crucial information about the bio-optical state of a water body.

Field-based sampling methods are traditional and accurate ways to measure bio-optical properties; however, they lack adequate spatial and temporal coverage for capturing short-period bio-optical alterations and estuarine-scale dynamics. In contrast, remote sensing platforms (e.g., satellite sensors) sense L_w (or R_rs) signal and are advantageous over field observations in providing better synoptic spatiotemporal coverage if the signal is successfully linked to in-water IOPs. As such, satellite-based remote sensing has been widely used to monitor harmful algal blooms (Carvalho et al., 2011; Hu et al., 2016), pollution events (Mishra et al., 2013; Zhao et al., 2014), suspended sediment dynamics (D'Sa et al., 2007), CDOM distribution and carbon flux (Joshi and D'Sa, 2015; Joshi et al., 2017a), and phytoplankton biomass and primary production (Uitz et al., 2010; Matsumoto et al., 2014), as well as to evaluate the effects of climate change on exotic marine biota (Liu et al., 2006; Castillo and Lima, 2010; Cavanaugh et al., 2011; Pu and Bell, 2017).

In recent years, both empirical and semi-analytical models have been frequently used to link satellite observations and in-water properties, such as IOPs, vertical diffuse attenuation coefficients (K_d), suspended particulate matter (SPM) concentrations, CDOM, pigment concentrations, phytoplankton cell counts and cell size, and particle size (D'Sa et al., 2003, 2006, 2007; Pan et al., 2010; Chen et al., 2013; Brewin et al., 2015; Joshi et al., 2017b). Empirical relationships are mathematical formulations (e.g., simple or multiple regressions) that directly link water-leaving measurements to the parameter of interest in surface waters. They are simple in nature, easy to implement, and do not require deep understanding of the underlying relationships between light and water properties. Because the performance of empirical relationships is uncertain outside the range of observations that are used to develop them, their applicability is doubtful and may cause significant errors if used in waters with different optical properties. In contrast, semi-analytical models, which are based on radiative transfer theory, invert R_rs using a suite of analytical and empirical relationships to derive the water IOPs (absorption and backscattering coefficients of water constituents; Lee et al., 2002). Because they solely depend on water-leaving radiance and require less information about in-water bio-optical properties, they have better applicability and accuracy compared to empirical methods in a variety of waters (IOCCG, 2006). However, a major drawback is that the retrieval of IOPs for individual water components is strongly dependent on the performance of the respective empirical models.

To some extent, this drawback is minimized in a multiband quasi-analytical algorithm (QAA) for optically deep waters (Lee et al., 2002). This algorithm analytically decomposes total non-water absorption coefficients (a_tnw) for combined CDOM and suspended sediments (a_dg) and phytoplankton (a_ϕ) using their spectral information. The QAA has been improved (e.g., QAA-v5 and QAA-v6) for better performance in turbid coastal waters (Lee et al., 2009). Several studies have also contributed to the standard QAA with regional and/or global modifications in a variety of waters such as the turbid waters of the Mississippi and Atchafalaya River system (Zhu et al., 2011), coastal waters of the South China Sea (Dong et al., 2013), and inland waters of the USA and China (Li et al., 2013). A large number of studies have evaluated the standard QAA in different regions (e.g., turbid inland waters of northeast China, shallow ponds of the northwestern Mississippi, Lake Taihu, Yellow Sea, East China Sea, Arctic, and low-latitude oceans) with acceptable performance in coastal and oceanic waters, but decreasing accuracy towards CDOM-rich and sediment-rich estuarine and inland waters (Lee et al., 2010; Qing et al., 2011; Zhu et al., 2011; Mishra et al., 2013; Mitchell et al., 2014; Zheng et al., 2014; Pitarch et al., 2016). Several factors could be responsible for the QAA's poor performance in shallow waters, including the following: (1) the empirical relationships of QAA were designed using field observations in coastal and oceanic environments and may not be suitable for the optically deep estuarine and nearshore waters; and (2) a majority of empirical models of the standard QAA use R_rs at blue wavelengths (e.g., 443 and 490 nm). However, it is well known that satellite products suffer from large errors at short wavelengths due to uncertainties in atmospheric correction, especially in coastal waters. This suggests the need for an estuarine-specific tuning of the QAA and its evaluation and application to newer ocean color satellite sensors.

In this study, we present a tuned multiband quasi-analytical algorithm (QAA-V) that is optimized primarily for the Visible and Infrared Imaging Radiometer Suite (VIIRS) ocean color sensor and calibrated for various other ocean color sensors, such as Sentinel-3 Ocean and Land Colour Instrument (Sentinel-3 OLCI), MODerate resolution Imaging Spectroradiometer (MODIS Aqua), MEdium Resolution Imaging Spectroradiometer (MERIS), Sea-viewing Wide Field-of-view Sensor (SeaWiFS), and land-observing sensors, such as Landsat 8 Operational Land Imager (Landsat 8 OLI) and Sentinel-2 Multispectral Instrument (Sentinel-2 MSI), to estimate IOPs in shallow estuarine and nearshore waters. First, synthetic data were generated using Hydrolight^® simulations (Mobley and Sundman, 2013) for highly absorbing and scattering waters and used collectively with estuarine in situ observations to update coefficients for the semi-analytical and empirical models of the standard QAA processing chain. The algorithm's performance is then evaluated on three datasets: (1) a Hydrolight-simulated dataset, (2) a subset of the well-known NASA bio-Optical Marine Algorithm Dataset (NOMAD), and (3) field observations that were obtained in various estuaries on the US East Coast and in the Gulf of Mexico. The QAA-V performance was compared to the QAA-v6, which was mainly tuned to improve QAA's performance in turbid coastal waters. Additionally, using a linear backscattering–SPM relationship, the QAA-V's applicability to VIIRS and estuarine waters is analyzed for SPM of various coastal sites, including Galveston Bay (USA). Finally, as a case study, VIIRS-derived SPM imagery of Galveston Bay was obtained following Hurricane Harvey to assess SPM dynamics in the bay and their impact on the coastal ocean.

Table 1 In situ estuarine and nearshore dataset (IES: IOPs and R_rs matchups) compiled from the SeaBASS with locations, possessing status, project instructors, and the purpose in this analysis. Apalachicola Bay, Barataria Bay, and Galveston Bay datasets were collected by the authors during various field surveys.

2 Materials and methods

2.1 Data for QAA-V

Three datasets were used for tuning and evaluating QAA-V's performance in a variety of waters ranging from highly turbid estuarine environments to relatively clear shelf waters. These datasets include (1) a synthetic dataset, (2) NASA's bio-optical marine algorithm dataset, and (3) an estuarine dataset. The availability of numerous observations representing the true state of natural systems is the primary requirement for any algorithm development and validation analysis. The Hydrolight^® radiative transfer model (Mobley and Sundman, 2013) was used to generate a large set of synthetic data (HL; N=561) for tuning and extending QAA-V's ability to perform in highly absorbing and highly scattering waters (e.g., turbid estuarine environments). The process of generating synthetic data using Hydrolight^® simulations was similar to the International Ocean-Colour Coordinating Group report (IOCCG Report 5; IOCCG, 2006) and is briefly described in Sect. S1 in the Supplement with necessary modifications based on in situ estuarine observations. The in situ estuarine and nearshore dataset (IES) included 340 concurrent water inherent optical properties (IOPs; e.g., absorption and backscattering coefficients) and above-water R_rs measurements at various locations on the US East Coast and in the northern Gulf of Mexico (Fig. 1a). Data were compiled from NASA's SeaBASS repository by applying a depth threshold of 10 m for obtaining measurements in estuarine and nearshore waters (Werdell et al., 2003). The IES dataset was further divided into a training set (EcoHAB and Tampa Bay; N=121) and a testing set (N=219) for tuning and validating QAA-V, respectively (Table 1). NOMAD (NASA bio-Optical Marine Algorithm Dataset) is a freely available, high-quality field dataset for ocean color algorithm development and validation (Werdell and Bailey, 2005). It includes IOPs and R_rs collected in waters ranging from oceanic to estuarine environments, but mostly in shelf waters around the world. We extracted a subset (N=547) containing complete observations of IOPs and R_rs (Fig. 1b). Data distributions of synthetic data clearly showed the representation of CDOM-rich and sediment-rich waters, whereas phytoplankton absorption was of secondary importance as generally observed in several estuarine environments (Fig. 2). The training and testing data (HL and IES datasets) ranged from approximately 0.1 to 7 m⁻¹ for the CDOM absorption coefficient (a_g443), 0.05 to 4.5 m⁻¹ (detritus + minerals or non-algal particles) for the absorption coefficient (a_NAP443), 0.05 to 2 m⁻¹ for the phytoplankton absorption coefficient (a_ϕ443), and 0.04 to 0.2 m⁻¹ for the particle backscattering coefficient (b_btnw532; Fig. 2).

Figure 1 (a) In situ dataset (IES) representing estuarine and nearshore waters (< 10 m) on the US East Coast and in the northern Gulf of Mexico (N=340), (b) NOMAD (N=547), and (c) Galveston Bay, Texas (USA). GIWW: Gulf Intracoastal Waterway.

Figure 2 Data statistics and distribution of water IOPs, (a) a_g443, (b) a_NAP443, (c) a_ϕ443, and (d) b_btnw532 for synthetic data (HL, Hydrolight^®; blue), in situ estuarine and nearshore data (IES; green), and the NOMAD (grey). The range of QAA-V indicates data that are used to update QAA-V for shallow waters (e.g., HL and IES).

2.2 Data for Galveston Bay

Galveston Bay, the seventh-largest estuary in the United States (area = ∼ 1600 km²; mean depth = ∼ 2 m), is located along the upper coast of Texas in the northern Gulf of Mexico (Fig. 1c). The Trinity River is the major source of fresh water (∼ 50 %) to the bay followed by the San Jacinto River (∼ 30 %) and local watersheds (∼ 20 %; Guthrie et al., 2012; Lucena and Lee, 2017). With the busiest petrochemical port in the US, Galveston Bay experiences frequent oil spills; ∼ 3500 oil spill (∼ 416 000 gallons) incidences were reported between 1998 and 2009 (Lester and Gonzalez, 2011). The bay is connected to the Gulf of Mexico via three passes: Bolivar Roads Pass, Rollover Pass, and San Luis Pass. Galveston Bay can be divided into four sections: (1) Trinity Bay (TB), (2) upper Galveston Bay (UGB), (3) lower Galveston Bay (LGB), and (4) East Bay (EB; Fig. 1c).

Surface water samples were collected at several stations during two field surveys on 29 September 2017 and 29–30 October 2017 as part of a larger study to investigate the aftereffects of Hurricane Harvey (25–29 August 2017) on the water quality of Galveston Bay. Available measurements of suspended particulate matter (SPM) concentration were utilized for evaluating the applicability of QAA-V in estuarine environments. Samples were filtered using precombusted and pre-weighed 47 mm, 0.7 µm porosity Whatman GF/F filters for SPM concentrations (Neukermans et al., 2012). An analytical scale with an accuracy of ±0.1 mg was used to measure the mass of SPM. Profiles of b_btnw were obtained at each station using the WETLabs VSF-3 (470, 530, 670 nm) and ECO BB (532 nm) backscattering sensors (D'Sa et al., 2006), and surface values were averaged for depth < 1 m. Above-water measurements of water-surface radiance (L_w, W m⁻² nm⁻¹ sr⁻¹, nadir = 40–50^∘, azimuth = 90–135^∘), sky radiance (L_sky, W m⁻² nm⁻¹ sr⁻¹, zenith = 40–50^∘, azimuth = 90–135^∘), and reference plate radiance (L_plate, W m⁻² nm⁻¹ sr⁻¹, nadir = 0^∘, azimuth = 90–135^∘) were collected using a GER1500 512iHR spectroradiometer under clear-sky conditions (Mobley, 1999). The spectroradiometer was set to provide an average of three internal scans for considering the variability in reference and target conditions. Hence, the final spectrum was an average of nine spectra (three replicates with three internal scans per measurement) at each station (Joshi et al., 2017a). The above-water remote sensing reflectance ($R_{\mathrm{rs}}^{\mathrm{0}+}$, unit: sr⁻¹) was obtained using the following equations (Mueller et al., 2003):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{Downwelling irradiance}\left(E_{\mathrm{d}}\right)\left(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}{\mathrm{sr}}^{-\mathrm{1}}\right)\\ \text{(1)}& {\displaystyle }& {\displaystyle }=\mathit{\pi }\times {\displaystyle \frac{L_{\mathrm{plate}}}{{\mathit{\rho }}_{\mathrm{plate}}}},{\displaystyle }& {\displaystyle }\text{Above-water remote sensing reflectance}\left(R_{\mathrm{rs}}^{\mathrm{0}+}\right)\left({\mathrm{sr}}^{-\mathrm{1}}\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{L_{\mathrm{w}}-\mathit{\rho }\times L_{\mathrm{sky}}}{E_{\mathrm{d}}}}-R_{\mathrm{rs}}\left(\mathrm{residual}\right),\end{array}$$

where ρ_plate is reference plate reflectance (99 %) and R_rs(residual) is attributed to residual sky radiance, which was taken as $R_{\mathrm{rs}}^{\mathrm{0}+}$(950 nm; Mobley, 1999).

Table 2 Processing steps of QAA-V for obtaining total absorption (a_t) and backscattering coefficients (b_bt). Levels 0 and 1C were adopted from Lee et al. (2002), whereas 1A and 1B were modified in this study^*. Level 1C1 is adopted from D'Sa et al. (2007).

^* Parameters g₀=0.0788 and g₁=0.2379 were derived with HL datasets. λ₀=551 or 555 nm; $R_{\mathrm{rs}}^{\mathrm{0}-}$: remote sensing reflectance just below water surface; a_tnw: total non-water absorption coefficient; b_btnw: total non-water backscattering coefficient; a_w: water absorption coefficient; b_bw: water backscattering coefficient; η: power-law exponent (D'Sa et al., 2007); SPM: suspended particulate matter concentration. Note: in the SPM model comparison, input backscattering values are obtained from QAA-V, whereas surface reflectance (ρ) is obtained by multiplying π and above-surface remote sensing reflectance ($R_{\mathrm{rs}}^{\mathrm{0}+}$).

In addition, discharge data were acquired from USGS river gauge sites for the Trinity River (Romayor site USGS 08066500 and Wallisville site USGS 08067252) and at the San Jacinto River (the eastern flank USGS 8070200 and the western flank USGS 8068090) for examining variations in freshwater flows after Hurricane Harvey. Wind speed and direction were acquired from NOAA Eagle Point station (ID-8771013) in Galveston Bay. Level-1 VIIRS (nine images) and MODIS Aqua (one image) products were obtained from NASA's Ocean Color data archive (OBPG, NASA), including three images during field surveys in Galveston Bay. VIIRS imagery was not available on 29 September 2017 corresponding to the first field survey; hence, the next-day image (30 September 2017) was used in this analysis. The Sentinel-3 OLCI Level-2 image was downloaded from the Earth Observation Portal (EUMETSAT) and the Landsat 8 OLI Level-1 image was downloaded from the USGS Earth Explorer for 29 October 2017.

2.3 QAA processing chain

The underlying structure of QAA-V is similar to the standard quasi-analytical algorithm (QAA; Lee et al., 2002) and its modifications (e.g., QAA-v5 and QAA-v6) suitable to coastal and open oceans (Lee et al., 2007, 2009). The processing pathway of QAA-V is illustrated in Table 2 and briefly mentioned here with justifications for necessary modifications. The QAA and its updated versions rely upon the principle that the spectral remote sensing reflectance just below the water surface ($R_{\mathrm{rs}}^{\mathrm{0}-}$) is a function of the spectral backscattering and absorption coefficients (Gordon et al., 1988), and it can be modeled using the following equation (Table 2, Level 1A):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{R}_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathit{\lambda }\right)\left({\mathrm{sr}}^{-\mathrm{1}}\right)=g_{\mathrm{0}}\times u\left(\mathit{\lambda }\right)+g_{\mathrm{1}}\times {\left[u\left(\mathit{\lambda }\right)\right]}^{\mathrm{2}},\\ \text{(3)}& {\displaystyle }& {\displaystyle }u\left(\mathit{\lambda }\right)={\displaystyle \frac{b_{\mathrm{bt}}\left(\mathit{\lambda }\right)}{a_{\mathrm{t}}\left(\mathit{\lambda }\right)+b_{\mathrm{bt}}\left(\mathit{\lambda }\right)}},\end{array}$$

where a_t and b_bt are the total absorption coefficient and total backscattering coefficient (m⁻¹), respectively. $R_{\mathrm{rs}}^{\mathrm{0}-}$ can be easily computed from above-surface remote sensing reflectance ($R_{\mathrm{rs}}^{\mathrm{0}+})$ using the following relationship (Lee et al., 1999; Table 2, Level 0).

$$\begin{array}{}\text{(4)}& {\displaystyle }{R}_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathit{\lambda }\right)\left({\mathrm{sr}}^{-\mathrm{1}}\right)={\displaystyle \frac{R_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathit{\lambda }\right)}{(\mathrm{0.52}+\mathrm{1.7}\times R_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathit{\lambda }\right))}}\end{array}$$

The coefficients g₀ and g₁ are empirically derived parameters related to the directional nature of the upwelling light field (Q) and f (well known as the f∕Q term). These coefficients depend on sun angle, viewing geometry, wind speed, and the bio-optical state of natural waters and vary with phase function (Morel and Gentili, 1991, 1993, 1996). Thus, appropriate coefficients are needed for different aquatic environments (Lee et al., 2002). For example, the values of g₀ and g₁ were previously suggested as 0.0949 and 0.0794, respectively, for oceanic waters (Gordon et al., 1988). Later, better approximations of g₀ (0.0895) and g₁ (0.1247) were proposed for reflective coastal waters using radiative transfer models on simulated data (Lee et al., 1999). Average values of g₀ and g₁ from Gordon et al. (1988) and Lee et al. (1999) were also used for both coastal and oceanic waters (Lee et al., 1999). The synthetic data for obtaining these historic values of g₀ and g₁ were generated using case-1 radiative transfer models. As the first modification, we updated these coefficients using Hydrolight^® simulations with a four-component case-2 model (Sect. S1) because a shallow water environment was the main focus of this study. In addition, it has been previously suggested that molecular scattering may primarily contribute to R_rs, especially in the blue and green wavelengths in oceanic waters. However, the phase function effect of water molecules could be much smaller than that of particles in nearshore and estuarine waters. Hence, we have avoided separating Eq. (3) as is usually done to account for the phase function effects of individual backscattering contributors (Lee et al., 2013; Zheng et al., 2014).

Next, QAA-v6 uses a set of empirical models (e.g., Eq. 5) based on above-surface R_rs threshold (0.0015 sr⁻¹) to estimate the total non-water absorption coefficient at a reference wavelength in coastal and oceanic waters (Table 2, Level 1B):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\text{if}}{R}_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathrm{670}\right)<\mathrm{0.0015}{\mathrm{sr}}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }{a}_{\mathrm{tnw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)={\mathrm{10}}^{\left(-\mathrm{1.146}-\mathrm{1.366}\times x-\mathrm{0.469}\times x^{\mathrm{2}}\right)},\\ {\displaystyle }& {\displaystyle }\mathrm{\text{where}}x={\log }_{\mathrm{10}}\left({\displaystyle \frac{R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{443}\right)+R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{490}\right)}{R_{\mathrm{rs}}^{\mathrm{0}-}\left({\mathit{\lambda }}_{\mathrm{0}}\right)+\mathrm{5}\times R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{670}\right)\times \frac{R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{670}\right)}{R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{490}\right)}}}\right),\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{\text{otherwise,}}{a}_{\mathrm{tnw}}\left({\mathit{\lambda }}_{\mathrm{1}}\right)=\mathrm{0.39}\times {\left[{\displaystyle \frac{R_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathrm{670}\right)}{R_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathrm{443}\right)+R_{\mathrm{rs}}^{\mathrm{0}+}\left(\mathrm{490}\right)}}\right]}^{\mathrm{1.14}},\end{array}$$

where a_tnw is the total non-water absorption coefficient (m⁻¹), λ₀=555, and λ₁=670 nm.

We avoided using blue wavelengths in our empirical models as blue bands, especially 443 nm or lower, suffer from large errors in atmospheric correction due to the high abundance of CDOM, suspended particles, and absorbing aerosols in a coastal environment. In contrast, the green to red band ratio (GRBR) can be used for estimating the absorption coefficient of an individual water constituent with the primary condition of their dominance in the study region. The GRBR has been used, for example, to monitor water constituents in various estuarine and coastal waters, e.g., a_g355 in Barataria Bay, USA (Joshi and D'Sa, 2015), a_g412 in Apalachicola Bay, USA (Joshi et al., 2017a), a_dg412 in Galveston Bay, USA (D'Sa et al., 2018), suspended particulate matter (D'Sa et al., 2007), and chlorophyll index (harmful algal bloom) in the northern Gulf of Mexico (Qi et al., 2015). Estuarine waters are generally characterized by a high abundance of CDOM, mineral particles, or both and are thus known to have strong light absorption towards shorter wavelengths, sometimes even in the green region. In contrast, the light absorption in the red region usually remains minimal for CDOM and mineral particles. As a result, the green band can be considered as a pilot band to capture variations in dissolved or mineral particle absorption, whereas the red band is a reference band. Hence, small variations in the GRBR are suitable to capture large variations in absorption at shorter wavelengths due to the exponential nature of CDOM and particle absorption. A similar band ratio (RGCI; red to green chlorophyll index) has been used in a semi-analytical approach to obtain chlorophyll-a concentrations for estimating phytoplankton absorption a_ϕ670 and subsequently a_tnw670 in the productive waters of Tampa Bay, USA (Le et al., 2013). The red to green band ratio (RGBR) works well to quantify variations in chlorophyll concentrations, especially in phytoplankton-dominated waters. Therefore, the GRBR (or RGBR) can be overall associated with the dominant water constituent in the study area. However, optically active water constituents collectively contribute to total light absorption, and hence these band ratios can also be used for the remote estimation of a_tnw at green wavelengths, especially in estuarine waters where variations in total absorption coefficients are often noticeable due to a high abundance of one or more (e.g., CDOM and particles) water components. A ratio of R_rs at the red and green wavelengths (e.g., 640∕645 and 555 nm for MODIS) was previously used in the standard QAA to improve the estimates of a_tnw at a reference wavelength (e.g., 555 nm for MODIS) in turbid waters, which worked reasonably in coastal waters (Lee et al., 2002; Chen and Zhang, 2015). The processing of the high-spatial-resolution I-1 band (640 nm; spatial resolution: 375 m) is not supported by NASA's SeaDAS tool; therefore, our analysis was limited to the available ocean color bands (M1 to M5; spatial resolution: 750 m).

Additionally, the selection of the reference wavelength in Eq. (5) is another important factor affecting the retrieval of IOPs in QAA. The use of a reference wavelength at red wavelengths (e.g., > 600 nm) was suggested for relatively turbid coastal waters where particulate and dissolved absorption is much lower than water absorption (Lee et al., 2002; Aurin and Dierssen, 2012). We used 555 nm as a reference wavelength for three reasons: (1) there is a lack of absorption measurements at red wavelengths (e.g., IES dataset) to tune empirical models (Table 2, Level 1B); (2) the use of a red reference wavelength would likely deteriorate the estimations in blue wavelengths due to errors in spectral extrapolation corresponding to the empirical nature of the backscattering power-law exponent, η; and (3) a strong relationship between the green to red band ratio and total non-water absorption at 555 nm is observed in this study.

Table 3 The calibration coefficients for sensor-specific QAA tuning. λ₀ is a sensor-specific reference wavelength.

Once u(λ₀) is obtained as the positive root of Eq. (3) using $R_{\mathrm{rs}}^{\mathrm{0}-}$ and the coefficients g₀ and g₁, the backscattering coefficient at a reference wavelength (b_btnw(λ₀)) can be easily obtained with u(λ₀), a_tnw(λ₀), and the following analytical model (Table 2, Level 1C0):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{b}_{\mathrm{btnw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)\left({\mathrm{m}}^{-\mathrm{1}}\right)=\left({\displaystyle \frac{u\left({\mathit{\lambda }}_{\mathrm{0}}\right)}{\mathrm{1}-u\left({\mathit{\lambda }}_{\mathrm{0}}\right)}}\right)\times (a_{\mathrm{tnw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)+a_{\mathrm{w}}\left({\mathit{\lambda }}_{\mathrm{0}}\right))\\ \text{(6)}& {\displaystyle }& {\displaystyle }-b_{\mathrm{bw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right),\end{array}$$

where a_w and b_bw are water absorption and backscattering coefficients, respectively. The spectral distribution of the particulate backscattering coefficient (b_btnw) can be modeled using the power-law model (Lee et al., 2002; Table 2, Level 1C2):

$$\begin{array}{}\text{(7)}& {\displaystyle }{b}_{\mathrm{bt}}\left(\mathit{\lambda }\right)\left({\mathrm{m}}^{-\mathrm{1}}\right)=b_{\mathrm{bw}}\left(\mathit{\lambda }\right)+b_{\mathrm{btnw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)\times {\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathit{\lambda }}}\right)}^{\mathit{\eta }},\end{array}$$

where η is the spectral shape of b_btnw distribution. The standard QAA-v6 uses the following empirical model to obtain η.

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathit{\eta }=\mathrm{2}\times \left(\mathrm{1}-\mathrm{1.2}\times e^{\left(-\mathrm{0.9}\times \frac{R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{443}\right)}{R_{\mathrm{rs}}^{\mathrm{0}-}\left(\mathrm{555}\right)}\right)}\right)\end{array}$$

The variation in the power-law exponent η depends on water properties and the size of particles according to Mie theory and it is extremely difficult to retrieve η from R_rs in nearshore waters (Aurin and Dierssen, 2012). Thus, we adopted a different approach of obtaining η from b_btnw555 with a linear relationship that was formulated using field observations in the turbid waters near the Mississippi River delta (D'Sa et al., 2007; Table 2, Level 1C1).

$$\begin{array}{}\text{(9)}& \mathit{\eta }=-\mathrm{0.566}-\mathrm{1.395}\times {\log }_{\mathrm{10}}\left(b_{\mathrm{btnw}}\mathrm{555}\right)\end{array}$$

The spectral distribution of total absorption coefficients was then obtained using b_btnw and u (Table 2, Level 1C3).

$$\begin{array}{}\text{(10)}& a_{\mathrm{t}}\left(\mathit{\lambda }\right)\left({\mathrm{m}}^{-\mathrm{1}}\right)=\left({\displaystyle \frac{\mathrm{1}-u\left(\mathit{\lambda }\right)}{u\left(\mathit{\lambda }\right)}}\right)\times b_{\mathrm{bt}}\left(\mathit{\lambda }\right)\end{array}$$

To extend and evaluate the applicability of estuarine-specific QAA tuning, it was further applied to various ocean color (Sentinel-3 OLCI, MODIS Aqua, MERIS, and SeaWiFS) and land-observing sensors (Landsat 8 OLI and Sentinel-2 MSI). The calibration coefficients for obtaining the total non-water absorption coefficient at a reference wavelength (a_tnw(λ₀); Level 1B in Table 2) are given in Table 3.

2.4 Atmospheric correction of satellite imagery

Level-1 satellite imagery was corrected for the atmosphere using the SeaDAS 7.4 image processing tool as described previously (Joshi et al., 2017a), and mean values of a 3×3 pixel box centered at a station location were considered as reasonable satellite matchups for the field measurements. In addition, in situ R_rs values were matched to the central wavelengths of spectral bands using the spectral response functions of respective satellite sensors prior to the sensor-specific tuning and validation:

$$\begin{array}{}\text{(11)}& {\displaystyle }{R}_{{\mathrm{rs}}_{\mathrm{RSR}}}\left({\mathrm{sr}}^{-\mathrm{1}}\right)={\displaystyle \frac{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}\mathrm{RSR}\left(\mathit{\lambda }\right)\times R_{{\mathrm{rs}}_{\mathrm{insitu}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}}{{\int }_{{\mathit{\lambda }}_{\mathrm{1}}}^{{\mathit{\lambda }}_{\mathrm{2}}}\mathrm{RSR}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\end{array}$$

where RSR is the relative spectral response for the satellite sensor, λ₁ is the lower bound of a spectral band, and λ₂ is the upper bound of a spectral band.

Table 4 The performance evaluation of the atmospheric correction procedure (MUMM) in Galveston Bay (USA) during two field campaigns (29 September and 29–30 October 2017). The MRE remained relatively lower than the blue for the green and red bands (bold fonts).

An iterative NIR atmospheric correction scheme was previously evaluated for estuarine environments (Bailey et al., 2010; Werdell et al., 2010; Joshi et al., 2017a); however, it yielded negative R_rs at blue wavelengths and atmospheric correction failure at several pixels in Galveston Bay. In comparison, the errors in atmospheric correction were considerably reduced with the MUMM NIR correction as it was designed for low to moderately turbid waters (Ruddick et al., 2006; Novoa et al., 2017). Furthermore, the validation of atmospheric-corrected VIIRS imagery showed reasonable performance of the MUMM atmospheric correction scheme during both field campaigns in Galveston Bay (Table 4). The MRE, which was relatively higher at the blue wavelengths, was greatly reduced towards the green and red wavelengths. Thus, the success of the atmospheric correction procedure was decided based on the green and red bands, since only these bands were used to tune QAA in this study. However, the observed large errors at other wavelengths (e.g., blue bands) were likely due to the high abundance of CDOM and particles in the study region and the aerosol model selection in the atmospheric correction procedure (Minu et al., 2014). The time difference between field and satellite measurements resulted in an error enhancement at longer wavelengths with relatively smaller errors in October (difference of 0 days) and larger errors in September (difference of +1 days; Table 4). Overall, low MRE in the green and red wavelengths indicated the usefulness of the MUMM atmospheric correction for investigating bio-optical properties with the QAA-V processing chain and VIIRS ocean color data in Galveston Bay (USA).

2.5 Statistical analysis

The algorithm's performance and atmospheric-corrected VIIRS imagery were evaluated using coefficients of determination (R²), root mean square error (RMSE_log 10), bias (Bias_log 10), and absolute mean relative error (MRE).

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathrm{Bias}}_{\log \mathrm{10}}={\displaystyle \frac{\mathrm{1}}{n}}\times \sum _{i=\mathrm{1}}^n\left[{\log }_{\mathrm{10}}\left(y_i\right)-{\log }_{\mathrm{10}}\left(x_i\right)\right]\text{(13)}& {\displaystyle }& {\displaystyle }{\mathrm{RMSE}}_{\log \mathrm{10}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\times \sum _{i=\mathrm{1}}^n{\left[{\log }_{\mathrm{10}}\left(y_i\right)-{\log }_{\mathrm{10}}\left(x_i\right)\right]}^{\mathrm{2}}}\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{MRE}\left(\mathit{\%}\right)={\displaystyle \frac{\mathrm{100}}{n}}\times \sum _{i=\mathrm{1}}^n\left[\left|{\displaystyle \frac{y_i-x_i}{x_i}}\right|\right]\end{array}$$

R software was used to generate synthetic IOPs for the Hydrolight^® case-2 model and the statistical analysis presented in this study.

Figure 3 (a) u (= $b_{\mathrm{b}}/(a+b_{\mathrm{b}})$) vs. $R_{\mathrm{rs}}^{\mathrm{0}-}$ for various studies (this study: orange; Lee et al., 1999: blue; Lee et al., 2002: red; Gordon et al., 1988: black) using HL datasets. (b) The relationship between green to red band ratio and a_tnw555 for different datasets (HL: blue, IES: green, NOMAD: grey). Black line shows a threshold to facilitate a smooth transition from in situ to synthetic data in modeling a_tnw555. A dashed line separates data with negative ρ.

Figure 4 Validation of $R_{\mathrm{rs}}^{\mathrm{0}-}$ vs. u model using HL data for (a) a_t555 and (b) b_bt555. Validation of QAA-V-modeled a_tnw555 with “true” a_tnw555 for (c) the HL synthetic dataset (bias${}_{\log \mathrm{10}}=-\mathrm{0.0208}$, RMSE_log10=0.0963 m⁻¹) and (d) IES dataset (bias${}_{\log \mathrm{10}}=-\mathrm{0.0294}$, RMSE_log10=0.190 m⁻¹). Validation of QAA-v6-modeled a_tnw555 with “true” a_tnw555 for (e) the HL synthetic dataset (bias${}_{\log \mathrm{10}}=-\mathrm{0.1180}$, RMSE_log10=0.1490 m⁻¹) and (f) IES dataset (bias${}_{\log \mathrm{10}}=-\mathrm{0.1252}$, RMSE_log10=0.278 m⁻¹). Note: the training set includes EcoHAB and Tampa Bay data (Table 1). IOP measurements of Galveston Bay (this study) are shown with orange hexagons. MRE: absolute mean relative error (%). TB: Tampa Bay, AB: Apalachicola Bay, BB: Barataria Bay, CB: Chesapeake Bay, FB: Florida Bay, GoM: Gulf of Mexico, HI: Horn Island, LE: Lake Erie, MS: Mississippi Sound, MB: Mobile Bay, SWFL: southwest Florida, GB: Galveston Bay (Table 1).

3 Results

3.1 Modifications to standard QAA

A new set of g₀ (0.0788) and g₁ (0.2379) obtained using Hydrolight-simulated IOPs and R_rs is proposed for the highly attenuating waters of this study. These coefficients showed a notable departure from the quadratic relationship, especially in highly scattering waters when compared to Gordon et al. (1988) and Lee et al. (1999, 2002; Fig. 3a). We used the historical values g₀ (0.0895) and g₁ (0.1245) for less reflective nearshore and shelf waters (i.e., threshold ρ ≥ 0.25) in the validation analysis (Fig. 3b). Overall, the threshold-based selection of the u vs. $R_{\mathrm{rs}}^{\mathrm{0}-}$ model showed a valid retrieval of u ($\approx b_{\mathrm{bt}}/(a_{\mathrm{t}}+b_{\mathrm{bt}})$) as seen when u was analytically used to obtain the total absorption coefficient (a_t555) and the total backscattering coefficients (b_bt555; HL dataset) and vice versa (Fig. 4a, b; Table 2, Level 1A).

The performance of QAA-V was largely dependent on the estimation of the total non-water absorption coefficient (a_tnw) in Level 1B (Table 2). In this study, the empirical model of standard QAA-v6 (Eq. 5) was replaced by a tuned empirical power-law relationship using a training set (IES; N=121; Table 1) of in situ observations and Hydrolight^® synthetic data (HL; N=561):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_{\mathrm{tnw}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)=\left\{\begin{array}{l}{\mathrm{10}}^{(\mathrm{0.139}-\mathrm{1.788}\times \mathit{\rho }+\mathrm{0.490}\times {\mathit{\rho }}^{\mathrm{2}})}\mathrm{if}\mathit{\rho }<\mathrm{0.25}\\ {\mathrm{10}}^{(\mathrm{0.406}-\mathrm{2.940}\times \mathit{\rho }+\mathrm{0.928}\times {\mathit{\rho }}^{\mathrm{2}})}\mathrm{if}\mathit{\rho }\ge \mathrm{0.25}\end{array}\right\},\\ \text{(15)}& {\displaystyle }& {\displaystyle }\mathit{\rho }={\log }_{\mathrm{10}}\left({\displaystyle \frac{R_{\mathrm{rs}}^{\mathrm{0}-}\left({\mathit{\lambda }}_{\mathrm{0}}\right)}{R_{\mathrm{rs}}^{\mathrm{0}-}\left({\mathit{\lambda }}_{\mathrm{1}}\right)}}\right),\end{array}$$

where a_tnw is the total non-water absorption coefficient, λ₀=551 or 555 nm, and λ₁=671 nm.

A threshold value of 0.25 was set for ρ to merge the HL and IES datasets excluding the NOMAD (Fig. 3b). The a_tnw555 (ρ>1.0) nearly reached the lower limit close to zero for oceanic waters, likely due to low concentrations of reflecting and absorbing materials. Thus, we suggest the upper threshold of ρ=0.65 beyond which Level 1B fails and overestimates water a_tnw(λ₀). Likewise, negative ρ values of synthetic data represented CDOM-rich waters with very strong absorption even at green wavelengths (e.g., 555 nm). For HL data, the modeled a_tnw555 showed a reasonable performance of the green to red band ratio model (MRE = 16.3 %, bias${}_{\log \mathrm{10}}=-\mathrm{0.0208}$, RMSE_log 10=0.0963 m⁻¹, N=561; Fig. 4c). The performance of the empirical model, however, showed a significant difference within the IES training data (R²=0.90, MRE = 21 %, bias${}_{\log \mathrm{10}}=-\mathrm{0.0015}$, RMSE_log 10=0.12 m⁻¹, N=120) and testing data (R²=0.72, MRE = 34 %, bias${}_{\log \mathrm{10}}=-\mathrm{0.0294}$, RMSE_log 10=0.19 m⁻¹, N=209; Fig. 4d).

Figure 5 Statistical assessment of QAA-V in comparison to QAA-v6 (updated on May 2015) using synthetic data (HL; a) and estuarine and nearshore field data (IES; b). Color scheme indicates a_tnw and b_btnw at 411 nm (dark blue), 443 nm (blue), 489 nm (cyan), 555 nm (green), and 532 nm (dark green) wavelengths. MRE: absolute mean relative error (%). RMSE and bias are in m⁻¹.

3.2 Comparison of QAA-V with the standard QAA-v6

The QAA-v6 algorithm was applied separately on HL synthetic data and IES field observations for a direct comparison with QAA-V in retrieving optical properties in estuarine waters. Figure 4e shows the performance of QAA-v6 in estimating the total non-water absorption coefficient at 555 nm (a_tnw555) with HL synthetic data. In comparison to QAA-V, the QAA-v6 showed an obvious underestimation at 555 nm, especially for a_tnw555 > ∼ 0.3 m⁻¹. Statistical assessment showed that QAA-V is more accurate than the QAA-v6 with approximately 83 % less bias, 35 % decreased RMSE, and 1.5-fold lower MRE (Fig. 5, Sect. S2). For IES data, the QAA-v6 showed a clear difference with a poor performance at several stations (Fig. 4f) compared to QAA-V (Fig. 4d). Furthermore, the retrieval errors were large towards the upper and lower ends (e.g., a_tnw555 > ∼ 0.3 and < ∼ 0.1 m⁻¹). Overall, the QAA-v6 had an obvious underestimation with approximately 75 % greater bias, 31 % increased RMSE, and 2-fold higher MRE than the QAA-V at 555 nm (Fig. 5, Sect. S2).

Figure 6 Validation of QAA-V on synthetic data (HL) for (a) a_tnw and (b) b_btnw at 411 nm (dark blue), 443 nm (blue), 489 nm (cyan), and 555 nm (green) wavelengths. MRE: absolute mean relative error (%).

Figure 7 Evaluation of QAA-V's performance on the estuarine and nearshore dataset (IES) for (a) a_tnw443 and (b) b_btnw532 and on the NOMADs for (c) a_tnw at 411 nm (dark blue), 443 nm (blue), and 488 nm (cyan) wavelengths. (d) SPM–b_btnw532 relationship (R²=0.89 without outliers) that was formulated based on field observations in Apalachicola Bay, Barataria Bay, and Galveston Bay. (e) A comparison of the SPM–b_btnw532 relationship in this study (grey triangles) with SPM models from Nechad et al. (2010; blue circles) and D'Sa et al. (2007; orange squares) and with field observations in AB, BB, and GB (black stars). TB: Tampa Bay, AB: Apalachicola Bay, BB: Barataria Bay, CB: Chesapeake Bay, FB: Florida Bay, GoM: Gulf of Mexico, HI: Horn Island, LE: Lake Erie, MS: Mississippi Sound, MB: Mobile Bay, SWFL: southwest Florida, GB: Galveston Bay (Table 1). Training set includes EcoHAB and Tampa Bay datasets. MRE: absolute mean relative error (%).

Figure 8 (a) Validation of QAA-V-derived SPM vs. in situ SPM for (b) 30 September 2017 (MRE = 39.9 %, N=13), (c) 29 October 2017 (MRE = 39.1 %, N=14), and (d) 30 October 2017 (MRE = 26.6 %, N=7). Stations for validation analysis (Fig. 8a) are also illustrated in corresponding SPM maps.

3.3 Evaluation of QAA-V on synthetic HL data, NOMAD, and IES data

Hydrolight-simulated case-2 water $R_{\mathrm{rs}}^{\mathrm{0}+}$ spectra were fed into QAA-V to derive a_tnw and b_btnw at 411, 443, 489, and 555 nm (Fig. 6). Although a negative bias indicated an overall underestimation of modeled absorption coefficient at 555 nm (Fig. 4c), QAA-V performed satisfactorily in the blue region (Fig. 6a). Furthermore, the modeled a_tnw at blue wavelengths showed relatively lower MRE < 13 % and RMSE_log 10 < 0.075 m⁻¹ despite the errors in modeled a_tnw555 at values less than 0.3 m⁻¹ (relatively low CDOM and mineral particle abundance), indicating the secondary importance of the green reference wavelengths in transferring errors to blue wavelengths in QAA-V's processing chain compared to the red reference wavelengths (Aurin and Dierssen, 2012). Similarly, QAA-V also estimated b_btnw with MRE < 16 % at four wavelengths (Fig. 6b). When compared to the QAA-v6, QAA-V performed with approximately 80–90 % lower bias, 30–40 % reduced RMSE, and 1–2-fold decreased MRE at blue wavelengths (Fig. 5, Sect. S2). QAA-V's performance on the standard NOMAD showed that the MRE was < 30 % for all blue wavelengths with a_tnw443 being the best-retrieved parameter (R²=0.94, MRE = 24.3 %, N=547). However, large errors were observed for a_tnw < 0.1 m⁻¹ at all wavelengths (Fig. 7a).

For field validation, $R_{\mathrm{rs}}^{\mathrm{0}+}$ spectra from 15 experiments of IES data were fed into QAA-V for obtaining a_tnw443 and b_btnw532 (Fig. 7b and c; Table 1). Modeled IOPs showed overall good results for a_tnw443 (MRE = 15 % for training data and 23 % for testing data; Fig. 7b) and b_btnw532 (MRE = 32 % for all data and 26 % for all data except Apalachicola Bay; Fig. 7c). A further analysis of individual datasets showed noticeable variations in MRE. For example, MRE varied from ∼ 11–24 % among the five best-performing datasets for estimating a_tnw443, whereas it varied from ∼ 40–75 % for three poorly performing datasets with preliminary processing status (Table 1). On the IES data, QAA-V showed fewer errors in retrieving water IOPs than QAA-v6 (Fig. 5, Sect. S2).

3.4 SPM–b_btnw532 relationship and validation

A linear relationship (R²=0.89; N=33) was observed between b_btnw532 and SPM (Eq. 16; Fig. 7d). It was used for evaluating the applicability of QAA-V to generate synoptic maps of SPM using the VIIRS ocean color sensor in estuarine environments.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{SPM}=\mathrm{103.07}\times \left(b_{\mathrm{btnw}}\mathrm{532}\right)+\mathrm{0.24};\text{where}{b}_{\mathrm{btnw}}\mathrm{532}\\ \text{(16)}& {\displaystyle }& {\displaystyle }=\text{backscattering coefficient at 532\hspace{0.17em}nm}\end{array}$$

It is important to note that the backscattering sensor had an upper threshold of 0.23 m⁻¹ at which the sensor saturated and failed to detect variability in SPM; however, the linearity in the SPM–b_btnw532 relationship was assumed beyond the threshold (0.23 m⁻¹) in this study. The SPM–b_btnw532 relationship was compared to a similar relationship developed in turbid coastal waters of the Mississippi River (D'Sa et al., 2007) and a generic multisensor algorithm (Nechad et al., 2010) using the IES dataset (Fig. 7e; Table 2). The three algorithms showed quite similar trends despite different water properties of the IES dataset. When compared to field observations of SPM in Apalachicola Bay, Barataria Bay, and Galveston Bay, these algorithms showed good performance with the lowest error for the single-wavelength generic algorithm of Nechad et al. (2010; MRE = 32.1 %, N=57) and the highest error for the SPM–b_btnw555 relationship of D'Sa et al. (2007; MRE = 38.8 %, N=57; Fig. 7e).

Figure 9 (a) River discharge of the Trinity River at Romayor, Texas (black line) and Wallisville, Texas (blue line); the west flank of the San Jacinto River (red line) and the east flank of the San Jacinto River (green line). The river discharge at the Wallisville site was not corrected for the tides. Black bars represent the days of VIIRS satellite imagery for SPM analysis in Galveston Bay. (b–g) Post-hurricane SPM maps of Galveston Bay. TB: Trinity Bay, UGB: upper Galveston Bay, LGB: lower Galveston Bay, EB: East Bay (Fig. 1c).

3.5 Evaluation of QAA-V in a turbid estuarine environment

The SPM–b_btnw532 relationship in this study was applied to VIIRS imagery for evaluating satellite-based application of QAA-V in a turbid estuarine environment (Fig. 8). The regional SPM–b_btnw532 relationship showed ∼ 35 % overall MRE in a satellite–field comparison (Fig. 8a) for Galveston Bay during the two field campaigns (Fig. 8b–d). Both field and estimated SPM concentrations followed a similar pattern of high to low values along the north-to-south transect. A large-to-small error trend from 30 September (MRE = 39.9 %) to 30 October (MRE = 26.6 %) was similar to the error trend in the atmospheric correction of the VIIRS imagery (Table 4). Based on these results, a sequence of cloud-free and atmospherically corrected VIIRS imagery was converted to SPM maps using the QAA-V processing chain (Table 2) for analyzing post-hurricane SPM dynamics in Galveston Bay (Fig. 9).

Figure 10 Application of sensor-specific QAA-V tuning to obtain maps of a_tnw443 using (a) VIIRS, (b) MODIS Aqua, and (c) Sentinel-3 OLCI on 29 October 2017 and (d) VIIRS and (e) Landsat 8 OLI on 30 September and 29 September 2017, respectively. The validation of these maps with the field observations along the transect (St. 1 to St. 14) is shown in (I) for Fig. 8a–c and in (II) for Fig. 8d and e. The maps of b_btnw470 were obtained similarly for (f) VIIRS (29 October 2017), (g) MODIS Aqua, (h) Sentinel-3 OLCI, (i) VIIRS (30 September 2017), and (j) Landsat 8 OLI (29 September 2017) with their validation results in (III) and (IV). Parameter values beyond the upper limit of tuned QAA (ρ > 0.65) are shown masked in white.

3.6 Extending the QAA-V tuning to various satellite sensors

The estuarine-specific green to red band tuning was further applied to evaluate and extend its applicability to past and present ocean color (e.g., SeaWiFS, MERIS, MODIS Aqua, and Sentinel-3 OLCI) and land-observing sensors (Landsat 8 OLI and Sentinel-2 MSI; Table 3). The validation analysis showed a promising performance of QAA-V tuning in obtaining the total non-water absorption coefficient (a_tnw443) and total-non water backscattering coefficient (b_btnw470) in the optically complex and shallow waters of Galveston Bay (Fig. 10). Overall, different satellite sensors showed similar trends of a_tnw443 and b_btnw470 along the transect despite having different spectral and spatial sensor resolutions (Fig. 10I–IV). The MREs were ∼ 15, 9, and 12 % for a_tnw443 retrievals from VIIRS, MODIS A, and Sentinel-3 OLCI sensors, respectively (Fig. 10a–c and I), whereas they were ∼ 26, 7, and 22 % for b_btnw470 retrievals on 29 October 2017 (Fig. 10f–h and III). For Landsat 8 OLI, the MREs were ∼ 20 and ∼ 10 % for a_tnw443 and b_btnw470, respectively, on 29 September 2017 (Fig. 10e, j, II, IV).

4 Discussion

4.1 QAA-V

The standard QAA (Lee et al., 2002) was first tuned with VIIRS sensor bands as QAA-V and then extended to other satellite sensors to obtain more accurate estimates of water IOPs (e.g., total absorption and backscattering coefficients) in shallow estuarine and nearshore waters. For this purpose, two modifications were applied to the QAA: (1) the coefficients g₀ and g₁ of a quadratic model were updated using Hydrolight^® simulations for highly absorbing and highly scattering waters (Eq. 3, Fig. 3a); and (2) an empirical model of QAA for obtaining the total absorption coefficient at a reference wavelength was replaced by a set of empirical models that were optimized for highly attenuating estuarine and nearshore waters (Eq. 15, Fig. 3b). The validity of these modifications is demonstrated in Fig. 4. The updated coefficients, g₀ and g₁, showed a significant departure from historical values (Gordon et al., 1998; Lee et al., 1999, 2002) especially at green and red wavelengths in highly attenuating waters. Likewise, these coefficients were updated using a synthetic dataset representing highly attenuating waters and they may not perform satisfactorily in less reflective nearshore and coastal environments. This limitation was addressed by adopting a water-type switching of the u vs. $R_{\mathrm{rs}}^{\mathrm{0}-}$ algorithm based on a green to red band ratio threshold (ρ) in the QAA-V processing chain (Table 2). The effectiveness of tuning g₀ and g₁ is also supported by a previous study showing an approximate 5-fold bias reduction in obtaining total non-water IOPs at 440 nm (Aurin and Dierssen, 2012).

Overall, the negative bias for the HL datasets indicated that the threshold-based empirical models underestimated a_tnw555; however, a major error in the model performance was observed at the lower end of the absorption of values, likely for two reasons: (1) the failure of the threshold (ρ) in providing a smooth switching of empirical relationships between the highly absorbing and scattering waters (synthetic HL dataset) and estuarine and nearshore waters (IES dataset; Fig. 4b); and (2) the unsuitability of the empirical relationship for some R_rs spectra in the HL synthetic data based on the green to red band ratio (e.g., a_tnw555 < 0.3 m⁻¹; Fig. 4a). In contrast, both empirical relationships worked reasonably well in providing a smooth transition from nearshore waters to the highly turbid and highly absorbing waters of the IES data (Fig. 4b). However, the large differences within an IES dataset (e.g., training data and test data) could be explained by the processing status of different individual experiments in the IES dataset (Table 1). Training data included R_rs measurements in the final processing status from two well-known experiments, namely the EcoHAB (N=74) and Tampa Bay monitoring programs (N=47; Table 1). Individual absolute mean relative errors (MREs) for these data are 24 and 16 %, respectively. In contrast, few datasets used in the testing set were in the preliminary processing stage with MRE > 35 % (e.g., Chesapeake Bay Light Tower, Horn Island, Cojet 7, and Lake Erie).

The validation of QAA-V in a variety of waters yielded reasonable performance as shown in Figs. 6 and 7. For the HL dataset, an error difference (∼ 2 %) between a_tnw (Fig. 6a) and b_btnw (Fig. 6b) could be associated with errors in Levels 0 and 1A of the QAA processing chain (Table 2). The effect of this systematic error (∼ 1–2 %) is also present in the estimates of backscattering coefficients at blue wavelengths; this could be due to uncertainties in various parameters of the semi-analytical power-law model such as the power-law exponent η (Table 2; Level 1C1). It has previously been shown that η is important for obtaining b_bt and its tuning can be responsible for an approximate 4-fold decrease in the percentage difference for b_bt440 in estuarine waters over QAA's η model (Eq. 8; Aurin and Dierssen, 2012). As we used a linear model based on b_btnw(λ₀) for η (D'Sa et al., 2007), the modeled backscattering coefficient may suffer from errors due to uncertainties in the empirical estimation of the power-law exponent. η does not play any role in retrieving b_btnw(λ₀), and hence errors associated with η should not affect modeled b_btnw(λ₀). However, the errors in η may affect the retrieval of b_btnw at other wavelengths (e.g., b_btnw470; Table 2). For the NOMAD, the observed large MRE in relatively clearer waters (a_tnw < 0.1 m⁻¹) can be due to the fact that the empirical relationships for modeling a_tnw(λ₀) were not designed to work with very low values of non-water absorption coefficients, as generally observed in shelf and oceanic environments (Fig. 7a; Eq. 5). For the IES datasets, the retrieval of a_tnw443 is obviously improved compared to a_tnw555 (Figs. 4d and 7b); however, a 1- to 2-fold increase in errors is observed in IOP retrieval when QAA-V is applied to the field data in contrast to synthetic data. This error enhancement can be due to several factors, including (1) failure to achieve a smooth switching between highly absorbing and reflecting waters and estuarine and nearshore waters in real data (Fig. 4a), (2) inclusion of preliminary data in validation analysis (Table 1), (3) the uncertainty propagation of the previous levels of the QAA-V processing chain (e.g., errors in various empirical, semi-analytical, and analytical relationships; Table 2), and (4) the difference in methodology of data collection, raw data processing, instrumentation, and measurement and instrument errors. For example, reflectance measurements, the primary input of QAA-V, can have large errors if measured in turbid waters and with large viewing angles (Lee et al., 1999), while bottom contamination may introduce an additional error since in situ data represent shallow water environments. The bottom effect could probably be reduced with sophisticated correction schemes (Lee et al., 1998, 1999; Li et al., 2017).

A statistical comparison showed that the QAA-V provided better estimations of IOPs than the standard QAA-v6 in estuarine and nearshore waters (Fig. 5). Also, a green to red band ratio is better suited for the satellite application of QAA in shallow waters due to fewer errors in these bands (Table 4). In contrast, several studies have demonstrated an overestimation in blue bands in turbid waters, likely due to errors in the atmospheric correction (Joshi et al., 2017a; Chen and Zhang, 2015, and reference therein). This overestimation could lead to the observed underestimation of a_tnw555 and subsequently the underestimation of a and b_b at different wavelengths in lower levels of the QAA processing chain (Fig. 5, Table 2).

4.2 The SPM optical model

Traditionally, suspended particulate matter (SPM) has been retrieved from remote sensing imagery using single or multiband empirical relationships between above-surface R_rs and SPM concentrations (Doxaran et al., 2002; Miller and McKee, 2004; D'Sa et al., 2007; Han et al., 2016). While empirical relationships are easy to implement, they are regionally limited and may cause large errors if applied to different waters due to differences in particle properties such as absorption, particle size and composition, and refractive index. Furthermore, above-surface R_rs, even at the red and NIR wavelengths at which particle attenuation controls the R_rs signal, is not a “true” representation of particles as the R_rs signal is also marginally contributed to by the absorption of water and other optically active materials. Hence, this approach has a limitation as it can have robust performance only in particle-dominated waters. Another way to estimate SPM concentration is to relate it to particle backscattering coefficients (b_bp). Several studies have reported the efficacy of this approach in estimating concentrations of total suspended material (TSM) and particulate organic carbon (POC) in a variety of waters ranging from estuarine to open oceans (Stramski et al., 1999; Loisel et al., 2001; Aurin and Dierssen, 2012).

When the SPM–b_btnw532 relationship was applied to the VIIRS imagery, both field- and satellite-estimated SPM concentrations showed a similar pattern of high to low along a north-to-south transect in Galveston Bay; however, the differences remained elevated in the turbid region of Trinity Bay (St. 2 to St. 6 in Fig. 8b and St. 16 to St. 23 in Fig. 8c). Factors contributing to these differences include (1) the error propagation from various steps of the QAA-V processing chain to b_btnw532 (e.g., 20–30 %, Fig. 7c) and hence further down to the SPM inversion, (2) the uncertainty in the atmospheric-corrected green and red R_rs (e.g., 5–20 %, Table 4), (3) the uncertainty in the SPM–b_btnw532 relationship due to limited observations, (4) the assumption of linearity in the SPM–b_btnw532 model beyond the instrument threshold, which may not hold well because the b_btnw532 to SPM ratio depends on the particle nature and it may not always be constant, especially in highly turbid waters, and (5) errors in SPM measurements. For example, SPM is usually measured with a 0.7 µm (average pore size) GF/F filter, which only represents the total concentration of particles greater than this size. However, while smaller particles may not have significant contributions to the mass-specific property of SPM (e.g., concentration), particles smaller than this size can contribute notably to the underestimation of SPM concentrations in highly turbid waters (Sect. S3). Furthermore, small particles and even colloidal particles (< 0.2 µm) are known to contribute significantly to total particle backscattering in coastal waters (Zhang and Gray, 2015; Zhang et al., 2011). For instance, high winds associated with the passage of a cold front on 28 October 2017 (not shown) could have resulted in the resuspension of smaller particles at the shallower stations on 29 October (Fig. 8c) and contributed to greater differences between satellite estimates and field SPM measurements. Interestingly, on the following day (30 October) under calmer conditions, the differences reduced substantially (Fig. 8d). Thus, the mismatch between mass-specific and optical properties could be a major source of error in the SPM–b_btnw532 relationship and hence the observed difference in field–satellite matchups.

4.3 Satellite application of QAA-V to Galveston Bay: post-Hurricane Harvey SPM dynamics

A sequence of SPM maps derived from VIIRS imagery using QAA-V revealed interesting patterns of SPM concentrations in Galveston Bay following Hurricane Harvey (Figs. 8 and 9). Hurricane Harvey, a category 4 hurricane, caused catastrophic flooding in the Houston metropolitan area and the surrounding regions of Galveston's Bay drainage basin with great potential to degrade the bay's water quality.

The dramatic increase in the discharge of water from the Trinity and San Jacinto rivers into Galveston Bay (Fig. 9a) following the record rainfall and flooding associated with Hurricane Harvey reveals different patterns of discharge, with the San Jacinto River peaking and retreating to pre-hurricane levels much quicker than the Trinity River. This reflects the differences of the two river basins, as well as the intensity and variability of the precipitation associated with the hurricane in the two basins. A sequence of post-hurricane SPM maps of Galveston Bay (Fig. 9b–g) reveals distinct spatial and temporal patterns of SPM variations within and outside the bay, including the offshore shelf waters that appeared to be strongly influenced by the river discharge and wind forcing. The immediate effect of Hurricane Harvey was clearly evident on 31 August when unusually high SPM concentrations (> 75 mg L⁻¹) were observed throughout the bay corresponding to high freshwater inputs from the Trinity and the San Jacinto rivers into Galveston Bay (Fig. 9a and b). The SPM-rich plume (> 75 mg L⁻¹) extended from the Bolivar Roads Pass to a large region of the coastal and shelf waters. In contrast, the bay experienced elevated SPM concentrations on 2 September, but the plume was limited by wind forcing and a reduction in fresh water to the bay (Fig. 9a and c).

The strong northeasterly winds (∼ 5–6 m s⁻¹) observed on 7 and 8 September appeared to restrict the SPM plume closer to the bay entrance and the inner shelf waters (Fig. 9d and e). Despite the reduced freshwater inflow into the upper Galveston Bay, sustained freshwater inputs from the Trinity River (∼ 500 m³ s⁻¹) and strong northeasterly winds (e.g., sediment resuspension) could have resulted in elevated SPM concentrations throughout the bay. Furthermore, wind-induced downwelling currents appeared to transport low SPM offshore waters nearshore, while the high SPM nearshore and plume waters were likely downwelled and could have eventually settled into the shelf sediments.

Despite the noticeable reduction in the Trinity River flow, SPM concentrations remained high within the bay on 12 September, indicating the importance of wind-induced sediment resuspension for SPM dynamics in the shallow water environments. However, outside the bay, a well-defined and elevated SPM plume extended offshore, likely associated with the southwesterly winds, which induced the offshore transport of the inner shelf waters (Fig. 9f). On 16 September, total freshwater inflow to Galveston Bay was reduced significantly (∼ 500 to ∼ 100 m³ s⁻¹ in the Trinity River) with SPM generally reduced throughout the bay. However, Trinity Bay and East Bay showed relatively higher SPM than the upper and lower Galveston Bay (Fig. 9g), likely due to a more delayed discharge through the Trinity River and the wetlands.

Overall, the QAA-V-based SPM maps of Galveston Bay showed distinct variations in SPM concentrations following Hurricane Harvey. Although the Trinity Bay and the upper Galveston Bay responded similarly a few days following the hurricane, distinct SPM patterns emerged (e.g., lower SPM in the western and higher in the eastern part of the bay) after 2 weeks, suggesting different influences in the eastern and western parts of the bay. For example, floodwaters from the Houston metropolitan and surrounding region appeared to have receded within a few days of the hurricane event (Fig. 9a), whereas the discharge of floodwaters was elevated through the Trinity River over the course of several weeks after Hurricane Harvey. These floodwaters could have accumulated first in wetlands and numerous water bodies in the eastern region and lower Trinity basin during the hurricane event and eventually emptied to the main Trinity River channel. SPM was also elevated in the East Bay, but concentrations were generally lower within the first 2 weeks, suggesting this region of the bay mostly remained isolated from the other regions of Galveston Bay. Nonetheless, this region receives discharge from the surrounding wetlands, which could have been elevated during this period. Overall, wind forcing was also important in controlling the extent and the dispersal of the sediment-rich plume waters into the shelf and in contributing to the SPM variability within the bay due to sediment resuspension and its transport into the shelf waters.

4.4 Application of the QAA tuning to various ocean color and land-observing sensors

Sensor-specific QAA tuning (e.g., QAA-V for VIIRS, MODIS, Landsat 8 OLI, and Sentinel-3 OLCI) showed an overall valid retrieval of absorption and backscattering coefficients with various ocean color and land observation sensors (Fig. 10). Although satellite-derived values and trends of a_tnw443 and b_btnw470 are similar to the field observations, the observed discrepancies could, in addition to the uncertainties in field measurements, be due to several sources of errors. For example, it is well known that satellite products suffer from large errors in the blue region, especially in coastal waters due to atmospheric correction (Table 4; Sect. S4). The large errors between field- and Landsat-8-OLI-derived a_tnw443 could have been due to the fact that the QAA processing chain uses these erroneous blue R_rs values to obtain a_tnw at blue wavelengths (Table 2; Level 3). However, the b_btnw retrievals at blue wavelengths are unaffected by the blue R_rs inputs (Table 2; Level 2). Likewise, the errors were relatively smaller at the reference wavelength (Table 4) because the proposed QAA tuning avoided using blue wavelengths in the primary step of getting a_tnw and b_btnw at a reference wavelength. Hence, the atmospheric correction procedure is an important step that would impact the performance of QAA-V in ocean color applications of shallow estuarine and nearshore waters.

Another important discrepancy among various sensor maps is the number of masked pixels in shelf waters (Fig. 10). We have used a VIIRS band-ratio-based threshold (ρ=0.65) to separate green waters (e.g., productive coastal waters) and blue waters (e.g., open ocean; Eq. 15) and applied it to various satellite sensors. Although this threshold worked well for the sensors with similar green and red bands (e.g., MODIS A and VIIRS), it did not perform as well for Sentinel-3 OLCI and Landsat 8 OLI because of notable differences between the green and red bands of these sensors and the VIIRS sensor (Table 3). Thus, while the proposed threshold works well to represent estuarine and nearshore waters for various ocean color and land-observing sensors, it could be further optimized for each satellite sensor.

5 Conclusions

A multiband quasi-analytical algorithm tuned for the VIIRS ocean color sensor (QAA-V) and for estuarine and nearshore waters was proposed. Two major changes were applied to the standard QAA (Lee et al., 2002): (1) the coefficients g₀ and g₁ of a semi-analytical quadratic relationship were updated to obtain u from the R_rs (Eq. 1), and (2) a threshold-based empirical model was proposed using the green to red band ratio (GRBR) to estimate the total absorption coefficient at a reference wavelength. The QAA-V-derived total absorption and backscattering coefficients showed a good relationship in a variety of waters ranging from highly turbid and highly absorbing (MRE < 17 %) to relatively clearer coastal waters (MRE < 30 %). Moreover, a reasonable performance (MRE < 25 %) using in situ estuarine and nearshore data indicated the usefulness of the GRBR in modeling total absorption coefficients in estuarine waters regardless of the dominance of one or more water constituents (e.g., CDOM, mineral particles, or phytoplankton). This band ratio needs to be explored further for various ocean color sensors and in different estuarine and nearshore environments around the world. The QAA-V may not perform satisfactorily in optically shallow waters as the empirical relationships were designed specifically for optically deep environments. This study showed good retrieval of backscattering coefficients (MRE = ∼ 25–30 %); however, the errors increased towards the lower levels of the QAA-V processing chain (e.g., MRE = ∼ 35 % in estimating SPM), likely due to the lack of b_btnw–SPM matchups in a variety of waters in formulating a robust SPM–b_btnw532 relationship corresponding to estuarine waters. This limitation suggests a great need for in situ backscattering measurements in various shallow environments and their availability for public use.

The QAA-V and a regional SPM–b_btnw532 relationship were applied to a sequence of VIIRS imagery for investigating post-hurricane SPM dynamics in Galveston Bay. Despite noticeable errors, especially in turbid regions of the bay, the application of QAA-V showed great potential in revealing the SPM patterns due to post-hurricane variations linked to freshwater inflow to the bay and wind forcing. A sequence of SPM maps after the passage of Hurricane Harvey showed that Galveston Bay received massive amounts of SPM due to large volumes of freshwater input from the two major rivers and the surrounding regions. However, while the freshwater inflow reduced within a few days in the western part of the bay, it remained high over the course of several weeks in the eastern part of the bay, providing evidence of the short-term storage capacity of wetlands and numerous freshwater reservoirs in the lower Trinity basin. This discharge pattern coupled with different residence times in the western and eastern parts of the bay resulted in distinct SPM patterns in the two regions of the bay. Furthermore, while winds appeared to have played an important role in resuspending sediments within the bay, they were critical in the transport and dispersion of sediments into the shelf waters of the Gulf of Mexico.

This study did not address bottom reflectance effects (e.g., during clear-water conditions), Raman scattering, and chlorophyll fluorescence that may have degraded the QAA-V's performance. Furthermore, the tuning of the semi-analytical and empirical relationships (e.g., η in Table 2) was not possible due to scarcity of field IOPs (e.g., backscattering measurements). Although further refinements and validation studies are necessary to improve the performance and applicability of QAA-V in spatially and temporally distinct shallow waters around the world, the promising results of this study suggest that the application of QAA-V to various ocean color and land observation satellites can be a useful tool to assess the bio-optical state and water quality dynamics in a variety of coastal systems around the world.