第2章 物理特性和原理

第2章 物理特性和原理

校正:高宏斌,骆慧斌,雷颖,解千里,吴小龙,张远,朱冬楠

统稿:吴小龙,朱冬楠

 

2.1达西定律

地下水水文学作为一门定量科学而诞生,其发展最早可以追溯到1856年。法国水力工程师亨利达西(Henry Darcy)于1856年发表了关于法国里昂市的供水报告。在报告中,达西介绍了他用于分析水流通过砂柱的实验室实验。实验结果被概括为经验公式,并以他的名字来命名。

实验设备如图2.1所示,向一个横截面面积为 A 的圆柱中注满砂,并将两端密封,从上下两个栓塞外分别插装入流管、出流管和一对流体压力计。水流通过圆柱,直到全部孔隙都充满水。从此刻起,入流量 Q 等于其出流量。如果我们在高程 z = 0 处设置一个任意基准点,那么两个压力计入口处的高程则为 z1z2,同时流体液面的高程分别为 h1h2,令两个压力计入口之间的距离为 Δl

图 2.1 达西定律实验的设备图示。

则我们可根据下式来计算通过圆柱体的单位面积流量,也即 v 值:

v = \frac{Q}{A} (2.1)

如果 Q 的量纲是 [L3/TΔl 的量纲是 [L2],那么 就具有一个速度的量纲 [L/T]。达西所做的这个实验表明:当 Δl 保持恒定时 vh1h2 成正比,而当 h1h2 保持恒定时,则与 Δl 成反比,如果我们确定了 Δh = h2h1 (这些特定符号适用于任意情况,并将为我们以后的实验带来便利),则有: v \propto 1/\Delta l。因此达西定律可写为如下形式:

v = -K\frac{\Delta h}{\Delta l} (2.2)

或者变为下面的微分形式:

v = -K\frac{\Delta h}{\Delta l} (2.2)

或者变为下面的微分形式:

v = -K\frac{dh}{dl} (2.3)

在方程(2.3)中,h 称为“水力水头”,dh/dl 是“水力梯度”。K 是一个比例常数。它是一个与圆柱体内土壤特征相关的参数。因为我们已保持水力梯度保持不变,不同土壤间的单位面积流量必定有所不同。换句话说,如果 dh/dl 保持恒定不变,则 vK。参数 K 为熟知的“水力传导率”。对于砂和砾石来说,它的值很大;对于粘土和大多数岩石来说,它的值较小。因为 ΔhΔl 都具有长度单位 [L],所以方程(2.2)的量纲分析指出,K 的量纲是流速 [L/T]。在2.3节中,我们将指出,K 不仅是介质特征的函数,也是通过介质而流动着的流体特性函数。

把方程(2.1)代入方程(2.3),我们可以得到达西定律的另一种表达形式

Q = -K\frac{dh}{dl}A (2.4)

它有时被进一步简化为如下的形式:

Q = KiA (2.5)

其中 i 是水力梯度。

达西定律对空间中任意方向的地下水流都是有效的。参照图2.1和方程(2.3),如果水力梯度 dh/dl 和水力传导率 K 保持不变,则 v 是相对角度 θ 的自变量。这一点,甚至在θ 角度大于90°,水流在压力作用下克服重力穿过圆柱体向上流动时,也是成立的。

我们已经注意到,单位面积流量 v 有速度或通量的量纲,因此,它有时被称为“达西速度”或“达西流”。单位面积流量是一个宏观概念,容易测量。必须把它与单个水分子在砂粒空隙中曲折前进的微观真实流速区分开(见图2.2)。 微观流速虽是真实的,但它却无法测量。在本章的剩余部分中我们将完全以宏观水流概念来进行讨论。不管它的量纲如何,我们不把 v 看作是一种流速,而将运用更为准确的表达—“单位面积流量”。

图 2.2 地下水流的宏观和微观概念。

最后一节可能不是必须掌握的,但它揭示了一个具有根本性重要意义的问题:当我们用达西方法来进行地下水流的分析时,它事实上意味着,我们通过代表性连续体(representative continuum)的概念定义了一系列宏观参数(诸如水力传导率),并利用宏观定律(诸如达西定律)将微观的现象进行平均化,以替代对多孔介质砂粒(或粘土颗粒或岩石碎块)真实情况的微观描述。概念性上看,这一范式简单而合理,但仍立足于复杂的理论基础之上。贝尔(Bear, 1972) 在其最近关于多孔介质和流体的著作中,详细讨论了这些基础理论。在2.12节中, 我们将进一步探讨地下水流中微观描述与宏观描述之间的相互联系。

达西定律是一种建立于经验数据基础上的定律。在使用基础物理学定律推导达西定律方面,有过很多研究和探索。贝尔在1972年的著作对这些研究有较详细的论述。最成功的方法,是将流体力学众所周知的纳维尔—斯托克(Navier—Stokes)方程应用于多孔介质中流体的理想化概念模型。哈伯特(Hubbert, 1956)和艾尔梅(Irmay, 1958)最早作了这一尝试。

达西定律在地下水流分析中具有根本性的重要意义,本书将为此提供依据,但值得注意的是,在很多多孔介质流体的其它应用中它也依然有很高的重要性。它可以描述土壤水的流动,并受到土壤物理学家、农业工程师和土壤力学专家的青睐。石油储层分析专家用它来描述深层地质构造中石油和天然气流动,化学工程师把它应用于过滤器的设计中,材料科学家把它应用于多孔陶瓷的设计中。它甚至曾被生物科学家用于描述流体通过体内孔隙薄膜的流动。

达西定律是一个强大的经验公式,它的分支值得我们投入更多的关注。下面两节将从物理意义上对水力水头 h 和水力传导率 K 给予更详细的说明。

2.2 水力水头和流体势能

涉及流动的物理过程分析,通常需要判明势能梯度。例如,固体中的热流是从高温部分向低温部分流动的,电路中的电流也是沿高电压向低电压方向流动的。对这些过程来说,温度和电压是势能的量,而热和电的流动速率与这些势能的梯度成正比例。我们的任务是确定多孔介质中控制水流的势能梯度。

哈伯特在其经典的地下水流论文(1940)中对该问题进行过论述。在本节开头部分我们将重温他定义的概念和推导过程。

哈伯特对流体势能的分析

哈伯特把“势能”定义为:“能够在流体系统中的每个位置测量到的一个物理量,其特性是,流动永远是从能量值较高处向能量值较低处流动,不管其在空间中的方向如何”。在达西实验(图2.1) 中,由测压计中水的液面所测定的水力水头 h 是满足这个定义的,但却如哈伯特所指出的:“不经过调查研究就采用这个经验,就如同从一个温度计上读取水银柱的长度而不知道温度才是它所要真正揭示的物理量”(795页)。

对于势能来说,两种潜在的表现形式是高程和流体压力。如使用垂直(θ = 0°) 圆柱体来进行达西实验 (图2.1),水流肯定是在重力的作用下通过圆柱体向下方运动 (即从高处向低处流动) 的。然而,如果圆柱体呈水平放置(θ = 90°),那么重力就不起作用,水流只能在一端增加压力而在另一端降低压力的情况下才能发生。无论高程还是压力,只有其中一个都不是准确的势能,但它们都是总势能的组成部分。

对于已经在基础物理或流体力学中接触到势能概念的读者来说,达到我们所追求目标最好的方法是去考察流动过程中的能量关系,这一模式已经变得不足为奇。事实上,势能的传统定义(数学家和物理学家经常提出的)是“在流动过程中所作的功”和“单位质量的流体在一个流动系统的任何两点间运动所作的功,它是这个单位质量能量损失的计量标准”。

流体通过多孔介质流动,是一个力学过程。要驱动流体向前,必须克服存在于流体和多孔介质颗粒间的摩擦力。所以,流动是机械能向热能转化的一种不可逆的过程(摩擦阻力机理)。因此,流体在空间中的流动方向,必然是从流体单位质量机械能较高的区域流向单位质量机械能较低的区域。一个水流系统中任一点处的单位质量机械能可以被定义为:从一个任意选择的标准状态中转移单位质量流体所需要作的功。显然,我们已经确定了一个既满足了哈伯特关于势能(以流动方向来表示)的定义,也满足了其传统定义(以所作的功来表示)的一种物理量。因此,通过多孔介质的“流体势能”就是“流体单位质量的机械能”。

现在来讨论势能与我们前面所说过的高程与压力间的关系。考虑一个在 z = 0 和压力 p = p0 (此处 p0 为大气压力) 时的任意标准状态 (图2.3),密度为 ρ0 的单位质量流体的体积为 V0V0 = 1/ρ0。我们要计算:把一个单位质量的流体从标准状态抬高到水流系统中的 P 点(该点的高程为 z,流体压力为P)所需要作的功。这里单位质量流体的密度为 p,体积为 V = 1/ρ。此外,我们将考虑该流体具有速度 v = 0 (标准状态),同时这个流速 v 是在点 P 处的。

图 2.3 单位质量流体的机械能计算。

这里,功的计算有三个组成部分。第一,把质量从高程 z = 0 抬高到高程 z 所需作的功,即

w_1 = mgz (2.6)

第二,把流体从速度 v = 0 加速到速度 v 所需作的功,即

w_2 = \frac{mv^2}{2} (2.7)

第三,把流体压力从 p = p0 提高到 P 所需作的功,即

w_3 = m \int^p_{p_0}\frac{V}{m}dp = m \int^p_{p_0}\frac{dp}{\rho} (2.8)

如果流体从点 P 向标准状态中的一点流动,方程(2.6)就代表势能损失,同时方程(2.7) 是动能的损失,方程(2.8)是弹性能损失或称为p-V功。
流体势能 Φ (即单位质量机械能)是 w1w2w3 之和。对于流体的一个单位质量来说,方程(2.6)、(2.7)和(2.8)中的 m = 1,可得

\Phi = gz + \frac{v^2}{2}+\int^p_{p_0}\frac{dp}{\rho} (2.9)

对于多孔介质中的水流来说,流速是极小的,因此第2项通常可以忽略不计。对于不可压缩流体(密度不变,因此 ρ 不是 p 的函数)来说,方程(2.9)可以进一步简化:

\Phi = gz + \frac{p-p_0}{\rho} (2.10)

因此,我们前面对流体势能构成的判断,现在看来是正确的。方程(2.10) 中的第一项包括高程 z,第二项包括流体压力 p

但是这些项如何与水力水头 h 发生关系呢?让我们回到达西压力计(图2.4)上来。在点 P,流体压力 p 由下式给出,即

p = \rho g\psi + p_0 (2.11)

其中 ψP 点上液柱的高度。p0 是大气压力或在标准状态时的压力。从图2.4及方程(2.11)可知

p = \rho g(h-z) + p_0 (2.12)

图 2.4 一个实验室测压计的水力水头 h,压力水头 ψ 和高程水头 z

把方程(2.12)代入方程(2.10)得:

\Phi = gz + \frac{[\rho g(h-z)+p_0]-p_0}{\rho} (2.13)

化简后可得:

\Phi = gh (2.14)

这一段冗长的演算可以得出一个最简单的结论,即在多孔介质中,任一点 P 处的流体势能 Φ,可简单通过该点的水力水头乘以重力加速度获得。因为 g 在地球表面的值几乎保持恒定,所以 Φh 就几乎完全是直接相关的,已知其中之一,也就可以求解另外一个,因此,水力水头 hΦ 一样,是势能合适的表征参数。回忆一下哈伯特的定义:“它是个物理量,并且是能够测量到的。永远从高 h 值区域向低 h 值区域流动”。事实上,通过方程(2.14)可知,如果 Φ 是单位质量的能量,那么 h 就是单位重量的能量。

在地下水水文学中,一般大气压力 p0 被设置为零,即用标准压力(以大气压力为起算点的压力)来计量。在这种情况下,方程式(2.10)和(2.14)变为:

\Phi = gz + \frac{p}{\rho} = gh (2.15)

两边同时除以 g,得

h = z + \frac{p}{\rho g} (2.16)

将读表压力带入方程(2.11)中,得,

p = \rho g \psi (2.17)

同时方程(2.16)也变为:

h = z + \psi (2.18)

因此,水力水头 h 就是所测点的高程(或称“高程水头”) z 与压力水头 ψ 两项之和。这个基本水头关系是研究地下水流的基础。图2.4代表了达西测压计中各种水头的相互关系。图2.5是在野外进行测量时的位置。

图 2.5 一个野外测压计的水力水头 h,压力水头 ψ 和高程水头 z

熟悉基础流体力学的读者可能已经认出方程(2.9)与伯努利方程很相似,它是流体流动期间能量损失的典型公式。一些研究者 (Todd, 1959; Domenico, 1972) 用伯努利方程作为他们研究和发展流体势能和水力水头概念的起点。如果把方程(2.9)都变为水头项并加以简单的脚注,得到:

h_T = h+z + h_p + h_v (2.19)

其中 hz 是高程水头,hp 是压力水头,hv 是流速水头,在我们前面的符号标记中,hz = z, hp = ψhv = v2/2ghT 项称为“总水头”。在 h0 = 0 的情况下,它等于水力水头 h,,与方程(2.18)相一致。

量纲和单位

hψz 等项中水头的量纲是长度 [L]。它们经常用“水米(meters of water)”或“水英尺(feet of water)”来表示,这种特殊标以“水”字的原因是强调水头测量要由方程(2.17)中的流体密度来决定。如果在图2.5中,在 P 点给以相同的压力 p,但在地质结构孔隙中的流体是石油而不是水,那么水力水头 h 和压力水头 ψ 就会有不同的值。在本书中,由于我们所讨论的一概是水,形容词短语将被省略掉,而将水头全部以米计之。

至于本节所介绍的其它项,在 SI 系统内是以 [M][L][T] 为基础的,即:压力的量纲为 [M/LT2],质量密度的量纲为 [M/L3],而流体势能,根据其定义,单位质量具有的能量 [L2/T2]。表2.1对此前已介绍的所有重要参数的量纲和常用单位作了分类。参考附录I,就可厘清所有的疑惑。在本书中,我们将采用 SI 制单位作为我们的基本单位系统,但是表2.1中包括 FPS 制的对应单位。附录 SI 中的表A1.3中提供有换算系统。

应该注意,在表2.1中,水的重量密度 γ 由下式来确定,即

\gamma = \rho g (2.20)

它对于 FPS 系统的单位来说,是比质量密度 ρ 更为适当的参数,可作为基本量纲。

表 2.1 地下水基本参数1 的量纲和常用单位

    SI 国际系统2   英尺-磅-秒系统 (FPS)3
参数 符号 量纲 单位   量纲 单位
水力水头 h [L] m   [L] ft
高程水头 ψ [L] m   [L] ft
高程水头 z [L] m   [L] ft
流体压力 p [M/LT2] N/m2 or Pa   [F/L2] lb/ft2
流体势能 Φ [L2/T2] m2/s2   [L2/T2] ft2/s2
质量密度 ρ [M/L3] kg/m3  
重量密度 γ   [F/L3] lb/ft3
单位排泄量 v [L/T] m/s   [L/T] ft/s
水力传导率 K [L/T] m/s   [L/T] ft/s
1见附录I中的表A1.1,A1.2和A1.3.
2基本量纲是长度[L],质量[M]和时间[T].
3基本量纲是长度[L],力[F]和时间[T].

潜水位计和潜水位计组

水头测量的基本设备是软管和硬管。利用它们可以确定出水面高程,在实验室内(图2.4),软管是一个“流体压力计”;在野外(图2.5),硬管称为“潜水位计”。潜水位计的硬管表面是密封的,但在底部开口以便进水,同时在顶部开口与大气相通。其进水部分经常是一段布有孔洞的硬管,或者商用过滤器。在任何情况下,进水部分都应设计成能够进水却不允许地层中的砂粒或粘土颗粒进入。这里要强调一下的是,潜水位计的测量点是水位计的底部,而不是液面。我们可以将潜水位计的作用温度计相比较。几乎可以说,它是在地下水层中的任意点 P 确定 h 值最简便的设备。在近几年的一些应用中,简单的竖管式潜水位计已被压力计,气压表和电子部件集成的设计所替代。

图 2.6 通过地下水位计的装设确定水力梯度。

潜水位计经常是成组装设的,以便能确定出地下水流的方向。在图2.6(a)中,有三支潜水位计被安装在一个含水的地质构造中。现在我们从图上移去测量设备而仅考虑其测量值[图2.6(b)]。地下水从 h 值较高处流向较低处的,也就是从图的右侧流向左侧的。如果各潜水位计间的距离是已知的,那么水力梯度 dh/dl 就可计算出来。同时,如果该地质结构的水力传导率 K 是已知的,就可用达西定律来计算单位流量(垂直于水流方向任意剖面面积的水流体积)。

有时遇到垂向势能梯度,在这种情况下,就要使用“潜水位计组”来测量。它在同一位置并排安装2个或多个潜水位计(或者可能就在同一个钻孔内),每个潜水位计分别安装于不同的深度,并且可能分布在不同的地质结构中。图2.6(c)和(d) 表示一个垂直向上的地下水流区域中的潜水位计组。

在一个地下水系统内,水力水头分布于整个三维空间中。图2.6中所示的潜水位计组,只能说明在所指示方向上有水流的一个“分量”存在。如果能有大量的潜水位计分布在整个三维水文地质系统内,我们就能够勾划出水力水头相等那些点的等值线(也即等势线—译者注)。在三维空间内,这些点的轨迹形成一个“等势面”。等势面与任何一个二维剖面(不论它是水平的、垂直的或其它状态的)相交所形成的迹线,就是“等势线”。如果在一个剖面内水力水头的形态是已知的,那么就可用垂直等势线勾划出“流线” (沿着最大势能梯度的方向)。等势线与流线交叉所形成的网,就称为“流网”。第5章将对建立这种流网进行详细介绍。第6章将对流网在区域性地下水流解读过程中的实用性进行介绍。

耦合流

现在已有的大量经验和理论证据表明,除水力水头之外,水还可以在其它动力梯度的影响之下在多孔介质中流动。例如,“温度梯度”可以使地下水在没有水力梯度的情况下发生流动 (就像热流一样)(Gurr et al., 1952; Philip & de Vries, 1957)。这一机制在土壤中低温楔(frost wedges)的形成上具有十分重要的意义(Hoekstra, 1966; Harlan, 1973)。

当在土壤中接通电流时,“电力梯度”可以使水从高电压向低电压处发生流动。这一水流流动的机理中涉及水中的带电离子与土壤粘土矿物所带的电荷之间的相互作用(Casagrande, 1952)。这个原理已经以电动力学的方式应用于土壤力学中实现土壤排水(Terzaghi & Peck, 1967)。

化学梯度 甚至可以在没有其它梯度的情况下使水从含盐量高的区域向含盐量低的区域发生流动(也可以使化学成分通过水进行移动)。化学梯度在水流形成方面并不十分重要,但是它对化学物质的运移具有直接影响,并在地下水污染分析中具有很显著的重要性。这些概念在第3,第7和第9章中逐渐介绍。

如果每种梯度都在水流的形成过程中起到一定作用,那么比方程(2.3)更具有普遍性的水流定律可表示为:

v = -L_1\frac{dh}{dl} - L_2\frac{dT}{dl} - L_3\frac{dc}{dl} (2.21)

其中 h 是水力水头,T 是温度,c 是化学浓度;L1L2L3 是比例常数。为了讨论研究,我们令 dc/dl = 0。所剩下的就是:水流在水力水头梯度和温度梯度两者的控制下而发生的情形,即:

v = -L_1\frac{dh}{dl} - L_2\frac{dT}{dl} (2.22)

一般情况下,L1 dh/dl \gg dT/dl

如果在多孔介质中,温度梯度能像控制热流那样促使流体发生流动,那么,水力梯度也就能像造成流体流动那样促成热流,这种共同的相互依赖关系反映的就是我们所熟知的“耦合流”这一热动力学概念。如果我们使 dh/dl = i1dT/dl = i2,那么我们就可以根据方程 (2.22)写出一新的方程:

v_1 = -L_{11}i_1 - L_{12}i_2 (2.23)

v_2 = -L_{21}i_1 - L_{22}i_2 (2.24)

其中 v1 是“流体”穿过介质时的单位流量,v2 是“热”穿过介质时的单位流量。各个L被称为是“唯象系数”。如果方程(2.23)中 L12 = 0,我们就只剩下地下水流的达西定律,而L_11就是水力传导率。如果方程(2.24)中 L21 = 0,我们就只剩下傅里叶(Fourier)定律,此时的 L22 就是热传导系数。

现在,已经可以表达一个完全的耦合方程组。方程组的形式如方程(2.23),还包括方程(2.21)中所有的梯度以及潜在的其它梯度。多孔介质中耦合流理论的发展起源于泰勒和卡里(Taylor & Cary, 1964)的探索。奥尔森(Olson, 1969)曾作了大量的实验研究。贝尔(Bear, 1972)对概念发展的情况作了更为详细的研究。通过热动力学对多孔介质水流物理过程的描述在概念上是非常强有力的,但实际上,在唯象系数 Lij 矩阵中的非对角线特性方面只有很少的几个数据。本书将假设地下水流能够完全通过达西定律[方程(2.3)]来解释,即水力水头[方程(2.18)]连同它的高程和压力能够用来表示总水头。同时水力传导率是方程(2.21)中唯一重要的唯象系数

2.3水力传导率和渗透率

如哈伯特(1956)所指出的,达西定律中被命名为“水力传导率”的比例常数不仅是多孔介质的函数,也是包括水在内的流体特性函数。回顾一下图2.1中的实验设备, 如果在两次操作中 ΔhΔl 保持不变,砂也是相同的砂,但是第一次通入的是水,而第二次通入糖蜜,那么在第二次操作中的单位流量 v 会毫无意外地比第一次要小得多。那么寻找一个与通过的流体无关,只描述多孔介质传导性的参数就很有必要。

作为一种实践,一些实验使用了理想的多孔介质(直径为 d 的均一玻璃珠)。当具有不同密度 ρ 和动力粘度µ的流体,在水力梯度 dh/dl 恒定的条件下通过实验试备时,我们观察到了下述比例关系:

v \propto d^2

v \propto \rho g

v \propto \frac{1}{\mu}

结合达西的原始观察,即 v \propto - dh/dl,上述三个比例关系可推导出达西定律的一个新的变式, 即

v = -\frac{Cd^2\rho g}{\mu}\frac{dh}{dl} (2.25)

参数 C 是另外一种比例常数。对于真实的土壤来说,它必须涵盖与平均粒径无关的其它介质特性对流体的影响,例如:颗粒大小的分布,颗粒的球度和圆度以及它们充填的密实度等。

对比方程(2.25)和原始达西方程[方程(2.3)]可知:

K = \frac{Cd^2\rho g}{\mu} (2.26)

在这个方程中 ρμ 只与流体相关,而 Cd2 则只与介质相关。如果我们定义:

k = Cd^2 (2.27)

那么则有:

K = \frac{k\rho g}{\mu} (2.28)

参数 k 被称为“比渗透率”或“固有渗透率”。如果将 K 称为水力传导率,那我们可以省略限定语根据实际意义将 k 称为“渗透率”,本书的讨论将使用这样的术语。但在较老的教科书和报告中这样的称呼会有些混淆,因为一些旧教科书有时会把水力传导率 K 称为“渗透系数”。

哈伯特(1940)通过微观层面流体通过多孔介质时驱动力与阻抗力之间关系的探索,从基本原理出发,通过方程(2.28)推导得出了方程(2.25)。在进行量纲分析时,他有创见地将常数 g 包含在比例关系中,从而得到方程(2.25)。通过这种方式,将 C 化作无量纲量。

渗透率 k 仅仅是介质的一个参数,并有量纲 [L2]。这个术语在石油工业中被广泛应用, 因为气体、石油和水共存的多相流动系统需要用一个不限于具体流体的传导性参数来描述。当以平方米(m2) 或平方厘米(cm2) 来测量时,K是非常小的。所以石油工程师将一个渗透率单位定义为“达西”(darcy)。把方程(2.28) 代入方程(2.3),则达西定律变为:

v = \frac{-k\rho g}{\mu}\frac{dh}{dl} (2.29)

根据这个方程,1达西被定义为:粘度为1 cp (即1厘泊—译者注)的流体在 ρg dh/dl 项等于1 atm/cm (即1个标准大气压/厘米—译者注) 的压力梯度下单位面积流量为1 cm/s (即1厘米/秒—译者注) 的渗透率。因此1达西就近似地等于 10-8 cm2

在地下水行业中,水力传导率广泛应用gal/day/ft2 作为单位。当达西定律用方程(2.4) 的各项来表示时,这个单位的相关性是显然的,即:

Q = -K\frac{dh}{dl}A

在美国地质调查局最初的定义中,实验室系数和野外系数之间在单位上是不同的。但是最新的定义 (Lohman, 1972) 已经摒弃了这个区别。值得注意的是,野外环境和实验室环境在测量时的温度差别会通过方程(2.28)中的粘度项影响水力传导率值。虽然这种差别一般很小,也很少引入修正系数,但是报告水力传导率的测量是在野外完成还是在实验室内完成还是十分有意义的。因为两者的测量方法截然不同,对于测量值的解释也取决于测量的形式。相比概念上的价值,这一点对于实践更加重要。

表2.2给出了在一个涵盖了很大范围内的多种地质材料在5个不同的单位系统内水力传导率和渗透性的取值范围。这个表部分是以戴维斯(David, 1969) 所总结的数据为基础的。从这些数据中可以得到的基本结论是,水力传导率有一个很大的变化范围。有少数的几个取值超过13个数量级,从实践的角度来说,这个特点意味着确认水力传导率的数量级是很有用的,而相对应数值中第三位以后的数位可能就意义不大了。

表2.3给出了 kK 各种常用单位之间的一组换算系数。比如说,以cm2 为单位的 k 值可以乘以 1.08 × 10–3 转换为以 ft2 为单位的 k值。反之单位从 ft2 转换成 cm2 应乘以 9.29 × 102

在实验室和野外进行水力传导率测量的各种方法详述于8.4至8.6节。

表 2.2 水力传导率和渗透系数的取值范围

表 2.3 渗透系数,水力传导率之间单位的换算系数

渗透系数1 水力传导率
cm2 ft2 darcy m/s ft/s U.S. 加仑/天/平方英尺
cm2 1 1.08 × 10–3 1.01 × 108 9.80 × 102 3.22 × 102 1.85 × 109
ft2 9.29 × 102 1 9.42 × 1010 9.11 × 105 2.99 × 106 1.71 × 1012
darcy 9.87 × 10–9 1.06 × 10–11 1 9.66 × 10–6 3.17 × 10–5 1.82 × 101
m/s 1.02 × 10–3 1.10 × 10–6 1.04 × 105 1 3.28 2.12 × 106
ft/s 3.11 × 10–4 3.35 × 10–7 3.15 × 104 3.05 × 10–1 1 6.46 × 105
U.S. 加仑/天/平方英尺 5.42 × 10–10 5.83 × 10–13 5.49 × 10–2 4.72 × 10–7 1.55 × 10–6 1
1为了获得以平方英尺为单位的 k,把平方厘米表示的 k 值乘以1.0 x 10–3即可.

2.4水力传导率的非均质性与各向异性

水力传导率在地质结构的整个空间内通常是有变化的,同时随着测量方向的不同,在地质结构同一点的水力传导率也会有变化。前一个特性称为“非均质性”,后一个称为“各向异性”。大量野外采样工作中的测量证实了这一常见的特性。对这两种普遍存在的现象进行地质学解释,主要还是依靠我们对造成各种地质环境的地质过程的了解。

均质性与非均质性

如果水力传导率 K 值的大小与地质结构中任一测点的位置无关。那么这个结构就是“均质的”。相反,如果与测点的位置有关,就是“非均质的”,即如果我们在均质结构中设置一个 xyz 坐标系统,则 K(x, y, z) = C, CC 为常数,而在非均质结构中,K(x, y, z) ≠ C

有多少种地质环境,就会有多少种非均质构造型。但值得注意的有三种广泛的类型。图2.7(a)是一个“层状非均质性”实例的垂向剖面,它一般分布在沉积岩和非固结的湖相和海相沉积中。在这里,形成结构的每个单独的地层结构都有一个均一的水力传导率 K1, K2, . . . , 等,但是整个系统却可以被设想为非均质的。有时地层非均质性可能使 K 值在表2.2的所有13个数量级范围内变动 (例如,在粘土和砂的交互沉积中)。同样悬殊的差别也可以在“非连续的非均质性”(由于存在断层,或多层巨厚地层而形成的)情况下产生。也许大多数普遍存在的非连续性地形是表面沉积—基岩接触型。图2.7(b)表示“趋势非均质性”的情况。趋势可能存在于任何类型的地质结构中,但它在形成三角洲,河流冲积扇和冰川沉积平原的沉积过程所形成的地层中特别普遍。就像各类岩石的水力传导率主要决定于节理和裂隙密度一样A, B, C土壤层位经常表示水力传导率的垂向趋势。在大的固结和非固结沉积结构中,趋势性的非均质性可能在几英里之内达到2—3个数量级的梯度。

许多水文地质学家和石油地质学家已用统计分布法提出地质结构中非均质程度的定量描述。现在有大量证据可以支持水力传导率概率密度函数是对数正态分布的说法。沃伦(Warren)和普赖斯(Price, 1961)与本宁(Benion) 和格里菲思(Griffiths, 1966) 也在油田贮油岩石中发现了这一情况。同时威拉德森(Willardson)和赫斯特(Hurst, 1965) 与戴维斯(1969) 支持关于非固结含水结构的结论。水力传导率概率密度函数是对数正态分布的具体情况之一是,K 的对数 Y (Y = log K) 是正态分布的。弗里泽(Freeze, 1975)参考以上资料得出了一张表,Y 上的标准差(不取决于测量单位的)的取值经常是在0.5—1.5范围内。这意味着大多数地质结构中的 K 值的内在非均质性变化在1—2个数量级内。一个地质结构中的趋势性非均质性可以被视为概率分布上“平均值”的一个趋势。在结构中的不同点测量会有相同的标准差,但是趋势的存在意味着结构中整个测量值范围的扩大。

图 2.7 层间非均质和趋势性非均质。

格林科恩(Greenkorn) 和凯斯勒(Kessler, 1969)提出了一种与观测值数理统计学概念相一致的关于非均质性的定义,即如果所有地质结构的空间变化都表现在 K 值中,那么对于一个地质结构来说,以统计学中的趋势来说明非均质性,就是可行的。所以他们把均质结构重新定义为:其水力传导率的“概率密度函数”是单一形态的。也就是说,其非均质性表现在 K 值中,而且这个 K 值在整个空间都保持一个不变的平均值。同样,一个非均质地质结构就被定义为:其概率密度函数是多形态的。也就是说,表现非均质性的、在空间上是有变化的 K 值。 描述一个满足均质性传统定义的多孔介质(就像在直径为 d 的玻璃珠试验中,K 在任何一点都为常数) 用的是“均一性”这个词。如果我们愿意把本节开头所给出的传统定义应用于这一组更合理的概念中,那么我们就能够通过加上形容词“平均”来实现,并把原始定义寓于“平均水力传导率”这一名词之中。

各向同性与各向异性

如果水力传导率 K 与地质结构中某测点上的测量方向无关,那么这个结构在该点上就是“各向同性”的。如果水力传导率 K 随着地质结构中某测点上的测量方向而变化,那么这个结构在该点上就是“各向异性”的。

考虑一个通过各向异性结构的二维垂向剖面,如果我们令结构中一些点位置上 K 值的测量方向与水平线的夹角为 θ,则 K = K(θ),与 K 值达到其最大与最小值时的 θ 角相应的空间方向则称为“各向异性的主方向”。这两个方向永远是相互垂直的。在三维系统中,只要一个平面垂直于主方向之一,那么其它两个主方向就是该平面内最大和最小 K值的方向。

如果建立一个 xyz 坐标系,使坐标的方向与各向异性的主方向相重合,则基本方向中的水力传导率可以被确定为Kx, Ky, Kz。在任一点 (x, y, z),一个各向同性的结构将有:Kx = Ky = Kz,而一个各向异性的结构将有:KxKyKz。如果 KxKyKz (在水平层状的沉积层中这是很普遍的),就把这个结构称为“横向各向同性的”。

为了充分描述地质结构中水力传导率的特性,需要使用两个形容词,一个描述非均质性,一个描述各向异性。例如,对于一个在二维上均质的,各向同性的结构来说,Kx(x, z) = Kz(x, z) = C,对于所有点 (x, z) 来说,C 是一个常数。在一个均质的各向异性系统中,对于所有点 (x, z) 来说,Kx(x, z) = C1,同时对于所有点(x, z) 来说,Kz(x, z) = C2C1C2。图2.8试图进一步阐释四种可能的组合。箭头向量的长度与在 (x1, z1) 和 (x2, z2) 两点处的 Kx 值和 Kz 值是成比例的。

小规模各向异性的根本形成原因是粘土矿物在沉积岩和非固结沉积中的取向。粘土和页岩的岩芯取样表明,水平与垂直的各向异性比值很少超过10:1,并且经常小于3:1。

图 2.8 均质、非均质,各向同性及各向异性可能存在的四种组合。
图 2.9 层间非均质性和各向异性间的关系。

在较大规模上可以表明(Maasland, 1957; Marcus & Evenson, 1961),在层状的非均质性与各向异性之间存在一种关系。图2.9中的地质结构,每一层都是均质的和各向同性的,具有水力传导率 K1, K2, . . . , Kn。我们将表示的这个系统从整体来说像一个简单的均质的各向异性地层;首先,考虑垂直于层面的水流。该系统中的单位面积流量 v 必须是流入量与流出量相等。即,事实上,在整个系统中它必须是个常数。令 Δh1 是穿过第一层时的水头损失,Δh2 是穿过第二层时的水头损失,等等类推,那么总水头损失就是Δh = Δh1 + Δh2 + . . . + Δhn,同时根据达西定律: (方程2.30)

v = \frac{K_1\Delta h_1}{d_1} = \frac{K_2\Delta h_2}{d_2} = ... = \frac{K_n\Delta h_n}{d_n} = \frac{K_z\Delta h}{d} (2.30)

其中 Kz 是这个系统中各层的等价垂向水力传导率。对于 Kz,运用 Δh1, Δh2, . . . , 等的内在含义,扩展方程(2.30)可得:

K_z = \frac{vd}{\Delta h} = \frac{vd}{\Delta h_1 + \Delta h_2 + ... \Delta h_n} =  \frac{vd}{vd_1/K_1+vd_2/K_2+...+vd_n/K_n}

由此可推导出

K_z = \frac{d}{\sum_{i=1}^n\frac{d_i}{K_i}} (2.31)

现在考虑平行于地层的水流。令 Δh 为沿水平距离 l 的水头损失。则通过系统单位厚度的流量 Q 是通过各层流量之和,因此单位面积流量 v = Q/d 由下式给出:

v = \displaystyle\sum_{i=1}^n \frac{K_id_i}{d}\frac{\Delta h}{l} = K_x\frac{\Delta h}{l}

其中 Kx 是一个等价的水平水力传导率。简化后得

K_x = \sum_{i=1}^n\frac{K_id_i}{d} (2.32)

对于一个均质,各向异性的地质结构,方程(2.31)和(2.32)给出了 KxKz 的值,该结构在水力学上等价于一个均质、各向同性层状地质结构,如图2.9。进一步的代数运算可知:对于 K1, K2, . . . , Kn 值所有可能的数组来说,Kx > Kz。 事实上,如果我们考虑一组 K1 = 104K2 = 102 的循环结构 K1, K2, K1, K2, . . . 则 Kx/Kz = 25。对于 K1 = 104K2 = 1 来说,Kx/Kz = 2500。在野外,层状非均质性导致区域各向异性值的数量级达100:1或更大的情况是不常见的。

斯诺(Snow, 1969)指出,裂隙岩也有各向异性的表现,因为它在节理裂隙和空隙上也有方向上的变化。在这种情况下,Kz > Kx 则是很常见的

三维中的达西定律

在可能是各向异性的介质中,对于三维流体来说,需要把前面提出的一维形式的达西定律[方程(2.3)] 加以归纳。在三维中,流速 v 是一个具有 vx, vyvz 分量的一个向量,一个简单和普遍适用的归纳形式是:

v_x =-K_x\frac{\partial h}{\partial x}

v_y =-K_y\frac{\partial h}{\partial y} (2.33)

v_z =-K_z\frac{\partial h}{\partial z}

其中 Kx, KyKz 是在 xyz 方向上的水力传导率,因为 h 现在是 xy 和z的函数, 所以求导必须以偏微分的形式进行。

在本书中,我们将假设这个最简形式是足够描述三维流体。但值得注意的是还可以写出如下的一组更具有普适性的方程形式:

v_x = -K_{xx}\frac{\partial h}{\partial x} - K_{xy}\frac{\partial h}{\partial y} - K_{xz}\frac{\partial h}{\partial z}

v_y = -K_{yx}\frac{\partial h}{\partial x} - K_{yy}\frac{\partial h}{\partial y} - K_{yz}\frac{\partial h}{\partial z} (2.34)

v_x = -K_{zx}\frac{\partial h}{\partial x} - K_{zy}\frac{\partial h}{\partial y} - K_{zz}\frac{\partial h}{\partial z}

这组方程指出了这样一种事实,即:在最一般的情况下,水力传导率实际上有9个分量。如果把这9个分量放在矩阵中,它们就构成一个称为“水力传导率张量”(Bear, 1972) 的二阶对称张量。对于 Kxy = Kyx = Kyz = Kzx = Kzy = 0 这一特殊情况,9个分量减少为3个,这时方程(2.33)就是达西定律在三维情况下的表达。在采用方程(2.34) 而不采用方程 (2.33) 时,所需要的条件是各向异性的主方向与坐标轴xyz一致。在大多数情况下,采用一个满足这个要求的坐标系统是有可能的。但是可以想象的是,在非均质各向异性系统中,不同结构的各向异性有不同的主方向,在这样的系统中选择适当的坐标轴是不可能的。

水力传导率椭圆

在一个均质、各向异性,主方向水力传导率为 KxKzxz介质平面上,我们设想一个任意的流线[图2.10(a)]。

图 2.10 任意水流方向上的单位流量 vs (b) 水力传导率椭圆。

在水流方向有

v_s = - K_s\frac{\partial h}{\partial s} (2.35)

方程中 Ks 是未知的,但可以预测它是在 KxKz 的范围内。我们可以把 vs 分离为 vxvz 两个分量,即

v_x = -K_x\frac{\partial h}{\partial x} = v_s \cos \theta
(2.36)

v_z = -K_z\frac{\partial h}{\partial z} = v_s \sin \theta

此时,由于 h = h(x, z),所以

\frac{\partial h}{\partial s} = \frac{\partial h}{\partial x} \cdot \frac{\partial x}{\partial s} + \frac{\partial h}{\partial z} \cdot \frac{\partial z}{\partial s}(2.37)

从几何学上可知,∂x/∂s = cos θ∂z/∂s = sin θ。把这些关系和方程(2.35)、(2.36) 一起代入方程(2.37)并加以简化,得

\frac{1}{K_s} = \frac{\cos^2 \theta}{K_x} + \frac{\sin^2 \theta}{K_z} (2.38)

这个方程使主方向水力传导率分量 KxKz 在任何角度 θ 的情况下都能够与合成矢量 Ks 联系起来。如果我们把方程 (2.38) 放入直角坐标系,并使 x = r cos θz = r sin θ,则:

\frac{r^2}{K_s} = \frac{x^2}{K_x} + \frac{z^2}{K_z} (2.39)

这是一个具有轴长 \sqrt{K_x}\sqrt{K_z} 的椭圆方程 [图2.10(b)]。在三维系统中,它变为一个具有轴长 \sqrt{K_x}\sqrt{K_y}\sqrt{K_z} 的椭球体。并被称为“水力传导率椭球体”,在图2.10(b)中,如果 KxKz 已知,则在一个各向异性的介质中,水流在任何方向上的水力传导率 Ks 值都可用图解法确定出来。

5.1节将讨论各向异性介质中的流网构图。在这里也要指出,与各向同性介质相反,在各向异性介质中,流线是不与等势线垂直的。

2.5孔隙度与空隙比

如果土壤或岩石单位体积的总量 VT 分为固相部分 Vs和空隙部分 Vv ,则孔隙度可定义为 n = Vv/VT,一般以一个小数或百分数表示。

图2.11表示各种岩石、土壤的纹理结构与孔隙度的关系。其中有必要区分出:由于岩石和土壤基质而形成的原生孔隙度[图2.11(a)、(b)、(c)和(d)]和由于一些次生溶解现象[图2.11(e)] 或在构造上受区域性裂隙[图2.11(f)]控制的“次生孔隙度”。

图 2.11 图 2.11 纹理和孔隙度间的关系。(a) 分选性良好,具有高孔隙度的沉积物; (b) 分选性差,具有低孔隙度的沉积物; (c) 分选性良好,由高孔隙度砾石组成的沉积物,因此该沉积物整体上的孔隙度也较高; (d) 分选性良好的沉积物,但它的孔隙被沉积矿物阻塞; (e) 由于溶解作用而重新生成孔隙的岩石; (f) 由于破碎而再度产生孔隙的岩石 (Meinzer, 1923) 。

表2.4以戴维斯(David, 1969)所收集的数据为基础列出了各种地质材料的代表性孔隙度范围。一般来说,岩石比土壤的孔隙度低;由尖角和圆颗粒组成的砾石、砂和淤泥,其孔隙度比富含扁平状粘土矿物的土壤孔隙度低。分选性不好的沉积物[图2.11(b)] 比分选性好的沉积物 [图2.11(a)] 孔隙度低。 

表 2.4 孔隙度的取值范围

n(%)
非固结沉积物
砾石 25–40
25–50
淤泥 35–50
粘土 40–70
岩石
裂隙玄武岩 5–50
喀斯特石灰岩 5–50
砂岩 5–30
石灰岩,白云岩 0–20
页岩 0–10
裂隙结晶岩 0–10
密室结晶岩 0–5
来源: Davis, 1969.

孔隙度 n 对于水力传导率 K是一个重要的控制性影响因素。从分选性好的沉积沙层或裂隙岩结构中采集的样本得知,具有较高孔隙度 n 的样品,一般具有较高的 K 值。然而,从区域规模上纵观可能遇到的岩石和土壤类型可知,并不总是存在这样的关系。例如富含粘土的土壤,一般比砂或砾石土有更高的孔隙度,但水力传导率却比较低。8.7节将介绍几种通过孔隙度和颗粒大小分析来估算水力传导率的技术。

孔隙度 n 与空隙比 e 有密切关系。空隙比被广泛应用在土壤力学中。e = Vv/Vsne 的关系如下:

e = \frac{n}{1-n} \hspace{1cm} \text{or} \hspace{1cm} n = \frac{e}{1+e} (2.40)

e 通常的取值范围是0—3。

实验室中对土壤样品的孔隙度的测量将在8.4节中谈到。

2.6非饱和流与潜水位

直到这里为止,达西定律、水力水头与水力传导率的概念,都是关于饱和介质的。就是说,其全部空隙都已被水充满。很明显,一些土壤,特别是接近地表的土壤,是很少处于饱和状态的,它们的孔隙通常只部分地被水充填。其余孔隙的空间则是被空气所占有。这种条件下的水流,被称为“非饱和流”或“部分饱和流”。历史上,非饱和流的研究曾是土壤物理学家和农业工程师的研究领域,但是近来土壤学家和地下水水文学家都已达成共识,需要把他们的聪明才智集中起来,发展统一的地下水研究方法,包括饱和与非饱和流。

本节的重点放在非饱和带中水的“液相”迁移水力学的问题上,即不讨论“气相”迁移,也不考虑“土壤水与植物的相互作用”。农业科学对这些课题特别感兴趣,与此同时这些课题对于土壤地球化学的解释占有重要地位。对于在非饱和带中水分转换的物理与化学的详细研究,可以在贝尔(Bear, 1972)等人的著作中看到较基础的介绍,更详细的介绍可查阅柯卡姆和鲍尔(Kirkham & Powers, 1972)以及蔡尔兹(Childs, 1969)的著作。

含水量

如果把土壤或岩石单位体积总量 VT 分为固相体积 Vs,水相体积 Vw 和气相体积 Va 三部分,则“体积含水量” θθ = Vw/VT 来确定。与孔隙度 n 一样,它是一个小数或百分数。对于饱和水流来说,θ = n; 对于非饱和水流来说,θ < n

潜水位

饱和与非饱和条件最简单的水文构造是近地表的“非饱和带”和在深处的“饱和带”[图2.12(a)]组成。我们一般将潜水位定义为两者之间的边界。除此以外,我们还观察到,在潜水位的上方,通常存在一个饱和的毛细边缘带。考虑到这种复杂的情况,我们在对各种饱和—非饱和概念进行定义时,必须注意其使用的一致性。

“潜水位”的最好定义是:在这个水位平面之上,多孔介质孔隙内的流体压力 p 正好是大气压力。这个表面的位置可以通过开放浅井钻凿时的水位来表示,需要注意的是,水井需要有足够的深度以形成积水。如果 p 是用压力计来测量的,那么在潜水位上,则 p = 0。这就意味着 ψ = 0,同时,因为 h = ψ + z,潜水位上任一点的水力水头必须等于该点处潜水位的高程 z。在图例中,我们常用一个小的反转三角形来表示它的位置 [图2.12(a)]。

图 2.12 地表附近的地下水条件。(a) 饱和和非饱和区; (b) 含水量与深度的剖面曲线; (c) 压力水头和水力水头关系;插入的小图:当压力水头低于顶部大气压力且高于底部的大气压力时的持水性; (d) 压力水头与深度的剖面曲线; (e) 水力水头与深度的剖面曲线。

负压力水头和张力计

我们已经知道,在饱和带中,ψ > 0 (通过压力计表示),而在潜水位处,则 ψ = 0。 随之,在非饱和带中,ψ < 0。它反映了这样一个事实:即在非饱和带中,水是在表面张力的作用下滞留在孔隙内的。微观的观察显示,颗粒与颗粒间的每个孔隙通道内分布着凹面的新月形水面[如图2.12(C)]。每个新月形凹面的半径,反映了单个气—水交界面的表面张力。根据滞水的物理机制,土壤物理学家常将该压力水头 ψ (当ψ < 0时)称为”张力水头”或”抽吸水头”。在本书中,立足于一个概念只应有一个术语的定义,所以我们将用”压力水头”来表示正负 ψ

不考虑 ψ 的正负号,水力水头 h 总是等于ψ 和 z 的代数和。但是,在潜水位以上,ψ < 0,所以压力计就不是测量 h 合适的工具了。作为一个替代方案,可以通过 ψ 的测量[使用”张力计”]间接地获得 h。柯卡姆 (Kirkham,1964) 和S.J.理查兹(S. J. Richards,1965)提供了这些设备的详细说明和使用方法。简单地说,张力计由一个孔隙杯与一根密闭的、充满水的软管相连。将孔隙杯插入至土壤在所需要的深度,它与土壤水接触并达到水力平衡。平衡过程中,水从软管中通过孔隙杯进入土壤。在密闭管顶部所形成的真空,可以作为土壤中压力水头的度量。它常通过地表上与软管联接在一起的真空压力计来测量,它所起的作用,可以被想象为把图2.12(c)所示的土壤点1处的压力计倒置过来。为获得水力水头 h,张力计的真空度量值所获得的 ψ 的负值必须加上测量点的高程 z。在图2.12(c)中,点1处的设备是一个张力计;点3处是一个压力计。这是一张示意图。 实践过程中,张力计由一根带度量的软管与一个孔隙杯相连;压力计则是底部带有一个尖端的开口硬管。

非饱和水力参数的特性曲线

在非饱和带中,问题的分析更为复杂。例如土壤含水量 θ 和水力传导率 K 两者都是压力水头 ψ 的函数。含水量与压力水头间存在的联系是显而易见的。由于土壤水分是通过表面张力存在于土壤颗粒之间的,同时这一现象也反映在毛细水每个新月形凹面的半径上。所以我们可以预见,土壤中的较高含水量会使毛细水面具有较大的凹面半径、较低的表面张力和较低的张力水头(即较小的负压力水头)。此外,从实验观察到,θ-ψ关系是滞后的;当土壤正在向潮湿方向变化与正在向干燥方向变化时,它们的关系走势是不同的。图2.13表示在砂质土壤中自然发生的 θψ 之间的滞后作用关系(Liokopoulos, 1965a)。如果这类土壤样品在一个大于零的压力水头下是饱和的,并且该水头逐步(分为若干个阶段)降低,直到显著低于大气压力的水平(即 ψ ≪ 0),那么,对每个阶段的含水量来说,都将表现为“干化曲线” (或称为“排水曲线”) [见图2.13(a)]的形状。反之,如果以较小的时间间隔向干燥土壤中加水,那么压力水头就将沿着返回的路线呈现为”湿化曲线”[或称为“吸入曲线”]的形状。重叠的线称为“扫描曲线”。它表示如果土壤仅仅是部分湿化了,那么在干化时,θψ 就将按照这一路线变化,反之也是一样(即如果土壤仅仅是部分干化了的,那么在湿化时,θψ 也将按照这一路线变化—译注)。

以目前介绍的程度,我们可以预想到,在粗粒土壤中,对于所有 ψ > 0 的情况来说,θ 将等于孔隙度 n。而对于细粒土壤来说,这种关系将保持在一个 ψ > ψa 较大的范围内。其中 ψa 是一个小的负压力水头[图2.13(a)],被称为“空气进入压力水头”。其相应的压力 pa 被称为“空气进入压力”或“气泡压力”。

图2.13(b)表示同一土壤中与压力水头 ψ 有关的水力传导率 K的滞后曲线。对于 ψ > ψa 来说,K = K0,其中 K0 现在称为“饱和水力传导率”。由于 K = K(ψ),θ = θ(ψ),故 K = K(θ)也成立。图2.13(b)的曲线反映了这样一个事实,即:非饱和土壤的水力传导率随着含水量的增大而增大,如果我们写出各向同性土壤中x方向上非饱和水流的达西定律

v_x = -K(\psi)\frac{\partial h}{\partial x} (2.41)

我们将看到,K(ψ) 关系的存在意味着,给定一个恒定的水力梯度,单位面积流速 v 就随着含水量的增加而增大。

图 2.13 一种天然砂质土壤中水力传导率和含水量与压力水头的特征曲线(根据Liakopoulos,1956a)。

在实际上,当含水量增加时,水力梯度保持不变将是不可能的。因为 h = ψ + zθ(ψ),水力水头 h 也受含水量的影响。换句话说,水力水头梯度可表示为压力水头梯度(纯重力作用水头不包括在内),同时这也表示为含水量梯度。图2.12是这三个变量的垂向曲线示意图,表示从地下水从地面向下渗透的情况。因为在图2.12(e)中的水力水头沿纵向下降,因此其水流必然是向下的。较大的 h 正值意味着 |z| \gg |\psi|。换句话说, z = 0 的基准面处于一定的深处。在真实情况中,这三个图象将通过该处土壤的 θ(ψ) 和 K(ψ) 曲线定量的相互关联。例如,如果该土壤的 θ(ψ) 曲线是已知的,同时在野外测量了 θ(z) 曲线,那么 ψ(K) 曲线和 h 曲线就可被计算出来。

如图2.13中所示, θ(ψ) 和 K(ψ) 曲线代表了任意给定土壤的特性。对同一均质土壤的不同样品所进行的一组测量表明,空间上各个分散采样点仅表现常见的统计变化。这些曲线经常被称为“特性曲线”。在饱和带中,我们有两个基本的水力参数 K0n;在非饱和带,这两个参数变为了 K(ψ) 和 θ(ψ) 的函数关系。更加简明的表示如下:

\theta = \theta(\psi) \hspace{1cm} \psi < \psi_a
(2.42)

\theta = n \hspace{1cm} \psi \geq \psi_a

K = K(\psi) \hspace{1cm} \psi < \psi_a
(2.43)

K = K_0 \hspace{1cm} \psi \geq \psi_a

图2.14表示了一些假设的单值特征曲线(即无滞后的),目的是表现土壤特性在曲线形状上的反映。关于对非饱和土壤持水性机理的更全面叙述,读者可参阅怀特(White, 1971)等人的著作。

图 2.14 三种假设土壤的单值特征曲线 (a) 均质砂;(b) 粉砂; (c) 粉砂质粘土。

饱和带、非饱和带和张力饱和带

由于前文经常提及饱和带与非饱和带,因此在这里有必要总结一下饱和带与非饱和带的特性。关于饱和带,可以概括如下:

  1. 它存在潜水位以下;
  2. 土壤孔隙被水所充满,且含水量 θ 等于孔隙度 n
  3. 流体压力 P 高于大气压力,所以压力水头 ψ (通过压力计测量)大于零;
  4. 水力水头 h 必须用压力计来测量;
  5. 水力传导率 K是一个常数,它不是压力水头 ψ 的函数。

关于非饱和带(或称为”包气带”或”渗流带”),概括如下,

  1. 它存在于潜水位以上,并在毛细带以上;
  2. 土壤空隙中仅部分充满了水,含水量 θ 比孔隙度 n 小;
  3. 流体压力 p 比大气压力低;压力水头 ψ 小于零;
  4. 水力水头 h 必须用张力计来测量;
  5. 水力传导率 K和含水量 θ 两者都是压力水头 ψ 的函数。

简而言之,对于饱和水流来说,ψ > ,θ = n,并且 K = K0;对于非饱和水流来说,ψ < 0,θ = θ(ψ) 并且 K = K(ψ)。

毛细带不属于上述任何一带,因为其土壤孔隙虽已被水充满,但压力水头却小于大气压力。现在一个慢慢被接受的表述性名词是”张力饱和带”。在图2.13中可以找到关于它的一个看似反常性质的解释。在特性曲线上存在一个使空气进入的压力水头 ψa < 0,这就是毛细带存在的原因。ψaψ 的一个值,它存在于张力饱和带的顶部,如图2.12(d)中以 ψa 表示的 A 点值。因为粘土质土壤的 ψA 负值比砂的大,所以细颗粒土壤的张力饱和带比粗颗粒土壤的更厚。

一些观点将张力饱和带看作是饱和带的一部分,但是在那种情况下,潜水位就不是两带之间的边界了。从物理学观点上来看,为了建立一套完整水文系统的概念,保留饱和、 张力饱和、非饱和这三个带可能是最恰当的。

对本节前面所讨论的内容中所直接引出的一个观点有必要作一个特别说明。在流体压力小于大气压力的非饱和带或张力饱和带的表面,不会有水流自然向大气一侧流动。水可以在蒸发和蒸腾作用下从非饱和带转移到大气中,但是像河岸泉水或井孔等的自然出流,必须来自于饱和带。5.5节介绍了饱和渗流面的概念,6.5节将着重强调其在山坡水文学中的重要性。

上层滞水和反向水位

到目前为止,我们所考虑的简单水文形态—由一个单一的非饱和带覆盖于一个主要的饱和地下水体之上—是一种常见的情况。这是均质地质沉积延伸到一定深度所表现出的规律。但是另一方面,复杂的地质环境可以造成更为复杂的饱和—非饱和条件,例如,一个存在于高渗透性砂层中的低渗透性粘土层,可以形成一个不连续的饱和透镜体,并且此透镜体的上面和下面都是非饱和的条件。如果我们考虑图 2.15中的 ABCDA 线是一条 ψ=0的等压线, 那么我们就把 ABC 这一部分称为“上层滞水水位”,同时把 ADC 这一部分称为“反向水位”。而 EF 则是真实的水位。

图 2.15 悬滞地下水位ABC,向下凸出的地下水位ADC,真实地下水位EF

饱和条件在时间和空间上都可以是不连续的。大雨可以在地表形成一个暂时的饱和带,这个带的下边界就是处于非饱和带上的反向水位。这种形式的饱和带将在向下渗透和表面蒸发的影响下随着时间而消逝。在第6章我们将更加详细地研究在饱和—非饱和系统中降雨与入渗之间的相互作用。

多相流

本节所概述的非饱和流分析方法,几乎是物理学家最普遍采用的方法之一,但从根本上说,它是一种近似的方法。实际上,非饱和流是多孔介质中“多相流”的 一个特殊情况,它在孔隙通道中有空气和水两相共存。含水量的体积用θw (前文是用 θ表示),θa 代表空气含量体积(与确定 θw 的方法相似),现在要考虑两个流体压力(液相的 pw,和气相的 pa)和两个压力水头 ψwψa。每种土质都有两条流体含量与压力水头的特性曲线,一个是关于水的 θw(ψw) 一个是关于空气的 θa(ψa)。

当涉及到传导率关系时,使用渗透率 k [方程(2.28)] 比使用水力传导率 j 意义更大。因为 K 与流体无关,而 K 与流体有关。流动参数 kwka 被称为介质对于水和空气的“有效渗透率”。每种土质都有两条有效渗透率与压力水头的特性曲线。一个是关于水的 kw(ψw) ; 另一个是对于空气的 ka(ψa)。

用于非饱和流的单相分析法,衍生出了几乎对所有实用目的都有足够精确性的分析技术。 但是也有一些非饱和流的问题,必须考虑到空气和水的多相流动。一些经常遇到的情况包括在湿润锋面前滞留的空气所积累的气体压力,会影响通过土质时湿润锋面的传播速率。威尔逊和鲁辛(Wilson & Luthin, 1963)在实验中遇到过这种影响。扬斯和佩克(Youngs & Peck, 1964)提供了理论上的探讨。同时麦克沃特(Mcwhorter, 1971)进行过全面的分析。正如6.8节将介绍的,空气的滞留也会影响潜水位的浮动,对此,安奇和哈斯克尔(Bianchi & Haskell, 1966)在野外研究了空气的滞留问题。格林(Green, 1970)等人描述了地下流动系统多相分析方法在野外的实际应用。

关于多孔介质中多相流的很多原始研究都是在石油工业中进行的。石油储层工程就涉及到了油、气和水的三相流分析。皮尔逊(Pirson, 1958) 和阿米克斯(Amyx, 1960)等人的著作是野外工作的经典参考文献。斯托尔曼 (Stallman, 1964)将石油多相(流)为地下水水文学做出的贡献做了说明性的综述。

非饱和流的两相分析是“非混相驱替”的一个实例,即,流体彼此驱替而不混合,同时在每个孔隙中存在清晰的流体—流体交界面。如果同时流动的两种流体是彼此可溶的,称为“混相驱替”。在这种情况下,不存在清晰的流体—流体交界面。贝尔(Bear, 1972) 对于多孔介质中混相与非混相驱替两者都提出了进一步的理论解释。在本书中,唯一非混相驱替的实例已在本小节中讨论。在本书的其余部分,非饱和流将用本节第一部分的概念和方法,被当作单相问题来处理。在地下水水文学中,最常发生的混相驱替发生在具有不同水化学成分的两种水(诸如海水和淡水,纯净水与污染水) 的混溶中。与混相驱替相关的迁移过程以及地下水污染分析技术将在第9章中讨论。

2.7含水层和滞水层

在水文学词汇的所有单词中,可能再没有任何一个比“含水层”这一名词在不同语境下拥有更多含义了。对它的理解因人而异,并且在不同时间的理解也不尽相同。它经常被看作是一个单一地层,或一个完整的地质结构,甚至有时是一组地质结构。所以这个词的含义必须永远由其应用时的场合和上下文来确定。

含水层、滞水层和隔水层

一个“含水层”的最佳定义为:在正常的水力梯度下,一个饱和的、能够传输大量水的透水性的地质单元。隔水层被定义为:在正常水力梯度下,不能够传输显著质量的水的饱和地质单元。

在地下水行业广泛使用的另一对定义是:含水层的渗透性较高,并足以使水井生产出有经济意义数量的水,而隔水层则相反。

近些年来,“滞水层”这一名词已经确立,被用来描述地质结构中低透水性地层。这种地层的所能传输的水量,足以影响区域性水流并值得对其进行研究,但是用作供应水井来说,其透水性就不足了。大多数地层都根据含水层或滞水层来分类。极少数地层适合于按隔水层的定义来分类。因此,目前的趋势是只用这些名词中的前两个,而舍弃第三个。

最常见的含水层,是那些水力传导率位于表1.1上半部的地质结构,即:非固结的砂和砾石;透水性的沉积岩,如砂岩和石灰岩;严重破碎了的火山岩和结晶岩。最常见的滞水层是粘土,页岩和致密结晶岩。在第4章中,将对基本的含水层与滞水层类型进行更充分的讨论,但这种讨论只局限在地质结构对地下水分布控制的范围内。

含水层和滞水层的定义,在水力传导率方面是有意未作交待的。这是为了能够在用词的相对性方面留下余地。例如,在砂和粉沙的交叠层序中,粉沙可以被考虑为滞水层,而在粉沙和粘土互层系统中,粉沙就是含水层了。

含水层经常以其地层学名词来称谓。例如,达科塔(Dakota)砂岩,主要是由于迈因策尔(Meinzer, 1923)把它作为含水层并根据其特性所作的评价而得名。其它两个众所周知的北美含水层是伊利诺斯州的圣彼得(St.Peter)砂岩和佛罗里达州的奥卡拉(Ocala)石灰岩。美国主要含水层系统的综述可以在麦吉尼斯(McGuinness, 1963)和马克西(Maxey, 1964)的著作中找到。它们是以迈因策尔(Meinzer, 1923)、托尔曼(Tolman, 1937)和托马斯(Thomas, 1951)早期复杂的资料为基础所定义的。布朗(Brown, 1967)则提供了有关于加拿大主要含水层的资料。

本书许多解释性章节中的理论性分析内容都基于一定的理想假设,其中趋向于把含水层看作是均质、各向同性、且具有不变厚度与简单几何形状的结构。我们希望读者在头脑中明确记住,现实条件与该假设的差别。水文地质学家经常面对的并不是书本上所勾画出的理想情况,而是由非均质、各向异性的结构所构成的复杂的含水层—滞水层系统。有时就像是各种地质过程故意为了使解释和分析更加困难而拼凑起来的。

承压含水层与非承压含水层

“承压含水层”,是指被上、下两个滞水层所封闭起来的含水层。“非承压含水层”或“自由水位含水层”,是一种上部边界为地下水位的含水层。承压含水层处于地层深处,非承压含水层靠近地表(图2.16)。以向下凸出的潜水位为边界的饱和透镜体(图2.15)是非承压含水层的一个特殊情况,这种透镜体有时称为“悬栖含水层”。

图 2.16 U非承压含水层及其地下水位;承压含水层及其潜势面。

在承压含水层内,井中的水位经常会高于含水层顶部的位置,这样的井称为“承压井”,并把这个含水层称为存在于“承压条件下的含水层”。在有些情况下,水位可能会上升到高出地面,这样的井就称为“自流井”,并把这个含水层称为存在于“自流条件下的含水层”。在6.1节中,我们将研究导致自流条件的地貌和地质环境。非承压含水层中水井的水位停留在潜水位处。

潜势面

对于已被供水井广泛开采了的承压含水层来说,业界已经形成了一个并不十分完美但在实用中却根深蒂固的传统概念,即:如果将承压含水层的那些井中的水位高程点绘在一张图上,即含水层中的水力水头的等高线,称为“潜势面”。从一个含水层的潜势面图,可以确定出该含水层中地下水流的方向。

潜势面这一概念,仅仅严格地对于水平含水层中的水平流动有效,而且水平流动也仅只在含水层的水力传导率大大高于其承压层的水力传导率时才能遇到。一些具有潜势面图的水文地质报告,其潜势面是以一组井的水位数据为基础的,而这些底部几乎处于同一高程的水井,通常不会检验它们是否打入了同一个承压含水层内。这种类型的潜势面图,本质上乃是一张从该地区地下所存在的三维水力水头场中所截取的二维水平剖面上的水力水头等高线图。此时,如果存在有垂直方向的流动成分(现实中常常存在),那么根据这种潜势面图所作的计算和解释就会导致明显的错误。

潜势面图也可能被同时存在有承压和非承压两种含水层地区的水位数据所误导,图2.16对这两者作了概念性的区分。一般来说,正如我们将在第6章中看到的那样,这两者是不一致的。

2.8 稳定流与非稳定流

“稳定流”是指水流场中任一点的流速大小和方向不随时间而变化的水流。“非稳定流”(或称“不稳定流”,“瞬变流”)则相反,它表示在水流场中任意位置的流速大小和方向随时间而变化的水流。

图2.17表示通过坝下透水性河流沉积层而流动的稳定地下水流形态(虚线为等势线,实线为流线)。沿着 AB 线,水力水头 hAB = 1000 米。它等于 AB 以上水库蓄水面的高程。同样,hCD = 900米(CD以上电站尾水池水面高程)。通过这个系统的水力水头下降值是100米。如果 AB 上面水库水位和CD上面电站尾水池水位不随时间而变,那么坝下的流网也将不随时间而变化。例如,在点E处的水力水头将是hE = 950米,并将保持为常数。在这种条件下,流速 v = -K \partial h/\partial l 也将随时间而保持为常数。在稳定流系统中不同位置的流速可以是不同的,但在任一给定点,它不随时间而变。

图 2.17 坝下稳态与非稳态地下水流。

我们现在来讨论图2.17(b)中的非稳定流问题,在时间t_0,坝下流网将与图2.17(a)中所示相同。同时 hE 为950米,如果允许水库的水位在 t0t1 这一期间内下降,直到大坝上下游的水位在时间 t1 时完全相等,则坝下流场的最终将达到静态,水将不会由坝的上游向下游流动。在点E,水力水头 hE 将经历一个随时间而下降的过程,即从 t0 时的hE = 950米到最终值的 hE = 900米。直到 t1 之后的某一时间,hE 才可能达到900米,这是由于系统内存在的时间滞后所致。

在稳定流与非稳定流系统中,地下水流的“流线”与“路径”意义不同,“流线”指的是在整个系统内水流的瞬间(在稳定流系统内,是全部时间;在非稳定流系统内,是一个给定时刻)方向。任何时候,它在整个区域性水流中都必须垂直于等势线。而“路径”则指的是在稳定流系统或非稳定流系统水流区域内,某个水质点的具体运动路径。在稳定流系统内,沿入流边界进入系统中的水质点,是沿着如图2.17(a)中所示的流线相一致的路径流向出流边界的,而在非稳定流系统内,其路径和流线就不是一致的。所以勾划流网只能描述某个非稳定流系统内一个给定时刻的具体水流条件,而流网中的流线也只代表这一时刻整个水流系统内水质点的运动方向。同时,流线的形状也是随时间而变化,从这个意义上说,要想确切地画出某个系统内水流的流线,也不是一件容易的事。但是,在地下水污染的研究中,流线的描绘具有显著的重要性。

地下水水文学家,必须懂得稳定流和非稳定流两者的分析技术。本章的最后一节将推导饱和与非饱和两种条件下的水流方程。接下来的章节中提出的实用方法通常都是以这些理论方程为基础的。但对于实践水文学家来说,数学知识的要求也并非十分严格。在地下水水文学中,稳定流基本应用在区域地下水流的分析中。而在水井水力学分析、地下水补给和许多地球化学与土力学的应用中,则需要对非稳定流有所了解。

2.9压缩性与有效应力

非稳定地下水流的分析要引入“压缩性”的概念。它描述的是材料在外加应力作用下所发生的体积变化或应变特性。在弹性材料力学的传统方法中,“弹性模量”是一种更为熟知的材料特性。它是根据应力变化 dσ 与最终应变 的比率来确定的。而压缩系数正是弹性模量简单的倒数,即应变与应力之比—/,而不是应力与应变之比—/。这个专有名词能够用于弹性材料与非弹性材料。对于通过多孔介质的水流来说,必须确定两个压缩系数,一个是水的,另一个是多孔介质的。

水的压缩性

应力是通过流体压力 p 施于流体的。一个压力的增量 dp 导致给定质量水的体积 Vw 产生一定的减小。因此,水的压缩系数 β 按下式确定:

\beta = \frac{-dV_w/V_w}{dp} (2.44)

如果我们要使 β 为一正值,这个负号是必须存在的。

方程(2.44)意味着体积变形 dVw/Vw 与流体应力产生的压力变化 dp 存在一种线性弹性关系。因此弹性系数 β 就是对于水来说受力形变的斜率,并且这个斜率在地下水水文学中所会遇到的流体压力范围内(包括在非饱和带中所遇到的小于大气压力的那些压力)都是不变的。对于经常遇到的地下水温度范围来说,温度对 β 有较小的影响。因此,以最实用的目的来说,我们可以认为 β 是一个常数。β 的量纲是压力或应力量纲的倒数。它的值可以被定为 4.4 × 10–10 m2/N (或Pa–1)。

对于给定质量的水来说,方程(2.44)可以写成如下形式:

\beta = \frac{dp/p}{dp} (2.45)

其中 ρ 是流体密度。方程(2.45)的积分可以获得水的“状态方程”

\rho = \rho_0 \text{exp}[\beta(p - p_0)] (2.46)

其中 ρ0 是在基准压力 p0 时的流体密度。当 p0 为大气压力时,方程(2.46)可根据下式写为标准压力

\rho = \rho_0 e^{\beta p} (2.47)

不可压缩流体指的是 β = 0 和 ρ = ρ0 = 常数的流体。

有效应力

现在我们考虑多孔介质的压缩性。假设一个应力被施加于单位质量的饱和砂上。有三种机理可以使其体积变小:(1)孔隙中的水受压缩;(2)单个砂颗粒受压缩;(3)砂粒重新排列成更为密实的状态。这些机理中的第一个是受流体压缩系数 β 控制的。我们假设第二个机理是忽略不计的,也就是说单个砂质颗粒是不可压缩的,我们的任务就是去确定影响第三个机理的压缩系数项。

这样一来,我们就将利用有效应力的原理。这个概念首先是太沙基(Terzaghi, 1925)提出的,并由斯肯普顿(Skempton, 1961)作过详细分析。大部分土壤力学参考书,诸如太沙基、佩克(1967)和斯科特(1963)等都对其作过充分的讨论。

我们的目的是考虑某饱和结构内,在特定深度的任一平面上的应力平衡(图2.18)。σT 是作用在平面上向下的总应力。它来自于上覆岩石和水的重量,即一部分由多孔介质的颗粒骨架所产生,另一部分由孔隙内水的流体压力 P 所产生。总应力中不是由流体所产生的那一部分被称为有效应力 σe,它是施加于多孔介质颗粒上的应力。土壤颗粒的重新排列与其所形成的颗粒骨架的压缩,是有效应力的变化所引起的,而不是由总应力变化所引起的。两者的关系可以通过下列方程来表示

\sigma_T = \sigma_e + p (2.48)

或以变化项表示为

d\sigma_T = d\sigma_e + dp (2.49)

图 2.18 通过饱和多孔介质任意平面上的总应力、有效应力和流体压力。

许多非稳定的地下水流问题都必须在不发生总应力变化的条件下去分析,也即要求系统中每一点处的上覆岩石和水的重量在整个时间内都基本上保持不变,即

d\sigma_e = -dp (2.50)

在这种条件下,如果流体压力增加,有效应力将以同样程度减小;如果流体压力减小,有效压力将以同样大小增加。在总应力不随时间而变的情况下,系统内任一点处有效应力及其所形成的体积变形,将受该点流体压力的直接影响。因为 p = ρgψψ = h – z (z 是常数,指代所研究点的高程),在某点有效应力的变化,是被该点水力水头的变化所控制的,即:

d\sigma_e = -\rho g \hspace{1mm} d\psi = -\rho g \hspace{1mm} dh (2.51)

多孔介质的压缩性

多孔介质的压缩系数由下式确定:

\alpha = \frac{-dV_T/V_T}{d\sigma_e} (2.52)

其中 VT 是某种土壤的质量总体积。e 是有效应力的变化。已知 VT = Vs + Vv,其中 Vs 是固相的体积,Vv 是被水填充所饱和的空隙体积。有效应力的增量 e 对土壤总体积中产生一个减少量 dVT。在颗粒材料中,该减少量几乎可以全部作为一种颗粒重新排列的结果而表现出来。诚然,单个颗粒本身也是可以压缩的,但其影响通常被认为可以忽略不计。因此,总体上dVT = dVS + dVv;但从我们的研究目的来考虑,我们将假设 dVS = 0,则 dVT = dVv

图2.19(a)表示一个被放入实验室测压器的饱和土壤样品。一个 σT = L/A 的总应力通过活塞施加于样品上,样品的横向被器壁所限制。样品中被截留住的水可以通过活塞上的孔口流到外部,连接到一个压力已知且不变的流体水槽内。土质样品在体积上的减少是在几个不同的 L 值情况下来测量的,而 L 值则是以逐步增大的方式发生变化。在每一步,增加的总应力 T 是由水增加了流体的压力所产生的。但是水从样品向外部水槽的排泄则慢慢地把应力从水转移到颗粒空隙中,这个瞬变过程称为“固结作用”。为使固结过程达到水力平衡,每个 L 值所对应的平衡时间可能很长。但是,它一但达到平衡,在样品以内就会形成 dp = 0。同时从方程(2.49)可知,e = T = dL/A。如果土壤样品有一个原始空隙比 e0 (其中 e = Vv/VT) 和一个原始高度b [图2.19(a)],同时假设 dVT = dVv,则方程(2.52)可以写为:

\alpha = \frac{-db/b}{d\sigma_e} = \frac{-de(1+e_0)}{d\sigma_e} (2.53)

压缩系数 α 通常是用 eσe 所绘出的应变—应力曲线斜率来确定。图2.19(b)中的曲线 AB 是加荷的(σe 增加),BC 是卸荷的(σe 减小)。一般来说,应变—应力关系既不是线性的,也不是弹性的。事实上,对于反复地加荷与卸荷来说,很多细颗粒土质会表现出滞后的特性[图2.19(c)]。土质的压缩系数与流体的压缩系数不同,它不是常数,而是所加应力的一个函数。这一函数也由以前所加负荷的历史来决定。

图 2.19 确定土壤压缩系数的实验用负荷箱(a);有效应力与空隙比之间关系的示意曲线(b)、(c)、(d)。

表2.5 压缩系数的取值范围1

  压缩系数 α (m2/N 或 Pa–1)
粘土 10–3–10–8
10–7–10–9
砾石 10–8–10–10
裂隙岩石 10–8–10–10
坚固岩石 10–9–10–11
水 (β) 4.4 × 10–10
1见附录I表A1.3的换算系数.

图2.19(d)表示粘土和砂的 eσe 曲线比较示意图。砂的曲线斜率较小意味着一个较小的α,同时它的线性意味着 α 值在一个较大的σe范围内都保持为常数。在地下水系统中,随时间而变化的 σe 的升降变化通常是很小的。所以甚至对于粘土来说,通常可以将 α 作为一个常数。表2.5是一个压缩系数表,它给出了各种地质材料的实测数值范围。压缩系数数据的最初来源是多门尼科和米夫林(Domenico & Mifflin, 1955)以及约翰逊(Johnson, 1938)等人的著作。α 的量纲和 β —样,是应力量纲的倒数。其值以SI制的m2/N或Pa–1来表示。应注意,水的压缩系数与弱压缩性地质材料的压缩系数是同一个数量级。

正如图2.19(b)和(c)所表明的,一些土壤在发生膨胀作用时,其膨胀系数与在压缩作用下的压缩系数相比要低得多。对于粘土来说,这两个 α 值的比例经常是10:1左右。对于均一的砂来说,可达1:1。对于膨胀系数远低于压缩系数的土壤来说,由有效应力增加[或者像方程(2.51)所表示的那样,由于水力水头的减小]所引起的体积变形基本上是不可逆的。即当最后去掉有效应力后,它们不能恢复原状。在粘土—砂质含水层—滞水层系统内,粘土滞水层中可能发生剧烈的压实作用(由于巨大的 α 值),这种压实作用是难于恢复的。而在沙质含水层中可能发生的小的形变(由于小的 α 值)却通常具有很大的弹性。

含水层的压缩性

方程(2.53)和图2.18、2.19中的压缩性概念是一维的。在野外场地的深处,如果土壤和岩石仅仅在垂直方向上受到应力时,这种一维概念是有意义的。在任一点处垂直方向总应力 σT 都是由上覆的岩石和水产生的。相邻的材料提供水平方向的约束。垂直方向有效应力 σe 等于 σTp。在这些条件下,含水层压缩系数 α 可以由方程(2.53)的第一个等式来确定。b 在此式中是含水层厚度,而不再是样品的高度。参数 α 是垂直方向的压缩系数。如果用图2.19(a)中所示的实验室设备来测量,土芯采样必须保持垂直,并且负荷必须对任何水平层垂直。在一个含水层内部,α 可能随水平位置而发生变化,也就是说,α 可以是非均一性的,即 α = (x, y)。

必须承认,在大多数的一般分析中,在地层深处存在的应力场不是一维的,而是三维的。在这种情况下,含水层的压缩性必须被看作是一个各向异性的参数。垂向压缩系数 α 就要使用有效应力垂向分量的变化。同时水平压缩性系数就要使用有效应力水平分量的变化。在通过多孔介质的流体流动研究中,三维应力分析概念的应用是一项高级课题而无法在这里详尽说明。幸运的是,对于许多实际情况来说,水平应力场中的变化是非常小的。并且在大多数分析中,它可以被假设为忽略不计的。把含水层压缩性系数 α 看作是一个各向同性的参数,对于我们的目的来说已经是足够的。但是必须记住,它实际上是垂直方向上的压缩系数,并且这是所预期的有效应力变化的唯一最大的方向。

为说明发生在可压缩性含水层中形变的性质,图2.20中展示了一个厚度为 b 的含水层。如果上覆材料的重量保持不变,并且含水层中的水力水头减少 -dh,有效应力 e 的增加可以根据方程(2.51)并由 ρg dh 得出。则根据方程(2.53),含水层的压实量为:

db = -\alpha b \hspace{1mm} d\sigma_e = -\alpha b \hspace{1mm} \rho g \hspace{1mm} dh (2.54)

其中的负号表明:水头减小,使含水层厚度 b 也减小。

图 2.20 由于抽取地下水而导致的含水层压。

使含水层中水力水头下降的途径之一,是从井中抽水。抽水使含水层中产生朝向水井的水平水力梯度,并导致靠近水井的每一点的水力水头的降低。相应地,在这些点处有效应力也増加了,并且导致含水层的压实。反过来,向含水层中注水,提高水力水头,减小有效应力,将使含水层膨胀,如果含水层—滞水层系统由于抽取地下水而导致的压实发展到地表面,其后果就是“地面沉陷”。在8.12节中,将详细讨论这一问题。

非饱和带中的有效应力

方程(2.51)中的第一个等式指出,有效应力 σe 和压力水头 ψ 之间应是线性关系。这个关系在饱和带中是成立的,图2.18所示的概念也是基于该关系。但有很多证据表明该关系在非饱和带中并不成立(Narasimhan, 1975)。对于非饱和流,毕晓普和布莱特(Bishop & Blight, 1963)认为方程(2.51)应修正为

d\sigma_e = -\rho g\chi \hspace{1mm} d\psi (2.55)

其中参数 χ 由饱和度、土质结构和土壤的湿化—干化历史共同决定。图2.21中的曲线ABC就表示了这样一种关系。对于ψ > 0 来说,χ = 1,对于ψ < 0, 来说,χ ≤ 1;对于 ψ ≪ 0 来说,χ = 0。

图 2.21 在饱和与非饱和带中有效应力与压力水头之间的关系(根据Narasimhan,1975)。

χ 方法是一种经验性的方法,与此同时,它的使用反映了一个事实,即在非饱和水流场中,对于小于大气压力的流体压力对总应力的贡献尚未得到充分了解。作为一个近似处理,如图2.21中曲线 ABD 所设想的那样,我们将其假设为没有相关影响也是合理的。在这样的假设下,对于 ψ < 0, 來说,χ = 0, e = T,同时非饱和带中压力水头(或湿度)的变化也不会导致有效应力的变化。

非饱和带中多孔介质压缩系数的确定,也正如它在饱和带中的那样,是由方程(2.52)确定的,但是压力水头对有效应力的影响现在则被认为不存在或可忽略不计。

2.10导水系数与贮水系数

为了描述饱和地下水流的水力状态,有六项流体与多孔介质的基本物理特性必须进行了解。这六项特征目前都已经得到介绍,即:水的密度 ρ,粘度 μ,和可压缩性 β;以及多孔介质的孔隙度 n(或空隙比 e),渗透率 k,和可压缩性 α。 所有其它用于描述地质结构的水文地质特性参数都可由这六项导出。例如,从方程(2.28)可以看出,饱和水力传导率 K是k、ρ 和μ的组合,在这一节中,我们将研究单位贮水量 Ss,贮水系数 S 和导水系数 T 等概念。

单位贮水量

某饱和含水层的“单位贮水量” Ss,是指该含水层水力水头下降一个单位时,单位体积含水层释放出来的水的体积。我们从2.9节中了解到,水力水头 h 的下降,会引起流体压力 p的下降和有效应力 σe 的增大。由于 h 下降从贮量中释放出的水,由两种机理产生的,即(1)因 σe 的增大而引起含水层的压实;(2)由于 P 减小而引发的水的膨胀。第一种机理是由含水层的压缩系数 α 来控制,而第二种机理是由流体的压缩系数 β 控制的。

我们首先考虑由含水层的压实而产生的水。这种在压实期间从单位体积的含水层中被挤压而释放的水的体积,等于单位体积含水层减小的体积。这种体积上的减小量 dVT 是负值,但是所产生的水量 dVw 是正值。所以,从方程(2.52)得到

dV_W = -dV_T = \alpha V_Td\sigma_e (2.56)

对于一个单位体积 VT = 1,根据方程(2.51),即 e = –ρg dh。对于一个单位水力水头的下降,dh = –1,同时可有

dV_W = \alpha \rho g (2.57)

现在来探讨由于水的膨胀所产生的体积。由方程(2.44)可知

dV_W = -\beta V_W dp (2.58)

在总单位体积 VT 中水的体积 Vw 等于 nVTn 是孔隙度。如果 VT = 1,dp = ρg dψ = ρg d(h – z) = ρg dh,当 dh = –1 时,方程 (2.58) 变为:

dV_W =\beta n \rho g (2.59)

单位贮水量 Ss 是方程(2.57 )和(2.59)两项之和:

S_s = \rho g (\alpha + n\beta) (2.60)

通过这个方程的量纲可知,Ss 具有一个独特的量纲 [L]-1。这也与 Ss 的定义相一致,即:单位体积含水层中单位水头下降时所产生的水体积。

承压含水层的导水系数与贮水系数

对于一个厚度为 b 的承压含水层来说,其导水系数(或导水能力)T 表示为:

T = Kb (2.61)

其贮水系数(或贮水率)S 由以下方程确定:

S = S_sb (2.62)

如果我们把(2.60)代入方程(2.62),S 的定义可扩展为

S = \rho gb(\alpha + n\beta) (2.63)

厚度为b的饱和承压含水层的贮水系数可以定义为:当垂直于含水层表面的水力水头下降一个单位时,从单位含水层表面面积上释放出来的贮水体积量。承压含水层的水力水头通常以潜势面形式的方程出现,图2.22(a)以此诠释了贮水系数的概念。

图 2.22 在承压含水层(a)和非承压含水层(b)中贮水性质的示意图(Ferris et al., 1962)。

其中水力传导率 K 的量纲是 [L/T],从方程(2.61)可以清楚地知道,导水系数 T 的量纲是 [L2/T]。SI制单位是m2/s。T和 S 是北美地下水行业中广泛使用的术语,经常用英尺·磅·秒(FPS)这一工程单位来表示。如果 K 以加仑/天/平方英尺来表示,那么 T 的单位就是加仑/天/英尺。T 的取值范围可以用合理的含水层厚度 (5—100米) 乘以表2.2中适当的 K 值获得。从地下水开发的角度来说,导水系数大于0.015平方米/秒 (或0.16平方英尺/秒,或10万加仑/天/英尺)的含水层较适于开发。贮水系数是无量纲量。在承压含水层中,其数值范围在0.005到0.00005之间。结合 S 的定义及其取值范围,可以清楚地了解到,要使承压含水层生产出相当的出水量,在较大的面积上会有显著的水头变化。

导水系数和贮水系数在滞水层和含水层中会被详细说明。在大多数应用中, 滞水层中纵向水力传导率比其导水系数具有更重要的意义。应注意到,在粘土滞水层中,αβ,贮水系数[方程(2.63)]和单位贮水量定义[方程(2.60)]中的 项可忽略不计。

可以把导水系数 TK 分别与贮水系数 SSs 结合起来,以确定某个结构的参数。这个参数,即“水力扩散系数” D 可表示为:

D = \frac{T}{S}=\frac{K}{S_s} (2.64)

这一术语目前的应用还不够广泛。

导水系数 T 和贮水系数 S 的概念,最初是从承压含水层中水井水力学分析中发展起来的。当地下水在厚度为 b 的承压含水层中向出水井进行二维水平流动时,各个术语都有明确的定义,但是,在很多其它地下水应用中它们就失去了意义。如果一个地下水问题具有三维的特征,最好使用水力传导率 K 和单位贮水量 Ss;或者最好去使用渗透率 k、孔隙度 n 和压缩系数 α 等基本参数。

非承压含水层的导水系数与单位产水量

在非承压含水层中,导水系数没有像承压含水层中那样的明确定义。但它仍可以被使用,并由同一个方程(2.61)来确定,但此时的 b表示含水层的饱和厚度或者下伏滞水层(约束含水层的)顶部以上潜水位的高度。

非承压含水层的贮水量这一术语以“单位产水量” Sy 被广泛熟知。它定义为:含水层单位表面面积在潜水位下降一个单位情况下,从贮存的水量中释放出来的水的体积。它有时被称为非承压贮水系数。图2.22(b)对它的概念进行了诠释。

单位产水量这一概念,清晰地解释了它所代表的饱和—非饱和相互作用这一问题。图2.23表示在非饱和带中,t1t2 两个时间的潜水位位置变化,以及随深度而变化的含水量纵向曲线图。斜线部分的面积代表从单位剖面柱状体的贮水量中释放出的水的体积。如果潜水位的下降代表单位降低值,那么阴影的面积就代表单位产水量。

图 2.23 从地下水位上方非饱和含水量剖面曲线来表示单位产水量概念。

非承压含水层的单位产水量比承压含水层的贮水系数大得多。Sy 通常的取值范围是0.01—0.30。较高值代表从非承压含水层的贮量中释放的水,它表示的时土壤孔隙的实际排水;而从承压含水层贮量中的释水,则是由于流体压力变化所引起的水膨胀和含水层压实等外部效应所产生的副效应。另外,非承压含水层易于循环补给的特性使它有利于通过水井来进行开采,而且与承压含水层相比较,获得同样的产水量只需要在较小的面积范围内,以较小的水头变化就可以实现。

非饱和带贮水量

在非饱和土质中,含水量 θ 的变化,如图2.23中所示,与压力水头 ψ 的同时发生变化,并通过图2.13(a)的特征曲线关系 θ(ψ) 来表示。特征曲线的斜率代表了土壤的非饱和贮水特性,称为“单位水分容量 C ”,并由下方程确定:

C = \frac{d\theta}{d\psi} (2.65a)

压力水头的增量 (譬如说,在图2.13中由 —200厘米增到 —100厘米),必然与非饱和土壤水分贮量的增量 dθ 相联系。由于 θ(ψ) 是非线性的,并且是滞后的,所以 C 也具有这些特征。C 不是一个常数;而是压力水头的一个函数,即 ψ: C = C(ψ)。事实上,在饱和带中,对于所有 ψa 来说,含水量 θ 都等于孔隙度 n(一个常数),所以 C = 0。对于C来说,与方程(2.42)相似的一个方程是

C = C(\psi) \hspace{1cm} \psi < \psi_a
(2.65b)

C = 0 \hspace{1cm} \psi \geq \psi_a

非饱和土壤的输水与贮水特性可以通过特征曲线 K(ψ) 以及 θ(ψ) 或 C(ψ) 所完全确定下来。与(2.64)中的方程相似,“土壤—水扩散系数”(Soil—water diffusivity)可被定义为

D(\psi) = \frac{K(\psi)}{C(\psi)} (2.66)

2.11地下水流方程

几乎在基础科学和工程学的每一个领域内,分析技术的基础都是对物理过程的理解。在大多数情况下,可以通过数学去描述这些过程,地下水流也不例外。地下水流动的基本定律是达西定律,同时,当它与描述流体流经多孔介质时的质量守恒连续性方程结合在一起时,其结果是一个偏微分方程。在本节中,我们将简要介绍(1)稳定饱和流;(2)瞬变饱和流;(3)瞬变非饱和流等方程的发展情况。对于这三种水流的方程,数学家们都非常熟悉。它们的演算技术广泛而有效,在科学和工程上的应用也很常见。一般情况下,水流方程出现在“边界值问题”中。所以在本节的最后,我们将对这个概念进行探讨。

在地下水水文学中,很多标准的分析技术都是以包括偏微分方程在内的边界值问题为基础。当我们着手学习各种技术时,先对这些方程的进行基本的了解是非常有益的。所幸的是,它不是一项绝对的要求。在大多数情况下,这些技术可以被解释和理解,而不需要回到数学每一步的基本推导中去。进行研究工作的水文地质学家必须每天与水流方程打交道;而参与实践的水文地质学家,如果愿意的话,可以避免使用高等数学。

稳定饱和流

以单位体积的多孔介质为例,如图2.24所示。这样的单元经常称为“控制性单元体积”。穿过饱和多孔介质的稳定流物质守恒定律,要求流入任一控制性单元体积的质量流速率等于流出这一控制性单元体积的质量流速率。把这一定律转换为数学方程的“连续性方程”可以表示为 (图2.24)

-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = 0 (2.67)

图 2.24 流经多孔介质的控制单元体积。

通过对 项的量纲检验可知,它们具有穿过控制性单元体积上单位横截面质量流速率的量纲。如果该流体是不可压缩的,则 ρ(x, y, z) = 常数,并且各个 ρ 都可从方程 (2.67) 中消去。 即使流体是可压缩的,即 ρ(x, y, z) ≠ 常数,也能看出,ρ ∂vx/∂x 类型的各项方程比 vx ∂ρ/∂x 类型的各项方程大得多。当展开方程 (2.67) 时,两者都会出现。在任何一种情况下,方程(2.67)都简化为

-\frac{\partial v_x}{\partial x} -\frac{\partial v_y}{\partial y} -\frac{\partial v_z}{\partial z} = 0 (2.68)

将达西定律代入方程(2.68),就可产生下列各向异性饱和多孔介质中稳态流的水流方程:

\frac{\partial}{\partial x}\left(K_x \frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(K_y \frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z}\left(K_z \frac{\partial h}{\partial z}\right) = 0 (2.69)

对于各向同性介质来说,Kx = Ky = Kz,如果介质也是均质的,则 K(x, y, z) = 常数。方程(2.69) 就简化为均质、各向同性介质中稳定流的水流方程:

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = 0 (2.70)

对于数学家们来说,方程(2.70)是最基本的几种偏微分方程之一的 “拉普拉斯方程”。方程的解 h(x, y, z) 描述了在任意一个三维空间中一点的水头 h。方程(2.70)的解让我们得以绘制 h 的等势线图,同时可以绘制流线和流网图。

对于二维水流场中(即 xz 平面内)的稳定饱和流来说,方程(2.70)的中间一项可以消去,得到解h (x, z)。

瞬变饱和流

饱和多孔介质中瞬变流的质量守恒定律要求进入任一控制性单元体积内的“净”流入率,等于单位时间内该单元内所贮流体质量的变化率。参考图2.24,其连续性方程将表示为:

-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = \frac{\partial(\rho n)}{\partial t} (2.71)

或将右边的表达进行展开,

-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = n\frac{\partial \rho}{\partial t} + \rho\frac{\partial n}{\partial t} (2.72)

方程(2.72)右边的第一项表示在密度 ρ 的变化下,由水的膨胀所产生的水的质量变化率。第二项表示由孔隙度 n 的变化所反映的水的质量变化率。第一项是受流体压缩系数 β 所控制,第二项受含水层压缩系数 α 控制。简化方程(2.72)右侧两项所需要的分析我们已经在2.10节中完成。我们知道,ρn 的变化,两者都是由水力水头 h 产生的。同时由于水头单位的降低,从两种机理中所产生的水的体积是 Ss。即,由 Ss = ρg(α + ) 所给出的得贮水量。故所产生的水的质量变化率(单位时间流体质量贮量的变化率)是 ρSS ∂h/∂t,同时方程(2.72)变为

-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = \rho S_s\frac{\partial h}{\partial t} (2.73)

应用链式法则展开左侧各项后,由于 ρ ∂vx/∂x 形式的各项远大于 vx∂ρ/∂x 形式的各项,因此可以从方程(2.73)的两侧消去 ρ。代入达西定律后得到

\frac{\partial}{\partial x}\left(K_x \frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(K_y \frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z}\left(K_z \frac{\partial h}{\partial z}\right) = S_s\frac{\partial h}{\partial t} (2.74)

这是在各向异性介质中,瞬变饱和流的方程。如果介质是均质、各向同性的,方程(2.74)就可以简化为

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = \frac{S_s}{K}\frac{\partial h}{\partial t} (2.75)

或将 Ss 展开,得到

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} + \frac{\partial^2h}{\partial z^2} = \frac{\rho g(\alpha+n\beta)}{K}\frac{\partial h}{\partial t} (2.76)

方程(2.76)是我们熟知的“扩散方程”。其解 h (x, y, z, t) 描述了某水流场内任一时刻在任一点处的水力水头值。求解时,需要知道 Kαn 等三个基本水文地质参数以及 ρβ 等流体参数。

对于厚度为 b 的水平、承压含水层来说,S = SsbT = Kb,方程(2.75)的二维方程变为

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2}  = \frac{S}{T}\frac{\partial h}{\partial t} (2.77)

其解 h(x, y, t) 描述了任一时刻通过水平含水层平面上任一点处的水力水头场。求解时,需要知道含水层参数 ST

瞬变饱和流 [通过方程(2.77)所给出的任一形方程(2.74)] 的水流方程的发展得益于达西(1856)所建立的流动方程,哈伯特(Hubbert, 1940)关于水力势能的阐述,迈因策尔(Meinzer, 1923)对含水层弹性概念的认识和太沙基(Terzaght, 1925)对有效应力问题的研究等。这一经典问题的发展,首先是由雅各布(Jacob,1940)所提出,在他1950 发表的论文那里可以找到最完整的方程。本节中所描述的发展情况与前几节所提到的贮水量这一概念,最初都来自雅各布。

最近几年,这一经典问题的研究出现了新的发展。拜欧特(Biot, 1955)认为,在压缩性含水层中,达西定律必须改用针对颗粒的相对流速来表示,同时库珀 (Cooper, 1966)指出,在变形的介质中,选取固定不变的控制性单元体积是不妥当的。库珀认为, 如果把流速看作是相对的,并且与其对应的系统也在发生形变,那么雅各布的经典理论是正确的。同时他也提出,德威斯特(De Wiest’s, 1966)对这个问题所作的说明并不正确(Davis & De Wiest, 1966)。关于雅各布—库珀研究的详细内容参见附录II。

使用了垂向含水层压缩性概念所进行的经典研究,基于一个假设:在压缩性含水层中,应力应变仅发生在垂直方向。因此,该方法与三维流场和一维应力场相耦合。更通用的三维流场与三维应力场耦合的方法首先由拜欧特(1941,1955)提出。沃瑞特(Verruijt, 1969)对这个方法进行了更完善的总结。

从实践的角度来说,不需要去考虑相对流速、变形的坐标系统或三维应力场。本节中所提出的经典水流方程方程已经足够。

瞬变非饱和流

我们将饱和度θ’ 定义为 θ’ = θ/n。其中 θ 是含水量,n 是孔隙度。对于非饱和的控制性单元体积中的水流来说,连续性方程必须揭示出含水量的单位时间变化率,以及由于水的膨胀和含水层承压而产生的贮水量的单位时间变化率。方程(2.71)中的 ρn 项必须变为 ρnθ’,与此同时,方程(2.72)变为

-\frac{\partial(\rho v_x)}{\partial x} -\frac{\partial(\rho v_y)}{\partial y} -\frac{\partial(\rho v_z)}{\partial z} = n\theta' \frac{\partial \rho}{\partial t} + \rho\theta' \frac{\partial n}{\partial t} + n\rho \frac{\partial \theta '}{\partial t} (2.78)

对于非饱和流来说,方程(2.78)右侧的前两项比第三项小得多。舍去这两项,并从两侧同时消去各个 ρ,代入非饱和形式的达西定律 [方程(2.41)],且 n dθ’ = ,得:

\frac{\partial}{\partial x}\left[K(\psi)\frac{\partial h}{\partial x}\right] + \frac{\partial}{\partial y}\left[K(\psi)\frac{\partial h}{\partial y}\right] + \frac{\partial}{\partial z}\left[K(\psi)\frac{\partial h}{\partial z}\right] = \frac{\partial \theta}{\partial t}
(2.79)

一般经常把方程(2.79)表示为 θ 或者 ψ 为自变量形式的方程。对于后者情况来说, 需要把右侧的分子分母都乘以 。然后,再根据单位含水量 C 的定义 [方程: (2.65)],且 h = ψ + z,得到:

\frac{\partial}{\partial x}\left[K(\psi)\frac{\partial \psi}{\partial x}\right] + \frac{\partial}{\partial y}\left[K(\psi)\frac{\partial \psi}{\partial y}\right] + \frac{\partial}{\partial z}\left[K(\psi)(\frac{\partial \psi}{\partial z}+1)\right] = C(\psi)\frac{\partial \psi}{\partial t}
(2.80)

方程(2.80)是在非饱和介质中,关于 ψ 的瞬时流方程。通常称为 “理查德方程(Richards equation)”,以纪念首先提出它的土壤物理学家(Richards, 1931)。它的解 ψ (x, y, z, t) 描述了某水流场中在任一时刻任一空间点的压力水头场。通过 h = ψ + z 这一关系,它能够被很容易地转换为一个水力水头解 h (x, y, z, t)。这一解需要首先熟知特征曲线 K(ψ) 和 C(ψ) 或 θ(ψ)。

弗里泽(Freeze, 1971a)和那拉西姆汉(Narasimhan, 1975)曾试图把非饱和流方程 [方程(2.80] 和饱和流方程 [方程(2.74)] 耦合起来。而欲改进寓于饱和—非饱和系统中的理论,还有待于对非饱和带中有效应力的原理进行更深入的理解。

边界值问题

边界值问题是一种“数学模型”问题。与其相关的分析技术由四个步骤组成,包括:(1)物理问题的研究;(2)以数学的形式来表达物理问题;(3)通过合理的数学方法求解;(4)通过物理问题的术语解释求解的结果。以流体物理学为基础的数学模型方程,通常采用势能场理论研究者所倡导的边界值问题形式的方程,如同它在物理学中应用的例子,固体中的传导问题(Carslaw & Jaeger,1959)。

为了准确定义地下瞬变流的边界值问题,我们需要知道(1)流动区域的大小和形状(2)区域的水流方程;(3)围绕区域的边界条件;(4)区域的初始条件;(5)控制水流的水文地质参数的空间分布;(6)求解的数学方法。如果这个边界值问题是针对稳定流系统的,则不需要条件(4)。

以图2.25(a) 中简单的地下水流问题为例。区域 ABCD 所包含的是一个均质、各向同性的多孔介质,其水力传导率为 K1,边界 ABCD 是不透水的。在 ADBC 上的水力水头分别是 h0h1。假设流场达到稳态,且 h0 = 100米,h1 = 0米,通过计算可知,点 E 处的水力水头是50米。显然,我们已经应用了上面所列出的(1)、(3)和(5)等条件;通过估算求解方法(6)也可以实现对结果的检查,即便对区域内的水流方程并不清楚。如果要求解图 2.25(b)中那样更为复杂的问题(坐落在一个斜坡上的填土坝),在点 F 处的水力水头值就不会这么容易地得到。这就要求我们对水流方程有清楚的认识,以获得相应的数学解。

图 2.25xy 平面内的两个稳态边界值问题。

求解方法可以被粗略地分为五个类别,即(1)目测求解,(2)绘图求解,(3)相似模型求解,(4)解析解,(5)数值解。我们在前面已看到了一个目测求解的实例。在第5章提出的流网构图法可以看作是边界值问题的绘图求解法。类似的电场模型将在5.2节 和8.9节中讨论。数值解是现代计算机模拟技术的基础,将在第5.3和8.8节中讨论。

求解边界值问题的最直接了当的方法是解析解。后面介绍的很多地下水分析技术都是以解析解为基础的,可以通过一个简单的实例来进行讨论。以图 2.25(a)中的边界值问题为例,其解析解为:

h(x, y) = h_0 - (h_0 - h_1) \frac{x}{x_L} (2.81)

该方程描述是一组分布在 ABCD 范围内,并平行于边界 ADBC 的一组等势线。由于这些等势线平行于 y 轴,由于 y 没有出现在方程(2.81)的右侧,所以 h 并不是一个 自变量为 y 的方程。在点E, x/xL 0.5 的位置,且 h0 = 100米,h1 = 0米,那么从方程(2.81)可知 hE = 50米,与估算的值相符。在附录II中,为获得方程(2.81)的解析解,应用了“分离变量法”,该解满足水流方程方程和边界条件的。

2.12达西方法的局限性

达西定律几乎为所有水文地质环境中地下水的流动提供了一种正确的解释。总的来说,达西定律适用于:(1)饱和流与非饱和流 (2)稳定流与瞬变流 (3)含水层中的水流与滞水层中的水流, (4)均质系统和非均质系统中的水流, (5)各向同性介质和各向异性介质中的水流, (6)岩石介质和颗粒介质中的水流。 在本书中,我们将达西定律作为我们定量分析的一个有效基础。

正因如此,研究达西方法在理论与实践上的局限性是极为重要的。这就需要看一下寓于我们关于连续体的定义中的假设;考察一下微观和宏观水流的概念;研究一下达西定律的上限和下限;考虑一下裂隙岩水流的特殊问题。

达西连续体和代表性单元体积

2.1节中已提到,达西定律要求用一个代表性的连续体以取代组成多孔介质颗粒的真实集合体。这个连续体的方法适用于宏观而不是微观层面。如果达西定律是一个宏观定律,那么就必然有—个定义达西定律有效性的多孔介质单元大小的下限。哈伯特(Hubbart, 1940)曾提出过这个问题,他借助于图2.26定义了“宏观”这个词。图2.26是多孔介质孔隙度的假定值,在多孔介质中某点 P 进行采样,体积按 V1, V2, . . . , 不断增加。贝尔(Bear, 1972)把图2.26中的体积 V3 定义为“代表性单元体积”。他提到这个体积必须比单个孔隙大。事实上,它必须包含有足够数量的孔隙,这样连续体方法中所要求的统计平均值才有意义。小于这个体积,就没有一个可以代表 P 点孔隙度的值。本书通篇提到的孔隙度、水力传导率和压缩系数都是在一个比代表性单元体积大的样品中进行测量得到的。实践过程中,它们通常是常见大小的岩芯土质样品中测出的数值。当分析的规模涉及到体积,如图2.26中的 V5,非均质介质有可能包含一个以上的地层,这种规模有时被称为“大规模(megascopic)”。

图 2.26 微观与宏观领域及其代表性单元体积 V3(Hubbert, 1956; Bear, 1972)。

2.11节中的每个水流方程式都引用了达西定律。因此,必须承认的是,以边界值问题为基础并且包含这些方程的分析方法,以达西连续体的方式被应用在宏观规模上。有一些地下水现象,诸如通过多孔介质中示踪剂的运动,是不能在这种规模下分析的。因此,很有必要研究一下宏观达西连续体定义下的达西流速(或单位面积流量)与多孔介质液相中微观流速之间的相互关系。

单位面积流量、宏观流速和微观流速

如果我们像贝尔(Bear, 1972)一样区别“体积孔隙度”(由2.5节所定义的)和“面积孔隙度” nA(可以通过某单位体积的任一横剖面根据 nA = Av/AT 确定,其中 Av 是孔隙所占的面积,AT 是总面积),那么我们的推导将会更严谨。如图2.27(a)所示,在一个给定的单位体积内,不同的横剖面有不同的面积孔隙度 nA1, nA2, . . . 体积孔隙度 n 是各种可能的面积孔隙度 nAi 的一个平均值。

图 2.27 面积孔隙度(a)与平均线性流速(b)。

对于任何横剖面 A 来说,其“单位面积流量” v 由方程(2.1)确定

v = \frac{Q}{A}

单位面积流量等于体积通量 Q 除以被整个横剖面面积(孔隙和基质)。这个流速是适于宏观连续体方法的。实际上,水流只是通过横剖面中孔隙所占据的部分。对于横剖面 A 来说,我们可以确定出一个流速 \bar{v}_1 = Q/n_{A1}A 代表被实际横剖面面积(水实际流过的面积)除得的体积流量,对于不同的横剖面 A1, A2, . . . 来说,我们可以确定出不同的\bar{v}_1, \bar{v}_2, . . . . ,如果我们以 \bar{v} 来表示它们的平均值,则

\bar{v} = \frac{Q}{nA} = \frac{v}{n} = -\frac{-K}{n}\frac{\partial h}{\partial l}
(2.82)

这个流速 \bar{v} 有许多其他的名称。我们将认为它是“平均线性流速”。其中 Q, n, A 以及 \bar{v} 是可以测量的宏观项。必须强调的是,\bar{v} 并不代表穿行于孔隙空间的水分子的平均流速,那些真正的微观流速一般比 \bar{v} 要大。因为水分子必须沿着比 \bar{v} 所代表的线性化的路径更长些的不规则路径来运动,这一点可通过图2.27(b)表示。所幸的是,对孔隙通道中真正的微观流速的研究兴趣并不大,因为它们很难被确定。在本书中所涉及的所有情况,达西流速 v 和平均线性流速 \bar{v} 已经足够。

作为进一步解释 \bar{v} 的基础,我们可以考虑一个试验,即用示踪剂来确定大量地下水沿着水流路径移动一个虽然较短,但也足够长的距离 AB 所需要的时间。v 由移动距离与时间之比来确定。移动距离为从 AB 的线性距离。移动时间是示踪剂从 AB 所需的时间。对于 \bar{v} 的定义,纳尔逊(Nelson, 1968)曾提出一个与上式(2.82)稍有区别的方程式

\bar{v} = \frac{Q}{\epsilon nA} = \frac{v}{\epsilon n}(2.83)

其中 ε 是一个由多孔介质特性的决定的经验常数。埃利斯(Ellis, 1968)等人在实验中得出的数据显示,当用相对均匀的砂样时,ε 值的范围是0.98—1.18。 对于非均匀的砂样和其它材料来说,目前为止 ε 值还没有得出。在地下水示踪剂和地下水污染的研究中,几乎通用的不须申明的假设是 ε = 1。对于颗粒介质来说,这可能引起细微差。在裂隙介质中,这个假设可能并不适用。

达西定律的上限和下限

如果我们仅从宏观规模上考虑通过达西连续体的单位面积流量,达西定律的适用性就有一定的局限。达西定律是一个线性定律。如果它是普遍适用的,那么对应于水力梯度 dh/dl 从0到∞,单位面积流量 v 与水力梯度的关系将是线性的 。对于通过颗粒物质的水流,在至少两种情况下,这种线性关系是有问题的;第一种是在很低的水力梯度下,通过低渗透率沉积物,第二种是高渗透率沉积物中大水流量的情况。换言之,应该有一个达西定律适用范围的上限和下限。因此,建议采用一个更具通用性的多孔介质水流规律的公式,即

v = -K\left(\frac{dh}{dl}\right)^m(2.84)

如果 m = 1, 就像它在所有普通情况中所表现的那样,其水流规律是线性的,并称之为达西定律。如果 m ≠ 1,其水流规律不是线性的,就不应称之为达西定律。

对于低渗透性的细粒物质来说,在已有的实验结论的基础上,人们认为这里可能存在一个临界水力梯度,即低于它就不产生水流。斯沃特曾德鲁伯(Swartzendruber, 1962),博尔特和格罗内威特(Bolt & Groenevelt, 1969)对这些实验结论进行了研究,并总结了被用来进一步解释这一现象的各种假设,只可惜在机理上还不能达到统一,并且实验证据上也还存在有一些疑问。无论如何,这一现象在实践上的重要性还较小。因为当水力梯度被认为可能是临界梯度的情况下,水流速度将非常小。

具有较大实践意义的是达西定律适用范围的上限。这个上限已经得到认可,并被接纳了多年(Rose, 1945; Hubbert, 1956)。对于高速水流,达西定律并不适用。托德 (Todd, 1959)和贝尔(Bear, 1972)两人的详细研究已经证实了这一结论。这个上限的确定,经常要借助于“雷诺数 Re ”来获得。它是一个无量纲的数,用以表示在水流动时惯性与粘滞力的比率。雷诺数广泛应用于流体力学中,以区分低速情况下的“层流”与高速情况下的“紊流”。通过多孔介质的水流的雷诺数由下式确定

R_e = \frac{\rho vd}{\mu}(2.85)

其中 ρμ 分别是流体的密度和粘度,v 是单位面积流量。d 是多孔介质的一个代表性长度量纲,根据不同情况它可代表:平均孔隙量纲,平均颗粒直径或者渗透率 k的平方根的一些函数。贝尔(Bear, 1972) 总结了试验结论,即:达西定律只有在以平均粒径为基础,且雷诺数不超出1—10这一范围时,才是有效的(126页)。对于雷诺数的这个范围来说,通过颗粒介质的所有水流都是层流。

超过达西定律这个上限的水流速度,一般出现在一些重要的岩石构造如溶洞石灰岩,白云岩和多孔洞的火山岩中。在非固结岩石和颗粒物质中,达西流几乎从不会超出该范围。裂隙岩 (我们将用这个名词去描述那些由于原始接触面、裂隙或开口而导致渗透率增大的岩石) 是需要单独注意的特殊情况。

裂隙岩中的水流

裂隙岩中的水流分析可使用本书已着重讨论过的连续体方法,或者单一裂隙中以水流水力学为基础的非连续体方法。在连续体方法中,与颗粒多孔介质类似,裂隙介质被连续体所替代。裂隙的水力学特征也被连续体空间的水力传导率、孔隙度和压缩系数所替代。但只有当裂隙岩体的密度足够大,使得它具有和颗粒多孔介质相似的水力传导率时这个方法才是有效的。虽然裂隙介质的代表性单元体积要比颗粒介质的大很多,但其概念的实质是相同的。如果裂隙间距在某个方向上是不同, 那么该介质呈现趋势性的非均质性。如果裂隙间距在某个方向上于其它方向都不同,则该介质将呈现为各向异性。斯诺(Snow, 1968, 1969)认为,任何裂隙水流问题,都可通过达西定律和一个各向异性的水力传导率张量所构成的标准多孔介质技术来解决。

如果裂隙的密度非常低,则需要进行单一裂隙中的水流分析。该方法已被用于岩土施工中,其中岩土力学分析指出,岩石的倾斜和开口会由于单个危险裂缝上的水流压力作用而遭到彻底破坏。这种分析方法建立于流体力学基本原理,纳维尔—斯托克斯方程之上,这里不再讨论。威特基(Wittke, 1973)撰写过一个介绍性的综述,可供参考。

即使我们应用连续体方法,仍有两个关于通过裂隙岩水流分析中的问题需加以论述。第一个问题是关于宽缝岩石裂隙中的非达西水流问题。夏普和梅尼(Sharp & Maini, 1972)提出了证实裂隙岩中非线性水流定律的实验数据。威特基(1973)建议在线性—层流范围(达西定律范围)、非线性层流范围和紊流范围内分别定义各自的流体定律。图2.28在单位面积流量与水力学梯度的曲线中囊括了这些概念。在宽的岩石缝隙中,单位面积流量和雷诺数相对较高,水力梯度通常小于1,且方程(2.84)中的指数 m 大于1。这些条件导致图 2.28中的曲线有一个向下的弯曲。

图 2.28 达西定律有效性的范围。

第二个问题是关于岩石中三维应力场和三维水流场的相互作用问题。关于两个场耦合的基本理论要求,已在第2.11节中作了简要讨论,并介绍了拜欧特(1941, 1955)关于多孔介质水流的经典研究工作。但是,对于裂隙介质来说,问题更复杂一些。因为裂隙岩的孔隙度很低,应力变化所导致的破碎裂隙的扩大和缩小,都会影响水力传导率 K 的值。因此,流体压力 p(x, y, z, t) 或水力水头 h (x, y, z, v) 与有效应力 σe(x, y, z, t) 之间的相互作用关系会变得更复杂,K必须用函数 K(σe) 来表示。这类系统的分析以及函数 K(σe) 本质的实验性研究,仍然是岩石力学和地下水水文学领域内需要进一步探索的内容。

许多将地下水理论应用到岩石力学中的学者提出了一个可以把裂隙的孔隙度 nf 和节理岩石的水力传导率 K,与节理的几何构造联系起来的公式。斯诺(1968)提到,对于垂直于岩石表面,单位距离内具有 N 个裂隙,裂隙的缝隙为 b 的平行裂隙组来说,nf = Nb,且

K = \left(\frac{\rho g}{\mu}\right) \left(\frac{Nb^3}{12}\right)(2.86)

k = \frac{Nb^3}{12}(2.87)

其中 k 是岩石的渗透率。Nb 的量纲分别是1/L和L,所以 k 的单位就变为 L2。方程(2.86)建立于一组平面裂隙的水动力学基础之上,并适用于达西定律有效的线性—层流范围内。它只能应用于具有足够大的可以认为是达西连续体的岩体中。用方程(2.87)计算的渗透率 k,可被看作是一个等同于多孔介质的渗透率;它的水力作用与裂隙岩相类似。

斯诺(1968)认为,该类裂隙岩体的立方体系统是一个孔隙度为 nf = 3Nb 的各向同性系统,其渗透率为该系统中任何一条裂隙的两倍,也就是说 k = Nb3/6。斯诺也提出了一个关于孔隙度与各向异性(裂隙间距或者孔径在各个方向都不同的三维节理几何构造)渗透率张量之间的预测性相关关系式。夏普和梅尼(Sharp & Maini, 1972)针对各向异性节理岩石的水力学性质进行了更深入的讨论。

2.13水动力弥散

在地下水流系统研究中,通过其迁移溶解物质(即熟知的溶质)的能力来评价水流势态问题已愈来愈普遍了。这些溶质可以是天然成分,人工示踪剂或者是污染物质。溶质借助地下水流的整体性运动得以迁移的过程,称为“平流”。由于平流的作用,(非反应性)溶质被水流以平均线性流速 \bar{v} 的所携带并发生迁移。但是溶质有离开平流水力学迁移路径的趋势。这个扩散现象被称为 “水动力弥散”,并导致了溶质的稀释。该现象的发生,是因为流体在平流过程中的机械性混合和溶质热动能所引起的分子扩散。扩散,是一种只在低速情况下才具有重要意义的弥散过程,这一点将在3.4节中加以叙述。本节的讨论着重于完全是被流体运动所控制的弥散。 也就是“机械弥散”或“水力弥散”。图2.29展示了均质颗粒介质中的弥散过程。

机械弥散是一种最容易被观察到的微观过程。在微观规模上,弥散的产生是源于三种机理(图2.30)。 第一种发生在单个孔隙通道内,因为分子从孔隙中不同位置流过时,由于孔隙表面的粗糙性导致流体速度发生变化而引起的。第二种过程是由水分子流经路径上的孔隙大小不同所引起的。由于单个孔隙通道内水体体积的表面面积以及通道粗糙度的不同,不同孔隙通道内的水流速度也不同。第三种弥散过程与孔隙通道的曲折性,分叉情况和重叠情况有关。溶质沿水流流动方向散开的过程被称为“纵向弥散”,沿垂直于流动方向散开的过程,称为“横向弥散”,纵向弥散一般比横向弥散强烈得多。

弥散是一种混合过程,从实质上说,它与地表水的紊流有相似的效应。对于多孔介质来说,平均线性流速与纵向弥散密切相关,纵向弥散是这样一种过程,一些水分子和溶质分子运动得比平均线性流速快,而另一些比平均线性流速慢。因此,溶质沿水流方向被散开,从而降低了浓度。

图 2.29 颗粒多孔介质中由于机械弥散作用产生的稀释过程。
图 2.30 微观过程中的弥散过程。

当在实验室内进行示踪试验时,唯一能够测量出来的弥散是那种在宏观上可以观察到的弥散。该过程假设这个宏观结果是由上述微观过程所导致的。一些研究者相信,宏观层面的非均质性可以在微观弥散的基础上产生额外的弥散效果。宏观弥散的机理尚未被充分认识,弥散过程将在第9章中继续探讨。

参考文献

BEAR, J. 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, pp. 15–24, 52–56, 85–90, 122–129, 136–148.

HUBBERT, M. K. 1940. The theory of groundwater motion. J. Geol., 48, pp. 785–822.

JACOB, C. E. 1940. On the flow of water in an elastic artesian aquifer. Trans. Amer. Geophys. Union, 2, pp. 574–586.

MAASLAND, M. 1957. Soil anisotropy and land drainage. Drainage of Agricultural Lands ed. J. N. Luthin. American Society of Agronomy, Madison, Wisc., pp. 216–246.

SKEMPTON, A. W. 1961. Effective stress in soils, concrete and rocks. Conference on Pore Pressures and Suction in Soils. Butterworth, London, pp. 4–16.

STALLMAN, R. W. 1964. Multiphase fluids in porous media-a review of theories pertinent to hydrologic studies. U.S. Geol. Surv. Prof Paper 411E.

VERRUIJT, A. 1969. Elastic storage of aquifers. Flow Through Porous Media, ed. R. J. M. De Wiest. Academic Press, New York, pp. 331–376.

习 题

  1. 下列野外记录是从并排安装在一个观测站的潜水位计中获得的(见下表)。 ABC 分别是潜水位计 abc 中测量的各点,试计算:

    潜水位计 a b c
    地面高程(平均海平面) 450 450 450
    潜水位计深度(m) 150 100 50
    潜水位埋深(m) 27 47 36
    1. BC 点的水力水头(米)
    2. BC 点的压力水头(米)
    3. BC 点的高程水头(米)
    4. B 点的流体压力(牛顿/平方米)
    5. AB, BC 间的水力梯度。你能否根据这些数据所指出的水流方向设想出一种水文地质情况?

  1. 绘制两个野外实际情况图,使每个图中的三个潜水位计并排装设,但是底部都在不同的深度,而水位却处在同一高程。
  1. 三个潜水位计彼此相距1000米,底部处在同一个水平含水层内,潜水位 A 在潜水位 B 的正南方,而潜水位 C 则在 AB 线的东面。ABC 的地面高程分别是95, 110和135米。A 中的水位埋深为5米,B 为30米,C 为35米,试确定通过三角形 ABC 的地下水流的方向,并计算其水力梯度。
  1. 对方程 \Phi = gz + p/\rho 进行量纲分析,以说明流体势能Φ是一个能量项。并用SI制和FPS制两种单位表示。
  2. 有三个地层,每个厚25米,互相叠置。如果在这个结构中设置一个不变流速的垂向水流场,使其顶部 h = 120米,底部 h = 100米,试计算其内部两个边界处的h值(设顶部地层的水力传导率为0.0001米/秒,中部地层为0.0005米/秒,底部为0.0010米/秒。
  1. 设某地质结构的渗透率为0.1达西(根据某石油公司为石油流动所确定的数值),那么对于水流来说,该结构的水力传导率是多少?请用米/秒和加仑/天/平方英尺两种单位来回答。同时推测这个结构像是何种岩石?
    1. 有四个水平、均质、各向同性的地质结构互相叠置,每层5米。如果其水力传导率分别是:10—4,10—6,10—4和10—6米/秒,试计算其等效的均质却各向异性结构的水力传导率的水平与垂向分量。
    2. 进行重复计算,设上述水力传导率为10—4,10—8,10—4和10—8米/秒;10—4,10—10 ,10—4和10—10米/秒。把这三组计算成果填入一张表格,使层状非均质性的数量级与其相对应的等效各向异性相结合。

    1. 通过孔隙度和空隙比的体积定义推导方程(2.40)中所给定的关系式。
    2. 当在同一个土质样品上测量时,孔隙度是否有可能比空隙比大?

  1. 设某土质水份测量点处地表面高程为300厘米,土质为砂,同时其未饱和的孔隙度可用图2.13中的干化曲线来代表,试绘出在下列条件下,在200厘米深处,对应于深度的含水量、压力水头和水力水头在数量上比较精确的一组垂向剖面图(如图2.12中所示):
    1. 整个剖面图的含水量都是20%。
    2. 整个剖面图的压力水头是—50厘米。
      整个剖面图(静态的)的水力水头是150厘米。

    对于情况(a)和(b),计算其水力梯度和剖面的水流速。对于情况(c), 确定到水位面的深度。

  1. 已知一个区域等势面的坡度为 7m/km,计算导水系数T = 0.002 m2/s 的承压含水层中的地下水自然流速。
  1. 使用量纲分析证明储水系数 S = \rho g b (\alpha + n\beta) 是无量纲量。
    1. 在一个水平含水层上方有厚度为50英尺的饱和粘土。粘土的比重(或单位干重)为120 磅/英尺3。水的比重为62.4 磅/英尺3。计算含水层上方的总应力。
    2. 如果含水层的压力水头为100英尺,计算含水层的有效应力。
    3. 如果从含水层中抽水,一些采样点的水力水头下降了10英尺,求压力水头的下降值,流体压力,有效应力和总应力。
    4. 如果厚度为25英尺的含水层的压缩性为10—6英尺2/磅,在经历(c) 的水头下降的过程中,含水层会有多大程度的压缩?
    5. 如果含水层的孔隙度和水力传导率为0.30和10加仑/天/英尺2

  1. 回顾以下经典定义衍生出的问题:等势面,渗透率,