基于完全非线性流动过程响应的非饱和土壤水分扩散度的测定方法

技术领域

本发明涉及土壤水分运动研究技术领域，特别是涉及一种非饱和土壤水分扩散度的测定方法。

背景技术

随着社会的进步与经济的快速发展，有关水资源的矛盾和问题日益尖锐，进一步加大了经济社会可持续发展面临水资源短缺的压力，如何保持农业的可持续发展和水资源的有效利用率成为亟待解决的问题。研究土壤水分运动与综合利用成为社会发展中的重要课题。土壤水分运动是一个复杂的过程，与灌溉方式、土壤条件和外界环境等众多因素密切相关。土壤水分扩散度是土壤水分运动方程的重要参数，现有的土壤扩散度的测定方法可以归纳为2类。

第一类方法：

以土壤含水率θ为变量的土壤水分运动方程为：

其中，θ为土壤含水率，x和t分别为空间和时间坐标，D(θ)为土壤扩散度，K为水力传导度。

现有的方法将方程(1)由二次偏微分方程变化为二次常微分方程后，获得解析解。然后在解析解的对应的定解条件(边界条件和初始条件)下进行测试，根据测试过程中各监测点水分变化，基于解析解反演D(θ)。

第二类方法：

土壤扩散度D(θ)与土壤非饱和水力传导度和单位含水率变化下的水势变化之间的关系可以表示为：

其中，h为土壤基质势，θ为土壤含水率，K(θ)为土壤非饱和水力传导度，为单位含水率变化下的水势变化。

基于这一性质，通过测定土壤水分特征曲线确定和非饱和水力传导度K(θ)，根据(2)式确定D(θ)。

需要指出，对于第一类方法，由于解析解得求解需要，对于偏微分方程转化为常微分过程中需要对(1)式的右侧两项分别进行线性化处理，也就是用两个含水率之间的D(θ)的割线代替测点位置的切线，也就是说基于解析解测定的D(θ)实际上是一个含水率条件下的平均值，而不是与θ的直接对应值，尽管在D(θ)～θ关系接近线性关系时，并不会产生严重的误差，然而，对于大部分土壤，特别是D(θ)～θ表现出显著的非线性关系的土壤，误差就非常显著，而且，就测试理论而言，采用分段线性的方法近似D(θ)～θ这一非线性关系，是没有办法根本解决这一测试问题的。

对于第二类方法，由于无法实现同一含水率(或基质势下)的土壤特征曲线和非饱和水力传导度的同步测定。因此更多的情况下是根据拟合确定的土壤水分曲线和非饱和水力传导度关系，在某一含水率情况下直接计算D(θ)，这样D(θ)就受到土壤水分特征曲线测量误差、非饱和水力传导度测量误差、和土壤水分特征曲线拟合误差的影响，一些情况下(特比是土壤中粘粒含量超过30％)，误差与D(θ)在同一个数量级。

此外，这两种方法还存在一个共性问题，均不能实现非饱和含水率区间的完全覆盖，其中的原因在于：在土壤含水率较小的情况下，土壤中水流运动主要在土壤固体颗粒的分子力的作用下发生移动，而随着土壤含水率的增加，土壤毛管力的作用增强。由于作用力不同，土壤中非线性流动关系随着土壤质地的变化差异性很大。现有的方法的含水率的测定范围很难覆盖整个非饱和含水率区间。

如上所述，一套完整的测定土壤水分扩散度的系统还没有建立起来，无法为解决农业水资源高效利用等问题提供科学支持。

发明内容

基于现有的方法存在的问题，本发明提出一种基于完全非线性流动过程响应的非饱和土壤水分扩散度的测定方法，在测试方法和解析方法上实现D(θ)～θ非线性关系的直接测定。相比现有的方法，本发明的先进性在于所提出的测试机理能够充分的反应D(θ)的物理性质，解决传统的方法对土壤非饱和流动这种非线性过程进行线性化处理这一重大测试机理问题。

本发明方法是一种基于水平非饱和渗吸过程和水平穿透过程连续测定的土壤非饱和水分扩散度的新技术，为更加精准地研究非饱和土壤水分运动、实现农业水资源的高效利用提供支持。

本发明所提供的技术方案具体如下：

一种基于完全非线性流动过程响应的非饱和土壤水分扩散度的测定方法，通过水平渗吸过程测定和水平穿透过程测定，实现覆盖整个非饱和含水率范围的非饱和土壤扩散度测定，包括以下步骤：

步骤1，准备水平土柱试验装置，该水平土柱试验装置包括依次连通的供水平水装置、水室、水平土柱、渗出通量测定装置，供水平水装置用于控制土柱的作用水头，水平土柱由左至均匀分隔成若干节子土柱，每节子土柱内布置土壤含水率及基质势测定传感器，各传感器与数据采集器连接；所述供水平水装置用于在控制水头高度不变的情况下，对实验土样进行供水；

步骤2，对非饱和土壤进行水平渗吸过程测定，具体实现如下：

①制备土样及装填土柱，准备好足够的试样，要求试样风干含水量；

②准备土柱，将试样按一定容重装填土柱，然后将土柱水平放置，并安装于水平土柱试验装置中；

③根据测试水头高度要求对供水平水装置进行水位设定；计时并记下供水平水装置初始水位读数；

④经过一定时间，且湿润锋面未达土柱末端之前，结束试验，此时停止供平水装置对土柱的供水，记录整个试验的历时，总入渗水量和各点含水率；

形成边界出流后，由水平渗吸过程转化为水平穿透过程；

步骤3，对非饱和土壤进行水平穿透过程测定，具体实现如下：

在土柱进水边界形成出流后，在试验过程中变换土柱出水边界水头，测定不同时刻t时刻的排水量；

步骤4，根据步骤2中记录的总入渗水量及各点含水率绘制θ～η关系曲线，从该曲线可以算得总入渗水量，其值应与实测入渗水量相当，如果相差过大，则应该检查是否有计算错误，若差别不大，则可对θ～η曲线作局部的修正，使其尽量接近实测值；

步骤5，根据步骤4绘制的θ～η关系曲线，按一定的θ值分割θ～η曲线，用图解分析法求D(θ)值，并绘制D(θ)～θ曲线；

步骤6，根据步骤3中的水平穿透过程对应的试验数据，绘制关系曲线，根据其斜率，确定相应非饱和含水率区间的扩散系数。

上述测定方法中，对非饱和土壤进行水平渗吸过程测定，具体数据解析方法如下：

非饱和水平渗吸过程以土壤水分在水平半无限边界均质土柱中发生水平运动的情况下进行分析计算的；已知一维土壤水平运动微分方程为：

式中：θ是入渗时刻为t时与输入端即进水边界的水平距离为x处的土壤含水率，D(θ)为土壤扩散度；

对方程式(1)，采用Boltmann变换法，得：

式(2)反映D(θ_x)与θ_x的关系，θ_x为x处的含水率；通过水平土柱试验，在某一时刻t，测定不同距离x处的含水率θ_x，即可绘制θ～η曲线，基于θ～η关系曲线，采用图解法确定值后，采用公式(1)确定D(θ_x)～θ_x关系；其中，θ是与输入端即进水边界的水平距离为x处的土壤含水率；为了求解方程式(1)，采用Boltmann变换法，以一个新的变量即含x和t两个分离变量代入方程使之变为常微分方程；设θ(x，t)＝f[η(x，t)]，而令η(x，t)＝xt^-1/2。

上述测定方法中，对非饱和土壤进行水平穿透过程测定，具体数据解析方法如下：

维持某一负压水头H₀达到稳定后，瞬时加大负压值即增加ΔH的负压，然后测定试样的出流过程，利用多步法中的解和所介绍的资料分析方法确定扩散系数，

h＝－H₀，t＝0，0≤x≤l(4)

其中，K为水力传导度，D为扩散度，K_p为多孔板的导水率；△H为在开始时刻的瞬时变压值；l为试样的厚度；l_p为多孔板的厚度；h为土壤基质势；t为时间，x为水平方向坐标；

式(3)为不考虑试样的重力的土壤水运动方程，式(4)(5)(6)为定解条件；在这些定解条件下，得到：

式中，W_t和W_∞分别为t时刻的测定排水量和累积排水量，l为试样的厚度，其中，K_p，l_p分别为多孔板的导水率和厚度，Rl＝a₁tga₁；a和a₁分别为a_n的零阶和一阶级数项，当时可取级数的第一项，其误差仅1％；其中i为级数，P_n为土壤水分运动方程进行拉普拉斯变化后的常数项，D为平均扩散度；

当Rl为定值时，则在半对数纸上呈直线关系，若直线与水平线夹角为β，其直线斜率为：

tanβ＝-0.4343a₁²(8)

因此，根据实测右侧边界的排水过程，据绘制关系曲线，根据其斜率，即可求得扩散度D值。

本发明非饱和土壤水平渗析过程和水平穿透过程测定机理具体如下：

1.非饱和土壤水平渗析过程测定机理

土壤水分在较长的(理论上说是水平半无限边界)均质土柱中、在以土壤分子力作用下发生水平运动的情况下进行分析计算的。已知一维土壤水水平运动(忽略重力作用)微分方程为：

式中：θ是与输入端(即进水边界)的水平距离为x处的土壤含水率。为了求解方程式(1)，采用Boltmann变换法。以一个新的变量(含x和t两个分离变量)代入方程使之变为常微分方程。

设

θ(x，t)＝f[η(x，t)]

而令η(x，t)＝xt^-1/2为新变量，这样初始条件及边界条件将变成：

θ＝θ_i，η＝0(η→∞)

θ＝θ_s，η＝0

对方程式(1)进行代换，由于

代入式(1)，得：

式(2)表明，当进行Boltmann变换后，式(1)变为常微分方程。式(2)两边乘以dη，并从θ_i积分到θ_x(x处的含水率)，得

最后得：

式(4)反映D(θ_x)与θ_x的关系，通过水平土柱试验，在某一时刻t，测定不同距离x处的含水率θ_x，即可绘制θ～η曲线，有了θ～η关系曲线，则值和均可用图解法求得，按式(4)进行计算，确定D(θ_x)～θ_x关系。

2.非饱和土壤水平穿透过程测定机理

在实验装置右侧出流后，由水平渗吸过程测定转换为水平穿透过程。

维持某一负压水头H₀，并已达到稳定后(即土壤含水率不变)，瞬时加大负压值(即增加ΔH的负压)，然后测定试样的出流过程，利用多步法中的解和所介绍的资料分析方法确定扩散系数。

h＝－H₀，t＝0，0≤x≤l(6)

其中，K为水力传导度，D为扩散度，K_p为多孔板的导水率；△H为在开始时刻的瞬时变压值；l为试样的厚度；l_p为多孔板的厚度。

式(5)为不考虑试样的重力的土壤水运动方程，式(6)(7)(8)为定解条件。

根据容水度的定义：

则式(5)可表示为：

即：

在右侧边界：从土样中排出的流量等于通过多孔处板的流量，即:

式中：Kp，lp分别为多孔板渗透系数和厚度；h为通过多孔板的压力差；l为土样长度。

对式(9)进行Laplace变换：

式(9)为：

上式通解为：

由式(12)得：

边界条件(式(7)和(8))，经laplace变换为：

由式(14)和(15)可知：

∴c₁＝c₂＝c(16)

将式(16)代入式(15),得：

令代入(17)式，经整理后得：

将式(18)代入式(16)，得：

式中：ψ(P)均为多项式。

将式(19)进行逆变换，根据留数定理得：

令：则式(20)为：

即rl＝atga(21)

由式(17)得：

将式(20)代入式(22)，且将P改写为P_n，得

由式(19)可知：

由式(22)、(24)得：

ψ′(0)＝1；

根据分解定理：

代入(25)式：

实验体的水断面面积为F，则水量Q_t为：

式中：V为土样体积。

由式(28)：

当t＝0时，Q_i＝AV+B·V·△h；

t→∞时，Q_f＝AV；

因此在负压增加某一△h值时，排水总量为：

W_∞＝Q_i-Q_f＝B·V·△h

即：

由式(29)可知

∴B实质上为平均的容水度(在负压为某一值时)。

在负压增加某一△h值时，t时刻的累计排水量为

当时可取级数的第一项，其误差仅1％。式(29)可写为

等式两边取对数

由式(31)可知，右端第一项为截距，即当时，仅与Rl有关。当Rl为定值时，则在半对数纸上呈直线关系。若直线与水平线夹角为β，其直线斜率为：

tanβ＝-0.4343a₁²(34)

因此，根据实测右侧边界的排水过程，据绘制关系曲线，根据其斜率，即可求得扩散度D值。

本发明中，供水平水装置的玻璃管管口深入到供水槽体的底部，这样供水槽体的静水压便等于此玻璃管口的液面高度，供水槽体内水体流出后，槽体内形成真空，空气只从玻璃管中进入，玻璃管口即为接触空气点，供水槽体内水体水位不低于玻璃管出口，玻璃管口以上供水槽体内水体的增减将不影响静水压，从而自动保持了流速的恒定。

与现有技术相比，本发明具有以下有益效果：

1.本发明通过水平渗吸过程测定和水平穿透过程测定，实现了非饱和含水率区间的完全覆盖的非饱和土壤水分扩散度D(θ)的测定；其中水平渗吸过程测定土壤含水率在田间持水率以下区段的的水分扩散度，而水平穿透过程测定土壤含水率在田间持水率～饱和含水率之间区段的土壤水分扩散度。

2.本发明测试原理能够有效的反映非饱和水流运动非线性流动这一物理特性，而不是现有方法中的线性近似，测定原理具有更为准确的物理机制。

3.本发明提出了一种具有完备土壤非饱和流动物理机制的非饱和扩散度D(θ)的测定机理，以及机理的具体实现。实际上，这种测试方法本身的测定原理的物理基础非常复杂，而我们为了使得本发明便于应用，对于两个过程都在原理的基础上，提出了较为容易实现的测定方法：采用水平土柱试验装置，在测试开始时，通过测定水平方向各个位置土壤含水率的变化，对水平渗吸过程进行测定；当渗透锋面达到右侧出水边界，形成出流后，改变左侧进水控制边界条件，通过测定水平方向各个位置的土壤基质势、以及边界渗出通量，对水平穿透过程进行测定，当渗出通量稳定时，土壤已经完全饱和，测定结束。这样就能够实现土壤从干土(含水率为0)到饱和含水率之间的D(θ)的测定，实现非饱和含水率区间的完全覆盖。

4.本发明实验方法更加简单，准确；采用Boltmann变换法变换方程，更加高效；克服了现有测定技术复杂、耗时、费用庞大的不足。

5.本发明的实验方法，提高了土壤水分扩散度D(θ)的计算精度，解决一般测定方法费时、费力、费资金的缺陷，信息明确直观，且在该领域具有独创性。

附图说明

图1为本发明中水平土柱试验装置结构示意图；

图2为本发明中供水平水装置结构示意图；

图3为本发明实施例的θ～η关系曲线；

图4为本发明实施例在半对数纸上绘制的右侧边界相对出流量关系图。

图中：1-供水平水装置；1a-控制水位的玻璃管；1b-供水槽体；1c-出水口；2-水室；3-排气孔；4-前过滤板；5-土柱；6-土壤含水率及基质势测定传感器；7-数据采集器；8-后过滤板；9-渗出通量测定装置。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

＜实施例＞

本实验通过两种实现方法来详细阐述本发明的方案，具体实施步骤如下：

①制备土样及装填土柱，准备好足够的试样(风干含水量)。按一定容重装填土柱，将螺杆旋紧然后水平放置。

②瞬时给水室充水并使供水平水装置定位。

③计时并记下供水箱初始水位读数。

④经过相当时间(湿润锋面未达土柱末端之前)后，即可结束试验，此时停止供水，松开坚固螺杆，按节取出土壤测定含水率，记下整个试验的历时及总水量，清洗器皿，试验结束。

绘制θ～η关系曲线，如图3所示，从该曲线可以算得总入渗水量，其值应与实测值W相当，如果相差过大，则应该检查是否有计算错误，若差别不大，则可对θ～η曲线作局部的修正，使其尽量接近实测值。

对于水平渗析过程：

有了θ～η关系曲线后，按一定的θ值分割θ～η曲线，见图3。设土柱初始含水率θ_i为3％，若求θ_x＝8％时的D(θ)，只要求得及两项(η_i和η_x为θ_i为3％，和θ_x＝8％对应的η值)，则D(θ)即可求得。在θ_x＝8％处求取θ～η曲线的斜率再用体积法量取abdc的面积(图3)，即可求得值(θ_i＝3％，θ_x＝8％)，将值和便算得θ＝8％时的D(θ)值，由于值为负数，故算出的D(θ)值为正数。按同样方法，可求得其他含水率值的D(θ)。

对于水平穿透过程：

右端第一项为截距，即当时，仅与Rl有关。当Rl为定值时，则在半对数纸上呈直线关系，见图4。若直线与水平线夹角为β，其直线斜率为：

tanβ＝-0.4343a₁²

根据根据实测右侧边界的排水过程，据绘制关系曲线，根据其斜率，即可求得扩散度D值。

本发明公开了一种土壤水分扩散度的测定方法，包括水平渗吸过程测定及水平穿透过程测定扩散度试验。首先按照水平土柱试验装置图准备试验，记录试验历时，总入渗水量及各点含水率等；然后绘制θ～η关系曲线并对θ～η曲线作局部的修正，使其尽量接近实测值；按一定的θ值分割θ～η曲线，用图解分析法求D(θ)值，并绘制D(θ)～θ曲线；最后利用砂性漏斗测定水分特征曲线的装置，在试验过程中观测出流过程，确定扩散系数。

本发明在特定水平土柱试验装置中进行，供水平水装置用以控制水平土柱的作用水头，

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。