Global Estimation of Suspended Particulate Matter From Satellite Ocean Color Imagery.

The suspended particulate matter (SPM) concentration (unit: mg l-1) in surface waters is an essential measure of water quality and clarity. Satellite remote sensing provides a powerful tool to derive the SPM with synoptic and repeat coverage. In this study, we developed a new global SPM algorithm utilizing the remote sensing reflectance (R rs (λ)) at near-infrared (NIR), red, green, and blue bands (NIR-RGB) as input. The evaluations showed that the NIR-RGB algorithm could predict SPM with the median absolute percentage difference of ∼35%-39% over a wide range from ∼0.01 to >2,000 mg l-1. The uncertainty is smaller (29%-37%) for turbid waters where R rs (671) ≥ 0.0012 sr-1 and slightly higher (41%-44%) for clear waters where R rs (671) < 0.0012 mg l-1. The algorithm was implemented with the global R rs (λ) data from the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi National Polar-orbiting Partnership (SNPP) satellite. We provided a brief characterization of the spatial distribution and temporal trends of the SPM products in global waters based on the monthly SPM composites. Case studies of the SPM time series in coastal and inland waters suggest that the satellite SPM estimations registered spatial and seasonal variation and episodic events in regional scales as well. The VIIRS-generated global SPM maps provide a valuable addition to the existing ocean color products for environmental and climate applications.

Introduction

Remote sensing of the suspended particulate matter (SPM) concentration (unit: mg l−1) in the surface water is driven by its close ties with the environment‐ and climate‐related problems. The suspended particles, living or nonliving, organic or inorganic, absorb and scatter sunlight penetrating through the upper water column, influencing the water turbidity and color (Kirk, 1994). Additionally, pollutants can attach to the suspended particles and disperse farther with them (James, 2002). Government‐led environmental agencies, such as the US Environmental Protection Agency (EPA), have designated SPM as a conventional water quality parameter (Bilotta & Brazier, 2008). In the land‐ocean interface zones, sediments from river discharge, shoreline erosion, and bottom resuspension can be transported by currents or settle down the water column. Satellite‐derived SPM products are useful for tracing ocean currents (Yang et al., 2014), for sediment transport modeling with satellite assimilation (Ouillon et al., 2004; Stroud et al., 2009) and for global land‐ocean flux (Milliman & Farnsworth, 2011). In open oceans, the suspended particles within the first optical depth are dominantly biogenic. As an important carbon pool, the particulate organic carbon (POC) is generated by phytoplankton photosynthesis, is recycled in the food web, and is eventually exported out of the euphotic layer by the biological pump (Ducklow et al., 2001; Mestre et al., 2018). Derivation of SPM concentrations over the global waters is thus valuable and needed to answer interdisciplinary questions.

Efforts to estimate SPM from ocean color satellites are hampered by the inherent complexity of the particle optics and the relationship between remote sensing reflectance (R (λ)) and SPM. The suspended particles in natural waters vary in size, shape, refractive index, and density (Babin et al., 2003; Shi & Wang, 2019; Twardowski et al., 2001; Xi et al., 2015). The optical properties of particles such as backscattering and absorption coefficients are not solely determined by SPM. On the other hand, R (λ) is a function of the total backscattering coefficient and total absorption coefficient of all water constituents (Gordon et al., 1988). According to previous studies, the R (λ) spectra vary with the SPM levels in a complex manner (Binding et al., 2005; Doxaran, Froidefond, Lavender, & Castaing, 2002; Shen et al., 2010). In Figure 1, we use in situ measurements to elaborate the variation of R (λ) with SPM. First, in extremely turbid environments, where SPM = 100 and 1,000 mg l−1, the R (λ) spectra are significantly different from each other at the near‐infrared (NIR) bands but appear hardly separable in the blue and green bands because of intense light attenuation by particles and colored dissolved organic matter (CDOM) at short wavelengths. For satellite observation, the blue‐green R (λ) values are also likely subject to saturation at bright pixels. In contrast, in relatively clear environments where SPM = 0.01, 0.1, and 1 mg l−1, the R (λ) spectra are much more distinguishable at the blue‐green bands than in the red‐NIR domain. Such spectral responses of R (λ) make it difficult to predict the SPM over global waters from satellite measurements.

Figure 1

Variation of R (λ) spectra as a function of the suspended particulate matter (SPM) concentration. The R (λ) spectra with SPM = 1,030.56 mg l−1 and SPM = 105.23 mg l−1 were measured in the Gironde Estuary, France, on June 14, 2012 (SeaSWIR project); the spectra with SPM = 10.48 mg l−1 and 1.04 mg l−1 were obtained from the northern Gulf of Mexico on September 12–13, 2013 (GEOCAPE project); the spectrum of SPM = 0.104 mg l−1 was retrieved from the Indian Ocean on February 15, 2008 (CLIVAR‐i06s project); and the spectrum with SPM = 0.012 mg l−1 was measured in the South Pacific Gyre on November 23, 2004 (BIOSOPE project). See Table 1 for a summary of field measurements.

Table 1

In Situ Data Used for Algorithm Development

Sampling areas	Number of data	Range of SPM (mg l−1)	Radiometric approach	Experiment	Satellite matchups
Yangtze Estuary a	130	0.7–2,068	Above‐water	Yangtze	–
Scheldt River, Gironde River, Rio de la Plata a	137	48–1,400	Above‐water	SeaSWIR	–
European coasts a	99	0.7–73	Above‐water	COASTLOOC	–
Massachusetts Bay a	31	0.7–12	On‐water	VIIRS CAL/VAL	18
Gulf of Mexico a	40	0.21–30	On‐water	GEOCAPE	10
Arctic/Alaska	19	0.5–8	In‐water	ICECAPE‐2010	–
Arctic/Alaska	22	0.15–100	In‐water	MALINA	–
Mid‐Atlantic Bight	46	0.09–2	In‐water	CLIVEC‐07	–
Indian/Southern Oceans	51	0.08–1.1	In‐water	CLIVAR‐i06s	–
Hawaii	8	0.04–0.15	In‐water	MURI_HI	6
South Pacific Ocean	31	0.01–0.6	In‐water	BIOSOPE	–
All	614	0.01–2,068	–	–	34

Data were used for the algorithm development by Yu et al. (2019).

The last two decades have witnessed dramatic progress in remote sensing of SPM. A number of algorithms were developed and tested with various satellite radiometric data. Two major categories of approaches exist: semi‐analytical and empirical algorithms. The semi‐analytical algorithms rely on the radiative transfer and, to varying degrees, on the relationships established between the measurements of water's inherent optical properties (IOPs) and SPM. Despite the favorable evaluations (Chen et al., 2013; Han et al., 2016; Tavora et al., 2020), the IOPs‐SPM relationships assumed by the semi‐analytical approaches still face challenges for general applicability. The empirical algorithms describe SPM directly as a function of R (λ), bypassing the need to define the IOPs‐SPM relationships. In practice, the satellite estimation of SPM was often pursued with empirical approaches. For instance, Dogliotti et al. (2015) and Shi et al. (2018) used the R (λ) at a single NIR band for the SPM retrieval in highly turbid waters. In parallel, studies also found that the R (λ) at a red band near 645 or 665 nm was closely related to SPM in moderately turbid waters (Binding et al., 2005; Miller & McKee, 2004; Ondrusek et al., 2012; Shi et al., 2015). Additionally, the R (λ) band ratios can also work as a proxy for SPM. In this regard, Doxaran, Froidefond, and Castaing (2002) recommended using R (850)/R (550) for SPM retrieval in sediment‐dominated coastal waters. In the northern Gulf of Mexico, D'Sa et al. (2007) demonstrated that R (670)/R (555) was a good indicator of SPM concentrations. Many other variants of the empirical algorithms adopt different bands and are discussed in the review articles (Matthews, 2011; Odermatt et al., 2012). Most of the existing empirical algorithms are designed for regional applications with simple formulations. They usually perform the best for a narrow range of SPM values.

The empirical algorithms applicable to more comprehensive ranges of waters often utilize a combination of empirical functions to arrive at the SPM concentrations. The conditional algorithms belong to a special type of empirical approach and are suitable for general purposes. In a study focused on the East China Sea, Shen et al. (2010) modeled SPM as a function of R (λ) at a single band of 560 nm for low SPM (<20 mg l−1), 620 nm for moderate SPM (20–80 mg l−1), and 708 nm or 778 nm for higher SPM. Similarly, Novoa et al. (2017) proposed an algorithm that switches among green, red, and NIR bands according to the R (λ) values at a reference band 655 nm. Balasubramanian et al. (2020) presented a multi‐conditional procedure for SPM retrieval: it first estimated the particle backscattering coefficient using semi‐analytical, machine‐learning, and empirical methods; then, it derived the SPM from the pre‐determined backscattering coefficients. Yu et al. (2019) proposed a weighted SPM function incorporating R (λ) at three visible bands and two NIR bands, with impressive results. Still, their weighted algorithm lays an emphasis on R (862) for SPM > ∼50 mg l−1 and on R (671) for SPM ≤ 50 mg l−1, resulting in estimations with high uncertainty in relatively clear waters (SPM < ∼2 mg l−1). Despite the promising performance of the conditional algorithms, none of them reviewed here have been applied to global satellite R (λ) data. According to Figure 1 and previous measurements, the SPM concentrations across the open ocean, coastal, and inland waters span more than five orders of magnitude, a range higher than those of existing algorithms.

The objective of the present study is to develop a robust algorithm for the estimation of SPM concentrations in global open ocean, coastal, and inland waters. The new algorithm takes account of the satellite R (λ) measurements at NIR, red, green, and blue bands (hereafter, NIR‐RGB) and predicts the SPM concentrations with two conditional formulations. The primary tasks are three‐fold: (a) to perform a validation analysis with field measurements and satellite‐field matchups; (b) to conduct a comparison with existing algorithms; and (c) to generate mission‐long SPM products from the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi National Polar‐orbiting Partnership (SNPP). The results of the study are presented in Section 3. For demonstration purposes, in Section 4, we analyze the VIIRS Level‐2 and Level‐3 SPM products from selected coastal and inland waters. Section 5 provides discussions on the SPM algorithm performance and satellite derivations. The study is concluded in Section 6.

Materials and Methods

VIIRS Satellite Images

The satellite ocean color images used in this study were from VIIRS onboard SNPP, captured between February 2012 and September 2020. The VIIRS sensor operates from a polar orbit and covers the global waters within ∼2–3 days. It is characterized by five visible bands at 410, 443, 486, 551, and 671 nm (plus an imaging band of 638 nm), two NIR bands at 745 and 862 nm, and multiple shortwave infrared (SWIR) bands. Such a spectral setting with both visible and NIR bands is optimal to represent the variation of SPM over a high dynamic range. The Level‐1 data were recorded in a moderate spatial resolution of 750 m. The Multi‐Sensor Level‐1 to Level‐2 (MSL12) ocean color data processing system was used to convert the VIIRS Level‐1 images to the Environmental Data Records (EDR) or Level‐2 ocean color products. Specifically, the MSL12‐SWIR utilized two SWIR bands of 1,238 and 1,601 nm to estimate aerosol types and contributions in turbid environments (Wang, 2007; Wang & Shi, 2007; Wang et al., 2009). The MSL12‐NIR procedure used the NIR bands at 745 and 862 nm (Gordon & Wang, 1994) for the atmospheric correction in non‐turbid environments. In addition, the present processing adopted an improved NIR reflectance correction (Jiang & Wang, 2014). The resultant science‐quality products include R (λ) at both visible and NIR bands and accompanying products such as the chlorophyll‐a (Chl‐a) concentration, quality measures, and attenuation coefficient for downwelling irradiance.

Data for Algorithm Development

In situ SPM and R (λ) data were used for the algorithm development (Table 1). Part of the in situ data set aligns with Yu et al. (2019). The Yu et al. (2019) data were collected mainly from environments with SPM > 0.2 mg l−1, including the Yangtze Estuary and adjacent waters, Río de La Plata, the European coasts, and the northern Gulf of Mexico, in which the SPM concentrations vary between ∼0.2 and 2,000 mg l−1, with about 90% of them higher than 1 mg l−1. The present study retrieved additional data from the Sea‐viewing Wide Field‐of‐view Sensor (SeaWiFS) Bio‐optical Archive and Storage System (SeaBASS) (Werdell & Bailey, 2002), covering the Mid‐Atlantic Bight, the Arctic waters off Alaska, the North Pacific Ocean off Hawaii, and the South Pacific Ocean. The corresponding SPM data vary from ∼0.01 mg l−1 to >100 mg l−1, with more than 95% of them below 5 mg l−1.

The SPM concentrations were determined with water‐filtering protocols by the principal investigators of every experiment. The methods involved in the SPM data processing, including types of filters, volumes of waters, and balance precision, are variable among experiments. Interested readers can find a more detailed description of the measurement protocols for each experiment from SeaBASS or refer to literature (e.g., Neukermans et al., 2012). Regardless, this study adopted the coastal SPM measurements either from zero‐meter depth or 5‐m depth, depending on availability. In rare cases, we obtained a mean concentration from measurements at zero and five meters. The open ocean SPM concentrations, especially off Hawaii and in the South Pacific Ocean, were measured at multiple discrete depths, including 3, 5, and/or 15 m, which were averaged for subsequent applications. Note that the SPM sampling depths are all limited within the upper layer of the first optical depth visible to satellites. Thus, the in situ SPM data are representative of the mean concentrations of particulates for surface waters.

Various radiometric approaches were exercised for the observation of water‐leaving radiance (L (λ)). The above‐water radiometry was used in the Yangtze Estuary experiments (Yu et al., 2019), the SeaSWIR experiments (Knaeps et al., 2018), and the COASTLOOC experiments (Babin et al., 2003). All the SeaBASS radiometric data used in this study were observed with in‐water radiometry. Some in‐water radiometric products were distributed in the format of R (λ). Others were only available in the format of the depth profiles of the upwelling radiance (L (λ,z)), downwelling irradiance (E (λ,z)), and above‐water downwelling irradiance (E (λ)). We processed these profiles following the NASA‐recommended protocols. Specifically, the L (λ,z) and E (λ,z) were first normalized by the corresponding E (λ). The diffuse attenuation coefficient for L (λ,z) (K (λ)) was determined with the least squares fit to log‐transformed L (λ,z) profiles. The K (λ) was then used to extrapolate the upwelling radiance profile back to right below the water surface to obtain L (λ,0−). With that, L (λ) was finally determined as L (λ) = 0.54 L (λ,0−), where the constant 0.54 is the upwelling radiance transmittance. The on‐water radiometry was used to determine R (λ) in the Massachusetts Bay and Gulf of Mexico. The protocol for on‐water data processing was provided elsewhere (Wei et al., 2021).

The R (λ) data were spectrally interpolated onto 410, 443, 486, 551, 671, 745, and 862 nm. In some rare cases, the R (λ) products had missing R (745) and/or R (862) values. For those data, we estimated R (NIR) following an IOP‐based procedure initially used by Yu et al. (2019). It should be emphasized that the R (NIR) data in these waters are often far smaller than R (671). Thus the errors associated with the estimated R (NIR) can be ignored in the present application (see discussion in Section 5.2).

We inspected the R (λ) data quality with the quality assurance (QA) model (Wei et al., 2016). The model assigned a QA score to each individual R (λ) spectrum based on the spectral similarity between the in situ R (λ) and the reference spectra at five nominal wavelengths (∼412, 443, 488, 551, and 670 nm) (Wei, Yu, et al., 2020). The highest QA score was set to one and the lowest to zero. The median QA score of the R (λ) spectra included in this study is equal to 0.8.

The NIR‐RGB SPM Algorithm

The newly developed NIR‐RGB SPM algorithm divides global waters into two groups demarcated by an empirical threshold value of R (671): “clear waters,” where R (671) < 0.0012 sr−1 and “turbid waters,” where R (671) ≥ 0.0012 sr−1. As shown in Figure 2a, when R (671) ≥ 0.0012 sr−1, the R (λ) spectra exhibit maximal spectral values in the green‐red‐NIR domain. When R (671) < 0.0012 sr−1, the highest spectral values of R (λ) are often present in the blue bands. According to the scatter plot in Figure 2b, R (671) increases with SPM over the full range of the SPM data (∼0.01–2,100 mg l−1). R (671) is relatively less sensitive to the variation of SPM when SPM is either extremely low (<∼0.2 mg l−1) or high (>50 mg l−1). The empirical threshold at R (671) = 0.0012 sr−1, where SPM ≈ 2 mg l−1, was somewhere in the middle of the R (λ)‐SPM relationship depicted in Figure 2b. This R (671) value used for the conditional procedure is close to that of Novoa et al. (2017), which defined R (655) < 0.0022 sr−1 for waters with SPM < 8.5 mg l−1.

Figure 2

In situ measurements used for algorithm calibration for (a) mean R (λ) spectra and the range of variation (denoted as shaded areas) and (b) SPM as a function of R (671).

We adopted the SPM algorithm of Yu et al. (2019) for turbid waters where R (671) ≥ 0.0012 sr−1. This algorithm was initially developed and tested with measurements from the Yangtze Estuary, SeaSWIR, COASTLOOC, Gulf of Mexico, and Massachusetts Bay (highlighted in Table 1). It estimates SPM with the two NIR bands of 745 and 862 nm and three visible bands of 486, 551, and 671 nm, with the following function: where B 1 = 20.43 mg l−1, B 2 = 2.15, and I is a band‐ratio term, which is modeled as where the empirical constants C 0, C 1, C 2, and C 3 are equal to 0.04, 1.17, 0.4, and 14.86, respectively; w represents a weighting factor, with where the subscript i varies from 1 to 3, with λ 1, λ 2, and λ 3 corresponding to 671, 745, and 862 nm, respectively. According to Yu et al. (2019), R (862) dominates the weighted function of Equation 2 for SPM ≥ 50 mg l−1, while R (671) matters for SPM < 50 mg l−1.

For R (671) < 0.0012 sr−1, or clear waters, we found that the following quadratic function can well characterize the variation of SPM where the model coefficients α 0 = 0.5192 mg l−1, α 1 = 0.9278 mg l−1, and α 2 = 0.4291 mg l−1 were determined through non‐linear optimization.

To ensure a smooth transition of the SPM derivations across R (671) = 0.0012 sr−1, we slightly adjusted the boundary of application for Equation 4 to R (671) < 0.0008 sr−1. Then we introduced a linear smoothing function to account for the data falling within 0.0008 ≤ R (671) < 0.0012 sr−1, with where β is a weighting factor, with β = 2,500 × [R (671) − 0.0008].

Data for Algorithm Validation

An independent in situ data set, including matched SPM and R (λ) measurements (N = 153), was acquired from SeaBASS for algorithm validation (Table 2). These measurements were observed in the San Francisco Bay, Chesapeake Bay, Arctic waters, East China Sea, and Great Lakes. The SPM data vary between 0.04 and 134 mg l−1, about 28% of which represent relatively clear waters with SPM ≤ 2 mg l−1. This validation data set is sufficient for the present purpose as it is primarily intended for algorithm evaluation in low‐SPM waters; the Yu et al. (2019) algorithm was already evaluated for relatively high‐SPM concentrations.

Table 2

In Situ Data Used for Independent Algorithm Validation

Sampling areas	Number of data	Range of SPM (mg l−1)	Radiometric approach	Experiment	Satellite matchups
California/San Francisco	82	2.3–134	In‐water	RIO‐SFE	–
Chesapeake Bay	43	3.6–36	In‐water	GEOCAPE	–
Arctic/Alaska	28	0.04–2	In‐water	ICESCAPE‐2011	–
East China Sea	9	0.9–2.1	In‐water	GOCI‐2013	–
Great Lakes	14	0.25–0.58	In‐water	LAKESUPERIOR	–
All	153	0.25–134	–	–	–

We also created matchups between the in situ observations and the VIIRS overpass. A time constraint of ±3 h was followed to perform the data search and matching. The VIIRS R (λ) data were derived as the mean values over a 3 × 3 pixel neighborhood, where the coefficient of variation of R (λ) measurements was small (Bailey & Werdell, 2006). The quality flags regarding land, high sun glint, high radiance, high sensor‐zenith angle, stray light, low radiance, and cloud/ice were considered. The final field‐satellite matchups contain 34 pairs of data, recovered from the Massachusetts Bay, Gulf of Mexico, South Atlantic Bight, and Hawaii (indicated in Tables 1 and 2).

Global SPM Level‐2 Products, Level‐3 Data Binning, and Trend Analysis

The NIR‐RGB SPM algorithm was implemented with the VIIRS Level‐2 ocean color imagery. For global analysis, the Level‐2 SPM products were merged into Level‐3 daily, 8‐day, monthly, yearly, and climatology products. The Level‐3 processing algorithm is generally the same as those used for the Moderate Resolution Imaging Spectroradiometer (MODIS) global ocean color products (Campbell et al., 1995; IOCCG, 2004). They were produced with a spatial resolution of 9‐km for subsequent applications on regional and global scales. The quality flags (high sun glint, high sensor‐zenith angle, high solar‐zenith angle, and cloud shadow/stray light) were taken into account for the spatial and temporal binning. The VIIRS Level‐2 and Level‐3 SPM products are currently accessible from the NOAA Ocean Color Team website (https://www.star.nesdis.noaa.gov/socd/mecb/color).

To understand the temporal variability, we evaluated the trends in the SPM concentration based on the monthly 9‐km composites. First, the monthly SPM residuals were computed by removing the seasonal climatology from the time series. Then, the statistical analysis for the trends was carried out with a MATLAB function “fitlm.”

Case Analyses

To demonstrate the SPM products at local scales, we analyzed the SPM distribution and variations over selected coastal and inland waters (see Table 3). The coastal sites of interest include the Massachusetts Bay, Long Island Sound, Delaware Bay, Chesapeake Bay, and Pamlico Sound. Among them, the Pamlico Sound is a lagoon separated from the Atlantic Ocean by barrier islands, while the others are semi‐enclosed estuaries or embayments. Riverine inputs of sediments remain important contributors to these water bodies. Also, bottom resuspension and erosion of shorelines are also significant sources. The studied inland waters consist of Lake Erie, Lake Ontario, Lake Winnipeg, and Lake Taihu. Among them, Lake Erie, Lake Winnipeg, and Lake Taihu are relatively shallow environments; the overall turbidity thereof is determined by wind‐induced resuspension and watershed inputs. In contrast, the mean depth of Lake Ontario is over 80 m. Lake Erie, Lake Ontario, and Lake Taihu were extensively investigated (Binding et al., 2010; Shi et al., 2015). Lake Winnipeg is the second‐largest watershed in Canada, which was rarely covered in previous SPM analyses. These selected water bodies span very clear, moderately turbid, and highly turbid environments.

Table 3

Satellite Level‐2 Images Selected for Case Studies

	Level‐2 images	Time series locations
Massachusetts Bay	V2015059180648_NPP_SCINSW_L2	42.3401°N, −70.6315°W
Long Island Sound	V2015059180648_NPP_SCINSW_L2	41.0982°N, −72.9735°W
Delaware Bay	V2015044174738_NPP_SCINSW_L2	39.0561°N, −75.1397°W
Chesapeake Bay	V2015044174738_NPP_SCINSW_L2	37.6512°N, −76.0842°W
Pamlico Sound	V2015016181228_NPP_SCINSW_L2	35.3867°N, −75.8171°W
Lake Erie	V2017267175535_NPP_SCINSW_L2	42.0662°N, −81.4691°W
Lake Ontario	V2017267175535_NPP_SCINSW_L2	43.6793°N, −77.8621°W
Lake Winnipeg	V2017251193913_NPP_SCINSW_L2	53.4300°N, −98.4327°W
Lake Taihu	V2017041045756_NPP_SCINSW_L2	31.2395°N, 120.1012°E

Note. VIIRS Level‐3 SPM time series data are extracted at the indicated sampling locations.

In this study, we chose the Level‐2 images under clear skies to illustrate the spatial variability of the satellite SPM products in these regional bodies of water. The SPM time series data from the Level‐3 9‐km mapped products were used to demonstrate the temporal variability. All time series data were extracted at a central location within each body of water. Table 3 lists the satellite Level‐2 files and the coordinates for the extraction of time series.

Uncertainty Measures

We used the median absolute percentage difference (MAPD), bias, and median absolute difference (MAD) to quantify the algorithm performance in estimating SPM from R (λ) data. Specifically, MAPD is defined as where E and M refer to the estimated and measured SPM values, respectively, with N equal to the total number of data points involved in the computation. Accordingly, the bias or relative percentage difference is described as

The median absolute difference is determined as

Note that the medians, instead of means, are calculated in Equations (6), (7), (8) to reduce the impact of data outliers. We also calculated the correlation coefficient (R 2) for log‐transformed E and M for performance evaluation.

Results

Algorithm Validation

We first present the result of the NIR‐RGB algorithm using the in situ data described in Table 1. This evaluation can provide first testing for the algorithm reliability. The scatter plot in Figure 3a compares the estimated SPM concentrations with in situ measurements. The agreement is reached between modeled and measured SPM data over five orders of magnitude, where MAPD = 39%, bias = −21%, and R 2 = 0.96 (see statistics in Table 4). As far as the waters where R (671) ≥ 0.0012 sr−1 are concerned, the estimated SPM values agree with the measurements where MAPD = 37%, bias = −20%, and R 2 = 0.94. For clear waters where R (671) < 0.0012 sr−1, the SPM values estimated by the NIR‐RGB SPM algorithm experienced a slight increase in MAPD (= 41%) and a decrease in R 2 (= 0.69). The decreased R 2 is partially explained by the fact that much of the data are clustered over a narrow range between 0.5 and 1 mg l−1. Noticeably, the clear‐water SPM data in this analysis have a small median value of 0.38 mg l−1; the absolute difference of the estimated SPM from in situ measurements is small, with MAD = 0.13 mg l−1. This difference corresponds to a relative percentage difference of ∼30%.

Figure 3

Comparison of the NIR‐RGB algorithm estimated and in situ measured SPM concentrations for (a) calibration of the NIR‐RGB algorithm with the field data described in Table 1 and (b) validation of the NIR‐RGB algorithm with the field measurements given in Table 2. The straight‐line refers to 1:1.

Table 4

Calibration Statistics for the NIR‐RGB SPM Algorithm Based on the In Situ Measurements Described in Table 1

	Range of SPM (median) (mg l−1)	N a	MAPD	MAD (mg l−1)	Bias	R 2
R rs (671) ≥ 0.0012 sr−1	0.10–2,069 (66)	406	37%	18	−20%	0.94
R rs (671) < 0.0012 sr−1	0.012–14 (0.38)	195	41%	0.13	−23%	0.69
All	0.012–2,069 (6.3)	601	39%	3.2	−21%	0.96

Note. Two subgroups of data are considered: R (671) ≥ 0.0012 sr−1 (representing turbid water) and R (671) < 0.0012 sr−1 (representing clear water).

The number of data used for comparison.

We further evaluated the algorithm performance with the independent field measurements. As shown in Figure 3b, the SPM values estimated from the NIR‐RGB algorithm are highly consistent with the measurements over the range of the data between 0.04–134 mg l−1. Two data sets are closely scattered about the 1:1 line, with a high correlation coefficient of R 2 = 0.96 and small errors with MAPD = 35% and bias = −19% (see Table 5). For waters of R (671) ≥ 0.0012 sr−1, the estimated SPM data are subjected to minor errors with MAPD = 29% and bias = −22%. Within clear waters where R (671) < 0.0012 sr−1, the estimated SPM experienced an increase in MAPD (= 44%) as a result of the small SPM values. It is emphasized that the absolute differences between estimated and measured SPM values in these clear waters are indeed minor, with MAD = 0.14 mg l−1. These results strongly suggest that the NIR‐RGB algorithm has a high performance for generating reliable SPM data over a wide range of waters.

Table 5

Validation Statistics for the NIR‐RGB Algorithm Based on In Situ Measurements Described in Table 2

	Range of SPM (median) (mg l−1)	N a	MAPD	MAD (mg l−1)	Bias	R 2
R rs (671) ≥ 0.0012 sr−1	1.13–134 (10)	126	29%	3.10	−22%	0.86
R rs (671) < 0.0012 sr−1	0.044–13 (0.38)	50	44%	0.14	15%	0.82
All	0.044–134 (8.2)	176	35%	1.74	−19%	0.96

Note. Two subgroups of data are considered: R (671) ≥ 0.0012 sr−1 (representing turbid water) and R (671) < 0.0012 sr−1 (representing clear water).

The number of data used for comparison.

The scatter plot in Figure 4 compares the estimated SPM values from satellite R (λ) and the in situ measured SPM data. A strong correlation is observed between these two data sets, with R 2 = 0.86 (see Table 6). The overall uncertainties of the estimated SPM are acceptably small, with MAPD = 35% and bias = −23%, indicative of robust retrieval of SPM. For R (671) < 0.0012 sr−1, our analysis shows the statistics with MAPD = 32% and bias = −23%. The MAD value is equal to 0.24 mg l−1, less than one‐third of the median SPM concentration (= 0.83 mg l−1). For R (671) ≥ 0.0012 sr−1, the satellite‐derived SPM data are subjected to a larger MAPD of 56% and a small correlation coefficient, partially due to the narrow range of the SPM data used therein. We note that the relatively high SPM values in Figure 4 were from dynamic nearshore environments in the northern Gulf of Mexico. The satellite and in situ observations with the ±3‐h time difference may not represent each other perfectly.

Figure 4

Comparison of Visible Infrared Imaging Radiometer Suite (VIIRS)‐derived SPM data with the in situ measurements. The straight‐line refers to 1:1.

Table 6

Validation Statistics for Satellite‐Derived SPM Data

	Range of SPM (median) (mg l−1)	N a	MAPD	MAD (mg l−1)	Bias	R 2
R rs (671) ≥ 0.0012 sr−1	1.1–12.9 (9)	10	56%	3.29	−22%	0.01
R rs (671) < 0.0012 sr−1	0.047–9.3 (0.83)	24	32%	0.24	−23%	0.87
All	0.047–12.9 (1.2)	34	35%	0.40	−23%	0.86

Note. Two subgroups of data are considered: R (671) ≥ 0.0012 sr−1 (representing turbid water) and R (671) < 0.0012 sr−1 (representing clear water).

The number of data used for comparison.

It is acknowledged that the global waters are composed of more than a dozen optical water types (OWTs) (Wei et al., 2016). The potential dependence of the SPM predictions on the OWTs constitutes an important aspect of the algorithm performance. To address this issue, we divided all in situ SPM and R (λ) measurements included in Tables 1 and 2 into a total of 23 OWTs following the scheme of Wei et al. (2016). In Figure 5a, OWTs 1–23 are arranged in descending order of R (443)/R (551) ratios, indicative of increasing Chl‐a from OWT 1 to OWT 23. According to the illustration, SPM concentrations increase from OWT 1 (median ∼0.03 mg l−1) to OWT 10 (median ∼10 mg l−1) and to OWT 20 (median ∼100 mg l−1). Figure 5b further shows the MAPD and MAD for the estimated SPM with regard to each OWT. It is conspicuous that the MAPD remains insensitive to the change of the water types or SPM. For instance, they are equal to ∼40% for OWTs 1–19 where the median SPM concentrations vary between 0.03 and 20 mg l−1, and only slightly decrease to below 35% for OWTs 20–23. This conservative variation of the MAPD supports the general applicability of the NIR‐RGB SPM algorithm to global waters. The MAD, as anticipated, generally increases from OWT 1 to OWT 23, a result of the increasing SPM from OWT 1 to OWT 23.

Figure 5

Data numbers and the corresponding median SPM values of the field measurements (upper panel) and the uncertainties for estimated SPM (lower panel) among different water types. The optical water types are described in Wei et al. (2016).

Comparison With State‐of‐the‐Art Algorithms

We conducted a comparative analysis of the NIR‐RGB algorithm, Yu et al. (2019), and other previously developed algorithms (Balasubramanian et al., 2020; Dogliotti et al., 2015; D'Sa et al., 2007; Han et al., 2016; Novoa et al., 2017; Tavora et al., 2020). This independent evaluation was based on the field measurements from Table 2 and the satellite‐field matchup data indicated in Table 1. The statistical results for the SPM estimation are given in Table 7. Compared with Yu et al. (2019), the NIR‐RGB algorithm has estimated SPM with much reduced MAPD and MAD values (by almost a factor of two) in clear waters where R (671) < 0.0012 sr−1; this uncertainty reduction has led to a significant increase in the overall performance of the NIR‐RGB algorithm over the SPM range of test data (0.044–134 mg l−1). The Han et al. (2016) and Dogliotti et al. (2015) algorithms were initially proposed for band settings different from VIIRS. Hence, our evaluation adopted new bands and new coefficients from Yu et al. (2019) to accommodate such differences. The estimated uncertainties for Han et al. (2016) and Dogliotti et al. (2015) are found to be much higher, especially in relatively clear waters. The Tavora et al. (2020) algorithm is a semi‐analytical approach; the relatively large uncertainties found here could be related to the default band setting considered (which only used the red and NIR bands) and the differences between the optical parameters used (e.g., the mass‐specific absorption coefficients) and their true values. Both the Novoa et al. (2017) and D'Sa et al. (2007) algorithms are subject to larger uncertainties than those from the NIR‐RGB algorithm. Such a phenomenon is partly because these empirical algorithms were designed for narrower ranges of SPM concentrations (see Table 7). The Balasubramanian et al. (2020) algorithm generated SPM values with more decent accuracy, except that about 20 data points were estimated as negative values and excluded for the present analysis. An updated version of the Balasubramanian et al. (2020) algorithm may solve this problem (Nima Pahlevan, personal communication). In conclusion, the comparisons confirm that the NIR‐RGB algorithm has the highest performance in either clear or turbid water or in the full range of the testing data.

Table 7

Performance Comparison of the NIR‐RGB Algorithm and Several Existing Algorithms

SPM algorithm	Evaluation data		Algorithm performance
	R rs (671) (sr−1)	SPM (mg l−1)	MAPD	MAD (mg l−1)	Bias	R 2	N
NIR‐RGB	≥0.0012	≥2	29%	3.03	−22%	0.72	130
	<0.0012	<2	41%	0.17	7%	0.71	74
	∼10−5–10−2	0.044–134	35%	1.18	−20%	0.90	204
Yu et al. (2019)	≥0.0012	≥2	29%	3.03	−22%	0.72	130
	<0.0012	<2	75%	0.30	−68%	0.43	74
	∼10−5–10−2	0.044–134	42%	1.43	−29%	0.84	204
Han et al. (2016)	≥0.0012	≥2	71%	6.48	−71%	0.51	130
	<0.0012	<2	85%	0.33	−81%	0.31	73
	∼10−5–10−2	0.044–134	74%	4.49	−73%	0.80	203
Balasubramanian et al. (2020)	≥0.0012	≥2	39%	3.69	−14%	0.33	109
	<0.0012	<2	47%	0.14	−29%	0.27	72
	∼10−5–10−2	0.044–134	40%	1.60	−21%	0.74	181
Tavora et al. (2020)	≥0.0012	≥2	83%	8.24	−83%	0.47	130
	<0.0012	<2	80%	0.35	−80%	0.49	67
	∼10−5–10−2	0.044–134	82%	5.98	−82%	0.81	197
Dogliotti et al. (2015)  a	≥0.0012	≥2	69%	5.95	−69%	0.52	130
	<0.0012	<2	85%	0.33	−81%	0.31	73
	∼10−5–10−2	0.044–134	73%	4.33	−73%	0.81	203
Novoa et al. (2017)  b	≥0.0012	≥2	50%	5.50	−47%	0.47	130
	<0.0012	<2	79%	0.33	61%	0.11	74
	∼10−5–10−2	0.044–134	58%	2.51	−25%	0.73	204
D'Sa et al. (2007)  b	≥0.0012	≥2	31%	2.39	−11%	0.70	130
	<0.0012	<2	68%	0.46	44%	0.24	73
	∼10−5–10−2	0.044–134	39%	1.57	−3%	0.78	203

Note. The evaluation data consisted of both the field measurements of Table 2 and satellite‐field matchups of Table 1. A total of 130 data points are included for turbid waters where R (671) ≥ 0.0012 sr−1 and 74 data points for clear waters where R (671) < 0.0012 sr−1.

The Dogliotti et al. (2015) algorithm was recalibrated by Yu et al. (2019).

These algorithms were designed for specific ranges of SPM concentrations: 2.6–1579 mg l−1 for Novoa et al. (2017) and 0.46–16 mg l−1 for D'Sa et al. (2007).

VIIRS‐Derived Global SPM Products

In Figure 6, four monthly 9‐km composites show the spatial and seasonal variability of the SPM concentration. Note that, for intuitive illustration, we set an upper limit of 500 mg l−1 to the color bar in Figure 6. The maximum SPM values derived from VIIRS are actually much higher than 500 mg l−1, such as over the Yangtze River Estuary and Amazon River Estuary. Nevertheless, on a global scale, the number of extremely high‐SPM pixels only accounts for a small percentage of the total pixels (Figure A1).

Figure 6

Global monthly SPM composites derived from VIIRS imagery in 2013 for (a) January, (b) April, (c) July, and (d) October. Black pixels indicate no valid SPM products.

The open oceans (>1,000 m depth) are characteristic of low‐SPM concentrations (<∼1 mg l−1), which is expected because they are far away from terrestrial sources. The lowest SPM values are found in the subtropical gyres. According to the brief statistics of example data (Table 8), SPM concentrations in the subtropical ocean gyres vary over a small dynamic range, a result of strong water‐column stratification and low primary production (Behrenfeld et al., 2006). For example, the sample data from the South Pacific Gyre (SPG) vary between 0.018 and 0.05 mg l−1, with a median value of 0.024 mg l−1. Away from the gyre centers, SPM concentrations slightly increase. In the Southern Ocean, for instance, the SPM concentrations can increase up to ∼0.5 mg l−1.

Table 8

Variability of the VIIRS‐Derived SPM Concentrations in the Open Ocean, Coastal, and Inland Waters

	Latitude (°)	Longitude (°)	SPM (mg l−1)			Trend (mg l−1 month−1)	p‐value	N a
			Median	min	max
Open ocean
North Pacific	20.0241	162.3093	0.042	0.025	0.057	−0.0018	3.0 × 10−19	104
South Pacific	−26.3822	−116.9538	0.024	0.018	0.050	−0.0003	5.3 × 10−2	104
North Atlantic	26.6466	−53.3904	0.041	0.026	0.066	−0.0015	2.9 × 10−17	104
South Atlantic	−14.1346	−25.2654	0.041	0.027	0.063	−0.0011	3.6 × 10−8	104
Southern Ocean	−49.848	22.3359	0.166	0.088	0.415	−0.0088	8.6 × 10−5	77
Coastal waters
Amazon River mouth	0.6380	−49.0616	789	63.3	2576	−13.5414	0.52	94
Baltic Sea	56.4514	19.2112	0.66	0.17	3.47	0.0289	0.15	79
Chesapeake Bay	37.6512	−76.0842	3.57	1.68	11.92	0.017	0.72	104
Delaware Bay	39.0561	−75.1397	6.58	3.69	14.56	−0.013	0.83	104
East China Sea	30.3587	125.4046	0.53	0.23	2.36	−0.0118	0.35	91
Gulf of Mexico	28.8932	−89.9349	3.62	0.12	9.76	0.1500	2.8 × 10−2	104
Long Island Sound	41.0982	−72.9735	3.07	0.56	7.31	0.022	0.50	104
Massachusetts Bay	42.3401	−70.6315	0.55	0.29	1.46	−0.022	4.0 × 10−5	104
Mekong River	9.3997	106.3468	10.1	2.47	764	−5.4293	6.5 × 10−2	100
Monterey Bay	36.8341	−121.9321	0.81	0.43	5.73	0.0059	0.82	104
North Sea	54.4794	2.8857	0.44	0.27	6.74	0.0364	0.24	80
Pamlico Sound	35.3867	−75.8171	8.80	4.40	30.22	0.109	0.38	104
Río de la Plata	−34.7309	−57.4003	76.2	23.9	430	−5.0551	1.5 × 10−2	104
Santa Barbara Channel	34.2496	−119.9320	0.48	0.23	1.09	−0.0061	0.15	104
Strait of Georgia	49.1621	−123.4295	5.85	0.96	20.1	−0.0415	0.65	96
Subei Bank	33.3237	121.1706	49.7	7.10	1311	−15.1376	6.1 × 10−3	102
Yangtze River Estuary	31.3045	122.0748	155	33.2	1069	−11.8353	0.18	87
Inland waters
Lake Erie	42.0662	−81.4691	2.04	0.47	18.5	−0.0436	0.52	93
Lake Huron	44.7434	−82.4852	0.27	0.12	0.54	−0.0060	3.6 × 10−3	96
Lake Michigan	42.6075	−87.0914	0.32	0.21	0.49	−0.0119	1.6 × 10−9	104
Lake Okeechobee	26.9606	−80.8	57.9	29.9	124	0.7151	0.27	102
Lake Ontario	43.6793	−77.8621	0.47	0.20	0.92	−0.0118	2.3 × 10−4	102
Lake Superior	47.9684	−87.2949	0.40	0.18	0.55	−0.0088	6.5 × 10−5	87
Lake Taihu	31.2395	120.1012	43.2	5.24	225	−1.9424	0.09	97
Lake Victoria	−0.9244	32.8522	2.30	0.73	9.28	−0.0753	0.15	103
Lake Winnipeg	53.4300	−98.4327	8.57	3.41	30.45	−0.34	0.09	55

Note. Results are based on monthly 9‐km composites (February 2012–September 2020). The trends in boldface represent significant results at the 95% confidence level (p < 0.05).

N refers to the total number of months with available SPM data.

The coastal oceans (≤ 1,000 m depth) are much more dynamic in terms of the SPM ranges and absolute values. Among the selected examples in Table 8, SPM concentrations vary between ∼7 and >1,000 mg l−1 in the Subei Shallow Bank (or Subei Shoal) on China's east coast and between ∼2 mg l−1 and >700 mg l−1 in the Mekong River Estuary, well above two orders of magnitude. Although the dynamic ranges of the sample data extracted for the Amazon River, the Yangtze River, and Río de la Plata estuaries are only slightly higher than two, they undoubtedly represent some of the most turbid environments (Shi & Wang, 2010). Note that it is not uncommon to observe low SPM values in coastal environments, such as the Santa Barbara Channel, where 0.2 < SPM < 1 mg l−1. As indicated by the example data, the SPM concentrations in coastal waters can span at least four orders of magnitude. The high contrast in the SPM concentrations entails using a “general” algorithm such as the NIR‐RGB algorithm for the SPM retrieval. Further demonstration of the SPM derivations in coastal waters is given in Section 4.1.

The SPM data from sizable inland water bodies are also exemplified in Table 8. Lake Erie, Lake Manitoba, and Lake Taihu demonstrate the most contrasting variation in the SPM concentrations. On average, Lake Okeechobee and Lake Taihu are among the most turbid environments. On the other hand, Lake Huron and Lake Michigan are relatively clear, with SPM < 0.5 mg l−1 in general. These examples show that the SPM concentrations in inland waters also vary over a high dynamic range from ∼0.10 mg l−1 to >300 mg l−1. Further demonstration of the SPM derivations in inland lakes is given in Section 4.2.

Figure 7 shows the monthly trends of the SPM concentrations between February 2012 and September 2020. Note that the values given in Figure 7 are all statistically significant (p < 0.05), with the red and blue hues for positive and negative trends, respectively. The white color refers to zero or insignificant trends. Overall, negative trends are prevalent in the SPM concentrations across the open oceans. The subtropical ocean gyres, including the North Pacific Gyre, the South Pacific Gyre, the North Atlantic Gyre, the South Atlantic Gyre, and the Indian Ocean Gyre, show decreasing trends around −0.001 mg l−1 month−1 (see statistics of example data in Table 8). In addition, the decreasing trends are also found in many coastal or marginal seas, such as the Gulf of Maine, Mid‐Atlantic Bight, Mediterranean Sea, and the Red Sea. In contrast, positive trends are rare and only present in a small number of regions, such as the northern Gulf of Mexico. In many coastal and inland waters, the SPM concentrations also decreased during this period.

Figure 7

Trends of the SPM concentration (mg l−1 month−1) derived from VIIRS monthly composites (from February 2012 to September 2020). White areas indicate that the trend is not significantly different from zero at the 95% confidence level (p < 0.05).

Case Studies

Along the US East Coast

In the following demonstration, example Level‐2 SPM images are shown in Figure 8. The corresponding Level‐3 SPM time series data are further depicted in Figure 9.

Figure 8

VIIRS Level‐2 SPM products along the US East Coast for (a) the Massachusetts Bay (February 28, 2015), (b) Long Island Sound (February 28, 2015), (c) Chesapeake Bay and Delaware Bay (February 13, 2015), and (d) Pamlico Sound (January 16, 2015). Symbols “∆” indicate the locations to extract SPM time series data.

Figure 9

VIIRS monthly SPM time series in (a) the Massachusetts Bay, (b) Long Island Sound, (c) Delaware Bay, (d) Chesapeake Bay, and (e) Pamlico Sound. The sampling locations for these time series data are indicated in Figure 8 and tabulated in Table 8.

As shown by Figure 8a, Massachusetts Bay is a clear coastal environment. More turbid pixels are often present in the west (e.g., Boston Harbor) and the south (e.g., Cape Cod Bay), where the riverine inputs and bottom suspension dominate. The SPM time series data (Figure 9a) reveal two seasonal maxima occurring in spring and fall, respectively, echoing the spring and fall blooms in this region (Hyde et al., 2007). The SPM values in February of both 2013 and 2015 exceeded the long‐term mean SPM (= 0.56 mg l−1) plus three standard deviations. These extremes are coincident with the historical winter storms that swept across New England in those 2 months. In addition, this SPM time series witnessed a significant decline from 2012 to 2020, with a trend of −0.023 mg l−1 month−1 (see the statistics in Table 8).

The image in Figure 8b recorded moderate SPM concentrations across the Long Island Sound. Some point sources showing abundant SPM are captured in the scene, which are likely related to dredging activities. The example time series data (Figure 9b) fluctuate between 0.5 and 8 mg l−1, with a mean of 3 mg l−1. This range of variation is comparable to historical in situ observations (Aurin et al., 2010). The maximum SPM values occurred between late fall and early winter. The two estuary systems of the Delaware Bay and the Chesapeake Bay are characteristic of high‐SPM values in the upper bay area, moderately high SPM values near eastern and western coasts, and relatively low‐SPM concentrations along the north‐south deep channel in the middle (Figure 8c). The Delaware Bay time series data span a range of 3.7–15 mg l−1 (with a mean of ∼7.2 mg l−1) (Figure 9c), and are comparable with previous in situ sampling results (e.g., McSweeney et al., 2017). The time series in the Chesapeake Bay (Figure 9d) varies from 1.73 to 8 mg l−1 (with a mean of 4 mg l−1). An extreme value was recorded in February 2015, possibly a result of increased freshwater discharge bringing sediments into the lower bay. The Pamlico Sound image revealed contrasting SPM values in the sound to neighboring waters (Figure 8d). This single scene encompassed a wide range of SPM concentrations (about three orders of magnitude). It also clearly captured the sediment plumes escaping into open waters to the east. The time series data exhibit characteristic seasonality (Figure 9e), with the maxima in winter and the minima in summer, which is similar to those in the Long Island Sound and Delaware Bay. The time series data in Figures 9b–9e show no significant trends (see statistics in Table 8).

Inland Lakes

Figure 10 presents the example Level‐2 images for selected inland waters. The monthly time series data extracted from each water body are subsequently plotted in Figure 11.

Figure 10

VIIRS Level‐2 SPM products over inland waters of (a) Lake Erie (September 24, 2017), (b) Lake Ontario (September 24, 2017), (c) Lake Winnipeg (September 8, 2017), and (d) Lake Taihu (February 10, 2017). Symbols of “∆” denote the locations where the SPM time series data were extracted.

Figure 11

VIIRS monthly SPM time series in (a) Lake Erie, (b) Lake Ontario, (c) Lake Winnipeg, and (d) Lake Taihu. The specific location of every time series is indicated in Figure 10.

Lake Erie is shown to be a turbid system, especially in its western basin, where SPM > 10 mg l−1 (Figure 10a). The extracted SPM time series reach annual maxima (>10 mg l−1) in early winter and minima (∼0.5 mg l−1) in summer (Figure 11a), consistent with the results of Binding et al. (2010). In contrast, Lake Ontario is a relatively clear environment, with SPM < 1 mg l−1 for most of the region (Figure 10b). The SPM time series in this clear water body picked up a bimodal seasonal pattern, with two SPM maxima (∼0.6 mg l−1) in spring and fall, respectively (Figure 11b). This unique seasonality is similar to that of the Chl‐a concentrations in this region (Watkins, 2009). The two time series data sets in Figures 11a and 11b both appear decreased, but only the one in Lake Ontario (by about −0.0118 mg l−1 month−1) is statistically significant (also see Table 8 and Figure 7).

According to Figure 10c, the SPM concentrations are more pronounced in the south basin of Lake Winnipeg due to its vicinity to tributary inputs. The time series data suggest that the monthly composite data are only available in summer and fall (Figure 11c); about half of the months have no valid observations due to the ice and cloudiness. Based on what we have obtained, the minimum SPM values (∼3.4 mg l−1) appear in late spring and the highest SPM concentrations (30 mg l−1) in late fall. Two exceptions existed in 2012 and 2014, which are likely related to melting snow/ice and intensified resuspension, and need further investigation. Lake Taihu is the most turbid environment of all selected regions (Figure 10d). The high SPM concentrations in the southwestern part of the lake are remarkable, where the rivers discharge a large amount of sediment. The time series from the central lake exhibits the maxima (up to ∼224 mg l−1) in winter and the minima (∼5.24 mg l−1) in summer, with a long‐term median value of ∼43 mg l−1. Episodic events with extremely high SPM values are recorded in the winter of 2013 and 2016, which could be related to the elevated rainfall, sediment discharge, and wind‐induced resuspension (Shi et al., 2015; Shi et al., 2018). The time series in Lake Winnipeg and Lake Taihu show negative but insignificant trends (see statistics in Table 8).

Discussion

Algorithm Robustness and Applicability

The newly developed NIR‐RGB SPM algorithm divides the global waters into “turbid” (R (671) ≥ 0.0012 sr−1) and “clear” (R (671) < 0.0012 sr−1) environments. With a combination of R (λ) at the NIR, red, green, and blue bands, the algorithm is capable of predicting the SPM concentrations over a dynamic range of about five orders of magnitude. This capability is essential as the SPM concentrations at a location can fluctuate dramatically from a small value to a high value as a result of sediment advection, resuspension, and primary production. Taking the VIIRS SPM products in Figure 11a as an example, Lake Erie may belong to the clear‐water regime in the summer months but change to the turbid‐water regime in winter. The NIR‐RGB SPM algorithm can cope with such dynamic situations and automatically switch from R (λ) at the NIR/red bands for the turbid regime and to R (λ) at the blue and green bands for the clear‐water regime. Subsequently, accurate SPM results can be achieved from the NIR‐RGB algorithm (Section 3.1). The better performance of the NIR‐RGB algorithm than the existing algorithms (Balasubramanian et al., 2020; Dogliotti et al., 2015; D'Sa et al., 2007; Han et al., 2016; Novoa et al., 2017; Tavora et al., 2020; Yu et al., 2019) can be ascribed to the fact that it incorporated optimal but different band settings for different water regimes.

We note that the turbid‐water pixels are located in coastal and inland regions (see Figure A2). The mid‐Atlantic area, the East China estuaries and coastal areas, the Amazon, the Arctic, the Ganges‐Brahmaputra Delta, and the Mississippi Delta all belong to this category. Although the turbid‐water pixels only account for about 5% of the global water surface area, they cover the primary land‐ocean interfaces for sediment flux to oceans (Milliman & Farnsworth, 2011). For these turbid‐water pixels, Equations (1), (2), (3) are employed to ensure the reliability of SPM retrievals. Specific to our algorithm, the open oceans and the vast coastal and much of the inland waters are treated as clear waters wherever R (671) < 0.0012 sr−1. This strategy has proved effective in generating SPM products from VIIRS. However, the NIR‐RGB algorithm does not separate optically shallow from deep waters. The R (λ) values in shallow environments such as the Great Bahamas can be high due to bottom reflection (Wei, Wang, et al., 2020). The resultant SPM products in such environments are not reliable. Although the SPM data in the Arctic and Antarctic coasts are of great value, they should be further evaluated because satellite R (λ) data are subject to increased uncertainties at high latitudes (Mikelsons et al., 2020). These and other limitations could be further addressed in the future.

Algorithm Uncertainty

Like all ocean color algorithms, the NIR‐RGB algorithm is subject to the satellite R (λ) measurement uncertainties. The function in Equation 4 utilizes the ratio of R (443)/R (551), the same blue‐green band combination with that of the Chl‐a algorithms (O'Reilly & Werdell, 2019; Wang & Son, 2016). One can deduce that the uncertainty distribution of SPM would be similar to that of the ocean Chl‐a concentrations (e.g., Moore et al., 2009). In coastal and inland waters, the atmospheric correction process may have more issues, where modeling of the aerosol properties is often more challenging than in open oceans. The resulting R (443) data from satellites often suffer substantial errors in coastal and inland waters (e.g., Wei, Yu, et al., 2020). As the NIR‐RGB algorithm involves this blue band, the quality of R (443) will play a role in the SPM uncertainty. However, we emphasize that the R (443) (and R (551)) measurement uncertainties are not a significant factor for the SPM estimation in turbid waters, where R (671) and R (862) matter. According to Yu et al. (2019), the uncertainties of R (862) and R (671) dominate the SPM estimations for SPM > 50 mg l−1 and 2 < SPM < 50 mg l−1, respectively; the weighted R (745) has a negligible contribution to the uncertainties of SPM overall. On the other hand, as R (671) is used for the separation of turbid and clear regimes, R (671) uncertainties could cause abrupt changes in the SPM estimations at the threshold value of R (671). The smoothing function of Equation 9 is devised for the NIR‐RGB algorithm to smooth transitions of SPM derivations across the threshold. The images and time series data (Figures 8–11) indeed suggest that the transitioning is a minor issue for SPM products.

SPM Products in Open Oceans

The spatial distribution of SPM in surface oceans (Figure 6) takes on a similar pattern to that of Chl‐a concentrations (O'Reilly & Werdell, 2019; Wang & Son, 2016). The similarity is as anticipated in that the suspended particles in clear ocean waters are primarily of biogenic origin and that the Chl‐a and SPM algorithms are both derived from R (443)/R (551) ratios in these waters. An implicit assumption with the use of R (443)/R (551) is that the absorption coefficient of CDOM (and detritus) covaries with that of phytoplankton. Studies showed that the variation of the CDOM and phytoplankton absorption coefficients tends to decouple more often than presumed (Siegel et al., 2005). This decoupling can partially explain the scatters and uncertainties revealed in the estimated SPM values in Figures 3 and 4. It is also noted that the low SPM values in open oceans, with 0.02–0.07 mg l−1 in the subtropical gyres, make the satellite inversion a challenging task. Indeed, in situ measurements of SPM in such environments require filtering very high volumes of waters to minimize the measurement errors. The large MAD values (>0.1 mg l−1) in our evaluations (Tables 4, 5, 6) should be interpreted as contributions from the algorithm itself and other sources, including the in situ measurement uncertainties. Nonetheless, the blue‐to‐green ratio function in Equation 4 provided a reasonable estimation of SPM over the clear waters.

Long‐Term Trends in the SPM Concentrations

The SPM trends in Figure 7 were estimated from the monthly SPM time series, on a per‐pixel basis. Note that some SPM data are missing in the time series due to invalid satellite R (λ) measurements (see example in Table 8). Missing data can cause unknown problems in the estimated trends. On the other hand, the average percentage of daily valid VIIRS R (λ) observations is only about 10% in global waters; it is higher in the subtropical gyres and in the coastal and inland regions at low‐mid latitudes (Hu et al., 2020). Some monthly SPM data were actually constructed with several or even fewer R (λ) measurements. Less representative data may introduce biases to the calculated trends. Addressing these potential problems is beyond the scope of the present study, which can be partially solved in the future by the use of merged ocean color products from multiple satellites (Liu & Wang, 2019). Nonetheless, several findings from the trend analysis are noteworthy.

First, the negative SPM trends in the open oceans resemble the previously reported chlorophyll trends (Beaulieu et al., 2013; Signorini et al., 2015). According to Signorini et al. (2015), the downward Chl‐a trends in the subtropical gyres are consistent with the biogeochemical response to the forcing factors, such as the warming surface temperature and shallowing mixed layer depth. We note that the controlling factors for the SPM variation are likely more complicated than those for Chl‐a. In addition to the forces related to phytoplankton growth, the composition and mass density of particles should play a role in determining the SPM. This complexity can help understand localized differences between the SPM trends in Figure 7 and the Chl‐a trends.

Second, the downward SPM trends in coastal and inland waters given in Figure 7 and Table 8 are also mentioned in previous regional analyses. For example, Shen et al. (2013) showed that the long‐term continuing decrease of river discharge caused an overall decline of the SPM in the Yangtze River Estuary. Binding et al. (2015) suggest a remarkable increase in water clarity in Lake Ontario, Lake Huron, and Lake Michigan. The increased water clarity could be a result of decreased SPM in these waters. It is understandable that the sources of suspended particles in coastal and inland waters are versatile and can include inputs from watersheds, bottom suspension, shoreline erosion, and primary production. Interested readers should refer to specific studies for further information.

Third, significant increasing trends do exist in the SPM monthly composites (Figure 7 and Table 8). The positive trends in the northern Gulf of Mexico, in particular, can be related to the increasing river discharges (Tao et al., 2014) and to the intensified wind forcing (Muller‐Karger et al., 2015) building up bottom suspension. The increasing trends in the SPM have not been reported before and merit follow‐up analysis.

Algorithm Compatibility With Landsat‐8

Many other multispectral satellite sensors, including the Ocean and Land Imager (OLI) onboard Landsat‐8, also measure R (λ) at the visible and NIR bands. It is possible to implement the NIR‐RGB SPM algorithm with R (λ) data from OLI. In particular, the OLI imagery allows for mapping in greater detail with 30‐m spatial resolution in coastal and inland waters (Luis et al., 2019; Pan et al., 2018; Wei et al., 2019). We tested two OLI images captured over the Mississippi River Delta and the Nantucket Sound off of Massachusetts. To account for the band differences between OLI and the NIR‐RGB SPM algorithm, we converted OLI's R (655) to R (671) by an empirical relationship:

The above function results from the least squares fit to a large set of hyperspectral measurements (Wei et al., 2016), with R 2 = 0.99. Similarly, R (561) of OLI was transformed to R (551) with the function given below (R 2 = 0.99):

The 745 nm band would not be needed in these cases as it only plays a secondary role to R (671) and R (862) (Yu et al., 2019). As shown in Figure 12, more than 90% of the pixels in the Mississippi River Delta show SPM varying between ∼2 mg l−1 and > 1,000 mg l−1. In the Nantucket Sound, about 25% of the pixels are characteristic of SPM higher than ∼2 mg l−1 (and <∼10 mg l−1). The SPM data derived from OLI and VIIRS are highly comparable in terms of their spatial patterns. Quantitatively, SPM differences between OLI and VIIRS are small, with MAPD = 26% and bias = −11% in the Mississippi River Delta, and with MAPD = 18% and bias = −8% in the Nantucket Sound.

Figure 12

The SPM concentrations in the Mississippi River delta (January 10, 2019) and the Nantucket Sound, Massachusetts (July 30, 2013). Atmospheric correction of the images in (a and c) was performed with Acolite (NIR‐SWIR method) and SeaDAS (standard method), respectively.

Conclusions

In this study, we developed a new conditional algorithm for the estimation of global SPM concentrations from satellite ocean color imagery. Validations indicated that the algorithm has a reliable performance in retrieving SPM products from ∼0.01 mg l−1 to >2,000 mg l−1. The NIR‐RGB SPM algorithm has been implemented with VIIRS data, and the mission‐long multi‐level global SPM products have been generated. We provided a brief characterization of the spatial distribution and temporal trend of the global SPM products. With case studies, we demonstrated the quality of the regional VIIRS SPM products across various waters. Despite some challenges, the new SPM data are an essential addition to the existing suite of ocean color products.

It is anticipated that the global SPM products will be of great value in many basic and applied applications. One important application is the use of SPM as a critical water quality indicator in environmental monitoring programs, especially where food security, including drinking water, is at risk. Additionally, the SPM data, particularly in the coastal and inland waters, could be incorporated into hydrodynamic models to predict the sediment transport as well as the dispersion of chemicals. In open oceans, the analysis of the SPM concentrations will convey information about the carbon cycles. Last but not least, satellite SPM products will support the classification of water resources according to their usefulness in agricultural, industrial, household, recreational and environmental activities.