CC BY
24
ТЕРМОГИДРОДИНАМИКА ОКЕАНА
УДК 551.465
Двухслойная вихреразрешающая модель ветровых течений в Черном море
© 2015 А.А. Павлушин, Н.Б. Шапиро, Э.Н. Михайлова,
Г.К. Коротаев
Морской гидрофизический институт РАН, Севастополь, Россия Поступила в редакцию 23.04.2015 г.
Процессы формирования и изменчивости циркуляции в Черном море исследуются с помощью двухслойной вихреразрешающей модели. В качестве возбуждающей силы рассматривается только действие ветра, сток энергии происходит за счет горизонтальной турбулентной вязкости, трения между слоями и придонного трения. Модель, в силу своей нелинейности, позволяет описывать, наряду с крупномасштабной циркуляцией, вихревые структуры, образование которых связано с гидродинамической неустойчивостью течений, особенностями береговой линии, рельефа дна и другими факторами. В данной работе обсуждаются результаты численных экспериментов, в которых для лучшего понимания процессов массопереноса и вихреобра-зования в море рассматривалось движение, возбуждаемое стационарным циклоническим ветром, и не учитывались рельеф дна, трение между слоями и ^-эффект. Исследована относительная роль горизонтальной вязкости и придонного трения. Показано, что модель воспроизводит основные черты циркуляции, наблюдаемой в Черном море, а главное - позволяет имитировать образование и трансформацию вихрей.
Ключевые слова: Черное море, вихреразрешающая модель, гидродинамическая неустойчивость, горизонтальное и придонное трение.
Введение. Цель настоящей работы - исследовать механизмы формирования течений и вихревых структур в Черном море. Для этого используется нелинейная двухслойная вихреразрешающая численная модель, аналогичная модели в одной из первых работ Холланда - Лина [1], посвященной моделированию синоптических вихрей в океане.
Исследованию и моделированию вихревой структуры циркуляции в Черном море посвящено большое число работ [2 - 7]. Удалось выявить физико-географические особенности синоптических и мезомасштабных вихрей, предложены возможные механизмы их формирования. Показано, что вихре-разрешающие численные модели неплохо воспроизводят вихревые структуры, которые по параметрам соответствуют реально наблюдаемым. Речь идет главным образом о наблюдениях с искусственных спутников Земли за уровнем моря (альтиметрия) и поверхностной температурой.
В то же время сложность наблюдаемой картины, как и сложность внешнего воздействия, существенно затрудняют интерпретацию полученных результатов. В частности, трудно отделить процессы, связанные с внутренней динамикой, обусловленной нелинейностью, неустойчивостью течений, рельефом дна, конфигурацией берегов, от процессов, вызванных непосредственно внешним воздействием, в данном случае ветром и его пространственно-временной изменчивостью.
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
3
Представляется поэтому целесообразным вернуться к более простым моделям для того, чтобы в «чистом» виде выявить роль отдельных факторов в образовании вихревых структур в Черном море. Предполагается начать с экспериментов, в которых число задаваемых факторов минимизировано, и продолжить исследование далее, усложняя модель за счет включения новых факторов.
Постановка задачи. Будем считать, что основным фактором, влияющим на сезонную и межгодовую изменчивость течений в Черном море, является ветер. Поскольку Черное море существенно бароклинное, то возникают сдвиги скорости не только по горизонтали, но и по вертикали, что может приводить как к баротропной, так и к бароклинной неустойчивости течений.
Двухслойная жидкость со слоями различной плотности является простейшей моделью бароклинного моря, которая в достаточной мере отражает стратификацию вод Черного моря, так что вполне естественным является ее использование для моделирования, в первом приближении, динамических процессов в море.
Итак, рассмотрим движение, возбуждаемое ветром в двухслойном море. Плотность воды в верхнем слое рх = const, в нижнем слое р2 = const, причем р2 > pi. Поверхность раздела слоев z = h(x, y, t) > 0 представляет собой поверхность тока, т. е. слои между собой не смешиваются; x, y, z - координаты по осям, направленным на восток, север и вертикально вниз соответственно; t - время. Толщина верхнего слоя h = h — ^, где £(x, y, z) - уровень, отсчитываемый вниз от невозмущенной поверхности моря, толщина нижнего слоя h 2 = H — h, где H(x, y) - глубина моря.
Используя приближения гидростатики, Буссинеска, yS-плоскости и предполагая, что скорости течения не изменяются по глубине в пределах слоя, уравнения движения и неразрывности запишем в следующем виде:
(uh)t + (M2h) x + (viu i h ) y + г(щ — м2) — f (v h ) = gh^x +rx + A1A(ufy),
(1)
(vh)t + (UiVih) x + (vhi)y + r (v — V2) + f (Uih ) = ghCy +Ty + A, A(vh),
(hi )t + (Uih ) x + (Vih ) y = 0, (2)
(U2h2 )t + (U22h2)x + (V2U2h2)y — r (u i — U2) — f Vh2) =
= gh2^x + g'h2 (hi)x — hU2 + Aj A(U2h2),
2
(3)
VК )t + (и2у2К )х + VК )у - Г (V - у2) + /(ы2н2 ) =
= ё^у + ё'К(И1 )у - Г2»2 + А1 Л^2h2),
(И2)1 + (и2к2) х + (у2И2) у = 0, (4)
где (и,у1 ),(и2,у2) - горизонтальные компоненты скорости течений в верхнем и нижнем слоях соответственно; г1, г2 - коэффициенты трения на поверхно-
4
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
сти раздела и у дна; / = /0 + /у - параметр Кориолиса, / = 10—4 1/с, / = 2 -10—131/(см-с); g = 980 см/с2 - ускорение свободного падения; g' = g(р2 — Рх)/ р2 ; - компоненты тангенциального напряжения ветра;
А - коэффициент горизонтальной турбулентной вязкости.
Следуя Холланду - Лину [1], используем приближение «твердой крышки», т. е. используем интегральное уравнение неразрывности
их + Гу = 0, (5)
что позволяет ввести интегральную функцию тока :
и = —¥у, V = ¥х, (6)
где и = Пуку + , V = у^ку + - компоненты полного потока.
Граничные условия на поверхности моря, поверхности раздела слоев и на дне учтены при выводе уравнений (1) - (4). На боковых границах задаются условия прилипания, т. е. равенства нулю обеих компонент скорости течения:
щ = Vк = 0, к = 1; 2. (7)
Сток рек в море и водообмен через проливы не учитываются. В начальный момент времени жидкость находится в состоянии покоя, поверхность раздела и поверхность моря горизонтальны.
Численная модель. Для численной аппроксимации используется конечно-разностная схема, основанная на бокс-методе с сеткой В (по терминологии Аракавы), двухслойной схеме интегрирования по времени, неявной аппроксимации трения у поверхности раздела и у дна и силы Кориолиса. Такая схема обычно применяется авторами в многослойной квазиизопикнической [8 -10] и многоуровенной [11 - 13] моделях.
Адвективные члены в уравнениях неразрывности и уравнениях движения аппроксимируются по-разному. В уравнениях неразрывности используется аппроксимация направленными разностями с первым порядком точности, в уравнениях движения для аппроксимации адвективных членов реализована схема Лакса - Вендроффа второго порядка точности. Учет схемной вязкости в уравнениях неразрывности обеспечивает устойчивость численной схемы и положительную определенность толщины слоев. Применение схемы Лакса -Вендроффа для вычисления толщины слоев приводит к неустойчивости численной схемы при расчете на большой срок.
В рассматриваемых экспериментах используется следующая последовательность вычислений. Предположим, что распределения всех полей в п-й момент времени известны. Вначале из уравнений неразрывности находятся толщины верхнего и нижнего слоев:
к"+1 к"
- = кк. — [(щкк )х + (Уккк )у ]", к = 1; 2. (8)
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
5
Затем вычисляются компоненты скорости течения и уровень моря в (п +1)-й момент времени. Конечно-разностный аналог уравнений движения запишем в виде:
-.и+1
(9)
(и1К1) чИ+1 . /- \И+1 1 П+1 П+1 , тП
-7--/^О + Г1(и1 - П2) = ёИ1 + Ъ ,
Аt
+/(пЛГ1 + ^ - у2)"+1 = ghln+1cy+1 + м;,
а
^Кт^ -/(УгКТ1 -ГП -и2)"+1 + г2п'П+1 = gh2¡+1a+1 + ^
Аt
+ /(и2к2)п+ -фх -v2)2+ + Т2^ = ghn+1Cy+1 + МП,
Аt
тП и /й тП и /й
где Ц,Мх,ь2,М2 - известные в данный момент выражения.
С помощью несложных выкладок можно получить соотношения, связывающие компоненты скорости течения с наклонами уровня:
пТ1 = и* +©к£2+1 -ЛкС1, у2+1 = V* +ЛкС2+1 +©к^;+1, (10)
* *
здесь и*,V* - результаты расчета скорости течения при нулевых наклонах уровня и ненулевых значениях Ь2к,М;; функции ©к, Лк - результаты расчета при наклонах уровня £х = 1, £у = 0 и нулевых значениях Ь2к,М; .
Умножая соотношения (10) на К;+1 и суммируя их по к, получим выражения для полных потоков
(11)
и;+1 = и* + &к£;+1 -\£;+1, у;+1 = V* + лкс;+1 + 4С+1,
где Uk,Vk*,^k,Лк - просуммированные по к функции и*,V*,©к,Лк .
Разрешая полученные соотношения относительно наклонов уровня и исключая затем уровень перекрестным дифференцированием, с учетом конечно-разностного представления производных £х, £ получим систему линейных алгебраических уравнений для интегральной функции тока на 9-точечном шаблоне. Эти уравнения решаются методом верхней релаксации с эмпирическим подбором оптимального параметра релаксации.
Вычислив \уп+, находим компоненты полного потока, наклоны уровня и затем компоненты скорости для каждого слоя. Зная наклоны уровня £х = А, £у = В, можно рассчитать сам уровень, используя уравнение Пуассона А£ = А + В при известной на границе нормальной производной уровня.
6 МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
Уравнение Пуассона также решается методом верхней релаксации. Поскольку уровень определяется с точностью до константы, а средний по площади уровень должен быть равен нулю, то из вычисленного уровня следует вычесть его среднее по площади значение.
Численные эксперименты. Данная постановка задачи подразумевает, что факторами, влияющими на формирование циркуляции в море, являются пространственно-временная изменчивость напряжения ветра, форма бассейна и особенности береговой черты, интенсивность сил трения, рельеф дна, в-эф-фект, разница плотностей в верхнем и нижнем слоях, начальная толщина верхнего слоя. В этой и следующих работах постараемся рассмотреть и проанализировать влияние указанных факторов на результаты численного моделирования как по отдельности, так и в совокупности друг с другом.
В серии экспериментов, о которых пойдет речь в этой статье, из модели исключено влияние рельефа дна, в-эффекта и трения на поверхности раздела слоев, дно считается горизонтальным (Н = 2200 м).
Для расчетов использовалась прямоугольная сетка с шагами Лх = 4 км, Лу = 3,5 км, шаг по времени Л! = 12 мин.
Приведенное ускорение свободного падения (ё' = 3,2 см/с2) и начальная толщина верхнего слоя (к0 = 175 м) выбирались согласно данным наблюдений в Черном море. В этом случае радиус деформации Россби (Л = у]1/) равен 23,4 км и почти в шесть раз превышает шаг сетки, что необходимо для удовлетворительного описания вихрей [14].
Движение возбуждалось стационарным циклоническим ветром, который рассчитывался следующим образом:
тх = Т + уТ -т%)/ 1у, ту =4 + х(туЕ-4)/4, (12)
х х у у X V
где т3,ты,Тф,тЕ - значения т ,ту в соответствии со сторонами света на границах прямоугольной области [ 1х х I ], в которую вписано Черное море. Завихренность ветра при таком распределении постоянна для всех точек бассейна. На рис. 1 изображено циклоническое поле ветра, которое использовалось в данной серии экспериментов, оно получено по формулам (12) при
Т =тЕ = 0,5 см2/с2, тхы =ту =-0,5 см2/с2.
Вычисления проводились на срок от 5 до 25 лет с различными значениями коэффициентов горизонтальной турбулентной вязкости А1 и придонного трения г2 . Расчеты прекращались, когда решение выходило на стационарный или квазистационарный режим. Коэффициенты А1 выбирались в диапазоне от достаточно малых значений (102 см2/с), при которых проявляются нелинейные эффекты, до сравнительно больших (107 см2/с), когда роль нелинейности незначительна. Коэффициенты придонного трения г2 выбирались в диапазоне от 0 до 1 см/с. Коэффициент г2, вообще говоря, характеризует
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
7
время приспособления к внешнему воздействию течений в нижнем слое ^* « ^ /г2, т. е. г2 = 1 см/с соответствует ^* « 2 сут при Л2 = 2 км.
Рис. 1. Поле тангенциального напряжения ветра
Результаты экспериментов. По результатам проведенных расчетов были построены графики зависимости от времени средних по площади доступной потенциальной и кинетической энергий, а также пространственные распределения полей \ и течений в верхнем и нижнем слоях щ, и2 соответственно. Наиболее интересные результаты размещены на сайте http://blacksea-model.ru.
Доступная потенциальная (БРЕ) и общая кинетическая (КЕ) энергии определялись как
БРЕ = (рг£(^ -\)2/^, КЕ = КЕХ + КЕ2,
КЕ1 =(РА(и1 + ^2)/2), КЕ 2 =[р2Ь2(и1 + у1)/2),
где КЕ2, КЕ2 - кинетическая энергия верхнего и нижнего слоев; Е = БРЕ+КЕ - суммарная энергия; угловые скобки означают осреднение по площади.
На рис. 2 приведены графики зависимости видов энергии от времени при фиксированном и достаточно малом коэффициенте придонного трения (г2 = 0,01 см/с) и различных значениях коэффициента горизонтальной вязкости Л1.
Как видно из рисунка, при Л = 107 см2/с решение имеет стационарный вид с близким к нулю значением кинетической энергии КЕ2 , что свидетельствует о практическом отсутствии течений в нижнем слое в стационарном случае. Энергии КЕ1, КЕ , БРЕ имеют наименьшие значения из рассматриваемых случаев.
8
МОРСКОЙ ГИДРОФИЗИЧЕСКИИ ЖУРНАЛ № 5 2015
12 8 4 О
- АЕ, кДж/м2 1,11111,
-------: 4 5/Х /К
/\ У>
: ! в
Г---
—1—1— —"—1—1— —1—1—1— —1—1—1— 1-^1-1 —1—
Рис. 2. Изменение во времени средних по площади составляющих энергии: Е(а), БРЕ(б), КЕ(в), КЕ1(г), КЕ2(д) при г2 = 0,01 см/с и различных значениях А1 (1 - Д = 107 см2/с; 2 -
А = 106 см2/с; 3 - А, = 5 • 105 см2/с; 4 - Л, = 105 см2/с; 5 - Л, = 104 см2/с; 6 - Д = 103 см2/с)
- ЭРЕ, к^ж/м2 1 1 Ы ^^Г : : 2
/ 1 1 : 1 б
^-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
9
При дальнейшем уменьшении А до 5-105 см2/с колебания энергии в равновесном режиме принимают более стохастический характер, что связано с возрастающей ролью нелинейности в уравнениях движения. Кинетическая энергия системы при этом увеличивается, а доступная потенциальная энергия, наоборот, уменьшается. Уменьшение БРЕ оказывает большее влияние на изменение суммарной энергии Е, чем увеличение КЕ, в результате суммарная энергия незначительно уменьшается по сравнению с экспериментом при А = 10б см2/с.
Графики энергий при значениях А, меньших 5-105 см2/с, в статистически равновесном режиме характеризуются наличием колебаний большой амплитуды, которые происходят в определенных пределах и также носят стохастический характер. Причем интервалы, в которых изменяется энергия, оказываются одинаковыми для различных А. Такое поведение модели при различных, сравнительно малых значениях коэффициента горизонтальной турбулентной вязкости, вероятнее всего, связано с влиянием схемной вязкости при расчете толщины слоев по используемой схеме, о чем уже говорилось ранее.
На рис. 3 представлены пространственные распределения толщины верхнего слоя \ и течений в верхнем слое щ , полученные в численных экспериментах при выходе на статистически равновесный режим в случаях А = 107 см2/с и А = 10б см2/с.
р1....................................................... Г I I I I I I I I I II I I I I I I I I I I ■ I I I I I I I I I I ■ I I I I I Н т г.....г 1-1 I I I I
О 100 200 300 400 500 600 700 800 900 1000 X. км о 100 200 300 400 500 600 700 800 900 1000 X км
Рис. 3. Пространственное распределение толщины верхнего слоя Н1, м (а), скорости течений в верхнем слое щ (б) при г2 = 0,01 см/с, А1 = 107 см2/с и толщины верхнего слоя Н1, м (в), скорости течений в верхнем слое щ (г) при г2 = 0,01 см/с, А1 = 106 см2/с (области с большими значениями толщины слоя обозначены более темным цветом)
При А = 107 см2/с (рис. 3, а, б) в верхнем слое имеет место стационарная циклоническая циркуляция с двумя круговоротами. Максимальные скорости
10
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
течения до 20 см/с отмечаются у западной границы Крымского п-ова. В нижнем слое течения практически отсутствуют. Поле \ характеризуется наличием двух куполов (областей с малыми значениями толщины слоя), которые соответствуют центрам циклонических круговоротов в поле течений. От центров круговоротов по направлению к берегу толщина верхнего слоя увеличивается. Такая ситуация находится в полном соответствии с линейной теорией течений в двухслойной жидкости, когда циркуляция определяется завихренностью ветра и конфигурацией бассейна.
При уменьшении А до 106 см2/с течения усиливаются и возникают возмущения в виде волн, которые перемещаются в направлении потока, но картина течений в основном сохраняется, а именно имеет место циклоническая циркуляция с несколькими ядрами. Максимальные скорости течений доходят в верхнем слое до 60 см/с, в нижнем слое, который также приходит в движение, - до 10 см/с. На рис. 3, в, г приведены мгновенные распределения \ и
щ , на которых видно, как трансформируются поля по сравнению с предыдущим экспериментом. На странице в интернете http://blacksea-model.ru/r0.01l6/hb17 можно посмотреть этот процесс в движении.
При дальнейшем уменьшении А до 5-10 см /с (рис. 4) возмущения в поле скорости становятся более значительными и начинают формироваться антициклонические вихри, причем как в верхнем, так и в нижнем слоях. Зоны образования этих вихрей привязаны к особенностям береговой линии, сами вихри хорошо видны в поле скорости (рис. 4, в, г). По сравнению с предыдущим случаем (рис. 3) течения в верхнем слое увеличиваются незначительно, в нижнем слое - более существенно. Антициклонические вихри в поле \ проявляются в виде локальных максимумов толщины слоя. Видно, что осреднен-ное за год поле \ (рис. 4, б) отличается от мгновенного поля (рис. 4, а), хотя не очень значительно.
Рис. 4. Пространственное распределение толщины верхнего слоя к1, м (а), среднегодовой толщины , м (б), скорости течений в верхнем слое щ (в) и скорости течений в нижнем слое и2 (г) при г2 = 0,01 см/с и А1 = 5105 см2/с для 05.06.0004 г. (области с большими значениями толщины верхнего слоя обозначены более темным цветом)
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
11
Наиболее интересная картина наблюдается в экспериментах с коэффициентом Л1 = 105 см2/с (рис. 5). Четко выделяются два крупных циклонических круговорота («Очки Книповича»), на периферии которых образуются довольно интенсивные антициклонические вихри. Антициклонам соответствуют локальные максимумы \ . Особенности циркуляции хорошо видны как в топографии верхнего слоя, так и в полях течений в верхнем и нижнем слоях. Образование антициклонических вихрей, связанное с гидродинамической неустойчивостью, носит нестационарный характер. Они спорадически возникают, перемещаются в направлении потока вдоль берега или отходят от берега вместе с потоком. На приведенных распределениях скорости течения можно увидеть проявление баротропизации процесса, а именно образование однонаправленных вихрей в одних и тех же местах одновременно в верхнем и нижнем слоях.
Рис. 5. Пространственное распределение толщины верхнего слоя к1, м (а), среднегодовой толщины , м (б), скорости течений в верхнем слое и (в) и скорости течений в нижнем слое и2 (г) при г2 = 0,01 см/с, А1 = 105 см2/с для 05.05.0004 г. (области с большими значениями толщины верхнего слоя обозначены более темным цветом)
О нестационарном характере циркуляции в статистически равновесном режиме можно судить, сравнивая мгновенное и осредненное за год распределения к1. Так, в осредненном за год поле антициклонические вихри, связанные с неустойчивостью течений, вообще не видны. Проявляются только крупные циклонические круговороты в западной и восточной частях моря.
Далее, чтобы оценить влияние придонного трения на решение задачи,
были проведены расчеты при фиксированном г2 > 0,01 см/с.
Л1 и разных значениях
12
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
Использование в расчетах малых значений г2 < 0,001 см/с приводит к появлению в нижнем слое нереально больших скоростей течений, что противоречит данным наблюдений и может вызывать неустойчивость схемы при длительных расчетах.
Результаты расчетов при А1 = 105 см2/с и различных г2 приведены на рис. 6. Видно, что при увеличении придонного трения увеличивается ОРЕ и уменьшается КЕ . Кинетическая энергия верхнего слоя КЕУ слабо зависит от Г, в противоположность кинетической энергии нижнего слоя КЕ2, которая резко уменьшается с ростом коэффициента придонного трения. При Г ^ 1,0см/с КЕ2 практически обращается в нуль. Таким образом, общая кинетическая энергия системы уменьшается за счет уменьшения кинетической энергии нижнего слоя.
На рис. 7 и 8 представлены распределения толщины верхнего слоя и скорости течений для экспериментов с г = 0,1 см/с иг = 1,0 см/с. Увеличение коэффициента приводит к существенному изменению картины циркуляции. Циклоническая циркуляция в верхнем слое представляет собой теперь один большой круговорот с двумя центрами в западной и восточной частях бассейна. Если интенсивность течений в верхнем слое не очень чувствительна к уменьшению придонного трения, то течения в нижнем слое изменяются значительно. Так, при г = 0,1 см/с скорости течения достигают 10 см/с, а при Г = 1,0 см/с - только 2 см/с. Также с увеличением коэффициента придонного трения уменьшается интенсивность образования антициклонических вихрей, происходит их локализация восточнее Синопа и западнее Крыма.
В мгновенных и осредненных полях толщины верхнего слоя \ и И\ крупномасштабный циклонический круговорот с двумя ядрами хорошо выражен в виде уменьшения толщины слоя. При увеличении придонного трения уменьшаются различия между мгновенными и осредненными полями ^ , что свидетельствует о более устойчивой циркуляции.
Пространственное распределение характеристик при г2 > 1,0 см/с практически не отличается от приведенного на рис. 8, за исключением скоростей в нижнем слое, которые оказываются близкими к нулевым.
Таким образом, проведенные численные эксперименты показали, что рассмотренная модель позволяет получить решение в большом диапазоне значений коэффициентов горизонтальной турбулентной вязкости и придонного трения, неплохо описывает основные черты циркуляции в Черном море, в том числе воспроизводит образование и трансформацию вихрей, связанных с гидродинамической неустойчивостью.
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
13
Рис. 6. Изменение по времени средних по площади составляющих энергии: Е(а), ВРЕ(б), КЕ(в), КЕ1(г), КЕ2(д) при Д = 105 см2/с и различных значениях г2 (1 - Г = 0,01 см/с; 2 -Г = 0,02 см/с; 3 - г2 = 0,1 см/с; 4 - г2 = 1,0 см/с; 5 - г2 = 10,0 см/с)
14
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
Рис. 7. Пространственное распределение толщины верхнего слоя ^, м (а), среднегодовой
толщины кх , м (б), скорости течений в верхнем слое щ (в) и скорости течений в нижнем слое
и2 (г) при г = 0,1 см/с, Л] = 105 см2/с для 10.03.0008 г. (области с большими значениями толщины верхнего слоя обозначены более темным цветом)
Рис. 8. Пространственное распределение толщины верхнего слоя ^, м (а), среднегодовой
толщины кг , м (б), скорости течений в верхнем слое щ (в) и скорости течений в нижнем слое
и2 (г) при г = 1,0 см/с, Л1 = 105 см2/с для 05.10.0004 г. (области с большими значениями толщины верхнего слоя обозначены более темным цветом)
Чтобы окончательно определиться с выбором коэффициентов, проверим, как выполняется энергетический баланс в модели. Дело в том, что в отличие МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015 15
от модели Холланда - Лина [1] используемая численная модель, строго говоря, не является энергосбалансированной, т. е. не гарантируется интегральный баланс энергии соответствующей аппроксимацией адвективных членов (согласно способу Аракавы) в уравнениях движения.
В рассматриваемой задаче изменение суммарной энергии системы по времени Et определяется совместной работой W сил тангенциального напряжения
ветра W, горизонтальной турбулентной вязкости Жж и трения о дно WR2 :
Et = W = WT + Wal + W2, WT=(Pl(u1rx + )),
(13)
WAL ={pA (uk Aukhk + vk Avkhk ))> k = 1; 2, WR2 =[p2r2 (u22 + v22 )).
С использованием результатов проведенных численных экспериментов были рассчитаны значения левой и правой частей уравнения энергии. На рис. 9 приведены графики зависимости от времени Et и W = WT+ WAL + WR2 при r2 = 0,01см/с и различных значениях Аг ■ Видно, что во всех случаях существует невязка S - различие между рассчитанными Et и W. Эта невязка всегда отрицательна и связана с особенностями используемой численной схемы.
Лучшие, на наш взгляд, результаты получаются при использовании коэффициентов Aj и r2 в достаточно узком диапазоне: At - в пределах 105 -1,5-105 см2/с, r2 - в пределах 0,005 - 0,01 см/с. При этом невязка оказывается наименьшей и близкой по величине к W и W . При r , меньших 0,005 см/с, значительно возрастают скорости в нижнем слое, а при больших значениях A хуже проявляются эффекты бароклинной неустойчивости течений.
Поля h ,u ,u , полученные в экспериментах с A и r из указанных диапазонов, оказываются подобными полям, приведенным на рис. 5. В море наблюдаются устойчивые крупномасштабные круговороты, на периферии которых спорадически образуются вихри синоптического масштаба. В крупномасштабных суббассейновых вихрях завихренность течений совпадает с завихренностью ветра, что соответствует общепринятым представлениям о течениях в Черном море [2, 15]. Отметим, что рассчитанные скорости течений в нижнем слое получаются больше, чем реально наблюдавшиеся. Согласно данным наблюдений, в глубоководной части Черного моря море скорости течений не превышают 10 см/с, а в численных экспериментах они могут достигать 20 - 40 см/с, разумеется, при интенсивном ветре и малом придонном трении.
Эксперименты показали, что интенсивность рассматриваемого циклонического ветра не влияет на качественную картину циркуляции в бассейне, а приводит лишь к пропорциональному изменению энергетических характеристик модели. На рис. 10 приведены временные графики составляющих энергии и работы сил (13) по результатам двух расчетов с разными значениями
16
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
тангенциального напряжения ветра (12): при тх3 =ту = 1,0 см2/с Тхн =туш =-1,0 см2/с2 и тх3 =туЕ = 0,5 см2/с2, тхн =хуш =-0,5 см2/с2.
Рис. 9. Графики производной по времени суммарной энергии у, совместной работы Ш сил тангенциального напряжения ветра, горизонтальной турбулентной вязкости, придонного трения и невязки 8 при г2 = 0,01 см/с и различных значениях Л1 (а - Лг = 103 см2/с; б -
Л = 105 см2/с; в - Л = 5 • 105 см2/с; г - Лг = 106 см2/с)
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
17
Рис. 10. Изменение по времени составляющих энергии Е, КЕ, БРЕ (а, в) и работы сил тангенциального напряжения ветра ^ , придонного трения , горизонтальной турбулентной вязкости ^ (б, г) при различных значениях тангенциального напряжения ветра: а, б -т* =те = 1,0 см2/с2, т* = ту = -1,0 см2/с2; в, г - т* = ту = 0,5 см2/с2, тхм = т£ = -0,5 см2/с2 (г = 0,007 см/с, А, = 1,4 -105 см2/с)
Видно, что характер поведения графиков одинаковый, за исключением того, что при более сильном ветре значения составляющих энергии и работы сил оказываются больше, чем при ветре меньшей интенсивности. Большими оказываются и амплитуды колебаний характеристик.
18
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
Наличие этих колебаний, очевидно, связано с гидродинамической неустойчивостью течений, которая приводит к меандрированию струйных течений и образованию вихрей. Этот процесс носит нестационарный характер и хорошо выражен на графиках работы сил горизонтальной вязкости и тангенциального напряжения ветра (рис. 10, б, г). Эти графики изменяются в проти-вофазе с небольшим отставанием работы силы ветра. Наименьшее влияние процесс вихреобразования оказывает на работу силы придонного трения.
Примечательно, что при стационарном ветре наблюдаются значительные изменения работы силы тангенциального напряжения ветра. Это происходит в результате рассогласования полей ветра и течений в верхнем слое вследствие образования синоптических и мезомасштабных вихрей противоположной ветру завихренности. Таким образом, в результате взаимного влияния работы силы ветра и процессов вихреобразования возникает механизм автоколебаний, при котором увеличивающийся поток энергии от ветра вызывает усиление интенсивности вихреобразования, вследствие чего деформируется поле течений. Рассогласование полей ветра и течений приводит к снижению притока энергии от ветра. Далее синоптические вихри затухают под действием сил трения, поле течений подстраивается под ветер, работа силы ветра начинает расти и т. д.
Из-за сильной нелинейности амплитуда автоколебаний и их период могут меняться. Эта изменчивость зависит от силы ветра. Так, при более сильном ветре (рис. 10, а, б) на фоне сравнительно малых колебаний спорадически возникали резкие падения значений характеристик («провалы»).
Для одного из таких «провалов», наблюдавшегося в эксперименте с более сильным ветром, приведены поля \ для нескольких последовательных моментов времени (рис. 11). Более подробно формирование и эволюцию поля кх для этого эксперимента можно посмотреть в интернете по ссылке http://blacksea-model.ru/experiment_a 1.html.
Рассмотрим, как выражен описанный выше механизм автоколебаний в изменчивости толщины верхнего слоя. В начале рассматриваемого периода (рис. 11, а) в поле Их наблюдаются два купола, соответствующие западному и восточному циклоническим круговоротам. Минимальная толщина в центрах циклонов составляет 90 и 115 м для западного и восточного круговоротов соответственно. На периферии циклонов возле берега наблюдаются антициклонические вихри. В это время кинетическая и доступная потенциальная энергии системы достаточно велики. Далее в течение 5 мес картина существенно не меняется, амплитуда колебаний энергии небольшая.
В начале июня возникает ситуация (рис. 11, б), при которой два крупномасштабных циклона начинают движение навстречу друг другу, а затем соединяются в центральной части моря (рис. 11, в). В районе Синопа и западнее Крыма в это время возникают антициклонические вихри, которые в дальнейшем усиливаются и перемещаются вглубь бассейна. Как видно на рис. 11, г, д, антициклонические вихри занимают положение, в котором раньше располагались центры циклонов. Область циклонической циркуляции при этом распадается на три части, наиболее интенсивный циклон наблюдается в центре бассейна. На графиках энергии (рис. 10, а, б) в это время наблюдается
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015 19
резкое падение значений, которое обусловлено уменьшением работы силы ветра.
v >Р1 ........................................................ .......................................................
О 100 200 300 400 500 600 700 800 900 1000 Х, кмО 100 200 300 400 500 600 700 800 900 1000 X. км
Рис. 11. Толщина верхнего слоя A¡, м для различных моментов времени в эксперименте со стационарным ветром при г2 = 0,007 см/с, Д = 1,4 -105 см2/с (области с б0льшими значениями толщины слоя обозначены более темным цветом)
В следующие месяцы (рис. 11, е, ж, з) поле h1 становится более равномерным, антициклонические вихри в центрах западной и восточной частей моря затухают, начинает увеличиваться приток энергии от ветра и происходит постепенная перестройка поля h1 к виду, когда в бассейне наблюдаются два крупных циклона с антициклонами на периферии.
Заключение. Проведенные эксперименты показали, что нелинейная двухслойная модель хорошо воспроизводит основные черты крупномасштабной
20 МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
циркуляции в Черном море [3, 4, 15], в частности струйное Основное Черноморское течение, а также формирование и трансформацию вихревых структур: крупномасштабных циклонических вихрей («Очки Книповича») и вихрей синоптического масштаба, обусловленных меандрированием Основного Черноморского течения и процессами гидродинамической неустойчивости.
Наблюдаемые в Черном море вихревые образования могут формироваться за счет внутренней нелинейной динамики даже при воздействии стационарного ветра без учета его мелкомасштабных особенностей. Выявлен и описан механизм автоколебаний, который связывает приток энергии от ветра с интенсивностью вихреобразования в море.
Определен диапазон возможных значений коэффициентов горизонтальной турбулентной вязкости и придонного трения, при которых в модели воспроизводятся эффекты гидродинамической неустойчивости течений.
Из полученных результатов следует, что в процессе вихреобразования участвуют оба слоя. Наиболее интенсивное и соответствующее данным наблюдений вихреобразование получается в случае, когда течения в нижнем слое достаточно развиты. Однако это противоречит данным наблюдений о течениях в глубоководной части Черного моря, которые являются сравнительно слабыми (их скорости не превышают 10 см/с). Это может быть связано с различными причинами, например со стационарностью ветра, отсутствием рельефа дна, наличием только двух слоев.
Использованная в данной версии модели численная схема имеет одно существенное достоинство - она позволяет производить расчеты в большом диапазоне задаваемых параметров и получать устойчивые решения. Недостатком схемы можно считать то, что она не обеспечивает выполнение энергетического баланса.
В дальнейшей работе планируется проведение как в рамках данной модели, так и с применением энергосбалансированной модели, основанной на схеме Аракавы, экспериментов, посвященных исследованию влияния ß-эф-фекта, формы бассейна, особенностей береговой черты, рельефа дна, сезонной изменчивости ветра и т. д.
Задачей планируемых исследований является анализ результатов моделирования вихревых структур в Мировом океане, полученных при использовании различных численных схем. Предполагается, что это позволит дать ответы на вопросы относительно адекватности воспроизводимой в численных экспериментах циркуляции, механизмов образования и временной эволюции получаемых при этом вихревых структур.
Работы по данной тематике выполнялись при финансовой поддержке РФФИ в рамках научного проекта № 14-45-01044 «р_юг_а». Результаты экспериментов, для лучшего их представления, размещены на сайте http://blacksea-model.ru в виде интерактивных графиков и видеофильмов.
СПИСОК ЛИТЕРАТУРЫ
1. Holland W.R., Lin L.B. On the generation of mesoscale eddies and their contribution to the oceanic general circulation. I. A preliminary numerical experiment. II. A parameter study // J. Phys. Oceanogr. - 1975. - 5, № 4. - P. 642 - 657.
2. Блатов А.С., Булгаков Н.П., Иванов В.А. и др. Изменчивость гидрофизических полей Черного моря. - Л.: Гидрометеоиздат, 1984. - 239 с.
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 5 2015
21
3. Коротаев Г.К., Саенко О.А., Коблинский Ч.Ж. и др. Точность, методология и некоторые результаты ассимиляции альтиметрических данных TOPEX/POSEIDON в модели общей циркуляции Черного моря // Исследование Земли из космоса. - 1998. - № 3. -С. 3 - 17.
4. Korotaev G.K., Saenko O.A., Koblinsky C.J. Satellite altimetry observation of the Black Sea level // J. Geophys. Res. - 2001. - 106, № C1. - P. 911 - 933.
5. Демышев С.Г., Кныш В.В., Коротаев Г.К. Численное моделирование сезонной изменчивости гидрофизических полей Черного моря // Морской гидрофизический журнал. -2002. - № 3. - С. 12 - 27.
6. Демышев С.Г. Численный прогностический расчет течений в Черном море с высоким горизонтальным разрешением // Там же. - 2011. - № 1. - С. 36 - 47.
7. Кныш В.В., Коротаев Г.К., Моисеенко В.А. и др. Сезонная и межгодовая изменчивость гидрофизических полей Черного моря, восстановленных на основе реанализа за период 1971 - 1993 гг. // Изв. РАН. Физика атмосферы и океана. - 2011. - 47, № 3. - С. 433 -446.
8. Михайлова Э.Н., Шапиро Н.Б. Квазиизопикническая слоистая модель крупномасштабной океанической циркуляции // Морской гидрофизический журнал. - 1992. - № 4. -С. 3 - 12.
9. Павлушин А.А., Шапиро Н.Б. Имитация сезонной изменчивости Атлантического океана в квазиизопикнической модели // Там же. - 1996. - № 4. - С. 12 - 25.
10. Шапиро Н.Б. Формирование циркуляции в квазиизопикнической модели Черного моря с учетом стохастичности напряжения ветра // Там же. - 1998. - № 6. - С. 26 - 42.
11. Андросович А.И., Михайлова Э.Н., Шапиро Н.Б. Численная модель и расчеты циркуляции вод северо-западной части Черного моря // Там же. - 1994. - № 5. - С. 28 - 42.
12. Михайлова Э.Н., Шапиро Н.Б. Метод вложенных сеток с улучшенным вертикальным разрешением в прибрежной зоне моря // Там же. - 2013. - № 4. - С. 3 - 12.
13. Михайлова Э.Н., Шапиро Н.Б. Трехмерная негидростатическая модель субмаринной разгрузки в прибрежной зоне моря // Там же. - 2014. - № 4. - С. 28 - 50.
14. Каменкович В.М., Кошляков М.Н., Монин А.С. Синоптические вихри в океане. - Л.: Гидрометеоиздат, 1982. - 264 с.
15. Иванов В.А., Белокопытов В.Н. Океанография Черного моря. - Севастополь: МГИ НАН Украины, 2011. - 212 с.
Two-layer eddy-resolving model of wind currents in the Black Sea A.A. Pavlushin, N.B. Shapiro, E.N. Mikhailova, G.K. Korotaev
Marine Hydrophysical Institute, Russian Academy of Sciences, Sevastopol, Russia
Processes of formation and variability of the Black Sea circulation are investigated using a two-layer eddy-resolving model. The wind effect only is considered as an inducing force; energy sink takes place due to horizontal turbulent viscosity, friction between the layers and the bottom friction. Being nonlinear, the model, alongside with large-scale circulation, permits to describe vortex structures formation of which is related to hydrodynamic instability of currents, coastline features, bottom topography and other factors. Discussed are the results of numerical experiments in which, for better understanding of mass transfer and vortex formation processes in the sea, the motion inducing by stationary cyclonic winds is considered, whereas bottom topography, friction between the layers, and ^-effect are not taken into account. Relative role of horizontal viscosity and bottom friction is studied. It is shown that the model reproduces the major features of the Black Sea circulation, and the most important particular is that it permits to simulate formation and transformation of vortices.
Keywords: Black Sea, eddy-resolving model, hydrodynamic instability, horizontal and bottom friction.
22
МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № S 2OiS
