写给初学者的话

简介

预测模型已经是一种地下水科学中的基本工具，绝大多数受过正规训练的地下水工作者都或多或少地接触过地下水溶质运移模型。初学者经常会误认为地下水模型的主要作用是预测水位变化或者污染运移所造成的影响，这至少是不全面的认识。地下水模型在更大程度上是用来认识含水层行为，测试水力控制和溶质运移假说的工具。也就是说，模型的主要功能是理解局部和区域含水层系统的工具。

本文简要介绍了反应和非反应情况下地下水溶质运移的基本理论。首先介绍了水流和溶质运移基本控制方程；其次是对解算方法的简要描述；最后简要介绍了将化学反应耦合进入模型的方法。

控制方程

建立溶质运移模拟的前提是优质的地下水野外勘察。由于地下水流场控制了污染烟羽的迁移和扩散，所以只有首先对地下水流进行模拟才能随后模拟溶质运移。同样道理，不深入了解含水层系统就想准确模拟地下水流，一样是不可能的任务。模拟地下水流的第一要务是定义水流系统，具体来说就是定义含水层性质、定义来自边界的影响、识别系统应力（如抽水和垂向补给）的位置和大小、量化含水层的水力特性和参数的分布。

尽管定义模型需要详尽的实地勘察资料，但不能因此忽视模型在勘察工作早期的作用。地下水工作者可以使用简易模型来检验关于含水层的一些关键假设，比如边界的真实属性、地下水流场形态、垂向补给率的指定、模型对关键物理化学参数的敏感性等。从一项地下水工作的最早期，模型就可以帮助研究人员确定数据的“含金量”，也就是某些特定信息在量化地下水系统和增强模型预测能力中所能发挥的真实作用。模型的构建和完善应该与数据的收集工作并行。

地下水模拟中最常见的关注点是饱和带中的水分运移。对于非均质、各向异性、二维结构的承压含水层，有如下的控制方程：

其中：Tij是导水系数，L2/T；h是水头，L；S是储水系数，无量纲；t为时间，T；W是单位面积上的体积通量，L/T；x，y是平面直角坐标。

所谓地下水流模型，就是使用解析方法或数值方法对上述控制方程求解的过程。水流模型的解算结果就是在给定的初始条件和边界条件下，水头在空间各处和时间各点上的分布。一般使用数值方法对控制方程求解，因为解析方法求解时对初始条件和边界条件有严格的限制，而客观世界中的情况极为复杂，只有数值方法能够处理与之相对应复杂性的模拟问题。地下水模拟的首要工作是准确地描述地下水系统。由于地下水数据存在长期而普遍的稀缺性，模型工作者一定会需要对模型结构和参数进行插值和调整。而对这一过程的校准常常通过试错法来进行，也就是反复调整模型中指定的输入参数，直到水头的模拟值与真实值达到一致。在数据和参数量较大的情况下，也可以使用专门为参数识别而设计的计算程序来实现自动反演。

完成一个经过校准的水流模型后，就可以着手构建溶质运移模型。与水流模型类似，在溶质运移模型里必须定义初始条件、边界条件、源汇项、以及控制溶质运移反应的关键参数。最为常见的溶质运移模型是基于物理过程（对流和弥散）的ADE模型（advective-dispersive model）。对流是指溶解态溶质在含水层中以地下水平均实际流速（即平均流速）迁移的现象。弥散描述的是溶质烟羽和周围地下水的混合过程，这一混合过程在很大程度上是由实际流速和平均流速的差异引起的。造成这些流速差异的来源很多，例如分子扩散、流线的扭曲、以及介质的非均质性。弥散作用的结果通常是溶质烟羽的扩大和边界的模糊。

举一个最简单的可溶性溶质的运移作为例子，考虑氚在地下水中的衰变和运移时，最少要指定三个参数：纵向弥散度、横向弥散度、以及有效孔隙度。两个弥散度参数有时也被称为含水层的“特征混合长度”，即溶质在复杂的非均质介质中运行时表现出的混合程度；而有效孔隙度与溶质运行的速度有关，因为溶质的实际运行速度是达西流速除以有效孔隙度。这三个参数在空间上都会变异，但一般而言，数据的稀缺性导致大多数模型中整个含水层使用一致的取值。

溶质与地质介质间的反应关系常常难以界定，最常见也是最常出现的错误就是把溶质简单地认为是保守的，既不与固相介质也不与其它液相组分发生反应，而且本身也不会发生任何形式的衰变。略微进步一些的是线性吸附模型，即溶质在液相和固相表面始终以平衡态维持线性的质量分配关系。过于简化的反应关系往往会为模拟带来许多问题，这可以通过细致的实验研究来解决，事实上全球范围内已经有很多成功将溶质运移模拟应用在复杂水岩作用中的先例。

溶质运移的模拟通常基于ADE模型的解算。如果忽略化学反应而仅考虑单一溶解因子，有如下的控制方程：方程左边是溶质浓度随时间的变化，方程右边分别是弥散项和对流项对溶质浓度场造成的影响。

其中： C为溶质的浓度，（ML-3）；t为时间，（T）；Xij为在直角坐标系下沿各方向上的距离，（L）；Dij为水动力弥散张量，（L2T-1）；Vi为孔隙水流实际速度，（LM-1）。

控制方程的解算

有时因为可用的数据过少，研究人员没有其它选择，只能假定含水层是各向同性的均质含水层，而且具有单向稳定的流场。这时可以对上述方程求取解析解，对溶质运移的情形进行非常粗略的估计。事实上，在一些情形下，这样级别的估算已经足够。然而大多数情况下，这样的处理是不够的，因为地下含水介质是地质年代间由各种沉积和构造活动形成的复杂地质体，不均一性和各向异性是常态。针对复杂水文地质问题的解析解是不存在的，此时唯一可用的办法是对控制方程求取数值解。

在求取ADE方程数值解的工具中，最常见的有三类：有限差分法（FD）、有限单元法（FE）、和特征线法（MOC）。这三类方法各有优缺点，互为对方的有机组成部分，其典型的剖分方式如下图所示。

有限差分法是最简单的方法，也是最常见的求取数值解的方法。在此方法中，微分方程被近似替换为相应的差分方程，而求解空间也从全维度空间压缩为一系列有限的单元格。对于每一个单元格，仅有一个水头值和浓度值。在有限差分法中，不仅空间被网格化离散，时间维度也被离散化，成为一个一个的时间步。这样，溶质浓度在时空中的变化就可以表达成一系列的代数方程，用通行的求解方法进行求解。

有限差分方法的优势在于它的概念非常直观，也便于使用编程的形式用计算机实现。有限差分方法的局限性之一在于，若想对客观世界做精细的模拟，则必须产生非常细碎的网格。虽然计算能力的飞速发展已经渐渐使得这一局限性变得不甚重要，但也不能否认，考虑溶质运移的有限差分模型，尤其是在运行自动调参或不确定性分析时，还是会消耗大量的计算资源。此方法的另一个局限性可能更为重要，即其理论局限性。当弥散项远小于对流项时，控制方程近似为一个双曲型差分方程（用来描述冲击波的一种数学物理方程）。有限差分法对于抛物线方程非常有效，而对双曲线方程则不稳定，会出现数值震荡和数值弥散。需要指出的是，单纯就弥散效果而言，数值弥散和真实弥散完全相同，只不过数值弥散是控制方程在解算过程中被动引入的弥散，无法人为控制。

为了部分解决有限差分方法中关于边界的几何形状的限制，地下水科学家在80年代引入了有限元方法。在这种方法中，模型区域被拆分为不规则的三角形单元，之后对每一个单元求取溶质浓度值。有限差分方法用近似的差分方程取代了原始的微分方程，然而有限元方法并不对控制方程进行近似，仅对其解进行近似。有限元方法中使用形函数（shape function）的概念，用形函数的线性组合来表征溶质浓度分布。伽辽金方法是最常见的一类有限元方法，它的基本原理是通过迭代将每个单元的残差均值降到最低，从而在整个有限单元网络内求取最佳近似解。

从上述简介可以看出，有限元方法在便于理解、便于数学表达、便于编程等方面不及有限差分方法，从而也在代码修改和传播中处于劣势。反观其优势，核心点在于当有限差分无法正常应用时，有限元方法有时可以取得令人满意的结果。

特征线法与前两类方法存在本质不同，因其使用了拉格朗日体系而非欧拉体系对控制方程进行求解。拉格朗日描述与欧拉描述乃描述流体运动的两种体系。前者是以移动中的单个质点为参照体系，而后者使用了固定的参照体系。在拉格朗日描述中，流体的任何物理量不仅是时间t的函数，也是初始坐标矢量R的函数，所以这是一种以流体质点初始位置R和时间t描述流体运动的方法。特征线法常见的应用方式是向流场中引入一系列随水流运动的质点，这些质点在运行中逐渐扩散，对这些质点的追踪可以提供对浓度变化的直观印象。这一方法的优点是非常直观，而缺点也很明显，即无法处理多组分和化学反应问题。

溶质运移中的化学反应

模拟溶质的运移远远复杂于单纯模拟地下水的流动，原因不仅包括上述的数值问题，还在于溶质之间以及溶质与介质之间会发生反应。对这一过程进行的模拟可能涉及到多种溶质，也就是说会需要解算多个包含反应项的偏微分方程。溶质的反应可以分为两大类：一类由反应速率控制，一类由化学平衡控制。平衡控制的反应很容易理解，此类反应是可逆的，并且可以在瞬间（理论上）达到平衡点。速率控制的反应在模拟中的主要难点是，人们很少能够知道此速率的真实取值，即便在实验室中也是如此。有两类反应必须在反应速率控制下的体系内进行模拟：一类是较为缓慢的反应（相对于对流和弥散而言）；另外一类是不可逆的反应。

化学反应的另外一个范畴是反应仅在液相中进行，还是会涉及固相。仅在同一相中进行的反应被称为均一性（Homogeneous）反应；而那些涉及在不同相中传质的反应称为异质性（Heterogeneous）反应，包括结晶、溶解、和类似于吸附和离子交换的表面反应。客观世界中发生的反应类型决定了模型中所使用的概念范畴。下面举一些具体实例来说明不同概念模型的数学表达。

在ADE模型中加入反应速率项是较为直观的操作。考虑最为简单的化合反应，组分P1和P2结合形成P3，三种组分的浓度分别为C1，C2，C3。

其中Kf和Kb分别是正向反应和逆向反应的速率常数。此时，三种组分的浓度值可分别表示为：

此时的控制方程变为以C1，C2，C3为未知数的多项式微分方程，在空间对其进行有限差或有限元离散后，可以使用非线性矩阵计算方法（比如预估-校正法或高斯-牛顿法）求解。

讨论

以上的论述说明了化学反应可以通过特定的数学手段与ADE方程联立起来，最终求取化学组分浓度在时空上的函数分布。控制方程表达方式要根据具体反应类型来确定，对于有些简单的反应类型（如吸附-解吸）单个偏微分方程即可作为控制方程。当考虑的反应类型变得复杂后，偏微分方程的数量就会增多，模型的非线性特征也会增强。

溶质运移模拟需要深入了解含水层的水力特征、运移参数、化学反应类型、化学参数、以及精确的污染现状。如果这些信息并不齐全，研究人员必须仔细考察模型的适用性。一个数据非常稀疏的模型如果用于污染预测，其可信度会大打折扣，这时模型的作用主要是帮助研究人员理解含水层的行为和特征，并指出数据缺口。一个完全没有实测数据的模型，其作用不会大于完全基于经验的猜测。