第二讲PHREEQC用户指导_地下水污染与防治

1、PHREEQC第二版是一个用C语言编写的计算机程序，对各种各样低温下的地球化学性质进行了演算。PHREEQC是以离子联系的水化学模型为基础的，可以推算（1）生成物和饱和系数；（2）涉及到可逆反应以及不可逆反应的批反应和一维（1D）运移计算，可逆反应包括水、矿物/无机溶液、气体、固体溶液、表面络合、离子交换平衡；涉及到的不可逆反应包括指定成分摩尔转换、动态控制反应、溶液混合和温度变化；（3）逆向模拟实验，其中多组的无机物和气体摩尔转换以解释在特定成分不确定范围水体之间混合物的不同。和第一版相比，PHREEQC第二版的新特点如下：具有在一维运移计算中模拟弥散（或扩散）和滞流区的能力，用用户确定的速

2、率表达式模拟分子反应，模拟标准的多种成分或非标准的两种成分的固体溶液的沉淀和溶解，模拟定体积气相和定压力气相，考虑表面系数或交换位置随着无机物的溶解和沉淀或者分子反应的变化而变化，自动采用多套收敛参数，打印用户指定量到原始输出文件和（或）适合输入到扩展表格的文件上，以一种与扩展表程序更兼容的形式确定溶解成分。这个版本报告中说明了化学平衡、动态平衡、运移计算以及逆向模拟计算基础方程式，描述了程序输入，以及举例说明了许多程序功能目的和范围这个报告的目的是说明PHREEQC程序的理论和操作，包括组成成分方程的确定，转换这些方程为数值计算方法的解释，补充数值方法计算机代码组织的描述，程序输入描述，以及

3、一系列数据组输入和许多论证的程序功能的数据输入和模拟结果的说明。生成方程和正向模拟在报告的这一部分，说明了用于确定水样热动力活动，离子交换物质、表面络合物质、气相成分、固体溶液和纯相的代数方程式。首先，说明了水样、交换种类和表面性质的热动力活动和质量作用方程。然后，定义一组函数，用f表示，他们必须能同时求解以确定给定条件下的平衡，许多这样的方程是各种元素或元素的价电子状态、交换位置和表面位置的摩尔平衡方程，或来自于纯相和固体溶液的质量作用方程。附加函数用于确定水的碱度、活度、电荷平衡、气相均衡、离子强度和表面络合均衡。每个函数包含最少量的变量，使函数的个数等于变量的个数。程序用一种修改过的牛顿

4、罗弗逊方法解答了伴随的非线性方程。这种方法利用了函数的偏差和每个函数关于主要未知数组或主要未知数的偏导数。为了清晰起见，微分计算中用的变量组指主要未知数，每个函数的全导数不用求导说明。在下列方程中，没有下标或下标“（aq）”指液相中的实体，“（e）”指交换，“（g）”指气体，“（s）”指表面，"（ss）“指固体溶液，“（p）”指相的状态。活度方程和质量作用方程在这一部分，定义了水的活度、交换过程和表层物质，并描述了每种物质的质量作用联系。质量作用方程是化学系统中根据主要未知数每种物质的摩尔数的质量作用表达式。然后，求这些方程关于主要未知数的差分。每种物质的摩尔数方程和偏导数方程将被组

5、成的摩尔平衡，电荷平衡和相均衡函数取代。水样PHREEQC允许单个液相的物质形成或物质平衡。然而多个液相可在一个运行过程中定义，一个液相可被定义为一种或多个液相的混合物（见数据输入说明”中的MIX关键词）。一个例外是液相中溶解物质假定处于热动力平衡状态。在最初的溶解计算中，允许氧化还原元素的化学价不平衡。每个水样i的未知数包括活度ai,活度系数，重量摩尔浓度mi,溶液中摩尔数ni。PHREEQC根据主要物质改写了所有的化学方程式。其中一个与每个元素（例如，钙离子Ca+2）或元素化学价态（例如，三价铁离子Fe+3）有关的主要水样，加上氢离子活度、水中电子活度和水的活度。一些程序，例如MINTEQ

6、A2（Allision和其他，1990）和MINEQL+（Schecher和Mcavoy,1991）,把这些物质用术语“成分”表示，但是，这个术语在这儿不能用，因为易与Gibbs'的相规则中的“成分”混肴。PHREEQC中，每个水样主要物质用SOLUTION-MASTER-SPECIES数据块确定（见“数据输入说明”）。数值方法将未知数的个数减少为主要未知数的最小个数，并且反复迭代这些未知数的值，直到得到数学方程组的解。水溶液的主要未知数是主要物质活度的自然对数、水的活度的自然对数Ct”门、离子H2O强度科、水溶液中有溶解能力的水的质量Waq。下列关系适用于所有水样（除了水成电子和水本

7、身）：ai=rimi和ni=miwaq.离子缔合模型中水样之间的均衡要求满足所有的水样质量方程。例如，水样CaSO；的缔合化学反应是Ca2+SO42-=CaSO4°。这个反应在25°C时的logK是2.3,形成质量作用方程：_02.3CaSO410二:2口2Ca2+SO4-(1)般情况下，质量作用方程式可被写为:Maq一._c.Ki=：i|I：-mmJ其中Ki是温度均衡常数，Cm是物质i中主要种类m的化学计量系数，Maq是水中主要物质总数，Cm，i的值可正可负。在PHREEQC中，缔合化学反应右边各项系数为负，左边各项系数为正。主要物质采用同样的形式，质量作用方程可简化为1

8、=am/am。一种水样的总摩尔数i可由质量作用表达式求得:MaqIMni=miWaq=KiWaq一i(3)牛顿.罗弗逊方法用摩尔数关于主要未知数的全微分。计算如下:Maq-dni=nid1n"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方程

9、中，a°是离子大小参数，反之在WATEQDebye-Huckel方程中，a和bi是适合于mean-salt的离子特定参数。否则，如果没有在数据库文件或输入数据组中说明，Davies方程用于带电物质。对于不带电物质，活度系数方程中第一项为零，WATEQDebye-Huckel方程简化为Setchenow方程（lnri=bw）（见Langmuir,讨论版1997）。否则，如果没有指定，对所有不带电物质，bi假定为0.1。Davies方程中，这些活度系数方程关于离子强度的偏微分为:ln;=-ln(10)1Azi21/出6+12-0.3(7)对于扩展的WATEQ或Debye-Huckel方程

10、为:22d-Az2(8)ln彳=-ln(10)zz-2+b护&何Bq0"+1)PHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式，logK和它的温度要求，每种水样的活度系数参数通过SOLUTION-SPECES数据块确定。元素的主要种类和元素的化合价态用SOLUTION-MASTER-SPECES数据块确定。一种溶液的组成用SOLUTION或SOLUTION-SPREAD数据块确定（见“数据输入说明”）。交换种类离子交换平衡通过交换位置的非齐次质量作用方程和摩尔平衡方程包括在模型中。PHREEQ比许多种交换与液相在均衡中共存，用“交换集合”表示。这种方法使用了质量作

11、用表达式，它以每个交换器中水样和虚构的空的交换位置之间的不完全反应为基础（Appeio和Postma1993）。这种空的交换点是交换器的重要种类，其活度的对数是一个附加的主要未知数。它由EXCHANGE-MASTER-SPECES数据块确定(见“数据输入说明”)。然而，交换器的主要物质不包括在摩尔平衡方程中，使其物质浓度为零。它的活度在物质上是无意义的，但是如此却使所有的交换位置充满了其他交换物质。交换计算的未知数是活度在PHREEQC中定义为当量比值乘以活度系数兀和交换器e中每种交换物质ie的摩尔数nieo当量比值是被一种交换物质占据的交换点的摩尔数比交换点的总数。一种交换物质的活度为6=丁

12、%上，其中be,ie是被交换物质ie占据-eTe的交换器e的当量数，Te是交换器中交换地点的当量总数。Te作为系统中交换器交换地点的当量总数，不必与每克水的当量数相等，因为系统中水的质量可能多或少于1kgo根据默认值，交换物质的活度系数为1.0,但是可以根据水的离子强度和被交换物质占据的交换地点当量数选择使用Davies,扩展的Debye-Huckel或WATEQDebye-Huckel活度系数。水和交换物质之间的平衡要求满足所有交换物质的质量作用方程。交换物质CaX2的缔合反应是Ca2+2X-=CaX2,其中X是默认数据库的交换主要物质。活度的当量比值和这种化学反应形式的使用统称为Gaine

13、s-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,工2Ca2X-(9)般情况写，质量作用方程可写为:Mcm,ieKie=5口«m,(10)m其中m随所有的主要物质变化，包括交换的主要物质，cm，ie是在物质ie不完全缔合反应中主

14、要物质m的化学计量系数，kie是不完全反应选择常数。cm，ie的值可正可负。PHREEQC中，缔合反应右半边各项系数为负，左半边各项系数为正。对于交换物质，物质ie总摩尔数方程为:(11)在数值方法中，交换器主要物质活度的自然对数是主要未知数。物质ie摩尔数关于主要未知数的全微分是;Mdni=niee(12)cmcm,iedln(o(m)-£ln(X)dN.mJPHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式，logK和它的温度要求以及每种交换物质的活度系数表达式通过EXCHANGE-SPECES数据块确定。交换的主要种类用EXCHANGE-MASTER-SPECES数据

15、块确定。交换位置数量和交换器组成用EXCHANGE数据块确定(见“数据输入说明”)。表面种类表面络合过程通过非齐次质量作用方程、表面位置的摩尔平衡方程和每个表面的电荷-电势联系包括在模型中。PHREEQC允许多种表面和表面位置类型液相在均衡中共存，用“表面集合”表示。PHREEQC表面种类质量作用方程有两个公式：(1)一个包括静电电势项；(2)另一个包括所有的静电电势项。如果使用模型静电电势项，由于存在表面电荷和表面静电电势，就需用到附加方程和质量作用术语。交换反应和表面反应计算之间的两个根本不同点在于，交换反应表示为不完全反应，使主要物质不能出现在任何摩尔平衡方程中，并且交换物质被期望是中性

16、的(或不带电的)。表面反应不是不完全反应，所以主要物质是物质上真正存在的物质，并且可出现于摩尔平衡方程中，表面物质可带正电，也可带负电，或不带电。在Dzombak和Morel(1990)中，表面络合反应的基本理论包括静电电势。该理论假定活性地点个数为Ts(eq),比面积As(m2/g),表面质量Ss(g)已知。两个附加的主要未知数是（1）参量lna乎s=lne2RT=，其中F是Faraday常数（96493.5JV-1eq-1）,2RT<JWs是表面电势（Volts）,R是气体常数（8.3/4）Jmol-1k-1）,T是温度（kelvin）o（2）主要表面物质活度的自然对数。参量lna.

17、s右侧项在术语标准中用a2定义。这是用于Dzombak和Morel（1990）中一个不同的未知数，但是由于所有的方程书写和这个主要未知数一致，所以结果相同。表面物质活度假定等于给出非空的表面位置类型的摩尔比值。换句话说，一种表面物质完全覆盖了一种给出的表面位置时是处于标准状态（活度为1）。这个协议不同于Dzombak和Morel（1990）,它假定一个表面物质（固相概念上）的活度数值上等于体积克分子浓度（溶液浓度）。假如只考虑monodeneatecomplexes（如Dzombak和Morel（1990所做）,消去质量作用方程中各项，不考虑标准状态协议，可得到一致的数值结果。然而，体积克分子

18、浓度（容摸）协议用于multideneatecomplexes（bidentate,tridentite,其他，cf.Appeloandpostma,1999）,表面位置浓度存在一个明显的不同。假如一个容器盛有表面包含multideneatespecies且处于均衡的溶液，并且增加了更多同样的溶液，体积克分子浓度不变，溶液的成分和表面不断改变。这种情况体积克分子浓度明显不正确。“Hfo”含水Fe2O3）在默认数据库中用作“_w；指明了较低的亲和力或很弱的位置，“揩明强的亲和力或强的位置。“Hfo_wOH”用于说明在软弱位置”一种中性表面物质，带负电软弱位置形成物的缔合反应（在一定意义上，在缔合

19、反应中定义的种类在方程的右侧）可以被写为：Hfo_wOH-Hfo_Wo-+H+.（13）质量作用表达式，包括静电电势项项，是ota+(14)int_Hfo_wOH十rtHfowOe-OfLXHfo_wOH其中是反应原有的平衡常数，eRT是一个因子，说明了将一种带电物质（H+）从带电表面移出。一般，表面物质i,）的质量作用方程为：K；ntiskaiskMCm,iskFPs【二me二；2mJRT(15)其中Kit是原有的平衡常数，是表面s中表层位置类型k（软弱的或坚固的Dzombak和Morel,1990中）的第i层物质；m随所有主要物质M变化，包括表层主要物质；cmi,表层物质）合、缔合反应主要

20、物质的化学计量系数,z,是形成表层物质时表面电荷的净变sk化。cmi可正可负。PHREEQC中，缔合反应中右边半边各项系数为负，左半边各项系数,sk为正。一种表层物质i（sk）的总摩尔数方程为:二:工二KReskiskbi飞kiskTc,M0=Kiiskbiskm(16)其中Tsk是表面位置类型总数，bi）是与物质相联系表面位置的数量。物质i（sk）的摩尔数关于主要未知数的全微分为:-Mdni=n.卜Cmidln：m-2zisksk1mskmi(17)表面物质质量作用方程的第二个公式在质量作用表达式中不包括静电电势项（-no-edlidentifier在SURFACE数据块中）。一种表层物质摩

21、尔数方程同方程16,除了涉及“收的因式没有显示。否则，摩尔数的全微分同方程17,除了没有最后一项。PHREEQC的数据输入，摩尔平衡化学方程式和质量作用表达式,logK和它的表面物质的温度要求通过SURFACE-SPECES数据块确定。表层主要物质或表层位置类型通过SURFACE-MASTER-SPECES数据块确定。表层的一致性，每个位置类型的当量数,表面组成，比表面积和表层质量用SURFACE数据块确定（见“数据输入说明”）。气相成分多成分气相和液相之间的均衡可用非齐次质量作用方程和一个总压力方程（仅为故乡压力气相）模拟。只有一个气相可与液相在均衡中存在，但气相可包括多种成分。所有的气体成

22、分假定为理想气体，气相假定为理想的气体混合物。如果气相体积固定，那么该体积内的压力将随着反应程度变化，但是气相中每种成分总是存在的。对于一个固定体积的气相，不需要另外的主要未知数，气体中一种成分的摩尔数可以从液体主要物质的活度中推算出来。如果气相压力固定，那么气相体积将随着反应程度变化，假如组成成分的部分压力之和小于指定的总压力，固定压力气相将不存在，并且气相中任一气体成分将不能存在。固定压力气相中，气体成分的总摩尔数方程中包括一个附加的主要未知数Ngas。根据理想假定，一种气体成分的逸度（或活度）等于其部分压力。PHREEQC使用溶解方程式，在一定程度上，气体成分假定处于化学反应的左边。CO

23、2溶解反应可写为CO2（g）=CO2（aq）。（18）亨利定律常数将气体成分的部分压力（理想气体数值上等于逸度）与液体物质的活度联系了起来。对于二氧化碳，亨利定律常数为10-1.468假定为理想气体，单位为大气压（atm）,下列质量作用方程适用于均衡状态：1.468,、PCO210"co2mqX，（19）22aq其中，Pco2是用液相的活度计算的部分压力（atm）。一般情况下，气体成分的部分压力可根据液相活度写为：1 Maqpg=n<m，g，（20）Kgm其中，Pg是气体成分g的部分压力，用液相的活度计算；Kg是气体成分的亨利定律常数；cm,g是溶解方程中水成主要m的化学计量系

24、数。Cm,g值可正可负。PHREEQC中，溶解反应方程中右半边各项系数为负，左半边各项系数为正。对于定体积气相，气相的总体积指定为Vtotal,但压力是变化的。在平衡状态，气体成分的摩尔数如下计算:_VtotalPg_Vtotal-cmng=ZTLZ-=mRTRTKgm,g(21)气相中该成分的全微分为MVaqdng詈工ngCm,gdln二mRTm(22)对于定压力气相，气相的总压力指定为Ptotal,但体积是变化的。在平衡状态，气体成分的摩尔数等于气体总压力的比值乘以气体的总摩尔数，如下计算:ng=NgasPgNgasPotalPtotalKMaqcnm-g：-mgm(23)气相中该成分摩尔

25、数的全微分为dngPMaq广dN(24)gas'、ngCm,gdln;PtotalmPHREEQC的数据输入、质量作用方程、亨利定律常数以及该常数的温度要求通过PHASE数据块确定。气相类型（定体积或定压力）、气相计算中包括的成分和初始气象成分通过CAS_PHASE数据块确定（见“数据输入说明”）。牛顿-罗弗逊方程用一组函数（用f表示）用于说明复杂的均衡，这些方程主要来源于用摩尔平衡方程和电荷平衡方程代换物质摩尔数方程（来源于前一部分的质量作用方程）。达到平衡时，有关待定平衡计算的所有函数值为零。函数的零点通过Tventm-Raohsm方法获得。据此每个函数关于每个主要未知数微分形成T

26、acobin矩阵。由Tacobin矩阵得到一组线性方程组，其解可近似为非线性方程的解。通过反复解答连续的线性方程组，可得到非线性方程的一个解。在这一部分说明了每个用于数值方法的函数f及其关于形成Jacobian矩阵的主要未知数的全导数水的活度水的活度用根据Raoult定律(Gaoult'sandChrist,1965P65)得出的近似值计算:-：H2ONaqn=1-0.017%.iWaq(25)函数fH2O定义为:fH2O-Waq-:H2ONaq-10.017%ni.(26)该函数的全微分为:NaqdfH2O=Waq：H20dln：h2O;H2O-1Waqdln叫0.017d.(27)

27、主要未知数是水的活度的自然对数ln0tH2O。离子强度水溶液的离子强度是一个主要未知数，定义如下:Naq2：斌,(28)函数f科定义为:1Naq2f-二5”,2i(29)该函数的全微分为:1Naq2df'-WaqdlnWaqNqdJ.三z2dn.2i(30)定体积多项成分气相平衡对于定体积气相，每种气体成分的摩尔数可由水中主要物质的活度计算得到，数值模型处理气相成分的方法与它处理水中物质的方法不同。每种气体成份的摩尔数ng用在摩尔平衡方程中用成分形式表示，dng在Jacobian矩阵中用摩尔平衡方程表示。定体积气相均衡计算不需要附加方程foPHREEQC数据输入，质量作用方程，Henr

28、y定律常数，以及其温度要求用PHASES数据块确定，气相类型（定体积或定压力），气相计算中包括的成分和初始气相成分由GAS_PHASE数据块确定（见“数据输入说明”）。定压力多成份气相平衡定体积气相中每种气体成份的摩尔数由水中主要物质的活度和气相中气体成份总摩尔数Ng决定。每种气体成份的摩尔数ng在摩尔平衡方程中用（elements）成份表示，dng在Jacobian矩阵中表示摩尔平衡方程。定压力多成份气相和液相平衡需要一新的方程气体成份部分压力之和等于总压力Ptotal。函数fR洞定义为:NgfPtotal=Pto-%Pg,g(31)其中，Ng是气相中气体成份的总摩尔数。fPtotal关于主

29、要未知数的全微分（其中正的dNgas表示溶液浓度增加）为：NgMaqdfMl""cm,gPgdln:gm(32)PHREEQC的数据输入，质量作用方程，Henry定律常数以及该常数的要求用数据块PHASES确定。气相类型(定体积或定压力)，气相计算中包括的成份和初始气相成份用数据块GAS-PHASES数据块定义(见“数据输入说明”)纯相均衡气相和液相之间的平衡，包括具有固定部分压力的气体，通过非齐次质量作用方程包括在模型中。PHREEQC允许多种纯相(用纯相集合表示)和液相在均衡中存在，受Gibbs纯相规则的限制，纯相的活度假定一致为1.0。每种纯相的附加未知数为系统中存在

30、的纯相的摩尔数np，其中p指第p个纯相，np表示每个纯相摩尔数的变化，在摩尔平衡方程中用成份(elements)表示。PHERRQC也允许如下平衡计算，通过增加或清除一种特定的反应物形成具有纯相的平衡(EQUIL旧RIUM_PHASES数据块中可供选择的公式或相)；计算形成纯相平衡必须的反应物摩尔换算。在这种计算模型中，摩尔平衡方程中各项来自反应物的化学计算而不是纯相的化学计算，未知数是进入或从溶液中析出的反应物摩尔数。与每个新未知数相应的新函数是每个纯相的质量作用表达式。PHREEQC中的溶解反应在一定意义上，纯相是化学方程式左半边的部分。方解石的溶解反应可写为：CaCO3=Ca2+CO32

31、-,(33)并且，用10-8.48的logK和纯固态活度1.0,质量作用表达式为Kca-=10'.48=七2二CO-(34)一般情况下，纯相均衡可用下述方程表示：Maq_cKp个一，.m(35)其中Cm,p是溶解反应中主要物质m的化学计量系数，其值可正可负。PHREEQC中溶解反应左边各项系数为负，右边各项系数为正。无机物SIp饱和指数定义为MaqSlp-logn.mm,p.(36)m数值方法中相均衡函数是Maqfp=QnKp+in(10虎骨吸卜£Cm,pln(«m)(37)m其中Slp.target是相的目标饱和指数，ln(10)将以10为底的对数转换为自然对数。

32、目标饱和指数由用户指定；正值、零、负值分别表示溶液中无机物过饱和、均衡、不饱和。对固定部分压力气体成分，SIp,target等于气体成分部分压力的对数，关于主要未知数的全微分是Maqdfp=-ZCm,pdlnc(m.(38)mPHKEEQC的数据输入，质量作用方程，均衡常数和纯相的温度要求用数据块PHASES确定。纯相集合的初始成分和目标饱和指数用数据块EQUILIBRIUM-PHASES确定。固体溶液均衡理想多成分或非理想，双成份固体溶液均衡的模拟是以Glynn的研究为基础的(Glynn和Readon,1990;Glynn和其他人，1990;Glynn,1991;GlynnandParklm

33、st,1992)。液相和固体溶液均衡通过非齐次质量方程包括在模型中。PHREEQC允许多种固体溶液和液相在均衡中存在，用固体溶液集合表示，受Gibbs相规则的限制，非理想固体溶液的模拟只限于双成份固体溶液；理想固体溶液可以有两种或多种成份。固体溶液的附加主要未知数是固体溶液中每种成份的摩尔数npss,其中ss指固体溶液ss。该项说明了在元素摩尔平衡Jacobian矩阵中每种成份摩尔数的变化。不象纯相，固体溶液中一种成份的活度不完全是1.0。其定义为a二九Xp,其中xppsspsspss,pss是固体溶液ss中成份p的摩尔比值，Kpss是活度系数，固体溶液中一种成份的摩尔比值定、，np义为xps

34、skA工np-psspss1,其中Nss是固体溶液ss中成份的摩尔数，对于理想固体溶液，活度系数为1.0,对于非理想双成份固体溶液中各成份的活度系数用2、入1=exp(：0-：1(4xi-1)x2)和入2=exp(X0+CC1(4x2-1)xi2),Gugenheim表达式确定:(39)(40)其中入1和入2是成份1和成份2的活度系数，”和a1是无因次的Gugenheim参数，无因次参数由过量自由能g0和g1(kJ/mol)的量纲参数根据方程ct0=g0/RT和/1=91小丁计算得到。过量自由能参数必和a2可以直接定义，也可能通过许多方法间接确定，包括界定可混合性范围的成份2的摩尔比值，界定s

35、pinodal范围的成份2的摩尔比值，临界点和临界温度处的成份2的摩尔比值、Thompson和Waldbaum参数、Mangules参数，成份2的摩尔比值和一个alyotropic点的所有溶解物质的对数，成份1和2微量浓度的固相活度系数，或者成份2的两个分布(配给)系数(glynn,1991)。相应与每个新未知数的新函数是每个固体溶液中每种成份的质量作用表达式。PHREEQC使用溶解反应，在一定意义上固定溶液成份指化学反应方程式左边各项。对于霞石(文石)的文石-菱铜矿的固体溶液，溶解反应可写为CaCO3=Ca2+CO32-,(41)KArag=10-8.34-Ca2"CO；一：Ca2

36、:CO-CaCO3(42)-,AragnAragArag1nAragnnStront)般情况下,固体溶液相均衡中，每种成份可用下列方程说明:使用10-8.43的logK和固体活度系数，形成质量作用表达式为：Maq(43)c-m,PssI1-mssKp:ppsspss其中Kd是纯相中成份p的均衡常数，cmD是固体溶液ss中成份p的溶解反应中，主要物pssmi,pss质m的化学计量系数，其值可正可负。PHREEQC中溶解反应左边各项系数为负，右边各项系数为正。固体溶液一种成份的溶解部分（商）可定义为Maqcm,pss(44)QpKpss：pss其中平衡时Qp等于1,inQp等于0。pssPss对于

37、非理想双成份固体溶液，其数值方法中每种成(45)(46)Maqdf1八Cm,1dln二mm222c23、+匚至+2口0乂2-6口诙+12。泾dn+、R口+n2)2：0x22：0x26：1x2-18：1x212：1x231dn2rhn2(47)分的函数为Maqf1="cm1ln二m-inK1-in-m'nn2Maqn2f2-%cm2in：m-inK2-inin'2.m,Rn2关于主要未知数的全微分是df2;加刈“m-2"2二。><2-6中18-6:-12二mnn2)dn2-2,-2一3一x142u0Xi+631x112%x1+1n2(48)理想固体

38、溶液中，在数值方法中每种成份的函数为ZMaqfPss=lnQd=lnpssCm，PssJ，,-mmk7pss一ln(49)其中NssNtotal=£nss,且jSS包括固体溶液ss中所有成份。jss其关于主要未知数的全微分是Maqdf='Cm,pssdIn二m-pssssmnNtotal-'nNtotaldnPssss，jss印ss1dnjjssNtotaljss(50)PHREEQC的数据输入，质量作用方程，平衡常数,没仲春相常数的温度要求用PHASES数据块定义。一种固体溶液的初始成份和非理想固体溶液的GUYYENHESIN参数用SOLID_SPLUTIONS数据

39、块确定（见“数据输入说明”表面层的摩尔平衡表层处摩尔平衡是一般摩尔平衡方程的特殊例子。表层集合是由一个或多个表层组成,每个表层有一个或多个位置类型。一个表层位置类型的总摩尔数由输入指定为如下之一:（1）固定值，（2）与一个纯相的摩尔数成正比，或（3）与动态反应物摩尔数成正比。下列函数来源于表面S的一个表面位置类型8k的摩尔平衡关系。一个位置类型的表层物质占据的表层位置的摩尔数总和必须等于表层位置类型的总摩尔数。卜列函数来源于表面S的一个表层位置类型Sk的摩尔平衡关系:NskfSk=Tsk_工bsk.ifknifk)，1sk(51)其中当达到平衡时，函数fsk的值为零，丁与是表层位置类型的摩尔数

40、，Nsk是某种位置类型表层物质的数量，bsk,i（%）是表层物质i（sk卢据的表层位置的数量，fsk的全微分为dfsk=ATsk-Zbsk/kM（52）kkkskskisk如果表面位置总数和纯相的摩尔数成正比例，那么ATsk=-Csk,pdnp,其中Csk,p是每摩尔p相表层位置的摩尔数。如果相溶解，dnp为正并且表面点减少。如果点的总数和动态反应物质成比例，在全微分方程ATSk=0中。点数的变化作为反应的一部分并与速率方程相结合，且Jacobia矩阵中不包括任何项。随着动态反应增加或减少，反应物质摩尔数、表面点数与此成比例调整。假如表面点数固定，ATSk=0。PHREEQC数据输入，表面点每

41、种类型的摩尔总数用SURFACE数据块确定，且可能是一个固定的量或者一个纯相的摩尔数或动态反应有关。表面点类型用surface_species数据块确定并且表层物质可用surface_species数据块确（见“数据输入说明”）交换处的摩尔平衡一个交换点的摩尔平衡是一般摩尔平衡方程的特例。一个交换点的总摩尔数由输入指定为如下之一：（1）固定值，（2）与一个纯相的摩尔数成正比，或（3）与一个动态反应物摩尔数成正比。被交换物质占据点的摩尔数之和必须等于交换点的总摩尔数。下列函数是一个交换点的摩尔平衡关系：Nefe=Tebe,ienie,ie（53）其中达到摩尔平衡时，函数fe值为零，Te是交换器e

42、中交换点的总摩尔数，口%是交换物质占据的交换点的数量。fe的全微分方程为Nedfe=ATe-Zbe,iednie.（54）ie如果点的总数与一个纯相的摩尔数成比例，那么ATe=-Ce,pdnp,其中Ce，p是每摩尔p相交换点的摩尔数，如果相溶解，那么dnp为正并且交换点数减少。如果点的总数与动态反应物摩尔数成比例，在全微分方程中ATe=0。点数的变化与作为反应部分与速率方程相结合，且没有相包括在Jacobian矩阵中。随着动态反应增加或减少，反应物摩尔数、交换点数成比例增减。如果交换点数固定,Ne&又=0。fe=Te-Zbe,ienieiePHREEQC数据输入，交换点摩尔数在EXCH

43、ANGE数据块中确定且可能是一个定量或者与一个纯相或一个动态反应物白摩尔数有关。交换点用EXCHANGE_SPECIES数据块确定且交换物质用EXCHANGE_SPECIES数据块确定(见“数据输入说明”)碱度摩尔平衡碱度摩尔平衡方程仅用于物质形成计算和逆向模拟中。它是一般摩尔平衡方程的特例，系数由每个水样的碱性确定。PHREEQC中碱性被确定为一种成分且一种主要物质与其相联系(见”数据输入说明"中SOLUTION_MASTER_SPECIES关键词)。在PHREEQC的默认数据库中，碱性主要物质是CO32-O碱度的主要未知数是lnaAlk，或默认数据库中lnaco32-.碱度当量总

44、数由模型输入确定，每种水样的碱性影响之和必须等于碱度当量总数。下列函数是碱度平衡方程：NaqfAlk=TAlk-ZbAlk,ini,(55)ie其中达到摩尔平衡时，函数fAlk的值为零，TAlk是溶液中碱度当量数，bAlk,I是水样i的碱度影响，fAlk全微分方程为NaqdfAlk二rbAlk,idn.(56)ie假定碳酸盐是主要碱性物质，TAlk的值必须为正。从概念上讲，用PHREEQC计算的碱度不同于实测的碱度。在默认PHREEQC的数据库文件中，bAlkj的值已选定，每种元素的所指状态(bAlki=0)或元素化合价态是PH=4.5时的主导物质。假定所有的元素或化学化合价态在碱性滴定中被转

45、化为这种主导物质。然而所指状态不存在的水样的有效浓度（具有非零碱度的物质）可能在滴定终点存在，这在一定程度上使用PHREEQC计算的碱度与所测碱度是一个不同的量。Fe和Al的氢氧化物络合物是最普通的例子，可能不被转换为所指状态。这样，一种溶液用PHREEQC计算的碱度尽管在数值上等于所测定碱度，实质上也是一个近似值，因为是假定滴定过程最终将元素和元素化合价态转变为它们的所指状态。在许多溶液中，碱度主要来源于碳酸盐，近似值是有效的。PHREEQC的数据输入，以每种物质的缔合反应计算的碱度在SOLUTION_SPECIES数据块中确定，主要物质的碱度由SOLUTION_MASTER_SPECIES

46、数据块确定，总碱度是溶液组成的一部分。用SOLUTIONorSOLUTION_SPREAD数据块确定（见“数据输入说明”）元素摩尔平衡包括纯相、固体溶液系统中一种元素的总摩尔数是下列各相中初始存在的摩尔数总和,集合、液相、交换集合、表面集合、气相和表面扩散层。卜列函数是总摩尔平衡方程。npfm1Tmbm,pnpSKsN%SSNssbm,pnssPssNgpssNaq-、bm,EiENefbm,ienieeieSNaq(57),sbm,isknisk,一bm,gng-bm,ini；k其中，达到摩尔平衡时，函数fm值为零，Tm是系统中成份的总摩尔数，Np是纯相集合中相的数量，SS是固体溶液集合中固

47、体溶液数量，Nss是固体溶液ss中成份的数量，Naq是水样数量，E是交换集合中交换器的数量，Ne是交换位置e的交换物质，S是表层集合中表层个数，Nsk是表层类型sk的表层物质数量，Ng是气相成份数量。系统中每个实体的摩尔数分别用np表示纯相集合中相的摩尔数，npss表示固体溶液的成分数，ni表示水样数，me表示交换位置e的交换物质数，门求示表层位置类型sk物质数，ng表示气体成份数，ni表sks示表面s扩散层中水样数。每个实体中元素m的摩尔数用bm表示，用一个附加的下标说明相关实体，bm通常但不总是等于Cm（质量作用方程中主要物质m的系数）。为了避免解决大数之间的微小差别，方程57中圆括弧中的

48、量没有明确包括在求解算法中，且Tm的值没有实际计算。函数fm中使用代替量NpSSNss*Tm=Tm-Zbmpnp-ZZbmpnp。最初，*Tm由液相m、交换集合、表面集合、气,pp,psspsspsspss相和表面扩散层中总摩尔数计算得到：NFNSKNQ1NSN。(58)aqFeSsskgSaq*Tm=£bm,iQ-Z工bm,ienie+£££bm,ifkbm,gng+ZZbm,in,sieieski%gsi在方程的反复迭代中，通过纯相和固体溶液成份的摩尔转换使*Tm得到更新：(59)NpSSNss*Tkmi="Tkm+£bm,pdn

49、p十££/加。,，psspss其中k指迭代数，在中间迭代中*Tm可能是负值，但达到均衡时必须为正。函数fm的全微分是NpSSNssNaqFNedfmbm,pdnp'二1口小小,-、bm,idn二二口仙e-(60)psspssieieSKsNskNgSNaq二.，二bm,isdni-'bm,gdngbm,idni,sskikgsiskPHRFFQC的数据输入，元素的总摩尔数最初用SOLUTIONorSOLUTION_SPRFAD数据块确定，对于交换集合，用EXCHANGE数据块确定，对于表面集合，用SURFACE数据块，对于气相，用GAS_PHASE数据块，纯

50、相集合中每个相的摩尔数用EQUILIBRIUM_PHASES数据块确定。在一个固体液体集合中每种固体溶液每种成份的摩尔数SOLID_SOLUTIONS数据块确定。元素的总摩尔数也可在批反应和迁移计算中改修动（见“数据输入说明”）。溶液中电荷平衡电荷平衡方程计算了水中阴阳离子的当量总和。在一些情况下，表面和交换器电荷不平衡。当具体指定时，在初始溶液计算中用电荷平衡方程调节PH值或一种主要物质的活度（继而一种元素的总浓度或元素的化合价态）使溶液达到电中性。电荷平衡方程对于在批量反应和迁移模拟中计算PH值是必要的。在现实的溶液中，阴阳离子的当量和必须为零。然而在化学分析中分析误差和没有分析的成份通常

51、引起溶液计算中电荷不平衡。如果一个初始溶液中计算到电荷不平衡，在后来的批量反应和运移模拟中需调节PH值，维持同样的电荷不平衡。将溶液混合，批反应的电荷不平衡是每种溶液其混合因素产生电荷不平衡偏差之和。如果模拟中用到一个表面，且没有指定显扩散层计算，那么带电的表层物质形成将产生表面电荷不平衡。相似地，如果交换物质不是电中性（在默认数据库中所有交换物质是电中性）,那么交换器将积累电荷。一般的电荷平衡方程包括表面和交换器电荷不平衡。在每个初始溶液和批反应步骤中，对于每个单元在运移模拟的每个时间区间中计算电荷不平衡使用方程NaqTz,q=zzn,(61)其中q代表液相，Tza是液相q的电荷不平衡，Zi

52、是水样i所带电荷。如果带电的表层和交z,q换器不存在，在溶液批反应中或运移模拟终点的电荷不平衡将和模拟始点相同。表面的电荷不平衡在初始表面组成计算和每个批反应步骤中计算，对于每个单元运移模拟的每个时间区间，使用方程KsNskNaqTz,sZzipfkJg4ni,s,kiski(62)其中Tz，s是指表面电荷不平衡以、是表面s的表面类型Sk的表层物质i所带电荷,方程中最后一项代表弥散层中积聚的电荷。最后一项仅用于弥散层组分明确包括在计算中时（在数据块SURFACE中的-diffuse-layer）。详细计算弥散层组分时，要求所有的溶液电荷平衡，并且Tz，s总为零。正常情况下，交换物质没有静电荷，

53、但是对于大部分，这不作要求。然而，如果带电物质的和不等于交换点（交换能力）的当量总数，交换物质（当量积分）的活度不能很好确定。如果带电物质存在，那么在每个批反应步骤中的初始交换组分计算中计算交换器的电荷不平衡，对于每个单元，运移模拟的每个时间区间用方程NeTz,e='、3eiNpSSNssNaqENedfm=A*pdnp;二bm,pssdnpss-均,idn=：b.,iednie-pssPssieieSKsNSkNgSNaq二：二二bm,idni-bm,gdng-'、bm,idni,sSkiskskskgsi(63)其中Tz，e是交换器的电荷不平衡,4是交换器e中交换物质i所带

54、电荷。系统电荷不平衡批反应步骤之初确定，并且对于每个单元，在运移模拟中每个时间段始点应该是:(64)QSETz=£«qTz,q+£Tz,s+£Tz,e,qse其中Tz是系统电荷不平衡，Q是在每个批反应步骤或一个运移步骤单元中混合的液相数量,气液相q的混合比值，S是表层的数量，E是交换器数量。电荷平衡函数是:NaqSKsNskfz=Tz-£zni-£££,*k)SNaqWS.ssiENe-Z£3，eie(65)其中达到电荷平衡时fz为零。如果计算弥散层组份，每个表面包括一个独立的电荷平衡方程,达到电荷平衡时，括号内项目的和将为零。如果弥散层组份不计算，括号内每二项为零。fz的全微分方程为：NaqSKsNskENedfz二八4ni二二二4sknisk-J4nieiskia,eiaske(66)其中，只有在弥散组份没有详细计算时，表面三项的总和存在