第三章 化学元素在地下水中的迁移

第三章  化学元素在地下水中的迁移 地下水中的溶解组分一方面可以随着地下水的运动而运动，同时它们也可以在自身浓度梯度的作用下进行迁移，研究化学元素在地下水中迁移规律的理论通常被称为水动力弥散理论。通过对地下水中所含溶质运移规律的研究，可对其在空间和时间上的分布与变化进行预测和控制，这就是水动力弥散问题研究的目的。 影响溶质运移的主要因素是：对流，机械弥散，分子扩散，固相—溶质间的相互作用（如溶解、吸附等），溶液内部的化学反应，溶质其他的源、汇作用（如放射性元素的衰减，作物根系对某些溶质的吸收等）。为达到研究弥散问题之目的，这些因素都是需要设法予以定量考虑和描述的。 3.1  水动力弥散现象及机理 3.1.1  水动力弥散现象 首先，我们来看几个实例。 【例1】 设有一地下水均匀流场，其宏观平均流速为 （真实速度），在某一过水断面上，同时开始均匀持续地注入某种示踪剂（例如NaCl溶液），观察示踪剂随水运移的情况（参见图3-1-1），会发现下游出现示踪剂的前缘部位，超前于按地下水宏观平均渗透速度 计算所预期能到达的部位。在与示踪剂推进锋面相邻的上游侧，存在一个过渡带，其中示踪剂的浓度由c0（地下水中原始的示踪剂浓度）向上游侧逐渐过度到c1。 图3-1-1  均匀流场中面状持续性注入示踪剂产生的水动力弥散现象 【例2】 在地下水均匀流场中的某点上瞬时注入示踪剂，并立即停止注入，观察下游不同时刻示踪剂浓度的分布状况。发现示踪剂所占据的区域逐渐向下游移动，并沿纵（水流方向）横（垂直水流方向）两个方向逐步扩展开来，范围愈来愈大。按宏观平均渗流速度计算的示踪剂在某时刻预期到达的位置，是该时刻示踪剂实际的似椭园形分布范围的形心。在形心处，示踪剂浓度最高，向四周逐渐变小。时间愈长，形心愈远离示踪剂注入点，其最大浓度值相对愈小，而示踪剂的似椭圆形分布范围则愈来愈大。示踪剂锋面出现的部位也超前于按宏观平均渗透速度所预计的部位，参见图（3-1-2）。 图3-1-2  均匀流场中点状瞬时注入示踪剂产生的水动力弥散现象 在该例条件下，如果在一点处示踪剂的注入不是脉冲式（瞬时注入、立刻终止）的，而是持续性的，则示踪剂将以似梨的形状逐渐扩展开来；同一时刻，距注入点不同距离处的横向浓度分布如图（3-1-3）所示。纵向浓度分布从注入点起是愈向下游愈小的。示踪剂锋面仍超前于按宏观平均渗流速度所预计的部位。 图3-1-3  均匀流场中点状持续性注入示踪剂产生的水动力弥散现象 【例3】 用含有示踪剂（浓度为c1）的水连续地驱替长度为L的均质园形砂柱中原来（示踪剂浓度为 c0）作稳定流动的水。示踪剂的浓度c1较低，它对流体密度产生的影响很小，从而对水流运动产生的影响可以忽略不计。在砂柱末端取水样测定示踪剂浓度随时间的变化 ，将结果绘成相对浓度  ( )与时间的关系曲线，参见图（3-1-4）。可见， ～t曲线呈s形，而不是按平均渗透速度所预期的阶梯形（虚线）。示踪剂锋面出现的仍超前于按宏观平均渗流速度所预期的部位。 以上各例中所表现出的示踪剂在流场中运移的共同特点是：示踪剂在注入地下水后，它就在流场中逐渐传播、扩展，占据的区域愈来愈大，超出了按宏观平均流动所预期的范围。这种现象就称为流体动力弥散。它是一种不稳定的、不可逆的过程。所谓不稳定，是指示踪剂 浓度随时间常常是变化的。不可逆是指弥散现象发生到一定程度时若改变渗流方向(即作逆转流动)，示踪剂不可能恢复到原来的初始分布状态。 3.1.2 弥散现象的机理与几个相关的概念 概括地说，弥散现象主要是由于两方面的因素综合作用而形成的。一是液体体系的微观运动速度(V)与宏观运动速度（ ）的不一至性；二是所考虑的α组分（某种示踪剂）的微观运动速度（Va）与流体体系的微观运动速度（V）的不一致性。这两个不一致性最终表现为 与 的不一致性。因此，可以更简捷地说，流体动力弥散现象是由于 与 的不一致性所引起的。下面对此作以具体的剖析。 关于V与 的不一致性，可以从以下几方面来理解。首先，我们知道 是以V为基础通过在表征体元（REV）上积分而得到的平均速度。显然，这样的 不论在大小上与方向上一般都不同于中各点水流的实际（微观）流速V，当然一般也不同于REV中心点处水流的实际流速V。有的地方V与 间的偏差可能是很大的，甚至是反方向的，参见示意图（3-1-5a）。 图3-1-5  水动力弥散现象的形成机制 其次，由于实际存在于孔隙空间中的水流是示踪剂的携带者和输运者，因此，孔隙水流的实际偏宏观平均流线必然要引起示踪剂分布范围偏离按宏观平均流动考虑所预期的范围。例如，按平均流动考虑，在稳定流条件下从上游投入图（3-1-5b）中所示流管内的示踪剂 ，不可能运移到该流管之外去。然而由于 与V的不一致性，示踪剂的分布范围会大大超出流管范围。例如，象 与VA那样的不一致性，就会使示踪剂运移出流管之外，发生横向弥散。在垂直流管方向上，若示踪剂的浓度梯度较大，则会加强其横向弥散性。 第三，在某些点上， 与V可能具有相同或大致相同的方向，但它们的大小却可能相差很大。例如在颗粒之间的窄长通中，即使 与V的方向大致相同，但由于流体自身的内摩擦力和流体与颗粒间的摩擦力相比要小得多，就可能使通中心部位的V的长度较 的长度大得多，参见图（3-1-5c）。在这种情况下示踪剂出现部位自然就会大超前于按 所预期的部位，发生了纵向弥散。若在 方向上示踪刘的浓度梯度较大，且是指向 反方向的，则会加强纵向弥散。 由以上可以看到，纵向弥散主要是当V与 的方向大体一致时由于 而引起的；横向弥散主要是当由于V方向大偏离 方向而引起的。纵向弥散与横向弥散通常总称为机械弥散或对流扩散。 关于V与Vα的不一致性，其成因可以从两方面来分析和认识。一方面，α组分（示踪剂）的浓度梯度力图驱使示踪剂沿着▽c的反方向运移；而▽c矢量与V矢量一般是不共线的，故▽c的存在就成为造成Vα偏离V的主要因素；另一方面，流体体系的运动速度V又力图驱使示踪剂随流体体系运动的总体方向运移。这两种影响示踪剂运移方向的因素综合作用的结果，一般使Vα即不取▽c的反方向，也不取V的方向，而使Vα与V保持某种程度的偏差。当V相对于▽c显得很小时，▽c对Vα的影响就显得很突出，以致处于控制地位，这常会形成Vα与V间很大的偏差；当V相对于▽c显得很大时，V对Vα影响就显得很突出，以致处于控制地位，这常会形成Vα与V接近一致。例如在本节的【例2】中，当示踪剂注入后，注入点示踪剂的浓度自然就大大高于附近地下水中这种示踪剂 和浓度（本底值），这就形成了浓度梯度（▽c），其方向显然具有向四周辐射的特点，它力图驱使示踪剂向注入点的四面八方扩散开来。而注入点周围V的方向一般来说是不具有向四周辐射之特点的，这就出现了▽c与V方向的不一致性，并导致Vα与V方向的偏差。                Vα与V的不一致性也导致Vα流线与V流线的不一致性，引起示踪剂的另一种弥散，通常称之为溶液中的分子扩散。由于示踪剂浓度在空间上分布的不均匀性是引起Vα与V偏离的主要原因，所以通常把由于示踪剂浓度分布的不均匀性而引起的弥散称为分子扩散。溶液中的分子扩散在多孔介质REV上平均所得和宏观表现，就是多孔介质中的分子扩散。 显然，分子扩散与机械弥散是同时发生并互相影响的。因为，一方面示踪剂浓度梯度的存在影响着Vα的大小和方向，而Vα的大小和方向又影响着V及 的大小和方向，影响着V与 的偏离程度；另一方面， 及V的大小和方向又影响着Vα的大小和方向，影响着Vα与V的偏离程度。机械弥散中蕴含着分子扩散的作用；分子扩散中也蕴含着机械弥散的作用。两者是一个有机的整体，很难将它们分割开来。机械弥散与分子扩散之分，仅在于各自的主导因素不同而已。 机械弥散与分子扩散，总称为流体（水）动力弥散。综上所述显然可以简捷地说，水动力弥散是由于 与 的不一致性所引起的。弥散现象的发生与流态无关。不论渗流是层流或是紊流，都会发生水动力弥散。 3.2  弥散通量与弥散系数 3.2.1  弥散通量 由于弥散作用所引起的单位时间通过单位溶液面积的溶质质量，称为弥散质量通量，简称为弥散通量。弥散通量可分为分子扩散通量与机械弥散通量两种。它们分别是由于分子扩散作用与机械弥散作用所引起的质量通量。 通过实验和描述多孔介质中溶质运移问题的简化（理想化）模型的研究，认为可以采用计算溶液中分子扩散通量的费克定律那种结构形式的公式，来计算多孔介质中的分子扩散与机械弥散通量。不过这时算式中系数的性质已发生了很大的变化，它一般不再是标量，而是二阶张量了。 1．  多孔介质中的分子扩散通量（J″）：可由下式计算： （3-2-1） 式中：J″— 由于分子扩散而引起的单位时间通过单位孔隙面积的溶质质量， 。这里假定含水层是被水饱和的，孔隙面积等于溶液面积。 － 多孔介质中的分子扩散系数，量纲为 ，在各向异性介质中为二阶张量，在各向异性介质中为标量。 c － 溶质（示踪剂）的浓度，即单位体积溶液中所含溶质的质量，量纲为 。 － 溶质浓度(c)的梯度，其各坐标的量纲为 。 2．多孔介质中的机械散弥散通量（J′）：可用下式计算： （3-2-2） 式中：J′－ 由于机械弥散作用而引起的单位时间通过单位孔隙面积的溶质质量。这里仍假定含水层是饱和的。 － 机械弥散系数，是二阶张量。 c、▽c的含义同上。 由式（3-2-1）、（3-2-2），多孔介质中的水动力弥散通量（J）可按下式计算： （3-2-3） 其中： － 称为多孔介质中的水动力弥散系数，是二阶张量。 3.2.2  弥散系数 1. 弥散系数的主轴和主值 弥散系数与各向异性介质中的渗透系数都是二阶张量，它们的物理意义虽然是完全不同的，但作为二阶张量，它们的某些性质则是完全类同的。渗透系数 有主轴与主值的概念，弥散系数 也有主轴与主值的概念，且它们的意义是完全类同的。在正交直角坐标系中，能使 成为对角型 （3-2-4） 的三个坐标轴称为弥散系数 的主轴，相应和三个方向称为弥散的主方向。对角型 中的三个非零元素 、 与 称为弥散系数的主值或主弥散系数。同时常把 称为纵向弥散系数，把 、 称为横向弥散系数。 Bachmat和，贝尔（Bear）证明了：在各向同性介质中， 的主轴中有一个轴处处都与地下水宏观平均渗流速度（V）方向一致（这里将 也简记为V，下同）。这样，其余两个主轴自然就处处都与V方向垂直。 当取坐标方向与弥散主轴方向一致时，水动力弥散通量J在三个坐标轴上的投影 、 、 可表为 （3-2-5） 2. 弥散系数的确定 （1）多孔介质中的分子扩散系数（ ）：该系数取决于溶液中的分子扩散系数（Dd）和多孔介质的弯曲率（ ），并按下式计算 （3-2-6） 其中 在各向异性介质中为二阶对称张量，在各向异性介质中退化为标量，无量纲。系数 是标量，其值的大小与溶液种类有关。对同种溶质，其 值还与溶液温度及溶质浓度有关。 经Bear分析，认为对于未固结的各向同性介质，可取2/3作为 的估计值。沙夫曼（Saffman）建仪用1/3作为 的估计值。T″＜1骨架颗粒的存在对分子扩散起了阻滞作用。在流体处于静止状态（V＝0）的各向同性介质中作弥散实验，可求得分子扩散系数D″（此时为标量，故不写作 ），再由（3-2-6）式即可求得弯曲率T″ 2．机械弥散系数（ ）：许多学者研究表明，机械弥散和系数（ ）取决于宏观渗透速度V、彼克来特（Peclet）数（Pe）和介质特性。Pe数是一个无量纲量，按下式计算： Pe＝                           （3-2-7） 式中V=︱V︱，d是介质颗粒的平均粒径。介质特性通过介质的（几何）弥散度（率）αijkm来表征，αijkm是四阶张量。 贝尔采用毛细管网络的简化模型，对弥散问题进行了定量研究，导出了 与相关量之间的下列函数关系（其中使用了爱因斯坦求和约规，k、m都是求和指标，k、m＝1，2，3） （3-2-8） 式中： （无量纲） δ－多孔介质单个通道的特征长度（L）与其横断面的水力半径（ ）之比 ，无量纲。 Vk、Vm－V在k、m坐标轴上的投影。 弥散度αijkm是介质骨架几何特性的物理量，具有长度量纲[L]。在各向同性介质中，αijkm是一个四阶对称张量，且其全部分量（元素）都可用两个参数αL（多孔介质的纵向弥散度）与αT（多孔介质的横向弥散度）的下列线性组合表出 （3-2-9） 其中δ ij是克罗拉格（Kronecker）δ函数，其定义为 （3-2-10） 当i、j、k、m分别取1、2、3（x、y、z）时，αijkm共计有81 个分量，其中只有36个非零元素，大量元素为零。 在Pe的表达式（3-2-7）中含有溶液中的分子扩散系数Dd，这反映出机械弥散与分子扩散的不可分性。当V很大因而Pe相当大时，有 ，这表明此时分子扩散对机械弥散的影响已很小。在法国，Pe数的经验范围为1＜Pe＜100。 将αijkm的表达式代入（3-2-8）式中，并取 ，就得到各向同性介质中计算机械弥散系数各分量的常用公式： （3-2-11） 因为此时 上式右端第一项 第二项 一、二两项结果相加就得到（3-2-11）式。 由（3-2-11）式可见，在各向同性介质中机械弥散系数 是二阶对称张量，在笛卡尔坐标系下，其六个独立分量为： （3-2-12） 对于{xoy}平面内的二维流，有 ，这时 简化为： （3-2-13） 即 上式表明，在{xoy}平面内的二维渗流中，虽然 ，但只要z方向浓度分布不均匀，致使 ，则在z 方向照样有机械弥散通量存在，其值为 。由此我们已看到，渗流维数与弥散维数常常是不一致的。二维渗流中可以发生三维弥散，一维渗流中也可以发生二维或三维弥散。 对于单向渗流，取流动方向为x方向是很方便的，这时 ，此时对于流动区域内任一点，处处都有 （3-2-14） 即 其中 、 分别称为纵向与横向机械弥散系数。由上式可见：（1）即使在单向渗流条件下，从某点注入示踪剂（或污染物由某点进入地下水）后，由于一般来说， 、 及 都不等于零，因而除了纵向（水流方向）弥散外，横向弥散也是存在的。（2）即使 ，由于 ，因而纵向弥散通量与横向弥散通量也是不相等的，这显示出弥散现象是各向异性的。（3） 、 都是V的线性函数。 显然，当 时，这里关于 、 是V的线性函数的结论就不再成立。许多实验和一些分析研究表明，在一般情况下， 、 具有下列形式： （3-2-15） 其中m1、m2为常数。 在各向异性介质中，对弥散度 研究得还不甚成熟，有的认为它仍是一个对称张量，有的则持怀疑态度。在非饱和区的弥散问题中， 、 、 及常数m1、m2都是含水率θ的函数，且一般要通过实验才能得出函数关系。 3.2.3  通过一维弥散实验综合成果分析弥散系数 上段中介绍的确定弥散系数的公式，是通过对理想化模型作数学物理分析得到的。它们对于实际情况符合程度如何，自然是人们关心的问题。 本段通过一维弥散实验的一个综合成果，来分析一下弥散系数。从中可看到，上段中关于V相当大（以致 ）时 与V 呈线性关系的结论大体上是符合实验结果的。 实验通常是这样进行的：在一长直圆筒中装满未固结的多孔介质（例如砂）。砂柱两端保持恒定的水头差，使其中的渗流达到承压稳定状态。测定这时渗透水流中所含某种溶质（例如Cl－）的原始浓度。然后在砂柱高水位一端突然改换为该种溶质浓度较高的同种溶液持续进水，去驱替砂柱中原先的水流。在砂柱的下游（低水位）一端不断地测定渗出液中所含溶质的浓度，直到与驱替液中溶质的浓度相等为止。在同一实验中，上下游水头差保持不变。调整不同的水头差进行类似的实验，就可以得到不同渗流速度下的多次实验。利用这些资料和相应弥散问题的解析解，反求出各次实验下的纵向水动力弥散系数（DL），可绘出 ～ 关系曲线。据许多学者的实验所得到的综合示意曲线如图（3-2-1）所示。 该图可分为五个区： Ⅰ区：V很小，从而Pe也很小， 。这表明此区内V的大小变化并不引起纵向水动力学弥散系数 的变化（这里 ）。由于机械弥散系数 是与V有关的，故上述情况说明在低速流动情况下 ，或 ，即在低速流动条件下的水动力弥散中分子扩散是起主导作用的，机械弥散是微不足道的。 Ⅱ区：流速较大，Pe数大致在0.4～5之间。曲线开始向上明显弯起，表明V的增大已使 明显增大， 的情况已不再成立， 的两个组成部分（ 与 ）那一个都不能忽略，机械弥散与分子扩散的作用大致相当。 Ⅲ区：V继续增大，机械弥散的作用愈显著，占据主要地位，同时也有一定的横向分子扩散，它们互相干扰，后者对前者往往起到削弱的作用。曲线以近乎直线的趋势上升，倾角略大于45°，实验得出 （3-2-16） 这表明 与V及 与V已近似呈直线关系。 Ⅳ区：V更大，但没有达到达西定律不适用的程度。机械弥散占居主导地位，分子扩散的作用可以忽略不计， 。曲线继续直线式上升，倾角约为45°，试验得出 （3-2-17） 这表明 与V及 与V呈直线关系。 Ⅴ区──V过大，这时机械弥散居于统治地位，但紊流和惯性力对弥散的影响已不可忽略，它们对机械弥散起到削弱的作用。曲线上升趋势逐渐缓和，倾角小于45°。 由上述实验综合结果可见：在V不过分大也不过分小的条件下（即在Ⅲ、Ⅳ两区中）， 与V及 与V基本上是呈直线关系的。这证明了前面关于 与V呈直线关系的理论基本上是正确的。 关于横向弥散系数 ，一般仍采用 这样的关系式。对于第Ⅲ区，常取 。注意到在该区中 这个前面已给出的结果，就得 即纵向弥散系数总大于横向弥散系数。一般文献中指出， 在5～24之间。关于 的实验资料还较少见。 3.3  水动力弥散方程 水动力弥散方程是一种习惯用语，其实质是水中某种溶质的质量守恒方程。它在溶质输运问题中起着基本的、控制性的作用。本节将导出该方程，并对源汇项予以讨论。 3.3.1  弥散方程的推导 过渗流区内任一点作一个长、宽、高分别为Δx、Δy、Δz的微元正六面体（参见图3-3-3），研究该控制体内溶质的质量守恒关系，就可得到微分形式的溶质质量守恒方程，即弥散方程。 按溶质质量守恒原理，当溶质无源无汇时，控制体内溶质质量的变化必等于通过控制面流入与流出控制体的溶质质量的代数和，或者说必等于通过控制面净流入控制体内的溶质质量。 由控制面流入或流出控制体的溶质质量有两类，一类是对流质量通量（cV），另一类是水动力弥散质量通量（J）。下面分别计算任一瞬时（t）由这两类质量通量所引起的控制体内溶质质量的净增量，建立质量守恒关系。 1. 由于水动力弥散所引起的控制体内溶质质量的增量 沿x方向通过控制体前后两个控制面流入与流出控制体的溶质质量及相应的质量增量分别为： 从后面流入      从前面流出      净流入（增量）  流出量的算式中略去了 等高阶无穷小量。这对最后结果并不带来什么误差。 沿 方向通过控制体左右（上下）两个控制面流入与流出控制体的溶质质量及相应的质量增量分别为： 左面            右面            增量            下面            上面            增量            综合以上结果，便得到由于水动力弥散作用所引起的控制体内溶质质量增量为： 2. 由于对流作用所引起的控制体内溶质质量的增量 作与上面类似的计算，易得由于x、y、z三个方向的对流作用所引起的控制体内溶质质量增量分别为： x方向        y方向        z方向        故由对流作用所引起的控制体内溶质质量增量为： 由于对流与水动力弥散的共同作用所引起的同一瞬时（t）控制体内溶质质量的变化显然为： 按质量守恒关系，该值必等于 ，约去 即得溶质无源无汇时的弥散方程为： （3-3-1） 或简写为 （3-3-2） 当研究区内部有溶质的源或/和汇时，设单位体积含水层单位时间所产生（源）或吸收（汇）的α组分溶质的质量为Wα，则在上式中补充上源汇项便可得到： （3-3-3） 这就是多孔介质中溶质的质量守恒方程，常称为水动力弥散方程，或较确切地称为对流-弥散方程。 假定n＝C（这意味着假定多孔介质骨架是均质且不可压缩的）， ，则（3-3-3）式可简化为： （3-3-4） 其中： （3-3-5） 当流体不可压缩时，有 ，从而（3-3-5）式可进一步简化为： （3-3-6） 这就是常见、常用的水动力弥散方程。将其右端两项展开，就得到： （3-3-7） 对上式采用爱因斯坦求和约规，可简化为 （3-3-8） 对于二维流（例如 的渗流）三维弥散的情况，（3-3-7）式可简化为 （3-3-8） 在所述的二维流中如果弥散也是二维（ ）的，则上式中的 ，其余不变。 在上述的二维渗流二维弥散中，常把弥散方程写成另一种形式。假定含水层基本呈水平分布，不可压缩，但含水层厚度（M）与孔隙率（n）在水平面内的变化不宜视为常数。这时可将溶质质量守恒方程扩大到一个具有单位水平面积的含水层铅直柱体中，得到如下的弥散方程： （3-3-9） 3.3.2  关于源汇项 除对流与水动力弥散作用外，凡存在于研究区域内部，能引起微元六面体内部某种溶质质量变化的其他一切因素，都称为源汇因素，均需要补充到弥散方程中去。这种补充项（例如3-3-3式中的Wα项）通常统称为源汇项。 源汇因素是相当多的，主要的有： （1）含水层颗粒对溶质的吸附； （2）溶液中溶质的沉淀； （3）固相颗粒表面与溶液间的离子交换； （4）固相颗粒中某些物质的溶解； （5）溶液中发生的化学变化； （6）放射性元素的衰变； （7）作物根系对溶质的吸收； 一般说来，源汇因素的具体内容不同，源汇项的具体表达式也就不同。下面给出几个例子。 1．放射性元素的衰变 若所研究的溶质是放射性元素，那么，即使没有对流与水动力弥散作用，其浓度c也在不断地变小的，其衰变规律为： （3-3-10） 其中：λ为常数，称为衰减常数；c0是初始时刻放射性元素的浓度。 据上式易得，由于衰变而引起的溶质浓度随时间的变化率为： （3-3-11） 显然，如果再没有其他源汇因素，则放射性元素浓度的时变率就是由对流、弥散与衰变这三方面因素共同引起的，因而相应的弥散方程为 （3-3-12） 由上式可见，即使等号右端第一、二两大项全为零， 也并不为零，而为 。按放射性元素的特性，正应如此。 若溶液中发生的某种化学反应使所研究的那种溶质的浓度减小，而减小的速度与该种溶质的浓度成正比，那么也可用 这个表达式作为源汇项。不过在这里 的意义已不是衰减常数，而是比例系数。 2．溶解与吸附 含水层骨架颗粒中某些物质成分的溶解，可能会增加溶液中所研究的溶质的浓度；溶液中所研究的某种溶质，也可能被颗粒吸附而减小其浓度。显然它们都属于源汇因素。这里分析一下这类源汇项的具体表达式。 设：F表示固相表面的溶质浓度，即固相表面单位体积中所含的我们所研究的那种溶质的质量； 表示单位时间从单位体积含水层的固相（体积为1－n）进入液相（体积为n）的溶质质量（固相溶解时是正的进入， ；固相吸附时是负的进入， ）。则对于单位体积溶液来说，单位时间内溶解或吸附所增加或减少的溶质质量就为： （3-3-13） 与此同时，对单位体积固相来说，单位时间内因溶解或吸附所减少或增加的溶质质量（即 ）就为： （3-3-14） 由式（3-3-14）便有： （3-3-15） 将上列源汇项的表达式加入到弥散方程右端，就得： （3-3-16） 即： 假定介质骨架不可压缩（n不随t变化），上式可写为： （3-3-17） 这就是考虑到溶解与吸附时的弥散方程。 若所研究的溶质是放射性元素，则在溶解与吸附作用的同时，还存在着它自身的衰减。这时的弥散方程显然应为 （3-3-18） 在（3-3-18、19）式中，c和F都是未知函数。这就要求有一个F～c间的关系式，才能与上列弥散方程一起组成封闭的方程组。若采用下列的线性吸附关系： （ ，为常数） 即假定F与溶液中溶质的浓度c成正比，则将其代入（3-3-18）式中就得 令： （3-3-19） 就有： （3-3-20） 由于 ，上式表明吸附所产生的影响等价于渗流速度和弥散系数的减小（ 减为 ， 减为 ）。这将使对流质量通量与弥散质量通量减小，故常将 称为阻滞因子。 3．越流补给和抽水消耗 在{xoy}平面内的二维渗流、二维弥散问题中，若含水层得到了越流补给，且含水层中有N个抽水井作用时，设越流补给强度为W，越流水中所含溶质的浓度为cW，第i口井的抽水流量为Qi，其中所含溶质的浓度为cQi，则可将溶质的这种补给与消耗作为源汇项加到（3-3-9）式中去。这时应采用的弥散方程就成为 （3-3-21） 其中 为狄拉克（Dirac） 函数。当采用上面这个方程研究问题时，水井就不能再作为边界对待了，这是需要注意的。 这里需要指出，除了较简单的少数几种源汇因素外，欲得到大多数源汇因素的源汇项的具体表达式常常是相当困难的事。因为它们有的需借助于总结实验资料，有的需借助于理论分析，更多的是需要借助于实验与理论的某种结合，才能获得较可靠的表达式。在实际问题中源汇因素总是存在的，而且常常是几种因素并存且相互影响的。因此，在溶质的对流与弥散规律已有了一定研究成果可用的今天，对源汇项的深入研究，实际上已成为能否较精确地预测溶质时空分布规律的关键和难点。对此应有清楚的认识，给予高度的重视。欲使多孔介质中溶质运移理论得到发展，走向成熟，这正是一个重要的生长点。 3.4  定解问题的类别与定解条件 3.4.1  溶质运移定解问题的类别 1. 溶质浓度（c）的变化对溶液密度（ρ）的影响较大，ρ不能作为常数处理的问题 研究多孔介质中溶质运移问题的最终目标，通常是求出人们所关心的某种溶质浓度的时空变化规律，从而为各种具体目的服务。溶质运移问题通常也简称为弥散问题（虽然这种说法是不准确的，但人们仍习惯于这种说法）。 从弥散方程来看，其中除了c之外还包含着Vx、Vy及Vz。这表明，如果不知道V的时空变化规律，则溶质浓度c的时空变化规律是无法确定的。 然而，c的变化会引起流体体系（溶液）密度ρ和粘滞系数μ的变化，因而也就会引起与之相关的渗流速度V、压强p等物理量的变化。而ρ、V、p等因素的变化反过来又会影响 c的变化（严格说来也会影响μ的变化，不过在等温流体中μ的变化很小，可以不予以考虑）。这就是说，c，ρ，V，p之间是互相联系、互相制约的，而不是互相独立的。这样，欲求得 就需要解一个以c， ρ，V，P为未知变量的方程组才能实现。在这种情况下，c，ρ，V，P的时空变化规律是同时联立求出的，这类问题简称为联合（或耦合）模型问题。 2. 溶质浓度（c）的变化对溶液密度（ρ）的影响不大，ρ可以作为常数处理的问题 例如溶质浓度相当小时，它的变化对ρ的影响就可忽略，ρ就可以作为常数处理。这时就可以把求 的问题与求 的问题作为两个彼此独立的问题分别去解决。其次序是先求得V，将其代入弥散方程中再求得c。而求V的问题这时已成为一个求水头H的问题；求得 后用达西定律即可求得V（若需知道 ，也易由 求得）。由此可见，这类问题最后被分解为两个彼此独立的问题：一个是求 的水流问题；另一个是求 的溶质运移问题。因此，这类问题可简称为独立模型问题。 在联合模型问题中，要给出所有因变量（ c， ρ，V，P）各自的定解条件才能求解。在独立模型问题中，只需给出关于水头H和溶质浓度c的定解条件即可。显然，独立模型较联合模型要简单得多。再要求精度不很高的情况下，都可采用独立模型。 关于水流问题的定解条件在第三章中已有详述。下面只讲在两类定解问题中都要涉及的关于c的定解条件。 3.4.2  关于溶质条件（c）的定解条件 关于c的定解条件仍包括初始条件与边界条件两类。若只需研究溶质运移最终所达到的稳定状态，则定解条件只需要边界条件，不需要初始条件。 1. 初始条件 在所选定的初始时刻（t=0或t=t0）研究区内某种溶质浓度的空间分布状况，即 称为该种溶质浓度的初始条件。 例如，t=0时在无界含水层中某点（x0，y0，z0）上瞬时注入质量为M的示踪剂，欲研究t＞0后这种示踪剂在含水层中的分布及随时间变化的情况，其初始条件可表为： （3-4-1） 这里假定质量M是注入到以（x0，y0，z0）为中心的单位体积含水层中的。单位体积饱水含水层中所含的溶液体积为n，故对注入点处单位体积溶液而言其所含示踪剂的浓度为（M/n）。乘以狄拉克（Dirac） 函数就表示除注入点外其他地方 。这种注入方式称为栓子式或脉冲式注入。 在一般情况下，初始条件需根据实测溶质浓度分布资料才能给出。 2. 边界条件 （1）第一类边界条件 已知（或给定）研究区域边界上溶质浓度分布状况的边界条件，称为第一边界条件。 例如，对于三维弥散问题，这类边界条件可表为： 其中A是所研究的（空间）区域（ ）的边界，一般是曲面，有时还可能包括个别属于边界的点或线段。 对于二维弥散问题，可表为： 其中 是所研究的（平面）区域（D）的部分或全部边界，一般是曲线，有时还可能包括个别属于边界的点。 对于无穷远处的边界，若所研究的溶质在那里没有源或汇时，可采用那里溶质的原始浓度（c0）作为边界条件。对于｛xoy｝平面内的二维弥散问题，可表为： （2）第二类边界条件 已知（或给定）研究区域边界上溶质通量的边界条件，称为第二类边界条件。 这类边界条件的基本物理关系为：在边界曲面的任一点处，界面两侧的溶质质量通量（Ia与Ib）在界面法线方向上的投影必须相等（参见图3-4-1），即： （3-4-2） 其中： （3-4-3） （3-4-4） 分别表示界面（a）侧与（b）侧的溶质质量通量（包括对流通量cV与弥散通量J）。 下面举几个第二类边界条件的例子。 ①两种不同空隙的多孔介质接触面边界 设研究区的边界是两种不同岩性的接触面，研究区内部多孔介质的空隙率为na，外部的空隙率为nb。在这种界面上的边界条件有两个： （a）从界面两侧逼近于界面上同一点时溶质浓度是相等的，即： （3-4-5） （b）界面两侧溶质通量的法向投影是相等的。 对多孔介质的单位面积而言，这种相等关系的一般表达式为 （3-4-6） 即： 将（3-4-6）式左、右两端展开，就得： 代入式（3-4-6）便得： （3-4-7） 作为边界条件，上式中等号右端的函数须是已知的。要注意，上式中等号两端的cos(n，ui)不能相约消掉；因为约去cos(n，ui)后的结果为： （3-4-8） 这已是完全不是原式 所表达的物理意义了。出问题的原因在于约去cos(n，ui)后的（3-4-8）式中，i已成为单指标（指定指标）而不是双指标（求和指标），从而就不能通过求和来保持（3-4-7）式的物理意义。 特别地，对于{xoy}平面内的二维渗流三维弥散问题，只要令（3-4-7）式中的Dxz= Dxz= Dzx= Dyz= Dzy=0及 ，得到的结果就是相应的边界条件的表达式。对于一维渗流三维弥散问题，当取x方向与V 方向一致时，则令（3-4-7）式中的Dxy= Dyx= Dxz= Dzx= Dyz= Dzy=0及Vy=Vz=0，就得到相应的边界条件的表达式。 对于{xoy}平面内的二维渗流二维弥散问题，其边界条件显然为： （3-4-9） 对于一维渗流一维弥散问题，若空隙率为na与nb的岩层分界面与渗流方向是正交的，这时取 x方向与V方向一致，则这样的边界条件就具有最简单的形式 考虑到这时Dxx=DL（=aLV+DdT″），naVa=va，nbVb=vb且va=vb，ca=cb，则上式可进一步简化为： （3-4-10） 上述（a）、（b）两个边界条件与不同K值岩层分界面上H1=H2与v1n=v2n，这两个条件是完全类似的。 ②隔水边界 隔水边界对于水流是隔绝的，即 vn=0（n是隔水边界的法线方向），但对于溶质不一定是隔绝的，即 。例如当隔水边界上存在着溶质的吸附、溶解或放射性衰变时就是如此。 假定隔水边界对于溶质也是隔绝的，则在这种边界上就有 。这时只要令（3-4-7）式等号右端为零，得到的结果 （3-4-11） 就是相应的边界条件。由于边界是隔水的（ ），所以 从而（3-4-11）式可简化为 （3-4-12） 当取坐标轴方向与弥散的主方向一致且n方向与某个坐标方向（例如z方向）一致时，（3-4-12）式可进一步简化为： 即： （3-4-13） 对于一维渗流一维弥散问题，若隔水边界与V方向是正交的，这时若取x方向与V方向一致，则由（3-4-10）式即得相应的边界条件为： （3-4-14） ③多孔介质与流体介质接触面边界 当含水层与河、湖、海洋、水库、渠道等水体毗邻时，含水层与水体的接触面就属于多孔介质与流体连续介质接触面。 设 为流体连续介质中溶质的浓度，c为多孔介质中同种溶质的浓度。一般来说， 。这里我们假定 是确定的且不依赖于 。 按照基本关系（3-4-2）式，由流体连续介质一侧通过对流与分子扩散两种方式输运到单位界面面积（包括空隙面积和颗粒面积）的溶质质量 必须通过对流与水动力弥散两种方式从该单位面积界面输运到多孔介质中去，即 必须等于 这样就得到了这类接触面上边界条件的表达式 或写作 （3-4-15） 若水体中溶质浓度是均匀分布的，或通过长时间充分混合达到了均匀分布时，水体中的分子扩散就不存在，这时（3-4-15）式简化为 （3-4-16） 一般说来，边界上 ，故由上式推知 ，即界面两侧同种溶质的浓度是不相等的，这一点是要注意的。当 保持常数并经过充分长的时间后，可能会出现 这种类型的边界条件。 在一维渗流一维弥散的情况下，若取x轴与V方向一致且假定边界与x轴垂直时，（3-4-16）式简化为： （3-4-17） 当 成立时，由上式可见在边界上也有 成立，这与（3-4-14）式是类似的。 ④渗出面边界 这类边界的内侧为多孔介质，外侧为空气，溶质随渗透水流由这类边界泄出。这里仍用 表示空气一侧靠近边界很薄的水层内溶质的浓度，则这类边界两侧显然有 成立。又由于空气一侧不存在溶质的分子扩散，故由（3-4-16）式可得到渗出面边界条件为： （3-4-18） 在一维渗流一维弥散的情况下，仍取x轴与V 方向一致，并假定边界与x轴垂直，则由（3-4-17）式可见，相应的边界条件仍为 。 以上（3-4-7、11、12、15、16、18）各式中都使用了求和约定。 3.5  简单问题的解析解 与渗流问题类似，弥散方程与其定解条件的总体称为弥散问题的数学模型。本节给出几个简单弥散问题数学模型的解析解。 3.5.1 半无限含水层中的一维问题 1. 假定条件及数学模型 在半无限含水层的始端连续稳定地注入浓度为c0的示踪剂，假定含水层中初始示踪剂浓度处处为零，渗流是稳定均匀流，流速为v=nV，弥散是一维的，含水层中示踪剂无源无汇。这种情况常见于室内，在相当长的砂柱中进行一维弥散研究，从砂柱的一端用含示踪剂浓度为从c0的水驱替原来不含示踪剂的水，就属于这种情况，参见图（3-5-1）。在野外，近似的条件出现在污染的河水切割含水层，而含水层中始终保持稳定均匀一维流的情形。 取x方向与渗流方向一致，这时 ， ， ，则对流一弥散方程（3-5-7）就简化为： 所论问题的数学模型可表为： （I） 2. 模型的解 使用Laplace变换可求得上述数学模型的解为： （3-5-5） 当x充分大，从而使得 时，式（3-5-5）右端的第二项与第一项相比可以忽略不计，若令：                          则有： （3-5-6） 式中，erfc为余误差函数，其定义式如下： 故有：                                                  （3-5-7） 当已知纵向弥散系数DL时，利用式（3-5-5）或（3-5-7）可以方便地求出示踪剂的时空分布c(x,t)。 3. 根据弥散实验资料计算弥散系数 当通过室内或野外弥散实验获得c(x, t)的观测资料时，也可利用上述的公式求得DL，进而计算出αL，具体计算如下： 对式（3-5-7）作代换： 。因为， ， ；同时，当 时， ；当 时， 。因此式（3-5-7）可表示为： （3-5-8） 可以证明：                                故有：                                          （3-5-9） 已知正态分布函数（数学期望μ为0，标准差σ为1）的密度函数为： （3-5-10） 相应的正态累积分布函数Ф0,1(x)为： （3-5-11） 与式（3-5-9）比较可知： （3-5-12） 令                                        查标准正态累积分布函数表有： 当 时：                  当 时：                  由此可见，与 和 相应的x值(x0.159和x0.841)之差为： 故：                                                    （3-5-13） 由式（3-5-13），只要根据弥散实验资料作出 ～x关系曲线（参见图3-5-2），并在该曲线上查得x0.159和x0.841，即可由 求得： （3-5-14） 当忽略分子扩散时，可求得多孔介质的纵向弥散度αL为： （3-5-15） 由于x=Vt时， ，故此时 。这说明相对浓度 出现的位置为x=Vt，由此可见示踪剂过渡带的中心部位是以平均速度V向前推进的。 利用距示踪剂注入处充分远的观测孔在不同时刻的观测资料，绘制 ～t关系曲线（图3-5-3），也可计算DL。设该观测孔出现 和 的时间分别为t0.159和t0.841，则根据前面的讨论显然有： 当 时：      ， ， 当 时：      ， ， 故：                            因此有：                                            （3-5-16） 3.5.2  瞬时点状注入示踪剂的二维弥散问题 设在平面上无界的均质等厚各向同性含水层中存在着流速v=nV的均匀流。取x方向与V方向一致，y方向与之垂直。 这一瞬时在通过坐标原点（0，0）的整个含水层厚度范围内注入示踪剂，单位厚度内注入的质量为m，注后即停，从而发生了二维弥散。假定含水层中原先不含所注入的示踪剂，源汇项为零（参见图3-5-4）。该问题的数学模型可表为： （II） 上述数学模型的解为（孙纳正，1988）： （3-5-21） 忽略多孔介质中的分子扩散作用，以DL=αL V、DT=αT V代入式（3-5-21）得： （3-5-22） 式中，αL 和αT 分别为纵向和横向弥散度。根据式（3-5-22）可方便地求得所给条件下示踪剂的时空分布，图3-5-5给出了不同时刻的计算结果。 3.5.3  连续点状注入示踪剂的二维弥散问题 基本假设与瞬时点状注入示踪剂的二维弥散问题相同，不同的是从t>0开始，以速率q向含水层连续注入示踪剂浓度为c0的流体。这时，只需要把数学模型（II）中的式（3-5-19）改为c(0, 0, t)=c0，便可得到该问题的数学模型。因为连续注入可看作为一系列的瞬时注入，故只需求出式（3-5-21）的时间积分，即可得到本问题的解为： （3-5-23） 令                                  ；                           （3-5-24） 同时进行变量代换 ，可将积分（3-5-23）变为： （3-5-25） 或                                （3-5-26） 其中                                                    （3-5-27） 即为我们熟知的Hantush越流井函数。 让 ，可得到这种情况下的渐进浓度分布为： （3-5-28） 其中K0为第二类零阶修正的贝塞尔函数。 若忽略分子扩散，即设DL=αL V、DT=αT V，则可将式（3-5-26）改写为： （3-5-29） 式中                                ； 把该问题的解表示为式（3-5-26）或（3-5-29）的形式使用起来是很方便的，因为在一般的地下水动力学著作中很容易找到Hantush越流井函数的表格，通过程序对该函数进行计算也不是一件很困难的事。 图3-5-6表示出了使用式（3-5-29）计算的示踪剂浓度随时间的变化。 图 3-5-6  连续点状注入示踪剂时等浓度曲线随时间的变化 3.5.4  径向弥散 设在水平、等厚、无限展布的均质各向同性承压含水层中有一口完整井，井半径为rw，通过它以定流量Q向含水层连续注入示踪剂浓度为c0的水。在初始承压水位水平、无天然流速的情况下，井附近会很快形成接近稳定的一维径向流。此时，通过以井为中心、半径为r的圆柱面的水量为： （3-5-30） 式中，K为含水层的渗透系数，M为含水层厚度。因此，含水层中地下水的平均流速为： （3-5-31） 这里，n为孔隙度， 。这样便可得到本问题的数学模型为：