如何用phreeqc软件来计算朗格利尔饱和指数数
点击联系发帖人
时间：2017-10-17 01:59
该用户的其他资料
(window.slotbydup = window.slotbydup || []).push({
id: '4540180',
container: s,
size: '250,200',
display: 'inlay-fix'
添加成功至
资料评价：
所需积分：5第二讲PHREEQC用户指导_甜梦文库
第二讲PHREEQC用户指导
摘要PHREEQC 第二版是一个用 C 语言编写的计算机程序， 对各种各样低温下的地球化学性 质进行了演算。PHREEQC 是以离子联系的水化学模型为基础的，可以推算（1）生成物和 饱和系数； （2）涉及到可逆反应以及不可逆反应的批反应和一维（1D）运移计算，可逆反 应包括水、矿物/无机溶液、气体、固体溶液、表面络合、离子交换平衡；涉及到的不可逆 反应包括指定成分摩尔转换、动态控制反应、溶液混合和温度变化； （3）逆向模拟实验，其 中多组的无机物和气体摩尔转换以解释在特定成分不确定范围水体之间混合物的不同。 和第一版相比， PHREEQC 第二版的新特点如下：具有在一维运移计算中模拟弥散（或 扩散）和滞流区的能力，用用户确定的速率表达式模拟分子反应，模拟标准的多种成分或非 标准的两种成分的固体溶液的沉淀和溶解， 模拟定体积气相和定压力气相， 考虑表面系数或 交换位置随着无机物的溶解和沉淀或者分子反应的变化而变化， 自动采用多套收敛参数， 打 印用户指定量到原始输出文件和（或）适合输入到扩展表格的文件上，以一种与扩展表程序 更兼容的形式确定溶解成分。 这个版本报告中说明了化学平衡、动态平衡、运移计算以及逆向模拟计算基础方程式， 描述了程序输入，以及举例说明了许多程序功能目的和范围这个报告的目的是说明 PHREEQC 程序的理论和操作，包括组成成分方程的确定，转 换这些方程为数值计算方法的解释，补充数值方法计算机代码组织的描述，程序输入描述， 以及一系列数据组输入和许多论证的程序功能的数据输入和模拟结果的说明。生成方程和正向模拟在报告的这一部分，说明了用于确定水样热动力活动，离子交换物质、表面络合物 质、气相成分、固体溶液和纯相的代数方程式。首先，说明了水样、交换种类和表面性质的 热动力活动和质量作用方程。然后，定义一组函数，用 f 表示，他们必须能同时求解以确 定给定条件下的平衡， 许多这样的方程是各种元素或元素的价电子状态、 交换位置和表面位 置的摩尔平衡方程， 或来自于纯相和固体溶液的质量作用方程。 附加函数用于确定水的碱度、 活度、电荷平衡、气相均衡、离子强度和表面络合均衡。每个函数包含最少量的变量，使函 数的个数等于变量的个数。 程序用一种修改过的牛顿? 罗弗逊方法解答了伴随的非线性方程。 这种方法利用了函数的偏差和每个函数关于主要未知数组或主要未知数的偏导数。 为了清晰 起见，微分计算中用的变量组指主要未知数，每个函数的全导数不用求导说明。在下列方程 中，没有下标或下标“ （aq） ”指液相中的实体， “（e） ”指交换， “（g） ”指气体， “（s） ” 指表面， “（ss） “指固体溶液， “（p） ” 指相的状态。活度方程和质量作用方程在这一部分，定义了水的活度、交换过程和表层物质，并描述了每种物质的质量作用联 系。 质量作用方程是化学系统中根据主要未知数每种物质的摩尔数的质量作用表达式。 然后， 求这些方程关于主要未知数的差分。 每种物质的摩尔数方程和偏导数方程将被组成的摩尔平 衡，电荷平衡和相均衡函数取代。水样PHREEQC 允许单个液相的物质形成或物质平衡。然而多个液相可在一个运行过程中 定义， 一个液相可被定义为一种或多个液相的混合物 （见D数据输入说明‖中的 MIX 关键词） 。 一个例外是液相中溶解物质假定处于热动力平衡状态。 在最初的溶解计算中， 允许氧化还原 元素的化学价不平衡。每个水样 i 的未知数包括活度 ai，活度系数 ri ，重量摩尔浓度 mi， 溶液中摩尔数 ni 。 PHREEQC 根据主要物质改写了所有的化学方程式。其中一个与每个元素（例如，钙离 子 Ca+2）或元素化学价态（例如，三价铁离子 Fe+3）有关的主要水样，加上氢离子活度、水 中电子活度和水的活度。一些程序，例如 MINTEQA2（Allision 和其他，1990）和 MINEQL+ （Schecher 和 Mcavoy，1991） ，把这些物质用术语“成分”表示，但是，这个术语在这儿不 能用，因为易与‘Gibbs’的相规则中的“成分”混肴。PHREEQC 中，每个水样主要物质 用 SOLUTION-MASTER-SPECIES 数据块确定（见“数据输入说明” ） 。数值方法将未知数 的个数减少为主要未知数的最小个数， 并且反复迭代这些未知数的值， 直到得到数学方程组 的解。水溶液的主要未知数是主要物质活度的自然对数、水的活度的自然对数 ? H2O、离子 强度μ 、水溶液中有溶解能力的水的质量 Waq。 下列关系适用于所有水样（除了水成电子和水本身） ：ai=ri mi 和 ni=miwaq. 离子缔合模 型中水样之间的均衡要求满足所有的水样质量方程。例如，水样 CaSO40 的缔合化学反应是 Ca2++SO42-=CaSO40 。 这 个 反 应 在 25 ° C 时 的 logK 是 2.3 ， 形 成 质 量 作 用 方 程 ：10 ?2.3CaSO 0 4 ? 2? ? 2? Ca SO 4?（1）一般情况下，质量作用方程式可被写为：M aq mKi ? ?i? ?m m ,i ,?c（2）其中 Ki 是温度均衡常数， cm，i 是物质 i 中主要种类 m 的化学计量系数，Maq 是水中主要物 质总数，cm，i 的值可正可负。在 PHREEQC 中，缔合化学反应右边各项系数为负，左边各项 系数为正。主要物质采用同样的形式，质量作用方程可简化为 1=am/am。 一种水样的总摩尔数 i 可由质量作用表达式求得：M aq mni ? miWaq ? KiWaq???ic m ,i m.（ 3）牛顿.罗弗逊方法用摩尔数关于主要未知数的全微分。计算如下：M aq mdni ? ni [d ln(Waq ) ? ? cm,i d ln(am ) ?? ln(? i )d? . ??（4）水样的活度系数由戴维斯方程确:? ? ? ? log ? i ? ? Azi2 ? ? 1 ? ? ? 0.3? ? , ? ?或者用扩展的或 Debye-Huckel 方程确定：（ 5） log ? i ? ?Azi2 ? ? bi? , 1 ? B?i0 ?（6）其中 zi 是水样 i 的离子电荷数，A 和 B 是关于温度的常数 。假如 bi 为零，方程 6 是扩展的 Debye-Huckel 方程；假如 bi 不为零，方程 6 是 Debye-Huckel 方程（见 Truesdell 和 Jones， 1974） 。在扩展的 Debye-Huckel 方程中，ai0 是离子大小参数，反之在 WATEQ Debye-Huckel 方程中，ai 和 bi 是适合于 mean-salt 的离子特定参数。否则，如果没有在数据库文件或输入 数据组中说明，Davies 方程用于带电物质。对于不带电物质，活度系数方程中第一项为零， WATEQ Debye-Huckel 方程简化为 Setchenow 方程 （lnri=biμ ） （见 Langmuir， 讨论版 1997） 。 否则，如果没有指定，对所有不带电物质，bi 假定为 0.1。 Davies 方程中，这些活度系数方程关于离子强度的偏微分为：? ? ?? ? 1 ?? ln ? i ? ? ln(10) ? Az i2 ? ? 0 . 3 2 ? ?? , ?? ? 2 ? ? ? 1 ? ?? ???（7）对于扩展的 WATEQ 或 Debye-Huckel 方程为：? ? ? Azi2 ? ln ? i ? ? ln(10) ? bi ? , 2 ? 2 ? B?0 ? ? 1 ? ?? i ? ???（8）PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的温度要 求，每种水样的活度系数参数通过 SOLUTION-SPECES 数据块确定。元素的主要种类和元 素 的 化 合 价 态 用 SOLUTION-MASTER-SPECES 数 据 块 确 定 。 一 种 溶 液 的 组 成 用 SOLUTION 或 SOLUTION-SPREAD 数据块确定（见“数据输入说明” ） 。交换种类离子交换平衡通过交换位置的非齐次质量作用方程和摩尔平衡方程包括在模型中。 PHREEQC 允许多种交换与液相在均衡中共存，用“交换集合”表示。这种方法使用了质量作 用表达式，它以每个交换器中水样和虚构的空的交换位置之间的不完全反应为基础 （Appeio 和 Postma 1993） 。这种空的交换点是交换器的重要种类，其活度的对数是一个附加的主要未 知数。它由 EXCHANGE-MASTER-SPECES 数据块确定（见“数据输入说明” ） 。然而，交 换器的主要物质不包括在摩尔平衡方程中， 使其物质浓度为零。 它的活度在物质上是无意义 的，但是如此却使所有的交换位置充满了其他交换物质。 交换计算的未知数是活度 aie，在 PHREEQC 中定义为当量比值乘以活度系数 ? ie 和交 换器 e 中每种交换物质 ie 的摩尔数 nie 。当量比值是被一种交换物质占据的交换点的摩尔数 比交换点的总数。一种交换物质的活度为 ?ie ? ? iebe,ie nie Te，其中 be，ie 是被交换物质 ie 占据的交换器 e 的当量数， Te 是交换器中交换地点的当量总数。 Te 作为系统中交换器交换地点的 当量总数，不必与每克水的当量数相等，因为系统中水的质量可能多或少于 1kg。根据默认 值，交换物质的活度系数为 1.0，但是可以根据水的离子强度和被交换物质占据的交换地点 当量数选择使用 Davies，扩展的 Debye-Huckel 或 WATEQ Debye-Huckel 活度系数。 水和交换物质之间的平衡要求满足所有交换物质的质量作用方程。交换物质 CaX2 的 缔合反应是 Ca2++2X-= CaX2，其中 X 是默认数据库的交换主要物质。活度的当量比值和这 种化学反应形式的使用统称为 Gaines-Thomas 协议（Gaines 和 Thomas，1953） ，并且在 PHREEQC 版本中用于数据库 phreeqc.dat 和 wateq4f.dat 文件。[在 PHREEQC 版本中，也可 能用 Gapon 协议，也使用当量比值，但是将交换反应写为 0.5 Ca2++X-= Ca0.5X。见 Appeio 和 Postma （1993）再次讨论。]在默认数据库文件中，钙离子交换的 logK 伪.8，形成质量 作用方程：100.8 ??CaX 2 ?Ca 2? ?2 X?,（9）一般情况写，质量作用方程可写为：MKie ? ?ie ? ? m m ,ie?c m,（10）其中 m 随所有的主要物质变化，包括交换的主要物质，cm，ie 是在物质 ie 不完全缔合反应中 主要物质 m 的化学计量系数， kie 是不完全反应选择常数。cm，ie 的值可正可负。PHREEQC 中，缔合反应右半边各项系数为负，左半边各项系数为正。 对于交换物质，物质 ie 总摩尔数方程为：nie ? Kie??mM? c m ,ie m? be,ie ? ie ? ? T ? e? ? ? ?.（11）在数值方法中， 交换器主要物质活度的自然对数是主要未知数。 物质 ie 摩尔数关于主要未知 数的全微分是；?M ? ? dnie ? nie ? ? ? cm,ie d ln(? m ) ? ?? ln(? i )d? ? ?. ? m ?（12）PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的温度要 求以及每种交换物质的活度系数表达式通过 EXCHANGE-SPECES 数据块确定。 交换的主要 种 类 用 EXCHANGE-MASTER-SPECES 数 据 块 确 定 。 交 换 位 置 数 量 和 交 换 器 组 成 用 EXCHANGE 数据块确定（见“数据输入说明” ） 。表面种类表面络合过程通过非齐次质量作用方程、 表面位置的摩尔平衡方程和每个表面的电荷电势联系包括在模型中。 PHREEQC 允许多种表面和表面位置类型液相在均衡中共存， 用 “表 面集合”表示。PHREEQC 表面种类质量作用方程有两个公式： （1）一个包括静电电势项； （2）另一个包括所有的静电电势项。如果使用模型静电电势项，由于存在表面电荷和表面 静电电势，就需用到附加方程和质量作用术语。 交换反应和表面反应计算之间的两个根本不同点在于，交换反应表示为不完全反应， 使主要物质不能出现在任何摩尔平衡方程中， 并且交换物质被期望是中性的 （或不带电的） 。 表面反应不是不完全反应， 所以主要物质是物质上真正存在的物质， 并且可出现于摩尔平衡 方程中，表面物质可带正电，也可带负电，或不带电。 在 Dzombak 和 Morel（1990）中，表面络合反应的基本理论包括静电电势。该理论假 定活性地点个数为 Ts（eq） ，比面积 AS（m2/g） ，表面质量 Ss（g）已知。两个附加的主要未 知数是（1）参量 ln ??s ? ln? e 2 RT ? ?? ? ?F?s? ? ?F?s ，其中 F 是 Faraday 常数（96493.5JV-1eq-1） ， 2 RTΨ s 是表面电势（Volts） ，R 是气体常数（8.3/4）Jmol-1k-1 ） ，T 是温度（kelvin） 。 （2）主要 表面物质活度的自然对数。参量 lnaΨ s 右侧项在术语标准中用 a2 定义。这是用于 Dzombak 和 Morel（1990）中一个不同的未知数，但是由于所有的方程书写和这个主要未知数一致， 所以结果相同。 表面物质活度假定等于给出非空的表面位置类型的摩尔比值。换句话说，一种表面物 质完全覆盖了一种给出的表面位置时是处于标准状态 （活度为 1） 。 这个协议不同于 Dzombak 和 Morel（1990） ，它假定一个表面物质（固相概念上）的活度数值上等于体积克分子浓度 （溶液浓度） 。假如只考虑 monodeneate complexes（如 Dzombak 和 Morel（1990 所做） ，消 去质量作用方程中各项，不考虑标准状态协议，可得到一致的数值结果。然而，体积克分子 浓度 （容摸） 协议用于 multideneate complexes （bidentate， tridentite， 其他， cf.Appelo and postma， 1999)，表面位置浓度存在一个明显的不同。假如一个容器盛有表面包含 multideneate species 且处于均衡的溶液，并且增加了更多同样的溶液，体积克分子浓度不变，溶液的成分和表面 不断改变。这种情况体积克分子浓度明显不正确。 “Hfo‖(含水 Fe2O3) 在默认数据库中用作D_w‖， 指明了较低的亲和力或很弱的位置， D_s‖指明强的亲和力或强的位置。 “Hfo_wOH‖用于说明D在软弱位置‖一种中性表面物质，带 负电软弱位置形成物的缔合反应（在一定意义上，在缔合反应中定义的种类在方程的右侧） 可以被写为： Hfo_wOH→Hfo_Wo-+H+. 质量作用表达式，包括静电电势项项，是 （13）Kint Hfo_wO -?? Hfo_wO - ? H ? ? Hfo_wOHe?F?s RT,（14）其中是反应原有的平衡常数， e?F? s RT是一个因子，说明了将一种带电物质（H+）从带电表面移出。一般，表面物质 i? s k ? 的质量作用方程为： M c m ,i? s ? ? F? ? s k Kiint ? ? ? ? i? sk ? ? ?m ? ?e RT ?zi? sk ? ? sk ? m ? ?（15）其中 K iint 是原有的平衡常数，是表面 s 中表层位置类型 k（软弱的或坚固的 Dzombak 和?sk ?Morel，1990 中）的第 i 层物质；m 随所有主要物质 M 变化，包括表层主要物质； cm , i? s ? 表k层物质 i? s k ? 缔合反应主要物质的化学计量系数， ?zi? s ? 是形成表层物质时表面电荷的净变k化。 cm , i? s ? 可正可负。PHREEQC 中，缔合反应中右边半边各项系数为负，左半边各项系数k为正。 一种表层物质 i(sk)的总摩尔数方程为：ni? s ? ? ? ikTs k bi? s ?k? sk ?? K i Ts k e? sk ? ?? F?s ? ? ? RT ?z i? sk ? ?? M ? ? ???mc m ,i? s ? k m? KiTs k bi? sk ?? sk ?? 2 ?z i? s ? k ?s??mM,（16）c m ,i? s ? k m其中 Tsk 是表面位置类型总数，bi（sk）是与物质相联系表面位置的数量。物质 i(sk)的摩尔数关 于主要未知数的全微分为：?M ? dni? sk ? ? ni? sk ? ?? cm,i? sk ? d ln ?m ? 2?zi? sk ? d ln ??s ? . ?M ?（17）表面物质质量作用方程的第二个公式在质量作用表达式中不包括静电电势项（ -no-edl identifier 在 SURFACE 数据块中） 。一种表层物质摩尔数方程同方程 16，除了涉及α 式没有显示。否则，摩尔数的全微分同方程 17，除了没有最后一项。 PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的表面物 质的温度要求通过 SURFACE-SPECES 数据块确定。表层主要物质或表层位置类型通过 SURFACE-MASTER-SPECES 数据块确定。表层的一致性，每个位置类型的当量数，表面组 成，比表面积和表层质量用 SURFACE 数据块确定（见“数据输入说明” ） 。Ψs的因 气相成分 多成分气相和液相之间的均衡可用非齐次质量作用方程和一个总压力方程（仅为故乡 压力气相）模拟。只有一个气相可与液相在均衡中存在，但气相可包括多种成分。所有的气 体成分假定为理想气体，气相假定为理想的气体混合物。 如果气相体积固定，那么该体积内的压力将随着反应程度变化，但是气相中每种成分 总是存在的。对于一个固定体积的气相，不需要另外的主要未知数，气体中一种成分的摩尔 数可以从液体主要物质的活度中推算出来。 如果气相压力固定，那么气相体积将随着反应程度变化，假如组成成分的部分压力之 和小于指定的总压力，固定压力气相将不存在，并且气相中任一气体成分将不能存在。固定 压力气相中，气体成分的总摩尔数方程中包括一个附加的主要未知数 Ngas。 根据理想假定，一种气体成分的逸度（或活度）等于其部分压力。PHREEQC 使用溶 解方程式，在一定程度上，气体成分假定处于化学反应的左边。CO2 溶解反应可写为 CO2（g）= CO2（aq） 。 （18）亨利定律常数将气体成分的部分压力 （理想气体数值上等于逸度） 与液体物质的活度联系了 起来。对于二氧化碳，亨利定律常数为 10-1.468[假定为理想气体，单位为大气压（atm）]， 下列质量作用方程适用于均衡状态：PCO2 ? 101.468 ?CO2 ? a q? ,（19）其中， P 。一般情况下，气体成分的部分压力可 CO2 是用液相的活度计算的部分压力（atm） 根据液相活度写为：1 aq cm , g Pg ? ? ?m , Kg mM（20）其中，Pg 是气体成分 g 的部分压力，用液相的活度计算；Kg 是气体成分的亨利定律常数； cm，g 是溶解方程中水成主要 m 的化学计量系数。cm，g 值可正可负。PHREEQC 中，溶解反应 方程中右半边各项系数为负，左半边各项系数为正。 对于定体积气相，气相的总体积指定为 Vtotal，但压力是变化的。在平衡状态，气体成 分的摩尔数如下计算：ng ?Vtotal pg RT?Vtotal RTKgM aq m??cm ,g m.（21）气相中该成分的全微分为dng ?Vtotal aq ? ng cm,g d ln ?m RT mM（22）对于定压力气相，气相的总压力指定为 Ptotal，但体积是变化的。在平衡状态，气体成 分的摩尔数等于气体总压力的比值乘以气体的总摩尔数，如下计算：ng ? N gasPg Ptotal?N gas Ptotal K g??mM aqcm ,g m.（23）气相中该成分摩尔数的全微分为dng ?Pg PtotaldNgas ? ? ng cm,g d ln ? m .mM aq（24）PHREEQC 的数据输入、质量作用方程、亨利定律常数以及该常数的温度要求通过 PHASE 数据块确定。气相类型（定体积或定压力） 、气相计算中包括的成分和初始气象成分 通过 CAS_PHASE 数据块确定（见“数据输入说明” ） 。牛顿-罗弗逊方程用一组函数（用 f 表示）用于说明复杂的均衡，这些方程主要来源于用摩尔平衡方程 和电荷平衡方程代换物质摩尔数方程（来源于前一部分的质量作用方程） 。达到平衡时，有 关待定平衡计算的所有函数值为零。函数的零点通过 Tventm-Raohsm 方法获得。据此每个 函数关于每个主要未知数微分形成 Tacobin 矩阵。由 Tacobin 矩阵得到一组线性方程组,其 解可近似为非线性方程的解。 通过反复解答连续的线性方程组， 可得到非线性方程的一个解。 在这一部分说明了每个用于数值方法的函数 f 及其关于形成 Jacobian 矩阵的主要未知数的全 导数水的活度水的活度用根据 Raoult 定律(Gaoult’s and Christ,1965,p.65-66)得出的近似值计算：N aq i? H 2 O ? 1 ? 0.017 ?ni . Waq（25）函数 fH2O 定义为：f H 2O ? Waq ? H 2O ? 1 ? 0.017? ni ,i??N aq(26)该函数的全微分为：dfH 2 O ? Waq ? H 2O d ln ? H 2O ? ? H 2O ? 1 Waq d ln?Waq ? ? 0.017? dni .i?? ??N aq（27）主要未知数是水的活度的自然对数 ln ? H 2 O 。离子强度水溶液的离子强度是一个主要未知数，定义如下：??1 aq 2 ni ? zi W , 2 i aqN（28）函数 f? 定义为：1 aq f? ? Waq? ? ? zi2 ni , 2 i该函数的全微分为：N（29） df? ? ?Waq d ln?Waq ? ? Waq d? ?1 aq 2 ? zi dni . 2 iN（30）定体积多项成分气相平衡对于定体积气相， 每种气体成分的摩尔数可由水中主要物质的活度计算得到， 数值模型 处理气相成分的方法与它处理水中物质的方法不同。 每种气体成份的摩尔数 ng 用在摩尔平 衡方程中用成分形式表示，dng 在 Jacobian 矩阵中用摩尔平衡方程表示。定体积气相均衡计 算不需要附加方程 f。 PHREEQC 数据输入，质量作用方程，Henry 定律常数，以及其温度要求用 PHASES 数据块确定，气相类型（定体积或定压力） ，气相计算中包括的成分和初始气相成分由 GAS_PHASE 数据块确定（见“数据输入说明” ） 。定压力多成份气相平衡定体积气相中每种气体成份的摩尔数由水中主要物质的活度和气相中气体成份总摩尔 数 Ng 决定。每种气体成份的摩尔数 ng 在摩尔平衡方程中用（elements）成份表示，dng 在 Jacobian 矩阵中表示摩尔平衡方程。 定压力多成份气相和液相平衡需要一新的方程――气体 成份部分压力之和等于总压力 Ptotal。函数 f Ptotal 定义为：Ng gf Ptotal ? Ptotal ? ? Pg ,其中，Ng 是气相中气体成份的总摩尔数。（31）f Ptotal 关于主要未知数的全微分（其中正的 dNgas 表示溶液浓度增加）为：dfPtotal ? ??? cm,g Pg d ln ? m .g m N g M aq（32）PHREEQC 的数据输入，质量作用方程，Henry 定律常数以及该常数的要求用数据 块 PHASES 确定。气相类型（定体积或定压力） ，气相计算中包括的成份和初始气相成份用 数据块 GAS-PHASES 数据块定义（见“数据输入说明” ）纯相均衡气相和液相之间的平衡， 包括具有固定部分压力的气体， 通过非齐次质量作用方程包括 在模型中。PHREEQC 允许多种纯相（用纯相集合表示）和液相在均衡中存在，受 Gibbs 纯 相规则的限制，纯相的活度假定一致为 1.0。每种纯相的附加未知数为系统中存在的纯相的 摩尔数 np ,其中 p 指第 p 个纯相，np 表示每个纯相摩尔数的变化，在摩尔平衡方程中用成份 （elements）表示。PHERRQC 也允许如下平衡计算，通过增加或清除一种特定的反应物形 成具有纯相的平衡（EQUILIBRIUM_PHASES 数据块中可供选择的公式或相） ；计算形成纯 相平衡必须的反应物摩尔换算。 在这种计算模型中， 摩尔平衡方程中各项来自反应物的化学 计算而不是纯相的化学计算，未知数是进入或从溶液中析出的反应物摩尔数。 与每个新未知数相应的新函数是每个纯相的质量作用表达式。PHREEQC 中的溶解反应 在一定意义上，纯相是化学方程式左半边的部分。方解石的溶解反应可写为： CaCO3=Ca2++ CO32-, 并且，用 10-8.48 的 logＫ和纯固态活度１.0，质量作用表达式为 （33）K calcite ? 10 ?8.48 ? ? Ca 2? ? CO 2? .3（34） 一般情况下，纯相均衡可用下述方程表示:M aq mK p ? ? ? mm , p .c（35） 其中 cm,p 是溶解反应中主要物质 m 的化学计量系数，其值可正可负。PHREEQC 中溶解反应 左边各项系数为负，右边各项系数为正。无机物 SIp 饱和指数定义为 SIp ? log? ? mm , p .c mM aq(36)数值方法中相均衡函数是f p ? ?ln K p ? ?ln?10??SI p, target ? ? ? cm, p ln?? m ? ,mM aq（37）其中 SIp，target是相的目标饱和指数，ln(10)将以 10 为底的对数转换为自然对数。目标饱和指数由用户指定；正值、零、负值分别表示溶液中无机物过饱和、均衡、不饱和。对固定部 分压力气体成分，SIp，M aq mtarget 等于气体成分部分压力的对数，关于主要未知数的全微分是df p ? ? ? ? cm, p d ln ?m .（38）PHKEEQC 的数据输入，质量作用方程，均衡常数和纯相的温度要求用数据块 PHASES 确定。纯相集合的初始成分和目标饱和指数用数据块 EQUILIBRIUM-PHASES 确定。固体溶液均衡理想多成分或非理想，双成份固体溶液均衡的模拟是以 Glynn 的研究为基础的（Glynn 和 Readon,1990; Glynn 和其他人，1990；Glynn,1991; Glynn and Parklmst,1992） 。液相和固体 溶液均衡通过非齐次质量方程包括在模型中。 PHREEQC 允许多种固体溶液和液相在均衡中 存在，用固体溶液集合表示，受 Gibbs 相规则的限制，非理想固体溶液的模拟只限于双成份 固体溶液； 理想固体溶液可以有两种或多种成份。 固体溶液的附加主要未知数是固体溶液中 每种成份的摩尔数 n p ss ,其中 ss 指固体溶液 ss。该项说明了在元素摩尔平衡 Jacobian 矩阵中 每种成份摩尔数的变化。 不象纯相， 固体溶液中一种成份的活度不完全是 1.0。 其定义为 ? pss ? ? pss x pss ,其中 x p ss 是固体溶液 ss 中成份 p 的摩尔比值， ? p ss 是活度系数，固体溶液中一种成份的摩尔比值定 义为 x p ssn p ssp ss ?1?nN ss, 其中 Nss 是固体溶液 ss 中成份的摩尔数，对于理想固体溶液，活度系p ss数为 1.0，对于非理想双成份固体溶液中各成份的活度系数用 Gugenheim 表达式确定： λ 1=exp(( ?0 - ?1 (4x1-1))x22) 和λ 2=exp(( ?0 + ?1 (4x2-1))x12)，1（39）（40）其中λ 1 和λ 2 是成份 1 和成份 2 的活度系数，?0 和 ? 是无因次的 Gugenheim 参数，无因次 参数由过量自由能 g0 和 g1(kJ/mol)的量纲参数根据方程 ?0 =g0/RT 和 ?1 =g1/RT 计算得到。 过 量自由能参数 ?1 和 ? 2 可以直接定义，也可能通过许多方法间接确定，包括界定可混合性范 围的成份 2 的摩尔比值，界定 spinodal 范围的成份 2 的摩尔比值，临界点和临界温度处的成 份 2 的摩尔比值、Thompson 和 Waldbaum 参数、Mangules 参数，成份 2 的摩尔比值和一个 alyotropic 点的所有溶解物质的对数，成份 1 和 2 微量浓度的固相活度系数，或者成份 2 的 两个分布（配给）系数（glynn,1991） 。 相 应 与 每 个新 未 知 数的新 函 数 是 每个 固 体 溶液中 每 种 成 份的 质 量 作用表 达 式 。 PHREEQC 使用溶解反应，在一定意义上固定溶液成份指化学反应方程式左边各项。对于霞 石（文石）的文石-菱锶矿的固体溶液，溶解反应可写为 CaCO3=Ca2++CO32-， 使用 10-8.43 logK 和固体活度系数，形成质量作用表达式为:的(41)K Arag ? 10?8.34 ??Ca 2? ? CO 2?3? Arag??Ca 2? ?CO 2?3? nArag ? Arag ? ?n ? Arag ? nStront? ? ? ?.（42）一般情况下，固体溶液相均衡中，每种成份可用下列方程说明：M aq mK p ss ???c m , pss mK p ss ? p ss,（43） 其中 K pss 是纯相中成份 p 的均衡常数，cm, p ss 是固体溶液 ss 中成份 p 的溶解反应中， 主要物 质 m 的化学计量系数，其值可正可负。PHREEQC 中溶解反应左边各项系数为负，右边各 项系数为正。固体溶液一种成份的溶解部分（商）可定义为M aqQ p ss ???mc m , pssmK p ss ? p ss,(44)其中平衡时 Qpss 等于 1，ln Qpss 等于 0。对于非理想双成份固体溶液，其数值方法中每种成 分的函数为M aq mf1 ? ? cm,1 ln ?m ? ln K1 ? lnM aq mn1 ? ln ?1 和 n1 ? n2n2 ? ln ?2 . n1 ? n2（45）f 2 ? ? cm, 2 ln ?m ? ln K2 ? ln关于主要未知数的全微分是（46）2 2 3 aq ? x2 2? 0 x2 ? ? 6?1 x2 ? 12?1 x2 ? df1 ? ? cm,1d ln ? m ? ? ? ? ? n ?dn1 ? n ? n m 1 2 ? 1 ? 2 2 3 2? 0 x2 ? 2? 0 x2 ? 6?1 x2 ? 18?1 x2 ? 12?1 x2 ? 1 dn2 n1 ? n2M（47）和df2 ? ? cm, 2 d ln ? m ?mM aq2 2 ? 2? 0 x1 ? 2? 0 x2 ? 6?1 x1 ? 18?1 x12 ? 6?1 x2 ? 12?1 x13 ? 1 dn1 ? n1 ? n2? x1 2? 0 x12 ? 6?1 x12 ? 12?1 x13 ? 1 ? ? ? ?? n ? ?dn2 n1 ? n2 ? 2 ?。（48） 理想固体溶液中，在数值方法中每种成份的函数为 f pss ? ln Q p ss? M aq c m , pss ? ?? m ? ? ln? m K p ss ? ? ?N ss j ss? ? ? n p ss ? ? ?N ? ? ? ? ln? ? total ? ? ? ?，（49）其中 Ntotal ??nss，且 jSS 包括固体溶液 ss 中所有成份。其关于主要未知数的全微分是M aq mdf pss ? ? cm , p ss d ln ? m ?1 n p ss? N total ? n p ss ? ? N total ?N ss, j ss ? p ss ? 1 ? dn ? dn j ss 。 ? ? p ss N j total ? ss（50）PHREEQC 的数据输入， 质量作用方程， 平衡常数， 没仲春相常数的温度要求用 PHASES 数据块定义。一种固体溶液的初始成份和非理想固体溶液的 GUYYENHESIN 参数用 SOLID_SPLUTIONS 数据块确定（见D 数据输入说明” ）表面层的摩尔平衡表层处摩尔平衡是一般摩尔平衡方程的特殊例子。表层集合是由一个或多个表层组成， 每个表层有一个或多个位置类型。一个表层位置类型的总摩尔数由输入指定为如下之一： （1）固定值， （2）与一个纯相的摩尔数成正比，或（3）与动态反应物摩尔数成正比。下列 函数来源于表面 s 的一个表面位置类型 sk 的摩尔平衡关系。一个位置类型的表层物质占据 的表层位置的摩尔数总和必须等于表层位置类型的总摩尔数。 下列函数来源于表面 s 的一个 表层位置类型 sk 的摩尔平衡关系：f s k ? Ts k ? ? bs k ,i? sk ? ni? sk ?i ? sk ?N sk，（51）其中当达到平衡时，函数 f s k 的值为零，Ts k 是表层位置类型的摩尔数， N s k 是某种位置类型 表层物质的数量， bsk ,i? sk?是表层物质 i? s k ? 占据的表层位置的数量， f s k 的全微分为 df s k ? ?Ts k ? ? bs k ,i? sk ? dni? sk ? 。i ? sk ?N sk（52）如果表面位置总数和纯相的摩尔数成正比例，那么 ?TS K ? ?csk , p dnp ，其中 csk , p 是每 摩尔 p 相表层位置的摩尔数。如果相溶解，dnp 为正并且表面点减少。如果点的总数和动态 反应物质成比例，在全微分方程 ?TS K ? 0 中。点数的变化作为反应的一部分并与速率方程 相结合，且 Jacobia 矩阵中不包括任何项。随着动态反应增加或减少，反应物质摩尔数、表 面点数与此成比例调整。假如表面点数固定， ?TS K ? 0 。 PHREEQC 数据输入，表面点每种类型的摩尔总数用 SURFACE 数据块确定，且可能是 一个固定的量或者一个纯相的摩尔数或动态反应有关。表面点类型用 SURFACE_SPECIES 数据块确定并且表层物质可用 SURFACE_SPECIES 数据块确（见“ 数据输入说明‖ ）交换处的摩尔平衡一个交换点的摩尔平衡是一般摩尔平衡方程的特例。 一个交换点的总摩尔数由输入指定 为如下之一： （1）固定值， （2）与一个纯相的摩尔数成正比，或（3）与一个动态反应物摩 尔数成正比。 被交换物质占据点的摩尔数之和必须等于交换点的总摩尔数。 下列函数是一个 交换点的摩尔平衡关系：Ne iefe ? Te ? ?be,ie nie ,（53） 其中达到摩尔平衡时，函数 fe 值为零，Te 是交换器 e 中交换点的总摩尔数， be ,ie 是交换物 质占据的交换点的数量。fe 的全微分方程为dfe ? ?Te ? ?be,ie dnie .ieNe（54）如果点的总数与一个纯相的摩尔数成比例，那么 ?Te ? ?ce, p dnp ，其中 ce，p 是每摩尔 p 相交换点的摩尔数，如果相溶解，那么 dnp 为正并且交换点数减少。如果点的总数与动态反 应物摩尔数成比例，在全微分方程中 ?Te ? 0 。点数的变化与作为反应部分与速率方程相结 合，且没有相包括在 Jacobian 矩阵中。随着动态反应增加或减少，反应物摩尔数、交换点数 成比例增减。如果交换点数固定， ?Te ? 0 。 f e ? Te ??bieNee , ie ienPHREEQC 数据输入， 交换点摩尔数在 EXCHANGE 数据块中确定且可能是一个定量或 者与一个纯相或一个动态反应物的摩尔数有关。交换点用 EXCHANGE_SPECIES 数据块确 定且交换物质用 EXCHANGE_SPECIES 数据块确定（见“ 数据输入说明‖ ）碱度摩尔平衡碱度摩尔平衡方程仅用于物质形成计算和逆向模拟中。它是一般摩尔平衡方程的特例， 系数由每个水样的碱性确定。 PHREEQC 中碱性被确定为一种成分且一种主要物质与其相联 系（见“数据输入说明”中 SOLUTION_MASTER_SPECIES 关键词） 。在 PHREEQC 的默认 数据库中，碱性主要物质是 CO32-。碱度的主要未知数是 lnaAlk ,或默认数据库中 lnaco32-. 碱度当量总数由模型输入确定， 每种水样的碱性影响之和必须等于碱度当量总数。 下列 函数是碱度平衡方程：N aq ief Alk ? TAlk ? ? bAlk,i ni ,（55）其中达到摩尔平衡时，函数 fAlk 的值为零，TAlk 是溶液中碱度当量数，bAlk,I 是水样 i 的碱度 影响，fAlk 全微分方程为df Alk ? ?? bAlk,i dni .ieN aq(56)假定碳酸盐是主要碱性物质，TAlk 的值必须为正。从概念上讲，用 PHREEQC 计算的碱度不 同于实测的碱度。在默认 PHREEQC 的数据库文件中， bAlk,i 的值已选定，每种元素的所指 状态（ bAlk,i ? 0 ）或元素化合价态是 PH=4.5 时的主导物质。假定所有的元素或化学化合价 态在碱性滴定中被转化为这种主导物质。 然而所指状态不存在的水样的有效浓度 （具有非零 碱度的物质）可能在滴定终点存在，这在一定程度上使用 PHREEQC 计算的碱度与所测碱 度是一个不同的量。Fe 和 Al 的氢氧化物络合物是最普通的例子，可能不被转换为所指状 态。这样，一种溶液用 PHREEQC 计算的碱度尽管在数值上等于所测定碱度，实质上也是 一个近似值， 因为是假定滴定过程最终将元素和元素化合价态转变为它们的所指状态。 在许 多溶液中，碱度主要来源于碳酸盐，近似值是有效的。 PHREEQC 的数据输入，以每种物质的缔合反应计算的碱度在 SOLUTION_SPECIES 数据块中确定，主要物质的碱度由 SOLUTION_MASTER_SPECIES 数据块确定，总碱度是 溶液组成的一部分。用 SOLUTION or SOLUTION_SPREAD 数据块确定（见“数据输入说 明” ）元素摩尔平衡系统中一种元素的总摩尔数是下列各相中初始存在的摩尔数总和， 包括纯相、 固体溶液 集合、液相、交换集合、表面集合、气相和表面扩散层。下列函数是总摩尔平衡方程。Np SS N SS E Ne ? ? N aq ? ? f m ? ? Tm ? ? bm, p n p ? ?? bm, p n p ss ? ? ? bm,i ni ? ?? bm,ie nie ? p ss PSS e ie ? ? i??? bs k i ? sk ?SK s N skm , i ? s k ? i ? sk ?n? ? bm, g ng ? ?? bm,i ni , sg s iNgS N aq,（57）其中，达到摩尔平衡时，函数 fm 值为零，Tm 是系统中成份的总摩尔数，Np 是纯相集合中相 的数量，SS 是固体溶液集合中固体溶液数量，Nss 是固体溶液 ss 中成份的数量，Naq 是水样 数量，E 是交换集合中交换器的数量，Ne 是交换位置 e 的交换物质，S 是表层集合中表层个 数，Nsk 是表层类型 sk 的表层物质数量，Ng 是气相成份数量。系统中每个实体的摩尔数分 别用 np 表示纯相集合中相的摩尔数，n p ss 表示固体溶液的成分数， ni 表示水样数， nie 表 示交换位置 e 的交换物质数， ni? s ? 表示表层位置类型 sk 物质数，ng 表示气体成份数， nis 表k示表面 s 扩散层中水样数。每个实体中元素 m 的摩尔数用 bm 表示，用一个附加的下标说明 相关实体，bm 通常但不总是等于 cm（质量作用方程中主要物质 m 的系数） 。为了避免解决大数之间的微小差别，方程 57 中圆括弧中的量没有明确包括在求解算 法 中 ， 且Np pTm的 值 没 有 实 际 计 算 。 函 数SS N ss ss p ssfm中 使 用 代 替 量? Tm ? Tm ? ? bm, p n p ? ?? bm, p ss n p ss 。最初，*Tm 由液相 m 、交换集合、表面集合、气相和表面扩散层中总摩尔数计算得到：? Tm ? ? bm ,i ni ? ?? bm,ie nie ? ??? bm ,i? sk ? ni? sk ? ? ? bm , g ng ? ?? bm, i ni , s 。 （58）i e ie s k i ? sk ? g s iN aqENeSK s N skNgSN aq在方程的反复迭代中，通过纯相和固体溶液成份的摩尔转换使*Tm 得到更新：?Tm k ?1? ?T ? ? bm, p dnp ? ?? bm, p ss dnp ss ，m k p ss p ssNpSS N ss（59）其中 k 指迭代数，在中间迭代中*Tm 可能是负值，但达到均衡时必须为正。函数 fm 的全微分是dfm ? ?? bm, p dnp ? ?? bm, p ss dnp ss ? ? bm,i dni ? ?? bm,ie dnie ?p ss p ss i e ieNpSS N ssN aqENe??? bs k i ? sk ?SKsN skm , i ? sk ?dni? sk ? ? ? bm, g dng ? ?? bm,i dni , sg s iNgS N aq（60）PHREEQC 的数据输入，元素的总摩尔数最初用 SOLUTION or SOLUTION_SPREAD 数据块确定，对于交换集合，用 EXCHANGE 数据块确定，对于表面集合，用 SURFACE 数 据 块 ， 对 于 气 相 ， 用 GAS_PHASE 数 据 块 ， 纯 相 集 合 中 每 个 相 的 摩 尔 数 用 EQUILIBRIUM_PHASES 数据块确定。在一个固体液体集合中每种固体溶液每种成份的摩 尔数 SOLID_SOLUTIONS 数据块确定。 元素的总摩尔数也可在批反应和迁移计算中改修动 （见“数据输入说明” ） 。溶液中电荷平衡电荷平衡方程计算了水中阴阳离子的当量总和。在一些情况下，表面和交换器电荷不 平衡。当具体指定时，在初始溶液计算中用电荷平衡方程调节 PH 值或一种主要物质的活度 （继而一种元素的总浓度或元素的化合价态） 使溶液达到电中性。 电荷平衡方程对于在批量 反应和迁移模拟中计算 PH 值是必要的。 在现实的溶液中， 阴阳离子的当量和必须为零。 然而在化学分析中分析误差和没有分析 的成份通常引起溶液计算中电荷不平衡。 如果一个初始溶液中计算到电荷不平衡， 在后来的 批量反应和运移模拟中需调节 PH 值，维持同样的电荷不平衡。将溶液混合，批反应的电荷 不平衡是每种溶液其混合因素产生电荷不平衡偏差之和。 如果模拟中用到一个表面， 且没有 指定显扩散层计算，那么带电的表层物质形成将产生表面电荷不平衡。相似地，如果交换物 质不是电中性（在默认数据库中所有交换物质是电中性） ，那么交换器将积累电荷。一般的 电荷平衡方程包括表面和交换器电荷不平衡。 在每个初始溶液和批反应步骤中， 对于每个单元在运移模拟的每个时间区间中计算电荷 不平衡使用方程Tz ,q ? ? zi ni ，iN aq（61）其中 q 代表液相， Tz , q 是液相 q 的电荷不平衡，zi 是水样 i 所带电荷。如果带电的表层和交 换器不存在，在溶液批反应中或运移模拟终点的电荷不平衡将和模拟始点相同。 表面的电荷不平衡在初始表面组成计算和每个批反应步骤中计算，对于每个单元运移 模拟的每个时间区间，使用方程Tz , s ? ?? zi? sk ? ni? sk ? ? ? zi ni , s ，k i ? sk ? iK s N skN aq（62）其中 Tz，s 是指表面电荷不平衡, zi? s ? 是表面 s 的表面类型 sk 的表层物质 i 所带电荷，方程中k最后一项代表弥散层中积聚的电荷。 最后一项仅用于弥散层组分明确包括在计算中时 （在数 据块 SURFACE 中的-diffuse-layer） 。详细计算弥散层组分时，要求所有的溶液电荷平衡，并 且 Tz，s 总为零。 正常情况下，交换物质没有静电荷，但是对于大部分，这不作要求。然而，如果带电 物质的和不等于交换点（交换能力）的当量总数，交换物质（当量积分）的活度不能很好确 定。 如果带电物质存在， 那么在每个批反应步骤中的初始交换组分计算中计算交换器的电荷 不平衡，对于每个单元，运移模拟的每个时间区间用方程Tz , e ? ? zie nieeiNe， dfm ? ?? bm, p dnp ? ?? bm, p ss dnp ss ? ? bm,i dni ? ?? bm,ie dnie ?p ss p ss i e ieNpSS N ssN aqENe??? bs k i ? sk ?SK s N skm , i ? sk ?dni? sk ? ? ? bm, g dng ? ?? bm,i dni , sg s iNgS N aq（63） 其中 Tz，e 是交换器的电荷不平衡， zie 是交换器 e 中交换物质 i 所带电荷。 系统电荷不平衡批反应步骤之初确定， 并且对于每个单元， 在运移模拟中每个时间段始 点应该是：Tz ? ? ?qTz , q ? ?Tz , s ? ?Tz ,e ，q s eQSE（64）其中 Tz 是系统电荷不平衡，Q 是在每个批反应步骤或一个运移步骤单元中混合的液相数量，?q 液相 q 的混合比值，S 是表层的数量，E 是交换器数量。电荷平衡函数是：S N aq ? S K s N sk ? E Ne ? f z ? Tz ? ? zi ni ? ??? zi? s ? ni? s ? ? ?? zi ni , s ? ? ?? zie nie ， k k ? s k i ? e i i s i ? sk ? e ? ? N aq（65）其中达到电荷平衡时 fz 为零。 如果计算弥散层组份， 每个表面包括一个独立的电荷平衡方程， 达到电荷平衡时，括号内项目的和将为零。如果弥散层组份不计算，括号内每二项为零。fz 的全微分方程为:df z ? ?? zi ni ? ??? zi? sk ? ni? sk ? ? ?? zie nie ，i s k i ? sk ? e ieN aqSK s N skENe（66）其中，只有在弥散组份没有详细计算时，表面三项的总和存在。 PHREEQC 的数据输入，用 SOLLUTION 或 SOLLUTION-SPREAD、EXCHANGE 和 SURFACE 数据块的数据输入结合物质组成，初始交换组份和初始表面组份计算确定。在 SOLLUTION-SPECIE、EXCHANGE-SPECIE 或 SURFACE-SPECIE 数据块中确定化学物质 反应的物质，电荷由平衡的化学反应确定（见“数据输入说明” ） 。不计弥散层组份的表面电荷电势方程根据缺省值，PHREEQC 应用 Dzombak 和 Morel（1990）说明的方法将表面电荷密度ζ s和表面电势Ψ s 相联系。表面电荷密度是每单位面积表面物质的电荷量，可由表面物质的分布计算而得：F ?s ? Asurf?? zk i ? sk ?K s N ski ? sk ? i ? sk ?n，（67）其中ζ s 是表面 s 的电荷密度（C/m2） ，F 是法拉第常数（96，493.5C/mol） ，Asurf 是物质的表 面积（m2） 。表面积由下列公式之一计算： （1）Asurf=ASSS，其中 As 是表面物质的比表面积 （m2/g） ，SS 是表面物质的质量（g） ，或（2）Asurf=Arnr，其中 Ar 是每摩尔纯相或热动力反 应物的表面积（m2/mol） ，nr 是纯相或热动力反应物的摩尔数。在 25°C 时，表面电荷密度 与表面电势如下：1 1 ? vF?s ? ?s ? ?8000 ?? 0 RT ?2 ???2 sinh? ?， ? 2 RT ?（68） 是自由空间的渗透性（8.854*10-12CV-1m-1 或其中ε 是水的介电常数（78.5，无量纲） ，ε0C2/m-j） ，v 是均质电解溶液的离子电荷，R 是气体常数（8.314Jmol-1K-1） ，T 是温度（K） ， μ 是离子强度，F 是法拉第常数（JV-1eq-1 或 C/mol） ，Ψ s 是用伏特计算的表面电势。在 25° C 时， ?8000 ?? 0 RT ?2 ? 0.1174。电解溶液的离子电荷假定为 1。 电荷-电势函数为：1 1 F ? vF ?s ? f ?s ? ?8000 ?? 0 RT ?2 ?? ?2 sinh ? ?? ? 2 RT ? Asurf1?? zk i ? sk ?K s N ski ? sk ? i ? sk ?n，（69）并且该函数的全微分为：df?s ??8000?? 0 RT ? ?? ?? 1 ? vF?s ? 2 sinh? ?d? ?2 ? 2 RT ?1 2 1 ? 21 2?8000?? 0 RT ? ?? ?F ? F?s ? cosh? ?d ln ? ?s ? ? Asurf ? 2 RT ??? zk i ? sk ?K s N sk。i ? sk ? i ? sk ?（70）nPHREEQC 的数据输入， 默认值不计弥散层组分。 比表面积 （As 或 Ar） 和表面质量 （Ss） 在 SURFACE 数据块中确定，表面位置摩尔数（1）如果数量固定在 SURFACE 数据块中确 定， （2）根据 SURFACE 数据块中的均衡系数，单相的摩尔数在 EQUILIBRIUM-PHASES 数据块中确定，或（3）根据 SURFACE 数据块中的均衡系数，热动力反应物的摩尔数在 KINETICS 数据块中确定。表面物质的电荷在平衡的化学反应中确定，并且该物质在 SURFACE-SPECIESE 数据块中确定（见“数据输入说明” ） 。 计算弥散层组份的表面电荷电势方程作为表面电荷 -电势关系先前模型的替换， PHREEQC 选择使用 Borkovec 和 Westall （1983）发展的方法。他们的发展在氧化电解液界面处求解 Poission-Boltzmann 方程来确定 弥散层中表面的过量离子。假定 1 升水为 1kg，利用下列微分计算。 表面过量离子为：??i , s ?x d ,s? ?c ?x ? ? c ?d xi, s 0 i，（71）其中 ?i , s 为表面 s 上每摩尔每平方米水样 i 的过量离子，xd，s 是外层 Helmholtz 平面的位置， ，作为表面的距离函数，ci0 是溶液的体积浓度。过量离子在 1.0kg ci, s ?x? 是浓度（mol/m3） 水的所指状态中与浓度有关： mi，s=Asurf ?i , s ， （72）其中 mi，s 是水样 i 的过量离子 （mol/kgw） ， 该表面过量离子浓度可在体积溶液中与浓度相联 系： mi，s= gi，smI， 其中 gi，s 是体积溶液中表面电势、所有离子的浓度和电荷的函数： （73）gi , s ? Asurf sign?X d , s ? 1??X d ,s?1?Xzi?1? ?1/ 2? 2 N aq ? z ? X ? ml X i ? 1 ? l ? ? ? ??dX ，（74）其中 X ? e?F?s RT，Xd，s 是外 Helmholtz 平面上 X 的值，Asurf 是表面积（m2） ， （Xd，s-1）为+1或-1，取决于括号中该项的符号，i 是所计算的表层过量离子的水样，zi 水样 i 所带电荷没， l 的范围为所有的水样，zl 是水样 l 所带电荷，并且 ? ? ??? 0 RT / 2?1/ 2。在 25°C 时，的值为 0.02931（l /mol）1/2Cm-2。Borkovec 和 Westall 使用的未知数 X 与 PHREEQC 所用的主 要未知数的关系为a?s ? X ?2 。Borkovec 和 Westall（1983）的发展仅计算了弥散层中每个水样的过量浓度。当一种 溶液从弥散层分离出时，在批反应和运移模拟中就会产生问题，例如，在平流模拟中，水从 一个水槽平流进另一个水槽。在这种情况下，需要知道保留在表面的总水样 i 一个经验假定 是弥散层为指定厚度并且所有的表面过量离子留在弥散层。 弥散层中一个水样总摩尔数等于 弥散层中表面过量离子的有效因素加上体积溶液： ni，s= ni，s，excess+ ni，s，aq= Wbulk+gi，sni n n ? Ws i ? gi , s ni ? Ws i Waq Waq Waq（75）其中 ni，s，aq 指由于体积溶液的存在弥散层中水样 i 的摩尔数， ni，s，excess 指表面过量离子的摩 尔数，Waq 是除了弥散层之外系统中水的质量，Ws 是表面 s 弥散层中水的质量，假定液相中 水量远大于弥散层中水量，使得 Wbulk≌ Waq， （版本 1 中，Wbulk= Waq+ 由弥散层的厚度和表面积计算而得，假定 1 升水为 1kg： Ws=tsAsurf （76） 其中 ts 为弥散层的厚度， 计量单位为米。 如果表层位置的摩尔数与一个纯相或热动力反应物 的摩尔数有关，那么 Asurf=Arnr，否则 Asurf 为常数并且由输入指定的比表面积和表面的质量 计 算而 得。根 据静 电理论 ，弥 散层的 厚度离 子强 度低 时较大 ，离子 强度 高时 较小。 PHREEQC 中弥散层的厚度的默认值为 1*10-8m，大约为由德拜理论计算的 0.001 克分子离 子强度的厚度，弥散层的德拜长计算为 1*10-7 m。如果表面积足够大，弥散层水量很小的假 定将水有效的；厚度为 1*10-7 m，表面积为 1000 m2 将造成弥散层 0.1 升的水量，这是 1 升 体积溶液有意义的一部分。 弥散层水样摩尔数的全微分方程为：? F?s ? ? ? ? ? 2 RT ? ? ?gi , s ? W ? s ? ? ?dni ? ni dni , s ? ? gi , s ? ? 2 e ? ? ? W ? X ? aq ? ? ??2?WsSs） 。水的质量? ?d ln ? ? n Ws d ln W ? n t s Ar dn ?s i aq i r， ? Waq Waq ? ?（77）其中第二项是关于表面电势主要未知数的偏微分，lnaΨ s 和最后一项仅在表面位置数量与一 纯相或热动力反应物有关时存在。偏微分?g i , s 等于方程 74 在 Xd，s 处的积分： ?Xzi d ,s?gi , s ?XX d ,s? Asurf sign?X d , s ? 1???X?1? ?1/ 2? 2 N aq ? z ? X d , s ? ml X d l, s ? 1 ? l ? ? ? ??,（ 78 ） 函数 gi，s 关于主要未知数的偏微分为：?2 ? F?s ? ? ? ? ? ? ?g i , s ?g i , s ? ? ? 2 RT ? ? ? ? 2e ? ? ? ln ? ?s ?X ? ? ? ? ? F?s ? ? 2 ? ? ? ? 2 RT ? ? X dz i, s ? 1 ? ? ? ? ? ? Asurf sign? X d , s ? 1??? 2e 1/ 2 ? ? ? ? 2 N aq ? ? ? ? X ? m X zl ? 1 ? d ,s l d ,s l ? ? ? ???.（79）??在数值方法中， gi，s 的计算很复杂， PHREEQC 版本使用和 Borkovec 和 Westall （1983） 相同的方法以减少函数计算的数量。 当弥散层明确包括在计算中时， 增加了一个新的迭代运 算。在每个弥散层迭代运算开始，函数和它们的偏微分必须明确求出。在弥散层的模拟迭代 中，函数值用下列方程迭代计算：?1 gik, s ? gik, s ??gi , s d ln ??s , ? ln ??s（80）其中 k 是模拟迭代次数，g0i，s 弥散层迭代开始求得的值。当牛顿方法收敛于一个解时，迭代 结束；然而，收敛是以函数 gi，s 的估计值为几基础的。这样，弥散层迭代一直继续到函数值 在指定允许范围内与逐次弥散层迭代相同。 当明确计算弥散层组分时，涉及电势未知数 sinh（方程 69）的函数用包括表面电荷和 弥散层电荷的电荷平衡函数代替：f z , s ? ?? zi? sk ? ni? sk ? ? ? zi ni , s ,k i ? sk ? iK s N skN aq（81）其中当达到电荷平衡时，fz，s 为零。fz，s 的全微分方程为：df z , s ? ?? zi? sk ? d ni? sk ? ? ? zi d ni , s .k i ? sk ? iK s N skN aq（82）PHREEQC 的数据输入，在 SURFACE 数据块中用-diffuse-layer 标识符计算弥散层。比表面 积（As 或 Ar）和表面质量（Ss）在 SURFACE 数据块中确定，表面位置摩尔数（1）如果数 量固定在 SURFACE 数据块中确定， （2）根据 SURFACE 数据块中的均衡系数，单相的摩尔 数在 EQUILIBRIUM-PHASES 数据块中确定，或（3）根据 SURFACE 数据块中的均衡系数， 热动力反应物的摩尔数在 KINETICS 数据块中确定。表面物质的电荷在平衡的化学反应中 确定，并且该物质在 SURFACE-SPECIESE 数据块中确定（见“数据输入说明” ） 。 非静电表面络合Davis 和 Kent（1990）说明了一个非静电表面络合模型。在这个模型中，表面络合质量 作用表达式中静电项被忽略。另外，不使用表面电荷平衡或表面电势关系；仅包括每个表面 位置类型的摩尔平衡方程。 PHREEQC 的数据输入， 在 SUQFACE 数据块中通过用-no-edl 标识符时产生非静电表面 模型（见“数据输入说明” ） 。物质形成和正向模拟的数值方法以前个部分用 PHREEQC 解决的任何化学均衡问题都来自用 f 表示的方程组。 包括 fAlk， fe，fg，fH，fH2O，fm`，fO，fPtotal，fp，fp，fpss，fsk，fz，fz，s，fu 和 fΨ s，其中 fH 和 fO 是氢和氧的 摩尔平衡函数，m`指所有水中的主要物质，除了 H+，e-，H2O 和碱性物质。相应的主要未 知数是， ln ? e ，ng， ln ?e ? ， ln ? H 2 O ， ln ?m` ，lnWaq，Ngas，np（或在物质形成计算中可 能为 ln ?m` ）nss， ln ?sk , ln ? H ? （或在物质形成计算中可能为 ln ?m` ） ， ln ??S （明确的 弥散层计算中） ， ? 和 ln ??S （隐含的弥散层计算中） 。当一给定计算的所有函数的偏差等 于零，即找到了非线性方程组的一个解，化学系统的均衡值已确定。 （注意：如果一纯相或 气相在均衡中不存在，出事包含在给定计算中的一些方程可能会遗漏。 ）解答方法赋予主要 未知数初始值，然后用牛顿-落弗逊方法不断修改直到在指定误差允许范围内找到方程的一 个解。 对于一组方程 fi=0， xj 是未知数， 牛顿-落弗逊方法需不断修改初始未知数组的值。 使 ? i =fi 方程未知数目前值的偏差。形成下列方程组：? i ? ??jJ?fi dxj , ?x j（83）其中 J 是所计算主要未知数的总数。方程组是线性的，可以同时求解未知数 dxj。未知数的 新值由公式计算：k ?1x j ?k x j ? dxj ，k 指迭代次数，然后计算偏差的新值。重复此过程直到偏差小于指定允许误差。 使用牛顿-落弗逊方法 解决化学均衡问题，需注意两点。第一，未知数的初始值必须充 分接近于均衡值，或该方法不收敛，第二，如果一组相的化学反应是非线性相关的，可能形 成 一个奇异矩阵。PHREEQC 使用 Barrodale 和 Roberts（1980）发展的优选技术避免出现 奇异矩阵。 该问题的优选技术也允许不等式约束， 对于限制相和能反应的固体溶液的总两是 有用的。 下面说明了每个模拟类型主要未知数初始值的选择问题。不顾赋予初始估计值的方法， 如果必要， 在牛顿-落弗逊方法产生近似的摩尔平衡之前修改元素和元素价态活度的估计值。 水成主要物质的程序如下。初始估计值确定后，计算每种元素（除了氢和氧）和初始溶液中 确定的各个原子价态的物质分配。然后，求出计算所得摩尔数和输入摩尔数的比率。如果一 种主要物质 m`的比率大于或小于 10-5，用下列方程矫正主要未知数的值：? N aq ? ? ? b ` ni ? ? i m ,i ? , ?1 k ?1 ln ? k ? ln ? ` ` ? w ln m m ? T ` ? m ? ? ? ? ? ?（84）其中如果比率大于 1.5，w 为 1.0，比率小于 10-5，w 为 0.3，Tm`是一种元素或元素原子价状 态的总浓度。交换和表面主要物质使用累似的方程。初始估计值矫正后，计算物质分配。继 续迭代直到比率在指定范围内，在该点，使用修改的牛顿-落弗逊技术。如果连续的修改不 能找到活度使得比率在指定范围内，那么，第二组迭代试图将比率减小到低于 1.5，这些比 率无更小的限制。无论第二组迭代是否成功，都用牛顿-落弗逊技术。 Barrodale 和 Roberts（1980）优选技术是简单线性算法的修正，使受相等和不等约束的 线性方程组偏差绝对值（L1 优选）的和最小。一般问题可用下列矩阵方程表示： AX=B CX=D EX≤F 第一个矩阵方程在一定意义上被减至最小， 使 （85）? b ? ??i i jIBJi, jx j 为一极小值，其中 IB 是 被优选方程的数量，受第二个矩阵方程相等条件约束和第三个矩阵方程不等条件约束。 PHREEQC 方法在优选方程的（AX=B）中将包含一些牛顿-落弗逊方程，而不包括所 有的牛顿-落弗逊方程，如等式（CX=D） ，在给定迭代中，包含在 A 矩阵中的方程可能不能 求得确切的等式，但是在一定意义上将得到优化。这样，在给定迭代中，即使无确切的等式 解存在，也能找到牛顿-落弗逊方程组的一个近似数学解，例如，所有方程的压力等式造成 一个无解的奇异矩阵。矩阵 A 包括碱度方程，气相中气体的总摩尔数，纯相和固体溶液组 分。矩阵 B 包括所有的摩尔平衡，电荷平衡和表面电势方程。矩阵 C 包括限制纯相的溶解、 固体溶液组分和气体成分在系统中存在程度的不等式。 如果相应未知数（在元素主要物质活度的对数中产生变化）的所有系数（一个直列矩阵 A 和 B）都小于 1e-10，为避免一些与矩阵 B 中小数字有关的数值问题，代表摩尔平衡方程 的一行矩阵将被改变（按比例缩放） 。在这种情况下，方程与 1e-10 除以最大系数的绝对值 成比例。当指定时（在 KNOBS 数据块-diagonal-scale 中） ，如果未知数的系数小于 1e-10， 可选择使用另一种方法，摩尔平衡方程将与 1e-10 除以相应未知数的系数成比例。 改变的矩阵用优选的解答器解答，返回值是一主要未知数值的变量。检查变量的值，确 保未知数变量小于限制变量最大值允许的范围的标准。这些标准在在程序中由默认值或在 KNOBS 数据块中的输入指定。如果变量太大，那么所有未知数的变量，除了纯相和固体溶 液组分的摩尔转换，按比例减少，使其满足所有标准。纯相和固体溶液摩尔转换不改变，除 了产生纯相和固体溶液组分总摩尔数的非负值。 计算未知数相适应的变量后， 更新主要未知 数，计算所有水样、交换物质和表面物质新的克分子浓度和活度以及所有函数的偏差 。偏 差用收敛验证（收敛判据在程序中确定，但可在 SELECTED-OUTPT 数据块 KNOBS 或 -high-precion 选择中-convergence-tolerance 转换为一替代组，如果没有收敛，则开始新的迭 代。水成物质计算水成物质计算用溶液的一种化学组分作为输入， 并计算水样的配比和相的饱和指数。 水 成物质计算包括方程 fH2O、fm`、和 fu，即元素或元素化合价态的摩尔平衡方程、水的活度和 离子强度方程。不包括氢和氧的摩尔平衡方程，因为一般不知道氢和氧的总质量。相反水的 质量假定为 1.0kg 或指定（在 SOLUTION 或 SOLUTION-SPREAD 数据块-water 中） ，氢和 氧的总质量在 水成物质计算中由水的质量和包含水样的氢和氧的浓度计算得到。 如果 PH 值、pe 或一种元素或元素化合价态的主要未知数指定为调整到溶液的电荷平 衡，fz 用于计算产生电荷平衡的主要未知数（ ln ? H ? 、 ln ?e ? 、 ln ?m` ）的值。在这种情况 下，计算得到的 PH、pe 或 m`的总浓度将不同于输入值。如果用 fz 计算主要未知数 ln ?m` ， 则不用方程 fm`。 如果输入中指定总碱度，用碱度摩尔平衡方程计算 ln ? Alk 和与碱度有关的元素（默认 数据库中为碳）的总克分子浓度。如果定义包括一个碳[碳（+4）]和碱度两者的摩尔平衡方 程，那么与这些方程有关的两个主要未知数是 ln ? Alk ? ln ?CO ?2 （默认数据库文件中）3ln ? H ? 。在这种情况下，PH 值将在水成物质计算得到并且不等于输入的 PH 值。对于物质形成计算， 如果在问题的阐述中包括碱度摩尔平衡方程， 则它为解答器唯一的 优选方程。所有其他的方程作为相等约束条件，物质形成计算不包括不等约束条件。 在 初 始 溶 液 计 算 中 允 许 部 分 氧 化 还 原 不 平 衡 反 应 ， 在 SOLUTION 或 SOLUTION-SPREAD 数据块中氧化还原反应的选择影响回程物质形成和饱和指数计算。根 据默认值，无论何时需要用电子活度计算一个水样的克分子浓度和活度，使用输入的 pe。 如果给出缺省的氧化还原对（-redox）或指定某元素的氧化还原对（或元素化合价态的组合） （见数据输入说明中的关键字 SOLUTION） ，然后重写每个水样氧化还原元素的质量作用表 达式，用氧化还原对的活度代替电子的活度。例如，使用硫酸盐-硫化物氧化还原对[S（+6） /S（-2）]铁离子将被改变，Fe+3 原始的化学反应： Fe+2 =Fe+3+e用硫的缔合反应写为， SO4-2+9H++8e-=HS-+4H2O， 生成下列不含电子的化学反应： （87） （86）1 ?2 9 ? 1 1 Fe ? 2 ? SO4 ? H ? Fe ? 3 ? HS ? ? H 2O . 8 8 8 2（88）最后的质量作用表达式将被用作 Fe+3 的质量作用表达式，Fe+3 摩尔数的微分 dnFe?3 ，也将以 这个质量作用表达式为基础。然而，原始的质量作用表达式（以方程 86 为基础）用于确定 含有 dnFe?3 的摩尔平衡方程，也就是说，Fe+3 将出现在铁的摩尔平衡方程中，而不出现在 S （+6）或 S（-2）的摩尔平衡氧化还原反应中。另一组方程中。这些变换的结果是 Fe+2 、 Fe+3、硫酸盐和硫化物处在平衡的氧化还原反应元素（例如，氧和氮）也可以形成在它们自 己的氧化还原反应平衡中，而不必与铁和硫一起生成氧化还原平衡反应。 根据默认值，如果饱和指数计算需要 pe（或电子的活度）的值，那么，用 pe 的输入值。 如果已确定氧化还原对（-redox） ，那么该相的溶解反映写法如上，消除 电子活度而用氧化 还原对的活度代替。 在一个计算中，主要未知数组可能随氧化还原元素的变化而变化。如果作为一个摩尔平 衡方程的主要未知数的主要物质活度变为小于同一摩尔平衡方程另一主要物质活度数量级 的十倍，基础转换过程将发生。在这种情况下，涉及目前主要未知数的质量作用表达式（包 括水样、交换物质、气体、表面物质和纯相）被写在具有更大活度的新的主要物质的项中。 这种过程的另一个离子为，如果氮在一个系统中减少，那么氮的主要未知数将从硝酸盐（在 氮减少的条件下，其量可忽略）转换为铵（此时将作为主要物质在） 。基础转换过程不影响 物质的最终平衡分配，但是在处理低浓度问题时，加快了计算速度并避免了数值问题。 估计主要未知数的初始值并且根据前一部分说明的方法修正。对于初始溶液计算，PH 值、pe 的输入值用作初始估计。水的质量为 1.0kg，除非有指定值，水的活度估计为 1.0。 假定主要物质是存在的唯一物质，并且它们的浓度等于输入浓度（转换为克分子浓度） ，估 算离子强度。元素（除了氢和氧）和元素化合价态主要物质的活度恒等于输入浓度（转换为 克分子浓度） 。如果用电荷平衡方程和相均衡方程代替一种元素和元素化合价态的摩尔平衡 方程，那么主要物质的初始活度恒等于输入浓度（转换为克分子浓度）的千分之一。 PHREEQC 的数据输入，物质形成计算使用的碱度方程、电荷平衡方程、相均衡方程和 氧化还原对在 SOLUTION 或 SOLUTION-SPREAD 数据块中确定（见“数据输入说明” ） 。交换器初始组分的计算如果一个交换器的组分没有明确给定，但是指出与一指定溶液组分在均衡中存在，那么 就需要计算初始交换组分。 在这种情况下， 交换器的组分未知， 只知道其与一溶液处于均衡。 初始交换组分计算方程是 fe、fm`、 f H 2 O 、和 fu，即每个交换器的摩尔平衡方程、每种元素和 元素化合价态的摩尔平衡方程、水的活度和离子强度方程。 对于初始交换组分计算，Tm`的值仅包括水样的浓度，fm 不包含交换器对总元素浓度的影 响。所有与液相有关的量和没有交换器存在时的溶液相同。只有交换集合未知数 ln ? e 的值 调整为达到交换器的摩尔平衡 。一旦达到摩尔平衡，即知每个交换器的组分。 所有初始交换组分计算方程作为解答器的相等条件约束。没有方程得到优选且不包括不 等条件约束。 只有交换器与指定溶液处于均衡时，才需计算初始交换组分。该溶液的物质配比已通过 初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关的未知数的值并且作为交 换计算的初始估计值。每个交换器的主要未知数的初始估计值恒等于该交换器交换的摩尔 数。 PHREEQC 的数据输入， 初始交换组分计算用 EXCHANGE 数据块 （见 “数据输入说明” ） 。表面初始组分的计算如果一个表面组分没有明确给定，但是指出与一指定溶液组分在均衡中存在，那么就需 要计算初始表面组分。在这种情况下，表面的组分未知，只知道其与一溶液处于均衡。初始 表面组分计算方程是 f s k 、 f ?s 或 fz，s、 fm`和 fu， 即表面集合中每种表面类型的摩尔平衡方程、 每个表面的电荷电势关系或每个表面的电荷平衡方程（非静电模型中不包括这些方程） 、每 种元素和元素化合价态的摩尔平衡方程、水的活度和离子强度方程。 对于初始表面组分计算，Tm`的值仅包括水样的浓度，相应的摩尔平衡方程 fm`不包含表 面对总元素浓度的影响。所有与液相有关的量和没有表面集合时的溶液相同。 计算弥散层组分时，每个表面使用电荷平衡方程 fz，s；表面集合每个表面类型主要未知 数 ln ?sk 的值调整为达到每个表面的摩尔平衡和电荷平衡 。如果不计算弥散层组分时，那 么用电荷电势方程 f ?s 代替表面电荷平衡方程。如果非静电模型用于表面集合，那么表面电 荷平衡方程和电荷电势方程都不包含在待解的方程组中。 所有初始表面组分计算方程都作为解答器的相等条件约束。没有方程得到优选且不包括 不等条件约束。 只有初始表面与指定溶液处于均衡时，才需计算初始表面组分。该溶液的物质配比已通 过初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关的未知数的值并且作为 表面计算的初始估计值。 每个表面主要主要物质活度的初始估计值恒等于该表面的表面位置 摩尔数的十分之一。 对于明确和隐含弥散层计算， 每个表面电势未知数 ln ??s 的初始估计值 为零，暗示了表面电荷为零。 PHREEQC 的数据输入，初始表面组分计算用 SURFACE 数据块（见“数据输入说明” ） 。固定体积气相初始组分计算如果气相组分没有明确确定，但是一直定体积气相组分与一特定溶液组分处于均衡中 时， 需要计算初始气相组分。 初始气相组分计算方程与初始溶液组分计算方程相同， 是 fm` 、 f H 2O 和 fu，即每种元素和元素化合价态的摩尔平衡方程、水的活度和离子强度方程。对于初始气相组分计算，Tm`的值仅包括水样的浓度，相应的摩尔平衡方程 fm`不包含气 体对总元素浓度的影响。 所有与液相有关的量和没有气相存在时的溶液相同。 一旦确定了液 相中组分的物质配比，气相中所有组分的部分压力即可计算。在气相中根据理想气体定 f ?s 律用部分压力和指定的体积计算每种成分的摩尔数。 所有初始气相组分组分计算方程都作为解答器的相等条件约束。没有方程得到优选且不 包括不等条件约束。 只有气相被确定具有固定体积且与指定溶液最初处于均衡时，才需计算初始气相组分。 该溶液的物质配比已通过初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关 的未知数的值并且作为初始气相组分组分计算的初始估计值。 PHREEQC 的数据输入， 初始气相组分计算用 GAS-PHASE 数据块 （见 “数据输入说明” ） 。批反应和运移模拟批反应和运移模拟计算要求计算液相和化学系统中存在的任何均衡相集合、表面集合、 交换器集合、固体溶液集合和气相的均衡。发生在均衡之前的不可逆反应包括混合、理想配 比反应、热动力反应和温度变化。批反应和运移模拟计算中完整的牛顿-罗弗逊方程组包括 fe， fH， f H 2 O ，fm`，fO， f ptotal ，fp， f p ss ， f s k ，fz，fz，s，fu 和 f ?s 。 总包括氢的摩尔平衡方程 fH、水的活度 f H 2 O 、氧的摩尔平衡方程 fO、电荷平衡方程 fz、 离子强度 fu，并且与总作为主要未知数的 ln ? H 2 O 、Waq（水的质量） 、 ln ? H ? 和 ? 有关。 摩尔平衡方程 fm 用于求元素的总浓度，而非各个化合价态或化合价态组合的浓度。不可 能包括碱度摩尔平衡方程；它用于初始溶液计算。 假如给定定压力气相且在均衡中存在，则包括方程 f ptotal 。假如给定交换集合，则包括方 程 fe。假如给定表面集合，则包括方程 f s k 。另外，假如给定每个表面的隐含弥散层计算， 则包括方程 f ?s 。假如给定明确的弥散层计算，则包括方程 fz，s。包括均衡中存在的每个纯 相的方程 fp。包括均衡中存在的每种固体溶液每种成分的方程 f p ss 。 在计算之初， 已知一纯相、 固体溶液或定体积气相是否在均衡中存在。 这样， 在每次迭 f p ss 代中，用下列推理确定均衡计算应包括哪个方程。如果一个相的摩尔数大于零，np&0, 或者 计算所得饱和指数大于目标饱和指数，则包括其方程。如果矩阵中不包括该方程，那么未知 数 dnp 所有的系数恒等于零。 对于一理想固体溶液，如果任一组分的摩尔数大于 1*10-13 或其和?p ssIAPp ss K p ss大于 1.0，则包括方程 f p ss 。对于一理想固体溶液,IAPp ss K p ss? x p ss 如果摩尔级分的和大于 1.0，则确定求和。如果固体溶液方程不包含在矩阵中，那么未知数 dnPss 所有的系数恒等于零。 对于非理想双成分固体溶液， 决定是否包括固体溶液方程的下列过程由 Glynn 和 Reardon 方程（1990，方程 37 到 48）发展而来。如果任一固体溶液组分的摩尔数大于 1*10 ，则包 括所有的固体溶液方程。否则，该组分水的活度级分由下式计算：-13x1, aq ?IAP IAP 1 2 和 x1, aq ? IAP ? IAP IAP ? IAP 1 2 1 2（89）其中 IAP 是离子活度积。下一步通过求解下列方程得到 x1 和 x2 的值从而确定将与水的活度 级分处于均衡的固体的摩尔级分：x1?1K1 ? x2? 2 K 2 ?1 x1, aq ?1K1 ? x2, aq ?2 K2，（90）其中 x1 和 x2 是在固相中的摩尔级分，K1 和 K2 是单成分的均衡常数， ?1 和 ? 2 是成分的活度 系数，据过量自由能的 Guggenheim 参数计算而来。该方程高度非线性且用第一个测试区间 [0，1]求得一满足方程包含摩尔级分的解，然后使区间缩小一半，使摩尔级分的估计值更精 确。一旦确定固体的摩尔级分， “总活度积”的两个值如下计算：??和aq? IAP1 ? IAP2solid（91） （92）??? x1?1K1 ? x2? 2 K2solid如果??? ?? aq ，那么包括固体溶液方程，否则不包括。如果固体溶液方程不包含在矩阵中，那么矩阵中未知数 dnp ss 所有的系数恒等于零。 在每次迭代中，对于固定压力气相，如果气相的摩尔数大于 1*10-14 或气相成分各部分 压力的和（根据水样的活度计算而来）大于总压力，则包括气相成分各部分压力和方程。如 果气相成分各部分压力和方程不包含在矩阵中，那么未知数 dNg 所有的系数恒等于零。f ptotal 、 f p 和 f p ss 作为优选方程包含在解答器中。所有其它方程作为解答器的相等约束条件。另外，也包括几个不等约束条件： （1）优选方程 f p 余项的值（等于）是非负值，并 作为无机物饱和或非饱和估计只值； （2） ）优选方程 f p ss 余项的值（等于 bpss ???jp ss , jxj ）是非负值，并作为固体溶液组分饱和或非饱和估计值； （3）优选方程 f ptotal 余项的值是非负 值，并作为总气体压力的非负估计值； （4）纯相质量的减少 dnp 小于或等于现存相的总摩尔 数 np； （5）固体溶液某种组分质量的减少 dnp ss 小于或等于存在组分的总摩尔数 np；和（6） 气相摩尔数的减少 dNgas 小于气相摩尔数 Ngas。 液相主要未知数的初始值为溶液组分最初的配比。如果涉及到两种或更多种溶液混合， 初始值为溶液的值加权它们的混合系数的和。 如果交换器或表面已经用一种溶液平衡， 初始 值即以前的平衡值。 如果它们没有用一种溶液平衡， 初始值与用于初始交换组分计算和初始 表面计算值相同。 纯相集合中每个相摩尔数的初始值、 固体溶液集合中每种成分摩尔数的初 始值、气相中每种气体成分摩尔数的初始值恒等于输入值或最后一次模拟中存入的值。 PHREEQC 的 数 据 输 入 ， 批 反 应 和 运 移 模 拟 计 算 依 靠 许 多 数 据 块 。 初 始 条 件 用 SOLUTION 或 SOLUTION-PPREAD 、 EXCHANGE 、 SURFACE 、 GAS-PHASE 、 EQUILIBRIUM CPHASES、SOLID-SOLUTIONS 和 USE 数据块确定。批反应用初始条件和 MIX、KINETICS、REACTION、REACTION-TEMPERATURE 和 USE 数据块确定。运移计 算用 ADVECTION 或 TRANSPORT 数据块确定（见“数据输入说明” ） 。化学动力学数值方法和速率表达式地球化学模型的不足之处是无机物、 有机物和其它反应物经常不能在一次实验或模拟期 间内反应到均衡。根据速率方程，均质或非均衡溶液的动态反应引起水样浓度的变化：dmi ? ci , k Rk ， dt（93）其中 ci , k 是动力反应中物质 i 的化学计量系数，Rk 是物质 k 的总反应速率（mol/kgw/s） 。一 般说来，数据块确定随反应程度变化，导致求解一常微分方程组。 对于不同的温度条件、压力和溶液组分，已经得到了许多动力反应速率。然而，不同研 究者使用不同的速率表达式以适合观察到的速率， 所以很难选择具有足够通用性的速率表达 式（一般很难编进程序） 。在 PHREEQC 版本中， 这个问题用根植于 BASIC 的翻译器得以 解决， 允许以一般的方式在输入文件中定义动力反应的速率表达式。 消除了在程序中硬行编 码速率表达式的需要。数值方法经过一段时间间隔，必须对速率求积分，这涉及到计算溶液中浓度的变化， ，同时，也 解释了其对反应速率的影响。 地球化学动力反应造成难解的方程组， 因为没有及时提取速率， 一些方程的速率（浓度变化对时间的微分）变化巨大，而另一些变化缓慢。 PHREEQC 用 Runge-Kutta（RK）算法解决了这个问题，求速率对时间的微分。Fehlberg（1969）的 RK 方 案，使用了达 6 个中间微分计算值。该方案包括用低阶 RK 方法求得一个误差估计值。用该 误差估计值和用户指定的容许误差比较，自动增加或减小时间间隔使误差在给定的范围内。 用户可在尝试最后的积分间隔前指定估计的中间 RK 的子间隔的数量（见数据输入说明） 。 该方案的系数取自 Cash 和 Karp（1990） 。速率表达式无机物和其它固体动力反应的总速率为：A ? mk Rk ? rk 0 ? V ? ? m0 k? ? ? ， ?n（94）其中 rk 是比速率（mol/m2/s） ，A0 为固体初始表面积（m2） ，V 为溶液的质量（kgw） ， m0 k 为 固体初始摩尔数， mk 为固体在一给定时间的摩尔数， ?mk / m0k ? 是一个系数，用于说明溶n解中 A0/V 的变化和固体的选择溶解和经过时间。对于统一的溶解 spheres andcubes n=2/3。 PHREEQC 中的所有计算用摩尔数，系数 A0/V 必须由用户指定得到合适的比率。 数据库的关键词 RATES 中包括一种物质的比速率表达式 rk。 这些比速率有不同的形式， 主要依赖经验资料的完善程度。资料缺乏时，经常使用的一个简单的速率表达式为：? ? IAP ?? ? ? ? rk ? kk ?1 ? ? ? ? ， ? ? K ? ? k ? ?（95） 其中 Kk 为一经验常数， IAP/ Kk 为饱和度 （SR） 。 这个速率方程可从 transition-state 理论得到， 当形成一活化络合物时（Aagaard 和 Helgeson，1982；Delany 和其它，1986） ，系数ζ 和反 应的化学计量有关。一般ζ =1。这个表达式的一个优点是该速率方程适用于超饱和和不饱 和两种情况，达到均衡时，速率为零。无论何时地球化学反应远离均衡时（IAP/ K&0.1） ，在 很长的时间范围内，速率为常数，并且 IAP/ K 趋近 1.0（均衡）时，速率趋近于零。 速率表达式也可以饱和指数为基础，采用下列形式：? IAP ? rk ? kk ? log? ? K ? ? 。 ? k ?该速率表达式已经成功用于大理石的溶解（Appelo 和其他，1984） 。 速率表达式经常包含与浓度有关的项。Monod 方程是其中之一：（96）? C ? rk ? rmax ? ?K ?C? ?， ? m ?（97）其中 rmax 是最大速率，Km 等于速率为最大值一半时的浓度。Monod 方程一般用于模拟有机 物氧化作用有序的步骤（Van Cappellen and Wang，1996） 。一系列速率表达式可由相应氧化 物的能量生成得到；首先消耗 O2，然后是 NO3-，接下来是其它，更慢反应的氧化物如 Fe2O3 和 SO4-2。Monod 方程的系数可由各个阶段的一级反应速率方程得到。土壤有机物降解的一 级反应速率方程为：dsC ? ? k1sC ， dt（98）其中 sC 是有机碳含量（mol/kg soil） ，K1 是一级衰变常数（s-1） 。在含氧环境下 K1 的值约等 于 0.025yr-1， 然而在荷兰砂质含水层中，NO3-是氧化物，k1≈5e-4 yr-1。在地下水甚至在有 机物氧化降解的氧化还原作用范围之外发现浓度高达 3μ Μ O2，而且可作为（依赖浓度的） 需氧降解的速率等于脱氮反应速率时的浓度。在 Monod 方程中利用系数 rmax= 1.57 e-9s-1 和 Km=294μ M 得到一级衰变反应（对于 0.3μ mO2，k1=0.025yr-1 和 3μ Μ O2，k1=5e-4 yr-1） ， 并且氧作为限制溶质。脱氮的类似估算是以 NO3-=3mM，k1=5e-4 yr-1 和 NO3-=3μ M，k1=1e-5 yr-为基础的，生成 rmax= 1.67 e-11s-1 和 Km=294μ M。在淡水含水层中，有机碳降解的综合 Monod 表达式为：?s RC ? 6sC ? C ? sC ? 01.67 *10?11 mNO ? ? ?? 1.57 *10?9 mO2 ? ? 3 ?? ? ? 4 ? 4 ?? 2.94 *10 ? mO 1.55*10 ? m ? ? NO3 ? 2 ?? ?（99）其中系数 6 是把土壤中的 sC 换算为空隙水中的 sC 时得到的。 有机物分解的另一重要方面是一部分很难溶而且特别不宜降解。已经提出一些模型解 释部分有机碳仍存在下去的趋势。 临时假定一系数 ?? sC ? sC ? 0? ? 组成总的速率二级反应。 该系数暗 ? ?示着浓度降为原始浓度的 1/10 将造成速率进一步降低为原来的/100。必须注意该系数用于 逼近一个很复杂的过程，并且需要严密分析，在 PHREEQC 版本中也可能给出确定速率的 其它灵活方法。 也有速率表达式是根据详细测量溶液中影响速率的浓度变化得出的。 例如 Willamson 和 Rimsidt（1994）给出的黄铁矿的速率表达式是：0.5 ?0.11 rpyrite ? 10?10.19 mO mH ? ， 2（100）式中含氧克分子浓度的平方根，并且说明了随 PH 值增大，速率有微小增加。该速率表达式 只用于溶解反应，且只用于溶液含氧时。当溶液接近均衡或氧被耗尽是，该速率表达式是不 适合的。 既适用于溶解又适用于沉淀的一个更完整的速率表达式是钙的速率方程：? ， rcalcite ? k1 H ? ? k2 ?CO2 ? ? k3 ?H2O? ? k4 Ca2? HCO3? ?????（101）其中括号内的化学符号表明为活度，系数 k1、k2 和 k3 已被 Plummer 等人（1978）在含带电 CO2 的钙溶解实验中确定为温度的函数。速率方程包括前半部分 rf（方程 101 的前三项）和 后半部分 rb（方程 101 的最后一项） ，并且既适用于溶解又适用于沉淀反应。后一部分包括 一个系数 k4，起值依赖于溶解组分。在单纯的方解石系统中，碳酸氢根离子浓度大约为钙离 子浓度的两倍，且后一部分速率表达式约等于：rb ? k4 Ca 2 ? HCO3? ? k4 2 Ca 2 ? 。均衡时，[Ca2+]为饱和时的活度[Ca2+]s。rcalcite=0，则??????2（102）2k4 ??Ca ?rf2? 2 s。（103）结合方程 101、102 和 103 得：2 ? ? 2? ? ? Ca ? ?。 rcalcite ? rf ?1 ? ? 2? 2 ? ? ? ? Ca s ? ? ? ?? ?? ?（104）CO2 压力为常数的纯系统 Ca-CO2 中，离子的活度积（IAP）为： IAPcalcite?Ca ??HCO ? ?2?? 2 3PCO2?Ca ? 和 K ?42? 3PCO 2calcite?Ca ? 。 ?42? 3 sPCO2（105）这样，对于一个方解石系统，方解石的速率近似为：2 ? ? ? IAP ? 3 ? ? ? rcalcite ? rf 1 ? ? ? ? Kcalcite ? ? ? ? ? ? ? ?（106）其中 rf 为方程 101 给出的前三项。数值方法和运移模拟方程PHREEQC 可 以 模 拟 几 个 一 维 运 移 过 程 ， 包 括 ：（ 1 ） 弥 散 ，（ 2 ） 平流（对流） ， （3）对流和弥散，和（4）对流和进入滞流区的弥散，这作为双空隙处理。所 有这些过程可和均衡及动力反应结合起来。平流-反应-弥散方程迁移化学物质（fig.1）的质量守恒生成对流-反应-弥散（ARD）方程：?C ?C ? 2C ?q ? ?v ? DL 2 ? ， ?t ?x ?x ?t（107）其中 C 是水中的浓度（mol/kgw） ，t 是时间（s） ，v 是水流速度（m/s） ，x 是距离（m） ，D L 是水力渗透系数[m2/s，DL = De+α Lv，De 是有效渗透系数，α L 是渗透率差（m）]，q 是固 相中的浓度（空隙中用 mol/kgw 表示） 。? v?C ? 2C 代表对流性运移， DL 代表弥散运移， ?x ?x 2?q 是由于反应在固相中浓度的变化（q 的单位同 C） 。一般假定 v 和 DL 等于所有的溶解物 ?t质，所以 C 可以是一种元素的总溶解浓度，包括所有的氧化还原物质。 图1对流-反应-弥散（ARD）方程中各项方程 107 的运移部分用先进行的、弥散中心的和对流运移迎风向的显有限差分法解答。 每种元素的化学交互作用项按每个时间段分别由运移部分计算， 并且是所有均衡和非均衡反 应速率的和。在 split-operator scheme（Press 和其他，1992；Yanenko，1971）中数值法随 ARD 方程中各项。对于每个时间步长，先计算对流运移，其次计算所有的均衡和动力控制 的化学反应， 再次计算弥散运移以及其所有的均衡和动力控制的化学反应。 这种方法不同于 其它多数的水文化学模型（Yeh 和 Tripathi，1989） ，因为动力和均衡化学反应是在对流计算 和弥散计算之后进行的。这减少了数值分散而且不必在物质组成和运移之间的重复计算。 split-operator scheme 的优点是通过对方程的各个部分的时间步长的调整，使其适合于 grid size（栅格尺寸） 。通过下列关系减小时间和距离离散之间的数值分散：??t ?A ? ?x ，v其中（108）是对流运移的时间步长（s） ，△x 是单元长度。弥散计算中数值的不稳定（波动）用下式估计：??t ?D ? ??x ?23DL，（109）其中 ??t ?D 是对流运移的时间步长（s） 。方程 108 和 109 的两个条件分别是对流运移计算的 Courant 和弥散运移计算的 Von Neumann 判据（例如，Press 和其他，1992） 。当 ?x ? ? L 时， 数值分散可忽略不计， 因为物质的对流运移和弥散运移一样重要甚至更重要。 当用一小栅格 间距减小数值分散时，对流运移计算的时间步长（方程 108）可能小于弥散运移计算的时间 步长（方程 108） ，因为前者包含 fine grid 的平方项用多个弥散时间步长解决这个矛盾，如? ??t ? ? ??t ?DA，并且在每个弥散时间步长计算化学反应。PHREEQC 的数据输入，一个弥散时间步长必须等于对流时间步长 ??t ?A ，或模拟弥散时等于弥散时间。此外，必须确定 变换次数，即计算的对流时间步长（弥散时间）的个数。 Central difference scheme 中的弥散运移必须是单元的混合。混合系数 mixf 定义为：m ixf ?DL ??t ?A ， 2 n??x ?（110）其中 n 为一正整数。限制条件为该单元和后面的不再混合，也就是说，mixf 必须小于 1/3， 如方程 109。当方程 110 中 n=1 时，mixf 大于 1/3，增加 n 的值使 mixf 小于或等于 1/3。弥 散时间为 ??t ?D ???t ?A ，并混合 n 次。n通过与简单线性交换的解析解相比较检查数值方法。 当交换不完全反应的交换系数等于 两个同价阳离子时形成线性交换结果。它给出了一个线性延迟值 R=1+CEC/C，其中 CEC 是 阳离子交换能力，用 mol/kgw 表示。在下列例子中，一个长 130m 的流水管中水的初始浓度 为 C（x，0）=Ci=0。排出的水浓度为 C=C0=1 mol/kgw，且空隙水流速为 v=15m/year。分散 率差为α L=5， 有效弥散系数为 De=0m2/s。 4 年后可给出两种化学成分的分布图， 一个为 R=1 （CL-） ，另一个为 R=2.5（Na+） 。 该问题需考虑两个边界条件。一个需要 x=0 时的流量或第三类边界条件：C ?0, t ? ? C0 ?DL ?C ? xend , t ? 。 v ?x（111）这个边界条件适用于输入管的横段面与实验管的横断面相比很小时。 ARD 方程的解 （Lindstrom 和其它，1967） ：C ? x, t ? ? Ci ?其中，DL=αL1 ?C0 ? Ci ?A 2v：（112）? x ? vt / R ? ? ?x ? vt / R ?2 ? x ?? A ? erfc? exp?? ?? ? 4? vt / R ? ?? 4 ? vt / R L L L ? ? ? ? ? x ? vt / R ? ? x ? 1? x vt / R ? ? ? ? ? ? ? 1 ? ? exp erfc ? ?? ? ? ? 2? ? ? 4 ? vt / R L L ? ? ? L? L ? ?。（113）插图 2 给出了具有不同（grid）栅格间距，△x=15、5 和 1.67m 时的三种模拟，相应的 和 和 和间距（m）图 2.具有离子交换和流通边界的一维运移解析解与 PHREEQC 计算在不同栅格间距时的比较时间分别为 ??t ?A =1、1/3 和 1/9 年。对于 CL- ，R=1，三种模拟的前部是不可区别的且与 解吸解完全一致。延迟离子 Na+ ，R=2.5，对于所有的（grid）栅格间距，漏过曲线的平均 位置是周围能正确的且与解析解一致。 最小的栅格间距结果与解析解最相近， 但需要最多的 计算时间。 计算时间主要依赖于 PHREEQC 中调用地球化学子程序的次数，没有动力反应时，调 用次数与单元个数*对流间隔个数*（1+弥散间隔个数）成比例。在这个例子中 DL= De+αLv=0+5*15m2/yr。这样，根据方程 110，对于逐渐减小的单元尺寸，mixf 分别等于 1/3、1 和 3。 对于 15m 的单元（mixf=1/3） ，对流步长取为弥散步长；对于 5m 的单元（mixf=1） ，三个对 流步长取为弥散步长；对于 1.67m 的单元（mixf=3） ，九个对流步长取为弥散步长。插图 2 说明了 CL（C/C0=0.5）四年旅行后对流的锋面曲线，当它达到 60m 时；对于 15m 的单元， 这需要 4 个对流步长。流管由 9 个单元组成，需对每个步长进行地球化学计算；因此，反应 计算次数为 9*4*（1+1）=72。大量的单元和对流步长适合于小 grid。另外两种情况反应计 算的调用次数为 27*12*（1+3）=1296；和 81*36*（1+9）=29160。 这儿给出的例子具有线性延迟使得能和解析解比较。然而线性延迟受到很大的数值分 散，而且在一定意义上，这个例子是数值分散糟糕的典型。在地球化学行业的许多方面，化 学反应帮助抵消数值分散，因为反应趋向于使波面曲度增大（sharpen fronts） ，如沉淀溶解 反应和沉降分离。 在其它例子中， 和不活泼的离子交换可能得到不受数值分散影响的真正的 物质弥散。 另一个边界条件是 Dirichlet 或者说是第一类边界条件， 它包括 x=0 时的常量浓度 C （0， t） ：C?0, t ? ? C0（114）当水从与下面土壤充分联系的蓄水池中渗漏时， 该边界掉件有效， 例如水从池塘或海水侵入 下覆沉积层。这种情况下 ARD 方程的解（Lapidus 和 Amundson，1952）为：C ? x, t ? ? Ci ?其中，1 ?C0 ? Ci ?B 2（115）? x ? vt / R ? ? ? ? ? ? ? exp? x ?erfc? x ? vt / R ? 。 A ? erfc? ? ? ? 4? vt / R ? ? 4? vt / R ? ? ?L ? L L ? ? ? ?（116）插图 3 说明了和以前运移例子相同的离散化模拟结果。 用三种栅格尺寸精确模拟守恒的 溶质（CL-，R=1） 。延迟物质（Na+，R=2.5）用大的网格说明了数值分散，但是，总的波面 是一致的。 在常浓度边界条件下， 因为指定的 x=0 处的条件弥散时间步长个数是流动时的两 倍。 第一类边界条件通过柱表面和外界溶液的接触会增大弥散。 通过边界的物质流动相应增 加，波锋比插图 2 中增加了几米高。例 11（临界点在柱的输出口）和例 12（弥散从恒定源 进行）给出了更多的分析溶液比较。 距离（m） 图 3. 具有离子交换和恒定边界条件的一维运移解析解与 PHREEQC 计算在不同栅格间距时的比较。热量传递由热量守恒的出热量传递方程， 或者说是温度变化。 一种物质的热量传递方程和其对流 -反应-弥散（ARD）方程是一致的：??wkw?T ?T ? 2T ? ?1 ? ???s k s ? ????wkw ?v ? ?L 2 ， ?t ?x ?x（117）其中 T 是温度（℃） ， ? 是空隙度（与总体积的比值，无单位） ，ρ 是密度（kg/m3） ，k 是比 热（kJ℃kg-1） ，κ 包含对流弥散和含水层的热传导（kJ℃m-1kg-1） ，下标 w 和 s 分别代表水 和固体。对于水和固体，温度 T 假定为一致的。 方程 117 两边同除以 ??wkw ：RT?T ?T ? 2T ? ?v ? ?L 2 ， ?t ?x ?x（118） 其中 RT ? 1 ??1 ? ???s ks 是温度延迟（无单位）系数， ???wkwL?? 热量耗散系数。该系数 ??wkw包括单纯弥散部分和由于对流耗散的部分： ?L ? ?e ? ?Lv ，与水力弥散系数相似。类似性}
(window.slotbydup = window.slotbydup || []).push({
id: '4540180',
container: s,
size: '250,200',
display: 'inlay-fix'
添加成功至
资料评价：
所需积分：5第二讲PHREEQC用户指导_甜梦文库
第二讲PHREEQC用户指导
摘要PHREEQC 第二版是一个用 C 语言编写的计算机程序， 对各种各样低温下的地球化学性 质进行了演算。PHREEQC 是以离子联系的水化学模型为基础的，可以推算（1）生成物和 饱和系数； （2）涉及到可逆反应以及不可逆反应的批反应和一维（1D）运移计算，可逆反 应包括水、矿物/无机溶液、气体、固体溶液、表面络合、离子交换平衡；涉及到的不可逆 反应包括指定成分摩尔转换、动态控制反应、溶液混合和温度变化； （3）逆向模拟实验，其 中多组的无机物和气体摩尔转换以解释在特定成分不确定范围水体之间混合物的不同。 和第一版相比， PHREEQC 第二版的新特点如下：具有在一维运移计算中模拟弥散（或 扩散）和滞流区的能力，用用户确定的速率表达式模拟分子反应，模拟标准的多种成分或非 标准的两种成分的固体溶液的沉淀和溶解， 模拟定体积气相和定压力气相， 考虑表面系数或 交换位置随着无机物的溶解和沉淀或者分子反应的变化而变化， 自动采用多套收敛参数， 打 印用户指定量到原始输出文件和（或）适合输入到扩展表格的文件上，以一种与扩展表程序 更兼容的形式确定溶解成分。 这个版本报告中说明了化学平衡、动态平衡、运移计算以及逆向模拟计算基础方程式， 描述了程序输入，以及举例说明了许多程序功能目的和范围这个报告的目的是说明 PHREEQC 程序的理论和操作，包括组成成分方程的确定，转 换这些方程为数值计算方法的解释，补充数值方法计算机代码组织的描述，程序输入描述， 以及一系列数据组输入和许多论证的程序功能的数据输入和模拟结果的说明。生成方程和正向模拟在报告的这一部分，说明了用于确定水样热动力活动，离子交换物质、表面络合物 质、气相成分、固体溶液和纯相的代数方程式。首先，说明了水样、交换种类和表面性质的 热动力活动和质量作用方程。然后，定义一组函数，用 f 表示，他们必须能同时求解以确 定给定条件下的平衡， 许多这样的方程是各种元素或元素的价电子状态、 交换位置和表面位 置的摩尔平衡方程， 或来自于纯相和固体溶液的质量作用方程。 附加函数用于确定水的碱度、 活度、电荷平衡、气相均衡、离子强度和表面络合均衡。每个函数包含最少量的变量，使函 数的个数等于变量的个数。 程序用一种修改过的牛顿? 罗弗逊方法解答了伴随的非线性方程。 这种方法利用了函数的偏差和每个函数关于主要未知数组或主要未知数的偏导数。 为了清晰 起见，微分计算中用的变量组指主要未知数，每个函数的全导数不用求导说明。在下列方程 中，没有下标或下标“ （aq） ”指液相中的实体， “（e） ”指交换， “（g） ”指气体， “（s） ” 指表面， “（ss） “指固体溶液， “（p） ” 指相的状态。活度方程和质量作用方程在这一部分，定义了水的活度、交换过程和表层物质，并描述了每种物质的质量作用联 系。 质量作用方程是化学系统中根据主要未知数每种物质的摩尔数的质量作用表达式。 然后， 求这些方程关于主要未知数的差分。 每种物质的摩尔数方程和偏导数方程将被组成的摩尔平 衡，电荷平衡和相均衡函数取代。水样PHREEQC 允许单个液相的物质形成或物质平衡。然而多个液相可在一个运行过程中 定义， 一个液相可被定义为一种或多个液相的混合物 （见D数据输入说明‖中的 MIX 关键词） 。 一个例外是液相中溶解物质假定处于热动力平衡状态。 在最初的溶解计算中， 允许氧化还原 元素的化学价不平衡。每个水样 i 的未知数包括活度 ai，活度系数 ri ，重量摩尔浓度 mi， 溶液中摩尔数 ni 。 PHREEQC 根据主要物质改写了所有的化学方程式。其中一个与每个元素（例如，钙离 子 Ca+2）或元素化学价态（例如，三价铁离子 Fe+3）有关的主要水样，加上氢离子活度、水 中电子活度和水的活度。一些程序，例如 MINTEQA2（Allision 和其他，1990）和 MINEQL+ （Schecher 和 Mcavoy，1991） ，把这些物质用术语“成分”表示，但是，这个术语在这儿不 能用，因为易与‘Gibbs’的相规则中的“成分”混肴。PHREEQC 中，每个水样主要物质 用 SOLUTION-MASTER-SPECIES 数据块确定（见“数据输入说明” ） 。数值方法将未知数 的个数减少为主要未知数的最小个数， 并且反复迭代这些未知数的值， 直到得到数学方程组 的解。水溶液的主要未知数是主要物质活度的自然对数、水的活度的自然对数 ? H2O、离子 强度μ 、水溶液中有溶解能力的水的质量 Waq。 下列关系适用于所有水样（除了水成电子和水本身） ：ai=ri mi 和 ni=miwaq. 离子缔合模 型中水样之间的均衡要求满足所有的水样质量方程。例如，水样 CaSO40 的缔合化学反应是 Ca2++SO42-=CaSO40 。 这 个 反 应 在 25 ° C 时 的 logK 是 2.3 ， 形 成 质 量 作 用 方 程 ：10 ?2.3CaSO 0 4 ? 2? ? 2? Ca SO 4?（1）一般情况下，质量作用方程式可被写为：M aq mKi ? ?i? ?m m ,i ,?c（2）其中 Ki 是温度均衡常数， cm，i 是物质 i 中主要种类 m 的化学计量系数，Maq 是水中主要物 质总数，cm，i 的值可正可负。在 PHREEQC 中，缔合化学反应右边各项系数为负，左边各项 系数为正。主要物质采用同样的形式，质量作用方程可简化为 1=am/am。 一种水样的总摩尔数 i 可由质量作用表达式求得：M aq mni ? miWaq ? KiWaq???ic m ,i m.（ 3）牛顿.罗弗逊方法用摩尔数关于主要未知数的全微分。计算如下：M aq mdni ? ni [d ln(Waq ) ? ? cm,i d ln(am ) ?? ln(? i )d? . ??（4）水样的活度系数由戴维斯方程确:? ? ? ? log ? i ? ? Azi2 ? ? 1 ? ? ? 0.3? ? , ? ?或者用扩展的或 Debye-Huckel 方程确定：（ 5） log ? i ? ?Azi2 ? ? bi? , 1 ? B?i0 ?（6）其中 zi 是水样 i 的离子电荷数，A 和 B 是关于温度的常数 。假如 bi 为零，方程 6 是扩展的 Debye-Huckel 方程；假如 bi 不为零，方程 6 是 Debye-Huckel 方程（见 Truesdell 和 Jones， 1974） 。在扩展的 Debye-Huckel 方程中，ai0 是离子大小参数，反之在 WATEQ Debye-Huckel 方程中，ai 和 bi 是适合于 mean-salt 的离子特定参数。否则，如果没有在数据库文件或输入 数据组中说明，Davies 方程用于带电物质。对于不带电物质，活度系数方程中第一项为零， WATEQ Debye-Huckel 方程简化为 Setchenow 方程 （lnri=biμ ） （见 Langmuir， 讨论版 1997） 。 否则，如果没有指定，对所有不带电物质，bi 假定为 0.1。 Davies 方程中，这些活度系数方程关于离子强度的偏微分为：? ? ?? ? 1 ?? ln ? i ? ? ln(10) ? Az i2 ? ? 0 . 3 2 ? ?? , ?? ? 2 ? ? ? 1 ? ?? ???（7）对于扩展的 WATEQ 或 Debye-Huckel 方程为：? ? ? Azi2 ? ln ? i ? ? ln(10) ? bi ? , 2 ? 2 ? B?0 ? ? 1 ? ?? i ? ???（8）PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的温度要 求，每种水样的活度系数参数通过 SOLUTION-SPECES 数据块确定。元素的主要种类和元 素 的 化 合 价 态 用 SOLUTION-MASTER-SPECES 数 据 块 确 定 。 一 种 溶 液 的 组 成 用 SOLUTION 或 SOLUTION-SPREAD 数据块确定（见“数据输入说明” ） 。交换种类离子交换平衡通过交换位置的非齐次质量作用方程和摩尔平衡方程包括在模型中。 PHREEQC 允许多种交换与液相在均衡中共存，用“交换集合”表示。这种方法使用了质量作 用表达式，它以每个交换器中水样和虚构的空的交换位置之间的不完全反应为基础 （Appeio 和 Postma 1993） 。这种空的交换点是交换器的重要种类，其活度的对数是一个附加的主要未 知数。它由 EXCHANGE-MASTER-SPECES 数据块确定（见“数据输入说明” ） 。然而，交 换器的主要物质不包括在摩尔平衡方程中， 使其物质浓度为零。 它的活度在物质上是无意义 的，但是如此却使所有的交换位置充满了其他交换物质。 交换计算的未知数是活度 aie，在 PHREEQC 中定义为当量比值乘以活度系数 ? ie 和交 换器 e 中每种交换物质 ie 的摩尔数 nie 。当量比值是被一种交换物质占据的交换点的摩尔数 比交换点的总数。一种交换物质的活度为 ?ie ? ? iebe,ie nie Te，其中 be，ie 是被交换物质 ie 占据的交换器 e 的当量数， Te 是交换器中交换地点的当量总数。 Te 作为系统中交换器交换地点的 当量总数，不必与每克水的当量数相等，因为系统中水的质量可能多或少于 1kg。根据默认 值，交换物质的活度系数为 1.0，但是可以根据水的离子强度和被交换物质占据的交换地点 当量数选择使用 Davies，扩展的 Debye-Huckel 或 WATEQ Debye-Huckel 活度系数。 水和交换物质之间的平衡要求满足所有交换物质的质量作用方程。交换物质 CaX2 的 缔合反应是 Ca2++2X-= CaX2，其中 X 是默认数据库的交换主要物质。活度的当量比值和这 种化学反应形式的使用统称为 Gaines-Thomas 协议（Gaines 和 Thomas，1953） ，并且在 PHREEQC 版本中用于数据库 phreeqc.dat 和 wateq4f.dat 文件。[在 PHREEQC 版本中，也可 能用 Gapon 协议，也使用当量比值，但是将交换反应写为 0.5 Ca2++X-= Ca0.5X。见 Appeio 和 Postma （1993）再次讨论。]在默认数据库文件中，钙离子交换的 logK 伪.8，形成质量 作用方程：100.8 ??CaX 2 ?Ca 2? ?2 X?,（9）一般情况写，质量作用方程可写为：MKie ? ?ie ? ? m m ,ie?c m,（10）其中 m 随所有的主要物质变化，包括交换的主要物质，cm，ie 是在物质 ie 不完全缔合反应中 主要物质 m 的化学计量系数， kie 是不完全反应选择常数。cm，ie 的值可正可负。PHREEQC 中，缔合反应右半边各项系数为负，左半边各项系数为正。 对于交换物质，物质 ie 总摩尔数方程为：nie ? Kie??mM? c m ,ie m? be,ie ? ie ? ? T ? e? ? ? ?.（11）在数值方法中， 交换器主要物质活度的自然对数是主要未知数。 物质 ie 摩尔数关于主要未知 数的全微分是；?M ? ? dnie ? nie ? ? ? cm,ie d ln(? m ) ? ?? ln(? i )d? ? ?. ? m ?（12）PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的温度要 求以及每种交换物质的活度系数表达式通过 EXCHANGE-SPECES 数据块确定。 交换的主要 种 类 用 EXCHANGE-MASTER-SPECES 数 据 块 确 定 。 交 换 位 置 数 量 和 交 换 器 组 成 用 EXCHANGE 数据块确定（见“数据输入说明” ） 。表面种类表面络合过程通过非齐次质量作用方程、 表面位置的摩尔平衡方程和每个表面的电荷电势联系包括在模型中。 PHREEQC 允许多种表面和表面位置类型液相在均衡中共存， 用 “表 面集合”表示。PHREEQC 表面种类质量作用方程有两个公式： （1）一个包括静电电势项； （2）另一个包括所有的静电电势项。如果使用模型静电电势项，由于存在表面电荷和表面 静电电势，就需用到附加方程和质量作用术语。 交换反应和表面反应计算之间的两个根本不同点在于，交换反应表示为不完全反应， 使主要物质不能出现在任何摩尔平衡方程中， 并且交换物质被期望是中性的 （或不带电的） 。 表面反应不是不完全反应， 所以主要物质是物质上真正存在的物质， 并且可出现于摩尔平衡 方程中，表面物质可带正电，也可带负电，或不带电。 在 Dzombak 和 Morel（1990）中，表面络合反应的基本理论包括静电电势。该理论假 定活性地点个数为 Ts（eq） ，比面积 AS（m2/g） ，表面质量 Ss（g）已知。两个附加的主要未 知数是（1）参量 ln ??s ? ln? e 2 RT ? ?? ? ?F?s? ? ?F?s ，其中 F 是 Faraday 常数（96493.5JV-1eq-1） ， 2 RTΨ s 是表面电势（Volts） ，R 是气体常数（8.3/4）Jmol-1k-1 ） ，T 是温度（kelvin） 。 （2）主要 表面物质活度的自然对数。参量 lnaΨ s 右侧项在术语标准中用 a2 定义。这是用于 Dzombak 和 Morel（1990）中一个不同的未知数，但是由于所有的方程书写和这个主要未知数一致， 所以结果相同。 表面物质活度假定等于给出非空的表面位置类型的摩尔比值。换句话说，一种表面物 质完全覆盖了一种给出的表面位置时是处于标准状态 （活度为 1） 。 这个协议不同于 Dzombak 和 Morel（1990） ，它假定一个表面物质（固相概念上）的活度数值上等于体积克分子浓度 （溶液浓度） 。假如只考虑 monodeneate complexes（如 Dzombak 和 Morel（1990 所做） ，消 去质量作用方程中各项，不考虑标准状态协议，可得到一致的数值结果。然而，体积克分子 浓度 （容摸） 协议用于 multideneate complexes （bidentate， tridentite， 其他， cf.Appelo and postma， 1999)，表面位置浓度存在一个明显的不同。假如一个容器盛有表面包含 multideneate species 且处于均衡的溶液，并且增加了更多同样的溶液，体积克分子浓度不变，溶液的成分和表面 不断改变。这种情况体积克分子浓度明显不正确。 “Hfo‖(含水 Fe2O3) 在默认数据库中用作D_w‖， 指明了较低的亲和力或很弱的位置， D_s‖指明强的亲和力或强的位置。 “Hfo_wOH‖用于说明D在软弱位置‖一种中性表面物质，带 负电软弱位置形成物的缔合反应（在一定意义上，在缔合反应中定义的种类在方程的右侧） 可以被写为： Hfo_wOH→Hfo_Wo-+H+. 质量作用表达式，包括静电电势项项，是 （13）Kint Hfo_wO -?? Hfo_wO - ? H ? ? Hfo_wOHe?F?s RT,（14）其中是反应原有的平衡常数， e?F? s RT是一个因子，说明了将一种带电物质（H+）从带电表面移出。一般，表面物质 i? s k ? 的质量作用方程为： M c m ,i? s ? ? F? ? s k Kiint ? ? ? ? i? sk ? ? ?m ? ?e RT ?zi? sk ? ? sk ? m ? ?（15）其中 K iint 是原有的平衡常数，是表面 s 中表层位置类型 k（软弱的或坚固的 Dzombak 和?sk ?Morel，1990 中）的第 i 层物质；m 随所有主要物质 M 变化，包括表层主要物质； cm , i? s ? 表k层物质 i? s k ? 缔合反应主要物质的化学计量系数， ?zi? s ? 是形成表层物质时表面电荷的净变k化。 cm , i? s ? 可正可负。PHREEQC 中，缔合反应中右边半边各项系数为负，左半边各项系数k为正。 一种表层物质 i(sk)的总摩尔数方程为：ni? s ? ? ? ikTs k bi? s ?k? sk ?? K i Ts k e? sk ? ?? F?s ? ? ? RT ?z i? sk ? ?? M ? ? ???mc m ,i? s ? k m? KiTs k bi? sk ?? sk ?? 2 ?z i? s ? k ?s??mM,（16）c m ,i? s ? k m其中 Tsk 是表面位置类型总数，bi（sk）是与物质相联系表面位置的数量。物质 i(sk)的摩尔数关 于主要未知数的全微分为：?M ? dni? sk ? ? ni? sk ? ?? cm,i? sk ? d ln ?m ? 2?zi? sk ? d ln ??s ? . ?M ?（17）表面物质质量作用方程的第二个公式在质量作用表达式中不包括静电电势项（ -no-edl identifier 在 SURFACE 数据块中） 。一种表层物质摩尔数方程同方程 16，除了涉及α 式没有显示。否则，摩尔数的全微分同方程 17，除了没有最后一项。 PHREEQC 的数据输入，摩尔平衡化学方程式和质量作用表达式，logK 和它的表面物 质的温度要求通过 SURFACE-SPECES 数据块确定。表层主要物质或表层位置类型通过 SURFACE-MASTER-SPECES 数据块确定。表层的一致性，每个位置类型的当量数，表面组 成，比表面积和表层质量用 SURFACE 数据块确定（见“数据输入说明” ） 。Ψs的因 气相成分 多成分气相和液相之间的均衡可用非齐次质量作用方程和一个总压力方程（仅为故乡 压力气相）模拟。只有一个气相可与液相在均衡中存在，但气相可包括多种成分。所有的气 体成分假定为理想气体，气相假定为理想的气体混合物。 如果气相体积固定，那么该体积内的压力将随着反应程度变化，但是气相中每种成分 总是存在的。对于一个固定体积的气相，不需要另外的主要未知数，气体中一种成分的摩尔 数可以从液体主要物质的活度中推算出来。 如果气相压力固定，那么气相体积将随着反应程度变化，假如组成成分的部分压力之 和小于指定的总压力，固定压力气相将不存在，并且气相中任一气体成分将不能存在。固定 压力气相中，气体成分的总摩尔数方程中包括一个附加的主要未知数 Ngas。 根据理想假定，一种气体成分的逸度（或活度）等于其部分压力。PHREEQC 使用溶 解方程式，在一定程度上，气体成分假定处于化学反应的左边。CO2 溶解反应可写为 CO2（g）= CO2（aq） 。 （18）亨利定律常数将气体成分的部分压力 （理想气体数值上等于逸度） 与液体物质的活度联系了 起来。对于二氧化碳，亨利定律常数为 10-1.468[假定为理想气体，单位为大气压（atm）]， 下列质量作用方程适用于均衡状态：PCO2 ? 101.468 ?CO2 ? a q? ,（19）其中， P 。一般情况下，气体成分的部分压力可 CO2 是用液相的活度计算的部分压力（atm） 根据液相活度写为：1 aq cm , g Pg ? ? ?m , Kg mM（20）其中，Pg 是气体成分 g 的部分压力，用液相的活度计算；Kg 是气体成分的亨利定律常数； cm，g 是溶解方程中水成主要 m 的化学计量系数。cm，g 值可正可负。PHREEQC 中，溶解反应 方程中右半边各项系数为负，左半边各项系数为正。 对于定体积气相，气相的总体积指定为 Vtotal，但压力是变化的。在平衡状态，气体成 分的摩尔数如下计算：ng ?Vtotal pg RT?Vtotal RTKgM aq m??cm ,g m.（21）气相中该成分的全微分为dng ?Vtotal aq ? ng cm,g d ln ?m RT mM（22）对于定压力气相，气相的总压力指定为 Ptotal，但体积是变化的。在平衡状态，气体成 分的摩尔数等于气体总压力的比值乘以气体的总摩尔数，如下计算：ng ? N gasPg Ptotal?N gas Ptotal K g??mM aqcm ,g m.（23）气相中该成分摩尔数的全微分为dng ?Pg PtotaldNgas ? ? ng cm,g d ln ? m .mM aq（24）PHREEQC 的数据输入、质量作用方程、亨利定律常数以及该常数的温度要求通过 PHASE 数据块确定。气相类型（定体积或定压力） 、气相计算中包括的成分和初始气象成分 通过 CAS_PHASE 数据块确定（见“数据输入说明” ） 。牛顿-罗弗逊方程用一组函数（用 f 表示）用于说明复杂的均衡，这些方程主要来源于用摩尔平衡方程 和电荷平衡方程代换物质摩尔数方程（来源于前一部分的质量作用方程） 。达到平衡时，有 关待定平衡计算的所有函数值为零。函数的零点通过 Tventm-Raohsm 方法获得。据此每个 函数关于每个主要未知数微分形成 Tacobin 矩阵。由 Tacobin 矩阵得到一组线性方程组,其 解可近似为非线性方程的解。 通过反复解答连续的线性方程组， 可得到非线性方程的一个解。 在这一部分说明了每个用于数值方法的函数 f 及其关于形成 Jacobian 矩阵的主要未知数的全 导数水的活度水的活度用根据 Raoult 定律(Gaoult’s and Christ,1965,p.65-66)得出的近似值计算：N aq i? H 2 O ? 1 ? 0.017 ?ni . Waq（25）函数 fH2O 定义为：f H 2O ? Waq ? H 2O ? 1 ? 0.017? ni ,i??N aq(26)该函数的全微分为：dfH 2 O ? Waq ? H 2O d ln ? H 2O ? ? H 2O ? 1 Waq d ln?Waq ? ? 0.017? dni .i?? ??N aq（27）主要未知数是水的活度的自然对数 ln ? H 2 O 。离子强度水溶液的离子强度是一个主要未知数，定义如下：??1 aq 2 ni ? zi W , 2 i aqN（28）函数 f? 定义为：1 aq f? ? Waq? ? ? zi2 ni , 2 i该函数的全微分为：N（29） df? ? ?Waq d ln?Waq ? ? Waq d? ?1 aq 2 ? zi dni . 2 iN（30）定体积多项成分气相平衡对于定体积气相， 每种气体成分的摩尔数可由水中主要物质的活度计算得到， 数值模型 处理气相成分的方法与它处理水中物质的方法不同。 每种气体成份的摩尔数 ng 用在摩尔平 衡方程中用成分形式表示，dng 在 Jacobian 矩阵中用摩尔平衡方程表示。定体积气相均衡计 算不需要附加方程 f。 PHREEQC 数据输入，质量作用方程，Henry 定律常数，以及其温度要求用 PHASES 数据块确定，气相类型（定体积或定压力） ，气相计算中包括的成分和初始气相成分由 GAS_PHASE 数据块确定（见“数据输入说明” ） 。定压力多成份气相平衡定体积气相中每种气体成份的摩尔数由水中主要物质的活度和气相中气体成份总摩尔 数 Ng 决定。每种气体成份的摩尔数 ng 在摩尔平衡方程中用（elements）成份表示，dng 在 Jacobian 矩阵中表示摩尔平衡方程。 定压力多成份气相和液相平衡需要一新的方程――气体 成份部分压力之和等于总压力 Ptotal。函数 f Ptotal 定义为：Ng gf Ptotal ? Ptotal ? ? Pg ,其中，Ng 是气相中气体成份的总摩尔数。（31）f Ptotal 关于主要未知数的全微分（其中正的 dNgas 表示溶液浓度增加）为：dfPtotal ? ??? cm,g Pg d ln ? m .g m N g M aq（32）PHREEQC 的数据输入，质量作用方程，Henry 定律常数以及该常数的要求用数据 块 PHASES 确定。气相类型（定体积或定压力） ，气相计算中包括的成份和初始气相成份用 数据块 GAS-PHASES 数据块定义（见“数据输入说明” ）纯相均衡气相和液相之间的平衡， 包括具有固定部分压力的气体， 通过非齐次质量作用方程包括 在模型中。PHREEQC 允许多种纯相（用纯相集合表示）和液相在均衡中存在，受 Gibbs 纯 相规则的限制，纯相的活度假定一致为 1.0。每种纯相的附加未知数为系统中存在的纯相的 摩尔数 np ,其中 p 指第 p 个纯相，np 表示每个纯相摩尔数的变化，在摩尔平衡方程中用成份 （elements）表示。PHERRQC 也允许如下平衡计算，通过增加或清除一种特定的反应物形 成具有纯相的平衡（EQUILIBRIUM_PHASES 数据块中可供选择的公式或相） ；计算形成纯 相平衡必须的反应物摩尔换算。 在这种计算模型中， 摩尔平衡方程中各项来自反应物的化学 计算而不是纯相的化学计算，未知数是进入或从溶液中析出的反应物摩尔数。 与每个新未知数相应的新函数是每个纯相的质量作用表达式。PHREEQC 中的溶解反应 在一定意义上，纯相是化学方程式左半边的部分。方解石的溶解反应可写为： CaCO3=Ca2++ CO32-, 并且，用 10-8.48 的 logＫ和纯固态活度１.0，质量作用表达式为 （33）K calcite ? 10 ?8.48 ? ? Ca 2? ? CO 2? .3（34） 一般情况下，纯相均衡可用下述方程表示:M aq mK p ? ? ? mm , p .c（35） 其中 cm,p 是溶解反应中主要物质 m 的化学计量系数，其值可正可负。PHREEQC 中溶解反应 左边各项系数为负，右边各项系数为正。无机物 SIp 饱和指数定义为 SIp ? log? ? mm , p .c mM aq(36)数值方法中相均衡函数是f p ? ?ln K p ? ?ln?10??SI p, target ? ? ? cm, p ln?? m ? ,mM aq（37）其中 SIp，target是相的目标饱和指数，ln(10)将以 10 为底的对数转换为自然对数。目标饱和指数由用户指定；正值、零、负值分别表示溶液中无机物过饱和、均衡、不饱和。对固定部 分压力气体成分，SIp，M aq mtarget 等于气体成分部分压力的对数，关于主要未知数的全微分是df p ? ? ? ? cm, p d ln ?m .（38）PHKEEQC 的数据输入，质量作用方程，均衡常数和纯相的温度要求用数据块 PHASES 确定。纯相集合的初始成分和目标饱和指数用数据块 EQUILIBRIUM-PHASES 确定。固体溶液均衡理想多成分或非理想，双成份固体溶液均衡的模拟是以 Glynn 的研究为基础的（Glynn 和 Readon,1990; Glynn 和其他人，1990；Glynn,1991; Glynn and Parklmst,1992） 。液相和固体 溶液均衡通过非齐次质量方程包括在模型中。 PHREEQC 允许多种固体溶液和液相在均衡中 存在，用固体溶液集合表示，受 Gibbs 相规则的限制，非理想固体溶液的模拟只限于双成份 固体溶液； 理想固体溶液可以有两种或多种成份。 固体溶液的附加主要未知数是固体溶液中 每种成份的摩尔数 n p ss ,其中 ss 指固体溶液 ss。该项说明了在元素摩尔平衡 Jacobian 矩阵中 每种成份摩尔数的变化。 不象纯相， 固体溶液中一种成份的活度不完全是 1.0。 其定义为 ? pss ? ? pss x pss ,其中 x p ss 是固体溶液 ss 中成份 p 的摩尔比值， ? p ss 是活度系数，固体溶液中一种成份的摩尔比值定 义为 x p ssn p ssp ss ?1?nN ss, 其中 Nss 是固体溶液 ss 中成份的摩尔数，对于理想固体溶液，活度系p ss数为 1.0，对于非理想双成份固体溶液中各成份的活度系数用 Gugenheim 表达式确定： λ 1=exp(( ?0 - ?1 (4x1-1))x22) 和λ 2=exp(( ?0 + ?1 (4x2-1))x12)，1（39）（40）其中λ 1 和λ 2 是成份 1 和成份 2 的活度系数，?0 和 ? 是无因次的 Gugenheim 参数，无因次 参数由过量自由能 g0 和 g1(kJ/mol)的量纲参数根据方程 ?0 =g0/RT 和 ?1 =g1/RT 计算得到。 过 量自由能参数 ?1 和 ? 2 可以直接定义，也可能通过许多方法间接确定，包括界定可混合性范 围的成份 2 的摩尔比值，界定 spinodal 范围的成份 2 的摩尔比值，临界点和临界温度处的成 份 2 的摩尔比值、Thompson 和 Waldbaum 参数、Mangules 参数，成份 2 的摩尔比值和一个 alyotropic 点的所有溶解物质的对数，成份 1 和 2 微量浓度的固相活度系数，或者成份 2 的 两个分布（配给）系数（glynn,1991） 。 相 应 与 每 个新 未 知 数的新 函 数 是 每个 固 体 溶液中 每 种 成 份的 质 量 作用表 达 式 。 PHREEQC 使用溶解反应，在一定意义上固定溶液成份指化学反应方程式左边各项。对于霞 石（文石）的文石-菱锶矿的固体溶液，溶解反应可写为 CaCO3=Ca2++CO32-， 使用 10-8.43 logK 和固体活度系数，形成质量作用表达式为:的(41)K Arag ? 10?8.34 ??Ca 2? ? CO 2?3? Arag??Ca 2? ?CO 2?3? nArag ? Arag ? ?n ? Arag ? nStront? ? ? ?.（42）一般情况下，固体溶液相均衡中，每种成份可用下列方程说明：M aq mK p ss ???c m , pss mK p ss ? p ss,（43） 其中 K pss 是纯相中成份 p 的均衡常数，cm, p ss 是固体溶液 ss 中成份 p 的溶解反应中， 主要物 质 m 的化学计量系数，其值可正可负。PHREEQC 中溶解反应左边各项系数为负，右边各 项系数为正。固体溶液一种成份的溶解部分（商）可定义为M aqQ p ss ???mc m , pssmK p ss ? p ss,(44)其中平衡时 Qpss 等于 1，ln Qpss 等于 0。对于非理想双成份固体溶液，其数值方法中每种成 分的函数为M aq mf1 ? ? cm,1 ln ?m ? ln K1 ? lnM aq mn1 ? ln ?1 和 n1 ? n2n2 ? ln ?2 . n1 ? n2（45）f 2 ? ? cm, 2 ln ?m ? ln K2 ? ln关于主要未知数的全微分是（46）2 2 3 aq ? x2 2? 0 x2 ? ? 6?1 x2 ? 12?1 x2 ? df1 ? ? cm,1d ln ? m ? ? ? ? ? n ?dn1 ? n ? n m 1 2 ? 1 ? 2 2 3 2? 0 x2 ? 2? 0 x2 ? 6?1 x2 ? 18?1 x2 ? 12?1 x2 ? 1 dn2 n1 ? n2M（47）和df2 ? ? cm, 2 d ln ? m ?mM aq2 2 ? 2? 0 x1 ? 2? 0 x2 ? 6?1 x1 ? 18?1 x12 ? 6?1 x2 ? 12?1 x13 ? 1 dn1 ? n1 ? n2? x1 2? 0 x12 ? 6?1 x12 ? 12?1 x13 ? 1 ? ? ? ?? n ? ?dn2 n1 ? n2 ? 2 ?。（48） 理想固体溶液中，在数值方法中每种成份的函数为 f pss ? ln Q p ss? M aq c m , pss ? ?? m ? ? ln? m K p ss ? ? ?N ss j ss? ? ? n p ss ? ? ?N ? ? ? ? ln? ? total ? ? ? ?，（49）其中 Ntotal ??nss，且 jSS 包括固体溶液 ss 中所有成份。其关于主要未知数的全微分是M aq mdf pss ? ? cm , p ss d ln ? m ?1 n p ss? N total ? n p ss ? ? N total ?N ss, j ss ? p ss ? 1 ? dn ? dn j ss 。 ? ? p ss N j total ? ss（50）PHREEQC 的数据输入， 质量作用方程， 平衡常数， 没仲春相常数的温度要求用 PHASES 数据块定义。一种固体溶液的初始成份和非理想固体溶液的 GUYYENHESIN 参数用 SOLID_SPLUTIONS 数据块确定（见D 数据输入说明” ）表面层的摩尔平衡表层处摩尔平衡是一般摩尔平衡方程的特殊例子。表层集合是由一个或多个表层组成， 每个表层有一个或多个位置类型。一个表层位置类型的总摩尔数由输入指定为如下之一： （1）固定值， （2）与一个纯相的摩尔数成正比，或（3）与动态反应物摩尔数成正比。下列 函数来源于表面 s 的一个表面位置类型 sk 的摩尔平衡关系。一个位置类型的表层物质占据 的表层位置的摩尔数总和必须等于表层位置类型的总摩尔数。 下列函数来源于表面 s 的一个 表层位置类型 sk 的摩尔平衡关系：f s k ? Ts k ? ? bs k ,i? sk ? ni? sk ?i ? sk ?N sk，（51）其中当达到平衡时，函数 f s k 的值为零，Ts k 是表层位置类型的摩尔数， N s k 是某种位置类型 表层物质的数量， bsk ,i? sk?是表层物质 i? s k ? 占据的表层位置的数量， f s k 的全微分为 df s k ? ?Ts k ? ? bs k ,i? sk ? dni? sk ? 。i ? sk ?N sk（52）如果表面位置总数和纯相的摩尔数成正比例，那么 ?TS K ? ?csk , p dnp ，其中 csk , p 是每 摩尔 p 相表层位置的摩尔数。如果相溶解，dnp 为正并且表面点减少。如果点的总数和动态 反应物质成比例，在全微分方程 ?TS K ? 0 中。点数的变化作为反应的一部分并与速率方程 相结合，且 Jacobia 矩阵中不包括任何项。随着动态反应增加或减少，反应物质摩尔数、表 面点数与此成比例调整。假如表面点数固定， ?TS K ? 0 。 PHREEQC 数据输入，表面点每种类型的摩尔总数用 SURFACE 数据块确定，且可能是 一个固定的量或者一个纯相的摩尔数或动态反应有关。表面点类型用 SURFACE_SPECIES 数据块确定并且表层物质可用 SURFACE_SPECIES 数据块确（见“ 数据输入说明‖ ）交换处的摩尔平衡一个交换点的摩尔平衡是一般摩尔平衡方程的特例。 一个交换点的总摩尔数由输入指定 为如下之一： （1）固定值， （2）与一个纯相的摩尔数成正比，或（3）与一个动态反应物摩 尔数成正比。 被交换物质占据点的摩尔数之和必须等于交换点的总摩尔数。 下列函数是一个 交换点的摩尔平衡关系：Ne iefe ? Te ? ?be,ie nie ,（53） 其中达到摩尔平衡时，函数 fe 值为零，Te 是交换器 e 中交换点的总摩尔数， be ,ie 是交换物 质占据的交换点的数量。fe 的全微分方程为dfe ? ?Te ? ?be,ie dnie .ieNe（54）如果点的总数与一个纯相的摩尔数成比例，那么 ?Te ? ?ce, p dnp ，其中 ce，p 是每摩尔 p 相交换点的摩尔数，如果相溶解，那么 dnp 为正并且交换点数减少。如果点的总数与动态反 应物摩尔数成比例，在全微分方程中 ?Te ? 0 。点数的变化与作为反应部分与速率方程相结 合，且没有相包括在 Jacobian 矩阵中。随着动态反应增加或减少，反应物摩尔数、交换点数 成比例增减。如果交换点数固定， ?Te ? 0 。 f e ? Te ??bieNee , ie ienPHREEQC 数据输入， 交换点摩尔数在 EXCHANGE 数据块中确定且可能是一个定量或 者与一个纯相或一个动态反应物的摩尔数有关。交换点用 EXCHANGE_SPECIES 数据块确 定且交换物质用 EXCHANGE_SPECIES 数据块确定（见“ 数据输入说明‖ ）碱度摩尔平衡碱度摩尔平衡方程仅用于物质形成计算和逆向模拟中。它是一般摩尔平衡方程的特例， 系数由每个水样的碱性确定。 PHREEQC 中碱性被确定为一种成分且一种主要物质与其相联 系（见“数据输入说明”中 SOLUTION_MASTER_SPECIES 关键词） 。在 PHREEQC 的默认 数据库中，碱性主要物质是 CO32-。碱度的主要未知数是 lnaAlk ,或默认数据库中 lnaco32-. 碱度当量总数由模型输入确定， 每种水样的碱性影响之和必须等于碱度当量总数。 下列 函数是碱度平衡方程：N aq ief Alk ? TAlk ? ? bAlk,i ni ,（55）其中达到摩尔平衡时，函数 fAlk 的值为零，TAlk 是溶液中碱度当量数，bAlk,I 是水样 i 的碱度 影响，fAlk 全微分方程为df Alk ? ?? bAlk,i dni .ieN aq(56)假定碳酸盐是主要碱性物质，TAlk 的值必须为正。从概念上讲，用 PHREEQC 计算的碱度不 同于实测的碱度。在默认 PHREEQC 的数据库文件中， bAlk,i 的值已选定，每种元素的所指 状态（ bAlk,i ? 0 ）或元素化合价态是 PH=4.5 时的主导物质。假定所有的元素或化学化合价 态在碱性滴定中被转化为这种主导物质。 然而所指状态不存在的水样的有效浓度 （具有非零 碱度的物质）可能在滴定终点存在，这在一定程度上使用 PHREEQC 计算的碱度与所测碱 度是一个不同的量。Fe 和 Al 的氢氧化物络合物是最普通的例子，可能不被转换为所指状 态。这样，一种溶液用 PHREEQC 计算的碱度尽管在数值上等于所测定碱度，实质上也是 一个近似值， 因为是假定滴定过程最终将元素和元素化合价态转变为它们的所指状态。 在许 多溶液中，碱度主要来源于碳酸盐，近似值是有效的。 PHREEQC 的数据输入，以每种物质的缔合反应计算的碱度在 SOLUTION_SPECIES 数据块中确定，主要物质的碱度由 SOLUTION_MASTER_SPECIES 数据块确定，总碱度是 溶液组成的一部分。用 SOLUTION or SOLUTION_SPREAD 数据块确定（见“数据输入说 明” ）元素摩尔平衡系统中一种元素的总摩尔数是下列各相中初始存在的摩尔数总和， 包括纯相、 固体溶液 集合、液相、交换集合、表面集合、气相和表面扩散层。下列函数是总摩尔平衡方程。Np SS N SS E Ne ? ? N aq ? ? f m ? ? Tm ? ? bm, p n p ? ?? bm, p n p ss ? ? ? bm,i ni ? ?? bm,ie nie ? p ss PSS e ie ? ? i??? bs k i ? sk ?SK s N skm , i ? s k ? i ? sk ?n? ? bm, g ng ? ?? bm,i ni , sg s iNgS N aq,（57）其中，达到摩尔平衡时，函数 fm 值为零，Tm 是系统中成份的总摩尔数，Np 是纯相集合中相 的数量，SS 是固体溶液集合中固体溶液数量，Nss 是固体溶液 ss 中成份的数量，Naq 是水样 数量，E 是交换集合中交换器的数量，Ne 是交换位置 e 的交换物质，S 是表层集合中表层个 数，Nsk 是表层类型 sk 的表层物质数量，Ng 是气相成份数量。系统中每个实体的摩尔数分 别用 np 表示纯相集合中相的摩尔数，n p ss 表示固体溶液的成分数， ni 表示水样数， nie 表 示交换位置 e 的交换物质数， ni? s ? 表示表层位置类型 sk 物质数，ng 表示气体成份数， nis 表k示表面 s 扩散层中水样数。每个实体中元素 m 的摩尔数用 bm 表示，用一个附加的下标说明 相关实体，bm 通常但不总是等于 cm（质量作用方程中主要物质 m 的系数） 。为了避免解决大数之间的微小差别，方程 57 中圆括弧中的量没有明确包括在求解算 法 中 ， 且Np pTm的 值 没 有 实 际 计 算 。 函 数SS N ss ss p ssfm中 使 用 代 替 量? Tm ? Tm ? ? bm, p n p ? ?? bm, p ss n p ss 。最初，*Tm 由液相 m 、交换集合、表面集合、气相和表面扩散层中总摩尔数计算得到：? Tm ? ? bm ,i ni ? ?? bm,ie nie ? ??? bm ,i? sk ? ni? sk ? ? ? bm , g ng ? ?? bm, i ni , s 。 （58）i e ie s k i ? sk ? g s iN aqENeSK s N skNgSN aq在方程的反复迭代中，通过纯相和固体溶液成份的摩尔转换使*Tm 得到更新：?Tm k ?1? ?T ? ? bm, p dnp ? ?? bm, p ss dnp ss ，m k p ss p ssNpSS N ss（59）其中 k 指迭代数，在中间迭代中*Tm 可能是负值，但达到均衡时必须为正。函数 fm 的全微分是dfm ? ?? bm, p dnp ? ?? bm, p ss dnp ss ? ? bm,i dni ? ?? bm,ie dnie ?p ss p ss i e ieNpSS N ssN aqENe??? bs k i ? sk ?SKsN skm , i ? sk ?dni? sk ? ? ? bm, g dng ? ?? bm,i dni , sg s iNgS N aq（60）PHREEQC 的数据输入，元素的总摩尔数最初用 SOLUTION or SOLUTION_SPREAD 数据块确定，对于交换集合，用 EXCHANGE 数据块确定，对于表面集合，用 SURFACE 数 据 块 ， 对 于 气 相 ， 用 GAS_PHASE 数 据 块 ， 纯 相 集 合 中 每 个 相 的 摩 尔 数 用 EQUILIBRIUM_PHASES 数据块确定。在一个固体液体集合中每种固体溶液每种成份的摩 尔数 SOLID_SOLUTIONS 数据块确定。 元素的总摩尔数也可在批反应和迁移计算中改修动 （见“数据输入说明” ） 。溶液中电荷平衡电荷平衡方程计算了水中阴阳离子的当量总和。在一些情况下，表面和交换器电荷不 平衡。当具体指定时，在初始溶液计算中用电荷平衡方程调节 PH 值或一种主要物质的活度 （继而一种元素的总浓度或元素的化合价态） 使溶液达到电中性。 电荷平衡方程对于在批量 反应和迁移模拟中计算 PH 值是必要的。 在现实的溶液中， 阴阳离子的当量和必须为零。 然而在化学分析中分析误差和没有分析 的成份通常引起溶液计算中电荷不平衡。 如果一个初始溶液中计算到电荷不平衡， 在后来的 批量反应和运移模拟中需调节 PH 值，维持同样的电荷不平衡。将溶液混合，批反应的电荷 不平衡是每种溶液其混合因素产生电荷不平衡偏差之和。 如果模拟中用到一个表面， 且没有 指定显扩散层计算，那么带电的表层物质形成将产生表面电荷不平衡。相似地，如果交换物 质不是电中性（在默认数据库中所有交换物质是电中性） ，那么交换器将积累电荷。一般的 电荷平衡方程包括表面和交换器电荷不平衡。 在每个初始溶液和批反应步骤中， 对于每个单元在运移模拟的每个时间区间中计算电荷 不平衡使用方程Tz ,q ? ? zi ni ，iN aq（61）其中 q 代表液相， Tz , q 是液相 q 的电荷不平衡，zi 是水样 i 所带电荷。如果带电的表层和交 换器不存在，在溶液批反应中或运移模拟终点的电荷不平衡将和模拟始点相同。 表面的电荷不平衡在初始表面组成计算和每个批反应步骤中计算，对于每个单元运移 模拟的每个时间区间，使用方程Tz , s ? ?? zi? sk ? ni? sk ? ? ? zi ni , s ，k i ? sk ? iK s N skN aq（62）其中 Tz，s 是指表面电荷不平衡, zi? s ? 是表面 s 的表面类型 sk 的表层物质 i 所带电荷，方程中k最后一项代表弥散层中积聚的电荷。 最后一项仅用于弥散层组分明确包括在计算中时 （在数 据块 SURFACE 中的-diffuse-layer） 。详细计算弥散层组分时，要求所有的溶液电荷平衡，并 且 Tz，s 总为零。 正常情况下，交换物质没有静电荷，但是对于大部分，这不作要求。然而，如果带电 物质的和不等于交换点（交换能力）的当量总数，交换物质（当量积分）的活度不能很好确 定。 如果带电物质存在， 那么在每个批反应步骤中的初始交换组分计算中计算交换器的电荷 不平衡，对于每个单元，运移模拟的每个时间区间用方程Tz , e ? ? zie nieeiNe， dfm ? ?? bm, p dnp ? ?? bm, p ss dnp ss ? ? bm,i dni ? ?? bm,ie dnie ?p ss p ss i e ieNpSS N ssN aqENe??? bs k i ? sk ?SK s N skm , i ? sk ?dni? sk ? ? ? bm, g dng ? ?? bm,i dni , sg s iNgS N aq（63） 其中 Tz，e 是交换器的电荷不平衡， zie 是交换器 e 中交换物质 i 所带电荷。 系统电荷不平衡批反应步骤之初确定， 并且对于每个单元， 在运移模拟中每个时间段始 点应该是：Tz ? ? ?qTz , q ? ?Tz , s ? ?Tz ,e ，q s eQSE（64）其中 Tz 是系统电荷不平衡，Q 是在每个批反应步骤或一个运移步骤单元中混合的液相数量，?q 液相 q 的混合比值，S 是表层的数量，E 是交换器数量。电荷平衡函数是：S N aq ? S K s N sk ? E Ne ? f z ? Tz ? ? zi ni ? ??? zi? s ? ni? s ? ? ?? zi ni , s ? ? ?? zie nie ， k k ? s k i ? e i i s i ? sk ? e ? ? N aq（65）其中达到电荷平衡时 fz 为零。 如果计算弥散层组份， 每个表面包括一个独立的电荷平衡方程， 达到电荷平衡时，括号内项目的和将为零。如果弥散层组份不计算，括号内每二项为零。fz 的全微分方程为:df z ? ?? zi ni ? ??? zi? sk ? ni? sk ? ? ?? zie nie ，i s k i ? sk ? e ieN aqSK s N skENe（66）其中，只有在弥散组份没有详细计算时，表面三项的总和存在。 PHREEQC 的数据输入，用 SOLLUTION 或 SOLLUTION-SPREAD、EXCHANGE 和 SURFACE 数据块的数据输入结合物质组成，初始交换组份和初始表面组份计算确定。在 SOLLUTION-SPECIE、EXCHANGE-SPECIE 或 SURFACE-SPECIE 数据块中确定化学物质 反应的物质，电荷由平衡的化学反应确定（见“数据输入说明” ） 。不计弥散层组份的表面电荷电势方程根据缺省值，PHREEQC 应用 Dzombak 和 Morel（1990）说明的方法将表面电荷密度ζ s和表面电势Ψ s 相联系。表面电荷密度是每单位面积表面物质的电荷量，可由表面物质的分布计算而得：F ?s ? Asurf?? zk i ? sk ?K s N ski ? sk ? i ? sk ?n，（67）其中ζ s 是表面 s 的电荷密度（C/m2） ，F 是法拉第常数（96，493.5C/mol） ，Asurf 是物质的表 面积（m2） 。表面积由下列公式之一计算： （1）Asurf=ASSS，其中 As 是表面物质的比表面积 （m2/g） ，SS 是表面物质的质量（g） ，或（2）Asurf=Arnr，其中 Ar 是每摩尔纯相或热动力反 应物的表面积（m2/mol） ，nr 是纯相或热动力反应物的摩尔数。在 25°C 时，表面电荷密度 与表面电势如下：1 1 ? vF?s ? ?s ? ?8000 ?? 0 RT ?2 ???2 sinh? ?， ? 2 RT ?（68） 是自由空间的渗透性（8.854*10-12CV-1m-1 或其中ε 是水的介电常数（78.5，无量纲） ，ε0C2/m-j） ，v 是均质电解溶液的离子电荷，R 是气体常数（8.314Jmol-1K-1） ，T 是温度（K） ， μ 是离子强度，F 是法拉第常数（JV-1eq-1 或 C/mol） ，Ψ s 是用伏特计算的表面电势。在 25° C 时， ?8000 ?? 0 RT ?2 ? 0.1174。电解溶液的离子电荷假定为 1。 电荷-电势函数为：1 1 F ? vF ?s ? f ?s ? ?8000 ?? 0 RT ?2 ?? ?2 sinh ? ?? ? 2 RT ? Asurf1?? zk i ? sk ?K s N ski ? sk ? i ? sk ?n，（69）并且该函数的全微分为：df?s ??8000?? 0 RT ? ?? ?? 1 ? vF?s ? 2 sinh? ?d? ?2 ? 2 RT ?1 2 1 ? 21 2?8000?? 0 RT ? ?? ?F ? F?s ? cosh? ?d ln ? ?s ? ? Asurf ? 2 RT ??? zk i ? sk ?K s N sk。i ? sk ? i ? sk ?（70）nPHREEQC 的数据输入， 默认值不计弥散层组分。 比表面积 （As 或 Ar） 和表面质量 （Ss） 在 SURFACE 数据块中确定，表面位置摩尔数（1）如果数量固定在 SURFACE 数据块中确 定， （2）根据 SURFACE 数据块中的均衡系数，单相的摩尔数在 EQUILIBRIUM-PHASES 数据块中确定，或（3）根据 SURFACE 数据块中的均衡系数，热动力反应物的摩尔数在 KINETICS 数据块中确定。表面物质的电荷在平衡的化学反应中确定，并且该物质在 SURFACE-SPECIESE 数据块中确定（见“数据输入说明” ） 。 计算弥散层组份的表面电荷电势方程作为表面电荷 -电势关系先前模型的替换， PHREEQC 选择使用 Borkovec 和 Westall （1983）发展的方法。他们的发展在氧化电解液界面处求解 Poission-Boltzmann 方程来确定 弥散层中表面的过量离子。假定 1 升水为 1kg，利用下列微分计算。 表面过量离子为：??i , s ?x d ,s? ?c ?x ? ? c ?d xi, s 0 i，（71）其中 ?i , s 为表面 s 上每摩尔每平方米水样 i 的过量离子，xd，s 是外层 Helmholtz 平面的位置， ，作为表面的距离函数，ci0 是溶液的体积浓度。过量离子在 1.0kg ci, s ?x? 是浓度（mol/m3） 水的所指状态中与浓度有关： mi，s=Asurf ?i , s ， （72）其中 mi，s 是水样 i 的过量离子 （mol/kgw） ， 该表面过量离子浓度可在体积溶液中与浓度相联 系： mi，s= gi，smI， 其中 gi，s 是体积溶液中表面电势、所有离子的浓度和电荷的函数： （73）gi , s ? Asurf sign?X d , s ? 1??X d ,s?1?Xzi?1? ?1/ 2? 2 N aq ? z ? X ? ml X i ? 1 ? l ? ? ? ??dX ，（74）其中 X ? e?F?s RT，Xd，s 是外 Helmholtz 平面上 X 的值，Asurf 是表面积（m2） ， （Xd，s-1）为+1或-1，取决于括号中该项的符号，i 是所计算的表层过量离子的水样，zi 水样 i 所带电荷没， l 的范围为所有的水样，zl 是水样 l 所带电荷，并且 ? ? ??? 0 RT / 2?1/ 2。在 25°C 时，的值为 0.02931（l /mol）1/2Cm-2。Borkovec 和 Westall 使用的未知数 X 与 PHREEQC 所用的主 要未知数的关系为a?s ? X ?2 。Borkovec 和 Westall（1983）的发展仅计算了弥散层中每个水样的过量浓度。当一种 溶液从弥散层分离出时，在批反应和运移模拟中就会产生问题，例如，在平流模拟中，水从 一个水槽平流进另一个水槽。在这种情况下，需要知道保留在表面的总水样 i 一个经验假定 是弥散层为指定厚度并且所有的表面过量离子留在弥散层。 弥散层中一个水样总摩尔数等于 弥散层中表面过量离子的有效因素加上体积溶液： ni，s= ni，s，excess+ ni，s，aq= Wbulk+gi，sni n n ? Ws i ? gi , s ni ? Ws i Waq Waq Waq（75）其中 ni，s，aq 指由于体积溶液的存在弥散层中水样 i 的摩尔数， ni，s，excess 指表面过量离子的摩 尔数，Waq 是除了弥散层之外系统中水的质量，Ws 是表面 s 弥散层中水的质量，假定液相中 水量远大于弥散层中水量，使得 Wbulk≌ Waq， （版本 1 中，Wbulk= Waq+ 由弥散层的厚度和表面积计算而得，假定 1 升水为 1kg： Ws=tsAsurf （76） 其中 ts 为弥散层的厚度， 计量单位为米。 如果表层位置的摩尔数与一个纯相或热动力反应物 的摩尔数有关，那么 Asurf=Arnr，否则 Asurf 为常数并且由输入指定的比表面积和表面的质量 计 算而 得。根 据静 电理论 ，弥 散层的 厚度离 子强 度低 时较大 ，离子 强度 高时 较小。 PHREEQC 中弥散层的厚度的默认值为 1*10-8m，大约为由德拜理论计算的 0.001 克分子离 子强度的厚度，弥散层的德拜长计算为 1*10-7 m。如果表面积足够大，弥散层水量很小的假 定将水有效的；厚度为 1*10-7 m，表面积为 1000 m2 将造成弥散层 0.1 升的水量，这是 1 升 体积溶液有意义的一部分。 弥散层水样摩尔数的全微分方程为：? F?s ? ? ? ? ? 2 RT ? ? ?gi , s ? W ? s ? ? ?dni ? ni dni , s ? ? gi , s ? ? 2 e ? ? ? W ? X ? aq ? ? ??2?WsSs） 。水的质量? ?d ln ? ? n Ws d ln W ? n t s Ar dn ?s i aq i r， ? Waq Waq ? ?（77）其中第二项是关于表面电势主要未知数的偏微分，lnaΨ s 和最后一项仅在表面位置数量与一 纯相或热动力反应物有关时存在。偏微分?g i , s 等于方程 74 在 Xd，s 处的积分： ?Xzi d ,s?gi , s ?XX d ,s? Asurf sign?X d , s ? 1???X?1? ?1/ 2? 2 N aq ? z ? X d , s ? ml X d l, s ? 1 ? l ? ? ? ??,（ 78 ） 函数 gi，s 关于主要未知数的偏微分为：?2 ? F?s ? ? ? ? ? ? ?g i , s ?g i , s ? ? ? 2 RT ? ? ? ? 2e ? ? ? ln ? ?s ?X ? ? ? ? ? F?s ? ? 2 ? ? ? ? 2 RT ? ? X dz i, s ? 1 ? ? ? ? ? ? Asurf sign? X d , s ? 1??? 2e 1/ 2 ? ? ? ? 2 N aq ? ? ? ? X ? m X zl ? 1 ? d ,s l d ,s l ? ? ? ???.（79）??在数值方法中， gi，s 的计算很复杂， PHREEQC 版本使用和 Borkovec 和 Westall （1983） 相同的方法以减少函数计算的数量。 当弥散层明确包括在计算中时， 增加了一个新的迭代运 算。在每个弥散层迭代运算开始，函数和它们的偏微分必须明确求出。在弥散层的模拟迭代 中，函数值用下列方程迭代计算：?1 gik, s ? gik, s ??gi , s d ln ??s , ? ln ??s（80）其中 k 是模拟迭代次数，g0i，s 弥散层迭代开始求得的值。当牛顿方法收敛于一个解时，迭代 结束；然而，收敛是以函数 gi，s 的估计值为几基础的。这样，弥散层迭代一直继续到函数值 在指定允许范围内与逐次弥散层迭代相同。 当明确计算弥散层组分时，涉及电势未知数 sinh（方程 69）的函数用包括表面电荷和 弥散层电荷的电荷平衡函数代替：f z , s ? ?? zi? sk ? ni? sk ? ? ? zi ni , s ,k i ? sk ? iK s N skN aq（81）其中当达到电荷平衡时，fz，s 为零。fz，s 的全微分方程为：df z , s ? ?? zi? sk ? d ni? sk ? ? ? zi d ni , s .k i ? sk ? iK s N skN aq（82）PHREEQC 的数据输入，在 SURFACE 数据块中用-diffuse-layer 标识符计算弥散层。比表面 积（As 或 Ar）和表面质量（Ss）在 SURFACE 数据块中确定，表面位置摩尔数（1）如果数 量固定在 SURFACE 数据块中确定， （2）根据 SURFACE 数据块中的均衡系数，单相的摩尔 数在 EQUILIBRIUM-PHASES 数据块中确定，或（3）根据 SURFACE 数据块中的均衡系数， 热动力反应物的摩尔数在 KINETICS 数据块中确定。表面物质的电荷在平衡的化学反应中 确定，并且该物质在 SURFACE-SPECIESE 数据块中确定（见“数据输入说明” ） 。 非静电表面络合Davis 和 Kent（1990）说明了一个非静电表面络合模型。在这个模型中，表面络合质量 作用表达式中静电项被忽略。另外，不使用表面电荷平衡或表面电势关系；仅包括每个表面 位置类型的摩尔平衡方程。 PHREEQC 的数据输入， 在 SUQFACE 数据块中通过用-no-edl 标识符时产生非静电表面 模型（见“数据输入说明” ） 。物质形成和正向模拟的数值方法以前个部分用 PHREEQC 解决的任何化学均衡问题都来自用 f 表示的方程组。 包括 fAlk， fe，fg，fH，fH2O，fm`，fO，fPtotal，fp，fp，fpss，fsk，fz，fz，s，fu 和 fΨ s，其中 fH 和 fO 是氢和氧的 摩尔平衡函数，m`指所有水中的主要物质，除了 H+，e-，H2O 和碱性物质。相应的主要未 知数是， ln ? e ，ng， ln ?e ? ， ln ? H 2 O ， ln ?m` ，lnWaq，Ngas，np（或在物质形成计算中可 能为 ln ?m` ）nss， ln ?sk , ln ? H ? （或在物质形成计算中可能为 ln ?m` ） ， ln ??S （明确的 弥散层计算中） ， ? 和 ln ??S （隐含的弥散层计算中） 。当一给定计算的所有函数的偏差等 于零，即找到了非线性方程组的一个解，化学系统的均衡值已确定。 （注意：如果一纯相或 气相在均衡中不存在，出事包含在给定计算中的一些方程可能会遗漏。 ）解答方法赋予主要 未知数初始值，然后用牛顿-落弗逊方法不断修改直到在指定误差允许范围内找到方程的一 个解。 对于一组方程 fi=0， xj 是未知数， 牛顿-落弗逊方法需不断修改初始未知数组的值。 使 ? i =fi 方程未知数目前值的偏差。形成下列方程组：? i ? ??jJ?fi dxj , ?x j（83）其中 J 是所计算主要未知数的总数。方程组是线性的，可以同时求解未知数 dxj。未知数的 新值由公式计算：k ?1x j ?k x j ? dxj ，k 指迭代次数，然后计算偏差的新值。重复此过程直到偏差小于指定允许误差。 使用牛顿-落弗逊方法 解决化学均衡问题，需注意两点。第一，未知数的初始值必须充 分接近于均衡值，或该方法不收敛，第二，如果一组相的化学反应是非线性相关的，可能形 成 一个奇异矩阵。PHREEQC 使用 Barrodale 和 Roberts（1980）发展的优选技术避免出现 奇异矩阵。 该问题的优选技术也允许不等式约束， 对于限制相和能反应的固体溶液的总两是 有用的。 下面说明了每个模拟类型主要未知数初始值的选择问题。不顾赋予初始估计值的方法， 如果必要， 在牛顿-落弗逊方法产生近似的摩尔平衡之前修改元素和元素价态活度的估计值。 水成主要物质的程序如下。初始估计值确定后，计算每种元素（除了氢和氧）和初始溶液中 确定的各个原子价态的物质分配。然后，求出计算所得摩尔数和输入摩尔数的比率。如果一 种主要物质 m`的比率大于或小于 10-5，用下列方程矫正主要未知数的值：? N aq ? ? ? b ` ni ? ? i m ,i ? , ?1 k ?1 ln ? k ? ln ? ` ` ? w ln m m ? T ` ? m ? ? ? ? ? ?（84）其中如果比率大于 1.5，w 为 1.0，比率小于 10-5，w 为 0.3，Tm`是一种元素或元素原子价状 态的总浓度。交换和表面主要物质使用累似的方程。初始估计值矫正后，计算物质分配。继 续迭代直到比率在指定范围内，在该点，使用修改的牛顿-落弗逊技术。如果连续的修改不 能找到活度使得比率在指定范围内，那么，第二组迭代试图将比率减小到低于 1.5，这些比 率无更小的限制。无论第二组迭代是否成功，都用牛顿-落弗逊技术。 Barrodale 和 Roberts（1980）优选技术是简单线性算法的修正，使受相等和不等约束的 线性方程组偏差绝对值（L1 优选）的和最小。一般问题可用下列矩阵方程表示： AX=B CX=D EX≤F 第一个矩阵方程在一定意义上被减至最小， 使 （85）? b ? ??i i jIBJi, jx j 为一极小值，其中 IB 是 被优选方程的数量，受第二个矩阵方程相等条件约束和第三个矩阵方程不等条件约束。 PHREEQC 方法在优选方程的（AX=B）中将包含一些牛顿-落弗逊方程，而不包括所 有的牛顿-落弗逊方程，如等式（CX=D） ，在给定迭代中，包含在 A 矩阵中的方程可能不能 求得确切的等式，但是在一定意义上将得到优化。这样，在给定迭代中，即使无确切的等式 解存在，也能找到牛顿-落弗逊方程组的一个近似数学解，例如，所有方程的压力等式造成 一个无解的奇异矩阵。矩阵 A 包括碱度方程，气相中气体的总摩尔数，纯相和固体溶液组 分。矩阵 B 包括所有的摩尔平衡，电荷平衡和表面电势方程。矩阵 C 包括限制纯相的溶解、 固体溶液组分和气体成分在系统中存在程度的不等式。 如果相应未知数（在元素主要物质活度的对数中产生变化）的所有系数（一个直列矩阵 A 和 B）都小于 1e-10，为避免一些与矩阵 B 中小数字有关的数值问题，代表摩尔平衡方程 的一行矩阵将被改变（按比例缩放） 。在这种情况下，方程与 1e-10 除以最大系数的绝对值 成比例。当指定时（在 KNOBS 数据块-diagonal-scale 中） ，如果未知数的系数小于 1e-10， 可选择使用另一种方法，摩尔平衡方程将与 1e-10 除以相应未知数的系数成比例。 改变的矩阵用优选的解答器解答，返回值是一主要未知数值的变量。检查变量的值，确 保未知数变量小于限制变量最大值允许的范围的标准。这些标准在在程序中由默认值或在 KNOBS 数据块中的输入指定。如果变量太大，那么所有未知数的变量，除了纯相和固体溶 液组分的摩尔转换，按比例减少，使其满足所有标准。纯相和固体溶液摩尔转换不改变，除 了产生纯相和固体溶液组分总摩尔数的非负值。 计算未知数相适应的变量后， 更新主要未知 数，计算所有水样、交换物质和表面物质新的克分子浓度和活度以及所有函数的偏差 。偏 差用收敛验证（收敛判据在程序中确定，但可在 SELECTED-OUTPT 数据块 KNOBS 或 -high-precion 选择中-convergence-tolerance 转换为一替代组，如果没有收敛，则开始新的迭 代。水成物质计算水成物质计算用溶液的一种化学组分作为输入， 并计算水样的配比和相的饱和指数。 水 成物质计算包括方程 fH2O、fm`、和 fu，即元素或元素化合价态的摩尔平衡方程、水的活度和 离子强度方程。不包括氢和氧的摩尔平衡方程，因为一般不知道氢和氧的总质量。相反水的 质量假定为 1.0kg 或指定（在 SOLUTION 或 SOLUTION-SPREAD 数据块-water 中） ，氢和 氧的总质量在 水成物质计算中由水的质量和包含水样的氢和氧的浓度计算得到。 如果 PH 值、pe 或一种元素或元素化合价态的主要未知数指定为调整到溶液的电荷平 衡，fz 用于计算产生电荷平衡的主要未知数（ ln ? H ? 、 ln ?e ? 、 ln ?m` ）的值。在这种情况 下，计算得到的 PH、pe 或 m`的总浓度将不同于输入值。如果用 fz 计算主要未知数 ln ?m` ， 则不用方程 fm`。 如果输入中指定总碱度，用碱度摩尔平衡方程计算 ln ? Alk 和与碱度有关的元素（默认 数据库中为碳）的总克分子浓度。如果定义包括一个碳[碳（+4）]和碱度两者的摩尔平衡方 程，那么与这些方程有关的两个主要未知数是 ln ? Alk ? ln ?CO ?2 （默认数据库文件中）3ln ? H ? 。在这种情况下，PH 值将在水成物质计算得到并且不等于输入的 PH 值。对于物质形成计算， 如果在问题的阐述中包括碱度摩尔平衡方程， 则它为解答器唯一的 优选方程。所有其他的方程作为相等约束条件，物质形成计算不包括不等约束条件。 在 初 始 溶 液 计 算 中 允 许 部 分 氧 化 还 原 不 平 衡 反 应 ， 在 SOLUTION 或 SOLUTION-SPREAD 数据块中氧化还原反应的选择影响回程物质形成和饱和指数计算。根 据默认值，无论何时需要用电子活度计算一个水样的克分子浓度和活度，使用输入的 pe。 如果给出缺省的氧化还原对（-redox）或指定某元素的氧化还原对（或元素化合价态的组合） （见数据输入说明中的关键字 SOLUTION） ，然后重写每个水样氧化还原元素的质量作用表 达式，用氧化还原对的活度代替电子的活度。例如，使用硫酸盐-硫化物氧化还原对[S（+6） /S（-2）]铁离子将被改变，Fe+3 原始的化学反应： Fe+2 =Fe+3+e用硫的缔合反应写为， SO4-2+9H++8e-=HS-+4H2O， 生成下列不含电子的化学反应： （87） （86）1 ?2 9 ? 1 1 Fe ? 2 ? SO4 ? H ? Fe ? 3 ? HS ? ? H 2O . 8 8 8 2（88）最后的质量作用表达式将被用作 Fe+3 的质量作用表达式，Fe+3 摩尔数的微分 dnFe?3 ，也将以 这个质量作用表达式为基础。然而，原始的质量作用表达式（以方程 86 为基础）用于确定 含有 dnFe?3 的摩尔平衡方程，也就是说，Fe+3 将出现在铁的摩尔平衡方程中，而不出现在 S （+6）或 S（-2）的摩尔平衡氧化还原反应中。另一组方程中。这些变换的结果是 Fe+2 、 Fe+3、硫酸盐和硫化物处在平衡的氧化还原反应元素（例如，氧和氮）也可以形成在它们自 己的氧化还原反应平衡中，而不必与铁和硫一起生成氧化还原平衡反应。 根据默认值，如果饱和指数计算需要 pe（或电子的活度）的值，那么，用 pe 的输入值。 如果已确定氧化还原对（-redox） ，那么该相的溶解反映写法如上，消除 电子活度而用氧化 还原对的活度代替。 在一个计算中，主要未知数组可能随氧化还原元素的变化而变化。如果作为一个摩尔平 衡方程的主要未知数的主要物质活度变为小于同一摩尔平衡方程另一主要物质活度数量级 的十倍，基础转换过程将发生。在这种情况下，涉及目前主要未知数的质量作用表达式（包 括水样、交换物质、气体、表面物质和纯相）被写在具有更大活度的新的主要物质的项中。 这种过程的另一个离子为，如果氮在一个系统中减少，那么氮的主要未知数将从硝酸盐（在 氮减少的条件下，其量可忽略）转换为铵（此时将作为主要物质在） 。基础转换过程不影响 物质的最终平衡分配，但是在处理低浓度问题时，加快了计算速度并避免了数值问题。 估计主要未知数的初始值并且根据前一部分说明的方法修正。对于初始溶液计算，PH 值、pe 的输入值用作初始估计。水的质量为 1.0kg，除非有指定值，水的活度估计为 1.0。 假定主要物质是存在的唯一物质，并且它们的浓度等于输入浓度（转换为克分子浓度） ，估 算离子强度。元素（除了氢和氧）和元素化合价态主要物质的活度恒等于输入浓度（转换为 克分子浓度） 。如果用电荷平衡方程和相均衡方程代替一种元素和元素化合价态的摩尔平衡 方程，那么主要物质的初始活度恒等于输入浓度（转换为克分子浓度）的千分之一。 PHREEQC 的数据输入，物质形成计算使用的碱度方程、电荷平衡方程、相均衡方程和 氧化还原对在 SOLUTION 或 SOLUTION-SPREAD 数据块中确定（见“数据输入说明” ） 。交换器初始组分的计算如果一个交换器的组分没有明确给定，但是指出与一指定溶液组分在均衡中存在，那么 就需要计算初始交换组分。 在这种情况下， 交换器的组分未知， 只知道其与一溶液处于均衡。 初始交换组分计算方程是 fe、fm`、 f H 2 O 、和 fu，即每个交换器的摩尔平衡方程、每种元素和 元素化合价态的摩尔平衡方程、水的活度和离子强度方程。 对于初始交换组分计算，Tm`的值仅包括水样的浓度，fm 不包含交换器对总元素浓度的影 响。所有与液相有关的量和没有交换器存在时的溶液相同。只有交换集合未知数 ln ? e 的值 调整为达到交换器的摩尔平衡 。一旦达到摩尔平衡，即知每个交换器的组分。 所有初始交换组分计算方程作为解答器的相等条件约束。没有方程得到优选且不包括不 等条件约束。 只有交换器与指定溶液处于均衡时，才需计算初始交换组分。该溶液的物质配比已通过 初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关的未知数的值并且作为交 换计算的初始估计值。每个交换器的主要未知数的初始估计值恒等于该交换器交换的摩尔 数。 PHREEQC 的数据输入， 初始交换组分计算用 EXCHANGE 数据块 （见 “数据输入说明” ） 。表面初始组分的计算如果一个表面组分没有明确给定，但是指出与一指定溶液组分在均衡中存在，那么就需 要计算初始表面组分。在这种情况下，表面的组分未知，只知道其与一溶液处于均衡。初始 表面组分计算方程是 f s k 、 f ?s 或 fz，s、 fm`和 fu， 即表面集合中每种表面类型的摩尔平衡方程、 每个表面的电荷电势关系或每个表面的电荷平衡方程（非静电模型中不包括这些方程） 、每 种元素和元素化合价态的摩尔平衡方程、水的活度和离子强度方程。 对于初始表面组分计算，Tm`的值仅包括水样的浓度，相应的摩尔平衡方程 fm`不包含表 面对总元素浓度的影响。所有与液相有关的量和没有表面集合时的溶液相同。 计算弥散层组分时，每个表面使用电荷平衡方程 fz，s；表面集合每个表面类型主要未知 数 ln ?sk 的值调整为达到每个表面的摩尔平衡和电荷平衡 。如果不计算弥散层组分时，那 么用电荷电势方程 f ?s 代替表面电荷平衡方程。如果非静电模型用于表面集合，那么表面电 荷平衡方程和电荷电势方程都不包含在待解的方程组中。 所有初始表面组分计算方程都作为解答器的相等条件约束。没有方程得到优选且不包括 不等条件约束。 只有初始表面与指定溶液处于均衡时，才需计算初始表面组分。该溶液的物质配比已通 过初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关的未知数的值并且作为 表面计算的初始估计值。 每个表面主要主要物质活度的初始估计值恒等于该表面的表面位置 摩尔数的十分之一。 对于明确和隐含弥散层计算， 每个表面电势未知数 ln ??s 的初始估计值 为零，暗示了表面电荷为零。 PHREEQC 的数据输入，初始表面组分计算用 SURFACE 数据块（见“数据输入说明” ） 。固定体积气相初始组分计算如果气相组分没有明确确定，但是一直定体积气相组分与一特定溶液组分处于均衡中 时， 需要计算初始气相组分。 初始气相组分计算方程与初始溶液组分计算方程相同， 是 fm` 、 f H 2O 和 fu，即每种元素和元素化合价态的摩尔平衡方程、水的活度和离子强度方程。对于初始气相组分计算，Tm`的值仅包括水样的浓度，相应的摩尔平衡方程 fm`不包含气 体对总元素浓度的影响。 所有与液相有关的量和没有气相存在时的溶液相同。 一旦确定了液 相中组分的物质配比，气相中所有组分的部分压力即可计算。在气相中根据理想气体定 f ?s 律用部分压力和指定的体积计算每种成分的摩尔数。 所有初始气相组分组分计算方程都作为解答器的相等条件约束。没有方程得到优选且不 包括不等条件约束。 只有气相被确定具有固定体积且与指定溶液最初处于均衡时，才需计算初始气相组分。 该溶液的物质配比已通过初始溶液或批反应或运移计算得到。 这样， 知道了所有与液相有关 的未知数的值并且作为初始气相组分组分计算的初始估计值。 PHREEQC 的数据输入， 初始气相组分计算用 GAS-PHASE 数据块 （见 “数据输入说明” ） 。批反应和运移模拟批反应和运移模拟计算要求计算液相和化学系统中存在的任何均衡相集合、表面集合、 交换器集合、固体溶液集合和气相的均衡。发生在均衡之前的不可逆反应包括混合、理想配 比反应、热动力反应和温度变化。批反应和运移模拟计算中完整的牛顿-罗弗逊方程组包括 fe， fH， f H 2 O ，fm`，fO， f ptotal ，fp， f p ss ， f s k ，fz，fz，s，fu 和 f ?s 。 总包括氢的摩尔平衡方程 fH、水的活度 f H 2 O 、氧的摩尔平衡方程 fO、电荷平衡方程 fz、 离子强度 fu，并且与总作为主要未知数的 ln ? H 2 O 、Waq（水的质量） 、 ln ? H ? 和 ? 有关。 摩尔平衡方程 fm 用于求元素的总浓度，而非各个化合价态或化合价态组合的浓度。不可 能包括碱度摩尔平衡方程；它用于初始溶液计算。 假如给定定压力气相且在均衡中存在，则包括方程 f ptotal 。假如给定交换集合，则包括方 程 fe。假如给定表面集合，则包括方程 f s k 。另外，假如给定每个表面的隐含弥散层计算， 则包括方程 f ?s 。假如给定明确的弥散层计算，则包括方程 fz，s。包括均衡中存在的每个纯 相的方程 fp。包括均衡中存在的每种固体溶液每种成分的方程 f p ss 。 在计算之初， 已知一纯相、 固体溶液或定体积气相是否在均衡中存在。 这样， 在每次迭 f p ss 代中，用下列推理确定均衡计算应包括哪个方程。如果一个相的摩尔数大于零，np&0, 或者 计算所得饱和指数大于目标饱和指数，则包括其方程。如果矩阵中不包括该方程，那么未知 数 dnp 所有的系数恒等于零。 对于一理想固体溶液，如果任一组分的摩尔数大于 1*10-13 或其和?p ssIAPp ss K p ss大于 1.0，则包括方程 f p ss 。对于一理想固体溶液,IAPp ss K p ss? x p ss 如果摩尔级分的和大于 1.0，则确定求和。如果固体溶液方程不包含在矩阵中，那么未知数 dnPss 所有的系数恒等于零。 对于非理想双成分固体溶液， 决定是否包括固体溶液方程的下列过程由 Glynn 和 Reardon 方程（1990，方程 37 到 48）发展而来。如果任一固体溶液组分的摩尔数大于 1*10 ，则包 括所有的固体溶液方程。否则，该组分水的活度级分由下式计算：-13x1, aq ?IAP IAP 1 2 和 x1, aq ? IAP ? IAP IAP ? IAP 1 2 1 2（89）其中 IAP 是离子活度积。下一步通过求解下列方程得到 x1 和 x2 的值从而确定将与水的活度 级分处于均衡的固体的摩尔级分：x1?1K1 ? x2? 2 K 2 ?1 x1, aq ?1K1 ? x2, aq ?2 K2，（90）其中 x1 和 x2 是在固相中的摩尔级分，K1 和 K2 是单成分的均衡常数， ?1 和 ? 2 是成分的活度 系数，据过量自由能的 Guggenheim 参数计算而来。该方程高度非线性且用第一个测试区间 [0，1]求得一满足方程包含摩尔级分的解，然后使区间缩小一半，使摩尔级分的估计值更精 确。一旦确定固体的摩尔级分， “总活度积”的两个值如下计算：??和aq? IAP1 ? IAP2solid（91） （92）??? x1?1K1 ? x2? 2 K2solid如果??? ?? aq ，那么包括固体溶液方程，否则不包括。如果固体溶液方程不包含在矩阵中，那么矩阵中未知数 dnp ss 所有的系数恒等于零。 在每次迭代中，对于固定压力气相，如果气相的摩尔数大于 1*10-14 或气相成分各部分 压力的和（根据水样的活度计算而来）大于总压力，则包括气相成分各部分压力和方程。如 果气相成分各部分压力和方程不包含在矩阵中，那么未知数 dNg 所有的系数恒等于零。f ptotal 、 f p 和 f p ss 作为优选方程包含在解答器中。所有其它方程作为解答器的相等约束条件。另外，也包括几个不等约束条件： （1）优选方程 f p 余项的值（等于）是非负值，并 作为无机物饱和或非饱和估计只值； （2） ）优选方程 f p ss 余项的值（等于 bpss ???jp ss , jxj ）是非负值，并作为固体溶液组分饱和或非饱和估计值； （3）优选方程 f ptotal 余项的值是非负 值，并作为总气体压力的非负估计值； （4）纯相质量的减少 dnp 小于或等于现存相的总摩尔 数 np； （5）固体溶液某种组分质量的减少 dnp ss 小于或等于存在组分的总摩尔数 np；和（6） 气相摩尔数的减少 dNgas 小于气相摩尔数 Ngas。 液相主要未知数的初始值为溶液组分最初的配比。如果涉及到两种或更多种溶液混合， 初始值为溶液的值加权它们的混合系数的和。 如果交换器或表面已经用一种溶液平衡， 初始 值即以前的平衡值。 如果它们没有用一种溶液平衡， 初始值与用于初始交换组分计算和初始 表面计算值相同。 纯相集合中每个相摩尔数的初始值、 固体溶液集合中每种成分摩尔数的初 始值、气相中每种气体成分摩尔数的初始值恒等于输入值或最后一次模拟中存入的值。 PHREEQC 的 数 据 输 入 ， 批 反 应 和 运 移 模 拟 计 算 依 靠 许 多 数 据 块 。 初 始 条 件 用 SOLUTION 或 SOLUTION-PPREAD 、 EXCHANGE 、 SURFACE 、 GAS-PHASE 、 EQUILIBRIUM CPHASES、SOLID-SOLUTIONS 和 USE 数据块确定。批反应用初始条件和 MIX、KINETICS、REACTION、REACTION-TEMPERATURE 和 USE 数据块确定。运移计 算用 ADVECTION 或 TRANSPORT 数据块确定（见“数据输入说明” ） 。化学动力学数值方法和速率表达式地球化学模型的不足之处是无机物、 有机物和其它反应物经常不能在一次实验或模拟期 间内反应到均衡。根据速率方程，均质或非均衡溶液的动态反应引起水样浓度的变化：dmi ? ci , k Rk ， dt（93）其中 ci , k 是动力反应中物质 i 的化学计量系数，Rk 是物质 k 的总反应速率（mol/kgw/s） 。一 般说来，数据块确定随反应程度变化，导致求解一常微分方程组。 对于不同的温度条件、压力和溶液组分，已经得到了许多动力反应速率。然而，不同研 究者使用不同的速率表达式以适合观察到的速率， 所以很难选择具有足够通用性的速率表达 式（一般很难编进程序） 。在 PHREEQC 版本中， 这个问题用根植于 BASIC 的翻译器得以 解决， 允许以一般的方式在输入文件中定义动力反应的速率表达式。 消除了在程序中硬行编 码速率表达式的需要。数值方法经过一段时间间隔，必须对速率求积分，这涉及到计算溶液中浓度的变化， ，同时，也 解释了其对反应速率的影响。 地球化学动力反应造成难解的方程组， 因为没有及时提取速率， 一些方程的速率（浓度变化对时间的微分）变化巨大，而另一些变化缓慢。 PHREEQC 用 Runge-Kutta（RK）算法解决了这个问题，求速率对时间的微分。Fehlberg（1969）的 RK 方 案，使用了达 6 个中间微分计算值。该方案包括用低阶 RK 方法求得一个误差估计值。用该 误差估计值和用户指定的容许误差比较，自动增加或减小时间间隔使误差在给定的范围内。 用户可在尝试最后的积分间隔前指定估计的中间 RK 的子间隔的数量（见数据输入说明） 。 该方案的系数取自 Cash 和 Karp（1990） 。速率表达式无机物和其它固体动力反应的总速率为：A ? mk Rk ? rk 0 ? V ? ? m0 k? ? ? ， ?n（94）其中 rk 是比速率（mol/m2/s） ，A0 为固体初始表面积（m2） ，V 为溶液的质量（kgw） ， m0 k 为 固体初始摩尔数， mk 为固体在一给定时间的摩尔数， ?mk / m0k ? 是一个系数，用于说明溶n解中 A0/V 的变化和固体的选择溶解和经过时间。对于统一的溶解 spheres andcubes n=2/3。 PHREEQC 中的所有计算用摩尔数，系数 A0/V 必须由用户指定得到合适的比率。 数据库的关键词 RATES 中包括一种物质的比速率表达式 rk。 这些比速率有不同的形式， 主要依赖经验资料的完善程度。资料缺乏时，经常使用的一个简单的速率表达式为：? ? IAP ?? ? ? ? rk ? kk ?1 ? ? ? ? ， ? ? K ? ? k ? ?（95） 其中 Kk 为一经验常数， IAP/ Kk 为饱和度 （SR） 。 这个速率方程可从 transition-state 理论得到， 当形成一活化络合物时（Aagaard 和 Helgeson，1982；Delany 和其它，1986） ，系数ζ 和反 应的化学计量有关。一般ζ =1。这个表达式的一个优点是该速率方程适用于超饱和和不饱 和两种情况，达到均衡时，速率为零。无论何时地球化学反应远离均衡时（IAP/ K&0.1） ，在 很长的时间范围内，速率为常数，并且 IAP/ K 趋近 1.0（均衡）时，速率趋近于零。 速率表达式也可以饱和指数为基础，采用下列形式：? IAP ? rk ? kk ? log? ? K ? ? 。 ? k ?该速率表达式已经成功用于大理石的溶解（Appelo 和其他，1984） 。 速率表达式经常包含与浓度有关的项。Monod 方程是其中之一：（96）? C ? rk ? rmax ? ?K ?C? ?， ? m ?（97）其中 rmax 是最大速率，Km 等于速率为最大值一半时的浓度。Monod 方程一般用于模拟有机 物氧化作用有序的步骤（Van Cappellen and Wang，1996） 。一系列速率表达式可由相应氧化 物的能量生成得到；首先消耗 O2，然后是 NO3-，接下来是其它，更慢反应的氧化物如 Fe2O3 和 SO4-2。Monod 方程的系数可由各个阶段的一级反应速率方程得到。土壤有机物降解的一 级反应速率方程为：dsC ? ? k1sC ， dt（98）其中 sC 是有机碳含量（mol/kg soil） ，K1 是一级衰变常数（s-1） 。在含氧环境下 K1 的值约等 于 0.025yr-1， 然而在荷兰砂质含水层中，NO3-是氧化物，k1≈5e-4 yr-1。在地下水甚至在有 机物氧化降解的氧化还原作用范围之外发现浓度高达 3μ Μ O2，而且可作为（依赖浓度的） 需氧降解的速率等于脱氮反应速率时的浓度。在 Monod 方程中利用系数 rmax= 1.57 e-9s-1 和 Km=294μ M 得到一级衰变反应（对于 0.3μ mO2，k1=0.025yr-1 和 3μ Μ O2，k1=5e-4 yr-1） ， 并且氧作为限制溶质。脱氮的类似估算是以 NO3-=3mM，k1=5e-4 yr-1 和 NO3-=3μ M，k1=1e-5 yr-为基础的，生成 rmax= 1.67 e-11s-1 和 Km=294μ M。在淡水含水层中，有机碳降解的综合 Monod 表达式为：?s RC ? 6sC ? C ? sC ? 01.67 *10?11 mNO ? ? ?? 1.57 *10?9 mO2 ? ? 3 ?? ? ? 4 ? 4 ?? 2.94 *10 ? mO 1.55*10 ? m ? ? NO3 ? 2 ?? ?（99）其中系数 6 是把土壤中的 sC 换算为空隙水中的 sC 时得到的。 有机物分解的另一重要方面是一部分很难溶而且特别不宜降解。已经提出一些模型解 释部分有机碳仍存在下去的趋势。 临时假定一系数 ?? sC ? sC ? 0? ? 组成总的速率二级反应。 该系数暗 ? ?示着浓度降为原始浓度的 1/10 将造成速率进一步降低为原来的/100。必须注意该系数用于 逼近一个很复杂的过程，并且需要严密分析，在 PHREEQC 版本中也可能给出确定速率的 其它灵活方法。 也有速率表达式是根据详细测量溶液中影响速率的浓度变化得出的。 例如 Willamson 和 Rimsidt（1994）给出的黄铁矿的速率表达式是：0.5 ?0.11 rpyrite ? 10?10.19 mO mH ? ， 2（100）式中含氧克分子浓度的平方根，并且说明了随 PH 值增大，速率有微小增加。该速率表达式 只用于溶解反应，且只用于溶液含氧时。当溶液接近均衡或氧被耗尽是，该速率表达式是不 适合的。 既适用于溶解又适用于沉淀的一个更完整的速率表达式是钙的速率方程：? ， rcalcite ? k1 H ? ? k2 ?CO2 ? ? k3 ?H2O? ? k4 Ca2? HCO3? ?????（101）其中括号内的化学符号表明为活度，系数 k1、k2 和 k3 已被 Plummer 等人（1978）在含带电 CO2 的钙溶解实验中确定为温度的函数。速率方程包括前半部分 rf（方程 101 的前三项）和 后半部分 rb（方程 101 的最后一项） ，并且既适用于溶解又适用于沉淀反应。后一部分包括 一个系数 k4，起值依赖于溶解组分。在单纯的方解石系统中，碳酸氢根离子浓度大约为钙离 子浓度的两倍，且后一部分速率表达式约等于：rb ? k4 Ca 2 ? HCO3? ? k4 2 Ca 2 ? 。均衡时，[Ca2+]为饱和时的活度[Ca2+]s。rcalcite=0，则??????2（102）2k4 ??Ca ?rf2? 2 s。（103）结合方程 101、102 和 103 得：2 ? ? 2? ? ? Ca ? ?。 rcalcite ? rf ?1 ? ? 2? 2 ? ? ? ? Ca s ? ? ? ?? ?? ?（104）CO2 压力为常数的纯系统 Ca-CO2 中，离子的活度积（IAP）为： IAPcalcite?Ca ??HCO ? ?2?? 2 3PCO2?Ca ? 和 K ?42? 3PCO 2calcite?Ca ? 。 ?42? 3 sPCO2（105）这样，对于一个方解石系统，方解石的速率近似为：2 ? ? ? IAP ? 3 ? ? ? rcalcite ? rf 1 ? ? ? ? Kcalcite ? ? ? ? ? ? ? ?（106）其中 rf 为方程 101 给出的前三项。数值方法和运移模拟方程PHREEQC 可 以 模 拟 几 个 一 维 运 移 过 程 ， 包 括 ：（ 1 ） 弥 散 ，（ 2 ） 平流（对流） ， （3）对流和弥散，和（4）对流和进入滞流区的弥散，这作为双空隙处理。所 有这些过程可和均衡及动力反应结合起来。平流-反应-弥散方程迁移化学物质（fig.1）的质量守恒生成对流-反应-弥散（ARD）方程：?C ?C ? 2C ?q ? ?v ? DL 2 ? ， ?t ?x ?x ?t（107）其中 C 是水中的浓度（mol/kgw） ，t 是时间（s） ，v 是水流速度（m/s） ，x 是距离（m） ，D L 是水力渗透系数[m2/s，DL = De+α Lv，De 是有效渗透系数，α L 是渗透率差（m）]，q 是固 相中的浓度（空隙中用 mol/kgw 表示） 。? v?C ? 2C 代表对流性运移， DL 代表弥散运移， ?x ?x 2?q 是由于反应在固相中浓度的变化（q 的单位同 C） 。一般假定 v 和 DL 等于所有的溶解物 ?t质，所以 C 可以是一种元素的总溶解浓度，包括所有的氧化还原物质。 图1对流-反应-弥散（ARD）方程中各项方程 107 的运移部分用先进行的、弥散中心的和对流运移迎风向的显有限差分法解答。 每种元素的化学交互作用项按每个时间段分别由运移部分计算， 并且是所有均衡和非均衡反 应速率的和。在 split-operator scheme（Press 和其他，1992；Yanenko，1971）中数值法随 ARD 方程中各项。对于每个时间步长，先计算对流运移，其次计算所有的均衡和动力控制 的化学反应， 再次计算弥散运移以及其所有的均衡和动力控制的化学反应。 这种方法不同于 其它多数的水文化学模型（Yeh 和 Tripathi，1989） ，因为动力和均衡化学反应是在对流计算 和弥散计算之后进行的。这减少了数值分散而且不必在物质组成和运移之间的重复计算。 split-operator scheme 的优点是通过对方程的各个部分的时间步长的调整，使其适合于 grid size（栅格尺寸） 。通过下列关系减小时间和距离离散之间的数值分散：??t ?A ? ?x ，v其中（108）是对流运移的时间步长（s） ，△x 是单元长度。弥散计算中数值的不稳定（波动）用下式估计：??t ?D ? ??x ?23DL，（109）其中 ??t ?D 是对流运移的时间步长（s） 。方程 108 和 109 的两个条件分别是对流运移计算的 Courant 和弥散运移计算的 Von Neumann 判据（例如，Press 和其他，1992） 。当 ?x ? ? L 时， 数值分散可忽略不计， 因为物质的对流运移和弥散运移一样重要甚至更重要。 当用一小栅格 间距减小数值分散时，对流运移计算的时间步长（方程 108）可能小于弥散运移计算的时间 步长（方程 108） ，因为前者包含 fine grid 的平方项用多个弥散时间步长解决这个矛盾，如? ??t ? ? ??t ?DA，并且在每个弥散时间步长计算化学反应。PHREEQC 的数据输入，一个弥散时间步长必须等于对流时间步长 ??t ?A ，或模拟弥散时等于弥散时间。此外，必须确定 变换次数，即计算的对流时间步长（弥散时间）的个数。 Central difference scheme 中的弥散运移必须是单元的混合。混合系数 mixf 定义为：m ixf ?DL ??t ?A ， 2 n??x ?（110）其中 n 为一正整数。限制条件为该单元和后面的不再混合，也就是说，mixf 必须小于 1/3， 如方程 109。当方程 110 中 n=1 时，mixf 大于 1/3，增加 n 的值使 mixf 小于或等于 1/3。弥 散时间为 ??t ?D ???t ?A ，并混合 n 次。n通过与简单线性交换的解析解相比较检查数值方法。 当交换不完全反应的交换系数等于 两个同价阳离子时形成线性交换结果。它给出了一个线性延迟值 R=1+CEC/C，其中 CEC 是 阳离子交换能力，用 mol/kgw 表示。在下列例子中，一个长 130m 的流水管中水的初始浓度 为 C（x，0）=Ci=0。排出的水浓度为 C=C0=1 mol/kgw，且空隙水流速为 v=15m/year。分散 率差为α L=5， 有效弥散系数为 De=0m2/s。 4 年后可给出两种化学成分的分布图， 一个为 R=1 （CL-） ，另一个为 R=2.5（Na+） 。 该问题需考虑两个边界条件。一个需要 x=0 时的流量或第三类边界条件：C ?0, t ? ? C0 ?DL ?C ? xend , t ? 。 v ?x（111）这个边界条件适用于输入管的横段面与实验管的横断面相比很小时。 ARD 方程的解 （Lindstrom 和其它，1967） ：C ? x, t ? ? Ci ?其中，DL=αL1 ?C0 ? Ci ?A 2v：（112）? x ? vt / R ? ? ?x ? vt / R ?2 ? x ?? A ? erfc? exp?? ?? ? 4? vt / R ? ?? 4 ? vt / R L L L ? ? ? ? ? x ? vt / R ? ? x ? 1? x vt / R ? ? ? ? ? ? ? 1 ? ? exp erfc ? ?? ? ? ? 2? ? ? 4 ? vt / R L L ? ? ? L? L ? ?。（113）插图 2 给出了具有不同（grid）栅格间距，△x=15、5 和 1.67m 时的三种模拟，相应的 和 和 和间距（m）图 2.具有离子交换和流通边界的一维运移解析解与 PHREEQC 计算在不同栅格间距时的比较时间分别为 ??t ?A =1、1/3 和 1/9 年。对于 CL- ，R=1，三种模拟的前部是不可区别的且与 解吸解完全一致。延迟离子 Na+ ，R=2.5，对于所有的（grid）栅格间距，漏过曲线的平均 位置是周围能正确的且与解析解一致。 最小的栅格间距结果与解析解最相近， 但需要最多的 计算时间。 计算时间主要依赖于 PHREEQC 中调用地球化学子程序的次数，没有动力反应时，调 用次数与单元个数*对流间隔个数*（1+弥散间隔个数）成比例。在这个例子中 DL= De+αLv=0+5*15m2/yr。这样，根据方程 110，对于逐渐减小的单元尺寸，mixf 分别等于 1/3、1 和 3。 对于 15m 的单元（mixf=1/3） ，对流步长取为弥散步长；对于 5m 的单元（mixf=1） ，三个对 流步长取为弥散步长；对于 1.67m 的单元（mixf=3） ，九个对流步长取为弥散步长。插图 2 说明了 CL（C/C0=0.5）四年旅行后对流的锋面曲线，当它达到 60m 时；对于 15m 的单元， 这需要 4 个对流步长。流管由 9 个单元组成，需对每个步长进行地球化学计算；因此，反应 计算次数为 9*4*（1+1）=72。大量的单元和对流步长适合于小 grid。另外两种情况反应计 算的调用次数为 27*12*（1+3）=1296；和 81*36*（1+9）=29160。 这儿给出的例子具有线性延迟使得能和解析解比较。然而线性延迟受到很大的数值分 散，而且在一定意义上，这个例子是数值分散糟糕的典型。在地球化学行业的许多方面，化 学反应帮助抵消数值分散，因为反应趋向于使波面曲度增大（sharpen fronts） ，如沉淀溶解 反应和沉降分离。 在其它例子中， 和不活泼的离子交换可能得到不受数值分散影响的真正的 物质弥散。 另一个边界条件是 Dirichlet 或者说是第一类边界条件， 它包括 x=0 时的常量浓度 C （0， t） ：C?0, t ? ? C0（114）当水从与下面土壤充分联系的蓄水池中渗漏时， 该边界掉件有效， 例如水从池塘或海水侵入 下覆沉积层。这种情况下 ARD 方程的解（Lapidus 和 Amundson，1952）为：C ? x, t ? ? Ci ?其中，1 ?C0 ? Ci ?B 2（115）? x ? vt / R ? ? ? ? ? ? ? exp? x ?erfc? x ? vt / R ? 。 A ? erfc? ? ? ? 4? vt / R ? ? 4? vt / R ? ? ?L ? L L ? ? ? ?（116）插图 3 说明了和以前运移例子相同的离散化模拟结果。 用三种栅格尺寸精确模拟守恒的 溶质（CL-，R=1） 。延迟物质（Na+，R=2.5）用大的网格说明了数值分散，但是，总的波面 是一致的。 在常浓度边界条件下， 因为指定的 x=0 处的条件弥散时间步长个数是流动时的两 倍。 第一类边界条件通过柱表面和外界溶液的接触会增大弥散。 通过边界的物质流动相应增 加，波锋比插图 2 中增加了几米高。例 11（临界点在柱的输出口）和例 12（弥散从恒定源 进行）给出了更多的分析溶液比较。 距离（m） 图 3. 具有离子交换和恒定边界条件的一维运移解析解与 PHREEQC 计算在不同栅格间距时的比较。热量传递由热量守恒的出热量传递方程， 或者说是温度变化。 一种物质的热量传递方程和其对流 -反应-弥散（ARD）方程是一致的：??wkw?T ?T ? 2T ? ?1 ? ???s k s ? ????wkw ?v ? ?L 2 ， ?t ?x ?x（117）其中 T 是温度（℃） ， ? 是空隙度（与总体积的比值，无单位） ，ρ 是密度（kg/m3） ，k 是比 热（kJ℃kg-1） ，κ 包含对流弥散和含水层的热传导（kJ℃m-1kg-1） ，下标 w 和 s 分别代表水 和固体。对于水和固体，温度 T 假定为一致的。 方程 117 两边同除以 ??wkw ：RT?T ?T ? 2T ? ?v ? ?L 2 ， ?t ?x ?x（118） 其中 RT ? 1 ??1 ? ???s ks 是温度延迟（无单位）系数， ???wkwL?? 热量耗散系数。该系数 ??wkw包括单纯弥散部分和由于对流耗散的部分： ?L ? ?e ? ?Lv ，与水力弥散系数相似。类似性}
我要回帖
更多关于 体重指数计算软件 的文章
更多推荐
- ·新注册的抖音号怎么才算养号成功会有播放量？
- ·抖音为什么零抖音发了半个小时播放量还是零？
- ·知乎2023年度盘点：如何通过“洞见对谈”展望留下一个你对未来的预言知乎？
- ·男子陌陌上警察会直接说约吗约美女学生酒店小聚被敲诈1.5万
- ·刘强东和奶茶生了几个为奶茶妹妹砸2000万注册公司，意欲何为？
- ·高清材质包会keyshot 怎么增加材质硬件要求吗
- ·IntelliJ IDEA 8ios 创建文件夹夹问题
- ·AIcohol Pad 遗精是什么原因造成的么
- ·如何和“小三卖子为男友”世界杯抢男友 网友爆笑支招
- ·messenger和kinky briefcasee的区别
- ·如何将torch.longtensor转换成intlong类型转换为int
- ·重庆php培训班哪家好多少钱
- ·windows10改密码 密码怎么改
- ·win7怎么win7无法删除打印机机驱动程序
- ·有人使用了APNs新的http 接口2接口了吗
- ·广东湖北交通运输管理局局网站
- ·ibatis 自定义sqlsql快定义的sql能在其他sql文件中引用吗
- ·大佬们 请问ps4千元ps4 pro显示器推荐有什么推荐吗
- ·如何用phreeqc软件来计算朗格利尔饱和指数数
- ·检测视频分辨率的真实分辨率
- ·海牛大大图标培训大数据好么？
- ·3.6里对arm64的编译器 支持 arm v7是怎么个情况
- ·s7s7总决赛小组赛赛b组fnc为什么打了8场
- ·中国女排吧江川张常宁现役国手颜值代表是张常宁吗 还是她
- ·芜湖羽毛接发一般多少钱，是步行街大万引进过来的吧
- ·光明轻巧包的机身，无论是手持还是夜间骑行都很方便
- ·体育产业发展存在问题教学及训练中存在哪些问题
- ·FIFAnba2konline配置要求3求一套最便宜好用的银卡套。各个位置都注明下多少钱。只求便宜好用，不想花人民币
- ·fifa fifaonline3好玩吗怎么换球队
- ·aj13伯爵爱与尊重有刮刮乐吗
- ·kpl王者荣耀耀kpl秋季赛 在哪里举办的最新相关信息
- ·18里球员2k18生涯模式能力值怎么提高球员能力值
