基于DEW模型的地球深部流体组成与水岩相互作用计算方法综述

Calculation methods for fluid composition and water-rock interaction in the deep Earth based on DEW model—a review

摘要

水岩相互作用会导致流体中元素的价态、赋存形式等发生改变,进而对元素的富集成矿、循环通量等产生影响。由于地球深部样品与实验数据有限,建立和使用地球深部流体模型可以有效地增加人们对深部流体及水岩相互作用的认知。Deep Earth Water(DEW)模型是一种描述地球深部流体热力学性质的数据库,其可以与矿物热力学数据库联用,实现对地球深部水岩相互作用过程的模拟研究。本文阐释了使用DEW模型描述深部流体的必要性,叙述了应用DEW模型进行深部流体物种和水岩相互作用计算的基本原理,介绍了一种基于DEW模型计算流体物种的软件——FluidsLab,列举了地球深部流体以及水岩相互作用的应用案例与研究现状,最后对DEW模型后续的应用与发展方向进行展望。

关键词

水岩相互作用 / 俯冲带 / 地球深部流体 / 热力学 / DEW模型

Key words

water-rock interaction / subduction zones / deep Earth fluids / thermodynamic model / DEW model

0 引言

地球上广泛分布的液态水是地球区别于太阳系其他行星的一个重要特征,水在地球各圈层之间的运移影响了地球内部的动力学过程以及储库之间物质与能量的交换^[1-2]。水在运移的过程中会与其流经的岩石体系发生相互作用,这一过程就被称作水岩相互作用,其对元素成矿、构造运动、生态系统演变等多尺度地质作用过程均有重要影响^[3]。同时水岩相互作用在现代工农业生产过程中被广泛应用,例如:矿物洗选加工、污水处理、土壤修复等^[4⇓⇓-7]。因此,水岩相互作用具有极其重要的研究价值。

水岩相互作用包含流体流经岩石过程中所产生的所有变化,其中既包括岩石对流体组分的改变,又有流体对岩石体系的影响。水岩相互作用是一个动态而复杂的过程,可以分为3个步骤:吸附、溶解、沉淀^[8]。在流体的影响下,矿物主体会发生蚀变,首先矿物的表面开始吸附溶剂分子,使得矿物外侧逐渐水化(图1)^[9];矿物外侧的元素会相继以离子、离子对或络合物的形式离开矿物进入流体^[10],即矿物发生溶解,同时也使得流体中溶质的组成发生改变(图1);随着溶解的持续进行,流体会达到饱和,若外部条件改变导致元素溶解度降低或溶剂减少,溶质则会发生沉淀,流体中溶质与固体矿物之间的交换过程始终保持动态平衡^[11]。水岩相互作用广泛存在于地下水径流^[12]、海底热液系统^[13]、俯冲带等环境中^[14-15]。目前对地球表层流体的组成以及演化过程可直接观测,因此人们对于浅部流体(海水、河流、地下水等)的地球化学性质有很好的限定^[16],其水岩相互作用过程可以通过PHREEQC等程序进行计算^[17]。而对于地球深部流体,特别是俯冲带中的流体,虽然其对地球的演变有非常重要的影响,但由于其温度和压力较高,人们对其认知相对匮乏。研究地球深部流体主要有以下3种方法:天然样品分析、高温高压实验、理论计算模拟。这3种研究方法是相辅相成的,通常来说,深部来源的流体包裹体样品数量少,并且会被后期地质过程改造^[18],而高温高压实验观测是对深部天然样品不足的有效补充,理论计算则是对已有高温高压实验结果与天然样品认知的拓展。因此建立有效的热力学模型是深部流体研究的一个重要目标。目前应用比较广泛的深部流体模型包括C-O-H模型以及DEW(Deep Earth Water)模型。

根据天然流体包裹体的认知,人们首先建立了C-O-H模型来描述深部地质流体^{[19⇓⇓⇓-23]},其在使用过程中可简化为6个流体物种:H₂O、CO₂、CH₄、CO、H₂和O₂^[24-25],并使用非线性方程组求解的方法来计算C-O-H流体的组成^[25]。虽然C-O-H模型可以解释一些地质流体的形成与演化,但后续的研究表明C-O-H模型在使用的过程中具有一定局限性:(1)高温高压实验以及天然流体包裹体的研究均表明矿物的溶解是普遍存在的现象^[26-27],即流体中存在由主要造岩元素、微量元素等组成的离子或者络合物,但C-O-H流体中不包含这些组分,因此其无法对这些情况进行量化^[28];(2)C-O-H流体组分求解过程中需要已知某些与体系组成有关的变量,例如:流体包裹体中常以石墨与流体的平衡共生作为流体中碳元素饱和的标志(此时碳元素活度为1),但是对于无定型石墨或碳元素未饱和时^[29],碳元素活度是由推测来建立的,因此计算结果具有不确定性^[25];(3)C-O-H流体模型与岩石体系的联系仅限于氧逸度与碳元素活度,因此本质上C-O-H流体模型不具备与岩石体系进行有效且完整交互的功能,考虑到实际矿物的溶解、蚀变,流体中主要元素的迁移、pH变化等因素,在描述流体的模型中引入离子与络合物是必要的^[30]。DEW模型是一种描述富水流体中流体物种热力学属性的数据库,其全称为Deep Earth Water模型,它包含了数百种分子、离子、络合物等流体物种,可以通过氧逸度(log f_{{\mathrm{O}}_2})、酸碱度(pH)、矿物溶解度等多方面实现与岩石体系的相互作用,因此DEW模型能更加真实地描述深部地质流体^[31],并已在地球深部流体及水岩相互作用的研究中有广泛的应用。目前有两种以DEW模型计算岩石与流体相互作用的方法:矿物溶解沉淀法(与EQs程序联合使用)^[32]和反向计算滞后流体法(与Perple_X程序联合使用)^[33]。本文将对深部流体组成以及水岩相互作用计算方法进行综述,之后介绍一种直接应用DEW模型计算流体物种的程序——FluidsLab,再结合当前的研究现状对地球深部流体以及深部水岩相互作用进行举例说明。

1 DEW模型概述

描述富水流体地球化学性质的理论由Helgeson^[34]在20世纪60年代建立,Helgeson通过4篇论文详细阐释了描述流体物种状态方程的HKF(Helgeson-Kirkham-Flowers)方法^{[35⇓⇓-38]}。后续研究对HKF方法进行改进并且开发了流体物种状态方程参数估计方法,从而使得改进的HKF方法可以对更多流体物种的热力学属性进行限定^{[39⇓⇓-42]}。水的介电常数是HKF方法中一个关键参数,其可以应用于玻恩方程计算溶解作用对流体物种吉布斯自由能的贡献^[43]。由于高压条件下水的介电常数未知,起初改进的HKF方法只能用于计算流体物种在小于0.5 GPa下的热力学参数。Pan等^[44]使用第一性原理分子动力学的方法,将水的介电常数拓展到10 GPa,这就促成了HKF方法向更高压力的拓展。最终在使用若干高温高压实验对流体物种热力学数据库进行校正后,Sverjensky等^[31]建立了适用温压范围达到1 000℃、6 GPa,用于描述流体热力学属性的DEW模型。Huang和Sverjensky^[45]又结合新的实验数据对DEW模型数据库中流体物种进行了改进。

DEW模型通过综合流体物种的非溶剂化与溶剂化作用得到任意温压下流体物种的标准偏摩尔吉布斯自由能^[43]。流体物种j在任意温度T以及压力p下的标准偏摩尔吉布斯自由能 G_{j;p,T}^0可以由公式(1)给出,其中p_r与T_r分别为标准压力(1 bar)与标准温度(298.15 K), S_j^0是流体物种j的标准偏摩尔熵、 C_j^0是流体物种j的标准偏摩尔热容、 V_j^0是流体物种j的标准偏摩尔体积。

(1) \begin{array}{l}{G}_{j;p,T}^0=G_{j;p_{\mathrm{r}}{T}_{\mathrm{r}}}^0-S_{j;p_{\mathrm{r}},T_{\mathrm{r}}}^0(T-T_{\mathrm{r}})+\\ {\int }_{T_{\mathrm{r}}}^T{C}_{j;p_{\mathrm{r}}}^0\mathrm{d}T-T{\int }_{T_{\mathrm{r}}}^T{C}_{j;p_{\mathrm{r}}}^0\mathrm{d}\mathrm{l}\mathrm{n}T+{\int }_{p_{\mathrm{r}}}^p{V}_{j;T}^0\mathrm{d}p\end{array}

离子在水中的标准偏摩尔体积由两部分构成(公式(2)):非溶剂化作用( \mathrm{\Delta } V_{n,j}^0)和溶剂化作用( \mathrm{\Delta } V_{s,j}^0)。

(2) \begin{array}{c}{V}_j^0=\Delta V_{n,j}^0+\Delta V_{s,j}^0\end{array}

相似的,焓也可以写成非溶剂化作用( \mathrm{\Delta } C p_{n,j}^0)与溶剂化作用( \mathrm{\Delta } C p_{s,j}^0)两部分组成的形式(公式(3))。

(3) \begin{array}{c}C{p}_j^0=\Delta C{p}_{n,j}^0+\Delta C{p}_{s,j}^0\end{array}

根据玻恩方程,溶剂化作用对标准偏摩尔体积以及标准偏摩尔热容的贡献如公式(4)与(5)所示,其中ω_j是玻恩溶剂化系数,ε是水的介电常数。因此要得到高温高压下流体物种的标准偏摩尔属性,需要高温高压下水的介电常数。

(4) \begin{array}{c}\Delta V_{s,j}^0=-{\omega }_{j;p_{\mathrm{r}},T_{\mathrm{r}}}\left(\frac{1}{\epsilon }\right){\left(\frac{\partial \mathrm{l}\mathrm{n}\epsilon }{\partial p}\right)}_T\end{array}

(5) \begin{array}{c}\Delta C{p}_{s,j}^0={\omega }_{j;p_{\mathrm{r}},T_{\mathrm{r}}}T\left(\frac{1}{\epsilon }\right)\left({\left(\frac{{\partial }^2\mathrm{l}\mathrm{n}\epsilon }{\partial T}\right)}_p-{\left(\frac{\partial \mathrm{l}\mathrm{n}\epsilon }{\partial p}\right)}_p^2\right)\end{array}

非溶剂作用对于流体物种热力学性质的贡献也可以用方程表示^[43]。同时考虑非溶剂化作用与溶剂化作用就得到完整的流体物种在溶液中的热力学属性,其中,公式(6)是流体物种的标准偏摩尔体积,公式(7)是标准偏摩尔热容。

(6) \begin{array}{l}{V}_j^0={\alpha }_{1,j}+\frac{{\alpha }_{2,j}}{p+\psi }+\left({\alpha }_{3,j}+\frac{{\alpha }_{4,j}}{p+\psi }\right)\left(\frac{1}{T-\theta }\right)-\\ {\omega }_j{Q}_{T,p}+\left(\frac{1}{\epsilon }-1\right){\left(\frac{\partial {\omega }_j}{\partial p}\right)}_T\end{array}

(7) \begin{array}{l}{C}_{p;j}^0=c_{1,j}+\frac{c_{2,j}}{{(T-\theta )}^2}-\\ \frac{2T}{{(T-\theta )}^3}\left[{\alpha }_{3,j}\left(p-p_{\mathrm{r}}\right)+{\alpha }_{4,j}\mathrm{l}\mathrm{n}\left(\frac{p+\psi }{p_{\mathrm{r}}+\psi }\right)\right]+\end{array} \begin{array}{c}{\omega }_{j}TX+2T{Y}_{T,p}{\left(\frac{\partial {\omega }_j}{\partial T}\right)}_p-T\left(\frac{1}{\epsilon }-1\right){\left(\frac{{\partial }^2{Y}_{T,p}}{{\partial }^{2}T}\right)}_p\end{array}

最终流体物种j在任意的温度T以及压力p下的标准偏摩尔吉布斯自由能如公式(8)所示。

(8)$ \begin{array}{c} \Delta G_{f, j ; p, T}^{0}=\Delta G_{f, j ; p_{\mathrm{r}}, T_{\mathrm{r}}}^{0}-S_{j ; p_{\mathrm{r}}, T_{\mathrm{r}}}^{0}\left(T-T_{\mathrm{r}}\right)- \\ c_{1, j}\left[T \ln \left(\frac{T}{T_{\mathrm{r}}}\right)-T+T_{\mathrm{r}}\right]- \\ c_{2, j}\left\{\left[\left(\frac{1}{T-\theta}\right)-\left(\frac{1}{T_{\mathrm{r}}-\theta}\right)\right]\left(\frac{\theta-T}{\theta}\right)-\frac{T}{\theta^{2}} \ln \left[\frac{T_{\mathrm{r}}(T-\theta)}{T\left(T_{\mathrm{r}}-\theta\right)}\right]\right\} \\ +\alpha_{1, j}\left(p-p_{\mathrm{r}}\right) \\ +\alpha_{2, j}\left[\ln \left(\frac{\psi+p}{\psi+p_{\mathrm{r}}}\right)\right]+\alpha_{3, j} \frac{1}{T-\theta}\left(p-p_{\mathrm{r}}\right) \\ +\alpha_{4, j} \ln \left(\frac{p+\psi}{p_{\mathrm{r}}+\psi}\right) \frac{1}{T-\theta}+\omega_{j}\left(\frac{1}{\varepsilon}-1\right) \\ +\omega_{j ; p_{\mathrm{r}}, T_{\mathrm{r}}}\left(\frac{1}{\varepsilon_{p_{\mathrm{r}}}, T_{\mathrm{r}}}-1\right)+\omega_{j ; p_{\mathrm{r}}, T_{\mathrm{r}}}\left[Y_{p_{\mathrm{r}}, T_{\mathrm{r}}}\left(T-T_{\mathrm{r}}\right)\right] \end{array}$

因此,根据以上的公式推导过程,在已知标准温压下的热力学参数后,还需要对改进的HKF状态方程中与温度以及压力变化有关的参数进行拟合,这些参数是α₁,α₂、α₃、α₄、c₁、c₂和ω。其余公式中的符号为常数^[43]。

2 流体物种的活度模型

水溶液中溶质物种的活度可由半经验的热力学模型计算得到,模型通常是基于Debye-Hückel方程的拓展形式,比如Helgeson形式^[35]与Davies形式^[46]等,二者的近似形式稍有不同,这里我们仅介绍Debye-Hückel方程的Helgeson形式。

对于某一电荷数为Z_j的带电离子j,其活度系数可以由公式(9)表示。

(9) \begin{array}{c}lg{\gamma }_{j;p,T}=-\frac{A_{\gamma }{Z}_j^2{\stackrel{-}{I}}^{0.5}}{1+{\alpha }_j{B}_{\gamma }{\stackrel{-}{I}}^{0.5}}+b_{\gamma ,j}\stackrel{-}{I}+{\Gamma }_{\gamma }\end{array}

其中,等号右侧的第一项为Debye-Hückel项,代表离子间的长程相互作用;A_γ和B_γ为两个系数,其与温度、流体密度、水介电常数有关^[35]。 \stackrel{-}{I}是溶液的离子强度,可以由方程(10)表示。α_j是流体中离子j的平均离子半径,可以用方程(11)计算得到,其中 m_i和 m_j为浓度。

(10) \begin{array}{c}\stackrel{-}{I}=\frac{1}{2}\sum _i{m}_i{Z}_i^2\end{array}

(11) \begin{array}{c}{\alpha }_j=2\times \stackrel{j}{\sum _i}\frac{m_i{r}_i}{m_j}\end{array}

公式(9)中等号右侧的第二项是离子之间短程相互作用的贡献, b_{\gamma ,j}描述溶液中第j个离子与体系的短程相互作用,其是温度以及压力的函数^[35]。公式(9)中的最后一项是转换项,将浓度单位从摩尔分数(mole fraction)转变为摩尔质量浓度(molality)。

溶液中电中性流体物种n的活度系数可以表示为:

(12) \begin{array}{c}lg{\gamma }_{n;p,T}=b_{\gamma ,n}\stackrel{-}{I}+{\Gamma }_{\gamma }\end{array}

其中 b_{\gamma ,n}为带电离子物种与电中性物种之间的短程相互作用。值得注意的是,Helgeson等^[35]并没有提出在不同温压下 b_{\gamma ,n}的数值或者其函数关系,其数值一般来自于实验经验值,而大多数研究认为可以将其定义为0,因此 {\gamma }_{n;p,T}一般认为等于1。

3 流体物种组成的计算

对于只有一个线性独立反应的体系,计算流体组成时,可以采用自由能最小化或者求解非线性方程组的方法,但是对于包含3个以上线性独立反应的体系,若采用自由能最小化方法,是在一个高维曲面上优化出自由能最小点,这种情况下很容易出现优化结果进入局部最小值盆地的情况,得到的自由能并非全局最小值,但这也是最优化研究领域的难点,因此对于具有多个独立反应的体系,推荐采用求解非线性方程组的方法。

3.1 非线性方程组联立法

非线性方程组联立法是指将所有流体物种的浓度作为未知数,通过求解包含体系中所有流体物种的平衡常数方程、元素质量守恒方程以及电荷守恒方程等组成的方程组,进而得到所有流体物种浓度的方法。这就要求所有的方程是线性无关的,同时方程的数目大于等于待求流体物种的数目,但是由于含有对数的方程或者指数的方程是高度非线性的,因此不能确保方程在实数空间内具有唯一解。根据构建质量守恒方程的需要,这种方法需要建立在热力学平衡以及体系封闭或者瞬态具有封闭特征的开放体系的基础之上。为了建立平衡常数方程组,我们需要给定基本流体物种,来构成其他流体物种的形成反应,通常来讲,每一种元素具有一个基本物种,其余流体物种全部写成基本物种组合的形式,下面我们通过一个应用DEW模型的例子来说明^[45]:一个岩石体系中包含的元素有Na、K、Mg、Ca、Al、Si、Fe、C、H和O。这些元素对应的基本流体物种分别如表1所示。

常用的其他流体物种如表2所示。表2中的每一个流体物种的形成反应都可以写成基本物种组合的形式,例如H {\mathrm{CO}}_3^{-},其形成反应可以写作方程(13),进而通过DEW模型获得其形成反应的平衡常数。

(13) \begin{array}{c}{\mathrm{H}}^{+}+C{\mathrm{O}}_3^{2-}=HC{\mathrm{O}}_3^{-}\end{array}

如果要考虑化合价的变化,则需要引入另外一个变量——氧逸度,其可以作为一个被外部限定的量,也可以通过引入氧元素质量守恒来求解。H₂O的含量与活度可以通过两种不同方式来引入:(1)体系为封闭体系,水的质量为1 kg,不随溶解过程而变化;(2)将水作为一个基本流体物种,此时需要额外引入H元素质量守恒方程。但是两种方式均可用方程(14)来计算水活度。

(14) \begin{array}{c}{\alpha }_{{\mathrm{H}}_2\mathrm{O}}=\frac{m_{{\mathrm{H}}_2\mathrm{O}}}{{\sum }_i{m}_i+m_{{\mathrm{H}}_2\mathrm{O}}}\end{array}

3.2 自由能最小化方法

我们以体系中仅包含一个独立化学反应的例子

来介绍自由能最小化方法。考虑化学反应 v_1 \mathrm{A} = v_2 \mathrm{B},流体物种A、B在此条件下的化学势分别为μ_A、μ_B,当这个反应达到热力学平衡时具有:

(15) \begin{array}{c}{v}_1{\mu }_{\mathrm{A}}=v_2{\mu }_{\mathrm{B}}\end{array}

随着反应的持续,物质A与B会发生转变,任一时刻的反应进度为ξ,A、B两个流体相起始含量分别为 n_{\mathrm{A}}^0、 n_{\mathrm{B}}^0,任意时刻,A、B的物质的量可表示为

(16) \begin{array}{c}{n}_{\mathrm{A}}=n_{\mathrm{A}}^0-v_1\xi \end{array}

(17) \begin{array}{c}{n}_{\mathrm{B}}=n_{\mathrm{B}}^0+v_2\xi \end{array}

A、B的化学势可以使用活度对标准化学势进行校正:

(18) \begin{array}{c}{\mu }_{\mathrm{A}}={\mu }_{\mathrm{A}}^0+RTln{\alpha }_{\mathrm{A}}={\mu }_{\mathrm{A}}^0+RTln{\gamma }_{\mathrm{A}}+RTln{n}_{\mathrm{A}}\end{array}

(19) \begin{array}{c}{\mu }_{\mathrm{B}}={\mu }_{\mathrm{B}}^0+RTln{\alpha }_{\mathrm{B}}={\mu }_{\mathrm{B}}^0+RTln{\gamma }_{\mathrm{B}}+RTln{n}_{\mathrm{B}}\end{array}

此时体系的总吉布斯自由能G可以写为物质A、B化学势之和:

(20) \begin{array}{c}G=n_1{\mu }_{\mathrm{A}}+n_2{\mu }_{\mathrm{B}}\end{array}

接下来分别把化学势与物质的量两组方程((15)~(19))代入总吉布斯自由能方程(20)得到一个只与反应进度ξ有关的方程(21),其余符号均为常量或通过其他方式可求的量:

(21) \begin{array}{l}G=n_1\left[{\mu }_{\mathrm{A}}^0+RT\mathrm{l}\mathrm{n}{\gamma }_{\mathrm{A}}+RT\mathrm{l}\mathrm{n}\left(n_{\mathrm{A}}^0-v_1\xi \right)\right]+\\ n_2\left[{\mu }_{\mathrm{B}}^0+RT\mathrm{l}\mathrm{n}{\gamma }_{\mathrm{B}}+RT\mathrm{l}\mathrm{n}\left(n_{\mathrm{B}}^0+v_2\xi \right)\right]\end{array}

体系的自由能会随着反应进度ξ而变化,因此可以得到总吉布斯自由能G(ξ)与ξ的关系,吉布斯自由能的最低点即对应体系热力学平衡时A、B的含量。由此可以推广到包含若干流体物种的情况。可以参考C-O-H模型^[47],以及Perple_X的热力学计算模型,他们均采用的是自由能最小化的方法^[48]。而对DEW模型可使用HighPGibbs软件提供自由能接口的HCh软件来实现热力学计算^[49-50]。

4 水岩相互作用的计算方法

在水岩相互作用过程中,岩石与流体之间会产生相互影响,矿物会被流体蚀变,同时岩石的溶解过程又会改变流体组分。当前,应用DEW模型计算流体与岩石间的相互作用一般有两种方式:(1)以矿物的溶解与沉淀过程为基础的矿物溶解沉淀法;(2)以岩石体系与流体中元素化学势关系为基础的反向计算滞后流体法。需要注意的是,滞后流体法目前无法对矿物相的水化与蚀变过程进行描述,所以本质上反向计算滞后流体法不能完整的计算水岩相互作用过程,而矿物沉淀法的计算流程更加接近于实际的水岩相互作用过程,其既可以量化矿物水化的结果,也可以计算矿物溶解对流体组分的改变。

4.1 矿物溶解沉淀法

矿物与流体的相互作用会使得矿物中的元素以各种流体物种的形式进入流体相,当溶液中各种流体物种的活度商达到矿物的溶解平衡常数时,流体中元素的浓度达到在此温度、压力、氧逸度下的最大值,再增加流体物种浓度时,则会产生沉淀。基于以上矿物溶解平衡的基本过程,可建立矿物溶解沉淀法用以计算水岩相互作用。矿物溶解沉淀法计算流体与岩石的相互作用可以通过EQs程序实现。其基本的计算过程与控制方程如下:

某一时刻溶液或者矿物的元素组成可以分别用一维向量来表示,向量中任意字母均代表一种元素的浓度:

(22) \begin{array}{c}{\stackrel{ \rightharpoonup }{e}}_{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{i}\mathrm{d}}=\left(\begin{array}{ll}{E}_1^f& E_2^f\end{array}\begin{array}{ll} \dots & E_n^f\end{array}\right)\end{array}

(23) \begin{array}{c}{\stackrel{ \rightharpoonup }{e}}_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{l}}=\left(\begin{array}{ll}{E}_1^m& E_2^m\end{array}\begin{array}{ll} \dots & E_n^m\end{array}\right)\end{array}

矿物发生溶解与沉淀的过程,就是对流体元素组成向量与矿物元素组成向量进行线性变换,矿物溶解就在流体组成向量与矿物组成向量相加,沉淀则是相减。对于具有任意元素组成的流体 {\stackrel{ \rightharpoonup }{e}}_{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{i}\mathrm{d}},在一定的温度、压力与氧逸度下,可以通过元素质量守恒方程、流体物种平衡常数方程、电荷守恒方程联立得到在此条件下的流体物种的种类与浓度,见本文3.1节。之后计算矿物沉淀反应的活度商并与其平衡常数进行比较,若过饱和则沉淀矿物,对于任意反应的活度商可以写为方程(24)。其中任意流体物种j的浓度与活度系数分别为m_j和γ_j;反应中的气体物种用k表示,其逸度为f_k。

(24) \begin{array}{c}lgQ=\sum _j{v}_j\left(\mathrm{l}\mathrm{g}{m}_j+\mathrm{l}\mathrm{g}{\gamma }_j\right)+\sum _k{v}_{k}log{f}_k\end{array}

具有n个元素组成的流体,其组成向量指向n维空间中的一点,如果矿物发生溶解与沉淀,流体元素组成就会沿着矿物元素组成向量 {\stackrel{ \rightharpoonup }{e}}_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{l}}发生变化,这样的变化一直进行,直到流体组分不再沉淀矿物,达到新的平衡为止。由此可知,流体组分的变化轨迹会因为沉淀矿物的不同而发生变化。沉淀的矿物通常来讲是在此条件下最稳定、最饱和的矿物,由此可以使用最大梯度下降法找到梯度最大的矿物用以计算矿物沉淀轨迹:首先筛选所有可能发生沉淀的矿物,需要计算所有可能发生沉淀的矿物对流体组成改变的梯度如方程(25)所示,其中 P_{\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}}为沉淀矿物物质的量。

(25) \begin{array}{c}{\nabla }_{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{l}}=\frac{\partial \mathrm{l}\mathrm{g}Q}{\partial P_{\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}}}\end{array}

矿物沉淀后,将会改变流体的元素组成,因此重复流体物种计算与最大梯度沉淀求解过程直到所有矿物不再沉淀为止,此时沉淀的矿物为蚀变产物,流体由于矿物的溶解与沉淀而产生变化,整体迭代计算过程如图2所示。

4.2 反向计算滞后流体法

反向计算滞后流体法^[33]是在简单滞后算法^[46]基础上改进得到的,概括来说,其可以分为两步:“反向”和“滞后”。第一阶段为反向计算电解质物种(back-calculated electrolyte speciation)方法:首先计算出岩石体系中所有元素的化学势,将元素的水解作为固体与流体交互的途径,通过元素的水解反应与流体中的离子或中性流体物种的活度建立联系,获得不同离子之间的活度比值,最后通过电荷守恒,解出所有流体物种的浓度。第二阶段为滞后流体算法(lagged speciation algorithm),其是建立在简单反向计算结果基础上的,使用流体体系部分元素的质量守恒与化学势关系增加了对流体中的溶质与溶剂的内恰,重复迭代直到流体与矿物体系化学势差值达到收敛范围。反向计算滞后流体法求解流体成分可以使用Perple_X软件实现,其具体计算过程如下:

自然界中的常见造岩矿物可以用12种氧化物的成分来描述,因此矿物在流体中的溶解可以用氧化物M_xO_y水解反应的形式来表示,如方程(26)所示:

(26) \begin{array}{c}{\mathrm{M}}_x{\mathrm{O}}_y+2y{\mathrm{H}}^{+}=x{\mathrm{M}}^{2y/x+}+y{\mathrm{H}}_{2}O\end{array}

在反应到达平衡的时候,氧化物的化学势与其对应的流体物种(M^2y/x+)之间的标准偏摩尔吉布斯自由能建立联系,如方程(27)所示:

(27) \begin{array}{l}\mu \left({\mathrm{M}}_x{\mathrm{O}}_y\right)+2y·g\left({\mathrm{H}}^{+}\right)=\\ x·g\left({\mathrm{M}}^{2y/x+}\right)+y·\mu \left({\mathrm{H}}_2\mathrm{O}\right)\end{array}

对于溶质j,其偏摩尔吉布斯自由能可以用方程(28)表示:

(28) \begin{array}{c}{g}_j=g_j^0+RTln{\alpha }_j\end{array}

其中 g_j^0是溶质j的标准偏摩尔吉布斯自由能,α_j是溶质j的活度。对于流体物种,H⁺的标准偏摩尔吉布斯自由能( g_j^0)为0,因此得到方程(29),用以表示离子之间比值关系。

(29) \begin{array}{l}{\mathrm{\Xi }}_{{\mathrm{H}}^{+}}^{{\mathrm{M}}^{2yx+}}=\mathrm{l}\mathrm{n}\left[\frac{\alpha \left({\mathrm{M}}^{2yx+}\right)}{\alpha {\left({\mathrm{H}}^{+}\right)}^{2yx}}\right]=\\ \frac{1}{RT}\left[\mu \left({\mathrm{M}}_x{\mathrm{O}}_y\right)-y·\mu \left({\mathrm{H}}_2\mathrm{O}\right)x-g^0\left({\mathrm{M}}^{2yx+}\right)\right]\end{array}

通过方程(29)计算出流体物种之间的相对含量,最后通过电荷守恒计算出所有流体物种的绝对含量。

虽然第一阶段的简单滞后算法已经获得了流体物种的浓度,但是其缺少溶质成分与流体全部成分之间的联系,也就是缺少了质量守恒。而滞后正向算法通过最小自由能迭代的方法来克服这个缺点:首先通过对全岩体系进行线性自由能最小化得到稳定的矿物相与C-O-H-N-S流体相组合,此时的流体被叫做溶剂。再通过简单反向算法计算与体系平衡的溶质,溶剂和溶质共同组成溶液,此时溶液体系各元素的化学势与溶剂不再相同,新的溶液体系中某一组分的化学势可以表示为方程(30)。

(30) g_{h+1} = \stackrel{t}{\sum _{i=1}} \left(m_h^i\left[g_{h+1}^{\mathrm{*},i}+RT\mathrm{l}\mathrm{n}\left(m_h^i{\gamma }_h^i\right)\right]\right) + \stackrel{t+r}{\sum _{i=t+r}} \left(m_h^i\left[g^{0,i}+RT\mathrm{l}\mathrm{n}\left(y_{h+1}^i{\gamma }_{0,h}^i\right)\right]\right)

溶液中任意元素的物质的量可以由方程(31)表示。

(31) \begin{array}{c}{n}_{j,h+1}=\stackrel{t}{\sum _{i=1}}{m}_h^i{n}_j^i+\stackrel{t+r}{\sum _{i=t+1}}{m}_{h+1}^i{n}_j^i j=1,\dots ,c\end{array}

溶液中的C-O-H-N-S元素浓度由于在第一步热力学相平衡计算中给出,它们在溶液中物质的量要与相平衡结果一致,因此以元素的质量守恒方程(31)限定流体体系的溶质组成变化,直到溶剂与溶质的自由能最小化,即方程(30)取得最小值,得到溶液体系的溶质流体物种浓度。

5 FluidsLab计算程序的使用

我们基于DEW模型开发了一个求解深部流体物种组成的软件——FluidsLab。其目前由两个模块组成:流体为主导的流体物种组成求解模块——FluidsLab_FE和可拓展工具——FluidsLab_Tools。

FluidsLab_FE通过读取用户输入的温度、压力、氧逸度、pH、元素含量,来计算该元素组成的流体在此条件下稳定存在的流体物种及其相应浓度,并可以对计算结果进行可视化。该软件操作界面如图3所示,左侧为输入界面,用以输入元素的浓度、温度、压力、氧逸度等;右侧为结果显示面板,用以显示计算结果中的流体物种浓度以及矿物的饱和状态。其通过求解非线性方程组,得到元素在流体中的存在形式及浓度。考虑到独立方程数与未知数(待求流体物种的数量),有时会出现超定方程组的情况,为了求解器的稳定性和普适性,我们采用了最小二乘法来优化求解非线性方程组。FluidsLab_FE中通常需要限定的外界条件是温度、压力、氧逸度,而当勾选限定pH的复选框后,才会对体系的pH进行外部控制,否则体系的pH由体系热力学平衡自恰得到。FluidsLab_FE可以应用于pH-log f_{O_2}等值线图的绘制、流体包裹体原始流体成分的恢复、高温高压实验流体物种探索等方面。目前仅开放使用部分已进行标定的元素(图3中浅蓝色的元素),更多元素将在实验充分标定后开放使用。

FluidsLab_Tools是FluidsLab为用户提供的独立开发工具,主要提供了FluidsLab_FE中计算流体物种求解器,用户可以直接调用求解器进行二次独立开发;此外我们还提供一些精简版的示例代码,例如计算不同条件下流体物种的存在形式(本文6.1部分)、计算矿物溶解度等等,用户也可以在这些代码的基础上进行修改。所有内容目前已上传至如下仓库 https://github.com/LANCY-universe/FluidsLab/。

6 地球深部流体模型的应用

6.1 地球深部流体中含碳物种的演变

地球深部流体中元素的存在形式,对于理解元素在储库之间的迁移机制和通量具有重要的意义。控制元素存在形式的因素为:温度、压力、氧逸度、酸碱度等,任何一个条件的改变,都将使得流体中碳的赋存形态发生转变。俯冲带流体中的含碳物种主要有CO₂、CH₄,以及以离子形式存在的 {\mathrm{CO}}_3^{2-}和 {{\mathrm{HCO}}_3^{-}}^{}^[26]。随着压力的升高,在5 GPa、600 ℃时,乙酸将作为一种稳定的流体物种出现^[2]。因此,我们使用以上的5个流体物种来探讨含碳流体在5 GPa、600 ℃的条件下碳元素的存在形式。在氧逸度(log f_{O_2})与酸碱度(pH)作为坐标轴的平面上,当体系中溶解的碳元素总浓度为0.2 mol/L时,求解不同pH-log f_{{\mathrm{O}}_2}条件下 {\mathrm{CO}}_3^{2-}、 {\mathrm{HCO}}_3^{-}、CO₂、CH₄和CH₃COO^-相对浓度的变化,5个流体物种之间的转变关系可以构建4个线性无关的方程((32)~(35)),再补充一个碳元素的质量守恒方程(36),可得5种含碳物种在不同pH以及log f_{{\mathrm{O}}_2}下的浓度变化,CH₄与四价碳之间具有与氧逸度相关的转变界面,这个氧逸度转变的区间同时受控于pH,酸性环境相对于碱性环境具有更高的氧逸度转换边界,这可以得出部分条件下,降低pH有利于CH₄的产生。氧化条件下,含碳相受控于pH的变化,随着pH的降低,主要含碳相依次从 {\mathrm{CO}}_3^{2-}变化为 {\mathrm{HCO}}_3^{-}以及CO₂。

(32) \begin{array}{c}C{\mathrm{O}}_3^{2-}+{\mathrm{H}}^{+}=HC{\mathrm{O}}_3^{-}\end{array}

(33) \begin{array}{c}HC{\mathrm{O}}_3^{-}+{\mathrm{H}}^{+}=C{\mathrm{O}}_2+{\mathrm{H}}_{2}O\end{array}

(34) \begin{array}{c}C{\mathrm{O}}_3^{2-}+2{\mathrm{H}}_{2}O+2{\mathrm{H}}^{+}=C{\mathrm{H}}_4+2{\mathrm{O}}_2\left(g\right)\end{array}

(35) \begin{array}{c}C{\mathrm{H}}_{3}CO{\mathrm{O}}^{-}+2{\mathrm{O}}_2\left(g\right)=2C{\mathrm{O}}_3^{2-}+3{\mathrm{H}}^{+}\end{array}

(36) \begin{array}{c}Tota{\mathrm{l}}_{\mathrm{C}}={\sum }_{i=1}^4{m}_i\times N_{\mathrm{C}i}\end{array}

我们的计算结果还表明,虽然乙酸根离子含量很低,在pH-log f_{O_2}浓度曲面上会被其他的流体物种完全遮盖(图4A),但是将Sverjensky等^[2]预测的乙酸活度稳定区域投影到我们计算的乙酸浓度等值线图上,可以发现二者边界形态一致(图4B)。因此在高温高压下乙酸根会作为流体物种稳定存在,但其含量要比其余含碳物种低一个数量级,这也就给出了虽然乙酸分子热力学稳定,但是在实验以及天然样品中很难发现乙酸分子的原因。同时DAC原位实验表明,在压力5 GPa、温度超过600 ℃时,乙酸分子将会转变为石墨质碳^[51]。

6.2 碳酸盐矿物的溶解度

矿物与流体在相互作用的过程中会产生溶解,俯冲带中流体溶解碳酸盐矿物造成的溶解脱碳是一种重要的俯冲板片脱碳机制^[52]。目前矿物与岩石溶解度的原始数据基本来自于高温高压实验,这其中的主要方法有:质量损失法^[53]、原位体积法^[54-55]、流体物种浓度法^[27]、流体金属元素含量法^[55-56]等。对于方解石而言,目前无论是通过活塞圆筒或者原位DAC实验都能对其溶解度有很好的限定,也有很多溶解度以及流体物种的实验数据用于流体热力学模型的校正。但是考虑到白云石在高压下的分解反应(白云石=文石+菱镁矿),DAC实验中白云石会分解为两种矿物,以及白云石的不一致溶解现象,原位体积法不能得到准确的高压下白云石溶解度。通过应用原位DAC多孔实验方法,可以得到文石与白云石溶解过程中相对碳酸根含量的变化^[57],并以此对DEW热力学数据库进行校正。该研究表明原始DEW模型在一些温压下高估了菱镁矿在纯水中的溶解度,通过抑制计算体系中水镁石及其流体物种Mg(OH)₂,可使得计算结果与实验结果符合的很好^{[54⇓⇓-57]}。当流体中Si等其他造岩元素浓度升高时,可能极大的促进碳酸盐的溶解度,并可能使得溶解脱碳成为俯冲带的最主要脱碳方式^[58]。

6.3 地球深部流体与岩石的相互作用

应用DEW模型对水岩相互作用进行计算,可以对深部流体演化与蚀变过程有更加深入的理解。当前关于地球深部流体的水岩相互作用的研究中,主要使用矿物溶解沉淀法(EQ3和EQ6),和反向计算滞后流体法(Perple_X),下面将列举一些应用DEW模型对深部流体以及水岩相互作用进行模拟计算的代表性研究。伴随俯冲深度的增加,洋壳型冷俯冲带榴辉岩中非生物成因烃类将作为主要含碳流体物种的形式存在,俯冲带中非生物烃生成的有利条件为450~560 ℃、1.5~3.5 GPa、氧逸度为FMQ-1到FMQ-3.5,其出现在冷俯冲超高压变质作用的峰期^[59⇓-61]。岩石学与DEW模拟计算结果表明,高压蛇纹石化过程中可以产生大量CH₄与H₂,这与古俯冲带中折返的高压蛇纹岩样品中存在的非生物成因甲烷与氢气的观察一致^[62-63]。俯冲带岩石退变质过程中榴辉岩以及蛇纹岩会发生碳酸盐化作用而固碳^[64-65],退变质时,沉积物产生的含碳流体会在与榴辉岩的相互作用过程中被固定成碳酸盐,榴辉岩的固碳量会受水岩比以及温压、氧逸度等条件影响;而俯冲蛇纹岩会在折返的过程中捕获大量流体中的碳,这些将对俯冲带的碳通量产生重要的影响。应用DEW模型还可以限定金刚石中常见流体包裹体成分端员的成因^[66]:卤水端员(富含Na、K和Cl)来自于橄榄岩的水岩相互作用、硅铝质端员(富含Si和Al)来自于榴辉岩的水岩相互作用、碳酸盐质端员(富含Ca、Mg和Fe)来自于碳酸盐化纯橄岩的水岩相互作用。不同端员的深部流体在迁移过程中经历混合或者再次发生水岩相互作用时,会由于pH降低以及氧逸度的变化而产生金刚石沉淀^[67]。同时DEW模型具有拓展性,可以通过对实验数据的拟合添加新流体物种,比如 {\mathrm{S}}_3^{-}的引入^[68-69],完善了深部流体中S元素的存在形式,进而对俯冲带硫元素的循环有了全新的认识^[70]。

7 结论与展望

DEW模型具有适用温压范围大、元素与流体物种丰富的特点,这使得其具有丰富的使用场景。目前应用DEW模型实现了对深部流体组分以及水岩相互作用过程的热力学计算,其已对深部水岩相互作用的研究产生了重要影响。本文作者也介绍了一种基于DEW模型开发的计算高温高压下流体物种的软件——FluidsLab。同时,作者认为在后续的工作中DEW模型仍需要被完善,遂在此针对研究现状提出如下4点展望:

(1)DEW模型中的热力学数据是通过相对低温低压下流体物种的热力学数据向高温高压下拓展得到的,DEW模型已通过当前的实验数据对大多数流体物种参数进行拟合,但是部分流体物种的适用性未知,因此还需要更加完备的原位高温高压实验来对流体物种的热力学性质进行校正。

(2)水的介电常数是DEW模型向高压下拓展的关键参数,但是对于非极性分子占主导的体系,比如二氧化碳流体等,其介电常数会与富水流体不同,同时溶质离子也会通过破坏水分子氢键网络的形式造成水的介电常数发生改变。因此,建立普适性的溶剂介电常数计算方法尤为重要。

(3)水岩相互作用计算的关键点是流体物种与固体矿物热力学数据的完备性,EQs采用Berman矿物热力学数据库,考虑到关键金属的研究是当前地质学研究的重点,DEW模型已包含多种金属元素的流体物种,但是含有关键金属元素的矿物却很少被包含在矿物热力学数据库中。因此,亟待对可应用于水岩相互作用的矿物热力学数据库进行完善。

(4)高温高压下水岩相互作用过程实际上是一个动态过程,也即很多情况下未达到热力学平衡,而当前的计算方法大多依赖于平衡态计算,为了量化完整的水岩相互作用过程,开发包含动力学过程的非平衡态高温高压水岩相互作用计算方法十分必要。

感谢黄放老师对本文写作思路以及计算方法的指导和帮助。

参考文献

[1] DONG J J, FISCHER R A, STIXRUDE L P, et al. Water storage capacity of the martian mantle through time[J]. Icarus, 2022, 385: 115113.

[2] SVERJENSKY D A, STAGNO V, HUANG F. Important role for organic carbon in subduction-zone fluids in the deep carbon cycle[J]. Nature Geoscience, 2014, 7(12): 909-913.

[3] WANI S P. Natural water remediation: chemistry and technology[J]. Current Science, 2020, 119(12): 2025-2026.

[4] TANG Y G, YANG C W, FINKELMAN R B, et al. Behavior of minerals and trace elements during cleaning of three coals with moderately high ash yields[J]. Energy and Fuels, 2020, 34(2): 2501-2515.

[5] YADAV V B, GADI R, KALRA S. Clay based nanocomposites for removal of heavy metals from water: a review[J]. Journal of Environmental Management, 2019, 232: 803-817.

[6] TONG S, RODRIGUEZ-GONZALEZ L C, PAYNE K A, et al. Effect of pyrite pretreatment, particle size, dose, and biomass concentration on particulate pyrite autotrophic denitrification of nitrified domestic wastewater[J]. Environmental Engineering Science, 2018, 35(8): 875-886.

[7] GODDERIS Y, ROELANDT C, SCHOTT J, et al. Towards an integrated model of weathering, climate, and biospheric processes[M]//OELKERS E H, SCHOTT J. Thermodynamics and kinetics of water-rock interaction. Toulouse: Mineralogical Society of America, 2009: 411-434.

[8] SCHOTT J, POKROVSKY O S, OELKERS E H. The link between mineral dissolution/precipitation kinetics and solution chemistry[M]//OELKERS E H, SCHOTT J. Thermodynamics and kinetics of water-rock interaction. Toulouse: Mineralogical Society of America, 2009: 207-258.

[9] SCHWEDA P, SJOBERG L, SODERVALL U. Near-surface composition of acid-leached labradorite investigated by SIMS[J]. Geochimica et Cosmochimica Acta, 1997, 61(10): 1985-1994.

[10] WESTRICH H R, CYGAN R T, CASEY W H, et al. The dissolution kinetics of mixed-cation orthosilicate minerals[J]. American Journal of Science, 1993, 293(9): 869-893.

[11] POKROVSKY O S, SCHOTT J. Experimental study of brucite dissolution and precipitation in aqueous solutions: surface speciation and chemical affinity control[J]. Geochimica et Cosmochimica Acta, 2004, 68(1): 31-45.

[12] BELKHIRI L, MOUNI L, TIRI A. Water-rock interaction and geochemistry of groundwater from the Ain Azel aquifer, Algeria[J]. Environmental Geochemistry and Health, 2012, 34(1): 1-13.

[13] STEELE-MACLNNIS M, MANNING C E. Hydrothermal properties of geologic fluids[J]. Elements, 2020, 16(6): 375-380.

[14] 李万财, 倪怀玮. 俯冲带脱水作用与板片流体地球化学[J]. 中国科学: 地球科学, 2020, 50(12): 1770-1784.

[15] GUO S, CHU X, HERMANN J, et al. Multiple episodes of fluid infiltration along a single metasomatic channel in metacarbonates (Mogok Metamorphic Belt, Myanmar) and implications for CO₂ release in orogenic belts[J]. Journal of Geophysical Research: Solid Earth, 2021, 126: e2020JB02098.

[16] JING X Y, YANG H B, CAO Y Q, et al. Identification of indicators of groundwater quality formation process using a zoning model[J]. Journal of Hydrology, 2014, 514: 30-40.

[17] CHARLTON S R, PARKHURST D L. Modules based on the geochemical model PHREEQC for use in scripting and programming languages[J]. Computers and Geosciences, 2011, 37(10): 1653-1663.

[18] MAFFEIS A, FERRANDO S, CONNOLLY J A D, et al. Thermodynamic analysis of HP-UHP fluid inclusions: the solute load and chemistry of metamorphic fluids[J]. Geochimica et Cosmochimica Acta, 2021, 315: 207-229.

[19] DUAN Z H, MOLLER N, WEARE J H. An equation of state for the CH₄-CO₂-H₂O style.1. Pure systems from 0 ℃ to 1000 ℃ and 0 to 8000 bar[J]. Geochimica et Cosmochimica Acta, 1992, 56(7): 2605-2617.

[20] DUAN Z H, MOLLER N, WEARE J H. An equation of state for the CH₄-CO₂-H₂O style.2. Mixtures from 50 ℃ to 1000^oC and 0 to 1000 bar[J]. Geochimica et Cosmochimica Acta, 1992, 56(7): 2619-2631.

[21] CONNOLLY J A D. Phase-diagram methods for graphitic rocks and application to the system C-O-H-FeO-TiO₂-SiO₂[J]. Contributions to Mineralogy and Petrology, 1995, 119(1): 94-116.

[22] OHMOTO H, KERRICK D. Devolatilization equilibria in graphitic systems[J]. American Journal of Science, 1977, 277(8): 1013-1044.

[23] FRENCH B M. Some geological implications of equilibrium between graphite and a C-H-O gas phase at high temperatures and pressures[J]. Reviews of Geophysics, 1966, 4(2): 223-253.

[24] HUIZENGA J M. Thermodynamic modelling of a cooling C-O-H fluid-graphite system: implications for hydrothermal graphite precipitation[J]. Mineralium Deposita, 2011, 46(1): 23-33.

[25] HUIZENGA J M. Thermodynamic modelling of C-O-H fluids[J]. Lithos, 2001, 55(1/2/3/4): 101-114.

[26] FREZZOTTI M L, SELVERSTONE J, SHARP Z D, et al. Carbonate dissolution during subduction revealed by diamond-bearing rocks from the Alps[J]. Nature Geoscience, 2011, 4(10): 703-706.

[27] FACQ S, DANIEL I, MONTAGNAC G, et al. In situ Raman study and thermodynamic model of aqueous carbonate speciation in equilibrium with aragonite under subduction zone conditions[J]. Geochimica et Cosmochimica Acta, 2014, 132: 375-390.

[28] 张志刚, 张驰, 耿明. 地幔条件下富水流体状态方程[J]. 中国科学: 地球科学, 2016, 46(5): 569-581.

[29] TUMIATI S, TIRABOSCHI C, MIOZZI F, et al. Dissolution susceptibility of glass-like carbon versus crystalline graphite in high-pressure aqueous fluids and implications for the behavior of organic matter in subduction zones[J]. Geochimica et Cosmochimica Acta, 2020, 273: 383-402.

[30] MANNING C E, FREZZOTTI M L. Subduction-zone fluids[J]. Elements, 2020, 16(6): 395-400.

[31] SVERJENSKY D A, HARRISON B, AZZOLINI D. Water in the deep Earth: the dielectric constant and the solubilities of quartz and corundum to 60 kb and 1200 ℃[J]. Geochimica et Cosmochimica Acta, 2014, 129: 125-145.

[32] WOLERY T J. EQ3/6: a software package for geochemical modeling of aqueous systems: package overview and installation guide (version 7.0)[M]. Livermore: Lawrence Livermore National Laboratory, 1992.

[33] CONNOLLY J A D, GALVEZ M E. Electrolytic fluid speciation by Gibbs energy minimization and implications for subduction zone mass transfer[J]. Earth and Planetary Science Letters, 2018, 501: 90-102.

[34] HELGESON H C. Thermodynamics of hydrothermal systems at elevated temperatures and pressures[J]. American Journal of Science, 1969, 267(7): 729-804.

[35] HELGESON H C, KIRKHAM D H, FLOWERS G C. Theoretical prediction of the thermodynamic behavior of aqueous-electrolytes at high-pressures and temperatures.4. Calculation of activity-coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600 ℃ and 5 kb[J]. American Journal of Science, 1981, 281(10): 1249-1516.

[36] HELGESON H C, KIRKHAM D H. Theoretical prediction of thermodynamic properties of aqueous electrolytes at high-pressures and temperatures.3. Equation of state for aqueous species at infinite dilution[J]. American Journal of Science, 1976, 276(2): 97-240.

[37] HELGESON H C, KIRKHAM D H. Theoretical prediction of thermodynamic behavior of aqueous electrolytes at high pressures and temperatures.2. Debye-Huckel parameters for activity-coefficients and relative partial molal properties[J]. American Journal of Science, 1974, 274(10): 1199-1261.

[38] HELGESON H C, KIRKHAM D H. Theoretical prediction of thermodynamic behavior of aqueous electrolytes at high pressures and temperatures.1. Summary of thermodynamic-electrostatic properties of solvent[J]. American Journal of Science, 1974, 274(10): 1089-1198.

[39] SHOCK E L, SASSANI D C, WILLIS M, et al. Inorganic species in geologic fluids: correlations among standard molal thermodynamic properties of aqueous ions and hydroxide complexes[J]. Geochimica et Cosmochimica Acta, 1997, 61(5): 907-950.

[40] SHOCK E L, HELGESON H C, SVERJENSKY D A. Calculation of the thermodynamic and transport-properties of aqueous species at high-pressures and temperatures - standard partial molal properties of inorganic neutral species[J]. Geochimica et Cosmochimica Acta, 1989, 53(9): 2157-2183.

[41] TANGER J C, HELGESON H C. Calculation of the thermodynamic and transport-properties of aqueous species at high-pressures and temperatures - revised equations of state for the standard partial molal properties of ions and electrolytes[J]. American Journal of Science, 1988, 288(1): 19-98.

[42] SHOCK E L, HELGESON H C. Calculation of the thermodynamic and transport-properties of aqueous species at high-pressures and temperatures: correlation algorithms for ionic species and equation of state predictions to 5 kb and 1000 ℃[J]. Geochimica et Cosmochimica Acta, 1988, 52(8): 2009-2036.

[43] SVERJENSKY D A. Thermodynamic modelling of fluids from surficial to mantle conditions[J]. Journal of the Geological Society, 2019, 176(2): 348-374.

[44] PAN D, SPANU L, HARRISON B, et al. Dielectric properties of water under extreme conditions and transport of carbonates in the deep Earth[J]. Proceedings of the National Academy of Sciences of the United States of America, 2013, 110(17): 6646-6650.

[45] HUANG F, SVERJENSKY D A. Extended deep Earth water model for predicting major element mantle metasomatism[J]. Geochimica et Cosmochimica Acta, 2019, 254: 192-230.

[46] GALVEZ M E, MANNING C E, CONNOLLY J A D, et al. The solubility of rocks in metamorphic fluids: a model for rock-dominated conditions to upper mantle pressure and temperature[J]. Earth and Planetary Science Letters, 2015, 430: 486-498.

[47] ZHANG C, DUAN Z H. A model for C-O-H fluid in the Earth's mantle[J]. Geochimica et Cosmochimica Acta, 2009, 73(7): 2089-2102.

[48] CONNOLLY J A D. Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation[J]. Earth and Planetary Science Letters, 2005, 236(1/2): 524-541.

[49] SHVAROV Y V. HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by windows[J]. Geochemistry International, 2008, 46(8): 834-839.

[50] ZHONG R C, LI Y X, ETSCHMANN B, et al. HighPGibbs, a practical tool for fluid-rock thermodynamic simulation in deep earth and its application on calculating nitrogen speciation in subduction zone fluids[J]. Geochemistry, Geophysics, Geosystems, 2020, 21(5). DOI: 10.1029/2020gc008973.

[51] SZLACHTA V, VLASOV K, KEPPLER H. On the stability of acetate in subduction zone fluids[J]. Geochemical Perspectives Letters, 2022, 21: 28-31.

[52] LAN C, TAO R, ZHANG L, et al. Carbon releasing mechanisms and flux estimation in subducting slabs: problems and progress[J]. Acta Petrologica Sinica, 2022, 38(5): 1523-1540.

[53] CACIAGLI N C, MANNING C E. The solubility of calcite in water at 6-16 kbar and 500-800 ℃[J]. Contributions to Mineralogy and Petrology, 2003, 146(3): 275-285.

[54] LI W C, WANG Q X. In situ determination of magnesite solubility and carbon speciation in water and NaCl solutions under subduction zone conditions[J]. Solid Earth Sciences, 2022, 7(3): 200-214.

[55] FARSANG S, LOUVEL M, ZHAO C S, et al. Deep carbon cycle constrained by carbonate solubility[J]. Nature Communications, 2021, 12: 4311.

[56] FARSANG S, LOUVEL M, ROSA A D, et al. Effect of salinity, pressure and temperature on the solubility of smithsonite (ZnCO₃) and Zn complexation in crustal and upper mantle hydrothermal fluids[J]. Chemical Geology, 2021, 578: 120320.

[57] LAN C Y, TAO R B, HUANG F, et al. High-pressure experimental and thermodynamic constraints on the solubility of carbonates in subduction zone fluids[J]. Earth and Planetary Science Letters, 2023, 603: 117989.

[58] TUMIATI S, TIRABOSCHI C, SVERJENSKY D A, et al. Silicate dissolution boosts the CO₂ concentrations in subduction fluids[J]. Nature Communications, 2017, 8: 616.

[59] ZHANG L J, ZHANG L F, TANG M, et al. Massive abiotic methane production in eclogite during cold subduction[J]. National Science Review, 2022, 10: nwac207.

[60] WANG C, TAO R B, WALTERS J B, et al. Favorable p-T-f_O2 conditions for abiotic CH₄ production in subducted oceanic crusts: a comparison between CH₄-bearing ultrahigh- and CO₂-bearing high-pressure eclogite[J]. Geochimica et Cosmochimica Acta, 2022, 336: 269-290.

[61] BROVARONE A V, MARTINEZ I, ELMALEH A, et al. Massive production of abiotic methane during subduction evidenced in metamorphosed ophicarbonates from the Italian Alps[J]. Nature Communications, 2017, 8: 14134.

[62] PENG W G, ZHANG L F, TUMIATI S, et al. Abiotic methane generation through reduction of serpentinite-hosted dolomite: implications for carbon mobility in subduction zones[J]. Geochimica et Cosmochimica Acta, 2021, 311: 119-140.

[63] BROVARONE A V, SVERJENSKY D A, PICCOLI F, et al. Subduction hides high-pressure sources of energy that may feed the deep subsurface biosphere[J]. Nature Communications, 2020, 11: 3880.

[64] HU H, BROVARONE A V, ZHANG L F, et al. Retrograde carbon sequestration in orogenic complexes: a case study from the Chinese southwestern Tianshan[J]. Lithos, 2021, 392/393: 106151.

[65] PENG W G, ZHANG L F, MENZEL M D, et al. Multistage CO₂ sequestration in the subduction zone: insights from exhumed carbonated serpentinites, SW Tianshan UHP belt, China[J]. Geochimica et Cosmochimica Acta, 2020, 270: 218-243.

[66] HUANG F, SVERJENSKY D A. Mixing of carbonatitic into saline fluid during panda diamond formation[J]. Geochimica et Cosmochimica Acta, 2020, 284: 1-20.

[67] SVERJENSKY D A, HUANG F. Diamond formation due to a pH drop during fluid-rock interactions[J]. Nature Communications, 2015, 6: 8702.

[68] POKROVSKI G S, KOKH M A, GUILLAUME D, et al. Sulfur radical species form gold deposits on Earth[J]. Proceedings of the National Academy of Sciences of the United States of America, 2015, 112(44): 13484-13489.

[69] POKROVSKI G S, DUBESSY J. Stability and abundance of the trisulfur radical ion S₃^- in hydrothermal fluids[J]. Earth and Planetary Science Letters, 2015, 411: 298-309.

[70] LI J L, SCHWARZENBACH E M, JOHN T, et al. Uncovering and quantifying the subduction zone sulfur cycle from the slab perspective[J]. Nature Communications, 2020, 11: 514.

基金资助

科学技术部国家重点研发计划项目(2019YFA0708501)