Модульный подход к расчету циркуляции и приливов в Черном море Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ТЕРМОГИДРОДИНАМИКА ОКЕАНА УДК 551.466.7+551.461

Модульный подход к расчету циркуляции и приливов

в Черном море

© 2016 С.Ф. Доценко*, В.Б. Залесный**, Н.К.В. Санникова*

*Морской гидрофизический институт РАН, Севастополь, Россия

E-mail: [email protected] **Институт вычислительной математики РАН, Москва, Россия E-mail:vzalesny@yandex. ru

Поступила в редакцию 20.10.2015 г. После доработки 10.12.2015 г.

Выполнен анализ циркуляции и приливов в Черном море. Даны оценки основных характеристик приливных колебаний. Помимо этого, с использованием модульного подхода проведены расчеты циркуляции вод Черного моря. Применена с-модель циркуляции океана, разработанная в Институте вычислительной математики РАН. Пространственное разрешение модели по долготе и широте составляет около 4 км. По вертикали задавалось 40 неравномерно распределенных с-уровней. Шаг по времени - 300 с.

В циркуляции Черного моря отчетливо проявляется вихревая структура. Воспроизводится Основное Черноморское течение, которое характеризует общую циклоническую циркуляцию по всему периметру Черного моря, образуя два заметных вихря. Для описания генерации приливов в с-модель циркуляции добавлен модуль, соответствующий приливообразующим потенциалам Луны и Солнца.

Ключевые слова: Черное море, общая циркуляция, приливы, математическая модель, модульный подход, вычислительные эксперименты.

Введение. Высота приливов в открытом океане составляет в среднем 1,2 м. Возникновение приливов и отливов обусловлено притяжением Луны (и в меньшей степени Солнца), действующим на вращающуюся Землю. Прили-вообразующая сила воздействует на всю поверхность Земли (включая сушу). Наиболее заметны приливы и отливы на побережьях океанов и морей.

Самые высокие приливы наблюдаются в Северной Америке в заливе Фанди и частично вдоль побережья штата Мэн (США). Максимальная измеренная здесь высота прилива - 18 м. Наивысший прилив у берегов Российской Федерации высотой 14 м был зарегистрирован в Пенжинской губе Охотского моря.

Характеристики черноморских приливов, как правило, устанавливаются из анализа данных мареографических измерений, отражающих колебания уровня в различных пунктах морского побережья, в частности в таких крупных портах как Одесса, Севастополь, Поти, Констанца, Самсун, Трабзон и др. Анализ этих данных показал, что амплитуда приливно-отливных колебаний поверхности Черного моря не превышает 17 см: в центральных частях западного и восточного участков побережья она равна 9 см, у Крымского п-ова — 2 - 3 см. Таким образом, приливы в Черном море относительно слабы, что и предопределило их недостаточную изученность.

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

3

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

Общая характеристика черноморских приливов. Черное море - внутреннее море Атлантического океана. Проливом Босфор оно соединяется с Мраморным морем, через пролив Дарданеллы - с Эгейским и Средиземным морями, а через Керченский пролив - с Азовским. Средиземноморские приливы затухают в проливах. Периоды черноморских сейш и приливов различаются, поэтому резонансные режимы в Черном и Азовском морях невозможны [4 - 10].

Параметры черноморских приливов устанавливаются по данным мареографов, отслеживающих колебания уровня в различных пунктах побережья Черного моря, в частности в крупных портах [4]. Анализ показал, что приливы и отливы в Черном море выражены слабо вследствие относительной малости площади бассейна. Черноморские проливы неглубокие и сравнительно узкие, что не способствует развитию приливов.

По данным измерений, средний период приливов в Черном море составляет около полусуток. Как уже упоминалось, приливы в Черном море - относительно слабое волновое явление, амплитуда приливно-отливных колебаний уровня Черного моря не превышает 17 см. Пространственная структура полусуточного прилива в Черном море показана на рис. 1.

Рис. 1. Полусуточный прилив в Черном море [11, 12]: а - котидальная карта; б - карта изоам-плитуд

а

б

4

МОРСКОЙ ГИДРОФИЗИЧЕСКИМ ЖУРНАЛ № 1 2016

Значения гармонических постоянных основных приливных гармоник (средняя амплитуда Н и угол положения у волны) для трех черноморских портов даны в таблице [11, 13].

Значения гармонических постоянных приливов Черного моря (Н - в сантиметрах, у - в градусах)

Порт М2 S2 N2 К2 К1 Р1 01 Средняя амплитуда квадратурного прилива, см Средняя амплитуда сизигийного прилива, см Нк1 + Но1 НМ 2 + Н82

Н 3,5 1,9 0,5 0,8 0,9 0,6 0,4

Одесса 10,8 3,2 0,24

У 142 148 136 155 79 109 112

Севас- Н 0,4 0,3 0,1 0,1 0,3 0,4 0,1

тополь У 110 90 126 100 79 23 61 0,57

Н 2,9 1,2 0,7 0,5 1,5 0,35 0,9

Поти 8,2 3,4 0,58

У 293 299 280 29'/ 299 291 277

Соотношение амплитуд суточных и полусуточных гармонических составляющих (Нк + Н^ )/(НМ2 + Н3 ) позволяет определить характер

приливов. Так, для Одессы они являются полусуточными, поскольку (Нк + Н0 )/(Нм + Н3 ) < 0,25, для Севастополя и Поти

0,5 < (Нк + Н01 )/(Нм + Н5) < 1,5, а поэтому в этих пунктах приливы являются смешанными, неправильными полусуточными [14]. В целом средний период черноморских приливов составляет 15 ч 25 мин (рис. 2) [15].

Рис. 2. Прилив в порту Констанца (5 - 7 марта 1962 г., полнолуние) [15]

Современные методы математического прогноза циркуляции. Остановимся на краткой характеристике с-модели термогидродинамики, предложенной Институтом вычислительной математики (ИВМ) РАН и применимой

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

5

для описания как приливов, так и циркуляции вод в ограниченных бассейнах. В последние годы эта модель получила дальнейшее развитие. В частности, она использовалась для расчета циркуляции с применением вихредопускаю-щего [16] и вихреразрешающего [17] подходов.

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

Остановимся на основных уравнениях динамики моря в рамках предлагаемого подхода. Модель относится к классу а-моделей морской среды. В качестве вертикальной координаты используется безразмерная переменная а е [0, 1], задаваемая в случае учета свободной поверхности жидкости соотношением

а = (г-ОН Н-О,

(1)

где г - декартова вертикальная координата; H - глубина океана в состоянии покоя, предполагаемая ограниченной функцией с ограниченными производными; О- смещение уровня моря от невозмущенного положения.

Приведем постановку задачи о моделировании циркуляции океана в обобщенной системе координат. Переход от декартовой системы координат к обобщенной задается прямым и обратным дифференциалами преобразова-

ния:

ОУ ОХ

дх / дХ дх / дХ2 ду / дХ ду / дХ2 дг / дХ дг / дХ2

дх / дХ3 ду / дХ3 дг / дХ

ОХ ОУ

дХ / дх дХ / ду дХ / дг' дХ2 / дх дХ2 / ду дХ2 / дг дХ / дх дХ / ду дХ / дг

(2)

где X = (X, Х2, Х3) - декартовы координаты с единичной матрицей метрики = сИа§(1, 1, 1), а У = (х,у, г) - произвольные обобщенные координаты. При этом в каждой точке пространства можно построить систему локальных базисных векторов (1, }, к) = (дХ / дх, дХ / ду, дХ / дг), направленных вдоль соответствующих обобщенных координат.

Если (ОХ/ ОУ)г (ОХ / ОУ) - диагональная матрица, то система локальных базисных векторов (1, к) является ортогональной. Тогда система координат У = (х, у, г) называется ортогональной, и матрица метрики для нее имеет вид

От = (ОХ / ОУ)гО( Х)(ОХ / ОУ) = ^(гД г22, г32), а метрические коэффициенты г могут быть вычислены по формуле

'дХ дХ дХ

г ,=

д/ ' д/ ' дг

(/' = х, у, г).

(3)

(4)

6

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

В основе математической модели океана лежит система примитивных уравнений в приближениях гидростатики и Буссинеска, записанная в обобщенных ортогональных координатах по горизонтали и в а-системе координат по вертикали. Уравнения модели имеют вид:

О и - (I + %)уН = -

H

1 0 1 дра дО —Рх +--— + ё

р0 р0 дх дх

Н

(

О.у + (I + %)иН =--

- Р + *

Ро Ро ду

1 дРа О ду

д V ди ^

+---+ Fv,

да Н да

д V дv ^

+-----Fu,

да Н да

1 (диНГу дуНг}

г г

х у

дх

- + -

ду

+

да дО да д1

О, в = ±*-в + Ов+™,

да Н да да

О,Б =^ + ОБ, ' да Н да

р = р (в, Б + 35%о, р„ ) - р (0,0, ро ёаН).

(5)

(6)

(7)

(8)

(9) (10)

Здесь ра - атмосферное давление; и = (и,у) - вектор горизонтальной скорости, и и V - его зональная и меридиональная проекции; а - вертикальная скорость в а-системе координат, связанная с вертикальной скоростью ^ в г-системе координат соотношением

а = ^ -

^ и д2 ^ V д2 гх дх г ду д1

в - потенциальная температура; Я - поток проникающей солнечной радиации; - соленость за вычетом константы 35%о; р - отклонение плотности воды от среднего распределения плотности, зависящего только от давления столба жидкости р0ёг при средней плотности в океане ро = 1,025 г-см"3 на глубине

2^ (, 1

1 + -

с 1 - угло-

г = аН; параметр Кориолиса I = 2О ътф, где О = -

86400^ 365-24

вая скорость вращения Земли с учетом годового вращения вокруг Солнца, ф - гео-

1 (дгу дг.

графическая широта; % = ■

гг

х у

Л

-V--и

дх ду

- слагаемое, описывающее до-

полнительный перенос импульса в криволинейных координатах; V и Vв, V.,- -коэффициент вертикальной турбулентной вязкости и коэффициенты вертикальной турбулентной диффузии тепла и соли соответственно, которые в случае устойчиво стратифицированного вертикального профиля потенциальной плотности рассчитываются согласно параметризации Пакановского -Филандера или Монина - Обухова, а в случае неустойчиво стратифициро-

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

г

х

7

ванного - полагаются большими для параметризации конвекции. Нелинейное уравнение состояния р = р(в, £ + 35 %о, рк ) для расчета плотности воды, учитывающее сжимаемость за счет давления столба воды р^, взято из работы [18].

Компоненты горизонтального градиента давления Рх и Ру в уравнениях

(5) и (6) рассчитываются с использованием уравнения гидростатики в специальной форме:

Р =1

х 2 Ы |дх

Н

-а\ —р-нр

дх дх ,

1

Р = т Е

д

2 I ду

Н\\Р-с

др ~да

ёа

(

-а

дН др -р-Н——

ду ду

м

(11)

которая позволяет уменьшить погрешности при их разностных аппроксимациях в а-системе координат. Так как Рх = Ру = 0, для линейного по вертикали

профиля плотности имеем р= constx аН Использование уравнения состояния в виде (10) также позволяет уменьшить эти погрешности, поскольку заранее вычитается та часть нелинейного по глубине профиля плотности, которая не дает вклада в горизонтальный градиент давления.

Оператор переноса, входящий в состав полной производной компонентов скорости в (5) и (6), используется в полудивергентной симметризованной форме:

а р = 1 гк ^ +

21 Я Я

+ -

2гхгу

УНи 1р + £ (ГУНир)+ГН5 + ('ИУр)

ду ду

+

(12)

1 Г др дар 21 да да

где р- переменная и или переменная V .

В новой версии модели оператор переноса, входящий в состав полной производной скалярных полей в уравнениях (8) и (9), используется в дивергентной форме:

дкр 1 ар = —-+

д ГхГу

V (ГуНи Р) +-д (.Гх^Р) дх ду

+

дар ~~да '

(13)

где р - это в или £, а также при необходимости любые другие скалярные поля.

Оператор боковой диффузии В тепла и соли выбирается одинаковым для в и £ в (8) и (9) и записывается в универсальном виде:

8

МОРЖОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

а

0

1

Г

Dp =

1 8

+ -

rxry 8x 1 8

8x

8a

1 8

rxry 8У

r

KyHrx-

8p 8p

У V

8y

- — к

8a

rxry 8a 1 8

KXH — к\д-* — кд-*

rx x V 8x x 8a ,

+

rxry 8a

r

KyH-x- к r

8p 8p

y

—к

a

(14)

где р есть либо в, либо 5", а Кх (х, у, аИ) и Ку (х, у, аИ) - коэффициенты горизонтальной диффузии 2-го порядка по х и у, выбираемые как некоторые функции от пространственных координат. Переменные Кх и Ку задают одну или комбинацию нескольких функций, вдоль изоповерхностей которых происходит боковая диффузия. В частности, это могут быть а-, Z- или р-поверх-ности.

Оператор боковой вязкости ¥ в уравнениях (5) и (6) представляет собой комбинацию операторов 2-го и 4-го порядков:

Fp = Hdiv h (AgradA )p—H (div h (b1/2gradA f p,

(15)

где p - либо u, либо gradh и div h - двумерные операторы боковых градиента и дивергенции, действующие на поверхностях a = const. Величины A и B - диагональные тензоры второго порядка:

A =

0

0 Ay

\

(

B

Bx 0

0

By

\

(16)

где Ах = Ах (х, у), Ау = Ау (х, у), Вх = Вх (х, у), Ву = Ву (х, у) - коэффициенты вязкости для операторов 2-го и 4-го порядков вдоль х и у , задаваемые как некоторые функции пространственных координат. Оператор 4-го порядка, по сравнению с оператором 2-го порядка, более эффективно подавляет высокочастотные пространственные гармоники и менее искажает основное крупномасштабное решение.

Граничные условия модели. В качестве граничных условий на поверхности океана (а = 0) для скорости задаются поток импульса от напряжений трения ветра тх ,ту и универсальное условие для с :

К ^ (Tx ,Ty )

— — la=0 =-> ®la=0 = 0

H 8a p0

(17)

а для температуры и солености - нормированные потоки тепла ^ и соли д3 :

-

H 8a 0='

vL 8S_ H 8a

|a=0 = 4S •

(18)

Поток ^ рассчитывается с учетом потоков явного и скрытого тепла, длинноволновой и коротковолновой радиации и потока, вызванного наличи-

r

x

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

9

ем льда, а поток - с учетом баланса пресной воды, обусловленного осадками, испарением, стоком рек и образованием или таянием льда.

На дне (а = 1) задаются условие непротекания, имеющее в а-системе координат весьма простой вид

а|а=0 = 0, (19)

и условие квадратичного придонного трения

1а=0 = С^и2 + V2 + в2ь и |а=1, (20)

И да

где Св = 2,5 -10-3 и = 5 см-с-1 - эмпирические константы.

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

В а-модель циркуляции океана может быть включена модель динамики -термодинамики морского льда [19].

Особенности численной реализации а-модели циркуляции ИВМ РАН. Главная особенность этой модели, которая отличает ее от известных моделей океана, заключается в том, что при численной реализации в ней используется метод расщепления [20] по физическим процессам и пространственным координатам.

Для этого уравнения термогидродинамики океана записываются в специальной, симметризованной форме. Она позволяет представить оператор дифференциальной задачи в виде суммы более простых операторов, каждый из которых является неотрицательным в норме, определяемой законом сохранения полной энергии. Это дает возможность расщепить оператор полной задачи на ряд операторов более простых задач [20] и построить пространственные аппроксимации соответствующих групп слагаемых (в разных уравнениях) так, чтобы закону сохранения энергии, выполняющемуся для исходной дифференциальной задачи, удовлетворяли все расщепленные дискретные задачи. Конечно-разностные аппроксимации уравнений строятся на сетке С.

Метод расщепления позволяет эффективно реализовывать неявные схемы интегрирования по времени для уравнений переноса - диффузии. В задаче геострофического приспособления компоненты ускорения Кориолиса аппроксимированы неявно.

Модульная структура модели. Перед решением уравнений (5) - (10) в модели производятся следующие вспомогательные расчеты, результаты которых используются при решении основной системы уравнений.

Исходные атмосферные данные заданы в обычной географической системе координат с пространственно-временным разрешением, отличным от модельного, поэтому они переводятся на модельную область внутри расчетного блока модели путем пространственной и временной интерполяции.

Расчет характеристик морского льда опирается на локально-одномерную модель термодинамики [19], модель переноса [21] и модель динамики льда [22].

10

МОРЖОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

Расчет потоков тепла, соли и импульса в океан производится с использованием как проинтерполированных на модельную область атмосферных данных, так и рассчитанных характеристик морского льда, а также характеристик поверхности океана из решения задачи на данный момент времени, который считается известным.

Расщепление системы уравнений (5) - (9) проводится на нескольких иерархических уровнях. Сначала используется расщепление по физическим процессам. На более высоких уровнях процесс расщепления доходит до выделения простейших локально-одномерных по пространству уравнений. На каждом временном интервале интегрирования (/■ процесс, описываемый частично линеаризованной системой уравнений (5) - (9), представляется в виде суперпозиции процесса переноса - диффузии для в, Б, и, V и процесса приспособления (адаптации) полей скорости и плотности. В дифференциальной постановке эти задачи описываются следующими уравнениями (решение исходной задачи на момент времени £. считается известным):

переноса - диффузии в и Б -

да Н да да

переноса - диффузии и и V -

ЦБ = ^ + Ш; ' да Н да

да Н да да Н да

Процесс адаптации (приспособления) гидрологических полей описывается системой трех уравнений:

Ро Рх +

ди , 1 ( 1 дрп

--IV =--~ О , га

д г

дv , 1 --+ 1и =--

д г„

Р0 Ру +

Ро дх 1 дРа

g

дх

Л

К

Ро ду " дУ

л

g

д£ _ 1 (дтуиН дг^НЛ

дл

г г

х у V

дх

+ -

ду

+ -

дю ~да

Процесс переноса - диффузии реализован с помощью расщепления по физическим процессам: перенос, боковая диффузия и вертикальная диффузия. Для решения задачи переноса по времени используется явная схема Адамса - Башфорта. Дивергентная форма оператора переноса обеспечивает сохранение тепла и соли в океане в случае отсутствия потоков на границах. Задача по времени для боковой диффузии решается по явной, а для вертикальной - по неявной разностной схеме.

Для описания процесса переноса - диффузии применяется расщепление по элементарным процессам переноса - диффузии вдоль координат. Это позволяет сделать полудивергентная форма, обладающая при условии непротекания на границах свойством кососимметрии (неотрицательности) для каждого направления отдельно.

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

11

Процесс приспособления гидрологических полей решается в три этапа. Сначала по значениям в и 5", полученным на первом этапе, рассчитывается плотность и далее по уравнению (11) находятся компоненты градиента давления Рх и Ру. По рассчитанным Рх и Ру вычисляется обусловленное ими изменение импульса:

ды

Р, — = --1-Р.

х -и У

р0Гх д Р0Гу

Остальная часть системы решается путем разделения на баротропную и бароклинную моды:

1 1 ы = Ы + ы' , V = V + V' , Ы = | ыёа, V = | \йа.

0 0

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

Система уравнений бароклинной адаптации имеет вид

--IV = 0, — + 1ы = 0 .

дГ д

При решении этой системы используется неявная схема с методикой диаго-нализации пространственного оператора для кориолисовых членов, возникающего при применении сетки С.

Вертикальная скорость находится путем интегрирования по глубине уравнения неразрывности (7) по рассчитанным горизонтальным составляющим бароклинной скорости с учетом условий непротекания и свободного скольжения на боковых границах:

с =

а 1 ^

•> г г

1 х У V

( ды Иг ду'ИкЛ

дх

■ + ■

ду

dа.

Граничные условия для вертикальной скорости на поверхности и дне

1 1

удовлетворяются автоматически, поскольку ^ ы 'с^а V^а = 0 .

0 0

Для реализации процесса баротропной адаптации требуется совместное решение трех уравнений, записанных с использованием неявной схемы по времени:

ды

1

--IV = —

д Гх

дv 1

--+ 1ы = —

г

1 дЕа

дх р0 дх

\

д_£

Ы

1 дЕа

_У Р0 _У

1

( дгыИ _гШ Л

г г

х'у V

■ + -

дх ду ^

Решение этой задачи может производиться как прямыми, так и итераци-

онными методами. 12

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

Реализация ст-модели для описания циркуляции Черного моря. Рассматриваемая версия ст-модели с процедурой расщепления охватывает акваторию Черного моря. Она допускает учет бассейна Азовского моря. Расчетная область расположена от 27°26'60" до 41°45'00" в. д. и от 40°54'36" до 47°16'12" с. ш. с пространственным разрешением (0°3') х (0°2'24") по долготе и широте соответственно. Сеточная область в горизонтальной плоскости содержит 287 х 160 узлов. По вертикали задано 40 неравномерно распределенных по глубине ст-уровней. Шаг по времени, исходя из условий устойчивости решения, был выбран 300 с.

Батиметрия дна получена по данным Севастопольского отделения Государственного океанографического института (ГОИН) им. Н.Н. Зубова, которые представляют собой отдельно взятые карты топографии акватории Черного моря с разрешением (0°0'36,00") х (0°0'36,00") по долготе и широте и Азовского моря с пространственным разрешением соответственно (0°0'35,06") х (0°0'35,06"). Исходные данные высокого разрешения сглажены для устранения локальных особенностей, а затем проинтерполированы на модельную область. Далее модельная батиметрия на сетке с разрешением по долготе и широте (0° 3') х (0°2'24") была еще раз сглажена для устранения изломов. Боковые границы области задаются в виде вертикальной стенки с минимальной глубиной 5 м. Сглаживание топографии и ненулевая глубина во всех точках области, включая береговые, необходимы для с-модели океана, поскольку здесь используется преобразование вертикальной координаты (1). Батиметрия бассейна представлена на рис. 3.

28° 30° 32° 34° 36° 38° 40°

Рис. 3. Батиметрия (м) Черного моря

Для задания начальных условий на поверхности океана, а также для построения распределений начальных данных по глубинам, были использованы материалы, предоставленные Севастопольским отделением ГОИН. Они содержат данные среднемесячных климатических полей по температуре и солености Черного моря. По этим данным рассчитывались начальные поля температуры и солености, а также среднемесячные климатические поля температуры и солености на поверхности Черного моря.

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

13

Поскольку мы имеем дело с моделью общей циркуляции океана, в которой используется ст-система координат, то возникает необходимость интерполяции 3-мерных полей данных. При передаче в модель данные переводят на ст-уровни. Напомним, что ст-координатами (а = г / И) называется система координат, в которой поверхность моря принята в качестве нулевой глубины, а дну (-Н соответствует максимальная глубина (в каждой точке) со значением координаты, равным единице. Уровни в такой системе координат не совпадают со стандартными горизонтами, на которых обычно располагаются исходные сеточные данные по температуре и солености. Пересчет на ст-уровни осуществляется следующим образом. Вначале производится двумерная интерполяция на модельную сетку по описанному выше алгоритму. Затем по вертикали используется кусочно-линейное представление данных, после чего расчет для нужных глубин осуществляется с помощью простой процедуры одномерной линейной интерполяции.

Приведем примеры реализации численных экспериментов по расчету циркуляции Черного и Азовского морей. Для верификации модели была проведена серия численных экспериментов за период с 1 января 2007 г. по 31 декабря 2008 г. с различными параметрами модели. В качестве характеристик атмосферного воздействия использованы данные реанализа ЕЕЛ-1п(впт.

Оказалось, что результаты расчетов, особенно воспроизведение холодного промежуточного слоя в поле температуры, существенно зависят от коэффициентов вязкости и диффузии как горизонтальной, так и вертикальной. После серии экспериментов были выбраны следующие параметры: привязка к поверхностной температуре - 0,0; привязка к поверхностной солености -2,010-3 смс-1; коэффициент горизонтальной диффузии 2-го порядка -5,0-105 см2•с-1; коэффициент боковой вязкости 4-го порядка - 1,0-1017 см4 с-1.

Использована вертикальная параметризация Филандера - Пакановского: коэффициент вертикальной диффузии для температуры — 0,5 - 0,0005 см2-с-1; коэффициент вертикальной диффузии для солености — 0,1 - 0,0001 см2 с-1; коэффициент вертикальной вязкости — 10,0 - 1,0 см2 с4.

28° 30° 32° 34° 36° 38° 40° в. д 1.8 3.6 5.4 7.2 9 10.8 12.6 14.4 16.2 18 19.8 см/с

Рис. 4. Поле течений Черного моря на глубине 10 м, рассчитанное для 7 октября 2007 г. (стрелки показаны для каждой третьей точки расчетной сетки)

14

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

На рис. 4 представлено рассчитанное описанным выше методом поле течений Черного моря на глубине 10 м для 7 октября 2007 г. Хорошо выделяется Основное Черноморское течение. Вихревая структура циркуляции Черного моря занижена вследствие большого коэффициента вязкости 2-го порядка в данном эксперименте (106 см2 с-1).

На рис. 5 показано изменение среднемесячной циркуляции на глубине 5 м для июня, июля и августа 2007 г.

28° 30° 32° 34° 36° 38° 40° в. д.

Рис. 5. Среднемесячная циркуляция Черного и Азовского морей на глубине 5 м, рассчитанная для июня, июля и августа 2007 г. (стрелки показаны для каждой третьей точки расчетной сетки)

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

15

По результатам расчетов прослеживается достаточно сильная изменчивость отклонений уровня в течение года, видны «очки Книповича». Батум-ский антициклон воспроизводится недостаточно четко, что связано с грубым заданием пространственного распределения атмосферного форсинга по данным реанализа ЕЯЛ-Шегт.

Добавление в модель блока с приливообразующим потенциалом. В целом современные глобальные модели приливов можно разделить на три группы: гидродинамические, эмпирические и гидродинамические с ассимиляцией данных наблюдений (мареографов, спутниковой альтиметрии) [14, 23, 24]. Прогноз приливов в Черном море с использованием последних двух типов моделей затруднен в связи с отсутствием непрерывных рядов наблюдений. Существующие данные разрозненны, что не позволяет провести интерполяцию данных между калибровочными станциями и получить картину приливных составляющих.

Для описания генерации приливов используются три основных подхода: в уравнения добавляется приливообразующий потенциал [25]; на открытых границах задается специальное граничное условие, имитирующее вхождение приливной волны в расчетную область [26, 27]; применяются оба подхода одновременно [28 - 31]. Черное море - внутреннее море Атлантического океана, большей частью закрытое от сообщения с Мировым океаном. Проливом Босфор оно соединяется с Мраморным морем и далее через пролив Дарданеллы - с Эгейским и Средиземным морями, а через Керченский пролив - с Азовским морем. Средиземноморские приливные волны затухают в проливах, поэтому приливы в Черном море образуются лишь под действием при-ливообразующих сил. Для описания генерации приливов следует использовать подход, подразумевающий добавление в исходную систему примитивных уравнений (5) - (9) в приближениях гидростатики и Буссинеска, записанную в обобщенных ортогональных координатах по горизонтали и в ст-сис-теме координат по вертикали, приливообразующих потенциалов Луны Ом или (и) Солнца О, определяемых следующим образом:

где G - гравитационная постоянная, Мм и Ms - массы Луны и Солнца, Ге -радиус Земли, 1м и ls - расстояния от центра Земли до центров Луны и Солнца, Zm и Zs - зенитные углы Луны и Солнца.

Зенитный угол выражается через склонение (Луны, Солнца) ё, широту наблюдателя Ф и часовой угол а (угловое расстояние вдоль небесного экватора от меридиана наблюдателя до меридиана небесного тела) по формуле cosZ = cos£cos0cosa + sin д sin Ф. Часовой угол рассчитывается так: a = Á + 36CPí / T, где X - долгота наблюдателя, T - период, через который небесное тело возвращается на меридиан наблюдателя [25, 29, 32].

Все перечисленные переменные необходимо знать в каждый момент расчетного времени. Их значения были получены с помощью специального функционала военно-морской обсерватории США [33, 34] и подготовлены в виде массива в качестве исходных данных.

16 МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

Таким образом, уравнения (5), (6) описанной выше ст-модели совместной циркуляции Черного и Азовского морей принимают вид:

n п ел zr Н ( 1 d 1 дРа д^Л д v ди „ да Dtu - (l + %)vH =---Px +--+ g-2-

T

+---+ Fv--

да H да дх

Dtv + (l + %)uH = -

H

r

y

VPo Po дх дх

( i 1 дра д£ л —Py +--— + g—

Po Po ду ду I да H да ду

д v дv „ да +-----Fu -

T

где Qt - гравитационный потенциал Луны (Ом) и Солнца (Os) или их сумма (Ом + Os). В отличие от исходных уравнений они приспособлены для расчета приливов.

Заключение. Таким образом, обобщены литературные материалы по приливным колебаниям гидрофизических полей в Черном море. Проведен анализ параметров приливов в Черном море по натурным данным для различных участков морского побережья. Дана характеристика основных параметров приливных колебаний в регионе, представлены примеры котидальной карты и карты изоамплитуд для полусуточного прилива.

Для расчета циркуляции вод в Черном и Азовском морях была использована ст-модель общей циркуляции океана, разработанная в ИВМ РАН. Представляемая версия имеет пространственное разрешение по долготе и широте около 4 км. По вертикали задавалось 40 неравномерно распределенных ст-уровней. Шаг по времени для обеспечения устойчивости решения принят 300 с.

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

Анализ существующих моделей описания приливных колебаний позволил выбрать наиболее оптимальный подход к моделированию приливов и сопутствующих им течений в бассейне Черного моря. Для описания генерации приливов в исходную ст-модель циркуляции Черного и Азовского морей был добавлен приливообразующий потенциал. Получены необходимые выражения для определения приливообразующих потенциалов Луны и Солнца.

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 14-45-01025 р_юг_а.

СПИСОК ЛИТЕРАТУРЫ

1. Каган Б.А. Гидродинамические модели приливных движений в море. - Л.: Гидрометео-издат, 1968. - 219 с.

2. Марчук Г.И., Каган Б.А. Динамика океанских приливов. - Л.: Гидрометеоиздат, 1991. -471 с.

3. Андросов А.А., Вольцингер Н.Е. Проливы Мирового океана. Общий подход к моделированию. - М.: Наука, 2005. - 187 с.

r

х

МОРСКОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

17

4. Горячкин Ю.Н., Иванов В.А. Уровень Черного моря: прошлое, настоящее и будущее. -Севастополь: МГИ НАН Украины, 2006. - 210 с.

5. Доценко С.Ф., Иванов В.А. Природные катастрофы Азово-Черноморского региона. -Севастополь: ЭКОСИ-Гидрофизика, 2010. - 174 с.

6. Marinescu A., Sclarin O. Les variations periodignes du nivean de la Mer Noire a longstanga // Trans. Museum histoire nature Gr. Antipa. - 1968. - № 1. - P. 531 - 535.

7. Endros A. Die Seiches des Schwarzen und Azowschen meers und die dortigen Hubhohen der Gezeiten // Ann. Hyd. Mar.Met. - 1932. - Bd. 60, Ht. 11. - S. 442 - 453.

8. Иванов В.А., Манилюк Ю.В., Черкесов Л.В. О сейшах Черного моря // Метеорология и гидрология. - 1996. - № 11. - С. 57 - 63.

9. Архипкин В.С., Иванов В.А., Николаенко Е.Г. Моделирование собственных колебаний в южных морях // Численное моделирование гидрофизических процессов и явлений в замкнутых водоемах / Под ред. А.С. Саркисяна. - М.: Наука, 1987. - С. 78 - 91.

10. Горячкин Ю.Н., Иванов В.А., Репетин Л.Н. , Хмара Т.В. Сейши в Севастопольской бухте // Тр. УкрНИГМИ. - 2002. - Вып. 250. - С. 342 - 353.

11. Defant A. Physical Oceanography. V. 2. - Oxford, New York: Pergamon Press, 1961. -598 p.

12. Ястреб Я.П., Хмара Т.В. Сравнительная характеристика приливных волновых движений в морях средиземноморского типа // Экологическая безопасность прибрежной и шельфовой зон и комплексное использование ресурсов шельфа. - Севастополь: ЭКОСИ-Гидрофизика, 2005. - Вып. 12. - С. 60 - 69.

13. Yuce H. Analysis of the water level variations in the eastern Black Sea // J. Coast. Res. -1993. - 9, № 4. - P. 1075 - 1082.

14. Boon J. Secrets of the tide: tide and tidal current analysis and applications, storm surges and sea level trends. - Woodhead Publishing, 2004. - 224 p.

15. Le Provost C.F., Lyard F., Molines J.M. et al. A hydrodynamic ocean tide model improved by assimilating a satellite altimeter-derived data set // J. Geophys. Res. - 1998. - 103, № C3. - P. 5513 - 5529.

16. Мошонкин С.Н. , Дианский Н.А., Гусев А.В. Влияние взаимодействия Атлантики с Северным Ледовитым океаном на Гольфстрим // Океанология. - 2007. - 47, № 2. -C. 197 - 210.

17. Дианский Н.А., Залесный В.Б., Мошонкин С.Н. и др. Моделирование муссонной циркуляции Индийского океана с высоким пространственным разрешением // Океанология. -2006. - 46, № 4. - C. 421 - 442.

18. Marchuk G.I., Rusakov A.S., Zalesny V.B. et al. Splitting numerical technique with application to the high resolution simulation of the Indian Ocean circulation // Pure Appl. Geophys. -2005. - 162. - P. 1407 - 1429.

19. Parkinson C.L., Washington W.M. A large scale numerical model of sea ice // J. Geophys. Res. - 1979. - 84, № C1. - P. 311 - 337.

20. Яковлев Н.Г. Воспроизведение крупномасштабного состояния вод и морского льда Северного Ледовитого океана в 1948 - 2002 гг. Часть I. Численная модель и среднее состояние // Изв. РАН. Физика атмосферы и океана. - 2009. - 45, № 3. - С. 383 - 398.

21. МарчукГ.И. Методы вычислительной математики. - СПб.: Лань, 2009. - 608 с.

22. Briegleb B.P., Hunke E.C., Bitz C.M. et al. The sea ice simulation of the Community Climate System Model, version two. - NCAR Tech. Note NCAR/TN-45+STR, 2004. - 34 p.

23. Hunke E.C., Dukowicz J.K. An elastic-viscous-plastic model for sea ice dynamics // J. Phys. Oceanogr. - 1997. - 27, № 9. - P. 1849 - 1867.

24. Arabelos D.N., Papazachariou D.Z., Contadakis M.E. et al. A new tide model for Mediterranean Sea based on altimetry and tide gauge assimilation // Ocean Sci. - 2011. - 7, № 3. -P. 429 - 444.

18

МОРЖОЙ ГИДРОФИЗИЧЕСКИЙ ЖУРНАЛ № 1 2016

25. Tierney C.C., Kantha L.H., Born G.H. Shallow and deep water global ocean tides from altim-etry and numerical modeling // J. Geophys. Res. - 2000. - 105, № C5. - P. 11259 - 11277.

26. Feistel R., Nausch G., WasmundN. State and evolution of the Baltic Sea, 1952 - 2005: a detailed 50-year survey of meteorology and climate, physics, chemistry, biology, and marine environment. - Hoboken: John Wiley & Sons, 2008. - 703 p.

27. MzadekR. Hydrodynamic tidal model of Cook Strait. - Ecole MATMECA, 2005. - 30 p.

28. Kowalik Z., Untersteiner N. A study of the M2 tide in the Arctic Ocean // Dt. Hydrogr. Z. -1978. - 31. - S. 216 - 229.

29. Tsimplis M.N., Proctor R., Flather R.A. A two dimensional tidal model for Mediterranean Sea // J. Geophys. Res. - 1995. - 100, № C8. - P. 16223 - 16239.

30. Einspigel D. Barotropnl oceanicky slapovy model. - Univerzita Karlova v Praze, 2012. -55 p.

31. Luettich R.A., Westerink Jr., Westerink J.J. Continental shelf scale convergence studies with a barotropic tidal model. Quantitative skill assessment for Coastal Ocean Models // Amer. Geophys. Union press. - 1995. - 48. - P. 349 - 371.

32. Kowalik Z., Proshutinsky A.Y. Diurnal tides in the Arctic Ocean // J. Geophys. Res. - 1993. -98, № C9. - P. 16449 - 16468.

33. Kowalik Z., Luick J. The Oceanography of Tides. - Fairbanks, 2013. - 157 p.

34. Kaplan G., Bangert J., Bartlett J. et al. User's Guide to NOVAS 3.0. - Washington: U.S. Naval Observatory, 2009. - 124 p.

Block approach to simulation of circulation and tides in the Black Sea

S. F. Dotsenko*, V. B. Zalesny**, N. K. V. Sannikova*

*Marine Hydrophysical Institute, Russian Academy of Sciences, Sevastopol, Russia e-mail: [email protected] **Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow, Russia e-mail:[email protected]

Circulation and tides in the Black Sea are analyzed. The basic characteristics of tidal oscillations are assessed. Besides, the Black Sea water circulation is simulated using the block approach. The ff-mo-del of the ocean circulation developed in the Institute of Numerical Mathematics of RAS is applied. The model spatial resolution over the longitude and latitude is about 4 km. 40 irregularly distributed ff-levels are preset over the vertical; the step in time is 300 s.

The vortex structure is distinctly manifested in the Black Sea circulation. The Rim Current which characterizes general cyclonic circulation along the Black Sea perimeter forming two evident vortices is reproduced. To describe generation of tides, the module corresponding to the Sun and Moon tide-generating potentials is introduced into the ff-model.

Keywords: Black Sea, general circulation, tides, mathematical model, block approach, computational experiments.

MOPCKOH rHflPOOH3HHECKHH ®yPHA.n № 1 2016 19