第二讲PHREEQC用户指导地下水污染与防治.docx
《第二讲PHREEQC用户指导地下水污染与防治.docx》由会员分享，可在线阅读，更多相关《第二讲PHREEQC用户指导地下水污染与防治.docx（83页珍藏版）》请在冰点文库上搜索。
第二讲PHREEQC用户指导地下水污染与防治
PHREEQC第二版是一个用C语言编写的计算机程序，对各种各样低温下的地球化学性质进行了演算。
PHREEQC是以离子联系的水化学模型为基础的，可以推算
（1）生成物和
饱和系数；
（2）涉及到可逆反应以及不可逆反应的批反应和一维（1D）运移计算，可逆反
应包括水、矿物/无机溶液、气体、固体溶液、表面络合、离子交换平衡；涉及到的不可逆反应包括指定成分摩尔转换、动态控制反应、溶液混合和温度变化；（3）逆向模拟实验，其
中多组的无机物和气体摩尔转换以解释在特定成分不确定范围水体之间混合物的不同。
和第一版相比，PHREEQC第二版的新特点如下：
具有在一维运移计算中模拟弥散（或扩散）和滞流区的能力，用用户确定的速率表达式模拟分子反应，模拟标准的多种成分或非
标准的两种成分的固体溶液的沉淀和溶解，模拟定体积气相和定压力气相，考虑表面系数或
交换位置随着无机物的溶解和沉淀或者分子反应的变化而变化，自动采用多套收敛参数，打
印用户指定量到原始输出文件和（或）适合输入到扩展表格的文件上，以一种与扩展表程序
更兼容的形式确定溶解成分。
这个版本报告中说明了化学平衡、动态平衡、运移计算以及逆向模拟计算基础方程式，描述了程序输入，以及举例说明了许多程序功能
目的和范围
这个报告的目的是说明PHREEQC程序的理论和操作，包括组成成分方程的确定，转
换这些方程为数值计算方法的解释，补充数值方法计算机代码组织的描述，程序输入描述，以及一系列数据组输入和许多论证的程序功能的数据输入和模拟结果的说明。
生成方程和正向模拟
在报告的这一部分，说明了用于确定水样热动力活动，离子交换物质、表面络合物
质、气相成分、固体溶液和纯相的代数方程式。
首先，说明了水样、交换种类和表面性质的热动力活动和质量作用方程。
然后，定义一组函数，用f表示，他们必须能同时求解以确定给定条件下的平衡，许多这样的方程是各种元素或元素的价电子状态、交换位置和表面位
置的摩尔平衡方程，或来自于纯相和固体溶液的质量作用方程。
附加函数用于确定水的碱度、
活度、电荷平衡、气相均衡、离子强度和表面络合均衡。
每个函数包含最少量的变量，使函数的个数等于变量的个数。
程序用一种修改过的牛顿・罗弗逊方法解答了伴随的非线性方程。
这种方法利用了函数的偏差和每个函数关于主要未知数组或主要未知数的偏导数。
为了清晰
起见，微分计算中用的变量组指主要未知数，每个函数的全导数不用求导说明。
在下列方程
中，没有下标或下标“（aq）”指液相中的实体，“（e）”指交换，“（g）”指气体，“（s）”
指表面，"（ss）“指固体溶液，“（p）”指相的状态。
活度方程和质量作用方程
在这一部分，定义了水的活度、交换过程和表层物质，并描述了每种物质的质量作用联系。
质量作用方程是化学系统中根据主要未知数每种物质的摩尔数的质量作用表达式。
然后，
求这些方程关于主要未知数的差分。
每种物质的摩尔数方程和偏导数方程将被组成的摩尔平衡，电荷平衡和相均衡函数取代。
水样
PHREEQC允许单个液相的物质形成或物质平衡。
然而多个液相可在一个运行过程中定义，一个液相可被定义为一种或多个液相的混合物（见数据输入说明”中的MIX关键词）。
一个例外是液相中溶解物质假定处于热动力平衡状态。
在最初的溶解计算中，允许氧化还原
元素的化学价不平衡。
每个水样i的未知数包括活度ai,活度系数「，重量摩尔浓度mi,
溶液中摩尔数ni。
PHREEQC根据主要物质改写了所有的化学方程式。
其中一个与每个元素（例如，钙离子Ca+2）或元素化学价态（例如，三价铁离子Fe+3）有关的主要水样，加上氢离子活度、水
中电子活度和水的活度。
一些程序，例如MINTEQA2（Allision和其他，1990）和MINEQL+（Schecher和Mcavoy,1991）,把这些物质用术语“成分”表示，但是，这个术语在这儿不
能用，因为易与‘Gibbs'的相规则中的“成分”混肴。
PHREEQC中，每个水样主要物质
用SOLUTION-MASTER-SPECIES数据块确定（见“数据输入说明”）。
数值方法将未知数的个数减少为主要未知数的最小个数，并且反复迭代这些未知数的值，直到得到数学方程组
的解。
水溶液的主要未知数是主要物质活度的自然对数、水的活度的自然对数Ct”门、离子
H2O
强度科、水溶液中有溶解能力的水的质量Waq。
下列关系适用于所有水样（除了水成电子和水本身）：
ai=rimi和ni=miwaq.离子缔合模
型中水样之间的均衡要求满足所有的水样质量方程。
例如，水样CaSO；的缔合化学反应是
Ca2++SO42-=CaSO4°。
这个反应在25°C时的logK是2.3,形成质量作用方程：
_0
2.3CaSO4
10二
:
•2口2
Ca2+SO4-
（1）
般情况下，质量作用方程式可被写为:
Maq
•一._c.
Ki=：
i|I：
-mmJ
其中Ki是温度均衡常数，Cm—是物质i中主要种类m的化学计量系数，Maq是水中主要物
质总数，Cm，i的值可正可负。
在PHREEQC中，缔合化学反应右边各项系数为负，左边各项
系数为正。
主要物质采用同样的形式，质量作用方程可简化为1=am/am。
一种水样的总摩尔数i可由质量作用表达式求得:
Maq
[IM
ni=miWaq=KiWaq」一i
（3）
牛顿.罗弗逊方法用摩尔数关于主要未知数的全微分。
计算如下:
Maq-
dni=ni[d1n"G;闻*）.1n（i）"
（4）
水样的活度系数由戴维斯方程确
10gii2——*
（5）
或者用扩展的或Debye-Huckel方程确定:
（6）
logi--
其中z是水样i的离子电荷数，A和B是关于温度的常数。
假如bi为零，方程6是扩展的
Debye-Huckel方程；假如bi不为零，方程6是Debye-Huckel方程（见Truesdell和Jones,1974）。
在扩展的Debye-Huckel方程中，a°是离子大小参数，反之在WATEQDebye-Huckel方程中，a和bi是适合于mean-salt的离子特定参数。
否则，如果没有在数据库文件或输入数据组中说明，Davies方程用于带电物质。
对于不带电物质，活度系数方程中第一项为零，WATEQDebye-Huckel方程简化为Setchenow方程（lnri=bw）（见Langmuir,讨论版1997）。
否则，如果没有指定，对所有不带电物质，bi假定为0.1。
Davies方程中，这些活度系数方程关于离子强度的偏微分为:
■
ln;=-ln（10）1Azi2
1
/出6+12
-0.3
（7）
对于扩展的WATEQ或Debye-Huckel方程为:
22
d-Az2
（8）
——ln彳=-ln（10）—zz-2+b
护&何Bq0"+1）
PHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式，logK和它的温度要
求，每种水样的活度系数参数通过SOLUTION-SPECES数据块确定。
元素的主要种类和元
素的化合价态用SOLUTION-MASTER-SPECES数据块确定。
一种溶液的组成用
SOLUTION或SOLUTION-SPREAD数据块确定（见“数据输入说明”）。
交换种类
离子交换平衡通过交换位置的非齐次质量作用方程和摩尔平衡方程包括在模型中。
PHREEQ比许多种交换与液相在均衡中共存，用“交换集合”表示。
这种方法使用了质量作
用表达式，它以每个交换器中水样和虚构的空的交换位置之间的不完全反应为基础（Appeio
和Postma1993）。
这种空的交换点是交换器的重要种类，其活度的对数是一个附加的主要未
知数。
它由EXCHANGE-MASTER-SPECES数据块确定（见“数据输入说明”）。
然而，交换器的主要物质不包括在摩尔平衡方程中，使其物质浓度为零。
它的活度在物质上是无意义
的，但是如此却使所有的交换位置充满了其他交换物质。
交换计算的未知数是活度在PHREEQC中定义为当量比值乘以活度系数兀和交
换器e中每种交换物质ie的摩尔数nieo当量比值是被一种交换物质占据的交换点的摩尔数比交换点的总数。
一种交换物质的活度为6=丁%上，其中be,ie是被交换物质ie占据-eTe
的交换器e的当量数，Te是交换器中交换地点的当量总数。
Te作为系统中交换器交换地点的
当量总数，不必与每克水的当量数相等，因为系统中水的质量可能多或少于1kgo根据默认
值，交换物质的活度系数为1.0,但是可以根据水的离子强度和被交换物质占据的交换地点
当量数选择使用Davies,扩展的Debye-Huckel或WATEQDebye-Huckel活度系数。
水和交换物质之间的平衡要求满足所有交换物质的质量作用方程。
交换物质CaX2的
缔合反应是Ca2++2X-=CaX2,其中X是默认数据库的交换主要物质。
活度的当量比值和这种化学反应形式的使用统称为Gaines-Thomas协议（Gaines和Thomas,1953）,并且在
PHREEQC版本中用于数据库phreeqc.dat和wateq4f.dat文件。
［在PHREEQC版本中，也可能用Gapon协议，也使用当量比值，但是将交换反应写为0.5Ca2++X-=Ca0.5X。
见Appeio
和Postma（1993）再次讨论。
］在默认数据库文件中，钙离子交换的logK伪.8,形成质量
作用方程：
100.8
工CaX2
.工2,工2
Ca2X-
（9）
般情况写，质量作用方程可写为:
Mcm,ieKie=5口«m,（10）
m
其中m随所有的主要物质变化，包括交换的主要物质，cm，ie是在物质ie不完全缔合反应中
主要物质m的化学计量系数，kie是不完全反应选择常数。
cm，ie的值可正可负。
PHREEQC
中，缔合反应右半边各项系数为负，左半边各项系数为正。
对于交换物质，物质ie总摩尔数方程为:
（11）
在数值方法中，交换器主要物质活度的自然对数是主要未知数。
物质ie摩尔数关于主要未知数的全微分是;
「M
dni=ni
ee
（12）
cmcm,iedln（o（m）-£ln（X）dN.
mJ
PHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式，logK和它的温度要
求以及每种交换物质的活度系数表达式通过EXCHANGE-SPECES数据块确定。
交换的主要
种类用EXCHANGE-MASTER-SPECES数据块确定。
交换位置数量和交换器组成用EXCHANGE数据块确定（见“数据输入说明”）。
表面种类
表面络合过程通过非齐次质量作用方程、表面位置的摩尔平衡方程和每个表面的电荷-
电势联系包括在模型中。
PHREEQC允许多种表面和表面位置类型液相在均衡中共存，用“表
面集合”表示。
PHREEQC表面种类质量作用方程有两个公式：
（1）一个包括静电电势项；
（2）另一个包括所有的静电电势项。
如果使用模型静电电势项，由于存在表面电荷和表面静电电势，就需用到附加方程和质量作用术语。
交换反应和表面反应计算之间的两个根本不同点在于，交换反应表示为不完全反应，使主要物质不能出现在任何摩尔平衡方程中，并且交换物质被期望是中性的（或不带电的）。
表面反应不是不完全反应，所以主要物质是物质上真正存在的物质，并且可出现于摩尔平衡
方程中，表面物质可带正电，也可带负电，或不带电。
在Dzombak和Morel（1990）中，表面络合反应的基本理论包括静电电势。
该理论假
定活性地点个数为Ts（eq）,比面积As（m2/g）,表面质量Ss（g）已知。
两个附加的主要未
知数是
（1）参量lna乎s=lne2RT=，其中F是Faraday常数（96493.5JV-1eq-1）,2RT
Ws是表面电势（Volts）,R是气体常数（8.3/4）Jmol-1k-1）,T是温度（kelvin）o
（2）主要
表面物质活度的自然对数。
参量lna.s右侧项在术语标准中用a2定义。
这是用于Dzombak
和Morel（1990）中一个不同的未知数，但是由于所有的方程书写和这个主要未知数一致，
所以结果相同。
表面物质活度假定等于给出非空的表面位置类型的摩尔比值。
换句话说，一种表面物
质完全覆盖了一种给出的表面位置时是处于标准状态（活度为1）。
这个协议不同于Dzombak
和Morel（1990）,它假定一个表面物质（固相概念上）的活度数值上等于体积克分子浓度
（溶液浓度）。
假如只考虑monodeneatecomplexes（如Dzombak和Morel（1990所做）,消去质量作用方程中各项，不考虑标准状态协议，可得到一致的数值结果。
然而，体积克分子
浓度（容摸）协议用于multideneatecomplexes（bidentate,tridentite,其他，cf.Appeloandpostma,1999）,表面位置浓度存在一个明显的不同。
假如一个容器盛有表面包含multideneatespecies
且处于均衡的溶液，并且增加了更多同样的溶液，体积克分子浓度不变，溶液的成分和表面
不断改变。
这种情况体积克分子浓度明显不正确。
“Hfo”含水Fe2O3）在默认数据库中用作“_w；指明了较低的亲和力或很弱的位置，
“—揩明强的亲和力或强的位置。
“Hfo_wOH”用于说明在软弱位置”一种中性表面物质，带
负电软弱位置形成物的缔合反应（在一定意义上，在缔合反应中定义的种类在方程的右侧）可以被写为：
Hfo_wOH-Hfo_Wo-+H+.（13）
质量作用表达式，包括静电电势项项，是
ota+
（14）
int_Hfo_wOH十rt
HfowO——e
-Of
LXHfo_wOH
其中是反应原有的平衡常数，eRT是一个因子，说明了将一种带电物质（H+）从带电表面
移出。
一般，表面物质i,）的质量作用方程为：
K；nt
isk
aisk
MCm,iskFPs
【二me二；2
mJRT
（15）
其中Kit是原有的平衡常数，是表面s中表层位置类型k（软弱的或坚固的Dzombak和
Morel,1990中）的第i层物质；m随所有主要物质M变化，包括表层主要物质；cmi,表
层物质）合、缔合反应主要物质的化学计量系数,
△z,是形成表层物质时表面电荷的净变sk
化。
cmi可正可负。
PHREEQC中，缔合反应中右边半边各项系数为负，左半边各项系数,sk
为正。
一种表层物质i（sk）的总摩尔数方程为:
「二:
工^二KRe
skiskbi飞k
isk
Tc,M0
=Ki
iskbiskm
（16）
其中Tsk是表面位置类型总数，bi⑼）是与物质相联系表面位置的数量。
物质i（sk）的摩尔数关
于主要未知数的全微分为:
-M
dni=n.卜Cmidln：
m-2z
isksk1mskmi
（17）
表面物质质量作用方程的第二个公式在质量作用表达式中不包括静电电势项（
-no-edl
identifier在SURFACE数据块中）。
一种表层物质摩尔数方程同方程16,除了涉及“收的因
式没有显示。
否则，摩尔数的全微分同方程17,除了没有最后一项。
PHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式,
logK和它的表面物
质的温度要求通过SURFACE-SPECES数据块确定。
表层主要物质或表层位置类型通过
SURFACE-MASTER-SPECES数据块确定。
表层的一致性，每个位置类型的当量数,
表面组
成，比表面积和表层质量用SURFACE数据块确定（见“数据输入说明”）。
气相成分
多成分气相和液相之间的均衡可用非齐次质量作用方程和一个总压力方程（仅为故乡压力气相）模拟。
只有一个气相可与液相在均衡中存在，但气相可包括多种成分。
所有的气
体成分假定为理想气体，气相假定为理想的气体混合物。
如果气相体积固定，那么该体积内的压力将随着反应程度变化，但是气相中每种成分总是存在的。
对于一个固定体积的气相，不需要另外的主要未知数，气体中一种成分的摩尔
数可以从液体主要物质的活度中推算出来。
如果气相压力固定，那么气相体积将随着反应程度变化，假如组成成分的部分压力之和小于指定的总压力，固定压力气相将不存在，并且气相中任一气体成分将不能存在。
固定
压力气相中，气体成分的总摩尔数方程中包括一个附加的主要未知数Ngas。
根据理想假定，一种气体成分的逸度（或活度）等于其部分压力。
PHREEQC使用溶
解方程式，在一定程度上，气体成分假定处于化学反应的左边。
CO2溶解反应可写为
CO2（g）=CO2（aq）。
（18）
亨利定律常数将气体成分的部分压力（理想气体数值上等于逸度）与液体物质的活度联系了
起来。
对于二氧化碳，亨利定律常数为10-1.468［假定为理想气体，单位为大气压（atm）］,
下列质量作用方程适用于均衡状态：
1.468,、
PCO2—10"co2mqX，（19）
22aq
其中，Pco2是用液相的活度计算的部分压力（atm）。
一般情况下，气体成分的部分压力可
根据液相活度写为：
1Maq
pg=nKgm
其中，Pg是气体成分g的部分压力，用液相的活度计算；Kg是气体成分的亨利定律常数；
cm,g是溶解方程中水成主要m的化学计量系数。
Cm,g值可正可负。
PHREEQC中，溶解反应方程中右半边各项系数为负，左半边各项系数为正。
对于定体积气相，气相的总体积指定为
Vtotal,但压力是变化的。
在平衡状态，气体成
分的摩尔数如下计算:
_VtotalPg_Vtotal-cm
ng=~ZTLZ-=■■m
RTRTKgm
g
（21）
气相中该成分的全微分为
M
Vaq
dng詈工ngCm,gdln二m
RTm
（22）
对于定压力气相，气相的总压力指定为
Ptotal,但体积是变化的。
在平衡状态，气体成
分的摩尔数等于气体总压力的比值乘以气体的总摩尔数，如下计算:
ng=Ngas
Pg
Ngas
PotalPtotalK
M
aq
cnm-g
：
-m
gm
（23）
气相中该成分摩尔数的全微分为
dng
PMaq
广dN
（24）
gas'、ngCm,gdln;
Ptotalm
PHREEQC的数据输入、质量作用方程、亨利定律常数以及该常数的温度要求通过
PHASE数据块确定。
气相类型（定体积或定压力）、气相计算中包括的成分和初始气象成分通过CAS_PHASE数据块确定（见“数据输入说明”）。
牛顿-罗弗逊方程
用一组函数（用f表示）用于说明复杂的均衡，这些方程主要来源于用摩尔平衡方程
和电荷平衡方程代换物质摩尔数方程（来源于前一部分的质量作用方程）。
达到平衡时，有
关待定平衡计算的所有函数值为零。
函数的零点通过Tventm-Raohsm方法获得。
据此每个
函数关于每个主要未知数微分形成Tacobin矩阵。
由Tacobin矩阵得到一组线性方程组，其
解可近似为非线性方程的解。
通过反复解答连续的线性方程组，可得到非线性方程的一个解。
在这一部分说明了每个用于数值方法的函数f及其关于形成Jacobian矩阵的主要未知数的全
导数
水的活度
水的活度用根据Raoult定律（Gaoult'sandChrist,1965P65）得出的近似值计算:
-：
：
H2O
Naq
n
=1-0.017%—^.
iWaq
（25）
函数fH2O定义为:
fH2O-Waq-:
^H2O
Naq
-10.017%ni.
（26）
该函数的全微分为:
Naq
dfH2O=Waq：
H20dln：
h2O;H2O-1Waqdln叫0.017d^.
（27）
主要未知数是水的活度的自然对数ln0tH2O。
离子强度
水溶液的离子强度是一个主要未知数，定义如下:
Naq
「2：
斌,
（28）
函数f科定义为:
1Naq2
f-二5”,
2i
（29）
该函数的全微分为:
1Naq2
df'-WaqdlnWaqNqdJ.三z2dn.
2i
（30）
定体积多项成分气相平衡
对于定体积气相，每种气体成分的摩尔数可由水中主要物质的活度计算得到，数值模型
处理气相成分的方法与它处理水中物质的方法不同。
每种气体成份的摩尔数ng用在摩尔平
衡方程中用成分形式表示，dng在Jacobian矩阵中用摩尔平衡方程表示。
定体积气相均衡计
算不需要附加方程fo
PHREEQC数据输入，质量作用方程，Henry定律常数，以及其温度要求用PHASES
数据块确定，气相类型（定体积或定压力），气相计算中包括的成分和初始气相成分由
GAS_PHASE数据块确定（见“数据输入说明”）。
定压力多成份气相平衡
定体积气相中每种气体成份的摩尔数由水中主要物质的活度和气相中气体成份总摩尔
数Ng决定。
每种气体成份的摩尔数ng在摩尔平衡方程中用（elements）成份表示，dng在
Jacobian矩阵中表示摩尔平衡方程。
定压力多成份气相和液相平衡需要一新的方程
气体
成份部分压力之和等于总压力Ptotal。
函数fR洞定义为:
Ng
fPtotal=Pto⑶-%Pg,
g
（31）
其中，Ng是气相中气体成份的总摩尔数。
fPtotal关于主要未知数的全微分（其中正的dNgas表示溶液浓度增加）为：
NgMaq
dfMl""cm,gPgdln:
gm
（32）
PHREEQC的数据输入，质量作用方程，Henry定律常数以及该常数的
要求用数据
块PHASES确定。
气相类型（定体积或定压力），气相计算中包括的成份和初始气相成份用
数据块GAS-PHASES数据块定义（见“数据输入说明”）
纯相均衡
气相和液相之间的平衡，包括具有固定部分压力的气体，通过非齐次质量作用方程包括
在模型中。
PHREEQC允许多种纯相（用纯相集合表示）和液相在均衡中存在，受Gibbs纯相规则的限制，纯相的活度假定一致为1.0。
每种纯相的附加未知数为系统中存在的纯相的
摩尔数np，其中p指第p个纯相，np表示每个纯相摩尔数的变化，在摩尔平衡方程中用成份
（elements）表示。
PHERRQC也允许如下平衡计算，通过增加或清除一种特定的反应物形
成具有纯相的平衡（EQUIL旧RIUM_PHASES数据块中可供选择的公式或相）；计算形成纯相平衡必须的反应物摩尔换算。
在这种计算模型中，摩尔平衡方程中各项来自反应物的化学计算而不是纯相的化学计算，未知数是进入或从溶液中析出的反应物摩尔数。
与每个新未知数相应的新函数是每个纯相的质量作用表达式。
PHREEQC中的溶解反应
在一定意义上，纯相是化学方程式左半边的部分。
方解石的溶解反应可写为：
CaCO3=Ca2++CO32-,（33）
并且，用10-8.48的logK和纯固态活度1.0,质量作用表达式为
Kca-=10'.48=七2二CO-
（34）
一般情况下，纯相均衡可用下述方程表示：
Maq
__c
Kp个一，.m
（35）
其中Cm,p是溶解反应中主要物质m的化学计量系数，其值可正可负。
PHREEQC中溶解反应
左边各项系数为负，右边各项系数为正。
无机物SIp饱和指数定义为
Maq
Slp-logn...mm,p.（36）
m
数值方法中相均衡函数是
Maq
fp=QnKp+in（10虎骨吸卜£Cm,pln（«m）（37）
m
其中Slp.target是相的目标饱和指数，ln（10）将以10为底的对数转换为自
升级会员 