張 賀

（遼寧潤(rùn)中供水有限責(zé)任公司,遼寧 沈陽(yáng) 110166）

1 研究區(qū)概況

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

2 耦合模型的構(gòu)建

2.1 構(gòu)建思路

流域內(nèi)的水循環(huán)是一個(gè)十分復(fù)雜的運(yùn)動(dòng)過(guò)程，而地表水和地下水兩種主要水資源形式的賦存和運(yùn)動(dòng)狀態(tài)之間也存在顯著差異[4]?；谒h(huán)的復(fù)雜性以及地表水和地下水循環(huán)上的差異，當(dāng)前的流域水循環(huán)的研究主要為地表水或地下水的單獨(dú)模擬。然而，上述研究模式中水循環(huán)過(guò)程缺失是不可避免的，研究成果的精度也難以貼合實(shí)際。因此，要獲得較高精度的模擬成果，需要將地表水與地下水這兩種主要水資源形式進(jìn)行綜合研究[5]。本文研究利用地下水模型MODFLOW模型，結(jié)合研究區(qū)地表水文資料進(jìn)行參數(shù)率定，在校驗(yàn)準(zhǔn)確并滿足模型精度要求的基礎(chǔ)上（限于篇幅這里不再詳述），將其輸出項(xiàng)作為Visual MODFLOW的輸入項(xiàng)對(duì)研究區(qū)地下水進(jìn)行數(shù)值模擬，試圖通過(guò)地表水和地下水的耦合模擬，提高模擬精度[6]。

2.2 模型準(zhǔn)備

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

2.3 網(wǎng)格剖分

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

2.4 模型輸入

由于MODFLOW是具有模塊化結(jié)構(gòu)的模型，而基于HYDRUS-MODFLOW地表水和地下水耦合模型的基本思路是將HYDRUS子程序嵌入到MODFLOW主程序中[9]。因此，在耦合模型中選擇最適合的模塊進(jìn)行模擬。結(jié)合本次研究的實(shí)際情況，增加河流子程序以及地表產(chǎn)匯流過(guò)程，以實(shí)現(xiàn)對(duì)研究區(qū)的地表水文過(guò)程的模擬。

2.5 參數(shù)的率定

根據(jù)研究區(qū)土壤的實(shí)際情況，對(duì)土壤類型進(jìn)行概化分為亞砂土和亞粘土兩個(gè)類型，其水力參數(shù)的初始值根據(jù)相關(guān)研究成果確定，然后利用人工試錯(cuò)法對(duì)參數(shù)進(jìn)行反復(fù)調(diào)整，直至模型的模擬和實(shí)測(cè)結(jié)果之間具有良好的擬合度，率定后的模型參數(shù)見表1。

表1 參數(shù)率定結(jié)果

3 計(jì)算結(jié)果與分析

3.1 地下水位擬合結(jié)果分析

由于設(shè)置的模型模擬應(yīng)力期長(zhǎng)度為1 d，共3650 個(gè)應(yīng)力期，數(shù)據(jù)量較大，因此在統(tǒng)計(jì)數(shù)據(jù)過(guò)程中選擇以月為時(shí)間單位，共120 個(gè)月。利用模型的模擬結(jié)果和觀測(cè)井的實(shí)測(cè)值，對(duì)研究區(qū)所有的125 眼觀測(cè)井的地下水位的平均值進(jìn)行對(duì)比，獲得圖1 和圖2 所示的模擬和實(shí)測(cè)地下水位過(guò)程和線性相關(guān)圖。由圖1 可知，研究區(qū)的地下水位的模型模擬結(jié)果和實(shí)測(cè)值具有基本一致的走向變化，具有較好的擬合度。由圖2 可知，研究區(qū)地下水位的模擬值與實(shí)測(cè)值的水位散點(diǎn)基本都落在直線附近，說(shuō)明兩者之間的相關(guān)性較好，計(jì)算兩者之間的相關(guān)系數(shù)，結(jié)果為0.79。由此可見，整個(gè)流域的實(shí)測(cè)值與模型的模擬值之間具有較好的吻合度。

圖1 地下水位模擬值與實(shí)測(cè)值關(guān)系圖

圖2 模擬和實(shí)測(cè)地下水位散點(diǎn)圖

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

圖3 嶺溝46#觀測(cè)井模擬和實(shí)測(cè)地下水位過(guò)程

圖4 嶺溝39#觀測(cè)井模擬和實(shí)測(cè)地下水位過(guò)程

對(duì)上述10 個(gè)觀測(cè)井的納什效率系數(shù)和模擬的相對(duì)誤差進(jìn)行統(tǒng)計(jì)，結(jié)果見表2。由表中的結(jié)果可知。各個(gè)典型觀測(cè)井的納什效率系數(shù)的值均在0 附近波動(dòng)，且波動(dòng)幅度不大，說(shuō)明各個(gè)觀測(cè)井的模擬結(jié)果和實(shí)測(cè)結(jié)果比較接近。從相對(duì)誤差的統(tǒng)計(jì)計(jì)算結(jié)果來(lái)看，各觀測(cè)井的相對(duì)誤差值在3%以內(nèi)，也說(shuō)明模型的模擬程度較好，可以用于研究區(qū)的地下水流動(dòng)過(guò)程的模擬計(jì)算。

表2 各觀測(cè)井的誤差情況分析結(jié)果

3.2 出口斷面流量過(guò)程計(jì)算結(jié)果分析

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

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

圖7 模擬與實(shí)測(cè)逐日流量過(guò)程

圖8 2000年模擬與實(shí)測(cè)逐日流量過(guò)程

4 結(jié)論

以遼寧省大洋河岫巖流域?yàn)槔?，?duì)基于HYDRUSMODFLOW的地表水和地下水耦合模型進(jìn)行研究，獲得如下結(jié)論：

（1）研究區(qū)的地下水位的模型模擬結(jié)果和實(shí)測(cè)值具有基本一致的走向變化，兩者之間的相關(guān)系數(shù)為0.79，具有較好的擬合度。

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

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

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