某污水池下游300m处有潜水监测井。该潜水含水层渗透系数为100m/d，有效孔隙度为25%，水力坡度为0.5%。若污水池发生泄漏，污水进入含水层后水平运移到达监测井处的时间约为（　）。

A、3D

B、150D

C、300D

D、600D

参考答案:

【正确答案：B】

根据达西定律，地下水渗流速度（u）=渗透系数（K）×水力坡度（i）=100×0.5％=0.5（m/d）。通过孔隙断面的水质点的实际平均流速（u）=0.5/砂的孔隙度（n）=0.5/25％=2（m/d）。则污水进入含水层后水平运移到达监测井处的时间=300／2=150（d）。

水量均衡法

（一）基本原理

水量均衡法是根据水量平衡原理，建立均衡方程计算水量的方法，表达式为

∑Q补-∑Q排=ΔQ储 （3-1）

式中：∑Q补为均衡期内地下水系统各种补给量的总和（m3）；∑Q排为均衡期内地下水系统各种排泄量的总和（m3）；ΔQ储为均衡期内地下水系统内部储存资源的变化量（m3）。

（二）一般步骤

1.确定均衡区

根据地下水系统理论的要求，均衡区应是地下水系统边界所界定的空间范围，一般要求以地下水系统天然边界作为划分依据。由于水量均衡法属于集中参数系统，为了提高区域地下水数量评价精度，在实际计算时可以根据不同水文地质条件划分为不同级别的子区，分别计算各均衡要素，然后进行综合。例如根据给水度、降水入渗系数、地下水埋藏深度等条件，将均衡区划分为若干子区，分别计算各子区的储变量、降水入渗量和潜水蒸发蒸腾量，然后求和。

2.确定均衡要素

确定式（3-1）中∑Q补和∑Q排的组成，即确定地下水系统三维空间区域边界上的输入和输出量。从外界进入地下水系统的各种水量统称为补给项，系统输出的各种水量统称为排泄项。

一般而言，补给项包括：大气降水入渗补给量、地表水体渗漏补给量（河流、湖泊、水库等）、地下侧向流入补给量、越流补给量、凝结水补给量、地表水灌溉入渗补给量、地下水灌溉回归补给量、渠系渗漏补给量、人工回灌补给量等。

排泄项包括：潜水蒸发蒸腾量、地下水侧向流出量、地下水开采量、泉水溢出量、越流排泄量、向河湖排泄量等。

需要指出的是，不同的地下水系统与外部环境之间的水量交换关系不同，所以均衡要素的组成因不同地下水系统而异。在实际工作中，需要与研究区具体条件紧密结合，确定均衡要素的组成。

3.确定均衡期

地下水均衡计算是针对某一特定时间段进行的，称为均衡期。如前所述，在地下水的资源功能评价中，要求地下水数量评价的时间尺度为5～12年，以此为均衡期进行水量均衡计算。为保证水量平衡，各均衡要素计算和相关的资料的选取应采用统一的时间序列。

（三）均衡项计算方法

1.降水入渗补给量

降水入渗补给量确定方法包括：直接测定法、零通量面法、包气带达西定律法、氯质量平衡法、示踪法等。

（1）直接测定法

通常利用测渗仪或通过包气带蒸渗试验直接测定不同岩性、不同地表覆盖情况下的降水入渗补给量（Young et al.，1996）。我国于20世纪70年代末期开始在华北地区和西北地区建立了许多包气带试验场，开展了大量的实验研究。

（2）零通量面法（ZFP）

零通量面是Richards于1956年首先提出的，是包气带水分运移的分界面，其上土壤水分向地表运移，其下水分向地下水运移，将该面以下的水分运移速率作为地下水补给速率，利用该法需要测定包气带垂直剖面土壤水势和含水量。我国于20世纪80年代后期引入ZFP法（张光辉，1988；张惠昌，1988），目前该方法仍在应用（程辉等，2000；周金龙等，2003；李茜等，2006）。

（3）包气带达西定律法

达西定律法是干旱、半干旱地区常用的方法，需要测定包气带水力梯度和不同含水量下的渗透系数，计算公式如下：

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/d）；K（θ）为包气带水渗透系数（m/d）；H为包气带水侧压水头（m）；Z为垂直位置高程（m）。

（4）氯质量平衡法

该方法主要应用了氯的化学稳定性，其应用前提是（Kinzelbach，2002）：①由于包气带溶质输入和向饱水带的输出存在时间滞后，所以必须假定在此期间没有重要的气候变化；

②没有额外的溶质加入，如肥料，同时也没有近期大气污染；

③在ZFP之上和之下没有溶质储存的净变化，这种变化可能由于动植物引起或矿物沉淀/溶解和吸附/解吸附。在满足以上条件的基础上，可采用下式计算降水入渗补给速率：

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/a）；P为多年平均降水量（m/a）；CP为降水中氯离子浓度（mg/L）；Fd为氯离子干沉降量（mg/m2·d）；CS为零通量面以下包气带水的氯离子浓度（mg/L）。

（5）示踪法

利用人工或环境示踪剂，通过失踪剂峰面移动来计算降水入渗补给速率。常用的人工示踪剂包括：氚、溴、碘、染色剂等（Athavale，1988；Kung，1990；Flury，1994；Aeby，1998；Forrer，1999），环境示踪剂包括氚、氯-36 等受核爆影响的放射性同位素（Scanlon，2002）。降水入渗速率计算公式为

区域地下水功能可持续性评价理论与方法研究

式中：q为降水入渗补给速率（m/a）；Δz为示踪剂浓度峰面的运移深度（m）；Δt为峰面运移时间（a）；θ为体积含水率（无量纲）。

由于受地形地貌、地表覆盖、包气带岩性及厚度、降水强度及频率、包气带水分状况、地下水埋深等条件的影响，不同地带的降水入渗速率不同。而采用以上方法获得的数据仅是某一或某些条件下的实验结果，所代表的空间尺度有限，且不同的方法所代表的时间尺度也不相同（表3-2）。因此，在区域地下水资源评价中，往往根据研究区实际条件进行适当分区，选用不同方法求得不同分区的降水入渗系数（a），然后采用下式计算降水入渗补给量：

Qp=a·P·F （3-5）

式中：Qp为降水入渗补给量（m3/a）；a为降水入渗系数（无量纲）；P为多年平均降水量（m/a）；F为计算区面积（m2）。

表3-2 不同方法确定的降水入渗速率范围及时空尺度对比

注：表中数据根据Scanlon等（2002）整理。

目前，我国北方大部分地区已经通过包气带入渗试验、水位动态分析等方法，建立起不同地区降水入渗系数与地下水位埋深和包气带岩性之间的关系。

2.地下水与河流之间交换量

（1）断面流量差法

若均衡区有河流穿过，则在均衡区上、下游边界处各选一个测流断面监测流量，并确定断面之间的距离、测流时间间隔、河流水面宽度和水面蒸发量，然后采用以下公式计算：

Qr=（Q1-Q2）·Δt-B·L·E （3-6）

式中：Qr为测流期间河道渗漏补给量（m3）；Q1，Q2分别为河流上、下游断面的平均流量（m3/s）；Δt为计算时段（s）；B为河流水面平均宽度（m）；L为河流两断面间的距离（m）；E为测流期间的水面蒸发量（m）。

（2）渗流断面法（达西定律）

当河水与地下水有直接水力联系时，采用达西定律计算河道侧渗量，公式为

Qr=K·L·I·h·Δt （3-7）

式中：K为含水层渗透系数（m/d）；L为河道渗漏段长度（m）；I为河渠一侧地下水水力梯度（无量纲）；h为过水断面的厚度（m）；Δt为计算时段（d）。

h的取值应根据河流与地下水的关系而定。当河流一侧接受地下水补给，另一侧补给地下水时（图3-2a），h取值为河床到地下水位（河水位）的距离；当河流两侧都补给地下水时（图3-2b），h取值为含水层的整个厚度。

（3）基流分割法

在地下水补给常年性河流的地区，在枯水期河水流量几乎全部由地下水补给维持，这时的河水流量被称为基流量。把河流流量过程线上的基流量分割出来，即为地下水对河流的补给量（房佩贤等，1987；曲焕林等，1991；徐恒力等，2001）。图3-3为典型的单峰流量过程线，由起涨部分、峰值和退落部分组成。起涨部分的起点称为起涨点（图3-3中的a点），在退落部分，当降水影响消失时，河流量由地下径流组成，其起点称为地下水退水点（图3-3中的d点）。起涨点很好确定，而确定地下水退水点比较困难，一般有3种方法：经验法、退水曲线法和作图法。

图3-2 河流与地下水补排关系示意图

经验法，就是在过程线上的退落部分，找到曲线曲率最大的点即视为地下水退水点。

退水曲线法，认为从退水点开始，流量变化满足布西涅斯克方程（退水曲线方程）：

Q=Q0·e-kt （3-8）

式中：Q为从d点开始的任一时刻的河流量（m3/s）；Q0为d点的流量（m3/s）；k为衰减系数；t为以d点为起点的时间（s）。

由式（3-8）可知，从退水点开始流量呈等比级数递减。利用这个规律，在退落部分找到流量大体成等比递减的开始时刻，即为地下水退水点。

作图法，在退落部分按相等的时段，选取一系列流量，计算流量差（ΔQ），然后以ΔQ为纵轴，以时间（t）为横轴绘制曲线，将曲线的拐点所对应的时刻作为退水曲线的初始时刻，然后在图3-3中找到该时刻所对应的点即为退水点。

图3-3 河流流量过程线

由于河流与地下水之间的水动力关系不同，基流分割方法也不同。一般有两种情况（图3-4）：一种情况是河流与地下水有直接水力联系；另一种情况是二者之间不存在直接水力联系。

当河流与地下水无直接联系时（图3-4（a）），若不考虑地下径流峰值，则直接连接ad，其下方阴影部分的面积即为基流量（图3-5（a））；若考虑地下径流峰值，则可分别计算流量过程线起涨部分（af段）的平均流量（Q1）和退落部分（df段）的平均流量（Q2）。然后，再计算出起涨时段内大气降水形成的平均流量（Q′1）：

图3-4 河流与地下水关系示意图

区域地下水功能可持续性评价理论与方法研究

式中：P为af段的总降水量（m）；F为测站所控制的流域面积（m2）；taf为从a点到f点经历的时间（d）；v为径流系数，等于径流深度与降水深度之比。

同理，求出退落时段内大气降水形成的平均流量（Q′2）之后，在流量过程线上找到分别与 和 相对应的点，过两点作垂直t轴的直线，在两直线上分别减去Q′1和Q′2，可得到b，c两点（图3-5（b）），连接a，b，c，d，其下的阴影部分面积即为基流量。

图3-5 地下水与河流无直接水力联系时的基流分割图示

当河流与地下水有直接水力联系时（图3-4（b）），枯水季节地下水补给地表水，丰水季节河水位高于潜水水位，地表水补给地下水。此种情况下，在枯水季节（起涨点之前和退水点之后），河流流量全部为基流量。将起涨点对应的时间记为ta，退水点对应的时间记为td。进入洪水期后，河水开始补给地下水，但在ta之前进入河道的地下水与洪水一起从上游向下游流动。河流源头到测站的距离很容易测定，则可以求出河水从源头流到测站所用的时间（记为Δt），也就是说，到ta+Δt时刻河水全部由洪水组成，其对应的点记为b点，则地下水径流量应按 ab 逐渐减少；同样，在洪峰过后，河流源头首先有地下水进入河道，起始时刻为td-Δt，其对应的点记为c，则地下水径流量应按 cd 逐渐减少。基流分割如图3-6中的阴影部分。

图3-6 地下水与河流有直接水力联系时的基流分割图示

（4）示踪法

水中的氢氧稳定同位素常用来示踪地下水与地表水的相互交换量，通过河道水量均衡方程和质量平衡方程的联合求解，计算地下水与河流交换量（Scanlon，2002）。水量均衡方程和质量平衡方程如下：

Qup+∑Qin+Qgi=Qdown+∑Qout+Qgo+Er （3-10）

Qup·δup+∑Qin·δin+Qgi·δgi=Qdown·δdown+∑Qout·δout+Qgo·δgo+Er·δEr　（3-11）

式中：Qup，Qdown分别为上、下游断面河流量（m3/s）；Qin，Qout分别为测流断面间各支流的流入、流出量（m3/s）；Qgi，Qgo分别为测流断面间地下水流入、流出量（m3/s）；Er为测流断面间河道水蒸发量（m3/s）；δup，δdown分别为上、下游断面河水氢氧稳定同位素δ值（‰）；δin，δout分别为测流断面间各支流的流入和流出水的氢氧稳定同位素δ值（‰）；δgi，δgo分别为测流断面间地下水流入、流出水的氢氧稳定同位素δ值（‰）；δEr为河道蒸发水的氢氧稳定同位素δ值（‰）。

（3-11）式中的δEr通常难以测定，Krabbenhoft（1990）给出了计算公式：

区域地下水功能可持续性评价理论与方法研究

式中：δL，δα分别为地表水和大气水汽同位素含量；h为相对湿度；α′为水-气界面温度下同位素平衡分馏因子，等于1/α；ε为总分馏因子，ε=1000（1-α′）+Δε；Δε为动力学分馏因子，对于δD，Δε=12.5（1-h），对于δ18O，Δε=14.2（1-h）（Gonfiantini，1986）。

平原区的河流往往在上游地带渗漏补给地下水，在下游地带接受地下水补给，因此，在实际应用时，通常需要首先确定该两种情况发生的分界面，然后分段进行计算。

3.地下水侧向流入流出量

一般采用达西定律计算，公式为

Qg=v·B·M=K·J·B·M （3-13）

式中：Qg为地下水侧向流入流出量（m3/d）；v为地下水渗流速度（m/d）；B为过水断面宽度（m）；M为含水层厚度（m）；J为地下水水力梯度（无量纲）。

计算时，需要注意：①采用测流计或地下水示踪技术测定地下水流速时，所测定的流速是地下水的实际流速（va）。

②采用达西定律计算时，需要将实际流速换算成渗透速度，即v=va·n（n为有效孔隙度）。

③不同地段含水层的孔隙度和地下水水力梯度不同，不同时段地下水的水力梯度也不相同。因此，实际应用时，应分地段、分时段分别进行计算。

4.地下水与湖泊或水库的交换量

（1）水量均衡法

计算公式为

Qlr=Qgi-Qgo=ΔV-P·F+E·F-Qsi+Qso （3-14）

式中：Qlr为地下水与水库或湖泊的净交换量（m3/a），Qlr＞0，则地下水向湖泊排泄量大于湖泊向地下水的渗漏量，湖泊接受地下水的净补给；Qgi为地下水向湖泊的排泄量（m3/a）；Qgo为湖泊向地下水的渗漏量（m3/a）；ΔV为水体体积的年变化量（m3/a）；P为年降水量（m/a）；F为水体水面面积（m2）；E为水面蒸发量（m/a）。Qsi为地表水年流入量（m3/a）；Qso为地表水年流出量（m3/a）。

（2）示踪法

采用式（3-14），一般只能获得地下水与湖泊之间的净交换量，为了分别求取Qgi和Qgo，可以利用水体中的天然示踪剂建立质量均衡方程（Sacks，1998）：

Qgi·Cgi-Qgo·CL=ΔV·CL-P·F·CP+E·F·CE-Qsi·Csi+Qso·CL　（3-15）

式中：Cgi为补给湖泊地下水的示踪剂浓度；CL为湖泊水的示踪剂浓度；CP为降水的示踪剂浓度；CE为湖泊蒸发水的示踪剂浓度；Csi为补给湖泊地表水的示踪剂浓度；其他符号同式（3-14）。

联合求解（3-14）和（3-15）两个方程，即可获得Qgi和Qgo。常用的天然示踪剂是水中的氢氧稳定同位素（Sacks，1998；Scanlon，2002），这时需要知道蒸发水的氢氧稳定同位素值，其计算方法见式（3-12）。

5.越流量（包括补给和排泄）

计算公式为

Qy=F·K·J （3-16）

式中：Qy为越流量（m3/d）；F为计算面积（m2）；K为弱透水层的垂直渗透系数（m/d）；J为弱透水层上下含水层间的水力梯度（无量纲）。

6.凝结水补给量

可根据均衡试验场地中渗透仪的观测资料求得，但是计算时应注意将观测中冬季潜水冻结层融化的水量扣除。

7.渠系渗漏补给量

根据渠系衬砌状况，选用实测或经验系数计算。若渠道没有任何衬砌，其渗漏补给量与河道渗漏补给量计算方法相同。若渠道有衬砌，则可采用如下公式计算：

Qci=r·（1-η）Qc·Δt （3-17）

式中：Qci为渠道渗漏量（m3）；r为渠道渗漏修正系数（无量纲）；η为渠系有效利用系数（无量纲）；Qc为渠道过水量（m3/s）；Δt为计算时段（s）。

8.田间灌溉入渗补给量

（1）入渗系数法

计算公式为

Qsi=β·Qs·F·N （3-18）

式中：Qsi为田间灌溉入渗量（m3）；β为入渗系数（无量纲）；Qs为灌溉定额（m3/m2）；F为灌溉面积（m2）；N为灌溉次数。

（2）水量均衡法

根据水均衡原理，用灌溉量减去排放量、蒸发量和其他消耗量计算。

（3）地中渗透仪法

在田间专门设置地中渗透仪，直接测定灌溉水渗漏补给量。

9.潜水蒸发蒸腾量

（1）蒸发系数法

计算公式为

Qe=E·c·F （3-19）

式中：Qe为潜水蒸发蒸腾量（m3/a）；E为水面蒸发量（m/a）；c为潜水蒸发系数（无量纲）；F为计算面积（m2）。

（2）经验公式法

通常利用经验公式求出潜水蒸发强度（ε），然后按下式计算：

Qe=ε·F （3-20）

式中：ε为潜水蒸发强度（m/a）；F为计算面积（m2）。

潜水蒸发强度一般采用柯达夫-阿维利扬诺夫公式计算，即

区域地下水功能可持续性评价理论与方法研究

式中：λ为植被修正系数（无量纲）；h为潜水水位埋深（m）；h0为潜水蒸发的临界深度（m）；θ为无量纲指数，因气候和土壤而异，取值1～3，一般可以取1；其他符号同前。

10.储变量计算

计算公式为

ΔS=μ·F·ΔH （3-22）

式中：F为计算面积（m2）；ΔH为水头变化（m）；μ为给水度（潜水）或储水系数（承压水）。

（四）若干问题说明与适用条件

1.参数的获取

以上介绍的方法中涉及各种参数（渗透系数、导水系数、给水度、储水率、储水系数、孔隙度、垂向渗透系数、越流系数、降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等），按“全国地下水资源及其环境问题调查评价技术要求系列（一）”的要求获取。

2.参数分区

水量均衡法是一种集中参数系统的方法，在地下水数量评价时往往难以满足计算精度要求，尤其是在区域地下水数量评价中。地下水系统是一个复杂的非均质系统，各种参数（渗透系数、导水系数、给水度、储水率、储水系数、孔隙度、垂向渗透系数、越流系数、降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等）是空间位置的函数，有些还是时间的函数（如降水入渗系数、灌溉入渗系数、潜水蒸发系数、渠系渗漏系数等），所以为了提高计算精度，需要综合考虑各种参数的时空变化特征，在空间上将地下水系统划分成若干子区块，在时间上划分成若干时段，在子区块和各时段可以认为各种参数是一个相对稳定的数值，然后分别计算各个区块和时段的水量，最后集成总量。理论上讲，划分的区块越小、时段越多，计算精度就越高，但是工作量就越大。在实际应用时，应根据评价区所拥有的资料状况和计算精度要求进行适当的划分。

3.点参数的区域化

通过各种方法获得的参数，大多是点源数据。在区域地下水数量评价时，需要将点源数据转化成区域数据。点源数据的区域化，常采用的方法是Kriging插值法。Kriging法是一种最佳空间估计法，其本质是最佳无偏估计，是对空间分布的数据求线性最优、无偏内插估计的一种方法（Gress N.A.C.，1990，1991；Dentsch C.V.，1992；Gelhar L.W.，1993）。常用的软件中都有Kriging插值功能模块，如Surfer、MapGIS等，也有一些文献中给出了计算程序（徐士良，1995；Dentsch C.V.，1992）。

4.适用条件

水量均衡法方法原理明确，计算公式简单，计算精度高低可调，适应性强。但是，在补给、排泄条件复杂的地区，涉及的均衡要素较多，某些均衡要素难以准确测定或求取成本和工作量较大，计算精度不如数值法高。

水量均衡法既可用于区域地下水数量评价，也可用于局域地下水数量评价或水源地评价；既可评价地下水补给资源量，又可评价可采资源量，是最常用、最基本的地下水数量评价方法，其成果是其他方法的验证依据之一。在补给、排泄条件简单，地下水系统边界比较清楚，水均衡要素容易确定的地区，应用效果较好，评价结果精度高。

模型的特色及改进

一、混合井孔的模拟

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

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

图5-2 混合抽水井

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

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

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

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

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

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

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

4）含水层的导水系数一般不随时间变化，即导水系数比T1/T2不会变化，但流量比Q1/Q2却不是一个常量。

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

6）从另一角度分析，混合观测孔是混合抽水井的特殊情况（Qw=0），对于两层混合的观测孔，其孔中水位（混合水位）必介于两含水层水位之间，即混合观测孔对于其中一含水层（例如1含水层）起抽水作用（Q1＞0），对于另一含水层（2含水层）起注水作用（Q2＜0）。如此，Q1/Q2＜0。而两含水层的导水系数的比值肯定是正值，即T1/T2＞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/d2。

当管流为层流时

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

式中：Re为雷诺数。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

则

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

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

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

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

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

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

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

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

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

D.Sokol（1963）采用Thiem“影响半径”稳定井流模型证明混合观测孔的水位等于以导水系数为权重的各层水头的代数平均值。该论文存在如下几个主要问题：①Thiem“影响半径”模型是不能形成稳定流的（陈崇希，1966，1975，1983）；

②Sokol隐含着假定，尽管两含水层的导水系数不等，流量也不同，但取其“影响半径”相等；

③Sokol还隐含着假定，井管的阻力可以忽略不计，即井管处处水头相等。在这些假定条件下得出的结论是不可信的。

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

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

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

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

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

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

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

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

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

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

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

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

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

三、泉流量动态的模拟

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

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

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

②用承压完整井的公式计算泉的流量，特别是柳林泉的流量，不妥；

③比例系数q在不同时间（不同条件）并不是常量。

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

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

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

式中：Q为泉流量；Kzsp为泉流区域的垂直渗透系数；Hb为泉口下层结点的水头值；Zsp为泉出露标高；Dzsp为泉流过程中的渗流长度（含水层厚度）；S为泉口结点控制面积；CA为泉颈修正系数；Cs为泉径修正系数；Tij，Tik为流段ij和ik的平均导水系数；，，，为线段的长度；hi，hj，hk为结点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）作用，则该月及后各月降雨入渗补给量之和∑REk（mm）就是该时段降雨量形成的全部补给量。该值∑REk与0月的降雨量RA之比就是（总）入渗补给系数α。各月份的入渗量与0月降雨量RA之比，则是0月降雨在各月的入渗补给系数αk。即

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

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

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

式中：RA为0月的降雨量，mm；REk为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个模拟层，即便是同一岩性（如亚粘土）而处于不同埋藏深度，由于其固结程度不同而导致其水文地质参数不相等。这种情况下如何给定初始参数估计值的空间分布又是一个难题。若按常规惯用的方法，则同一岩性根据不同层位（不同埋藏深度）要给出不同参数分区，如此的参数分区数过多，给调参带来很大的工作量。

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

六、自流井的模拟

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

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

七、初始水头的形成

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

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