North to South Variations in the Suspended Sediment Transport Budget within Large Siberian River Deltas Revealed by Remote Sensing Data

Abstract

This study presents detailed suspended sediment budget for the four Siberian river deltas, representing contrasting conditions between Northern and Southern environments. Two of the studied rivers empty their water and sediments into the marine located in the permafrost zone in the Arctic region (Lena and Kolyma), and the other two (Selenga and Upper Angara) flow into Lake Baikal located in the steppe and forest-steppe zone of Southern Siberia. For the first time, these poorly monitored areas are analyzed in terms of the long-term and seasonal changes of spatial patterns of suspended sediment concentrations (SSC) over distributaries systems. Remote sensing reflectance is derived from continuous time series of Landsat images and calibrated with the onsite field measurements of SSC. Seasonal variability of suspended sediment changes over deltas was captured for the period from 1989 to 2020. We identify significant variability in the sedimentation processes between different deltas, which is explained by particularities of deltas networks and geomorphology and the existence of specific drivers—continuous permafrost impact in the North and abundant aquatic vegetation and wetland-dominated areas in the South. The study emphasizes that differences exist between Northern and Southern deltas regarding suspended sediments transport conditions. Mostly retention of suspended sediment is observed for Southern deltas due to sediment storage at submerged banks and marshlands located in the backwater zone of the delta during high discharges. In the Northern (arctic) deltas due to permafrost impacts (melting of the permafrost), the absence of sub-aquatic banks and river to ocean interactions of suspended sediment transport is mostly increased downwards, predominantly under higher discharges and along main distributary channels. These results shine light on the geochemical functions of the deltas and patterns of sequestering various metals bound to river sediments.

1. Introduction

Most of the World’s rivers deltas are hot spots in productive land construction [1], where the main river distributes water among comprehensive channel networks and connected wetlands [2]. Deltas with complex channel networks influence terrestrial sediment flux [3,4,5]. The stored material contributes to the sub-aerial growth of a delta, affects the morphology of channel networks, and therefore, influences the overall distribution of floating particles. According to Passalacqua et al. [6], river deltas should have a high sediment supply and a long-term net profit in the land area. The deltas are recognized as an important lateral geochemical barrier, where over 90% of suspended and 40% of dissolved carbon associated with global geochemical fluxes is trapped [7]. At the same time, knowledge about the spatiotemporal distribution of sediment properties within the deltas is limited due to the infrequent single character of sediment monitoring, and does not reflect long-term and seasonal variations. Influence on storage functions and processes of deltas by various hydroclimatic processes and channel morphology remains unclear, e.g., climatic change effects on sediment supply along deltas [8,9,10]. An increased sea level results in the flooding of the terrestrial parts of the deltas [6], which drives sinking of the deltas and impacts sediment provenance in the deltas. It should be noted that in the northern latitudes, in the Arctic Ocean Basin, the impact of climate change is a direct cause of the intensification of the arctic hydrologic cycle. Here, river flow experiences approximately a 7% increase over the period from 1936 to 1999 [11], i.e., increase in spring by 49%, winter by 31% [12], and a further growth is expected [13]. Ahmed et al. [14] recently proved that of all the great Arctic catchment areas especially in the Lena basin during 1980–2000, peak freshet magnitudes were largely increasing. This may have a potential shift from a distinctly nival river regime to a more hybrid or even a pluvial regime [15].

The recent detailed studies on suspended sediment budget over some large deltas such as Volga [16], Terek [17], Atabaska [18], and Niger [19] revealed uncertainties in understanding patterns of sediment concentrations along deltas. Deltas with complex channel networks can efficiently distribute and store sediment loads along the channels and within floodplain water bodies [3,4,5]. Here, submerging of the banks and wetlands during floods was considered at the main driver of sediment retention. At the same time, increase of sediment transport has been reported based on grab sampling campaigns. For example, for the Lena delta, both increase [20] and decrease [21] of sediment load have been reported. Different environments might represent contrasting patterns of suspended sediment concentration (SSC) due to the contradictive nature of flow–channel interactions and channel–floodplain exchange.

Due to mentioned limitation in knowledge, there is an open question “Will deltas be sustainable in the near future?” [6]. Applications of remote sensing data are the only plausible and efficient tool for identifying long-term sediment concentration patterns over extended river reaches [18,22]. Few studies [23,24] successfully demonstrated the capability of various satellite images to capture the spatial and temporal variability of suspended sediments over large rivers as well as smaller streams and plumes over river–ocean interface [25,26]. The present study is based on the application of Landsat images to study spatio-temporal variation of suspended sediment transport over large deltas located in the South (steppe and forest-steppe zone and discontinuous/sporadic permafrost—52.0–52.4 N for the Selenga delta and 55.6–56.0 N for the Upper Angara delta) and North (tundra, continuous permafrost area, 190 km from north to south—72.0–73.8 N for the Lena delta and 68.5–69.4 N for Kolyma). These deltas are considered as typical for the particular natural zones and thus represent possible North to South variations in sedimentation processes. Based on Landsat images collected, obtained, and processed for the last 30 years, we aimed at (1) identifying typical suspended sediment variations along studied deltas; (2) identifying environmental drivers of sediment transport changes within particular deltas and possible changes of the existing patterns; and (3) providing a comparative study of suspended sediment patterns between Northern and Southern Siberian deltas. The standardized methodology and abundant ground-truth observations both with applied remote sensing model verification enable to provide a novel comparative study.

2. Materials and Methods

2.1. The Study Area

The Selenga River delta (Russia) is a large (680 km²) fluvially-dominated fresh water system that transfers water and sediment from an undammed drainage basin into Lake Baikal, a UNESCO (United Nations Educational, Scientific and Cultural Organization) World Heritage Site. This is almost twice bigger then the Upper Angara delta (

Table 1. Main characteristics of the studied deltas.

Delta	Deltas Area, km2	Catchment Size, km2	Average Discharge at the Delta Topset, m3 s−1 (for Period)	Peak Discharge at the Delta Topset, m3 s−1 (Date)	Number of Main Distributaries	Gauging Station (Distance to the Delta Topset, km)
Selenga	680 [30]	447,000 [30]	865 (1934–2015) [22]	6420 (1993) [22]	7	Mostovoy (150 km)
Upper Angara	400 [31]	21,400 [31]	269 (1946–2016) [31]		2	Verkhnyaya Zaimka (10 km)
Kolyma	350 [32]	526,000 [32]	3210 (1977–2004) [32]	42,900 (1999) [32]	2	Kolymskoe (150 km)
Lena	32,000 [33]	2,430,000 [33]	16,800 (1935–2005) [32]	220,000 (1989) [32]	4	Kyusyur (380 km)

), which is the second largest delta among Lake Baikal tributaries. The area is located in southeastern Siberia, Russia (Figure 1) and is formed by a half-graben–styled rift basin, which extends over 700 km in length, has an average width of 60 km, and has a maximum depth in excess of 1600 m [27]. The elevation of both deltas is less than 2 m above the average water level of Lake Baikal. The Selenga delta contains three main distributary systems (sectors): the Lobanovskiy channel system (Right), the Selenga (Kharauz) channel system (Central), and the Levoberezhnaya channel system (Left), named after the largest channels of the dispersal system (Figure 1). Seven main distributaries are averaged between three sectors: P1 (Levoberezhnaya branch), P2 (Main channel), and P3 (Galutaj branch) are located within the left channel system; P4 (Srednyaya) and P5 (Kolpinnaya)—the central channel system; and P6 (Lobanovskiy) and P7 (Dologan))—right channel system. Altogether, the delta comprises more than 30 channels of various sizes. Upper Angara delta consists of two main distributary systems. The left system is fed by few tributaries, whereas the right is dominant in terms of water discharge. Over 49% of the water discharge enters a large pool separated from the recipient (Lake Baikal) by extended bar. A significant part of both deltas is swamped and possesses a great number of tiny lakes. Both deltas experience significant changes in water regime—decline in the case of Selenga river [28] and increase in the case of Upper Angara River [29].

The Lena River delta is the largest Arctic delta with an area of about 32 000 km² [32]. It consists of more than 800 branches with a total length of 6500 km. Here four main distributaries (Figure 1) represent the main sectors of the delta: The biggest is the Trofimovskaya branch (II); the second largest channel system by water discharge is the Bykovskaya branch (I) that turns sharply to the east after Sardakh Island and flows into Buor Khaya Gulf. The third largest Olenekskaya branch (IV) flows westward into the Kuba Gulf, whereas the large sector of small channels is called the Tumatskaya branch (III) [34]. Significantly smaller then Lena delta, Kolyma delta consists of two main distributaries. Both deltas were formed by the accumulation of sediments from the river and from the erosion and abrasion of ancient landforms due to sea-level fluctuations and the tectonic movements of the Earth’s crust. The Lena delta is underlain by continuous permafrost 500–600 m thick [35]. The active (annually thawed) layer is 30–50 cm thick. The entire deltas lie in the tundra zone except for the southern part.

2.2. Methods

The estimation of sediment concentrations is based on the sensitivity of remote sensing reflectance to the floating matter. Most of the existing works consider site-specific empirical relationships between SSC and reflectance at different satellite wavebands by fitting in situ measurements of SSC (mg/L) with satellite-derived reflectance.

Sediment concentration is mainly estimated using the red [36,37,38] (and near-infrared) bands; the green band [39] is used rarely, and panchromatic images are considered only at single cases. In general, the red band reflectances have the most robust correlations with SSC over the range of values. Our algorithm of SSC assessment is based on using visible red Landsat satellite band (630–680 nm for Landsat 8 and 630–690 nm for Landsat 5).

Our approach used in this study relies on regression relationships between SSC and either raw pixel values of satellite images, so called DN—digital numbers, which depend on the remote sensing system, satellite sensor parameters, and approach only one remote sensing system [39,40,41], or reflectance coefficients (absolute units, ρ). The latter is independent from the remote sensing instrument and sensing parameters [36,42,43].

The image processing workflow to estimate suspended sediment budget and extract surface concentration maps consisted of several continuous steps (Figure 2). The calibration of the SSC = f(ρ, DN) relationship was done based on the combination of the SSC field data with the available Landsat image taken during the period of the in situ observations. Next step is model development, which is done as empirical relationships between these two parameters, plotted on a graph. On the total available dataset, we carried out statistical analyses to generate least square regression models to predict surface sediment concentration as a function of surface reflectance for each selected river. The equation of the approximating line presents the relationships of image brightness and SSC. The highest accuracy is achieved using a single technique and preparation of the satellite data on the image brightnesses to eliminate radiometric and atmospheric influence.

All Landsat 5 and Landsat 8 images were pre-processed to a radiometric correction that converted raw digital numbers (DN) values into a reflectance coefficient. The radiometric correction of the Landsat images was derived by [44] in the ArcGIS 10 software. The use of this correction eliminates the systematic differences in the lighting due to the differences in the heights of the Sun above the horizon and in the distances from the Earth to the Sun on different survey dates. It includes a radiometric calibration and correction for lighting conditions.

Radiometric calibration converts the DN values into physical units of the spectral radiance at the sensor’s aperture—$L_{\mathsf{\lambda }}$ (W/m² sr μm) recorded by the sensor on the satellite. Next, the obtained values of $L_{\mathsf{\lambda }}$ are converted to the TOA (top-of-atmosphere) reflectance values at the upper boundary of the atmosphere—${\mathsf{\rho }}_{\mathsf{\lambda }}$.

Three advantages to use reflectance instead of at-sensor spectral radiance were revealed by Chander [44]. Firstly, it removes the cosine effect of different solar zenith angles due to the time difference between data acquisitions. Secondly, TOA reflectance compensates for different values of the exoatmospheric solar irradiance arising from spectral band differences. Thirdly, the TOA reflectance corrects for the variation in the Earth–Sun distance between different data acquisition dates. These variations can be significant geographically and temporally.

The atmospheric correction eliminates relative or absolute influence of haze and clouds that creates variation in the values of pixel brightness of the same objects. The “Dark Object Method” implies the presence of objects on the image with the minimum brightness in a red spectral range. This method is based on selection of one constantly dark object on the reference image and applying the pixel value of this object to calibrate all the images in the series [45]. For relative atmospheric correction, we reduced a series of images “to the common denominator”, which is called “reference” image. The remaining images of the series are reduced to first reference values. On the series of images taken on different dates, dark objects have the same brightness and any differences on the images are the influence of the atmosphere. In case of the research reported here, the “dark object” was the deep part of Baikal Lake in the case of Selenga and Upper Angara images. For the Arctic deltas, clean thermokarst and floodplain lakes were used as a “dark object”. Then each satellite image calculated differences between reflectance coefficients of “dark objects”.

2.3. Materials

For calibration of the SSC model, routine measurements of suspended sediment concentrations were performed. Surface water samples (0.5 m from the surface) were collected with a bottle sampler at various locations of the river channel. Water samples were filtered through pre-weighed membrane filters (with pore size 0.45 µm) from the “Millipore” filtration system. The samples were then oven-dried (105 °C) and re-weighed to determine suspended sediment concentrations (SSC, mg/L = g/m³). Grain size analysis of suspended sediments was conducted with the Laser Particle Sizer Fritsch Analysette.

The field campaigns with SSC measurements were done in the delta of Selenga River and the adjacent water area of Baikal in 6–23 July 2011 (98 measurements)) and compared to Landsat 8 image from 23 July 2013 (Figure 3a). The field measurements at Upper Angara river delta were taken on 3–5 August 2017 and compared with Landsat 8 image from 3 August 2017 (Figure 3a). For the generation of regression models for the Lena River delta, we used field data on the measured SSC of water during expedition in the middle of the Lena River from Yakutsk to the Viluy River in 25 June–9 July 2020 (26 measurements) and 43 images from Landsat 5 and Landsat 8 for the period from 2000 to 2020 (Figure 3b). For the Kolyma River delta, we generated a single model of the SSC of water on the reflection coefficient and used a date of field campaign in lower reaches of the Kolyma River near the Cherskiy (12 measurements) during 3–18 August 2019 (Figure 3c). We used for correction a reference image from Landsat 8 (9 July 2016) for Lena and Kolyma Rivers. All images were compared to this image by the method of relative atmospheric correction. In all cases, time delays between the field sampling and satellite imagery were less than 5 days and were related to homogenous conditions of water discharges, thus minimizing the possible effects of changed sediment flow properties on the model.

A linear regression model built from a combination of data for each river displayed a statistically significant relationship (

Table 2. Obtained regression models for SSC and number and period of Landsat 5 and Landsat 8 images analyzed in the study.

Delta	Equation	Number and Period of Images
Selenga	$SSC=429\times \rho -17.5$ R = 0.87	82 images for the period from 1989 to 2015
Upper Angara	$SSC=1.29\times DN+4.08$ R = 0.84	87 images for the period from 1986 to 2018
Kolyma	$SSC=4.32\times e^{35.4\mathsf{\rho } }$ R = 0.79	92 images for the period from 1999 to 2020
Lena	$SSC=12 345\times {\rho }^{1.8}$ R = 0.84	43 images for the period from 2000 to 2020

$SSC$—suspended sediment concentration of water, $\rho$—the reflectance coefficient, and $DN$—digital numbers of the Landsat 5 and Landsat 8 image in the red spectral region; R—correlation coefficient; $e$—exponential regression model.

, Figure 4). In all cases, the robust empirical models were derived between field surface sediment concentration and surface reflectance data (0.79 to 0.92, slopes significant at 99% confidence level). The verification of the revealed models was performed for each of the studied rivers based on the independent measured SSC datasets. For the Selenga River, the whole measured SSC dataset was divided into two parts: Two of three of the points were used for model calibration, and one out of three for model verification. The accuracy of the modelled SSC calculated at the difference between modelled and measured values is less than 2 mg/L, which is below 20%. Verification of the model for Lena and Kolyma was also based on the independent measured SSC datasets, obtained for Lena Middle reach and Lena delta by Sankt-Petersburg State University team. In the Lena delta, five measurements of surface SSC in the Trofimovskaya branch (II, Figure 1) were done between 19 July to 27 July 2016 [46]. These observations were compared with modeled SSC from Landsat 8 imaged on 17 July 2016. The accuracy of the modelled SSC calculated at the difference between modelled and measured values is less than 8.3 mg/L, which is below 25% error threshold. The difference between field measurements and satellite-derived estimates should be attributed to SSS concentration heterogeneity across the river reach. These results generally agree with the error estimates for the remote sensing SSC models for other rivers of the World. For the Amazon river, J. Martinez et al. [24] reported 36% relative error, whereas E. Park and E.M. Latrubesse [23] highlighted statistically significant relations between SSC and SR at the three gauge stations of Amazon River and its tributaries at the 99% confidence level. In the study of subglacial plume of Svalbard tidewater glaciers [26], the root mean square error (RMSE) between Landsat-8 and in situ measurements of surface reflectance was found to be 3.4%, whereas for the Mekong River [47], RMSE estimated for different seasons of Landsat TM reflectance and in situ measurements also confirmed a high applicability of satellite image for monitoring of suspended sediment loads.

The significance of the obtained relationships in our study is proved by high spatial resolution remote-sensing data (up to 30 m) compared to river size (the width of the river channels is between 200 and 2000 m), and the regionalization (river-specific) of the surface reflectance, which considers the effects of various sediment sizes, concentrations, and material types, including both mineral and biological components [48,49]. It is also explained by relatively low values of sediment concentration at the downstream reaches of the large river systems, which are below 80 mg/L, and which are reported as the threshold level for the high accuracy of determining SSC from surface reflectance (e.g., [50]) above which the surface reflectance saturates. Theoretical analysis [51] indicates that the algorithm’s sensitivity to variable mass-specific particulate absorption has a negligible impact on SSC in the low reflectance regime, but becomes an additional source of error for high reflectance values.

The obtained models (Table 2) were compared to each other and to the ones reported for other river systems of the large rivers of the World [22,23,38,39,52]. Different relationships between SSC and ρ are explained by differences in sediment properties, mainly grain size. All models were obtained for average discharge conditions of the studied rivers, thus the effects of discharge cycle on regression models were negligible.

Eighty-two images of Landsat 5 TM and Landsat 8 OLI (spatial resolution of 30 m in the red and near-infrared (NIR) channels) for the period from 1989 to 2015 for Selenga River delta, and 87 images for the period from 1986 to 2018 for Upper Angara River delta were taken [8]. Empirical models were applied to 50 images of Lena River delta taken between 1999 to 2019 and 70 images of Kolyma River delta.

Representative distributaries were selected for each delta to calculate longitudinal change in SSC (Figure 1). Further, the estimates of suspended sediment budget were averaged for each sector. For Selenga delta were selected three sectors with seven representative distributaries for calculation of longitudinal change in SSC. Distributaries P1 (Levoberezhnaya branch), P2 (Main channel), and P3 (Galutaj branch) are located within the left channel system (first sector); P4 (Srednyaya) and P5 (Kolpinnaya) the central channel system (second sector); and P6 (Lobanovskiy) and P7 (Dologan) the right channel system (third sector). In the Upper Angara delta were selected one channel sector from the top of the delta to the lake through the main channel. For Lena delta were selected five sections: the main channel and four branches: Bykovskaya, Trofimovskaya, Tumatskaya, and Olenekskaya channels. In the Kolyma delta: main channel (Cherskiy) and two branches: Pokhodskaya and Kamennaya channels.

The SSC balance (${\mathsf{\Delta }\mathrm{S}}_0$) was calculated for each distributary on the basis of the difference between SSC at the upstream ($SS{C}_1$) and downstream ($SS{C}_2$) ends [8]:

$${\mathsf{\Delta }\mathrm{S}}_0=SS{C}_2-SS{C}_1,$$

(1)

Values for $SS{C}_1$ and $SS{C}_2$ were estimated as an average for every 500 m of each distributary. Relative sediment retention $\mathsf{\Delta }\mathrm{S}$ (%) was calculated as:

$$\mathsf{\Delta }\mathrm{S} = {\mathsf{\Delta }\mathrm{S}}_0/SS{C}_1,$$

(2)

Three possible ratios between sediment input and output were analyzed: $SS{C}_2$>$SS{C}_1$, $SS{C}_2$ < $SS{C}_1$, and $SS{C}_2$≈$SS{C}_1$, which correspond to either erosional ($\mathsf{\Delta }\mathrm{S}$ > 0) or deposition ($\mathsf{\Delta }\mathrm{S}$ < 0) patterns.

Relative sediment retention ($\mathsf{\Delta }\mathrm{S}$, mg/L) for each sector were compared with the discharges in the upper part of the delta: for Selenga—Mostovoy gauging station; for Upper Angara—Verkhnyaya Zaimka gauging station; for Lena—Kyusyur hydrological station; and for Kolyma—Srednekolymsk hydrological station. We have selected the best available hydrological stations where monitoring is carried out by the Russian Hydrometeo Service.

The comparisons of deltas are based on the conversion average daily water discharge to relative values ($Kq,i$) by average long-term water discharge (for each delta separately)

$$Kq,i=\raisebox{1ex}{$Q_o$}\!\left/ \!\raisebox{-1ex}{$Qav$}\right.$$

(3)

where $Q_o$ is observed water discharge, and $Qav$ islong-term average water discharge.

$$Qav=\raisebox{1ex}{${\sum }_{i=1}^n{Q}_o$}\!\left/ \!\raisebox{-1ex}{$n$}\right.$$

(4)

where $Qav$ is therelative water discharge.

The delta topset discharge data used in Equation (3) relied on Roshydromet gauging stations listed in Table 1 and were downloaded from Arctic Great River Observatory archive [53].

Image processing in ArcGIS and QGIS software resulted in the format of ESRI Shapefile format layers, which contain information about the SSC value with an accuracy of 30 m. SSC data interpolation and construction of the SSC model were performed using the Multilevel B-spline algorithm for spatial interpolation of scattered data as proposed by Lee, Wolberg, and Shin [54]. The algorithm makes use of a coarse-to-fine hierarchy of control lattices to generate a sequence of bicubic B-spline functions.

3. Results

The Equations (1) and (2) were used to estimate the suspended sediment budget along river deltas. The calculations for the four deltas show that both erosion ($\mathsf{\Delta }\mathrm{S}$ > 0) and deposition ($\mathsf{\Delta }\mathrm{S}$ < 0) patterns can be observed under changing hydrological conditions at various parts of particular deltas. At the same time, in some systems, the specific role of sedimentation or erosion processes was found.

The most evident deposition patterns were observed at Selenga and Upper Angara river deltas. In the Upper Angara delta, 49% of the sediments, which enter separately from the recipient estuary (Figure 1), are stored. The remaining 51% of the flow is delivered to the Lake by two main distributaries that have a clear negative downward trend of SSC along the river channel from the delta top to the lake edge. Along the main channel, on average, only 18% of sediments are delivered to Lake Baikal, whereas the remaining 82% are stored in the channel and near-channel wetlands. The abundant aquatic vegetation increases the hydraulic flow resistance here, which increases the depositional role of this section. Downstream of the channel, dense macrophytes form permanent sub-aqueous channels fully determine the planform of the channel (Figure 5b).

In the Selenga River, $\mathsf{\Delta }\mathrm{S}$ (2) is controlled by water discharge. Longitudinal increases in SSC were observed only during low-water periods (discharges less than 1500 m³ s⁻¹). Relative values of sediment retention ($\mathsf{\Delta }\mathrm{S}$, %) generally similar between various branches were higher in the left or right delta sectors than in the central delta sector. Average sediment retention over 82 observational patterns (Landsat images) was estimated to be 6% in the branches of the left sector (I), 3% in the central (II), and 6% in the right (III) sectors. The highest deposition rate was observed in the Selenga delta sector ($\mathsf{\Delta }\mathrm{S}$= 58%, when Q_o = 575 m³ s⁻¹). In the other delta sectors, the highest deposition rate was 17% (II) and 39% (III). The same data ordered by months of the captured conditions in the delta show that in May and June (the beginning of the growing season in the delta [55]), the sediment concentrations mostly increase, which implies net erosion along the channels. The most intensive sediment storage occurs in wetland-dominated areas of the lower part of the delta [56], in particular small ponds and backswamps, connected with the main channel over the year or in particular floods (Figure 5a). The average $\mathsf{\Delta }\mathrm{S}$ for all seven delta distributaries can be approximated as a function of Selenga water discharge:

$$\mathsf{\Delta }\mathrm{S}={10}^{-6}{Q}^2-0.091Q+4.654$$

(5)

The analysis of satellite images on the Lena delta section revealed contrasting patterns $\mathsf{\Delta }\mathrm{S}$ in the left and right sectors of the delta: Erosion is dominant in the Bykovskaya and Trofimovskaya channels. The predominance of a positive balance ($\mathsf{\Delta }\mathrm{S}$ > 0) is associated with thermal erosion processes. The satellite images reveal a significant increase of sediment concentrations along the left bank, which suggests the role of average daily air temperatures in bank erosion. Based on a comparison with air temperature datasets (according to the Tiksi weather station), the increase in the average daily air temperatures influences suspended sediment budget along Bykovskaya and Trofimovskaya distributaries [57]. This is seen especially in the range of air temperatures from 5 °C to 15 °C. The influence of a positive influx of total solar radiation and an increase in air temperatures contributes to the degradation of permafrost rocks and the activation of the thermal erosion [58]. Summer thawing and associated activation of bank erosion raise sediment delivery into the hydrological network and further increase the concentration of solids in the river (Figure 5c). Such effects are mostly marked along the left river banks of the southern exposure of sublatitudinal sections of the Bykovskaya and Trofimovskaya channels. The predominance of the erosion patterns is typical for the largest channels (more than 75% of the total number of situations).

The accumulative regime mostly dominates in the small branches of the river (left part of the delta—Figure 1). Thus, there is a longitudinal decrease in sediment concentrations (accumulation) for the Tumatskiy channel system. Similar processes are typical for the Olenekskaya channel, which is distinguished by the maximal channel length. This part of the delta is also characterized by high sedimentation rates. These differences are explained by the different morphology of the channel of the river. The Tumatskiy and Olenekskiy sectors are characterized by shallow and narrow branches with a significant drop of transporting capacity due to the longer length and associated surface slopes decline.

The erosion patterns are also dominated within Kolyma delta distributaries (Figure 5d). Here, the river flow almost entirely passes through two large main branches (Figure 1). In general, the delta is distinguished by a lower number of small channels. In the structure of the delta, there are almost no channels that lose their connection with the main riverbed network during the year; $\mathsf{\Delta }\mathrm{S}$ < 0 are noted only in conditions of low water content, when water consumption is less than 15,000 m³ s⁻¹. This pattern is reversed in comparison with the situations that are typical for the deltas of the southern regions (Selenga, Upper Angara); accumulation prevails only during the period of increased water content. The situation in which erosion processes prevail in the Kolyma delta (a positive balance of suspended sediments) is observed both under high and low waters. Under discharges higher then 15,000 m³ s⁻¹, longitudinal increase in suspended sediment concentration (SSC) is observed, which is opposite to depositional effects in the Southern deltas. This can be associated to the scarce vegetation on the surface of the floodplain, and thermal erosion of frozen permafrost banks. The influence of destructive processes is also observed during the spring flood due to an increase in the impact of the thermal and mechanical components of the water flow. During summer, an increase in both air and water temperatures leads to permafrost melting [59]. The ratio of positive and negative balances does not differ much between the two main branches of the river—Pokhodskaya and Kamennaya.

4. Discussion

4.1. Drivers of Sediment Transport Changes along Deltas

The observed results indicate that sedimentation processes over Northern and Southern deltas are driven by opposite factors. The intensity of deposition and erosion processes within deltas firstly depends on the relations between accommodation space provided by pre-existing drowned river-valley topographies, and hydrodynamic processes. Deltas are seen as a river transitions areas, where the energy available to transport suspended sediment decreases [60]. Submerging of banks and decrease of water slopes and transport capacity is particular important for suspended-sediment transport because water has no turbulent energy to suspend sediment, favoring deposition. Additionally, aquatic vegetation (macrophytes) dampens currents and turbulence [61], and thus also impacts sedimentation processes.

One of the aspects of the functioning of river deltas is the structure of the network itself—hydraulic connectivity, which, together with the resulting geometry of the islands, reflects the main processes shaping its system [62,63]. In Southern deltas, high flow events, which increase water levels, favor the hydraulic connectivity between the main channels and the adjacent wetlands [64] and, therefore, enhance sediment delivery from the main channel to the deltaic plain. Both Selenga and Upper Angara deltas submerged bank locations can accumulate various size fractions of sediment, even during moderate water discharge conditions [65]. In the Selenga delta, 10-fold increase of the water surface is observed under maximum discharge (4700 m³ s⁻¹), which leads to flooding over 60% of the delta surface (Figure 6b). Due to hindered connectivity between closed lakes and channels, the input of fine sediment to these lakes is facilitated during flooding events [55], particularly if a tie channel connects the lake and distributary channel [5]. This statement support observations at other deltas, such as Mekong River delta, which on average stores 28% of the incoming sediment load [66] due to high net accumulation rates during flooding dynamics, or Ganges-Brahmaputra Rivers Delta that stores 40% of total fluvial sediment discharge [6]. The Selenga delta can store up to 34% of the incoming suspended sediment load (67% of the total sediment load) during flow conditions considered as high (Q ~ 3000 m³ s⁻¹) [67]. Changes in the height of the delta surface in the years 1956–1998 [67] confirm the lack of a clear relationship between the main deposition zones and the vicinity of the main distributaries. Taking the Ganges-Brahmaputra Mega Delta as an example, the deposition mainly takes place close to the major distributors [68] and redirects ~10 ± 5% of the discharge to smaller distant regions of the delta, serving to compensate for slow subsidence [69]. On the Selenga delta, under a recent low-water period that started in the region after 1995 [70], due to the decreased frequency of high discharge events (Q > 1350 m³ s⁻¹), as well as the decreased magnitude of annual high flow events over the past two decades [9], floodplain lakes could remain disconnected from the delta’s main channels for longer time intervals. In particular, decreasing frequencies of high flow events (Q > 1350 m³ s⁻¹) were observed for the Selenga River System [71].

Submerging of the bank and near-shore areas is not typical for the Northern deltas. In particular, even high flow events in the Kolyma delta are not associated with bank submerging. Here along the main distributaries under maximum discharge (30,000 m³ s⁻¹) due to confined valleys, the flow is not extended over near-channel areas (Figure 4). This lead to the increase of stream energy during floods. The rapid rate of water level decrease in summer (even 50–60 cm per day) can have a significant impact on fluvial transport. In Lena, larger flooding is observed especially along smaller left distributaries where up to 5-fold increase of water surface under maximum discharge (120,000 m³ s⁻¹) occurs (Figure 7).

These differences are explained by the evolutionary history of the deltas. Both Lena and Kolyma deltas are much older drowned river-valley topographies, which thus provide smaller accommodation space for sediment deposition compared to actively formed Selenga and Upper Angara deltas. The Selenga and Upper Angara deltas are 500,000 years old [72] and recently experienced active tectonic subsidence [73]. Lena delta consists mainly of organomineral sediments, commonly called peat [74]. Lena and Kolyma deltas terraces were not formed as deposits made by the flooding river. i.e., river terraces [33]; rather, they were formed by Eurasia Basin of Arctic Ocean sea-level fluctuations over a long period of time, and thus, the surface is significantly higher, hampering possible submerging.

These differences are the primary explanation for the contrasting relationships between $\mathsf{\Delta }\mathrm{S}$ (%) и Q (m³ s⁻¹) of particular deltas (Figure 8). Deltas where the bank submerging is associated with increased sediment deposition demonstrate a mostly negative relationship between $\mathsf{\Delta }\mathrm{S}$ and water discharge Q. In both cases (Selenga and Upper Angara), the threshold value of discharge associated with bank flooding separates deposition processes that are typical for high flow events and both erosion patterns, which are observed only during low flow periods (Figure 8). These effects are also increased by summer vegetation growth in the Southern deltas, which marks an important change in sedimentary processes on estuarine margins during periods of high flow [61]. An increase of sediment concentrations during a low water period is caused by in-channel erosion and resuspension of former deposited sediments. Another factor influencing the increase of sediment load in the lower part of the Delta is Baikal wind surge, which effect sediment load during low water conditions. In the Northern deltas where the sediment transfer is driven by permafrost thaw (peak freshet) and unchanged by aquatic vegetation, the events associated with longitudinal increase of sediment concentrations are observed especially under higher discharges (Figure 8). Absence of sub-aqueous channel banks hampers the effects of hydraulic connection between channel and floodplain during both higher flood water discharge conditions. Thus, the effects of high discharge events, which result in shortened backwater lengths [73] and increased transport capacities, lead to increased erosion additionally to permafrost thaw effects.

The scale effects of particular distributaries play an important role in the mentioned differences between certain deltas (Figure 9). In particular, this is important for the largest Lena delta, which consists of distributaries of various size and morphology. As it was discussed in the results section, the erosion is dominant along the largest Eastern distributaries (Bykovskaya and Trofimovskaya), whereas in the smaller Tumatskaya and Olenekskaya distributaries, deposition (ΔS < 0) is more important compared to erosion (ΔS > 0). Following this, it is important to note that Figure 6 is related to the sediment concentrations and not to river flow (sediment fluxes). In such complex channel networks as the studied deltas, the general sediment budget of the whole system (in terms of volume of stored or eroded sediments) mostly depends on the processes within the main distributaries, which distribute the main part of water flow. In particular, in the Lena delta, even though the number of deposition patterns ∆S < 0 equals the number of ∆S > 0 patterns, the latter are observed in the largest channels. Thus, the total sediment budget will be mostly determined by the Bykovskaya and Trofimovskaya channels where ∆S > 0, and they contribute over 90% of water discharge.

Additional effects are associated with tides and backwater flows effects. A considerable area of the Selenga River delta (i.e., within ~ 9 km from the delta outlets; [73]) is impacted by backwater flow, where the range of water stage variability is limited and the entire channel perimeter is persistently wetted even during low water discharge conditions, thus favoring hydraulic connectivity between the main channels and the surrounding water bodies. Here, connectivity can be facilitated by tie channels and/or exchange of water at the channel–wetland interface, where levees of the main channel are developing and remain below the water line. In the Kolyma River, these effects can be seen under low discharges when there are abrupt declines in sediment concentrations in the upper part of the delta (Figure 10).

Our observations suggest that during wind events, there are significant backwater effects in the delta topset and opposite flows in the Kolyma tributaries (e.g., −180 m³ s⁻¹ were measured at Panteleiha river on 15.07.2021 with total loss of discharge at Kolyma delta topset (Cherskiy) 3500 m³ s⁻¹). Temporary loss of sediment transport capacity downstream and an additional increased re-suspension of sediments in the shallow parts of the cross-section can be attributed to a relationship between bed shear stress and water depth [75]. The transient increases and decreases of discharge (without expressive water level changes) mentioned above may be responsible for sudden changes in the transport conditions of the suspension [76]. The described phenomenon can be related to a kind of flood wave diffusion, where the influence of the cross-sectional shape on the spread of the single-peak flood wave is marked [77]. Sudden changes in the base flow observed in the estuary section of Kolyma may intensify towards the sea. Additional impacts on the sediment budget in the Arctic deltas are related to coastal upwelling events, which are regularly increased sea surface turbidity at certain areas adjacent to the Lena delta in the Laptev Sea and the Kolyma deltas in the East-Siberian Sea [78]. These events are formed under strong easterly and southeasterly wind forcing and are estimated to occur during up to 10%–30% of ice-free periods. Both with intensive coastal erosion [79], they increase sediment concentrations in the most downstream parts of the Lena distributaries (Figure 11).

4.2. Sedimentation and Biogeochemical Impacts

The obtained results are important also regarding understanding the biogeochemical role of the deltas. The Northern deltas due to active thermoerosional processes and the dominant role of large channel distributaries are primary erosional ones (

Table 3. Occurrence of various sediment budget patterns along studied deltas (number of events/% from total number of events).

Delta		ΔS < −20%	−20% < ΔS < 0%	0% < ΔS < 20%	0% < ΔS < 20%
Selenga		13/5.4	184/76.0	34/14.0	11/4.6
Upper Angara		91/52.3	60/34.5	10/5.7	13/7.5
Lena	all distributaries	24/23.8	34/33.6	23/22.8	20/19.8
	distributaries 1,2	5/8.9	14/25.0	21/37.5	16/28.6
	distributaries 3,4	19/42.3	20/44.5	2/4.4	4/8.8
Kolyma		5/5.6	18/20.0	40/44.4	27/30.0

The intensity of the color corresponds to the frequency of events, ΔS < 0 and ΔS > 0 are either depositional and erosional states, respectively.

). Over 70% of the observed events in the largest distributaries of Lena River delta (Bykovskaya and Trofimovskaya channels) and Kolyma River delta represent situations when ΔS > 0 in %. The Southern deltas might be considered as primary depositional systems due to sub-aqueous banks, abundant aquatic vegetation, and submerging of the wetlands. Over 80% of flow events of both Selenga and Upper Angara deltas are related to retention of suspended sediments along distributaries. Additionally, the total retention of the sediments is even more pronounced due to their capturing in the ponds located at the outlets of the deltas (Figure 1). This explain the downstream sediment fining of sediments over tens of kilometers of deltas, ranging from predominantly gravel (coarse pebble) and sand near its apex to silt and sand at the delta–lake interface. Both deltas are characterized by significant hydrodynamic control over sediment transfer along deltas, which is characterized by downstream elimination of gravel in distributary channels and is caused by declining boundary shear stress as a result of water discharge partitioning among the bifurcating channels [71] and decline of water slopes [65].

The observations within Southern deltas showed that the largest decrease in total metal concentrations occurred along smaller wetland-dominated channels. In the Selenga delta, reduction of the heavy metals concentration by 77–99%, during both herein considered moderate and high flow conditions (Q ~ 1000 m³ s⁻¹ and Q ~ 3000 m³ s⁻¹), were reported [65,67]. Submerged bank and connected marshland locations in the backwater zone may act as active sinks for all the various metals depending on the ambient geochemical conditions within the delta, such as the water pH, sediment mineralogy, OM and oxygen availability, as well as the related speciation of the metals [80,81,82].

More complicated systems are observed along Northern deltas, which reflect contrasting patterns of sediment transport. In the Kolyma river delta, specific observations were done in 2019 and relied on three sampling points at the three locations along right Kamennaya branch (Figure 1). The locations were placed in the delta topset, central part of the delta, and delta outlet [83]. At each location, the suspended sediment samples were taken from three depths (subsurface, middle, and near-bottom). Most of the metals concentrations (Al, V, Mn, Fe, Ni, Cu, Zn, As, Sr, Pb, and Mo) during the period of sampling demonstrated 10–30% decrease along the upper part of the delta, whereas in the downstream part of the delta, the concentration mostly increased up to 10–20% compared to the central part of the delta. The herein studied sedimentation processes, which were independent on the ambient flow conditions, can possibly explain these changes of metals concentrations.

Our proposed remote sensing and modeling approach to capture fluvial flux distribution and network dynamics along deltaic channels and coastal transport paths could be of major importance in integrating these methods into adaptive delta management strategies. The Landsat record is long enough to today quantify the changes that occur in the water distribution network, deltas, and dissolved particles within large river deltas.

5. Conclusions

Deltas are usually classified as depositional zones where, due to bifurcation, loss of stream power, and existence of abundant wetlands, large amounts of sediments are stored. Our study, which benefits from application of remote sensing data, reveals comprehensive spatio-temporal variations in suspended sediment budget within deltas of large rivers. In particular, we conclude:

• Patterns of sediment sinks and sources in the deltas significantly varied between various deltas and in general change due to latitudinal gradients.
• Zoning and seasonality might be considered in future studies as fundamental parameters of river delta sedimentation. Permafrost in the North and aquatic vegetation in the South are dominant drivers that determine zonal North-to-South variation in suspended sediment budget.
• Lena and Kolyma deltas represent mostly erosional patterns where suspended load increases along the main distributaries; here, an increase of sediment discharge in a permafrost-affected region may be predominantly driven by the rate of summer thawing and associated activation of erosion features. Upper Angara and Lena deltas mostly capture suspended sediments.

The performed study, due to a remote sensing approach, shines a light on previously unknown phenomena of river sedimentation related to suspended sediment behavior in deltas. Analyzing the hydrological connectivity and accumulation of sediment within river deltas is particularly important due to measured climatic multi-decadal trends in water discharges [84,85]. Understanding these processes is important when trying to find ways to maintain, restore, or even enhance a delta’s capacity to cope with the effects of climate change.