Numerical model of two-region flow by constant-head pumping in an unconfined aquifer
-
摘要:
为了揭示在非承压含水层中定水头抽水试验引起的达西−非达西两区流动机理，提出基于有限差分法的地下水定水头抽水井流数值模型。该模型根据抽水的流态特征将含水层分为2个区域：靠近抽水井的有限非达西渗流区域和远离抽水井的半无限达西渗流区域，其中非达西流区域流态的模拟基于Izbash方程实现。通过与COMSOL Multiphysics的有限元数值解进行比较，验证了所提出数值解的可靠性。最后，研究有限非达西流效应对水头和抽水井抽水速率的影响以及井内水头对抽水井抽水速率的影响。研究表明：在抽水试验中非达西区域的影响不可忽略，湍流会分别导致两区流中水头较纯非达西流和纯达西流的水头偏大和偏小，且随抽水时间的增加逐渐变大；通过减小抽水井井内水头或增大非达西系数可提高抽水速率，但该影响会随抽水时间的增加而逐渐减弱；断面流量随径向距离的增大而不断减小，断面流量与径向距离曲线下降速率不断减小，且在转换界面处会出现转折点。该模型为定量研究在非达西流和达西流耦合作用下抽水井附近的井流水头特征提供了一种简洁的方法，并为调查定水头抽水测试期间的抽水速率提供理论依据。
Abstract:In order to reveal the mechanism of two-region Darcian and non-Darcian flow induced by constant-head pumping tests in unconfined aquifers, a numerical model of well flow by constant-head pumping underground based on finite difference solution was proposed. In the model, it is assumed that the aquifer is divided into two regions according to the pumping flow characteristics: the finite non-Darcian flow region near the pumping well and the semi-infinite Darcian flow region far away from the pumping well. Specifically, the flow regime in the non-Darcian flow region was simulated by the Izbash’s equation, and the reliability of the proposed solution was verified by comparison with finite element numerical solutions by COMSOL Multiphysics. Finally, the influence of finite non-Darcian effect on hydraulic head and pumping rate, as well as the influence of hydraulic head in the pumping well on pumping rate, was especially studied. The results show that the influence of non-Darcian flow region in the pumping test cannot be ignored. The turbulent flow makes the hydraulic head of the two-region flow larger than that of the pure non-Darcian flow and smaller than that of pure Darcian flow. Besides, such difference of the hydraulic head is increased with the pumping time. The pumping rate can be increased by reducing the hydraulic head in the pumping well or increasing the non-Darcian coefficient, but the effect is gradually decreased as the pumping continues. The flow rate at the cross-section is decreased with the increase of radial distance, the gradient of the flow rate-radial distance curve is decreased with time, and a turning point will appear at the conversion interface between the two regions. The proposed model provides a simple method for quantitative studies on the characteristics of hydraulic head near the pumping well under the coupling of non-Darcian and Darcian effects, and provides a theoretical basis for determining the pumping rate during constant-head pumping test.
-
定水头抽水试验是一种特殊的抽水试验，其井中水头在抽水时需保持稳定[1]，常被应用于渗透系数较低的地层以及抽水量难以维持稳定的工况[2]。定水头抽水试验对于确定水力参数或解决垃圾填埋场引起的环境破坏问题具有重要意义[2-3]。自1950年以来，水文地质学家一直在研究由定水头抽水试验引起的抽水井流问题，C. E. Jacob等[4]第一个提出了用于研究定水头抽水试验期间的达西流流动特性的解析解。此后，水文学家们在定水头抽水井流的理论和物理模型方面进行了大量有价值的工作[2,5-8]。
一般情况下，抽水引起的地下水流动遵循达西定律，符合线性流动假设[9]。然而，许多学者发现，当水力梯度较大或较小时，流量与水力梯度会呈非线性变化关系。这种非线性抽水井流被定义为非达西流[10-11]，通常由Izbash方程或Forchheimer方程描述，并根据现场条件在2个方程之间进行选择[12-13]。在过去10年中，大量学者基于上述2种方法，通过原位试验研究了定水头抽水引起的非达西流流动特性问题[5,8,12,14-15]。例如，P. M. Quinn等[15]在加拿大安大略省圭尔夫市附近的裂缝性白云岩和砂岩中进行了定水头原位试验，并指出使用线性特征的达西流来拟合实验数据可能会产生很大的误差。Dan Hancheng等[6]基于定水头试验的降深数据，对无黏结级配骨料的参数进行评估，发现即使在低水力梯度下非达西流现象也很普遍。Liu Zhongyu等[8]使用改进的渗透率测试装置在饱和黏土含水层中进行定水头测试，发现抽水试验会导致水流具有明显的非线性特性。这些研究认为，忽略非达西效应会使抽水井流模拟和确定含水层参数产生较大的误差。然而，现有的定水头抽水模型大多是基于现场或室内试验数据拟合得出，适用地质范围较为狭窄，且不具备理论依据。到目前为止，关于定水头抽水引起的非达西流数值模型研究很少，提出可通用于常见含水层系统的相关数值模型是具备现实价值的。
由文献评述笔者发现，现有定水头抽水试验研究通常基于抽水含水层是承压层的假设上提出的。然而，世界上许多抽水井都处于非承压含水层中[16-17]，并且关于非承压含水层定水头抽水井流的研究一直较为缺乏。到目前为止，现有的非承压含水层抽水井流模型主要分为考虑和不考虑垂向流分量2类。在考虑垂直流分量模型中，重力释水不是瞬间释放的，这进一步导致水头上方的非饱和带会对下方饱和带进行补给[18]。基于含水层重力储水量远大于弹性储水量的事实，在抽水井内降深较小或抽水时间足够长的定水头试验中，垂直流的影响有限，可以合理使用不考虑垂向流分量的Dupuit假设，该假设认为重力储存为瞬时释放[19-21]。根据Dupuit的调查，当非承压含水层降深变化较缓时，可以忽略垂直流的影响，并假设所有补给均来自水平方向。因此，它可用于简化无约束流动问题，使其在解析上易于处理。在非承压含水层中，当水位降深相对初始水位较小时，可以优先考虑使用Dupuit假设来模拟流动，以满足实际工程的需要[2,13,20,22-24]。
从理论上说，由于抽水井附近的流速足够大，会对相应区域的流态造成影响。足够大的流速可以使雷诺数大于临界雷诺数，从而导致抽水井附近水流呈现非达西流态。而随着径向距离的增加，雷诺数随流速逐渐减小。因此，它可以进一步引起非达西流向达西流的转变，导致含水层中呈现出两区流特性[25-29]。迄今为止，在非承压含水层中通过定水头试验引发的这种两区流动机制仍不清楚，相关研究不足。
针对上述问题，笔者提出一种非承压含水层定水头抽水试验导致的非达西−达西两区渗流数值模型。当含水层储水系数较小时，通过Izbash方程描述比流量和水力梯度之间的关系，以进行数学建模。基于Dupuit假设，利用有限差分法推导出易于实际操作的简洁数值模型。所提出模型通过与COMSOL Multiphysics的数值解进行比较，验证渗流模型的可靠性，以期为两区流特性分析奠定基础。
1 两区流数学模型
图1为在非承压含水层中定水头抽水试验导致的两区渗流模型示意图，模型考虑了非达西和达西2个区域的水流流动，其数学模型的假设包括：(1) 该含水层是非承压、均质、各向同性的非承压含水层，初始水头为一定值；(2) 抽水井完全穿透非承压含水层，抽水井内的水头自抽水开始后保持不变；(3) 非达西区域(简称N区)为有限区域，主要发生在抽水井周围；达西区域(简称D区)为半无限区域，位于N区外围；(4) 抽水井半径与有效半径相同，为一个固定正数[23]。
![]()
图 1 考虑两区转换的非承压含水层定水头抽水模型rw—抽水井有效半径，m；hw—抽水井中水头，m；R—N区和D区转换界面到抽水井中心的距离，m；b—含水层厚度，mFigure 1. Schematic of the constant-head test in an unconfined aquifer considering the two-region transform基于质量守恒定律和Dupuit假设，根据Boussinesq方程将N区和D区中的抽水流量分别用下式描述。
$$ \frac{\partial {q}_{1}}{\partial r}+\frac{{q}_{1}}{r}=\frac{{S}_{1}}{{h}_{1}}\frac{\partial {h}_{1}}{\partial t} $$ (1) $$\frac{\partial {q}_{2}}{\partial r}+\frac{{q}_{2}}{r}=\frac{{S}_{2}}{{h}_{2}}\frac{\partial {h}_{2}}{\partial t} $$ (2) 式中：
${q}_{1}$ 和${q}_{2}$ 分别为N区和D区中的比流量，m/s；${h}_{1}$ 和${h}_{2}$ 分别为N区和D区中的水头，m；上述参数均为时间$ t $ 和水平距离$ r $ (指某一点与抽水井中心的水平距离)的函数；$ {S}_{1} $ 和$ {S}_{2} $ 分别为N区和D区的储水系数。假设N区和D区之间的转换界面位于
$ r=R $ 处，根据流体连续性假设得到边界处的边界条件如下：$$ \underset{r\to R}{\mathrm{lim}}{q}_{1}=\underset{r\to R}{\mathrm{lim}}{q}_{2} $$ (3) $$ {h}_{1}{|}_{r\to R}={h}_{2}{|}_{r\to R} $$ (4) 达西定律描述了水力坡度(
$ I $ )和比流量($ q $ )之间的线性关系。对于多孔介质流，研究表明其线性关系在以下2种情况下是不成立的。一种是在较低水力梯度下通过低渗介质的流动，另一种则是在较高水力梯度下通过粗粒高渗介质的流动。J. Bear[30]提到，基于平均直径的雷诺数($ Re $ )在$Re < 1 \sim 10$ 时，达西定律是适用的。在此范围内，渗流流动可以被归纳为层流，称为达西流，否则即可被称为非达西流。此外，$ Re $ 较高时($ Re $ >10)的流动为湍流。对于抽水井附近的流动，水力梯度通常很高，流速随之增大[31]。1.1 非达西区域
Izbash方程可用于描述非达西流的
$I $ 和$ q $ 之间的关系。使用该非线性函数的一个假设是，至少在某些流动阶段，所涉及的非达西系数$ n $ 和准渗透系数$k_{\rm{c}}$ 与空间和时间无关[32]。其公式表示为：$I=c{q}^{n}$ ,${c}$ 为常数，当参数$ 1 < n < 2 $ 时，该公式用于表征因显著惯性效应导致的非线性渗流特性；当$ n=1 $ 时，公式用于表示线性达西定律[33]；$ 0 < n < 1 $ 时，该公式可用于表征低渗介质中因显著的固液界面效应导致的非线性渗流。当
$ Re > 10 $ 时，可认为该区域处于N区，其抽水井流为非达西渗流[25]，N区的比流量可用Izbash方程表示。$$ {q}_{1}=\left({K}_{1}\frac{\partial {h}_{1}}{\partial r}\right)^{\tfrac{1}{{n}_{1}}} $$ (5) 式中：
$ {K}_{1} $ 为N区中的准渗透系数，$({{{\rm{m}}}}/{{{\rm{s}}}})^{n_1}$ (n1条件下)；$ {n}_{1} $ 为N区中的非达西湍流系数。当$ {n}_{1} < 1 $ 时渗流被称为前达西流；当$1 < {n}_{1}\leqslant 2$ 时为后达西流；当$ {n}_{1}=1 $ ，式(5)即为达西定律，此时$ {K}_{1} $ 为渗透系数[10]。模型的初始条件是：
$$ {h}_{1}\left(r, 0\right)=b $$ (6) 在定水头抽水试验中，代表完整井的第一类边界条件如下：
$$ \underset{r\to {r}_{{\rm{w}}}}{\mathrm{l}\mathrm{i}\mathrm{m}}{h}_{1}={h}_{{\rm{w}}} $$ (7) $$ \underset{r\to {r}_{{\rm{w}}}}{\mathrm{lim}}2\text{π} r{h}_{{\rm{w}}}\left({K}_{1}\frac{\partial {h}_{1}}{\partial r}\right)^{\tfrac{1}{{n}_{1}}}=-Q $$ (8) 式中：
$ Q $ 为抽水流量，$ {\mathrm{m}}^{3}/\mathrm{s} $ 。将式(5)代入式(1)可得：
$$ \frac{{\partial }^{2}{h}_{1}}{\partial {r}^{2}}+\frac{{n}_{1}}{r}\frac{\partial {h}_{1}}{\partial r}=\frac{{S}_{1}}{{h}_{1}}\frac{{n}_{1}}{{{K}_{1}}^{\tfrac{1}{{n}_{1}}}}{\left(\frac{\partial {h}_{1}}{\partial r}\right)}^{\tfrac{{n}_{1}-1}{{n}_{1}}}\frac{\partial {h}_{1}}{\partial t} $$ (9) 已知N区中的流量等于比流量乘以过流断面面积，如下式所述：
$$ -Q\approx A{q}_{1}=2\pi r{h}_{1}\left({K}_{1}\frac{\partial {h}_{1}}{\partial r}\right)^{\tfrac{1}{{n}_{1}}} $$ (10) 式中：
$ A=2\pi r{h}_{1} $ 为井壁处过流断面面积，$ {\mathrm{m}}^{2} $ ；在抽水井中$ {h}_{1} $ 为已知参数。假设非承压含水层的降深远远小于
$ b $ ，则可以基于Dupuit假设，在式(8)和式(9)中近似引入$\dfrac{{S}_{1}}{{h}_{1}}\dfrac{{n}_{1}}{{{K}_{1}}^{\tfrac{1}{{n}_{1}}}}\approx \dfrac{{S}_{1}}{b}\dfrac{{n}_{1}}{{{K}_{1}}^{\tfrac{1}{{n}_{1}}}}$ 来线性化一小部分推导过程。研究表明采用该近似方法可使建模过程大为简化，同时该近似方法所引起的计算误差相对较小，可以忽略不计[34-35]。因此，式(9)和式(10)可以近似改写为：$$ \frac{{\partial }^{2}{h}_{1}}{\partial {r}^{2}}+\frac{{n}_{1}}{r}\frac{\partial {h}_{1}}{\partial r}=\frac{{S}_{1}}{b}\frac{{n}_{1}}{{{K}_{1}}^{\tfrac{1}{{n}_{1}}}}{\left(\frac{\partial {h}_{1}}{\partial r}\right)}^{\tfrac{{n}_{1}-1}{{n}_{1}}}\frac{\partial {h}_{1}}{\partial t} $$ (11) $$ \frac{\partial {h}_{1}}{\partial r}\approx -\dfrac{{\left(\dfrac{Q}{2\pi rb}\right)}^{{n}_{1}}}{{K}_{1}} $$ (12) 结合式(11)与式(12)，N区中的控制方程可表示如下：
$$ \frac{{\partial }^{2}{h}_{1}}{\partial {r}^{2}}+\frac{{n}_{1}}{r}\frac{\partial {h}_{1}}{\partial r}={\varepsilon }_{1}{{r}}^{1-{n}_{1}}\frac{\partial {h}_{1}}{\partial t} $$ (13) 其中：
${\varepsilon }_{1}=\dfrac{{S}_{1}}{b}\dfrac{{n}_{1}}{{K}_{1}}{\left(\dfrac{{Q}}{2\pi b}\right)}^{{n}_{1}-1}$ 。1.2 达西区域
当
$ 1 < Re < 10 $ 时，可认为该区域处于D区，其抽水井流满足达西定律[25]，D区的比流量可用Darcian方程表示。$$ {q}_{2}={K}_{2}\frac{\partial {h}_{2}}{\partial r} $$ (14) 式中：
$ {K}_{2} $ 为D区中的渗透系数，m/s。初始条件为：
$$ {h}_{2}\left(r, 0\right)=b $$ (15) 无穷远处的边界条件为：
$$ {h}_{2}\left(\infty ,t\right)=b $$ (16) 将式(14)代入式(2)可得：
$$ \frac{{\partial }^{2}{h}_{2}}{\partial {r}^{2}}+\frac{1}{r}\frac{\partial {h}_{2}}{\partial r}=\frac{{S}_{2}}{{h}_{2}{K}_{2}}\frac{\partial {h}_{2}}{\partial t} $$ (17) 同理，为了简化建模过程，假设非承压含水层的降深远远小于
$ b $ ，基于Dupuit假设，引入近似关系$\dfrac{{S}_{2}}{{h}_{2}{K}_{2}}\approx \dfrac{{S}_{2}}{b{K}_{2}}$ 来线性化式(17)，可得D区中的控制方程可表示如下：$$ \frac{{\partial }^{2}{h}_{2}}{\partial {r}^{2}}+\frac{1}{r}\frac{\partial {h}_{2}}{\partial r}={\varepsilon }_{2}\frac{\partial {h}_{2}}{\partial t} $$ (18) 其中：
${\varepsilon }_{2}=\dfrac{{S}_{2}}{b{K}_{2}}$ 。2 基于有限差分的数值解
有限差分法简称FDM，常用于求解大型非线性地下水流运动问题，FDM在不规则含水层边界和含水层内区域参数的紧密空间近似方面具有灵活性[36]。因此，本研究采用FDM基本思想按时间步长和空间步长对定解区域进行网格划分[37]，把从无穷远处到抽水井处的空间区域及从抽水开始到抽水结束的时间区域用有限个网格节点代替。
采用FDM推求上述方程的数值解，需先对抽水井流模型进行如下两步预处理。第一步，针对无穷远边界，本文拟采用恒定不变且取值较大的半径值
${r}_{{\rm{e}}}$ 来代替，使抽水井流模型可被认为是一个因变量为半径$ r $ 和时间$ t $ 的一维非稳定流模型。第二步，将时间区间$ [0,t] $ 离散为$ N $ 个步长，径向空间$[{r}_{{\rm{w}}},{r}_{{\rm{e}}}]$ 离散为$ J $ 个步长。其中，任何一个时间坐标上节点$ m $ 处的时间值$ {t}_{m} $ 必满足$0\leqslant {t}_{m}\leqslant t$ ,$m=0, 1, \cdots ,N$ ，任何一个径向坐标上节点$ j $ 处的径向距离$ {r}_{j} $ 需满足${r}_{{\rm{w}}}\leqslant {r}_{j}\leqslant {r}_{{\rm{e}}}, j= 0, 1,\cdots ,J$ 。其步长大小的具体值表示为：$ \Delta t=t/N $ ，$\Delta r= \left({r}_{{\rm{e}}}-{r}_{{\rm{w}}}\right)/J$ 。由于在抽水井附近水位降深变化较大，因此在离抽水井较近的地方空间步长应尽可能小，而离抽水井较远的地方空间步长可以适当加大。对N区和D区的控制方程及其所有边界条件进行离散化，采用向后差分的方法，依次得出N区和D区的水头。2.1 非达西区域
N区控制方程式(13)为一个二阶齐次微分方程，基于FDM的基本理论，将控制方程离散化，得到离散化后的控制方程如下：
$$ {A}_{1j}{{h}_{1}}_{j-1}^{m+1}+{B}_{1j}{{h}_{1}}_{j}^{m+1}+{C}_{1j}{{h}_{1}}_{j+1}^{m+1}={{h}_{1}}_{j}^{m} $$ (19) 其中：
$$ {A}_{1j}=-\frac{\Delta t{({r}_{{\rm{w}}}+j\Delta r)}^{{n}_{1}-1}}{{\left(\Delta r\right)}^{2}{\varepsilon }_{1}} $$ (20) $$ {B}_{1j}=\frac{2\Delta t{({r}_{{\rm{w}}}+j\Delta r)}^{{n}_{1}-1}+{n}_{1}\Delta r\Delta t{({r}_{{\rm{w}}}+j\Delta r)}^{{n}_{1}-2}+{\varepsilon }_{1}{\left(\Delta r\right)}^{2}}{{\left(\Delta r\right)}^{2}{\varepsilon }_{1}} $$ (21) $$ {C}_{1j}=-\frac{\Delta t{({r}_{{\rm{w}}}+j\Delta r)}^{{n}_{1}-1}+{n}_{1}\Delta r\Delta t{({r}_{{\rm{w}}}+j\Delta r)}^{{n}_{1}-2}}{{\left(\Delta r\right)}^{2}{\varepsilon }_{1}} $$ (22) 对于N区的边界条件式(7)和式(8)，离散化后的边界条件如下：
$$ {{h}_{1}}_{j}^{0}=b $$ (23) $$ 2\pi b{r}_{{\rm{w}}}{\left({K}_{1}\frac{{{h}_{1}}_{1}^{m}-{{h}_{1}}_{0}^{m}}{\Delta r}\right)}^{\tfrac{1}{{m}_{1}}}=-Q+\pi {{r}_{{\rm{w}}}}^{2}\frac{{{h}_{1}}_{0}^{m+1}-{{h}_{1}}_{0}^{m}}{\Delta t} $$ (24) 式中：
${{h}_{1}}_{j}^{m}$ 为N区在时间节点为$m(m=0, 1,\cdots ,N)$ 、径向节点为$ j(j=0, 1,\cdots ,J) $ 处的水头。2.2 达西区域
同理，D区控制方程式(18)为一个二阶齐次微分方程，基于FDM的基本理论离散化后可得：
$$ {A}_{2j}{{h}_{2}}_{j-1}^{m+1}+{B}_{2j}{{h}_{2}}_{j}^{m+1}+{C}_{2j}{{h}_{2}}_{j+1}^{m+1}={{h}_{2}}_{j}^{m} $$ (25) 其中：
$$ {A}_{2j}=-\frac{\Delta t}{{\left(\Delta r\right)}^{2}{\varepsilon }_{2}} $$ (26) $$ {B}_{2j}=\frac{2\Delta t\left({r}_{{\rm{w}}}+j\Delta r\right)+\Delta r\Delta t+{\varepsilon }_{2}\left({r}_{{\rm{w}}}+j\Delta r\right){\left(\Delta r\right)}^{2}}{{\left(\Delta r\right)}^{2}{\varepsilon }_{2}({r}_{{\rm{w}}}+j\Delta r)} $$ (27) $$ {C}_{2j}=-\frac{\Delta t\left({r}_{{\rm{w}}}+j\Delta r\right)+\Delta t\Delta r}{{\left(\Delta r\right)}^{2}{\varepsilon }_{2}\left({r}_{{\rm{w}}}+j\Delta r\right)} $$ (28) 对于D区的边界条件式(15)和式(16)，离散化后的边界条件如下：
$$ {{h}_{2}}_{j}^{0}=b $$ (29) $$ {{h}_{2}}_{j}^{m}=b $$ (30) 式中：
${{h}_{2}}_{j}^{m}$ 为D区在时间节点为$m(m=0, 1,\cdots ,N)$ 、径向节点为$ j(j=0, 1,\cdots ,J) $ 处的水头。利用Matlab软件，将式(19)—式(24)与式(25)—式(30)分别导入，输入合适的步长，即可由软件自动计算并分别得出N区及D区各个节点上的流量及水头值。通过在
${r}_{j}\leqslant R$ 的径向节点上利用式(19)以及在$ {r}_{j} > R $ 的径向节点上利用式(25)计算得出的结果，调用Matlab中名为Equations and System Solver的内置模块，即可得到在假设的定水头试验条件下任意点的流量及水头情况，N区与D区求解矩阵分别如下：$$ \left(h_{10}^{m+1} \; h_{11}^{m+1} \; h_{12}^{m+1} \; \cdots \; h_{1 J}^{m+1} \; h_{1 J+1}^{m+1}\right)\left(\begin{array}{cccc} A_{11} & 0 & \cdots & 0 \\ B_{11} & A_{12} & \cdots & 0 \\ C_{11} & B_{12} & \cdots & 0 \\ 0 & C_{12} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & A_{1 J} \\ 0 & 0 & \cdots & A_{1 J} \\ 0 & 0 & \cdots & A_{1 J} \end{array}\right)=\left( h_{10}^m \; h_{11}^m \; \cdots \; h_{1J}^m\right) $$ (31) $$ \left( h_{20}^{m+1} \; h_{21}^{m+1}\; h_{22}^{m+1}\; \cdots\; h_{2 J}^{m+1}\; h_{2 J+1}^{m+1} \right)\left(\begin{array}{cccc} A_{21} & 0 & \cdots & 0 \\ B_{21} & A_{22} & \cdots & 0 \\ C_{21} & B_{22} & \cdots & 0 \\ 0 & C_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & A_{2 J} \\ 0 & 0 & \cdots & A_{2 J} \\ 0 & 0 & \cdots & A_{2 J} \end{array}\right)=\left( h_{20}^m \; h_{21}^m \; \cdots \; h_{2 J}^m \right) $$ (32) 定水头抽水试验可以使抽水井附近的水头在抽水早期迅速下降，并在短时间内接近
$ {h}_{\mathrm{w}} $ 。在此之后，水头变化的幅度会大大减小。因此，除了在抽水刚开始时的短时间内，其余阶段基于Dupuit假设提出的解是可以接受的。在定水头抽水试验中，由于流动的湍流程度会随着与泵井径向距离的增加而减弱，因此可将N区定义为
$ Re $ 大于临界雷诺数$ {Re}_{\mathrm{c}} $ 的区域，并且通过判断$ Re={Re}_{\mathrm{c}}{|}_{r\to R} $ 来确定该时刻的$ R $ 值。此外，还应该注意的是，所提出的解也可通过假设$ R={r}_{\mathrm{w}} $ 来评估纯达西流或设定$ R=\infty $ 调查纯非达西流的降深。计算雷诺数的公式[25]为$ Re=qd/\upsilon $ ，其中$ d $ 为介质的特征直径，m；$ \upsilon $ 为水的运动黏度，$ {\mathrm{m}}^{2}/\mathrm{s} $ ，其值通常为$1\times {10}^{6}\;{\mathrm{m}}^{2}/\mathrm{s}$ 。可以得到达西区和非达西区转换界面处的$ {Re}_{\mathrm{c}} $ 表达式如下：$$ {Re}_{\mathrm{c}}=\frac{d}{\upsilon }{K}_{2}\frac{{h}_{{2}{j+1}}-{h}_{{2}{j}}}{\Delta r} $$ (33) 根据相关研究，通常认为
$ {Re}_{\mathrm{c}}=10 $ 。基于计算得到的水头$ h $ 值，可以通过Matlab进而推求两区流界面位置$ R $ 值。拟通过计算得出的达西区水头$ {h}_{2} $ 进行迭代计算，其具体步骤如下：(1) 选取初始值。假设在抽水前，井流为达西流，即
${t}_{m}=0$ 时，时间步长数$m=0$ ，两区流交界面位置为$R={R}_{m}=0$ 。(2) 计算流速。计算第
$m+1$ 个时间步长${t}_{m+1}$ 时的水头${{h}_{2}}_{j}^{m+1}$ ，并计算每个节点处的流速及雷诺数$ Re $ ，然后根据临界雷诺数${Re}_{{\rm{c}}}$ 确定两区流交界面位置${R}_{m{\rm{c}}}$ 。(3) 收敛性判断。
$ \delta $ 为一个给定的判别标准，若$\left|{R}_{m}\right.\left.-{R}_{m{\rm{c}}}\right| < \delta$ ，则认为第$m$ 个时间步长的两区流交界面位置${R=R}_{m}$ 。然后计算下一时间步长两区流交界面位置：令$m=m+1$ ，若$m < N$ ，返回步骤(2)计算下一时刻的流速，反之则计算终止。若$\left|{R}_{m}\right.\left.-{R}_{m{\rm{c}}}\right| > \delta$ ，令${{R}_{m}=R}_{m{\rm{c}}}$ 。(4) 在新交界面处计算流速。根据计算得出的水头值，计算每个时间节点处的流速与雷诺数，由
${Re}_{{\rm{c}}}$ 确定两区流交界面位置${R}_{m{\rm{c}}}$ ，返回步骤(3)进行判断，进入下一步迭代计算，直至满足$ \delta < 5\text{%} $ 的精度要求，迭代终止。3 模型验证与对比
鉴于目前缺少非承压含水层定水头抽水计算模型的相关研究，难以找到相关工况完全一致的模型研究或实测案例，为验证本文提出有限差分数值模型的可靠性，本章将提出的非承压定水头抽水井流模型退化为承压定水头达西抽水井流模型，并与现有的经典承压定水头达西抽水井流解析模型(Jacob和Lohman模型)进行对比分析[4]。同时，通过假设案例，对所提出有限差分模型与基于COMSOL有限元模型的计算结果进行对比分析，以验证所提出的有限差分模型的可靠性。为方便计算，本文基于Matlab软件实现对提出模型的自动化求解。
3.1 与Jacob和Hantush解析解比较
在抽水过程中，通过控制含水层压力水头始终大于含水层厚度，将模型退化成承压含水层定水头抽水井流模型。然后，再令
$ {n}_{1}=1 $ 将模型退化成达西流，采用Jacob和Lohman模型的参数，将比流量$ q $ 及抽水时间$ t $ 无量纲处理并代入本模型中，可得到对比结果如图2所示。从图中可以看出，本文在定水头抽水的承压含水层达西流情况下的有限差分解与Jacob和Lohman模型解高度吻合，说明在模型推导过程中由数值差分带来的计算误差可忽略不计[38]。![]()
图 2 有限差分解与Jacob和Lohman解对比Figure 2. Comparison between the finite difference solution with Jacob and Lohman solution3.2 与COMSOL有限元解比较
非承压含水层非达西定水头抽水径流模型采用有限元模型基于COMSOL Multiphysics软件建立，该软件已被证明是最可靠的多物理场仿真软件之一[39]。COMSOL模型是基于N区的式(6)—式(9)及D区的式(15)—式(17)推导而来。
需要注意的是，式(16)中无穷远边界的距离在COMSOL解中用
$r=1\;000\;\mathrm{m}$ 表示，假设非承压含水层是砂质介质，基于现有研究可知非承压含水层的储水系数为$0.01 \sim 0.03$ ，非承压含水层中细砂的渗透系数为$1 \sim 5\;\mathrm{ }\mathrm{m}/\mathrm{d}$ [20]。因此，案例A参数选取见表1。考虑到Dupuit假设的要求，抽水井内部降深需要取一个相对较小的值，故取$b=12\;\mathrm{ }\mathrm{m}$ ,${h}_{\mathrm{w}}=9\;\mathrm{ }\mathrm{m}$ 。表1中的案例B、C和D用于下文讨论部分的参数选取。表 1 假设案例的参数值Table 1. Parameters of each hypothetical cases案例 $ {h}_{\mathrm{w}}/\mathrm{m} $ ${K}_{1}/{(\mathrm{m}\cdot\mathrm{d}^{-1})}^{ {{n} }_{1} }$ ${K}_{2}/(\mathrm{m}\cdot \mathrm{d}^{-1})$ $ {n}_{1} $ $ {S}_{1} $ $ {S}_{2} $ $ b/\mathrm{m} $ ${r}_{{\rm{w}}}/\mathrm{m}$ $ r/\mathrm{m} $ $ d/\mathrm{m}\mathrm{m} $ A 9 1.1 1.1 1.1 0.01 0.01 12 0.2 3 0.3 B 7,8,9 1.1 1.1 1.2 0.01 0.01 10 0.2 5 0.3 C 9 1.1 1.1 1.0,1.2,1.4 0.01 0.01 12 0.2 2 0.3 D 9 1.1 1.1 1.2 0.01 0.01 12 0.2 2 0.3 图3展示了两区流(
$ R $ 为变量)、纯达西流($R= {r}_{\mathrm{w}}=0.2\;\mathrm{m}$ )和纯非达西流($ R=\infty $ )的水头−时间曲线对比。纯非达西流COMSOL解近似为$R=1\;000\;\mathrm{m}$ 。应该注意到，$ {K}_{1} $ 和$ {S}_{1} $ 用于纯非达西流动的模拟，而$ {K}_{2} $ 和$ {S}_{2} $ 用于纯达西流动的模拟。这一对比清楚地表明，在所有3种流动条件下，模拟结果与有限差分解之间的偏差都很小。在抽水开始后第一个小时内出现一定偏差，可认为是忽略了所提出的有限差分解的垂直流动而引起，这种偏差随着时间的推移逐渐消失，模拟结果进一步验证了所提出的有限差分模型用于两区流动的水头模拟的有效性。![]()
图 3 不同R值下有限差分解与数值解的水头−时间曲线对比Figure 3. Comparison of head-time curves between finite difference solution and numerical solution for different R此外，两区流的水头明显小于纯非达西流的水头，而大于纯达西流的水头(图3)，说明非达西流对水头下降有负影响。这是因为
$ R $ 值越大表示非达西范围越大，抽水井附近的非承压含水层湍流范围更广。相比于纯达西流，纯非达西流的湍流程度更高。紊流会导致过水断面流速变大，相应比流量变大，含水层从外界得到水的能力增强，无穷远处对抽水井附近含水层补给变大，进而延缓因抽水导致的含水层水位下降[40]。因此，在抽水后期，纯非达西流的水头能够明显大于纯达西流的水头。图4展示了基于不同模型
$R\text{-}t$ 曲线的对比。结果表明，不同时间下本文提出的有限差分模型解与COMSOL有限元模型解拟合程度较高，验证了通过迭代法计算出的$ R $ 值的可靠性。而$ R $ 值随时间推移不断减小，说明非达西流对两区流动的影响随着抽水的继续而减弱。在抽水开始后，非承压含水层水头逐渐下降，接近于$ {h}_{\mathrm{w}} $ ，它进一步导致水力梯度和比流量的下降。比流量越小，流速越慢，雷诺数越小，最终导致N区消失。当抽水时间足够长时，N区可以完全消失，两区流可以转化为纯达西流。![]()
图 4 有限差分解界面距离−时间曲线与数值解的比较Figure 4. Interface distance-time curve of finite difference solution vs numerical solution3.3 与实测数据比较
L. Jones等[41]在美国威斯康星州开展了定水头抽水试验，所调查的非承压含水层由厚度为2.5 m的含泥量低、含砂量高的浅风化冰碛层组成。抽水井有效半径为0.051 m，且完全贯穿非承压含水层，其底部由砾石充填。抽水共进行24 h，抽水井内定水头为1.5 m[42]。两组降深数据分别在径向距离为0.87、3.62 m的观测井中采集。根据C. S. Chen等[2]的研究，水平方向的给水度和渗透系数分别为
$0.014 \sim 0.042$ 和$0.241 \sim 0.362\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}/\mathrm{d}$ 。因此，基于式(19)、式(25)经对比分析，选取参数$ {S}_{1}={S}_{2}=0.031 $ ，${K}_{1}=0.285\;\mathrm{ }\mathrm{m}/\mathrm{d}$ ，${K}_{2}=0.320\;\mathrm{m}/\mathrm{d}$ ，$ {n}_{1}=1.3 $ 时，可得到实测数据与本文模型解的对比曲线拟合度良好，如图5所示。但与图3的结果类似，有限差分解模拟曲线在抽水期间与实测数据的偏差很小，仅在抽水中期会出现一定程度的误差，且这种偏差随着时间的增加逐渐消失。图5进一步验证了所提出的有限差分模型在实际工程中的适用性。4 讨 论
在实际抽水试验中，一般认为抽水附近区域流态为达西流，以方便概化模型及预测井内降深。本研究提出的定水头抽水两区流数值解，将抽水井附近的有限区域视为非达西区域，远离抽水井的半无限区域视为达西区域，从而将非承压含水层的地下水渗流根据流态区分为2个流域。相比于纯达西流模型，该两区流模型更符合实际工况，计算结果精度更高，但相应计算过程也较复杂。该模型适用于沙土或碎石比例较大、渗透性较强的含水层，其渗流流速相对较快，
$ Re $ 相对较高，在抽水井附近极易形成非达西流。但在黏土等渗透能力弱的含水层，渗流流速较小，非达西范围相对小，该模型将不再适用。在实际工程应用中，通常会着重观测井内水头及抽水量的动态变化。因此，该模型有必要讨论偏离达西程度对流态的影响，进而判断不同非达西湍流系数$ {n}_{1} $ 对水头及抽水井抽水量的影响，并讨论井内水头$ {h}_{\mathrm{w}} $ 对抽水井抽水速率的影响。准确评估抽水流量$ Q $ 的动态变化对于设计定水头抽水试验至关重要。在下文中，基于假设案例B、C和D(参数见表1)揭示非达西流的$ {n}_{1} $ 、两区流转换界面位置$ R $ 和井内水头$ {h}_{\mathrm{w}} $ 对抽水速率动态变化的影响，以及径向距离对断面流量的影响。图6展示了案例B中不同时刻
$ {h}_{\mathrm{w}} $ 的$Q\text{-}t$ 曲线的模拟结果。结果表明，固定时刻下$ Q $ 随着$ {h}_{\mathrm{w}} $ 的增加不断降低，在抽水早期不同曲线对应的$ Q $ 之间偏差较大，随着抽水的继续，流量偏差逐渐减小。$ Q $ 的变化与$ {h}_{\mathrm{w}} $ 呈负相关关系，在抽水早期$Q\text{-}t$ 曲线下降速度较快，随着抽水的继续$ Q $ 下降的趋势越来越平缓，符合$ {h}_{\mathrm{w}} $ 的小幅降低会使流速大幅度升高的事实。流量变化实际上是由抽水井内部水头下降而引起的水力梯度变化造成的，因此，较大的$ {h}_{\mathrm{w}} $ 意味着较小的比流量，从而导致抽水流量较小。![]()
图 6 泵井内不同水头下流量的动态变化Figure 6. Dynamic development of the flow rates with different hydraulic head inside the pumping well图7展示了案例C中非达西湍流系数对泵送流量的影响，其中图7a和图7b分别为固定径向距离
$r=2\;\mathrm{m}$ 时的$Q\text{-}t$ 曲线和相同抽水时间$ t=1\;\mathrm{d} $ 时的$h\text{-}r$ 曲线。从图7a中可以发现，在抽水早期，固定时间点的$ Q $ 随着$ {n}_{1} $ 的增加而增加。而随着抽水的继续，不同$ {n}_{1} $ 对应的流量偏差不断减小。从理论上讲，这是可以理解的。抽水井附近含水层补给水量以重力释水为主，如式(10)所示，假设其等于比流量乘以过流断面面积，比流量越大，抽水流量越大。由此可见，在抽水早期，抽水井附近区域含水层的重力储水逐渐释放并补充到抽水井内，$ {n}_{1} $ 越大，湍流程度越大，导致非承压含水层储水的释放越快。随着抽水的继续，N区重力储水的释放逐渐趋于稳定，抽水井中的水主要通过D区的水流进行补给，导致$Q\text{-}t$ 曲线逐渐趋同。![]()
图 7 不同非达西湍流系数对泵送流量的影响Figure 7. Effects of different non-Darcy coefficients on pumping flow图7b体现了
$ t=1\;\mathrm{d} $ 时不同径向距离的非达西湍流的发展特征。结果表明，固定$ r $ 下$ h $ 随着$ {n}_{1} $ 的增大而增大，在$ r $ 较小时不同曲线对应的$ h $ 之间的偏差较大，随着$ r $ 的增大，$ h $ 偏差逐渐减小。非达西系数与水头之间存在正相关关系，在$ r $ 较小时$h\text{-}r$ 曲线上升的速度较快，随着$ r $ 的增大，曲线斜率逐渐减小。$ {n}_{1} $ 越大，表示水力梯度越大，对应流量越大。在抽水井附近，水头波动较大，但没有出现不合理的突变，证明了所提有限差分解的适用性。图8展示了径向距离对断面流量的影响，其中图8a和图8b分别用于表示抽水前期 (
$t=0.1\;\mathrm{d}$ ) 和抽水后期 ($t=100.0\;\mathrm{d}$ ) 的${Q}_{1}\text{-}r$ 曲线(Q1为断面流量)，其水力参数见表1的案例D。对比图8a和图8b可以发现，$t=0.1\;\mathrm{d}$ 时$ {Q}_{1} $ 随$ r $ 的增加而快速下降，其$ {Q}_{1} $ 较$t=100.0\;\mathrm{d}$ 时更大。随着抽水时间的增加，$ R $ 逐渐变小，该结论符合图4中$R\text{-}t$ 曲线的趋势。$ {Q}_{1} $ 随着$ r $ 的增大不断减小，在N区部分$ {Q}_{1} $ 曲线随距离下降的速度较快，并在经过转换界面$ R $ 后出现突变，在D区$ {Q}_{1} $ 曲线随距离下降的速率明显变慢。从理论上来说，定水头抽水过程中，在井壁处存在$ {Q}_{1} $ 变化的最大值，随着$ r $ 的增大，距离抽水井越远，对抽水井补给量越小，断面流量越小。在抽水前中期，含水层中距离抽水井较远的部分无法及时对抽水井附近进行充足的补给，抽水井附近区域含水层的重力储水会逐渐释放并补充到抽水井内。随着抽水时间的增长，在抽水后期，远处的水流逐渐补充到抽水井附近，抽水井周边区域的水位和断面流量的下降速度逐渐变小。相比于N区，D区水力梯度更小，其比流量随径向距离的下降速率更慢，其${Q}_{1}\text{-}r$ 曲线的转折点出现在$ r=R $ 处。![]()
图 8 不同时间下径向距离对断面流量的影响Figure 8. Influence of radial distance on cross-section flow at different time5 结 论
a. 所提出的有限差分解可用于在由定水头抽水试验引起的两区流、纯达西流和纯非达西流的情况下进行降深模拟。在两区流动的情况下，所提出的解也可用于评价抽水流量和非达西区域的动态发展。
b. 忽略非达西区域的影响可能会给两区流动的水头模型带来误差。有限非达西区域的湍流会分别导致两区流中水头较纯非达西流和纯达西流的水头偏大和偏小，且随抽水时间的增加逐渐变大。
c. 利用Izbash方程，发现瞬态抽水速率−时间变化关系可以由抽水井内水头
${h}_{{\rm{w}}}$ 和非达西滞流系数$ {n}_{1} $ 显著控制。在抽水过程中，减小${h}_{{\rm{w}}}$ 或增大$ {n}_{1} $ 可以提高抽水速率，这种影响随着抽水时间的增加会不断减弱并在抽水结束时消失。d. 在抽水过程中，
$ {Q}_{1} $ 随着$ r $ 的增大不断减小，并在$ r=R $ 处出现转折点，具体表现为：在N区中$ {Q}_{1} $ 曲线随距离下降的速度较快，在经过转换界面$ R $ 处后出现突变，在D区中$ {Q}_{1} $ 曲线随距离下降的速度明显变慢。 -
![]()
图 1 考虑两区转换的非承压含水层定水头抽水模型
rw—抽水井有效半径，m；hw—抽水井中水头，m；R—N区和D区转换界面到抽水井中心的距离，m；b—含水层厚度，m
Fig. 1 Schematic of the constant-head test in an unconfined aquifer considering the two-region transform
![]()
图 2 有限差分解与Jacob和Lohman解对比
Fig. 2 Comparison between the finite difference solution with Jacob and Lohman solution
![]()
图 3 不同R值下有限差分解与数值解的水头−时间曲线对比
Fig. 3 Comparison of head-time curves between finite difference solution and numerical solution for different R
![]()
图 4 有限差分解界面距离−时间曲线与数值解的比较
Fig. 4 Interface distance-time curve of finite difference solution vs numerical solution
![]()
图 6 泵井内不同水头下流量的动态变化
Fig. 6 Dynamic development of the flow rates with different hydraulic head inside the pumping well
![]()
图 7 不同非达西湍流系数对泵送流量的影响
Fig. 7 Effects of different non-Darcy coefficients on pumping flow
![]()
图 8 不同时间下径向距离对断面流量的影响
Fig. 8 Influence of radial distance on cross-section flow at different time
表 1 假设案例的参数值
Table 1 Parameters of each hypothetical cases
案例 $ {h}_{\mathrm{w}}/\mathrm{m} $ ${K}_{1}/{(\mathrm{m}\cdot\mathrm{d}^{-1})}^{ {{n} }_{1} }$ ${K}_{2}/(\mathrm{m}\cdot \mathrm{d}^{-1})$ $ {n}_{1} $ $ {S}_{1} $ $ {S}_{2} $ $ b/\mathrm{m} $ ${r}_{{\rm{w}}}/\mathrm{m}$ $ r/\mathrm{m} $ $ d/\mathrm{m}\mathrm{m} $ A 9 1.1 1.1 1.1 0.01 0.01 12 0.2 3 0.3 B 7,8,9 1.1 1.1 1.2 0.01 0.01 10 0.2 5 0.3 C 9 1.1 1.1 1.0,1.2,1.4 0.01 0.01 12 0.2 2 0.3 D 9 1.1 1.1 1.2 0.01 0.01 12 0.2 2 0.3 
下载: 导出CSV
-
[1] MARKLE J M,ROWE R K,NOVAKOWSKI K S. A model for the constant–head pumping test conducted in vertically fractured media[J]. International Journal for Numerical and Analytical Methods in Geomechanics,1995,19(7):457−473. DOI: 10.1002/nag.1610190702
[2] CHEN C S,CHANG C C. Well hydraulics theory and data analysis of the constant–head test in an unconfined aquifer with the skin effect[J]. Water Resources Research,2003,39(5):1121−1135.
[3] 柯瀚,胡杰,吴小雯,等. 竖井抽水下垃圾填埋场渗滤液运移规律研究[J]. 岩土工程学报,2018,40(5):786−793. KE Han,HU Jie,WU Xiaowen,et al. Investigation into leachate transport in MSW landfills under pumping of vertical wells[J]. Chinese Journal of Geotechnical Engineering,2018,40(5):786−793.
[4] JACOB C E,LOHMAN S W. Nonsteady flow to a well of constant drawdown in an extensive aquifer[J]. Eos,Transactions American Geophysical Union,1952,33(4):559−569. DOI: 10.1029/TR033i004p00559
[5] CHANG Y C,YEH H D,CHEN G Y. Transient solution for radial two–zone flow in unconfined aquifers under constant–head test[J]. Hydrological Processes,2010,24(11):1496−1503. DOI: 10.1002/hyp.7610
[6] DAN Hancheng,HE Linhua,XU Bo. Experimental investigation on non−Darcian flow in unbound graded aggregate material of highway pavement[J]. Transport in Porous Media,2016,112(1):189−206. DOI: 10.1007/s11242-016-0640-z
[7] 王军辉,王峰. 论抽水的降落漏斗范围、影响半径与环境影响范围[J]. 水利学报,2020,51(7):827−834. WANG Junhui,WANG Feng. Discussion on the range of groundwater depression cone,radius of influence and scope of environmental impacts during pumping[J]. Journal of Hydraulic Engineering,2020,51(7):827−834.
[8] LIU Zhongyu,XIA Yangyang,SHI Mingsheng,et al. Numerical simulation and experiment study on the characteristics of non−Darcian flow and rheological consolidation of saturated clay[J]. Water,2019,11(7):1385. DOI: 10.3390/w11071385
[9] ZHUANG Chao,ZHOU Zhifang,ZHAN Hongbin,et al. New graphical methods for estimating aquifer hydraulic parameters using pumping tests with exponentially decreasing rates[J]. Hydrological Processes,2019,33(17):2314−2322. DOI: 10.1002/hyp.13470
[10] LI Yabing,ZHOU Zhifang,ZHUANG Chao,et al. Non–Darcian effect on a variable–rate pumping test in a confined aquifer[J]. Hydrogeology Journal,2020,28(8):2853−2863. DOI: 10.1007/s10040-020-02223-w
[11] 刘元会,常安定,邓秋霞. 非线性井流问题水头降深的分区研究[J]. 中国农村水利水电,2007(10):11−12. LIU Yuanhui,CHANG Anding,DENG Qiuxia. Study on sub–area water head drawdown of nonlinear well flow[J]. China Rural Water and Hydropower,2007(10):11−12.
[12] YEH H D,CHANG Yachi. Recent advances in modelling of well hydraulics[J]. Advances in Water Resources,2013,51:27−51. DOI: 10.1016/j.advwatres.2012.03.006
[13] XIAO Liang,YE Ming,XU Yongxin,et al. A simplified solution using Izbash’s Equation for non–Darcian flow in a constant rate pumping test[J]. Groundwater,2019,57(6):962−968. DOI: 10.1111/gwat.12886
[14] 刘元会,常安定,邓秋霞. 线性非线性流并存区域井流问题的水头降深研究[J]. 西北农林科技大学学报(自然科学版),2005,33(3):113−115. LIU Yuanhui,CHANG Anding,DENG Qiuxia. Drawdown of well flow in the co−existed linear and nonlinear exponents[J]. Journal of Northwest Sci−Tech University of Agriculture and Forestry (Natural Science Edition),2005,33(3):113−115.
[15] QUINN P M,PARKER B L,CHERRY J A. Validation of non−Darcian flow effects in slug tests conducted in fractured rock boreholes[J]. Journal of Hydrology,2013,486:505−518. DOI: 10.1016/j.jhydrol.2013.02.024
[16] ALI R,MCFARLANE D,VARMA S,et al. Potential climate change impacts on the water balance of regional unconfined aquifer systems in south−western Australia[J]. Hydrology and Earth System Sciences Discussions,2012,9(5):6367−6408.
[17] NARASIMHAN T N,ZHU Ming. Transient flow of water to a well in an unconfined aquifer:Applicability of some conceptual models[J]. Water Resources Research,1993,29(1):179−191. DOI: 10.1029/92WR01959
[18] BOULTON N S. Analysis of data from non–equilibrium pumping tests allowing for delayed yield from storage[J]. Proceedings of the Institution of Civil Engineers,1963,26(3):469−482. DOI: 10.1680/iicep.1963.10409
[19] DUPUIT J. Études théoriques et pratiques sur le mouvement des eaux dans les canaux découverts et á travers les terrains perméablesm (2ème edition)[M]. Libraire des corps impérious des ponts et chaussées et des mines, 1863.
[20] KRUSEMAN G P, RIDDER N A D, VERWEIJ J M. Analysis and evaluation of pumping test data[J]. ILRI, 1990, 47: 1−377.
[21] 陈崇希, 林敏, 成建梅. 地下水动力学[M]. 北京: 地质出版社, 2011. [22] BRESCIANI E,DAVY P,DREUZY J R. Is the Dupuit assumption suitable for predicting the groundwater seepage area in hillslopes?[J]. Water Resources Research,2014,50(3):2394−2406. DOI: 10.1002/2013WR014284
[23] DAGAN G. A method of determining the permeability and effective porosity of unconfined anisotropic aquifers[J]. Water Resources Research,1967,3(4):1059−1071. DOI: 10.1029/WR003i004p01059
[24] NEUMAN S P. Analysis of pumping test data from anisotropic unconfined aquifers considering delayed gravity response[J]. Water Resources Research,1975,11(2):329−342. DOI: 10.1029/WR011i002p00329
[25] WEN Zhang,HUANG Guanhua,ZHAN Hongbin,et al. Two–region non–Darcian flow toward a well in a confined aquifer[J]. Advances in Water Resources,2008,31(5):818−827. DOI: 10.1016/j.advwatres.2008.01.014
[26] LI Yabing,ZHOU Zhifang,SHEN Qi,et al. Two–region Darcian and non–Darcian flow towards a well with exponentially decayed rates considering time–dependent critical radius[J]. Journal of Hydrology,2021,601:126712. DOI: 10.1016/j.jhydrol.2021.126712
[27] 文章,黄冠华,刘壮添,等. 越流含水层中抽水井附近非达西流两区模型近似解析解[J]. 地球科学(中国地质大学学报),2011,36(6):1165−1172. WEN Zhang,HUANG Guanhua,LIU Zhuangtian,et al. An approximate analytical solution for two–region non–Darcian flow toward a well in a leaky aquifer[J]. Earth Science (Journal of China University of Geosciences),2011,36(6):1165−1172.
[28] CHANG Yachi,YEH H D. Skin effect in generalized radial flow model in fractured media[J]. Geophysical Journal International,2011,185(1):78−84. DOI: 10.1111/j.1365-246X.2011.04943.x
[29] 王全荣,唐仲华,文章,等. 越流含水层抽水井附近非达西流与达西流区界面位置变化规律研究[J]. 水利学报,2012,43(10):1171−1178. WANG Quanrong,TANG Zhonghua,WEN Zhang,et al. Numeric simulation for flow to a pumping well with moving boundary of the non–Darcian flow region in a leaky aquifer[J]. Journal of Hydraulic Engineering,2012,43(10):1171−1178.
[30] BEAR J. Hydraulics of groundwater[M]. New York: Dover Publications, 1979.
[31] MATHIAS S A,BUTLER A P,ZHAN Hongbin. Approximate solutions for Forchheimer flow to a well[J]. Journal of Hydraulic Engineering,2008,134(9):1318−1325. DOI: 10.1061/(ASCE)0733-9429(2008)134:9(1318)
[32] WEN Zhang, HUANG Guanhua, ZHAN Hongbin. Non–Darcian flow in a single confined vertical fracture toward a well[J]. Journal of Hydrology, 2006, 330(3/4): 698–708.
[33] BORDIER C, ZIMMER D. Drainage equations and non−Darcian modelling in coarse porous media or geosynthetic materials[J]. Journal of Hydrology, 2000, 228(3/4): 174–187.
[34] MALAMA B,KUHLMAN K L,REVIL A. Theory of transient streaming potentials associated with axial–symmetric flow in unconfined aquifers[J]. Geophysical Journal International,2009,179(2):990−1003. DOI: 10.1111/j.1365-246X.2009.04336.x
[35] MOENCH A F,PRICKETT T A. Radial flow in an infinite aquifer undergoing conversion from artesian to water table conditions[J]. Water Resources Research,1972,8(2):494−499. DOI: 10.1029/WR008i002p00494
[36] 卢洪健,王卓然. 地下水模拟方法与应用软件研究进展[J]. 地下水,2022,44(6):49−52. LU Hongjian,WANG Zhuoran. Research progress in groundwater simulation methods and application softwares[J]. Groundwater,2022,44(6):49−52.
[37] 孙广才. 水力学模型及其数值求解[J]. 渭南师范学院学报,2015,30(22):17−20. SUN Guangcai. Hydraulics model and its numerical solution[J]. Journal of Weinan Normal University,2015,30(22):17−20.
[38] 文章,黄冠华,李健,等. 越流含水层中抽水井附近非达西流动模型的数值解[J]. 水动力学研究与进展,2009,24(4):448−454. WEN Zhang,HUANG Guanhua,LI Jian,et al. A numerical solution for non–Darcian flow toward a well in a leaky aquifer[J]. Chinese Journal of Hydrodynamics,2009,24(4):448−454.
[39] CHUI T F M, FREYBERG D L. The use of COMSOL for integrated hydrological modelling[R]. COMSOL Conference, 2007.
[40] 刘凯. 非完整井附近非达西渗流规律研究[D]. 武汉: 中国地质大学(武汉), 2014. LIU Kai. Investigation on non–Darcian flow toward a partially penetrating well[D]. Wuhan: China University of Geosciences (Wuhan), 2014.
[41] JONES L,LEMAR T,TSAI C–T. Results of two pumping tests in Wisconsin age weathered till in Iowa[J]. Groundwater,1992,30(4):529−538. DOI: 10.1111/j.1745-6584.1992.tb01529.x
[42] WEN Zhang,HUANG Guanhua,ZHAN Hongbin. An analytical solution for non−Darcian flow in a confined aquifer using the power law function[J]. Advances in Water Resources,2008,31(1):44−55. DOI: 10.1016/j.advwatres.2007.06.002
-
期刊类型引用(2)
1. 史鹏钰，韦宇杰，肖良，陈立华. 非达西效应下瞬态承压-无压转换数值解. 广西大学学报(自然科学版). 2024(03): 517-528 . 
百度学术
2. 滕开庆，刘健军，史鹏钰，宗一杰，肖良. 非达西效应下不稳定承压-无压井流的有限差分解. 安全与环境工程. 2024(06): 179-188 . 百度学术
其他类型引用(0)
计量
- 文章访问数: 183
- HTML全文浏览量: 7
- PDF下载量: 30
- 被引次数: 2
