Моделирование фазовых равновесий при кристаллизации базальтовых магм

Заключительная глава книги посвящена проблемам ЭВМ-моделирования дифференциации основных магм в плоскопараллельной закрытой магматической камере. Хотя исторически эта часть исследования была выполнена раньше, чем получены данные геохимической термометрии и проанализированы результаты фракционирования в природных системах, мы рассматриваем их в последнюю очередь, учитывая, что здесь представлены наиболее сложные модели, сочетающие расчеты реалистичных фазовых диаграмм с динамикой остывания и дифференциации магмы. Работы по моделированию процессов внутрикамерной дифференциации были начаты в первой половине 70-х годов группой сотрудников ГЕОХИ РАН и кафедры геохимии МГУ по инициативе А.А.Ярошевского. Эти исследования шли сразу по нескольким направлениям, включая работы по моделированию термики остывания пластовых интрузивов с учетом перемещения продуктов кристаллизации на границах камеры и в основном объеме магмы (М.Я.Френкель), построение реалистичной модели фазовой диаграммы для кристаллизующихся базальтов (А.А.Арискин) и тестирование разработанных моделей на природных объектах (Г.С.Бармина, Б.С.Киреев, Е.В.Коптев-Дворников). В начале 80-х годов результаты этих исследований удалось совместить в рамках единого алгоритма - программы ИНТРУЗИВ. Теоретические основы построения этой сложной физико-химической модели были сформулированы в работах (Френкель, 1982; Френкель и др., 1985).

Наша задача состояла в реализации этих идей, совершенствовании процедуры расчета фазовых равновесий, отладке методических приемов работы с программой и анализе результатов ЭВМ-моделирования. Во второй половине 80-х годов при помощи программы ИНТРУЗИВ были проведены расчеты по динамике дифференциации трапповых магм Сибирской платформы - эти исследования изложены в серии публикаций (Френкель и др., 1985; Френкель, Арискин, 1990; Frenkel et al., 1989) и обобщены в коллективной монографии (Френкель и др., 1988). По этой причине мы считаем возможным только коротко описать важнейшие результаты, уделив главное внимание анализу деталей динамического ('интрузивного") блока модели КОМАГМАТ и новым данным, касающихся использования этой программы применительно к проблемам генезиса комплекса изверженных пород Садбери (Ariskin et al., 2000).

Среди важнейших достижений в понимании процессов, протекающих в магматических камерах, особое место занимают исследования конвекционных режимов, ответственных за дифференциацию магмы и формирование структурно-химических особенностей пород. Вслед за пионерской работой (Bartlett, 1969) петрологи акцентировали внимание на оценках числа Рэлея, как главного показателя эффективности термической конвекции в гомогенной магме. В последнее десятилетие было показано, что важным фактором конвекции является скорость зарождения и распределение кристаллов, а также градиент температуры в гетерогенном пограничном слое, примыкающем к верхнему фронту полного затвердевания (Френкель и др., 1988; Morse, 1988; Marsh, 1988; 1989; Worster et al., 1990; Mangan, Marsh, 1992; Simakin et al., 1994; Френкель, 1995; Jaupart, Tait, 1995). Между тем результаты этих работ оказались противоречивы.

В частности, Marsh (1988) пришел к выводу, что для интрузивных тел, охлаждающихся преимущественно через верхний контакт, кристаллизация вблизи кровли должна препятствовать проявлениям глобальной термической конвекции в камере. Этот вывод основан на постулате, что основной объем магмы находится в изотермических условиях, отвечающих температуре ликвидуса на границе переходного гетерогенного слоя. Утверждалось, что высокая температура магмы "не позволяет" кристаллам перемещаться (оседать) в объеме расплава и достигать дна камеры, что ставит под сомнение возможность внутрикамерной дифференциации по седиментационному механизму (Marsh, 1989).

Противоположная точка зрения была высказана в работе (Worster et al.,1990), главный вывод которой состоит в том, что "общая конвекция, значительное внутреннее охлаждение и дифференциация могут происходить в магматической камере, которая остывает в результате кондуктивного отвода тепла в перекрывающие породы, даже при условии отсутствия начального перегрева (магмы)". М.Я.Френкель и А.А.Ярошевский с коллегами еще в 70-е годы исследовали эти проблемы путем построения численной модели внутрикамерной дифференциации основных магм, сопоставляя свойства модельных систем с наблюдаемыми характеристиками природных объектов.

В ГЕОХИ РАН была разработана серия ЭВМ-программ, моделирующих процессы кристаллизационной дифференциации магмы в плоскопараллельной камере при наличии градиентов температуры и фазовых переходов, вызванных оседанием кристаллов в объеме магматического расплава (Френкель, Ярошевский, 1976, 1978; Коптев-Дворников и др., 1979). В этих алгоритмах была задействована сильно упрощенная модель фазовой диаграммы базальтов (типа Di-An-Fo - Рис. 6.1), но задача тепломассопереноса решалась максимально корректно, путем разбиения слоя магмы (мощностью 200 м) на ряд интервалов (не менее 50-100), для каждого из которых на заданных временных отрезках контролировался баланс тепла и вещества. Это позволило проследить эволюцию температуры, фазового и химического состава магмы вдоль вертикального разреза камеры по мере затвердевания интрузивного тела в зависимости от исходного состава, степени кристалличности и краевых условий.

Рис. 6.1. Модельная диаграмма кристаллизации, использованная при разработке ранних моделей внутрикамерной дифференциации магматических расплавов (Френкель и др., 1978) Точками обозначены составы исходных расплавов при исследованиях свойств седиментационной и конвекционно-кумуляционной модели.

Отличительная особенность этих моделей - возможность варьирования скоростей оседания кристаллов Ol, Pl и Px. Используя имеющиеся в то время средства вычислительной техники (БЭСМ-4 и БЭСМ-6), были проведены десятки расчетов, моделирующих затвердевание интрузивных базитов разного состава. Напомним, что возможности этих машин и реального доступа к "времени" требовали не менее 2-3 недель для завершения одной серии вычислений, включающих полную историю остывания и дифференциации магмы. Значение этих результатов трудно переоценить.

Во-первых, было показано, что растворение оседающих кристаллов в нижележащем относительно высокотемпературном расплаве действительно имеет место, но вследствие общего остывания системы этот процесс не в состоянии предотвратить продвижение "облака кристаллической фазы" в направлении дна магматической камеры и появление кумулятов. Во-вторых, установлено, что седиментация кристаллов является гораздо более мощным механизмом теплообмена по сравнению с кондуктивным переносом. Это приводит к тому, что температура во внутренних частях камеры достаточно быстро выравнивается, достигая значений, отвечающих эвтектике Ol-Pl-Px. Конечное распределение минералов по разрезу (порядок смены кумулятивных парагенезисов) определяется последовательностью кристаллизации в переходном слое и скоростями оседания кристаллов (Френкель, Ярошевский, 1978). При этом удалось впервые воспроизвести прогрессивное обогащение пород оливином вглубь интрузивного тела от нижнего контакта и обеднение от верхнего - так называемый S-образный профиль модальных вариаций, часто наблюдаемый в интрузивных телах, см. напр. Рис. 1 в работе (Mangan, Marsh, 1992). Детальный анализ петрологических следствий седиментационной модели дифференциации представлен в работе (Коптев-Дворников и др., 1979), где структура модельных тел была сопоставлена с данными по строению силлов Сибирской платформы.

Результаты этой работы показали, что многие особенности строения реальных интрузивов воспроизводятся в модельных объектах: это общее преобладание в породах кумулятивного материала, порядок смены кристаллизующихся ассоциаций, строение верхней и нижней приконтактовой зоны, поведение микроэлементов. Все эти признаки указывали на определяющую роль оседания кристаллов как фактора внутрикамерной дифференциации. Проблема состояла в том, что степень дифференцированности модельных интрузивов (контрастность минерального состава пород) оказалась гораздо выше по сравнению с данными о строении дифференцированных силлов мощностью 100-200 м. Этот результат проявлялся независимо от принятой скорости оседания кристаллов: даже в случае низких значений порядка первых м/год границы между модельными кумулятами носили резкий, скачкообразный характер. Анализ этого вопроса привел к выводу, что в реальных интрузивах проявляется дополнительный процесс, дестабилизирующий термические и фазовые границы внутри камеры - общая конвекция.

В разделе 1.5 мы уже отмечали, что главной физической причиной общей конвекции (глобального перемешивания магмы) является инверсия плотности, которая возникает как результат кристаллизации расплава в прикровельной части интрузива. Это приводит к гидродинамической неустойчивости системы и появлению разнонаправленных конвекционных токов магмы, относительно обогащенной и обедненной твердой фазой. Теоретический анализ этих явлений был представлен в конце 80-х - первой половине 90-х годов (Френкель и др., 1988; Френкель, 1995), тогда как на раннем этапе работ по ЭВМ-моделированию больше внимания уделялось исследованию свойств конвекционно-кумуляционной модели дифференциации, допускающей полное перемешивание основного объема магмы. Технически эта задача решалась просто - на каждом временном отрезке теплосодержание и валовый состав магмы выравнивался по всему объему; в результате происходило выравнивание температуры и фазового состава магматической смеси (Рис. 6.2). При этом постулировалось, что возможность перемешивания магмы не препятствует интегральному оседанию (или всплыванию) кристаллов данного минерального вида: у каждого индивида в дополнение к стоксовской появляется конвективная составляющая скорости.

Рис. 6.2. Динамика внутрикамерной дифференциации в соответствии с конвекционно-кумуляционной моделью Френкеля и др. (1988) А - эволюция фазового состава магмы и положение фронтов кумуляци и полного затвердевания в зависимости от времени; Б - вертикальный разрез магматической камеры в произвольный момент 1 - Ol+l, 2 - Ol+Pl+l, 3 - Ol+Pl+Px+l Рисунок показывает, что главный процесс химического разделения компонентов системы прекращается после заполнения камеры осадком кристаллических фаз (кумулусом).

Подобный подход представляет крайнее упрощение реальной ситуации в камере, но здесь можно заметить следующее. Если эффективность конвективного перемешивания в магматическом процессе была действительно высока, то можно ожидать, что при достаточно грубом сеточном разбиении по времени характеристики модельных тел могут соответствовать основным особенностям строения природных объектов. Именно это соответствие продемонстрировали первые результаты расчетов с использованием конвекционно-кумуляционной модели, основанной на примитивной фазовой диаграмме: фактор перемешивания проявился в постепенных переходах от одних кумулятов к другим (Коптев-Дворников и др., 1979) и плавных изгибах эволюционных линий для содержаний микроэлементов (Френкель, 1982).

Таким образом, к середине 80-х годов сложились предпосылки для построения более реалистичной модели дифференциации трапповых магм. С одной стороны был накоплен опыт моделирования динамики дифференциации и остывания пластовых интрузивов, выявлены главные механизмы перераспределения тепла и вещества в магматической камере. С другой - разработаны модели кристаллизации базальтовых расплавов, дающие хорошие результаты для толеитовых систем. Задача состояла в том, чтобы "связать" эти алгоритмы или представить некоторую аналитическую аппроксимацию динамики конвекционно-кумуляционной модели, совместимую с расчетами траекторий кристаллизации базальтовой магмы. Было решено пойти по второму пути.

Дело в том, что возможность построения строгой модели остывания и дифференциации магмы, основанной на сеточном разбиении и учитывающей потоки тепло- и вещества внутри камеры, была обусловлена линеаризованной формой диаграммы кристаллизации (Рис. 6.1), допускающей вывод простых алгебраических соотношений, связывающих температуру и фазовый состав системы с теплосодержанием - см. описание блока ДИАГРАММА в книге (Френкель и др., 1988). Это позволяло максимально эффективно использовать возможности доступной вычислительной техники при многократных обращениях (десятки и сотни тысяч раз на каждом уровне) к процедуре расчета фазового состава. Описанный в этой книге алгоритм моделирования траекторий кристаллизации построен на других принципах и позволяет решать задачу равновесия при известном соотношении кристаллов и расплава, но не заданном запасе тепла в системе.

Используя результаты ранних работ по ЭВМ-моделированию, М.Я.Френкель предложил ряд аналитических зависимостей, аппроксимирующих процессы остывания и дифференциации магмы по конвекционно-кумуляционному механизму. Перед нами была поставлена задача разработки ЭВМ-программы, объединяющей эту систему динамических уравнений с термодинамической моделью кристаллизации в единый алгоритм, который можно использовать при моделировании расслоенности пластовых интрузивов конкретного состава и мощности. Основные соотношения и структура этого алгоритма, интегрированного в модели КОМАГМАТ в виде подпрограммы ИНТРУЗИВ, рассмотрены ниже (Ariskin et al., 1993).

Возможность разработки новой, петрологически ориентированной версии конвекционного-кумуляционной модели связана с тем, что учет конвекции как фактора выравнивания температуры и состава по объему магмы позволяет отвлечься от решения задач, связанных с теплопереносом внутри слоя интрузивного расплава и ограничиться рассмотрением ситуации на его контактах с вмещающими породами. Это эквивалентно замене разбиения камеры на серию слоев (интервалов) только одним слоем, который характеризуется однородным распределением фазового состава и температуры. Вновь обратимся к Рис. 6.2. Из этого достаточно типичного графика видно, что процесс формирования строения пластового интрузива разбивается на три стадии, включающие: (1) образование зон закалки, (2) встречное движение верхнего и нижнего фронта кристаллизации с образованием приконтактовых зон и (3) остановку верхнего фронта кристаллизации с одновременным нарастанием мощности кумулуса в нижней части камеры. Интересно, что остановка верхнего фронта не относится к заранее заданным условиям, а проявляется как "ответ модели" на возможность оседания кристаллов. В этой асимметрии камеры (и температурного поля) заключается одно из важнейших отличий моделей Френкеля от других термических моделей интрузивов, не учитывающих перенос кристаллической фазы (Marsh, 1988; Worster et al., 1990; Mangan, Marsh, 1992) или связывающих процесс седиментации кристаллов с кинетикой зарождения твердой фазы в переходном слое (Simakin et al., 1994).

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

Взаимосвязь термики и динамики

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

, , (6.1)

где t - время, - теплоемкость (кал/г* град), - плотность (г/см³) и - температуропроводность (см²/сек) вмещающих пород, T₀ - температура на верхнем фронте кристаллизации. Выражения (6.1) представляют строгое решение задачи о прогреве полубесконечной среды, на границе которой поддерживается постоянная температура (Тихонов, Самарский, 1966). Возможность применения этих соотношений для описания зависимости теплового потока от времени обусловлена динамикой верхнего фронта кристаллизации, положение которого практически неизменно на главном этапе магматической дифференциации (Рис. 6.2). Для ранних моделей с упрощенной фазовой диаграммой и фиксированной температурой эвтектики (Рис. 6.1) эти уравнения являются довольно точной аппроксимацией теплового потока, рассчитанного при помощи сеточного разбиения. Для реальных базальтовых систем величина T₀ не является постоянной и понижается с течением времени.

На данном этапе мы не предпринимали попыток учесть этот фактор и скорректировать временную зависимость (равно как возможные вариации теплофизических параметров) и умышленно пошли на некоторое "искажение" истории остывания интрузива, сконцентрировав внимание на химизме процессов внутрикамерной дифференциации (см. ниже). Это замечание справедливо и в отношении ситуации на нижнем фронте, где используется более грубое приближение:

, , (6.2)

где - коэффициент пропорциональности, связывающий тепловые потоки на верхнем и нижнем фронте кристаллизации. Последующий опыт работы с программой ИНТРУЗИВ показал, что в ряде случаев величину необходимо корректировать в зависимости от времени (степени фракционирования расплава). Эти эмпирические поправки вводились на основании ранних (более строгих) термических моделей (Френкель, Ярошевский, 1978). Важно также заметить, что согласно конвекционно-кумуляционной модели тепловой поток на нижнем фронте обусловлен выделением теплоты кристаллизации захваченного (интеркумулятивного) расплава (кал/г):

(г/сек). (6.3)

Принципиально иная ситуация в прикровельной части камеры. Суммарный поток твердой фазы на верхнем фронте также зависит от теплоты кристаллизации расплава:

(г/сек), (6.4)

но для разных стадий дифференциации (Рис. 6.2) разделяется на поток направленной () и объемной () кристаллизации:

. (6.5)

Физический смысл этого разделения вытекает из специфики главных стадий: величина определяет скорость направленного роста верхней приконтактовой зоны (за счет кристаллов, механически связанных с кровлей); - скорость поступления кристаллического материала в объем конвектирующей магмы.

В качестве условного рубежа, отделяющего период формирования зон закалки от стадий направленного роста и эффективного оседания кристаллов было предложено ввести понятие времени закалки t^*(Френкель, 1982). На интервале 0 < t t* кристаллическая фаза в объем магмы не поступает и в процессах массопереноса не участвует (кристаллизация in situ), так что:

. (6.6)

При t > t*  доля материала, кристаллизующегося путем направленного роста, уменьшается пропорционально величине , соответственно нарастает поток объемной кристаллизации:

. (6.7)

Для нижнего фронта кристаллизации на любом временном интервале справедливо

. (6.8)

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

(6.9)

и общей долей кумулятивных кристаллов в нижней зоне ()

, (6.10)

где - плотность расплава (г/см³), - объемная доля взвешенных в магме зерен различных минералов, - доли кумулятивных минералов, S - площадь единичного сечения (см²). Соотношения (6.1-6.10) связывают термическую историю с динамикой продвижения фронтов кристаллизации в камере.

Взаимосвязь динамических и термодинамических параметров. Значения и средняя кристалличность являются экстенсивными параметрами, отражающими P-T условия существования смеси кристаллов и расплава в конвектирующем объеме магмы. Они определяются фазовой диаграммой системы для заданной фигуративной точки состава. Доля захваченных нижним фронтом кумулятивных кристаллов рассчитывается при помощи соотношений вида:

, (6.11)

где - скорость оседания кристаллика данного минерального вида. Уравнение (6.11) легко выводится из условия баланса масс для кристаллов, оседающих на продвигающейся вверх горизонтальной поверхности. Учитывая, что значения и по определению имеют разный знак, величина > , если > 0. Таким образом, соотношения (6.9-6.11) дают возможность рассчитывать фазовый состав пород, отвечающих этапу встречного движения верхнего и нижнего фронта кристаллизации.

Эти породы содержат интрателлурическую фазу (верхний фронт) или кумулятивные кристаллы (нижний фронт), однако называть их кумулятами не совсем корректно, поскольку на ранних стадиях количество кристаллизующегося in situ расплава >> (). По мере остывания интрузивного расплава в нем нарастает количество взвешенной фазы и замедляется скорость продвижения фронтов кристаллизации. Это приводит к увеличению концентрации кумулятивных кристаллов на нижнем фронте - именно эта особенность отвечает за формирование нижней "ветви" S-образной структуры некоторых минеральных и геохимических распределений в интрузивах.

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

Выражение для скорости продвижения нижнего фронта кумуляции кристаллов может быть получено суммированием уравнений (6.11) для всех кристаллических фаз, учитывая условие =:

. (6.12)

Очевидно, что при условии накопления в магматическом расплаве взвешенной кристаллической фазы (например за счет менее эффективной отсадки плагиоклаза) значение знаменателя в уравнении (6.12) понижается, а числителя - возрастает. В результате происходит ускорение фронта кумуляции, которое объясняет форму пунктирной линии на Рис. 6.2. Понятно также, что условие означает окончание процессов седиментации в системе. При расчете фазового состава "истинных" кумулятов значение в уравнении (6.11) заменяется на :

. (6.13)

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

. (6.14)

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

Моделирование внутрикамерной дифференциации с использованием подпрограммы ИНТРУЗИВ начинается с обычного расчета траектории равновесной кристаллизации расплава заданного состава (отвечающего исходной магме) по модели КОМАГМАТ - Рис. 2.28. Эти вычисления проводятся до заданной степени кристаллизации системы , которая соответствует предполагаемой доли интрателлурических вкрапленников в магме в момент внедрения . На этой стадии термодинамическая информация о равновесии (T, фазовые пропорции, составы кристаллов и расплава) впервые поступает в динамический блок ИНТРУЗИВ (Рис. 6.3): исходный объем 1 моля расплава увеличивается пропорционально заданной мощности интрузива h₀, определяются значения теплофизических характеристик расплава и вмещающих пород, начальная температура магмы на верхней границе T_o принимается равной рассчитанной температуре равновесия T.

Это дает возможность рассчитать тепловой поток из слоя магмы на ее верхнем и нижнем контакте и соответствующие суммарные потоки твердой фазы - см. уравнения (6.1-6.4). Отсюда для заданного времени закалки t^o можно рассчитать мощность закалочных зон (Френкель, 1982):

. (6.15)

Образование зон закалки приводит к сокращению объема магмы за счет "удаления" части расплава и интрателлурических фаз. После коррекции исходного объема информация из блока ИНТРУЗИВ поступает в главный цикл модели КОМАГМАТ для последующей кристаллизации (Рис. 6.3). Степень фракционирования системы смещается на величину .

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

. (6.16)

Более точное выражение учитывает также охлаждение магматического расплава по мере кристаллизации (Арискин, 1985).

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

. (6.17)

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

. (6.18)

Зная фазовый и химический состав кумулятов, не трудно рассчитать изменение валового состава магматической системы при формировании новых горизонтов пород. С учетом этой коррекции информация о текущем состоянии магмы после ее "частичной дифференциации" вновь поступает в блок решения задачи равновесия. Заметим, что система уравнений конвекционно-кумуляционной модели (6.1-6.18) допускает не только оседание, но также всплывание кристаллов с образованием кумулятов верхнего фронта - при этом применимы все приведенные выше соотношения, но с поправкой на изменение знака скорости седиментации и возможность "убегания" одной или нескольких фаз от фронта кристаллизации.

Таким образом, для каждого текущего значения времени положение верхней (h_U) и нижней (h_L) границы конвектирующего расплава задается соотношениями:

, (6.19)

которые описывают формирование внутренней структуры интрузива в зависимости от времени. Последовательное повторение всех действий, начиная с нулевого t при h_U (0) = h_o и h_L (0) = 0 есть решение динамической части построения модельного разреза интрузива. Эти итерации продолжаются до тех пор, пока поверхность кумулуса не достигнет верхнего фронта () или содержание кристаллов в магме не превысит заданного значения предельной кристалличности . При выполнении этих условий расчеты прекращаются. Рассчитанные вариации кумулятивных фаз и расплава в сочетании с эволюцией их составов определяют геохимическую структуру модельного интрузива.

Рис. 6.3. Схема включения динамического блока (программы ИНТРУЗИВ) при расчетах траектории кристаллизации магмы в комплексе КОМАГМАТ-3.0 (Ariskin et al., 1993) Детали термодинамической процедуры показаны на Рис. 2.28.

При помощи программы ИНТРУЗИВ можно строить модели строения конкретных интрузивных тел, включающие распределение по разрезу кумулятивных минералов, захваченного расплава, а также содержаний 10 главных и 20 примесных элементов. Эта работа является типичным примером "прямого" моделирования (forward modeling), при котором изучается "ответ" системы на принятые начальные условия и подгоночные параметры модели. Результаты расчетов сопоставляются с природными данными, затем в начальные условия и прогнозные параметры вносятся коррекции и вычисления повторяются снова до тех пор, пока не будет достигнуто (если это возможно в принципе) соответствие модельных распределений наблюдаемым.

К начальным условиям относятся мощность интрузива, состав исходной магмы, теплофизические свойства вмещающих пород (температура и плотность расплава рассчитываются по мере кристаллизации), давление и режим f_O2. Подгоночные параметры определяют полуэмпирический характер модели и включают: [1] долю интрателлурических вкрапленников в момент внедрения (), [2] значения "скорости теплопотерь" и , в том числе коэффициент пропорциональности - см. уравнения (6.1, 6.2), [3] время закалки t*, [4] предельную плотность кумулуса на нижнем и верхнем фронте ( - для случая всплывания кристаллов) и [5] эффективные скорости оседания минеральных индивидов .

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

Методические аспекты. В случае, если петрологу удалось подобрать набор начальных условий и параметров, при которых модель ИНТРУЗИВ реалистично воспроизводит строение интрузивного тела в комплексе минералогических и геохимических данных, есть основания говорить о правомерности приближения эффективной конвекции в камере, а некоторые из подгоночных параметров рассматривать как реальные физические величины.

Главный вопрос, который возникает в этой связи, касается однозначности решения задачи ЭВМ-моделирования. Может показаться, что число варьируемых в модели параметров слишком велико, так что одинаковые распределения главных и примесных элементов можно получить при различных сочетаниях начальных условий и подгоночных коэффициентов. Иными словами, в каждом случае необходимо доказывать устойчивость решения в зависимости от исходных данных. Методы анализа таких систем включают внесение случайных вариаций в исходные параметры и исследование воспроизводимости результатов моделирования (Шапкин, 1998). Применение данного подхода для модели ИНТРУЗИВ практически не возможно, поскольку при этом необходимо не только многократно варьировать 10-20 параметров, но также исследовать "отклик" модели одновременно по содержаниям 20-30 элементов.

Тем не менее, несколько моментов указывают на устойчивость решения задачи моделирования строения интрузивных тел. Во-первых, при поиске подгоночных параметров путем сопоставления модельных распределений с природными надо учесть, что петрохимические и геохимические данные представляют содержания порядка 30 элементов для нескольких десятков образцов. Таким образом, система аналитических данных сильно переопределена по отношению к числу варьируемых параметров модели. Второй момент - чисто практический. Личный опыт авторов по моделированию строения дифференцированных силлов показал, что независимо от прогнозных значений оптимальная модель распределения элементов по разрезу удовлетворяет устойчивой комбинации подгоночных параметров, которые варьруют в пределах 10-20 отн.%. Наконец, существует методика независимого контроля результатов прямого моделирования - геохимическая термометрия, которая позволяет определить некоторые из параметров конвекционно-кумуляционной модели, напр. среднюю кристалличность магмы и кумулуса на каждом уровне (см. раздел 3.4), и сравнить их с данными вариационных расчетов.

Первый опыт прямого моделирования строения дифференцированных силлов Кузьмовка (р. Подкаменная Тунгуска) и Палисейд (Нью-Джерси, США) был представлен в диссертации (Арискин, 1985). Параллельно Е.В.Коптев-Дворников и Б.С.Киреев провели работы по расчетам динамики формирования расслоенности силлов В-304 и Вавукан с верховьев р. Вилюй (Френкель и др., 1985). В последующем конвекционно-кумуляционная модель была использована для моделирования строения Талнахского интрузива (Днепровская и др., 1985). Вторая "волна" динамических вычислений захватывает конец 80-х годов, когда данные для Кузьмовского, Вилюйского (В-304) и Вавуканского интрузивов были пересчитаны с использованием скорректированной модели фазовой диаграммы (Арискин и др., 1986). Эти данные подробно описаны в коллективных публикациях (Френкель и др., 1988; Frenkel et al., 1989). Недавно были представлены результаты применения программы ИНТРУЗИВ для моделирования строения нижней зоны плутона Партридж Ривер (Chalokwu et al., 1996).

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

Вавуканский интрузив вскрывается в ряде обнажений по берегам р. Вилюй и его правого притока Вавукан. Ни в одном из обнажений интрузив не вскрыт полностью, однако в двух наиболее полных скальных выходах, отстоящих на расстоянии 1.5 км, были опробованы протяженные по вертикали фрагменты разреза (с шагом 5 м). Результаты геохимических исследований показывают, что эти фрагменты с перекрытием охарактеризовывают строение интрузивного тела. Это позволило составить полный сводный разрез Вавукана мощностью 100 м. По петрографическим данным верхняя и нижняя приконтактовые зоны представлены мелкозернистыми долеритами. В расслоенной серии снизу вверх сменяются горизонты пойкилоофитовых (Ol-Pl "кумулят") и такситоофитовых долеритов, габбродолеритов и феррогаббро (Ol-Pl-Aug± Mt).

В результате перебора более 50 вариантов начальных условий и подгоночных параметров, были найдены оптимальные условия вычислений, при которых модельные характеристики объекта во многих деталях совпали с наблюдаемыми (Рис. 6.4). Наиболее впечатляет соответствие данных для микроэлементов (Ni, Cr, Sc, V, Ba), которые представляют некогерентные и совместимые элементы. При этом каждый изгиб эволюционных линий или скачок содержаний находит объяснение как результат эволюции фазового состава кумулятов, включая долю захваченного в интеркумулусе расплава (Рис. 6.5). Важный результат этих работ заключается в установлении низкой исходной кристалличности вавуканской магмы ~ 2% и предельной кристалличности кумулуса 42-49 об.%. Оцененные эффективные скорости оседания минералов варьировали от 10 м/год (Pl) до 100 м/год (Aug).

Мы провели серию термометрических расчетов для пород Вавуканского интрузива по методике, апробированной ранее на долеритах Вельминского силла (Глава 3). Данные этой работы, полученные для основной части расслоенной серии, показаны на графиках Рис. 6.5 отдельными значками. Как видно, результаты решения прямой и обратной задачи физико-химического моделирования практически идентичны. Это сильный аргумент в пользу реалистичности конвекционно-кумуляционного механизма дифференциации и внутренней согласованности термодинамических и динамических параметров модели ИНТРУЗИВ. Заметим, что методики прямого и обратного моделирования независимы - общей их особенностью является только использование одной и той же модели фазовых равновесий.

Анализ показанных на Рис. 6.5 данных приводит к принципиальному выводу: однозначность (устойчивость) решения задачи моделирования проявляется в том, что существует единственное первичное распределение кумулятивных фаз по разрезу, которое описывает поведение всей совокупности главных и примесных элементов. Любые отклонения от этого распределения (скажем для Ol) немедленно приведут к искажениям концентрационных линий для других элементов, причем не только тех, которые совместимы с оливином (Mg, Ni, Fe, Co, Mn). Данный вывод представляется важным с методической точки зрения: при поиске оптимальной комбинации подгоночных параметров (прежде всего ) имеет смысл ориентироваться на фазовый состав кумулятов, установленный по данным геохимической термометрии. Это должно существенно ускорить поиск других оптимальных параметров модели.