автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Численное моделирование тепловых конективных течений в верхней мантии Земли

РТВ ОН

И

На правах рукописи

Рычкова Блена Владимировна

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕПЛОВЫХ КОНВЕКТИВНЫХ ТЕЧЕНИЙ В ВЕРХНЕЙ МАНТИИ ЗЕМЛИ

05.13.16 — применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Новосибирск - 1997

Работа выполнена в Институте вычислительных технологий СО РАН (г. Новосибирск)

Научные руководители:

доктор физико-математических наук, профессор Г.Г. Черных, кандидат геолого-минералогических наук, с.н.с. С.А. Тычков

Официальные оппоненты:

доктор физико-математических наук, профессор А.Ф. Воеводин, кандидат геолого-минералогических наук, с.н.с О.П. Полянский.

Ведущая организация — Институт теоретической и прикладной механики СО РАН (г. Новосибирск).

Зашита диссертации состоится " 24" февраля 1998 года в " 14м" часов на заседании диссертационного совета Д 002.10.02 при Институте вычислительной математики и математической геофизики СО РАН по адресу: 630090, Новосибирск-90, проспект академика Лаврентьева, 6.

С диссертацией можно ознакомиться в читальном зале библиотеки ИВМиМГ СО РАН, проспект акад. Лаврентьева, 6.

Автореферат разослан "20" января 1998 г.

Ученый секретарь диссертационного совета

к.т.н.

Г.И. Забиняко

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность проблемы. Интенсивное использование математических методов в естественных пауках п, в частности, в геологии началось с середины 60-х годов нашего столетия, когда произошло бурное разлитие вычислительной техники и расширение возможностей современных ЭВМ. Существенный вклад в развитие этого направления внесли A.C. Алексеев. М.М. Лаврентьев, Э.Э. Фотиадп, А.Ф. Белоусов, Ю.А. Воронин, А.Н. Дмитриев н другие. Многие направления наук о Земле попугали ускоренное развитие с привлечением математической обработки данпых, а некоторые из них поднялись на нозую качественную ступень в своем развитии. Это относится прежде иссго к изучению глубоких недр Земли, где главную роль, как известно, играют сейсмические исследования, поставляющие информацию о составе и состоянии вещестаа, а также о процессах, протекающих во всем земном шаре, вплоть до ядра. Благодаря современным математическим методам1 обработка данных сейсмологии и сейсморазведки мы сейчас имеем представление о распределении физических параметров с глубиной во всей Земле, а развитие сейсмической томографии в последнее десятилетие позволяет говорить о возникновении нового направления в пауках о Земле — глубинной геодинамики, которое синтезирует всю доступную информацию о недрах с целью выявления принципиальных механизмов, отвечающих за структуру и динамику недр. Специфика глубинной геодинамики состоит в том, что в отличие от рада других естественных наук, здесь не удается моделировать в лабораторных условиях большинство процессов, определяющих тектонический режим и эволюцию нашей планеты. Поэтому математическое моделирование этих процессов является весьма актуальной задачей. Математические модели дпнампкп недр обязательно должны включать в себя расчет геологических п геофизических параметров, что позволяет при сравнении их с наблюдениями, составить впечатление о степени адекватности модели реальным процессам. Кроме того, при формулировке граничных и начальных условии необходимо максимально пепользовать

'Алексеев A.C., Рябой В.З. - Модель строения верхней мантии по объемным сейсмическим волнам //Строение земной коры и верхней маитпи по данным сейсмических исследований. Киев: Наукоиа дунка. 1977. С. 67-82.

всю имеющуюся информацию об объекте исследования. Настоящая работа посвящена моделированию мантийной динамики внутриконтинентальных областей. Здесь важно отметить, что большая часть созданных ранее математических моделей процессов в недрах относится к океаническим областям, а появляющиеся сейчас подобные модели для континентов зачастую механически переносят отработанные для океанов технологии и подходы создания моделей недр. В существующих моделях динамики континентальных недр (Flcitout., Yuen2, Schmcling, Maiquart3,4) областям относительного погружения поверхности соответствует литосфера в 200 км толщины. Это означает, что до сих пор нет модели, геофизические характеристики которой даже на качественном уровне соответствовали бы наблюдениям. Развиваемый в настоящей работе подход обеспечивает подобное соответствие (относительно пониженный рельеф поверхности участков с тонкой литосферой, что подтверждается данными сейсмических наблюдений), что несомненно является достижением п создании моделей такого класса.

Цель работы состоит в разработке численных моделей мантийной динамики внутриконтинентальных областей, в частности, тепловой конвекции в верхней мантии северной половины Евразийской плиты.

Научная новизна изложенных в диссертационной работе результатов заключается в следующем:

• литосфера рассматривается как жесткий кондуктивный слой, не принимающий участия в конвективном перемешивании вещества;

о реализован учет концентрации радиоактивных изотопов в коре континентов (исходя из современных работ по теплозьш моделям региона, созданным на основе многолетних наблюдений за тепловым потоком Сибири);

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

•'l'Vitoui |„, Уши D.A. Secondary convod.ion and the growth of the oceanic lithosphere //I'iiysifs of 1 If; Ivirtli and Planetary Interiors. H)M. V. 36. P. 181-212.

'S*liinclini" II., Mar*|Uart. (j. Tlio influence of scr-oiid-scalo conversion on the thickncss of

• onl Illental Iii husplxTc and crust//'lb.toiiopliy.sirx li»l. V. 18!). I'. 281 .'!0(i.

1S< ImnTiri^ II., Marrjuarl (I. Maritl'1 (low and evolution of the lithosphere //Physics of the brill an<l I'hnmtary Interiors. 1!Ш. V. 7!). P. 211 267.

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

Достоверность полученных результатов подтверждается проведением

многочисленных тестовых расчетов с применением мер контроля точности получаемых решений н удовлетворительным согласованием результатов с теоретическими и экспериментальными данными других авторов (для модельных задач).

Научная и практическая ценность работы. Разработанная численная модель может быть использована для численного моделирования широкого круга задач геодинамики, а также при прогнозировании землетрясений а в сейсморайонироваяии.

Представленные в диссертации исследования проводились а рамках программы СО РАН "Новые поколения вычислительной техники, математическое моделирование и информационные технологии": "Численное моделирование течений вязкой жидкости и турбулентных течений" (Лз гос. регистрации 01.9.40 000839), "Разработка и исследование математических моделей п численных методов решения задач аэро-гпдродинамики" (№ гос. регистрации 01960011628). Работа поддержана Российским фондом фундаментальных исследований (гранты № 94-05-16544, № 96-05-65970, № 96-05-66227). Международным научным фондом (грант ЛКЮО).

На защиту выносятся:

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

• разработка эффективного алгоритма, на основе апробированных конечно-разностных методов, для численного решения двумерных уравнений Навье-Стокса в приближении ОберСека-Буссинеска при больших числах Рэлея (Ю4 — 106) и бесконечном числе Прандтля, позволяющего варьировать большое число входных геофизических данных;

• результаты числешгого .моделирования тепловой конвекции под литосферой постоянной п переменной мощности;

• результаты сопоставления вычисленных геофизических характеристик (тепловой поток, рельеф, гравитационное поле) с реально наблюдаемыми на поверхности Земли для верхнемантийной конвекции под континентальной плитой Северной Евразии.

Апробация работы. Основные научные результаты диссертации до-

кладывались на Межреспубликанской Школе-семинаре по численным методам механики вязкой жидкости (Новосибирск, 1994), Пятой Международной конференции памяти Л.П.Зоненшайна по тектонике плит (Москва, 1995), Международной конференции "Математические модели и численные методы механики сплошных сред" (Новосибирск, 1996), Международной Летней школе "Конвекция п геофизике и астрофизике" (Аоста, Италия, 199G), Восьмой Международной конференции по методам аэрофизпческих исследований (Новосибирск, 1996), Научной конференции РФФИ "Геодинамика и эволюция Земли" (Новосибирск, 199G), Сибирской школе-семинаре "Математические проблемы механики сплошных сред" (Новосибирск, 1997), обсуждались на семинарах в Объединенном институте геологии, геофизики и минералогии СО РАН под руководством акад. Н.Л. Добреиова, в Институте вычислительных технологий СО РАН под руководством акад. Ю.И. Шокина и проф. D.M. Ковснн, в Институте, вычислительной математики и математической геофизики СО РАН под руководством акад. A.C. Алексеева.

Публикации. Основное содержание диссертации отражено в 7 работах, список которых приведен в конце автореферата.

Структура и объем диссертации. Работа состоит из введения, трех глав, заключения ц списка литературы из 121 наименования. Полный объем диссертации 121 страницы, включая 30 рисунков и 10 таблиц.

СОДЕРЖАНИЕ РАБОТЫ

Но Введении обосновывается актуальность темы. Приводен обзор работ. посвященных 'пиленному моделированию естественно-конвективных процессов п пязкой жидкости (прежде всего, работы В.И. Полежаева и

его учеников5,6, E.J1. Тарушгна7), изучению процессов тепловой конвекции в недрах Земли методами математического моделирования (работы В.П.Трубицына8,9 ц других). Сформулированы цели работы и приведено краткое содержание глав.

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

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

Численное интегрирование уравнения теплопереноса осуществлялось методом предиктор-корректор со схемой расщепления10 на этапе предиктора и аппроксимацией конвективных слагаемых направленными разностями. На этапе корректора применялась центрально-разностная аппроксимация конвективных слагаемых. При интегрирования уравнеппя переноса вихря использовалась конечно-разяоегная схема переменных направлений с центральными разностями по пространству. Уравнение Пуассона для функции

5Пасконов В.М., Полежаев В.И., Чудов Л. А. Численное моделирование процессов тепло п массообмена. М.: Наука, 1981.

'Полежаев В.И., Бунэ A.B., Всрезуб H.A. и др. Матсмыпческое моделирование тепломас-собмеяа на основе уравнений Назье-Стокса. М.: Наука, 1987.

7Таруяип Е.Л. Вычислительный эксперимент в задачах свободной конвехпии. Иркутск: Нзд-во Иркут. ун-та, 1990.

sTrubitsyn V.P., Rykov V.V. 3D numerical model of the Wilson cycle //Geodynamics. 1995. V. 20. P. 63-75.

9Трубицьга В.П., Рыков B.B., Трубииын А.П. Конвекция и распределение вязкости в мантии //Физика Земли. 1997. № 3. С. 3-10.

10Я,ненко H.H. - Метод дробных шагов решешм многомерных задач математической физики. Новосибирск: Наука, 1967.

тока решалось итерационным методом неременных направлений.

В случае зависимости вязкости от температуры и давления для численного решения уравнения четвертого порядка с переменными коэффициентами для функции тока использовался итерационный метод стабилизирующей поправки10,11.

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

Во втором пункте приведены тестовые расчеты конвективных течений для различных значении чисел Прандтля Рг л Рэлся Ra. Число Ra определяется по отношению к критическому числу Rac и 657.5, при котором равновесие подогреваемой снизу жидкости теряет устойчивость13. Рассматривается двумерная конвекция в квадратной области между свободными границами (при подогреве снизу) жидкости с постоянной вязко-пью. В ходе расчетов определялось число Нуссельта Nu, характеризующее отношение «сродненного по горизонтали потока тепла при наличии конвекшш. к потоку тепла, переносимому теплопроводностью. В табл. 1 приведены значения числа Нуссельта в зависимости от числа Прандтля (Rii/Ibv « 50). Во второй строке табл. 1 представлены значения числа Нуссельта из работы"'^, в третьей строке — результаты тестирования численной модели на равномерной сетке (hT = hy = 1/60). Использовалась постановка задачи для (?;), îj, Т) системы. При вычислении Nu применялась

11 Houston M.П., Jr., lie Hreiniif-rkcr ,1'1. Л1)1 solution of frnoconvcctiori in a variable viscosity llni.l //.I. of Computational Physics. I!)74. V. Jfi. 1'. 221-2-4!).

'-Ilnp'on К.Г.. Численный меч с« д.'1я мпач никого потока //"Механика", 1965. №6. I'. 77.

Moon' U.K.. \V< iss X.O. 'I4vo-rliiiicri.sioija! Uaylcigli I!éuar<I «mvcx.tion //.!. Fluid Moch. I')7Л. V. .VS. I'.wl >. I». Ж1 :ш.

квадратурная формула Спмпсона.

Таблица 1 Таблица 2

Pr 0.1 1. 100.

Тест 7.92 7.81 7.30

61 x 61 7.14 7.13 7.12

Ra/Rac 20 100 200 400 1000

Тест 5.37 9.25 11.52 14.45 19.47

61 x 61 5.34 9.20 11.54 14.41 19.31

В случае бесконечного числа Прандтля численные эксперименты проводились для (-ф,Т)-системы при следующих числах Рэлея: 1.315 ■ 1С1. 6.575 • 10\ 1.315 • 10®, 2 .63 • 105, 6.575 • 105. Их результаты представлены в табл. 2 в виде значений числа Нуссельта в зависимости от Ra/Ra Таблицы 1, 2 свидетельствуют о достаточно хорошем соответствии расчетов настоящей работы и данных^'3'.

Третий пункт посвящен сопоставлению с данными физических экспериментов (которые были любезно предоставлены автору д. т.н. Кирдяшкп-ным А.Г.), выполненных при Рг = 28.7 и Ra = 1.1 • Ю4. В экспериментах, проведенных B.C. Берднпковым, А.Г. Кпрдяшкиным1,1 обнаружено, что прп и = const в области значений 5 • 103 < Ra < 2 ■ 104 (Рг = 14) образуется квазидвумерная стационарная валпковая структура течения. Следовательно, применение двумерной математической модели для сравнения с результатами физического эксперимента является допустимым. Из представленных на рис. 1 значений горизонтальной компоненты вектора скорости в вертикальном сечении конвектирующей ячейки (сечение проходит через точку минимума функции тока) видно' вполне удовлетворительное соответствие экспериментальных н расчетных данных. Вычисления были реализованы в переменных у!>, и\ в качестве граничного условия для температуры на дне области задавалась экспериментальная кривая зависимости от времени.

В четвертом пункте рассматривается модельная постановка задачи для вязкости, зависящей от температуры: v(T) — ехр[С(0.5 — Т)], С = const (Torrance, Turcotte15).

"Бердкиков B.C., Ккрдишкия А.Г. - Структура свободпо-конвективпых течений в горизонтальном слое жидкости при различных граничных условиях //Структура пристенного пограничного слоя. Новосибирск: ИГФ СО АН СССР. 1978. С. 5-48.

"Torrance К.Е, Turcotte D.L. - Thermal convection with large viscosity variations //J. Fluid Mech. 1971. V. 47. Part 1. P. 113-125.

so o.s мм/с

Рис. 1.

Профиль горизонтальной компоненты скорости в вертикальном сечении ячейки:

• - эксперимент (Кирдяшкин А.Г.), — - численное моделирование.

В данном случае используется система уравнений Навье-Стокса в переменных ф, и\

а / > ^ öt „ г

Д(м) + Ra— + 2

Д ф = -ui,

<Pv du д^и dv d2v dx* ду d,ß Ох дхду

(dv ЗгЛ] _

д~у ~ Ш/1 ~

о.

(1) (2)

Алгоритм определения значений температуры, вихря и функции тока на новом шаге по времени состоит в проведении совместных итераций по нелинейности. Здесь также было получено удовлетворительное соответствие с результатами'15'. Однако, при решении реальных задач геодинамики вязкость имеет существенные изменения своих значений в полосе малой толщины (большие градиенты), вследствие этого, наличие вторых производных от вязкости в уравнении для завихренности (2) может привести к расходимости итерационного алгоритма при проведении расчетов (что подтверждено численными экспериментами). Неприемлемой такая постановка является также в случае разрывного коэффициента вязкости. Поэтому более предпочтительным представляется другой подход, а именно, отказ от введения вихря скорости. Тогда для функции тока получается уравнение четвертого порядка с переменными коэффициентами.

Сопоставление с системой тестой Blankenbach et al.lc, ориентированных

'Hlfitikctiharh Ii. г! al. Л bciirhmark comparisott for (nantlc convcction Codes //Gcophys. J. Im. im v. 98. I'. z< :»*.

н

на задачи мантийной конвекции осуществлено в пятом пункте. Данное сопоставление проведено для (ф, Т)-системы с уравнением четвертого порядка для ф. Расчеты настоящей работы выполнялись на сетках, сгущающихся в области пограничных слоев (у верхней и нижней границ области): 41 х 69 п 81 х 69.

В результате численных экспериментов определялись следующее »сличили: число Нуссельта ¡\ги, среднеквадратичная скорость Кггм. безразмерный градиент температуры в углах ячейки д, экстремумы Ге температуры около нижней и верхней границ на центральной линии и их положение ус. Табл. 3 составлена для случая постоянной вязкости: Па = 10е, 1/(1 = 1. Данные табл. 4 соответствуют случаю и = ¡/„ ехр [— В ■ Т/9 + С ■ (1 — ;/)/(!}, В - 1п(16384), С = 1п(64), и0 = 2.5Лй19м2[с, 6. - 106 м, в = 1000 К, 1(й = 2.5.

Таблица 3, Па = 106, и = 1

Nu Vr^i q-2 % У'

Тест 21.972 833.989 45.9G4 0.877 0.432 0.058

41 х 69 20.283 733.209 41.251 0.756 0.426 0.065

Таблица 4, Ra = 10\ v = и(Т, у)

Nu vrms Я1 «2 9з 94

Тест 6.929 171.755 18.48 0.177 14.168 0.618

81 х 69 6.674 158.794 15.15 0.209 13.719 0.779

Результаты вычислений демонстрируют удовлетворительное соответствие с данными работы^16).

В шестом пункте второй главы представлены результаты моделирования конвективных процессов в жидкости с сильно меняющейся нелинейной вязкостью. Вид зависимости вязкости от температуры и давления взят из работы Fleitout, Yuen'2'. Численное интегрирование уравнения теп.топереноса осуществлялось с применением конечно-разностного алгоритма предиктор-корректор с продольно-поперечной прогонкой на этапе предиктора и аппроксимацией конвективных слагаемых направленными разностями. Прп решении уравнения четвертого порядка для функции тока использовался итерационный метод стабилизирующей поправки. Получено достаточно хорошее

соответствие результатов расчетов настоящей работы и данных'2' но распределению вязкости. А. в случае применения уравнений (1), (2) получить положительный результат не удалось.

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

В первом параграфе приведен обзор современных представлений о режимах тепловой конвекции в мантии Земли (D.P. McKenzie, J.M. Roberts, N.O. Weiss17; В.П. Трубицын, B.B. Николаичик18; H.J1. Добрецов, А.Г. Кирдяшкин19; В. Travis, P. Olson20 и другие). Рассматривается влияние числа Рэлея и начальных условий на эволюцию течения.

Второй параграф содержит описание постановки задачи. "В приближении Оборбека -Буссинеска для бесконечного числа Прандтля уравнения движения и переноса тепла задаются следующим образом в безразмерном виде:

от иот .

Ot Эх

от а (.эт\ д ,

1 а2 Г . д'Ф 1 __ Клдт

ду1 с)!1] \ \ду'1 Ох'1)) дхду\ дхду ) Ох'

Оф Оф .

где х, у оси декартовой системы координат, х = у = 0 в левом верхнем углу, ось у направлена вверх; и, у - компоненты вектора скорости; X время; Т температура; а - коэффициент температуропроводности: А скорос ть генерации радиоактивного тепла; ф - функция тока; ¡/(.г. у) кинематическая вязкость; 11а - число Рэлея.

В системе (3) (о) используется линейное уравнение состояния: о = о,\\ - п(Т - Т„)].

,7.М< Krnzit! IJ.J'., Roberts J..VI, Weiss N.O. Convection in Earth's mantle: towards a numerical simulation //.!. Fluid Mech. НЖ. V. 62. I'art 3. P. 465-538.

'"Трубицын 11.11., Никшлй'/ик li.ll. Режимы тсиловой конвекции //Физика З'.'мли. ¡991. .V (,. (,'. .'{ 12.

'Vloi/piW/H П Л.. Кирляшкич Л.Г. Глубинна* «.-«динамика. Новосибирск, 19.91. ' I)iivis IV. Olson I'. C<;nv.T.iioii wi1.li intimia! Iiral swirri-s anil thermal turbulence in the Earth's iii,.iill<- //С.-.,pliys. .). In1«;r. I!)!M. V. 1IK. I>. 1 19.

Обезразмерпвание величин выполнялось следующим образом (штрихом обозначены безразмерные величины, в (3)-(5) знак обезразмеривания опу-щеп):

(х,у) = d{x,y)\ v=v0-i/, Г = dfi0T', А = A'Qo/d (u,v) = (и, v)'a0/d, x¡¡ — a0 • ф', t = t'd1/a0.

Здесь d - характерный размер области, в данном случае он равен ^верхней мантии ^литосферьч ^о -V« ¡Qoср - характерное значение коэффициента температуропроводностии (термодиффузии), ср - удельная теплоемкость прп постоянном давлений, д0 - плотность при постоянной средней температуре Та, \м - коэффициент теплопроводности в мантии, в коре литосферы à равно отношению коэффициентов теплопроводности коры (Лк) и мантии, i/0 = \ín¡ о о - характерное значение вязкости (fia - динамическая вязкость), Д, = —дТ/ду - вертикальный градпенг температуры, Qa = Хм/30 - тепловой поток на дне верхней мантии.

Значения параметров среды в модели: d = 660 км - 120 км = 540 км, а0 = 8 • 10~7д(2/с, Qa = 27 • 10~3Втп/м2, Хм = 4Вт/(м°С), Хк = 2.2Вт/(м°С), Д, = 6.75 • 10"3oC/.ti, а = 3 ■ 10~51/°С, ju„ = 1021Яа • с, g = 9.8-w/c2, £>а = 3.3 • 103кг/л3; Ra = gago0od4f(aoii 0) = 7.103 • 105.

Геометрия расчетной области выбирались следующим образом: ¿z = 6000 к.и - горизонтальный и = 660 к.и - вертикальный размеры области, Lx/Ly =9:1. Граничные условия приведены на рис. 2.

-200-

s

•i-

(0 X -400-

f

с -600-

|120к

2160 км I

3660 Ш Г= 0. и = 0, и = 0 6000 км

литосфера

160 - 280 км

u — 0t v=0

Конвектирующая мантия

ах и дх V

Рис. 2. Геометрия области и граничные условия.

Начальное распределение температуры было выбрано в виде Некой!;, Уиеп'2|(г = —у, т.е. ось г направлена вниз):

где Г„ = 1350°С, Ть = 1800СС, гь = 660 км, \а = 14.04 М лет

(1 М лет = 10б лет и 3.15 • 1013 с). Тогда дТ/дг = ,% при г = 660 км.

После обезразмеривания Тн к нему добавляется малое возмущение в (СЬ^евдеп21):

^iflsin^V Ю Ьк \LyJ

( кжх | cos ——

*=1 - V '

Третий параграф посвящен изучению эволюционного поведения конвективных течений под литосферной плитой постоянной мощности. Анализ эволюции течения показал, что в качестве начального условна для расчетов течения под литосферной плитой переменной мощности следует выбирать неустановившуюся структуру течения с развитыми пограничными слоями, отвечающей абсолютному времени процесса в 2 3 • 103 М лет (рис. 3). Стационарное состояние достигается после ~ 7000 М лет и имеет вид одной вытянутой ячейки (рис. 4),

В четвертом параграфе рассматривается зависимость структуры течения от мощности литосферной плиты (в центральной части области). Тол-шина литосферы менялась от 160 до 280 км. Показано, что утолщенная до 160 км континентальная литосфера не влияет на структуру верхнемантийного течения — после 5000 М лет остается одна ячейка, длина которой соответствует длине расчетной области (как и для ровной: литосферы, см. рис. 4), т.е. амплитуда латеральных вариаций мощности литосферы в 40 км недостаточна для возмущения течения.

Во втором эксперименте мощность толстой литосферы была увеличена до Лк = 170 км. Такое небольшое увеличение, менее чем на 10%, привело к принципиальным изменениям в истории развития структуры течения. Начальный этап до ~ 2 000 М лег как и ранее отличался уменьшением количества ячеек с 6 до 3. Его продолжительность возросла почти вдвое по

2,C'hristensen U. - Convection with pressure- and temperature-dependent non-Newtonian rlie-ology //Geophys. J. Roy. Astr. Soc. 1984. V. 77. P. 343-384.

сравнению с предыдущими расчетами. Другой его отличительной особенностью является достаточно быстрое смещение восходящего потока конвекции, который находился первоначально на удалении в ~ 300 км, под область утолщенной литосферы. Это смещение произошло за t ~ 500 М лет. А через t ~ 800 М лет от начала счета дпссипировал первоначально существовавший в центральной частя толстой плиты нисходящий поток конвекции (рис. 5). К локальным, но принципиальным изменениям в структуре тепловой конвекции при росте мощности литосферы следует отнести развитие конвективной неустойчивости у краев толстой плиты (рис. 6). Впервые развитие подобной неустойчивости наблюдалось при толщине плиты в 190 км. Причина образования локальных конвективных ячеек у краев плиты состоит в существовании выраженного горизонтального градиента температуры при движении от холодной толстой плиты в сторону конвектирующей мантии. Таким образом, пижний горизонтальный поток локальной ячейки выносит мантийный материал из-под утолщенной литосферы к подошве тонкой литосферы на расстояние от края толстой плиты, равное горизонтальному размеру ячейки, что составляет первые сотни километров. Средняя скорость перемешивания вещества в локальной ячейке не превосходит средней скорости глобальной конвекции, т.е. ~ 1 — 2 см ¡год.

Изучение роли толстой литосферы в эволюции структуры течения состояло в последовательном увеличении толщины плиты. Были рассчитаны варианты для 180 км, 190 км, 200 км, 220 км, 240 км и 280 км толщины (рис. 7). По мере увеличения толщины плиты устойчивый восходящий поток под толстой литосферой мигрировал к центру утолщенного участка а при hK = 280 км практически достиг его. С увеличением толщины сокращалось также время появления восходящего потока под толстой литосферой. Так для hK = 170 км оно составляло ~ 500 М лет, для 200 км снизилось до 400 М лет, а при толщине в 280 км восходящий поток смещался под толстую плиту за 200 М лег. Уменьшалось и время диссипации изначально существующего под толстой плитой нисходящего потока. Оно снизилось почти в четыре раза и для hK = 280 км составляло всего 280 М лет.

В пятом параграфе представлена численная модель тепловой конвекции

в верхней мантии северной половины Евразийской плиты, в которой реа-

лпзован учет концентрации радиоактивных изотопов в коре континентов, исходя из работ Дучковаи др.22'23 по тепловым моделям региона, созданным на основе многолетних наблюдений за тепловым потоком Сибири. Генерация внутренних источников тепла задавалась в виде двух- и трехслойной модели. Для Западно-Сибирской плиты в первом слое Ьц = 14 км теплоге-нерация составляла А\ = 0.7 • 10~6Вт/м3, в следующем слое, расположенным ниже Лг = 17 км, Аг = 0.3- 10~6Вт/м3. Кора Сибирской платформы была представлена трехслойной моделью: Ь.х = 14 км, А\ = 0.6 • 10~^Вт/мг; к2 = 21 км, А2 = 0.2 ■ 10~6Вт/м3; /г3 = 10 км, А3 = 0.04 • 10~6Вт/м3. В результате вычислений наблюдается формирование восходящего конвективного потока под утолщенной лптосферой. Эта структура течения обес-печшзает пониженные значения гравитационного поля и теплового потока на поверхности, и относительно повышенный рельеф поверхности, что характерно для докембрийских континентальных платформ (рис. 8).

В Заключении сделаны следующие выводы:

1. Построена численная модель конвективных процессов в верхней ман-тип Земли. Модель основана на уравнениях Навье-Стокса в приближении Обербека-Буссинеска. Осуществлено подробное тестирование этой модели.

2. Выполнено численное моделирование тепловой конвекции под лито-сферной плитой постоянной мощности. Анализ эволюции течения позволил определить состояние структуры, которое было выбрано в качестве начального при изучении конвекции под плитой переменной мощности.

3. Проведено систематическое изучение влияния кондуктивной плиты на пространственно-временную эволюцию конвекции при постепенном увеличении толщины кратона от 160 до 280 км. Показано, что латеральные неоднородности стабилизируют структуру течения, что является принципиальным для понимания причин унаследованности и периодичности теологических процессов на континентах.

4. Представлена физико-математическая модель тепловой конвекции в верхней мантии северной половины Евразийской плиты, позволяющая удо-

?2Дучков А .Д., Соколова Л.С. - Геотермические исследования в Сибири. Новосибирск: Наука, 1974.

"Дучков А .Д., Балобаев В.Т., Володько Б.В. и др. - Температура, криолитозона и радиогенная генерация в земной коре Северной Азии. Новосибирск: ОИГГМ СО РАН, 1994.

влетворнтельно рассчитывать геофизические характеристики, соответствующие наблюдаемым на поверхности.

Данная модель континентальных недр принципиально отличается от ранее предложенных в следующем:

— литосфера рассматривается как жесткий коядуктивнын слой повышенной прочности, имеющий химические отличия по веществу от остальной конвектирующей мантии. Поэтому вещество литосферы в данной модели не принимает участия в конвективном перемешивании вещества;

— астеносфера моделируется слоем, имеющем зависимость величины вязкости от глубины;

— реализован учет концентрации радиоактивных изотопов в коре континентов (исходя из современных работ по тепловым моделям региона, созданным на основе многолетних наблюдений за тепловым потоком Сибири).

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

Автор выражает глубокую признательность научным руководителям д.ф.-м.н. Геннадию Георгиевичу Черных и к.г.-м.н. Сергею Анатольевичу Тычкову. Автор также благодарит д.т.н. Анатолия Григорьевича Кирдяш-кина (ОИГГиМ СО РАН) за предоставленные экспериментальные данные, Александра Николаевича Василевского (ОИГГиМ СО РАН) за предоставленные результаты расчетов ряда геофизических характеристик. Отдельная благодарность к.ф.-м.н. Жанне Леонидовне Коробицыной за ценные советы, высказанные в процессе работы над диссертацией.

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ:

1. Рычкова Е.В. Численные эксперименты в задачах свободной конвекции //Материалы XXXII Международной научной студенческой конференции "Студент и научно-технический прогресс": Математика. Новосибирск: Нзд-во Новости, ун-та, 1994. С. 52.

2. Мошкин Н.П., Рычкова Е.В., Тычков С.А., Черных Г.Г. Тестирование некоторых численных моделей конвективных течений применительно к задачам геодинамики //Вычислительные технологии. Новосибирск: ИВТ СО РАН. 1995. Т. 4, № 13. С. 224-231.

3. Рычкова Е.В., Тычков С.А. Численное моделирование тепловой конвекции в слое жидкости с сильной зависимостью вязкости от температуры и давления //Математические модели и численные методы механики сплошных сред. Тезисы докладов международной конференции. Новосибирск, 1996. С. 454.

4. Chernykh G.G., Moshkin N.P., Rychkov«i E.V., Tychkov S.A. Comparison of some numerical algorithms for two-dimensional convection of fluid with nonlinear viscosity //Proceedings of 8th International Conference on the Methods of Aerophysical Research. Novosibirsk, 1996. Part I. P. 79-84.

5. Тычков C.A., Рычкова Е.В. Континентальная геодинамика: роль литосферы в эволюции мантийных течений //Геодинамика и эволюция Земли. Материалы к научной конференции РФФИ. Новосибирск, 1996.

С. 40-41.

6. Рычкова Е.В., Тычков С.А. Численная модель тепловой конвекции в верхней мантии Земли под литосферой континентов //Вычислительные технологии. 1997. Т. 2, № 5 С. 66-81.

7. Тычков С.А., Рычкова Е.В., Василевский А.Н. Взаимодействие плюма и тепловой конвекции в верхней мантии под континентом //Геология и геофизика. 1998. № 4 (в печати).

ПОЛЕ СКОРОСТИ

3.98 см/гоб

% £ * 9 » * « « 1 4 * т 4

! 1} $ } И

-- г г ! г £

пшШш!

И!

о

-200 -400 -С00

1000 ' 2000 ' 3000 ' '.ООО 5000 6000

Рис. 3. Мощность литосферы Нл = 120 км; ? = 2500 М лет.

ПОЛЕ СКОРОСТИ -. 4.42 см/год

1000 - 2000 ФУНКЦИЯ Т0КЙ

Г-| «Г Ь«т Г,-.

ЗООО 4 ООО 5000

6000

1000 2000

0 750 1200 1300 1400 1500 1600 1700"С

4000 5000

т-г I I

Рис. 4. Мощность литосферы кл = 120 км; I = 8000 М лет.

о

ПОЛЕ СКОРОСТИ

4.22 см/гоЗ

-200 *

-400

-600

ШШШШ1ШЖН П/Ш

ft«*«-» . > — *->-.-.•• * - . i ' • ....

О 1000

ФУНКЦИЯ TOKfi

T»it f

2000

зооо

4000

5000

6000

1000 2000 3000 4000 ' '5000 '......6000

Рис. 5. Мощность кратона hK — 170 км; t = 840 М лет.

функция тока температура

О Т60 1200 1300 1400 1500 1600 1700*С

Рис. 6. Мощность кратона hK = 240 км; t = 4 000 М лет. Средняя часть области (в рамке) представлена вверху в увеличенном виде.

¡H-r-f

5000 6000 (П

íooo то' зооо' '.ooo •то

t*mw д**,-** iiii

О 750 1200 1300 1400 1500 1600 1700'C

Рис. 7. t = 4 000 M лет. Мощность кратона: (1)-180 км, (2)—190 км, (3)-200 км, (4)-220 км, (5)-240 кя, (6)-280 км

Рис. 8. Мощность кратона кк = 200 км; £ = 10000 М лет. Иа = 8.82 ■ 105, тепловой поток на дне 28 • 10~*Вт/м2. Тепловой поток на поверхности, аномалии гравитации и топография:

----наблюдения, Дучков22,23;--расчет.

Изолинии температуры (в градусах Цельсия).