Коллектив авторов - Строение и история развития литосферы

Для решения задачи о распределении температур в разрезе используется численный метод конечных элементов с квадратичной аппроксимацией функции температуры между узлами прямоугольной сетки. В программе предусматривается сетка 41х41 узел (т. е. решается двумерная задача), линейные размеры узла по осям X и Z возможно изменять по требованию оператора. Внутри области моделирования с помощью программы ввода и оцифровки данных задается конфигурация контрастных сред и их теплофизические свойства: температуропроводность а ( м 2/с ), теплопроводность k ( Вт/(мК) ) и нормированная плотность тепловых источников Q/c∙ ρ ( К/с ). В расчетной части комплекса задаются линейные размеры области моделирования ( Lx и Lz , в км), которые определяют линейные размеры узла ( Lx/41 и Lz/41 ), а также временной интервал дискретизации решения (в млн лет). Временной шаг итерационного процесса автоматически выбирается программой и рассчитывается как τ = 10 -7 ( Z 2 /4a ), где Z – толщина области моделирования.

В результате численного решения уравнения теплопроводности:

,(1)

где k, c, ρ – соответственно, теплопроводность, теплоемкость и плотность слоев литосферы, А ( х, z ) – плотность источников тепла в слое τ, – время, мы получаем распределение температур и тепловых потоков q ( z ) и q ( x ) для принятой теплофизической среды в конечный момент временного этапа дискретизации. Полученный файл результатов переименовывается в файл начальных температур и на следующем этапе начинается расчет с конечного момента предыдущего этапа. Возможность дискретизации решения удобна, если есть необходимость изменить теплофизическую среду в связи со структурно-вещественными перестройками геологического разреза, задать распределение новых источников и стоков тепла, а также просмотреть результаты расчета палеотемпературного поля. Если в модели требуется задать тепломассоперенос, то его можно имитировать установкой граничных температур и/или адиабатического градиента в интервале глубин, на которые распространяется адвекция.

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

q гр. =q набл . – q р .

Последний рассчитывается на основании сейсмической информации о мощности слоя ( z i ) и его составе, а также из общепринятых традиционных величин удельной теплогенерации ( А(x,z) ) для соответствующего типа пород: ( q р .) I =А(х, z)z i ( Смыслов и др.,1979 ).

Точность расчетов оценивалась по двум критериям: во-первых, по совпадению модельного и измеренного в скважинах теплового потока; во-вторых, по совпадению температур на пересечении профилей.

Построение трехмерной региональной геотермической модели производится с помощью пакета трехмерной графики «TECPLOT v.9.0» (Amtec Engineering Inc.), который позволяет провести объемную интерполяцию наблюденного поля (в нашем случае температуры, теплового потока, а также структурных сейсмотомографических границ) в координатах: «широта-долгота-глубина». Для подготовки файлов данных в формате «TECPLOT v.9.0» нами была написана специальная программа, которая при задании координат начала и конца профиля, а также интервала разбиения по глубине производит трансформацию текстового файла, содержащего результаты термического моделирования, в формат базы данных «TECPLOT». Программа «TECPLOT» предусматривает трехмерную интерполяцию по сетке любой конфигурации. Во всех случаях нами используется неравномерная сетка, «привязанная» к простиранию сейсмических профилей, вдоль которых выполнялись двумерные расчеты глубинных температур. Таким образом, уравнение (1) дополняется третьей недостающей производной и «превращается» в трехмерное уравнение теплопроводности. Результаты расчета температурного поля с помощью интерполяционного приближения будут отличаться от результатов, получаемых с помощью решения прямой задачи для трехмерной среды, только на величину искажений, обусловленных пертурбацией теплового потока вдоль оси Y. Однако многолетний опыт двухмерного термического моделирования показывает, что эта величина ничтожно мала. Например, для типичного осадочного бассейна тепловой поток вдоль оси Y составляет не более 0,04 мВт/м 2, т. е. на три порядка ниже фоновой величины глубинного теплового потока. Следовательно, получение трехмерного распределения температур с помощью объемной интерполяции практически адекватно решению прямой задачи для объемного теплового поля.