Петрология, 2024, № 1

Российская академия наук
ПЕТРОЛОГИЯ
Том 32     № 1     2024     Январь–Февраль
Основан в 1993 г. 
Выходит 6 раз в год 
ISSN 0869-5903
Журнал индексируется 
в Current Contents
Журнал издается под руководством 
Отделения наук о Земле РАН
Главный редактор
В.В. Ярмолюк
Редакционная коллегия:
Акинин В. В., Аранович Л. Я., Арискин А. А., 
Дубинина Е. О., Каменецкий В. С.,
Каргин А. В., Копылова М. Г., Котов А. Б., 
Латыпов Р. М., Носова А. А. (ответственный секретарь),
Плечов П. Ю., Портнягин М. В., Пухтель И. С., 
Самсонов А. В. (заместитель главного редактора), 
Сафонов О. Г., Силантьев С. А., Симакин А. Г., 
Скляров Е. В., Соболев А. В.
Зав. редакцией И.И. Невская
Адрес редакции: 119017 Москва, Старомонетный пер., 35 
e-mail: petrolog@igem.ru
Москва
ФГБУ «Издательство «Наука»
© Российская академия наук, 2024
© Редколлегия журнала 
     “Петрология” (составитель), 2024

СОДЕРЖАНИЕ
Том 32, номер 1, 2024
“В кильватере большого корабля”: плаванье продолжается!
Л.Я. Аранович
3
Моделирование фильтрации многокомпонентного флюида в деформируемых 
вмещающих породах с учетом химического взаимодействия
Л.А. Хакимова, Ю.Ю. Подладчиков
5
Тектониты приенисейской сдвиговой зоны (Енисейский кряж): свидетельства 
и термомеханическая численная модель генерации сверхлитостатического давления
О.П. Полянский, И.И. Лиханов, А.В. Бабичев, П.С. Козлов, С.В. Зиновьев, В.Г. Свердлова
19
Растворимость хлора в силикатных расплавах: новые эксперименты 
и термодинамическая модель смешения
Л.Я. Аранович, М.А. Голунова, Дж. А.Д. Коннолли, М.В. Иванов
46
Метаультрамафиты максютовского комплекса, Южный Урал: высокобарный 
Si-Al метасоматоз и карбонатизация на границе кора–мантия в зоне субдукции
А.Л. Перчук, Н.Г. Зиновьева, А.В. Сапегина, П.М. Вализер, В.М. Козловский, 
В.М. Григорьева, С.Т. Подгорнова
59
Метасоматизм в докембрийской коре Сибирского кратона: результат изучения 
ксенолитов гранат(±ортопироксен)-биотит-полевошпатовых пород из кимберлитовых 
трубок Юбилейная и Сытыканская, Якутия
Н.Е. Селютина, О.Г. Сафонов, В.О. Япаскурт, Д.А. Варламов, И.С. Шарыгин, 
К.М. Константинов, В.М. Козловский
91
Геохимическая термометрия рудоносных габброноритов из апофиза йоко-довыренского 
массива: состав, количество оливина и условия сульфидного насыщения исходной магмы
И.В. Пшеницын, А.А. Арискин, С.Н. Соболев
119

ПЕТРОЛОГИЯ,  2024, том 32, № 1,  с.  3–4
“В КИЛЬВАТЕРЕ БОЛЬШОГО КОРАБЛЯ”: 
ПЛАВАНЬЕ ПРОДОЛЖАЕТСЯ!
Специальный выпуск журнала “Петрология” 
посвящен памяти выдающегося отечественного 
петролога – Леонида Львовича Перчука, 90-летний юбилей которого отмечался 20 ноября 2023 г. 
в  рамках уже ставшей традиционной научной 
конференции “В кильватере большого корабля”. 
Помимо огромного собственно научного вклада, 
важной заслугой Леонида Львовича было создание 
научной школы, работы представителей которой 
предлагаются вниманию читателей.
В статье Л.А. Хакимовой и Ю.Ю. Подладчикова предложена гидро-хемо-механическая модель, 
позволяющая в рамках унифицированного подхода проводить расчеты, связанные с фильтрацией 
многокомпонентного флюида в деформируемых 
химически активных вмещающих породах с учетом 
изменения плотностей сосуществующих фаз и их 
химического состава.
В работе О.П. Полянского с коллегами развиваются представления о вкладе тектонических процессов в давление при метаморфизме и показана 
возможность превышения давления над литостатическим в локальном масштабе в породах, попавших в условия сдвиговых деформаций.
Л.Я. Аранович с соавторами получили новые 
экспериментальные данные по растворимости хлоридов в модельных базальтах и представили новую 
термодинамическую модель, описывающую энергию смешения солевых частиц в базальтовом и гранитном расплавах.
А.Л. Перчук с коллегами детально описали минеральные ассоциации и структуру Ti-клиногумитовых ультрамафитов максютовского комплекса 
и на основании термодинамического моделирования показали роль высокобарного Si-Al метасоматоза и карбонатизации на границе кора–мантия 
в образовании этих пород.
Н.Е. Селютина с соавторами реконструировали Р-Т и флюидные условия метаморфизма гранат-биотит-полевошпатовых и ортопироксен-гранат-биотит-полевошпатовых ксенолитов из кимберлитовых трубок Юбилейная и Сытыканская. 
Им впервые удалось установить декомпресионный 
тренд метаморфизма.
В статье И.В. Пшеницына с коллегами приводятся оценки параметров исходной магмы рудоносного апофиза Йоко-Довыренского массива, 
основанные на результатах термодинамического 
3

ПЕТРОЛОГИЯ,  2024, том 32, № 1,  с.  4–4
гранулитовой фации рассматривается роль растворенной во флюидах серы в реакциях дегидратации, 
в состоянии окисленности пород и в подвижности 
редких элементов.
Этот краткий обзор показывает, что “выпускники” и последователи научной школы Л.Л. Перчука 
продолжают развивать главные направления исследований, проводившихся любимым учителем: петрология и эволюция условий метаморфизма, экспериментальное и теоретическое моделирование 
метаморфических реакций и процессов плавления, 
термодинамические свойства породообразующих 
растворов и расплавов, численное моделирование 
геодинамических процессов и общетеоретические 
проблемы физико-химической петрологии.
Плавание в  “кильватере большого корабля” 
продолжается!
Л.Я. Аранович,
приглашенный редактор
моделирования равновесной кристаллизации расплавов по развиваемому А.А. Арискиным методу геохимической термометрии. Авторам удалось 
установить начальную температуру исходной магмы, ее состав и порядок кристаллизации, а также 
близость к насыщению магмы сульфидной серой.
Еще две статьи приглашенных авторов публикуются только в английской версии журнала.
Т.В. Геря (T.V. Gerya) предложил новое определение крупномасштабной долговременной прочности литосферы как меры ее механического сопротивления необратимым деформациям. В соответствии с этим новым определением прочность 
является отношением интегральной (по объему 
литосферы и  времени) диссипации механической энергии к интегральному необратимому вязко-пластическому напряжению.
В  обзорной статье Д. Харлова (D. Harlov) 
на примерах хорошо изученных гранитоидов 
4

ПЕТРОЛОГИЯ,  2024, том 32, № 1,  с.  5–18
УДК 550.41
МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОГО 
ФЛЮИДА … В  ДЕФОРМИРУЕМЫХ ВМЕЩАЮЩИХ ПОРОДАХ 
С УЧЕТОМ ХИМИЧЕСКОГО ВЗАИМОДЕЙСТВИЯ
© 2024 г.   Л.А. Хакимова, a, b, c, *, Ю.Ю. Подладчиковa, c
aUniversity of Lausanne, Lausanne, Switzerland
bСколковский институт науки и технологий, Центр добычи углеводородов, Москва, Россия
cМосковский государственный университет имени М.В. Ломоносова, Механико-математический факультет, 
Москва, Россия
*e-mail: liudmila.khakimova@unil.ch
Поступила в редакцию 13.06.2023 г.
После доработки 02.08.2023 г.
Принята к публикации 16.09.2023 г.
Мы предлагаем сопряженную гидро-хемо-механическую модель и ее численную 1D-реализацию, 
позволяющую в рамках унифицированного подхода проводить расчеты, связанные с фильтрацией многокомпонентного флюида в деформируемых химически активных вмещающих породах 
с учетом изменения плотностей сосуществующих фаз и их химического состава.
Ключевые слова: гидро-хемо-механическая сопряженная модель, реактивный транспорт, численное 
моделирование, термодинамическая согласованность
DOI: 10.31857/S0869590324010021
ВВЕДЕНИЕ
CH4 и др. – вызывают химические и фазовые превращения как в поровых флюидах, так и в породах, через которые фильтруются поровые жидкости и с которыми они химически взаимодействуют.
Более 70-ти лет назад, на основании детального 
парагенетического анализа высоко метаморфизованных пород и гранитоидов Восточной Сибири, 
Д.С. Коржинский сформулировал принцип подвижности щелочей при метаморфизме и гранитизации (Коржинский, 1946, 1961). В соответствии 
с этим принципом химические потенциалы K2О 
и Na2О могут являться такими же интенсивными 
факторами глубинного минералообразования, как 
и химические потенциалы летучих компонентов.
Идеальная термодинамическая модель смешения предполагает простейшее и  универсальное 
равенство химической активности любого компонента его мольной доли в растворе. Однако термодинамические свойства смешения компонентов 
значительно отличаются от идеальных как в смесях 
Н2О-неполярный газ (СО2, N2, CH4, H2 и др.), так 
и в смесях Н2О-сильный электролит, таких как рассолы H2O-NaCl, H2O-KCl и H2O-(Na, K)Cl. Активность воды в смесях вода-неполярный газ больше, 
чем ее мольная доля, причем положительное отклонение от идеальности возрастает с ростом давления (например, Jacobs, Kerrick, 1981; Aranovich, 
Моделирование фильтрационных течений 
многофазных многокомпонентных смесей с учетом химических и фазовых превращений, а также термомеханических эффектов взаимодействия 
порового флюида с деформируемыми и химически активными вмещающими породами является 
ключевой и основной задачей при количественном и предсказательном исследовании ряда важнейших техногенных и природных геодинамических процессов, а также при решении прикладных 
задач, связанных с этими процессами. Успешное 
решение этих задач требует развития сопряженных математических моделей физических процессов и их эффективной программной реализации, 
учитывающей особенности современных вычислительных архитектур. Одним из ключевых компонентов такого подхода является самосогласованная термодинамическая модель, позволяющая 
описывать как фазовые равновесия многокомпонентного флюида и вмещающих пород и плотности сосуществующих фаз, так и теплофизические 
и реологические свойства фаз.
Общепринято, что изменения температуры, 
давления и концентрации летучих компонентов 
в поровых флюидах или расплавах – Н2О, СО2, 
5

ХАКИМОВА, ПОДЛАДЧИКОВ
Newton, 1999; Аранович, 2013), вплоть до появления второй критической точки.
В противоположность смесям вода-неполярный газ, в соответствии с экспериментальными 
данными (Aranovich, Newton, 1996, 1997), активность воды a(H2O) в рассолах H2O-NaCl, H2O-KCl 
и H2O-(Na, K)Cl сильно зависит от давления (при 
постоянной температуре и концентрации солей), 
резко уменьшаясь от значений, близких к идеальному молекулярному раствору, при сравнительно 
низком давлении ~2 кбар до значений, близких 
к идеальному ионному раствору, aH2O (XH2O)2, при 
10 кбар. Тенденция уменьшения активности воды 
в рассолах с ростом давления сохраняется по крайней мере до 45 кбар (Tropper, Manning, 2004).
Как следствие, в  присутствии рассолов 
кварц-полевошпатовые породы плавятся при гораздо более высокой температуре, чем в присутствии флюида Н2О-СО2 с той же концентрацией 
Н2О. Реакции дегидратации, напротив, протекают в рассолах при более низкой температуре, чем 
в смесях Н2О-неполярный газ. В качестве примера на изобарической (Р = 7 кбар) Т–XH2O диаграмме Л.Я. Аранович (Aranovich, 2017) сопоставил 
кривые плавления простого гранита и реакции 
дегидратации флогопита с кварцем в присутствии 
рассола H2O-(Na, K)Cl и смеси Н2О-СО2 по экспериментальным данным (Aranovich et al., 2013) 
и (Aranovich, Newton, 1998) соответственно. При 
одной и той же XH2O = 0.5 плавление в присутствии 
рассола происходит при температуре (~800°C) 
примерно на 80°C выше, чем в Н2О-СО2 (~720°C), 
а реакция дегидратации – почти на 100°C ниже. 
Этот пример показывает, что в присутствии рассолов при фиксированном давлении появляется 
окно значений температуры, внутри которого дегидратация минералов (и метасоматические преобразования, связанные с привносом/выносом 
щелочей и Са) может протекать без плавления 
пород (Аранович и др., 1987), причем, как показывают экспериментальные данные (Aranovich et 
al., 2013), это окно расширяется с увеличением 
давления.
Важным фактором взаимодействия флюид–порода являются также значительные вариации растворимости породообразующих минералов в зависимости от химической природы флюида, что 
в свою очередь может существенно влиять на пористость/проницаемость пород (Aranovich, 2017; 
Aranovich et al., 2020).
Таким образом, огромное разнообразие эффектов взаимодействия сложных флюидов с кристаллическими породами необходимо учитывать 
при моделировании потоков поровых флюидов, 
поскольку именно ими во многом определяются 
фазовый состав и, соответственно, плотность пород, и, как следствие, их пористость, а значит, и их 
проницаемость (Malvoisin et al., 2015; Plumper et al., 
2017; Beinlich et al., 2020).
Данные наблюдений о сосуществующих минералах и исследование их локальных равновесий позволили оценить Р-Т эволюцию глубинных метаморфических комплексов (Перчук, 2006; Perchuk, 
2011, а также ссылки в них). Всеобщее признание 
надежности этой количественной информации 
(исторический обзор в Green, 2005), естественно, 
привело к возможности и необходимости моделирования геодинамических процессов (Brown, 2014, 
а также ссылки в ней). Первые геодинамические 
модели были основаны на термо-механических 
моделях несжимаемой вязкой жидкости, учитывающих фазовые превращения только для оценок 
плавучести в так называемом приближении Буссинеска. К удивлению и сожалению, моделирование 
Р-Т эволюции глубинных метаморфических комплексов столкнулось с серьезными проблемами 
по всем направлениям. Полученные оценки изменения температуры в природных процессах привели к пониманию слишком быстрого остывания 
мантийных плюмов (Weinberg, Podladchikov, 1994 
и ссылки в ней) и потерянном тепле (missing heat) 
в орогенезе (Hartz, Podladchikov, 2008 и ссылки 
в ней) и к необходимости объяснения свыше тысячеградусного роста температуры в коровых псевдотахилитах. С первых установленных оценок давления появились подозрения о динамо-метаморфизме неопределенной природы (Petrini, Podladchikov, 
2000; Moulas et al., 2013, а также ссылки в них). 
Объяснение температур, превышающих консервативные оценки на основе простейшего приближения Буссинеска, потребовало введения механизмов диссипативного разогрева при необратимой деформации метаморфизуемых пород (Hartz, 
Podladchikov, 2008; Perchuk et al., 1992; Перчук, 
Подладчиков, 1993; Perchuk, Gerya, 2011). Отклонение давления от простейшего литостатического, 
равного весу столба вышележащих пород, обычно 
связывают с тектоническими стрессами (Petrini, 
Podladchikov, 2000; Moulas et al., 2013; Schmalholz, 
Podladchikov 2014, а также ссылки в них). Однако 
следы тектоники часто отсутствуют при записи пиковых значений давления. Это приводит к необходимости рассмотрения “автоклавных” процессов 
по аналогии с ростом давления пара в плотно закрытой пароварке (Добрецов, 1974).
Учет взаимного влияния химических превращений и термо-механических процессов в открытых 
системах при возможной фильтрации флюидов или 
расплавов требует развития термо-гидро-хемо-механических (Thermo-Hydro-Mechanical-Chemical, 
THMC) моделей. Основной причиной необходимости сопряжения термодинамических и  механических процессов являются сильные изменения плотности, важные для механики и надежно 
предсказываемые из термодинамики. В настоящее 
ПЕТРОЛОГИЯ
том 32
№ 1
2024

МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОГО ФЛЮИДА… 
7
соотношений для каждой из рассматриваемых 
фаз и компонентов, а также на основе принципов 
классической необратимой термодинамики. Предполагая гипотезу о локальном термодинамическом 
равновесии в  слабой формулировке (Yarushina, 
Podladchikov, 2015), можно получить выражения 
для определяющих соотношений, гарантирующих термодинамическую согласованность системы. Данная методология вывода термодинамически согласованных систем уравнений механики 
сплошной среды представлена в (Landau, Lifshitz, 
1987), а  также в  таких работах, как (Yarushina, 
Podladchikov, 2015; Rass et al., 2018, 2019; Duretz 
et al., 2019; Malvoisin et al., 2021), где приводится вывод систем уравнений самосогласованного 
сопряжения.
Закон сохранения массы всей системы или 
уравнение неразрывности имеет вид:
(
)
∂
∂ρ = −∂
∂
−
ρ φ +
ρ
t
x
v
v
v
(
)
,
t
i
i
f
i
s
f
i
s t
(1)
где rf, rs – массовые плотности флюида (или расплава) и твердой матрицы соответственно; φ – пористость; vi
f , vi
s – компоненты скоростей флюида и твердой матрицы в приближении сплошной 
среды соответственно; φρ
φ ρ
ρ
f
s
t
+
−
(
)
=
1
 – плотность всей системы.
Закон сохранения массы однофазного флюида:
∂
i
f
i
s
f
i
s
f
ρ φ
ρ φ
ρ φ
(
)
.
(2)
∂(
) = −∂
∂
−
+
(
)
t
x
v
v
v
f
i
Закон сохранения импульса флюида (или расплава) имеет форму закона Дарси:
f
φ
η
ρ
v
v
k
P
x
g
i
f
i
s
i
i
f
−
(
) = −
∂
∂
+

f




,
(3)
время нами ведутся разработки полностью сопряженных геодинамических моделей (Schmalholz 
et al., 2020; Malvoisin et al., 2021; Bessat et al., 2022).
В  рамках настоящей работы мы предлагаем 
один из промежуточных шагов в создании полностью сопряженных моделей – сопряженную гидро-хемо-механическую модель для расчета процессов, связанных с фильтрацией флюидов в деформируемых химически активных вмещающих 
породах, учитывающую изменения плотности 
минеральной матрицы, обусловленные фазовыми 
превращениями.
Предлагаемый подход включает в себя два принципиальных блока: расчет P-T-X многофазных равновесий в многокомпонентных системах с учетом 
химических превращений как части композиционного симулятора и его интеграция с гидро-хемо-механическим солвером. Для расчета P-T-X
равновесий мы используем метод прямой минимизации энергии Гиббса системы (Gordon, McBride, 
1994; Vrijmoed, Podladchikov, 2022; Исаева и др., 
2021; Khakimova et al., 2021). В этом случае константы фазового равновесия не используются в рабочих формулах алгоритма, но их значения можно 
легко вычислить, используя результаты алгоритма 
минимизации. Результаты работы алгоритма представляют собой табличные значения равновесных 
концентраций компонентов в равновесных фазах 
в зависимости от P-T условий и валового состава 
многофазной смеси, таким образом, они параметризуются в виде кусочно-линейных функций или 
гладких сплайнов для последующей интеграции 
с гидро-хемо-механическим солвером путем внутреннего сопряжения. Демонстрацию численной 
реализации мы провели для задачи фильтрации 
H2O-CO2 в деформируемой и реагирующей вмещающей породе в 1D-постановке на примере реакций 
дегидратации и карбонатизации. Исследуется формирование магнезита и талька из дегидратирующегося серпентинита при фильтрации водного флюида с низким содержанием CO2.
где P f – давление флюида, hf – динамическая вязкость флюида, gi – ускорение свободного падения, k = k(φ) – проницаемость, которая нелинейно зависит от пористости, согласно соотношению 
Кармана–Козени:
n
=

k
k
bg


0
φ
φ
,
(4)



МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА 
ЗАДАЧИ И ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ
Законы сохранения и определяющие соотношения
Мы рассматриваем процесс фильтрации однофазного сверхкритического флюида (или расплава) 
в реагирующей и деформирующей вязкоупругой 
среде. Рассматриваемая двухфазная система порода–флюид состоит из N компонентов, из которых 
Nim < N компонентов не могут переходить во флюид (или расплав) при интересующих P-T-X условиях, а N – Nim = Nm компонентов могут быть как 
в составе твердой фазы, так и в составе флюида.
Система уравнений, описывающая данный процесс, может быть получена на основе балансовых 
где k0 – референсная проницаемость, φbg – референсная пористость среды, n – показатель нелинейности. Проницаемость сильно нелинейных или 
трещиноватых пород можно воспроизвести, используя большой показатель степени и, в этом случае, данное соотношение можно интерпретировать 
как линейную аппроксимацию экспериментально 
измеренной зависимости проницаемости от пористости в логарифмических координатах (Yarushina 
et al., 2021).
Закон сохранения полного импульса обеих фаз, 
где предполагается, что инерциальные силы пренебрежимо малы, имеет вид:
ПЕТРОЛОГИЯ
том 32
№ 1
2024

ХАКИМОВА, ПОДЛАДЧИКОВ
∂
t
ij
ij
t
i
t
δ
τ
ρ
0,
(5)
∂
−
+
(
) −
=
x
P
g
j
где полное давление системы, P
P
P
f
s
t
φ
φ
+
−
(
) =
1
и компоненты девиатора напряжений, τ φ
τ
φ
τ
ij
f
ij
s
ij
t
+
−
(
) =
1
φ
τ
φ
τ
ij
f
ij
s
будет использовано для интегрирования по времени закона сохранения массы.
Закон сохранения массы для Nm – 1 мобильных 
компонентов системы записывается в следующем 
виде в предположении отсутствия диффузионных 
потоков:
t k
,
t C
x
v
v
C
v
C
,
,
,
(
)
(
)
∂
∂ρ
= −∂
∂
−
ρ φ
+
ρ
i
i
f
i
s
f
f k
s
t k
(9)
m
k
N
1...
1,
=
−
ij
t
+
−
(
) =
1
введены согласно формулировкам, представленным в (Lopatnikov, Cheng, 2002); dij – дельта 
Кронекера.
Определяющее соотношение, связывающее 
компоненты девиатора напряжений и градиенты 
компонентов скорости твердого скелета, удовлетворяющее вязкоупругой модели Максвелла, имеет 
форму (Beuchert, Podladchikov, 2010):
ρ φ
ρ
φ
ρ
f
f k
s
s k
t k
C
C
C
,
,
, ,
+
−
(
)
=
1
(10)
s
1
τ
τ
∂
∂
+
∂
∂
=
+
x v
x v
G
D
Dt
i
j
s
j
i
s
ij
t
ij
t
η
,
(6)
где G – модуль сдвига; hs – вязкость твердого скелета; D
Dt  – производная Яумана.
Определяющее соотношение для скорости деформации насыщенной поровязкоупругой среды, 
удовлетворяющее термодинамической согласованности, выведено в (Yarushina, Podladchikov, 2015) 
и имеет вид:
P
P
где C f,k, Cs,k  – массовая концентрация компонента k во флюидной и твердой фазах соответственно; 
Ct,k – полная массовая концентрация компонента k
в системе.
Далее мы рассматриваем уравнение, описывающее закон сохранения массы нереагирующей части твердой матрицы. Мы рассматриваем систему, 
в которой Nim компонентов не могут переходить во 
флюид при интересующих P-T-X условиях. В таком 
случае закон сохранения массы всех компонентов 
(k = 1…Nim), составляющих нереагирующую часть 
твердого скелета, записывается следующим образом:
∂
∂
= −
−

k
k
s
d
im
s
t
f
f
f
t
1
α
ηφ
, (7)

+
−
x v
K
d p
dt
d
p
dt



∂
1
(
)
,
1
s
s k
k
N
s
=
t
C
∂
ρ
−φ





=
f
,
  – субстан(11)
i
f
im
∑
где ∂
∂+
∂
∂
=
∂
∂+
∂
∂
=
t
v
x
d
dt
t
v
x
d
dt
i
s
i
i
1
,
(
)
∑
,
1
=
x
v
C
= −∂
∂
ρ
−φ






i
i
s
s
s k
k
N
где Cs,k  – массовая концентрация компонента k 
(k = 1…Nim) в твердой матрице.
Термодинамическая модель
циональная производная по отношению к системе отсчета, связанной с vs и v f соответственно, 
Kd  – объемный модуль упругости насыщенной 
пористой среды, a – коэффициент Био–Виллиса, hφ – эффективная объемная вязкость. В работе 
(Yarushina, Podladchikov, 2015) также показано, что 
для случая вязкой деформации определяющее соотношение (7) сводится к виду
f
t
.
∂
∂
=
−
x v
P
P
k
k
s
ηφ
Именно этот случай и будет рассмотрен в настоящей работе. Стоит отметить, что параметр hφ
отражает характер течения флюида в  пористой 
среде и может быть оценен как экспериментально 
(Zimmerman, 1991; Dong et al., 2010), так и теоретически с использованием методов теории эффективных сред (Yarushina, Podladchikov, 2015).
В то же время дивергенция поля скорости vs, 
кинематически связанная со скоростью изменения объемной деформации твердой матрицы, зависит от изменения объема системы V следующим 
образом:
∂
∂
=
=
(
)
x v
V
d V
dt
d
dt lnV
i
i
s
s
s
1
.
(8)
Это соотношение можно считать определяющим соотношением для объема системы V, и оно 
Для замыкания вышеописанных уравнений 
необходимо получить выражение для плотности 
результирующих фаз, а  также для равновесной 
термодинамической зависимости концентрации 
компонентов рассматриваемой системы в сосуществующих фазах.
В рамках предлагаемого подхода искомые функции представляют собой табулированные значения 
равновесных концентраций компонентов в сосуществующих фазах, а также плотностей сосуществующих фаз в зависимости от набора P-T-X условий, в которых система находится в разные моменты времени. Для расчета термодинамических 
равновесий используется метод прямой минимизации энергии Гиббса системы, реализованный 
в пакете ThermoLab (Vrijmoed, Podladchikov, 2022). 
Пакет ThermoLab разработан на языке MATLAB 
и уже включает в себя большой набор современных 
термодинамических баз данных для расчета свойств минералов, расплавов и флюидов, а также фазовых равновесий (Holland, Powell, 1998; Aranovich, 
ПЕТРОЛОГИЯ
том 32
№ 1
2024

МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОГО ФЛЮИДА… 
9
Newton, 1999; White et al., 2003; Padrόn-Navarta 
et al., 2013). На практике результаты серии расчетов 
алгоритма минимизации для разных интересующих диапазонов P-T-Х условий сохраняются в виде 
отдельного многомерного массива данных с целью 
последующей обработки и формирования табулированных таблиц данных, поступающих в транспортный солвер.
Для сведения изначально нелинейной задачи 
минимизации энергии Гиббса, нелинейно зависящей от концентраций в растворах с линейными 
ограничениями по валовому составу системы к задаче линейного програмирования, применяется 
дискретизация непрерывных растворов по композиционному пространству, или, другими словами, 
все смеси переменного состава (твердые растворы, 
жидкие растворы, расплавы) приближенно заменяются на большой набор фаз, имеющих фиксированный состав (Connolly, 2005; Коннолли, 2017). 
Алгоритм минимизации определяет, какие из этих 
дискретных фаз находятся в равновесии. Технически результатом расчета равновесия является вектор, который содержит общее количество молей 
каждой дискретной фазы из рассматриваемых в системе. После получения результатов алгоритма минимизации возможно вычисление необходимых наборов табулированных данных, характеризующих 
свойства результирующих фаз, согласно формулам, 
представленным в (Vrijmoed, Podladchikov, 2022).
Рис. 1. Цикл по времени 1D-численной имплементации гидро-хемо-механической модели, реализованной в MATLAB.
Численная имплементация
твердого скелета, т.е. в данном случае для C, записывается следующим образом (см. уравнение (11) 
при  Nim = 1):
∂
i
s
s
s C
ρ
φ
ρ
φ
1
1
,
,
. (14)
∂
−
(
)
(
) = −∂
∂
−
(
)
(
)
t
C
x
v
C
s
s C
i
Используя уравнение (8) для дивергенции скорости твердой фазы и субстанциональную произs
, мы можем привести уравнение (11) 
водную d
dt
к следующему:
d
s
s C
ρ
φ
1
0
−
(
)
(
) =
,
,
(15)
dt
VC
s
Стоит отметить, что уравнение (9) представляет собой закон сохранения полной массы компонентов системы, которое распадается на несколько 
уравнений в зависимости от количества мобильных Nm (могут присутствовать как в твердой, так 
и в жидкой фазах) и немобильных Nim (могут присутствовать только в твердой фазах) компонентов 
системы. Для пояснения приведем пример, где будем рассматривать систему, состоящую из Nc = 3 
компонентов, A, B и C, предполагая, что N m = 2 
компонентов, а именно A и B, могут переходить во 
флюид (или расплав), а Nim = 1 компонент, а именно C, при рассматриваемых P-T-Х условиях не может переходить во флюид (или расплав). В данном 
случае закон сохранения массы для N m = 1 мобильных компонентов системы (например, для компонента A) записывается в следующем виде (см. уравнения 9 и 10 при N m = 2):
которое является обыкновенным дифференциальным уравнением вдоль траекторий твердой фазы 
и может быть проинтегрировано по времени с учетом начальных условий:
∂
i
f
i
s
f
f A
i
s
t A
ρ
φρ
ρ
,
,
,
)
,(12)
∂
= −∂
∂
−
+
(
)
t C
x
v
v
C
v
C
t A
i
ρ
φ
ρ
φ
s
s C
s
s C
VC
V C
1
1
0
0
0
0
−
(
)
=
−
(
)
,
, ,
(16)
где ρ0
s , φ0, V0, C s C
φρ
φ ρ
ρ
f
f A
s
s A
t A
C
C
C
,
,
, .
+
−
(
)
=
1
(13)
Далее закон сохранения массы всех компонентов Nim , составляющих нереагирующую часть 
0
,  – плотность твердого скелета, 
пористость, объем и концентрация компонента C
в твердой матрице в начальный момент времени 
соответственно.
ПЕТРОЛОГИЯ
том 32
№ 1
2024