基于Visual MODFLOW的黑龍江省巴彥縣地下水流數值模擬

婁安宇，陳 末，臧珊珊，吳志剛

(黑龍江大學 a.水利電力學院; b.寒區地下水研究所, 哈爾濱 150080)

0 引 言

隨著社會的發展，巴彥縣的城鎮化建設也迅速發展，地表水資源已經供不應求，地下水資源的利用在該地區已達到不可替代的地位。不過，地下水不是無窮無盡的，所以要對該地區的地下水進行分析研究。地下空間錯綜復雜，地下水含水層中有著各種形狀復雜的空隙[1]。因此，分析研究區地下水的補給來源、排泄途徑，進行邊界條件的概化，構建研究區地下水流概念模型，從而得到相關的水文地質參數是一個極為簡單的方法[2]。

數值模擬方法在地下水研究方面的效率、靈活性和低成本管理等方面已成為一種重要的方法，并日益引起注意和廣泛的應用[3]。數值模擬方法是一種數值模擬程序的微分方程，根據邊界和初始基本條件，來描述地下水水位隨著時間的推移而產生連續變動得數學模型,并且通過電腦軟件中的離散化方法來解決這個數學問題[4]。

本文應用數值模擬方法，通過使用Visual MODFLOW 軟件建立巴彥縣地下水流數值模型，得到該地區相應的水文地質參數，并對這一地區地下水水位進行預測。

1 研究區概況

巴彥縣位于黑龍江省西部偏南，隸屬于哈爾濱市，巴彥縣境南北長85 km、東西寬72.7 km，幅員3 137.7 km2。境內主要有“一江三河”，其中“一江”是松花江，“三河”分別為泥河、少陵河和五岳河。見圖1。

圖1 研究區示意圖

巴彥縣地區屬于中溫帶大陸性季風氣候，每個季節的溫差很大，溫度最高可以達到36.5℃，最低可以到-40.2℃，年平均氣溫約為2.9℃[5]。地下水開采量達12 734.19 m3，主要為農業灌溉用水，部分用于鄉鎮生活用水、工業生產和生態環境改善等方面。全縣地下水資源開發利用量約為7 400×104m3[6]。

地下水貯存條件及分布規律為，一級階地區含水層分布均勻，厚度穩定，堆積物松散，水量較豐富，地下水類型為階地松散層孔隙承壓水。漫灘區孔隙潛水在2月下旬至3月初為枯水期，進入雨季后水位逐漸上升，至8-9月初上升到最高值，水位年變幅1～4.5 m。階地區孔隙承壓水在4月下旬至5月上旬為枯水期，9月末至10月初出現最高值，地下水動態類型為降水入滲-徑流型[7]。

綜合全區地貌成因和形態，將巴彥縣劃分為3種成因類型，3個形態單元見表1。

表1 地貌單元分區表

2 地下水流數值模型

2.1 水文地質概念模型

以相對完整的水文地質單元為數值模擬計算區，盡量將計算區邊界設置在自然邊界處，或者設置在容易確定流量或地下水位的人為邊界處。本次研究區位于哈爾濱市巴彥縣，北部邊界為泥河，西部邊界為漂河和少陵河，南至松花江干流。位于松嫩平原東部，松花江中游左岸，依據地形地貌特征、地層結構及水文地質特征，將區內地下水劃分為河漫灘松散層孔隙潛水、一級階地松散層孔隙承壓水和高平原(松花江二級階地)松散層孔隙承壓水。隨著系統的變化，三維空間參數向量的一般地下水流動速度反映了輸入和輸出系統的異質性，輸入和輸出系統隨時間和空間的變化而變化，來顯示地下水的流動情況。

2.2 含水層系統的概化

該研究區的地下水主要是以水平運動為水流的運動方向，而垂直方向的運動則很少，水文地質參數伴隨空間分布的變化而改變，這反映了地下水系統的均質性。地下水系統輸入和輸出隨時間和空間的變化而變化，地下水流是不穩定的[8]。綜上所述，將該模擬區的含水系統概化為均質的各向同性的三維非穩定流。

2.3 模型邊界條件的概化

地下水流向主要從泥河到松花江，泥河、大荒溝、少陵河以及漂河的流向切割了該區域的含水層并且地下水與它們有補排關系。當河流或湖泊切割含水層，兩者有直接水力聯系時，可作為第一類邊界條件處理[9]，因此將泥河、大荒溝、少陵河以及漂河作為定水頭邊界。由于其他邊界和地下水流方向平行，所以把其定為地質零通量邊界，沒有側向流量的補給，所以作為隔水邊界。

2.4 數學模型

數學模型描述了一個研究領域即地下流動的水流，并確定數學方程的解決方案，這些條件構成一個真正的數學結構問題。含水系統概化為均質的各向同性的三維非穩定流，其相應的偏微分方程以及定解條件如下：

式中：h為水位，m；t為時間，d；W為源匯項，1/d；K為滲透系數，m/d；Ω為模型模擬區；Γ1為第一類邊界；Γ2為第二類邊界；n為邊界面的外法線方向；H0為地下水初始水位函數；q(x,y,z,t)為第二類邊界上已知流量函數。

3 模型的識別與驗證

3.1 研究區的剖分

對于研究區的剖分采用有限差分方法求解，并且采用強隱式法聯立迭代求解代數方程組[10]，在水平方向上對潛水含水層用相互垂直的平行線對研究區進行網格剖分，對研究區剖分成100×100的矩形單元格，其中有效單元計算點共3 164個。

3.2 計算時限與步長

識別期為2015年1月1日至2016年12月31日，驗正期是2017年1月1日至2017年12月31日。由于地下水位受地下水開采與利用的影響，所以在對模型進行輸入參數時將每個季度作為一個地下水開采期。

3.3 參數的選取

3.3.1 滲透系數K

根據巴彥縣當地的水文地質條件以及含水層巖性，將其劃分為3個不同滲透系數值的區域，見圖2。通過查閱資料以及查看巴彥縣的水文地質圖可以發現，巴彥縣的含水層巖性基本為亞黏土、粉土質砂及粉砂等地質；根據水文地質手冊可以得到，所有分區的滲透系數經驗值的區間為0.5～2 m/d。不同的分區輸入相應的滲透系數初始值，由于此次模型屬于各向同性，所以水平滲透系數KX、KY與垂向滲透系數KZ相同。

圖2 研究區滲透系數分區圖

3.3.2 降雨入滲系數

通過查找資料可知，該地區是含水層研究的重要組成部分，橫向滲透、邊境供應和其他形式的補給作用相對較弱，且觀測井的埋深都在1～5 m之間，所以判斷巴彥縣的地下水補給主要為大氣降水補給。根據巴彥縣水文地質條件、流域水系、行政分區現狀、地貌和巖性等特征，把研究區的大氣降水滲透性分為6個區。見圖3。

圖3 研究區大氣降水滲透性分區圖

3.3.3 初始水頭

此次模型中初始水頭為各個觀測井中2015年1月1日的水頭值。觀測井位置見圖4。

圖4 觀測井位置圖

3.4 模型的識別

模型采用間接求參法達到識別目的, 通過不斷調整參數誤差的方法對測井計算的初值，使其計算水頭和測量水頭符合擬合要求, 即計算水頭和實測水頭之間相差很小。若計算水頭和實測水頭相差較大，說明設定參數或邊界條件不合適，需調整，只能通過模型識別、調整、再識別、再調整的多次反復過程，一直到計算水位與實測水位誤差滿足要求為止，此時所選用的參數即為校正后的參數。觀測井校準水位圖見圖5，誤差分析見表2，觀測井水頭時間圖見圖6，單個觀測井水位擬合圖見圖7(以5號觀測井為例)。

圖5 觀測井水位校準圖

圖6 觀測井水頭時間曲線

圖7 5號觀測井水位擬合曲線

表2 誤差分析表

通過以上圖表可以發現，各個觀測井的計算水位與觀測水位變化已基本一致，除個別點處的水位變化步調有偏差外，大部分點的變化基調基本上是同步的。通過該圖可知，該模型具有良好的擬合效果，說明建立的數學模型、各類邊界條件概化、水文地質參數的選擇、源匯條件的處理是合理的，基本反映了研究區內地下水真實的變化特征[11]。

3.5 模型的驗證

為使建立的數值模型更為準確地描述出地下水流系統的真實情況，需要利用現有的另一個時段的地下水動態觀測資料對已經建立的模型進行檢驗，以此來進一步證明所建立的數值模型和所得模型參數的可靠性，本次驗證的初始時間為2017年1月1日。觀測井校準水位圖見圖8，誤差分析見表3，觀測井水頭時間圖見圖9，單個觀測井水位擬合圖見圖10(以5號觀測井為例)。

圖8 觀測井水位校準圖

表3 誤差分析表

圖9 觀測井水頭時間曲線

圖10 5號觀測井水位擬合曲線

從模型識別結果和擬合效果可以看出，地下水位擬合程度非常高，建立的模型可以真實地反映該地區實際的地下水情況，同時也論證其數學模型、邊界條件、水文地質參數和源匯條件建立的正確性。

3.6 水文地質參數

經過模型的識別與驗證可以得到，研究區各個分區的滲透系數和降雨入滲系數值見表4、表5。

表4 分區滲透系數表

表5 分區降雨入滲系數表

4 地下水位預測

本文建立的數值模型通過不同監測井的結果經過分析與計算得到含水層的含水量變化情況，并根據線性回歸方法對研究區地下水位進行預測。大多數情況下，僅有流經方向和降雨量等一些自然環境影響地下水位。將模型初始時間設定為2015年1月1日，不改變已經建立的模型的任何條件與參數，對未來10年的地下水位進行預測。研究區未來10年地下水位預計變化情況見表6,地下水位折線圖見圖11，單個觀測井地下水位折線圖見圖12-圖14。

表6 地下水位預計變化情況

圖11 地下水位預測折線圖

圖12 4號觀測井水位預計變化圖

圖13 5號觀測井水位預計變化圖

圖14 3號觀測井水位預計變化圖

根據表6、圖11總體來看，巴彥縣地下水位將會逐年下降，造成這一現象的根本原因是地下水開采量大于補給量。巴彥縣經濟主要以農業生產為主，是黑龍江省的糧食主產區，由于地表水的供水能力有限，需要不斷開采地下水來維持農業的發展。地下水開采主要用途是農業灌溉，部分用于鄉鎮生活用水、工業生產和生態環境改善等方面。

根據圖12-圖14可以看出，巴彥縣不同地區地下水位的變化情況并不相同，大部分地區地下水位有下降趨勢。4號井位于洼興鎮附近，地下水不僅用于農業灌溉，附近的城鎮居民生活用水以及工業用水也需要通過地下水進行補給，所以位于城鎮附近的地下水位下降很明顯。5號井位于村莊附近，地下水基本只用于農業灌溉，地下水位只是略有下降。3號井位于少陵水庫和漂河附近，能夠得到足夠的地表水補給，這一地區的地下水位會有些上升。

針對地下水過度開采的情況，應該嚴格控制地下水的開采量，可以通過增加工廠企業的中水回用，取締新水，逐漸降低對地下水的依賴。還可以考慮通過建設各種大中小型蓄水工程，在汛期時，對過境洪水進行攔蓄使用。同時，相關部門應加強監測站網的建設，進一步積累和完善監測資料，提高地下水資源評價成果的可靠性和精度。

5 結 論

本文以巴彥縣為研究區域，根據前期水文地質資料的搜集、水文地質勘察與試驗結果，建立了水文地質概念模型， 并運用Visual MODFLOW軟件進行三維地下水流模擬，經過模型識別與校正，建立了符合目標的地下水流模型。得到以下結論：

1) 農業用水主要來源是地下水，同時隨著城市規模的不斷擴大，人口增長和經濟建設的迅猛發展使得工業和生活用水量日益增大，區域地下水資源供需矛盾日益突出。

2) 通過建立巴彥縣地下水流數值模型，經過模型的識別與驗證，可以確定巴彥縣地區的滲透系數為：河漫灘地區0.8 m/d；山區1.4 m/d；山前地區1.7 m/d。