Mathematical modeling of Elbrus glaciers in the 21st century. Part 1. Glaciological model and setup of numerical experiments.

Full Text

ВВЕДЕНИЕ

Исследования оледенения Эльбруса и климатических условий в регионе имеют длительную, более чем столетнюю, историю (Золотарев, Харьковец, 2012). Накопленные к настоящему времени данные позволяют детально восстановить фактическую сторону современного состояния трёх десятков ледников Эльбруса, климатических условий на них, а также оценить масштабы и параметры изменения оледенения Эльбруса на протяжении десятилетий (Золотарев, 2009; Золотарев, Харьковец, 2012; Торопов и др., 2016; Kutuzov et al., 2019). Особо стоит отметить две фундаментальные монографии, опубликованные с интервалом в половину столетия (Оледенение Эльбруса, 1968; Ледники и климат Эльбруса, 2020), в которых собраны и обобщены сведения об этом уникальном природном объекте.

Объём оледенения Эльбруса составлял в 2017 г. 5.03±0.85 км³, 68% которого находилось ниже 4000 м над ур. моря со средней толщиной снежно-ледовой толщи 44.6±7.3 м (Kutuzov et al., 2019). Его площадь по состоянию на тот же год оценивалась в 112.20±0.58 км², что составляло немногим менее 10% площади оледенения всего Большого Кавказа (1132±0.58 км² по состоянию на 2014 г. (Tielidze, Wheat, 2018)).

Реки, берущие начало на склонах Эльбруса, относятся к бассейнам Чёрного и Каспийского морей. Учитывая, что доля ледникового стока в общем стоке рек северного склона Большого Кавказа оценивается в 1–9%, однако в тёплый период года на отдельных реках может достигать 40–45% (Лурье, Панов, 2013), можно предположить, что значительный по региональным меркам масштаб оледенения Эльбруса должен оказывать определённое влияние на формирование режима речного стока в регионе и косвенным образом на экосистемы.

В одном ледниковом массиве наблюдаются существенные региональные различия: ледники Эльбруса имеют разную экспозицию, оканчиваются в разных высотных зонах. В результате они теряют массу с разной скоростью (Золотарев, 2009; Золотарев, Харьковец, 2012; Kutuzov et al., 2019). Так, на южном склоне сокращение массы льда в среднем больше, чем на северном. Однако самые большие потери массы наблюдаются на ледовом плато Джикаугенкёз на северо-восточном участке, которое является частью ледников Чунгурчатчиран и Бирдждалычиран (Золотарев, Харьковец, 2012), а самые маленькие – на леднике Кюкюртлю, который выделяется среди прочих благодаря большой роли лавинного питания в аккумуляции (Kutuzov et al., 2019).

Исторические изменения параметров оледенения Эльбруса соответствуют общей тенденции сокращения горного оледенения мира в целом и на Кавказе в частности, хотя темпы этого сокращения, по-видимому, несколько различались. Однозначные выводы сделать трудно, поскольку приводимые оценки относятся по большей части к несовпадающим периодам времени. Так, площадь оледенения Эльбруса сократилась с 1997 по 2017 г. на 10.8% (0.51% в год). В то же время в 1971–2011 гг. на всём северном склоне Большого Кавказа площадь ледников сокращалась со скоростью 0.41% в год (Лурье, Панов, 2014), однако в 2000–2020 гг. скорость сокращения увеличилась до 1.15% в год (Котляков и др., 2023). Последняя оценка относится, правда, ко всему Большому Кавказу, однако, даже учитывая региональные различия, цифра сокращения площади превышала 1% в год. В тот же период площадь ледников Эльбруса сокращалась гораздо меньшими темпами, на 0.57% год (Котляков и др., 2023).

Начиная с 30-х годов XX века и до недавнего времени выполнялись фрагментарные метеонаблюдения разной продолжительности на склонах Эльбруса. Только в последнее десятилетие они стали носить относительно регулярный характер (Торопов и др., 2016), в частности на леднике Гарабаши, а также охватывали незначительные промежутки времени на Западном плато и Восточной вершине. В совокупности с длительными регулярными наблюдениями на сетевой метеостанции Терскол, расположенной у подножия Эльбруса в долине р. Баксан, появилась возможность реконструировать текущие климатические условия в регионе. Наличие данных наблюдений – это необходимое условие для коррекции региональных климатических прогнозов, основанных на математическом моделировании (Корнева и др., 2023). Что касается палеоклиматических реконструкций, то исходные данные для них были получены в ходе бурения толщи ледников Эльбруса. Среди кернов различной длины особо выделяется несколько кернов, полученных на Западном плато; длина самого глубокого из них превышает 180 м. Физико-химический анализ результатов глубокого бурения позволил восстановить ряды приземной температуры воздуха и скорости аккумуляции на протяжении нескольких последних столетий, а также и особенности циркуляции атмосферы в регионе (Ледники и климат Эльбруса, 2020).

Актуальность исследования возможной эволюции оледенения Эльбруса в меняющихся климатических условиях имеет важный аспект. Рецессия ледников приводит к активизации опасных явлений, связанных с ледниками. Отступание ледников сопровождается образованием большого объёма мёртвых и погребённых льдов, что зачастую служит источником селей вследствие прорыва заполненных водой термокарстовых каверн (Черноморец, 2005) или подпруженных мореной с погребённым под ней стагнирующим льдом прогляциальных озёр (Докукин, Хаткутов, 2016). На Эльбрусе наиболее опасным с точки зрения селей гляциально-ливневого генезиса считается северо-восточный склон. Так, на месте отступившего ледника Бирджалычиран были образованы озёра, большая часть которых расположена поверх мёртвого льда и подпружена последним (Petrakov et al., 2007).

Проблемы, обусловленные отступанием фронтов ледников и формированием областей мёртвого льда, обуславливают необходимость решения задачи моделирования формирования и трансформации поверхностной морены. Прогноз эволюции моренного покрова важен по нескольким причинам. Cлой морены толщиной более 5–7 см на леднике, выполняя изолирующую роль, задерживает его поверхностное таяние (Østrem, 1959; Поповнин и др., 2015) и, таким образом, замедляет скорость отступания фронта. Под моренным покровом возможно сохранение областей мёртвого льда, оставшегося после отступания ледника, которые могут служить подпрудой для формирующихся гляциальных озёр (Лаврентьев и др., 2020; Rounce et al., 2021).

Моренный покров на Эльбрусе сконцентрирован главным образом на ледниках Уллучиран, Большой Азау, Кюкюртлю, Битюгтюбе и ледовом плато Джикаугенкёз (рис. 1, а). В целом доля покрытого мореной льда мала (менее 5% в 2014 г.), однако с 1986 г. она увеличилась практически в 2.5 раза, что с учётом интервала неопределённостей соответствует скорости увеличения поверхностной морены на Большом Кавказе в целом (Tielidze et al., 2020). Максимальный рост доли заморененного льда наблюдался на восточном склоне Эльбруса.

 

Рис. 1. Эльбрус. (а) Ледники Эльбруса (1 – Уллучиран, 2 – Карачаул, 3 – Уллукол, 4 – Микелчиран, 5 – Джикаугенкёз, 6 – Ирикчат, 7 – Ирик, 8 – Теркол, 9 – Гарабаши, 10 – Малый Азау, 11 – Большой Азау, 12 – Кюкюртлю, 13 – Битюгтюбе); 1 – области, покрытые моренным чехлом по состоянию на 2000 г., 2 – очертания ледников в 2000 г. согласно RGI6.0. (б) Снимок со спутника Sentinel-2 (дата съёмки 01.09.2023). 1 – очертания ледников в 2000 г. согласно RGI6.0. (в) Толщина ледников на 2000 год (м) (данные Huss, Farinotti et al., 2012, обновлены в 2019 авторами); (г) Изменение высоты поверхности ледников с 2000 по 2019 г. (м) (данные Hugonnet et al., 2021)

Fig. 1. Elbrus. (a) Elbrus glaciers (1 – Ulluchiran, 2 – Karachaul, 3 – Ullukol, 4 – Mikelchiran, 5 – Dzhikiugankoz, 6 – Irikchat, 7 – Irik, 8 – Terkol, 9 – Garabashi, 10 – Maliy Azau, 11 – Bolshoy Azau, 12 – Kyukyurtlyu, 13 – Bityugtyube); 1 – areas covered by debris cover as of 2000, 2 – outlines of glaciers in 2000 according to RGI6.0. (б) Sentinel-2 satellite image (imaging date 01.09.2023). 1 – outlines of glaciers in 2000 according to RGI6.0. (в) Glacier thickness as of 2000 (data from Huss, Farinotti et al., 2012, updated in 2019 by the authors); (г) Change in glacier surface elevation from 2000 to 2019 (m) (data from Hugonnet et al., 2021)

 

Если исторические изменения параметров оледенения Эльбруса прямыми или косвенными методами изучены с большей или меньшей детализацией, начиная с периода формального окончания малого ледникового периода в середине XIX века (Золотарев, Харьковец, 2012), то долгосрочные прогностические количественные оценки его к настоящему времени отсутствуют. Были проведены исследования более глобального масштаба, относящиеся, в частности, к оледенению всего северного склона Большого Кавказа в целом (Лурье, Панов, 2014), но прогноз эволюции самих ледников Эльбруса до настоящего времени не выполнялся. Оценки составляющих баланса массы ледников Эльбруса в конце XIX столетия были получены в работе (Рыбак и др., 2019). Однако в ней не учитывались потенциальные изменения конфигурации ледников. Полноценных же количественных модельных оценок будущих очертаний, площадей и объёмов ледников Эльбруса в условиях различных сценариев климатических изменений до настоящего времени сделано не было. Цель настоящей работы заключается в том, чтобы заполнить этот пробел. Кроме того, на основе расчётов темпов отступания фронтов ледников нами решается задача определения времени вероятного формирования прогляциальных озёр и областей мёртвого льда, покрытых мореной. Численные эксперименты выполнялись на глобальной гляциологической модели GloGEMflow-debris, особенности архитектуры которой, настройка и валидация разбираются в первой части работы. В ней же приводится детальное описание использованных климатических сценариев.

Результаты моделирования и их обсуждение содержатся во второй части.

МЕТОДИКА ИССЛЕДОВАНИЯ

Модель GloGEMflow-debris. Для реализации прогностических расчётов мы использовали глобальную гляциологическую модель (ГГМ) GloGEMflow-debris (Postnikova et al., 2023). Эта модель представляет собой модификацию модели GloGEMflow (Zekollari et al., 2019), которая была дополнена блоком, описывающим трансформацию моренного покрова. ГГМ, их типичная архитектура, особенности настройки и применения детально разобраны в (Постникова, Рыбак, 2021, 2022), в связи с чем ниже приводятся только минимально необходимые для понимания настоящего исследования сведения, связанные, главным образом, с описанием работы блока моренного чехла.

В модели горный ледник представлен как одномерный объект, и все динамические процессы рассматриваются в 10 высотных поясах. При таком подходе Эльбрус представляется как совокупность тридцати отдельных ледников (Kutuzov et al., 2019), из них тринадцать крупных (см. рис. 1, а). Модель состоит из модулей (рис. 2), предназначенных для инициализации и прогностических расчётов. Инициализация модели служит для обеспечения внутренней согласованности геометрии ледника, его динамики, поверхностного баланса массы и моренного покрова. Она сводится к калибровке трёх групп ключевых параметров: параметра, который объединяет реологические свойства льда и условия базального (глыбового) скольжения (1), поправки поверхностного баланса массы (2) и параметров в описании динамики моренного покрова (3).

 

Рис. 2. Структура модели GloGEMflow-debris. Моренный блок выделен тонкой штриховой линией. Серым выделены основные блоки модели. Прямоугольниками обозначены используемые или рассчитанные данные, овалами – процедуры

Fig. 2. Structure of the GloGEMflow-debris model. The debris-cover module is highlighted with a thin dashed line. The main modules of the model are highlighted in gray. Rectangles indicate data used or calculated, ovals indicate procedures

 

Расчёт поверхностного баланса массы. Поверхностный баланс массы для исторического периода (1980–1990 гг.), периода валидации (1991–2020 гг.) и прогностического (2021–2100 гг.) рассчитывался с использованием индексно-температурного метода (Huss, Hock, 2015), с тем отличием, что параметры масс-балансового блока были откалиброваны индивидуально для каждого ледника для соответствия набору данных изменения объёма ледников в 2000–2019 годах полученных геодезическим методом (Hugonnet et al., 2021).

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

$M_{\text{total}}= M_{ice}\cdot \left(1-\frac{A_{debris}}{A}\right)+M_{debris}\cdot \frac{A_{debris}}{A}$, (1)

где $M_{debris}$ – абляция льда, покрытого мореной; $M_{ice}$ – абляция чистого льда; A – площадь льда; $A_{debris}$ – площадь морены в данной ячейке сетки. Для вычисления абляции под мореной используется функция зависимости таяния льда от толщины морены (Postnikova et al., 2023):

$M_{debris}=M_{ice}\cdot f_{debris}$, (2)

где $f_{debris}$ – коэффициент уменьшения или увеличения таяния в зависимости от толщины моренного покрова.

Для расчёта изменения скорости таяния под мореной в зависимости от её толщины используются так называемые кривые Острёма (Østrem, 1959). Это функции зависимости величины таяния от толщины моренного покрова, для которых характерно усиление таяния при росте толщины морены от нуля до так называемой эффективной толщины морены $h_{debris}^{eff}$, при которой достигается максимальное усиление таяния, как правило 1–3 см (Mattson et al., 1993; Verhaegen et al., 2024). Для ледников Эльбруса это неизвестная величина, которую мы принимаем за 2 см. При толщине морены больше $h_{debris}^{eff}$, скорость таяния льда под моренным чехлом экспоненциально уменьшается с увеличением толщины моренного слоя (Østrem, 1959). При этом толщина моренного чехла, начиная с которой таяние под мореной становится меньше, чем таяние чистого льда, называется критической, $h_{debris}^{crit}$ . Для Кавказа, по оценкам из работы (Поповнин и др., 2015), $h_{debris}^{crit}$ составляет 5–7 см. Поэтому коэффициент изменения скорости таяния $f_{debris}$ (3) при толщине моренного слоя не менее 7 см (Поповнин и др., 2015) рассчитывается в зависимости от толщины моренного покрова (Verhaegen et al., 2020) следующим образом:

$f_{debris}^{thick}=e^{\frac{-h_{debris}}{h_{debris}^{*}}},$ (3)

где $h_{debris}^{\text{*}}$ – характерная толщина моренного покрова, при которой фактическое таяние под мореной уменьшается до $e^{-1}$ или приблизительно 37% от скорости таяния чистого льда. Для ледников Эльбруса было взято значение $h_{debris}^{*}$= 44 см, рассчитанное для ледника Джанкуат (Verhaegen et al., 2024). Для более точной оценки величин $h_{debris}^{eff}$ , $h_{debris}^{crit}$, $h_{debris}^{*}$ для каждого ледника в отдельности требуется энергобалансовая модель для расчёта теплообмена с участием моренного покрова (Verhaegen et al., 2024), применение которой в рамках данного исследования нецелесообразно в силу малого распространения морены и недостатка входных характеристик моренного чехла на разных ледниках Эльбруса. Чувствительность модели к изменению параметра $h_{debris}^{*}$, определяющего форму кривой уменьшения таяния под мореной, описана в разделе «Оценка чувствительности модели».

Если же толщина моренного покрова в рассматриваемом высотном поясе не превышает 7 см, то применяется следующая расчётная схема для $f_{debris}$ : 1) при нулевой толщине морены $f_{debris}=1$ , то есть берётся исходный баланс массы, рассчитанный для чистого льда; 2) при толщине морены $h_{debris}^{eff}$= 2 см $f_{debris}^{eff}$ рассчитывается исходя из данных (Rounce et al., 2021), где для каждого ледника представлена величина таяния чистого льда и льда, покрытого мореной толщиной 2 см. Характерная величина коэффициента усиления таяния составляет около $f_{debris}^{eff}=1.38$; 3) при толщине морены от 0 до 2 см $f_{debris}$ линейно растёт от 1 до $f_{debris}^{eff}$; 4) при толщине морены от 2 до 7 см $f_{debris}$ линейно убывает от $f_{debris}^{eff}$ до $f_{debris}^{thick}\left(0.07\right)=0.94$.

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

 

Рис. 3. Ледник Уллучиран в 2024 г. (сценарий SSP2-4.5). (а) Баланс массы: 1 – чистого льда; 2 – при учёте моренного покрова и $h_{debris}^{*}=44$ см; 3 – при учёте моренного покрова и $h_{debris}^{*}=115$ см. (б) Коэффициент изменения таяния под моренным покровом: 4 – при $h_{debris}^{*}=44$ см; 5 – при $h_{debris}^{*}=115$ см. (в) Профиль ледника: 6 – моренный покров (толщина которого увеличена в 30 раз для наглядности); 7 – ледник при $h_{debris}^{*}=44$ см; 8 – ложе; 9 – поверхность ледника при $h_{debris}^{*}=115$ см; 10 – модельная площадь покрытости льда мореной $\frac{A_{derbis}}{A}$ ; 11 – толщина морены (м). Правая вертикальная ось относится, во-первых, к толщине морены (м), во-вторых, к $\frac{A_{derbis}}{A}$ (безразмерная величина)

Fig. 3. Ulluchiran Glacier in 2024 (SSP2-4.5 scenario). (a) Mass balance: 1 – clean ice; 2 – when debris cover is taken into account and $h_{debris}^{*}$ = 44 cm; 3 – when debris cover is taken into account and $h_{debris}^{*}$ = 115 cm. (б) Coefficient of melt correction under debris cover: 4 – when $h_{debris}^{*}=$ 44 cm; 5 – when $h_{debris}^{*}=$115 cm. (в) Glacier profile: 6 – debris cover (whose thickness is increased 30 times for clarity); 7 – glacier when $h_{debris}^{*}=$ 44 cm; 8 – bed; 9 – glacier surface when $h_{debris}^{*}=$ 115 cm; 10 – modeled area of debris-covered ice $\frac{A_{derbis}}{A}$; 11 – debris cover thickness (m). The right vertical axis refers firstly to the debris cover thickness (m), secondly to $\frac{A_{derbis}}{A}$ (dimensionless value).

 

Динамика ледника. Динамика ледника описывается законом сохранения массы и реологической зависимостью скорости ледника от напряжения:

$\left\{\begin{array}{c}\frac{\partial H}{\partial t}=-\nabla \left(\overline{u}H\right)+b\\ \overline{u}=\frac{2A}{n+2}{\tau }^{n}H\\ \tau =-\rho gH\frac{\partial s}{\partial x}\end{array}\right.$, (4)

где H – толщина ледника, $\overline{u}$ – вертикально усреднённая скорость, $\nabla \left(\overline{u}H\right)$ – дивергенция локального потока льда, b – поверхностный баланс массы, τ – движущее напряжение, n – экспонента в законе Глена, $\frac{\partial s}{\partial x}$ – наклон поверхности, A – коэффициент деформации-скольжения. Последний параметр объединяет эффекты базального скольжения и деформации льда, так как оба они определяются наклоном поверхности и толщиной льда, и, таким образом, имеют очень похожие пространственные характеристики (Zekollari et al., 2013). Ключевой пункт в настройке модели заключается в калибровке именно параметра А (см. ниже).

Моренный блок. Описание трансформации моренного чехла кардинально отличает модель GloGEMflow-debris от предыдущей версии GloGEMflow (Zekollari et al., 2013), поэтому целесообразно остановиться на этой модификации более подробно. Алгоритм описания эволюции моренного покрова был разработан на основании материалов работ (Anderson, Anderson, 2016) и (Verhaegen et al., 2020) и подключён к ГГМ GloGEMflow, как описано в статье (Postnikova et al., 2023).

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

$\frac{\partial h_{debris}}{\partial t}=-\frac{c_{debris}\min \left(\mathrm{0,}{b}_a\right)}{\left(1-{\varphi }_{debris}\right){\rho }_{debris}}-\frac{\partial \left(u_{surf}{h}_{debris}\right)}{\partial x}+I_{debris}\text{.}$ (5)

(I) (II) (III)

В уравнении (5) $h_{debris}$ – толщина моренного покрова, $t$ – время, $c_{debris}$ – концентрация моренного материала во льду, ${\varphi }_{debris}$ – пористость моренного материала, ${\rho }_{debris}$ – плотность пород моренного покрова, $b_a$ – годовой поверхностный баланс массы, $u_{surf}$ – скорость поверхности ледника, $I_{debris}$ – поступление/расход моренного материала. Цифрами I, II, III обозначены соответственно члены, ответственные за вытаивание, адвекцию и приходно-расходную компоненты. Плотность ${\rho }_{debris}$ = 2600 кг м^–3 характерна для гранитов и гнейсов. Значение ${\varphi }_{debris}$ на Кавказе было измерено только на леднике Джанкуат в Приэльбрусье (Bozhinsky et al., 1986) и задаётся как константа для ледников Эльбруса: ${\varphi }_{debris}$ = 0.43. Чувствительность модели к этому параметру низка по сравнению с геометрическими характеристиками моренного покрова (см. раздел «Оценка чувствительности модели») Поступление моренного материала на ледники Эльбруса происходит по большей части в результате экзарации (в результате взбрасывания донной морены по плоскостям сдвига – в этом случае можно наблюдать тонкую приязыковую морену). Также часть моренного материала поступает в верхних поясах некоторых ледников в результате эрозии склонов (в результате образуются срединные морены, как, например, на леднике Уллучиран). Эти процессы у нас задаются при помощи параметра концентрации моренного материала в толще льда $c_{debris}$, который настраивается для каждого ледника отдельно. Чем больше этот параметр, тем больше моренного материала вытаивает на поверхность ледника.

Моренный покров выносится в прифронтальную зону со скоростью $F_{debris}^{out}$ м в год:

$F_{debris}^{out}\left(x_{front}\right)=-c_{front}{h}_{debris}\left(x_{front}\right),$ (6)

где $c_{front}$ – константа, описывающая скорость удаления моренного материала с ледника в прифронтальную зону.

Изменение площади моренного чехла. Доля покрытой мореной площади вдоль линии потока параметризуется в зависимости от расстояния от фронта ледника $(x-x_{front})$:

$\frac{A_{debris}\left(x\right)}{A\left(x\right)} = G_A \cdot \left(a \cdot e^{\left(x-x_{front}\right)\cdot b}\right)$ (7)

где $A_{debris}$ – площадь льда, покрытого мореной, A – площадь ледника. $G_A={\alpha }_{debris}{\left(H_{debris}^{front}\right)}^{{\beta }_{debris}}$ – это коэффициент роста относительной площади морены, который обновляется каждый модельный год. $G_A$ зависит от толщины морены во фронтальной зоне $H_{debris}^{front}$, а параметры ${\alpha }_{debris}$ и ${\beta }_{debris}$ определяются путём процедуры калибровки, описанной ниже. Параметры a и b получаются из данных распределения моренного покрова в зависимости от расстояния от фронта ледника на 2001 год, которое приближается экспоненциальной функцией.

ИСХОДНЫЕ ДАННЫЕ

Конфигурация и толщина ледников Эльбруса. Данные гипсометрии ледников были получены из (Huss, Farinotti, 2012, обновлены авторами в 2019 г.). В качестве границ ледников за 2000 г. используются данные RGI 6.0 (RGI Consortium, 2017) (см. рис. 1, в). В целом указанные данные хорошо согласуются с данными, полученными на основе радиозондирования (Kutuzov et al., 2019)

Моренный покров. На первый взгляд влияние моренного покрова на внешний теплообмен ледников Эльбруса не может быть значительным, тем не менее, будучи сосредоточен во фронтальных частях ледников (см. рис. 1, а), он должен замедлять скорость таяния и, соответственно, тормозить их отступание. Для калибровки моренного блока использовались данные распространения моренного чехла по состоянию на 2001 г., а для его валидации – по состоянию на 2018 г. Соответствующие данные были получены из анализа спутниковых снимков Landsat 7 ETM+ и Sentinel-2. Относительная площадь покрытого мореной льда была рассчитана путём деления площади моренного покрова на площадь ледника в каждой ячейке пространственной сетки. После этого распределение относительной площади покрытого мореной льда в каждом высотном поясе вдоль профиля ледника было аппроксимировано экспоненциальной функцией в соответствии с уравнением:

$\frac{A_{debris}\left(x\right)}{A\left(x\right)}=\zeta e^{\left(x-x_{front}\right)\xi },$ (8)

где $x_{front}$ – положение фронта ледника, $A_{debris}\left(x\right)$ – площадь моренного покрова в ячейке сетки x, A(x) – площадь ледника в ячейке сетки x, а ζ и ξ – коэффициенты, которые были рассчитаны для каждого ледника отдельно.

Толщина моренного покрова была оценена методом инверсии температуры поверхности и таяния под мореной (см. рис. 1, а; Rounce et al., 2021). Предположительно, область моренного покрова на Эльбрусе, представленная в (Rounce et al., 2021), меньше, чем в реальности, так как в этой работе рассматривался моренный чехол на стагнирующих областях льда, где скорость льда не превышает 7.5 м в год. Из той же работы для каждого ледника были взяты значения величины таяния чистого льда и льда под слоем морены толщиной 2 см, для того чтобы примерно оценить эффект усиления таяния при наличии тонкого слоя морены.

Поверхностный баланс массы. Данные поверхностного баланса массы рассчитывались индексно-температурным методом (Huss, Hock, 2015). Параметры, используемые в расчётах, были откалиброваны в соответствии с набором данных изменений объёма ледников в 2000–2019 гг. геодезическим методом (Hugonnet et al., 2021).

Климатический форсинг в исторический период. Инициализация и валидация блока расчёта поверхностного баланса массы осуществляется с использованием данных приповерхностной температуры и количества осадков с месячным разрешением в 1979–2020 гг. из реанализа Европейского центра среднесрочного прогноза погоды (ERA-5; Hersbach, 2019). Пространственное разрешение данных составляет 31 × 31 км.

Климатические сценарии SSP. Для прогностического периода используются результаты тринадцати климатических моделей для пяти сценариев SSP (Shared Socioeconomic Pathways – общие социально-экономические пути развития) из проекта межмодельного сравнения CMIP6 (Eyring et al., 2016), с пространственным разрешением от 0.7° × 0.7° до 1.5° × 2°, в зависимости от климатической модели (табл. 1). Среди участников проекта CMIP6 были выбраны те модели, для которых были доступны данные месячного разрешения по крайней мере для сценариев historical, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5 (рис. 4, см. табл. 1). Исторический сценарий (historical) распространяется на 1850–2014 гг., доступен для всех климатических моделей и используется для сравнения данных ре-анализа и CMIP6 и коррекции последних.

 

Таблица 1. Климатические модели из проекта CMIP6, результаты которых были использованы для моделирования. Вертикальное разрешение (L) определяется количеством уровней

Table 1. Climate models from the CMIP6 project whose results were used for modeling. Vertical resolution (L) is determined by the number of levels

№	Название	Организация	Разрешение атмосферного блока (широта, долгота, количество уровней)
1	BCC-CSM2-MR	Пекинский климатический центр, Китай	1.125° × 1.125°, L46
2	CAMS-CSM1-0	Китайская академия метеорологических наук, Китай	1.13° × 1.12°, L31
3	CESM2	Национальный центр исследования атмосферы, США	1.25° × 0.94°, L32
4	CESM2-WACCM		1.25° × 0.94°, L70
5	EC-Earth3	Европейский консорциум национальных метеорологических служб и научно-исследовательских институтов, Европа	0.7° × 0.7°, L91
6	EC-Earth3-Veg		0.7° × 0.7°, L91
7	FGOALS-f3-L	Китайская Академия Наук, Китай	2.25° × 2°, L32
8	GFDL-ESM4	Геофизическая лаборатория гидродинамики NOAA, США	1° × 1.25°, L49
9	INM-CM4-8	ИВМ РАН, Россия	1.5° × 2°, L21
10	INM-CM5-0		1.5° × 2°, L73
11	MPI-ESM1-2-HR	Институт меторологии Макса Планка, Германия	0.9375° × 0.9375°, L95
12	MRI-ESM2-0	Метеорологический исследовательский институт, Япония	1.125° × 1.125°, L80
13	NorESM2-MM	Норвежский климатический центр, Норвегия	0.9375° × 1.25°, L32

 

Рис. 4. Аномалии климата в окрестностях Эльбруса (координаты 43.2 с.ш., 42.26 в.д.) по данным CMIP6 (без поправок в соответствии с высотным градиентом)

Fig. 4. Climate anomalies in the vicinity of Elbrus (coordinates 43.2 N, 42.26 E) according to CMIP6 data (without corrections according to the altitudinal gradient)

 

Согласованность между прошлыми климатическими данными и будущими климатическими сценариями из проекта CMIP6 была достигнута введением поправок для температуры и осадков (Huss, Hock, 2015; Compagno et al., 2022). Исторические данные CMIP6 на 1980–2010 гг. в каждой точке сетки были сравнены с данными реанализа в 1980–2010 гг. Для коррекции данных климатических моделей из CMIP6 в соответствии с данными реанализа, были рассчитаны, во-первых, аддитивные поправки для температуры: $T_{cor{r}_{fut}}^{m,Y}=T_{c{m}_{fut}}^{m,Y}+\left( \overline{T_{ERA}^m}-\overline{T_{c{m}_{hist}}^m}\right)$, где $\overline{T_{ERA}^m}$ – средняя температура за 1980–2010 гг. за месяц под номером m по данным реанализа, $\overline{T_{c{m}_{hist}}^m }$ – средняя температура для месяца m за 1980–2010 гг. по данным исторического сценария CMIP6, $T_{c{m}_{fut}}^{m,Y}$ – прогнозная температура на месяц m года Y.

Во-вторых, к данным климатических моделей из CMIP6 были применены дополнительные поправки температуры воздуха, чтобы учесть разницу в межгодовой изменчивости по сравнению с данными реанализа. Для каждого месяца m = [1, ..., 12] стандартное отклонение температуры за 1980–2010 гг. было рассчитано для данных реанализа $\left({\sigma }_{ERA}^m\right)$ и для данных каждой из климатических моделей $\left({\sigma }_{GCM}^m\right)$. Поправки к межгодовой изменчивости были заданы как ${\text{Ô}}^m=\frac{{\sigma }_{ERA}^m}{{\sigma }_{GCM}^m}$.

Таким образом, для каждого месяца m и года Y прогнозируемого периода данные межгодовой изменчивости температуры воздуха были скорректированы следующим образом (Huss, Hock, 2015): $T_{corr}^{m,Y}= \overline{T_{25}^{m,Y}}+{\text{Ô}}^m\cdot \left(T_{c{m}_{fut}}^{m,Y}-\overline{T_{25}^{m,Y}}\right)$, где $\overline{T_{25}^{m,Y}}$ – средняя температура за 25-летний период в окрестностях года Y по данным CMIP6. Суточные колебания были наложены на прогностические данные температуры месячного разрешения, с использованием статистики температурных колебаний из данных реанализа. Для коррекции осадков были вычислены мультипликативные поправки (аналогично Huss, Hock, 2015), средние за период 1980–2010 гг., на которые были домножены ряды из CMIP6: $P_{cor{r}_{fut}}^{m,Y}=P_{c{m}_{fut}}^{m,Y}\cdot \left(\overline{P_{ERA}^m}/\overline{P_{c{m}_{hist}}^m}\right)$ , где $P_{c{m}_{fut}}^{m,Y}$ – прогнозные осадки из CMIP6 на месяц m и год Y, $\overline{P_{ERA}^m}$ – средние осадки в период 1980–2010 гг. по данным реанализа, $\overline{P_{c{m}_{hist}}^m}$ – средние осадки в период 1980–2010 гг. по данным исторического сценария CMIP6. Все перечисленные поправки были наложены на проекции CMIP6 до 2100 г., в предположении, что в будущем они не изменятся.

Для вычисления изменения температуры с высотой были использованы данные температуры воздуха на уровнях с разным давлением (500 гПа, 600 гПа, 700 гПа, 850 гПа, 925 гПа, 1000 гПа) и соответствующие им геопотенциальные высоты для расчёта месячных значений температурного градиента в свободной атмосфере для ближайшего узла сетки. Осадки для каждого высотного пояса i для каждого ледника за месяц m корректировались следующим образом (Huss, Hock, 2015): $P_i^m=P_{cell}^m\cdot c_{prec}\cdot \left(1+\left(z_i-z_{ref}\right)\cdot \frac{dP}{dz}\right)$, где $P_{cell}^m$ (м) – осадки в ближайшей точке сетки климатической модели за месяц m; $c_{prec}$ (безразмерная величина) – поправочный коэффициент на количество осадков, который вычисляется для каждого ледника отдельно в ходе калибровки; $\frac{dP}{dz}$ (1/м)– высотный градиент осадков, равный 1% на 100 м; $z_{ref}$ – это высота над уровнем моря, соответствующая узлу сетки климатической модели. Для ледников с высотным диапазоном >1000 м осадки уменьшаются экспоненциально в верхних высотных поясах ледника (верхние 25%), чтобы учесть типично пониженное содержание влаги в воздухе и усиление ветровой эрозии по мере увеличения высоты (Huss, Hock, 2015).

НАСТРОЙКА МОДЕЛИ

Перед постановкой и проведением численных экспериментов математическая модель должна пройти через процедуры калибровки ключевых параметров и валидации. Калибровка ГГМ особенно важна в силу самой природы этого класса моделей, использующих значительные обобщения и упрощения энергообменных и динамических процессов. В нашем случае калибровка затрагивает масс-балансовый, динамический и моренный блоки ГГМ. Валидация – это, по сути дела, проверка того, насколько адекватно модель воспроизводит реальную динамику моделируемых объектов по заданным входным данным.

Калибровка масс-балансового блока ГГМ. Масс-балансовый блок модели был откалиброван для каждого ледника в отдельности в соответствии с данными изменения высоты поверхности ледников (Hugonnet et al., 2021) (рис. 1г). Допустимая точность – 0.01 м в.э. год^–1. Процедура калибровки была применена к трём ключевым параметрам: мультипликативный поправочный коэффициент на количество осадков $c_{prec}$ (варьирует между 0.6 и 2), который выполняет функцию адаптации климатических данных к специфическим особенностям конкретного ледника (местные топографические эффекты, дождевая тень и т.д.) (1); коэффициенты градусо-дней (degree-day factor – DDF) для снега $DD{F}_{snow}$ (варьирует между 1.75 и 4.5 мм ∙ К/день) и льда $DD{F}_{ice}$ (в два раза больше, чем $DD{F}_{snow}$), которые переводят количество дней с положительной температурой в количество растаявшего снега или льда (2); температурная поправка, нивелирующая недостаточное пространственное разрешение входных климатических данных (3). Сама процедура состоит из трёх шагов: сначала калибруется параметр коррекции осадков; затем, если отклонения от данных баланса массы остаются большими, калибруется параметр DDF (по умолчанию $DD{F}_{snow}=3 \text{мм}\cdot К/\text{день}$ и $DD{F}_{snow}=6 \text{мм}\cdot К/\text{день}$); если второй шаг не даёт достаточно хорошего результата, систематически смещается параметр коррекции температуры $t_{off}$. Более подробно процедура калибровки разобрана в (Huss, Hock, 2015).

В ходе проведения численных экспериментов значения калибруемых параметров (представленные в табл. 2) корректируются для соответствия с эволюционирующим моренным чехлом.

 

Таблица 2. Параметры масс-балансового блока, полученные в результате калибровки, для основных ледников

Table 2. Mass-balance module parameters obtained as a result of calibration for the main glaciers

№	Название ледника	Параметры масс-балансового блока				Параметры динамического блока	Параметры моренного блока
		$DD{F}_{snow}$ (мм день–1K–1)	$DD{F}_{ice}$(мм день–1K–1)	$c_{prec}$	$t_{off}$ K	$A$ (Па–3год–1)	$c_{debris}$(кг м–3)
1	Уллучиран	4.8	9.6	1.6	2.83	1.00e–16	10.85
2	Карачаул	4.5	9	1.6	3.56	1.70e–16	5.07
3	Уллукол	4.49	8.97	1.6	3.37	1.91e–16	19.89
4	Микелчиран	4.49	8.98	1.6	2.17	2.39e–16	26.2
5	Джикаугенкёз	4.72	9.44	1.6	1.62	9.93e–17	4.23
6	Ирикчат	4.49	8.99	1.6	1.57	1.17e–18	1.05
7	Ирик	4.55	9.11	1.6	1.86	1.99e–16	1.05
8	Терскол	4.50	9.00	1.6	2.00	2.12e–16	–
9	Гарабаши	4.49	8.98	1.6	2.72	2.35e–16	13.71
10	Малый Азау	4.80	9.59	1.6	3.37	2.14e–16	4.58
11	Большой Азау	4.51	9.03	1.6	1.37	2.17e–16	7.74
12	Кюкюртлю	4.8	9.649	1.6	2.98	2.21e–16	1.05
13	Битюгтюбе	4.44	8.87	1.6	1.37	3.57e–16	17.18

 

Калибровка моренного и динамического блоков. Динамика морены и ледника в модели связаны, поэтому их калибровка проводится одновременно и проходит в три этапа (рис. 5).

 

Рис. 5. Этапы калибровки модели GloGEMflow-debris

Fig. 5. Stages of the GloGEMflow-debris model calibration

 

Этап I. Численные эксперименты начинаются с условия отсутствия ледникового покрова. Далее возникающий ледник эволюционирует под действием рассчитанного поверхностного баланса массы, который соответствует климату 1981–1990 гг., пока не будет достигнуто стационарное состояние, на что требуется приблизительно от ста до шестисот модельных лет. Затем ледник развивается с условного 1990 г. до даты инвентаризации RGI (2000 год для Эльбруса) без эволюционирующей морены. Коэффициент A (см. выше) калибруется так, чтобы объём модельного ледника соответствовал объёму ледника из RGI с отклонением не более 1%. Одновременно корректируется расчётный поверхностный баланс массы так, чтобы длина модельного ледника не отклонялась от зафиксированной в RGI более чем на 1%.

Этап II. На этом этапе к модели подключается моренный блок, параметры которого калибруются для соответствия реальному распределению моренного покрова по состоянию на 2001–2004 гг. Для каждого ледника, на котором он присутствует (Большой Азау, Малый Азау, Уллучиран, Гарабаши, Битюгтюбе, Ирик, Ирикчат, Кюкюртлю, Терскол, Карачаул, плато Джикаугенкёз), настраиваются три параметра: $F_{debris}^{input}$ – скорость поступления моренного материала (1); параметры ${\alpha }_{debris}$ (2) и ${\beta }_{debris}$ (3), регулирующие скорость роста относительной площади покрытого мореной льда $G_A$ в зависимости от толщины моренного чехла на языке ледника:

$G_A={\alpha }_{debris}{\left(H_{debris}^{front}\right)}^{{\beta }_{debris}}$. (9)

Параметры ${\alpha }_{debris}$ (2) и ${\beta }_{debris}$ подбираются таким образом, чтобы среднеквадратическая ошибка (с.к.о.) воспроизведения относительной площади покрытого мореной льда во всех диапазонов высот, где присутствует морена, не превышала 0.01.

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

Калибровка параметра концентрации моренного материала в толще ледника. Поступление моренного материала, как уже упоминалось, происходит в результате вытаивания в зоне абляции. Параметр концентрации морены в толще ледника был откалиброван с использованием модельных данных о толщине морены на 2008 год (см. рис. 1, а; Rounce et al., 2021). Цель процедуры калибровки заключается в достижении соответствия средней по всему леднику толщины моренного покрова в 2008 году.

Этап III. На этом этапе производится повторная калибровка параметров динамического и масс-балансового блоков, которая обусловлена введением моренного чехла. В динамическом блоке корректируется параметр А, до достижения соответствия площади и объёма модельных ледников данным RGI6.0 и (Huss, Farinotti et al., 2012, обновлены авторами в 2019 г.), соответственно.

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

ВАЛИДАЦИЯ МОДЕЛИ

Валидация (рис. 6–8) необходима для проверки слаженности работы настроенных блоков модели. В ходе проведения численного эксперимента происходит тесное взаимодействие различных модельных блоков. Общая картина выглядит следующим образом. Во-первых, изменение климатических условий преобразуется в изменения в удельном балансе массы для каждого 10-метрового высотного пояса ледника, что, в свою очередь, влияет на расчёт толщины льда. Во-вторых, на поверхностный баланс массы влияет моренный покров в зависимости от его площади и толщины, эволюция которых воспроизводится в моренном блоке и которые зависят от и динамики ледника. В-третьих, в динамическом блоке рассчитывается вертикальное (вдоль линии тока льда) перераспределение объёма льда и моренного покрова, а также обновляется их распределение во времени.

 

Рис. 6. Валидация. (а) Наблюдаемые (в 2000 г. по данным RGI6.0 и в 2018 г. по данным Хромовой и др., 2021) и модельные высоты фронтов ледников Эльбруса в 2000 г. (год калибровки) и в 2018 г. (год валидации), м над ур. моря 1 – Большой Азау, 2 – Ирик, 3 – Терскол. (б) Наблюдаемая (Kutuzov et al., 2019) и модельная потеря массы ледников Эльбруса в 1997–2017 гг. 4 – Джикаугенкёз, 5 – Кюкюртлю. (в) Измеренный (Kutuzov et al., 2019) и модельный объём ледников в 2017 и в 1997 г., а также разница между ними. (г) Измеренное (Kutuzov et al., 2019) и (д) модельное сокращение массы ледников в 1997–2017 гг. Границы ледников показаны по состоянию на 2000 г.

Fig. 6. Validation. (a) Observed (in 2000 based on RGI6.0 data and in 2018 based on data from Khromova et al., 2021) and modeled heights of Elbrus glacier fronts in 2000 (calibration year) and in 2018 (validation year), m a.s.l. 1 – Bolshoi Azau, 2 – Irik, 3 – Terskol. (б) Observed (Kutuzov et al., 2019) and model mass loss of Elbrus glaciers in 1997–2017. 4 – Dzhikiugankyoz, 5 – Kyukyurtlyu. (в) Measured (Kutuzov et al., 2019) and modeled glacier volume in 2017 and in 1997, and the difference between them. (г) Measured (Kutuzov et al., 2019) and (д) modeled glacier mass reduction in 1997–2017. Glacier boundaries are shown as of 2000.

 

Рис. 7. Баланс массы ледника Гарабаши: 1 – модельный, 2 – наблюдённый (WGMS, 2022)

Fig. 7. Mass balance of the Garabashi Glacier: 1 – modelled, 2 – observed (WGMS, 2022)

 

Рис. 8. (а) Модельные и наблюдённые скорости движения основных ледников Эльбруса: 1 – модельные, 2 – наблюдённые. (б) Скорости движения льда в ледниках Эльбруса на 2018–2019 гг., в соответствии с данными (Millan et al., 2022)

Fig. 8. (a) Modeled and observed velocities of the main Elbrus glaciers: 1 – modeled, 2 – observed. (б) Ice velocities of Elbrus glaciers for 2018–2019, according to (Millan et al., 2022)

 

Как уже было упомянуто выше, масс-балансовый блок был откалиброван по данным (Hugonnet et al., 2021) за период 2000–2019 гг., в связи с чем его валидация заключалась в проверке модельных результатов условного модельного периода 1997–2017 гг. на соответствие оценкам (Kutuzov et al., 2019), которые также были получены геодезическим методом, но по другим данным (изображения Ames Stereo Pipeline Обсерватории ледников Pléiades на 2017 г. и стереопары, полученные из данных аэрофотосъёмки географического факультета МГУ, которая была выполнена в 1997 г.). В целом, сравнение модельных результатов с рассчитанными, исходя из данных наблюдений, показывают удовлетворительное соответствие: модельный баланс массы всего ледникового массива сокращается в среднем со скоростью около –0.51 м в.э. в год в обоих случаях. Некоторые различия наблюдаются при сравнении динамики баланса массы отдельных ледников (рис. 6, б, г, д), хотя и в этом случае расхождения некритичны (с.к.о. = 0.3 м в.э. год^–1, R² = 0.71, см. рис. 6, б). Они обусловлены тем, что процессы энергообмена в ГГМ не могут, естественно, охватить весь спектр реальных процессов, и региональные различия неизбежны. Заметим, что результаты модельных реконструкций показывают, что в начале 2000-х годов поверхностный баланс массы оледенения Эльбруса был в среднем положителен, что качественно соответствует упомянутому в работе Золотарева (2013) наступанию некоторых ледников на Эльбрусе.

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

Работа динамического блока была оценена с помощью спутниковых данных (Millan et al., 2022) за 2018 год (см. рис. 8, б). Учитывая упрощения динамического блока (моделирование вдоль осевой линии), модель удовлетворительно передаёт динамику основных ледников Эльбруса (см. рис. 8, а). Исключения представляют, во-первых, ледники Уллукол и Микелчиран, для которых границы RGI6.0 не соответствуют направлению потока льда, а во-вторых, ледник Уллучиран, который в реальности разделён на два потока.

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

Значения высоты фронта ледника над уровнем моря сравнивались в 2018 г. (модельные на 2018 год и наблюдённые, согласно данным (Хромова и др., 2021; см. рис. 6, а). Первая из проблем валидации по данным (Хромова и др., 2021) о состоянии ледников Кавказа на 2018 г. состоит в том, что эти данные для некоторых ледников могут быть несопоставимы с данными RGI6.0 за 2000 г. К примеру, сравнение положения фронтов ледников за 2000 г. (RGI6.0) и 2018 г. (Хромова и др., 2021), показывает, что ледники Ирик и Джикаугенкёз наступали в 2000–2018 гг., но это не соответствует действительности (Shahgedanova et al., 2014), а Большой Азау отступал со скоростью 10 м/год, в то время как модельный результат составляет 28 м/год, что более реалистично по сравнению с прошлыми исследованиями (22 м/год до 2010 г. по результатам из (Shahgedanova et al., 2014). Однако в целом модельная высота фронта на год валидации хорошо соответствует наблюденной (коэффициент детерминации R²=0.89, среднеквадратичное отклонение RMSE=103 м; см. рис. 6, а).

Модельная конфигурация ледников в 2017 г. сравнивалась с границами ледников, в работе (Kutuzov et al., 2019). Учитывая сложность корректного определения границ ледников, полностью покрытых мореной у фронта, можно предположить, что погрешности в оконтуривании ледников как в год инициализации (2000) по RGI6.0, так и в год валидации (2017) по (Kutuzov et al., 2019), сопоставимы с погрешностями моделирования до 2017 г. Несмотря на то, что объём некоторых ледников (например, Джикаугенкёз, Большой Азау) недооценен моделью, по сравнению с данными (Kutuzov et al., 2019), как в 1997 г, так и в 2017 г. (см. рис. 6, в), в целом достигается хорошее соответствие (коэффициент детерминации не меньше 0.9, среднеквадратичное отклонение не больше 0.4 км³, см рис. 6, в). Причина некоторого расхождения модельных результатов и оценок (Kutuzov et al., 2019), очевидно, заключается в особенностях процедуры инициализации. В частности, толщина ледников, используемая нами (Huss, Farinotti et al., 2012, обновлены в 2019 г.), немного завышена, в особенности в области аккумуляции выше 4600 м. В этой части оледенения Эльбруса расхождение в объёмах льда может сохраняться до конца века и приводить к несколько завышенным модельным оценкам. Модельное сокращение объёма ледников хорошо соответствует наблюдённому (коэффициент детерминации 0.95, среднеквадратичное отклонение модельного изменения объёма от измеренного составляет 0.064 км³, см рис. 6, в). Таким образом, результаты валидации модели можно считать удовлетворительными.

ОЦЕНКА ЧУВСТВИТЕЛЬНОСТИ МОДЕЛИ

Параметры моренного покрова. Характерная толщина моренного покрова $h_{debris}^{*}$ определяет скорость убывания таяния с ростом толщины морены. Чем больше $h_{debris}^{*}$, тем сильнее изолирующие свойства морены и тем меньше толщина моренного чехла, требующаяся для полного прекращения таяния. Для оценки чувствительности модели к параметру $h_{debris}^{*}$ были произведены расчёты при $h_{debris}^{*}=115$ см (Verhaegen et al., 2020; Postnikova et al., 2023). Результаты показали, что уменьшение $h_{debris}^{*}$ в 2.6 раза приводит к задержке отступания фронта ледника (см. рис. 3) и увеличению времени существования мёртвых льдов, образованных при отступании ледника, однако в среднем в 2030-x годах различия в модельной геометрии ледника при разных значениях $h_{debris}^{*}$ нивелируются. Влияние этого параметра на объём ледников не превышает 4%.

Так как для моделирования ледников Эльбруса была выбрана модификации модели трансформации моренного покрова, в которой основной член уравнения (5) – вытаивание, то константы, входящие в компоненту вытаивания $({\rho }_{debris}, {\varphi }_{debris})$ , могут иметь существенное значение. Для сравнения, в том случае, когда основным источником моренного материала на поверхности ледника принимается эрозия склонов и перераспределение этого материала ниже по течению (Postnikova et al., 2023), компонента вытаивания и, соответственно, плотность и пористость морены играют незначительную роль по сравнению с толщиной морены, накопленной в результате сноса обломочного материала.

В силу того, что плотность моренного покрова на Эльбрусе неизвестна, были проведены эксперименты на чувствительность результатов к этому параметру, путём увеличения или уменьшения его на 15%: ${\rho }_{debris}=3000 кг/{м}^3$ и ${\rho }_{debris}=2200 кг/{м}^3$. Чувствительность модели к этому параметру может быть существенной: до 5% разницы в объёме льда $(100\%\cdot \frac{V_{modif}}{V})$ до 2070 г. и до 10–15% после 2070 г. (в частности в силу того, что общий объём ледников V значительно уменьшается к концу века). Схожие результаты получились при изменении пористости морены ${\varphi }_{debris}$ на 15%.

В модель трансформирующегося моренного чехла из работы (Postnikova et al., 2023) мы ввели эффект усиления таяния в результате наличия тонкого слоя морены. Заметим, что толщина моренного покрова и его площадь относительно «чистого» льда убывают от фронта ледника к границе питания (см. рис. 3, в). Таким образом, зона «тонкой» морены находится в тех высотных поясах, где доля покрытости льда мореной минимальна. Поэтому в формуле (1) таяние на участках с тонкой мореной $M_{debris}$ усилено, но площадь распространения таких участков относительно чистого льда в том же высотном поясе $\frac{A_{debris}}{A}$ невелика. Поэтому в результате градиент баланса массы при учёте усиления таяния и без него почти одинаков. Влияние эффекта усиления таяния на объём льда составляет не более 2%.

ЗАКЛЮЧЕНИЕ

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

Модель GloGEMflow-debris была настроена для условий ледников Эльбруса. В частности, модуль эволюции морены был адаптирован так, чтобы источник поступления отложений на поверхность ледника был не коллювиальным (обвалы, денудация склонов), а экзарационным (взбрасывание донной морены по плоскостям сдвига).

Валидация модели показала, что в начале XXI века в численных экспериментах систематически несколько недооценивалась потеря массы по сравнению с данными 1997–2017 гг. из работы (Kutuzov et al., 2019). Тем не менее, модель успешно воспроизводит общие закономерности потери массы, а также кинематику ледников, подтверждая её применимость в контексте исследований изменений ледникового покрова на Эльбрусе.

Благодарности. Исследование выполнено при финансовой поддержке РНФ в рамках научного проекта 23-27-00050. Авторы выражают искреннюю признательность С.С. Кутузову (Институт географии РАН) за предоставленные данные по конфигурации и толщине ледников Эльбруса. Авторы очень благодарны рецензентам за ценные замечания, которые позволили значительно улучшить статью.

Acknowledgements. This work was supported by the Russian Science Foundation, grant 23-27-00050. The authors express sincere gratitude to S.S Kutuzov (Institute of Geography of RAS) for providing data on configuration and thickness of Elbrus glaciers. The authors are grateful to the reviewers, whose valuable comments allowed improving the original paper.

About the authors

T. N. Postnikova

Water Problems Institute of RAS; Department of Geography, Lomonosov Moscow State University

Author for correspondence.
Email: tasinidze@gmail.com
Russian Federation, Moscow; Moscow

O. O. Rybak

Water Problems Institute of RAS; Institute of Natural and Technical Systems; Vrije Universiteit Brussel

Email: tasinidze@gmail.com

Earth System Science and Departement Geografie, Vrije Universiteit Brussel

Russian Federation, Moscow; Sevastopol; Brussels, Belgium

A. S. Gubanov

Department of Geography, Lomonosov Moscow State University

Email: tasinidze@gmail.com
Russian Federation, Moscow

H. Zekollarie

Vrije Universiteit Brussel

Email: tasinidze@gmail.com

Department of Water and Climate

Belgium, Brussels

M. Huss

ETH Zürich; Swiss Federal Institute for Forest, Snow and Landscape Research (WSL); University of Fribourg

Email: tasinidze@gmail.com

Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zürich; Department of Geosciences, University of Fribourg

Switzerland, Zürich; Birmensdorf; Fribourg

References

Supplementary files