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

: 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 inﬂuence 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 reﬂectance by a single radiometer and the concentration of water constituents by an ultraviolet ﬂuorescence 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 ﬁrst 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.


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.

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. 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  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  (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. 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]. 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 ( rs R ) and water sampling and require another approach for in situ measurements. 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].

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;

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. 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 minutes, 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 Spectrometer LiDAR Figure 3. The spectrometer and LiDAR position on a board of the high-speed gliding motorboat.
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 Remote Sens. 2019, 11, 1215 6 of 29 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.
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 29 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.

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 hour in a refrigerator at a temperature of −4 °С. 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].

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].

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): Chl a = 0.019x 3.109 where x 685 and x 355 are the LiDAR signals at 685 nm and 355 nm in Raman units, respectively, and R 2 is the coefficient of determination.

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): where 685 x and 355 x are the LiDAR signals at 685 nm and 355 nm in Raman units, respectively, and 2 R is the coefficient of determination.

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 ≈ 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. For receiving the water-leaving radiance from above-water measurements it is necessary to measure the total radiance u L and the sky radiance reflected by the water surface r L [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 u L above the water surface and upwelling radiance cuv L 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, cuv r L L ≈ . Measurements of u L and r L are performed consequently within a few minutes. After that, the water-leaving radiance w L can be obtained as a difference u r L L − .
Measurements of the downwelling irradiance d E are usually made by the third spectrometer [57,58]. But in our case, d E was estimated through the radiance of the Lambertian surface which was a horizontal plaque p L with known reflection coefficient p R close to Spectralon reflectance standard ( Figure 6c). This method is well described in the NASA protocols [38]. The irradiance measurement was performed immediately after r L by the same spectrometer and took a similar time.
Summing up, the proposed method [56] consists in consequent measurement of u L , r L , and After that these spectra were converted into rs R similarly to Mobley [54]: Figure 6. Schematic explanation of radiometric measurements: (a) the total upwelling radiance, (b) the surface-reflected radiance, and (c) the plaque radiance. Section (d) presents a photo of a water-filled cuvette at field measurements of the surface-reflected radiance.
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]: where 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 L r and L p by data from two mentioned points. After that we obtained R rs spectra using the next equation: where E d = πL p /R p .

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.

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).

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.

Accuracy Assessment
Evaluating the models performance were done based upon three statistical metrics: where y i and y m i are the predicted and measured values.

Statistics of Chl a and TSM Variations
Descriptive statistics of 21502 Chl a and TSM measurements per two days are given in Table 2. 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 3 . 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.

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 3 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 3 , 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.

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 rs R spectra were excluded, too. The resulting dataset including 4087 Chl a and TSM concentrations (Table 3), as well as rs R 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 minutes). 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 ) (λ rs R as: Figure 7. Examples of R rs spectra measured in the Gorky reservoir for two days: (a) spectra included in the temporal-synchronous dataset for calibration/validation of empirical models; and (b) spectra excluded from calibration/validation dataset.

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), 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 (λ) as: where SRF(λ) is the spectral response function of Sentinel-2B/MSI; R meas rs (λ) are in situ hyperspectral radiometric data; and λ 1 and λ 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.

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. . 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 3 ; MAPE = 39.0%; Bias = 0.02 mg/m 3 ) than the linear and polynomial 3B models due to more accurate fitting for Chl a > 80 mg/m 3 . 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. 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%). Despite the fact that most of the models explain about 70% of the variation of Chl a concentration (R 2~0 .7), and RMSE~9-10 mg/m 3 , 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 2 for the spectrometer and 0.02 m 2 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). retrieved Rrs 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 Rrs was one order less than the measured one. Such difference from [5] may be related to different water optical properties of the Estonian lakes ( . Two other plugins, ACOLITE and Sen2Cor, provided similar Rrs 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 Rrs values at B1-B5 bands were similar. In addition, iCOR had no pixel geolocation, which made it impossible for usage.  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 2 = 0.80) obtained in [19] were quite close to those obtained in our research.

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.
The single-band algorithms in Table 5 showed a high correlation (R 2 > 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 2 = 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. 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).   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).

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 (ρ w (560) = 0.003 − 0.02, R rs = 0.001-0.0065 sr −1 ) and the Gorky Reservoir (ρ w = 0.04-0.11; R rs = 0.013-0.035 sr −1 ). 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 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 in situ data at 671 match-up points for these two bands, was 0.0013 and 0.0032 sr -1 , while Bias was Figure 10. Comparison of Sentinel-2-retrieved R rs by means of atmospheric correction processors (ACOLITE, C2RCC, Sen2Cor) with in situ measurements in pixels depicted as stars in Figure 13a: (a) pixel #1, (b) pixel #2, (c) pixel #3, and (d) pixel #4. Red lines are hyperspectral measurement of R rs , blue lines are measured R rs simulated on spectral bands of Sentinel-2/MSI according [14].
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 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 R rs,ACOLITE /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 −1 , while Bias was 0.0013 и 0.0032 sr −1 , 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 (R rs,ACOLITE /R rs,insitu = 0.98, Bias = −0.0001 sr −1 ) 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 R rs,ACOLITE /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 2 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 3 . In such circumstances, satellite data will have lower values.
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 3 . In such circumstances, satellite data will have lower values.

Chl a and TSM Retrieval Models Using Sentinel-2 Data.
The high correlation of rs R 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), but we chose the polynomial NDCI algorithm as the best ( Figure 12).

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), but we chose the polynomial NDCI algorithm as the best (Figure 12). 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.  2B index algorithm was used for TSM retrieval (Table 8). 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). Table 8. Remote sensing algorithms for TSM retrieval and error analysis using RMSE, MAPE, and 2B index algorithm was used for TSM retrieval (Table 8). 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).  2B index algorithm was used for TSM retrieval (Table 8). 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).

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

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 3 . At the same time, in the upstream area near left bank maximum Chl a concentrations reached 70 mg/m 3 . Pixels with Chl a > 70 mg/m 3 were masked because ρ 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.
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.

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

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.
Relationships between Chl a concentration and 2B and PH indices were close to linear. R 2 , 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 3 ) and Fremont Lake (Chl a = 2-120 mg/m 3 ) 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 3 , and in [19] for the range 2-120 mg/m 3 . Additionally, the calibration coefficients of the linear PH model (y = 2231x + 12.7; R 2 = 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 3 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 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.
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.

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 2 = 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 R rs,ACOLITE /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

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.

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 2 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 3 ). 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 2 = 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 2 < 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 3 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.
Author Contributions: A.A.M. realized the conceptualization, methodology, in situ measurements, data analyses, as well as wrote the original manuscript. S.V.F. performed data processing, analyses, and validation, and wrote the original manuscript. They contributed towards editing and reviewing the manuscript. V.V.P. performed LiDAR measurements and data processing. E.N.K. supported the radiometric in situ measurements and data processing.
Funding: This research was funded by the Russian Science Foundation (Project RSF 17-77-10120) regarding the development of regional bio-optical algorithms and by the Ministry of Science and Education of Russian Federation (theme no. 0149-2019-0003 and Agreement 14.W03.31-0006) regarding the field ground-truth LiDAR data collection.