6.2.2.1　模型

已有的研究（郑西来等，1990；陈芸等，1994；王锦国等，1999）表明，水—岩系列间相互作用的规律遵从最小自由能原理、质量守恒原理、电子及电荷守恒原理。对于由多个反应构成的水—岩系列间的相互作用系统，可用如下模型来描述。

1.目标函数

由化学热力学的基本原理知，一个化学反应的自由能变化ΔGr 越小，正反应趋势则越大。对于由多个反应构成的水—岩相互作用系统，可建立如下的目标函数：

式中：Xi为第i 个物理化学反应发生的量（以mmol/L表示）。当Xi＞0 时，表示反应正向进行；当Xi＜0 时，表示反应逆向进行；当Xi＝0 时，则表示反应在所研究的系统中可能未发生。

2.约束条件

由质量守恒原理知，系统中地下水从A 到B 点，其化学组分的变化必须满足：

上游（A点）化学组分＋反应相＝下游（B 点）化学组分＋生成相

因此有：

式中：αik为第i 个物理化学反应、反应相矿物中第k 种组分的化学计量数；βik为相应生成相矿物中第k 种组分的化学计量数。令bik＝αik－βik，式（6.2.32）可写为

对于系统中所发生的氧化还原反应，必须满足电子守恒。即有

而对系统中的阳离子交换反应，必须满足如下电中性方程

式中：ri为第i 种发生阳离子交换作用的离子化合价。

可将式（6.2.32）～式（6.2.35）合并为如下形式：

式中：m为约束条件方程个数。

这样，由目标函数式（6.2.31）和约束条件式（6.2.36），构成了水—岩相互作用的基本模型。

上述数学模型对于描述水—岩系列间相互作用的复杂系统仍然具有普遍的意义。并且可以得出，运用化学热力学方法来研究水—岩作用具有两大优点：①可以不考虑参与化学反应的物质结构；②可以不考虑化学反应的途径和过程，而是通过反应过程的始态和到达平衡后的终态，进行热力学推导和计算，来探求化学反应的程度和方向。因此，可以根据对系统的认知来推知系统状态的变化。

6.2.2.2　求解方法

1.确定性方法

如上述，已有多位学者应用上述水—岩作用的线性规划模型进行地下水地球化学模拟，并取得了一些富有意义的成果。其共同的特点是，根据单次样本分析结果或根据多次样本分析的结果采用将数据平均化后，输入水—岩作用模型，应用最优化方法（如单纯形方法等）求解，作为输出结果就是目标函数（反映体系中化学反应的速度及反应方向）及有关矿物的溶蚀量等具有确定性的数值。此类求解方法就是所谓的确定性方法。显然，把此类方法应用于求解水—岩作用模型，忽略了水—岩作用系统具有不确定性，即具有随机性的特点。

2.随机性方法

人们在实践中逐步认识到，地下水系统具有随机性，水—岩作用系统亦具有随机性或不确定性。这种不确定性反映在多个方面：如水—岩系列间水文地球化学作用过程大多具有非线性，可能出现化学波以及混沌等现象；与溶质运移转化有关的如弥散性能、阻滞吸附性能、降解性能以及氧化还原环境变化的不确定性；与初始条件有关的如溶质初始浓度、衰减系数及其变化的不确定性；与溶质之间相互作用有关的如溶出产生阶段和运移转化阶段的不确定性；与尺寸效应有关的如在小范围内测试的数据应用在大范围中产生的不确定性等。这种不确定性除了反映在客观方面，还体现在主观方面，如建立概念模型中出现的不确定性、模型求解误差，以及参数测试过程产生的不确定性等（宋汉周，2003）。

针对水—岩作用系统存在的如上所述的不确定性，宜采用随机分析方法进行探讨。目前，在水文地质领域中常用的随机分析方法，有原序分析法、扰动分析法、可靠度分析法以及蒙特卡罗分析法等。后一种方法在实际工作中得到了比较广泛的应用。

实际上，蒙特卡罗（Monte Carlo）分析法是一种统计试验的方法，即为一种求取工程技术问题近似解的数值方法（Smith等，1981）。其求解过程一般包括如下步骤。(www.daowen.com)

（1）随机变量的抽样试验。根据基本随机变量的已知概率分布进行随机抽样，即输入随机变量进行数字模拟。

（2）样本反应求解。对每个抽取的样本，按照问题的性质采用确定性的控制数学、物理方程求解样本反应。

（3）计算反应量的统计量估计。对于所有的样本反应，按待求解的类型分别求出输出随机变量的均值、方差或概率分布。

研究表明，蒙特卡罗分析法既适用于求解确定性问题，也适用于求解随机性问题。

当求解确定性问题时，其基本思路可归纳为：首先，根据所提出的问题构造一个简单且适用的概率模型，使问题的解对应于该模型中随机变量的某些数字特征（如概率、数学期望及方差等）；然后，在计算机上生成随机数，并对生成的随机数进行统计分析试验；最后，利用试验所获结果求出具有统计特征的估计值作为问题的近似解。

当求解随机性问题时，其基本思路亦可归纳为：首先，根据问题的物理性质建立随机模型；然后，根据模型中各个随机变量的分布，在计算机上产生随机数，进行大量的统计试验，以取得所求问题的大量试验值；最后，根据这些试验结果求其统计特征值，从而获得所求问题的解。由此可见，采用蒙特卡罗方法求解随机性问题的步骤与求解确定性问题的步骤基本一致。

蒙特卡罗法也存在缺点。其中之一，就是计算成本高，其理论误差为ε～σ/n1/2。因此，通常需要生成大量的样本数n 才能保证一定的精度。对于计算高阶矩和多维概率分布，必须产生大量的数据才能使结果有意义。

蒙特卡罗分析法仍是目前模拟输入变量不确定性对输出变量不确定性最有效的模拟方法，只要已知输入随机变量的概率函数，即可通过随机模拟计算输出随机变量的概率分布。蒙特卡罗法能够模拟输出变量的整体概率特征，它适用性广，既可用于解析模型，又可用于数值模型；既可用于有限域，又可用于无限域，还可用于刻画任意程度的输入不确定性。

在采用随机性分析方法求解水—岩作用模型之前，需要论述随机数的生成技术，因为随机数的生成技术是进行随机模拟的基础。

设X 为随机变量，概率分布为F（·）。由于F（·）是非减函数，且值域为［0，1］，其逆函数为F－1（·），在［0，1］上有定义。设u 是［0，1］上均匀分布的随机变量，则有

那么变量x ＝F－1（u）的概率分布函数必为F（·）。由此，要产生服从分布为F（·）的随机数，就只要产生［0，1］上的均匀的随机数u，然后由F－1（·）求得。这种产生随机数的方法称为逆变换法。例如，设F 为参数为1 的指数分布，即

对于求出的F 的反函数F－1（y）＝－ln（1－y），若u是［0，1］上均匀的随机变量，则F－1（u）＝－ln（1－u）服从参数为1 的指数分布。

从原理上讲，逆变换法可以求得服从任何连续分布的随机数，然而在有些情况下反函数无解析表达式以至于计算很复杂，此时可采用如下方法。

对于均匀分布：u（a，b），设随机变量x 服从［a，b］上的均匀分布，其概率密度函数为：

记为u（a，b）。均匀分布的随机数是构造具有其他形式分布的随机数的基础，可以由计算机附加的硬件来生成，但不具有再生性。通用的方法是由计算机按一定的算法生成从统计意义而言满足独立性及均匀性的一串数字，这样的数称为伪随机数。

用算法生成伪随机数，至少应满足：算法容易执行，且生成的速度快；产生的数字序列能通过独立性及均匀性的统计检验。

目前，同余法常用于产生伪随机数。令

式中：a 为正整数，称为乘子；c 为非负整数，称为增量；0≤x1≤m，x1称为种子，m称为模数，同时也是伪随机数序列的长度。这样，对任何一个初始值x1，可由上式产生一个序列｛x1，x2，…，xn｝。通过下式，就可以产生区间［0，1］上的均匀分布的随机数

等价地，［a，b］上的均匀分布的随机数可以由

产生。由以上各式可知模m越大，随机数的均匀性越好。从计算的角度来讲，取m＝2l最方便，其中l 为正整数，例如取a＝27＋1，c＝1，m＝235。

对于正态分布：N（u，σ2），设随机变量x 服从正态分布，其概率密度函数

对于式（6.2.36）中呈随机变化的变量bj的每一个实现，水—岩作用模型应该总是成立的。这里，采用Monte Carlo方法，随机地从变量bj的分布中抽取一个实现［b′1，b′2，…，b′m］，利用单纯形法产生符合水—岩作用模型｛式（6.2.31）和式（6.2.36）｝的一组解X′i＝［x′1，x′2，…，x′n］i。如此用同样的方法，产生N 组解，当N 足够大时，可以统计得出各矿物溶蚀量的特征值。其均值和标准差分别为：