A Quantitative Comparison of Total Suspended Sediment Algorithms: A Case Study of the Last Decade for MODIS and Landsat-Based Sensors

Abstract

A quantitative comparative study was performed to assess the relative applicability of Total Suspended Solids (TSS) models published in the last decade for the Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat-based sensors. The quantitative comparison was performed using a suite of statistical tests and HydroLight simulated data for waters ranging from clear open ocean case-1 to turbid coastal case-2 waters. The quantitative comparison shows that there are clearly some high performing TSS models that can potentially be applied in mapping TSS concentration for regions of uncertain water type. The highest performing TSS models tested were robust enough to retrieve TSS from different water types with Mean Absolute Relative Errors (MARE) of 69.96%–481.82% for HydroLight simulated data. The models were also compared in regional waters of northern Western Australia where the highest performing TSS models yielded a MARE in the range of 43.11%–102.59%. The range of Smallest Relative Error (SRE) and Largest Relative Error (LRE) between the highest and the lowest performing TSS models spanned three orders of magnitude, suggesting users must be cautious in selecting appropriate models for unknown water types.

1. Introduction

The health of coastal waters not only determines the health of marine habitats in the region but also signifies the health of the nearby human inhabitants with nearly 60% of the earth’s population settled in the coastal zones of our oceans and seas, and over 90% of the world’s fish caught for consumption being sourced from coastal waters [1]. The health of water systems is typically determined from a key indicator, the water clarity (turbidity) which is influenced by the amount of dissolved matter and total suspended solids (TSS) comprising organic matter such as algae and other micro-organisms and inorganic particulate matter from minerals [2]. Monitoring TSS along with other water quality parameters is crucial for coastal ecology because TSS can directly affect the turbidity and color of water [3] and turbidity determines the amount of light availability at depth for primary production [4,5,6].

Monitoring the temporal and spatial distribution of TSS in the coastal environment can be a huge undertaking and nearly impossible in terms of financial and time resources if performed using traditional in situ water sampling methods [7] unless coupled with satellite-based remote sensing. Since the early space-borne sensors of the 1970s there has been, and continues to be, a great improvement in the spectral, spatial and temporal resolutions [8]. For example, the Landsat-based series of sensors has evolved over the years from three (red, green and blue) spectral bands with spatial resolutions of 185 m and a revisit time of 18 days to the newest Landsat-8 with 11 spectral bands (433–12,500 nm) with spatial resolutions of 30 m (and 15 m panchromatic) and a revisit time of 16 days. The shortcoming of the long revisit time for Landsat can be filled by the readily available MODIS-Aqua and Terra sensors which have shorter revisit times of one day, and with 36 spectral bands (405–14,385 nm) and spatial resolutions from 250 m to 1000 m.

Since the launch of the early remote sensing satellites in the early 1970s many studies have been conducted in remotely mapping TSS, driven in part by the capability of satellite remote sensing to cover large spatial domains in near real time [7]. Considering the past decade, remote sensing studies of the spatial and temporal mapping of TSS have utilized moderate resolution sensors including Landsat [9,10,11,12,13,14,15], MERIS [6,16,17,18,19,20], MODIS [17,21,22,23,24,25,26,27], and high resolution sensors including SPOT [28], IKONOS [29], and THEOS [30]. Further, the TSS mapping studies encompassed waters with diverse optical and physical properties, from inland lakes and river systems [21,25,31,32,33] to coastal waters [34,35,36,37,38] and from different geographical locations including America [3,32,39], Africa [26,40], Asia [16,34,41], Australia [42,43], and Europe [17,23,27,37].

The majority of the models developed in retrieving TSS by remote sensing methods are typically locally tuned to a regional water or waters with similar optical properties. Regional tuning of a TSS model is necessary because of the potentially large variation in the inherent optical properties (IOPs) of the water constituents. The theoretical basis of ocean color remote sensing has shown that sensor-measured reflectance of the water is related to the IOPs of the water—absorption and scattering coefficients. IOPs vary with the types and amounts of the water’s constituents, such as sediments, phytoplankton, detrital matter and CDOM [44] which may be different for different sediment types and phytoplankton types in different regions. In addition, factors such as water depth, viewing geometry, and atmospheric conditions all add to the complexity of the relationship between the measurement of reflectance of the water surface and the IOPs and concentrations of constituents [45].

TSS models are generally classified into three categories, (1) an empirical model where TSS is modelled directly using a statistical analysis to relate the apparent optical properties (AOPs); (2) an analytic model that relates the IOPs and AOPs of water through radiative transfer theory to derive TSS; and (3) a semi-analytic model that is partly based on the empirical analysis and grounded on the radiative transfer theory [7]. Individual TSS model designs have their own limitations and advantages. An empirical model is often sought for its simplicity and explanatory power because unique properties of local waters are tuned to each model, but it may lack general applicability. An analytic model is potentially applicable to other water bodies because it is not dependent on the in situ water constituents, but it requires accurate knowledge of water column properties which is often difficult to acquire. The semi-analytic model has both the limitations and advantages associated with the first two models, and it is generally preferred because it has higher explanatory power and is more convenient than the analytic model [7,8].

In the last decade, various TSS models have been developed [35,37,38,46,47] and applied to their respective regions with a wide range of success with reported retrieval errors ranging from lows of ~18% to highs of ~61%. Considering each model is developed and tuned for a specific region, water type and its associated IOPs, the application or transferability of the models to other regions is limited, and the likely accuracy of the results unknown. Even when an existing TSS model is applied to waters in similar regions it is often first re-calibrated before being applied. The availability of many TSS algorithms for different regions and sensors warrants one to ask if we can use someone else’s algorithm to estimate TSS in regions where we do not have any in situ observations? For the cross applicability of TSS models between different regions the design of a TSS algorithm has to either be based on analytic methods and grounded on theoretical functions of radiative transfer theory, or the waters must be assumed to have similar optical and physical properties. However, considering the vast number of TSS models that have been developed across different geographical regions with different optical and physical properties we can seek to establish the robustness in the applicability of these existing TSS algorithms for different regions.

A recent study by Brewin, et al. [44] developed an objective methodology where comparison of different bio-optical algorithms are quantitatively and qualitatively considered for use in climate studies. Following the methods of Brewin, et al. [44] and their quantitative methodology to rank the algorithms, in this study we objectively compare the performance of TSS algorithms for MODIS and Landsat sensors developed during the last decade using HydroLight simulated data for different water and sediment types. If shown to be robust, these algorithms would provide marine remote sensing scientists and coastal managers some level of confidence in their ability to assess the quality of water with minimal resource for coastal monitoring of optically unexplored waters. Specifically, this study aims to quantitatively assess the applicability of established TSS algorithms to different water types and quantify the variability in retrieving TSS when using off the shelf TSS algorithms for MODIS and Landsat sensors.

2. Materials and Methods

2.1. Dataset

2.1.1. HydroLight Simulation

A set of ocean reflectance spectra were derived using the radiative transfer numerical model HydroLight 4.2 (Sequoia Scientific, Inc., Bellevue, WA 98005, United States of America) in the four component case-2 waters mode. Using a forward model HydroLight solves radiance distributions and derives reflectance and radiance for water bodies with specific inherent optical properties (SIOPs) for given sky and water state conditions [48]. Sub-surface remote sensing reflectance’s (r_rs) were computed for infinitely deep water using a range of SIOPs, sea-state, and sky conditions. The spectral range for r_rs from HydroLight was simulated for wavelengths (λ) in the range of 400 nm–800 nm at a nominal bandwidth of 4 nm.

For all the HydroLight simulations the sea state was chosen to have a wind speed of 5 m·s⁻¹ and the sky radiance computed using the Harrison and Coombes (1988) normalized radiance model for a clear sky. The diffuse and direct sky irradiances were computed using the Gregg and Carder (1990) irradiance model for a solar zenith angle of 30° [49]. The four components, pure water, chlorophyll (CHL), colored dissolved organic matter (CDOM), and mineral (TSS) were modelled in varying concentrations, presented in

Table 1. Concentration of colored dissolved organic matter (CDOM), chlorophyll (CHL), and total suspended solids (TSS) used in HydroLight modelling. The pure water component in all the HydroLight runs remains unchanged.

CHL (mg/m3)	CDOM (m−1)	TSS (mg/L)
0.01, 3.0, 20.0	0.001, 1.0, 10.0	0.01–1.00 at 0.01 interval
		1.00–10.00 at 0.1 interval
		10.00–50.00 at 1.0 interval
		50.00–100.00 at 2.0 interval
		100.00–250.00 at 5.0 interval
		250.00–500.00 at 10.0 interval
		500.00–2000.00 at 50.0 interval
		2000.00–7000.00 at 250.0 interval

, to be representative of open ocean to turbid coastal waters. For the TSS component, five different sediment types were used, namely (1) brown earth; (2) calcareous sand; (3) yellow clay; (4) red clay; and (5) Bukata from the default database of HydroLight. The phase functions for the components were modelled as Rayleigh like phase function for pure water, Fournier-Forand phase function with b_b(λ)/b(λ) of 0.01 for CHL, and Petzold “average particle” phase function for TSS for all the aforementioned HydroLight simulations. In addition to the aforementioned parameters for HydroLight simulations, we further carried out additional simulations using the parameters outlined above but with solar zenith angles of 15°, 30°, 45°, and 60° and b_b(λ)/b(λ) ratios of 0.001, 0.01, 0.018, 0.05, and 0.1 for calcareous sand to study the robustness of TSS models to changes in solar angles and the backscattering ratios.

The IOP models used in this HydroLight simulation are described by Equations (1) and (2). The total absorption coefficient (a(λ)) is the sum of absorption coefficients of pure water (a_w(λ)), CHL (a_ϕ(λ)), CDOM (a_cdom(λ)) and TSS (a_p(λ)):

$$a(\lambda )=a_{\mathrm{w}}(\lambda )+a_{\phi }(\lambda )+a_{\mathrm{cdom}}(\lambda )+a_{\mathrm{p}}(\lambda )$$

(1)

The total scattering coefficient (b) is the sum of scattering coefficients of pure water (b_w(λ)), CHL (b_ϕ(λ)), and TSS (b_p(λ)):

$$b(\lambda )=b_{\mathrm{w}}(\lambda )+b_{\phi }(\lambda )+b_{\mathrm{p}}(\lambda )$$

(2)

The total backscattering coefficient is expressed as the sum of backscattering coefficients for pure sea water (b_bw(λ)), particulates (b_bp(λ)), and phytoplankton pigments (b_bφ(λ)).

$$b_b(\lambda )=b_{\mathrm{bw}}(\lambda )+b_{\mathrm{bp}}(\lambda )+b_{\mathrm{b}\mathsf{\phi }}(\lambda )$$

(3)

The SIOP models allow the scaling of the IOP of each component with concentration (X):

$$a_i(\lambda )={a_{\mathrm{i}}}^{*}(\lambda )\times X_i$$

(4)

$$b_i(\lambda )={b_{\mathrm{i}}}^{*}(\lambda )\times X_i$$

(5)

where i is the component and a_i*(λ) and b_i*(λ) are component specific absorption and scattering coefficients.

The SIOP of each component was either obtained from HydroLight’s default dataset or modeled using established models. For the specific absorption and scattering coefficients: the absorption coefficient for pure water was obtained from Pope and Fry (1997) [50] and mass-scattering coefficient from Smith and Baker (1981) [51], the CHL mass-specific absorption coefficient (a_ϕ*(λ)) from Prieur-Sathyendranath (1981) [52] and the CHL mass-specific scattering coefficient modeled using Equation (6), the CDOM mass-specific absorption was modeled using Equation (7) and CDOM was considered to be a non-scattering component, and the mineral mass-specific absorption and scattering coefficients were obtained from HydroLight’s default dataset for brown earth, calcareous sand, yellow clay, red clay, and Bukata. Figure 1a,b shows the mass-specific absorption and scattering coefficients of the five different minerals used in the HydroLight modelling of water reflectance.

$${b_{\phi }}^{*}(\lambda )=0.407CH{L}^{0.795}\left(\frac{600}{\lambda }\right)$$

(6)

$${a_{cdom}}^{*}(\lambda )=0.1\exp (-0.014\times (\lambda -440))$$

(7)

2.1.2. Extrapolation of Simulated Dataset

The IOP data output by HydroLight do not extend beyond 800 nm, however some of the TSS algorithms for MODIS and Landsat utilize bands beyond the 800 nm reflectance data generated by the HydroLight simulations. To include algorithms which utilize bands in the NIR region of the electromagnetic spectrum, we extrapolated the r_rs(λ) data from HydroLight to 1300 nm using Equation (1) of the quasi-analytical model of Lee, et al. [53] at a nominal wavelength of 1.0 nm:

$$r_{\mathrm{rs}}(\lambda )=g_0\left(\frac{b_{\mathrm{b}}(\lambda )}{a(\lambda )+b_{\mathrm{b}}(\lambda )}\right)+g_1{\left(\frac{b_{\mathrm{b}}(\lambda )}{a(\lambda )+b_{\mathrm{b}}(\lambda )}\right)}^2$$

(8)

where g₀ and g₁ are assigned either g₀ = 0.0949 and g₁ = 0.0794 for oceanic case-1 water [54], g₀ = 0.084 and g₁ = 0.17 for coastal water, or averaged values of g₀ = 0.0895 and g₁ = 0.1247 for coastal and case-1 waters [53]. The selection of values for g₀ and g₁ were based on the condition that the selected values provided the minimum Mean Absolute Relative Error (MARE) as defined in Equation (C1) in the Appendix C between HydroLight and Equation (8) r_rs (λ) spectra.

To model the r_rs (λ) spectra to 1300 nm using Equation (8), we used the following IOPs—the total absorption coefficient was computed using Equation (1) while the total backscattering coefficient was computed using Equation (2). Equations (4) and (5) were used to compute individual component-specific absorption and scattering coefficients using the respective component concentration and the phase function used in the HydroLight simulations as mentioned in Section 2.1.1. The total backscattering coefficient in Equation (8) was computed from the respective backscattering components in Equation (3) which in turn were computed using respective scattering components from Equation (2) and scattering phase functions and backscattering ratios discussed in Section 2.1.1. The mineral specific absorption and backscattering coefficients were spline extrapolated to 1300 nm to compute the mineral-specific absorption and backscattering coefficients required in Equations (4) and (5). The r_rs (λ) spectra generated using HydroLight and modelled using Equation (8) had MARE of 1.6% to 13.73%. The higher relative error was toward the blue end of the spectral region.

2.1.3. Grouping of Datasets

Using the extrapolation methods discussed in Section 2.1.2, in total 2.2 × 10⁴ r_rs(λ) spectra were generated for the spectral range of 350 nm to 1300 nm at the nominal wavelength of 1.0 nm for the parameters discussed in Section 2.1.1. The water, from the point of view of remote sensing, can be classified into case-1 and case-2 water types: case-1 waters are optically dominated by phytoplankton (CHL) while case-2 waters are more optically complex with varying concentrations of CHL, CDOM and TSS that are region specific [28,45]. With respect to modelling the water types, it is not feasible to model each water type that is optically similar to the optical properties of the water where each individual TSS model was developed. The TSS models that are robust enough in one region can often fail when applied to other regions because each TSS model is typically tuned to a specific region where the waters are optically unique. Thus, due to the problem of accurately modelling the waters to suit any specific TSS model, and acknowledging the fact that we cannot simulate all the conditions and compositions of ocean constituents for different regions, we resorted to five different classes (shown in

Table 2. Five different water classes.

CLASS	CDOM (m−1)	CHL (mg/m3)
I	0.01	20.0
II	10.0	0.1
III	10.0	20.0
IV	0.01	0.1
V	1.0	5.0

) to represent varying cases of water where concentrations of one ocean constituent might dominate the others or there are different degrees of contributions from each constituent. CLASS I from the water classification in Table 2 represents high CHL and low CDOM concentration which in a physical world would be associated with high phytoplankton blooms in eutrophic lakes where concentration of CHL dominates other optically active substances [55]. CLASS II with high CDOM and low CHL represents water where CDOM dominates other optically active substances, which is the case in lakes where CHL is generally low, for example as in the case in lakes in boreal regions and waters off the coast in the Baltic Sea [55]. CLASS III and IV represent the extreme cases where both CDOM and CHL are either high or low, which can be associated with high phytoplankton blooms in coastal waters for CLASS III and open ocean water with low CHL for CLASS IV. CLASS V represents a general case of coastal waters where CHL and CDOM are moderate. For all the classes of water discussed above, the TSS is varied in its concentration independent of different water cases considered.

2.1.4. HydroLight-Derived Reflectance to Sensor Equivalent Reflectance

The TSS retrieval algorithms developed by various researchers use different types of reflectance measurements to relate to TSS concentrations. The most common choice among all the TSS algorithms considered here is the remote sensing reflectance (Rrs(λ)), which is defined by Equation (9).

$$R_{\mathrm{rs}}(\lambda )=\frac{L_w(0^{+},\lambda )}{E_d(0^{+},\lambda )}$$

(9)

where Lw (0⁺, λ) is the water leaving radiance and Ed (0⁺, λ) is the downwelling irradiance evaluated above the water surface. The HydroLight generated r_rs(λ) was converted to R_rs(λ) following [53] as defined by Equation (10).

$$R_{\mathrm{rs}}(\lambda )=\frac{0.52r_{\mathrm{rs}}(\lambda )}{1-1.7r_{\mathrm{rs}}(\lambda )}$$

(10)

After converting r_rs(λ) to R_rs(λ), depending on the sensor and the bands used by particular TSS algorithms, we convolved R_rs(λ) from Equation (10) to each sensor’s respective band reflectance using the spectral response function of the sensor in their respective bands using Equation (11).

$$R_{\mathrm{rs}}^{\mathrm{k}}=\frac{{\displaystyle \underset{\Delta \mathrm{k}}{\int }{R}_{\mathrm{rs}}(\lambda )s(\lambda )}d\lambda }{{\displaystyle \underset{\Delta \mathrm{k}}{\int }s(\lambda )}d\lambda }$$

(11)

where $R_{\mathrm{rs}}^{\mathrm{k}}$ is the band averaged R_rs for each band, k, with band width Δk and spectral response function s(λ) of the sensor.

The next common reflectance type used in TSS algorithms is a normalized water-leaving reflectance which is related to R_rs(λ) as follows:

$${{\rho }_w(\lambda )|}_N=\pi R_{\mathrm{rs}}(\lambda )$$

(12)

There are also algorithms which employ normalized water leaving radiance which is calculated using Equation (13).

$${L_w(\lambda )|}_N=\frac{{\rho }_w(\lambda )\times F_o(\lambda )}{\pi }$$

(13)

where F_o(λ) is the extraterrestrial solar irradiance band averaged to each sensor’s band using their respective band spectral response functions.

2.2. TSS Models

This section lists the available TSS algorithms from 2000–2015 that are empirical and semi-analytic in their design for MODIS and Landsat-based sensors. We made an effort to select all the available TSS algorithms for the sensors considered in this study using a search database ‘Scopus’ (https://www.scopus.com/), but we acknowledge that some of the literature for TSS algorithms, which were not present in the database, might have been missed. However, within the limitation of our search capability we made an effort to use other science databases and discovered 42 MODIS empirical models and 7 semi-analytical models, 22 Landsat empirical models and 5 semi-analytical models. The summaries of each TSS algorithm are provided in

Table A1. The summary of TSS algorithms mentioned in Section 2.2.

Algorithm	Reference	Location	TSS Range (mg/L)	Bands/Algorithms	Regression Coefficient (R2)	Error	N
MOD-E1	Kumar, et al. (2016) [71]	Chilika Lagoon, India	3.9–161.7	$TSS=13181R_{rs}{(B1)}^2-1408.6R_{rs}(B1)+44.15$	0.915	RMSE = 2.64 mg/L	54
MOD-E2	Ayana, et al. (2015) [40]	Gumera catchment, Lake Tana, Ethiopia	~5–255	$TSS=2371\rho (B2)-62.8$	0.95	SE = 10.77 mg/L	54
MOD-E3	Chen, et al. (2015) [22]	Estuary of Yangtze River and Xuwen Coral Reef, China	5.8–577.2	$\begin{array}{c}log(TSS)=\frac{-b+\sqrt{b^2-4a(c-y)}}{2a}\}{R}_{rs}(B1)<0.025\\ \log (TSS)=\frac{-b+\sqrt{b^2-4a(c-y)}}{2a}\}{R}_{\mathrm{rs}}(B1)>0.025\\ a=-0.334,b=1.0046,c=0.8251,(b2-4a(c-y))\ge 0\\ y=\log (Rrs(B2))/\log (Rrs(B1))\end{array}$	0.752	RMSE = 2.1 mg/L RMSE = 38.6 mg/l	40
MOD-E4	Zhang, et al. (2016) [72] and Shi, et al. (2015) [21]	Lake Taihu, China	1.7–343.9	$TSS=9.65\exp (58.81R_{rs}(B1))$	0.70	RMSE = 14.0 mg/L	150
MOD-E5	Choi, et al. (2014) [34]	Mokpo coastal area, Korea	1.03–193.10	$TSS=1.545\exp (179.53R_{rs}(B1))$	0.92	-	96
MOD-E6	Feng, et al. (2014) [35]	Yangtze estuary	4.3–1762.1	$\begin{array}{l}TS{S}_{low}=2.49\exp \left(97.19R_{rs}(B1)\right)\}TSS<10(mg/l)\\ TS{S}_{high}=57.58\exp \left(3.48\left(R_{rs}(B2)/R_{rs}(B1)\right)\right)\}TSS>150(mg/l)\\ TS{S}_{mid}=\alpha \times TS{S}_{low}+\beta \times TS{S}_{high}\}50<TSS<150(mg/l)\\ \alpha =\frac{\ln \left(TS{S}_{low}/50\right)}{\ln \left(150/50\right)},\beta =\frac{\ln \left(150/TS{S}_{low}\right)}{\ln \left(150/50\right)}\end{array}$	0.88 (low) 0.93 (high)	RMSE = 27.7%	78
MOD-E7	Hudson, et al. (2014) [23]	Fjord in Southwest Greenland	1.2–716	$TSS=1.80\exp (19.11(R_{rs}(B1)+R_{rs}(B2))$	0.84	-	143
MOD-E8	Kaba, et al. (2014) [31]	Lake Tana, Ethiopia	~5–255	$TSS=2371\times \rho (B2)-62.8$	0.95	RMSE = 16.5 mg/L	54
MOD-E9	Lu, et al. (2014) [73]	Bohai Sea, China	~<160	$TSS=\exp \left(\frac{R_{rs}(B1)-0.0123}{0.0038}\right)$	0.75	RE ≤ 20%	627
MOD-E10	Park and Latrubesse, (2014) [32]	Amazon River system	30–150	$TSS=27.05\exp (7.83\rho (B1))$	0.88	RMSE = 6.2 mg/L	232
MOD-E11	Sokoletsky, et al. (2014) [74]	Yangtze river estuary	0–2500	$TSS=10^{\left(2\times \left[R_{rs}(B1)/R_{rs}(B4)\right]\right)}$	-		361
MOD-E12	Chen, et al. (2014) [61]	Bohai Sea	4–106.4	$\begin{array}{c}TSS=3\times {10}^8{x}^3-2\times {10}^6{x}^2+5453.3x+3.8825\\ x={\rho }_w(B15)-1.107{\rho }_w(B16)\end{array}$	0.954	RMS = 30.12%	48
MOD-E13	Cui, et al. (2013) [75]	Ponyang lake, China	0–141.9	$TSS=1.063\exp (27.859\rho (B1))$	0.91	SE = 11.20 mg/L	54
MOD-E14	Kazemzadeh, et al. (2013) [76]	Bahmanshir River, Iran	30–500	$TSS=22.92\rho {(B1)}^{0.82}$	0.63	RMSE = 261.84	23
MOD-E15	Raag, et al. (2013) [17]	Pakri Bay, Gulf of Finland	0–10	$TSS=592.54\rho (B1))+1.1136$	0.52		77
MOD-E16	Qui (2013) [46]	Yellow River Estuary, China	1.9–1896.5	$TSS=\exp \left(1.1932{\left(\frac{Rrs(B12)}{Rrs(B14)}\right)}^{0.875}\right)$	0.95	MAE = 24.5 mg/L	81
MOD-E17	Villar, et al. (2013) [77]	Maderia River	25–622	$TSS=1020{\left(\frac{Rrs(B2)}{Rrs(B1)}\right)}^{2.94}$	0.62	-	282
MOD-E18	Min, et al. (2012) [78]	Saemangeum coastal area, Korea	0.1–55	$TSS=0.24\exp (188.3R_{rs}(B12))$	0.90	-	88
MOD-E19	Ondrusek, et al. (2012) [62]	Chesapeake Bay	4.5–14.92	$\begin{array}{l}TSS=3.8813{(n{L}_w(645))}^3-13.822{(n{L}_w(645))}^2\\ +19.61(n{L}_w(645))\end{array}$	0.95	MPD = 4.2%	35
MOD-E20	Son and Wang, (2012) [39]	Chesapeake Bay	1.0–20	$TSS=1.7+5.263K_{\mathrm{d}}(490)$	0.77	STD = 0.48	15,720
MOD-E21	Wang, et al. (2012) [60]	Hangzhou Bay, China	133–1,950	$TSS=1.4599{\left(\frac{Rrs(B1)}{Rrs(B2)}\right)}^{2.3874}$	0.82		35
MOD-E22	Chen, et al. (2011) [24]	Apalachicola Bay, USA	1.29–208	$\log (\rho (B2)/\log (\rho (B1))=-0.1325\log {(TSS)}^2+0.7429\log (TSS)+0.6768$	0.86	RMSE = 4.76 mg/L	32
MOD-E23	Chen, et al. (2011) [36]	Apalachicola Bay, USA	1.29–208	$\log (\rho (B2)/\log (\rho (B1))=0.4339\log (TSS)+0.8288$	0.8	RMSE = 4.79	25
MOD-E24	Jiang and Liu(2011) as cited in [22]	Poyang Lake, China	0–40	$\begin{array}{l}TSS=1365.5x^2-369.08x+27.216\\ x=Rrs(B3)+Rrs(B4)\end{array}$	0.81	-	27
MOD-E25	Siswanto, et al. (2011) [79]	Yellow and East China Sea	0.04–340.07	$TSS={10}^{\left(0.649+25.623(Rrs(B12)+Rrs(B13))-0.646\left(\frac{Rrs(B10)}{Rrs(B12)}\right)\right)}$	0.92	RPD = 15.7%	223
MOD-E26	Zhao, et al. (2011) [80]	Mobile Bay estuary, Alabama	0–87.8	$TSS=2.12\exp (42.92\rho (B1))$	0.781	RMSE = 5.42	63
MOD-E27	Petus, et al. (2010) [81] and Petus, et al. (2014) [37]	Bay of Biscay, France	0.3–145.6	$TSS=12450Rrs{(B1)}^2+666.1Rrs(B1)+0.45$	0.97	RMSE = 61%	74
MOD-E28	Wang and Lu (2010) [25]	Yangtze River, China	45–909	$\begin{array}{l}TSS=61.369\exp (0.2623x)\\ x(\%)=\rho (B2)-\rho (B5)\end{array}$	0.78	RRMSE = 36.5%	35
MOD-E29	Wang, et al. (2010) [82]	Apalachicola Bay, USA	1–64	$\log (TSS)=1.5144\left(\log (\rho (B2))/\log (\rho (B1))\right)-0.5755$	0.72	-	16
MOD-E30	Wang, et al. (2010) [47]	Middle and Lower Yangtze River, China	75–881	$\begin{array}{l}TSS=60.24x-23.03\\ x(\%)=\rho (B2)-\rho (B5)\end{array}$	0.73	RMSE = 29.7%	153
MOD-E31	Zhang, et al. (2010) [83]	Yellow and East China Sea	0.68–27.2	$\begin{array}{l}\log (TSS)=s_0+s_1(Rrs(B4)+Rrs(B1))\\ +s_2(Rrs(B10)/Rrs(B4))\\ s_0=0.6311,s_1=22.2158,s_2=-0.5239\end{array}$	0.87	ARE = 26%	81
MOD-E32	Chen, et al. (2009) [63]	Apalachicola Bay, USA	1.29–208	$\log (\rho (B2)/\log (\rho (B1))=-0.1356\log {(TSS)}^2+0.7402\log (TSS)+0.6836$	0.853	RMSE = 5.5 mg/L	25
MOD-E33	Chu, et al. (2009) [84]	Kangerlussuaq Fjord, Greenland	~500	$TSS={10}^{\left(\frac{Rrs(B1)-1.6}{7.5}\right)}$	-	-	-
MOD-E34	Doxaran, et al. (2009L) [85]	Gironde Estuary, France	77–2182	$TSS=12.996\exp \left(Rrs(B2)/(0.189Rrs(B1))\right)$	0.89	RMSE: 18%–22%	204
MOD-E35	Jiang, et al. (2009) [86]	Taihu Lake, China	0–170	$TSS={10}^{\left(0.3568\ln (Rrs(B2)+3.3431\right)}$	0.81	ARE = 20.5%	56
MOD-E36	Liu and Rossiter (2008) as cited in [22]	Poyang Lake, China	15.6–518.8	$TSS=7167\rho (B1)-42.0$	0.91	-	25
MOD-E37	Wang, et al. (2008) [87]	Hangzhou Bay, China	17–6949	$TSS=\exp \left(43.233Rrs(B2)+1.396\right)$	0.76	RMSE = 424 mg/L	25
MOD-E38	Wu and Cui (2008) as cited in [22]	Poyang Lake, China	0-142	$\begin{array}{l}TSS=86236.23{(Rrs(B1))}^3-15858.70{(Rrs(B1))}^2\\ +1005.29(Rrs(B1))-15.67\end{array}$	0.92	-	42
MOD-E39	Kutser, et al. (2007) [26]	Muuga and Sillmae Port, Estonia	2–8	$TSS=349.83\rho (B1)+2.9663$	0.86	-	11
MOD-E40	Liu, et al. (2006) [58]	Middle Yangtze River, China	23.4–61.2	$\begin{array}{l}TSS=\exp (2.495x+1.81)\\ x=\frac{\rho (B1)-\rho (B2)}{\rho (B1)+\rho (B2)}\end{array}$	0.72	RE = 34.7%	41
MOD-E41	Sipelgas, et al. (2006) [27]	Parki Bay, Finland	3–10	$TSS=110.3\rho (B1)+2.0$	0.58	-	48
MOD-E42	Miller and Mckee, (2004) [3]	Northern Gulf of Maxico, USA	1.0–55.0	$TSS=1140.25\times R_{rs}(B1)-1.91$	0.89	RMSE = 4.74 mg/L	52
MOD-A1	Dorji, et al. (2016) [67]	Onslow, Western Australia	2.4–69.6	$\begin{array}{l}TSS=\frac{23.47\left(\frac{x(B1)}{1-x(B1)}\right)}{1-0.69\left(\frac{x(B1)}{1-x(B1)}\right)}|x(B1)=\frac{-g1+\sqrt{{\left(g1\right)}^2+4g2r_{rs}(B1)}}{2g2}\\ g1=0.084,g2=0.17,and{r}_{rs}(B1)=Rrs(B1)/(0.52+1.7Rrs(B1))\end{array}$	0.85	MARE = 33.33%	48
MOD-A2	Han, et al. (2016) [88]	Europe, French Guiana, Vietnam, North Canada, and China	0.154–2627	$\begin{array}{l}TSS=\frac{W_L\cdot TS{S}_L+W_H\cdot TS{S}_H}{W_L+W_H}\\ W_L=\{\begin{array}{l}1,ifRrs(B1)\le 0.03s{r}^{-1}\\ 0,ifRrs(B1)\ge 0.04s{r}^{-1}\\ {\log }_{10}(0.04)-{\log }_{10}(Rrs(B1)),otherwise\end{array}|TS{S}_L=\frac{404.4{\rho }_w(B1)}{1-{\rho }_w(B1)/0.5}\\ W_H=\{\begin{array}{l}0ifRrs(B1)\le 0.03s{r}^{-1}\\ 1,ifRrs(B1)\ge 0.04s{r}^{-1}\\ {\log }_{10}(Rrs(B1))-{\log }_{10}(0.03),otherwise\end{array}|TS{S}_H=\frac{1214.669{\rho }_w(B1)}{1-{\rho }_w(B1)/0.3394}\end{array}$	-	MRAD = 51.9-59%	TSSL = 366 TSSH = 46
MOD-A3	Shen, et al. (2014) [89]	Yangtze estuary, China	-	$\begin{array}{l}Rrs=\frac{\alpha \beta \times TSS}{1+\beta \times TSS+\sqrt{1+2\beta \times TSS}}\\ Rrs(B2)\}\alpha =0.1038,\beta =1.8042\end{array}$	0.91	RMSE = 0.0048 (sr−1)	144
MOD-A4	Vanhellemont and Ruddick (2014) [11]	Southern North Sea, UK	0.5–100	$\begin{array}{l}TSS=\frac{258.85\rho (B1)}{1-\rho (B1)/0.1641}\\ \rho (B1)=0.529\times \pi \times r_{rs}(B1)\end{array}$	-	-	-
MOD-A5	Chen, et al. (2013) [56]	Changjiang River Estuary, China	70–710	$TSS=1.7492{(Rrs{(B16)}^{-1}-Rrs{(B15)}^{-1})}^{-1}+0.0912$	0.89	MRE = 28.99%	20
MOD-A6	Katlane, et al. (2013) [90]	Gulf of Gabes	0.7–30	$\begin{array}{l}TSS=62.86\left(\raisebox{1ex}{$\rho (B1)$}\!\left/ \!\raisebox{-1ex}{$(0.1736-\rho (B1))$}\right.\right)\\ \rho (B1)=\pi Rrs(B1)\end{array}$	-	-	56
MOD-A7	Nechad, et al. (2010) [38]	Southern North Sea	1.24–110.27	$\begin{array}{l}TSS=\frac{400.75\rho (B1)}{1-\rho (B1)/0.1774}+1.02\\ \rho (B1)=0.539\times \pi r_{rs}(B1)\end{array}$	0.80	RMSE = 11.23 mg/L MRE = 38.9%	72
LAN-E1	Cai, et al. (2015) [91]	Hangzho Bay, China	203–481	$\ln (TSS)=0.01113115\times L(TM4)+4.794229$	0.951	-	35
LAN-E2	Cai, et al. (2015) [92]	Hangzho, Bay	179–389.58	$TSS=314.435\rho (TM3)+3805.982\rho (TM4)+28.54$	0.976	-	27
LAN-E3	Kong, et al. (2015) [7]	Gulf of Bohai Sea	2.1–208.7	$TSS=5.184+1349.63R_{rs}(TM4)+614561.673{\left[R_{rs}(TM4)\right]}^2$	0.844	RMSE = 5.59	70
LAN-E4	Kong, et al. (2015) [93]	Caofeidian, Bohai Sea	4.3–104.1	$\begin{array}{l}TSS=296.29x^2-272.62x+70.939\\ x=R_{rs}(TM3)/R_{rs}(TM2)\end{array}$	0.977	RMSE = 7.22 mg/L MRE = 25.35
LAN-E5	Lim and Choi (2015) [94]	Nakdong River, South Korea	~3–14	$\begin{array}{l}TSS=11.80-50.608\rho (OLI2)+14.58\rho (OLI5)\\ -4.764\rho (OLI5)/\rho (OLI3)\end{array}$	0.74	RMSE = 1.40	48
LAN-E6	Wu, et al. (2015) [9]	Dongting Lake, China	0–63.2	$TSS=1.1034\times \exp (23.61\times \rho (OLI4))$	0.91	RMSE = 4.41 mg/L	52
LAN-E7	Zheng, et al. (2015) [95]	Dongting Lake, China	4.0–101	$TSS=4616.4\times R_{rs}(TM4)-4.362$	0.82	MAPE = 21.3% RMSE = 7.01 mg/L	42
LAN-E8	In-Young, et al. (2014) [96]	Old Women Creek Estuary, Ohio, US	1.0–278	$\begin{array}{l}Turb=14.44\times \rho (\%,ETM4)+7.61\\ Turb=1.03TSS-3.54\end{array}$	0.65		11
LAN-E9	Zhang, et al. (2014) [10]	Yellow river estuary	1.0–1500	${\log }_{10}(TSS)=44.072\times R_{rs}(TM3+0.1591)$	0.9672	MRE = 26.1%	44
LAN-E10	Hao, et al. (2013) [97]	Yangtze Estuary, China	~40.0–750	$\begin{array}{c}R(B4)=0.0444+32.7736\times \left[TSS/\left(51.6753+TSS\right)\right]\\ -32.7260\times \left[TSS/\left(51.6753+TSS\right)\right]\times \exp (-0.0000012081\times TSS)\end{array}$	0.8175	ARE = 36.83	17
LAN-E11	Hicks, et al. (2013) [98]	Waikato River, New Zealand	2.0–962	$TSS=-52.817+1449.4\rho (ETM4)$	0.939	RMSE = 21.3	35
LAN-E12	Min, et al. (2013) [78]	Saemangeum coastal area, Korea	0.1–55	$TSS=0.24\exp (188.3R_{rs}(TM2))$	0.90	-	88
LAN-E13	Miller, et al. (2011) [99]	Albemarle-Pamlico Estuarine System, North Carolina, USA	~5.0–30	$TSS=1.7+684.76\rho (\%,ETM3)$	0.87	-	599
LAN-E14	Li, et al. (2010) [100]	Changjiang Estuary	~1.5–560	$\begin{array}{l}1o{g}_{10}(TSS)=0.892+6.2244x\\ x=\frac{\left[\rho (ETM2)+\rho (ETM3)\right]}{\left[\rho (ETM2)/\rho (ETM3)\right]}\end{array}$	0.915	-	21
LAN-E15	Wang, et al. (2009) [12]	Yangtze river, China	22–2610	$\ln (TSS)=3.18236\times \ln (\%,{\rho }_w(ETM4))$	0.88	MRE = 14.83%	24
LAN-E16	Onderka and Pekarova (2008) [101]	Danube River, Slovakia	19.5–57.5	$TSS=4.17\times L(ETM4)-43.22$	0.93	SE = 3.2 mg/L	10
LAN-E17	Teodoro, et al. (2008) [102]	Douro River and Mira Lagoon, Portugal	14–449	$TSS=15.483-12.688\rho (TM1)+44.495\rho (TM3)$	0.995	RMSE = 25.3 mg/L	11
LAN-E18	Aparslan, et al. (2007) [103]	Omerli Dam, Turkey	0.4–2.9	$\begin{array}{l}TSS=42.2671-0.8694\rho (ETM1)-0.3716\rho (ETM2)\\ +1.05\rho (ETM3)+0.1247\rho (ETM4)\end{array}$	0.99	SE = 0.0085 mg/L	6
LAN-E19	Wang, et al. (2007) [104]	Yangtze River, China	0–900	$TSS=71.392\times \rho (\%,ETM4)-272.48$	0.92	MAE = 68.9 RMSE = 83.2	14
LAN-E20	Doxaran, et al. (2006) [105]	Gironde Estuary, France	10–2000	$\begin{array}{l}TSS=29.022\times \exp (0.0335x)\\ x=R_{rs}(ETM4)/R_{rs}(ETM2)(\%)\end{array}$	0.88	SD = 21%	132
LAN-E21	Wang, et al. (2006) [33]	Lake Reelfoot, USA	11.5–33.5	$\begin{array}{l}TSS=22.74+156.202\rho (TM2)-147.62\rho (TM3)\\ -45.66\rho (TM4)\end{array}$	0.52	-	18
LAN-E22	Zhou, et al. (2006) [15]	Lake Taihu, China	48.32–120.80	$TSS=-122.4+304.1\left(\raisebox{1ex}{$\rho (TM4)$}\!\left/ \!\raisebox{-1ex}{$\rho (TM1)$}\right.\right)+10.37$	0.74	MPE = 65.40%
LAN-A1	Dorji, et al. (2016) [67]	Onslow, Western Australia	2.4–69.6	$\begin{array}{l}TSS=\frac{25.34\left(\frac{x(B1)}{1-x(B1)}\right)}{1-0.69\left(\frac{x(B1)}{1-x(B1)}\right)}|x(B1)=\frac{-g1+\sqrt{{\left(g1\right)}^2+4g2r_{rs}(B1)}}{2g2}\\ g1=0.084,g2=0.17,and{r}_{rs}(B1)=Rrs(B1)/(0.52+1.7Rrs(B1))\end{array}$	0.85	MARE = 33.36%	48
LAN-A2	Han, et al. (2016) [88]	Europe, French Guiana, Vietnam, North Canada, and China	0.154–2627	$\begin{array}{l}TSS=\frac{W_L\cdot TS{S}_L+W_H\cdot TS{S}_H}{W_L+W_H}\\ W_L=\{\begin{array}{l}1,ifRrs(B1)\le 0.03s{r}^{-1}\\ 0,ifRrs(B1)\ge 0.04s{r}^{-1}\\ {\log }_{10}(0.04)-{\log }_{10}(Rrs(B1)),otherwise\end{array}|TS{S}_L=\frac{346.353{\rho }_w(OLI4)}{1-{\rho }_w(OLI4)/0.5}\\ W_H=\{\begin{array}{l}0ifRrs(B1)\le 0.03s{r}^{-1}\\ 1,ifRrs(B1)\ge 0.04s{r}^{-1}\\ {\log }_{10}(Rrs(B1))-{\log }_{10}(0.03),otherwise\end{array}|TS{S}_H=\frac{1221.390{\rho }_w(OLI4)}{1-{\rho }_w(OLI4)/0.3329}\end{array}$	-	MRAD = 51.9%–59%	TSSL = 366 TSSH = 38
LAN-A3	Zhang, et al. (2016) [106]	Xinánjiang Resevoir, China	0.67–5.66	$\begin{array}{c}TSS=38.08\times [0.42\times ({(Rrs(OLI2)+Rrs(OLI3))}^{-1}\\ -{(Rrs(OLI3)+Rrs(OLI4))}^{-1})\times (Rrs(OLI4)+0.24]-3.39\end{array}$	>0.8	MRE = 24.3%	45
LAN-A4	Kong, et al. (2015) [7]	Gulf of Bohai Sea	2.1–208.7	$TSS=8.602+\text{}1805.26Rrs(TM4)+900713.14Rrs{(TM4)}^2$	0.844	RMSE = 4.53	70
LAN-A5	Vanhellemont and Ruddick (2014) [11]	Southern North Sea, UK	0.5–100	$\begin{array}{l}TSS=\frac{289.29{\rho }_w(OLI4)}{1-{\rho }_w(OLI4)/0.1686}\\ {\rho }_w(OLI4)=0.529\times \pi \times r_{rs}(OLI4)\end{array}$	-	-	-

MRE = Mean Relative Error, MARE = Mean Absolute Relative Error, MAE = Mean Absolute Error, MPE = Mean Percentage Error, SD = Standard Deviation, MRAD = Mean Relative Absolute Difference, RMSE = Root Mean Square Error, SE = Standard Error, ARE = Absolute Relative Error, RPD = Relative Percentage Difference, APD = Absolute Percentage Difference, RRMSE = Relative Root Mean Square Error.

. Semi-analytical models described in this section encompass all the semi-analytical models from MODIS (MOD-A) and Landsat (LAN-A). Models are considered semi-analytic because they are derived based on a physical form [56] or one or more parameters in the TSS algorithms are either parameterized using site-specific or global in-water bio-optical properties [38]. Semi-analytic algorithms for the two sensors considered here consist of algorithms that are based on radiative transfer modelling to relate the dependence of geo-physical properties of the water, TSS in our case, to the reflectance via IOPs of the water.

Empirical models consist of TSS algorithms that are directly related with in situ AOPs of water and the TSS using linear or non-linear regression methods. For the two optical sensors considered here the empirical algorithms from MODIS (MOD-E) and Landsat (LAN-E) will be collectively known as empirical algorithms unless otherwise stated explicitly. The form of the equations used in the empirical methods ranged from simple linear [17,26,57,58], exponential [9,21,34,59], power [10,46,60] and other polynomial relationships [61,62,63] using single, multiple or combinations of different bands in band ratio or self-formulated indexes. To differentiate the algorithms within each sensor, algorithms will be labeled with a respective number following each sensor’s name, MOD-A1 and MOD-E1 will represent MODIS semi-analytic algorithm 1 and MODIS empirical algorithm 1 respectively; likewise, a similar naming convention is followed for TSS algorithms for Landsat-based sensors.

2.3. Statistical Tests and Scoring System

The statistical tests used to evaluate the performance of each TSS algorithm for different types of water described in Section 2.1.2 are based on the statistical tests used by Brewin et al. [44]. Further, to objectively rank the TSS algorithms we used the point scoring system of Brewin et al. [44]. The details of each statistical test and scoring system of each test adopted from [44] are described in the following sections. Further, to contain the effect of spurious TSS generated by some of the TSS models being applied outside their range, we only included TSS estimations that were between a lower bound available in each TSS model (zero for the TSS models which did not contain the lower bound) and an upper bound of twice the highest TSS concentration reportedly used to calibrate each TSS model.

2.3.1. Pearson Correlation Coefficient (r) Test

The point scoring system for the r test involves determining if the r-value for each TSS algorithm is statistically significant when compared with the mean r-value for all TSS algorithms. The statistical significance is determined through z-scores and the z-score is computed through Fisher’s r-to-z transformation using relationships between the r-values of two models and the total number of samples used to determine the r-values, described in [44] as:

$$z_1=0.5\log \left(\frac{1+r_1}{1-r_1}\right)$$

(14)

$$z_2=0.5\log \left(\frac{1+r_2}{1-r_2}\right)$$

(15)

$$z_{score}=\frac{z_1-z_2}{\sqrt{\frac{1}{n_1-3}+\frac{1}{n_2-3}}}$$

(16)

where r₁ is the r-value of a specific TSS algorithm and r₂ is the mean of all r-values from all the TSS algorithms. Similarly, n₁ is the number of samples in a specific TSS algorithm and n₂ is the mean number of samples from all TSS algorithms. In the event that the TSS model fails to produce a reasonable estimate of TSS within the accepted bounds of each TSS model when tested for a particular water type then in such cases the value of n₁ can be different between two different water conditions, similarly, the value of n₂ also changes as it is the average number of samples of all TSS models in that particular water type.

For algorithm comparison, a two-tailed test was performed using the z-score to determine the p-value. If the p-value was less than 0.05 then the r-values were considered as statistically significant and for each TSS algorithm that were statistically significant the following scores were assigned comparing the r-value and the mean r-value $\left(\overline{r}\right)$ of all TSS algorithms:

$$\begin{array}{}\\ r-\mathrm{test}\\ \end{array}\}\begin{array}{c}0\mathrm{points}ifr<\overline{r}\\ 1\mathrm{point}ifr=\overline{r}\\ 2\mathrm{points}ifr>\overline{r}\end{array}$$

(17)

2.3.2. Root Mean Square Error (ψ) Test

The Root Mean Square Error (ψ) of a model estimate, y_i, with respect to a true value, x_i, can be computed using Equation (18):

$$\psi =\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(y_i-x_i\right)}^2}}$$

(18)

The 95% confidence intervals were also calculated for each TSS algorithm and the mean of all TSS algorithms. For each TSS algorithm, the following scoring points were assigned according to the conditions in Equation (19):

$$\begin{array}{}\\ \psi -test\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array}\begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}\psi -{\psi }_{95\%CI}>\overline{\psi }+{\overline{\psi }}_{95\%CI}\\ \overline{\psi }-{\overline{\psi }}_{95\%CI}\le \psi -{\psi }_{95\%CI}\le \overline{\psi }+{\overline{\psi }}_{95\%CI}\begin{array}{cc}\mathrm{or}& \overline{\psi }-{\overline{\psi }}_{95\%CI}\le \psi +{\psi }_{95\%CI}\le \overline{\psi }+{\overline{\psi }}_{95\%CI}\end{array}\\ \psi +{\psi }_{95\%CI}<\overline{\psi }-{\overline{\psi }}_{95\%CI}\end{array}$$

(19)

where ${\psi }_{95\%\mathrm{CI}}$ and ${\overline{\psi }}_{95\%\mathrm{CI}}$ is the 95% confidence interval of ψ and mean—ψ ($\overline{\psi }$) of all TSS algorithms respectively.

Figure 2 shows an example of scoring point classification for Landsat algorithms used in retrieving TSS concentration for the ψ—test.

2.3.3. The Bias (δ) Test

The bias (δ) of model estimate y_i and true x_i is calculated using Equation (20):

$$\delta =\frac{1}{N}{\displaystyle \sum _{i=1}^N}(y_i-x_i)$$

(20)

For each TSS algorithm, following score points were awarded according to the conditions in Equation (21):

$$\begin{array}{}\\ \delta -test\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array} \begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}{\delta }_{95\%\mathrm{CI}}>{\overline{\delta }}_{95\%\mathrm{CI}}\mathrm{AND}\delta -{\delta }_{95\%\mathrm{CI}}>0+{\delta }_{95\%\mathrm{CI}}\mathrm{or}\delta +{\delta }_{95\%\mathrm{CI}}>0-{\delta }_{95\%\mathrm{CI}}\\ {\delta }_{95\%\mathrm{CI}}<{\overline{\delta }}_{95\%\mathrm{CI}}\mathrm{or}0-{\overline{\delta }}_{95\%\mathrm{CI}}\le \delta +{\delta }_{95\%\mathrm{CI}}\le 0+{\overline{\delta }}_{95\%\mathrm{CI}}\mathrm{or}0-{\overline{\delta }}_{95\%\mathrm{CI}}\le \delta -{\delta }_{95\%\mathrm{CI}}\le 0+{\overline{\delta }}_{95\%\mathrm{CI}}\\ {\delta }_{95\%\mathrm{CI}}<{\overline{\delta }}_{95\%\mathrm{CI}}\mathrm{AND}0-{\overline{\delta }}_{95\%\mathrm{CI}}\le \delta -{\delta }_{95\%\mathrm{CI}}\le 0+{\overline{\delta }}_{95\%\mathrm{CI}}\mathrm{or}0-{\overline{\delta }}_{95\%\mathrm{CI}}\le \delta +{\delta }_{95\%\mathrm{CI}}\le 0+{\overline{\delta }}_{95\%\mathrm{CI}}\end{array}$$

(21)

where ${\delta }_{95\%\mathrm{CI}}$ and ${\overline{\delta }}_{95\%\mathrm{CI}}$ is the 95% confidence interval of mean—δ $(\overline{\delta })$ of all TSS algorithms respectively. Further, to score one point only one conditions must be satisfied while to score two points both the conditions must be satisfied.

2.3.4. The Center-Pattern Root Mean Square Error (Δ) Test

The center-pattern Root Mean Square Error (Δ) is calculated using Equation (21):

$$\Delta =\sqrt{\frac{1}{N}{{\displaystyle \sum _{i=1}^N\left\{\left(y_i-\overline{y}\right)-\left(x_i-\overline{x}\right)\right\}}}^2}$$

(22)

The 95% confidence intervals were also calculated for each TSS algorithm and the mean of all TSS algorithms. For each TSS algorithm, the following scores were assigned according to the conditions in Equation (23):

$$\begin{array}{}\\ \Delta -test\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array}\begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}\Delta -{\Delta }_{95\%\mathrm{CI}}>\overline{\Delta }+{\overline{\Delta }}_{95\%\mathrm{CI}}\\ \overline{\Delta }-{\overline{\Delta }}_{95\%\mathrm{CI}}\le \Delta -{\Delta }_{95\%\mathrm{CI}}\le \overline{\Delta }+{\overline{\Delta }}_{95\%\mathrm{CI}}\mathrm{or}\overline{\Delta }-{\overline{\Delta }}_{95\%\mathrm{CI}}\le \Delta +{\Delta }_{95\%\mathrm{CI}}\le \overline{\Delta }+{\overline{\Delta }}_{95\%\mathrm{CI}}\\ \Delta +{\Delta }_{95\%\mathrm{CI}}<\overline{\Delta }-{\overline{\Delta }}_{95\%\mathrm{CI}}\end{array}$$

(23)

where ${\Delta }_{95\%\mathrm{CI}}$ and ${\overline{\Delta }}_{95\%\mathrm{CI}}$ is the 95% confidence interval of Δ and mean—Δ ($\overline{\Delta }$) of all TSS algorithm respectively.

2.3.5. The Slope (S) and Intercept (I) of a Type-2 regression Test

The Slope (S) and Intercept (I) of a type-2 regression [64] were calculated using Equation (24):

$$Y=X\times S+I$$

(24)

where Y is the TSS estimates derived from the TSS algorithms and X the true TSS. The following scores were assigned by comparing the S—value of each TSS algorithm and mean—S ($\overline{s}$) value of all TSS algorithms.

$$\begin{array}{}\\ S-\mathrm{test}\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array}\begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}{\sigma }_s>{\sigma }_{\overline{s}}\mathrm{AND}S-{\sigma }_s>1+2{\sigma }_{\overline{s}}\mathrm{or}S+{\sigma }_s<1-2{\sigma }_{\overline{s}}\\ {\sigma }_s<{\sigma }_{\overline{s}}\mathrm{or}1-2{\sigma }_{\overline{s}}\le S-{\sigma }_s\le 1+2{\sigma }_{\overline{s}}\mathrm{or}1-2{\sigma }_{\overline{s}}\le S+{\sigma }_s\le 1+2{\sigma }_{\overline{s}}\\ {\sigma }_s<{\sigma }_{\overline{s}}\mathrm{AND}1-2{\sigma }_{\overline{s}}\le S-{\sigma }_s\le 1+2{\sigma }_{\overline{s}}\mathrm{or}1-2{\sigma }_{\overline{s}}\le S+{\sigma }_s\le 1+2{\sigma }_{\overline{s}}\end{array}$$

(25)

where ${\sigma }_{\overline{s}}$ is the standard deviation of $\overline{s}$ from all TSS algorithms.

For the I parameter, for each TSS algorithm, the following scores were assigned according to the conditions in Equation (26).

$$\begin{array}{}\\ I-\mathrm{test}\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array}\begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}{\sigma }_I>{\sigma }_{\overline{I}}\mathrm{AND}I-{\sigma }_I>0+2{\sigma }_{\overline{I}}\mathrm{or}I+{\sigma }_I<0-2{\sigma }_{\overline{I}}\\ {\sigma }_I<{\sigma }_{\overline{I}}\mathrm{or}0-2{\sigma }_{\overline{I}}\le I-{\sigma }_I\le 0+2{\sigma }_{\overline{I}}\mathrm{or}0-2{\sigma }_{\overline{I}}\le I+{\sigma }_I\le 0+2{\sigma }_{\overline{I}}\\ {\sigma }_I<{\sigma }_{\overline{I}}\mathrm{AND}0-2{\sigma }_{\overline{I}}\le I-{\sigma }_I\le 0+2{\sigma }_{\overline{I}}\mathrm{or}0-2{\sigma }_{\overline{I}}\le I+{\sigma }_I\le 0+2{\sigma }_{\overline{I}}\end{array}$$

(26)

where ${\sigma }_{\overline{I}}$ is the standard deviation of mean—I ($\overline{I}$) from all TSS algorithms. Further, in the S and I—test in Equations (25) and (26), to score one point only one of the two conditions must be satisfied while to score two points both the conditions must be satisfied.

2.3.6. Percentage of Possible Retrievals (η)

The percentage of possible retrievals (η) was calculated using Equation (27):

$$\eta =\frac{N^E}{N^M}\times 100\%$$

(27)

where N^E is the total number of TSS retrieved using each TSS algorithm from the total number of TSS concentrations (N^M) considered in the study. For the point scoring system the following basis was followed:

$$\begin{array}{}\\ \eta -test\\ \end{array}\}\begin{array}{c}0\mathrm{points}\\ 1\mathrm{point}\\ 2\mathrm{points}\end{array}\begin{array}{c}if\\ if\\ if\end{array}\begin{array}{c}\eta <\overline{\eta }-{\sigma }_{\overline{\eta }}\\ \overline{\eta }-{\sigma }_{\overline{\eta }}\le \eta \le \overline{\eta }+{\sigma }_{\overline{\eta }}\\ \eta >\overline{\eta }+{\sigma }_{\overline{\eta }}\end{array}$$

(28)

where $\overline{\eta }$ and ${\sigma }_{\overline{\eta }}$ is the mean η-value and its standard deviation for all TSS algorithm in η-test.

2.3.7. Total Points

For objective comparison the performance of each TSS algorithm with respect to different water types, all points from each statistical test were added and normalized by the mean score of all TSS algorithms. Thus, a score of zero indicates that the TSS algorithm is performing lower than the mean of all TSS algorithms, a score of one indicates that the TSS algorithm is at par with the average of all TSS algorithms, and a score of 2 means the TSS algorithm is better than the mean of all TSS algorithms.

Further, to test the stability of the scoring systems used in this study, we used a bootstrap method [65], with 1000 runs and each time generating a new dataset by resampling via the replacement method for model generated TSS and HydroLight TSS. Each new dataset was tested using the statistical tests and scored using the scoring system described in Section 2.3. The results reported in the Section 3 are the mean values of the total points from the bootstrapping method with 2.5 and 97.5 percentiles reported as uncertainty estimates for 95% confidence limits.

2.3.8. Mean of Total Points

The mean of total points is achieved by averaging the score of each TSS model across different sediment types or solar zenith angles or backscattering ratios for the particular water classes described in Section 2.3.7. For example, in the case of different sediment types in CLASS-I water, the mean of total points in CLASS-I is an aggregate of total scores of each TSS model for different sediment types. For the case of different backscattering ratios and solar zenith angle, the mean of total points is an aggregate of each TSS model for different backscattering ratios and solar zenith angles, respectively, for a specific sediment type in a particular water class. The error bars in the mean of total points are the mean of uncertainty estimates of the total points obtained from the 95% confidence limit from the bootstrapping method.

2.3.9. Final Score

The final score is the aggregate of the mean of total points across all water classes for different sediment types, backscattering ratios and solar zenith angles. For example, the final score for MOD-E1 is derived as the mean from the aggregate score of MOD-E1 at five different sediment types, backscattering ratios, and solar zenith angles across all five different water classes. The error bars are the standard deviation of errors from the mean of total points across all five different water classes. Figure 3 shows an illustration of the point score system adapted from [44] and used in comparing TSS models in this study. The error bars in the Final score are the mean of uncertainty estimates from the mean of total points.

3. Results

3.1. TSS Model Comparisons

Figure 4 and Figure 5 show the quantitative comparison between the models using the final scores which are aggregates of the total scores from different sediment types, backscattering ratios, and solar zenith angles across all five different water classes for MODIS and Landsat-based models respectively. The final results presented in Figure 4 and Figure 5 are indications of the overall performance of the TSS models when weighted across different water types, sediment types and backscattering ratios. The detailed results of individual model performance in respective sediment types, backscattering ratios, and water types are presented in Supplementary Material, S11. In addition, the Supplementary Materials S1–S10 also provide the detailed statistical test results for each TSS model.

From the final scores displayed in Figure 4 and Figure 5 we can visually observe that there are clearly high and low performing models. The high performing MODIS TSS models with final scores greater than 1.5, in the order of highest to lowest final score, are MOD-E6, MOD-A1, MOD-E28, MOD-A4, MOD-E10, and MOD-E42 and low performing MODIS TSS models with scores less than 0.5 are MOD-E8, MOD-E2, MOD-E24, MOD-E22 and MOD-E32. For the Landsat TSS models, LAN-E3, LAN-A4, LAN-E9, LAN-A5, and LAN-A1 have final scores greater than 1.5 while LAN-E11, LAN-E22, LAN-E16, and LAN-E18 have final scores less than 0.5. In the final scores of low performing TSS models, the LAN-E18 model has scores of zeros which shows that LAN-E18 failed to derive TSS within the acceptable TSS bounds of 0.4–5.8 mg/L. We suspect the published algorithm includes an error. The overall ranking of the TSS models using the final scores for each TSS model is also presented in

Table B1. Mean of Total Point and Final Scores of MODIS TSS models across different water classes as derived from different sediment types, backscattering ratios and Solar Zenith Angles. The top five and bottom five scores from each water types and the final scores are in bold (top) and bold italics (bottom).

MODEL	Mean Total Score from Sediment					Mean Total Score from Backscattering Ratio (bb/b)					Mean Total Score from Solar Zenith Angles					Final Score	Error Final Score
	I	II	III	IV	V	I	II	III	IV	V	I	II	III	IV	V			Lower Bound	Upper Bound
MOD-E6	1.69	1.61	1.66	1.61	1.63	2.00	1.72	1.98	1.75	1.72	1.71	1.60	1.67	1.53	1.59		1.70	1.64	1.76
MOD-A1	1.46	1.53	1.50	1.56	1.46	1.54	1.71	1.57	1.82	1.67	1.54	1.67	1.55	1.73	1.65	1.60	1.55	1.63
MOD-E28	1.53	1.51	1.53	1.51	1.51	1.71	1.71	1.71	1.74	1.71	1.52	1.55	1.52	1.51	1.56	1.59	1.56	1.63
MOD-A4	1.47	1.55	1.48	1.42	1.43	1.57	1.71	1.57	1.59	1.51	1.61	1.62	1.60	1.54	1.57	1.55	1.51	1.60
MOD-E10	1.54	1.54	1.54	1.54	1.54	1.57	1.57	1.57	1.57	1.57	1.59	1.50	1.55	1.47	1.50	1.54	1.54	1.56
MOD-E42	1.48	1.49	1.46	1.42	1.47	1.57	1.16	1.57	1.76	1.57	1.61	1.17	1.60	1.51	1.62	1.50	1.40	1.63
MOD-E21	1.57	1.50	1.58	1.49	1.50	1.73	1.46	1.76	1.51	1.53	1.68	1.24	1.37	1.20	1.29	1.49	1.40	1.60
MOD-E31	1.45	1.46	1.43	1.42	1.42	1.55	1.60	1.52	1.46	1.46	1.55	1.51	1.51	1.48	1.55	1.49	1.38	1.58
MOD-A6	1.47	1.46	1.49	1.42	1.40	1.43	1.57	1.43	1.57	1.43	1.47	1.54	1.49	1.50	1.47	1.48	1.44	1.53
MOD-A7	1.50	1.47	1.54	1.47	1.44	1.44	1.53	1.57	1.55	1.43	1.53	1.31	1.59	1.28	1.25	1.46	1.39	1.51
MOD-E44	1.32	1.30	1.31	1.26	1.30	1.57	1.56	1.57	1.51	1.55	1.58	1.47	1.57	1.44	1.56	1.46	1.39	1.49
MOD-E27	1.38	1.42	1.37	1.41	1.41	1.46	1.57	1.47	1.57	1.54	1.49	1.33	1.49	1.27	1.35	1.44	1.38	1.50
MOD-E4	1.47	1.41	1.47	1.40	1.42	1.57	1.43	1.57	1.43	1.45	1.49	1.36	1.47	1.32	1.39	1.44	1.41	1.49
MOD-E34	1.43	1.43	1.43	1.43	1.43	1.43	1.43	1.43	1.45	1.43	1.49	1.43	1.50	1.44	1.43	1.44	1.43	1.46
MOD-E41	1.41	1.40	1.41	1.40	1.41	1.43	1.43	1.43	1.43	1.43	1.46	1.34	1.45	1.32	1.46	1.41	1.41	1.43
MOD-E20	1.36	1.40	1.33	1.46	1.46	1.33	1.47	1.30	1.57	1.52	1.37	1.30	1.33	1.43	1.53	1.41	1.34	1.49
MOD-E35	1.15	1.52	1.15	1.56	1.23	1.29	1.58	1.29	1.68	1.39	1.29	1.53	1.28	1.56	1.40	1.39	1.33	1.51
MOD-E39	1.31	1.31	1.31	1.31	1.31	1.28	1.29	1.29	1.29	1.29	1.31	1.26	1.30	1.25	1.31	1.29	1.29	1.30
MOD-E25	1.15	1.19	1.14	1.32	1.20	1.31	1.40	1.24	1.15	1.34	1.39	1.36	1.31	1.32	1.39	1.28	1.10	1.49
MOD-E3	0.99	1.21	0.83	1.25	1.10	1.39	1.75	1.09	1.67	1.53	1.29	1.33	1.03	1.23	1.53	1.28	1.09	1.48
MOD-E19	1.39	1.22	1.42	1.26	1.24	1.30	1.12	1.40	1.36	1.14	1.38	0.90	1.43	0.98	1.22	1.25	1.05	1.36
MOD-E40	1.14	1.20	1.14	1.23	1.20	1.14	1.29	1.14	1.29	1.29	1.16	1.25	1.15	1.29	1.24	1.21	1.20	1.22
MOD-E11	1.15	1.19	1.12	1.21	1.18	1.23	1.26	1.16	1.18	1.26	1.26	1.19	1.15	1.18	1.25	1.20	1.11	1.28
MOD-E37	1.13	1.09	1.13	1.09	1.10	1.24	1.22	1.25	1.27	1.23	1.14	1.11	1.14	1.11	1.14	1.16	1.08	1.23
MOD-E36	1.18	1.17	1.19	1.17	1.16	1.14	1.14	1.14	1.14	1.14	1.16	1.11	1.15	1.10	1.11	1.15	1.12	1.17
MOD-A5	1.14	1.12	1.14	1.12	1.13	1.14	1.14	1.14	1.14	1.14	1.14	1.14	1.13	1.14	1.14	1.14	1.13	1.14
MOD-E5	1.30	1.16	1.32	1.18	1.19	1.19	0.90	1.20	1.08	0.96	1.33	0.94	1.22	1.01	0.97	1.13	1.02	1.26
MOD-E30	1.12	1.09	1.12	1.08	1.09	1.14	1.14	1.14	1.14	1.14	1.16	1.14	1.16	1.13	1.14	1.13	1.10	1.16
MOD-E17	1.22	1.06	1.24	1.06	1.07	1.27	1.00	1.29	1.00	1.00	1.30	0.81	1.31	0.75	1.10	1.10	1.07	1.12
MOD-E18	1.23	1.12	1.08	1.05	1.11	1.32	1.16	1.12	0.84	0.92	1.47	1.01	1.28	0.71	0.88	1.09	0.88	1.30
MOD-E7	1.17	1.13	1.19	1.15	1.15	0.98	1.04	1.01	1.13	1.02	0.99	1.05	0.99	1.09	1.05	1.08	0.93	1.19
MOD-E14	1.18	1.03	1.21	1.03	1.03	1.16	1.00	1.22	1.00	1.00	1.20	0.80	1.24	0.75	1.07	1.06	1.02	1.11
MOD-E13	0.85	1.02	0.82	1.08	1.03	0.84	1.25	0.81	1.41	1.28	0.97	1.11	0.93	1.17	1.16	1.05	0.87	1.22
MOD-E16	1.11	1.04	1.06	1.16	1.12	1.00	1.00	1.00	1.14	1.05	1.04	0.80	1.05	0.85	1.15	1.04	1.01	1.07
MOD-E12	1.00	1.14	1.04	1.14	1.13	0.88	1.00	0.89	1.07	1.09	0.99	0.97	0.99	1.01	0.96	1.02	0.88	1.16
MOD-E33	0.94	1.03	0.91	1.07	1.03	0.90	1.10	0.82	1.14	1.12	0.93	1.02	0.88	1.03	1.07	1.00	0.90	1.09
MOD-E29	1.00	0.91	1.03	0.87	0.92	1.14	1.00	1.14	1.00	1.00	1.10	0.81	1.10	0.78	0.94	0.98	0.92	1.03
MOD-E45	0.94	1.09	0.94	1.09	0.93	0.87	1.13	0.87	1.22	0.93	0.82	0.98	0.80	0.98	0.86	0.96	0.88	1.08
MOD-E1	0.91	0.85	0.92	0.84	0.86	0.76	0.72	0.78	0.72	0.72	0.82	0.82	0.85	0.83	0.86	0.82	0.77	0.93
MOD-E26	0.85	0.79	0.86	0.83	0.78	0.62	0.50	0.64	0.72	0.55	0.80	0.65	0.82	0.80	0.76	0.73	0.55	0.92
MOD-E15	0.45	0.86	0.44	0.85	0.58	0.35	0.96	0.35	0.98	0.67	0.49	0.85	0.49	0.82	0.77	0.66	0.52	0.88
MOD-E9	0.60	0.71	0.59	0.75	0.73	0.48	0.49	0.47	0.58	0.56	0.68	0.57	0.68	0.67	0.60	0.61	0.49	0.80
MOD-E38	0.45	0.52	0.42	0.61	0.56	0.50	0.64	0.50	0.75	0.63	0.62	0.62	0.58	0.73	0.66	0.59	0.47	0.89
MOD-E23	0.75	0.44	0.79	0.44	0.44	0.75	0.30	0.80	0.53	0.30	0.80	0.27	0.86	0.34	0.57	0.56	0.43	0.69
MOD-E8	0.51	0.41	0.51	0.35	0.48	0.62	0.60	0.59	0.24	0.55	0.60	0.44	0.56	0.24	0.40	0.47	0.18	0.67
MOD-E2	0.51	0.42	0.51	0.34	0.48	0.63	0.59	0.60	0.23	0.54	0.61	0.45	0.56	0.24	0.40	0.47	0.17	0.67
MOD-E24	0.46	0.45	0.45	0.48	0.46	0.44	0.44	0.43	0.55	0.47	0.45	0.49	0.44	0.57	0.51	0.47	0.43	0.58
MOD-E22	0.44	0.23	0.54	0.32	0.24	0.31	0.11	0.42	0.37	0.16	0.46	0.09	0.58	0.25	0.40	0.33	0.17	0.52
MOD-E32	0.38	0.31	0.49	0.42	0.32	0.04	0.14	0.36	0.45	0.18	0.29	0.15	0.53	0.40	0.39	0.32	0.18	0.57

and

Table B2. Mean of Total Point and Final Scores of Landsat TSS models across different water classes as derived from different sediment types, backscattering ratios and Solar Zenith Angles. The top five and bottom five scores from each water types and the final scores are in bold (top) and bold italics (bottom).

MODEL	Mean Total Score from Sediment					Mean Total Score from Backscattering Ratio (bb/b)					Mean Total Score from Solar Zenith Angles					Final Score	Error
																	Final Score
	I	II	III	IV	V	I	II	III	IV	V	I	II	III	IV	V		Lower Bound	Upper Bound
LAN-E3	1.66	1.69	1.67	1.69	1.70	1.85	1.81	1.86	1.82	1.84	1.79	1.61	1.74	1.61	1.64	1.73	1.60	1.77
LAN-A4	1.54	1.63	1.54	1.63	1.64	1.62	1.69	1.64	1.71	1.67	1.59	1.50	1.56	1.48	1.50	1.60	1.46	1.69
LAN-E9	1.28	1.38	1.24	1.36	1.39	1.74	1.94	1.65	1.97	1.98	1.47	1.54	1.39	1.55	1.62	1.57	1.47	1.64
LAN-A5	1.38	1.51	1.39	1.52	1.43	1.52	1.60	1.53	1.59	1.49	1.52	1.54	1.52	1.52	1.52	1.51	1.44	1.59
LAN-A1	1.33	1.53	1.34	1.58	1.45	1.54	1.69	1.46	1.80	1.63	1.52	1.27	1.48	1.33	1.59	1.50	1.43	1.56
LAN-E14	1.47	1.32	1.46	1.33	1.45	1.76	1.37	1.78	1.42	1.56	1.75	1.09	1.74	1.12	1.56	1.48	1.35	1.60
LAN-E20	1.56	1.57	1.53	1.61	1.60	1.52	1.45	1.52	1.48	1.48	1.54	1.17	1.51	1.16	1.49	1.48	1.37	1.59
LAN-E4	1.53	1.42	1.42	1.33	1.58	1.41	1.54	1.46	0.91	1.45	1.51	1.21	1.46	0.68	1.50	1.36	1.23	1.48
LAN-E1	1.36	1.36	1.36	1.36	1.36	1.37	1.34	1.36	1.35	1.35	1.29	1.33	1.29	1.32	1.30	1.34	1.28	1.40
LAN-E8	1.31	1.35	1.32	1.35	1.35	1.35	1.36	1.36	1.41	1.36	1.29	1.26	1.27	1.28	1.28	1.33	1.18	1.42
LAN-E13	1.36	1.39	1.38	1.35	1.37	1.30	1.28	1.35	1.27	1.30	1.35	1.20	1.35	1.12	1.35	1.31	1.28	1.37
LAN-E2	1.33	1.33	1.32	1.34	1.33	1.34	1.30	1.33	1.32	1.30	1.16	1.26	1.21	1.26	1.23	1.29	1.27	1.35
LAN-A2	1.18	1.08	1.19	1.11	1.12	1.38	1.04	1.43	1.20	1.23	1.42	1.08	1.41	1.16	1.26	1.22	1.13	1.38
LAN-E21	1.16	1.11	1.20	1.10	1.11	1.28	1.13	1.43	1.01	1.08	1.25	1.16	1.40	1.05	1.15	1.17	1.07	1.24
LAN-E7	1.11	0.93	1.09	0.93	0.85	1.38	1.04	1.39	1.09	0.89	1.47	0.79	1.46	0.74	0.91	1.07	0.83	1.31
LAN-E17	1.09	1.04	1.10	1.08	1.09	1.00	0.99	1.01	1.02	1.00	0.89	0.97	0.89	0.97	0.94	1.01	0.98	1.03
LAN-E12	1.13	1.02	0.96	1.12	1.24	1.09	0.74	1.11	0.61	1.05	1.25	0.71	1.19	0.50	0.75	0.96	0.73	1.20
LAN-E15	0.98	0.91	0.97	0.99	0.97	1.04	0.92	1.02	0.97	0.99	1.06	0.71	1.09	0.69	0.99	0.95	0.83	1.04
LAN-E5	0.97	0.95	0.95	0.99	0.99	0.94	0.88	0.90	0.89	0.96	0.97	0.70	0.94	0.72	1.03	0.92	0.76	1.05
LAN-A3	0.93	0.93	0.90	0.93	0.89	0.92	0.85	0.88	0.66	0.61	0.92	0.91	0.90	0.72	0.69	0.84	0.68	1.02
LAN-E19	0.66	0.67	0.67	0.69	0.64	0.64	0.73	0.80	0.87	0.76	0.60	0.65	0.66	0.70	0.67	0.69	0.45	1.07
LAN-E6	0.59	0.68	0.57	0.73	0.66	0.61	0.61	0.56	0.76	0.63	0.68	0.58	0.62	0.69	0.65	0.64	0.53	0.81
LAN-E10	0.42	0.45	0.39	0.45	0.45	0.65	0.59	0.61	0.65	0.48	0.66	0.44	0.60	0.44	0.36	0.51	0.28	0.78
LAN-E11	0.40	0.46	0.40	0.48	0.41	0.45	0.37	0.46	0.52	0.38	0.42	0.30	0.40	0.36	0.27	0.41	0.23	0.67
LAN-E22	0.99	0.84	1.02	0.67	0.75	0.56	0.00	0.19	0.00	0.00	0.47	0.05	0.34	0.00	0.00	0.39	0.31	0.51
LAN-E16	0.29	0.20	0.30	0.30	0.31	0.43	0.30	0.43	0.42	0.45	0.36	0.27	0.41	0.32	0.44	0.35	0.16	0.62
LAN-E18	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00

for MODIS and Landsat respectively. Further, Table B1 and Table B2 provides mean total scores for different sediment types, backscattering ratios, and solar zenith angles in all five water classes for respective TSS model.

With respect to the results displayed in Figure 4 and Figure 5, without the inclusion of error bars the distinction between the high performing TSS models is clear and we can easily compare the scores of each TSS model to obtain a ranking. For instance, in Figure 4 and Figure 5, the MOD-E6 and LAN-E3 are the highest scoring models with final scores of 1.70 and 1.73 respectively. However, on inclusion of the error bars, all high performing TSS models may be considered comparable and difficult to separate in terms of robustness, thus may all be ranked equally. Likewise, the case is similar for low performing TSS models where their error bars overlap. Further, we observe that two and three of the top five high scoring TSS models in MODIS and Landsat respectively are semi-analytic while none of the semi-analytic models were in the bottom five low scoring models.

3.2. Evaluation of Models

3.2.1. Model Evaluation Using HydroLight Data

The five high and low scoring models from MODIS and Landsat TSS models were selected to further evaluate their performance. From all available HydroLight data discussed in Section 2.1.2, the aforementioned high scoring TSS models were evaluated for their Relative Error (RE) between model-derived and HydroLight TSS. From the results presented in

Table 3. Relative Error and ΔR_rs Uncertainty Tolerance results for the highest and lowest scoring models’ evaluation using HydroLight Data. The highest scoring models are in bold text and the lowest scoring models are in regular italic text. The results provided in parenthesis represent the +ΔR_rs and ‘-’ indicates the model failed to provide TSS estimation within acceptable bounds. SRE: Smallest Relative Error. LRE: Largest Relative Error. MARE: Mean Absolute Relative Error. ARE: Absolute Relative Error.

Model	Relative Errors from HydroLight Data Validation			ARE from RRS Uncertainty (%)
	SRE (%)	MARE (%)	LRE (%)	−(+) 10% Δ Rrs	−(+) 20% Δ Rrs	−(+) 50% Δ Rrs
MOD-E6	59.35	94.30	139.35	70.46 (113.02)	44.59 (129.11)	91.94 (170.65)
MOD-A1	15.00	75.56	151.14	39.24 (126.59)	38.89 (182.84)	97.92 (294.93)
MOD-E28	51.61	148.62	191.97	97.96 (211.76)	49. 89 (271.68)	53.30 (497.21)
MOD-A4	63.14	257.59	386.87	157.51 (346.27)	68.10 (410.35)	96.13 (530.23)
MOD-E10	32.17	92.42	171.47	53.64 (149.97)	33.54 (242.01)	49.85 (396.29)
MOD-E8	189.55	220.69	344.16	244.77 (197.29)	268.89 (180.18)	341.16 (164.68)
MOD-E2	189.55	220.69	344.16	244.77(197.29)	268.89 180.18()	341.16 (164.68)
MOD-E24	77.87	141.49	218.80	10824.61 (9960.40)	11278.06 (9549.92)	12747.84 (8416.88)
MOD-E22	42.31	1832.79	5403.47	2461.87(1149.55)	1369.44 (1306.50)	187.31 (1206.94)
MOD-E32	39.90	1717.85	6778.93	2575.05(1067.58)	1381.65 (1385.73)	184.20 (288.28)
LAN-E3	59.31	120.37	166.68	69.03 (170.14)	33.14 (220.15)	76.58 (387.62)
LAN-A4	57.05	197.26	266.40	134.36 (262.03)	72.73 (331.63)	74.29 (541.89)
LAN-E9	23.52	481.82	1109.80	171.42 (857.00)	51.00 (1167.00)	92.43 (1974.47)
LAN-A5	62.86	244.28	362.44	149. 20 (341.63)	66. 53 (414.85)	95.90 (543.85)
LAN-A1	16.07	69.96	141.53	38.02 (115.85)	39.00 (169.17)	97. 78 (286.31)
LAN-E10	76.17	106.43	118.16	88.74 (126.91)	82.69 (161.62)	−(357.92)
LAN-E11	213.54	241.28	337.58	260.07 (22.48)	278.86 (203.89)	335.21 (177.52)
LAN-E22	19.41	110.69	164.56	110. 70 (110.688)	110. 64 (110.72)	196.66 (110.60)
LAN-E16	77.55	135.45	222.93	150.00 (109.18)	151.20 (103.59)	223.24 (85.67)
LAN-E18	-	-	-	-	-	-

we observe that there is high variability in the RE results amongst the respective MODIS and Landsat TSS models. The differences in the Smallest Relative Error (SRE) for high scoring TSS models were not as large as the differences within the MARE and Largest Relative Error (LRE). The MARE ranged from a low of 69.96% to a high of 481.82% while the SRE and LRE ranged from 15% to 63.14% and 139.35% to 1109.80% respectively. In the low scoring models, the high variability in the RE was observed with the MARE for low performing models ranging from 106.43% to 1832.79% while the SRE and LRE ranged from 39.90% to 213.54% and 118.16% to 6778.93% respectively. In both MODIS and Landsat high scoring models, the LRE results were for backscattering ratios of 0.001 and for Bukata type sediment. The SRE results were for backscattering ratios of 0.01 and calcareous sand sediment. Further, for the SRE in both the high and low performing TSS models, we observe that the high and low performing TSS models scored reasonably well in either one of the categories in sediment types, backscattering ratios, solar zenith angle and water classes. For instance, the low performing LAN-E22 scored higher than most of the high scoring TSS models in SRE results which indicated that LAN-E22 retrieves better in one of the water types.

The TSS derived using real satellite-data are bound by uncertainty related to observational, instrumental, measurement and data processing errors, the latter largely associated with the atmospheric correction procedure [66]. Thus, to assess the tolerance of high and low performing TSS models to the uncertainties in R_rs, which is the key input in derivation of the TSS concentration, we simulated the effect of R_rs uncertainty (󠄢ΔR_rs) by varying the R_rs by ±10%, ±20% and ±50% of the HydroLight generated R_rs at each of the MODIS and Landsat bands. The R_rs ± ΔR_rs was used in deriving TSS concentration and compared with HydroLight input TSS to calculate the Absolute Relative Error (ARE) of the TSS model. Table 3 reports the ARE and the MARE of HydroLight Data Validation as defined in Equation C2 in Appendix C. In general, we observe that with the increase in ΔR_rs the ARE also increases and the errors are higher for +ΔR_rs than −ΔR_rs. The ARE for high scoring TSS models ranged from 33.14% to 1974.47% while for low scoring TSS models it ranged from 82.69% to 12747.84% which shows both high and low performing TSS models are not impervious to uncertainty in R_rs measurements. However, high scoring TSS models show better tolerance to ΔR_rs than the low scoring TSS models. The details of the TSS models deviation in estimating TSS concentration from the error-free HydroLight data with ΔR_rs are shown in Table 3.

3.3.2. Model Evaluation Using In situ Data

As part of a regional water monitoring program, in situ reflectance and TSS measurements were carried out for the waters off the coast of northern Western Australia to develop regional TSS models (see MOD-A1 and LAN-A1 in Appendix A) [67]. The details of the in situ measurements and regional TSS model developed using in situ data can be obtained from [67]. A set of high scoring models (MOD-E10, MOD-A4, LAN-E9, and LAN-A5) and low scoring models (MOD-E1, MOD-E38, LAN-E6, and LAN-A3) were selected to compare with MOD-A1 and LAN-A1 in the context of in situ data comparisons. These subsets of models were selected because the reflectance bands used by other high scoring models were beyond the available reflectance bands in the in situ data.

Table 4. The MARE for high and low scoring models for in situ data. The high scoring models are in bold text and the low scoring models are in italics.

Error/Model	MOD-E10	MOD-A1 *	MOD-A4	MOD-E1	MOD-E38	LAN-E9	LAN-A1 *	LAN-A5	LAN-E6	LAN-A3
Mare (%)	46.20	33.33	100.85	341.04	256.00	43.11	33.36	102.59	55.23	35.62

* MARE was obtained using the leave-one-out cross validation method discussed in [67].

shows the Mean RE results obtained from each of the model evaluations against in situ data. Table 4 displays a high variability in the Mean RE for model comparisons for high scoring models with in situ data, from a low of 43.11% for LAN-E9 to a high of 102.59% for LAN-A5. When compared with the regional model’s MOD-A1 and LAN-A1 MARE results, we see that both the high scoring TSS models MOD-E10 and LAN-E9 and low scoring TSS models LAN-E6 and LAN-A3 were comparable. However, the results presented in Table 4 also show the extreme variability observed in the Mean RE for the low scoring models with a low of 35.62% and a high of 256%.

4. Discussion

4.1. Data and Methodological Limitation

The data used in this study to quantitatively compare TSS models have been generated using the widely used [68,69] in-water radiative transfer model HydroLight 4.2. The simulated data do not encompass all different water types in which each TSS model was developed to be used, however, it does provide us with a dataset that is independent of the data that has been used to parameterize the models to avoid biases in the results. To include all the models in comparisons, the simulated data were extrapolated to the NIR region of the spectrum using the methods discussed in Section 2.1.2. The extrapolation of reflectance data can introduce unrealistic values if the underlying assumptions of the spline extrapolation methodology does not hold true for the NIR regions. The extrapolation of the data is not ideal when used in modelling remote sensing products but the error for extrapolation had a MARE of 4.0% which was considered to be acceptable for this study. The ideal case for data for model comparisons would be to use a real global water data base, which is currently not available. The NOMAD dataset (http://seabass.gsfc.nasa.gov/) that is currently the most extensive dataset of in situ reflectance measurement and in-water variables did not contain the TSS measurements essential for this study.

The use of the objective methodology [44] of comparing models, used in this study to compare TSS models, can aid users in selection of TSS models that are best suited for waters of regional interest in the absence of means and a method to produce their own regionally tuned TSS algorithms. However, the objective methodology used here is not without limitations, as discussed by [44] with respect to using average performance to classify between high and low performing models. The very low performance of one particular model would affect the average of all other models to the extent that it becomes difficult to differentiate scores between models. For example, in Figure S11.1 for the score of MODIS TSS models in yellow clay, MOD-E1-2, E8-9, E15, E22-24, E32, and E38 all have low scores which increases the score of other TSS models making it difficult to differentiate among high scoring models. This problem is further exacerbated when the majority of TSS models score low which makes the few remaining high scoring models to appear similar in score, which is the case in Figure S11.22 for b_b/b of 0.001.

The objective classification was conducted on a case by case basis for different water types, sediment types, solar zenith angles, and backscattering ratios. The overall low performance of models in the final scores in Figure 4 and Figure 5 does not necessarily mean that low performing TSS models scored less in all the categories used in deriving the final score. For example, in Figure 5, LAN-E22 scored a very low final score when compared to other TSS models, but when considering specific results as presented in Figures S11.16–S11.20, LAN-E22 received a score at least comparable with most of the best performing models in all water classes for the red clay sediment type. Likewise, similar cases can be ascertained for all the respective TSS model’s scores for specific water classes, sediment types, backscattering ratios and solar zenith angle (Results provided as Supplementary Material S11 for other overall low scoring TSS models in Figure 4 and Figure 5). An additional disadvantage of the objective methodology used here is that the final score does not necessarily show the performance of all models in different categories considered, it shows only the relative performance of models in comparison with the mean scores of TSS models. In Figure S11.1, we observe that almost all TSS models score relatively higher total points for brown earth and lower for Bukata sediment types when compared with other sediments.

To account for the methodological uncertainties from the range of univariate statistical tests described in Section 2.3, we used a bootstrapping method [44,65] which generates the confidence limit in the final score. The results from the 1000 bootstrap runs presented in all the score charts shows that the mean score of models did not vary significantly for each different run, the ranges of 95% confidence limits were smaller for most of the models. Further, to limit the effect of spurious TSS values derived by some of the models, especially models with exponential and power functions, we filtered out any derived TSS value below a minimum of 0.001 mg/L and greater than a twice the maximum TSS concentration of each TSS model. Filtering out the spurious results can artificially inflate the final scores because only values that are within the upper and lower bounds would be considered for statistical tests. However, the possible percentage retrieval test discussed in Section 2.3.6 negate such an effect because filtering out spurious results would result in lower possible percentage retrieval and lower score in the percentage retrieval test.

4.2. TSS Model Selection Guidelines

Even though there were clearly distinct higher and lower performing TSS models from the final score chart presented in Figure 4 and Figure 5, the performance of individual models varied widely when viewed against respective water types, sediment types, and backscattering ratios. The results presented in Figure 4 and Figure 5 can be of use to the end-users who are clearly interested in TSS models that are robust enough to be used in waters for which they have little or no information of their optical and physical properties to generate TSS products. Figure 4 and Figure 5 indicate that the MODIS TSS models MOD-E6, MOD-A1, MOD-E28, MOD-A4 and MOD-E10 and the Landsat TSS models LAN-E3, LAN-A4, LAN-E9, LAN-A5 and LAN-A1 are ranked the highest in terms of likely suitability for estimating TSS concentration of unknown water types. An example of the selection of high performing TSS models using a real water dataset was demonstrated in Section 3.3.2 and it can be seen that the results varied widely among the high scoring TSS models, with MOD-E10 and LAN-E9 producing results within a MARE of 46.20% and 43.11% and other higher scoring models producing results as high as 102.59%. Considering the retrieval error of TSS concentrations from MODIS algorithms are typically reported as in the range of ~18.0% to ~61% for many studies conducted in the last decade [35,37,38,46,47], we consider the regional TSS models MOD-A1 and LAN-A1, and the empirical models MOD-E10 and LAN-E9 as being the most appropriate for the waters in the north of Western Australia.

However, readers with prior information of water and sediment types can use information provided in S11, and Table B1 and Table B2 as a guideline in selecting the model that is best suited for that particular water type. The difference in Relative Error between the high and low scoring models validated using HydroLight data and the in situ data showed that there is a huge difference between the two. The best performing model from the high scoring models shows that TSS can be estimated with a Mean RE between 69.96% and 481.82% (for different water conditions), but the low scoring model’s results can vary dramatically within a Mean RE ranging from 106.43% to 1832.79%. The high Mean RE for low scoring models does not necessarily mean that the low scoring model performs low for all waters types. The low scoring TSS model’s performance in one category or more can be significantly better than other models, but overall on average the model performs poorly when compared with high scoring models across all water types. For example, the low scoring model LAN-E22 displays the Smallest RE of 19.41% which is certainly better than the Smallest RE of most of the high scoring model’s Smallest RE. Thus, with prior knowledge of water types and bio-geochemical properties of the region, we can select a TSS model from both high and low performing TSS models presented in Table B1 and Table B2 that have higher scores in the water that are similar to the region where TSS model would be applied.

The results also showed that semi-analytic models were generally higher in ranking when compared with empirical models. The reason for most semi-analytic models performing better than empirical models can be attributed to the fact that semi-analytic models, by design, were based on radiative transfer theory [38,70] and one or more parameters were calibrated using general in situ bio-optical properties representative of a wide range of global waters [7,38].

5. Conclusions

In summary, in this study we have applied an objective methodology to compare the TSS models and their suitability in use for retrieving TSS in the absence of a regionally tuned TSS model. From the study we have identified the MODIS TSS models MOD-E6, MOD-A1, MOD-E28, MOD-A4 and MOD-E10 and the Landsat TSS models LAN-E3, LAN-A4, LAN-E9, LAN-A5 and LAN-A1 as suitable for estimating TSS concentration in waters with no prior knowledge of bio-optical or bio-geochemical properties. The results from this study highlighted the impact of “local tuning” of algorithms, showing that some low scoring models performed better than the high scoring models in one or more specific sediment, backscattering, solar zenith and water types. The results from this study can be used to ascertain which TSS models perform well in particular water types, sediment types and backscattering ratios for use in aiding the selection of a TSS model suited for use in a particular water type. In addition, the results also show that the semi-analytic TSS models are generally better than empirical TSS models in deriving TSS estimation in unknown water types.