基于改進的SWAT 模型的山地流域吸附態(tài)總磷負荷模擬

龍?zhí)煊澹娚贅s，李業(yè)盛，鄭岐全，劉 敏

（重慶大學 三峽庫區(qū)生態(tài)環(huán)境教育部重點實驗室，重慶 400045）

隨著點源污染管控力度的加大，非點源污染已成為水體污染的主要原因[1]。2015年，國家正式發(fā)布《水污染防治行動計劃》，提出要全面控制工業(yè)、農業(yè)和城鎮(zhèn)污染物排放，深化污染物總量控制。明確非點源污染的時空效應，研究非點源污染的產生和運移機制，識別關鍵污染源區(qū)是開展非點源污染治理的關鍵。非點源污染從形式上可劃分為吸附態(tài)和溶解態(tài)兩種類型。吸附態(tài)非點源污染是指污染物以陸面上被侵蝕泥沙為載體，在徑流的作用下遷移進入水體所形成的污染。磷主要以吸附態(tài)的形式進入水體，在三峽水庫，80%以上的磷是以吸附態(tài)形式存在[2]。

非點源污染具有起源分散、隨機性強、污染范圍廣等特點，因此，監(jiān)測識別非常困難，特別是大型流域或區(qū)域。應用數(shù)學模型進行模擬和分析是研究非點源污染最有效、最直接的方法。通過構建數(shù)學模型，從景觀生態(tài)學和系統(tǒng)動力學角度出發(fā)，研究自然特性和人類活動對水環(huán)境的影響已成為當前非點源污染研究的熱點之一。在流域尺度的非點源污染模型中，SWAT模型被廣泛采用[3-7]。吸附態(tài)污染是陸面上被侵蝕土壤在降雨形成的地表徑流的遷移作用下伴隨著入河泥沙所形成的污染，因此，產沙量的模擬是吸附態(tài)污染負荷模擬的基礎。由于被侵蝕的土壤在遷移過程中會產生部分沉積，因此最終進入水體的泥沙量總是小于被侵蝕的土壤量。在SWAT模型中，產沙量的模擬采用的是修正的土壤流失方程（MUSLE）。在該方程中，涉及到陸面泥沙輸移比（陸面上被侵蝕的土壤進入水體的比例）的因子“峰值流量”和“地表徑流量”，其通常采用小流域出口值，出口值不能反映泥沙輸移比的空間異質性或空間分布，僅能反映整個流域形成的地表徑流隨時間變化所引起的泥沙輸移比的變化，因此模擬結果無法確定陸面上產生泥沙以及吸附態(tài)非點源污染的關鍵區(qū)域，這顯然不利于從源頭控制和有效治理產沙和吸附態(tài)非點源污染。

為此，本文針對山地流域，提出坡面（山地流域的陸面）泥沙輸移比的時空分布式，對SWAT模型中的泥沙模型進行改進，并在此基礎上，構建吸附態(tài)總磷負荷模型，以山地流域，即三峽庫區(qū)包含的御臨河流域、小江流域、大寧河流域和香溪河流域為例，對它們的吸附態(tài)總磷負荷進行模擬。由于三峽庫區(qū)為行政區(qū)域，而不是完整的流域，因此為模擬三峽庫區(qū)產生的吸附態(tài)總磷負荷，采用了小流域擴大法，將庫區(qū)分為御臨區(qū)、小江區(qū)、大寧區(qū)和香溪區(qū)4個小區(qū)，將其中包含的流域負荷模型及其模型參數(shù)擴大到整個小區(qū)進行模擬。

1 研究區(qū)域與方法

1.1 研究區(qū)域及數(shù)據(jù)來源

三峽庫區(qū)（105°44′～110°39′E，28°32′～31°44′N）處于長江上游向中游過渡的平原和四川盆地的結合地區(qū)，面積廣闊，總面積約為5.8萬km2。庫區(qū)包括御臨河流域、小江流域、大寧河流域和香溪河流域。庫區(qū)內地形以丘陵、山地為主，丘陵約占總面積的21.7%，山地約占74.0%[8]。該區(qū)域屬于亞熱帶季風性濕潤氣候，多年平均氣溫17～19℃。降水豐富，年平均降雨量為1 150.26 mm，但其空間分布不均，不同年份相差頗大，降雨主要集中在夏季，多為暴雨或大暴雨，導致庫區(qū)水土流失嚴重。

模擬所需資料包括空間數(shù)據(jù)和屬性數(shù)據(jù)等，具體來源見表1。數(shù)字高程圖（DEM）和重分類的土地利用類型圖見圖1和圖2。

1.2 坡面輸沙模型

1.2.1 坡面產沙模型構建

圖1 三峽庫區(qū)數(shù)字高程圖Figure 1 DEM diagram of the Three Gorges Reservoir area

表1 數(shù)據(jù)及來源Table 1 Data and sources

圖2 三峽庫區(qū)土地利用圖Figure 2 Land use map of the Three Gorges Reservoir area

在陸面上被侵蝕的土壤在降雨徑流的作用下遷移進入水體，形成入河泥沙。由于被侵蝕的土壤在遷移過程中將產生部分沉積，因此最終進入水體的泥沙量小于被侵蝕量。對于山地流域，坡面被侵蝕的土壤進入水體的比例通常稱為坡面泥沙輸移比HSDR（The hillslope sediment delivery ratio），即

式中：HSDR為坡面某單元泥沙輸移比，無量綱；Y為產沙量，t；E為坡面某單元的土壤侵蝕量，t。顯然，不同時刻、不同空間位置的HSDR不同，其隨時空的變化而變化。對于陸面上某空間單元的土壤侵蝕量E，通常采用修正的通用土壤流失方程（RUSLE）進行計算：

式中：R為降雨侵蝕力因子，MJ·mm·hm-2·h-1；K為土壤可蝕性因子，t·h·MJ-1·mm-1；LS為坡長坡度因子，無量綱；C為植被覆蓋度因子，無量綱；P為水土保持措施因子，無量綱。由式（2）和式（3）可得：

在SWAT模型中，采用改進的通用土壤流失方程（MUSLE）計算產沙量。

式中：Q為地表徑流量，m3；qp為峰值流量，m3·s-1。在現(xiàn)有文獻中，qp和Q基本上采用流域出口的實測值。

對比式（4）與式（5）可知，在MUSLE中與泥沙輸移比HSDR以及降雨侵蝕力因子R有關的因子為地表徑流量Q和峰值流量qp，因常采用出口值計算，因此代表的是流域空間上的平均值，不能反映產沙量的空間變異性或空間分布，以及對應的吸附態(tài)非點源污染負荷的空間分布。為此，本文針對山地流域，構建具有時空分布特性的坡面泥沙輸移比HSDR，并在此基礎上形成時空分布的產沙模型以及吸附態(tài)總磷負荷模型，模擬產沙量和吸附態(tài)磷負荷的時空分布，識別產沙和形成污染負荷的“關鍵源區(qū)”。

1.2.2 泥沙輸移比的時空分布式

針對山地流域，在分析和研究了現(xiàn)有有關坡面泥沙輸移比計算式的基礎上，Broselli等[9]提出了“流動連通性”因子IC（Flux connectivity），并由其形成了HSDR的計算式。Broselli等認為，坡面某空間單元上被侵蝕土壤進入水體的比例取決于該單元與坡上單元和坡下單元的水力連通狀況。坡上單元指匯流于該單元的所有單元，而坡下單元指從該單元進入到水體所經過的所有單元。假定坡上單元水力連通的程度與坡上各單元的植被覆蓋率、坡度和總匯流面積成正比，而坡下單元水力連通的程度與坡長成正比、與坡下各單元的植被覆蓋率和坡度成反比，得到：

式中：Dup為坡上單元水力連通分量；Ddn為坡下單元水力連通分量；----Cup為坡上各單元平均植被覆蓋率，無量綱；Sup為坡上各單元平均坡度，無量綱；Aup為坡上單元的總匯流面積，m2；lw為坡下單元w的長度，m；Cw為坡下單元w的植被覆蓋率，無量綱；Sw為坡下單元w的坡度，無量綱。由于0≤HSDR≤1，且HSDR應隨IC非線性地單調增加，因此可近似假定HSDR與IC之間滿足玻爾茲曼S形關系，即

式中：IC0和k IC為參數(shù)，由實測數(shù)據(jù)率定；HSDRmax為HSDR最大值，取決于表層土壤的類型，其數(shù)值范圍為0.7～0.9[10]。研究表明，采用式（7）的坡面泥沙輸移比與式（3）的RUSLE結合模擬的產沙量精度明顯提高[10]。

分析式（6）和式（7）可以發(fā)現(xiàn)，HSDRIC表征了坡面泥沙輸移比的空間變化，即在坡面不同位置被侵蝕的土壤進入水體的比例不同，但無法表征降雨產生的徑流等隨時間變化而導致的變化。鑒于徑流是遷移被侵蝕土壤的動力，且徑流隨時間的變化是引起坡面泥沙輸移比隨時間變化的最重要因素，定量表示徑流強度和徑流隨時間變化的單元累積地表降雨徑流深能綜合反映坡上各單元的植被覆蓋率、坡度和總匯流面積對“流動連通性”的作用，因此，采用徑流深來表征與坡上各單元的水力連通程度應比式（6）更精確，這也解決了“HSDRIC無法表征降雨產生的徑流隨時間變化引起的坡面泥沙輸移比變化”的問題。為此，假定坡上單元水力連通程度與該單元的累積地表降雨徑流深成正比，而坡下單元水力連通程度與式（6）相同，重新定義“流動連通性”因子。新定義的“流動連通性”因子稱為“徑流連通性”RC（Runoff connectivi?ty），RC以及相應的HSDRRC的計算式為：

式中：H為某單元的累積地表降雨徑流深，m；RC0和kRC為參數(shù)，由實測數(shù)據(jù)率定。式（9）為具有時空分布特性的HSDRRC計算式。

1.2.3 模型的驗證

為對比HSDRRC與HSDRIC的模擬效果，應用式（4）和式（9）或式（7）分別對三峽庫區(qū)中的山地流域——大寧河流域和小江流域的產沙量進行模擬，并將模擬結果與實測值進行對比。對于大寧河流域，根據(jù)巫溪站2006—2009年各年中4—10月泥沙觀測值（1—3月為枯水期，產沙量近似為0）對模型參數(shù)進行率定，采用巫溪站2010—2013年各年中4—10月的泥沙觀測值來檢驗HSDR的合理性和精確性。對于小江流域，選用寶塔窩站2002—2003年各年3—10月的泥沙觀測值對模型參數(shù)進行率定，采用2004年3—10月的泥沙觀測值進行檢驗。使用納什系數(shù)（Nash-Sutcliffe effi?ciency，NS）和決定系數(shù)R2評價模型的合理性和精度。

本文式（3）中各因子的計算方法與現(xiàn)有研究成果[11-12]相同，式（8）中累積地表降雨徑流深H的計算采用分布式徑流曲線（SCS-CN）的方法[11]，并借助ArcGIS中加權累積流函數(shù)得出其值，HSDRmax依據(jù)土壤數(shù)據(jù)分別賦值。由于三峽水庫蓄水后支流的流速不大，因此不計河道內河床侵蝕產生的泥沙，坡面入河泥沙在河道內的沉積采用改進的Bagnold泥沙輸移方程進行計算[13]。

表2為HSDRRC和HSDRIC的模擬精度對比。由表2可知，基于HSDRRC的產沙量模擬效果比基于HS?DRIC的模擬效果更佳。對于大寧河流域，在模型檢驗期應用HSDRRC模擬產沙量的NS值和R2值相比HSDRIC提高了14.1%和7.8%。各流域HSDRRC中參數(shù)率定結果見表3。

1.3 吸附態(tài)總磷負荷的時空分布模型

在改進了SWAT模型的產沙量計算方法后，采用SWAT模型中由產沙量計算吸附態(tài)總磷負荷的方法，用式（10）和式（11）計算該單元產生的吸附態(tài)總磷負荷[14]。

式中：LP為吸附態(tài)總磷負荷，t；A為單元面積，hm2；Cs為土壤中磷背景含量，g·kg-1；η為富集比，無量綱。由于改進的產沙模型為時空分布模型，由此形成的吸附態(tài)總磷負荷模型也變?yōu)闀r空分布模型。

1.4 小流域分區(qū)推廣法

三峽庫區(qū)為行政區(qū)域，而非完整的流域，為模擬三峽庫區(qū)的產沙量以及吸附態(tài)總磷負荷的時空分布，參考有關文獻[15-16]，選用小流域分區(qū)推廣法對三峽庫區(qū)進行模擬。首先，對庫區(qū)所包含的4個小流域，即御臨河流域、小江流域、大寧河流域和香溪河流域，分別根據(jù)各自流域出口實測的泥沙數(shù)據(jù)率定其模型參數(shù)（HSDRRC中的RC0和kRC），并對其模擬效果進行驗證（表3），然后將4個小流域模型模擬的區(qū)域分別推廣到相應的御臨區(qū)（YL）、小江區(qū)（XJ）、大寧區(qū)（DN）和香溪區(qū)（XX）4個區(qū)域（圖3），4個區(qū)域覆蓋了整個三峽庫區(qū)。

表2 模擬精度對比Table 2 Comparison of simulating precision

表3 各小流域產沙模型參數(shù)的率定值及模擬效果Table 3 HSDR parameter rate setting and accuracy of each zone

2 結果與分析

2.1 三峽庫區(qū)輸沙量模擬

以2016年為例，4個區(qū)域的產沙量如圖4所示。4個區(qū)域產沙量的空間差異明顯，產沙最多的區(qū)域為大寧區(qū)，其東北部為產沙的關鍵源區(qū)，該區(qū)域單位面積產沙量峰值為2 154.26 t·hm-2，平均值為18.03 t·hm-2；產沙最少的區(qū)域為香溪區(qū)，其峰值為1 245.27 t·hm-2，平均值為5.98 t·hm-2。

2.2 三峽庫區(qū)吸附態(tài)總磷負荷模擬

2.2.1 負荷的空間分布

圖3 三峽庫區(qū)模擬分區(qū)Figure 3 The Three Gorges Reservoir area division

依據(jù)小流域擴大法匯總獲得的三峽庫區(qū)產沙量分布，由式（10）可以模擬出三峽庫區(qū)吸附態(tài)總磷負荷，2016年的模擬結果如圖5所示。從圖5中可以看到，總體上負荷的空間變化較大，小江區(qū)北部負荷最高，最高達到144.6 kg·hm-2，此外，大寧區(qū)部分地域和御臨區(qū)東部的負荷也較高，這些地區(qū)是負荷形成的關鍵源區(qū)。表4為各分區(qū)總磷負荷量和平均負荷強度（單位面積上的負荷量）?？傮w來說，大寧區(qū)產生的污染負荷最多，其次為小江區(qū)。整個庫區(qū)平均吸附態(tài)總磷負荷強度為0.49 kg·hm-2。

圖4 4個區(qū)域的產沙量分布Figure 4 Distribution of slope sediment transport capacity in four regions

2.2.2 負荷隨時間的變化

由于庫區(qū)各年以及同一年中各月的降雨量不同，導致產沙量和吸附態(tài)總磷負荷隨時間發(fā)生變化。圖6為2007—2016年庫區(qū)產沙量和總磷負荷的年際變化，圖7為2007年4月—2017年6月庫區(qū)總磷各月間的變化。從圖6中可以看到，產沙量和吸附態(tài)總磷年際差異顯著，產沙量峰值出現(xiàn)在2016年，約6 020萬t，最小的年份為2012年，約3 610萬t；吸附態(tài)總磷負荷峰值出現(xiàn)在2014年，約2 890 t，最小為2012年，約1 750 t。吸附態(tài)總磷負荷不僅年際差異大，由于同一年內各月降雨量不同，各月的差異也十分明顯。由于每年1—3月及11—12月降雨少、雨量低，基本不產生地表徑流，產沙量和總磷負荷可認為是0。從圖7中可以看到，全年的污染負荷主要來源于夏季。

3 結論

圖5 三峽庫區(qū)2016年吸附態(tài)總磷負荷分布Figure 5 Distribution of adsorbed phosphorus loads in the Three Gorges Reservoir area in 2016

表4 各分區(qū)吸附態(tài)總磷負荷量及負荷強度Table 4 Adsorbed phosphorus load and load intensity of each region

圖6 三峽庫區(qū)產沙量和吸附態(tài)總磷年際變化Figure 6 Interannual variation of slope sediment transport capacity,adsorbed phosphorus

圖7 三峽庫區(qū)各月吸附態(tài)總磷的變化Figure 7 Monthly variation of adsorbed phosphorus in the Three Gorges Reservoir area

（1）針對山地流域，提出“徑流連通性”因子，基于“徑流連通性”因子形成坡面泥沙輸移比的時空分布式，對SWAT模型中產沙量的計算方法進行了改進，構建了時空分布的山地流域坡面產沙模型以及吸附態(tài)總磷負荷模型。產沙量的模擬結果表明，改進后的模型不僅能模擬產沙量的時空分布，而且模擬精度高。

（2）應用所構建的吸附態(tài)總磷負荷模型，采用小流域推廣法，對三峽庫區(qū)進行分區(qū)模擬，實現(xiàn)了對吸附態(tài)總磷負荷的模擬。模擬結果顯示三峽庫區(qū)的吸附態(tài)總磷負荷的時空差異性大，在空間上，負荷強度從0到144.6 kg·hm-2，在大寧區(qū)和小江區(qū)，平均負荷強度較大，分別為 0.81 kg·hm-2和 0.69 kg·hm-2；同一年內各月的吸附態(tài)總磷負荷差異明顯，全年的污染負荷主要來源于夏季。

（3）本文基于“徑流連通性”因子所形成山地流域產沙量和吸附態(tài)負荷時空分布的計算方法可用于其他山地流域，通過模擬，獲得產沙和產生吸附態(tài)負荷的關鍵區(qū)域，有利于吸附態(tài)負荷的有效治理。