楊妮娟, 王曉云, 李建柱, 張 婷, 馮 平

(1.天津大學(xué) 水利工程仿真與安全國家重點實驗室, 天津 300350; 2.中水北方勘測設(shè)計研究有限責(zé)任公司, 天津 300222)

1 研究背景

非點源污染是目前我國水環(huán)境治理的主要對象[1-3],由于其過程具有隨機(jī)性和分散性,污染的形成與流失機(jī)理較為復(fù)雜,難以對其進(jìn)行有效的管理與控制[4-5]。非點源污染受眾多因素影響,其中土地利用方式直接影響著流域的水文循環(huán)與物質(zhì)遷移過程,成為影響非點源污染的關(guān)鍵因子[6-7],由于近年來極端氣候事件頻發(fā)和人類活動加劇,土地利用方式的時空變化顯著,因此明晰不同土地利用方式對污染流失的影響以及流域環(huán)境對土地利用時空變化的響應(yīng)具有重要意義。

國內(nèi)外學(xué)者采用不同的方法對非點源污染進(jìn)行了研究。鄧華等[8]采用田間試驗法研究在不同的土地利用方式下石盤丘小流域非點源污染的流失變化,并針對總氮、總磷的流失特征分別提出了相應(yīng)的治理措施,為流域的污染治理提供了科學(xué)依據(jù)。但試驗法需要花費(fèi)大量的時間和精力,通常只適合于面積較小的流域。付紹桐等[9]利用“源-匯”理論分析了太原市景觀空間特征對非點源污染的影響,在“源-匯”理論的基礎(chǔ)上,Huang等[10]結(jié)合最小阻力累積(minimum cumulative resistance, MCR)模型對黃河三角洲的非點源污染進(jìn)行了動態(tài)模擬,指出耕地的污染風(fēng)險等級為中度以上的面積占耕地總面積的47.57%。Srinivas等[11]采用多元回歸方法探究了印度Ganges流域的土地利用類型與關(guān)鍵水質(zhì)參數(shù)之間的關(guān)系。也有學(xué)者利用輸出系數(shù)模型[12-14]、潛在非點源污染指數(shù)(potential non-point pollution index, PNPI)模型[15-16]研究了流域非點源污染對不同土地利用方式的響應(yīng),均指出耕地和城鎮(zhèn)用地的污染流失較為嚴(yán)重。經(jīng)驗統(tǒng)計模型對數(shù)據(jù)要求不高,操作簡便[17-18],但沒有反映真實的水文過程和物質(zhì)的遷移轉(zhuǎn)化規(guī)律,不能及時響應(yīng)水文要素的動態(tài)變化[19-20]。因此,探究土地利用方式的時空動態(tài)變化對非點源污染的影響,使用機(jī)理模型更為合適[21]。Shi等[22]、李明濤等[23]分別采用SWAT(soil and water assessment tool)模型、HSPF(hydrological simulation program fortran)模型模擬了景觀格局變化對非點源污染的影響,發(fā)現(xiàn)農(nóng)田的高度集中在一定程度上會加劇非點源污染負(fù)荷。宗敏等[24]以遼寧中部城市群為研究對象,利用SWMM(storm water management model)研究了城市化的快速發(fā)展帶來的城市非點源污染變化。天目湖流域因茶園開發(fā)使土地利用方式發(fā)生了改變,賴正清等[25]利用APEX(air pollutants exposure)-SWAT集成模型研究了該變化對流域非點源污染負(fù)荷的影響,結(jié)果表明茶園用地的污染負(fù)荷貢獻(xiàn)率逐年上升。

目前,研究土地利用方式對流域非點源污染的影響主要是通過對比分析污染負(fù)荷在單一土地利用類型下產(chǎn)出量的差異,而實際中流域的計算單元往往包含多種土地利用類型,并不都以某一類型為主,但考慮土地利用類型的不同組合對流域非點源污染的影響研究還較少。本研究以桃林口水庫流域為研究對象,通過構(gòu)建SWAT非點源污染負(fù)荷模型研究不同土地利用結(jié)構(gòu)對各子流域氮、磷流失率的影響,并定量分析研究區(qū)在2014和2019年兩種土地利用模式下非點源污染產(chǎn)量的變化,以期為桃林口水庫流域合理的土地利用規(guī)劃、水環(huán)境保護(hù)及秦皇島市水資源管理提供科學(xué)依據(jù)。

2 研究區(qū)概況

選擇位于秦皇島市青龍滿族自治縣青龍河上游的桃林口水庫流域為研究區(qū)。青龍河發(fā)源于承德市平泉縣境內(nèi),是灤河的第二大支流。桃林口水庫是國家和河北省“八五”“九五”期間的重點建設(shè)工程,其水質(zhì)情況對當(dāng)?shù)氐纳詈桶l(fā)展都有著重大影響。桃林口水庫流域控制面積為4 949 km2,占青龍河流域總面積的79%,其地理位置為118°37′～119°37′E,39°51′～41°07′N,海拔在92～1 841 m之間,地形總體上呈東北高西南低。桃林口水庫流域地理位置及水系分布如圖1所示。

圖1 桃林口水庫流域地理位置及水系分布

研究區(qū)屬北溫帶半濕潤大陸性季風(fēng)氣候,四季特征變化明顯,晝夜溫差大,無霜期較長,多年平均氣溫為24.7 ℃,年平均降水量為500～700 mm,降水主要集中在7—9月份,約占全年降水量的60%～70%,多年平均徑流量為2.55×108m3。流域內(nèi)的土地利用方式以林地為主,棕壤土和褐土是主要的土壤類型,主要農(nóng)作物為小麥和玉米。根據(jù)實測資料,桃林口水庫流域內(nèi)多年平均總磷(TP)、總氮(TN)濃度分別為0.022和3.934 mg/L,根據(jù)《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》(GB 3838—2002),水庫TP負(fù)荷滿足Ⅱ類水質(zhì)標(biāo)準(zhǔn),TN負(fù)荷為Ⅴ類水質(zhì)標(biāo)準(zhǔn)。桃林口水庫是秦皇島市重要的飲用水源地,流域內(nèi)點源污染得到了有效控制,污染來源主要以非點源形式存在。

3 數(shù)據(jù)來源與研究方法

3.1 數(shù)據(jù)來源與處理

構(gòu)建桃林口水庫流域SWAT模型需要的基礎(chǔ)數(shù)據(jù)包括數(shù)字高程模型(digital elevation model,DEM)、土壤類型、土地利用類型等空間數(shù)據(jù),以及氣象、土壤物理及化學(xué)屬性、污染源以及水文水質(zhì)等屬性數(shù)據(jù),數(shù)據(jù)來源及時空尺度見表1。

表1 模型數(shù)據(jù)庫所需基礎(chǔ)數(shù)據(jù)及其來源

研究區(qū)非點源污染來源主要分為農(nóng)村生活、農(nóng)藥化肥和禽畜養(yǎng)殖3種類型。農(nóng)村生活污染與畜禽養(yǎng)殖污染負(fù)荷產(chǎn)量較為固定,因此根據(jù)兩種污染源的TN、TP和NH3—N的比例定義新的化肥類型,并以連續(xù)施肥的形式輸入SWAT模型[26-27],施肥持續(xù)時間設(shè)置為1 a,施肥頻次設(shè)置為每日1次[28],同時綜合考慮各子流域的人口與面積來確定子流域的污染排放量。對于農(nóng)藥化肥污染,查閱年鑒發(fā)現(xiàn)桃林口水庫流域主要農(nóng)作物為小麥和玉米,小麥在3—4月份播種,7月份收獲;玉米在6—7月份播種,10月份收獲,因此結(jié)合兩種作物的生長時間及所需施肥量,將施肥時間設(shè)置為每年的4—9月,施肥頻次設(shè)置為每月1次。

3.2 研究方法

3.2.1 SWAT模型基本原理 SWAT模型可以模擬流域在不同情況下的水循環(huán)過程,能較好地預(yù)測人類活動對流域內(nèi)水、泥沙以及化學(xué)物質(zhì)的長期影響。水量平衡方程是SWAT模型的基礎(chǔ)[29-30],其表達(dá)式如下:

(1)

式中:SWt為第i天的土壤最終含水量,mm;SW0為土壤初始含水量,mm;t為時間,d;Rday為第i天的降水量,mm;Qsurf為第i天的地表徑流量,mm;Ea為第i天的蒸發(fā)量,mm;Wseep為第i天存在于土壤剖面地層的滲透量和側(cè)流量,mm;Qgw為第i天進(jìn)入主河道的地下水量,mm。

SWAT模型的水質(zhì)模塊來源于QUAL2E模型,能對流域中多種形式的氮和磷進(jìn)行遷移轉(zhuǎn)化模擬,氮和磷的轉(zhuǎn)化分別由氮循環(huán)和磷循環(huán)控制。在SWAT模型模擬污染物的過程中首先對子流域進(jìn)行離散,將其劃分為屬性一致的水文響應(yīng)單元(hydrologic response unit, HRU),基于HRU對污染物的流失與遷移進(jìn)行分析,再匯總到各個子流域,然后模擬污染物在河道及水庫等水體中的遷移轉(zhuǎn)化。

3.2.2 模型參數(shù)率定與驗證 SWAT模型的率定與驗證采用“先徑流后營養(yǎng)物質(zhì)”的順序,選取桃林口水庫上游雙山子水文站2004—2020年實測逐月數(shù)據(jù)對徑流模擬參數(shù)進(jìn)行校準(zhǔn),其中2004—2005年為預(yù)熱期,2006—2016年為率定期,2017—2020年為驗證期。對水質(zhì)參數(shù)的率定和驗證采用桃林口水庫壩上逐月實測水質(zhì)數(shù)據(jù),基于率定好的SWAT水量模塊模擬2019—2020年TN和TP負(fù)荷,2019年為率定期,2020年為驗證期。使用SWAT-CUP軟件的SUFI-2算法對模型參數(shù)進(jìn)行率定和驗證,該算法計算效率高[31],并可以同時分析多個參數(shù)的敏感性。

采用決定系數(shù)(R2)和納什效率系數(shù)(Ens)作為模型適用性評價指標(biāo),其表達(dá)式分別為[32]:

(2)

(3)

3.2.3 土地利用方式對非點源污染影響分析 考慮土地利用結(jié)構(gòu)以及不同研究年份的土地利用變化研究土地利用方式對非點源污染的影響?；诹饔?019年的土地利用數(shù)據(jù)將土地利用類型重分類為林地、耕地、草地、住宅用地和其他5種類型,并利用率定好的SWAT模型對流域的非點源污染進(jìn)行模擬,分別統(tǒng)計各子流域2019年TN、TP負(fù)荷輸出量和5種土地利用類型面積占比,根據(jù)子流域內(nèi)5種土地利用類型的面積占比情況,將子流域劃分為8種不同的土地利用結(jié)構(gòu),見表2。

表2 土地利用結(jié)構(gòu)分類

以2019年為研究年份分析不同土地利用結(jié)構(gòu)對桃林口水庫流域TN、TP負(fù)荷的影響?？紤]到降雨徑流對非點源污染負(fù)荷的影響較大[35-37],選取降雨徑流特征相似的2014和2019年作為研究年份,分析不同年份的土地利用變化對桃林口水庫流域TN、TP負(fù)荷的影響。

4 結(jié)果與分析

4.1 子流域劃分

基于桃林口水庫流域DEM提取的河道水系,將研究流域劃分為29個子流域,如圖2所示。將模型中土地利用類型面積、土壤類型面積以及坡度閾值均設(shè)置為10%,按照相同的土地利用、土壤屬性和坡度組合將研究區(qū)進(jìn)一步劃分為2 124個水文響應(yīng)單元(HRUs)。

圖2 研究區(qū)子流域劃分

4.2 模型率定和驗證

首先對模型參數(shù)進(jìn)行敏感性分析,選取對徑流和水質(zhì)有重要影響的31個參數(shù)進(jìn)行參數(shù)率定,率定參數(shù)敏感性及最優(yōu)值見表3。SWAT模型模擬的率定期和驗證期桃林口水庫流域雙山子水文站徑流量和水庫TN、TP負(fù)荷與實測值對比見圖3,雙山子水文站徑流量模擬效果及桃林口水庫TN、TP負(fù)荷模擬評價結(jié)果見表4、5。

表3 模型參數(shù)敏感性及最優(yōu)值率定結(jié)果

表4 研究區(qū)雙山子水文站徑流量模擬評價結(jié)果

圖3 研究區(qū)徑流量及庫內(nèi)TN、TP負(fù)荷模擬值與實測值對比

由表4可知,率定期桃林口水庫流域徑流量模擬結(jié)果的R2=0.93,Ens=0.92;驗證期R2=0.82,Ens=0.80,說明構(gòu)建的SWAT模型能較好地模擬該流域的水文特征,從圖3(a)也可以看出,模擬的各月平均徑流量能夠反映出該流域徑流的年內(nèi)豐枯變化,并且峰值的模擬效果較好。表5顯示,率定期桃林口水庫TN、TP負(fù)荷模擬結(jié)果的R2分別為0.95、0.94,驗證期R2也均在0.85及以上,表明校準(zhǔn)后的SWAT模型能較好地反映研究流域TN和TP的遷移轉(zhuǎn)化情況。因此,構(gòu)建的SWAT非點源污染模型對研究區(qū)的水文特征及TN和TP負(fù)荷的模擬具有較強(qiáng)的適用性。

表5 桃林口水庫TN、TP負(fù)荷模擬評價結(jié)果

4.3 土地利用方式對非點源污染的影響

4.3.1 不同土地利用結(jié)構(gòu)下子流域的非點源污染特征差異 利用SWAT模型對桃林口水庫流域非點源污染進(jìn)行模擬,得到研究區(qū)2019年各子流域的TN、TP負(fù)荷量,由于各子流域在空間上的污染物施用量存在較大差異,本研究以TN、TP負(fù)荷流失率(流失量/施用量)作為污染指標(biāo),以此反映土地利用結(jié)構(gòu)對子流域污染流失的影響,研究區(qū)不同土地利用結(jié)構(gòu)下各子流域TN、TP負(fù)荷流失率見圖4。

圖4 研究區(qū)不同土地利用結(jié)構(gòu)下各子流域TN、TP負(fù)荷流失率(2019年)

從圖4中可以看出,研究區(qū)土地利用結(jié)構(gòu)為林地(Ⅰ)的子流域數(shù)量最多,TN、TP負(fù)荷在不同土地利用結(jié)構(gòu)下的變化特征基本一致。在單一土地利用類型(Ⅰ、Ⅱ、Ⅲ、Ⅷ)中,子流域的TN、TP負(fù)荷流失率按從大到小排序為耕地(Ⅱ)>草地(Ⅲ)>住宅用地(Ⅷ)>林地(Ⅰ);在多種土地利用類型的組合(Ⅳ、Ⅴ、Ⅵ、Ⅶ)中,土地利用結(jié)構(gòu)為林地-耕地-草地(Ⅶ)的子流域TN、TP負(fù)荷流失率最低,土地利用結(jié)構(gòu)為耕地-草地(Ⅴ)的子流域負(fù)荷流失率最高;與相對應(yīng)的單一土地利用類型的土地利用結(jié)構(gòu)相比,林地的組合能有效降低子流域的污染流失率。

4.3.2 土地利用類型的年際變化對流域非點源污染的影響 2014、2019年研究區(qū)的土地利用轉(zhuǎn)移矩陣見表6。分析表6可知,與2014年相比,2019年的林地和住宅用地面積有所增加,耕地和草地面積相應(yīng)減小,2014年的土地利用類型中分別有25.26 km2的耕地、153.10 km2的草地轉(zhuǎn)化成2019年的林地。表7和圖5為2014、2019年研究區(qū)各土地利用類型面積及TN、TP負(fù)荷量。由表7可以看出,2019年桃林口水庫流域TN、TP總負(fù)荷量較2014年明顯減少,分別減少了34.43%、23.75%,說明林地面積的增加以及耕地和草地面積的減少對研究流域的TN、TP污染負(fù)荷有一定的削減作用。圖5顯示,研究區(qū)不同土地利用類型的單位面積產(chǎn)污量差異較大,林地的產(chǎn)污量最小,耕地的產(chǎn)污量最大;在研究年份(2014、2019年)中林地的平均單位面積TN、TP負(fù)荷量分別為10.73和0.40 kg/km2,而耕地的平均單位面積TN、TP負(fù)荷量分別高達(dá)87.69和4.28 kg/km2。

表6 2014、2019年研究區(qū)土地利用轉(zhuǎn)移矩陣 km2

表7 2014、2019年研究區(qū)各土地利用類型面積與TN、TP負(fù)荷量

圖5 2014、2019年研究區(qū)各土地利用類型面積及單位面積TN、TP負(fù)荷量

5 討 論

桃林口水庫流域SWAT模型對非點源污染的模擬具有較強(qiáng)的適用性,但從水質(zhì)模擬結(jié)果中發(fā)現(xiàn)模擬值較實測值偏小,其主要原因可能為:(1)SWAT模型中需要對子流域土地利用和土壤類型的面積閾值進(jìn)行設(shè)定,因此在模擬過程中會忽視對某些小面積區(qū)域的分析計算,但這些區(qū)域污染輸出量可能會較大;(2)未考慮點源污染排放、大氣沉氮等情況,因而造成了誤差。

隨著近幾年環(huán)境保護(hù)力度的加大,退耕還林措施效果明顯,2019年桃林口水庫流域TN、TP輸出量與2014年輸出量相比有所減小。在單一的土地利用類型下,流域內(nèi)林地的污染物流失量最小,草地次之,耕地最大[38-39],林地中樹木的根系能固持土壤,具有較強(qiáng)的土壤保持能力,且能較好地過濾上方徑流中的淤泥,通過改變徑流對TN、TP負(fù)荷的運(yùn)輸過程,可以在一定程度上抑制污染物的遷移與轉(zhuǎn)化[40-42]。桃林口水庫流域的草地多為牧草地,所涉及的行政區(qū)養(yǎng)殖業(yè)較為發(fā)達(dá),存在較高的牲畜養(yǎng)殖污染,因此以草地為主的區(qū)域污染流失也較嚴(yán)重,這一結(jié)論與Zeiger等[43]在James河流域的研究結(jié)果一致。與相對應(yīng)的單一土地利用類型的土地利用結(jié)構(gòu)相比,林地的組合能有效降低子流域的負(fù)荷流失率,說明林地對非點源污染具有較強(qiáng)的調(diào)控作用。

6 結(jié) 論

根據(jù)桃林口水庫流域水文氣象、水質(zhì)、土地利用、土壤類型等資料,構(gòu)建了流域SWAT非點源污染模型,對TN和TP負(fù)荷進(jìn)行了模擬,分析了不同土地利用方式對流域非點源污染負(fù)荷的影響,主要結(jié)論如下:

(1)采用雙山子水文站2004—2020年實測徑流量和桃林口水庫2019—2020年實測TN、TP負(fù)荷量數(shù)據(jù),對構(gòu)建的桃林口水庫流域SWAT非點源污染模型進(jìn)行參數(shù)率定和驗證。率定期徑流決定系數(shù)R2=0.93,納什效率系數(shù)Ens=0.92,TN、TP負(fù)荷決定系數(shù)R2分別為0.95、0.94,納什系數(shù)Ens分別為0.84、0.89;驗證期徑流決定系數(shù)R2=0.82,納什系數(shù)Ens=0.80,TN、TP負(fù)荷決定系數(shù)R2也均在0.85及以上,說明構(gòu)建的SWAT模型適合桃林口水庫流域氮、磷負(fù)荷特征的模擬。

(2)不同土地利用結(jié)構(gòu)下的子流域非點源污染特征差異較大,TN、TP負(fù)荷的分布規(guī)律相似,土地利用結(jié)構(gòu)為林地的子流域污染流失率最低,耕地的污染流失率最高;與相對應(yīng)的單一土地利用類型相比,林地的組合能有效降低污染流失率。

(3)不同土地利用方式的產(chǎn)污能力差異顯著,林地的單位面積產(chǎn)污量最低,耕地最高;土地利用方式的年際變化對流域非點源污染也有較大影響,與2014年相比,2019年桃林口水庫流域林地面積增加而耕地與草地面積減少,TN、TP總負(fù)荷量分別減少了34.43%、23.75%。