，，

(1.江蘇省地質(zhì)調(diào)查研究院，江蘇 南京 210018；2.南京市測繪勘察研究院股份有限公司，江蘇 南京 210005；3.河海大學(xué)水文水資源學(xué)院/河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098)

自由滲流面與地下水潛水水面特征相似，具有復(fù)雜的非線性，是河流等地表水入滲補(bǔ)給地下水、大壩滲流等定量化計算的難點，目前，尚未發(fā)現(xiàn)確定自由滲流面的解析法，在實際工程應(yīng)用中常用數(shù)值方法求解。求解該類問題的數(shù)值方法可分為兩類：變網(wǎng)格法和固定網(wǎng)格法。由于變網(wǎng)格法操作復(fù)雜，容易使計算網(wǎng)格畸變，難以處理存在水平介質(zhì)分層及各種復(fù)雜夾層的情況等缺點,因而研究較少。為此，國內(nèi)外學(xué)者提出了許多與固定網(wǎng)格法相關(guān)的處理技術(shù)，主要有剩余流量法、初流量法、節(jié)點虛流量法、單元滲透矩陣調(diào)整法、復(fù)合單元法及截止負(fù)壓法等[1-6]。但這些處理技術(shù)大多方法復(fù)雜，操作困難，應(yīng)用不便。

本文利用開源地下水?dāng)?shù)值模擬軟件MODFLOW[7]及SUTRA[8]，提出了采用干濕單元轉(zhuǎn)化、緩變滲透系數(shù)矩陣法推求自由滲流面的方法，成功推求了自由滲流面，有效地提高了目前河流與地下水交互作用計算的求解精度。

1 典型自由滲流面問題

典型自由滲流面問題可概化為圖1所示。為簡化求解難度，一般假設(shè)：(1)土壤均質(zhì)各向同性；(2)毛管作用及蒸散發(fā)作用忽略；(3)水流服從Darcy定律；(4)穩(wěn)定水流。基于以上假設(shè)可建立滲流問題的數(shù)學(xué)模型：

(1)

式中：h為水頭/m；y為相對高程/m；Ks為飽和滲透系數(shù)/(m·s-1)。

為便于建立數(shù)值模型，并不失一般性，假定上、下河道水頭y1=6 m，y2=1 m，滲流區(qū)距離x1=4 m。對于均質(zhì)各向同性土壤介質(zhì)中穩(wěn)定流問題，滲透系數(shù)對自由面計算無影響，故可設(shè)定為任意值，如采用Ks=10-5m/s。

圖1 滲流概化圖

2 數(shù)值法求解

2.1 MODFLOW求解方法

MODFLOW是美國地質(zhì)調(diào)查局開發(fā)的模擬三維地下水運動的數(shù)值計算軟件。該軟件采用有限單元差分法離散地下水水流方程，廣泛應(yīng)用于地下水?dāng)?shù)值模擬及地下水資源管理。MODFLOW將單元格分為四類：無效單元格、完全飽和單元格、部分飽和單元格、定水頭單元格。前三種類型可通過輸干及再濕潤等處理方法相互轉(zhuǎn)化。當(dāng)計算單元(i,j,k-1)水頭hi,j,k-1大于底板高程hb(hi,j,k-1>hb)時，無效單元格轉(zhuǎn)為有效計算單元格，初始水頭hi，j，k=C(hi,j,k-1-hb)+hb，其中C為水力傳導(dǎo)系數(shù)，該過程稱為浸潤；當(dāng)hi，j，k

自由滲流面附近的計算單元格為部分飽和單元， MODFLOW通過將這些單元的導(dǎo)水系數(shù)設(shè)置為T=Δh×Ks(其中Δh=hi，j，k-hb)，再利用調(diào)和平均、對數(shù)平均等處理技術(shù)計算單元格間的傳導(dǎo)能力，很好地解決了這些單元格的水流計算問題[7]。

在MODFLOW推求自由面問題方面，可采用排水溝方法模擬穩(wěn)定自由滲流面，但排水溝水位設(shè)置十分重要，如果排水溝水位設(shè)置不當(dāng)，引起干濕單元轉(zhuǎn)化混亂，溢出點附近水位有奇異變化，從而致使計算不穩(wěn)定[9]。

本文將圖1中滲流區(qū)域劃分為300×200個計算單元，設(shè)定MODFLOW為穩(wěn)定流計算模式，式(1)中邊界條件設(shè)置如下：af、bc概化為定水頭邊界；ab不透水邊界(無需直接設(shè)置，MODFLOW默認(rèn)的邊界即為不透水邊界)；fd為自由滲流面，無法直接給定，可在計算迭代中求出。具體計算方法為：利用MODFLOW中計算單元格的有效、無效及定水頭屬性設(shè)置功能，將自由面右側(cè)的單元格設(shè)置為無效單元格，這是因為自由滲流面為單調(diào)函數(shù)，其右側(cè)無水流運動，通過迭代不斷調(diào)整計算單元格活動、非活動屬性直至收斂；cd滲流邊界，將該處單元格設(shè)置為定水頭，其水頭設(shè)置為單元格頂板高程，由于MODFLOW采用單元格水深計算單元格導(dǎo)水系數(shù)，這種設(shè)置可保證溢流面排水通暢;在初始設(shè)置中可將整個bg設(shè)置為滲流邊界，迭代過程中非滲流單元格會被剔除(被設(shè)置成了無效單元格)。按以上方法設(shè)置的數(shù)值模型求解得到的滲流區(qū)自由水面線及流場水頭分布分別見圖2和圖3。

2.2 SUTRA求解方法

SUTRA是美國地質(zhì)調(diào)查局研發(fā)的用于模擬變密度、變飽和度水流運動及溶質(zhì)運移的數(shù)值模擬軟件。該軟件采用Galerkin有限元法離散非飽和水流運動方程(Richard方程)及溶質(zhì)運移方程，廣泛應(yīng)用于海水入侵模擬研究[8][10]。

利用SUTRA構(gòu)建圖1中滲流區(qū)域的數(shù)值模型，網(wǎng)格剖分及邊界條件處理與MODFLOW相近：af、bc概化為定水頭邊界，ab為不透水邊界；fd為自由滲流面，無需直接給定。但對于滲流邊界cd，由于SUTRA不能剔除非計算單元格，因此不能采用定水頭模擬，而邊界cg段計算節(jié)點具有如下特征：(1)溢流段cd，計算節(jié)點壓強(qiáng)為0；(2)非溢流段dg，溢流量為0，因此可在SUTRA中將cg設(shè)置為變水頭邊界，依據(jù)臨近計算節(jié)點判定cg段計算節(jié)點壓強(qiáng)值：若臨近節(jié)點壓強(qiáng)大于0(飽和)，則該節(jié)點位于滲流面處，壓強(qiáng)為0，反之，可判定該節(jié)點位于非滲流面外，壓強(qiáng)設(shè)置為臨近節(jié)點壓強(qiáng)。

SUTRA程序未對計算節(jié)點分類，不能隨意剔除非水流區(qū)計算節(jié)點。通過設(shè)置滲透系數(shù)來區(qū)分水流區(qū)及非水流區(qū)計算節(jié)點的單元滲透矩陣調(diào)整法處理技術(shù)，具有方便簡潔優(yōu)點。黨發(fā)寧等人[11]提出了變單元滲透系數(shù)法，該方法將自由面以下的計算單元設(shè)置為飽和滲透系數(shù)，自由面以上單元近似設(shè)置為0，重新建立總體滲透矩陣，試圖消除非飽和區(qū)對自由面滲流求解的影響，但是該處理方法會造成自由面附近的高斯點求解結(jié)果出現(xiàn)振蕩、穩(wěn)定性差等問題。馬淑芝等[5]、付延玲等[12]提出并改進(jìn)了復(fù)合單元滲透矩陣調(diào)整法，他們將罰函數(shù)引入自由面附近單元格，減少離散矩陣的奇異性。裴利華[13]研究認(rèn)為單元滲透矩陣調(diào)整法求出的單元滲透系數(shù)矩陣不能真實反映滲流區(qū)透水特性，且矩陣的主系數(shù)不占優(yōu)，影響計算精度和穩(wěn)定性，并且對于三維問題，自由面穿越的單元形式很復(fù)雜，數(shù)學(xué)上不易處理，程序工作量很大。

本文先按傳統(tǒng)的單元滲透矩陣調(diào)整法將計算結(jié)點分為兩類：飽和節(jié)點與非飽和節(jié)點，非飽和節(jié)點滲透系數(shù)近似為0，飽和節(jié)點滲透系數(shù)為飽和滲透系數(shù)。這種設(shè)置求出的滲流區(qū)自由水面線及流場水頭分布分別見圖2和圖3。在溢出點附近計算水頭出現(xiàn)奇異變化，計算過程中出現(xiàn)不穩(wěn)定，且較難收斂，這主要是由于滲透系數(shù)矩陣突變型設(shè)置導(dǎo)致雅克比矩陣不能滿足嚴(yán)格對角占優(yōu)，采用迭代法求解線性方程會引起數(shù)值震動。因此為提高求解結(jié)果的穩(wěn)定性，需引入緩變型滲透系數(shù)矩陣。由于非飽和滲透系數(shù)函數(shù)描述從飽和到干涸過程滲透系數(shù)的緩變過程，因此，通過引入非飽和滲透系數(shù)函數(shù)可解決傳統(tǒng)單元滲透矩陣調(diào)整法不穩(wěn)定性問題。

SUTRA采用Van Genuchten-Mualem(VGM)模型描述非飽和滲透系數(shù)[14]：

(2)

式中：p為水壓/pa；K為土壤非飽和滲透系數(shù)/(m·s-1)；α，n，m為VGM模型經(jīng)驗參數(shù)，其中m=1-1/n。

對于非飽和區(qū)滲透系數(shù)K，由式(2)可知，當(dāng)α→∞時，K→0。因此，如采用α取較大值方法，進(jìn)行飽和、非飽和水流計算，求解自由面問題，可避免在益處點附近計算水頭奇異變化及數(shù)值解不穩(wěn)定問題。當(dāng)α=0.002/pa、α=0.000 05/pa、n=2時計算的滲流區(qū)自由水面線及流場水頭分布分別見圖2和圖3。由圖可以看出，當(dāng)值α較大時，自由水面線推求結(jié)果雖然仍存在水頭震蕩變化，但計算穩(wěn)定性有較大提高。

圖2 三種方法求得的自由水面線(a)及其局部放大圖(b)

表1 邊界出入流量表 ×10-5/m2·s-1

(a)MODFLOW(b)SUTRA變滲透率法(c)SUTRA改進(jìn)變滲透率法

MODFLOW、SUTRA傳統(tǒng)單元滲透矩陣調(diào)整法和SUTRA緩變滲透系數(shù)矩陣法求得的溢出點(d點)高程分別為3.25 m、3.22 m、3.24 m，與甘油試驗[15]結(jié)果(3.25 m)對比可知，MODFLOW求解精度最高。圖2和圖3中比較三種方法表明，緩變滲透系數(shù)矩陣法相對于傳統(tǒng)單元滲透矩陣調(diào)整法穩(wěn)定性有較大提高。三種方法求得的邊界出af段入流量、bc段出流量及cd段溢流量見表1，MODFLOW和SUTRA緩變滲透系數(shù)矩陣法求得的af段入流量一致，MODFLOW、SUTRA傳統(tǒng)單元滲透矩陣調(diào)整法和SUTRA緩變滲透系數(shù)矩陣法求得的bc段出流量無顯著差異，但SUTRA模擬的cd溢出段出流量比MODFLOW模擬結(jié)果偏小。

3 結(jié)語

MODFLOW、SUTRA傳統(tǒng)單元滲透矩陣調(diào)整法及SUTRA緩變滲透系數(shù)矩陣法求解結(jié)果表明采用MODFLOW運用干濕轉(zhuǎn)化技術(shù)求解自由滲透面的方法穩(wěn)定性最好、精度最高，但MODFLOW中的有限單元差分方法要求計算網(wǎng)格為矩形網(wǎng)格，矩形網(wǎng)格難于刻畫復(fù)雜形體，致使MODFLOW應(yīng)用方面存在較大局限性，而采用緩變滲透系數(shù)矩陣法的SUTRA程序，改善了傳統(tǒng)單元滲透矩陣調(diào)整法的不穩(wěn)定性，提高了數(shù)值計算精度，避免了MODFLOW必須矩形網(wǎng)格的局限性，是一種計算自由滲流面，估算地下水與河流水量交換量的實用方法。