Resuspension of Bottom Sediments in a Shallow Lagoon by Currents and Waves Based on the Numerical Modeling Data (Using the Example of Sivash Bay, the Sea of Azov)

Abstract

Purpose. The work is purposed at studying the intensity of resuspension of silty bottom sediments in Eastern Sivash Bay (the Sea of Azov) during an extreme storm, as well as at assessing the contribution of currents and wind waves to the resuspension processes.

Methods and Results. The current fields are calculated using a three-dimensional σ-coordinate water circulation model of the POM type supplemented with a block of silty sediments resuspension. The SWAN spectral model is applied to calculate wind waves. In both models a rectangular computational grid with the horizontal resolution 300 m is involved. The ERA-Interim atmospheric reanalysis data corresponding to the extreme storm situation in November 10–13, 2007 are used as a forcing. The performed calculations constituted a base for analyzing the structure of the fields of waves, currents, bottom shear stresses and suspended matter concentration in Eastern Sivash for different phases of the storm. A technique for assessing the resuspension model sensitivity to the variations in the input parameter values is proposed.

Conclusions. The applied resuspension model is most sensitive to the variations in the parameter values that condition intensity of the silt particles vertical flow from the basin bottom. During the period of the storm maximum development, conditions for forming resuspension zones arise on 80 % of the total area of Eastern Sivash Bay. If, while modeling, the contribution of the waves is not taken into account, the total area of resuspension is reduced by four times. This fact testifies to a decisive contribution of the bottom wave stresses in formation of the resuspension zones in bottom sediments in the bay.

Full Text

Introduction

The need to study water dynamics of Sivash Bay of the Sea of Azov is associated with active anthropogenic transformations of this water area of the Crimean coastal zone. The bay is a recoverable deposit of mineral salts, it is characterized by high biodiversity and included in the list of wetlands of international importance. A special economic zone “Sivash” with a developed chemical industry and agriculture is located on the coast of the bay [1].

The present-day Sivash is a vast shallow bay of the Sea of Azov (a lagoon-type sea bay) with an indented coastline. Sivash is usually divided into two large regions – Western and Eastern Sivash (Fig. 1, a). Western Sivash is currently an almost completely isolated basin, its water level is regulated by a dam, its depth usually does not exceed 0.3–0.4 m. Eastern Sivash is a deeper (up to 3 m) section of the bay, occupying over 60% of the total area of Sivash.

F i g. 1. Location of East Sivash Bay (a), model bottom relief of East Sivash Bay (b)

Several publications are devoted to the study of hydrodynamic processes in Eastern Sivash. The work [2] presents calculations of surface currents and level oscillations in the bay under winds of various directions. The characteristics of a freshwater plume in the area where the Salgir River flows into Sivash Bay are studied in [3]. In [4], current fields in Eastern Sivash are modeled for different wind conditions. It is demonstrated that the structure of water circulation in the bay is determined by the general direction of the wind. The field of depth-average currents consists of several eddy structures. These structures are expressed especially well in the southern part of the water area. When winds are directed along the water area of the bay, the circulation is most intense. In these cases, in the narrowness connecting the northern and southern parts of Eastern Sivash, a jet current with a speed of more than 1 m/s occurs. In the deepest parts of the water area, a bottom countercurrent is formed near the bottom.

In [5], modeling of the wind wave characteristics in Eastern Sivash was carried out. It is revealed that due to the shallowness of the water area and relatively small effective fetches, the main factor affecting the intensity of waves is wind speed. The waves reach a steady state 3 h after the start of the wind effect. After the wind stops, the waves completely cease after 2 h. Estimates of the maximum values of wave parameters for different wind speed gradations showed that even during the periods of strong storm, the heights and periods of waves in the waters of Eastern Sivash do not exceed 1 m and 3 s, respectively. In these works, the problems of the effect of currents and waves on the dynamics of bottom sediments in the waters of Eastern Sivash remained beyond the scope of the studies.

The problems of mathematical modeling of suspended matter evolution in shallow waters are considered in [6, 7]. The work [6] simulates the processes of bottom sediment resuspension and sedimentation at the Black Sea northwestern shelf during the passage of an atmospheric cyclone with no regard to wind waves. Bottom sediments are assumed to consist of a single silt fraction. The areas of the most intense rise of suspended matter and vertical profiles of suspended bottom sediment concentration in individual areas are analyzed. A similar study for the eastern part of the Gulf of Finland, but with regard to wind waves, was carried out in [7]. Here, the calculations were carried out for two fractions of bottom sediments: fine-grained sand and silt. The contribution of currents and waves to the sediment resuspension intensity is studied. The importance of taking into account the contribution of wind waves to resuspension processes is shown.

This work is purposed at mathematical modeling of the process of bottom sediment resuspension in Eastern Sivash. Mathematical modeling is an important part of environmental monitoring of the water area; it allows one to calculate various parameters of the bay state, which are not always possible to obtain experimentally. The paper presents the results of numerical experiments on modeling the resuspension of silty sediments during an extreme storm. A description of the mathematical models we applied is given below. The main attention is paid to the analysis of the contribution of currents and waves to the resuspension intensity.

The major problem of mathematical modeling of bottom sediment resuspension is the lack of the required amount of initial data. As a rule, the granulometric composition of bottom sediments is known fragmentarily; there is not enough field observation data to validate models; the selection of critical values of bottom shear stresses determining resuspension processes is a separate task. Currently available field studies in Eastern Sivash [8, 9] do not provide a complete understanding of suspended matter dynamics in the bay, since they are episodic and carried out directly offshore in relatively calm weather. For this cause, the bottom sediment resuspension model sensitivity to changes in input parameters was analyzed.

Materials and methods

Wave model. To calculate wind wave fields in Eastern Sivash Bay, we apply the Simulating Waves Nearshore (SWAN) model [10], based on the numerical solution of the wave energy balance equation in spectral form

$\frac{\partial }{\partial t}N+\frac{\partial }{\partial x}\left(c_{X}N\right)+\frac{\partial }{\partial y}\left(c_{Y}N\right)+\frac{\partial }{\partial \omega }\left(c_{\omega }N\right)+\frac{\partial }{\partial \theta }\left(c_{\theta }N\right)=S/\omega$, (1)

where $N=E/\omega$ is wave action density; E is wave energy spectrum; x, y, t are spatial coordinates and time; $\omega ,\theta$ are frequency and angular coordinates; $c_X,c_Y,c_{\omega },c_{\theta }$ are rates of wave action density transfer along spatial and frequency-angular coordinates.

The source function on the right side of equation (1) has the following form:

$S=S_{\text{in}}+S_{\text{nl}}+S_{\text{wc}}+S_{\text{bf}}+S_{\text{dib}}$, (2)

where $S_{\text{in}}$ is a source of wind-induced waves; $S_{\text{nl}}$ are nonlinear interactions of spectral harmonics; $S_{\text{wc}}$ is energy dissipation due to wave crests breaking; $S_{\text{bf}}$ is energy dissipation due to bottom friction; $S_{\text{dip}}$ is wave breaking at critical depths. Details of parameterization of individual source function terms (2) are given in the manual ¹.

The model bottom relief of Eastern Sivash Bay is represented in Fig. 1, b. To calculate wave characteristics, we use a rectangular grid with a horizontal resolution of Δx = Δy = 300 m (284 × 334 nodes). The discreteness of the model along the angular coordinate is 10°. An uneven grid with 31 nodes is applied along the frequency coordinate, varying within 0.04–1.0 Hz range. We use a non-stationary version of the SWAN model. Integration over time is carried out using an implicit difference scheme with 30 min step.

The outputs of the SWAN model are fields of significant wave height $h_{\text{s}}$, mean wave direction ${\text{θ}}_{\text{w}}$, peak wave period $T_{\text{p}}$ and amplitudes of near-bottom orbital motion velocity $U_{\text{w}}$. These parameters are applied to calculate bottom shear stresses in the model of currents and the bottom sediment resuspension model.

Model of currents. In order to calculate current fields, we used a barotropic version of the three-dimensional hydrostatic water circulation model proposed in [11]. The model is based on three-dimensional $\sigma$-coordinate equations of hydrodynamics in the hydrostatic approximation (hereinafter, summation is carried out over repeating indices $\alpha$ and $\beta$ from 1 to 2):

$\frac{\partial }{\partial t}\left(D{u}_{\text{α}}\right)+\Lambda u_{\text{α}}+{\text{ε}}_{\text{αβ}}fD{u}_{\text{β}}+gD\frac{\partial \text{η}}{\partial x_{\text{α}}}=\frac{\partial }{\partial x_{\text{β}}}\left(D{\text{τ}}_{\text{αβ}}\right)+\frac{\partial }{\partial \text{\hspace{0.17em}σ}}\left(\frac{K_{\text{M}}}{D}\frac{\partial u_{\text{α}}}{\partial \text{σ}}\right),$ (3)

$\frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x_{\alpha }}\left(D{u}_{\alpha }\right)+\frac{\partial w_{\ast }}{\partial \sigma }=0$, (4)

$\Lambda \varphi =\frac{\partial }{\partial x_{\beta }}\left(D{u}_{\beta }\varphi \right)+\frac{\partial }{\partial \sigma }\left(w_{*}\varphi \right)$, ${\text{τ}}_{\text{αα}}=2A_{\text{M}}\frac{\partial u_{\text{α}}}{\partial x_{\text{α}}}$, ${\text{τ}}_{\text{αβ}}={\text{τ}}_{\text{βα}}=A_{\text{M}}\left(\frac{\partial u_{\text{β}}}{\partial x_{\text{α}}}+\frac{\partial u_{\text{α}}}{\partial x_{\text{β}}}\right)$,  (5)

where $(x_1,x_2)$ = $(x,y)$; $\sigma$ is dimensionless vertical coordinate varying from –1 to 0; (u₁, u₂) = (u, v) are current velocity components along the axes $x_1,x_2$; w is current velocity component directed along the normal to the surfaces $\sigma$ = const; $D=h+\eta$; h is basin depth; $\eta$ is free surface coordinate; f is the Coriolis parameter; ${\text{ε}}_{\text{αβ}}=0$ at $\text{α}=\text{β}$ and ${\text{ε}}_{12}=-1$, ${\text{ε}}_{21}=1$; ${\text{τ}}_{\text{αβ}}$ are components of the turbulent stress tensor; $A_{\text{M}},K_{\text{M}}$ are coefficients of turbulent viscosity; g is gravitational acceleration.

At the solid lateral boundaries of the computational domain, no-slip conditions are set for velocities and conditions for turbulent fluxes to be equal to zero. To simplify the problem, water exchange through the Tonkiy Strait is not taken into account.

On the free surface of the basin ($\text{σ}=0$), the boundary conditions have the following form:

$w_{*}=0$, $\text{ρ}\frac{K_{\text{M}}}{D}\frac{\partial u_{\text{α}}}{\partial \text{σ}}=\text{ρ}{}_0с_{0}W\cdot W_{\text{α}}$, (6)

where ${\text{ρ}}_0$ is air density; $c_0=(0.49+0.065W)\cdot {10}^{-3}$ is surface friction coefficient; $W=\sqrt{W_1^2+W_2^2}$ is modulus of near-water wind speed; $W_1$, $W_2$ are zonal and meridional wind speed components.

At the bottom ($\text{σ}=-1$), the boundary conditions are set by the following expressions:

$w_{*}=0$, $\text{ρ}\frac{K_{\text{M}}}{D}\frac{\partial u_{\text{α}}}{\partial \text{σ}}={\text{τ}}_{\text{с}}\frac{u_{\text{α}}}{U_{\text{c}}}$, (7)

where $U_{\text{c}}=\sqrt{u_1^2+u_2^2}$ is modulus of horizontal velocity of currents; ${\text{τ}}_{\text{с}}$ is near-bottom shear stress determined by the formula

${\text{τ}}_{\text{с}}=\text{ρ\hspace{0.05em}}{c}_{\text{b}}{U}_{\text{c}}^2$, (8)

where $c_{\text{b}}=0.16/{\text{ln}}^{\text{2}}\text{(δ\hspace{0.05em}}z\text{/}{z}_{\text{b}}\text{)}$ is bottom friction coefficient; $\text{δ\hspace{0.05em}}z$ is vertical distance from the bottom to the point where the friction coefficient is determined; $z_{\text{b}}=d_{50}/30$ is bottom surface roughness parameter; $d_{50}$ is average diameter of bottom soil particles corresponding to silt deposits.

To determine the coefficient of vertical turbulent exchange $K_{\text{M}}$, the Mellor–Yamada model is applied [12]. The horizontal turbulent exchange coefficient $A_{\text{M}}$ is calculated using the Smagorinsky formula [13].

The boundary value problem (3)–(8) is solved numerically on the basis of explicit difference schemes for horizontal coordinates and implicit difference schemes for vertical coordinates. To approximate  advection operator, monotonic difference schemes are used. A detailed description of numerical algorithm is given in [11]. For horizontal coordinates, the same computational grid as in the SWAN model is applied. Along $\sigma$-coordinate, 11 uniform computational levels are specified.

Model of bottom sediment resuspension. The model of bottom sediment resuspension is incorporated into the water circulation model and is based on numerical solution of transport-diffusion equation of the following form:

$\frac{\partial }{\partial t}\left(DC\right)+\Lambda C-\frac{\partial \left(w_{\text{s}}C\right)}{\partial \text{σ}}=\frac{\partial }{\partial x_{\text{β}}}\left(A_{\text{C}}D\frac{\partial C}{\partial x_{\text{β}}}\right)+\frac{\partial }{\partial \text{\hspace{0.17em}σ}}\left(\frac{K_{\text{C}}}{D}\frac{\partial C}{\partial \text{σ}}\right)$, (9)

where C is volumetric concentration of suspended particles, m³/m³; $w_{\text{s}}$ is rate of gravitational settling of particles, depending on their density ${\text{ρ}}_{\text{c}}$ and average diameter $d_{50}$; $A_{\text{C}},K_{\text{C}}$ are coefficients of horizontal and vertical diffusion determined in the water circulation model.

At the solid lateral boundaries of the computational domain, conditions for the turbulent fluxes of substance are set $C$ to be equal to zero.

On the free surface, the condition of the suspended flow absence is specified [6, 14]

$-w_{\text{s}}C-\frac{K_{\text{C}}}{D}\frac{\partial C}{\partial \text{σ}}=0$. (10)

The resuspended matter flow from the bottom is determined by the difference between F_e erosion and F_d sedimentation flows [6, 14]:

$-w_{\text{s}}C-\frac{K_{\text{C}}}{D}\frac{\partial C}{\partial \text{σ}}=F_{\text{e}}-F_{\text{d}}$. (11)

F_e and F_d flows are functions of bed shear stress ${\text{τ}}_{\text{b}}$ caused by the combined contribution of currents and waves. In the context of the problem under consideration, the value ${\text{τ}}_{\text{b}}$ plays a key role in modeling the process of bottom sediment resuspension.

The erosion flow is non-zero in the case when bottom shear stresses exceed the critical value t_ce:

$F_{\text{e}}=\left\{\begin{array}{l}{M}_0\left(\frac{{\text{τ}}_{\text{b}}}{{\text{τ}}_{\text{ce}}}-1\right),{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}τ}}_{\text{b}}\ge {\text{τ}}_{\text{ce}},\\ \mathrm{0,}{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}τ}}_{\text{b}}<{\text{τ}}_{\text{ce}},\end{array}\right.$ (12)

where M₀ is empirical coefficient characterizing the intensity of sediment resuspension and varying within the range of 10^–6–10^–2 kg/(m²×s) [15].

The sedimentation flow is determined by the concentration of sediments at the bottom $C_{\text{a}}$ and is different from zero when bottom shear stresses are less than a certain critical value t_cd [6, 14]:

$F_{\text{d}}=\left\{\begin{array}{l}{w}_{\text{s}}{C}_{\text{a}}\left(1-\frac{{\text{τ}}_{\text{b}}}{{\text{τ}}_{\text{cd}}}\right),{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}τ}}_{\text{b}}<{\text{τ}}_{\text{cd}},\\ \mathrm{0,}{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}τ}}_{\text{b}}\ge {\text{τ}}_{\text{cd}},\end{array}\right.$ (13)

where value $C_{\text{a}}$ is assessed by the following formulas [16, p. 673]:

${С}_{\text{a}}=0.015\left(\frac{d_{50}}{a}\right)\frac{{({\text{τ}}_{\text{b}}/{\text{τ}}_{\text{s}}-1)}^{1.5}}{d_{*}^{0.3}}$, $d_{\ast }=d_{50}{\left[\frac{g\text{\hspace{0.17em}Δ}}{{\text{ν}}^2}\right]}^{1/3}$, $\text{Δ}=\frac{{\text{ρ}}_{\text{c}}}{\text{ρ}}-1$, (14)

here $a=0.05D$ is reference level; $d_{\ast }$ is dimensionless diameter of particles; v = 10^-6 m²/s is molecular viscosity coefficient; ${\text{τ}}_{\text{s}}$ is bottom shear stress determined by the Shields criterion.

To evaluate ${\text{τ}}_{\text{b}}$, the following ratios are used ² [7, p. 37]:

${\text{τ}}_{\text{b}}=\sqrt{{({\text{τ}}_{\text{m}}+{\text{τ}}_{\text{w}}\left|\text{cos}\phi \right|)}^2+{\left({\text{τ}}_w\left|\sin \phi \right|\right)}^2}$, ${\text{τ}}_{\text{m}}={\text{τ}}_{\text{c}}\left[1+1.2{\left(\frac{{\text{τ}}_{\text{w}}}{{\text{τ}}_{\text{c}}+{\text{τ}}_{\text{w}}}\right)}^{3.2}\right]$, (15)

where φ is an angle between the direction of current and waves.

Bottom shear stress caused by wave action is defined as

${\text{τ}}_{\text{w}}=\frac{1}{2}\text{ρ\hspace{0.05em}}{f}_{\text{w}}{U}_{\text{w}}^2$, (16)

where $f_{\text{w}}$ is wave drag coefficient; $U_{\text{w}}$ is velocity of bottom orbital movements from the SWAN model. A semi-empirical dependence is used to assess $f_{\text{w}}$ [5, p. 100; 17]

$f_{\text{w}}=\min \left\{\exp \right[5.5{\left(k_{\text{b}}/A\right)}^{0.2}-6.3],0.3\}$, (17)

where $k_{\text{b}}=2.5d_{50}$; $A=U_{\text{w}}{T}_{\text{p}}/2\text{π}$.

Research results and discussion

Dynamics of the waters of Eastern Sivash Bay is largely determined by the local features of atmospheric processes on a synoptic scale. According to ERA-Interim global atmospheric reanalysis data ³ over the period of 1979–2020, statistical characteristics of surface wind speed for the water area under study are as follows[5]: long-term average value is 7.3 m/s; standard deviation is 3.6 m/s. The greatest occurrence frequency (23%) has the northeastern wind and the least occurrence frequency (5%) has the wind of southeastern direction. The frequency occurrence of winds of other directions does not exceed 10–12%.

From the point of view of studying bottom sediment resuspension in Eastern Sivash Bay, the greatest interest is in cases where the water area of the bay is affected by intense non-stationary atmospheric disturbances such as cyclones. Therefore, as an atmospheric forcing from the ERA-Interim reanalysis data, we selected a synoptic situation (from 00:00 November 10, 2007 to 00:00 November 13, 2007) when an anomalously deep cyclone formed in the western part of the Black Sea. The cyclone moved eastward towards the Sea of Azov and during its movement crossed the waters of Eastern Sivash Bay creating significant atmospheric disturbances over it.

Resuspension was modeled with the following values of input parameters: ${\text{ρ}}_{\text{c}}$ = 2000 kg/m³; $d_{50}$ = 0.01 mm; t_ce = 0.13 N/m²; t_cd = 0.1 N/m²; M₀ = 10^–5 kg/m²/s.

These constants are taken from [6, p. 9; 7, p. 40]. Sedimentation rate of particles $w_{\text{s}}$ is determined by predetermined values ${\text{ρ}}_{\text{c}}$ and $d_{50}$. Calculations were carried out for the silt fraction of bottom particles, since the bottom of Eastern Sivash Bay is a silt layer of up to 5 m thickness or more ⁴.

The numerical modeling was carried out in two stages. At the first stage, the wave parameters included in the formulas for calculating bottom stresses (12)–(17) were calculated using the SWAN model. Wave fields were stored with 1 h discreteness. At the second stage, current fields and concentrations of suspended particles were calculated using a water circulation model.

Figure 2, a represents the temporal variation of surface wind speed W over the water area of Eastern Sivash Bay for the synoptic situation under consideration (solid curve – average W value over the water area, dashed curve – maximum W value over the water area). The dependence of average wind direction θ_A over the bay area on time is shown in Fig. 2, b. Time moment t = 0 corresponds to the initial date of synoptic situation in this and other graphs.

The most significant W variations occur in the time interval from 0 to 36 h. Until t = 18 h, wind speed monotonically increases from 4 m/s to its maximum value. At t ≤ 12 h, the prevailing wind is from the south and southeast. After 12 o'clock, wind changes its direction to southwestern. At t > 18 h, the storm fades and the wind speed decreases.

F i g. 2. Wind speed W (a), wind direction θA (b), significant wave heights hs (c) (solid curve shows average values over the water area, dashed line – maximum values over the water area), and integral indicator of bottom sediments resuspension $S_1$ (solid line – with regard to wave stress, dashed line – with no regard to it) (d)

When an atmospheric cyclone passes through the bay, a non-stationary wave field, caused by temporal variability of wind speed and direction, is generated. The dependences of significant wave heights $h_{\text{s}}$ on time for the considered synoptic situation are given in Fig. 2, c (solid curve – average $h_{\text{s}}$ value over the water area, dashed curve – maximum $h_{\text{s}}$ value over the water area). It can be seen that wave height is sensitive to changes in both wind speed and direction. When the wind direction changes, the wave height decreases. A phase shift of 1 h is observed between the time of wind speed maximum occurrence and wave height. The maximum wind speed is reached at t = 18 h and the maximum wave height is reached at t = 19 h.

To quantify the intensity of bottom sediment resuspension in Eastern Sivash Bay, an integral resuspension indicator was applied

$S_n=100\%\cdot A_n/A_0$, (18)

where $A_n$ is bay surface area where the condition for the occurrence of bottom sediment erosion ${\text{τ}}_{\text{b}}\ge n{\text{τ}}_{\text{ce}}$ is fulfilled, $n\ge 1$ is integer; $A_0$ is area of the entire bay.

We performed two calculations of integral resuspension index. In the first calculation of indicator (18), wave stress ${\text{τ}}_{\text{w}}$ was taken into account when determining bottom stresses ${\text{τ}}_{\text{b}}$. $S_1$ dependence on t (solid curve in Fig. 2, d) has two peaks; they correspond to the maximum values: $S_1$ = 50% at t = 15 h; $S_1$ = 80% at t = 20 h. Thus, during the period of the maximum storm development, conditions for resuspension of bottom sediments are created on 80% of the bay water area. $S_1$ curve peaks are shifted relative to the maximum wave heights by 1 h and they are quite sharp, which indicates a rapid course of erosion and deposition processes in bottom sediments. The strongest resuspension occurs during the period of maximum development of wind waves and lasts 3–4 hours. At $n>1$ (zones of intense resuspension), the following values of the maximum resuspension index were obtained: $S_2$ = 35%; $S_3$ = 10%; $S_4$ = 5%.

The second calculation of the resuspension index (dashed curve in Fig. 2, d) was carried out with no regard to the wave stress (${\text{τ}}_{\text{b}}={\text{τ}}_{\text{c}}$). In this case, only one maximum $S_1$ = 21% at t = 18 h arises, which corresponds to the occurrence time of the maximum wind speed. From a comparison of both options for calculating $S_1$ indicator, the following conclusion can be drawn: bottom wave stress makes a decisive contribution to the formation of areas of bottom sediment resuspension, which is explained by the bay shallowness.

Now, let us move on to considering the spatial structure of model fields. As a typical example, Fig. 3 indicates spatial distributions of wind speed W, significant wave height $h_{\text{s}}$, velocity of bottom wave currents $U_{\text{w}}$ and velocity of wind currents in the bottom layer $U_{\text{c}}$.

F i g. 3. Wind speed and direction (a), significant wave heights (b), velocity of the bottom wave (c) and bottom wind (d) currents in Eastern Sivash Bay at t = 19 h

These fields correspond to the time moment t = 19 h, when the wind waves in the bay reached their maximum intensity. The wind has a predominant western-southwestern direction and its speed is maximum in the southern part of the bay (Fig. 3, a). Here, we can also see the area of maximum waves with wave heights exceeding 0.8 m (Fig. 3, b). Wave intensification occurs in the hollows. The areas of maximum waves are shifted to the windward coast. The velocity of bottom wave currents reaches 0.4 m/s. In this case, the configurations of the velocity field of near-bottom wave currents (Fig. 3, c) and the wave height field are well-corresponded to each other. The areas of maximum bottom wind currents with up to 0.8 m/s velocities are localized in shallow water near the leeward shores of the basin and in the narrowness connecting the northern and southern parts of the bay (Fig. 3, d).

F i g. 4. Spatial distribution of bottom shear stresses in Eastern Sivash Bay and wind current velocities in the bottom layer at t = 19 h

Figure 4 demonstrates the fields of bottom shear stresses ${\text{τ}}_{\text{с}}$, ${\text{τ}}_{\text{w}}$ and ${\text{τ}}_{\text{b}}$ for t = 19 h normalized to the critical value t_ce. In general, the configuration of these fields is similar to the one of the bottom current fields (Fig. 4, d). The component of shear stresses caused by wind currents increases near the leeward shores and the wave component of shear stresses, on the contrary, increases near the windward shores of the basin.

In some cases, the values ${\text{τ}}_{\text{с}}$ are noticeably higher than ${\text{τ}}_{\text{w}}$. However, strong bottom stresses caused by currents are concentrated in small areas, while waves generate bottom stresses over ${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$ a much larger area of the bay. Thus, when calculating the total bottom stress ${\text{τ}}_{\text{b}}$, it is necessary to take into account both the bottom wave-generated stress and the stress generated by currents.

Analysis of C vertical structure of the concentration field showed that the concentration profiles of silty suspended matter depend weakly on depth during a storm. The resuspension of silty sediments occurs very quickly. The maximum differences between C values at the surface and at the bottom are no more than 0.1 mg/l. This consistency of profiles is due to bay shallowness and a vertical mixing effect.

In Fig. 5, the fields of depth-average suspended matter concentration for characteristic moments of time are shown:

$C_{\text{m}}=\underset{-1}{\overset{0}{\int }}Cd\text{σ}$.

Local turbidity sources, which increase in size and intensify over time, appear near the coast at the initial stage (t = 8 h). By t = 20 h time point, the bottom sediment resuspension reaches its maximum intensity. Further, sedimentation processes begin – suspended matter concentration in the water column decreases quite quickly and the turbidity areas begin to shrink (t = 24 h).

During a storm, in areas where the bottom wind current velocities are maximum, the calculated values $C_{\text{m}}$ at individual points reach 250 mg/l. In numerical experiments, this result is determined primarily by the bay shallowness, the silty nature of bottom and low values of particle sedimentation rate. The question of whether such large concentration values can exist in Eastern Sivash remains open, since no data on suspended matter concentrations under storm conditions is currently available. It can be noted that calculations for the completely silty bottom of Neva Bay provide the same order of magnitude values for the suspended matter concentration [7].

In the bottom sediment resuspension model (9)–(17), as well as in other similar models, the values of input parameters have a large degree of uncertainty. Therefore, it is of interest to evaluate the sensitivity of the resuspension model under consideration to changes in the values of these parameters.

F i g. 5. Field of the depth-average concentration of suspended matter ${С}_m$ in Eastern Sivash Bay for different time points

For quantitative assessment of resuspension model sensitivity to variations in input parameters, 12 numerical experiments with different combinations of parameters were performed. The list of parameter values is given in Table 1. Model solution deviations from the base case (k = 1) were analyzed.

T a b l e 1

Parameters of bottom sediment resuspension model

k	${\text{ρ}}_{\text{c}}$, kg/m3	$d_{50}$, mm	tce, Н/m2	tcd, Н/m2	M0, kg/m2/s
1	2000	0.010	0.130	0.10	10−5
2	2000	0.010	0.156	0.10	10−5
3	2000	0.010	0.104	0.10	10−5
4	2000	0.015	0.130	0.10	10−5
5	2000	0.020	0.130	0.10	10−5
6	2000	0.010	0.130	0.10	10−4
7	2000	0.010	0.130	0.10	10−6
8	2000	0.010	0.130	0.10	10−3
9	2000	0.010	0.130	0.12	10−5
10	2000	0.010	0.130	0.08	10−5
11	2100	0.010	0.130	0.10	10−5
12	1900	0.010	0.130	0.10	10−5

N o t e: k is the experiment number.

For each k we calculated integral indicator

$B_{kp}=100\%\cdot Q_{kp}/Q_0$,  (19)

where index p is threshold value of average suspended matter concentration $C_{\text{m}}$, characterizing resuspension intensity; $Q_{kp}$ is total surface area of the bay where the condition $C_{\text{m}}\ge p$ is satisfied; $Q_0$ is total surface area of the entire bay.

Indicator (19) provides the total relative size of bottom sediment resuspension sources of different intensity at a specific time point. Five threshold p values (5, 10, 25, 50 and 100 mg/l) were considered. It was revealed that the greatest increase in B_kp index occurs during the period of maximum storm development (t = 19–22 h); therefore, B_kp maxima sensitivity to changes in input parameters was analyzed.

Values of $\max B_{kp}$, as well as the difference ${\text{γ}}_{kp}=\max B_{kp}-\max B_{1p}$ (Table 2), which is a deviation of B_kp maximum from $B_{1p}$ maximum of the base experiment, are given below for each numerical experiment.

T a b l e 2

Values of $\max B_{kp}$  (%) and ${\gamma }_{kp}$  (%) for different options of specifying the resuspension model parameters

k	p, mg/l
	5		10		25		50		100
	$\max B_{kp}$	${\text{γ}}_{kp}$	$\max B_{kp}$	${\text{γ}}_{kp}$	$\max B_{kp}$	${\text{γ}}_{kp}$	$\max B_{kp}$	${\text{γ}}_{kp}$	$\max B_{kp}$	${\text{γ}}_{kp}$
1	94	0	91	0	82	0	70	0	49	0
2	91	−3	86	−5	73	−9	58	−12	36	−13
3	97	3	95	4	88	6	80	10	65	16
4	90	−4	86	−5	76	−6	62	−8	37	−12
5	87	−7	82	−9	70	−12	53	−17	27	−22
6	96	2	94	3	90	8	87	17	83	34
7	76	−18	57	−34	24	−58	9	−61	2	−47
8	92	−2	90	−1	85	3	82	12	78	29
9	94	0	91	0	81	1	68	2	47	2
10	95	0	92	1	83	0	72	0	51	0
11	94	0	91	1	81	0	70	0	49	0
12	94	0	92	1	82	0	70	0	49	0

In all experiments, p increases with B_kp index decrease. The behavior of ${\text{γ}}_{kp}$ parameter demonstrates the opposite trend: with increasing p, the deviations of maxima from the base experiment in absolute value increase, i.e. the impact of variations in the model parameters is more pronounced for areas with large suspended matter concentrations.

In experiments 2 and 3, the influence of critical shear stress t_ce on variations in resuspension areas was assessed. The t_ce value changed by ±20% relative to the base value. With an increase in t_ce by 20%, depending on p values, the resuspension areas decrease by 3–13%. When t_ce decreases by 20%, the resuspension areas are reduced by 3–16%.

In experiments 4 and 5, the average particle diameter $d_{50}$ varied. The increase of $d_{50}$ by 50 and 100% results in a decrease in resuspension areas by 4–12 and 7–22%, respectively.

In experiments 6, 7 and 8, the value of M₀ parameter varied. When M₀ decreases by an order of magnitude, the resuspension areas are reduced by 18–61%. An increase in this parameter by two orders of magnitude increases the resuspension areas by 2–34%.

In experiments 9 and 10, t_cd stress effect on variations in resuspension areas was assessed. A weak dependence of resuspension areas on the values of this parameter was revealed (variations do not exceed 2%).

As shown by experiments 11 and 12, variations in the density of bottom sediment particles ${\text{ρ}}_{\text{c}}$ by ±5% relative to the base value also do not have a noticeable effect on variations in resuspension areas (variations do not exceed 2%).

Thus, the most significant parameters in the bottom sediment resuspension model we applied, are the parameters t_ce and M₀, which determine the intensity of vertical flow of particles from the bottom.

It seems possible to use an integral indicator of B_kp type to clarify the values of the resuspension model input parameters based on remote sensing data for certain areas of the Sea of Azov. To do this, it is necessary to obtain the resuspension areas from satellite images using a regional algorithm for estimating the suspended matter concentration [18] and calculate the B_kp index. Next, based on a series of model calculations, we are to select the value of this indicator, which is close to the indicator value obtained from remote sensing data. The problem of estimating the parameters of a model for bottom sediment resuspension in shallow waters is of independent interest and the authors considered it as a continuation of the research presented in this paper.

Conclusion

Based on a combination of numerical models of wind waves and currents, the process of bottom silty sediment resuspension in the eastern part of Sivash Bay during the extreme storm on November 10–13, 2007 was studied. Analysis of numerical modeling results revealed the following.

The strongest bottom sediment resuspension in the bay occurs during the period of wind wave maximum intensification and lasts 3–4 h. During the storm, vertical profiles of silty suspended matter C concentration weakly depend on depth. The maximum differences between C values at the surface and at the bottom are no more than 0.1 mg/l.

Bottom wave stresses ${\text{τ}}_{\text{w}}$ make a decisive contribution to the formation of bottom sediment resuspension areas in the bay. When they are taken into account, the conditions for forming resuspension areas (${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$) are created on 80% of the bay area. If wave stresses are not taken into account during modeling, the total resuspension area is reduced by ~ 4 times.

Currents can also generate bottom stresses ${\text{τ}}_{\text{c}}\ge {\text{τ}}_{\text{ce}}$, but this occurs in small parts of the water area near the shore, while waves generate bottom stresses ${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$ throughout most of the bay. Thus, when calculating the total bottom stress ${\text{τ}}_{\text{b}}$, it is necessary to take into account both wave stresses and stresses caused directly by currents.

Based on multivariate calculations, it was found that the silt resuspension model we applied is most sensitive to variations in the values of t_ce and M₀ parameters, which determine intensity of the vertical flow of particles from the basin bottom.

The results of the work can be useful in planning and interpreting field experiments in Sivash Bay and other shallow areas of the Sea of Azov.

¹ USER MANUAL SWAN Cycle III Version 41.45. Delft University of Technology. [online] Available at: https://swanmodel.sourceforge.io/online_doc/swanuse/swanuse.html [Accessed: 03 June 2024].

² Soulsby, R.L., 1997. Dynamics of Marine Sands: A Manual for Practical Applications. London: Tomas Telford Services, 249 p.

³ ECMWF. Archive Catalogue. [online] Available at: apps.ecmwf.int [Accessed: 03 June 2024].

⁴ Stashchuk, M.F., Suprychev, V.A. and Khitraya, M.S., 1964. [Mineralogy, Geochemistry and Conditions for the Formation of Bottom Sediments of the Sivash]. Kiev: Naukova Dumka, 174 p. (in Russian).