基于FEFLOW的岸坡地下水三維滲流場(chǎng)模擬研究

李慶偉，晏鄂川，楊 廣，崔學(xué)杰

(1.中國(guó)地質(zhì)大學(xué)(武漢)工程學(xué)院，湖北 武漢 430074；2.中交第二航務(wù)工程勘查設(shè)計(jì)院有限公司，湖北 武漢 430060)

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

基于以上認(rèn)識(shí)，本文在三峽水庫(kù)蓄水條件下，通過構(gòu)建涉水岸坡的滲流模型，基于FEFLOW軟件對(duì)岸坡地下水滲流場(chǎng)的變化情況進(jìn)行了三維數(shù)值模擬研究，從而揭示水庫(kù)正常運(yùn)行期間庫(kù)岸內(nèi)部地下水滲流場(chǎng)的發(fā)展動(dòng)向與變化趨勢(shì)。

1 涉水岸坡工程地質(zhì)概況

圖1 研究區(qū)地層分布圖Fig.1 Formation distribution map

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

2 計(jì)算原理與方法

2.1 計(jì)算原理

當(dāng)不考慮水的密度的變化條件下，在孔隙介質(zhì)中地下水在三維空間的流動(dòng)可以用下面的偏微分方程來(lái)表示[13]：

(1)

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

公式(1)構(gòu)成了對(duì)于一個(gè)實(shí)際問題地下水流動(dòng)的定解問題，可采用數(shù)值計(jì)算方法進(jìn)行求解，如有限差分法、有限單元法等，求解結(jié)果即為地下水水位的分布值。

2.2 計(jì)算方法

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

3 數(shù)值模擬方案設(shè)計(jì)

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

圖2 三峽工程正常運(yùn)行期水庫(kù)調(diào)度圖Fig.2 Reservoir operation chart of Three Gorges Project during normal operation period注：(1)為防破壞線；(2)為限制供水線；(Ⅰ)為防洪區(qū)；(Ⅱ)為裝機(jī)預(yù)想出力區(qū)；(Ⅲ)為保證出力區(qū)；(Ⅳ)為降低出力區(qū)(資料來(lái)源于《三峽庫(kù)區(qū)地質(zhì)災(zāi)害防治工程地質(zhì)勘查技術(shù)要求》)

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

3.1 數(shù)值模型的構(gòu)建

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

圖3 數(shù)值模型范圍Fig.3 Domain of numerical model注：zk22〇為鉆孔編號(hào)及位置

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

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

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

3.2 邊界條件、參數(shù)選取和初始條件

3.2.1 邊界條件

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

3.2.2 參數(shù)選取

本次調(diào)查期間，庫(kù)區(qū)處于降水期，但庫(kù)水位保持在166 m以上，鉆孔顯示地下水水位基本與其持平。隨著庫(kù)水位的下降，部分單元將變?yōu)楦蓡卧?，在這個(gè)過程中，利用三維滲流模型計(jì)算地下水在這些單元的變化時(shí)所用的參數(shù)也由儲(chǔ)水系數(shù)S轉(zhuǎn)變?yōu)榻o水度μ，因此三維滲流模型所涉及的計(jì)算參數(shù)包括滲透系數(shù)K(水平方向Kx、Ky，垂直方向Kz)、給水度μ和儲(chǔ)水系數(shù)S。限于所采用的手段和可用的資料，本文將巖體近似作為一種連續(xù)、均質(zhì)的介質(zhì)。

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

表2 三維滲流模型的計(jì)算參數(shù)

3.2.3 初始條件

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

4 模擬結(jié)果與分析

岸坡地下水三維滲流場(chǎng)數(shù)值模擬分兩個(gè)水文年進(jìn)行，第一個(gè)水文年的模擬結(jié)果用以確定初始水位，重點(diǎn)分析在第二個(gè)水文年中舌狀涉水岸坡地下水三維滲流場(chǎng)特征及其變化規(guī)律。

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

第二個(gè)水文年模擬水庫(kù)正常運(yùn)行期間，庫(kù)水位對(duì)地下水滲流場(chǎng)的影響，其中庫(kù)水位變化與第一個(gè)水文年相同，其模擬結(jié)果見圖5至圖8。

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

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

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

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

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

由圖8可見，當(dāng)水庫(kù)保持145 m低庫(kù)水位運(yùn)行時(shí)，舌狀涉水岸坡地下水水位逐漸下降，地下水水位等值線整體上表現(xiàn)為“西疏東密”，即地下水水位在西部區(qū)域變化最快，往東逐漸減緩；保持145 m低庫(kù)水位運(yùn)行105 d后，岸坡內(nèi)地下水滲流場(chǎng)仍未達(dá)到穩(wěn)態(tài)，地下水仍然在向長(zhǎng)江和大寧河排泄，但新一輪的庫(kù)水位上升即將開始，地下水滲流場(chǎng)將重新進(jìn)行調(diào)整。

圖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Ⅰ

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

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

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

(3) 庫(kù)水位的變化對(duì)紅梁子溝—小狗架子溝以西區(qū)域地下水滲流場(chǎng)的影響較大，地下水水位等值線的變化更為明顯，而東側(cè)區(qū)域地下水滲流場(chǎng)則較為穩(wěn)定。

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

5 結(jié) 論

本文運(yùn)用有限元方法，通過三維數(shù)值模擬軟件FEFLOW，對(duì)簡(jiǎn)化條件下的涉水岸坡地下水滲流場(chǎng)進(jìn)行了三維數(shù)值模擬研究，并得出了以下結(jié)論：

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

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

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