□陳 鵬 □孫 剛 □周子東(河南省水利勘測有限公司)
預測水庫蓄水后庫周水位分布情況是研究浸沒發生范圍、程度的關鍵一步,在查清地塊水文工程地質條件的前提下運用MODFLOW可以幫助我們較準確的預測蓄水后水位的變化,從而對比臨界埋深得出可靠的結論。
MODFLOW是地下水模擬系統 (GroundwaterModeling Systems)的一個計算模塊,是美國地質調查局在20世紀80年代開發出來的專門用于孔隙介質中三維有限差分地下水流數值模擬的軟件,問世以來,已經在全世界范圍內科研、生產、環境保護和水資源利用等領域得到了廣泛應用,成為最為普及的地下水運動數值模擬的計算軟件。應用MODFLOW進行地下水流模擬的步驟如下:
滲流區的離散化→建立地下水流動問題的差分方程組→求解差分方程組。
擬建水庫位于淮河干流上游。因水庫左岸地形平緩開闊,第四紀沖洪積物分布廣泛,地下水位埋深較淺,初判將遭受較大范圍浸沒災害的影響。選擇了龍井~下土城地塊為作為浸沒勘察研究的典型地塊。本區氣候屬華中副熱帶氣候區。多年平均降雨量為1059mm。多年平均水面蒸發量790mm,陸面蒸發量650mm。
該地塊地形總體由東北向西南河谷緩傾,地塊以北為由中更新世沖洪積物組成的崗地,向南依次降為二級階地、一級階地及河床漫灘,地表坡降多在2.00×10-3~2.00×10-2之間。
據現場調查及鉆探資料可知,本區分布的地層中有第三系和第四系,由老到新簡述如下:一是第三系毛家坡組(Em)以紫紅色粘土巖為主,夾有泥質填充的砂礫巖;二是中更新統(alplQ2)為褐紅色低液限粘土;三是上更新統(alplQ3)上部為低液限粘土,下部為砂礫石;四是全新統(alQ4)上部為低液限粘土,下部為細砂土,構成一級階地的主體。水庫水邊線位于一級階地。
表1中的水文地質參數是通過現場和室內試驗,結合類似工程和經驗確定。底部第三系粘土巖和砂礫巖,透水性小,可作為本區相對隔水層。場區等水位線近似平行于河床,其補給來源以后緣中更新統組成的崗地以潛流形式補給,其次為降水補給。地下水的排泄方式有兩種:一是以潛流和泉的形式向游河排泄;二是通過民井排泄。

表1 龍井~下土城地塊水文地質參數建議值表
龍井~下土城地塊可能浸沒范圍內主要為一級階地Q4低液限粘土,通過野外試驗確定土壤毛管水上升帶的高度取62.50 cm。
研究區的地下水主要儲存于砂礫層中,存在統一的地下水面,地下水流動系統遵從水均衡原理和達西定律。研究區域地下水流從空間上看以水平運動為主,含水介質的巖性特征在垂向上變化較顯著,為保證計算準確性將水庫蓄水前后的地下水運動均概化為非均質各向同性三維穩定流。
根據上述的概念模型,結合研究區的邊界條件,可得出研究區地下水流動數學模型:

式中:h0—初始水頭;H—含水層厚度;n—含水層邊界上內法向單位矢量;Ω—模擬區范圍;B1—模擬區上下游定流量邊界;B2—模擬區兩側零流量邊界;q—含水層邊界單寬流量,流入為正,流出為負。
龍井~下土城地塊水文地質結構主要為非均質多層水文地質結構。模擬區域長1980m,寬1200m的地塊,模擬區面積為1.90 km2,見圖1。模型兩側邊界分別取在2個勘探線外側,并與地下水等水位線垂直,為零流量邊界;下游邊界取已知水頭邊界,其值為各處的河水位;上游邊界取在Q2和Q3地層界線附近,為定流量邊界,其值可根據達西定律確定:

依據實測的地下水等水位線,分別計算出弱透水層和強透水層流量,最后疊加得到模擬區上游邊界流量約2000m3/d。
利用GMS(Groundwater Modeling System)中Modflow模塊建立柵格模型。按模擬區地層巖性將模擬區分成4層,每層33×20個單元格,其中每層有效單元格為577個,總共有效單元格為2308個,見圖2。分別給出4層的頂底板高層,然后賦給各單元格水文地質參數,并進行地下水流動模擬。實測水位時間為3月份,因降雨較少,所以模型識別時不考慮降雨。

圖1 龍井~下土城地塊模擬區邊界條件示意圖
由對水文地質結構的認識到水文地質概念模型再到地下水數值模型的整個過程經歷了一系列概化處理,這些因素導致模型不可能細致地、完全真實地刻畫出各處的地下水運動,而只能對模擬區內地下水運動進行宏觀的、趨勢性的擬合。為了判別所概化的模型與實際地質原型的符合程度,必須對模型進行識別與檢驗。

采用試算法進行校驗,不斷調整含水層參數,使模擬得到的水位與實測水位吻合。模擬的水位與實測水位有一定的差別,這是因為:①概化的水文地質概念模型是理想化的表述,例如巖層中存在的透鏡體難以一一刻畫;②邊界條件都是人為邊界,與實際的邊界條件也有一些差別;實測的等水位線是根據離散的水位值經插值而得,其精度有一定局限。
盡管如此,模擬所得的水位和觀測水位形狀相似,趨勢相近,兩者誤差最大不超過0.50m,地下水流動方向和觀測的流動方向基本一致。模型可以用于預測典型地塊水庫蓄水至90m時地塊的等水位線變化情況。計算參數采用經模型識別調整后的計算值,如表2所示。

表2 龍井~下土城地塊浸沒預測計算采用水文地質參數表
水庫蓄水后,在龍井~上土城形成三面臨水的半島,其地下水流動方向將發生變化。模擬區下游邊界為庫水位定水頭邊界,上游邊界偏安全考慮,仍然為定流量邊界,但其水位將隨庫水位上升有所壅高。取平行于流線的一個剖面,下游邊界水位初始高程h為85m,蓄水后至90m;上游邊界初始水位h1為90.50m,距離庫水邊界約900m,可以近似計算得到庫水位上升后上游邊界的水位。當庫水位上升至90m時,上游邊界水位上升至90.89m,偏安全考慮取91m。蓄水后龍井~下土城地塊模擬區面積為1.60 km2。
降雨入滲作為模擬區上部邊界條件,在模擬初始水位時,由于實測水位的時間為3月下旬至4月上旬,降雨量小,所以在模型識別與校驗時未考慮降雨影響。在預測水庫蓄水后,典型地塊地下水位壅高時,應當考慮降雨入滲的影響。
模型計算時,偏安全考慮取6-8月降雨集中時間,平均日降雨量為2.50mm/d,兩典型地塊地表都以低液限粘土為主,其入滲系數取經驗值0.01,得到降雨入滲量為2.50×10-5m/d。
如前所調整的模型邊界及計算參數,代入模型計算,可以得到地塊模擬區在水庫蓄水以后的等水位線圖,進而作出模擬區的地下水埋深等值線圖。在龍井~下土城地塊模擬區,在蓄水后模擬區的滲流場特征比蓄水前發生了明顯變化。其一是隨著庫水位上移至高程90m,從宏觀上看,全區地下水位上升后的運動總趨勢大致仍由東北向西南流動。但在上土城和下土城將變成三面臨水的半島狀庫岸,地下水呈扇狀流動,分別向淮河和半島兩側的沖溝或低地(都是庫區)排泄。其二是隨著庫水位上升,地下水位發生不同程度的壅高,升幅多在0.40~4.50m。其中以沿庫周邊線附近上升較明顯,向后緣升幅逐漸減小,至后緣上游邊界附近升幅僅40~50 cm。其三是地下水埋深發生不同程度的變淺,最明顯的是沿庫邊線附近地下水埋深<1m的范圍有所擴大。
依據《水利水電工程地質勘察規范》(GB50487-2008),評價浸沒的地下水臨界埋深Hcr=土壤毛管水上升帶的高度(m)+安全超高值(m)。
庫區浸沒影響主要有農作物、居民,農作物主要為小麥、水稻、油菜等,最大根系發育深度約0.60m,一般民房基礎砌置深度為0.50m左右。故安全超高值取0.60m。龍井~下土城地塊Q4低液限粘土毛細上升高度取0.63m,可以得到浸沒的地下水臨界埋深為1.23m,偏安全取1.25m。
通過預測的地下水埋深Hm與浸沒的地下水臨界埋深Hcr進行比較,可判斷模擬區的浸沒程度。當地下水埋深Hm>Hcr時,不會產生浸沒。當Hm≦Hcr時,將產生浸沒危害,即建筑物地基承載力下降,威脅建筑物安全穩定;因過度潮濕使當地群眾生活環境惡化;因鹽漬化危及農作物的正常生長等。其中當Hm<0m時,將產生沼澤化,危害更大。
龍井~下土城地塊蓄水后模擬區地下水埋深在0~3m,地下水臨界埋深Hcr為1.25m,以此可以得到該地塊的浸沒范圍。當庫水位達到90m時,該地塊模擬區浸沒寬度(90米庫水位至浸沒邊界)在60~780m之間不等,90m庫邊線長3300m,浸沒面積約為0.60 km2。受災對象以農田為主,其次沿庫岸邊線的村民和農田也受到一定程度的危害。
影響庫區浸沒的因素較多,如地層巖性、透水性、厚度空間分布、第四系潛水埋深、補給徑流排泄、第四系含水層與其下覆含水層的相互關系、氣候特征、庫水位高低、庫岸地形坡度等。庫區浸沒范圍大小受以上各種要素的綜合作用的影響。根據敏感度系數(SAF)的研究成果,認為庫岸地形坡度是影響浸沒范圍的主要因素。
在準確掌握模擬地塊水文地質條件前提下,MODFLOW可以幫助我們通過三維數值模擬預測水庫蓄水后庫周的地下水位,從而分析研究浸沒發生的范圍和程度。某水庫龍井~下土城地塊通過以上方法對浸沒問題進行了評價,取得了可靠的結果。
[1]賀國平,張彤,趙月芬等.GMS數值建模方法研究綜述[J].地下水,2007(03).
[2]武強,朱斌,徐華等.MODFLOW在淮北地下水數值模擬中的應用[J].遼寧工程技術大學學報,2005(04).
[3]陳全禮,張志敏.庫區蓄水浸沒預測評價[J].資源環境與工程,2008(B10).
[4]張穎,許模,王子忠.應用三維地下水滲流數值模擬方法分析水庫浸沒問題[J].四川水利,2008(04).
[5]GB50487-2008,水利水電工程地質勘察規范[S].
