Вилюйская синеклиза по высокой концентрации углеводородных природных ресурсов и их экономической значимости принадлежит к числу важнейших регионов России. В настоящее время в пределах Вилюйской синеклизы открыты 11 газовых и газоконденсатных месторождений, наиболее крупные из которых находятся в промышленной разработке. Помимо углеводородных месторождений, синеклиза относится к Ленскому угольному бассейну с большим количеством месторождений угля, а западная часть структуры богата месторождениями каменной соли [1].
Тепловое состояние горных пород и изменения криогенной толщи во многом определяют особенности формирования месторождений и освоения природных ресурсов в северных регионах. В связи с этим изучение особенностей и закономерностей формирования геотемпературного поля криолитозоны рассматриваемой территории является важной научно-практической задачей.
Поздний неоплейстоцен, включая его финальную эпоху голоцен, является важнейшим историческим интервалом, последовательность событий которого полностью определила современный рельеф, ландшафт и особенности криолитозоны не только рассматриваемого региона, но и всей Земли. Поэтому развитие фундаментальных представлений об эволюции криолитозоны неразрывно связано с базовыми палеогеографическими гипотезами – последовательностью ландшафтно-климатических изменений в четвертичное время с акцентом на последние 125 тыс. лет. Палеогеографические гипотезы носят универсальный характер, сквозь их призму происходит формирование современных представлений о структуре и генезисе современных климата, ландшафтов, гидросферы, криосферы и биосферы.
Представления о палеогеографических изменениях территории всегда были и остаются теоретической основой для интерпретации геотемпературного поля современной криолитозоны, наблюдаемой в различных геологических и орографических структурах. Тепловое поле криогенной толщи является в значительной мере реликтом – результатом ее длительного охлаждения во время последних криохронов, что позволяет использовать этот факт для анализа состояния природной среды.
Известно, что скорость и глубина промерзания грунтов определяются многими причинами, среди которых главными всегда считались климат, условия теплообмена на поверхности, теплопроводность и льдистость пород. Эти же факторы, предположительно, оказывают влияние на инертность криолитозоны, ее устойчивость к климатическим колебаниям различных циклов.
При интерпретации аномалий температурного поля, которое унаследовано от предыдущих криохронов и является в целом нестационарным современному климату, необходимо учитывать не только вышеперечисленные признаки, но локальные палеогеографические обстановки его формирования. К таким условиям рассматриваемого региона относятся микроклимат, растительные и почвенные покровы, озерно-аллювиальные и эолово-криогенные процессы. Причем последние оказывали значительно большее влияние в регионе, нежели имеющиеся на сегодняшний день представления большинства исследователей.
Вопросам изучения геотемпературного поля различных тектонических структур посвящены работы: Е.А. Любимовой, Я.Б. Смирнова, Б.Г. Поляка, М.Д. Хуторского, Lachenbruch, Sass, Vitorello, Pollack, и многих других, а по районам Сибири – А.Д. Дучкова, В.Т. Балобаева, С.В. Лысак, В.Н. Девяткина, А.И. Левченко, Б.В. Володько, М.Н. Железняка и др. Оценка полученных данных по тепловому потоку и теплофизическим свойствам горных пород говорит о специфичности этих параметров в рамках отдельных структур.
Геотемпературное поле и мощность многолетнемерзлой толщи определяется тремя основными факторами – температурой на поверхности пород, теплофизическими свойствами разреза и внутриземным тепловым потоком. Именно эти параметры, наряду с климатом и условиями теплообмена на поверхности, необходимы для выполнения достоверных прогнозов и палеореконструкций при решении научных фундаментальных и прикладных задач освоения территории.
Изучение глобальных изменений природной среды в будущем производится на основе ретроспективных схем развития ландшафтов в прошлом. Этапы развития ландшафтов в позднем плейстоцене и голоцене служат аналогами прогнозируемых изменений в будущем.
Метод реконструкции температурной истории земной поверхности по данным геотермии скважин прочно вошел в арсенал методов палеоклиматических исследований. Анализируя современное распределение температур горных пород с глубиной по замерам в скважинах, можно оценить колебания среднегодовой температуры земной поверхности в прошедшие геологические периоды. Эти колебания, связанные с климатическими вариациями, распространяются через толщу горных пород и вызывают возмущение стационарного геотемпературного поля, формируемого глубинным тепловым потоком и тепловыми свойствами среды (температурную аномалию). Годовые колебания не проникают глубже 20÷30 м, колебания столетнего масштаба оставляют аномальный след в интервале первых сотен метров, а потепление в конце последней ледниковой эпохи (~10 тыс. лет назад) фиксируется температурной аномалией до глубины 1,5÷2 км. Плейстоцен-голоценовое потепление представляет особый интерес. Оно завершает ближайший к нам цикл ледниково-межледниковых колебаний, характерных для климата четвертичного периода. Несмотря на более чем столетнюю историю изучения этого вопроса, механизм ледниковых колебаний во многом остается неясным [2].
Основная группа геокриологов России до недавнего времени считала, что многолетнее промерзание пород (особенно в южных районах Сибири) началось только в раннем неоплейстоцене, так как оно было обусловлено общепланетарным похолоданием климата на границе неогена и четвертичного периода. В то же время на палеоклиматической кривой, построенной для Центральной Якутии на основе комплексного изучения опорного разреза Мамонтова гора, показано, что в конце плиоцена средняя годовая температура воздуха Тв была порядка +12 °С, а многолетнее промерзание пород началось только в середине раннего плейстоцена [3]. Некоторые ученые считают, что многолетнее промерзание пород началось только в среднем плейстоцене, синхронно с эпохой максимального самаровского оледенения, а другие полагают, что во многих районах криогенная толща мерзлых пород имеет позднеплейстоценовый возраст [4].
Материалы и методы исследования
Определение геотемпературных полей и криогенной толщи конкретного региона, помимо сведений о сценарии палеоклимата в данном регионе, о теплофизических характеристиках горных пород геологического разреза включает в себя математическую постановку соответствующей задачи и выбор адекватного этой постановке метода ее решения. Эта постановка должна обязательно включать в себя описание фазового перехода «вода – лед» и различие теплофизических свойств пород в талом и мерзлом состояниях.
Выбор метода решения соответствующей обратной задачи будет определяться граничными условиями, то есть изменениями температуры на поверхности земной коры, которые зависят от выбранного сценария палеоклимата, но в любом случае будут циклическими функциями времени. Это сразу заставляет отказаться от методов с явным выделением неизвестной границы фазового перехода, так как при этом таких границ может быть несколько, их движение будет во времени немонотонно, а некоторые из них могут сливаться друг с другом или исчезать.
Понятно, что такие задачи могут быть решены только численными методами, из которых наиболее подходящим является метод Самарского – Моисеенко, или, иначе, метод сквозного счета со сглаживанием разрывных коэффициентов уравнения теплопроводности в окрестности температуры фазового перехода. Этот метод, в котором граница раздела фаз явно не выделяется, позволяет использовать однородные разностные схемы. При этом скрытая теплота фазового перехода вводится с применением d-функции Дирака в коэффициент теплоемкости – как сосредоточенная теплоемкость. Получаемая таким образом разрывная функция затем «размазывается» по температуре и не зависит от размерности задачи и количества фазовых границ. Это равносильно предположению, что фазовый переход происходит не при одной определенной температуре, а в некотором интервале температур, который определяется величиной параметра сглаживания («размазывания»). Выбранный метод лег в основу математической модели для расчетов температурного поля и динамики криогенной толщи при заданной температурной истории в конкретных геологических условиях.
Целью данной работы является моделирование геотемпературного поля и динамики многолетнемерзлой толщи для различных сценариев изменения климата и условий теплообмена применительно к различным областям Вилюйской синеклизы. В модели используется двухфазная задача Стефана [5], с помощью которой выполняются расчеты температурного поля и динамики криогенной толщи при заданной температурной истории в конкретных геологических условиях. На первом этапе выбирается и тестируется наиболее надежный метод решения прямых задач, так как существующие алгоритмы решения обратных задач, то есть задач реконструкции температурной истории земной поверхности по результатам геотермических исследований, сводятся к последовательности их решения.
На верхней границе массива горных пород в качестве сценариев изменения климата рассмотрены: 1) кривая изменения относительной температуры на границе льдообразования в Антарктиде [6] и 2) геокриологические интерпретации палеоклиматической летописи [4, 7] на основе комплексных исследований содержания створок диатомовых водорослей и биогенного кремнезема в донных осадках озера Байкал. На нижней границе массива приняты два условия: постоянство температуры или постоянство внутриземного теплового потока.
Расчеты выполнялись при усредненных значениях теплофизических свойств осадочных горных пород, соответствующих Средневилюйской площади. Теплопроводность горных пород в талом и мерзлом состояниях с учетом их влагонасыщения, а также их плотность были взяты из работ [1, 8]. Влажность песчанистых пород принята по работе [9], удельные теплоемкости в сухом состоянии и плотность внутриземного теплового потока – по [1].
Результаты исследования и их обсуждение
В начальном приближении был проведен вычислительный эксперимент, в котором считалось, что промерзание горных пород началось 1,3 млн лет назад (середина эпохи плейстоцена). Расчет был проведен при следующих исходных данных: температура дневной поверхности до начала ее изменения T0 = 0 °C, среднегодовая температура в наше время Tt = –10 °C, температура горных пород TL = 83 °C на глубине L = 3000 м; геотермические градиенты в основных типах горных пород отложений Вилюйской синеклизы Гr1 = 0,025 °С/м, Гr2 = 0,030 °С/м, Гr3 = 0,027 °С/м. На верхней границе задается линейная функция температурной истории 
, где (Tt – T0) – полная амплитуда изменения температуры на момент tmax = 1,3 млн лет. Получено, что при задании граничного условия I-го рода на нижней границе за 500 тыс. лет образовалась мерзлота мощностью 149 м, а при задании граничного условия II-го рода – 151 м. За это время температура дневной поверхности стала T0 = –3,85 °C. Тогда при задании граничного условия I-го рода геотермические градиенты мерзлых и талых горных пород равны Г1 = 0,02584 °С/м и Г2 = 0,02911 °С/м соответственно. При задании граничного условия II-го рода получим Г1 = 0,02550 °С/м, Г2 = 0,02829 °С/м, температура на нижней границе составит TL = 80,6 °C.
Далее была использована геокриологическая летопись неоплейстоцена Сибири (эпоха Брюнес – последние 800 тыс. лет) [4], которая составлена на основе интерпретации байкальской записи биогенного кремнезема [10]. С.М. Фотиев в своей работе [4] рекомендует рассматривать эту летопись в качестве опорной для всех континентальных районов России. Для интерпретации геокриологической летописи голоцена (11 тыс. лет назад – до современности) используем данные из кривой изменения относительной температуры на границе льдообразования в Антарктиде [6].
Следующий вычислительный эксперимент выполнен с использованием характеристик климата и геокриологических условий хронов неоплейстоцена (табл. 1–3 из [4]) и голоцена [6] при tmax = 800 тыс. лет и ранее найденных значениях H, T0, Г1, Г2. Получено, что при задании граничного условия I-го рода на нижней границе за 600 тыс. лет криогенная толща увеличилась до 431 м, а при задании граничного условия II-го рода – до 477 м. За это время температура поверхности уменьшилась до T0 = –11,47 °С. Тогда при задании граничного условия I-го рода геотермические градиенты мерзлых и талых горных пород будут равны Г1 = 0,02661 °С/м и Г2 = 0,03339 °С/м соответственно. При задании граничного условия II-го рода получим Г1 = 0,02405 °С/м, Г2 = 0,02925 °С/м, температура на нижней границе TL = 73,8 °С.
На рис. 1 представлена динамика мощности многолетнемерзлой толщи для Средневилюйской площади. Результаты расчетов показали, что нижняя граница многолетнемерзлой толщи за период 200 тыс. лет могла варьироваться от 100 до 140 м. При t = 0 и условии постоянного значения температуры на глубине 3000 м мощность мерзлоты достигает 450 м в случае изменения климата по байкальской летописи, 477 м – по Антарктиде. При условии постоянного значения геотермического теплового потока на этой же глубине мощность мерзлоты в случае изменения климата по байкальской летописи и по Антарктиде составляют соответственно 516 м и 560 м. Эти значения соответствуют геотермическим наблюдениям современности, так как по результатам геотермических исследований для данной площади глубина залегания нулевой изотермы колеблется от 480 до 630 м [9, 11, 12].
Рис. 1. Динамика подошвы многолетнемерзлой толщи для Средневилюйской площади при различных сценариях изменения климата: 1, 2 – по байкальской летописи; 3, 4 – по Антарктиде; сплошные кривые – при условии постоянства температуры горных пород (83 °С) на нижней границе (3000 м), точечные кривые – при условии постоянства внутриземного теплового потока (61 мВт/м2) на нижней границе
Рис. 2. Распределение температуры горных пород для Средневилюйской площади при различных сценариях изменения климата: пунктир 1 – начальное (200 тыс. лет назад); 2, 4 – по байкальской летописи; 3, 5 – по Антарктиде; сплошные кривые 2 и 3 – при условии постоянства температуры горных пород (83 °С) на нижней границе (3000 м), точечные кривые 4 и 5 – при условии постоянства внутриземного теплового потока (61 мВт/м2) на нижней границе; 6 – геотермические данные скважины № 3 (1971 г.)
На рис. 2 приведены начальная и конечные (через 200 тыс. лет) термограммы, рассчитанные для описанной выше модели изменения температуры земной поверхности по байкальской летописи и по Антарктиде. Видно, что термограмма, полученная при задании граничного условия I-го рода, в области талых горных пород практически совпадает с современной термограммой (ср. кривые 2 и 3 с кривой 6) до подошвы мерзлоты.
Заключение
Расчеты показали, что положение нижней границы многолетнемерзлой толщи зависит и от выбранного сценария изменения палеоклимата, и от условий задания величины внутриземного теплового потока. Первый фактор смещает ее глубину примерно на 20–40 м, второй – на 70–90 м, однако во всех случаях вычисленные значения практически соответствуют современным результатам геотермических измерений на данной площади.
Результаты таких вычислений позволят моделировать и прогнозировать состояние природной среды, в том числе развитие экзогенных процессов, контролировать устойчивость инженерных сооружений, обеспечивая обоснованные проектные решения, в области развития нестационарных многолетнемерзлых толщ и в первую очередь для региона Вилюйской синеклизы.
Работа выполнена в рамках Комплексной программы фундаментальных исследований СО РАН № II.1 «Междисциплинарные интеграционные исследования» на 2018–2020 гг. (проект – Палеореконструкция теплового поля и криолитозоны Вилюйской синеклизы в позднем плейстоцене-голоцене).