李 乐, 马 巍, 勾蒙蒙, 王 娜, 刘常富,2, 肖文发,2
(1.中国林业科学研究院森林生态环境与自然保护研究所,国家林业和草原局森林生态环境重点实验室,北京 100091;2.南京林业大学南方现代林业协同创新中心,南京 210037;3.中国水利水电科学研究院水生态环境研究所,北京 100038)
农业面源输出的氮素是造成地表水富营养化的主要原因。硝态氮(NO—N)通常是进入水体氮的主要形式,这是由于NO—N具有高度的溶解性和稳定性。有研究指出,气候变化、土地利用/覆被变化及集约化农业管理之间的协同作用放大氮素流失的负效应,并且破坏当前的缓解效果。为了制定有效缓解策略,不仅要了解NO—N从径流输送到地表水的过程,还要系统了解NO—N的来源和驱动因素。
目前,基于观测试验和模型模拟开展了大量土地利用类型对面源污染影响的研究,但养分传输还受气候、水文特征、下垫面属性等多种因素及其相互作用的影响。然而,以往的研究仅考虑土地利用类型或很少结合其他因素,这可能导致研究结果出现偏差。径流是NO—N向受纳水体运移的主要途径,但径流组成(地表径流和基流)在不同流域间具有显著差异,并随降雨、生产活动等因素的变化发生季节性变化。由于基流难以直接观测和量化,目前的研究多关注地表径流或总径流,很少关注基流的贡献。将不同土地利用类型的NO—N负荷与水文过程相结合,量化不同水文过程的贡献有助于确定NO—N来源和实现NO—N污染的有效控制。另外,土壤条件决定径流分布,影响土壤水和溶质中污染物的运移。例如,土壤渗滤性是影响养分淋溶损失的主要因素之一;地形特征也在一定程度上决定污染物由扩散源向河流的输送路径。气候、地形、土壤和其他自然属性往往决定某一地区的土地利用活动,且这些因素以复杂的方式相互作用,但综合考虑这些因子对面源污染的影响仍未得到充分认识。流域自然地理特征对面源污染的影响通常被误解为是土地利用变化所造成影响的一部分。Varanka等研究发现,地貌(土壤和基岩属性)对芬兰北方河流水质的独立贡献高于土地利用,这一发现强调地貌等自然属性在预测水质空间变化方面的重要作用。若在分析面源污染的影响因素过程中单独考虑土地利用,则可能低估或高估受纳水体实际受到的外界影响程度。此外,土地利用和面源污染的关系往往是区域特有的。因此,在研究土地利用对面源污染的影响时,有必要全面考虑自然和人为因素并确定各因素的贡献。
通过大量控制试验分离影响面源污染的多种因素具挑战性,这不仅因为众多因素之间的相互作用,也源于水文过程的高度可变性以及面源污染的间歇性、随机性和难于监测的特性。基于物理机制的水文模型将模型参数直接与物理可观测的地表特征联系起来,为面源污染的定量研究提供了有效途径。在这些机理模型中,SWAT(soil and water assessment tool)模型较全面考虑氮循环过程(矿化、吸收、挥发、淋溶等过程),适用于不同下垫面特征的复杂流域(包括喀斯特流域),在数据不充足流域也具有较高模拟精度。
三峡库区位于保障长江流域生态安全的关键地带,是我国乃至东亚典型生态敏感区和重要的淡水水源地。面源污染(尤其是库区腹地)是威胁水库及长江中下游生态安全的主要问题,但是库区腹地在土壤、地形等自然和社会经济背景方面存在明显空间差异。因此,选择具有不同地理背景的流域揭示NO—N来源及驱动因素的一致性和差异性,具有典型性和代表性。综上,本文以三峡库区腹地的大宁河流域(喀斯特流域)和梅溪河流域(喀斯特—非喀斯特组合流域)为研究区,利用SWAT模型模拟水文和NO—N流失过程,解析流域NO—N输出的季节变化特征和来源贡献,通过随机森林回归评估影响径流和NO—N负荷因素的影响程度,以期为三峡库区面源污染有效控制和治理提供科学依据。
梅溪河流域(108.93°—109.55°E,31.03°—31.53°N)和大宁河流域(108.73°—110.18°E,31.06°—31.73°N)位于三峡库区腹地,地理位置毗邻,水系都起源于巫溪县,向南分别经奉节县和巫山县汇入长江。2个流域均位于大巴山弧形褶皱带,岩石发育广泛,以山地地貌为主,地形条件复杂,属亚热带湿润季风气候,4—9月降水量占60%~80%。2个流域土地利用类型均以林地为主(60%和72%),其次为旱地(31%和18%)。耕作方式相似,属传统农耕区,耕作历史悠久,坡耕地分布范围较广,也是典型的水土流失区。虽然2个流域自然地理环境具有一定的相似性,但各自特点鲜明。梅溪河全长112.8 km,流域面积1 898 km,为喀斯特—非喀斯特组合流域,砂岩和石灰岩分布广泛,土壤以紫色土为主。大宁河全长162 km,流域面积4 420 km,为喀斯特流域,石灰岩广布,土壤以黄壤为主。此外,梅溪河流域海拔相对较低,人口众多,农业活动较大宁河流域更加发达。因此,本研究所选区域具有典型性和代表性。
SWAT模型输入数据主要包括数字高程模型(DEM)、土地利用/覆被数据、土壤数据、气象数据、水文水质数据以及当地农业管理措施数据等,数据来源见表1。
表1 梅溪河和大宁河流域基础数据及来源
DEM用于划分流域,在ArcSWAT 2012平台中,梅溪河和大宁河流域先被划分为27,37个子流域,随后根据土地利用、土壤类型和坡度的分布特征进一步将流域细分为2 005和2 812个水文响应单元(HRU)。HRU具有唯一的土地利用类型、土壤类型和坡度等级,是模型基本计算单元。SWAT模型中径流模拟分地表径流、地下水和河道汇流进行。NO—N被模拟在地表水中运移,也可通过渗漏和地下水一起运移。本文选择soil conservation’s curve(SCS)曲线计算地表径流,地下水则分别通过对浅层和深层含水层的水量平衡方程来计算。由于CMADS数据集可提供全部气象要素,因此选择Penman-Monteith法计算潜在蒸散发。SWAT模型更多细节见文献[20]。
本研究利用SWAT自动校准与不确定性分析程序SWAT-CUP 2012内置的连续不确定率定算法(SUFI-2)对月径流量和NO—N进行校准。根据数据可用性,2008—2010年为模型预热期,径流校准期为2011—2015年,验证期为2016—2018年,梅溪河流域NO—N校准期为2013年和2017年,验证期为2018年。由于梅溪河流域无水文站,本研究选用水利工程中较为常用的水文比拟法推求得到梅溪河流域逐日流量数据。水文比拟法是指将气候、下垫面等条件相似的参证流域断面水文资料移用至待求流域,考虑流域面积、降雨量来计算待求断面径流量,已被广泛应用于少资料、无资料地区的径流量计算,计算方法见文献[21]。本文选择大宁河流域巫溪水文站作为参考断面,计算2008—2018年梅溪河库湾的逐日来流过程,其结果已在之前的研究中验证。
通过全局敏感性分析确定模拟结果对特定参数的敏感性程度,SWAT-CUP中统计量(-Stat)和值作为参数的敏感性指标。选取决定系数()和Nash-Sutcliffe Efficiency(NSE)效率系数作为模型评估指标。和NSE分别反映模拟值与观测值之间相关程度和拟合程度,被广泛应用于模型性能评估。综合SWAT模型在国内外的应用情况,一般认为,当>0.5,NSE>0.5时模拟结果是可接受的。本文采用Moriasi等在综合大量文献基础上提出的评价标准,即模型月尺度模拟结果NSE≥0.65,说明模型模拟效果是可接受的。
鉴于无法获得大宁河流域水质监测数据,因此选择参数扩展的方法进行大宁河NO—N模拟,即在数据充足的流域进行模型校准和验证,然后将这些典型流域的参数扩展到周边地区,从而获得整个流域的面源污染负荷。该方法是水文资料匮乏地区进行水文模拟的有效方法,已在三峡库区取得良好效果。在先前研究基础上,本文基于空间临近法将梅溪河流域校准参数扩展至大宁河流域。
1.5.1 基流分割 SWAT模型是一种常用的基流分离工具,特别是对于无资料流域。在SWAT中,总径流量是一定时间步长内离开HRU进入主河道的总水量,计算公式为:
WYLD=SURQ+LATQ+GWQ-TLOSS
式中:WYLD为总径流量(mm);SURQ为地表径流(mm);LATQ为壤中流或侧向流(mm);GWQ为地下径流(mm);TLOSS为传输损失(mm)。
NO—N的计算则为通过地表径流、壤中流和地下径流输送到河流的NO—N的总和减去传输损失。模型将地下水分为浅层和深层地下水,并将深层地下水作为流域损耗排泄到流域外。因此,SWAT模拟的基流对应于在某个时间步长内从浅含水层进入主河道的水。本研究借鉴Zhu等的计算方法,将模拟的地表径流从总径流中减去得到基流,即基流包括地下径流和壤中流,这也是国际上普遍采用的径流划分概念。最终通过HRU输出文件统计得到不同土地利用类型的径流和NO—N负荷,以及不同季节下不同径流组成对NO—N负荷的贡献。根据当地气象记录和水文季节分离标准,雨季定义为4—9月,旱季定义为10月至翌年3月。
1.5.2 随机森林回归算法 本文研究土地利用、气候、土壤、地形对径流和NO—N的影响,考虑到变量之间的多重共线性和非线性关系,利用随机森林回归量化不同因素的相对贡献。随机森林是一种非参数统计方法,相对于其他回归方法具有的优点为:(1)无需预先设定函数形式,模型结构灵活、预测精度高;(2)不受线性回归假设条件限制,可处理变量多元共线性问题;(3)引入随机性使模型不易过拟合,抗噪声能力强;(4)袋外数据的存在使模型运算无需单独留出验证集。随机森林算法具体原理见文献[24]。
以土地利用、土壤类型、坡度、坡长、降水和潜在蒸散发为自变量,通过R语言“randomForest”程序包构建径流和NO—N负荷的随机森林模型。SWAT模型划分的每个HRU作为1个样本。随机森林模型中回归树数量()和分裂次数()是2个重要的自定义参数。随着值增大,袋外误差显著降低后基本保持稳定。值通常选取于自变量个数的1/3处。本研究通过tuneRF函数调优模型获得最佳参数。最终取值为3,取值为600。模型通过计算变量排列时均方误差(mean square error, MSE)的增加来估计每个预测变量的相对重要性,每个变量的相对重要性与MSE的增加成正比。模型显著性和变量显著性通过“A3”和“rfPermute”程序包进行(1 000次)。基于MSE的增加将预测变量的相对重要性归一化为100%。
本文选择18个参数来校准流量,11个参数来校准NO—N,29个参数的率定值见表2。梅溪河流域径流和NO—N在校准期和验证期最佳模拟值的均大于0.73,NSE均大于0.68(图1a~图1d),大宁河流域径流模拟和NSE均大于0.69(图1e~图1f),表明模拟结果较好,可进一步用于解析流域径流和NO—N输出特征。将梅溪河流域参数率定值直接扩展至大宁河流域,用模拟流量与观测流量进行比较。结果显示,两者之间存在较强的线性关系(=0.87,<0.01,图1g),说明相邻流域之间的参数扩展是可行的。但就2个流域径流参数率定结果来看,参数敏感性等级和取值仍具有一定差异。这些结果既表明特定于区域的属性,又表明2个流域模型参数存在“异参同效”现象。为保证径流过程的准确性,本文最终以大宁河实际率定参数模拟径流过程,仅将梅溪河流域NO—N参数扩展至大宁河流域。
图1 径流与硝态氮的校准与验证
表2 参数敏感性分析及校准值
总径流量在不同土地利用类型之间相差不大(图2a),多年平均值波动于333~677 mm。建设用地径流量最高,旱地和水田次之,林地、草地、园地相对较低。由于建设用地不透水表面大幅削弱其持水蓄水能力,因此产水量最高。其他透水类型中,旱地和水田裸露表面较多,且与林地相比表现较低的蒸散,因此产水量相对较高。
NO—N负荷强度在不同土地利用类型间差异显著(图2b)。负荷强度排序为园地(20.41 kg/hm)>旱地(12.51 kg/hm)>水田(10.31 kg/hm)>建设用地(7.09 kg/hm)>林地(0.62 kg/hm)>草地(0.46 kg/hm)。从输出总量(图2c、图2d)来看,旱地、水田和园地等农业用地是NO—N主要来源,贡献了梅溪河流域90%和大宁河流域74%的NO—N输出,其中又以旱地对NO—N的贡献最大,占梅溪河和大宁河流域总负荷的80%和67%。此外,由于具有最高的负荷强度,仅约占2个流域面积2%的园地贡献6%的NO—N输出。
注:(a)(b)分图中不同字母表示不同土地利用类型间差异显著(p<0.05)。
梅溪河和大宁河流域水文状况具有一定相似性(表3)。2个流域多年径流深分别为560,540 mm,总径流系数分别为45%和48%。2个流域地表径流处于较低水平,基流较丰富。梅溪河流域基流系数(67%)略高于大宁河流域(62%)。
表3 梅溪河和大宁河流域产流总体特征
从图3可以看出,径流补给来源表现出明显的季节性和不均匀性。旱季,基流是梅溪河和大宁河流域径流量的主要补给来源(图3a、图3c)。除建设用地外,其他5种透水类型基流占总径流量的70%以上。流域尺度上,旱季基流占2个流域总径流量的92%(梅溪河)和84%(大宁河)。雨季,地表径流是径流量的重要组成部分(图3b、图3d)。雨季地表径流约占2个流域全年地表径流的89%(梅溪河)和88%(大宁河)。雨季各土地利用类型地表径流对总径流的贡献明显增加,由旱季的6%(林地)~66%(建设用地)增长至22%~90%。地表径流和基流在不同土地利用类型间差异较明显,建设用地地表径流最大,其次为旱地和水田,草地、园地、林地较小;而基流恰好相反,雨季草地、园地、林地基流仍占总径流的65%以上,其中林地基流占比最高,为78%。这是由于林地土壤物理结构良好、透水性强,且地表凋落物具有较强持水性能。
图3 不同土地利用类型间径流组分变化
雨季不同土地利用类型的NO—N负荷与旱季相比增长2~7倍(图4),雨季大约贡献2个流域全年77%(梅溪河)和86%(大宁河)的NO—N负荷。旱季,基流是NO—N的主要输出途径,输出了2个流域85%(梅溪河)和72%(大宁河)的NO—N。建设用地NO—N主要通过地表径流输出,其他土地利用类型中基流对NO—N的贡献在67%(草地)~99%(水田)。雨季,随着地表径流量增加,地表径流对NO—N的贡献明显增加,尤其是旱地和水田。梅溪河流域旱地和水田地表径流对NO—N的贡献由旱季的16%和1%增长至雨季的43%和19%;大宁河流域由19%和6%增长至61%和45%。此期间,地表径流共输出36%和42%的NO—N。
图4 不同径流组成对硝态氮的贡献
随机森林模型结果(图5、图6)表明,所选的6个变量可以很好地解释径流和NO—N负荷(<0.001)。2个流域中不同因素对径流和NO—N负荷的影响具有一致性。降雨量和潜在蒸散发是影响总径流的主要因子,二者相对重要性达78%以上,其中降雨量相对重要性达58%(图5)。对于地表径流和基流,土壤类型是其主要影响因素,相对重要性为45%~54%。土地利用、降雨量和潜在蒸散发对地表径流的相对重要性大致相等,为14%~17%。
土地利用和土壤类型是影响不同来源NO—N的主要因子,二者相对重要性达70%以上(图6)。其中,土地利用的影响最大,相对重要性为43%~58%。地形因子对径流和NO—N的相对影响较小。
图6 随机森林模型中影响硝态氮自变量的相对重要性
本研究结果表明,基流是梅溪河和大宁河流域径流的主要补给来源,分别占全年总径流量的67%和62%,这与Zhu等在大宁河流域(66%)研究结果一致。研究流域处于喀斯特地貌区,受到独特的地表地下二元水文结构影响,使得地表水易漏失至地下。此外,土壤是影响基流的重要因素(图5)。梅溪河和大宁河流域土壤水文分组大多为A组和B组,其均有团聚体结构少、孔隙度大、土壤压实度低和入渗能力强的特点,易冲刷流失,也易造成雨水渗漏。特别是紫色土区域底部页岩隔水层的存在使土壤水分易于下渗,从而造成壤中流的普遍存在。熊子怡等观测到三峡库区紫色土坡耕地壤中流流量高达总径流量的60%~89%。虽然研究流域旱季降雨量小,但基流量仍较高,这可能与壤中流的滞后性有关。再者,降雨也是影响水文过程的主要因素之一,这不仅体现在降雨总量上,也与降雨强度和降雨历时密切相关。
图5 随机森林模型中影响径流自变量的相对重要性
三峡库区腹地全年雨强>1.00 mm/min的降雨占总降雨次数的比例不足10%,且绝大多数<0.5 mm/min。这种小雨强、长历时的降雨特征产流方式主要以蓄满产流为主,更有利于降雨向基流的分配。因此,地质结构、土壤性质及降雨特征决定研究区基流系数较高。
由于流域径流输出方式以基流为主,因此基流是研究流域NO—N输出的主要途径,特别是旱季。在其他由基流主导水文过程的流域中均发现一致的结论。这一现象除与丰富的基流流量有关外,也与基流中养分浓度较高有关。贾海燕等研究表明,无论施肥与否,紫色土壤中流中NO—N浓度均高于地表径流。壤中流形成在土壤的浅层且在降雨结束后还具有较长持续时间,最终沿风化裂隙继续向下运动,特别是大孔隙流的流速较快,造成大量氮素来不及与土壤发生作用而随水进入地下水,导致地下水中NO—N含量显著增加。陈成龙等在三峡库区农业小流域观测到不同土地利用类型坡面浅层地下水NO—N对TN的贡献为68%~79%。目前,三峡库区经过水土保持治理后地表径流和土壤侵蚀已得到有效控制。李红颖等在三峡库区砂质土梯田柑橘园中研究发现,坡改梯和生草覆盖管理模式下渗漏是养分流失的主要途径,分别占TN和TP流失量的99%和77%。此外,壤中流和地下水循环的滞后性导致基流养分输出不与地表径流同源同步,这无疑加大面源污染的治理难度。因此,当前面源污染防控策略亟待纳入基流传输路径削减策略。
本研究表明,土地利用是影响不同径流组分中NO—N的主要因素。陈成龙等研究发现,三峡库区浅层地下水中NO—N受下渗坡面土地利用类型的显著影响。旱地是研究流域乃至三峡库区面源污染的主要来源地。一方面,大部分旱地位于结构较差、有机质含量较低的陡坡地区;另一方面,人地矛盾突出客观上推动剩余农田化肥施用量进一步增大以应对巨大的粮食压力,因为增加无机肥料投入是提高区域作物产量的主要因素。从负荷强度来看,3种农业用地呈现出园地>旱地>水田的特点,这与朱波等在三峡库区石盘丘小流域研究结果一致。虽然之前的研究表明,经济果林保证土地持续的植被覆盖,从而提高水分利用效率,减少土壤侵蚀同时兼顾经济效益。但在实地调查中发现,园地(主要是柑橘)化肥施用量是农田的3.3倍。鉴于经济林种植的日益扩大和较高的负荷强度,园地将是三峡库区面源污染治理中另一关键的土地利用类型。此外,土壤类型也是影响NO—N产生的重要因素(图6)。土壤性质是影响近地表地球物理过程的关键因素。研究流域2种主要土壤类型紫色土和黄壤在整个三峡库区土壤类型中具有较高的污染负荷强度,在负荷总量方面也占主导地位,这不仅与2种土壤类型覆盖范围广有关,也与长期施肥以及土壤和地下水中养分停留时间较长有关。相关研究表明,氮素的淋失量和残留量与施肥量呈显著正相关。NO—N通过淋溶离开根区后大量NO—N积累储存在包气带中,导致NO—N进入地下水的时间延长,这些遗留在包气带中的氮素可以通过基流过程的不断缓慢输出作为地表和地下水污染的长期来源。由此可见,在高渗透性土壤上的农业土地利用是造成河流NO—N增加的主要原因。这种土地利用不符合土壤性质(如质地、有机质含量)或环境条件(如坡度)所规定的土地能力,可定义为环境土地利用冲突(environmental land use conflicts)。Pacheco等调查了葡萄牙大陆环境土地利用冲突表明,环境土地利用冲突是造成农业流域NO—N富集的主要原因,特别是在葡萄种植园大量侵占陡坡区域的流域。而消除土地利用冲突则是缓解淡水NO—N污染的有效途径。SWAT模拟结果显示林地等自然土地类型NO—N负荷显著低于其他土地利用,可见退耕还林是消除环境土地利用冲突的途径之一。然而,在人地矛盾突出的山区,通过牺牲农用地为代价的调控模式已日趋困难。因此,源头控制仍然是三峡库区面源污染治理的关键。
水文模型的参数选择和调参是水文建模中的重要环节。本文所选参数(表2)与三峡库区相关研究较为一致,其中大宁河部分参数取值(ALPHA_BF、GW_DELAY、CH_K2、SLSUBBSN、REVAPMN)与Shen等在该流域率定结果相似,NO—N部分参数(CDN、SDNCO、RCN、BC1~3)与石荧原在梅溪河率定结果相似。总体上,参数敏感性分析结果显示出与其他研究的共同特征。但是,参数敏感等级和部分参数取值与前人研究具有差异。例如,CN2是地表径流模拟中的重要参数。然而,CN2在本研究中敏感性等级相对较低,这一结果与Shen等在大宁河流域研究不一致,但与刘伟等在大宁河流域研究及袁江等在其他喀斯特流域研究结果相似。这些不确定性源自不同的敏感性分析方法、不同数据源(土地利用、土壤、气象数据等)造成的背景偏差以及不同参数之间的交互作用和参数敏感性的时变特征。因此,参数敏感性分析需因地制宜,而不能普遍化。本研究所选参数中SOL_BD、SOL_K、SOL_AWC决定水分在土壤层的分布及下渗过程,ALPHA_BNK、GWQMN、GW_DELAY、REVAPMN、ALPHA_BF、CH_K2、CH_N2等决定了地下水分布及地表水—地下水转化过程。以上参数体现了研究区壤中流现象普遍和地下水丰富的特点,与实际水文过程较为符合。
从模拟结果来看,参数扩展取得较满意的效果,但2个流域径流控制参数具有一定差异。这是2个流域不同的土地利用和土壤组成及地理环境异质性导致的。因此,应用参数区域化方法进行无资料流域水文预测并提升预测准确性仍是未来研究的重大挑战。另外,2个流域模型参数值不一致却产生相似输出现象表明研究区存在较高“异参同效”性。这一现象主要源自水文模型的非线性、参数相关性和模型公式的设计缺陷,是当前水文模型参数率定中普遍存在且亟待解决的难题。建议未来应用模型时使用实测流域特征值约束模型以降低“异参同效”风险,提高流域模型对水文过程的模拟性能。
(1)NO—N负荷在不同土地利用类型间差异显著,年负荷强度大小顺序为园地(20.41 kg/hm)>旱地(12.51 kg/hm)>水田(10.31 kg/hm)>建设用地(7.09 kg/hm)>林地(0.62 kg/hm)>草地(0.46 kg/hm),其中旱地对2个流域NO—N输出贡献最高(80%和67%)。从负荷强度和输出总量来看,旱地和园地是三峡库区面源污染防控的关键土地利用类型。
(2)基流是梅溪河和大宁河径流主要输出方式和NO—N的主要输出途径,分别占2个流域总径流的67%和62%,输出NO—N负荷分别占总负荷的68%和60%。但径流和NO—N负荷具有明显的季节变化和不均匀性。旱季径流量主要由基流补给,地表径流是雨季径流量的重要组成部分。在集约化水土流失生态治理背景下,基流养分流失已经成为三峡库区面源污染的一个不可忽视的污染源。当前面源污染治理策略不应仅限于严格控制地表径流的养分流失,应同时考虑地表径流和基流途径,研发流域综合管理体系。
(3)降雨量是控制总径流的主要因素,而土壤类型是影响地表径流和基流的主要因素。不同径流来源NO—N负荷主要受土地利用影响,其次为土壤类型,二者相对重要性在70%以上。综合来看,在高渗透性土壤上的农业土地利用,即环境土地利用冲突是造成NO—N流失的主要原因。因此,以减少化肥投入和提高养分利用效率为主的各种源头控制措施仍需广泛推广。