基于地表水與地下水耦合模型模擬計(jì)算

王祥偉

（東港市防汛抗旱指揮部辦公室，遼寧 東港118300）

1 研究區(qū)概況

青龍河是灤河重要支流，其干流總長(zhǎng)度246km，發(fā)源于河北省平泉縣臺(tái)頭山鄉(xiāng)，流經(jīng)遼寧省凌源市，向南流入河北省青龍縣，最終在灤縣匯入灤河入渤海［1］。青龍河在凌源市境內(nèi)流經(jīng)劉杖子鄉(xiāng)、三道河子鄉(xiāng)、刀爾登本鄉(xiāng)及河坎子鄉(xiāng)等4個(gè)鄉(xiāng)鎮(zhèn)，流域面積589km2。研究區(qū)處于遼西低山丘陵地帶，受到燕山構(gòu)造運(yùn)動(dòng)影響，地質(zhì)環(huán)境比較復(fù)雜。凌源市青龍河流域?qū)儆谶|西半干旱地區(qū)，該地區(qū)屬于中溫帶內(nèi)陸地區(qū)，屬于溫帶大陸性季風(fēng)氣候。流域內(nèi)年平均氣溫較低，降水量小但是蒸發(fā)量較大［2］。流域內(nèi)具有典型的大陸性氣候特征，春秋降水較少風(fēng)力較大，夏季短促多短時(shí)強(qiáng)降雨、秋季干爽降雨較少，冬季漫長(zhǎng)寒冷干燥。流域內(nèi)多年平均溫度為5.6 ℃，極端最低氣溫-42.4 ℃，極端最高氣溫39.2 ℃。降雨主要集中于夏季的7，8月份，約占全年降水量的60%～75%，流域最大年降雨量683mm，最小年降雨量221mm。流域內(nèi)的潛水主要是松散巖孔隙水和基巖裂隙水，承壓水主要為碎屑巖裂隙水與基巖裂隙水，在天然狀態(tài)下的補(bǔ)給主要依靠大氣降水和基巖裂隙水的側(cè)向補(bǔ)給，并主要表現(xiàn)為丘陵區(qū)向山谷流動(dòng)。

2 地表水與地下水耦合模型的構(gòu)建

2.1 模型構(gòu)建思路

水循環(huán)是一個(gè)比較復(fù)雜的運(yùn)動(dòng)過程，地表水和地下水等不同類型水資源的賦存條件及運(yùn)動(dòng)狀態(tài)也存在諸多差異?；谒h(huán)的復(fù)雜性，當(dāng)前對(duì)流域水循環(huán)的研究主要表現(xiàn)為不同類型水資源的單獨(dú)模擬，而這種研究方式容易導(dǎo)致水循環(huán)過程缺失，進(jìn)而影響研究成果，因此應(yīng)該將地表水與地下水這兩種主要水資源形式進(jìn)行綜合研究。本文利用寧吉才等人的研究成果［3-4］，對(duì)改進(jìn)后的SWAT模型，結(jié)合研究區(qū)地表水文資料進(jìn)行參數(shù)率定，在校驗(yàn)準(zhǔn)確并滿足模型精度要求的基礎(chǔ)上（限于篇幅這里不再詳述），將其輸出項(xiàng)作為Visual MODFLOW的輸入項(xiàng)對(duì)進(jìn)行數(shù)值模擬，以凌源市青龍河流域地表水和地下水耦合模擬方法進(jìn)行研究，通過考慮地表水影響提高對(duì)地下水的模擬精度。

2.2 Visual MODFLOW軟件與耦合原理

Visual MODFLOW軟件由加拿大開發(fā)的地下水流與溶質(zhì)運(yùn)移模擬軟件，主要由MODFLOW、MODPATH、MT3D 3個(gè)地下水模塊組成。該軟件利用有限差分法實(shí)現(xiàn)對(duì)地下水流的數(shù)值模擬，同時(shí)可將不同模擬過程進(jìn)行串聯(lián)，使模型輸入、運(yùn)行和結(jié)果輸出過程比較便捷和清晰。

雖然利用Visual MODFLOW軟件對(duì)區(qū)域地下水運(yùn)動(dòng)規(guī)律進(jìn)行建模分析具有諸多優(yōu)勢(shì)，但該軟件弱點(diǎn)十分明顯，因?yàn)槟Ｐ托枰揽刻囟〝?shù)據(jù)的輸入進(jìn)行運(yùn)行和模擬［5］，當(dāng)模型通過參數(shù)輸入的方式進(jìn)行數(shù)值模擬時(shí)，參數(shù)估值造成模擬計(jì)算精度降低。要實(shí)現(xiàn)地表水與地下水模型的耦合，首要問題是解決兩模型在最小單元上的空間差異。改進(jìn)SWAT模型以HRU （水文響應(yīng)單元） 為最小計(jì)算單元，而Visual MODFLOW模型以CELL（有限差分網(wǎng)格）為最小計(jì)算單元。解決思路是在構(gòu)建起研究區(qū)的SWAT模型，并進(jìn)行參數(shù)率定之后，輸入研究期（2016年6月～2018年6月）氣象數(shù)據(jù)進(jìn)行徑流模擬，并基于模擬結(jié)果獲得土壤蒸發(fā)量和降雨入滲補(bǔ)給兩個(gè)主要參數(shù)。然后，利用Arcgis軟件將HRU轉(zhuǎn)化為ASCII格式的文件，在賦予其空間屬性后實(shí)現(xiàn)與Visual MODFLOW 模型的CELL的對(duì)應(yīng)關(guān)系，進(jìn)而實(shí)現(xiàn)兩模型的耦合［6-7］。

2.3 模型離散

在模型構(gòu)建過程中，利用Arcgis軟件從研究區(qū)的DEM資料中提取到地表高程的數(shù)據(jù)，在對(duì)部分異常值進(jìn)行剔除之后，對(duì)獲得數(shù)據(jù)利用Surfer軟件進(jìn)行插值處理獲得地表高程高柵格數(shù)據(jù)，然后，利用隔水基巖的具體位置，確定隔水底板，最終將研究區(qū)在垂直方向上劃分為一層。在對(duì)研究區(qū)進(jìn)行研究邊界導(dǎo)入后，將整個(gè)研究區(qū)劃分為縱橫各50行，共2500個(gè)網(wǎng)格單元，其中有效計(jì)算單元為1405個(gè)。此外，為了驗(yàn)證模型的擬合度，在研究區(qū)內(nèi)選擇6個(gè)地下水位觀測(cè)孔，作為模型校準(zhǔn)和驗(yàn)證依據(jù)。

2.4 模型參數(shù)與初始條件

水文地質(zhì)參數(shù)可以對(duì)地下水含水層的儲(chǔ)水、輸水特征進(jìn)行定量表征，是進(jìn)行模型建構(gòu)和計(jì)算的重要指標(biāo)［8］。鑒于本次研究的主要目的是通過模型計(jì)算對(duì)研究區(qū)內(nèi)的潛水含水層的地下水運(yùn)動(dòng)規(guī)律進(jìn)行模擬，因此，不同地層滲透系數(shù)與給水度數(shù)值是否準(zhǔn)確，將直接影響到模型表征地下水運(yùn)動(dòng)實(shí)際的契合度。另一方面，在耦合模型構(gòu)建中需要將SWAT模型的輸出數(shù)據(jù)輸入Visual MODFLOW，由于研究區(qū)的土

表1 水文地質(zhì)參數(shù)經(jīng)驗(yàn)值 單位：m/d

表2 研究區(qū)水文地質(zhì)參數(shù)初始值 單位：m/d

2.5 模型的匯源項(xiàng)處理

研究區(qū)主要地貌為丘陵低山，耕地不多且主要分布在河谷地帶，因此大氣降水入滲及側(cè)向流入是流域內(nèi)地下水補(bǔ)給的主要形式，潛水蒸發(fā)、側(cè)向排泄、河流排泄及人工開采則是地下水的主要排泄形式。在研究中，將利用SWAT模型模擬獲取的大氣降雨入滲補(bǔ)給與潛水蒸發(fā)進(jìn)行基于HRU和CELL的空間耦合，并將結(jié)果作為補(bǔ)給和蒸發(fā)邊界輸入到Visual MODFLOW中。由于青龍河凌源段屬于丘陵河谷區(qū)，河流水位一般低于周邊的地下水位，屬于淺層地下水排泄區(qū)。因此，基于河流和周邊地區(qū)地下水水位之間的變化關(guān)系，就能對(duì)排泄量進(jìn)行相應(yīng)計(jì)算。此外，研究區(qū)內(nèi)的地下水開采主要是居民生活用水，可根據(jù)統(tǒng)計(jì)資料對(duì)開采井進(jìn)行日概化分配。最終得到研究區(qū)地下水匯源項(xiàng)統(tǒng)計(jì)結(jié)果，如表3。壤類型單一，在參數(shù)分區(qū)時(shí)可將青龍河凌源段的每個(gè)子流域作為一個(gè)水文地質(zhì)區(qū)，這樣共獲得7個(gè)水文地質(zhì)區(qū)。在參數(shù)初始值選取中，結(jié)合表1，相關(guān)領(lǐng)域的研究經(jīng)驗(yàn)值，在通過抽水試驗(yàn)獲取如表2所示的水文地質(zhì)參數(shù)初始值。

表3 研究區(qū)地下水匯源項(xiàng)統(tǒng)計(jì)結(jié)果 單位：106m3/a

3 模型識(shí)別和驗(yàn)證

3.1 模型的識(shí)別

在模型構(gòu)建完成之后，選擇2016年6月1日的研究區(qū)地下水流場(chǎng)作為初始流場(chǎng)，將2017年5月30日作為末端時(shí)刻，將一個(gè)完整的水文年作為模擬期。在模型運(yùn)行結(jié)束之后，對(duì)6個(gè)選定的地下水位觀測(cè)孔的實(shí)測(cè)值進(jìn)行對(duì)比，結(jié)果顯示，2#觀測(cè)井的水位殘差最小，為0.88m，而4#觀測(cè)井的殘差最大，為30.16m，所有6個(gè)觀測(cè)孔的均值為13.27m。因此，模型的模擬結(jié)果與實(shí)測(cè)值之間存在較大差別，需要對(duì)水文地質(zhì)參數(shù)進(jìn)行不優(yōu)化調(diào)整。

利用Visual Modflow 軟件自帶的校準(zhǔn)程序進(jìn)行自動(dòng)擬合，并反復(fù)調(diào)整水文地質(zhì)參數(shù)，直至獲得最佳擬合狀態(tài)，概化和優(yōu)化后水文地質(zhì)參數(shù)如表4。

表4 水文地質(zhì)參數(shù)識(shí)別結(jié)果 單位：m/d

3.2 模型驗(yàn)證

研究中選取2017年6月1日至2018年5月30日這一水文年作為驗(yàn)證期，將該水文年的地下水初始流場(chǎng)及匯源項(xiàng)等相關(guān)數(shù)據(jù)輸入和運(yùn)行模型，并與6個(gè)觀測(cè)井的實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，進(jìn)而對(duì)模型的模擬結(jié)果的精確度進(jìn)行驗(yàn)證。驗(yàn)證期間，6個(gè)觀測(cè)井的實(shí)測(cè)與模擬計(jì)算結(jié)果的擬合效果如圖1～圖6。

圖1 A觀測(cè)井?dāng)M合結(jié)果

圖2 B觀測(cè)井?dāng)M合結(jié)果

圖3 C觀測(cè)井?dāng)M合結(jié)果

圖4 D觀測(cè)井?dāng)M合結(jié)果

圖5 E觀測(cè)井?dāng)M合結(jié)果

圖6 F觀測(cè)井?dāng)M合結(jié)果

由圖1～圖6擬合曲線可知，6個(gè)水位觀測(cè)井實(shí)測(cè)值和模擬值之間具有基本一致的變化趨勢(shì)，且兩者之間差別不大，擬合度良好。因此，模型符合凌源市青龍河流域的水文地質(zhì)條件，能夠用于研究區(qū)地下水運(yùn)動(dòng)變化特征的相關(guān)研究。

4 結(jié)語

（1）利用改進(jìn)后的SWAT模型，將其輸出項(xiàng)作為Visual Modflow的輸入項(xiàng)對(duì)凌源市青龍河流域的地下水進(jìn)行數(shù)值模擬，建立起地表水和地下水的耦合模型。

（2） 通過模擬期模型運(yùn)行結(jié)果，利用Visual Modflow 軟件自帶的校準(zhǔn)程序進(jìn)行自動(dòng)擬合，反復(fù)調(diào)整水文地質(zhì)參數(shù)，獲得符合流域地下水流運(yùn)動(dòng)特征的水文地質(zhì)參數(shù)。

（3）模型驗(yàn)證結(jié)果顯示，模擬值和實(shí)測(cè)值之間存在良好的擬合度，耦合模型可用于研究區(qū)地下水運(yùn)動(dòng)變化特征的相關(guān)研究。