Modelling salt water intrusion into main branches of the Don delta depending on wind situation

Full Text

Introduction

The Don mouth area is a key region of the basin of the Sea of Azov, where complex processes of interaction between river and sea waters occur. Here, surge level fluctuations, associated with the flow of sea transformed waters into the branches of the Don delta, are especially pronounced [1, 2].

A serious threat of flooding is posed by extreme downsurges with easterly winds and upsurges with westerly winds. The situation is particularly dangerous when the water surges down before the upsurge. This occurs when the easterly wind suddenly changes its direction to the west. Thus, flooding occurs faster, and it is greater in scale than with a constant westerly wind. A similar picture was observed on 23 March 2013, when the easterly wind (3–11 m/s) changed to southwesterly (15 m/s with gusts of 20–22 m/s). This phenomenon is confirmed by the performed numerical experiments, in which the wind direction changes to the opposite one [3]. Such extreme upsurges are only possible with strong southwesterly winds. In addition, surges can be influenced by seiche fluctuations, which are significant in the Sea of Azov [4], but are not taken into account in the presented model.

Under certain hydrometeorological conditions, salt waters of Taganrog Bay intrude into the Don River delta, where the water intakes of the Rostov agglomeration largest cities, such as Azov and Taganrog, are located. Moreover, a new water pipeline was put into operation in the Donetsk region. During a strong upsurge in June 2014, the sea level rose to 1.7 m, and the salinity at the Don mouth reached 5 PSU [5, 6]. In September 2014, water with increased salinity (5–9 PSU) intruded along the Don to Azov and entered the city water supply system. Subsequently, the similar situation was observed several times, being the most significant in February 2021 [7].

In recent decades, we have observed large-scale climate changes, which led to significant changes in environmental conditions in the Don basin and the Sea of Azov. The Don low water has been continuing for 17 years and is the longest one in the last 100 years. This led to a significant decrease in freshwater flow into Taganrog Bay of the Sea of Azov and an increase in salinity both in the sea and in Taganrog Bay. The average sea water salinity increased to 14 PSU (in the 1970s, such salinity was considered critical for the marine ecosystem), and in Taganrog Bay it reached 10 PSU [8]. Analysis of field data from the oceanological research of the Southern Scientific Center of the Russian Academy of Sciences indicates a continuing increase in salinity. The decrease of the Don flow to historical minimums, as well as the anomalous advection of the Black Sea waters led to the abnormally high salinity observed in Taganrog Bay during 2014–2016 (up to 12 PSU) [2]. It should be noted that the current rate of increase in the average annual sea water salinity exceeds the rate observed in the mid-1970s, during the previous period of the Don low water levels [9]. In this regard, the study of the intrusion of salt sea water into the Don River delta is becoming increasingly relevant.

Currently, many papers are aimed at modeling the intrusion of salt water into river mouths. These studies are mainly related to tidal processes. Thus, the model of salt inflow into the mouth of the San Francisco River is described in [10]. In this case, a numerical solution of advection–diffusion problems in hydrodynamics on an unstructured grid is used [11].

A numerical model for studying seasonal variability of currents and salinity at the Indus River mouth in Pakistan is given in [12]. The model is calibrated using observed water level, current velocity, and salinity data. The modeling results show that salt water intrudes far upstream, approximately 65 km.

The Don mouth is not subject to any tidal processes. The intrusion of salt water into the Don main channels occurs during upsurges, when the water level rises significantly in Taganrog Bay. Therefore, the intrusion of salt water into the Don main branches is directly related to the salinity of water at the mouth of the branches, as well as the movement of water in the channels. Here, the flow of water in the river channel is influenced primarily by the water level in the receiving basin, in this case in Taganrog Bay.

The paper is aimed at computational study of this phenomenon. Use of a mathematical model based on equations for incompressible liquid motion and equations for convection–diffusion (transfer) allows studying saline water intrusion into the Don delta during its flooding.

For these reasons, when modeling the process of intrusion of salt water into the Don main channels, the following steps are necessary to calculate the required parameters.

1. Determination of the water level at the mouths of the branches for a given wind situation over the Sea of Azov and Taganrog Bay.
2. Determination of the salinity of water at the mouths of the branches.
3. Determination of the velocity of water movement and distribution of salinity in the Don main branches.

The first two problems are solved using a two-dimensional model of the hydrodynamics of wind currents in the Sea of Azov, as well as a model of salinity distribution in it [13]. The third problem is solved using one-dimensional models of water movement in channels and substance transport [14].

Materials and methods

Since the beginning of the 2000s, during marine expeditions, the Southern Scientific Center of the Russian Academy of Sciences has been conducting systematic observations of the thermohaline structure of the waters of the Don River coastal estuary and delta. Since the summer of 2014, studies of the dynamics of salinity and water temperature changes, together with measurement of the current direction and velocity, have been regularly carried out using an integrated current meter Aanderaa RSM-9LW (https://www.aanderaa.com/media/pdfs/Seaguard_RCM-TD262b_001.pdf) at a stationary buoy station on the seaside 5 km from the edge of the delta (the first leading mark of the Azov-Don Seaway Canal) and at the mouth of the Don (the khutor of Donskoy, the village of Kagalnik), as well as in the Don delta at the network of ship oceanological stations. In parallel with this, meteorological parameters and water level have been observed at the hydrological station of Donskoy and water level station of Taganrog.

To study the spatial distribution of sea water temperature and salinity in the Sea of Azov and Taganrog Bay, continuous recording of data on the temperature and electrical conductivity of surface water was performed with thermosalinograph SBE21 SEACAT (Available at: https://www.seabird.com/sbe-21-seacat-thermosalinograph/product?id=60762467702) during the voyages of the R/V Deneb.

The area of calculation is the mouth section of the Lower Don from the stanitsa of Razdorskaya to Taganrog Bay, including its eastern part. This section consists of the main channel of the Don and its main branches, the Old Don, the Bolshaya Kalancha turning into the Mokraya Kalancha and the Bolshaya Kuterma (Fig. 1). White circles indicate points in Taganrog Bay where calculated water level values were taken, black triangles indicate the hydrological stations where observations were made.

 

Fig. 1. The area of calculation of the Don River mouth section. The inset shows the diagram of the delta part: 1, 3, 5, 6 – end nodes; 2, 4 – branch nodes

 

The salinity transport in the Sea of Azov is described by a system containing equations for long waves in a homogeneous incompressible fluid in a Coriolis force field and a transport equation under the assumption that the propagated substance is conservative [15]:

 $\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}-f\cdot v=-g\frac{\partial \zeta }{\partial x}+\frac{{\tau }_{sx}}{H}-\frac{{\tau }_{bx}}{H},$

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+f\cdot u=-g\frac{\partial \zeta }{\partial y}+\frac{{\tau }_{sy}}{H}-\frac{{\tau }_{by}}{H},$ (1)

$\frac{\partial \zeta }{\partial t}+\frac{\partial \left(Hu\right)}{\partial x}+\frac{\partial \left(Hv\right)}{\partial y}=\mathrm{0,}$

$\frac{\partial c}{\partial t}+\frac{\partial \left(uc\right)}{\partial x}+\frac{\partial \left(vc\right)}{\partial y}={\epsilon }_{xy}\left(\frac{{\partial }^{2}c}{\partial x^2}+\frac{{\partial }^{2}c}{\partial y^2}\right),$

where H = h + ζ; h = h(x, y) – depth; f – Coriolis parameter; ζ – water level difference; u = u(x, y, t), v = v(x, y, t) – velocities; c – concentration; ε_xy – horizontal turbulent diffusion coefficient; τ_sx, τ_sy – projections on the OX and OY axes of the force of wind friction on the surface of the basin; τ_bx, τ_by – projections on the OX and OY axes of the force of wind friction on the bottom. These values depend on the velocity of wind W_B = {W_x; W_y} and current W_T = {u; v} and are determined as follows [16]:

${\tau }_s=\gamma \left|{\text{W}}_B\right|{\text{W}}_B$, ${\tau }_b=\beta \left|{\text{W}}_T\right|{\text{W}}_T$,

where $\left|{\text{W}}_B\right|=\sqrt{W_x^2+W_y^2}$, $\left|{\text{W}}_T\right|=\sqrt{u_{}^2+v_{}^2}$, ϒ – coefficient of the friction of wind on a free surface, β – coefficient of the friction of liquid on the bottom.

Sliding conditions are set along the solid boundary or the rates of inflow or outflow of water are specified (for example, for river mouths). It is assumed that there is no flow through the free surface and the side boundary.

In accordance with [17], the Don delta is represented as a graph (inset in Fig. 1), consisting of five edges and six vertices. The edges correspond to sections of open channels, i. e. the Don main channel and its branches. Four vertices correspond to end nodes (1, 3, 5, 6), and two – to branch nodes (2, 4).

The intrusion of salt water into the Don branches is described by one-dimensional equations of water movement in an open channel and the conservative substance transport [18]. No zones of sharp changes in the channel cross section, as well as distributed lateral inflow due to its insignificance, are assumed to be in the studied section of the Don channel. In the case when the channel cross section has a parabolic profile, this system can be rewritten in the following form:

$\frac{\partial Q}{\partial t}+gW\left(\frac{\partial z}{\partial x}+\frac{Q\left|Q\right|}{K^2}\right)=0$,

$b\frac{\partial z}{\partial t}+\frac{\partial Q}{\partial x}=0$, (2)

$\frac{\partial S}{\partial t}+v\frac{\partial S}{\partial x}-\mu \frac{{\partial }^{2}S}{\partial x^2}=0$,

where x – coordinate; t – time; Q – water consumption; z – water level; W – cross section area; K – rate of discharge; g – acceleration of gravity; b – channel width; S – concentration; v – velocity of water movement in the channel; m – turbulent diffusion coefficient. Rate of discharge К is calculated by formula K = ω·C$\sqrt{R}$, where R – hydraulic radius; С – Chezy’s velocity factor found using Manning formula $C=\sqrt[6]{R}/n$, n – bottom roughness.

The first two equations of system (2) (hydrodynamic component) are reduced to the following characteristic equations:

$\frac{gW}{c_{*}}\frac{\partial z}{\partial t}+\frac{\partial Q}{\partial t}+gW\frac{\partial z}{\partial x}+c_{*}\frac{\partial Q}{\partial x}=-gW\frac{Q\left|Q\right|}{K^2}$,

$-\frac{gW}{c_{*}}\frac{\partial z}{\partial t}+\frac{\partial Q}{\partial t}+gW\frac{\partial z}{\partial x}-c_{*}\frac{\partial Q}{\partial x}=-gW\frac{Q\left|Q\right|}{K^2}$,

where ${с}_{*}=\sqrt{gW/b}$.

As boundary conditions, discharge of incoming water Q₀ (0, t) is set at the initial point (node 1), and level in the receiving basin z_k (X_k, t) is set at the end points (nodes 3, 5, 6). At branching nodes, the boundary conditions are set as follows: the sum of inflows and outflows is equal to zero ${\sum }_i^{}{Q}_i=0$ (i – number of branches coming to the branching node), and the water surface levels are equal to each other: z_i = z.

The boundary conditions are set at the ends of the branches depending on the sign of the current velocity. If water enters the branch, the following condition is set: $\frac{\partial S}{\partial x}=\frac{v}{\mu }\left(S-S*\right)$, where S * – salt concentration value at the end of the branch. In the case when water leaves the branch, the following condition is set: $\frac{\partial S}{\partial x}+S\frac{v}{\mu }=0$, which corresponds to the disposal of salt from the branch.

Thus, during strong upsurges at the boundary of the branches coinciding with the mouths in Taganrog Bay, the following condition is set: $\frac{\partial S}{\partial x}=\frac{v}{\mu }\left(S-S_{zal}\right)$, where S_zal – salt concentration value in Taganrog Bay. If the end of the branch contains lumped parameters, in particular concentration, S_zal should be replaced with S * – concentration at the branch point, which is calculated as follows:

$S*=\frac{{\displaystyle \sum S_i{Q}_i}}{{\displaystyle \sum Q_i}}.$

Summing-up is performed only over the branches flowing into the branch point.

The problem is solved by finite-difference methods using implicit schemes. In the “marine” model (1), a grid with steps Δx = 660 m, Δy = 685 m is constructed, and in the “channel” model (2), the step is Δx = 1000 m. Stable counting is observed at time step Δt = 30 s. At each time step, hydrophysical parameters are first calculated for the entire Sea of Azov, then the values of salinity and water level are selected at the points corresponding to the mouths of the Don branches. These points are indicated by circles in Fig. 1. Next, the calculation continues in the Don channel area. Then, a transition to a new temporary layer is performed.

Results of calculation and discussion

Based on the salinity observations, calculations were performed under corresponding wind situations. Coastal weather stations located around the Sea of Azov provided the wind data. For each calculated point, three nearest weather stations were determined, and then, using linear interpolation, zonal and meridional coordinates of the wind velocity were determined in this point.

From 23 September 2014 to 27 September 2014, an extreme upsurge accompanied by salt water intrusion into the Don branches, was observed [6]. At the beginning of the period, an easterly wind of 2–4 m/s took place. From 10 a. m. on 24 September 2014, the wind changed its direction to southwesterly and increased to 20 m/s. This led to a sharp rise in water and increased salinity. The maximum values of water level and salinity in the area of the port of Azov occurred on September 24–25 and amounted to 3.82 m and 5.59 PSU, respectively. Fig. 2 shows changes in water salinity and water level in the water area of the port of Azov. Calculations confirm the fact that changes in salinity are directly proportional to changes in water level. The calculation error is 19.1 % for salinity and 22.8 % for water level.

 

Fig. 2. Water salinity and water level in the water area of the port of Azov from 24 September 2014 to 27 September 2014. The triangles are observed concentration; the solid line is calculated concentration. The diamonds are the observed level; the dashed line is the calculated level

 

One of the last relatively strong upsurges with the salt water intrusion into the Don branches was observed on 12–16 February 2021. During this period, water rise was observed twice. From 03:00 p. m. on 12 February 2021, influenced by a westerly wind of 13–16 m/s, the first rise of water took place. Then, as the wind weakened to 4–6 m/s, the level decreased. But starting from 12:00 p. m. on 13 February 2021, the wind again rose up to 12–16 m/s, which led to a new water rise. Then the wind weakened to 1–2 m/s and took a northwesterly direction.

For a second time, sampling was performed in Rostov-on-Don (city beach) and Azov (port). In such points of Rostov-on-Don as Nakhichevan Duct, khutor of Koluzaevo, khutor of Dugino, khutor of Rogozhkino, sampling was performed only once. Fig. 3 shows changes in salinity in the area of the Rostov beach and the port of Azov, as well as in the current velocity in the area of the port of Azov.

The table shows the time of water sampling, calculated and observed salinity, as well as calculation error.

The change in salinity depends directly on fluctuations in the water level at the mouths of the branches, which is quite natural. It can be seen that the change in salinity in the port of Azov, which is located 15 km from the mouth, occurs with larger amplitudes than in the area of the Rostov beach, located 50 km from the same mouth (Fig. 3). It can also be seen that salinity increases almost immediately after the current velocity becomes negative (reverse current), and decreases when the direct current returns.

 

Fig. 3. Salinity and velocity of current in the area of the port of Azov from 13 February 2021 to 15 February 2021. The triangles are the Rostov beach (observed). The dashed line is the Rostov beach (calculated). The diamonds are the port of Azov (observed). The solid line is the port of Azov (calculated). The dotted line is the velocity of current (the port of Azov)

 

From 22 November 2022 to 23 November 2022, measurements of salinity and current velocity were performed at the hydrological station of Donskoy. At the beginning of this period, a southwesterly wind was observed with a force of 4–5 m/s, which then intensified to 12–14 m/s, resulting in an upsurge with the salt water intrusion into the Don delta. A comparison of the results of calculating salinity and current velocity with the observed data is shown in Fig. 4. A noticeable increase in salinity begins at the moment when the current changes direction to the opposite, from the sea to the river, and the velocity becomes negative (Fig. 4).

 

Fig. 4. Salinity and velocity of current at the hydrological station of Donskoy from 22 November 2022 to 23 November 2022. The black lines are salinity; the grey lines are velocity; the solid lines are calculations; the dotted lines are observations

 

The calculation errors are 27.7 % for water level, 16.8% for salinity, 92.6 % for current velocity. The large current velocity calculation error is explained by the assumption that the bottom profile is parabolic, which does not fully correspond to reality in this hydrological section.

When comparing the observed and calculated values, attention is drawn to the large error in the calculations of salinity in the khutors of Dugino and Rogozhkino located in the Bolshaya Kalancha branch, while at other points located in the Don main channel and in the Old Don branch, this error is acceptable (table).

 

Comparison of calculated salinity values (PSU) with observations

Sampling location	Distance from the mouth, km	Date, time	Salinity		Error, %
			Observed	Calculated
Rostov-on-Don (Nakhichevan Duct)	52.3	13.02.2021 14:54	0.79	0.756	4.30
Rostov-on-Don (city beach)	49.4	13.02.2021 15:41	0.79	0.768	2.78
Rostov-on-Don (city beach)	49.4	14.02.2021 12:27	0.83	0.760	8.43
khutor of Koluzaevo	34.5	14.02.2021 16:35	0.85	0.694	18.35
khutor of Dugino	22.3	13.02.2021 17:48	3.24	0.768	76.30
khutor of Rogozhkino	11.7	13.02.2021 17:22	3.97	1.023	74.23
city of Azov (port)	15.0	13.02.2021 14:36	2.56	2.534	1.02
city of Azov (port)	15.0	14.02.2021 12:15	5.28	4.72	10.61

 

To clarify the reasons of the above mentioned, additional hydrodynamic calculations were performed in the Don delta region in November 2022 based on the wind situation for a moderate upsurge. At that time, measurements of the velocity and direction of currents and water salinity were performed from the board of the R/V Deneb in the Old Don branch (hydrological station of Donskoy, 7.2 km from the mouth) and in the Bolshaya Kalancha branch (down the hydrological station of Dugino, 18 km from the mouth). The change of the water level in the delta during the upsurge was measured by a level gauge at the hydrological station of Donskoy. Fig. 5 shows the observed and calculated values of the water level for the observation period from 17 November 2022 to 21 November 2022.

 

Fig. 5. Change in the water level at the hydrological station of Donskoy from 17 November 2022 to 21 November 2022. The dotted line is observations. The solid line is calculations

 

Fig. 6 shows changes in calculated and observed current velocities and salinity at the hydrological station of Donskoy for the period from 17 November 2022 to 21 November 2022. Moreover, the observations here are presented only from 07:00 p. m. on 17 November 2022 to 07:00 a. m. on 18 November 2022, and then the values obtained from the calculation based on the wind situation for the specified period, are given.

 

Fig. 6. Calculated and observed values of the current velocity and salinity at the hydrological station of Donskoy from 17 November 2022 to 21 November 2022. The dashed line is observed salinity; the solid grey line is calculated salinity. The dotted line is observed velocity; the solid black line is calculated velocity

 

Thus, a clear connection between the direction of the current and the change in salinity is observed: at a negative velocity, when there is a reverse current directed from the sea to the river, salinity increases, and at a positive velocity, salinity decreases. This connection is confirmed by observations (Fig. 6).

From 10:00 p. m. on 17 November 2022, the current velocity decreased, remaining positive until approximately 02:30 a. m. At that time, the salinity values were close to 0.60 PSU. Then, a reverse current began to develop, and at 05:00 a. m. on 18 November 2022, a sharp increase in salinity occurred. The time difference between the beginning of the reverse current and the jump in salinity is explained by the fact that Donskoy is located 6.5 km from Taganrog Bay, and it takes time for salt water to reach the observation point. It can be seen (Fig. 7) that during the reverse current, starting from 05:00 a. m. on 18 November 2022, the salinity increases sharply. This figure is an enlarged copy of the initial part of the graph presented in Fig. 6.

 

Fig. 7. The dependence of salinity at the hydrological station of Donskoy on the current velocity from 17 November 2022 to 18 November 2022 (observed values). The solid line is observed salinity. The dashed line is observed velocity

 

Fig. 8 shows change in calculated and observed values of current velocity during the observation period from 10:00 p. m. on 17 November 2022 to 07:00 a. m. on 18 November 2022. The dynamics of the calculated and observed values agree well with each other. However, there is a discrepancy between observation data and modeling results at the night period, when the observed level values are oscillatory. This can be due to an error in the wind field interpolation or the occurrence of proper oscillations.

 

Fig. 8. Calculated and observed values of the current velocity at the hydrological station of Donskoy from 10 p. m. on 17 November 2022 to 7 a. m. on 18 November 2022. The dotted line is observed velocity. The solid line is calculated velocity

 

At the observation point in the area of Kostina Yama, the relationship between the current direction and the change in salinity is not as clear as in Donskoy. In fact, in this section of the Bolshaya Kalancha channel, changes in the current velocity and direction during upsurge conditions are instrumentally recorded and coincide satisfactorily with the calculated values, but almost no change in salinity is observed (Fig. 9).

 

Fig. 9. The dependence of salinity on the velocity of current in the area of Kostina Yama from 18 November 2022 to 19 November 2022 (observed values). The dotted line is observed salinity. The solid line is observed velocity

 

This is explained by the fact that observation points Donskoy and Kostina Yama are located at different distances from the mouth of the corresponding branches of the delta. In addition, the distribution of flow along the branches is also uneven, as the Old Don accounts for 40% of the Don flow, and the Bolshaya Kalancha – 60%. And finally, the longitudinal profile of the Old Don channel is more uniform, as the Azov-Don Seaway Canal passes through this branch, and dredging works are regularly performed resulting in favorable conditions for the salt water intrusion. At the same time, the profile of the Bolshaya Kalancha channel is characterized by shallower average depths and alternating river reaches and rifts, which complicates salt water the intrusion into this branch.

Conclusions

Correct setting of the initial salinity field in Taganrog Bay is a challenging part of the problem under consideration. Further distribution of salinity values at the mouth coastal area and then in the Don branches, depends on it, which significantly affects the calculation results. Based on the data obtained during spring, summer, and autumn surveys, salinity values were calculated along the route of the R/V Deneb from the Don mouth in Taganrog Bay to the Kerch Strait. Using the obtained data, the Sea of Azov and Taganrog Bay were divided into the areas where salinity was considered constant. In the absence of such data, it is possible to use corresponding salinity maps for a certain time of year. This will make it possible to obtain the initial distribution of water salinity with some approximation. However, this approach cannot ensure an acceptable initial distribution of salinity during the time period when the calculation is performed, especially for a future forecast.

This model makes it possible to determine the current velocity, water surface level, and salt concentration in the Don delta main branches depending on the wind situation in the Sea of Azov and Taganrog Bay.

The comparison of the calculated water level values with the values obtained at observation stations shows identical dynamics of their change, which indicates the adequacy of the presented model. The proposed methodology can be applied to other water areas in the joint calculation of currents in mouth channels and wind currents in a receiving basin.

About the authors

Aleksey L. Chikin

Institution of Science "Federal Research Center Southern Scientific Center of the Russian Academy of Sciences"

Author for correspondence.
Email: chikin1956@gmail.com
ORCID iD: 0000-0002-4065-010X
SPIN-code: 3845-9760
Scopus Author ID: 8240627300

Chief Research Associate, Dr.Sci. (Phys.-Math.), Senior Research Associate

Russian Federation, 344006, Rostov-on-Don, Chekhova Ave., 41

Aleksey V. Kleshchenkov

Institution of Science "Federal Research Center Southern Scientific Center of the Russian Academy of Sciences"

Email: chikin1956@gmail.com
ORCID iD: 0000-0002-7976-6951
Scopus Author ID: 57016697100
ResearcherId: E-6619-2014

Leading Research Associate, Ph.D. (Geogr.)

Russian Federation, 344006, Rostov-on-Don, Chekhova Ave., 41

Lyubov G. Chikina

Federal State Autonomous Educational Institution of Higher Education "Southern Federal University"

Email: chikin1956@gmail.com
SPIN-code: 4295-3950

Professor, Dr.Sci. (Phys.-Math.), Associate Professor

Russian Federation, Rostov-on-Don

References

Supplementary files