大沽河流域土壤水-地下水流耦合模擬及補(bǔ)給量估算*

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

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

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

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

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

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

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

1 材料與方法

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

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

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

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

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

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

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

式中，K為滲透系數(shù)(m·d-1)；h為潛水層的厚度（m）；H為含水層任一點的水頭標(biāo)高（m）；W為單位時間單位面積上的垂直水量交換（m·d-1）；μ為給水度；t為時間（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），屬華北暖溫帶沿海濕潤季風(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ā)量和植物蒸騰量作為模型的主要源匯項；耦合模型的下邊界條件設(shè)定為隔水邊界，這主要是由于大沽河流域地下含水層下伏地層為滲透能力極差的基巖；將研究區(qū)的東西兩側(cè)視為隔水邊界，主要原因是研究區(qū)內(nèi)的含水砂層沿大沽河周圍分布、向東西兩側(cè)尖滅，且東西兩邊也是流域地表分水嶺邊界；研究區(qū)北部的含水層在大、小沽河入境處與外界的含水層相互連通，故將研究區(qū)的北部邊界視為透水邊界；南部邊界根據(jù)課題組多年收集的地下水位實際觀測資料，取多年平均水位作為定水頭邊界。

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

1.4 模型離散化

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

流域內(nèi)土壤類型根據(jù)本課題組多年采樣的測定結(jié)果，并通過kriging插值獲得黏粒和砂粒在整個流域內(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），將流域土地利用類型簡化為耕地（3 461 km2）和城鎮(zhèn)居民用地（1 289 km2）；由于大沽河流域地下水埋深大都小于6 m，將地下水埋深劃分為4個等級：0～2 m、2～4 m、4～6 m和＞6 m；鑒于流域內(nèi)降雨量空間分布的不一致性，考慮到研究區(qū)內(nèi)有藍(lán)村、嵐西頭、郭莊等8個氣象站，采用泰森多邊形方法進(jìn)行分區(qū)，每一分區(qū)內(nèi)的降雨條件視為均一，分區(qū)內(nèi)的降雨由所屬氣象站數(shù)據(jù)代表。

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

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

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

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

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

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

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

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

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

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

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

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

由圖7和圖8可知：驗證期耦合模擬的地下水位、土壤剖面含水量與實測數(shù)據(jù)的擬合度較高：地下水位的模擬值與實測值的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。冬小麥生長期內(nèi)土壤剖面含水量的擬合結(jié)果明顯比夏玉米生長期的擬合結(jié)果好，這是因為在夏玉米生長期，降水和植物蒸散發(fā)作用強(qiáng)烈，土壤剖面含水量變化劇烈，導(dǎo)致擬合結(jié)果較差。驗證期內(nèi)地下水位和土壤剖面含水量的擬合結(jié)果表明，土壤水-地下水耦合模型能夠準(zhǔn)確的模擬研究區(qū)土壤水和地下水的運動過程，模擬所得的地下水補(bǔ)給量準(zhǔn)確度和可信度較高，可以用于計算不同降雨水平年條件下地下水的補(bǔ)給量。

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

3 結(jié) 論

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

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

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

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