天目湖流域茶园非点源污染多尺度模拟

赖正清，肖雪纯，李硕*，费国松，朱立国，闫浩，夏玉林

（1.南京师范大学地理科学学院，教育部虚拟地理环境重点实验室，南京 210023；2.南京师范大学海洋科学与工程学院，南京210023；3.江苏省水文水资源勘测局常州分局，江苏 常州 213000；4.江苏省水文水资源勘测局无锡分局，江苏 无锡 214000）

土地利用变化是非点源污染的重要影响动因，会对流域的水量与水质产生强烈的影响[1-3]。科学地分析土地利用变化对非点源污染的影响，对提高流域水环境质量以及实现流域生态环境持续发展具有重要的研究意义[4]。我国东南丘陵山区优质耕地缺乏，适宜建设的土地资源具有局限性，地方政府出于缓解人地矛盾的考虑，相继出台了引导和鼓励丘陵山区农业开发的政策[5-6]。江苏溧阳天目湖流域是溧阳市近80万人口的饮用水源地和国家级旅游度假区，自21 世纪初开始，当地政府实施天目湖5A级景区的打造，在流域内部实施了以乡村旅游和鼓励茶果园种植生产为主的丘陵农业开发，迄今逐渐形成了集供水、旅游、农业于一体的区域发展模式[7-10]。天目湖流域丘陵山区茶园的规模化开发显著促进了该区域的社会经济发展，但大范围茶园种植伴随着的水土流失加剧、施肥量剧增形成的非点源污染对下游水源地水质安全造成极大威胁[11-12]。2015—2018 年周监测资料显示，沙河水库和大溪水库共5 条入库支流的总氮（TN）都达到劣Ⅴ类水质标准，大溪水库水质TN 指标仅为Ⅳ类水质，而沙河水库在2017 年和2019 年均有超过Ⅳ类水质标准的情况出现[13]。两个水库现状水质与《常州市天目湖保护条例》提出的全湖Ⅱ类水质目标有较大差距。2015 年以来，地方政府逐渐认识到现有开发模式的环境负面效应，因此陆续出台了一系列的政策，实施了退耕还林、茶果园收储等相关措施[14]。从前期的大力支持开发，到现阶段的补救措施，政策必然会对当地经济发展产生很大程度的影响。如何对规模化茶园开发产生的长期环境效应进行科学评价，水源地保护、茶园开发和流域生态文明建设如何协同发展，成为该区域及类似发展模式地区生态文明建设中迫切需要解决的现实问题。

得益于计算机技术和“3S”技术的发展，国内外学者将地理过程模型与遥感技术、地理信息技术、计算机技术相结合，研发出了多种流域分布式非点源污染过程模型[15-18]，这些模型成为流域非点源污染时空过程定量评估的重要工具。各类非点源污染过程模型开发的背景和应用对象不同，导致其适用性和功能存在一定差异。SWAT 模型（Soil and Water Assessment Tool model）是由美国农业部农业研究中心开发的流域尺度非点源污染模型，主要用于预测较大流域复杂多变的土壤类型、土地利用方式和管理措施条件下土地管理对水分、泥沙和化学物质的长期影响[15]。该模型在计算中采用流域-子流域-水文响应单元的空间离散方式，子流域内部的水文响应单元采用统计叠加的方式生成，因此没有明确的地块划分。该模型的模拟产出是将各个水文响应单元计算结果进行统计汇总后作为子流域的单独产出，进行河道传输的演算，因此在子流域内部无法进行地块尺度的空间确定性模拟输出。APEX 模型（Agricultural Policy/Environ⁃mental eXtender model）由美国德州农工大学在EPIC模型（Environmental Policy Integrated Climate model）的基础上开发，具备对中小型流域或复杂结构农场的土地管理影响的模拟能力，可配置灌溉、排水、沟渠、缓冲带、梯田、施肥、水库、轮作、农药施用、放牧和耕作等多种土地管理策略[19-20]。APEX 模型以子区即同一性质的地块单元为最小离散单元，不同子区通过流域内的水流方向相互联系，可从子流域层面对地块产出进行详尽模拟[21]，但该模型缺乏子流域内部和子流域之间河道过程的模拟，因而不具备对较大流域处理的能力。

模型集成是实现不同非点源污染模型之间的优势互补、提高模拟能力与效率的有效手段[22-25]。因此，本研究首先利用遥感手段揭示天目湖流域土地利用的时空变化特征，解析土地利用的数量结构和空间格局，再通过集成SWAT 模型和APEX 模型，构建基于“源-汇”的一系列要素和过程的多尺度集成模拟方法，定量解析天目湖流域低山丘陵区以茶园为主要土地利用类型的非点源污染产出对水环境的影响，探讨并制订丘陵山地茶园开发影响下的水源地水生态保护调控方案，推动区域经济和生态河湖建设的和谐发展。

1 材料与方法

1.1 研究区概况与基础数据

天目湖流域（图1）属江苏省常州市溧阳市，位于苏浙皖三省交界，属天目山余脉的丘陵地区，介于31°07′~31°20′N、119°21′~119°29′E 之间，流域面积约245 km2，属北亚热带季风气候，雨量充沛，四季分明，年平均气温16.6 ℃，最高气温38.7 ℃，最低气温-4.6 ℃，年平均降水量1 251 mm。天目湖流域主要由沙河水库流域和大溪水库流域构成，拥有沙河、大溪两座大（二）型水库，是太湖流域上游重要的水源涵养区。

图1 天目湖流域概况Figure 1 An overview of the Lake Tianmu basin′s topography with stream networks and reservoirs

天目湖流域两大水库上游的低山丘陵主要来源于天目山余脉，丘陵区从东、南、西3 个方向环绕水库，具有阶梯式地貌特征。天目湖流域的地带性土壤为黄棕壤，其呈微酸性，主要分布在低山山顶以及缓坡，保水性较差，易造成水土流失现象。为了改善当地经济，缓解人地矛盾，20 世纪末当地政府采取了鼓励大规模种植茶叶、果树和以水库为重点的旅游景点建设等政策，流域内越来越多的丘陵坡地被改造用于农业生产。大规模的茶园开发以及茶叶种植中大量的氮肥施用，导致天目湖流域两大水库的TN 始终处于Ⅳ类至劣Ⅴ类标准，难以达标，严重威胁了流域生态环境的可持续发展和居民饮用水的安全。

本研究使用的基础数据主要包括地形、土地利用和土壤等，具体见表1。根据模型输入要求分别构建天目湖流域气象、土壤类型、土地利用类型和农业管理措施参数库。

表1 模型输入数据Table 1 Model input data

1.2 天目湖流域土地利用变化分析

对天目湖流域 2000、2005、2010、2015 年以及2019年5期土地利用类型数据进行统计，获得各期不同土地利用类型面积、占比和变化面积，见表2。从表2 中可以看出，2000 年到2019 年以来天目湖流域的土地利用面积呈现出较大的变化。2000 年时，溧阳市政府尚未对茶、果园等经济作物进行规模化的开发种植，两者仅呈零星状分布，分布在盆地边缘的林地占据了天目湖流域较大的土地利用面积（110.86 km2），其次为农田（104.01 km2），两者的面积占比分别为47.10%和44.19%。2005 年开始，该流域的土地利用发生了较为显著的变化，林地、农田开始向茶、果园转变，林地面积减少了15.39 km2，农田更是大面积减少至55.39 km2，原来分布在水库周边的部分林地也已经被逐渐开发转变为可持续带来经济收益的茶、果园，茶园和果园面积分别达到了9.91 km2和47.10 km2。自2005 年以后，天目湖流域内的农田、林地、草地面积呈现出不断减少的趋势，分别由2005 年的55.39、95.47、1.75 km2下降至 2019 年的 21.84、85.18、0.32 km2。果园和茶园的面积持续增加，果园面积由2005 年的 47.10 km2增加至 2019 年的 57.74 km2，茶园面积也由9.91 km2扩张至21.72 km2，果园以及茶园15 a 的面积扩张速度分别达到了每年近0.7 km2和0.8 km2。此外，随着人口激增带来的用地需求的增加，使居民区也呈现出上升的趋势。纵观2000—2019 年的土地利用变化数据，农田、林地呈现出下降的趋势，面积占比分别由2000 年的44.19%和47.10%下降至2019 年的8.95%和34.90%；而果园和茶园呈现出不断上升的趋势，面积占比分别由2000 年的0.03%和0.16%上升至2019年的23.66%和8.90%。

表2 天目湖流域土地利用结构变化Table 2 Land use change in Lake Tianmu basin

纵观2000—2019 年天目湖流域近20 a 茶园面积变化的统计分析，发现其时间变化特征呈现出由快速扩张向缓慢增长转变的趋势，开发模式由规模化开发逐渐向限制保护性开发转变。

从空间尺度上看（图2），2000 年时天目湖流域的茶园仅零星分布在洙漕河附近，开始规模化开发建设茶园后，茶园主要分布在大溪水库沿湖东岸、洙漕河流域、沙河水库沿湖西岸以及平桥河流域。2005—2019 年，由于政策的驱动和茶园高效益产出所吸引的企业投资，位于沙河西岸、洙漕河流域、平桥河流域的茶园在原有基础上迅猛扩张。洙漕河流域的茶园面积由 2005 年的 3.13 km2增长至 2019 年的 8.03 km2，沙河西岸的茶园面积由2005 年的2.84 km2增长至2019年的4.59 km2，平桥河流域的茶园面积由2005年的1.85 km2增长至2019年的4.39 km2。这三大茶园的聚集区均为居民区，道路较为集中，可达性较高，可见茶园在天目湖流域空间分布以及变化上的差异受到了社会经济以及基础设施完备性的极大影响。

图2 2000—2019年天目湖流域茶园空间分布图Figure 2 Spatial distribution of tea plantations in Lake Tianmu basin

1.3 SWAT-APEX多尺度集成模拟

1.3.1 模型选择与集成方法

本研究采取SWAT 模型和APEX 模型集成模拟的方法，即：在SWAT 模型划分的子流域中选择茶园分布的子流域，利用APEX 模型进行空间子区细分，将子流域模拟结果导入SWAT 模型，进行大流域的河道演算。通过集成方式分别发挥两个模型的优势，完成天目湖流域以茶园为主要土地利用类型的非点源污染产出对水环境影响的评价目标。

在本研究构建的SWAT 模型中，根据地形、土地利用类型、土壤类型和河道位置，将天目湖流域划分为359 个子流域。为实现对流域内各茶园地块更高的空间分辨率定量模拟，在359个子流域中筛选出82个茶园聚集的子流域，针对这些子流域进行子流域-子区的空间离散化，作为APEX模型的基本模拟单元。为保证子区划分的科学性与合理性，子流域-子区的空间离散化对地形、土地利用以及土壤类型进行了综合考虑，共划分出139 个茶园子区。最终的天目湖流域SWAT-APEX模型空间离散化结果如图3所示。

图3 天目湖流域空间离散化结果Figure 3 Subbasin and subarea of the Lake Tianmu basin

1.3.2 模型率定与验证

针对SWAT-APEX 集成模型的模拟需要分别完成SWAT 模型和APEX 模型相关参数的更新和校正。SWAT模型和APEX模型提供的土地覆盖类型参数数据库中不包含茶园类型的参数，本研究基于SWAT 模型已有作物参数库，根据文献在作物参数库中添加了茶园类型的参数[26]。在此基础上，首先根据参数的敏感性完成SWAT 模型参数的校正，参数校正结果见文献[27]；再针对APEX 模型控制数据库、子区数据库、参数数据库中相关参数，如对径流量模拟影响较大的SCS 曲线系数、对氮磷产出模拟影响较大的氮磷富集率系数等进行率定[20，28-29]。2010—2012年作为模型预热期，利用鲢鱼桥站2013—2018年月均径流观测数据对模型径流模拟效果进行评价，利用鲢鱼桥站、洙漕坝桥站、潘村桥站和徐家园桥站2013—2017 年月均总氮、总磷监测数据对模型水质模拟效果进行评价。精度评价指标选取决定系数（R2）和纳什系数（Ens）[28，30]，评价函数的计算公式如下：

式中：Oi和Om分别代表实测值及其平均值；Si和Sm分别表示模拟值及其平均值；n为实测或模拟数据的个数。

2 结果与讨论

2.1 SWAT-APEX模型的率定与验证结果

鲢鱼桥站的月均径流模拟和实测结果对比如图4所示。2013—2015年为模型径流模拟校正期，2016—2018 年为模型径流模拟验证期，校正期和验证期的R2和Ens均大于0.6，表明SWAT-APEX 集成模型得到了较好的径流模拟效果。

图4 鲢鱼桥站月均径流量模拟值与实测值对比Figure 4 Calibration and validation results of monthly average flow from 2013 to 2018 at the Lianyuqiao station

洙漕坝桥站、潘村桥站、鲢鱼桥站和徐家园桥站的月均总氮、总磷模拟和实测结果对比以及精度评价如图5 和图6 所示。洙漕坝桥站、潘村桥站和徐家园桥站3 个站点在校正期（2013—2015 年）和验证期（2016—2017 年）的总氮和总磷模拟结果R2都大于0.5、Ens都大于 0.4，且其中大部分R2大于 0.6、Ens大于0.5，模拟精度满足要求[26-27]。鲢鱼桥站在2013—2015年的模拟精度偏低，但是在2016—2017年的模拟精度较高。原因是当地在2013—2015 年期间对中田舍河实施了大量河道清淤和整治工程，位于中田舍河的鲢鱼桥断面该阶段的水质监测数据受到强烈的人为影响，超出了模型的模拟能力。而在后期不受河道工程影响时，模型则较好地模拟了河流水质过程。总体来看，SWAT-APEX集成模型获得了较好的氮、磷模拟精度，其中总磷模拟精度整体优于总氮模拟精度。

图5 月均总氮模拟值与实测值对比Figure 5 Calibration and validation results of monthly TN from 2013 to 2017

图6 月均总磷模拟值与实测值对比Figure 6 Calibration and validation results of monthly TP from 2013 to 2017

2.2 茶园规模开发前后氮磷产出变化模拟与分析

基于天目湖流域SWAT-APEX 集成模型参数率定结果，在SWAT 模型中分别输入2000、2005、2010、2015、2019 年 5 期土地利用类型数据，在 APEX 模型中分别设定不同的情景，对天目湖流域2000—2019年的总氮、总磷产出分别进行流域整体、子流域、茶园子区的多尺度分析。

在整个流域层面，根据模型模拟结果对2000、2005、2010、2015 年和 2019 年不同用地类型总氮和总磷的贡献率进行计算，如图7所示，从2000—2019年，茶园用地类型对天目湖流域的氮磷负荷贡献持续上升。至2019年，茶园的年均总氮和总磷负荷贡献率分别达到了30%和20%，是天目湖流域年均总氮负荷的第一贡献者、年均总磷负荷的第二贡献者。

图7 天目湖流域不同用地类型的总氮与总磷产出占比Figure 7 Percentage of TN and TP contribution from different land use types in Lake Tianmu basin

在子流域层面，利用2019 年与2000 年天目湖流域两期土地利用模拟的子流域多年平均总氮、总磷产量相减，获得了2000—2019 年总氮、总磷产量变化的空间分布图，如图8 所示。年均总氮、总磷产量增加的子流域面积分别占55%和77%，子流域总氮产出最大增加了3 891.7 kg，子流域总磷产出最大增加了326 kg。经分析，2000—2019 年总氮、总磷产出明显增加的子流域与其内部茶园用地类型面积的增加具有正相关关系。

图8 天目湖流域2000—2019年总氮与总磷产量变化空间分布Figure 8 Spatial distribution of changes in TN and TP production from 2000 to 2019 in Lake Tianmu basin

在茶园子区层面，结合SWAT-APEX 集成模型的模拟结果制作出2019年茶园子区和2000年对应范围子流域的年均总氮与总磷产出空间分布图，如图9 所示。2019年天目湖流域的茶园子区在2000年时的土地利用类型主要为耕地和林地。根据土地利用变化分析，近20 a中，近56%的茶园是由林地转换而来，其余茶园则是由耕地转换而来。在2000 年，茶园子区对应范围子流域的年均总氮产出为0.87~710.63 kg，年均总磷产出为0.19~151.62 kg。至2019 年，茶园子区的年均总氮产出增加到3.10~3 982.81 kg，年均总磷产出增加到2.98~316.07 kg。从2000 年到2019 年，经过茶园的规模化开发，茶园子区的年均总氮与总磷产量均得到大幅增加。

图9 2000年与2019年各茶园子区年均总氮与总磷产出空间分布Figure 9 Spatial distribution of TN and TP production from tea plantation subareas in 2000 and 2019

通过茶园规模化开发后的2019 年茶园子区总氮、总磷产出空间分布图可以看出，茶园子区年均氮磷产出较为严重的区域主要聚集在洙漕河流域和平桥河流域，这与茶园集中种植区相吻合。其中，洙漕河流域既属于大溪水库西部的水源涵养区，也属于十里长山-南山山水生态保护带，平桥河流域主要属于沙河水库东部的水源涵养区。2000 年部分位于大溪水库西部水源涵养区的洙漕河流域以及沙河水库东部水源涵养区平桥河流域的水源涵养林，到2019 年时已被聚集成片的茶园所替代。茶园植株高度远低于涵养林，且在采摘茶叶时需经常对其除草，这必然导致两大重要水源涵养区的入湖污染拦截能力大幅降低。一旦该区域发生降雨，将产生严重的氮磷流失，继而引发水库富营养化问题，此种风险在茶园开发初期时尤为显著。因此，就茶园子区现状的污染程度而言，有必要对水源涵养区的茶园采取优化种植管理、设置拦污设施，甚至退茶还林等措施。

2.3 针对茶园管理措施的污染削减效果评估

根据天目湖流域不同用地类型对养分负荷的贡献，明确了茶园是自然用地中造成氮磷污染的最重要用地类型。参考《天目湖保护规划（2018—2035）》草案，分别从农业生产措施管理（施肥削减）、工程管理（过滤带设置）和土地资源管理（退茶还林）3 方面因素设计3 种情景来研究不同管理措施下茶园污染产出的削减效果。情景1 设计为将茶园化肥施用量削减至原来的80%；情景2 设计为在茶园的集中种植区设置10 m的植被过滤带；情景3设计为将所有茶园转换为林地。各情景下茶园子区的氮磷削减效果如表3所示。

表3 3种茶园管理措施下平均总氮、总磷削减率Table 3 Average TN and TP reduction rate of three scenarios of tea plantation managements

情景1 中，各茶园子区的年均总氮、总磷产出削减率分别为6.23%~24.31%、1.75%~13.22%，平均削减率分别为14.47%和6.77%。在化肥施用量适当削减的情况下，茶园子区的年均氮磷产出削减效果并不显著。茶园的土壤氮素呈现出输入量高、利用率低、流失量大的特点[31]，在茶园施肥过程中也常采用高氮含量的肥料。尽管削减施肥量对茶园年均总氮、总磷产出的削减效果有限，但仍可通过适当减少施肥量、避免雨前施肥或者施用稀释肥料等措施来减少氮磷流失。

情景2 中，各茶园子区的年均总氮、总磷产出削减率分别为25.18%~53.85%、31.58%~84.69%，平均削减率分别为29.46%和42.37%，年均总磷产量的削减效果优于总氮。整体来说，设置10 m 植被过滤带的效果要优于控制化肥施用量的效果，然而也需要较高成本。

情景3 中，各茶园子区的年均总氮、总磷产出削减效率分别为14.57%~98.68%、17%~96.76%，平均削减效率分别为71.97%和74.42%。将茶园转换为林地是保护流域水环境最直接、最高效的方法，然而，这种方式将大幅降低丘陵山区农民的经济收益，与地区经济发展相矛盾。根据前文分析，天目湖流域的水源涵养区聚集了大量茶园，是氮磷流失的高风险区。因此，可适当将位于水源涵养区的茶园转换为涵养林，加强水源涵养区的污染物入湖拦截能力。

综合看来，3 种情景中，退茶还林的污染削减效果最优，设置10 m 植被过滤带的污染削减效果次之，控制施肥量的污染削减效果较为有限。然而考虑到现实情况与经济成本，可将3 种方法适当结合，如将水源涵养区中造成高污染负荷的茶园区用林地替代，在其余茶树种植聚集区设置5~10 m 的植被过滤带，同时避免雨前施肥。

3 结论

（1）2000 年，茶园仅零星分布在洙漕河流域。2005—2019 年规模化开发茶园后，沙河水库西岸、洙漕河流域、平桥河流域的茶园迅速扩张，成为天目湖流域主要的茶园聚集区。茶园开发模式呈现出水库水源涵养区集中成片、水库周边向上游持续蔓延的特征。

（2）经过APEX-SWAT 集成模型的参数率定与结果验证，模型月尺度径流产出的确定性系数和纳什效率系数在校正期达到0.6 以上，验证期达到0.8 以上；月尺度的总氮和总磷产出的确定性系数在率定期与验证期基本达到0.6以上，纳什效率系数基本达到0.5以上。SWAT-APEX 集成模型在天目湖流域具有较好的适用性。

（3）通过模型模拟对比了茶园规模化开发前后流域总氮、总磷的时空产出特征，即：茶园用地类型对天目湖流域的氮磷负荷贡献持续上升，至2019 年，茶园总氮产出贡献超出其他土地利用类型成为第一贡献者，总磷产出仅次于耕地成为第二贡献者，洙漕河流域和平桥河流域成为两大重点氮磷产出的茶园聚集区。

（4）设定施肥削减（将茶园化肥施用量削减至原来的80%）、过滤带设置（在茶园的集中种植区设置10 m 植被过滤带）和退茶还林（将所有茶园转换为林地）3 种不同情景对茶园氮磷产出的削减效果进行评估，3 种情景的平均总氮削减率分别为14.47%、29.46%、71.97%，平均总磷削减率分别为6.77%、42.37%、74.42%。建议根据当地实际情况与经济效益，将水源涵养区高氮磷流失的茶园区转换为林地，结合过滤带设置与避免雨前施肥等措施，实现丘陵山地茶园开发影响下的水源地水生态保护调控。