*Article* **Monitoring Lakes Surface Water Velocity with SAR: A Feasibility Study on Lake Garda, Italy**

**Marina Amadori 1,2,† , Virginia Zamparelli 3,† , Giacomo De Carolis <sup>1</sup> , Gianfranco Fornaro <sup>3</sup> , Marco Toffolon <sup>2</sup> , Mariano Bresciani <sup>1</sup> , Claudia Giardino <sup>1</sup> and Francesca De Santi 1,\***


**Abstract:** The SAR Doppler frequencies are directly related to the motion of the scatterers in the illuminated area and have already been used in marine applications to monitor moving water surfaces. Here we investigate the possibility of retrieving surface water velocity from SAR Doppler analysis in medium-size lakes. ENVISAT images of the test site (Lake Garda) are processed and the Doppler Centroid Anomaly technique is adopted. The resulting surface velocity maps are compared with the outputs of a hydrodynamic model specifically validated for the case study. Thermal images from MODIS Terra are used in support of the modeling results. The surface velocity retrieved from SAR is found to overestimate the numerical results and the existence of a bias is investigated. In marine applications, such bias is traditionally removed through Geophysical Model Functions (GMFs) by ascribing it to a fully developed wind waves spectrum. We found that such an assumption is not supported in our case study, due to the small-scale variations of topography and wind. The role of wind intensity and duration on the results from SAR is evaluated, and the inclusion of lake bathymetry and the SAR backscatter gradient is recommended for the future development of GMFs suitable for lake environments.

**Keywords:** SAR; Doppler Centroid Anomaly; inland waters; physical limnology; hydrodynamics

#### **1. Introduction**

Microwave methodologies based on the use of Synthetic Aperture Radar (SAR) sensors have been widely developed for many Earth surface monitoring applications. In the last decade, images acquired by SAR have been increasingly exploited by the scientific community for deformation surveys [1] even at single facilities scale, that is, buildings and infrastructures [2]. More recently, increasing attention has been paid to the observation of sea surface, complementing the traditional use of optical or multispectral images.

Differently from optical sensors, SAR allows the obtaining of all weather and day– night 2D images of the illuminated area of the Earth's surface [3]. SAR remote sensing for ocean, seas and coastal applications mostly exploits the amplitude of the backscattered signal for, for example, monitoring oil-spills [4] and sea-ice [5], ship detection [6] and high-resolution wind fields retrieval [7]. However, SAR uses coherent radiation and the complementary information carried in the phase of the received complex signal can be also exploited. An example from land applications is the use of the complex information from a single SAR image to reconstruct infrastructures micro-motion [8]. By analysing the complex backscattered (received) signal, it is possible to measure the Doppler properties

**Citation:** Amadori, M.; Zamparelli, V.; De Carolis, G.; Fornaro, G.; Toffolon, M.; Bresciani, M.; Giardino, C.; De Santi, F. Monitoring Lakes Surface Water Velocity with SAR: A Feasibility Study on Lake Garda, Italy. *Remote Sens.* **2021**, *13*, 2293. https://doi.org/10.3390/ rs13122293

Academic Editor: Joong-Sun Won

Received: 30 March 2021 Accepted: 7 June 2021 Published: 11 June 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/).

of the scatterers. The latter are directly related to the motion of the scatterers along the radar line of sight (LOS) in both land [9] and marine [10] contexts. In ocean applications, these properties are currently used to retrieve near surface wind speed [11] and surface current [12].

Several factors contribute to the measured Doppler frequency, such as the bulk surface current, the wind waves and the wave–current interactions at different scales [13]. The drift current induced by the wind on the sea surface (wind drift) and the sea state (more precisely the high-frequency waves) depend on the near-surface wind field. In this regard, [14] performed a global ocean analysis and and found that an increase of the wind speed (projected along the LOS) corresponds to a first-order increase of the Doppler anomaly. Such a wind-wave bias is empirically estimated and removed to retrieve the radial component of the bulk surface current. To this end, empirical Geophysical Model Functions (GMF) have been determined [12].

This analysis is going toward a consolidated state for oceanic applications, such that surface velocity maps in open waters are distributed by the European Space Agency (ESA) among the second level products of Sentinel 1 [15]. On the contrary, the extension of a similar approach to coastal zones is more delicate. At least two main assumptions of the analysis of the Doppler Centroid Anomaly (DCA) in the open ocean do not hold for coastal areas: (i) uniform surface flow field and (ii) fully developed wind-waves. In the near-shore areas, small- scale variations of topography and wind can cause higher variability in the currents and the waves to develop in fetch-limited condition.Thus, the application to other contexts of the above mentioned methods for the removal of the wind waves bias may not be straightforward as in the open ocean case, although examples along this line are given in the recent literature [10,13,14,16–19]. At the smaller spatial scale of enclosed basins (e.g., lakes), no studies have yet explored the possibility of inferring surface currents from SAR Doppler frequencies.

More generally, SAR images of lakes have been seldom exploited for several reasons, among which the minor intensity of hydrodynamic processes (e.g., smaller wave heights and currents than in open seas and oceans) and the effect of surrounding topography (e.g., foreshortening and layover effect [20]).

Currently, SAR applications for lakes are limited to the analysis of the backscattered signal amplitude. SAR amplitude images have been used in support of optical or thermal imagery for identifying and mapping inland waters surface features, for example, ice [21] and cyanobacteria blooms [22,23]. A few tests also showed the potentialities of using SAR images for retrieving the size of surface eddies in large lakes [24] and spatially distributed information on wind intensity [25].

The SAR Doppler analysis could represent a useful way to retrieve surface velocity in lakes and would compensate for the lack of synoptic measurements of the key physical quantities determining lake hydrodynamics. Indeed, while the use of remote sensing products is well consolidated to infer lake surface water temperature (LSWT) [26], water velocity, as well as all wind velocity, is traditionally monitored with in-situ point measurements. In most cases, traditional instrumentation requires direct access to the lake, with strong operational limitations, and provide data at single locations. Remote sensing techniques provide instead synoptic coverage, fine spatial detail and repeated regular sampling.

Periodic snapshots of lake surface currents, at the present day, cannot be achieved alternatively than with remote sensing instruments, for example, airborne or space, the latter being traditionally exploited in marine applications [27–29]. Few attempts towards the reconstruction of the surface transport from thermal infrared imagery [30] showed encouraging results in this direction, but such approach requires high temporal resolution, which is not always guaranteed. Moreover, the use of thermal infrared remote sensing is hindered by cloud cover. In this regard, the use of SAR images combined with the DCA technique could overcome the technical issues typically encountered for reconstructing the surface flow field in lakes, asides from being rather innovative.

This study aims to investigate the feasibility of extracting the surface velocity from DCA in lakes by considering as test case Lake Garda, a large and deep lake in northern Italy. The test case is a clear oligotrophic lake, thus allowing for neglecting turbidity and algal blooms which might affect the SAR signal.

We analyze and process images acquired by the ENVISAT C-Band sensor. Results are compared with the outputs of a three-dimensional (3D) numerical model validated for Lake Garda [31]. Water surface temperature products from the MODIS (Moderate Resolution Imaging Spectroradiometer) Terra sensor are also considered to further validate the simulation results. The numerical results are used to interpret the SAR retrieved signal and to evaluate the different factors influencing the Doppler shift analysis in the case study. Based on this, we draw some guidelines on the variables to be considered for a future definition of GMFs for lakes.

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

In this section we describe the test site (Section 2.1), the sensors and datasets for SAR and thermal infrared images (Section 2.2) , the procedure for the analysis of the DCA (Section 2.3) and the adopted numerical models (Section 2.4).

#### *2.1. Case Study*

Lake Garda (Figure 1a) is one of the large European perialpine lakes and is the largest lake in Italy by surface extension and volume. It covers an area of 370 km<sup>2</sup> , with a water volume of 50 km<sup>3</sup> and a maximum depth of 346 m. The lake represents an essential water supply for many sectors of the local economy (e.g., agriculture, industry, fishing and drinking [32]). Thanks to its attractive landscape, mild climate and water quality, Lake Garda is also an important resource for recreation and tourism.

**Figure 1.** (**a**) Location of Lake Garda in northern Italy; orography of the surrounding area in m above sea level (m a.s.l.); lake bathymetry in m below water surface (m b.w.s.). (**b**) Scheme of the SAR acquisition geometry. (**c**) Example of an amplitude focused SAR image of an ENVISAT acquisition centered on lake. The SAR image is geocoded and superimposed on the corresponding Google Earth image.

The lake has a heterogeneous bathymetry (Figure 1a). A narrow (average width 4 km) and deep (maximum depth 346 m) northern trunk, enclosed between steep mountains, is connected to a southern large (maximum width 18 km) and shallower (average depth 65 m) basin, which is characterized by more gentle slopes and lies in a flat plain.

The typical winds in Lake Garda are thermal breezes [33]. In the northern region, the most frequent provenance direction of high speed winds is along the longitudinal axis of the lake, where winds are channeled by the steep lateral mountains (see the orography in Figure 1a). During summer, northerly and southerly breezes alternate in the night/morning and in the afternoon/evening, respectively [34]. Such breezes are characterized by moderate to high velocities (≈5–10 m/s), alternating directions and not necessarily uniform spatial distribution. During winter and mid-seasons, winds are usually weaker. However, it is not uncommon for this lake to be subject to long-lasting storm winds at a synoptic scale, for example, Föhn winds. These winds frequently come from North-East, have almost uniform spatial distribution on the lake surface, reach speeds of more than 10 m/s and often last for more than one day [35].

#### *2.2. Sensor and Dataset*

In the last years, there has been a consistent proliferation of SAR systems with different technical characteristics, that is, resolution, operational frequencies and revisiting time [36]. Among them we selected the ENVISAT, following the experience gained in marine environment with C-Band sensors for the estimation of the surface current velocity. In this work, we adopt a methodology validated in the coastal area of the Gulf of Naples [16,27,37] and the Gulf of Trieste [38]. ASAR ENVISAT operated in C Band (wavelength *λ* equal to 5.6 cm), in multiple polarization mode (HH/VV,HH/HV,VV/VH). In stripmap mode it covered up to 100 km with a single swath. The spatial resolution of the sensor's Single Look Complex (SLC) product is 5 m in azimuth and 20 m ground range (see Figure 1b for SAR acquisition geometry). A considerable archive of data exists from the ENVISAT operative period, that is, from 2002 to 2012, on ascending and descending orbits. Within the broad landscape of SAR sensors currently available, data acquired by Sentinel 1 (S1) could also be exploited. S1 operates in C-Band and regularly provides very wide swaths of about 250 km by using the TOPSAR (Terrain Observation with Progressive Scans SAR) operation mode. However, due to the steering of the antenna along the azimuth, this operation mode complicates the estimation of the Doppler Centroid. The presence of discontinuities in the estimated Doppler Anomaly, which in turn affects the retrieved radial velocity, has been already observed [39]. Based on these considerations and on the experience we gained in marine environments, ENVISAT data were chosen for this exploratory study on Lake Garda. The ENVISAT archive was screened to eliminate images characterized by very low backscattering coefficients. We selected six images showing amplitude patterns, that is, sufficiently intense and spatially varying backscattering, potentially associated with interesting hydrodynamic features and fulfilling the wind criteria specified in detail in Section 3.2. The selected images are acquired over ascending orbit, corresponding to evening acquisitions (around 21:00 UTC), and span from 2008 to 2010 in summer, autumn and winter seasons. The choice of ascending acquisitions is due to the favourable acquisition geometry related to the orientation of the lake. In fact, the ground projection of the ascending sensor LOS is more aligned with the longitudinal axis of the lake than the descending one. We note that the longitudinal axis of the lake is expected to be the main direction of lake currents and waves propagation, as it also corresponds with the typical direction of winds (see Section 2.1). The images have been classified based on the wind velocity registered at the acquisition time. In Table 1 the selected dates are listed and their simulated geophysical characteristics are summarized. For a more detailed description of the investigated dates, we refer to Section 3.2.

For the selected dates, thermal images from MODIS were also analyzed to validate the simulations. In particular, we preferred nighttime Level 2 MODIS Terra Sea Surface Temperature (SST) products, because they are closer in time to SAR acquisitions. MODIS

images are cloud free for four out of six investigated dates. The four dates are indicated with an asterisk in Table 1.

**Table 1.** Mean, standard deviation and maximum of the simulated key geophysical quantities for each date considered. Module of wind speed from WRF (*W*SIM); module of the surface water velocity from Delft3D + SWAN (*v* SIM); significant wave height (*H*s) and peak wavelength (*λ*peak) from SWAN.


∗ MODIS product available.

#### *2.3. SAR Doppler Centroid Anomaly Estimation*

For the retrieval of information about the ocean surface velocity from SAR data, two methodologies have been developed: the Along-Track Interferometry (ATI) [18] and the DCA method [14]. Here we exploit the technique of the DCA, which was successfully applied in open oceans [14,40,41] and ocean coastal areas [14,18,42]. The same approach was also adopted for the estimation of water currents in the Mediterranean sea, a semienclosed sea, where they show a smaller range of intensity, [16,27,38] . The DCA technique estimates the surface velocity of the illuminated area by measuring the variation that such a velocity produces in the spectrum of the SAR image [14]. The principle is rather well known: let us consider a SAR system acquiring data on a scene where a target is in relative motion with respect to the sensor. When the transmitted radiation hits the target, it is backscattered toward the radar along the antenna pointing direction, also called Line of Sight (LOS) or radial direction (see Figure 1b). If the velocity associated with the motion has a non-zero component along the radial direction (*v*r), this component determines a shift in the azimuth spectrum of the received signal. Such a shift is called Doppler Centroid frequency *fDC* and is proportional to the radial component of the velocity associated to the relative motion between target and sensor.

Some remarks on the formation of the Doppler shift, which is influenced by the target and sensor motion along the LOS direction, are in order. First of all, the radial direction is influenced by the platform attitude as well as by possible electronic steering introduced to point the beam in a direction different from the broadside one (squinted acquisition). Secondly, the measurement of *fDC* from the power spectral density allows the estimation of the radial component of the (relative) target motion [43,44]. For each pixel, say with azimuth *x* and slant range *r*, *fDC*(*x*,*r*) can be considered the sum of two contributions: the Doppler Centroid *fDC*0(*x*,*r*) and Doppler Centroid Anomaly *fDCA*(*x*,*r*). The first one (*fDC*<sup>0</sup> ) corresponds to a stationary scene and is referred to as "background DC", due to a scene in which all the scatterers move at the same velocity (background motion). The second one (*fDCA*) is the Doppler component associated with the target motion with respect to the background one.

From what is stated above, it follows that for the satellite case the background DC, *fDC*0(*x*,*r*), is contributed by both the platform movement and the Earth rotation, in a way that depends on the platform attitude and beam orientation. This term is generally a slow varying function along the space, typically well approximated by a polynomial in *r*. It is worth noting that, strictly speaking, ionospheric variations could also impact *fDC*0(*x*,*r*): this phenomenon is, however, more relevant for low frequency systems, for example, Land P-Band SAR sensors [45]. The latter are usually not exploited for water surface velocity retrieval, due to the reduced sensitivity associated with the increase of the wavelength.

As for ASAR-Envisat, the term *fDC*0(*x*,*r*) can be retrieved from the ancillary data through annotated polynomial coefficients, which are computed knowing the platform attitude, the antenna pointing and the scene latitude. Previous experiences carried out with other platforms showed that this information is frequently not sufficiently accurate, especially for the range variation. This fact might be due to inaccurate attitude information, to unexpected beam bending as well as to other spurious effects. Following the strategy adopted in [16] for marine applications, we estimated *fDC*0(*x*,*r*) directly from the data via a polynomial Least Square fitting along *r*. More details on the estimation of the background Doppler component associated with the stationary scene are given in [16]. In the same work it is also shown how to handle unwanted azimuth oscillating components, frequently observed for ASAR-ENVISAT acquisitions and evidently not associated with any geophysical process.

The velocity at which the target proceeds toward or away from the radar, with respect to the background scene motion, can be thus obtained by removing from the measured Doppler Centroid its component associated to the stationary scene, as follows:

$$\begin{split} v\_{\mathbf{r}}(\mathbf{x},r) &= \frac{\lambda}{2} [f\_{\text{DC}}(\mathbf{x},r) - f\_{\text{DC0}}(\mathbf{x},r)] \\ &= \frac{\lambda}{2} f\_{\text{DCA}}(\mathbf{x},r), \end{split} \tag{1}$$

where *λ* is the sensor wavelength. The scatters' surface velocity along the ground range direction *v*gr can be finally obtained from *v*<sup>r</sup> by considering that the latter is its projection along the LOS (see also Figure 1b):

$$v\_{\rm gr} = \frac{v\_{\rm r}}{\sin \theta'} \tag{2}$$

where *θ* is the incidence angle.

In order to retrieve the actual lake ground range surface velocity, any possible spurious effect, such as the wind-wave bias [12], must be precisely identified, estimated and removed from *v*gr. Thus, *v*gr represents a raw variable to be processed with an appropriate GMF to obtain the ground range surface velocity, hereinafter surface velocity. The identification of the spurious effects on *v*gr is the goal of this study, where no GMFs are used, and it is performed with the help of the numerical model.

In Figure 2a ,we summarize the procedure for the generation of *v*gr maps. After the SAR focusing operation on the raw data, the measurement of *fDC* is performed pixel by pixel by estimating the azimuth self-correlation function along the azimuth direction [43]. The Doppler spectrum is evaluated on patches, hereafter referred to as Doppler Resolution Cell (DRC), sliding along the whole SAR image. The DRC has been chosen to guarantee a good balance between spatial and spectral resolutions. In particular, we considered patches of width 512 azimuth by 128 range pixels, corresponding to a DRC of 2.5 km × 2.5 km. In this condition the spectral resolution is 2.85 Hz, which corresponds to a velocity resolution of 7.98 cm/s for *v*r.The Doppler Centroid estimation accuracy can be estimated through the calibration of the radial velocity over the land area (i.e masking out Lake Garda). For all the acquisitions here analyzed (see Table 1) *v*<sup>r</sup> has zero mean and a standard deviation ranging from 27 cm/s to 30 cm/s. Note that these values account for the variability of the residual signal all over the scene with the lake only covering a small portion.

Subsequently, the compensation of the stationary component of the Centroid is performed and the final DCA is estimated (Equation (1) and then *v*gr from Equation (2)). The DCA map is converted to surface velocity and geocoded through a SRTM DEM on a geographical grid with a spacing of about 92.7 m N and 64.9 m E [46]. At this stage, the geocoded surface velocity map is comparable with geo-referred simulations results and MODIS products.

The output product is supplemented by the basckscattering information. A standard calibration procedure is applied to the focused amplitude data to get the normalised measure of the radar cross-section per unit area on the ground, that is, the backscatter coefficient *σ* 0 . A multi-look algorithm is then applied to reduce the noise (or "speckle") of SAR images by averaging adjacent pixels. The resulting *σ* <sup>0</sup> with resolution 100 m <sup>×</sup> 100 m is finally geo-referred (Figure 1c).

**Figure 2.** (**a**) Procedure for the generation of *v*gr and *σ* <sup>0</sup> maps from SAR data; (**b**) modeling chain for the simulation of surface water velocity, temperature and waves field.

#### *2.4. Numerical Modeling*

For simulating the flow field, we use a modeling chain (Figure 2b) composed by three numerical models: an atmospheric, a hydrodynamic and a wave model. In this work, we start from the outputs of an existing model for the case study [35,47] in the configuration presented in [31], which consists of the hydrodynamic model Delft3D [48] fed by the results of the Weather Research and Forecasting (WRF) atmospheric model [48].

The WRF model provides the atmospheric variables as time and space varying fields of wind velocity, air pressure, temperature and relative humidity, incoming shortwave radiation and cloudiness at the lake surface. The resolution of the simulated weather forcing is 2 km in space and 1 h in time.

From the WRF outputs, the hydrodynamic model Delft3D computes the heat fluxes between air and water and wind stress at the lake surface, which are necessary for simulating the flow field, turbulence, heat and mass transport within the lake. The nominal horizontal resolution of the hydrodynamic model is ≈100 m, with smaller cells in the northern and larger cells in the southern part. The vertical layers are 1 m thick near the surface (over the upper 40 m) and smoothly increase their thickness up to 25 m at the deepest point. The outputs of the described model are available over the period 2004–2018 and will be henceforth referred to as "long-run data". The long-run data are saved on a daily basis at 10:00 UTC.

For the present work, the available long-run data were used to restart ad hoc simulations of the selected dates (Table 1) with the same setup as in [31]. The restart was necessary to produce outputs synchronous with the ENVISAT acquisitions (21:00 UTC) and to couple the Delft3D model with the wave model. We adopted a third-generation wind-wave model known as Simulating Waves Nearshore (SWAN) [49]. The wave spectrum employed in the simulations was composed of 101 logarithmically spaced frequencies in the range of 0.15–3 Hz. The frequency range was specifically chosen to be representative of wave field expected in Lake Garda and therefore included higher frequencies than those typically employed in coastal applications, but neglected very low frequencies that hardly occur in medium-size lakes. Since no data were available for the calibration of the model, the parameters were chosen in consistency with previous experiences in other large and deep lakes (e.g., Lake Constance [50] and Michigan [51]).

The results of the six Delft3D + SWAN simulations were saved at hourly intervals, over the whole computational domain, and covered a period of two days (the SAR image acquisition day plus one day before).

#### **3. Results**

In this work we use the model outputs with the aim of interpreting the results obtained from the SAR images. Thus, we rely here on the assumption that the model is capable of reproducing wind fields and surface transport. The reliability of the wind fields simulated by the WRF model has been extensively demonstrated by [31,52]. As for the surface transport, we base our work on [31], who showed that the Delft3D model results are coherent with in-situ and remotely sensed water temperature observations. We stress that [31] provided an indirect verification of the flow field, since in-situ measurements of water velocity are not available. In Section 3.1 we perform a similar verification by comparing the spatial patterns of surface temperature simulated by the model and reconstructed from the MODIS images. In Section 3.2, we introduce the main features of the investigated dates and in Sections 3.3 and 3.4, we evaluate the outputs of our procedure for the extraction of lake surface currents from SAR.

#### *3.1. Model Verification Against MODIS*

We provide here the comparison between the Level 2 MODIS Terra SST products (*T* MODIS) in Celsius degrees and the simulated lake surface water temperature (*T* SIM). Figure 3 shows that the model correctly reproduces the remotely sensed patterns of lake surface temperature and the spatial variability of this quantity in the different seasons. This confirms that the heat fluxes between atmosphere and lake are well simulated, and that the numerical flow field is coherent with the transport patterns responsible for horizontal anomalies in lake surface temperature. The difference between the spatial mean values of *T* SIM and *T* MODIS (reported in each panel of Figure 3) can be due to model errors as well as to the skin-bulk bias [53]. However, in all investigated dates the mean error is below 1 ◦C, well in line with the performances of the long-run model in [31].

**Figure 3.** Top plots: maps of SST from MODIS Terra on cloud-free investigated dates. Bottom plots: maps of simulated lake surface water temperature. The spatial mean of MODIS and simulated water temperature is displayed together with the existing bias (*T*SIM − *<sup>T</sup>*MODIS). Note: in each panel, colorbar limits are set differently to highlight the spatial patterns. The colorbar limits of bottom plots are set as those from top plots at the corresponding date + the computed bias.

#### *3.2. Main Features of Investigated Dates*

The wind field has been recognized as playing a major role for the use of SAR in the open ocean [54]. In particular, in order generate capillary waves detectable from C-band sensors, the wind forces have to overcome the viscous ones. A wind speed threshold value for this process to occur is estimated to be about 3.25 m/s (at 10 m above the surface, [55]).

Therefore, we classify the six cases considered in this study (and listed in Table 1) based on the wind intensity. In Figure 4 the daily cycle of wind speed and direction simulated in a mid-lake point (point P in Figure 1) is displayed for the selected dates. The SAR acquisition time is indicated as a red vertical line in the figure panels.

We define "low wind" dates (Figure 4a–c) the three acquisitions in which the spatial mean of the simulated wind speed is equal or smaller than this threshold value (29 October 2009, 5 August 2010, 14 October 2010). In these dates mild evening breezes blew at 21:00 UTC over the lake surface (≈1–3 m/s). This condition is frequent in sunny and warm days, when thermal breezes develop with alternating directions during daylight hours and calm down after the sunset. Under low wind conditions, simulations predict mean surface velocities below 10 cm/s, significant wave heights *H*<sup>s</sup> of the order of few cm, and peak wavelengths *λ*peak of the order of 1 m (Table 1). Such values suggest that gravity effects do not completely dominate the dynamics as surface tension still plays a significant role [56].

**Figure 4.** Simulated wind direction and speed on the investigated dates a in-lake point (point P in Figure 1a). The red arrow in each panel corresponds to the wind simulated in P at the SAR acquisition time (i.e., the red line at 21:00 UTC). The sketch at the right displays the location of P and the direction of the wind with respect to the lake geography.

The remaining three acquisitions are classified as "high wind" dates (7 February 2008, 13 November 2008, 11 June 2009, Figure 4d–f). On those days, the mean value of wind speed is above 7 m/s and its maximum value exceeds 15 m/s in large regions of the lake. Both the mean and the maximum wind speed values are significantly higher than the above-mentioned threshold. Such high speeds are often due to large scale synoptic winds, blowing almost uniformly over the lake's surface for more than 12 h. Figure 4d,e shows that the latter is indeed the case for the dates 7 February 2008 and 13 November 2008, when a northerly Föhn wind persistently blew for one day before the SAR acquisition time. On the contrary, the date 11 June 2009 (Figure 4f) coincides with the onset of a similar wind, which started one hour before the acquisition time and lasted all night long (not shown).

The simulated gravity waves on high wind speed days have *H*<sup>s</sup> of the order of 50 cm on average (Table 1), overcoming 1 m in some areas of the lake, and *λ*peak of the order of 10 m. These values are fully consistent with the gravity waves field that can be theoretically expected [57] and observed [58] in a medium-size lake.

In the next two sections, we will provide for all investigated dates the amplitude SAR backscatter image (*σ*0), the ground range surface velocity retrieved from SAR DCA (Figure 2a, *v*gr) and simulated by the numerical model (Figure 2b, *v* SIM gr ), and the simulated wind field (direction and ground range component magnitude *W*SIM gr ). The represented quantities are obtained as follows:


#### *3.3. Low Wind Dates*

The results for the "low wind" dates are reported in Figure 5. In marine/ocean observations, a low-wind regime corresponds to a low mean amplitude and a limited spatial variation of the backscatter. In our case study, Figure 5 clearly shows a different behaviour. The SAR backscatter retrieved on these dates (first column of the Figure 5, panels a, e, i) shows recurring patterns in all examined images. Higher values are detected in the northwestern basin, lower values in the south-eastern region. A high intensity feature located at the NE edge of the lake is present in all the ascending SAR images analysed. It consists in a cross-track radiometric compression, known as foreshortening effect, which occurs when the radar beam impinges on the foreslope of a mountainous area [20]. The component of surface motion retrieved from SAR DCA *v*gr is represented in the second column of Figure 5, panels b, f, j. Red and blue colors (redshift and blueshift) indicate a motion towards and away from the radar antenna, respectively. We recall here that the surface velocity has been estimated with a DRC of 2.5 km × 2.5 km (see Section 2.3) and that any variation of the resulting *v*gr at a smaller spatial scale has to be interpreted with care. While this is not a critical issue along the longitudinal axis of the lake and in the southeastern basin, where both the length and the width of the lake are significantly larger than the DRC size, difficulties arise in the northern trunk of the lake. Here the DRC size is comparable to the width of the lake and the estimation might be affected by the presence of portions of land inside the DRC. For such a reason we recommend caution in the interpretation of the variations of *v*gr along the crosswise direction.

The spatial patterns of *v*gr resemble those from *σ* 0 . In the northern trunk of the lake, *v*gr shows in all dates a recurring pattern: it overcomes 2 m/s in the central and western region, and is almost null along the eastern shore. This behaviour suggests a south-westward surface transport which is not confirmed by the numerical model (third column of Figure 5, panels c, g, k). In fact, the simulated field of *v* SIM gr shows significant differences from one date to another, mostly due to the differences in the spatial distribution of the wind field (fourth column of Figure 5, panels d, h, l). Moreover, the simulations show significantly lower surface velocities, below 0.1 m/s, as it is reasonably expected in a lake of the dimensions of Lake Garda under low wind conditions [59].

In the southern sub-basin, similar patterns are found in *v*gr, retrieved on 29 October 2009 and 14 October 2010 (Figure 5b and j respectively), with high and positive velocities at the center of the sub-basin and almost null or negative velocities close to the shores. The spatial distribution of *v*gr in this region is similar to that of *v* SIM gr on the same dates (Figure 5c,k). On 5 August 2010 *v*gr (Figure 5f) is strongly positive over the entire southern region, while *v* SIM gr (Figure 5g) is strongly negative. The spatial distribution of *v* SIM gr on this date reveals a predominant eastward surface transport in the southern basin and an

opposite tendency in the northern basin, which is the direct consequence of two divergent wind fields with respect to the ground range direction (Figure 5h).

**Figure 5.** Maps from low wind days (**a**–**d**) 29 October 2009, (**e**–**h**) 5 August 2010, (**i**–**l**) 14 October 2010 of: first column (**a**,**e**,**i**) SAR bakcscatter amplitude; second column (**b**,**f**,**j**) *v*gr estimated from SAR DCA (Equation (2)); third column (**c**,**g**,**k**) *v* SIM gr simulated by the numerical model; fourth column (**d**,**h**,**l**) ground range component of simulated wind speed *W*SIM gr . Black arrows in fourth column indicate the wind vector. The ground range direction is displayed by the blue and red arrow at the bottom-right corner of panels.

#### *3.4. High Wind Dates*

The results for the "high wind" dates are reported in Figure 6. The SAR backscatter images (Figure 6, panels a,e,i) and the surface velocity *v*gr (Figure 6b,f,j) present similar patterns, which are in good agreement with the simulated surface velocity *v* SIM gr (c,g,k), especially in the elongated region.

**Figure 6.** Maps from high wind days (**a**–**d**) 7 February 2008, (**e**–**h**) 13 November 2008, (**i**–**l**) 11 June 2009 of: first column (**a**,**e**,**i**) SAR backscatter amplitude; second column (**b**,**f**,**j**) *v*gr estimated from SAR DCA (Equation (2)); third column (**c**,**g**,**k**) *v* SIM gr simulated by the numerical model; fourth column (**d**,**h**,**l**) ground range component of simulated wind speed *W*SIM gr . Black arrows in fourth column indicate the wind vector. The ground range direction is displayed by the blue and red arrow at the bottom-right corner of the panels.

On all dates *σ* 0 shows a pattern characterized by increasingly high values as moving south-west in the northwestern basin (Figure 6a,e,i). Such pattern can be observed also in the *v*gr images (Figure 6b,f,j), with a redshift located mostly along the western shore indicating a south-westerly surface current moving towards the sensor. The value of *v*gr in this area is higher than 2 m/s and slightly decreases in the southeastern part of the region. This behavior finds qualitative correspondence with the patterns of *v* SIM gr (Figure 6c,g,k), where the signature of wind waves can be detected from localized peaks of velocity up to 0.5 m/s. The black arrows displayed in Figure 6d,h,l indicate the presence of a wind blowing from north-east to south-west, which is responsible for the surface transport observed both from SAR and the model.

#### **4. Discussion**

In all cases examined, the surface ground range velocity retrieved from SAR (*v*gr) is much higher than the simulated *v* SIM gr . Regardless of the wind intensity, *v*gr always overcomes ±2 m/s, which is too high for a medium-size lake environment. Despite the empirical relations linking lake surface velocity and wind speed are often site-specific, the order of magnitude of the first quantity is generally assumed to depend on the second. A wind factor can be computed as the ratio between water surface velocity and wind speed. In open waters the wind factor is typically of the order of 3%, but it can also be lower for the nearshore areas [60]. In the investigated dates, the simulated water velocity was found to be, on average, between 2% and 3.5% of the wind speed. Even assuming that the observed *v*gr is the only component of the surface current (i.e., assuming that the water is moving precisely along the ground range direction), a value of 2 m/s would require a wind speed of more than 60 m/s, which is an absolutely unrealistic wind speed. Such wind speed would be even higher if we take into account also the other component of surface motion, which is likely to exist and contributes to the module of surface current, despite it can not be retrieved with the DCA technique. Hence, the velocity retrieved from SAR is a large overestimation of the actual lake surface velocity.

Such overestimation is also found in open ocean applications [61], where *v*gr is affected by the surface bulk velocity and by a term of the order of the 30% the ground range component of the wind velocity *W*SIM gr [14]. The latter term represents a bias to be removed for the effective estimate of the surface ground range velocity. Such bias is typically predicted and removed using an ad hoc geophysical model function (GMF), for example, CDOP [12], which exploits the correlation between DCA and the ground range component of the wind velocity *W*SIM gr . The key geophysical process underlying this function is represented by wind waves, due to their orbital velocities.

Following this approach, the existence of a correlation between *v*gr and *W*SIM gr is investigated. Scatter plots of *v*gr vs *W*SIM gr are shown in the first row of Figure 7 (panels a, b from the low and high wind dates respectively). The colour scale reflects the probability density estimation from data computed via a Kernel Density Estimation (KDE) algorithm. The two variables appear not to be correlated in all cases examined here. This was reasonably expected, as in small to medium-size lakes the wind–waves field is typically underdeveloped, with complex dynamics often subject to local effects (such as topography) [57]. Thus, the validated approaches operationally used in the open ocean [10,15] cannot be applied to our case study and other quantities possibly affect *v*gr.

**Figure 7.** Scatter plots of *v*gr vs *W*SIM gr (**a**,**b**); *v*gr vs depth (**c**,**d**); *σ* <sup>0</sup> vs *λ*peak (**e**,**f**); *v*gr vs *σ* 0 (**g**,**h**). Panels on left column (**a**,**c**,**e**,**g**) refer to low wind cases, right column (**b**,**d**,**f**,**h**) to high wind cases. Color scale reflects the probability density estimate (PDE) computed via the Kernel smoothing function, which evaluates the spatial density of nearby points.

Patterns observed in Figures 5 and 6 for *σ* <sup>0</sup> and *v*gr suggest a possible contribution from bathymetry to the SAR products. Figure 7c,d shows that a correlation between *v*gr and bathymetry exists in all the investigated dates: higher correlation between *v*gr and water depth is found the shallower areas (0 to 100 m deep), and in the deepest area (below 300 m). It is recognized that underwater bottom topography can be sensed indirectly via surface effects [62], although the electromagnetic waves emitted by the radar penetrate water only to a depth that is small in comparison to the radar wavelength. In fact, the underwater bottom topography generates short-scale surface roughness variations, which in turn modulate the radar reflectivity. As a complement to the analysis of the correlation between *v*gr and water depth, we note that the deepest region of the lake (depth > 300 m, where highly correlated values are visible in Figure 7c,d) coincide with the northern narrow area of the lake. In this area, the width of the lake ranges between 3 and 5 km, that is, between one and two times the DRC size (2.5 km). Hence, we assume that an effect of the mixed signals from water and land shall also be considered. To what extent such an effect combines with the bathymetry and the wind waves field and alters the measured *v*gr is beyond the aims of this work and we recommend further investigation in this sense.

Given the specific geometry of Lake Garda, the observed significance of bathymetry could, in principle, also be ascribed to the presence of longer (meters scale) wind waves. As shown in Figure 7c,d there is a peak of correlated values around ≈170 m depth and *v*gr ≈ 1.5 m/s. Such a feature, in the case of high wind, is likely related to the wave field development. In this regard, in Figure 8a–c we report the peak wavelength map in the high wind dates. As it is clear from the first two days (panels a,b), *λ*peak reaches its maximum at the end of the northern trunk, where water depth is below the two 150 m isobaths. This area also corresponds to maximum effective fetch for a wind blowing from North-East. The relative weight of the short and long waves also affects *v*gr and is at the base of the applications of the DCA method to oceans [11,61]. Nevertheless, it cannot be easily quantified in lake applications, despite the fact that it can be inferred from *σ* 0 variations. According to the Bragg scattering model, the radar backscatter is proportional to the surface wave spectral density at the Bragg wavelength. As these centimeters-scale waves are tilted by longer (meters-scale) waves, *σ* <sup>0</sup> varies along with such longer wave profiles. This leads to a correlation between *σ* <sup>0</sup> and the orbital velocities of the peak waves [63]. Theoretical findings predict an average increase in the mean spectral density of the short waves due to the interaction between short "Bragg" waves and meters-scale waves [64]. Such a relationship also holds in the high wind cases analyzed here, as shown in Figure 7f, where a clear correlation exists between *σ* <sup>0</sup> and *λ*peak. However, no correlation is found in the case of low wind dates (Figure 7e), when the wave field is not developed.

Figure 7g,h clearly shows that a tight relationship exists between *v*gr and *σ* 0 . A remark about the compensation of the effects of the azimuth variation of the backscattering on the DCA is now in order. The Doppler shift is estimated by exploiting the spectral tapering of the azimuth antenna beam pattern. The algorithms for the DC estimation assume a backscattering spectral flatness, that is, a flat backscattering spectral power density [43,44]. Evidently, the presence of azimuth variations of the backscattering within the Doppler Resolution Cell (DRC) can negatively impact this assumption leading to a spatially varying Doppler shift bias, that is, a disturbance in the measured DCA. In [12], this phenomenon was physically explained by referring to the unbalancing, appearing in the presence of a backscattering azimuth gradient, between the contributions to the Doppler Centroid associated with the scatterers ahead and behind the zero Doppler scatter. In [12], the authors proposed a solution to this problem, which is based on the subtraction of the contribution correlated to the azimuth backscattering gradient from the estimated DCA. A detailed analysis of the effect of azimuth backscattering variations on the estimated DCA, as well as its mitigation, is beyond the scope of the present work and is certainly worth future research, especially for the specific application to lakes. The compensation procedure proposed in [12] was not adopted here due to the possible consequences on the cancellation of the useful signal.

**Figure 8.** Maps of peak wavelength for high wind dates 7 February 2008 (**a**), 13 November 2008 (**b**), 11 June 2009 (**c**).

#### **5. Conclusions**

In recent years, the literature has reflected increasing awareness of the opportunities of satellite-based monitoring of lake surfaces. In this study, we investigate the possibility of inferring lake surface water velocity from existing methods, exploiting the SAR Doppler Centroid Anomaly. To our knowledge, this is the first application of such methodologies to inland waters. The analysis of ENVISAT SAR images and the related DCA maps for Lake Garda returns evident structures with intense and spatially varying values of backscattering and Doppler frequency. These are compared with the results of a numerical model. In the three dates when a high wind blows uniformly over the lake, the spatial patterns retrieved from SAR resemble those from the numerical model, while this is not the case for the low wind dates. In all cases, the surface velocity retrieved from DCA appears to always be significantly overestimated.

The factors potentially influencing the extraction of surface velocity on the investigated dates are identified as are the intensity of the wind forcing, the bathymetry, and the wind waves' wavelength. The simulated wind waves are also found to be well correlated with the backscatter amplitude. Some hypotheses on the potential biases to be removed are addressed, to preliminarily evaluate the applicability of existing GMFs.

Lake dynamics can be significantly different to those in open seas, especially in medium-size lakes where the response to the wind forcing varies locally. In our case study, wind space and time variability is mainly responsible for the surface flow field [59]. The cases analyzed so far suggest that the GMFs derived for inferring the surface velocity from DCA in open seas [12] cannot be applied to our case study, as they account for the wind vector solely and assume a fully developed wind waves field. Indeed, none of the cases analyzed show a clear correlation between ground range velocity measured from SAR and the correspondent wind component. On the other hand, the data processed suggest that a specific GMF could be defined for our case study, and potentially for the more general case of medium-size lakes. This function should include parameterizations handling the small scale variations associated with wind forcing, bathymetry and waves. The possibility of including the information carried by the SAR backscatter should be evaluated, and the interpretation of the environmental variables signature on the SAR backscatter should be further improved.

The application to lakes of Doppler methods commonly exploited in the open sea is not only challenging for the site related features, for example, specific environmental conditions, hydrodynamic and topographic aspects, but also from the data processing perspective. Further developments shall be carried out in the processing chain to account for possible size limitations of the lake compared to the Doppler resolution cell. Adaptive strategies based on multi-resolution Doppler estimation could be investigated, as well as approaches for the correction of the bias induced by a non-uniform backscattering within the Doppler resolution cell. Such data processing developments are considered priorities for improving the application of SAR-based current estimation methods to lakes.

The results and the interpretations provided in this work must be deeply validated and massive research activity is required, as well as dedicated field campaigns. In this sense, we hope that our contribution will pave the way for the development of a procedure for the extraction of the current signature from SAR that is adequately re-adapted for lakes.

**Author Contributions:** Conceptualization, M.A., V.Z., G.D.C., G.F. and F.D.S.; Data curation, M.A., V.Z. and F.D.S.; Formal analysis, M.A., V.Z. and F.D.S.; Investigation, M.A., V.Z., G.D.C., G.F. and F.D.S.; Methodology, M.A., V.Z., M.T., M.B. and C.G.; Project administration, F.D.S.; Resources, M.A. and V.Z.; Software, M.A. and V.Z.; Supervision, G.F. and F.D.S.; Visualization, M.A.; Writing—original draft, M.A., V.Z. and F.D.S.; Writing—review & editing, G.D.C., G.F., M.T., M.B. and C.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study has been supported the EU Horizon 2020 project Water-ForCE (grant nr. 101004186).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is no applicable to this article.

**Acknowledgments:** We thank Lorenzo Giovanni for the WRF simulations. We are also thankful to Marco Pilotti, Giulia Valerio and Francesco Serafino for insightful discussions on the potential impacts of this preliminary work.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


## *Article* **CDOM Optical Properties and DOC Content in the Largest Mixing Zones of the Siberian Shelf Seas**

**Anastasia N. Drozdova 1,\* , Andrey A. Nedospasov <sup>1</sup> , Nikolay V. Lobus <sup>2</sup> , Svetlana V. Patsaeva <sup>3</sup> and Sergey A. Shchuka <sup>1</sup>**


**Abstract:** Notable changes in the Arctic ecosystem driven by increased atmospheric temperature and ice cover reduction were observed in the last decades. Ongoing environmental shifts affect freshwater discharge to the Arctic Ocean, and alter Arctic land-ocean fluxes. The monitoring of DOC distribution and CDOM optical properties is of great interest both from the point of view of validation of remote sensing models, and for studying organic carbon transformation and dynamics. In this study we report the DOC concentrations and CDOM optical characteristics in the mixing zones of the Ob, Yenisei, Khatanga, Lena, Kolyma, and Indigirka rivers. Water sampling was performed in August–October 2015 and 2017. The DOC was determined by high-temperature combustion, and absorption coefficients and spectroscopic indices were calculated using the seawater absorbance obtained with spectrophotometric measurements. Kara and Laptev mixing zones were characterized by conservative DOC behavior, while the East Siberian sea waters showed nonconservative DOC distribution. Dominant DOM sources are discussed. The absorption coefficient *a*CDOM (350) in the East Siberian Sea was two-fold lower compared to Kara and Laptev seawaters. For the first time we report the DOC content in the Khatanga River of 802.6 µM based on the DOC in the Khatanga estuary.

**Keywords:** CDOM absorbance; spectroscopic indices; DOC; Arctic; shelf seas; estuarial and coastal areas

#### **1. Introduction**

The oceanic dissolved organic matter (DOM) pool is one of Earth's large organic carbon reservoirs [1], and, therefore, represents an important component of marine ecosystems and the carbon cycle [2]. The main sources of DOM in the world ocean are primary production of phytoplankton and ice algae, as well as the secondary production of zooplankton. Additional DOM sources include more refractory, compared to the autochthonous material, terrestrial-derived DOM supplied by river runoff (more than 80% of terrigenous DOM), aeolian dust, and coastal abrasion [3].

The Arctic Ocean represents a unique ecosystem, which, on the one hand, is highly sensitive to climate changes occurring during previous decades, and on another, is an important feedback component of the global climate system [4]. An essential feature of the Arctic region is its exposure to large river discharge. Representing only 1% of the global ocean volume, the Arctic Ocean receives more than 10% of the global freshwater and, therefore, the vast amounts of riverine DOM [5,6]. This is a key factor in regulating biogeochemical cycles in the area. In estuarine and coastal zones, fresh waters and terrestrial material control the distribution of flora and fauna, their productivity, and consumption [7]. Here, the limitation of photosynthetic activity due to humic substances, absorbing sunlight in the blue spectral range, where chlorophyll and photosynthetic carotenoids have

**Citation:** Drozdova, A.N.; Nedospasov, A.A.; Lobus, N.V.; Patsaeva, S.V.; Shchuka, S.A. CDOM Optical Properties and DOC Content in the Largest Mixing Zones of the Siberian Shelf Seas. *Remote Sens.* **2021**, *13*, 1145. https://doi.org/10.3390/ rs13061145

Academic Editor: Giacomo De Carolis

Received: 26 January 2021 Accepted: 14 March 2021 Published: 17 March 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/).

absorption maxima [8], is the most significant. Conversely, humic substances absorb in the UVB (280–320 nm) and/or UVA (320–400 nm) ranges [9], diminishing the negative effects of ultraviolet radiation on plankton populations [10]. High concentrations of terrigenous colored DOM (CDOM) rich in humic substances have a significant impact on penetration of light into the water column [11], water color, and its spectral features (see [12] and references therein). Thus, humic substances providing an essential component of the remotely sensed optical signal, affect estimates of chlorophyll concentration using satellite imagery [13–15] or shipborne lidar measurements [16,17].

Numerous studies have contributed to understanding biogeochemical cycling of organic carbon in the Arctic Ocean, see, for example, the PANGAEA database (https: //www.pangaea.de/, accessed on 15 March 2021), the monographs and reviews [18–21], etc. At the same time, seasonal variability of DOC distribution is not quite clear: significant changes, often comparable with multiyear variations, may occur on time scales of a few weeks [19], while most sampling expeditions are restricted to a few months during summer. Thus, for example, Dai et al. [22] estimated the uncertainty of global river DOC discharge to the coastal seas as 30%. Climate change in the Arctic [23–26] has already resulted in reduced ice cover and increased flows of terrigenous DOM due to permafrost thawing and coastal erosion. These changes require monitoring of the concentration and quality of DOM [27] to better understand the Arctic ecosystem and its response to changing conditions and anthropogenic stress. The development of forecasting models for predicting river export of DOM also involves the acquisition of new field data. Thus, a great need for additional DOC data, which can only be obtained in the field, was emphasized by Harrison et al. [28]. Taking into account the recent progress in development of quantitative CDOM and DOC determination with the use of satellite remote sensing [29,30], new field studies of DOC and optical properties of natural waters might be useful for validation of remote sensing models, as well as regional algorithms for estimation of absorption coefficient of colored organic matter [31].

In this paper, we present new data on the content of DOC in the Kara, Laptev, and East Siberian seas obtained during expeditions in summer and fall of 2015 and 2017. We analyze the data on the optical characteristics of the colored DOM fraction [32] to reveal its quality and the sources of DOM input. For the first time, we present DOM concentration at the section from the mixing zone of the Khatanga River to the continental slope obtained with the high temperature combustion technique, which allowed us to estimate the content of DOC in the Khatanga River. In previous studies (see, for example, the paper of Wheeler at al. [33]), the content of DOC in the Khatanga River was evaluated on the basis of total organic carbon concentration (TOC) and the ratio DOC/TOC = 0.9, typical for the Russian rivers [34].

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

#### *2.1. Sampling*

Samples were collected during the 63rd (28 August–6 October 2015) and 69th (22 August–26 September 2017) cruises aboard the R/V *Akademik Mstislav Keldysh* (Figure 1). Exact sampling dates for each sample are given in Supplementary Table S1. In 2015, the sampling was performed along a cross-slope transect in the Laptev Sea, starting from the Lena River Delta region (station 5216) and moving along the 130◦ E towards the continental slope. Absorbance was measured spectrophotometrically for five surface water samples, and DOC concentration was determined in 36 water samples from different depths. In 2017, 141 samples for DOC and 135 ones for absorbance were taken in the Kara, Laptev and East Siberian seas. We consider four shelf-crossing transects from the mixing zones of the Lena, Khatanga, Indigirka, and Kolyma rivers, as well as individual sampling sites in the Kara and Laptev seas. In the Kara Sea, the samples were collected in the estuarine zone of the Ob and Yenisei rivers (station 5588), in the central part of the Kara Sea to the north of Novaya Zemlya Trough (station 5587), and in the Blagopoluchiya Bay (the entrance—5641\_2 and inner part of the bay water area—5642). One more station (5586) was

located close to Novaya Zemlya Trough. In previous years, this region was characterized as less affected by Ob and Yenisei riverine waters [35].

**Figure 1.** Sample site locations during the 63rd (black crosses) and 69th (blue dots) cruises of the R/V *Akademik Mstislav Keldysh*. The data from the study of Amante et al. [36] were used for the map plotting.

Water samples were taken using Niskin bottles of 5 L volume mounted on the CTD/rosette system at surface and discrete depths, associated with boundaries of large gradients of hydro-physical parameters. All the samples were filtered through precombustion at 450 ◦C Whatman GF/F filters with a nominal pore size of 0.7 µm. The filtrate was collected into the acid-cleaned 10 mL glass vials and stored under dark conditions at 4 ◦C until further analysis. For the DOC measurements, the filtrate was acidified up to pH = 2 before storage. Water temperature and practical salinity were derived from CTD measurements.

#### *2.2. Flow Measurements*

During the 69th cruise, we performed continuous observations of conductivity and temperature of subsurface water layer at 2.7 m depth. The measurements were carried out with a temporal resolution of 3 s (at a speed of 10 knots—15 m) using a SBE 21 SeaCAT Thermosalinograph (Sea-Bird Scientific) and a pump to supply outboard water.

#### *2.3. DOC*

DOC concentration was measured onshore by high-temperature combustion with a Shimadzu TOC-VCPH/CPN analyzer. Precision and accuracy of our measurements were determined relative to external laboratory standards, namely solutions of potassium hydrogen phthalate and sodium hydrogen carbonate diluted to different concentrations according to estimated DOC content, and amounted to ±5% and 1%, respectively.

#### *2.4. Optical Measurements and Spectroscopic Indices*

Absorbance A(λ) of water samples have been registered in the laboratory conditions at room temperature 22 ± 2 ◦C. Measurements were performed within the spectral range from 200–700 nm at 1 nm increments using double-beam scanning spectrophotometer Solar PB2201 with 3 or 5 cm quartz cuvettes depending on the DOC concentration and

Milli-Q water as a blank. CDOM absorbance spectra are available online as supplementary information to the study of Drozdova et al. [32]. The blank-corrected absorbance spectra were converted into the Napierian absorption coefficients *a*CDOM(λ) by multiplying CDOM absorbance by 2.303 and dividing by the cuvette path length taken in meters. In accordance with the study of Helms et al. [37], the spectral slope for the 275–295 nm range (S275–295) and the ratio of S275–295 and S350–400 (SR) were determined using linear regression of the log-transformed functions of absorption coefficients *a*CDOM(λ) defined as:

*aCDOM*(*λ*) = *aCDOM*(*λ*0)*e* −*S*(*λ*−*λ*0) , (1)

where λ<sup>0</sup> is a reference wavelength [38]. S<sup>R</sup> index was reported by Helms et al. [37] to correlate with molecular weight and the degree of photochemical degradation of CDOM. Thus, S<sup>R</sup> <1 is typical for terrestrial CDOM, while S<sup>R</sup> > 1.5 indicate the presence of oceanic and photodegraded terrestrial CDOM, see also [12,39]. The value of S275–295 is also used for the CDOM source discrimination so that S275–295 > 20 is typical for marine waters, while S275–295 < 16 is a characteristic of riverine waters. A specific UV absorbance (SUVA) was calculated by normalizing the decadic absorption at 254 nm to the DOC concentration in milligrams per liter (mg/L). It was shown to be a useful parameter for estimating the dissolved aromatic carbon content in aquatic systems [40]. The value of SUVA < 1.8 is an evidence of algae and bacteria CDOM predominance, and SUVA >3 suggests the terrestrial CDOM origin. Spectroscopic indices are given in our recent studies [32,41].

#### **3. Results**

#### *3.1. Subsurface Water Temperature and Salinity*

At the end of August 2017, an influx of warm and salty waters that originated from the Barents Sea, was recorded in the Kara Sea (Figure 2). The eastern type of freshened water distribution [42] was observed in the Kara Sea at the beginning of the cruise. It is characterized by the transfer of low-salinity waters, formed under the influence of Ob and Yenisei river runoff, along the coast towards the Wilkitsky Strait. To the east of 94◦ E and up to the Laptev Sea, the drifting ice was constantly met, therefore, the thermosalinograph was switched off. On the way back, ice did not occur, but in the northern part (north of 76◦ N), the subsurface temperature was still below zero. Freshened waters were shifted westward a month later, which is typical for the central type of freshened water distribution.

In the northwestern part of the Laptev Sea, to the east of the Vilkitsky Strait, surface water temperature was below zero (Figure 2). The highest values of about 5 ◦C were observed for the freshened waters formed under the influence of the Lena River discharge. On the way back, the temperature and salinity distribution in the central and eastern parts of the Laptev Sea did not change significantly. In the northeast, the temperature dropped by 1–1.5 ◦C and negative water temperature values were observed between the station 5634 and the Vilkitsky Strait.

According to the thermosalinograph data, the water temperature of the subsurface layer in the East Siberian Sea decreased from south to north (Figure 2). Negative values were observed near the ice massif, the edge of which was located above the 70 m isobath. Water freshening was observed in the shelf region. At the transect from the Kolyma River estuary to the continental slope, a pronounced frontal zone to the south of station 5615 was observed. Salinity did not change significantly to the north, while to the south it dropped sharply from 28 to 17 within 125 km distance. On the contrary, no large horizontal gradients were observed in the Indigirka section. Similar to the Laptev Sea, the distribution of river waters reaches an isobath of 25 m.

**Figure 2.** Water temperature and salinity at the depth of 2.7 m (thermosalinograph data) along the route of the 69th cruise of the R/V *Akademik Mstislav Keldysh* in the Kara, Laptev and East-Siberian seas.

#### *3.2. Vertical Sensing of Hydrophysical Parameters*

#### 3.2.1. Kara Sea

In the first part of the cruise, the surface waters of the Kara Sea warmed up to 2.5– 5.5 ◦C (Figure 3). With the exception of the stations located in the Blagopoluchiya Bay (5642, 5641, and 5644), a pronounced cold intermediate layer (CIL), formed during winter convection, was observed.

**Figure 3.** Temperature—salinity diagram for CTD profiles observed in the Kara Sea in Autumn, 2017. Color coding indicates DOC concentration (µM).

The core of CIL was located at depths of about 50 m. This agrees well with the data reported by Pavlov and Pfirman [42], who states that summer warming extends to depths of 20–60 m. The upper heated layer thickness in different areas varied from 5 m at station 5641 to 25 m at the Ob-Yenisei coast (5588), and was absent completely in the inner part of Blagopoluchiya Bay (5642). It was reported previously that the depth of the surface layer in the Kara Sea was 6–8 m in shallow parts of the sea, ranging up to 20–30 m over deeper regions [42]. Under the CIL core, the temperature profiles were quite different. The increase of temperature with depth in the south (5586) is much sharper than the one in the northern part of the Novaya Zemlya basin (5587, 5649). In the Novaya Zemlya Trough at depths above 120 m, we suggest the active admixing of waters, which flow from the Barents Sea along the St. Anna Trough. A month later (stations 5586\_2, 5587\_2, and 5588\_2), the surface water temperature dropped by 1–2 ◦C. The remaining profiles have changed insignificantly (Figure 3). Salinity in the studied areas of the Kara Sea varied from 20, measured for the Ob-Yenisei coast surface waters, to 34.6 (deep waters). Below 50 m, all profiles are similar, except for the station 5587, which is characterized by fresher deep waters. High salinity >31 was measured in the Bay and at the southern station 5586, where the waters of the Ob and Yenisei did not extend. The influence of Ob and/or Yenisei riverine waters was reported in 2015 in the Oga and Tsivol'ki Bays of the Severny Island of Novaya Zemlya [35,43], when the western distribution of low-salinity waters [44] was observed. The temperature in the

CIL core at station 5587 was −1.8 ◦C. This roughly corresponds to the freezing point of water with a salinity of 32.8. Taking into account that the formation of CIL is related to vertical mixing during the winter period, we suggest that the salinity in the upper layer of the Kara Sea in winter is quite high—above 32. The waters with lower salinities distributed in the Kara Sea in the summer and fall periods is the result of current-year freshening caused by the Ob and Yenisei rivers, the Baydaratskaya Bay, and melt waters. In the study of Olsson and Anderson [45], such a threshold salinity value in case of the Siberian shelf seas was reported to be 24.

#### 3.2.2. Laptev Sea

Along the transects in the Lena Delta region, water was warmed up to 10–15 m (Figure 4). Closer to the coastline, the temperature of the low-salinity water layer was 5 ◦C and decreased to 1–2 ◦C away from the Lena Delta region. The bottom water layer was characterized by temperatures below zero. The mixing of riverine and sea waters occurred throughout the entire water column of 20 m depth in the southern part of the sections. Compared to 2015, the spreading of freshened waters along the transect was weaker. In 2015, the salinity value of 25 was observed 200 km farther to the north. The waters in the area of the shallow shelf were freshened up to the 20 m depth. A frontal zone was observed above the sloping shelf brow at about 20 m depth, were the waters with salinity of 15 or less were distributed. In opposite, salinity of 15 was measured 50 km southerly in 2015 (Figure 4). We assume that a more active horizontal mixing of sea and river waters took place in 2015. Consequently, low-salinity waters spread far to the north, but at the same time low salinity values were observed closer to the Lena Delta. The halocline was located at depths of 10–15 m along the entire transect, and isohaline concentration within the frontal zone increased above the brow of the slope. This result does not completely support a schematic diagram of the interactions occurring in the Laptev Sea continental shelf close to the Lena River delta region reported by Gonçalves-Araujo et al. [46], since it implied reducing the thickness of plume-influenced freshened water layer, but agrees well with the data on salinity and temperature vertical distribution along the 130◦ E transect discussed by Bauch et al. [47]. A small lens of fresher and warmer waters was observed at station 5595. A similar lens was noted in 2015, but it was larger and located further from the coast. It very likely was formed due to the "offshore" atmospheric forcing [47].

The vertical salinity structure and distribution of freshwater fraction at the transect from the Khatanga River estuary to the continental slope are discussed in detail by Osadchiev et al. [48]. Shortly, the Khatanga plume was weakly-stratified and occupied the whole water column in the shallow inner part of the estuary (stations 5627–2629) due to intense tidal mixing in the Khatanga Gulf. Tidal-induced dilution caused an increase of surface salinity and depth of the plume from 4–7 m (station 5628) to 17–25 m (station 5630) at 120 km along the transect. In the outer part of the estuary, the plume detached from sea bottom and its depth steadily decreased to 11 m, while surface salinity increased to 21 (station 5632). Penetration of marine waters (with salinity above 30) into the bay was observed up to Bolshoy Begichev Island. Closer to the mouth in the bottom layer, salinity gradually decreased from 25 to 7 (Figure 5). The impact of the Khatanga River fresh waters decreased sharply not far from the entrance to the open sea (station 5633) due to primarily moderate Khatanga River flow, which is about five times smaller than that of the Lena River [49]. The boundary of freshened water lies above the 25 m isobath. The surface water temperature varied from 3 ◦C in the Khatanga River estuary to −1 ◦C near the continental slope (station 5635). Deep waters (>30 m) were characterized by temperatures below zero (Figure 5).

**Figure 4.** Temperature and salinity across the mixing zones of the Lena River in fall 2015 (**upper**) and 2017 (**lower**).

**Figure 5.** Distribution of temperature, salinity, and DOC along the Khatanga, Indigirka, and Kolyma transects.

#### 3.2.3. East Siberian Sea

In the East Siberian Sea, surface water temperature decreased from south to north (Figure 5). At the section from the Indigirka River to the continental slope, the surface water temperature near the river mouth reached 6.2 ◦C and gradually decreased to −1.4 ◦C at station 5607 located near the ice edge. Negative temperatures of bottom water were observed on the inner shelf (station 5601) at the depth of 26 m. Low-salinity waters supplied by the river runoff extend to a depth of 10–15 m. At the transect from Kolyma River estuarine region to the continental slope, the surface water temperature on the inner shelf was 6.7 ◦C (stations 5619 and 5620) and decreased to 0.5 at the northern station 5612. In most of the Kolyma section, the total temperature difference in the water column did not exceed 1 degree. Salinity of surface waters increased from southwest to northeast (Figure 5). In the shelf area adjacent to the estuaries of the Kolyma and Indigirka, the minimum salinity was 17 and 15, respectively. The maximum salinity values were recorded at the northernmost stations 5607 (30 and 32.5 at the surface and bottom, respectively) and 5612 (29.2 and 31.2). The vertical distribution of salinity at the Indigirka section showed a pronounced freshened surface water layer, typical for the Arctic shelf under the influence of continental runoff. A unique feature of the eastern part of the East Siberian Sea (Kolyma transect) was a region of nearly 150 km of practically homogenous vertically mixed water column. This area may provide vertical convection down to the bottom during autumn water cooling.

#### *3.3. DOC*

The concentration of DOC in seawater of Arctic seas varied in a wide range between 82.8 and 886.7 µM. The complete DOC dataset is given in Supplementary Table S2. In the Kara and Laptev seas, the higher DOC concentrations were measured for the upper fresher water layer in the mixing zones formed under the influence of Ob, Yenisei, Lena, and Khatanga runoff (Figures 3 and 5, Tables 1 and 2). In the Kara Sea, the mean DOC content in the upper 25-m water layer was approximately two-fold higher compared to deep waters below 25 m. In the Blagopoluchiya Bay, the DOC was lower, compared with the data obtained recently for the Oga and Tsivolki Bays in the case of the western distribution of low-salinity riverine waters [35]. It confirms that terrestrial-derived DOM, supplied by Ob and Yenisei rivers, represents the main DOM source in the bays of Severny Island of Novaya Zemlya archipelago. At the transect from the Kolyma River mixing zone to the East Siberian Sea continental slope, DOC varied between 125.8 and 505.0 µM for the salinity rage 17.0–31.5. The Indigirka transect covered a larger salinity gradient from 15.2 to 33.4. The values of DOC varied there from 165.0 to 526.7 µM. The surface waters of the transect were characterized by moderate DOC concentrations of 236.7–393.3 µM with its local increase at the station 5606 up to 520.0 µM. In contrast to the Kara and Laptev seas, the linkage between DOC and hydrological parameters in the East Siberian Sea was not observed. DOC was distributed rather randomly, showing that significant DOC concentrations of 300–526.7 µM were typical for both Indigirka and Kolyma mixing zones and continental slope region. The DOC values measured in 2017 in the East Siberian Sea were higher than the ones reported by Alling et al. [50]. Thus, in 2008, the upper water layer of 15 m to the west of 160◦ E was characterized by mean DOC of 170 µM (mean salinity S = 22), and to the east of 160◦ E 93 µM (S = 28). In the present study, corresponding mean DOC concentrations were 285.1 µM (S = 26.0) and 168.1 µM (S = 29.7).


**Table 1.** Variation of DOC and optical characteristics of the Kara Sea waters.

**Table 2.** Variation of DOC and optical characteristics of the Laptev Sea and East Siberian Sea shelf waters (2017).


#### *3.4. Optical indices*

3.4.1. Kara Sea

CDOM concentrations, depicted as *a*CDOM(375) [51], followed similar trends to that observed for DOC. The surface of mostly freshened waters had the highest *a*CDOM(375) (up to 3.64 and the mean value 1.43 m−<sup>1</sup> ), while *a*CDOM(375) of Blagopoluchiya Bay and Kara Sea deep waters did not exceed 0.64 m−<sup>1</sup> (Table 1).

Spectral slope ratio S<sup>R</sup> varies between 0.9–3.4. It strongly correlates with *a*CDOM(375) showing exponential decrease to the values, typical for the Ob (~0.87) and Yenisei (~0.91) freshwaters reported by Stedmon et al. [26], see also Discussion Section. Similar dependence was obtained for the S275–295 spectral slope, indicating the predominance of terrestrial material in a fewer number of samples (S275–295 < 20) and mixed or mostly autochthonous DOM character for the others. Interestingly, the data from Blagopoluchiya Bay are grouped separately in the S275–295—*a*CDOM(375) plot and differ by the lower S275–295 values (Figure 6). It apparently reflects local input of terrestrial-derived material from Novaya Zemlya island.

**Figure 6.** Spectral slope S275–295 plotted against *a*CDOM(375) for the Kara Sea waters.

SUVA varied between 0.39 and 2.30 m<sup>2</sup> g C−<sup>1</sup> with the highest values observed for the upper water layer. Since SUVA has been shown to be positively correlated to molecular weight [52] and also to be one of the most reliable parameters for the DOM source discrimination with regard to its changes during bio- and photodegradation [53], we suggest the upper water layer DOM to have higher molecular weight, which is explained by a larger fraction of humic acids [54] supplied by Ob and Yenisei rivers. SUVA values of the Kara Sea waters were found to be lower than the ones of the Ob (2.58–2.75 m<sup>2</sup> g C−<sup>1</sup> ) and Yenisei (1.95–2.97 m<sup>2</sup> g C−<sup>1</sup> ) end members [26].

#### 3.4.2. Laptev Sea

CDOM optical properties suggest the dominance of terrigenous humic substances in surface and thermocline waters (~10 m layer) at the transect along the 130◦ E (see Table 2). Absorption coefficient *a*CDOM(375) varied between 0.34 and 7.20 m−<sup>1</sup> . The lowest absorption coefficient *a*CDOM(350) was measured for the bottom waters from the station 5592 and amounted to 0.5 m−<sup>1</sup> , which is close to the ones reported for the Polar Waters of the East Greenland Current [55,56]. In the study of Pugach et al. [57] a comparable spatial variability of *a*CDOM(350) was demonstrated. The maximal *a*CDOM(350) of 11.2 measured for the surface waters with salinity 6.6 (station 5596\_2) was slightly lower than the one reported for a mid-flow regime of the Lena River of about 13.1 m−<sup>1</sup> [58]. The values of *a*CDOM(350) were found to be lower compared to the data reported by Gonçalves-Araujo et al. for the Lena Delta region [46] (0.9< *a*CDOM(350) < 15.7 m−<sup>1</sup> ), which is likely explained by a smaller terrestrial CDOM contribution. Spectral slope *S275–295* varied between 16.38– 20.35 µm−<sup>1</sup> . Spectral slope ratios *S<sup>R</sup>* obtained in the present study (0.9< *S<sup>R</sup>* < 1.2) are generally consistent with the results of Pugach et al. [57] and Gonçalves-Araujo et al. [46] (2015) (~0.87–1.00) for the Lena Delta—sea mixing zone. For comparison, *S<sup>R</sup>* values of the Lena River water samples were reported to vary between seasons and estimated by Stedmon et al. [26] as 0.81–0.89. Higher *S<sup>R</sup>* values were typical for deep waters at depths 15–44 m, see Table 2. Values of SUVA of the upper 10 m water layer fall in the range 1.72 < SUVA < 2.52 m<sup>2</sup> g C−<sup>1</sup> , which is comparable with the results 1.33 < SUVA < 4.80 m<sup>2</sup> g C−<sup>1</sup> reported by Gonçalves-Araujo et al. [46] for salinities of 0.90–32.63. For deep waters with salinities 30.1–33.9, SUVA varied from 0.38 to 1.28 m<sup>2</sup> g C−<sup>1</sup> , indicating a lower impact of terrestrial-derived DOM.

The Khatanga discharge experienced intense estuarine tidal mixing and therefore was distributed from surface to the bottom in the inner estuary and over the 20–25 m deep water column in the outer estuary, see Section 3.2.2 and the study of Osadchiev et al. [48]. The difference in DOC concentration and CDOM absorption between the upper and deep waters was, therefore, not as pronounced as for the Kara Sea and the Lena Delta region Tables 1 and 2. The data for salinities above and below 25 are summarized in Supplementary Table S4 for convenience. At the beginning of the transect (stations 5627, 5629 and 5630), the entire water column was characterized by maximal along the transect values of absorption (*a*CDOM(375) was 2.29–7.10 m−<sup>1</sup> ) and specific absorbance SUVA (1.44–2.48 m<sup>2</sup> g C−<sup>1</sup> ), while *S275–295* (15.73–16.40 µm−<sup>1</sup> ) and *S<sup>R</sup>* (0.92–1.09) were low. This indicates the predominance of terrigenous CDOM in this location [12,37]. Further north (stations 5631–5633), the contribution of terrigenous CDOM decreases and becomes significant in the upper 10–20 m water layer only (Figure 7). In contrast, deep waters had lower absorption and SUVA (0.32–0.93, mean 0.75 m<sup>2</sup> g C−<sup>1</sup> ). The increase of spectral slope *S275–295* was observed for the bottom waters. At the northernmost stations of the section, the waters of different optical characteristics were observed at depths 10–40 m. They were characterized as being higher compared to oceanic waters [12] absorption at 375 nm (*a*CDOM(375) was 3.0 m−<sup>1</sup> ) and *S275–295* values typical for estuarine and coastal waters with strong humic character (14.1–18.5 µm−<sup>1</sup> ). At the same time, high salinities and spectral slope ratio *S<sup>R</sup>* varying from 2.23–2.39 clearly indicates the presence of oceanic and/or photodegraded terrestrial CDOM. We, therefore, assume that increased absorption at 375 nm is related to the recently produced CDOM [59,60].

**Figure 7.** Distribution of *a*CDOM(375), S275–295, S<sup>R</sup> and SUVA along the Khatanga, Indigirka, and Kolyma transects.

#### 3.4.3. East Siberian Sea

At the beginning of the transect from the Indigirka River mixing zone to the continental slope, the values of *a*CDOM(375) were typical for estuaries and coastal waters [12] and accounted for 0.7–2.4 m−<sup>1</sup> . The waters farthest from the coast (stations 5605–5607) exhibited lower *a*CDOM(375) values of 0.38–0.56 m−<sup>1</sup> , which is a characteristic of oceanic waters [59]. Spectral slope *S275–295* and *S<sup>R</sup>* were distributed along the transect rather irregularly (Figure 7), resulting in similar mean values for the upper 10 m water layer and deep waters. Spectral slope ratio *S<sup>R</sup>* varied between 0.9 and 1.4, indicating that CDOM had a weak humic character and intermediate or strong autochthonous component. The mean values of the specific UV absorbance were 1.45 m<sup>2</sup> g C−<sup>1</sup> for the upper water layer and 1.0 m<sup>2</sup> g C−<sup>1</sup> for the deep waters and testify the predominance of algae and bacterial CDOM. The surface waters of the inner plume (stations 5598–5602) exhibited SUVA up to 2.6 m<sup>2</sup> g C−<sup>1</sup> caused by the presence of terrestrial-derived DOM supplied by the Indigirka River.

Pronounced plume was not seen from the Kolyma River. The impact of the Kolyma River waters was notable for over 150 km from the beginning of the transect (stations 2619 and 5617). Absorption at 375 nm was about 1.5 times lower compared to the Indigirka transect. The maximal *a*CDOM(375) values were measured for the inner plume (1–2.1 m−<sup>1</sup> ), as well as for the entire water column at the northernmost station 5612. Similar to the Khatanga transect, an increase in absorption at 375 nm, observed at salinities >29, was accompanied by a decrease in *S275–295* and by an increase in *SR*, which indicates the recently produced CDOM. The values of *a*CDOM(350) and *S<sup>R</sup>* obtained in the present study are consistent with the data reported by Pugach et al. [57]. The values of SUVA were below 1.8 m2g C−<sup>1</sup> for all the samples, showing little freshwater input.

#### **4. Discussion**

#### *4.1. Conservative DOC Behavior in the Kara and Laptev Seas*

The DOC versus salinity plot (Figure 8) for the Kara Sea showed that the data fell close to the DOC = 927.3 − 23.0 × *Salinity* regression line, characterized by the coefficient of determination *R* <sup>2</sup> = 0.91. No difference in conservative DOM distribution was found between surface (open circles) and depth profile (filled circles) samples. This agreed well with the observations reported by Stein et al. [7], suggesting similar vertical and horizontal mixing in the Ob and Yenisei estuaries. The conservative DOM behavior in the Kara Sea was demonstrated recently by several studies [19,35,61–63]. We summarized the coefficients *a* and *b* for the regression line DOC = *a* − *b*× *Salinity,* obtained during late summer and fall periods 1997–2017 in Supplementary Table S3.

In September 2015, the DOC at the transect along 130◦ E (Laptev Sea) was distributed conservatively. A negative correlation of DOC with respect to salinity is described by the following equation: DOC = 545.5 − 10.4× *Salinity*, *R* <sup>2</sup> = 0.74. In September 2017, four stations were examined for DOC (stations 5592 and 5596 on the way there and back, 5592\_2 and 5596\_2). Higher DOC values were observed and accounted to 887 µM at salinity 6.6. Assuming conservative DOM behavior, our estimates for the DOC in fresh water are 545.5 µM C and 1015.4 µM C in 2015 and 2017, respectively. These results are in a good agreement with the data reported previously (506–1252 µM C) [26,46,64–68]. The mean DOC concentrations for the upper and deep waters are consistent with the data given by Alling et al. [50].

At the transect from the Khatanga River to continental slope, a linear correlation DOC − salinity was described as DOC = 802.6 − 19.3 × *Salinity*, *R <sup>2</sup>* = 0.97, for salinities varying from 3.5–31.5 (Figure 8). Considering conservative DOC behavior within the Khatanga River mixing zone, our evaluation of DOC in the Khatanga River is 802 µM. It is almost twofold higher compared to the DOC value of 472 µM reported by Wheeler at al. [33]. This assessment was based on the mean TOC in the Khatanga River and typical for Russian rivers mean fraction of DOC/TOC = 0.9 [34]. According to the data on annual variation of DOC in six major Arctic rivers [26], DOC concentration in September can exceed the mean annual value no more than 10% (Ob and Kolyma rivers). For the Mackenzie, Yenisei, Yukon, and Lena rivers, the mean annual DOC was found to be even higher than DOC measured during the period from the end of August to the beginning of October. We, therefore, assume that the mean annual DOC concentration of 472 µM in the Khatanga River waters might be underestimated.

**Figure 8.** Distribution of DOC along the salinity gradient in the Kara, Laptev, and East Siberian seas. The data on surface and subsurface (1 m depth) waters are shown by open circles.

#### *4.2. Nonconservative DOC Behavior in the East Siberian Sea*

In the study of Alling et al. [50], nonconservative DOM behavior was revealed in the East Siberian Sea. It was related to the DOM removal, which explained the net DOC deficit. The field studies conducted in 2017 have also demonstrated nonconservative DOC distribution along the salinity gradient in the East Siberian Sea (Figure 8). While in the Arctic region, the maximum DOC content is usually a characteristic of low-salinity waters of the mixing zones affected by river runoff, in this study the regions of the Arctic shelf remote from the estuaries and deltas showed DOC concentrations that were comparable with the ones observed at lower salinities. In order to identify the samples with a high DOC content, which was not caused by newly released terrestrial-derived material, we plotted DOC against absorption coefficient at 350 nm (Figure 9A). The substantial DOC concentrations of >300 µM and *a*CDOM(350) < 2 m−<sup>1</sup> were found at the farthest from the coast stations of the Indigirka section (5605–5607) as well as in the upper 15 m water layer throughout the entire Kolyma section (stations 5613, 5615 and 5617). The possible mechanisms of formation of high salinity waters exhibiting high DOC concentration is an autochthonous DOM production, which is one of the major DOM sources to the marine environment with limited continental influence [69,70]. This assumption is supported by the increase in the above areas the spectral slope *S275–295* (Figure 7). A similar increase in the DOC content, accompanied by an increase in absorption in the short-wavelength spectral range, was observed in the region of the continental slope in the Khatanga section (station 5634).

**Figure 9.** (**A**)—DOC against absorption coefficient at 350 nm for the studied samples; (**B**)—spectral slope ratio plotted against the absorption coefficient at 375 nm; (**C**)—*a*CDOM(350) against salinity for the individual mixing zones; (**D**)—*a*CDOM(440) plotted against *a*CDOM(350).

#### *4.3. CDOM Sources*

A criterion suggested by Helms et al. [37] for identifying the sources of CDOM shows that the dominant contribution of terrigenous OM (*S<sup>R</sup>* < 1) is a characteristic of the mixing

zones of the Khatanga and Lena rivers (salinities 3.5–21.5) (Figure 9B). *S<sup>R</sup>* values varied between 0.91 and 1. For comparison, *S<sup>R</sup>* values for river water samples were reported to vary between seasons and estimated by Stedmon et al. [26] as 0.82–0.92 (Kolyma), 0.81–0.89 (Lena), 0.83–0.92 (Ob), and 0.79–0.93 (Yenisei). CDOM of the Indigirka transect was of mixed autochthonous-allochthonous character, while the stations east of 160◦ E (Kolyma section) are distinguished by the presence of autochthonous CDOM in seawater. This is consistent with the results, published by Semiletov et al. [71], demonstrating that a significant component of freshwater from Siberian river inflows into the coastal East Siberian Sea, extending to approximately 160◦ E, where the long-term average position of the Pacific frontal zone is located.

The predominance of autochthonous CDOM was also demonstrated in the area of the Novaya Zemlya Trough, in the Blagopoluchiya Bay, and the northern part of the Khatanga transect. Expectedly, the decrease of the influence of the Lena, Khatanga, Indigirka, and Kolyma river runoff farther seaward along the studied transects was accompanied by an increase of salinity, *S275–295* and *SR*, while the water temperature, CDOM absorption at 350 and 375 nm, and SUVA decreased. At the northern stations of the Khatanga and Kolyma sections, however, an increase of *a*CDOM(375) was observed and very likely was associated with the recently produced DOM [59].

#### *4.4. CDOM Absorption at 350 nm and 440 nm*

The absorption coefficient *a*CDOM(350) was repeatedly used earlier as a quantitative measure for CDOM concentrations, see for example [58], due to its ability to estimate lignin concentrations and inputs of terrestrial DOM to the Arctic Ocean [26,72]. While no correlation between DOC and salinity was found in the case of mixing zones of the Kolyma and Indigirka rivers, *a*CDOM(350) plotted against salinity (Figure 9C) showed good correlations described separately for each studied water area in Table 3. We suggest that Indigirka and Kolyma river waters were characterized by similar *a*CDOM(350) values of ~6.2 m−<sup>1</sup> . This is about twofold smaller compared to the Khatanga, Lena, and Ob/Yenisei rivers (~12.3 m−<sup>1</sup> ). The obtained results are consistent with data on the Ob, Lena, and Kolyma rivers during mid-flow [58]. Lower *a*CDOM(350) values in the Kolyma River were explained by lower vascular plant inputs during freshet and its more extensive microbial degradation in the Kolyma watershed.

**Table 3.** Coefficients *a* and *b* of the *a*CDOM(350) = *a* +*b* × *Salinity* regression line and corresponding coefficients of determination obtained for the Kara, Laptev, and East Siberian seas during August– September 2017.


Another important optical characteristic of seawater is a CDOM absorption coefficient at 440 nm *a*CDOM(440). CDOM represents an essential constituent affecting ocean color. Thus, in the Arctic Ocean, the contribution of *a*CDOM(443) to the total non-water absorption can reach ~50% [73]. It was shown that systematic differences in chlorophyll retrievals resulting from different ocean color models are related to each model's ability to account for the absorption of light by CDOM [14]. In the present study, most of the *a*CDOM(440) data showed a good negative correlation with salinity, similar to the ones reported for the *a*CDOM(350) absorption coefficients (Figure 9C). Some of the water samples, however, were failed to be described by the linear dependence on salinity due to high absorption coefficients *a*CDOM(440); they were taken in the Blagopoluchiya Bay, at the station 5586 in the Kara Sea and northern parts of the Khatanga (stations 5633 surface waters, 5590\_2

24 m, and 5634 surface waters and 18 m) and Kolyma (5612) transects. In the *a*CDOM(350) against *a*CDOM(440) plot (Figure 9D), such data points are grouped separately, they are characterized by *a*CDOM(350) below 4 m−<sup>1</sup> and higher *a*CDOM(440) values. This group of points that are aligned on a regression line characterized by different slopes represents a group of CDOM absorption spectra characterized by a shallower slope coefficient. The increase of *a*CDOM(440) may be caused by autochthonous CDOM. Thus, in the study of the South Brazilian Bight, a strong correlation between Chl-a and *a*CDOM(440) was revealed and described with a regression line lying close to the one observed for the global pelagic oceans reported by Bricaud et al. [74]. The CDOM—Chl-a correlation allowed suggesting the presence of an autochthonous source of CDOM to the region driven mostly by the phytoplankton community over the shelf domain [75]. Our assumption is also supported by the study of phytoplankton of the Khatanga transect [76]. It was found that the area of the continental slope in the western Laptev Sea represents a specific local biotope. Phytoplankton in the area of the continental slope was characterized by high abundance and biomass, dominance of diatoms, and the formation of the deep maximum formed by actively growing algae. At the station 5633, maximum of phytoplankton biomass was observed in surface waters, while at station 5635 (about 40 km north from 5634), it was found at 45 m depth.

Obtained results on CDOM absorption can be valuable in remote sensing and modeling issues. For example, calculated from satellite data, CDOM absorption coefficients may be used as an effective indicator of the Kara Sea surface desalinated layer distribution and dynamics [31]. As the values of light absorption in this layer are significantly higher than in surrounding seawaters [56,77], its characteristics must be taken into account in heat budget models.

#### **5. Conclusions**

The complex field studies, conducted in fall 2015 and 2017, covered a large area of the eastern Arctic shelf of the Kara, Laptev and East Siberian seas. Analysis of DOC concentration and CDOM optical properties, supported by CTD data, allowed us to consider DOM distribution and its quality in the mixing zones of the Ob/Yenisei, Khatanga, Lena, Indigirka, and Kolyma rivers. It was demonstrated that the Kara and Laptev mixing zones were characterized by conservative DOC and *a*CDOM(350) behavior, while the East Siberian sea waters showed nonconservative DOC distribution. We provide the first estimates on DOC content, based on the high-temperature combustion technique, in the Khatanga River during mid-flow regime, it accounted for 802.6 µM (9.6 mg/L). Assuming conservative DOM behavior, our estimates for the DOC in fresh water are 545.5 µM C and 1015.4 µM C in 2015 and 2017, respectively, which is consistent with the results of previous studies. Despite the individual watershed characteristics of the rivers flowing into the eastern Arctic shelf seas, variation of the absorption at 350 nm along the salinity gradient was found to be similar for the Laptev and Kara seas. Absorption of the East Siberian Sea waters was found to be two-fold smaller, which is explained by lower CDOM content in the Indigirka and Kolyma rivers, as well as degradation of humic substances supplied by the Lena River during the transport to the East Siberian Sea through the Dmitry Laptev Strait.

Estuarine and delta regions were characterized by the predominance of terrestrialderived DOM supplied by river runoff. The increase of DOC content was observed at the most distant from the shore stations in the area of the continental slope. It was frequently accompanied by growth of absorption at short-wave spectral range (*S275–295*), *SR*, and *a*CDOM(440), which indicates the production of autochthonous DOM by marine biota to be the dominant CDOM source at those locations. The literature overview also demonstrated the correlation between high DOC values and the increase of phyto- or zooplankton populations.

The OLCI ocean color scanners launched in February 2016 (Sentinel-3A) and in April 2018 (Sentinel-3B) should provide satellite data in the next decade. It was shown by Glukhovets et al. [31] that the standard OLCI algorithm for estimating the CDOM absorption coefficient ADG443\_NN gives high errors in the Arctic seas. The dataset presented in this work may be used to improve existing standard and regional [78] algorithms and to create new ones.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2072-429 2/13/6/1145/s1, Table S1: Sampling dates during the 63rd and 69th cruises of R/V *Akademik Mstislav Keldysh*, Table S2: Salinity and DOC concentration of water samples, Table S3: Coefficients *a* and *b* of the DOC = *a* +*b* × *Salinity* regression line and corresponding coefficients of determination obtained for the Kara Sea waters during August–September periods 1997–2017, Table S4: DOC and optical characteristics of the Khatanga transect waters (Laptev Sea).

**Author Contributions:** Conceptualization, A.N.D.; data curation, A.N.D., A.A.N., N.V.L. and S.A.S.; formal analysis, A.N.D.; funding acquisition, A.N.D.; investigation, A.N.D., A.A.N., N.V.L., S.V.P. and S.A.S.; supervision, A.N.D.; visualization, A.A.N. and A.N.D.; writing—original draft preparation, A.N.D. and A.A.N.; writing—review and editing, A.N.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study of the Kara Sea waters was performed in the framework of the state assignment of IO RAS (theme 0128-2021-0016). The study of DOC along the Khatanga and Lena Delta transects was supported by the RFBR (project 18-05-60214). Investigation of CDOM of the East Siberian Sea was supported by the RSF Grant (project 18-77-00053).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data on DOC and salinity of are given in Supplementary Materials to this study. Spectroscopic indices and CDOM absorbance spectra are provided in [32].

**Acknowledgments:** The data are included as supporting information. The authors are grateful to Mikhail V. Flint for the opportunity to take part in the impressive complex expeditions to the Arctic shelf seas and the crew of the R/V *Akademik Mstislav Keldysh* for their contribution during field studies. We thank Dmitry I. Glukhovets, Alexander A. Osadchiev and Marina D. Kravchishina for fruitful discussion, as well as anonymous reviewers for their comments on an earlier version of this paper.

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

#### **References**

