phreeqc实例练习

《phreeqc实例练习》由会员分享，可在线阅读，更多相关《phreeqc实例练习（102页珍藏版）》请在金锄头文库上搜索。

1、PHREEQCPHREEQC 实例分析实例分析例例 1 1物种形成分析物种形成分析这个例子计算了海水中矿物质的分布以与一组有关矿物在海水中的饱和程度。 为了证明如何在这个模型中应用新的元素，将元素铀添参加由phreeqc.dat定义的液相模型中wateq.dat是包含于程序分类中的一个数据库文件， 它来自于 WATEQ4F Ball and Nordstrom,1991 ，并包含铀。物质形成计算所需要的数据包括温度、Ph、元素的浓度和/或其元素的化合价。海水中的这些数据见表10。这个例子计算中输入的数据组见表11。在模拟中所运用的有关计算的注释包含在TITLETITLE关键字中。SOLUTIO

2、NSOLUTION数据块定义了海水的成分。注意：元素的化合价用元素化学符号后面圆括号中的数字表示S(6),N(5), N(-3)和O(0)。表表 1010海水的成分海水的成分未指定浓度时，其浓度的单位为 ppm分析的组分分析的组分钙镁钠钾铁锰硅石，SiO2氯化物-碱度，HCO32-硫酸盐，SO4-硝酸盐，NO3+铵，NH4铀pH，标准单位pe，无单位温度，密度，千克/升PHREEQCPHREEQC 符号符号CaMgNaKFeMnSiClAlkalinityS(6)N(5)N(-3)UpHpetemperaturedensity浓度浓度412.31291.810768.0399.1.002.00

3、024.2819353.0141.6822712.0.29.03.00338.228.45125.01.023用于分配氧化复原元素和计算饱和指数的pe由redoxredox标识符所指定。在这个例子中，用氧化复原电对O(-2)/O(0) 计算的pe值相对应于溶解氧/水， 并且这个pe适用于需要pe值的所有的计算。如果redoxredox没有指定，那么缺省的值将会是所输入的pe。缺省的氧化复原标识符可被任何氧化复原元素代替，如输入元素锰时，那么输入的 pe被用来表示各1 / 102种化合价状态的锰；输入铀时，这里是氮/铵电对将会用来计算所形成各种价态铀的pe值。数据组中缺省的单位为ppmunits

4、units标识符。这个缺省值可以替换为任何浓度单位，如指定铀的浓度为ppb来代替ppm。因为ppm是一个质量单位，而不是一个摩尔单位，这个程序必须用分子量来将浓度单位转化为摩尔单位。每一种主要物质缺省的分子量在SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES输入中指定 缺省数据库phreeqc.dat的值列在表4和附录B中 。如果提交的分子量数据不同于其缺省值， 必须在输入数据的设置中指定适当的分子量。 这可以用gfwgfw标识符来完成，在这里输入真正的分子量，转化硝酸盐的分子量为62.0 g/mol，或是更简便的是以asas标识符来完成，在这里输入

5、所使用的化学分子式的单位，正如在这个例子中输入的碱和铵是一样的。注意最后给定的溶解氧 O(0)的浓度是1ppm的初始估计值，但它的浓度将会得以调整，直到氧气分压的对数达到-0.7。O2(g)的定义是在缺省数据库文件中在PHASESPHASES输入(附录B)。当使用相均衡来指定初始浓度正如这个例子中的O(0)，那么仅有一种浓度是得以调整。例如，例如石膏被用来调整钙的浓度，钙的浓度会改变，而硫酸盐的浓度却保持不变。表 11例 1 的输入数据TITLE Example 1.-Add uranium and speciate seawater.SOLUTION 1 SEAWATER FROM NORD

6、STROM ET AL. (1979)units ppmpH 8.22pe 8.451density 1.023temp 25.0redox O(0)/O(-2)Ca 412.3Mg 1291.8Na 10768.0K 399.1Fe0.002Mn 0.0002 peSi 4.28Cl 19353.0Alkalinity 141.682 as HCO3S(6) 2712.0N(5) 0.29 gfw 62.0N(-3) 0.03 as NH4U3.3 ppb N(5)/N(-3)O(0) 1.0 O2(g) -0.7SOLUTION_MASTER_SPECIESU U+4 0.0 238.02

7、90 238.0290U(4) U+4 0.0 238.02902 / 102U(5) UO2+ 0.0 238.0290U(6) UO2+2 0.0 238.0290SOLUTION_SPECIES#primary master species for U#is also secondary master species for U(4)U+4 = U+4log_k 0.0U+4 + 4 H2O = U(OH)4 + 4 H+log_k -8.538delta_h 24.760 kcalU+4 + 5 H2O = U(OH)5- + 5 H+log_k -13.147delta_h 27.5

8、80 kcal#secondary master species for U(5)U+4 + 2 H2O = UO2+ + 4 H+ + e-log_k -6.432delta_h 31.130 kcal#secondary master species for U(6)U+4 + 2 H2O = UO2+2 + 4 H+ + 2 e-log_k -9.217delta_h 34.430 kcalUO2+2 + H2O = UO2OH+ + H+log_k -5.782delta_h11.015 kcal2UO2+2 + 2H2O = (UO2)2(OH)2+2 + 2H+log_k -5.6

9、26delta_h -36.04 kcal3UO2+2 + 5H2O = (UO2)3(OH)5+ + 5H+log_k-15.641delta_h -44.27 kcalUO2+2 + CO3-2 = UO2CO3log_k 10.064delta_h 0.84 kcalUO2+2 + 2CO3-2 = UO2(CO3)2-2log_k 16.977delta_h 3.48 kcalUO2+2 + 3CO3-2 = UO2(CO3)3-4log_k 21.397delta_h -8.78 kcalPHASESUraniniteUO2 + 4 H+ = U+4 + 2 H2Olog_k -3.

10、490delta_h -18.630 kcal3 / 102END程序的的数据库文件phreeqc.dat中不包含铀。这样，当应用这个数据库文件时，输入文件中一定得包括描述热动力学和液相中含铀组分的数据。 需要两个关键字来定义铀的形态，即SOLUTION_MASTER_SPECIES和SOLUTION_SPECIES。通过把这两个数据块加到输入文件中，将会在程序运行中确定液相中含铀组分。为把铀稳定地加到列出的元素中，那么这些数据块应参加到数据库文件中。这里铀的数据是说明性的，而不是铀物质的完整描述。使用SOLUTION_MASTER_SPECIES输入来定义含铀的主要物质成分是必要的。因为铀是

11、活泼的氧化复原元素，所以定义具有不同化合价的次要含铀物质也是很有必要的。SOLUTION_MASTER_SPECIES(表11)数据块定义了U+4为主要的含铀物质，同时+4价的铀也是次级主要物质。 UO2是化合价为+5的次级主要含铀物质， UO2是化合价为+6的次级主要含铀物质。 定义这些液相和其它任何铀络合物的方程必须通过SOLUTION_SPECIES输入来进展。在数据块SOLUTION_SPECIES (表11)中，主要的和次要的物质均附有注释。首要的主要物质总是以恒等反响U+4 = U+4的形式来定义的。次主要物质是在化学反响中仅有的含有电子的液相。另外的氢氧化物和碳酸盐络合物定义为+

12、4和+6价，无+5价。最后，在PHASES输入中定义一种新的含铀矿物。在物质形成模拟中该物质将会被用来计算饱和指数， 在计算机运行中的批反响、 运移或是反向模拟中， 如果没有重新定义，那么不能使用。表 12-例 1 的输出Input file: ex1Output file: ex1.outDatabase file: ./phreeqc.dat-Reading data base.-SOLUTION_MASTER_SPECIESSOLUTION_SPECIESPHASESEXCHANGE_MASTER_SPECIESEXCHANGE_SPECIESSURFACE_MASTER_SPECIES

13、SURFACE_SPECIESRATESEND-Reading input data for simulation 1.-4 / 102+2SOLUTION 1 SEAWATER FROM NORDSTROM ET AL. (1979) units ppm pH 8.22 pe 8.451 density 1.023 temp 25.0 redox O(0)/O(-2) Ca 412.3 Mg 1291.8 Na 10768.0 K 399.1 Fe 0.002 Mn 0.0002 pe Si 4.28 Cl 19353.0 Alkalinity 141.682 as HCO3 S(6) 27

14、12.0 N(5) 0.29 as NO3 N(-3) 0.03 as NH4 U 3.3 ppb N(5)/N(-3) O(0) 1.0 O2(g) -0.7SOLUTION_MASTER_SPECIES U U+4 0.0 238.0290 238.0290 U(4) U+4 0.0 238.0290 U(5) UO2+ 0.0 238.0290 U(6) UO2+2 0.0 238.0290SOLUTION_SPECIES U+4 = U+4 log_k 0.0 U+4 + 4 H2O = U(OH)4 + 4 H+ log_k -8.538 delta_h 24.760 kcal U+

15、4 + 5 H2O = U(OH)5- + 5 H+ log_k -13.147 delta_h 27.580 kcal U+4 + 2 H2O = UO2+ + 4 H+ + e- log_k -6.432 delta_h 31.130 kcal U+4 + 2 H2O = UO2+2 + 4 H+ + 2 e- log_k -9.217 delta_h 34.430 kcal UO2+2 + H2O = UO2OH+ + H+ log_k -5.7825 / 102 delta_h 11.015 kcal 2UO2+2 + 2H2O = (UO2)2(OH)2+2 + 2H+ log_k

16、-5.626 delta_h -36.04 kcal 3UO2+2 + 5H2O = (UO2)3(OH)5+ + 5H+ log_k -15.641 delta_h -44.27 kcal UO2+2 + CO3-2 = UO2CO3 log_k 10.064 delta_h 0.84 kcal UO2+2 + 2CO3-2 = UO2(CO3)2-2 log_k 16.977 delta_h 3.48 kcal UO2+2 + 3CO3-2 = UO2(CO3)3-4 log_k 21.397 delta_h -8.78 kcalPHASES Uraninite UO2 + 4 H+ =

17、U+4 + 2 H2O log_k -3.490 delta_h -18.630 kcalEND-TITLE- Example 1.-Add uranium and speciate seawater.-Beginning of initial solution calculations.-Initial solution 1. SEAWATER FROM NORDSTROM ET AL. (1979)-Solution composition-Elements Molality MolesAlkalinity 2.406e-03 2.406e-03Ca 1.066e-02 1.066e-02

18、Cl 5.657e-01 5.657e-01Fe 3.711e-08 3.711e-08K 1.058e-02 1.058e-02Mg 5.507e-02 5.507e-026 / 102Mn 3.773e-09 3.773e-09N(-3) 1.724e-06 1.724e-06N(5) 4.847e-06 4.847e-06Na 4.854e-01 4.854e-01O(0) 3.746e-04 3.746e-04 Equilibrium with O2(g)S(6) 2.926e-02 2.926e-02Si 7.382e-05 7.382e-05U 1.437e-08 1.437e-0

19、8-Description of solution- pH = 8.220 pe = 8.451 Activity of water = 0.981 Ionic strength = 6.748e-01 Mass of water (kg) = 1.000e+00 Total carbon (mol/kg) = 2.180e-03 Total CO2 (mol/kg) = 2.180e-03 Temperature (deg C) = 25.000 Electrical balance (eq) = 7.936e-04 Percent error, 100*(Cat-|An|)/(Cat+|A

20、n|) = 0.07 Iterations = 7 Total H = 1.110147e+02 Total O = 5.563047e+01-Redox couples-Redox couple pe Eh (volts)N(-3)/N(5) 4.6737 0.2765O(-2)/O(0) 12.3893 0.7329-Distribution of species- Log Log LogSpecies Molality Activity Molality Activity GammaOH- 2.674e-06 1.629e-06 -5.573 -5.788 -0.215H+ 7.981e

21、-09 6.026e-09 -8.098 -8.220 -0.122H2O 5.551e+01 9.806e-01 -0.009 -0.009 0.000C(4) 2.180e-03HCO3- 1.514e-03 1.023e-03 -2.820 -2.990 -0.170MgHCO3+ 2.195e-04 1.640e-04 -3.658 -3.785 -0.127NaHCO3 1.667e-04 1.948e-04 -3.778 -3.710 0.0677 / 102MgCO3 8.913e-05 1.041e-04 -4.050 -3.982 0.067NaCO3- 6.718e-05

22、5.020e-05 -4.173 -4.299 -0.127CaHCO3+ 4.597e-05 3.106e-05 -4.337 -4.508 -0.170CO3-2 3.821e-05 7.959e-06 -4.418 -5.099 -0.681CaCO3 2.725e-05 3.183e-05 -4.565 -4.497 0.067CO2 1.210e-05 1.413e-05 -4.917 -4.850 0.067UO2(CO3)3-4 1.255e-08 1.183e-10 -7.901 -9.927 -2.025UO2(CO3)2-2 1.814e-09 5.653e-10 -8.7

23、41 -9.248 -0.506MnCO3 2.696e-10 3.150e-10 -9.569 -9.502 0.067MnHCO3+ 6.077e-11 4.541e-11 -10.216 -10.343 -0.127UO2CO3 7.429e-12 8.678e-12 -11.129 -11.062 0.067FeCO3 1.952e-20 2.281e-20 -19.709 -19.642 0.067FeHCO3+ 1.635e-20 1.222e-20 -19.786 -19.913 -0.127Ca 1.066e-02Ca+2 9.504e-03 2.380e-03 -2.022

24、-2.623 -0.601CaSO4 1.083e-03 1.265e-03 -2.965 -2.898 0.067CaHCO3+ 4.597e-05 3.106e-05 -4.337 -4.508 -0.170CaCO3 2.725e-05 3.183e-05 -4.565 -4.497 0.067CaOH+ 8.604e-08 6.429e-08 -7.065 -7.192 -0.127Cl 5.657e-01Cl- 5.657e-01 3.528e-01 -0.247 -0.452 -0.205MnCl+ 9.582e-10 7.160e-10 -9.019 -9.145 -0.127M

25、nCl2 9.439e-11 1.103e-10 -10.025 -9.958 0.067MnCl3- 1.434e-11 1.071e-11 -10.844 -10.970 -0.127FeCl+2 9.557e-19 2.978e-19 -18.020 -18.526 -0.506FeCl2+ 6.281e-19 4.693e-19 -18.202 -18.329 -0.127FeCl+ 7.786e-20 5.817e-20 -19.109 -19.235 -0.127FeCl3 1.417e-20 1.656e-20 -19.849 -19.781 0.067Fe(2) 6.909e-

26、19Fe+2 5.205e-19 1.195e-19 -18.284 -18.923 -0.639FeCl+ 7.786e-20 5.817e-20 -19.109 -19.235 -0.127FeSO4 4.845e-20 5.660e-20 -19.315 -19.247 0.067FeCO3 1.952e-20 2.281e-20 -19.709 -19.642 0.067FeHCO3+ 1.635e-20 1.222e-20 -19.786 -19.913 -0.127FeOH+ 8.227e-21 6.147e-21 -20.085 -20.211 -0.127FeHSO4+ 3.0

27、00e-27 2.242e-27 -26.523 -26.649 -0.127Fe(3) 3.711e-08Fe(OH)3 2.841e-08 3.318e-08 -7.547 -7.479 0.067Fe(OH)4- 6.591e-09 4.924e-09 -8.181 -8.308 -0.127Fe(OH)2+ 2.118e-09 1.583e-09 -8.674 -8.801 -0.127FeOH+2 9.425e-14 2.937e-14 -13.026 -13.532 -0.506FeSO4+ 1.093e-18 8.167e-19 -17.961 -18.088 -0.127FeC

28、l+2 9.557e-19 2.978e-19 -18.020 -18.526 -0.506FeCl2+ 6.281e-19 4.693e-19 -18.202 -18.329 -0.1278 / 102Fe+3 3.509e-19 2.796e-20 -18.455 -19.554 -1.099Fe(SO4)2- 6.371e-20 4.760e-20 -19.196 -19.322 -0.127FeCl3 1.417e-20 1.656e-20 -19.849 -19.781 0.067Fe2(OH)2+4 2.462e-24 2.322e-26 -23.609 -25.634 -2.02

29、5FeHSO4+2 4.228e-26 1.318e-26 -25.374 -25.880 -0.506Fe3(OH)4+5 1.122e-29 7.679e-33 -28.950 -32.115 -3.165H(0) 0.000e+00H2 0.000e+00 0.000e+00 -44.436 -44.369 0.067K 1.058e-02K+ 1.041e-02 6.494e-03 -1.982 -2.187 -0.205KSO4- 1.639e-04 1.225e-04 -3.785 -3.912 -0.127KOH 3.137e-09 3.664e-09 -8.504 -8.436

30、 0.067Mg 5.507e-02Mg+2 4.742e-02 1.371e-02 -1.324 -1.863 -0.539MgSO4 7.330e-03 8.562e-03 -2.135 -2.067 0.067MgHCO3+ 2.195e-04 1.640e-04 -3.658 -3.785 -0.127MgCO3 8.913e-05 1.041e-04 -4.050 -3.982 0.067MgOH+ 1.084e-05 8.100e-06 -4.965 -5.092 -0.127Mn(2) 3.773e-09Mn+2 2.171e-09 4.982e-10 -8.663 -9.303

31、 -0.639MnCl+ 9.582e-10 7.160e-10 -9.019 -9.145 -0.127MnCO3 2.696e-10 3.150e-10 -9.569 -9.502 0.067MnSO4 2.021e-10 2.360e-10 -9.695 -9.627 0.067MnCl2 9.439e-11 1.103e-10 -10.025 -9.958 0.067MnHCO3+ 6.077e-11 4.541e-11 -10.216 -10.343 -0.127MnCl3- 1.434e-11 1.071e-11 -10.844 -10.970 -0.127MnOH+ 2.789e

32、-12 2.084e-12 -11.555 -11.681 -0.127Mn(NO3)2 1.375e-20 1.606e-20 -19.862 -19.794 0.067Mn(3) 5.993e-26Mn+3 5.993e-26 4.349e-27 -25.222 -26.362 -1.139N(-3) 1.724e-06NH4+ 1.648e-06 9.272e-07 -5.783 -6.033 -0.250NH3 7.507e-08 8.769e-08 -7.125 -7.057 0.067N(5) 4.847e-06NO3- 4.847e-06 2.846e-06 -5.315 -5.

33、546 -0.231Mn(NO3)2 1.375e-20 1.606e-20 -19.862 -19.794 0.067Na 4.854e-01Na+ 4.791e-01 3.387e-01 -0.320 -0.470 -0.151NaSO4- 6.053e-03 4.522e-03 -2.218 -2.345 -0.127NaHCO3 1.667e-04 1.948e-04 -3.778 -3.710 0.067NaCO3- 6.718e-05 5.020e-05 -4.173 -4.299 -0.127NaOH 3.117e-07 3.641e-07 -6.506 -6.439 0.067

34、O(0) 3.746e-04O2 1.873e-04 2.188e-04 -3.727 -3.660 0.0679 / 102S(6) 2.926e-02SO4-2 1.463e-02 2.664e-03 -1.835 -2.574 -0.740MgSO4 7.330e-03 8.562e-03 -2.135 -2.067 0.067NaSO4- 6.053e-03 4.522e-03 -2.218 -2.345 -0.127CaSO4 1.083e-03 1.265e-03 -2.965 -2.898 0.067KSO4- 1.639e-04 1.225e-04 -3.785 -3.912

35、-0.127HSO4- 2.089e-09 1.561e-09 -8.680 -8.807 -0.127MnSO4 2.021e-10 2.360e-10 -9.695 -9.627 0.067FeSO4+ 1.093e-18 8.167e-19 -17.961 -18.088 -0.127Fe(SO4)2- 6.371e-20 4.760e-20 -19.196 -19.322 -0.127FeSO4 4.845e-20 5.660e-20 -19.315 -19.247 0.067FeHSO4+2 4.228e-26 1.318e-26 -25.374 -25.880 -0.506FeHS

36、O4+ 3.000e-27 2.242e-27 -26.523 -26.649 -0.127Si 7.382e-05H4SiO4 7.110e-05 8.306e-05 -4.148 -4.081 0.067H3SiO4- 2.720e-06 2.032e-06 -5.565 -5.692 -0.127H2SiO4-2 7.362e-11 2.294e-11 -10.133 -10.639 -0.506U(4) 1.040e-21U(OH)5- 1.040e-21 7.773e-22 -20.983 -21.109 -0.127U(OH)4 1.662e-25 1.941e-25 -24.77

37、9 -24.712 0.067U+4 0.000e+00 0.000e+00 -46.994 -49.020 -2.025U(5) 1.627e-18UO2+ 1.627e-18 1.216e-18 -17.789 -17.915 -0.127U(6) 1.437e-08UO2(CO3)3-4 1.255e-08 1.183e-10 -7.901 -9.927 -2.025UO2(CO3)2-2 1.814e-09 5.653e-10 -8.741 -9.248 -0.506UO2CO3 7.429e-12 8.678e-12 -11.129 -11.062 0.067UO2OH+ 3.386

38、e-14 2.530e-14 -13.470 -13.597 -0.127UO2+2 3.019e-16 9.410e-17 -15.520 -16.026 -0.506(UO2)2(OH)2+2 1.780e-21 5.547e-22 -20.750 -21.256 -0.506(UO2)3(OH)5+ 2.908e-23 2.173e-23 -22.536 -22.663 -0.127-Saturation indices-Phase SI log IAP log KTAnhydrite -0.84 -5.20 -4.36 CaSO4Aragonite 0.61 -7.72 -8.34 C

39、aCO3Calcite 0.76 -7.72 -8.48 CaCO3Chalcedony -0.51 -4.06 -3.55 SiO2Chrysotile 3.36 35.56 32.20 Mg3Si2O5(OH)4CO2(g) -3.38 -21.53 -18.15 CO2Dolomite 2.41 -14.68 -17.09 CaMg(CO3)2Fe(OH)3(a) 0.19 -3.42 -3.61 Fe(OH)310 / 102Goethite 6.09 -3.41 -9.50 FeOOHGypsum -0.63 -5.21 -4.58 CaSO4:2H2OH2(g) -41.22 1.

40、82 43.04 H2Hausmannite 1.57 19.56 17.99 Mn3O4Hematite 14.20 -6.81 -21.01 Fe2O3Jarosite-K -7.52 -42.23 -34.71 KFe3(SO4)2(OH)6Manganite 2.39 6.21 3.82 MnOOHMelanterite -19.35 -21.56 -2.21 FeSO4:7H2OO2(g) -0.70 -3.66 -2.96 O2Pyrochroite -8.08 7.12 15.20 Mn(OH)2Pyrolusite 6.95 5.29 -1.66 MnO2:H2OQuartz

41、-0.11 -4.06 -3.96 SiO2Rhodochrosite -3.27 -14.40 -11.13 MnCO3Sepiolite 1.16 16.92 15.76 Mg2Si3O7.5OH:3H2OSepiolite(d) -1.74 16.92 18.66 Mg2Si3O7.5OH:3H2OSiderite -13.13 -24.02 -10.89 FeCO3SiO2(a) -1.35 -4.06 -2.71 SiO2Talc 6.04 27.44 21.40 Mg3Si4O10(OH)2Uraninite -12.67 4.40 17.06 UO2-End of simulat

42、ion.-模拟中的输出表12包含标题所描绘的几个信息块。首先，是运行的输入、输出、数据库文件的名字。其次，在标题“Reading data base以下出了在读数据库中碰到的所有关键字。后面，输入数据在标题“Reading input data for simulation 1下进展重复输出，不包括注释和空行。所有的输入数据以与END关键字构成了该模拟。在这个模拟中，TITLE关键字后面所碰到的任何注释都将打印在后面。名称后面是标题，“Beginning of initial solution calculations ，在它的下面是海水物质形成计算的结果。浓度数据，转化成重量摩尔数，在子标题

43、“Solutioncomposition下所给出。对初始溶液计算而言，溶液中的摩尔数目在数字上等于它的重量摩尔数，这是因为假定的是1kg的水。标识符-water是用来定义溶液中不同质量的水。在批反响计算中，水的质量可发生变化，液相中的摩尔数不会准确地等于组分的重量摩尔数。注意，产生分压力对数为-0.7的溶解氧的重量摩尔数已被计算，并且在输出中会加以解释。在子标题“Description of solution之后，在输出的第一个数据块中所列出的一些属性等于它们的输入值，另一些是计算值。在这个例子中， pH，pe和温度都等于它们的输入值。离子强度，总碳碱度是输入数据，总无机碳“ Total CO

44、2，电子平衡，和百分误差在模拟计算中得到。在子标题“Redox couples下面打印的是可获得的每个氧化复原电对的pe和Eh，11 / 102在该例中，为铵/硝酸盐，以与水/溶解氧。在子标题“Distribution of species下面，列出的是每种元素所有形态和价态的摩尔浓度、活度与活度系数。它的顺序是根据元素的字母顺序或是按每种元素的浓度或是元素化合价递减的顺序列出的。 除了每种元素或元素的化合价， 也给出了总摩尔数。最后，在子标题“Saturation indices下，适于给出分析数据的矿物的饱和指数，在输出的末尾局部以物质的名字按子母顺序列出。 饱和指数是在标有“SI的栏中给

45、出的，跟在后面的是的一栏是离子活度积的对数 “log IAP和溶解常数的对数“logKT。每种物质的化学分子式都在右边栏中打出。注意，例如没有包括含铝的矿物，这是因为在分析数据中不包括铝。同样也注意，输出中也没有包括四方硫铁矿FeS和其它的硫化矿物，这是因为没有指定S-2的分析数据。如输入了S代替S6或S-2的浓度，那么这个S-2的浓度将会被加以计算，四方硫铁矿和其它硫化矿物的饱和指数也将得到计算。例例 2 2矿物相的溶解平衡矿物相的溶解平衡这个例子测定了最稳定相石膏或是无水石膏在一定温度围的溶解性。 输入数据组见表 13。仅用 pH 和温度来定义纯水溶液。缺省单位是millimolal，但没

46、有指定浓度。缺省状态下，pe 为 4.0，缺省的氧化复原计算使用pe，且水的密度为1.0这是没有必要的，因为没有浓度为“每升 。在批反响期间，所有允许反响达到指定的饱和指数的反响都列在EQUILIBRIUM_PHASES 中，无论开始时它们是否存在。输入数据包含相的名字之前通过PHASES 输入在数据库或输入文件中定义 、指定的饱和指数和以摩尔数表示的当前存在的相的数量。如果一种相在当前不存在，那么在纯相集合中给定其浓度为0.0 mol。在这个例子中，石膏和无水石膏允许反响来达到平衡饱和指数等于 0.0 ，初始相的集合中每种矿物都为 1mol。每一种矿物或是反响达到平衡或是被耗尽。在大多数的情

47、况下， 1 mol 的单纯相足够达到反响平衡。表表13-13-例例2 2输入数据的设置输入数据的设置TITLE Example 2.-Temperature dependence of solubility of gypsum and anhydrite SOLUTION 1 Pure water pH 7.0 temp 25.0 EQUILIBRIUM_PHASES 1 Gypsum 0.0 1.0 Anhydrite 0.0 1.0 REACTION_TEMPERATURE 112 / 102 25.0 75.0 in 51 steps SELECTED_OUTPUT file ex2.s

48、el si anhydrite gypsumEND在 REACTION_TEMPERATURE数据块中， 指定反响温度计算步长为1， 从 25开始， 75完毕，共计算温度变化 51 次。温度的每一度变化都详细指定输入数据，在可能的情况下，在 EQUILIBRIUM_PHASES石膏和无水石膏中定义的物质将会反响达到平衡，或是直到两种物质都完全溶解。 最后在每一次计算之后， 用 SELECTED_OUTPUT 来输出石膏和无水石膏的饱和指数到文件 ex2.sel 中。然后用这个文件来生成图5。初始溶液计算和第一次批反响的结果见表14。纯水中物质的分布列于标题 “开始初始溶液计算下。在 25下，与

49、所给定石膏和无水石膏量达到平衡的体系是第一次批反响，反响结果列在标题“开始批反响计算之后。 紧接着这个标题，定义了批反响阶段，接下来是计算中确定关键字数据的列表。在这个例子中，计算时所应用的反响组分保存为序号1，纯相集合保存为序号 1，反响温度保存为序号 1。概念上，溶液和纯相放在一个烧杯中，并调整温度为 25，使其反响达到体系平衡。图 5-在 25 到 75的温度围，溶液中石膏和硬石膏稳定时的饱和指数，在子标题“Phase assemblage下，列出了由 EQUILIBRIUM_PHASES 所定义的每一种相的饱和指数和数量。在第一次的批反响阶段， 最后相的集合不含无水石膏， 对溶液而言这

50、并没有达到平衡 饱和指数为-0.22 ， 1.985mol 的石膏， 与溶液达到平衡 饱和指数等于 0.0 。所有的无水石膏得以溶解， 且大多数的钙和硫酸盐是生成石膏沉淀下来。“溶液组分说明，15.64mmol/kgw 的钙和硫酸盐留在溶液中，这确定了纯水中石膏的溶解度。然而，在液相中每一成分的摩尔总数为 15.08，这是因为水的质量仅为0.9645kg “溶液的定义 。在形成石膏CaSO4.2H2O时，水从溶液中移除。因此，在批反响计算中，溶剂水的质量并不是常13 / 102数；在溶解和沉淀相中，反响和水合作用可以增加或减少溶剂水的质量。所有批反响阶段的饱和指数绘于图5。在每一阶段，纯水与相

51、在不同的温度下反响反响并不累加 。PHREEQC 缺省的数据块说明：在温度低于 57时，石膏是稳定相饱和指数等于 0.0 ；在高于这个温度时，无水石膏是稳定相。表 14例 2 的选择性输出-Beginning of initial solution calculations.-Initial solution 1.Pure water-Solutioncomposition-Elements Molality MolesPure water-Descriptionsolution- pH = 7.000 pe = 4.000 Activity of water = 1.000 Ionic st

52、rength = 1.001e-07 Mass of water (kg) = 1.000e+00 Total alkalinity (eq/kg) = 1.082e-10 Total carbon (mol/kg) = 0.000e+00 Total CO2 (mol/kg) = 0.000e+00 Temperature (deg C) = 25.000 Electrical balance (eq) = -1.082e-10 Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.05 Iterations = 0 Total H = 1.110124

53、e+02 Total O = 5.550622e+01-Distributionspecies- Log Log LogSpecies Molality Activity Molality Activity GammaOH- 1.002e-07 1.001e-07 -6.999 -6.999 -0.00014 / 102ofofH+ 1.001e-07 1.000e-07 -7.000 -7.000 -0.000H2O 5.551e+01 1.000e+00 0.000 0.000 0.000H(0) 1.416e-25H2 7.079e-26 7.079e-26 -25.150 -25.15

54、0 0.000O(0) 0.000e+00O2 0.000e+00 0.000e+00 -42.080 -42.080 0.000-Saturationindices-Phase SI log IAP log KTH2(g) -22.00 -22.00 0.00 H2O2(g) -39.12 44.00 83.12 O2-Beginning of batch-reaction calculations.-Reaction step 1.Using solution 1.Pure waterUsing pure phase assemblage 1.Using temperature 1.-Ph

55、aseassemblage- Moles in assemblagePhase SI log IAP log KT Initial Final DeltaAnhydrite -0.22 -4.58 -4.36 1.000e+00 -1.000e+00Gypsum 0.00 -4.58 -4.58 1.000e+00 1.985e+00 9.849e-01-Solutioncomposition-Elements Molality MolesCa 1.564e-02 1.508e-02S 1.564e-02 1.508e-02-Descriptionofsolution-15 / 102 pH

56、= 7.062 Charge balance pe=10.691Adjusted toredox equilibrium Activity of water = 1.000 Ionic strength = 4.178e-02 Mass of water (kg) = 9.645e-01 Total alkalinity (eq/kg) = 1.122e-10 Total carbon (mol/kg) = 0.000e+00 Total CO2 (mol/kg) = 0.000e+00 Temperature (deg C) = 25.000 Electrical balance (eq)

57、= -9.766e-11 Percent error, 100*(Cat-|An|)/(Cat+|An|) = -0.00 Iterations = 23 Total H = 1.070728e+02 Total O = 5.359671e+01-Distributionspecies- Log Log LogSpecies Molality Activity Molality Activity GammaOH- 1.402e-07 1.155e-07 -6.853 -6.937 -0.084H+ 1.006e-07 8.665e-08 -6.997 -7.062 -0.065H2O 5.55

58、1e+01 9.996e-01 -0.000 -0.000 0.000Ca 1.564e-02Ca+2 1.045e-02 5.176e-03 -1.981 -2.286 -0.305CaSO4 5.191e-03 5.242e-03 -2.285 -2.281 0.004CaOH+ 1.192e-08 9.910e-09 -7.924 -8.004 -0.080H(0) 4.363e-39H2 2.181e-39 2.202e-39 -38.661 -38.657 0.004O(0) 1.701e-15O2 8.504e-16 8.587e-16 -15.070 -15.066 0.004S

59、(-2) 0.000e+00HS- 0.000e+00 0.000e+00 -117.650 -117.734 -0.084H2S 0.000e+00 0.000e+00 -117.807 -117.803 0.004S-2 0.000e+00 0.000e+00 -123.278 -123.590 -0.312S(6) 1.564e-02SO4-2 1.045e-02 5.075e-03 -1.981 -2.295 -0.313CaSO4 5.191e-03 5.242e-03 -2.285 -2.281 0.004HSO4- 5.141e-08 4.276e-08 -7.289 -7.36

60、9 -0.080-Saturation16 / 102ofindices-Phase SI log IAP log KTAnhydrite -0.22 -4.58 -4.36 CaSO4Gypsum 0.00 -4.58 -4.58 CaSO4:2H2OH2(g) -35.51 -35.51 0.00 H2H2S(g) -116.81 -158.45 -41.64 H2SO2(g) -12.11 71.01 83.12 O2Sulfur -87.18 -122.94 -35.76 S例例 3 3混合反响混合反响这个例子证明了 PHREEQC 进展一系列地球化学模拟的能力， 在同一次运行中， 最

61、后的模拟依赖于以前模拟的结果。这个例子研究了海水与碳酸盐地下水混合区域中发生的成岩反-2.0响。该例子分五次模拟，在表 15 中标为 A 到 E。 A在 PCO2为 10atm 时，通过纯水与方解石的溶解平衡定义碳酸盐地下水。 B用在表 10 中给定的主要离子数据定义海水。 C将 70%的地下水与 30%的海水溶液进展混合。 D混合物与方解石和白云石达到平衡。 最后，E假定白云石沉淀可以忽略的情况下，研究混合物仅与方解石达到平衡的化学演化。例 3. 输入的数据组TITLE Example 3, part A.-Calcite equilibrium at log Pco2 = -2.0 and

62、 25.SOLUTION 1 Pure water pH 7.0 temp 25.0EQUILIBRIUM_PHASES CO2(g) -2.0 Calcite 0.0SAVE solution 1ENDTITLE Example 3, part B.-Definition of seawater.SOLUTION 2 Seawater units ppm pH 8.22 pe 8.451 density 1.023 temp 25.0 Ca 412.3 Mg 1291.8 Na 10768.0 K 399.1 Si 4.28 Cl 19353.0 Alkalinity 141.682 as

63、HCO317 / 102 S(6) 2712.0ENDTITLE Example 3, part C.-Mix 70% ground water, 30% seawater.MIX 1 1 0.7 2 0.3SAVE solution 3ENDTITLE Example 3, part D.-Equilibrate mixture with calcite and dolomite.EQUILIBRIUM_PHASES 1 Calcite 0.0 Dolomite 0.0USE solution 3ENDA 局部 表 15 的输入包括在 SOLUTION 输入中纯水的定义， 以与在 EQUIL

64、IBRIUM_PHASES输入中纯相组分的定义。在定义这些相时，仅给定每一种相的饱和指数。如果没有指定，那每一种相的缺省量为 10.0mol，对大多数相而言，这根本上是无限的。暗中定义的批反响是一个平衡，即在模拟中第一次定义的溶液与在模拟中第一次定义的纯相的集合所达到的平衡。 批反响的显性定义是用 USE 关键字来完成的 。SAVE 关键字命令这个程序保存最终的批反响中形成的批反响溶液组分为溶液1。这样，当模拟开始时，溶液1 是纯水。在模拟的批反响计算完成后，批反响溶液水与方解石和CO2达到平衡保存为溶液1。B 局部定义了海水的成分， 保存为溶液 2。C 局部是地下水溶液1与海水溶液2在封闭体

65、系中相混合，其中 Pco2应计算得到，而不指定。MIXMIX 关键字是用来定义混合物中每种溶液所占的比例大致的混合体积 。SAVESAVE 关键字将混合物保存为溶液 3。MIXMIX 关键字允许以指定的任何比例混合。比例分数体积不必是1.0。如果用7.0 和 3.0 来代替 0.7 和0.3，那么溶液 1 中每一种元素包括氢和氧摩尔数的数目将乘以7.0，在溶液 2 中每一种元素摩尔数的数目将乘以3.0，元素的总摩尔数是两者的和。在混合物中水的质量将大约为 10kg 来自溶液 1 的为 7.0， 来自溶液 2 的为 3.0 ， 代替比例分别为 0.7 和 0.3 时的 1kg。混合物的浓度会与任

66、一混合比例相一致，因为溶液1 和溶液 2 的相比照例是一样的。然而，在后面的反响中， 对于 7.0 和 3.0 的混合分数， 将会消耗比表 16 中所给定数值大 10 倍的摩尔转换，这是因为在体系中具有10 倍的水。D 局部是与方解石和白云石的平衡。USEUSE 关键字指定了达到平衡后的溶液为溶液 3，这是来自 C 局部的混合物。通过使用“EQUILIBRIUM_PHASES1来定义相的集合，这个相的集合将代替了前面在 A 局部所定义的集合 1。E 局部进展与 D 局部相似的计算，但使用相集合2，集合 2 中白云石不作为反响物。表表 16.16. 例例 3 3 的选择输出结果的选择输出结果模拟

67、 A 产生了碳酸盐地下水；B 定义了海水；C 进展没有其它摩尔转移的混合；D 平衡了方解石和白云石的混合物；E 仅仅平衡了方解石的混合物。摩尔转移与相集合中的摩尔数有关；正数说明当前相数量的增18 / 102长，也就是说沉淀；负数说明了当前相数量的降低，即溶解。饱和指数：“-说明饱和指数的计算是不可行的，这是由于组成元素的一种在溶液中不存在。摩尔转移： “-说明这种矿物在模拟中不允许进展摩尔转换模拟ABCDEpH7.2978.2207.3517.0567.442logPCO2-2.00-3.38-2.23-1.98-2.31饱和指数方解石0.00.76-.10.00.00白云石-2.41.52

68、.00.73CO2-1.976-摩尔转移方解石-1.646-15.71-.040白云石-7.935-在表 16 中所列出了例 3 的选择输出结果。由 A 局部的地下水与方解石达到平衡，且在输入中指定二氧化碳分压的对数为-2.0。在相集合中的 CO2摩尔数会降低大约为 2.0mmol，这意味着大约有 2.0mmol CO2溶解到溶液中。同样地，大约有 1.6mmol 的方解石被溶解。B局部定义了海水，它的计算值比大气中的二氧化碳要大一些 -3.38 与-3.5 比照 ，与方解石饱和指数为0.76和白云石2.41将会达到过饱和。在B 局部不允许有矿物的转移。C 局部进展混合，并计算了混合物中各种物

69、质的平衡分布，同样不允许有矿物的摩尔转移。二氧化碳分压的对数为-2.23，方解石没有达到饱和，白云石达到了过饱和。该饱和指数说明热动力学状态下白云石化作用应该发生，也即是方解石应该溶解而白云石应当沉淀。 D 局部计算了应当发生反响的方解石和白云石的量。为了达到平衡， 15.71mmol 的方解石应当溶解，7.935mmol 的白云石应当沉淀。在当前区域环境中，没有观测到白云石化作用，即使白云石是热动力学稳定相。无白云石化作用是由于白云石形成时缓慢的反响动力学。因此，E局部模拟了如果白云石不沉淀时将会发生的反响。如果白云石没有沉淀，对于该混合比率，将仅一小局部的方解石会溶解0.040mmol 。

70、例例 4 4蒸发和类似的氧化复原反响蒸发和类似的氧化复原反响蒸发是通过除去化学系统中的水完成的， 水能够以三种方法除去： 1在REACTION关键字输入中， 可指定水为具有负反响系数的不可逆反响物， 2在MIX关键字输入中，将溶液与给定负混合分数的纯水相混合，或3在EQUILIBRIUMPHASES关键字输入中，指定“H2O为可替代的反响，在这种情况下，从液相中增减水以达到纯相的指定饱和指数。这个例子应用的是第一种方法； REACTION数据块用来模拟大约20倍雨水的浓度，它是通过移除95%的水实现的。产生的结果溶液仅包含有0.05kg的水。在后来的模拟中，用MIX关键字来生成与蒸发后的溶液具

71、有一样浓度的溶液，但所具有的水的质量大约为1kg。19 / 102第一种模拟数据输入组表17包含四个关键字：1TITLE用来描述包含在输出文件中的模拟，2SOLUTION是用来定义来自中部Oklahoma地区的雨水的组分，3REACTION是用来指定从液相中移除水的量，以摩尔数表示，和4SAVE是用来保存批反响计算的结果，存储为溶液2。表表 17.17. 例例 4 4 的输入数据组的输入数据组TITLE Example 4a.-Rain water evaporationSOLUTION 1 Precipitation from Central Oklahoma units mg/L pH 4

72、.5 # estimated temp 25.0 Ca .384 Mg .043 Na .141 K .036 Cl .236 C .1 CO2(g) -3.5 S(6) 1.3 N(-3) .208 N(5) .237REACTION 1 H2O -1.0 52.73 molesSAVE solution 2END除非使用了标识符-water，否那么，通过SOLUTION输入定义的所有溶液都准确地具有1kg大约为55.5mol的水。为了将溶液浓缩20倍，大约需要移除52.8mol的水55.5x 0.95。应用MIX关键字关键字的第二次模拟是将溶液中包含的所有元素的摩尔数乘以20，包括氢和氧。

73、这个过程有效地增加了液相的总质量或是体积，但保持了同样的浓度。为了便于识别，来自于MIX模拟的溶液用SAVE关键字保存为溶液3。溶液3与溶液2前面模拟中产生的具有一样的浓度，但溶液3中水的质量大约为1kg。模拟的选择输出结果见表18。在中部的Oklahoma，根据蒸散过程中水的均衡，20倍的浓度因子是合理的Parkhurst和其他, 1996。PHREEQC模拟假定蒸发和蒸散具有一样的效果，且蒸散不影响离子比率。这个假定并没有被证明，也可能是不正确的。蒸发作用之后， 所模拟溶液的组分相对于方解石、 白云石和石膏仍不饱和。 正如所想的一样，当水从反响中移除时，1 kg雨水的质量溶液1大约降低为溶

74、液2中的0.05 kg。一般而言，反响之后，所保存水的量是大约的，因为可以通过类似的水解反响、外表配位反响以与纯相的溶解和沉淀反响所消耗或生成水。减少的水对氯化物数目没有影响，因为水量的降少，氯化物的浓度20 / 102mol的摩尔mol/kgw会有所增大。第二种混合模拟增加水的质量和氯化物的摩尔数20倍。因此，氯化物的摩尔数升高了，然而氯化物的浓度在混合模拟之前溶液2和之后溶液3是一样的，这是因为水的质量按比例增大了。这些模拟结果解释了有关类似氧化复原反响的一个重点表18。批反响计算和运移计算总是保持每种氧化复原元素的液相平衡。雨水分析包括铵和硝酸盐的数据，但并不含有溶解氮。雨水的pe对初始

75、溶液中离子的分布是没有影响的，这是因为氧化复原元素各个氧化复原形态是指定的C, N, 和S。尽管在热力学平衡中硝酸盐与铵不能共存，物种形成计算允许氧化复原作用不均衡，并承受在输入数据中定义的氮的两种氧化复原态浓度，而不管热动力学平衡。在批反响蒸发步骤中，达到液相的氧化复原平衡， 会引起铵的氧化和硝酸盐的复原， 而生成溶解氮在PHREEQC中的符号为N2(aq),或N(0)。最初的批反响溶液溶液2包含有氮的平衡分布，由硝酸盐和溶解氮组成，但没有铵表18。在批反响模拟计算中发生的铵的氧化和硝酸盐的复原，是为将定义雨水组分时固有的氧化复原不均衡调整为氧化复原平衡。在模拟中会产生氮的氧化复原，即使RE

76、ACTION关键字已指定溶液中没有水分移除。防止氮氧化复原态完全平衡的唯一 方 法 是 定 义 各 个 氧 化 复 原 态 作 为 单 独 的 SOLUTION_MASTER_SPECIES 和SOLUTION_SPECIES， 例如在SOLUTION_MASTER_SPECIES过定义新成分 “Amm， 并根据AmmAmmH3, AmmH4+和其它定义NH3和其它的N(-3)AmmH3, AmmH4+等。在这种情况下，将会在含N物质和含Amm物质之间取得平衡，但在N和Amm间不存在均衡。表表 1818 例例 4 4 的选择结果的选择结果kg，千克，mol，微摩尔组分组分水的质量，kgCl，m

77、olCl，mol/kg 水硝酸盐N5，mol/kg 水溶解氮N0，mol/kg 水铵N-3，mol/kg 水方解石饱和指数白云石饱和指数石膏饱和指数溶液溶液 1 1雨水雨水1.0006.6576.65716.9014.8-9.21-19.02-5.35溶液溶液 2 2浓缩浓缩 2020 倍倍0.050026.657133.1160.1475.10-9.37-19.35-2.91溶液溶液 3 3乘以因数乘以因数 2020 混合混合1.000133.1133.1160.1475.10-9.37-19.35-2.91例例 5 5不可逆反响不可逆反响在这个例子中，通过用 PHREEQC 模拟黄铁矿氧化

78、反响证明了其模拟不可逆反响的能力。21 / 102将氧气O2和 NaCl 以五种不同的量不可逆地参加到纯水中0.0, 1.0, 5.0, 10.0,50.0mmol ；在不可逆反响中，O2和 NaCl 的相比照例分别为 1.0 和 0.5。黄铁矿、方解石和针铁矿允许溶解以达到平衡，二氧化碳分压力的保持为10大气分压力 。另外，达到过饱和时，允许生成石膏沉淀。-3.5表表 19.19. 例例 5 5 输入数据组输入数据组TITLE Example 5.-Add oxygen, equilibrate with pyrite, calcite, and goethite.SOLUTION 1 PU

79、RE WATER pH 7.0 temp 25.0EQUILIBRIUM_PHASES 1 Pyrite 0.0 Goethite 0.0 Calcite 0.0 CO2(g) -3.5 Gypsum 0.0 0.0REACTION 1 O2 1.0 NaCl 0.5 0.0 0.001 0.005 0.01 0.05SELECTED_OUTPUT file ex5.sel totals Cl si Gypsum equilibrium_phases pyrite goethite calcite CO2(g) gypsumEND表表 20.20. 例例 5 5 的选择输出结果的选择输出结果摩

80、尔转移与相集合中的摩尔数有关；正数示当前相数量的增长，也即是，沉淀；负数表示当前相数量降低，也即溶解 反响物量，毫摩尔O20.01.05.0NaClpHpe黄铁矿-.27-1.33-2.67-13.33摩尔转移，毫摩针铁矿.271.332.6713.33方解石CO2(g)石膏0.0.0.00.09.00-.93-2.94-5.56.142.405.11石膏的饱和指数-6.13-2.02-1.06-.65.00.08.28-4.94-0.0000320.000011-0.49-0.490.58.17-4.292.57.98-3.9710.05.07.88-3.8250.025.07.72-3.5

81、7-26.8426.49定义纯水在 SOLUTION 关键字中输入 表 19 。 纯相组分在 EQUILIBRIUM_PHASES中定义。缺省状态下，在纯相集合中存在 10 mol 的黄铁矿，针铁矿，方解石和二氧化碳；而石膏在纯相集合中为 0.0 mol。石膏达到了过饱和时将会沉淀；由于在当前石膏摩尔数为零，所以不能溶解。REACTION 数据块定义了所要模拟的不可逆反响。在这个例子中，氧“O222 / 102气和 NaCl 以 1:0.5 的相比照例参加。 反响量分别为 0.0, 0.001, 0.005, 0.01, 和 0.05mol。反响物以化学分子式来定义，正如在这种情况下的O2或是

82、在 PHASES 输入中所定义的相的名字。因此，缺省数据块文件中名字为“O2(g)或是“Halite的相可代替“O2或“NaCl得到一样的结果。 在每一步反响中， 参加水中的氧元素的摩尔数等于分子式“O2(2.0)中的氧的化学计量系数乘以相对系数 1.0 ， 并乘以由反响步骤中 0.0, 0.001, 0.005,0.01, 或 0.05定义的反响摩尔数。加到每一步中氯化物的摩尔数是在分子式“NaCl(1.0)中的氯化物计量系数乘以相对系数 0.5并乘以反响步骤中定义的摩尔数。 在每一步平衡计算中，SELECTED_OUTPUT 用来打印氯化物的总浓度，石膏的饱和指数，黄铁矿、针铁矿、二氧化碳

83、和石膏的总数量和转移的摩尔数到文件ex5.sel 中。例 5 中的结果概括在表 20 中。当没有氧和氯化钠加到体系中时，很少量的方解石和二氧化碳溶解，痕量的黄铁矿和针铁矿发生反响； 由于与黄铁矿达到平衡， 石膏未饱和达六个数量级饱和指数为-6.13 ，故pH 为 8.28，pe 很低-4.94 。当参加氧和氯化钠，黄铁矿氧化，针铁矿相对不溶，发生沉淀。反响产生了硫磺酸，pH 值降低，pe 轻微增大，引起了方解石溶解并释放出二氧化碳。在参加10 和 50 mmol 氧之间的某点上，方解石达到了饱和，且开始沉淀。当已经参加了50 mmol 的氧和 25 mmol 的氯化钠时，共有 9.00 mmo

84、l 的石膏沉淀出来。例例 6 6反响路径模拟反响路径模拟在这个例子中，研究了相的沉淀，它是不同的钾长石微斜长石不全等溶解的结果。在这个例子中考虑了有限的相钾长石、水铝矿、高岭石、云母白云母。这种相设置的路径最初是由Helgeson和其他人提出的 1969 。 在这个例子中， 这个相的热动力学数据 表21，PHASESPHASES关键字是来自于Robie和他人的成果1978，这个成果与在 PHREEQC是用手册中实验问题5的一样Parkhurst和他人，1980。PHREEQC能够用三种方法来解决这个问题：在相图表中，反响路径和相边界的单一交集能够被计算例子6A，反响路径能够被用来逐一的计算 6

85、B，或是反响路径能够作为动力学过程进展计算6C。在第一种方法中，没有必要知道有关反响的数量，但需要大量的模拟来找到适宜的相边界交点。 在第二种方法中， 仅要一次模拟就足够了， 但必须提前知道一定数量的反响。在第三种方法中，用动力学反响速率表达式来计算反响路径， 运用步长的大小来调整运算法那么， 这将通过在必要的时候自动地降低时间间隔处理相边界的国度。 仅仅要求总时间达到钾长石的平衡点。在这个例子中证明了所有这三种方法。PHREEQC暗中包含了完全反响路经程序的逻辑性 例如Helgeson和他人, 1970,Wolery, 1979,Wolery和他人,1990。此外，直接计算相边界交集的能力提

86、供了在相图上描述反响路径，逐步增加反响的选项和自动地寻找固定相的集合允许容易和快速计算在相边界之间的反响路径上的点。 对每23 / 102种相转移而言， PHREEQC中的动力学方法和Basic编译程序能够用来保存和打印所耗的时间和液相的组分。在概念上，这个例子考虑了将钾长石放入烧杯中， 并且允许缓慢反响时，将会发生的反响。当钾长石溶解时，其它的相便开始沉淀。在这个例子中，假定了仅有水铝矿、高岭石或钾云母能形成，如果这些相达到饱和，那么这些相将会可逆地沉淀。在反响开始时，沉淀的矿物相随着反响的进展也可能再次溶解。表表 2121例例 6 6 输入数据的设置输入数据的设置TITLE Example

87、 6A.-React to phase boundaries.SOLUTION 1 PURE WATER pH 7.0 charge temp 25.0PHASES 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 + 10 H+ = 3 Al+3 + 3 H4SiO4 + K+

88、 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 kcalSELECTED_OUTPUT -file ex6A-B.sel -activities K+ H+ H4SiO4 -si Gibbsite Kaolinite K-mica K-feldspar -equilibrium Gibbsite Kaolinite K-mica K-feldsparENDTITLE Example 6A1.-Find

89、amount of K-feldspar dissolved to reach gibbsite saturation.USE solution 1EQUILIBRIUM_PHASES 1 Gibbsite 0.0 KAlSi3O8 10.0 Kaolinite 0.0 0.0 K-mica 0.0 0.0 K-feldspar 0.0 0.024 / 102ENDTITLE Example 6A2.-Find amount of K-feldspar dissolved to reach kaolinite saturation.USE solution 1EQUILIBRIUM_PHASE

90、S 1 Gibbsite 0.0 0.0 Kaolinite 0.0 KAlSi3O8 10.0 K-mica 0.0 0.0 K-feldspar 0.0 0.0ENDTITLE Example 6A3.-Find amount of K-feldspar dissolved to reach K-mica saturation.USE solution 1EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 0.0 K-mica 0.0 KAlSi3O8 10.0 K-feldspar 0.0 0.0ENDTITLE Example 6A4

91、.-Find amount of K-feldspar dissolved to reach K-feldspar saturation.USE solution 1EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 0.0 K-mica 0.0 0.0 K-feldspar 0.0 KAlSi3O8 10.0ENDTITLE Example 6A5.-Find point with kaolinite present, but no gibbsite.USE solution 1EQUILIBRIUM_PHASES 1 Gibbsite 0

92、.0 KAlSi3O8 10.0 Kaolinite 0.0 1.0ENDTITLE Example 6A6.-Find point with K-mica present, but no kaoliniteUSE solution 1EQUILIBRIUM_PHASES 1 Kaolinite 0.0 KAlSi3O8 10.0 K-mica 0.0 1.0ENDTITLE Example 6B.-Path between phase boundaries.USE solution 125 / 102EQUILIBRIUM_PHASES 1 Kaolinite 0.0 0.0 Gibbsit

93、e 0.0 0.0 K-mica 0.0 0.0 K-feldspar 0.0 0.0REACTION 1 K-feldspar 1.0 0.04 0.08 0.16 0.32 0.64 1.0 2.0 4.0 8.0 16.0 32.0 64.0 100 200 umolENDTITLE Example 6C.-kinetic calculationSOLUTION 1 units mol/kgw Al 1.e-13 K 1.e-13 Si 3.e-13EQUILIBRIUM_PHASES 1 Gibbsite 0.0 0.0 Kaolinite 0.0 0.0 K-mica 0.0 0.0

94、KINETICS 1K-feldspar parms 1.36e-11 m0 2.16 m 1.94 step_divide 1e-6 steps 1e2 1e3 1e4 1e5 1e6 1e7 1e8INCREMENTAL_REACTIONS trueRATESK-feldspar-start 10 REM store the initial amount of K-feldspar 20 IF EXISTS(1) = 0 THEN PUT(M, 1) 30 REM calculate moles of reaction 40 SR_kfld = SR(K-feldspar) 50 mole

95、s = PARM(1) * (M/M0)0.67 * (1 - SR_kfld) * TIME 60 REM The following is for printout of phase transitions 80 REM Start Gibbsite 90 if ABS(SI(Gibbsite) 1e-3 THEN GOTO 150 100 i = 2 110 GOSUB 1500 150 REM Start Gibbsite - Kaolinite 160 if ABS(SI(Kaolinite) 1e-3 THEN GOTO 20026 / 102 170 i = 3 180 GOSU

96、B 1500 200 REM End Gibbsite - Kaolinite 210 if ABS(SI(Kaolinite) 1e-3 OR EQUI(Gibbsite) 0 THEN GOTO 250 220 i = 4 230 GOSUB 1500 250 REM Start Kaolinite - K-mica 260 if ABS(SI(K-mica) 1e-3 THEN GOTO 300 270 i = 5 280 GOSUB 1500 300 REM End Kaolinite - K-mica 310 if ABS(SI(K-mica) 1e-3 OR EQUI(Kaolin

97、ite) 0 THEN GOTO 350 320 i = 6 330 GOSUB 1500 350 REM Start K-mica - K-feldspar 360 if ABS(SI(K-feldspar) 1e-3 THEN GOTO 1000 370 i = 7 380 GOSUB 1500 1000 SAVE moles 1010 END 1500 REM subroutine to store data 1510 if GET(i) = M THEN RETURN 1520 PUT(M, i) 1530 PUT(TOTAL_TIME, i, 1) 1540 PUT(LA(K+)-L

98、A(H+), i, 2) 1550 PUT(LA(H4SiO4), i, 3) 1560 RETURN-endUSER_PRINT 10 DATA A: Gibbsite, B: Gibbsite- Kaolinite ,C: Gibbsite - Kaolinite , D: Kaolinite - K-mica , E:Kaolinite- K-mica , F: K-mica - K-feldspar 20 PRINT Transition Time K-feldspar LA(K/H)LA(H4SiO4) 30 PRINT reacted 40 PRINT (moles) 50 FOR

99、 i = 2 TO 7 60 READ s$ 70 IF EXISTS(i) THEN PRINT s$, GET(i,1), GET(1) - GET(i), GET(i,2), GET(i,3) 80 NEXT iSELECTED_OUTPUT file ex6C.sel reset false27 / 102USER_PUNCH heading pH+logK logH4SiO4 10 PUNCH LA(K+)-LA(H+) LA(H4SiO4)END输入数据的设置表21，首先在SOLUTIONSOLUTION输入中定义纯水，在PHASESPHASES输入中定义相的热动力学数据。一些矿物

100、在数据库文件中已经定义phreeqc.dat， 但是在运行过程中，在输入数据设置中包含的相代替以前的定义数据库文件没有改变。SELECTED_OUTPUTSELECTED_OUTPUT是用来产生在表22中列出的所有数据的文件，这些数据被用来会制图6。指定了钾离子、氢离子和硅酸活度系数的对数；水铝矿、 高岭石、钾云母、钾长石的饱和指数以与在相集合中的总数量，并且在每一次计算之后，水铝矿、高岭石、钾云母、和钾长石的摩尔转换都将写到文件ex6A-B.sel中。SELECTED_OUTPUTSELECTED_OUTPUT的定义在运行所有模拟中都保持有效，直到读入新的SELECTED_OUTPUTSEL

101、ECTED_OUTPUT数据块，或是直到在PRINTPRINT数据块中，用-selected_output-selected_output标识符将写到文件暂停。表表 22.22. 例例 6A6A 的选择输出的选择输出例6A输入数据组中模拟的参考标志。负摩尔转移标明溶解，正摩尔转移标明沉淀。图形上的点标明的是例6中的标志点钾长石钾长石模模拟拟的摩尔的摩尔转移转移, ,微微摩尔摩尔6A16A26A36A46A56A6+ +活度对数活度对数摩尔转移，微摩尔摩尔转移，微摩尔水水铝铝矿矿高高岭岭石石水水铝铝矿矿饱和指数饱和指数高高岭岭石石钾钾云云母母钾钾长长石石图图形形上上的的点点ABDFCEH HKH

102、 H4 4SiOSiO4 4H钾钾 云云母母-0.03-2.18-7.01-0.57-7.100.000.000.000.0-3.8-10.7-14.7-8.212.55-5.201.78.00.00.0-.7.0.0-1.9-5.9.0.0-2.5.0-20.02-9.114.41-4.47.009.71.00-190.9-9.395.49-3.55.00.0063.61-2.0-.7-3.02-8.352.83-5.20.001.24.00.0.0.0-1.6-5.6.0-2.1-32.68-9.074.41-4.25.00.0010.78-.9模 拟 6A1 允 许 钾 长 石 参 加

103、反 响 ， 直 到 达 到 水 铝 矿 矿 物 的 溶 解 平 衡 。 这 在EQUILIBRIUM_PHASESEQUILIBRIUM_PHASES输入中设置的，它是通过指定水铝矿平衡饱和指数为0.0和一个选择的反响达到其平衡，KAlSi3O8钾长石的分子式。大量的钾长石10.0mol可确保能达到水铝矿的平衡。高岭石，钾云母，和钾长石在达到饱和之后，可允许生成沉淀，但不允许其溶解，这是因为在相集合中，所给定它们的初始摩尔数为0。在这个模拟中所计算的反响的数量是刚好溶解足够的钾长石以达到与水铝矿的平衡， 这也可能包含一种或是几种其它的矿物沉淀。没有水铝矿将会溶解或是沉淀； 选择的反响物KAlS

104、i3O8将会在它们本来的位置溶解或沉淀。模拟6A2-6A4中对高岭石，钾云母、和钾长石进展了一样的计算。在其它28 / 102的温度下或是应用其它的矿物时， 可能的情况是，不管所增加的替代反响的数量有多少， 目标相仍不饱和，这是因为相对于其它相而言，该相不稳定。在这种情况下， 这种数字化的方法将会找到产生最大饱和指数的选择反响的数量。图图6 625下，纯水中钾长石溶解的相图表微斜长石说明了稳定相边界的交叉点例 6A和穿过稳定区域的反响路径例6B。图表的建立用到了来自Robie和他人1978所做的水铝矿，高岭石，钾云母白云母，和微斜长石的热动力学数据。模拟6A1-6A4的选取结果见表22， 并且

105、在图6中以点A、 B、 D和F标出。 基于热动力学数据，在图上标出了相的稳定区域，但没有建立这些模拟来进展计算。在点B，水铝矿开始形成高岭石，这个反响消耗Si。这个反响路径必须沿着水铝矿-高岭石相的边界达到一些中间点C，直到水铝矿被转换，然后这个路径穿过高岭石区域到D点。相似的，在高岭石-K-云母相边界会形成点E，在那里反响路径开始穿过K-云母区域到点F。模拟6A5和6A6表21解释了这两点。在模拟6A5中，通过允许钾长石溶解，并在某一点上，高岭石达到饱和且出现在相集合中，而水铝矿饱和却不出现在相集合中计算得到点C。同样，模拟6A6解答了钾云母处于饱和态且出现在相集合中，而高岭石达到饱和，但不

106、出现在相集合中的点。 在6A5和6A6中，为高岭石和钾云母指定1 mol的初始量是任意的，这个数量仅仅要求必须足够达到矿物平衡。决定反响路径的一种简单方法是增加钾长石的反响，允许在沿路径的每一点形成水铝矿，高岭石，钾云母，和钾长石之间稳定的相集合。在这种方法中，唯一的难点是知道所加反响的适当的数量。在表22中，从点A到点F钾长石的溶解围从0.03到190.0 mmol。在6B局部中表21，反响增量的对数围是用来定义路径固体线，这个路径穿过相图，从水铝矿的平衡A点开始到钾长石的平衡F点完毕。然而，A点到F点的准确位置不是由B局部中应用的任意的反响增量决定的。 通过6B局部中计算的反响路径标出在图

107、6中的相图中， 从A点到在6A局部F点。29 / 102最后，在动力学方法中， 钾长石的动力学分解是遵循不同的时间变量而决定的， 随着动力学反响的进展，允许水铝矿、高岭石和钾云母沉淀或再次溶解。 SOLUTIONSOLUTION 1 1定义为溶解有少量的钾长石1e-13 moles，那么这种溶液包含有与EQUILIBRIUM_PHASESEQUILIBRIUM_PHASES中的相相关的所有元素，虽然程序的成功运行不要求，但可以排除一些警告信息。在反响速率的积分期间，基于过渡态理论，假定了一个简单的溶解速率定律。RK feldspar k1这里，k = 1e-16 mol/cm /s。2AmIA

108、P()(1()K feldspar)V m0K这个KINETICS数据块是用来输入动力学模拟的指定数据。 动力学反响的化学计量是钾长石的化学分子式；缺省状态下，速率的名字是假定在PHASES数据块中所定义的相， 且相分子式用作反响的化学计量系数。 据假定，在0.1 立方毫米的古土壤中含有10%的钾长石，其中，b/ 6g /cm3，因此，A/V 136/cm。用-parms-parms标识符在 KINETICSKINETICS数据块中输入k1A1.36e 11V的值假定1kgw=1升，且能够在RATESRATES数据块中的Basic速率定义中作为“PARM(1)调用。假定土壤在一定程度上风化，最

109、初的钾长石仅余下了 90%-m0-m0 2.16 和-m-m1.94，这里m0m0表示初始质量1 kg soil x 0.1= 100 g / 278.3 g/mol = 0.359 mol/kg x 6kg/L = 2.16 mol/L，m为剩余的质量2.16的90%为1.94mol/L。在任何时间间隔反响的最大数量限制在1e-6mole(-step_divide-step_divide 1e-61e-6)。时间步长秒以 -steps-steps标识符来定义。INCREMENTAL_REACTIONS trueINCREMENTAL_REACTIONS true使得模拟的总时间为所有时间步长

110、的总和1.111111e8秒；每一次新的时间步长都在前面时间步长的完毕后开始。钾长石分解速率在RATESRATES数据块中以Basic声明的形式定义。为了证明Basic编译的一些特点，Basic程序同样在相转换时确定和保存一些信息，这在运行完毕时通过USER_PRINTUSER_PRINT来打印。确定相转换的准确性是通过用户可定义的积分的准确性决定的。小的公差-tol-tol，大于1的大的-step_divide-step_divide初始时间间隔将会通过这个数字来剖分，以与小于1的-step_divide-step_divide 反响所指定的最大数量 将会产生小一些的时间间隔膜更准确的相转换

111、证明。在这个例子中，-step_divide-step_divide设置为1e-6，这将会限制任何的时间段反响的最大数量小于1micromole。这样，达到相转换的反响的数量应当为的精度1micromole。然而，在积分反响期间，限制反响数量需要更小的时间间隔，因此，需要更多的时间间隔用来完成积分反响，这将会增加CPU的运行时间。Basic程序不同局部函数的描述见表23。函数PUT, GET,和EXISTS用作操纵静态，和动态存贮的数据。应用在PUT声明中的下标指明唯一的数据。EXISTS能够用来决定是否具有给定下标的数据已被存贮，GET是用来检索已被存贮的数据。 一旦一种数据已经用PUT保存

112、，那么将以后的运行中存在， 除非它被另外具有一样下标的PUT声明所覆盖。 用PUT储存的数据可被Basic程序检索，包括那些在RATES, USER_PRINT,RATES, USER_PRINT, 和USER_PUNCHUSER_PUNCH数据块中定义的数据。在这个例子中，数据由RATESRATESBasic程序储存，USER_PRINTUSER_PRINTBasic程序检这些数据，并打印相转换的总和。在一个时间步长的动力学积分期间 为达到所需求的准确性， 可能需要对许多的30 / 102时间间隔进展积分 ， RATESRATES程序可运行屡次， 在每个时间步长的积分完毕时运行USER_PR

113、INTUSER_PRINT程序。表表 23.23. 钾长石溶解动力学和相转移标识的钾长石溶解动力学和相转移标识的 BasicBasic 程序描述程序描述行号行号2040-5090-110160-180200-230250-280300-330350-38010001010作用作用保存钾长石的初始摩尔数1.94mol在由TIMETIME给定的时间间隔，钾长石溶解速率的积分在水铝矿达饱和时，确定当前钾长石的最大量参加反响的最小数量在高岭石达饱和时，确定当前钾长石的最大量在高岭石达饱和，但水铝矿没有达到饱和时，确定当前钾长石的最大量钾云母达饱和时，确定当前钾长石的最大量在钾云母达饱和，但高岭石没有达

114、到饱和时，确定当前钾长石的最大量在钾长石达饱和时，确定当前钾长石的最大量保存积分反响Basic程序的完毕1500-1560保存相转移值的子程序。对序号i而言，如果钾长石的数量比当前保存的值大，保存钾长石的数量，累积时间，由氢离子所划分的钾长石活度比率的对数，和硅酸的活度对数。表24列出了在例6C中最后一个时间步长完毕时所遇到的相转换。对每一种相的转移而言，给出了发生相转移的时间、参加反响的钾长石的总数量以与图6中向转移的坐标。尽管表24中所给定的值是近似的， 钾长石的数量和相转换的坐标可与表22中的值相比拟。 正如所期望的一样，达到相转换的近似转移的摩尔数是与表22中列出值相差在1微摩尔之。表

115、表 24.24. 在例在例 6C6C 中模拟了钾长石的动力学溶解，相转移由中模拟了钾长石的动力学溶解，相转移由 RATESRATES BasicBasic 程序确程序确定，打印到输出文件中由定，打印到输出文件中由 USER_PRINT BasicUSER_PRINT Basic 程序控制。程序控制。-User print- Transition Time K-feldspar LA(K/H) LA(H4SiO4) reacted (moles)A: Gibbsite 1100 1.4048e-007 3.5755e-001 -6.3763e+000B: Gibbsite - Kaolinite

116、 1.7434e+005 2.2064e-006 2.5609e+000 -5.1950e+000C: Gibbsite - Kaolinite 2.3929e+005 3.0284e-006 2.8352e+000 -5.1943e+000D: Kaolinite - K-mica 1.5869e+006 2.0070e-005 4.4080e+000 -4.4659e+000E: Kaolinite - K-mica 2.5972e+006 3.2791e-005 4.4103e+000 -4.2509e+000F: K-mica - K-feldspar 4.7840e+007 1.90

117、72e-004 5.4879e+000 -3.5540e+000SELECTED_OUTPUTSELECTED_OUTPUT数据块指定了应用于这个模拟中的新的选择输出文件，ex6C.sel，并且所有打印到选取的输出文件都将停止-reset-reset为假。USER_PUNCHUSER_PUNCH数据块使得写到选取输出文件中的每一行中有两栏， 即钾离子与氢离子活数比率的对数和硅酸活度的对数。 在每31 / 102一个时间步长已被模拟之后-steps,KINETICS-steps,KINETICS数据块，这个数据将会被写入。下表列出了写到ex6C.sel中的结果。表表 25.25. 模拟了钾长石的

118、动力学溶解的例模拟了钾长石的动力学溶解的例 6C6C 中由中由 USER_PUNCH BasicUSER_PUNCH Basic 程序写到程序写到选取输出文件中的结果选取输出文件中的结果“t标明的是制表符 pH+logKt-6.0002e+00t-1.8975e+00t-8.5316e-01t 3.5755e-01t 2.1823e+00t 4.1360e+00t 5.2614e+00tlogH4SiO4t-1.2524e+01t-8.4212e+00t-7.3798e+00t-6.3763e+00t-5.3818e+00t-4.6005e+00t-3.7187e+00t5.4884e+00t

119、-3.5536e+00t例例 7 7气相计算气相计算该例中运用 PHREEQC 模拟了定压或定体积条件下， 气相与液相均衡时其成分的变化。 对于固定压力气相，当气相成分偏压力的总和大于气相的指定压力时， 将形成气泡。一旦气泡形成，其成分和体积将随着反响程度的变化而变化。 对于固定体积气相， 水溶液与定体积的顶部空间相联系。气相总是存在于顶部，并且压力和组分随反响程度的变化而变化。气液反响可以用 PHREEQCPHREEQC 以三种方法模拟： 1可以用EQUILIBRIUM_PHASESEQUILIBRIUM_PHASES 数据块定义一种气体在固定的偏压力下进展反响， 2定压力多成分气相可以用具

120、有-fixed-pressure-fixed-pressure 标识符默认的 GAS_PHASEGAS_PHASE 数据块模拟，或3定体积多成分气相可用具有-fixed-volume-fixed-volume 标识符的 GAS_PHASEGAS_PHASE 数据块模拟。从概念上讲，定偏压力是假定气相的容器为无限大，就象在大气或非饱和带中一样。 无论反响程度如何，容器中气体成分的偏压力不变。如果气体成分的容器是无限的， 那么气相的压力恒定， 就象河口和湖底的沉积环境中的气泡一样，那么适宜选择定压力气相。 如果气体成分的容器是有限的， 且气相体积恒定，就象实验中固定的顶部空间一样，那么应选择定体积

121、气相。在这个例子中， 用 GAS_PHASEGAS_PHASE 数据块模拟在固定压力或固定体积条件下纯水中有机物的分解，按照分解反响，假定碳、氮、氢和氧按化学方程式 CH2O(NH3).07释放。在纯水中没有其它的电子承受器， 相关的微生物分解反响是甲烷的生成。 通过微生物分解释放的碳和氮假定为发生了氧化复原反响，并达到气液均衡。溶液中碳C(+4) 和 C(-4) (甲烷)由默认的SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES 或 SOLUTION_SPECIESSOLUTION_SPECIES 数据块定义；没有定义中间价态的碳。溶液中氮元素的化合

122、价可能为 +5、+3、0 和-3。考虑的气体成分为二氧化碳(CO2)、甲烷(CH4)、氮气(N2)和氨气(NH3)。32 / 102图7. 对气相而言的，在定容积和定压力下的纯水中含有组分CH2O(NH3)0.07的有机质分解期间的气相组分。氨气的分压力不小于大气中的10 没有标出 。-7在第一次模拟中，初始水定义为二氧化碳(CO2)的偏压力为 10-1.5时与方解石达到均衡。纯水通过 SOLUTIONSOLUTION 数据块中的默认值(pH = 7, pe = 4, temperature = 25)定义；方解石和 CO2由 EQUILIBRIUM_PHASESEQUILIBRIUM_PHA

123、SES 数据 块定义；用SAVE 保存均衡溶 液(表 26)。用SELECTED_OUTPUTSELECTED_OUTPUT 数据块定义文件(ex7.sel)，用于写入每次计算中指定的数据。所有的打印到文件通过标识符-reset-reset 设置为 falsefalse 来停止。在每个模拟的每次计算中，附加的标识符使得指定的数据项被写到选择的输出文件： -simulation-simulation， 模拟数， -state-state， 计算类型 初始溶液、反响、运移和其它 ；-reaction-reaction，每次计算中增加的反响量如REACTIONREACTION 数据块所定义 ；-si

124、-si，指定每种矿物的饱和指数或每种气体偏压力的对数； 和-gas-gas, 指定的每种气体成分的摩尔数。在第二次模拟中，将碳氮比大约为 15:1 的有机物分解反响不可逆地从 1 到 1000mmol增加到溶液中 (REACTION 关建字)。当偏压力大于1.1 atm 时，允许生成气相中最初不存在的气体成分；在气相中只允许有CO2、CH4、N2和 NH3(GAS_PHASE 数据块)气体。在第三次模拟中，使用和第二次模拟同样的初始溶液和反响。 气相开始无气体成分存在， 但是定义固定体积为 22.5 L，即 1.0 mol 有机物反响之后定压力气相的体积。1.0 mol 有机物反响后，定压力和

125、定体积气相将有同样的压力、体积和成分；在所有其它的反响中， 增加压力、体积和成分时，两个气相将有所不同。表表26.26. 例例7 7输入的数据组输入的数据组TITLE Example 7.-Organic decomposition with fixed-pressure and fixed-volume gas phasesSOLUTION 1EQUILIBRIUM_PHASES 133 / 102 Calcite CO2(g) -1.5SAVE solution 1SELECTED_OUTPUT -reset false -file ex7.sel -simulation true -st

126、ate true -reaction true -si CO2(g) CH4(g) N2(g) NH3(g) -gases CO2(g) CH4(g) N2(g) NH3(g)END#Simulation 2: Decomposition of organic matter, CH2O(NH3).07,#at fixed pressure of 1.1 atmUSE solution 1GAS_PHASE 1 Fixed-pressure gas phase fixed_pressure pressure 1.1 CO2(g) 0.0 CH4(g) 0.0 N2(g) 0.0 NH3(g) 0

127、.0REACTION 1 CH2O(NH3)0.07 1.0 1. 2. 3. 4. 8. 16. 32 64. 125. 250. 500. 1000. mmolEND#Simulation 3: Decomposition of organic matter, CH2O(NH3).07,#at fixed volume of 22.5 LUSE solution 1USE reaction 1GAS_PHASE 1 Fixed volume gas phase fixed_volume volume 22.5 CO2(g) 0.0 CH4(g) 0.0 N2(g) 0.0 NH3(g) 0

128、.0END对于定压力气相，约增加3 mmol 反响时，形成气泡见图7 。初始气体成分CH4占 90%以上， CO2不足 10%， 以与很少量的 N2和 NH3(整个批反响计算中 NH3的偏压力少于 10 atm ，图上没有表示)。产生气体体积的围从3 mmol 反响时的 1 mL到 1 mol 反响时的 22.5 L。反响增加到 1 mol 后，几乎所有的碳和氮都以气相存在。对于定体积气相，在所有的有机物分解过程中都存在气相见图 7 。初始气体成分主34 / 102-7要是 CO2，以与一定量的 CH4和少量 N2。随着气相成分的变化，CO2和 CH4的偏压力变得几乎相等。N2的偏压力总是比

129、CO2和 CH4小一个数量级；NH3的偏压力总是很小图中没有表示 。反响达到 1 mol 之前，定体积气相的偏压力小于定压力气相的偏压力。如果继续反响超过1.0 mol，定体积气相的偏压力将继续增加，并将大于定压力气相的偏压力，并保持不变。相反， l 反响达到 1.0 mo 之前，定压力气相的体积小于定体积气相的体积。 如果反响持续，超过 1.0 mol 时，定压力气相的体积将大于定体积气相的体积。？例？例 8 8外表络合作用外表络合作用PHREEQC 包括三个外表络合模型： 1根据默认值，使用没有明确计算扩散层组分的双层模型。 (2) 也可选择有明确计算扩散层组分的静电双层模型(-diffu

130、se_laye-diffuse_layer)。 (3) 最后，可能选择使用非静电模型 (-no_edl-no_edl)。根据下面的修改，在 Dzombak and Morel (1990)将静电模型归纳为双层模型： (1) 外表层可能有超过两种的键联位置， 2不包括外表析出， 3根据 Borkovec 和 Westall (1983)，可选择使用电荷势联系公式，它明确计算了扩散层组分 -diffuse_layer-diffuse_layer 。 非静电模型不考虑外表络合形成中产生的外表电荷的影响，其外表络合的数学处理结果很象没有活度系数项的水溶液络合。 以下归纳除来的没有明确计算扩散层组分的例

131、子来自于Dzombak 和 Morel (1990, 第 8 章)。在氧化物外表， 假定有强和弱两种类型的位置进展水合氧化铁吸附锌的模拟。 质子和锌离子两种键型竞争， 可用质量作用描述。 外表物质的活度大小取决于由外表电荷的变化引起的外表电势。该例中，描述了在 0.1m 硝酸钠电解液中锌浓度较低(10-7 m)和较高(10-4 m)时，以 PH值为函数的水合氧化铁吸附锌的变化。要求三个关建字数据块确定模拟的外表络合数据：SURFACE_MASTER_SPECIESSURFACE_MASTER_SPECIES ， SURFACE_SPECIESSURFACE_SPECIES ， 和 SURFAC

132、ESURFACE 。 在 默 认 数 据 块 文 件 中SURFACE_MASTER_SPECIESSURFACE_MASTER_SPECIES 数据块用两个键位定义了以名为“Hfo 水合氧化铁 的化学键。在数据库文件中，化学键的名字由外表物质的名字“Hfo组成，之后可加上一个下划线合一个小写的化学键名表示， “Hfo_w和“Hfo_s分别表示强和弱。 只有单个外表物质存在两个或多个化学键时， 必须用。 下划线符号中， 一个化学键位用一个摩尔平衡方程表示， (该例中为 Hfo_w 和 Hfo_s)， 外表用单个电荷电势方程和电荷平衡方程表示(该例中为 Hfo)。这样， 每个变面层物质中， 每个

133、化学键上的电荷将需要单个电荷电势方程或电荷平衡方程表示。表表 2727例例 8 8 输入数据的设置输入数据的设置TITLE Example 8.-Sorption of zinc on hydrous iron oxides.SURFACE_SPECIES35 / 102 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

134、_wO- + H+ log_k -8.82 Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+ log_k -2.32SURFACE 1 Hfo_sOH 5e-6 600. 0.09 Hfo_wOH 2e-4SOLUTION 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+ l

135、og_k 0.0END#Zn = 1e-7#SELECTED_OUTPUT file ex8.sel molalities Zn+2 Hfo_wOZn+ Hfo_sOZn+USE solution 1USE surface 136 / 102EQUILIBRIUM_PHASES 1 Fix_H+ -5.0 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.25 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.5 N

136、aOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.75 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.0 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.25 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.5 NaOH 10.0EN

137、DUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.75 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.0 NaOH 10.0ENDUSE solution 137 / 102USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.25 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.5 NaOH 10.0ENDU

138、SE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.75 NaOH 10.0ENDUSE solution 1USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -8.0 NaOH 10.0END#Zn = 1e-4#USE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.0 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.25 NaOH 10.0END

139、USE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.5 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -5.75 NaOH 10.0ENDUSE solution 2USE surface 138 / 102EQUILIBRIUM_PHASES 1 Fix_H+ -6.0 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.25 NaOH 10.0ENDUS

140、E solution 2USE surface 1EQUILIBRIUM_PHASES 2 Fix_H+ -6.5 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -6.75 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.0 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.25 NaOH 10.0ENDUSE solution

141、 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.5 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -7.75 NaOH 10.0ENDUSE solution 2USE surface 1EQUILIBRIUM_PHASES 1 Fix_H+ -8.0 NaOH 10.0ENDPHREEQC 中，Dzombak 和 Morel (1990)总结的外表络合反响在默认数据库文件中由39 / 102SURFACE_SPECIES 数据块确定。但是，在 Dzombak a

142、nd Morel (1990, 第 8 章) 例中，使用的固有稳定常数与它们的和不同，因此，在输入文件中由SURFACE_SPECIES 数据块表27专门确定。Dzombak 和 Morel (1990, p. 259)中的质量作用方程在输入数据组中表27给出。其中活度系数或不包括在质量作用方程中；程序本身包括电势项在。外表层集合的组分和其它他特征由SURFACE 数据块确定。 具有多个结合位置的多外表层组分可能也在这个数据块确定。 对于每个外表，位置每种类型的摩尔数、 外表层的初始成分和外表积必须确定。 外表层的组分随着反响程度的变化而变化。 结合位置和外表积的数量可能保持不变，如果外表层与

143、一个均衡相或热动力反响有关，那么可能变化。在这个例子中，确定了一个外表层(Hfo) 和两个结合位置(Hfo_w and Hfo_s) ，并且结合位置和了外表积的数量是固定不变的。强键Hfo_s 5*10-6 ，弱键Hfo_w 的摩尔数是 2*10-4。最初，所有的外表位置是不带电的，以质子形式存在。整个外表 Hfo 的外表积必须用两个数确定，即外表物质的质量比外表 thearea per mass ( 该例中为 600 m2/g)和外表物质的总质量 (该例中为0.09 g)。习惯上用这两个数来确定外表积，但是，在这个模型中用这两个数的积来定义外表积；不分别使用单个数。一个外表的任何键位的外表积

144、可能由这个数来输入；在该例中，外表积用 Hfo_s 来输入。Figure 8总锌浓度为10-7和10-4克分子浓度水合氧化铁的强合弱外表位置作为PH值的函数时，氧化锌在液相中的分布。40 / 102为了确定模拟的初始条件， 用不同浓度的氧化锌来确定两个硝酸钠溶液(SOLUTION 1 和2 数据块)。用 PHASES 数据块来确定假相 “Fix_H+。这个相不是真实的，但在每个批反响计算中都使用，以调整PH 值到固定值。最后，“USE surface none 行在最初的模拟中自行消除了定义的用批反响计算。默认情况下。如果在模拟中定义了一个 SOLUTION 和SURFACE 数据块，那么在模

145、拟中定义的第一个溶液 (在该例中为 SOLUTION 1)和在模拟中定义的第一个外表结合起来可能和其它集合以与一个气相允许平衡。 对于批反响计算，具有“solution none的 USE 关建字从系统中排除溶液。而且不执行批反响计算。 (无论何时在模拟中定义了一种溶液或混合物， 即暗含确定了一个批反响， 并且在同一个模拟中确定了任意一个关建字数据块 EXCHANGE, QUILIBRIUM_PHASES,GAS_PHASE, KINETICS, REACTION,REACTION_TEMPERATURE, SOLID_SOLUTIONS, 或 SURFACE)。对于 PH 值为 5 到 8

146、时，输入数据集中与剩余的溶液具有溶液1 或 2 的外表集合平衡。在单一的模拟中，可能用 REACTION 数据块增加不同量的 NaOH 到一种溶液中，但是反响增量不能产生均匀间隔的 PH 值，并且增量的大小预先不知。在这个例子中，采取不同的方法产生在要求的 NaOH 量时均匀间隔的 PH 值，但是需要许多模拟来达到预期的PH 值。在具有变化饱和指数的 EQUILIBRIUM_PHASES 数据块中，每额模拟用相“Fix_H+ 调节 PH 值。在每个模拟中增加或排除反响NaOH 以达到指定的饱和指数， 根据反响的定义， “Fix_H+ 数值上等于氢活度的对数或负的PH 值。注意，尽管在所有的这些

147、模拟中，可以达到预期的PH值，但太小的 PH 值将使得程序出错，因为即使除去溶液中所有的钠，也不能达到如此低的PH 值。模拟的结果表示在图 8 中，并与 Dzombak 和 Morel (1990, 图 8.9)中结果一致。PH 值较高时比拟低时对锌的吸附更强烈。另外，锌浓度较低时，在整个PH 值围，强健对锌的吸附均可大于弱键， 且 PH 值较高时， 大局部留在强键处。 锌的浓度较大， 且仅在 PH 值较低时，强健站主导地位。因为在PH 值较高时，所有的强健都被填充，PH 值较高，且锌的浓度较大时，大局部锌被吸更多的弱键位。例例 9 9 溶解铁离子与氧的动力学氧化复原反响溶解铁离子与氧的动力学

148、氧化复原反响动力学速率表达式的定义完全可以以一种一般的方式来定在，在PHREEQC中是在RATESRATES数据块中应用Basic声明。速率表达式可以在KINETICSKINETICS数据块中应用在多组反响或是运移反响中。对运移反响而言(ADVECTIONADVECTION或TRANSPORTTRANSPORT)，在关键字KINETICSKINETICS(KINETICSKINETICSm-n)之后，动力学反响能够以数字围的形式区域来定义。速率表达式是含第四次序和第五次序的Runge-Kutta-Fehlberg算法的综合。 在动力学计算被初始化之前和当动力学反响增量被增加时，平衡才被计算。

149、平衡的计算包括溶液中离子的种类、 各种交换物质、 平衡相、 固体溶液、溶液物质外表的集合， 和以与已经定义的气相。 应执行检查以确保在第四次序和第五次序经过一段时间间隔的速率评估差异不变化这是为不超过用户所定义的误差。 如果这个误差并不41 / 102满以，那么对时间的积分会自动地以更小的时间间隔来开始。在固体相和水相之间的动力学反响能够不修改任何的数据库文件进展计算。PHREEQC也能够计算正常情况下认为是在平衡相下的水中离子的动力学反响， 但这需要数据库来重新定义。动力学反响的水中的离子用SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES来作为独立

150、的元素进展定义是必要的。 这个例子证明了一种元素 铁 的两种化合价的状态， 以与说明如何在水中用PHREEQC来计算Fe2+到Fe3+的动力学氧化复原反响。由(Singer and Stumm, 1970)给定的在水中的Fe2+被氧化的氧化复原反响的速率是：dmFe2dt (2.91e 91.33e12a2OHPO2)mFe2这里t是时间，以秒表示，aOH-是氢氧离子的活度，是溶液中的mFe2+总摩尔数，是氧气的分压力atm。铁离子完全氧化复原反响时间是当pH值高于7.0时在溶液入二氧化碳的时间问题。 然而，Fe3+形式溶解OH-的配合物，它同样也作为铁离子的氢氧化合物来沉淀，因此，在氧化复原

151、反响期间，pH会降低。因为氧化复原速率对OH-的活度具有二次独立性，那么，这个速率随着pH的降低而迅速的变小。 这个速率等式在非缓冲溶液中是高度非线性的， 且一定是需要以数字来综合。这个例子模拟了反响容器中具有10mmolNaCl/kgw和0.1mmolFeCl2/kgw的溶液，其中pH=7.0，通过其中的空气是缓缓流入的；在溶液组分随时间的变化是计算得来的。这个计算需要分开Fe(2) 和Fe(3)离子之间平衡。有两种新的“元素定义在SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES中“Fe_di，这与Fe(2)相关；和“Fe_tri，这与Fe(3)相

152、关。这些元素的主要离子是被定义为Fe_di+2和Fe_tri+3，以与所有溶液中的离子、相、交换的离子和溶液外表的离子一定要使用这些新的元素的主要的离子来重新打印。 一些转换是列在表28中，这给出了这个例子的局部输入文件。表表2828例例9 9局部输入数据的设置局部输入数据的设置TITLE Example 9.-Kinetically controlled oxidation of ferrous iron. Decoupled valence states of iron.SOLUTION_MASTER_SPECIESFe_di Fe_di+2 0.0 Fe_di 55.847Fe_tri

153、Fe_tri+3 0.0 Fe_tri 55.847SOLUTION_SPECIESFe_di+2 = Fe_di+2 log_k 0.0Fe_tri+3 = Fe_tri+3 log_k 0.0# Fe+2 species#Fe_di+2 + H2O = Fe_diOH+ + H+42 / 102 log_k -9.5 delta_h 13.20 kcal#. and also other Fe+2 species#Fe_di+2 + Cl- = Fe_diCl+ log_k 0.14Fe_di+2 + CO3-2 = Fe_diCO3 log_k 4.38Fe_di+2 + HCO3- =

154、 Fe_diHCO3+ log_k 2.0Fe_di+2 + SO4-2 = Fe_diSO4 log_k 2.25 delta_h 3.230 kcalFe_di+2 + HSO4- = Fe_diHSO4+ log_k 1.08Fe_di+2 + 2HS- = Fe_di(HS)2 log_k 8.95Fe_di+2 + 3HS- = Fe_di(HS)3- log_k 10.987Fe_di+2 + HPO4-2 = Fe_diHPO4 log_k 3.6Fe_di+2 + H2PO4- = Fe_diH2PO4+ log_k 2.7Fe_di+2 + F- = Fe_diF+ log_

155、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 kcalFe_tri+3 + 3 H2O = Fe_tri(OH)3 + 3 H+ log_k -12.56 delta_h 24.8 kcalFe_tri+3 + 4 H2O = Fe_tri(OH)4- + 4 H+ log_k -21.6

156、delta_h 31.9 kcal43 / 1022 Fe_tri+3 + 2 H2O = Fe_tri2(OH)2+4 + 2 H+ log_k -2.95 delta_h 13.5 kcal3 Fe_tri+3 + 4 H2O = Fe_tri3(OH)4+5 + 4 H+ log_k -6.3 delta_h 14.3 kcalFe_tri+3 + Cl- = Fe_triCl+2 log_k 1.48 delta_h 5.6 kcalFe_tri+3 + 2 Cl- = Fe_triCl2+ log_k 2.13Fe_tri+3 + 3 Cl- = Fe_triCl3 log_k 1.

157、13Fe_tri+3 + SO4-2 = Fe_triSO4+ log_k 4.04 delta_h 3.91 kcalFe_tri+3 + HSO4- = Fe_triHSO4+2 log_k 2.48Fe_tri+3 + 2 SO4-2 = Fe_tri(SO4)2- log_k 5.38 delta_h 4.60 kcalFe_tri+3 + HPO4-2 = Fe_triHPO4+ log_k 5.43 delta_h 5.76 kcalFe_tri+3 + H2PO4- = Fe_triH2PO4+2 log_k 5.43Fe_tri+3 + F- = Fe_triF+2 log_k

158、 6.2 delta_h 2.7 kcalFe_tri+3 + 2 F- = Fe_triF2+ log_k 10.8 delta_h 4.8 kcalFe_tri+3 + 3 F- = Fe_triF3 log_k 14.0 delta_h 5.4 kcalPHASESGoethite Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O log_k -1.0ENDSOLUTION 1 pH 7.0 pe 10.0 O2(g) -0.67 Fe_di 0.144 / 102 Na 10. Cl 10. chargeEQUILIBRIUM_PHASES 1 O2(g) -0.

159、67RATESFe_di_ox-start10 Fe_di = TOT(Fe_di)20 if (Fe_di = 0) then goto 20030 p_o2 = 10(SI(O2(g)40 moles = (2.91e-9 + 1.33e12 * (ACT(OH-)2 * p_o2) * Fe_di * TIME200 SAVE moles-endKINETICS 1Fe_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.728e51.72

160、8e5INCREMENTAL_REACTIONS trueSELECTED_OUTPUT -file ex9.sel -reset falseUSER_PUNCH-headings Days Fe(2) Fe(3) pH si_goethite10PUNCHSIM_TIME/3600/24TOT(Fe_di)*1e6,TOT(Fe_tri)*1e6,-LA(H+),SI(Goethite)END这个 SOLUTIONSOLUTION 数据块定义了具有 0.1mmol/kgw 铁离子Fe-di的氯化钠溶液， 以与它与大气中的氧达到平衡。这个EQUILIBRIUM_PHASESEQUILIBRIU

161、M_PHASES 相数据块指定了所有的多组反响也将与大气中的氧达到平衡；这样，对于铁离子的氧化而言保持了充足的氧供给。在 RATESRATES 数据块中，速率表达式指定的名字为“Fe_di_ox，与根据等式159 来定义。注意应用特定的 Basic 函数“TOT来取得铁离子 line 10的总浓度molality， “SI是得到其饱和指数，或，在气体情况下，气体分压力的对数氧， line 30，和“ACT是取得 OH-line 40的活度。Line 40 定义了反响的摩尔数。注意，变量moles 的计算是通过速率乘以当前的时间间隔TIME， 与速率的定义是以 SAVE 声明来完毕的，SAVE

162、和 TIME声明一定要包含于速率的定义中；它们指定了经过时间间隔sub-的反响摩尔数。在KINETICSKINETICS数据块中，名为的速率表达式被调用和其参数被定义。在KINETICSKINETICS数据块中，速率的名字等于定义在PHASESPHASES下的矿物的名字，矿物的计量将会应用在反响中。然而，没有矿物与这个例子中速率的名字相关，与标识符- formula formula一定用来指定反响计量，这个反响涉与到溶液中的Fe_di等于Fe(2)的损失， 这正如计量系数为-1.0所说明的一样。 在溶45 / 102液中损失平衡的取得是通过计量系数为1.0的Fe_tri等于Fe(3)。 注意方

163、程式所包含的元素即为体系中质量的变化。这样，这个例子整个的动力学反响是Fe2+H+ 0.25O2 = Fe3+0.5H2O，但是形成水的反响的质子和氧并不改变体系中总的氢和氧的质量。 氢和氧因此并不包含于其分子式中。在这个例子中，氧是通过与大气中的O2g相平衡来补充的，转移1摩尔的氧发生于来自于EQUILIBRIUM_PHASESEQUILIBRIUM_PHASES中O2g相到溶液中。在接近于氧的体系中，溶解氧将会局部地被消耗。KINETICSKINETICS数据块-steps-steps标识符给予的是动力学反响必须被积分的时间步。当INCREMENTAL_REACTIONS trueINCR

164、EMENTAL_REACTIONS true被用到时，每一次时间步将会增加到所模拟的总的时间上与来自于前面时间步的结果将会作为当前时间步的起点来使用。数据块指定了所选取的输出文件的名字并且消除了所有打印到文件中的缺省的打印值- reset reset为假。在这个例子中所选取输出文件的唯一输出结果是定义在USER_PUNCHUSER_PUNCH数据块中。在所用到的Basic程序指定了在每一次动力学反响之后的接着的打印结果-steps-steps定义了11动力学时间步：模拟的时间累积是以天表示；总铁和铁离子以mol/kgw表示；与pH；和针铁矿的饱和指数。当运行输入文件时，在积分期间将会产生两个警

165、告信息， 如果积分时间间隔太大， 动力学反响增量将会产生负的溶液浓度的初始估计是可能的。 当发生这种情况时， 这个程序会打印出警告信息，降低其时间步的大小，和重新开始积分。这些信息是警告信息，不是错误，这个程序能够顺利地完成模拟。 通过减低初始时间间隔来减少这些警告信息是可能的。 如果标识符-step_divide-step_divide 应用为 100KINETICSKINETICS ，这表示为把初始整个时间步为 100 平均，将不会打印警告信息。否那么，如果标识符-step_divide-step_divide 应用为 1e-7，那将不打印其警告信息，这将引起反响增量小于1e-7mol。前

166、一种方法，-step_divide-step_divide 100 通常是更可取的，因为，尽管初始反响增量编译起来是很小的，后面在其积分的反响增量是可能的。应用-step_divide-step_divide 1e-7 会致使反响增量在整个积分期间留在更小的围上，运行时间大约比使用-step_divide-step_divide 长 5 倍，如果根本没有使用时，运行时间大约为其10 倍。图 9 说明了总 Fe(2)，总 Fe(3)的浓度，与在反响容器中经过去 10 天模拟反响的 pH。在反响开始时，能够看到 pH 值迅速的降低。Fe(2)随时间性的下滑是非常快的， 但随着反响的进展， 将逐渐的

167、松懈， 这与等式 159 相一致。 当在实际情况下的非缓冲溶液中操作实验时，需注意的是 pH 最初是升高的。pH 的这种提升将与缓慢形成的 Fe(3)络合物相一致的。因为氧化复原反响的自身消耗质子， 如果氢氧根络合物降低了pH 缓慢形成， 那 pH 最初是升高的。这样水中络合物的动力学信息最初也包含于PHREEQC 的模拟中，但它需要Fe(3)的氢氧根络合物也使用独立的 SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES 来定义，与其定义的络合物动力学信息的速率表达式。46 / 102图 9总 Fe(2)， 总 Fe(3)的浓度， 和 pH 随亚铁离子

168、Fe(2)氧化成三价铁离子Fe(3)的变化。例例 1010文石和菱锶矿的固体溶液文石和菱锶矿的固体溶液PHREEQC有能力模拟理想的、多组分的或是非理想的二元固体溶液。对理想固体溶液而言，每个最后一种固体成分的活度等于它的摩尔分数。 对非理想的固体溶液而言， 每个最后一种固体成分的活度是摩尔分数和活度系数的乘积， 这是由摩尔分数和Guggenheim过量自由能参数所决定的。下面的例子考虑了文石CaCO3-菱锶矿SrCO3的固体溶液，并证明了将碳酸锶参加到初始的纯相碳酸钙系统中后固体溶液和水相的组分发生的变化。这个例子来自于Glynn和Parkhurst (1992)所提出的曲线图。平衡常数为在

169、25时，KSrCO3109.271和KCaCO3108.336, 和Guggenheim参数,0 3.43与1 1.82(Plummer和Busenberg，1987)。输入数据见表29。PHASESPHASES数据块定义了文石和菱锶矿的logKs，并且覆盖了任何可能出现在数据库文件中的这些矿物的数据。SOLID_SOLUTIONSSOLID_SOLUTIONS数据块定义了无单位的Guggenheim过量自由能参数，与初始固体溶液组47 / 102分，其中文石和菱锶矿都为0摩尔。初始溶液1定义了碳酸氢钙溶液。然后， 这个溶液在二氧化碳分压大约为1大气压下与文石达到平衡，并作为新的溶液1存贮。图

170、 10. A固体溶液中菱锶矿和文石的摩尔分数， B液相中钙和锶的摩尔分数，C固体溶液中菱锶矿和文石的摩尔数，和D固体溶液中可混合区间的摩尔数，是碳酸锶增量的函数。虚线表示在可混合区间中的组分。在下一次模拟中，溶液1将与固体溶液相混合USEUSE关键字， 分500步参加5millimoles的碳酸锶REACTIONREACTION数据块。PRINTPRINT关键字数据块排除了所有缺省的打印到输出文件中，仅打印在USER_PRINTUSER_PRINT数据块所定义的打印项。USER_PRINTUSER_PRINT数据块指定了在每次反响后将要打印到输出文件中的固体溶液的信息：模拟的序号，反响步骤编号

171、，所参加碳酸锶的量，log，菱锶矿和文石的摩尔分数，液相中钙和锶的摩尔分数，和存在于混和性区间的48 / 102两种固体溶液的组分。这个SELECTED_OUTPUTSELECTED_OUTPUT数据块定义了输出文件为ex10.sel，取消了任何打印到所选取的输出文件中的缺省项 -reset-reset为假，并要求每一步中所增加的反响的数量定义在REACTIONREACTION数据块中写到所选取的输出文件中-reaction-reaction为真。USER_PUNCHUSER_PUNCH数据块打印到选择输出文件中的附加信息栏包括： 绘制图10所需要的所有的信息。 两次额外的模拟需要连续地把大量

172、的碳酸锶加到系统中，总的增加量为10摩尔。表 29例 10 输入数据的设置TITLE Example 10.-Solid solution of strontianite and aragonite.PHASES Strontianite SrCO3 = CO3-2 + Sr+2 log_k -9.271 Aragonite CaCO3 = CO3-2 + Ca+2 log_k -8.336ENDSOLID_SOLUTIONS 1 Ca(x)Sr(1-x)CO3 -comp1 Aragonite 0 -comp2 Strontianite 0 -Gugg_nondim 3.43 -1.82EN

173、DSOLUTION 1 -units mmol/kgw pH 5.93 charge Ca 3.932 C 7.864EQUILIBRIUM_PHASES 1 CO2(g) -0.01265 10 AragoniteSAVE solution 1END# Total of 0.00001 to 0.005 moles of SrCO3 added#USE solution 1USE solid_solution 1REACTION 1 SrCO3 1.0 .005 in 500 stepsPRINT -reset false -user_print trueUSER_PRINT49 / 102

174、-start 10 sum = (S_S(Strontianite) + S_S(Aragonite) 20 if sum = 0 THEN GOTO 110 30 xb = S_S(Strontianite)/sum 40 xc = S_S(Aragonite)/sum 50 PRINT Simulation number: , SIM_NO 60 PRINT Reaction step number: , STEP_NO 70 PRINT SrCO3 added: , RXN 80 PRINT Log Sigma pi: , LOG10 (ACT(CO3-2) * (ACT(Ca+2) +

175、 ACT(Sr+2) 90 PRINT XAragonite: , xc 100 PRINT XStrontianite: , xb 110 PRINT XCa: , TOT(Ca)/(TOT(Ca) + TOT(Sr) 120 PRINT XSr: , TOT(Sr)/(TOT(Ca) + TOT(Sr) 130 PRINT Misc 1: , MISC1(Ca(x)Sr(1-x)CO3) 140 PRINT Misc 2: , MISC2(Ca(x)Sr(1-x)CO3)-endSELECTED_OUTPUT -file ex10.sel -reset false -reaction tr

176、ueUSER_PUNCH-headlg_SigmaPi X_Arag X_Stront X_Ca_aq X_Sr_aq mol_Misc1 mol_Misc2mol_Arag mol_Stront-start 10 sum = (S_S(Strontianite) + S_S(Aragonite) 20 if sum = 0 THEN GOTO 60 30 xb = S_S(Strontianite)/(S_S(Strontianite) + S_S(Aragonite) 40 xc = S_S(Aragonite)/(S_S(Strontianite) + S_S(Aragonite) 50

177、 REM Sigma Pi 60 PUNCH LOG10(ACT(CO3-2) * (ACT(Ca+2) + ACT(Sr+2) 70 PUNCH xc # Mole fraction aragonite 80 PUNCH xb # Mole fraction strontianite 90 PUNCH TOT(Ca)/(TOT(Ca) + TOT(Sr) # Mole aqueous calcium 100 PUNCH TOT(Sr)/(TOT(Ca) + TOT(Sr) # Mole aqueous strontium 110 x1 = MISC1(Ca(x)Sr(1-x)CO3) 120

178、 x2 = MISC2(Ca(x)Sr(1-x)CO3) 130 if (xb x2) THEN GOTO 250 140 nc = S_S(Aragonite) 150 nb = S_S(Strontianite) 160 mol2 = (x1 - 1)/x1)*nb + nc 170 mol2 = mol2 / ( (x1 -1)/x1)*x2 + (1 - x2) 180 mol1 = (nb - mol2*x2)/x1 190 REM #Moles of misc. end members if in gap 200 PUNCH mol150 / 102 210 PUNCH mol2

179、220 GOTO 300 250 REM # Moles of misc. end members if notin gap 260 PUNCH 1e-10 270 PUNCH 1e-10 300 PUNCH S_S(Aragonite) # Moles aragonite 310 PUNCH S_S(Strontianite) # Moles Strontianite-endEND# Total of 0.001 to 0.1 moles of SrCO3 added#USE solution 1USE solid_solution 1REACTION 1 SrCO3 1.0 .1 in 1

180、00 stepsEND# Total of 0.1 to 10 moles of SrCO3 added#USE solution 1USE solid_solution 1REACTION 1 SrCO3 1.0 10.0 in 100 stepsEND过量的自由能参数描述了一种具有可混和性区间的非理想固体溶液。 对于落在可混合性区间的组分，液相中钙和锶的活度都将保持为常数， 且与两种组分的固体达到平衡， 一种固体中锶的摩尔分数为0.0048和另一种固体锶的摩尔分数为0.8579。 对于这个例子中的模拟而言，碳酸锶的每一次增量都会增加固体中碳酸锶的摩尔分数，直到参加其中的碳酸锶达到0.001

181、摩尔图10A。那一点是可混合性区间的开始图 10和固体组分中锶的摩尔分数是0.0048。下一次碳酸锶的增量达到所加到的碳酸锶为0.005mol在溶液中图10B产生钙和锶的常摩尔分数，并与可合混性区间的两个端点达到平衡。 然而，在固体相中碳酸钙和碳酸锶的数量图10C以与每个可混合性区间的端点的量 图10D随着碳酸锶的增加的而变化。最后，大约0.005mol的碳酸锶参加后达到可混合性区间的终点。在此处，这种溶液与锶的摩尔分数为0.8579的单一固体处于平衡。 碳酸锶的增加提高了液相中锶的摩尔分数， 直到增加了10摩尔的碳酸锶之后，固体溶液中的两种摩尔分数接近于1.0。51 / 102例例 1111

182、运移和阳离子交换运移和阳离子交换在存在阳离子交换的情况下，以下有关平流运移的例子来自于PHREEQM程序的简单计算Appelo和Postma，例10.13, p. 431-434。模拟了从含有阳离子交换剂的柱子中流出物的化学组分。最初， 柱子中是钠-钾-氮溶液，并与阳离子交换剂达到平衡。然后用三倍空隙体积的氯化钙溶液冲洗柱子。钙、钾和钠持续与交换剂反响， 并达到平衡。这个问题可以用两种方法运行，即使用ADVECTIONADVECTION数据块，仅模拟平流；使用TRANSPORTTRANSPORT数据块，将模拟平流和弥散相混合。输入数据组见表30，柱子中有40个单元，与Appelo和Postma

183、 (1993)所描述的一次运行相一致。这需要定义40种溶液，序号为1到40，溶液的序号与柱子中的单元对应。在这个例子中，所有的单元包含的溶液一样，但这不是必要的。每个单元中溶液的定义可以不同，也可根据当前或前面模拟中的反响定义 使用SAVESAVE关键字。每一单元中溶液的定义是强制性的，但每一单元换剂的定义是可选择的。 交换剂的序号与柱子中的单元号对应， 如果交换剂用单元号定义，它用于该单元的计算。在这个例子中，所有的单元使用同样的交换剂。SOLUTIONSOLUTION 1-40数据块定义柱子中40个单元中填充的溶液。柱子的填充溶液必须以SOLUTIONSOLUTION 0来进展定义，它是一

184、种氯化钙溶液。在这40个单元中，交换剂的量和组分是以EXCHANGEEXCHANGE 1-40数据块来进展定义的。在每一单元换点的数目为 1.1mmol，计算交换剂的初始组分，使其与溶液1相平衡。注意初始交换成分的计算是假定溶液1的组分是固定的；溶液1的组分在初始交换组分计算期间不发生变化。ADVECTIONADVECTION数据块仅需要包含单元号和模拟转变号。计算仅解决了解流入单元中孔隙体积的数目。没有用到时间和距离的显定义。标识符-punch_cells-punch_cells和-punch_frequency-punch_frequency指定了在每次转换中40单元的数据将来写到所选取的

185、输出文件中。-punch_cells-punch_cells和-punch_frequency-punch_frequency标识符说明40单元每二十次转换的数据将写到输出文件中。数据块指定了转换平流步数和钠、氯、钾与钙的总溶解浓度，它们将会被写到文件ex11adv.sel中。孔隙的容积可以通过转换的数目进展计算；一次转换移动溶液到下一个单元中，最后的溶液将会溢出柱子。PHREEQC计算了单元中心的浓度，因此，最后单元的浓度在柱子中的终端将会达到半次转换。 在这个例子中， 一次转换代表了1/40的柱子的空隙体积。因此，从柱子中流入的孔隙的数目PVPV= (转换的次数+ 0.5) / 40。使用

186、USER_PUNCHUSER_PUNCH数据块计算孔隙的容积并打印到所选取的输出文件中。表表3030例例1111输入数据的设置输入数据的设置TITLE Example 11.-Transport and ion exchange.SOLUTION 0 CaCl2 units mmol/kgw temp 25.0 pH 7.0 charge52 / 102 pe 12.5 O2(g) -0.68 Ca 0.6 Cl 1.2SOLUTION 1-40 Initial solution for column units mmol/kgw temp 25.0 pH 7.0 charge pe 12.5

187、 O2(g) -0.68 Na 1.0 K 0.2 N(5) 1.2EXCHANGE 1-40 equilibrate 1 X 0.0011ADVECTION -cells 40 -shifts 120 -punch_cells 40 -punch_frequency 1 -print_cells 40 -print_frequency 20SELECTED_OUTPUT -file ex11adv.sel -reset false -step -totals Na Cl K CaUSER_PUNCH -heading Pore_vol 10 PUNCH (STEP_NO + .5) / 40

188、.ENDSOLUTION 1-40 Initial solution for column units mmol/kgw temp 25.0 pH 7.0 charge pe 12.5 O2(g) -0.68 Na 1.0 K 0.2 N(5) 1.2EXCHANGE 1-40 equilibrate 1 X 0.0011TRANSPORT -cells 40 -length 0.00253 / 102 -shifts 120 -time_step 720.0 -flow_direction forward -boundary_cond flux flux -diffc 0.0e-9 -dis

189、persivity 0.002 -correct_disp true -punch 40 -punch_frequency 1 -print 40 -print_frequency 20SELECTED_OUTPUT -file ex11trn.sel -reset false -step -totals Na Cl K CaEND图 11. 在离子交换器中参加氯化钙溶液来置换钠和钾的运移模拟结果。线代表在 PHREEQC 中仅有平流ADVECTIONADVECTION 关键字和同时有平流与弥散TRANSPORTTRANSPORT 关键字时所进展计算的柱子出口处的浓度。平流计算ADVECTIO

190、NADVECTION之后，初始条件用第二组数据块SOLUTIONSOLUTION和EXCHANGEEXCHANGE设置为平流和弥散计算TRANSPORTTRANSPORT。根据ADVECTIONADVECTION模拟，SOLUTIONSOLUTION 0不改变，所以不需要重新定义。TRANSPORTTRANSPORT数据块相比ADVECTIONADVECTION数据块有更为明显的运移过程的描述。每个单元的54 / 102长度 -length-length ， 柱子末端的边界条件 -boundary_cond-boundary_cond ， 流动方向 -flow_direction-flow_d

191、irection ，弥散度-dispersivity-dispersivity，以与扩散系数-diffc-diffc都能够被指定。当用流动边界条件模拟柱子中流出物时， -correct_disp-correct_disp标识符应设置为真。 -punch-punch, -punch_frequency-punch_frequency, -print-print,和-print_frequency-print_frequency标识符在ADVECTIONADVECTION数据块中起一样的作用。第二个SELECTED_OUTPUTSELECTED_OUTPUT数据块指定了运移步转换的数目，以与将被写

192、到ex11trn.sel文件中的钠、氯、钾和钙的总溶解浓度。 来自于平流模拟的USER_PUNCHUSER_PUNCH数据块仍旧有效， 每次运移步长的孔隙体积也会得以计算，并写到选取的输出文件中。使用ADVECTIONADVECTION and TRANSPORTTRANSPORT关键字计算的例11的结果可用曲线表示，见图11。40单元的浓度，即末端单元，根据孔隙容积绘图。在这两种运移模拟中的计算具有一样的特征。氯化物是一种守恒溶液，在大约1个孔隙的体积时达到流出物。柱子中最初存在的钠与参加的钙发生交换，并被洗堤，只要交换器中含有钠。钠离子临界曲线的中点出现在大约1.5个孔隙体积处。由于钾的交

193、换能比钠的交换能更强在交换反响中 logK的值更大一些， 钾将在钠交换完毕才得以释放。最后，当所有的钾都释放出后，钙的浓度升高至一个稳定的值，此时等于填充溶液的浓度。在流出物中钾和钠的浓度变化形成了色谱法模式，常可以用简单的方法进展计算Appelo, 1994b。达到钠降低锋所需的孔隙体积数可用公式Pv1VS进展计算，这里Vs=q/c，q是吸附浓度的变化mol/kgw，c是溶质浓度对于锋的变化。最初填充柱子的溶液中钠的浓度是1.0mmol/kgw，且钠的初始吸附浓度为0.55；在填充溶液中钠的浓度为0，这最终导致吸附钠为0。这样，(Vs)1 q/c (0.550)/(10) 0.55，且Pv

194、1.55,这说明了钠锋的中点应当在达到1.55孔隙容积后出现于柱子的末端。下一步，从交换器中置换钾。溶液中的浓度升至1.2 mmol/kgw以平衡Cl 的浓度，当交-换物质消耗尽时，便降低到0。当钾为溶液中仅有的阳离子时，它也是交换剂中仅有的阳离子。对钾而言，(V )2 q/c (1.10)/(1.20) 0.917，且Pv=1.917孔隙容积，可以看出， (V)1和(V)2的锋面位置与列于图11中浓度变化的中点严密匹配。这两次模拟的差异完全是由于在TRANSPORTTRANSPORT计算中包括了弥散。在TRANSPORTTRANSPORT计算中，氯化物的临界曲线与守恒溶质的平流弥散方程的分析

195、溶液一致Appelo and Postma, 1993,p.433。没有弥散，ADVECTIONADVECTION计算产生的氯化物的临界曲线为方波。 在其他元素的锋面计算中也缺少弥散的拖尾效应特性，尽管由于交换反响的影响也存在着一定的曲率。在ADVECTIONADVECTION模拟中钾的最大浓度相对要大一些，因为无视了弥散的作用。sss例例 1212热与溶质的平流弥散流动热与溶质的平流弥散流动下面这个例子证明了 PHREEQC 软件用来计算一维环境中的热与溶质的运移。 这一列的开55 / 102始最初是为冲淡的 KCl 溶液，这种溶液是在 25下，与阳离子达到平衡条件中的。KNO3 溶液然后平

196、流进入到这一栏中， 建立新的 0。后来是，在 24的氯化钠溶液将允许从这一栏的两端进展弥散，假定在通过这一栏时， 没有热的损失。在一端施加的是常边界的情况，而在另一端最后的剖分区域那么充满的是氯化钙溶液， 且假定的是封闭边界条件。 对具有常边界条件的一端而言， 所分析的溶液与 PHREEQC 的结果进展比照， 这适用于延迟R= 1.0 (Cl-)和R= 3.0 (Na+和温度)。最后，数字化方法的第二次序的准确度是通过增加剖分区域的数目来确定的，这是通过三个因素来确定的，同时，证明了这些数目的溶液出错情况的降低，这是通过与所分析溶液比照而具有近似的数量级得出的。表 31例 12 输入数据的设置

197、TITLE Example 12.-Advective and diffusive transport of heat and solutes. Constant boundary condition at one end, closed at other. The problem is designed so that temperature should equal Na-conc (in mmol/kgw) after diffusion.EXCHANGE_SPECIES Na+ + X- = NaX log_k 0.0 -gamma 4.0 0.075 H+ + X- = HX log

198、_k -99. -gamma 9.0 0.0 K+ + X- = KX log_k 0.0 -gamma 3.5 0.015SOLUTION 0 24.0 mM KNO3 units mol/kgw temp 0 # Incoming solution 0C pH 7.0 pe 12.0 O2(g) -0.67 K 24.e-3 N(5) 24.e-3SOLUTION 1-20 0.001 mM KCl units mol/kgw temp 25 # Column is at 25C pH 7.0 pe 12.0 O2(g) -0.67 K 1e-6 Cl 1e-6EXCHANGE 1-20

199、KX 0.048TRANSPORT # Make column temperature 0C, displace Cl -cells 20 -shifts 1956 / 102 -flow_d forward -bcon flux flux -length 1.0 -disp 0.0 # No dispersion -diffc 0.0 # No diffusion -thermal_diffusion 1.0 # No retardation for heatPRINT -reset falseENDSOLUTION 0 Fixed temp 24C, and NaCl conc (firs

200、t type boundary cond) at inlet units mol/kgw temp 24 pH 7.0 pe 12.0 O2(g) -0.67 Na 24.e-3 Cl 24.e-3SOLUTION 20 Same as soln 0 in cell 20 at closed column end (second type boundarycond) units mol/kgw temp 24 pH 7.0 pe 12.0 O2(g) -0.67 Na 24.e-3 Cl 24.e-3EXCHANGE 20 NaX 0.048TRANSPORT # Diffuse 24C, N

201、aCl solution from column end -shifts 1 -flow_d diffusion -bcon constant closed -thermal_diffusion 3.0 # heat is retarded equal to Na -diffc 0.3e-9 # m2/s -timest 1.0e+10 # 317 years, 19 substeps will be usedSELECTED_OUTPUT -file ex12.sel -high_precision true -reset false -dist true -temp trueUSER_PU

202、NCH -head Na_mmol K_mmol Cl_mmol 10 PUNCH TOT(Na)*1000, TOT(K)*1000, TOT(Cl)*1000END57 / 102第12个例子中输入数据的设置列于表31中。应用了EXCHANGE_SPECIESEXCHANGE_SPECIES数据块1使KX的交换常数等于KNa的log_k = 0.0，2有效且尽可能地移除氢离子的交换离子，与3设置交换种类的交换系数等于他们水离子的配对物-gamma-gamma标识符，因此，Na+ 和K+之间的交换是线性的，其延迟为常数。运移的注入溶液，即溶液0，它的定义在在0下的24 mmol/kgw KN

203、O3。为了使pe稳定，氧浓度的定义是与大气中氧的分压力相平衡。这一栏离散为20个剖分区域，这最初注入的是为a1mol/kgw KCl 溶液(SOLUTIONSOLUTION 1-20)。每一剖分区域都具有48mmol的交换位置，它的定义仅包含了EXCHANGEEXCHANGE 1-20数据块中的钾。TRANSPORTTRANSPORT数据块定义了1米长的区域长度-length-length 1，没有扩散-disp-disp 0，没有弥散-diffc-diffc 0，且没有温度的延迟-thermal_diffusion-thermal_diffusion 1。SOLUTIONSOLUTION 0

204、转换了19次进入了这一栏中-shifts-shifts 19, -flow_d-flow_d向前，在最后的平流弥散运移模拟中到达了区域19中，这一栏两端的边界条件都是流动的类型-bcon-bcon flux flux。在这个初始平流-弥散运移模拟中，没有交换发生溶液仅包含有钾阳离子，交换点完全由钾所填满，在 19区域中开始物质的浓度和温度分别为24 mmol/kgw的 KNO3和0oC。这个结果能够用数据块更为直接容易的取得如附录C，这个模拟证明了多么短暂的边界和流动条件能够被描绘。在平流-弥散运移模拟中，每一区域中的溶解物质、固体组分和温度能够在每一次转换之后自动地保存，PRINTPRINT

205、关键字用来排除所有的打印到输出文件中-reset-reset为假。在下面的平流-弥散运移模拟中， 弥散从这一栏的两端开始计算， 以这一栏的组分开始，并且温度为存在于下面的第一次平流扩散运移模拟中， 这当然是除了现在的NaCl溶液是作为溶液0来定义的，它将扩散到这一栏的上部，至于溶液20，那么是从这一栏的底部开始弥散的。新SOLUTIONSOLUTION 0的定义是在24oC下的24 mmol/kgw NaCl。最后的区域SOLUTIONSOLUTION 20也定义为具有这种溶液的组分和温度。区域20的交换物质定义为与日俱增0区域的新的溶液组分相平衡EXCHANGEEXCHANGE 20。TRA

206、NSPORT数据块定义了一段弥散运移时期-shifts-shifts 1; -flow_d-flow_d扩散。在第一区域的边界条件为常浓度，且在最后的区域这一栏是封闭的-bcon-bcon常封闭。有效的平流扩散系数设置-diffc-diffc为0.3e-9 m2/s，时间步-timest-timest定义为1.e10 s。然而，这个时间步将会自动地化分为一系列的时间子步， 以满足数字化方法的稳定的标准。 热延迟因素设为3.0-thermal_diffusion-thermal_diffusion 3.0。对Na+而言的可交换的浓度最大NaX为48 mmol/kgw与溶质浓度最大a+ =24 m

207、mol/kgw相对所有浓度的比率为2.0，因此，这个延迟为，它在数字上等于温度的延迟。SELECTED_OUTPUTSELECTED_OUTPUT数据块指定了所选择输出文件的名字为ex12.sel。标识符“-high_precision-high_precision true是用来获得在打印中所增加的阿拉伯数字的序号，而标识符-dist-dist和-temp-temp指定了将会打印到文件中的每一区域的距离和温度。USER_PUNCHUSER_PUNCH数据块是用来打印钠、钾和氯的浓度到所选择的输出文件中，单位为mmol/kgw。58 / 102图12热和钠离子延迟R=3与Cl-R=1柱子末端弥

208、散的结果与常边界条件分析溶液的比照在模拟中，温度的计算是以线性延迟分子式进展的；然而，钠离子的浓度的计算是通过阳离子交换反响进展的。 即使钠离子的浓度和温度的计算是以不同的方法进展的， 但在数值上是应该相等的， 这是因为初始条件和瞬时条件在数值上是一样的， 并且延迟因素也是一样的。在PHREEQC所选择输出文件中， 温度和钠离子浓度至少有6位是一样的，这说明了在这个例子中考虑简单的化学情况， 其化学运移计算算法是正确的。 更进一步准确度的取得是通过与所分析溶液比照模拟结果来取得的。对明确的具有Cx = 0 fort= 0,和扩散从x= 0，以Cx=0 =C0 对t 0 而言的分析解决方法为：C

209、x,tC0erfc(x2 Det/ R)这里De是有效扩散系数。PHREEQC的结果是在图12中与Cl-的分析溶液和的Na+温度相比照，说明极好的一致性。注意，来自这一栏中Cl-的弥散仍没有“达到中间的局部，因此， 这一栏仍是有效明确的，并且分析溶液是适当的。尽管这一栏的两端都是以一样的温度和浓度开始，x=0仍是由于常边界条件而在一样的温度和浓度条件下取得的。 温度和浓度已在20区域中降低 划分为区域的中点，x= 19.5，这是因为封闭的边界条件适用于x= 20；没有热和质量通过这个边界是允许的， 并且温度和浓度由于弥散到这一栏中而导致降低。 钠的浓度不如氯的浓度消耗的那样快，这是由于交换点一

210、定在沿扩散路径注入。由于这个例子是分析溶液，在PHREEQC中使用的数字化算法中确定第二次序的准确性是可行的。 对第二种方法而言， 通过三个因素之一降低剖分区域的大小应当能够改善结果中的九个因素之一。在附录C中所给出的输入文件执行了在这个例子中的20个剖分区域和60剖分区域的计算。在每一时间步完毕时的分析溶液偏离的计算，其距离是从0.5到8.5米以0.5米增加的。这些结果列于表32中。第二种方像所期望的一样， 分析溶液的偏离的降低大约是作59 / 102为降低结果的三个因素之一产生的。表表3232例例1212的的2020区域和区域和6060区域模拟解析解法的数字误差区域模拟解析解法的数字误差距

211、离距离ClCl浓度的误差浓度的误差2020区域模拟区域模拟0.51.52.53.54.55.56.57.58.53.32E-058.17E-059.18E-057.15 E-054.24 E-052.00 E-057.81 E-052.55 E-055.58 E-056060区域模拟区域模拟3.03 E-067.66 E-069.09 E-067.65 E-064.98 E-062.61 E-061.12 E-063.97 E-077.65 E-08NaNa浓度的误差浓度的误差2020区域模拟区域模拟5.75 E-045.54 E-048.29 E-05-5.07 E-05-2.54 E-05

212、-5.44 E-06-7.20 E-07-6.77 E-08-4.90 E-096060区域模拟区域模拟4.42 E-056.08 E-051.43 E-05-5.64 E-06-3.26 E-06-6.27 E-07-6.15 E-08-3.48 E-09-1.21 E-10？例？例 1313阳离子交换双多孔柱的阳离子交换双多孔柱的 1 1 维运移维运移这个例子证明了PHREEQC计算双孔隙度媒介中发生交换反响物质的流动的能力，交换反响通过弥散在流动的孔和不流动孔。PHREEQC中输入格式灵活、附加溶液和过程的模块化定义允许在1维柱中包含的不同成分和各种各样的配合物。 该例子使用直径为2厘米

213、柱子， 且充满了小粘土颗粒中充当滞留区。 该例应用了一次交换的近似值和有限差分， 并分别考虑了守恒和延迟通过离子交换的化学运移。 它也更进一步说明了能够通过指定混合系数， 用关键字MIXMIX来解释流动单元和不流动单元的弥散来模拟不同种类的柱子。使用具有在混合系数的一使用具有在混合系数的一次交换近似值进展停滞区计算交换近似值进展停滞区计算这个例子的输入文件表33适用于停滞多孔沿柱统一分布的柱子。 这20个活动单元的序号为1-20。每一个活动单元，n，与一个不活动的单元发生交换，其序号为20+1+n单元22-41是不活动单元。所有的单元给定一样的初始溶液和交换配合物，但每个单元的定义可以不同。也

214、可能不统一分配柱子中的不流动单元，仅仅忽略掉不存在的滞留单元的溶液。流动带和滞流带之间的联系和滞流带单元之间的关系沿柱子可以变化，这需要用关键字MIXMIX来指定流动和非流动单元之间混合的系数。正如在表33中所定义的一样，在流动和滞流带溶液 1-41的柱子中最初包含1 mmol/LKNO3溶液。NaCl-NO3流进到柱子中溶液0。每一单元定义1mmol吸附点的交换配合物交换1-41，并且交换系数也适用于Na EXCHANGE_SPECIESEXCHANGE_SPECIES交换种类的线性延迟系数R=2。第一个TRANSPORTTRANSPORT数据块是用来定义第一次运移模拟的物理和流动特性。 柱

215、子长2米， 剖分为60 / 102+20个0.1m长-length-length的单元-cells-cells。溶液注入柱子中为溶液0一个脉冲5一次转换 -shifts-shifts 。 每次转换时间为3600秒 -timest-timest ， 这使得在流动空隙中的速率为：vm= 0.1/ 3600 = 2.78e-5 m/s。所有单元的弥散度(-disp-disp)设置为0.015米。扩散系系设置为0.0-diffc-diffc。使用一阶交换近似值定义滞流/流动交换。根据表1见“在双孔介质中的运移，滞流单元由半径为r=0.01m的球体组成，扩散系数De= 3.e-10 m2/s，形状参数f

216、S1 0.21。变量所给定的交换因素为= 6.8e-6s-1。流动的孔隙率为m = 0.3 (-stag-stag) 与不流动的孔隙率为= 0.1。对PHREEQC中的一阶交换近似值而言，每个单元的不流动区和其参数,、m和im都是以指定的。滞流单元的定义会导致每一剖分单元的流动单元序号1-20都与相关单元的不流动单元相关 序号为22-41。数据块用来消除所有的打印到输出文件中。下面连续的注入NaCl溶液，以与10次转换1mmol的KNO3/L第二个溶液0都将参加到这一栏中。第二个TRANSPORTTRANSPORT数据块没有重新定义任何的柱子或是流动的特征，但却指定了1到20单元(-punch

217、_cells-punch_cells)的结果在10次转换(-punch_frequency-punch_frequency)之后都将写到所选取的输出文件中。数据块SELECTED_OUTPUTSELECTED_OUTPUT和USER_PUNCHUSER_PUNCH指定了将要写到输出文件中的数据。表表3333例例13A13A输入数据的设置输入数据的设置TITLE Example 13A.-1 mmol/l NaCl/NO3 enters column with stagnant zones. Implicit definition of first-order exchange model.SO

218、LUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, stagnant exchange N(5) 1.0 # NO3 is conservative# charge imbalance is no problem .ENDSOLUTION 1-41 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0

219、EXCHANGE 1-41 -equil 1 X 1.e-361 / 102EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX K+ + X- = KX log_k 0.0 -gamma 3.5 0.015ENDTRANSPORT -cells 20 -shifts 5 -flow_d forward -timest 3600 -bcon flux flux -diffc 0.0 -length 0.1 -disp 0.015 -stag 1 6.8e-6 0.3 0.1# 1 stagnant l

220、ayer, alpha, theta(m), theta(im)PRINT -reset falseENDSOLUTION 0 # Original solution reenters units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-20SELECTED_OUTPUT -file ex13a.sel -reset false -solution -distance trueUSER_PUNCH -head Cl_mmol Na

221、_mmol10 PUNCH TOT(Cl)*1000, TOT(Na)*1000END这个例子的一阶交换近似值的混合因素mixfm和mixfim都是来自于下面的等式 121 和123：mixfimmt(mim)(1exp()和mixfm mixfimimmimmimm延迟因素Rm和Rim都没有包含在这里的mixfim和mixfm公式中，这是因为在PHREEQC中延迟是62 / 102化学反响的结果。 根据等式161和162， 这个例子的混合因素计算为mixfim= 0.20886和mixfm= 0.06962。流动单元和不流动单元的混合因素是不同的，它们代表了流动的水和不流动的水体积的不同。在

222、PHREEQC中， 流动的和滞流水混和是在每次扩散/弥散步完成后才完成的。 这意味着当剖分单元变小或是划分了更多的扩散步“mixruns执行时时间步会变小。在例13中20个单元模拟完成了一个混和运行。100个单元将会完成有3个混和运行等式110要求mixf1/3的n=3，mixfim计算的时间步将会为(3600/5) / 3 = 240 s。时间步t= 240 s将会导致在100单元模拟中mixfim= 0.01614。利用具有显性混合因素的一阶交换近似值计算滞流单元利用具有显性混合因素的一阶交换近似值计算滞流单元滞流单元统一分布的显性混合因素的输入文件是在表34中所给出的。 SOLUTION

223、SOLUTION数据块等于前面的输入文件。没有进一步信息的滞流层将会被得以定义-stag-stag 1, 在TRANSPORTTRANSPORT数据块中。流动/非流动的交换是由数据块所给定的因素进展设置的。这个输入文件的结果与前面所使用的具有简单定义概念的输入文件是一样的。 然而， 混合因素的显性定义证明了滞流单元的非正式分布， 或是滞流单元的其它物理特性， 它们都能够通过定义流动和非流动区交换的混合因素来包含于PHREEQC的模拟中。表表3434例例13B13B输入数据的设置输入数据的设置TITLE Example 13B.-1 mmol/l NaCl/NO3 enters column w

224、ith stagnant zones. Explicit definition of first-order exchange factors.SOLUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, stagnant exchange N(5) 1.0 # NO3 is conservative# charge imbalance is no problem .ENDSOLUTION 1-4

225、1 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0EXCHANGE 1-41 -equil 1 X 1.e-363 / 102EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX K+ + X- = KX log_k 0.0 -gamma 3.5 0.015ENDMIX1;1 .93038;22 .06962;MIX2;2 .93038;23 .06962;MIX3;3 .93038;24 .06962;

226、MIX4;4 .93038;25 .06962;MIX5;5 .93038;26 .0696227 .06962;MIX7;7 .93038;28 .0696229 .06962;MIX9;9 .93038;30 .0696231 .06962;MIX 11; 11 .93038; 32 .06962 ;MIX 12; 12MIX 13; 13 .93038; 34 .06962 ;MIX 14; 14MIX 15; 15 .93038; 36 .06962 ;MIX 16; 16MIX 17; 17 .93038; 38 .06962 ;MIX 18; 18MIX 19; 19 .93038

227、; 40 .06962 ;MIX 20; 20#MIX 22; 1 .20886; 22 .79114 ;MIX 23; 2MIX 24; 3 .20886; 24 .79114 ;MIX 25; 4MIX 26; 5 .20886; 26 .79114 ;MIX 27; 6MIX 28; 7 .20886; 28 .79114 ;MIX 29; 8MIX 30; 9 .20886; 30 .79114 ;MIX 31; 10MIX 32; 11 .20886; 32 .79114 ;MIX 33; 12MIX 34; 13 .20886; 34 .79114 ;MIX 35; 14MIX 3

228、6; 15 .20886; 36 .79114 ;MIX 37; 16MIX 38; 17 .20886; 38 .79114 ;MIX 39; 18MIX 40; 19 .20886; 40 .79114 ;MIX 41; 20TRANSPORT -cells 20 -shifts 5 -flow_d forward -timest 3600 -bcon flux flux -diffc 0.064 / 102;MIX6;MIX8;MIX 10;.93038; 33.93038; 35.93038; 37.93038; 39.93038; 41.20886; 23.20886; 25.208

229、86; 27.20886; 29.20886; 31.20886; 33.20886; 35.20886; 37.20886; 39.20886; 416 .93038;8 .93038;10 .93038;.06962;.06962;.06962;.06962;.06962;.79114;.79114;.79114;.79114;.79114;.79114;.79114;.79114;.79114;.79114; -length 0.1 -disp 0.015 -stag 1PRINT -reset falseENDSOLUTION 0 # Original solution reenter

230、s units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-20SELECTED_OUTPUT -file ex13b.sel -reset false -distance true -solutionUSER_PUNCH -head Cl_mmol Na_mmol10 PUNCH TOT(Cl)*1000, TOT(Na)*1000END应用有限差分近似值进展滞流单元的计算应用有限差分近似值进展滞流单元的计算滞流单元所包含的围为半径

231、为r= 0.01 m的圆周。扩散到这一圆周单元会致使迅速平衡C2C2C De()浓度的变化，这是根据下面的微分等式进展的：trrr因而在有限差分中的计算能够简化为一维放射状的情况。表表3535在圆周弥散中有限差分的混合因素在圆周弥散中有限差分的混合因素单元单元102826242221r r,m,m0.001.003.005.007.009-v vj j,m,m3 33.35e-82.35e-76.37e-71.24e-62.04e-61.26e-5A Aijij,m,m2 25.30e-52.01e-44.52e-48.04e-41.26e-3-h hijij,m,m0.002.002.002

232、.002.002-f fbcbc11111.72-mixmixf fijij0.81.463.384.350.571-mixmixf fijij0.19.421.446.453.217.907MixMixf fijij-0.116.170.197.212.09365 / 102这个计算是依赖于在这个手册 “在双孔隙度介质中的运移局部所概括的理论进展的。 在这个例子中，半径被划分为r= 0.002 m的五个等距的层。这五个滞流层的定义是在关键字TRANSPORTTRANSPORT中以-stagnant-stagnant 5来进展定义的表36 。在滞流单元中相邻单元间混合确实定是以MIXMIX数据

233、块进展的，混合因素是通过等式128和129进展的。对计算而言，单元cellj(m3)的体积Vj是需要的。共享的有单元celli和j(m2)的交接单元Aij，单元cellsi和j(m)中点之间的距离，与其边界单元fbc无量纲的修正因素。流动单元1和和相关的不流动单元的值都在表35中给出。在不流动层中单元给定的序号为n+lxcells+ 1，这里n流动单元的数目，l是滞流单元的数目，以与cells为流动单元的总数目。这个例子中，滞流单元边界单元的编号是从单元序号1到单元序号22， 与其它四个滞流单元的编号为42， 62， 82和102，编号102是这个围中的最里面的单元，它只与一个别的单元相关单元

234、82。流动单元的表示是相对于半径为0.01m的圆的围而言的，通过乘以m /im的体积(4.19e-6 0.3 / 0.1 =1.26e-5)。在表35中，fbc的值为1.72，正如等式127所计算的一样。 值得一提的是应用fbc=1.81明显地证明了Crank (1975)所给出的扩散到密闭容器中有溶液和粘土层的烧杯分析溶液的合理性。然而，改变fbc到1.81对柱子中的影响是很小的，这在图13中有所说明表表3636例例13C13C的输入数据的输入数据Table 36. -Input data set for example 13C: Stagnant zone with diffusion c

235、alculatedby finite differences (partial listing)TITLE Example 13C.-1 mmol/l NaCl/NO3 enters column with stagnant zones. 5 layer stagnant zone with finite differences.SOLUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, sta

236、gnant exchange N(5) 1.0 # NO3 is conservative# charge imbalance is no problem .ENDSOLUTION 1-121 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0EXCHANGE 1-121 -equil 1 X 1.e-3EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX66 / 102 K+ + X- = KX log_k

237、 0.0 -gamma 3.5 0.015ENDMIX 1; 1 0.90712; 22 0.09288MIX 22; 1 0.57098; 22 0.21656; 42 0.21246MIX 42; 22 0.35027; 42 0.45270; 62 0.19703MIX 62; 42 0.38368; 62 0.44579; 82 0.17053MIX 82; 62 0.46286; 82 0.42143; 102 0.11571MIX 102; 82 0.81000; 102 0.19000# MIX definitions omitted for mobile cells# 2-19

238、 and associated immobile cells#MIX 20; 20 0.90712; 41 0.09288MIX 41; 20 0.57098; 41 0.21656; 61 0.21246MIX 61; 41 0.35027; 61 0.45270; 81 0.19703MIX 81; 61 0.38368; 81 0.44579; 101 0.17053MIX 101; 81 0.46286; 101 0.42143; 121 0.11571MIX 121; 101 0.81000; 121 0.19000TRANSPORT -cells 20 -shifts 5 -flo

239、w_direction forward -timest 3600 -tempr 3.0 -bcond flux flux -diffc 0.0 -length 0.10 -disp 0.015 -stag 5PRINT -reset falseENDSOLUTION 0 # Original solution reenters units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-2067 / 102SELECTED_OUTPUT

240、-file ex13c.sel -reset false -solution -distance trueUSER_PUNCH -head Cl_mmol Na_mmol10 PUNCH TOT(Cl)*1000, TOT(Na)*1000END注意在表36中定义了121种溶液， 1-20是流动单元的溶液， 其他的为不流动单元的溶液。输入文件与前面的一样，只除了-stag 5和单元之间的混合因素。并不是所有的混合因素都列于表36中，对后面的单元和它们邻近的单元， 它们是相等的。在这个具有粘土颗粒的例子中，仅考虑到了1D的弥散，并且仅在不同层中的不流动单元被定义。然而，包含所有的邻近单元中不流动

241、单元中的混合是可能的。图13扩散至球形滞流单元的运移模拟的结果，运用有限差分和一价交换近似值进展模拟的。图13比照了在例13A和B与例13C中所得到的流动单元剖面。这两个模拟的根本的特征是一样的。峰值的位置，正像所计算的一样是相似的。这Cl-峰接近于1.2m，但在缺少滞流单元时将可能达到1.45m。在可变孔隙度的单元整个的浓度大约等于一阶交换和有限差分模拟。一阶交换近似值的交换因素fs1 = 0.21的出现提供了这个模拟的足够的准确度。然而， 一阶交换近似值产生了较低的峰值和与准确溶液相比用有限差分取得了更长的谱尾。 差68 / 102异和偏离同样也出现在曲线上Van Genuchten, 1

242、985。当应用模拟这一围的运移模拟时，一阶交换模拟可能至少是准确的； 滞流单元的其它形状能够给定一个更好的圆周。 因此线性的交换模型更容易适用，这是因为显性的模拟需要混合因素扩展序列的准备注意应用USER_PUNCHUSER_PUNCH的单独模拟能够达到这个目的 ，当采用离散化时它将会发生变化。 有限差分模拟的计算时间乘以不流动单元层也能够比拟的长， 这当然是对单个不流动单元层的一阶交换近似值而言的。？例？例 1414平流运移，平流运移，阳离子交换，阳离子交换，外表络合反响和矿物均外表络合反响和矿物均衡衡这个例子应用了 PHREEQC 中的相的平衡， 阳离子交换和溶液外表络合反响的能力， 这与

243、在平流运移中用来模拟在中部 Oklahoma 含水层中的水的演化是相比照的。含水层的地球化学特性已由 Parkhurst 和他人(1996)所描述。 在这一层中有两种主要的水的类型， 在含水层中承压水层中的重碳酸钙锰水的 pH 值的围为 7.0 到 7.5，且在含水层中的承压水局部重碳酸钠水的 pH 值的围在 8.5 到 9.2。另外，来自于海水中的氯化钠盐水存在于含水层之下，假定在含水层中为流动的含物和末端为有孔的空间。砷几乎特定地是与高pH 值的重碳酸钠水的类型相关。这个例子计算的概念模型假定盐水最初是灌注到含水层中去的。这个含水层包含有钙，白云石，具有阳离子交换能力的粘土，和氧化水铁离子

244、的外表；最初， 阳离子的交换和溶液的外表是与盐水相平衡。 含水层假定是由雨水重新灌注， 同时通过蒸发和在渗流区通过与钙和白云石相平衡而达到的。 这些水然后进入到饱和区与出现阳离子交换和水相氧化铁离子外表的钙和白云石相反响。这个计算使用 PHREEQC 软件的平流运移能力， 它是由仅由单个的剖分区域来代表饱和区域。总共有 200 个的孔隙区域灌注水平流到这些剖分区域中， 在每一个孔隙容积中，水都将与矿物， 阳离子交换和这区域外表中的水质相平衡。 在这些区域中的水化学的演化代表了距含水层饱和区顶部点处的水的化学演化。初始条件初始条件Parkhurst 和他人1996提供了如下的数据：在含水层中的1

245、 升水中估计大约含有的钙离子，白云石和阳离子交换点的摩尔数。钙的重量摩尔围是从 0 到 2%，且白云石是从 0到 7%，当然，白云石的重量可能是更多。孔隙度声明为0.22。粘土阳离子交换的围是从20到 50 meq/100 g，与其 30%粘土的平均含量。在这些例子的计算中，钙离子假定为0.1 的重量百分点和白云石假定为3%的重量百分点， 这是通过假定岩石的密度为2.7kg/L 而得出的，69 / 102同时它对应于0.1mol/L的钙与其1.6mol/L的白云石。 阳离子交换点的数量假定为1.0eq/L。在溶液外表的砷估计是从核心样品Mosier 和他人, 1991中连续萃取而得出的数据。在

246、固体相中砷的浓度的变化围是10 to 到 20ppm，这与 1.3 到 2.6 mmol/L 的砷相对应。外表点的估计是与来自于沉淀物中可萃取的铁离子的数量相关，它的围是由1.6 到 4.4%Mosier 和他人, 1991。沉淀物中2%铁离子的含量是与 3.4mol/L 的铁相对应。然而，大多数的铁离子是在针铁矿和赤铁矿中， 这比水相氧化铁离子具有更少的外表粘合点。 在水相氧化铁离子中的铁局部的含量武断地假定为0.1。这样，总数为0.34mol 的铁离子是假定在水相中的氧化铁离子， 并且应用值为 0.2 表表示每摩尔铁离子的数量， 在计算中每升水中总共有 0.07mol 的点。应用 89 的

247、克分子式量来估计水相氧化铁离子量为 30g/L。这个指定的外表区域是指定为 600m /g。2表表3737例例1414输入数据的设置输入数据的设置TITLE Example 14.-Transport with equilibrium_phases, exchange, and surfacereactions*PLEASE NOTE: This problem requires database file wateq4f.dat!*SURFACE_SPECIES Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+ log_k -15. Hfo_wOH + Ca+2 = Hfo_w

248、OCa+ + H+ log_k -15.SOLUTION 1 Brine pH 5.713 pe 4.0 O2(g) -0.7 temp 25. units mol/kgw Ca .4655 Mg .1609 Na 5.402 Cl 6.642 charge C .00396 S .004725 As .05 umol/kgwENDUSE solution 1EQUILIBRIUM_PHASES 1 Dolomite 0.0 1.6 Calcite 0.0 0.1SAVE solution 1# prints initial condition to the selected-output f

249、ileSELECTED_OUTPUT70 / 102 -file ex14.sel -reset false -stepUSER_PUNCH -head m_Ca m_Mg m_Na umol_As pH10 PUNCH TOT(Ca), TOT(Mg), TOT(Na), TOT(As)*1e6, -LA(H+)ENDPRINT# skips print of initial exchange and initial surface to the selected-output file -selected_out falseEXCHANGE 1 -equil with solution 1

250、 X 1.0SURFACE 1 -equil solution 1# assumes 1/10 of iron is HFO Hfo_w 0.07 600. 30.ENDSOLUTION 0 20 x precipitation pH 4.6 pe 4.0 O2(g) -0.7 temp 25. units mmol/kgw Ca .191625 Mg .035797 Na .122668 Cl .133704 C .01096 S .235153 chargeEQUILIBRIUM_PHASES 0 Dolomite 0.0 1.6 Calcite 0.0 0.1 CO2(g) -1.5 1

251、0.SAVE solution 0ENDPRINT -selected_out trueADVECTION -cells 1 -shifts 200-print_frequency 20END71 / 102最初注入到含水层中的盐水是来自于Parkhurst和他人1996的成果且是在这个例子的输入数据设置中是作为溶液1表37给定的。纯相集合含有钙和白云石是以EQUILIBRIUM_PHASESEQUILIBRIUM_PHASES 1数据块所定义的。阳离子交换点的数目是以EXCHANGEEXCHANGE 1所定义的且外表点的数目是以SURFACESURFACE 1所定义的。初始的交换和初始外表组

252、分的定义是与盐水相平衡而达到的，这是在与钙和白云石相平衡之后达到的 注释交换物质和外表的平衡， 是在矿物平衡之前达到，通过缓冲注入吸附元素会产生不同的结果 。在盐水中砷的浓度是通过屡次试验得到的，且给定在溶液外表上的误差砷大约为2mmol，这与随后的萃取数据相一致。数据库wateq4f.dat它包含有砷元素和来自于Dzombak和Morel的外表络合物组分适用于除两种外表反响的热动力学的数据。在每一次运行之后，它决定了砷浓度更好结果的取得，在这种情况下，钙和锰络合反响都将排除。 这个数据块是用来降低平衡常数， 这个平衡常数会降低每丙种反响的大约有10个数量级。这将有效地排除了钙和锰的溶液外表的

253、络合反响。这个SURFACE_SPECIESSURFACE_SPECIES数据块是用来降低每两种反响的平衡常数大约有10个数量级。这也将有效地排除了钙和锰的外表络合反响。 可选性的， 这些反响能够从缺省的数据库中排除 。如果阳离子和阴离子没有为同一点而竞争， 这将得到调整；一般而言，在阳离子和阴离子之间的竞争吸附是不为人知的。重新灌注水重新灌注水进入到含水层饱和区域中的水假定是在渗流区PCO2 为10-1.5时的钙和白云石相平衡。这输入设置中的第四个模拟这个模拟紧跟在第三个ENDEND声明之后产生了这样的水的组分并且它是应用SAVESAVE表37作为溶液0来保存的。平流运移计算这个ADVECT

254、IONADVECTION数据块表37提供了用来平流灌注水到代表饱和区域中区域的必要信息。总共指定了200的转换点，这等于200的孔隙容积，这是因为在这个计算中仅有一种单个的区域。计算的结果是列于表14中。在初始5个孔隙体积中，高浓度的钠，钙和锰都会得以降低以致使钠是主要的阳离子与钙和锰的浓度是非常小的。这个pH值升高至超过9.0。且砷的浓度升高至接近于2 mol/kgw。后面45个孔体积中pH逐渐降低且砷的浓度也降低至可忽略的浓度。在大约100个孔隙体积中，钙和锰会变成主要的离子与其稳定的pH值，这出现在后面灌注水的pH值中。72 / 102图 14地下水化学演化运移模拟的结果，主要是由于钙、

255、镁、重碳酸盐水流入到最初含有卤水、方解石、白云石、离子交换剂和含有砷配合物外表的含水层中。平流运移反响产生了三种类型的水， 这与在含水层中所观测到的水的类型相似： 初始盐水， 一种重碳酸盐水， 和重碳酸钙锰水。 所计算的pH值是与所观测到的含水层中的水相一致。在主要含有钠离子的水域中， 所计算的pH值一般大于8.0且有时甚至高至9.2； 在重碳酸钙锰水域中， pH稍微大于7.0。 更为敏感的计算说明了pH的最大值是依赖于当前交换物质的数量。阳离子交换点数目的减少降低了最大的pH。 所模拟砷的浓度的值是与在含水层中所观测到的值相类似，这里的最大浓度为1到2mol/kgw。PH最大值的降低会使砷浓

256、度最大值降低。溶液外表络合反响稳定常数是直接取自于文献中； 占主要砷络合反响的logK的降低也倾向于降低砷的最大浓度。在结论中，模拟的结果，绝大局部是依赖于所测量的值和文献的热动力学数据，这提供了在在含水层中主要离子化学反响中的pH，和砷的浓度变化的满意解释。73 / 102？例例 1515一维运移：动力学生物降解，单元增长与吸附一维运移：动力学生物降解，单元增长与吸附平流弥散反响运移的测试问题是经由Tebes-Steven and Valocchi(1997, 1998)开展起来的。 尽管是基于相对简单的物种形成反响， 对于这个问题的解答证明了几种相互的化学反响进程，这根本上是与许多环境问题

257、相关： 以细菌为媒介的有机物的降介； 细菌区域的增长和其腐烂；和其水相物种的形成包括金属络合物。在这个例子中，所测试问题的解决是用PHREEQC来进展解决的，这产生的结果几乎与那些由Tebes-Steven和Valocchi (1997, 1998)所提出的结果相类似。然而，在使用PHREEQC中的平流弥散运移模拟时， 细心是很有必要的，这正如许多反响运移所模拟的一样，用来确保准确数字化溶液的取得。所测试的问题模拟了运移的进程， 这是当含有NTAnitrylotriacetate的水和其钴注入到这一栏中而发生的。这一问题包含有平流和这一栏中的弥散反响与其水相平衡反响和NTA降介的动力学反响，细

258、菌的生长和钴的吸附。运移参数运移参数这些维数和这一栏的水动力学属性列于表38中。表38例15柱子的水动力学和物理学性质性质柱子长度孔隙度体积密度值10.0m.41.56e6g/m3每升水中沉淀物 来自于孔隙度和体积密度 3.7533g/L孔中水流速度纵向弥散度水相模型水相模型Tebes-Steven和Valocchi (1997)定义了水相模型用来对包含有物质种类和其水相离子特性和其物质的 logKs 值进展测试；其活度系数假定为1.0。列于表 39 中的数据库文件的建立是基于它们的水相模型建立的。对PHREEQC 模拟而言，NRA 是作为一种新的“元素来定义的，在 SOLUTION_MAST

259、ER_SPECIESSOLUTION_MASTER_SPECIES 数据块中命名为“Nta。来自于“NTA中的这一点都将指定为“Nta，这是为了与 PHREEQC 中的概念相一致。在 SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES中的 Nta 的克分子式量是非物质性的，这是当其在 SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES 数据块中输入单位为摩尔，且简单地设为 1。这个 Nta 的水水相络合是定义在 SOLUTION_SPECIESSOLUTION_SPECIES 数据块中。注意水相中所有物质的活度

260、系数都定义了a参数的较大值1x107以gammagamma 标识符标识，这将使得其活度系数非常的接近于1.0。表 39例 15 的数据库74 / 1021.0m/hr.05mSOLUTION_MASTER_SPECIESC CO2 2.0 61.0173 12.0111Cl Cl- 0.0 Cl 35.453Co Co+2 0.0 58.93 58.93E e- 0.0 0.0 0.0H H+ -1. 1.008 1.008H(0) H2 0.0 1.008H(1) H+ -1. 1.008N NH4+ 0.0 14.0067 14.0067Na Na+ 0.0 Na 22.9898Nta N

261、ta-3 3.0 1. 1.O H2O 0.0 16.00 16.00O(-2) H2O 0.0 18.016O(0) O2 0.0 16.00SOLUTION_SPECIES2H2O = O2 + 4H+ + 4e- log_k -86.08; -gamma 1e7 0.02 H+ + 2 e- = H2 log_k -3.15; -gamma 1e7 0.0H+ = H+ log_k 0.0; -gamma 1e7 0.0e- = e- log_k 0.0; -gamma 1e7 0.0H2O = H2O log_k 0.0; -gamma 1e7 0.0CO2 = CO2 log_k 0

262、.0; -gamma 1e7 0.0Na+ = Na+ log_k 0.0; -gamma 1e7 0.0Cl- = Cl- log_k 0.0; -gamma 1e7 0.0Co+2 = Co+2 log_k 0.0; -gamma 1e7 0.0NH4+ = NH4+ log_k 0.0; -gamma 1e7 0.0Nta-3 = Nta-3 log_k 0.0; -gamma 1e7 0.0Nta-3 + 3H+ = H3Nta log_k 14.9; -gamma 1e7 0.0Nta-3 + 2H+ = H2Nta- log_k 13.3; -gamma 1e7 0.0Nta-3

263、+ H+ = HNta-2 log_k 10.3; -gamma 1e7 0.0Nta-3 + Co+2 = CoNta-75 / 102 log_k 11.7; -gamma 1e7 0.02 Nta-3 + Co+2 = CoNta2-4 log_k 14.5; -gamma 1e7 0.0Nta-3 + Co+2 + H2O = CoOHNta-2 + H+ log_k 0.5; -gamma 1e7 0.0Co+2 + H2O = CoOH+ + H+ log_k -9.7; -gamma 1e7 0.0Co+2 + 2H2O = Co(OH)2 + 2H+ log_k -22.9;

264、-gamma 1e7 0.0Co+2 + 3H2O = Co(OH)3- + 3H+ log_k -31.5; -gamma 1e7 0.0CO2 + H2O = HCO3- + H+ log_k -6.35; -gamma 1e7 0.0CO2 + H2O = CO3-2 + 2H+ log_k -16.68; -gamma 1e7 0.0NH4+ = NH3 + H+ log_k -9.3; -gamma 1e7 0.0H2O = OH- + H+ log_k -14.0; -gamma 1e7 0.0END初始和边界条件初始和边界条件背景值的浓度是列于表 40 中这一栏中。这一栏中最初不

265、包含有 Nta 或是钴，但含有1.36x10-4 g/L 生物数量。流运的边界条件适用于这一栏的入口处，并且在开始的 20 个小时会遇到在溶液中含有 Nta 和钴；这些数量同样也列于表40 中。在 20 小时之后，溶液将会介绍至入口中至达实验达到 75 小时之后。Na 和 Cl 也不出现在最初问题的定义中，但将被PHREEQC 加到电荷平衡吸附反响中看下面的“吸附反响 。表表4040例例1515的浓度数据的浓度数据组分组分H+总碳NH4+O2Nta3-Co2+NaCl单位面积生物数量ClNta(ads)水相水相水相水相水相水相水相水相不流动不流动类型类型脉冲浓度脉冲浓度10.0e-6mol/L

266、4.9e-7mol/L.03.125e-5mol/L5.23e-6mol/L5.23e-6mol/L1.0e-3mol/L1.0e-3mol/L-76 / 102背景浓度背景浓度10.0e-6mol/L4.9e-7mol/L.03.125e-5mol/L.0.01.0e-3mol/L1.0e-3mol/L1.36e-4g/L.0Co(ads)不流动-.0NtaNta 的动力学降解和剖分区域的增长的动力学降解和剖分区域的增长Nta假定在出现生物数量和氧的情况下以以下的反响来进展降解的：HNta21.62O21.272H2O 2.424H 0.576C5H7O2N 3.12H2CO30.424NH

267、4PHREEQC需要单独定义其动力学反响物，在进入或是从溶液中去除的每一种元素的反响物都是以摩尔数表示的。进一步而言，反响物应当达到电荷平衡 没有净电荷应当进入或是离开溶液 。这个Nta转换1mol的HNta2- (C6H7O6N)到0.576 mol C5H7O2N，这里的后者是不活泼的， 因此它的浓度可以被忽略。 在这两种反响物中所包含的元素数量的不同提供了反响中元素C，H，O和N的化学计量。这个化学计量等于等式右边元素的总和，减去除了C5H7O2N，的等式左边的元素的总和。在水相中每摩尔的HNta2-反响中的元素的相关变化是列于表41中正的系数标明了水相浓度的增高，而负的系数说明了水相浓

268、度的降低。表表4141NtaNta氧化复原反响反响的计量氧化复原反响反响的计量组分组分NtaCHON下面的这个乘法Monod速率表达式是用来描述Nta的速率降解：-1.03.121.9684.848.424系数系数RHNTA2 qmXm(cHNTA2kscHNTA2)(cO2KacO2)这里RHNT A2-是HNta2-降解的速率(mol/L/hr),qm是使用酶mol/g cells/hr时所指定的最大速率，Xm是单位面积的生物数量g cells/L，是酶Ntamol/L的半饱和常数，Ka是电子承受者O2 (mol/L)的半饱和常数，与ci说明的是浓度mol/L。所产生生物数量是依赖于使用酶

269、的速率和生物数量的腐烂速率：Rcells YRHNTA2bXm这里的Rcells是剖分区域增长的表达式 g cells/L/hr ， Y是微生物产量系数 g cells/molNta，和b是生物数量的腐烂速率hr-1。这些等式参数的值是列在表42中。表表4242例例1515中所用到动力学速率参数中所用到动力学速率参数参数参数描述描述77 / 102参数值参数值KsKaQm给体半饱和常数承受体半饱和常数使用培养基的最大指定速率7.64e-7mol/L6.25e-6mol/L1.418e-3mol Nta/gcells/hrYb吸附反响吸附反响微生物生产系数一阶微生物腐烂系数65.14 g cel

270、ls/mol Nta0.00208 hr-1Tebes-Steven和Valocchi (1997)以下面的等式定义了Co2+ 和CoNta-的动力学吸附反响：Ri km(cisi)kd这里i或者是Co2+或CoNta- (mol/L)，si是吸附浓度mol/g沉淀物，km是质量转移系数hr-1，与Kd是分布系数L/g。这个系数的值是列于表43中。表表4343Co2+Co2+和和CoNta-CoNta-的吸附系数的吸附系数种类种类Co2+CoNta-1 hr-11 hr-1kmkmKdKd5.073-3 L/g5.33e-4 L/gKd值的定义是给定了Co2+和CoNta-各自的滞流系数为20

271、和3。由于吸附反响定义的是动力学的，这些反响物的初始摩尔数和反响速率是以KINETICSKINETICS和RATESRATES数据块所定义的；没有溶液外表的物质SURFACESURFACE,SURFACE_MASTER_SPECIESSURFACE_MASTER_SPECIES, or SURFACE_SPECIESSURFACE_SPECIES是明确需要的。另外，所有的动力学反响物都是静止的，因此没有吸附的种类需要转移。当使用PHREEQC来模拟时，动力学反响物一定得达到电荷平衡。对于Co2+和CoNta-的吸附而言，1 mmol的NaCl被加到所定义的溶液中，以达到吸附进所需要的配对离子。

272、动力学吸附反响然后被定义移除或是进入依赖于所转移摩尔数到CoCl2 和NaCoNta溶液中，这是达到电荷平衡。为了从每千克的沉淀物si的吸附转移到到每升水中所吸附的摩尔数，乘以在每升水中的每克沉淀物是非常有必要的，3.75e3 g/L。输入数据的设置输入数据的设置表44说明了来自前面问题所定义的输入数据的设置。 尽管速率是以单位mol/L/hr来给定的， 在PHREEQC中速率通常是mol/s， 并且所有的速率在输入数据设置中都是以秒来定义速率的。假定在每一剖分区域中水的体积为1升，这在当前的问题中是合理的，其原因在于在每一种溶液中水的质量都接近于1kg，并且溶液都是相对冲淡的。如果在溶液中水

273、的质量别离为1kg，这可能打破常容积的假定。这个10米的栏可以别离1米的长度为10个剖分区域。开始两个灌注区域都定义为可灌注78 / 102的溶液且初始溶液编号为1到10。数据块定义了四种动力学反响表达式：HNta-2, 生物数量, Co_吸附, 和CoNta_吸附。这个速率表达式是以-start-start开始的， 以数字化定义为Basic语言的声明， 且是以-end-end来中止。每一表达式最后的声明是跟在变量之后的SAVE。 这个变量是经过一定时间间隔的反响的数字摩尔数并且它的计算是来自于即时速率表达式mol/s乘以时间间隔的长度 s，这是由变量“TIME来给定的。在第一和第二速率表达式

274、中的Line 30和20的速率表达式与其在第三和第四速率表达式中的第10行调整参数小时单位为秒单位。函数“MOL返回了物种的浓度 mol/kgw 与函数 “M返回了反响物的摩尔数， 这个速率表达式是这时被定义的。 “KIN返回的是动力反响物的摩尔数， 这可定义了一个区域中的任何反响物。 函数 “PUT和 “GET是用来保存和排除这样一个术语，这是Hnta-2和生物数量速率表达式是很普遍的。KINETICSKINETICS数据块定义了适用于每一个剖分区域的速率表达式；在这个例子中区域1到区域10是同时定义的。 对适用于每一个区域的每一个速率表达式， 反响物的分子式 -formula-formul

275、a和反响物的初始摩尔数在开始时就是被定义-m-m，如果需要的。定义一个误差是可行的-tol-tol， 以摩尔表示，这是对于速率表达式的数字综合和准确性而言的。 注意HNta-2速率表达式产生了一个负的速率， 因此，从溶液中移除元素分子式的系数是正数， 而元素加到溶液中时，系数为负数。一般而言，速率和系数的积是正数，表示元素进入到溶液中，但其积如果是负数，那么表示，元素离开溶液。生物质量反响加了“ H 0.0，或是没有其摩尔数，这指定了生物质量增长的动力学反响， 没有或是从溶液中移除。 碳的同化作用和与生物生长有关的营养物质在这个模拟中都将忽略。SELECTED_OUTPUTSELECTED_O

276、UTPUT数据块使水中物种Nta-3, CoNta-, HNta-2和Co+2写到数据块文件ex15.sel中。对文件中的每一行而言，这个数据块所填加的时间以时间表示，与其吸附浓度转换为mol/g沉淀物，和生物数量。TRANSPORTTRANSPORT数据块定义了实验开始的前20个小时， 这是对Nta和钴加到这一栏输入期间而言的。如果这一栏定义有一米的长度-length-length有10个区-cells-cells。这个平流弥散运移模拟的持续时间为3600秒-time_step-time_step的20步-shifts-shifts。流动的方向是向前的-flow_direction-flow

277、_direction。每一栏的最后都定义有一个流动的边界条件-boundary_condition-boundary_condition。这个弥散度为0.05米弥散度和弥散系数设置为0-diffusion_coef-diffusion_coef 。 在每一次转换之后 -punch_frequency-punch_frequency 仅有10个区 -punch_cells-punch_cells的数据写到了所选取的输出文件中。在每五个转换-print_frequency-print_frequency之后仅有10个区的数据写到输出文件中-print_cells-print_cells。在第一个平

278、流-弥散运移模拟之后，一种新的填入溶液会被定义SOLUTIONSOLUTION 0，这没有包含有Nta或钴。对相关的初始溶液计算而言，打印到所选取的输出文件将会去除，然后重新声明-selected_out-selected_out为假和-selected_outselected_out为真，在PRINTPRINT数据块中。最后，第二个TRANSPORTTRANSPORT数据块定义了实验的最后55个小时，这时，在灌注溶液中不出现Nta和钴。所有的参数与前面TRANSPORTTRANSPORT数据块中的参数是一致的，只是平流的步数79 / 102-shifts-shifts会增加到55。表表444

279、4例例1515输入数据的设置输入数据的设置TITLE Example 15.-1D Transport: Kinetic Biodegradation, Cell Growth, and Sorption*PLEASE NOTE: This problem requires database file ex15.dat!*SOLUTION 0 Pulse solution with NTA and cobalt units umol/L pH 6 C .49 O(0) 62.5 Nta 5.23 Co 5.23 Na 1000 Cl 1000ENDSOLUTION 1-10 Backgroun

280、d solution initially filling column # 1-20 units umol/L pH 6 C .49 O(0) 62.5 Na 1000 Cl 1000ENDRATES Rate expressions for the four kinetic reactions# HNTA-2 -start10 Ks = 7.64e-720 Ka = 6.25e-630 qm = 1.407e-3/360040 f1 = MOL(HNta-2)/(Ks + MOL(HNta-2)50 f2 = MOL(O2)/(Ka + MOL(O2)60 rate = -qm * KIN(

281、Biomass) * f1 * f270 moles = rate * TIME80 PUT(rate, 1) # save the rate for use in Biomass rate calculation90 SAVE moles -end# Biomass -start10 Y = 65.1480 / 10220 b = 0.00208/360030 rate = GET(1) # uses rate calculated in HTNA-2 rate calculation40 rate = -Y*rate -b*M50 moles = -rate * TIME60 if (M

282、+ moles) 0 then moles = -M70 SAVE moles -end# Co_sorption -start10 km = 1/360020 kd = 5.07e-330 solids = 3.75e340 rate = -km*(MOL(Co+2) - (M/solids)/kd)50 moles = rate * TIME60 if (M - moles) 0 then moles = M70 SAVE moles -end# CoNta_sorption -start10 km = 1/360020 kd = 5.33e-430 solids = 3.75e340 r

283、ate = -km*(MOL(CoNta-) - (M/solids)/kd)50 moles = rate * TIME60 if (M - moles) 0 then moles = M70 SAVE moles -endKINETICS 1-10 Four kinetic reactions for all cells # 1-20 HNTA-2 -formula C -3.12 H -1.968 O -4.848 N -0.424 Nta 1. Biomass -formula H 0.0 -m 1.36e-4 Co_sorption -formula CoCl2 -m 0.0 -to

284、l 1e-11 CoNta_sorption -formula NaCoNta -m 0.0 -tol 1e-11SELECTED_OUTPUT81 / 102 -file ex15.sel -mol Nta-3 CoNta- HNta-2 Co+2USER_PUNCH -heading hours Co_sorb CoNta_sorb Biomass -start10 punch TOTAL_TIME/3600 + 1800/3600 # TOTAL_TIME/3600 + 900/360020 punch KIN(Co_sorption)/3.75e330 punch KIN(CoNta_

285、sorption)/3.75e340 punch KIN(Biomass) -endTRANSPORT First 20 hours have NTA and cobalt in infilling solution -cells 10 # 20 -length 1 # 0.5 -shifts 20 # 40 -time_step 3600 # 1800 -flow_direction forward -boundary_condition flux flux -dispersivity .05 -correct_disp true -diffusion_coef 0.0e-9 -punch_

286、cells 10 # 20 -punch_frequency 1 # 2 -print_cells 10 # 20 -print_frequency 5 # 10ENDPRINT -selected_out falseSOLUTION 0 New infilling solution, same as background solution units umol/L pH 6 C .49 O(0) 62.5 Na 1000 Cl 1000ENDPRINT -selected_out trueTRANSPORT Last 55 hours with background infilling so

287、lution -shifts 55 # 110END网格收敛网格收敛利用平流-弥散-反响运移模拟来检查结果的数字准确性是必要的。 一般而言， 对这些复杂的模拟而言没有分析溶液， 因此数字准确性的测试重新定义了网格和时间步， 重新运行模82 / 102拟与比照结果， 如果在两种不同的网格上给出了两种相似的结果， 那么数字化的一些确定性是相对小的。 如果基于两种不同网格上的模拟给出了重要的不同结果， 这个网格必须重新定义与这个过程必须重复。 不幸的时， 使网格尺寸扩大两倍， 至少有四倍数目的溶液必须制取，这是因为区域的数目进展翻倍和时间步减半。 如果区域的大小接近弥散度的大小， 这可能需要更多的溶

288、液计算，这是因为弥散计算中混合的步数也将会增高。为了测试这个例子中的网格收敛， 对总共这20个区域而言， 柱子中区域的数目将是翻倍的。1-10围的所有定义组分的数据块都将变为1-20。另外，平流-弥散运移的参数将调整至与这些新剖分区域的数目相一致。 区域和转换的数目都将翻倍； 区域长度和时间步都将减半。为了打印与这10个区域模拟中一样位置的信息柱子的末端，-punch_cells-punch_cells和-print_cells-print_cells设为区域20。 为了打印这10个区域模拟中一样时间的信息， -punch_frequency-punch_frequency设置为每两个转换，-

289、print_frequency-print_frequency设置为每10个转换，以与从区域中点到柱子完毕时间步的运行在第10行以USER_PUNCHUSER_PUNCH进展减半。 这20个区域所有的变化同样也标注在表44中以行完毕的注释进展。结果结果T柱子中水相和不流动要素的分布是在75个小时完毕之后，列于10和20区域模拟的图15和16中，在这个实验中，含有Nta和钴的两种孔隙水在开始的20小时中参加到柱子中，后面55小时参加的是5.5体积的背景水。在10个小时时，HNta2-开始在柱子出口处出现，此时伴着pH值的升高图15。如果Nta和钴是保守的，并且弥散是可无视的，那么这个图说明的是在

290、10小时升高和30小时降低的矩形脉冲。然而，Nta和钴相对于吸附反响的保守运动而言是停止的。 Nta和钴浓度的顶峰发生在30和40小时之间的CoNta混合物中。 Co2+浓度的顶峰对它的吸附反响而言更是延迟的，且至到实验完毕时才会出现的。图 15. 分别在 10 和 20 区域 Nta 和钴运移模拟柱子出口处的水中物质浓度和 pH 值。在图16中， 在柱子最后区域的浓度的固相浓度是概括时间来进展划分的。 吸附的CoNta-浓度顶峰的出现在30和40小时之间出现的，稍稍落后于CoNta-配合物溶解浓度顶峰的出现。83 / 102开始， 在柱子中没有Nta出现， 并且单位面积的生物数量在开始的10

291、个小时之有轻微的降低，这是由于生物数量第一次序的腐烂速率。 随着，Nta移进这个区域，底层Nta单位面积生物数量的增长变为可能。在柱子中出现Nta峰之后，单位面积生物数量的浓度开始变得平稳，然后，由于腐烂而开始降低。钴吸附的Kd比CoNta-吸附的Kd具有更大一些的延迟系数，Co2+吸附浓度出现在实验完毕的停止增长上。图 16. 分别在 10 和 20 区域 Nta 和钴运移模拟柱子出口处吸附离子的浓度和单位面积生物数量。10-区域和20-区域模拟给出了一样的结果，这说明了，平流-弥散运移模拟的数字上的错误是相对很小的，这个结果与Tebes-Steven和Valocchi (1997, 199

292、8)所给出的结果相类似。然而，Tebes-Steven和Valocchi (1997)包括了他们测试问题的其它局部，也就是从1到1000hr-1的吸附反响速率常数的升高。 所增加的速率常数产生了一组偏微分等式， 速率限制的进程发生在不同的时间衡量上。具有快速吸附反响的问题证明了PHREEQC显性算法的难以处理性，但当动力学吸附的计算是作而有效组成的平衡过程时，能够被成功的加以解决。然而，即使有平衡吸附，网格收敛的计算具有更大的强度；它需要使用100区域或是更多来达到满意的溶液。相关CPU次数的估计，10-区域的模拟使用270秒，20-区域的模拟需要732秒，这是在Pentium I, 133

293、MHz的计算机上而言的。一个200区域的模拟大约需要600倍的10-区域模拟的时间。例例 1616SierraSierra 泉水的逆向模拟泉水的逆向模拟该例对Sierra Nevada地区泉水组分的化学演化进展了重复的反向模拟计算，这些计算在Garrels和Mackenzie (1967)的文章中已经描述过。在反向模拟程序ETPATH (Plummer等，1991和1994)的操作手册中描述过同样的例子。这个例子使用了两种泉水的组分，一种是暂时泉水，它经历的化学演化较少； 而另一种那么是长年存在的泉水， 它在下层土中可能已存在很长的时间了。 这两种泉水组分的差异是假定与它们与相接触的水、 矿物

294、以与气体发生了84 / 102反响。在这个模拟中反向模拟的目的是， 当参与反响的量适宜时， 找到矿物和气体可定量地解释这两种溶液组分之间的不同。表表4646例例1616泉水的分析数据泉水的分析数据单位为毫摩尔/升，来自于Garrels和Mackenzie (1967)泉泉pHpHSiOSiO2 2CaCa2+2+MgMg2+2+NaNa+ +K K+ +HCLHCL3 3- -SOSO4 42-2-ClCl- -暂时泉 6.20.2730.0780.0290.1340.0280.3280.0100.014长年泉6.8.410.260.071.259.040.895.025.030NETPATH

295、 (Plummer等，1991, 1994)和PHREEQC都具有能够进展反向模拟运算的能力。NETPATH相对PHREEQC而言有两个明显的优点： 1 NETPATH提供了完整的处理同位素的方法，包括同位素的摩尔平衡，同位素分馏，以与 C 示踪，然而PHREEQC仅具有同位素摩尔平衡的能力，和2NETPATH提供了数据输入和模拟进展的完全交互环境，然而PHREEQC版本2主要是批反响计算的程序。另一方面，在PHREEQC版本1中也是可应用完全的图形用户界面，不过它缺少同位素摩尔平衡能力Charlton等，1997， PHREEQC版本2也可应用完全的图形用户界面PHREEQC forWind

296、ows, V.E.A. Post, writtencommun.,1999,.geo.vu.nl/users/posv/phreeqc.html。PHREE相对于NETPATH而言，进展反向模拟的主要优点是， 反向模拟计算中包括分析数据不确定度的计算。 这种能力使反向模拟具有更强的生命力， 也就是说， 在输入数据中的微小的变化不会导致模拟计算的摩尔转换产生大的变化。PHREEQC的另一种优点是所设定的任何元素都可以包含在反向模拟计算中，然而，NETPATH那么局限于所选取的一组元素，这相对要全面一些。表表47.47. 由由Garrels and Mackenzie (1967)Garrels

297、and Mackenzie (1967)给出的反响物组分和摩尔转换给出的反响物组分和摩尔转换摩尔转换单位为毫摩尔/千克水，正数表示溶解，负数表示沉淀-14反响物反响物“岩盐“硬石膏高岭石Ca-蒙脱石CO2气体方解石硅黑云母斜长石组分组分NaClCaSO42H2OAl2Si2O5(OH)4Ca0.17Al2.33Si3.67O10(OH)2CO2CaCO3SiO2KMg3AlSi3O10(OH)2Na0.62Ca0.38Al1.38Si2.62O8摩尔转移摩尔转移0.016-0.015-0.033-0.810.4270.1150.00.0140.175这两种泉水的分析数据见表46。由Garrel

298、s和Mackenzie (1967)给出的假定参与反响的矿物和气体的化学组分以与它们的摩尔转换见表47。 反响相的特性和组分的选择是反向模拟最难的局部。一般而言，这局部是根据流动系统和沿流动路径的矿物学的知识得出的； 含水层物质成分的微观和化学分析、 水和矿物的同位素组分为这局部的反响物的选择提供了附加的信息。 没有必要准确地知道是哪一种矿物参加了反响， 但却有必要全面知道潜在的反响物。这个例子的输入数据组见表48。SOLUTION_SPREADSOLUTION_SPREAD 数据块用来定义这两种泉水。85 / 102INVERSE_MODELINGINVERSE_MODELING 数据块用来

299、定义所有反向模拟计算的特性，包括所用到的溶液和相，摩尔平衡方程，不确定围，以与是否所有的或仅是“最小数量 的模拟将会被打印，是否与不确定限一致的围的摩尔转换都将会被计算。一系列的标识符在连字符之后的子关键字指定了反向模拟的特性。表表48.48. 例例1616的数据输入的数据输入TITLE Example 16.-Inverse modeling of Sierra springsSOLUTION_SPREAD -units mmol/LNumber pH Si Ca MgNa K Alkalinity S(6)Cl1 6.2 0.273 0.0780.0290.1340.0280.328 0.

300、010.0142 6.8 0.41 0.260.0710.2590.04 0.8950.0250.03INVERSE_MODELING 1 -solutions 1 2 -uncertainty 0.025-balancesCa 0.05 0.025 -phases Halite Gypsum Kaolinite precip Ca-montmorillonite precip CO2(g) Calcite Chalcedony precip Biotite dissolve Plagioclase dissolve-rangePHASESBiotite KMg3AlSi3O10(OH)2 +

301、 6H+ + 4H2O = K+ + 3Mg+2 + Al(OH)4- + 3H4SiO4 log_k 0.0 # No log_k, Inverse modeling onlyPlagioclase Na0.62Ca0.38Al1.38Si2.62O8 + 5.52 H+ + 2.48H2O = 0.62Na+ + 0.38Ca+2 + 1.38Al+3 + 2.62H4SiO4 log_k 0.0 # No log_k, inverse modeling onlyEND-solutions-solutions 标识符选择应用的溶液的序号，两种或是更多的溶液序号必须在其标识符之后86 / 1

302、02列出。如果仅仅是给出了两种溶液的序号， 第二种溶液假定是从第一种溶液演化来的。 如果给出了更多的溶液序号， 最后列出的溶液假定是由前面溶液的混合物演化来的。 在反向模拟中所用到溶液的定义与在PHREEQC 中模拟所用到任何溶液的定义是一样的。 通常分析数据是在 SOLUTIONSOLUTION 或 SOLUTION_SPREADSOLUTION_SPREAD 数据块中输入的，但在当前或是前面模拟中批反响计算中所定义的溶液也可使用，条件是它们须用SAVESAVE 关键字进展保存。-uncertainty-uncertainty 标识符设定了每一分析数据的缺省不确定限。在这个例子中，相对不确定

303、限为 0.025(2.5%)，它适用于除 pH 之外的所有分析数据。在缺省状态下，pH 值的不确定限为 0.05 个单位。pH 的不确限也可以用-balances-balances 标识设置为绝对值 标准单位 。任何溶液任何数据的不确定限， 能够用-balances-balances 标识符明确地设置为小数值或是绝对值 摩尔数；碱度当量 。在缺省状态下，每个反向模拟包括包含于-phases-phases 中除了氢和氧任何相的每一种元素的摩尔平衡方程。如果需要不包括在相中的元素的摩尔平衡方程， 也就是说，元素没有源或汇守恒混合 ，-balances-balances 标识符能够用来包含那些在反向

304、模拟方程中的公式中的元素见例 17 。另外，-balances-balances 标识符能够用来指定每一种溶液中每一种元素的不确定限。在这个例子中为了证明这个目的，钙离子的不确定限在溶液 1 中设置为 0.055%和在溶液 2 中设置为 0.0252.5%在反向模拟计算中用到的相用-phases-phases 标识符进展定义。另外，这个标识符可以用来约束某相仅被溶解，或是仅为沉淀。在这个例子中，高岭石，钙-蒙脱石和玉髓(SiO2)都要求仅为沉淀。 这意味着在包含有高岭石相的任何模拟中高岭石都将作为沉淀生成 摩尔转移为负 ；同样，也适用于钙-蒙脱石和玉髓。如果在反向模拟中出现黑云母和斜长石， 那

305、么要求它们仅为溶解摩尔转移为正 。在反向模拟中用到所有的相必须在数据块文件中或是在输入文件中用PHASESPHASES 或EXCHANGE_SPECIESEXCHANGE_SPECIES 数据块来定义。这样在所有默认数据库文件中所定义的所有相，phreeqc.dat或wateq4f.dat在反向模拟中都可用。 黑云母和斜长石不在缺省的数据库文件phreeqc.dat中，因此，在输入文件设置中的以 PHASESPHASES 数据块显性定义。为简单起见，这些相的 logKs 设置为 0，由于仅是应用到了矿物的化学计量，这样设置不会影响反向模拟。然而，这些相饱和指数的计算将会是不合实际的。 在反向模

306、拟中用到的相均须具有电荷平衡反响。 这是因为每一种溶液都必须处于电荷平衡状态。 在用最正确化方法减少目标函数的过程中见“反向模拟方程和数值方法 ，通过在不确定限围调整元素的浓度，保证每个模型中的每一种溶液均达到电荷平衡。 如果一种溶液在给定的不确定限围不能调整到电荷平衡， 那么在输出文件中将会标记出该溶液和找不到模型。 由于在模拟过程中的所有溶液都处于电荷平衡，相也必须达到电荷平衡，或在模型中没有该相。注意，斜长石的反响表48有两行，由于在第一行的最后以反斜杠“完毕，程序将这两行解释为单逻辑行。-range-range 标识符说明，除了发现所有的反向模拟，在不确定限的围之，发现的每一种模型都需

307、要进展额外的计算以决定摩尔转移值的围。87 / 102以下方程也包含于每一种反向模拟之中： 系统中每种元素的摩尔数以与每种元素化合价如在-phases-phases 相中所定义的元素和在-balances-balances 中所列出的元素的摩尔平衡方程，每一种溶液的电荷平衡，体系的碱度平衡，体系的电子平衡和体系的水平衡。 在这些方程中的未知数包括相转移的摩尔数， 氧化复原反响转移的摩尔数， 和每一种溶液中每一种元素的不确定度不包括氢和氧 。不确定性包括每一种溶液的碱度和pH 值。优化方法解决了满足所有方程的未知数，并满足所有的不确定限， 同时最小化了目标函数， 这是不确定未知数的加权和见“反向

308、模拟方程和数值方法 。该例中发现的两个反向模拟的结果见表 49。模型中每一种溶液的每种成分用三列描述。所有的列除了 pH如果有同位素，也包括在 ，单位都是 mol/kgw。第一栏包括溶液的初始分析数据输入 。第二栏包括模拟中所计算的分析数据的修正(Delta)。这些修正必须在指定的不确定围。第三栏包含溶液修订后的分析数据，它等于原始数据加上修正值(Input+Delta)。在溶液列表之后，将会打印出反向模拟中每一种溶液的相比照例溶液分数 。在模型中仅有两种溶液，通常情况下每种溶液的分数都为1.0。在反向模拟中如果超过两种溶液，除了最后一种溶液外，其它溶液所占的百分数的和通常等于1.0。这个分数

309、实际上是来自于水的摩尔平衡，因此，如果消耗了含水矿物或是产生了一定数量的水， 或是模拟蒸发见例17 ，分数的和那么可能不是 1.0。在这个例子中，所有的分数都等于1.0；来自于石膏溶解水的量很小， 不影响混合分数的四个有效数字。 这个数据块的第二栏和第三栏给出的溶液分数的最小值和最大值，这是在指定的不确定限围取得的。只有使用了-range-range 标识符，这两栏才为非 0。在列表中的下一个数据块包括三列，分别描述了相的摩尔转移相摩尔转移 。第一栏包含与打印到溶液列表中修正的浓度相一致的反向模拟。 在这个例子中， 修正的溶液 1 加上在第一栏中的摩尔转移恰好等于修正后的溶液 2。摩尔转移为正

310、表示溶解；为负表示沉淀。注意：在多组反响计算中相集合的摩尔转移是相对于相， 而不是相对于溶液的第二和第三列的摩尔转移是在指定不确定限的围每种相摩尔转移的最小值和最大值。只有使用了-range-range 标识符，这两列为非0。一般而言，这个最小值和最大值不是独立的，也就是说，取得一种相最大的摩尔转移将会影响到摩尔平衡模拟中其它相的摩尔转移。在这个反向模拟中没有计算氧化复原摩尔转移。 如果计算了某个氧化复原摩尔转移， 每种元素不同化合价之间的摩尔转移将会被打印在标题“氧化复原摩尔转移之下。下一个数据块打印了在某种程度上与这个模拟修正的分析数据有关的结果； 如果没有进展任何修正，打印的这三个数将都

311、为0。首先，打印偏差的和，这是反向模拟不确定围的总和( uqmm,qm,qqmqm,qum,q)，其次，是反向模拟不确定围加权并修正的每种元素浓度和同位素组分的总和 uqmm,qm,q，delta/uncertainty不确定围的和。最后，88 / 102将会打印每种溶液中任何元素或是同位素组分 元素浓度的最大局部 的最大分数修正值 元素浓度的最大修正误差。所有的这三个值应用于打印在左栏中的模型。对于给定的摩尔平衡模型，如果没有发现包含适当的溶液和相子集的更为简单的模型，将会为给定的模型打印“模拟包含最少的相。所有模拟都打印之后是简短的计算总结， 列出了发现的模型数量、 最小模型的数量包含最少

312、相的模拟 ，所验证的不可实行模拟的数目， 不等式计算调用的数次cl1计算的时间通常与其调动的数目成比例 。例子的结果说明，使用Garrels和Mackenzie (1967)提出的相可发现两种反向模拟。主要反响是方解石和斜长石的溶解，消耗了二氧化碳；在第一次模拟中生成高岭石和钙-蒙脱石沉淀，在第二次模拟中生成高岭石和玉髓沉淀。 在模拟中需要溶解少量的岩盐， 石膏和黑云母。在除了二氧化碳之外，Garrels和Mackenzie (1967)得出的结果与其它所有相在第一个PHREEQC模型计算出的摩尔转移围。第一次模拟的二氧化碳转移摩尔数不同于Garrels和Mackenzie (1967)所得出

313、的值，这是因为他们没有说明泉水中二氧化碳的溶解。Garrels和Mackenzie (1967)也无视了钾离子摩尔平衡的微小差异。PHREEQC通过调整两种溶液中元素的浓度防止了钾离子的摩尔不平衡问题。PHREEQC计算说明，能够在不超过指定不确定限2.5%的围调整浓度得到两种反向模拟。如果不用PHREEQC进展计算，不考虑不确定限的数量级，就不清楚被Garrels和Mackenzie 无视的钾的差异是否重要。PHREEQC的结果与NETPATH的结果是相一致，但其条件是NETPATH也不考虑钾摩尔转移的差异。表表4949例例1616的选择输出结果的选择输出结果Solution 1: Inpu

314、t Delta Input+Delta pH 6.200e+000 + 1.246e-002 = 6.212e+000 Al 0.000e+000 + 0.000e+000 = 0.000e+000 Alkalinity 3.280e-004 + 5.500e-006 = 3.335e-004 C(-4) 0.000e+000 + 0.000e+000 = 0.000e+000 C(4) 7.825e-004 + 0.000e+000 = 7.825e-004 Ca 7.800e-005 + -3.900e-006 = 7.410e-005 Cl 1.400e-005 + 0.000e+000

315、 = 1.400e-005 H(0) 0.000e+000 + 0.000e+000 = 0.000e+000 K 2.800e-005 + -7.000e-007 = 2.730e-005 Mg 2.900e-005 + 0.000e+000 = 2.900e-005 Na 1.340e-004 + 0.000e+000 = 1.340e-004 O(0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+000 + 0.000e+000 = 0.000e+000 S(6) 1.000e-005 + 0.000e+000 = 1.000e-0

316、05 Si 2.730e-004 + 0.000e+000 = 2.730e-004Solution 2:89 / 102 Input Delta Input+Delta pH 6.800e+000 + -3.407e-003 = 6.797e+000 Al 0.000e+000 + 0.000e+000 = 0.000e+000 Alkalinity 8.951e-004 + -1.796e-006 = 8.933e-004 C(-4) 0.000e+000 + 0.000e+000 = 0.000e+000 C(4) 1.199e-003 + 0.000e+000 = 1.199e-003

317、 Ca 2.600e-004 + 6.501e-006 = 2.665e-004 Cl 3.000e-005 + 0.000e+000 = 3.000e-005 H(0) 0.000e+000 + 0.000e+000 = 0.000e+000 K 4.000e-005 + 1.000e-006 = 4.100e-005 Mg 7.101e-005 + -8.979e-007 = 7.011e-005 Na 2.590e-004 + 0.000e+000 = 2.590e-004 O(0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+00

318、0 + 0.000e+000 = 0.000e+000 S(6) 2.500e-005 + 0.000e+000 = 2.500e-005 Si 4.100e-004 + 0.000e+000 = 4.100e-004Solution fractions: Minimum Maximum Solution 1 1.000e+000 1.000e+000 1.000e+000 Solution 2 1.000e+000 1.000e+000 1.000e+000Phase mole transfers: Minimum Maximum Halite 1.600e-005 1.490e-005 1

319、.710e-005 NaCl Gypsum 1.500e-005 1.413e-005 1.588e-005 CaSO4:2H2O Kaolinite -3.392e-005 -5.587e-005 -1.224e-005 Al2Si2O5(OH)4Ca-Montmorillon-8.090e-005-1.100e-004-5.154e-005Ca0.165Al2.33Si3.67O10(OH)2 CO2(g) 2.928e-004 2.363e-004 3.563e-004 CO2 Calcite 1.240e-004 1.007e-004 1.309e-004 CaCO3 Biotite

320、1.370e-005 1.317e-005 1.370e-005 KMg3AlSi3O10(OH)2 Plagioclase1.758e-0041.582e-0041.935e-004Na0.62Ca0.38Al1.38Si2.62O8Redox mole transfers:Sum of residuals (epsilons in documentation): 5.574e+000Sum of delta/uncertainty limit: 5.574e+000Maximum fractional error in element concentration: 5.000e-002Mo

321、del contains minimum number of phases.=Solution 1:90 / 102 Input Delta Input+Delta pH 6.200e+000 + 1.246e-002 = 6.212e+000 Al 0.000e+000 + 0.000e+000 = 0.000e+000 Alkalinity 3.280e-004 + 5.500e-006 = 3.335e-004 C(-4) 0.000e+000 + 0.000e+000 = 0.000e+000 C(4) 7.825e-004 + 0.000e+000 = 7.825e-004 Ca 7

322、.800e-005 + -3.900e-006 = 7.410e-005 Cl 1.400e-005 + 0.000e+000 = 1.400e-005 H(0) 0.000e+000 + 0.000e+000 = 0.000e+000 K 2.800e-005 + -7.000e-007 = 2.730e-005 Mg 2.900e-005 + 0.000e+000 = 2.900e-005 Na 1.340e-004 + 0.000e+000 = 1.340e-004 O(0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+000 +

323、0.000e+000 = 0.000e+000 S(6) 1.000e-005 + 0.000e+000 = 1.000e-005 Si 2.730e-004 + 0.000e+000 = 2.730e-004Solution 2: Input Delta Input+Delta pH 6.800e+000 + -3.407e-003 = 6.797e+000 Al 0.000e+000 + 0.000e+000 = 0.000e+000 Alkalinity 8.951e-004 + -1.796e-006 = 8.933e-004 C(-4) 0.000e+000 + 0.000e+000

324、 = 0.000e+000 C(4) 1.199e-003 + 0.000e+000 = 1.199e-003 Ca 2.600e-004 + 6.501e-006 = 2.665e-004 Cl 3.000e-005 + 0.000e+000 = 3.000e-005 H(0) 0.000e+000 + 0.000e+000 = 0.000e+000 K 4.000e-005 + 1.000e-006 = 4.100e-005 Mg 7.101e-005 + -8.980e-007 = 7.011e-005 Na 2.590e-004 + 0.000e+000 = 2.590e-004 O(

325、0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+000 + 0.000e+000 = 0.000e+000 S(6) 2.500e-005 + 0.000e+000 = 2.500e-005 Si 4.100e-004 + 0.000e+000 = 4.100e-004Solution fractions: Minimum Maximum Solution 1 1.000e+000 1.000e+000 1.000e+000 Solution 2 1.000e+000 1.000e+000 1.000e+000Phase mole tr

326、ansfers: Minimum Maximum Halite 1.600e-005 1.490e-005 1.710e-005 NaCl Gypsum 1.500e-005 1.413e-005 1.588e-005 CaSO4:2H2O Kaolinite-1.282e-004-1.403e-004-1.159e-00491 / 102Al2Si2O5(OH)4 CO2(g) 3.061e-004 2.490e-004 3.703e-004 CO2 Calcite 1.106e-004 8.680e-005 1.182e-004 CaCO3 Chalcedony -1.084e-004 -

327、1.473e-004 -6.906e-005 SiO2 Biotite1.370e-0051.317e-0051.370e-005KMg3AlSi3O10(OH)2 Plagioclase1.758e-0041.582e-0041.935e-004Na0.62Ca0.38Al1.38Si2.62O8Redox mole transfers:Sum of residuals (epsilons in documentation): 5.574e+000Sum of delta/uncertainty limit: 5.574e+000Maximum fractional error in ele

328、ment concentration: 5.000e-002Model contains minimum number of phases.=Summary of inverse modeling:Number of models found: 2Number of minimal models found: 2Number of infeasible sets of phases saved: 20Number of calls to cl1: 62例例 1717蒸发作用的反向模拟蒸发作用的反向模拟蒸发模拟与其它多相反响的反向模拟一样。为了模拟蒸发或稀释，H2O中必须包含一种物质。 模拟蒸发

329、的重要原理是包括在每一个反向模拟问题公式中的水摩尔平衡方程 见“反向模拟方程和数值方法 。初始溶液中水的摩尔数乘以它们的混合分数， 再加上通过相溶解或沉淀得到或是失去的水以与由于氧化复原反响得到或是失去的水之后必须等于最终溶液中水的摩尔数。 这个等式是近似的， 因为它不包括多相水解或配位反响中得到或是失去水的摩尔数。表表5050例例1717输入数据的设置输入数据的设置TITLE Example 17.-Inverse modeling of Black Sea water evaporationSOLUTION 1 Black Sea water units mg/L density 1.01

330、4 pH 8.0 # estimated92 / 102 Ca 233 Mg 679 Na 5820 K 193 S(6) 1460 Cl 10340 Br 35 C 1 CO2(g) -3.5SOLUTION 2 Composition during halite precipitation units mg/L density 1.271 pH 5.0 # estimated Ca 0.0 Mg 50500 Na 55200 K 15800 S(6) 76200 Cl 187900 Br 2670 C 1 CO2(g) -3.5INVERSE_MODELING -solution 1 2

331、-uncertainties .025 -range -balances Br K Mg -phases H2O(g) pre Calcite pre CO2(g) pre Gypsum pre Halite preEND这个例子使用到的是Carpenter1978提出的黑海水蒸发的数据。选取了两个分析，即初始黑海水，以与在岩盐沉淀过程中蒸发阶段的黑海水组分。假设方解石、石膏、和岩盐的蒸发和沉淀， 与二氧化碳的溶解足够用来解释水组分中的所有主要离子和溴化物的变化。 输入数据组表50包含在SOLUTIONSOLUTION数据块的溶液组分中。溶液中的总碳未知，但可通过假设这两种溶液都与大气中的二氧

332、化碳平衡来估计。INVERSE_MODELINGINVERSE_MODELING关键字定义了这个例子的反向模拟。溶液2，岩盐沉淀期间的溶液，93 / 102是从溶液1黑海水演化来的。所有数据的不确定围均为2.5%。指定水、方解石、二氧化碳、石膏和岩盐为潜在的反响物-phases-phases。每一种相都必须为沉淀，也就是说，它必须从任何有效的反向模拟的液相中除去。缺省状态下， 反向模拟公式包含水、 碱度和电子的摩尔平衡方程。 另外， 在缺省状态下，包括指定相中所有元素的摩尔平衡方程。 在该例中，缺省值包含了钙、 碳、硫和氯的摩尔平衡方程。 balancesbalances标识符用来指定附加的溴

333、、 镁和钾的摩尔平衡方程。 在缺乏碱度数据时，这些溶液所计算的碱度通过pH值的选择来控制， 并且假定， 溶液与大气中的二氧化碳处于均衡状态。对合理的pH值而言，碱度对电荷平衡的影响很小。在反向计算中仅发现一种模型。这个模型说明黑海水溶液 1必须被浓缩 88 倍才能生成溶液 2，在反向模拟输出中列出了这两种溶液的比值分数表 51 。这样，大约 88kg的黑海水能够浓缩成 1kg 的溶液 2 水。 在蒸发过程中， 生成岩盐 19.75mol 和石膏 0.48mol沉淀。注意，这些摩尔转移是相对于88kg 的水。为了得到每千克黑海水的损失，有必要除以溶液 1 所占的比值。 结果是从每千克黑海水中除去

334、了54.9mol 的水， 0.0004mol 的方解石，0.0004mol 的二氧化碳，0.0054mol 的石膏，和 0.22mol 的岩盐。 这个计算能够从溶液 2得到溶液 1 来，只需要将矿物的约束条件为从沉淀改为溶解。 所有的其他离子镁、钾和溴都保持在所指定的 2.5%的不确定围。反向模拟说明，在给定的不确定围，蒸发失去水 ，二氧化碳脱气，与方解石，岩盐和石膏的沉淀足以解释这两种溶液中主要离子组分发生的变化。表表5151例例1717选择输出结果选择输出结果Solution 1: Black Sea water Input Delta Input+Delta pH 8.000e+000

335、+ 0.000e+000 = 8.000e+000 Alkalinity 8.625e-004 + 0.000e+000 = 8.625e-004 Br 4.401e-004 + 0.000e+000 = 4.401e-004 C(-4) 0.000e+000 + 0.000e+000 = 0.000e+000 C(4) 8.284e-004 + 0.000e+000 = 8.284e-004 Ca 5.841e-003 + 0.000e+000 = 5.841e-003 Cl 2.930e-001 + 7.845e-004 = 2.938e-001 H(0) 0.000e+000 + 0.0

336、00e+000 = 0.000e+000 K 4.959e-003 + 1.034e-004 = 5.063e-003 Mg 2.806e-002 + -7.016e-004 = 2.736e-002 Na 2.544e-001 + 0.000e+000 = 2.544e-001 O(0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+000 + 0.000e+000 = 0.000e+000 S(6) 1.527e-002 + 7.768e-005 = 1.535e-002Solution 2: Composition during ha

337、lite precipitation Input Delta Input+Delta94 / 102 pH 5.000e+000 + -9.369e-014 = 5.000e+000 Alkalinity -9.195e-006 + 0.000e+000 = -9.195e-006 Br 3.785e-002 + 9.440e-004 = 3.880e-002 C(-4) 0.000e+000 + 0.000e+000 = 0.000e+000 C(4) 7.019e-006 + 0.000e+000 = 7.019e-006 Ca 0.000e+000 + 0.000e+000 = 0.00

338、0e+000 Cl 6.004e+000 + 1.501e-001 = 6.154e+000 H(0) 0.000e+000 + 0.000e+000 = 0.000e+000 K 4.578e-001 + -1.144e-002 = 4.463e-001 Mg 2.353e+000 + 5.883e-002 = 2.412e+000 Na 2.720e+000 + -4.500e-002 = 2.675e+000 O(0) 0.000e+000 + 0.000e+000 = 0.000e+000 S(-2) 0.000e+000 + 0.000e+000 = 0.000e+000 S(6)

339、8.986e-001 + -2.247e-002 = 8.761e-001Solution fractions: Minimum Maximum Solution 1 8.815e+001 8.780e+001 8.815e+001 Solution 2 1.000e+000 1.000e+000 1.000e+000Phase mole transfers: Minimum Maximum H2O(g) -4.837e+003 -4.817e+003 -4.817e+003 H2O Calcite -3.802e-002 -3.897e-002 -3.692e-002 CaCO3 CO2(g

340、) -3.500e-002 -3.615e-002 -3.371e-002 CO2 Gypsum -4.769e-001 -4.907e-001 -4.612e-001 CaSO4:2H2O Halite -1.975e+001 -2.033e+001 -1.901e+001 NaClRedox mole transfers:Sum of residuals (epsilons in documentation): 1.947e+002Sum of delta/uncertainty limit: 7.804e+000Maximum fractional error in element co

341、ncentration: 2.500e-002Model contains minimum number of phases.=Summary of inverse modeling:Number of models found: 1Number of minimal models found: 1Number of infeasible sets of phases saved: 6Number of calls to cl1: 2295 / 102例例 1818MadisonMadison 含水层的反向模拟含水层的反向模拟在这个例子中，反向模拟包括同位素摩尔平衡模拟，应用于模拟蒙大纳州美国

342、州名Madison地区含水层中水的演化。 Plummer等1990应用摩尔平衡模拟来定量计算整个含水层的脱白云石作用。在脱白云石作用过程中，硬石膏的溶解引起方解石沉淀和白云石溶解。由摩尔平衡模拟所确定的附加反响包括硫酸盐复原，阳离子交换，以与岩盐和钾盐的溶解Plummer等, 1990。 C和 S数据用来证实摩尔平衡的模拟，C-14用来估计地下水的年龄Plummer等, 1990。初始和最终水样从沿怀俄明州中北部向东北穿过蒙大纳地区的流动路径中选取Plummer等, 1990，流动路径3。这对水样是经过特别选取的，因为只有很少几对水样能表现出以前的摩尔平衡方法和包括不确定围的PHREEQC摩尔

343、平衡法之间相对的较大差异；许多样品对使用这两种方法模拟的结果差异不大。 另外，选取这对样品主要是由于Plummer等已对其进展了详细模拟，并得出摩尔平衡模拟对不同模拟假想的敏感性，且被作为NETPATH手册的例子Plummer等, 1994，例6。PHREEQC的计算结果与NETPATH的计算结果进展了比照。同时Parkhurst1997也讨论了这个例子。表表52.52. 例例1818中溶液的分析数据中溶液的分析数据电荷平衡单位是毫克当量/千克水。其它数据的单位是毫摩尔数/每千克水，除了 C、 S和 C。Fe(2)1413341334亚铁离子，TDIC，总溶解无机碳。 C是相对PDBPee D

344、ee Belemnite标准以千分率表示的TDIC的组分碳-13。 S(6)是相对CDTCaon Diablo Troilite标准以千分率表示的硫酸盐中的硫-34。 S(-2)343413是相对CDT以千分率表示的总硫化物中的硫-34。 C是以百分比目前碳的碳-14的组分。，表示在逆向模拟中所指定的不确定围。pH值的不确定度为0.1，其它数据的不确定度为5%，但铁为100%。-，表示没有测定14分析物分析物温度，PHCaMgNaKFe(2)TDICSO4H2SCl13C34S(6)9.97.551.201.01.02.02.0014.30.160.02溶液溶液1 163.06.61溶液溶液2

345、 211.284.5431.892.54.00046.8719.86.2617.85-2.30.216.31.5-7.01.49.70.996 / 10234S(-2)14C电荷平衡水质成分和反响物水质成分和反响物-52.3+0.11-22.17.0.8+3.24摩尔平衡模拟的初始水溶液1，表52是流动路径3的补给水Plummer等, 1990，重碳酸钙镁水是典型的方解石和白云石地层的补给水。最终的水溶液2，表52是硫酸钠钙水含有较高浓度的氯Plummer等, 1990“Mysse FlowingWell。这种水电荷不平衡，为+3.24 meq/kgw。最终的水也含有适量的硫化物。除了铁离子，

346、指定初始水溶液和最终水溶液所有的化学数据的不确定度围为5%。 指定初始水的不确定度为5%， 主要是因为与最终的水具有一样的流动路径的补给水位置的空间不确定， 对最终的水而言， 由于它接近于最小的不确定限，因此，取得电荷平衡是非常有必要的。由于铁离子浓度较低，其不确定度指定为100%。 pH值的不确定度指定为0.1个单位， 由于在这个取样点存在潜在的二氧化碳脱气，所以是一个保守的估计L.N. Plummer，美国地质调查局，written commun., 1996。从初始水到最终水， C值升高-7.0到-2.3， S值也有所增加9.7到16.3。初始溶液同位素值的不确定度设置为流动路径3和4中

347、4个补给水样中同位素组分围的一半Plummer等, 1990表52。同样，最后水中同位素值的不确定限设置为流动路径 3样品中同位素组分围的一半Plummer等, 1990表52。Plummer等1990考虑的反响物为白云石，方解石，硬石膏，有机物质CH2O，针铁矿， 黄铁矿， Ca/Na2阳离子交换， 岩盐， 钾盐， 和二氧化碳气体。 在其敏感性计算中， Mg/Na2阳离子交换和甲烷是潜在的反响物。这个含水层对于二氧化碳来说是封闭体系，也就是说，不能从外界得到二氧化碳，也不会从相中失去，也不可能得失甲烷Plummer等, 1990。因此，二氧化碳和甲烷气体在PHREEQC摩尔模拟中不作为反响物

348、考虑。二氧化碳气体包含在NETPATH模拟中，但是通过调整硬石膏中的 S，其摩尔转移将会减少为0。溶解一样位素组分的不确定度由Plummer等1990提出，并做了很小的改动，如下：白云石的C，1到5， 有机碳的 C， 为-30到-20； 硬石膏的 S， 11.5到15.5。 沉淀方解石中的 C取决于溶液中同位素的演化，且受溶液同位素分馏的影响。分馏方程不包含在PHREEQC中，因此，有必要来假定方解石组分的围，它代表了方解石沉淀的平均同位素组成。NETPATH计算的方解石沉淀的平均同位素组成大约为-1.5Plummer等, 1990，选取1.0的不确定度来表示分馏系数的不确定度。所有C-14的

349、模拟都是使用PHREEQC模拟中的摩尔转移用NETPATH完成的。沉淀黄铁矿的 S估计为-22，不确定度为2Plummer等, 1990；灵敏度分析说明黄铁矿沉淀的同位素值对摩尔转移具有很少的影响。PHREEQC的输入数据见表53。 注意钾盐， CH2O， 和Ca0.75Mg0.25/Na2交换反响的logK值在PHASESPHASES和EXCHANGE_SPECIESEXCHANGE_SPECIES数据块中都设为0。每一种反响物的化学计量都是正确的， 在摩尔平衡模拟中也需要； 然而，任何饱和指数或是正向模拟在使用这些反响时将是不正确的，因为logK值没有正确定义。97 / 102341334

350、1313341334表表 5353例例 1818 输入数据的设置输入数据的设置TITLE Example 18.-Inverse modeling of Madison aquiferSOLUTION 1 Recharge number 3 units mmol/kgw temp 9.9 pe 0. pH 7.55 Ca 1.2 Mg 1.01 Na 0.02 K 0.02 Fe(2) 0.001 Cl 0.02 S(6) 0.16 S(-2) 0 C(4) 4.30 -i 13C -7.0 1.4 -i 34S 9.7 0.9SOLUTION 2 Mysse units mmol/kgw t

351、emp 63. pH 6.61 pe 0. redox S(6)/S(-2) Ca 11.28 Mg 4.54 Na 31.89 K 2.54 Fe(2) 0.0004 Cl 17.85 S(6) 19.86 S(-2) 0.26 C(4) 6.87 -i 13C -2.3 0.2 -i 34S(6) 16.3 1.5 -i 34S(-2) -22.1 7INVERSE_MODELING 1 -solutions 1 2 -uncertainty 0.05 -range -isotopes 13C 34S -balances98 / 102 Fe(2) 1.0 ph 0.1 -phases D

352、olomite dis 13C 3.0 2 Calcite pre 13C -1.5 1 Anhydrite dis 34S 13.5 2 CH2O dis 13C -25.0 5 Goethite Pyrite pre 34S -22. 2 CaX2 pre Ca.75Mg.25X2 pre MgX2 pre NaX Halite SylvitePHASES Sylvite KCl = K+ + Cl- -log_k 0.0 CH2O CH2O + H2O = CO2 + 4H+ + 4e- -log_k 0.0EXCHANGE_SPECIES 0.75Ca+2 + 0.25Mg+2 + 2

353、X- = Ca.75Mg.25X2 log_k 0.0END摩尔平衡计算包括反响相中所有元素 列于-phases-phases标识符之下 与 S的方程。 NETPATH计算包括同位素分馏方程，用来计算最终水的 C，而PHREEQC计算包括 C的摩尔平衡方程。PHREEQC结果中的修正浓度原始数据加上S使用NETPATH再次运行以取得C-14的年龄，并考虑方解石沉淀对其分馏的影响。NETPATH计算使用电荷平衡选项以确定电荷平衡误差的影响。 NETPATH电荷平衡选项通过分数f修正所有阳离子元素的浓度， 通过分数131334I修f正阴离子的浓度以取得溶液的电荷的平衡。NETPATH的电荷平衡的选

354、项在版本2.13中得到改良以取得准确的电荷平衡；以前的版本仅能取得近似的电荷平衡。对所有的NETPATH计算而言包括使用PHREEQC中修正后浓度，二氧化碳作为潜在反响相， 但是调整硬石膏的 S为不产生二氧化碳摩尔转移。 在其不确定度调整白云石和有机物质中的 C，使其值尽可能达到最终溶液中的 C。MadisonMadison 含水层的结果和讨论含水层的结果和讨论由摩尔平衡模拟所决定的主要反响是脱白云石作用， 离子交换，岩盐溶解，和硫酸盐的99 / 102131334复原反响， 正如下面在表中讨论的不同模拟选项。 脱白云石化作用的主要驱动力是硬石膏的溶解大约为20 mmol/kgw，表54，这将

355、导致方解石沉淀和白云石溶解。一些来自硬石膏溶解的钙和 或 白云石溶解的镁发生离子交换， 释放钠到溶液中。 大约有岩盐溶解。有机质引起的硫酸盐和氢氧化铁复原反响导致产生黄铁矿沉淀。表表54.54. MadisonMadison含水层的摩尔平衡结果含水层的摩尔平衡结果如果没有注明，结果以毫摩尔每千克水数表示。 C以碳-14现在碳的百分比pmc表示， C，碳-13的千分率PDB， S是硫-34的千分率CDT。CH2O代表有机物质。矿物质量转移为正数表示溶解；负数表示沉淀。对交换反响而言，正数表示溶液中钙和或镁浓度降低以与钠浓度升高。-，表示模拟中不包括该反响物。在所有模拟中黄铁矿中的 S大约为-22

356、。为了比照计算的同位素值；测量的 C，-2.3；测量的 S总硫酸盐加硫化物，15.8；测定的 C，0.8pmc14341334341413mmol/kgwCa/Na2结果Mg/Na2Ca0.75Mg0.25/NaNEPATH C电荷平衡-7.6-5.3-12.322.54.31.0-1.015.82.5.03.813,.00012.55.0-20.0-4.315.9PHREEQC-7.7-5.4-12.122.53.5.8-.815.32.5-3.812,90013.45.0-20.0-3.316.0NEPATH ANEPATH BPHREEQC BNEPATH C8.3-3.5-5.320.

357、1.8.1-.115.32.5.012.522,70015.63.6-25.0-2.315.8-8.311.8-21.820.1.8.1-.115.32.5.0.6-2,20015.61.0-30.0-2.215.8-7.711.2-23.922.94.11.0-1.015.32.5-.4-540012.83.0-21.4-3.016.1-8.3-5.6-9.420.1.8.1-.115.32.5-05.916,50015.61.9-25.0-2.315.8Ca/Na2交换Ca0.75Mg0.25/Na交换Mg/Na2交换白云石CaMa(CO3)2方解石(CaCO3)硬石膏(CaSO4)CH2

358、O针铁矿FeOOH黄铁矿(FeS2)盐岩NaCl钾盐KCl二氧化碳CO214C，反响调整外观上年龄年 S，硬石膏 C，白云石 C，CH2O计算的 C，最终水计算的 S，3413131334100 / 102最终水Plummer等1990意识到交换反响的化学计量没有很好地定义，并考虑了在的摩尔平衡模拟的灵敏度分析中这些反响的两个变化。认为纯Ca/Na2交换和纯Mg/Na2交换是潜在的反响物NETPATHA和B，表54。当PHREEQC用这些反响物运行时，将会发现Mg/Na2交换模型PHREEQCB，但是没有发现纯Ca/Na2交换模型。NETPATH和PHREEQC结果的差异来自于溶液中电荷的不平

359、衡。溶液2表52的电荷不平衡为3.24 meq/kgw，这相对于阳离子和阴离子当量浓度和，超过了3%。这不是异常的百分比误差， 但是这个毫克当量的绝对数量与摩尔平衡模拟的一些摩尔转移相比，其值相对较大。当包含电荷平衡约束条件时，通过使用PHREEQC中修订的摩尔平衡方程，使用纯Ca/Na2交换作为仅有的交换反响，不可能使每一种溶液同时取得元素和同位素的摩尔平衡， 达到电荷平衡， 并在指定的不确定度围保存不确定度项。发现的含有大量钙组分交换反响的模拟是关于Ca0.75Mg0.25/Na2(PHREEQCC)。然后这个交换反响应用在NETPATH中以找到NETPATHC。NETPATHC是通过使用

360、NETPATH中的电荷平衡选项以与来与NETPATH C中是一样约束条件来计算的。无电荷平衡选项(NETPATHA,B,和C)的NETPATH模拟与PHREEQC模拟的一致性差异是有机物质的氧化量与针铁矿和黄铁矿的摩尔转移比在PHREEQC模拟中的大。这种差异是由于摩尔转移对电荷平衡的影响。 电荷平衡误差经常表现为单组分反响物的摩尔转移不正确， 例如二氧化碳或有机物质Plummer等，1994。摩尔转移的差异除了有机质，针铁矿，和黄铁矿外，Mg/Na2模拟是相似的NETPATHB和PHREEQCB。然而，这两种模拟都暗含了C-14的年龄为负，这是不切实际的，这一点Plummer等1990也注意

361、到了。PHREEQC模拟与纯Ca/Na2交换模拟(NETPATHA)最相似的是关于Ca0.75Mg0.25/Na2模拟(PHREEQCC)。这个模拟具有比Ca/Na2模拟较大的碳酸盐矿物和有机质的摩尔转移，这降低了修正反响种C-14的活度，产生的地下水年龄为12900PHREEQCC，与22700NETPATHA相比更为年轻。在年龄计算中的变化是由于反响中涉与到的碳的不同。 有两种影响，交换中的反响变化以与电荷平衡差的修正。 交换反响中变化的影响能够从下面两种模拟的差异中估计得到，NETPATHA包含纯Ca/Na2交换，NETPATHC包含Ca0.75Mg0.25/Na2交换，但在溶液组分中没

362、有任何一种模拟矫正了电荷不平衡。 交换反响中Mg的升高可导致方解石和白云石中较大的摩尔转移，计算的年龄从22700降低到16500年。电荷平衡误差的估计是根据NETPATHC和C之间的差异得出，这仅仅是因为在NETPATHC中使用了NETPATH电荷平衡选项。电荷平衡溶液产生了大量的有机质和方解石摩尔转移， 计算的年龄从16500年降低到13000年。 摩尔转移和NETPATHC计算的年龄与NETPATHC得出的相似，仅有些微的不同，因为PHREEQC模拟中的不确定项用来计算得到的不仅是电荷平衡，同样也产生了最接近的最终溶液中的 C。在PHREEQC中所修订的摩尔平衡方程的一个优点是，以前通过

363、建立和运行屡次模拟完成的敏感性分析，现在能够同时包括所有化学和同位素数据的不确定度来实现。 例如，修订一次运行摩尔平衡方程就发现了没有纯Ca/Na2模拟， 即使是任何或所有的化学数据围都加上或13101 / 102是减去10%。这类信息建立前面的摩尔平衡公式非常困难，并且非常消耗时间。另外一个改良是直接包括了电荷平衡条件。 在这个例子中， 包括电荷平衡条件需要改变一次交换反响和调整溶液组分，两者的综合影响是将估计的最下水年龄降低了大约10000年。如果Mg/Na2交换是唯一的交换反响，那么这个年限将是现代的。这样，年限的估计围是非常大的，从0到13000年。然而，由于溶液中钙与镁的比率大约为2.5:1，钙和镁的阳离子交换常数大约是相等的Appelo and Postma, 1993，那么含有占优势的钙的混合交换反响似乎更合理，这将使得更信任老的年龄。 而且， 含水层中C-14年龄和地下水中流动模拟年龄的比照也说明了更老的年龄围具有更合理性。102 / 102