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

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

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

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

(2.1)

при выполнении условий баланса вещества, где - фазы системы. Предложенные алгоритмы минимизации свободной энергии разделяются на две главные группы (Френкель и др., 1988; Борисов, Шваров, 1992): (I) программы первого типа строятся как процедуры непосредственной минимизации функции G на многограннике ограничений; (II) программы второго типа реализуют численные схемы решения систем нелинейных алгебраических уравнений, построенных по закону действующих масс.

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

Метод прямой минимизации термодинамического потенциала предпочтителен по многим причинам - он использует наиболее систематично организованную информацию о термодинамических свойствах веществ, допускает применение высоко развитой теории нелинейного математического программирования и обеспечивает максимальную универсальность ЭВМ-программ в смысле их применимости к разнообразным химическим системам, способам и диапазонам задания внешних условий. К этой группе алгоритмов относятся большинство программ, разработанных для термодинамического моделирования атмосферных, гидротермальных, метасоматических и метаморфических процессов (White, 1967; Волков, Рузайкин, 1974; Карпов и др., 1976; Шваров, 1976; Карпов, 1981; Harvie et al., 1987; Борисов, Шваров, 1992; Mironenko et al., 1992; Powell et al., 1997).

Алгоритмы этого вида нашли геохимическое применение в системах, для которых имеются надежные теоретические методы расчета химических потенциалов компонентов в зависимости от состава - не слишком плотные газы, подчиняющиеся теории Дебая-Хюккеля растворы сильных электролитов и т.д. В магматической петрологии по этому пути пошли авторы ЭВМ-программы MELTS (Ghiorso, 1985, 1994; Ghiorso, Carmichael, 1985), которые дважды перекалибровывали термодинамические параметры, описывающие свойства смешения силикатных расплавов на базе модели регулярных растворов (Ghiorso et al., 1983; Ghiorso, Sack, 1995). В нашей стране для расчета равновесий в магмах предпринимались попытки использовать алгоритмы И.К.Карпова (Феоктистов, 1983, 1987).

"Альтернативный" метод решения систем уравнений равновесия основан на использовании констант равновесия и стехиометрических коэффициентов для протекающих в системе реакций кристаллизации (1.25). При этом уравнения закона действующих масс выступают в качестве формального критерия равновесия и решаются совместно с уравнениями материального баланса. Способы решения этих уравнений основаны на разнообразных итерационных методах - метод Ньютона-Рафсона, метод оценки сдвига реакции и др. Алгоритмы, построенные на "константах", получили распространение в магматической петрологии (Langmuir, Hanson, 1981; Арискин, Френкель, 1982; Nielsen, Dungan, 1983; Френкель, Арискин, 1984аб; Ariskin et al., 1993; Nielsen, 1990; Weaver, Langmuir, 1990; Рябчиков, 1993; Camur, Kilinc, 1995; Ariskin et al., 1997)⁷ . В значительной мере это связано с тем, что для силикатного расплава, как "главной" растворной фазы, до сих пор отсутствуют надежные методы предсказания активностей компонентов, что затрудняет построение корректной поверхности свободной энергии в системе расплав-кристаллы (Гиорсоу, 1992).

Частично разрешить эту ситуацию позволяет использование эмпирических моделей активности компонентов в минералах и жидкости, которые определяют расчет "эффективных" констант равновесия минерал-расплав (Арискин, Френкель, 1982; Nielsen, Dungan, 1983; Weaver, Langmuir, 1990), см. уравнение (1.26). Замена активностей на концентрации в геотермометрах минерал-расплав (типа приведенных в Табл. 1.1) компенсируется тем, что заложенная в них связь между температурой и составами фаз ликвидусной ассоциации "снимается" непосредственно с экспериментальных данных и не включает никаких дополнительных пересчетов, связанных с использованием данных по термодинамике краевых компонентов силикатных систем - форстериту Mg₂SiO₄, анортиту CaAl₂Si₂O₈, хромиту FeCr₂O₄ и т.д.

Дело в том, что задача описания свободной энергии минеральных фаз в классической постановке минимизации функции G системы требует привлечения независимых термохимических данных для чистых миналов - стандартные значения энтальпий их образования, зависимости теплоемкости от температуры и т.п. В этом смысле модели решения задачи равновесия по алгоритмам первого типа должны быть согласованы с имеющимися базами данных по термодинамическим свойствам чистых веществ - например упомянутые выше модели Гиорсоу внутренне согласованы с базой данных (Berman, 1988). При отсутствии общей теории строения силикатных расплавов попытка такого согласования имеет и негативные последствия, поскольку в обработку экспериментальных данных по фазовым равновесиям с участием многокомпонентной силикатной жидкости извне привносится информация по краевым системам (Ghiorso et al., 1983; Ghiorso, Sack, 1995). Понятно, что при анализе непосредственной экспериментальной информации имеется больше шансов на удачное согласование данных по фазовым равновесиям, чем при условии наложения дополнительных, пусть даже термодинамически обоснованных, ограничений. В этом состоит определенное преимущество метода решения систем уравнений закона действующих масс: по сравнению с прямой минимизацией термодинамического потенциала он хотя и представляется излишне эмпиричным (Борисов, Шваров, 1992) или даже "примитивным" (Гиорсоу, 1992), но зато не требует детального знания зависимости химических потенциалов компонентов от составов минералов и расплава.

Забегая вперед, заметим, что специалисты в области магматической петрологии неоднократно проводили сравнительное тестирование алгоритмов кристаллизации, построенных на использовании "термодинамических" (MELTS - Ghiorso, Sack, 1995), полуэмпирических (COMAGMAT - Ariskin et al., 1993) и эмпирических моделей. Результаты такого тестирования оказываются обычно не в пользу программ, основанных на непосредственной минимизации термодинамического потенциала (Yang et al., 1996).

Таблица 2.1. Характеристики ЭВМ-моделей равновесной кристаллизации, для которых дано термодинамическое обоснование алгоритмов

Разработчики	Критерии	Активности компонентов		Примечания
программ	равновесия	Минералы	Расплав	;
Френкель, Арискин, (1984аб) Ariskin et al. (1993)	Закон действующих масс для составов + соответствие фазовых пропорций минимуму G при заданной степени кристаллизации системы	Однопозиционные модели идеального смешения. Минералы: Ol, Pl, Aug, Pig, Opx, Mt, Ilm	Двухпозиционные модели  смешения для компонентов: SiO2, TiO2, AlO1.5, FeO, FeO1.5, MnO, MgO, CaO, NaAlO2, KAlO2, PO2.5.	Программы серии КОМАГМАТ. Отклонения от идеальности в расплаве учитываются при помощи эмпирических коэффициентов при стуктурно-химических параметрах в уравнениях равновесия.
Ghiorso (1985, 1994) Ghiorso, Sack (1995)	Минимум свободной энергии G при заданной температуре	Модели регулярных растворов. Компоненты расплава: SiO2, TiO2, Al2O3, Fe2O3, MgCr2O4, Fe2SiO4, Mg2SiO4, CaSiO3, Na2SiO3, KAlSiO4, Ca3(PO4)2, H2O. Минералы: Ol, Pl, Aug, Opx, Sp, Ilm, Qtz, Lc, Ap, корунд и витлокит		Программы серии MELTS. Отклонения от идеальности в расплаве и минералах учитываются за счет параметров смешения регулярных растворов .
Langmuir, Hanson (1981) Weaver, Langmuir (1990)	Индекс насыщения (Reed, 1982) при сохранении баланса масс	Однопозиционные модели идеального смешения. Компоненты расплава: CaSiO3, NaAlO2, CaAl2O4, TiO2, FeO, MgO. Минералы: Ol, Pl, Aug		Отклонения от идеальности в расплаве никак не учитываются. Обновленные версии этой модели представлены в работах (Langmuir et al., 1992; Reynolds, 1995)
Camur, Kilinc (1995)	Ясно не сформулированы. В основе алгоритма - решение систем уравнений равновесия и баланса масс (минимизация G подразумевается).	Одно- и двухпозиционные модели идеального смешения. Минералы: Ol, Pl, Aug, Opx, Ne, Lc	Однопозиционная модель смешения. Компоненты расплава по модели Ghiorso et al., 1983	Отклонения от идеальности в расплаве учитываются путем введения для каждой реакции избыточной свободной энергии в виде функционала состава расплава, температуры и летучести кислорода.

При анализе данных Табл. 2.1, очевидно, следует пояснить отсутствие в ней информации по моделям кристаллизации, разработанным Нильсеном с коллегами (Nielsen, Dungan, 1983; Nielsen, 1990). Несмотря на популярность этих программ среди петрологов, в алгоритмах Нильсена отсутствует блок решения задачи термодинамического равновесия и построены они по принципу расчета траекторий фракционной кристаллизации, подробно рассмотренному в Главе 1 (Арискин, Френкель, 1982).

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

Задача ставится следующим образом. Предположим, что силикатная система известного валового состава (для простоты примем моль) произвольно разделена на некоторое количество минералов (j) и расплава, так что состав силикатной жидкости по каждому компоненту отвечает условиям материального баланса

, (2.2)

где - краевые компоненты твердых растворов (миналы), - стехиометрические коэффициенты реакций образования миналов (1.25), - количества отдельных миналов в смеси.

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

. (2.3)

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

. (2.4)

Откуда можно оценить влияние реакции образования каждого минала на изменение свободной энергии прогнозной расплавно-минеральной ассоциации

. (2.5)

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

,   , (2.6)

где значения и отвечают стандартным условиям чистых компонентов, и - составы исходного расплава и j-минерала. Соответственно, вычитаемая часть уравнения (2.5) может быть представлена как

. (2.7)

Из уравнений (2.5-2.7) следует, что

. (2.8);

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

. (2.9)

Два крайних правых члена в уравнении (2.8) представляют концентрационную константу , рассчитанную исходя из прогнозного состава силикатной жидкости и j-минерала:

. (2.10)

В результате, выражение (2.8) преобразуются к виду

, (2.11)

напоминающему уравнение изотермы реакции Ван-Гоффа, которое обычно применяется для анализа условий протекания самопроизвольных процессов в гомогенной системе (Николаев, 1979).

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

Например, учитывая выражения для температурной зависимости константы равновесия (1.26), уравнение (2.11) можно записать в виде

, (2.12)

или

, (2.13)

где T - температура искомого равновесия, а - температура метастабильного (прогнозного) равновесия с расплавом для данного минала. Уравнение (2.13) позволяет с терминов "химических потенциалов" (2.5) и "констант равновесия" (2.11) перейти к соотношениям температур, определяющих условия минимизации свободной энергии. Напомним, что R - газовая постоянная, а значения кратны энтальпии плавления чистого компонента и всегда являются величиной положительной, см. уравнения (1.19-1.24).

В случае, когда температура искомого равновесия T известна заранее, при помощи уравнения (2.13) можно предсказать, что при необходимо выделять из расплава r-минал j-минерала, а при переводить его в расплав (если количество этого компонента в смеси отлично от нуля), чтобы обеспечить понижение свободной энергии системы dG<0. Таким образом, расчет температур метастабильного равновесия с расплавом для каждого минала представляет важнейший элемент процедуры минимизации свободной энергии расплавно-кристаллической системы.

Если бы минералы представляли собой идеальные растворы компонентов, то эту задачу можно было бы решать абстрагируясь от того, что тот или иной минал входит в конкретную минеральную фазу, раздельно рассчитывая значения, например, и , а также молекулярные количества миналов ( , ). В действительности, отклонения от идеальности установлены для всех минеральных фаз, а использование эмпирических геотермометров (1.26), калиброванных на фазах переменного состава (1.19, Табл. 1.1), подразумевает необходимость оперировать при термодинамическом описании единым значением температуры для всех компонентов данного минерала .

Предложенный в работе метод оценки температуры метастабильного равновесия минерал-расплав имеет общие особенности с методикой расчета псевдоликвидусных температур (1.28), но отличается тем, что учитывает присутствие в системе помимо расплава определенных количеств минералов. Важная особенность этого алгоритма состоит в том, что он ориетирован на решение задачи равновесия не при заданой температуре, а при фиксированной степени кристаллизации системы, когда суммарное количество молей компонентов расплава не меняется в процессе кристаллизации и растворения минералов ; (Френкель, Арискин, 1984аб). Это несколько осложняет процедуру минимизации функции G, зато дает возможность разрабатывать ЭВМ-программы c фиксированным шагом кристаллизации (Ariskin et al., 1993) или плавления (Ariskin et al., 1997) в системе. Значение этой возможности смогут оценить специалисты, заинтересованные в расчетах равновесного, фракционного и критического плавления мантийных материалов, когда температура солидуса системы заранее не известна.

Вернемся к начальному прогнозному состоянию: силикатная система включает молей расплава состава и определенные количества минералов состава :

, (2.14а)

и

, (2.14б)

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

. (2.15)

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

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

, (2.16)

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

Система уравнений вида (2.16) становится разрешимой относительно температуры при использовании дополнительного условия стехиометрии данной минеральной фазы (1.27). При этом уравнение "смещенного" равновесия для одного их миналов (2.16) необходимо заменить на очевидное условие материального баланса в твердой фазе:

или , (2.17)

которое гарантирует постоянство количества минерала в смеси и, соотвественно, .

В качестве примера решения подобной задачи можно привести систему уравнений равновесия и баланса масс, использованных для уравновешивания заданного количества плагиоклаза и комплементарного расплава в программе РТРК (Арискин, 1985)⁸ :

, ,
,
,
,    (2.18)
,
,
,

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

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

, (2.19)

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

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

Дело в том, что каждый шаг кристаллизации или плавления минералов фиксированного состава (2.16), хотя и приближает систему к равновесному распределению фазовых пропорций, но неизбежно нарушает выполнение условий (1.29). Это нарушение проявляется в том, что изменение пропорций минералов не сопровождается изменением их состава, хотя состав жидкой фазы при этом эволюционирует (представим выделение оливина и комплементарное растворение плагиоклаза). Разуравновешивание составов вынуждает повторно обращаться к процедуре поиска температур метастабильного равновесия и новых составов минералов. Таким образом, основной алгоритм минимизации свободной энергии частично расплавленной системы (2.13) строится как совокупность двух главных итерационных циклов (Рис. 2.1): цикл Термометрии расплавно-минеральных равновесий переводит систему в метастабильное состояние, где состав каждой твердой фазы связан с составом жидкости законом действующих масс; цикл Расчета фазовых пропорций минимизирует разницу между температурами равновесия отдельных минералов путем их кристаллизации или плавления.

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

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

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

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

На ВХОДЕ в программу РТРК необходимо также инициализировать параметры геотермометров минерал-расплав и (1.26). Эти параметры могут быть постоянны или переопределяться в процессе вычислений, как это будет показано на примере моделей полибарического фракционирования (см. Главу 4). При известном составе исходного расплава и заданных значениях и в блоке ТЕРМО рассчитываются псевдоликвидусные температуры равновесия для всех минералов (1.28) и соотвествующие составы минералов (1.29). После этого программа находит минеральную фазу j1 , отвечающую температуре ликвидуса и в соотвествии с начальным условием считает задачу равновесия решенной при (Рис. 2.2). Выход из блока ИНФОРМ при условии < переводит выполнение программы в блок расчета пропорций кристаллизации и точности итерационных циклов на следующей стадии вычислений.

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

, (2.20)

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

Точность расчета температур и составов необходимо переопределять на каждом этапе решения задачи равновесия, поскольку от этих параметров зависит сходимость комбинированных итерационных циклов (Рис. 2.2). Эту математическую точность вычислительной схемы следует отличать от точности решения задачи равновесия, которая определяется заданием интегральных внешних параметров и (Рис. 2.1 и 2.2): при шаге кристаллизации =1-2 мол.% оптимальные значения составляют =0.5-1.0^оС; для точности расчета состава жидкой фазы мы обычно принимаем =0.1-0.2 мол.%¹⁰ . Эти неопределенности удовлетворяют решению большинства петрологических и геохимических задач, на которые ориентировано использование данных моделей. Проблема состоит в том, что интегральные погрешности не могут варьироваться произвольным образом, поскольку имеют ограничения, связанные с точностью решения задач метастабильного равновесия в отдельных итерационных циклах.

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

, , . (2.21)

Соотношения (2.15, 2.21) дают уравнение

, (2.22)

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

Точность расчета состава минерала определяется точностью оценки температуры для каждой твердой фазы - : эти параметры связаны через константу равновесия (1.29). На начальных стадиях кристаллизации эта зависимость может быть оценена при помощи уравнения (1.26):

, (2.23)

. (2.24)

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

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

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

Примечания: