Вдоль каждого из профилей выполнялся расчет глубинных температур с помощью программного пакета «TERMGRAF», разработанного нами (Хуторской, 1996).
Для решения задачи о распределении температур в разрезе используется численный метод конечных элементов с квадратичной аппроксимацией функции температуры между узлами прямоугольной сетки. В программе предусматривается сетка 41х41 узел (т. е. решается двумерная задача), линейные размеры узла по осям X и Z возможно изменять по требованию оператора. Внутри области моделирования с помощью программы ввода и оцифровки данных задается конфигурация контрастных сред и их теплофизические свойства: температуропроводность а (м2/с), теплопроводность k (Вт/(мК)) и нормированная плотность тепловых источников Q/c∙ρ (К/с). В расчетной части комплекса задаются линейные размеры области моделирования (Lx и Lz, в км), которые определяют линейные размеры узла (Lx/41 и Lz/41), а также временной интервал дискретизации решения (в млн лет). Временной шаг итерационного процесса автоматически выбирается программой и рассчитывается как τ = 10-7(Z2/4a), где Z – толщина области моделирования.
В результате численного решения уравнения теплопроводности:
,(1)
где k, c, ρ – соответственно, теплопроводность, теплоемкость и плотность слоев литосферы, А(х, z) – плотность источников тепла в слое τ, – время, мы получаем распределение температур и тепловых потоков q(z) и q(x) для принятой теплофизической среды в конечный момент временного этапа дискретизации. Полученный файл результатов переименовывается в файл начальных температур и на следующем этапе начинается расчет с конечного момента предыдущего этапа. Возможность дискретизации решения удобна, если есть необходимость изменить теплофизическую среду в связи со структурно-вещественными перестройками геологического разреза, задать распределение новых источников и стоков тепла, а также просмотреть результаты расчета палеотемпературного поля. Если в модели требуется задать тепломассоперенос, то его можно имитировать установкой граничных температур и/или адиабатического градиента в интервале глубин, на которые распространяется адвекция.
Для каждого профиля при моделировании задается краевая температура на верхней границе в соответствии с метеорологическими данными и тепловой поток – на нижней границе (qгр.), соответствующий измеренному в ближайших скважинах значению (qнабл.) за вычетом теплового потока, генерируемого в слое земной коры выше нижней границы области моделирования при спонтанном распаде долгоживущих радиоизотопов (qр.), т. е.
qгр.=qнабл. – qр.
Последний рассчитывается на основании сейсмической информации о мощности слоя (zi) и его составе, а также из общепринятых традиционных величин удельной теплогенерации (А(x,z)) для соответствующего типа пород: (qр.)I=А(х, z)zi (Смыслов и др.,1979).
Точность расчетов оценивалась по двум критериям: во-первых, по совпадению модельного и измеренного в скважинах теплового потока; во-вторых, по совпадению температур на пересечении профилей.
Построение трехмерной региональной геотермической модели производится с помощью пакета трехмерной графики «TECPLOT v.9.0» (Amtec Engineering Inc.), который позволяет провести объемную интерполяцию наблюденного поля (в нашем случае температуры, теплового потока, а также структурных сейсмотомографических границ) в координатах: «широта-долгота-глубина». Для подготовки файлов данных в формате «TECPLOT v.9.0» нами была написана специальная программа, которая при задании координат начала и конца профиля, а также интервала разбиения по глубине производит трансформацию текстового файла, содержащего результаты термического моделирования, в формат базы данных «TECPLOT». Программа «TECPLOT» предусматривает трехмерную интерполяцию по сетке любой конфигурации. Во всех случаях нами используется неравномерная сетка, «привязанная» к простиранию сейсмических профилей, вдоль которых выполнялись двумерные расчеты глубинных температур. Таким образом, уравнение (1) дополняется третьей недостающей производной и «превращается» в трехмерное уравнение теплопроводности. Результаты расчета температурного поля с помощью интерполяционного приближения будут отличаться от результатов, получаемых с помощью решения прямой задачи для трехмерной среды, только на величину искажений, обусловленных пертурбацией теплового потока вдоль оси Y. Однако многолетний опыт двухмерного термического моделирования показывает, что эта величина ничтожно мала. Например, для типичного осадочного бассейна тепловой поток вдоль оси Y составляет не более 0,04 мВт/м2, т. е. на три порядка ниже фоновой величины глубинного теплового потока. Следовательно, получение трехмерного распределения температур с помощью объемной интерполяции практически адекватно решению прямой задачи для объемного теплового поля.
2. Термотомография Западно-Арктического бассейна
Нами проведено два «этапа» геотермического моделирования вдоль сети геотраверсов в Баренцевом и Карском морях с целью определения глубинных температур в земной коре, для оценки глубины залегания температурного интервала катагенетического преобразования органического вещества, а также для исследования характерных проявлений геотемпературного поля, контролирующих локализацию известных месторождений газа и газоконденсата.
На первом этапе нами использованы 123 сейсмических геотраверса, на которых строение геологического разреза определено только по геофизическим данным, т. е. на основе интерпретации положения отражающих (МОВ-ОГТ) и преломляющих (ГСЗ) границ. На втором этапе у нас была возможность по тому же региону построить термотомографическую модель по семи геотраверсам, на которых учтены не только сейсмические данные, но и результаты бурения глубоких поисково-разведочных скважин.
Геотермические исследования на шельфе Евразийского бассейна Арктики начались в 70-е годы с измерений погружными зондами ПТГ-3М в Баренцевом море (Методические…, 1983). Однако эти измерения были сделаны на шельфе при глубинах моря до 300 м и глубине внедрения одноканального зонда в осадки максимально на два метра. Полученные результаты продемонстрировали ожидаемое очень мощное влияние экзогенных термических полей на глубинный тепловой поток – сказывалось воздействие сезонных периодических колебаний температуры дна моря за счет инсоляции, а также придонных течений с большим дебитом, приносящих массы воды с контрастной относительно фоновой температурой. Особенно этот фактор заметен в южной и западной частях моря, где проявляется влияние Гольфстрима и его ответвления – Нордкапского теплого течения. В связи с этим зондовые измерения не позволяли получить кондиционные оценки фонового глубинного теплового потока, т. к. диапазон измеренных значений был очень велик – от 0 до более 500 мВт/м2, и не адекватно отражал распределение глубинных термических источников.
Начавшаяся в 80-е годы интенсивная разведка нефтегазовых месторождений на шельфе сопровождалась бурением на акватории и на островах глубоких скважин, в которых проводились каротажные исследования, в том числе и термокаротаж. К этому времени относятся и первые скважинные измерения теплового потока в южной части Карского моря. Обработка термических измерений позволила оценить значения градиентов температуры, а теплофизические исследования керна скважин – теплопроводность пород. Таким образом были получены первые кондиционные измерения теплового потока в регионе (Цыбуля, Левашкевич, 1992), которые после применения методики учета сезонных колебаний температуры были уточнены (Левашкевич, 2005) (рис. 1).
Рис. 1. Геотермическая изученность Баренцевского региона.
Тем не менее, скважинные геотермические измерения для такой обширной территории были редки и не позволяли составить картину распределения температур и тепловых потоков, а тем более рассчитать глубинные температуры в литосфере региона.
Обобщение имеющихся скважинных и зондовых измерений тем не менее позволяет говорить о тенденции повышения теплового потока в северо-восточном и северо-западном направлениях. Так, в зоне сочленения Кольской микроплиты и Балтийского щита среднее значение теплового потока составляет 54 мВт/м2, а в районах Северо-Баренцевской впадины и Центрально-Баренцевского поднятия – 70 мВт/м2. В первом приближении такую тенденцию тренда теплового потока можно объяснить приближением к Северо-Атлантическому центру спрединга, в котором возрастает термическая активность астеносферы. Однако полученные нами в 25-м рейсе НИС «Академик Николай Страхов» данные (подробное описание которых приведено в нашей другой статье) позволяют говорить о новейших геодинамических проявлениях в земной коре Баренцевской плиты. Мы ранее уже высказывали предположение о рифтогенном характере геодинамической активности на основании интерпретации геотермических данных (Хуторской и др., 2000).
