Regional Satellite Algorithms to Estimate Chlorophyll-a and Total Suspended Matter Concentrations in Vembanad Lake

Abstract

A growing coastal population is leading to increased anthropogenic pollution that greatly affects coastal and inland water bodies, especially in the tropics. The Sustainable Development Goal-14, ‘Life below water’ emphasises the importance of conservation and sustainable use of the ocean and its resources. Pollution management practices often include monitoring of water quality using in situ observations of chlorophyll-a (chl-a) and total suspended matter (TSM). Satellite technology, including the MultiSpectral Instrument (MSI) sensor onboard Sentinel-2, enables the continuous monitoring of these variables in inland waters at high spatial and temporal resolutions. To improve the monitoring of water quality in the tropical Vembanad-Kol-Wetland (VKW) system, situated on the southwest coast of India, we present two regionally tuned satellite algorithms developed to estimate chl-a and TSM concentrations. The new algorithms estimate the chl-a and TSM concentrations from the simulated reflectance values as a function of the inherent optical properties using a forward modelling approach. The model was parameterised using the National Aeronautics and Space Administration (NASA) bio-Optical Marine Algorithm Dataset (NOMAD) and in situ measurements collected in the VKW system. To assess model performance, results were compared with in situ measurements of chl-a and TSM and other existing satellite-based models of chl-a and TSM. For satellite application, two different atmospheric correction methods (ACOLITE and POLYMER) were tested and satellite matchups were used to validate the new chl-a and TSM algorithms following standard validation procedures. The results demonstrated that the new algorithms were in good agreement with in situ observations and outperform existing chl-a and TSM algorithms. The new regional satellite algorithms can be used to monitor water quality within the VKW system to support the sustainable management under natural (cyclones, floods, rainfall, and tsunami) and anthropogenic pressures (industrial effluents, agricultural practices, recreational activities, construction, and demolishing concrete structures) and help achieve Sustainable Development Goal 14.

1. Introduction

Coastal, estuarine, and inland waters are major carbon reservoirs [1,2,3,4]; support diverse species in multiple habitats; provide a wide range of ecosystem services; and function as natural barriers for the protection of coastal areas from extreme climate events, including floods, cyclones, and tsunamis [5]. However, anthropogenic activities, such as the overexploitation of resources and pollution, have adversely affected many coastal aquatic ecosystems and the livelihood of people living in the surrounding areas [6,7,8,9,10]. The United Nations, through its Sustainable Development Goals-6 and 14, emphasise the importance of clean water and sanitation, and the sustainable use of aquatic resources [11]. In this context, monitoring water quality becomes an important step in our efforts to understand the stresses on coastal and estuarine ecosystems that are associated with anthropogenic activities, and to move towards sustainable management of their resources [12,13,14]. Some of the most vulnerable of coastal and estuarine ecosystems are in the tropics [15,16,17,18,19], which merit particular attention.

In this paper, we study the water quality of a tropical coastal, estuarine and lake system, the Vembanad-Kol-Wetland (VKW) system, situated on the southwest coast of India. VKW is a large ecosystem, adjoining some of the major towns of the state Kerala, including Alappuzha, Kochi, Kottayam, and Thrissur. The population surrounding the VKW system relies on the waterbody for many aspects of their life: provision of drinking water, fishing, agriculture, transport, tourism, and aquaculture [20]. The shipyard and port at Kochi benefit from the waterways of VKW. Despite the importance of the VKW system, heavy loads of suspended particulate and dissolved matter and nutrients (which determine the concentration of chlorophyll-a, the major pigment in phytoplankton), and the prevalence of pathogenic organisms responsible for water-borne diseases have degraded the water quality of the system [21,22,23,24]. The VKW is one of the most polluted wetland systems in India, exposed to anthropogenic stresses that include the inflow of domestic waste and sewage causing faecal contamination, and heavy metal contamination from nearby industries [25,26,27,28,29]. Construction activities along its shorelines in violation of the Indian coastal zone management regulations present yet another threat to the integrity of the VKW ecosystem [30,31]. The ecological and economic importance of the VKW system, and the multiple stresses it is currently facing, raise the need for regular monitoring of its water quality, to enable recovery and restoration, and to achieve the sustainable development goals. In situ observations of water quality are essential in this context, but they are labour-intensive. Furthermore, the complex nature of the wetland system (about 100 km long, with many inlets and meanders) makes it costly to cover the entire lake by local observations alone. Citizen science can help address some of these limitations [31,32]. However, another important avenue that we explore here is the use of satellite data as a complement to in situ observations, for monitoring water quality on a sustained and continuous basis.

Chlorophyll-a (chl-a) and total suspended matter (TSM) are two major water constituents that affect water quality, and which are amenable to remote sensing. Many models have been developed to estimate chl-a and TSM concentrations from satellite observations. Among them, the model by O’Reilly et al. [33], which produces accurate results in the open-ocean waters, has been widely employed to estimate chl-a concentrations. To estimate high concentrations of chl-a (>100 mg m⁻³), especially in turbid, productive, and inland waters, models based on the ratio of near-infrared (NIR) reflectances [34,35,36,37] have been found to be more suitable. For the estimation of TSM, models based on reflectance at red and infrared wavelengths [38,39] have performed well in sediment-dominated waters. The VKW is an optically complex system with diverse types of water appearing in different areas of the system [25]. Therefore, we cannot assume that the above-mentioned algorithms would perform equally well throughout the system.

In this paper, we evaluate the performance of existing chl-a and TSM satellite retrieval algorithms, and use a forward modelling technique to develop improved algorithms in the study area. Though earlier studies have used satellite data to investigate the dynamics of chl-a and TSM in Vembanad lake using existing algorithms [40,41,42], to our knowledge, this is the first time that such algorithms have been validated using in situ-satellite matchup data. Furthermore, satellite-derived water-leaving reflectance (${\rho }_w$) values retrieved using two atmospheric-correction algorithms are fine-tuned using in situ observations. Using Sentinel-2 MultiSpectral Instrument (MSI) data, we demonstrate that the performance of the modified chl-a and TSM algorithms is better than those of existing algorithms for the VKW system.

2. Materials and Methods

This study estimates reflectance as a function of the absorption and backscattering coefficients based on the in situ observations from Vembanad Lake and the NOMAD dataset, where the estimated reflectance was used to develop regionally tuned chlorophyll-a ($B$) and total suspended matter ($S$) algorithms suitable for Vemabanad Lake (Figure 1, discussed in detail in the following sections). The developed algorithms were applied on the MSI data after the atmospheric signal was removed based on the two different atmospheric correction methods, ACOLITE and POLYMER, and corrected using the satellite matchups to of estimated reflectance.

2.1. Study Site

Vembanad Lake (VL) is a part of the Vemband–Kol–Wetland (VKW) system, which is the second largest Ramsar site in India. The lake is ~100 km long and up to ~14 km wide and is open to the Arabian Sea in the north, near Kochi (Figure 2). Its depth varies from 1.5 to 6 m except near the port of Kochi, where a depth of 13 m is maintained by dredging to make it suitable for shipping [25]. The lake receives freshwater input from six major rivers [20]. The region experiences two monsoons: the south-west monsoon (June–September) and the north-east monsoon (October–November), which form the wet season. The dry season is during the winter (December–February) and the spring inter-monsoon (March–May) [43].

Based on salinity and hydrological conditions, VL can be divided into a brackish-water-dominated northern region (Stations VL01 to VL06, Figure 2) and a fresh-water-dominated southern region (Stations VL08 to VL13, Figure 2) [22]. The northern region also becomes dominated by freshwater during the monsoons. The Thanneermukkom Bund isolates the southern region (stations–VL10, VL11, VL12, and VL13) of the lake from the incursion of saltwater from the sea during the dry season.

2.2. In Situ Data from Vembanad Lake

In situ data, used for the development of regional satellite-algorithms (VL dataset,

Table 1. Details of datasets used. Note that the $\lambda$ represents the wavelengths—443, 490, 560, 665, and 683 or 705. The 683 nm was the longest band available for the reflectance from the NOMAD dataset, whereas the fifth band of MSI on Sentinel 2 is 705. There were $N=12$ satellite matchups available for ACOLITE (${\rho }_w{}^A$) and $N=14$ for POLYMER (${\rho }_w{}^P$) data used for the reflectance model from Vembanad Lake (VL) dataset. The superscript ‘$’ denotes that the $S$ values were estimated using Equation (11).

Vembanad Lake (VL) Dataset
$\mathit{B}$	${\mathit{a}}_{\mathit{S}}\left(400\right)$	${\mathit{a}}_{\mathit{Y}}\left(412\right)$	$\mathit{S}$	${\mathit{a}}_{\mathit{B}}\left(\mathit{\lambda }\right)$	${\mathit{a}}_{\mathit{S}}\left(\mathit{\lambda }\right)$	${\mathit{a}}_{\mathit{Y}}\left(\mathit{\lambda }\right)$	${\mathit{b}}_{\mathit{b}\mathit{t}}\left(\mathit{\lambda }\right)$	$\mathit{R}\left(\mathit{\lambda }\right)$
$228$	$228$	$228$	$162$	$228$	$228$	$228$	-	-
NOMAD dataset
$839$	$839$	$839$	${111}^{ \$}$	$839$	$839$	$839$	$111$	${228}^{R_i}$
Satellite matchup dataset
$12$	$12$	$12$	$12$	$12$	$12$	$12$	-	${12}^{R_s{}^A}$
$14$	$14$	$14$	$14$	$14$	$14$	$14$	-	${14}^{R_s{}^P}$

), were collected in VL at 13 stations (VL01–VL13, Figure 2) during 15 field campaigns between March 2018 and May 2019. Samples were collected for the analysis of bio-optical properties, including the absorption coefficients ($a$), the concentrations of chlorophyll-a biomass ($B$), and total suspended matter ($S$). Samples were filtered through 25 mm Whatman glass-fibre filters (GF/F) with a 0.7 μm pore size to collect the algal and non-algal particulate matter to determine the optical density ($D$) of phytoplankton biomass ($D_B$), non-algal suspended particles ($D_S$) and coloured dissolved organic matter (CDOM), or yellow substances ($D_Y$). The filters were analysed using a Shimadzu Ultraviolet-Visible (UV) Spectrophotometer (UV-2600) following standard methods [44,45,46]. The filter, with all material retained on it, was used to determine the optical density of particulate matter (phytoplankton biomass + non-algal suspended particles) ($D_{B+S}$). Phytoplankton pigments were removed from the filter using a methanol extraction, and the optical density $D_S$ of the residual material was measured. Then, $D_B$ was calculated by subtracting $D_S$ from $D_{B+S}$ [45,46]. For measuring $D_Y$, water samples were filtered through 0.45 and 0.2 μm nitrocellulose filters [44,46], and the optical density of the filtrate was measured [47]. The $D_B$, $D_S$ and $D_Y$ values were used to estimate their respective absorption coefficients ($a_B, a_S$, and $a_Y$) [48,49,50].

To estimate the chlorophyll-a biomass ($B$), samples were filtered through 47 mm, 0.7 μm Whatman GF/F filters and phytoplankton pigments were extracted in 90% acetone for 24 h, and the optical density $D_B$ of the extract was measured spectrophotometrically at 630, 647, and 664 nm [51,52]. Then, the chlorophyll-a biomass was calculated using the equation, $B=\left[\left(11.85\times D_B\left(664\right)\right)-\left(1.54\times D_B\left(647\right)\right)-\left(0.08\times D_B\left(630\right)\right)\right]\times \left(v/V\right)\times 1000$, where $v$ is the volume of acetone (units in ml) and $V$ is the volume of water (units in ml) filtered. Concentrations of total suspended matter ($S$) were measured gravimetrically by filtering the sample through pre-weighed 25 mm, 0.7 μm Whatman GF/F filters, which were then dried before weighing following standard techniques [53,54].

Table A1. Notations, their descriptions, and units used in the present study.

Notations	Description	Units
$a$	Absorption coefficient	m−1
$a_B$	Absorption coefficient of phytoplankton	m−1
$a_m$	Absorption coefficient of microphytoplankton	m−1
$a_m{}^{*}$	Chlorophyll-specific absorption coefficient of microphytoplankton	m2 (mg chl-a)−1
$a_n$	Absorption coefficient of nanophytoplankton	m−1
$a_n{}^{*}$	Chlorophyll-specific absorption coefficient of nanophytoplankton	m2 (mg chl-a)−1
$a_p$	Absorption coefficient of picophytoplankton	m−1
$a_p{}^{*}$	Chlorophyll-specific absorption coefficient of picophytoplankton	m2 (mg chl-a)−1
$a_S$	Absorption coefficient of non-algal suspended particles	m−1
$a_{B+S}$	Absorption coefficient of particulate matter (phytoplankton biomass + non-algal suspended particles)	m−1
$a_t$	Total absorption coefficient	m−1
$a_W$	Absorption coefficient of water	m−1
$a_Y$	Absorption coefficient of coloured dissolved organic matter, or yellow matter	m−1
$A$	ACOLITE	-
$b_b$	Back-scattering coefficient	m−1
$b_{bS}$	Back-scattering coefficient of non-algal suspended particles	m−1
$b_{bS}{}^{*}$	Specific-back-scattering coefficient of non-algal suspended particles	m2 g−1
$b_{bt}$	Total back-scattering coefficient	m−1
$b_{bW}$	Back-scattering coefficient of water	m−1
$B$	Phytoplankton biomass, in units of chlorophyll-a	mg m−3
$B_m$	Microphytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$B_n$	Nanophytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$B_p$	Picophytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$B_{p,n}$	Combined pico- and nanophytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$B_{p,n}{}^m$	Asymptotic maximum value of combined pico- and nanophytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$B_p{}^m$	Asymptotic maximum value of picophytoplankton biomass, in units of chlorophyll-a	mg chl-a m−3
$D$	Optical density	Dimensionless
$D_B$	Optical density of phytoplankton biomass	Dimensionless
$D_S$	Optical density of non-algal suspended particles	Dimensionless
$D_{B+S}$	Optical density of particulate matter (phytoplankton biomass + non-algal suspended particles)	Dimensionless
$D_Y$	Optical density of coloured dissolved organic matter	Dimensionless
$E_s$	Downwelling surface irradiance	μW cm−2 nm−1
$f$	A proportional constant for IOP-based reflectance	-
$I$	Intercept values estimated between measured and modelled data	-
$k$	Spectral slope of back-scattering coefficient of non-algal suspended particles	Dimensionless
$L_w$	Water-leaving radiance	μW cm−2 nm−1 sr−1
$m$	Fitted coefficients	Dimensionless
$m^S$	Spectral slope of absorption coefficient of non-algal suspended particles	nm−1
$m^Y$	Spectral slope of absorption coefficient of coloured dissolved organic matter	nm−1
$N$	Number of samples	-
$P$	POLYMER	-
$Q$	Bi-directional factor	sr
$R$	Reflectance	Dimensionless
$R_i$	In situ irradiance reflectance	Dimensionless
$R_e$	Estimated/Modelled reflectance	Dimensionless
$R_{rs}$	Remote-sensing reflectance	sr−1
$r^2$	Determination coefficient	-
$S$	Total suspended matter	g m−3
$s$	Slope values estimated between measured and modelled data	-
$S_p$	Slope to estimate picophytoplankton biomass	Dimensionless
$S_{p,n}$	Slope to estimate combined pico- and nanophytoplankton biomass	Dimensionless
$W$	Notation for water	-
$Y$	Notation for coloured dissolved organic matter, or yellow substances	-
$\delta$	Bias	-
$\mathsf{\Delta }$	Mean relative error	-
$\lambda$	Wavelength	nm
${\rho }_w$	Water-leaving reflectance (referred to as ‘satellite reflectance’ when it is derived from satellite)	Dimensionless
${\rho }_w^A$	Uncorrected ACOLITE-based satellite reflectance	Dimensionless
${\rho }_w^P$	Uncorrected POLYMER-based satellite reflectance	Dimensionless
${\rho }_w{}^{\prime }$	Corrected satellite reflectance	Dimensionless
${\rho }_w^{A^{\prime }}$	Corrected ACOLITE-based satellite reflectance	Dimensionless
${\rho }_w^{P^{\prime }}$	Corrected POLYMER-based satellite reflectance	Dimensionless
$\psi$	Root mean square error	-

provides details on the notations used in this study.

2.3. In Situ Data from NOMAD

In addition to the in situ measurements from VL, samples from the NASA bio-Optical Marine Algorithm Dataset (NOMAD) v2.0 (NOMAD dataset, Table 1) [55] were also used in this study for model parameterizations. NOMAD is a high-quality in situ dataset covering global waters and is used for developing ocean-colour models. For this study, the water-leaving radiance ($L_w$), downwelling surface irradiance ($E_s$), back-scattering coefficient ($b_b$), absorption coefficient of particulate matter (phytoplankton biomass + non-algal suspended particles, $a_B+a_S$); non-algal suspended particles ($a_S$); coloured dissolved organic matter ($a_Y$), and chlorophyll-a biomass ($B$) were extracted from this dataset. The remote-sensing reflectance ($R_{rs}$) was calculated as the ratio of $L_w$ to $E_s$ and converted into irradiance reflectance ($R_i$) using $R_i=Q\times R_{rs}$, where Q is set to 3.14.

2.4. Satellite Dataset

Data acquired by the MultiSpectral Instrument (MSI) onboard Sentinel-2A and -2B between 2018–2019 were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu, accessed on 6 January 2021). The two tiles of Sentinel-2 that cover Vembanad Lake were downloaded, processed, and mosaiced to generate images of the entire study area. Two processors, ACOLITE and POLYMER, were used to correct for the atmospheric effects. ACOLITE (v20190326.0; [56,57]) is based on the ‘dark spectrum fitting’ approach and POLYMER developed by Steinmetz et al. [58] removes the atmospheric signal and estimates the water-leaving reflectance using a polynomial function of the wavelength. The water-leaving reflectances (${\rho }_w$), derived from ACOLITE and POLYMER, were used to retrieve $B$ and $S$. For simplicity, hereafter, the water-leaving reflectance (${\rho }_w$) from satellites will be referred to as ‘satellite reflectance’.

2.5. In Situ Satellite Matchup Dataset

In situ data from VL was matched with the corresponding Sentinel-2 data to construct a dataset, called the Satellite matchup dataset (Table 1), which was used for the validation of regional satellite-retrieval algorithms of $B$ and $S$. Comparison between the modelled ($R_e$) (described in the later Section 2.6 and Section 2.7) and satellite (${\rho }_w$) reflectance demonstrated some systematic biases. Therefore, a correction is calculated as: ${\rho }_w{}^{\prime }=\left(l+{\rho }_w\right)\times m$, where $l$ and $m$ are fitted coefficients (Figure A1) that were applied to the processed satellite data. For both ACOLITE- and POLYMER-processed data, the intercept $l$ was wavelength independent and equal to 0.0 and 0.01, respectively. For the five MSI wavebands—443, 490, 560, 665, and 705 nm—the values of the slope $m$ were, respectively, 1.6, 2.0, 2.2, 1.6, and 1.2 for ACOLITE, and 3.5, 3.6, 3.8, 3.7, and 3.6 for POLYMER. The MSI products for the VL waters demonstrated that reflectances at bands 665 and 705 nm were frequently negative for POLYMER data before the bias correction, but were pushed to positive values with the addition of the offset. ACOLITE did not have the problem of negative values occurring at the long wavebands. Due to the small number of satellite matchups, all of them ($N=12$ and $14$ for ACOLITE, and POLYMER, respectively) were used to correct satellite reflectances.

Standard quality control methods were followed for the in situ satellite matchup [59]: (i) the covariance of the 3 × 3 pixel box (central pixel plus surrounding eight pixels) was less than 0.15 for the five Sentinel-2 reflectance bands, to ensure the homogeneity of the pixel box; (ii) more than fifty percent of the pixel box must contain valid data (i.e., was not ‘Not A Number’ or negative); (iii) the time difference between the in situ and satellite overpass was within $\pm 4$ h. For the matchup data that passed the quality criteria, the mean value of reflectances at each MSI waveband of the pixel box were computed and used as the input to the $B$ and $S$ algorithms. Based on these criteria, 12 and 14 satellite matchups were available for the ACOLITE- and POLYMER-processed satellite data, respectively.

2.6. Inherent Optical Properties

This section presents models to estimate the inherent optical properties (IOPs), and the absorption ($a$) and back-scattering ($b_b$) coefficients, of the constituents present in the water column. The absorption is contributed by water ($a_W$), phytoplankton as indexed by chlorophyll-a biomass ($a_B$), non-algal suspended particles ($a_S$), and CDOM ($a_Y$). Thus, the total absorption coefficient ($a_t$) can be expressed as:

$$a_t\left(\lambda \right)=a_W\left(\lambda \right)+a_B\left(\lambda \right)+a_S\left(\lambda \right)+a_Y\left(\lambda \right)$$

(1)

where the individual contribution of each constituent to $a_t$ is represented by the subscripts $W$, $B$, $S,$ and $Y$. The $a_W$ values were taken from Pope and Fry [60].

In this study, the method of Brewin et al. [61,62] was used to estimate $a_B\left(\lambda \right)$ values with the chlorophyll-a concentration as the input to the computation. In this approach, $a_B$ is the sum of the absorptions of three different size classes of phytoplankton, namely picophytoplankton ($a_p$ in m⁻¹), nanophytoplankton ($a_n$ in m⁻¹), and micro-phytoplankton ($a_m$ in m⁻¹):

$$a_B\left(\lambda \right)=a_p\left(\lambda \right)+a_n\left(\lambda \right)+a_m\left(\lambda \right)$$

(2)

Furthermore, each of the components of $a_B$ can be expressed as the product of the specific absorption coefficient of each phytoplankton size class and its respective biomass (in chlorophyll-a units), such that Equation (2) can be rewritten as:

$$a_B\left(\lambda \right)=\left(a_p{}^{*}\left(\lambda \right)\times B_p\right)+\left(a_n{}^{*}\left(\lambda \right)\times B_n\right)+\left(a_m{}^{*}\left(\lambda \right)\times B_m\right)$$

(3)

where $a_p{}^{*}$,$a_n{}^{*}$, and $a_m{}^{*}$ are the chlorophyll-specific absorption coefficients (in m² (mg chl-a)⁻¹) and $B_p$, $B_n$, and $B_m$ are the specific concentrations (in mg chl-a m⁻³) of the pico-, nano-, and micro-phytoplankton size classes, respectively. The specific absorption coefficients of each size class can be found in Table 2 of Brewin et al. [61]. The concentration of each size class was estimated using the following set of equations, which are all functions of the total chlorophyll-a biomass concentration, B:

$$B_p=B_p{}^m\left[1-\exp \left(-S_{p}B\right)\right]$$

(4)

$$B_{p,n}=B_{p,n}{}^m\left[1-\exp \left(-S_{p,n}B\right)\right]$$

(5)

$$B_n=B_{p,n}-B_p$$

(6)

$$B_m=B-B_{p,n}$$

(7)

where the asymptotic maximum values $B_p{}^m$ and $B_{p,n}{}^m$ are 0.13, and 0.77, respectively, and the slopes $S_p$ and $S_{p,n}$ are 6.15, and 1.22, respectively [62]. Note that $B_n$ was estimated by subtracting B_p from the combined concentration of pico- and nano-phytoplankton ($B_{p,n}$) biomass, as shown in Equation (5).

The spectral form of $a_S$ and $a_Y$ can be described by exponential functions with negative slopes, as their absorptions decrease with increasing wavelength:

$$a_S\left(\lambda \right)=a_S\left(400\right)\times \exp \left[-m^S\times \left(\lambda -400\right)\right]$$

(8)

$$a_Y\left(\lambda \right)=a_Y\left(412\right)\times \exp \left[-m^Y\times \left(\lambda -412\right)\right]$$

(9)

where $a_S\left(400\right)$ and $a_Y\left(412\right)$ are the absorption coefficient of the non-algal suspended particles at 400 nm and the CDOM absorption coefficient at 412 nm; and $m^S$ and $m^Y$ are the spectral slopes of $a_S$ and $a_Y$, for which we assumed values of 0.012 [63] and 0.015 [64,65], respectively. The values of $a_t$ can be computed once $a_B$, $a_S$ and $a_Y$ are estimated from the Equations (2), (8) and (9) in Equation (1), respectively.

When it came to modelling particle back-scattering, we faced a conundrum: neither the NOMAD nor the VL datasets included all the elements necessary for parameterising particle back-scattering as a function of total suspended matter. The VL dataset included measurements of total suspended matter ($S$) concentration and absorption coefficient of non-algal suspended particle ($a_S$), but not non-algal suspended particle back-scattering ($b_{bS}$). On the other hand, the NOMAD dataset included non-algal suspended particle back-scattering ($b_{bS}$) and absorption coefficient of non-algal suspended particle ($a_S$), but not the total suspended matter ($S$) concentration. The path we followed, therefore, was to use the common measurement of the non-algal suspended particle absorption coefficient ($a_S$) to link back-scattering by non-algal suspended particles ($b_{bS}$) and the concentration of total suspended matter ($S$), as described below.

First, the back-scattering coefficient of non-algal particles at 665 nm $b_{bS}\left(665\right)$ was expressed as a function of $a_s\left(400\right),$ using the NOMAD dataset ($N=91$, Figure 3a) as:

$$b_{bS}\left(665\right)=0.083\ast a_S{\left(400\right)}^{0.873}$$

(10)

Then, the modelled $b_{bS}\left(665\right)$ was used to estimate $S$ using the relationship between $b_{bS}\left(665\right)$ and $S$ proposed by Balasubramanian et al. [66] (Figure A2), based on data from a variety of locations and environmental conditions:

$$S=53.736\times b_{bS}{\left(665\right)}^{0.856}$$

(11)

Next, the total back-scattering coefficient ($b_{bt}$) in the NOMAD dataset was treated as the sum of the contributions to back-scattering by water ($b_{bW}$) and non-algal suspended particle ($b_{bS}$), and the contribution of phytoplankton is assumed to be relatively small compared with that of non-algal suspended particles for our study area where the non-algal suspended particle scattering dominates that of phytoplankton, such that we have:

$$b_{bt}\left(\lambda \right)=b_{bW}\left(\lambda \right)+b_{bS}\left(\lambda \right)$$

(12)

The spectral values of $b_{bW}$ were taken from Morel [67]. The spectral variation in $b_{bS}$ derived from the equation above was then represented as [68]:

$$b_{bS}\left(\lambda \right)=b_{bS}\left(665\right)\times {\left(\frac{\lambda }{665}\right)}^{-k}$$

(13)

where the exponent $k$ determines the spectral variations in $b_{bS}$, which was estimated as 1.32 by fitting Equation (13) to $b_{bS}\left(\lambda \right)$. For the NOMAD dataset, $b_{bS}\left(665\right)$ was computed using Equation (10), and tested against in situ data (Figure 3b). Due to the lack of in situ $b_b$ from VL, the $b_b$ model was only compared with the NOMAD dataset. For the Vembanad Lake dataset, we estimated the back-scattering coefficient of non-algal suspended particle ($b_{bS}\left(665\right)$) at 665 nm from measured S, by inverting Equation (11):

$$b_{bS}\left(665\right)={\left(\frac{S}{53.736}\right)}^{\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$0.856$}\right.\right)}$$

(14)

The computed $b_{bS}$ values were then used in Equation (12) to estimate $b_{bt}$.

Table 2. Table showing the different satellite-retrieval algorithms and their coefficients used in the present work to estimate chlorophyll-a ($B$) and total suspended matter ($S$) concentrations.

Algorithms	Equations	Variables	References
B1	$\begin{array}{c}B={10}^{\left[0.19-\left(1.81X\right)+\left(1.97X^2\right)-\left(2.56X^3\right)\right]}\\ -0.722\end{array}$	$X=lo{g}_{10}\left[R_{rs}\left(443\right)/R_{rs}\left(560\right)\right]$	O’Reilly et al. [33] with parameters from Franz et al. [69] and Vanhellemont and Ruddick [70]
B2	$B={\left[\left(35.8X\right)-19.3\right]}^{1.12}$	$X=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	Gilerson et al. [36]
B3	$B=40.6-\left(138X\right)+\left(122X^2\right)$	$X=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	This study
B4	$\begin{array}{c}B=-55.0-\left(56.1X_1\right)+\left(50.3X_2\right)\\ -\left(20.8X_3\right)+\left(89.0X_4\right)\end{array}$	$X_1=R_{rs}\left(443\right)/R_{rs}\left(560\right)$ $X_2=R_{rs}\left(490\right)/R_{rs}\left(560\right)$ $X_3=R_{rs}\left(665\right)/R_{rs}\left(560\right)$ $X_4=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	This study
S1	$S=-1.91+\left(1140X\right)$	$X=R_{rs}\left(665\right)$	Miller and McKee [38]
S2	$S=1.74+\frac{\left(356X\right)}{1-\left(X/0.1728\right)}$	$X=R\left(665\right)$	Nechad et al. [39]
S3	$S=1.1+\left(332X\right)+\left(7188X^2\right)$	$X=R_{rs}\left(665\right)$	This study

Table 2. Table showing the different satellite-retrieval algorithms and their coefficients used in the present work to estimate chlorophyll-a ($B$) and total suspended matter ($S$) concentrations.

Algorithms	Equations	Variables	References
B1	$\begin{array}{c}B={10}^{\left[0.19-\left(1.81X\right)+\left(1.97X^2\right)-\left(2.56X^3\right)\right]}\\ -0.722\end{array}$	$X=lo{g}_{10}\left[R_{rs}\left(443\right)/R_{rs}\left(560\right)\right]$	O’Reilly et al. [33] with parameters from Franz et al. [69] and Vanhellemont and Ruddick [70]
B2	$B={\left[\left(35.8X\right)-19.3\right]}^{1.12}$	$X=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	Gilerson et al. [36]
B3	$B=40.6-\left(138X\right)+\left(122X^2\right)$	$X=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	This study
B4	$\begin{array}{c}B=-55.0-\left(56.1X_1\right)+\left(50.3X_2\right)\\ -\left(20.8X_3\right)+\left(89.0X_4\right)\end{array}$	$X_1=R_{rs}\left(443\right)/R_{rs}\left(560\right)$ $X_2=R_{rs}\left(490\right)/R_{rs}\left(560\right)$ $X_3=R_{rs}\left(665\right)/R_{rs}\left(560\right)$ $X_4=R_{rs}\left(705\right)/R_{rs}\left(665\right)$	This study
S1	$S=-1.91+\left(1140X\right)$	$X=R_{rs}\left(665\right)$	Miller and McKee [38]
S2	$S=1.74+\frac{\left(356X\right)}{1-\left(X/0.1728\right)}$	$X=R\left(665\right)$	Nechad et al. [39]
S3	$S=1.1+\left(332X\right)+\left(7188X^2\right)$	$X=R_{rs}\left(665\right)$	This study

2.7. Forward Reflectance Model

The IOPs were used to model reflectance $R_e$ (Morel and Prieur [71]; see also Sathyendranath and Platt [72]) as:

$$R_e\left(\lambda \right)=0.33\times \left[\frac{b_{bt}\left(\lambda \right)}{\left(a_t\left(\lambda \right)+b_{bt}\left(\lambda \right)\right)}\right].$$

(15)

2.8. Chl-a and TSM Satellite Algorithms

The computed $R_e$ values were used to develop regional satellite algorithms for retrieval of $B$ and $S$ in VL. Several existing and newly-developed algorithms (Table 2) were tested to evaluate the efficiency of $B$ and $S$ retrieval in VL using the VL dataset ($N=162$, Table 1).

To estimate $B$, the algorithms of O’Reilly et al. [33], based on a blue-green reflectance ratio (Algorithm B1; Table 2) and Gilerson et al. [36], based on a NIR reflectance band ratio (Algorithm B2; Table 2) were examined, as they are both used extensively in the literature for open-ocean waters [33], and for more complex coastal waters [36]. The model coefficients used by Gilerson et al. [36] were tuned to improve the performance of this algorithm (Algorithm B3; Table 2) using the VL dataset. In addition, by experimenting with different combinations of reflectance bands, a new algorithm based on multi-linear regression (BGR, Algorithm B4, Table 2) was introduced that used all five reflectance bands from Sentinel-2.

For the estimation of $S$ in VL, the algorithms of Miller and McKee [38] (Algorithm S1, Table 2) and Nechad et al. [39] (Algorithm S2, Table 2) were tested. It has been reported that reflectance in the red bands increases with $S$, especially in sediment-dominated turbid waters [73]. In VL, a polynomial relationship was observed (Algorithm S3, Table 2) between $R_e\left(665\right)$ and in situ $S$:

$$S\propto \exp \left(R_e\left(665\right)\right)$$

(16)

and this relationship was also explored, to estimate $S$ from satellite data.

The performance of all $B$ and $S$ algorithms (Table 2) were assessed based on the following statistical metrics: determination coefficient ($r^2$), root mean square error ($\psi$), bias ($\delta$), mean relative error ($∆$), and the slope ($s$), and the intercept ($I$) of a linear fit between the model and observations.

3. Results

3.1. Inherent Optical Properties

The comparison between the in situ and model-derived values of $a_B$, $a_S$, $a_Y$, and $a_t$ at the five key wavelengths demonstrated that the modelled absorption values were consistent with in situ measurements, both from the VL ($N=228$) and NOMAD ($N=839$) datasets (Figure 4). The model-derived $a_Y$ showed more scatter at 665 and 683 nm compared with the $a_B$ and $a_S$ estimates at the same wavelengths (Figure 4i,j), but the effect on the total absorption values ($a_t$) is insignificant, as the magnitude of $a_Y$ is small compared with the absorption coefficients of the other constituents ($a_W$, $a_B$, and $a_S$) at these wavelengths (Figure 4s,t). The estimated $b_{bt}$ values are in close agreement with in situ values from the NOMAD dataset ($N=111$) (Figure 5).

3.2. Reflectance Model

The model-derived $a_t$ and $b_{bt}$ were used to estimate the $R_e$ using a forward reflectance model. To assess the performance of the $R$ model, the $R_e$ values were compared with in situ ($R_i$) and satellite (${\rho }_w$) measured reflectances. The results demonstrated that the $R_e$ values are in close agreement with $R_i$ from the NOMAD dataset (Figure 6, grey circles) demonstrating the quality of the $a$ and $b_b$ values that were used to compute the in situ reflectance values. Spectral comparison of the data demonstrated that the satellite-derived ${\rho }_w$ (${\rho }_w^A$ and ${\rho }_w^P$) in both ACOLITE- and POLYMER-processed datasets underestimated $R_e$ (Figure A3a,b,d,e), and these data thus required a correction. The application of the empirical correction term reduced the differences between the $R_e$ and ${\rho }_w$ (Figure 6, bright coloured stars and triangles). After correction, the spectral comparison between ${\rho }_w{}^{\prime }$ and $R_e$ demonstrated a better agreement between modelled and satellite-derived reflectances, and their spectral curvature matched that of $R_e$ reasonably well (Figure A3c–f).

3.3. Chl-a and TSM Satellite Algorithms

Existing algorithms for chlorophyll-a retrieval that were tested (Algorithms B1 and B2, Table 2; Figure 7a,b) did not appear to be suitable for the study area, even after regional tuning. The Gilerson et al. [36] algorithm with tuned coefficients (Algorithm B3) showed improvement at concentrations > 10 mg m⁻³, but high variance remained at lower concentrations (Figure 7c). The locally tuned multi-linear regression algorithm (Algorithm B4, Table 2, and Figure 7d) performed best in comparison with the other algorithms tested in this study. Hence, algorithm B4 was chosen as the chl-a algorithm for satellite application in VL.

For the estimation of total suspended matter in VL, good correlations ($r^2$ > 0.85) were obtained for the algorithms of Miller and McKee [38] (Algorithm S1, Table 2) and Nechad et al. [39] (Algorithm S2, Table 2) when compared with in situ measurements, but the tuning of both algorithms was required to improve the estimation of $S$ (Figure 7e,f). The tuned Algorithm S3, in which $S$ is estimated as a function of $R\left(665\right)$, demonstrated better performance compared with the other $S$ algorithms (Figure 7g). Therefore, Algorithm S3 was selected for the satellite retrieval of $S$ in VL.

3.4. Application of Regionally Tuned Satellite Retrieval Algorithms

To identify the best atmospheric correction technique suitable for estimations of chlorophyll-a biomass ($B$) and total suspended matter ($S$) in VL, the concentrations retrieved using the empirically corrected satellite reflectance values were compared with the satellite matchup dataset (Table 1) (Figure 1). The regionally tuned Algorithms B4 and S3 were used for the retrieval of $B$ and $S$, respectively. Statistical analysis demonstrated that the satellite-based $B$ from the ACOLITE-processed data and the satellite-based $S$ from the POLYMER-processed data compared best with in situ measurements (Figure 8a,d). We note that the estimates of $B$ and $S$ in both ACOLITE- and POLYMER-processed data demonstrated four outliers, all of them for samples collected on the 4th of January 2019 (Figure 8a,b). For these matchups, the corresponding in situ $S$ values are much higher than for other samples. The reason for this anomaly remains unclear.

As examples of products that can be derived from the regionally tuned algorithms developed in this study, we demonstrate the spatial distribution of satellite-derived $B$ and $S$ in VL for the 4th of January and the 25th of March 2019 (Figure 9). We observe that $B$ ranged from 3–15 mg m⁻³ across the lake, with occasionally higher concentrations up to 60–80 mg m⁻³ in the north (Figure 9a,b), and concentrations between 6–9 mg m⁻³ near the Thanneermukkom Bund. Concentrations of $S$ ranged from 9–15 g m⁻³ on the 4th of January 2019 (Figure 9c). A maximum concentration of $S$ of nearly 45 g m⁻³ was observed for some pixels, especially near Kochi and the south of the lake on the 25th of March 2019 (Figure 9d). High values of $B$ and $S$ in the southeast side of VL may be related to nutrient and sediment inputs from the adjacent paddy fields.

4. Discussion

4.1. Satellite Products for Sustainable Development Goals

For decades, Vembanad Lake has been impacted by anthropogenic stresses that affect its ecological health and socioeconomic status [25,74,75]. To move towards the sustainable use of Vembanad Lake, a detailed water quality monitoring programme is required [8,76,77]. In this study, we addressed the estimation of two key indicators of water quality, chlorophyll-a and total suspended matter, from remote sensing observations to support such a monitoring programme. As an indicator of ecosystem health, chlorophyll-a concentration can be used to study phytoplankton dynamics in general as well as to monitor harmful algal blooms resulting from eutrophication that affect ecosystem health [78,79]. The total suspended matter affects the turbidity of the lake with high concentrations leading to reduced light penetration, and hence decreasing the light available for photosynthesis by phytoplankton. Earlier studies in other regions have reported that high levels of $S$ can also affect zooplankton and fish survival [80]. Satellite-based routine monitoring of these variables for the assessment of water quality can be efficient and cost-effective compared with in situ observations. For instance, there are studies that assessed the variability in water quality during the COVID-19 lockdown using satellite data in Vembanad Lake and other aquatic systems in India [41,42,81]. However, the assessment of water quality using satellite data can be improved when complementary in situ observations are available, as was the case in this study.

4.2. Challenges with Quantitative Water Quality Measurements from Satellites over Vembanad Lake

Although our study relied on an extensive field campaign in Vembanad Lake that measured various water quality indicators over an annual cycle, only a small number of matchup data points were available. Matching satellite and in situ data in Vembanad Lake was difficult because: (i) logistical issues and other conflicting requirements affected the scheduling of field campaigns during satellite overpasses, even though Sentinel-2 has a high overpass frequency of 5 days; (ii) carrying out simultaneous in situ measurements at 13 stations spread across the 100 km length of the lake is difficult; in fact, it took 2 days to complete one transect; and (iii) frequent cloud cover, especially during the monsoon seasons, limited satellite retrievals. Therefore, we employed a modelling approach to estimate the reflectance values based on IOPs that were either observed or estimated from in situ observations. The simulated reflectances constituted a large dataset (N = 162, VL dataset) that could be used to improve satellite algorithms. Using this approach, we demonstrated that the locally tuned algorithms performed better than previously established models [33,36,38,39] for the satellite retrievals of the concentrations of $B$ and $S$ (discussed later in this section).

The geometry of Vembanad Lake, which is long and narrow, with a highly indented shoreline and multiple waterways; the proximity to land and the potential for water-pixel contamination by neighbouring pixels; and above all, the optically complex nature of the waters; all rendered the study area particularly problematic from a remote-sensing point of view. The spatial resolution of satellite sensors such as the Ocean Colour and Land Imager (OLCI) on Sentinel-3 (300 m) was insufficient to resolve the spatial distribution of water quality properties in the lake. The solution was to make use of the MultiSpectral Instrument (MSI) on Sentinel-2, which is designed primarily for land applications, but has a spatial resolution of 10m in the visible domain. While the spatial resolution was ideal for the purposes of the study, we had to compromise on the spectral quality and resolution, and adapt in-water algorithms for the application with MSI data.

Phytoplankton-size-based models have been demonstrated to be useful in reproducing the spectral values of the absorption coefficient of phytoplankton [61,62,82,83,84,85]. Although the size-class model of Brewin et al. [61,62] worked well in the estimation of $a_B\left(\lambda \right)$ values for our study site, we recognise that the absorption values of non-algal suspended particles vary with sediment type, composition, and concentration [86,87] and, consequently, with aquatic environments [88]. For $a_Y$, the model values demonstrated some scatter in the longer bands, especially for samples from the VL waters (Figure 4). The use of a constant slope for the absorption spectrum of yellow substances ($Y$) is known to produce high errors at long wavelengths when concentrations of $Y$ are high or when different sources contribute to $Y$ [89]. Previous studies have reported that the slope of $a_Y$ can vary between 0.01 and 0.02 with the composition and source of the substance [64,65,90,91,92,93]. This is also potentially true for VL, where phytoplankton, macrophytes, and materials of terrestrial origin from river runoff could all contribute to the dissolved organic matter in the water. However, the overall effect of the deviation in $a_Y$ on total absorption is minimal, since concentrations of $Y$ are low in VL compared with $B$ and $S$.

The atmospheric correction of Sentinel-2 data was another challenge. In addition to the complex nature of the water body itself, problems associated with poor air quality, high aerosol content, and high humidity, which could also have affected the aerosol properties, eventually makes it difficult to perform the atmospheric correction. We therefore implemented two well-known atmospheric correction procedures that are known to perform well in complex situations [66,94,95]. Reflectance values retrieved using both procedures—ACOLITE and POLYMER—demonstrated a significant offset when compared with simulated in situ reflectances, and the data needed a spectral bias correction to bring the magnitudes and spectral shapes to within reasonable values. Though the underlying cause responsible for the offset remains unknown, one might speculate that the algorithms were unable to distinguish fully between high aerosols in air and scattering particles in water.

4.3. Chlorophyll-a and Total Suspended Matter in Vembanad Lake

For estimating chl-a, one of the algorithms tested here was the blue-green algorithm of O’Reilly et al. [33], which is designed for use in open-ocean, Case-1 waters, where phytoplankton and chl-a can be treated as the primary agents responsible for variations in optical properties, with other substances, when present, covarying with chl-a. The poor performance of this algorithm is therefore expected in Case-2 waters, such as those of Vembanad Lake, where high absorption and scattering by particulate and dissolved organic matter, varying independently of phytoplankton, can degrade algorithm performance.

Mittenzwey et al. [34] and Gitelson [96], among others, proposed using the NIR reflectance bands, which contain information on phytoplankton, to estimate chl-a in turbid, but productive waters. The NIR chlorophyll-a algorithm of Gilerson et al. [36], after tuning for Vembanad Lake (Algorithm B3), performed well at high concentrations. At low chl-a concentrations, the phytoplankton signal is generally weak in the NIR region, introducing noise in the algorithm estimates. The conditions in Vembanad Lake called for an algorithm that works well in Case 2 environments, at both low and high chl-a concentrations. The multi-linear regression algorithm (Algorithm B4), which uses a combination of the five reflectance bands, performed effectively in the retrieval of $B$ with high accuracy in VL waters. However, the use of short (blue and green) wavebands in satellite retrieval algorithms can also be problematic, as the reflectances at these bands are often erroneous after atmospheric corrections [97]. Hence, the performance of the multi-band algorithm should also be tested frequently and modified accordingly, as more validation data become available.

Sathyendranath et al. [68], Miller and McKee [38] and Nechad et al. [39], and others have demonstrated that reflectances in the red wavebands are suitable for the retrieval of the total suspended matter concentration. Following these earlier studies, we chose to estimate $S$ as a function of reflectance at 665 nm, because there is a proportional increase in reflectance at this band, with an increase in $S$. However, the model coefficients needed to be tuned for Vembanad Lake to improve the quantitative retrieval. The need for re-parameterising the model may be attributed to the impact of the size, shape, and nature of the inorganic suspended particles on their optical properties [86,87]. Classifying the pixels into water types and applying class-specific algorithms [66,98] could be a potential avenue, to avoid the need for re-parameterising algorithms for different regions. However, the lack of $b_b$ measurements from the study site did not allow us to explore such algorithms in our study.

5. Conclusions

In this study, we presented site-specific products for chlorophyll-a and total suspended matter concentrations in Vembanad Lake, derived using high-spatial-resolution Sentinel-2 data. To this end, remote-sensing reflectances were simulated using a forward modelling approach that uses absorption, and back-scattering coefficients as inputs. Then, the simulated reflectances were used to test algorithms for estimating the chlorophyll-a biomass ($B$) and total suspended matter ($S$) concentrations. We used the NOMAD and Vembanad Lake in situ datasets for developing the reflectance model and for validation. Two different atmospheric correction techniques—ACOLITE and POLYMER—were examined for the application of regionally tuned chlorophyll-a and total suspended matter algorithms in Vembanad Lake. The satellite-retrieved products were validated against in situ matchups. When compared with other commonly-used algorithms, in situ data and values estimated from the site-specific algorithms proposed in the current study were in good agreement despite the fact that the number of matchups was small. Although ACOLITE and POLYMER atmospheric correction procedures produced similar results, we recommend the use of ACOLITE-based reflectances for the estimation of chlorophyll-a and POLYMER-based reflectances for the estimation of total suspended matter concentrations. The statistical performance of both ACOLITE and POLYMER were comparable for the retrieval of total suspended matter concentration; POLYMER had a small advantage with respect to spatial coverage. This study indicates that the algorithms developed here can be used for the routine monitoring of water quality in the management of Vembanad Lake, which is valuable in the assessment of the effects of anthropogenic activities on human and aquatic life. In particular, the satellite-derived chlorophyll-a, when combined with photosynthetically available radiation (PAR), can support studies related to primary production. Systematic studies of chlorophyll-a concentration are also useful for monitoring phytoplankton blooms, including harmful algal blooms. A recent study also demonstrated [22] that there is a complex relationship between chlorophyll-a concentration and the presence of Vibrio cholerae bacteria in Vembanad Lake. Sediment dynamics, critical for the local harbour, can now be monitored from the satellite, and could help investigate the impact of dredging activities to maintain the water depth in shipping channels. Furthermore, such observations can be used to infer the flushing rate of water in the lacustrine–estuarine system, a critical component among various factors that maintain the health of the lake ecosystem. Such observations can also help to monitor erosion and deposition along the shoreline. Thus, this study supports the United Nation’s sustainable development goals—3, 6, and 14 by monitoring the water quality and by helping to maintain good human and animal health both in water and along the shores of Vembanad Lake. Future efforts include calibrating the model coefficients with more satellite matchups; implementing an improved atmospheric correction method; using in situ back-scattering coefficients and eliminating uncertainties associated with using a constant bidirectional factor (f/Q), which could help in enhancing the model performance.