PHREEQC实例1

第四节 PHREEQC 一、PHREEQC的特点 PHREEQC的先前版本是PHREEQE。它可以模拟较宽范围的水岩反应，然而，PHREEQE有许多不足。 PHREEQC是在PHREEQE的源程序基础上用C语言重写而成的，它消除了PHREEQE的缺陷与限制。例如，PHREEQC可以定任何价态或价态组合而进行离子配分的摩尔平衡计算；也可以计算基于特定的pe或任意的氧化还原条件对各种价态间的氧化还原元素的分配；可以通过调整元素的浓度而达到与特定相态间的平衡，如达到特定的饱和指数或气体压力。溶液的浓度可以很容易地计算成各种所需的单位。 新版本增加的功能包括： （1）模拟一维溶质运移中的弥散（或扩散）和滞流区内的运移 （2）使用用户定义的速率表达式（或方程）模拟动力学反应，使模拟能完成与时间有关的动态化学反应； （3）模拟理想溶液、多组分溶液、非理想溶液及二元固体溶液的形成或溶解，同时也能计算固体间的化学反应； （4）除已有的固定压力气相反应，引入了固定体积的气相反应； （5）表面络合与离子交换模拟中加入了其交换与络合格点随溶解与沉淀的矿物或反应物的量改变而变化的功能； （6）在反向模拟计算中加入了同位素摩尔平衡，使反向模拟的解更加精确，同位素的计算可以用于校正地下水流动的时间或年龄； （7）自动使用多组收敛参数，用户自定义输出文件，以及多种输入文件格式，使输入与输出十分方便。 二、PHREEQC的局限性 （一）水溶液模型 PHREEQC用相关的离子和Debye Hückel方程计算非理想溶液，这种类型的水溶液模型在较低的离子强度下比较合适，但在高离子强度下（海水范围或更高）则会崩溃。 （二）离子交换 离子交换模型假设交换组分的热力学活度等于它的摩尔浓度，或者用摩尔浓度乘以Debye - Hückel 活度系数，以此定义成交换组分的活度（Appelo, 1994）。 （三）表面络合 PHREEQC加入了Dzombak 和 Morel （1990）的双层模型。其它模型，如三层、四层模型在PHREEQC中当前仍无法处理。表面络合位置数目、表面积、被吸附相成分、合适的logK等的确定也是不确定的，实际模拟时需要用实际的野外介质进行实验室研究，通过获取的实验数据来确定合适的表面络合模型，加入PHREEQC中进行模拟。 （四）固溶体 PHREEQC使用Guggenheim方法来确定非理想二元相固溶体。尚无法计算非理想的三元固溶体。 （五）溶质运移模拟 模拟一维对流弥散或滞流区中的扩散时PHREEQC使用的是显式有限差分算法。当网格划分较粗时，此算法出现数值弥散。推荐的做法是逐步模拟。首先用粗网格尽快获得结果并研究水化学反应过程，最后用细划网格评价数值弥散对反应组分和结果组分的影响。 三、PHREEQC计算实例 例1、矿物相的溶解平衡模拟 这个例子定义了最稳定相石膏或是硬石膏在一定温度范围的溶解度变化。输入的数据集见表9-1。定义纯水溶液时仅用pH和温度两个参数。默认单位为毫摩尔（millimolal）。缺省状态下，pe为4.0，氧化还原计算缺省使用pe值，水的密度为1.0。在多组反应期间，所有列于EQUILIBRIUM_PHASES中平衡相在批反应计算中不管其存在与否，都允许达到指定的饱和指数。相的输入按相名称、期望达到的饱和指数和当前该相矿物的总量，用摩尔数表示。如果一种矿物在开始时不存在，那么其总量为0.0 mol。该例中，石膏和硬石膏允许反应达到平衡（饱和指数0.0指矿物与处于平衡状态），矿物相初始量都为1.0 mol。每一种矿物或是反应达到平衡，或者在矿物分配时消耗完。在大多数的情况下，1mol的单矿物相对反应平衡计算而言是足够的。 表9-1 例1溶解平衡模拟的数据输入格式 SOLUTION 1 Pure water temp 25.0 pH 7.0 EQUILIBRIUM_PHASES 1 Anhydrite 0.0 1.0 Gypsum 0.0 1.0 REACTION_TEMPERATURE 1 25.0 75.0 in 51 steps SELECTED_OUTPUT file ex2.sel si anhydrite gypsum END 在REACTION_TEMPERATURE数据块中，指定温度计算步长为1℃，开始于25℃，结束于75℃，共计算温度变化为51次。计算结果见图9-1所示。 图9-1 在25到75℃的温度范围内溶液中石膏和硬石膏的饱和指数（据David et al., 1999） 由图9-1可以看出，在初始温度25℃时计算两矿物相的平衡饱和指数是硬石膏处于不饱和状态，其饱和指数为-0.22，说明其已溶解完毕。而石膏则已达到饱和状态（饱和指数为0.0），因而石膏在此种状态时处于向沉淀进行的趋势。计算结果显示，硬石膏已转化为石膏沉淀达0.985摩尔。由于石膏沉淀带走了溶液中的水，于是溶液中的水不再为1kg。随温度升高，硬石膏饱和指数升高，而石膏不变化，直至温度达到57℃后，石膏变为不稳定相，饱和指数小于0.0，溶解完毕，而硬石膏则处于过饱和，开始沉淀。 例9-2、不可逆反应模拟——黄铁矿氧化 这个例子反映了PHREEQC模拟不可逆反应的能力，以黄铁矿氧化反应为例。氧气（O2）和NaCl被不可逆地加入到纯水之中，加入的量为五个不同浓度（0.0, 1.0, 5.0, 10.0, 50.0 mmol）；在不可逆反应中O2 和NaCl的相对比例分别为1.0和0.5。黄铁矿、方解石和针铁矿允许溶解以达到平衡，二氧化碳分压力保持在10-3.5（大气分压力）。另外，在石膏达到过饱和时，允许其沉淀。   表9-2 例9-2数据输入文件 TITLE EXAMPLE 2. Add oxygen, equilibrate with pyrite, calcite, and geothite SOLUTION 1 PURE WATER pH 7.0 temp 25.0 EQUILIBRIUM_PHASES 1 Pyrite 0.0 Goethite 0.0 Calcite 0.0 CO2(g) -3.5 Gypsum 0.0 0.0 REACTION 1 O2 1.0 NaCl 0.5 0.0 0.001 0.005 0.01 0.05 SELECTED_OUTPUT file ex5.sel totals Cl si Gypsum equilibrium_phases pyrite goethite calcite CO2(g) gypsum END 表9－3 例9-2的结果选择如下 [摩尔转移与相分配的摩尔数有关；正值表示当前相数量的增加，也即沉淀；负值表示当前相数量的减少，也即溶解] 反应物加入量 毫摩尔 pH pe 摩尔转移（毫摩尔） 石膏的饱和指数 O2 NaCl 黄铁矿 针铁矿 方解石 CO2(g) 石膏 0.0 0.0 8.28 -4.94 -0.000032 0.000011 -0.49 -0.49 0.0 -6.13 1.0 0.5 8.17 -4.29 -0.27 0.27 -0.93 0.14 0.0 -2.02 5.0 2.5 7.98 -3.97 -1.33 1.33 -2.94 2.40 0.0 -1.06 10.0 5.0 7.88 -3.82 -2.67 2.67 -5.56 5.11 0.0 -.65 50.0 25.0 7.72 -3.57 -13.33 13.33 -26.84 26.49 9.00 0.0 表9－2中定义了纯水溶液，平衡相以EQUILIBRIUM_PHASES标识，有黄铁矿，针铁矿，方解石和二氧化碳，缺省量均为10mol，而只是石膏为0.0mol。石膏项的定义指示石膏如达到了过饱和，它将会沉淀；而在初始状态下不可能有溶解的可能，因为其初始值为0.0。反应标识REACTION定义参加不可逆反应物质和数量，该例子中，氧（“O2”）以1.0的相对系数加入，而NaCl则会以0.5的相对系数加入。反应量分别定义为0.0, 0.001, 0.005, 0.01, 和0.05 mol。在每一步平衡计算中SELECTED_OUTPUT用来保存氯化物的总浓度，石膏的饱和指数，黄铁矿、针铁矿、二氧化碳和石膏的总数量和转移的摩尔数，文件名为ex5.sel。 例9-2中的部分结果列于表9-3中。当没有氧和氯化钠加到体系中时，很少量的方解石和二氧化碳溶解，以及微量的黄铁矿和针铁矿发生了反应；由于与黄铁矿达到平衡，pH为8.28，pe很低（-4.94），方解石低于饱和状态的6个数量级（饱和指数为-6.13）。当加入氧和氯化钠，黄铁矿被氧化，针铁矿由于相对难溶解而沉淀。反应中产生了硫酸， pH值降低了，pe轻微上升，导致方解石溶解及二氧化碳的析出。在加入10和50mmol氧之间的一些点上，方解石达到了饱和，开始沉淀。当加入了50mmol的氧和25mmol的氯化钠后，总共有9.00mmol的总石膏沉淀出来了。 例9-3、反应路径的模拟 PHREEQC中有三种方法可用于求解反应路径问题：求解反应路径和相边界的交点；不断增加反应物质量来求解反应路径；反应路径按动力学过程求解。在第一种方法中，不必知道反应数量，但需要一系列的模拟以便寻找到相关相边界。在第二种方法中，只需一次模拟便足够了，但必须先知道反应合适量。在第三种方法中，求解反应路径需要动力学反应速率表达方程式，运用时间步长来调整算法。 在下面这个例子中，钾长石被放进烧杯中，慢慢地反应。当钾长石溶解时，其他矿物相便可能沉淀。该例中，假定了仅能形成水铝矿、高岭石或是钾云母，如果这些相达到饱和，那么这些相将发生沉淀。在反应开始时沉淀的矿物相，随着反应的进行也可再次被溶解。下面的模拟用于求解水铝矿、高岭石、钾云母与钾长石系统中的反应路径。 表9-4 例9-3的数据输入 TITLE Example 3A.--React to phase boundaries. SOLUTION 1 PURE WATER pH 7.0 charge temp 25.0 PHASES Gibbsite Al(OH)3 + 3 H+ = Al+3 + 3 H2O log_k 8.049 delta_h -22.792 kcal Kaolinite Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3 log_k 5.708 delta_h -35.306 kcal K-mica KAl3Si3O10(OH)2 + 10H+ = 3Al+3 + 3H4SiO4 + K+ log_k 12.970 delta_h -59.377 kcal K-feldspar KAlSi3O8 + 4 H2O + 4 H+ = Al+3 + 3 H4SiO4 + K+ log_k 0.875 delta_h -12.467 kcal SELECTED_OUTPUT -file ex3A-B.sel -activities K+ H+ H4SiO4 -si Gibbsite Kaolinite K-mica K-feldspar -equilibrium Gibbsite Kaolinite K-mica K-feldspar END TITLE Example 3A1.--Find amount of K-feldspar dissolved to reach gibbsite saturation. USE solution 1 EQUILIBRIUM_PHASES 1 Gibbsite 0.0 KAlSi3O8 10.0 Kaolinite 0.0 0.0 K-mica 0.0 0.0 K-feldspar 0.0 0.0 END TITLE Example 3A2.--Find amount of K-feldspar dissolved to reach kaolinite saturation. USE solution 1 EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 KAlSi3O8 10.0 K-mica 0.0 0.0 K-feldspar 0.0 0.0 END TITLE Example 3A3.--Find amount of K-feldspar dissolved to reach K-mica saturation. USE solution 1 EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 0.0 K-mica 0.0 KAlSi3O8 10.0 K-feldspar 0.0 0.0 END TITLE Example 3A4.--Find amount of K-feldspar dissolved to reach K-feldspar saturation. USE solution 1 EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 0.0 K-mica 0.0 0.0 K-feldspar 0.0 KAlSi3O8 10.0 END TITLE Example 3A5.--Find point with kaolinite present,but no gibbsite. USE solution 1 EQUILIBRIUM_PHASES 1 Gibbsite 0.0 KAlSi3O8 10.0 Kaolinite 0.0 1.0 END TITLE Example 3A6.--Find point with K-mica present,but no kaolinite USE solution 1 EQUILIBRIUM_PHASES 1 Kaolinite 0.0 KAlSi3O8 10.0 K-mica 0.0 1.0 END 表9-4定义的数据，SOLUTION数据块定义的是纯水，其中只需温度和pH值；PHASES数据块中定义的是水铝矿、高岭石、含钾云母、钾长石的溶解化学方程及其化学热力学数据； SELECTED_OUTPUT指定数据的输出是钾离子、氢离子和硅酸活度系数的对数；水铝矿、高岭石、钾云母、钾长石的饱和指数及在矿物相分配过程中的转移总量，其计算结果保存于表9－5之中。EQUILIBRIUM_PHASES数据块则定义矿物平衡参数，如3A1中第一行定义10.0 mol KAlSi3O8（钾长石的分子式）溶解反应足以使水铝矿达到饱和状态，而高岭石、钾云母和钾长石在此反应中初始量为0.0，饱和指数为0.0，此设置指示这三种矿物达到饱和时可以沉淀，但开始时不可能溶解。3A2~3A6定义的平衡相数据块与3A1类似。   表9-5 例9-3选择性输出 [输入设置模拟的编号。负摩尔转移指示溶解，正摩尔转移指示沉淀。 模拟编号 钾长石的摩尔转移,微摩尔 活度对数 摩尔转移，微摩尔 饱和指数 图形上的点 H+ H4SiO4 水铝矿 高岭石 钾云母 水铝矿 高岭石 钾 云 母 钾 长 石 3A1 -0.03 -7.01 -0.57 -7.10 0.00 0.00 0.00 0.0 -3.8 -10.7 -14.7 A 3A2 -2.18 -8.21 2.55 -5.20 1.78 .00 .00 .0 .0 -1.9 -5.9 B 3A3 -20.02 -9.11 4.41 -4.47 .00 9.71 .00 -.7 .0 .0 -2.5 D 3A4 -190.9 -9.39 5.49 -3.55 .00 .00 63.61 -2.0 -.7 .0 .0 F 3A5 -3.02 -8.35 2.83 -5.20 .00 1.24 .00 .0 .0 -1.6 -5.6 C 3A6 -32.68 -9.07 4.41 -4.25 .00 .00 10.78 -.9 .0 .0 -2.1 E 将模拟结果投点于这些相的热动力学稳定场中，得A、B、C、D、E点，如图9-2所示。通过A点与表9-5可以看出，A点属于水铝矿的稳定存在区，溶液中水铝矿已达到饱和，开始沉淀。而高岭石、钾云母和钾长石均不饱和。在B点水铝矿即将开始转化成高岭石，这一转化过程是一个消耗Si的过程；但在B点，水铝矿与高岭石同处于饱状态，若要转化为高岭石，该反应路径需沿着该两相边界到达一个C点，然后进入高岭石相区，即生成区，在此区域中，只有高岭石沉淀，其它几类矿物均处于溶解状态下；反应路径穿过高岭石稳定区后到达D点，此点位于钾云母与高岭石两相的边界上，此时从表9-5可以看出，钾云母已达到饱和状态，但要沉淀尚需沿两相边界达到一个E点后才会进入钾云母的生成区，直至达F点。F点是钾长石与钾云母的饱和点。在模拟3A5中，计算允许钾长石溶解到达C点，该点高岭石呈饱和状态且存在于物相分配之中，但此时当水铝矿虽处在饱和状态，并不出现物相的分配，即转移的摩尔数为0.0（见表9-5）。同样，模拟3A6求解出E点，此点钾云母为饱和态且参与物相分配，此时高岭石虽达到饱和，但未沉淀，即转移的摩尔数为0.0（见表9-5）。在3A5和3A6中，均假定高岭石和钾云母有1mol的初始量，该量足以达到矿物平衡状态。 图9-2 25℃下，纯水中钾长石（微斜长石）溶解相图（据David et al., 1999） 该例通过模拟来确定反应路径，方法简单。只需通过增加钾长石反应量即可同时允许水铝矿，高岭石，钾云母和钾长石物相间发生相态转变。同时计算得出从点A到点F钾长石的溶解量介于0.03到190.0mmol范围。 例9-4、表面络合作用模拟 PHREEQC的表面络合模型有三种：（1）默认的是一般性的双层模型，扩散层组分的计算使用非显式算法。（2） 静电双层模型，使用显式计算扩散层组分，通过设置-diffuse_layer项，可选用此模型。（3）第三种是非静电模型，通过设置-no_edl来选用。静电双层模型是在一般性双层模型的基础上修改的模型，修改包括：（1）表面可能有两种以上的成键位置；（2）不包括表面沉淀，（3）可使用一种代替分子式，该分子式可表示电荷与势能之间的关系，显式地计算扩散层的组分，此种形式用-diffuse_layer来设置。非静电模型不考虑表面电荷对表面络合作用形成的影响，结果是表面络合的数学处理方式如同水溶液中的络合反应，只是没有活度系数项。 本例中使用一般性双层模型、扩散层非显式计算。在氢氧化铁表面，假定有强和弱两种成键类型位置，模拟锌在其表面的吸附。H离子和锌离子竞争两种成键位置，用质量作用方程描述该过程。表面组分的活度大小取决于由表面电荷的变化引起的表面电势。该例中，考虑锌吸附的变化是pH的函数，在0.1M硝酸钠电解液中锌浓度从低（10-7 M）到高（10-4 M）变化。 表9-6 例9-4数据输入 SURFACE_SPECIES Hfo_sOH + H+ = Hfo_sOH2+ log_k 7.18 Hfo_sOH = Hfo_sO- + H+ log_k -8.82 Hfo_sOH + Zn+2 = Hfo_sOZn+ + H+ log_k 0.66 Hfo_wOH + H+ = Hfo_wOH2+ log_k 7.18 Hfo_wOH = Hfo_wO- + H+ log_k -8.82 Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+ log_k -2.32 SURFACE 1 Hfo_sOH 5e-6 600. 0.09 Hfo_wOH 2e-4 SOLUTION 1 units mmol/kgw pH 8.0 Zn 0.0001 Na 100. charge N(5) 100. SOLUTION 2 units mmol/kgw pH 8.0 Zn 0.1 Na 100. charge N(5) 100. USE solution none # # Model dedinitions # PHASES Fix_H+ H+ = H+ log_k 0.0 END # # Zn = 1e-7 # SELECTED_OUTPUT file ex8.sel molalities Zn+2 Hfo_wOZn+ Hfo_sOZn+ USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.0 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.25 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.5 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.75 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.0 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.25 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.5 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.75 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.0 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.25 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.5 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.75 NaOH 10.0 END USE solution 1 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -8.0 NaOH 10.0 END # # Zn = 1e-4 # USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.0 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.25 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.5 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -5.75 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.0 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.25 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 2 Fix_H+ -6.5 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -6.75 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.0 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.25 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.5 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -7.75 NaOH 10.0 END USE solution 2 USE surface 1 EQUILIBRIUM_PHASES 1 Fix_H+ -8.0 NaOH 10.0 END 表9-6中，使用了三个关键字：SURFACE_MASTER_SPECIES，SURFACE_SPECIES和SURFACE。在缺省状态，SURFACE_MASTER_SPECIES数据块定义了两个成键位置，名字为“Hfo” “Hfo_w”和“Hfo_s”，分别表示强和弱成键位置。 图9-3 总锌浓度为10-7和10-4摩尔浓度水合氧化铁的强和弱表面位置作为pH值的函数时，氧化锌在液相中的分布。（据David et al., 1999） 模拟的初始条件是定义两个不同锌浓度的硝酸钠溶液(SOLUTION 1 和2 数据块)。“Fix_H+”是一个假设的物相，用于定义固定的pH条件，其pH通过加入不同量NaOH来调节。模拟目的是研究两种不同的锌溶液中在pH从5至8的变化范围条件下，表面络合反应的情况。 模拟的结果见图9-3。图9-3显示模拟结果与Dzombak 和 Morel（1990）中结果一致。锌在高pH条件下强烈被吸附，低pH条件下则吸附较弱。此外，锌浓度较低时，在整个pH值范围内，强成健位置对锌的络合作用造成的吸附量均大于弱成键位置的吸附量，且pH值较高时，大部分锌被强成键位置吸附。在锌的浓度较大的情况下，且仅在低pH值时，强成键位置吸附占据主导，高pH值时，由于所有的强成键位置都占满，因此在pH值较高，锌的浓度较大时，更多的锌被弱成键位吸附。 例9-5 溶解铁离子与氧的动力学氧化反应 动力学化学反应过程的模拟最重要的是描述动力学的参数与方程。PHREEQC中处理这个问题使用一种非常简单、通用的形式。即使用速率关键字RATES来描述有关的参数与方程，在这个由RATES描述的数据块中用一种简单Basic语言来定义一个与时间相关速率算法。在模拟过程，程序不断地调用这个算法以获取所需时段的物质反应量。同时这个速率算法块也可以用在KINETICS数据块中进行批反应或运移计算。模拟溶质运移（对流或运移）过程的动力学反应时，可用关键字KINETICS来定义一个一个的动力学反应。动力学速率表达式中结合一个4级和5级的Runge-Kutta-Fehlberg算法，平衡反应在动力学计算被初始化之前计算一次，在加入动力学反应增量之后再计算一次。平衡反应的计算包括所有已经定义的反应，包括溶液中离子配分、各种离子交换、平衡相、固溶液、表面配分与气相等。 本例模拟计算Fe2+被氧化成Fe3+的的动力学反应。数据输入见表9-7。 在水中的Fe2+被O2氧化的反应速率可写成为： （9-1） 式中：t —— 时间（s）； aOH- —— 氢氧根离子的活度； mFe2+ —— 溶液中二价铁的总摩尔浓度； Po2 —— 氧气的分压力（atm）。 表9-7. 例5数据输入 SOLUTION_MASTER_SPECIES Fe_di Fe_di+2 0.0 Fe_di 55.847 Fe_tri Fe_tri+3 0.0 Fe_tri 55.847 SOLUTION_SPECIES Fe_di+2 = Fe_di+2 log_k 0.0 Fe_tri+3 = Fe_tri+3 log_k 0.0 # # Fe+2 species # Fe_di+2 + H2O = Fe_diOH+ + H+ log_k -9.5 delta_h 13.20 kcal # #... and also other Fe+2 species # Fe_di+2 + Cl- = Fe_diCl+ log_k 0.14 Fe_di+2 + CO3-2 = Fe_diCO3 log_k 4.38 Fe_di+2 + HCO3- = Fe_diHCO3+ log_k 2.0 Fe_di+2 + SO4-2 = Fe_diSO4 log_k 2.25 delta_h 3.230 kcal Fe_di+2 + HSO4- = Fe_diHSO4+ log_k 1.08 Fe_di+2 + 2HS- = Fe_di(HS)2 log_k 8.95 Fe_di+2 + 3HS- = Fe_di(HS)3- log_k 10.987 Fe_di+2 + HPO4-2 = Fe_diHPO4 log_k 3.6 Fe_di+2 + H2PO4- = Fe_diH2PO4+ log_k 2.7 Fe_di+2 + F- = Fe_diF+ log_k 1.0 # # Fe+3 species # Fe_tri+3 + H2O = Fe_triOH+2 + H+ log_k -2.19 delta_h 10.4 kcal # #... and also other Fe+3 species # Fe_tri+3 + 2 H2O = Fe_tri(OH)2+ + 2 H+ log_k -5.67 delta_h 17.1 kcal Fe_tri+3 + 3 H2O = Fe_tri(OH)3 + 3 H+ log_k -12.56 delta_h 24.8 kcal Fe_tri+3 + 4 H2O = Fe_tri(OH)4- + 4 H+ log_k -21.6 delta_h 31.9 kcal 2 Fe_tri+3 + 2 H2O = Fe_tri2(OH)2+4 + 2 H+ log_k -2.95 delta_h 13.5 kcal 3 Fe_tri+3 + 4 H2O = Fe_tri3(OH)4+5 + 4 H+ log_k -6.3 delta_h 14.3 kcal Fe_tri+3 + Cl- = Fe_triCl+2 log_k 1.48 delta_h 5.6 kcal Fe_tri+3 + 2 Cl- = Fe_triCl2+ log_k 2.13 Fe_tri+3 + 3 Cl- = Fe_triCl3 log_k 1.13 Fe_tri+3 + SO4-2 = Fe_triSO4+ log_k 4.04 delta_h 3.91 kcal Fe_tri+3 + HSO4- = Fe_triHSO4+2 log_k 2.48 Fe_tri+3 + 2 SO4-2 = Fe_tri(SO4)2- log_k 5.38 delta_h 4.60 kcal Fe_tri+3 + HPO4-2 = Fe_triHPO4+ log_k 5.43 delta_h 5.76 kcal Fe_tri+3 + H2PO4- = Fe_triH2PO4+2 log_k 5.43 Fe_tri+3 + F- = Fe_triF+2 log_k 6.2 delta_h 2.7 kcal Fe_tri+3 + 2 F- = Fe_triF2+ log_k 10.8 delta_h 4.8 kcal Fe_tri+3 + 3 F- = Fe_triF3 log_k 14.0 delta_h 5.4 kcal PHASES Goethite Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O log_k -1.0 END SOLUTION 1 pH 7.0 pe 10.0 O2(g) -0.67 Fe_di 0.1 Na 10. Cl 10. charge EQUILIBRIUM_PHASES 1 O2(g) -0.67 RATES Fe_di_ox -start 10 Fe_di = TOT("Fe_di") 20 if (Fe_di <= 0) then goto 200 30 p_o2 = 10^(SI("O2(g)")) 40 moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME 200 SAVE moles -end KINETICS 1 Fe_di_ox -formula Fe_di -1.0 Fe_tri 1.0 -steps 100 400 3100 10800 21600 5.04e4 8.64e4 1.728e5 1.728e5 1.728e5 1.728e5 INCREMENTAL_REACTIONS true SELECTED_OUTPUT -file ex9.sel -reset false USER_PUNCH -headings Days Fe(2) Fe(3) pH si_goethite 10 PUNCH SIM_TIME/3600/24 TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite") END 该例中定义溶液（用SOLUTION关键字）为氯化钠溶液，其中有0.1mmol/kgw二价铁离子（Fe-di）的，并与大气中的氧达到平衡。EQUILIBRIUM_PHASES关键字定义所有的批次反应都将与大气中的氧达到平衡；这样，对于铁离子的氧化而言提供了充足的氧气供应。 在RATES数据块中，速率表达式指定名字为“Fe_di_ox”的氧化反应速率按上面式（1）计算，一共用了5行Basic语言来描述。 在KINETICS数据块中，名为Fe-di_ox的速率表达式被调用，并定义了有关参数。在KINETICS数据块中，如速率的名字与PHASES下的矿物的名字相同时，则矿物分子表达式用于化学反应计算。 KINETICS数据块-steps标识符给定动力学反应积分的时间步长。当设置INCREMENTAL_REACTIONS为true时，每一次时间步长将会加到所模拟的总的时间上进行模拟，前一时段计算的结果作为当前时间的计算的起点来使用。 图9－4为上述数据的模拟结果。图9-4反映出在10天模拟反应期间，反应容器中总Fe(II)，总Fe(III)的浓度，及pH的变化。在反应开始时，能够看到pH值迅速地降低。Fe(II)随时间下降得非常快，但随着反应的进行，趋于稳定，这与上面的式（1）相一致。当在实际情况的非缓冲溶液中进行实验时，会发现pH最初是升高的。pH的这种升高与Fe(III)的氢氧络合物缓慢形成有关。因为氧化反应的消耗氢离子，pH开始时是升高的。 PHREEQC中这种定义反应速率的通用方法，使任何动力学反应的速率表达都能实现，换言之，PHREEQC中可以计算任何动力学过程。 图9-4 总Fe(II)，总Fe(III)的浓度和pH随亚铁离子[Fe(II)]氧化成三价铁离子[Fe(III)]的变化。（据David et al., 1999）