ВВЕДЕНИЕ

Потепление климата в высоких широтах в течение 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].

Для представляемого в данной статье исследования были использованы три основных источника полевых данных. Первый – это отчет службы «Парки Канады» о ресурсах НПИ [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)) оценивалась по следующей формуле:

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

Среднемесячные значения температуры воздуха и общие месячные количества осадков были преобразованы в суточные данные с использованием вышеописанного метода (см. формулы (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% (на основе уравнения Стефана).

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

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

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

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

С прогрессирующим увеличением глубины летнего протаивания глубокий слой, оттаявший летом, может не замерзнуть следующей зимой и образовать талик (круглогодично незамерзаюший слой) между многолетней мерзлотой и слоем сезонного протаивания. В настоящее время многолетняя мерзлота подстилает почти все участки суши в парке и лишь очень небольшая доля площади содержит талики (в 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) даже без учета данных для упомянутых выше семи участков с наибольшими различиями. Низкая корреляция в основном обусловлена вариациями грунтовых условий в пределах экосистемы того или иного типа. Поэтому для картирования многолетней мерзлоты с высоким разрешением необходимы более детальные в пространственном отношении карты грунтовых условий.

На основе полевых наблюдений автор работы [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).

В работе [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).

Во-вторых, хотя при некоторых исследованиях была закартирована многолетняя мерзлота в горных регионах с высоким пространственным разрешением и было рассмотрено влияние топографических особенностей на инсоляцию (например, [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]:

где 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). Редактировал статью (исходную статью на английском языке. – Ред.) Т. Чжан.