*Article* **Improving Water Leaving Reflectance Retrievals from ABI and AHI Data Acquired Over Case 2 Waters from Present Geostationary Weather Satellite Platforms**

#### **Bo-Cai Gao \* and Rong-Rong Li**

Remote Sensing Division, Code 7232, Naval Research Laboratory, Washington, DC 20375, USA; rong-rong.li@nrl.navy.mil

**\*** Correspondence: gao@nrl.navy.mil

Received: 2 September 2020; Accepted: 4 October 2020; Published: 7 October 2020

**Abstract:** The current generation of geostationary weather satellite instruments, such as the Advanced Baseline Imagers (ABIs) on board the US NOAA GOES 16 and 17 satellites and the Advanced Himawari Imagers (AHIs) on board the Japanese Himawari-8/9 satellites, have six channels located in the visible to shortwave IR (SWIR) spectral range. These instruments can acquire images over both land and water surfaces at spatial resolutions between 0.5 and 2 km and with a repeating cycle between 5 and 30 min depending on the mode of operation. The imaging data from these instruments have clearly demonstrated the capability in detecting sediment movements over coastal waters and major chlorophyll blooms over deeper oceans. At present, no operational ocean color data products have been produced from ABI data. Ocean color data products have been operationally generated from AHI data at the Japan Space Agency, but the spatial coverage of the products over very turbid coastal waters are sometimes lacking. In this article, we describe atmospheric correction algorithms for retrieving water leaving reflectances from ABI and AHI data using spectrum-matching techniques. In order to estimate aerosol models and optical depths, we match simultaneously the satellite-measured top of atmosphere (TOA) reflectances on the pixel by pixel basis for three channels centered near 0.86, 1.61, and 2.25 µm (or any combinations of two channels among the three channels) with theoretically simulated TOA reflectances. We demonstrate that water leaving reflectance retrievals can be made from ABI and AHI data with our algorithms over turbid case two waters. Our spectrum-matching algorithms, if implemented onto operational computing facilities, can be complimentary to present operational ocean versions of atmospheric correction algorithms that are mostly developed based on the SeaWiFS type of two-band ratio algorithm.

**Keywords:** remote sensing; sensors; ocean color; sediment; turbid water; chlorophyll; geostationary satellite

#### **1. Introduction**

The present geostationary weather satellite instruments, such as the Advanced Baseline Imagers (ABIs) on board the US GOES 16 and 17 satellites and the Advanced Himawari Imagers (AHIs) on board the Japanese Himawari-8/9 satellites, have six channels located in the visible to shortwave IR (SWIR) spectral range. These instruments can take images over land and water areas at spatial resolutions between 0.5 and 2 km and with a repeating cycle between 5 and 30 min depending on the mode of operation. The imaging data from these instruments and from other geostationary weather satellite instruments have clearly demonstrated the capability in detecting sediment movements over coastal waters and major chlorophyll blooms over deeper oceans [1–3].

During the early planning stage of GOES-R (now renamed as GOES-16), NASA, and NOAA envisioned the Hyperspectral Environment Suite (HES) instrument [4] to fulfill the future needs

and requirements of the NESDIS (National Environmental Satellite, Data, and Information Service) Office [4]. NOAA was once very interested in the coastal waters (CW) task. HES was supposed to have an instrument like an imaging spectrometer to have hourly coverages over the US east coast, west coast, Gulf of Mexico, and Hawaii coastal waters. Later on, HES was canceled by the GOES-R Project because of cost constraints, while ABI remained on GOES-R. Due to the lack of a green band on ABI, NOAA has not generated operational ocean color data products, including suspended sediments, from ABI data. On the other hand, through visual inspection of GOES-R visible band images and movies provided on a web site at the Colorado State University (CIRA), it is easy to find interesting cases of ocean color observations with ABI. Figure 1 shows such an example. Figure 1A–D are false color images (Red: 0.64 µm; Green: 0.865 µm, Blue: 0.47 µm) acquired at UTC 1332, 1532, 1732, and 1932, respectively, on September 14, 2017 over the coastal areas of Florida. The silt stirred up by Hurricane Irma was captured in the Gulf Stream and advected northward while spinning off eddies. The original movie images at a time interval of five minutes can be found from http://rammb.cira.colostate.edu/ramsdis/online/loop\_of\_the\_day/ for the day of 2017/09/14.

**Figure 1.** (**A**–**D**): time series of Advanced Baseline Imagers (ABI) images acquired at UTC 1332, 1532, 1732, and 1932 on September 14, 2017. The silt stirred up by Hurricane Irma was captured in the Gulf Stream and advected northward while spinning off eddies.

It should be pointed out that the silt features in Figure 1 are absent in images acquired under normal atmospheric wind conditions. With polar orbiting ocean color instruments, such as the aqua moderate resolution imaging spectroradiometer (MODIS) [5] and visible infrared imaging radiometer suite (VIIRS) [6], having a two-day global coverage, it was impossible to capture the rapid silt movement after the major hurricane event. Figure 1 also demonstrates partial coastal watch capabilities with ABI originally intended by NOAA using the HES instrument. At present, due to the lacking of a green channel on ABI, no operational ocean color data products have been produced from ABI data at NOAA. If operational ocean color data products were produced, our ability in monitoring water discharging from major US river systems to oceans and tracking turbid water movements over coastal waters could have been greatly improved because of rapid repeating cycles of ABI measurements.

The Japanese AHI instruments on board the Himawari-8/9 geostationary satellite platforms are equipped with a green band. Operational ocean color data products have been generated at the Japan Aerospace Exploration Agency (JAXA). The atmospheric correction algorithm for processing AHI data was described by Murakami [7]. This algorithm was adapted from an earlier version of the atmospheric correction algorithm [8] developed for processing the Japanese Global Imager (GLI) data, while the GLI algorithm shared the same basis with the SeaWiFS algorithm [9]. During the atmospheric correction processes using SeaWiFS-types of algorithms, the pure Rayleigh scattering effects for two bands are first subtracted out. The ratio of the two bands is then calculated. From the ratio value, an aerosol model for a given pixel is estimated. One problem with the ratio calculation for optical bands with additive noises is that it magnified the noise by a factor of 2 [10]. The noise is then propagated down to the derived water leaving reflectances of visible bands. In order to improve signal to noise ratios of measured AHI data, temporal averaging on the hourly basis is required for operational atmospheric corrections to AHI data [7]. Figure 2 shows sample ocean color results downloaded from a public Himawari web site (eroc.jaxa.jp/ptree/index.html). Figure 2a is an RGB image covering eastern part of China–Korea Peninsula, and part of Japan. The image was acquired on May 28, 2017 at UTC 0200. Land/water boundaries are drawing in yellow colored lines. Turbid coastal waters in eastern part of China and around the Korean Peninsula are seen obviously. Figure 2b is the corresponding false color chlorophyll concentration image overlaid on the RGB image. Over portions of turbid coastal water areas, such as those areas pointed by arrows, no retrievals are made. Figure 2c is another case of false color chlorophyll concentration image overlaid on the RGB image. The AHI data was acquired on April 26, 2020 at UTC 0300. In this specific case, no retrievals were made over most water surfaces in eastern coastal areas of China.

Historically, the first geostationary ocean color sensor, the geostationary ocean color imager (GOCI) [11], was launched into space in 2010. This sensor has eight SeaWiFS-like bands. Several atmospheric correction algorithms [12–16] sharing the same base as the SeaWiFS 2-band ratio algorithm but with additional assumptions about absorption and scattering properties of coastal waters in the 0.66–0.865 µm spectral range have been developed. The performance of these algorithms has been objectively evaluated [17]. The majority of multi-spectral and hyperspectral ocean versions of atmospheric correction algorithms have recently been reviewed in a lengthy paper by Frouin et al. [18]. Because GOCI lacks SWIR bands above 1 µm for atmospheric correction purposes, the spatial coverage of atmosphere-corrected GOCI data products over very turbid coastal waters is also frequently lacking.

In view of the absence of ABI ocean versions of atmospheric correction algorithms and the needs for improving spatial coverage over turbid coastal water areas with AHI ocean color retrievals, we have recently developed spectrum-matching versions of atmospheric correction algorithms using SWIR bands above 1 µm, where the turbid coastal water surfaces are often very dark, for retrieving water leaving reflectances from ABI and AHI data. In Section 2 of this article, we describe the spectral characteristics of ABI and AHI instruments and the spectrum-matching technique for water leaving reflectance retrievals. We present sample retrievals from ABI and AHI data in Section 3. We give a short discussion in Section 4. Finally, we provide a brief summary in Section 5.

**Figure 2.** (**A**): A true color RGB image acquired with the Japanese Advanced Himawari Imagers (AHI) instrument over part of Japan, Korean Peninsula, and eastern part of China at UTC 0200 on May 28, 2017; (**B**): Similar to (**A**) but with a false chlorophyll concentration image overlaid; and (**C**): Another case of an AHI false color chlorophyll concentration image overlaid on an RGB image that was acquired on April 26, 2020 at UTC 0300.

#### **2. Data and Methods**

#### *2.1. The ABI and AHI Instrument Characteristics*

The ABI instruments, currently on board the NOAA GOES 16 and 17 geostationary weather satellite platforms, have six bands in the solar spectra range. These bands are centered near 0.47, 0.64, 0.865, 1.378, 1.61, and 2.25 µm. The positions and widths of these bands are illustrated in Figure 3. The bands are plotted above a green algae reflectance spectrum, which was retrieved from an AVIRIS (airborne visible infrared imaging spectrometer) [19] data set acquired over a salt pond in the San Francisco Bay in August of 2009, with a version of our hyperspectral atmospheric correction algorithm (ATREM) [20]. It is obvious to see that ABI does not have a green band to center near the 0.55 µm. To improve the ocean color capability, the Japanese AHI instruments on board the Himawari-8/9 geostationary satellite platforms are equipped with a green band centered at 0.51 µm, which is illustrated in Figure 3. The basic designs of AHI are the same as those of ABI, except filter changes. AHI replaced ABI's 1.378-µm cirrus band with the 0.51-µm green band. The addition of the

green band enhances the water color information that can be obtained from the AHI data in comparison with that from the ABI data.

**Figure 3.** Bandpasses of ABI and AHI instruments. An airborne visible infrared imaging spectrometer (AVIRIS) green algae reflectance spectrum is also shown.

#### *2.2. Atmospheric Corrections for ABI and AHI Data*

At present, most ocean color data products have been generated from polar orbiting satellite data at various data centers. Atmospheric correction algorithms used for data product generations are mostly based on the Gordon and Wang [9] (hereafter referred as GW94) SeaWiFS types of 2-band algorithms that were designed mainly for clear deep ocean waters (case 1 waters). In the GW94 algorithm, the SeaWiFS 0.75-µm and 0.86-µm bands are used as atmospheric correction bands assuming water leaving reflectances from the two bands are close to zero. Over coastal waters, suspended sediments can contribute non-negligible water leaving reflectances to the 0.75- and 0.86-µm bands. Several investigators have tried to remove the NIR water leaving reflectance contributions from the measured top-of-atmosphere NIR apparent reflectances with different techniques [21,22], so that a "black pixel" could be provided to the SeaWiFS 2-band ratio type of atmospheric correction algorithms.

Previously, we developed an ocean version of hyperspectral atmospheric correction algorithm [20] mainly for supporting the Navy COIS (coastal ocean imaging spectrometer) project [23]. For the turbid coastal environment, the water leaving radiances in the 0.75–0.865 µm spectral range are typically not close to zero mainly because of scattering by suspended materials, particularly for the 0.75-µm band. Under these conditions, the channels in this spectral region have limited use for retrieving information on atmospheric aerosols. Because the SeaWiFS types of algorithms derive aerosol information from channels in the 0.75–0.865 µm spectral range, these algorithms cannot be easily adapted for operational retrieval of water leaving radiances over coastal waters. In view of this limitation, we designed a spectrum-matching algorithm [20] that allowed the use of channels in shortwave IR (SWIR) spectral regions between 1 and 2.5 µm, where the turbid waters are much darker, for the estimates of aerosol models and optical depths. Our algorithm was based on Robert Fraser's radiative transfer formulation and algorithm [24,25]. The adoption of the Fraser formulation permits the simultaneous matching

between measured radiances of several bands centered at different wavelengths with those from theoretical simulations and results in more stable estimates of aerosol models and optical depths.

In our algorithm, we express radiances in reflectance units. We adopt the standard definition of apparent reflectance ρ *\* obs* at a satellite level for a given wavelength as [20,26]

$$
\stackrel{\*}{\rho}\stackrel{\*}{\text{obs}} = \pi \text{ L}\_{\text{obs}} / (\mu\_o \text{ E}\_o) \,\tag{1}
$$

where *Lobs* is the radiance of the ocean-atmosphere system measured by a satellite instrument, µ*<sup>o</sup>* the cosine of solar zenith angle, and *E<sup>o</sup>* the downward solar irradiance at the top of the atmosphere when the solar zenith angle is equal to zero. Neglecting the interactions between atmospheric gaseous absorption and molecular and aerosol scattering, ρ *\* obs* can be expressed as [20,26]

$$\stackrel{\ast}{\rho}\stackrel{\ast}{\text{obs}} = T\_{\mathcal{S}} \left[\stackrel{\ast}{\rho}\stackrel{\ast}{\text{atm}+\text{sfc}} + \rho\_{\text{w}} \stackrel{\ast}{\text{t}}\stackrel{\ast}{\text{t}}\!\_{\text{ul}}\!/(\text{1}-\text{s}\;\!\!\rho\_{\text{w}})\right] \tag{2}$$

where *T<sup>g</sup>* is the total atmospheric gaseous transmittance on the Sun-surface-sensor path, ρ \* *atm*+*sfc* the reflectance resulted from scattering by the atmosphere and specular reflection by ocean surface facets, *t<sup>d</sup>* the downward transmittance (direct + diffuse), and *t<sup>u</sup>* the upward transmittance, *s* the spherical albedo that takes account of reflectance of the atmosphere for isotropic radiance incident at its base, and ρ*<sup>w</sup>* the water leaving reflectance.

Solving Equation (2) for ρ*<sup>w</sup>* yields

$$\boldsymbol{\rho}\_{w} = (\boldsymbol{\rho}^{\star}\_{\text{obs}}/T\_{\text{g}} - \boldsymbol{\rho}^{\star}\_{\text{ atm} + \text{sfc}}) \text{[} \mathbf{t}\_{d} \,\, \mathbf{t}\_{u} + \text{s} \,\, (\boldsymbol{\rho}^{\star}\_{\text{obs}}/T\_{\text{g}} - \boldsymbol{\rho}^{\star}\_{\text{ atm} + \text{sfc}}) \text{]} \tag{3}$$

Given a satellite measured radiance, the water leaving reflectance can be derived according to Equations (1) and (3) provided that the other quantities in the righthand side of Equation (3) can be modeled theoretically. Because of the availability of the vector radiative transfer code, the proper atmospheric layering structure in this code (up to 80 layers with proper mixing of aerosol particles and atmospheric molecules in each layer), and the treatment of wind-roughened water surfaces, we have used a modified version of Ahmad and Fraser code [27] to generate lookup tables used in our retrieving algorithm. Specifically, we use the code to generate the quantities ρ *\* atm*+*sfc*, *td, tu*, and *s* in Equation (2). Lookup tables for 14 wavelengths between 0.39 and 2.5 µm in atmospheric "window" regions, sets of aerosol models, optical depths, solar and view angles, and surface wind speeds have been generated. Aerosol models, similar to those used in the SeaWiFS algorithm, are used during our table generation. A sophisticated line-by-line based atmospheric transmittance code is used in our algorithm to calculate contiguous atmospheric gaseous transmittance spectra (*Tg*).

We match the satellite measured "apparent reflectances" of selected NIR and SWIR bands with those pre-computed apparent reflectances for different aerosol models and optical depths. The simulated lookup table apparent reflectances (after interpolation) have step sizes of 0.001 or smaller. The fine step sizes are smaller than the noise levels of satellite measured data. During the retrieving process for a given pixel, we search through all possible combinations of aerosol models and optical depths in lookup tables, and find the best match of the aerosol model and optical depth that minimizes the sum of the squared differences between the measured apparent reflectances (two or more NIR and SWIR bands) and those from theoretical simulations. Because of large computer memory sizes for storing lookup tables and because of the array processing capability of our codes written in Fortran 90, our table searching and minimization process is very fast. We did not need to implement an additional numerical optimization method, such as the commonly used gradient descent method, to speed up the table searching and minimization process. For example, in one test case, we processed one hyperspectral data set having 614 pixels in the cross-track direction, 1024 lines in the along-track direction, and 224 spectral bands, using a Macbook computer having a 2.9 GHz Dual-Core Intel Core i5 processor. It took about 500 megabytes of memory and 140 seconds of computer time to finish the data processing. Our hyperspectral atmospheric correction algorithm has been tested with images

acquired with a number of instruments, including the airborne visible infrared imaging spectrometer (AVIRIS) [19] from an ER-2 aircraft at an altitude of 20 km.

We also adapted the hyperspectral ocean atmospheric correction algorithm for processing multi-channel imagery, such as those acquired with the Terra and Aqua MODIS instruments [28]. In the MODIS algorithm, the lookup tables corresponding to MODIS wavelengths are obtained through linear interpolation of the tables generated for the hyperspectral atmospheric correction algorithm. The lookup table quantities, ρ *\* atm*+*sfc*, *td, tu*, and *s*, are functions of wavelength, solar zenith angle, view zenith angle, relative azimuth angle, aerosol model, optical depth, relative humidity, and surface wind speed, respectively. For MODIS data, each pixel's solar zenith angle, view zenith angle, and relative azimuth angle are variable because of the large swath width (~2800 km for one scan line). In order to speed up the retrieving process, the storage order for different MODIS lookup table quantities (such as wavelengths, particle size distributions, solar and view angles, etc.) is optimized. Because MODIS has water vapor bands centered near 0.94 µm, we derive water vapor values from MODIS data on the pixel by pixel basis, and then use the derived water vapor values to calculate water vapor transmittances for other bands used for retrievals. Climatological atmospheric ozone, carbon dioxide, and methane amounts are used for modeling relevant gaseous transmittances. Similarly, we developed a VIIRS ocean version of algorithm. The VIIRS channels centered at 1.24, 1.61, and 2.25 µm with proper modeling of atmospheric CO<sup>2</sup> and CH<sup>4</sup> absorption effects and a spectrum-matching technique were used for atmospheric corrections.

Based on our experience in developing MODIS and VIIRS versions of atmospheric correction algorithms, we have recently developed spectrum-matching versions of ABI and AHI multi-channel atmospheric correction algorithms. The main tasks in the ABI and AHI algorithm development were to generate proper transmittance tables of atmospheric gases, including those of ozone, water vapor, carbon dioxide, and methane, for ABI and AHI bands in the 0.4–2.3 µm, as illustrated in Figure 2. The spectrum-matching technique is also used for atmospheric corrections and for the subsequent retrieval of water leaving reflectances from ABI and AHI data. Figure 4 illustrates the spectrum-matching technique.

Figure 4a is a false color RGB image (Red: 0.637; Green: 0.864; Blue: 0.471 µm) acquired at UTC 1100 on December 17, 2018 with ABI on board the GOES-16 satellite platform. The image covered algae blooming areas off the coast of Argentina in the Atlantic Ocean. Figure 4b illustrates the spectrum-matching technique for water leaving reflectance retrievals. The red line (marked with the plus symbol) in Figure 4b shows the apparent reflectances for 5 ABI bands for a water pixel located at the center of the red square in the right side of the Figure 4a image. The 0.86-, 1.61-, and 2.24-µm bands were used for the estimation of aerosol model and optical depth by minimizing the sum of the squared differences between the measured data and calculated data during the spectrum-matching process. The estimated aerosol model and optical depth were extracted back to the ABI blue and red bands. The green line (marked with "\*") shows the estimated atmospheric path radiances in reflectance unit, ρ *\* atm*+*sfc* (as defined in Equation (2)), for all the 5 ABI bands. The differences between the red line and the green line are proportional to water leaving reflectances (see Equation (3)). The blue line (marked with the triangle symbol) shows the water leaving reflectances of ABI bands for the algae pixel. The water leaving reflectance value for the blue band is close to 0.1, indicating that the pixel is highly reflecting in the blue spectral region.

**Figure 4.** (**A**)—a false color ABI image acquired at UTC 1100 on December 17, 2018 off the coast of Argentina covering algae blooming regions; (**B**)—spectral plot for a bright pixel for top of atmosphere (TOA) apparent reflectances (red line), atmospheric reflectances (green line), and water leaving reflectances (blue line); (**C**)—the retrieved false color water leaving reflectance image for the UTC 1100 scene; and (**D**)—the derived false color water leaving reflectance image for the UTC 1800 scene.

#### **3. Results**

We have selected three examples to demonstrate the capabilities in deriving water leaving reflectances from ABI and AHI data. The selection of the ABI and AHI cases was based on ABI and AHI time sequences of movie images available from a public web site associated with the Cooperative Research Program at Colorado State University. The first example is for the ABI scene covering a major blooming event off the coast of Argentina (see Figure 4a). The second example is for an ABI scene over the Gulf of Mexico where large amounts of sediments were discharged from rivers into the Gulf of Mexico. The third case is for a Japanese AHI scene for the coastal areas of Japan, Korea Peninsula, and eastern part of China, where very turbid coastal waters were present.

#### *3.1. ABI Scene Over Coastal Area of Argentina, December 17, 2018*

The ABI instrument on GOES-16 captured a major chlorophyll blooming event occurred in the middle of December of 2018. As described in Section 2, Figure 4a is a false color RGB image acquired at UTC 1100 on December 17 over algae blooming areas off the coast of Argentina in the Atlantic Ocean. Figure 4b shows the results of matching technique for water leaving reflectance retrievals from ABI data. Because ABI didn't have a true green band, it is not possible to observe the green reflectance peak centered near 0.55 µm in the Figure 4b water leaving reflectance curve. Figure 4c is the false color RGB image processed from the retrieved water leaving reflectances of the Figure 4a UTC 1100 scene. In principle, water leaving reflectance images, similar to the one illustrated in Figure 4c, at a time interval of approximately 10 minutes can be obtained from ABI data because of the rapid repeating cycle of ABI measurements. Figure 4d shows another false color ABI water leaving reflectance image for UTC 1800. Although the NASA MODIS instrument onboard the Aqua Spacecraft acquired image at UTC 1755 over the same area and covered the same chlorophyll patches, the NASA operational MODIS algorithm didn't make ocean color retrievals over areas containing the chlorophyll features in the MODIS scene. As a result, it is not possible to make quantitative comparisons between our retrieval results with those from the NASA operational algorithm. Nevertheless, the images and line plots in Figure 4 demonstrate that water leaving reflectance images can be retrieved from ABI data with our ABI version of multi-channel atmospheric correction algorithm.

#### *3.2. ABI Scene Over Coastal Area of Gulf of Mexico, March 20, 2019*

Figure 5 shows another example of water leaving reflectance retrieval from ABI data. The data sets were acquired on March 20, 2019 off the coastal area of Louisiana and Texas. At the time, large amounts of sediments were discharged from rivers into the Gulf of Mexico. Figure 5a,b is a false color RGB image (R: 0.637; G: 0.864; B: 0.471 µm) acquired at UTC 1730 and 1930, respectively. The brownish-colored areas contained large amount of sediments. Figure 5c,d is the corresponding false color water leaving reflectance images. The sediment-dominated areas were picked up very nicely in the Figure 5c,d. The ABI 1.61-, and 2.24-µm bands were used during atmospheric corrections. Figure 5e illustrates the atmospheric correction for one sediment-dominated pixel located near the center of the rectangle outlined in the Figure 5b,d images. The six points marked with "\*" are the top of atmosphere (TOA) apparent reflectances for the six ABI bands. The five points marked with small green-colored triangles are the retrieved water leaving reflectances for the pixel. The atmosphere-corrected 0.637-µm red band reflectance is about 0.12, which indicates that the sediment-dominated pixel is highly reflecting in the red band. Although the Aqua MODIS instrument had a good overpass over the same area at UTC 1925, the sediment-dominated areas were masked out by the NASA MODIS atmospheric correction algorithm and no water leaving reflectance retrievals were made. Figure 5 demonstrates that ABI data can be used for water leaving reflectance retrievals over sediment-dominated areas. The data products from ABI can be complimentary to the existing NASA polar-orbiting dominated ocean color data products, particularly over turbid coastal waters.

For comparison and validation purposes, we have also retrieved water leaving reflectances from a MODIS data set acquired at UTC 1925 on the same day over similar areas. Figure 6a is a true color MODIS RGB image. Figure 6b is a false color image (R: 0.645 µm; G: 0.845 µm; B: 0.466 µm) of the same scene. Figure 6b covers nearly the same area as that of Figure 5b. Major land and water surface features in both images are approximately the same. However, because the Aqua spacecraft has a declination angle of 98.5 degree and because no spatial registration is made between the MODIS image and the ABI image, the small turbid water areas in the upper right portion of Figure 6b are not covered in Figure 5b. Figure 6c is the false color water leaving reflectance image (R: 0.645 µm; G: 0.845 µm; B: 0.466 µm). Based on visual inspection of images in Figures 6c and 5c, we selected nearly same spatial areas and extracted spectral reflectances. Figure 6d shows water leaving reflectances derived from ABI (dashed green curve) and from MODIS data (solid line) using narrow bands near 0.86, 1.61, and 2.25 µm. By comparing the two curves in Figure 6d, it is seen that the ABI band reflectance values agreed quite well with those of the corresponding MODIS band values. The good agreement partially validated our ABI version of the atmospheric correction algorithm. It should be pointed out that, because MODIS has more narrow bands, the MODIS curve provides more spectral information than the ABI curve. For example, MODIS has a band centered near 0.52 µm and this band has much larger

reflectance value than the two nearby bands, while ABI does not have a band centered near 0.52 µm and has less capability in charactering the spectral reflectance properties of water surfaces.

**Figure 5.** (**A**,**B**): false color ABI images acquired at UTC 1730 and 1930, respectively on March 20, 2019 near the Gulf of Mexico; (**C**,**D**): the corresponding derived water surface reflectance images; (**E**): spectral plot for a pixel before and after atmospheric corrections.

**Figure 6.** (**A**,**B**)—true color and false color moderate resolution imaging spectroradiometer (MODIS) images acquired near the Gulf of Mexico at UTC 1925 on March 20, 2019; (**C**)—the corresponding derived water surface reflectance images; and (**D**)—water leaving reflectances derived from ABI (dashed green curve) and from MODIS data (solid line).

#### *3.3. AHI Scene Over Eastern Asia, May 28, 2017*

Figure 7 shows an example of water leaving reflectance retrievals from the Japanese AHI data measured on May 28, 2017 over part of Japan, Korea, and eastern China. Figure 7a is the RGB image (R: 0.64-; G: 0.51-; B: 0.47-µm) acquired at UTC 0200. It can be difficult to distinguish between bright coastal waters in these visible band images, particularly in the left portion of the scene. Figure 7b shows the AHI image for the 1.61-µm band. A large inland lake (Lake Taihu located near the lower left corner), rivers, and coastal waters are quite dark in this image because of strong liquid water absorption at 1.61 µm. Through several trials and tests, we decided to select a 1.61-µm band threshold. Pixels with the 1.61-µm band TOA apparent reflectance values greater than 0.03 were classified as land or cloudy pixels. These pixels were masked out in our water leaving reflectance data product. Figure 7c shows the retrieved water leaving reflectance image for the scene. With the SWIR spectrum-matching algorithm, we are able to do retrievals over the deep ocean waters, bright coastal waters, and inland lakes and rivers. By comparing Figure 7c with Figure 2b, it seen that we have greatly expanded the retrievals over very turbid water areas in eastern part of China and around Korea Peninsula. In particular, we are able to do retrievals over Lake Taihu, which is seen in the lower left part of Figure 7c. Figure 7d illustrates the atmospheric correction for one turbid water pixel in west coastal area of the Korea Peninsula using two narrow channels centered near 1.61 and 2.25 µm for aerosol retrievals. The six points marked with the "+" symbol are the TOA apparent reflectances for the six AHI bands. The six points marked with small diamonds are the retrieved water leaving reflectances for the pixel. The near zero surface reflectance values for the AHI bands centered near 0.865, 1.61, and 2.25 µm demonstrate that proper minimization is made during the spectrum-matching process.

**Figure 7.** (**A**)—a true color AHI image acquired at UTC 0200 on May 28, 2017 over part of Japan, Korea, and eastern China; (**B**)—the 1.61-µm AHI B/W image; (**C**)—the retrieved true color water leaving reflectance image; and (**D**)—spectral plot for a bright pixel for TOA apparent reflectances marked with the "+"symbols and for the corresponding surface reflectances marked with small diamonds.

#### **4. Discussions**

It should be pointed out that water leaving reflectance retrievals, such as those showing in Figures 4, 5 and 7, are made on the pixel by pixel basis with our spectrum-matching versions of ABI and AHI atmospheric correction algorithms. No spatial and temporal averaging of ABI and AHI data is required during the retrievals. Similar retrievals with the SeaWiFS-type of algorithm adopted for AHI data processing by Murakami [7] required temporal averaging on the hourly basis. In the original operational SeaWiFS algorithm, spatial averaging of the 0.75- and 0.86-µm bands was made for improving signal to noise ratios of the two atmospheric correction bands. The calculation of two band ratios magnifies the noise by a factor of two [10]. The spectrum-matching algorithms described here do not need the calculation of two-band ratios. In fact, if two bands are used in the spectrum-matching process (e.g., Figures 5 and 7), the signal to noise ratio is increased by a factor of square root of two. If three bands are used in the matching process (e.g., Figure 4), the signal to noise ratio is increased by square root of three [10]. Therefore, the use of spectrum-matching algorithms for atmospheric corrections and for derivation of water leaving reflectances in the visible requires less signal to noise

ratios of the original satellite-measured data than the use of SeaWiFS-type of two-band ratio algorithms. This fine point is very important to atmospheric corrections for the present generation of geostationary weather satellite data, because the relevant instruments were not geared for ocean color measurements with exceedingly high signal to noise ratios. In addition, performing spatial and temporal averaging to improve signal to noise ratios can sometimes be difficult. For example, during a time interval of one hour for the day May 28, 2017 (see Figure 6), the aerosol properties over a given area in the Figure 6 scene changed quite a bit because of the movement of a dust storm through the area. Temporal averaging is not quite suited for improving signal to noise ratios for atmospheric correction purposes for this particular day.

The details on our radiative transfer simulations [27], generation of lookup tables, and spectrum-matching processes were previously described in our hyperspectral ocean version of atmospheric correction algorithm [20] and multi-channel MODIS algorithm [28]. The water leaving reflectance retrieving results from MODIS data were verified with measurements with a portable field spectrometer over surface stations in Bahamas Banks and Florida Bay area [28]. The results were also compared with those from the NASA operational MODIS algorithm. Because the lookup tables used in our ABI and AHI algorithms were obtained through interpolation of the same hyperspectral lookup tables simulated with the Ahmed and Fraser radiative transfer code [27] and because our hyperspectral and MODIS multi-spectral retrieving results were previously validated, we feel that it is not absolutely necessary to have additional validation of our several sample retrieving cases from ABI and AHI data. In addition, we are now having practical difficulties in obtaining additional ABI and AHI data sets for more case studies, because the public ftp access is no longer permitted at our computing facility.

In the commonly used SeaWiFS types of ocean version of atmospheric correction algorithms, the atmospheric path radiances in reflectance unit after subtraction of Rayleigh scattering effects are stored in lookup tables, such as the tables used in the NASA operational SeaWiFS, MODIS, and VIIRS algorithms. Under very clear atmospheric conditions, the Rayleigh-subtracted atmospheric path reflectances for bands centered near 1.24, 1.61, and 2.25 µm can be very small and close to zero. The extension of the SeaWiFS types of 2-band ratio algorithms to SWIR band pairs, such as the 1.24- and 1.61-µm band pair and the 1.24- and 2.25-µm band pair, can encounter the problems of zero divided by zero or zero divided by negative numbers due to the presence of noises in measured SWIR band data. As a result, the use of 2-SWIR band ratio technique for the estimates of aerosol models and optical depths during atmospheric correction processes can be problematic. The multi-channel ocean color research community does not seem to well understand the intrinsic problems [29].

#### **5. Summary**

We have developed atmospheric correction algorithms for retrieving water leaving reflectances from ABI and AHI data using spectrum-matching techniques. Our spectrum-matching algorithms reduced significantly the signal to noise ratio requirements for atmospheric correction channels in comparison with those required with the widely used SeaWiFS types of two-channel ratio algorithms by the current multi-channel ocean color research community. We have demonstrated using three case studies that water leaving reflectance retrievals can be made from satellite-measured ABI and AHI data with our spectrum-matching versions of atmospheric correction algorithms. Our algorithms, if implemented onto operational computing facilities, can be complimentary to current operational ocean versions of atmospheric correction algorithms in the area of improving spatial coverage over turbid coastal waters.

**Author Contributions:** B.-C.G. made initial observations of chlorophyll blooming features and suspended sediment features from ABI and AHI data sets, and developed spectrum-matching versions of algorithms for retrieving water leaving reflectances from ABI and AHI data. R.-R.L. carried out retrievals and detailed data analysis. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is partially supported by the US Office of Naval Research and by a research grant managed by Jared K. Entin at the Science Mission Directorate of National Aeronautics and Space Administration.

**Acknowledgments:** The authors are grateful to Weile Wang affiliated with California State University in Monterey and NASA Ames Research Center for providing the ABI and AHI data in HDF4 format used in this study. The authors are also grateful to Pengwang Zhai of University of Maryland at Baltimore County for useful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Synergy between Satellite Altimetry and Optical Water Quality Data towards Improved Estimation of Lakes Ecological Status**

**Ave Ansper-Toomsalu 1,\*, Krista Alikas <sup>1</sup> , Karina Nielsen <sup>2</sup> , Lea Tuvikene <sup>3</sup> and Kersti Kangro 1,3**


**Abstract:** European countries are obligated to monitor and estimate ecological status of lakes under European Union Water Framework Directive (2000/60/EC) for sustainable lakes' ecosystems in the future. In large and shallow lakes, physical, chemical, and biological water quality parameters are influenced by the high natural variability of water level, exceeding anthropogenic variability, and causing large uncertainty to the assessment of ecological status. Correction of metric values used for the assessment of ecological status for the effect of natural water level fluctuation reduces the signal-to-noise ratio in data and decreases the uncertainty of the status estimate. Here we have explored the potential to create synergy between optical and altimetry data for more accurate estimation of ecological status class of lakes. We have combined data from Sentinel-3 Synthetic Aperture Radar Altimeter and Cryosat-2 SAR Interferometric Radar Altimeter to derive water level estimations in order to apply corrections for chlorophyll a, phytoplankton biomass, and Secchi disc depth estimations from Sentinel-3 Ocean and Land Color Instrument data. Long-term in situ data was used to develop the methodology for the correction of water quality data for the effects of water level applicable on the satellite data. The study shows suitability and potential to combine optical and altimetry data to support in situ measurements and thereby support lake monitoring and management. Combination of two different types of satellite data from the continuous Copernicus program will advance the monitoring of lakes and improves the estimation of ecological status under European Union Water Framework Directive.

**Keywords:** ecological status class of lakes; European Union Water Framework Directive (2000/60/EC); water quality parameters; water level; Sentinel-3; Cryosat-2; shallow lakes; synergy; altimetry data; optical data

#### **1. Introduction**

Consistent monitoring and estimation of the ecological status of lakes is essential for having an overview of current conditions and prediction of water quality (WQ) [1,2]. As lakes regulate the climate and carbon cycle, provide habitat for biota, drinking water, and other ecosystem services, it is important to monitor their ecological status and indications of their decline. To implement water management for sustainable environment, it is essential to evaluate the influence of human activity to lake ecosystems and also take into account variability of natural factors [3,4]. The European Union (EU) has established Water Framework Directive (2000/60/EC) [5] (WFD) for all European countries to achieve at least "good" ecological status class in lakes through consistent monitoring and water policy by the latest 2027 [6,7]. In total, five ecological quality classes from "high" to "bad" have been determined for every defined natural water body type, according to different biological (e.g., phytoplankton), physico-chemical (e.g., transparency) and hydro-morphological (e.g.,

**Citation:** Ansper-Toomsalu, A.; Alikas, K.; Nielsen, K.; Tuvikene, L.; Kangro, K. Synergy between Satellite Altimetry and Optical Water Quality Data towards Improved Estimation of Lakes Ecological Status. *Remote Sens.* **2021**, *13*, 770. https://doi.org/ 10.3390/rs13040770

Academic Editor: Giacomo De Carolis

Received: 30 December 2020 Accepted: 17 February 2021 Published: 19 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

depth variations) quality elements [5]. The main emphasis is given to the biological quality elements while the physico-chemical and hydro-morphological quality elements only affect the decision indirectly through their effect on the state of biological elements [5].

Traditional in situ measurements (e.g., collection of water samples, laboratory analyses, gauging stations, etc.) tend to be money and time consuming, and tend to take a lot of resources [8]. It is difficult to collect data from remote areas, where infrastructure is missing and access to the water body is complicated [9]. In situ gauging stations provide more frequent data than traditional measurements, but transmission of the data to the further workflow is rather slow and time consuming [10–12]. These kinds of systems are usually expensive to buy and difficult to maintain, therefore the number of gauging stations has decreased in recent years [11–14]. As lakes are mostly categorized as heterogeneous bodies, the biggest disadvantage of traditional in situ measurements is providing only point measurements, which do not represent spatial and temporal variations of the entire water body [1,11,15]. The combination of Earth Observation (EO) and in situ measurements complement each other by fulfilling their shortages and also helps to optimize and improve the selection of sampling locations and surveying times [16–18]. Earth Observation gives opportunity to have large scale data for spatial analysis and synchronized view of group of lakes [8,19]. In addition, digital data forms accelerates data reading and processing time and helps to produce historical and long time series [15]. The European Space Agency's (ESA) Copernicus program [20] provides full, free, open access, high quality near-realtime EO data on a global level for all users. Since the launch of the first Sentinel mission satellite in 2014, the purpose of the program is to launch nearly 20 satellites by the year of 2030 to provide long-term time series for various applications, e.g., environmental protection, climate change monitoring, and civil security [21]. The continuity of data allows to develop and advance methods to better integrate EO data into environmental monitoring, management, and policy.

As the satellite era has lasted almost 40 years by now, many algorithms and methods have been developed for different domains in lake research [10,22–31]. As the WFD requirements are obligatory to every EU country, it has been indicated that the in situ methods are not feasible to cover the monitoring expectations, therefore many EO-based applications have been developed for implementation of requirements of WFD [32]. According to the WFD, phytoplankton is the key parameter of estimating ecological status of water [33–36], which also indicates eutrophication status of the lake [37]. Estimation of the concentration of the photosynthetic pigment chlorophyll a (chl a) has been one of the main goals in all ocean color missions starting from the first National Aeronautics and Space Administration (NASA) Coastal Zone Color Scanner on board Nimbus-7 in 1978. Additionally, water transparency is one of the key indicators for lake's ecological status class according to the WFD [38–40], where also EO-based assessment methods have shown the suitability to obtain additional information [7,25,32,41]. Covering various ocean color missions, the approaches and algorithms have been developed to provide trophic status based estimations [31,42–48] to obtain chl a and transparency estimates. However, estimation is based only on certain parameters (e.g., chl a and transparency), which do not take into account any external factors [34,35,38,49]. Ranges of different ecological status classes of lakes may be very narrow, which is the main reason why accurate WQ parameter estimation is essential. Water level (WL) estimation, using satellite altimetry data over inland water bodies, has become more advanced during recent years [10,16,17,50–52], which allows the monitoring of its influence on the physical environment, which in turn changes the biota and whole lake ecosystem [53], especially in the case of shallow lakes.

In general, the ecological status of lakes is affected by anthropogenic factors, such as agricultural pollution by nutrients and industrial waste and various chemical elements [54], nevertheless natural factors, such as diurnal and seasonal changes, and prolonged changes in atmospheric circulations should be taken into account as well [55–57]. Water level is one of the natural factors, which affects WQ directly, worsening WQ parameters in low WL periods [58,59]. Especially in large, shallow lakes, the influence of lake depth variations on

biota is remarkable, due to large seasonal and inter-annual WL fluctuations. The periods of extremely low WL are the most critical for the ecological status in such lakes, strengthening sediment resuspension, increasing concentrations of ions, phytoplankton biovolume and chl a concentration [60]. In contrast, the phytoplankton biomass (FBM) is significantly lower due to light limitation in years of high WL [61]. Such high variation in biological quality indicators, caused by natural factors, may over-mask human induced effects on a water body [62]. Based on the long-term in situ data, Tuvikene et. al. [62] showed the significant influence of WL while estimating the ecological status class according to the chl a, FBM and transparency in case of large shallow lakes. It was shown that, while phytoplankton species composition and traditional WQ parameters indicated different ecological status class, correction of metrics for WL enabled more consistent assessment of the ecological status of the lake.

As assessment of ecological status of lakes is currently parameter-based, thus accounting for the external natural factors is needed for more precise estimation. In the increasing constellation of ocean color satellites, carrying passive optical sensor usable for mapping various optical WQ parameters, and altimeters for estimating WL, the synergy of data could potentially improve the ecological status assessment over large shallow lakes. Therefore, the aim of this paper is to (1) create synergy between the optical and altimetry data in order to, (2) derive chl a, total FBM and transparency based ecological status classes taking into account the dependence of biological status indicators on the changing WL over large shallow lakes.

#### **2. Materials and Methods**

#### *2.1. Study Sites Description*

Data from two large and shallow lakes are used for this study. Võrtsjärv is situated in the central part of Estonia (Figure 1) with surface area 270 km<sup>2</sup> (length 35 km, width ~14 km). With the average depth of 2.8 m (the maximum depth is 6 m), the lake has strong WL fluctuation as the annual mean amplitude is 1.4 m and absolute range is 3.1 m, causing up to 3.2-fold difference in the lake volume [61]. Water level is generally characterized by a high period in spring, due to snow melting until end of April or beginning of May, a low period in winter and summer, and increasing WL from middle of the October due to autumn rainfalls. Võrtsjärv is a eutrophic lake with low water transparency (<1 m) [63]. According to the Estonian regulation of Water Act [64], Võrtsjärv belongs to independent ecological status class Type 6 (Table 1) (non-stratified, light water color, low amount of chlorides, moderate hardness, surface area 100–300 km<sup>2</sup> ). Annually, national monitoring is performed at one station once in a month and from additional 2 or 9 stations in August [65].

Peipsi is a transboundary lake situated in the border of Estonia and Russia (Figure 1) with surface area 3555 km<sup>2</sup> (length 152 km, width 48 km). The average depth is 7.1 m (the maximum depth is 15.3 m) and the annual mean amplitude of WL change is 1.5 m. Water level of Peipsi is characterized the highest in spring during the ice melting period and is the lowest in October [66]. It is a eutrophic lake, which consists of three parts: Peipsi *sensu stricto* (*s.s.*, 2611 km<sup>2</sup> average depth 8.3 m), Lämmijärv (236 km<sup>2</sup> average depth 2.5 m) and Pihkva (708 km<sup>2</sup> , average depth 3.8 m). According to the Estonian regulation of Water Act [8], Peipsi belongs also to the independent ecological status class Type 7 (Table 1) (non-stratified, light water color, low amount of chloride, moderate hardness, surface area >1000 km<sup>2</sup> ) and is monitored 7 times per year (in March, when the lake is ice-covered, and from May to October) from 7 stations [65].

The both lakes are eutrophic and shallow, with similar WL fluctuation in terms of the annual mean (Võrtsjärv 1.4 m and Peipsi 1.5 m), but with different average depth (Võrtsjärv 2.8 m and Peipsi 7.1 m) leading to differences in the effect of WL. In addition, both lakes were considered unique, having no undisturbed reference sites for comparison, therefore, pressure-based approach, morphoedaphic index, and historical data were used to derive site-specific reference conditions and classification criteria for them [67].

**Figure 1.** Locations of in situ stations in Võrtsjärv (Vx) and Peipsi (Px). Peipsi consists of three parts: Peipsi *s.s.* (light purple), Lämmijärv (pink), Pihkva (light green); Võrtsjärv is colored as orange. Buffers of 1 and 2 km from the shore are created for Võrtsjärv and Peipsi respectively. Sentinel-3- altimetry tracks are shown in light blue, whereas Cryosat-2 tracks are shown in brown. Locations of hydrometric stations are numbered as 1 Rannu-Jõesuu (blue), 2 Mustvee (green), 3 Praaga (yellow), 4 Mehikoorma (red).

**Table 1.** Thresholds for chlorophyll a (chl a) (mg/m<sup>3</sup> ), phytoplankton biomass (FBM) (g/m<sup>3</sup> ) and Secchi disk depth (SDD) (m) estimating ecological status class in Võrtsjärv and Peipsi according to Estonian regulation of Water Act [64]. Each color represents ecological status class defined by Water Framework Directive (WFD) [5] and the range of measured values for each parameter during the study period from both lakes are shown in the last column, whereas the number in the brackets indicates average of the measured values.



**Table 1.** *Cont.*

#### *2.2. In Situ Data*

#### 2.2.1. Water Level

Water level was measured by the Vaisala automatic weather station named as MAWS110, which sensor was located near the bottom of the lake and it measured the pressure of the water column, which is converted to the WL in centimeters [68]. Measurements were performed every hour, but WL was averaged over a day for the study. To derive WL above sea level (in Estonian Height system (EH2000), which corresponds to European Vertical Reference System), height of the sensor from sea level (national geoid named as EST-GEOID2017 (Estonian national geoid model)) [69] was added. Hydrometric stations were located near the shore or in the mouth of the river (Figure 1). Water level data was available in Rannu–Jõesuu (hereafter Rannu) hydrometric station from 1927, in Mustvee and Praaga from 1921 and in Mehikoorma from 1949 until 2019.

#### 2.2.2. Chl a, FBM, and SDD

Chl a and FBM samples of integral water at discrete depths (0.5 m, 1 m, etc., until 0.5 m from the bottom) were collected from different stations and analyzed in a laboratory from the 1960s. Chl a samples were filtered through Whatman glass-microfiber filters (GF/F) ø 25 mm filter and pigment extraction was done with 5 mL 96% ethanol [70]. Chl a was measured spectrophotometrically with Hitachi U-3010 (430–750 nm) and the concentrations of chl a were derived according to Jeffrey and Humphrey [71]. For determination of FBM, up to 400 counting units were counted and measured with an inverted microscope based on Utermöhl technique [72]. SDD values were measured with white Secchi disc from the 1950s.

#### *2.3. Satellite Data*

Sentinel-3A/3B (S3A/S3B) and CryoSat-2 (C2) data from 2016 to 2019 covering the period from May to October were used. Buffers of 1 km (Võrtsjärv) and 2 km (Peipsi) were created to reduce the possible errors in the vicinity of land in case of optical and altimetry data (Figure 1). Depending on the cloud cover, 300 m resolution data from S3 Ocean and Land Colour Instrument (OLCI) is available daily in the study location. Surface elevations from the Synthetic Aperture Radar Altimeter (SRAL) sensor onboard S3A and S3B are available every 27 days for a given track. Lake Peipsi is measured by both S3A and S3B, hence, data are available from 2016 and onwards. Lake Võrtsjärv is only measured by S3B, hence data are available from the launch in 2018. C2, carrying Synthetic Aperture Radar/Interferometric Radar Altimeter (SIRAL) sensor, is operating in a geodetic orbit with a repeat period of 369 days. S3 and C2 tracks during the study period are shown on the Figure 1 in blue and gray, respectively.

#### 2.3.1. Optical Data

S3 OLCI Level-1 (L1B) Full Resolution "Non Time Critical" images were used from Estonian National Satellite images database named as ESTHub [73]. Data were used from

the beginning of S3 operation, respectively, S3A since 2016, and additional data from S3B since 2018. L1B data was filtered by quality flags ("bright", "sun\_glint\_risk", "duplicated", "invalid") and phytoplankton related parameters were derived via lake specific empirical relationship (Equations (2) and (3)) [74] applied on Maximum Chlorophyll Index (MCI) (Equation (1)) [75]:

$$M\text{CI} = L\_{709} - L\_{681} - \frac{709 - 681}{753 - 681} \times (L\_{753} - L\_{681}) \tag{1}$$

where *L<sup>x</sup>* corresponds to top-of-atmosphere radiance at certain wavelength.

$$\text{Chl } a = 10.9 \times \text{MCI} + 15.3 \tag{2}$$

$$FBM = 5.8 \times MCI + 5.4\tag{3}$$

For water transparency estimates, OLCI L1B data was processed with Polymer v4.10 for retrieving water-leaving reflectance (*ρ*w). Processed data were filtered by quality flag Bitmask ("0", "1024", "1023"). For the SDD estimates, first the weighting function (Equation (4)) [39] was applied on *ρ*<sup>w</sup> which defined the appropriate empirical algorithm to derive diffuse attenuation coefficient of downwelling irradiance (*K<sup>d</sup>* (490)) (Equation (5)).

$$\mathcal{W} = 5.098 - 2.2099 \left( \frac{\rho(560)}{\rho(709)} \right) \tag{4}$$

$$K\_{\rm f}(490) = (1 - W) \times \exp(-0.9 \ast (\log\_{10} \left(\frac{\rho\_{490}}{\rho\_{799}}\right) + 0.4) + 0.02) + W \times \exp(-1.2 \ast (\log\_{10} \left(\frac{\rho\_{560}}{\rho\_{799}}\right) + 1.4) + 0.02) \tag{5}$$

Calculated *W* ≤ 0 corresponds to more transparent waters and algorithm (*ρ*490/*ρ*709) was used by assigning *W* = 0. Calculated *W* ≥ 1 corresponds to more turbid waters and the algorithm *ρ*560/*ρ*<sup>709</sup> was used by assigning *W* = 1. When the *W* is between 0 and 1, then the algorithm merges the retrievals from both (*ρ*490/*ρ*<sup>709</sup> and *ρ*490/*ρ*709) algorithms [39]. Secondly, the SDD was derived based on the empirical algorithm from *K<sup>d</sup>* (490):

$$SDD\_{\text{Vortsjarv}} = -0.038 \times K\_d (490) + 0.9448 \tag{6}$$

SDD for Peipsi *s.s.* (Equation (7)) was calculated using band ratio *K<sup>d</sup>* (490) algorithm (Equation (5)) [38]:

$$SDD\_{\text{Peips}\,s.s.} = 3.42 \times K\_d (490)^{-0.82} \tag{7}$$

SDD for Lämmijärv and Pihkva was calculated using empirical algorithm (Equation (8)) [38]:

$$SDD\_{\text{Lammijarv}, \text{Pihva}} = \left(\frac{\rho(490)}{\rho(709)}\right)^{0.7} \times 2.14 \tag{8}$$

#### 2.3.2. Altimetry Data

To construct water surface elevations from S3A and S3B the ESA "Non Time Critical" Level-2 product "Enhanced measurements" was used. This product consists of 20Hz alongtrack measurements of, e.g., tracker range, satellite altitude, atmospheric and geophysical corrections, and waveforms. To identify the nadir water surface, the waveforms were retracked with a narrow primary peak threshold retracker [76]. For C2 we applied the ESA L1B product from baseline D. This product also contains 20Hz measurements of, e.g., waveforms, tracker range, satellite altitude and corrections. The same retracker was applied to the C2 waveforms.

The water surface elevations are derived from the following relation (Equation (9))

$$H = h - R - N\_\prime \tag{9}$$

where *h* is the altitude of the satellite, and *N* is the geoid height, which here is the EST-GEOID2017 national geoid model. *R* is the range, i.e., the distance from the satellite to the surface which can be expressed as (Equation (10))

$$R = R\_{tracker} + R\_{retrack} + R\_{atm} + R\_{geo}.\tag{10}$$

Here *Rtracker* is the distance to the presumed surface measured by the onboard tracker and *Rretrack* is the retracking correction. *Ratm* represents the atmospheric corrections consisting of the wet and dry troposphere and ionosphere. *Rgeo* represents the geophysical corrections including solid earth tide, pole tide, and ocean loading tide. All corrections are here taken from the altimetry products.

The altimetry-based WL time series based for Peipsi and Võrtsjärv were reconstructed using the "R" package "tsHydro" available from https://github.com/cavios/tshydro (accessed on 14 October 2020). The software package is based on an implementation of a simple state-space model where the observations are assumed to be described by a random walk plus noise. The observation noise is described by a mixture of a Normal and the heavy tailed Cauchy distributions making it more robust against erroneous observation. In the solution, bias between the different missions is also estimated. The model is described in detail in [10]. Here, one time series based on all available altimetry data is constructed for each lake.

#### *2.4. Data Correction for the Effects of Water Level*

EO-derived WQ variables (chl a, FBM, SDD) were corrected for the effects of changing WL by the methodology developed by Tuvikene et al. [62]. The methodology is based on the long-term in situ data, where statistically significant relationships (*p*-value < 0.05) between the specific WQ variables and WL were found (Table A1). Corrected WQ variables were calculated using regression equations according to the long-term monthly mean WL and daily measured WL. If the WL of the sampling day corresponded to the long-term average for that period, no correction was needed, but if it differed, the correction either up or down was proportional to the deviation of the WL from the mean.

#### **3. Results**

#### *3.1. Inter-Annual Dynamics in Water Level*

Comparison of long-term WL measured in three hydrometric stations (Mustvee, Mehikoorma, Praaga) in Peipsi (Figure 2) show similar dynamics and very small spatial differences, e.g., the mean difference between Praaga and Mustvee is 0.02 m, and between Mehikoorma and Mustvee is 0.04 m. In both lakes, the higher and lower WL periods vary similarly and long-term monthly mean WL is highest in May and decreases towards October. While the annual difference can exceed the long-term monthly mean WL (Figure 2) in some years more than 1 meter, its effect is more pronounced in shallower Võrtsjärv (average depth 2.8 m) compared to Peipsi (average depth 7.1 m).

Water levels based on altimetry data and WL measured at the Mustvee hydrometric station show very similar dynamics (Figure 3) with a bias of 0.4 m. Bias corrected Root Mean Standard Error (RMSE) is 0.07 m, whereas the Pearson correlation between WL measured in Mustvee hydrometric station, and derived from altimetry data, is 0.99. Comparison between altimetry data and Rannu hydrometric station in Võrtsjärv shows similar results, with a bias 0.4 m. Bias corrected RMSE is 0.08, whereas the Pearson correlation between WL measured in Rannu hydrometric station and WL, derived from altimetry data, is 0.98. The combination of altimetry data from C2 and S3 missions allows to obtain much more frequent time series, especially since 2018 due to the launch of S3B (Figure 3).

**Figure 2.** Long-term water level (WL) in Rannu, Mustvee, Praaga, and Mehikoorma hydrometric stations from May (5) to October (10). Long-term monthly mean WL is showed in solid line.

#### *3.2. Inter-Annual Dynamics of Water Quality Parameters*

Long-term in situ time series of chl a, FBM, and SDD show variation in seasonal dynamics and pronounced inter-annual differences in Võrtsjärv (Figure 4A). For example, chl a shows increasing trend from 1995 until 2008, then it started to decrease again (Figure 4A). FBM shows increasing trend from the early 2000′ s. SDD shows seasonal variation from less than 0.5 to 1.5 m. As WQ parameters are in situ measured once during the vegetation period in Võrtsjärv station named as V2 (Figure 4B), the seasonal dynamics are better obtained with more frequent satellite data resulting on average 50 points, which doubled with S3B launch. This allows to acquire more representative data to analyze the WQ parameters over specific time period and have better knowledge of the seasonal and inter-annual differences.

The dynamics of satellite data corresponded well to in situ measured WQ parameters, showing similar dynamics during vegetation period with a bias for chl a 4.9 <sup>±</sup> 11.8 mg/m<sup>3</sup> , for FBM 2.7 <sup>±</sup> 14.2 g/m<sup>3</sup> , and for SDD was 0.3 ± 0.8 m. Chl a derived from OLCI corresponds to in situ data with Pearson correlation 0.7, whereas FBM 0.6 and SDD 0.4.

Depending on the monitoring system, some areas are not covered with sufficient monitoring frequency (Figure A1 Pihkva A) or include gaps (Figure A1 Peipsi s.s A., Figure A1 Lämmijärv A) in its long-term time series. Satellite data covers these areas 2–3 times per week and with two satellites the amount of data is doubled from 2018.

**Figure 3.** Water level (WL) measured in Mustvee and Rannu hydrometric stations from May to October. WL derived from altimetry data (S3A, S3B, C2).

**Figure 4.** Long-term time series of in situ measured chl a, FBM, and SDD in Võrtsjärv in station named as V2 (**A**) and the same parameters measured from S3 Ocean and Land Color Instrument (OLCI) (**B**). Dots indicate in situ measured values. Triangles and squares indicate S3A and S3B values, respectively. Different colors represent ecological status classes according to European Union (EU) WFD, while FBM is not considered as important parameter to assess the ecological status in Võrtsjärv according to the WFD.

#### *3.3. Correction of chl a, FBM, and SDD Based on the Water Evel*

Relationships between WL and WQ parameters vary in both lakes (Table A1). In Võrtsjärv, chl a was significantly related to WL (*p* < 0.05) in July, August, and October, FBM was significantly related in July, whereas SDD was related in all months (from May to October). In Peipsi *s.s.*, WL was significantly related only with FBM in June and September and in Pihkva, FBM was related in August. In Lämmijärv, chl a was significantly related in July and August, FBM was related in June and SDD in June, July, and August. When the relationship was significant between WL and specific WQ parameter, correction was applied on a daily basis according to the WL value in respect to long-term monthly mean WL. When WL was higher than long-term monthly mean WL, WQ parameter values increased after correction, if the WL was lower, then WQ parameters were corrected to lower values. This changed the mean value of the month, as well as minimum and maximum values.

In Võrtsjärv, WL was higher than long-term monthly mean WL in 2016 (up to 0.3 m) and lower in 2017 until 2019, especially low in 2018 and 2019 towards autumn, when WL was lower by more than 0.5 m (Figure 5). As WL was high in 2016, the mean value of chl a and SDD increased after applying the WL correction, whereas, in 2017 until 2019, mean values decreased. Due to the decrease in mean values, ecological status class of WQ parameters also changed. In 2017 and 2018, the ecological status class of chl a changed from "moderate" to "good", and in 2019 from "bad" to "moderate", whereas the ecological status class remained the same in 2016, although mean value of chl a increased. The ecological status class of SDD, changed from "good" to "very good" in 2018 and 2019, when the WL was the lowest compared to long-term monthly mean WL (Figure 5).

In Peipsi, there was less significant relationships between WL and WQ parameters than in Võrtsjärv (Table A1). In Peipsi *s.s*., there was no statistically significant relationships between WL and WQ parameters, except for FBM in June and September, which only leads to a small change of mean values of FBM (Figure A2 Peipsi *s.s.* FBM). As WL was lower than long-term monthly mean WL in 2016, 2018, and 2019, mean FBM value slightly decreased without changing the ecological status class. More significant relationships between WL and WQ parameters appeared in Lämmijärv (Figure A2 Lämmijärv), which also changed the ecological status class of SDD. As WL was more than 0.4 m lower than long-term monthly mean WL in 2019, then ecological status class of SDD changed from "moderate" to "good", and the mean value of SDD in 2016 also decreased. In Pihkva, the amount of in situ data was mostly insufficient for statistical analysis (significant relationship was found only in August for FBM) and therefore, there were no changes in ecological status class estimations of WQ parameters (Figure A2 Pihkva).

As Võrtsjärv and Lämmijärv are shallower than Peipsi *s.s.*, statistically significant relationships between WL and specific WQ parameters tend to occur more frequently there on a monthly basis. More pronounced changes in WL (difference approximately from 0.4 m) from long-term monthly mean WL were influencing not only the mean values of WQ parameters, but even changed the ecological status class from one to another. Smaller changes in WL in respective to long-term monthly mean influenced only the mean values of WQ, not the ecological status class. However, it is still important, considering long-term monitoring of the lake.

**Figure 5.** Water level change from long-term monthly mean WL during study period in Võrtsjärv (**1**). Value 0 indicates long-term monthly mean WL. Mean, maximum and minimum values of EO-derived WQ parameters are represented on the graphs (**2**–**4**), where dark gray indicates values before and light gray indicates values after applying WL correction. Different colors represent ecological status classes according to EU WFD.

#### **4. Discussion**

Long-term time series of WQ parameters in lakes are important to obtain information about the seasonal dynamics, their past and present situation and form the basis for future predictions [77,78]. The variety of measured parameters allow the partition of the influences by anthropogenic, as well as natural factors. Moreover, long-term time series help to understand the necessity for current management improvements for the future [79]. Water level as a natural factor, has an impact on lakes ecosystem and on WQ parameters, which in turn have a high impact on the ecological status of lakes [80,81]. Changes in the lake ecosystem influence biota existence and survival complexity in the lake, which affect ecosystem structure and future perspective [82–84]. Water level plays an important role in the lake ecosystem, especially in large shallow lakes [85], where monthly WL differences can be highly different from long-term observations, therefore removing WL influence from WQ parameters is beneficial for fair status assessment [62]. This study shows the possibility to create synergy between optical and altimetry data for lakes research, especially in terms of WFD reporting needs. The removal of one natural aspect will improve the estimation of the ecological status class of a lake, which will allow to keep the focus on the anthropogenic changes in the lake ecosystem.

The study shows high variations between seasonal and long-term WL in large and shallow lakes. Low WL periods vary with higher WL periods in both Võrtsjärv and Peipsi, which are caused by variation in precipitation, hydrological inflows, and evaporation in different years [86]. Long-term monthly mean WL decrease, from May towards October, characterizing higher WL after snow and ice melting, and lower WL in summer due to less inflow and higher rate of evaporation during the warm period [86,87]. Not only seasonal temperature and precipitation variability influences WL, but climate change has also direct effect to lake ecosystem [88–91], which has been even pointed out by The Intergovernmental Panel on Climate Change [92] referring to a global problem. The decrease in WL in some areas is often not compensated by the precipitation, which is caused by the higher annual mean air temperature, which, in turn, causes intensive evaporation and higher water temperature, causing a decrease in ice cover thickness and duration [93,94]. On the other hand, extreme events (flooding, rainfalls, storms, etc.) due to climate change

are causing an abrupt increase in the WL in some areas, causing also more apparent WL fluctuations [81,86].

As WL is sensitive to external factors, it is important to monitor it continuously [86,90,95–97]. As the number of in situ gauging stations is decreasing due to their cost and maintenance needs, the use of altimetry data becomes an important supplement to measure lake level variations. The use of altimetry has already provided promising results over lakes [10,11,16,50,98–105]. This is confirmed in this study, where a good agreement between the altimetry-based WL and WL measured at hydrometric stations with a correlation up to 0.99 in both lakes, with RMSE for Võrtsjärv 0.07 and for Peipsi 0.08, respectively. By applying measurements from different altimetry satellites the amount of data is increased, hence, providing a better understanding of WL dynamics and its changes [106–108]. However, combining data located at different tracks may introduce errors due to unresolved geoid signals, or other effects that causes the lake surface so be spatially different such as wind or land signals near the shoreline, but could be combined when the spatial differences in WL are small inside the water body. For example, comparing WL in Peipsi hydrometric stations, seasonal dynamics of the WL were similar and differences remained between 0.02 and 0.04 m, which shows possible WL homogeneity in the water body. Here, S3 SRAL and C2 SIRAL data were combined to increase the data frequency, but other missions could also have been used. Additionally, S3 and C2 provide SAR mode data in the study area, which give higher along-track resolution than conventional altimetry, with dense C2 spatial tracks suitable for lakes. The continuation of S3 missions, and the recently launched S6 mission, ensure data for future monitoring, which allows to continue developing EO-based applications to monitor WL and its implication for ecosystems and their services.

The possibility to derive certain WQ parameters from satellites helps to collect data more frequently and from places, which are not easily accessible. Satellite data corresponded to in situ data showing correlation 0.7 for chl a and lower than 0.6 for FBM and SDD, whereas bias for chl a was 4.9 <sup>±</sup> 11.8 mg/m<sup>3</sup> , for FBM 2.7 <sup>±</sup> 14.2 g/m<sup>3</sup> , and for SDD 0.3 ± 0.8 m. More data from satellites gives an extra information to in situ data, which enables to detect seasonal dynamics and rapid changes in different parts of a lake, especially in terms of estimating chl a, FBM, SDD based ecological status class. In some cases, satellite data overestimated (Figure 4 Võrtsjärv B SDD, Figure A1 Lämmijärv B SDD) or underestimated (Figure A1 Peipsi *s.s.* B chl a) in comparison of in situ data. This can be due to spatial differences (point measurement vs. 300 m × 300 m pixel), measurement technique (integral sample vs upper water column captured by satellite), but still, both in situ and EO data show similar trends during the vegetation period and complement each other. Different WQ parameters have been used for estimating ecological status of lakes from satellite data by many authors already for a while [7,25,33,36,41,109], which confirms necessity of this work for WFD.

Based on the long-term in situ data, the study shows significant relationships (*p* < 0.05) between WL and WQ more frequently in Võrtsjärv than in Peipsi. It is caused by the fact, that Võrtsjärv is shallower than Peipsi and high seasonal WL amplitude affects its biota in greater extent. In addition, the relationships between WL and WQ parameters were more often significant in the shallower parts of Peipsi (e.g., in Lämmjärv). Correcting WQ parameters with WL results in the change of mean annual WQ parameter values and even in changes of their ecological status class. The impact was more pronounced in Võrtsjärv, where WL was much lower than long-term monthly mean WL (approximately from 0.5 m), which did not only decrease the mean value of WQ parameters, but even changed the ecological status estimation to a higher class. As WL has a high impact to WQ parameters and the whole ecosystem, it is encouraging to see that EO data can be used to remove WL effect from data of WQ estimates in large shallow lakes, thus it could be one of the research topics deserving more attention due to its applicability to other large and shallow lakes. Moreover, the variety of past, present, and future satellite missions allows to apply the demonstrated methodology to extend the temporal and also spatial extent of the results. This could be done by combining satellite sensors for WQ, e.g., Environmental Satellite

(ENVISAT) Medium Resolution Imaging Spectrometer (MERIS), Moderate Resolution Imaging Spectroradiometer (MODIS) to estimate WQ parameters [35,40,110–114], with Topex Poseidon and Jason mission satellites to estimate WL data [51,104,105,115,116]. Additionally, several databases like Hydroweb [105,117], Database for Hydrological Time Series of Inland Waters (DAHITI) [98,118], Global Reservoirs and Lakes Monitor database (G-REALM) of United States Department of Agriculture (USDA) [119] provide water level data of lakes, whereas Globolakes [120,121] and Climate Change Initiative (CCI) [122] provide data of WQ parameters from different satellites. However, Sentinel missions will provide data at least for 20 years, which gives a possibility to gather continuous data from now.

As sustainable lake ecosystems are essential in the present and in the future, it is important to ensure frequent monitoring system and to develop EO-based approaches to improve monitoring and assessment of the ecological status class of lakes, leading, finally, to the better management of lakes. In situ monitoring is still essential to gather data to monitor the ecosystem, but also it is basis for developing and improving new tools and approaches, while satellite data would support it with better spatial and temporal resolution for certain parameters. Synergy between altimetry and optical data, gives possibility to estimate WL, as well as WQ parameters for more precise analysis. Results of this study show that, WL can have a strong impact on WQ parameters in large shallow lakes, where amplitude in WL change is high. This effect can be corrected with the proposed method to improve the estimation of ecological status class of such lakes.

#### **5. Conclusions**

EU WFD obligates European countries to monitor WQ of lakes and estimate their ecological status according to reference condition. The ecological status of lakes is influenced by anthropogenic, as well as by natural factors. To decrease the effect of natural factors, WQ parameters should be corrected for WL in case of shallow lakes with changing WL for more accurate identification of human impact and estimation of ecological status. Traditional monitoring methods give valuable long-term time series, however often with constrained spatial and temporal coverage, while satellites provide more frequent data. Many water quality parameters (e.g., chl a, FBM, and SDD) and WL can be derived using EO based data. The development and validation of EO algorithms allow to create synergy between different types of satellite data which enables to analyze and assess the lake ecosystems and provide complementary information to improve the monitoring and management decisions.

Long-term in situ data show inter-annual and seasonal dynamics of indicators characterizing lake ecosystems. In large and shallow lakes such as Võrtsjärv and Peipsi, amplitude of seasonal WL changes is high (more than 1 m), which influences WQ parameters, especially in shallower Võrtsjärv. The results support previous studies [10,99,101,104] in the fact that WL can be estimated with high accuracy (bias 0.4 m, R 0.99) over lakes with EO data. Correcting WQ parameters with the WL data derived from satellite, the mean values of the parameters and resulting ecological status classes changed in some cases. The highest changes occurred in Võrtsjärv during the years with lower WL compared to the long-term monthly mean. In this case the correction of the values of WQ parameters changed the resulting ecological status to better quality class. Therefore, taking into account the variability of WL and correcting its effect in case of shallow lakes is important to decrease the influence of the natural factors on the assessment of the lakes ecological status class in terms of EU WFD.

**Author Contributions:** Conceptualization, A.A.-T. and K.A.; methodology, A.A.-T. and K.A.; validation, A.A.-T. and K.A.; formal analysis, A.A.-T. and K.A.; resources, A.A.-T., K.A., K.N., L.T., and K.K.; writing—original draft preparation, A.A.-T.; writing—review and editing, A.A.-T., K.A., K.N., L.T., and K.K.; visualization, A.A.-T.; supervision, K.A., K.N., and L.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by (1) Estonian Research Council research and development program grant number RITA1/02-52-01; (2) EU's Horizon 2020 research and innovation program (grant agreement no. 730066, EOMORES) and; (3) Estonian Research Council personal start-up grant named as PSG10 and team grant named as PRG709.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request.

**Acknowledgments:** We are grateful to anonymous reviewers for their feedback and suggestions that helped to improve the manuscript. This study was financially supported by (1) the European Regional Development Fund within National Program for Addressing Socio-Economic Challenges through R&D (project named as RITA1 grant number RITA1/02-52-01); (2) European Union's Horizon 2020 research and innovation program under grant agreement no. 730066; (3) Estonian Research Council grants PSG10 and PRG709. Authors thank people from Centre for Limnology providing in situ water quality data and Estonian Weather Service providing in situ water level data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

#### **Figure A1.** *Cont*.

**Figure A1.** Long-term time series of in situ measured chl a, FBM, and SDD in Peipsi *s.s.*, Lämmijärv, and Pihkva (**A**) and same parameters measured from S3 OLCI (**B**). Dots indicate in situ measured values. Triangles and squares indicates S3A and S3B values, respectively. Different colors represent ecological status classes according to EU WFD.

− − − − − − − − − − − − − − − − − − − − − − − − − − − − −

#### **Appendix B**

**Table A1.** Relationship between in situ measured WL (cm) and water quality parameters from May (5) to October (10), where bold values indicate statistically significant relationships (*p* < 0.05), N indicates number of observations during the study period and "-" indicates not enough observations for statistical analysis.




#### **Appendix C**

**Figure A2.** *Cont*.

− −

− − − − − − − − − −

**Figure A2.** Water level change from long-term monthly mean WL during the study period in Peipsi *s.s.*, Lämmijärv, and Pihkva (**1**). Value 0 indicates long-term monthly mean WL. Mean, maximum, and minimum values of water quality parameters are represented on the graphs (**2**–**4**), where dark gray indicates values before and light gray indicates values after applying WL correction. Different colors represent ecological classes according to EU WFD.

#### **References**

