Взмучивание донных осадков в мелководной лагуне течениями и волнами по данным численного моделирования (на примере залива Сиваш, Азовское море)

Полный текст

Введение

Необходимость изучения динамики вод залива Сиваш Азовского моря связана с активными антропогенными преобразованиями данной акватории прибрежной зоны Крыма. Залив является восстанавливаемым месторождением минеральных солей, отличается высоким биологическим разнообразием и входит в перечень водно-болотных угодий международного значения. На побережье залива функционирует специальная экономическая зона «Сиваш» с развитой химической промышленностью и сельским хозяйством [1].

Современный Сиваш представляет собой обширный мелководный залив Азовского моря (морской залив лагунного типа) с изрезанной береговой линией. Сиваш принято разделять на два больших района – Западный и Восточный Сиваш (рис. 1, а). Западный Сиваш в настоящее время является почти полностью изолированным водоемом, уровень воды в котором регулируется дамбой, его глубина обычно не превышает 0,3–0,4 м. Восточный Сиваш – более глубоководный (глубины до 3 м) участок залива, занимающий свыше 60 % общей площади Сиваша.

 

Рис. 1. Расположение залива Восточный Сиваш (а), модельный рельеф дна залива Восточный Сиваш (b)

Fig. 1. Location of the East Sivash Bay (a), model bottom relief of the East Sivash Bay (b)

 

Изучению гидродинамических процессов в Восточном Сиваше посвящено несколько публикаций. В работе [2] приведены расчеты поверхностных течений и колебаний уровня в заливе при ветрах различных направлений. В работе [3] изучены характеристики пресноводного плюма в районе впадения реки Салгир в залив Сиваш. В статье [4] моделируются поля течений в акватории Восточного Сиваша для разных ветровых условий. Показано, что структура циркуляции вод в заливе определяется генеральным направлением ветра. Поле средних по глубине течений состоит из нескольких вихревых структур. Особенно хорошо эти структуры выражены в южной части акватории. При ветрах, направленных вдоль акватории залива, циркуляция наиболее интенсивна. В этих случаях в узости, соединяющей северную и южную части Восточного Сиваша, возникает струйное течение со скоростью более 1 м/с. На самых глубоких участках акватории у дна формируется придонное противотечение.

В [5] выполнено моделирование характеристик ветровых волн в Восточном Сиваше. Показано, что вследствие мелководности акватории и сравнительно небольших эффективных разгонов основным фактором, влияющим на интенсивность волнения, является скорость ветра. Волнение выходит на установившийся режим через 3 ч после начала действия ветра. После прекращения действия ветра волнение полностью затухает через 2 ч. Оценки максимальных значений параметров волнения для разных градаций скорости ветра показали, что даже в период сильных штормов высоты и периоды волн в акватории Восточного Сиваша не превышают 1 м и 3 с соответственно. В указанных работах за пределами проведенных исследований остались проблемы воздействия течений и волн на динамику донных осадков в акватории Восточного Сиваша.

Задачи математического моделирования эволюции взвеси в мелководных акваториях рассмотрены в [6, 7]. В работе [6] моделируются процессы взмучивания и осаждения донных осадков на северо-западном шельфе Черного моря при прохождении атмосферного циклона без учета ветрового волнения. Донные отложения предполагаются состоящими из одной илистой фракции. Анализируются области наиболее интенсивного подъема взвеси и вертикальные профили концентрации взвешенных донных осадков в отдельных районах. Подобное исследование для восточной части Финского залива, но с учетом ветрового волнения, проведено в работе [7]. Здесь расчеты были выполнены для двух фракций донных осадков – мелкозернистого песка и ила. Исследуется вклад течений и волнения в интенсивность взмучивания осадков. Показана важность учета вклада ветрового волнения в процессы взмучивания.

Цель настоящей работы – математическое моделирование процесса взмучивания донных осадков в Восточном Сиваше. Математическое моделирование является важной частью экологического мониторинга акватории, с его помощью можно осуществить расчет различных параметров состояния залива, которые не всегда удается получить в результате натурных измерений. В статье приведены результаты численных экспериментов по моделированию взмучивания илистых осадков в период экстремального шторма. Ниже дается характеристика используемых математических моделей. Основное внимание уделяется анализу вклада течений и волнения в интенсивность взмучивания.

Основной проблемой математического моделирования взмучивания донных осадков является отсутствие необходимого количества исходных данных. Как правило, гранулометрический состав донных отложений известен фрагментарно, нет достаточного количества данных натурных наблюдений для валидации моделей, выбор критических значений придонных сдвиговых напряжений, определяющих процессы взмучивания, представляет собой отдельную задачу. Имеющиеся в настоящее время натурные исследования в Восточном Сиваше [8, 9] не дают полного представления о динамике взвешенного вещества в заливе, поскольку они эпизодичны, выполнены непосредственно у берега в относительно спокойную погоду. По этой причине был проведен анализ чувствительности модели взмучивания донных осадков к изменению входных параметров.

Материалы и методы

Волновая модель. Для расчета полей ветрового волнения в заливе Восточный Сиваш используется модель Simulating Waves Nearshore (SWAN) [10], основанная на численном решении уравнения баланса волновой энергии в спектральной форме

$\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)

где $N=E/\omega$ – плотность волнового действия; E – энергетический спектр волнения; x, y, t – пространственные координаты и время; $\omega ,\theta$ – частотная и угловая координаты; $c_X,c_Y,c_{\omega },c_{\theta }$ – скорости переноса плотности волнового действия вдоль пространственных и частотно-угловых координат.

Функция источника в правой части уравнения (1) имеет вид

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

где $S_{\text{in}}$ – источник генерации волн ветром; $S_{\text{nl}}$ – нелинейные взаимодействия спектральных гармоник; $S_{\text{wc}}$ – диссипация энергии вследствие обрушения волновых гребней; $S_{\text{bf}}$ – диссипация энергии, обусловленная донным трением; $S_{\text{dib}}$ – обрушение волн на критических глубинах. Детали параметризации отдельных слагаемых функции источника (2) даны в руководстве ¹.

Модельный рельеф дна залива Восточный Сиваш показан на рис. 1, b. Для расчета волновых характеристик используется прямоугольная сетка с горизонтальным разрешением Δx = Δy = 300 м (284 × 334 узла). Дискретность модели по угловой координате составляет 10°. По частотной координате, изменяющейся в диапазоне 0,04–1,0 Гц, применяется неравномерная сетка с 31 узлом. Используется нестационарный вариант модели SWAN. Интегрирование по времени проводится по неявной разностной схеме с шагом 30 мин.

Выходными данными модели SWAN являются поля высоты значительных волн $h_{\text{s}}$, среднего направления волнения ${\text{θ}}_{\text{w}}$, пикового периода волн $T_{\text{p}}$ и амплитуды скорости придонных орбитальных движений $U_{\text{w}}$. Эти параметры используются для расчета придонных сдвиговых напряжений в модели течений и модели взмучивания донных осадков.

Модель течений. Для расчета полей течений используется баротропный вариант трехмерной гидростатической модели циркуляции вод, предложенный в работе [11]. Модель основана на трехмерных σ-координатных уравнениях гидродинамики в гидростатическом приближении (здесь и далее проводится суммирование по повторяющимся индексам α и β от 1 до 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)

где $(x_1,x_2)$ = $(x,y)$; σ – безразмерная вертикальная координата, изменяющаяся в пределах от –1 до 0; (u₁, u₂) = (u, v) – компоненты скорости течения по осям ;– компонента скорости течения, направленная по нормали к поверхностям σ = const; $D=h+\eta$; $h$ – глубина бассейна; $\eta$ – координата свободной поверхности; f – параметр Кориолиса; ${\text{ε}}_{\text{αβ}}=0$ при $\text{α}=\text{β}$ и ${\text{ε}}_{12}=-1$, ${\text{ε}}_{21}=1$; ${\text{τ}}_{\text{αβ}}$ – компоненты тензора турбулентных напряжений; $A_{\text{M}},K_{\text{M}}$ – коэффициенты турбулентной вязкости; g – ускорение силы тяжести.

На твердых боковых границах расчетной области ставятся условия прилипания для скоростей и условия равенства нулю турбулентных потоков. Для упрощения задачи водообмен через пролив Тонкий не учитывается.

На свободной поверхности бассейна ($\text{σ}=0$) граничные условия имеют вид

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

где ${\text{ρ}}_0$ – плотность воздуха; $c_0=(\mathrm{0,49}+\mathrm{0,065}W)\cdot {10}^{-3}$ – коэффициент поверхностного трения; $W=\sqrt{W_1^2+W_2^2}$ – модуль скорости приводного ветра; $W_1, W_2$ – зональная и меридиональная компоненты скорости ветра.

У дна ($\text{σ}=-1$) граничные условия задаются следующими выражениями:

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

где $U_{\text{c}}=\sqrt{u_1^2+u_2^2}$ – модуль горизонтальной скорости течений; ${\text{τ}}_{\text{с}}$ – придонное сдвиговое напряжение, определяемое формулой

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

где $c_{\text{b}}=\mathrm{0,16}/{\text{ln}}^{\text{2}}\text{(δ\hspace{0.05em}}z\text{/}{z}_{\text{b}}\text{)}$ – коэффициент донного трения; $\text{δ\hspace{0.05em}}z$ – расстояние по вертикали от дна до точки, в которой определяется коэффициент трения; z_b = d₅₀/30 – параметр шероховатости донной поверхности; d₅₀ – средний диаметр частиц донного грунта, соответствующий иловым отложениям.

Для определения коэффициента вертикального турбулентного обмена К_М используется модель Меллора – Ямады [12]. Коэффициент горизонтального турбулентного обмена А_М вычисляется по формуле Смагоринского [13].

Краевая задача (3)–(8) решается численно на основе явных разностных схем по горизонтальным координатам и неявных разностных схем по вертикальной координате. Для аппроксимации оператора адвекции $\Lambda$ применяются монотонные разностные схемы. Детальное описание численного алгоритма дано в [11]. По горизонтальным координатам используется та же расчетная сетка, что и в модели SWAN. По σ-координате задается 11 равномерных расчетных уровней.

Модель взмучивания донных осадков. Модель взмучивания донных осадков включена в модель циркуляции вод и основана на численном решении уравнения переноса – диффузии вида

$\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)

где С – объемная концентрация взвешенных частиц, м³/м³; $w_{\text{s}}$ – скорость гравитационного оседания частиц, зависящая от их плотности ${\text{ρ}}_{\text{c}}$ и среднего диаметра d₅₀; $A_{\text{C}},K_{\text{C}}$ – коэффициенты горизонтальной и вертикальной диффузии, определяемые в модели циркуляции вод.

На твердых боковых границах расчетной области ставятся условия равенства нулю турбулентных потоков субстанции C.

На свободной поверхности задается условие отсутствия потока взвеси [6, 14]

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

Поток взмученного вещества со дна определяется разностью потоков эрозии F_e и седиментации F_d [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 и F_d являются функциями придонного сдвигового напряжения ${\text{τ}}_{\text{b}}$, обусловленного совместным вкладом течений и волн. В контексте рассматриваемой задачи величина ${\text{τ}}_{\text{b}}$ играет ключевую роль в моделировании процесса взмучивания донных осадков.

Поток эрозии отличен от нуля в случае, когда придонные сдвиговые напряжения превышают критическое значение τ_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)

где M₀ – эмпирический коэффициент, характеризующий интенсивность взмучивания осадков и изменяющийся в пределах 10^–6–10^–2 кг/(м²×с) [15].

Поток седиментации определяется концентрацией осадков у дна C_a и отличен от нуля, когда придонные сдвиговые напряжения меньше некоторого критического значения τ_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)

где величина C_a оценивается по формулам [16, c. 673]

${С}_{\text{a}}=\mathrm{0,015}\left(\frac{d_{50}}{a}\right)\frac{{({\text{τ}}_{\text{b}}/{\text{τ}}_{\text{s}}-1)}^{\mathrm{1,5}}}{d_{*}^{\mathrm{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)

здесь $a=\mathrm{0,05}D$ – отсчетный уровень; $d_{\ast }$ – безразмерный диаметр частиц; $\text{ν}$ = 10^-6 м²/с – коэффициент молекулярной вязкости; ${\text{τ}}_{\text{s}}$ – придонное сдвиговое напряжение, определяемое по критерию Шильдса.

Для оценки ${\text{τ}}_{\text{b}}$ используются соотношения ² [7, c. 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+\mathrm{1,2}{\left(\frac{{\text{τ}}_{\text{w}}}{{\text{τ}}_{\text{c}}+{\text{τ}}_{\text{w}}}\right)}^{\mathrm{3,2}}\right]$, (15)

где φ – угол между направлением течения и волн.

Придонное сдвиговое напряжение, обусловленное действием волн, определяется в виде

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

где $f_{\text{w}}$ – коэффициент волнового сопротивления; $U_{\text{w}}$ – скорость придонных орбитальных движений из модели SWAN. Для оценки $f_{\text{w}}$ применяется полуэмпирическая зависимость [5, с. 100; 17]

$f_{\text{w}}=\min \left\{\exp \right[\mathrm{5,5}{\left(k_{\text{b}}/A\right)}^{\mathrm{0,2}}-\mathrm{6,3}],\mathrm{0,3}\}$, (17)

где $k_{\text{b}}=\mathrm{2,5}{d}_{50}$; $A=U_{\text{w}}{T}_{\text{p}}/2\text{π}$.

Результаты исследования и их обсуждение

Динамика вод залива Восточный Сиваш в значительной мере определяется локальными особенностями атмосферных процессов синоптического масштаба. Согласно данным глобального атмосферного реанализа ERA-Interim ³ за период 1979–2020 гг., статистические характеристики скорости приземного ветра  для исследуемой акватории следующие [5]: среднее многолетнее значение – 7,3 м/c; среднеквадратическое отклонение – 3,6 м/c. Наибольшую повторяемость (23 %) имеет ветер северо-восточного направления, а наименьшую (5 %) – ветер юго-восточного направления. Повторяемость ветров других направлений не превышает 10–12 %.

С точки зрения исследований взмучивания донных осадков в заливе Восточный Сиваш наибольший интерес представляют случаи воздействия на акваторию залива интенсивных нестационарных атмосферных возмущений типа циклонов. Поэтому в качестве атмосферного форсинга из данных реанализа ERA-Interim была выбрана синоптическая ситуация (с 00:00 10 ноября 2007 г. до 00:00 13 ноября 2007 г.), когда в западной части Черного моря сформировался аномально глубокий циклон. Циклон смещался в восточном направлении в сторону Азовского моря и при своем движении пересек акваторию залива Восточный Сиваш, создав над ним существенные атмосферные возмущения.

Взмучивание моделировалось при следующих значениях входных параметров: ${\text{ρ}}_{\text{c}}$ = 2000 кг/м³; d₅₀ = 0,01 мм; τ_ce = 0,13 Н/м²; τ_cd = 0,1 Н/м²; M₀ = 10^–5 кг/м²/с.

Эти константы взяты из работ [6, c. 9; 7, c. 40]. Скорость оседания частиц  определяется по заданным значениям ${\text{ρ}}_{\text{c}}$ и d₅₀. Расчеты проводились для илистой фракции донных частиц, поскольку дно залива Восточный Сиваш представляет собой слой ила толщиной до 5 м и более ⁴.

Численное моделирование выполнялось в два этапа. На первом этапе с помощью модели SWAN рассчитывались волновые параметры, входящие в формулы расчета придонных напряжений (12)–(17). Волновые поля сохранялись с дискретностью 1 ч. На втором этапе с использованием модели циркуляции вод рассчитывались поля течений и концентрации взвешенных частиц.

На рис. 2, а показано изменение во времени скорости приводного ветра  над акваторией залива Восточный Сиваш для рассматриваемой синоптической ситуации (сплошная кривая – среднее по акватории значение , штриховая кривая – максимальное по акватории значение ). Зависимость от времени среднего направления ветра θ_A над акваторией залива дана на рис. 2, b. На этом и других графиках момент времени t = 0 соответствует начальной дате синоптической ситуации.

Наиболее значительные изменения  происходят на интервале времени от 0 до 36 ч. До t = 18 ч скорость ветра монотонно увеличивается от 4 м/с до своего максимального значения. При t ≤ 12 ч преобладает ветер южного и юго-восточного направлений. После 12 ч ветер меняет свое направление на юго-западное. При t > 18 ч происходит затухание шторма и уменьшение скорости ветра.

 

Рис. 2. Скорость ветра W (а), направление ветра θ_A (b), высоты значительных волн h_s (c) (сплошная кривая – средние по акватории значения, штриховая – максимальные по акватории значения) и интегральный показатель взмучивания донных осадков  (сплошная линия – с учетом волнового напряжения, штриховая – без учета) (d)

Fig. 2. Wind speed W (a), wind direction θ_A (b), significant wave heights h_s (c) (solid curve shows the average values over the water area, dashed line – the maximum values over the water area), and integral indicator of the bottom sediments resuspension  (solid line – with regard for a wave stress, dashed line – with no regard for it) (d)

 

При прохождении атмосферного циклона в заливе генерируется нестационарное поле волнения, что обусловлено изменчивостью во времени скорости и направления ветра. Зависимости от времени высоты значительных волн $h_{\text{s}}$ для рассматриваемой синоптической ситуации приведены на рис. 2, c (сплошная кривая – среднее по акватории значение $h_{\text{s}}$, штриховая кривая – максимальное по акватории значение $h_{\text{s}}$). Видно, что высота волн чувствительна к изменениям как скорости ветра, так и его направления. При смене направления ветра высота волн уменьшается. Прослеживается сдвиг по фазе в 1 ч между временем возникновения максимумов скорости ветра и высоты волн. Максимум скорости ветра достигается при t = 18 ч, а максимум высоты волн – при t = 19 ч.

Для количественной оценки интенсивности взмучивания донных осадков в заливе Восточный Сиваш использовался интегральный показатель взмучивания

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

где $A_n$ – площадь поверхности залива, для которой выполняется условие возникновения эрозии донных осадков ${\text{τ}}_{\text{b}}\ge n{\text{τ}}_{\text{ce}}$, $n\ge 1$ – целое число; $A_0$ – площадь поверхности всего залива.

Было выполнено два расчета интегрального показателя взмучивания. В первом расчете показателя (18) при определении придонных напряжений ${\text{τ}}_{\text{b}}$ учитывалось волновое напряжение ${\text{τ}}_{\text{w}}$. Зависимость $S_1$ от t (сплошная кривая на рис. 2, d) имеет два пика, им соответствуют максимальные значения:  $S_1$= 50 % при t = 15 ч; $S_1$ = 80 % при t = 20 ч. Таким образом, в период максимального развития шторма на 80 % акватории залива создаются условия для взмучивания донных осадков. Пики кривой смещены относительно максимумов высоты волн на 1 ч вперед и они достаточно острые, что указывает на быстрое протекание процессов эрозии и осаждения донных осадков. Наиболее сильное взмучивание возникает в период максимального развития ветрового волнения и длится в течение 3–4 ч. При $n>1$ (зоны интенсивного взмучивания) получены следующие значения максимумов показателя взмучивания: $S_2$ = 35 %; $S_3$ = 10 %; $S_4$ = 5 %.

Второй расчет показателя взмучивания (штриховая кривая на рис. 2, d) выполнялся без учета волнового напряжения (${\text{τ}}_{\text{b}}={\text{τ}}_{\text{c}}$). В этом случае проявляется только один максимум $S_1$ = 21 % при t = 18 ч, что соответствует времени наступления максимума скорости ветра. Из сопоставления обоих вариантов расчетов показателя $S_1$ следует вывод: придонное волновое напряжение вносит определяющий вклад в формирование областей взмучивания донных осадков, что объясняется мелководностью залива.

Перейдем к рассмотрению пространственной структуры модельных полей. В качестве характерного примера на рис. 3 показаны пространственные распределения скорости ветра W, высоты значительных волн $h_{\text{s}}$, скорости придонных волновых течений $U_{\text{w}}$ и скорости ветровых течений в придонном слое $U_c$.

 

Рис. 3. Скорость и направление ветра (а), высота значительных волн (b), скорость придонных волновых (c) и придонных ветровых (d) течений в заливе Восточный Сиваш при t = 19 ч

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

 

Эти поля соответствуют моменту времени t = 19 ч, когда ветровое волнение в заливе достигло максимальной интенсивности. Ветер имеет преобладающее направление западо-юго-запад и его скорость максимальна в южной части залива (рис. 3, а). Здесь же прослеживается область максимального волнения с высотами волн, превышающими 0,8 м (рис. 3, b). В котловинах происходит интенсификация волнения. Области максимумов волнения смещены к наветренному берегу. Скорости придонных волновых течений достигают 0,4 м/с. При этом конфигурации поля скоростей придонных волновых течений (рис. 3, c) и поля высот волн хорошо соответствуют друг другу. Области максимумов придонных ветровых течений со скоростями до 0,8 м/с локализованы на мелководье у подветренных берегов бассейна и в узости, соединяющей северную и южную части залива (рис. 3, d).

 

Рис. 4. Пространственное распределение придонных сдвиговых напряжений в заливе Восточный Сиваш и скорости ветровых течений в придонном слое при t = 19 ч

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

 

На рис. 4 для t = 19 ч приведены поля придонных сдвиговых напряжений ${\text{τ}}_{\text{с}}$, ${\text{τ}}_{\text{w}}$ и  ${\text{τ}}_{\text{b}}$, нормированные на критическое значение τ_ce. В целом конфигурация этих полей подобна конфигурации полей придонных течений (рис. 4, d). Компонента сдвиговых напряжений, обусловленная ветровыми течениями, усиливается у подветренных берегов, а волновая компонента сдвиговых напряжений, напротив, усиливается у наветренных берегов бассейна.

В отдельных случаях значения ${\text{τ}}_c$ заметно превосходят значения ${\text{τ}}_w$. Тем не менее сильные придонные напряжения, вызванные течениями, сосредоточены в небольших по площади районах, в то время как волны генерируют придонные напряжения ${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$ на гораздо большей площади залива. Таким образом, при расчете суммарного придонного напряжения ${\text{τ}}_{\text{b}}$ необходимо учитывать как придонное напряжение, создаваемое волнами, так и напряжение, генерируемое течениями.

Анализ вертикальной структуры поля концентрации C показал, что профили концентрации илистой взвеси во время шторма слабо зависят от глубины. Взмучивание илистых осадков происходит очень быстро. Максимальные перепады между значениями C на поверхности и у дна составляют не более 0,1 мг/л. Такое постоянство профилей обусловлено мелководностью залива и действием вертикального перемешивания.

На рис. 5 для характерных моментов времени показаны поля средней по глубине концентрации взвешенного вещества

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

На начальной стадии (t = 8 ч) у берегов возникают локальные очаги взмучивания, которые со временем увеличиваются в размерах и усиливаются. К моменту времени t = 20 ч взмучивание донных осадков достигает максимальной интенсивности. Далее начинаются процессы осаждения – концентрация взвеси в толще воды достаточно быстро уменьшается, а области взмучивания начинают сокращаться (t = 24 ч).

В период шторма в зонах, где максимальны скорости придонных ветровых течений, расчетные значения $C_{\text{m}}$ в отдельных точках достигают 250 мг/л. В численных экспериментах такой результат определяется в первую очередь малыми глубинами залива, илистым характером дна и малыми значениями скорости оседания частиц. Вопрос о том, могут ли быть в Восточном Сиваше такие большие значения концентрации, остается открытым, поскольку в настоящее время нет данных о концентрациях взвеси в штормовых условиях. Можно отметить, что расчеты для полностью илистого дна Невской губы дают такие же по порядку значения концентрации взвеси [7].

В модели взмучивания донных осадков (9)–(17), как, впрочем, и в других подобных моделях, значения входных параметров имеют большую степень неопределенности. Поэтому представляет интерес оценить чувствительность рассматриваемой модели взмучивания к изменениям значений этих параметров.

 

Рис. 5. Поле средней по глубине концентрации взвешенного вещества  в заливе Восточный Сиваш для разных моментов времени

Fig. 5. Field of the depth-average concentration of suspended matter  in the Eastern Sivash Bay for different time points

 

Для количественной оценки чувствительности модели взмучивания к вариациям входных параметров было выполнено 12 численных экспериментов с разными комбинациями параметров. Список значений параметров приведен в табл. 1. Анализировались отклонения модельного решения от базового варианта (k = 1).

 

Таблица 1

Table 1

Параметры модели взмучивания донных осадков

Parameters of the bottom sediments resuspension model

k	${\text{ρ}}_{\text{c}}$, кг/м3 / ${\text{ρ}}_{\text{c}}$, kg/m3	d50, мм / d50, mm	τce, Н/м2 / τce, Н/m2	τcd, Н/м2 / τcd, Н/m2	M0, кг/м2/с / 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

Примечание: k – номер эксперимента.

Note: k is the experiment number.

 

Для каждого k рассчитывался интегральный показатель

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

где индекс p – пороговое значение средней концентрации взвешенного вещества $C_{\text{m}}$, характеризующее интенсивность взмучивания; $Q_{kp}$ – общая площадь поверхности залива для которой выполняется условие $C_{\text{m}}\ge p$; $Q_0$ – общая площадь поверхности всего залива.

Показатель (19) дает суммарный относительный размер областей взмучивания донных осадков разной интенсивности в конкретный момент времени. Рассматривалось пять пороговых значений p (5, 10, 25, 50 и 100 мг/л). Было выявлено, что наибольшее увеличение показателя B_kp происходит в период максимального развития шторма (t = 19–22 ч), поэтому анализировалась чувствительность максимумов B_kp к изменениям входных параметров.

Ниже для каждого численного эксперимента дан $\max B_{kp}$, а также разность ${\text{γ}}_{kp}=\max B_{kp}-\max B_{1p}$ (табл. 2), представляющая собой отклонение максимума B_kp от максимума $B_{1p}$ базового эксперимента.

 

Таблица 2

Table 2

Значения  (%) и  (%) для разных вариантов задания параметров модели взмучивания

Values of  (%) and  (%) for different options of specifying the resuspension model parameters

k	p, мг/л / 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

 

Во всех экспериментах с увеличением p происходит уменьшение показателя B_kp. В поведении параметра ${\text{γ}}_{kp}$ прослеживается обратная тенденция – с увеличением p отклонения максимумов от базового эксперимента по абсолютной величине увеличиваются, т. е. влияние вариаций параметров модели сильнее сказывается для областей с большими значениями концентрации взвеси.

В экспериментах 2 и 3 оценивалось влияние критического сдвигового напряжения τ_ce на изменения площадей взмучивания. Значение τ_ce менялось на ±20 % относительно базового значения. При увеличении τ_ce на 20 % происходит, в зависимости от значений p, сокращение площадей взмучивания на 3–13 %. При уменьшении τ_ce на 20 % площади взмучивания увеличиваются на 3–16 %.

В экспериментах 4 и 5 варьировался средний диаметр частиц d₅₀. Увеличение d₅₀ на 50 и 100 % дает уменьшение площадей взмучивания на 4–12 и 7–22 % соответственно.

В экспериментах 6, 7 и 8 менялось значение параметра M₀. При уменьшении M₀ на порядок площади взмучивания сокращаются на 18–61 %. Увеличение этого параметра на два порядка дает увеличение площадей взмучивания на 2–34 %.

В экспериментах 9 и 10 оценивалось влияние напряжения τ_cd на изменения площадей взмучивания. Выявлена слабая зависимость площадей взмучивания от значений этого параметра (изменения не превышают 2 %).

Как показали эксперименты 11 и 12, изменения плотности частиц донных наносов  на ±5 % относительно базового значения также не оказывают заметного влияния на изменения площадей взмучивания (изменения не превышают 2 %).

Таким образом, наиболее значимыми в используемой модели взмучивания донных осадков являются параметры τ_ce и M_0, определяющие интенсивность вертикального потока частиц со дна.

Интегральный показатель типа B_kp представляется возможным использовать для уточнения значений входных параметров модели взмучивания по данным дистанционного зондирования для отдельных районов Азовского моря. Для этого необходимо по космическим снимкам с использованием регионального алгоритма оценки концентрации взвешенного вещества [18] получить площади взмучивания и рассчитать показатель B_kp. Далее на основе серии модельных расчетов подобрать значение этого показателя, которое близко к значению показателя, полученного по данным дистанционного зондирования. Задача оценки параметров модели взмучивания донных осадков в мелководных акваториях представляет самостоятельный интерес и рассматривается авторами как продолжение исследований, представленных в данной статье.

Заключение

На основе комбинации численных моделей ветрового волнения и течений исследован процесс взмучивания донных илистых осадков в восточной части залива Сиваш в период экстремального шторма 10–13 ноября 2007 г. Анализ результатов численного моделирования показал следующее.

Наиболее сильное взмучивание донных осадков в заливе происходит в период максимального усиления ветрового волнения и длится в течение 3–4 ч. В период шторма вертикальные профили концентрации илистой взвеси С слабо зависят от глубины. Максимальные перепады между значениями С на поверхности и у дна составляют не более 0,1 мг/л.

Придонные волновые напряжения ${\text{τ}}_{\text{w}}$ вносят определяющий вклад в формирование областей взмучивания донных осадков в заливе. При их учете на 80 % площади залива создаются условия для формирования областей взмучивания (${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$). Если при моделировании волновые напряжения не учитываются, то общая площадь взмучивания сокращается в ~ 4 раза.

Течения также могут формировать придонные напряжения ${\text{τ}}_{\text{c}}\ge {\text{τ}}_{\text{ce}}$, однако это происходит на небольших участках акватории вблизи берега, в то время, как волны генерируют придонные напряжения ${\text{τ}}_{\text{w}}\ge {\text{τ}}_{\text{ce}}$ на большей части залива. Таким образом, при расчете суммарного придонного напряжения ${\text{τ}}_{\text{b}}$ необходимо учитывать как волновые напряжения, так и напряжения, вызванные непосредственно течениями.

На основе многовариантных расчетов установлено, что используемая модель взмучивания илистых осадков наиболее чувствительна к вариациям значений параметров τ_ce и M₀, которые определяют интенсивность вертикального потока частиц со дна бассейна.

Результаты работы могут быть полезны при планировании и интерпретации натурных экспериментов в заливе Сиваш и других мелководных районах Азовского моря.

 

¹ USER MANUAL SWAN Cycle III version 41.45. Delft University of Technology. URL: https://swanmodel.sourceforge.io/online_doc/swanuse/swanuse.html (date of access: 03.06.2024).

² Soulsby R. L. Dynamics of marine sands: a manual for practical applications. London : Tomas Telford Services, 1997. 249 p.

³ URL: apps.ecmwf.int (дата обращения: 03.06.2024).

⁴ Стащук М. Ф., Супрычев В. А., Хитрая М. С. Минералогия, геохимия и условия формирования донных отложений Сиваша. Киев : Наукова думка, 1964. 174 с.

Об авторах

Владимир Владимирович Фомин

Морской гидрофизический институт РАН

Email: v.fomin@mhi-ras.ru
ORCID iD: 0000-0002-9070-4460
SPIN-код: 6992-2688
ResearcherId: H-8185-2015

заведующий отделом вычислительных технологий и математического моделирования, доктор физико-математических наук

Россия, Севастополь

Елена Васильевна Иванча

Морской гидрофизический институт РАН

Автор, ответственный за переписку.
Email: ev_ivancha@rambler.ru
ORCID iD: 0000-0001-8328-4915
SPIN-код: 9981-7991
ResearcherId: AAG-9265-2020

младший научный сотрудник, отдел вычислительных технологий и математического моделирования

Россия, Севастополь

Антон Александрович Полозок

Морской гидрофизический институт РАН

Email: polozok.umi@gmail.com
ORCID iD: 0000-0002-0825-8669
SPIN-код: 5745-1021
ResearcherId: ADJ-1790-2022

старший инженер-программист, отдел вычислительных технологий и математического моделирования

Россия, Севастополь

Список литературы

Дополнительные файлы