基于有限元法計算三江平原地下水開發潛力研究

宮興龍，付 強，李 衡，郭 佳(東北農業大學水利與建筑學院，哈爾濱 150030)

0 引 言

從20世紀90年代以后隨著墾區興起打井種稻,水稻的種植技術取得了長足的發展，產量穩步提高，同時由于打井大量開采地下水，降低了地下水位，澇漬堿危害減小。但人們擔心,抽取的地下水水溫低，影響著影響水稻米質與產量，在不斷擴大水井稻的同時需要不斷加大抽取地下水，導致地下水位不斷下降，使這些地區可能面臨超采，經濟資源環境發展難可持續[1-3]。故本文僅重點論證地下水資源開發能否支撐本農場糧食生產、經濟發展和生態環境可持續問題,即地下水資源開發潛力分析。很多人進行了三江平原地下水開發潛力研究，但大部分都是基于水量平衡法[4-9]，基于此本文嘗試通過分析三江平原灌區典型區域853農場的水文地質條件，采用Galerkin法推導出該地區承壓水的有限元方程組，并采用了改進平方根法求解該方程組，進而計算出853農場的地下水開發潛力。

1 地下水流數學模擬

承壓水地下水運動采用的數學模擬模型如下[8]:

(1)

H(x,y,0)=H0(x,y)t=0,(x,y)∈Q

(2)

H(x,y,t)|Γ1=φ1(x,y)t=0,(x,y)∈Q

(3)

(4)

式中:T為導水系數；S為貯水系數；H0為水頭初始值；φ1為Γ1(第一類邊界)的已知水位函數；q為Γ2(第二類邊界)上的側向補給量；n為邊界Γ1的外法線方向；Q為方程建立區域。

1.1 基于Galerkin法構建地下水運動有限元方程

對式(1)～式(4)采用Galerkin法得到有限元方程為：

(5)

式中：φi為基函數；n為除了第一類邊界上的節點外的節點數。

應用分布積分法(5)式得：

(6)

式(6)左端第一項應用Green公式:

(7)

利用式(7)對式(6)整理：

(8)

以三角元離散式(8)：

(9)

在單元上式(9)變形為：

(10)

式中:e表示三角單元上的變量。

1.2 有限元方程的構建

為了構建有限元方程對(10)變形為：

(11)

式中:Δe為三角元面積；qai,qik為邊界段的流量；Lai,Li為邊界段長度；a,b,c采用下式計算：

(12)

式中:(xi,yi)、(xj,yj)、(xk,yk)為三角元三個節點i,j,k的坐標。

1.3 方程組求解

式(12)利用矩陣形式為：

(13)

式中：導水系數矩陣D、貯水系數矩陣P、常數項矢量矩陣F。上述方程是具有對稱，正定系數矩陣的線性方程，因此采用改進平方根法求解[9,10]。

1.4 抽水井有限元方程

“853”農場的取得地下水主要方式是利用抽水井進行抽水，如果抽水井在三角單元節點上應在三角單元方程(13)式上加一項：

(14)

式中：rwi為水井的半徑。

如抽水井在單元的其他段，單位時間，單位面積上的垂直流出(入)含水層的量為ω=-Q/Fw(Fw單元面積)，單元e上的節點i,j,k依次有：

(15)

2 計算實例

2.1 研究區地質概況

853農場位于三江平原東部，寶清縣境內，屬于溫暖半濕潤農業氣候區，多年平均降水量557.2 mm。在地理位置上，東南部緊靠完達山，西北及北部以平原河流撓力河為界，西南的以生產蛤蟆的蛤蟆通河為邊界。本區有基巖裂隙孔隙含水層、第四系孔隙承壓含水層、第三系裂隙承壓含水層,第四系孔隙含水層系統在區內最重要，含水層較厚供水量充沛,從山前向河谷含水層逐漸增厚，地下水埋深逐漸變淺。第三系裂隙孔隙水系統內流動較弱，承壓水的賦存、運輸開采條件差于第四系，富水性差于第四系。該區農業灌溉用水主要來自于第四系含水層，該層地下水易于補給和排泄，因該區內地下水的開采還可能深入到第三系含水層,但開采量比較少，基于此本次模擬以第四系含水層為主。該區的主要補給方式是河流滲漏、降水入滲和側向地下徑流；人工抽水井開采是地下水排泄的主要方式,在開采量小的時間里還有一部分河流排泄和側向徑流排泄。

2.2 研究區域離散

為建立研究區的地下水流運動介質中的水文地質模型，基本步驟是要對實際水文地質條件加以概化，建立水文地質概念模型。根據上述水文地質特殊條件,本文主要致力于研究第四系含水層，山區第四系覆蓋較薄，僅有少量孔隙裂隙潛水可供利用，其貯存于完達山麓，本次模擬不計。第三系裂隙孔隙承壓水儲存于巨厚的第四系含水層之下，在目前的開采條件和經濟條件下尚難以對其進行開采，且收集資料困難，所以本次模擬也不考慮。由于第四系含水層其覆蓋有10～35 m厚的亞黏土層,最厚達42 m，具有良好的隔水能力，地下水具有承壓性質。故將含水層在垂向上概化為1層，即第四系孔隙承壓含水層，其上腹亞黏土，下部為不透水基巖。水平三條河流概化為流量邊界。采用有限元法數值求解式(7)對上述數學模型進行求解，對比例尺為1：10萬的八五三農場綜合水文地質圖及地形圖進行掃描，然后用數字化工具進行矢量化處理為柵格文件,將柵格文件導入River Tool工具，提取出撓力河，蛤蟆通河流，七里沁河和圖1右下方的邊界。模型作為計算模擬區的剖分底圖計算區面積為1 182 km2，采用1 890個三角元進行離散見圖1。圖1上邊的左邊邊線為撓力河，右邊的線為七里沁河，下邊的線為不透水邊界，左下邊的線為蛤蟆通河。

圖1 853農場三角元離散圖Fig.1 Discrete map of Triangular element in 853 farm

2.3 水位觀測點

853農場有長期地下水觀測井，本文選取7個長期觀測井，分別為853農場1分場6小隊的觀井編號32號，2分場1小隊的觀井編號6號3分場7小隊的觀井編號10號，4農場1小隊的觀井編號25號，5農場2小隊的觀井編號20號，6農場5小隊的觀井編號14號，6農場9小隊的觀井編號13號。

2.4 降水量

降雨通過分析1995-2001年的月降水量數據，分析得到最大降水年份是1959年降水量為729.6m，多年平均降水量523.35。最小降水年份2004年458 mm。其過程線見圖2。可見本地區降水量主要集中在6-8月份。本地區的降水入滲系數空間分布難于確定，因此本文根據清河灌區的設計值采用0.06。

圖2 月降水量數據圖Fig.2 Monthly rainfall data map

2.5 地下水開采控制水位判斷

文獻[4]還指出，根據水文地質條件，開采技術條件和社會生產需要，利用地下水變化直接反映開采程度的特點，提出合理的“地下水開采控制水位”概念，防止超量開采地下水導致地下水枯竭。筆者暫提出“理想水位(埋深)”概念，853農場承壓開采為主區的理想水位本文控制在10～15 m。

2.6 目前開采條件下地下水開發程度

采地下水一直存在爭議，本文以“853”農場的1998年到2012的觀測井數據進行分析。以“853”農場的5農場2小隊的觀井編號20，因為其在清河灌區內，存在開采地下水進行水稻田灌溉的量大于其他各觀測井所在的地區。從圖可以看出地下水下降最大的一年是2012年，地下水為45.712 7 m。20號觀測井的地面高程為58.912 7，其差值為13.5，由于其為承壓區不超過15 m，因為本位認為853農場目前開采地下水不存在超采，而且還有一定開采的空間。由八五三農場五分場2隊地下水位高程曲線圖(圖3)可以看出其地下水恢復地比較快，可見地下水很豐富,故有很大開發潛力。

圖3 5分場2隊地下水位高程曲線圖Fig.3 The underground water level elevation curve in fifth parvial field second team

2.7 基于數值法推求地下水開發潛力研究

“853”農場從近年地下水動態分析，每年4月恢復水位埋深大約3.0 m，尚小于理想埋深最高值(10 m)。“853”農場屬于以承壓水為主的地下水開采由。撓力河河道寬廣，一年時間里水深變化不大，在56.5 m波動，本文據此給出撓力河水位為56.5 m。蛤蟆通河水位較淺，本文采用的水位為河道高程，同理對于七里沁河也采用同樣處理。對7個水位觀測點的多年1月份的地下水水位分析得出其多年平均初始水位，采用距離倒數法插值得到計算所用的初始水位(見圖4)。第四系承壓含水層由于資料較少，其參數的選取主要是根據含水層水文地質特征,及前人的抽水試驗,參照經驗值給出初始值彈性釋水系數0.003 5，導水系數1 560 m2/d。控制水位采用10～15 m。為了計算地下水開發潛力本例按照承壓水理想埋深10～15 m計算地下水開發潛力，地下水為控制在水深10 m是計算的開發潛力為8 328.09萬m3/a、地下水水位控制在15 m時，可得其開發潛力為9 456.09萬m3/a。

圖4 853農場多年春季平均水位圖Fig.4 The average water level diagram of spring for many years in 853 farm

2.8 激發出江河對地下水補給

“853”農場周圍三面鄰河，北及北部以撓力河為界，西南的蛤蟆通河，西北的撓力河，東北的越嶺河，因此在控制開采地下水水位降到一定程度，就發生江河加大補給地下水[6]，基于此本位在控制承壓區10～15 m，計算其激發的河流補給量。由于蛤蟆通河和七里沁河水位相對于承壓水比較高，激發補給量不明顯。因此本文的激發江河補給主要考慮撓力河，由計算可得，在控制降深水位取15 m，可得其激發的補給量為256.09 萬m3/a。

3 結 語

綜上所述，由“853”農場從1998-2013的觀測井數據及地下水水位控制數據分析得出，目前“853”農場地下水開采不僅沒有超采，而且還有很大開發潛力，地下水位下降只是7、8月份，之后隨著地下水開采量減少，地下水水位迅速的恢復。由基于有限元法構建的地下水運動數值模型計算得到“853”農場可以開發9 456.09 萬m3/a水量，“853”農場的北部靠近撓力河和西部的蛤蟆通的地區地下水開發量明顯大于東部的靠近七里沁河的開采量。在合理控制地下水情況下，可以激發出江河更多的補給，撓力河中部補給量為256.09 萬m3/a。

□

[1] 王立權,馮建維.水稻熱的隱憂[J].水利天地,2004,(6):18-19.

[2] 張惠斌,于 東,姚章村.論“打井種稻”與“循環經濟”[J].水利科技與經濟,2006,12(12):820-821.

[3] 王立權,姚章村.高寒地區井灌水田井水增溫技術的探索與實踐[J].地下水,2006,28(3):46-49.

[4] 何 璉.中國三江平原[M].哈爾濱:黑龍江科學技術出版社,2000.

[5] 姚章村,安瑞強,董經財,等.三江平原建三江地下水激發補給初步分析[J].黑龍江水專學,2001,28(3):20-23.

[6] 李宏勛,陳 杰,李 偉,等.291農場地下水資源開發潛力初步分析[J].水利科技與經濟,2008,14(1):49-51.

[7] 王貴玲,劉志明,劉花臺,等.地下水潛力評價方法[J].水文地質工程地質,2003,(1):63-66.

[8] 薛禹群,謝春紅.地下水數值模擬[M]. 北京：科學出版社,2007.

[9] 張德良.計算流體力學教程[M]. 北京：高等教育出版社,2009.

[10] 吳頌平,劉趙淼 譯.計算流體力學基礎及其應用[M]. 北京：高等教育出版社,2013.