CC BY
17
GEODYNAMICS & TECTONOPHYSICS
PUBLISHED BY THE INSTITUTE OF THE EARTH'S CRUST SIBERIAN BRANCH OF RUSSIAN ACADEMY OF SCIENCES
ISSN 2078-502X
2019 VOLUME 10 ISSUE 1 PAGES 21-41
https://doi.org/10.5800/GT-2019-10-1-0402
Mathematical models simulating the formation
of the stress-strain state of epiplatform orogens D. S. Myagkov, Yu. L. Rebetsky
O.Yu. Schmidt Institute of Physics of the Earth of RAS, Moscow, Russia
Abstract: The sources of the natural stress-strain state (SSS) of epiplatform orogens are investigated by tectonophy-sical methods based on seismological data. According to the available data, the horizontal axes of the main deviatoric extension are dominant in depressions, while in the ridges of the orogens, the axes of the main deviatorial compression are dominant. Our comparative analysis is focused on SSS of the orogenic crust. It is generally accepted that the sources of such SSS are geodynamic processes, including the pressure on the Eurasian Plate from the Indian Plate, and the small-scale thermogravitational asthenospheric convection. In the mathematical (analytical) simulation technique used in our study, the main criterion for the correctness of models in terms of tectonophysics is the correspondence between the orientation pattern of the principal stress tensor axes in the crust model to the natural data. According to Model I, the lithospheric SSS under lateral compression is less consistent with the sought-for SSS. Model II also gives the results that do not fully correspond to the stress data from tectonophysical reconstructions. However, additional analysis suggests that asthenospheric convection is a more promising (from the point of view of tectonophysics) geodynamic process for explaining epiplatform orogenesis. In our opinion, more complex and probably non-analytical mathematical models should consider this source of loading of the lithosphere as one of the most significant factors in the formation of the orogenic crust SSS in Central Asia.
Key words: analytical modeling; epiplatform orogenesis; stress state of mountain ranges
RESEARCH ARTICLE Received: August 3, 2017
Revised: February 8, 2019 Accepted: February 20, 2019
For citation: Myagkov D.S., Rebetsky Yu.L., 2019. Mathematical models simulating the formation of the stress-strain state of epiplatform orogens. Geodynamics & Tectonophysics 10 (1), 21-41. doi:10.5800/GT-2019-10-1-0402.
Математическая модель формирования
напряженно-деформированного состояния
эпиплатформенных орогенов
Д. С. Мягков, Ю. Л. Ребецкий
Институт физики Земли им. О.Ю. Шмидта РАН, Москва, Россия
Аннотация: Рассматривается вопрос об источниках формирования природного напряженно-деформированного состояния (НДС) эпиплатформенных орогенов, исследуемого тектонофизическими методами на основе сейсмологических данных. Такого рода данные свидетельствуют о преобладании горизонтальной ориентации осей главного девиаторного растяжения во впадинах и осей главного девиаторного сжатия в хребтах орогенов. Проводится сравнительный анализ НДС коры орогена, источником которого выступают «общепринятые» геодинамические процессы: давление на Евразийскую плиту со стороны Индийской и действие мелкомасштабной термогравитационной астеносферной конвекции. Исследование проводилось методом математического (аналитического) моделирования, основным критерием тектонофизической корректности модели считалось соответствие распределения ориентаций главных осей тензора напряжений в коровой части моделей природным данным. Моделирование показало, что НДС литосферы, формирующееся в обстановке латерального сжатия, менее соответствует искомому. Вторая модель также показала результаты, не вполне соответствующие данным тектонофизической реконструкции напряжений. Однако дополнительный анализ позволил установить, что астеносферная конвекция является более перспективным, с точки зрения тектонофизики, геодинамическим процессом для объяснения эпиплатформенного орогенеза. Мы считаем, что в рамках более сложных и, вероятно, неаналитических математических моделей этот источник нагруже-ния литосферы должен рассматриваться как один из наиболее существенных факторов формирования НДС коры орогенов Центральной Азии.
Ключевые слова: аналитическое моделирование; эпиплатформенный орогенез; напряженное состояние горных массивов
1. Введение
Проблема происхождения эпиплатформенных орогенов представляет особый интерес, так как, с одной стороны, по данным структурам накоплено большое количество и геологических [Khain, Lomize, 1995; Belousov, 1989; Burtman, 2012], и тектонофизи-ческих данных. С другой стороны, простые схемы, объясняющие механизм их формирования [Cloetingh et al., 2002] обычным латеральным сжатием, не подтверждаются ни современными данными о закономерностях распределения напряжений в коре [Rebetsky, Alekseev, 2014], ни численными моделями [Mikhailov, 1999; Mikhailov et al., 1999; Timoshkina et al., 2010], которые зачастую также не согласуются с этими данными [Liu et al., 2007] либо создавались до их получения [England, Houseman, 1986]. Распространенность эпиплатформенных орогенов на территории России (Алтай, Саяны) и постсоветского пространства (Тянь-Шань) придает созданию подобной модели практическую важность. Такая модель могла бы лечь в основу более детальных региональных тектонофизических моделей, на основе которых, в первую очередь, можно было бы выполнять прогноз
сейсмической опасности. Отметим, что данное исследование ориентировано на приложение к проблемам тектонофизики орогенов Центральной Азии, однако рассматриваемые модели, при некотором расширении, актуальны и для континентальных орогенов в целом, в первую очередь - для орогенов Высокой Азии (Памир, Тибет) и активных континентальных окраин (Анды и пр.), по которым также имеются данные тектонофизических реконструкций природных напряжений и исследованию проблем генезиса которых отводятся ведущие роли в современной геодинамической науке.
На текущий момент существуют самые разнообразные модели формирования эпиплатформенных орогенов. Говоря о геодинамических источниках данного процесса, следует выделить внешние по отношению к участку литосферы, деформирование которой приводит к возникновению орогена, и внутренние, связанные с перераспределением энергии в самом рассматриваемом участке литосферы. В качестве первых рассматривается в основном термогравитационная мелкомасштабная конвекция в астеносфере, формирующая вынужденное течение в вышележащем искомом участке
литосферы, и коллизионное давление на окраины литосферной плиты со стороны иных плит, вызывающее деформацию внутренней части искомой плиты в обстановке латерального сжатия. При рассмотрении орогенов Центральной Азии большое значение придается, конечно, взаимодействию Индийской и Евразийской плит. Источники второй группы - это деформации, вызванные инверсивным распределением плотности и упругих модулей в коре и литосфере [Makarov, 2005; Dobretsov et al., 1995], гравитационной неустойчивостью в системе литосфера - астеносфера, существенной латеральной неоднородностью плотности литосферы и перераспределением внутренней энергии, связанной, к примеру, с метаморфическими и геохимическими процессами, протекающими под воздействием тепловой энергии, поступающей из астеносферы. В работах [Rebetsky, 2012; Rebetsky, Pogorelov, 2013] на основе результатов решения простых аналитических моделей обсуждалась проблема несоответствия некоторых моделей процесса формирования внутриконтинентальных орогенов имеющимся тектонофизическим данным [Rebetsky, 2015]. С целью уточнения сделанных в этих работах выводов были проведены дополнительные исследования по созданию набора моделей «геодинамический процесс - НДС литосферы» как первого этапа в решении обратной задачи о механизме формирования внутриплитных орогенов.
В текущем исследовании рассматриваются процессы первого типа: 1) мелкомасштабная астено-сферная конвекция - модель I; 2) связанная с коллизией деформация внутриплитных областей (которая, с плейт-тектонических позиций, вызывается действием общемантийной конвекции [Trubitsyn et al., 2006]) - модель II. Исследование ведется математическими методами, позволяющими получить решение в аналитическом виде (конечные формулы). Допустимость применения аналитических методов, являющихся более простыми по реализации, но ограниченными по сложности используемых моделей среды, по отношению к численным методам, позволяющим использовать сложные модели, объясняется вышеуказанным отсутствием даже приблизительной начальной модели и точного знания реологических свойств литосферы (вязкость, к примеру, известна с точностью лишь до порядка), поэтому численный расчет заранее предполагает наличие последующего перебора пара-
метров, в то время как возможность проведения математически простого исследования зависимости аналитической модели от входных параметров (в первую очередь - вязкости, скорости мелкомасштабной конвекции или горизонтального сокращения плиты в обстановке латерального сжатия) является большим преимуществом.
Отправным пунктом выбора моделей являлся установленный по геологическим данным факт, что до начала активной стадии горообразования основные границы структурных этажей коры будущих Центрально-Азиатских орогенов были относительно плоскими и согласными, а также имели несколько повышенную мощность [Krestnikov et al., 1979]. Учитывая этот факт и определенную латеральную периодичность строения большей части структур (один цикл включает, соответственно, одно горное поднятие и соседнюю межгорную или крупную внутригорную впадину) в одном направлении и протяженность структур в другом (вдоль горных хребтов), было принято решение использовать двумерные латерально-периодические реологические дискретные модели, границы которых в начальный момент времени считались плоскими. Именно характер латеральной периодичности считался гармоническим.
Основными эмпирическими данными, с которыми соотносились результаты моделирования, являлись распределения ориентаций главных осей тензора напряжений, полученные для рассматриваемых регионов. Их источник - результаты текто-нофизических реконструкций, выполненных методом катакластического анализа разрывных смещений на основе данных о механизмах очагов землетрясений [Rebetsky et al., 2012, 2013, 2016] для Тянь-Шаня и Алтая. На рис. 1 приведены подобные данные для Алтайского горно-складчатого сооружения. Основной полученной закономерностью является преобладание субгоризонтальной ориентации осей максимального сжатия в коре хребтов и субгоризонтальной ориентации осей минимального сжатия (максимального девиаторного растяжения) в коре межгорных и крупных внутригорных впадин. Так как основным источником данных для реконструкций являлись относительно неглубоко-фокусные землетрясения (10-20 км), рассматриваемая закономерность наиболее актуальна для верхней коры и верхней части средней, т.е. до глубин 15-20 км.
2. ПОСТАНОВКА ЗАДАЧ. ОБЩЕЕ РЕШЕНИЕ
Модели литосферы рассчитывались в рамках приближения сплошной среды. Структура моделей в начальный момент времени для рассматриваемых задач приведена на рис. 2 и рис. 3. В модели I литосфера представлена двумя слоями (кора и мантия) с вязкостями, отличающимися на два порядка, а воздействие
Рис. 1. Результаты тектонофизической реконструкции для Алтая (по [Rebetsky et al, 2013]). Fig. 1. Tectonophysical reconstruction for the Altai area (after [Rebetsky et al., 2013]).
t M 0 Д ' ■ Дневная поверхность (подвижная граница) Ре. Лс Земная кора г
e ■ л ь - 1 Граница Мохоровичича (подвижная граница) Рт'Т,т Литосфера (мантия) г
Y X Подошва литосферы (неподвижная граница) ' Астеносфера
Рис. 2. Схема модели I (задача о влиянии астеносферной конвекции).
Fig. 2. Diagram of Model I (problem of the influence of asthenospheric convection).
Рис. 3. Схема модели II (задача о горизонтальном сокращении литосферы).
Fig. 3. Diagram of Model II (problem of the influence of horizontal shortening of the lithosphere).
со стороны астеносферы, в которой происходит термогравитационная конвекция, задается в виде скоростей течения мантийного вещества через подошву литосферы или вдоль нее. В модели II кора представлена двумя слоями, верхний из которых упругий (10-15 км), а нижний - вязкий. Кора лежит на менее вязкой мантии (разделение мантии на литосферную часть и астеносферную не проводилось).
Особо следует отметить, что реализация стремления системы кора - мантия к изостатическому равновесию обеспечивалась за счет выдерживания стандартного скачка плотности на подошве литосферы 0.4 г/см3 и использования вязкой модели мантии и всей коры или большей ее части. Отсутствие фактической границы литосфера-астеносфера (модель I) или скачка плотности порядка 0.02 г/см3 на ней (модель II) несущественно искажало амплитуды изостатических сил.
В текущем исследовании моделируются сверхмедленные процессы в литосфере (характерное время которых первые десятки миллионов лет), поэтому уравнения движения и связи напряжений и скоростей деформаций решались в линейной и квазистатической постановке (вместо уравнений движения выписывались уравнения равновесия). Определяющая система уравнений решалась методом разделения переменных Фурье. Реология основных тел моделей бралась в виде несжимаемого тела Ньютона (линейно-вязкая жидкость). Исключением являлась верхняя кора в плейт-тектонической модели (модель II), о которой подробнее будет сказано отдельно. Такой выбор реологии тел моделей обусловлен, с одной стороны, необходимостью использовать максимально простые реологические тела и с другой - тем фактом, что для большей части земной коры и, тем более, мантии, находящихся в закритическом и относительно разогретом состоянии на рассматриваемых временах все равно будут преобладать диссипативные процессы. Далее будет представлено общее решение краевой задачи, одинаковое для тел с ньютоновской реологией обеих моделей, после чего краевые задачи будут рассматриваться уже раздельно.
Итак, уравнения движения записываются в квазистатическом приближении, используется реология несжимаемого тела Ньютона, и задача решается в линейной и двумерной постановке [Landau, Lifshits, 1986]:
°ij,j = 0, (1) Оц= a8ij+T]{viij + vjii), (2)
pVi i = 0, i, j = x,y. (3)
Здесь и далее для удобства оператор дифференцирования по пространственным координатам записывается в индексе через запятую, как это принято в тензорном анализе. р - плотность, ^ - динамическая вязкость, 5 - символ Кронекера. Уравнение (1) - уравнение равновесия для дополнительного напряженного состояния, основным состоянием системы считается то, которое характеризуется наличием только лито-статического давления, т.е. нулевыми девиаторными напряжениями и нормальными напряжениями, равными:
= ~9У$1<1Р> 1 = х,у. (4)
Далее все выкладки и результаты приведены именно для дополнительного НДС системы. Под номером (2) записаны линейные соотношения Коши, под (3) - уравнение неразрывности. Поставленную задачу можно сформулировать в ином виде, например посредством замены дифференциальной формы записи соответствующих законов механики на интегральную, что часто используется при создании численных моделей.
Система (1)-(3) стандартными выкладками сводится к единому уравнению (Навье - Стокса), в данном случае принимающему следующий вид:
Чхки = 0, ¿,7, к = х, у, (5)
представляющему собой бигармоническое дифференциальное уравнение относительно любой из компонент вектора скорости смещения. Поясним, что, благодаря (3), уравнение (5) является набором двух эквивалентных уравнений, переход между которыми с легкостью осуществляется посредством (3). Решение (5) будет находится методом Фурье, применение которого базируется на вышеупомянутой латеральной гармонической периодичности моделей. Введем пространственную частоту к = 2п/А, где А - длина волны системы (латеральный размер цикла хребет - впадина). Гармоничность периодичности означает возможность представления базовых механических функций системы (компоненты вектора скорости, тензора деформации и напряжений) в виде:
Ъу(х,у) = ру(у)соз(кх), (6)
для прочих величин - аналогично. Подставляя (6) в (5), получим биквадратное обыкновенное дифференциальное уравнение 2-й степени относительно Уууу\
^уууууу + = ° í=x,y,
решение которого дает следующие выражения для компонент вектора скорости:
рх = С1скх + С2зкх + С3 хскх + ру = (С3 ~С2)скх + (С4 ~С1)зкх - С^кх - С4хскх, (7)
где С; - произвольные постоянные, а х = ку. Сходные выражения для компонент тензора напряжений сразу могут быть получены из (2) и (7):
дхх = 2г]к(С1скх + С^кх + С3(25Л* + *сЛ*)+С4(2сЛ* +
аху = 2г]к(С2скх + С1Бкх + САхс^х + С3^х) , (8)
дуу = 2г]к(-С1скх - С^кх + С3(25^ - хсЬх)+СА{2скх - Х^Ьх)), д = 2г]к(С3зкх + С4скх).
На основе (8) вычисляются требуемые для тектонофизической интерпретации инварианты тензора напряжений - главные значения (о^), максимальное касательное напряжение (ттах) и углы ориентации главных осей от горизонта (аг и а2).
3. Задача о влиянии мелкомасштабной астеносферной конвекции. Эволюционный метод получения решения
После получения общего решения (7) и (8) требуется найти удовлетворяющие краевым условиям частные решения поставленных краевых задач. При том что задача решается в двумерной и линейной постановке, деформация границ модели и структура формирующегося НДС характеризуются достаточно сложным и нелинейным поведением во времени даже с учетом латеральной гармонической симметрии модели, что требует поиска специальных и нестандартных подходов для получения аналитического решения. В рамках текущего исследования будет применяться подход, который будет далее называться эволюционным и структура которого будет детально раскрыта в применении к рассматриваемой в этом разделе задаче о формировании НДС литосферы под влиянием мелкомасштабной термогравитационной конвекции в астеносфере (модель I). В плейт-тектонической задаче (модель II) будет использоваться тот же подход, поэтому последовательность построения решения уже не будет так подробно поясняться.
Модель литосферы, используемая в задаче о влиянии конвекции, уже была выше кратко обрисована и представлена на рис. 2. Она состоит из двух однородных линейно-вязких тел, верхнее из которых соответствует земной коре, а нижнее - мантийной литосфере (далее, учитывая структуру модели, будем применять термин «двуслой»). Сама астеносфера напрямую в модели не участвует, в рамках текущего исследования было достаточно учесть только ее влияние на литосферу. Таким образом, здесь представлено именно решение задачи о вынужденном течении через подошву литосферы, а не сама задача об астеносферной конвекции. В обзорной работе [Dobretsov et al., 2013] приведены характерные скорости конвективного процесса, составляющие первые см/год.
В нашей модели влияние астеносферной конвекции на литосферу задается в виде кинематического краевого условия на подошве двуслоя. В статье [Myagkov, Rebetsky/ 2016] уже рассматривалась подобная задача, однако в качестве соответствующего кинематического условия задавалось распределение горизонтальной компоненты вектора скорости смещения. Такой вид краевых условий в той работе был использован для сопоставления получаемого решения с результатами, представленными в работе [Mikhailov et al., 1999; Timoshkina et al., 2010]. Здесь же будет рассмотрено действие вертикальной компоненты течения, связанное, очевидно, с сообщением подошвенной части литосферы направленного к дневной поверхности импульса в области восходящей части конвективного потока и обратно направленного в области нисходящего. Отметим, что в работах [Bobrov, Trubitsyn, 2003; Bobrov, Baranov, 2011] показано, что
--------— ......*......
м о д___ ♦ * Земная кора
е л ь тг * Литосфера (мантия)
W Ж
Астеносфера
Рис. 4. Схема распределения возмущений в модели I. Тонкими стрелками в астеносфере показано направление течения, связанного с конвекцией, толстыми стрелками в литосфере изображены действующие на границы модели дополнительные гравитационные силы.
IFig. 4. Disturbance pattern in Model I. Thin arrows - direction of the convective flow in the asthenosphere. Thick arrows -additional gravitational forces in the lithosphere at the boundary of Model I.
конвекция в мантии способна создать на подошве литосферы дополнительное давление в области восходящего потока до 40 МПа. В данной работе мы рассматриваем именно влияние вертикальной компоненты скорости течения конвективного процесса на подошву литосферы.
Схема задачи влияния астеносферной конвекции (базового возмущения) представлена на рис. 4. Расчет в модели I, так же как и в модели II, ведется в эйлеровых координатах, поэтому нижнюю границу модели, где задается кинематическое граничное условие (ГУ), нет нужды задавать подвижной, она будет оставаться плоской на всем протяжении формирования НДС системы. Само кинематическое условие запишем в виде:
г£озм(х,у) = 0, ^озм(х,у) = v0cos(kx). (9)
Теперь опишем структуру эволюционного метода построения решения. Наибольшую сложность при этом представляет необходимость записывать ГУ на искривленных деформациями границах модели (подчеркнем, что, в отличие от нижней, прочие границы деформируемы), а именно - на дневной поверхности и границе Мохо. Поскольку при задаваемой мощности коры 40 км и характерной длине цикла хребет-впадина для орогенов Центральной Азии более 100 км (в данном исследовании будет взята величина 120 км) амплитуда формирующегося рельефа на подвижных границах не будет превышать первые километры, имеется возможность использовать подход Лява [Love, 1911]. Этот подход подразумевает необходимость записывать ГУ на плоских границах, а эффект от их искривления - учитывать в виде динамической поправки к ГУ. Однако остается еще открытым вопрос о самом характере искривления границ, без знания которого нельзя ввести корректно упомянутую поправку к ГУ (поскольку речь идет о сложном процессе деформирования границ во времени, знания современной амплитуды рельефа явно недостаточно для введения такой поправки). Необходимо заранее вычислить функции изменения величины амплитуды рельефа во времени (далее, для краткости, они будут называться «эволюционные кривые»), на основе которых уже получить окончательное нестационарное решение краевой задачи.
Методика получения эволюционных кривых базируется на том, что, в силу линейной постановки общей задачи, результат влияния нескольких факторов на общее НДС системы представим в виде суперпозиции НДС, формирующихся за счет отдельных факторов. В нашем случае имеются три основных фактора: базовое возмущение, влияние искривления внутренней границы модели и поверхности, что можно записать как
<щее = ^ + о* + = xy, (10)
где верхние индексы указывают на принадлежность текущей функции состояния системы к НДС, формирующемуся за счет перечисленных выше факторов (отдельно). Так как оператор пространственного диф-
ференцирования, также в силу линейности постановки, не нарушает подобную (9) закономерность, из (2) следует, что и для полей вектора скорости смещения, соответствующих тем же НДС, формирующихся от отдельных факторов, можно записать:
Далее, несложно установить характер влияния, оказываемого на систему искривлением отдельной границы. В силу исходной симметрии базового возмущения, границы также будут искривляться по закону (9). Плотность нижележащих слоев модели выше, поэтому на области близ положения изначально плоской границы будут действовать положительные и отрицательные (архимедовы) гравитационные силы, где формируются соответственно положительные и отрицательные формы рельефа (рис. 4). В рамках подхода Лява добавим в записанный для нормальных напряжений ГУ дополнительный член:
где Ар - разность плотности между разграниченными телами, g - ускорение свободного падения, ^ - амплитуда рельефа.
Теперь, в соответствии с намеченным подходом, получим решения краевых задач, для которых НДС формируется за счет действия только одного из намеченных факторов: базового возмущения, искривления внутренней границы модели и искривления поверхности модели. Эти подзадачи рассматриваются в стационарном приближении, нестационарным предполагается уже финальное решение. ГУ для первой подзадачи, связанной с базовым возмущением, имеет, с учетом реологии тел модели и (9), вид:
рх = 0 ру = р0 (нижняя граница)
[рх] = 0 [ру] = 0 [&Ху\ = 0 №уу] = 0 (внутренняя граница)
аху = 0 дуу = 0 (верхняя граница).
ГУ для второй подзадачи (возмущение - искривление внутренней границы) вытекает из (12), здесь, очевидно, Ар = рт— рс. Имеем:
рх = 0 ру = 0 (нижняя граница)
[рх] = 0 [ру] = 0 [дху\ = 0 [ауу] = (рт - рс)д(т (внутренняя граница)
аху = 0 ®уу= 0 (верхняя граница).
Аналогично, ГУ для третьей, связанной с искривлением поверхности, подзадачи:
рх = 0 ру = 0 (нижняя граница)
[рх] = 0 [ру] = 0 [&Ху\ = 0 [&уу\ = 0 (внутренняя граница)
дху = 0 &уу= Рсд^с (верхняя граница).
Для выписанных ГУ гармоническая х-составляющая уже опущена, для определенности здесь и далее будем считать, что оставшиеся ^-составляющие относятся к центру нисходящего конвективного потока по оси x. Подставляя (7), (8) в ГУ, получим систему из восьми алгебраических уравнений для определения восьми констант С]т и С]с, ¿=1.3, у'=1.4 (/' - индекс подзадачи). Тривиальное, но весьма громоздкое решение этих систем для подзадач здесь не приводится.
Полученные решения зависят от параметров первый из которых (скорость астеносферной
конвекции) является входным параметром модели и считается постоянным (обратная связь между вынужденной вторичной литосферной конвекцией и астеносферной здесь, естественно, не учитывается). Для рельефа подвижных границ известно лишь начальное значение (ноль), характер его изменения как раз требуется выяснить - это и будут искомые эволюционные кривые. Для этого рассмотрим, для определенности, поверхность модели и запишем для нее соотношение (11) (будем добавлять нижний индекс «с» для обозначения принадлежности к этой границе и, далее, «т» - для внутренней):
общее
= pf + pf + vf, i = x,y.
(11)
/доп = kpgt при <= < cos(fcx),
(12)
общее
Vе
Для поля стуу, полученного из решения третьей подзадачи, имеет место, как зависящего от параметра:
0уу — оуу/дош или же, с учетом (12):
<7уу — ОууРсд^с,
откуда, аналогично переходу от (10) к (11), следует
^у — %Рс9$с . (14)
Но, так как
9 ^ ^общее Ьс — ■'ус ,
уравнение (13), учитывая (14), можно переписать в следующем виде:
(с = Яу^щ + р$сАрд(т + у^Рсд^. (15а)
В последнем выражении было введено обозначение £у, аналогичное прочим, произведены полностью аналогичные преобразования члена ^у!с и для краткости введено обозначение Ар = рт— рс. Аналогично получим для внутренней границы:
<Гш = + $утАрд$т + $утрсд(с. (15б)
Решение системы (15), дополненное вышеуказанными начальными условиями (плоские границы в начальный момент времени), дает выражения для амплитуд рельефа в явном виде. Это и есть искомые эволюционные кривые:
& = С0с + Ард г)у3с(С1сехр(у1£) + С2свхр(у21)) ^т С0т ^ (/1 - Рс^с^т^РМ + (72 - Рс^с^т^Р^),
(16)
где /1(2) = 0.5д( ЩтАр + ЩсРс ± ^($утАр + х>^срОг -4рсАр{ д^ - Щт)), а выражения для констант Ст и Сс опущены из-за громоздкости.
На рис. 5 изображены эволюционные кривые, рассчитанные для значений параметров, которые будут использованы и в дальнейшем: Ьс—40 км, ^т—60 км, рс—2.7 г/см3, рт—3.1 г/см3, ^с—1023 Па-с, ^т—1021 Па-с (далее - базовый набор параметров). Структура деформации модели, соответствующая этим кривым, схематически изображена на рис. 6. Полученные эволюционные кривые полностью определяют ход деформации модели и требуют отдельного пояснения.
Отметим, что приведенные далее данные о времени изменения геометрии коры, так же как и скорости внутрилитосферного течения и напряжения, существенно зависят от принятых выше значений вязкости коры и мантии. Однако если уменьшить вязкость, сохраняя соотношение между корой и мантией в два порядка, то напряжения остаются неизменными и при этом пропорционально увеличиваются скорости и уменьшается время формирования характерных структурных форм коры.
Будем, как и раньше, рассматривать формирование рельефа над центром нисходящего конвективного тока. Граница Мохо здесь монотонно прогибается конформно направлению астеносферного течения. Дневная поверхность первое время проявляет сходную динамику, с меньшей, конечно, амплитудой прогиба. В какой-то момент (стадия 2 на рис. 6, максимум см. на рис. 5) этот прогиб достигает максимума. Во всей соответствующей области модели вертикальная компонента скорости течения положительна. Однако далее силы, связанные с гравитационной нескомпенсированностью, приводят к обращению процесса, стремясь достигнуть изостазии. Система сначала полностью нивелирует сформированный ранее прогиб
Рис. 5. Функции эволюции амплитуды рельефа дневной поверхности (синий цвет, построена относительно левой вертикальной шкалы) и подошвы коры (красный цвет, построена относительно правой вертикальной шкалы) для модели I. Пунктирной линией обозначены кривые, построенные без учета влияния экзогенных процессов. Для кривых, обозначенных сплошной линией, Л=0.1
Fig. 5. Model I. Evolution of the amplitudes of the ground surface relief (blue - function built relative to the left vertical scale) and the crust bottom (red -function built relative to the right vertical scale). Dashed lines - curves plotted without taking into account the effect of exogenous processes. A=0.1 for the curves shown by solid lines.
Земная кора
Мантия
Земная кора
Мантия
Земная кора
Мантия
—- Земная кора Мантия
Земная кора Мантия
I Рис. 6. Основные этапы деформации моде ли I.
Fig. 6. Model I. Main stages of deformation.
дневной поверхности (рис. 6, стадия 3), а далее начинает формироваться инверсионная форма рельефа (рис. 6 - стадии 4 и 5), амплитуда которой в некоторый момент (около 8-10 млн лет) выходит на асимптоту. В результате формируется классическая геодинамическая пара - корни под горными поднятиями и антикорни под впадинами.
Отдельно стоит заметить, что даже ко времени выхода на асимптоту эволюционных кривых система не достигает полной изостазии, что объясняется незамкнутостью используемой модели: астеносфера явно не входит в ее состав, а ее влияние перманентно присутствует в составе соответствующего ГУ. Вследствие этого нисходящий астеносферный ток формирует дополнительную положительную добавку к полю вертикальной компоненты течения в области инверсионного поднятия, которое в результате не достигает амплитуды, соответствующей компенсации. При данных параметрах «коэффициент изостатической нескомпенсированности» (отношение амплитуд рельефа на границах, нормированное на то, которое было бы при полной скомпенсированности) примерно равен 3. На рис. 7 приведен результат решения отдельной задачи об эволюции этой же модели при «выключении» конвекции (в качестве начального состояния - асимптотическое состояние для модели с включенной конвекцией). Слева -эволюционные кривые, а справа - значение коэффициента изостатической нескомпенсированности. Видно, что система быстро достигает изостазии, после чего крайне медленно (сотни млн лет) релаксирует. Только отсюда уже может быть сделан вывод, что при объяснении влиянием астеносферной конвекции эпиплатформенного орогенеза для изостатически скомпенсированных структур разумно предположить, что интенсивность астеносферной конвекции заметно ниже таковой при активной стадии горообразования, что может быть связано, к примеру, с деградацией конвективного процесса или смещением конвективной структуры в иные части астеносферы по латерали.
Подставляя полученные выражения (16) в решения 2-й и 3-й стационарных подзадач, мы получим окончательное нестационарное решение исходной краевой задачи. Поскольку
-400-
-800-
-1200-
0.8-1 1.21.62.32.8-
3.2
0
100
млн лет 200
Рис. 7. Функции амплитуды рельефа (слева) и график коэффициента изостатической нескомпенсированности (справа) (определение в тексте) для модели I после прекращения воздействия астеносферной конвекции. Точка отсчета соответствует 50 млн лет рис. 5.
Fig. 7. Model I. Functions of the relief amplitude (left), and the curve of isostatic non-compensation coefficient (right) (see the definition in the text) after termination of asthenospheric convection effects. Calculation reference point: 50 Ma. See Fig. 5.
при этом НДС системы меняется настолько медленно, что позволило использовать уравнения равновесия (1) вместо уравнений движения, в каждый момент такое решение близко к стационарному, поэтому для него корректно использовать термин «квазистационарное». Однако перед непосредственным получением решения рассмотрим еще один важный фактор.
Экзогенные процессы. Важным фактором, влияющим на формирующееся НДС системы, является наличие процессов эрозии и денудации осадочного материала из областей поднятий и последующей аккумуляции в депрессиях, протекающих на поверхности земной коры [Ollier, 1981; Makarov et al., 2010, 2011]. Эрозионно-аккумуляционные процессы не обязательно сводятся к релаксационным [Rebetsky, 2012; Rebetsky, Pogorelov, 2013]. При наличии активно формирующегося рельефа они также выполняют активную функцию, замедляя, конечно, рост амплитуды рельефа со временем, однако повышая скорости движения материала в геодинамической системе, особенно близ поверхности, что, в свою очередь, приводит и к росту напряжений в среде. Отсюда вытекает необходимость учета экзогенных процессов в рамках используемой модели.
Вопрос о способе математической формализации данных процессов на текущий момент достаточно хорошо изучен как с позиций согласованности с геологическими данными [Culling, 1960], так и с физико-теоретических позиций [Blatt et al., 1972; Ahnert, 1970]. Существует большое количество численных моделей связи горообразования и экзогенных процессов [Thieulot et al., 2014; García-Castellanos et al., 2012]. Наиболее часто применяемая при решении геодинамических задач, подобных рассматриваемой в данной статье, математическая модель аккумуляционно-денудационных процессов предполагает зависимость между скоростью изменения амплитуды рельефа за счет сноса/накопления осадков и формой рельефа в виде:
(=ЛД(, (17)
где £ - амплитуда рельефа, Д - лапласиан, Л - числовой коэффициент, характеризующий интенсивность процесса денудации поднятий и осадконакопления во впадинах. Коэффициент Л имеет обратную ко времени размерность (обычно исчисляется в [ед/год], а дифференцирование ведется по латеральным пространственным координатам. Подобная модель учета поверхностных процессов успешно применялась, например, в работах [Mikhailov, 1983; Horowitz, 1976].
Специально отметим, что форма уравнения (17) определяет именно денудацию горных поднятий, а не эрозию их склонов. Форма уравнений удобна тем, что она позволяет задать столь разные поверхностные процессы в виде простой гармонической функции. На самом деле в исследуемых нами центрально-
азиатских орогенах на вершинах хребтов сохранился пенеплен, что говорит о необходимости использования модели процесса в виде эрозии склонов. Однако данная модель сложно реализуема в рамках используемого нами аналитического подхода.
В рамках текущей задачи применение подхода, базирующегося на уравнении (17), крайне удобно, так как дифференцирование ведется только по х-координате, и, в силу гармонической латеральной симметрии модели, оно эквивалентно домножению дифференцируемой функции (в данном случае £с) на множитель -к2. Учет влияния экзогенных процессов на амплитуду рельефа должен быть произведен, очевидно, в виде добавления дополнительного члена ЯА^с в эволюционном уравнении дневной поверхности (15а), который, соответственно, будет иметь вид —к2Л^с. Тогда вместо (15) мы получим систему:
поэтому решение (17) будет даваться также формулами (16), в которых произведена замена (18).
В разделе «Обсуждение результатов» данные, рассчитанные без учета влияния экзогенных процессов, отдельно не приводятся. Укажем только, что при использовании базового набора параметров экзогенного процесса увеличение абсолютной величины вектора скорости в среднем в коре за счет этого фактора происходит более чем в пять раз.
Важно отметить, что учет экзогенных процессов приводит к размыканию конвективных ячеек близ поверхности, формируемому на постинверсионной стадии эволюции модели (стадии 4 и 5 на рис. 10). Естественно ожидать сильное влияние денудации на амплитуду формирующегося рельефа, для сравнения, на рис. 5 изображены (пунктиром) эволюционные кривые, рассчитанные без учета влияния экзогенных процессов. Помимо двукратного, на начальной стадии, и более чем четырехкратного, на асимптотической, снижения рельефа, денудационные процессы приводят также и к изменению темпа эволюции модели, а именно - несколько ускоряют достижение системой основных структурных точек эволюции (так, выравнивание доинверсионных форм рельефа происходит за 2.9 вместо 5.6 млн лет). В общем, характер влияния экзогенных процессов математически достаточно сложный и существенно нелинейный, поэтому его детальное исследование требует отдельного исследования.
Итак, окончательное решение поставленной задачи, которое учитывает влияние экзогенных процессов, зависит от эволюционных кривых, в выражении (16) для которых необходимо сделать замену (18). Само решение строится следующим образом: в любой момент времени мы имеем значения амплитуды рельефа и £с. Для этих значений вычисляется решение соответствующих подзадач, описанных в предыдущем разделе. Далее строится суперпозиция (10) или (11) для функций НДС системы. Эти функции (компоненты вектора скорости смещения и тензора напряжений) составляют искомое решение. Результаты расчетов представлены на рис. 10, их обсуждение ведется в разделе «Обсуждение результаты».
4. Задача о формировании НДС литосферы в обстановке горизонтального сжатия
Как уже было сказано, в рамках текущей статьи вязкая реология аппроксимирует реально упругопла-стическое деформирование литосферы. При этом верхняя кора [М1ко!аеузку, 1983] находится в закритиче-ском состоянии только в зонах разломов, разделяющих ее на блоки, которые могут быть представлены упругим телом. В предыдущей задаче (модель I) земная кора была представлена единым слоем, в результате чего корректнее было рассматривать ее в виде вязкого тела (так как верхняя кора составляет относительно небольшую (меньше четверти) часть всей коры). В текущей задаче рассматривается НДС, формирующееся в коре (области формирования эпиплатформенных орогенов Центральной Азии) в обстановке горизонтального сжатия, как результат влияния коллизии плит. Во внутренней части платформы, вдали от зоны коллизии, основная часть силового воздействия передается именно через упругую верхнюю часть коры, поэтому в модели II (см. рис. 3) она фигурирует отдельно.
Структура решения текущей задачи строится сходным с предыдущей образом, поэтому здесь мы только кратко обрисуем последовательность построения решения. Деформации верхней коры будут рассмат-
(18)
которая, в общем-то, формально совпадает с исходной при проведении замены:
С — (к2А/Ард),
(19)
риваться в рамках балочной теории. Считается, что за счет давления со стороны соседней плиты происходит потеря устойчивости балки при превышении эйлеровой (первой критической) силы. Для оценки значения этой силы с учетом влияния вязкой коры и мантии можно воспользоваться решением, приведенным в монографии [Schubert et al., 2001].
Исследуя эпиплатформенный орогенез, будем считать, что потеря устойчивости балки происходит с той же длиной волны, что и в предыдущей задаче (равной размеру цикла хребет - впадина). Образующееся НДС самой верхней коры будет рассмотрено в следующем разделе, здесь же покажем, как с помощью эволюционного подхода получается решение для прочей части модели. Поскольку, как будет видно далее, основная часть энергии концентрируется в коре, текущая модель не содержит мантийной литосферы в качестве отдельного тела. Ниже коры в модель включается аппроксимирующее мантию полупространство.
В отличие от предыдущей задачи (модель I), базовым возмущением модели II является движение потерявшей устойчивость верхней коры. Положим, что кора сокращается с постоянной скоростью деформации ¿хх. Из элементарных геометрических соотношений следует:
= ßL2l&l (20)
где - амплитуда вертикальной скорости движения осевой поверхности балки, а ß для условий рассматриваемой геомеханической задачи в приложении к Алтаю может быть принято равным 15, что дает скорости порядка 1-2 мм/год. В рамках балочного приближения корректно считать, что это и скорость верхней поверхности вязкой части модели. Условие (20) полностью аналогично (9) для предыдущей задачи. Однако теперь базовое возмущение задается на поверхности модели, а не на подошве. Не дублируя аналогичные предварительные рассуждения, выпишем ГУ для соответствующих подзадач. Для основной подзадачи, связанной с давлением балки:
vy = öxy = 0 (верхняя граница)
[vx] = 0 [vy] = 0 [<7ху] = 0 [дуу] = 0 (внутренняя граница)
vx = 0 vy= 0 (на бесконечности).
В качестве второго условия на верхней границе было выбрано условие отсутствия проскальзывания между верхней и средней корой. Для задачи, связанной с искривлением границы Мохо:
vy = 0 дху = 0 (верхняя граница)
[vx] = 0 [vy] = 0 [аху] = 0 [öyy] = (рт - pc)g^m (внутренняя граница)
vx = 0 vy= 0 (на бесконечности).
Так как литосферная мантия здесь не вводится как отдельная часть модели, мы решаем только две подзадачи. Однако эволюционных кривых будет при этом по-прежнему две, в качестве первой будет выступать кривая для рельефа поверхности коры, определяемая тем же соотношением (20). По аналогии с предыдущей задачей, выпишем уравнение для определения второй эволюционной кривой (граница Мохо):
<fm = ^ут^у + Цт^Рд^т . (21)
Также рельеф в начальный момент времени полагается плоским. Решение (21) запишем в виде:
V1 V0
Ст = Ст(ехр( Щт^рдЬ) ~ 1) при Ст = JT^g- (22)
Следует, однако, учесть и тривиальное напряженное состояние, формирующееся чисто за счет сокращения самой вязкой части модели (до этого анализировался только вклад давления упругой балки на вязкую часть коры). Добавка к напряженному состоянию вязкой средней и нижней коры будет состоять только из горизонтальной нормальной компоненты тензора напряжений и равна:
о^ = (23)
Однако важнее не только учесть (23), но и тот факт, что горизонтальное сокращение вязкой части коры приводит к ее утолщению (упругая часть коры реагирует потерей устойчивости). При постоянной скорости деформации толщина вязкой части коры будет утолщаться следующим образом:
I Рис. 8. Функции эволюции амплитуды рельефа границы Конрада (синий цвет) и подошвы коры (красный цвет) для модели II.
IFig. 8. Model II. Evolution of the Conrad boundary relief amplitude (blue) and the crust bottom (red).
кс = йс0ехр(|£ххЮ,
откуда, дифференцируя по времени, получаем добавку к скорости:
£уоп = кхх^соехрО^ххю,
которую необходимо учесть в виде дополнительного члена в (20). Так как этот член - функциональный, он несколько усложняет получение решения. Опустив для краткости выкладки, выпишем лишь конечный результат:
Çm
'Ст ^
1£хх|"с0
Рут1-
expG^XXlO + (С г
|ёХХ|^со
дут1-
-)ехр( tfmApgt).
(24)
Земная кора
Мантия
Земная кора
Кривые, определяемые выражениями (20) и (24), представлены на рис. 8. В качестве значения ¿хх была взята величина 10-9 ед/год.
Существенным моментом, определяющим различие между структурой эволюции в текущей и предыдущей
- задачах, является отсутствие инверсионных процессов
для рельефа, формирующегося на границах модели (рис. 9). Этот результат никак не связан с относительной простотой модели, необходимой для проведения аналитического расчета. По изложенной в статье методике были рассчитаны и более сложные модели, включающие, к примеру, среднюю кору как отдельное тело. Однако и они дали сходные результаты. Кроме того, существуют и численные модели, согласные с рассматриваемой здесь. Например, в статье [Cloetingh et al., 2002] получены сходные результаты для коры Иберии. Таким образом, следует считать, что НДС литосферы, формирующееся в результате потери устойчивости упругой части коры в обстановке горизонтального сжатия, не приводит к формированию корней/антикорней в горно-складчатой системе, что одно уже говорит против использования этого процесса для объяснения эпиплатформенного орогенеза.
Балочное решение для упругой коры. Для решения общего вопроса о структуре НДС, возникающего в обстановке горизонтального сжатия литосферы, достаточно было рассмотреть эволюцию системы, литосфера которой испытывала горизонтальное сокращение с постоянной скоростью деформации ¿х^. В этих условиях влияние упругой верхней коры (балки) на вязкую часть модели могло быть учтено в виде кине-
Мантия
Земная кора
Мантия
Рис. 9. Основные этапы деформации модели II. Fig. 9. Model II. Main stages of deformation.
матического ГУ, определяемого формулой (20); никакого дополнительного знания о процессах, происходящих в самой балке, не требовалось. Для определения же НДС упругой коры, представляющего самостоятельный интерес, необходимо учитывать взаимовлияние упругой и вязкой части модели друг на друга. Ниже кратко представлена последовательность построения подобного решения.
Во-первых, определим величину первой критической силы, достижение которой необходимо для потери устойчивости балки. Она рассчитывается по формуле [Turcotte, Schubert, 1982]:
Рэ = (25)
где V - коэффициент Пуассона балки, Е - модуль Юнга. Удобнее, правда, пересчитать Рэ в терминах напряжений. Соответствующее критическое напряжение получается в результате нормировки Рэ в (25) на мощность верхней коры hcl. Полученное значение приблизительно равно 30 МПа и, очевидно, достигается.
Теперь получим непосредственное решение для напряжений в балке. Для этого необходимо учесть все действующие на балку силы, это:
1) сила продольного сжатия P (нормированная на площадь сечения балки), соответствующая деформации ^ ;
2) силы реакции со стороны вязкого слоя qr(x) = ауу (на границе вязкая/упругая кора);
3) массовые силы гравитационной природы, аналогичные рассмотренным во вспомогательных подзадачах для основных задач ранее, qm(x) = рсд^с, где £с - амплитуда рельефа поверхности вязкой части модели, которая в балочном приближении равна смещению осевой поверхности балки £с = и®.
Соответственно, уравнение изгиба упругой балки [Rabotnov, 1988] запишем в виде:
+ PhcluliXiX + qr+ qm = 0. (26)
Решение (26) находится посредством разделения переменных, при постоянной скорости деформации ¿х^ решение для продольной силы будет следующим:
Р = ((Е1хк4 + рсд)и°у + асуу)/ hclk2u°y, (27)
где определено (20) и равно Î7°t. На основе (27) рассчитываются компоненты тензора напряжений, здесь мы приведем только конечные формулы:
&бх = Ek2UyW + Р,
обу = 2- EIxk*u0y(1 - QVci, (28)
öly = l EIxk3Uy(w — + \hcl+ acyy),
в которые была введена вспомогательная переменная w, соответствующая измерению по вертикали отдельно в балке. Результаты расчета по формуле (28) строились для того же сечения, что и эволюционные кривые ранее. Приведем результаты расчета по формулам (27) и (28) для значения Е=5 ГПа и амплитуды рельефа =0.5 км, считая Ьс1=10 км. Получим значение отношения продольной силы к площади сечения балки Р=29 ГПа, тогда как связанные с изгибом (без учета Р) напряжения колеблются в пределах от 340 до -340 МПа. Это означает, что горизонтальные нормальные напряжения не меняют свой знак в вертикальном сечении балки. Отметим, однако, что в геологии существуют ошибочные представления, что изгиб упругой литосферы должен в соответствующих областях сопровождаться растяжением (верхняя часть балки в антиклинали и нижняя часть балки в синклинали). Это совершенно неверно в том случае, если вся литосфера в целом находится под действием продольных сил сжатия. Даже там, где изгиб приведет к снижению величины а , оно все же останется того же знака, что и в остальных частях системы, т.е. сжатием. Таким образом, верхняя кора целиком находится в обстановке горизонтального сжатия. Кроме того, отметим необходимость достижения колоссального уровня горизонтальных напряжений (десятки ГПа) для реализации вышеописанной потери устойчивости верхней коры, что означает, что в действительности она в подобной геодинамической обстановке достаточно быстро деструктирована посредством формирования шарьяжей и надвигов, что само по себе говорит против соответствия модели II природным данным.
5. Обсуждение результатов
Результаты моделирования (рис. 10, 11) будут показаны для значений входных параметров, соответствующих базовому набору (представлен выше). Отметим, что большинство параметров определены: для континентальной коры известна как плотность, так и мощность основных структурных этажей с высокой точностью, период возмущения по латерали определен из географических данных, хотя и несколько варьируется. Амплитуды возмущений и вязкости подбирались так, чтобы рассчитываемые девиаторные напряжения в среднем соответствовали по уровню имеющимся в среде, а время эволюции модели было подходящим с геологических позиций (первые млн лет/десятки млн лет).
Приступая непосредственно к анализу, отметим, что наиболее важными является структура деформации модели, включающая в себя эволюцию формирования рельефа и структуру формирования поля течения (конвективных ячеек) в модели, и структура формирующегося напряженного состояния модели, т.е. ориентации главных осей. Сами напряжения здесь не так важны, так как их детальное распределение в коре, в общем, неизвестно (напомним - здесь не идет речь о литостатическом напряжении), а выполненные в работах [ЯвЪвЬзку вЬ а!., 2012, 2013, 2016] оценки величин природных напряжений в большей степени являются априорным ориентиром для подбора адекватных значений амплитуд базовых возмущений. Все же для репрезентативности на рис. 10 и рис. 11 во второй колонке ориентации указаны поверх поля максимального касательного напряжения.
Первая модель (о влиянии астеносферной конвекции) продемонстрировала многоэтапную структуру эволюции НДС. На первом этапе (см. рис. 10, строка 1 и 2), продолжающемся до 12-15 млн лет, область литосферы, находящаяся над нисходящим астеносферным конвективным током (далее для определенности будем рассматривать именно эту область, для литосферы над восходящим током верны все те же рассуждения с поправкой на смену знака рельефа) испытывает погружение и на поверхности рельефа формируется прогиб, амплитуда которого достигает максимума в 8.7 млн лет (см. рис. 10, строка 2), после чего начинается медленное воздымание изначального прогиба. Максимальная амплитуда погружения составляет 700 м, уровень напряжений в коре достигает 100 МПа на подошве, но в среднем составляет 40-50 МПа. Ось ст3 направлена субгоризонтально, что прямо противоположно требуемому. Далее постепенно воздымание приводит к полному выравниванию форм рельефа (см. рис. 10, строка 3) на времени 60 млн лет (вид-
но, что скорость эволюции системы снижается со временем, см. рис. 5). При этом сглаживаются абсолютно все особенности в напряженном состоянии: к примеру, поле Ттах не меняется по латерали (одномерно). Уровень девиаторных напряжений в коре сильно снижается. При этом ориентации главных осей пока сохраняют прежнюю структуру. Однако далее воздымание не сменяется на погружение, оно продолжается (см. рис. 10, строка 4 и 5).
Переходим к постинверсионной стадии эволюции модели, которая характеризуется тем, что, во-первых, в коре формируется замкнутая конвективная ячейка (объединяющая области поднятия и впадины), а ориентации ст3 в верхней части коры сменяют направление с субгоризонтального на субвертикальное. При этом на поверхности над центром конвективной ячейки формируются локальные максимумы напряженного состояния (около 15 МПа для гтах). После инверсии на исследуемом участке коры исходное поле ориентаций соответствует уже воздыманию, что совпадает с нашими представлениями. Однако формирующаяся постинверсионная смена структуры ориентаций на обратную вновь приводит к несоответствию с данными реконструкций. Следовательно, именно начало постинверсионной стадии эволюции, когда область инверсии режима напряженного состояния в верхней коре достаточно маломощная (до 5-10 км), будет наиболее удовлетворять природным данным.
Вторая модель показала несоответствие природным данным еще на стадии получения эволюционных кривых: отсутствие инверсии (в разделе, где давалось построение решения, этот вопрос детальнее освещался) является существенным недостатком с точки зрения соответствия природным данным. Тем не менее рассмотрим эволюцию напряженного состояния, полученного для данной модели. Структура течения здесь не представляет особого интереса, на рис. 11 она достаточно явно показана и не требует отдельных пояснений. Напряженное же состояние имеет достаточно сложную структуру, связанную с тем, что составляющая (23) нарушает исходную латеральную симметрию НДС. Результаты моделирования представлены для времени 0.3 млн лет, соответствующего моменту формирования рельефа с перепадом высот 1 км. В области прогиба мы наблюдаем преобладание субвертикальной ориентации оси ст3 на глубинах ниже 15 км. Однако в верхней части модели под погружающейся верхней корой, под впадинами, все же преобладают субгоризонтальные ориентации оси ст3, как и в области поднятий. Отметим, что речь идет сейчас только о средней и нижней коре, для верхней коры было использовано балочное решение (подробнее см. в соответствующем разделе).
Рис. 10. Сводная таблица полученных данных (для указанного в тексте базового набора параметров) для модели I. В левой колонке - абсолютные значения вектора скорости течения с нанесенными поверх ориентация-ми вектора скорости. Абсолютные значения скорости в коре увеличены в 100 раз. В правой колонке - максимальные касательные напряжения, с нанесенными поверх ориентациями оси стз. Данные приведены последовательно (сверху вниз) для значений времени 0; 8.7; 61.5; 93.0; 150 млн лет.
Fig. 10. Summary table for Model I. The base set of parameters is specified in the text. Left column - absolute values of the flow velocity vector; orientations of the velocity vector are marked above. The absolute velocity values in the crust are increased 100 times. Right column - maximum tangential stress values; orientations of axis аз are marked above. The data are given consecutively (top to bottom) for time values of 0; 8.7; 61.5; 93.0; and 150 Ma.
0.5 km 0-0.1-
k -0.2-] о
p -0.3-a
0.7 км -0.5-
M -0.6-a
» -0.7-
И-0.8-
-0.9-
1 Г
0..5.KM 0
мм/год
хЮОкм -0.2
0.2
0.4
0.6
МПа
—г 280 -260 -240 -220 -200 -180 -160 -160 -120 -100 -80 -60 -40 .-20
хЮОкм -0.2
Рис. 11. Сводная таблица полученных данных (для указанного в тексте базового набора параметров) для модели II. Слева - абсолютные значения вектора скорости течения с нанесенными поверх ориентациями вектора скорости. Справа - максимальные касательные напряжения с нанесенными поверх ориентациями оси аз. Данные приведены для времени 0.28 млн лет.
Fig. 11. Summary table for Model II. The base set of parameters is specified in the text. Left column - absolute values of the flow velocity vector; orientations of the velocity vector are marked above. Right column - maximum tangential stress values; orientations of axis аз are marked above. The data are given for the time value of 0.28 Ma.
Кроме несоответствия ориентаций осей главных напряжений для модели II также получен крайне высокий (200-300 МПа) уровень напряжений, что на порядок выше реальных значений.
Резюмируя, можно утверждать, что результаты моделирования выявили несоответствие «плейт-тектонической» модели II и имеющихся для эпи-платформенных орогенов тектонофизических данных. Геодинамическое происхождение этого феномена стоит искать «на месте» - в подстилающей литосферу мантии или в самой литосфере ороген-ных областей. Моделирование показало, что именно вынужденная конвекция в литосфере, формирующаяся под влиянием мелкомасштабной астено-сферной конвекции, приводит к формированию НДС, наблюдаемого у эпиплатформенных орогенов Центральной Азии. Именно эта модель, с незначительными оговорками, показывает соответствие результатам тектонофизических реконструкций.
Следует отметить, что в нашем решении получены слишком большие времена формирования орогенов, превышающие реально наблюдаемые на порядок и даже более. Это связано с используемой для коры реологической моделью в виде вязкого тела. В реальности верхняя часть коры (5-10 км) упруга в зонах консолидированных блоков и вяз-копластична в зонах разломов, разделяющих эти блоки. Средняя и нижняя кора в большей степени отвечает упругопластическому телу [М1ко1авУБку, 1983]. Таким образом, в реальности приращения деформации пород коры определяются интенсивностью нагружения (идет мгновенная реакция на приращение нагружения), в то время как скорость деформаций в теле с вязкой реологией зависит от
скорости нагружения. Использование упругопла-стической модели коры и упруговязкой модели литосферы и астеносферы должно дать более реалистическое значение длительности процесса.
6. Заключение
В рамках исследования были созданы две математические модели эпиплатформенного орогенеза, отвечающие рассматриваемым чаще всего в качестве обусловливающих орогенез геодинамическим процессам. Анализ результатов моделирования проводился, в первую очередь, с учетом формирующейся в земной коре структуры напряженного состояния, которая сопоставлялась с данными тектонофизических реконструкций, обобщенных в работе [ЯвЪвЬБку, 2015].
Результаты моделирования показали, что под влиянием астеносферной конвекции происходит изостатическая инверсия образующихся форм рельефа в течение первых 50 млн лет, после чего система постепенно переходит в состояние частичной изостазии (100-120 млн лет), причем в случае прекращения действия астеносферной конвекции система достаточно быстро (8-10 млн лет) достигает изостазии. До инверсии в области поднятий наблюдается девиаторное растяжение, которое сменяется на сжатие на постинверсионной стадии, что соответствует природным данным. Однако в верхней части коры (10 км) формируется зона локального растяжения. Уровень напряжений (в среднем ~45 МПа) зависит от поверхностных процессов, учет которых был произведен и привел
к снижению амплитуды рельефа в 1.5 раза, но почти в 5 раз к повышению средней амплитуды скорости течения в коре и, соответственно, напряжений.
Вторая («плейт-тектоническая») модель не показала наличия инверсионных процессов, соответственно не привела к формированию корней/антикорней для поднятий и впадин, что явно противоречит природным данным, как и чрезмерно высокий уровень девиаторных напряжений. Кроме того, на зрелых этапах эволюции модели обстанов-
ка горизонтального сжатия формируется как во впадинах, так и в поднятиях, что не вполне соответствует природным данным о НДС эпиплатфор-менных орогенов Центральной Азии.
7. Благодарности
Работа выполнена при поддержке РФФИ (проекты № 18-35-00482 и № 16-05-01115) и в рамках Государственного задания ИФЗ РАН.
8. Литература / References
Ahnert F, 1970. Functional relationships between denudation, relief, and uplift in large, mid-latitude drainage basins. American Journal of Science 268 (3), 243-263. https://doi.Org/10.2475/ajs.268.3.243.
Belousov V.V., 1989. Fundamentals of Geotectonics. Nedra, Moscow, 382 p. (in Russian) [Белоусов В.В. Основы геотектоники. М.: Недра, 1989. 382 с.].
Blatt H., Middleton G., Murray R., 1972. Origin of Sedimentary Rock. Prentice-Hall, Inc., New Jersey, 634 p.
Bobrov A.M., Baranov A.A., 2011. Horizontal stresses in the mantle and in the moving continent for the model of two-dimensional convection with varying viscosity. lzvestiya, Physics of the Solid Earth 47 (9), 801-815. https:// doi.org/10.1134/S1069351311090023.
Bobrov A.M., Trubitsyn V.P., 2003. Evolution of viscous stresses in the mantle and in moving continents in the process of the formation and breakup of a supercontinent. lzvestiya, Physics of the Solid Earth 39 (12), 963-973.
Burtman V.S., 2012. Tian Shan and High Asia. Geodynamics in Cenozoic. GEOS, Moscow, 188 p. (in Russian) [Бурт-ман В.С. Тянь-Шань и Высокая Азия. Геодинамика в кайнозое. М.: ГЕОС, 2012. 188 с.].
Cloetingh S., Burov E., Beekman F., Andeweg B., Andriessen P.A.M., García-Castellanos D., de Vicente G., Vegas R., 2002. Lithospheric folding in Iberia. Tectonics 21 (5), 1041. https://doi.org/10.1029/2001TC901031.
Culling W.E.H., 1960. Analytical theory of erosion. The Journal of Geology 68 (3), 336-344. https://doi.org/10.1086/ 626663.
Dobretsov N.L., Berzin N.A., Buslov M.M., Ermikov V.D., 1995. General aspects of the evolution of the Altai region and the interrelationships between its basement pattern and the neotectonic structural development. Geologiya i Geofizika (Russian Geology and Geophysics) 36 (10), 5-19 (in Russian) [Добрецов Н.Л., Берзин Н.А., Буслов М.М., Ермиков В.Д. Общие проблемы эволюции Алтайского региона и взаимоотношения между строением фундамента и развитием неотектонической структуры // Геология и геофизика. 1995. Т. 36. № 10. С. 5-19].
Dobretsov N.L., Koulakov l.Y., Polyansky O.P., 2013. Geodynamics and stress-strain patterns in different tectonic settings. Russian Geology and Geophysics 54 (4), 357-380. https://doi.org/10.1016Zj.rgg.2013.03.001.
England P., Houseman G., 1986. Finite strain calculations of continental deformation: 2. Comparison with the India-Asia collision zone. Journal of Geophysical Research: Solid Earth 91 (B3), 3664-3676. https://doi.org/10.1029/ JB091iB03p03664.
García-Castellanos D., Cloetingh S.I.E.R.D., Busby C., Azor A., 2012. Modeling the interaction between lithospheric and surface processes in foreland basins. In: C. Busby, A. Azor (Eds.), Tectonics of sedimentary basins: recent advances. Wiley-Blackwell, Oxford, p. 152-181. https://doi.org/10.1002/9781444347166.ch8.
Horowitz D.H., 1976. Mathematical modeling of sediment accumulation in prograding deltaic systems. In: D.F. Merriam (Ed.), Quantitative techniques for the analysis of sediments. Pergamon, Oxford, p. 105-119. https://doi.org/ 10.1016/B978-0-08-020613-4.50015-X.
Khain V.E., Lomize M.G., 1995. Geotectonics with Fundamentals of Geodynamics. Publishing House of Moscow State University, Moscow, 480 p. (in Russian) [Хаин В.Е., Ломизе М.Г. Геотектоника с основами геодинамики. М.: Изд-во Московского госуниверситета, 1995. 480 с.].
Krestnikov V.N., Belousov T.P., Ermilin V.l., Chigarev N.V., Shtange D.V., 1979. Quaternary Tectonics of the Pamirs and Tien Shan. Nauka, Moscow, 116 p. (in Russian) [Крестников В.Н., Белоусов Т.П., Ермилин В.И., Чигарев Н.В., Штанге Д.В. Четвертичная тектоника Памира и Тянь-Шаня. М.: Наука, 1979. 116 с.].
Landau L.D., Lifshits E.M., 1986. Hydromechanics. Theoretical Physics. Volume 6. 3rd Edition. Nauka, Moscow, 736 p. (in Russian) [Ландау Л.Д., Лифшиц Е.М. Гидромеханика. Теоретическая физика. Том 6. 3-е изд., перераб. М.: Наука, 1986. 736 с.].
Liu J., Liu Q.-Y., Guo B., Yuen D.A., Song H.-Z., 2007. Small-scale convection in the upper mantle beneath the Chinese Tian Shan Mountains. Physics of the Earth and Planetary Interiors 163 (1-4), 179-190. https://doi.org/10.1016/j.pepi. 2007.04.019.
Love A.E.H., 1911. Some Problems of Geodynamics. Cambridge University Press, London, 180 p.
Makarov V.I. (Ed.), 2005. Recent Geodynamics of Areas of Intracontinental Collision Mountain Building (Central Asia). Nauchny Mir, Moscow, 400 p. (in Russian) [Современная геодинамика областей внутриконтинентального коллизионного горообразования (Центральная Азия) / Ред. В.И. Макаров. М.: Научный мир, 2005. 400 с.].
Makarov V.I., Alekseev D. V., Batalev V. Y., Bataleva E.A., Belyaev I. V., Bragin V.D., Dergunov N.T., Efimova N.N., Leonov M.G., Munirova L.M., Pavlenkin A.D., Roecker S., Roslov Yu.V., Rybin A.K., Shchelochkov G.G., 2010. Underthrusting of Tarim beneath the Tien Shan and deep structure of their junction zone: Main results of seismic experiment along MANAS Profile Kashgar-Song-Kol. Geotectonics 44 (2), 102-126.
Makarov V.I., Rybin A.K., Matyukov V.E., Pushkarev P.Yu., Shcherbina F.A., 2011. Features of the deep structure of depression areas of the Central Tien Shan. Inzhenernye Izyskaniya (Engineering Surveys) (1), 42-51 (in Russian) [Макаров В.И., Рыбин А.К., Матюков В.Е., Пушкарев П.Ю., Щербина Ф.А. Особенности глубинной структуры депрессионных областей Центрального Тянь-Шаня // Инженерные изыскания. 2011. № 1. С. 42-51].
Mikhailov V.O., 1983. Mathematical model of the evolution of structures formed as a result of vertical movements. Izvestiya AN SSSR, Seriya Fizika Zemli (6), 3-18 (in Russian) [Михайлов В.О. Математическая модель эволюции структур, образующихся в результате вертикальных движений // Известия АН СССР, серия Физика Земли. 1983. № 6. С. 3-18].
Mikhailov V.O., 1999. Modeling the Extension and Compression of the Lithosphere by Intraplate Forces. Izvestiya, Physics of the Solid Earth 35 (3), 228-238.
Mikhailov V.O., Timoshkina E.P., Polino R., 1999. Foredeep basins: the main features and model of formation. Tectono-physics 307 (3-4), 345-359. https://doi.org/10.1016/S0040-1951(99)00052-9.
Myagkov D.S., Rebetsky Y.L., 2016. The evolution of structure of crustal flow and relief of epi-platformal orogens effected by small-scale convection in the asthenosphere. Bulletin of Kamchatka Regional Association Education-Science Centre. Earth Sciences (1), 89-100 (in Russian) [Мягков Д.С., Ребецкий Ю.Л. Эволюция структуры течения и рельефа коры эпиплатформенных орогенов под воздействием мелкомасштабной астеносферной конвекции // Вестник КРАУНЦ. 2016. № 1. С. 89-100].
Nikolaevsky V.N., 1983. Mechanics of geomaterials and earthquakes. Results in Science and Technique. Mechanics of De-formable Solid Body Series 15, 817-821 (in Russian) [Николаевский В.Н. Механика геоматериалов и землетрясения // Итоги науки и техники, серия Механика деформируемого твердого тела. 1983. Т. 15. С. 817-821].
Ollier C.D., 1981. Tectonics and Landforms (Geomorphology Texts). Addison-Wesley Longman, London, 324 p. [Русский перевод: Оллиер К. Тектоника и рельеф. М.: Недра, 1984. 460 с.].
Rabotnov Yu.N., 1988. Mechanics of a Deformable Solid. Nauka, Moscow, 712 p. (in Russian) [Работнов Ю.Н. Механика деформируемого твердого тела. M.: Наука, 1988. 712 с.]
Rebetsky Y.L., 2012. About one new form of instability of the continental crust. In: Sedimentary basins and geological prerequisites for forecasting new objects promising for oil and gas. Proceedings of the XLIV Tectonic Meeting. Vol. II. GEOS, Moscow, p. 355-359 (in Russian) [Ребецкий Ю.Л. Об одной новой форме неустойчивости континентальной коры // Осадочные бассейны и геологические предпосылки прогноза новых объектов, перспективных на нефть и газ: Материалы XLIV тектонического совещания. Т. II. М.: ГЕОС, 2012. С. 355-359].
Rebetsky Y.L., 2015. On the specific state of crustal stresses in intracontinental orogens. Geodynamics & Tectonophysics 6 (4), 437-466 (in Russian) [Ребецкий Ю.Л. Об особенности напряженного состояния коры внутриконтинен-тальных орогенов // Геодинамика и тектонофизика. 2015. Т. 6. № 4. С. 437-466]. https://doi.org/10.5800/ GT-2015-6-4-0189.
Rebetsky Y.L., Alekseev R.S., 2014. The field of recent tectonic stresses in Central and South-Eastern Asia. Geodynamics & Tectonophysics 5 (1), 257-290 (in Russian) [Ребецкий Ю.Л., Алексеев Р.С. Поле современных тектонических напряжений Средней и Юго-Восточной Азии // Геодинамика и тектонофизика. 2014. Т. 5. № 1. С. 257-290]. https://doi.org/10.5800/GT-2014-5-1-0127.
Rebetsky Y.L., Kuchai O.A., Marinin A.V., 2013. Stress state and deformation of the Earth's crust in the Altai-Sayan mountain region. Russian Geology and Geophysics 54 (2), 206-222. https://doi.org/10.1016Zj.rgg.2013.01.011.
Rebetsky Y.L., Kuchai O.A., Sycheva N.A., Tatevossian R.E., 2012. Development of inversion methods on fault slip data: Stress state in orogenes of the Central Asia. Tectonophysics 581, 114-131. https://doi.org/10.1016/j.tecto.2012. 09.027.
Rebetsky Y.L., Pogorelov V.V., 2013. Tectonophysical model of the loading mechanism and the evolution of the stressstrain state of the lithosphere of continental mountain-folded areas. In: Geological history, possible mechanisms, and the problem of formation of depressions with suboceanic and anomalously thin crust in continental lithosphere provinces. Proceedings of the XLV Tectonic Meeting. Vol. II. GEOS, Moscow, p. 181-185 (in Russian) [Ребецкий Ю.Л., Погорелов В.В. Тектонофизическая модель механизма нагружения и эволюции напряженно-деформированного состояния литосферы континентальных горно-складчатых областей // Геологическая история, возможные механизмы и проблема формирования впадин с субокеанической и аномально тонкой корой в провинциях с континентальной литосферой: Материалы XLV тектонического совещания. Т. II. М.: ГЕОС, 2013. С. 181-185].
Rebetsky Y.L., Sycheva N.A., Sychev V.N., Kuzikov S.I., Marinin A.V., 2016. The stress state of the northern Tien Shan crust based on the KNET seismic network data. Russian Geology and Geophysics 57 (3), 387-408. https://doi.org/ 10.1016/j.rgg.2016.03.003.
Schubert B., Turcotte D.L., Olson P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press, Cambridge, 940 p.
Thieulot C., Steer P., Huismans R.S., 2014. Three-dimensional numerical simulations of crustal systems undergoing orogeny and subjected to surface processes. Geochemistry, Geophysics, Geosystems 15 (12), 4936-4957. https:// doi.org/10.1002/2014GC005490.
Timoshkina E.P., Leonov Y.G., Mikhailov V.O., 2010. Formation of the orogen-foredeep system: A geodynamic model and comparison with the data of the northern Forecaucasus. Geotectonics 44 (5), 371-387. https://doi.org/10.1134/ S0016852110050018.
Trubitsyn V.P., Simakin A.G., Baranov A.A., 2006. The effect of spatial variations in viscosity on the structure of mantle flows. lzvestiya, Physics of the Solid Earth 42 (1), 1-12. https://doi.org/10.1134/S1069351306010010.
Turcotte D.L., Schubert G., 1982. Geodynamics: Application of Continuum Physics to Geological Problems. John Wiley & Sons, New York, 464 p.
СВЕДЕНИЯ ОБ АВТОРАХ | INFORMATION ABOUT AUTHORS
Дмитрий Сергеевич Мягков н.с.
Институт физики Земли им. О.Ю. Шмидта РАН
123242, ГСП-5, Москва Д-242, ул. Большая Грузинская, 10, Россия
И e-mail: [email protected]
Dmitriy S. Myagkov
Researcher
O.Yu. Schmidt Institute of Physics of the Earth of RAS
10 Bol'shaya Gruzinskaya street, Moscow D-242 123242, GSP-5, Russia
Юрий Леонидович Ребецкий
докт. физ.-мат. наук, зав. лабораторией
Институт физики Земли им. О.Ю. Шмидта РАН
123242, ГСП-5, Москва Д-242, ул. Большая Грузинская, 10, Россия
e-mail: [email protected]
© https://orcid.org/0000-0003-3492-2452
Yuri L. Rebetsky
Doctor of Physics and Mathematics, Head of Laboratory
O.Yu. Schmidt Institute of Physics of the Earth of RAS
10 Bol'shaya Gruzinskaya street, Moscow D-242 123242, GSP-5, Russia
