基于VFSMOD 模型的植被缓冲带对千岛湖地区农田面源污染磷负荷削减效果模拟

2022-02-20李芸王斌袁静储昭升金春玲

环境工程技术学报 2022年1期

李芸，王斌，袁静，储昭升*，金春玲

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 所示。建设植被缓冲带时既要保证对污染物的削减效果也要考虑土地资源利用的合理性，所以建议根据不同汇水区长度设计不同的植被缓冲带宽度，做到“一案一策”。可根据期望的削减目标和径流区实际情况确定相应的缓冲带宽度。根据图9，在径流区长度为300 m 的农田下方，需要达到TP 削减率目标为30%时，植被缓冲带的设计宽度建议为25～30 m。该设计曲线可为千岛湖地区针对农田面源磷净化的缓冲带设计提供参考。

图9 丰、平、枯水文年植被缓冲带宽度设计曲线Fig.9 Design curve of buffer zone during the wet,flat and dry years

由于农田地表径流具有流动的不确定性，存在不可控的径流以及未收集径流，因此需要建设布水设施收集更多的农田径流，此外模型模拟时由于参数条件更加理想化，模拟结果会比实际削减率要高，所以对于面积较大的径流区，可以考虑在缓冲带之前设置塘库等调蓄型湿地滞留初期雨水，以提高植被缓冲带对污染物的截留能力，还可以建设碎石床、下凹式绿地等透水渗滤设施，促进来水的下渗[31]。