气候变化下乌伦古河流域农业面源污染负荷响应

邹凯波， 张玉虎， 刘晓伟， 薛淑慧， 杨博文， 崔艳欣

（1.首都师范大学资源环境与旅游学院，北京 100048；2.生态环境部华南环境科学研究所，广东 广州 510530）

近年来，随着人类社会经济的高速发展，由此引发的全球气候变暖及其对生态环境产生的负面影响，越来越受到社会各界的重视，成为学术界关注的热点问题。在我国西北内陆干旱区，水资源是制约社会经济发展，影响生态安全的关键因素，对未来社会经济可持续发展起着至关重要的作用。农业面源污染由于其污染物的广域性、分散性、污染物质迁移途径的无序性，具有机理模糊、潜伏周期长、危害性大等特点［1］，相对而言，点源污染易于识别和治理，随着工业和城市生活点源污染控制措施逐步取得成效，目前面源污染成为水环境污染重要来源［2］。而我国西北地区是农业面源污染增长速率和排放强度较高的主要区域［3-4］，排污来源主要包括畜禽养殖和化肥养分流失。2018 年新疆单位播种面积的化肥总施用量为5.42 kg·hm-2，是安全化肥施用标准的1.7倍［5］。

我国西北干旱区地处中纬度地带的欧亚大陆腹地，是对全球气候变化响应最为敏感的地区之一［6］。近些年来，西北部分地区暖干向暖湿转变趋势明显。李哲等［7］研究表明1960—2019 年西北地区年均气温整体呈显著上升趋势，增温速率为0.26 ℃·（10a）-1。施雅风等［8］利用西北地区气温、降水、径流等大量资料，研究分析西北地区温度、降水变化趋势，认为该地气候由暖干向暖湿转型。李江等［9］对1956—2016 年新疆气象资料分析表明气温增长速率约0.3 ℃·（10a）-1，降水增长速率约10 mm·（10a）-1。而乌伦古河流域由于周边土地利用资源过度开发和农牧业迅速发展使得流域面源污染逐渐严重［10-11］。因此，探讨气候变化对流域农业面源污染影响，对于应对气候变化情景下流域农业面源污染防控与治理有重要实践意义。

SWAT 是美国农业部农业研究局（USDA-ARS）历经30 多年开发出来的一套长时段流域分布式水文模型，其具有很强的物理基础，适用于具有不同的土壤类型、土地利用方式和管理条件下的复杂大流域，并能在资料缺乏的地区建模，目前在国内外农业面源污染研究中得到了广泛的应用［12-17］。Wang 等［18］利用SWAT 评价未来气候变化对辽河源区融雪期面源污染的影响。国内已有不少学者将SWAT 模型应用到气候变化对面源污染的影响，刘吉开等［19］采用SWAT模型对渭河流域陕西段面源污染负荷影响进行研究。赵雪松等［20］应用SWAT模型定量分析气候变化对营口地区农业面源污染的影响，研究表明气候变化对区域农业面源污染有显著影响。

综合以上SWAT 模型的应用研究，主要集中于气候变化对面源污染负荷的总量变化方面，系统分析径流和面源污染对气候变化响应的研究不多，且将SWAT模型应用到探究气候变化下西北地区流域面源污染响应的文献尚不多见。因此，本研究以新疆乌伦古河流域为研究区域，构建乌伦古河流域径流及面源污染SWAT 模型，重点模拟分析不同气候变化情景对流域径流、总氮、总磷变化影响及氮磷空间分布特征。以期为乌伦古河流域在应对气候变化下水资源的保护与利用、土地利用规划、农业生产活动方面提供科学依据。

1 研究区概况

乌伦古河流域位于新疆阿勒泰地区，地理位置位于46°10′～47°28′N，87°05′～90°45′E，由发源于阿尔泰山东段南坡的大青河、小青河、查干郭勒河、强罕河和发源于阿尔泰山北坡蒙古国境内的布尔根河等组成。自东向西流经青河县、富蕴县、福海县，汇入新疆第二大湖——乌伦古湖。乌伦古河流域总面积37882 km2，其中，国外面积为10310 km2，河流全长811 km，主要依靠山地降水和夏季融雪补给，用水以农牧业灌溉为主。流域年均气温2.3 ℃，年均降水量129.8 mm，年均蒸发量867.2 mm，无霜期110～130 d，多年平均日照时数3157.5 h。乌伦古河流域包含青河、富蕴、福海3 个县级行政单位，涉及17个乡镇、2个农场和1个兵团，流域地理位置如图1所示。

图1 乌伦古河流域位置示意图Fig.1 Location of Ulungur River basin

2 模型构建与研究方法

2.1 SWAT模型构建与校准

2.1.1 SWAT 模型构建 SWAT 模型是美国农业部开发的适用于较大流域尺度的面源污染计算模型［21］，不仅能够模拟长时间序列水文过程，而且在资料缺乏地区同样能建立水文过程模拟［22］。因此，被广泛应用在非点源污染的管理与控制过程中。构建SWAT模型的基本数据分为空间数据和属性数据，空间数据包括数字高程模型（DEM）、土地利用图、土壤类型图等，属性数据包括土壤属性、气象、水文水质等数据［23］。

在SWAT 模型数据构建过程中，需要将所有数据图层的坐标系和投影类型保持一致，按照SWAT模型的一般要求，采用Albers 等圆锥投影进行图层和数据的投影变换［24］。本文根据研究区所处地理位置，投影均选择WGS_1984_UTM_Zone_45N，基准面D_WGS_1984，数据格式为栅格型。构建模型采用的数据如表1所示。

表1 模型构建采用的数据Tab.1 Data required for model construction

（1）数字高程模型（DEM）数据

数字高程模型（DEM）数据来自于地理空间数据云平台所提供的30 m 分辨率数据，并利用Arc-GIS 10.2 对其进行投影转换及一系列预处理操作，生成SWAT模型所需要的DEM，如图2a所示。

（2）土地利用类型分布

土地利用数据来源于清华大学地球系统科学系宫鹏教授科研团队发布的1:250 000 土地利用栅格数据，参照SWAT 自带的各种土地利用类型参数库，最终将其重分为7类土地利用类型，重分类结果及空间分布如表2和图2b所示。

图2 乌伦古河流域DEM和土地利用Fig.2 DEM and land use of Ulungur River basin

表2 土地利用SWAT重分类结果Tab.2 Land use SWAT reclassification results

（3）土壤类型分布

土壤类型数据来源于联合国粮食及农业组织（FAO）创立的分辨率为1:1000 000的世界土壤数据库（HWSD）。中国境内数据是在第二次全国土地调查中由南京土壤所提供的1:1000 000 的土壤数据。根据SWAT 模型输入要求，将土壤类型进行重分类（表3），其空间分布如图3所示。

图3 乌伦古河流域土壤类型Fig.3 Soil type map of Ulungur River basin

表3 乌伦古河流域土壤类型Tab.3 Soil types in the Ulungur River basin

（4）土壤属性数据库

构建SWAT模型中的土壤属性数据库主要包括以下相关参数，参数名称和含义如表4所示，土壤的物理属性直接决定了土壤剖面中水和气的运动情况，对SWAT 模型中水文响应单元（HRU）的水循环有重要作用。大部分参数可以直接在HWSD 土壤数据集查询到，少数需要借助土壤特性计算工具SPAW（Soil Plant Atmosphere Water）进行计算，结合HWSD 数据库中查询获取的有机物质（Organic Matter）和不同粒径土壤值，则可计算获取以下相关参数的值。

表4 土壤属性数据库相关参数Tab.4 Soil attribute database related parameters

（5）气象数据与天气发生器

SWAT模型所需的实测气象数据包括气象站点的位置、高程和气象站逐日观测数据，气象数据包括日降水量、日气温（最高气温、最低气温）、日相对湿度、日太阳辐射、日平均风速。由于研究区内气象站点较少，本研究采用来自中国水利水电科学研究院发布的SWAT 模型中国大气同化数据集CMADS V1.0，数据集空间分辨率为1/3°×1/3°，时间尺度为2008—2018 年的逐日数据。乌伦古河流域CMADS数据集气象站点分布如图4所示。

图4 乌伦古河流域内CMADS数据集气象站点分布Fig.4 Distribution map of meteorological stations in the CMADS dataset in the Ulungur River basin

（6）划分子流域与水文响应单元

基于DEM、土地利用类型、土壤类型数据将乌伦古河流域划分为23 个子流域（图5），分别设定土地利用、坡度分布和土壤类型的最小面积比率阈值为10%、15%、10%，如果低于该阈值，则子流域内所属类型忽略不计，最终SWAT 模型在该研究区域内一共生成了276 个水文响应单元（HRUs）。经计算提取流域面积为19572 km2。

图5 乌伦古河流域子流域划分Fig.5 Sub-basin division map of the Ulungur River basin

2.1.2 模型率定与验证 采用决定系数R2（Coefficient of Determination）和 纳 什 系 数Ens（Nash-Suteliffe）这两个指标来评价模型的适用性。一般认为，R2＞0.75 且Ens＞0.75 时，模型的适用性非常好；R2≤0.50、且Ens≤0.50时，模型模拟结果不可信；当结果值介于上述范围之间时，认为模拟达到了可用性，模型模拟结果是可以接受的。

根据输入气象数据时间范围和实测水文水质数据时长，本研究模拟时间段确定为2008—2018年，选取流域中段二台水文水质站点数据进行月尺度的径流、总氮、总磷的模拟。其中2008—2009 年为模型运行的预热期，利用2010—2018年径流实测数据对径流参数进行率定验证，利用2017—2018年总氮、总磷实测数据对水质参数进行率定验证。由率定期及验证期月径流量和营养物质模拟效果评价指标可知（表5），径流和营养物质在率定期与验证期的R2、Ens均大于0.75，满足了模型的精度要求。因此，该模型能够用于乌伦古河流域农业面源污染负荷模拟。

表5 二台站月径流及营养物质模拟评价指标Tab.5 Simulated evaluation index of monthly runoff and nutrient substance in Ertai hydrological station

2.1.3 主要参数 依次完成流域径流、总氮和总磷负荷率定验证后，SWAT 模型相关参数经过一系列调整，最后确定的参数值如表6所示。

表6 SWAT模型参数取值Tab.6 SWAT model parameter value

2.2 研究方法

2.2.1 气候变化特征分析 为了探明乌伦古河流域的气候变化特征，对1962—2018年乌伦古河流域福海、阿勒泰、富蕴3 个气象站的气温、降水数据进行分析（图6）。乌伦古河流域近57 a气候变化特征分析结果表明，该流域气温增长率范围为0.191～0.653 ℃·（10a）-1，降水增长率范围为7.290～16.531 mm·（10a）-1，说明乌伦古河流域气候已向暖湿的趋势发展。

图6 1962—2018年乌伦古河流域气象站点年均气温和降水量变化趋势Fig.6 The trend of annual average temperature and precipitation at meteorological stations in the Ulungur River basin from 1962 to 2018

2.2.2 气候变化情景设置 以气候变化特征分析结果为依据，采用任意情景设置的方法［25］分别设置3种降水和3 种气温变化互相组合情景，在乌伦古河流域2008—2018 年逐日气温、降水基础上进行改变，将温度分别增加0 ℃、1.0 ℃、2.0 ℃，降水量分别增加0%、10%、20%，将变化的气温与降水情景组合构成9种气候变化情景（表7），带入构建完成的乌伦古河流域SWAT 模型进行模拟，设置2008—2009年为预热期，对2010—2018年的径流和水质模拟结果进行分析。

表7 气候变化情景组合Tab.7 Climate change scenario portfolio

3 模拟结果及分析

3.1 气候变化情景影响

3.1.1 径流对气候情景变化的响应 通过将9 种气候变化情景分别带入SWAT 模型，模拟得到乌伦古河流域9 种气候变化情景下的模拟结果，将2010—2018年的模拟结果进行整理，计算得到不同气候变化情景下的年平均径流量及变化百分比（表8）。从多年平均径流量结果可知，降水量变化对径流量影响较大。在气温保持不变的情况下，增加10%的降水量，径流量增加23.07%，增加20%的降水量，径流量增加48.01%。这主要因为径流产生是由于强降水，增加全年整体降水量就从整体上增加了强降水次数，从而增加径流量［25］。在降水量保持不变的情况下，气温上升1 ℃，径流量减少3.22%，气温上升2 ℃，径流量减少6.50%。气温上升会导致冰雪融化，使得流域上游径流量有所增加。但同时由于气温的上升也导致区域总的蒸发量增加，使得径流量有所减少。气温主要通过冰雪融化、流域水面蒸发等因素间接影响径流。而降水是流域径流形成的基础，是流域水资源得到补充最重要的一个方面，降雨量会直接影响该流域的水资源变化和水文过程。通过对比可以发现，气温变化对径流量的影响要明显小于降水量的变化对径流量的影响。

表8 不同情景下流域多年平均径流量及变化百分比Tab.8 Annual average runoff and change percentage of the basin under different scenarios

3.1.2 水质对气候情景变化的响应 水质对气候情景变化的响应主要是通过气候变化使得流域径流量发生改变，而径流又是面源污染的主要驱动力。依次对9 种气候情景下2010—2018 年的总氮和总磷输出负荷进行计算，可以得到每种气候变化情景下总氮、总磷的年平均负荷量。每种气候情景与S11基准情景（即气温、降水无变化时，年平均总氮负荷为10.17×103t、总磷负荷为1.88×103t）进行对比，可以得到不同情景下年平均总氮负荷、总磷负荷变化百分比（表9）。

表9 不同情景下流域多年平均总氮、总磷负荷量及变化百分比Tab.9 Multi-year average total nitrogen and total phosphorus loadings and percentage changes in the basin under different scenarios

从总氮、总磷负荷变化百分比来看，其负荷变化主要受降水影响较大，并且变化趋势与径流相似，呈现较为明显的变化特征。当气温保持不变时，降水量增加10%，总氮负荷增加11.25%、总磷负荷增加14.32%；降水量增加20%，总氮负荷增加23.19%、总磷负荷增加29.65%。当降水量保持不变时，气温升高1 ℃，总氮负荷减少1.86%、总磷负荷减少3.21%；气温升高2 ℃，总氮负荷减少5.66%、总磷负荷减少7.11%。

3.2 污染负荷空间分布情况

3.2.1 总氮空间分布特征 9 种不同气候变化情景下，2010—2018年各子流域年均单位面积总氮负荷输出情况如图7 所示，为方便比较单位面积污染负荷变化情况，图例采用同一限值范围。首先对情景S11下的单位面积总氮负荷输出空间分布特征进行分析可知，1 号子流域单位面积总氮负荷输出量最大，为60.87～164.79 kg·hm-2·a-1；其次是子流域15、19，单位面积总氮负荷为26.83～60.87 kg·hm-2·a-1；子流域4、5、7、8、9、12、14、17、18、21、23 次之，为1.67～26.83 kg·hm-2·a-1。子流域2、3、6、10、11、13、16、20、22 单位面积总氮负荷最小，为0.01～1.67 kg·hm-2·a-1。

图7 不同气候变化情景下年均单位面积总氮负荷输出空间分布Fig.7 Spatial distribution of annual average total nitrogen load output per unit area under different climate change scenarios

当由情景S11变为情景S12时，即气温保持不变，降水量增加10%，子流域8、15、20、21单位面积总氮负荷升高；当由情景S12变为S13时，即气温保持不变，降水量增加20%，可以发现，子流域4、19 的单位面积总氮负荷在情景S12的基础上有所上升。

当由情景S11变为S21时，即降水量保持不变，气温升高1 ℃，可以发现子流域20单位面积总氮负荷升高；由情景S21变为S22时，即气温保持不变，径流量增加10%，子流域8、15、21、23单位面积总氮负荷增加；当由情景S22变为S23时，即气温保持不变，降水量增加20%，可以发现子流域19单位面积总氮负荷在情景S22的基础上有所增加。

当由情景S11变成S31时，即降水量保持不变，气温升高2 ℃，可以发现仅有子流域20单位面积总氮负荷增加；当由情景S31变为S32时，即气温保持不变，径流量增加10%，子流域8、15、21、23单位面积总氮负荷上升；当由情景S32变为S33时，即气温保持不变，降水量增加20%，仅有子流域19单位面积总氮负荷范 围 由26.83～60.87 kg·hm-2·a-1增 加 为60.87～164.79 kg·hm-2·a-1。

综上所述，1号子流域单位面积总氮负荷最大，当气温保持不变时，子流域8、15、19、20、21、23会随着降水量的增加而不同程度的增加，且子流域19只有在降水量增加20%时才会出现显著上升；当降水量不变时，随着气温的升高，子流域20 单位面积总氮负荷增加。子流域2、3、6、10、11、13、16、22 受气候变化影响时单位面积总氮负荷未发生明显改变，通过叠加空间分布图层到乌伦古河流域遥感影像图中发现以上子流域范围内荒漠化程度较高，由此推断在荒漠化地区污染负荷受气候变化影响较小。

3.2.2 总磷空间分布特征 9 种不同气候变化情景下，2010—2018年各子流域年均单位面积总磷负荷输出情况如图8所示，同总氮一样，对情景S11的单位面积总磷负荷输出空间分布特征进行分析可知，1号子流域单位面积总磷负荷输出量最大，为10.59～45.10 kg·hm-2·a-1；其次是子流域15、19，单位面积总磷负荷为4.53～10.59 kg·hm-2·a-1；子流域4、5、7、8、9、12、14、17、18、21、23 次之，为0.36～4.53 kg·hm-2·a-1。子流域2、3、6、10、11、13、16、20、22单位面积总磷负荷最小，为0.01～0.36 kg·hm-2·a-1。

图8 不同气候变化情景下年均单位面积总磷负荷输出空间分布Fig.8 Spatial distribution of annual average total phosphorus load output per unit area under different climate change scenarios

当由情景S11变为情景S12时，即气温保持不变，降水量增加10%，子流域8、12、14、15、17、20单位面积总磷负荷升高；当由情景S12变为S13时，即气温保持不变，降水量增加20%，可以发现，子流域19单位面积总磷负荷在情景S12的基础上有所上升。

当由情景S11变为S21时，即降水量保持不变，气温升高1 ℃，可以发现子流域20单位面积总磷负荷升高；由情景S21变为S22时，即气温保持不变，径流量增加10%，子流域8、12、14、15、17 单位面积总磷负荷增加；当由情景S22变为S23时，即气温保持不变，径流量增加20%，可以发现子流域19单位面积总磷负荷在情景S22的基础上有所增加，区间范围由4.53～10.59 kg·hm-2·a-1增加为10.59～45.10 kg·hm-2·a-1。

当由情景S11变成S31时，即降水量保持不变，气温升高2 ℃，可以发现仅有子流域20单位面积总磷负荷增加；当由情景S31变为S32时，即气温保持不变，降水量增加10%，子流域8、12、14、15单位面积总磷负荷上升；当由情景S32变为S33时，即气温保持不变，降水量增加20%，仅有子流域19单位面积总磷负荷升高。

综上所述，1号子流域单位面积总磷负荷最大，当气温保持不变时，子流域8、12、14、15、17、20会随着降水量的增加出现不同程度的增加，且子流域19只有当降水量增加20%时才出现显著上升；当降水量不变时，随着气温的升高，子流域20 单位面积总磷负荷增加，推断可能随着温度的升高流域上游部分地区的融雪过程可能会导致土壤侵蚀的加剧，进而导致总磷负荷增加。与总氮一样，在子流域2、3、6、10、11、13、16、22 范围内单位面积总磷负荷也未因气候变化改变而出现明显改变。从子流域单位面积总磷负荷分布特征来看，可以发现与总氮基本保持一致，但单位面积总磷负荷强度要小于单位面积总氮负荷强度。

3.3 气候变化对面源污染影响

不同气候变化情景对流域径流及污染负荷都有不同程度影响。面源污染主要是由于降水带来的径流引起，随着径流量的增加，总氮、总磷年均负荷量增加较为显著，这表明径流量是乌伦古河流域农业面源污染物流失的主要载体。而随着温度的增加，子流域20总氮、总磷负荷量随之增加，表明流域上游部分地区的融雪过程可能会导致土壤侵蚀的加剧，进而导致污染负荷增加。

根据总氮和总磷单位面积污染负荷空间特征分析结果，考虑到以综合防治为出发点。将子流域1、15、19 作为流域农业面源污染重点防治区域；子流域8、12、14、17、20、21、23作为气候变化下农业面源污染重点影响区域（图9），可以看出这些区域主要分布在福海县的阿尔达乡、富蕴县的吐尔洪乡、库尔特乡、恰库尔图镇、克孜勒希力克乡，及青河县的阿热勒乡、阿热勒托别镇、塔克什肯镇、萨尔托海乡。结合实地调研结果分析，阿尔达乡所处范围内畜禽养殖场分布较多，且阿尔达乡处于流域末端，汇聚了部分来自上游的污染负荷，使得其污染负荷较为严重。针对这一现状，应加快推进建设堆粪场，减少畜禽养殖带来的粪污污染，使得大型规模养殖场粪污处理设施配套率达标。气候变化主要通过降水、温度对流域径流量带来影响，从而对流域面源污染负荷产生影响。因此，应合理规划放牧时间，在降水时期，应限制大规模放牧活动。同时为削减由于降水导致的农田径流作用，应对河流沿岸实施绿化带修复工程，能极大程度降低氮磷等营养物质进入水体。

图9 流域农业面源污染区县图Fig.9 Map of the area and counties of the watershed agricultural non-point source pollution

4 结论

（1）构建了适用于乌伦古河流域的面源污染SWAT模型。利用径流和水质数据对模型相关参数进行率定与验证。径流、总氮和总磷在率定期和验证期的决定系数R2均在0.75以上，纳什系数Ens均在0.55 以上。模拟结果精度均达到模型模拟的要求，说明建立的SWAT模型能够用于乌伦古河流域农业面源污染负荷模拟。

（2）气温变化对径流量的影响小于降水量变化对径流量的影响。总氮负荷及总磷负荷随着降水量的增加而增加，且总磷负荷的变化幅度要高于总氮负荷的变化幅度。增加10%的降水量，径流量、总氮负荷、总磷负荷分别增加23.07%、11.25%、14.32%。增加20%的降水量，径流量、总氮负荷、总磷负荷分别增加48.01%、23.19%、29.65%。

（3）乌伦古河流域各子流域年均单位面积总氮负荷为0.01～164.79 kg·hm-2·a-1，单位面积总磷负荷为0.01～45.10 kg·hm-2·a-1。1 号子流域单位面积总氮及总磷负荷最大。气候变化对面源污染影响较大的区域主要分布在福海县的阿尔达乡，富蕴县的吐尔洪乡、库尔特乡、恰库尔图镇、克孜勒希力克乡，青河县的阿热勒乡、阿热勒托别镇、塔克什肯镇、萨尔托海乡。这些乡镇是面源污染防治的关键区域，乌伦古河流域未来面源污染防治应重点控制畜禽养殖和农田径流这两类污染源。