Scientific journal
Advances in current natural sciences
ISSN 1681-7494
"Перечень" ВАК
ИФ РИНЦ = 0,791

MODELING GEOTEMPERATURE FIELD AND THICKNESS OF CRYOGENIC ROCKS IN THE VILYUY SYNECLYSE

Rozhin I.I. 1 Argunova K.K. 1
1 Institute of Oil and Gas Problems SB RAS – a separate subdivision of the FGBUN Federal Research Center «Yakut Scientific Center SB RAS»
The work is devoted to numerical modeling of the dynamics of the temperature field and permafrost for various climate scenarios and heat transfer conditions in relation to different areas of the Vilyuy syneclise. It is based on the mathematical model of the two-phase Stefan problem, and in the computational algorithm for solving it – on the Samarsky-Moiseenko method. Using this model, the temperature field and dynamics of the cryogenic stratum were calculated for a given temperature history in specific geological conditions. Climate changes are realized by setting the following scenarios at the upper boundary of the rock massif: the relative temperature change curve at the ice formation boundary in Antarctica and geocryological interpretations of the paleoclimatic records based on comprehensive studies of the content of diatom valves and biogenic silica in the bottom sediments of Lake Baikal. At the lower boundary of the array, two conditions are set: the constancy of temperature or the constancy of the underground heat flux. The results of a computational experiment with averaged values of the thermophysical properties of sedimentary rocks corresponding to the Srednevilyuyskoe field are presented. Calculations showed that the position of the lower boundary of the permafrost layer depends on the selected scenario of paleoclimate change and on the conditions for setting the value of the intraterrestrial heat flux. The first factor shifts its depth by about 20-40 m, the second – by 70-90 m, however, in all cases, the calculated values practically correspond to modern results of geothermal measurements on this area. The results of such calculations will allow us to simulate and predict the state of the environment for other areas, including development of exogenous processes, control the stability of engineering structures, providing sound design decisions in the field of the evolution of non-stationary cryogenic rocks, and especially for the Vilyuy syneclise region.
geotemperature field
permafrost
paleoclimate reconstruction
Stefan problem
numerical modeling

Вилюйская синеклиза по высокой концентрации углеводородных природных ресурсов и их экономической значимости принадлежит к числу важнейших регионов России. В настоящее время в пределах Вилюйской синеклизы открыты 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 °С/м. На верхней границе задается линейная функция температурной истории rogin01.wmf, где (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].

rogin1.tif

Рис. 1. Динамика подошвы многолетнемерзлой толщи для Средневилюйской площади при различных сценариях изменения климата: 1, 2 – по байкальской летописи; 3, 4 – по Антарктиде; сплошные кривые – при условии постоянства температуры горных пород (83 °С) на нижней границе (3000 м), точечные кривые – при условии постоянства внутриземного теплового потока (61 мВт/м2) на нижней границе

rogin2.tif

Рис. 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 гг. (проект – Палеореконструкция теплового поля и криолитозоны Вилюйской синеклизы в позднем плейстоцене-голоцене).