基于FEFLOW的岸坡地下水三維滲流場模擬研究

李慶偉，晏鄂川，楊 廣，崔學杰

(1.中國地質大學(武漢)工程學院，湖北 武漢 430074；2.中交第二航務工程勘查設計院有限公司，湖北 武漢 430060)

三峽庫水位的周期性波動引起岸坡內地下水滲流場不斷地發生相應的變化，使得坡體與水的相互作用在很大程度上得以加強，導致岸坡穩定性的改變，甚至可能引發岸坡的局部或整體變形破壞。國內外許多學者對岸坡地下水滲流特征進行了一系列研究，取得了豐碩的科研成果。如Sevaldson[1]、Henkel[2]、Liu[3]和Wang等[4]對庫水位升降作用下，岸坡地下水滲流場和失穩機理進行了相關研究，認為水是導致坡體變形破壞的主控因素；Szabo等[5]利用有限差分模型,模擬了瞬時水流的非穩定滲流；榮冠等[6]在研究降雨入滲機理和模擬方法的基礎上，編寫了非飽和滲流程序SUSC，并模擬了某邊坡降雨過程中地下水滲流場的變化；謝瑾榮等[7]基于軟巖邊坡非飽和滲流-應力計算原理，將滲流與軟化相結合，采用間接耦合法分析軟巖邊坡的降雨入滲，建立了軟巖邊坡滲流效應的分析方法；楊廣等[8]利用數值模擬方法研究了地下水滲流和岸坡的穩定性，認為岸坡不同位置和高程處地形的改變因地表水力特征的差異而有所不同；還有一些學者也進行了大量的相關研究[9-12]。但由于不同庫岸段工程地質條件的差異，部分特殊岸坡地下水的滲流特征仍沒有得到合理的解釋。

基于以上認識，本文在三峽水庫蓄水條件下，通過構建涉水岸坡的滲流模型，基于FEFLOW軟件對岸坡地下水滲流場的變化情況進行了三維數值模擬研究，從而揭示水庫正常運行期間庫岸內部地下水滲流場的發展動向與變化趨勢。

1 涉水岸坡工程地質概況

圖1 研究區地層分布圖Fig.1 Formation distribution map

研究區位于重慶巫山縣長江與大寧河交匯處江東咀地區，平面上呈“舌狀”，為長江與大寧河的階地，屬于構造侵蝕、剝蝕河谷地貌。據現場調查和鉆孔揭露，研究區內地層主要為第四系人工填土層、沖洪積層、崩坡積層和三疊系下統嘉陵江組第四段灰巖地層。第四系覆蓋物以粉質黏土夾碎石、碎塊石土、砂卵石層為主。其中，粉質黏土夾碎石在研究區分布較為廣泛[見圖1(a)]，主要位于大和尚包、小和尚包及小狗架子溝附近，厚度分布不均勻，最厚可達70 m；碎塊石土主要分布在紅梁子溝—小狗架子溝以北，厚度一般為50～60 m，局部地區可達到80 m[見圖1(b)]；紅梁子溝—小狗架子溝以南的部分鉆孔及以北的鉆孔SXZK10中揭露到砂卵石夾碎石層[圖1(c)]，砂卵石層一般直接覆蓋在基巖之上，上覆碎塊石土或粉質黏土夾碎石,其厚度變化較大，最厚處達5.2 m。基巖以三疊系下統嘉陵江組第四段灰巖為主，主要在研究區東北部出露地表，即神女廟周邊地區，基巖的埋深變化較大，基巖面陡峭[見圖1(d)]，基巖頂面高程為90～240 m,但在紅梁子溝附近基巖頂面高程急劇降低(90～100 m)，推測是巖溶塌陷造成；此外，基巖在神女廟附近埋深較淺，基巖頂面高程為210～240 m，覆蓋物厚度為10～20 m。研究區內及其附近無斷層通過。

2 計算原理與方法

2.1 計算原理

當不考慮水的密度的變化條件下，在孔隙介質中地下水在三維空間的流動可以用下面的偏微分方程來表示[13]：

(1)

式中：Kxx、Kyy、Kzz分別為x、y、z主方向的滲透系數；ω為源匯項，包括蒸發、降雨入滲補給、井的抽水量和泉的排泄量等；μs為儲水率；H0為初始地下水水位；Ω為地下水滲流區域；H1為河流水位；S1為第一類邊界(給定水位邊界條件)；S2為第二類邊界(給定流量邊界條件)；q(x,y,z,t)表示在邊界不同位置上不同時間的流量。

公式(1)構成了對于一個實際問題地下水流動的定解問題，可采用數值計算方法進行求解，如有限差分法、有限單元法等，求解結果即為地下水水位的分布值。

2.2 計算方法

目前地下水滲流數值計算主要采用有限單元法和有限差分法。其中，美國地質調查局在20世紀80年代開發的MODFLOW(Modular three-dimensional finite difference groundwater flow model)為有限差分法的典型代表[14]，而本文研究所使用的FEFLOW(Finite Element subsurface FLOW system)則是基于有限元法的地下水數值模擬軟件。該軟件是一種求解偏微分方程定解問題的有效數值方法，成功地解決了眾多領域的許多計算問題[15-17]，如地下水流動、水動力彌散問題等。

3 數值模擬方案設計

在現場調查的基礎上，結合研究區的地質資料，建立了水庫正常運行條件下岸坡地下水運動的三維地質模型，并利用數值計算方法對不同工況下岸坡地下水三維滲流場進行模擬。依據三峽工程正常運行期水庫的調度特征(見圖2)，本次數值模擬研究采用如下4種工況：①145 m低庫水位(即工況Ⅱ)；②145 m庫水位上升至175 m(即工況Ⅱ)；③175 m高庫水位(即工況Ⅲ)；④175 m庫水位下降至145 m(即工況Ⅳ)。為了簡化模擬過程，每月按30 d計，對較小的庫水位波動忽略不計，以該期間的特征庫水位為準，設計的涉水岸坡地下水三維滲流場數值模擬具體方案，見表1。

圖2 三峽工程正常運行期水庫調度圖Fig.2 Reservoir operation chart of Three Gorges Project during normal operation period注：(1)為防破壞線；(2)為限制供水線；(Ⅰ)為防洪區；(Ⅱ)為裝機預想出力區；(Ⅲ)為保證出力區；(Ⅳ)為降低出力區(資料來源于《三峽庫區地質災害防治工程地質勘查技術要求》)

工況時間庫水位/m歷時/dⅠ6月中旬至9月145105Ⅱ10月145～17530Ⅲ11月至次年4月末175180Ⅳ5月至6月中旬175～14545

3.1 數值模型的構建

數值計算范圍根據江東咀庫岸段研究范圍和庫水位變化特征綜合確定。根據三峽庫區的庫水位調度特征，本次岸坡地下水三維滲流場的數值模擬庫水位變動范圍設定為145～175 m，選取130 m等高線在水平面上的投影作為模型邊界，頂面高程則依據實際地形確定。為了便于建模，將地形整體順時針旋轉了43°，模型在平面上的投影近似于一個梯形，見圖3。數值計算模型長約726 m，寬約592 m，平面面積約為0.42 km2，地表高程范圍為130～268 m，相對高差為138 m。

圖3 數值模型范圍Fig.3 Domain of numerical model注：zk22〇為鉆孔編號及位置

研究區內土層的分布情況較為復雜，且存在尖滅情況，完全再現地質結構是很難實現的。另外，研究區粉質黏土夾碎石與碎塊石土滲透系數的差異遠小于庫水位的日升降幅度，因此在進行岸坡地下水三維滲流場數值模擬時可將其視為同一層；砂卵石層分布范圍小且厚度薄，對岸坡地下水滲流場的影響十分有限，因此本次模擬忽略砂卵石層，將研究區地層簡化為上覆第四系土層和下伏基巖的結構。

建立的舌狀涉水岸坡三維滲流模型，見圖4。該模型分為兩層，即第四系覆蓋層和基巖。地表使用實測地形圖生成，其他各層則根據鉆孔資料并結合surfer平臺，采用克里金插值法生成初步模型，再進行微調生成。模型共包括節點14 820個，單元18 790個。研究區單元的劃分不均一，在145～175 m高程范圍內進行了網格加密。

圖4 舌狀涉水岸坡三維滲流模型Fig.4 Three-dimensional seepage numerical model of groundwater in the tongue-shaped bank slope

3.2 邊界條件、參數選取和初始條件

3.2.1 邊界條件

模型的邊界條件是指確定模型中各單元的水頭性質，用以判定是定水頭單元、變水頭單元或是無效水頭單元。江東咀庫岸段斜坡整體呈舌狀，處于三面涉水的環境中，僅向東北側延伸，由于庫水位周期性升降，漲落高差達30 m，故認為模型東北側是變水頭較符合實際情況；模型東側邊界鉆孔基本無水，庫水位的變化對其影響不大，可作為不透水邊界；其他三面的變水頭邊界則視表1中的工況予以確定。

3.2.2 參數選取

本次調查期間，庫區處于降水期，但庫水位保持在166 m以上，鉆孔顯示地下水水位基本與其持平。隨著庫水位的下降，部分單元將變為干單元，在這個過程中，利用三維滲流模型計算地下水在這些單元的變化時所用的參數也由儲水系數S轉變為給水度μ，因此三維滲流模型所涉及的計算參數包括滲透系數K(水平方向Kx、Ky，垂直方向Kz)、給水度μ和儲水系數S。限于所采用的手段和可用的資料，本文將巖體近似作為一種連續、均質的介質。

經現場鉆孔抽水和注水試驗，粉質黏土夾碎石和碎塊石土的滲透系數分別為2.64×10-6m/s(即0.228 m/d)、3.63×10-6m/s(即0.314 m/d)，第四系覆蓋層的滲透系數綜合取值為0.3 m/d；基巖巖體裂隙發育，本身較破碎，但普遍存在泥質、鈣質填充、膠結的現象，其滲透性較弱，可將其視為隔水層，滲透系數取值為10-5m/d。給水度μ和儲水系數S按經驗取值，具體取值詳見表2。

表2 三維滲流模型的計算參數

3.2.3 初始條件

初始條件即初始水位的確定。目前共有三種方法可以用來確定庫區地下水水位的分布情況：一是根據鉆孔中地下水水位的埋深資料，利用插值方法或趨勢面分析得到；二是先利用類比方法確定該地區的水力坡度，利用該水力坡度和現今的庫水位推算出各單元的地下水水位埋深；三是對模型進行穩定流模擬，利用模擬的結果作為模型的初始水位。由于野外鉆探工作在庫區降水位期間進行，且鉆孔水位隨庫水位變化明顯，實時記錄存在誤差，無法直接用來確定初始水位，而類比方法精度較差，故本次模擬采用了第三種方法，即將第一個水文年的穩態模擬結果作為下一個水文年的初始水位。

4 模擬結果與分析

岸坡地下水三維滲流場數值模擬分兩個水文年進行，第一個水文年的模擬結果用以確定初始水位，重點分析在第二個水文年中舌狀涉水岸坡地下水三維滲流場特征及其變化規律。

第一個水文年模擬三峽庫區首次蓄水至175 m庫水位的工況：首先，庫水位在30 d內由145 m上升至175 m，升幅30 m，平均上升速度為1 m/d；水庫蓄水至175 m后，開始進入正常運行階段；然后，經過180 d(11月初至次年4月末)后，在汛期前大幅度緩慢下降，即在大約45 d(每年5月初至6月中旬)內庫水位從175 m降至防洪限制水位145 m，降幅30 m，平均下降速度約為0.67 m/d；最后，保持145 m低庫水位運行105 d，地下水水位分布見圖5(a)，并且該水位將作為第二個水文年的初始水位。

第二個水文年模擬水庫正常運行期間，庫水位對地下水滲流場的影響，其中庫水位變化與第一個水文年相同，其模擬結果見圖5至圖8。

圖5 工況Ⅱ下岸坡地下水水位分布圖Fig.5 Distribution map of groundwater level of the bank slope in working condition Ⅱ

由圖5可見，在庫水位由145 m上升至175 m期間，舌狀涉水岸坡內地下水水位變化有很大的區別：在大和尚包以西范圍內地下水水位值的變化最快，當庫水位上升至175 m，地下水水位幾乎同步上升至175 m，這主要是由于該處地形狹窄，滲流路徑比較短，另外該區域內地表高程(174～176 m)與庫水位十分接近；大和尚包以東至紅梁子溝一帶，地下水水位變化也相對較快，同時在大、小和尚包一帶出現類似“地下水漏斗”的滲流特征，即中部水位低而周邊水位高，175 m地下水水位等值線平面范圍逐漸向內擴展;在紅梁子溝—小狗架子溝東側區域，由于相對遠離庫水且靠近基巖面，地下水水位的變化較緩慢。

由圖6可見，當水庫保持175 m高庫水位運行時，舌狀涉水岸坡地下水水位逐漸上升。在運行的前90 d內，地下水水位的變化較明顯；而90 d后，地下水水位等值線的空間范圍和數值基本都不再發生變化，除墳院子和神女廟附近地下水水位依然保持145 m左右，其余各處地下水水位均達到175 m。分析原因認為：紅梁子溝—小狗架子溝以西，基巖埋深較深，第四系覆蓋物較厚，且地形較狹窄，滲流路徑較短，地表高程普遍低于190 m，地下水可以很快與庫水產生聯系，并達到平衡；而紅梁子溝—小狗架子溝以東則剛好相反，不僅海拔高，而且地形寬闊，導致滲流路徑變長，加之作為相對隔水層的基巖埋深淺且基巖面陡峭，因此當保持175 m高庫水位運行180 d時，高程175 m處基巖面附近的第四系覆蓋層部分可以達到175 m的地下水水位，而基巖內仍然保持穩定流模擬的145 m地下水水位，其變化梯度和范圍在理論上存在，但限于數值模擬的精度，很難在圖形中予以區別。

圖6 工況Ⅲ下岸坡地下水水位分布圖Fig.6 Distribution map of groundwater level of the bank slope in working condition Ⅲ

由圖7可見，庫水位由175 m下降至145 m，舌狀涉水岸坡內地下水水位變化與庫水位上升過程有所不同。庫水位下降過程中，地下水水位值的變化慢慢減緩，小和尚包至墳院子一帶在庫水位下降至145 m過程中，始終保持170 m以上的地下水水位，這主要是由于地形坡度導致滲流路徑的變長。

由圖8可見，當水庫保持145 m低庫水位運行時，舌狀涉水岸坡地下水水位逐漸下降，地下水水位等值線整體上表現為“西疏東密”，即地下水水位在西部區域變化最快，往東逐漸減緩；保持145 m低庫水位運行105 d后，岸坡內地下水滲流場仍未達到穩態，地下水仍然在向長江和大寧河排泄，但新一輪的庫水位上升即將開始，地下水滲流場將重新進行調整。

圖7 工況Ⅳ下岸坡地下水水位分布圖Fig.7 Distribution map of groundwater level of the bank slope in working condition Ⅳ

圖8 工況Ⅰ下岸坡地下水水位分布圖Fig.8 Distribution map of groundwater level of the bank slope in working conditionⅠ

通過對比分析圖5至圖8可知，舌狀涉水岸坡地下水三維滲流場還具有以下特征：

(1) 地下水水位分布與地形有極大的相似性，沿舌狀岸坡山脊線呈對稱分布，且地下水水位隨地形的起伏而產生相應的變化，但變化的幅度略小于地形坡度。地下水水位等值線整體上表現為“西疏東密”，主要受地形坡度的控制，與基巖(相對隔水層)的埋深也有一定的關系。

(2) 地下水水位上升和下降的速度明顯滯后于庫水位的變化。這主要是由于庫水位的日變化幅度遠大于岸坡的入滲量，即巖土體的滲透系數太小。另外，通過對比分析圖5(a)和圖8(f)可知，相同時間節點，第二個水文年的地下水水位略高于第一個水文年的地下水水位，證明地下水在坡體內有一定的“滯留積累”效應，部分區域地下水水位在水庫保持145 m低庫水位運行105 d后仍然保持在170 m以上。

(3) 庫水位的變化對紅梁子溝—小狗架子溝以西區域地下水滲流場的影響較大，地下水水位等值線的變化更為明顯，而東側區域地下水滲流場則較為穩定。

本次數值模擬研究能基本反映在庫水調動的條件下，該舌狀涉水岸坡體內地下水三維滲流場的變化特征。但也存在些許誤差，如東部區域基巖范圍內一直保持著穩定流模擬的145 m地下水水位，地下水水位的變化梯度和范圍難以看出。這種誤差產生的原因包括兩個方面：①建模誤差。建模過程中，簡化了岸坡的物質組成結構，且限于計算機硬件的限制，網格的劃分也不能太密，因此無法保證岸坡的三維數值模型與其真實的地質結構保持完全一致；②參數誤差。滲透系數采用了野外抽水試驗測得的參數，而給水度和貯水系數無法在野外測得，則采用了經驗系數，存在一定的誤差。

5 結 論

本文運用有限元方法，通過三維數值模擬軟件FEFLOW，對簡化條件下的涉水岸坡地下水滲流場進行了三維數值模擬研究，并得出了以下結論：

(1) 岸坡地下水水位變化與庫水位呈正相關性，表現為地下水水位隨庫水位的升高而升高，且這種變化具有明顯的“滯后效應”，存在地下水滲流中的“漏斗”現象。

(2) 岸坡地下水水位分布與地形密切相關，紅梁子溝—小狗架子溝以西區域地下水滲流場受庫水的影響較大，而東側區域地下水滲流場則較為穩定，且坡體內的地下水水位沿舌狀涉水岸坡山脊線呈對稱分布。在地形狹窄的區域，由于滲流路徑短且緊臨庫水，地下水水位變化較快，表現為“西疏東密”的分布特征；庫水位下降階段，地形坡度的存在導致滲流路徑變長，因此岸坡地下水水位的變化也比庫水位上升階段要緩慢。

(3) 岸坡地下水水位的分布與基巖面的位置(相對隔水層)有關。該舌狀涉水岸坡東部區域基巖埋深淺且基巖面陡峭，其附近的覆蓋層區域地下水水位等值線分布密集，而基巖內部由于相對隔水，其庫水位值一直保持初始穩態分析的145 m低庫水位。