Regional Models for High-Resolution Retrieval of Chlorophyll a and TSM Concentrations in the Gorky Reservoir by Sentinel-2 Imagery

Abstract

The possibilities of chlorophyll a (Chl a) and total suspended matter (TSM) retrieval using Sentinel-2/MSI imagery and in situ measurements in the Gorky Reservoir are investigated. This water body is an inland freshwater ecosystem within the territory of the Russian Federation. During the algal bloom period, the optical properties of water are extremely heterogeneous and vary on scales of tens of meters. Additionally, they vary in time under the influence of currents and wind forcing. In this case, the usage of the traditional station-based sampling to describe the state of the reservoir may be uninformative and not rational. Therefore, we proposed an original approach based on simultaneous in situ measurements of the remote sensing reflectance by a single radiometer and the concentration of water constituents by an ultraviolet fluorescence LiDAR from a high-speed gliding motorboat. This approach provided fast data collection including 4087 synchronized LiDAR and radiometric measurements with high spatial resolutions of 8 m for two hours. A part of the dataset was coincided with Sentinel-2 overpass and used for the development of regional algorithms for the retrieval of Chl a and TSM concentrations. For inland waters of the Russian Federation, such research was performed for the first time. The proposed algorithms can be used for regular environmental monitoring of the Gorky Reservoir using ship measurements or Sentinel-2 images. Additionally, they can be adapted for neighboring reservoirs, for example, for other seven reservoirs on the Volga River. Moreover, the proposed ship measurement approach can be useful in the practice of limnological monitoring of inland freshwater ecosystems with high spatiotemporal variability of the optical properties.

1. Introduction

The largest cities, industrial centers, and agricultural lands are often located on the banks of inland waters. Their overall prolonged impact on the aquatic environment and specific hydrological characteristics of inland waters lead to significant declining of water quality [1,2]. Self-purification mechanisms existing in the seas and oceans are weak [3] in inland waters. Due to the importance of freshwater bodies for human activities, much attention is paid to the monitoring of their quality. For these reasons, regional bio-optical algorithms for retrieval of water quality parameters from satellite images [4,5,6] are developed for many lakes and reservoirs all over the world. In the Russian Federation, unlike Europe, Asia or America, such algorithms were developed for significantly fewer water bodies: the Ivankovo [7], the Krasnoyarsk [8], and the Gorky [9] Reservoirs, and the Lakes Baikal [10,11] and Ladoga [12]. These studies were limited by retrieval of Chl a concentration by regression to Level 2 (L2) satellite data. The radiometric measurements and atmospheric correction were not performed. Meanwhile, the amount of the Russian inland waters which are important for ensuring human activities and protection of unique flora and fauna is significant. Among them there a great deal of highly productive reservoirs with different sizes from tens to hundreds of kilometers. The presence of regular channel currents and variable wind forcing leads to significant spatial-temporal variability of their water optical properties on scales of tens of meters and some hours. These factors limit the usage of the MODIS and Landsat-7,8 data due to insufficient resolution and low revisit interval (16 days), respectively.

MSI (MultiSpectral Instrument) optical sensors installed on Sentinel-2 satellites are a promising tool for studying of inland freshwater ecosystems [13]. Despite the fact that the main purpose of Sentinel-2 is global land monitoring, there are some studies demonstrating the possibility of its usage for assessment of inland water quality. In particular, in [14] the signal-to-noise ratio of MSI bands relative to water bodies was estimated, the atmospheric correction (AC) in SeaWiFS Data Analysis System (SeaDAS) was carried out, its accuracy according to AERONET-OC was estimated, and vicarious calibration was performed. In [15], different AC methods for MSI images by in situ radiometric measurements in four Amazon lakes were evaluated. Authors showed significant accuracy variations for all spectral bands depending on the water optical properties of studied lakes. Additionally, the water reflectance above 705 nm is sensitive to the adjacency effect. Similar results were obtained for Lake Starnberg, Germany [16].

The possibility of MSI bands usage for constituent estimation and the development of appropriate bio-optical algorithms were studied in [17,18,19,20,21,22,23]. The accuracy of the Chl a retrieval algorithm for eutrophic waters developed for MODIS, MERIS, and other ocean color sensors was evaluated using MSI data by Grendaitė [17]. In [18], calibration and validation of well-proven algorithms [24,25] were performed by in situ measurements in the Barra Bonita Reservoir, Brazil. The suitability of the commonly used green to red band ratio was tested for estimation of the colored dissolved organic matter (CDOM), dissolved organic carbon (DOC), and water color in Estonian lakes with different trophic state indices in [19]. Empirical relation with CDOM absorption coefficient was obtained by the usage of the model radiance spectrum [20]. At the same time, Kutser showed that information on Chl a and TSM concentrations in CDOM-rich lakes can be obtained from peaks at 710 and 810 nm using B5 (705 nm) and B7 (783 nm) bands [21]. In [22], the single-band algorithm for TSM concentration retrieval was developed for Lake Poyang, China. In [23], the TSM retrieval algorithm [26] was successfully applied to B4 band and the consistency between Sentinel-2A/2B and Landsat-8 data was shown.

The most part of algorithms for the processing of the MSI data is developed on the basis of in situ measurements for inland waters with smooth heterogeneity of optical properties. The last one strongly differs for the waters of the Gorky Reservoir during the algal bloom period. The spatial patterns of Chl a concentration vary in time so rapidly that satellite image may not correspond to in situ data after a few hours or faster. Another problem is the regular presence of cumulus and cirrocumulus clouds. In this case, some stations for sampling can appear under clouds or its shadow in a satellite image. Therefore, the usage of the traditional station-based measurements may not be rational. On the one hand, it will require a large number of expeditions to obtain statistically valid results. On the other hand, it will require additional significant financial and time costs for laboratory analysis of water samples, which may be useless due to above-mentioned reasons. Thus, another approach of in situ measurements may be required in order to minimize these risks and conservation of the results reliability. One of them is proposed in the present paper. The obtained in situ data were used to investigate practical possibilities of Sentinel-2/MSI for inland waters with high spatiotemporal variability of the water optical properties on the example of the Gorky Reservoir.

2. Materials and Methods

2.1. Study Area

The Volga River is the main waterway of central Russia and the longest European river. Its current ecological state is under the close attention of environmentalists due to shoaling, pollutions and phytoplankton blooms [27,28,29]. The last factor is the most pronounced in reservoirs of the upper and middle Volga, for example, in the Gorky Reservoir. This reservoir (56.65°–58.08°N, 38.83°–43.37°E) has 427 km long and covers 1590 km² (Figure 1). Its volume is 8.71 km³. Average and maximum depths are 3.65 m and 26.6 m respectively. A shipping channel of the Volga River runs along the reservoir (Figure 1a). The last 100 km forms a lake (Figure 1b). The lake part is well suited for studying different hydrophysical processes due to the presence of unique features of hydrology, variety of winds and current velocities, intensive phytoplankton bloom [30,31,32,33,34,35]. Specific parameters that quantitatively characterize the hydrological and bio-optical regimes of the Gorky Reservoir are presented in

Table 1. Hydrological and bio-optical characteristics of the lake part of the Gorky Reservoir according to authors own measurements of 2016–2018 years [34,35] and previous studies [30,31,32,33].

Parameter	Value
Maximum depth, m	26.6
Average depth, m	3.65
Maximum current velocity, cm/s	12.0
Average current velocity, cm/s	3.0
Average wind speed, m/s	1.7
Average wind gust speed, cm/s (for a 10 min interval)	3.3
Average wind gust speed, cm/s (for a 1-min interval)	14.3
Prevailing wind direction	SW
Maximum water temperature, °C	32
Chlorophyll a concentration, mg/m3	0.5–460
TOC concentration, mg/L	9–21
TSM concentration, mg/L	5–20
Secchi depth, m	0.2–3.5
Photic zone depth, m	1.0–4.1
Trophic state index	eutrophic

.

According to [30,31,34], the intensive phytoplankton bloom (Figure 1c) is observed from June through October and basically presented by algae of green, blue-green, cryptophytes, and diatom species. During this period there is quantitative domination of the blue-green algae, represented by Aphanizomenon flos-aquae (L.) Ralfs ex Born. et Flah., Anabaena sp., Anabaena flos-aquae Bréb. ex Born. et Flah., Merismopedia tenuissima Lemm., Microcystis aeruginosa (Kütz.) Kütz., Microcystis wessenberdgii (Kom.) Kom. In calm weather, the algal spatial distribution is highly heterogeneous, as demonstrated by the Sentinal-2A image of 15 July 2016 (Figure 1d). There are various structures in this image: sharp fronts with the spatial scale of tens of meters, large-scale areas of quasi-uniform distribution, and vortex structures with increased concentrations of algae. These structures move with average speeds of 3–5 cm/s under the action of channel current [35]. This leads to algal shifting of more than 100–200 m per 1 h and more than 2.5–5 km per day. At windy weather, shifting speed may increase and reach values more than 7 cm/s. As an example, large-scale algal shifting from the right bank to the left one under west wind proceeded four days is shown by a series of MODIS images in Figure 2. According to archival meteorological data [36], west and northern winds are regular and continual for this region. Generated wind waves have maximum fetch and effect on the vertical mixing of algae. As a result, it leads to a decrease of phytoplankton concentration in the near surface water layer where the satellite signal is formed. These factors limit the applicability of traditional station-based measurements of remote sensing reflectance (R_rs) and water sampling and require another approach for in situ measurements.

2.2. Field Measurements

2.2.1. Approach Description

The necessity of a suitable approach for in situ measurements in the Gorky Reservoir was determined by the following reasons:

1)

strong heterogeneity of phytoplankton on scales from tens to thousands of meters;

2)

short lifetime of phytoplankton distribution relative to the moment of satellite overpass due to river current and wind forcing; and

3)

potential possibility to study small-scale patterns with sizes of tens of meters using Sentinel-2/MSI radiometer.

According to conditions #1 and #3, it is required to perform field measurements with resolution equal to the minimum MSI radiometer resolution or less. At the same time, with accordance to condition #2, the measurements should be performed within time interval equaled to the lifetime of phytoplankton distribution or faster. Simultaneous radiometric and water constituents’ measurements from a high-speed gliding motorboat completely satisfy these conditions. We used the single Ocean Optics USB2000+ spectrometer to measure the remote sensing reflectance and the ultraviolet fluorescent LiDAR UFL-9 [37] to assess Chl a and TSM concentrations. Both devices were installed on the bow deck (Figure 3). Motorboat length was about 9 m which, several times, exceeded the length of dominant long waves. In this case, the movement was stable: noticeable pitching and rolling were absent, and ship waves and splashes near the bow did not affect registered signals. While moving with a cruise speed of 8 m/s, we continuously registered field data with frequencies of 1 Hz for the spectrometer and 2 Hz for the LiDAR. This, the spatial data resolutions were equal to 8 m and 4 m, respectively.

In situ measurements were carried out from the south side of the reservoir (Figure 4) under Sentinel-2B and MODIS Aqua/Terra overpass on 21 and 22 September 2018, respectively (the results of MODIS images processing are not included in this study). They were performed from 8:00–9:00 UTC. The sky was clear, the weather was sunny. During this time, the Sun azimuth angles and solar elevations varied within 162.8°–180.8° and 33.3°–34.5°, respectively (Figure 4). Wind waves were smooth (height was about 0.3–0.5 m) according to WMO Sea State Code on the first day and calm on the second one. The motorboat route began at the start point, passed along four tracks and ended at the Finish point, which coincides with the start point (Figure 4). Each track was about 6–8 km and took about 10–15 min, meanwhile, the Sun position changed by azimuth and elevation for 4.5° and 0.3°, respectively. Therefore, we considered that lighting conditions change slightly. However, to control this assumption, the downwelling irradiance was measured by the spectrometer at the Start and Finish points for each track. Simultaneous measurements of the upwelling radiance and fluorescence signals were performed by the spectrometer and LiDAR, respectively, along each track (details are shown below).

The spectrometer with field-of-view (FOV) of 20° was installed on the bow railing in the center of the motorboat at the zenith angle of 30°. These conditions were chosen based on the proximity to the requirements of the NASA protocols [38] and practical feasibility. Its azimuth angles were mechanically changed to keep the angle of 90° to the Sun when motorboat changed its track. It was necessary to ensure the constancy of the observation geometry and illumination. The second optical instrument LiDAR with FOV of 1° was also installed on the bow, slightly behind the spectrometer. It was oriented at an angle of 30° to the zenith and at 45° to the motion direction. Such a position of both optical devices was necessary to perform passive optical observations and active laser sensing of the unperturbed water surface in front of motorboat and minimized the falling of splashes and sun glints to the FOV. The motorboat position was registered by onboard Chartplotter Garmin EchoMap 721.

2.2.2. Water Sample

To recalculate fluorescence LiDAR signals to Chl a, and TSM concentration, surface water samples were collected from depth of 0–30 cm using clear polyethylene bottles at 10 stations. They were delivered to the shore within 1–2 h in a refrigerator at a temperature of −4 °C. Filtration was realized through 47 mm Whatman GF/F fiberglass filters with a pore size of 0.7 μm with a low vacuum (~0.2 bar). Filtered volume was 2 L. Filters were frozen at −16 °C and stored under dark conditions for one week. In the laboratory, Chl a concentration was determined using the spectrophotometric method [39] and calculated according to the equation for mixed phytoplankton [40]. Chlorophyll was extracted in 10 ml of 90% aqua acetone solution twice during an hour. The extracts were clarified twice by centrifugation for 10 min at 8000 r/min speed. Chl a concentration was measured by SF-14 spectrophotometer (Russia), previously calibrated using pure chlorophyll (Sigma) as a standard. Despite the fact that the spectrophotometric method does not satisfy the NASA protocols [41], it is often used to retrieve the concentration of Chl a as the most accessible method providing reliable accuracy (for example [24,39,42,43]). Intercomparison of the spectrophotometric method with two others valid by NASA protocols, fluorimetric and high-performance liquid chromatography methods were performed, for example, in [44,45].

TSM concentrations were determined gravimetrically by weight following the drying of filtered samples of known volume on pre-dried and weighed GF/F filters (pore size 0.7 μm), wherein the organic and mineral suspended matter concentrations were obtained by their spectral absorption in accordance with the procedure described in [46].

2.2.3. LiDAR Measurements

Fluorescence LiDAR systems are widely used around the world for detecting bio-optical parameters of the oceans, seas, and inland waters. Various studies have been conducted on the comparison of LiDAR measurements and satellite data of MODIS, MERIS, and SeaWiFS radiometers in the open seas [47,48,49,50,51]. Especially, LiDAR systems are relevant for small inland waters where spatial patterns of water quality parameters change very fast.

Mentioned LiDAR UFL-9 has been involved in field measurements worldwide: in the Atlantic Ocean, in the Black, the Kara, the Aral, the Caspian, the Baltic, the South China, the Barents, the North, and the Mediterranean Seas, on Lakes Balaton and Issyk-Kul, in the Ikshinsky and the Gorky Reservoirs. Recently, it was used on the Lake Balaton in Hungary, whose geometric dimensions, shape and the water optical properties are similar to the Gorky Reservoir. As a result, high-quality ground-truth LiDAR data were obtained and used for calibration L2 MODIS data [52].

The high quality of the LiDAR data is reached due to its physical principles and technical characteristics. The ultraviolet fluorescence LiDAR UFL-9 analyses returned signal from dual excitation (355 and 532 nm) Nd:YAG laser pulses emitted at 2 Hz with the energy of 2 mJ. Detection is carried out consistently across 11 bands (355, 385, 404, 424, 440, 460, 499, 532, 620, 651, and 685 nm) on stations simultaneously with water sampling for the instrument calibration, and across four bands (355, 404, 440, and 685 nm) simultaneously in transect mode while the motorboat moves. Fluorescence intensities at 440 nm (CDOM) and 685 nm (Chl a) and the backscattering signal at 355 nm (TSM) are normalized to the Raman scattering at 404 nm and then calibrated using a set of laboratory-measured concentrations of CDOM, Chl a, and TSM.

Tested on a great number of water bodies, it allows measuring bio-optical properties with high accuracy for non-contact and express methods. According to [37,52] the total relative measurement error of UFL-9 is 10% for TSM and CDOM and 16% for Chl a. LiDAR signals processing and its calibration are exhaustively described in [53].

Based on the results of laboratory analysis of water samples, a good correlation between couples “fluorescence signal at 685 nm–Chl a concentration” and “backscattered laser signal–TSM concentration” was established (Figure 5):

$$\mathrm{Chl} a={0.019\mathrm{x}}_{685}^{3.109}, R^2=0.91,$$

(1)

$$\mathrm{TSM}={0.06\mathrm{x}}_{355}+5.37, R^2=0.64,$$

(2)

where ${\mathrm{x}}_{685}^{}$ and ${\mathrm{x}}_{355}$ are the LiDAR signals at 685 nm and 355 nm in Raman units, respectively, and R² is the coefficient of determination.

2.2.4. Radiometric Measurements

For receiving the water-leaving radiance from above-water measurements it is necessary to measure the total radiance L_u and the sky radiance reflected by the water surface L_r [54]. Usually, two separate spectrometers are used for registration of these components [55]. However, an approach using a single spectrometer was presented in [56]. The total upwelling radiance L_u above the water surface and upwelling radiance L_cuv above water-filled cuvette are consequently measured (Figure 6a,b). This cuvette with the sizes of 200 mm × 100 mm × 75 mm (length × width × height) is filled by water to the top. Its walls and bottom absorb 98% of the incident light. Therefore, the upwelling underwater radiance is considered to be zero, so, $L_{cuv}\approx L_r$. Measurements of L_u and L_r are performed consequently within a few minutes. After that, the water-leaving radiance L_w can be obtained as a difference $L_u-L_r$.

Measurements of the downwelling irradiance E_d are usually made by the third spectrometer [57,58]. But in our case, E_d was estimated through the radiance of the Lambertian surface which was a horizontal plaque L_p with known reflection coefficient R_p close to Spectralon reflectance standard (Figure 6c). This method is well described in the NASA protocols [38]. The irradiance measurement was performed immediately after L_r by the same spectrometer and took a similar time.

Summing up, the proposed method [56] consists in consequent measurement of L_u, L_r, and L_p immediately one by one, within a few minutes (Figure 6a–c). This short time interval lets us assume that illumination conditions do not change during measurements.

The obtained time series of L_u, L_r, and L_p were averaged and smoothed by the median filter. After that these spectra were converted into R_rs similarly to Mobley [54]:

$$R_{rs}=\frac{L_u-L_r}{E_d},$$

(3)

where $E_d=\mathsf{\pi }{L}_p/R_p$.

The proposed method was used in sea expeditions when the radiometric measurements were performed at ship stops (for example, [59]) or from a stationary oceanographic platform [56]. At the same time, this method does not require underwater measurements allowing to apply it while the ship moves.

We modified this approach for continuous ship measurements as follow. The spectral radiances L_r and L_p were measured in the spectral range of 378-760 nm with a resolution of 1 nm by the zenith angle of 30° with the help of Ocean Optics USB 2000 spectrometer with FOV of 20°. Such observation geometry goes beyond the NASA protocol [38], but not significantly. Therefore, it is often used in ship measurements (for example, [42]). These radiances were registered at the start and finish points of each track (Figure 6d). In our case, it was enough to measure L_r and L_p at one point for R_rs retrieval. However, each track took about 10–15 min on average. Therefore, in order to be sure that illumination conditions are constant, we used averaged spectra of ${\overline{L}}_r$ and ${\overline{L}}_p$ by data from two mentioned points. After that we obtained R_rs spectra using the next equation:

$$R_{rs}=\frac{L_u-{\overline{L}}_r}{{\overline{E}}_d},$$

(4)

where $E_d=\mathsf{\pi }{\overline{L}}_p/R_p$.

2.3. Empirical Models

2.3.1. Chl a Models

Blue-green ratio algorithms [60] are often used to derived Chl a concentration in the open ocean (Case 1 waters), where absorption of phytoplankton pigments mainly affect the water-leaving reflectance. However, in optically-complex water (Case 2 waters) due to the absorption of CDOM and scattering of TSM, they have a poor accuracy [61]. To retrieve Chl a concentration in turbid high-productive waters (eutrophic waters with phytoplankton domination) the NIR-red ratio algorithms are often used [24,25,61,62]. These algorithms use band combinations in the range of 665–680 and 700–710 nm, where water-leaving reflectance maximally and minimally sensitive to the second absorption maximum of phytoplankton pigments, and CDOM influence are negligible.

To develop an empirical Chl a model for the Gorky reservoir, we considered the semi-empirical NIR-red edge algorithms that are popular for optically complex waters: two-band (2B) and three-band (3B) algorithms [24], the Normalized Difference Chlorophyll Index (NDCI) algorithm [25], and the peak height algorithm (PH) [21]. Applied to the central wavelengths of the MSI spectral bands, the indices (band ratio) are as follows:

$$2B=R_{rs}(705)/R_{rs}(665),$$

(5)

$$3B=\left(R_{rs}{(665)}^{-1}-R_{rs}{(705)}^{-1}\right)R_{rs}(740),$$

(6)

$$NDCI=\left(R_{rs}(705)-R_{rs}(665)\right)/\left(R_{rs}(705)-R_{rs}(665)\right),$$

(7)

$$PH=R_{rs}(705)-\left(R_{rs}(740)+R_{rs}(665)\right)/2.$$

(8)

To relate results of usage indices from Equations (5)–(8) with constituent concentrations, linear (linear), 2nd polynomial (poly), exponential (exp), and power (power #1, #2) approximations were used:

$$\mathrm{Chl} a, \mathrm{TSM}=\mathrm{a}\cdot Index+b, \left(\mathrm{linear}\right),$$

(9)

$$\mathrm{Chl} a, \mathrm{TSM}=\mathrm{a}\cdot Inde{x}^2+\mathrm{b}\cdot Index+c, \left(\mathrm{poly}\right),$$

(10)

$$\mathrm{Chl} a, \mathrm{TSM}=\mathrm{a}\cdot \exp \{b\cdot Index\}, \left(\exp \right)$$

(11)

$$\mathrm{Chl} a, \mathrm{TSM}=\mathrm{a}\cdot Inde{x}^b, (\mathrm{power} \#1)$$

(12)

$$\mathrm{Chl} a, \mathrm{TSM}={(\mathrm{a}\cdot Index+b)}^c, (\mathrm{power} \#2)$$

(13)

where Index denotes the MSI indices from Equations (5)–(8).

2.3.2. TSM Models

TSM concentration is an important indicator of water quality, as well as Chl a, especially for reservoirs with significant river discharges, coastal erosion, and sediment resuspension. To retrieve TSM concentration in optically complex waters, satellite bands in the red and NIR ranges are often used [26,63,64,65]. The widespread single-band algorithms use features of the reflectance spectrum in the red region and near the reflectance peaks at 560 and 810 nm. Most of these algorithms were obtained for sediment-rich waters, where water constituents significantly differ from those in highly productive waters of the Gorky Reservoir. It was found that TSM concentration strongly correlates with Chl a concentration (coefficient of correlation r~0.85) in our study area. Consequently B4 band (665 nm) of Sentinel-2/MSI will not be optimal for retrieval of TSM concentration due to strong absorption of phytoplankton pigments and phycocyanin fluorescence in the red region. Therefore, various single-band and band-ratio models of TSM in the form of Equations (9)–(13) were considered. B3 and B5 bands associated with reflectance maxima at 560 and 705 nm were used for single-band models, as well as B6 (740 nm).

2.4. Satellite Data

One Sentinel-2B (L1C) MSI image was downloaded from the Sentinels Scientific Data Hub [66]. The time of satellite acquisition 09/21/2019 was 08:29 UTC and coincided with the time of field measurements (from 08:00–09:00 UTC). To obtain L2 products we considered the next software for the processing of Sentinel-2 images: SeaDAS, ACOLITE, C2RCC, iCOR, and Sen2Cor.

The SeaDAS [67] is a reliable atmosphere correction (AC) tool suitable for many ocean color sensors. The usage of the AC and vicarious calibration for Sentinel-2 is fully described in [14] and opens up the possibility to retrieve R_rs spectra with high accuracy. The AC for OLI “lite” (ACOLITE) is a processing software developed for aquatic applications of high-resolution Landsat (5/7/8) and Sentinel-2 (A/B) data [68]. It allows the retrieval of the water reflectance on the bottom of the atmosphere and was used to solve various tasks [69,70,71]. The SeNtinel Application Platform (SNAP) software package provides several tools for processing Sentinel-2 images and supports the installation of external plugins. We considered C2RCC [72], iCOR [73] and Sen2Cor [74], used for study of inland waters in [17,19].

Sentinel data were resampled on the resolution of 20 m. On the average, two of three in situ measurements were fallen in one 20 m pixel. To reduce the difference between inhomogeneous of in situ data and averaged satellite product, the match-up data were smoothed with a cosine filter.

2.5. Accuracy Assessment

Evaluating the models performance were done based upon three statistical metrics:

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{(y_i-y_i^m)}^2}},$$

(14)

$$MAPE=\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{\left|y_i-y_i^m\right|}{y_i^m}}\cdot 100,$$

(15)

$$Bias=\frac{1}{N}{\displaystyle \sum _{i=1}^N(y_i-y_i^m)}.$$

(16)

where $y_i$ and $y_i^m$ are the predicted and measured values.

3. Results

3.1. Statistics of Chl a and TSM Variations

Descriptive statistics of 21502 Chl a and TSM measurements per two days are given in

Table 2. Statistics of Chl a and TSM concentrations measured on 21 September 2018 and 22 September 2018.

Parameter	N	Min	Max	Mean	Median	STD
21 September 2018
Chl a (mg/m3)	15,545	0.45	96.94	15.57	13.25	10.88
TSM (mg/L)	15,545	5.00	26.77	5.87	5.88	1.10
22 September 2018
Chl a (mg/m3)	5957	4.26	463.42	57.41	38.77	52.48
TSM (mg/L)	5957	5.38	38.76	8.42	6.54	4.39

. Here. N is a number of measurements; Min, Max, and Mean are the minimum, maximum and average values respectively; Median is the median; STD is the standard.

Despite the fact that measurements were performed on closed transects for 2 days (Figure 4), the daily variability of Chl a and TSM concentrations was great. On the second day, the concentration of Chl a and TSM increased by factors of 4 and 1.5, respectively. In some places, the maximum of Chl a concentration reached a value of 463.4 mg/m³. These variations were associated with the changes in the water surface state. According to the meteorological archive [36], the first day of measurement was the last day of prolonged wind forcing of the western direction. In that morning, surface roughness was characterized by waves with an average height of 0.3–0.5 m (in opposition to 0.5–1.0 m in the previous days). In the evening the roughness was completely damped. As a result, in the morning of the second day of measurements, phytoplankton concentrated in the thin near-surface water layer and formed strong heterogeneous structures with different scales, starting from a few meters. These structures were observed during all day long due to calm weather.

3.2. Reflectance Spectra

A total of 6020 R_rs spectra were obtained per two days. Its examples for Chl a concentration varied from 0–100 mg/m³ are shown in Figure 7a. These spectra are typical for inland water with dominant cyanobacteria bloom [60,61,62,63]. The absorption maxima of phytoplankton and CDOM in the blue region well expressed. Two peaks at 560–570 nm and 710–720 nm are associated with scattering on suspended particles (organic and mineral) and with total absorption minimum. Local minimum at 620–630 nm and 670–680 nm are caused by absorption of phycocyanin (PC) and phytoplankton, respectively. A peak at 640–660 nm is associated with PC fluorescence.

In some points on the reservoir with Chl a > 100 mg/m³, reduction of R_rs at 740–760 nm due to water absorption was weak or absent (Figure 7b). Such spectra corresponded to floating algae or surface accumulation (surface scum, dense mats) showed in Figure 1c. Empirical models for these areas should be developed separately. Therefore, 1993 spectra approaching to the vegetation spectrum were excluded. The ratio Rrs(755)/Rrs(705) ≥ 0.9 was used as a criterion for exclusion.

3.3. In Situ Dataset

Two datasets, including 21,502 LiDAR and 6020 radiometric measurements, were used for the composition of the combined dataset.

LiDAR continuously registered fluorescent signals, including the motorboat stops between tracks, while the continuous recording of the upwelling radiance by the spectrometer was interrupted for measurements of water surface radiance and downwelling irradiance. In addition, these optical instruments had different temporal resolutions: 2 Hz for the LiDAR and 1 Hz for the spectrometer. As a result, we obtained two datasets with different dimension. Using averaging of the LiDAR data over one second, and temporal synchronization of two datasets, the combined dataset including 6020 measurements per two days with a spatial resolution of 8 m was collected. After that, 1933 values corresponding to excluded R_rs spectra were excluded, too. The resulting dataset including 4087 Chl a and TSM concentrations (

Table 3. Statistics of temporal-synchronous dataset used for the development of empirical models.

Parameter	N	Min	Max	Mean	Median	STD
Chl a (mg/m3)	4087	1.06	99.32	29.85	27.11	18.37
TSM (mg/L)	4087	5.00	19.48	6.55	5.88	1.65

), as well as R_rs spectra were used for the development of regional empirical models for Chl a and TSM estimation. For this purpose, 60% of this dataset were randomly selected for models calibration, and the remaining 40% for its validation. For the creation of Chl a and TSM retrieval algorithms by Sentinel-2 imagery, we used only a certain part of this dataset, including 671 measurements which were maximally close to the satellite overpass (within 10 min). These data corresponded to measurements at the first and half second tracks marked by an arrow line in Figure 4. As before, 60% of this dataset were randomly selected for model calibration and the remaining 40% for its validation.

To develop empirical models for estimation of Chl a and TSM concentration, and assess AC accuracy, the hyperspectral radiometric data were resampled in Sentinel-2B/MSI bands to obtain the simulated remote sensing reflectance $R_{rs}(\lambda )$ as:

$$R_{rs}(\lambda )={\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}SRF(}\lambda )R_{rs}^{meas}(\lambda )d\lambda /{\displaystyle \underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}SRF(}\lambda )d\lambda ,$$

(17)

where $SRF(\lambda )$ is the spectral response function of Sentinel-2B/MSI; $R_{rs}^{meas}(\lambda )$ are in situ hyperspectral radiometric data; and ${\lambda }_1$ and ${\lambda }_2$ are the lower and upper wavelengths of the MSI spectral band. Further, they were overlapped on the satellite image and averaged pixel-by-pixel.

3.4. Calibration and Validation of Empirical Models

3.4.1. Chl a Models

Accuracy assessment for different empirical models on the basis of the calibration dataset (N = 2452, 60%) and the validation dataset (N = 1635, 40%) are presented in

Table 4. Calibration and validation results of Chl a models. All parameters were derived from the calibration dataset (N = 2452, 60%) and the validation dataset (N = 1635, 40%).

	Calibration					Validation
Index	R2	a	b	c		RMSE (mg/m3)	MAPE (%)	Bias (mg/m3)	Fits
2B	0.71	24.463	−8.356	-		9.86	44.1	0.03	linear
	0.71	−2.155	33.034	−15.59		9.76	41.6	0.02	poly
	0.62	12.546	0.526	-		11.19	60.5	0.59	exp
	0.70	17.346	1.198	-		9.98	47.9	0.27	power #1
	0.71	96.806	−61.71	0.764		9.76	40.6	−0.05	power #2
3B	0.68	35.812	19.24	-		10.26	53.6	0.06	linear
	0.73	−10.227	52.867	16.91		9.54	46.6	0.07	poly
	0.49	25.114	0.554	-		12.72	74.4	0.77	exp
	0.74	4373.068	296.6	0.479		9.33	39.0	0.02	power #2
NDCI	0.69	90.101	13.75	-		10.13	37.1	−0.02	linear
	0.71	88.757	52.715	15.06		9.79	41.9	0.01	poly
	0.71	16.365	2.754	-		9.83	45.2	0.15	exp
	0.71	6.251	3.359	2.217		9.80	40.4	−0.01	power #2
PH	0.69	2340.123	9.185	-		10.12	47.4	−0.02	linear
	0.69	−17567.0	2814.313	7.11		10.06	45.4	−0.03	poly
	0.62	17.808	52.313	-		11.21	60.8	0.41	exp
	0.69	890.568	0.708	-		10.06	42.9	−0.13	power 1
	0.69	7834.173	8.051	0.789		10.05	44.8	0.02	power 2

. The exponential Chl a models for all indices gave the most inaccurate prediction of Chl a concentration (MAPE = 45–74% and Bias = 0.15–0.77 mg/m³). The power 2B and PH models showed the highest values of Bias (0.27 and –0.13 mg/m³). Among other models, it was difficult to select the best one because of the closeness of their efficiency parameters. The power and polynomial 2B models gave approximately the same estimates (RMSE = 9.76 mg/m³; MAPE = 40.6% and MAPE = 41.6%; Bias = 0.05 mg/m³ and Bias = 0.02 mg/m³). The linear 2B model had slightly lower accuracy (MAPE = 44.16%), but could also be used. The power 3B model [61] had better parameters (RMSE = 9.33 mg/m³; MAPE = 39.0%; Bias = 0.02 mg/m³) than the linear and polynomial 3B models due to more accurate fitting for Chl a > 80 mg/m³. In addition, this model showed the lowest RMSE among all models. The power #2 and the polynomial NDCI models had the same RMSE and MAPE as similar 2B models, however, their Bias is smaller. At the same time, the polynomial NDCI model overestimates forecasts, but the power model underestimates it. RMSE and Bias of the linear NDCI models are slightly higher but its MAPE of 37.1% is the lowest among all investigated models. For PH, the polynomial and the power #2 models gave the best results, but their efficiency parameters are worse than NIR-red edge models. The best models for each index are marked in bold in Table 4.

Despite the fact that most of the models explain about 70% of the variation of Chl a concentration (R²~0.7), and RMSE~9–10 mg/m³, all of them cannot be considered as models with high prediction accuracy (MAPE ≥ 40%). Perhaps this result is due to insufficient consistency of radiometric and LiDAR data (Figure 8a) produced by two reasons. The first one is different footprints on the water surface: 0.4 m² for the spectrometer and 0.02 m² for the LiDAR. The second one is its shifting from each other to 1.5 m. Thus, these instruments scan different parts of the water surface. As a result, the difference between radiometric and LiDAR signals (Figure 8a) could lead to a significant scattering in the synchronous dataset used for calibration and validation of models (Figure 8b,c).

Furthermore, it should be noted that the calibration coefficients for all models were significantly lower in comparison with [24,25,61,62], possibly due to features of the hydro-optical water properties of the Gorky Reservoir. However, the calibration coefficients of the linear PH model (y = 2231x + 12.7; R² = 0.80) obtained in [19] were quite close to those obtained in our research.

3.4.2. TSM Models

Calibration and validation of these models were carried out using the same dataset (Table 3). Their efficiency was evaluated according to Equations (14)–(16). The main results are shown in

Table 5. Calibration and validation results of TSM models. All parameters were derived from the calibration dataset (N = 2452) and the validation dataset (N = 1635).

	Calibration					Validation
Index	R2	a	b	c		RMSE (mg/L)	MAPE (%)	Bias (mg/L)	Fits
Single band algorithms
B3	0.69	141.339	3.949	-		0.68	7.20	−0.008	linear
	0.71	2277.23	42.359	4.85		0.66	6.49	−0.001	poly
	0.71	4.437	20.386	-		0.66	6.69	−0.004	exp
	0.64	26.494	0.348			0.75	7.86	−0.057	power #1
B5	0.71	90.735	4.646	-		0.66	6.48	−0.006	linear
	0.71	149.813	82.16	4.734		0.76	7.76	−0.376	poly
	0.70	4.98	12.356	-		0.67	6.39	0.002	exp
	0.65	18.071	0.257			0.74	7.62	−0.060	power #1
B6	0.74	122.89	5.27	-		0.63	6.21	0.006	linear
	0.75	−871.689	157.24	5.103		0.60	6.13	−0.001	poly
	0.70	5.509	15.105	-		0.70	6.89	0.016	exp
	0.64	14.838	0.170	-		0.73	7.85	−0.056	power #1
Band ratio algorithms
B3/B4	0.66	3.127	1.857	-		0.72	6.93	−0.026	linear
	0.69	1.495	−1.631	5.449		0.69	6.16	−0.018	poly
	0.68	3.213	0.466	-		0.70	6.51	−0.020	exp
B5/B4	0.68	1.683	3.751	-		0.69	6.48	−0.022	linear
	0.68	0.143	1.128	4.213		0.69	6.16	−0.019	poly
	0.68	4.341	0.239	-		0.69	6.09	−0.015	exp
B6/B5	0.62	8.032	3.119	-		0.74	8.93	−0.008	linear
	0.71	16.955	−7.345	6.306		0.64	6.78	0.014	poly
	0.67	3.79	1.251	-		0.69	8.25	−0.006	exp
PH	0.66	160.207	4.96	-		0.72	6.62	−0.016	linear
	0.66	788.31	139.449	5.048		0.71	6.45	−0.015	poly
	0.66	5.167	22.503	-		0.72	6.40	−0.011	exp

.

The single-band algorithms in Table 5 showed a high correlation (R² > 0.7). B3, B5, and B6 bands provided good efficiency. Among the models based on R_rs maximum, B3 polynomial and B5 exponential ones were more accurate. However, the B6 polynomial model (marked bold in Table 5) showed the highest coefficient of determination (R² = 0.75) and the prediction accuracy (RMSE = 0.60 mg/L; MAPE = 6.13%; Bias = −0.001 mg/L) (Figure 9a,b). Probably, this is due to the absorption by CDOM and phytoplankton pigments at 740–760 nm is negligible and water-leaving radiance is mainly determined by the absorption of pure water and scattering on suspended particles. For this reason, we think that B6-based models are more reliable.

We also considered band ratio models using bands corresponding to the minimum and maximum absorption, and the peak at 705 nm. Despite the fact that B2 and B4 single-band models showed the same accuracy of TSM retrieval, B3/B2 and B3/B4 models showed different efficiency. In particular, coefficients of correlation $r$ were 0.45 (not shown) and 0.82, respectively, which confirms the unsuitability of the blue-green band ratio algorithms for optically complex waters even in the presence of covariance between TSM and Chl a concentrations. At the same time, it was difficult to determine the best algorithm because the all band ratios presented in Table 5 have similar statistical parameters. We have marked italic one algorithm for each index. Finally, PH models were slightly worse relative to other band-ratio algorithms and, therefore, it can also be used to predict TSM concentrations (Figure 9c,d).

3.5. Atmospheric Correction

The aerosol radiance due to scattering by aerosols has been greatly overestimated when the Sentinel-2 image was processed in l2gen (SeaDAS). As a result, retrieved R_rs spectra were negative at all wavelengths. Therefore, we were unable to use this processing tool in our study. Examples of retrieved R_rs spectra obtained by ACOLITE, C2RCC, and Sen2Cor are shown in Figure 10. In our case, the C2RCC algorithm, using neural network technologies, could not predict the shape of the spectral curve. At the same time, the magnitude of R_rs was one order less than the measured one. Such difference from [5] may be related to different water optical properties of the Estonian lakes (${\rho }_w(560)=0.003-0.02$, R_rs = 0.001–0.0065 sr⁻¹) and the Gorky Reservoir (${\rho }_w$ = 0.04–0.11; R_rs = 0.013–0.035 sr⁻¹). Two other plugins, ACOLITE and Sen2Cor, provided similar R_rs spectra. Usually, they were lower than measured ones, excepting blue bands where Sen2Cor gave much higher values. The adjacency effects at 740 nm also appeared the same. The iCOR (not shown) did not reproduce the spectral shape: there was no decrease at 665 nm and R_rs values at B1–B5 bands were similar. In addition, iCOR had no pixel geolocation, which made it impossible for usage.

Based on the analysis results, we chose ACOLITE because it better reproduces the shape of the reflectance spectra. Additionally, ACOLITE is more flexible in configuration than SNAP plugins. It has the ability to apply coefficients for vicarious calibration, to choose AC algorithms (SWIR-SWIR or dark spectrum), to select bands for AC, etc.

Table 6. Atmospheric correction accuracy of ACOLITE at Sentinel-2/MSI bands based on 671 match-up points.

λ	440	490	560	665	705	740
R2	0.28	0.09	0.72	0.56	0.87	0.44
RMSE, sr−1	0.0014	0.0008	0.0041	0.0017	0.0054	0.0039
MAPE, %	33.6	9.4	17.5	12.1	20.3	71.5
Bias, sr−1	0.0013	−0.0001	−0.0030	−0.0015	−0.0037	0.0032
Slope	0.362	0.367	3.320	1.289	3.440	1.888
Intercept, sr−1	0.002	0.004	−0.027	−0.002	−0.027	−0.012
${\overline{R_{}}}_{rs,ACOLITE}/{\overline{R_{}}}_{rs,insitu}$	1.33	0.98	0.80	0.88	0.77	1.46

shows the results of AC accuracy assessment by 671 match-up point.

The retrieved reflectance at B1 and B6 bands are greater than in situ measured ones. As one can see from Table 5, the ratio ${\overline{R_{}}}_{rs,ACOLITE}/{\overline{R_{}}}_{rs,insitu}$ of the average R_rs calculated by ACOLITE and in situ data at 671 match-up points for these two bands, was 0.0013 and 0.0032 sr⁻¹, while Bias was 0.0013 и 0.0032 sr⁻¹, respectively. Perhaps, this overestimation at B6 band was associated with the adjacency effect (Figure 10a,b). When this effect did not manifest, the retrieved values are close to the measured ones (Figure 10c,d). Reflectance at B2 band almost coincides (${\overline{R_{}}}_{rs,ACOLITE}/{\overline{R_{}}}_{rs,insitu}$ = 0.98, Bias = −0.0001 sr⁻¹) which was mainly connected with uniform distribution of retrieved R_rs relatively to its true values, i.e., about half of the points were overestimated while the other half were underestimated (Figure 10). At B3-B5 bands we observed regular underestimation (Figure 10c,d) but for some pixels where R_rs (560) = 0.011–0.014 we had the opposite situation (Figure 10a,b). Generally, the obtained estimation of ${\overline{R_{}}}_{rs,ACOLITE}/{\overline{R_{}}}_{rs,insitu}$ was close to presented in [15]. Since the estimates of AC accuracy made on the basis of independent measurements coincide, we can conclude that they are not associated with measurement errors, but with atmospheric correction. At the same time, there was a significant correlation with measurements (r = 0.85, 0.75, and 0.93) so they can be used for the retrieval of Chl a and TSM concentrations.

Correlation of the band ratio was better than for individual bands (Figure 11) because AC errors were partly compensated. For 2B and NDCI indices R² were 0.93, while the slopes were 2.91 and 2.47. According to Figure 11a, the time series of 2B ratio along the ship track had common features with ratio calculated by measurements, but their magnitudes differed on 50%. Probably, this difference was caused not only by AC errors, but also a difference between in situ point measurements and area-averaged satellite data. According to Figure 8a, spatial heterogeneity of Chl a concentration on scales of tens of meters can vary by more than 50 mg/m³. In such circumstances, satellite data will have lower values.

3.6. Chl a and TSM Retrieval Models Using Sentinel-2 Data

The high correlation of R_rs at B4 and B5 bands with in situ measurements allows us to use them for satellite algorithm development. Linear and polynomial algorithms based on 2B and NDCI indices showed similar results (

Table 7. Remote sensing algorithms for Chl a retrieval and error analysis using RMSE, MAPE, and Bias on the basis of 671 match-up points. A total of 401 points were used for algorithm calibration, and the remaining 270 points for algorithm validation.

	Calibration					Validation
Index	R2	a	b	c		RMSE (mg/m3)	MAPE (%)	Bias (mg/m3)	Fits
2B	0.84	64.536	−57.800	-		3.27	13.89	0.003	linear
	0.86	−73.669	252.808	−176.68		3.01	13.15	0.031	poly
NDCI	0.86	167.293	4.756	-		3.11	13.21	−0.000	linear
	0.86	−300.260	235.556	1.586		3.02	13.14	−0.001	poly

), but we chose the polynomial NDCI algorithm as the best (Figure 12).

2B index algorithm was used for TSM retrieval (

Table 8. Remote sensing algorithms for TSM retrieval and error analysis using RMSE, MAPE, and Bias on base 671 match-up points. A total of 401 points were used for algorithm calibration, and the remaining 270 points for algorithm validation.

	Calibration					Validation
Index	R2	a	b	c		RMSE (mg/L)	MAPE (%)	Bias (mg/L)	Fits
2B	0.78	2.208	3.169			0.13	1.56	0.0003	linear
	0.78	−0.406	3.246	2.513		0.13	1.56	0.0008	poly

). Since the dependence of 2B index on TSM concentration is close to linear, there was not noticeable difference with other regression models. Therefore we chose it as the best one (Figure 13).

3.7. Chl a and TSM Mapping Using Sentinel-2 Data

The obtained algorithms were used to map Chl a and TSM concentrations (Figure 14). The white line in Figure 14a shows the track that was close to satellite overpass. Corresponding to it, in situ data were chosen for the match-up dataset and for calibration and validation of satellite algorithms. The stars show pixels corresponding to four spectra in Figure 10. In the study area, the maximum observed Chl a concentrations were about 40−50 mg/m³. At the same time, in the upstream area near left bank maximum Chl a concentrations reached 70 mg/m³. Pixels with Chl a > 70 mg/m³ were masked because ${\rho }_w(1600)$ exceeded the threshold value of 0.0215, which was used to separate water from the clouds. Areas with high Chl a concentration corresponding to phytoplankton clusters are shown in Figure 1c.

In accordance with the results of the laboratory analyses and earlier study [75], the organic suspended matter was the main component (~75%) of suspended matter in the Gorky Reservoir in that period. Therefore, TSM and Chl a maps greatly correlate with each other (Figure 14b). The maximum and average TSM concentrations were approximately 7–8 mg/L and 5.8 mg/L, respectively, and agreed with in situ measurements (Table 2). A large number of suspended constituents, likely phytoplankton Chl a, were concentrated along the left bank.

4. Discussion

4.1. Weakness of Station-Based Measurements

Basic knowledge about the relation between remote sensing reflectance and water constituents in inland waters are obtained on traditional station-based in situ measurements in conditions close to calm weather. To correct Chl a retrieval in wind-wave conditions, corrective coefficients based on wind speed and wave fetch can be applied [76]. However, it is not a solution for water bodies when the optical properties of water change rapidly due to heterogeneous currents. On the example of the Azov Sea, significant variations in fluorescence of Chl a within every 300-m and 1-km lengths along the transect were found [24]. According to the comparison with the MERIS data, the authors made two important conclusions: 1) within-pixel spatial heterogeneity of Chl a distribution resulting in the station-based in situ measurements is not representative of the satellite pixel; and 2) an actual biophysical change in the water body between the times of in situ data collection and satellite acquisition should be taken into account. As a result, they offered to account within-pixel spatial heterogeneity and temporal variation through taking multiple measurements around each station so as to characterize the spatial variation within the satellite pixel area. In the Gorky Reservoir, water property varied in smaller scales, i.e., from tens to hundreds of meters. Additionally, daily variability of TSM and Chl a was more than 150% and 400%, respectively. Thus, for inland waters, like reservoirs with channel currents, this problem is more acute. Due to these examples, it becomes obvious that traditional station-based measurements are not effective for such water bodies. Therefore, the usage of remote sensing methods with known accuracy seems more reliable.

The approach proposed in the present study has not been described in the literature before. However, its components, like radiometric and LiDAR measurements, are widely used in limnological practice. It should be noted, those continuous radiometric measurements were repeatedly used for validation of atmospheric correction. In [77], for the AERONET OC program, a consequent measurement of the sky radiance and water leaving radiance from an oceanographic tower in the Adriatic Sea were realized using a modified photometer CE-318, Cimel. Again, in [78] synchronous measurements of such parameters were performed during transects. The authors indicated the real possibility of obtaining a large number of R_rs spectra and their applicability for the validation of MODIS and VIIRS data. Such investigations have a basic goal: the changeover from manual long in situ measurements to automatic systems. To automate measurements, LiDAR systems are widely used in practical oceanology as express methods for assessing water quality [47,48,49,50,51]. They provided a high spatial resolution and a satisfied accuracy (according to some sources, like [52,79], the measurement error lies within 5–16%).

Therefore, our approach for in situ measurement seems justified. Moreover, it can be considered as a small, but significant, step towards automated remote sensing of the optical properties of water and have particular importance for productive inland water.

4.2. Chl a Retrieval Models

To develop empirical models, we considered two band (2B) and three band (3B) NIR-red algorithms, as well as models based on NDCI and peak height at 705 nm (PH), which are well suited for Chl a evaluation in eutrophic basins [5,17,18,24,25,61,62,63]. Chl a concentrations varied from 1–100 mg/m³.

Relationships between Chl a concentration and 2B and PH indices were close to linear. R², RMSE, Bias, and MAPE were similar for linear, polynomial, and power models (Table 3). On the contrary, the exponential model describing the nonlinear concentration growth from the index showed the worst result. 2B model for the Azov Sea (Chl a = 0.6–60 mg/m³) and Fremont Lake (Chl a = 2–120 mg/m³) were also linear in [24,61]. The linear relationship between Chl a concentration and PH was also obtained in [21] for the range 7–203 mg/m³, and in [19] for the range 2–120 mg/m³. Additionally, the calibration coefficients of the linear PH model (y = 2231x + 12.7; R² = 0.80) [20] were quite close to those received in our research.

Nevertheless, despite the almost linear relationship between Chl a and 2B and PH, there was a weak nonlinearity above 50 mg/m³ in our studies. Thus, the polynomial and power model y = (a × Index + b)^c were slightly better (Table 3).

The NDCI model is a nonlinear and is equally well described by both the 2nd order polynomial and the power model y = (a × Index + b)^c, which corresponds to [25]. Non-linear regression was also received for 3B index (Table 3, Figure 7b). These results mismatch [24,55], which showed that the 3B model with sufficient accuracy can be described by linear fit, and [18,25] showed that nonlinearity appears after 100 mg/m³.

4.3. TSM Retrieval Models

Most of TSM algorithms were developed for sediment-rich waters, where water constituents significantly differ from those in highly productive waters of the Gorky Reservoir. We found that TSM concentration strongly correlates with Chl a concentration (r~0.85) during cyanobacteria bloom. Therefore, we investigated various models based on B3 (centered at 560 nm) and B5 (centered at 705 nm) bands associated with reflectance maxima in eutrophic waters and B6 band (centered at 740 nm), which are minimally affected by the absorption of phytoplankton pigments.

It is known that TSM models in the red and NIR bands become strongly nonlinear at high concentrations, but could be linear in SWIR [26,65]. The analysis of regression models between TSM concentration and the simulated R_rs at Sentinel-2/MSI bands [22] also showed that nonlinearity is more pronounced at shorter wavelengths. In sediment-rich waters (17–295 mg/L) on B1–B5 bands, the TSM concentration increases exponentially, but the dependencies at B6–B8 bands become a less steep and is described by power functions with the exponent 1.22–1.38.

We also found that the regressions become more linear with increasing wavelength (Table 5). Thus, the exponential and polynomial B3-based models (R² = 0.71, RMSE = 0.66 mg/L, MAPE = 6.7%, Bias = −0.004 mg/L) are more accurate than the linear B3-based model (R² = 0.69, RMSE = 0.68 mg/L, MAPE = 7.2%, Bias = 0.008 mg/L). For B5 band, the difference between the exponential and linear models is not so obvious (Table 5). But the linear B6-based model (R² = 0.74, RMSE = 0.63 mg/L, MAPE = 6.21%, Bias = 0.006 mg/L) has a minor advantage on the exponential one (R² = 0.70, RMSE = 0.70 mg/L, MAPE = 6.89%, Bias = 0.016 mg/L). However, for B6 band at concentrations above 9–10 mg/L, nonlinearity was also significantly manifested. In this regard, the polynomial B6 model are more accurate (Table 5, Figure 9a,b).

A few words should be said about the models based on B1, B2 and B4 bands, which are not shown in Table 5, because they showed the worse results. R_rs at B1 band (centered at 440 nm) did not correlate with TSM concentration (coefficient of correlation r = 0.17) due to strong absorption of CDOM and phytoplankton pigments. R_rs at B2 (centered at 490 nm) and B4 (centered at 665 nm) bands showed a moderate correlation (r is equal to 0.7, and 0.66, respectively). At the same time, R_rs at B2 band had better correlation with TSM and provided more accurate predictions than R_rs at B4 band. Perhaps, this can be explained by the non-optimal location of B4 band due to it partially coinciding with phycocyanin PC fluorescence and the second absorption maximum of phytoplankton pigments. We found that relationships between Chl a concentration and R_rs at B4. The regressions of Chl a and TSM concentrations for the B4 band were worse than those for the B3 and B5 bands. Therefore, the usage of water reflectance in the red B4 band (centered at 665 nm), which is often used to determine the TSM concentration in sediment-rich waters, is not suitable for highly productive inland waters. Among the single-band algorithms, B6-based models are more reliable, due to the absorption by CDOM and phytoplankton pigments at 740–760 nm is negligible and water-leaving radiance is mainly determined by the absorption of pure water and scattering on suspended particles.

Analysis of the MSI/Sentinel-2 bands sensitivity to changes in water constituents allowed to determine optimal bands for assessing their concentration, as well as confirmed previous conclusions about the large potential of Sentinel-2 data for monitoring highly productive internal waters [5,16,17,18,19,20,21,22]. Empirical models developed on the basis of 4087 in situ measurements reflect the regional characteristics of the optical characteristics in the Gorky Reservoir and can further be used for rapid estimation of Chl a and TSM concentration by radiometric measurements.

4.4. Atmospheric Correction

Considered Sentinel-2 image processing tools were not allowed to retrieve remote sensing reflectance with the accuracy required for water applications. Only ACOLITE and Sen2Cor could reproduce the spectrum features at B3–B5 bands (R² = 0.56–0.87), but significantly underestimated remote sensing reflectance: average and maximum values are less by 12–23% and 100–125%, respectively (Table 6). The performance of blue bands (B1, B2) and the NIR band (B5) was considerably worse. It was noted in [68] that high performance of ACOLITE-retrieved reflectances in green (B3) and red (B4) channels mainly related to the higher reflectance range than in the blue and NIR bands. Large overestimation in NIR bands could be caused by the calibration performance at low radiances, adjacency effects, or surface effects. This explains the difference in the accuracy of the MSI bands that were also noted in [5,15], but does not explain the large difference from the in situ measurements.

We understand that these differences may be due to both measurement inaccuracies and atmospheric correction uncertainties. Therefore, we calculated the ratio ${\overline{R_{}}}_{rs,ACOLITE}/{\overline{R_{}}}_{rs,insitu}$ in 671 match-up points (Table 6) and compared them with an independent study [15]. Since the results of this comparison were almost similar, we concluded that differences between satellite-retrieved and in situ measurements are associated with the AC uncertainties. Unfortunately, no other measurements have been made in the Gorky Reservoir. The two nearest AERONET stations are 600 km and 3600 km away, so we cannot use any other data to verify this assumption. In our opinion, this difference may also be caused by a discrepancy between the point in situ measurements and area-averaged satellite data in the case of high spatial heterogeneity of water optical properties.

The atmospheric correction for turbid waters is a great challenge. It was shown in [17,19] that L1C top-of-atmosphere radiances correlate better with water quality parameters than retrieved bottom-of-atmosphere radiances (L2A products). These, and our results of AC accuracy assessment, shows that there is a strong need for further improvement of the existing algorithms of Sentinel-2 image processing for aquatic applications.

Due to the poor accuracy of the retrieved R_rs data, we cannot use them in empirical models. However, we used the B4 and B5 bands, which are well correlated with in situ measurements, to develop a remote sensing algorithm on the basis of a dataset of 671 matchups. The developed NDCI algorithm for retrieval of Chl a concentration (Figure 12) and the 2B algorithm for retrieval of TSM concentration (Figure 13) have good accuracy (R² = 0.86, RMSE = 3.02 mg/m³; MAPE = 13.14%; Bias = −0.001 mg/m³; slope = 0.864) and (R² = 0.78, RMSE = 0.13 mg/L, MAPE = 1.56%, Bias = 0.0003 mg/L, slope = 0.785) accordingly. These algorithms can be further used for mapping of Chl a and TSM concentration in the Gorky Reservoir and neighboring inland waters.

4.5. Future Plans

After completing this study, the following ways to improve measurement accuracy can be taken: (i) LiDAR calibration by a larger number of water samples [37]; (ii) the usage of a spectrometer with a collecting lens that provides measurements from a footprint similar to that of the LiDAR; and (iii) orientation of both devices at the same point of the water surface. At the same time, we realize that the usage of three synchronized spectrometers for radiometric measurements [57,58] will provide greater accuracy, guarantee clear result comparability, and eliminate the measurements of the sky radiance and the radiance of the Lambertian surface between transects. Additionally, future efforts should be focused on approaches of validation and analysis of its uncertainties. For this, the intercomparison with traditional radiometric protocols will be planned. To assess the propriety of aerosol model selection for atmospheric correction, we are going to measure the aerosol optical thickness (AOT) using a MICROTOPS portable sun photometer.

5. Conclusions

We proposed the approach of remote rapid in situ measurements with the resolution of 8 m suitable for inland waters with spatial heterogeneity and temporal variation of its optical properties. It is based on the simultaneous measurements of the remote sensing reflectance by the single spectrometer Ocean Optics USB 2000 and the concentration of water constituents by the ultraviolet fluorescent LiDAR UFL-9, both installed on a high-speed gliding motorboat. This approach was tested on the Gorky Reservoir of the Russian Federation under Sentinel-2B acquisition on 21 September 2018. As a result, the dataset including 4087 values of Chl a, TSM, and R_rs were collected from the area of 150 km² per two hours.

The regional Chl a and TSM models were developed based on 2B, 3B, NDCI and PH indices. To assess their accuracy we investigated linear, polynomial, exponential, and two power functions for approximations. All Chl a models, except the exponential one, showed approximately the same efficiency. Power models showed slightly better results due to a more accurate description of high Chl a concentrations (Chl a > 80 mg/m³). According to TSM retrieval, the single-band TSM models were more effective than the models using band ratios. The peak height at 705 nm (PH models) was also a representative parameter for Chl a and TSM prediction.

Atmospheric correction of the satellite image was realized using SeaDAS, ACOLITE, C2RCC, iCOR, and Sen2Cor processing software. All of them did not allow retrieval of Rrs spectra with sufficient accuracy. At the qualitative level ACOLITE and Sen2Cor reproduce the spectrum features at B3-B5 bands (R² = 0.56–0.87), but significantly underestimate remote sensing reflectance: average and maximum values are less by 12–23% and 100–125%, respectively. The retrieved remote sensing reflectance at B1, B2, and B6 bands are weakly correlated with measurements (R² < 0.44). The high correlation of remote sensing reflectance at B4 and B5 bands with in situ measurements allowed us to develop the algorithms to predict Chl a and TSM concentration using ACOLITE processing software and applied them for mapping in the Gorky Reservoir.

Both empirical models and remote-sensing algorithms can be applied for Chl a and TSM retrieval in next ranges: 1 < Chl a < 100 mg/m³ and 5 < TSM < 20 mg/L. These models can be adapted for neighboring reservoirs, for example, for the other seven reservoirs on the Volga River. At the same time, the proposed approach of in situ measurements can be useful in limnological monitoring of highly productive water bodies with sharp spatiotemporal variability of their optical properties.