Численная реконструкция термической истории осадочных бассейнов

счета на элемент dZ каждый промежуток времени dt, консолидации среды и

теплопереноса в ней

Cp – теплоемкости на единицу объёма и массы породы, ρ

- плотность, A - генерация тепла в единице обьема, Z - глубина и t - время.
Скорости конвективных движений V и само решение уравнения рассматриваются в системе координат, связанной с “неподвижным” фундаментом,

сохранения энергии для прямоугольного параллелепипеда ΔZ×ΔS.
Уравнение существенно нелинейно, так

как все его коэффициенты изменяются с изменением пористости (глубины), температуры и времени. Область решения наращивается при отложении осадков и уменьшается при эрозии.

(2а)

его основании z=z2 поддерживается температура Т1 и Т2, соответственно, и

Vw – скорость движения грунтовых вод (Vw > 0 для нисходящего движения), то стационарное распределение температур имеет вид (Огильви, 1959):



где ξ=(z - z1) / (z2 - z1) и число Пекле, Pe, (или число Нуссельта, Nu)

Число Пекле характеризует теплоперенос в конвектирующей среде и равно отношению конвективной составляющей тепло-вого потока с перепадом температуры ΔТ=Т2 – Т1 к кондук-тивной (Bredehoeft and Papadopulos, 1965).

0 - восходящее течение повышает температуру пород слоя над линейным

распределением. q > q0 в верхней половине слоя и q < q0 в его нижней половине. Pe > 0 - нисходящее течение - ситуация обратная.

где .

При Pe<<1 получаем:

течения грунтовых вод имеют пренебрежимое влияние на температурный режим осадков.

Поэтому для скоростей осадконакопления, типичных в осадочных бассейнах конвективное движение выжатых поровых вод не оказывает заметного влияния на термический режим бассейнов.

Исключение составляют бассейны с лавинным осадконакоплением (V > 0.5 км/млн.лет), например, в Южно-Каспийский бассейе.

ростом солености вод с глубиной (Kukkonen, 1995).
Другой пример: бассейн

Уэд эль-Миа, где скорости осадконакопления не превосходили 145 м/млн.лет, а скорости восходящего движения выжатых - 3 мм/год и число Пекле Pe < 0.05.

сотни км практически не сказывается на вертикальном (глубинном) профиле температур

при условии достаточного удаления изучаемых областей от мест инфильтрации метеорных вод на поверхности бассейна.
Причина – вода успевает достичь теплового равновесия с вмещающими породами на расстояниях x >> [Vw⋅(Δh)2/κrock], где Vw – скорость воды в водоносном горизонте, Δh – мощность горизонта и κrock – термическая диффузия пород. Для Vw=10 м/год, Δh=100 м и κrock=3×10-7 м2/сек термическое равновесие водоносного горизонта с вмещающими породами достигается уже на расстояниях 300 – 500 м от входного сечения горизонта.

изменение коэффициентов с глубиной и временем, определяемое зависимостью теплофизических параметров

от литологического состава пород, их пористости ϕ и температуры Т


Плотность: ρs(Z) = ρm· [1 - ϕ (Z)] + ρw ·ϕ (Z)
Теплоёмкость: Cvs(Z) = Cvm ·[1 - ϕ (Z)] + Cvw ·ϕ (Z)
Теплопроводность: Ks(Z) = Km (1 - ϕ (Z)) ·Kwϕ(Z)
Теплогенерация: As(z)=Am⋅[1-ϕ(z)]

этих пород и “w” – поровые воды,
ϕ - пористость

пород, Z – глубина.
Теплопроводность минерального скелета уменьшается с ростом температуры:
Km = Km(T=0°C) / (1+Al⋅T)
(Km(галит) падает в 1.5 раза при росте Т от 0° до 100°С)

Термофизические параметры осадочных пород

ρw = ρw0⋅[1 – 0.000317⋅(T-T0) – 0.00000256⋅(T-T0)2]

ρw0 = 998.2 и

1030 кг/м3 - плотности пресной и

морской воды при Т= 10°С

Ттеплопроводность

Кw = 0.565 + 0.00188⋅T - 0.00000723⋅T2 (0 < T < 137°C)

Kw = 0.602 + 0.00131⋅T - 0.00000514⋅T2 (137 < T < 300°C)

Kw в Вт/м⋅°К и Т в °С.

зависимости от пористости, температуры и литологии пород
Рис. K(z) и A(z)

в разрезе скв. Ахмерова (Западный Башкортостан) (Пунктир - пределы изменения измеренных величин).
Породы верхнего венда - максимум песчаной фракции и максимум Кm.
.Для z < 6 км – эффект пористости на К, для z > 6 км – эффект температу-
ры на Кm. Породы среднего рифея – максимум глин и максимальное А.

температуры dT/dz = 30°C/км; глины: пористость - ϕ=0.70⋅exp(-z/1.8(км)), Km(Т=0°С)=2.09 Вт/м°К,

α=0.0005°С-1; песчаник: ϕ=0.40⋅exp(-z/3.0(км)), Km(Т=0°С)=5.44 Вт/м°К, α=0.0030°С-1; теплопроводность воды бралась постоянной Kw=0.565 Вт/м°К; теплопроводности пород вычислялись по формуле (5-7) и (5-11), пористость и матричная теплопроводность смеси пород вычислялись по формулам (4-10) и (5-12).

пород - ключевой параметр в моделировании бассейнов.

Корреляционная связь теплопроводности со

значениями пористости, ϕ, плотности, ρ, и скоростью продольных волн Vp (Blackwell and Steele, 1989):

для пород Калифорнийской долины (коэффициент корреляции 0.926):
K = 0.84 – 0.040×ϕ + 0.000695×Vp
и для пород Северного моря (коэффициент корреляции 0.9):
K = 3.60×ρ - 0.049×ϕ - 0.160×Vp –5.50
Выше ϕ - пористость в %, Vp – в м/сек, ρ - плотность в г/см3 и К в Вт/м°К. Последнее соотношение проверено на измерениях теплопроводности образцов ряда пород (Blackwell and Steele, 1989): алевролитов и аргиллитов, сланцев, песка, песчаника, кварцита, литофицированных песков, грауваков, известняков, доломитов, солей, ангидритов, углей, воды, гранита, базальта и андезита, риолитового стекла, риолитового пепла и риолитового консолидированного туфа.

осадочных породах, может вносить заметный вклад в формирование температурного режима

бассейна

Её значения сильно зависят от свойств осадочныx пород как адсорбентов, и соответственно от настоящего и прошлого гидрогеологического режима осадочной толщи. По этой причине отсутствует прямая связь генерации тепла с современной пористостью осадков. Предполагая, что радио-активные элементы концентрируются в основном минеральной частью породы (её скелетом), следует ожидать роста объёмной генерации тепла по мере сокращения размеров пор по закону, близкому к:
A(Z) = A(0) * [ 1 - ϕ(Z) ] / [ 1 - ϕ(0) ]
Тогда генерация радиогенного тепла для глин возрастает в 2.5 раза при погружении породы на 4 - 6 км. Движение флюидов через матрицу пород, в том числе и поровыx вод, может нарушать эту зависимость, становясь важным фактором для накопления радиоактивныx элементов в объёме осадков. Не только гидрологический режим, но и состав пород и флюидов в значительной степени влияют на этот процесс.

9.79 Сu + 2.64 СТh) 10-5
ρ (кг/м3), CU

и CTh (1 ppm=10-6 г/г); CK (wt% = г К/ г породы)%

покрова и может достигать 30% от теплового потока на поверхности

в бассейнах с мощным осадочным заполнением (Баренцево море (Вержбицкий, 2002) и Мексиканский залив (McKenna et al., 1998)).

Материнские породы имеют высокое содержание U, Th, K Cилурийские глины (Алжир; ≈70 м): A≈6.7 мкВт/м3 (ΔТ≈0.45°С) Баженовские глины (Уренгой; ≈ 15-50 м): A≈7 мкВт/м3 (ΔТ≤0.5°С). Т.е. тепловой эффект выделения радиогенного тепла в узких слоях материнских пород (≤100-150 м) не имеет заметного влияния на температуру и зрелость этих пород.

В моделировании полагают, что по мере погружения
осадков генерация тепла в них увеличивается по закону: A(z)=Am⋅[1-ϕ(z)]
(см. слайд 10)

в соот-
ветствии с глобальным цик-
лом эволюции литосферы
Земли.
Более того:
разные площади одного

бас-
сейна могут отвечать раз-
ным типам литосферы (пас-
сивные окраины, слайд 17).
Тип фундамента бассейна
может меняться также и в
процессе его развития.

океаническому
Бассейн Сантос, Бразилия

4 типа литосферы : 1) Континентальная; 2) Океаническая; 3) Задугового

(краевого) моря и 4) с произвольными толщинами и характеристиками слоёв.
Вклад радиогенной составляющей пород литосферы в поверхностный тепловой поток увеличивается от Qрад = 10.5 мвт/м2 для океанической и Qрад = 13.-15. мвт/м2 для литосферы окраинных морей до Qрад = 23-25.мвт/м2 для континентальной литосферы.
Но различие в величине Qрад в различныx тектоническиx провинцияx может быть значительнымв

от Qрад = 10.5 мвт/м2 для океанической и Qрад =

13.-15. мвт/м2 для литосферы окраинных морей до Qрад = 23-25.мвт/м2 для континентальной литосферы.

Такие значения генерации тепла Qрад означают, что, например, для потока тепла на поверхности 40 - 50 мВт/м2 температура пород океанической мантии будет расти с глубиной в полтора - два раза быстрее, чем в континентальной мантии.

1999; McKenzie et al.,2005) :
С ростом температуры увели-чивается плотность фононов,

но падает длина их свободного
пробега, что ведёт к падению теплопроводности. Но одновре-менно с ростом Т растёт интен-сивность радиационного пере-носа тепла фотонами.

Здесь Т в °С, b = 5.3, c = 0.0015, d0 = 1.753×10-2, d1 = - 1.0365×10-4,
d2 = 2.245×10-7, d3 = - 3.4071×10-11. Вычисленные k(T) показаны на рисунке сплошной линией.

1972), заметно превышает измеренные значения при температурах 900 – 1000ºС

и выше.

Вклад радиационной составляющей в теплопроводность мантийных пород оказался заметно меньше предполагаемого в (Schatz and Simmons, 1972) из теоретических соображений.

смесь 0.89 форстерита и 0.11 фаялита) с температурой, Ср (Т),

определяются согласно (Berman, Aranovich, 1996) и (McKenzie et al.,2005) по формуле:

Ср = 1618.27 – 12538.9 / SQRT(T°K) – 173.7×107 / T(°K)3 (Дж/кг°K)

Они согласуются с литературными данными для диопсида (Cawthorn, 1975; Robie et al,1978; Bottinga, Steinmetz, 1979). Изменение Ср от Т согласно этой формуле:
 
Т (°К) = 298.15 600 1000 1600
Т (°С) = 25 226.85 726.85 1326.85
Ср (Дж/кг°K) = 826.55 1098.3 1220.0 1304.4

или в периоды его тепловой реактивизации характеризуются заметными изменениями толщины

литосферы в ответ на плавление пород в её основании или на затвердевание пород прилегающей астеносферы. Выделение или поглощение скрытой теплоты плавления, заметно увеличивают «тепловую инерционность литосферы бассейна», замедляя процессы её остывания или разогревания.
Энтальпийное приближение: (Карслоу, Егер, 1964).
Доля расплавленной фракции, f, в породе линейно увеличивается с ростом температуры, T, в пределах интервала температур между солидусом и ликвидусом породы:

при Ts < T < Tl

В этом интервале температур в уравнении теплопроводности в члене c ( ∂/∂t) теплоёмкость Ср заменялась на Cp‘ :

перидотитовых пород), Ts и Tl - температура солидуса и ликвидуса

пород, соответственно. Кривая Тs для перидотита ( < 0.2% H2O; Wyllie, 1979) имеет вид:
Ts =1.03 + 0.00391×P (0 ≤ P ≤ 12.8)
Ts = 1.236 - 0.0122P × (12.8 ≤ P ≤ 20.16)
Ts = 0.742 + 0.0122 × P (P ≤ 20.16)
Ts в 1000°C и давление P в Кбар. Кривая солидуса Ts(P) определяет глубины подошвы литосферы.
Температура ликвидуса Тl определяется как:
Tl = Ts + 450 - 600oC для перидотитовых и
Tl = Ts + 75 - 150oC для базальтовых пород.

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

(McKenzie, Bickle, 1988).
2 и 7 – TL и Ts

влажного перидотита с 0.5% Н2О
(Wyllie, 1979).
4, 6 и 8 - солидус для сухого оливина (4), оливина с 810 H/106Si (=0.00547% (wt H2O) – типичное содержаие летучих в породах мантии на глубинах 90 – 120 км, генерирующих базальт под осевыми зонами СОХ) (6)) и влажного оливина (8) (Hirth&Kohlstedt, 1996).

5 – Ts перидотита океанической литосферы (Asimov et al.,2001).
9 и 10 – Ts влажного и сухого перидотита мантии под осевыми зонами СОХ (Braun et al.,2000),.
11 - солидус нормального сухого перидотита (Hirschmann M.M., 2000).

пород мантии Ts(z) с текущей геотермой литосферы Т(z,t).

(Z=0) определяется палеоклиматическими условиями на поверхности бассейна. Она может быть

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

Европе, Азии, Западной Африке, Гренландии,
Австралии.

ПАЛЕОЗОЙСКАЯ ЛЕДНИКОВАЯ ЭРА (460

– 230 млн.лнт)

Поздний ордовик – ранний силур (460 – 420 млн.л. назад)\\
Восточный Ньюфаундленд, север Аравийского п-ва, Централь-
ная Сахара.
Поздний девон – ранний карбон (370 – 355 млн.л. назад)
Бразилия, Африка (?)
Каменоуглльно-пермский период (360 – 230 млн.л. назад)
Пик - (310 – 270 млн.л. назад)
Гондвана – площадь оледенения в 2-3 раза выше площади
современной Антарктиды

КАЙНОЗОЙСКАЯ ЛЕДНИКОВАЯ ЭРА (30 мл назад – 0)

затухать с глубиной. Так для решения T(0,t)=T0+ΔT⋅cos(ω⋅t) на поверxности (Карслоу

Егер 1964)




Для К=0,0082 кал/см°Cсек, Cр=0,25 кал/г°C, ρ=3,3 г/см3, т.е. χ =0,01 см2/сек (фундамент) xарактерная глубина затуxания поверxностныx колебаний температуры будет около 0,17 м для суточныx 2.25 м - для годовыx и около 3.2 км для колебаний с периодом 1 млн.лет.

Карта геоизотерм «нейтрального слоя» Прикаспийской впадины и обрамляющих её территорий (Котровский, 1986).
1 – станции, где определялись температуры «нейтрального слоя», 2 – температуры, °С.
Глубина нейтрального слоя - от 15 до 25 м.

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

с учётом уменьшения температуры воды с глубиной. Для многих современных

открытых океанов температуры воды убывает со скоростью около 4°C/100 м в верхних 200 м воды и около 2°C/100 м в более глубокиx слояx. Минимальная температура воды на больших глубинах около 1-2°C. В Полярных морях: на поверхности Т(воды (льда)) = Т (воздуха) < 0°C. Т растёт до 0°C до глубины 2 м, затем слегка увеличивается чуть выше 0оС, затем уменьшается, второй раз переходит через 0оС при z=7-8 м и понижается до –1.7 - -1.8оС при z=35-40 м. Вода остается при такой температуре до глубин 150 – 200 м и не подвержена здесь сезонным колебаниям температуры. При z > 150-200 м температура растёт и достигает –0.8оС на глубине z = 500 м.

системе моделирования бассей-нов ГАЛО включает в себя осадочную толщу бассейна

и его фундамент. Последний, кроме литосферы, может включать и верхнюю часть астеносферы.
Основное требование – достаточ-ная глубина нижней границы, слабое влияние её выбора на результаты расчётов.

Нижняя граница области счёта – это и уровень изостазии, поэтому она располагается в пределах реологически слабого слоя мантии, залегая либо глубже подошвы литосферы, либо в её нижних слоях. Это предполагает возможность перемещения вещества мантии под действием минимальныx разностей напряжений, что Такое движение будет способствовать выравниванию температуры в горизонтах мантии на глубинах порядка ZM и установлению здесь термического режима со сравнительно слабым изменением температуры, характеризующимся градиентами близкими к адиабатическим (около 0,3°С/км).

Qmin)
Для определения ZМ и ТМ расcчитывается распределение температуры T(z,Q) с

глубиной в фундаменте для предполагаемого минимального теплового потока в истории бассейна на его поверхности Q = Qmin и затем пересечение этой геотермы с кривой солидуса пород мантии. Если такого пересечения нет, то ZM=200 км и TM = T(ZM). (ZМ ≈ 90 км, TМ ≈ 1050°С в сл. 11)
Тепловая активизация литосферы бассейна численно воспроизводится подъёмом одной из изотерм (1100, 1050, 1000, 950 или 900°С < TM), ближайшей к ТМ, до определённого, предполагаемого в модели уровня (слайд 11).
В термически активныx районах (рифтогенез или тепловая реактивизация) поддержание на границе высокого теплового потока приводит к завышению температурных градиентов в мантии и появлению здесь аномально высоких температур, достигающих 1600-1800°C. Такие значения противоречат геофизическим оценкам верхнего предела 1400-1500°C для температур мантийных пород на глубинах от 100 до 200 км (Anderson 1979, 1980).
В пакетах MATOIL, GENEX, TEMISPACK, PDI тепловой поток задаётся в виде ступенеобразной функции геологического времени. Высоты ступеней подбираются из условия совпадения вычисленных и наблюдённых значений отражательной способности витринита, Ro%.

температур - ввод произвольного кусочно-линейного распределения температуры с глубиной (тепловая

релаксация надвиговыx структур или бассейнов с межпластовыми интрузиями в осадочной толще или фундаменте).
Для большинства осадочныx бассейнов начальные условия задаются через значение начального теплового потока на поверхности бассейна Q0 (или градиента температуры), характеризующего термическое состояние бассейна на стартовом этапе его развития.
При заданном Qo и известном распределении термофизических характеристик литосферы бассейна с глубиной начальное распределение температур определялось решением стационарного уравнения (5) T(z, t=0) = T(z, Q0).

To и Qo – температура и тепловой поток на повер-хности бассейна в момент его зарождения (t=0).

ТМ (в области Ts < T < TМ) предполагается линейное,

квазиадиабатическое распределение Т(z).
Эволюционный ряд осадочных бассейнов даёт возможность сопоставить каждому из бассейнов тектоническую обстановку его зарождения. При выборе значений начального теплового потока Qo (или градиента температур) они принимаются близкими к современным потокам, наблюдаемым в районаx с тектонической обстановкой аналогичной или близкой к обстановке в период зарождения рассматриваемого бассейна.

Qo ≅ 100-110 мвт/м2 - осевые зоны континентальных рифтов (озеро Байкал, Африканская рифтовая система)
Qo ~ 65 - 85 мвт/м2 - плечи рифтовых структур (Смирнов, 1980).
- бассейны на океанической коре
где Qo в мВт/м2 и t - возраст коры в
млн.лет. Она может использоваться как
при анализе классических океанических бассейнов, так и "океаническиx" бассейнов - останцев, возникающиx на нерегулярностяx континентальныx границ при закрытии океанов.

на следующем, (n+1)-ом, шаге времени.
Неявная конечно-разностная схема уравнения теплопроводности

значений 0.5 - 10 м на поверхности осадочной толщи до

1000 -3000 м в основании области счета на глубине 100-200 км. Корректное численное воспроизведение термической истории бассейна требовало достаточно мелкого разбиения осадочной толщи по глубине на шаги Δz, что приводило к увеличению числа узлов в процессе счёта и одновременно времени счёта. Чтобы держать этот процесс в разумных границах, в программе предусмотрена процедура переписывания (нового разбиения) области счёта на шаги по глубине всякий раз как число шагов превосходило 800.
Размер шагов Δz рос линейно с глубиной в пределах каждой литологически однородной свиты и оставался непрерывной функцией во всей области счёта.
Число шагов по времени достигало 3000, а размер шагов менялся от 1 000 до 1 000 000 лет.

конвективныx членов неявная конечно-разностная сxема слайда 31 уравнения теплопроводности обладает

безусловной устойчивостью.
Зависимость теплофизических параметров от глубины и времени как функций пористости, литологии пород, истории их погружения, выделение и поглощение тепла плавления, движение выжатых вод, консолидация осадков и т.д. - всё это затрудняет явную оценку сходимости разностной схемы.
Корректность применяемой разностной схемы проверяется тогда сравнением численных решений с имеющимися аналитическими и полуаналитическими аналогами, а также сопоставлением решений, полученных с использованием различных шагов Δt и Δz.

используемые для проверки корректности разностной схемы

1) Отложение несжимаемых однородных осадков

на однородный фундамент с идентичными свойствами

- функция вероятности,

без-

размерное время, κ= K/ ρ⋅Cp - термическая диффузия, ρ - плотность, К - теплопроводность, Ср - теплоёмкость осадков и фундамента, t - время, V - скорость осадконакопления, q∞ - невозмущённый тепловой поток на большой глубине, - безразмерная глубина, z - глубина,

Z1 = (Z - P)/2 , Z2 = (Z + P)/2, - температура на

глубине тепловой волны √κt, соответствующая тепловому потоку q∞.

схеме слайда 31. Рассматривалось отложение осадков со скоростью V =

1 км /млн. лет в течении 20 млн. лет на фундамент со следующими термическими характеристиками: К = 0,005 кал/см сек°C, Cр = 0,25 кал/г°С, κ = 0,008 см2/сек. В процессе численного моделирования основание области счёта принималось на глубине 200 км и здесь поддерживалась температура 1200°C. При воспроизведении процесса с отложением на поверхности порции осадков Δz=20 м каждые 20 тысяч лет значения теплового потока на поверхности и распределения температур на глубине, рассчитанные по разностной схеме, отличались от вычисленных по аналитической формуле не более чем на 0,1% во всем интервале времени 0 ≤ t ≤ 20 млн. лет и глубин 0 ≤ z ≤ 200 км.

Z2, Ф(y) те же, что на слайде 34. Сравнение проведено

для эрозии с постоянными скоростями V = 500 м / млн. лет (t ≤ 20 млн. лет) и V = 1000 м / млн. лет (t ≤ 40 млн. лет) и значениями термической диффузии κ= 0.008 см2/сек и 0,0095 см2/сек, соответственно. При ZM=200 км T=1200°C. Расчёт по разностной схеме с шагами эрозии Δz=10 м за 20 и 10 тысяч лет, соответственно, воспроизводит аналитические значения Т и q c точностью не менее 0,1 %.

q1 - аналитические решение;
q2 - по разностной схеме с Δz = 10 м;
q3 - то же с шагом эрозии Δz = 50 м..

Эрозия 10 км/ 20 млн. лет при фоновом потоке q∞ = 0.3 ЕТП

с V=const на фундамент с другими характеристиками пород
5) Периодически меняющаяся

температура на поверхности однородного полупространства
Во всех случаях адекватный подбор шагов Δz и Δt позволялобеспечить совпадение вычисленных температур и тепловыx потоков с результатами аналитическиx и полуаналитическиx решений с точностью до десятыx долей процента