黄土区裸露坡地径流养分流失模型的建立与验证*

邵凡凡，吴军虎，李玉晨

（西安理工大学省部共建西北旱区生态水利国家重点实验室，西安 710048）

溶质从表层土壤进入地表径流的过程非常复杂。雨滴击溅、土壤中的养分浓度以及对流、扩散作用和土壤颗粒吸附均会对这一过程产生影响[1-3]。Ahuja等[4]1981年通过在饱和土层不同深度位置放置一定量的32P，发现表层土壤的溶质进入径流的可能性最高，随着土层深度的增加，进入径流的概率呈指数递减。然后，Ahuja和Lehman[5]在1983年发现，观察到的交换层深度较通过拟合模型和实测数据获得的深度小得多。结合交换层理论和Rose土壤侵蚀模型，Gao等[6]在2003年建立了基于物理过程并考虑雨滴击溅和扩散作用的溶质运移模型；该模型的所有参数均可通过实验测量获得。鉴于黄土高原地区养分流失的特点，王全九等[7]提出了一种新的方法来改进等效对流传质模型。在此基础上，Dong等[8]假设交换层被混合层代替，且交换率被雨滴诱导水分转移率所代替。Yang等[9]结合质量守恒方程和降雨对表层土壤的剥蚀过程，建立了可预测黄土高原地表养分随径流流失的数学模型。但是，这些模型无法描述径流发生之前土壤溶质浓度的变化，并且仅能通过拟合曲线来获得初始土壤溶质浓度。Tong等[10]基于质量守恒方程和水平衡方程建立了一维两层的溶质运移模型。该模型结合了入渗和扩散作用，并由不完全混合参数来描述。基于该模型，Tong等[11]使用集合卡尔曼滤波数据同化方法（EnKF）来校准参数并更新可溶性化学物质从土壤至地表径流的转移过程，并消除了实验观测数据的误差。因此，基于混合层理论的模型由于具有明晰的物理意义而被广泛用于预测斜坡上的溶质运移。

养分流失的模拟是在径流模拟的基础上进行的。通常使用 Saint-Venant方程（即连续性方程和动量方程）来描述地表径流过程[12]。然而，由于Saint-Venant方程是高度非线性的，很难获得解析解，这意味着仅能使用数值方法对其进行求解。但当忽略 Saint-Venant方程的加速度项时，可使用扩散波方程来对其进行简化[13]。若同时忽略Saint-Venant方程的加速度和压力项时，Saint-Venant方程可表示为运动波方程。Luce和Cundy[14]通过使用菲利普（Philip）入渗方程修改了运动波方程来预测超渗降雨条件下的产流过程。Yang等[15]通过假设水深与入渗率之间的线性关系简化运动波模型中的水深项，并结合Philip入渗方程得到了运动波模型的近似解析解。该模型因其参数简单易获取而被广泛使用于坡面径流的模拟。在养分随径流流失的研究中，Gao等[16]在2004年提出了一个基于溶质守恒方程的模型，该模型考虑了雨滴飞溅和径流冲刷作用。但是，该模型仅用于模拟积水条件下饱和土壤的养分流失过程，这与黄土区初始非饱和土壤条件下的流失过程存在较大差异。因此，本文以 Yang等[15]建立的坡面径流近似解析解为基础，进一步延伸于养分随地表径流流失过程的模拟中，并修改了Gao等[16]的模型以适应本文的测试条件。通过模型参数分析揭示了不同因素对养分流失的贡献作用，提出了防止养分流失的有效措施。该研究可为防治农田退化和农业面源污染提供有力基础。

1 材料与方法

1.1 理论与模型

1.1.1 坡面径流运动过程 采用运动波模型来描述次降雨条件下的坡面水流流动过程[15]，其中超渗净雨可用降雨强度与入渗率的差值来表示，见式（1）：

式中，h为径流深，cm；t为径流时间，min；q为单宽流量，cm2·min–1；x为坡面任一位置距离入流口的长度，cm；p为降雨强度，mm·h–1；i为土壤入渗率，cm·min–1。

由于坡面水深与入渗率之间存在关联关系[9]，径流水量为超渗净雨所产生的，Yang等[15]用Philip公式表示降雨条件下的入渗过程，进一步求解了单宽流量和坡面水深，见式（2）和式（3）：

式中，c为入渗率参数[15]；S为土壤吸渗率，cm·min–0.5；Δt=3S2/（16p2）。

式中，n为曼宁糙率系数，；S0为水力梯度，本文中坡度为15°，故S0为Sin 15°。

1.1.2 径流养分流失过程 降雨条件下土壤表层养分在雨滴击溅和水分入渗的作用下随径流迁移并在土壤中重新分配，因此土壤剖面的水和养分运移系统自上而下可分为3层：径流积水层、养分交换层和交换层以下土壤，如图1所示。

养分交换层是径流积水层与土壤剖面交界面以下厚度较薄的土层。交换层中化学物质的传输主要受入渗、水动力弥散和雨滴飞溅侵蚀控制[16-17]。

式中，de为交换层深度，cm；Ce为交换层中溶质浓度，mg·L–1；Cw为径流中溶质浓度，mg·L–1；is为径流水进入交换层的入渗速率，cm·min–1；ix为交换层水分进入更深土层的入渗速率，cm·min–1；er为雨滴诱导水分转移速率，cm·min–1；λCw为径流层进入交换层的溶质浓度（0 ≤λ≤ 1，Gao等[6]研究表明计算模型对该参数不敏感，取λ=0），mg·L–1；J为较深土壤层与交换层内部的溶质扩散通量，mg·cm–2·min–1。

为了简化计算过程，Gao等[16]对J进行近似求解如下：

式中，Ds为养分在土壤中的扩散性，cm2·min–1；Cs为更深土层的溶质浓度，mg·g–1；γ为土壤容重，g·cm–3；K为土壤吸附系数，mL·g–1；β=er/（αde）。

从降雨开始，将整个降雨过程划分为3个阶段。

第一阶段：从降雨开始t0至交换层完全饱和tsa。在这一阶段，土壤入渗率为降雨强度，土壤表层未产生径流，故i=p，q=0，交换层完全饱和的时间tsa可以表示为：

式中，tsa为交换层完全饱和所需时间，min；θs为饱和含水率，cm3·cm–3；θ0为初始含水率，cm3·cm–3。

第二阶段：从交换层完全饱和tsa至土壤表层出现径流tp。在这个阶段，径流层溶质浓度Cw和雨滴诱导水分转移速率er的取值为0，ix = p。交换层中溶质浓度见式（9）：

式中，C0为初始溶质浓度，mg·L–1。

当式（9）中t=tp时，开始产流时刻交换层中的溶质浓度见式（10）：

式中，A=ix/（αde）。

第三阶段：从开始产流至降雨结束。这一过程中，径流中的养分浓度远低于交换层中的养分浓度，因此为了简化计算过程，忽略了径流养分在入渗作用下对交换层养分的微小补给作用，ix=0.01 cm·min–1。结合这一阶段的起始产流时间，即：t=tp，可求解得到交换层中溶质浓度表示如下：

式中，B=（ix+er）/（αde）。

产流过程中，径流中化学溶质的质量守恒关系可表示为：

结合式（1）和式（12），可得到：

为了简化计算，忽略了入渗和扩散作用，式（13）可表示为：

对式（14）进行积分得到径流中化学溶质的浓度，见式（15）：

径流中化学溶质的浓度与径流量的乘积便为径流中化学溶质的流失速率，结合式（2）、式（3）和式（15），径流中溶质流失速率可表示为：

式中，Mw为养分随径流流失的速率，mg·min–1。

1.2 试验区概况

试验于2019年5月在中国科学院水利部水土保持研究所长武黄土高原农业生态试验站担水沟流域野外模拟降雨小区（35°12′N，107°10′E）进行，试验区平均海拔为1 200 m，气候属暖温带半湿润大陆性季风气候，年平均气温为 9.1℃，年平均降水量580 mm，地下水位50～80 m，无灌溉条件，属典型的旱作雨养农业区。该流域内典型土壤为粉砂质壤土，母质为深厚的中壤质马兰黄土，具体土壤物理特性见表1，流域塬面面积占35%，梁坡占35.6%，沟谷占29.4%，各约占1/3；流域地貌属典型的黄土高原沟壑区[18-20]。

表1 试验区土壤物理化学特性Table 1 Physical and chemical properties of the soil tested

1.3 试验方法

试验用地为3年闲置坡耕地，模拟降雨试验的小区设置尺寸为1.0 m×1.0 m，并根据当地典型坡耕地坡度和坡面侵蚀的临界坡度，设置小区坡度为15°；使用针孔式人工模拟降雨装置进行降雨试验（图2），其主要由主体支架（可调节高度）、底板布设有针孔的水槽（可根据降雨强度更换不同孔径的针孔）和供水装置 3部分组成，有效降雨面积为1.0 m2。经过测试：该套人工降雨器的平均雨滴直径为2 mm，降雨均匀度在80%以上，雨滴终速符合天然降雨特征[3]。

为了消除土壤前期含水率对试验结果的影响，每次开始降雨24 h前以25 mm·h–1的降雨强度在试验小区进行预降雨，直至开始产流时停止降雨。开始正式降雨试验前测定小区内土壤表层 0～20 cm剖面的初始含水率，采用烘干法测定质量含水率为0.11± 0.003 7 g·g–1（ 即 体 积 含 水 率 为 0.15±0.005 cm3·cm–3）时开始试验。为了提高土壤初始养分浓度值，使养分在土壤中均匀分布，预降雨结束后，在小区土壤表面均匀喷洒氯化铵和硝酸钾混合溶液，其具体操作方式为：将预先配置好的5.0 g·L–1的氯化铵溶液和10.0 g·L–1的硝酸钾溶液各取1 L混合均匀（为消除溶液喷施次序对其分布的影响），并用压力喷壶在每个小区分别定量喷洒2 L氯化铵和硝酸钾的混合溶液，为尽可能减小喷壶压力对表层土壤造成的压实作用，将喷嘴调节至雾化度最强位置处，喷嘴雾化半径为5 cm，采用左右往复的方式将混合溶液均匀喷洒在小区表层，使其在2 cm范围内均匀分布。并在小区坡面上、中、下三个部位分别取土样来测定表层 2 cm土壤中的溶质浓度作为模型计算所需的初始土壤养分浓度（表2）。根据研究区暴雨实测资料及降雨分级标准[20]，设计30、45、60、75、90 mm·h–1的 5种降雨强度，设计总降雨历时为120 min，按照0～10 min之间，每隔2 min承接1次径流，10～120 min之间，每隔5 min承接1次径流的频率用量杯承接出口处径流，并用量筒进行精确测量，通过沉淀过滤除去径流中的泥沙，用50 mL的塑料瓶收集径流水样并存放于实验室冰箱中，用全自动高通量间断分析仪（SmartChem450，AMS Allinace公司，意大利）测定径流和土壤中的养分浓度。

1.4 模型基本参数

通过现场测定和文献查阅等方式获取了模型计算所需参数[20-21]（表2）：

表2 模型基本参数Table 2 Basic parameters of the calculation model

1.5 数据处理

所有试验实测数据均为 3次重复试验的平均值，使用Matlab 2015b进行参数求解和模型模拟；使用SPSS 24.0进行数据分析，使用Origin 2018进行图表绘制和函数拟合。

2 结果与讨论

2.1 产流过程分析及模拟

2.1.1 产流过程分析 土壤表层的养分通常会随着地表径流而流失。因此，探究产流规律是模型准确预测养分随径流流失过程的基础。在5种降雨强度下，开始产流的时间点分别为20.5、8.5、4.8、3.0和 1.8 min（图3）；90 mm·h−1较 30 mm·h−1提前19 min产流；说明随降雨强度的增大，开始产流时间开始显著缩短。降雨强度与产流时间的关系可用幂函数来描述，R2=0.997 6。单宽流量开始产流后快速增大，而后进入稳定产流阶段（图4）。这是由于表层土体中的黏粒分散堵塞了土壤的孔隙，并伴随着雨滴的飞溅使表层土壤变得密实，降低了土壤的入渗能力[13]。5种降雨强度下实测单宽流量的标准差分别为 0.01～0.16、0.01～0.16、0.10～0.41、0.13～0.71和0.23～0.69。这可能是由于土壤的非均质性、蚁穴和植物根系对小区土壤入渗过程的影响，以及雨滴飞溅和径流发育过程中微地形的形成，可能导致径流滞后。进一步分析发现，90 mm·h−1在稳定产流阶段的单宽流量分别较其他降雨强度依次增加6.3倍、2.7倍、1.6倍和1.2倍，这表明降雨强度的增加显著增大了坡面径流率。

2.1.2 产流过程模拟 通过将已知参数S、p（表2）和实测单宽流量代入式（2）中推求入渗率参数c，采用R2、均方根误差（RMSE）和纳什效率系数（NSE）对模拟结果进行评价。从表3可看出，入渗率参数c随着降雨强度的增大呈减小趋势，且分布在0.003 1～0.006 0之间；R2均在0.89以上，随着降雨强度的增大，RMSE也随之增大，取值分布在0.406～1.052之间，NSE均大于0.397，而当降雨强度大于等于 60 mm·h–1时，NSE 则进一步增大至0.783以上；说明降雨强度越大，模型计算值和实测值的匹配度也随之提高。指数函数可很好地拟合参数c与降雨强度之间的关系（图5），决定系数R2为0.976 5，表达式为：

表3 入渗率参数c的最佳拟合值Table 3 Optimal fitting values of c，R2，RMSE and NSE relative to rainfall intensity

图6分别显示了5种降雨强度下的单宽流量模拟过程。可以看出，产流模型能够较好地模拟地表径流过程，且随着降雨强度的增大，模拟趋势变得更加准确。在产流初期，计算值的上升趋势均慢于实测值；在稳定产流阶段，30 和45 mm·h−1下的模拟值均大于实测值，而在 60、75、90 mm·h−1下，实测值与计算值的匹配程度较好。这可能是由于降雨强度较小时，雨滴动能的溅蚀增加了前期土壤表面粗糙度；雨滴击溅形成的微地形和洼地拦截部分地面径流，从而削减了连续径流的冲刷作用。这表明产流模型可准确模拟大于等于 60 mm·h−1的产流过程。此外，雨滴的击溅使表层土壤被压实，容重增大，土壤表层形成密封层，降低了土壤的入渗能力。但该模型未考虑地表土壤容重和孔隙率的变化，导致模型计算的土壤入渗能力明显大于实测值。

2.2 养分流失过程分析及模拟

2.2.1 养分随径流流失过程分析 养分从土壤至径流的传输是通过雨滴击溅作用和径流溶解作用来完成的[16]。不同降雨强度下硝态氮和铵态氮随时间的流失过程可用单峰形式来描述；即：径流初期养分流失速率迅速增大，到达峰值后开始减少，最后进入稳定流失阶段的趋势（图7）。这是由于土壤水和土壤颗粒表面吸附的养分在雨滴击溅作用下进入径流所引起的；随着土壤结皮厚度和径流深度的增大，土壤表层形成“坚实的保护壳”削弱了雨滴动能，延缓土壤水和径流的交换作用，使得进入径流的土壤水和溶解态氮显著减少。此外，径流中氮素浓度的降低也是由于随着降雨过程的推移，表层土壤中氮含量的逐渐减少所引起的。同时，降雨强度、入渗能力、养分浓度和径流率对养分的峰值流失速率及其发生时间均有一定的影响。简言之，硝态氮和铵态氮的峰值损失率随降雨强度的增大而增大。以硝态氮流失过程为例，当降雨强度为 30 mm·h–1时，硝态氮流失速率在 25 min左右达到峰值5.74 mg·min–1，而当降雨强度为 45、60、75 和90 mm·h–1时，分别在 13、8、6和 5 min达到硝态氮流失速率的峰值：35.21、121.3、280.4和468.4 mg·min–1。因此，降雨强度对硝态氮的峰值流失速率具有较大贡献。通过对比养分的峰值流失速率出现时间和稳定产流时间，各降雨强度下养分流失速率峰现时间分别为25、13、8、6、5 min，而稳定产流时间分别为40、28、18、16、14 min；由此看出，硝态氮流失速率的峰值出现时间要早于径流速率达到稳定阶段所需的时间，这可能是随着产流时间的推移，交换层土壤中硝态氮浓度的不断减小和径流量的增大共同作用所造成的。5种降雨强度下硝态氮损失率的标准误差分别分布在0.03～0.53、0.05～5.13、0.03～6.9、0.1～21.0 和 1.1～31.7。

2.2.2 养分流失过程模拟 交换层深度de和雨滴诱导水分转移速率er是养分流失模型中的两个重要参数。由于受室外实验条件的限制，交换层深度de很难通过实地测量得到，因此需借助模型拟合实测的养分流失速率来反推交换层深度de。研究[4-5]发现交换层的深度de在2～3 mm的范围内。Tong等[11]指出，交换层深度随入渗率的增加而减小。有研究[21-23]指出，交换层深度随着初始含水量的增大而增大。关于雨滴诱导水分转移速率er，Gao等[16]在2004年提出了适用于初始饱和土壤的雨滴诱导水分转移率er的计算方法，由于黄土区坡耕地在降雨前为非饱和土壤，因此，借助 Matlab非线性拟合的方法，将式（3）所计算出的坡面径流深度代入式（16）来计算养分流失速率，并结合硝态氮流失速率的实测值推求出了式（16）中的交换层深度de和雨滴诱导水分转移速率er（表4），并进一步模拟了铵态氮流失过程。可以看出：de和er均随着降雨强度的增大而增大，其分别从0.68增至1.32、从0.006增至0.023。这与前述单宽流量随着降雨强度的增大而增加是一致的[24]。硝态氮和铵态氮流失速率的R2值分别分布于 0.834～0.922和 0.800～0.921之间，RMSE值分别分布在1.188～58.50和0.974～58.37之间，NSE值分别分布在0.653～0.881和0.546～0.775之间。从图8中可以看出，计算出的曲线可很好地模拟养分流失过程。当降雨强度为30 mm·h–1时，初始增大阶段的测量值与计算值之间的差异较大，但随着降雨强度的增大，差异逐渐减小。而在养分流失的稳定阶段，实测数据大于计算值。这可能是由于本文建立的养分流失模型近似求解了对流弥散项所造成的，这使得稳定减小阶段养分流失速率衰减得过快[9]。同时模型忽略了径流层养分对交换层的微弱补给作用，然而在降雨开始时可能存在从交换层至径流层的扩散过程[16]。从交换层完全饱和（tsa）到坡面开始产流（tp）的时间段内，模型假设交换层中的养分随入渗水向土壤深处迁移的速率大小即为对应的降雨强度，从而使得计算出的开始产流时土壤表层的浓度Ce（tp）小于理论值。

表4 不同降雨强度下参数de和er的最佳拟合值Table 4 Optimal fitting values of de，er，R2，RMSE and NSE relative to rainfall intensity

以上分析表明本文建立的模型可以很好地模拟裸露坡面径流和养分随径流迁移过程。但该模型未用于模拟不同坡度、坡长和土壤初始含水量条件下的养分流失过程，本研究所获得的参数是否具有普遍适用性，需要在以后的研究中加以验证。可以预见的是，坡长的增大将显著增加径流量和泥沙量，初始含水量的增大将提前产流时间并增加养分的峰值流失速率[25，26]。同时，土壤中植物根系的生长和土壤生物活性可能形成连通的土壤孔隙结构，导致优先流的出现，这将对模型的模拟精度产生较大影响。泥沙颗粒中通常吸附有大量养分，但本研究的模型并未将泥沙考虑在内，使得模型并不能完整模拟径流过程所带走的养分总量。简言之，该近似解析模型充分考虑了非饱和土壤水分入渗对交换层中养分运移过程的影响，因此，该模型可用于预测干旱和半干旱气候条件下裸露坡耕地的养分流失过程。但是，径流过程的精确计算是进行养分流失模拟的基础，应根据土壤质地、养分类型和雨水中养分浓度选择合适的入渗公式和溶质吸附系数。

3 结 论

本研究以交换层理论为基础，根据黄土区降雨量少，降雨前土壤通常为非饱和状态，其产流需要较长时间的实际情况对降雨过程进行划分，建立了基于坡面径流养分迁移理论的机理模型，并通过 5个降雨强度的模拟降雨试验对模型进行了验证。3组实测重复试验间的标准误差均较小，试验结果具有可靠性。模型验证的结果表明，本文建立的养分流失近似解析模型能够准确描述不同降雨强度下的坡面流和养分流失特征（R2> 0.8，NSE > 0.347）。参数c（入渗率参数）、de（交换层深度）、er（雨滴诱导水分转移速率）均随降雨强度的增大而增大。养分流失模型对交换层的深度de较雨滴诱导水分转移速率er更敏感，de可显著影响可交换溶质的量。因此，在施肥过程中应采取一些措施，如施肥后覆盖坡面土壤、暴雨前避免施肥等，以达到防止土壤贫瘠化并控制农业面源污染的目的。