等效渗透系数计算方法

A. 模型的特色及改进

一、混合井孔的模拟

为了简化成井工艺，增加出水量和降低成本，本区地下水开采多采用混合井（图5-2）。

混合井在地下水流模型中的刻画，目前大都采用由美国地调局1988年推出，并在国际上广泛应用的三维有限差分地下水流模型MODFLOW软件处理。

图5-2 混合抽水井

MODFLOW建议：多层井的流量必须以某种方式人为地分配给每一单层，……把井流量按每层的导水系数大小分配，即Q_i/Q_w=T_i/∑T（1988，1996，2000年）。其中Q_i和Q_w分别为第i层流量和总（井口）流量，T_i和∑T分别为第i层导水系数和总导水系数。

为了便于讨论，又不失其一般，我们以贯穿两个含水层的混合井孔（图5-2）为例进行讨论。如此，上式可表述为Q₁/Q₂=T₁/T₂。

MODFLOW对混合井的这种处理方法没有给出理论上任何的分析、说明，缺乏理论依据，应用中也与实际不符（陈崇希等，1998；Chen Chong-xi和Jiao J J，1999）。这是因为：

1）含水层导水系数对井孔流量的影响，不会如此简单。例如，混合井附近岩性（渗透系数）发生变化，甚至混合井打在岩性透镜体上（图5-3），怎样影响各层流量的分配?含水层厚度发生变化又如何改变流量的分配?MODFLOW无法解答这些经常遇到的实际问题。

图5-3 混合井管贯穿岩性透镜体

2）含水层的参数影响混合井流量的分配，导水系数只是其中一个因素，含水层的弹性给水度（储水系数）也应该起作用。

3）井管的流量分配不仅与含水层的水文地质参数分布有关，还与外边界条件、井径（有效井径）、水泵吸水管的位置及其他抽水井的干扰等因素有关，MODFLOW仅依导水系数对混合井的流量进行预分配，未能考虑其他因素的作用。

4）含水层的导水系数一般不随时间变化，即导水系数比T₁/T₂不会变化，但流量比Q₁/Q₂却不是一个常量。

5）按MODFLOW的方法，一个混合抽水井的流量Q₁/Q₂始终是个常量。然而在模拟过程中，特别是预测时，周围随时可加入或关闭混合的或非混合的抽水井，在这种井群干扰下，原混合井的流量比还会保持常量吗?显然是不可能的。

6）从另一角度分析，混合观测孔是混合抽水井的特殊情况（Q_w=0），对于两层混合的观测孔，其孔中水位（混合水位）必介于两含水层水位之间，即混合观测孔对于其中一含水层（例如1含水层）起抽水作用（Q₁＞0），对于另一含水层（2含水层）起注水作用（Q₂＜0）。如此，Q₁/Q₂＜0。而两含水层的导水系数的比值肯定是正值，即T₁/T₂＞0。如此，两个比值怎能相等?有人认为：MODFLOW并没有说上式可用于混合观测孔。对此，很容易论证：当混合抽水井的抽水量足够小，足以保持混合井中的水位介于上下两含水层的水头之间，上述论证一样成立。

7）按各层导水系数T的比例来预先人为划分各层的流量。这是缺乏理论依据的，因为这种做法要求各分层的有效井径相等和井壁处的水力坡度上下处处相等，这两个条件不可能人为控制，预先也不得而知。

显然，用导水系数的比值预先给定各分层的流量是不妥的。实质上，这种方法不是模拟，而是“处理”，一种未考虑机理的“处理”。

“防止地下水模拟失真，提高仿真性”是水文地质模拟工作者的核心任务，而地下水流系统中普遍存在的混合井，又是当前国内外模拟失真的主要问题之一，应当引起我们的足够重视。

研究区的混合井，既有抽水井，也有观测孔。三维流场中，即使属于均质含水层（无需多层含水系统），常规（理论上非点状滤管）抽水井和观测孔都属于混合井孔，因为滤管中不同深度处的水头是不相等的，因此滤管中的水要发生垂直流动，即使在不抽水的观测孔中也一样。这就是地下水流对混合井孔响应的本质所在。

本研究采用“渗流-管流耦合模型”（陈崇希等，1992，1996）来刻画混合井孔，以“渗流”刻画地下水的运动，“管流”刻画井孔中的水流，以解决混合井孔的模拟问题，大大地提高了模型的仿真性。“渗流-管流耦合模型”的基本思路如下（图5-4）：

图5-4 计算框图

1）将产生管流的混合井视为渗透系数很大的“圆柱形透镜体”，含水层通过这个具有很大导水性的圆柱体（井筒）强烈地交换水量，那么混合抽水问题就可以视为一个特殊的“越流系统”。

2）求越流系统“圆柱形透镜体”的“渗透系数”。需要强调的是，“圆柱形透镜体”的“渗透系数”不是常量，而是随流态的变化而变化。

具体的计算方法如下：

由流体力学知识可知，当管流呈层流状态时，其水头损失可依Darcy-Weisbach方程计算：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ΔH为水头损失；f为摩擦系数；l为管长，m；d为管内直径，m；u为管内平均流速，m/d；g为重力加速度，m/d²。

当管流为层流时

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：Re为雷诺数。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ν为流体的运动黏度；μ为流体的动力黏度；ρ为流体密度。

将式（5-2）至（5-4）代入式（5-1），得

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：J为水力坡度；γ为流体重度。

将上式写成渗透流速的形式，对于管流，其空隙率n=1，则渗透流速v=nu=u，故有

河西走廊疏勒河流域地下水资源合理开发利用调查评价

将此式与达西定律v=KJ对比，可得出层流状态下管流的等效渗透系数K_L的表达式

河西走廊疏勒河流域地下水资源合理开发利用调查评价

此式由陈崇希（1966）和J.Bear（1972）分别获得。

根据流体力学知识可知，当雷诺数大于100000时，摩擦系数f与雷诺数Re无关，而取决于管内壁的相对粗糙度（），这时的水头损失ΔH与流速的平方u²成正比。在3000＜Re＜100000的范围内，其中存在f=，即ΔH∝u^1.7⁵的区段。在紊流条件下，除了上述两个区段外，还存在两个过渡区，共分为4个区段。可见，紊流是个比较复杂的问题。为了解决此问题，这里提出等效渗透系数K_N的概念。

当管流呈紊流状态时，式（5-1）可改写为

河西走廊疏勒河流域地下水资源合理开发利用调查评价

如果定义紊流状态下管流的等效渗透系数

河西走廊疏勒河流域地下水资源合理开发利用调查评价

则

河西走廊疏勒河流域地下水资源合理开发利用调查评价

该式也具有达西定律的形式。

本来，紊流的运动规律有别于达西定律，而且不同流态区具有不同的形式。引入等效渗透系数后，将5个流态分区（1个层流区和4个紊流区）的运动规律统一为达西定律形式，且与地下水渗流定律一致。即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

其中K_e为渗流-管流耦合模型的等效渗透系数。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

需注意的是，紊流条件下的等效渗透系数K_N与流速v及摩擦系数f有关，它是个随雷诺数而变化的量。另外，雷诺数Re的确定与流速v有关，而流速v又依赖于摩擦系数f，摩擦系数f的确定反过来又取决于雷诺数Re，因此必须采用迭代法来确定三者，进而确定等效渗透系数K_e。具体计算流程如图5-4。

二、混合观测孔水位的模拟

由于种种原因，几乎没有无混合观测孔的地区。什么是混合观测孔，混合观测孔中的水位（混合水位）如何形成，如何根据混合水位来确定初始水头的分布，如何利用混合水位求取各分层水文地质参数，等等，是“防止模拟失真，提高仿真性”（陈崇希，2003）的重要问题之一。然而MODFLOW（McDonald等，1988）对混合观测孔没有做任何分析、讨论，也没有混合观测孔的模块。简单地放弃混合观测孔的信息是不对的，因为“混合观测孔”只是孔口的流量为零，混合观测孔内却存在“抽水”与“注水”。

D.Sokol（1963）采用Thiem“影响半径”稳定井流模型证明混合观测孔的水位等于以导水系数为权重的各层水头的代数平均值。该论文存在如下几个主要问题：①Thiem“影响半径”模型是不能形成稳定流的（陈崇希，1966，1975，1983）；②Sokol隐含着假定，尽管两含水层的导水系数不等，流量也不同，但取其“影响半径”相等；③Sokol还隐含着假定，井管的阻力可以忽略不计，即井管处处水头相等。在这些假定条件下得出的结论是不可信的。

Hantush（1961）和Бочевер等（前苏联）（1961）在对承压含水层非完整井流的研究中，都认为观测孔中的水头降深s（r，l′，d′，t）反映该孔滤水管中各点降深s_p（r，z，t）的平均降深，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：l′，d′是观测孔滤管顶点、底点的标高。本书称之为Hantush-Бочевер方程。

对于潜水三维流的研究，Neuman（1972）也认为，观测孔中的水头降深可视为滤管内各点降深的平均值，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：r是观测孔至抽水井的距离；z₁、z₂是观测孔滤管顶点、底点的标高。

上述诸学者的见解与确定观测孔中水位的方法在水文地质界一直沿用至今，40多年来未见到异议。然而，上述观测孔中水位的确定方法，缺乏对形成机理的基本分析。简单分析之，三维流中的垂直观测孔，由于观测孔滤管中的水头不等，会发生垂向流动，而水的流动则会导致水头的变化，并引起井孔周围地下水的运动与水头再分布。因此，与测压计式观测孔不同，一般观测孔实质上是一部分滤管进水（抽水）而另一部分滤管出水（注水），即观测孔并不简单地反映含水层地下水的水头，而是它兼有“抽水”与“注水”作用的井孔，只是孔口的流量为零，井管内的“抽水”量与“注水”量绝对值相等而已（如果忽略微小的井筒储存量的变化）。观测孔非“观测”孔也！

显然，上述Hantush、Бочевер等以及Neuman提出的“积分平均水位”都未涉及水流的机理，是缺乏物理基础的纯数学方法。作者曾用“渗流-管流耦合模型”来模拟混合观测孔中的水位（混合水位）。算例的条件是：

一个均质、等厚、水平无限延伸的承压含水层，于顶面处有一微小的半球状井定流量抽水，抽水延续时间t=10.02d，在不同径距r处有不同孔径的完整观测孔（此即混合观测孔），模拟这些观测孔中的水头分布和流量变化。含水层的厚度M=100m，渗透系数K=100m/d，单位弹性给水度μ_s=0.00001m^-1，抽水流量Q=36000m³/d。模拟结果示于图5-5和图5-6。由图5-5可见，孔径≤0.1m的观测孔中，不同深度处的水头值明显存在差异，与Hantush-Бочевер的积分平均相差也很大。

图5-5 不同径距r观测孔中水头降深分布及Hantush-Бочевер公式计算值s_ga

图5-6 不同径距观测孔流量分布

由于抽水井（点汇）位于承压含水层顶面处，因此所有径距r处的顶部处的水头降深s大于底部，即底部的水头H大于顶部的水头，于是观测孔中的水从下而上流动。基于水流连续性原理，底部观测孔要吸水（抽水）而顶部观测孔要排水（注水），图5-6表示观测孔中的流量分布。图5-6a表明：径距21.19m，直径0.20m的观测孔中，大约在z=35m以下为“抽水”，35m以上为“注水”；观测孔中垂向最大流量超过500m³/d。这就是我们所说的，观测孔非“观测”孔也！

基于上述分析，本项研究采用“渗流-管流耦合模型”来模拟混合观测孔中的水位，使混合水位成为有用的信息，即用于初始水头分布和拟合求参。

三、泉流量动态的模拟

泉是地下水转化为地表水的主要形式之一，不同类型泉的出现及其流量动态，是水文地质工作者要重点分析的问题，因为这些数据为认识水文地质条件提供了极为重要的信息。泉是数值模型中主要的模拟要素；泉流量的预测是地下水资源评价、管理不可缺省的内容。

美国地质调查局（2001）建议采用MODFLOW软件的Drain模块迭代计算泉流量，实质上将泉作为第二类边界来刻画，即先不考虑泉的存在求解水头分布，再用泉所在格点的地下水水头与泉口标高之差乘以某比例系数（该系数缺乏物理意义）计算出泉流量，再以此流量置于该格点（第二类边界）重新求解水头分布，再次重新计算出新的泉流量，如此迭代直至收敛为止。当地下水水位低于泉口标高时，流量为零；当地下水水位高于泉口标高时，泉流量与水头差成正比，其比例系数必须通过流量的拟合来反算。

我国以往数值模型中关于泉流的刻画，采用的像抽水井一样输入（给定）流量。如此处理，预测怎么办?泉流量是未知的。近年来，我国也有类似上述MODFLOW建议的迭代算法。例如有的研究报告将鄂尔多斯大向斜的东侧山西省柳林泉（实为单斜自流斜地上升泉）刻画为平面二维流模型，进而用承压完整井的公式计算泉的流量，即Q=qs。这里存在几个问题：①柳林泉是上升泉，泉附近的钻孔愈深其地下水头也愈高，这是单斜自流斜地上升泉的特点，这种条件具有明显的三维流特征；②用承压完整井的公式计算泉的流量，特别是柳林泉的流量，不妥；③比例系数q在不同时间（不同条件）并不是常量。

“泉，必需将其放在水文地质体中去研究其成因类型，才能正确模拟仿真”。本模型提出的方法不同于MODFLOW等，是将泉口标高作为第一类边界条件（当泉流存在时），并将泉口格点与同层周边格点及下格点分别依达西定律（如果地下水流属于线性流，例如陈崇希等（1995），Chen Chongxi等（1996），陈崇希等（2003）；或采用非达西定律，例如陈崇希（1995），成建梅等（1998），ChengJianmei等（1998））和水均衡原理建立关系，整个计算方法建立在流动机理之上，因此由模型运转直接产生泉流量，无需通过迭代求解（且有的情况下迭代并不收敛）。正因为此，在模型识别阶段才有可能将泉流量作为拟合对象。这一点十分重要。

具体地说（图5-7），根据达西定律，泉流量Q为水平向的流量与垂向上的流量之和，如果以水平向的流量为主，则为下降泉，否则为上升泉。具体可通过下式计算：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：Q为泉流量；K_zsp为泉流区域的垂直渗透系数；H_b为泉口下层结点的水头值；Z_sp为泉出露标高；D_zsp为泉流过程中的渗流长度（含水层厚度）；S为泉口结点控制面积；C_A为泉颈修正系数；C_s为泉径修正系数；T_ij，T_ik为流段ij和ik的平均导水系数；，，，为线段的长度；h_i，h_j，h_k为结点i，j，k的水头值；e为结点i周围的均衡单元，如ipbq四边形。

图5-7 泉的流量计算示意图

在以往的地下水资源评价数值模拟中，没有将泉的流量作为拟合要素。实际上，象MODFLOW等的做法也不可能将泉的流量作为拟合要素，因为他们在模型识别其间要求解一个缺乏物理意义的系数。

本模型选择了3个典型点泉和22条线泉（泉沟），模拟其流量动态，并将其实测流量作为一个重要的拟合对象用以求取含水系统的水文地质参数和预测泉流量的动态。如此将会提高模型的仿真性。

四、大气降雨、地表水等入渗补给潜水滞后性的刻画

降水、地表水、渠系等入渗补给潜水是地下水的重要补给来源，若处理不当，则会导致模型失真。因此，正确的刻画入渗补给的滞后性对于模型仿真性的提高具有重要的意义。

目前地下水资源评价中降水、地表水、渠系等入渗补给量的计算主要采用下列两种方法：

1）未考虑滞后的入渗系数法。该法对于小埋深区域大体上可以用来刻画入渗对潜水的补给，但有些地区潜水埋深较大，若不考虑降雨、地表水等入渗补给的滞后特征，则会导致模型失真。

2）平移滞后入渗系数法。有的研究报告把潜水位出现高峰时间与降雨高峰出现时间的时间差，作为统一的滞后时间。例如若该时间差为3个月，则1月份的降雨量都在4月份入渗补给，2月份的降雨量则在5月份入渗补给；若3月份无降雨，则6月份无降雨入渗补给。如此处理滞后补给也存在明显的不合理性。因为潜水位高峰出现的时间，是所有影响地下水位动态诸因素作用的总和，除降雨补给外，还有地下水蒸发，灌溉入渗，越流和地下水开采等因素作用下产生地下水不平衡流动的结果，特别是地下水开采这个人为因素在本区已经成为地下水动态的主要控制因素，它已经完全掩盖了地下水的天然动态，因此用该法研究降雨入渗补给的滞后性是不恰当的。即使这一地区完全处于未开采状态，甚至忽略除降雨以外的其他因素对地下水动态的影响，仅就降雨入渗补给单一因素对潜水位动态的影响而言，该法也存在明显的问题。降雨过程是间断的，不连续的，但它对潜水的补给则是连续的，这是因为地下水流动，特别是非饱和流动，其介质具明显的储存性而有调节功能，从而使补给具有滞后性及连续性。一年内虽然降雨仅有几十天，其余时间均无降雨，但在潜水位深埋区，补给却天天作用着（浅埋区由于存在强蒸发作用而使情况复杂化了），这种连续补给作用是不可能用一个时间差的平移来刻画其滞后性的。

上述“入渗系数法”在潜水小埋深区有一定的适用性。但由于本区部分地段潜水埋深较大，如昌马洪积扇顶部潜水水位最大埋深达290m上下。一个月的降雨量乘以入渗系数的水量不可能在当月内全部补给其下的潜水含水层。因此，不考虑降雨补给的滞后特征，则会导致模型的失真。

本研究采用“降雨补给滞后权系数法”（陈崇希等，1991，1998）来刻画降雨入渗补给。该方法既能基本反映客观实际情况，又很实用。

根据地下水非饱和流动理论，一次降雨过程对潜水补给量的分布曲线如图5-8所示。它是一单峰曲线。曲线的形态取决于潜水位的埋深和包气带的岩性，当然还与包气带的初始含水率的分布有关。我们将后面的因素固定起来，专门讨论前面因素的作用。

图5-8 一次降雨对补给量的分布曲线图

我们以包气带地下水非饱和流动理论为基础，考察不同水位埋深和岩性条件下，某时段（例如月、旬等）降雨对该时段及其后各时段补给强度的分布特征。为了使问题的讨论能与地下水资源评价数值法相匹配，将连续的特征曲线改为离散点的形式（因为数值法在时间上是离散的）。不失一般性，我们将时段取为1个月。

对于潜水位埋深较大和（或）包气带渗透系数较小的条件，入渗强度分布是一单峰函数，当潜水位埋深很小时，可能是单调减函数（图5-9）。若在0月有某降雨量RA（mm）作用，则该月及后各月降雨入渗补给量之和∑RE_k（mm）就是该时段降雨量形成的全部补给量。该值∑RE_k与0月的降雨量RA之比就是（总）入渗补给系数α。各月份的入渗量与0月降雨量RA之比，则是0月降雨在各月的入渗补给系数α_k。即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

图5-9 降雨对不同埋深（岩性）条件的入渗补给历时分布曲线图

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：RA为0月的降雨量，mm；RE_k为0月的降雨在k月的入渗补给量，mm；RE为（总）补给量，mm；α_k为0月降雨在k月的入渗补给系数；α为（总）入渗补给系数。

如果令图5-9上0月的降雨量为100mm/月，则该图纵坐标转变为以百分数表示的月入渗系数α_k（k=0，1，2，3，…）。

月入渗补给系数组已能很好地反映降雨入渗补给滞后性的特征，为了使刻画滞后性的特征数组更通用，我们定义：

河西走廊疏勒河流域地下水资源合理开发利用调查评价

因此

河西走廊疏勒河流域地下水资源合理开发利用调查评价

该式表明ω_k满足权的定义，所以我们称ω_k（k=0，1，2，…）为滞后补给权系数（组）。

引入滞后补给权系数（组）的优点是，不同入渗系数的诸分区（不同区，同月的入渗补给量不等），有可能取用同一滞后补给权系数来计算各月的补给量，也就是说，滞后补给权系数更具通用性。于是我们的任务转到寻找一个能满足上述条件的特征函数，即：

1）它们是单峰函数或单调减函数，并随峰值（或最大值）出现时间的后延，其峰值（或最大值）减小；

2）各权系数之和为1；

3）适应性较强，能较好地与不同潜水位埋深，不同岩性条件下滞后入渗补给量相拟合。

经反复分析研究，我们确定采用下面的离散函数来表征滞后入渗补给权系数的分布。

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：ω（j-i）为i月降雨，j月补给的权系数；h为潜水位埋深，m；B为与包气带岩性等因素有关的经验系数，m；可称为埋深区间系数。该系数是包气带渗透系数的增函数。

该式中h是已知的，仅有一个与岩性等因素有关的经验系数B是待求的。这个系数可以通过模型识别来确定，这样就可以利用上式计算出滞后入渗补给权系数。各区（总）入渗系数也是通过模型识别来确定的，由此可计算出在某月i降雨量RA（i）作用下，该月及其后各月j的补给量RE（j-i），即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

至此，我们分析了某月降雨量对该月及其后各月入渗补给的贡献。对于实际问题——多月降雨对某月j的补给量问题，为了简化计算，我们近似用该月及其以前各月i降雨量对潜水补给贡献的代数和来计算该月j获得的总入渗补给量，即

河西走廊疏勒河流域地下水资源合理开发利用调查评价

式中：RES（j）为j月及其以前各月的降雨量在j月的总入渗补给量，mm；RE（j-i）为j月降雨对j（j＞i）月的入渗补给量，mm；α为入渗补给系数；ω（j＞i）为i月降雨对j（j＞i）月的滞后补给权系数；RA（i）为i月的降雨量，mm。

依式（5-17），获得的降雨入渗补给滞后权系数曲线共75条，它对应于75个潜水位埋深区间，每个区间埋深差Bm（即埋深区间系数）即每增加埋深Bm，换一条权系数曲线的序号）。图5-10中仅绘出12条典型的权系数曲线，而各两条曲线之间，尚有5～7条曲线图中并未绘出。

降雨补给滞后权系数法的提出和使用，不仅大大提高了模拟的仿真性，而且使用十分方便。地表水入渗补给滞后性的机理与降雨入渗的类似，也采用上述方法刻画。

图5-10 降雨入渗补给滞后权系数曲线

五、参数的分区

本模型含水介质厚度大，垂向上共分6个模拟层，即便是同一岩性（如亚粘土）而处于不同埋藏深度，由于其固结程度不同而导致其水文地质参数不相等。这种情况下如何给定初始参数估计值的空间分布又是一个难题。若按常规惯用的方法，则同一岩性根据不同层位（不同埋藏深度）要给出不同参数分区，如此的参数分区数过多，给调参带来很大的工作量。

对此问题，本模型采用如下方法，处于不同埋藏深度同一岩性的K_h值作为同一分区，但其参数值随埋深的负指数衰减。这样，每个参数分区只要1个零埋深的水平渗透系数K_h0和1个衰减系数γ刻画，即K_h=K_h0×e^-γh（h为埋深值）。其他参数，如K_z值、μ_d值和μ_s值的空间分布也按同样方法给定。如此处理，就大大地减少了模型参数的数量，相应地也减少了模型调参的工作量。

六、自流井的模拟

自流井普遍地存在于我国西北内陆盆地，冲洪积扇的前缘洼地，大型岩溶上升泉的周围以及平原区的深部承压含水层等区段。我国至今还存在大面积的地下水自流区（尽管由于大量开采地下水而缩小了其自流面积）。因此模型中如何刻画自流井是一个重要的水文地质问题。然而，美国地调局推出的MODFLOW软件尚无模拟自流井流量的模块，这似乎不是疏忽，而是存在一定的难度。常规的地下水流模型已经解决不了这类问题。至今，对于数值模型中的自流井，仍是无可奈何地采用人为地给定其自流量。然而，自流井与抽水井不同，后者人们一般可以人为给定抽水流量，而井中水位是作为其响应；当井中水位降至井底时，人们就不能自由地抽取了（尽管如此，大多含有抽水井的数值模型，未能客观地将“井中水位作为抽水量的约束”来认识）。然而，自流井就不能人为地给定自流量了，因为自流量是自流井对周围环境（气象、水文、地下水的开采等因素）的响应，自流井的流量是属于预测范畴，是输出信息，而不是输入信息。

本项目采用“渗流-管流耦合模型”模拟自流井，只要在自流井的井口处设置水头已知边界，即水头等于井口标高，即可模拟出自流井的流量。该法已在“陕西渭北东部岩溶水资源数值模拟研究”中获得成功的应用（陈崇希等，2003）。

七、初始水头的形成

初始水头的分布是地下水不稳定流数值模拟不可缺少的条件。通常利用一定数量观测孔的水位通过插值获得各结点（格点）的初始水头值。然而，当研究区范围较大，观测孔较少，观测孔所在的层位不一，而且含有混合孔时，难以用插值的办法获得每一模拟层的初始水头分布。

本项目采用“参数-初始水头迭代法”（陈崇希、林敏等，1991；陈崇希、裴顺平，2001），来确定各层的初始水头分布。具体处理办法是：取模拟的初始时间为2003年7月，用水文地质参数初始估计值从2001年7月开始模拟，至2003年7月将模拟所得的该月水头值与实测值拟合，粗略确定模拟时段的初始水头（2003年7月）；再以此初始水头按常规方法求参，如此迭代求解。当拟合达到一定程度之后，必须全面地检验模型是否符合基本规律，否则要再次重新调整，直至满足要求。上述方法虽然工作量很大，但效果良好。

B. 土石混合体渗透性能的正交试验研究

周中¹ 傅鹤林¹ 刘宝琛¹ 谭捍华² 龙万学² 罗强²

（1.中南大学土木建筑学院 湖南 长沙 410075

2.贵州省交通规划勘察设计研究院 贵州 贵阳 550001）

摘要 土石混合体作为土和石块的介质耦合体，具有非均质性、非连续性及试样的难以采集性等独特的性质，从而给研究带来极大的困难。土石混合体属于典型的多孔介质，其渗透特性与颗粒的大小、孔隙比及颗粒形状关系密切。本文采用室内正交实验，利用自制的常水头渗透仪，研究了砾石含量、孔隙比和颗粒形状三个因素在不同水平下对土石混合体渗透系数的影响。通过正交试验确定了三种因素对土石混合体渗透系数的影响顺序及各因素的显着性水平。提出了土石混合体渗透系数计算公式，并通过试验结果验证了计算公式的正确性，为土石混合体渗透系数的理论计算提供了一个简明有用的计算工具。

关键词 土石混合体 多孔介质 渗透性能 计算公式 正交试验

土石混合体一般由作为骨料的砾石或块石与作为充填料的粘土或砂组成，是介于土体与岩体之间的一种特殊的地质体，是土和石块的介质耦合体^［1］。因为土石混合体具有物质组成的复杂性、结构分布的不规则性以及试样的难以采集性等独特的性质，从而给研究带来极大的困难，目前人们对于它的研究仍处于探索之中^［2］。渗透与强度和变形特性，都是土力学中所要研究的主要力学性质，其在土木工程的各个领域中都有重要的作用^［3］。土石混合体属于典型的非均质多孔介质^［4］，其渗透特性与颗粒的大小、颗粒组成、孔隙比及颗粒形状关系密切。土的渗透系数可以通过室内试验由达西定理计算得出，然而土石混合体的渗透系数却难以确定，主要原因是：取样困难；难以进行常规的渗透试验；大尺度的渗透试验不仅造价高准确性差，而且试验结果离散度大，难以掌握其规律性。迄今为止，国内还没有对土石混合体渗透性能进行研究的资料，现有研究成果局限于利用物理和数值模拟试验对其变形和力学性质进行研究，而对渗透性还未涉及。因此，能够求出土石混合体渗透系数的计算公式具有重要的理论意义和工程应用价值。

本文研究土石混合体中砾石含量、孔隙比（压实度）和颗粒形状三个因素在不同水平下对土石混合体的渗透系数的影响，找出三因素与土石混合体渗透系数之间的关系，并提出土石混合体渗透系数计算公式。

1 土石混合体渗透性能的正交试验

1.1 正交试验方案设计

在室内试验中考虑砾石含量、孔隙比（压实度）和颗粒形状三个因素对土石混合体渗透系数的影响，就每种因素拟考虑3个水平。对于这种3因素3水平的试验，如果考虑每一个因素的不同水平对基材的影响，则根据组合可得有3³组试验，这对人力、物力与时间来说都是一种浪费，因此采用正交试验设计来研究这一问题更为合理。本试验所选取的正交表为L₉（3⁴），考虑试验误差的影响，但不考虑各因素间的交互作用（即假定他们之间相互没有影响）。共需9组试验，每组作平行试验3次，共27次渗透试验。本试验中采用的因素与对应的水平数如表1所示，其中粗粒形状分为球形体、六面体和三棱锥3个水平，分别由卵石、强风化石块和新打碎的碎石来近似替代。

表1 正交试验的因素水平

1.2 试样的基本物理力学性质

试验所取土样为正在修建的上瑞高速公路贵州段晴隆隧道出口处典型性土石混合体，其天然状态土的物理指标及颗粒级配曲线见表2和图1。由图1可知现场取回土样的不均匀系数C_u为12.31，说明土样中包含的粒径级数较多，粗细粒径之间差别较大，颗粒级配曲线的曲率系数C_c为1.59，级配优良。

表2 天然状态土的基本物理指标

图1 天然状态土的颗粒级配曲线

1.3 大型渗透仪的研制

《土工试验规程》（SL237—1999）规定粗粒土的室内渗透系数需由常水头渗透仪测试，国内常用的常水头渗透仪是70型渗透仪。70型渗透仪的筒身内径为9.44cm，试验材料的最大粒径为2cm，规范^［5］要求筒身内径应为最大粒径的8～10倍，因此70型渗透仪的筒身内径过小，有必要研制大尺寸的渗透仪。自制渗透仪的内径和试样高度至少应为最大颗粒粒径的8倍，即至少应为16cm，另外，考虑到边界效应，试样的上下两头分别增加2cm，因此，自制渗透仪的内径和试样高分别取为16cm和20cm。考虑到土石混合体的渗透性较强，选取进排水管的口径为2cm。自制的大型常水头渗透仪的如图2 和图3所示。

图2 常水头渗透仪示意图

数据单位为cm

图3 自制渗透仪

2 试验结果分析

2.1 试验结果

按正交试验表L₉（3⁴）的安排，共需作9组试验，每组试验作平行试验3次，取3次测量的平均值，并乘以温度校正系数

，即可求出每组试验20℃时的渗透系数，渗透系数的测量结果见表3。

表3 渗透试验测定结果

续表

2.2 试验分析

运用正交试验的直观分析法和方差分析法，分析各因素对土石混合体渗透系数影响的主次顺序，绘出因素水平影响趋势图，求出各因素的显着性水平。

2.2.1 直观分析

对试验所得的土石混合体的渗透系数进行正交试验的极差分析，并画出各因素的水平影响趋势图。正交试验的极差分析表见表4，3个因素与渗透系数的关系见图4。

表4 极差分析表

图4 各因素与渗透系数的关系

A—砾石含量；B—孔隙比；C—粗粒形状

由正交试验的极差分析表可以看出，对土石混合体渗透系数影响的主次顺序为A→B→C，即砾石含量→孔隙比→颗粒形状。由各因素与渗透系数的关系图可以看出砾石含量越多渗透系数越大，孔隙比越大渗透系数越大，颗粒磨圆度越大渗透系数越小。在路基工程及大坝工程中，可以通过调节粗颗粒的含量、压实度及颗粒形状以获得工程所需的渗透系数。

2.2.2 方差分析

为了确定因素各水平对应的试验结果的差异是由因素水平不同引起的，还是由试验误差引起的，并对影响土石混合体渗透系数的各因素的显着性水平给予精确的数量评估，需采用正交试验的方差分析法对试验数据进行分析，分析结果如表5所示。

表5 方差分析结果

方差分析结果表明：

（1）因素各水平对应的试验结果的差异是由因素水平不同引起的，而不是由试验误差引起的；

（2）砾石含量对土石混合体渗透系数的影响高度显着，孔隙比对土石混合体渗透系数的影响显着，颗粒形状土石混合体渗透系数的影响不显着。

3 土石混合体渗透系数

3.1 渗透系数与砾石含量之间的关系

众所周知，土石混合体的渗透系数与颗粒的大小及级配有关，本文选择等效粒径d₂₀和曲率系数C_c来表示土的颗粒大小和颗粒级配，原因是文献［3］认为等效粒径d₂₀比其他粒径特征系数更能准确地表示颗粒的大小，而与颗粒级配有关的系数是不均匀系数C_u和曲率系数C_c，不均匀系数C_u只反映土粒组成的离散程度，曲率系数C_c能在一定程度上反映颗粒组成曲线的特性，因而曲率系数C_c更适合于评价土的颗粒级配。不同砾石含量的颗粒级配曲线如图5所示。由图5可以求出各曲线的粒径特征系数，见表6。

图5 试样的颗粒级配曲线

表6 不同粗粒含量时的粒径特征

由图6可知，其他条件相同时，土石混合体的渗透系数k与函数f（d₂₀，C_c）呈线性关系，其中

。

图6 k₂₀-f（d₂₀，C_c）关系曲线

3.2 渗透系数与密实度之间的关系

由正交试验的方差分析可知，孔隙率e对渗透系数的影响虽不如粗粒含量大，但也是很显着的。在其他条件相同时，k与

呈线性关系，如图7所示。

土石混合体

3.3 渗透系数与颗粒形状之间的关系

狄凯尔与海阿特（Tikell and Hiatt）于1938年探讨了颗粒的“棱角性”与“圆度”对渗透系数的影响，并指出颗粒的棱角性越大，渗透系数越大^［6］。由正交试验分析表可知C_s1∶C_s2∶C_s3＝0.9∶1∶1.2，并且将试验数据进行回归分析，当形状系数C_s1＝0.18，C_s2＝0.2，C_s3＝0.24时与试验结果最为接近，此结论与卡门（Carmen）的研究成果^［7］相近。

3.4 土石混合体的渗透系数

由以上分析可知土石混合体的渗透系数与颗粒大小、颗粒级配、颗粒形状及孔隙比有关，同时渗透流体对渗透性也有一定的影响，主要是受液体的动力粘滞度η的影响，大量研究成果表明渗透系数k 与g/η 成正比^{［3，4，7］}。因此，土石混合体的渗透系数计算公式为

土石混合体

式中：k为土石混合体的渗透系数，cm/s；C_s为颗粒的形状系数，m^-3；d₂₀为等效粒径，小于该粒径的土重占总土重的20%，m；C_c为颗粒级配曲率系数，

；e为孔隙比；g为重力加速度，9.8 N；η 为液体的动力粘滞度，kPa · s（10^-6），η₂₀＝1.01×10^-6kPa·s。

由公式（1）计算出20℃时土石混合体的渗透系数k₂₀列于表7。与其他物理力学参数相比，土石混合体的渗透性变化范围要大得多。同时，受宏观构造和微观结构复杂性的影响，其渗透性具有高度的不均匀性^［8］。为进一步验证公式（1）的正确性，将实测值与由公式（1）得出的计算值进行对比分析，见图8。由图8可知由公式（1）计算出的渗透系数值与实测值基本吻合，9组试样的平均相对误差为21%，这对于离散性很强的土石混合体的渗透系数来说已经具有足够的精确性。

表7 计算值与实测值对应关系

图8 计算值与实测值关系

4 结论

（1）通过正交试验获取了砾石含量、孔隙比和颗粒形状对土石混合体渗透系数影响的主次顺序，并得出各因素的显着性水平，工程设计中可以通过合理调整土石混合体的砾石含量、孔隙比（压实度）和颗粒形状，以达到控制其渗透能力的目的。

（2）土石混合体的渗透系数与等效粒径d₂₀和曲率系数C_c组成的函数

成正比，并与孔隙比函数

成正比。

（3）提出了土石混合体渗透系数的计算公式，并通过试验结果验证了计算公式的正确性，为土石混合体渗透系数的定量预测提供了一个简明有用的计算工具。

参考文献

［1］油新华.土石混合体随机结构模型及其应用研究.北方交通大学博士论文，2001：1～18

［2］油新华，汤劲松.土石混合体野外水平推剪试验研究.岩石力学与工程学报，2002，21（10）：1537～1540，60～129

［3］刘杰.土的渗透稳定与渗流控制.北京：水利电力出版社，1992：1～20

［4］薛定谔A E.多孔介质中的渗流物理.北京：石油工业出版社，1984：141～173

［5］中华人民共和国水利部.土工试验规程（SL237—1999）.北京：中国水利水电出版社，1999：114～120

［6］ Tickell FG，Hiatt WN.Effect of angularity of grains on porosity and permeability of unconsolidated sands.AAPG Bulletin，1938，22（9）：1272～1274

［7］黄文熙.土的工程性质.北京：水利电力出版社，1984：60～129

［8］邱贤德，阎宗岭，刘立等.堆石体粒径特征对其渗透性的影响.岩土力学，2004，25（6）：950～954

C. 计算方法、步骤

（一）建立水文地质概念模型

解析法对水文地质条件限制较多，有严格的理想化要求，而实际水文地质条件往往十分复杂，为了能够用解析法计算，必须对水文地质条件进行合理的简化和概化，经过简化和概化后的水文地质条件称水文地质概念模型，它是对地下水系统的定性描述。

1.分析疏干流场的水力特征

矿床的疏干流场，是在天然流场背景下，叠加人为开采因素演变而成的，因此分析疏干流场各种水力特征时，均应以天然条件为基础，充分考虑开采的影响。

（1）区分非稳定流与稳定流

一般，疏干排水时，矿区地下水多为非稳定状态，但当疏干排水量小于地下水补给量时，可出现稳定状态。

矿山开采初期（开拓阶段），开拓井巷不断发展变化，疏干漏斗的外边界不断扩展，矿坑涌水量以消耗含水层储存量为主，该阶段疏干场一般为非稳定流，矿山开采后期（回采阶段），疏干流量主要受流场外边界的补给条件所控制，在补给条件不充分的矿区，疏干流场以消耗含水层储存量为主，仍为非稳定流，在补给条件充足的矿区，或具定水头补给边界的矿区，矿坑涌水量（或疏干量）被补给量平衡，一般出现相对的稳定流，矿坑涌水量预测可以稳定井流理论为基础。

（2）区分层流与紊流

矿区地下水在疏干条件下与天然运动状态相比，在大面积内仍为层流，仅在疏干工程附近常出现紊流，故达西定律（直线渗透定律）仍然是建立确定性模型的基础。

一般，常以抽（放）水试验为依据，用单位涌水量（q_i）法对层流、紊流进行判别，计算式为：

承压水

图13-7 水位降深为S_k的Q-t曲线

D. 地下水流基本问题的计算原理

一、地下水向河渠的运动

（一）河渠间地下水的稳定运动

1.潜水的稳定运动

河渠间潜水的运动由于受大气降水入渗补给或蒸发消耗的影响，应该属于非稳定运动。但为了简化计算，当入渗在时间和空间分布上都较均匀时，可以把潜水运动看作稳定运动。

假设条件：

（1）含水层为均质各向同性，隔水底板水平；

（2）河渠间潜水有垂向均匀入渗补给或蒸发消耗，设其强度W为常数；

（3）河渠基本上平行展布，潜水流可视为一维、渐变流并趋于稳定。

基于上述假设条件，取垂直于河渠的单位宽度进行研究并按图4-24取坐标，根据式（4-62），可写出上述问题的数学模型为：

图4-24 河渠间潜水的运动

续表

E. 层状岩层的等效渗透系数

在自然界中很常见的非均质岩层多是由许多透水性各不相同的薄层相互交替组成的层状岩层。每一单层的厚度比其延伸长度小得多（图1—20）。其平行于层面的渗透系数K_p和垂直于层面的渗透系数K_v不等。当每一分层的渗透系数K_i和厚度M_i已知时，可求出K_p和K_v。当水流平行于层面时（图1—20），通过层状含水层的总的单宽流量q等于各分层的单宽流量之和，总厚度M等于各分层厚度之和。对于每一分层而言，水力坡度J均为△H/l。因此，每一分层的流量为：

地下水动力学（第二版）

如果我们用一等效的均质含水层代替层状岩层，显然等效层的厚度等于层状岩层的总厚度，并且在同一水力坡度△H/l作用下应当有相同的流量q，因而有：

地下水动力学（第二版）

由此得

地下水动力学（第二版）

因而求得平行层面方向的等效渗透系数为：

地下水动力学（第二版）

等效导水系数为：

地下水动力学（第二版）

类似地，如果渗透系数在垂直方向变化，且没有明显的分层界线，而是逐渐连续过渡的，即：

K=K（z）

地下水动力学（第二版）

下面考虑水流方向垂直于岩层层面的情况（图1—21）。该情况下通过各分层的流量相同，即：

q₁=q₂＝…＝q_i＝…＝q_n=q

但水头降落和水力坡度不同，总的水头降落△H等于各分层水头降落△H_i之和。因此，对每一层都有：

地下水动力学（第二版）

用类似方法可得垂直于层面方向的等效渗透系数为：

地下水动力学（第二版）

由（1—46）式可以发现一个有趣的现象：垂直于层面的等效渗透系数主要取决于渗透系数最小，即阻力最大的分层。如有一层K_i=0，为不透水层，则K_v=0。

图1—21层状岩层中垂直于层面的渗流（据J.Bear修改）

（a）—水在承压含水层中流动；（b）—通过层状土的垂直下渗，水在一个大气压（101325Pa）下流出

平行层面的等效渗透系数K_p总是大于垂直于层面的等效渗透系数K_v。

思考题：

1.假设有一层状岩层，由三层组成，K₁=K₃。当K₂大于K₁和K₃，水流斜向穿过该层状岩层时，应如何折射？当K₂小于K₁和K₃时，应如何折射？画出示意图。

2.试自行证明，由二层组成的层状岩层，K_p＞K_v。

F. 作者———陈崇希教授简介

陈崇希 教授，博士生导师，男，1933年10月30日出生，浙江温州人，1956年毕业于北京地质学院水文地质及工程地质专业，中国地质大学（武汉）环境地质研究所前所长。

从事渗流理论，地下水数值模拟，地下水资源-环境评价与管理等方面的教学与科研工作。主持国家重点科技攻关项目专题（3项），国家自然科学基金（3项）等科研项目38项，以第一作者身份获奖的有：国家科技进步三等奖1项；省部级科技进步二等奖3项；部优秀教材二等奖1项。发表着作12部，论文100篇。是中国地质大学第一、二期211工程建设“地质环境保护及地质灾害防治”和“地下水资源与地质环境保护”重点学科群的学术带头人及首席科学家。

陈崇希教授主要学术成果简介如下：

1.地下水资源评价理论方面

地下水可持续开采量的评价是水文学与水资源工程等专业最重要的研究方向之一。然而长期以来，国内外多以地下水补给量为其评价基础。1966～1978年陈崇希教授论证：“井孔中抽水引起的影响范围随时间延续是发展的，因为井孔抽水破坏了原有的水动力平衡。假如没有因抽水而得到补给量的增加或排泄量的减少，那么影响范围将无限制地发展，”此时抽取的全部地下水为储存量；“只有当地下水抽水得到地下水的补给增量和排泄减量之和等于抽水量时，地下水才能形成稳定井流”。这是评价地下水可持续开采量的基本准则。2002年美国着名学者Bredehoeft证明的“可持续开采量”评价准则与陈崇希教授1966年论证的结论完全相同。多年来，“地下水的补给增量和排泄减量”的评价准则逐渐得到推崇，尚未见到非议。

2.地下水动力学解析理论方面

（1）1966～1981：纠正了“影响半径稳定井流”模型的错误，恢复裘布依“圆岛模型”的原貌。Dupuit模型是地下水稳定井流最基本的模型，20世纪50年代该模型由前苏联以“影响半径”模型传入我国，作为地下水资源评价的主要基础沿用了数十年。后来了解到，当时欧美同样也有将Dupuit模型刻画为“影响半径”模型，直至1972年，着名学者Bear J在其名着《多孔介质流体动力学》中也还是如此提法。

1966年，陈崇希教授对“影响半径稳定井流”模型提出质疑，从理论上证明了该模型不仅在原假定的初始水头水平、无入渗补给（地下水静止）条件下不可能形成稳定井流，而且在初始地下水流动（存在补给的径流场）条件下开采地下水也不都能形成稳定井流。陈崇希先生依据水均衡原理严格证明了形成稳定井流的条件。

1974年3月，陈崇希教授又进一步指出：Dupuit稳定井流方程的应用条件应是“园岛模型”，而非“影响半径”模型，Dupuit公式中的“影响半径”应改为“园岛半径”。这就从理论上纠正了20世纪50年代由国外传入、在我国影响深远的“影响半径模型”的错误。

陈崇希教授的科学见解在当时迅速得到响应。1975年，陈崇希先生应邀为中国地质科学院水文地质工程地质研究所以及来自国内各省市的同行阐述对Dupuit稳定井流模型的认识、系统介绍地下水非稳定井流的理论和方法。为此，陈崇希教授撰写了我国最早的一本系统阐述地下水不稳定井流理论的着作———《地下水不稳定井流计算方法》。

此后几年间，陈崇希教授跑遍国内各大城市图书馆、情报所寻找Dupuit的原着，未果。直至1981年，托人从法国巴黎图书馆复印到1863年Dupuit的原着，才见“庐山真面目”。原着所表达的正是陈崇希教授1974年提出的“圆岛模型”。陈崇希教授十多年执着追踪这一问题，不仅纠正了水文地质文献中长期的误传误用，证明了“影响半径稳定井流”模型的理论错误，而且在此基础上建立了地下水开采条件下的质量守衡方程，为正确评价地下水可持续开采量确定了基本准则。

（2）拓宽Theis公式和Hantush公式等的应用条件。Theis模型和越流系统的Hantush模型等模型都假定初始水头水平分布，此条件自然界很难满足。1975年，陈教授证明Theis降深公式和Hantush降深公式可拓宽用于初始稳定的流场；对于不稳定的初始流场，必需经过天然水头动态校正之后才能使用。Theis模型要求井径趋于零，1975年证明：对于实际的有限井径，只要经过很短的抽水时间即可使用Theis公式。2002和2003年，与研究生合作提出了抽水导致承压含水层顶板变形的修正Theis井流模型和Hantush井流模型，并得到解析解，从理论上证明该两模型抽水初期水头降深偏小的现象。

（3）1966年建立了孔隙介质和裂隙介质渗透系数的水力学模型/方程。该成果揭示了水文地质学最主要参数———渗透系数的物理实质。1972年Bear提出了与陈崇希先生相同的多孔介质渗透系数水力学方程，1967年法国学者Louis也提出了与陈先生类似的裂隙介质水力学模型。

（4）1974年，提出考虑含水层导水系数随时间变化的地下水承压-无压不稳定井流的解析解，与Moench（1972）把热传导中导管周围的冻结或融解问题的解直接移植用于承压-无压不稳定井流问题相比，各自有优缺点。2006～2008年又与研究生合作在该领域获得新成果。

（5）在岩溶泉流量的衰减方面，国际上多用含有衰减系数α的负指数函数进行分析和预测，但该系数一直作为经验参数来使用。1988年，林敏、陈崇希两教授合作，获得层流和紊流状态衰减系数的水力学方程，揭开了该系数与含水层基本参数之间的关系。

（6）2011年，建立了混合井流的某些基本解析解，获得若干重要认识。

（7）把不稳定流解析法与数值法的优势相结合———提出“数值-解析法”，解决了毛里塔尼亚首都伊迪尼水源地由于含水层存在岩性“天窗”获得补给造成纯解析法预测地下水开采动态的困难（1978～1980），数值模型中边界处理的困难（如韩城水源地，1989）以及传统有限元法和有限差分法，对抽水井水位的模拟存在明显的理论缺陷问题（1989）。

3.地下水流数值模拟方面

（1）改进滨海地区含水层系统的模型，提出确定海向边界的理论和方法。在众多的滨海地区的数值模型中，国内研究者都把海岸线作为直立边界；国外大部分同上，个别研究者无依据地取某个延伸距离为其边界。1986～1988年，陈崇希教授认为绝大多数含水系统向海底延伸，对不同排泄类型概括为“等效排泄边界”的概念，同时提出运用地下水的潮汐效应信息确定滨海含水层的“等效排泄边界”和水文地质参数的方法，并分别用于北海市几个水源地、海南岛洋浦港地区的地下水可持续开采量评价和烟台市的海水入侵模型等实例。这一系统而完整的确定海底边界的理论和方法在建立滨海水文地质模型上具有重要的理论意义和实践价值。（获国家科技进步三等奖）

（2）把传统的基于线汇的井孔-含水系统模型提升为“渗流-管流耦合模型”，解决了长期来井孔边界刻划的水文地质难题。为解决混合抽水井的模拟问题，陈崇希教授提出“渗流-管流耦合模型”。该模型己成功地解决北海市两层混合抽水试验确定分层水文地质参数（1992）、郑州市三层混合抽水试验确定分层水文地质参数（1995）以及我国西北地区4项地下水混合井开采动态预测（1999～2005）等实例。新模型的提出及混合井实例的应用，为“井孔-含水系统模型”开创了高一级的平台。

1995年，进一步把“渗流-管流耦合模型”发展为“岩溶管道-裂隙-孔隙”三重介质的地下水“线性非线性流”模型，突破了国外双重介质模型和线性管流模型的局限性。该模型成功地用于模拟广西环江岩溶强发育地区十分复杂的泉流量动态。

2003～2005年又将“渗流-管流耦合模型”用于自流井和水平井，将原模型提高了一个台阶。2003年，又向观测孔水位形成的传统观念提出质疑。国外着名学者（Hantush，Бочевер和Neuman等）都认为，常规观测孔中的水头降深可视为滤管内各点降深的平均值。然而，陈崇希教授指出：“这些计算观测孔水位降深的方法是缺乏物理基础的纯数学方法”，并用“渗流-管流耦合模型”从机理上仿真了观测孔水位的形成。这一问题的解决具有重要的理论意义和实用价值。

2004年完成了“渗流-管流耦合模型”的砂槽物理模拟，以进一步检验之。结果表明，数值模拟的水头、流量动态相当好地再现了物理模拟结果。为此，陈崇希教授从提出“渗流-管流耦合模型”和“等效渗透系数”的理论，物理模拟的验证以及不同条件下应用中的进一步检验等各环节，前后花去了10多年的时间。

1992年陈教授提出“渗流-管流耦合模型”，克服了传统基于线汇的井孔-含水系统模型的不实，将井孔的边界从滤管壁面处移至井口，从而无需人为给定滤管壁面处的边界条件，只需真实地给定井口处的流量或水头。该成果将长期以来水文地质界惯用的井孔-含水系统的模型提高到一个新的水平。陈教授卓有成效的工作，使我国相关领域的科学研究走在世界前列。

（3）“九五”国家科技攻关：发展地下水开采-地面沉降模型（以苏州地区为例）。土层固结与地下水流如何耦合？陈教授批评了国际上的“分两步走耦合”的错误，提出单步求解方法的理论依据；国际上多采用准三维流模型，陈教授指出“千层饼状”弱透水层的多层含水系统不宜用准三维流模型（其误差不是小于5％，而是超过30％），应采用真三维模型；建立了考虑土层固结引起孔隙度、渗透系数和给/储水系数变化的水流模型；刻画了地面沉降滞后于地下水开采层水头动态的普遍规律，也客观地仿真了苏州市地面沉降中心与地下水开采漏斗不一致的事实。

（4）2003年提出：“防止模拟失真，提高仿真性是数值模拟的核心”。在数值模拟上，陈崇希教授强调正确分析水文地质条件，研究水流机理和提高仿真技术，完成了30多个实际模拟模型，从科学性和实用性两方面完善地下水数值模拟技术。2005年与研究生一起，开发了总体上优于MODFLOW的PGMS软件，该软件含混合抽水井、常规观测孔、自流井、水平井、降雨/渠道等入渗补给地下水的滞后性、抽水井中水位、泉流量、地下水的非线性蒸发量、地下水-地表水相互作用等模块。

4.结语

回顾陈崇希教授从业的50年，处处可以感触到他发现科学问题的敏锐性、解决问题的智慧和敢于纠正前人错误、直面真理的勇气。陈崇希教授一向推崇学术争鸣，认为只有积极开展学术争鸣和交流，才能更好地促进学科的发展。陈崇希先生始终如一地坚持“求真、求实”的学术思想，锲而不舍、精益求精地寻求科学问题的解答。他时常教导学生“不仅要学会解决问题的方法，而且要培养发现问题和科学地提出问题的能力”。

陈崇希先生性格开朗、襟怀坦白，他对科学的热爱之情、对工作一丝不苟的科学家风范和对难题的攻关精神都深深地影响着他的学生和工作伙伴。

（摘引自《陈崇希教授的学术思想和成就综述》（焦赳赳等，2003，略加补充）

G. 等效渗透系数的名词解释

天然沉积土往往是由渗透性不同的土层所组成。 对于与土层层面平行和垂直的简单渗流情况，当各土层的渗透系数和厚度为已知时，我们可求出整个土层与层面平行和垂直的平均渗透系数，作为渗流计算的依据。 需要对水平渗流情形和竖直渗流情形分别单独考虑。

H. 土石混合体渗透性能的试验研究

周中¹ 傅鹤林¹ 刘宝琛¹ 谭捍华² 龙万学² 罗强²

（1.中南大学土木建筑学院 湖南 长沙 410075

2.贵州省交通规划勘察设计研究院 贵州 贵阳 550001）

摘要 土石混合体属于典型的多孔介质，其渗透特性与砾石的百分含量关系密切。通过自制的常水头渗透仪，测定了不同含砾量时土石混合体渗透系数值，研究发现含砾量与土石混合体渗透系数之间存在指数关系；基于幂平均法，提出了土石混合体复合渗透系数的计算公式，并通过试验结果验证了该式的正确性，为土石混合体渗透系数的理论计算提供了一个简明有用的计算工具。

关键词 土石混合体 多孔介质 渗透性能 复合渗透系数 经验公式

土石混合体一般是由作为骨料的砾石或块石与作为充填料的粘土或砂组成，它是介于土体与岩体之间的一种特殊的地质体，是土和石块的介质耦合体^［1］。因为土石混合体具有物质组成的复杂性、结构分布的不规则性以及试样的难以采集性等特殊的性质，从而给研究带来极大的困难，目前人们对于它的研究仍处于探索之中^［2］。渗透与强度和变形特性都是土力学中所要研究的主要力学性质，其在土木工程的各个领域都有重要的作用^［3］。土石混合体属于典型的非均质多孔介质^［4］，其渗透系数是由高渗透性的砾石和低渗透性的土体复合而成的。土的渗透系数可以通过室内试验由达西定理计算得出，然而土石混合体的渗透系数却难以确定，主要原因是：取样困难；难以进行常规的渗透试验；大尺度的渗透试验不仅造价高、准确性差，而且试验结果离散度大，难以掌握其规律性。因此能够求出土石混合体复合渗透系数的计算公式具有重要的理论意义和工程应用价值。

土石混合体中土与砾石粒径的界限值为5mm，即将粒径小于5mm的颗粒称为土、大于5mm的颗粒称为石，砾石含量用P₅表示^［1］。利用自制的常水头渗透仪，研究砾石体积百分含量P₅从0%逐步过渡到100%（间隔10%）时土石混合体的渗透系数，每种配比作平行试验3次，共33次渗透试验。

1 土石混合体渗透性能试验

1.1 试样的基本物理力学性质

试验所取土样为正在修建的上瑞高速公路贵州段晴隆隧道出口处典型性土石混合体，其天然状态土的物理指标及颗粒级配曲线见表1和图1。由图1可知现场取回土样的不均匀系数C_u为12.31，说明土样中包含的粒径级数较多，粗细粒径之间差别较大，颗粒级配曲线的曲率系数C_c为1.59，级配优良。

表1 天然状态土的基本物理指标

图1 天然状态土的颗粒级配曲线

1.2 大型渗透仪的研制

《土工试验规程》（SL237—1999）规定粗粒土的室内渗透系数需由常水头渗透仪测试，国内常用的常水头渗透仪是70型渗透仪。70型渗透仪的筒身内径为9.44cm，试验材料的最大粒径为2cm，规范^［5］要求筒身内径应为最大粒径的8～10倍，因此70型渗透仪的筒身内径过小，有必要研制大尺寸的渗透仪。自制渗透仪的内径和试样高度至少应为最大颗粒粒径的8倍，即至少应为16cm，另外，考虑到边界效应，试样的上下两头分别增加2cm，因此，自制渗透仪的内径和试样高分别取为16cm和20cm。考虑到土石混合体的渗透性较强，选取进排水管的口径为2cm。自制的大型常水头渗透仪如图2和图3所示。

图2 自行研制的渗透仪

图3 常水头渗透仪示意图

数据单位为cm

1.3 试验步骤

首先，将由现场取回的土样烘干、过筛，并根据粒径的大小分为0～5 mm的土和5～20mm的砾石两部分。然后，按照试验要求的砾石体积百分含量P₅，以10%的初始含水量配制试样，静置24 h。试验时，将配制好的试样分层装入圆桶中，每层装料厚度30mm左右，分层压实，记录每层的击实数。按上述步骤逐层装样，至试样顶部高出测压孔约3cm为止。测出装样高度，准确至0.1cm。在试样顶部铺一层2cm厚的细砾石作缓冲层。之后，由进水管注入蒸馏水，直至出水孔有水流出，静置24 h使试样充分饱和。用量筒从渗透水出口测定渗透量，同时用温度计测量水温，用秒表测记经一定时间的渗水量，共测读6次，取其平均值，6次结果相差不得超过7%，否则需重新测定。

1.4 试验数据

按照试验设计的各种砾石体积百分含量P₅共需作11组试验，每组试验作平行试验三次，取3次测量的平均值，并乘以温度校正系数

，即可求出每组试验20℃时的渗透系数，渗透系数的测量结果见表2。

表2 渗透系数测定结果

2 试验结果分析

2.1 渗透系数与砾石含量的关系

不同含砾量的颗粒级配曲线如图4所示，由图4可以求出各曲线的粒径特征系数及不均匀系数C_u和曲率系数C_c。

图4 试样的颗粒级配曲线

图5为土石混合体砾石含量P₅与20℃时渗透系数的关系曲线。从图5可以看出，随着含砾量的增加，渗透系数急剧增加，可见，在设计中可以通过调节砾石的含量来控制土石混合体的宏观渗透性能。

图5 粗粒含量与渗透系数的关系

从图5还可以发现，土石混合体中砾石的含量P₅与渗透系数k之间存在指数关系，与文献［6］的研究成果相似，即

土石混合体

式中：k₀为P₅＝0时土的初始渗透系数；n为与土石混合体本身性质相关的常数。对于文中试验值，k₀与n分别为0.0006cm/s和8.82。在工程中可以通过少量试验来确定k₀，n值，以此来预测不同级配土石混合体的渗透性。

2.2 土石混合体的复合渗透系数

近几十年来，许多学者在揭示影响和决定土的渗透系数内在因素及其相互关系方面进行了大量工作，并取得了有益的成果^［7～12］，被认为依然有效且目前常用的确定渗透系数的半经验、半理论公式有：

（1）水利水电科学研究院公式^［7］：

土石混合体

式中：k₁₀，k₂₀分别为温度为10℃和20℃时的渗透系数（cm/s）；η₁₀/η₂₀为温度为10℃和20℃的粘滞系数比；n为孔隙率；d₂₀为等效粒径（mm）。

（2）泰勒（Taylor）^［9］用毛管流的哈根-伯努力（Hange-Poiseuille）方程导出渗透系数的表达式：

土石混合体

式中：d_s为当量圆球直径，可以用等效粒径d₂₀代替；γ_w为液体容重；μ为液体粘滞度；e为孔隙比；C为形状系数，通常取C＝0.2。

式（2）和式（3）均是针对土体的渗透特性提出的半经验、半理论公式，然而对于非均质性更强、粒径差别更大的土石混合体来说，其适用性不是很强。土石混合体中砾石形成骨架，细颗粒充填孔隙，其渗透系数是由低渗透介质土体的渗透系数k_S和高渗透性介质砾石的渗透系数k_G复合而成。土石混合体复合渗透系数不是按体积百分含量的简单复合，而是高低渗透性介质的耦合。在参考相关文献^{［10～12］}的基础上，基于幂平均法，本文提出的土石混合体复合渗透系数k复合的表达式为

土石混合体

式中：P₅为砾石的体积百分含量，%；k_G为砾石的渗透系数，cm/s；k_S为土的渗透系数，cm/s；f为系数。

砾石的体积百分含量P₅可以由筛分法求出；土的渗透系数k_S和砾石的渗透系数k_G可以由室内试验直接求出或参考相关资料确定；系数f可以通过少量试验回归分析确定，因此可以说（4）式是一个简明实用的土石混合体复合渗透系数计算公式。

图6 不同计算方法结果比较

为进一步验证（4）式，我们将试验测得的k值与用（2），（3），（4）式计算得到的k值进行对比分析。结果见图6，具体数值见表3。由图6和表3可知据水利水电科学研究院公式和泰勒公式计算结果均高于实测值，尤其是当P₅≤30%时，（2）式计算结果和（3）式计算结果比实测值大2～3个数量级，与实测值相差较大。而用本文方法得到的土石混合体的渗透系数最接近实测值，平均相对误差仅为0.6%，能够作为土石混合体渗透系数定量预测的有效工具。在工程设计中，可以根据工程对土石混合体渗透性的要求，依据本文提供的经验公式，调整土石混合体中砾石的含量，达到控制土石混合体渗透能力的目的。

表3 土石混合体渗透系数及相关参数

3 结论

（1）利用自制的常水头渗透仪，测定了不同含砾量时土石混合体的渗透系数值，并指出含砾量与土石混合体渗透系数之间存在指数关系。在工程设计中可以通过合理调整土石混合体中砾石的含量，达到控制其渗透性能的目的。

（2）指出土石混合体的渗透系数是一种由高渗透性的砾石和低渗透性的土体复合而成的，给出了土石混合体复合渗透系数的计算公式，并通过试验结果验证了计算公式的正确性，为土石混合体渗透系数的定量预测提供了一个简明有用的计算工具。

参考文献

［1］油新华.土石混合体随机结构模型及其应用研究.北方交通大学博士论文.2001：1～18

［2］油新华，汤劲松.土石混合体野外水平推剪试验研究.岩石力学与工程学报.2002，21（10）：1537～1540，60～129

［3］黄文熙.土的工程性质.北京：水利电力出版社.1984：60～129

［4］薛定谔 A E.多孔介质中的渗流物理.北京：石油工业出版社.1984：141～173

［5］中华人民共和国水利部.土工试验规程（SL237—1999）.北京：中国水利水电出版社，1999：114～120

［6］邱贤德，阎宗岭，刘立等.堆石体粒径特征对其渗透性的影响.岩土力学，2004，25（6）：950～954

［7］刘杰.土的渗透稳定与渗流控制.北京：水利电力出版社，1992：1～20

［8］ Wen X H，Gomez-Hernandez J J.Upscaling hydraulic conctivities in heterogeneous media：An overview.Journal of Hydrology，1996，183：ix～xxxii

［9］ Taylor D W.Fundamentals of soil mechanics.John Wiley ＆ SONS，Inc.，1948

［10］ Brown W F.Solid mixture permitivities.Journal of Chemical Physic，1955，23（8）：1514～1517

［11］ Dagan G.Analysis of flow through heterogeneous random aquifers by the method of embedding matrix—1：Steady flow.Water Resources Research，1981，17（1）：107～122

［12］ Noetinger B.The effective permeability of a heterogeneous porous medium.Transport in Porous Media，1994，15：99～127