Modelling of Multicomponent Fluid Flow in Deforming and Reacting Porous Rock

Full Text

ВВЕДЕНИЕ

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

Общепринято, что изменения температуры, давления и концентрации летучих компонентов в поровых флюидах или расплавах – H₂O, CO₂, CH₄ и др. – вызывают химические и фазовые превращения как в поровых флюидах, так и в породах, через которые фильтруются поровые жидкости и с которыми они химически взаимодействуют.

Более 70-ти лет назад, на основании детального парагенетического анализа высоко метаморфизованных пород и гранитоидов Восточной Сибири, Д.С. Коржинский сформулировал принцип подвижности щелочей при метаморфизме и гранитизации (Коржинский, 1946, 1961). В соответствии с этим принципом химические потенциалы K₂O и Na₂О могут являться такими же интенсивными факторами глубинного минералообразования, как и химические потенциалы летучих компонентов.

Идеальная термодинамическая модель смешения предполагает простейшее и универсальное равенство химической активности любого компонента его мольной доли в растворе. Однако термодинамические свойства смешения компонентов значительно отличаются от идеальных как в смесях H₂O-неполярный газ (CO₂, N₂, CH₄, H₂ и др.), так и в смесях H₂O-сильный электролит, таких как рассолы H₂O-NaCl, H₂O-KCl и H₂O-(Na, K)Cl. Активность воды в смесях вода-неполярный газ больше, чем ее мольная доля, причем положительное отклонение от идеальности возрастает с ростом давления (например, Jacobs, Kerrick, 1981; Aranovich, Newton, 1999; Аранович, 2013), вплоть до появления второй критической точки.

В противоположность смесям вода-неполярный газ, в соответствии с экспериментальными данными (Aranovich, Newton, 1996, 1997), активность воды a(H₂O) в рассолах H₂O-NaCl, H₂O-KCl и H₂O-(Na, K) Cl сильно зависит от давления (при постоянной температуре и концентрации солей), резко уменьшаясь от значений, близких к идеальному молекулярному раствору, при сравнительно низком давлении ~2 кбар до значений, близких к идеальному ионному раствору, a_H₂O (X_H₂O)², при 10 кбар. Тенденция уменьшения активности воды в рассолах с ростом давления сохраняется по крайней мере до 45 кбар (Tropper, Manning, 2004).

Как следствие, в присутствии рассолов кварц-полевошпатовые породы плавятся при гораздо более высокой температуре, чем в присутствии флюида H₂O-CO₂ с той же концентрацией H₂O. Реакции дегидратации, напротив, протекают в рассолах при более низкой температуре, чем в смесях H₂O-неполярный газ. В качестве примера на изобарической (Р = 7 кбар) Т–X_H₂O диаграмме Л.Я. Аранович (Aranovich, 2017) сопоставил кривые плавления простого гранита и реакции дегидратации флогопита с кварцем в присутствии рассола H₂O-(Na, K)Cl и смеси H₂O-CO₂ по экспериментальным данным (Aranovich et al., 2013) и (Aranovich, Newton, 1998) соответственно. При одной и той же X_H₂O = 0.5 плавление в присутствии рассола происходит при температуре (~800°C) примерно на 80°C выше, чем в H₂O-CO₂(~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) моделей. Основной причиной необходимости сопряжения термодинамических и механических процессов являются сильные изменения плотности, важные для механики и надежно предсказываемые из термодинамики. В настоящее время нами ведутся разработки полностью сопряженных геодинамических моделей (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 условий и валового состава многофазной смеси, таким образом, они параметризуются в виде кусочно-линейных функций или гладких сплайнов для последующей интеграции с гидро-хемо-механическим солвером путем внутреннего сопряжения. Демонстрацию численной реализации мы провели для задачи фильтрации H₂O-CO₂ в деформируемой и реагирующей вмещающей породе в 1D-постановке на примере реакций дегидратации и карбонатизации. Исследуется формирование магнезита и талька из дегидратирующегося серпентинита при фильтрации водного флюида с низким содержанием CO₂.

МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ И ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ

Законы сохранения и определяющие соотношения

Мы рассматриваем процесс фильтрации однофазного сверхкритического флюида (или расплава) в реагирующей и деформирующей вязкоупругой среде. Рассматриваемая двухфазная система порода–флюид состоит из N компонентов, из которых N^im < N компонентов не могут переходить во флюид (или расплав) при интересующих P-T-X условиях, а N – N^im = N^m компонентов могут быть как в составе твердой фазы, так и в составе флюида.

Система уравнений, описывающая данный процесс, может быть получена на основе балансовых соотношений для каждой из рассматриваемых фаз и компонентов, а также на основе принципов классической необратимой термодинамики. Предполагая гипотезу о локальном термодинамическом равновесии в слабой формулировке (Yarushina, Podladchikov, 2015), можно получить выражения для определяющих соотношений, гарантирующих термодинамическую согласованность системы. Данная методология вывода термодинамически согласованных систем уравнений механики сплошной среды представлена в (Landau, Lifshitz, 1987), а также в таких работах, как (Yarushina, Podladchikov, 2015; Rass et al., 2018, 2019; Duretz et al., 2019; Malvoisin et al., 2021), где приводится вывод систем уравнений самосогласованного сопряжения.

Закон сохранения массы всей системы или уравнение неразрывности имеет вид:

$\frac{\partial }{\partial t}{\rho }^t=-\frac{\partial }{\partial x_i}\left(\left(v_i^f\mathit{-}{v}_i^s\right){\mathrm{\rho }}^f\mathrm{\varphi }+v_{\mathit{i}}^{\mathit{s}}{\mathrm{\rho }}^f\mathrm{\varphi }\right)$                                                                  (1)

где ${\rho }^f, {\rho }^s$ – массовые плотности флюида (или расплава) и твердой матрицы соответственно; $\mathrm{\varphi }$ – пористость; $v_i^f$, $v_i^s$ – компоненты скоростей флюида и твердой матрицы в приближении сплошной среды соответственно; ${\mathrm{\varphi \rho }}^f+\left(1-\mathrm{\varphi }\right){\mathrm{\rho }}^{\mathrm{s}}={\mathrm{\rho }}^{\mathit{t}}$ – плотность всей системы.

Закон сохранения массы однофазного флюида:

$\frac{\partial }{\partial t}\left({\rho }^f\mathrm{\varphi }\right)=-\frac{\partial }{\partial x_i}\left(\left(v_{\mathit{i}}^{\mathit{f}}\mathit{-}{v}_{\mathit{i}}^{\mathit{s}}\right){\mathrm{\rho }}^f\mathrm{\varphi }+v_{\mathit{i}}^{\mathit{s}}{\mathrm{\rho }}^f\mathrm{\varphi }\right)$                                                              (2)

Закон сохранения импульса флюида (или расплава) имеет форму закона Дарси:

$\mathrm{\varphi }\left({\mathrm{v}}_i^f-{\mathrm{v}}_{\mathit{i}}^{\mathit{s}}\right)=-\frac{k}{{\mathrm{\eta }}_f}\left(\frac{\partial P^f}{\partial {\mathrm{x}}_{\mathrm{i}}}+{\mathrm{g}}_i{\mathrm{\rho }}^f\right)$                                                                       (3)

где P^f – давление флюида, η_f – динамическая вязкость флюида, g_i – ускорение свободного падения, k = k(φ) – проницаемость, которая нелинейно зависит от пористости, согласно соотношению Кармана–Козени:

$k=k_0{\left(\frac{\mathrm{\varphi }}{{\mathrm{\varphi }}_{bg}}\right)}^n,$                                                                                       (4)

где k₀ – референсная проницаемость, ${\mathrm{\varphi }}_{\mathit{ }bg }$– референсная пористость среды, n – показатель нелинейности. Проницаемость сильно нелинейных или трещиноватых пород можно воспроизвести, используя большой показатель степени и, в этом случае, данное соотношение можно интерпретировать как линейную аппроксимацию экспериментально измеренной зависимости проницаемости от пористости в логарифмических координатах (Yarushina et al., 2021).

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

$\frac{\partial }{\partial x_i}\left(-P^t{\delta }_{ij}+{\tau }_{ij}^t\right)-g_i{\rho }^t=0,$                                                                    (5)

где полное давление системы, $P^f\mathrm{\varphi }+P^{\mathrm{s}}\left(1-\mathrm{\varphi }\right)\mathit{=}{P}^{\mathit{t}}$ и компоненты девиатора напряжений, ${\tau }_{ij}^f\mathrm{\varphi }+{\mathrm{\tau }}_{ij}^{\mathrm{s}}\left(1-\mathrm{\varphi }\right)={\mathrm{\tau }}_{\mathit{i}\mathit{j}}^{\mathit{t}}$ введены согласно формулировкам, представленным в (Lopatnikov, Cheng, 2002); δ_ij – дельта Кронекера.

Определяющее соотношение, связывающее компоненты девиатора напряжений и градиенты компонентов скорости твердого скелета, удовлетворяющее вязкоупругой модели Максвелла, имеет форму (Beuchert, Podladchikov, 2010):

$\frac{{\partial }^{}}{\partial x_i}{v}_j^s+\frac{{\partial }^{}}{\partial x_j}{v}_i^s=\frac{1}{G}\frac{D{\tau }_{ij}^t}{Dt}+\frac{{\tau }_{ij}^t}{{\mathrm{\eta }}^s}$                                                               (6)

где G – модуль сдвига; η^s – вязкость твердого скелета; $\frac{D}{Dt}$ – производная Яумана.

Определяющее соотношение для скорости деформации насыщенной поровязкоупругой среды, удовлетворяющее термодинамической согласованности, выведено в (Yarushina, Podladchikov, 2015) и имеет вид:

$\frac{\partial }{\partial x_k}{v}_k^s=-\frac{1}{K_d}\left(\frac{d^s{p}^t}{dt}-\alpha \frac{d^f{p}^f}{dt}\right)+\frac{P^f{P}^t}{{\mathrm{\eta }}_{\mathrm{\varphi }}},$ (7)

где $\frac{\partial }{\partial t}+v_k^s\frac{\partial }{\partial x_i}=\frac{d^s}{dt}, \frac{\partial }{\partial t}+v_i^f \frac{\partial }{\partial x_i}=\frac{d^f}{dt}$ – субстанциональная производная по отношению к системе отсчета, связанной с v^s и v^f соответственно, K_d – объемный модуль упругости насыщенной пористой среды, α – коэффициент Био–Виллиса, ${\mathrm{\eta }}_{\mathrm{\varphi }}$ – эффективная объемная вязкость. В работе (Yarushina, Podladchikov, 2015) также показано, что для случая вязкой деформации определяющее соотношение (7) сводится к виду

$\frac{\partial }{\partial x_k}{v}_k^s=\frac{P^f-P^t}{{\mathrm{\eta }}_{\mathrm{\varphi }}}$.

Именно этот случай и будет рассмотрен в настоящей работе. Стоит отметить, что параметр ${\mathrm{\eta }}_{\mathrm{\varphi }}$ отражает характер течения флюида в пористой среде и может быть оценен как экспериментально (Zimmerman, 1991; Dong et al., 2010), так и теоретически с использованием методов теории эффективных сред (Yarushina, Podladchikov, 2015).

В то же время дивергенция поля скорости v^s, кинематически связанная со скоростью изменения объемной деформации твердой матрицы, зависит от изменения объема системы V следующим образом:

$\frac{\partial }{\partial x_i}{v}_i^s=\frac{1}{V}\frac{d^s}{dt}=\frac{d^s}{dt}\left(lnV\right)$                                                               (8)

Это соотношение можно считать определяющим соотношением для объема системы V, и оно будет использовано для интегрирования по времени закона сохранения массы.

Закон сохранения массы для N^m – 1 мобильных компонентов системы записывается в следующем виде в предположении отсутствия диффузионных потоков:

$\begin{array}{c}\frac{\partial }{\partial t}\rho C^{t,k}=-\frac{\partial }{\partial x_i}\left(\left(v_i^f-v_i^s\right){\rho }^f\mathrm{\varphi }{C}^{\mathit{f}\mathit{,}\mathit{k}}+{\mathrm{v}}_{\mathrm{s}}\mathrm{\rho }{C}^{\mathit{t}\mathit{,}\mathit{k}}\right),\\ k=1...N^m-1,\end{array}$                                                  (9)

${\rho }^f\mathrm{\varphi }{C}^{f,k}+{\rho }^s\left(1-\mathrm{\varphi }\right)C^{s,k}=\rho C^{t,k},$                                                                   (10)

где C^f,k, C^s,k – массовая концентрация компонента k во флюидной и твердой фазах соответственно; C^t,k – полная массовая концентрация компонента k в системе.

Далее мы рассматриваем уравнение, описывающее закон сохранения массы нереагирующей части твердой матрицы. Мы рассматриваем систему, в которой N^im компонентов не могут переходить во флюид при интересующих P-T-X условиях. В таком случае закон сохранения массы всех компонентов (k = 1…N^im), составляющих нереагирующую часть твердого скелета, записывается следующим образом:

$\frac{\partial }{\partial t}\left({\rho }^s\left(1-\mathrm{\varphi }\right)\sum _{k=1}^{N^{im}}{C}^{s,k}\right)=\frac{\partial }{\partial x_i}\left(v_i^s{\rho }^s\left(1-\mathrm{\varphi }\right)\sum _{k=1}^{N^{im}}{C}^{s,k}\right),$                                       (11)

где C^s,k – массовая концентрация компонента k (k = 1…N^im) в твердой матрице.

Термодинамическая модель

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

В рамках предлагаемого подхода искомые функции представляют собой табулированные значения равновесных концентраций компонентов в сосуществующих фазах, а также плотностей сосуществующих фаз в зависимости от набора P-T-X условий, в которых система находится в разные моменты времени. Для расчета термодинамических равновесий используется метод прямой минимизации энергии Гиббса системы, реализованный в пакете ThermoLab (Vrijmoed, Podladchikov, 2022). Пакет ThermoLab разработан на языке MATLAB и уже включает в себя большой набор современных термодинамических баз данных для расчета свойств минералов, расплавов и флюидов, а также фазовых равновесий (Holland, Powell, 1998; Aranovich, Newton, 1999; White et al., 2003; Padrόn-Navarta et al., 2013). На практике результаты серии расчетов алгоритма минимизации для разных интересующих диапазонов P-T-Х условий сохраняются в виде отдельного многомерного массива данных с целью последующей обработки и формирования табулированных таблиц данных, поступающих в транспортный солвер.

Для сведения изначально нелинейной задачи минимизации энергии Гиббса, нелинейно зависящей от концентраций в растворах с линейными ограничениями по валовому составу системы к задаче линейного програмирования, применяется дискретизация непрерывных растворов по композиционному пространству, или, другими словами, все смеси переменного состава (твердые растворы, жидкие растворы, расплавы) приближенно заменяются на большой набор фаз, имеющих фиксированный состав (Connolly, 2005; Коннолли, 2017). Алгоритм минимизации определяет, какие из этих дискретных фаз находятся в равновесии. Технически результатом расчета равновесия является вектор, который содержит общее количество молей каждой дискретной фазы из рассматриваемых в системе. После получения результатов алгоритма минимизации возможно вычисление необходимых наборов табулированных данных, характеризующих свойства результирующих фаз, согласно формулам, представленным в (Vrijmoed, Podladchikov, 2022).

Численная имплементация

Стоит отметить, что уравнение (9) представляет собой закон сохранения полной массы компонентов системы, которое распадается на несколько уравнений в зависимости от количества мобильных N^m (могут присутствовать как в твердой, так и в жидкой фазах) и немобильных N^im (могут присутствовать только в твердой фазах) компонентов системы.

 

Рис. 1. Цикл по времени 1D-численной имплементации гидро-хемо-механической модели, реализованной в MATLAB.

 

Для пояснения приведем пример, где будем рассматривать систему, состоящую из N^c = 3 компонентов, A, B и C, предполагая, что N^m = 2 компонентов, а именно A и B, могут переходить во флюид (или расплав), а N^im = 1 компонент, а именно C, при рассматриваемых P-T-Х условиях не может переходить во флюид (или расплав). В данном случае закон сохранения массы для N^m = 1 мобильных компонентов системы (например, для компонента A) записывается в следующем виде (см. уравнения 9 и 10 при N^m = 2):

$\frac{\partial }{\partial t}\rho C^{t,A}=-\frac{\partial }{\partial x_i}\left(\left(v_i^f-v_i^s\right){\mathrm{\varphi \rho }}^f{C}^{f,A}+v_i^s{\mathrm{\rho }}^f{C}^{f,A}\right)$                                             (12)

${\mathrm{\varphi \rho }}^f{C}^{f,A}+\left(1-\mathrm{\varphi }\right){\mathrm{\rho }}^s{C}^{s,A}=\mathrm{\rho }{C}^{\mathit{t}\mathit{,}\mathit{A}}\mathit{.}$                                                                    (13)

Далее закон сохранения массы всех компонентов N^im , составляющих нереагирующую часть твердого скелета, т. е. в данном случае для C, записывается следующим образом (см. уравнение (11) при N^im = 1):

$\frac{\partial }{\partial t}\left({\rho }^s\left(1-\varphi \right)C^{s,C}\right)=-\frac{\partial }{\partial x_i}\left(v_i^s{\rho }^s\left(1-\mathrm{\varphi }\right)C^{s,C}\right).$                                               (14)

Используя уравнение (8) для дивергенции скорости твердой фазы и субстанциональную производную $\frac{d^s}{dt}$, мы можем привести уравнение (11) к следующему:

$\frac{d^s}{dt}\left({\rho }^s\left(1-\mathrm{\varphi }\right)V{C}^{s,C}\right)=0$                                                                                      (15)

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

${\rho }^s\left(1-\mathrm{\varphi }\right)V{C}^{s,C}={\rho }_0^s\left(1- {\mathrm{\varphi }}_0\right)V_0{C}_0^{s,C}$                                                                    (16)

где  ${\rho }_0^s, {\mathrm{\varphi }}_0, V_0, C_0^{s,C}$  – плотность твердого скелета, пористость, объем и концентрация компонента C в твердой матрице в начальный момент времени соответственно.

Для решения данной системы уравнений мы используем метод конечных разностей на разнесенной сетке и ускоренный метод установления (Rass et al., 2022) для решения статических подзадач (подпроцессов), который технически заключается в рассмотрении итеративных шагов по времени наряду с физическими. На рис. 1 представлен типичный цикл по времени, который соответствует 1D-численной имплементации гидро-хемо-механической модели, реализованной в MATLAB.

РЕЗУЛЬТАТЫ

Пример расчета задачи фильтрации глубинных химически активных флюидов

Для демонстрации мы остановились на рассмотрении случая образования талька и магнезита из серпентинита при 300°C и 0.3 ГПа (Beinlich et al., 2020). Для этого была проведена 1D-численная реализация предложенной гидро-хемо-механической модели в MATLAB. Изначально всю расчетную область занимает серпентинит, находящийся в равновесии с поровым H₂O-CO₂ флюидом. Минеральный состав представлен на рис. 2в. Начальная пористость в модели составляет 0.02 и имеет аномалию с повышенным значением на нижней границе расчетной области (зеленая линия на рис. 2а и 2б. Также на нижней границе расчетной области происходит закачка флюида с некоторой концентрацией CO₂ (рис. 2г), не находящегося в равновесии с исходной породой. Флюид вступает в реакцию с исходной породой с образованием новых минеральных комплексов, постепенно трансформирующих породу.

 

Рис. 2. Начальные условия для численной 1D-гидро-хемо-механической модели: на нижней границе происходит закачка H₂O-CO₂.

(а) Начальное распределение пористости. (б) Начальное распределение концентрации CO₂ в твердой матрице. (в) Минеральный состав породы. Показана объемная доля соответствующего минерала в каждой точке расчетной области. Левая граница данного графика соответствует нижней границе графиков (а) и (б), т. е. на данном графике изменение минерального состава будет происходить слева направо. (г) Красной линией показано равновесное соотношение между концентрацией CO₂ во флюиде и в твердой матрице, ромбами показаны комбинации C^s,CO₂ и C^f,CO₂, которые формируются в ходе расчета. Соотношение C^s,CO₂ и C^f,CO₂, соответствующее закачиваемому флюиду, отмечено стрелкой.

 

В данном случае система состоит из четырех основных компонентов MgO, SiO₂, H₂O и CO₂, из которых два компонента считаются мобильными, а именно CO₂ и H₂O, и могут находиться как во флюиде, так и быть структурно связанными в минеральной матрице. Остальные компоненты в рассматриваемой системе не могут переходить во флюид и формируют нереагирующую (немобильную, с точки зрения транспортного кода) часть твердой матрицы. Таким образом, уравнения (12)–(16) в конкретном случае могут быть преобразованы к виду:

$\frac{\partial }{\partial t}\rho C^{t,C{O}_2}=\frac{\partial }{\partial x}\left(\left(v^f-v^s\right)\mathrm{\varphi }{\rho }^f{C}^{f,C{O}_2}+v^s\rho C^{t,C{O}_2}\right),$                                              (17)

${\mathrm{\varphi \rho }}^f{C}^{f,C{O}_2}+\left(1-\mathrm{\varphi }\right){\mathrm{\rho }}^s{C}^{s,C{O}_2}=\mathrm{\rho }{C}^{t,C{O}_2}$                                                                 (18)

$\mathrm{\varphi }=1-\frac{{\rho }_0^1 \left(1-{\mathrm{\varphi }}_0\right) V_0 C_0^{s,\mathrm{MgO}}}{{\mathrm{\rho }}^s V C^{s,\mathrm{MgO}}}$                                                                             (19)

где C^f,CO₂, C^s,CO₂, C^t,CO₂ – массовая концентрация CO₂ во флюиде, в твердой фазе и во всей системе соответственно; C^s,MgO – массовая концентрация MgO в твердом скелете. Что касается гидро-механической части вышеописанных уравнений, то она остается неизменной.

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

$\left\{\begin{array}{l}{C}^{f,C{O}_2}=C^{f,C{O}_2}\left(C^{f,C{O}_2}, T, P^f\right)\\ C^{s,\mathrm{MgO}}=C^{s,\mathrm{MgO}}\left(C^{s,\mathrm{MgO}} T, P^f\right)\\ {\rho }^f={\rho }^f\left(C^{s,\mathrm{MgO}} T, P^f\right)\\ {\rho }^s={\rho }^s\left(C^{s,\mathrm{MgO}} T, P^f\right)\end{array}\right.$                                                                     (20)

Расчет данных функций происходит в результате интерполяции предварительно рассчитанного и табулированного набора термодинамических данных на каждом шаге по времени транспортного гидро-хемо-механического солвера. Расчет набора термодинамических данных вынесен в препроцессинг, осуществляется с помощью linprog минимизации по циклу для каждой точки из интересующего диапазона условий P^f – T – C^f,CO₂ (рис. 3).

 

Рис. 3. Визуализация табулированных данных для системы MgO, SiO₂, H₂O и CO₂ при 300°C и 0.3 ГПa, которые используются для интерполяции в транспортном гидро-хемо-механическом коде. Результат серии расчетов c использованием linprog минимизации. (а) –минеральный состав, (б) – плотность твердой матрицы и флюида, (с) – массовая доля CO₂ во флюиде, (г) – массовая концентрация MgO в твердой матрице. Все графики построены в зависимости от C^s,CO₂.

 

На рис. 4 и 5 приведены результаты расчетов при закачке флюида с низким содержанием растворенного CO₂ (C^f,CO₂ = 0.44 мас. %). На рис. 4 представлены распределения пористости, эффективного давления, потока Дарси и концентрации CO₂ в твердой матрице в разные моменты времени. Вначале происходит распространение аномалии пористости по типу солитона и формирование стабильной зоны с более высокой пористостью и проницаемостью. На рис. 4г видно, что в данном случае не происходит формирования резкого фронта реакции. Данный факт также виден на рис. 5, где представлено распределение минерального состава в системе в начальный и конечный моменты времени. В данном случае антигорит не исчезает из системы, происходит постепенное формирование фронта реакции из серпентинита в магнезит-тальк-серпентинит. На рис. 5в красной линией показано равновесное соотношение между концентрацией CO₂ во флюиде и в твердой матрице. Ромбами показаны комбинации C^s,CO₂ и C^f,CO₂, которые формируются в ходе расчета.

 

Рис. 4. Результаты расчетов гидро-хемо-механической численной модели при закачке флюида с низким содержанием растворенного CO₂ (Cf,CO₂ = 0.44 мас. %).

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

 

Рис. 5. Результаты расчетов гидро-хемо-механической численной модели при закачке флюида с низким содержанием растворенного CO₂ (C^f,CO₂ = 0.44 мас. %). Количественный минеральный состав породы на (а) момент начала расчетов, т. е. закачки флюида и (б) финальный момент расчетов, соответствующий образованию фронта реакции. (в) Равновесное соотношение между концентрацией CO₂ во флюиде и твердой матрицей (красная линия), ромбами показаны комбинации C^s,CO₂ и C^f,CO₂, которые формируются в ходе расчета. В данном случае не происходит формирование резкого фронта реакции, антигорит не исчезает из системы, происходит постепенное формирование фронта реакции из серпентинита в магнезит-тальк-серпентинит.

 

Концептуально другой режим реализуется при закачке флюида с чуть более высоким содержанием CO₂ (C^f,CO₂ = 1.38 мас. %). Вначале, аналогично предыдущему случаю, происходит волнообразное распространение аномалии пористости по механизму вязкой (де)компакции, приводящее к увеличению скорости Дарси потока (рис. 6), поле чего начинается формирование и распространение резкого фронта реакции (рис. 6, 7).

 

Рис. 6. Результаты расчетов гидро-хемо-механической численной модели при закачке флюида с повышенным содержанием растворенного CO₂ (C^f,CO₂ = 1.38 мас. %). (а) – пористость, (б) – эффективное давление, (в) – поток Дарси и (г) – массовая концентрация CO₂ в твердой матрице в разные моменты времени.

 

Наблюдается волнообразное распространения аномалии пористости по механизму вязкой (де)компакции, сопровождающееся увеличением скорости потока Дарси. Также происходит формирование и продвижение резкого фронта реакции.

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

 

Рис. 7. Результаты расчетов гидро-хемо-механической численной модели при закачке флюида с повышенным содержанием растворенного CO₂ (C^f,CO₂ = 1.38 мас. %). Количественный минеральный состав породы на (а) момент начала расчетов, т. е. закачки флюида и (б) финальный момент расчетов, соответствующий образованию фронта реакции. (в) Равновесное соотношение между концентрацией CO₂ во флюиде и твердой матрицей (красная линия), ромбами показаны комбинации C^s,CO₂ и C^f,CO₂, которые формируются в ходе расчета. В данном случае задействован сильный скачок на изотерме, что отражается в формировании резкого фронта реакции. Антигорит полностью исчезает из системы, серпентинит трансформируется в магнезит и тальк.

 

ДИСКУССИЯ

В современной литературе представлены многочисленные примеры моделирования сопряженных процессов течения, деформаций, тепло-массопереноса с учетом химических реакций как в одномерном, так и в двух- и трехмерном случаях, в том числе и не использующие приближение Буссинеска (Roded et al., 2018; Xu et al., 2005; Ghorbani et al., 2016; Pesavento et al., 2016; Katz, 2008; Keller, Suckale, 2019; Gerya, 2019 и ссылки в них них). В работах (Lacinska et al., 2017; Klein, Garrido, 2011; Picazo et al., 2020; Malvoisin et al., 2015; Plumper et al., 2017; Beinlich et al., 2020; Vrijmoed, Podladchikov, 2022) изучались такие же или близкие к рассмотренным реакции дегидратации и карбонатизации в качестве примера. Новизна представленного подхода в настоящей работе заключается в добавлении большой вязкой объемной деформации к методу, представленному в нашей серии работ, учитывающих большие изменения плотности и пористости в реагируюшей системе минеральная матрица–поровый флюид (Malvoisin et al., 2015; Plumper et al., 2017; Beinlich et al., 2020; Vrijmoed, Podladchikov, 2022). В то время как в текущей литературе плотности часто для простоты предполагаются постоянными, в этой серии работ плотности берутся из термодинамических расчетов, делая их согласованными с изменяющимся химическим и фазовым составом системы. Например, в недавно опубликованной работе, посвященной моделированию магматических процессов, связанных с неизбежными эффектами больших изменений плотностей как минералов твердой фазы, так и расплава, авторам пришлось признать, что в приближении постоянства плотностей требуется замена закона сохранения массы на несуществующий закон сохранения объема (Hu et al., 2022), что приводит к термодинамической несогласованности исследуемой модели. В опубликованной нами серии работ пористость вычисляется из безусловно корректного закона сохранения массы химически инертного компонента системы без упрощающего предположения о постоянстве плотности, но предполагая для простоты отсутствие деформации твердого скелета:

$\mathrm{\varphi }=1-\frac{{\rho }_0^s\left(1-{\mathrm{\varphi }}_0\right)C_0^{s,\mathrm{MgO}}}{{\rho }_0^s{C}_{}^{s,\mathrm{MgO}}}$                                                                 (21)

В работе (Bessat et al., 2022) рассмотрен случай отсутствия химически инертного элемента и деформирующийся твердый скелет, но потеряна возможность аналитически проинтегрировать закон сохранения массы какого-либо компонента для получения аналогичной алгебраической формулы для вычисления пористости. В настоящей работе мы вернулись к случаю наличия инертного элемента, но сохранили большую вязкую деформацию по сравнению с (Bessat et al., 2022), необходимую для спонтанного формирования каналов в двух- и трехмерных задачах (Rass et al., 2018; Bessat et al., 2022). Обобщением алгебраической формулы для пористости стало уравнение (19).

Два рассмотренных нами примера с небольшой разницей в составе закачиваемого флюида, но с максимально возможной разницей в морфологии зоны реакции, меняющейся от диффузионной зоны до очень резкого фронта, находятся в полном соответствии с классическими результатами о важности отклонения зависимости концентраций флюида и твердого скелета от линейной (Fletcher, Hoffman, 1974; Vrijmoed, Podladchikov, 2022). Линейная зависимость или постоянные константы равновесия между фазами, характерные для идеальных растворов или простых примесей, приведут к линейной адвекции любых неоднородностей химических составов без изменения их формы в пространстве. Смена знака кривизны в функциональной зависимости концентрации компонента во флюиде от концентрации этого же компонента в твердом скелете, характерная для главных компонентов и при смене фаз в равновесии, определяет морфологию фронта реакции. На рис. 3в кривизна меняется при 5 мас. % CO₂ в системе или около 1 мас. % CO₂ во флюиде. В представленных численных примерах происходит смена морфологии фронта с диффузионного (волна разрежения в газовой динамике) до резкого фронта (ударная волна или скачки в нелинейных кинематических волнах, см. Уизем, 1977; Бхатнагар, 1983; Orr, 2007; Panfilov, 2018) при увеличении концентрации CO₂ в закачиваемом флюиде с 0.44 до 1.38 мас. %.

ВЫВОДЫ

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

Мы представили 1D-численную реализацию на примере карбонатизации серпентинита при фильтрации флюида H₂O-CO₂ в комбинации с вязкой деформацией минеральной матрицы, рассмотрев систему, состоящую суммарно из MgO-SiO₂-H₂O-CO₂.

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

В результате непрерывного процесса фильтрации флюида H₂O-CO₂ с низкой концентрацией растворенного CO₂ (начиная примерно с 1 мас. % CO₂ во флюиде) наблюдается карбонатизация гидратированного серпентинита, сопровождающаяся исчезновением антигорита и появлением магнезита и талька.

Благодарности. Авторы выражают благодарность Х. Фримоду, а также Т.В. Гере, В.И. Мальковскому и Л.Я. Арановичу за тщательное рецензирование, ценные замечания и рекомендации по оформлению, которые позволили повысить качество работы.

Источники финансирования. Работа выполнена при поддержке гранта Министерства науки и высшего образования РФ № 075-15-2022-1106.

About the authors

L. A. Khakimova

University of Lausanne; Skolkovo Institute of Science and Technology; Lomonosov Moscow State University

Author for correspondence.
Email: liudmila.khakimova@unil.ch
Russian Federation, Lausanne, Switzerland; Moscow; Moscow

Y. Y. Podladchikov

University of Lausanne; Lomonosov Moscow State University

Email: liudmila.khakimova@unil.ch

Механико-математический факультет

Russian Federation, Lausanne, Switzerland; Moscow

References

Supplementary files