越流含水层系统中渗透系数随深度衰减的抽水试验解析模型

Analytical Model of Pumping Tests with Depth-Dependent Hydraulic Conductivity in Leakage Aquifer System

摘要

传统的抽水试验模型假定含水层的渗透系数是恒定的，然而，由于地质作用等因素的影响，含水层渗透系数存在随深度衰减的现象.建立考虑越流、井储效应和随深度衰减的渗透系数的抽水试验数学模型，并推导了半解析解，其中，假设渗透系数随深度呈指数衰减.系统分析随深度衰减的渗透系数对抽水试验结果的影响.结果表明：含水层渗透系数随深度衰减系数（A）越大，井筒内地下水的降深越大；当含水层渗透系数随深度衰减时，滤水管的位置对抽水试验结果的影响比较显著，滤水管在含水层顶部时的井筒内降深比滤水管在含水层底部时的井筒水位降深要小；传统的抽水试验模型采用井筒水位降深数据反演的常数渗透系数与随深度衰减的渗透系数平均值近似.

关键词

抽水试验 / 数值模拟 / 渗透系数 / 越流 / 水文地质.

Key words

pumping test / numerical modeling / hydraulic conductivity / leakage / hydrogeology

0 引言

含水层水动力学参数对水资源评价、环境污染预测，以及地热开发利用具有重要的意义（万力等，2010；陈崇希，2020；陈崇希和唐仲华，2021；罗颖等，2023）.抽水试验是获取野外含水层水动力学参数最广泛的方法之一.1935年，Theis（1935）提出了第一个非稳定流抽水试验模型，然而该模型包含许多假设条件.随后，许多学者对Theis模型进行了拓展，例如，考虑表皮效应的抽水试验模型（Novakowski，1989；Yang and Yeh，2005； Tsai and Yeh，2012； Lin et al.，2016；Chen et al.，2020）；考虑变流量的抽水试验模型（Lai and Su，1974； Mathias and Todman，2010；Chen et al.，2020）；考虑非达西流的抽水试验模型（Wen et al.，2008； Mathias and Todman，2010； Wen et al.，2011； 文章等，2011；Wang et al.，2021；魏世毅等，2022）；考虑非完整井的抽水试验模型（Hantush，1960， 1964； Cassiani et al.，1999； Yang and Yeh，2005；Yang et al.，2006；周振方，2018）等.

相比沙质越流含水层系统，许多场地研究已经发现淤泥质和粘土质含水层往往可压缩性更大，厚度更厚，含水层的渗透系数在含水层应力下呈现出随深度衰减的现象，例如Achtziger-Zupančič et al.（2017）和Chen et al.（2018）收集的很多场地数据发现含水层渗透系数随深度衰减.最近，Zhuang et al.（2022）的研究发现：忽略弱透水层中的渗透系数随深度衰减时，会导致计算的越流量结果偏低.然而，现有的抽水试验模型假定含水层是均质的，即含水层渗透系数是常数.事实上，大量的野外试验表明：含水层的渗透系数呈现随深度衰减的现象，而且满足如下式（Louis，1972； Zhang and Franklin，1993； Saar and Manga，2004；Wang et al.，2009）：

K \left(z\right) = K_{\mathrm{0}} \mathrm{e} \mathrm{x} \mathrm{p} \left(-Az\right),(1)

式（1）中， K表示渗透系数； K_{\mathrm{0}}表示含水层顶部位置的渗透系数； A表示衰减系数； \mathrm{e} \mathrm{x} \mathrm{p} \left(\bullet \right)表示指数函数； z表示垂向距离.Louis （1972）基于大量的试验得出 A的范围在0.003 4~0.270 0 {\mathrm{m}}^{-\mathrm{1}}，Zhang and Franklin （1993）反演了加拿大某页岩含水层中的 A为0.01 {\mathrm{m}}^{-\mathrm{1}}，Sakata and Ikeda （2013）在地表30 m以内的砾石含水层中获取的 A为0.11 {\mathrm{m}}^{-\mathrm{1}}，Wang et al.（2009）获取了中国西南地区玄武岩中的 A为0.022 {\mathrm{m}}^{-\mathrm{1}}.

由于现有的抽水试验模型无法反演深度衰减的渗透系数，为此，本研究将构建考虑越流、井储效应和随深度衰减的渗透系数的抽水试验数学模型，并采用拉普拉斯变换、Goldstein-Weber变换和格林函数法推导解析解，为反求含水层水动力学参数提供了理论模型.

1 考虑越流和变渗透系数抽水试验模型

1.1　数学模型

考虑越流和变渗透系数抽水试验模型如图1所示，坐标原点位于含水层顶部井筒中央位置.基于前人模型（Wen et al.，2011； Sakata and Ikeda，2013），建立越流条件下变渗透系数的抽水试验数学模型：

\frac{\mathrm{1}}{r} \frac{\partial }{\partial r} \left[rK\left(z\right)\frac{\partial h}{\partial r}\right] + \frac{\partial }{\partial z} \left[K\left(z\right)\frac{\partial h}{\partial z}\right] = S_{\mathrm{s}} \frac{\partial h}{\partial t}，r \ge r_{\mathrm{w}}，,(2)

{\left.h\left(r,z,t\right)\right|}_{t=\mathrm{0}} = h_{\mathrm{0}}，r \ge r_{\mathrm{w}}，,(3)

{\left.h\left(r,z,t\right)\right|}_{r\to \mathrm{\infty }} = h_{\mathrm{0}}，,(4)

\mathrm{2} \pi \left(l-d\right) r_{\mathrm{w}} K \left(z\right) {\left.\frac{\partial h}{\partial r}\right|}_{r=r_{\mathrm{w}}} = \left[-Q+\pi r_{\mathrm{w}}^{\mathrm{2}}\frac{\mathrm{d}h}{\mathrm{d}t}\right] H_{\mathrm{v}} \left(z\right)，,(5)

H_{\mathrm{v}} \left(z\right) = H \left(z-d\right) - H \left(z-l\right)，,(6)

K \left(z\right) {\left.\frac{\partial h}{\partial z}\right|}_{z=\mathrm{0}} = K_a \frac{h_{\mathrm{0}}-h}{m}，r \ge r_{\mathrm{w}},(7a)

{\left.\frac{\partial h}{\partial z}\right|}_{z=B} = \mathrm{0}，r \ge r_{\mathrm{w}},(7b)

式（2~7）中， h是承压含水层水头［ L］； h_{\mathrm{0}}是潜水含水层水头［ L］； B是承压含水层的厚度［ L］； S_s是弹性贮水或释水率［ \mathrm{1} / L］； r表示径向距离［ L］； r_w表示井半径［ L］； t表示时间［ T］； K_a是弱透水层的渗透系数［ L / T］， m表示弱透水层的厚度［ L］； Q表示抽水流量［ L^{\mathrm{3}} T］； l和 d分别表示滤水管的顶部和底部距离坐标原点的距离［ L］； H \left(\bullet \right)表示海维塞德函数.为了求解上述数学模型的解析解，式（5）中的进水段的渗透系数 K \left(z\right)假定为常数 K_{\mathrm{0}}.

1.2　解析解的推导

公式（2~7）构成了考虑井储效应和越流作用下的变渗透系数抽水试验数学模型，本研究将采用拉普拉斯变换法、格林函数法和Goldstein-Weber变换法推导拉普拉斯空间下的解析解（Kabala and Cassiani，1997）.假设 s表示降深［ L］，且 s = h_{\mathrm{0}} - h，将式（1）带入式（2）得：

\frac{\mathrm{1}}{r} \frac{\partial }{\partial r} \left[r\frac{\partial s}{\partial r}\right] + \frac{{\partial }^{\mathrm{2}}s}{\partial z^{\mathrm{2}}} - A \frac{\partial s}{\partial z} = \frac{S_{\mathrm{s}}}{K_{\mathrm{0}}} \mathrm{e} \mathrm{x} \mathrm{p} \left(Az\right) \frac{\partial s}{\partial t}，r \ge r_{\mathrm{w}}， .(8)

对式（3~8）进行拉普拉斯变换和Goldstein-Weber变换，得：

\frac{{\partial }^{\mathrm{2}}\check{\overline{s}}}{\partial z^{\mathrm{2}}} - A \frac{\partial \check{\overline{s}}}{\partial z} - \left[\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(Az\right)+{\alpha }^{\mathrm{2}}-\frac{r_{\mathrm{w}}p{H}_{\mathrm{v}}\left(z\right)}{\pi \alpha \left(l-d\right)K_{\mathrm{0}}}\right] \check{\overline{s}} = \frac{Q{H}_{\mathrm{v}}\left(z\right)}{{\pi }^{\mathrm{2}}\alpha \left(l-d\right)r_{\mathrm{w}}{K}_{\mathrm{0}}p},(9)

{\left.\frac{\partial \check{\overline{s}}}{\partial z}\right|}_{z=\mathrm{0}} - \frac{K_a}{m{K}_{\mathrm{0}}} \check{\overline{s}} = \mathrm{0},(10)

{\left.\frac{\partial \check{\overline{s}}}{\partial z}\right|}_{z=B} = \mathrm{0},(11)

式（9~11）中， p和 \alpha分别表示拉普拉斯变换参数和Goldstein-Weber变换参数，上标“ \vee”和“-”分别表示Goldstein-Weber域和拉普拉斯域.

类似于Wang et al.（2018），式（9~11）将采用格林函数法进行求解.式（9）是非齐次变系数偏微分方程，其拉普拉斯域下的通解为：

\overline{s} = {\int }_{\mathrm{0}}^{\mathrm{\infty }} \check{\overline{s}}\left(\alpha \right) \frac{J_{\mathrm{0}}\left(\alpha r\right)\left[F\alpha Y_{\mathrm{1}}\left(\alpha r_{\mathrm{w}}\right)+Y_{\mathrm{0}}\left(\alpha r_{\mathrm{w}}\right)\right]-Y_{\mathrm{0}}\left(\alpha r\right)\left[F\alpha J_{\mathrm{1}}\left(\alpha r_{\mathrm{w}}\right)+J_{\mathrm{0}}\left(\alpha r_{\mathrm{w}}\right)\right]}{{\left[F\alpha J_{\mathrm{1}}\left(\alpha r_{\mathrm{w}}\right)+J_{\mathrm{0}}\left(\alpha r_{\mathrm{w}}\right)\right]}^{\mathrm{2}}+{\left[F\alpha Y_{\mathrm{1}}\left(\alpha r_{\mathrm{w}}\right)+Y_{\mathrm{0}}\left(\alpha r_{\mathrm{w}}\right)\right]}^{\mathrm{2}}} \mathrm{\alpha } \mathrm{d} \mathrm{\alpha },(12)

\check{\overline{s}} = \frac{\mathrm{2}Q}{{\pi }^{\mathrm{2}}{r}_{\mathrm{w}}^{\mathrm{2}}{p}^{\mathrm{2}}} {\int }_{\mathrm{0}}^B g\left(z,p,\alpha ,\lambda \right) \mathrm{d} \lambda,(13)

g \left(z,p,\alpha ,\lambda \right) = \left\{\begin{array}{c}{g}_{\mathrm{1}}=N_{\mathrm{1}}{F}_{\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{Az}{\mathrm{2}}\right)I_{-v}\left[\frac{\mathrm{2}}{A}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(Az\right)}\right]+N_{\mathrm{2}}{F}_{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{Az}{\mathrm{2}}\right)I_v\left[\frac{\mathrm{2}}{A}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(Az\right)}\right],\mathrm{0}\le z\le \lambda \\ g_{\mathrm{2}}=N_{\mathrm{3}}{F}_{\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{Az}{\mathrm{2}}\right)I_{-v}\left[\frac{\mathrm{2}}{A}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(Az\right)}\right]+N_{\mathrm{4}}{F}_{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{Az}{\mathrm{2}}\right)I_v\left[\frac{\mathrm{2}}{A}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(Az\right)}\right],\lambda \le z\le B\end{array}\right.,(14)

式（12~14）中， \lambda是积分变量； f \left(\lambda \right) = \frac{Q{H}_{\mathrm{v}}\left(\lambda \right)}{{\pi }^{\mathrm{2}}\alpha \left(l-d\right)K_{\mathrm{0}}{r}_{\mathrm{w}}p}； g \left(z,p,\alpha ,\lambda \right)表示格林函数； F = \frac{\mathrm{2}\left(l-d\right)K_{\mathrm{0}}}{r_{\mathrm{w}}p{H}_{\mathrm{v}}\left(z\right)}； v = \frac{\sqrt[]{A^{\mathrm{2}}+\mathrm{4}{\alpha }^{\mathrm{2}}}}{A}； F_{\mathrm{1}} = \Gamma \left(\mathrm{1}-v\right) i^{-v}； F_{\mathrm{2}} = \Gamma \left(\mathrm{1}+z\right) i^v； i是虚数单位； \Gamma \left(\bullet \right)是伽马函数； I_v \left(\bullet \right)是 v阶第一类修正贝塞尔函数； N_{\mathrm{1}}， N_{\mathrm{2}}， N_{\mathrm{3}}和 N_{\mathrm{4}}是待求解的系数.

由于式（13）中包含4个待求的系数，仅由式（10）和式（11）难以求解，因此，需要再补充两个额外的条件.根据格林函数的性质，补充的两个方程为：

g_{\mathrm{1}} \left(z=\lambda ,p,\alpha ,\lambda \right) = g_{\mathrm{2}} \left(z=\lambda ,p,\alpha ,\lambda \right)，(15)

{\left.\frac{\mathrm{d}{g}_{\mathrm{2}}}{\mathrm{d}z}\right|}_{z=\lambda } - {\left.\frac{\mathrm{d}{g}_{\mathrm{1}}}{\mathrm{d}z}\right|}_{z=\lambda } = - \mathrm{1} .(16)

将式（14）带入式（10~11）和式（15~16），有：

N_{\mathrm{1}} H_{\mathrm{1}} + N_{\mathrm{2}} H_{\mathrm{2}} = \mathrm{0}，(17)

N_{\mathrm{1}} H_{\mathrm{3}} + N_{\mathrm{2}} H_{\mathrm{4}} = N_{\mathrm{3}} H_{\mathrm{3}} + N_{\mathrm{4}} H_{\mathrm{4}}，(18)

N_{\mathrm{3}} H_{\mathrm{5}} + N_{\mathrm{4}} H_{\mathrm{6}} - N_{\mathrm{1}} H_{\mathrm{5}} - N_{\mathrm{2}} H_{\mathrm{6}} = - \mathrm{1}，(19)

N_{\mathrm{3}} H_{\mathrm{7}} + N_{\mathrm{4}} H_{\mathrm{8}} = \mathrm{0} .(20)

求解式（17~20）可得：

N_{\mathrm{2}} = \frac{H_{\mathrm{1}}}{W{H}_{\mathrm{3}}+H_{\mathrm{4}}} N_{\mathrm{3}} + \frac{H_{\mathrm{2}}}{W{H}_{\mathrm{3}}+H_{\mathrm{4}}} N_{\mathrm{4}},(21b)

N_{\mathrm{3}} = \frac{H_{\mathrm{6}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)-H_{\mathrm{2}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)}{H_{\mathrm{1}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)-H_{\mathrm{3}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)} N_{\mathrm{4}} + \frac{W{H}_{\mathrm{3}}+H_{\mathrm{4}}}{H_{\mathrm{1}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)-H_{\mathrm{3}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)},(21c)

N_{\mathrm{4}} = \frac{-\frac{W{H}_{\mathrm{3}}+H_{\mathrm{4}}}{H_{\mathrm{1}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)-H_{\mathrm{3}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)}{H}_{\mathrm{7}}}{\frac{H_{\mathrm{6}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)-H_{\mathrm{2}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)}{H_{\mathrm{1}}\left(W{H}_{\mathrm{5}}+H_{\mathrm{6}}\right)-H_{\mathrm{3}}\left(W{H}_{\mathrm{3}}+H_{\mathrm{4}}\right)}{H}_{\mathrm{7}}+H_{\mathrm{8}}},(21d)

式（21）中，

W = - \frac{H_{\mathrm{2}}}{H_{\mathrm{1}}}，(22)

\begin{array}{l}{H}_{\mathrm{1}}=\frac{A}{\mathrm{2}}{F}_{\mathrm{1}}{I}_{-v}\left[{\phi }_{\mathrm{1}}\right]+\frac{F_{\mathrm{1}}\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\left\{I_{-\mathrm{1}-v}\left[{\phi }_{\mathrm{1}}\right]+I_{\mathrm{1}-v}\left[{\phi }_{\mathrm{1}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)}}-\\ \frac{K_{\mathrm{a}}}{m{K}_{\mathrm{0}}}{F}_{\mathrm{1}}{I}_{-v}\left[{\phi }_{\mathrm{1}}\right]\end{array}，(23)

H_{\mathrm{2}} = \frac{A}{\mathrm{2}} F_{\mathrm{2}} I_v \left[{\phi }_{\mathrm{1}}\right] + \frac{F_{\mathrm{2}}\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\left\{I_{-\mathrm{1}+v}\left[{\phi }_{\mathrm{1}}\right]+I_{\mathrm{1}+v}\left[{\phi }_{\mathrm{1}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)}} - \frac{K_{\mathrm{a}}}{m{K}_{\mathrm{0}}} F_{\mathrm{2}} I_{-v} \left[{\phi }_{\mathrm{2}}\right]，(24)

H_{\mathrm{3}} = F_{\mathrm{1}} e x p \left(\frac{A\lambda }{\mathrm{2}}\right) I_{-v} \left[{\phi }_{\mathrm{2}}\right]，(25)

H_{\mathrm{4}} = F_{\mathrm{2}} e x p \left(\frac{A\lambda }{\mathrm{2}}\right) I_v \left[{\phi }_{\mathrm{2}}\right]，(26)

H_{\mathrm{5}} = \frac{A}{\mathrm{2}} F_{\mathrm{1}} \mathrm{e} \mathrm{x} \mathrm{p} \left(\frac{A\lambda }{\mathrm{2}}\right) I_{-v} \left[{\phi }_{\mathrm{2}}\right] + \frac{F_{\mathrm{1}}\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{3}A\lambda }{\mathrm{2}}\right)\left\{I_{-\mathrm{1}-v}\left[{\phi }_{\mathrm{2}}\right]+I_{\mathrm{1}-v}\left[{\phi }_{\mathrm{2}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(A\lambda \right)}}，(27)

H_{\mathrm{6}} = \frac{A}{\mathrm{2}} F_{\mathrm{2}} \mathrm{e} \mathrm{x} \mathrm{p} \left(\frac{A\lambda }{\mathrm{2}}\right) I_v \left[{\phi }_{\mathrm{2}}\right] + \frac{F_{\mathrm{2}}\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{3}A\lambda }{\mathrm{2}}\right)\left\{I_{-\mathrm{1}+v}\left[{\phi }_{\mathrm{2}}\right]+I_{\mathrm{1}+v}\left[{\phi }_{\mathrm{2}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(A\lambda \right)}}，(28)

H_{\mathrm{7}} = \frac{A}{\mathrm{2}} F_{\mathrm{1}} \mathrm{e} \mathrm{x} \mathrm{p} \left(\frac{AB}{\mathrm{2}}\right) I_{-v} \left[{\phi }_{\mathrm{3}}\right] + \frac{F_{\mathrm{1}}\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{3}AB}{\mathrm{2}}\right)\left\{I_{-\mathrm{1}-v}\left[{\phi }_{\mathrm{3}}\right]+I_{\mathrm{1}-v}\left[{\phi }_{\mathrm{3}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(AB\right)}}，(29)

H_{\mathrm{8}} = \frac{A}{\mathrm{2}} F_{\mathrm{2}} \mathrm{e} \mathrm{x} \mathrm{p} \left(\frac{AB}{\mathrm{2}}\right) I_v \left[{\phi }_{\mathrm{3}}\right] + \frac{F_{\mathrm{2}}\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{3}AB}{\mathrm{2}}\right)\left\{I_{-\mathrm{1}+v}\left[{\phi }_{\mathrm{3}}\right]+I_{\mathrm{1}+v}\left[{\phi }_{\mathrm{3}}\right]\right\}}{\mathrm{2}\sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(AB\right)}},(30)

{\phi }_{\mathrm{1}} = \frac{\mathrm{2}}{A} \sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)},(31)

{\phi }_{\mathrm{2}} = \frac{\mathrm{2}}{A} \sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(A\lambda \right)},(32)

{\phi }_{\mathrm{3}} = \frac{\mathrm{2}}{A} \sqrt[]{\left(\frac{S_{\mathrm{s}}p}{K_{\mathrm{0}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(AB\right)} .(33)

解析解（12~13）是拉普拉斯空间下的解析解，由于包含复杂的贝塞尔函数，很难通过解析逆变换的方法获取实数空间下的解析解.Wang and Zhan （2015）指出Stehfest数值逆变换的方法在求解径向流问题上具有较高的精度，因此，本研究将采用Stehfest数值逆变换获取该模型在时间空间下的结果.

1.3　研究解析解特点

本文建立的考虑越流、井储和渗透系数随深度衰减的抽水试验解析模型是许多前人模型的进一步拓展，正如引言所述，大量的研究表明很多情况下存在含水层的渗透系数呈现随深度衰减的现象，本文提出的模型旨在依托抽水试验精确获取该这类地质含水层系统中水力参数，从而为地下水资源的开发利用与保护提供理论支撑.该模型在特殊条件下能够还原为先前的解析模型：当 A = \mathrm{0}时，本文模型还原为Hantush （1964）的解；当 A = \mathrm{0}且消除井储效应时，本文建立的模型还原为现有不考虑越流和井储效应的均质含水层中的抽水试验模型，例如Yang and Yeh （2005）和Yang et al.（2006）；当 A = \mathrm{0}， K_{\mathrm{a}} \to \mathrm{0}， d = \mathrm{0}和 l = B时，模型变为Papadopulos and Cooper （1967）的完整井抽水的情况；当 A = \mathrm{0}， K_{\mathrm{a}} \to \mathrm{0}， d = \mathrm{0}， l = B和 r_{\mathrm{w}} \to \mathrm{0}，本文的模型还原为Theis（1935）的模型.

1.4　模型的验证

本研究以Yang et al.（2006）模型为标准，验证解析解（12~13）的推导过程以及Stehfest数值逆变换的精度，选取3个不同位置（ r = r_{\mathrm{w}} = \mathrm{0.15} \mathrm{m}， r = \mathrm{0.25} \mathrm{m}和 r = \mathrm{0.35} \mathrm{m}）作为观测点进行对比.解析解模型中采用的参数除了 A = \mathrm{0}和 K_{\mathrm{a}} = \mathrm{0}外，其他参数值见表1、表1中的参数值源于Zheng et al.（2019）的野外试验数据.对比结果如图2所示，本研究的解与前人模型结果拟合效果较好，表明新模型的推导过程无误，数值逆变换的精度可以接受.

2 结果与讨论

2.1　衰减系数对流场的影响

衰减系数（ A）是描述渗透系数随深度衰减的主要参数.当 A = \mathrm{0}，渗透系数为常数，即含水层变为均质含水层.为了研究 A对流场的影响，参考Louis （1972）和Wang et al.（2009），选择3组不同的 A值（0、0.011 {\mathrm{m}}^{-\mathrm{1}}和0.022 {\mathrm{m}}^{-\mathrm{1}}），当含水层的厚度为8 m时，3组不同的 A对应的渗透系数变化量（ ∆ K）分别为：0、0.337 0和0.645 5 \mathrm{m} / \mathrm{d}，其他的参数详见表1.

图3和图4分别表示不同 A值下井筒内水位降深情况和不同 A值下含水层中水位降深的分布.对比发现，当 A越大，井筒内的降深值越大；一个更大的 A导致一个更大的影响范围.图3和图4的结果表明：当含水层渗透系数随深度衰减越快，抽水试验过程中井筒内的降深越大，导致含水层中存在一个更大的影响范围（降落漏斗更大），这是因为含水层渗透系数随深度衰减越快，含水层渗透系数变化量越大，整个含水层厚度上的平均渗透系数越小，导致抽水过程中井筒内降深越大.

2.2　越流补给对抽水试验的影响

由于抽水的作用，导致目标含水层水头降低，相邻含水层中的地下水在水力梯度的作用下，通过弱透水层进入目标含水层，这个现象被称为越流补给.类似于Hantush （1964），本研究采用越流系数刻画越流补给量，这里采用 R表示越流系数［ L T］，其定义为 B m / K_{\mathrm{a}}，为了研究越流补给对抽水试验的影响，选择3组不同的 R值，分别是 R = \mathrm{1} \mathrm{000} \mathrm{m} \mathrm{d}， R = \mathrm{2} \mathrm{000} \mathrm{m} \mathrm{d}和 R = \mathrm{3} \mathrm{000} \mathrm{m} \mathrm{d}.这些参数都来自于Zheng et al.（2019）场地研究.其他参数与表1一致. R值越小，表明越流量越大，当 R = \mathrm{0}，表示没有越流量，含水层的顶板是不透水的.

图5表示不同越流系数 R下的井筒内水位降深，为了对比研究，将没有越流条件下的结果也放在该图中.图5表明，越流系数 R对抽水试验的结果影响显著，因此，在抽水试验过程中越流补给不能忽视.同时， \mathrm{越流} \mathrm{系数} R越小，井筒内抽水试验降深越小.根据 \mathrm{越流} \mathrm{系数} R的定义，越流系数 R与含水层和弱透水层的厚度的乘积成正比，而与弱透水层的渗透系数成反比，对于含水层和弱透水层厚度固定的情况下，弱透水层的渗透系数越大， \mathrm{越流} \mathrm{系数} R越小，进而导致更多的水能通过越流补给含水层，最终表现为抽水试验过程中井筒内水位降深越小.

2.3　井储效应对结果的影响

井储效应是影响抽水试验的主要因素之一.在抽水试验初期，抽出的水主要来自井筒内原有的地下水.根据式（5）可知，影响井储效应的因素主要是井内地下水的体积，本研究采用3个不同的井半径来研究井储效应对抽水试验结果的影响： r_{\mathrm{w}} = \mathrm{0.1} \mathrm{m}， r_{\mathrm{w}} = \mathrm{0.15} \mathrm{m}和 r_{\mathrm{w}} = \mathrm{0.2} \mathrm{m}.图6表示不同井半径下的井筒内降深随时间的变化，其他参数与表1一致.

由图6可知，井储效应对抽水试验结果影响显著，具体表现为：井径越小，抽水井内的水位降深越大.这是由于井径越小，抽水井中地下水的储量越少，井储作用所释放水量越少，抽水井中所抽水量主要来源于含水层释水，导致井筒内水位降深越大.

2.4　滤水管位置对结果的影响

当含水层比较厚时，井筒往往不是完全穿透含水层，即非完整井.对于非均质介质的含水层，滤水管的位置可能对抽水试验结果造成较大的影响，因此，本文设置3组不同的滤水管位置进行研究： d = \mathrm{0} \mathrm{m}， l = \mathrm{4} \mathrm{m}； d = \mathrm{2} \mathrm{m}， l = \mathrm{6} \mathrm{m}和 d = \mathrm{4} \mathrm{m}， l = \mathrm{8} \mathrm{m}； 3种情况分别代表的是滤水管位于含水层的顶部、中央和底部位置，除了 R = \mathrm{1} \mathrm{000} \mathrm{m} \mathrm{d}外，其他参数采用表1中的值.图7a和图7b分别表示均质含水层（ A = \mathrm{0} {\mathrm{m}}^{-\mathrm{1}}）和渗透系数随深度衰减的含水层（ A = \mathrm{0.011} {\mathrm{m}}^{-\mathrm{1}}）中抽水井井筒内降深随时间的变化图.

由图7a可知，对于均质含水层，当滤水管的长度相同时，滤水管位于含水层顶部和底部对抽水试验结果没有影响，而滤水管在含水层中央时，抽水井井筒内的降深比滤水管在 d = \mathrm{0} \mathrm{m}， l = \mathrm{4} \mathrm{m}和 d = \mathrm{4} \mathrm{m}， l = \mathrm{8} \mathrm{m}情况下要小，这主要是由于滤水管在含水层中央位置时，含水层顶部和底部的水能够流到井筒中.与图7a不同的是，对于渗透系数随深度衰减的含水层，滤水管在 d = \mathrm{0} \mathrm{m}， l = \mathrm{4} \mathrm{m}和 d = \mathrm{4} \mathrm{m}， l = \mathrm{8} \mathrm{m}两种情况下的井筒内降深不一致，滤水管在 d = \mathrm{4} \mathrm{m}， l = \mathrm{8} \mathrm{m}具有更大的降深，如图7b所示.这主要是因为含水层底部的渗透系数比顶部要小.图7a和7b表明，在抽水试验过程中，滤水管的位置对抽水试验结果影响比较显著，尤其是针对渗透系数随深度衰减的含水层中的抽水试验.

3 理论抽水试验研究

目前，大部分的恒定渗透系数抽水试验模型被耦合到Aquifer Test 软件中，被广泛用于反求含水层水动力学参数.为了研究恒定渗透系数模型的精度，首先，基于给定的含水层参数，采用本研究模型产生一系列理论抽水试验数据，即井筒水位降深.其次，采用传统的抽水试验模型（Aquifer Test软件）来反演含水层参数.通过对比已知的参数与反演的参数，研究传统模型的精度.基于上述研究思路，本研究首先选取4个不同的衰减系数（ A = \mathrm{0.01}，0.10，0.20和0.30 \mathrm{1} / \mathrm{m}）用于模型计算，为了方便对比，本研究也计算出不同 A值计算的含水层渗透系数的平均值：

K_{\mathrm{a}\mathrm{v}\mathrm{g}} = \frac{\mathrm{1}}{B} {\int }_{\mathrm{0}}^B K\left(z\right) \mathrm{d} z，(34)

其他的含水层参数值采用表1中给定的，采用本研究模型计算的不同 A下的井水位降深如图8所示.

通过Aquifer Test软件进行渗透系数反演的结果如表2所示，其中 K_{\mathrm{e}\mathrm{s}\mathrm{t}}表示采用Aquifer Test软件反演的渗透系数， E表示反演的渗透系数 K_{\mathrm{e}\mathrm{s}\mathrm{t}}和含水层平均渗透系数 K_{\mathrm{a}\mathrm{v}\mathrm{g}}之间的相对误差（相对误差计算公式： \left|\frac{K_{\mathrm{e}\mathrm{s}\mathrm{t}}-K_{\mathrm{a}\mathrm{v}\mathrm{g}}}{K_{\mathrm{a}\mathrm{v}\mathrm{g}}}\right| \times \mathrm{100} \mathrm{\%}）.由表2可知，随着 A值增加，常数渗透系数模型采用井筒水位降深数据反演的渗透系数呈现出减小的趋势，并且反演的渗透系数值与含水层平均渗透系数比较接近，即当含水层的渗透系数随深度衰减时，传统的常数渗透系数抽水试验模型采用井筒水位降深数据反演的渗透系数与随深度衰减的渗透系数平均值近似.

4 结论

抽水试验被广泛用于反求含水层水动力学参数，其参数反演的精度不仅仅取决于野外试验操作，还取决于参数反演所采用的模型.然而，现有的模型往往假设含水层渗透系数是常数.事实上，大量的野外试验表明实际含水层中的渗透系数呈现出随深度衰减的现象，该现象无法用现有的抽水试验模型来刻画.为此，本文建立考虑越流、井储效应和随深度衰减的渗透系数的抽水试验数学模型，并采用拉普拉斯变换、Goldstein-Weber变换和格林函数法推导半解析解，作为现有的模型的拓展.主要结论如下.

（1）在抽水试验过程中，随深度衰减的渗透系数、越流补给和井储效应对结果的影响不能忽视：当含水层渗透系数随深度衰减越快，抽水试验过程中井筒内的降深越大，导致含水层中存在一个更大的影响范围（降落漏斗更大）；越流系数 R越小，井筒内抽水试验降深越小；井筒内地下水的储量越小，抽水过程中井筒内降深越大.

（2）相比均质介质含水层中，当含水层渗透系数随深度衰减时，滤水管的位置对抽水试验结果的影响比较显著，滤水管在含水层顶部时的井筒内降深比滤水管在含水层底部时的井筒水位降深要小.

（3）理论抽水试验结果表明，当含水层的渗透系数随深度衰减时，恒定渗透系数模型采用井筒水位降深数据反演的渗透系数是实际含水层平均渗透系数的近似结果.

参考文献

[1] Achtziger-Zupančič,P.,Loew,S.,Mariéthoz,G.,2017.A New Global Database to Improve Predictions of Permeability Distribution in Crystalline Rocks at Site Scale.Journal of Geophysical Research: Solid Earth,122(5): 3513-3539.https://doi.org/10.1002/2017jb014106

[2] Cassiani,G.,Kabala,Z.J.,Medina,M.A.,1999.Flowing Partially Penetrating Well: Solution to a Mixed-Type Boundary Value Problem.Advances in Water Resources,23(1):59-68.https://doi.org/10.1016/s0309-1708(99)00002-0

[3] Chen,C.,Wen,Z.,Zhou,H.,et al.,2020.New Semi-Analytical Model for an Exponentially Decaying Pumping Rate with a Finite-Thickness Skin in a Leaky Aquifer.Journal of Hydrologic Engineering,25(8).https://doi.org/10.1061/(asce)he.1943-5584.0001956

[4] Chen,C.X.,2020.Improvement of Dupuit Model:With Infiltration Recharge.Hydrogeology & Engineering Geology,47(5):1-4 (in Chinese with English abstract).

[5] Chen,C.X.,Tang,Z.H.,2021.Improvement of the Theis Unsteady Well Flow Model with Infiltration Recharge in a Phreatic Aquifer.Hydrogeology & Engineering Geology,48(6):1-12 (in Chinese with English abstract).

[6] Chen,Y.F.,Ling,X.M.,Liu,M.M.,et al.,2018.Statistical Distribution of Hydraulic Conductivity of Rocks in Deep-Incised Valleys,Southwest China.Journal of Hydrology,566:216-226.https://doi.org/10.1016/j.jhydrol.2018.09.016

[7] Hantush,M.S.,1964.Advances in Hydroscience.Elsevier,Amsterdam,281-432.https://doi.org/10.1016/b978-1-4831-9932-0.50010-3

[8] Hantush,M.S.,1960.Modification of the Theory of Leaky Aquifers.Journal of Geophysical Research,65(11):3713-3725.https://doi.org/10.1029/jz065i011p03713

[9] Kabala,Z.J.,Cassiani,G.,1997.Well Hydraulics with the Weber-Goldstein Transforms.Transport in Porous Media,29(2):225-246.https://doi.org/10.1023/A:1006542203102

[10] Lai,R.Y.S.,Su,C.W.,1974.Nonsteady Flow to a Large Well in a Leaky Aquifer.Journal of Hydrology,22(3-4):333-345.https://doi.org/10.1016/0022-1694(74)90085-7

[11] Lin,Y.C.,Yang,S.Y.,Fen,C.S.,et al.,2016.A General Analytical Model for Pumping Tests in Radial Finite Two-Zone Confined Aquifers with Robin-Type Outer Boundary.Journal of Hydrology,540: 1162-1175.https://doi.org/10.1016/j.jhydrol.2016.07.028

[12] Louis, C. ，1972. Rock Hydraulics. In: Müller, L. ，ed.，Rock Mechanics. International Centre for Mechanical Sciences. Springer, Vienna.

[13] Luo,Y.,Liu,P.,Chen,W.X.,et al.,2023.Pumping Test to Interpret Karst Aquifer Types and Hydrogeological Parameters.Journal of China Hydrology,43(5):39-44 (in Chinese with English abstract).

[14] Mathias,S.A.,Todman,L.C.,2010.Step-Drawdown Tests and the Forchheimer Equation.Water Resources Research,46(7):W07514.https://doi.org/10.1029/2009wr008635

[15] Novakowski,K.S.,1989.A Composite Analytical Model for Analysis of Pumping Tests Affected by Well Bore Storage and Finite Thickness Skin.Water Resources Research,25(9):1937-1946.https://doi.org/10.1029/wr025i009p01937

[16] Papadopulos,I.S.,Cooper,H.H.Jr.,1967.Drawdown in a Well of Large Diameter.Water Resources Research,3(1):241-244.https://doi.org/10.1029/wr003i001p00241

[17] Saar,M.O.,Manga,M.,2004.Depth Dependence of Permeability in the Oregon Cascades Inferred from Hydrogeologic,Thermal,Seismic,and Magmatic Modeling Constraints.Journal of Geophysical Research:Solid Earth,109(B4).https://doi.org/10.1029/2003jb002855

[18] Sakata,Y.,Ikeda,R.,2013.Depth Dependence and Exponential Models of Permeability in Alluvial-Fan Gravel Deposits.Hydrogeology Journal,21(4):773-786.https://doi.org/10.1007/s10040-013-0961-8

[19] Theis,C.V.,1935.The Relation between the Lowering of the Piezometric Surface and the Rate and Duration of Discharge of a Well Using Ground-Water Storage.Transactions,American Geophysical Union,16(2): 519-524.https://doi.org/10.1029/tr016i002p00519

[20] Tsai,C.S.,Yeh,H.D.,2012.Wellbore Flow-Rate Solution for a Constant-Head Test in Two-Zone Finite Confined Aquifers.Hydrological Processes,26(21): 3216-3224.https://doi.org/10.1002/hyp.8322

[21] Wan,L.,Jiang,X.W.,Wang,X.S.,2010.A Common Regularity of Aquifers: The Decay in Hydraulic Conductivity with Depth.Geological Journal of China Universities,16(1):7-12 (in Chinese with English abstract).

[22] Wang,Q.R.,Shi,W.G.,Zhan,H.B.,et al.,2018.Models of Single-Well Push-Pull Test with Mixing Effect in the Wellbore.Water Resources Research,54(12):10155-10171.https://doi.org/10.1029/2018wr023317

[23] Wang,Q.R.,Zhan,H.B.,2015.On Different Numerical Inverse Laplace Methods for Solute Transport Problems.Advances in Water Resources,75:80-92.https://doi.org/10.1016/j.advwatres.2014.11.001

[24] Wei,S.Y.,Wan,J.W.,Sun,W.,et al.,2022.Identification of Hydrogeological Parameter in Double-Porosity Fractured Aquifer Based on Single-Well Pumping Test.Earth Science(in press)(in Chinese with English abstract).

[25] Wang,X.S.,Jiang,X.W.,Wan,L.,et al.,2009.Evaluation of Depth-Dependent Porosity and Bulk Modulus of a Shear Using Permeability-Depth Trends.International Journal of Rock Mechanics and Mining Sciences,46(7):1175-1181.https://doi.org/10.1016/j.ijrmms.2009.02.002

[26] 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

[27] Wen,Z.,Huang,G.H.,Zhan,H.B.,2011.Non-Darcian Flow to a Well in a Leaky Aquifer Using the Forchheimer Equation.Hydrogeology Journal,19(3):563-572.https://doi.org/10.1007/s10040-011-0709-2

[28] 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

[29] Wen,Z.,Huang,G.H.,Liu,Z.T.,et al.,2011.An Approximate Analytical Solution for Two-Region Non-Darcian Flow Toward a Well in a Leaky Aquifer.Earth Science,36(6):1165-1172 (in Chinese with English abstract).

[30] Yang,S.Y.,Yeh,H.D.,Chiu,P.Y.,2006.A Closed Form Solution for Constant Flux Pumping in a Well under Partial Penetration Condition.Water Resources Research,42(5).https://doi.org/10.1029/2004wr003889

[31] Yang,S.Y.,Yeh,H.D.,2005.Laplace-Domain Solutions for Radial Two-Zone Flow Equations under the Conditions of Constant-Head and Partially Penetrating Well.Journal of Hydraulic Engineering,131(3):209-216.https://doi.org/10.1061/(asce)0733-9429(2005)131:3(209)

[32] Zhang,L.,Franklin,J.A.,1993.Prediction of Water Flow into Rock Tunnels:an Analytical Solution Assuming an Hydraulic Conductivity Gradient.International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,30(1):37-46.https://doi.org/10.1016/0148-9062(93)90174-c

[33] Zheng,G.,Ha,D.,Loaiciga,H.,et al.,2019.Estimation of the Hydraulic Parameters of Leaky Aquifers Based on Pumping Tests and Coupled Simulation/Optimization:Verification Using a Layered Aquifer in Tianjin,China.Hydrogeology Journal,27(8):3081-3095.https://doi.org/10.1007/s10040-019-02021-z

[34] Zhuang,C.,Li,Y.B.,Zhou,Z.F.,et al.,2022.Effects of Exponentially Decaying Aquitard Hydraulic Conductivity on Well Hydraulics and Fractions of Groundwater Withdrawal in a Leaky Aquifer System.Journal of Hydrology,607:127439.https://doi.org/10.1016/j.jhydrol.2022.127439

[35] Zhou,Z.F.,2018.Application of the Monoexponential Decay Function on Analyzing Water-Inflow in Underground Boreholes.Journal of Mining & Safety Engineering,35 ( 3 ): 642 - 648 (in Chinese with English abstract).

基金资助

国家重点研发计划项目(2017YFC0804108)

中煤科工西安研究院(集团)有限公司科技创新基金项目(2021XAYKF02)

湖北省自然科学基金(2021CFA089)