Аналитическое исследование малоразмерной модели динамики углерода в биосфере Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 17, № 3, 2012

Аналитическое исследование малоразмерной модели динамики углерода в биосфере*

А. М. Федотов1, С. Б. Медведев1, А. И. Пестунов1, И. А. Пестунов1, С. И. Барцев2, А. Г. Дегерменджи2 1 Институт вычислительных технологий СО РАН, Новосибирск, Россия 2Институт биофизики СО РАН, Красноярск, Россия e-mail: fedotov@ict.nsc.ru, medvedev@ict.nsc.ru, alexpi@ngs.ru, pestunov@ict.nsc.ru, bartsev@yandex.ru, ibp@ibp.ru

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

Ключевые слова: круговорот углерода, парниковый эффект, малоразмерная модель динамики биосферы, устойчивость по Ляпунову, триггерный режим, прогнозные сценарии.

Введение

Проблема круговорота углерода в биосфере привлекает внимание исследователей ввиду наблюдаемых в последние десятилетия изменений климата [1] и существующей гипотезы "парникового эффекта" [2]. Вряд ли в ближайшее время удастся доказать или опровергнуть эту гипотезу, поскольку нельзя поставить глобальный эксперимент над биосферой. Даже если такой эксперимент удастся провести, то он будет весьма дорог и его последствия могут быть катастрофическими. Более подробно аргументы за и против парникового эффекта приведены в работе [3]. Тем не менее антропогенное воздействие на биосферу продолжается, и остаётся лишь руководствоваться принципом предосторожности [2] и прибегать к математическому моделированию как фактически к единственному способу изучения глобального круговорота углерода.

Математические модели глобальной динамики углерода развиваются в двух основных направлениях [4]: 1 — усложнение моделей за счёт ввода дополнительных уравнений и параметров, что предположительно повышает адекватность описания, но усложняет их анализ, 2 — упрощение моделей, т. е. построение так называемых минимальных или малоразмерных моделей [5], использующих небольшое число интегральных переменных. Минимальные модели не претендуют на количественное описание динамики

* Работа выполнена в рамках Интеграционного проекта СО РАН № 50 "Модели изменения биосферы на основе баланса углерода (по натурным и спутниковым данным и с учётом вклада бореальных экосистем)" при поддержке президентской программы "Ведущие научные школы РФ" (грант № НШ-6293.2012.9).

системы биосфера — климат, определяемой совокупным действием множества известных климатических и биосферных процессов, однако исследование таких моделей позволяет понять качественные механизмы основных процессов, происходящих в биосфере, и оценить их возможные последствия. В данной работе исследуется малоразмерная модель, представленная в [6].

1. Малоразмерная модель динамики биосферы

Модель [6] при условии отсутствия антропогенных выбросов углерода в атмосферу можно записать в виде системы трёх уравнений

dx

-37 = Vxx(xmax - x)V(A)/mg(T(A)) - VdX, dt

dt = Vdx - VSyfMD(T(A)), (1)

C = A + x + y.

Здесь /MG(T) = Td(TMG — T)0(T)0(TMG — T) (MG — MaxGrowth) — зависимость скорости прироста биомассы растений от среднегодовой глобальной температуры, /md(T) = Td(TMD — T)0(T)0(TMD — T) (MD — MaxDecay) — зависимость скорости почвенного дыхания от температуры, 0(x) — функция Хевисайда

^ I 1 при x > 0,

0(x) = <

10 при x < 0.

Зависимость вида Моно представляет собой функцию от переменной углерода атмо-A

сферы —-- и описывает скорость прироста биомассы растений. Первое уравнение

Kо + A

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

A

круговороте углерода. Эта зависимость имеет вид T(A) = T0 + Tdel log2 -г-. Параметр

A0

Tdei показывает число градусов Цельсия, на которое увеличится глобальная температура при удвоении углерода в атмосфере. В случае отсутствия парникового эффекта глобальная температура предполагается постоянной: T(A) = T0. Сведения о параметрах модели представлены в табл. 1. Термин "глобальная среднегодовая температура поверхности" используется межправительственной группой экспертов по изменению климата (Intergovernmental Panel on Climate Change — IPCC) для описания средней по Земному шару температуры.

Из закона сохранения массы углерода следует, что количество углерода в атмосфере A = C — x — у. Подставим это выражение в первые два уравнения системы (1). Учтём начальные данные: C = A0 + x0 + у0. Рассмотрим полученную систему как систему двух дифференциальных уравнений от переменных x и у. Относительно стационарных

Таблица 1. Оценки параметров системы

Параметр Минимально Максимально

возможное возможное Описание

значение значение

А0 600 760 Количество углерода в атмосфере, последние десятилетия, Гт

Х0 500 850 Количество углерода в биомассе растений, последние десятилетия, Гт

У0 1080 2011 Количество углерода гумуса, последние десятилетия, Гт

( 1.4 1.6 Показатель степени, приближающий вид температурной зависимости к более реалистичному

Тмс 30 45 Температура прекращения роста биомассы растений Т, °С

Т0 13 17 Глобальная среднегодовая температура поверхности в настоящее время Т, °С

Тае1 2 11 Чувствительность климата Т, °С

К0 837 1023 Параметр в зависимости вида Моно

Р0 50 60 Стационарная скорость круговорота углерода, Гт/год

TMD 35 55 Максимальная температура распада органики Т, °С

С 1.1 1.6 Коэффициент ёмкости среды, Сх0 = хтах

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

¿х

— = Р (х,у),

¿у _ (2)

■ И = Q(x'y)■

Существование начального стационарного состояния системы определяется выбором масштабных множителей Ух, Д/^, V в соответствии с условиями Уххо(хтах — х0)У(А0)х /мс(Т(А0)) = Уа,х0 = Vу0/Мв(Т(А0)) = Р0, обеспечивающими равные скорости прироста чистой первичной продукции (ЧПП), отмирания биомассы растений и интенсивности почвенного дыхания.

2. Анализ малоразмерной модели

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

2800

н и

и о

>.2000

га Ч О Рч

и

0

03

1

с; о

1200 -

400-

Изоклиг горизон га талей Иг юкли на

вертикалей

0 200 400 600 800

Количество углерода в биомассе растений, Гт

£ 2800

и и >.

>.2000 -

га

се «

о 1200 -

Сц

и

Ч -

>. О

п

ч

о

400

II 1 : 1 1

11 : )

1 -- 1 ^ Г 1 1 г г

1 -

Изоклина

горизонталей ___________ 1 1 Изоклина

_______— \ 1 вертикален

0 200 400 600 800

Количество углерода в биомассе растений, Гт

а

Рис. 1. Главные изоклины системы при следующих параметрах: а — Ао = 680, Хо = 675, уо = 1546, й = 1.5, Тис = 38, То = 15, ТЛе1 = 7.5, Ко = 930, Ро = 55, Тмв = 43, Хтах = 912; б — Ао = 600, хо = 650, уо = 1100, й = 1.5, Тмс = 30, То = 15, ТЛеЛ = 8, Ко = 930, Ро = 60, ТмБ = 55, хтах = 715. Вставка показывает детальное устройство фазовой плоскости

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

Замечание 1. Стационарная точка (х = 0, у = 0) названа точкой деградации, поскольку она соответствует исчезновению растительности и почвы.

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

2.1. Возможность реализации триггерного режима в биосфере

Кратко результаты анализа, представленного в данном разделе, описаны в [8]. Проведём анализ устойчивости по Ляпунову начального стационарного состояния.

Утверждение 2 (о характеристическом уравнении в начальной стационарной точке модели с парниковым эффектом). Характеристическое уравнение в точке (х = хо, у = уо) является приведённым квадратным уравнением А2 + рА + q = 0, где

= / 1 + Ко___ТЛе1

\ хтах хо Ао (Ко + Ао) Ао (Тмс - То) 1п 2 уо

Т4в1

+

Ао (Тмв - То) 1п 2

Р 2

РП

1

+

Ко

+

¿г

<1е1

(Хтах - Хо) уо Аоуо (Ко + Ао) АоУоТо 1п 2

т

(е1

¿Т-

(е1

Аоуо (Тмо - То) 1п 2 Ао ( Хтах Хо) То 1п 2

(е1

+

Ко

Ао ( Хтах - Хо) (Тмб - То) 1п 2 АоХо (Ко + Ао)

+

ЛТ-

(1е1

Т

(е1

АоХоТо 1п 2 АоХо (Тмо - То) 1п 2^

В соответствии с критерием Рауса — Гурвица [9] оба корня характеристического уравнения имеют отрицательную действительную часть тогда и только тогда, когда р > 0 и д > 0. Эта ситуация по теореме Ляпунова соответствует устойчивому случаю. Если р < 0 или д < 0, то начальная точка (х = хо, у = уо) неустойчива. При определённых параметрах, попадающих в диапазоны своих оценок, начальное состояние биосферы является неустойчивым по Ляпунову, причём неустойчивая точка имеет тип седло (см. ниже раздел 2.4).

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

$р

дР (х,у)

+

(хо,Уо)

дУ

(хо,Уо)

где £р — след матрицы линейного приближения системы в окрестности особой точки. Неявный вид линии нейтральности для начальной стационарной точки —

Т(е1 (ТМБ - ТМО)

1

1

1

1

- - = 0.

1п2Ао (Тмо - То) (Тмо - То) Ко + Ао Ао х(С - 1) у

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

А

дР (х,у)

(хо,Уо)

ду

(хо,Уо)

дР (Х,у)

ду

(хо,Уо)

(хо,Уо)

д

+

0

0

где А — определитель матрицы линейного приближения системы в окрестности особой точки.

Граница устойчивости находится в области отрицательных значений количества органического вещества в почве (рис. 2). Это означает, что одиночная особая точка всегда устойчива при реальных начальных значениях параметров. При переходе значений параметров через границу моностационарности появляются сразу три особых точки, одна из которых неустойчива, что приводит к возможности триггерных переключений между двумя устойчивыми состояниями. Термин триггерный режим означает, что система может находиться в двух состояниях, между которыми осуществляется относительно быстрый переход.

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

Рис. 2. Параметрический портрет при Ао = 600, хо = 650, уо = 1100, й = 1.5, Тмс = 35, То = 15, Ко = 930, То = 60, Тмб = 55, хтах = 715. Область моностационарности расположена ниже линий моностационарности, область устойчивости — над линией нейтральности

Рис. 3. Два сценария триггерного режима при А0 = 600, х0 = 650, у0 = 1100, й = 1.5, Тмс = 35, То = 15, Ко = 930, То = 60, Тмб = 55, Хтах = 715, ТЛе1 = 8

2.2. Критические параметры, способствующие возникновению триггерного режима

Коэффициенты характеристического уравнения из утверждения 2 — вещественные, поэтому его корни лежат на вещественной прямой и являются комплексно-сопряжённы-

.. + Р /Р2

ми. Решения уравнения даёт формула А± = — ± у —— д. Рассмотрим случай веще-

2 V 4

р2

ственных корней при —— д > 0. Нас интересует пограничная ситуация возникновения неустойчивости, которая может иметь место, когда первое собственное значение корня отрицательно, а второе — равно нулю: А- = —р < 0, А+ = 0 и д = 0, причём д = 0 задаёт некоторое многообразие в пространстве параметров. Обозначим набор параметров, при котором достигается это условие, через Го = (г0,г0, ). Рассмотрим малые возмущения параметра гг = Г0 + Аг* в окрестности точки г0. Разложим выражение для большего собственного значения А+ в ряд Тейлора по линейному приближению:

\ + дА+^Д * Р(г0 К, /Р2(г0) / ч,

А+(го) + "зТ"(го)Аг =--¡Г" + v — д(го) +

+

Аг*

( 1 ( дР дд \

_1ар (г) + 2р(го)дг*(го) — дг*(го)

2дг*( 0) + 2 IР2(го) ( )

v 2у-4--д(го) )

Учитывая, что р(г0) > 0, А+(г0) = 0 и д(г0) = 0, получим

дА+ . 1 Зд

д^<г»)Аг> = — ддд(го)Аг'- (3)

1

Ввиду того, что множитель ——- является величиной постоянной, для сравнения скоро-

р(го)

сти изменения собственного значения А+ в окрестности нулевого значения достаточно

Зд Зд Зд

сравнить производные от д по различным параметрам: ——, ——, ... , ——-. Параметр,

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

Р2

В случае комплексных корней имеем —— д < 0. Условие р(г0) = 0 вытекает из того, что рассматривается граница областей устойчивости, а комплексная часть собственного значения на устойчивость не влияет (т. е. Р = 0 также задаёт некоторое многообразие в пространстве параметров). Результаты численного моделирования методом Монте-Карло (см. ниже раздел 2.4) показали, что случай комплексных корней с положительной вещественной частью нереализуем, поэтому ситуация неустойчивого фокуса невозможна, и в этом случае мы не будем производить разложение корней характеристического уравнения в ряд Тейлора. Отметим также, что если корни вещественны, то они имеют различные знаки (см. раздел 2.4). Следовательно, неустойчивая особая точка является седлом.

Наиболее вероятное положение системы (назовём его типичной точкой) имеет координаты, приведённые в табл. 2. Эти координаты примерно соответствуют серединам

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

q{Aq,Xo, Уо, d, TMG,To,Tdei ,К0, Po,TMD,Xmax) = 0.

Чтобы это расстояние отражало процентное отклонение параметра от своего типичного значения, необходимо перейти к безразмерным переменным. Для этого введём новые

Рг „ ,

переменные хг = —, где рг — исходный параметр, ti — соответствующая координата

ti

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

Далее решалась задача минимизации евклидовой нормы

11

V^ (хг — 1)2 —> min

' J X1,X2,...,X11

г=1

от типичной точки до многообразия Mq = {х : q(x1, х2,..., х11) = 0}. В этом случае уравнение для многообразия можно записать в новых переменных.

\

Таблица 2. Формулировка задачи оптимизации. Результаты оптимизации*

Переменная x1 X2 Х3 x4 X5 X6 X8 X9 X10 X11

оптимизации

Левое 0.882 0.74 0.699 0.8 0.789 0.8 0.267 0.8 0.9 0.814 0.815

ограничение

Правое 1.118 1.259 1.301 1.2 1.184 1.2 1.467 1.2 1.1 1.163 1.184

ограничение

Результат 0.989 1.057 1.015 1.024 0.994 0.972 1.017 0.998 1 1.016 0,932

оптимизации

Исходный Ао xo Уо d TMG To Tdel Ko Po TMD Xmax

параметр

Типичное 680 675 1546 1.5 38 15 7.5 930 55 43 912

значение

Значение 673 713 1569 1.54 37.7 14.6 7.7 928 55 43.6 850

в точке

минимума

Первый показатель 1 5.3 1.2 0.3 0.8 2.7 2.6 0.2 0 1.4 7.3

критичности, %

Второй показатель 0.01 0.065 0.018 0.027 0.007 0.032 0.019 0.002 0 0.018 0.077

критичности

* Подчеркнуты значения безразмерных переменных, соответствующие критическим параметрам.

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

Вторым показателем критичности параметра является модуль производной в

Зд

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

др.

ного значения в ряд Тейлора по линейному приближению (см. (3)).

Результаты оптимизации (см. табл. 2) показывают, что состояние триггерного переключения находится не более чем в 8-процентном изменении параметров от типичной точки. Если рассмотреть данные по критическим параметрам, то оба показателя выделяют два параметра жтах — максимально возможное количество углерода в биомассе растений и ж0 — начальное количество биомассы растений. Кроме того, движение в область неустойчивости (первый критерий) определяется главным образом увеличением ж0 и уменьшением жтах, т. е. снижением коэффициента ёмкости среды С. Этот же параметр критичен с точки зрения перехода границы устойчивости (второй критерий).

2.3. Возможность деградации биосферы при массовых лесных и торфяных пожарах

Утверждение 3 (о собственных значениях в точке деградации модели с парниковым эффектом). Если выполнено условие

Т0 + ТЛе11о§2 0 + л° + У° < тт (Тыс, Тиб) , (4)

Ао

то собственные значения, отвечающие за устойчивость точки деградации (ж = 0, у = 0), записываются в виде

(Ао + ж0 + у0\

Т0 + Т<Ы 1о§2 -1- / „ N

_ _Ао / гт гр гр ! А0 + ж0 + Уо .п

А1 =--/лр Ы ,ГТ-- Т МБ — Т 0 — Т йе1 1о§2 -1- < 0,

Уо(То)" (ТмБ — То) у Ао

Ро Р (Ко + Ао)

жтах

(Ао + жо + уо)

2 Жо Жо (жтах — жо) Ао(То)^ (Тмо — То) (Ко + Ао + жо + уо)

Ао + ж0 + у0\ ( Ао + ж0 + уо х (То + ТЛеЛ 1о§2-а-) I Тмо — То — ТЛеЛ 1og2-а-) . (5)

Замечание 2. Условие (4) необходимо для того, чтобы ступенчатая функция 0(ж) в выражениях для /ыс и /иб была равна единице, иначе рост биомассы растений либо почвенное дыхание прекратятся и мы будем иметь дело с вырожденной системой.

Преобразуем формулу (5) и получим эквивалентное условие устойчивости

Ао + ж0 + у0 \

А2 < 0 ^ (Ко + Ао) жтах (Ао + жо + уо) ( То + ТЛех 1og2 -А- I х

! ч 3600 о о

<В rv

£ » ° я ч о ЬИ

ч я

Рч

о

О

н

О U

з-S

ч о

U .iri'

Я о в

Ч о s

2 Ssh

в< 2 о L

Н О н

Ч S о

О н

а ~ —

н ч

о о о

W rv о

£ о >

3400 3200

200 100

300 -200100

о

50 100 150 200 250

50 100 150 200 250

ч__

------

50 100 150 200 250

Время, годы

Рис. 4. Деградация системы с парниковым эффектом при сильном возмущении, вызванном массовыми лесными и торфяными пожарами. Параметры модели: Ao = 760, xo = 850, yo = 1985, d = 1.5, Tmg = 30.7, To = 15, Tdel = 7, Ko = 930, Po = 60, Tmd = 50, xmax = 1359. Возмущение переводило начальное положение (xo = 850, yo = 1985, Ao = 760) в положение (x = 178, y = 356, A = 3061)

x Tmg - To - Tdel log

Ao + xo + yo

A

- (Xmax - xo) Ao(To)d (Tmg - To) (Ko + Ao + xo + yo) < 0.

(6)

Численно найдены значения параметров, при которых условие (6) выполнено. Таким образом, рассматриваемая точка является притягивающей (рис. 4).

2

2.4. Оценка условной вероятности реализации сценариев методом Монте-Карло

Если предположить, что вероятность принять какое-либо значение из области изменения параметра имеет равномерное распределение в этом диапазоне, то можно смоделировать вероятность того или иного сценария, например, вероятность возникновения триггерного режима. Пусть 0 < A < B. Выберем p = A + random (B - A), где random — псевдослучайная величина с равномерным распределением на интервале (0,1). Тогда p имеет равномерное распределение на интервале (A,B). Каждый из одиннадцати рассматриваемых параметров разыгрывался независимо. В результате получался набор, задающий конкретную модель, что соответствовало одному испытанию. Было проведено более 107 испытаний для каждого из сценариев.

После задания набора параметров рассматривалась вещественная часть собственных значений характеристического уравнения, отвечающего за устойчивость начального состояния из утверждения 2. Вероятность неустойчивости составляет 12%, причём ситуация комплексных корней с положительной вещественной частью является нереализуемой. Кроме того, если начальная стационарная точка неустойчива, то она является седлом (оба корня действительны и имеют разные знаки). Вероятность неустойчивости существенно уменьшается (до 4 %), если выбрать максимальное значение Tdei, равное 8. Условие (4) выполнено в 87 % испытаний, а одновременное выполнение условий (4)

и (6) происходит в 2.3% испытаний, т.е. ситуация, когда точка деградации устойчива, достаточно редка. Кроме того при Т^ < 5 °С точка деградации неустойчива.

Замечание 3. В действительности параметры не распределены равномерно по интервалам, а имеют,, вероятно, большую плотность около середин интервалов. Так, поскольку при задании набора на середине интервалов оценок начальное стационарное состояние устойчиво, вероятность триггерного режима меньше полученной оценки. Меньше данной оценки и ситуация устойчивости точки деградации, так как при типичных значениях параметров эта точка неустойчива.

2.5. Различия в свойствах моделей с наличием и отсутствием парникового эффекта

Утверждение 4 (о модели без парникового эффекта). В модели без парникового эффекта всегда существуют три стационарных состояния. Первая точка — неустойчивые деградации (х = 0, у = 0), вторая — устойчивое начальное положение системы (х = хо, у = уо), третья — недостижимое без нарушения закона сохранения массы углерода положение (х = х^3, у = у3ы), где

КоХтах + Коуо х^ + Ко А + (Л)2 + Лхо + Л уо х0

(Ко + А) ( 1 + -х0

хтах

Кохтахуо + Ко^о)2^^ + КЛуо + (А)2уо + Ахоуо + Л (уо)2 уш = °(Ко + А) (хо + уо) . (8)

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

Заключение

Анализ системы (1) показал возможность реализации в биосфере триггерного режима, способного привести к кардинальному изменению состояния биосферы даже без дополнительного сжигания ископаемых топлив. Этот режим возможен при значениях параметров биосферы, входящих в интервалы их существующих оценок. Отсюда следует потенциальная опасность любых резких изменений состояния биосферы, в том числе вызываемых хозяйственной деятельностью человека. Так, уменьшение биомассы в пересчете на углерод (например, в результате вырубки лесов) на 4 ГтС/год в течение четырёх лет приводит к последующему повышению среднегодовой температуры поверхности до 17 °С, а равное по величине, но противоположное по знаку воздействие (например, массовые лесопосадки) вызывает понижение температуры до 5 °С. При этом процесс перехода биосферы в новое состояние происходит примерно в течение 100 лет. Переход в состояние неустойчивости из типичной точки осуществляется

главным образом при уменьшении ёмкости среды G. Вероятность триггерного режима составляет 12% (при равномерном распределении параметров по диапазонам).

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

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

Список литературы

[1] IPCC. Climate Change: The Physical Science Basis. Summary for Policymakers. Fourth Assessment Report. Intergovernmental Panel on Climatic Change. Geneva, Switzerland, 2007.

[2] Тарко А.М. Можем ли мы затормозить глобальное потепление? Россия в окружающем мире — 2008. Устойчивое развитие: Экология, политика, экономика. Аналитический ежегодник / Отв. ред. Н.Н. Марфенин; под общей редакцией Н.Н. Марфенина, С.А. Степанова. М.: МНЭПУ, 2008. 328 с.

[3] Федотов А.М., Медведев С.Б., Пестунов А.И., Пестунов И.А. О нестандартном поведении минимальной модели углеродного цикла // Вестник НГУ. Информационные технологии. 2011. Т. 9, вып. 1. С. 82-88.

[4] Кондратьев К.Я., Крапивин В.Ф. Моделирование глобального круговорота углерода. М.: Физматлит, 2004. 336 с.

[5] Svirezhev Y.M., von Bloh W. A minimal model of interaction between climate and vegetation: qualitative approach // Ecolog. Model. 1996. Vol. 92, is. 1. P. 89-99.

[6] Барцев С.И., Дегерменджи А.Г., Ерохин Д.В. Глобальная минимальная модель многолетней динамики углерода в биосфере // Докл. АН. Геофизика. 2005. Т. 401(2). С. 233-237.

[7] ИвАницкий RP., Кринский В.И., Сельков Е.Е. Математическая биофизика клетки. М.: Наука, 1978. 308 с.

[8] Барцев С.И., Дегерменджи А.Г., Федотов А.М. и др. Биосферный триггер в минимальной модели углеродного цикла // Докл. АН. География. 2012. Т. 443(4). С. 500-503.

[9] Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). М.: Наука, 1973. 832 с.

Приложение (доказательства утверждений)

Доказательство утверждения 1 (о стационарных состояниях модели с парниковым эффектом). Нетрудно видеть, что точка (x = 0, y = 0) обращает в ноль

оба уравнения системы системы (1), т. е. является стационарной. Точка (ж = жо, у = уо) в силу выбора масштабных множителей

Ро(Ко + Ао) Ро Ро

К =-—-¡г-, УЛ = —, V =-- (9)

Хо(жтах - Хо)Ао (То) (Тмо - То) жо уо (То) (Тмб - То)

также стационарна. Для проверки этого подставим условия (9), А = Ао, Т(А = Ао) = Ао

То + Т^ log2 —- = То в систему дифференциальных уравнений и также получим, что Ао

аж ^у

-Г= -г = 0. М М

Доказательство утверждения 2 (о характеристическом уравнении в начальной стационарной точке модели с парниковым эффектом). Учтём, что в стационарной точке (ж = жо, у = уо)

0(То) = 1, ©(Тмо - То) = 1, 0(Тмв - То) = 1

10)

Используя условия (10), получим систему, описывающую динамику углеродного круговорота в окрестности стационарного состояния:

аж ^ж(жтах - ж)А / А \ / А

Я = Ко + А + ^ log2 А^ (^Тмо - То - ^ log2 Ао • - ^

ау А ¡ А

— = V«ж - У8у I То + Tde1 log2 I ТмБ - То - Тм log2 А

Подставим выражение для А из закона сохранения массы углерода в первое и второе уравнения этой системы и вычислим матрицу частных производных:

аж „ ч ау /

Л = Р<ж-у»- i = « м

^ 5Р(ж, у) 5Р(ж, у) \

5ж 5у

д^(ж,у) д^(ж,у)

\ 5ж

5у /

Для вычисления

5Р (ж, у) 5ж

учтём, что

(/1/2/3/4/5)' = (Л)' /2/3/4/5 + /1 (/2)' /3/4/5 + /1/2 (/а)' /4/5 + /1/2/3 (/4)' /5 + /1/2/3/4 (/5)'

дж 5

(/1)' = =1, (/2)' = дж (жтах - ж) = -1,

5

(/3)' = ТТ"

Ао - жо - уо - ж - у

-Ко

5ж \ Ко + Ао + жо + уо - ж - у / (Ко + Ао + жо + уо - ж - у)

7Х

х

5

дж

То + Т^ log2

Ао + жо + уо - ж - у

ао

¡

,т ! Ao + x0 + yo — x — у\ d ^

d ( To + Tdel log2 -A- ) Tdel

(Ao + Xo + У0 — X — y) ln 2

д Ao + xo + yo — x — y Tdel

"TT" TMG — To — Tdel log2

д^ \ Ao / (Ao + x0 + y0 — x — y) ln 2

д—(x,y) ч Ao + x0 + y0 — x — y ,

( ,У) = — Vd + Vx(Xmax — x)—0-y-Td(TMG — T) —

дx Ko + Ao + xo + yo — x — y

Ao + xo + yo x y

— VxT ° . 0 У0- T d(TMG — T ) —

Ko + Ao + xo + yo — x — y

— VxX(Xmax — x) —---0--¿T d(TMG — T ) —

(Ko + Ao + Xo + yo — x — y)

1 dT¿ T d— l

— Vxx(xmax — x) ^ , , .-¡--f^-(TMG — T) +

Ko + Ao + xo + yo — x — y ln 2

1 d Tdel

+VXx(xmax — x)^ : i : : T ï ^.

Ko + Ao + xo + yo — x — y ln 2

Аналогично

д—(x,y)

ду Vxx(xmax x) Х

x

— KoTd(TMG — T ) — dTdel Td—1(Tmg — T ) + TdT

MG — T ) + — dT delT (T MG — T ) + T T del

(Ko + Ao + Xo + yo — X — y)2 (Ko + Ao + xo + yo — x — y) ln 2

d-1

дО^, y)

Ao + xo + yo — x — y 1 d ( To + Tdel log2 -A- ) Tdel I — ^

Ao + xo + yo — x — y

-;-ln 2

Ao

Ao + xo + yo — x — y Х ( Tmd — To — Tdel log2 -A- I +

d — Tdel I--1-

+ To + Tdel log2

Ao + xo + yo — x — y Ao

Ao Ao + xo + yo — x — y

A

ln 2

o

+ (Ao + x, x — y) ln 2 <dTd—1 (Tmd — T) — Td)

(glg2g3)/ = (gl)/g2g3 + gl (g2)/g3 + glg2 (g3)/ ,

/ \ d , V дУ -, rp ^rp ! A0 + x0 + yo — x — y

(gl) = TT = 1, g2 = T0 + Tdel log2

ду - i — —-2 a,

x

o

#3

ТмБ - То - Тйе1 log2

Ао + жо + уо - ж - у

ао

В итоге получим

д^(ж,у) 5у

-V,

С - ж - у\ / С - ж - у' То + Ти log2 --:- ТмБ - То - Т^ log2 -;- | +

уЦ То + Ти log2

А0

Ао + жо + уо - ж - у

ао

А0

¡1

т(ы '- ао

-

Ао + жо + уо - ж - у

-х

А0

х ТмБ - То - Ти log2

1п 2

Ао + жо + уо - ж - у

АО

-

+у То + Те log2

Ао + жо + уо - ж - у

АО

Т(И 1 - АО

Ао + жо + уо - ж - у

А0

1п 2

) \

К

-Тл (Тмв - Т) +

уаТ ¡-1Т«е (С - ж - у) 1п 2

(ТмБ - Т) -

уТ «¡Те

(С - ж - у) 1п 2

Подставим (ж = ж0, у = у0) в компоненты матрицы частных производных М:

5Р (ж,у)

5ж

Ро

РоКо

Р0аТ«е1

РоТ

0Т4е1

(хо,Уо)

5Р (ж,у)

жтах - жо Ао (Ко + Ао) АоТо 1п 2 Ао (Тмо - То) 1п 2

5у

(хо,Уо)

5^(ж,у)

5ж

5^(ж,у)

-

det

5Р (ж,у)

5ж

5у

М 1(хо,Уо)

Ро Ао

(хо,уо)

= Ро

Ко

¡е1

-

Т

¡е1

Ко + Ао То 1п 2 (Тмо - То)1п2

РО + РоТ«е1 (±__

жо Ао 1п 2 ^ То ТмБ - То

аТае1 Тйе1 1

(хо,уо)

- АЕ

тоао 1п 2 ао (тмв - то) 1п 2 уо

А2 А

5Р (ж,у)

5ж

-

5^(ж,у)

(хо,Уо)

5у

5Р (ж,у)

(хо,Уо)

5у

(хо,уо)

5^(ж,у)

(хо ,Уо)

5ж

5у

(хо ,Уо),

-

(хо,уо),

А2 + рА + д = 0.

«

«

После преобразований и сокращений получим искомые выражения для р и д.

Доказательство утверждения 3 (о собственных значениях в точке деградации модели с парниковым эффектом). Подставим (х = 0, у = 0) в компоненты матрицы частных производных М, найденные в утверждении 2:

дР (х,у)

Зх

Ро (Ко + Ло) хтах (Ао + хо + уо)

(о,о)

Ро

---1--

хо хо (хтах - хо) Л(То)а (Тмо - То) (Ко + Ао + хо + уо)

X

х То + ТЛе1 1с§

Ао + хо + уо

ао

ТМО — То — Тйе1 1с®

2

Ао + хо + уо

дР (х,у)

Зу

зд(х,у)

(о,о)

Зх

(хо ,Уо)

Ро

хо'

Зу

Ро То + Ти lcg.

Ао + хо + уо

ао

(о,о)

уо(То)^ (Тмв - То)

ТМБ — То — Тйе1 lcg:

Ао + хо + уо

ао

|м(о,о) - ле | = л2 -л

дР (х,у)

Зх

+

(о,о)

Зу

дР (х,у)

(о,о)

Зх

(о,о)

Зу

0.

(о,о)

Решая последнее уравнение, найдём искомые выражения для собственных значений:

Л-1

зд(х,у)

Зу

< 0, Л

дР (х,у)

2=

(о,о)

Зх

(о,о)

Доказательство утверждения 4 (о модели без парникового эффекта). При

рассмотрении этого доказательства следует помнить, что все выкладки проводятся при условии Т^е1 = 0, что соответствует постоянной глобальной температуре То.

Существование трёх стационарных состояний. Для нахождения стационарных состояний подставим выражения для масштабных множителей (формула (9) в доказательстве утверждения 1) в систему (1) и учтём условия (10) в доказательстве утверждения 2. Произведя сокращения, получим систему

х (Ко + Ао) (хтах - х) (Ао + хо + уо - х - у) = 0

ч (хтах - хо) Ао (Ко + Ао + хо + уо - х - у) ,

ху

--— = 0.

хо уо

В этом случае имеем решение (х = 0, у = 0), т. е. получим первую стационарную точку. хуо

Выразим у = -из второго уравнения и учтём, что (Ко + Ао + хо + уо - х - у) > 0.

хо

Получим уравнение

(Ко + Ао) (хтах - х) ( С - х--—

хо

хуо

(хтах - хо) Ао Ко + С - х--

хо

а

2

0

а

2

где С = Ао + ж0 + у0. Вспомним, что масштабные множители выбирались так, чтобы начальное состояние было стационарным, следовательно, (ж = ж0) есть одно из решений этого уравнения. Заметим, что данное уравнение квадратное и один корень известен. Для нахождения третьего решения продолжим преобразования, вынося множитель (ж0 — ж):

( Жтах Жуо (Ко + Ао) АоЖтах + ЖоЖтах + уоЖтах — ЖЖтах---АоЖо — ЖоЖ — уоЖ+

Жо

(

Ж Уо

+Ж +--

Жо

КоЖтахАо + (Ао)2 Жтах + Ж0Жтах А0 + УоЖтахАо ЖЖтах А0

Жу0ЖтахА0

Жо

+

+ (—КоЖоАо — (Ао)2Жо — УоЖоАо + жжоАо + жуоАо) ^

^ КоЖтах(жо — ж) +--0У0 тах(жо — ж) + КоАо(жо — ж) — Кож(жо — ж) + (Ао)2(жо — ж) —

Жо

Коуож Аоуо ж

(ж0 — ж) + А0ж0(ж0 — ж)--(ж0 — ж) — А0ж(ж0 — ж) + Аоуо(жо — ж) = 0.

Жо

Жо

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

Устойчивость начального состояния. Для получения характеристического уравнения в точке (ж = Жо, у = у0) модели без парникового эффекта достаточно подставить значение = 0 в уравнение из утверждения 2. Получим

Л2 + Ар + д = 0,

где

р = —о

1

Ко

(Жтах — Жо) Ао(Ко + Ао) уо

> 0,

д = —о2

1

Ко

+

+

Ко

> 0.

(Жтах — Жо)уо Аоуо(Ко + Ао) АоЖо(Ко + Ао)

В силу положительности р и д по критерию Рауса — Гурвица и теореме Ляпунова система всегда устойчива в начальной стационарной точке при реальных значениях параметров.

Неустойчивость точки деградации. Матрица частных производных в точке (0, 0)

( ( РоЖтах (Ко + Ао) (Ао + Жо + уо) _ —о

Жо (Жтах — Жо) Ао (Ко + Ао + Жо + уо) Жо —

Жо

А1 = — — < 0, у0

\

0

—

уо )

PoXmax (Ко + Ao) (Ao + Xo + yo) Po

Ao

Хо (хтах - Хо) Ао (Ко + Ао + Хо + Уо) Хо _ Ро (Хтах (Ко + Ао) (Ао + Хо + уо) - (Хтах - Хо) Ао (Ко + Ао + Хо + уо)) Хо (Хтах - Хо) Ао (Ко + Ао + Хо + Уо) '

Чтобы доказать положительность корня Л2 (Л2 > 0), достаточно показать, что

Хтах (Ко + Ао) (Ао + Хо + Уо) - (Хтах - Хо) Ао (Ко + Ао + Хо + Уо) > 0 ^

^ (КоХтах + АоХтах) (Ао + Хо + уо) - (АоХтах - АоХо) (Ко + Ао + Хо + Уо) > 0 ^ ^ КоХтах Ао + КоХтахХо + КоХтах Уо + (Ао) Хтах + АоХтахХо + АоХтах Уо-

АоХтах Ко - (Ао)2Хтах - АоХтахХо - АоХтахУо + АоХо Ко + (Ао)2Хо +

+Ао(Хо)2 + АоХоУо > 0 ^

^ КоХтахХо + КоХтахУо + АоХоКо + (Ао)2Хо + Ао(Хо)2 + АоХоУо > 0.

Последнее условие выполнено при любых реальных значениях параметров, следовательно, Л2 > 0 и по теореме Ляпунова стационарное состояние неустойчиво.

Недостижимость состояния (х _ х^з, У _ Уш) без нарушения закона сохранения массы углерода. Вычтем из углерода биомассы растений и гумуса в третьей стационарной точке (формулы (7), (8) в утверждении 4) весь углерод системы: С _ Ао + хо + Уо:

x

KoXmax + Koyo—max + KoAo + (Ao)2 + AoXo + Ao yo

Xst3 + yst3 - C =-—

(Ko + Ao) ( 1 + yo Xo

Xmax

KoXmaxyo + ^(yo)2^^ + KoAoyo + (Ao)2yo + Ao Xoyo + Ao(yo)2

Xo__C =

(Ko + Ao) (Xo + yo)

2 Xmax

KoXo (Xmax - Xo) + 2Koyo (Xmax - Xo) + Ko (yo)--1

Xo

'11)

(Ко + Ао) (Хо + Уо)

Учитывая, что Хтах > хо, выражение (11) положительно, следовательно, в третьей стационарной точке закон сохранения массы углерода нарушен.

Поступила в 'редакцию 15 марта 2012 г., с доработки — 27 апреля 2012 г.