东北黑土区是我国重要的粮食主产区，也是我国重要的商品粮基地，被视为中国粮食安全的“压舱石”。该地区纬度较高，雨热同期，属寒温带大陆性季风气候。黑土具有有机质含量高、质地黏重、入渗能力较弱、易产生地表径流的特征。在近百年的开垦过程中，由于缺乏保护，黑土区水土流失导致的土壤退化问题越来越明显，并得到广泛关注。地表径流是坡面土壤侵蚀发生的动力基础。径流是由降雨或灌溉水到达地表之后超过土壤入渗能力的水分形成的。因此，土壤入渗是影响侵蚀过程的重要环节。研究不同侵蚀程度黑土的入渗特征对治理坡面土壤侵蚀具有重要意义。

不同类型的土壤入渗过程差异很大。每种入渗模型都不可能适用于所有土壤。依据土壤的特性选择适用的入渗模型并加以改进，是目前主流的土壤入渗模拟方法。HYDRUS-1D是由美国国家盐改中心(US Salinity Laboratory)于1991年开发的商业化软件。该软件是一种用于模拟变饱和带多孔介质中的水分、热量和溶质运移的数值模型，能够较好地模拟水分、溶质和能量在土壤中的时空分布变化及运移规律^[1]。郭雯等^[2]在内蒙古高原西部腰坝绿洲利用该模型模拟土壤水分入渗过程，模拟结果与实测含水率拟合较好，模型平均R_E、R_MSE值分别为6.1%、0.015 cm³/cm³，误差较小，精度较高；李琦等^[3]以华北平原区典型冬小麦农田土壤为研究对象, 结合野外观测与室内土柱试验，借助该软件建立水盐运移模型，模拟结果的纳什效率系数平均值为0.826，变异系数为0.056 0，能较好且稳定地模拟土壤内部水分的运移过程；王国帅等^[4]利用该模型模拟沙丘、沙丘—荒地交界和荒地土壤水分和盐分运移动态，R_MSE为0.01~0.03 cm³/cm³，R²为0.85~0.92，其结果能较好地反映出土壤水分的动态变化。我国HYDRUS-1D模型的应用研究多集中于西北干旱地区，该模型是否能模拟黑土入渗需要进一步研究。

笔者将基于室内土柱试验，观测典型黑土区不同侵蚀程度黑土土壤入渗过程，并利用HYDRUS-1D模型模拟，通过比较模型模拟结果与实测数据，评价该模型对东北黑土区不同侵蚀程度土壤的模拟效果，为黑土土壤入渗模拟和土壤侵蚀预报提供基础。这对黑土土壤水蚀过程模拟和水土流失治理有重要的实际意义。

室内土柱试验于2020年在黑龙江省黑河市嫩江县鹤山农场境内的北京师范大学九三水土保持试验站(E 125°16′~125°21′，N 48°59′~49°03′)完成。供试土壤选择轻度、中度及重度3种不同侵蚀程度的土壤(表 1)，测定其一维垂直入渗过程，每种土壤类型设置3个重复。不同侵蚀程度黑土根据表层的黑土层(A层)厚度和其下淀积层(B层)厚度确定：轻度侵蚀的土壤A、B层分别为30和25 cm，中度侵蚀的土壤A、B层分别为30和10 cm，重度侵蚀的土壤A、B层分别为0和15 cm(表 2)。

表 1(Tab. 1)

表 1 3种侵蚀程度土壤密度及机械组成 Tab. 1 Soil bulk densities and mechanical compositions at three eroded degrees 侵蚀程度 Eroded degree 土层 Soil horizon/cm 土粒密度 Bulk density/(g·cm-3) 砂粒 Sand content/% 粉粒 Silt content/% 黏粒 Clay content/% 轻度Slight 0~10 1.06 12.45 45.56 42.00 10~20 1.28 11.27 47.69 41.04 20~30 1.24 9.47 48.99 41.54 30~40 1.49 10.23 48.17 41.60 40~50 1.52 11.75 47.28 40.96 中度Moderate 0~10 1.18 19.47 46.94 33.59 10~20 1.39 28.50 48.97 22.53 20~30 1.47 29.46 48.24 22.30 30~40 1.46 31.40 38.07 30.53 40~50 1.45 30.08 38.13 31.79 重度Severe 0~10 1.43 63.94 15.70 20.36 10~20 1.53 66.34 14.62 19.04 20~30 1.81 64.33 10.39 25.28 30~40 2.07 68.98 12.79 18.23 40~50 1.71 74.67 11.27 14.06

表 1 3种侵蚀程度土壤密度及机械组成 Tab. 1 Soil bulk densities and mechanical compositions at three eroded degrees

表 2(Tab. 2)

表 2 3种侵蚀程度土壤取土位置和土壤剖面 Tab. 2 Soil sampling position and soil profile at three eroded degrees 侵蚀程度 Eroded degree 取土位置 Soil location 地貌部位 Landform position 黑土(A)层深度 Black soil(A)depth/cm 淀积(B)层深度 Deposition(B) layer depth/cm 母质出现深度 Appearing parent material depth/cm 母质类型 Parent material type 轻度 Slight E125°18′21″N49°01′58″ 坡中Middle slope 0~30 30~55 55 黄黏土 Yellow clay 中度 Moderate E125°19′40″N49°00′39″ 坡中上 Middle and upper slope 0~30 30~40 40 黄黏土 Yellow clay 重度 Severe E125°19′24″N49°00′42″ 坡中上 Middle and upper slope 0 0~15 15 黄砂 Yellow sand

表 2 3种侵蚀程度土壤取土位置和土壤剖面 Tab. 2 Soil sampling position and soil profile at three eroded degrees

试验由土柱系统和供水装置组成。

土柱系统：土柱装土50 cm高，直径28.4 cm，为透明有机玻璃圆柱。每次试验准备3个土柱，分别按照轻、中、重3种侵蚀程度的原状土壤分层，每层10 cm分层填装进土柱，装土密度与原状土密度一致，层间用毛刷打毛。将土柱挂上铁架，每个土柱每层中间位置(即5、15、25、35和45 cm处)安装水分温度传感器探头，探头长度为15 cm，数据采集器每5 min输出相应位置的土壤含水量数据。为确保试验的准确性，再设置2次重复试验。

供水装置：马氏瓶，保证土面上方有3 cm的稳定水头。

试验开始，马氏瓶中水位每下降1 cm手动记录1次历时。数据采集器每隔1 min记录土柱质量的变化，每隔5 min记录剖面含水量变化。若连续5个时刻记录点的入渗速率基本相同，将第1个时刻记作达到稳渗的时刻，并把连续5个记录点的平均入渗速率做为稳渗速率。同时记录湿润锋运移情况，湿润峰到达土柱底部后停止马氏瓶供水。用塑料薄膜密封土柱上表面防止蒸发，待所有土柱排水口不再有水流出时，撤下塑料薄膜，试验结束。

综上，本研究进行3次试验，共9组数据，分别计算：1)累积入渗量与入渗速率。土柱试验过程中需要手动记录入渗过程中马氏瓶中水位每下降1 cm所需要的时间。测量马氏瓶水位每下降1 cm的水流通量，根据记录表中的时间间隔计算入渗过程中的累积入渗量与入渗速率。

式中：I_n为每段时间的累积入渗量，cm；Δl为马氏瓶水位差；Q₀为马氏瓶每cm的水流通量，mL；S为土柱的横截面积，cm²；I为总累积入渗量，cm；f为入渗速率，cm/h；t为时间，h。

2) 土壤含水量。数采每5 min输出相应位置的土壤含水量数据。

3) 湿润锋运移。土柱上贴有刻度尺，根据湿润锋运移的速率选择记录间隔。在试验刚开始时入渗速率较快，每3 min记录1次湿润锋位置，随着试验的进行湿润锋移动速率逐渐减慢，记录时间间隔增大，最大时间间隔不超过2 h。将记录下的刻度转化为入渗深度，得到湿润锋运移情况。

笔者使用软件的模拟水流程序模拟土壤水分入渗。该模型仅考虑一维垂向运移时的土壤水分运动，并采用Richards方程来描述一维非饱和土壤水流的运动^[5]。HYDRUS-1D模型中水分特征曲线模型有Van Genuchten、Modified Van Genuchten、Brooks-Corey和Kosugi等。笔者采用运用最为广泛的、也更适合东北黑土^[6]的Van Genuchten模型来拟合土壤水力性质。HYDRUS-1D中VG模型表达式为：

式中：θ(h)为土壤含水量随土壤水势的变化函数；h为土壤的有压水头，cm；θ_r和θ_s分别为残余含水量和饱和含水量, cm³/cm³；K(h)为非饱和导水率，cm/h；K_s为饱和导水率，cm/h；α为进气吸力值的倒数；n为孔隙度连通性参数，m为经验参数；l为形状参数，l=0.5；S_e为有效水分含量。

入渗模拟过程中需要设置的指标有剖面形态、时间步长、土壤水力参数和边界条件等。其中土壤水力参数是进行模拟的主要指标，包括：残余含水量θ_r、饱和含水量θ_s、进气吸力值α、孔隙度连通性参数n、饱和导水率K_s和形状参数l等，这些参数主要通过土壤水分特征曲线模拟来获得。本实验中θ_r、θ_s和K_s均有实测数据，α、n和l 3个参数依靠模型内置函数模拟，共设置5种模拟情况：1)θ_r、K_s、θ_s均为模拟值；2)θ_r、K_s、θ_s均为实测值；3)θ_r、K_s为模拟值，θ_s为实测值；4)K_s为模拟值，θ_r、θ_s为实测值；5)θ_r为模拟值，K_s、θ_s为实测值。采用Excel 2010和SPSS 24.0软件对数据进行分析，用Pearson法对累积入渗量、入渗速率的模拟值与实测值进行相关分析，利用Origin 2010软件作图。

由图 1可见，重、中、轻度侵蚀土壤达到稳渗的时间分别为8、26和36 h，稳渗时的入渗速率为0.72、0.25和0.16 cm/h，完成整个入渗的历时分别为10.2、31.5和50 h。针对一次完整的入渗过程，轻度侵蚀土壤每层累积入渗量分别为7.3、4、3.9、3.4和3.6 cm；中度侵蚀土壤每层累积入渗量分别为7、3.6、3.4、4.5和2.5 cm；重度侵蚀土壤每层累积入渗量分别为5、3.2、2.6、2和3.7 cm。通常情况下，黏粒含量越高的土壤对水的吸附能力越强，相同势能下的水分通量会相应变小，因此侵蚀程度越高的土壤在相同的入渗历时内累积入渗量会越大(图 1)。

图 1(Fig. 1)

图 1 不同侵蚀程度黑土土壤入渗速率及累积入渗量 Fig. 1 Infiltration rates and cumulative infiltrations of black soils with different erosion degrees

湿润锋的运移可以反映出土壤的入渗能力^[7]。轻、中、重3种侵蚀程度的土壤均是越向下层入渗速率越慢，湿润锋运移速率相应减慢(图 2)。轻度侵蚀的土柱在1 h处完成第1层10 cm土壤的入渗，并开始第2层的入渗，5 h开始第3层入渗，14.5 h开始第4层入渗，28.7 h开始第5层入渗，50 h的时候完成整个入渗过程。每层入渗耗时分别为1、4、9.5、14.2和21.3 h。中度侵蚀的土柱每层入渗耗时分别为1、4.4、7、10.1和14 h；重度侵蚀的土柱每层入渗耗时分别为0.25、0.6、1.7、3.5和5.2 h。轻、中、重3种侵蚀程度湿润锋运移的平均速度分别为：2.52、3.32和9.85 cm/h。

图 2(Fig. 2)

图 2 轻、中、重度侵蚀土壤入渗湿润锋运移变化 Fig. 2 Changes of wetting front migrations of infiltrations in slight, moderate and severe eroded soils

将入渗过程中土壤含水量最大并不再发生明显变化的数值作为土壤饱和含水量。3种土壤在入渗过程中的土壤水分含量随时间变化的趋势存在一定的差异。土壤饱和含水率随着侵蚀程度的增大而减少，轻度侵蚀土壤比中度侵蚀土壤饱和含水率平均高4.9%左右，比重度侵蚀土壤饱和含水率平均高9.4%左右。每种土壤的饱和含水率随土壤层数的加深而减少，3种土壤按侵蚀强度由弱到强，其表层土壤饱和含水率比最下层分别大7.8%、8.0%和12.8%(图 3)。

图 3(Fig. 3)

图 3 轻、中、重度侵蚀土壤10~50 cm分层土壤含水率随时间变化 Fig. 3 Variation of soil moisture contents with time in 10-50 cm divided layers of slight, moderate and severe eroded soils

通过HYDRUS-1D内置函数所得到的模拟值如表 3所示，其他模拟情况均根据表中参数进行校正设置。

表 3(Tab. 3)

表 3 3种侵蚀程度土壤土柱所有模拟参数数值 Tab. 3 All simulated parameters of soil column with three eroded degrees 侵蚀程度Eroded degree 土层Soil horizon/cm θr θs α/(cm-1) n Ks/(cm·h-1) 轻度Slight 0~10 0.102 8 0.566 0 0.014 3 1.389 4 2.496 7 10~20 0.096 8 0.499 7 0.011 7 1.419 6 0.721 7 20~30 0.098 6 0.515 1 0.012 1 1.413 7 0.902 9 30~40 0.091 0 0.442 7 0.011 2 1.397 5 0.211 7 40~50 0.089 0 0.430 9 0.011 2 1.392 0 0.179 2 中度Moderate 0~10 0.091 3 0.503 9 0.009 5 1.490 8 1.366 2 10~20 0.068 9 0.411 6 0.006 9 1.586 7 0.539 6 20~30 0.066 1 0.393 4 0.007 5 1.558 4 0.375 4 30~40 0.077 8 0.417 8 0.010 8 1.453 1 0.298 3 40~50 0.080 0 0.424 2 0.010 9 1.449 5 0.313 3 重度Severe 0~10 0.062 9 0.419 5 0.022 8 1.422 8 1.457 5 10~20 0.058 1 0.391 3 0.025 8 1.397 5 1.096 3 20~30 0.050 6 0.313 6 0.014 3 1.311 4 0.090 0 30~40 0.044 9 0.258 0 0.038 9 1.250 8 0.115 8 40~50 0.047 4 0.337 6 0.036 3 1.398 1 0.969 6 注：θr：残余含水量，cm3/cm3；θs：饱和含水量，cm3/cm3；α：进气吸力值的倒数；n：孔隙度连通性参数；Ks：饱和导水率，cm/h。下同。Notes: θr: Residue water, cm3/cm3; θs: saturation capacity, cm3/cm3; α: the inverse of air-entry value; n: porosity connectivity parameters; Ks: saturated hydraulic conductivity, cm/h.The same below.

表 3 3种侵蚀程度土壤土柱所有模拟参数数值 Tab. 3 All simulated parameters of soil column with three eroded degrees

根据入渗过程拟合可以看出，总体上轻度侵蚀土壤的累积入渗量被高估，重度侵蚀土壤的累积入渗量被低估，中度侵蚀土壤的拟合情况相对最好，模拟值与实测值分布于1∶1线两侧，且较为集中，所以HYDRUS-1D对于相对极端情况的黑土拟合情况较差(图 4)。

图 4(Fig. 4)

S：模拟值；M：实测值。下同。 S: Simulated value. M: Measured value. The same below. 图 4 轻度、中度和重度侵蚀土壤累积入渗量模拟值与实测值相关性分析 Fig. 4 Correlation analysis between simulated and measured cumulative infiltration in slight, moderatel and severe eroded soils

入渗速率由土壤累积入渗量计算得到，所以其拟合结果的特征与累积入渗量的规律相似(图 5)。除重度侵蚀土壤能够看出实际入渗速率明显高于模拟结果，其余土壤实际入渗与拟合结果的差异并不明显。轻、中、重度3种土壤的稳定入渗率模拟值与实测值也存在差异(表 4)，轻中度侵蚀土壤的稳定入渗率被高估，模拟值分别为实测值的2.5倍和1.6倍，而重度侵蚀土壤的稳渗率被低估，模拟值为实测值的0.4倍。

图 5(Fig. 5)

图 5 轻度、中度和重度侵蚀黑土入渗速率模拟值与实测值相关性分析 Fig. 5 Correlation analysis between simulated and measured values of infiltration rates of slight, moderate and severe eroded black soil

表 4(Tab. 4)

表 4 3种侵蚀程度土壤稳定入渗率及模拟值 Tab. 4 Soil stable infiltration rate and simulated values at three eroded degrees 侵蚀程度 Eroded degree 实测稳定入渗率 Measured stable infiltration rate/(cm·h-1) 5种模拟情况下的稳定入渗率 Stable infiltration rate under five simulated conditions/(cm·h-1) S：θr、Ks、θs M：θr、Ks、θs S：θr、Ks；M：θs S：Ks；M：θr、θs S：θr；M：Ks、θs 轻度Slight 0.16 0.398 0.005 7 0.398 0.398 0.007 3 中度Moderate 0.25 0.398 — 0.446 0.446 — 重度Severe 0.72 0.216 0.268 0.202 0.217 0.278 RMSE — 0.367 — 0.383 0.390 — 注：—：模拟过程提前结束，无数值的情况。Notes: —: The simulation process ends early, and there is no numerical case.

表 4 3种侵蚀程度土壤稳定入渗率及模拟值 Tab. 4 Soil stable infiltration rate and simulated values at three eroded degrees

本研究通过室内土柱试验，分别分析不同侵蚀程度黑土的初始入渗速率、稳定入渗率、湿润锋运移速率以及土壤剖面含水量随时间变化特征，可知初始入渗速率主要受供水装置影响，入渗过程刚开始时风干土壤的含水量很小，土壤的基质势梯度很大，土壤的入渗能力较强，入渗速率主要取决于供水速率。但随着入渗过程的不断进行，土壤的基质势梯度会逐渐趋近于0，使入渗速率逐渐保持在一个较小的稳定值，即进入稳定入渗阶段。稳定入渗阶段是指土壤入渗经过瞬变和渐变阶段之后达到的一种入渗速率处于相对稳定的阶段。稳渗速率是反映土壤渗透性能的主要指标^[8]。不同侵蚀程度黑土土壤机械组成结构不同，侵蚀程度强的土壤中砂粒含量越高，土壤孔隙度越大，导水率相应变大，入渗过程加快，稳渗率越高。因此土壤达到稳渗的时间越短，完成入渗的时间也越短，最终的稳定入渗率较高，反之递减。不同侵蚀程度黑土下不同剖面含水量其分布规律与土壤累积入渗量的规律相似，这是由于土壤中黏粒含量较高会吸附更多的水分导致土壤饱和含水率更高。

3种侵蚀程度土壤的土壤含水量的模拟结果与入渗速率和累积入渗量的规律相符，均为轻中、度侵蚀土壤模拟的入渗过程快于实测过程，而重度侵蚀土壤模拟的入渗过程慢于实测过程。这表明HYDRUS-1D模拟入渗过程会低估砂粒含量较高的土壤入渗，高估黏粒含量较高的土壤入渗，黑土较为黏重的质地影响HYDRUS-1D中内置函数对于土壤入渗过程的模拟，导致入渗模型的适用性稍差。在所模拟的6个水力参数中，针对轻度侵蚀地块土壤，K_s是主要的误差来源，对拟合结果的影响最明显：拟合时表层K_s值有所低估，而在深层K_s值被高估；中度侵蚀地块的拟合情况较好，有完整模拟过程的模拟情况1)、3)和4)有相同的趋势；重度侵蚀地块的所有模拟情况均明显少于实际入渗。这是由于重度侵蚀土壤底层密度较大，模拟出的K_s值偏小，实测出的K_s值并没有明显大于其他2种侵蚀程度土壤所导致的。水势梯度是由水头、重力、基质力和毛管力决定的，在实测入渗与模拟入渗设置的前期条件相同的时候，两者的水势梯度是没有太大差别的。所以入渗速率的差异主要是由饱和导水率值的差异造成的，饱和导水率是影响模拟结果最关键的因子。

通过室内土柱实验发现，不同侵蚀程度黑土的入渗性能不同：3种侵蚀程度黑土的入渗速率在入渗开始时均较大，而后逐渐减小，最后趋于稳定，每种土壤上层的入渗速率均大于下层。随着侵蚀程度加剧，土壤入渗速率、累积入渗量及湿润锋运移速度加大，土壤达到稳渗的时间越短，完成入渗的时间也越短，最终的稳定入渗率也较高。

土壤饱和含水率随着侵蚀程度的增大而减少。3种不同侵蚀程度土壤入渗过程的含水量模拟结果与实测结果差异较大，轻中度侵蚀土壤模拟的入渗过程快于实测过程，而重度侵蚀土壤模拟的入渗过程慢于实测过程。HYDRUS-1D模拟入渗过程会低估砂粒含量较高的土壤入渗而高估黏粒含量较高的土壤入渗。

应用HYDRUS-1D模型模拟不同侵蚀程度黑土入渗过程的结果发现，轻、中度侵蚀土壤的入渗能力被高估，而重度侵蚀土壤的入渗能力被低估。其中，中度侵蚀的土壤入渗的拟合结果好于其他2种侵蚀程度土壤。可尝试进一步修正残余含水量和饱和含水量的参数数值，优化拟合结果。