Моделирование и картирование влияния изменений климата на многолетнюю мерзлоту в регионе со сложным рельефом с высоким пространственным разрешением

ВВЕДЕНИЕ

Потепление климата в высоких широтах в течение XX века было примерно в два раза выше, чем в среднем по миру  [1]. Наблюдения показали увеличение температуры приповерхностного слоя грунта и толщины слоя сезонного протаивания/промерзания (деятельного, активного, сезонно-талого слоя), а в некоторых местах и исчезновение многолетней мерзлоты (ММ) (например, [2, 3]). Большинство климатических моделей прогнозирует, что потепление климата в высоких северных широтах будет продолжаться со скоростью, превышающей среднюю мировую, в течение XXI века [1], что вызовет еще бОльшую деградацию ММ. Таяние ММ влияет на инфраструктуру, экосистемы, места обитания диких животных и имеет сильные обратные связи с климатической системой [1].

Эффективным подходом к пониманию распределения многолетней мерзлоты и его изменений с изменениями климата является пространственное моделирование. Пространственные модели ММ можно в целом разделить на равновесные и переходные [4, 5]. Из-за медленного реагирования термических режимов грунта на изменения климата состояние многолетней мерзлоты находится и будет продолжать находиться в XXI веке в неравновесии с климатом [6, 7]. Поэтому необходимо количественно оценивать изменения ММ на основе переходных моделей. Большинство исследований по переходному моделированию и картированию многолетней мерзлоты проводилось с использованием пространственного разрешения в полградуса по широте и долготе или хуже (например, [7–10]). Исследования с грубым разрешением не могут точно учитывать влияние топографии на инсоляцию, поскольку топографические условия обычно меняются на более коротких расстояниях. Кроме того, результаты таких исследований трудно проверить путем сравнения с данными полевых наблюдений и они не подходят для планирования освоения земель и экологических оценок.

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

Так, авторы статьи [11] составили карту термических условий грунтов и их изменений с изменениями климата на Аляске с разрешением 2 км, используя неявную конечноразностную численную модель. Влияние снежного покрова учитывалось в явной форме, а влияние топографии на инсоляцию не учитывалось. Авторы доклада [12] составили карту многолетней мерзлоты и ее изменений с климатом в бассейне реки Маккензи с разрешением 30 м на основе модели процессов теплопроводности. В этой модели использовались сезонные n-факторы для оценки приповерхностной температуры грунта по температуре воздуха.

Авторы статей [13, 14] закартировали ММ северо-западной части Низменности Гудзонова залива (Hudson Bay Lowlands) с разрешением 30 м, для чего использовалась более подробная модель, в которой были интегрированы эффекты других климатических переменных (например, количества атмосферных осадков, солнечной радиации и давления пара) и изменений в снежном покрове и влажности грунтов. Данный регион в основном представляет собой равнину, сложенную с поверхности торфом, где толщина торфяного слоя может быть оценена на основе высотных отметок.

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

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

МЕТОДЫ И ДАННЫЕ

Исследуемая территория и источники полевых данных

Исследуемая территория – это национальный парк «Иввавик» (НПИ), расположенный в северной части провинции Юкон (территории на северо-западе Канады) (рис. 1). Этот парк занимает площадь 10,17 тыс. км² и состоит из двух основных зон: высокогорного внутреннего региона, который относится к Британским горам, и прибрежной равнины у моря Бофорта. Самая высокая гора достигает высоты 1655 м. НПИ характеризуется необычайным природным разнообразием и включает арктическую тундру, горную тундру, лесные массивы и арктические прибрежные местообитания. БОльшая часть территории парка не была покрыта льдом, поэтому содержит древние формы рельефа, которые за миллионы лет эволюционировали до нынешнего вида. Парк является частью района отёла и летнего пребывания стада поркьюпайнских карибу (северных оленей), насчитывающего около 165 тыс. животных [15]. В настоящее время в НПИ распространена сплошная многолетняя мерзлота [16].

Рис. 1. Карта рельефа национального парка «Иввавик» (а) и расположение парка (серый цвет) в провинции Юкон (б). Указаны точки размещения климатических станций (зеленые квадратики) и мест наблюдений за глубиной летнего протаивания (красные кружки), данные которых использовались в исследовании. Красный кружок на рисунке «б» – место наблюдений на озере Иллисарвик на острове Ричардс

Для представляемого в данной статье исследования были использованы три основных источника полевых данных. Первый – это отчет службы «Парки Канады» о ресурсах НПИ [15], в котором предоставлена общая информация о климате, геологии, геоморфологии, грунтах, гидрологии, гидрогеологии, растительности и дикой природе парка. Второй источник – детальная классификация экосистем вдоль бассейна реки Ферт, который охватывает около трети НПИ  [17]. На основе этой классификации авторы расширили карту экосистем на весь парк и определили грунтовые и дренажные условия для каждого типа. Третий источник данных – полевые наблюдения авторов, проводившиеся в течение четырех летних сезонов с 2008 по 2011 год с определением местоположения, топографии, дренажных и грунтовых условий, состава и состояния растительности, типа экосистемы, глубины летнего протаивания. Последняя измерялась с помощью стальных зондов и выкапывания в грунте шурфов на 162 участках по всему парку, причем на 62 участках из них не был достигнут мерзлый слой из-за высокого содержания в грунте камней. На рисунке 1, а показано распределение 100 участков, на которых удавалось определять глубину летнего протаивания.

Модель

Для расчета и картирования состояния многолетней мерзлоты на территории парка была использована модель Northern Ecosystem Soil Temperature model (NEST – «Модель температуры грунтов северных экосистем»). Это одномерная переходная модель, в которой учитываются воздействия климата, растительности, снежного покрова и грунтовых условий на динамику температурного режима грунтов на основе переноса энергии и массы в системе «грунт – растительность – атмосфера» [18]. Температура грунта рассчитывается путем решения одномерного уравнения теплопроводности. Учитываются динамика глубины снежного покрова, его плотность и их влияние на температуру грунта. Динамика влажности грунта моделируется с учетом поступления воды (из атмосферных осадков и при таянии снега), ее ухода (при испарении и транспирации), а также распределения по слоям грунта. Протаивание и промерзание грунта и связанные с ними изменения в содержании льда и воды определяются на основе закона сохранения энергии. Подробные описания модели и путей ее валидации можно найти в работах [10, 18, 19]. Латеральные водные потоки и перенос снега ветром параметризуются упрощенным образом [13, 20].

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

Входные данные и их обработка

Климатические данные

В пределах или поблизости от НПИ имеется семь метеорологических станций. На трех станциях около северной береговой линии наблюдения ведутся около 50 лет, а на четырех остальных, в том числе на трех удаленных от моря, – менее 15 лет. Δ

В большинстве имеющихся данных есть некоторые пробелы. В нескольких работах представлены разработанные методы пространственной интерполяции среднемесячных климатических показателей на основе наблюдений на станциях (например, [21, 22]).

Авторы настоящей статьи оценили временнЫе закономерности по данным наблюдений на репрезентативных метеорологических станциях, а также усредненное пространственное распределение температуры воздуха и количества атмосферных осадков на основе среднемесячных пространственных данных из работы [21].

Результаты наблюдений показали тесную корреляцию хода температур на разных станциях в НПИ. Например, коэффициент корреляции (R) между среднесуточными температурами воздуха с 1968 по 2010 год на станциях «Комакук-Бич» (69,62 ° с.ш., 140,20° з.д.) и «Олд-Кроу» (67,34° с.ш., 139,50° з.д.) составляет 0,92 (при количестве данных N=15706), хотя эти две станции находятся примерно в 250 км друг от друга (первая – на побережье, а вторая – на внутренней территории (см. рис. 1)). Корреляция между суточными количествами осадков между ними слабая (R=0,16; N=15706). Однако осадки, в сумме выпавшие за месяц на этих станциях, коррелируют (R=0,50; N=516). Поэтому авторы использовали данные наблюдений на климатических станциях для представления как долгосрочных среднемесячных климатических закономерностей, так и ежесуточных колебаний в течение месяца в соответствующих районах.

Среднемесячная температура воздуха в ячейке сетки g в месяце M года Y (T_m,g(Y,M)) оценивалась по следующей формуле:

где T_0m,g(M) – среднемесячная температура воздуха в ячейке сетки g в месяце M, усредненная за 30 лет (1971–2000 гг.), которая была интерполирована на основе данных метеорологических станций из работы [21] с использованием пространственного разрешения 30 м на 30 м (рис. 2); T_m,s(Y,M) – отклонение среднемесячной температуры воздуха в месяце M года Y от среднего значения для этого месяца за 30 лет, оцененного на основе наблюдений на станции.

Рис. 2. Распределение в национальном парке «Иввавик»: а – среднегодовых температур воздуха, усредненных за период 1971–2000 гг. на основе данных из работы [21]); б –общих годовых количеств осадков, усредненных за период 1971–2000 гг. на основе данных из работы [21]); в – индексов листовой поверхности (на основе оценок по снимкам со спутников Landsat); г – типов экосистем (карта экосистем была взята из работы [23], названия типов экосистем приведены в таблице 1). Для удобства интерпретации все это приведено на фоне оттененного рельефа. Акватории, в том числе подо льдом, или территории вне парка везде показаны черным цветом

Таблица 1. Биофизические группы и типы экосистем в национальном парке «Иввавик» (более подробное описание особенностей рельефа, растительности и грунтовых условий можно найти в брошюре [17])

Среднесуточная температура воздуха в ячейке сетки оценивалась на основе ежедневных наблюдений на климатической станции по формуле:

где T_d,g(Y,M,D), T_d,s(Y,M,D) – суточная температура воздуха (максимальная или минимальная температура за сутки) в день D месяца M года Y для ячейки сетки g и климатической станции s соответственно; ΔT_m,gs(Y,M) – разница между среднемесячными температурами воздуха в ячейке сетки g и на климатической станции s для месяца M года Y.

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

Для прогнозирования климатических условий в будущем (2012–2100 гг.) авторы выбрали два сценария, которые были сгенерированы с помощью двух глобальных климатических моделей (Global Climate Models – GCM), таких как CCCma (CGCM3.1) и MPI ECHAM5. Эти два сценария будем называть CCCma и ECHAM соответственно. Они были выбраны потому, что две использованные для их создания модели лучше подходили для моделирования климата XX века в Северной Америке [24] и два полученных на их основе сценария представили общий диапазон изменений температуры воздуха в соответствии со сценариями промежуточных уровней выбросов (A1B) для этого региона. Среднемесячные данные для CCCma и ECHAM были взяты в апреле 2011 года с сайта Всемирного центра климатических данных (http://mud.dkrz.de/ wdcc-for-climate/).

Сначала помесячно рассчитали значения относительной разницы между величинами для спрогнозированного будущего климата и усредненными величинами за прошлый 30-летний период (1971–2000 гг.), полученными с помощью одной и той же модели. Затем эти относительные изменения применили для помесячной оценки будущего климата для всех ячеек сетки НПИ, поскольку парк в основном охватывается одной ячейкой сетки климатической модели. В результате между 1971–2000 и 2090-ми годами прогноз показал увеличение среднегодовой температуры воздуха на 4,2 и 8,8 °C, а общего годового количества осадков – на 26 и 37% по сценариям CCCma и ECHAM соответственно (рис. 3).

Рис. 3. Изменения на климатических станциях «Комакук-Бич» и «Олд-Кроу»: а, б – среднегодовой температуры воздуха: в, г – общего годового количества осадков. Наблюдения на этих двух станциях использовались для представления временнЫх климатических закономерностей в прибрежных и внутренних районах соответственно

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

Парциальное давление водяного пара в воздухе оценивалось на основе минимальной температуры воздуха [13]. Общая суточная инсоляция без учета топографических эффектов (входные данные для учета влияния топографических особенностей на инсоляцию описаны в приложении к статье) оценивалась на основе широты, дня года, суточного диапазона температур и давления пара [13]. Эти параметры определялись на основе наблюдений на климатической станции «Инувик» (68,32° с.ш., 133,52° з.д.).

Прибрежная территория парка, граничащая с Северным Ледовитым океаном, характеризуется морским климатом, а внутренняя территория – континентальным [15]. Авторы настоящей статьи очертили южную границу прибрежного региона по контурной линии 300 м над уровнем моря. Чтобы представить временнЫе климатические закономерности в прибрежных и внутренних регионах, использовались результаты ежедневных наблюдений соответственно на станциях «Комакук-Бич» (за 1958–2011 гг.) и «Олд-Кроу» (1968–2011 гг.), поскольку там проводились самые продолжительные непрерывные наблюдения в пределах или вблизи парка. Пробелы в данных заполнялись на основе измерений на соседних метеорологических станциях.

Классификация усредненных климатических и топографических воздействий на инсоляцию

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

Состояние многолетней мерзлоты в основном зависит от сезонных и годовых климатических условий и не очень чувствительно к ежедневным колебаниям климатических параметров. Например, толщина деятельного слоя грунта в основном определяется общим годовым количеством градусо-дней, когда среднесуточная температура воздуха больше нуля градусов по Цельсию (Total Degree-Days at T>0 °C – TDD_T>0), в соответствии с уравнением Стефана [25]. Поэтому для группирования (кластеризации) средних климатических показателей авторы использовали TDD_T>0 и общее годовое количество осадков. Значения TDD_T>0 с разрешением 30 м и общие годовые количества осадков были повторно дискретизированы до разрешения 10 м с использованием такого метода кластерного анализа, как метод ближайшего соседа.

Топографическими атрибутами, влияющими на инсоляцию, являются крутизна склона, его экспозиция (ориентация) и зона видимости. Крутизна и экспозиция были рассчитаны с использованием цифровой модели рельефа (ЦМР) с разрешением 10 м, которая была повторно дискретизирована на основе данных ЦМР с разрешением 30 м с использованием метода билинейной интерполяции. Авторы сделали эту повторную дискретизацию, чтобы достичь соответствия разрешению карты экосистем. Данные ЦМР были взяты с сайта Центра топографических данных Министерства природных ресурсов Канады.

Зона видимости для ячейки сетки представляет собой угловое распределение видимости неба в зависимости от препятствий, аналогичное видам на полусферических фотографиях, снятых с помощью объектива типа «рыбий глаз», направленного в небо из центра ячейки сетки [26]. Поскольку большинство долин в парке имеют ширину менее 3 км (между вершинами), авторы рассчитали зоны видимости для каждой ячейки сетки размером 30 м на 30 м, используя окно размером 201 на 201 ячейку сетки ЦМР (или 3 км от ячейки сетки до сторон окна). Углы видимости были рассчитаны для каждого азимута через 22,5° (всего для 16 азимутов) для каждой ячейки сетки. Затем углы видимости были интерполированы до пространственного разрешения 10 м.

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

Авторы выполнили кластеризацию средних климатических и топографических показателей на основе TDD_T>0, общего годового количества осадков, среднегодовой инсоляции, крутизны склона и его экспозиции с использованием метода неконтролируемой классификации Isodata с помощью программного обеспечения PCI Isodata от канадской компании PCI Geomatics Enterprises Inc. Поскольку эта программа допускает максимум 255 кластеров, авторы сначала получили 9 больших кластеров для внутреннего региона, а затем каждый из них дополнительно разбили на 254 кластера. Данные для прибрежного региона были напрямую разбиты на 254 кластера, причем ошибки кластеризации были меньше, чем для внутреннего региона. На рисунке 4 показаны примеры климатических параметров, инсоляции и топографических условий для кластеров прибрежного региона. В таблице 2 перечислены средние ошибки кластеризации. Влияние ошибок кластеризации на толщину смоделированного деятельного слоя в целом было менее 0,6% (на основе уравнения Стефана).

Рис. 4. Условия в кластерах прибрежного региона: а – общая годовая инсоляция; б – уклон и экспозиция кластеров, соответствующих той или иной общей годовой инсоляции (если экспозиция больше 180°, то она показана на графике как 360° минус экспозиция); в – общее годовое количество градусо-дней с температурой воздуха выше 0 °C и общее годовое количество осадков; г – зоны видимости для некоторых кластеров с указанной в легенде общей годовой инсоляцией (для склона, обращенного на север, экспозиция равна нулю, а увеличивается она в направлении против часовой стрелки)

Таблица 2. Количество кластеров и средние ошибки кластеризации для прибрежного и внутреннего регионов НПИ

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

Аналогичным образом среднемесячная температура воздуха и общее месячное количество осадков для кластера определялись с использованием широты той ячейки сетки в его пределах, в которой разница между TDD_T>0 для этой ячейки и средним значением по кластеру была наименьшей.

Типы экосистем и их распределение

На рисунке 2, г показано распределение типов экосистем (см. таблицу 1) на территории НПИ. Эту карту составили авторы статьи [23], используя снимки со спутников SPOT с разрешением 10 м, сделанные в период с 13 июля по 22 августа 2006 года. Она была разработана с применением нового метода прогнозного картирования экосистем на основе снимков – с объединением карты растительного покрова по данным дистанционного зондирования и прогнозных атрибутов рельефа по ЦМР  [23]. Классификатор дерева решений был обучен на основе карты экосистем бассейна реки Ферт, разработанной с использованием аэрофотоснимков и детальных полевых наблюдений на 367 участках авторами работы [17]. Общая точность классификации составила 85% [23].

Грунтовые условия и гидрогеологические параметры

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

Грунтовые условия, определенные для разных типов экосистем, приведены в таблице 3.

Таблица 3. Грунтовые условия, определенные для разных типов экосистем

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

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

Содержание органического вещества в верхней части минерального грунта оценивалось в соответствии с типичными результатами полевых измерений. Изменения в содержании органического вещества в грунте с глубиной оценивались на основе полевых данных и общих закономерностей, описанных в работе [27].

Грунты большинства территорий парка, особенно в более глубоких слоях, содержат некоторое количество крупнообломочного материала от гравия до валунов. Его доля оценивалась авторами на основе результатов полевых наблюдений, биофизических и ландшафтных особенностей типов экосистем, описанных в работе [17]. Этот материал уменьшает объем грунта, способный удерживать воду и обмениваться ей. Поэтому мощность подстилающего грунта оценивалась с допущением того, что доля крупнообломочного материала увеличивается с глубиной на 2% каждые 10 см до тех пор, пока не достигнет 100%.

Параметры бокового притока и оттока оценивались на основе класса дренируемости и класса увлажнения (водного режима) грунта для типов экосистем в соответствии с работой [17]. Авторы настоящей статьи приняли, что наименьшая глубина уровня грунтовых вод (УГВ) для бокового оттока, составляющая от 0 до 0,4 м, пропорциональна номеру класса увлажнения грунта (в диапазоне от 0 до 8, то есть от очень сухого до водонасыщающего соответственно), а скорость бокового оттока экспоненциально зависит от номера класса дренируемости (в диапазоне о  0 до 6, то есть от очень слабой до очень быстрой соответственно). Авторы также приняли, что нет поверхностного бокового притока и вода будет постепенно утекать, когда УГВ будет выше поверхности земли (если последнее случится, то УГВ будет становиться на 5% глубже каждый день).

Параметр переноса снега ветром рассчитывался с учетом влияния экспозиции поверхности, форм и плотности растительного покрова.

Экспозицию поверхности в ячейке сетки определяли на основе типичных особенностей рельефа для типов экосистем, описанных в работе [17].

Растительность подразделили на шесть форм: хвойные деревья; лиственные деревья; кустарники от средней и большой высоты (0,5 м); низкие кустарники (ниже 0,5 м); осоки или злаки; низкие осоки или злаки. Влияние плотности растительности оценивали на основе индексов листовой поверхности (Leaf Area Indices, LAI).

Поскольку геотермические наблюдения в арктических регионах были редки, было принято, что все ячейки сетки территории парка имеют одинаковую плотность геотермального теплового потока (0,08 Вт/м²) и одинаковую теплопроводность коренных пород (2,28 Вт /(м·°С), оцененные на основе результатов полевых наблюдений [30–32].

Индексы листовой поверхности

Значения индекса листовой поверхности (Leaf Area Index – LAI) были доведены до пространственного разрешения 10 м с использованием оценок биомассы листвы по снимкам со спутников Landsat с разрешением 30 м с помощью следующей формулы:

где где L_10m,i - расчетное значение LAI для 10-метрового i-го пикселя (см² листовой поверхности на см² поверхности грунта); E_i, E_j – типы экосистем для 10-метровых i-го и j-го пикселей соответственно (j варьирует от 1 до 9 для девяти 10-метровых пикселей в 30-метровом пикселе снимка со спутника Landsat); B_30m - биомасса листвы в пределах 30-метрового пикселя (г/см²) которая рассчитывается с использованием уравнения регрессии простого соотношения вегетационного индекса (полоса-4/полоса-3) по снимкам, полученным со спутника Landsat 30 августа 2001 года [33]; B₀(E_j) - средняя биомасса листвы (г/см²) для экосистемы типа Ej, которая была рассчитана как среднее значение для всех 30-метровых пикселей (в пределах территории парка), где E_j – доминирующий тип экосистемы; S(E_i) - удельная листовая поверхность (см²/г) для типа экосистемы E_i (122 см²/г для осок, 154 см²/г для листопадных кустарников, 74 см²/г для вечнозеленых кустарников, 50  см²/г для ели черной) [34].

На снимках, сделанных со спутников Landsat, некоторые крутые участки были скрыты тенью. Поэтому авторы настоящей статьи оценивали LAI на таких участках как среднее значение LAI для экосистемы типа B₀(E_i)S(E_i). Распределение LAI на территории парка показано на рисунке 2, в. При моделировании значения LAI были разделены на 18 классов, которые экспоненциально увеличивались с 0 до 4,0.

Процедуры расчетов

Полученная модель запускалась для всех комбинаций классов входных данных (типов экосистем, климатотопографических кластеров, классов LAI) для прибрежного и внутреннего регионов. Одна комбинация типов/классов/кластеров входных данных обычно представляла условия для более чем одной ячейки сетки. За исключением поверхностей воды и наледей, для прибрежного и внутреннего регионов парка есть соответственно 49,441 тыс. и 258,633 тыс. различных комбинаций. Общее количество комбинаций классов входных данных составило 0,3% от общего количества ячеек сетки сухопутной части парка с разрешением 10 м (в парке насчитывается 107 млн ячеек сетки с разрешением 10 м, не считая водных объектов и наледей). В результате при таком пространственном моделировании с высоким разрешением время вычислений значительно сократилось.

Посуточные климатические данные, полученные на станциях «Олд-Кроу» и «Комакук-Бич», относятся к 1968 и 1958 годам соответственно. Для пространственной согласованности авторы моделировали многолетнюю мерзлоту в парке с 1968 года. Для инициализации (задания начальных условий) модели были выбраны посуточные климатические данные наблюдений на станциях «Олд-Кроу» и «Комакук-Бич» в 1972 и 1964 годах соответственно, когда климатические условия были близки к таковым при экстраполяции на 1967 год. Модель прогоняли итеративно, используя климатические данные для указанных лет, до тех пор, пока смоделированная температура грунта не стала стабильной. Затем ее прогоняли по данным за 1968–2011 годы на основе ежедневных наблюдений. После этого модель запускалась для двух рассмотренных ранее климатических сценариев до 2100 года.

ПОЛУЧЕННЫЕ РЕЗУЛЬТЫ И ИХ АНАЛИЗ

Пространственное распределение толщины деятельного слоя

На рисунках 5, а, б показано смоделированное пространственное распределение толщины слоя сезонного оттаивания/промерзания в 2000-х годах (2000–2009 гг.). В основном она уменьшалась на территории парка с юга на север. Деятельный слой был тоньше в большинстве долин, на прибрежных территориях. Более толстым он был выше в горах и вдоль дренирующих русел в некоторых долинах. Полученные результаты также показали значительные различия для территорий, занятых разными типами экосистем (рис. 6, в). Горные склоны и скальные участки, поросшие лишайниками, имели мощные деятельные слои, поскольку на их поверхности не было органического слоя и в грунте обычно содержалось большое количество крупнообломочного материала. В действующих поймах (с типами экосистем 23–26) и дренирующих руслах (с типами экосистем 11–13) деятельный слой также был более мощным из-за очень влажных условий и высокого содержания крупнообломочного материала. Склоны педиментов, то есть слабонаклонные предгорные/подгорные равнины (с типами экосистем 14–17), и заболоченные территории, поросшие злаками (экосистемы типа 21), имели тонкий активный слой из-за мощного слоя торфа и очень влажных условий.

Рис. 5. Результаты моделирования: а – средняя толщина деятельного слоя в 2000-х годах (2000–2009 гг.); б – увеличенное изображение части модели, ограниченной черным квадратом на рисунке «а»; в, г – относительные изменения толщины деятельного слоя с 2000-х по 2090-е годы в соответствии с климатическими сценариями CCCma и ECHAM соответственно. Для удобства интерпретации все это приведено на фоне оттененного рельефа. Акватории, в том числе подо льдом, или территории вне парка везде показаны черным цветом

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

Рис. 6. Средняя толщина деятельного слоя грунта в 1970-х годах в климатотопографических кластерах парка с разными: а – средними летними (с июня по сентябрь) температурами воздуха; б – общегодовой инсоляцией; в – типами экосистем в соответствии с таблицей 1 (показаны также средние абсолютные отклонения в пределах классов); г – индексами листовой поверхности (показаны также средние абсолютные отклонения в пределах классов)

Чтобы сравнить относительную важность входных переменных для смоделированного пространственного распределения толщины слоя сезонного протаивания/промерзазния в парке, были рассчитаны средние абсолютные отклонения внутри каждого типа/кластера/класса (далее – класса) и между разными классами на основе входных переменных для типов экосистем, климатотопографических кластеров и классов индексов листовой поверхности (таблица 4). Среднее отклонение между классами является мерой различий между ними, в то время как среднее отклонение внутри каждого класса является мерой изменчивости внутри него.

Таблица 4. Средние абсолютные межклассовые и внутриклассовые отклонения для толщины активного слоя на территории парка в 1970-х годах и соотношения между ними

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

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

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

Изменения толщины деятельного слоя и распространения многолетней мерзлоты

Результаты моделирования показали, что средняя толщина слоя сезонного протаивания/промерзания в НПИ с 1970-х по 2000-е годы (2000–2009 гг.) выросла на 0,1 м (13%). Это увеличение было более значительным на возвышенностях в южной части парка. С 2000-х годов до конца XXI века смоделированная средняя толщина активного слоя выросла на 0,3 м (34%) и на 0,9 м (99%) в соответствии со сценариями изменений климата CCCma и ECHAM соответственно (рис. 7, а), то есть разница между этими двумя сценариями оказалась существенной. Согласно CCCma относительные изменения толщины деятельного слоя в основном составляли от 20 до 40% по всему парку. Однако по сценарию ECHAM они были больше и в более широком диапазоне (в основном от 70 до 150%) (см. рис. 5, г).

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

Рис. 7. Изменения средней толщины слоя сезонного протаивания/промерзания (а) и изменения долей площади парка (%), находящихся на разных стадиях деградации многолетней мерзлоты согласно сценариям изменений климата CCCma и ECHAM (б, в соответственно). Стадии деградации приведены на основе данных из работы [14]

С прогрессирующим увеличением глубины летнего протаивания глубокий слой, оттаявший летом, может не замерзнуть следующей зимой и образовать талик (круглогодично незамерзаюший слой) между многолетней мерзлотой и слоем сезонного протаивания. В настоящее время многолетняя мерзлота подстилает почти все участки суши в парке и лишь очень небольшая доля площади содержит талики (в 2000-х годах общие площади участков с таликами и вообще без многолетней мерзлоты составляли в парке всего 0,0007 и 0,004% соответственно). По модельному прогнозу, с середины XXI века талики будут формироваться на большем количестве учатков, а в некоторых южных районах парка многолетняя мерзлота вообще исчезнет. К концу XXI века по сценариям CCCma и ECHAM талики образуются на 0,1 и 13,7% площади парка соответственно, а многолетняя мерзлота исчезнет на 2,1 и 5,2% площади соответственно. Развитие таликов будет происходить в основном в южных районах парка (см. рис. 5 в, г). Талики могли бы формироваться в большинстве типов экосистем, за исключением склонов педиментов и заболоченных территорий, поросших злаками. Многолетняя мерзлота могла бы исчезнуть в лесах и в некоторых южных долинах – в основном в поймах и на неактивных аллювиальных террасах.

Стадии деградации многолетней мерзлоты

В работе [14] процесс деградации многолетней мерзлоты, вызванный потеплением климата, был подразделен на пять стадий: постепенного оттаивания, повышенного оттаивания, часто встречающихся таликов, изотермической ММ и отсутствия ММ. Авторы настоящей статьи определили эти стадии для каждой ячейки сетки территории НПИ, основываясь на результатах моделирования.

С 1967 по 2011 год 8,7% площади распространения многолетней мерзлоты в парке перешли из стадии постепенного оттаивания в стадию повышенного оттаивания с очень небольшой долей (<0,01%) площади, находящейся на стадии множественных таликов и отсутствия ММ (см. рис. 7, б). Согласно модельному прогнозу к концу XXI века все бОльшая доля территории будет находиться на более высоких стадиях деградации многолетней мерзлоты, хотя между результатами для двух рассматриваемых климатических сценариев имеются значительные различия.

Для сценария CCCma доля площади с постепенно оттаивающей ММ сократится к концу XXI века до 84,8%, в то время как доли территории парка на стадии повышенного оттаивания и без ММ увеличатся до 12,8 и 2,1% соответственно (см. рис. 7, б). В соответствии с тем же сценарием зоны с многолетней мерзлотой на стадии повышенного оттаивания и зоны без ММ в основном будут располагаться в некоторых долинах и на склонах, обращенных на юг.

По сценарию ECHAM прогнозируется более значительная деградация ММ: к концу XXI века только 19,9% площади парка будут содержать постепенно оттаивающую многолетнюю мерзлоту. Доля площади парка с ММ на стадии повышенного оттаивания увеличится до 53,0%. Территории на других стадиях деградации также расширятся (см. рис. 7, в). Северная часть парка в основном будет находиться на стадии повышенного оттаивания, в то время южная часть – главным образом на стадиях множественных таликов и изотермической многолетней мерзлоты, а некоторые зоны будут вообще свободны от ММ.

Сравнение с результатами измерений и других исследований

На рисунке 8 показана корреляционная связь между смоделированной и измеренной глубиной летнего оттаивания на 100 участках. Они были сопоставимы для большинства участков, но для некоторых были большие различия, особенно для тех семи, данные по которым обведены на рисунке 8 круговой штриховой линией (где модельные данные значительно завышены по сравнению с измеренными). Эти семь участков находятся на горных склонах, где на поверхности земли и в грунте имеется много каменного материала. Наблюдаемые небольшие глубины летнего оттаивания на этих участках, вероятно, обусловлены локальными вариациями в грунтовых условиях или тем, что скальные породы ошибочно принимаются за мерзлый грунт. На некоторых горных склонах было обнаружено летнее протаивание глубиной более 1 м, однако при большинстве измерений мерзлый слой не был достигнут (поэтому такие данные не были отражены на рисунке 8). Модельные тесты тоже показали, что деятельный слой на горных склонах был толстым – из-за отсутствия органического слоя, высокого содержания крупнобломочного материала и редкой растительности [35]. Коэффициент корреляции R между смоделированной и наблюдаемой глубиной летнего протаивания оказался низким (R=0,56 при количестве участков измерений N=93) даже без учета данных для упомянутых выше семи участков с наибольшими различиями. Низкая корреляция в основном обусловлена вариациями грунтовых условий в пределах экосистемы того или иного типа. Поэтому для картирования многолетней мерзлоты с высоким разрешением необходимы более детальные в пространственном отношении карты грунтовых условий.

Рис. 8. Корреляционная связь между смоделированной и измеренной глубиной летнего протаивания. Модель значительно завысила данные по семи участкам (обведены штриховой окружностью), обладающим поросшей лишайниками скальной поверхностью и находящимся на горных склонах, – вероятно, из-за пространственной неоднородности грунтовых условий или ошибочного принятия скальных пород за мерзлый грунт. Сплошная прямая линия соответствует уравнению линейной регрессии, приведенному в левом верхнем углу. Штриховая прямая – это линия несуществующего соответствия 1:1 (для справки)

На основе полевых наблюдений автор работы [36] создал полигональную карту толщины деятельного слоя для территории НПИ. По этой карте медианная толщина слоя сезонного протаивания составляет 0,4 м (0,3–0,5 м) в большинстве северных прибрежных районов парка и 0,75 м (0,7–1,0 м) – в северных дельтах. Медианный деятельный слой вдоль крупных долин составляет 0,7 м (0,5–0,9 м), а во внутренних горных районах – 1,35 м (0,5–2,0 м). Эта пространственная картина и диапазон толщины деятельного слоя аналогичны результатам моделирования, выполненного авторами настоящей статьи, хотя на карте из работы [36] на территории парка было представлено только 15 полигонов. Смоделированная многолетняя мерзлота в этом регионе является сплошной, что согласуется с картой ММ Канады [16]. Смоделированная толщина ММ составила 150–300 м, что соответствует результатам наблюдений в этом регионе [37].

Авторы работы [38] сообщили об изменениях толщины деятельного слоя и температуры грунта на острове Хершел недалеко от прибрежного региона. Откалибровав использованную ими модель с помощью детальных результатов полевых измерений, они подсчитали, что с начала XX века (1899–1905 гг.) по 2006 год среднегодовая температура грунта на глубине 1 и 20 м увеличилась на 2,6 и 1,9 °C соответственно, а изменение температуры грунта (хотя и менее чем на 0,1 °C) достигло примерно 120 м глубины.

Авторы настоящей статьи проверили смоделированные профили температуры грунта для схожих условий почвенно-растительного покрова (для осоковых кочкарников) в прибрежном регионе НПИ. С 1967 года (при условии равновесия) по 2006 год смоделированная температура грунта на глубине 1 и 20 м увеличилась на 3,4 и 2,5 °C соответственно. Эти значения получились больше, чем в работе [38], хотя рассматривалось более короткое время. Вероятно, это было связано с тем, что в климатических данных, выбранных для задания начальных условий модели, было меньше атмосферных осадков, а тренды к увеличению количества осадков и толщины снежного покрова за последние 45 лет не могли соответствовать тенденциям за весь XX век, поскольку в 1958–1967 годах было зафиксировано бОльшее количество осадков (см. рис. 3, в). Результаты моделирования, выполненного авторами настоящей статьи, показали, что температура грунта на высотной отметке 120 м увеличилась на 0,04 °C, а это аналогично результату из статьи [38].

В работе [38] также представлены данные измерений глубины летнего протаивания грунта на острове Хершел в период с 2003 по 2007 год. Ее величина и межгодовые колебания оказались сопоставимыми с толщиной деятельного слоя, смоделированной авторами настоящей статьи, за исключением 2007 года, когда наблюдения показали значительное увеличение (эти результаты показаны точками в виде черных квадратиков на рисунке 9).

Рис. 9. Корреляционная связь между смоделированной толщиной деятельного слоя в прибрежной зоне НПИ и глубиной летнего протаивания грунта в районе станции «Иллисарвик» на острове Ричардс, измеренной в период 1983–2012 гг. (белые кружки) и на острове Хершел в 2003–2007 гг. (черные квадратики). Наблюдения в районе станции «Иллисарвик» на острове Ричардс в 1983–2008 годах проводились авторами работы [39], а неопубликованные данные измерений в 2009–2012 годах были любезно предоставлены Кристофером Берном (Christopher Burn). Наблюдения на острове Хершел выполнялись авторами статьи [39]. Смоделированная толщина активного слоя в 2012 году – это среднее модельное значение, рассчитанное на основе двух климатических сценариев для этого года

В работе [39] приведены данные долговременных измерений глубины летнего протаивания с 1983 года в районе станции «Иллисарвик» на острове Ричардс (69,48° с.ш., 134,59° з.д.) (см. красную точку на рисунке 1, б). Эти наблюдения, вероятно, были самыми продолжительными из непрерывных измерений глубины летнего протаивания в этом регионе. И с ними тесно коррелирует (R=0,71, N=30) толщина деятельного слоя по данным моделирования, выполненного авторами настоящей статьи, хотя абсолютные значения несколько отличаются из-за различий в климатических и грунтовых условиях между местами наблюдений и прибрежным регионом НПИ (см. рис. 9). Смоделированный тренд увеличения мощности активного слоя составил 2,9 см за десятилетие с 1983 по 2012 год, что было похоже на тренд по результатам наблюдений, который показал увеличение на 2,3 см за десятилетие.

ОБСУЖДЕНИЕ ПОЛУЧЕННЫХ РЕЗУЛЬАТОВ

Авторами настоящей статьи была составлена карта распространения многолетней мерзлоты и изменений климата в арктическом регионе со сложным рельефом с пространственным разрешением 10 м на основе модели, основанной на процессах.

Стоит отметить несколько особенностей представленного исследования по сравнению с предыдущими. Здесь использовалось более высокое пространственное разрешение, чем в ранее опубликованных работах по моделированию и картографированию, например с разрешением 30 м [12, 13], 2 км [11], в половину градуса по широте и долготе или даже с более грубым разрешением [7, 10]. Данные дистанционного зондирования Земли со спутников были использованы более эффективно. Результаты влияния топографии на инсоляцию были рассчитаны более точно. Размер каждого участка полевых наблюдений был сопоставим с размером ячейки сетки (поскольку большие однородные участки в регионах со сложным рельефом обычно найти трудно), поэтому данные наблюдений можно было непосредственно использовать для калибровки и валидации модели. Полученные результаты выявили больше пространственных деталей и, следовательно, больше подошли для планирования освоения землель, экологического мониторинга и оценок.

Для проверки влияния пространственного разрешения на моделируемые условия многолетней мерзлоты была выбрана территория размером 20 км на 20 км (сопоставимая по площади с размером ячеек сетки в полуградуса по широте и долготе, который в рассматриваемом регионе составляет около 20 км на 55 км) и было выполнено моделирование с разрешением 100 м, 1 км, 5 км и 20 км на основе доминирующего типа почвенно-растительного покрова в каждой ячейке сетки. Более грубое разрешение не только выявило меньше пространственных деталей (рис. 10), но и дало результаты со смещением характеристик ММ в сторону показателей, свойственных доминирующему типу почвенно-растительного покрова На указанной тестовой территории доминирующим типом экосистемы является каменистая поверхность, покрытая лишайниками, где слой сезонного протаивания более мощный, чем в большинстве других типов. С увеличением размера ячеек сетки смоделированная средняя толщина деятельного слоя на этой территории увеличилась на 22%, а пространственные вариации между ячейками сетки уменьшились с 0,4 до 0 м при ухудшении разрешения с 10 м до 20 км (рис. 11). ВременнЫе закономерности смоделированной средней мощности активного слоя для данной территории при различных пространственных разрешениях были схожи (коэффициенты корреляции R варьировали в пределах 0,97–1,00 при количестве данных N=134).

Рис. 10. Смоделированное пространственное распределение толщины деятельного слоя (ее средних значений за 2000-е годы) при различных пространственных разрешениях для тестовой территории размером 20 км на 20 км (а–ж). Расположение тестовой территории в парке показано красным квадратом (з). Чтобы показать разницу между изображениями «а» и «б», небольшая область, обведенная черные квадратами на рисунках «а» и «б», была увеличена и показана на рисунках «г» и «д».

Рис. 11. Влияние пространственного разрешения на толщину смоделированного деятельного слоя на тестовой территории размером 20 км на 20 км (см. рис. 10, з). Средние значения мощности активного слоя были рассчитаны для всех ячеек сетки в пределах этой территории за период 1967–2100 гг. Пространственные стандартные отклонения были вычислены на основе средних значений для этих ячеек за все эти годы (1967–2100). данные для будущих лет (2012–2100) были смоделированы на основе сценария изменений климата ECHAM (при использовании сценария CCCma результаты получились очень похожими)

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

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

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

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

Еще одним важным источником неопределенности являются представления о будущем климате. Результаты выполненного исследования показали существенные различия в прогнозируемых условиях ММ для двух климатических сценариев. При более «теплом» сценарии (ECHAM) прогнозируемые изменения в толщине деятельного слоя характеризовались не только в целом бОльшим увеличением, но и более широким диапазоном флуктуаций, а территории с таликами или даже с исчезновением многолетней мерзлоты имели гораздо бОльшую общую площадь, чем при более «холодном» сценарии (CCCma). На точность модельных результатов также влияют сильная территориальная разбросанность и относительно короткие периоды сбора климатических данных в арктических регионах. Климатические параметры первого года, выбранного для задания начальных условий модели, могут повлиять на смоделированную толщину многолетней мерзлоты и долговременные тренды температурных режимов грунтов.

Хотя авторы учитывали перенос снега ветром и латеральные водные потоки, созданная модель NEST все же является одномерной при допущениях, что каждая ячейка сетки однородна, а боковым тепловым потоком можно пренебречь. Поэтому полученные результаты не отражают условия для территорий с сильными боковыми тепловыми потоками, таких как островерхие горные пики, гребни хребтов или районы, очень близкие к водным объектам.

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

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

ЗАКЛЮЧЕНИЕ

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

Результаты моделирования с высоким разрешением показали как крупномасштабные изменения по всему канадскому национальному парку «Иввавик», так и изменения в пределах того или иного ландшафта, обусловленные топографией, экосистемами и характеристиками растительного покрова.

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

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

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

ПРИЛОЖЕНИЕ

Расчет влияния топографических особенностей на инсоляцию

Рельеф влияет на инсоляцию поверхности земли в зависимости от крутизны склона, локальных топографических особенностей и затеняющего воздействия окружающих гор и холмов. Эффекты затенения рассчитывались на основе зон видимости [26]. Зона видимости для ячейки сетки – это угловое распределение видимости неба, частично загороженного окружающими горами и холмами, аналогичное изображениям на полусферических фотографиях, снятых с помощью объектива типа «рыбий глаз», направленного в небо из центра ячейки сетки. Максимальный угол закрывающего небо препятствия в направлении азимута может быть определен на основе расстояния и высотных отметок ячеек в этом направлении.

Зона видимости влияет на эффекты как прямой, так и рассеянной солнечной радиации. Интенсивность рассеянной радиации, полученной ячейкой сетки (R_diff), может быть рассчитана по следующей формуле:

где R_diff0 - рассеянная солнечная радиация, полученная поверхностью без загораживающих препятствий вокруг (Вт/м₂); F_diff,topo – относительная рассеянная радиация, полученная ячейкой сетки, по сравнению с ячейкой без окружающих препятствий (отношение).

Если принять, что рассеянная солнечная радиация поступает только с неба и что ее распределение является изотропным по всей полусфере [41], а также если пренебречь эффектами крутизны и экспозиции склона, то F_diff,topo можно оценить так:

где α – азимут (азимутальный угол); β – угол возвышения (относительно горизонтальной поверхности); ν(α) – максимальный угол возвышения, закрываемого окружающими горами и холмами в направлении α.

Прямая солнечная радиация, получаемая ячейкой сетки за раз, может быть рассчитана по формуле:

где α, θ – азимутальный и зенитный углы солнца в данный момент времени, (они могут быть определены на основе широты, дня года и времени суток); G_s, G_a – крутизна и экспозиция поверхности ячейки сетки; R_dir0 - прямая солнечная радиация, полученная поверхностью без загораживающих препятствий вокруг (Вт/м²).

Относительное воздействие особенностей рельефа на получаемую прямую солнечную радиацию (F_dir,topo) можно выразить так:

Общая инсоляция ячейки сетки в любое время суток, представляет собой сумму прямой и рассеянной радиации, попавшей на эту ячейку.

Для пользователей сложно и обычно не нужно предоставлять почасовые климатические данные для долгосрочного моделирования многолетней мерзлоты, хотя они необходимы для расчета энергетического баланса. Поэтому авторы настоящей статьи оценили внутрисуточные вариации как прямой, так и рассеянной солнечной радиации для каждой даты, используя функции косинуса на основе работы [42]:

где θ _n – зенитный угол солнца в полдень (рад) в рассматриваемый день, F_1,θ_n – отношение прямой солнечной радиации в полдень к среднесуточной прямой солнечной радиации, попавшей на горизонтальную поверхность; F_2,θ_n – отношение рассеянной солнечной радиации в полдень к среднесуточной рассеянной солнечной радиации, попавшей на горизонтальную поверхность; R_day0 – среднесуточная инсоляция горизонтальной поверхности; (Вт/м²); F_diff0,day – доля суточной рассеянной солнечной радиации в общей суточной радиации, попавшей на горизонтальную поверхность (общей инсоляции).

R_day0 можно рассчитать по формуле:

где R'_day0 – суммарная суточная инсоляция горизонтальной поверхности (МДж/(м²•сут)), которая относится к входным данными для модели; D – длительность светового дня (с), определяемая на основе широты и дня года.

Авторы работы [42] разработали регрессионные зависимости для определения F_1,θ_n и F_2,θ_n на основе θ_n для той или иной локации. Численное моделирование, выполненное авторами настоящей статьи, показало, что эти зависимости меняются с широтой, особенно когда продолжительность светового дня достигает 24 ч. Поэтому F_1,θ_n и F_2,θ_n были численно определены для каждого дня года на основе исходных уравнений из работы [42]:

где t – время (с).

Суточная доля рассеянной солнечной радиации в общей суточной инсоляции горизонтальной поверхности (F_diff0,day) может быть оценена с помощью логистической функции по работе [43]. Авторы настоящей статьи немного модифицировали параметры, основываясь на данных по суточной инсоляции, полученных на станции «Инувик» (133,52 Вт, 68,32° с.ш.):

где F₀ – отношение суточной инсоляции горизонтальной поверхности, к суточной инсоляции верхней границы атмосферы; R² – коэффициент детерминации (квадрат коэффициента корреляции); N – количество валидных суточных данных станции «Инувик» за период 1959–2005 гг., опубликованных Министерством окружающей среды Канады и Национальным исследовательским советом Канады [44].

Суточную инсоляцию верхней границы атмосферы (R_0,atm, МДж/(м²•сут)) можно рассчитать по следующей формуле [45]:

[an error occurred while processing this directive]

где S_c — солнечная постоянная (Вт/м²).

Дополнительный материал по исходной программе (на языке C++), использованной для реализации этих алгоритмов, доступен онлайн (по адресу: http://www.the-cryosphere.net/7/1121/2013/tc-7-1121-2013-supplement.pdf ).

-

Авторы хотели бы поблагодарить Ин Чжана (Ying Zhang) и Сергея Самсонова (Sergey Samsonov) за критический анализ статьи. Также выражается благодарность Кристоферу Берну (Christopher Burn) за его замечания и полевые данные для валидации результатов. Для улучшения качества статьи также были очень полезны замечания двух анонимных рецензентов. Представленное в статье исследование было поддержано Канадским космическим агентством в рамках проекта ParkSPACE GRIP и Министерством природных ресурсов Канады в рамках Программы дистанционного зондирования (вклада сектора наук о Земле № 20120211). Редактировал статью (исходную статью на английском языке. – Ред.) Т. Чжан.

ИСТОЧНИК ДЛЯ ПЕРЕВОДА

Zhang Yu., Wang X., Fraser R., Olthof I., Chen W., Mclennan D., Ponomarenko S., Wu W.Modelling and mapping climate change impacts on permafrost at high spatial resolution for a region with complex terrain // The Cryosphere. 2013. Vol. 7. P. 1121–1137. DOI:10.5194/tc-7-1121-2013.

СПИСОК ЛИТЕРАТУРЫ, ИСПОЛЬЗОВАННОЙ АВТОРАМИ ПЕРЕВЕДЕННОЙ СТАТЬИ

1. ACIA. Impact of a Warming Arctic: Arctic Climate Impact Assessment. Cambridge University Press, 2005. 1042 p.
2. Vitt D.H., Halsey L.A., Zoltai S.C. The changing landscape of Canada’s western boreal forest: the current dynamics of permafrost // Can. J. Forest Res. 2000. Vol. 30. P. 283–287.
3. Smith S.L., Romanovsky V.E., Lewkowicz A.G., Burn C.R., Allard  M., Clow G.D., Yoshikawa K., Throop J. Thermal state of permafrost in North America: a contribution to the International Polar Year // Permafrost Periglac. 2010. Vol. 21. P. 117–135. DOI:10.1002/ppp.690,
4. Riseborough D., Shiklomanov N., Etzelmuller B., Gruber S., Marchenko S. Recent advances in permafrost modelling // Permafrost Periglac. 2008. Vol. 19. P. 137–156. DOI:10.1002/ppp.615.
5. Shiklomanov N.I., Anisimov O.A., Zhang T., Marchenko S., Nelson F.E., Oelke C. Comparison of model produced active layer fields: results for northern Alaska // J. Geophys. Res. 2007. Vol. 112. Paper F02S10. DOI:10.1029/2006JF000571.
6. Osterkamp T.E. The recent warming of permafrost in Alaska // Global Planet. Change. 2005. Vol. 49. P. 187–202. DOI:10.1016/j.gloplacha.2005.09.001.
7. Zhang Y., Chen W., Riseborough D.W. Disequilibrium response of permafrost thaw to climate warming in Canada over 1850–2100 // Geophys. Res. Lett. 2008. Vol. 35. Paper L02502. DOI:10.1029/2007GL032117.
8. Anisimov O., Reneva S. Permafrost and changing climate: the Russian perspective // Ambio. 2006. Vol. 35. P. 169–175.
9. Marchenko S., Romanovsky V., Tipenko G. Numerical modeling of spatial permafrost dynamics in Alaska // Proceedings of the 9th International Conference on Permafrost (edited by Kane D. L., Hinkel K.M.). US: Institute of Northern Engineering, University of Alaska Fairbanks, 2008. P. 190–204.
10. Zhang Y. Chen W., Riseborough D.W. Temporal and spatial changes of permafrost in Canada since the end of the Little Ice Age // J. Geophy. Res. 2006. Vol. 111. Paper D22103. DOI:10.1029/2006JD007284.
11. Jafarov E.E., Marchenko S.S., Romanovsky V.E. Numerical modeling of permafrost dynamics in Alaska using a high spatial resolution dataset // The Cryosphere. 2012. Vol. 6. P. 613–624. DOI:10.5194/tc- 6-613-2012.
12. Duchesne C., Wright J.F., Ednie M. High-resolution numerical modeling of climate change impacts on permafrost in the vicinities of Inuvik, Norman Wells, and Fort Simpson, NT, Canada // Proceedings of the 9th  International Conference on Permafrost (edited by Kane D.L., Hinkel K.M.). US: Institute of Northern Engineering, University of Alaska Fairbanks, 2008. P. 385–390.
13. Zhang Y., Li J.,Wang X., Chen W., Sladen W., Dyke L., Dredge L., Poitevin J., McLennan D., Stewart H., Kowalchuk S., Wu W., Kershaw G.P., Brook R.K. Modelling and mapping permafrost at high spatial resolution in Wapusk National Park, Hudson Bay Lowlands // Can. J. Earth Sci. 2012. Vol. 49. P. 925–937. DOI:10.1139/E2012-031.
14. Zhang Y. Spatio-temporal features of permafrost thaw projected from long-term high-resolution modeling for a region in the Hudson Bay Lowlands in Canada // J. Geophys. Res.-Earth. 2013. Vol. 118. P. 1–11. DOI:10.1002/jgrf.20045.
15. Canadian Parks Service: Northern Yukon National Park Resource Description and Analysis: RM REPORT 93-01/INP. Natural Resources Conservation Section, Canadian Parks Service, Prairie and Northern Region. Winnipeg, 1993.
16. Heginbottom J.A., Dubreuil M.A., Harker P.A. Canada Permafrost // National Atlas of Canada (5th edn.). Ottawa, Canada: Natural Resources Canada, 1995.
17. Mackenzie W., MacHutchon A.G. Habitat Classification for the Firth River Valley, Ivvavik National Park, Yukon, Western Arctic District. Inuvik, Parks Canada, 1996. 81 p.
18. Zhang Y., Chen W., Cihlar J. A process-based model for quantifying the impact of climate change on permafrost thermal regimes // J. Geophys. Res. 2003. Vol. 108. Paper 4695. DOI:10.1029/2002JD003354.
19. Zhang Y., Chen W., Smith S.L., Riseborough D.W., Cihlar J. Soil temperature in Canada during the twentieth century: complex responses to atmospheric climate change // J. Geophys. Res. 2005 . Vol. 110. Paper D03112. DOI:10.1029/2004JD004910.
20. Zhang Y., Li C., Trettin C.C., Li H., Sun G. An integrated model of soil, hydrology and vegetation for carbon dynamics in wetland ecosystems // Global Biogeochem. Cy. 2002. Vol. 16. Paper 1061. DOI:10.1029/2001GB001838.
21. Wang T., Hamann A., Spittlehouse D., Aitken S.N. Development of scale-free climate data for western Canada for use in resource management // Int. J. Climatol. 2006 . Vol. 26. P. 383–397. DOI:10.1002/joc.1247.
22. McKenney D.W., Hutchinson M.F., Papadopol P., Lawrence K., Pedlar J., Campbell K., Milewska E.M., Hopkinson R.F., Price D., Owen T. Customized spatial climate models for North America // B. Am. Meteorol. Soc. 2011. Vol. 92. P. 1611–1622. DOI:10.1175/2011BAMS3132.1.
23. Fraser R., McLennan D., Ponomarenko S., Olthof I. Image-based predictive ecosystem mapping in Canadian Arctic Parks // Int. J. Appl. Earth Obs. 2012. Vol. 14. P. 129–138. DOI:10.1016/j.jag.2011.08.013.
24. Walsh J.E., Chapman W.L., Romanovsky V., Christensen J.H., Stendel M. Global climate model performance over Alaska and Greenland // J. Climate. 2008. Vol. 21. P. 6156–6174. DOI:10.1175/2008JCLI2163.1.
25. Lunardini V. Heat Transfer in Cold Climates. New York: Van Nostrand Reinhold, 1981. 731 p.
26. Fu P., Rich P.M. The Solar Analyst 1.0 Manual. USA: Helios Environmental Modeling Institute (HEMI), 2000.
27. Hossain M.F., Zhang Y., Chen W., Wang J., Pavlic G. Soil organic carbon content in northern Canada: a database of field measurements and its analysis // Can. J. Soil Sci. 2007. Vol. 87 . P. 259–268.
28. Agriculture Canada Expert Committee on Soil Survey: The Canadian System of Soil Classification, 2nd edn. Ottawa, Canada: Agriculture Canada, 1987. Publication 1646. 164 p.
29. Walmsley M., Utzig G., Vold T., Moon D., Van Barneveld J. (eds.). Describing Ecosystems in the Field. RAB Technical Paper 2. Victoria, BC Canada: Resource Analysis Branch, British Columbia Ministry of Environment and British Columbia Ministry of Forests, 1980. Rep. № 7.
30. Pollack H.N., Hurter S.J., Johnson J.R. Heat flow from the earth’s interior: analysis of the global data set // Rev. Geophys. 1993. Vol. 31. P. 267–280.
31. Jessop A.M., Lewis T.J., Judge A.S., Taylor A.E., Drury M.J. Terrestrial heat flow in Canada // Tectonophysics. 1984. Vol. 103. P. 239–261.
32. Majorowicz J.A., Skinner W.R., Safanda J. Large ground warming in the Canadian Arctic inferred from inversions of temperature logs // Earth Planet. Sc. Lett. 2004. Vol. 221. P. 15–25.
33. Chen W., Zorn P., Chen Z., Latifovic R., Zhang Y., Li J., Quirouette J., Olthof I., Fraser R., Mclennan D., Poitevin J., Stewart H.M., Sharma R. Propagation of errors associated with scaling foliage biomass from field measurements to remote sensing data over a Canada’s northern national park // Remote Sens. Environ. 2013 . Vol. 130. P. 205–218. DOI:10.1016/j.rse.2012.11.012.
34. Kattge K., Diaz S. Lavorel S., Prentice I.C., Leadley P., et al. TRY-a global database of plant traits // Glob. Change Biol. 2011. Vol. 17. P. 2905–2935. DOI:10.1111/j.1365-2486.2011.02451.x.
35. Wang X., Zhang Y., Fraser R., Chen W. Evaluating the major controls on permafrost distribution in Ivvavik National Park based on process-based modelling // Proceeding of The 63rd Canadian Geotechnical Conference and the 6th Canadian Permafrost Conference (GEO2010) . Calgary, Canada, 2010. P. 1235–1241.
36. Tarnocai C. Soil Landscapes of the Firth River Area, Northwest Territories – Yukon Territory. Ottawa, Canada, Research Branch, Agriculture Canada, 1986. 1:1 million map.
37. Smith S.L. Burgess M.M. A Digital Database of Permafrost Thickness in Canada. Ottawa, Canada: Geological Survey of Canada, 2002. Open File № 4173. P. 38.
38. Burn C.R., Zhang Y. Permafrost and climate change at Herschel Island (Qikiqtaruq), Yukon Territory, Canada // J. Geophys. Res. 2009. Vol. 114. Paper F02001. DOI:10.1029/2008JF001087.
39. Burn C.R., Kokelj S.V. The environment and permafrost of the Mackenzie Delta area // Permafrost Periglac. 2009 . Vol. 20. P. 83–105. DOI:10.1002/ppp.655.
40. Bonnaventure P.P., Lewkowicz A.G., Kremer M., Sawada M.C. A permafrost probability model for the southern Yukon and northern British Columbia, Canada // Permafrost Periglac. 2012. Vol. 23. P. 52–68. DOI:10.1002/ppp.1733.
41. Goudriaan J. Crop Micrometeorology: a Simulation Study // Simulation Monographs. Wageningen, Netherlands: PUDOC, 1977. 236 p.
42. Wang S., Chen W., Cihlar J. New calculation methods of diurnal distribution of solar radiation and its interception by canopy over complex terrain // Ecol. Model. 2002. Vol. 155. P. 191–204.
43. Boland J., Ridley B., Brown B. Models of diffuse solar radiation // Renew. Energ. 2008. Vol. 33. P. 575–584.
44. Environment Canada and the National Research Council of Canada. Canadian Energy and Engineering Data Sets (CWEEDS) and Canadian Weather for Energy Calculations (CWEC). 2007. CD-ROM.
45. Ruth D.W., Chant R.E. The relationship of diffusive radiation to total radiation in Canada // Sol. Energy. 1976. Vol. 18. P. 153–154.