Formation Mechanism of Water Level and Its Determination Method in Conventional Observation Wells for Three-Dimensional Groundwater Flow
-
摘要: 对地下水三维流中常规观测孔中水位的传统计算方法提出质疑, 认为计算观测孔中水位的Hantush Бочевер方程是缺乏物理基础的纯数学方法.分析了形成观测孔中水位的机理, 提出三维地下水流中常规观测孔中只是孔口的流量为零, 而孔内却存在"抽水"和"注水"的井孔; 多层井(multilayerwell) 不一定要求"多层"的条件, 在均质单一含水层中的井孔可以具有多层井的基本特征; 混合井孔的水位并不"混合", 混合观测孔中存在符合机理的水头分布和流量分布规律等观点.普遍而言, 三维流中的观测孔不能用通常所说的线汇/线源刻画, 也不能用近几年提出的孔内为层流(线性流) 的假定来研究该问题, 然而可用笔者于1993年提出的"渗流-管流耦合模型"和"等效渗透系数"方便、有效地模拟.就说明性算例而言, Hantush Бочевер方程只能近似用于孔径大于0.2m且径距大于10~20m的条件.Abstract: The paper questions the traditional method for the calculation of the water level in the conventional observation wells for three-dimensional groundwater flow, and proposes that the Hantush -Бочевер equation, a pure mathematical method, is lacking in the basis of physics. The analysis of the formation mechanism of water level in the observation wells shows that the flux is zero only at the mouth of three-dimensional groundwater observation well, but the interior of the well contains actually "discharge" and "recharge" from the aquifers; that the water level in one homogeneous aquifer may have the same essential features as the multilayer well has; that the water level in the multilayer well is not mixed, and that the distribution of the head and flow in the multilayer well conforms to its own mechanism. Generally speaking, the observation well in the three-dimensional groundwater flow may not be described with the line sinks and line sources, nor can it be assumed that the flow in the observation well is determined by the laminar flow (linear flow) that has been proposed in recent years. However, the "seepage-pipe coupling model" and "equivalent hydraulic conductivity" (EHC) presented by the author in 1993 can be used to simulate very conveniently and effectively the three-dimensional groundwater flow. For the further illustration, the Hantush-Бочевер equation can only be approximately used in the conditions where the diameter of observation well is more than 0.2 m and that the distance between the well and the observation well is more than 10-20 m.
-
0. 引言
观测孔是用以监测地下水水头动态的主要手段, 是数值模型反演求取水文地质参数拟合的主要对象, 因此观测孔在水文地质学中具有极重要作用.
观测孔分两类[1]: 一是测压计式观测孔, 理论上它只有一个点与含水层相连通, 反映该点的水头值; 另一类是常规观测孔, 它有一定的长度与含水层相连通.通常所谓观测孔是指后者.美国著名学者Hantush[2]、前苏联著名学者Бочевер和Веригин[3]在对承压含水层非完整井流的研究中, 都认为观测孔中的水头降深s (r, l′, d′, t) 反映该孔滤水管中各点降深sp (r, z, t) 的平均值, 即
(1) 式中, l′、d′是观测孔滤管顶点、底点的标高.本文称之为Hantush-Бочевер方程.
对潜水三维流的研究, 美国著名学者Neuman[4]也认为, 观测孔中的水头降深可视为滤管内各点降深的平均值, 即
(2) 式中, r是观测孔至抽水井的距离; z1、z2是观测孔滤管顶点、底点的标高.上述诸学者的见解与确定观测孔中水位的方法在水文地质界一直沿用至今, 近半个世纪来未见到异议.然而, 上述观测孔中水位的确定方法, 缺乏对形成机理的基本分析.
简单分析之, 三维流中的垂直观测孔, 由于观测孔滤管中的水头不等, 会发生垂向流动, 而水的流动则会导致水头的再分布, 并引起井孔周围地下水的运动与水头再分布.因此, 与测压计式观测孔不同, 一般观测孔实质上是一部分滤管进水(抽水) 而另一部分滤管出水(注水), 即观测孔并不只简单地反映含水层地下水的水头, 而是兼有抽水与注水作用的井孔, 只是孔口的流量为零, 井管内的“抽水”量与“注水”量绝对值相等而已(如果忽略微小的井筒储存量的变化) (图 1).观测孔非“观测”孔也!
显然, 上述Hantush、Бочевер和Веригин以及Neuman提出的“积分平均水位”都未涉及水流的机理, 是缺乏物理基础的纯数学方法.
如果一个观测孔中部有粘性土夹层, 此即通常所谓的混合观测孔(多层观测孔), 因此混合观测孔也归属于三维流中的观测孔.而混合观测孔上下层水头差达几十cm甚至数m者, 是常见现象.
地下水的水位(水头) 靠观测孔监测; 观测孔的水位(水头) 又是地下水抽水试验和开采动态用以确定含水系统水文地质参数的最基本、最重要的数据, 也是地下水资源评价、地下水开采动态预测及管理的最重要的要素.因此, 对地下水三维流中观测孔水位的形成机理及确定方法的研究, 不管是在地下水资源评价还是地下水开采引发地质环境问题(地面沉降、海水入侵、地面塌陷、生态环境退化……), 以及地下水作为边坡稳定性中最活跃的因素作用的评价、治理等中都具有极其重要的应用前景.
1. 水流机理与模拟方法
如图 1所示, 在地下水三维流场中, 穿入一根不抽水的滤管——观测孔, 原来流场中滤管位置处的水头值不相等(如图 1所示下部水头值高于上部水头值), 滤管穿入后, 由于滤管的水流阻力远远小于原来的孔隙介质, 因此滤管中的原始水头差导致比原来孔隙介质条件下大得多的垂向流速和井筒垂向流量.如此, 依水流连续性原理, 井管下部必须从含水层中进水——抽水; 井管上部必须向含水层中出水——注水.滤管中水流流速(流量) 的增量(与设置观测孔之前的孔隙介质比较) 和滤管与孔隙介质间水量的交换(抽水、注水) 必导致观测孔及其周围含水介质水头的再分布.总的趋势是滤管内的水头差将会缩小, 但水头值不会变成相等, 否则, 滤管内水流将会停滞, 违背了部分进水部分出水的水流连续性原理.
由此可见, 即使是均质单含水层中的井孔, 只要存在垂直水力坡度分量, 就使其具有多层井(multilayer well) (混合井) 的基本性质.多层井不一定要求“多层”的条件, 均质单含水层也可存在“多层井”!
有了上述关于三维流观测孔中水位形成的基本认识后, 我们意识到, 若将观测孔视为边界, 用线汇(线源) 的理论方法处理观测孔是难以奏效的, 因为我们无法事先给出滤管段的流量或水头分配.
为了模拟三维地下水流观测孔中的水位及观测孔的存在对周围水头分布的影响, 本文采用陈崇希等[5]提出的“渗流-管流耦合模型”来模拟该问题, “渗流”刻画地下水的运动, “管流”刻画观测孔中的水流.笔者提出的方法是将管流部分视为渗透系数很大的圆柱形“透镜体”, 并对管流引入形式上服从达西定律的“等效渗透系数”, 从而将渗流-管流模型视为含有圆柱形“透镜体”的新的“渗流”系统.
对于这个新“渗流”系统, 关键问题是如何确定井管的“等效渗透系数”.
从流体力学知识可知, 当圆柱管中的水流呈层流状态时, 其水头损失可依Darcy-Weisbach方程计算[6], 即
(3) 式中: ΔH是水头损失(m); f是摩擦系数; l是管长(m); d是管内直径(m); u是管内平均流速(m/d); g是重力加速度(m/d2).
当管流为层流时,
(4) 式中: Re是雷诺数.而
(5) (6) 式中: ν是流体运动粘度; μ是流体动力粘度; ρ是流体密度.
将(4) 至(6) 代入(3) 式, 得
(7) 式中: J是水力坡度; γ是流体重度.
将上式写成渗透流速的形式, 对于管流, 其空隙率n=1, 则渗透流速V=nu=u, 故有
(8) 将此式与达西定律V=KJ对比, 写为
(9) 可得圆管“渗透系数”的表达式
(10) 其中, Kl定义为层流状态下管流的等效渗透系数, 方程(10) 已分别由陈崇希和林敏[7]以及Bear[8]获得.
当管流呈紊流状态时, 方程(3) 可改写为
(11) 如果定义紊流状态下管流的等效渗透系数Kn为
(12) (13) 该式也具有达西定律的形式.
本来, 紊流的运动规律有别于达西定律, 而且不同流态具有不同的形式.引入等效渗透系数后, 将5个流态分区(1个层流区和4个紊流区) 的运动规律统一为达西定律形式, 且与地下水渗流定律一致.即
(14) 其中:
(15) 需注意的是, 紊流条件下的等效渗透系数Kn与流速V及摩擦系数f有关, 它是随雷诺数而变化的量.另外, 雷诺数Re的确定与流速V有关, 而流速V又依赖于摩擦系数f, 摩擦系数f的确定反过来又取决于雷诺数Re, 因此必须采用迭代法来确定三者, 进而确定等效渗透系数.有了上述建立的等效渗透系数的概念和方程, 就可以利用常规的地下水流控制方程来刻画渗流-管流问题, 只要用等效渗透系数代替原控制方程中的渗透系数.
另外要注意: 从测压计原理出发, 混合观测孔中的水位反映的是当时流动条件下(包括井管中的流动) 观测孔(滤管) 揭露含水层最上点处的水头值.它是在与孔外(多孔介质) 流量交换条件下使得水头再分布而形成的[5].
2. 说明性算例
1993年笔者对井孔-含水系统提出通用的“渗流-管流耦合模型”之后, 已得到砂槽物理模拟[9]和现场试验[5]等的验证, 并成功地用于地下水混合井流[10, 11]、岩溶管道-裂隙-孔隙三重介质地下水流[12, 13]和水平井流[14]等问题.本文进一步将其用来研究三维地下水流中常规观测孔中水位形成机理及其某些规律, 并用下述理想模型加以说明.
本文提出的理想模型是: 均质、各向同性、等厚(厚度100 m) 的扇形(扇角为10°) 承压含水层, 有一个刚好打在扇角处穿入含水层顶面的点状非完整井, 由此抽水形成地下水三维流, 在此抽水井附近(扇形含水层的一侧面上) 打一完整观测孔.
上述理想的井孔-含水系统可用下面的数学模型刻画, 其中
(16) (17) 依上述,
(18) (19) 式中: Ke是等效渗透系数(m/d); s是水头降深(m); Qw是位于原点处抽水井的流量(m3/d); Vw是含水层中抽水井滤管所控制的体积(m3); Ss是单位储水系数(1/m); B2是零通量边界(所有的边界).其他符号的意义同上.
本文对该数学模型垂向上采用隐式有限差分法、平面上采用Galerkin有限元法使其离散化.当观测孔中的水流是紊流时, 这是一非线性问题.为了减少运算量而又能基本说明所研究的问题, 应该适当减少其结点数, 为此扇形含水层的扇角取为10°, 且只在扇边线上设置结点.为了说明基本数值解(观测孔的直径为零) 的精度, 将其与解析模型作对比, 为此全部边界取为不透水的, 且取模型的径距R足够大, 即R=868 146.5 m, 以保证扇外缘的边界对内部流场不产生明显的影响(模拟结果表明这一点).
关于结点的设置, 在抽水井的井心处设置结点, 并设为坐标原点.平面上采用不等距的结点设置, 近抽水井处结点间距小, 向外逐渐增大.平面上结点数计151个.对于本问题感兴趣的小于径距100 m的地段, X轴线上设置了26个结点.垂向上每5 m设置一个结点, 计21个结点.因此, 整个模型为151×21=3 171个结点.另外, 观测孔置于结点处.
设含水层的渗透系数为100 m/d, 单位储水系数为0.000 01 m-1, 抽水流量为1 000 m3/d, 观测孔的内径do和至抽水井的径距r取不同值, 以分析其完整观测孔的水位和其他效应.初始时间步长为0.001 d, 并以1.05的因子变步长地增加.运行步长数127个, 计10.188 99 d.当观测孔的直径为零时, 该问题为线性问题, 存在解析解.利用无限空间上一点汇作用下的水头降深方程, 通过无限次反映使含水层的上、下隔水边界同时获得满足而得解析解.无限空间上一点汇作用下的水头降深方程为
(20) 式中, s是水头降深; Q是点状抽水井的流量; K是渗透系数; ρ是观测点至点状抽水井的距离; a=K/Ss, 是含水层压力传导系数, 即扩散系数; t是抽水延续时间.
为了对比数值解与解析解, 设置6个测压计式观测点(表 1), 对比情况见图 2.由图 2可见, 拟合情况良好.
表 1 用于有限元解与解析解对比的观测点位置Table Supplementary Table Position of observation spot used to compare defined element result with analytical result对于具有一定直径的完整观测孔, 本文采用“渗流-管流耦合模型”模拟其水流运动规律.上文在理论上已作了必要的分析, 下面对其作具体应用.
本文要回答的主要问题是, 承压含水层三维地下水流中的完整观测孔中的水位(降深s) 是否由Hantush-Бочевер方程确定?对此可以根据“渗流-管流耦合模型”, 采用不同内径、不同径距的完整观测孔模拟试验得以证明.图 3a—3i分别表示径距r=2.583, 5.192, 11.839, 21.186, 30.948, 45.005, 65.247, 94.396和113.475 m的不同观测孔内径do (分别为0.00, 0.02, 0.05, 0.10和0.20 m, 其中do=0.00者即为测压计式观测点), M=100 m, K=100 m3/d, μs=0.000 01 m-1, t=10.019 6 d条件下的模拟结果和对比.由图 3可以看出:
(1) 观测孔中不同深度z处的水头降深s并不相等, 其顶、底层水头降深差Δs随观测孔径距的增大而减小.例如, 当do=0.02 m时, r=5.192 m和r=45.005 m的Δs分别为6.812 m和0.253 m, 有其明显的变化.其原因是这个水头降深差Δs的产生起因于三维流水力坡度的垂直分量, 而随着观测孔径距的增大, 其水力坡度的垂直分量会减小.此点说明, 不能简单地用等水头线刻画观测孔(和抽水井).
(2) 观测孔中顶、底层的水头降深差Δs随孔径的增大而变小, 例如, 当r=2.583 m时, do=0.02 m和do=0.20 m的Δs分别为13.758 m和0.03 m.这是因为存在do > 0的常规观测孔时, 其井孔滤管的等效渗透系数大于含水层的渗透系数而加大了垂向流量(可称为常规观测孔的“通道效应”), 为适应新情况的水均衡, 水头分布必定会自行调整.就本例而言, 观测孔的下部必须降低水头值(加大降深) 以吸收周围更多的来水量, 提供“加大了的垂向流量”; 观测孔的顶层必须提高水头值(减小降深) 以排泄“加大了的垂向流量”到周围的含水介质去, 结果导致顶、底模拟层水头差的减小.显然, 观测孔的孔径愈大, 孔滤管的等效渗透系数也愈大, 为满足水均衡, 观测孔顶、底模拟层水头的调整也愈大, 此即导致观测孔顶、底模拟层水头降深差Δs随孔径的增大而变小.
(3) 观测孔中的水位与按Hantush-Бочевер方程确定的值一般不相等.以do=0.02 m的观测孔为例, 当r=2.583 m时, 观测孔中的水位降深为17.199 m, 而按Hantush-Бочевер方程的计算值却仅有5.236 m, 两者相差达11.963 m之多.当r=45.005 m时, 观测孔中的水位降深为3.455 m, 而按Hantush-Бочевер方程的计算值为3.316 m, 两者相差为0.14 m, 仍然是个不可忽视的差别.只有当r=113.475 m时, 观测孔中的水位降深为2.792 m, 而按Hantush-Бочевер方程的计算值为2.788 m, 两者相差为0.004 m.这个差别, 一般情况下可以忽略不计.然而我们注意到, 这个径距已大于承压含水层的厚度100 m, 地下水流的垂向水头差非常之小, 通常可视为二维流了. (5) 沿观测孔井管中垂直向上流量的分布: 根据模拟获得观测孔内两结点的水头差(降深差)、井管等效渗透系数Ke (方程18或19) 和孔径, 可以计算两结点间井管段(中点) 的流量, 为此可获得流量沿井管的分布.图 4a—4d表示不同内径、不同径距的完整观测孔内抽水末期流量分布的数值模拟试验结果.图中, 若流量沿井管向上增大, 表示此结点起“抽水”作用; 若流量沿井管向上减小, 表示此结点起“注水”作用.下、上相邻两段(中点) 的流量差(Q上-Q下) 即为两段(中点) 间结点的“抽水”量, 负值为“注水”量.
(4) 表 2表示不同内径、不同径距的完整观测孔内水位降深数值模拟试验结果及与Hantush-Бочевер方程的对比.如果以0.01 m作为可以接受的水位误差的话, 那么在上述水文地质条件下, 只有当孔径do≥0.20 m且径距r大约大于10~20 m条件下, 观测孔水位可用Hantush-Бочевер方程近似确定.由于观测孔的孔径通常为0.02~0.05 m, 仅这一点而言, Hantush-Бочевер方程在大部分情况下不适用.
(21) 表 2 完整观测孔内水位降深与Hantush-Бочевер降深对比Table Supplementary Table Contrast of decline of complete observation well and decline of Hantush-Бочевер对于同一径距、同一标高处不同内径的观测孔, 内径大的流量大于内径小的, 这由井管等效渗透系数和孔径(过水断面积) 所决定.对于同一内径的观测孔, 同一标高处径距大的流量小于径距小的, 这是由垂直水力坡度随径距增大而减小所控制.
(6) 沿观测孔井管水流雷诺数的分布: 我们有兴趣地看一下不同径距、不同内径观测孔内水流雷诺数(方程5) 的分布, 图 5a—5d表示了这样一组曲线.对比相同的径距、相同的内径的雷诺数分布和流量分布, 其曲线形态是相似的(图 4a, 5a), 这不难理解; 但是上述“对于同一径距、同一标高处不同内径的观测孔, 内径大的流量大于内径小的”的规律, 对于雷诺数分布就不适用了, 因为若将Q与u的关系代入方程(5), 则雷诺数表示为
(22) 可见, 雷诺数不仅与流量呈正比, 还与孔径呈反比.因此对于同一径距、同一标高处, 内径大者的雷诺数可以大于、也可以小于内径小者的(图 5b).值得注意的另一现象是, 对于上述图 5的各种条件下, 没有一种条件能使观测孔中水流全段保持层流(线性流) 状态, 即使像图 5d中r=65.25 m和do=0.02 m、雷诺数偏小的条件, 全孔20段中也只有下部6段和上部4段为层流状态(Re≤2 300), 其他段已为非线性流.由此可见, 若采用Therrien和Sudicky[15]孔内为层流(线性流) 的假定来研究该问题是不妥的.尽管这种条件下直径2 cm的观测孔中的流量仅有3.192 m3/d, 已经出现非线性流了.
(7) 观测孔的“通道效应” (三维流中在垂直水力坡度作用下沿井管产生的水流效应) 不仅涉及观测孔自身的水头变化, 理论上要影响整个流场, 在本例中, 最上的第一模拟层地下水水头降落漏斗的反应最为明显(图 6a-6d)——可称为降落漏斗的畸变.其中, 随着观测孔径距的减小和孔径的加大, 笫一层的降落漏斗的提高值(与无观测孔相比较) 愈加明显(图 6a).
由于受比例尺的限制, 图 6在降落漏斗的改变量上表现得不是很明显, 特别是当观测孔的径距rg较大时.为此表 3中分别表示不同径距rg存在观测孔条件下, 不同内径do和径距r处降落漏斗的上升值ΔH (等于存在观测孔条件下的水头-无观测孔条件下的水头).由表 3可见, 观测孔的通道效应在改变降落漏斗方面是相当可观的, 位于径距2.583 m的观测孔, 对径距r=0.000 m处的水头影响, 观测孔直径do=0.10 m者可使水头上升11.374 m; do小到0.02 m者也可使水头上升1.289 m; 即使对于径距r大到54.206 m处的水头影响, 一般也可达到数cm.对于三维流特征已明显衰减的观测孔径距rg达到45.005 m的情况, 径距r分别为0.000和54.206 m处的降落漏斗的上升值ΔH尚有0.052和0.077 m, 而观测孔处(r=45.005 m) 则达0.142 m.如果抽水试验水位观测值的误差如此之大, 那是不允许的.降落漏斗的畸变现象应当引起人们的注意.
表 3 径距不同时观测孔的通道效应ΔH (do, r)Table Supplementary Table Channel effect of observation well in different radius distances3. 结论
(1) 1961年美国学者Hantush[2]和前苏联学者Бочевер和Веригин[3]分别认为承压含水层三维流观测孔中的水头反映该孔滤水管中各点水头的平均值, 这种提法未涉及水流的机理, 是缺乏物理基础的纯数学方法. (2) 通过观测孔中水流的机理分析, 三维流中的观测孔实质上是一部分滤管进水(抽水) 而另一部分滤管出水(注水), 它是兼有“抽水”与“注水”作用的井孔, 只是孔口的流量为零.观测孔非“观测”孔也! (3) 就普遍而言, 三维流中的观测孔不能用通常所说的线汇、线源刻画, 也不能用近几年提出的孔内为层流(线性流) 的假定来研究该问题, 然而可用笔者1993年所提出的“渗流-管流耦合模型”方便、有效地模拟. (4) 多层井(multilayer well) 不一定要求“多层”的条件, 均质单一含水层也可存在多层井的特征. (5) 采用“渗流-管流耦合模型”刻画观测孔——形成新的含强透水圆柱形透镜体的非均质含水层模拟表明, 混合井并不“混合”, 观测孔中存在符合机理的水头分布和流量分布规律. (6) 观测孔的“通道效应”不仅涉及观测孔自身的水头变化, 理论上还影响整个流场.就本文算例而言, 若在45.0 m范围内存在直径do≥0.05 m的完整观测孔, 对径距r≤54.0 m范围内的水头都要产生cm~m级(部分达十多m) 的畸变; 对于直径do=0.02 m的完整观测孔, 则产生mm~m级的畸变.
在上述说明性实例的水文地质条件下, 当孔径do≥0.20 m且径距r大约大于10~20 m条件下, 观测孔水位才可以用Hantush-Бочевер方程近似确定.
致谢: 博士研究生黎明为论文绘制图件和表格, 表示感谢! -
表 1 用于有限元解与解析解对比的观测点位置
Table 1. Position of observation spot used to compare defined element result with analytical result
表 2 完整观测孔内水位降深与Hantush-Бочевер降深对比
Table 2. Contrast of decline of complete observation well and decline of Hantush-Бочевер
表 3 径距不同时观测孔的通道效应ΔH (do, r)
Table 3. Channel effect of observation well in different radius distances
-
[1] 陈崇希. 地下水不稳定井流计算方法[M]. 北京: 地质出版社, 1983.222.CHEN C X. Calculation methodsof groundwater unsteadily flowing to well[M]. Beijing: Geological Publishing House, 1983.222. [2] Hantush M S. Aquifer test on partially penetrating wells [J]. Proc Am Soc Civil Engrs, 1961, 87 (5): 171 - 195. [3] Бочевер Ф М, Веригин Н Н. Методическое пособие по расчетам зксплуатащонных запасов подземных вод для водоснабжения[M]. Москва: Игостроииздат, 1961. [4] Neuman S P. Theory offlow in unconfined aquifers considering delayed response of the water table[J]. Water Resour Res, 1972, 18(4): 1031-1045. [5] 陈崇希, 林敏, 叶善士, 等. 地下水混合井流理论及其应用[M]. 武汉: 中国地质大学出版社, 1998.43.CHEN C X, LIN M, YE S S, et al. Groundwater flow model of mixed well and its application[M]. Wuhan: China University of Geosciences Press, 1998.43. [6] Streeter V L. Handbook of fluid dynamics[M]. New York: McGraw-Hill Book Co, 1961. [7] 陈崇希, 林敏. 地下水动力学[M]. 武汉: 中国地质大学出版社, 1999.CHEN C X, LIN M. Groundwater hydraulis[M]. Wuhan: China University of Geosciences Press, 1999. [8] Bear J. Dynamics of fluids in porous media[M]. New York: Elsevier, 1972. [9] Chen C X, Wan J W, Zhan H B. Theoretical and experimental studies of coupled Darcian-pipe flow to a horizontal well[J]. J of Hydro, 2003, (5)(in press). [10] 陈崇希. 地下水不稳定混合抽水的模型及模拟方法[J]. 地球学报, 1996, 17(增刊): 222-228. https://cpfd.cnki.com.cn/Article/CPFDTOTAL-DQXB199612001033.htmCHEN C X. The model and simulated method for the mixpumping in the unstable groundwater[J]. Earth Journal, 1996, 17(Suppl): 222-228. https://cpfd.cnki.com.cn/Article/CPFDTOTAL-DQXB199612001033.htm [11] Chen C X, Jiao J J. Numerical simulation of pumping test in multilayer wells with non-Darcian flow in the wellbore [J]. Ground Water, 1999, 37(3): 465-474. doi: 10.1111/j.1745-6584.1999.tb01126.x [12] 陈崇希. 岩溶管道-裂隙-空隙三重空隙介质地下水流模型及模拟方法研究[J]. 地球科学——中国地质大学学报, 1995, 20(4): 361-366.CHEN C X. Groundwater flow model and simulation method in triple medium of karst tube-fissure-pore[J]. Earth Science-China University of Geosciences, 1995, 20(4): 361-366. [13] Cheng J M, Chen C X. Preliminary numerical study of karst groundwater flow in Beishan area[A]. Proceedings of 2nd international conference of future groundwater at risk[C]. Changchun: [s. n. ], 1998.65-66. [14] 陈崇希, 万军伟. 地下水水平井流的模型及数值模拟方法——考虑井管内不同流态[J]. 地球科学——中国地质大学学报, 2002, 27(2): 135-140. https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX200202004.htmCHEN C X, WAN J W. A new model of groundwater flowing to horizontal well and the numerical simulation approach[J]. Earth Science-Journal of China University of Geosciences, 2002, 27(2): 135-140. https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX200202004.htm [15] Therrien R, Sudicky E A. Well bore boundary conditions for variably saturated flow modeling[J]. Adv Water Resour, 2001, 24(2): 195-201. -