刘 昭，徐燕星，聂小飞，胡 皓，郑海金

(1.江西省土壤侵蚀与防治重点实验室，江西 南昌 330029；2.江西省水土保持科学研究院，江西 南昌 330029； 3.江西水利职业学院，江西 南昌 330013)

坡耕地作为南方红壤区的重要农业生产资源被高强度开发利用，加之南方红壤区易出现集中强降雨情况，因此极易发生水土流失，造成大量氮素流失，而该区域特殊的地形和土壤理化性质，使得驱动土壤氮素运移的地表-土壤水文过程和土壤侵蚀过程具有十分复杂的特征。地表径流和土壤水运动过程相互作用特征明显，氮素迁移随地表径流、泥沙及土壤水文过程在时间和空间上呈动态变化，是地表径流、土壤水运动和泥沙输移等多种过程耦合作用的结果。目前国内外对坡地氮素流失的研究主要集中在径流流失方面，对地表径流氮素迁移过程的影响因素和相关的调控措施进行了深入的分析[1-2]，但对氮素流失与水文泥沙过程(尤其是土壤水文过程)的相互关系尚未深入研究。因此，构建氮素随水文泥沙过程迁移的数学模型，将有助于防治坡地氮素流失和减少水体污染。

土壤水流、地表径流和泥沙输移具有明确的动力学机制，目前已有大量模型考虑了这些动力学过程。在土壤水流运动方面，以达西定律和质量守恒定律推导的Richards方程具有坚实的物理基础，对于水分通量计算具有较大的优势[3]，在各种水文模型中得到了广泛应用，如HYDRUS[4]、SWAP[5]等，但是这些模型以土壤水流模拟为主，很少详细考虑地表水沙动力学过程。由经典的Navier-Stokes流体运动方程，可以得到在坡面地表径流模拟中广泛使用的二维浅水方程、一维Saint Venant方程，以及更为简化的运动波近似方程；对于泥沙输移模拟，一般结合地表水动力学模拟结果，采用泥沙输移方程求解。目前，也有许多基于这类地表水沙动力学方程的模型，如国外的WEPP模型[6]、EROSEM模型[7]等，国内的雷廷武模型[8]、龙满生模型[9]等。但是上述这类地表水沙模型往往忽略土壤水分过程或以经验模型进行描述。因此，坡地土壤水沙动力学模型构建的关键点之一是如何耦合土壤水流、地表径流和泥沙输移的动力学过程。

本研究利用江西水土保持生态科技园大型土壤水分渗漏装置的地表径流、泥沙、壤中流及其中的氮素含量实测数据，构建氮素输出与水文泥沙过程的数学关系，耦合土壤水分运移模型HYDRUS-2D和坡面土壤侵蚀模型WEPP，构建可以动态模拟水文泥沙过程的水沙动力学模型，最终形成一套同时可定量描述地表径流、泥沙、壤中流、氮素输出过程的综合模型，以期为红壤坡地水土流失预测和农业面源污染防治提供计算工具。

1 材料与方法

1.1 试验装置及数据采集

试验场地设置在江西水土保持生态科技园，园区气候、地貌、土壤等自然条件在南方红壤区均具有代表性。园区属于亚热带季风气候区，多年平均降水量1 469 mm，多年平均气温16.7 ℃；地貌为浅丘岗地，坡度5°～25°，海拔30～100 m；土壤主要为第四纪红黏土发育的红壤，呈酸性至微酸性。试验装置为大型土壤水分渗漏装置(裸地)，如图1所示。小区坡度为14°，水平投影长15 m、宽5 m。试验装置的周围及底板用钢筋混凝土浇筑，地板上设砂砾反滤层，坡脚修筑梯形钢筋混凝土挡土墙，挡土墙高出地表30 cm以阻止水分进出小区，从而形成一个封闭排水式的土壤入渗装置。小区土壤为第四纪红黏土发育的红壤，土层深105 cm，剖面为三层，其中0～30 cm为Ah层(淋溶层)，60 cm以下为Bsv层(淀积层)，30～60 cm为Bs层(过渡层)。坡底设置地表径流和30、60、105 cm深壤中流的收集装置。

图1 土壤水分渗漏装置示意

2015年5月22日在小区施尿素300 kg/hm2，施肥后进行逐场次降雨的径流、泥沙和氮素输出试验观测。降雨量采用试验区内设置的虹吸式自记雨量计和全自动气象采集系统进行监测，可以获得每次降雨的降雨量和降雨历时；地表径流量通过地表径流池池壁水尺读数计算；不同层次壤中流量通过自记水位计数据读出；产流结束后，从地表径流池取样并分析径流中泥沙和总氮含量，从不同层次壤中流收集装置中取等量水样充分混匀后分析其总氮含量，泥沙量测定采用称重法进行，总氮的测定采用碱性过硫酸钾消解紫外分光光度法。至2016年4月22日，共收集到21场次降雨的地表径流数据，由于地表径流产流历时短，因此记录数据为每场降雨发生的累积地表径流量；而壤中流持续时间很长(尤其是底层壤中流)，故壤中流数据为日记录数据，共337组。因人工采样受限于时间和设备等因素，故泥沙量、地表氮素输出量和壤中氮素输出量共有10、12、14场次降雨的数据。

1.2 模拟评价相关说明

本研究将2015年5月22日至2016年1月31日间的数据用于经验模型构建及机理性模型参数校正，这段时期称为模型校正期，期间地表径流量、壤中流量、泥沙量、地表氮素输出量和壤中氮素输出量数据分别有16、255、6、8和9组；而将2016年2月1日至2016年4月22日间的数据用于模型验证，这段时期称为模型验证期，期间地表径流量、壤中流量、泥沙量、地表氮素输出量和壤中氮素输出量分别有5、82、4、4和5组样本数据。

研究中为了判别模拟值与观测值之间的拟合精度，需要引入统计学指标。常用于判别拟合精度的统计学指标有决定系数R2(R2∈[0,1])、均方根误差RMSE(RMSE∈[0,∞))、偏差BIAS(BIAS∈[0,∞))、一致性指标IA(IA∈[0,1])等。其中：BIAS表示的是平均误差与观测均值之比，与R2、IA都是相对指标，无量纲；RMSE的量纲与统计对象相同。各指标计算公式为

(1)

(2)

(3)

(4)

(5)

(6)

2 模型构建

2.1 氮素输出经验模型

如图2所示，分析模型校正期有效观测数据发现，壤中氮素输出量与壤中流量显著相关(相关系数为0.95)。在实际中，若无壤中流产生，则氮素不会从土壤中输出，因此本研究采用无截距的二次函数表示壤中氮素输出量与壤中流量的关系，即

图2 壤中氮素输出经验模型

(7)

式中：Ng为壤中氮素输出量，g；Qg为壤中流量，mm。

相关性分析表明，地表氮素输出量与泥沙量的相关系数达到0.92，而地表氮素输出量与地表径流量的相关系数为0.78，可见地表氮素输出量与土壤侵蚀关系更为密切。通过分析，本研究采用幂函数表示地表氮素输出量与泥沙量的关系，即

(8)

式中：Ns为地表氮素输出量，g；Me为泥沙量，kg。

图3 地表氮素输出经验模型

由表1可知，式(7)和式(8)均具有较高的模拟精度。

表1 氮素输出经验模型评价

2.2 壤中流和地表径流耦合下氮素输出模型构建

Richards方程由达西定律和质量守恒定律推导而来，对于水分通量计算具有较大的优势，可用于描述饱和-非饱和土壤水流动力学过程。本研究假设垂直于坡向上的水流过程是相同的，则土壤水流动力学过程以二维Richards方程表示，即

(9)

式中：t为时间，h为土壤压力水头；K(h)为水力传导度；C(h)为容水度；β为参数，在土壤非饱和时为0，饱和时为1；μs为弹性释水系数；x为坡向水平投影方向坐标；z为距离地表的深度，向下为正。

式(9)的定解条件可根据实际情况进行设置，下边界根据地下水位情况设为定水头边界或者根据不透水层位置设为0通量边界；在降雨条件下，上边界可设为给定通量边界，即降雨量与地表径流量之差。以混合型Richards方程为基础，采取有限单元法进行时空离散，采取Picard迭代方法进行数值求解，能够模拟土壤水分、压力和边界通量动态过程，具有较高的模拟精度，在农田水利、水文水资源及地下水污染防治方面得到了广泛应用，因此本研究选取HYDRUS-2D模型对式(9)进行求解。

对于径流泥沙过程的模拟，可参考美国农业部开发的WEPP模型。其单次降雨模块以一维运动波方程模拟地表径流，描述为

(10)

式中：l为某点沿下坡方向的距离，h为地表水深，r为净雨强，i为入渗速率，θ为坡面与水平面夹角，q为地表径流通量。

泥沙输移方程可表示为

(11)

式中：G为输沙量；DL为从相邻坡面流入的泥沙量；DF为水流对沟道的剥蚀量或水流中泥沙的沉积量。对于裸土而言，当确定了土壤相关参数、土壤含水量和气象因素后，即可借助WEPP模型对式(10)和(11)进行求解。

为保证上边界通量和土壤水分计算的一致性，本研究所采取的耦合方法如图4所示，即：将HYDRUS-2D模型土壤水分模拟结果作为WEPP模型的输入条件，然后以WEPP模型计算的径流量分割降雨量(扣除蒸发量后)，从而得到HYDRUS-2D模型的上边界通量；以FORTRAN 95程序编写WEPP模型和HYDRUS-2D模型的输入输出文件交互接口，并控制程序运行，模拟出壤中流量和侵蚀泥沙量后，根据式(7)和式(8)即可得到壤中和地表氮素的输出量。

图4 综合模型框架

3 参数校正和模型检验

采用图4所示模型进行模拟，还需要确定机理性模型的参数。对于模型参数初值，一部分采用实测法获取，其余则通过模型自带的估算方法(如HYDRUS-2D中的土壤水动力参数、WEPP中的土壤侵蚀参数等)获取。其中实测参数见表2。

表2 实测的模型参数

研究中发现将上述参数初值直接用于模拟，会使得地表径流、壤中流和泥沙的模拟值与实测值相差较大，因此需要进一步校正参数。通过敏感性分析可知，对地表径流量、壤中流量和泥沙量影响最大的为饱和水力传导度(Ks)，其余参数影响相对较小，因此主要是对各层饱和水力传导度进行了校正。如前所述，模型校正期选为2015年5月22日至2016年1月31日，参数校正前后模拟值与观测值的比较结果见表3，校正后的参数模拟精准度得到了很大程度的提高(R2和IA增加，BIAS和RMSE减少)。以校正后的参数Ks进行模拟，壤中流的模拟精度最高，泥沙量的模拟精度相对较差。值得注意的是，各土层饱和水力传导度实测值与校正值可相差1.4～4.7倍(见表4)，这是导致模拟值与实测值相差较大的根本原因。

进行参数校正后，即可对模型验证期(2016年2月1日至2016年4月22日)的地表径流量、 壤中流量、泥沙量、地表氮素输出量、壤中氮素输出量进行模拟。如图5所示，模拟结果不仅可以较好地匹配地表径流量、壤中流量、泥沙量、地表氮素输出量、壤中氮素输出量等观测值，还可以估计出未能实测到的数据。精度评价指标计算结果(表5)表明，壤中流量和壤中氮素输出量的模拟效果(决定系数达到0.84以上，一致性指标达到0.96以上)要强于地表径流量、泥沙量和地表氮素输出量(决定系数为0.41～0.68，一致性指标为0.58～0.87)。分析模拟结果(图5)发现，在降雨产流期壤中流的氮素输出量显著高于地表径流的氮素输出量，且在未降雨期间氮素随着壤中流持续输出，说明壤中流是红壤坡地氮素的主要流失途径，这与前人研究结论一致[1]。此外，壤中流产流高峰较降雨峰值延后，且一次降雨的壤中流会持续数天，说明土壤起到了一定的缓冲作用。

表3 参数修正前后模拟评价

表4 各土层饱和水力传导度校正前后比较

4 结 论

以实测数据分析为基础，建立了氮素输出经验模型，并将该经验模型与土壤水运动模型(HYDRUS-2D)和坡面土壤侵蚀模型(WEPP)耦合，形成了一套模拟径流-泥沙-壤中流-氮素输出的综合模型。以实测数据对模型参数进行校正并进行了模型检验，校正后的参数对径流量、泥沙量、壤中流量、氮素输出量均具有较好的模拟效果(决定系数在0.41以上)，其中壤中流量和壤中氮素输出量模拟效果更好， 决定系数分别达0.88和0.84。分析模拟结果表明，在降雨产流期壤中氮素输出量显著高于地表氮素输出量，且在未降雨期间氮素随着壤中流持续输出，进一步印证了壤中流是红壤坡地氮素的主要流失途径。此外，壤中流产流高峰较降雨峰值延后，且一次降雨的壤中流会持续数天，说明土壤起到了一定的缓冲作用。

图5 模型验证期动态模拟结果

数据类别样本数量评价指标R2RMSEBIASIA地表径流量50.5140.9750.6130.585泥沙量40.6850.0060.3870.869地表氮素输出量40.4100.0710.5720.648壤中流量820.8840.7650.2690.972壤中氮素输出量50.8400.0640.1440.962

本研究已基本建立了具有一定物理意义、适合红壤坡地产流的数学模型，并可定量描述壤中流和地表径流耦合下的土壤氮素迁移输出过程。这将为红壤坡地水土流失预测和农业面源污染防治提供计算工具。受限于时间，模型研究中还需要进一步加强和完善的有以下方面：一是今后应加强观测，减少漏测现象，提高观测精度(如土壤水分、泥沙量等)，因为观测数据的质量在很大程度上决定了模型及参数的准确程度；二是模型自身缺陷还需修补，例如HYDRUS-2D模型遇强降雨时不收敛的问题，可以通过修改其源代码数值格式进行修补。