第⒉章 达西定律与对流迁移

        本章简要介绍达西定律与地下水流方程，并引入对流迁移的基本概念。

        首先建立对流迁移方程，然后介绍对流迁移问题的两种求解方法——欧拉法与拉格朗日法。

接着通过一个简单的例子说明欧拉（质量守恒）法引起的数值问题的性质及其产生原因，并由此引出拉格朗日（粒子追踪）法。

2.1粒子平均速度与运移时间

        

图2.1示意一个导管,

其中一段填满砂；

具有均匀密度与黏滞性的水在一定压力下流过该管，

导管始终充满水，

砂处于饱和态。

设单位时间内流过导管的水的体积流量为Q。

在砂段的两端分别安装测压管测量水头h,

即基准面以上测压水位的高度；

管的横截面积记为A，

砂段长度记为L。

则通过砂段的流量为：

                                           

        其中K为砂的渗透系数；

        h1为上游断面的水头；

        h2为下游断面的水头。

        (2.1)式是达西定律的一种表达形式,

刻画了多孔介质中水流运动的基本规律(Darcy，1856)。

        表达为导数或水力梯度，图2.1对应的达西定律可写为：   

     

图2.1 流量为Q的砂管
 

        其中l为沿该管的距离；

为水力梯度。

        水头h表示测量点处单位重量水体所具有的势能。

砂段内任意两点间的水头损失表示单位重量的水在这两点之间运动所需的能量或功。

注意砂段中任意一点的水头包括两项:

其中，

z为该点相对于基准面的高度

，P为压强，

p为水的密度，

g为重力加速度，

P/(ρg)为连通到该点测压管中的水柱上升到的高度。

熟悉水力学的读者可看出(2.3)式中不包括动力项v^2/2g，

其中v为水的速度。

这是因为孔隙介质中水流速度很小，

可以将此项忽略。

达西定律还可以表达为更一般的形式（见2.2节与附录A），

伯(2.1)式与(2.2)式完全可以描述图2.1所示系统的水流运动。

现在我们不仅要知道流过砂段的流量，

还要确定从一个断面运动到另一断面的时间。

假设开始时一部分孔隙中的速度为零，

即其中的水完全静止；

例子可包括“死端”孔隙或完全封闭的孔隙。

我们采用有效孔隙度一词表示余下的孔隙体积在总体积中所占的比率，

即含有可流动水的孔隙体积与砂的总体积之比。

为简化分析，我们假设孔隙内的速度分布完全均匀。

在此假设条件下，

一定时间间隔内砂中自由流动的水会被砂段上游断面注入的水完全替换。

任意时刻砂段中流动水的体积为θAL，θ为有效孔隙度。

计算运移时间时把θAL作为“孔隙体积”，它表示发生流动的孔隙空间。

如图2.2所示，

我们可以把水流视为一个接一个沿着导管运动的孔隙体积，

每一个孔隙体积被与其相邻的下一个孔隙体积替换。

（注意，

在进入砂段前及流出砂段后每个孔隙体积沿导管占据一定的距离$L'$，

其中$L'=Lθ$。）

图2.2 图2.1中的管段及砂段中的孔隙体
 

从砂段中替换一个孔隙体积需要一定的时间，

记为$\Delta t$；

根据定义，

被替换的液体体积除以该时间为流量$Q$:

或:
 

如果把图2.2中连续的孔隙体积看做被液体界面相分隔，

明显的就是其中某个界面完全流过该砂段的时间；

因此，

可以解释为一个流体粒子穿过长度为的砂段所需要的时间。

更具体地说，

若给定体积为的含水层中流量为，

则为某流体界面沿流动路径通过体积$V$所需的时间。

我们定义通过图2.2中砂段的平均渗流速度，或孔隙速度为：

将(2.2)式代入(2.6)式，得：

在达西定律表达式中频繁地用到达西速度,q或Q/A:
 

注意将达西速度与(2.6)式或(2.7)式中的渗流速度区分开。

达西速度为砂段中单位截面面积的过水量，常称为单位流量。

它表示图2.2管中砂段以外导管中的粒子平均速度，其值比砂段对应的渗流速度小0倍。

从另一方面来看，砂段中水流的截面积(0A)小于空管的截面积(A)。

由于各管段的流量相同，砂中的速度必然大一些。

还应该注意，

虽然我们假设孔隙里有效孔隙中的水流速度均匀，

但渗流速度v并非砂段中的真实速度。

流体粒子在砂段中会绕砂粒移动，

其在砂段中的实际移动距离大于L。

因此，实际粒子速度必然大于L/△t。

所以渗流速度是沿着水流运动总方向线性距离上的表观速度；

其值可以通过野外和实验测得或算出。

在结束这部分内容之前，

还应该认识到孔隙介质中的流动并不像前面假设的那样简单。

明显地，

不可能像前面的假定一样，

一部分孔隙中有理想均一的速度，

而其余部分速度为0。

实际上，孔隙对溶质迁移的影响比上述情况复杂得多。

但是人们发现若将实测总孔隙度（孔隙总体积与样本总体积之比）代入(2.5)式，

得到的运移时间与示踪实验结果不同，

因此必须采用一个较小的孔隙度调整计算结果。

这个孔隙度就是有效孔隙度。

下一章将说明孔隙中的实际速度变化范围很大。

当“孔隙体积”（用有效孔隙度来衡量）中的液体发生替换时，孔隙空间中一些部位的间隙水会发生多次替换，而其余部分只发生部分替换或根本未替换。

因此有效孔隙度是一个协调运移时间计算值与观测值的孔隙度。

通过示踪剂的观测运移时间计算出渗流速度，

达西速度与计算所得的渗流速度之比为有效孔隙度。

这种方法有多个意义。

在讨论一般溶质迁移过程时，

我们会更全面地讨论该问题。

对溶质迁移计算中孔隙度的进一步讨论在第3章进行。

若无特别说明，本书中的孔隙度及字符0均指本节所介绍的有效孔隙度。

回过头来讨论运移时间、流体体积与流量的关系。

图2.2中的导管可类比于地下水稳定流流管，与以河面或没有水流穿过的界面为界限的流场相似。

如果流管系统处于稳定状态，就可以求出运移时间。

如图2.2所示，流管中发生替换的水体积除以通过该流管的流量可得出运移时间（关于流管与流函数更详尽的讨论见附录B）。

 

2.2 广义达西定律与地下水流方程

实际情况中，

如图2.1所示的达西速度与渗流速度只在一个方向的情形很少。

(2.7)式给出的渗流速度与(2.8)式给出的达西定律形式简单，但是不满足常见问题的需要。

本节给出适用于三维问题的两种达西定律表达式。

有关达西定律更详细的讨论可参阅地下水教材（例如：Freeze和Cherry,1979；de Marsily,1986；Fetter,1994；Domenico和Schwartz,1998）。

达西速度的三维矢量表达式为：

其中q为达西速度矢量;

i,j ,k 分别为z , y,方向的单位矢量。

4.,qy 4。为r , y,z方向的达西速度分量。

若孔隙介质采用渗透系数的三个主方向描述,且主方向的轴与坐标轴重合,对具有均匀密度与均匀黏性的水,达西速度的三个分量可表示为:   

  修改标点符号

其中，Kx，Ky，Kz表示各个方向上的渗透系数。

若水的密度、黏性是不均匀的，上述关系式可用固有渗透率和压强梯度表示为：

其中，u为水动力黏度

；kx，ky，kz为三个方向上的固有渗透率；

z是纵坐标，向上为正。固有渗透率只与孔隙介质本身性质有关，

而渗透系数还受流体性质影响。两者有以下关系：

上式为α方向上渗透率与渗透系数的关系，

同理可得y、z方向上的关系式。

由水头h与压强P的关系方程(2.3)式可知，

当密度与黏性为常数时，(2.11)式可简化为(2.10)式。

应该指出，渗透系数和渗透率为二阶张量。

若渗透系数和渗透率的主分量与坐标轴不重合，每个方向上的速度分量就是所有三个方向上水头与压强梯度的函数，而不仅仅是一个方向的函数。

此时以水头作为变量，(2.10)式写成：

其中，

Kx，Ky，Kz为渗透系数张量的主分量；

Ky，Kz，Kxy，Kyz，Kzx为该张量的交叉项。

实际中很少直接用(2.13)式。

通常假定渗透系数张量的主分量与水平及垂直坐标轴相重合，此时渗透系数张量的交叉项为0。

 

(2.9)式算出的达西速度为矢量点函数，

平均渗流速度亦为点函数，

即将达西速度除以（标量）有效孔隙度：

或:

为得出(2.15)式中的渗流速度分量，

首先要求解地下水流微分方程，

给出水头与压强的空间及时间分布；

再由(2.10)式或(2.11)式求出研究区域与研究时段内的达西速度分量；

最后利用孔隙度将达西速度换算为渗流速度。

若水的密度与黏性均匀不变，可以写出以水头表示的最简单的地下水流微分方程：

其中，

q为流体的源/汇项，

即流进或流出单位体积含水层的体积流量；

S为单位储水量，

即水头下降单位值时单位体积含水层释放出的水的体积。

含水层有三个渗透系数主分量，它们的排列方向与坐标轴相吻合。

若流体的黏性与密度发生变化，

水流方程采用压强项表示，

并用质量流量代替体积流量更为方便。

对流体黏性与密度变化条件下的地下水水流方程与溶质迁移的进一步讨论可参考第15章与附录A。

2.3对流迁移

流体粒子的平均速度与运移时间在2.1节已讨论。

若地下水中流体粒子的物理、化学特性相同，本节研究的内容将没有意义。

但是流体粒子的性质一般不会完全相同。

虽然一些粒子与地下水能充分混合，但仍然可以被“标记”出来，如溶解的示踪剂、高盐分的水、溶解的有机污染物等。

这些“标记”出的粒子的运动即溶质迁移

。

本节根据2.1节的假设

，

即所有的溶质粒子以平均渗流速度运动，

推导溶质迁移方程。

所得方程即对流迁移方程。

其他过程对迁移方程的影响将在第3章与第4章讨论。

 

2.3.1质量守恒定律与求解对流迁移的欧拉法

首先根据质量守恒定律建立对流迁移方程，

并以溶质浓度作为基本变量。

设图2.2的砂段中溶解盐或示踪剂的浓度为C，

即单位体积水中溶质的质量。

在图2.2的砂段中，通过横截面的溶质质量为：

 

其中，

Qm为单位时间通过截面的溶质质量；

Q为通过砂的体积流量（单位时间的体积）。

若以达西速度矢量q代替体积流量Q，则对流质量通量的矢量式可以记为：

 

其中，

qm为矢量，

方向与达西速度一致，物理意义为单位时间单位截面积上的对流迁移溶质质量。

图2.3的设置与图2.1、图2.2相同，

图中所示两断面与流动方向正交，间距为△l；

两截面间的体积为A△l，A为管的截面积。

假设盐分或示踪剂的浓度沿管轴方向发生变化（在垂直于管轴的截面上浓度为定值）。

设图2.3上游截面的浓度为C1，下游截面的浓度为C2，那么进入A△l段的溶质质量流量Qm1等于QC1，流出该段的溶质质量流量Qm2为QC2。

因为这两个溶质质量流量不同，A△l段内的溶质质量必然随时间变化，变化率等于流入与流出质量之差。以M表示A△l段的溶质质量，dM/dt为质量随时间的变化率，可得：

其中，q为达西速度标量值。
 

图2.3 图2.1中的砂管及溶质进入或离开体积元的迁移过程

为简化分析，

设对流迁移过程中溶质始终保留在流动的水体内，

即溶质不会因为扩散进入静止或近于静止的孔隙中，

或从其中扩散出来

。

即假定沿孔隙壁与死端孔隙扩散进入静止水体中的溶质可以忽略，

下一章将介绍如何从其他角度来进行考虑。

在此假设下，

A△l段内含溶质的水体积为0A△l，

根据定义，

0为有效孔隙度；

C为该体积内的平均浓度，

两横截面间任意时刻的溶质质量为0A△lC。

因此，溶质质量随时间的变化率可以用浓度随时间的变化率表示：

 

对(2.19)式与(2.20)式的dM/dt项建立等式，得到：

沿管轴浓度梯度为dC/dl，因此浓度差C2-C1可以记为：
 

其中，l沿流动方向取正。联立(2.21)式与(2.22)式得到：

或:

进一步考虑更一般的情况，

设达西速度在图2.3的两个横截面间发生变化。

这表明A△l段内的水量发生变化，

即产生弹性储水效应，

这种情况与承压水系统的储水效应相似。

为简化起见，

设弹性释水效应不会对A△l段内的孔隙体积产生显著影响，

即弹性储水效应仅仅对应于因压力变化导致的刚性孔隙结构内的水密度变化。

在此条件下，(2.24)式可记为：

此时单元体界面间(qC)值的变化导致溶质质量累积变化。

还可以建立更一般的对流迁移方程，

如图2.1砂段中所示的小体积元，

其任意一边与流动方向不相正交，

但与垂直-水平笛卡尔坐标系的方向相同。

该单元体△x△y△z如图2.4所示。

单元体的任一坐标轴都不与管轴平行或正交；

因此，

达西速度矢量q在三个坐标轴上都有分量，

溶质浓度在三个方向也相应发生变化。

为简化分析，

我们再次假设弹性储水效应只引起刚性孔隙结构中水密度的变化。

体积元中垂直于x方向的截面积为△y△z；

根据推导(2.25)式的方法，

可得x方向上流入减去流出的净溶质质量变化量为：

其中，

qx为达西速度在x方向上的分量。

同理可得y方向与z方向上的表达式，

体积元的入流减去出流的净质量变化总量为：

图2.4 笛卡儿直角坐标系中建立对流迁移方程的体积元
 

若图2.4的体积元中有源使得水进入流场，

或有汇使水流出流场，

还应该在表达式中加入源/汇项，

表示流入、流出单元体的质量。

以Qs表示源/汇的体积流量，

设源为正，汇为负，

Cs表示源/汇的溶质浓度，

则由源/汇进入或流出体积元的溶质质量为QsCs，

表示每单位时间通过源/汇的质量。

将该项加进表示体积元的流入减去流出量的(2.27)式中，得到：

 

体积元内的质量累积变化率为0Ax△y△zC/dt。

与(2.28)式建立等式关系，得到：

此时设体积元△x△y△z内水的体积变化量可忽略，

并且q1=Qs/(△z△y△x)是由于源或汇引起的单位体积含水层的体积流量。

将(2.29)式除以孔隙度，

则达西速度分量变为渗流速度分量v1、v2、v3，得到上式的另一种表达形式：

 

(2.30)式的矢量形式为：

下标形式为：
 

其中，

下标i表示坐标方向，

并表示要对所有坐标方向(α,y,z)求和。

以上(2.30)式～(2.32)式为对流迁移方程的一般表达式，

适用于有源/汇项、流速及浓度在三个坐标方向上变化的对流迁移。

由推导过程可以看出，

若坐标系本身与水流方向不重合，

即使流速和浓度梯度是单向的且线性相关的，仍需要采用一般表达式。

建立(2.30)式～(2.32)式时假定含水层为刚性，

即储水量的变化是由于间隙水的膨胀或压缩引起的。

承压储水量的变化通常会引起水密度的变化以及含水层变形(Cooper,1966)。

根据质量守恒定律分析变形含水层比研究刚性含水层复杂得多。

如果含水层发生形变，孔隙度0则随时间变化。

对固定空间体积元建立质量守恒方程时，孔隙度将出现在溶质质量累积变化率对时间的导数中。

例如，在推导(2.29)式的过程中，体积元△x△y△z中的溶质质量变化率应写为△x△y△zα(0C)/dt而不是△x△y△zC/dt，因此有：

若此时(2.33)式左边的单位给水量q1、q2、q3表示流过变形含水层体积元各个方向的真实流量，可以得出溶质质量守恒方程的完全表达式。

然而，

由于假设体积元△x△y△z的空间位置确定，

且因为孔隙介质的固体颗粒本身基本上不可压缩，

只有在体积元边界处固体颗粒以有限的粒子速度u移动才会引起含水层变形及孔隙度变化。

这反过来意味着在变形过程中，

孔隙介质中的一些孔隙空间及该孔隙中的间隙水会以速度u穿过体积元边界。

Cooper(1966)指出，

达西速度及渗流速度与孔隙介质中的颗粒或固体基质有关。

因此，相对于固定参考系，变形含水层中间隙水的运动总速度是渗流速度与颗粒速度之和；

同样，

总的单位流量，

即每单位截面流量为达西速度与ug之和，

后者是颗粒速度与孔隙度之积。

在净溶质流入量的表达式如(2.30)式、(2.33)式的左边，

必须用这些总速度或总单位流量导出变形系统的质量守恒方程式。

实际上，承压储水量的变化对溶质质量守恒的影响很小，在刚性含水层假定下推导的迁移方程式满足一般情况。在不能近似表示的情况下，可以在控制迁移方程中增设源/汇项，模拟承压储水作用，令源的强度等于与含水层变形有关的储水累积速率或水量释放的速率，也就是说，等于刚性含水层储量变化之外的储水率或释水率。

野外最具代表性的储水过程是非承压含水层的水量存储，即潜水位或自由水位变动引起的储水效应。在地下水流偏微分方程中，把潜水储量的变化表示为随时间变动的边界条件。但是在实际计算中，通常把潜水效应作为含水层上部的储水量变化来处理，而不作为边界条件。因此，把它们放在(2.16)式的右边，但此时单位储水量较大，通常取单位给水率除以潜水面出现的水文地质区间的厚度。

(2.30)式～(2.33)式是针对含水层中完全饱和的含水单元导出的，不能反映含水层单元有自由水面时溶质的质量守恒关系，也不能解释自由水面上升或下降引起的溶质质量变化。迁移模拟中，含自由水面的单元或节点的近似方程采用含水层的饱和厚度而不是固定厚度△z。当水位变化时，需要对不同时间段更新饱和厚度；这样在计算中才能考虑到储水量的体积变化，以及相应的溶质质量变化。

 

在(2.30)式～(2.33)式中，偏微分式∂C/∂t表示空间某点溶质浓度C的变化率，这些公式称为欧拉表达式(例如Daily和Harleman,1966)。这类方程式可以直接采用以质量守恒定律为基础的常规数值法求解。然而，欧拉型对流迁移方程直接求解时会遇到严重的困难，称为数值弥散。此处只示意数值弥散，第6章、第7章将详细讨论这类问题以及与迁移模拟有关的数值方法遇到的问题。

再回到图2.3的设置，设初始时刻流动系统内盐的浓度为10 mg/L。设砂段的达西速度为10 cm/h，管的横截面积为100 cm^2，那么砂段中的流量为1,000 cm^3/h。设△l长为100 cm，砂的孔隙度为0.1；因此体积元A△l的孔隙体积为1,000 cm^3，即1 L，流量可表示为每小时一个孔隙体积。最后，设t0时刻瞬时注入上游入口的水突然变成浓度为100 mg/L的盐溶液。

现在用质量守恒定律对A△l进行分析。在t0时刻之前，通过上游断面的盐浓度为10 mg/L，下游断面处流出的盐浓度相同；因此，该段盐的质量累积率为零。在t0时刻，上游的浓度变为100 mg/L，流入的盐量增至100 mg/h。根据(2.25)式并将偏微分式∂C/∂t写为有限差分式△C/△t，其中△C表示有限差分时间步长△t上体积元A△l内盐的平均浓度变化量，△t取6 min，即0.1 h；假设任一时段流出该单元下游断面的浓度等于每个时段开始时对应的平均浓度值。因此，在t0之后第一个0.1 h进入体积元的溶质质量流量为100 mg/h，流出的溶质质量流量为10 mg/h。结果得到第一步长的质量累积变化率为90 mg/h；该值表示起始步长段体积元A△l中的盐增加了9 mg，浓度增加9 mg/L。由此可以计算出第一步长结束时平均浓度为19 mg/L。第二个步长段流入的质量流量仍为100 mg/h，这时流出的质量流量为19 mg/h；因此，体积元内质量累积率为81 mg/h。在该累积率下，第二步长结束时体积元A△l内的浓度增加8.1 mg/L，达27.1 mg/L。表2.1为采用该法得出的体积元在每个时间步长结束时的平均浓度计算值。
 

上述计算过程中采用有限差分法求解(2.25)式并运用向前差分近似表示时间导数。因为假定下游面的浓度为体积元A△l的平均浓度，表2.1可以看做是下游面浓度对时间的函数。把表2.1中的浓度-时间关系绘成曲线，结果得到图2.5中的一条光滑曲线，出流浓度由10.0 mg/L渐增，最终达到入流浓度100.0 mg/L。
 

理论上，假定只有对流迁移时，在t0时刻后1.0 h，下游面浓度应该由10.0 mg/L突增至100.0 mg/L(图2.5)。可以清楚地看出，用质量守恒定律计算引进了出流面的浓度-时间曲线，也称穿透曲线的抹平效应。在下一章将看到实际上穿透曲线根本不会呈现完善的突变形态；天然流速变化过程中会发生传播或弥散。这里的计算效果显示了这一情况，因而称之为数值弥散。应该指出，数值弥散是计算上的人为效应，其本质与物理机制无关。

图2.5 纯对流迁移例子里的计算出流浓度与理论出流浓度(假定入流浓度为常数)

如果将图2.3中的整个体积元A△l作为我们计算的基础，将其分为一系列垂直于管轴的薄片，根据质量守恒定律计算，会得到较好的计算结果。如果对体积元A△l采用1 h作为计算步长，或用一个接近于流体粒子通过该计算体积的运动时间值，也可以减小数值弥散。这说明数值弥散效应受时间与空间离散方案的影响。

2.3.2 求解对流迁移的粒子追踪法

因为根据质量守恒定律进行计算引起了数值弥散，因此常用粒子追踪法解决对流迁移问题。在粒子追踪法中，每个流体粒子对应一个浓度值，通过流场分析可以计算这些粒子在研究区域内的运动。粒子追踪法是根据拉格朗日的观点建立的(例如，Daily和Harleman,1966；Bear,1972)，该法中浓度不与空间固定点或单元体对应，而是以优势流动速度运动的流体元或粒子的浓度。

对流迁移的拉格朗日方程根据质量守恒方程的下标形式[即(2.32)式]建立。注意式中的C现在表示随流体流动的流体元或一个粒子的浓度。将(2.32)式左边第一项按乘法定律展开可得：

 

由于下标记法表示各坐标方向的和，∇·v实际上表示流体散度或单位孔隙体积的净出流量。对于稳定流，该项必然等于源/汇项的变化量q0，即单位孔隙体积由于源、汇引起流体增加或损失的体积变化率：
 

将(2.34)式与(2.35)式代入(2.32)式，并整理得：

(2.36)式的左边表示流场中沿流线(也称为特征线)流动的粒子或流体的浓度变化率。此两项和称为随体导数，记为DC/Dt：
 

由此，(2.36)式可写为：
 

(2.38)式可以进一步简化。根据拉格朗日法的原理，浓度C对应于运动单元或粒子。注意只有对流迁移时，由源点流出的粒子其浓度必然与源的浓度相同，同时进入汇的粒子其浓度必然等于该粒子在汇处的浓度。因此按照拉格朗日法原理，(2.38)式中的C与C'相等。所以(2.38)式可简化为：

(2.39)式表明，纯对流运动流线上的运动粒子其浓度不随时间变化。粒子沿流线运动，其浓度保持在流线起点处的浓度。粒子由源出发沿流线运动，在整个流线上维持着源处的浓度；类似地，进入汇点的粒子沿流线运动到汇，进入汇点的浓度与沿流线的浓度相同。因此，以拉格朗日法处理纯对流迁移问题就是确定流动系统中溶质粒子的流线。

回到上一节讨论的图2.2中盐溶液的运动问题，现在用拉格朗日粒子追踪法进行研究。开始时在某一界面放很多粒子，表示浓度突然从原来的10.0 mg/L突变至100.0 mg/L。那么我们来计算该界面在砂段中的推进，从初始时刻t0开始，渗流速度为100 cm/h，界面之后的浓度为100.0 mg/L，界面之前的浓度为10.0 mg/L。根据此法计算得到出流面在t时刻起一直至1.0 h浓度仍为10.0 mg/L，而在1.0 h时突增至100.0 mg/L。注意在计算过程中不产生数值弥散，对于对流占优的迁移问题，由该法得到的计算结果比质量守恒法(欧拉法)的结果更符合观测情形。

在本例中，图2.2的砂管有两个区域，一个匀速前进的平面状界面将它们分隔开来，每个区域中的浓度为定值。粒子追踪法的应用范围十分广泛。例如，对于浓度分布可以采用一系列三维等浓度面描述的问题，可以把每个浓度面作为一个界面。利用渗流速度与运移时间的概念计算出不同时刻每个等浓度面的位置，获得每个计算时刻新的浓度分布。虽然用此法能获得新的浓度分布，但是该法的计算变量不是浓度；计算中只涉及用流体粒子表示的不同界面的位置、速度及运移时间。

界面的计算工作量取决于界面的几何形状及水流形式的复杂程度。图2.2中界面始终是与管轴垂直的平面，并且以均一速度沿已知方向前进。然而野外多为三维流场并且流场随时间变化，界面位置可能更加复杂；图2.6示意了一个简单例子。计算一个界面的移动涉及确定界面上大量粒子的位置。图2.2中所有粒子的速度不随时间变化并且为定值。粒子速度只需要计算一次，不需要对大量粒子进行跟踪。而对于野外复杂的三维界面，界面上单个粒子的速度不同，并随时间变化。因此通常采用(2.9)式和(2.10)式的达西定律三维矢量形式分别计算每个粒子的流速，并且需要计算粒子在不同时刻的分布情况。时间步长增加后，必须使用新算出来的值确定每个粒子的位置，该过程中要用到流速、时间增量、粒子在前一时间的位置；通常需要占用大量的计算机内存来存储粒子位置的信息。虽然该计算方法十分烦琐，但与计算图2.2中界面在管中移动的方法本质上相同。

粒子追踪法以及解决纯对流迁移问题的计算程序将在第6章详细讨论。

 

图2.6 纯对流条件下逐渐远离连续点源的溶质晕的三维锋面
 

深入阅读与思考题：

2.1 Bear (1972, 第70页) 在讨论欧拉法与拉格朗日法的异同点时进行了这样的描述：“描述流体系统中的流动有两种方法。第一种方法是描述单个粒子的历程，为拉格朗日法；而第二种方法着眼于空间上的固定点，为欧拉法。”这就像在大海中漂浮的一艘船，为描述它的运动，可以在船上记录它的运动轨迹（拉格朗日法），也可以在海上布置固定的观测网，监测它的位置（欧拉法）。对欧拉法与拉格朗日法更严谨的分析比较可参考 Bear (1972) 的著作。

2.2 附录A 给出达西定律的推导过程以及密度变化条件下的水流方程。附录B 介绍了流函数的概念以及计算地下水流迹线与对流迁移的方法。这部分内容可以参考这两个附录。

2.3 验证表2.1的结果。并取时间步长为0.2 h 进行验算。讨论步长值的变化对浓度曲线的影响。

2.4 假设在流场中心有一个连续点源（见图2.6），并形成三维球形稳定流场。若该点源的流量为1 m3/d，含水层孔隙度为0.2，根据关系式 dr = v(r)dt，其中v(r)为半径r处的径向速度，求第100天时注入流体的运动距离。