曹彭強(qiáng)，束龍倉，魯程鵬

河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實驗室，南京 210098

0 前言

河流－地下水系統(tǒng)事關(guān)河流與含水層之間水量交換評價、農(nóng)田灌排水渠設(shè)計、河道基流衰減過程研究、潛水含水層參數(shù)評估、河流與地下水水質(zhì)交換等方面的問題［1－5］，是地下水滲流力學(xué)中的研究熱點(diǎn)和難點(diǎn)之一。從20世紀(jì)40年代開始，學(xué)者們［6－8］就河流－地下水系統(tǒng)提出了不同的模型，并對模型進(jìn)行了求解，就河流入滲對潛水的影響進(jìn)行了分析?，F(xiàn)有文獻(xiàn)［9－10］主要針對隔水底板水平的條件下假設(shè)初始自由面水平，或隔水底板傾斜條件下初始自由面平行于傾斜隔水底板，這些假設(shè)簡化了描述初始面的數(shù)學(xué)方程，從而方便模型的求解；但在自然狀態(tài)下，初始自由面非水平更具普遍性，因此，如何描述與計算非水平初始面成為滲流計算的關(guān)鍵。

筆者以河流附近具水平隔水底板的潛水含水層中的潛水滲流模型為例，構(gòu)造初始面非水平條件下的一維潛水非穩(wěn)定滲流模型；利用Boussinesq第一線性化方法，通過Laplace變換，再結(jié)合迭加原理，提出非水平初始面和參數(shù)的簡易計算方法，導(dǎo)出模型的解，并根據(jù)誤差分析討論解的適用條件，為初始面非水平條件下的潛水滲流計算提供方法。

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

假設(shè)一順直河流完整切割含水層（圖1），潛水含水層均質(zhì)各向同性，具有水平隔水底板，側(cè)向無限延伸；潛水初始自由面非水平hx，0；河流初始水位為h1，無限遠(yuǎn)處含水層水位為h∞，0；河流水位瞬時上升至h2，水位升幅為h2－h(huán)1，且上升后保持不變；潛水流可視為一維流。

上述水文地質(zhì)概念模型是在經(jīng)典的Ferris模型［11］基礎(chǔ)上，增加了面平均入滲或蒸發(fā)條件，且假設(shè)潛水初始自由面hx，0非水平。把hx，0看作原河流水位及水平面水位都為h∞，0（稱為最初初始面），河流水位變化h1－h(huán)∞，0，歷經(jīng)時間tN，含水層中形成的自由面。河流附近初始面非水平的半無限潛水含水層滲流控制方程為

邊界條件及初始條件為

圖1 河流附近潛水含水層Fig.1 Phreatic aquifer near river

其中：μ為給水度；K為滲透系數(shù)（m/d）；h為水位（m）；w為入滲強(qiáng)度（m/d，補(bǔ)給為正，排泄為負(fù)）；x為距河邊距離（m）；t為河流水位波動時間（d）；tN為初始面形成時間（d）。

控制方程與邊界條件及初始條件構(gòu)成了河流附近初始面非水平的半無限潛水含水層滲流數(shù)學(xué)模型。其計算坐標(biāo)如圖1所示（對于寬深比較小的河流，將坐標(biāo)原點(diǎn)定于水位上升后的水土交界處）。初始自由面hx，0非水平，模型（1）難以直接進(jìn)行線性化求解。

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

2.1 數(shù)學(xué)模型求解

模型（1）求解的關(guān)鍵點(diǎn)及難點(diǎn)在于初始自由面hx，0的確定。形成初始自由面hx，0的過程是個滲流過程，滿足地下水流非穩(wěn)定流規(guī)律。

河流水位常有一定漲落，呈階梯變化或連續(xù)變化，為了便于計算，把河流水位變化曲線概化為階梯狀線段［12］。當(dāng)水位變化相對于含水層厚度較小時，將水位連續(xù)變化過程看作數(shù)個瞬時變化的疊加［13］，參照初始自由面水平條件下河流附近潛水含水層地下水流非穩(wěn)定流計算公式求解方法［14］，當(dāng)hx，t－h(huán)x，0≤0.1hm（hm為潛水含水層的平均厚度）時，利用Boussinesq第一線性化方法，通過Laplace變換，得河流水位連續(xù)變化情況下河流附近的初始面計算公式：

其中：Δhi為第i次水位變化ti－1為第i－1次水位變化后的時刻；t0＝0；a＝Khm/μ，為導(dǎo)壓系數(shù)；erf（）和erfc（）分別為誤差函數(shù)和余誤差函數(shù)。

河流水位瞬時上升Δh2并保持不變，由迭加原理，含水層中水位可由下式計算：

當(dāng)h1＝h∞，0，即初始面水平時，Δhi＝0，式（3）變?yōu)?/p>

式（4）為初始面水平潛水含水層非穩(wěn)定流計算公式。

式（2）中的hi與ti－1在河流水位變化過程中不易取得，難以計算，若將n次水位變化過程簡化為1次水位變化過程，即式（2）變?yōu)?/p>

式（3）變?yōu)?/p>

式（6）即為簡化后的初始面非水平潛水含水層非穩(wěn)定流計算公式。

2.2 參數(shù)的確定

應(yīng)用式（6）計算初始面非水平含水層水位時，需獲得參數(shù)x、t、a和tN的取值。其中，x和t是確定的，對于某一次水位的變化，還需要確定參數(shù)a和tN。a和tN的具體求解過程有2種方法：查表法和配線法。

2.2.1 查表法

現(xiàn)已知河流初始水位h1、河流水位上升后水位h2、距離河流無窮遠(yuǎn)處含水層水位h∞，0。h∞，0即是離河流無窮遠(yuǎn)處的潛水水位，但在實際應(yīng)用中只需取足夠遠(yuǎn)、含水層水位影響較小處的水位即可。

式（6）減式（5）得

因此，式（7）簡化為

即

式（9）中的未知數(shù)a，可選取離河流較近觀測孔中的數(shù)據(jù)hx，t，查erfc（）數(shù)值表求得。再將a代入式（5），求得tN。

查表法計算過程簡單，所需水位資料較少，只需要河流水位變化初期（t→0）的一個離河流較遠(yuǎn)的水位觀測孔和一個離河流較近的觀測孔的水位數(shù)據(jù)即可，適用于河岸地下水水位觀測資料不足的情況。但由于erfc（）數(shù)值表中數(shù)據(jù)有限，在查表求解a時需進(jìn)行插值計算。

2.2.2 配線法

令

對式（9）和式（10）同時取對數(shù)。對于同一個觀測孔而言，lg（x2/4a）為一常數(shù)。則lg［（hx，t－h(huán)x，0）/（h2－h(huán)1）］－lg（t）曲線（實際資料曲線）和lg［erfc（λ）］－lg（1/λ2）曲線（標(biāo)準(zhǔn)曲線）的圖形是相同的，只是橫坐標(biāo)平移了x2/4a而已。如果已知坐標(biāo)值及水位觀測值，就可以計算參數(shù)a。

具體步驟為：根據(jù)erfc（λ）數(shù)值表繪制lg［erfc（λ）］－lg（1/λ2）標(biāo)準(zhǔn)曲線（圖2）；根據(jù)觀測孔的實際水位觀測數(shù)據(jù)在另一張有相同坐標(biāo)的透明紙上繪制lg［（hx，t－h(huán)x，0）/（h2－h(huán)1）］－lg（t）實際資料曲線；將實際資料曲線疊置在標(biāo)準(zhǔn)曲線上，在保持對應(yīng)坐標(biāo)軸平行的條件下，移動坐標(biāo)紙，直至2條曲線重合為止；之后在圖上任取一點(diǎn)作為匹配點(diǎn)。讀出該點(diǎn)的坐標(biāo)值t、1/λ2，則

再將a代入式（5），計算不同點(diǎn)、不同時刻潛水水位。配線法對于只有河流水位未發(fā)生變化前觀測孔的潛水水位數(shù)據(jù)的情況也適用，但在進(jìn)行配線時會受個人主觀因素影響，對于觀測孔水位資料序列較短的情況，不同人給出的配線數(shù)值差別較大。

圖2 lg［erfc（λ）］－lg（1/λ2）標(biāo)準(zhǔn)曲線Fig.2 lg［erfc（λ）］－lg（1/λ2）standard curve

2.3 簡化后解的誤差分析

式（6）減式（3）得x處潛水初始面計算誤差ε為

其中：ε是x，a，tN，t，ti－1的 函 數(shù)，ε＝f（x，a，tN，t，ti－1）；0≤ti－1≤tN。當(dāng)Δhi＞0時，ε≤0，即用1次水位變化過程代替多次水位變化過程計算所得的初始潛水水位小于實際觀測值。

式（12）中ε分別對x，a，tN，t，ti－1求導(dǎo)，得ε與各參量的關(guān)系（圖3）。由圖3可知：隨x增大，ε先增后減，在x＝0和x＝∞處，ε＝0，即在離河流較近距離和足夠遠(yuǎn)處，采用式（6）計算的誤差較小；隨a增大，ε先增后減，當(dāng)a較小時，ε較小，也即當(dāng)含水層K和hm較小、μ較大時，計算的誤差較??；隨tN增大，ε先增后減，當(dāng)tN＝0時，ε＝0；隨t增大，ε先增后減，即水位變化初期和末期，計算誤差較小；ε隨ti－1增大而增大，形成初始面的時間越長，后續(xù)水位變化計算的水位誤差越大。

由以上分析可得，在河流水位變化初期，水位變幅相對于潛水含水層厚度較小、觀測孔距河岸較近，即t較小、Δhi與x較小時，計算的ε較小，才能用式（6）代替式（3）進(jìn)行計算。

3 算例

圖3 ε與參數(shù)關(guān)系Fig.3 Relationship ofεand parameters

某河流在某時刻的水位從23.00m上升為25.00m，并保持不變。垂直河岸布置了2個觀測孔（圖1），觀測水位資料見表1，含水層入滲w為0。

3.1 查表法

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

由表1可知，河流水位發(fā)生變化后，離河岸較遠(yuǎn)的觀測孔（K446.4）的水位變化極小，可以把該點(diǎn)t＝0時刻的水位作為最初初始水位，即H∞，0＝22.16m。

選取t值較?。?.1d）時、離河較近的觀測孔（K55.4）的水位數(shù)據(jù)代入式（5），即x＝55.4m，t＝1.1d，hx，0＝22.87m，hx，t＝23.02m，得erfc（λ1）＝7.50×10－3。

查erfc（λ1）數(shù)值表，插值得λ1≈1.26，則a1＝

將a1值代入式（5），得erfc（λ1）＝0.85，查erfc（λ）數(shù)值表，插值得λ2≈0.14，

表1 水位觀測數(shù)據(jù)Table 1 Water table observed date m

3.2 配線法

根據(jù)觀測孔K55.4的水位觀測資料繪制實際資料曲線，與標(biāo)準(zhǔn)曲線疊置直至2條曲線重合（圖4）；在橫軸上任取一點(diǎn)作為匹配點(diǎn)A（圖4中箭頭所指）。讀出該點(diǎn)的坐標(biāo)值t＝14、1/λ2＝9.633，代入便可算出a＝528.0。將a代入式（5）求得tN＝76.2d，確定河岸無垂向入滲（或蒸發(fā)）地下水流非穩(wěn)定流自由面位置計算方程，此處不再重述。

3.3 誤差分析

分別用查表法、配線法及式（4）對觀測孔K55.4的水位進(jìn)行計算，計算結(jié)果與觀測值對比見表2。

由表2可知：采用初始面非水平的解（查表法和配線法）較好地描述了含水層水位的變化規(guī)律，其計算結(jié)果較準(zhǔn)確，水位計算值與觀測值誤差在0.02m以內(nèi)；而采用初始面水平的解（公式（4））忽略了初始面非水平的影響，誤差較大，水位計算值與觀測值誤差在0.0d時刻達(dá)到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 結(jié)論

1）把非水平初始面看作是最初初始面水平情況下，河流水位變化h1－h(huán)∞，0，歷經(jīng)時間tN，含水層中形成的自由面；結(jié)合初始條件與邊界條件建立河流附近初始面非水平含水層非穩(wěn)定滲流模型，提出模型的簡化解及其適用條件；誤差分析表明簡化的解適用于河流水位變化初期、水位變幅相對于潛水含水層厚度較小、觀測孔距河岸較近時的情況。

2）查表法和配線法的計算結(jié)果表明，參數(shù)a和tN的計算值不同，水位計算誤差在0.02m以內(nèi)，較原公式精度有明顯提高。該法可用于缺少含水層參數(shù)、觀測數(shù)據(jù)較少條件下的自由面計算。

3）在求解非水平初始面時，將n次水位變化過程簡化為1次水位變化過程，便于求解，但在水位變化較大情況下會產(chǎn)生較大誤差，有待進(jìn)一步的研究。

（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］潘世兵，王忠靜，邢衛(wèi)國.河流含水層系統(tǒng)數(shù)值模擬方法探討［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］陶月贊，席道瑛.垂直與水平滲透作用下潛水非穩(wěn)定滲流運(yùn)動規(guī)律［J］.應(yīng)用數(shù)學(xué)和力學(xué)，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］王文科，李俊亭，王釗，等.河流與地下水關(guān)系的演化及若干科學(xué)問題［J］.吉林大學(xué)學(xué)報：地球科學(xué)版，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］薛禹群.地下水動力學(xué)原理［M］.北京：地質(zhì)出版社，1986.Xue Yuqun.Principle of Underground Water Dynamics［M］.Beijing： Geological Publishing House，1986.

［13］陶月贊，姚梅，張炳峰.時變垂向入滲影響下河渠－潛水非穩(wěn)定滲流模型的解及應(yīng)用［J］.應(yīng)用數(shù)學(xué)和力學(xué)，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］陶月贊，曹彭強(qiáng)，席道瑛.垂向入滲與河渠邊界影響下潛水非穩(wěn)定流參數(shù)的求解［J］.水利學(xué)報，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.