桃林口水庫流域氮磷污染負荷對土地利用的響應分析

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

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

1 研究背景

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

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

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

2 研究區(qū)概況

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

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

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

3 數據來源與研究方法

3.1 數據來源與處理

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

表1 模型數據庫所需基礎數據及其來源

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

3.2 研究方法

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

(1)

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

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

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

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

(2)

(3)

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

表2 土地利用結構分類

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

4 結果與分析

4.1 子流域劃分

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

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

4.2 模型率定和驗證

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

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

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

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

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

表5 桃林口水庫TN、TP負荷模擬評價結果

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

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

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

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

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

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

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

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

5 討 論

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

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

6 結 論

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

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

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

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