潜水稳定井流的剖面二维数值模拟方法

王旭升; 谢永桦; 陈崇希

doi:10.19509/j.cnki.dzkq.tb20230024

潜水稳定井流的剖面二维数值模拟方法

doi: 10.19509/j.cnki.dzkq.tb20230024

1a.
中国地质大学(北京)水利部地下水保护重点实验室(筹), 北京 100083
1b.
中国地质大学(北京)水资源与环境学院, 北京 100083
2.
中国地质大学(武汉)环境学院, 武汉 430078

基金项目:

国家自然科学基金项目 41972263

国家自然科学基金项目 41772249

详细信息

作者简介:
王旭升(1974—), 男, 教授, 博士生导师, 主要从事地下水动力学和水文模型研究工作。E-mail: wxsh@cugb.edu.cn

中图分类号: P641
计量
- 文章访问数: 989
- PDF下载量: 95
- 被引次数: 3
出版历程
- 收稿日期: 2023-01-13
- 录用日期: 2023-04-19
- 修回日期: 2023-04-14

Sectional 2D numerical modelling method for steady state well-flow in an unconfined aquifer

1a.
Key Laboratory of Groundwater Conservation of Ministry of Water Resources(in preparation), China University of Geosciences(Beijing), Beijing 100083, China
1b.
School of Water Resources & Environment, China University of Geosciences(Beijing), Beijing 100083, China
2.
School of Environment Studies, China University of Geosciences(Wuhan), Wuhan 430078, China

摘要

摘要:
经典Dupuit井流模型与考虑入渗的改进Dupuit井流模型, 都受到Dupuit假定的影响, 可能存在系统误差。建立反映三维流或轴对称二维流特性的潜水井流数值模型, 是检验Dupuit模型可靠性的重要手段。提出一种模拟潜水稳定井流的剖面二维数值模拟方法, 通过参数转换把柱坐标系的渗流方程变换为等效的直角坐标系方程, 利用MODFLOW方体网格有限差分模型实现剖面二维流场模拟。数值模型以抽水井的井中水位为已知条件, 渗出面的排水量按照Darcy定律差分公式计算, 潜水面则通过MODFLOW对干-湿单元的处理加以识别, 抽水流量经水均衡计算得到。通过采用精细化网格建立典型案例模型, 获得模拟精度很高的结果, 使反算抽水井流量的相对误差不超过0.2%。以此检验Dupuit井流模型, 发现解析公式得到的水位线总体与数值模拟结果一致, 仅在抽水井附近由于没有考虑水跃而偏低, 且误差受到含水层渗透系数各向异性的影响。在有入渗的情况下, 分水岭附近的渗流违反Dupuit假定。然而, 改进的Dupuit井流公式计算的分水岭水位相对误差低于1%。这一数值模拟方法简单实用, 但也受到MODFLOW本身局限性的约束。
- Dupuit假定 /
- 抽水井 /
- 入渗补给 /
- 渗出面 /
- 潜水面 /
- 各向异性 /
- 有限差分法
Abstract: Objective
Both the classical Dupuit model and the modified Dupuit model including infiltration are influenced by Dupuit's assumption and have potential systematic errors. Building a numerical model of well flow in an unconfined aquifer by characterizing the three-dimensional or axisymmetric two-dimensional (2D) flow is an essential approach to verify the performance of Dupuit-type models.
Methods
In this study, a 2D numerical model is proposed for the steady state well-flow in an unconfined aquifer, in which the control equation of seepage in the cylindrical coordinates is transformed equivalently to the Cartesian coordinates through parameter transformation, and the sectional 2D modelling is implemented via the MODFLOW finite-difference grid of cubic blocks. In the numeric model, the water level in the pumping well is a given condition, the flux across the seepage face is estimated by difference equation according to Darcy's law, the phreatic surface is identified by the treatment of dry and wet cells in MODFLOW, and the pumping rate is determined from the water debug calculation.
Results
Fine grids are constructed innumerical models of typical cases to obtain high-precision results, in which the relative error of the backwards estimated pumping rate is no more than 0.2%. This numerical model is used to check the Dupuit-type well-flow models. As indicated, the groundwater level estimated from the analytical formulas generally agrees well with the numerical modelling results, except that near the well, where the analytical solution underestimates the groundwater level due to ignoring the waterjump and the errors depend on the anisotropic permeability of the aquifer. When infiltration exists, the flow in the vicinity of watershed does not follow Dupuit's assumption. However, the estimated groundwater level on the watershed by modified Dupuit well-flow equation has a low level of relative error, which is less than 1%.
Conclusion
This numerical method is simple and practical however it is also influenced by limitations in MODFLOW.
- Dupuit's assumption /
- pumping well /
- infiltration recharge /
- seepage face /
- phreetic surface /
- anisotropic /
- finite difference method
注释:

(所有作者声明不存在利益冲突)

HTML全文

地下水动力学中的井流模型具有轴对称特征，使用柱坐标系进行描述。经典Dupuit井流模型^[1]引入Dupuit假定，把地下水流简化为单纯的径向一维流问题，便于获得控制方程的解析解。Dupuit井流模型忽略了降水入渗的影响。为此，有研究者提出了改进模型^[2-4]，将入渗补给考虑在径向一维流问题的控制方程中。这些模型能够提供抽水井周围潜水面形态的理论公式。

轴对称状态的潜水井流实际属于柱坐标系下的二维渗流过程，采用Dupuit假设是将其简化为径向一维流问题的手段。在含水层的局部地段(如地下分水岭处)，地下水的垂向流速可能显著大于水平流速，不满足Dupuit假定。那么，如何检验Dupuit假定对潜水稳定井流问题的有效性？这是地下水动力学中一个受到关注的问题。2020年陈崇希^[4]在提出考虑入渗的改进Dupuit模型时，也指出所建立的方程包含Dupuit假定，可能在地下分水岭和抽水井附近带来误差，但只做了简单定性讨论，未定量分析Dupuit假定对解析结果的影响。

检验Dupuit假定的影响，需要求解柱坐标系的二维潜水井流问题。不过，由于包含潜水面这样的自由面边界，获得解析解的难度很大。潜水井流问题还涉及另一个需要考察的实际因素，就是渗出面。地下水流向抽水井时，井壁水位高于井中水位(即水跃)，介于两者之间的一段井壁就是渗出面。这是地表水等势体附近流线和等水头线弯曲状态(不满足Dupuit假设)所产生的必然结果，其存在性与流场是否稳定无关^{[2, 5-7]}。前苏联学者已经证明^[5-6]：即使存在水跃，Dupuit模型的流量公式仍然是可靠的(无入渗补给条件)，尽管在抽水井附近实际地下水位高于Dupuit公式计算的结果。对于潜水井流的水跃现象，Boulton^[7]利用Laplace方程的有限差分解法进行了系统的研究。张有龄^[8]在阐述潜水井流问题和水跃现象时，详细引用了杨式德所做的有限差分模型，其方法与Boulton所用方法类似。其他学者也用有限差分法^[9]或有限单元法^[10]对此做了数值模拟研究。由于非线性带来的困难，目前尚未研究得到水跃的严格解析解。

前人所做的潜水井流数值模拟没有考虑降水入渗补给，而且需要采用柱坐标系进行网格离散，不容易在现有的通用数值模拟工具中实现。如果能够在直角坐标系建立与柱坐标系等价的剖面二维渗流数值模型，将为模型网格离散以及处理非均质含水层、入渗补给、渗出面等水文地质要素带来便利。目前国际上流行的三维地下水流有限差分模拟程序MODFLOW^[11]采用的是直角坐标系，可在一些专业软件进行快速建模。笔者拟利用控制方程的等效原理，在直角坐标系下采用等效方法模拟潜水稳定井流剖面二维流场，与Dupuit模型或改进Dupuit模型解析解进行对比，对有关问题进行讨论。

1. 直角坐标系模拟轴对称渗流方法

1.1 两种坐标系的渗流方程

在直角坐标系中(图 1-a)，二维地下水稳定渗方程^[12]的一般形式为：

图 1 直角坐标系(a)和柱坐标系(b)的潜水面示意图(z为相对高度; x为横向距离; r为径向距离；后同)

Figure 1. Schematics of the phreetic surface in the cartesian (a) and cylindrical (b) coordinate systems

下载: 全尺寸图片幻灯片

$\frac{\partial}{\partial x}\left(K_x \frac{\partial H}{\partial x}\right)+\frac{\partial}{\partial z}\left(K_z \frac{\partial H}{\partial z}\right)=0$

(1)

式中：H为水头；x，z分别为水平与垂直坐标，而K_x和K_z分别为对应坐标方向的渗透系数。该方程成立的条件是渗透张量的主轴与坐标轴一致，且流场内部没有源汇项。潜水面作为二维流场的上边界，其稳态约束条件可表示为^[12]：

$K_x \frac{\partial H}{\partial x} \frac{\partial z_{\mathrm{wt}}}{\partial x}-K_z \frac{\partial H}{\partial z}+w=0$

(2)

式中：z_wt为潜水面的高度(随x坐标变化)；w为直角坐标系入渗补给强度。

在柱坐标系中(图 1-b)，轴对称渗流方程^[13]的一般形式为：

$\frac{1}{r} \frac{\partial}{\partial r}\left(K_{\mathrm{h}} r \frac{\partial H}{\partial r}\right)+\frac{\partial}{\partial z}\left(K_{\mathrm{v}} \frac{\partial H}{\partial z}\right)=0$

(3)

式中：r为径向距离；K_h和K_v分别为水平和垂直方向的渗透系数。这种柱坐标系下的潜水面稳态边界条件表示为：

$K_{\mathrm{h}} r \frac{\partial H}{\partial r} \frac{\partial z_{\mathrm{wt}}}{\partial r}-K_{\mathrm{v}} r \frac{\partial H}{\partial z}+\varepsilon r=0$

(4)

式中：ε为柱坐标系下的入渗补给强度。

1.2 等效变换方法

已有一些研究者提出了在直角坐标系等效模拟轴对称地下水流的方法。Samani等^[14]将直角坐标系的水平距离变换为x=ln(r)建立MODFLOW模型，然后按照这种对数变换，保持水平渗透系数不变，而对垂向渗透系数实施变换，产生等效结果。Langevin^[15]提出另外一种在MODFLOW中转换参数的方法，他发现MODFLOW用对数加权平均法计算的块体单元界面等效渗透系数正好与Dupuit公式的效果一致。Louwyck等^[16]建立x方向均匀划分的MODFLOW模型，然后通过参数变换获得等效模拟结果，将水平渗透系数变换为：

$K_x=\frac{2 \pi K_{\mathrm{h}}}{\ln \left(r_{j+1 / 2} / r_{j-1 / 2}\right)}$

(5)

式中：j为单元体在x方向的离散坐标；r_j+1/2，r_j-1/2分别为单元外侧和内侧的径向距离。这种方法实际上也是取x与lnr成正比。

本研究从控制方程等价的角度建立变换关系，即保持直角坐标系的水平距离x与柱坐标系的径向距离r相同，通过参数变换，使式(1)与式(3)等价，而式(2)与式(4)等价。满足这种等价关系的变换式为：

$K_x=\frac{2 \pi x K_{\mathrm{h}}}{\Delta y}, K_z=\frac{2 \pi x K_{\mathrm{v}}}{\Delta y}, w=\frac{2 \pi x \varepsilon}{\Delta y}, x=r$

(6)

式中：乘积因子2πx为径向距离r处完整过水断面的周长；Δy为模型在x，z正交方向的宽度。在式(6)中变量与Δy相除，可以保持参数的量纲不变。从物理意义上来讲，Δy是剖面模型所代表的过水断面单位宽度。MODFLOW模型默认单元都是三维方体，因此必然会有y方向的宽度Δy，将其设定为一个任意单位宽度(例如1 m)代入式(6)即可。Langevin^[15]提出的方法本质也是如此，但他没有从控制方程等价的角度进行解释，也没有处理渗出面。

1.3 渗出面的处理

渗出面是地下水流的特殊边界，在忽略流速水头的情况下，其水头与渗出点的位置高度相等。对于有入渗补给的潜水稳定井流模型，渗出面既可能出现在抽水井的井壁，也可能出现在外周含水层与大气的交界面上。假设井壁和含水层外周均沿铅直方向延伸，分别位于r=r_w(井半径)和r=R(含水层外周半径)处，则渗出面的边界条件可表示为：

$H(x, z)=z, x=r_{\mathrm{w}}, h_{\mathrm{w}}<z<\left(h_{\mathrm{w}}+L_{\mathrm{w}}\right)$

(7)

$H(x, z)=z, x=R, h_{\mathrm{R}}<z<\left(h_{\mathrm{R}}+L_{\mathrm{R}}\right)$

(8)

式中：h_w为抽水井中的水位；h_R为含水层外周水位，而L_w和L_R分别为井壁与含水层外周渗出面产生的水跃值。

式(7，8)表明渗出面属于已知水头的边界。然而，在实际的模拟过程中，预先并不知道L_w和L_R的数值，因此难以在模型中直接设置某个高度段为已知水头边界。为了能让模型识别出渗出面，可以改用以下描述方式：

$V_x=-K_x \frac{\partial H}{\partial x}, H\left(r_{\mathrm{w}}+\delta_x, z>h_{\mathrm{w}}\right)>z$

或

$\begin{aligned} & H\left(R-\delta_x, z>h_{\mathrm{R}}\right)>z \\ & V_x=0, H\left(r_{\mathrm{w}}+\delta x, z>h_{\mathrm{w}}\right) \leqslant z \end{aligned}$

(9)

或

$H\left(R-\delta_x, z>h_{\mathrm{R}}\right) \leqslant z$

(10)

式中：δ_x>0，表示偏移一小段距离。采用偏移小段距离的水头大小来判断是否处于渗出面附近，符合数值模型中边界单元的特性。以井壁附近的单元为例，如图 2-a所示，第k层和k+1层的相邻单元(宽度为Δx)计算水头分别为H(k)和H(k+1)，它们均与x=r_w+0.5Δx处的潜水面高度z_w近似相等，但是只有第k层的水头满足H(k)>z(k)而处于渗出面附近。这里有δ_x=0.5Δx。

图 2-a中，第k层单元与渗出面接触，并且排泄的侧向流量Q_x按照Darcy定律计算，即将式(9)的差分格式与单元界面的面积相乘，计算为：

图 2 含水层网格单元与井孔的关系示意图

a.井壁渗出面相邻单元；b.井孔饱水带相邻单元；r_w为抽水井半径；Δx为单元长度；Δz为单元厚度；Δy为单元宽度；Q_x为侧面流量；z_w为井中水位；z(k)为单元中点高度；H(k)为单元格点水头

Figure 2. Schematics of the relationship between grid blocks of the aquifer and the well

下载: 全尺寸图片幻灯片

$Q_x=-K_x \frac{2\left(H_k-z_k\right)}{\Delta x} \Delta z \Delta y, H_{\mathrm{k}}>z_{\mathrm{k}}$

(11)

Q_x为负值表示流向与x轴相反。图 2-a中第k+1层单元格点高于潜水面，符合式(10)所述条件，即Q_x=0。在含水层外周，可以证明存在类似关系。因此，可以将渗出面的边界条件替换为第三类边界，即：

$\begin{aligned} & Q_{\mathrm{d}}=C_{\mathrm{d}}\left(H_{\mathrm{k}}-z_{\mathrm{k}}\right), H_{\mathrm{k}}>z_{\mathrm{k}} ; \\ & Q_{\mathrm{d}}=0, H_{\mathrm{k}} \leqslant z_{\mathrm{k}} \end{aligned}$

(12)

式中：Q_d为模型单元向外部排泄的流量；C_d称为模型单元的排水系数。根据式(11)，并把式(6)代入，可得：

$C_{\mathrm{d}}=\frac{2 K_x \Delta z \Delta y}{\Delta x}=\frac{4 \pi x K_{\mathrm{h}} \Delta z}{\Delta x}, x=r_{\mathrm{w}} \text { 或 } x=R$

(13)

显然，排水系数的取值与模型网格剖分有关，并非自然界含水层的固定参数。

实际上，井中水位以下的井壁(水头值为h_w)，也可以用类似的方法处理。如图 2-b所示，设第k层有一个模型单元与井孔单元相邻，其宽度为Δx，2个单元交界面上的水头固定为h_w，则侧向流量的差分公式为：

$Q_x=-K_x \frac{2\left(H_{\mathrm{k}}-h_{\mathrm{w}}\right)}{\Delta x} \Delta z \Delta y$

(14)

式(14)与处理渗出面的式(11)是类似的。因此也可使用类似式(12)的公式计算单元的排泄流量：

$Q_{\mathrm{d}}=C_{\mathrm{d}}\left(H_{\mathrm{k}}-h_{\mathrm{w}}\right)$

(15)

其中排水系数C_d的计算方法与式(13)相同。

2. 模型设计与MODFLOW用法

2.1 含水层模拟情景

本研究模拟均质含水层中的潜水井流问题，以便与Dupuit模型和改进的Dupuit模型^[4]解析解进行对比。模型如图 3所示，有一半径为R=2 000 m的圆岛状均质含水层，外周边界水位为h_R=50 m，圆岛中心处有一个半径为r_w=0.1 m的抽水井，以恒定流量Q_w抽水，最终形成井中水位h_w。含水层从隔水底板到地面的总厚度为70 m，地下水接受大气降水入渗补给，强度为ε，最大可达1 mm/d。地下水流场达到稳定，形成非线性的潜水面形态，并可能在抽水井的井壁发育长度为L_w的渗出面(即水跃值为L_w)，在含水层外周边界发育长度为L_R的渗出面。

图 3 潜水井流模型示意图

x, r为距离；R为圆岛半径；L_R为外侧边界水跃；H(r, z)为水头分布函数；h_w为井中水位；L_w为井壁水跃

Figure 3. Diagrammatic map of the model for the wellflow in an unconfined aquifer

下载: 全尺寸图片幻灯片

按照含水层渗透系数与入渗补给强度的变化，设置不同情景进行模拟。情景编号和参数见表 1。

表 1 模拟情景编号和参数

Table 1. Code and parameters of modelling scenarios

模拟情景	水平方向渗透系数K_h/(m·d^-1)	各向异性比值K_v/K_h	入渗补给强度ε/(mm·d^-1)	井中水位h_w/m
A	10	1	0.0	35
B	10	1	1.0	35
C	10	0.1, 10.0	1.0	35
D	1	0.1	1.0	20

下载: 导出CSV

| 显示表格

2.2 离散网格

模拟工作是在Aquaveo公司发布软件Groundwater Modeling System (GMS)中完成的，构建三维方体有限差分网格(y方向只取一个单元宽度)，调用MODFLOW模块处理水文地质要素。

模型网格的初始剖分方式为：在垂向上分为70层，每层厚度1 m；水平方向也进行密集划分，总体每个网格单元的宽度为2 m，抽水井本身占一排宽度为0.1 m的单元，其附近含水层单元加密，宽度缩小到0.01 m。模型中井孔所在的这一列单元设置为无效单元，即并不直接参与有限差分计算，而是以边界的方式作用到模型中(具体处理见2.3节)。模型最右侧一排单元与含水层外周边界相邻。在模拟过程中，根据渗出面和潜水面分水岭的结果进行模型网格的微调，提高精度。每个含水层单元的水平渗透系数K_x和垂向渗透系数K_v，根据式(6)计算，其中x为单元形心与抽水井中心线的距离。

2.3 抽水井和外周边界的处理

本次模拟并不像解析模型那样把抽水井处理为已知流量的边界，因为抽水井本身也被离散为很多个单元格(类似混合井)，事先并不能准确掌握每个单元分配的流量。为此，本研究将抽水井的水位处理为已知条件，即给定h_w的数值，然后在模拟形成稳定流之后反求出抽水井的流量。在井壁处，z坐标高于h_w的点可能位于渗出面，其水头值与z坐标一致。与井壁相邻的一排模型单元，采用MODFLOW中的Drain模块^[11]处理地下水向抽水井的排泄。高于井中水位的边界单元，排泄量用式(12)计算，以z_k为Drain模块的排水高程。低于井中水位的边界单元，排泄量用式(15)计算，以h_w为Drain模块的排水高程。Drain模块的排水系数统一采用式(13)计算，取x=r_w。

调用MODFLOW运行稳定流模拟之后，再根据模拟结果确定抽水井的流量Q_w。它实际上由2个部分构成的：一部分是从渗出面排出的地下水量，一部分是从井水位下方井壁排出的地下水量。因此，抽水井流量的统计可以表示为:

$Q_{\mathrm{w}}=\sum\limits_{k=1}^N C_{\mathrm{d}}\left(H_{\mathrm{k}}-h_{\mathrm{w}}\right)+\sum\limits_{k=N+1}^M C_{\mathrm{d}}\left(H_{\mathrm{k}}-z_{\mathrm{k}}\right)$

(16)

式中：N为井中水位以下的最高井壁单元层号；M为对应渗出面的最高井壁单元层号。从模拟结果计算Q_w的方法，是将井壁相邻的饱和带单元组合设置为一个均衡区，进行编码，调用MODFLOW的Budget模块计算均衡区流进和流出的水量。没有计算误差的情况下，两者的绝对值应相等，并且与Q_w一致。如果潜水面存在分水岭，也可以将分水岭半径范围以内的全部入渗补给水量计为Q_w。这2种计算方法的结果不一定相等，只要控制相对误差不超过1%，就可以认为数值模拟结果是可靠的。

在模型的外周边界处，渗出面的情况与井壁渗出面是类似的，也用Drain模块处理。对高于水位h_R的边界单元，排泄量用式(12)计算，以z_k为Drain模块的排水高程，取x=R计算排水系数C_d。低于h_R的边界单元则将水头固定为h_R。

2.4 潜水面的处理

潜水面属于自由水面，各点的水头值与其高度相等，这可以用来判断潜水面的位置。与此同时，潜水面还需要满足式(2)，作为渗流场的边界条件。MODFLOW并不直接求解式(2)，而是以入渗补给量作为源、建立潜水面所在网格单元的水均衡方程，与饱和带的方程组一起求解。从数值计算的角度来讲，这与求解式(2)是等价的。如果计算水头高于单元的底部，则该单元为“湿单元”，是有效的。否则，就判断为“干单元”，相当于无效单元，不再参加下一次计算，而是将潜水面移入下部的相邻单元^{[11, 17]}。MODFLOW的Recharge模块可以自动将入渗补给加入到最高处的湿单元上。模型入渗强度w根据式(6)计算，其中x为单元形心的r值。

模拟过程中，潜水面的形成实际上表现为网格单元的逐渐疏干过程。先将所有不属于边界的网格单元水头设置为模型顶部高程(70 m)。然后启动渗流场计算和“干”、“湿”单元判断，反复迭代，逐渐将水头高于潜水面的单元疏干，而余下的网格单元构成饱和带的渗流空间。迭代过程中可能会出现疏干过快或振荡现象，导致计算不稳定，因此需要对迭代算法进行控制，使疏干过程平稳进行。

2.5 求解算法和收敛精度

MODFLOW形成有限差分方程组后，可采用多种方法进行求解。本次模拟采用预共轭梯度法(PCG2模块^[18])，这是一种迭代计算效果较高的算法。为保证数值模拟精度，PCG2模块中水头计算的收敛标准设置为0.001 m，水均衡计算的收敛标准设置为0.01 m³/d。由于潜水面的形成是一个将包气带单元逐步疏干的过程，不宜太快，将PCG2模块中迭代计算的松弛因子(relaxation parameter)从默认的1.0调整为0.5。

3. 典型模拟结果分析

3.1 对比Dupuit模型解析解

表 1中的情景A是无入渗、均质各向同性含水层条件，剖面二维渗流模拟结果可以与Dupuit公式进行对比。图 4给出了抽水井附近的模型网格与模拟结果。显然，潜水面的模拟水位高于解析解，在井壁处模拟得到渗出面的顶部高为40.5 m，比井中水位35.0 m高出了5.5 m，这就是水跃值。圆岛外周的水跃值小于1 m。

图 4 情景A模拟水位与经典Dupuit模型解析解的局部对比

Figure 4. Local comparison between the modelled groundwater level of the scenario-A and the analytical solution of the classical Dupuit model

下载: 全尺寸图片幻灯片

根据模拟结果，对抽水井流量进行了计算，得到Q_w=4 051.8 m³/d，其中从渗出面排泄的地下水量为364.2 m³/d，占9.0%。模型的总体水均衡误差(流入、流出之差)绝对值小于0.01 m³/d，说明有很高的差分计算精度。解析方法计算流量的Dupuit公式^[1]为：

$Q_{\mathrm{w}}=\frac{\pi K_{\mathrm{h}}\left(h_{\mathrm{R}}^2-h_{\mathrm{w}}^2\right)}{\ln \left(R / r_{\mathrm{w}}\right)}$

(17)

把含水层参数和边界水位代入式(17)，得到Q_w=4 044.6 m³/d。与之相比，数值模型得到的流量偏高0.2%。Dupuit模型的流量公式对三维流也是适用的，而且不受水跃的影响，这一点已经得到了严格的证明^[5]。因此，上述0.2%的相对误差代表了数值模型在流量计算上的误差水平。

3.2 对比改进的Dupuit模型解析解

表 1中的情景B是有入渗、均质各向同性含水层条件，剖面二维渗流模拟结果可以与改进的Dupuit公式^[4]进行对比(图 5)。宏观上看，模拟的潜水面曲线与解析解很接近，在分水岭处几乎重合(图 5-a)。在抽水井附近，解析解则明显低于模拟结果(图 5-b)，渗出面产生的水跃值L_w≈7.1 m。

对于情景B，抽水井流量与井中水位的解析关系为^[4]：

$Q_{\mathrm{w}}=\frac{\pi K_{\mathrm{h}}\left(h_{\mathrm{R}}^2-h_{\mathrm{w}}^2\right)}{\ln \left(R / r_{\mathrm{w}}\right)}+\frac{\varepsilon \pi R^2}{2 \ln \left(R / r_{\mathrm{w}}\right)}$

(18)

把参数和边界条件代入式(18)，得到抽水井流量的解析解为Q_w=4 679.0 m³/d，高于无入渗的情景。数值模型反算的结果为Q_w=4 687.9 m³/d，其中渗出面流量510.7 m³/d，占10.9%。模型的水均衡误差仍然小于0.01 m³/d，数值模拟得到的抽水流量比解析解偏大0.2%，相对误差很小。这种一致性也反映在分水岭上，解析解确定的分水岭在x=1 220.4 m处，数值模拟的分水岭在x=1 221 m处，相对误差仅有0.05%，而且数值模拟的分水岭水位是50.52 m，与解析解一致。分水岭处地下水主要呈现垂向流动，明显违反Dupuit假设。然而，从本例对比结果来看，引入Dupuit假设并没有影响解析解计算分水岭水位的准确性。

图 5 情景B模拟水位与改进Dupuit模型解析解的对比

a.整体范围；b.抽水井附近

Figure 5. Comparison between the modelled groundwater level of scenario-B and the analytical solution of the modified Dupuit model

下载: 全尺寸图片幻灯片

3.3 有入渗的各向异性含水层情景

表 1中的情景C和D都是各向异性含水层，可以考察含水层各向异性特征对潜水井流的影响。改进Dupuit模型^[4]并没有包含对各向异性的处理。

情景C与情景B的基本条件是相同的，唯一区别在于情景C的垂向渗透系数与水平渗透系数不同。根据模拟结果(图 6)，各向异性对水位的影响主要表现在抽水井附近(图 6-a)，而距离大于100 m时的水位线对各向异性比(K_v/K_h)不敏感(图 6-b)。分水岭位置和水位的模拟结果与情景B几乎一致，说明各向异性对抽水井流量也几乎没有影响。从图 6-a可以看出，K_v/K_h越小(即垂向渗透系数越低)则模拟水位线越高于改进Dupuit模型的解析解。当K_v/K_h=0.1时，水跃值为L_w≈9.5 m。当K_v/K_h=10时，水跃值只有L_w≈4.5 m。这说明：垂向渗透系数越大，解析公式得到的水位越逼近实际水位。

情景D的水平渗透系数为其他情景的1/10，而且井中水位比其他情景低了10 m。模拟结果见图 7。与情景B和C相比，分水岭的位置(在r=643 m)更加靠近抽水井，但水位更高(62.0 m)，相应的抽水流量也降低到Q_w=1 287.0 m³/d，这主要是水平渗透性降低导致水流阻力增加的结果。这种情况下的水跃值也很大，达到L_w≈33.5 m，使得在抽水井附近模拟水位远高于改进Dupuit公式的解析解。当r>117 m时，解析公式计算水位的相对误差低于1%。在分水岭处，解析公式计算水位为61.8 m，相对误差约0.4%。

图 6 情景C的模拟水位与改进Dupuit模型解析解的对比

a.抽水井附近；b.全部范围

Figure 6. Comparison between the modelled groundwater level of scenario-C and the analytical solution of the modified Dupuit model

下载: 全尺寸图片幻灯片

图 7 情景D的模拟结果

Figure 7. Numerical modelling results of scenario-D

下载: 全尺寸图片幻灯片

4. 讨论

4.1 流速的换算

如果要根据上述剖面二维渗流模拟结果绘制流线，需要注意2个问题。第1个问题是等效模型计算的流速与轴对称模型是否具有相同的方向？如果方向相同，那么流线的形态是一样的。第2个问题是两者计算的流速大小是否相同？如果有差异，则意味着不能直接计算流线上质点的运移时间。下面对此进行讨论。

轴对称坐标系与直角坐标系的Darcy流速分别为：

$V_{\mathrm{r}}=-K_{\mathrm{h}} \frac{\partial H}{\partial r}, V_{\mathrm{v}}=-K_{\mathrm{v}} \frac{\partial H}{\partial z}$

(19)

$V_{\mathrm{x}}=-K_x \frac{\partial H}{\partial x}, V_z=-K_z \frac{\partial H}{\partial z}$

(20)

式中：V_r和V_v分别为轴对称坐标系的径向流速与垂向流速；V_x和V_z分别为直角坐标系的横向流速与垂向流速。把1.2节所提供的参数等效变换方法代入流速计算式，保持(∂H/∂r)=(∂H/∂x)，可建立2种坐标系Darcy流速的转换关系：

$V_{\mathrm{r}}=\frac{\Delta y}{2 \pi x} V_{\mathrm{x}}, V_{\mathrm{v}}=\frac{\Delta y}{2 \pi x} V_{\mathrm{z}}$

(21)

式(21)表明，2种坐标系Darcy流速是不同的，但V_r/V_v与V_x/V_z的数值是相同的，即流动方向不变。因此，2种坐标系的流线形态是相同的，而计算质点运移时间需要用式(21)换算流速。GMS软件提供了MODPATH模块与MODFLOW对接，可追踪粒子绘制流线，采用的是Pollock算法^[19]。该算法计算粒子对流运移速度用到介质的孔隙度，即Darcy流速与孔隙度的比值。因此，在调用MODPATH模块时，把实际孔隙度与式(21)中的比例系数(2πx/Δy)相乘作为等效孔隙度输入，即可得到符合轴对称坐标系的运移时间。

图 7展示了MODPATH模块所绘制情景D中的流线，且流线上相邻箭头之间表示一段相同的运移时间。潜水面上粒子的投放间距随径向距离的增大而增加，以保证任意2条相邻流线之间获得的补给流量为固定值。在分水岭下方，可以追踪出一根几乎垂直的流线，即为分水线(三维空间为分水面)，其内侧区域的流线较为稀疏(反映抽水井流量占补给流量较小)，流线上的箭头也较少(反映地下水流速较快)。

4.2 为何Dupuit假设未导致显著的分水岭水位误差

在有入渗的情况下，只要抽水流量小于圆岛范围的入渗补给总量，就会出现分水岭。在分水岭处，地下水存在垂直向下的流速，而水平流速为零，显然违反Dupuit假定。如果按照直觉进行推测，Dupuit假定应该在分水岭处造成较大的误差。然而，在本研究的有入渗情景B、C和D中，可用改进的Dupuit井流公式计算分水岭的水位，与数值模拟结果对比，其相对误差均低于1%，说明并没有受到Dupuit假定的显著影响。下面对于这种“意外”或“反常”现象做一个理论上的探讨。

在剖面二维的潜水井流模型中，底部为隔水边界，因此垂向流速为零，而顶部潜水面处的边界条件由式(4)描述。如果忽略二阶小项，则流速的上、下边界条件可近似为：

$\left.V_{\mathrm{v}}\right|_{z=h}=-\left.K_{\mathrm{v}} \frac{\partial H}{\partial z}\right|_{z=h} \approx-\varepsilon,\left.V_{\mathrm{v}}\right|_{z=0}=0$

(22)

式中：V_v为垂向Darcy流速，以向上为正；h为潜水面相对底板的高度，即以底板为零高程时式(4)中的z_wt。尽管垂向流速随z的变化是非线性的，我们可以根据边界条件用线性化公式来近似评估水力梯度，即：

$\frac{\partial H}{\partial z} \approx \frac{\varepsilon}{K_{\mathrm{v}}} \frac{z}{h}$

(23)

因此水头随z变化的近似式为：

$H(z) \approx H_{z=0}+\frac{\varepsilon}{2 K_{\mathrm{v}} h} z^2$

或

$H(z) \approx h-\frac{\varepsilon}{2 K_{\mathrm{v}} h}\left(h^2-z^2\right)$

(24)

设分水岭处的水位为h_d，若采用Dupuit假定忽略垂向水力梯度，则底部水头的假设值为H_z=0=h_d，而根据式(24)其实际值(取z=h_d)近似为：

$H_{z=0} \approx h_{\mathrm{d}}-\frac{\varepsilon h_{\mathrm{d}}}{2 K_{\mathrm{v}}}$

(25)

因此Dupuit假定计算含水层底部水头的相对误差为:

$e=\frac{\left|h_{\mathrm{d}}-H_{z=0}\right|}{h_{\mathrm{d}}} \approx \frac{\varepsilon}{2 K_{\mathrm{v}}}$

(26)

该误差水平取决于入渗强度与垂向渗透系数的比值。在案例情景B，C和D中，ε=0.001 m/d，K_v=0.1~100 m/d，因此理论上近似估计的相对误差不超过0.5%。

如果要分析分水岭处水位的计算误差，则式(26)应改为：

$e=\frac{\left|h_{\mathrm{d}}-h_{\mathrm{D}}\right|}{h_{\mathrm{d}}}$

(27)

式中：h_D是采用改进Dupuit模型^[4]计算的分水岭水位，即：

$h_{\mathrm{D}}^2=h_{\mathrm{R}}^2+\frac{\varepsilon}{K_{\mathrm{h}}}\left(\frac{R^2-r_{\mathrm{d}}^2}{2}-r_{\mathrm{d}}^2 \ln \frac{R}{r_{\mathrm{d}}}\right)$

(28)

其中r_d为分水岭的横向坐标。实际水位h_d目前还没有严格的解析解，所以与数值模拟结果对比是检验Dupuit假定的主要方法。为了掌握更普遍的情况，我们做理论上的近似分析。首先，从水均衡角度可以得到：

$\begin{array}{*{20}{l}} {Q(r){\text{ }} = 2\pi \int\limits_0^{\text{h}} {{V_{\text{h}}}rdz = \pi \varepsilon \left( {{r^2} - r_{\text{d}}^2} \right)} ,} \\ {\;\;\;\;\;\;\;{r_{\text{d}}} \leqslant r \leqslant R} \end{array}$

(29)

式中：Q(r)为穿过距离r处圆环过水断面的流量(向外为正)；V_h为水平方向的Darcy流速。把Darcy定律形式V_h=-K_h∂H/∂r代入式(29)，得到：

$\begin{array}{*{20}{l}} {\int\limits_0^{h} {\frac{{\partial H}}{{\partial r}}{\rm{d}}z = \frac{{\rm{d}}}{{{\rm{d}}r}}\left( {\int_0^h H {\rm{d}}z} \right) - h\frac{{{\rm{d}}h}}{{{\rm{d}}r}} = \frac{\varepsilon }{{2{K_{\rm{h}}}}}\left( {\frac{{r_{\rm{d}}^2}}{r} - r} \right),} }\\ \ \ \ \ \ \ {{r_{\rm{d}}} \le r \le R} \end{array}$

(30)

在式(30)等号两边对r积分，利用外侧边界条件但忽略渗出面，可得以下近似解：

$\begin{array} \mathop \smallint \limits_0^{\text{h}} H \mathrm{~d} z-\frac{1}{2} h^2 \approx \\ \ \ \ \ \frac{1}{2} h_{\mathrm{R}}^2+\frac{\varepsilon}{2 K_{\mathrm{h}}}\left(\frac{R^2-r^2}{2}-r_{\mathrm{d}}^2 \ln \frac{R}{r}\right), r_{\mathrm{d}} \leqslant r \leqslant R \end{array}$

(31)

把式(24)代入式(31), 有：

$\begin{array}{*{20}{l}} {\left( {\frac{1}{2} - \frac{\varepsilon }{{3{K_{\text{v}}}}}} \right){h^2} \approx } \\ \ \ \ \ {\quad \frac{1}{2}h_{\text{R}}^2 + \frac{\varepsilon }{{2{K_{\text{h}}}}}\left( {\frac{{{R^2} - {r^2}}}{2} - r_{\text{d}}^2\ln \frac{R}{r}} \right),{r_{\text{d}}} \leqslant r \leqslant R} \end{array}$

(32)

把r=r_d代入式(32)得到分水岭处实际水位的近似解，即:

$\begin{array}{*{20}{l}} {\left( {\frac{1}{2} - \frac{\varepsilon }{{3{K_{\text{v}}}}}} \right)h_{\text{d}}^2 \approx } \\ \ \ \ \ {\quad \frac{1}{2}h_{\text{R}}^2 + \frac{\varepsilon }{{2{K_{\text{h}}}}}\left( {\frac{{{R^2} - {r^2}}}{2} - r_{\text{d}}^2\ln \frac{R}{{{r_{\text{d}}}}}} \right) = \frac{1}{2}h_{\text{D}}^2} \end{array}$

(33)

由式(33)得到h_D/h_d，再代入式(27)，有：

$e \approx 1-\sqrt{1-\frac{2 \varepsilon}{3 K_{\mathrm{v}}}}, \varepsilon<\frac{3}{2} K_{\mathrm{v}}$

(34)

式(34)从理论上评估了采用Dupuit假定计算分水岭水位的相对误差，也取决于入渗强度与垂向渗透系数的比值。在案例情景B、C和D中，ε/K_v的数值小于0.01，误差应不超过0.34%，这与情景D数值模拟检验的误差不超过0.4%是基本一致的。因此，只要ε/K_v很小，即使采用Dupuit假设计算分水岭水位，也不会导致显著误差。

4.3 方法特点与局限

数值模拟是研究地下水运动普遍采用的方法，模型设计会直接影响模拟过程的仿真性与可靠性^[20]，模拟方法随含水层介质和水流过程的特点而变化^[21-22]。用于检验解析理论的地下水流数值模型不一定要做得很复杂，但需要抛弃一些理论模型的简化假设，考虑更符合实际的因素。本研究所提出的潜水稳定井流模拟方法能够在直角坐标系下完成建模，得到与轴对称模型等价的模拟结果，同时又能够处理渗出面与潜水面这样的自由面边界，比经典Dupuit模型或改进的Dupuit模型更加具有普遍适用性。式(6)这样的变换式较为简单，比前人提出的对数变换公式^{[14, 16]}更实用，而且能够保证不同坐标系下的流速具有式(21)所示的正比关系。该方法便于利用MODFLOW方体网格建立潜水稳定井流模型。地下水与井孔、外周渗出面的水量交换用Darcy定律差分公式等效处理为Drain边界，如式(12, 15)所示，能够与MODFLOW中的Drain模块兼容。对二维渗流场自由面边界的处理，则利用了MODFLOW处理“干”、“湿”单元的技术。采用精细化网格能够提高数值模拟的精度。根据情景A和情景B模拟结果与解析解的对比，抽水井流量和分水岭的模拟结果具有很高的精度，相对误差不超过0.2%，说明数值模型是可靠的。基于这一点，我们判断水跃值的模拟结果也是可靠的，但目前尚未研究得到水跃的精确解，因此不能对其计算误差做出准确的评估。

本研究提出的方法也能够处理非均质含水层，只要把非均质条件下随坐标变化的K_h和K_v代入式(6)即可。不过，这样的变换式对强烈非均质含水层可能还是有缺陷的。数值模型需要计算相邻单元界面的等效渗透系数或导水系数，当K_h和K_v随坐标强烈变化时，式(6)可能会增加差分公式的误差。目前，渗出面和自由面的算法处理对MODFLOW有依赖性，也会把MODFLOW本身的缺陷带入到模拟过程中。例如Drain模块需要输入的排水系数C_d并非含水层的真实物理参数，而是随着模型网格剖分情况而变化，使用不慎可能导致系统误差或错误的理解。用干-湿单元区分潜水面边界也会产生截断误差，并影响迭代计算的收敛性，如果垂向网格剖分太粗，难以得到光滑的潜水面形态。采用节点法构建模型网格，并把潜水面处理为可移动的节点，效果可能会比MODFLOW更好。此外，如果抽水井附近存在高速非Darcy流或井管内的流动对渗流场有较大影响，本研究方法也是不适用的。这些问题的处理需要引入非Darcy渗流公式^[23-24]以及渗流-管流耦合模型^[25]，有待进一步研究。

5. 结论

(1) 本研究提出了一种模拟潜水稳定井流的剖面二维数值模拟方法，根据柱坐标系渗流方程与直角坐标系方程的等效关系，建立了参数转换关系式，从而能够利用方体网格有限差分模型(如MODFLOW)实现剖面二维流场的数值模拟。该模拟方法中的数值模型考虑了二维的渗流场，也考虑了抽水井和外侧边界可能出现的渗出面(水跃)，并且把潜水面作为自由面边界，渗出面的排水量采用Darcy定律的差分公式计算，潜水面则通过MODFLOW对干-湿单元的处理加以识别。抽水井的井中水位作为已知条件配置在模型中，而抽水流量是在模拟结果基础上通过水均衡计算得到的。流量计算的误差水平可以用于检验数值模型的精度。本研究采用精细化网格建立的数值模型，典型算例中抽水井流量的相对误差不超过0.2%，模型精度很高。

(2) 数值模型考虑了二维流、渗出面和自由面边界，不受Dupuit假定的限制，因此可以用来检验包含Dupuit假定的井流理论模型。通过把无入渗情景的经典Dupuit模型和有入渗的改进Dupuit模型解析解与案例中的模拟水位进行对比，发现解析公式得到的水位线总体上与数值模拟结果一致，仅在抽水井附近由于没有考虑水跃而偏低。其误差受到含水层渗透系数各向异性的影响，垂向渗透性越低，解析解计算的水位越偏低。在有入渗的情况下，分水岭附近的渗流违反Dupuit假定。尽管如此，改进的Dupuit井流公式计算的分水岭水位相对误差低于1%，仍然具有很高的可靠性。本研究的数值模拟方法可用于模拟非均质各向异性含水层潜水稳定井流的二维流场，操作方式明确，但也受到MODFLOW本身局限性的约束。

(所有作者声明不存在利益冲突)

图 1 直角坐标系(a)和柱坐标系(b)的潜水面示意图(z为相对高度; x为横向距离; r为径向距离；后同)

Figure 1. Schematics of the phreetic surface in the cartesian (a) and cylindrical (b) coordinate systems

下载: 全尺寸图片幻灯片

图 2 含水层网格单元与井孔的关系示意图

Figure 2. Schematics of the relationship between grid blocks of the aquifer and the well

下载: 全尺寸图片幻灯片

图 3 潜水井流模型示意图

x, r为距离；R为圆岛半径；L_R为外侧边界水跃；H(r, z)为水头分布函数；h_w为井中水位；L_w为井壁水跃

Figure 3. Diagrammatic map of the model for the wellflow in an unconfined aquifer

下载: 全尺寸图片幻灯片

图 4 情景A模拟水位与经典Dupuit模型解析解的局部对比

Figure 4. Local comparison between the modelled groundwater level of the scenario-A and the analytical solution of the classical Dupuit model

下载: 全尺寸图片幻灯片

图 5 情景B模拟水位与改进Dupuit模型解析解的对比

a.整体范围；b.抽水井附近

Figure 5. Comparison between the modelled groundwater level of scenario-B and the analytical solution of the modified Dupuit model

下载: 全尺寸图片幻灯片

图 6 情景C的模拟水位与改进Dupuit模型解析解的对比

a.抽水井附近；b.全部范围

Figure 6. Comparison between the modelled groundwater level of scenario-C and the analytical solution of the modified Dupuit model

下载: 全尺寸图片幻灯片

图 7 情景D的模拟结果

Figure 7. Numerical modelling results of scenario-D

下载: 全尺寸图片幻灯片

表 1 模拟情景编号和参数

Table 1. Code and parameters of modelling scenarios

模拟情景	水平方向渗透系数K_h/(m·d^-1)	各向异性比值K_v/K_h	入渗补给强度ε/(mm·d^-1)	井中水位h_w/m
A	10	1	0.0	35
B	10	1	1.0	35
C	10	0.1, 10.0	1.0	35
D	1	0.1	1.0	20

下载: 导出CSV

参考文献(25)

[1]	Dupuit A J E J. Etudes theoretiques et pratiques sur le mouvement des eaux[M]. Paris: Dunod, 1863.
[2]	Каменский Г И. Основы дина-мики подземных вод[M]. Москва: Госгеолиздат, 1943(in Russian).
[3]	Haitjema H M. Analytic element modeling of groundwater flow[M]. San Diego: Academic Press, Inc., 1995.
[4]	陈崇希. Dupuit模型的改进: 具入渗补给[J]. 水文地质工程地质, 2020, 47(5): 1-4. https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG202106001.htm Chen C X. Improvement of Dupuit model: With infiltration recharge[J]. Hydrogeology & Engineering Geology, 2020, 47(5): 1-4(in Chinese with English abstract). https://www.cnki.com.cn/Article/CJFDTOTAL-SWDG202106001.htm
[5]	Чарный И A. Строгое доказателъство формулм дюпюидля безнапорной филътрации с промеҗутком высачнвания[J]. Докл. АН СССР, 1951, 79(6): 937-948. https://www.cnki.com.cn/Article/CJFDTOTAL-NWYJ202105015.htm
[6]	陈崇希, 林敏. 地下水动力学[M]. 武汉: 中国地质大学出版社, 1999. Chen C X, Lin M. Groundwater hydraulics[M]. Wuhan: China University of Geosciences Press, 1999(in Chinese).
[7]	Boulton N S. The flow pattern near a gravity well in a uniform water-bearing medium[J]. Journal of the ICE, 1951, 36(10): 534-550.
[8]	张有龄. 河床地下水运动的供水理论分析[M]. 北京: 科学出版社, 1958. Zhang Y L. Theoretical analysis on water yield from groundwater flow beneath riverbed[M]. Beijing: Science Press, 1958(in Chinese).
[9]	Taylor G S, Luthin J N. Computer methods for transient analysis of water-table aquifers[J]. Water Resources Research, 1969, 5(1): 144-152. doi: 10.1029/WR005i001p00144
[10]	Neuman S P, Witherspoon P A. Finite element method of analyzing steady seepage with a free surface[J]. Water Resources Research, 1970, 6(3): 889-897. doi: 10.1029/WR006i003p00889
[11]	McDonald M G, Harbaugh A W. A modular three-dimensional finite-difference ground-water flow model[R]. Denver: Techniques of Water-Resources Investigations of the U.S. Geological Survey, Chapter A1, Book 6, 1988.
[12]	王旭升, 万力. 地下水运动方程[M]. 北京: 地质出版社, 2011. Wang X S, Wan L. Equations of groundwater movements[M]. Beijing: Geological Publishing House, 2011(in Chinese).
[13]	陈崇希, 林敏, 成建梅. 地下水动力学[M]. 北京: 地质出版社, 2011. Chen C X, Lin M, Cheng J M. Groundwater hydraulics[M]. Beijing: Geological Publishing House, 2011(in Chinese).
[14]	Samani N, Kompani-Zarea M, Barry D. MODFLOW equipped with a new method for the accurate simulation of axisymmetric flow[J]. Advances in Water Resources, 2004, 27(1): 31-45. doi: 10.1016/j.advwatres.2003.09.005
[15]	Langevin C D. Modeling axisymmetric flow and transport[J]. Groundwater, 2008, 46(4): 579-590. doi: 10.1111/j.1745-6584.2008.00445.x
[16]	Louwyck A, Vandenbohede A, Bakker M, et al. MODFLOW procedure to simulate axisymmetric flow in radially heterogeneous and layered aquifer systems[J]. Hydrogeology Journal, 2014, 22(5): 1217-1226. doi: 10.1007/s10040-014-1150-0
[17]	董佩, 王旭升. MODFLOW模拟自由面渗流的应用与讨论[J]. 工程勘察, 2009(7): 27-30. https://www.cnki.com.cn/Article/CJFDTOTAL-GCKC200907009.htm Dong P, Wang X S. Application and discussion of MODFLOW's simulation to the seepage of free surface[J]. Journal of Geotechnical Investigation & Surveying, 2009(7): 27-30(in Chinese with English abstract). https://www.cnki.com.cn/Article/CJFDTOTAL-GCKC200907009.htm
[18]	Hill M C. Preconditioned conjugate-gradient 2 (PCG2), a computer program for solving, ground-water flow equations[R]. Water-Resources Investigations Report 90-4048, Denver: U.S. Geological Survey, 1990.
[19]	Pollock D W. User's Guide for MODPATH/MODPATH-PLOT, Version 3: A Particle tracking post-processing package for MODFLOW, the U.S. Geological survey finite difference groundwater flow model[R]. Open-file report 94-464, Denver: U.S. Geological Survey, 1994.
[20]	陈崇希, 唐仲华, 胡立堂. 地下水流数值模拟理论方法及模型设计[M]. 北京: 地质出版社, 2014. Chen C X, Tang Z H, Hu L T. Theory, method and model design for numerical simulation of groundwater flow[M]. Beijing: Geological Publishing House, 2014 (in Chinese).
[21]	成建梅, 罗一鸣. 岩溶多重介质地下水模拟技术及应用进展[J]. 地质科技通报, 2022, 41(5): 220-229. doi: 10.19509/j.cnki.dzkq.2022.0220 Chen J M, Luo Y M. Overview of groundwater modeling technology and its application in karst areas with multiple-void media[J]. Bulletin of Geological Science and Technology, 2022, 41(5): 220-229(in Chinese with English abstract). doi: 10.19509/j.cnki.dzkq.2022.0220
[22]	郑小康, 杨志兵. 岩溶含水层饱和-非饱和流动与污染物运移数值模拟[J]. 地质科技通报, 2022, 41(5): 357-366. doi: 10.19509/j.cnki.dzkq.2022.0211 Zheng X K, Yang Z B. Numerical simulation of saturated-unsaturated groundwater flow and contaminant transport in a karst aquifer[J]. Bulletin of Geological Science and Technology, 2022, 41(5): 357-366(in Chinese with English abstract). doi: 10.19509/j.cnki.dzkq.2022.0211
[23]	Wen Z, Liu Z T, Jin M G, et al. Numerical modeling of Forchheimer flow to a pumping well in a confined aquifer using the strong-form mesh-free method[J]. Hydrogeology Journal, 2014, 22(5): 1207-1215. doi: 10.1007/s10040-014-1136-y
[24]	Wang Q R, Zhan H B, Tang Z H, et al. Forchheimer flow to a well-considering time-dependent critical radius[J]. Hydrology and Earth System Sciences, 2014, 18(6): 2437-2448. doi: 10.5194/hess-18-2437-2014
[25]	陈崇希, 林敏, 叶善士, 等. 地下水混合井流的理论及应用[M]. 武汉: 中国地质大学出版社, 1998. Chen C X, Lin M, Ye S S, et al. Theory of multi-layer mixed well flow and its application[M]. Wuhan: China University of Geosciences Press, 1998(in Chinese).