考虑井储效应的双重介质单孔抽水非稳定流解析模型

Identification of Hydrogeological Parameter in Double⁃Porosity Fractured Aquifer Based on Single⁃Well Pumping Test

摘要

为探究双重介质含水层中单孔抽水试验水文地质参数求取方法，在不考虑井储效应的双重介质非达西流模型基础上，提出了考虑井储效应的双重介质非达西流模型，并以某野外试验场中实地进行的单孔抽水试验为例，比较了不同求参方法对单孔抽水试验参数识别的效果.结果表明，使用考虑井储的双重介质非达西流模型，反演结果与实测水位降深的吻合度最高；考虑非达西流态，更加适用于抽水流量大，地下水流场呈现非达西流的状态；考虑井储效应，能够更好地刻画出抽水过程前期地下水位迅速下降的现象；考虑双重介质结构，能够更好地刻画出抽水过程前‒中期孔隙‒裂隙介质间水量交换作用由弱变强的现象.说明本文提出的考虑井储效应的双重介质非达西流模型具有一定的适用性.

Abstract

To find out a suitable equation for single-well pumping test in double-porosity aquifer, it takes wellbore storage effect into consideration and present the double-porosity non-Darian model with wellbore storage effect based on double-porosity non-Darcian model in this study. Taking a single-well pumping test as an example, the results show that the simulation curve has the highest accuracy when using the double-porosity non-Darcian model proposed in this study; taking non-Dacian effect into consideration is more suitable when the pumping rate is relatively large and the flow in aquifer shows non-Darcian property. When considering non-Darcian effect, we can have a better inversion of the instantaneous drawdown in early time; when regarding the aquifer as double-porosity, we can portray the weak to strong process of the water exchange in fracture and porosity media. In that case, the double-porosity non-Darcian model with wellbore storage effect presented in this study is suitable for single-well pumping test in double-porosity aquifer.

Graphical abstract

关键词

单孔抽水 / 双重介质 / 非达西流 / 井储效应 / 工程地质 / 水文地质.

Key words

single⁃well pumping test / double⁃porosity / non⁃Darcian flow / wellbore storage effect / engineering geology / hydrogeology

0 引言

在工程勘察中，地下水井流试验，是确定含水层水力联系与水文地质参数最直接的方法，同时也为地下水资源的评价和开发、污染物运移预测以及地下水修复等提供有效的科学依据（李娜， 2020； 李旭等，2022），根据地下水井流试验获取的含水层水位随时间变化曲线，并结合水文地质条件选择相应的井流模型即可求取含水层的水文地质参数（李丞等， 2020）.然而，面对多变的野外现场环境和错综复杂的水文地质条件，地下水井流试验常面临着多重挑战，例如，由于条件所限，无法布设观测孔，只能进行单孔抽水试验；含水层呈现非均质态，而传统抽水试验求解含水层渗透系数时，将含水层介质等效为均质介质（黄康等，2022）；受多因素扰动的影响，抽水井附件呈现非达西流状态等；同时，抽水孔中原先存储水量的释放也可能对抽水试验造成影响，因此，为尽可能准确的对抽水试验结果进行拟合，获取符合含水层结构特征的水文地质参数，探究综合考虑了井储效应和非达西流效应的双重介质含水层井流模型是有必要的.

1856年，法国工程师Darcy通过砂柱试验提出了水在孔隙介质中渗透的定律，即Darcy定律.以此为开端开始出现对井流理论的相关研究，并先后提出了在圆形透水层中心抽水的潜水井和承压水井稳定井流Dupuit模型；在无界承压含水层中用一个影响半径R来表达抽水过程影响范围的Thiem模型；在无限承压含水层中定流量抽水的地下水完整井非稳定井流Theis模型.随着研究的深入，地下水流态、抽水井类型及含水层性质对井流模型的影响也逐渐被研究者所重视，例如，Forchheimer（1901） 、Izbash （1931）认为，随着渗透流速的增大，地下水运动存在偏离达西线性渗流方程的现象，而万军伟等（2013）、Li et al.（2019）则通过室内实验证明，无论是低速流或高速流，均不存在达西线性渗流；Fenske（1977）在Theis模型的基础上，考虑了井储效应对水位降深的影响，给出了考虑井储效应的地下水呈线性渗流状态下的井流模型；Wen et al. （2008）、Wen and Wang（2013）则考虑了井储效应、表皮效应与相邻含水层的越流作用，使用Izbash公式将抽水过程中的地下水流场概化为达西‒非达西流两区模型，探究了水位降深与抽水时间的关系.

值得注意的是，上述模型大都认为主含水层为单一介质且呈现均质各项同性的特征，对于受构造作用或成岩作用影响而呈现高度非均质性特征的裂隙含水层抽水试验的适用性有限.以火成岩地区碎屑岩含水层为例，由于受到成岩作用和后期构造作用的影响，储水体往往由成岩裂隙‒孔洞、一般构造裂隙及大型构造裂隙或剪切破碎带等空间分布不连续、含（导）水性能差异大的非均匀双重含水介质组成.针对这类含水层，Barenblatt et al.（1960）、Warren and Root（1963）、Moench（1984）等将其概化为具有低渗透性高储水性的孔隙介质和高渗透性低储水性的裂隙介质共同构成的双重介质模型，并给出了各自的解析解表达式；何林青等（2015）使用分段Theis模型的方法，对双重介质含水层的释水机理及不同时期的水位降深曲线所反映的含水介质特性进行了定性描述和定量分析；Wang et al.（2021）则考虑到抽水过程中非达西流的影响，提出了不考虑井储条件的双重介质达西流模型.

然而，由于未考虑井储效应，Wang模型对多孔抽水试验的适用性存在一定限制.实际勘察工作中，由于受多因素影响，常常只具备开展单孔抽水试验的条件.受井储效应的影响，单孔抽水试验初期抽水孔水位极速下降，亟需寻找一种能够描述该情况的解析解模型.本文则在Wang模型的基础上，考虑了井储效应对水位降深的影响，为双重介质含水层单孔抽水试验水文地质参数的求取提供了参考.

1 考虑井储效应的双重介质非达西流模型

1.1　数学模型

如附图1所示，认为含水层具有承压性，且水平等厚、侧向无限延伸；孔隙介质和裂隙介质连续体均为均质各向同性；抽水井为完整井，抽水井处地下水无垂向流动，抽水井具有一定的半径r_w 且套管半径为r_c （注：r_c ≥r_w ）；抽水过程流量恒定为Q（m³/h），抽水前，承压含水层水头面保持水平静止，且抽水期间地下水的垂向流速可忽略.

基于以上假设条件，可写出双重介质含水层中裂隙介质的流场控制方程：

\frac{\partial q\left(r,t\right)}{\partial r}+\frac{q\left(r,t\right)}{r}=\frac{S_f}{b}\frac{\partial s_f\left(r,t\right)}{\partial t}+q_a,(1)

式中，q（r，t）为径向渗流速度 ［L/T］；S_f 为裂隙介质储水系数［-］；为裂隙介质含水层水位降深［L］；q_a 为孔隙‒裂隙介质间的水量交换量［L/T］；b为承压含水层厚度［L］；t为抽水时间［T］；r为距抽水井径向距离（r≥r_w ）［L］.

在双重介质模型假设中，认为孔隙块体介质地下水流呈现达西流，因此孔隙块体介质有如下流场控制方程：

K_m{\nabla }^{\mathrm{2}}{s}_m\left(r,t\right)=\frac{S_m}{b}\frac{\partial s_m\left(r,t\right)}{\partial t}-q_a,(2)

式中，K_m 为孔隙块体介质渗透系数［L/T］；s_m （r，t）为孔隙块体介质水位降深［L］.由于K_m 极小，因此式（2）左端可视为0（De Smedt， 2011）.

孔隙‒裂隙介质间的流量交换符合拟稳态渗流假定，因而有控制方程：

q_a=\alpha K_m\left[s_f\left(r,t\right)-s_m\left(r,t\right)\right],(3)

式中，α为形状系数［L^-2］，该系数与孔隙介质的几何形状有关，且控制着孔隙介质和裂隙介质间的流量交换.

由于式（2）左端为0，得到q_a 的近似表达式，将该表达式代入式（1）、式（3）得到：

\frac{\partial q\left(r,t\right)}{\partial r}+\frac{q\left(r,t\right)}{r}=\frac{S_f}{b}\frac{\partial s_f\left(r,t\right)}{\partial t}+\frac{S_m}{b}\frac{\partial s_m\left(r,t\right)}{\partial t},(4)

S_m\frac{\partial s_m\left(r,t\right)}{\partial t}=C\left[s_f\left(r,t\right)-s_m\left(r,t\right)\right],(5)

式中，C=αK_mb，为孔隙‒裂隙介质间的水量交换系数［T^-1］.

初始情况下，孔隙块体介质和裂隙介质的水位降深均为0：

s_f\left(r,\mathrm{0}\right)=s_m\left(r,\mathrm{0}\right)=\mathrm{0},(6)

抽水过程中，无限远处孔隙块体介质和裂隙介质的水位降深始终为0：

s_f\left(\mathrm{\infty },t\right)=s_m\left(\mathrm{\infty },t\right)=\mathrm{0},(7)

考虑井储效应时，抽水井处的控制方程为：

\underset{r\to r_w}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{2}\mathrm{\pi }rbq\left(r,t\right)-\mathrm{\pi }{r}_c^{\mathrm{2}}\frac{\partial s_w\left(t\right)}{\partial t}=-Q,(8)

式中，r_w 为抽水井半径［L］，r_c 为套管半径［L］，s_w 为抽水井水位降深［L］，Q为抽水流量（流出为负，流入为正）［L³/T］.

根据Bordier and Zimmer（2000）的研究成果，由于与达西公式在形式上相似，使用Izbash公式描述抽水过程中的非达西流状态较为便捷，有：

q\left(r,t\right)=K_q^{\mathrm{1}/n}{\left[\frac{\partial s_f\left(r,t\right)}{\partial r}\right]}^{\mathrm{1}/n},(9)

式中，K_q 和n均为常量，K_q 为准渗透系数［（L/T） ⁿ ］，n为经验参数，取值范围在1~2，当n为1时式（9）变为达西流，K_q 变为渗透系数.

将式（9）代入式（4）、式（8），有：

\begin{array}{l}\frac{{\partial }^{\mathrm{2}}{s}_f}{\partial r^{\mathrm{2}}}+\frac{n}{r}\frac{\partial s_f}{\partial r}=\frac{S_f}{b}\frac{n}{K_q^{\mathrm{1}/n}}{\left(\frac{\partial s_f}{\partial r}\right)}^{\frac{n-\mathrm{1}}{n}}\frac{\partial s_f}{\partial t}+\\ \frac{S_m}{b}\frac{n}{K_q^{\mathrm{1}/n}}{\left(\frac{\partial s_f}{\partial r}\right)}^{\frac{n-\mathrm{1}}{n}}\frac{\partial s_m}{\partial t}\end{array},(10)

\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q^{\mathrm{1}/n}{\left[\frac{\partial s_f\left(r,t\right)}{\partial r}{\left|\right.}_{r\to r_w}\right]}^{\mathrm{1}/n}-\mathrm{\pi }{r}_c^{\mathrm{2}}\frac{\partial s_w\left(t\right)}{\partial t}=-Q,(11)

式（10）、式（11）中存在非线性项，根据Wen et al.（2008）提出的与时间变量相独立的线性化方法替代式（10）中由\partial s_f/\partial r构成的非线性项，有：

\frac{\partial s_f\left(r,t\right)}{\partial r}=\frac{{\left[q\left(r,t\right)\right]}^n}{K_q}\approx -\frac{{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }rb}\right)}^n}{K_q},(12)

将式（12）代入式（10）、式（11），有：

\begin{array}{l}\frac{{\partial }^{\mathrm{2}}{s}_f}{\partial r^{\mathrm{2}}}+\frac{n}{r}\frac{\partial s_f}{\partial r}\approx \frac{S_f}{b}\frac{n}{K_q}{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right)}^{n-\mathrm{1}}{r}^{\mathrm{1}-n}\frac{\partial s_f}{\partial t}+\\ \frac{S_m}{b}\frac{n}{K_q}{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right)}^{n-\mathrm{1}}{r}^{\mathrm{1}-n}\frac{\partial s_m}{\partial t}\end{array},(13)

\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b{r}_w}\right]}^{\left(\mathrm{1}-n\right)}\frac{\partial s_f\left(r,t\right)}{\partial r}{\left|\right.}_{r\to r_w}-\\ \mathrm{\pi }{r}_c^{\mathrm{2}}\frac{\partial s_w\left(t\right)}{\partial t}\approx -Q\end{array}.(14)

1.2　模型求解

为对以式（13）为流场控制方程，式（7）、式（14）为边界条件的地下水井流数学模型进行求解，需对相关等式进行Laplace变换，流场控制方程变换结果如式（15）所示.

\frac{{\mathrm{d}}^{\mathrm{2}}\overline{s_f}}{\mathrm{d}{r}^{\mathrm{2}}}+\frac{n}{r}\frac{\mathrm{d}\overline{s_f}}{\mathrm{d}r}=A{r}^{\mathrm{1}-n}p\overline{s_f}+B{r}^{\mathrm{1}-n}p\overline{s_m},(15)

式中，A=\frac{S_f}{b}\frac{n}{K_q}{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right)}^{n-\mathrm{1}}，B=\frac{S_m}{b}\frac{n}{K_q}{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right)}^{n-\mathrm{1}}.

边界条件变换结果如式（16）、式（17）所示：

\overline{s_f}\left(\mathrm{\infty },p\right)=\mathrm{0},(16)

\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b{r}_w}\right]}^{\left(\mathrm{1}-n\right)}\frac{\mathrm{d}\overline{s_f}\left(r,t\right)}{\mathrm{d}r}{\left|\right.}_{r\to r_w}-\\ \mathrm{\pi }{r}_c^{\mathrm{2}}p\overline{s_f}\left(r_w,p\right)=-\frac{Q}{p}\end{array},(17)

对式（5）进行Laplace变换，有：

p{S}_m\overline{s_m}=C\left(\overline{s_f}-\overline{s_m}\right),(18)

对式（18）进行移项，得到：

\overline{s_m}=C\overline{s_f}/\left(p{S}_m+C\right),(19)

将式（19）代入式（15）得到：

\frac{{\mathrm{d}}^{\mathrm{2}}\overline{s_f}}{\mathrm{d}{r}^{\mathrm{2}}}+\frac{n}{r}\frac{\mathrm{d}\overline{s_f}}{\mathrm{d}r}=\left[Ap+Bp\frac{C}{\left(p{S}_m+C\right)}\right]r^{\mathrm{1}-n}\overline{s_f},(20)

式（20）有通解

\overline{s_f}\left(r,p\right)=r^{\frac{\mathrm{1}-n}{\mathrm{2}}}\left[\begin{array}{l}{M}_{\mathrm{1}}{I}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{\left(p{S}_m+C\right)}}\right)+\\ M_{\mathrm{2}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{\left(p{S}_m+C\right)}}\right)\end{array}\right],(21)

式中，I_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(x\right)和K_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(x\right)分别为改进第一类和第二类贝塞尔\frac{\mathrm{1}-n}{\mathrm{3}-n}阶函数，根据贝塞尔函数的性质，I_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}=\mathrm{\infty }，因此，将边界条件式（16）代入通解式（21）可知M₁ =0，将边界条件式（17）代入通解式（21）可得：

M_{\mathrm{2}}=\frac{Q}{p\left\{\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b{r}_w}\right]}^{\left(\mathrm{1}-\mathrm{n}\right)}{r}_w^{\mathrm{1}-n}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\\ K_{\frac{\mathrm{2}}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\\ +\mathrm{\pi }{r}_c^{\mathrm{2}}p{r}_w^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\end{array}\right\},}(22)

因此，考虑井储效应条件下，抽水过程中裂隙地下水位降深在Laplace空间内的表达式为：

\overline{s_f}\left(r,p\right)=\frac{Q{r}^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{\left(p{S}_m+C\right)}}\right)}{p\left\{\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b{r}_w}\right]}^{\left(\mathrm{1}-\mathrm{n}\right)}{r}_w^{\mathrm{1}-n}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\\ K_{\frac{\mathrm{2}}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\\ +\mathrm{\pi }{r}_c^{\mathrm{2}}p{r}_w^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\end{array}\right\}}.(23)

2 模型讨论

2.1　与其他模型的关系

在不考虑井储效应（即r_w →0，r_c =0）时，双重介质含水层任意位置处裂隙介质水位降深表达式如式（24）所示：

\overline{s_f}\left(r,p\right)=\frac{Q{r}^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)}{p\left\{\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right]}^{\left(\mathrm{1}-\mathrm{n}\right)}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\\ K_{\frac{\mathrm{2}}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\end{array}\right\}}{\left|\right.}_{r_w\to \mathrm{0}},(24)

考虑到贝塞尔函数的性质K_v\left(x\right)\approx \frac{\mathrm{\Gamma }\left(v\right)}{\mathrm{2}}{\left(\frac{x}{\mathrm{2}}\right)}^{-v}，有：

\begin{array}{l}\overline{s_f}\left(r,p\right)=\frac{\mathrm{2}{\left(\frac{Q}{\mathrm{2}\mathrm{\pi }b}\right)}^n{\left(\frac{\mathrm{1}}{\mathrm{3}-n}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}-n}}}{K_{q}p\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\mathrm{\Gamma }\left(\frac{\mathrm{2}}{\mathrm{3}-n}\right)}{r}^{\frac{\mathrm{1}-n}{\mathrm{2}}}\\ K_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap+Bp\frac{C}{p{S}_m+C}}\right)\end{array},(25)

与Wang et al.（2021）提出的不考虑井储效应的双重介质非达西流模型一致.

当C=0或S_m =0时，含水层为单一介质，任意位置处水位降深表达式如式（26）所示：

\overline{s_f}\left(r,p\right)=\frac{Q{r}^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap}\right)}{p\left\{\begin{array}{l}\mathrm{2}\mathrm{\pi }{r}_{w}b{K}_q{\left[\frac{Q}{\mathrm{2}\mathrm{\pi }b{r}_w}\right]}^{\left(\mathrm{1}-\mathrm{n}\right)}{r}_w^{\mathrm{1}-n}\sqrt[]{Ap}{K}_{\frac{\mathrm{2}}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap}\right)+\\ \mathrm{\pi }{r}_c^{\mathrm{2}}p{r}_w^{\frac{\mathrm{1}-n}{\mathrm{2}}}{K}_{\frac{\mathrm{1}-n}{\mathrm{3}-n}}\left(\frac{\mathrm{2}}{\mathrm{3}-n}{r}_w^{\frac{\mathrm{3}-n}{\mathrm{2}}}\sqrt[]{Ap}\right)\end{array}\right\}},(26)

与Wen and Wang（2013）提出的考虑井储效应的均质含水层非达西流模型一致.

取Q=50 m³/d、n=1.5、K_q =0.1（m/d） ⁿ 、r= 0.5 m、r_w =r_c =0.1 m、b=50 m、A=1、B=10、C=5、S_m =0. 1，得到水位降深s⁃时间t曲线如图1所示.

由图1可知，本文和Wang模型均考虑了双重介质假设，抽水过程进行到一定时间后，裂隙介质与孔隙介质间的水头差增大，孔隙介质释水作用显著，从而有裂隙水位下降速率减缓、水位降深‒时间曲线出现拐点的过程；Wen模型则未考虑双重介质条件，曲线整体较为平滑；Wang模型和本文模型在抽水过程后期水位降深逐渐趋于一致，而在抽水过程前期，本文模型由于受到井储效应的影响，水位降深整体比Wang模型大，且在抽水过程中期，拐点的出现较Wang模型提前，说明抽水过程中，抽水井中地下水储量的释放，将对抽水影响范围内地下水位造成扰动.

2.2　参数敏感性分析

2.2.1　非达西系数

在Izbash公式中，当n=1时，地下水为非达西流，而当地下水流速较大时，出现非达西流效应（n>1），该效应对抽水过程中含水层的水位降深产生一定影响.取Q=50 m³/d，K_q =0.1（m/d） ⁿ，r=0.5 m，r_w =r_c =0.1 m，b= 50 m，A=1，B=10，C=5，S_m =0.1，n=1、1.2、1.5、1.7、2，其水位降深s⁃时间t曲线如图2所示.

由图2可知，其他条件相同时，n值越大，非达西效应越明显，地下水流速越大，所表现出的含水岩组释水能力越强.从而，当抽水时间相同时，n越大，裂隙含水层受到的地下水补给量越多，水位降深越小.

2.2.2　孔隙‒裂隙水量交换系数

孔隙‒裂隙水量交换系数（C值）代表其他条件相同时，孔隙‒裂隙介质间水量交换作用的强度，该值体现出双重介质含水层储水介质（孔隙）与导水介质（裂隙）间的地下水联通情况.取Q=50 m³/d，K_q =0.1（m/d） ⁿ，r=r_w =r_c =0.1 m，b=50 m，A=1，B=10，S_m =0.1，n=1.5，C=0、1、10、100、 1 000，其水位降深s⁃时间t曲线如图3所示.

由图3可知，孔隙‒裂隙水量交换系数对水位降深的影响主要表现在抽水过程前期.当C=0时，模型为均质孔隙含水层，不存在孔隙‒裂隙介质间的水量交换，其他条件相同时，水位降深最大；当C≠0时，模型为双重介质含水层，水位降深s⁃时间t曲线表现出双重介质所特有的“S”形曲线特征，且孔隙‒裂隙介质水量交换系数C越大，孔隙‒裂隙介质水量交换越明显，抽水过程前期导水介质（裂隙）受到储水介质（孔隙）的补给作用越强，抽水过程前期水位降深越小.

2.2.3　井储效应

当考虑井储效应时，抽水过程中由抽水井释放的水量由两部分构成.一为含水层本身释水，这部分水量主要来源于含水层中的储水介质（孔隙）；二为井储作用释水，即钻孔中存储的水量随抽水过程中地下水位下降而释水.因此，井储效应的强弱将对抽水过程产生影响.为探究其他条件相同时，不同钻孔半径对抽水井中水位降深的影响，取Q= 50 m³/d，K_q =0.1（m/d） ⁿ，n=1.5，b=20 m，A=1，B=10，C=5，S_m =0.1，r=r_c =r_w =0.1 m、0.3 m、0.5 m、0.7 m、0.9 m，则水位降深s⁃时间t曲线如图4所示.

由图4可知，井储效应对抽水试验全过程均产生影响，具体表现为：井径越小，抽水井处的水位降深越大，而该现象在抽水试验前期最为明显，抽水试验后期，不同井径抽水井地下水水位下降速率趋于一致.对其原因进行分析，发现井径越小，抽水井中地下水的储量越少，抽水试验前期，释放相同水量时抽水井中水位下降幅度越大.抽水试验后期，含水层水位逐渐趋于稳定，释出水量主要来源于含水层释水，抽水井中地下水储量已基本释出，占总抽水流量的比例逐渐减小，井储效应对抽水试验的影响逐渐减弱.

3 案例分析

3.1　研究区概况

研究区位于江西省赣州市寻乌县境内，属中低山丘陵及山间河谷地貌.区内出露第四系人工填土（Q₄^ml）和残坡积层（Q^el+dl）粉质黏土、粗角砾土，下伏基岩以碎屑岩为主，地层岩性由老到新依次为侏罗系上统菖蒲群下亚群上组（J₃ch^1b）凝灰熔岩夹凝灰岩，寒武系下组（∈a）长石石英细砂岩夹绢云母粉砂板岩，震旦系上组（Z_Ⅱ）变质砂岩夹板岩，其中，第四系人工填土和粉质粘土则较为致密，可视作区内隔水顶板，而碎屑岩地层则在成岩过程中发育裂隙‒孔洞，形成主要含水层，因此，区内地下水具有承压性的特点.研究区地处华南褶皱系之武夷隆起带，于加里东旋回和燕山旋回期发生多次构造运动，导致研究区构造发育，断层密布，大型断层带及破碎带受剪切作用的影响，岩性极为破碎，渗透性极强，形成主要的导水介质.因此，本区地下水赋存介质包括成岩裂隙‒孔洞和大型构造裂隙带与剪切破碎带，两者含水性和导水性差异明显，具有双重介质特征.

3.2　抽水试验

为分析研究区双重介质含水层的释水过程，获取含水层水文地质参数，本次研究开展了单孔抽水试验，抽水孔半径0.037 5 m，抽水共历时30 000 s，使用加拿大Solinst公司的Levelogger自动水位计监测抽水孔的水位降深情况，水位计分辨率达 1 mm.根据水文地质钻探结果，抽水孔深520 m，其中，0~39.5 m深度为凝灰岩风化裂隙带，以风化裂隙水为主，含水量极为有限，306~320 m深度为凝灰岩破碎带，为钻探过程中发生涌水的主要层位，而破碎带上部和下部分别为弱风化凝灰岩，分别可视为含水顶板和含水底板.抽水孔代表性岩心如表1所示.

3.3　结果与讨论

3.3.1　抽水试验结果

抽水过程中地下水位降深s⁃时间t曲线如图5所示.由图5可见，抽水过程前期地下水位下降迅速，而当抽水过程进行到500 s后，地下水位下降逐渐放缓，出现“S”型曲线的特征.根据Kazemi et al.（1969）对双重介质含水层水位降深过程的分析，抽水前期，主要为裂隙介质释水，孔隙介质中存储的水量释水程度弱，因此，地下水位下降速度快，与拟合初期降深的Theis曲线基本吻合；抽水后期，孔隙介质中的存储的水量逐渐释放，裂隙介质释放水量占总释水量的比例逐渐降低，在总释水量恒定的情况下，地下水位下降速率逐渐放缓.因此，抽水过程反映出含水层所具备的双重介质特征.

3.3.2　参数反演

使用Stehfest方法（Stehfest， 1970a，1970b）将Laplace空间解转换为实域解，并利用Matlab自带的配线工具箱进行参数率定，可实现对水位降深s⁃时间t曲线的拟合.分别采用本文提出的考虑井储的双重介质非达西流模型、考虑井储效应的双重介质达西流模型、Wang et al.（2021）提出的无井储效应的双重介质非达西流模型、Wen and Wang（2013）提出的考虑井储效应的单一介质非达西流模型和目前广泛使用的Theis模型对野外抽水试验实测数据进行拟合，其结果如表2及图6所示.

由表2及图6可知，采用本文提出的考虑井储效应的双重介质非达西流模型，拟合值与预测值的均方根误差（RMSE）最小，与实测数据吻合程度最高.对比考虑井储效应的双重介质非达西流模型和考虑井储效应的达西流模型，使用达西流模型反演所得的含水层给水度较非达西流模型小、而达西流模型渗透系数较非达西流模型大，且达西流模型拟合曲线水位降深在抽水过程前期、后期偏大，说明采用达西流模型反演所得给水度较实际情况偏小，而渗透系数大、地下水流速快，则更易出现非达西流现象，因此，采用非达西流模型进行参数反演更为合理；对比考虑井储效应的双重介质非达西流模型和不考虑井储效应的双重介质非达西流模型（Wang模型），不考虑井储效应时，抽水过程前期拟合降深较实际降深大，对抽水过程前期抽水孔原有水量大量释放，含水层释水比例占总释水量较小，水位降深变幅较缓的现象刻画能力有限；对比考虑井储的双重介质非达西流模型和考虑井储效应的均质含水层非达西流模型（Wen模型），二者对于抽水过程后期的参数拟合程度较好，而采用双重介质模型，则能够较好地刻画出抽水过程前‒中期裂隙含水层储量释放，孔隙‒裂隙介质水头差增大，孔隙含水层释水逐渐占主导作用这一过程中，裂隙含水层逐渐接受孔隙含水层补给，地下水水位下降受抑制的现象；传统Theis模型仅能够对抽水过程后期，含水层流场逐渐趋于稳定时的抽水孔水位降深情况进行拟合，无法刻画抽水过程前期裂隙介质快速释水，地下水位急剧下降的过程.

4 结论

本文在Wang模型的基础上，提出了考虑井储效应的双重介质非达西流/达西流模型，并将该模型与考虑井储效应的均质含水层非达西流模型、不考虑井储效应的双重介质非达西流模型和Theis模型进行对比，同时，讨论了不同参数对水位降深s⁃时间t曲线的影响，最后，将本文提出的模型运用于工程实际中的单孔抽水试验中，得到如下结论.

（1）当抽水流量较大、地下水流速较快时，抽水井附近存在非达西流效应.非达西流效应越明显，则反映出含水层地下水补给、径流、排泄条件越好，因此，当抽水流量相同时，非达西流效应越强的含水层，往往更容易接受补给，从而在相同时间内，水位降深越小.

（2）双重介质含水层中，孔隙‒裂隙介质地下水的相互作用强弱可通过孔隙‒裂隙介质水量交换系数来定量表征，该系数对抽水过程的影响主要体现在抽水过程前‒中期，因为该时期孔隙介质与裂隙介质水头差较小，孔隙‒裂隙介质水量交换系数越大，则起到导水作用的裂隙介质更易接受起到储水作用的孔隙介质含水层的补给，其他条件相同时，抽水孔水位降深相对较小；抽水过程后期，含水层水位逐渐趋于稳定，双重介质含水层中的地下水逐渐形成新的平衡态，其他条件相同时，不同孔隙‒裂隙介质水量交换系数的双重介质含水层水位降深也逐渐趋于一致.

（3）抽水过程中，井储效应对抽水孔水位降深的影响主要表现于抽水过程前期抽水孔地下水位急剧下降的阶段，上述阶段，由于地下水位变动迅速，抽水孔中的地下水储量也相应明显得到释放，抽水孔地下水储量释放量占总抽水流量的比例较大，抽水过程后期，抽水孔和含水层地下水位逐渐趋于稳定，抽水孔地下水储量释放量占总抽水流量的比例逐渐减小，井储效应对抽水孔水位降深的影响也相应减弱.

（4）本文提出的考虑井储效应的双重介质非达西流模型，能够较好地适用于本文出现的由导水断裂控制的裂隙含水层，特别是在仅允许开展单孔抽水的情况时，能够较好地刻画出抽水过程初期，由于抽水孔释出的水量占总抽水流量比例较大，而导致抽水孔中地下水位快速下降的现象，对工程实践有一定参考价值.

（5）本文提出的考虑井储效应的双重介质单孔抽水非稳定流解析解模型将孔隙介质和裂隙介质概化为两组渗透系数、储水系数不同，而彼此之间存在水量交换的含水介质，无法对裂隙发育程度、发育位置等因素做更为详细的刻画，建议在后续研究中考虑上述因素，对模型加以改进与优化.

附图见https：//doi.org/10.3799/dqkx.2022.053.

参考文献

[1] Barenblatt, G. I., Zheltov, Y. P., Kochina, I. N., 1960. Basic Concepts in the Theory of Seepage of Homogeneous Liquids in Fissured Rocks. Journal of Applied Mathematics, 24(5): 1286-1303.

[2] Bordier, C., Zimmer, D., 2000. Drainage Equations and Non⁃Darcian Modelling in Coarse Porous Media or Geosynthetic Materials. Journal of Hydrology, 228(3/4): 174-187. https://doi.org/10.1016/S0022⁃1694(00)00151⁃7

[3] De Smedt, F., 2011. Analytical Solution for Constant⁃Rate Pumping Test in Fissured Porous Media with Double⁃Porosity Behaviour. Transport in Porous Media, 88(3): 479-489. https://doi.org/10.1007/s11242⁃011⁃9750⁃9

[4] Fenske, P. R., 1977. Radial Flow with Discharging⁃Well and Observation⁃Well Storage. Journal of Hydrology, 32(1-2): 87-96. https://doi.org/10.1016/0022⁃1694(77)90120⁃2

[5] Forchheimer, P. H., 1901. Wasserbewegun Durch Boden. Zeitschrift des Vereines Deutscher Ingenieure, 50: 1781-1788. https://doi.org/10.1016/0022⁃1694(77)90120⁃2

[6] He, L. Q., Huang, K., Wan, J. W., et al., 2015. High⁃Resolution Identification of Hydrogeological Parameters in Dual⁃Aqueous Medium. Geological Science and Technology Information, 34(1): 216-220 (in Chinese with English abstract).

[7] Huang, K., Sun, R. L., Yuan, S. Q., et al., 2022. Effect of Number of Pumping Tests and Prior Information on Hydraulic Conductivity Estimation of Three⁃Dimensional Heterogeneous Aquifer. Earth Science, 47(2): 689-699 (in Chinese with English abstract).

[8] Izbash, S., 1931. O Filtracii Krupnozernistom Materiale. Izv Nauchnoissled Inst. Gidrotechniki (NIIG), Leningrad.

[9] Kazemi, H., Seth, M. S., Thomas, G. W., 1969. The Interpretation of Interference Tests in Naturally Fractured Reservoirs with Uniform Fracture Distribution. Society of Petroleum Engineers Journal, 9(4): 463-472.

[10] Li, C., Wang, N. T., Hu, C., et al., 2020. Hydrogeological Parameter Determination of the Confined Aquifer under Confined⁃Unconfined Conditions by Means of Multiple⁃Well Pumping Test. Safety and Environmental Engineering, 27(4): 1-7 (in Chinese with English abstract).

[11] Li, N., 2020. The Formation Mechanism and Hydrogeological Parameters of Karst Thermal Water: A Case Study of Nanyang Thermal Spring in Xingshan County of Hubei Province, Xiangxi River Basin (Dissertation). China University of Geosciences, Wuhan (in Chinese with English abstract).

[12] Li, X., Su, S. L., Wen, Z., et al., 2022. Numerical Analysis of Estimating Groundwater Velocity through Single⁃Well Push⁃Pull Test. Earth Science, 47(2): 633-641 (in Chinese with English abstract).

[13] Li, Z. X., Wan, J. W., Zhan, H. B., et al., 2019. Particle Size Distribution on Forchheimer Flow and Transition of Flow Regimes in Porous Media. Journal of Hydrology, 574: 1-11. https://doi.org/10.1016/j.jhydrol.2019.04.026

[14] Moench, A. F., 1984. Double⁃Porosity Models for a Fissured Groundwater Reservoir with Fracture Skin. Water Resources Research, 20(7): 831-846. https://doi.org/10.1029/WR020i007p00831

[15] Stehfest, H., 1970a. Algorithm 368: Numerical Inversion of Laplace Transforms. Communications of the ACM, 13(1): 47-49. https://doi.org/10.1145/361953.361969

[16] Stehfest, H., 1970b. Remark on Algorithm 368: Numerical Inversion of Laplace Transforms. Communications of the ACM, 13(10): 624. https://doi.org/10.1145/355598.362787

[17] Wan, J. W., Huang, K., Chen, C. X., 2013. Reassessing Darcy’ Law on Water Flow in Porous Media. Earth Science, 38(6): 1327-1330 (in Chinese with English abstract).

[18] Wang, Y. L., Zhan, H. B., Huang, K., et al., 2021. Identification of Non⁃Darcian Flow Effect in Double⁃Porosity Fractured Aquifer Based on Multi⁃Well Pumping Test. Journal of Hydrology, 600: 126541. https://doi.org/10.1016/j.jhydrol.2021.126541

[19] Warren, J. E., Root, P. J., 1963. The Behavior of Naturally Fractured Reservoirs. Society of Petroleum Engineers Journal, 3(3): 245-255. https://doi.org/10.2118/426⁃pa

[20] Wen, Z., Huang, G. H., Zhan, H. B., et al., 2008. Two⁃Region Non⁃Darcian Flow toward a Well in a Confined Aquifer. Advances in Water Resources, 31(5): 818-827. https://doi.org/10.1016/j.advwatres.2008.01.014

[21] Wen, Z., Wang, Q. R., 2013. Approximate Analytical and Numerical Solutions for Radial Non⁃Darcian Flow to a Well in a Leaky Aquifer with Wellbore Storage and Skin Effect. International Journal for Numerical and Analytical Methods in Geomechanics, 37(11): 1453-1469. https://doi.org/10.1002/nag.2091

基金资助

国家自然科学基金项目(42172281)