CC BY
64
МОДЕЛИРОВАНИЕ КЛИМАТА ¡
УДК 551.509.313.41
В. П. П а р х о м е н к о
МОДЕЛЬ КЛИМАТА С УЧЕТОМ ГЛУБИННОЙ ЦИРКУЛЯЦИИ МИРОВОГО ОКЕАНА
Рассмотрена модель климата, включающая блоки океана, атмосферы и морского льда, взаимодействующие между собой. Модель описывает глубинную термохалинную циркуляцию Мирового океана и основные характеристики остальных элементов климатической системы. В работе рассматривается функционирование модели в режиме сезонного хода солнечной радиации. Рассчитаны изменения температуры атмосферы в XXI веке для различных сценариев изменения концентрации С02.
E-mail: parhom@ccas.ru
Ключевые слова: климатическая модель, термохалинная циркуляция.
Глобальная модель климата включает блоки океана, атмосферы, морского льда. Описывает термохалинную циркуляцию Мирового океана и основные характеристики остальных элементов климатической системы. В настоящей работе исследуется функционирование модели в режиме сезонного хода солнечной радиации.
Система уравнений модели океана рассматривается в геострофическом приближении [1] с фрикционным членом в уравнениях импульса по горизонтали. Значения температуры T и солености S удовлетворяют адвекционно-диффузионным уравнениям, что позволяет описать термохалинную циркуляцию океана. Учитывается также процедура конвективного приспособления.
Таким образом, система основных уравнений в безразмерной форме, записанных в сферических координатах (ф, s, z), где ф — долгота, s = sin в, в — широта и z — высота, направленная вверх, имеет следующий вид [2]:
уравнения импульса по горизонтали
-sv = -1 ^ - Xu + дт\ (1)
c дф dz
su = -cdp - Xv + — ts ; (2)
ds dz
уравнение гидростатики
уравнение неразрывности
d (u^ д , N dw
дф
Vе/
+ ^(vc ) + ~Г = 0; (4)
ds dz
уравнение состояния морской воды
Р = р(?, Т); (5)
уравнение переноса и диффузии трассеров X (температуры и солености)
dX = к V2 X + А
dt dz
' дХЛ Kv Т"
V dz
+ C. (6)
Здесь u, v, w — компоненты вектора скорости; Л — переменный в пространстве фрикционный член, увеличивающийся к береговым границам и экватору; c = cos 0, T, S,p — температура, соленость, давление
соответственно; т = (Т ,ts ) — безразмерное напряжение трения ветра; р и р0 — плотность воды и ее характерное значение; к Kh — коэффициенты турбулентной диффузии трассеров по вертикали и горизонтали соответственно.
Характерный масштаб для расстояния по горизонтали определяется радиусом Земли r0, а по вертикали — максимальной глубиной океана D. Горизонтальные компоненты скорости (u, v) в направлениях (ф, s) выражаются через характерную скорость U0, а вертикальная компонента w — через U0D/r0. Характерные масштабы для давления p и плотности р получаются из геострофических и гидростатического соотношений, соответственно. Следовательно, градиент плотности
G = др/dz выражается через р0fU0r0/(gD2), гдеf— удвоенная угловая скорость суточного вращения Земли. Характерное время определяется выражением r0/U0. Величина d/dt в уравнении (6) — материальная (полная) производная. Масштабные множители для T и S не используются. Величина C определяется из процедуры конвективного приспособления для устранения статической неустойчивости с учетом консервативности T и S. Уравнение состояния для размерной плотности р* имеет вид
р* = 1000 + 0,7968S - 0,0559T - 0,0063T2 + 3,731510-5 T3. (7)
Условие отсутствия нормального потока требуется на всех границах. На границах материков также принимаются равными нулю нормальные составляющие потоков теплоты и солей. Океан подвергается воздействию напряжения трения ветра т на поверхности. Потоки T и S у дна полагаются равными нулю, а на поверхности определяются
взаимодействием с атмосферой. Поверхностный температурный поток ¥т связан с поверхностным тепловым потоком Qв соотношением Qв = р0ер0Гт, где ср0 — удельная теплоемкость морской воды.
Уравнения дискретизируются на сетке Аракавы [3] с использованием простых центральных разностей по пространству для диффузии и схемой с весами вверх по потоку для адвекции. Простые явные конечные разности по времени обеспечивают требуемую точность и, хотя шаг по времени численно ограничен, являются более эффективными, чем центральные разности по времени с большим шагом по времени. Неявный алгоритм [4] также может быть использован в программе, но для стандартных параметров он является менее эффективным. На каждом шаге по времени поле скоростей определяется диагностически из поля плотностей.
Вертикальные уровни модели равномерно распределены в логарифмических координатах £ = 1о§(1 - 2 + 0,1) так, что верхние слои тоньше, чем нижние. Горизонтальная сетка является равномерной в (<, э)-координатах (долгота и синус широты), определяя при этом ячейки одинаковой площади в пространстве. В настоящей модели используется восемь вертикальных уровней для плотности. Максимальная глубина принимается равной 5 км.
Для описания процессов, протекающих в атмосфере, используется энерго- и влагобалансовая модель или модель общей циркуляции атмосферы. Для первой из них прогностическими переменными являются температура воздуха Та и удельная влажность qa на подстилающей поверхности. В модели решается вертикально проинтегрированное уравнение для Та, определяющее баланс приходящего и уходящего радиационных потоков, явных (турбулентных) обменов потоками теплоты с подстилающей поверхностью, высвобождения скрытой теплоты из-за осадков и простой однослойной параметризации горизонтальных процессов переноса. Источники в уравнении переноса для удельной влажности q определяются осадками, испарением и сублимацией с подстилающей поверхности. Уравнения для баланса теплоты и влаги атмосферы (на единицу площади) имеют следующий вид:
ph
' a q
Pahtc
гдТ Л
-f + ß-W(uTa )-V(vVTa) = Qta; (8)
öt j
+ ßqV(uqa )-V(*?qa )) = Po (E - P), (9)
dt
где к и к — толщины атмосферных пограничных слоев для темпера-
I q
туры (8,4 км) и влажности (1,8 км) соответственно; V и к— коэффициенты турбулентной диффузии для температуры и влажности соответ-
ственно; Qta — суммарный поток теплоты в атмосферу; Е — скорость испарения или сублимации; Р — скорость выпадения осадков; ра и р0 — плотность воздуха и воды; сра — удельная теплоемкость воздуха
при постоянном давлении. Параметры вТ и в являются масштабны-
Т Я
ми множителями для оценки вклада адвективного переноса. Они могут быть необходимы вследствие однослойного представления атмосферы, особенно если при счете используются данные о поверхностных скоростях взамен вертикально усредненных данных. В работе [5] вТ = 0 и вд = 0,4 или 0.
Полный поток теплоты в атмосферу
Qaa = QsWCA + Q ьш - + Q 5Н + Q ьН
(10)
Здесь QSW — приходящая коротковолновая солнечная радиация, задающаяся соотношением
Qsw = 5с1 (в)(1 - а), (И)
где 5с — солнечная постоянная; I — множитель, учитывающий широтное распределение солнечной радиации; а — планетарное альбедо, которое над океаном и сушей определяется косинусом широты в соответствии с данными наблюдений; СА — коэффициент поглощения, параметризующий поглощение солнечной радиации водяным паром, пылью, озоном, облаками и т. д. Над океаном СА = 0,3 (т. е. океан поглощает 70 % инсоляции), над сушей СА = 1,0.
Суммарный поверхностный поток длинноволновой радиации QLW в атмосферу равен разности между длинноволновым излучением подстилающей поверхности и атмосферы:
QLW = срт: - еаоТ:, (12)
где е и £а — коэффициенты излучательной способности поверхности и атмосферы; Т — температура подстилающей поверхности (над сушей Т = Та); о— постоянная Стефана - Больцмана. Уходящая планетарная длинноволновая радиация QPLW должна учитывать парниковый эффект из-за наличия атмосферного водяного пара. В связи с отсутствием в модели схем для описания облачности и радиации используется полиномиальная функция [5, 6], кубическая по температуре Т и квадратичная по относительной влажности г = я /я (я — удельная влажность насыщения). Также учитывается член, отвечающий за нагрев атмосферы в результате увеличения концентрации С углекислого газа по отношению к значению исходной концентрации С0:
QPLW = Е Е с/ТА С, (13)
1=0 ]=0 С0
где ЛР2 = 4/1п2 — показатель радиационного воздействия при удвоении концентрации углекислого газа.
Поток явной теплоты QSH параметризует вертикальньный турбулентный обмен с подстилающей поверхностью и задается соотношением
Qsн = РаСнСраи ( - Та ), (14)
где и — скорость ветра у поверхности (получена из напряжения трения ветра).
Число Стэнтона, характеризующее интенсивность диссипации энергии в потоке жидкости или газа,
Сн = 0,9Се, (15)
где СЕ — число Дальтона:
С =
Е (16)
= тах ^ 6 • 10-5, шт^ 2,19 • 10-3, 10-3 (1,0022 - 0,0822 (Та -Т5) + 0,026би))). Латентное выделение теплоты при выпадении осадков
QLH = Р0 К P, (17)
где Lv — латентная теплота парообразования.
Значение Р вычисляют в предположении, что вся избыточная влага мгновенно удаляется прямо за один шаг по времени, когда относительная влажность г превышает пороговое значение гшах. Скорость испарения или сублимации
Е = ВШ-( - яа), (18)
р0
где qs — удельная влажность насыщения:
dT
ls
С
qo exp
\ s s
s = a,0, i. (19)
При значении индекса ^ = а используется температура атмосферы для расчета скорости выпадения осадков Р. Индекс ^ = 0 относится к температуре океана, необходимой для расчета испарения с поверхности океана, а ^ = 1 определяет температуру поверхности льда, требуемую для подсчета сублимации морского льда.
Температура на поверхности принимается равной температуре атмосферы Та, и испарение приравнивается к нулю; таким образом, источники атмосферной теплоты упрощаются, так как QLW = Qsн = QLH = = 0. Осадки над поверхностью земли прибавляются к стоку в океан в соответствующих ячейках.
В модели эволюции морского льда часть океанической поверхности, покрытой морским льдом в любой заданной ячейке (сплоченность льда),
обозначается A. Динамические уравнения решаются для A и для средней толщины льда H. Отметим, что H представляет собой усредненную толщину льда в ячейке с учетом как открытой поверхности океана, так и покрытой льдом. Среднее значение толщины льда на покрытой льдом части площади ячейки определяется выражением H/A. Для температуры поверхности льда T. решается диагностическое уравнение.
Рост и таяние льда в модели зависят только от разности между потоком теплоты из атмосферы в морской лед и потока теплоты изо льда в океан. Так как A — это часть ячейки, покрытая льдом, то поток любой субстанции F между атмосферой и океаном (или морским льдом) определяется взвешенной суммой составляющих потоков Fi и F0 над открытой и покрытой льдом частями:
F = AFx + (1 - A)Fo. (20)
Потоки теплоты из атмосферы в океан и область океана, покрытую льдом, соответственно равны:
0,0 = (1 - Ca )Qsw - Qlwo - Qsho - PooLE (21)
Qi=(l - Ca )QsW - QlW - QsHi - PoLE. (22)
Члены в правой части этих уравнений определены соотношениями (11) - (19) с параметрами для океана в (21) и для льда — в (22). Потери теплоты из-за испарения или сублимации заменяют приток теплоты от выпадения осадков. В уравнении (22) Ls является скрытой теплотой сублимации.
Скорость роста Gi морского льда в части океана, покрытой льдом, определяется из разности тепловых потоков в морской лед и обратно, минус латентные тепловые потери из-за сублимации. Образование снега не рассматривается в модели, все осадки над океаном или морским льдом добавляются непосредственно в поверхностный слой океана.
Отсюда можем подсчитать суммарную скорость роста G и скорость изменения средней толщины льда, которая подвергается также влиянию адвекции поверхностными течениями океана и диффузии, имеющими место при детальном представлении подсеточных процессов реологии и адвекции морского льда: dH
— + H = AGi + (1 - A) Go = G, (23)
где Khi — эффективный коэффициент горизонтальной диффузии.
Скорость изменения доли A площади ячейки океана, покрытой льдом, равна
dA ( G Л ( A Л
— + KhVlA = max 0, (1 - A)G + min 0, AG — . (24) dt H r.
2 H )
Первый член в правой части уравнения (24) описывает возможный рост льда на открытых поверхностях океана. Влияние этого члена заключается в том, что если G0 положительно, то доля поверхности без льда убывает экспоненциально со скоростью 00/Н0, где Н0 — минимально допустимая толщина льда. Второй член этого уравнения описывает возможное таяние льда и отвечает за скорость, с которой площадь А будет уменьшаться, если весь лед будет равномерно распределен по толщине от 0 до 2Н/А в части ячейки А, покрытой льдом.
Окончательно на каждом шаге необходимо убедиться, что вычисленная толщина льда положительна. Практически игнорируется присутствие тонкого льда и считается Н = 0 всюду, где Н < Н0. Это означает, что потоки теплоты и пресной воды должны быть модифицированы согласно новому состоянию ледового покрова. Для положительных или отрицательных значений Н, если Н становится равным нулю, соответствующее количество тепла (-НрХу) добавляется в океан; при этом соответствующее добавленное количество пресной воды дается выражением -Нр/р0. Кроме того, необходимо контролировать, что ошибки аппроксимации не приводят к нефизичным значениям А (т. е. А < 0 или А > 1).
Все блоки климатической модели связаны между собой обменом импульсом, теплотой и водой. Обмен импульсом состоит только в использовании скорости верхнего слоя океана для адвекции морского льда. Всеми другими обменами импульсом пренебрегают.
Потоки теплоты между «смежными» блоками могут быть модифицированы фазовыми переходами на границах (испарением, таянием и т. д.). Потоки из одного блока в другой могут отличаться на величину, определяемую латентными тепловыми эффектами. При этом материковые стоки воды Я добавляются в океанические ячейки на каждом шаге по времени.
Полный обмен потоками теплоты и воды над океаном определяется соотношением (20) и зависит от сплоченности А морского льда. Модель морского льда является связующим звеном всех трех компонент как в теории, так и в программной реализации.
Поток теплоты в атмосферу
Qta = AQ1a +(1 - ,
где Qia и Qoa определяются из уравнения (10). В этом соотношении член Qsw отвечает за распределение падающей солнечной радиации между атмосферой и подстилающей поверхностью, т. е. связывает атмосферу и морской лед посредством зависимости альбедо морского льда от температуры воздуха. Отметим, что выпадение осадков, испарение и сублимация являются причиной нетривиальной сшивки блоков модели из-за пресноводных потоков, а не из-за тепловых, потому
что соответствующие тепловые источники и стоки в каждой из компонент зависят только от ее внутренних переменный. Из условий сохранения энергии члены QLW и QSH из (22) и (21) должны иметь противоположные знаки в выражениях для потоков теплоты из атмосферы в морской лед и океан. Поток теплоты изо льда в океан задается соотношением (27), но суммарный поток распределяется между морским льдом и открытыми частями океана.
Поток пресной воды в атмосферу (9) равен Е - Р с учетом испарения с поверхности Земли и сублимации морского льда. Считаем, что осадки выпадают непосредственно в океан, без учета присутствия льда, а испаренная или сублимированная вода удаляется из океана или льда соответственно. В формулировке модели океана в приближении «твердой крышки», используемой здесь, модель представляет океан как неисчерпаемый источник пресной воды для морского льда и атмосферы. Пресная вода в них сохраняется, но конвертируется в соленость на поверхности океана. В модели океана соленость является консервативной переменной вследствие использования постоянного множителя перехода S0 от количества пресной воды к солености.
Для расчета одного шага по времени для океана, морского льда и поверхностных потоков требуется несколько шагов по времени для атмосферы. Для атмосферы шаг примерно равен одним, а для океана — нескольким суткам. Как отмечено выше, расчеты для морского льда и поверхностных потоков тесно связаны между собой и играют роль связующего звена. Все потоки между компонентами вычисляются в одни и те же моменты времени для гарантии консервативности, но их значения берутся с предыдущего шага по времени, чтобы избежать сложностей неявной схемы. Диагностическая поверхностная температура морского льда удовлетворяет неявному уравнению. Дальнейшее усложнение состоит в том, что определенные ограничения на толщину льда и его площадь могут быть оценены после обновления характеристик морского льда, так как зависят от потока пресной воды в океан. По этой причине обновление характеристик морского льда происходит после определения других потоков.
Прогностические уравнения модели (6), (8), (9), (23) и (24) решаются методом центральных разностей второго порядка по пространству и простык разностей вперед по времени. Как альтернативный, можно применить метод решения уравнений (6), (8) и (9) с помощью неявной схемы по времени, что в принципе создает возможность использования больших шагов по времени. Неявная схема включает в себя неявный шаг — предиктор с переменным числом итераций [3], за которым следует шаг — корректор. Обычно используется четыре итерации. Большее количество итераций повышает точность, но может повлиять на устойчивость системы. Неявная схема является сильно ограничен-
ной для океана и не имеет преимуществ перед более простым явным методом. Для атмосферных уравнений она приводит к значительным улучшениям в эффективности модели.
В модели используется равномерная по долготе и синусу широты конечно-разностная сетка 36x36 ячеек (рис. 1). Разрешение модели по долготе составляет 10°, а по широте оно изменяется приблизительно от 3° у экватора до 20° у полюсов. Глубина океана представляется в виде восьмиуровневой логарифмической шкалы до 5000 м (рис. 1).
Рис. 1. Расчетная сетка и восьмиуровневая шкала глубин (м). Серым цветом выделена суша
Численные эксперименты показывают, что модель выходит на равновесие за период около 2000 лет. Начальное состояние системы характеризуется постоянными температурами океана, атмосферы и нулевыми скоростями течений океана.
В целях исследования чувствительности термохалинной циркуляции и Северо-Атлантического океанского переноса теплоты к коэффициентам диффузии к и кь и неявному переносу пресных вод из Атлантического океана в Тихий океан проведены серии экспериментов. Наиболее близкие к данным наблюдений результаты получаются при следующих значениях параметров: к = 1,0.10-4 м-2с-1, К = 2000 м-2. с-1. У
На рисунках представлены некоторые результаты расчетов. Они показывают достаточно хорошее воспроизведение основных характеристик глобальной климатической системы.
На рисунках далее показан выход средней глобальной температуры атмосферы (рис. 2, а) и удельной влажности (рис. 2, б) на современные
1
1
/
1
0,012 0,011 0,01
л
§ 0,009 г
| 0,008 ш
0,007 0,006 0,005
400 800 1200 1600 2000 Годы
v:
400 800 1200 1600 2000 Годы
Рис. 2. Выход средних глобальных температуры атмосферы (а) и удельной влажности (б) на стационарные значения:
1 — среднегодовой режим; 2 — режим с сезонным ходом инсоляции
значения этих показателей, т. е. температуру около 14,3 °С и удельную влажность около 11,3 г/кг. Как видно, эти глобальные характеристики при расчете в среднегодовом режиме и с учетом сезонного хода практически не различаются. То же относится и к площади морского льда (рис. 3), которая стабилизируется приблизительно через 1000 расчет-
Рис. 3. Выход площади морского льда на стационарные значения: 1 — среднегодовой режим; 2 — режим с сезонным ходом инсоляции
ных лет. Однако средняя толщина морского льда сильно различается для этих двух вариантов расчета.
На рис. 4 и 5 показано распределение некоторых основных характеристик климата для зимних и летних условий.
Межправительственная группа экспертов по изменению климата (МГЭИК) разработала долгосрочные сценарии эмиссии парниковых газов и аэрозоля в атмосферу в XXI в., которые опубликованы в Специальном докладе о сценариях выбросов (СДСВ) [10]. По сценарию А2, развитие мира проходит при постоянном росте общей численности
■260 -220 -180 -140 -100 -60 -20 20 60 100
Рис. 4. Поток теплоты в атмосферу для января (Вт/м2). Отрицательные значения выделены серым цветом
Рис. 5. Температура верхнего слоя (175 м) океана для июля (°С)
населения в мире. Экономическое развитие имеет региональную направленность, а экономический рост в расчете на душу населения и технологические изменения происходят медленнее по сравнению с другими сценариями. В результате такого развития ожидается значительное увеличение концентрации основных парниковых газов в атмосфере. К 2100 г. концентрации основных парниковых газов в атмосфере увеличатся по сравнению с показателями 1990 г.: С02 — в 2,42, СН4 — в 2,19 и К20 — в 1,45 раза. Сценарий В1 содержит описание мира с глобальным населением, которое достигнет максимальной числен-
Рис. 6. Временной ход средней годовой аномалии температуры атмосферы (°С) в XXI в. по отношению к базовому периоду (1980-1999 гг.) с учетом неизменных условий (1), сценариев В1 (2) и А2 (3)
Рис. 7. Изменение температуры атмосферы в конце XXI в. при увеличении концентрации С02 согласно сценарию А2. Лето
ности к середине XXI в., а затем сократится при быстрых изменениях в экономических структурах с уменьшением материальной интенсивности и внедрением экологически чистых и ресурсосберегающих технологий. К 2100 г. концентрация С02 и К20 в атмосфере увеличится по сравнению с данными 1990 г. соответственно в 1,53 и 1,22 раза, а концентрация СН4 уменьшится на 6 %.
Рассчитанное по описанной модели изменение средней глобальной температуры атмосферы в конце XXI в. для сценария А2 составляет +2,2 °С, для сценария В1 — +1,4 °С, при неизменной концентрации — +0,35 °С (рис. 6). Глобальное распределение изменения температуры для июля при увеличении концентрации СО2, согласно сценарию А2, показано на рис. 7.
Работа выполнена в рамках Программы фундаментальных исследований президиума РАН № 14, проектов РФФИ 11-01-93003 и 11-01-00575.
СПИСОК ЛИТЕРАТУРЫ
1. К о ч е р г и н В. П. Теория и методы расчета океанических течений. М.: Наука, 1978. 128 с.
2. M a r s h R., E d w a r d s N. R., S h e p h e r d J. G. Development of a Fast Climate Model (C-GOLDSTEIN) for Earth System Science // SOC. 2002. No. 83. 54 p.
3. A r a k a w a A., L a m b V Computational Design of the Basic Dynamical Processes of the Ucla General Circulation Model. In Methods in Computational Physics. Academic Press, 1977. Vol. 17. P. 174-207.
4. S h e p h e r d J. G. Overcoming the CFL Time-Step Limitation: a Stable Iterative Implicit Numerical Scheme for Slowly Evolving Advection-Diffusion Systems // Ocean Modelling. 2002. Vol. 4. P. 17-28.
5. W e a v e r A. J., E b y M., W i e b e E. C., et al. The UVic Earth System Climate Model: Model Description, Climatology, and Applications to Past, Present and Future Climates // Atmos-Ocean. 2001. Vol. 39. P. 361-428.
6. T h o m p s o n S. L., and W a r r e n S. G. Parametrization of Outgoing Infared Radiation Derived from Detailed Radiative Calculations // J. Atmos. Sci. 1982. Vol. 39. P. 2667-2680.
7. H o l l a n d D. M., M y s a k L. A., M a n a k D. K., O b e r h u b e r J. M. Sensitivity Study of a Dynamic Thermodynamic Sea Ice Model // J. Geophys. Res. 1993. Vol. 98. P. 2561-2586.
8. M i l l e r o F. J. Annex 6, Freezing Point of Seawater // Unesco Tech. Papers in the Marine Sciences. 1978. Vol. 28. P. 29-35.
9. M c P h e e M.G. Turbulent Heat Flux in the Upper Ocean Sea Ice // J. Geophys. Res. 1992. Vol. 97. P. 5365-5379.
10. Nakic enovic N., et al., IPCC Special Report on Emission Scenarios. United Kingdom and New York, NY, USA, Cambridge University Press, 2000.
Статья поступила в редакцию 27.10.2011.
Пархоменко Валерий Павлович, в 1974 г. окончил физический факультет МГУ им. М.В. Ломоносова, канд. физ.-мат. наук. В настоящее время зав. сектором моделирования климатических и биосферных процессов Вычислительного центра РАН (ВЦ РАН), доцент МГТУ им. Н.Э. Баумана. Имеет более 90 печатных работ. Научные исследования посвящены нестационарной газовой динамике, численным методам динамики сплошных сред, проблемам математического моделирования климата. За последний период им выполнены разработка и адаптация глобальных климатических численных моделей, используемых в исследованиях в ВЦ РАН, для проведения и обработки численных расчетов при моделировании и прогнозировании глобальных климатических и биосферных процессов на базе параллельных технологий. Разработаны варианты моделей циркуляции атмосферы, гидротермодинамики океана и морского льда в Арктическом регионе. Проведено множество компьютерных экспериментов с моделями.
