基于HYDRUS-MODFLOW的大洋河岫巖流域地表水與地下水耦合模型

張 賀

（遼寧潤中供水有限責任公司,遼寧 沈陽 110166）

1 研究區概況

大洋河發源于岫巖縣偏嶺鄉境內的一棵樹嶺南側，流經岫巖縣、東港市，在東溝縣的黃土坎入黃海[1]，干流全長230 km，主要支流有雅河、牤牛河、連河、哨子河等，全部流域面積6004 km2，多年平均徑流量31 億m3。其中，大洋河岫巖縣境內長度180.2 km，流域面積1968.4 km2，是岫巖縣的重要飲用水水源地[2]。大洋河岫巖流域位于遼東半島東部，地貌以低山丘陵為主，呈現出西北高，東南低的特點。研究區屬于溫帶濕潤性季風氣候，多年平均氣溫為8 ℃，多年平均降水量為767 mm，多年平均蒸發量為1066 mm。由于研究處于距離黃海較近的東部迎風坡，因此受季風影響比較顯著，汛期的7 月、8 月集中了全年6 成以上的降水[3]。研究區的地下水主要是松散巖類孔隙水，且主要富集于全新統地層中，沿河谷和山間谷地均有分布，成因主要是沖積、沖洪積、坡洪積等。根據地質鉆孔和民井調查資料，含水層分布在階地和河漫灘中，上覆黃褐色亞砂土層厚度為1.0 m～1.5 m，粗砂和砂礫石互層厚度為2 m～3 m，底部的砂礫卵石層厚3 m～5 m，含水層總厚度為4.5 m～6.5 m。研究區的地下水主要是側向徑流補給、大氣降水補給、農田灌溉水滲入補給以及開采條件下的河流入滲補給。隨著經濟社會的發展，研究區的水資源需求量逐年攀升，而氣候變遷的大背景下，研究區的水資源總量呈遞減狀態。因此，地下水資源的不斷開發利用將成為研究區的必然趨勢。所以，加強地下水資源研究，對研究區水資源的合理開發和利用具有重要意義。

2 耦合模型的構建

2.1 構建思路

流域內的水循環是一個十分復雜的運動過程，而地表水和地下水兩種主要水資源形式的賦存和運動狀態之間也存在顯著差異[4]。基于水循環的復雜性以及地表水和地下水循環上的差異，當前的流域水循環的研究主要為地表水或地下水的單獨模擬。然而，上述研究模式中水循環過程缺失是不可避免的，研究成果的精度也難以貼合實際。因此，要獲得較高精度的模擬成果，需要將地表水與地下水這兩種主要水資源形式進行綜合研究[5]。本文研究利用地下水模型MODFLOW模型，結合研究區地表水文資料進行參數率定，在校驗準確并滿足模型精度要求的基礎上（限于篇幅這里不再詳述），將其輸出項作為Visual MODFLOW的輸入項對研究區地下水進行數值模擬，試圖通過地表水和地下水的耦合模擬，提高模擬精度[6]。

2.2 模型準備

大洋河岫巖流域的孔隙水在天然狀態下的水力坡度不大，滲流狀態也基本符合達西定律。因此，地下水的水流形式可以概化為平面二維流，同時在計算時段內表現為非穩定流特征。研究區的地下含水層可以概化為非均質各向同性的潛水層，且下部直接與不透水黏土層相接。基于此，研究區的南北方向邊界設定為零流量邊界條件，東西方向為GHB通用水頭邊界[7]。模型的計算時間為2006 年1 月～2015 年12 月，地下水模型MODFLOW的計算步長和應力期均設定為1 d，全部計算時段共3650 個應力期。非飽和帶模型HYDRUS的初始時間步長設定為0.5 d，最大和最小時間步長分別為1 d和0.3 d。研究區地下水的初始水位采用2006 年1 月1 日的觀測井實測數據。

2.3 網格剖分

在HYDRUS-MODFLOW地表水和地下水耦合模型時，首先利用MODFLOW模型對大洋河岫巖流域進行空間離散。網格單元在X方向和Y方向上的尺寸分別為974.55 m和913.05 m，最終劃分為2678 個計算單元。在豎直方向上，根據觀測井的實測地下水位數據，確定出含水層的飽和帶與非飽和帶的厚度。然后，利用HYDRUS模型，在研究區的水平方向上，根據下墊面以及土壤特性，研究區劃分為150 個土壤單元體，每個土壤單元按照不同深度部位的特征劃分為3 個～38 個計算節點，進而實現土壤單元體的有限元網格劃分[8]。

2.4 模型輸入

由于MODFLOW是具有模塊化結構的模型，而基于HYDRUS-MODFLOW地表水和地下水耦合模型的基本思路是將HYDRUS子程序嵌入到MODFLOW主程序中[9]。因此，在耦合模型中選擇最適合的模塊進行模擬。結合本次研究的實際情況，增加河流子程序以及地表產匯流過程，以實現對研究區的地表水文過程的模擬。

2.5 參數的率定

根據研究區土壤的實際情況，對土壤類型進行概化分為亞砂土和亞粘土兩個類型，其水力參數的初始值根據相關研究成果確定，然后利用人工試錯法對參數進行反復調整，直至模型的模擬和實測結果之間具有良好的擬合度，率定后的模型參數見表1。

表1 參數率定結果

3 計算結果與分析

3.1 地下水位擬合結果分析

由于設置的模型模擬應力期長度為1 d，共3650 個應力期，數據量較大，因此在統計數據過程中選擇以月為時間單位，共120 個月。利用模型的模擬結果和觀測井的實測值，對研究區所有的125 眼觀測井的地下水位的平均值進行對比，獲得圖1 和圖2 所示的模擬和實測地下水位過程和線性相關圖。由圖1 可知，研究區的地下水位的模型模擬結果和實測值具有基本一致的走向變化，具有較好的擬合度。由圖2 可知，研究區地下水位的模擬值與實測值的水位散點基本都落在直線附近，說明兩者之間的相關性較好，計算兩者之間的相關系數，結果為0.79。由此可見，整個流域的實測值與模型的模擬值之間具有較好的吻合度。

圖1 地下水位模擬值與實測值關系圖

圖2 模擬和實測地下水位散點圖

鑒于研究區內地表水和地下水的補排關系十分密切，是研究區地下水位的重要影響因素，因此距離河道較近區域的地下水位受河流徑流影響比較顯著，變化也比較劇烈。因此，研究中選擇位于河道附近的嶺溝46#、39#、12#、23#、17#觀測井，哨子6#、14#、17#、24#、31#觀測井等10個代表性觀測井進行地下水位模擬，并與實測值進行對比。其中，嶺溝46#、39#觀測井和哨子6#、14#觀測井的地下水位過程模擬和實測結果見圖3 ～圖6 。由模擬結果可知，大部分觀測井的模擬結果和實測結果之間具有良好的擬合度，僅有個別監測井部分時段的模擬誤差較大。究其原因，主要是部分地區的土壤類型和下墊面變化過于復雜，地下水的開采量難以進行準確估算，同時枯水年份地下水位受到灌溉用水的影響比較大等。

圖3 嶺溝46#觀測井模擬和實測地下水位過程

圖4 嶺溝39#觀測井模擬和實測地下水位過程

對上述10 個觀測井的納什效率系數和模擬的相對誤差進行統計，結果見表2。由表中的結果可知。各個典型觀測井的納什效率系數的值均在0 附近波動，且波動幅度不大，說明各個觀測井的模擬結果和實測結果比較接近。從相對誤差的統計計算結果來看，各觀測井的相對誤差值在3%以內，也說明模型的模擬程度較好，可以用于研究區的地下水流動過程的模擬計算。

表2 各觀測井的誤差情況分析結果

3.2 出口斷面流量過程計算結果分析

為了進一步驗證模型的模擬能力，對研究區出口斷面的流量進行模擬計算，以便對模型的模擬效果進行驗證。圖7為根據模擬值和實測值繪制的逐日流量過程圖，由圖7 可知，模擬值和實測值之間具有良好的擬合度。具體來看，在參數的率定期和驗證期，納什效率系數分別為0.81 和0.79，徑流量相對誤差分別為3.66%和3.85%，整個模擬期的納什效率系數為0.79，徑流量相對誤差為0.44%，各個階段和整個模擬期的納什效率系數和徑流量相對誤差均在允許范圍內。

鑒于模擬期間的1999 年、2000 年以及2001 年，由于降水量偏大造成徑流量變幅較大。因此對選取其中變幅最大的2000 年，進行逐日流量的模擬值和實測值對比，結果見圖8。由圖8 可知，模型的模擬值和實測值之間具有良好的擬合度，納什效率系數為0.88，相對誤差為6.23%，均在允許范圍內，說明模型在出口斷面流量變幅較大的情況下，亦可以獲得良好的模擬效果。

圖7 模擬與實測逐日流量過程

圖8 2000年模擬與實測逐日流量過程

4 結論

以遼寧省大洋河岫巖流域為例，對基于HYDRUSMODFLOW的地表水和地下水耦合模型進行研究，獲得如下結論：

（1）研究區的地下水位的模型模擬結果和實測值具有基本一致的走向變化，兩者之間的相關系數為0.79，具有較好的擬合度。

（2）研究中選擇位于河道附近的10 個代表性觀測井進行地下水位模擬，并與實測值進行對比。由模擬結果可知，各個典型觀測井的納什效率系數值均在0 附近波動，相對誤差值在3%以內，因此模擬結果和實測結果之間具有良好的擬合度。

（3）對研究區出口斷面的流量進行模擬計算，結果顯示整個模擬期的納什效率系數為0.79，徑流量相對誤差為0.44%；對徑流量變幅較大的2001 年的逐日流量的模擬值和實測值進行對比，結果顯示模擬值和實測值之間具有良好的擬合度，納什效率系數為0.88，相對誤差為6.23%，均在允許范圍內。

（4）綜合上述，基于HYDRUS-MODFLOW的地表水和地下水耦合模型具有較高的模擬精度和適應性，可以用于研究區的相關研究。