有限差分法模拟潜水含水层井流的水位校正与流线识别

Water Level Correction and Streamline Tracing in Phreatic Aquifer Pumping Well Flow Simulation Using Finite Difference Method

摘要

中心网格有限差分法模拟水位与流线是开展众多水文地质工作的基础.为了解决潜水含水层中抽水井处水位模拟误差问题并精确识别流线.本研究基于复势理论，推导了中心网格有限差分法抽水井处水位校正公式，耦合数值模型量化分析了水位校正对流线精度的影响.结果表明，校正公式可显著提升抽水井处水位模拟精度；同时明确了当离散网格尺寸约为井半径的5倍时，无需校正即可同时保证水位与流线的模拟精度.该研究为解决数值模拟中离散化带来的模拟误差提供了理论依据与实用判据.

Abstract

Block-centered finite difference method for simulating water levels and streamlines is the basis of numerous hydrogeological works. To address the simulation error of water levels at pumping wells in a phreatic aquifer and accurately identify streamlines, this study derives a water level correction fomula for pumping wells in the block-centered finite difference method based on complex potential theory. A coupled numerical model is developed to quantitatively analyze the impact of water level correction on the accuracy of streamlines. The results demonstrate that the correction formula significantly improves the simulation accuracy of water levels at pumping wells. Furthermore, it is clarified that when the discrete grid size is approximately 5 times the well radius, both water level and streamline simulation accuracy can be ensured without the need for correction. This study provides a theoretical basis and practical criterion for mitigating simulation errors caused by discretization in numerical modeling.

Graphical abstract

关键词

中心网格有限差分 / 潜水含水层 / 抽水井 / 水位 / 流线 / 环境地质学.

Key words

block-centered finite difference / phreatic aquifer / pumping well / water level / streamline / environmental geology

0 引言

精确模拟抽水井附近地下水位与流线是开展地下水资源管理、地下水源地保护区划分、地下水污染修复等众多水文地质工作的基础（刚什婷和邓英尔， 2015； 汪洋等， 2018； 李旭等， 2022； 王新港， 2023； 苏世林等， 2024）.其方法有解析与数值模拟两种.解析法通常是用复势函数表示地下水等水位线与流线，但使用条件苛刻，难以刻画复杂边界条件的地下水流动（Zhang， 2024）.数值模拟法适用于复杂水文地质条件，是绘制地下水流场最常用的方法（Poeter and Hsieh， 2020）.在模拟地下水位与流线的众多模型中，中心网格有限差分模型因求解速度快，成熟度高而被广泛应用（杨建宏， 2008），例如GMS（groundwater modeling system）软件中的MODFLOW与MODPATH模型，MODFLOW是通过中心网格有限差分方式模拟地下水位空间分布；MODPATH是基于MODFLOW水位模拟结果追踪粒子运移轨迹，识别出地下水流线.

中心网格有限差分法模拟已知井周半径大小的抽水井处水位往往存在误差（夏强等， 2007； MacMillan and Schumacher， 2015）.在正方形中心网格有限差分模型中，网格中心水位相当于0.2倍正方形网格边长处的水位（Prickett， 1967； Lerner， 1989）.中心网格有限差分法模拟承压含水层中抽水井与自流井处水位的校正公式已经被前人提出（陈崇希等， 2007； 王旭升， 2008），但潜水含水层中抽水井处模拟水位的校正问题尚未解决.

有限差分法识别流线是基于有限差分法模拟水位的水力梯度并以一定时间或空间步长的正向或逆向粒子追踪识别粒子轨迹运移的过程（Pollock， 2016），识别的流线精度受抽水井处水位模拟误差的影响程度尚不明确.识别地下水流线的首要任务是释放粒子追踪起点.例如，在GMS中，MODPATH模型无法输入井周半径参数，在抽水井上释放粒子追踪起点时，MODPATH仅可以设置粒子数量，而默认粒子按照某种规律规则分布，但抽水井上粒子追踪起点分布规则现有文献资料并未报道.

本文将以地下水模拟软件GMS为例，通过复势函数推导矩形中心网格有限差分法模拟平面二维潜水含水层抽水井处水位校正公式，明确MODPATH中抽水井上粒子追踪起点分布规则，给出抽水井处水位与流线模拟误差的解决办法，为精确绘制地下水流场提供有力的理论依据和技术支持.

1 理论方法

1.1　复势函数

根据复势理论，平面二维稳态地下水流场中的等水位线与流线可用流量势（后文简称势函数）与流函数来表示（Campos， 2010； Strack， 1989）.复势函数（Ω）由实部（Φ）与虚部（Ψ）两部分组成，形如：

\Omega =\Phi +i\Psi，(1)

式（1）中，Φ、Ψ分别为势函数和流函数，L³/T.

在含水层水平、垂向无补给排泄、Dupuit假设条件下，均质各向同性介质潜水含水层中势函数与地下水位之间的关系可以表示为（Strack， 1989）：

\Phi =\frac{\mathrm{1}}{\mathrm{2}}K{h}^{\mathrm{2}}，(2)

式（2）中，K为含水层渗透系数，L/T；h为潜水含水层中相对含水层底板的地下水位高度，L.

在笛卡尔坐标系中，均匀侧向流和抽水井共同扰动的流场符合线性叠加关系，势函数与流函数可分别表示为（Strack， 1989）：

\Phi =-q(X\mathrm{c}\mathrm{o}\mathrm{s}\alpha +Y\mathrm{s}\mathrm{i}\mathrm{n}\alpha )+\frac{Q}{\mathrm{4}\pi }\mathrm{l}\mathrm{n}\left[{\left(X-X_p\right)}^{\mathrm{2}}+{\left(Y-Y_p\right)}^{\mathrm{2}}\right]+{\Phi }_{\mathrm{0}},(3a)

\Psi =q(X\mathrm{s}\mathrm{i}\mathrm{n}\alpha -Y\mathrm{c}\mathrm{o}\mathrm{s}\alpha )+\frac{Q}{\mathrm{2}\pi }\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{2}(Y-Y_p,X-X_p)+{\Psi }_{\mathrm{0}},(3b)

式（3）中，q为均匀侧向流的单宽流量，L²/T；α为均匀侧向流方向，沿X正半轴方向为0基准线，沿逆时针方向α逐渐增大，α取值范围为［0， 2π］；Q为抽水井流量，L³/T，Q为正数（Q为负数表示注水井）；（X_p， Y_p ）为抽水井位置；Φ₀、Ψ₀是常量，表示基准势函数与流函数，L³/T.值得注意的是，多井扰动的流场可由叠加单个抽水井得到.若单宽流量q沿X、Y方向的分量记为q_X、q_Y，则潜水含水层中的q_X、q_Y 可分别表示为：

q_X=-Kh\frac{\partial h}{\partial X}=Kh{J}_X，

q_Y=-Kh\frac{\partial h}{\partial Y}=Kh{J}_Y，(4)

式（4）中，J_X 、J_Y 是X、Y方向上的水力梯度，当潜水面可近似水平时，h、J（J_X 或J_Y ）可近似为常数.

1.2　矩形中心网格有限差分法抽水井处校正公式

1.2.1　抽水井处水位模拟原理

在中心网格有限差分模型中，抽水井位于网格中心的格点位置.网格大小用ΔX、ΔY表示（ΔX、ΔY的量纲为L），字母i、j和数字表示格点位置与网格编号，如图1所示.

网格（i， j）与周边相邻网格水量均衡方程为（Ramadhan， 2015）：

Q=\frac{K}{\mathrm{\Delta }X}\left[A_{\mathrm{1}}\left(h_{i-\mathrm{1},j}-h_{i,j}\right)+A_{\mathrm{2}}\left(h_{i+\mathrm{1},j}-h_{i,j}\right)\right]+\frac{K}{\mathrm{\Delta }Y}\left[A_{\mathrm{3}}\left(h_{i,j-\mathrm{1}}-h_{i,j}\right)+A_{\mathrm{4}}\left(h_{i,j+\mathrm{1}}-h_{i,j}\right)\right],(5)

式（5）中，A_n（n=1，…，4）为过水断面.需要注意的是，为了保证现有中心网格有限差分模型中Q与复势函数中Q符号一致，式（5）中的Q取正数表示抽水井，取负数表示注水井.

对于潜水含水层，过水断面由饱和带厚度表示为：

A_{\mathrm{1}}=\mathrm{\Delta }Y{h}_{i-\mathrm{1},j}，A_{\mathrm{2}}=\mathrm{\Delta }Y{h}_{i+\mathrm{1},j}，

A_{\mathrm{3}}=\mathrm{\Delta }X{h}_{i,j-\mathrm{1}}，A_{\mathrm{4}}=\mathrm{\Delta }X{h}_{i,j+\mathrm{1}} .(6)

在计算均质各向同性介质模型的网格间流量时，现有中心网格有限差分法将过水断面处理为（Mcdonald and Harbaugh， 1988）：

A_{\mathrm{1}}=\mathrm{\Delta }Y\frac{h_{i-\mathrm{1},j}+h_{i,j}}{\mathrm{2}}，A_{\mathrm{2}}=\mathrm{\Delta }Y\frac{h_{i+\mathrm{1},j}+h_{i,j}}{\mathrm{2}}，

A_{\mathrm{3}}=\mathrm{\Delta }X\frac{h_{i,j-\mathrm{1}}+h_{i,j}}{\mathrm{2}}，A_{\mathrm{4}}=\mathrm{\Delta }X\frac{h_{i,j+\mathrm{1}}+h_{i,j}}{\mathrm{2}} .(7)

将式（7）代入式（5），整理得

Q=K\mathrm{\Delta }Y\left(\frac{h_{i+\mathrm{1},j}^{\mathrm{2}}-h_{i,j}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }X}+\frac{h_{i-\mathrm{1},j}^{\mathrm{2}}-h_{i,j}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }X}\right)+K\mathrm{\Delta }X\left(\frac{h_{i,j+\mathrm{1}}^{\mathrm{2}}-h_{i,j}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }Y}+\frac{h_{i,j-\mathrm{1}}^{\mathrm{2}}-h_{i,j}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }Y}\right) .(8)

这就是中心网格有限差分方程.

1.2.2　校正公式

在均质各向同性潜水含水层中，不考虑垂向补给与排泄、井损与表皮效应，基于Dupuit假设条件下的平面二维潜水稳定流中，抽水井所在网格格点与周边相邻网格格点的水位二次方满足线性关系（Bear， 1979），由式（2）与式（3a）可得：

h_{i+\mathrm{1},j}^{\mathrm{2}}-{h_w}^{\mathrm{2}}=\frac{Q}{\pi K}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }X}{r_w}\right)，h_{i,j+\mathrm{1}}^{\mathrm{2}}-{h_w}^{\mathrm{2}}=\frac{Q}{\pi K}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }Y}{r_w}\right)，

h_{i-\mathrm{1},j}^{\mathrm{2}}-{h_w}^{\mathrm{2}}=\frac{Q}{\pi K}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }X}{r_w}\right)，h_{i,j-\mathrm{1}}^{\mathrm{2}}-{h_w}^{\mathrm{2}}=\frac{Q}{\pi K}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }Y}{r_w}\right).(9)

对于一般矩形网格，ΔX≠ΔY，将式（9）代入式（7），得到矩形中心网格有限差分法潜水含水层抽水井处水位校正公式为：

h_w^{\mathrm{2}}=h_{i,j}^{\mathrm{2}}-\frac{Q}{\pi K}\frac{\mathrm{1}}{\mathrm{\Delta }{X}^{\mathrm{2}}+\mathrm{\Delta }{Y}^{\mathrm{2}}}\left(\mathrm{\Delta }{Y}^{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }X}{r_w}\right)+\mathrm{\Delta }{X}^{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }Y}{r_w}\right)-\pi \mathrm{\Delta }X\mathrm{\Delta }Y\right) .(10)

若网格大小ΔX=ΔY，由式（10）整理可得正方形中心网格有限差分法潜水含水层抽水井处水位校正公式为：

h_w^{\mathrm{2}}=h_{i,j}^{\mathrm{2}}-\frac{Q}{\pi K}\left(\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }X}{r_w}\right)-\frac{\pi }{\mathrm{2}}\right) .(11)

1.3　MODPATH流线识别原理

1.3.1　核心算法原理

MODPATH识别地下水流线的核心算法是基于中心网格有限差分的粒子追踪法，该算法最早由David W. Pollock于1988年开发（Pollock， 1988），也称为Pollock法.Pollock算法的关键思想是假设两点间地下水速度呈线性分布（Ramadhan， 2015），具体原理如下.

在平面二维地下水流场中，如果粒子从点（X_p， Y_p ）运移到点（X_p+1， Y_p+1）的时间步长为ΔT，并且以（X_p， Y_p ）处的速度V（X_p， Y_p ）运移，则X，Y方向的位移分量可分别计算为：

\mathrm{\Delta }X=V_X\left(X_p,Y_p\right)\mathrm{\Delta }T，\mathrm{\Delta }Y=V_Y\left(X_p,Y_p\right)\mathrm{\Delta }T，(12)

式（12）中，V_X，V_Y 分别为X，Y方向的速度分量.假设点（X_p， Y_p ）周围流速的空间分布是线性的，则其变化梯度可表示为：

A_X=\frac{V_X\left(X+\mathrm{\Delta }X,Y+\mathrm{\Delta }Y\right)-V_X\left(X,Y\right)}{\mathrm{\Delta }X}，

A_Y=\frac{V_Y\left(X+\mathrm{\Delta }X,Y+\mathrm{\Delta }Y\right)-V_Y\left(X,Y\right)}{\mathrm{\Delta }X}，(13)

式（13）中，A_X 和A_Y 在一个步长计算中为常数，分别表示X和Y方向上的流速梯度.按照这种流速变化规律进行粒子追踪，则下一步的实际位置为：

X_{p+\mathrm{1}}=X_p+V_X(X_p,Y_p)\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(A_X\mathrm{\Delta }T\right)-\mathrm{1}}{A_X}，

Y_{p+\mathrm{1}}=Y_p+V_Y(X_p,Y_p)\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(A_Y\mathrm{\Delta }T\right)-\mathrm{1}}{A_Y} .(14)

在计算过程中，A_X 或A_Y 可能趋近于0，此时应注意以下极限公式：

\underset{A\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{e}\mathrm{x}\mathrm{p}(A\mathrm{\Delta }T)-\mathrm{1}}{A}=\mathrm{\Delta }T,(15)

式（15）中，A可以取A_X 或A_Y .这就是MODPATH粒子追踪的核心算法，此外，MODPATH还通过粒子进出网格单元位置优化粒子追踪步长ΔT，提高流线识别精度.

1.3.2　井上粒子追踪起点分布规则

合理布设粒子追踪起点位置是识别地下水流线的首要任务.在GMS软件中，粒子追踪起点布设方式有3种，分别为：（1）在网格单元空间内布设粒子追踪起点；（2）导入外部散点坐标数据作为粒子追踪起点；（3）在抽水井上布设粒子追踪起点（Pollock， 2016； 李若怡， 2022）.第（1）种粒子追踪起点的分布规律是在网格空间内按行和列空间均匀分布，第（2）种粒子追踪起点分布规则取决于散点坐标数据，第（3）种粒子追踪起点数量最大限值为100.笔者通过大量模拟实验最终确定，在MODPATH模型中，抽水井上释放的粒子追踪起点（P _n ）从X轴正方向出发沿逆时针均匀分布于椭圆上，如图2所示，椭圆大小与抽水井所在网格大小相关，即椭圆的轴长为抽水井所在网格长宽的一半：

\frac{\mathrm{16}{\left(X-X_{i,j}\right)}^{\mathrm{2}}}{\mathrm{\Delta }{X}^{\mathrm{2}}}+\frac{\mathrm{16}{\left(Y-Y_{i,j}\right)}^{\mathrm{2}}}{\mathrm{\Delta }{Y}^{\mathrm{2}}}=\mathrm{1}.(16)

当ΔX=ΔY时，椭圆方程简化为圆形方程，

\frac{{\left(X-X_{i,j}\right)}^{\mathrm{2}}}{{\left(\frac{\mathrm{\Delta }X}{\mathrm{4}}\right)}^{\mathrm{2}}}+\frac{{\left(Y-Y_{i,j}\right)}^{\mathrm{2}}}{{\left(\frac{\mathrm{\Delta }X}{\mathrm{4}}\right)}^{\mathrm{2}}}=\mathrm{1}.(17)

而抽水井的井周圆形方程为：

\frac{{\left(X-X_{i,j}\right)}^{\mathrm{2}}}{{r_w}^{\mathrm{2}}}+\frac{{\left(Y-Y_{i,j}\right)}^{\mathrm{2}}}{{r_w}^{\mathrm{2}}}=\mathrm{1} .(18)

联立方程（17）与（18）得到，ΔX=ΔY=4r_w 时，粒子追踪起点分布于抽水井井周上.

2 案例分析

2.1　案例设计

为了验证中心网格有限差分模拟潜水含水层中抽水井处水位校正公式的合理性、定量分析抽水井处水位模拟误差对MODPATH识别流线的影响，本文设计2个典型案例A和B，剖面示意如图3.案例A：在均质各向同性介质的潜水含水层中，含水层可近似为圆岛；地面近似水平且距离隔水底板的高度h_g =50 m，潜水面水平；圆岛边界为定水头边界，定水头为25 m；初始水位h₀=25 m；隔水底板水平，且被视为基准面，水位h_r =0 m；坐标原点（0 m， 0 m）处有一口抽水井，井滤管贯穿整个含水层，抽水流量Q=100 m³/d，井半径r_w =0.1 m，忽略井损与表皮效应；不考虑降水和蒸发的垂向补给作用，地下水流动可近似为平面二维稳定流；渗透系数K=1 m/d.案例B：考虑平面无限大含水层存在侧向流，侧向流方向水平，抽水井开采前的背景潜水面水力梯度J=1/200 000，渗透系数K=10 m/d.

在MODFLOW中，模型范围设置为100 km × 100 km.抽水井位于模拟区中心，网格大小一般为1 000 m × 1 000 m.为了分析抽水井处水位模拟精度及其对流线模拟结果的影响，在案例A中，以圆周为常水头边界，水头值为25 m，对抽水井所在网格进行局部加密，加密后抽水井所在网格保持为正方形，边长ΔX取值为0.32~1 000.00 m，如图4a所示.在案例B中，MODFLOW模型上下两侧边界默认为零通量边界，左右两侧边界设置为定水头边界，水头值分别为25.25 m、24.75 m，如图4所示.

2.2　井中水位模拟结果

通过案例A分析，得到抽水井处水位模拟值，真实解（基于复势理论的解析解），潜水含水层抽水井处水位校正公式（11）修正的水位校正值，结果如图5所示.对比抽水井处水位模拟值与真实值，得到以下结果：模拟值随网格尺寸ΔX减小而减小，模拟值与ΔX呈对数关系，如图5所示.当ΔX 4.81r_w 时，抽水井处水位模拟值大于真实值，如：ΔX=1 000 m时，抽水井处水位模拟值为21.19 m，比抽水井处水位真实值14.40 m偏大6.80 m；当ΔX=4.81r_w =0.48 m时，抽水井处水位模拟值约14.35 m，与真实值14.40 m基本一致；当ΔX4.81r_w 时，抽水井处水位模拟值小于真实值；如：ΔX=0.32 m时，抽水井处水位模拟值为13.67 m，模拟值比真实值14.40 m偏小0.72 m.

2.3　流线识别结果

在案例B中，由式（3b）计算得到真实地下水流线，地下水流线在空间分布上关于Y=0对称，如图6a所示.流函数值等于0的流线是上游地下水汇入抽水井的最外围流线，即地下水能否被抽水井抽出的分界线.井周上位于分界线的点是（0.1 m， 0 m），井周外，当Y=0时，分界线上X取最大值12 700 m；当X=-10 000 m时，分界线上的纵坐标为Y=25 000 m.

在抽水井井周上均匀布设17个粒子追踪起点，若以角度β表示井周上粒子追踪起点位置，则相邻两个起点与井点连线间的夹角为π/8.通过GMS软件中MODFLOW与MODPATH相结合模拟得到地下水流线，如图6b和6c所示.当ΔX=200 m时，水位模拟误差使MODPATH识别的流线与真实流线明显存在区别.如：粒子追踪起点位于β=0的流线，X=-10 000 m，真实流线的纵坐标为25 000 m，而MODPATH识别流线的纵坐标仅为0.7 m；粒子追踪起点位于β=π/2的流线，X=-10 000 m，真实流线的纵坐标为10 000 m，而MODPATH识别的流线是分界线，纵坐标为1 247 m，如图6b.当ΔX=0.32 m时，MODPATH识别的所有流线与真实流线基本一致，X=-10 000 m，MODPATH识别流线的纵坐标为24 122 m，相对误差为4%，如图6c.

以X=-10 000 m时流线的纵坐标为参数，分析MODPATH模拟流线结果随ΔX变化特征.对比流线MODPATH模拟解与真实解，如图6d结果表明：ΔX≤50 m时，流线识别结果可靠；ΔX50 m时，流线识别结果严重偏离真实流线.

3 讨论

3.1　抽水井处水位校正与网格大小

当抽水井处水位与网格格点模拟水位相同时，无需再对抽水井处水位模拟结果进行校正.根据潜水含水层抽水井处水位校正公式（11），ΔX与r_w 满足关系式：

\mathrm{l}\mathrm{n}\frac{\mathrm{\Delta }X}{r_w}-\frac{\pi }{\mathrm{2}}=\mathrm{0},(19)

求得，

\mathrm{\Delta }X={\mathrm{e}}^{\pi /\mathrm{2}}{r}_w=\mathrm{4.81}{r}_w .(20)

这就是潜水含水层抽水井处水位校正值与真实值一致的理论依据，也为中心网格有限差分法精确模拟三维地下水流场中抽水井处水位时，抽水井所在网格剖分精度提供了依据.

抽水井处水位校正值精度受模拟精度直接影响，与模拟值偏大或偏小无关.抽水井周围网格局部加密区越大，网格尺寸越小，抽水井处水位模拟值精度越高，相应的校正值精度越高，如图7所示.

3.2　非均质与非稳定流潜水含水层抽水井处水位校正

式（10）与式（11）仅适用于一般矩形与正方形网格剖分的平面二维均质各向同性潜水含水层模型.推导式（10）与式（11）的核心在于潜水含水层有限差分模型求解的关键问题在于过水断面饱和带厚度的处理.在MODFLOW中，有且仅有隔水底板水平的均质各向同性含水层，中心网格有限差分方程将过水断面饱和带厚度处理成相邻两个格点处水位的平均值（Mcdonald and Harbaugh， 1988）.本文中的潜水含水层抽水井处水位校正公式不适用于沿X或Y方向渗透系数呈现非均质性或隔水底板明显起伏的潜水含水层.

对于均值各向异性含水层，K_X¹K_Y，将K_X 、K_Y 分别代入式（8）与（9），可得均值各向异性潜水含水层抽水井处水位校正公式为：

h_w^{\mathrm{2}}=h_{i,j}^{\mathrm{2}}-\frac{Q}{\pi }\frac{\mathrm{1}}{K_Y\mathrm{\Delta }{X}^{\mathrm{2}}+K_X\mathrm{\Delta }{Y}^{\mathrm{2}}}\left(\mathrm{\Delta }{Y}^{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }X}{r_w}\right)+\mathrm{\Delta }{X}^{\mathrm{2}}\mathrm{l}\mathrm{n}\left(\frac{\mathrm{\Delta }Y}{r_w}\right)-\pi \mathrm{\Delta }X\mathrm{\Delta }Y\right) .(21)

此外，本文抽水井处水位校正公式仅适用于中心网格有限差分法的平面二维均质潜水含水层抽水井处水位校正，对于其他有限差分求解法（如layer property flow）及非稳定流条件的适用性有待进一步研究.

3.3　中心网格有限差分模型与MODPATH模型抽水井所在网格适宜大小

中心网格有限差分模型与MODPATH精确模拟结果的网格大小与抽水井井周半径关系不一致.据抽水井处水位校正公式（式（10）或式（11）），中心网格有限差分法模拟抽水井处精确水位，网格大小ΔX与抽水井井周半径r_w 满足关系ΔX=4.81r_w；而MODPATH井周上粒子追踪起点分布规则（式（18））为：网格大小ΔX与抽水井井周半径r_w 满足关系ΔX=4r_w .两者不一致的原因是抽水井处水位校正公式主要考虑的是修正网格剖分造成的模拟误差，未充分考虑有限差分法求解与模型边界刻画等误差因素.

据式（11），抽水井处水位模拟值无需校正时的网格大小为4.81r_w，水位为：

h_{w\mathrm{1}}^{\mathrm{2}}=h_{i,j}^{\mathrm{2}}-\frac{Q}{\pi K}\left(\mathrm{l}\mathrm{n}\left(\mathrm{4.81}\right)-\frac{\pi }{\mathrm{2}}\right) .(22)

如果抽水井所在网格大小为4r_w，则有：

h_{w\mathrm{2}}^{\mathrm{2}}=h_{i,j}^{\mathrm{2}}-\frac{Q}{\pi K}\left(\mathrm{l}\mathrm{n}\mathrm{4}-\frac{\pi }{\mathrm{2}}\right).(23)

将网格大小为4r_w 时校正水位相对初始水位模拟误差控制在0.1内，则水位平方相对误差约0.01，则有：

\frac{h_{w\mathrm{2}}^{\mathrm{2}}-h_{w\mathrm{1}}^{\mathrm{2}}}{h_{\mathrm{0}}^{\mathrm{2}}}=\frac{\mathrm{0.184}\mathrm{4}Q}{\pi K{h}_{\mathrm{0}}^{\mathrm{2}}}\le \mathrm{0.01}.(24)

在实际应用中，若将抽水井所在网格大小加密至4.81r_w，则无需再用水位校正公式（11）；如果近似将抽水井所在网格大小剖分为4r_w，相对误差控制在0.1内，则抽水井流量Q不宜超过0.054πKh₀².对于K=1 m/d、h₀=25 m的案例，Q≤106 m³/d.ΔX=4r_w 与4.81r_w 时，MODPATH识别的流线精度基本一致，如图6d所示.

4 结论

本文基于复势函数等数学理论，推导了中心网格有限差分模型的潜水含水层抽水井处水位校正公式；识别了MODPATH中抽水井上粒子追踪起点按照椭圆分布，椭圆轴长等于抽水井所在网格长宽的一半，讨论了潜水含水层抽水井处水位校正公式的适用性，以及在GMS软件应用中精确模拟抽水井处水位与流线的方法，得到以下结论.

（1）本文提出的中心网格有限差分法抽水井处水位校正公式能系统性地将潜水含水层抽水井处水位误差从米级降至厘米级，其精度由局部网格尺度控制.

（2）研究确定当抽水井所在网格尺寸约为井半径的5倍时，可同时满足MODFLOW直接输出高精度水位与MODPATH精确追踪下游流线的要求，为建立可靠、高效的地下水流与运移模型提供了明确的实践指南.

参考文献

[1] Bear,J.,1979.Hydraulics of Ground Water.McGraw⁃Hill Inc.,New York.

[2] Campos,L.M.B.C.,2010.Complex Analysis with Applications to Flows and Fields.CRC Press,Boca Raton.

[3] Chen,C.X.,Wang,X.S.,Hu,L.T.,2007.Emendation of Drawdown in Pumping Wells for Numerical Modeling of Groundwater Flow.Journal of Hydraulic Engineering,38(4):481-485 (in Chinese with English abstract).

[4] Gang,S.T.,Deng,Y.E.,2015.Comprehensive Summary of GMS Application for Groundwater Resource Evaluation and Management.Ground Water,37(2):33-36 (in Chinese with English abstract).

[5] Lerner,D.N.,1989.Predicting Pumping Water Levels in Single and Multiple Wells Using Regional Groundwater Models.Journal of Hydrology,105(1-2):39-55.https://doi.org/10.1016/0022⁃1694(89)90095⁃4

[6] Li,R.Y.,2022.Study on the Three⁃Dimensional Shape and Recharge Sources for Catchment Zone of a Pumping Well near Rivers (Dissertation).China University of Geosciences,Beijing (in Chinese with English abstract).

[7] 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).

[8] MacMillan,G.J.,Schumacher,J.,2015.Correction of Discretization Errors Simulated at Supply Wells.Ground Water,53(4):651-657.https://doi.org/10.1111/gwat.12254

[9] Mcdonald,M.G.,Harbaugh,A.W.,1988.A Modular Three⁃Dimensional Finite⁃Difference Ground⁃Water Flow Model.US Geological Survey,U.S.A.,83-875.

[10] Poeter,E.,Hsieh,P.,2020.Graphical Construction of Groundwater Flow Nets.The Groundwater Project:Guelph,Ontario,Canada.https://doi.org/10.21083/978⁃1⁃7770541⁃3⁃7

[11] Pollock,D.W.,1988.Semianalytical Computation of Path Lines for Finite⁃Difference Models.Groundwater,26(6):743-750.https://doi.org/10.1111/j.1745⁃6584.1988.tb00425.x

[12] Pollock,D.W.,2016.User Guide for MODPATH Version 7—A Particle Tracking Model for MODFLOW.Open⁃File Report 2016-1086.U.S.Geological Survey,Reston,VA.https://doi.org/10.3133/ofr20161086

[13] Prickett,T.A.,1967.Designing Pumped Well Characteristics into Electric Analog Models.Groundwater,5(4):38-46.https://doi.org/10.1111/j.1745⁃6584.1967.tb01625.x

[14] Ramadhan,M.,2015.A Semi⁃Analytical Particle Tracking Algorithm for Arbitrary Unstructured Grids.University of Waterloo,Waterloo.

[15] Strack,O.D.,1989.Groundwater Mechanics.Prentice⁃Hall,Englewood Cliffs.

[16] Su,S.L.,Li,X.,Guo,Q.,et al.,2024.Effect Mechanistic of Partially Penetrating Well on Single⁃Well Push⁃Pull Tests for Groundwater Velocity Estimation.Earth Science,49(1):288-298 (in Chinese with English abstract).

[17] Wang,Y.,Li,J.,Xi,B.D.,et al.,2018.Research on the Division Technology of Karst Groundwater Source Protection Areas Based on Numerical Simulation.Carsologica Sinica,37(6):799-809 (in Chinese with English abstract).

[18] Wang,X.G.,2023.Experimental Study on Hydraulic Control and Remediation of Trichloroethylene Pollution in Single Well Pumping and Infiltration Groundwater (Dissertation).Hebei University of Engineering,Handan (in Chinese with English abstract).

[19] Wang,X.S.,2008.Revised Model for Finite⁃Difference Modeling of Flowing Artesian Wells.Earth Science,33(1):112-116 (in Chinese with English abstract).

[20] Xia,Q.,Wang,X.S.,Lu,L.B.,2007.Characteristics of Error in Numerical Modeling of Well Flow with MODFLOW.Geotechnical Investigation Surveying,35(10):29-32,37(in Chinese with English abstract).

[21] Yang,J.H.,2008.Comparative Analysis of Three Fundamental Numerical Computation Methods.Science Technology Information,35:19 (in Chinese).

[22] Zhang,S.,Wang,X.S.,Han,P.F.,et al.,2024.New Particle Tracking Method for 2D Steady⁃State Groundwater Flow around Wells:A Modified Algorithm Using the Stream Function.Journal of Hydrology,632:130928.https://doi.org/10.1016/j.jhydrol.2024.130928

基金资助

水利部地下水保护重点实验室开放课题资助项目

河南省自然科学基金项目(252300423266)

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

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