CC BY
42
GEODYNAMICS & TECTONOPHYSICS
PUBLISHED BY THE INSTITUTE OF THE EARTH'S CRUST SIBERIAN BRANCH OF RUSSIAN ACADEMY OF SCIENCES
2011 VOLUME 2 ISSUE 4 PAGES 356-377 DOI:10.5800/GT-2011-2-4-0050
ISSN 2078-502X
New information technology and its capacities for modeling of geosystems
V. I. Gunin
MoGeos Geosystems Simulation Center, Ulan-Ude, pr. 50 let Oktyabrya, 38-18, Russia
Abstract: Numerical solution of the Navier-Stokes equations for free and forced convection in solid and porous media is challenging as an evolution equation for pressure is lacking. An information technology is proposed on the basis of a new system of hydrodynamic equations. It provides for distinguishing between free and forced convection components and allows estimation of parameters of the cumulative convective flow by calculating values of its two components. A classical system of equations of the elliptic-parabolic type is developed for solving the problems of heat and mass transfer. The system describes vortex structures which occur in the gravitational field in all media in case of density stratification. A majority of the available computation methods and schemes can be applied for numerical solutions of the proposed system; therefore modeling can be simplified, while the scope of the system's application can be expanded. Possible applications of the proposed information technology are demonstrated by examples showing how problems of low-mantle plume and diapir formation and rifting in the lithosphere-crust system can be solved.
Key words: information technology, convection, vortex structures, heat and mass transfer, mathematical model, numerical experiment, mantle, lithosphere, crust, mantle plume, diapir, rift.
Recommended by V.A. San'kov, 10 October 2011
Citation: Gunin V.I. New information technology and its capacities for modeling of geosystems // Geodynamics & Tectonophysics. 2011. V. 2. № 4. P. 356-377. doi:10.5800/GT-2011-2-4-0050.
Новая информационная технология и ее возможности при
моделировании геосистем
В. И. Гунин
Центр моделирования геосистем «МоГеос», Улан-Удэ, пр. 50 лет Октября, 38-18, Россия
Аннотация: При численном решении системы уравнений Навье-Стокса в задачах со свободной и вынужденной конвекцией в сплошных и пористых средах возникают сложности по причине отсутствия эволюционного уравнения для давления и из-за ряда других особенностей. Предложена информационная технология, основанная на новой системе уравнений гидродинамики, позволяющая разделить конвекцию на две составляющие - свободную и вынужденную, и, вычислив их значения, получить параметры суммарного конвективного потока. Для решения задач тепломассопе-реноса получена классическая система уравнений эллиптическо-параболического типа, описывающая вихревые структуры, которые формируются в гравитационном поле во всех средах при возникновении плотностной стратификации. При численном решении этой системы можно использовать большинство разработанных расчетных методов и схем, что упрощает моделирование и расширяет область ее применения. На примерах решения задач формирования нижнемантийного плюма, диапира и рифта в системе литосфера-кора показаны возможности этой информационной технологии.
Ключевые слова: информационная технология, конвекция, вихревые структуры, тепломассоперенос, математическая модель, численный эксперимент, мантия, литосфера, кора, плюм, диапир, рифт.
1. Введение
Самым ценным товаром на мировом рынке становятся новые знания, полученные на основе современных информационных технологий (ИТ). Новые знания приобретаются посредством преобразования с помощью передовых ИТ старых знаний (данных) более низкого уровня, накопленных в огромном количестве. Одним из наиболее эффективных методов получения новых знаний является математическое моделирование различных природных систем (геосистем). Полученные таким образом новые знания являются знаниями более высокого уровня, позволяющими понять причинно-следственные связи сложных процессов, явлений и на основе этих понятий вначале научиться их прогнозировать, а затем вводить элементы управления ими.
Основной движущей силой переноса вещества в природе является конвекция, которая, в свою очередь, зависит от сил, вызванных какими-то внешними воздействиями (вынужденная конвекция), и внутренних сил, вызванных плотностными неоднородностями, возникающими за счет градиентов температуры и концентрации вещества (свободная конвекция). Для описания гидродинамики в задачах со свободной и вынужденной конвекцией в сплошных и пористых средах широко используется система уравнений Навье-Стокса, в которой эти силы описываются совместно -одним уравнением. Это приводит к трудностям при ее решении, так как она не является системой типа Коши-Ковалевской, то есть нет эволюционного уравнения для давления. Существуют и другие особенности, возникающие при численном решении данной системы уравнений. Например, наличие малого параметра при второй производной может приводить к существенному изменению гладкости решения, а также вызывать ряд других изменений. К тому же с помощью данной системы уравнений сложно описать вихревые структуры, которые формируются в гравитационном поле во всех средах при возникновении плотностной стратификации.
Для избавления от этих недостатков существует подход, основанный на исключении давления из уравнений за счет применения операции ротор к уравнению для вектора скорости и переходе описания скоростей через функцию тока и вихрь. Такой переход по-
вышает порядок дифференциального уравнения по пространственным переменным, но за счет введения функции тока дает автоматическое удовлетворение уравнения неразрывности, при этом система уравнений усложняется и обычно решается раздельным двухполевым методом. Приравнивание операции ротор скорости к вихрю, которые по определению тождественно равны, не позволяет вычислить функцию тока через другую независимую функцию, что приводит к новым, такого же порядка, трудностям, связанным с заданием граничных условий для вихря [Тару-нин, 1990]. Возникает необходимость в построении новой математической модели, позволяющей избежать указанных затруднений. Для вычисления функции тока и затем определения через нее вектора скорости можно использовать независимую векторную функцию - генератор вихрей. Похожий метод использовался некоторыми исследователями для решения задач в пористых средах, но только для двумерных уравнений [Кирюхин, Сугробов, 1987].
Предлагаемый автором подход основан на расщеплении конвекции на свободную и вынужденную, что позволяет записать новую систему уравнений гидродинамики и в результате ее решения получить значения общего конвективного потока. В работе дано обоснование предлагаемого подхода, а на примерах решения задач формирования нижнемантийного плю-ма, диапира и рифтинга, в системе литосфера - кора, показаны возможности информационной технологии, разработанной на основе новой трехмерной системы уравнений тепломассопереноса.
2. Обоснование нового подхода
И МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Из теории поля известно, что любое векторное поле M можно разложить на две составляющие: M=Mi+M2.
Поле M1=grad U - источниковое или дивергентное поле, где U - произвольная функция, имеющая первую производную, и выполняется равенство rot grad U=0. Функцию U называют скалярным потенциалом векторного поля M. Тогда div M1=div grad U=k15,^ V2U=AU=k15, где k1 - константа пропорциональности, 5 - плотность источников.
Поле M2=rot A - вихревое или циркуляционное поле,
тогда, div rot A=0. Функцию A называют векторным потенциалом векторного поля M. Существует бесконечное число векторных функций A, описывающих одно и то же поле M, такие функции могут отличаться друг от друга с точностью до градиента некоторой скалярной функции ф, так как rot grad ф=0. Это приводит к следующему соотношению: rot rot A=k2W ^ grad div A-V2A=k2W. Так как существует бесконечное множество функций A, описывающих одно и то же поле M, можно выбрать из них такое, которое удовлетворяет равенству div A=0 (калибровка кулона) ^ V2A=AA=k2W, где k2 - константа пропорциональности, W - вектор плотности вихревых линий (линий тока).
Все процессы, протекающие в природе, в той или иной степени связаны с гидродинамическими потоками, вызванными перераспределением вещества в гравитационном поле. Известно, что генератором свободной конвекции является неоднородность плотности среды, а интенсивность конвекции пропорциональна градиенту этой плотности [Гебхарт и др., 1991] Значит, векторную функцию W можно представить в виде градиента скалярной функции, например р (плотность среды), тогда W=grad р. Свободно конвективные потоки вещества, вызванные градиентом плотности, формирующим гидродинамический диполь, удобно описывать с помощью векторной функции тока ¥ - векторный потенциал, а потоки вынужденной конвекции, вызванные источниками, через функцию U - скалярный потенциал, давление или напор.
Тогда перенос среды (вещества) можно записать в виде векторных уравнений, AU=ki5, описывающих источниковое или дивергентное поле, записанное через скалярную функцию U, и A¥=k2Vp, описывающее вихревое или циркуляционное поле, записанное через векторную функцию которое в декартовой системе координат распадается на три скалярных уравнения, инвариантных относительно системы координат:
A¥x=k2 др/дх, A¥y=k2 др/ду, A¥z=k2dp/dz.
(1)
В общем случае векторное поле скоростей, вызванное источниками и генераторами вихрей, можно записать так:
V=Vi+V2=grad U+rot
(2)
Это равенство представляет собой одну из форм теоремы Гельмгольца.
На этом основании, автором предложена новая трехмерная сопряженная математическая модель, позволяющая разделить конвекцию на две составляющие - свободную и вынужденную - и, вычислив их значения, получить параметры суммарного конвективного потока. Данный подход хотя и увеличивает порядок системы уравнений в частных производных, но позволяет более точно вычислить конвективную состав-
ляющую, которая играет основную роль при вычислении полей температуры и концентрации вещества. Добавив к полученной системе уравнений гидродинамики уравнения теплопроводности, массопереноса, уравнение состояния, граничные и начальные условия, получим следующую полную систему уравнений тепло-массопереноса в пористых и вязких средах.
Эта модель в терминах функции тока, давления, температуры и концентрации выглядит так:
(A д^И + д^ ( A д^ 1 + д^ дх [ x dx J ду [ у ду J dz [ z dz
дРф дх
+hit,x, y,z);
дх
д^2 ( A
дРф ду
дх J ду + fг(t,x,У,z) '
д¥2 ( A
ду
8z
д^2 ( A
8z
^ Г Ax ^ | +
дх I 8x J ду
дЩ+д^(a
ду J 8z [ z 8z
дРф ' 8z
+ f зx, У, z);
(3)
(4)
(5)
1+8U (• ■ да
8z ~8z
д T дТ (ах дТ
~д7 _ дх дх
+ дТ ( дТ1 дТ( дТЛ ,,
гу ["у Ту J+& ^ J +Ut,x,у,z);
C+v: C + Vу C + V; C _
8i дх у ду 8z
= C (ядС 1+дС(^+
дх [ дх J ду [ ду
+°ci_ ( scL 1- х, у,
8z [ z 8z J Л ^ ''
-ди. _a .
x ду 8z дх у 8z дх ду
(7)
(8)
V, _
(9)
дУ 2 ду дЦ_; дх ду дz Рф =Ро - Ров(т - Т0 )+ Ров (СХ - о )+Рое(Р - Ро),
С, = Е С,, (10)
где У 12,3 - проекция вектора функции тока на координатные оси; Т - температура; и - потенциальная функция; С1 - концентрация 1-го компонента в системе, 1=1, т, т - количество компонент; t - время; х, у, ъ - координаты, У„, Уу, Уг, Щ, ах, ау, а2, Ах, Ау, Аг, Вх, Ву, Вг, ах, ау, аъ - проекции векторов скоростей конвекции (фильтрации) и тепловой волны, коэффициентов текучести (фильтрации), гидродинамического сопротивления, гидродисперсии, температуропроводности на координатные оси; в, /в, 8, Ро, Рф- термический и концентрационный коэффициенты объемного расширения, коэффициент сжимаемости, исходная и
+
У
у
текущая плотность среды (флюида); T0, C0, T, C - фоновые и текущие температура и концентрация; P0, P -гидростатическое давление, или напор, фоновое и текущее; D - обобщенный коэффициент дисперсии, или коэффициент гидродисперсии; Fi - обобщенная функция кинетики, превращения веществ; f (t, x, y, z) - функция источников или генераторов,
формирующих внешние (сторонние) поля, например транзитный поток или внутренние поля, сформированные в предыдущие моменты времени i=1^5.
Для решения выписанной системы уравнений задаются граничные и начальные условия:
г =Ф; y, z )г=о = F(( У, z);
T\ г = G; t(( y, z =о = F2 (x, y, z);
C г = R; C (x, У, z )t=о = f3 (x, У, z);
P| Г =н; P(, У, z) t=о = F4 (x, y, z); (11)
где Ф, G, R, H, F1-4 - известные функции.
Это система уравнений шестого и более порядка, в зависимости от количества веществ в рассматриваемой системе. Первые четыре уравнения описывают гидродинамику, причем первые три, записанные через функцию тока, - свободную конвекцию, вызываемую внутренними силами (градиентом плотности), а четвертое, записанное через скалярный потенциал, - вынужденную конвекцию, вызываемую любыми внешними силами. Пятое уравнение - теплопроводности, шестое - массопереноса. Прямая связь между этими двумя блоками осуществляется через скорости конвекции (фильтрации) (уравнение (9), обратная - через плотность среды (флюида) (уравнение состояния (10)). Для решения данной системы задаются начальные и граничные условия (11). Любые начальные состояния и внешние воздействия (транзитные потоки и т.д.) или взаимодействия разных веществ между собой и их фазовые переходы задаются с помощью функций источника, записанных в правой части для каждого уравнения.
Представленная модель эллиптическо-парабо-лического типа является классической системой уравнений, для численного решения которой подходят большинство разработанных и используемых расчетных методов и схем. Повышенный порядок системы уравнений практически не влияет на методику и технологию разработки программного кода. Это позволило автору с применением экономичных конечно-разностных схем, основанных на расщеплении уравнений по координатам осей и физическим процессам, используя абсолютно устойчивый локально-одномерный метод решения [Тихонов, Самарский, 1975] разработать пакет программ на языке Fortran, с помощью которого можно решать широкий круг трехмер-
ных задач тепломассопереноса в пористых и вязких средах на персональных компьютерах. При решении задач используются три безразмерных комплекса. Первые два - это тепловое и концентрационное числа Фурье Ро(=ат/Ь2, Рос=Бт/Ь2. В уравнениях, описывающих свободную и вынужденную конвекцию, возникающую за счет градиента плотности среды (флюида) или внешних сил, коэффициент пропорциональности к=(^ц,)х(р/Ъ), определяющий плотность потока линий тока и характеризующий его интенсивность, в общем случае является тензором второго ранга и имеет размерность 1/(м2сек). Физический смысл этого коэффициента можно интерпретировать как текучесть, а его обратное значение - как гидродинамическое сопротивление. Третий безразмерный комплекс Б1с=^/ц,)х(р/Ь)хЬ2т можно назвать критерием плотности распределения линий тока. В уравнениях этих трех комплексов а, Б, ц, g - коэффициенты температуропроводности, диффузии, динамической вязкости, ускорения свободного падения, т, Ь2, Ь - характерные время, площадь и размер в решаемой задаче. Более подробно с методом решения данной системы уравнений можно ознакомиться в работе [Гунин, 2003], в которой сделана и оценка точности решения. Основную нагрузку при интерпретации огромного объема результатов вычислений несет пакет графических программ визуализации, разработанный автором, позволяющий в динамике просматривать значения любой функции, вычисленной при решении поставленной задачи.
Выписанная выше система уравнений универсальна, поэтому ее можно использовать не только для задач гидродинамики. Так как все напряжения возникают за счет неоднородности среды, через перемещения (мгновенные скорости) можно определить поле деформаций, дилатацию, девиатор (тензор напряжений) и решать задачи прочности. Можно использовать первые четыре уравнения для решения задач электромагнитной динамики, через векторный и скалярный потенциалы определять напряженность электрического или магнитного поля, используя их коэффициенты сопротивления или проводимости.
Ниже будут представлены результаты численных решений трех связанных в пространстве и последовательных во времени задач, показывающих возможности новой информационной технологии. Следует подчеркнуть, что рассматриваемые задачи имеют теоретический характер. Они основаны на усредненных данных и параметрах рассматриваемых геосистем, не привязаны ни к каким реальным объектам, хотя могут быть совпадения в тех или иных деталях с объектами, известными читателю, и показывают только возможности новой ИТ. Тем не менее представленные результаты являются новыми знаниями о зарождении и развитии сложных процессов в различных средах, так как получены при решении уравнений, основанных на фундаментальных законах физики.
3. Задача № 1
3.1 Оценка условий формирования и эволюции
нижнемайтийного плюма
Понятие о мантийных плюмах было введено в 70-х годах прошлого века и впоследствии получило широкое распространение. Под этим термином понимают глубинные восходящие столбообразные, грибовидные и близкие к ним по форме структуры (струи), характеризующиеся перемещением мантийного вещества, обособленного от окружающей среды, повышенной температурой, пониженной плотностью и геохимическими особенностями [Красный, 2003]. Плюмы, как это установлено сейсмографическими и петролого-геохимическими исследованиями, связаны с процессами в переходном слое Б", разделяющем мантию и ядро, где идет преобразование вещества нижней мантии с разделением его на высокоплотное, опускающееся во внешнее ядро, и низкоплотное, скапливающееся в подошве нижней мантии [Лобковский и др., 2004].
Есть основание считать, что на ядро-мантийной границе формируются термохимические плюмы, определяющие тектонику горячих полей. Термохимический плюм может образоваться при наличии теплового потока из внешнего ядра и локальном поступлении химической добавки, понижающей температуру плавления вблизи подошвы нижней мантии. В области теплового пограничного слоя при понижении температуры плавления ниже температуры границы происходит плавление нижней мантии и подъем плюма до подошвы литосферы. Локальными источниками термохимических плюмов могут быть выступы на границе ядро - мантия высотой 10-20 км, а источником добавки, понижающей температуру плавления, - реакции железосодержащих фаз нижней мантии с водородом и метаном, выделяющимися на ядро-мантийной границе [Dobretsov et al.., 2003]. Эволюция мантийного плюма зависит от физических свойств области, в которой он зарождается, и свойств окружающей его мантии. Возникает вопрос, при каких условиях идет зарождение и подъем плюма и какие основные факторы влияют на этот процесс? В этой задаче на основе результатов численного эксперимента сделана оценка условий формирования и эволюции нижнемантийного плюма.
3.2. Схематизация задачи
Для расчета взята трехмерная область с декартовой системой координат в виде параллелепипеда размером 5000x5000, высотой 3000 км и разбита сеткой с шагом от 50 до 200 км на 41x41x40 объемных ячеек. Расчеты проводились на 100 - 200 млн лет с шагом по времени 10000 лет.
Предполагалось, что вещество мантии, при характерных для этого процесса времени и нагрузке, ведет себя подобно ньютоновской несжимаемой вязкой
жидкости. В рассматриваемой задаче внешних сил, возбуждающих вынужденную конвекцию, нет, поэтому источниковое поле не рассматривалось. В начальных условиях задавалось распределение температуры Т(х, у, ъ, t=0) = Т0(ъ), где Т0(ъ) - распределение температуры по вертикали в соответствии с температурным градиентом, и плотности среды, выраженной через концентрацию вещества С0(х, у, ъ, t=0)=C0(ъ), где С0(ъ) - вертикальное распределение концентрации вещества в зависимости от плотности, ъ - вертикальная координата, для функции тока задавались нулевые значения ^(х, у, ъ, t=0)=0. Вязкость в литосфере принималась
23 21 20
10 П-с, для верхней и нижней мантии 10 -10 П-с [Trubitsyn et al., 2006], коэффициенты температуропроводности и диффузии для всей области были одинаковые: 5=1х10-6 м2/сек, А,=1х10-8 м2/сек, а коэффициент температурного расширения а=5х10-5 °С-1.
На всех границах области для всех скоростей задавалось условие непротекания и проскальзывания
г=0. На верхней и нижней границах задавались постоянные значения температуры и плотности Т(х,у) | гв =Т>(0), Т(х,у) | гн=Т0(Ь+1), С(х,у) | гв=С>(0), С(х,у)[гн= =С0(Ь+1), а на боковых их распределение Т(х,ъ) | гп= =Т(х,ъ) | гз=Т(у,ъ) | гл=Т(у,ъ) | гп=Т0(ъ), С(х,ъ) | гп=С(х,ъ) | | гз=С(у,ъ) | гл=С(у,ъ) | гп=С0(ъ), Ь - количество ячеек по вертикальной координате.
На подошве нижней мантии, в ограниченной области, с помощью дельта-функции задавался источник тепла с температурой Т0=3200 °С и химической добавкой с концентрацией С1=1-20 мас. %, верхний предел концентрации выбран для оценки условий подъема низковязких плюмов. При этом максимальное отношение между плотностью нижней мантии и плотностью расплава (плавучесть) составляло Др=0.2, 0.35, 0.5 г/см3, а размер источника принимался 300x300 и 150x150 км, со временем существования от 40 до 80 млн лет.
Изменение плотности с глубиной задавалось с учетом ее скачков на границах фазовых переходов пород [Burmin, 2006]. В расчетах использовалась эффективная плотность, которая определялась через концентрацию литосферного или мантийного вещества C0=cxp'(ъ), где р'(ъ) - скачок плотности на заданной глубине, с - переводной безразмерный коэффициент (с=10-100). Считалось, что в точке с концентрацией флюида С1>1 % вещество среды приобретает свойства расплава - пониженную плотность и вязкость, т.е. С=С*, (Т,С1), где С* - значение концентрации вещества среды, при которой плотность снижается на Др, С1 - концентрация флюида. Отклонение плотности (плавучесть) определялось в приближении Буссинеска р(ъ)=рo(ъ)x(1-а(To-Т)-P(Co-C)), где р - переводной коэффициент от концентрации вещества к его плотности Р=1/с, а - коэффициент температурного расширения, С0, Т0, С, Т - исходные и текущие значения концентрации вещества и температуры в рассматриваемой
точке, ро(г) - плотность среды на определенной глубине (рис. 1). Наиболее важным, но плохо известным параметром для мантии является вязкость и ее зависимость от температуры, давления и концентрации химической добавки. Есть данные, что при повышении температуры среды на 100 °С вязкость понижается в 10 раз, т.е. на порядок [ТгиЬквуп еЬ а!., 2006]. В этом случае при изменении температуры в точке на 2000 °С, что может быть при подъеме плюма, вязкость должна измениться на 20 порядков. Тогда с учетом того, что при нисходящем потоке на границах тела плюма с окружающей средой (мантией) температура понижается, перепад вязкости может составлять более 20 порядков. В таких жестких условиях плюм подниматься не сможет (см. ниже разрушение плюма). Автором был выбран более мягкий вариант, заключающийся в том, что на начальном этапе, при изменении температуры на 100 °С, изменение вязкости шло со скоростью около порядка, но затем постепенно затухало, а максимальное ее снижение составляло не более десяти порядков. С учетом повышения вязкости за счет холодного нисходящего мантийного потока перепад вязкости не превышал 15 порядков. Тот же закон использовался для определения вязкости при изменении на 1 % концентрации химической добавки. Такую затухающую функцию можно описать с помощью гиперболического тангенса TANH(T, С). Эта зависимость выглядит так: ц=ц,0х10-А, где ц - новое значение вязкости, ц0 -исходное значение вязкости, А=тах( | АТ |, АС), Ат=TANH(ЛT/1000)хB, Ac=TANH(Сl/10)хB, АТ - перепад температуры в точке относительно исходной (заданной), С1 - концентрация химической добавки в точке, коэффициент В>1 характеризует степень изменения параметра А. Тогда максимальное снижение вязкости при концентрации химической добавки С1>20 % или при перепаде температуры Л^! 2000 °С |, для В=10, будет составлять не более десяти порядков. Температура плавления тугоплавкого слоя на глубине 100, 200 км принималась равной 1500, 1800 °С соответственно [БоЬгеЬвоу еЬ а!., 2006].
3.3. Анализ результатов
Результаты расчетов показали следующее. При наличии ограниченного источника химической добавки вблизи подошвы нижней мантии и теплового потока из внешнего ядра идет плавление вещества со снижением плотности и вязкости. При этом образуется канал, соизмеримый с площадью источника и увеличенной головной частью, заполненный расплавом. Внутри канала и в окружающем его массиве, за счет горизонтального градиента плотности, формируется конвективный поток с восходящей ветвью в канале и нисходящей в массиве. Это способствует интенсивному прогреву массива и проникновению в него химической добавки. При повышении концентрации химической
добавки в источнике (>10 %) канал расширяется за счет плавления боковых частей и нисходящая ветвь потока частично проникает в канал. Конвективный поток имеет торообразную форму. За счет вертикального градиента плотности возникают тангенциальные силы, которые в горизонтальном сечении закручивают верхнюю часть потока против часовой (циклонический вихрь), а нижнюю по часовой стрелке (антициклонический вихрь) (рис. 2). Скорость восходящего потока внутри канала достигает 50-100 см/год, а нисходящего в окружающем массиве в 20-100 раз меньше. В центральной части канала восходящий поток формирует высокотемпературную струю (до 3200 °С) с максимальной концентрацией химической добавки, за счет которой идет плавление пород в головной части плюма и его равномерный подъем. На периферии канала эти параметры расплава гораздо ниже. Нисходящяя ветвь конвективного потока снижает температуру в окружающих канал плюма породах мантии, что приводит к росту их вязкости. При этом, чем выше скорость конвекции, тем интенсивнее идет снижение температуры и рост вязкости. Это приводит к формированию в мантии областей с повышенной вязкостью вдоль всего пути подъема плюма, а в процессе образования линзы в подошве литосферы может возникнуть область с повышенной вязкостью кольцевой формы, сдерживающей ее рост. Области с повышенной вязкостью могут иметь различную форму в зависимости от динамики подъема плюма и места их расположения на его пути. Так, в нижней части образуются крестообразные, в средней - вихревые и в верхней - кольцевые формы высоковязких областей (рис. 3). Такие высоковязкие области вокруг канала плюма и сформировавшейся линзы могут служить одним из критериев для выделения его местонахождения и пути подъема при интерпретации сейсмических данных.
Если породы мантии имеют особую структуру (например, крупные зерна) [ТгиЬквуп еЬ а!., 2006], вязкость образовавшегося расплава понижается незначительно. Даже при максимальном повышении температуры ее значения в 10-100 раз меньше вязкости окружающего массива (высоковязкий расплав, в формуле для вязкости В=1-2), формируется ламинарный поток и подъем плюма идет равномерно. При концентрации в источнике химической добавки С1>5 %, поступающей постоянно в течение 80 млн лет, расплав начинает всплывать со скоростью около одного сантиметра в год. Постепенно скорость увеличивается и через 20 млн лет (для Ар=0.35, 0.5 г/см3) достигает 8-10 см/год и в дальнейшем меняется в пределах 5-10 %. Через 30-40 млн лет, в зависимости от плавучести, головная часть плюма, размером до 300-500 км, достигает литосферы, вязкость которой на 2-3 порядка выше, чем в мантии, и начинает растекаться по ее подошве, постепенно проникая вовнутрь. При уменьшении концентрации химической добавки в источнике до 1-2 % время подъема плюма возрастает до 60-70 млн лет.
Рис. 1. Начальные условия для температуры, плотности вещества, концентрации химической добавки и распределения вектора скорости Vmax=0.25 м/год (жёлтый цвет), Уш1п=0.013 м/год в расчетной области: а - распределение температуры Tmax=3200 °С, Tmin=0 °С; b - распределение плотности вещества среды max=3.4, min=2.6 г/см3; с - распределение концентрации химической добавки С1>1 %.
Fig. 1. Initial conditions for temperature, density and concentration of chemical agents, and the distribution of velocity vector Vmax=0.25 m/year (yellow), Vmin=0.013 m/year in the computation domain: a - distribution of temperature, Tmax=3200 °С, Tmin=0 °С; b - distribution of density, max =3.4, min=2.6 g/cm3; c - distribution of chemical agents concentration, C1>1 %.
tti-tttttt
+
i. -к- 4 4
\ It t t t
i T * » + *
t * * ?
М» //* t /
у- ^ y* jt
rf ^ л _ж- -—Л s-*
Ч 4 ■--^.'--t -* 4 **
У '+ '-t '+ —I- -+
- + --!■ --!■ -* ■+ ■+
H| t
* <k. 4. 1
i- V +
4 1 t t * Л .-+ *
V T ^ -t +
' r- r~ Г ^J J- r f
I
ну**
+ * * + ^и
M и il + ч ■
■+ ч. ч ■"-+ I ч- ч ■-+■-+ ■+
- Ч -+ % --t Ч Ч ■-+ ■"ч. Ч ■-I- "-!■ "-+■+■+
■ ч 41 N- ч ч
1'+ 41 \ 4 4
■■> + 1 4
* 4- +
Щ * 4 4 *
К- ^ s * *
Ш ^ v *
г -г" У г
L " '
L У г- У
Рис. 2. Распределение вектора скорости конвекции Vmax=0.56 м/год и концентрации вещества среды (голубым цветом выделен расплав) в увеличенных сечениях: а - горизонтального в головной части; b - вертикального по центру; с - горизонтального в тыловой (хвостовой) части плюма.
Fig. 2. Distribution of convention velocity, Vmax=0.56 m/year, and substance concentration (melt is shown in blue) in enlarged sections: a - horizontal head; b - vertical center, and c - horizontal rear (tail) parts of the plume .
Для Др=0.2 г/см максимальная скорость подъема плюма уменьшается до 5 см/год, характер конвективного потока остается прежним, а время подъема до подошвы литосферы увеличивается до 65 млн лет. Максимальная температура вещества, выносимая
плюмом к подошве литосферы, может достигать 3100 °С. При уменьшении времени работы источника с 80 до 20 млн лет время подъема плюма увеличивается в 1.5 раза, а температура вещества, вынесенного к подошве литосферы, уменьшается в 1.2-1.5 раза. Основное влияние на скорость подъема плюма, при рассмотренных выше параметрах, оказывает разность между плотностью расплава и окружающего массива - «плавучесть». При этом источник химической добавки, понижающей температуру плавления, может работать только на начальном этапе формирования плюма для его разгона, затем его движение вверх будет осуществляться за счет сухого плавления вещества мантии при возрастающем градиенте температуры (рис. 4).
При резком снижении вязкости расплава с ростом температуры (до 1 порядка на 100 °С, В=3-10) скорость восходящего потока в канале плюма возрастает, а нисходящего в окружающем массиве убывает. В этом случае режим плавления окружающего массива переходит из конвективного в кондуктивно-диф-фузионный, а при малой концентрации химической добавки подъем плюма осуществляется за счет кон-дуктивной теплопередачи и «сухого» плавления окружающего массива. Такой механизм движения может работать только при большом градиенте температуры на его кровле, что резко снижает вязкость в канале плюма и замедляет скорость его подъема, так как высокотемпературная струя упирается в высоковязкие породы мантии, с трудом проникая в них.
При дальнейшем подъеме температурный градиент продолжает расти (так как снижается температура окружающих пород), вязкость расплава - уменьшаться, скорость восходящего потока в канале возрастает, а нисходящего в окружающих породах снижается. Высокотемпературная струя начинает пульсировать, равномерный (ламинарный) поток переходит в турбулентный, «размазывая» температуру и химическую добавку по каналу, расширяя его. Температура и концентрация химической добавки, размазанные по каналу плюма, уменьшаются, скорость плавления мантии в его головной части снижается, что замедляет подъем плюма, он разрушается и, ослабленный, поднимается только до границы верхней мантии. При этом температура вынесенного вещества снижается на 600-800 °С, что может привести к формированию в подошве литосферы только горячего выступа, на 300-400 °С выше температуры верхней мантии (рис. 5).
Для опережающего проникновения химической добавки из головной части в окружающий массив необходимо значительно повысить ее концентрационный градиент. Увеличение концентрации химической добавки в источнике способствует росту ее градиента в головной части плюма и увеличению скорости ее проникновения в окружаюший массив, при этом в головной части образуется высоковязкий расплав. Дальнейшее уменьшение вязкости расплава за счет высокого градиента температуры, относительно вмещающих
b
4100 0 410 820 1230 1640 2050 2460 2870 3280 3690 4100g'
3280 ..........~T ,nf
■ ■;« ■ ■ 'Ж :
3000 4100
4100 KM
Рис. 3. Распределение вязкости и вектора скорости конвекции Vmax=0.92 м/год после формирования линзы расплава. Максимальный перепад вязкости 1010 (10 порядков). Срез сделан по вертикальной координате на глубине 500 км. Верхняя часть блока срезана по подошве литосферы.
Fig. 3. Distribution of viscosity and convection velocity vector, Vmax=0.92 m/year after formation of the melt lens. Maximum difference of viscosity 1010 (10 orders). The vertical profile goes at a depth of 500 km. The upper part of the block is cut off at the base of the lithosphere.
пород, может быть большим, то есть высоковязкий расплав формируется за счет опережающего поступления химической добавки, а затем снижение его вязкости идет за счет поднимающегося вслед высокотемпературного расплава, что способствует стабильному подъему плюма. Так, при максимальной разности вязкости расплава и окружающих пород на 5 порядков концентрация химической добавки должна составлять С1>12 %, а на 9 и более порядков - С1>18 %. Увеличение размера источника на эволюцию плюма оказывает незначительное влияние, но при этом возрастает количество выносимого к подошве литосферы вещества и его температура, что способствует более глубокому проникновению легкого вещества в литосферу и интенсивному ее прогреву, при этом у высоковязкого вещества скорость и глубина проникновения больше, чем у низковязкого.
После достижения головной частью плюма литосферы в ее подошве на глубине 200 км формируется
линза расплава, размер которой постепенно растет (обратно пропорционально вязкости расплава), а температура снижается. Через 80 млн лет линза, при минимальной вязкости, имеет размер до 4000 км в диаметре, мощность около 100 км и температуру расплава до 1800-2000 °С. Чем выше вязкость расплава (на 3-5 порядков ниже окружающего массива), тем меньше размер линзы и выше его температура. Так как высокая вязкость литосферы препятствует проникновению конвективных потоков, прогрев ее вышележащих слоев и перенос химической добавки идут только за счет медленных процессов кондуктивного теплообмена и диффузии, поэтому литосфера на глубине 100 км может прогреться только до 1200-1300 °С и расплав не появляется (рис. 6). Через 80-100 млн лет, в зависимости от размеров источника, формирующего плюм, температура падает ниже температуры плавления (1800 °С), но линза расплава может сохраниться, пока концентрация химической добавки не снизится до 1 %.
5200 о 820 1400 1810 2300 2710 3200 3610 4100 4820
км
5200 q
I 300 I 640 j. 1000 L 1340 I. 1700 I 2040
Конц-ция вещества н .000 м .100 .200 .300 .400 .500
I Э >00
I 3 ^900
L 3100 5200
5200 0 8,20 1400 1810 2300 2710 3200 3610 4100 4820 5200р
1 300
L 1340 [ 1700 I 2040
L 3100 5200
0 820 1400 1810 2300 2710 3200 3610 4100 4820
км 5200 q
3100 5200
Рис. 4. Динамика подъема плюма на три момента времени с распределением концентрации вещества среды (белым цветом выделен расплав) и вектора скорости конвекции (максимальная скорость - стрелки жёлтого цвета): а - через 24 млн лет, Vmax=0.28 м/год; b - через 32 млн лет, Vmax=0.42 м/год; с - через 48 млн лет, Vmax=0.8 м/год, после начала подъема плюма. Верхняя часть блоков срезана по подошве литосферы.
Fig. 4. Dynamics of plume rise at three time points, distribution of substance concentration (melt is shown in white), and convection velocity vector (yellow arrows show maximum velocities): a - 24 million years, Vmax=0.28 m/ year; b - 32 million years, Vmax=0.42 m/year; c - 48 million years, Vmax=0.8 m/year, from plume rise. The upper part of the block is cut off at the base of the lithosphere.
2460
1640
24&0
Скорость конвекции м/гол
Unax= .9068 Unax= 4534
Unax= 0906E Unax=.04533
Вязкость
среды
.00000
.00001
.00010
.00050
.00100
.00500
.01000
.05000
.10000
.50000
1.00000
Unax= . 9068 Unax= 4534
Unax= 0906E Uriax=. 04533 Uriax=.00906
Скорость конвекции м/гол
Unax= .9068 Unax= 4534
Unax= 0906E Unax=.04533 Unax=.00906
Конц-ция ким-добавки б 9
Контд-ция вещества в среде-
Рис. 5. Увеличенные вертикальные срезы по центральной части плюма и его окрестности (турбулентный режим) с распределением вектора скорости конвекции Vmax=0.98 м/год и: а - вязкости (максимальный перепад 9 порядков); b - химической добавки (концентрация в источнике 15 %); с - температуры; d - концентрации вещества среды (голубым и белым цветом выделен расплав); через 80 млн лет после начала подъема.
Fig. 5. Vertical profiles of the central part of the plume and its surroundings (turbulent regime), and distributions of convection velocity vector, Vmax=0.98 m/year: a - viscosity (maximum difference of 9 orders); b - chemical agents (concentration of 15% in the source); c - temperature, d -substance concentration (melt is shown in blue and white); 80 million years from commencement of plume rise.
L 1340 ¡_ 1700 L 2040 2400
[" 2740 Конц-ция J. 3100Еещ-ва E op^г
Рис. 6. Распределение: а - концентрации вещества среды (белым цветом выделен расплав); b - температуры Тш1п=200 °С, Tmax=3100 °С и вектора скорости конвекции Vmax=0.8 м/год (жёлтый цвет стрелок) в линзе расплава у подошвы литосферы через 100 млн лет.
Fig. 6. Distributions of a - substance concentration (melt is shown in white); b - temperature, Tmin=200 °C, Tmax=3100 °C and convection velocity vector, Vmax=0.8 m/year, in the melt lens at the base of the lithosphere in 100 million years.
В области расположения линзы, после исчезновения расплава, за счет повышенной температуры еще длительное время может сохраняться зона с пониженной вязкостью.
При наличии в литосфере ослабленных зон, в виде линейных структур с пониженной вязкостью, картина меняется. Прогрев и распространение химической добавки в пределах ослабленных зон идет за счет восходящих конвективных потоков. Это способствует в границах ослабленной зоны формированию на глубине от 100 км и выше линз расплава, которые при определенных условиях могут внедряться в горизонты земной коры. В нисходящих потоках холодный материал верхней части литосферы может проникать в верхнюю мантию. При этом очаги расплава возникают в восхо-
дящей и исчезают в нисходящей ветви конвективной ячейки в разных местах этих линейных структур, проникая в верхнюю часть земной коры. Время существования очагов расплава достигает 5-10 млн лет. Минимальная вязкость расплава в этих очагах может снижаться на 8 порядков относительно вязкости литосферы, а максимальная концентрация легкой химической добавки - достигать 6 % при исходной концентрации в источнике 12-20 % (рис. 7).
3.4. Заключение
Результаты численного эксперимента при использовании трехмерной математической модели показали, что за счет потоков тепла и химической добавки из
Рис. 7. Распределение концентрации вещества среды (белым цветом выделен расплав) после формирования линз расплава вдоль линейных ослабленных зон через 100 млн лет. Срез сделан по вертикальной координате на глубине 100 км.
Fig. 7. Distribution of substance concentration (melt is shown in white) after formation of melt lenses along linear weakened zones in 100 million years. The vertical profile goes at a depth of 100 km.
внешнего ядра идет плавление пород нижней мантии с формированием восходящего торообразного потока ее вещества. Если образуется высоковязкий расплав на 12 порядка ниже вязкости окружающего массива, то подъем плюма идет равномерно (ламинарный поток) и достигает подошвы высоковязкой литосферы, имея температуру до 3100 °С, а источник химической добавки может работать только на начальном этапе его формирования с концентрацией не более 5 %. Если вязкость расплава с ростом температуры снижается быстро (до 1 порядка на 100 °С), то при низкой (менее 10 %) концентрации химической добавки в источнике за счет большого перепада вязкости между расплавом и окружающим массивом плюм может разрушиться. Для его стабильного, равномерного подъема необходима высокая (до 15-20 %) концентрация химической
добавки в источнике.
Плавление литосферы и внедрение ультрабазит-базитовых расплавов в верхнюю часть литосферы и земную кору возможны при прогреве и распространении химической добавки в пределах ослабленных зон, происходящих за счет восходящего конвективного потока. При этом нисходящий поток холодного материала верхней части литосферы может проникать в верхнюю мантию. При отсутствии ослабленных зон плавление пород литосферы не происходит, так как высоковязкая литосфера, которая на 5 порядков и более выше вязкости вещества, вынесенного плюмом, препятствует конвективному переносу вещества, а за счет медленных процессов тепловой кондукции и концентрационной диффузии ультрабазит-базитовые расплавы не формируются.
Рис. 8. Начальные условия для температуры и плотности вещества в расчетной области: а - распределение температуры Tmax=1500 °С, Tmin=0 °С, температура вещества линзы T0=2000 °С; b - распределение плотности вещества среды max=3.4, min=2.6, плотность вещества линзы р0=3.0 г/см3.
Fig. 8. Initial conditions for temperature and substance density in the computation domain; a - distribution of temperature Tmax=1500 °С, T^^O^, temperature of the lens substance, T0=2000 °С; b - distribution of substance density, max=3.4, min=2.6; density of the lens substance, р0=3.0 g/cm3.
4. Задача № 2
4.1. Оценка условий формирования и развития диапиров
в системе литосфера - кора
Всплывание менее плотных (легких) пород через более плотные (тяжелые) в поле силы тяжести называется диапиризмом. К диапирам относятся различные восходящие куполообразные структуры, в том числе гранитогнейсовые купола и гранитоидные батолиты. Для формирования диапира необходим легкий питающий слой, перекрытый более тяжелыми слоями пород, и наличие неоднородности (неравномерность толщины питающего или перекрывающего слоя) или нарушения сплошности [Рамберг, 1985]. Питающим слоем может являться линза расплава, сформированная нижнемантийным плюмом в подошве высоковязкой литосферы,
а неоднородностью - вогнутость этой подошвы. В этой задаче на основе результатов численного эксперимента рассмотрено влияние линзы расплава с различной вязкостью, плотностью и концентрацией флюида на формирование диапиров в земной коре и литосфере.
4.2. Схематизация задачи
Для расчета была выбрана трехмерная область с декартовой системой координат в виде параллелепипеда с основанием 2000x2000 км, высотой 250 км и разбита сеткой с шагом по времени 10000 лет. В рассматриваемой задаче внешних сил, возбуждающих вынужденную конвекцию, нет, поэтому источниковое поле не рассматривалось.
Все допущения, ограничения, коэффициенты тем-
пературопроводности и диффузии, уравнение для вязкости принимались такие же, как в первой задаче, а начальные и граничные условия в таком же виде, но с учетом расположения рассматриваемой расчетной области (верхняя мантия - кора).
Значения плотности составляли 2.6 г/см3 в верхней и 3.4 г/см3 в нижней части расчетной области. В нижней части литосферы задавалась линза расплавленного мантийного вещества диаметром 1600, мощностью 100 км, с температурой Т0=2000 °С и концентрацией флюида С1=2-6 мас. % (см. результаты примера № 1). При этом максимальная разность между плотностью среды и плотностью расплава (плавучесть) составляла ЛР=0.2-0.4 г/см3 (рис. 8).
4.3. Анализ результатов
Результаты расчетов показали следующее. На контакте линзы расплава с областью вогнутости литосферы за счет прогрева литосферы и проникновения флюида формируется восходящий конвективный поток вихревой структуры. Со временем менее интенсивные потоки формируются по всему контакту линзы с литосферой. Поток образует куполообразную структуру вещества с пониженной плотностью, которая начинает всплывать в верхнюю часть рассматриваемой области со скоростью 2-3 см/год, максимальная скорость может достигать 6-8 см/год. За счет пониженной плотности и вязкости формируется конвективный поток торообразного вида с закручиванием в горизонтальной плоскости во фронтальной части диапира против часовой и в тыловой по часовой стрелке (аналогичный потоку задачи № 1). Конвективный поток высокотемпературного и флюидонасыщенного вещества частично захватывает окружающий массив, что способствует его частичному плавлению, за счет прогрева и внедрения флюида, и перемешиванию с веществом линзы. Чем ниже вязкость вещества диапира, тем интенсивнее идет смешение. При определенных условиях (сравнительно низкой вязкости пород литосферы) может сформироваться диапир с несколькими симметрично расположенными восходящими куполообразными структурами (рис. 9).
Если при подъеме диапира составляющие его породы (расплав) претерпевают фазовый переход на границах скачков плотности вмещающих пород, то его плавучесть (исходный перепад плотности) сохраняется, и он может всплыть до уровня верхней коры. Дальнейшему его подъему будет препятствовать повышенная вязкость пород верхней коры, но при наличии в ней ослабленных зон (разломов) он может выйти на поверхность земли. Если породы диапира не претерпевают фазового перехода на границах скачков плотности, то диапир поднимается только до уровня, на котором плотность его вещества становится равной или выше (в зависимости от вязкости) плотности вмещающих пород, и затем начинает оседать (растекать-
ся), перемешиваясь с окружающими породами.
Существенное влияние на формирование и подъем диапира оказывает вязкость. Если вещество диапира претерпевает фазовые переходы, то относительно высокая его вязкость (на 1-3 порядка ниже вязкости окружающих пород) замедляет его подъем, но при этом диапир имеет компактную (аккуратную) форму с равномерным распределением температуры и вещества. Если вещество диапира не претерпевает фазовых переходов, то при этой вязкости диапир может внедриться в породы с меньшей плотностью, чем его вещество, и подняться в них тем выше, чем выше его вязкость. После внедрения диапир начинает оседать (растекаться) и перемешиваться с окружающими породами, при этом, чем выше вязкость его вещества, тем меньше скорость оседания и перемешивания (рис. 10).
Низкая вязкость (4 порядка и более ниже окружающих пород) вещества диапира, претерпевающего фазовые переходы, ускоряет его подъем, а диапир принимает форму струи с повышением значений температуры и концентрации флюида в центральной части. Максимальная температура вещества диапира и концентрация флюида при уменьшении вязкости снижаются за счет более высокой скорости конвекции. Диапир при низкой вязкости вещества без фазовых переходов поднимается только до уровня равных по плотности окружающих его пород. В этом случае скорость растекания вещества диапира и перемешивания с окружающими породами больше, чем у диапира с высокой вязкостью. Это способствует относительно быстрому изменению формы диапира из куполообразной в линзовидную.
При подъеме диапира происходит смешение его вещества с окружающими породами (ассимиляция). Интенсивность смешения без фазовых переходов вещества диапира при его подъеме мало зависит от вязкости.
Максимальная величина смешения концентрируется на периферии диапира. При оседании диапира на контакте с его головной частью идет затягивание верхнего более легкого слоя вмещающих пород вниз. Чем ниже вязкость, тем сильнее расплывается диапир и глубже затягивает верхние более легкие породы. Если при подъеме диапира происходят фазовые переходы, то процесс смешения идет тем интенсивнее, чем меньше вязкость его вещества. При высокой вязкости диапир движется медленнее, раздвигая окружающие его породы с незначительной их ассимиляцией, и более интенсивное затягивание всех слоев окружающих пород идет на периферии диапира (рис. 11).
4.4. Заключение
При наличии в вогнутой подошве литосферы линзы расплава с повышенной температурой, содержащей флюид, может возникнуть вихревой поток, формирующий восходящую куполообразную структуру -
Рис. 9. Распределение температуры после подъёма диапира при различной вязкости окружающих пород: а - вязкость окружающих пород равна значению, заданному в исходных данных, Tmax=1730 °С; b - вязкость окружающих пород в десять раз меньше, Tmax=1690 °С.
Fig. 9. Distribution of temperatures after diapir rise for different viscosities of the surrounding rocks: a - viscosity of the surrounding rocks is equal to the value specified in the initial data, Tmax=1730 °С; b - viscosity of the surrounding rocks is lower by a factor of 10, Tmax=1690 °С.
диапир. Если породы (расплав) диапира на границах скачков плотности вмещающих пород претерпевают фазовый переход, то при наличии в верхней коре разломов он может выйти на поверхность земли. Если фазового перехода нет, то диапир может подняться только до уровня, на котором плотность его вещества становится выше плотности вмещающих пород и затем начинает оседать (растекаться).
Высокая вязкость вещества диапира (на 1-3 порядка ниже вязкости окружающих пород) замедляет его подъем при фазовом переходе вещества и его оседание (растекание) без фазового перехода. Относительно низкая вязкость его вещества (на 4 порядка и более ниже окружающих пород) ускоряет подъем в первом случае и оседание диапира во втором.
При подъеме диапира происходит смешение (ассимиляция) его вещества с окружающими породами. Этот процесс идет более интенсивно при условии фа-
зовых переходов вещества диапира и меньшей его вязкости.
5. Задача № 3
5.1. Оценка условий формирования рифтовых зон
Рифт - линейно вытянутая щелевидная или ровооб-разная структура глубинного происхождения. Существуют разные гипотезы образования рифтовых зон. Это, например, идея пассивного рифтинга, которая объясняет растяжение земной поверхности напряжениями в литосфере, вызванными межплитными взаимодействиями. При этом плиты могут приходить в движение при воздействии на них внешних сил. Модели пассивного рифтинга противопоставляется концепция активного рифтинга, согласно которой раздвижение лито
a
Рис. 10. Распределение плотности вещества среды после подъема диапира при относительно высокой вязкости диапира (1-3 порядка ниже вязкости окружающих пород): а - вещество диапира претерпевало фазовые переходы; b - вещество диапира не претерпевало фазовые переходы.
Fig. 10. Distribution of substance density after diapir rise for relatively high viscosity of the diapir (1-3 orders of below the viscosity of the surrounding rocks): a - diapir substance was subject to phase transitions; b - diapir substance was not subject to phase transitions.
сферы происходит в результате воздействия мантийного диапира. В подтверждение этой идеи существует ряд геофизических наблюдений, указывающих на расположение у границы Мохо линзы низкоскоростного вещества, интерпретируемого как альтернативная мантия. Существует еще одна модель образования рифтов, основанная на идее подъема верхнемантийного диапи-ра с последующей частичной кристаллизацией вынесенного им тяжелого мантийного вещества и погружением его в легких, коровых породах [Рамберг, 1985]. В данной задаче представлена оценка условий, при которых возможно формирование рифтовых зон.
5.2. Схематизация задачи
Для расчета в декартовой системе координат была взята трехмерная область в виде параллелепипеда размером от 2000x2000 до 4000x4000 и высотой от 250 до 300 км. Она была разбита сеткой с шагом от 5 до 100
км на 44x41x41 объемных ячеек. Расчеты проводились на 80-800 млн лет с шагом по времени 10000-100000 лет. Верхние четыре точки с нулевой плотностью имитировали слой атмосферы, в котором вязкость принималась на 5-6 порядков ниже, а коэффициенты температуропроводности на порядок выше. Вязкость в литосфере принималась 1022 П-с, для верхней и нижней коры 1024-1023 П-с.
На всех боковых границах области для нормальных скоростей задавалось условие протекания, а для касательных (тангенциальных) скоростей прилипания | г=0, для температуры и плотности их распределение по вертикали. На верхней и нижней границе для скоростей задавались условия непротекания и проскальзывания г=0, а для температуры и плотности постоянные значения. Поверхность верхней коры считалась свободной, а ее граница меняла свое положение в зависимости от концентрации среды (плотности). При увеличении концентрации вещества в четырех
Рис. 11. Смешение вещества диапира с окружающими породами при условии фазового перехода: а - начальное распределение пород по плотности, вертикальный срез по центру блока; вязкость пород диапира ниже вязкости окружающих пород: b - на 5 порядко в 105 раз; c -на 3 порядка в 103 раза; d - на 1 порядок в 10 раз.
Fig. 11. Mixing of the diapir substance with the surrounding rocks in case of phase transition; a - initial distribution of rocks by density; the vertical profile in the center of the block; the viscosity of the diapir rocks is lower than the viscosity of the surrounding rocks: b - by 5 orders (by a factor of 105); c - by 3 orders (by a factor of 103); d - by 1 order (by a factor of 10).
верхних точках (вздымание коры) вязкость в них росла по линейному закону пропорционально коэффициенту отношения концентрации вещества в рассматриваемой точке и коре под ней. При уменьшении концентрации в коре (погружение коры) вязкость по тому же закону снижалась. Такое изменение вязкости можно записать в виде уравнения ц^ц^к, где ц - значение вязкости в корректируемой точке, ц - исходное значение вязкости в коре под рассматриваемой точкой, к=С/С0, где С - значение концентрации вещества в корректируемой точке, С0 - исходное значение концентрации вещества в той же точке коры при погружении или в коре под корректируемой точкой при вздымании, с ограничением на коэффициент 1>к>0.
В расчетах использовалось три варианта форм ли-тосферной плиты: сплошная однородная плита, эта же плита с линейной ослабленной зоной в ее центральной части, вязкость которой была в 50-100 раз меньше вязкости основной части плиты, и область, состоящая из двух плит с разной плотностью и (или) вязкостью, имеющих Б-образный контакт.
Для сценария пассивного рифтинга задавалось дви-
жение плит навстречу друг другу с постоянной скоростью за счет внешнего источника, а контакт между плитами имел Б-бразную форму. При линейном контакте между плитами движение задавалось расходящееся и навстречу друг другу. Для модели активного рифтинга в нижней части литосферы задавалась линза расплавленного вещества, вынесенного нижнемантийным плюмом к подошве литосферы, диаметром 15002000, мощностью 70-100 км, с температурой Т0=1800-2000 °С и концентрацией флюида С1=2-6 мас. % (см. пример №1). При этом максимальное отношение между плотностью среды и плотностью расплава (плавучесть) составляло Др=0.2-0.5 г/см3. Линза имела изо-метричную и продолговатую (вытянутую в одном направлении) форму (рис. 12).
5.3. Анализ результатов
Результаты расчетов показали следующее. Для реализации сценария пассивного рифтинга необходимы высокие скорости движения плит (десятки см/год), так как в гравитационном поле при более низких скоро-
Рис. 12. Начальные условия задачи: a - распределение температуры с изометричной формой линзы расплавленного мантийного вещества; b - распределение концентрации вещества (плотности) с продолговатой формой линзы расплавленного мантийного вещества; c - распределение вещества и горизонтальный срез верхней части коры с S-образной формой контакта между плитами и заданными скоростями движения за счет внешнего источника для сценария пассивного рифтинга.
Fig. 12. Initial conditions of the problem: a - temperature distribution for an isometric-shape lens of mantle melt; b - distribution of substance density of an elongated lens of mantle melt; c - substance distribution, the horizontal profile of the upper crust with the S-shaped contact between plates, preset movement velocities due to an external source for the passive rifting scenario.
стях движения формирующийся рифт заполняется веществом коры. Окружающие породы за счет гидростатических сил выдавливаются с подошвы и бортов в
Рис. 13. Распределение концентрации вещества (плотности) при формировании рифта за счет внешних сил при S-образной форме контакта двух плит (чем светлее окраска, тем меньше концентрация вещества, значит, больше глубина погружения). Скорости движения плит навстречу друг другу: a - 4 м/год; b - 0,7 м/год; c - 0,2 м/год.
Fig. 13. Distribution of substance density during rifting due to external forces at the S-shaped contact between two plates (lighter colours show lower densities, i.e. larger depths). Velocities of plate movements towards each other: a - 4 m/year; b - 0.7 m/year; c - 0.2 m/year.
зону формирования рифта, что не позволяет ему развиваться (рис. 13).
В модели активного рифтинга на контакте линзы расплава с областью вогнутости литосферы идет прогрев и проникновение флюида. При этом возникает градиент плотности, за счет которого формируется восходящий конвективный поток вихревой структуры - диапир. Это конвективный поток торообразного вида с закручиванием в горизонтальной плоскости во фронтальной части диапира против часовой (циклонический вихрь), а в тыловой по часовой стрелке (антициклонический вихрь). Он захватывает окружающие породы, что способствует их частичному плавлению за счет прогрева и внедрения флюида и перемешива
Рис. 14. Распределение: а - концентрации вещества (плотности), светлой окраской выделяются области погружения; b - температуры (на подошве самого глубокого рифта область повышенной температуры продолговатой формы) при подъеме диапира.
Fig. 14. Distributions of a - substance density (sinking areas are shown by light colour) and b - temperature (there is an elongated area of higher temperatures at the bottom of the deepest rift) during diapir rise.
нию с веществом линзы, которое при низкой вязкости идет более интенсивно, чем при высокой. Со временем менее интенсивные потоки могут сформироваться по всему контакту линзы с литосферой. Эти потоки могут образовать несколько куполообразных структур вещества, имеющих пониженную плотность, вязкость и всплывающих с различной скоростью. Скорость подъема вещества может меняться от 0.1 до 5.0 см/год в зависимости от его плавучести и вязкости вмещающих пород. Время подъема диапира может варьироваться в широких пределах и составлять от 60 до 200 млн лет.
При подъеме диапира происходит вздымание коры в ее центральной части и опускание на периферии. В верхней части коры идет раздвижение пород с формированием прогиба, который постоянно заполняется породами, выдавливаемыми диапиром, поэтому фор-
мирование рифта из прогиба не происходит. Температура поднимающегося вещества может превышать значения вмещающих пород на 100-800 °С в зависимости от скорости подъема диапира и температуры его источника (линзы исходного мантийного вещества) (рис. 14).
Чем ниже вязкость вещества диапира, тем интенсивнее оно перемешивается с окружающими породами, во время его подъема. При высокой вязкости смешение происходит только в краевых частях диапира, в этом случае он в основном будет состоять из мантийного вещества. При высокой плавучести вещество диапира может выйти на поверхность, а при низкой завершить подъем в коре. По мере снижения температуры и частичной кристаллизации вещества формируются пластичные породы, имеющие зоны различной
Рис. 15. Рифтовые зоны (выделены серым и голубым цветом), образованные при погружении тяжёлого вещества, вынесенного диапиром: а - литосфера неоднородная с S-образным контактом двух плит; b - литосфера однородная с ослабленной линейной зоной в центре.
Fig. 15. Rift zones (shown in grey and blue) that formed during sinking of heavy materials drawn by diapirs; a - the lithosphere is not homogeneous, and there is an S-shaped contact between two plates; b - the lithosphere is homogeneous, and there is a weakened linear zone in the center.
плотности. В тех зонах, где преобладают породы, вынесенные из мантии и имеющие высокую плотность, может возникнуть нисходящий поток за счет их погружения в более легких породах. Погружение может идти неравномерно, более плотное вещество будет погружаться быстрее менее плотного, при этом верхние более легкие слои вмещающих пород будут затягиваться нисходящим потоком вниз, за счет чего будет формироваться рифтовая зона (рис. 15).
Форма линзы играет определенную роль при формировании диапира и рифтовой зоны. В сплошной однородной литосфере форма диапира и рифтовой зоны примерно соответствуют форме линзы. При введении вытянутой ослабленной зоны с пониженной вязкостью над продолговатой линзой формируется рифт примерно такой же формы, а над изометричной линзой - изогнутой неоднородной формы. В литосфере, состоящей
из разных плит с Б-образным контактом, из симметричной линзы может подняться диапир вдоль этого контакта и иметь продолговатую форму, но может сформироваться симметричный с несколькими головными частями диапир, а рифтовые зоны будут иметь примерно такие же формы. В этом случае формы диа-пира и рифтовой зоны будут в основном зависеть от параметров мантии (плотность, вязкость), степени плавучести вещества диапира при подъеме и степени его кристаллизации при погружении.
Чем выше вязкость тяжелого вещества, тем дольше оно будет погружаться (растекаться) и глубже затягивать верхние более легкие породы. Форма рифтовой зоны будет зависеть от размеров и формы погружающихся частей диапира. На бортах формирующейся рифтовой зоны будет возникать обратный восходящий поток вещества, который приведет к вздыманию ее
краевых частей и подъему изотермы. При наличии в коре ослабленной области с вязкостью, в 50-100 раз меньшей окружающих пород, погружение их в ее центральной части будет идти более интенсивно, чем на периферии, и может сформироваться рифт, ограниченный бортами с более жесткими породами.
Можно выделить три этапа развития рифтовой зоны: 1) подъем верхнемантийного диапира с возможным выходом его на поверхность, 2) остывание с частичной кристаллизацией вынесенного им мантийного вещества, 3) погружение тяжелого, пластичного вещества с затягиванием за собой верхних более легких слоев вмещающих пород. Продолжительность каждого из этапов зависит от параметров коры, мантии, линзы расплавленного вещества (плотность, вязкость, концентрация флюида) и может меняться от 60 до 250 млн лет.
5.4. Заключение
Для реализации сценария пассивного рифтинга необходимы высокие скорости движения плит (десятки см/год). В модели активного рифтинга при раздвиже-нии пород верхней коры образуется прогиб, который заполняется выдавливаемыми диапиром породами, и развития из него рифта не происходит. Формирование рифта возможно при погружении частично закристаллизованных пластичных мантийных пород, вынесенных диапиром.
6. Выводы
На основе изложенного материала можно сделать следующие заключения и выводы.
1. Предложена новая информационная технология, с помощью которой можно решать сложнейшие задачи теоретического и прикладного характера по изучению и прогнозированию развития процессов в геосистемах.
2. Ядром технологии является трехмерная сопряженная математическая модель тепломассопереноса, основанная на новой системе уравнений гидродинамики, позволяющей разделить конвекцию на две составляющие - свободную и вынужденную, и, вычислив их значения, получить параметры суммарного конвективного потока. Такой подход дает бо-
лее адекватное приближение модели к реальному объекту.
3. Основную нагрузку при интерпретации огромного объема результатов численного эксперимента несет пакет графических программ, позволяющий в динамике просматривать значения любой вычисленной функции, наблюдать за развитием взаимосвязанных процессов на изучаемом объекте и при необходимости определять точки ввода элементов управления ими.
4. Применение данной технологии помогает понять причинно-следственные связи сложных природных явлений и на этой основе генерировать новые знания.
7. Литература
Гебхарт Б., Джалурия Й, Махаджан Р., Саммакия Б. Свободно-конвективные течения, тепло- и массообмен. В 2-х книгах. Пер. с англ. М: Мир, 1991. 678 с.
Гунин В.И. Новая трехмерная математическая модель тепломассопереноса в пористых средах и ее возможности // Геоэкология. 2003. № 4. С. 355-370.
Кирюхин А.В., Сугробов В.М. Модели теплопереноса в гидротермальных системах Камчатки. М.: Наука, 1987. 150 с.
Красный Л.И. Эволюция тектонических идей от середины XIX столетия до современности. СПб. 2003.
Лобковский Л.И., Никишин А.М., Хаин В.Е. Современные проблемы геотектоники и геодинамики. М.: Научный мир, 2004. 610 с.
Рамберг Х. Сила тяжести и деформации в земной коре. Пер.с англ. М: Недра, 1985. 399 с.
Тарунин Е.Л. Численный эксперимент в задачах свободной конвекции, Иркутск: Изд-во Иркутского государственного ун-та, 1990. 223 с.
Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1975. 735 с.
Burmin V. Yu. Distributions of Density and Elastic Parameters in the Earth // Izvestiya, Physics of the Solid Earth. V. 42. № 7. 2006. P. 608-620. doi: 10.1134/S106935130607007X.
Dobretsov N.L., Kirdyashkin A.A., Kirdyashkin A.G. Physicochemical Conditions at the core-mantle boundary and formation of thermochemical plumes // Doklady Earth Sciences. V. 393A. № 9. 2003. P. 1319-1322.
Dobretsov N.L., Kirdyashkin A.A., Kirdyashkin A.G. Diameter and formation time of plume head at the base of refractory lithospheric layer // Doklady Earth Sciences. 2006. V. 406. № 1. P. 56-60. doi:10.1134/S1028334X06010144.
Trubitsyn V.P., Baranov A.A., A.N. Evseev A.N., Trubitsyn A.P. Exact analytical solutions of the stokes equation for testing the equations of mantle convection with a variable viscosity // Izvestiya, Physics of the Solid Earth. V. 42. № 7. 2006. P. 537-545. doi: 10.1134/S1069351306070019.
Гунин Владимир Иванович, директор
Независимый коммерческий центр моделирования геосистем «МоГеос» 670034, Улан-Удэ, пр. 50 лет Октября, 38-18, Россия H e-mail: [email protected]
Gunin Vladimir I., Director
MoGeos Independent Commercial Centre for Modeling of Geological Systems 670034, Ulan-Ude, pr. 50 let Oktyabrya, 38-18, Russia H e-mail: [email protected]
