河流附近非水平初始面潜水渗流分析

曹彭强，束龙仓，鲁程鹏

河海大学水文水资源与水利工程科学国家重点实验室，南京 210098

0 前言

河流－地下水系统事关河流与含水层之间水量交换评价、农田灌排水渠设计、河道基流衰减过程研究、潜水含水层参数评估、河流与地下水水质交换等方面的问题［1－5］，是地下水渗流力学中的研究热点和难点之一。从20世纪40年代开始，学者们［6－8］就河流－地下水系统提出了不同的模型，并对模型进行了求解，就河流入渗对潜水的影响进行了分析。现有文献［9－10］主要针对隔水底板水平的条件下假设初始自由面水平，或隔水底板倾斜条件下初始自由面平行于倾斜隔水底板，这些假设简化了描述初始面的数学方程，从而方便模型的求解；但在自然状态下，初始自由面非水平更具普遍性，因此，如何描述与计算非水平初始面成为渗流计算的关键。

笔者以河流附近具水平隔水底板的潜水含水层中的潜水渗流模型为例，构造初始面非水平条件下的一维潜水非稳定渗流模型；利用Boussinesq第一线性化方法，通过Laplace变换，再结合迭加原理，提出非水平初始面和参数的简易计算方法，导出模型的解，并根据误差分析讨论解的适用条件，为初始面非水平条件下的潜水渗流计算提供方法。

1 潜水初始面非水平模型的建立

假设一顺直河流完整切割含水层（图1），潜水含水层均质各向同性，具有水平隔水底板，侧向无限延伸；潜水初始自由面非水平hx，0；河流初始水位为h1，无限远处含水层水位为h∞，0；河流水位瞬时上升至h2，水位升幅为h2－h1，且上升后保持不变；潜水流可视为一维流。

上述水文地质概念模型是在经典的Ferris模型［11］基础上，增加了面平均入渗或蒸发条件，且假设潜水初始自由面hx，0非水平。把hx，0看作原河流水位及水平面水位都为h∞，0（称为最初初始面），河流水位变化h1－h∞，0，历经时间tN，含水层中形成的自由面。河流附近初始面非水平的半无限潜水含水层渗流控制方程为

边界条件及初始条件为

图1 河流附近潜水含水层Fig.1 Phreatic aquifer near river

其中：μ为给水度；K为渗透系数（m/d）；h为水位（m）；w为入渗强度（m/d，补给为正，排泄为负）；x为距河边距离（m）；t为河流水位波动时间（d）；tN为初始面形成时间（d）。

控制方程与边界条件及初始条件构成了河流附近初始面非水平的半无限潜水含水层渗流数学模型。其计算坐标如图1所示（对于宽深比较小的河流，将坐标原点定于水位上升后的水土交界处）。初始自由面hx，0非水平，模型（1）难以直接进行线性化求解。

2 潜水初始面非水平模型的求解

2.1 数学模型求解

模型（1）求解的关键点及难点在于初始自由面hx，0的确定。形成初始自由面hx，0的过程是个渗流过程，满足地下水流非稳定流规律。

河流水位常有一定涨落，呈阶梯变化或连续变化，为了便于计算，把河流水位变化曲线概化为阶梯状线段［12］。当水位变化相对于含水层厚度较小时，将水位连续变化过程看作数个瞬时变化的叠加［13］，参照初始自由面水平条件下河流附近潜水含水层地下水流非稳定流计算公式求解方法［14］，当hx，t－hx，0≤0.1hm（hm为潜水含水层的平均厚度）时，利用Boussinesq第一线性化方法，通过Laplace变换，得河流水位连续变化情况下河流附近的初始面计算公式：

其中：Δhi为第i次水位变化ti－1为第i－1次水位变化后的时刻；t0＝0；a＝Khm/μ，为导压系数；erf（）和erfc（）分别为误差函数和余误差函数。

河流水位瞬时上升Δh2并保持不变，由迭加原理，含水层中水位可由下式计算：

当h1＝h∞，0，即初始面水平时，Δhi＝0，式（3）变为

式（4）为初始面水平潜水含水层非稳定流计算公式。

式（2）中的hi与ti－1在河流水位变化过程中不易取得，难以计算，若将n次水位变化过程简化为1次水位变化过程，即式（2）变为

式（3）变为

式（6）即为简化后的初始面非水平潜水含水层非稳定流计算公式。

2.2 参数的确定

应用式（6）计算初始面非水平含水层水位时，需获得参数x、t、a和tN的取值。其中，x和t是确定的，对于某一次水位的变化，还需要确定参数a和tN。a和tN的具体求解过程有2种方法：查表法和配线法。

2.2.1 查表法

现已知河流初始水位h1、河流水位上升后水位h2、距离河流无穷远处含水层水位h∞，0。h∞，0即是离河流无穷远处的潜水水位，但在实际应用中只需取足够远、含水层水位影响较小处的水位即可。

式（6）减式（5）得

因此，式（7）简化为

即

式（9）中的未知数a，可选取离河流较近观测孔中的数据hx，t，查erfc（）数值表求得。再将a代入式（5），求得tN。

查表法计算过程简单，所需水位资料较少，只需要河流水位变化初期（t→0）的一个离河流较远的水位观测孔和一个离河流较近的观测孔的水位数据即可，适用于河岸地下水水位观测资料不足的情况。但由于erfc（）数值表中数据有限，在查表求解a时需进行插值计算。

2.2.2 配线法

令

对式（9）和式（10）同时取对数。对于同一个观测孔而言，lg（x2/4a）为一常数。则lg［（hx，t－hx，0）/（h2－h1）］－lg（t）曲线（实际资料曲线）和lg［erfc（λ）］－lg（1/λ2）曲线（标准曲线）的图形是相同的，只是横坐标平移了x2/4a而已。如果已知坐标值及水位观测值，就可以计算参数a。

具体步骤为：根据erfc（λ）数值表绘制lg［erfc（λ）］－lg（1/λ2）标准曲线（图2）；根据观测孔的实际水位观测数据在另一张有相同坐标的透明纸上绘制lg［（hx，t－hx，0）/（h2－h1）］－lg（t）实际资料曲线；将实际资料曲线叠置在标准曲线上，在保持对应坐标轴平行的条件下，移动坐标纸，直至2条曲线重合为止；之后在图上任取一点作为匹配点。读出该点的坐标值t、1/λ2，则

再将a代入式（5），计算不同点、不同时刻潜水水位。配线法对于只有河流水位未发生变化前观测孔的潜水水位数据的情况也适用，但在进行配线时会受个人主观因素影响，对于观测孔水位资料序列较短的情况，不同人给出的配线数值差别较大。

图2 lg［erfc（λ）］－lg（1/λ2）标准曲线Fig.2 lg［erfc（λ）］－lg（1/λ2）standard curve

2.3 简化后解的误差分析

式（6）减式（3）得x处潜水初始面计算误差ε为

其中：ε是x，a，tN，t，ti－1的 函 数，ε＝f（x，a，tN，t，ti－1）；0≤ti－1≤tN。当Δhi＞0时，ε≤0，即用1次水位变化过程代替多次水位变化过程计算所得的初始潜水水位小于实际观测值。

式（12）中ε分别对x，a，tN，t，ti－1求导，得ε与各参量的关系（图3）。由图3可知：随x增大，ε先增后减，在x＝0和x＝∞处，ε＝0，即在离河流较近距离和足够远处，采用式（6）计算的误差较小；随a增大，ε先增后减，当a较小时，ε较小，也即当含水层K和hm较小、μ较大时，计算的误差较小；随tN增大，ε先增后减，当tN＝0时，ε＝0；随t增大，ε先增后减，即水位变化初期和末期，计算误差较小；ε随ti－1增大而增大，形成初始面的时间越长，后续水位变化计算的水位误差越大。

由以上分析可得，在河流水位变化初期，水位变幅相对于潜水含水层厚度较小、观测孔距河岸较近，即t较小、Δhi与x较小时，计算的ε较小，才能用式（6）代替式（3）进行计算。

3 算例

图3 ε与参数关系Fig.3 Relationship ofεand parameters

某河流在某时刻的水位从23.00m上升为25.00m，并保持不变。垂直河岸布置了2个观测孔（图1），观测水位资料见表1，含水层入渗w为0。

3.1 查表法

已知河流原水位h1＝23.00m；河流水位上升后水位h2＝25.00m。

由表1可知，河流水位发生变化后，离河岸较远的观测孔（K446.4）的水位变化极小，可以把该点t＝0时刻的水位作为最初初始水位，即H∞，0＝22.16m。

选取t值较小（1.1d）时、离河较近的观测孔（K55.4）的水位数据代入式（5），即x＝55.4m，t＝1.1d，hx，0＝22.87m，hx，t＝23.02m，得erfc（λ1）＝7.50×10－3。

查erfc（λ1）数值表，插值得λ1≈1.26，则a1＝

将a1值代入式（5），得erfc（λ1）＝0.85，查erfc（λ）数值表，插值得λ2≈0.14，

表1 水位观测数据Table 1 Water table observed date m

3.2 配线法

根据观测孔K55.4的水位观测资料绘制实际资料曲线，与标准曲线叠置直至2条曲线重合（图4）；在横轴上任取一点作为匹配点A（图4中箭头所指）。读出该点的坐标值t＝14、1/λ2＝9.633，代入便可算出a＝528.0。将a代入式（5）求得tN＝76.2d，确定河岸无垂向入渗（或蒸发）地下水流非稳定流自由面位置计算方程，此处不再重述。

3.3 误差分析

分别用查表法、配线法及式（4）对观测孔K55.4的水位进行计算，计算结果与观测值对比见表2。

由表2可知：采用初始面非水平的解（查表法和配线法）较好地描述了含水层水位的变化规律，其计算结果较准确，水位计算值与观测值误差在0.02m以内；而采用初始面水平的解（公式（4））忽略了初始面非水平的影响，误差较大，水位计算值与观测值误差在0.0d时刻达到0.71m，随时间增大，误差逐渐减小，在60.0d为0.06m。

图4 叠置曲线Fig.4 Overlay curve

表2 水位计算值与观测值对比Table 2 Water table comparison of calculation value and observed value m

4 结论

1）把非水平初始面看作是最初初始面水平情况下，河流水位变化h1－h∞，0，历经时间tN，含水层中形成的自由面；结合初始条件与边界条件建立河流附近初始面非水平含水层非稳定渗流模型，提出模型的简化解及其适用条件；误差分析表明简化的解适用于河流水位变化初期、水位变幅相对于潜水含水层厚度较小、观测孔距河岸较近时的情况。

2）查表法和配线法的计算结果表明，参数a和tN的计算值不同，水位计算误差在0.02m以内，较原公式精度有明显提高。该法可用于缺少含水层参数、观测数据较少条件下的自由面计算。

3）在求解非水平初始面时，将n次水位变化过程简化为1次水位变化过程，便于求解，但在水位变化较大情况下会产生较大误差，有待进一步的研究。

（References）：

［1］Zlotnik V A，Huang Huihua.Effect of Shallow Penetration and Streambed Sediments on Aquifer Response to Stream Stage Fluctuations（Analytical Model）［J］.Ground Water，1999，37（4）：599－605.

［2］Moench A F，Barlow P M.Aquifer Response to Stream－Stage and Recharge Variations：I：Analytical Step－Response Functions［J］.Journal of Hydrology，2000，230（3）：192－210.

［3］Girard P，Silva C J，Abdo M.River－Groundwater Interactions in the Brazilian Pantanal：The Case of the CuiabáRiver［J］.Journal of Hydrology，2003，283（1）：57－66.

［4］O’Driscoll M，Johnson P，Mallinson D.Geological Controls and Effects of Floodplain Asymmetry on River－Groundwater Interactions in the Southeastern Coastal Plain， USA［J］.Hydrogeology Journal，2010，18（5）：1265－1279.

［5］Page R M，Lischeid G，Epting J，et al.Principal Component Analysis of Time Series for Identifying Indicator Variables for Riverine Groundwater Extraction Management［J］.Journal of Hydrology，2012，432：137－144.

［6］潘世兵，王忠静，邢卫国.河流含水层系统数值模拟方法探讨［J］.水文，2002，22（4）：19－21.Pan Shibing， Wang Zhongjing， Xing Weiguo.Discussion on Numerical Simulation for Stream－Aquifer System［J］.Hydrogeology，2002，22（4）：19－21.

［7］陶月赞，席道瑛.垂直与水平渗透作用下潜水非稳定渗流运动规律［J］.应用数学和力学，2006，27（1）：53－59.Tao Yuezan，Xi Daoying.Rule of Transient Phreatic Flow Subjected to Vertical and Horizontal Seepage［J］.Applied Mathematics and Mechanics，2006，27（1）：53－59.

［8］王文科，李俊亭，王钊，等.河流与地下水关系的演化及若干科学问题［J］.吉林大学学报：地球科学版，2007，37（2）：231－238.Wang Wenke， Li Junting， Wang Zhao，et al.Evolution of the Relationship Between River and Groundwater and Several Scientific Problems［J］.Journal of Jilin University：Earth Science Edition，2007，37（2）：231－238.

［9］Zissis T S，Teloglou I S，Terzidis G A.Response of a Sloping Aquifer to Constant Replenishment and to Stream Varying Water Level［J］.Journal of Hydrology，2001，243（3）：180－191.

［10］Singh S K.Explicit Estimation of Aquifer Diffusivity from Linear Stream Stage［J］.Journal of Hydraulic Engineering，2003，129（6）：463－469.

［11］Chang Yachi，Yeh H D.Analytical Solution for Groundwater Flow in an Anisotropic Sloping Aquifer with Arbitrarily Located Multiwells［J］.Journal of Hydrology，2007，347（1）：143－152.

［12］薛禹群.地下水动力学原理［M］.北京：地质出版社，1986.Xue Yuqun.Principle of Underground Water Dynamics［M］.Beijing： Geological Publishing House，1986.

［13］陶月赞，姚梅，张炳峰.时变垂向入渗影响下河渠－潜水非稳定渗流模型的解及应用［J］.应用数学和力学，2007，28（9）：1047－1053.Tao Yuezan，Yao Mei，Zhang Bingfeng.Solution and Its Application of Transient Stream/Groundwater Model Subjected to Time Dependent Vertical Seepage［J］.Applied Mathematics and Mechanics，2007，28（9）：1047－1053.

［14］陶月赞，曹彭强，席道瑛.垂向入渗与河渠边界影响下潜水非稳定流参数的求解［J］.水利学报，2006，37（8）：913－917.Tao Yuezan，Cao Pengqiang，Xi Daoying.Parameter Estimation for Semi－Infinite Phreatic Aquifer Subjected to Vertical Seepage and Bounded by Channel［J］.Journal of Hydraulic Engineering，2006，37（8）：913－917.