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

Представленный в предыдущих разделах алгоритм моделирования кристаллизации силикатных расплавов, база данных по фазовым равновесиям и процедуры калибровок уравнений равновесия составляют систему ЭВМ-моделирования кристаллизующихся магм. Важным элементом этой системы являются поисковые опции программ ИНФОРЭКС, которые объединены в едином интерфейсе с процедурами расчета регрессионных коэффициентов для геотермометров минерал-расплав. Эти коэффициенты записываются в выходной файл, формат которого позволяет непосредственно стартовать программу расчета траекторий кристаллизации. Такая организация работы дает возможность калибровать модели для заданного диапазона составов магм и P-T-f_O2 уcловий.

Для верификации этих моделей программа ИНФОРЭКС предоставляет возможность создавать файлы вида минерал-расплав (Ol-l, Pig-l, Mt-l и т.д.), включающие данные по температурам и составам стекол, закаленных в присутствии данной минеральной фазы, даже если отсутствуют данные по составам этого минерала. Этот массив независимой от калибровок информации может быть использован для дополнительного тестирования и коррекции интегральных расплавно-минеральных геотермометров. Другой важный способ "настройки" и тестирования моделей фазовых равновесий задействует данные по валовым составам экспериментальных материалов. Используя эти составы в качестве исходных для моделирования, не трудно рассчитать траектории равновесной кристаллизации в условиях проведения экспериментов и сравнить их с наблюдаемыми. При наличии систематического сдвига температур появления на ликвидусе одной из фаз, модели серии КОМАГМАТ позволяют ввести температурную поправку, которая составляет обычно не более 10^оС и находится в рамках точности разработанных геотермометров.

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

Программный комплекс КОМАГМАТ предназначен для расчета фазовых равновесий и моделирования процессов кристаллизационной дифференциации магм широкого диапазона составов, включающих магнезиальные и высокоглиноземистые базальты, разнообразные толеиты и андезит-базальтовые расплавы. Первая версия этой системы (КОМАГМАТ-1.n) была разработана в 1989 г. для использования на машинах серии ЕС. В начале 90-х годов появились DOS-программы, предназначенные для IBM совместимых персональных компьютеров (версии КОМАГМАТ-2.n). В 1993 г. они получили законченный облик и были представлены в виде системы КОМАГМАТ-3.0, которая в настоящее время широко распространена среди специалистов в области магматической петрологии (Ariskin et al., 1993; Ariskin, Nielsen, 1993).

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

1. Термометрия равновесий минерал - расплав (Thermometry of Mineral-Melt Equilibria). Этот режим позволяет рассчитывать температуры насыщения магматических расплавов минеральными фазами при заданных внешних условиях (P-T-logf_O2 и состав расплава, включая содержание H₂O). Рекомендуется для оценки температур кристаллизации природных закалочных стекол (расплавных включений) и тестирования встроенных в КОМАГМАТ геотермометров при помощи файлов минерал-расплав (Ol-l, Pig-l, Mt-l), созданных в базе ИНФОРЭКС.

2. Моделирование равновесной кристаллизации (Simulating Equilibrium Crystallization) предполагает расчет линий эволюции состава расплава и минералов, а также изменение фазовых пропорций при кристаллизации в условиях сохранения валового состава системы. Этот режим используется для дополнительной настройки и коррекции модели КОМАГМАТ путем сравнения рассчитанных и экспериментальных порядков кристаллизации минералов (см. ниже). Расчет траекторий равновесной кристаллизации расплавов природных образцов является основой метода Геохимической Термометрии, который будет рассмотрен в Главе 3.

3. Моделирование фракционной кристаллизации (Simulating Fractional Crystallization). Этот режим позволяет рассчитывать линии эволюции химического и фазового состава системы при идеальном фракционировании, обеспечивающем полное отделение кристаллизующихся фаз от расплава. Предназначен для исследования гипотезы формирования ассоциации магматических пород путем фракционирования некоторого источника - первичной или родоначальной магмы. В ряде случаев, сравнение линий фракционной кристаллизации и наблюдаемых в породах геохимических трендов приводит к оценкам условий протекания кристаллизационного процесса, в частности - давления и летучести кислорода (см. разделы 4.3 и 4.4).

4. Моделирование формирования расслоенных интрузивов (Simulating Layered Intrusion Formation). Этот блок вычислений включает алгоритм моделирования динамики тепломассопереноса и химической дифференциации в магматических камерах, которые могут быть аппроксимированы формой пластового интрузива (Френкель и др., 1988). Примеры использования этой программы будут даны в Главе 6. Отметим также, что данный режим позволяет моделировать случаи частичного и промежуточного фракционирования магм (см. раздел 1.3).

Условия вычислений. Задание внешних условий в программе КОМАГМАТ сводится к определению барического и окислительно-восстановительного режима (Рис. 2.27). Для указания барических условий необходимо выбрать (а) режим изобарической кристаллизации P=const (Isobaric crystallization), (б) режим увеличения давления при проведении серии изобарических расчетов P=P+D P (Increasing pressure routine) или режим декомпрессионной кристаллизации (Decompression crystallization), моделирующий условия сброса давления по мере кристаллизации исходного расплава P(j )=P(j )-D P.

При задании редокс-условий предусмотрено 2 главные ситуации: (а) закрытость системы по кислороду (Closed system (Fe2+/FeO) ratios) и (б) условия буферированности по давлению O₂ - Open system (oxygen buffers). В первом случае пользователю необходимо задать общую степень восстановленности железа в расплаве Fe(II)=Fe²⁺/(Fe³⁺ + Fe²⁺); обычно эта величина лежит в диапазоне 0.80 =< Fe(II) =< 0.95. Во втором - указать один из 12 кислородных буферов, традиционно используемых при петрологических интерпретациях (QFM, WM, NNO и др.) или задать постоянное давление кислорода logf_O2=const.

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

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

Сетевая версия программы КОМАГМАТ. Широкое распространение программы КОМАГМАТ-3.0 среди петрологов, обновление локальных DOS-версий и внедрение современных коммуникационных технологий вызвали необходимость создания сетевой версии программы, доступной пользователям через ИНТЕРНЕТ. Для решения этой задачи мы объединили усилия с коллегами кафедры петрологии МГУ, занимающимися обеспечением сетевого доступа к информационным ресурсам в области наук о Земле (рук-ль работ П.Ю.Плечов). В результате этой работы из локальной версии КОМАГМАТ-3.0 удалось вычленить расчетное ядро, которое получает из интерфейсной программы поток данных, обрабатывает их и выдает результат обратно в интерфейсную программу (Рис. 2.29). Сетевая версия ЭВМ-модели КОМАГМАТ-3.5 (1997) доступна теперь по электронному адресу http://www.geo.web.ru/users/comagmat и позволяет проводить расчеты фазовых равновесий различных магматических систем, работая в режиме удаленного доступа.

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

Представление об эффективности комплекса КОМАГМАТ и точности численного моделирования фазовых равновесий дает сравнение расчетных и экспериментальных траекторий равновесной кристаллизации. Для этой цели могут быть использованы данные melting-экспериментов с образцами известного состава, когда смеси выдерживаются при температурах, отвечающих преобладанию расплава над твердой фазой. В этих опытах температура последовательно понижается от ликвидусной до значений, при которых еще удается получить гомогенные закалочные стекла (обычно до 50-70% кристаллов).

К сожалению, у этого метода существует серьезный недостаток - при проведении экспериментов практически никогда не удается сохранить неизменным валовый состав системы. Большинство опытов при атмосферном давлении проводится "методом платиновой петли" в токе газовой смеси, задающей давление кислорода. В условиях высоких температур эта методика приводит к потери летучих элементов (Na, K, P), которая при длительности плавок в сотни часов может достигать десятков относительных процентов (напр., Grove et al., 1982). В ряде случаев не удается избежать и потери FeO в платину (до нескольких отн.%). Это обстоятельство оказывается весьма существенно, поскольку условие сохранения баланса масс в программе КОМАГМАТ является принципиальным моментом, определяющим температуры кристаллизации, пропорции и составы минералов по мере увеличения степени кристалличности системы.

Учитывая этот фактор, при сопоставлении результатов расчетов и экспериментов необходимо акцентировать внимание на ранних и средних стадиях кристаллизации, уделяя главное внимание последовательности появления минералов на ликвидусе. Первые результаты такого сравнительного анализа были даны в работе (Арискин и др., 1986), где рассмотрены экспериментальные и модельные последовательности кристаллизации для 9 расплавов, преимущественно толеитового состава. При проведении этих вычислений использовалась программа, которая позднее явилась основой версии КОМАГМАТ-1.n (см. выше). В этом разделе монографии будут приведены результаты новых расчетов, основанных на использовании одной из последних версий КОМАГМАТ-3.51, включающей уточненные геотермометры Mt-расплав, Ilm-расплав и некоторые дополнительные модификации (Ariskin, Barmina, 1999).

Исходные данные. Для моделирования траекторий равновесной кристаллизации было выбрано 10 исходных составов из 10 различных экспериментальных исследований, проведенных за период 1979-1996 гг. (Табл. 2.9 и 2.10). Данные (Walker et al., 1979; Grove et al., 1982; Grove, Bryan, 1983; Tormey et al., 1987) были задействоваваны ранее при калибровке модели КОМАГМАТ-3.0, остальные работы представляют независимые эксперименты. 6 выбранных исходных составов относятся к различным толеитовым базальтам (Табл. 2.9), 4 состава включают высокоглиноземистые и субщелочные базальтоиды (Табл. 2.10).

Табл. 2.9. Составы экспериментально изученных толеитовых базальтов, которые использовались при тестировании программы КОМАГМАТ

Условия вычислений. Расчет траекторий равновесной кристаллизации проводился при помощи режима Simulating Equilibrium Crystallization (Рис. 2.27) в условиях, максимально приближенных экспериментальным: давление - 1 атм, летучесть кислорода соответсвует буферному равновесию QFM. Шаг кристаллизации составлял 1 мол.%, при этом максимальное значение степени кристаллизации расплава варьировало от 70 до 80%. Рассчитанные и экспериментальные последовательности кристаллизации минералов приведены на Рис. 2.30 и 2.31. Эти графикы построены в диапазоне температур 1050-1300^оС, который отвечает условиям калибровки базовых геотермометров минерал - расплав (см. раздел 2.3).

Рис. 2.30 Сравнение экспериментальных и расчетеных   температур кристаллизации минераорв при равновесном  затвердевании расплавов толеитовых базальтов. Номера графиков отвечают номерам исходных составов в  табл. 2.9. Расчеты проводились с использованием программы   КОМАГМАТ-3.51 (1997).

Сравнение порядка кристаллизации. Данные Рис. 2.30 показывают, что модель КОМАГМАТ с хорошей точностью воспроизводит последовательность выделения минералов в толеитовых системах: Ol Pl-> Aug-> Pig-> оксиды (Mt Ilm). При этом отклонения расчетных температур начала кристаллизации минералов от экспериментальных обычно не выходят за пределы 10-15^оС. Этот результат отражает тот факт, что калибровка базовой модели проводилась на данных, большинство которых получено в системах толеитового состава. Особое внимание заслуживает 5-й график, показывающий, что программа адекватно "реагирует" на насыщение остаточного расплава Fe-Ti оксидами.

Большее разнообразие демонстрируют данные для высокоглиноземистых и субщелочных базальтов (Рис. 2.31). В 3-х случаях здесь ликвидусная фаза представлена Pl, шире поля кристаллизации оксидных минералов. Модельные температуры кристаллизации Aug и Ilm несколько завышены, но в целом рассчитанная последовательность появления и набор ликвидусных фаз отвечают результам экспериментов.

Табл. 2.10. Составы высокоглиноземистых и субщелочных базальтоидов, использованных при тестировании программы КОМАГМАТ

Для последующего использования программы КОМАГМАТ важно отметить, что корректность модельных линий эволюции состава расплава определяется не столько абсолютными, сколько относительными температурами кристаллизации минералов. В этом смысле данные экспериментов и расчетов на Рис. 2.30(2) практически идентичны и могут давать на вариационных диаграммах хорошее соответствие модельных и экспериментальных трендов - см. Рис. 2.34 и 2.35.

Пропорции кристаллизации. При анализе результатов расчетов с использованием программы РТФК отмечалось (см. разделы 1.3 и 1.4), что направленность петрохимических и геохимических трендов определяется соотношениями между кристаллизующимися минералами; погрешности в оценках составов минералов оказываются менее существенны. Учитывая, что в модели КОМАГМАТ пропорции кристаллизации не задаются, а рассчитываются по разности количества фаз в предыдущем и текущем состоянии (2.20), контроль корректности этих соотношений имеет принципиальное значение. При разработке ранних версий модели КОМАГМАТ эмпирическая база для подобных сравнений была весьма ограничена: нам удалось найти лишь одну работу, где при помощи оптических методов в продуктах экспериментов прослеживалось изменение количеств Ol и Pl (Bender et al., 1978). Сопоставление этих экспериментальных данных и расчетных фазовых пропорций убедило, что разработанные модели с хорошей точностью предсказывают соотношения между кристаллизующимся Ol и Pl - примерно 3:7 (Френкель, Арискин, 1984б; Френкель и др., 1988).

Рис. 2.32. Изменение фазового состава системы при равновесной кристаллизации  расплава толеитового базальта Результаты моделирования отвечают данным графика 4, показанном на Рис. 2.30 (исходный состав 4 в Табл. 2.9).

В 80-х и 90-х годах появилось большое количество новых экспериментальных данных, позволяющих оценивать пропорции кристаллизации главных породообразующих минералов. Они представляют результаты математической обработки информации по составам минеральных фаз и расплава, закаленных в системах известного валового состава. Основой для этих оценок являются расчеты по методу МНК, которые проводятся для каждого эксперимента аналогично масс-балансовым вычислениям типа источник - дифференциат (Wright, Doherty, 1970) - см. раздел 1.3. Пример эволюции фазового состава системы, рассчитанной таким образом по экспериментальным данным (Grove et al., 1990), показан на Рис. 2.32.

Данные Рис. 2.32 построены в зависимости от степени кристаллизации ,которая комплементарна доли остаточного расплава, рассчитанной одновременно с оценками относительного количества Ol, Pl и Aug (). Такое построение удобно для иллюстрации метода расчета пропорций кристаллизации. Очевидно, что в выбранных координатах пропорции кристаллизации представляют наклон зависимости для каждого j -минерала. Учитывая неопределенности, связанные с погрешностями оценок составов синтезированных фаз и проблемами сохранения валового состава системы, более информативным является расчет не "точечной", а средней пропорции кристаллизации , где отвечает интервалу устойчивости данной 2-х, 3-х или 4-х фазовой минеральной ассоциации.

Используя этот подход, мы обработали данные по фазовым пропорциям, представляющим данные опытов на 18 образцах из 10 независимых исследований. Некоторые из этих результатов суммированы на Рис. 2.33. Здесь показано только две серии данных, включающих "экспериментальные" пропорции кристаллизации для ассоциаций Ol-Pl и Ol-Pl-Aug. Для сравнения приведены результаты расчетов по модели КОМАГМАТ для образца океанического толеита (Grove et al., 1990) - см. график 4 на Рис. 2.30 и Рис. 2.32.

Показанные на Рис. 2.33 данные позволяют сделать два вывода. Во-первых, модельные пропорции кристаллизации попадают в диапазон экспериментально установленных значений, что еще раз подтверждает адекватность модели КОМАГМАТ¹³ . С другой стороны, заметные колебания "точечных" значений указывают на неустойчивость вычислительной схемы, которая возрастает при смене кристаллизующейся ассоциации Ol-Pl на Ol-Pl-Aug. Причина этой неустойчивости заключается в сильном влиянии математической точности вычислений (см. раздел 2.1) на фазовые пропорции. При анализе этого вопроса выяснилось, что рекомендуемая нами обычно точность итерационных циклов порядка 1^оС (Рис. 2.27) хотя и гарантирует стабильность интегрального результата ЭВМ-моделирования, далеко не всегда дает устойчивое решение на небольших интервалах кристаллизации.

Таблица 2.11. Средние пропорции кристаллизации минералов при равновесном затвердевании природных силикатных расплавов (по результам ЭВМ-моделирования)

В Табл. 2.11 сведены данные по средним пропорциям кристаллизации для минеральных ассоциаций, представляющих траектории, рассчитанные при помощи программы КОМАГМАТ и показанные на Рис. 2.30 и 2.31. Не трудно заметить, что для одних и тех же котектик эти значения несколько варьруют, но разделить определенно эффект влияния состава расплава и неустойчивости математического решения в настоящее время не представляется возможным.

Эволюция состава расплава. Хорошее воспроизведение модельных пропорций кристаллизации позволяет надеяться, что при условии корректного расчета последовательности появления минералов зависимость состава остаточного расплава от температуры может быть смоделирована достаточно точно. Сравнение модельных и экспериментальных линий эволюции состава расплава показано на Рис. 2.34. Для тестирования были выбраны результаты экспериментов на 2-х образцах толеитовых базальтов - см. составы 2 и 3 в Табл. 2.9 (Рис. 2.30, графики 2 и 3). В обоих случаях результаты вычислений хорошо воспроизводят составы экспериментальных стекол, причем все изгибы эволюционных линий находят объяснение как результат смены фазового состава кристаллизующихся ассоциаций и состава минералов. Расчет температурно-композиционных зависимостей для жидкой фазы лежит в основе методик Геохимической Термометрии и Редоксиметрии магматических пород, которые будут рассмотрены в последующих главах.

Другая форма представления эволюционных линий - вариационные диаграммы. На Рис. 2.35 показаны те же результаты расчетов и экспериментов, однако составы жидкой фазы приведены в зависимости от содержания в расплаве MgO. Структура модельных трендов при этом принципиально не изменилась - см. графики MgO-T на рис. 2.34. Эти данные свидетельствуют, что программа КОМАГМАТ может быть использована в традиционном петрологическом анализе, включая исследование вопроса о возможности формирования петрохимических трендов в кристаллизационных процессах.

Рис. 2.35. Сравнение экспериментальных и модельных линий эволюции состава остаточного расплава при равновесном затвердевании расплавов толеитовых базальтов Исходные составы отвечают номерам 2 и 3 в Табл. 2.9. Расчеты проводились при помощи программы КОМАГМАТ-3.51 (1997). Порядок кристаллизации - см. Рис. 2.30.

В целом, сравнительный анализ экспериментальных данных и результатов расчетов при атмосферном давлении показывает, что низкобарная версия модели КОМАГМАТ реалистично предсказывает фазовые равновесия для различных сочетаний породообразующих минералов. Более проблематично моделирование высокобарной кристаллизации.

Иллюстрация "высокобарных возможностей" ЭВМ-модели КОМАГМАТ проведена на примере построения P-T диаграммы фазовых отношений магнезиального толеита (Рис. 2.36). Для расчетов был выбран исходный состав, представляющий базальты Срединно-Атлантического хребта с содержанием ~ 10 мас.% MgO. Этот магнезиальный толеит был использован в серии экспериментов при давлениях 1 атм - 15 кбар (Bender et al., 1978), которые послужили основой для калибровки высокобарной версии модели КОМАГМАТ - см. раздел 2.3 (Рис. 2.24).

Условия вычислений. Расчеты траекторий равновесной кристаллизации проводились при помощи режима Simulating Equilibrium Crystallization с использованием опции последовательного увеличения давления в серии изобарических расчетов (Increasing pressure routine) - Рис. 2.27. Окислительно-восстановительные условия отвечали экспериментальным (буфер железо-вюстит), а шаг кристаллизации составлял 1 мол.%. Вычисления прекращались при степени кристаллизации расплава 60%. В интервале давлений 1 атм - 10 кбар изобарические расчеты проводились с шагом = 2 кбар; выше 10 кбар сетка вычислений была сужена до 0.5 кбар, чтобы более определенно зафиксировать ожидаемую смену порядка появления Ol и Aug на ликвидусе. Такая процедура позволила сконструировать модель P-T диаграммы, которая наглядно воспроизводит общие особенности смены фазовых ассоциаций базальтов при изменении давления.

Контроль точности. Данные фазовой диаграммы на Рис. 2.36 показывают, что в условиях умеренных давлений 3-7 кбар предложенная нами модель удовлетворительно предсказывает порядок кристаллизации Ol-> Pl+Aug (или Aug+Pl). Это подтверждают и результаты расчета кристаллизации расплава другого толеитового базальта, проведенные по программе КОМАГМАТ для давления 4 кбар (Yang et al., 1996). Такое соответствие представляется важным в связи с практическим использованием модели высокобарных фазовых равновесий для целей барометрии условий образования магматических расплавов, представленных среди закалочных стекол Мирового океана (см. раздел 4.4). При давлениях выше 10-12 кбар смещение модельных фазовых полей относительно экспериментальных данных может достигать 3-4 кбар (Рис. 2.36). Этот фактор необходимо учитывать при интерпретации данных высокобарических расчетов, в частности - касающихся декомпрессионного фракционирования островодужных высокомагнезиальных магм ( Глава 5).

Главные проблемы применения высокобарной версии программы КОМАГМАТ связаны с оценкой стабильности ортопироксена. На данном этапе разработки модели появление этого минерала и перитектические соотношения с оливином даже при 1 атм предсказываются не слишком надежно. В случае высокобарных расчетов, проведенных с магнезиальными базальтами, эти неопределенности усиливаются расхождениями разных авторов в идентификации поля устойчивости Opx при давлениях до 10-12 кбар (Bender et al., 1978; Green et al., 1979; Fujii, Bougault, 1983; Draper, Johnston, 1992). По этой причине при проведении высокобарных вычислений для природных объектов в дальнейшем мы ограничимся рассмотрением систем, для которых фракционирование Opx существенной роли не играло.

Попытки перекалибровки программы КОМАГМАТ, анализа и обработки данных высокобарных опытов приводят к пессиместичному выводу, что качество экспериментальных данных, полученных при высоких давлениях, заметно уступает данным при 1 атм. Это касается степени равновесности продуктов закалочных опытов, точности оценки давления и температуры "равновесия", а также неопределенности режима f_O2. По-видимому, эти данные позволяют разработать ЭВМ-модели, средняя точность которых не превышает 20-30^оС. Определенный скепсис в этом отношении высказывают и другие разработчики моделей фазовых равновесий (Ghiorso, Sack, 1995; Ghiorso, 1997).

При моделировании фазовых равновесий в водосодержащих системах мы рассматриваем воду как аналог некогерентного компонента, полностью растворенного в расплаве. Вычисления "в присутствии воды" организованы таким образом, что H₂O выступает в качестве фиктивного компонента, концентрация которого не входит в расчеты констант равновесия реакций кристаллизации и не оказывает влияния на активности минералообразующих компонентов. Такой подход предполагает, что в процессах кристаллизации "вода" монотонно накапливается в жидкости, силикатная часть которой эволюционирует по законам сухой системы. Накопление "воды" в расплаве может происходить вплоть до концентраций, отвечающих пределу растворимости при данном давлении. При этом влияние H₂O на фазовые равновесия учитывается опосредованно - через изменения относительных температур кристаллизации минералов - см. уравнение (2.49) и Рис. 2.26.

Алгоритм моделирования. Содержание воды в исходном расплаве задается при помощи основного интерфейса (Рис. 2.27) и постоянно находится в памяти машины. В соответствии с предложенным алгоритмом, на каждой стадии процесса сначала рассчитываются прогнозные температуры кристаллизации минералов для "сухого" расплава, а затем они корректируются в сторону уменьшения - пропорционально коэффициентам, приведенным в Табл. 2.5, и модельной "концентрации" H₂O (2.49). При проведении вычислений с использованием "водной" версии программы КОМАГМАТ расчетное содержание воды в жидкой фазе после каждого шага кристаллизации сравнивается с предельной растворимостью, оцененной по уравнению (2.51). Если модельное содержание H₂O выше растворимости, вычисления могут быть остановлены или продолжены, но при условии отвода из расплава на каждом шаге 3-5% относительных процентов воды. Последняя ситуация моделирует дегазацию водонасыщенного магматического расплава и эволюцию фазового состава системы под давлением паровой фазы. Этот примитивный подход к расчету расплавно-минеральных равновесий в присутствии воды следует рассматривать в качестве первого приближения, которое, однако, не плохо воспроизводит известные закономерности смещения фазовых полей в водосодержащих системах (Альмеев, Арискин, 1996).

На Рис. 2.37 показаны данные ЭВМ-моделирования фазовых равновесий при кристаллизации расплава высокоглиноземистого базальта в сухих и водосодержащих условиях. Исходный состав отвечают образцу 79-35g, который был использован в серии водонасыщенных экспериментов при давлении 2 кбар (Sisson, Grove, 1993). Моделирование условий насыщения расплава H₂O проводилось следующим образом: к исходному составу ВГБ "добавлялось" заданное количество воды (1.0, 1.5, 2.0 мас.% и т.д.) и для каждой системы рассчитывалась траектория равновесной кристаллизации с шагом =1 (P=2 кбар, буфер QFM). Таким образом, для составов с различным исходным содержанием воды были найдены параметры равновесия расплав-минералы-пар, включающие степень кристаллизации расплава , температуру равновесия T_j, пропорции минералов на ликвидусе и состав жидкой фазы.

В результате было установлено, что при исходных содержаниях воды 2-2.5 мас.% модельные траектории кристаллизации реалистично воспроизводят экспериментальные температуры (1020-1070^оС) и концентрацию H₂O (5.5-6.0 мас.%) в водонасыщенных расплавах (Альмеев, Арискин, 1996). При этом на фоне общего понижения ликвидусной температуры наблюдается заметное расширение поля фемических минералов и сокращение поле кристаллизации плагиоклаза относительно данных, полученных для сухих систем. Таким образом, модельные линии эволюции состава высокоглиноземистого расплава отвечают известным представлениям о задержке кристаллизации плагиоклаза и накоплении Al₂O₃ в остаточной жидкости в присутствии воды.

Примечание:

¹³ Недавние результаты независимого тестирования программы КОМАГМАТ и сравнения с другими моделями показали, что она наиболее точно воспроизводит пропорции кристаллизации силикатных минералов при атмосферном давлении (Yang et al., 1996).