庫水波動下庫岸邊坡地下水浸潤線解析解及誤差分析

胡 慶，李 漪，李 恒，胡凱鋒，鄭德賓

(1.中國地震局地震研究所 地震大地測量重點(diǎn)實(shí)驗(yàn)室，武漢 430071;2.武漢地震工程研究院，武漢 430071;3.湖北省電力勘測設(shè)計(jì)院，武漢 430040)

1 研究背景

庫水波動是影響庫岸邊坡變形破壞最主要的動態(tài)因素，庫水波動下庫岸邊坡地下水浸潤線的合理計(jì)算對庫區(qū)滑坡穩(wěn)定性分析和預(yù)測預(yù)報(bào)都有重要意義。

近幾十年來，很多學(xué)者采用不同方法對邊/滑坡地下水滲流問題進(jìn)行了研究，如文獻(xiàn)［1－2］分別用有限元方法研究了降雨和庫水波動下邊坡地下水滲流場的變化趨勢;文獻(xiàn)［3］用有限差分法研究了降雨和庫水波動聯(lián)合作用下的滑坡地下水浸潤線。數(shù)值方法能較好解答復(fù)雜初始條件和邊界條件下的邊/滑坡地下水浸潤線問題，但目前還未在工程實(shí)際中廣泛使用，解析法是比數(shù)值法更容易推廣應(yīng)用的實(shí)用方法。文獻(xiàn)［4－6］基于含水層均質(zhì)、側(cè)向無限延伸并具水平不透水層，庫岸垂直，初始地下水位水平和潛水流為一維流等假定，用解析解法研究了不考慮降雨和蒸發(fā)作用時(shí)庫水位等速變化下的滑坡地下水浸潤線。文獻(xiàn)［7－10］在上述假定的基礎(chǔ)上研究了考慮降雨和蒸發(fā)影響時(shí)的情況，文獻(xiàn)［11］用解析解法研究了庫水位呈正弦半波曲線變化時(shí)的滑坡地下水浸潤線。上述解析解都是基于一定假定并對方程線性化處理后得到的，然而，用解析法計(jì)算的邊/滑坡地下水浸潤線經(jīng)常與實(shí)際情況不符，這主要是因?yàn)闆]有系統(tǒng)地研究各種假定和線性化過程帶來的誤差并分析其適用條件。

結(jié)合前人研究成果，本文基于一定假定建立了不同的庫岸邊坡地下水滲流數(shù)學(xué)模型，并求出了庫水位呈任意函數(shù)波動下的庫岸邊坡地下水浸潤線解析解。通過對比分析不同數(shù)學(xué)模型的數(shù)值解和解析解，研究了庫岸邊坡垂直處理、潛水運(yùn)動基本方程線性化處理、假定庫水位等速變化、不考慮非飽和滲流和Dupuit假定帶來的誤差及其規(guī)律，從而分析計(jì)算庫岸邊坡地下水浸潤線解析法的適用范圍。

2 解析解答

2.1 基本假定

由于潛水運(yùn)動基本方程(即布西涅斯克方程)是一個(gè)二階非齊次偏微分方程，只有通過一定的假定并對方程線性化處理后才能給出浸潤線的解析解。本文在建立用于解析解答的數(shù)學(xué)模型時(shí)，有以下基本假定:①含水層均質(zhì)、各向同性，側(cè)向無限延伸并具有水平隔水層;②潛水流為一維流，不考慮垂直潛水面方向的滲流(dupuit假設(shè));③ 庫水波動范圍的岸坡垂直;④不考慮非飽和滲流，即巖土介質(zhì)的滲透性與其含水量無關(guān);⑤由于降雨入滲和蒸發(fā)作用下巖土體中水的滲流過程十分復(fù)雜，目前的處理方法都有較大的誤差［12］，本文不考慮降雨入滲和蒸發(fā)情況。

2.2 計(jì)算模型

基于上述假設(shè)，可建立如圖1所示的幾何模型(將隔水底板高程定為0 m)。

圖1 幾何模型Fig.1 Geometric model

不考慮降雨入滲和蒸發(fā)作用時(shí)，潛水一維非穩(wěn)定運(yùn)動的基本方程為

式中:潛水面上升時(shí)μ為重力給水度，下降時(shí)μ為飽和差;t為時(shí)間;K為滲透系數(shù);x為位置坐標(biāo);h為含水層厚度;H為總水頭(隔水底板水平時(shí)，h和H相等)。式(1)為一非線性方程，須對其線性化后才能解答，通常采用的線性化方法是用平均含水層厚度hm近似代替h。設(shè)水位傳導(dǎo)系數(shù)a=Khm/μ，有

設(shè) s(x，t)=h(x，t)－ h(x，0)，可得到以含水層厚度變化函數(shù)s(x，t)為未知函數(shù)的數(shù)學(xué)模型:

式中f(t)為庫水位變動函數(shù)，可由庫水位－時(shí)間曲線經(jīng)過函數(shù)擬合得到。

2.3 解析解

對于上述數(shù)學(xué)模型，通過對t進(jìn)行拉普拉斯正變換，解答常微分方程，可得

式中:S(x，p)為 s(x，t)關(guān)于 t的拉普拉斯變換，F(xiàn)(p)為f(t)關(guān)于t的拉普拉斯變換。運(yùn)用拉普拉斯變換的微分性質(zhì)和卷積定理對上式進(jìn)行拉普拉斯逆變換，有

所以有:

若假定庫水位等速變化，即f(t)=vt，v為庫水變化速度，上升為正，下降為負(fù)。則

對上式(7)進(jìn)行拉普拉斯逆變換，有

由于M(0)=1，h(x，t)可寫成下列統(tǒng)一格式:

其中

M(z)為一特殊函數(shù)，可查表獲得對應(yīng)不同z值的函數(shù)值［13］，本文對M(z)函數(shù)用多項(xiàng)式擬合成如下公式:

具體解答時(shí)，平均含水層厚度hm涉及到滲流結(jié)束時(shí)的含水層厚度，需迭代才能確定。迭代初值取為h0+1/3H'，h0為初始狀態(tài)的平均含水層厚度，H'為庫水上升或下降的幅度，上升為正，下降為負(fù)。用初值代入浸潤線計(jì)算函數(shù)，可得滲流結(jié)束時(shí)的平均含水層厚度hmn，當(dāng)hmn與上一步計(jì)算結(jié)果hmn－1相差足夠小時(shí)，即可停止迭代計(jì)算。

3 誤差分析

針對上述建立數(shù)學(xué)模型時(shí)的諸多假定和解答過程中的線性化方法，建立不同的“參照”模型，用有限元程序解答上述“參照”模型并與對應(yīng)的解析解做比較分析，就可得到各種誤差的大小和規(guī)律。本文用表示浸潤線計(jì)算誤差，式中(如圖2所示):s1為x軸、“參照”模型的浸潤線、過2條浸潤線交點(diǎn)且平行于h軸的直線三者圍成的面積;s2為2條浸潤線圍成的面積。本文以三峽庫區(qū)巴東新縣城東部的趙樹嶺滑坡為例進(jìn)行研究。趙樹嶺滑坡為一基巖順層滑坡，滑體南北長約1 000 m，東西寬平均350 m，滑體平均厚度約35 m?；w主要由塊裂、碎裂巖、含泥碎塊石及碎塊石組成;基巖則以T2b2紫紅色粉砂質(zhì)泥巖、泥質(zhì)粉砂巖為主。

圖2 岸坡垂直處理誤差分析Fig.2 Analysis of errors caused by assuming that the reservoir bank is vertical

3.1 岸坡垂直處理帶來的誤差

解析解答中建立滲流數(shù)學(xué)模型時(shí)假定庫水波動范圍內(nèi)的岸坡(圖1中的A'B段)坡度為90°。分別建立對應(yīng)圖1中實(shí)際剖面和假定剖面的有限元數(shù)值模型，取平均滲透系數(shù)為0.3 m/d，庫水位自34 m高程以1 m/d速度經(jīng)30 d上升至64 m時(shí)的地下水浸潤線計(jì)算結(jié)果如圖2所示。

2種數(shù)值模型的浸潤線計(jì)算結(jié)果較為接近，P=6.9%。實(shí)際岸坡剖面的浸潤線要比垂直處理后岸坡的浸潤線計(jì)算結(jié)果更陡一些。如圖1所示，岸坡剖面的垂直化處理相當(dāng)于“削”掉了一部分岸坡巖土體，使得庫水在補(bǔ)給岸坡地下水時(shí)，其滲流路徑更短，這就導(dǎo)致岸坡內(nèi)地下水浸潤線隨庫水位變化的響應(yīng)更快。因此，庫水上升時(shí)，垂直處理后岸坡的浸潤線形狀相對緩一些。當(dāng)平均滲透系數(shù)取值變大，2種模型計(jì)算的浸潤線都會越來越接近水平，這時(shí)P將逐漸變小至0。

3.2 方程線性化處理帶來的誤差

有限元法不必對方程線性化處理即可得到浸潤線計(jì)算結(jié)果，因此，通過數(shù)值解和解析解的比較就可得到方程線性化處理過程帶來的誤差。建立k=0.3 m/d，初始水位為20 m高程并以1 m/d速度上升的地下水滲流數(shù)學(xué)模型。不同庫水位波動幅度下的數(shù)值解和解析解如圖3所示，P值如表1所示?？梢钥闯?，隨著庫水位波動幅度ΔH'的增大，解析解的誤差越來越大。若以P=15%為容許誤差值，解析解法只適用于ΔH'/hm小于30%的情況。

3.3 假定庫水位等速變化帶來的誤差

本文雖然給出了庫水位呈任意函數(shù)變化時(shí)的庫

圖3 不同庫水位波動幅度下的數(shù)值解和解析解Fig.3 Arithmetic and analytic solutions of the equation with different water-level fluctuations

表1 不同庫水上升幅度下的P值Table 1 Values of P in the presence of different water level rise

岸邊坡浸潤線計(jì)算公式，但是其結(jié)果為函數(shù)的變限積分形式，由于實(shí)際庫水波動曲線比較復(fù)雜，很多情況下不能求得便于計(jì)算的表達(dá)式。因此，常假定庫水位呈等速變化，對于實(shí)際水位變化復(fù)雜的情況可分段考慮，但這就要求每一段都要有一個(gè)初始地下水位;因而在實(shí)際運(yùn)用中往往僅考慮水位變化的初始值和末尾值而將其簡化為等速變化。經(jīng)過不同庫水變化函數(shù)下浸潤線計(jì)算結(jié)果的比較，發(fā)現(xiàn)假定庫水位等速變化帶來的誤差有以下特點(diǎn):當(dāng)最后時(shí)間段庫水實(shí)際變化速度和假定的平均變化速度相差越大時(shí)，誤差越大;當(dāng)滲透系數(shù)取值變大時(shí)，2種模型計(jì)算的浸潤線都會越來越接近水平，同時(shí)P將逐漸變小至0。

k=0.3 m/d時(shí)，圖4中2種庫水位變化曲線下的浸潤線計(jì)算結(jié)果如圖5所示。

對于三峽水庫而言，水位變化速度一般不超過1 m/d，且快速蓄水和泄水時(shí)，其變化曲線與直線很接近，所以，假定庫水位等速變化僅對三峽庫區(qū)中滲透性很差的庫岸邊坡浸潤線計(jì)算有影響。

3.4 不考慮非飽和滲流帶來的誤差

如果考慮非飽和滲流，則潛水非穩(wěn)定運(yùn)動方程中不僅含水層厚度是一個(gè)變量，而且滲透系數(shù)也是一個(gè)變量，其方程如下:

圖4 庫水位變化曲線Fig.4 Curves of water level fluctuation

圖5 庫水位變化函數(shù)簡化誤差分析Fig.5 Errors caused by assuming that the changes in water level is uniform

這是一個(gè)非線性化程度更高的方程，解答起來將更加困難。因此，目前的解析解答都沒有考慮非飽和滲流，常用一個(gè)平均的滲透系數(shù)k代替非飽和滲透函數(shù)k(θ)。在有限元數(shù)值模型中可以用非飽和滲透函數(shù)K(θ)=k(k為常數(shù))代表解析解答中不考慮非飽和滲流的情況。通過有限元數(shù)值方法解答上述數(shù)值模型和實(shí)際情況(其非飽和滲透函數(shù)按文獻(xiàn)［5］取值)下的數(shù)值模型。庫水位自34 m高程以1 m/d速度經(jīng)5 d上升至39 m高程時(shí)的浸潤線計(jì)算結(jié)果如圖6所示，P=17.3%。通過不同的滲透系數(shù)取值可以發(fā)現(xiàn)，滲透性越好，其誤差越小。需要說明的是，不考慮非飽和滲流所引起的誤差很大程度上受控于滲流介質(zhì)的非飽和滲透函數(shù)，顯然，滲透系數(shù)在不同含水量狀態(tài)下差值越大，誤差越大。

圖6 不考慮非飽和滲流的誤差分析Fig.6 Analysis of errors caused by the neglect of unsaturated seepage

3.5 Dupuit假設(shè)帶來的誤差

目前廣泛應(yīng)用的非穩(wěn)定滲流基本方程是基于Dupuit假設(shè)的，即假定地下水不發(fā)生垂直于浸潤線(面)方向的滲流。這種假定在滲出面等水頭變化較大的地方附近有一定誤差，因?yàn)檫@些地方垂直方向的滲流速度較大。其誤差大約為［14］

4 結(jié)論

(1)由于潛水非穩(wěn)定滲流基本方程為二階非線性偏微分方程，須基于一定假定并對方程線性化處理后才能求得庫岸邊坡地下水浸潤線的解析解，這些假定和線性化處理過程必然帶來一定的誤差。

(2)庫岸垂直處理帶來的誤差要小于方程線性化處理和不考慮非飽和滲流帶來的誤差，而且它們都隨滲透系數(shù)變大而減小;方程的線性化處理帶來的誤差隨庫水位波動幅度與含水層厚度之比增大而增大;在三峽庫區(qū)庫水調(diào)度情況下，假定庫水位等速變化帶來的誤差可以忽略。

(3)求解庫岸邊坡地下水浸潤線時(shí)，解析法有一定的適用范圍，針對不同的情況要具體分析。一般而言，解析法適用于庫岸邊坡巖土結(jié)構(gòu)簡單、幾何邊界規(guī)則、滲透性較好且?guī)焖▌臃刃∮?0%含水層厚度時(shí)的情況。

［1］吳宏偉，陳守義，龐宇威.雨水入滲對非飽和土坡穩(wěn)定性影響的參數(shù)研究［J］.巖土力學(xué)，1999，3(1):1 －14.(WU Hong-wei， CHEN Shou-yi， PANG Yu-wei.Parametric Study of Effects of Rain Infiltration on Unsaturated Slopes［J］.Rock and Soil Mechanics，1999，3(1):1 － 14.(in Chinese))

［2］章廣成，唐輝明，胡 斌.非飽和滲流對滑坡穩(wěn)定性的影響研究［J］.巖土力學(xué)，2007，28(5):965 －970.(ZHANG Guang-cheng，TANG Hui-ming，HU Bin.Study of Influence of Unsaturated Seepage on Stability of Landslide［J］.Rock and Soil Mechanics，2007，28(5):965 －970.(in Chinese))

［3］李 曉，張年學(xué)，廖秋林，等.庫水位漲落與降雨聯(lián)合作用下滑坡地下水動力場分析［J］.巖石力學(xué)與工程學(xué)報(bào)，2004，23(21):3714 － 3720.(LI Xiao，ZHANG Nianxue，LIAO Qiu-lin，et al.Analysis of Hydrodynamic Field Influenced by Combination of Rainfall and Reservoir Level Fluctuation［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23(21):3714 －3720.(in Chinese))

［4］鄭穎人，時(shí)衛(wèi)民，孔位學(xué).庫水位下降時(shí)滲透力及地下水浸潤線的計(jì)算［J］.巖石力學(xué)與工程學(xué)報(bào)，2004，23(18):3203 － 3210.(ZHENG Ying-ren，SHI Wei-min，KONG Wei-xue.Calculation of Seepage Forces and Phreatic Surface under Drawdown Conditions［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23(18):3203 －3210.(in Chinese))

［5］馬崇武，劉忠玉，苗天德，等.江河水位升降對堤岸邊坡穩(wěn)定性的影響［J］.蘭州大學(xué)學(xué)報(bào)(自然科學(xué)版)，2000，36(3):56 － 60.(MA Chong-wu，LIU Zhong-yu，MIAO Tian-de，et al.The Influence of Water Level Changing on the Stability of River Embankment［J］.Journal of Lanzhou University(Natural Sciences)，2000，36(3):56 －60.(in Chinese))

［6］楊紅坡，謝新宇，張繼發(fā)，等.潛水一維非穩(wěn)態(tài)運(yùn)動的解析理論及應(yīng)用［J］.水科學(xué)進(jìn)展，2004，15(1):82－86.(YANG Hong-po，XIE Xin-yu，ZHANG Ji-fa，et al.Analytical Solution of One-dimensional Transient Phreatic Flow and Its Application［J］.Advances in Water Science，2004，15(1):82 －86.(in Chinese))

［7］WARRICK A W，WIERENGA P J，PAN L.Downward Water Flow Through Sloping Layers in the Vadose Zone:Analytical Solutions for Diversions［J］.Journal of Hydrology，1997，192:321 －337.

［8］馮文凱，石豫川，柴賀軍，等.降雨及庫水升降作用下地下水浸潤線簡化求解［J］.成都理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，2006，33(1):90－94.(FENG Wen-kai，SHI Yuchuan，CHAI He-jun，et al.The Simplified Solution of Phreatic Saturation Line under the Actions of Rainfall and Reservoir Water Level Fluctuation［J］.Journal of Chengdu University of Technology(Science＆ Technology Edition)，2006，33(1):90 －94.(in Chinese))

［9］張友誼，胡卸文.庫水位等速上升作用下岸坡地下水浸潤線的計(jì)算［J］.水文地質(zhì)工程地質(zhì)，2007，(5):46 －49.(ZHANG You-yi，HU Xie-wen.Calculation of Saturation Line of Groundwater under Reservoir Water Table Uniform Rising［J］.Hydrogeology and Engineering Geology，2007，(5):46 －49.(in Chinese))

［10］吳 瓊，唐輝明，王亮清，等.庫水位升降聯(lián)合降雨作用下庫岸邊坡中的浸潤線研究［J］.巖土力學(xué)，2009，30(10):3025 －3031.(WU Qiong，TANG Hui-ming，WANG Liang-qing，et al.Analytic Solutions for Phreatic Line in Reservoir Slope with Inclined Impervious Bed under Rainfall and Reservoir Water Level Fluctuation［J］.Rock and Soil Mechanics，2009，30(10):3025 － 3031.(in Chinese))

［11］李相依，晏鄂川.庫水位周期性升降作用下滑坡堆積體潛水位計(jì)算［J］.四川大學(xué)學(xué)報(bào)(工程科學(xué)版)，2007，39(增刊):78 － 81.(LI Xiang-yi，YAN E-chuan.Calculation of Phreatic Surface of Bank Talus under Periodic Fluctuating Reservoir Water Level［J］.Journal of Sichuan University(Engineering Science Edition)，2007，39(Sup.):78 －81.(in Chinese))

［12］朱 偉，陳學(xué)東，鐘小春.降雨入滲規(guī)律的實(shí)測與分析［J］.巖土力學(xué)，2006，27(11):1873 － 1879.(ZHU Wei，CHEN Xue-dong，ZHONG Xiao-chun.Observation and A-nalysis of Rainfall Infiltration［J］.Rock and Soil Mechanics，2006，27(11):1873 －1879.(in Chinese))

［13］張蔚榛.地下水非穩(wěn)定滲流計(jì)算和地下水資源評價(jià)［M］.北京:科學(xué)出版社，1983.(ZHANG Wei-zhen.Calculation of Unsteady Groundwater Seepage and Evaluation of Groundwater Resources［M］.Beijing:Science Press，1983.(in Chinese))

［14］薛禹群.地下水動力學(xué)［M］.北京:地質(zhì)出版社，1997.(XUE Yu-qun.Groundwater Dynamics［M］.Beijing:Geological Press，1997.(in Chinese))