An algorithm to estimate suspended particulate matter concentrations and associated uncertainties from remote sensing reflectance in coastal environments

: Suspended Particulate Matter (SPM) is a major constituent in coastal waters, involved in processes such as light attenuation, pollutant propagation, and waterways blockage. The spatial distribution of SPM is an indicator of deposition and erosion patterns in estuaries and coastal zones and a necessary input to estimate the material ﬂuxes from the land through rivers to the sea. In-situ methods to estimate SPM provide limited spatial data in comparison to the coverage that can be obtained remotely. Ocean color remote sensing complements ﬁeld measurements by providing estimates of the spatial distributions of surface SPM concentration in natural waters, with high spatial and temporal resolution. Existing methods to obtain SPM from remote sensing vary between purely empirical ones to those that are based on radiative transfer theory together with empirical inputs regarding the optical properties of SPM. Most algorithms use a single satellite band that is switched to other bands for different ranges of turbidity. The necessity to switch bands is due to the saturation of reﬂectance as SPM concentration increases. Here we propose a multi-band approach for SPM retrievals that also provides an estimate of uncertainty, where the latter is based on both uncertainties in reﬂectance and in the assumed optical properties of SPM. The approach proposed is general and can be applied to any ocean color sensor or in-situ radiometer system with red and near-infra-red bands. We apply it to six globally distributed in-situ datasets of spectral water reflectance and SPM measurements over a wide range of SPM concentrations collected in estuaries and coastal environments (the focus regions of our study). Results show good performance for SPM retrieval at all ranges of concentration. As with all algorithms, better performance may be achieved by constraining empirical assumptions to specific environments. To demonstrate the flexibility of the algorithm we apply it to a remote sensing scene from an environment with highly variable sediment concentrations.


Introduction
Suspended Particulate Matter concentration (SPM) is a major constituent in coastal waters that is involved in a variety of processes (e.g., carrying adsorbed pollutant, reflecting and absorbing light modulating its availability to planktonic and benthic organisms, clogging waterways). SPM is a necessary input in models solving the sub-surface light field and is a state variable in sediment transport and biogeochemical algorithms of coastal seas. In addition, the geographical distribution of SPM is needed to analyze the deposition and erosion patterns in estuaries and coastal zones and to estimate the material fluxes from land, through rivers, to sea. Depending on the composition of SPM (organic or inorganic) it may be indicative of availability of food of interest to the aquaculture industry (organic) or particles that are detrimental to bi-valve growth (inorganic). Estimates of the distribution of SPM are thus valuable for coastal management.
Methods to monitor SPM concentrations in coastal waters include both in-situ and remote sensing approaches. In-situ measurements, apart from being expensive, provide data that are limited in space and/or time and therefore do not always represent well the temporal and spatial dynamics of a river, estuary or coastal system. Remote sensing techniques, such as satellite ocean color remote sensing, provide spatial distributions of surface SPM concentration in natural waters not possible with in-situ tools, with spatial resolution as high as 10 m (e.g., as with the spatial resolution of the Sentinel 2a and 2b satellites) and temporal resolution as high as one hour (e.g., Geostationary Ocean Color Imager, GOCI). Both approaches, however, are most useful when done in synergy, as remote sensing estimates of SPM require ground-truth data to insure they are unbiased, and their estimates are confined to surface water (though sub-surface dynamic may be inferred [1]).
Algorithms to obtain SPM from remote sensing reflectance in coastal environments include purely empirical approaches (e.g., statistical regressions [2,3]; band ratios [4]; neural network [5,6]) and, in the past decades, semi-analytical approaches based on analytical relations between reflectance and the optical properties of SPM and on the empirical knowledge of their spectral backscattering and absorption coefficients (hence referred to as semi-analytical). Semi-analytical SPM inversions are typically performed using a single satellite band (e.g., References [7][8][9]), that may be switched for different ranges of turbidity (e.g., Reference [10]), or using the full spectrum like attempted by References [11,12]. The necessity to switch bands is due to saturation of reflectance as SPM concentration increases [13]. Under such circumstances near-infra-red (NIR) and short-wave-infra-red (SWIR) bands are usually used [14]. Those bands, however, usually perform poorly for moderate to low SPM concentrations due to lower reflectance as a result of water absorption and lower signal-to-noise ratio (SNR) than bands at shorter wavelengths.
Any approach for SPM retrievals, however, has method-specific limitations. Empirical methods, while easy to implement, depend strongly on the SPM ranges and sediment characteristics with which they have been developed. These algorithms are usually site-specific and their coefficients need to be adapted for a defined coastal domain. The performance of the semi-analytical algorithms based on the theoretical relationship between SPM and the absorption and backscattering coefficients [15][16][17][18] is limited by the accuracy of the SPM-IOP (Inherent Optical Properties) relationship, which vary with sediment characteristics ( [9], for example, size, shape, and mineralogy composition).
Here we propose a multi-wavelength semi-analytical algorithm (herein after MW algorithm) that, in addition to SPM concentration, also provides an associated uncertainty estimate and can be applied to any sensor. The uncertainty in estimated SPM concentration is based on uncertainties in both measurement and empirical inputs used. This study uses all the existing wavelengths that could provide information regarding SPM that do not exhibit significant saturation. We avoid bands where CDOM (Colored Dissolved Organic Matter) and phytoplankton are likely to be significant contributors to the signal. The algorithm is designed to be used in (and is tested for) rivers, estuaries and coastal environments that are optically deep (the bottom does not contribute to reflectance) and devoid of floating vegetation. We caution the readers that if, in their applications, phytoplankton and CDOM dominate the remote sensing reflectance signal at the bands we have retained, our algorithm may not work well as it was not tested in such waters. The approach used here could be generalized to other environments but appropriate datasets would be necessary for testing and possibly modification of the algorithm (e.g., explicitly solving for phytoplankton and CDOM), though this is out of the scope of this paper.
Assigning an uncertainty to SPM enhances user confidence and defines the range of possible applications of data products [19]. It is thus important that uncertainties in ocean color remote sensing data and in the products derived from it (e.g., SPM) are documented. Several techniques have been proposed [19][20][21], mostly for Chlorophyll or IOPs which are derived from Remote Sensing Reflectance (R rs ). Despite the importance of uncertainty estimates, no algorithm-derived SPM, to date, has been provided with uncertainties, other than can be estimated from match-ups with in-situ data. Those uncertainties, however, represent only the environment where the in-situ data comes from and may not work similarly at other locations. Here we propose a way to provide uncertainties in products derived from remote sensing reflectance.

Data Sources
Data used to evaluate our algorithm performance were collected worldwide ( Figure 1) representing a variety of water and SPM types. The datasets include in-situ measurements from different coastal and estuarine environments. The field data measurements are characterized by a wide range of SPM concentrations (four orders of magnitude of SPM concentration from approximately 0.4-3980 g.m −3 - Table 1) as well as spanning a range of particle composition and characteristics (or sources). Water absorption is of significant influence to R rs (λ) at the range of wavelengths (red, NIR, and SWIR) used to retrieve SPM and temperature contributes to it significantly in that range. For some of the datasets (e.g., Yangtze19, MCR13) we had concurrent temperature data. When in-situ temperature measurements were not available, climatological monthly sea surface temperature was used instead (https://seatemperature.info). This source of clmatological data was chosen because it has an extensive database of current and historical water temperature data around the world. The websites database uses both satellite data and in-situ observations to get reliable information of surface water temperatures.

Field Data
We have assembled a total of 420 in-situ data points with both R rs (λ) and SPM, spanning nine sites (Table 1). About 39% of these data were provided from Reference [22] (hereon Nechad15); 32% from Reference [23], (hereon Knaeps18); 12% from the RIVERCOLOR project (hereon Rivercolor14); 8% from the Mouth of the Columbia River dataset (hereon MCR13); 4% from Reference [24] (Yangtze19); and 4% from the RIVET project (NewRiver12). While SPM measurement protocol and data processing are similar across the datasets, the spectral radiance and irradiance measurements follow different approaches depending on the equipment available to the operator ( Table 2). Measurements were also made in different sky conditions for which we do not have complete information. Under such variable circumstances the assessment of quality for acquired data becomes necessary. Ruddick et al. [25] suggested the application of a quality control of above-water reflectance measurements named by the authors as 'similarity NIR reflectance spectrum'. This proposed approach is based on the theory that the shape of the reflectance spectra rarely changes in the range of 710-900 nm given the strong influence light absorption by water molecules , a w (λ), has on the water-leaving reflectance signal, ρ w (λ). Because the spectral shape is almost invariant it could be used to assess data quality. The method was applied to our six datasets resulting in the total removal of sixteen spectra not meeting the similarity spectrum criteria. This resulted in the removal of one spectra for Nechad15 for Indonesian waters, fifteen spectra for MCR13, and the removal of parts of the spectrum unusual shapes (Knapes18 for the Scheldt Estuary, Nechad15 for the North Sea region, and MCR13 for the Mouth of Columbia River).

The MW algorithm
A schematic of the MW algorithm is provided in Figure 2 (see Table 3 for symbols, definitions, and units). This Supplementary Material provides a package of Matlab functions to estimate Suspended Particulate Matter (SPM) and associated uncertainties from remote sensing reflectance at any spectral resolution or radiometer sensor.  Table 3 16 SPM uncertainty range g.m −3 σ 84, 16 Weighted SPM uncertainty range Nadir viewing angle rad ∆ ϕ Azimuth angle rad ρ w water-leaving reflectance −

Approach
The present study follows an approach linking analytically derived water-leaving reflectance, ρ w (λ), (Equation (1)), and IOPs. The link is based on the semi-analytical quasi-single scattering approximation (QSSA) (see review in Reference [26]). R rs (λ) is computed from the ratio of water-leaving radiance (L w ) and downwelling irradiance (E d ) measured above the water surface (Equation (2)): We convert the above water surface R rs (λ) to the in-water r rs (λ), for which the IOP relations have been derived following [27]: r rs (λ) is expressed as a function of IOPs (e.g., absorption and scattering) based on the relationship established by [28]: where G 1 , and G 2 (0.0949 and 0.0794, respectively) were estimated in Reference [28] and where b b (λ) and a(λ) are the backscattering and absorption coefficients, respectively. u(λ) is found by solving the quadratic Equation (4) for its only positive root. For our coastal application, we construct an algorithm under the assumption of optically deep waters capable of using all available bands from the red (630 nm) to the SWIR (1300 nm), from both multispectral and hyperspectral radiometers. We exclude the region between 670-700 nm to avoid potential contamination by chlorophyll fluorescence at sites with low SPM concentration. The motivation behind our approach is the expectation that using more bands should, in principle, provide additional pertinent information compared to a single band, resulting in a better SPM retrieval. We do not extend to wavelengths shorter than 630 nm to avoid the influence of phytoplankton absorption, which is highly variable in spectral shape, as well as the influence of CDOM (Colored Dissolved Organic Matter).

Parameterization of Inherent Optical Properties
The total absorption coefficient consists on the contribution of water, phytoplankton, CDOM, and non-algal particles-a(λ) = a w (λ) + a φ (λ) + a CDOM (λ) + a N AP (λ). Here, for wavelengths from 630 nm and longer, especially in sediment-laden waters, we neglect phytoplankton absorption (a φ ) assuming its contribution to total absorption is negligible [29]. We also neglect CDOM absorption as it decreases exponentially with increasing wavelength [30] at nearly twice the rate of absorption by non-lagal-particles [15]. In addition, CDOM absorption is assumed to have negligible influence on the SPM retrieval for sediment-dominated waters in the range from 700 to 1300 nm [8].
Water absorption, a w (λ), is a significant absorber in the range of interest [25,[31][32][33][34][35][36] and is temperature dependent (Equation (6)). In order to treat all of our datasets similarly, we use the results from recent studies (e.g., References [34,36]) covering the wavelength range of interest. We account for the temperature dependence using the Taylor series of Reference [37], applying the temperature coefficients by Reference [38]. Temperature changes by one degree Kelvin result in relative differences in a w ranging from −0.4 to 0.7% depending on wavelength.
Non-algal particles (NAP) can also contribute significantly to the IOPs and reflectance in coastal regions. NAP absorption (a N AP ) is modeled as SPM multiplied by mass-specific absorption (Equation (7)). Mass-specific absorption varies significantly less than SPM. The spectrum of particulate absorption has been found to be well approximated by an exponential shape (Equation (8)) in the visible wavelengths and a range of values for the exponent have been reported (see Table 2). Here we add a constant value at 750 nm consistent with recent reports [38]: where: where a * N AP (λ 0 ) is the mass specific absorption at the reference wavelength λ 0 . The NIR absorption offset, a * N AP (750), is added based on mounting observations using integrating spheres suggesting NAP absorption is non-zero in the NIR while not varying significantly with wavelength [38].
Backscattering (b b ) is typically assumed to be dominated by particles, b bp (λ), with contribution from seawater, b sw (λ), such that b b (λ) = b bp (λ) + b sw (λ). Here we neglect the backscattering by seawater relative to that of particles. The spectrum of particle backscattering is assumed to follow a power-law shape with a power-law exponent γ (Equation (9)). Refer to variations of γ in Table 4.
where b * bp (λ 0 ) is the mass specific backscattering at a reference wavelength.

Sensitivity of Solution to the Assumed Spectral Shape of Particulate IOPs
The solution employed here follows that of Reference [7], where at each wavelength we solve an algebraic equation for SPM, but with three important differences. (1) We solve this equation at all available wavelengths specified above (see Section 2.1); (2) We explicitly weigh solutions at wavelengths differently, dependent on uncertainties in inputs; (3) We solve the equation multiple times, changing the assumptions about IOP spectral shape and mass specific coefficients within their observed range of variability (Table 4); (4) we filter out solutions that exhibit significant saturation (see below). For each input shape parameter we use equal intervals between their observed maximum and minimum, resulting in multiple different inversion computations for each measured R rs (λ) spanning all the possible combinations of the shape parameters. The solutions are then used to obtain the most likely SPM solution and its uncertainty (e.g., References [39,40]). If the shape descriptors of absorption and backscattering are known for a given water body, this range could be significantly narrowed resulting in site-specific SPM solutions with smaller uncertainties. R rs (λ) measurements have uncertainties that are due to calibration uncertainties and methodological uncertainties (e.g., for those measured from satellites the processes of atmospheric and glint corrections). Accounting for these uncertainties in the estimation of SPM is necessary to remove compromised data as well as to weight less those wavelengths for which r rs (λ) has larger uncertainties than those with small uncertainties [45].
Uncertainties in reflectance include signal-to-noise and calibration uncertainties, and are computed as the maximum between a relative and absolute uncertainty (denoted as W below). Here, for the in-situ reflectance measurements, the absolute uncertainty, δ 1 r rs (λ) (Equation (10)), is computed based on the standard deviation of r rs (λ) of measurement replicates. When replicates are not available, the absolute uncertainty is computed based on the standard deviation of the 'noise' in the r rs (λ) spectra. The noise is computed from the difference between the measured spectrum and the same spectrum smoothed using a ten-point moving average filter. The relative uncertainty, δ 2 r rs (λ), is based on the SNR of the sensor. For in-situ sensors we assume a 5% uncertainty as stated by one of the sensors manufacturers as a typical upper bound of variability in radiative quantities (downwelling irradiance and upwelling radiance) between successive calibration (see also Reference [20] p. 44).
where √ 2 is the error propagation from the 5% uncertainty of downwelling and upwelling radiances.
where G 1 and G 2 are coefficients described in Section 2.2.1, where P 50 is the 50th percentile of solutions (SPM and mass-specific coefficients), and SPM(λ) is described in Section 2.2.5.
As SPM concentration increases, the water-leaving reflectance signal at a given wavelength tends towards an asymptotic value (for given wavelength and mass-specific IOPs, r rs (λ) is a monotonically increasing function of SPM) resulting in reflectance becoming less sensitive to SPM. This, called saturation, occurs first at visible wavelengths then progressively in the red and finally the NIR and SWIR spectral regions (when SPM reaches values > 1000 g.m −3 ) [13]. Saturation at a given wavelength occurs when the contribution of SPM to light absorption becomes predominant (i.e., the relative contributions of water to absorption is diminished). The asymptotic reflectance in this case becomes independent of SPM (e.g., Reference [13]). To avoid using wavelengths and/or choices of IOP shape parameters that result in reflectances exhibiting a significant level of saturation in reflectance, we take it into account in our inversion scheme as explained below. We compute a saturation parameter at each wavelength, λ 0 : At a specific wavelength, λ 0 , the theoretical value for a saturated u is a constant and equals To remove solutions where saturation may contribute, we remove those for which Q(λ 0 ) ≥ 0.5. The removal of solutions where the saturation threshold surpasses 50% is designed to avoid SPM underestimation while still obtaining solutions for all cases for which we have data.
IOPs (a w , a N AP , and b bp ) are not constant functions of wavelength (refer to Equations (7)- (9)). Because some sensors have varying bandwidths, for bands wider than 10 nm the band-weighted reflectance should be used: where the band-average IOP is computed for each band using its band-response function, k(λ), of bandwidth ∆k. While the spectral width of all in-situ sensors used here is narrower than 10 nm, Equation (14) should be used when using wide satellite bands (e.g., Reference [46]).

The SPM Retrieval Approach
At each wavelength and for each given spectral shape chosen for backscattering and absorption (Equations (7) and (9)), we obtain a SPM solution from u(λ 0 ) (e.g., a specific value of a N AP ,b bp , S ap and γ; [39,40]), Varying the spectral shape parameters results in a large number of solutions for which, after removal of solution for which Q(λ 0 ) > 0.5, we obtain the median solution as well as the 16th and 84th percentiles of the solutions. The difference between these two percentiles is equivalent to ±1 standard deviation for a normal distribution. To obtain our 'best' SPM solution we take the median solutions obtained at each wavelength and compute their uncertainty-weighted average: where N is the number of wavelengths for which we have SPM solutions and W the estimated uncertainty for each wavelength ( refer to Section 2.2.4). The same procedure is applied to the 16th and 84 th percentiles to provide us with the uncertainty estimate (Equations (17) and (18)).
where P 16,84 are respectively the 16th and the 84th percentiles of SPM solutions. The ( √ M) −1 factor in Equation (18) account for the fact that we are looking for an uncertainty in SPM and not the spread around it (similar to the computation of the standard error of the mean), with M the number of degrees of freedom. While we have datasets with hundreds of wavelengths of reflectance these measurements are not truly independent requiring us to investigate the information content they exhibit. Principal Component Analysis (PCA) is one method one could use to compute it (e.g., Reference [47]). we performed a PCA to establish the likely number of degrees of freedom in the input reflectances using the 'broken stick' method [48]. We applied the PCA to each reflectance dataset (normalized by the respective area under the curve) to estimate the degrees of freedom based on the number of PCA modes needed to explain the spectral variance in the data. For each dataset we computed the number of spectral modes describing >98% of the variance of the r rs spectra, and used in M in Equation (18) (varied from two to four).

The IOP Retrieval Approach
In the retrieval of SPM above, the solutions at the different wavelengths were derived independently of each other. For IOPs we look for solutions having the same spectral shapes (S ap , and γ) and mass normalized amplitudes at given wavelengths (a * N AP (443), a * N AP (750) and b * bp (700)) that provide solutions in the vicinity of the chosen SPM solution at as many wavelengths as possible. For each reflectance we find the twenty most common combinations of IOP providing SPM(λ) solutions (Equation (15)) within ±10% of the SPM. We report the median and respective ranges (maximum and minimum) for each parameter in Section 3.2.

Comparison to State-of-the-Art Algorithms
Many algorithms are available to estimate SPM from R rs (λ), especially at high SPM concentrations. Although an exhaustive comparison of all available SPM inversion algorithms is not the goal of this paper, we provide a comparison with two popular inversion models [7,10]-1. semi-analytical single band algorithm and 2. a semi-analytical switching-band algorithm. A brief description of these published algorithms is provided below:

Nechad et al., 2010 Algorithm
Reference [7] is a single band algorithm (Equation (20)), broadly applicable in turbid waters due to its semi-analytical basis and with coefficients provided for MERIS satellite data. The algorithm was developed with a range of SPM from 1.24 g.m −3 to 110.27 g.m −3 .
where A ρ and C ρ are dimensionless variables given by Reference [7] in Tables 2 and 6, respectively. Nechad10 states that λ > 680 nm shows lower relative and absolute errors for high SPM concentrations, therefore the wavelength chosen here to apply this algorithm to field data is: λ = 710 nm.
Nechad10 algorithm (Equation (20)) and the MW algorithm (Equation (16)) use the same SAA equation. The essential differences between the approaches are that the MW algorithm formulation considers a variety of shape parameters for IOPs, that it uses multiple wavelengths, and that it provides an uncertainty.

Novoa et al., 2017 Algorithm
The Reference [10] algorithm was developed based on two regional field datasets: Gironde Estuary and Bourgneuf Bay. The algorithm developed based on Gironde Estuary was chosen for our study. The dataset used to develop this algorithm ranged from of 2.6 g.m −3 to 1579 g.m −3 .
Novoa17 is a multi-conditional algorithm applicable to a broad range of low to high SPM concentrations using band-switch weights to ensure a smooth transition between the different SPM algorithms. The wavelength switching algorithm as well as the switching criteria for SPM are described in Table 5.

SPM Estimates
The MW algorithm performed well when compared to the other algorithms (see Figure 3, Tables 6 and 7 for the performance metrics) with the advantage of also providing uncertainties.
Not surprisingly, Nechad10 does not perform as well as the other algorithms in highly turbid waters for which it was not designed (as saturation will decrease its performance); Table 6. The best performance among SPM algorithms for high SPM concentration was obtained from the MW algorithm (Table 6-Yangtze19, Kaneps18, and Rivercolor14 with BIAS as low as 3%).  While the algorithms exhibit similar performance (not surprising given the similarities in their construction), the MW algorithm provides uncertainties with the SPM estimate. These uncertainty estimates vary between solutions in the range of ±7.3 to 41.9% (with a median uncertainty of about ±26.9%, Figure 4). Uncertainty estimates show similar ranges despite the differences in each environment (sediment properties), the radiometric sensor characteristics, and data acquisition method.

Estimates of IOP Parameters
The ranges of IOP parameters are constrained based on SPM solutions ( Figure 5) and reported in Table 8. Among all IOP parameters , backscattering mass-specific coefficient (b * bp (700)) is most correlated with SPM relative to the other optical parameters. The Yangtze dataset exhibited the most constrained range in IOP parameters. The table below reports the range of IOP parameters and its median solution by dataset, region, and year of data acquisition (when applicable). For example, Gironde Estuary data reported on Knaeps18 dataset, were acquired at two different dates, therefore it was split here (Gironde Estuary, FRA (2012) and (2013).

Validation of IOP Estimates
The only datasets with optical parameters available for comparison are NewRiver12 and MCR13. The two datasets do not represent a large dynamic range of SPM but are adequate for testing of the MW algorithm assessment of IOPs. The NewRiver12 and MCR13 datasets allowed for comparison between derived and measured b bp (λ), and γ. In order to compare information, b * bp (650) was selected from the SPM solution pool, and field γ was calculated based on the slope between measured b * bp (λ) at available wavelengths (see Figure 6 for b bp (650) and γ).
There is a good agreement between derived and measured backscattering for both datasets. The number of matchup measurements is small (N = 10 for NewRiver12; N = 13 to MCR13) but metrics show promising results. The shape parameter γ (Figure 6), however, had agreement. For the MCR13 data, retreived estimates of γ were lower than the measured ones while for the NewRiver12 dataset retrievals were closer to the measured of γ.

Implications and Limitations of the MW Algorithm
The MW algorithm includes several improvements over previous algorithms-1. the temperature dependence of water absorption is taken into account [25,37,38]; 2. all available wavelengths from the visible (from 630 nm to 670 nm) to the NIR and SWIR (from 700 nm to 1300 nm) are used [7,10,49]; 3. an allowance is made for an offset in absorption of non-algal particle coefficient in the NIR [38] (Equations (6) and (7)), 4. saturation of reflectance is avoided to restrict solutions [13], and 5. the retrieved SPM includes uncertainties [19]. Some previously published algorithms perform best under a relatively narrow suite of conditions [50]. The MW algorithm, however, was set-up to be used as a global algorithm allowing for a range of IOP shapes. Those ranges allow the MW algorithm to be applied to a broader range of water and sediment conditions (the MW algorithm is however not designed for locations considered optically shallow (e.g., coral reefs) or where CDOM contributes significantly to absorption in the red wavelengths included in this study). Improvement in performance for a specific body of water can be achieved by choosing spectral shapes of IOPs known to be characteristic of a given location if such information is available.
To express IOPs as a function of r rs (λ), we applied the expression developed by Reference [28]. In addition to this widely applied radiative transfer approximation, other numerical approximations exist to relate Rr rs (λ) to IOPs (see review in Reference [26]). In practice, these methods for estimating IOP cannot, individually, represent all illumination conditions and geometries, sea surface properties, and SPM properties. Therefore, the assumption of using one model to our algorithm might imply a bias. It is suggested that future work test the sensitivity to the specific R rs -IOP models chosen.
Numerous water absorption datasets exist [31][32][33][34][35][36][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66]. Some, however, show considerable differences as a function of wavelength. The observed differences may come from differences in measurement methodology (i.e., temperature control, equipment accuracy) as well as the method employed to obtain a "pure" water sample [67]. Temperature is also known to affect water absorption up to 0.7% per degree kelvin (e.g., 725-750 nm; 1145-1155 nm). Pure water absorption increases by almost a factor of 10 from the visible to the NIR (changes in temperature therefore potentially have significant impact on inversions) [13]. In conditions with high SPM concentrations, wavelengths in the visible region are expected to suffer from saturation. If saturation surpasses the suggested threshold of 50%, these saturated wavelengths are removed from providing solution for the MW algorithm which then relies on the NIR wavelengths. It is therefore expected that when temperature measurements are available, a more accurate SPM will be retrieved. We performed a sensitivity test with solutions obtained at two temperature extremes, 0 and 30 o C. SPM retrievals changed by up to 30%, with 90% of the solutions affected by 0.04% to 4.3%. It is thus recommended to take this effect into account in SPM retrievals. When such data are not available, a site-specific temperature climatology or space-based estimates of sea surface temperature could be used.
Absorption by particles is generally assumed negligible in the NIR part of the spectrum (e.g., Reference [13]). The study by Reference [38] however, shows that this might not be the case, therefore applying the suggested NIR offset would improve SPM estimates. Sensitivity tests were realized with and without the introduction of a NIR offset for non-algal particles (a * N AP (750)). The offset suggested in the study of Reference [38] for the Elbe River was added (Equations (7) and (8)) with the expectation of improving SPM estimates specially in highly turbid waters. Our sensitivity analysis demonstrate that the addition of a NIR offset lead to an overall improvement on SPM retrievals of about 2-3%. With the NIR offset, the dataset with highest SPM concentrations (Yangtze19) showed an improvement of more than 30% while the retrievals of datasets of low SPM concentration (MCR13, NewRiver12) worsen by 1-2%. Also, the addition of the offset to the total absorption by particles directly implies a higher saturation estimate (Equation (14)) leading to less solution with different combinations of IOP parameters passing the saturation threshold of 50%.
The use of spectral weights ensure that noisy NIR wavelengths do not introduce noise to solutions when SPM is low, while allowing these wavelengths to contribute more to the solution as SPM increases (and as shorter wavelengths begin to be affected by saturation). Typically, satellites have SNR of NIR wavelengths lower than the SNR in the visible. Using this weighting, the MW algorithm does not need to use discrete switching boundaries between wavelengths that could generate sharp transitions observed with some state-of-the-art switching schemes [4].
Previous studies have shown that NIR-based algorithms provide improved SPM retrieval in coastal and estuarine waters compared to visible bands only [9]. Here we use both visible, NIR, and SWIR wavelengths (up to 1300 nm) given the limitations of the data available. SWIR wavelengths (≥1300 nm), if available, could be added to our scheme, which is expected to further improve SPM retrievals [14] if weighted correctly.
In fact, the application of the MW algorithm to a Landsat8/OLI scene using a band ≥1300 nm shows promising results. Figure 7 presents the application of the MW algorithm using the red, NIR and SWIR bands (655, 865, and 1609 nm, respectively). Although the MW algorithm was only applied here to field acquired data, the MW algorithm shows great potential for applicability to ocean color sensors in conditions of large range of SPM concentration. Note, however, that the computation performed here does not include accounting for uncertainties associated with the process of computing reflectance from top-of-the-atmosphere measurements. In routine applications we recommend that such uncertainties should be taken into account (e.g., References [68,69] on how to compute them for a specific satellite). Finally, the MW algorithm may help to obtain information regarding the composition and size of SPM. Both NewRiver12 and MCR13 datasets, show high uncertainties regarding the estimate of backscattering spectral shape. The low correlation between measured and modeled backscattering spectral shapes could be related, at a first order, to the high number of degrees of freedom provided by the many different combination of mass-specific coefficients and shape parameters. Reducing these with in-situ data could help better constrain these parameter. If the uncertainty in the spectral shape of backscattering, γ, were sufficiently small, it could provide a size estimate as in Reference [17].

Performance of the SPM Algorithms
We compared the performance of the MW algorithm to two commonly used SPM algorithms using an extensive in-situ dataset covering about four orders of magnitude in SPM concentration. While all algorithms performed relatively well over the SPM range for which they were designed, they performed less well outside of these ranges. The MW algorithm had the lowest percentage error for the dataset with highest SPM concentration, 22.6%. Some of the uncertainty may be due to variability in the methods used to measure R rs (λ) among the six in-situ datasets, and different processing procedures (see different methods in Table 2). Additionally, it is not possible to have uncertainties in retrieved SPM significantly lower than the ranges observed of mass-specific IOPs of SPM, as these will directly propagate into the possible solutions. Future algorithms could take into account the observed probability distribution in mass-specific IOPs (assumed uniform here).
While the variety of data acquisition methods, spectral resolution and SPM concentration may account for differences between predicted and measured SPM values, it also highlights the robustness of the MW algorithm. The MW algorithm performed best in our comparisons which could be attributed to the fact that the MW algorithm takes advantage of the full information available in the R rs (λ) spectra, as the SPM retrieved with the MW algorithm is based on many SPM solutions.
Our calculations, however, show that the uncertainties derived by the MW algorithm (±σ SPM ) underestimate the average differences in the match-ups by about a factor of about two. The difference may be due to uncertainties in the R rs -IOP relationship used which we did not study here. Some of the dispersion may be due to the uncertainties in measured SPM which requires significant handling. Nonetheless, the order of magnitude of uncertainty is reasonable and provide a robust basis for an error estimate.

Conclusions
A semi-analytical multi-wavelength SPM algorithm using reflectance spanning the range from 630 nm to 1300 nm was developed. We demonstrated that using all available bands is useful in environments with large ranges of SPM. Comparison with state-of-the-art algorithms shows similar performance. The MW algorithm has the advantage of being easily adapted to any range of wavelengths available. This enhances applicability for extremely turbid waters such as the world's most turbid estuaries like Yangtze River and the Gironde Estuary. Application of the MW algorithm to satellite sensors shows great promise (despite the lack of a validation dataset) and has the advantage of providing a map of uncertainty (see Figure 7). It takes into account the different spectral ranges, bandwidths, signal to noise, and others characteristic parameters of different space-borne sensors. Additionally, it can incorporate uncertainties associated with computation of reflectance from space (e.g., atmospheric correction), though we did not pursue it here. Our approach may help in obtaining information regarding SPM size characteristics if the uncertainty in the spectral shape of backscattering were sufficiently small [17]. In any case, proof of this concept will require size data at all SPM ranges. Estuary Cruise. The RIVERCOLOR project for making the dataset available. The data efforts from the SOMLIT and MAGEST projects, specially Sabine Schimit. Special thank you to Tim Milligan and Jing Tao for data acquisition and processing of NewRiver12 and MCR13 datasets.Finally, we would like to thank two anonymous reviewers whose comments have significantly improved this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.