Журнал вычислительной математики и математической физики, 2021, T. 61, № 10, стр. 1693-1703

1. ВВЕДЕНИЕ

Исследования по газовой динамике, инициированные и поддерживаемые Ю.Д. Шмыглевским, были ценны не только сами по себе, но и как тематика для подготовки научных кадров, в частности, учеником Юрия Дмитриевича был В.В. Александров, которому полученные компетенции позволили применить их в области моделирования глобальных процессов в атмосфере и климата. Развитие этих исследований представлено и в настоящей работе.

При исследовании долгосрочных изменений климата требуется рассматривать всю атмосферу, глобальный океан (с морским льдом) и верхний слой суши (почва и растительность) как взаимодействующие части единой системы, называемой климатической системой (см. [1]). Она имеет глобальный характер со значительно отличающимися временными и пространственными масштабами. Математическое моделирование является мощным инструментом для исследования климатической системы и прогнозирования (см. [2]).

В настоящей работе представлена глобальная трехмерная гидродинамическая модель климата, объединяющая модель общей циркуляции атмосферы (ОЦА), модель термохалинной циркуляции океана и модель эволюции морского льда. До этого модель океана использовалась с достаточно сильно агрегированной энерго-влаго-балансовой моделью атмосферы, описывающей характеристики приземного слоя (см. [3]). Модель ОЦА значительно более сложная и позволяет более точно описывать процессы в атмосфере. Функционирование совместной климатической модели рассматривается в режиме сезонного хода солнечного излучения.

Рассмотрена процедура проведения совместных расчетов модели океана и модели ОЦА. Для этого требуется обеспечить синхронизацию ряда параметров и процедур в обеих моделях. В связи с этим применяется процедура двумерной интерполяции данных, заданных на сетке модели океана и данных модели атмосферы и обратно. Особенностью этой процедуры является несовпадение конфигураций материков в моделях. В работе приведены результаты численных экспериментов с климатической моделью по воспроизведению климата и его чувствительности к увеличению концентрации углекислого газа для двух сценариев.

2. МОДЕЛЬ ТЕРМОХАЛИННОЙ ЦИРКУЛЯЦИИ ОКЕАНА

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

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

Учитывая эти замечания, система уравнений модели океана записывается в геострофическом приближении с фрикционным членом в уравнениях импульса по горизонтали (см. [3], [7], [8]). Значения температуры T и солености S подчиняются адвекционно-диффузионным уравнениям, что дает возможность описать крупномасштабную термохалинную циркуляцию океана. Приближенным образом учитываются также конвективные процессы (см. [7]).

Таким образом, система основных уравнений, для наглядности записанных в локальных декартовых координатах (x, y, z), где x, y – горизонтальные координаты и z – высота, направленная вверх, имеет следующий вид:

уравнения импульса по горизонтали

уравнение неразрывности

уравнение гидростатики

уравнение состояния морской воды

уравнение переноса и диффузии субстанции X (температуры и солености)

где u, $v$, w – компоненты вектора скорости течений, λ – переменный в пространстве фрикционный член, увеличивающийся к материкам и экватору; T, S, p – температура, соленость, давление соответственно; ${{\tau }_{x}},\;{{\tau }_{y}}$ – компоненты напряжения трения ветра; ρ – плотность воды; l – параметр Кориолиса, g – ускорение свободного падения, ${{\kappa }_{{v}}},\;{{k}_{h}}$ – коэффициенты турбулентной диффузии температуры и солености по вертикали и горизонтали соответственно, С – источники субстанции X. Уравнение состояния для плотности морской воды имеет вид

Условие отсутствия нормального потока принимается на всех границах. На границах с сушей также полагаются равными нулю нормальные компоненты потоков тепла и солей. Океан подвергается воздействию напряжения трения ветра на поверхности. Потоки T и S у дна принимаются нулевыми, а на поверхности определяются взаимодействием с воздухом.

Уравнения дискретизируются на горизонтальной сетке Аракавы (см. [3], [7]) c применением простых центральных разностей по пространству для диффузионных членов и схемой с весами вверх по потоку для адвективных. Явные конечные разности по времени для соответствующих уравнений дают требуемую точность, и хотя шаг по времени численно ограничен, являются более эффективными, чем центральные разности по времени с большим шагом по времени. Возможно также использование в программе неявного алгоритма (см. [5]). На всех шагах по времени поле скоростей вычисляется диагностически из поля плотностей с использованием уравнений импульса по горизонтали.

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

В настоящей реализации модели используется восемь вертикальных слоев для плотности в растянутой логарифмической шкале. Поверхность океана соответствует уровню с номером восемь. Максимальная глубина океана принимается равной 5000 м.

В термодинамической модели (см. [3], [7]) морского льда уравнения, описывающие его эволюцию, решаются для сплоченности льда и для средней толщины льда. Нарастание и таяние льда в модели определяются разностью между атмосферным потоком тепла в морской лед и тепловым потоком изо льда в океан. Температура поверхности льда в каждый момент времени определяется из уравнения для распределения температуры по толщине льда.

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

где K_hi – эффективный коэффициент горизонтальной диффузии. Скорость роста G_i толщины морского льда в части океана, уже покрытой льдом, вычисляется из разности тепловых потоков в морской лед и обратно, с учетом латентных тепловых потерь из-за сублимации. Образование снега в модели явно не учитывается, все осадки над океаном или морским льдом попадают непосредственно в верхний слой океана. В области расчетной ячейки океана, свободной ото льда, определенный дефицит тепла приводит к росту льда в открытой области ячейки и скорость роста льда в открытой области ячейки задается G₀.

Скорость изменения сплоченности льда A, т.е. доли площади ячейки океана, занятой льдом, равна

Первый член в правой части уравнения определяет потенциальный рост льда на свободных поверхностях океана. Значение этого члена заключается в том, что если G₀ положительно, то доля поверхности без льда убывает экспоненциально со скоростью G₀/H₀, где H₀ – минимально допустимая толщина льда. Второе слагаемое в правой части уравнения определяет потенциальное таяние льда и соответствует скорости, с которой площадь A будет уменьшаться, если весь лед будет однородный по толщине от 0 до 2H/A в части ячейки А, занятой льдом.

Все блоки модели связаны условиями непрерывности потоков импульса, тепла и влаги. Применяются реальная конфигурация материков и структура глубин мирового океана (см. [3], [8]). Уравнения в сферической системе координат решаются численным конечно-разностным методом. Глубина океана представляется в виде восьмиуровневой логарифмической шкалы до 5000 м.

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

При расчете современного состояния климата начальные условия системы задаются постоянными температурами мирового океана, атмосферы и нулевыми скоростями течений океана. Численные эксперименты показывают, что модель достигает равновесия за период около 2000 лет (см. [8]).

3. МОДЕЛЬ ОБЩЕЙ ЦИРКУЛЯЦИИ АТМОСФЕРЫ

Модель ОЦА описывает тропосферу, расположенную ниже предполагаемого уровня изобарической тропопаузы (см. [9], [10]). Как указано ниже, гидростатическое приближение по вертикали позволяет использовать вертикальные координаты, связанные с давлением атмосферы. В данном случае используется безразмерная σ-система координат (см. [11]):

где p – давление, ${{p}_{T}}$ – постоянное давление на уровне тропопаузы, $pS$ – переменное давление у поверхности Земли. По определению, на тропопаузе (верхняя граница тропосферы) $\sigma = 0$ и у поверхности земли $\sigma = 1$.

Уравнения горизонтального движения (в σ-системе координат) могут быть записаны в векторной форме следующим образом:

где для вектора $A = ({{A}_{\lambda }},{{A}_{\varphi }})$, λ – долгота и φ – широта точки. Здесь V – вектор горизонтальных компонент скорости, $\pi = pS - {{p}_{T}}$, $\dot {\sigma } = d\sigma {\text{/}}dt$, f – параметр Кориолиса, k – единичный вектор в вертикальном направлении, α – удельный объем, F – горизонтальная сила трения, Φ – геопотенциал.

Термодинамическое уравнение энергии имеет вид

где c_p – удельная теплоемкость сухого воздуха, T – температура воздуха, $\dot {H}$ – скорость выделения тепла в единице массы воздуха.

Уравнения неразрывности и переноса влаги, соответственно,

где q – отношение смеси водяного пара, равное отношению массы водяного пара (в граммах) к массе сухого воздуха (в килограммах) в единице объема, $\dot {Q}$ – скорость генерации влаги в единице массы воздуха.

Приведенные выше уравнения являются прогностическими (т.е. эволюционными). К ним добавляются уравнение состояния воздуха $\alpha = RT{\text{/}}p$, где R – газовая постоянная для влажного воздуха, и диагностическое уравнение гидростатики, в которое вырождается вертикальная проекция уравнения движения:

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

Для численного решения задачи атмосфера разбивается на слои в вертикальном направлении пропорционально массе воздуха (давлению) (см. [9], [12]). Количество слоев варьировалось от двух до 18 в разных экспериментах.

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

Для определения источников водяного пара $\dot {Q}$ и тепла $\dot {H}$ применяются точечные модели, описывающие гидрологический цикл и процессы распространения теплового и солнечного излучения (см. [13]). Явная численная схема, используемая в модели, накладывает ограничение сверху на шаг интегрирования по времени, который при данном пространственном разрешении модели не превышает 1 ч (см. [14]).

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

Испарение, конденсация и конвективные процессы зависят от термического состояния атмосферы, которое, в свою очередь, является функцией обмена теплом, имеющим место в этих процессах. Вместо получения одновременного решения для влаги и термического состояния, в модели рассчитываются испарение и компоненты конденсации последовательно. На каждом шаге определяется термическое состояние атмосферы, новые значения температуры используются на следующем шаге. При этом для описания конвективных процессов применяется так называемая процедура конвективного приспособления (см. [10]). Она состоит в следующем. Во-первых, температурный градиент между слоями приводится к сухоадиабатическому градиенту, если обнаруживается сухоадиабатическая неустойчивость. Во-вторых, если воздух в верхних слоях перенасыщен, здесь происходят крупномасштабная конденсация и приведение температуры и отношения смеси к устойчивому состоянию. В-третьих, градиенты температуры между уровнями и влажность проверяются на существование влажноконвективной неустойчивости. Если обнаруживается неустойчивость, происходит конденсация, а температуры и отношения смеси адаптируются в результате конвективных процессов. После конвективного приспособления в модели предусматривается переход к расчету влаги в слоях и конденсации.

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

4. ОРГАНИЗАЦИЯ И РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

В модели океана используется по горизонтали расчетная сетка 72 × 72 ячеек, неравномерная по широте, в то время как модель ОЦА рассчитывается на сетке 72 × 46 ячеек, равномерной по широте (46 ячеек) и долготе (72 ячейки).

Расчетная сетка модели океана – равномерная по долготе с тем же шагом 5°, что и для модели атмосферы. Сетка модели океана имеет переменный шаг по широте (сгущается от полюсов к экватору). Для определения значений климатических характеристик в этих точках развита и применяется процедура кусочно-линейной интерполяции функции двух переменных (см. [15], [16]).

В случае интерполяции массива 72 × 72 из модели океана на сетку 72 × 46 атмосферы необходима дополнительная корректировка массива 72 × 46, связанная с тем, что в модели атмосферы используется более точная карта суши, определяющая тип поверхности в данной точке (океан, суша, лед и т.п.).

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

Некоторые результаты совместных расчетов по модели общей циркуляции атмосферы и модели термохалинной циркуляции океана показаны на фиг. 2–8. На фиг. 2 приведены графики средней глобальной и приземной температуры атмосферы в зависимости от времени. Эти кривые демонстрируют выход на установившийся режим и наличие межгодичной изменчивости температуры атмосферы в установившемся состоянии климатической системы.

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

Поля температуры атмосферы вблизи подстилающей поверхности и давления на уровне моря (мбар) для января из модели ОЦА при совместном расчете с моделью термохалинной циркуляции океана обладают сильной изменчивостью над материками, что видно на фиг. 4, 5.

На фиг. 6 и 7 представленo распределение температур океана и меридиональной функции тока, соответственно, в вертикальном сечении до глубины 5000 м в среднем для Атлантического океана. По глубине использована логарифмическая расчетная шкала, сгущающаяся к поверхности. Распределение температуры демонстрирует большие вертикальные градиенты в термоклине вблизи поверхности океана в слое толщиной около 1000 м и более однородное поле в глубоких слоях. На фиг. 7 изолинии определяют направление и величину меридионального потока в Свердрупах (1 Св = 10⁶ м³/с). В верхних слоях океана он направлен от экватора к северному полюсу, где холодные тяжелые массы воды опускаются вниз и более медленно смещаются в обратном направлении, постепенно замещая верхние слои в экваториальной области.

По комплексу моделей на первом этапе проведены расчеты прогнозирования климата до 2100 г. с использованием сценариев роста СО₂ (фиг. 8) под названием РТК8.5 (концентрация 860 ppm в 2100 г.) и РТК4.5 (концентрация 560 ppm в 2100 г.), предложенных Межправительственной группой экспертов по изменению климата (МГЭИК) (см. [17]) и отличающихся прогнозами развития мировой энергетики.

Для первого из них среднеповерхностная глобальная температура атмосферы к 2100 г. выросла на 2.7°С, влажность атмосферы – на 11.5%, уменьшение толщины морского льда – на 25% (фиг. 9). Арктический морской лед летом практически полностью растаял. Площадь вечной мерзлоты на территории России сократилась на 22%. Повышение приземной температуры атмосферы больше над материками, в средних и высоких широтах, достигая величины 5.2°С в северных областях Евразии (фиг. 10). В южном полушарии потепление не превышает 2°С. При реализации второго сценария среднеповерхностная глобальная температура атмосферы к 2100 г. выросла на 1.4°С, влажность атмосферы – на 8.0%, уменьшение толщины морского льда составило 15%.

Изменение температуры верхнего слоя океана для сценария РТК8.5 ожидаемо составляет меньшие значения (фиг. 11), чем для атмосферы, не превышая 2°С. В низких и средних широтах океан прогрет достаточно равномерно на величину около 1.8°С.

К концу века согласно проведенным расчетам прогнозируется существенное уменьшение мощности меридионального потока воды в Атлантическом океане при реализации сценария РТК8.5. На фиг. 12 показан прогноз изменения среднего меридионального потока в Атлантическом океане для 2100 г. при реализации сценария РТК8.5 роста CO₂ по сравнению с его значением при современном климате (фиг. 7). Наблюдается значительное уменьшение потока максимально на 27%, что означает уменьшение потока теплых масс воды из зоны экватора в северные области Атлантики.

5. ЗАКЛЮЧЕНИЕ

В настоящей работе описана организация взаимодействия модели ОЦА и модели термохалинной циркуляции океана и рассматривается функционирование полученной глобальной климатической модели в режиме сезонного хода солнечной радиации. Предложена схема их взаимодействия через граничные условия на поверхности океана. Реализована процедура интерполяции данных на расчетных сетках моделей с учетом конфигурации материков и океанов. Проведены долговременные расчеты на период более 400 лет по совместной модели, которые показали ее устойчивую работу. По комплексу моделей проведены расчеты прогнозирования климата до 2100 г. с использованием различных сценариев роста СО₂, предложенных Межправительственной группой экспертов по изменению климата и отличающихся прогнозами развития мировой энергетики.