一维非均质条件下农田土壤水分运移模拟及应用

王佩浩，张 茜，王国芳，张吴平

（1.山西农业大学资源环境学院，山西太谷030801；2.山西农业大学软件学院，山西太谷030801）

土壤水分作为四水转化的中枢、植物生长发育的关键元素，无论是在自然降水条件还是人为的农事活动等的影响下，都有其重要的研究意义。国内外学者在设定均质土壤条件下已对不同初边条件下土壤的入渗、再分配过程进行了大量的理论研究[1-5]。但在实际的田间土壤中，水力特征差异显著的层状土是田间土壤普遍存在的土壤质地结构，会导致土壤水分在均质土壤和层状土壤中的运移具有很大的差异，进而导致模拟结果与真实的土壤水分运移过程存在较大差异，影响了土壤水分运移模拟模型在SPAC 系统和农业精准灌溉中的应用[6]。为此，国内外有学者采用HYDRUS 系列模型以及SWMS 等经典的模拟软件模拟研究了非均质多层土壤中水分的入渗、再分配和排水等过程，得到了理想的模拟结果[7-8]。然而，在实际的土壤水分定量化运移过程中，需要与特定的研究目标结合起来，如研究犁底层特征对大气降雨入渗和土壤再分配的影响，既需要把犁底层的存在对土壤水分运移过程实时地反映出来，又需要采用不同的关系模式。这种情景下，采用已存在的经典模型限制较多，不能灵活地改变不同变量及参数间的关系，限制了土壤水动力学模型在农业生产实践中的应用。

本研究采用Richards 方程的基质势形式描述了田间非均质土壤水分的运移过程，并采用Visual Studio 2017 面向对象的程序语言利用数值差分法求解了多层基质势方程的数值解，在设定犁底层存在的条件下，模拟了土壤水分的运移过程，揭示出多层质地下农田土壤水分运移的规律特征，以期为研究区域农田土壤水分入渗变化提供依据。

1 材料和方法

1.1 研究区选取

大田土壤水分入渗试验于2018 年8 月23 日至9 月8 日在太谷县侯城乡孟家庄村东北部试验地进行。试验地土壤类型为褐土性土，地理位置为东经112°35′，北纬37°26′。

1.2 模型初始值测定

试验地土壤数据采集深度180 cm，每10 cm 为一层，分层用环刀取适量样品进行室内试验测定。测试指标包括土壤粒径组成、容重、饱和含水量和分层土壤含水量（表1），其测定方法依据《土壤物理实验指导》[9]。气象数据取自太谷县气象局。

表1 各层土壤粒径组成、容重及饱和含水量

1.3 模型原理与方法

1.3.1 土壤水分运动方程 本研究忽略土壤中的温度与气相对土壤水分的影响，借助一维饱和-非饱和Richards 控制方程描述多层质地下土壤水分运移[10]。

1.3.2 边界条件的设定

1.3.2.1 初始条件 本研究设定时段初的土壤基质势剖面为一维垂直土壤水分入渗的初始条件。

1.3.2.2 边界条件 考虑到研究区域实际的降雨与蒸发过程，设定上边界条件为二类边界条件，在不同时段设定相应的土壤入渗速率，即在降雨时设定地表入渗速率，蒸发时设定蒸发速率[11]。

式中，f（t）在降雨时刻表示土壤入渗速率，蒸发时刻表示蒸发速率。

下边界条件设定为常数边界条件。

1.3.3 土壤水分特征曲线 通过粒径组成并借助土壤传递函数得到本研究中各层土壤采用Van-Genuchten 模型表示的水分特征曲线初值[12]。

1.4 数值模型参数敏感性分析与确定

采用Visual Studio 2017 面向对象的程序语言，利用数值差分法求解多层基质势方程的数值解。

1.4.1 敏感性参数的确定 将数值模型模拟结果与实测数据进行对比分析，通过对表征模型精度的主要参数进行敏感性分析，确定参数的调整范围。土壤质地状况的差异会引起土体孔隙状况的不同，导致参数α 和n 发生变化，进而影响土壤水分特征曲线的变化[13-16]。因此，本研究确定的主要参数为α和n。

1.4.2 参数优化 将模型推求的各项土壤物理参数初值代入数值模型中进行模拟，不断调整参数α和n 的大小，直至模型输出含水量值与实测含水量（8 月23—30 日）达到最优拟合。参数最优调整依次为：0～10 cm 范围n 值不变、α 值增加10%；10～50 cm 范围n 值与α 值同时减少10%；50～120 cm范围n 值与α 值同时增加10%；120～140 cm 范围n 值与α 值同时增加30%；140～150 cm 范围n 值与α值同时增加20%；150～170 cm 范围n 值与α值同时增加30%；170～180 cm 范围n 值与α值同时增加20%。参数α和n 的初值和优化值列于表2。

表2 关键参数初值与优化值

1.4.3 模型精度的验证 利用参数优化后的模型进行土壤水分运动变化模拟，模拟的土壤含水量值与实测值（9 月1—7 日）通过均方根误差[17]（RMSE）和纳什系数[18]（NSE）2 个指标进行验证。二者统计分析结果如图1 所示。

图1 结果显示，本次模拟结果表层0～10 cm模拟效果误差大，其余各层的模拟效果都很好，表明模型模拟可对实测值进行预测，进而准确地反映出农田土壤水分的运移变化情况。

2 模型应用与分析

土壤的黏粒含量和容重是影响土壤水分入渗的重要因素，二者的变化会使土壤孔隙发生变化，进而引起参数α和n 的变化[13，19-20]。为了更好地检验模型关键参数α和n 对水分运移模拟的影响以及模型的适用性，本研究进行了如下模拟设定（图2）。

据太谷县多年降雨数据，设定模拟初期降雨强度为0.208 cm/h，降水时长24 h，后期设置为自然条件下的蒸发速率（此处输入的是实测的蒸发速率数据）。下边界条件设定为常数边界条件，土壤水分自由下渗。代入数值模型进行模拟，模拟结果如图3所示。由图3 可知，总体来看，除第1 天外，在降水后的一段时间内，0～70 cm 土层土壤含水量随时间增加基本呈逐渐降低趋势，70～120 cm 土层含水量随时间增加逐渐增大，120～180 cm 土层土壤含水量随时间增加逐渐降低。120～180 cm 土层，2 种类型土壤各层含水量基本随时间增加而降低，且变化趋势基本相同。考虑该段土壤质地类型主要为粉黏壤和黏壤土，粉黏粒含量较大，表明土壤深度超过120 cm 土层处含水量变化受粒径组成影响持水稳定，且受上层土壤影响较小，随时间增加，下层土壤水分在重力和土壤本身吸力作用下入渗，因此，土壤含水量变化基本稳定。

图3-a 与图3-b 相比，0～50 cm 处含水量值图3-a 各层含水量均低于图3-b，该范围内αa＞αb、na＞nb，使得图3-a 的土壤进气值小于图3-b，因此，图3-a 的土壤含水量较图3-b 低。图3-b 在30 cm处由于犁底层存在的影响，土壤含水量曲线较图3-a 出现明显的弯折，且较无犁底层的土壤含水量曲线分布更为紧密，此处na＜nb，因此，同一点处土水势图3-b 较图3-a 降低，土壤吸力增大，受αa＞αb的影响，含水量图3-b 高于图3-a。50～70 cm 处模拟图3-a 与图3-b 的各层含水量均随时间增加而降低，但图3-a 的变化幅度明显大于图3-b，图3-b 的变幅较小，这也表明具有较小空隙和较多黏粒组成的犁底层能够很好地降低土壤水分入渗速率和流失速率，具有良好的吸持水能力。

3 结论与讨论

本研究首先通过对研究区域进行实地采集与测量得到所需各项土壤物理数据，结合数值理论模型，搭建基于Visual Studio 的一维土壤水分运移模型。其次，借助实测数据对数值模型进行了参数敏感性分析及模拟精度验证。最后，为更好地检验模型关键参数α 和n 对水分运移模拟的影响以及模型本身的适用性，通过分设不同性质土层进行模拟。研究结果表明，模型参数敏感性主要受土壤水力特征参数α 和n 的影响，参数值越大，土壤进气值越低，含水量越小；n 值减小，土水势降低，土壤吸力增大，含水量增加；在调整数值模型精度时，α 和n 的单个误差一般为10%～30%；模拟值和测量值的变化过程比较吻合与一致，表明本研究构建的基质势方程能够实现对田间一维多层质地土壤水分运移过程的准确模拟。

本研究土壤水分运移模型可以较好地反映土壤水分变化特征，具有较高的应用价值。与其他模型相比，数据输入与操作相对简便，整体过程可控，能够满足对特定问题的专门研究。模型应用也较好地体现了土壤水力参数变化对模型精度的适用性。本研究后期将考虑在不同时空状态下土壤水分的运移变化，同时结合植被生长、根系吸水等地上植物情况，完善土壤水分运移模型，以期为更好地研究区域农田土壤水分入渗状况提供相应的参考依据。