Sensitivity of the ocean circulation model to the k–omega vertical turbulence parametrization

Full Text

1. ВВЕДЕНИЕ

В настоящее время интенсивно развиваются численные модели общей циркуляции океана (МОЦО), основанные на примитивных уравнениях. Модели применяются для изучения изменчивости климатической структуры гидрофизических полей и для их гидродинамического прогноза. Кардинальными направлениями повышения адекватности моделей являются три: увеличение пространственного разрешения, улучшение физических параметризаций и ассимиляция данных наблюдений [1–4]. В этой работе мы описываем развитие МОЦО ИВМ РАН по пути усовершенствования параметризации вертикального турбулентного обмена.

Модель представлена в эволюционной форме, а алгоритм ее решения основан на методе многокомпонентного расщепления [1, 2]. Основным является расщепление по физическим процессам — это внешний контур расщепления модели [1, 2]. Решение новой модельной параметризации также проводится в рамках метода расщепления. Встраиваемая в МОЦО новая подсистема записывается в эволюционной форме, а ее оператор представляется в виде суммы подоператоров более простой структуры [1, 2]. Этому условию удовлетворяет рассматриваемая в работе k–ω модель вертикальной турбулентности.

Важным фактором развития МОЦО является улучшение физического описания вертикального турбулентного обмена импульсом, теплом и солью. По сравнению с крупномасштабной циркуляцией пространственно-временные масштабы вертикального турбулентного обмена значительно меньше: порядка 1–10 км по горизонтали, 1–10 м по вертикали и от минут до нескольких часов по времени. В рамках традиционной модели, основанной на «примитивных» уравнениях, это — подсеточный процесс, подлежащий параметризации. Параметризация подсеточных процессов остается одной из актуальнейших задач моделирования и прогноза гидрофизических и метеорологических полей [4–14].

Часто турбулентное перемешивание описывают в МОЦО оператором второго порядка, с коэффициентами турбулентного обмена импульсом ν_u, теплом ν_T и солью ν_S [1–5]. Для определения коэффициентов обмена используются модели, основанные на двух уравнениях [7]. Первое уравнение описывает поведение кинетической энергии турбулентности k (КЭТ). Второе уравнение записывается либо для масштаба турбулентности l (k – kl-замыкание), либо для скорости удельной диссипации КЭТ ε (k – ε-замыкание), либо для частоты ω (k – ω-замыкание). Характеристики l, ε, ω связаны алгебраическими соотношениями. Эти модели описывают развитую турбулентность и разделяют слои с развитой и слаборазвитой перемежающейся турбулентностью. Перемешивание в слоях слаборазвитой турбулентности осуществляется за счет двойной диффузии, разрушения внутренних и приливных волн и других эффектов [8, 9].

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

2. ОСОБЕННОСТИ k – ω МОДЕЛИ И МЕТОД РЕШЕНИЯ УРАВНЕНИЙ ТУРБУЛЕНТНОГО ОБМЕНА

Сравним особенности часто используемой для параметризации вертикального турбулентного обмена k – ε модели [7, 9], с особенностями k – ω модели. Первая основана на системе двух уравнений для кинетической энергии турбулентности (КЭТ) k и ее диссипации ε. Вторая — на уравнениях для КЭТ и ее частоты диссипации ω. Запишем уравнения k – ε и k – ω моделей в σ-системе координат [7] (предполагая k ≠ 0, ε ≠ 0):

 (1)

 (2)

Здесь   — геопотенциальная вертикальная координата, уровень и глубина океана; G, N — частоты сдвига скорости и плавучести, рассчитываемые в МОЦО (N ² > 0 соответствует устойчивой стратификации)

,

;

 – оператор горизонтальной диффузии, ν_ρ — коэффициент вертикального турбулентного обмена потенциальной плотностью; σ_k, σ_ε, σ_ω — безразмерные турбулентные числа Шмидта для k, ε и ω.

Для функций и параметров моделей турбулентности выполняются соотношения:

 (3)

где l — характерный масштаб турбулентных возмущений,  — безразмерные функции устойчивости для вектора и скаляра. Значения параметров  можно найти в [7]:

 (4)

К (1)–(2) присоединяются соответствующие граничные и начальные условия [13].

Учитывая (3), системы (1), (2) можно переписать в виде

 (5)

 (6)

Видно, что в отличие от k – ε, система k – ω уравнений (6) имеет более простой вид, поскольку уравнение для ω не зависит от k!

Применим к (5), (6) метод расщепления по физическим процессам. Процедуру решения уравнений расщепим на два этапа, описывающих процессы переноса–диффузии и генерации–диссипации. На каждом этапе по времени решаются более простые расщепленные подсистемы, полученное на текущем этапе решение используется в качестве начального условия на последующем этапе.

Процессы переноса–диффузии и генерации–диссипации имеют разные характерные времена. Перенос–диффузия — медленная трехмерная эволюция полей k, ε, ω, аналогичная эволюции фоновых полей МОЦО. Генерация–диссипация — быстрый локальный по пространству процесс, описывающий динамику турбулентных возмущений.

На этапе трехмерного переноса–диффузии имеем:

 (7)

На верхней границе океана, при σ = 0, поставим условия:

 (8)

где  — скорость трения в воде у поверхности океана, c_g — безразмерный параметр, зависящий от ветра и волн (параметр ветровой генерации):  [11] или  [12], ,  — «потоки диссипации и частоты диссипации» на поверхности. Об одном способе задания «потока диссипации» дана информация на с. 22–23 работы [18]. В экспериментах далее мы полагали  = 0.

На дне океана, при σ = 1, поставим условия отсутствия потоков k, ε и ω по нормали.

На этапе генерации–диссипации имеем уравнения:

(9)

 (10)

 (11)

 (12)

3. АНАЛИТИЧЕСКОЕ РЕШЕНИЕ k – ω МОДЕЛИ НА ЭТАПЕ ГЕНЕРАЦИИ–ДИССИПАЦИИ

Перепишем (11), (12) в виде

 (13)

 (14)

 (15)

(16)

 (17)

Отметим, что B > 0, C > 0, D > 0, в то время как A может менять знак. Предположим, что A, B, C, D не зависят от времени на временном шаге модели. Тогда, поскольку уравнение для ω не зависит от k, его несложно решить аналитически, а затем найти из (13) k. Аналитическое решение (14), (13) имеет вид:

или

 (18)

 (19)

Здесь

ω⁰ и k⁰ — значения в начальный момент времени на этапе генерации–диссипации.

Отметим, что для коэффициентов A, B с учетом (3) можно использовать две формы записи (см. (15), (16)):

 (20)

или

 (21)

Алгоритм, использующий аналитическое решение (18), (19), мы называем «алгоритмом расщепления турбулентности» (АРТ) [10, 13, 14].

Важными факторами повышения адекватности моделирования крупномасштабной циркуляции океана является подбор модельных параметров и использование в расчетах данных наблюдений. В наших работах [13, 14] изучена возможность повышения адекватности МОЦО с помощью учета в A и B данных наблюдений о среднегодовой климатической частоте плавучести. С этой точки зрения алгоритм является достаточно гибким. Во-первых, он позволяет варьировать A и B, изменяя функции устойчивости в (20) или число Прандтля в (21), сохраняя простоту и устойчивость вычислений. Во-вторых, — решать уравнения турбулентности с большим шагом по времени, равным шагу МОЦО. В-третьих, — алгоритм можно эффективно использовать при решении задач четырехмерной вариационной ассимиляции данных наблюдений [15, 16].

Замечание 1. Используя соотношение   можно выписать аналитическое выражения для ν_u. Из (13), (14) следует уравнение для ν_u

 (22)

 (23)

Уравнение (22), с точностью до коэффициентов (23), совпадает с (13), и его решение имеет вид

. (24)

Замечание 2. Из вида решений (18), (19), (24) следует, что функции ω, k и ν_u на этапе генерации–диссипации остаются положительными. Аналитические решения для ω, k и ν_u позволяют оценить их характерную изменчивость по времени. Например, переписывая (18), имеем:

При  функция  а при  ее значение будет  Отсюда характерную изменчивость ω, или время стационирования решения, можно оценить как

Наши численные эксперименты показывают, что характерное время процессов вертикальной турбулентности в верхнем перемешанном слое Северной Атлантики в зимний период изменяется от 2–3 ч в тропической зоне до 1–2 мин в средних широтах.

4. ПАРАМЕТРЫ МОДЕЛИ ОБЩЕЙ ЦИРКУЛЯЦИИ ОКЕАНА И k – ω ПАРАМЕТРИЗАЦИИ

Модель общей циркуляции океана ИВМ РАН формулируется в сигма-системе координат со свободной поверхностью [1, 2, 15]. Она основана на примитивных уравнениях, записанных в приближениях Буссинеска и гидростатики. Уравнения записываются в биполярной ортогональной системе координат на сфере, полюсы расположены на географическом экваторе в точках 120° з.д. и 60° в.д. [14]. Прогностическими переменными модели являются горизонтальные скорости течения, высота уровня поверхности океана  потенциальная температура и соленость. Морской лед рассчитывается по модели [17].

Система уравнений МОЦО расщепляется на две подсистемы: перенос–диффузия субстанций и адаптация полей течений и плотности. При решении подсистем переноса–диффузии и адаптации повторно применяется метод расщепления по физическим процессам и по отдельным пространственным координатам. Разностная аппроксимация модельных уравнений по пространству осуществляется на сетке «C». Процессы вертикального турбулентного обмена описываются в рамках «диффузионного» подхода, коэффициенты обмена рассчитываются с помощью встроенной k – ω модели. Более подробное описание МОЦО приводится в [2, 15].

Опишем основные модельные параметры. При вычислении коэффициентов вертикальной турбулентной вязкости и диффузии в МОЦО используются соотношения подобия:

 (25)

В слоях развитой турбулентности, где k > k_min = 0.03 cм²/с² полагается ν_S = ν_T [18]. В слоях слабой турбулентности используются фоновые значения ν_u = 1, ν_S = ν_T = 0.05 cм²/с.

В численных экспериментах для МОЦО шаг по времени τ_ост равняется 1 ч. На этапе переноса–диффузии уравнения k – ω модели также решаются с шагом по времени 1 ч. На этапе генерации–диссипации уравнения решаются с шагом τ_T, который менялся в пределах 1 мин ≤ τ_T ≤ 1 ч.

Для задания коэффициентов A, B используются два подхода. В первом используем представление (21) и полагаем, следуя [18], в каждом слое:

 (26)

B (21) также полагаем:

 (27)

В этом случае, объединяя (21) с (26), (27), имеем

 (28)

Коэффициенты турбулентного обмена, с учетом (3) и (26), имеют вид:

 (29)

где число Прандтля, согласно [19], есть функция от числа Ричардсона Ri = N ²/G ² в виде:

 (30)

Во втором подходе выбираем A, B в форме (20). В этом случае, — в предположении о локальном равновесии напряжений Рейнольдса и турбулентных потоков тепла [18, 20], — функции устойчивости в каждом слое определяются в виде:

 (31)

где

Система (31) разрешима, поскольку ее определитель d ≠ 0. В качестве ограничений на развитие свободной конвекции и вырождение турбулентности можно использовать ограничения [18]:

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Цель численных экспериментов — изучить чувствительность решения МОЦО с встроенной в нее k – ω моделью к вариациям коэффициентов A и B. Расчеты проводятся для акватории, включающей Атлантический океан к северу от 30° ю.ш., Северный Ледовитый океан и Берингово море. Открытые границы области расположены на 30° ю.ш. и в проливах Алеутских островов. Область включает также Средиземное, Черное и Балтийское моря. Шаг сетки по широте и долготе равен 0.25°. По вертикали задаются сорок сигма-уровней со сгущением к поверхности океана. Рельеф дна океана сглажен в соответствии с горизонтальным разрешением модели так, чтобы отсутствовали резкие градиенты дна. Модельная глубина ограничена минимальной величиной 10 м.

Граничные условия на поверхности океана рассчитываются с использованием атмосферных характеристик по данным CORE-2 [21]. Потоки явного и скрытого тепла, влаги и напряжение ветра рассчитываются с дискретностью в один час по данным CORE-2 о температуре воздуха, влажности, компонентах скорости ветра и давлении на уровне моря с использованием модельной температуры воды. Потоки длинноволновой и коротковолновой радиации задаются с дискретностью одни сутки [21]. Для коротковолновой радиации учитывается ее проникающая способность. Интенсивность атмосферных осадков и речного стока задается с дискретностью один месяц. Сток основных рек учитывается косвенно в граничном условии для солености.

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

В качестве начальных условий выбираются климатические январские поля температуры и солености океана, отсутствие движения и морского льда [14]. В граничном условии (8) положено

Численные эксперименты проведены на период 20 лет с различными входными параметрами, указанными в табл. 1.

Таблица 1. Список вычислительных экспериментов

Эксперимент EP40. Коэффициенты A и B выбираются в форме (28). Параметр ветровой генерации в (8) c_g = 40. Коэффициенты турбулентной диффузии и вязкости рассчитываются по (29), где число Прандтля есть функция числа Ричардсона (30). Уравнения турбулентности решаются аналитически согласно (18), (19) с шагом по времени τ_T = 1 ч. Заметим, что шаг модели турбулентности по времени τ_T равен шагу МОЦО: τ_T = τ_ост.

Эксперимент EC40. Коэффициенты A и B выбираются в форме (20) с функциями устойчивости в форме (31), коэффициенты вертикальной турбулентной вязкости и диффузии в форме (25), параметр c_g = 40 в (8). Уравнения турбулентности решаются аналитически согласно (18), (19) с шагом по времени τ_T = 5 мин.

Эксперимент EC10. Он отличается от EC40 только тем, что в граничном условии (8) c_g = 10.

Эксперимент EС1h. В отличие от EC40 уравнения турбулентности решаются аналитически с шагом по времени равным шагу МОЦО: τ_T = τ_ост = 1 час. Функции устойчивости (31) на расчетном интервале t ^j < t < t ^{j + 1} берутся с предыдущего шага по времени t = t ^j. Функции устойчивости пересчитываются далее при t = t ^{j + 1}.

6. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Рассмотрим вопрос о чувствительности толщины верхнего квазиоднородного слоя (ВКС) к вариациям коэффициентов аналитического решения (18), (19). Сравнение проведем для периода максимального развития свободной конвекции. На рис. 1 показана толщина ВКС для Северной Атлантики и Арктики в феврале 20-го года модельного расчета, полученная в экспериментах EP40 и EС1h. Отметим, что данные распределения типичны для всех зимних периодов. Видна высокая чувствительность толщины ВКС к изменению коэффициентов A и B. Разница между двумя расчетами составляет десятки метров в средних широтах и сотни, — а в отдельных районах интенсивной конвекции до 1–2 тысяч, — метров в Лабрадорском, Норвежском и Гренландском морях. В эксперименте EP40 (рис. 1a) в Лабрадорском и Норвежском морях толщина ВКС заметно больше по сравнению с данными наблюдений [22]. В экспериментах EС1h и EС40, в которых толщины ВПС близки друг другу, этот недостаток устраняется (рис. 1б). В них более реалистично воспроизводится также толщина ВКС в восточной Атлантике. В целом, результаты EС1h и EС40 лучше соответствуют данным наблюдений [22], чем результаты EP40.

Рис. 1. Толщина верхнего квазиоднородного слоя в Северной Атлантике и Арктике, в феврале 20-го года модельного расчета в экспериментах EP40 (а) и EС1h (б). В пределах ВКС потенциальная плотность воды отличается от плотности на поверхности океана менее чем на 0.15 кг/м³. Координаты модельные (см. текст).

На рис. 2 показаны типичные для летнего периода профили температуры Т(z) и солености S(z) на Северном географическом полюсе при различных значениях параметров. Для сравнения приведены соответствующие климатические профили по данным Атласа океанов [23, 24]. Устойчивость стратификации верхнего слоя океана здесь в основном связана с распределением солености по вертикали. Летний период выбран для иллюстрации наибольшей чувствительности вертикальной структуры расчетных полей к вариациям параметра c_g. В этот период зона развитой турбулентности у поверхности океана относительно невелика. Уменьшение в четыре раза притока КЭТ в океан у поверхности в эксперименте EС10 по сравнению с EС40 уменьшает перемешивание. Опресненные (за счет таяния льдов и речного стока) воды проникают на меньшую глубину, и в верхнем 15-метровом слое в EС10 воспроизводится соленость на 0.6–0.9 PSU ниже, чем в EС40. Ниже 15 м, в галоклине, из-за недостатка КЭТ и слабого перемешивания (с верхними опресненными слоями), соленость в EС10 выше, чем в EС40 (рис. 2а). Указанные причины обуславливают наличие в EС10 более теплой воды в верхнем 25-метровом слое, в июне, по сравнению с EС40 (рис. 2б). В целом S(z) и T(z) в эксперименте EС40 ближе к климатическим данным наблюдений, чем в EС10.

Рис. 2. (a) Соленость (в PSU), (б) температура (в °С) на Северном полюсе. Июль, слой 0–50 м. Эксперимент EС40 – черная линия и полые кружки; EС1h — черная линия, штрихи; EС10 — красная линия, кружки; EP40 – зеленая линия, квадраты; соответствующие климатические данные для июля из Атласа океанов — синяя линия, крестики [23, 24].

При использовании функций устойчивости в виде (31) модельное решение лучше соответствует данным наблюдений, чем при задании числа Прандтля как функции от числа Ричардсона в эксперимент EP40 (рис. 2). В эксперименте EP40, как и в EС10, приповерхностный слой более пресный, чем в EС40 и EС1h. При этом в EP40 в 20-метровом слое заметно завышаются градиент солености и температура воды, и возникает завышенная по сравнению с наблюдениями термическая инверсия (рис. 2б). В целом можно сказать, что в EP40 уменьшен приток КЭТ. Распределение КЭТ в слое 0–20 м показывает, что в эксперименте EС40 величина КЭТ равна ~30 см²/с², а в EP40 ~5–15 см²/с².

Следует особо отметить малую чувствительность вертикальной структуры верхнего слоя океана при переходе от эксперимента EС40 к EС1h (рис. 2).

Эксперименты показывают заметное изменение полей температуры, солености и плотности в верхнем слое океана при вариации коэффициентов A и B (18), (19), что отражается и в поле течений. На рис. 3а показано среднее за 20 февралей поле течений в слое 10–300 м в Субполярном циклоническом круговороте в эксперименте EС1h. Хорошо выражены все известные течения региона. Пространственное распределение течений в эксперименте EP40 имеет ту же структуру, что и в EС1h. Чтобы показать различия в скоростях между этими экспериментами, на рис. 3б приведена их разность: «EС1h минус EP40». Относительные изменения модулей течений в отдельных районах могут достигать величины 20–30%. Так, на южной периферии Северо-Атлантического течения скорости в EС1h меньше, чем в EP40. Скорости в EС1h возрастают по сравнению с EP40 на северной периферии Северо-Атлантического течения и на западной периферии Лабрадорского течения (рис. 3). В целом в экспериментах EС40 и EС1h Северо-Атлантического течение приобретает более выраженный струйный характер, чем в EP40.

Рис. 3. (a) Линии тока и модуль скоростей течений (см/c, цветной фон) в слое 10–300 м. Среднее за 20 февралей (1948–1967 гг.): (а) — эксперимент EС1h, (б) — векторная разность течений «EС1h минус EP40». Координаты модельные: вверху — очертания Гренландии и Исландии, слева — Лабрадор и Ньюфаундленд.

7. ВЫВОДЫ

Описывается развитие модели общей циркуляции океана ИВМ РАН, используемой в различных исследованиях динамики морей и океанов [1–4, 15, 16]. Модель циркуляции включает k – ω параметризацию вертикального турбулентного обмена. Уравнения турбулентности расщепляются на два этапа, описывающие процессы трехмерного переноса–диффузии и локальной генерации–диссипации k и ω. На этапе генерации–диссипации уравнения решаются аналитически, что обуславливает эффективность и устойчивость расчетов. С вычислительной точки зрения алгоритм является достаточно гибким. Он позволяет варьировать коэффициенты аналитического решения, изменяя функции устойчивости или число Прандтля и сохраняя при этом простоту и устойчивость вычислений.

Изучена чувствительность модельных крупномасштабных полей Северной Атлантики — Северного Ледовитого океана к изменению коэффициентов турбулентного обмена. Изменение коэффициентов обмена производится путем изменения коэффициентов аналитического решения k – ω модели. Эксперименты показывают, что, варьируя коэффициенты аналитического решения на стадии генерации–диссипации k – ω модели, можно повышать адекватность моделирования. Приводятся оценки чувствительности совместной модели циркуляции и турбулентности к вариациям притока кинетической энергии турбулентности на поверхности океана. Даны также оценки времени выхода на стационар турбулентных процессов, показавшие принципиальную необходимость использования нестационарных уравнений турбулентных характеристик.

Источники финансирования. Работа выполнена в ИВМ РАН при поддержке Российского научного фонда (грант № 18-11-00163) и РФФИ (грант 18-05-00177).

About the authors

V. B. Zalesny

Author for correspondence.
Email: atarexm@himki.net
Russian Federation, Gubkina str, 8, Moscow, 119333

S. N. Moshonkin

Email: atarexm@himki.net
Russian Federation, Gubkina str, 8, Moscow, 119333

References

Supplementary files