基于VFSMOD 模型的植被緩沖帶對千島湖地區農田面源污染磷負荷削減效果模擬

李蕓，王斌，袁靜，儲昭升*，金春玲

1.核資源與環境國家重點實驗室,東華理工大學

2.中國環境科學研究院湖泊生態環境研究所

3.東華理工大學水資源與環境工程學院

湖泊是地球生態環境的重要組成部分，湖泊生態健康受到威脅會影響到人類的健康以及社會可持續發展。研究表明，近年來農業面源污染已成為威脅湖庫水質安全的重要因素之一，并導致湖庫水體富營養化加劇[1]。其中，農業面源磷流失對湖庫富營養化具有重要影響。研究認為，磷是藻類生長的第一限制因素[2]。植被緩沖帶是攔截入湖污染的最后一道生態屏障，因其對地表徑流攜帶的泥沙、氮與磷、農藥以及其他無機污染物有著良好的攔截效果，已經被歐美國家作為控制農業非點源污染的有力措施之一[3]。因此在湖濱區域構建植被緩沖帶已成為入湖徑流污染控制技術研究的熱點。

VFSMOD 模型是目前比較先進的、可全面模擬植被緩沖帶對地表徑流污染物攔截以及泥沙削減的工具，其受到美國農業部《緩沖帶、廊道和綠色通道設計指南》的推薦[4]。Mu?oz-Carpena 等[5]比較了27 場自然降水條件下植被緩沖帶對泥沙的截留試驗結果與VFSMOD 模型模擬結果，發現精確度最高達87%，模擬效果較好。Dosskey 等[6]通過VFSMOD模型對泥沙和徑流進行模擬，從泥沙和水的模型關系解釋了泥沙結合和溶解污染物的關系，為植被緩沖帶設計與建設提供了參考。楊寅群等[7]通過對比野外試驗植被緩沖帶凈化效果及模型模擬情況，發現VFSMOD 模型對國內緩沖帶模擬具有較高的精度與較強的適應性，植被緩沖帶出流流量及泥沙濃度模擬與實測系數分別為0.995 和0.889，可用于國內植被緩沖帶寬度的劃定。目前關于VFSMOD 模型的研究主要集中在植被緩沖帶對徑流泥沙攔截與削減效率的模擬和分析，但對徑流污染負荷削減效益研究不足。

筆者以千島湖地區為例，模擬了植被緩沖帶寬度、坡度以及降水量變化對緩沖帶削減入湖徑流磷污染負荷效果的影響，并根據千島湖地區不同長度徑流區的模擬結果，提出植被緩沖帶設計工作曲線，以期為湖濱區植被緩沖帶的規劃與設計提供參考。

1 研究區概況

千島湖（118°34′E～119°15′E，29°22′N～29°50′N）位于浙江省西北部的杭州市淳安縣。區域年平均氣溫為16.3～17.2 ℃，年平均降水天數為156 d，多年平均降水量為1 515.1 mm，降水量年內分布不均衡，多雨季（4—9 月）降水量占全年的68.0%。土壤種類以紅壤、黃壤、巖性土、水稻土為主，其中紅壤面積最大（412.1 hm2），占淳安縣總面積的62.1%[8]，黏質土壤占22.56%，砂質占28.41%，礫質占20.77%。土壤肥力中等，主要種植作物為水稻、小麥、玉米、大豆等。千島湖西南部汾口鎮、大墅鎮、界首鄉等地集中有農田、茶園、村落等典型農業污染源，且地形較為典型，因此選擇該處為典型片區（圖1）進行模型模擬。

圖1 千島湖典型片區范圍及匯水區入湖分布概況Fig.1 Typical area of Qiandao Lake and distribution of catchment area into the lake

2 研究方法

2.1 VFSMOD 模型

VFSMOD 模型是一款用于預測植被緩沖帶對地表徑流泥沙及污染物凈化的大規模尺度機理模型[9]，該模型系統基于研究區緩沖帶坡度、植被等地表條件以及復雜的降水類型和降水強度，模擬泥沙沉積所導致的徑流中污染物濃度的變化過程[6]。VFSMOD 模型由模擬植被緩沖帶內運移的水文模型(模擬水流)和泥沙輸移模型(模擬泥沙)構成，具體包括4 個模塊：1）滲透模塊，主要用來計算表層土壤水量平衡；2）地表徑流模塊，主要用來計算滲透性土壤表面徑流深和徑流量；3）泥沙沉淀過濾模塊，模擬入流泥沙量在植被緩沖帶內的輸移與沉積；4）污染物轉移模塊。模型的組成結構如圖2 所示[10]。

圖2 VFSMOD 模型結構示意Fig.2 Schematic diagram of VFSMOD model structure

2.2 VFSMOD 模型主要參數確定

VFSMOD 模型需要輸入的參數分為5 類：地表徑流模擬參數、降水和入流參數、入滲模型土壤參數、植被緩沖帶性能參數以及泥沙特征參數。

2.2.1地表徑流模擬參數

地表徑流參數主要包括緩沖帶長度、寬度、坡度，各段曼寧糙率系數等。通過對千島湖典型片區的匯水區進行GIS 分析，徑流區長度主要集中在100～600 m，極少數沿湖匯水區以及一級支流匯水區的徑流長度小于100 m 或大于600 m，因此在對比不同徑流區削減效果時，徑流區長度選擇100～600 m。在0～30 m 內選取6 種緩沖帶寬度進行模擬。模擬過程中，根據緩沖帶表面性質不同，將緩沖帶分為若干小段，每段長度設為1 m。研究坡度對削減效果的影響時，考慮到淳安縣耕地主要分布在坡度小于25°的低洼平坦地區，坡度大于25°的耕地面積僅占全縣耕地面積的4%[11-12]，結合楊寅群等[13]提出的坡度大于60°的地表徑流難以保持薄層水流狀態，將設計坡度定為3%～30%。根據現場踏勘結果，匯水區坡度主要集中在2%～4%（圖3），因此，在研究降水量變化對削減效果的影響及不同長度徑流區模擬對比時，匯水區坡度均取3%。緩沖帶各段曼寧糙率系數與緩沖帶內植被狀況有關，本研究模擬緩沖帶內植被為濃密草地，糙率系數取0.24 s/cm1/3。

圖3 典型片區匯水區坡度分布Fig.3 Slope distribution of the catchment area in typical areas

2.2.2降水和入流參數

通過國家氣象科學數據中心（http://data.cma.cn/）獲取淳安縣1989—2018 年降水、氣象數據，根據降水量等級標準[14]對典型平水年日降水狀況進行強度分級，結果如表1 所示。

由表1 可知，年降水量超過30 mm 的天數為15 d，占年降水天數的10%，其累計降水量占年降水量的44.46%。日降水量可能為一場降水或多場降水事件的降水量總和，根據千島湖地區降水特征，假設日降水量為30 mm 的降水事件是6 h 的單場降水事件，其降水量過程如圖4 所示（徑流區長度為100 m，入流峰值為0.027 m3/s）。

圖4 VFSMOD 模型降水量過程線Fig.4 Precipitation hydrograph of VFSMOD model

表1 典型平水年不同強度降水天數分布情況Table 1 Distribution of precipitation days with different intensities in typical flat water years

2.2.3入滲模型土壤參數

入滲模型土壤參數主要包括飽和導水率（VKS）、濕潤鋒平均吸力（SAV）、初始含水量（OI）和飽和含水量（OS），其中VKS 被認為是最敏感的參數。Mu?oz-Carpena 等[15]研究發現，VKS 是控制徑流量、緩沖帶滲透性以及泥沙傳輸的主要因素。如前所述，千島湖地區主要以壤土為主，因此模擬過程中，參照Rawls 等[16]研究給出的VFSMOD 模型參考參數，選擇VKS 為3.67×10-6m/s，SAV 為0.089 9 m。初始含水量和飽和含水量使用采集的過濾土樣實測平均值，分別為0.200 和0.499 m3/m3。

綜上所述，雖然Logistic模型是針對于二值變量的回歸，但不論是從現實意義，還是統計檢驗甚至整體預測擬合，逐步回歸模型明顯更貼近實際生活．

2.2.4植被緩沖帶性能參數

植被緩沖帶內植被類型以及植被分布密度也是影響緩沖帶削減效果的關鍵因素。本研究中緩沖帶內植被選用混合草本，參考美國典型緩沖帶植被參數[10]，將植物莖干間距設為2.15 cm，植物高度設為18 cm，緩沖介質修正糙率（VN）和泥沙淤滿緩沖帶后裸露表面的糙率（VN2）對模型輸出的結果不敏感，VN 取0.012，VN2 取裸露黏壤土的糙率0.02。

2.2.5泥沙特性參數

入流泥沙顆粒分級數（NPART）參考美國農業部（USDA）泥沙粒徑分級標準[10]。模擬中根據入流泥沙類型設定NPART，由于泥沙特性為細顆粒聚合體，NPART 取3。入流泥沙濃度采用農田徑流實測濃度，為0.002 g/cm3。模型具體參數匯總見表2。

表2 VFSMOD 模型模擬參數Table 2 VFSMOD model simulation parameters

2.3 徑流磷污染負荷削減估算

目前VFSMOD 模型的輸出機制主要模擬出水及泥沙輸移，為了評估磷污染負荷的削減，需要結合試驗和經驗公式進行估算。磷在地表徑流中流失按物理狀態可分為顆粒態磷和溶解態磷[17]，植被過濾帶徑流中顆粒態污染物使其隨泥沙沉淀而得以削減。袁溪等[18]發現徑流中磷素流失量與SS 流失量具有顯著線性關系，約占36%。Schmitt 等[19]研究發現，溶解態磷削減率約占徑流削減率的41%～60%。

根據已有研究經驗，結合千島湖實地情況，本研究在模擬磷負荷削減的估算中，泥沙中磷削減率確定為泥沙削減率的40%；徑流中磷削減率確定為徑流削減率的60%。估算TP 削減效率時，選用以下公式進行計算：

式中：LTP為地表徑流中TP 的削減率，%；Q泥沙為泥沙的削減率，%；Q徑流為徑流的削減率，%。

3 結果與分析

3.1 緩沖帶寬度對削減效果的影響

經過模型模擬和估算，不同寬度的植被緩沖帶對泥沙以及TP 負荷削減效果如圖5 所示。由圖5可知，植被緩沖帶對泥沙負荷的削減率隨著緩沖帶寬度的增大而增大。TP 削減率的變化趨勢與泥沙負荷削減率變化規律一致，這是因為顆粒態磷在TP 中的占比較高。此外，植被緩沖帶對泥沙和TP 的削減更多集中在前端區域，當緩沖帶寬度為10 m 時，對泥沙削減率達85%，表現出良好的泥沙削減效果。Wanyama 等[20]在維多利亞湖濱帶削減試驗中同樣發現，在自然降水條件下超過70%的泥沙會在植被緩沖帶的前5～10 m 被攔截。隨著寬度的增加，植被緩沖帶對泥沙的削減率增幅逐漸變緩：當緩沖帶寬度從5 m 增至10 m 時，泥沙削減率提高了20 個百分點；而當緩沖帶寬度從10 m 增至15 m 時，泥沙削減率僅提高了6 個百分點。國內許多學者在研究植被緩沖帶對面源污染影響時同樣發現，緩沖帶對徑流污染物削減能力在前段的10～15 m 較強[21-23]，這與本研究結果相似。然而對比TP 和泥沙削減率發現，相同削減率下TP 所需的削減距離比泥沙大，當植被緩沖帶寬度為5 m 時，泥沙削減率達到67.1%，而TP 削減率僅為30%；TP 削減率達到60%時，緩沖帶寬度為30 m，此時泥沙削減率為97.8%。這是因為磷素結合的顆粒物更為細小，顆粒態磷中有32%是通過小于0.045 mm團聚體流失的[24]。

圖5 不同寬度植被緩沖帶對泥沙和TP 的削減效果Fig.5 Reduction effects of different width vegetation buffer zones on sediment and TP

3.2 坡度對削減效果的影響

坡度對植被緩沖帶削減效果的影響如圖6 所示。由圖6 可知，總體上植被緩沖帶對泥沙和TP 的削減率均隨著坡度的增大而逐步減小。這是由于坡度越小，徑流的流速越小[25]，徑流通過植被緩沖帶的時間則越長，緩沖帶削減泥沙和污染物的效果也就越好[26]。坡度為3%～30%時，植被緩沖帶削減凈化率降幅較小，泥沙削減率降幅不超過1 個百分點，TP 削減率降幅不超過1.76 個百分點。這可能是因為模擬中未改變污染物上源區坡度，污染物負荷不變，植被緩沖帶在此坡度范圍內仍能保證徑流處于薄層漫流的狀態，致使植被緩沖帶坡度改變對削減性能的影響較弱，這與國內其他研究者對植被緩沖帶坡度研究結論相似[26-28]。

圖6 不同坡度植被緩沖帶對泥沙和TP 的削減效果Fig.6 Reduction effect of different slope vegetation buffer zones on sediment and TP

3.3 降水量對削減效果的影響

根據千島湖年降水量分布特征，降水量選取10～90 mm，降水量對泥沙和TP 削減效果的影響如圖7 所示。由圖7 可知，植被緩沖帶的削減效果與降水量的變化趨勢相反，即隨著降水量的增加，植被緩沖帶對泥沙和TP 削減率均減小，且降水量變化對植被緩沖帶削減效果的影響較大。降水量為30 mm 時，30 m 寬度的植被緩沖帶泥沙削減率最高（97.8%），5 m 寬度的植被緩沖帶泥沙削減率最低（67.1%）；當降水量增至50 mm 時，30 m 寬度的植被緩沖帶泥沙削減率最高（84.3%），但泥沙削減率降低了13.5 個百分點，5 m 寬度的植被緩沖帶對泥沙削減率最低（33.3%），泥沙削減率降低了33.8 個百分點。這是因為降水量的變化會影響入流流量的變化，進而影響到緩沖帶內徑流流速和徑流挾沙能力的變化，降水量增大，緩沖帶內徑流流速和徑流挾沙能力均變大，植被對徑流的削減效率降低導致僅少部分的泥沙和顆粒態污染物通過沉積作用被削減[20]。鄧娜等[29]研究發現，將進水流量 從7.10×10-3～7.80×10-3m3/s 降 至1.00×10-3～2.90×10-3m3/s 時，植被緩沖帶對SS、TP、顆粒態磷的凈化效果顯著增強。

圖7 不同降水量條件下植被緩沖帶對泥沙和TP 的削減效果Fig.7 Reduction effects of vegetation buffer zone on sediment and TP under different rainfall conditions

對比6 種寬度的植被緩沖帶發現，當降水量小于15 mm 時，植被緩沖帶對泥沙的削減效果基本相同，這是因為降水量小，植被緩沖帶可在前5～10 m對徑流和泥沙進行削減。而當降水量增至30 mm時，不同寬度的植被緩沖帶對泥沙削減率的降幅差異顯著，5 m 寬度的植被緩沖帶對泥沙的削減率下降了32.6 個百分點，10 m 寬度的下降了13.8 個百分點，而30 m 寬度的僅下降2.2 個百分點。所以在植被緩沖帶設計過程中，應分析不同類型的降水產生的徑流量，徑流量增大時可以考慮優化徑流進入緩沖帶的布水方式，通過增加緩沖帶寬度和調整植被密度等方式保證植被緩沖帶對泥沙和TP 的削減效果。

3.4 不同長度徑流區模擬削減效果對比

在研究區內選擇100、300、600 m 3 種長度的徑流區模擬對比削減效果，根據千島湖典型平水年降水特征，降水量取15～90 mm，模擬結果如圖8 所示。由圖8 可知，徑流區越長，緩沖帶對泥沙和TP 的削減效果越差。徑流區長度為300 m 時，5～10 m 寬度的緩沖帶在大于15 mm 的降水量下泥沙削減率最高達53.9%，比相同情況下徑流區長度為100 m 時降低32.2 個百分點；TP 的削減率最高為24.02%，比相同情況下徑流區長度為100 m 時降低16.6 個百分點。而相同情況下徑流區長度為600 m時泥沙與TP 的削減率最高僅為34.00%和15.46%，最低不足3%。徑流區長度較大的農田地表徑流往往集中在農田內部的自然排水道中，并以集中流的形式穿過植被緩沖帶，在這種集中流動的情況下入流流量變大，且植被緩沖帶的有效削減徑流的面積變小，所以植被緩沖帶的削減效果并不理想[30]，需要輔以其他攔截凈化措施。

圖8 不同徑流區長度下植被緩沖帶對泥沙和TP 的削減率Fig.8 Reduction rate of silt and TP in the buffer zone under different runoff zone lengths

此外，徑流區長度為100 m 時，30 mm 降水量下植被緩沖帶寬度從5 m 增至10 m 時，泥沙削減率增加了20 個百分點左右；而當植被緩沖帶寬度從10 m增加到15 m 時，泥沙削減率僅增加6 個百分點；繼續增加緩沖帶寬度，泥沙削減率幾乎不變。徑流區長度為300 m時，相同情況下植被緩沖帶寬度從5 m增加到10 m，泥沙削減率增加了20 個百分點；緩沖帶寬度從10 m增至15 m 時，泥沙削減率增加了14 個百分點；繼續增加緩沖帶寬度至30 m 時，泥沙削減變化速率逐漸趨于平緩。徑流區長度為600 m時，相同情況下植被緩沖帶寬度從5 m 增至10 m時，泥沙削減率增加了14 個百分點，緩沖帶寬度從10 m 增至15 m 時，泥沙削減率增加了12 個百分點，繼續增加緩沖帶寬度至45 m 時，泥沙削減變化速率逐漸趨于平緩。綜上，徑流區長度增加時，應考慮增加緩沖帶的寬度來提高植被緩沖帶對泥沙和TP 的削減效率。

3.5 研究區植被緩沖帶設計方案建議

根據研究區單場降水緩沖帶設計曲線特征，結合千島湖年內降水量分布情況（表1），模擬年度植被緩沖帶設計曲線，具體公式如下：

式中：P為TP 的年削減率，%；Q0為典型水文年的日降水量中不同降水區間的降水量之和，mm；P0為典型水文年的日降水量中不同降水區間對應的TP 削減率，%；Q年為年降水量，mm；β為徑流收集率，本研究取60%。

在農田徑流收集率為60%條件下模擬千島湖地區緩沖帶寬度設計曲線，結果如圖9 所示。建設植被緩沖帶時既要保證對污染物的削減效果也要考慮土地資源利用的合理性，所以建議根據不同匯水區長度設計不同的植被緩沖帶寬度，做到“一案一策”?？筛鶕谕南鳒p目標和徑流區實際情況確定相應的緩沖帶寬度。根據圖9，在徑流區長度為300 m 的農田下方，需要達到TP 削減率目標為30%時，植被緩沖帶的設計寬度建議為25～30 m。該設計曲線可為千島湖地區針對農田面源磷凈化的緩沖帶設計提供參考。

圖9 豐、平、枯水文年植被緩沖帶寬度設計曲線Fig.9 Design curve of buffer zone during the wet,flat and dry years

由于農田地表徑流具有流動的不確定性，存在不可控的徑流以及未收集徑流，因此需要建設布水設施收集更多的農田徑流，此外模型模擬時由于參數條件更加理想化，模擬結果會比實際削減率要高，所以對于面積較大的徑流區，可以考慮在緩沖帶之前設置塘庫等調蓄型濕地滯留初期雨水，以提高植被緩沖帶對污染物的截留能力，還可以建設碎石床、下凹式綠地等透水滲濾設施，促進來水的下滲[31]。

4 結論

（1）VFSMOD 模型模擬結果表明，緩沖帶寬度與徑流中泥沙和TP 的削減率呈正相關，對于污染物削減集中在前段10 m 處，當緩沖帶寬度達到15 m時，削減效果趨于平穩。

（2）植被緩沖帶坡度、降水量與徑流中泥沙及TP 的削減率呈負相關，因本研究中緩沖帶設置在沿湖平坦區域，模擬過程中未改變污染物上源區坡度，使得緩沖帶坡度從3%增至30%時，污染物削減率降幅為1%～1.76%，坡度對緩沖帶削減效果影響不大。

（3）根據不同長度徑流區提出千島湖地區基于豐、平、枯水文年植被緩沖帶寬度的設計工作曲線，可為千島湖地區針對農田面源磷凈化的植被緩沖帶設計提供參考。

（4）當徑流區長度較大、植被緩沖帶削減效果不理想時，可以考慮在緩沖帶之前建設調蓄型濕地等輔助措施，截留凈化初期雨水，以提高緩沖帶對污染的截留能力。