MODFLOW在庫區(qū)浸沒預(yù)測上的應(yīng)用

□陳 鵬 □孫 剛 □周子?xùn)|（河南省水利勘測有限公司）

0 引言

預(yù)測水庫蓄水后庫周水位分布情況是研究浸沒發(fā)生范圍、程度的關(guān)鍵一步，在查清地塊水文工程地質(zhì)條件的前提下運(yùn)用MODFLOW可以幫助我們較準(zhǔn)確的預(yù)測蓄水后水位的變化，從而對比臨界埋深得出可靠的結(jié)論。

1 MODFLOW簡介

MODFLOW是地下水模擬系統(tǒng) (GroundwaterModeling Systems)的一個計算模塊，是美國地質(zhì)調(diào)查局在20世紀(jì)80年代開發(fā)出來的專門用于孔隙介質(zhì)中三維有限差分地下水流數(shù)值模擬的軟件，問世以來，已經(jīng)在全世界范圍內(nèi)科研、生產(chǎn)、環(huán)境保護(hù)和水資源利用等領(lǐng)域得到了廣泛應(yīng)用，成為最為普及的地下水運(yùn)動數(shù)值模擬的計算軟件。應(yīng)用MODFLOW進(jìn)行地下水流模擬的步驟如下：

滲流區(qū)的離散化→建立地下水流動問題的差分方程組→求解差分方程組。

2 研究地塊的工程水文地質(zhì)條件

擬建水庫位于淮河干流上游。因水庫左岸地形平緩開闊，第四紀(jì)沖洪積物分布廣泛，地下水位埋深較淺，初判將遭受較大范圍浸沒災(zāi)害的影響。選擇了龍井～下土城地塊為作為浸沒勘察研究的典型地塊。本區(qū)氣候?qū)偃A中副熱帶氣候區(qū)。多年平均降雨量為1059mm。多年平均水面蒸發(fā)量790mm，陸面蒸發(fā)量650mm。

該地塊地形總體由東北向西南河谷緩傾，地塊以北為由中更新世沖洪積物組成的崗地，向南依次降為二級階地、一級階地及河床漫灘，地表坡降多在2.00×10-3～2.00×10-2之間。

據(jù)現(xiàn)場調(diào)查及鉆探資料可知，本區(qū)分布的地層中有第三系和第四系，由老到新簡述如下：一是第三系毛家坡組（Em）以紫紅色粘土巖為主，夾有泥質(zhì)填充的砂礫巖；二是中更新統(tǒng)（alplQ2）為褐紅色低液限粘土；三是上更新統(tǒng)（alplQ3）上部為低液限粘土，下部為砂礫石；四是全新統(tǒng)（alQ4）上部為低液限粘土，下部為細(xì)砂土，構(gòu)成一級階地的主體。水庫水邊線位于一級階地。

表1中的水文地質(zhì)參數(shù)是通過現(xiàn)場和室內(nèi)試驗，結(jié)合類似工程和經(jīng)驗確定。底部第三系粘土巖和砂礫巖，透水性小，可作為本區(qū)相對隔水層。場區(qū)等水位線近似平行于河床，其補(bǔ)給來源以后緣中更新統(tǒng)組成的崗地以潛流形式補(bǔ)給，其次為降水補(bǔ)給。地下水的排泄方式有兩種：一是以潛流和泉的形式向游河排泄；二是通過民井排泄。

表1 龍井～下土城地塊水文地質(zhì)參數(shù)建議值表

龍井～下土城地塊可能浸沒范圍內(nèi)主要為一級階地Q4低液限粘土，通過野外試驗確定土壤毛管水上升帶的高度取62.50 cm。

3 地塊的地下水流動數(shù)學(xué)模型

研究區(qū)的地下水主要儲存于砂礫層中，存在統(tǒng)一的地下水面，地下水流動系統(tǒng)遵從水均衡原理和達(dá)西定律。研究區(qū)域地下水流從空間上看以水平運(yùn)動為主，含水介質(zhì)的巖性特征在垂向上變化較顯著，為保證計算準(zhǔn)確性將水庫蓄水前后的地下水運(yùn)動均概化為非均質(zhì)各向同性三維穩(wěn)定流。

根據(jù)上述的概念模型，結(jié)合研究區(qū)的邊界條件，可得出研究區(qū)地下水流動數(shù)學(xué)模型：

式中：h0—初始水頭；H—含水層厚度；n—含水層邊界上內(nèi)法向單位矢量；Ω—模擬區(qū)范圍；B1—模擬區(qū)上下游定流量邊界；B2—模擬區(qū)兩側(cè)零流量邊界；q—含水層邊界單寬流量，流入為正，流出為負(fù)。

4 地塊的空間離散化與邊界條件的確定

龍井～下土城地塊水文地質(zhì)結(jié)構(gòu)主要為非均質(zhì)多層水文地質(zhì)結(jié)構(gòu)。模擬區(qū)域長1980m，寬1200m的地塊，模擬區(qū)面積為1.90 km2，見圖1。模型兩側(cè)邊界分別取在2個勘探線外側(cè)，并與地下水等水位線垂直，為零流量邊界；下游邊界取已知水頭邊界，其值為各處的河水位；上游邊界取在Q2和Q3地層界線附近，為定流量邊界，其值可根據(jù)達(dá)西定律確定：

依據(jù)實測的地下水等水位線，分別計算出弱透水層和強(qiáng)透水層流量，最后疊加得到模擬區(qū)上游邊界流量約2000m3/d。

利用GMS(Groundwater Modeling System)中Modflow模塊建立柵格模型。按模擬區(qū)地層巖性將模擬區(qū)分成4層，每層33×20個單元格，其中每層有效單元格為577個，總共有效單元格為2308個,見圖2。分別給出4層的頂?shù)装甯邔?，然后賦給各單元格水文地質(zhì)參數(shù)，并進(jìn)行地下水流動模擬。實測水位時間為3月份，因降雨較少，所以模型識別時不考慮降雨。

圖1 龍井～下土城地塊模擬區(qū)邊界條件示意圖

5 典型地塊數(shù)值模型的識別與校驗

由對水文地質(zhì)結(jié)構(gòu)的認(rèn)識到水文地質(zhì)概念模型再到地下水?dāng)?shù)值模型的整個過程經(jīng)歷了一系列概化處理，這些因素導(dǎo)致模型不可能細(xì)致地、完全真實地刻畫出各處的地下水運(yùn)動，而只能對模擬區(qū)內(nèi)地下水運(yùn)動進(jìn)行宏觀的、趨勢性的擬合。為了判別所概化的模型與實際地質(zhì)原型的符合程度，必須對模型進(jìn)行識別與檢驗。

采用試算法進(jìn)行校驗，不斷調(diào)整含水層參數(shù)，使模擬得到的水位與實測水位吻合。模擬的水位與實測水位有一定的差別，這是因為：①概化的水文地質(zhì)概念模型是理想化的表述，例如巖層中存在的透鏡體難以一一刻畫；②邊界條件都是人為邊界，與實際的邊界條件也有一些差別；實測的等水位線是根據(jù)離散的水位值經(jīng)插值而得，其精度有一定局限。

盡管如此，模擬所得的水位和觀測水位形狀相似，趨勢相近，兩者誤差最大不超過0.50m，地下水流動方向和觀測的流動方向基本一致。模型可以用于預(yù)測典型地塊水庫蓄水至90m時地塊的等水位線變化情況。計算參數(shù)采用經(jīng)模型識別調(diào)整后的計算值，如表2所示。

表2 龍井～下土城地塊浸沒預(yù)測計算采用水文地質(zhì)參數(shù)表

6 水庫蓄水后地塊邊界條件的確定

水庫蓄水后，在龍井～上土城形成三面臨水的半島，其地下水流動方向?qū)l(fā)生變化。模擬區(qū)下游邊界為庫水位定水頭邊界，上游邊界偏安全考慮，仍然為定流量邊界，但其水位將隨庫水位上升有所壅高。取平行于流線的一個剖面，下游邊界水位初始高程h為85m，蓄水后至90m；上游邊界初始水位h1為90.50m，距離庫水邊界約900m，可以近似計算得到庫水位上升后上游邊界的水位。當(dāng)庫水位上升至90m時，上游邊界水位上升至90.89m，偏安全考慮取91m。蓄水后龍井～下土城地塊模擬區(qū)面積為1.60 km2。

降雨入滲作為模擬區(qū)上部邊界條件，在模擬初始水位時，由于實測水位的時間為3月下旬至4月上旬，降雨量小，所以在模型識別與校驗時未考慮降雨影響。在預(yù)測水庫蓄水后，典型地塊地下水位壅高時，應(yīng)當(dāng)考慮降雨入滲的影響。

模型計算時，偏安全考慮取6-8月降雨集中時間，平均日降雨量為2.50mm/d，兩典型地塊地表都以低液限粘土為主，其入滲系數(shù)取經(jīng)驗值0.01，得到降雨入滲量為2.50×10-5m/d。

7 水庫蓄水后典型地塊地下水位預(yù)測

如前所調(diào)整的模型邊界及計算參數(shù)，代入模型計算，可以得到地塊模擬區(qū)在水庫蓄水以后的等水位線圖，進(jìn)而作出模擬區(qū)的地下水埋深等值線圖。在龍井～下土城地塊模擬區(qū)，在蓄水后模擬區(qū)的滲流場特征比蓄水前發(fā)生了明顯變化。其一是隨著庫水位上移至高程90m，從宏觀上看，全區(qū)地下水位上升后的運(yùn)動總趨勢大致仍由東北向西南流動。但在上土城和下土城將變成三面臨水的半島狀庫岸，地下水呈扇狀流動，分別向淮河和半島兩側(cè)的沖溝或低地（都是庫區(qū)）排泄。其二是隨著庫水位上升，地下水位發(fā)生不同程度的壅高，升幅多在0.40～4.50m。其中以沿庫周邊線附近上升較明顯，向后緣升幅逐漸減小，至后緣上游邊界附近升幅僅40～50 cm。其三是地下水埋深發(fā)生不同程度的變淺，最明顯的是沿庫邊線附近地下水埋深<1m的范圍有所擴(kuò)大。

8 臨界埋深的確定

依據(jù)《水利水電工程地質(zhì)勘察規(guī)范》（GB50487-2008），評價浸沒的地下水臨界埋深Hcr=土壤毛管水上升帶的高度（m）+安全超高值（m）。

庫區(qū)浸沒影響主要有農(nóng)作物、居民，農(nóng)作物主要為小麥、水稻、油菜等，最大根系發(fā)育深度約0.60m，一般民房基礎(chǔ)砌置深度為0.50m左右。故安全超高值取0.60m。龍井～下土城地塊Q4低液限粘土毛細(xì)上升高度取0.63m，可以得到浸沒的地下水臨界埋深為1.23m，偏安全取1.25m。

9 水庫蓄水后典型地塊浸沒預(yù)測

通過預(yù)測的地下水埋深Hm與浸沒的地下水臨界埋深Hcr進(jìn)行比較，可判斷模擬區(qū)的浸沒程度。當(dāng)?shù)叵滤裆頗m>Hcr時，不會產(chǎn)生浸沒。當(dāng)Hm≦Hcr時，將產(chǎn)生浸沒危害，即建筑物地基承載力下降，威脅建筑物安全穩(wěn)定；因過度潮濕使當(dāng)?shù)厝罕娚瞽h(huán)境惡化；因鹽漬化危及農(nóng)作物的正常生長等。其中當(dāng)Hm<0m時，將產(chǎn)生沼澤化，危害更大。

龍井～下土城地塊蓄水后模擬區(qū)地下水埋深在0～3m，地下水臨界埋深Hcr為1.25m，以此可以得到該地塊的浸沒范圍。當(dāng)庫水位達(dá)到90m時，該地塊模擬區(qū)浸沒寬度（90米庫水位至浸沒邊界）在60～780m之間不等，90m庫邊線長3300m，浸沒面積約為0.60 km2。受災(zāi)對象以農(nóng)田為主，其次沿庫岸邊線的村民和農(nóng)田也受到一定程度的危害。

10 典型地塊浸沒長度影響因素的敏感性分析

影響庫區(qū)浸沒的因素較多，如地層巖性、透水性、厚度空間分布、第四系潛水埋深、補(bǔ)給徑流排泄、第四系含水層與其下覆含水層的相互關(guān)系、氣候特征、庫水位高低、庫岸地形坡度等。庫區(qū)浸沒范圍大小受以上各種要素的綜合作用的影響。根據(jù)敏感度系數(shù)（SAF）的研究成果，認(rèn)為庫岸地形坡度是影響浸沒范圍的主要因素。

11 結(jié)語

在準(zhǔn)確掌握模擬地塊水文地質(zhì)條件前提下，MODFLOW可以幫助我們通過三維數(shù)值模擬預(yù)測水庫蓄水后庫周的地下水位，從而分析研究浸沒發(fā)生的范圍和程度。某水庫龍井～下土城地塊通過以上方法對浸沒問題進(jìn)行了評價，取得了可靠的結(jié)果。

[1]賀國平，張彤，趙月芬等.GMS數(shù)值建模方法研究綜述[J].地下水,2007(03).

[2]武強(qiáng)，朱斌，徐華等.MODFLOW在淮北地下水?dāng)?shù)值模擬中的應(yīng)用[J].遼寧工程技術(shù)大學(xué)學(xué)報,2005(04).

[3]陳全禮，張志敏.庫區(qū)蓄水浸沒預(yù)測評價[J].資源環(huán)境與工程,2008(B10).

[4]張穎，許模，王子忠.應(yīng)用三維地下水滲流數(shù)值模擬方法分析水庫浸沒問題[J].四川水利,2008(04).

[5]GB50487-2008,水利水電工程地質(zhì)勘察規(guī)范[S].