張旭洋 林 青 黃修東 徐紹輝?

（1 青島大學(xué)環(huán)境科學(xué)與工程學(xué)院，山東青島 266071）

（2 青島市水文局，山東青島 266071）

大沽河水源地是青島市主要的供水水源地之一，由于大規(guī)模開采地下水，造成地下水位大幅下降、海水入侵等一系列環(huán)境問題，給當(dāng)?shù)厣鐣?huì)經(jīng)濟(jì)發(fā)展造成很大影響［1-2］。準(zhǔn)確計(jì)算大沽河流域的地下水補(bǔ)給量，對(duì)保障該地區(qū)經(jīng)濟(jì)可持續(xù)性發(fā)展和地下水資源的合理開發(fā)利用具有重要的理論和實(shí)際意義。土壤水是大氣降水、地表水、土壤水和地下水“四水轉(zhuǎn)化”中的一個(gè)非常重要的環(huán)節(jié)［3-5］，與地下水之間存在著密切聯(lián)系，且在一定條件下可以相互轉(zhuǎn)化，綜合考慮水分在包氣帶和地下含水層的運(yùn)動(dòng)過程，可以更加準(zhǔn)確地估算地下水補(bǔ)給量。

估算地下水補(bǔ)給的方法大致可以分為物理法、化學(xué)法和數(shù)值模擬方法。物理方法主要包括水均衡法［6-7］、零通量面法、地中滲透儀法［8］、地下水位波動(dòng)法［9］等；化學(xué)方法主要有人工示蹤劑法［10］和氯離子質(zhì)量平衡法［11-12］；數(shù)值模擬方法可以分為：基于土壤水動(dòng)力學(xué)過程的滲流帶模型［13］（Richards 方程）和地下水流模型［14］。由于研究對(duì)象和原理的限制，傳統(tǒng)方法大多只單獨(dú)從地下水或者包氣帶的角度考慮問題，忽略了土壤水和地下水之間的相互聯(lián)系，導(dǎo)致計(jì)算結(jié)果存在一定誤差。

將土壤水與地下水作為一個(gè)整體來進(jìn)行研究，可以更加全面地認(rèn)識(shí)土壤水、地下水的運(yùn)動(dòng)規(guī)律。韓雙平等［15］通過調(diào)節(jié)試驗(yàn)區(qū)的潛水埋深，研究了在不同水位埋深下，冬小麥和夏玉米生長(zhǎng)期內(nèi)土壤水和地下水的轉(zhuǎn)化關(guān)系；Adomako等［16］通過測(cè)量和分析降水、土壤水和地下水中δD和δ18O的含量差異，表明降水、土壤水和地下水之間關(guān)系密切；牛赟等［17］選取黑河荒漠過渡帶作為實(shí)驗(yàn)地點(diǎn)，分析了該地區(qū)降水、土壤水和地下水的相關(guān)性，并建立三者之間的回歸模型；竇超銀等［18］通過觀測(cè)土壤水分、鹽分以及地下水位在一年內(nèi)的變化情況，對(duì)地下水淺埋條件下土壤水分的運(yùn)動(dòng)進(jìn)行了研究；盧小慧等［19］以中國(guó)科學(xué)院欒城農(nóng)業(yè)生態(tài)試驗(yàn)站的地下水位觀測(cè)資料及氣象資料為基礎(chǔ)，綜合運(yùn)用降水、蒸發(fā)、土壤水、地下水動(dòng)態(tài)觀測(cè)資料，利用EARTH 模型計(jì)算了河北平原地下水垂向入滲補(bǔ)給量。Chen和Hu［20］在考慮水分在包氣帶和地下水運(yùn)動(dòng)的基礎(chǔ)上，開發(fā)出一款土壤水文模型，并利用Nebraska Sand Hills的實(shí)測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證，研究表明受到地下水位影響的土壤含水量可以提高21%，蒸發(fā)量也相應(yīng)提高7%～21%。

已有研究未從水動(dòng)力學(xué)角度將土壤水和地下水真正作為一個(gè)整體看待，也未建立流域尺度上的土壤水-地下水耦合模型，并以此計(jì)算大氣降水和灌溉水等對(duì)地下水的補(bǔ)給。本研究參照Seo等［21］的HYDRUS package for MODFLOW軟件的原理，結(jié)合GIS技術(shù)，建立流域尺度土壤水-地下水流耦合模型，利用研究區(qū)夏玉米生長(zhǎng)期內(nèi)的降水、蒸發(fā)、地下水位和土壤水分含量等數(shù)據(jù)對(duì)模型進(jìn)行校正，并用2012—2013年的資料進(jìn)行驗(yàn)證。

1 材料與方法

1.1 大沽河流域土壤水-地下水耦合模型的原理

土壤水-地下水耦合模型包含兩大部分：HYDRUS子模塊和MODFLOW主模塊。HYDRUS子模塊基于?im?nek等［22］編譯的HYDRUS-1D源程序，考慮了土壤水分運(yùn)動(dòng)的主要過程和影響因素，如降水、蒸發(fā)、滲透、毛細(xì)上升、植物根系吸水、地表積水和徑流、土壤水儲(chǔ)量等。MODFLOW模塊基于Harbaugh等［23］編譯的MODFLOW-2000源程序，主要用來模擬各種條件下水流在地下含水層中的運(yùn)動(dòng)。

HYDRUS子模塊和MODFLOW模塊耦合的原理主要基于Seo等［21］的HYDRUS package for MODFLOW軟件，將對(duì)土壤水分運(yùn)動(dòng)的模擬作為一個(gè)子程序耦合到地下水運(yùn)動(dòng)模型MODFLOW模塊中。地下水流區(qū)域的空間離散參照Harbaugh等［23］的離散方法，將含水層離散成若干網(wǎng)格，將地下水運(yùn)動(dòng)簡(jiǎn)化為平面二維流動(dòng)；土壤包氣帶根據(jù)植被類型、地下水埋深、土壤類型以及降水分布等因素的不同，將其劃分為N個(gè)不同的區(qū)域（N≤地下水網(wǎng)格單元），每個(gè)區(qū)域包含一個(gè)或多個(gè)地下水網(wǎng)格單元，該分區(qū)內(nèi)所有地下水網(wǎng)格單元的地下水位平均值將作為土壤底部邊界的水頭值，土壤包氣帶中水分的運(yùn)動(dòng)則簡(jiǎn)化為垂向一維流動(dòng)，這樣整個(gè)研究區(qū)從地表至隔水底板就被劃分為N個(gè)柱體，其中的水流運(yùn)動(dòng)被簡(jiǎn)化為準(zhǔn)三維流動(dòng)。每個(gè)土壤柱以地下水位為界，地下水位以下，用MODFLOW模塊求解；地下水位以上，用HYDRUS-1D模塊求解，在MODFLOW的每一個(gè)時(shí)間步長(zhǎng)內(nèi)，HYDRUS子模塊通過多次迭代求解Richards方程得到土壤剖面底部流量，MODFLOW模塊獲得該流量經(jīng)計(jì)算得到一個(gè)新的水位，作為下一計(jì)算步長(zhǎng)內(nèi)土壤剖面底部邊界的水頭值。HYDRUS模塊和MODFLOW模塊耦合計(jì)算的過程如圖1所示。

HYDRUS-MODFLOW耦合模型中，水在土壤和地下含水層中的流動(dòng)分別采用不同的方程，土壤水分運(yùn)動(dòng)采用修正的Richards方程，忽略土壤水側(cè)向流動(dòng)，僅考慮一維垂向運(yùn)移和根系吸水項(xiàng)，方程求解采用伽遼金線性有限元法。方程具體形式如下：

式中，C(h)為比水容量（cm-1），C(h)=dθ/dh，θ為土壤體積含水量（cm3·cm-3）；h為壓力水頭（mm）；K(h)為非飽和土壤導(dǎo)水率（cm·d-1）；S(z, t)為t時(shí)刻z深度處耗水速率，取該處作物根系吸水率（cm3·d-1）；t為時(shí)間（d）；z為土壤深度（cm），坐標(biāo)向下為正。

圖1 土壤水-地下水耦合模型流程圖Fig. 1 Schematic diagram of soil water/groundwater coupled model

地下水流動(dòng)數(shù)值模擬的理論基礎(chǔ)是孔隙介質(zhì)中地下水二維運(yùn)動(dòng)方程、定解條件及數(shù)值計(jì)算方法。運(yùn)動(dòng)方程根據(jù)質(zhì)量守恒定律導(dǎo)出，定解條件由研究區(qū)域地下水的初始條件與計(jì)算區(qū)域的邊界條件給定，地下潛水非穩(wěn)定流的二維運(yùn)動(dòng)偏微分方程為：

式中，K為滲透系數(shù)(m·d-1)；h為潛水層的厚度（m）；H為含水層任一點(diǎn)的水頭標(biāo)高（m）；W為單位時(shí)間單位面積上的垂直水量交換（m·d-1）；μ為給水度；t為時(shí)間（d）；q(x, y, z)為側(cè)向補(bǔ)給量；φ為Γ2-1上的已知函數(shù)；Γ2-1為第一類邊界，定水頭邊界；Γ2-2為第二類邊界，隔水邊界；Γ2-3為二類邊界，側(cè)向補(bǔ)給邊界。

1.2 研究區(qū)概況

大沽河是山東半島主要河流之一，位于東經(jīng)120°03′—120°25′，北緯36°10′—37°12′之間。青島市境內(nèi)流域總面積為4 781 km2（見圖2），屬華北暖溫帶沿海濕潤(rùn)季風(fēng)氣候，多年平均氣溫12.3℃，多年平均降水量685.4 mm，降水主要集中在6—9月份，約占全年降水量的69.2%。多年平均蒸發(fā)量為983.9 mm，是平均降雨量的1.45倍。流域內(nèi)農(nóng)作物主要為冬小麥—夏玉米輪作。

1.3 大沽河流域水文地質(zhì)概念模型

圖2 大沽河流域的地理位置Fig. 2 Geographical location of the Dagu River Basin

土壤水-地下水流耦合模型的上邊界條件為大氣邊界，以大氣降水量、農(nóng)田灌溉水量、地表蒸發(fā)量和植物蒸騰量作為模型的主要源匯項(xiàng)；耦合模型的下邊界條件設(shè)定為隔水邊界，這主要是由于大沽河流域地下含水層下伏地層為滲透能力極差的基巖；將研究區(qū)的東西兩側(cè)視為隔水邊界，主要原因是研究區(qū)內(nèi)的含水砂層沿大沽河周圍分布、向東西兩側(cè)尖滅，且東西兩邊也是流域地表分水嶺邊界；研究區(qū)北部的含水層在大、小沽河入境處與外界的含水層相互連通，故將研究區(qū)的北部邊界視為透水邊界；南部邊界根據(jù)課題組多年收集的地下水位實(shí)際觀測(cè)資料，取多年平均水位作為定水頭邊界。

大沽河流域地下水的補(bǔ)給來源主要包括大氣降水和農(nóng)田灌溉；含水層的排泄項(xiàng)主要是潛水蒸發(fā)及人工開采。

1.4 模型離散化

研究區(qū)總面積為4 781 km2，將滲流區(qū)剖分為113行×77列，每個(gè)單元格為1 000 m×1 000 m，共劃分為8 701個(gè)正方形單元，其中有效單元4 750個(gè)，非活動(dòng)單元3 951 個(gè)。根據(jù)以往對(duì)大沽河水源地，即地下水庫(kù)的勘察結(jié)果［24-25］，將整個(gè)流域的含水層參數(shù)分為4個(gè)區(qū)，其中地下水庫(kù)區(qū)域由于含水層給水度大且差異明顯，細(xì)分為3個(gè)區(qū)；地下水庫(kù)以外的區(qū)域，因透水性差，故概化為一個(gè)區(qū)。給水度和滲透系數(shù)分區(qū)一致，如圖3a所示。含水層垂向上簡(jiǎn)化為一層，為潛水含水層。

流域內(nèi)土壤類型根據(jù)本課題組多年采樣的測(cè)定結(jié)果，并通過kriging插值獲得黏粒和砂粒在整個(gè)流域內(nèi)不同深度的分布情況，得到相應(yīng)的土壤質(zhì)地；流域內(nèi)的土地利用類型依據(jù)歐盟聯(lián)合研究中心(Joint Research Centre，JRC)空間應(yīng)用研究所(Institute of Space Applications，SAI)2000年全球土地覆蓋數(shù)據(jù)（2000 Global Land Cover Data，GLC2000），將流域土地利用類型簡(jiǎn)化為耕地（3 461 km2）和城鎮(zhèn)居民用地（1 289 km2）；由于大沽河流域地下水埋深大都小于6 m，將地下水埋深劃分為4個(gè)等級(jí)：0～2 m、2～4 m、4～6 m和＞6 m；鑒于流域內(nèi)降雨量空間分布的不一致性，考慮到研究區(qū)內(nèi)有藍(lán)村、嵐西頭、郭莊等8個(gè)氣象站，采用泰森多邊形方法進(jìn)行分區(qū)，每一分區(qū)內(nèi)的降雨條件視為均一，分區(qū)內(nèi)的降雨由所屬氣象站數(shù)據(jù)代表。

針對(duì)流域內(nèi)土壤類型、土地利用類型、地下水埋深以及降水分布等具體情況，將整個(gè)流域劃分為47個(gè)區(qū)域如圖3b所示，每個(gè)區(qū)域包含一個(gè)或多個(gè)網(wǎng)格單元，土壤包氣帶中的水分運(yùn)動(dòng)簡(jiǎn)化為一維流動(dòng)。

圖3 模型空間離散示意圖Fig. 3 Schematic diagrams of spatial discreteness

根據(jù)模擬期間地下水埋深的變化情況，為使不同分區(qū)土壤剖面的底部邊界在模擬期間始終處于地下水位以下，以含水層的隔水底板作為土壤剖面的底部邊界。

1.5 模型相關(guān)參數(shù)及初始條件

以地下水的埋深為界，將土壤剖面簡(jiǎn)化為2層：下層土壤質(zhì)地為砂土，代表地下水位波動(dòng)帶的土壤質(zhì)地；上層根據(jù)實(shí)測(cè)數(shù)據(jù)細(xì)分為2～3種不同的質(zhì)地，土壤水力學(xué)參數(shù)以課題組多年實(shí)測(cè)數(shù)據(jù)為準(zhǔn)。

表1 不同類型土壤水力學(xué)參數(shù)Table 1 Soil hydraulic parameters relative to type of the soil

研究區(qū)地下水開采量參照青島市水利局發(fā)布的2012、2013年《青島市水資源公報(bào)》的統(tǒng)計(jì)結(jié)果。青島地區(qū)的地下水開采主要集中在地下水庫(kù)范圍內(nèi)，地下水庫(kù)以外由于沒有好的含水層，其滲透系數(shù)和給水度較小，開采量較低。因此，根據(jù)不同行政區(qū)淺層淡水的總開采量求得不同行政區(qū)在地下水庫(kù)內(nèi)、外每平方千米的開采量，將其作為耦合模型中每個(gè)單元格的地下水開采量，以抽水井的形式加入到模型中，流域內(nèi)每個(gè)單元井的開采量如表2所示。

表2 大沽河流域年均供水量統(tǒng)計(jì)結(jié)果Table 2 Statistics of annual water supply of the Dagu River Basin

根系吸水模型采用van Genuchten宏觀吸水模型，根系吸水參數(shù)選取HYDRUS-1D數(shù)據(jù)庫(kù)中的玉米（wesseling，1991）和小麥(wesseling，1991)經(jīng)驗(yàn)值［22］，由于夏玉米生長(zhǎng)期內(nèi)，大約91.6 %的根系分布在0～1 m的土壤層內(nèi)，故模型中僅在0～1 m的土層設(shè)置根系吸水項(xiàng)。降水、灌溉水和田間蒸散發(fā)構(gòu)成作物生長(zhǎng)期內(nèi)土壤水分運(yùn)動(dòng)模型的上邊界條件和源匯項(xiàng)，農(nóng)業(yè)灌溉量根據(jù)當(dāng)?shù)剞r(nóng)民的實(shí)際灌溉時(shí)間和強(qiáng)度以降雨的方式加入到模型中，在夏玉米整個(gè)生長(zhǎng)期內(nèi)，灌溉主要集中在播種時(shí)；冬小麥生育期內(nèi)灌溉一次，主要集中在拔節(jié)期前后（生育期第140天左右），每次的灌水量80 mm［26-27］。降水和蒸發(fā)數(shù)據(jù)主要來自于青島市水利局2012年和2013年的氣象統(tǒng)計(jì)結(jié)果。

大沽河流域土壤水-地下水耦合模型的初始流場(chǎng)，以模擬時(shí)間開始時(shí)的地下水位的監(jiān)測(cè)資料作為依據(jù)，在流域的上、中、下游一共選取35口水位監(jiān)測(cè)井的實(shí)測(cè)數(shù)據(jù)，采用內(nèi)插法形成耦合模型的初始地下水流場(chǎng)。

由于流域內(nèi)地勢(shì)平坦，地表導(dǎo)水率較大，故模擬中不考慮地表徑流。模型模擬期內(nèi)，大沽河流域降水量偏小，河流大部分時(shí)間處于干涸狀態(tài)，故不考慮河流入滲對(duì)地下水的補(bǔ)給。

2 結(jié)果與討論

2.1 校正期土壤剖面含水量和地下水位的擬合結(jié)果

利用土壤水-地下水流耦合模型和上述參數(shù)對(duì)2013年6月16日—2013年9月29日夏玉米生長(zhǎng)期內(nèi)（共105 d）土壤水和地下水的運(yùn)動(dòng)過程進(jìn)行模擬，將土壤剖面含水量、地下水位的模擬結(jié)果分別與實(shí)測(cè)數(shù)據(jù)進(jìn)行擬合分析，采用均方根誤差（RMSE）和決定系數(shù)（R2）來評(píng)價(jià)耦合模型的模擬效果。

根據(jù)已有的監(jiān)測(cè)資料，隨機(jī)選取大沽河流域上中下游6個(gè)土壤水分監(jiān)測(cè)點(diǎn)（5#、3#、2#、8#、1#、10#）在7月12日的實(shí)測(cè)數(shù)據(jù)與耦合模擬所得的土壤剖面含水量進(jìn)行擬合分析，除個(gè)別點(diǎn)外，土壤剖面含水量的R2在0.65～0.91之間，RMSE在0.005～0.01之間，因篇幅所限只給出5#、2#和10#三個(gè)監(jiān)測(cè)點(diǎn)的擬合結(jié)果，如圖4所示。由于模擬時(shí)間處于夏玉米生長(zhǎng)旺期，玉米蒸散發(fā)量較大，而且7月9日至7月12日，大沽河流域經(jīng)歷了較強(qiáng)的降水過程，導(dǎo)致整個(gè)土壤剖面含水量的變化比較劇烈，模擬結(jié)果較差。5#和10#觀測(cè)點(diǎn)0.4～0.6 m 處的實(shí)測(cè)土壤含水量明顯偏高，這是因?yàn)樵搶油寥勒沉：枯^高，持水性好，造成水分在該層集聚，從而對(duì)水分的下滲補(bǔ)給過程產(chǎn)生重要影響。土壤剖面含水量的擬合結(jié)果表明：耦合模型模擬出的土壤剖面含水量與研究區(qū)實(shí)測(cè)點(diǎn)的實(shí)測(cè)數(shù)據(jù)的變化趨勢(shì)基本相同，模擬得到的土壤含水量與實(shí)測(cè)值總體上差別不大，說明耦合模型能較準(zhǔn)確地模擬玉米生長(zhǎng)期內(nèi)水分在包氣帶中的運(yùn)動(dòng)過程。

圖4 校正期土壤剖面含水量的擬合結(jié)果Fig. 4 Fitting of soil profile water content in the calibration phase

將耦合模擬所得的地下水位和實(shí)測(cè)水位進(jìn)行擬合分析，隨機(jī)選取大沽河上中下游有地下水位觀測(cè)井的三個(gè)網(wǎng)格，繪制地下水位隨時(shí)間的變化趨勢(shì)線，如圖5所示，地下水位的R2在0.52～0.68之間，RMSE在0.12～0.38之間。從圖5中可以看出，實(shí)測(cè)地下水位在7月16日前后有明顯的升高，根據(jù)青島市水利局統(tǒng)計(jì)的2013年降水?dāng)?shù)據(jù)顯示，在7月9日至7月12日，青島市有明顯的降雨過程，降雨量達(dá)到100 mm，由于大沽河流域平均地下水埋深在4～6 m左右，而且可能存在大孔隙優(yōu)先流，較強(qiáng)程度的降水可以迅速入滲補(bǔ)給到地下水，在短時(shí)間內(nèi)導(dǎo)致地下水位明顯上升。而模型中沒有考慮優(yōu)先流，對(duì)于地下水位升高的現(xiàn)象，模擬結(jié)果滯后于實(shí)際觀測(cè)結(jié)果，并且模擬曲線比實(shí)測(cè)曲線更加平緩。地下水位的擬合結(jié)果表明：土壤水-地下水流耦合模型模擬得出的地下水位與研究區(qū)實(shí)測(cè)地下水位的變化趨勢(shì)基本相同，耦合模擬得到的地下水位與實(shí)測(cè)地下水位總體上差別不大，說明耦合模型能較好地模擬玉米生長(zhǎng)期內(nèi)水分在地下含水層中的運(yùn)動(dòng)過程。

圖5 校正期地下水位的擬合結(jié)果Fig.5 Fitting of groundwater table in the calibration phase

圖6為2013年8月10日研究區(qū)范圍內(nèi)地下水位實(shí)測(cè)值和模擬值的等值線圖。從圖中可以看出，在整個(gè)研究區(qū)范圍內(nèi)，模擬地下水位和實(shí)測(cè)地下水位大體相同，證明耦合模型可以準(zhǔn)確地模擬流域范圍內(nèi)地下水位的變化趨勢(shì)。

圖6 地下水位模擬值和實(shí)測(cè)值等值線圖Fig. 6 Contour map of simulated and measured groundwater tables

基于對(duì)土壤剖面含水量和地下水位擬合結(jié)果的分析，可以證明土壤水-地下水流耦合模型的概化及選取的土壤水力學(xué)參數(shù)和水文地質(zhì)參數(shù)是合理的，校正后的土壤水-地下水流耦合模型能夠比較準(zhǔn)確地模擬大沽河流域夏玉米生長(zhǎng)期內(nèi)水分在包氣帶和地下含水層中的運(yùn)動(dòng)情況，因此耦合計(jì)算所得的地下水補(bǔ)給量具有較高的準(zhǔn)確性。識(shí)別后的含水層給水度μ和滲透系數(shù)K如表3所示。

2.2 驗(yàn)證期土壤剖面含水量和地下水位的擬合結(jié)果

利用校正后的耦合模型對(duì)大沽河流域土壤水和地下水的運(yùn)動(dòng)過程進(jìn)行模擬。模擬時(shí)間從2012年6月16日至2013年6月16日，模擬時(shí)間根據(jù)夏玉米和冬小麥的生長(zhǎng)過程分為兩個(gè)應(yīng)力期，夏玉米期從2012年6月16日至2012年9月29日；冬小麥期從2012年9月30日至2013年6月16日。因篇幅限制，隨機(jī)選取大沽河流域上中下游（K094029A、K094007A和K0920430）三個(gè)水位觀測(cè)點(diǎn)的實(shí)測(cè)數(shù)據(jù)說明地下水位的擬合情況；隨機(jī)選取上中下游（5#、2#和1#）三個(gè)土壤水分監(jiān)測(cè)點(diǎn)在2012年7月14日和2013年3月28日的實(shí)測(cè)數(shù)據(jù)說明土壤剖面含水量的擬合結(jié)果。

表3 不同分區(qū)給水度μ和滲透系數(shù)K值表Table 3 Specific yield ( μ ) and permeation coefficient (K ) relative to zone

由圖7和圖8可知：驗(yàn)證期耦合模擬的地下水位、土壤剖面含水量與實(shí)測(cè)數(shù)據(jù)的擬合度較高：地下水位的模擬值與實(shí)測(cè)值的R2在0.56～0.81，RMSE為0.17～0.19，擬合結(jié)果較好。2012年7月14日土壤剖面含水量的R2為0.53～0.86，RMSE為0.006～0.011；2013年3月28日土壤剖面含水量的R2在0.73～0.92，RMSE為0.008～0.011。冬小麥生長(zhǎng)期內(nèi)土壤剖面含水量的擬合結(jié)果明顯比夏玉米生長(zhǎng)期的擬合結(jié)果好，這是因?yàn)樵谙挠衩咨L(zhǎng)期，降水和植物蒸散發(fā)作用強(qiáng)烈，土壤剖面含水量變化劇烈，導(dǎo)致擬合結(jié)果較差。驗(yàn)證期內(nèi)地下水位和土壤剖面含水量的擬合結(jié)果表明，土壤水-地下水耦合模型能夠準(zhǔn)確的模擬研究區(qū)土壤水和地下水的運(yùn)動(dòng)過程，模擬所得的地下水補(bǔ)給量準(zhǔn)確度和可信度較高，可以用于計(jì)算不同降雨水平年條件下地下水的補(bǔ)給量。

根據(jù)耦合模型的輸出文件，得到含水層的水量平衡關(guān)系，進(jìn)而求取地下水的補(bǔ)給量，地下含水層的水量平衡關(guān)系如表4所示。校正期大沽河流域大氣降水量為328 mm，灌溉量為80 mm，耦合模擬計(jì)算出的地下水補(bǔ)給量為3.15×109m3；驗(yàn)證期大沽河流域大氣降水量為653 mm，灌溉量為160 mm，土壤水-地下水耦合模型計(jì)算出的地下水補(bǔ)給量為4.77×109m3，可以為制定流域水資源優(yōu)化配置方案提供科學(xué)依據(jù)。

3 結(jié) 論

基于HYDRUS package for MODFLOW軟件的原理，結(jié)合GIS技術(shù)，建立大沽河流域土壤水-地下水流耦合模型。根據(jù)2013年夏玉米生長(zhǎng)期內(nèi)的降水、蒸發(fā)、地下水位以及土壤剖面含水量等資料，構(gòu)建2013年夏玉米生長(zhǎng)期內(nèi)土壤水-地下水耦合運(yùn)動(dòng)模型，經(jīng)過模型校正，地下水位、土壤剖面含水量的模擬值和實(shí)測(cè)值吻合度較好，在夏玉米生長(zhǎng)期內(nèi)，整個(gè)流域范圍內(nèi)地下水補(bǔ)給量為3.15×109m3。校正期的擬合結(jié)果表明模型的概化及選取的土壤水力學(xué)參數(shù)和水文地質(zhì)參數(shù)是比較合理的，校正后的耦合模型基本可以反映大沽河流域土壤水和地下水的運(yùn)動(dòng)狀況，耦合計(jì)算所得的地下水補(bǔ)給量可信度較高。利用校正后的模型及參數(shù)，對(duì)2012年6月16日—2013年6月16日期間的土壤水和地下水運(yùn)動(dòng)情況進(jìn)行模擬，結(jié)果表明：土壤水和地下水流耦合模型能夠較好的反映研究區(qū)土壤水和地下水的時(shí)空變化，耦合計(jì)算所得的地下水補(bǔ)給量為4.77×109m3。土壤水-地下水流耦合模型能夠較好地模擬大沽河流域土壤水和地下水的時(shí)空變化，并且能夠準(zhǔn)確地計(jì)算地下水垂向入滲補(bǔ)給量，可以為流域尺度制定更加合理的水資源優(yōu)化配置方案提供科學(xué)依據(jù)。

圖7 驗(yàn)證期地下水位的擬合結(jié)果Fig. 7 Fitting of groundwater table in the validation phase

圖8 驗(yàn)證期土壤剖面含水量的擬合結(jié)果Fig. 8 Fitting of g soil profile water content in the validation phase

表4 驗(yàn)證期和校正期地下水水量平衡表Table 4 Groundwater balance in the validation and calibration phases /(×106 m3)