Abstract
Regular monitoring of inland water bodies is essential to identify the areas of deteriorating water quality and to take measures to curb the impairment. The present study aims to develop a semi-analytical model for estimating Trophic State Index (TSI) from Sentinel 2 (S2) imagery. A semi-analytical algorithm, QAA-v6, is parameterized for S2 data (referred to as QAA-S2) based on the correlation between the retrieved inherent optical properties (IOPs) from S2 imagery and the measured TSI data of Milford Lake in the USA. The accuracy of estimation of TSI by this modified model from three different lakes located in USA and from one lake in India has increased by almost 50%, when compared with that of QAA-v6. A correlation analysis of the retrieved IOPs from QAA-S2 using the outputs of three atmospheric correction processors (ACPs) (namely C2RCC, Acolite and Sen2Cor) was carried out and C2RCC gave the least mean absolute percentage error (MAPE<8%) for TSI estimation. TSI estimation using single-scattering albedo (u) at B5 and B6 bands of S2 was reasonable (MAPE<12%) to mark them as computationally efficient estimators of TSI. These results indicate the usefulness and transferability of the QAA-S2 to the various parts of the globe for estimating TSI.
HIGHLIGHTS
A parameterized semi-analytical model (QAA-S2) is proposed for the estimation of TSI from Sentinel 2 imagery.
C2RCC-retrieved reflectance data provide the best results for the model; three atmospheric correction processors, namely C2RCC, Acolite and Sen2Cor, are compared.
Single-scattering albedo at B5 and B6 of S2L1C and indices derived from backscattering coefficients are identified as good estimators for TSI.
Graphical Abstract
INTRODUCTION
Remote sensing has been of interest to researchers and environmentalists owing to its ability to get a synoptic view and a real-time source of information about the topic of interest (García 2021). The potential of satellite remote sensing has been effectively used for water quality monitoring of oceans. Retrieving of bathymetric information, water colour, temperature, salinity and concentration of optically active constituents (OACs) from oceans is carried out competently with the aid of ocean monitoring satellites. Satellite remote sensing has been used on coastal waters and estuaries for wetland mapping, shoreline changes (Cham et al. 2020), bathymetric studies and submerged vegetation studies. However, its application in the area of water pollution monitoring is still in the developing stage and has not reached satisfactory levels.
The models for retrieving information about water quality from satellite imagery can be grouped into empirical, semi-analytical and analytical. Semi-analytical algorithms lie midway between the empirical models and analytical models. The Quasi-Analytical Algorithm (QAA), developed by Lee et al. (2002), is a semi-analytical model that estimates the concentrations of OACs (namely Chlorophyll-a (Chl-a), Coloured Dissolved Organic Matter (CDOM) and non-algal particles (NAPs)) present in a water body from the above-surface remote sensing reflectance (Rrs) sensed from it. In that algorithm, the radiative transfer equations (RTEs) are used to retrieve the inherent optical properties (IOPs) of OACs from Rrs (Zhu & Yu 2013). Ocean monitoring satellites like MODIS and MERIS have been developed with the spectral resolutions necessary for these algorithms and they have been effectively used as a real-time water quality monitoring tool for oceanic water bodies. But the spatial resolution of these satellite imagery is inadequate to capture the variability in inland water bodies. The other challenges in using the QAA for coastal and inland water bodies are the complex uncorrelated relationships between the constituents like sediments, organic and inorganic matter present in them and their unpredictable temporal and spatial variabilities.
Another major limitation in using satellite imagery as input in the QAA is their atmospheric interference that accounts for 80–90% of the received signal. The errors in the retrieval of the corrected spectrum greatly influence the accuracy of the model. The inability of the spatial resolution of ocean colour sensors to capture the variability in inland water bodies and the absence of fine spectral bands required for the QAA in multispectral land monitoring satellite missions have limited the employment of this algorithm for inland quality monitoring. Most of the studies used in situ spectral measurements for the development and validation of the QAA for inland water bodies (Warren et al. 2019). Hence, it is the need of the hour to parametrize the QAA to suit the resolutions of multispectral satellite imagery without compromising its efficiency.
The QAA is a set of analytical and empirical steps developed from in situ measurements and online databases (International Ocean Colour Coordinating Group (IOCCG) synthetic data and NASA Bio-Optical Marine Algorithm Data Set (NOMAD) databases) of Rrs and IOPs. The algorithm generalizes the spectral behaviour of the OACs and any deviation from this behaviour affects the accuracy of the model. A thorough knowledge of the spectral behaviour, composition and origin of the OAC is necessary for the accurate application of the algorithm. In other words, the concentration, composition and inter-relationship of OACs are to be well understood. Hence, the partitioning of the IOPs from the total IOP is a major uncertainty in this method (Watanabe et al. 2016; Li et al. 2017). Instead of using the individual concentrations of OAC as the parameter for quality monitoring, the Trophic State Index (TSI), an overall index representing the impact of all the constituents present in the water body, can be utilized as the quality indicator. Theoretically, it should be correlated to the total IOPs of the water body. The TSI was developed by Carlson (1977) and has been adopted by ecologists and limnologists as a tool for monitoring and managing inland water bodies (Papoutsa et al. 2014). Using the TSI as the water quality parameter helps in evading the drawback of the extensive data required for the portioning of the total IOP into that of the respective OACs. The correlation between the total IOP of the water body and the TSI may be utilized for parametrizing QAA.
The TSI can be computed from log-transformed values of Secchi depth (SD), Total Nitrogen (TN), Total Phosphorus (TP) or Chlorophyll-a (Chl-a). Though all the four parameters are correlated, only SD and Chl-a are optically active and detectable by remote sensing techniques. Retrieval of Chl-a by QAA requires Rrs at 412 nm (unavailable in S2) and 443 nm (blue band). Studies on the efficiency of atmospheric correction processors (ACPs) have mentioned the higher uncertainties in the retrieval of the blue band (Warren et al. 2019; Renosh et al. 2020). Hence, it is better to avoid these bands while using satellite imageries for QAAs on water bodies. SD is highly correlated to turbidity and suspended matter, both of which are detectable from the red bands (Zheng & DiGiacomo 2017). Hence, TSI(SD) (TSI derived from SD) would be correlated to IOPs at these bands. Thus, the possibility of a relation between TSI(SD) (hereafter mentioned as TSI) and IOPs (namely absorption coefficient ‘a’, backscattering coefficient ‘b’ and single scattering albedo ‘u’) is worth exploring. A few studies have attempted to use remote sensing for quantifying TSI with airborne sensors, Landsat imagery and IRS LISS imagery (Papoutsa et al. 2014). However, they were focused on developing multiple regression models between the TSI and top of atmosphere (TOA) reflectance data of the study area, and these models had limited transferability to other areas. Moreover, the possibility of deriving TSI by semi-analytical methods has been rarely explored.
Sentinel 2 (S2), commissioned in 2016, has been of interest to researchers because of its improved spatial, spectral and temporal resolutions suitable for monitoring of inland water bodies. The red-edge bands (B5-704 and B6-740 nm), absent in most other terrestrial satellites, is an added advantage of S2 and the potential of these bands for water quality monitoring has not been adequately explored.
From these discussions, it is revealed that a semi-analytical algorithm for estimating the TSI of an inland water body from S2 imagery has not been previously developed and/or validated. Also, the QAA has not been parametrized for S2 imagery. Hence, the objective of the present study is to develop a model for mapping the trophic status (TSI) of an inland water body by using the semi-analytically derived IOP from Sentinel 2 Level 1C (S2L1C) imagery. The study focuses on parameterizing QAA for S2L1C data so that the capabilities of S2 may be utilized for estimating the quality of inland water bodies. As the QAA requires bottom of atmosphere (BOA) reflectance and the S2L1C gives the TOA reflectance, a suitable ACP that would provide the closest above-surface (BOA) reflectance data needs to be identified. BOA reflectance (Rrs) from three ACPs, namely Case 2 Regional Coast Colour (C2RCC), Acolite and Sentinel 2 Correction (Sen2Cor), are used as inputs in the model and the difference in the retrieved IOPs using these three datasets is analysed to explore the effect of the ACP on the efficiency of the algorithm.
MATERIALS AND METHODS
Study area
The area chosen for study is Milford Lake, Kansas, USA (Figure 1(a)) which is a man-made reservoir, fed by the Republican River. It has a surface area of 64 km2 and a catchment area of 6,200 km2. The catchment consists of cropland (51%), grassland (36%) and urban area (13%). The lake has been categorized by the Kansas Department of Health and Environment (KDHE) as a water body that needs restoration of water quality (Nejadhashemi et al. 2009). Nutrients from fertilizer runoff are the major cause of the water quality deterioration of the lake.
Three lakes in the USA (Conroe Lake, Canyon Lake and Texoma Lake) (Figure 1(b)–1(d)) and a lake in India (Vembanad Lake, Kerala) (Figure 1(e)) were chosen for validation. The lakes in the USA were chosen based on the availability of turbidity measurements from the USGS dataset and cloud-free S2 L1C imagery.
The area chosen for study from India is the Vembanad Lake located in Kerala. The lake is the longest (96 km) in India, the largest in Kerala and one of the three protected Ramsar sites in Kerala. Various studies conducted on the environmental status of the lake have revealed its deteriorating condition due to anthropogenic activities. In situ measurements were taken from two regions of this lake: Kochi Kayal in the north and Punnamada Kayal in the south (Figure 1(e)).
The details of these lakes with the dates of sample collection, satellite sensing date and the number of sample points are given in Table 1.
Details of study area, dates of sample collection and satellite sensing date
Study area . | Surface area of lake (km2) . | No. of sample points . | Date of sample collection . | Sensing date of image . |
---|---|---|---|---|
Milford Lake, Kansas | 64 | 26 | 10 July 2018 | 9 July 2018 |
26 | 9 August 2018 | 3 August 2018 | ||
13 | 16 October 2018 | 22 October 2018 | ||
13 | 17 October 2018 | |||
Texoma Lake, Texas | 360 | 6 | 3–4 August 2020 | 2 August 2020 |
Conroe Lake, Texas | 84 | 7 | 4 February 2020 | 6 February 2020 |
Canyon Lake, California | 1.55 | 3 | 12 June 2020 | 1 June 2019 |
Kochi Kayal | 2,033a | 24 | 26 January 2021 | 26 January 2021 |
Punnamada Kayal | 16 | 13 March 2021 | 14 March 2021 |
Study area . | Surface area of lake (km2) . | No. of sample points . | Date of sample collection . | Sensing date of image . |
---|---|---|---|---|
Milford Lake, Kansas | 64 | 26 | 10 July 2018 | 9 July 2018 |
26 | 9 August 2018 | 3 August 2018 | ||
13 | 16 October 2018 | 22 October 2018 | ||
13 | 17 October 2018 | |||
Texoma Lake, Texas | 360 | 6 | 3–4 August 2020 | 2 August 2020 |
Conroe Lake, Texas | 84 | 7 | 4 February 2020 | 6 February 2020 |
Canyon Lake, California | 1.55 | 3 | 12 June 2020 | 1 June 2019 |
Kochi Kayal | 2,033a | 24 | 26 January 2021 | 26 January 2021 |
Punnamada Kayal | 16 | 13 March 2021 | 14 March 2021 |
aSurface area of the Vembanad Lake.
Satellite data
S2L1C imagery was used for the study. S2 is a multispectral satellite mission commissioned by the European Space Agency (ESA) in July 2015. The imagery provides radiometrically and geometrically corrected TOA reflectance in 13 bands: four in the visible region, six in the near-infrared (NIR) region and three in the short-wave infrared (SWIR) region. The bands in the visible region, B2–B4 and the NIR band B8 have a spatial resolution of 10 m. Bands B5–B7 have 20 m spatial resolution, whereas B1, B9 and B10 have 60 m spatial resolution. The bands B8A and B9 are water vapour retrieval bands. Three bands in red-edge B5–B7 (central wavelength at 704, 740 and 780 nm) with a spatial resolution of 20 m is an advantage of S2L1C imagery. The high temporal resolution of 5 days (2–3 days for mid-latitudes) combined with its better spatial resolution and free availability has made it an extremely useful resource for research purposes. Sentinel data are freely downloadable at Copernicus Open Access Hub (https://scihub.copernicus.eu; accessed 21 June 2021).
The pre-processing and analysis of the imagery were carried out by using SeNtinel Application Platform (SNAP) software (version 8.0.2) developed by the ESA. SNAP consists of several open-source toolboxes and bundled packages for processing MERIS, Sentinel and Landsat imageries. The pre-processing steps include resampling the imagery to 10 m spatial resolution and clipping the area of interest from the downloaded images.
Three ACPs were used in this study: C2RCC, Acolite and Sen2Cor. Sen2Cor and C2RCC are available as plugins in SNAP. Acolite is a stand-alone processor.
Water quality data
Water quality data of all the lakes in the USA were downloaded from the USGS National Water Information System (https://nwis.waterdata.usgs.gov; accessed 21 December 2020). The in situ measurements of Milford Lake were used for the development of the model. The data from other lakes are used for validation.
METHODS
Atmospheric correction processor
C2RCC developed by Doerffer & Schiller (2007) was available in S2 toolbox under the thematic water processing menu of SNAP. The C2RCC processor uses artificial neural networks (ANNs) based on a large database of simulated water-leaving radiances and related TOA radiances. Around 5 million cases generated from in situ measurements obtained from NOMAD database, Coastcolour database and simulated results from the radiative transfer model (HYDROLIGHT) were used by Doerffer & Schiller (2007) for developing the C2RCC model (Brockmann et al. 2016). The water-leaving reflectance in six bands (B1–B6) is available as the output from the processor.
Acolite is a stand-alone processor for S2 and Landsat 5/7/8. In this model, Rayleigh correction is based on a lookup table generated by 6SV (Second Simulation of the Satellite Signal in the Solar Spectrum, Vermote) and the aerosol correction is carried out by the Dark Spectrum Fitting (DSF) or Exponential (EXP) method (Vanhellemont & Ruddick 2016). The output from Acolite can be viewed and analysed in SNAP.
The Sen2Cor processor is based on Dense Dark Vegetation (DDV) and it is developed for land products. The algorithm estimates the aerosol thickness using the dark vegetation pixels at 550 nm based on the correlation between B12 (SWIR) and visible bands. L2A images provided by Sentinel use the Sen2Cor algorithm for atmospheric correction (Main-Knorn et al. 2017).
Deriving IOP from satellite reflectance
IOPs (absorption coefficient a and backscattering coefficient b) at various bands were derived by a semi-analytical approach from the BOA reflectance (Rrs) of S2L1C imagery atmospherically corrected by the three processors. Three separate models, one for each ACP, were developed by inputting the Rrs obtained from each processor. The basic steps of QAA version 6 (Lee et al. 2002) were chosen for the model. The above-surface reflectance (Rrs) was converted to below-surface reflectance (rrs) (Equation (2)). rrs has a robust non-linear relationship with single-scattering albedo, u(λ) represented by Equation (3) (Gordon 1988). The constants go and g1 used in Equation (3) were that for inland water bodies (Lee et al. 2002).
The next step was the retrieval of a(λo), where λo is the reference wavelength at which the contribution of absorption coefficient of water (aw) is greater than that of the other constituents present in water. 560 nm was used as λo for oceans, but longer wavelengths of 665,740 or 778 nm were found to be suitable for inland and coastal waters (Watanabe et al. 2016). Lee et al. (2002) considered 670 nm as λo if Rrs(670) > 0.0015 sr−1. In this study, all the red bands, B4 (665 nm), B5 (704 nm) and B6 (740 nm) were tested as λo in the model. a(λo) was derived by using the empirical relation proposed in QAA-v6 (Equation (5)).

Indices for trophic status
Development of model
The methodology adopted in this study for the development of the model is shown in Figure 2 and summarized below.
S2L1C imagery of the study area (Milford Lake) was downloaded and pre-processed. It was atmospherically corrected by three ACPs, namely C2RCC, Acolite and Sen2Cor.
Independent models were developed for each ACP. The output (Rrs) from each ACP was used as the input in the algorithm. QAA-v6 is the basic algorithm used. This algorithm was parametrized for S2L1C imagery (QAA-S2) based on the correlation of derived IOPs (a, b and u) and their indices with TSI.
- The correlation analysis was carried out to study the relationship between the IOPs (a, b and u) derived from the corrected S2L1C imagery of Milford Lake and TSI. The correlation of the indices, a/b, a/(a + b), ratios (BR) of a at various bands (aBR), ratios of b at various bands (bBR), the normalized difference (ND) of a between various bands (aND) and normalized difference of b between various bands (bND) were also studied. Normally, ND and BR indices enhance the characteristic spectral variations in an imagery and are hence used in this study. λo and η were parameterized based on the correlation of the indices and TSI. TSI derived from the turbidity readings was used for this study. The Pearson correlation coefficient (r) given by Equation (10) was used to assess the degree of correlation.
where x and y are variables, whose correlation is studied.
The BOA reflectance spectrum obtained from each ACP was compared with the simulated spectrum derived from the in situ spectral measurements to assess the efficiency of each ACP. The ACP whose spectrum was closest with the simulated spectrum was identified.
The regression relations with r2 > 0.7 were validated with data from validation sites.
- The accuracy of the derived relations was measured by mean absolute error (MAE) (Equation (11)) and mean absolute percentage error (MAPE) (Equation (12)). Indices suitable for estimating TSI were identified based on MAE and MAPE.where ym and yp are the measured and estimated variables, respectively.
RESULTS AND DISCUSSION
The primary step for using any satellite data is to remove the atmospheric interference. Three ACPs were used in the present study to obtain the BOA reflectance (Rrs) from TOA reflectance data of S2L1C imagery. The basic algorithm used for the development of the model was QAA-v6 which was parametrized by modifying the reference wavelength, λo, and the exponent, η, in the power relation of backscattering coefficients. Lee et al. (2002) recommend 670 nm for λo, if Rrs(670) > 0.15 sr−1. In the present study, the bands B4, B5 and B6 with central wavelengths at 665, 704 and 740 nm, respectively, were tested as λo. Shifting λo to B5 and B6 did not show significant variation in the performance of the model. Hence, the B4 band was fixed as λo. For the computation of b(λ), ratios of rrs at various bands were attempted in the equation of η (Equation (7)). η is represented in the body of the paper as and in the figures as
, where λ1 and λ2 are the central wavelengths of the bands whose rrs is considered for computing η. The parameterized algorithm is named as QAA-S2.
Rrs obtained from the three ACPs were used as inputs in the model and their performance was assessed by the correlation of the derived IOPs and their indices with TSI. For the sake of simplicity of representation, the models are referred to with the name of the ACP whose data are used as input in the model, namely C2RCC model, Acolite model and Sen2Cor model. The correlation results for the three models with the various η functions are shown in Figure 5 (C2RCC), Figure 6 (Acolite) and Figure 7 (Sen2Cor) that are discussed in subsequent sections.
Single-scattering albedo (u)
Single-scattering albedo denotes the scattering efficiency of particulate matter to the total extinction of electromagnetic energy, which is equivalent to . However, in a semi-analytical model, it is computed from rrs (Equation (3)). A correlation of r > 0.8 for u versus TSI is observed at all the red bands (B4, B5 and B6) in the C2RCC model, and at B4 and B5 in the Acolite model (Figure 3).
The reflectance from the red bands is mainly contributed by the scattering of light by particulate matter present in the water body and is proportional to the concentration of suspended matter (Zheng & DiGiacomo 2017). This good correlation of red bands with TSI indicates the retrieval efficiency of the red bands of these ACPs. The good linear correlations between u(B5) and TSI for C2RCC and Acolite models are shown in Figure 4. The slope of u(B5) versus TSI for Acolite model is observed to be steeper than that of the C2RCC model. This implies a larger window of spread for u(B5) to represent the variations in TSI in the C2RCC model and hence, a better depiction of the fine variations of TSI.
(a–f) Correlation coefficients of different indices* (derived from the C2RCC model) with TSI. *The different indices are as follows: aBR: ratios (BR) of a at various bands, bBR: ratios of b at various bands; aND normalized difference (ND) between a at various bands and bND normalized difference (ND) between b at various bands.
(a–f) Correlation coefficients of different indices* (derived from the C2RCC model) with TSI. *The different indices are as follows: aBR: ratios (BR) of a at various bands, bBR: ratios of b at various bands; aND normalized difference (ND) between a at various bands and bND normalized difference (ND) between b at various bands.
(a–f) Correlation coefficients of different indices (derived from the Acolite model) with TSI.
(a–f) Correlation coefficients of different indices (derived from the Acolite model) with TSI.
(a–f) Correlation coefficients of different indices (derived from the Sen2Cor model) with TSI.
(a–f) Correlation coefficients of different indices (derived from the Sen2Cor model) with TSI.
C2RCC model
The correlation analysis of a and b with TSI at various bands (Figure 5(a) and 5(b)) reveals that a at green (B3) and red bands (B4, B5 and B6) are well correlated to TSI (the red bands having a negative correlation). The correlation coefficient of b is around 0.8 for all the bands and for all values of η tested in this study. However, no such good correlations are observed for a/b or for all bands with TSI.
Lee et al. (2002) recommended a reflectance ratio at 443 and 560 nm for η in QAA-v6. However, a poor correlation is observed when η(443/560) is used in the C2RCC model. Among the indices developed from a, namely aBR24, aBR25, aBR26, aBR35, aBR36, aBR45 (Figure 5(c)) and aND35, aND36, aND45 (Figure 5(d)) show a correlation (r > 0.8) with TSI. In the case of backscattering coefficient, the correlations of all BR and ND indices with TSI (Figure 5(e) and 5(f)) are observed to be good (r > 0.85) when η is calculated with red bands. The C2RCC model gives the best correlation (r > 0.9) for all ND and BR indices except for bBR12 and bND12 when η(665/740) is used.
Acolite model
When Rrs values obtained from the Acolite processor were used as input in the model, the correlations of a and b with TSI are similar to that of the C2RCC model (Figure 6(a) and 6(b)).
The correlations of a with TSI at the green and red bands are good (r > 0.7) and that of b is at approximately 0.8 for all the bands irrespective of the tested η values. However, the correlation of the indices with TSI is not as good as the C2RCC model. A correlation above 0.7 is visible only for aBR25,aBR35, aND25 and aND35 when η(665/704) and η(665/740) are used in the model (Figure 6(c) and 6(d)). Both ND and BR indices of b fail to give any good correlation with TSI (r < 0.5) (Figure 6(e) and 6(f)).
Sen2Cor
When Rrs from Sen2Cor processor is used as input in the model, the correlation of a with TSI is above 0.6 only for B2 and B3 bands (Figure 7(a)), whereas for the red bands the correlation of a is observed to be poor. However, the correlation of b at all the red bands is found to be good (r ≅ 0.8) (Figure 7(b)). When aBR is tested for correlation with TSI, only those with B2 and B3 show r > 0.8. The correlations of all aND,bBR and bND indices are found to be poor (Figure 7(d)–7(f)).
The poor correlation of all the indices developed from the Sen2Cor model may be an indication of the inaccuracy in the corrected BOA spectrum. Sen2Cor is a land-based ACP which uses DDV for identifying the aerosol contribution to the reflectance sensed by the sensor. A lack of sufficient dark vegetation pixels in the scene for estimation of aerosol thickness may have produced the errors in the BOA reflectance spectrum for the study area. Studies on the Sen2Cor performance for water bodies have not shown very dependable results (Pereira-Sandoval et al. 2019; Warren et al. 2019). Now, a question arises as to why this ACP was included in the comparison. The Sen2Cor algorithm is tested here because the Level 2A products of Sentinel 2 (S2L2A) are the S2L1C imagery atmospherically corrected by the Sen2Cor algorithm.
The inferences from the correlation analysis of the three models are summarized in the following.
One of the objectives of this study is the parametrization of η used for deriving b. Lee et al. (2002) recommended the ratio in Equation (7), but this study did not give good results with the above ratio. Instead, the use of reflectance at the red bands gave better results. Correlations of IOP and their indices with TSI is best when η is derived using the reflectance ratio
for both C2RCC and Acolite models. Pereira-Sandoval et al. (2019) had evaluated various ACPs for S2 imagery and reported the inaccuracy of all the algorithms in retrieving the B1 band (443 nm). This limitation of the ACPs may be the reason for the poor correlation of the models when the B1 band is used for deriving η. Hence, this study recommends
for deriving η in this semi-analytical model (QAA-S2) that uses S2L1C data.
The observations from this study also establish that the ACP is as important as the algorithm for the accurate performance of any retrieval model based on the satellite imagery. The results from the C2RCC model are found to be the best in this study, with bND and bBR indices displaying the potential to be considered as tools for estimating TSI from S2L1C imagery. The efficiency of C2RCC has been validated in various studies (Sherjah & Sajikumar 2021). The Acolite model shows good correlation with TSI only for aBR index with B2 and B3 with B5 bands. However, the correlation of u(B5) with TSI is fairly good for both the models to be used as an estimator for TSI. The indices from the Sen2Cor model fail to show any good correlation with TSI.
Spectrum comparison
The inconsistency in the behaviour of the three models when QAA-S2 is used can only be explained by the difference in the efficiency of the ACPs in retrieving Rrs from TOA reflectance. To evaluate the performance of each ACP, the simulated spectrum derived from the in situ measurements taken at Kochi Kayal and Punnamada Kayal were compared to the BOA spectrum of S2L1C imagery obtained from each ACP for these regions. In situ spectrum measured from the two areas are shown in Figures 8 and 9.
S2L1C imagery of 26 and 28 January 2021 (in situ measurements done on 26 January 2021) had to be merged as the Kochi Kayal area lies in two tiles of imagery. S2L1C imagery for Punnamada Kayal was available for 14 March 2021 (in situ measurements were done on 13 March 2021). Atmospherically corrected spectrums from the three ACPs, the simulated S2 spectrum and the TOA spectrum for Kochi Kayal and Punnamada Kayal are shown in Figures 10 and 11, respectively.
The BOA spectrum obtained from the C2RCC and the Acolite processors match with the simulated spectrum from Kochi Kayal (Figure 10(a)–10(d)) even though Rrs obtained from the corrected imagery is lower than the in situ values. Rrs for all the bands in the Sen2Cor corrected spectrum is higher than that of the simulated spectrum (Figure 10(c)). The typical peak at B3 is not observed in the Sen2Cor corrected spectrum.
Atmospherically corrected spectrum from three processors (a–c), simulated spectrum (d) and TOA spectrum (e) – Kochi Kayal.
Atmospherically corrected spectrum from three processors (a–c), simulated spectrum (d) and TOA spectrum (e) – Kochi Kayal.
The BOA spectrum obtained from the C2RCC and the Acolite processors for some of the sample points on Punnamada Lake showed deviations from the simulated spectrum (Figure 11(a), 11(b) and 11(d)). In situ spectrum measured from the Kochi and Punnamada areas (Figures 8 and 9) did not vary much, so the variations between the simulated and the BOA spectrum can be explained only by the difference in the TOA reflectance spectrum captured from the two areas (Figures 10(e) and 11(e)).
Atmospherically corrected spectrum from three processors (a–c), simulated spectrum (d) and TOA spectrum (e) – Punnamada Kayal.
Atmospherically corrected spectrum from three processors (a–c), simulated spectrum (d) and TOA spectrum (e) – Punnamada Kayal.
High TOA reflectance is observed at B8A and B9 (water vapour retrieval bands) and also at B11 (SWIR) for the points which show erroneous values for the BOA spectrum from Punnamada Kayal. These sampling points are located on narrow streams. The intended purpose of the B11 band (1,610 nm) is snow and cloud separation and it is sensitive to lignin, starch and forest (ESA Standard Document 2015). The adjacency effect of the neighbouring land pixels may be the cause of this high TOA reflectance and it may also be the reason for errors in the corrected BOA spectrum (Pereira-Sandoval et al. 2019). It is observed that the C2RCC processor shows errors in the corrected spectrum when TOA reflectance at B8A goes above 0.04 sr−1 and that at B11 is above 0.01 sr−1. The C2RCC corrected spectrum obtained for such pixels display low Rrs values at B3, B4 and B5 (Figure 11(a)). The peak of the BOA spectrum shifts from B3 to B2 for those points.
The corrected spectrum from the Acolite shows very low Rrs values for all points (10 times less) (Figure 11(b)). These erroneous Rrs values would lead to inaccuracy in the results from any retrieval algorithm from the imagery. The BOA spectrums from the Sen2Cor processor for both Kochi and Punnamada areas (Figures 10(c) and 11(c)) do not match with the simulated spectrum. The deviations in the Sen2Cor corrected spectrum from the simulated spectrum explain the weak correlation of IOPs and their indices (derived from Sen2Cor data) with TSI. S2L2A are level 2 products of S2 that are atmospherically corrected by the Sen2Cor algorithm and caution should be exercised while using them for any water body analysis.
Validation
The correlation analysis of the three models with the data from Milford Lake has helped us to identify the indices which have the potential for estimating TSI. The indices which have a correlation of r2 > 0.7 with TSI were used to estimate TSI from the validation sites (Table 2 for the C2RCC model and Table 3 for the Acolite model).
Indices and their regression equations for the C2RCC model
aBR | BR15 | BR16 | BR24 | BR25 | BR36 | BR45 | |||||||
R2 | 0.77 | 0.78 | 0.82 | 0.79 | 0.78 | 0.78 | |||||||
Slope | 3.95 | 10.27 | 8.24 | 6.69 | 6.69 | 68.32 | |||||||
Intercept | 51.12 | 52.83 | 48.81 | 52.18 | 52.18 | 2.41 | |||||||
aND | ND35 | ND36 | ND45 | u(B5) | u(B6) | ||||||||
R2 | 0.73 | 0.78 | 0.75 | 0.71 | 0.74 | ||||||||
Slope | 30.39 | 41.22 | 126.83 | 41.16 | 76.59 | ||||||||
Intercept | 71.37 | 90.55 | 71.11 | 56.49 | 58.39 | ||||||||
bBR | BR13 | BR14 | BR15 | BR23 | BR24 | BR25 | BR26 | BR34 | BR35 | BR36 | BR45 | BR46 | BR56 |
R2 | 0.86 | 0.85 | 0.85 | 0.71 | 0.80 | 0.81 | 0.81 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 |
Slope | −253.3 | −104.8 | −82.1 | −367.8 | −143.0 | −110.4 | −90.3 | −375.4 | −251.9 | −190.1 | −1376 | −684.4 | −1673 |
Intercept | 461.5 | 293.1 | 264.9 | 542.7 | 323.2 | 288.1 | 265.1 | 586.9 | 456.7 | 389.3 | 1604.2 | 904.9 | 1900.8 |
bND | ND13 | ND14 | ND15 | ND16 | ND24 | ND25 | ND26 | ND34 | ND35 | ND36 | ND45 | ND46 | ND56 |
R2 | 0.86 | 0.86 | 0.86 | 0.86 | 0.81 | 0.82 | 0.82 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 |
Slope | −819.4 | −511.8 | −465.9 | −437.5 | −555.9 | −494.4 | −456.2 | −1060 | −809.2 | −682.2 | −3077 | −1685 | −3668 |
Intercept | 245.5 | 254.4 | 259.4 | 264.0 | 224.5 | 231.7 | 237.3 | 238.0 | 240.6 | 242.9 | 236.9 | 236.9 | 235.2 |
aBR | BR15 | BR16 | BR24 | BR25 | BR36 | BR45 | |||||||
R2 | 0.77 | 0.78 | 0.82 | 0.79 | 0.78 | 0.78 | |||||||
Slope | 3.95 | 10.27 | 8.24 | 6.69 | 6.69 | 68.32 | |||||||
Intercept | 51.12 | 52.83 | 48.81 | 52.18 | 52.18 | 2.41 | |||||||
aND | ND35 | ND36 | ND45 | u(B5) | u(B6) | ||||||||
R2 | 0.73 | 0.78 | 0.75 | 0.71 | 0.74 | ||||||||
Slope | 30.39 | 41.22 | 126.83 | 41.16 | 76.59 | ||||||||
Intercept | 71.37 | 90.55 | 71.11 | 56.49 | 58.39 | ||||||||
bBR | BR13 | BR14 | BR15 | BR23 | BR24 | BR25 | BR26 | BR34 | BR35 | BR36 | BR45 | BR46 | BR56 |
R2 | 0.86 | 0.85 | 0.85 | 0.71 | 0.80 | 0.81 | 0.81 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 |
Slope | −253.3 | −104.8 | −82.1 | −367.8 | −143.0 | −110.4 | −90.3 | −375.4 | −251.9 | −190.1 | −1376 | −684.4 | −1673 |
Intercept | 461.5 | 293.1 | 264.9 | 542.7 | 323.2 | 288.1 | 265.1 | 586.9 | 456.7 | 389.3 | 1604.2 | 904.9 | 1900.8 |
bND | ND13 | ND14 | ND15 | ND16 | ND24 | ND25 | ND26 | ND34 | ND35 | ND36 | ND45 | ND46 | ND56 |
R2 | 0.86 | 0.86 | 0.86 | 0.86 | 0.81 | 0.82 | 0.82 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 | 0.85 |
Slope | −819.4 | −511.8 | −465.9 | −437.5 | −555.9 | −494.4 | −456.2 | −1060 | −809.2 | −682.2 | −3077 | −1685 | −3668 |
Intercept | 245.5 | 254.4 | 259.4 | 264.0 | 224.5 | 231.7 | 237.3 | 238.0 | 240.6 | 242.9 | 236.9 | 236.9 | 235.2 |
Indices and their regression equations for the Acolite model
a . | BR25 . | BR35 . | ND25 . | ND35 . | u (B5) . |
---|---|---|---|---|---|
R2 | 0.81 | 0.78 | 0.76 | 0.74 | 0.74 |
Slope | 14.18 | 31.50 | 76.07 | 74.99 | 128.71 |
Intercept | 34.42 | 30.55 | 38.61 | 62.57 | 41.62 |
a . | BR25 . | BR35 . | ND25 . | ND35 . | u (B5) . |
---|---|---|---|---|---|
R2 | 0.81 | 0.78 | 0.76 | 0.74 | 0.74 |
Slope | 14.18 | 31.50 | 76.07 | 74.99 | 128.71 |
Intercept | 34.42 | 30.55 | 38.61 | 62.57 | 41.62 |
‘Measured TSI’ was computed from the turbidity readings of water samples collected from the sample locations. ‘Estimated TSI’ was computed using QAA-S2 derived indices from the corrected S2L1C imagery of the respective regions. The validation result for each lake is shown separately to comprehend the results and figure out the possible cause for nonconformity. The MAPEs for the estimated TSI for the validation sites are shown in Figures 12 and 13 for the C2RCC model and Figure 15 for the Acolite model.
(a–d) Mean absolute percentage error for TSI estimation using indices developed from the C2RCC model.
(a–d) Mean absolute percentage error for TSI estimation using indices developed from the C2RCC model.
Validation results of the C2RCC model
TSI could be estimated from Kochi Kayal with a reasonable accuracy using the aND and aBR indices (MAE < 8; MAPE ≅ 9%) and using the bND and bBR indices with an MAE of 5 (MAPE ≅ 7.5%). TSI estimation from Punnamada could not achieve the same accuracy because of inaccuracy in the corrected spectrum from C2RCC. It is observed that deviations from actual values were high for the sample points on narrow streams. When the sample points from the narrow streams were excluded, the MAE and MAPE were reduced to 3 and 6%, respectively. This indicates that the cause of error in the estimation of TSI is due to the adjacency effect of the neighbouring land pixels.
TSI estimation from the Conroe Lake is the most accurate (MAE ≅ 2; MAPE < 2%). The MAPE of around 8% (MAE < 5) is observed for the Texoma Lake and around 13% (MAE ≅ 9) for the Canyon Lake. It may be noted that the Canyon Lake is a small lake of surface area of 1.55 km2 only. The adjacency of land may have influenced the radiance sensed from the water body by the satellite. The error of TSI estimation for the Punnamada Lake and the Canyon Lake is the highest when the TSI is estimated with aND45 and aBR45 (Figure 12(a) and 12(b)), and with bBR23 and bND23 (Figure 12(c) and 12(d)) indicating that the error in the corrected spectrum is greater in these bands. From Figure 11(a) and 11(d), it is observed that Rrs at B3, B4 and B5 bands of the C2RCC retrieved spectrum are lesser than that of the simulated spectrum. Hence, it can be inferred that the lesser accuracy of estimation from some of these lakes is because of the inability of the ACP to consider the adjacency effect of the land pixels. The accuracy of estimation with the backscattering coefficients is better in the C2RCC model (MAE ≅ 5; MAPE < 7%) (Figure 12(c) and 12(d)).
To evaluate the influence of parameterization of η on the accuracy, the TSI is estimated with the regression relations derived from the original QAA-v6. The MAPE for the estimated TSI (using bND and bBR indices) reduces when QAA-S2 is used instead of QAA-v6 (Table 4). The higher MAPE in TSI estimation for Punnamada Kayal is also reduced when QAA-S2 is used, though it is not as good as for the other lakes, owing to the reasons stated earlier. These results confirm the improved performance of the parametrized algorithm.
Comparison of the MAPE for the TSI estimated using different indices
Sl. No . | Site . | BND and bBR (QAA-v6) (%) . | bND and bBR (QAA-S2) (%) . | u (%) . |
---|---|---|---|---|
1 | Kochi Kayal | 9 | 7.5 | 6.7 |
2 | Punnamada Kayal | 68 | 34 | 9.5 |
3 | Conroe Lake | 7.8 | 1.7 | 3.2 |
4 | Texoma Lake | 13 | 7.7 | 10.6 |
5 | Canyon Lake | 42.5 | 13 | 12 |
Average | 28 | 13.5 | 8.4 |
Sl. No . | Site . | BND and bBR (QAA-v6) (%) . | bND and bBR (QAA-S2) (%) . | u (%) . |
---|---|---|---|---|
1 | Kochi Kayal | 9 | 7.5 | 6.7 |
2 | Punnamada Kayal | 68 | 34 | 9.5 |
3 | Conroe Lake | 7.8 | 1.7 | 3.2 |
4 | Texoma Lake | 13 | 7.7 | 10.6 |
5 | Canyon Lake | 42.5 | 13 | 12 |
Average | 28 | 13.5 | 8.4 |
TSI was estimated using u(B5) and u(B6) indices with an average MAE of 5 (MAPE < 9%) from all the sites (Figure 13 and Table 4). The reasonable accuracy of estimation using u(B5) is depicted in Figure 14, which displays the measured and estimated TSI plotted against a 1:1 line. u(B5) and u(B6) could estimate TSI with good accuracy (MAPE <12%) from all the lakes. Though the accuracy of estimation is better with bND and bBR indices derived from QAA-S2, u is computationally simpler than bND and bBR indices. Hence, u can be marked as the computationally efficient estimators of TSI.
Validation results of the Acolite model
The Acolite model gave good correlations with TSI only for aBR and aND indices with B2, B3 and B5 bands. However, the accuracy of estimation of TSI was not as good as the C2RCC model (Figure 15).
The u(B5) index is a better estimator for TSI using the Acolite model also. TSI estimation from Punnamada Lake is also not good using the Acolite-derived indices. The BOA reflectance spectrum for the sample points could not be obtained for the Texoma Lake and the Canyon Lake from the Acolite processor. The Acolite processor extracts the water pixels from the imagery by considering a threshold reflectance at 1,600 nm (B11) as 0.0215. Most of the sample points from these two lakes were excluded from water pixels due to high reflectance at the B11 band.
The better accuracy of TSI estimation using the C2RCC-derived indices establishes the C2RCC processor as a better ACP for S2 imagery containing inland water bodies. The relations developed from u(B5), u(B6), bND and bBR have demonstrated their effectiveness and transferability to different areas of the globe with reasonable accuracy. The u(B5) and u(B6) are identified as simple and effective indices for estimating TSI. These indices can be utilized to obtain a synoptic view of the TSI of the area of study from their S2 imagery. The synoptic view of spatial and temporal variations in TSI would help in identifying the critical areas of poor water quality and make their regular monitoring simpler and economically feasible.
CONCLUSION
The real-time monitoring of inland water bodies is possible only if the information about their eutrophic status can be obtained from satellite data. This study proposes a semi-analytical model for estimating the eutrophic status of an inland water body (based on the TSI) from S2L1C imagery of the area.
The exponent η (power relation for backscattering coefficients) in QAA-v6 is parameterized by modifying wavelengths used in the ‘below-surface reflectance’ ratio. The parametrized model, QAA-S2, has increased the accuracy of estimation by almost 50% for most of the lakes and hence, it can be considered as an appropriate tool for TSI estimation from S2 imagery.
Based on the comparative analysis of three ACPs, namely C2RCC, Acolite and Sen2Cor, the C2RCC is recommended for QAA-S2.
Single-scattering albedo (u) at the red-edge bands of S2L1C are identified as computationally efficient estimators for TSI from their ability to estimate TSI within a MAPE of 12% for all the lakes.
Normalized Difference (ND) and ratios (BR) of backscattering coefficient (b) could also estimate the TSI with good accuracy (MAE < 5 and MAPE < 8%) except from the regions with close proximity to land.
These indices can be utilized to estimate TSI from S2L1C imagery to get a synoptic view of water quality of an inland water body. The spatial and temporal synoptic views of variations in TSI can make the identification of the critical areas with poor water quality easier and make their regular monitoring economically viable.
TSI estimation from the pixels with close proximity to land shows deviations from the measured values. The adjacency effect of the neighbouring land pixels is found to have a strong influence on the retrieval accuracy of the ACPs, which in turn impacts the TSI estimation from these pixels. ACPs with the ability to detect and correct the adjacency effect of neighbouring pixels may be tried in the model to improve the accuracy for pixels close to land.
Future scope for research
This study used IOP derived semi-analytically from corrected S2 reflectance data for estimating TSI. However, these derived IOPs were not validated with their in situ measurements.
This model could be further improved to estimate the TSI based on Chlorophyll-a if blue bands could be retrieved with more accuracy and also 412 nm band could be included in the terrestrial satellite missions. More research may be dedicated to explore the utilization of backscattering coefficients and the red-edge bands for water quality monitoring.
DATA AVAILABILITY STATEMENT
All relevant data are available from an online repository or repositories.