河流附近非水平初始面潛水滲流分析

曹彭強，束龍倉，魯程鵬

河海大學水文水資源與水利工程科學國家重點實驗室，南京 210098

0 前言

河流－地下水系統事關河流與含水層之間水量交換評價、農田灌排水渠設計、河道基流衰減過程研究、潛水含水層參數評估、河流與地下水水質交換等方面的問題［1－5］，是地下水滲流力學中的研究熱點和難點之一。從20世紀40年代開始，學者們［6－8］就河流－地下水系統提出了不同的模型，并對模型進行了求解，就河流入滲對潛水的影響進行了分析?，F有文獻［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）計算的誤差較?。浑Sa增大，ε先增后減，當a較小時，ε較小，也即當含水層K和hm較小、μ較大時，計算的誤差較?。浑StN增大，ε先增后減，當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.