CN114397427A - 一种基于环刀土样上渗过程的土壤水力学特性预测方法 - Google Patents

Abstract

本发明涉及一种基于环刀土样上渗过程的土壤水力学特性预测方法，采用全新设计方法，基于恒定水头边界条件下均质干土一维水分上渗解析解的逆过程，通过记录环刀土样上渗过程中累积入渗量与时间的数据、湿润锋到达土样上表面的时间，测定土样饱和导水率，快速获取土壤水力学特性，克服了现有技术测定土壤水力学特性耗时长、效率低，且不适于田间大面积快速获取等问题，避免了参数多解、不收敛等问题，提高了土壤水力学特性的获取效率，为田间大面积的土壤水力特性调查提供了有效手段。

Description

一种基于环刀土样上渗过程的土壤水力学特性预测方法

技术领域

本发明涉及一种基于环刀土样上渗过程的土壤水力学特性预测方法，属于土壤物理和土壤水文学技术领域。

背景技术

土壤水分特征曲线是土壤孔隙结构的量化和表征，吸渗率是土壤孔隙结构吸收和解吸收液体能力的量度。如何准确、快速和低成本的获取土壤水力学特性参数(土壤水分特征曲线、饱和/非饱和导水率、吸渗率)，是帮助我们理解包气带中土壤水肥运动的关键。环刀采样，室内分析是当前获取完整土壤水力学特性参数最为有效的方法，但传统的测量方法需要使用昂贵的设备(如：压力膜仪、吸力平板仪、沙箱等)，通过缓慢的排水、蒸发过程来获取参数，并且各装置都会受到水势测定范围的影响，耗时费力，成本高昂，无法满足田间大面积快速获取的需求。

水分上渗过程因在短距离入渗中也可以获得详细的过程测量，已经被发展成为一种快速、有效的间接方法来估计环刀土样的水力学特性。Moret-Fernández(2017；2019；2021)等通过对数值算法性能的提升，提出了一系列基于上渗过程的数值反演方法来估计土壤水力学特性，然而该方法受到参数估计的非唯一性和不收敛的影响，极大地限制了上渗方法在田间快速获取土壤水力特性上的应用。解析反演方法可以部分克服这些问题，但由于上渗问题的复杂性，现存的上渗解析模型要么形式过于复杂(Swartzendruber，2002)，亦或是简化了太多的过程细节(Prevedello，2009；Ogunmokun and Wallach,2021)，均难以发展出准确的解析反演方法。并且，申请号为202122414262.0的专利申请文件公开了一种野外便携式土壤水力学特性快速测定装置，其提供了一种在田间详细获取环刀土样水分上渗、稳渗过程的测量手段。但由于当前反演技术在理论方法上的限制，难以满足大面积土壤普查即时和高效的调查需求。

发明内容

本发明所要解决的技术问题是提供一种基于环刀土样上渗过程的土壤水力学特性预测方法，通过环刀土样的水分上渗过程，快速预测土壤水力学特性，有效提高了土壤水力学特性的获取效率。

本发明为了解决上述技术问题采用以下技术方案：本发明设计了一种基于环刀土样上渗过程的土壤水力学特性预测方法，执行如下步骤，实现对环刀土样的土壤水力学特性预测；

步骤A.针对环刀土样对应恒定下边界水头h_p下均质干土一维水分上渗过程，记录获得环刀土样累计入渗量I随累计时长t的变化关系，土壤中湿润锋前进距离到达环刀土样表面时的累计时长t_e，以及测量获得环刀土样的饱和导水率K_s，然后进入步骤B；

步骤B.根据环刀土样累计入渗量I随累计时长t的变化关系，获得累计时长t_e所对应的累计入渗量I_te、以及土壤中湿润锋前进距离z_fte，并获得z_fte与I_te的比值，构成环刀土样土壤湿润区含水量的平均增长量U，然后进入步骤C；

步骤C.根据环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r，结合环刀土样土壤湿润区含水量的平均增长量U，获得环刀土样土壤水分剖面的形状系数a，然后进入步骤D；

步骤D.根据环刀土样土壤水分剖面的形状系数a，获得环刀土样的土壤孔隙尺寸分布指数n。

作为本发明的一种优选技术方案：还包括步骤E如下，执行完步骤D之后，进入步骤E；

步骤E.根据环刀土样的土壤孔隙尺寸分布指数n、环刀土样的饱和导水率K_s，结合环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r、以及环刀土样土壤水分剖面的形状系数a，获得环刀土样的有效饱和导水率K_e，然后进入步骤F；

步骤F.根据环刀土样的有效饱和导水率K_e、以及环刀土样累计入渗量I随累计时长k的变化关系，获得环刀土样的吸渗率s。

作为本发明的一种优选技术方案：所述步骤E中，根据环刀土样的土壤孔隙尺寸分布指数n、环刀土样的饱和导水率K_s，结合环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r、以及环刀土样土壤水分剖面的形状系数a，按如下公式：

获得环刀土样的有效饱和导水率K_e，其中，m＝3n+2。

作为本发明的一种优选技术方案：所述步骤F中，根据环刀土样的有效饱和导水率K_e、以及环刀土样累计入渗量I随累计时长t的变化关系，按如下公式：

执行拟合操作，获得环刀土样的吸渗率s。

作为本发明的一种优选技术方案：还包括步骤G如下，执行完步骤F之后，进入步骤G；

步骤G.根据环刀土样的吸渗率s、环刀土样的土壤孔隙尺寸分布指数n、环刀土样的土壤饱和含水量θ_s、环刀土样的土壤初始含水量θ_i、环刀土样的饱和导水率K_s、以及恒定下边界水头h_p，按如下公式：

获得环刀土样的土壤进气吸力值h_d，其中，m＝3n+2。

作为本发明的一种优选技术方案：所述步骤C中，根据环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r，结合环刀土样土壤湿润区含水量的平均增长量U，按如下公式：

获得环刀土样土壤水分剖面的形状系数a。

作为本发明的一种优选技术方案：所述步骤D中，根据环刀土样土壤水分剖面的形状系数a，结合如下公式：

应用不动点迭代算法，获得环刀土样的土壤孔隙尺寸分布指数n。

本发明所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，采用以上技术方案与现有技术相比，具有以下技术效果：

(1)本发明所设计基于环刀土样上渗过程的土壤水力学特性预测方法，采用全新设计方法，基于恒定水头边界条件下均质干土一维水分上渗解析解的逆过程，通过记录环刀土样上渗过程中累积入渗量与时间的数据、湿润锋到达土样上表面的时间、以及测定土样饱和导水率，快速获取土壤的水力学特性，克服了现有技术测定土壤水力学特性耗时长、效率低，且不适于田间大面积快速获取等问题，避免了参数多解、不收敛等问题，有效提高了土壤水力学特性的获取效率，为田间大面积的土壤水力特性调查提供了有效手段。

附图说明

图1是本发明所设计基于环刀土样上渗过程的土壤水力学特性预测方法的流程图；

图2a、图2b、图2c是在三种土壤中上渗解析模型模拟的累积入渗量曲线与数值模拟曲线的对比图，其中，图2a为黏土的累计入渗量对比图，图2b为粉沙土的累积入渗量对比图，图2c为砂壤土的累积入渗量对比图。

具体实施方式

下面结合说明书附图对本发明的具体实施方式作进一步详细的说明。

本发明所设计一种基于环刀土样上渗过程的土壤水力学特性预测方法中，首先分析理查德水流控制方程及其边界条件如下：

并分析Brooks-Corey模型描述土壤非饱和导水率、土壤基质势和土壤含水量之间关系的表达式分别如下：

式中：h为土壤基质势(cm)。

通过最小作用量原理，通量假定以及积分中值定理，获得了描述土壤水分上渗的剖面方程，其中，采用的通量假定与积分中值定理是用以简化剖面方程的函数形式，其中通量假定的表达式为：

式中：J_w0是进水口处的水流通量(cm min^-1)，J_w是土壤水流通量(cm min^-1)，K_i是对应于土壤初始含水量时的非饱和导水率(cm min^-1)；其他参数意义同上。

再基于质量守恒原理与剖面方程，推导得出完整的描述恒定水头边界条件下均质干土一维水分上渗的解析模型(6)–(8)。

其中，采用的质量守恒原理是用以获得累积入渗量I与时间t之间的函数关系，其中质量守恒原理的具体物理表现形式如下：

I＝Uz_f (8)

其中，

s²＝2K_e(θ_s-θ_i)(h_p-H_f) (10)

基于上述公式(6)–(13)，本发明设计了一种基于环刀土样上渗过程的土壤水力学特性预测方法，在实际应用当中，如图1所示，执行如下步骤A至步骤G，实现对环刀土样的土壤水力学特性预测。

步骤A.针对环刀土样对应恒定下边界水头h_p下均质干土一维水分上渗过程，记录获得环刀土样累计入渗量I随累计时长t的变化关系，土壤中湿润锋前进距离到达环刀土样表面时的累计时长t_e，以及测量获得环刀土样的饱和导水率K_s，然后进入步骤B。

步骤B.根据环刀土样累计入渗量I随累计时长t的变化关系，获得累计时长t_e所对应的累计入渗量I_te、以及土壤中湿润锋前进距离z_fte，并获得z_fte与I_te的比值，构成环刀土样土壤湿润区含水量的平均增长量U，然后进入步骤C。

步骤C.根据环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r，结合环刀土样土壤湿润区含水量的平均增长量U，按如下公式：

获得环刀土样土壤水分剖面的形状系数a，然后进入步骤D。

步骤D.根据环刀土样土壤水分剖面的形状系数a，结合如下公式：

应用不动点迭代算法，获得环刀土样的土壤孔隙尺寸分布指数n，然后进入步骤E。

步骤E.根据环刀土样的土壤孔隙尺寸分布指数n、环刀土样的饱和导水率K_s，结合环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r、以及环刀土样土壤水分剖面的形状系数a，按如下公式：

获得环刀土样的有效饱和导水率K_e，其中，m＝3n+2，然后进入步骤F。

步骤F.根据环刀土样的有效饱和导水率K_e、以及环刀土样累计入渗量I随累计时长t的变化关系，按如下公式：

执行拟合操作，获得环刀土样的吸渗率s。，然后进入步骤G。

步骤G.根据环刀土样的吸渗率s、环刀土样的土壤孔隙尺寸分布指数n、环刀土样的土壤饱和含水量θ_s、环刀土样的土壤初始含水量θ_i、环刀土样的饱和导水率K_s、以及恒定下边界水头h_p，按如下公式：

获得环刀土样的土壤进气吸力值h_d，其中，m＝3n+2。

为了验证上渗解析模型的准确性，如图2a、图2b、图2c分别展示了在三种典型土壤质地中(土壤水力学性质见表1)利用解析解模拟的累计入渗量曲线与对公式(1)直接进行数值求解结果的差异比较。从图中可以看出解析解结果与数值模拟结果几乎是相一致的，二者间的相对误差分别为2.17％、1.12％和1.56％，决定系数分别为0.996、0.999和0.998，因此该上渗解析模型可以准确模拟不同质地土壤的上渗过程，为获取相应的水力学参数提供了理论基础。

将本发明设计应用于实际当中，采用以下三种典型质地土壤(土壤水力学性质见表1)、不同环刀高度的土样试验案例进一步说明本发明公布的一种基于水分上渗过程快速估计土壤水力学参数的方法。

表1

注：表中参数意义同上。

①对于环刀高度为5cm，质地为黏土的土样，在一维水分上渗试验中下边界水头设置为-37.31cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的黏土土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝7800s，以及测定土样的饱和导水率Ks＝0.001cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝1.536cm，湿润锋前进距离为z_fte＝5cm，根据公式(8)，计算得到U＝0.307，再根据公式(14)得到α＝0.139，最后根据公式(15)获得土样的孔隙分布指数为0.120。

步骤E-G：根据公式(16)计算得到Ke＝0.00031cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.1371cm min^-0.5，再根据公式(18)计算得到进气吸力值为36.98cm。

②对于环刀高度为5cm，质地为粉沙土的土样，在一维水分上渗试验中下边界水头设置为-20.75cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的粉沙土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝1913s，以及测定土样的饱和导水率Ks＝0.011cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝1.981cm，湿润锋前进距离为z_fte＝5cm，根据公式(8)，计算得到U＝0.396，再根据公式(14)得到α＝0.193，最后根据公式(15)获得土样的孔隙分布指数为0.186。

步骤E-G：根据公式(16)计算得到Ke＝0.0037cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.3593cm min^-0.5，再根据公式(18)计算得到进气吸力值为19.58cm。

③对于环刀高度为5cm，质地为砂壤土的土样，在一维水分上渗试验中下边界水头设置为-14.66cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的粉沙土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝694s，以及测定土样的饱和导水率Ks＝0.043cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝1.522cm，湿润锋前进距离为z_fte＝5cm，根据公式(8)，计算得到U＝0.304，再根据公式(14)得到α＝0.255，最后根据公式(15)获得土样的孔隙分布指数为0.290。

步骤E-G：根据公式(16)计算得到Ke＝0.0156cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.4714cm min^-0.5，再根据公式(18)计算得到进气吸力值为12.94cm。

④对于环刀高度为10cm，质地为黏土的土样，在一维水分上渗试验中下边界水头设置为-37.31cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的黏土土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝539min，以及测定土样的饱和导水率Ks＝0.001cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝3.059cm，湿润锋前进距离为z_fte＝10cm，根据公式(8)，计算得到U＝0.306，再根据公式(14)得到α＝0.143,最后根据公式(15)获得土样的孔隙分布指数为0.123。

步骤E-G：根据公式(16)计算得到Ke＝0.00031cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.1371cm min^-0.5，再根据公式(18)计算得到进气吸力值为36.36cm。

⑤对于环刀高度为10cm，质地为粉沙土的土样，在一维水分上渗试验中下边界水头设置为-20.75cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的粉沙土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝138min，以及测定土样的饱和导水率Ks＝0.011cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝3.926cm，湿润锋前进距离为z_fte＝10cm，根据公式(8)，计算得到U＝0.393，再根据公式(14)得到α＝0.204最后根据公式(15)获得土样的孔隙分布指数为0.203。

步骤E-G：根据公式(16)计算得到Ke＝0.0038cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.3638cm min^-0.5，再根据公式(18)计算得到进气吸力值为19.77cm。

⑥对于环刀高度为10cm，质地为砂壤土的土样，在一维水分上渗试验中下边界水头设置为-14.66cm。

由上述的计算参数，土壤水力学参数的计算过程如下。

步骤A：在恒定下边界水头下的粉沙土样一维水分上渗试验中，记录累积入渗量与时间的数据，记录湿润锋到达土样上表面时间t_e＝54min，以及测定土样的饱和导水率Ks＝0.043cm min^-1。

步骤B-D：根据湿润锋到达土样上表面的时间，查询累计入渗量与时间的数据，得到此时刻下的累积入渗量I_te＝2.992cm，湿润锋前进距离为z_fte＝10cm，根据公式(8)，计算得到U＝0.299，再根据公式(14)得到α＝0.275，最后根据公式(15)获得土样的孔隙分布指数为0.333。

步骤E-G：根据公式(16)计算得到Ke＝0.0161cm min^-1，通过Ke与累积入渗量的时间序列数据，通过公式(17)拟合得到吸渗率为0.4758cm min^-0.5，再根据公式(18)计算得到进气吸力值为13.51cm。

从以上六个具体计算结果(表2所示基于本发明计算得到的土壤水力学参数与被试土壤参数之间的比较)可以看出，本发明计算得到的土壤水力学参数与被试土壤参数间的差距不大。对于土壤孔隙分布指数而言，二者间最大的误差仅0.032，发生在环刀高度为5cm的砂壤土中；对于吸渗率而言，二者是几乎相一致的；对于进气吸力值而言，二者间最大的误差仅1.72cm，发生在环刀高度为5cm的砂壤土中。因此，本发明可以基于恒定水头边界条件下均质干土一维水分上渗过程，准确的预测出环刀土样的水力学特性，尤其是对于质地偏黏的土样。

表2

注：表中参数意义同上，其中n′、s′和h_d′表示对于环刀土样高度为5cm的估计参数，n″、s″和h_d″表示对于环刀土样高度为10cm的估计参数。

上述技术方案所设计基于环刀土样上渗过程的土壤水力学特性预测方法，采用全新设计方法，针对恒定水头边界条件下均质干土一维水分上渗解析解的逆过程，通过记录环刀土样上渗过程中累积入渗量与时间的数据、湿润锋到达土样上表面的时间、以及测定土样饱和导水率，快速获取土壤水力学特性，克服了现有技术测定土壤水力学特性耗时长、效率低，且不适于田间大面积快速获取等问题，避免了参数多解、不收敛等问题，提高了土壤水力学特性的获取效率，为田间大面积的土壤水力特性调查提供了有效手段。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (7)

1.一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于，执行如下步骤，实现对环刀土样的土壤水力学特性预测；

步骤A.针对环刀土样对应恒定下边界水头h_p下均质干土一维水分上渗过程，记录获得环刀土样累计入渗量I随累计时长t的变化关系，土壤中湿润锋前进距离到达环刀土样表面时的累计时长t_e，以及测量获得环刀土样的饱和导水率K_s，然后进入步骤B；

步骤B.根据环刀土样累计入渗量I随累计时长t的变化关系，获得累计时长t_e所对应的累计入渗量I_te、以及土壤中湿润锋前进距离z_fte，并获得z_fte与I_te的比值，构成环刀土样土壤湿润区含水量的平均增长量U，然后进入步骤C；

步骤C.根据环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r，结合环刀土样土壤湿润区含水量的平均增长量U，获得环刀土样土壤水分剖面的形状系数a，然后进入步骤D；

步骤D.根据环刀土样土壤水分剖面的形状系数a，获得环刀土样的土壤孔隙尺寸分布指数n。

2.根据权利要求1所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：还包括步骤E如下，执行完步骤D之后，进入步骤E；

步骤E.根据环刀土样的土壤孔隙尺寸分布指数n、环刀土样的饱和导水率K_s，结合环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r、以及环刀土样土壤水分剖面的形状系数a，获得环刀土样的有效饱和导水率K_e，然后进入步骤F；

步骤F.根据环刀土样的有效饱和导水率K_e、以及环刀土样累计入渗量I随累计时长t的变化关系，获得环刀土样的吸渗率s。

3.根据权利要求2所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：所述步骤E中，根据环刀土样的土壤孔隙尺寸分布指数n、环刀土样的饱和导水率K_s，结合环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r、以及环刀土样土壤水分剖面的形状系数a，按如下公式：

获得环刀土样的有效饱和导水率K_e，其中，m＝3n+2。

4.根据权利要求2所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：所述步骤F中，根据环刀土样的有效饱和导水率K_e、以及环刀土样累计入渗量I随累计时长t的变化关系，按如下公式：

执行拟合操作，获得环刀土样的吸渗率s。

5.根据权利要求2所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：还包括步骤G如下，执行完步骤F之后，进入步骤G；

步骤G.根据环刀土样的吸渗率s、环刀土样的土壤孔隙尺寸分布指数n、环刀土样的土壤饱和含水量θ_s、环刀土样的土壤初始含水量θ_i、环刀土样的饱和导水率K_s、以及恒定下边界水头h_p，按如下公式：

获得环刀土样的土壤进气吸力值h_d，其中，m＝3n+2。

6.根据权利要求1所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：所述步骤C中，根据环刀土样的土壤饱和含水量θ_s、土壤初始含水量θ_i、土壤残余含水量θ_r，结合环刀土样土壤湿润区含水量的平均增长量U，按如下公式：

获得环刀土样土壤水分剖面的形状系数a。

7.根据权利要求1所述一种基于环刀土样上渗过程的土壤水力学特性预测方法，其特征在于：所述步骤D中，根据环刀土样土壤水分剖面的形状系数a，结合如下公式：

应用不动点迭代算法，获得环刀土样的土壤孔隙尺寸分布指数n。

Images