Тепломассоперенос в континентальной коре в условиях гравитационной неустойчивости: геологический анализ и многопроцессорное моделирование (стр. 1 )

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

О Т Ч Е Т

по междисциплинарному интеграционному проекту

фундаментальных исследований

Сибирского отделения РАН №2

за 2009 год

Тепломассоперенос в континентальной коре

в условиях гравитационной неустойчивости:

геологический анализ и многопроцессорное моделирование

Организации – исполнители:

Институт геологии и минералогии им. СО РАН

Институт гидродинамики им. СО РАН

Институт вычислительных технологий СО РАН

Научный координатор проекта:

Академик РАН , советник РАН

Новосибирск

2009 г.

Цели проекта.

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

Результаты исследований. I. Геологические примеры гравитационно-неустойчивых систем в земной коре.

Среди многообразия проявлений гранитоидного магматизма в континетальной коре мы сосредоточились в рамках проекта на диапиризме и пока не рассматриваем альтернативный механизм – гранитные интрузии. Диапиризм гранитной магмы широко развит в архее и нижнем протерозое. Кислый расплав в нижней коре может сформироваться в следующих случаях: 1) в процессе термической релаксации утолщенной коры в орогенических областях; 2) за счет регионального увеличения теплового потока из мантии над поднятиями астеносферы; 3) при глубинных интрузиях базитов под основание коры (андерплейтинга); 4) ввиду внутрикорового разогрева вследствие аномального распределения радиоактивных элементов в существенно гранитной коре. Дальнейшее движение расплава определяется перераспределением вещества по плотности в гравитационном поле, происходит путем интрузий или диапиризма и контролируется рядом факторов, главными из которых являются температура и реология вмещающей среды.

Геологическим примером диапировой структуры может являться Тейский гранито-гнейсовый купол, расположенный в заангарской части Енисейского кряжа [Ножкин и др., 1983; 1999]. Купольная структура площадью 1500 км2 является самым верхним элементом структурно и орографически приподнятого сиалического блока в пределах центральной зоны Енисейского кряжа. Ядерная часть купола образована гнейсами, гранито-гнейсами, порфиробластическими гнейсо-гранитами, а также интрузивными магматическими гранитами, средняя плотность которых составляет 2600 кг/м3 [13]. Вмещающие породы – метаморфические толщи тейской и сухопитской серии протерозойского возраста, состоящие из высокоглиноземистых гранат-ставролит-биотитовых кристаллических сланцев (плотностью 2960 кг/м3) и ортоамфиболитов (2920 кг/м3). Контакты гранитоидного массива с вмещающими породами преимущественно согласные, но на отдельных участках в плане и на глубине (по геофизическим данным) они становятся секущими. Тейский массив выделяется как крупное поднятие первого порядка, осложненное куполами второго порядка (Итуйский, Индольский, Тырыдинский купола площадью 6-10 км2). Последние характеризуются кольцевыми и полукольцевыми разломами, залегающими с крутым падением в сторону периферии купола. Геохронологические исследования U-Pb методом по цирконам, отобранным во внешней периферийной зоне Итуйского купола, показали возраст одного из этапов формирования Тейского гранито-гнейсового купола – 866±16 млн. лет. Этот возрастной рубеж подтвержден позднее SHRIMP-анализом цирконов [Верниковский, Верниковская 2006], отобранных из гранитоидов северной части Тейского массива (864±9–868±10 млн. лет). Достаточно близкие оценки возраста говорят о том, что структурные элементы разного порядка Тейского гранитогнейсового купола формировались практически одновременно. Наши модельные расчеты (см. раздел II) указывают на более быструю скорость всплывания диапира относительно скорости кондуктивного нагрева, а затем - охлаждения гранитоидного массива. Так как температура закрытия U-Th-Pb системы составляет около 940ºС, а всплывание возможно только при высоких температурах, возраст, определенный по цирконам соответствует тому моменту времени, когда диапир полностью поднялся.

Другими примерами гранито-гнейсовых диапиров являются купол Фангшан (Северо-Китайский кратон) (He et al., 2009) и купол г. Тор-О’дин (Канадские Кордильеры) (Norlander et al., 2002).

Параметры формирования купольной структуры Фангшан следующие. Возраст кристаллизации – 130 млн. л. (данные по датированию цирконов методом SHRIMP (Davis, 2001)), возраст охлаждения 128.7 млн. л. (данные по Ar/Ar-датированию (Ma, 1989)), глубина внедрения – 6-7 км (данные по минеральному геобарометру 2 кб). Ширина метаморфической зональности варьирует в диапазоне м при размере гранодиоритового ядра купола 8-10 км. Наблюдается обрамляющая синклинальная структура по периферии ядра купола, что является признаком совместной деформации гранитоидов и вмещающих пород при подъеме диапира. Эрозия вмещающих метаосадков над головной частью диапира составляет более 4 км, рассчитанная по сравнению с мощностью неэродированных осадочных отложений на удалении от купольной структуры. Это означает, что минимальная высота подъема диапира была не менее толщины эродированных пород. Гранодиоритовое ядро окружено высокотемпературным концентрическим ореолом сдвигово-деформированных пород.

Мигматитовый купол г. Тор-Один представляет собой один из совокупности куполов плутонического пояса Роки Маунтин (Западные Кордильеры, Канада) (Norlander et al., 2002). Он представляет собой овал размером 5-7 на 20 км. Ядро составляют мигматито-гнейсы, лейкограниты; обрамление купола сложено грубообломочными метабазитами зеленосланцевой фации. По реконструкции термохронологической истории формирования купола выделяются следующие этапы формирования купола: I этап – погружение пород при коллизии орогена и анатексис, II этап – изотермическая декомпрессия от 10 до 4-5 кбар (17-21 км) в процессе всплывания дипира, III этап – растяжение и эксгумация ядер диапиров и метаморфических пород обрамления. Длительность всплывания диапира оценивается приблизительно в ~4 млн. лет на основе датирования зерен цирконов методом SHRIMP в диапазонемлн. лет (Vanderhaege, 1999). Полученные из геологического анализа параметры процесса диапирового подъема гранитной магмы применялись для построения теплофизической и компьютерной моделей, описанных в разделе II и III.

II,а. Моделирование диапиризма гранито-гнейсов в земной коре с учетом плавления.

В разделах II, а-в приводятся результаты двух подходов в моделировании диапиризма: 1) на основе механики деформируемого твердого тела (МДТТ); 2) с применением гидродинамики вязкой ньютоновой жидкости. По каждому направлению разработаны самосогласованные модели, выполнена их численная реализация в виде программных приложений, решены тестовые задачи. Произведена закупка коммерческого пакета программ MSC. Marc2008, реализующего метод конечных элементов.

В результате выполнения этапа проекта разработан новый подход, описывающий процессы частичного плавления и развития гравитационной неустойчивости в рамках механики деформируемого твердого тела. Постановка задачи состоит в следующем. Рассматривается двумерная прямоугольная область земной коры размером 38х60 км (глубина/ширина) (рис.1). Решаются уравнения механического равновесия в «слабой» форме (уравнение принципа возможных скоростей перемещений или уравнение баланса виртуальных мощностей) и уравнение теплопроводности с переменными коэффициентами и учетом фазового перехода при плавлении. Классическая постановка задачи Стефана заменяется учетом повышенной теплоемкости при фазовом переходе. Уравнения, лежащие в основе математического моделирования квазистатического деформирования, решались численно в рамках приближения задачи о плоской деформации.

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

Для дискретизации уравнений МДТТ использовался метод конечных элементов. Для численного моделирования использовался пакет программ MSC. Marc 2005, в котором предусмотрен учет всех типов нелинейности уравнений МДТТ. Уравнения квазистатического движения и теплопроводности решались в верхней области моделирования (верхних 30 км), а в нижней области (8 км) решалось только уравнение теплопроводности. Результаты расчетов приводятся только для верхней части модели. Ввиду больших деформаций сетка в верхней области перестраивалась по условию достижения критической деформации элементов: при изменении угла между ребрами элемента более 20° при средней длине ребра 300 м.

Граничные условия выбирались следующими: верхняя граница – свободная поверхность с постоянной температурой Т=0°С; боковые границы изолированы для передачи тепла и неподвижны в отношении механического движения в горизонтальном направлении; нижняя граница фиксирована в вертикальном направлении и допускает свободное скольжение вдоль нее. Начальные условия принимались такими: отсутствие начальных перемещений, литостатическое распределение напряжений, температура соответствовала геотермическому градиенту 18°С/км при экспоненциальном распределении по глубине коры радиоактивных элементов. В середине нижней коры задается блок размером 8х20 км, состоящий из базитового расплава при начальной температуре 1200°С. На остальной части нижней границы принималась постоянная температура 545°С или постоянный мантийный тепловой поток 0.03 Вт/м2 (рис. 1). Предполагается, что 8-километровый слой базитовых интрузий внедряется в основании нижней коры в течение 2 млн. лет, что обусловило прогрев и частичное плавление вышележащей гранитной коры. Используется диаграмма плавления водосодержащего гранита, для которого температура солидуса выражается аналитически в виде:

где Tsol ‑ температура в K, P ‑ давление в МПа. Принималось, что при плавлении в гранитной коре в интервале температур 650-700°С происходит уменьшение плотности частично расплавленных пород до 2600 кг/м3, при плотности исходных пород 2800 кг/м3. Коэффициенты теплопроводности и теплоемкости расплавленной и вмещающей породы задавались одинаковыми и составляли 2 Вт/(м∙К) и 1250 Дж/(кг∙К), соответственно, скрытая теплота плавления гранита принималась равной 300 кДж/кг. Земная кора принимается однородной и изотропной с упругими характеристиками: модуль Юнга , коэффициент Пуассона . Минимальная температура перехода от хрупкого к пластическому поведению для разных типов пород зависит от скорости деформации, всестороннего давления, водонасыщенности и составляет 250°С для влажного кварцита, 280-300°С для сухого кварцита, 500-550°С для водонасыщенного дунита. Экспериментальные данные говорят о том, что, начиная с некоторых глубин и температур, материал коры подчиняется законам пластической деформации. Поэтому была выбрана упруго-пластическая модель материала среды, в которой пластичность описывалась уравнением Хубера-Мизеса , где – предел текучести (МПа); – второй инвариант тензора–девиатора напряжений , где ‑ компоненты тензора напряжений Коши.

Реология частично-расплавленного материала, как и вмещающей породы, предполагалась температурно-зависимой, подчиняющейся закону идеальной пластичности с пределом текучести, увеличивающимся от 1 МПа в расплаве до 10 МПа в окружающей среде. Температурная зависимость предела текучести по аналогии с вязкостью принималась в форме закона Аррениуса σy = 106 (Па)+A exp(-E/RT), где A и E – эмпирические константы. Модель строилась для описания процесса плавления и вызванного им всплывания легкого вещества в результате андерплейтинга базитовой магмы в основании континентальной коры. Цель моделирования состояла в определении структуры течения всплывающей гранитной магмы и предсказании возможной формы гранито-гнейсовых диапировых тел. Были проведены 2D численные эксперименты с разными пределами текучести материала коры. Важным параметром является степень плавления и, следовательно, разность плотностей расплавленной и твердой породы (50-200 кг/м3). Для упрощения анализа принято постоянное значение 200 кг/м3; таким образом, считалось, что степень затвердевания охлаждающегося расплава не учитывается. Решения получены в виде двумерных картин распределения полей напряжений, деформаций, температуры в разные моменты времени. Первая серия тестов проводилась при постоянной величине предела текучести σy =1 МПа, не зависящего от температуры. Результаты моделирования представлены на рис. 2 в виде картин температурного поля промежуточных состояний всплывающего расплава. В начальный этап (~ 2 млн. лет) происходит распространение фронта частичного плавления над тепловым источником. Кондуктивный прогрев и формирование области частичного плавления, которое происходит без заметного движения материала коры, при заданных параметрах длится около 2 млн. лет. За это время формируется область с поверхностью в форме купола высотой 6.7 км и шириной, определяемой заданным размером базитовой интрузии (около 20 км), поставляющей тепло. Далее начинается всплывание расплава. Конвективное движение преобразует форму поверхности в грибовидную (рис. 2, а-в), затем снизу формируется канал, который сужается в процессе всплывания (рис. 2в-г). Ширина канала составляет около 3.5 км, а ширина головной части диапира – около 13 км. Верхняя свободная поверхность испытывает перемещения порядка сотни метров с максимальным воздыманием над «головой» диапира. На поздних временах (более 2 млн. лет) подъем диапира прекращается и следует стадия охлаждения (затвердевания) на месте. Этот момент фиксируется, когда происходит «пережатие» канала и отрыв головной части диапира (рис. 2д). По результатам расчетов сделаны оценки скорости проплавления нижней коры при кондуктивном (начальном) этапе теплопереноса и скорости диапирового всплывания при конвективном (последующем) этапе. По модельным оценкам фронт плавления перемещается со скоростью 4.7 мм/год, а совместное движение фронта плавления с учетом всплывания происходит со скоростью 10 мм/год. Из-за большой скорости всплывания плавление новых порций вещества на верхних уровнях коры не происходит, и всплывает лишь объем расплава, сформированного на начальной стадии.