CC BY
18
11. Allen, M.P. Computer simulation of liquids |Text| / M.P. Allen, D.J. Tildesley.- Oxford: Clarendon Press, 1987,- 387 p.
12. Tersoff, J. Modeling solid-state chemistry: Interatomic potentials for multicomponent systems
I Text I / J. Tersoff// Phys. Rev. B. -1989. -Vol. 39,-P. 5566-5568.; 1990,- Vol. 41,- P. 3248.
13. Ziegler, J.F. The stopping and range of ions in solids [Text] / J.F. Ziegler, J.P. Biersack, U. Littmark / / The Stopping and Range of ions in Matter.— 1995.— Vol. 1,— New York : Pergamon, 1985,— 321 p.
УДК 519.71 1.3
И.А. Судаков
ДИНАМИКА ПРОТАИВАНИЯ МЕРЗЛОТНЫХ ОЗЕР И ИЗМЕНЕНИЯ КЛИМАТА
Деградация вечной мерзлоты и дестабилизация мегапула газовых гидратов, захороненных в ней, представляет собой единственный возможный механизм быстрого вовлечения огромных количеств метана в современный цикл углерода, что может привести к климатической катастрофе.
Установлено, что выброс метана в атмосферу из вечной мерзлоты, как минимум, соизмерим с эмиссией метана из морей всего Мирового океана. Летом 2010 года в отдельных районах Арктического региона было обнаружено резкое увеличение атмосферного содержания метана, по сравнению с предыдущими наблюдениями (2006-2009 гг.) до 3,5 мд (3,5-10"5 %), что в два раза выше среднепланетарного значения. Первоочередная задача — это выявление масштабов деградации вечной мерзлоты и количественная оценка потоков метана в атмосферу
Для решения проблемы прежде всего необходимо разработать эффективные модели процессов протаивания вечной мерзлоты на основании которых становится возможной оценка эмиссии метана в атмосферу и прогноз поведения климатической системы в будущем. Очевидно, что при этом возникает сложная задача: связать микроскопические процессы, приводящие к таянию вечной мерзлоты (например фазовые переходы), с макроскопическими, происходящими в климатической системе (например явление положительной обратной связи, приводящее к «парниковому эффекту»).
В зависимости от вида мерзлотных образований существуют различные географические типы источников мерзлотного метана: мерзлые
грунты, торфяники, мерзлотные озера, подводная мерзлота на шельфе.
К мерзлотным относятся озера, имеющие в своем основании многолетнемерзлые грунты и расположенные преимущественно в зоне тундры; они слабо изучены в плане моделирования эмиссии мерзлотного метана [1]. В настоящее время существует ряд моделей, описывающих процессы таяния мерзлоты в таких озерах на основе уравнений переноса тепла. Однако эти модели реализуются только с использованием трудоемких и времязатратных численных методов и не содержат связи микроскопических процессов протаивания мерзлоты с макроскопическими процессами, происходящими в климатической системе планеты.
Основная цель данной статьи — представить модель динамики радиуса мерзлотного озера, разработанную автором на основе теории фазовых переходов Гинзбурга — Ландау с применением асимптотических методов. Кроме того, в статье показана возможность применения данной модели к задаче прогнозирования эмиссии мерзлотного метана в атмосферу
Модель динамики радиуса мерзлотного озера
При моделировании многолетнемерзлых фунтов (например, в мерзлотных озерах) часто предполагают, что они состоят всего из двух фаз — твердой и жидкой, на границе которых происходит фазовый переход «промерзание — оттаивание» [1]. Термический режим многолетнемерзлых грунтов в таких моделях описывается моделями на основе задач типа задачи Стефана.
Рассмотрим цилиндрическую область В = = Пх[0,й], где О — подобласть двухмерной сферы соответствующей поверхности вечной мерзлоты; (х,у) — координаты в подобласти 0;ге[0,Щ — вертикальная координата, соответствующая глубине озера А. Неизвестная функция температуры определяется как и = м(х, у, г, г). Предположим, что /г«Д где Х = ^ат(Х) и граница подобласти О — гладкая кривая. Определим положение фронта протаивания как /(х, >\ /).
Таким образом, уравнение теплопроводности соответственно для жидкой и твердой (ледяной) фаз, имеет вид:
и,=К{Ащ 0<г</(г);
щ = К2А щ т<кк,
граничное условие на глубине
йг(х, у, К О + уЩх,у, К О = У(х, у, (); (2) на поверхности —
и(х,у,0,1) = и(х,у,1); (3)
условие на границе раздела фаз —
п( У(и(х, .у, г, 0 - й(х,у,г,^) |(Л. ^ г)еГ) = (4)
Здесь п — вектор нормали к фронту протаивания; V — скорость фронта; 0 — безразмерная скрытая теплота; К(х,у,/), ¿/(х,у,/) — граничные функции температуры; К1,К2 — коэффициенты теплопроводности жидкой и твердой фаз соответственно; Ь — феноменологический параметр.
Поверхность раздела фаз Г описывается следующим уравнением:
и(х,у,1(х,у,(),{) = д, (5)
где 9 — температура фазового перехода.
Жидкую фазу будем определять как
{м(х, у, г, г) > 0}, а твердую — {м(х, у, г, г) <0}.
Обратимся теперь к неклассической схеме моделирования движения границы раздела фаз твердое тело — жидкость, впервые предложенной в работе [2]. Будем рассматривать границу раздела фаз как область конечной толщины. Далее, как было ранее сделано в [2], представим единую для всех точек рассматриваемой физической системы «твердое тело — граница раздела фаз (конечнойтолщины) — жидкость» систему двух дифференциальных уравнений в частных производных для двух неизвестных полей:
ut = KM-^wiy
а^2фг='2Лф + а_1£(ф) + 2( и-0), (7) где ф = ф(х, y,z, t) — параметр порядка; g =
о
= 0,5(ф-ф )— производная симметричного двухъямного потенциала с минимумами при ф = = +1 (подробнее см. в статье [2]), а,а,Ь,К — безразмерные постоянные, которые должны быть физически интерпретированы в ходе экспериментов [2].
Неизвестные поля в указанной системе урав-
=
и (гладкое) поле параметра порядка ф = ф(х, у, z, t), интерпретируемое в [2] как «локальная средняя фаза». Уравнение (6) является строгим в приближении бесконечно тонкой границы раздела фаз, предполагающемся в классической задаче Стефана. Следуя работе [2], будем использовать уравнение (6) и для описания границы раздела фаз
конечной толщины при гладком поле параметра фф
правую часть уравнения (7), то получится уравнение Эйлера — Лагранжа, которое можно вывести в рамках решения вариационной задачи о минимизации функционала свободной энергии системы, если выбрать для указанного функционала форму, аналогичную предложенной в теории фазовых переходов Гинзбурга — Ландау [2]. Переход к неравновесному приближению обуславливает появление в левой части уравнения (7) частной производной по времени от параметра порядка. Упомянем, что теория фазовых переходов Гинзбурга — Ландау была первоначально развита для фазовых переходов нормальный металл — сверхпроводник, хотя впоследствии она применялась и для других физических систем, например фазовых переходов «порядок — беспорядок» в ферромагнетиках.
В работе [2] было показано, что в ряде предельных случаев описываемый подход, основанный на теории Гинзбурга — Ландау, сводится к трем известным формулировкам задачи о движении границы раздела фаз, а именно, к классической задаче Стефана и двум ее модификациям. Остановимся на одной из таких модификаций.
Пусть в уравнении (7) параметр а фиксирован, параметры и 'а= const; тогда при
Кх= К2 = К приходим к формулировке задачи Стефана в следующем виде [2]:
(Ум(х, у, г, 0 - Vйь(х, у, г, ¿))-п|г =
= Ьу{х,у,1,гу, (8)
Л5((м(х, у, ж, /) |г -6)) = -ок - аор,
(9)
где Д5 — разность энтропий жидкой и твердой фаз; а — поверхностное натяжение; к{1,х) — параметр, описывающий полную кривизну в точке на межфазной границе.
Первое условие (уравнение (8)) совпадает с условием непрерывности теплового потока на границе раздела фаз в классической задаче Стефана. Второе (уравнение (9)) описывает эффект переохлаждения на границе раздела фаз, обусловленный двумя вкладами различной природы: эффектом кривизны границы, описанным еще в работах Гиббса и отвечающим за первое слагаемое в правой части условия (8);
добавочным переохлаждением, описываемым эмпирической поправкой, предложенной в [2], второе слагаемое в правой части уравнения (9).
Применим вышеописанную модель к задаче о динамике радиуса озера в тундре. Рассматриваемая модель математически эквивалентна модели вытаивания «тонкого жидкого диска» на поверхности ледяного полупространства (рис. 1).
Рассмотрим случай, когда в условии (9) разность энтропий жидкой и твердой фаз Д5" = 0. При этом оказывается, что нормальная составляющая скорости движения границы раздела фаз V пропорциональна средней кривизне к указанной границы:
•\>(х,у,г,1) = -\к(х,у,г,1), (10)
где \ — положительная феноменологическая константа.
Горизонтальный радиус, о. е.
Рис. 1. Фронт протаивания грунта
в мерзлотном озере: 1 — исходный; 2— через 300 лет
Для дальнейшего анализа напомним, как вычисляется средняя кривизна произвольной поверхности в какой-либо ее точке [3—6]. Для каждой точки поверхности существует два нормальных сечения с взаимно перпендикулярными плоскостями, которым соответствует наибольшая кри-
=
=
нормальные сечения; , Л2 — соответствующие радиусы кривизны). По определению, средняя кривизна дается формулой
2 1 2 2
Я1 +Я
(Н)
и
Чтобы получить модель расширения (или сужения) озера в горизонтальном направлении, будем исходить из следующих соображений. Поскольку глубина озера много меньше его радиуса, озеро можно моделировать квазидвумерным кругом (точнее, тонким жидким диском) радиуса Я, на образующей которого (вертикальный берег озера) в соответствии с формулой (10) будем иметь выражения: кх =1/Я, к2=$,Н = 1/2 Я, а на плоском дне — кх= к2 = Н = 0.
Соответственно, граница раздела фаз движется в горизонтальном направлении (движение вертикального берега), радиус озера увеличивается (или уменьшается). Используем результаты работы [7]; тогда можно получить следующее дифференциальное уравнение для радиуса озера:
Л
(12)
По аналогии с работами [7—9] модифицируем уравнение динамики изменения радиуса озера (12), добавляя в его правую часть феноменологическую положительную константу 5 (это
Д
ропий жидкой и твердой фаз):
Л
(13)
Возвращаемся к граничному условию (9) модифицированной задачи Стефана и предполагаем теперь, что ДБ ф 0 при большой кривиз-
5
величиной. Физический смысл этого параметра состоит в том, что он описывает плоский фронт протаивания при ненулевой разности энтропий в жидкой и твердой фазах и определяется через микроскопические параметры а, 6, которые мо-
гут быть найдены в ходе экспериментальных исследований динамики протаивания озер в тундре. Фактически теперь соотношение (9) связывает микроскопические процессы, происходящие при фазовом переходе в озере, с макроскопическими процессами роста озер.
Для более реалистичного описания процесса роста озера (в действительности в тундре часть озер испытывает разрастание, а часть — схло-пывание) учтем стохастические эффекты при помощи уравнения Фоккера — Планка:
д(ДЛ)р) = ^
(14)
и количество метана, поступающего из озера в атмосферу, в зависимости от его площади.
Скорость эмиссии метана определяется соотношением
V
= Ро |ехР
кви
ёхёуёг, (16)
дЯ дЯА где р = р(Л,г) — функция, определяющая изменение радиуса озера с течением времени (Я > 0); й — параметр; /=-т/Я (т — постоянная).
Модифицированное уравнение Фоккера — Планка в случае (I = 0 сводится к детерминистическому уравнению (13) динамики радиуса озера. Из уравнения (14) находим вид функции
р(Л) = ^Ягде и — константа.
р
блюдаемому распределению озер в тундре, которое можно определить на основании распределения Парэто:
р(Л) = ^пЛ-т+1\ Ле[Лт{п,«>). (15)
Здесь — минимальный радиус озера, к — положительная константа.
Численные расчеты по стандартным алгоритмам для распределения Парэто соответствуют результатам наблюдений. Последние показывают, что величина максимального наблюдаемого радиуса озера в районе Центральной Якутии за период 1976—2000 гг. лежит в интервале 1000 — 2000 м [10], а по результатам численных расчетов по формулам (8), (9), (13), (14) эта величина составляет в среднем 1600 м.
Применение модели динамики мерзлотного озера к задаче об эмиссии метана из вечной мерзлоты
Теперь на основе полученных результатов построим полуэмпирическое соотношение, описывающие эмиссию метана из системы озер. Действительно, из известных данных наблюдений [11] система озер в тундре претерпевает изменение за счет различных процессов, прежде всего из-за протаивания вечной мерзлоты, следовательно, радиус таких озер изменяется, как
где и = м(х, у, ж, г) — абсолютная температура в точке в данный момент времени; р0 — плотность метана в объеме, изменяющемся стечением времени; ка — постоянная Больцмана; ¿70 — энергия активации.
Уравнение (16) основано на классических представлениях о том, что вероятность продуцирования молекулы метана в процессе реакции разложения метана определяется как
/> = ехр
.На.
кви
Используя уравнение (13), описывающее динамику изменения радиуса озера, и предполагая, что средний радиус системы озер определяется как
N
;=1
N
а общая площадь озер — как $ = , получим
/=1
соотношение для потоков метана из озера:
Гтег = С^ = 2СЩ-^ + ЬЯ). (17) аТ
Здесь коэффициент С зависит от температуры и и определяется как больцмановская экспонента; N — количество озер на всей исследуемой площади.
Усредняя температуру и по некоторому периоду времени с учетом того, что рост озер происходит гораздо медленнее, чем изменения в атмосфере, получим простое соотношение
Ута= р(ЯЛ-1)ехр(-^/"«)> <18)
где р — некоторая постоянная, определяемая из экспериментальных данных, 60, В — положительные константы.
Феноменологические уравнения для эмиссии метана, полученные в работах [12, 13] как из результата эксперимента, так и путем моделирования, имеют вид:
я I
? 7
Щ
Л
О О СЗ
S 6
5 11:1-1-1-1-1-
0 100 200 300 400 500
Время, годы
Рис. 2. Прогнозируемое количество метана в атмосфере Земли на 500 лет с учетом (1) и без учета (2) эмиссии из мерзлотных озер
Гтег = ехр(0,492 + 0,126(Г-273)-0,057И^), (19)
где Гте1 — поток метана с поверхности озера в районе станции Бочкаревская, Г — температура на глубине озера 12 см, ]¥— уровень воды в озере.
При разложении подэкспоненциального выражения формулы (18) в ряд Тейлора по отклонениям температуры от некоторого усредненного ее значения, полученное выражение имеет ту же структуру, что и уравнение (19). Это подтверждает соответствие полуэмпирического соотношения (18), полученного в рамках теории о динамике изменения радиуса озера, соотношениям, полученным экспериментально.
Применим развитый здесь подход для описания динамики радиуса озера к поиску решения важнейшей задачи об увеличении количества метана в атмосфере. Известно, что резкое увеличение концентрации метана в атмосфере за короткий период времени может привести к «парниковому эффекту» и усилить положительную обратную связь в климатической системе [14,15].
Учитывая влияние внешних условий на поверхности озера (глобального климата) на температуру воды в озере, предположим, что усредненная температура озера и линейно зависит от общего содержания метана в атмосфере; тогда на
основании численных расчетов по уравнениям (13), (18) получим, что при типичном значении коэффициента положительной обратной связи у = 0,8-Ю-18 К/мг в климатической системе количество метана в атмосфере с учетом эмиссии из мерзлотных озер будет увеличиваться в ближайшие триста лет с большей скоростью, по сравнению с ситуацией, когда метан из озер практически не выделяется (рис. 2).
Таким образом, в данной статье предлагается модель динамики радиуса мерзлотного озера, которая может быть использована для прогнозирования эмиссии метана в атмосферу при таянии вечной мерзлоты. Представленная модель является аналитической и не требует сложных численных расчетов. Она базируется на методах теории фазовых переходов Гинзбурга — Ландау и асимптотической теории радиуса средней кривизны, т. е. основывается на фундаментальных физических теориях, что делает ее более надежной, по сравнению с простыми численными моделями протаивания мерзлотных озер. Кроме того, в рассматриваемой модели используется уравнение Фоккера — Планка, позволяющее учитывать случайные эффекты.
В результате последовательного применения аналитических методов получено явное соотношение, описывающие динамику радиуса озера, которое в дальнейшем использовано для расчетов потоков метана из мерзлотного озера. Показано, что это соотношение согласуется с имеющимися данными наблюдений и простыми эмпирическими моделями, предложенными ранее. Следует отметить, что в полученное соотношение входит ряд эмпирических параметров, которые могут быть найдены из данных наблюдений или эксперимента, поэтому модель может использоваться для анализа таких параметров, получаемых при мониторинге состава атмосферы. В статье рассмотрено использование модели для описания положительной обратной связи в климатической системе. Из представленных расчетов следует, что резкое увеличение концентрации метана в атмосфере может произойти в ближайшие триста лет.
СПИСОК ЛИТЕРАТУРЫ
1. Комаров, И.А. Термодинамика и тепломассообмен в дисперсных мерзлых породах [Текст] / И.А. Комаров,— М.: Научный мир, 2003,— 608 с.
2. Caginalp, G. Stefan and Hele — Shaw type problems as asymptotic limits of the phase field equations |Text| / G. Caginalp // Phys. Rev. A.- 1989,- Vol. 39,- P 128-146.
3. Angenent, S. A computed example of non-uniqueness of mean curvature [Text] / S. Angenent, D.L. Chopp, T. llmanen // Comm. Partial Differential Equations.- 1995,- Vol. 20,- № 11-12,- P. 1937-1958.
4. Angenent, S. Some recent results on mean curvature flow | Text| / S. Angenent // RAM Res. Appl. Math. -1994,- Vol. 30,- P. 1-18.
5. Ecker, K. Mean curvature evolution of entire graphs |Text| / K. Ecker, G. Huisken // Ann. of Math. -1989,- Vol. 130,- № 3,- P. 453-71.
6. Ecker, K. Interior estimates for hyper surfaces moving by mean curvature |Text| / K. Ecker, G. Huisken // Invent. Math.- 1991,- Vol. 105,- № 3,-Р. 547-569.
7. Молотков, И.А. Сосредоточенные нелинейные волны [Текст] / И.А. Молотков, С.А. Вакулен-ко,— J1.: Изд-во Ленингр. ун-та, 1988,— 240 с.
8. Gage, М. The heat equation shrinking convex plane curves |Text| / M. Gage, R. S. Hamilton // J. Differential Geom.- 1986. -Vol. 23,- № 1,- P. 69-96.
9. De Mottoni, P. Geometrical evolution of developed interfaces [Text] / P. de Mottoni, M. Schatzman // Trans. Amer. Math. Soc.— 1995,— Vol. 347,- № 5,- P. 1533-1589.
10. Downing, J.A. The global abundance and size
distribution of lakes, ponds, and impoundments [Text] / J.A. Downing, Y.T. Prairie, J.J. Cole, |et al.| // Eimnol. Oceanogr.- 2006,- Vol. 51,- № 5,- P. 2388-2397.
11. Matthews, E. Methane emissions from natural wetlands: global distribution, area, and environmental characteristics of sources [Text] / E. Matthews, l.Y. Fung // Global Biogeochemical Cycles.— 1987.— Vol. 1,- P. 61-86.
12. Frolking, G. Climate controls on temporal variability of methane flux from a poor fen in southeastern New Hampshire: Measurement and modeling [Text] / G. Frolking, P. Grill //Global Biogeochem. Cycles.— 1994,- Vol. 8,- P. 385-387.
13. Cao, M.K. Global carbon exchange and methane emissions from natural wetlands: Application of a process-based model |Text| / M.K. Cao, S. Marshall, K. Gregson // J. Geophys. Res.- 1996,- Vol. 101.— R 14399-14414.
14. Kerr, R. Arctic 'Armageddon' needs more science less hype |Text| / R. Kerr // Science.— 2010.— Vol. 8,- P. 620-621.
15. Kiehl, J.T. Earth's annual global mean energy budget | Text | / J.T. Kiehl, K.E. Trenberth //Bulletin of the Am. Meteorol. Soc.- 1997,- Vol. 78,- No. 2,-P. 197-208.
