CN101788449B - 一种用于测量土壤垂向饱和渗透系数的方法 - Google Patents

Abstract

本发明是一种用于测量土壤垂向饱和渗透系数的方法，属于水文测验领域，该方法使用单环、自记水位计、土壤含水量测定仪的设备测量，然后根据观测的单环水位变化过程，利用在Philip与Nestingen土壤水运动计算方法基础上，推导出的描述单环变水头入渗试验中水深变化过程的计算公式，采用最小二乘原理计算出土壤垂向饱和渗透系数。该方法可用于野外迅速准确的测定土壤垂向饱和渗透系数。该仪器需水量少，设备简易，操作简单，野外携带方便，并且测量时入渗单环的半径不受限制，可用于较大尺度测量土壤渗透系数。适用于水利行业、水文地质行业及农田水利行业用于测量土壤饱和渗透系数，易于推广使用。

Description

一种用于测量土壤垂向饱和渗透系数的方法

技术领域

本发明一种用于测量土壤垂向饱和渗透系数的方法，属于水文测验领域，尤其是一种野外便携式设备测量土壤垂向饱和渗透系数的方法。

技术背景

野外试验测定土壤渗透系数等流域水文地质参数对于研究降水入渗过程、流域产汇流规律具有重要意义。目前野外常用测量渗透系数的方法有：单/双环定水头法，井流法，吸力盘法等。但是上述方法存在如下缺陷：(1)需要大量注水以维持定水头；(2)测量的尺度都比较小。

在单环变水头入渗试验推求土壤渗透系数等水力特征参数的计算方法方面，Philip于1992年基于Green-Ampt入渗模型(简称GA模型)导出了一维变水头GA入渗公式，1993年基于土壤水球面扩散假说及GA模型导出了三维井流变水头入渗公式，并提出了两点法推求土壤吸力(C)与垂向渗透系数(K)的方法。Regalado运用Philip三维井流变水头入渗方法的室内外试验成果，得出了推求C、K的经验公式。Nestingen在Philip研究基础上进一步考虑单环对水流影响，提出了球缺扩散假说，导出了单环变水头入渗试验推求C、K的方法。然而Nestingen方法存在如下概念上的偏差及解法上的缺陷：(1)仅考虑土壤水三维入渗，实际上由于单环插入土壤一定深度，因此土壤水进入三维入渗前需先经历一维入渗；(2)球缺扩散假说只有当作用在湿润峰上的力大小相等且都沿球面法线方向的情况下才成立，而三维土壤入渗水在受重力作用下，湿润峰“受力”不均，因此使用球缺扩散假说需从土壤入渗水中除去重力入渗水；(3)采用差分离散法求解微分方程，其数值计算结果受限于实测资料的精度，容易造成计算结果失真。为此，本发明借鉴Philip及Nestingen有关土壤水运动的研究方法，将非饱和一维入渗与三维入渗水运动相结合，推导出在单环变水头入渗试验中描述单环水深变化过程的公式，提出了推求土壤垂向渗透系数的改进Nestingen方法(Modified Nestingen方法，简称MN方法)。

发明内容

本发明一种用于野外测定土壤垂向饱和渗透系数的方法，该方法具有需水量少，设备简易，操作简单，野外携带方便等优点，并且其半径可变，适用于较大尺度测量土壤渗透系数。

一种用于测量土壤垂向饱和渗透系数的方法，其特征在于该方法包括以下步骤：

(1)将单环缓慢均匀插入土壤，记录单环插入土壤的深度；

(2)将土壤含水量测量仪插入单环内的土壤中，测量记录土壤初始含水量；

(3)将自记水位计放入单环内，然后往单环内注水，记录单环水位变化，即时间序列{t_i}及单环内水深变化序列{H_i}；

(4)在单环中水入渗结束时，再用土壤含水量测量仪测定单环内土壤含水量，此时所测得单环内土壤含水量作为土壤饱和含水量；

(5)根据单环变水头入渗试验中水深变化过程的计算原理，将单环内水流入渗分为两个阶段，单环中的一维入渗及出单环后的三维入渗。单环中一维入渗采用一维变水头Green-Ampt入渗模型计算：

$t=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(\frac{H_0-j}{\mathrm{\Δ\θ}}\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+(H_0-h)\×(1-\mathrm{\Δ\θ})/\mathrm{\Δ\θ}}\right))$ H₀-h≤LΔθ   (1)

式中：Δθ为土壤缺水量；H₀为单环初始水深；C为土壤吸力；h为单环内水深；t为入渗时间；K为土壤垂向渗透系数；L为单环插入土壤深度；

当湿润峰到达单环出口(H₀-h＞LΔθ)后，土壤水进入三维入渗阶段。出单环后的三维入渗包含两部分：重力引起的下渗土壤水以及单环水压力与土壤吸力引起的土壤水扩散，重力入渗水量等于渗透系数乘入渗时间，土壤水扩散采用Nestingen提出的球缺扩散模型计算。对于土壤水扩散过程，假设土壤均质且各向同性，由球缺面积与单环面积相等，可得初始球缺半径r₀：

$r_0=\sqrt{\frac{L^2}{4}+\frac{r_1^2}{2}}-\frac{L}{2}---\left(2\right)$

式中r₁为单环半径。

在球缺扩散中，由水量平衡原理，当土壤水扩散半径为R时，在球缺初始半径r₀处土壤水扩散速度v为：

$v=\frac{\mathrm{\Δ\θ}}{\mathrm{\π}{r}_1^2}\frac{\mathrm{dV}}{\mathrm{dt}}=\frac{2\mathrm{\Δ\θ}}{r_1^2}R(R+L)\frac{\mathrm{dR}}{\mathrm{dt}}---\left(3\right)$

式中V为球缺扩散中土壤水扩散体积；π为圆周率。

根据Green-Ampt模型原理，每个球缺面上水通量相等，则当球缺扩散半径为R时，在半径r(r₀≤r≤R)处，土壤水扩散速度v_r为：

$v_r=\frac{r{S}_0}{S_r}=\frac{2\mathrm{\Δ\θ}}{r_1^2}R(R+L)\frac{\mathrm{dR}}{\mathrm{dt}}\frac{r_0(L+r_0)}{r(L+r)}---\left(4\right)$

式中S₀为半径为r₀时的球缺表面积；S_r为半径为r时的球缺表面积。

根据达西定理，从球缺半径r₀到R处的土壤水扩散的势能P为：

$P={\∫}_{r_0}^R\frac{v_r}{K}\mathrm{dr}=\left(\frac{2\mathrm{\Δ\θ}}{K{r}_1^2}R(R+L)\frac{\mathrm{dR}}{\mathrm{dt}}\right)r_0(L+r_0){\∫}_{r_0}^R\frac{1}{r(L+r)}\mathrm{dr}---\left(5\right)$

利用球缺扩散与实际单环扩散之间的几何转换系数，可推导出从单环三维入渗出口到土壤水实际扩散锋面的土壤水扩散的势能P为：

$P=\frac{{\mathrm{\π}}^2}{8}\frac{2R(R+L)\mathrm{\Δ\θ}}{K{r}_1^2}\frac{\mathrm{dR}}{\mathrm{dr}}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)---\left(6\right)$

势能P又等于单环三维入渗出口处水压力与土壤吸力之和：

$P=h+L+\frac{L}{K}\frac{\mathrm{dh}}{\mathrm{dt}}+C---\left(7\right)$

单环入渗水量应等于重力入渗水量与单环水压力及土壤吸力引起的扩散水量之和，由重力入渗水速度K，可得h与R的关系式：

$h=H_0-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{2R^3+3{\mathrm{LR}}^2-2r_0^3-3{\mathrm{Lr}}_0^2}{3r_1^2}\mathrm{\Δ\θ}---\left(8\right)$

式中t₀为一维土壤水入渗结束时间，由式(1)得：

$t_0=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(L\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+L\×(1-\mathrm{\Δ\θ})}\right))---\left(9\right)$

式(8)求导得： $-\frac{\mathrm{dh}}{\mathrm{dt}}=\frac{2R(R+L)\mathrm{\Δ\θ}}{r_1^2}\frac{\mathrm{dR}}{\mathrm{dt}}+K---\left(10\right)$

将式(6)、(8)和(10)代入式(7)可得描述单环三维扩散的微分方程：

$\left\{\begin{array}{c}\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{H_0+C-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{2R^3+3{\mathrm{LR}}^2-2r_0^3-3L{r}_0^2}{3r_1^2}\mathrm{\Δ\θ}}{(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)+L)\frac{2R(R+L)\mathrm{\Δ\θ}}{K{r}_1^2}}\\ {R|}_{t=t_0}=r_0\end{array}\right.---\left(11\right)$

将式(10)代入式(11)，得C的表达式：

$C=(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)+L)\left(\right|\frac{\mathrm{dh}}{\mathrm{Kdt}}|-1)-h---\left(12\right)$

(6)整理试验数据，根据实测时间序列{t_i}及水头序列{H_i}计算出入渗末期单环水位入渗率dH/dt|_end；然后假设土壤饱和渗透系数K，再联立求解下列方程组(13)、(14)和(15)，得出土壤吸力C：

$C=(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R_{\mathrm{end}}}{R_{\mathrm{end}}+L}\frac{r_0+L}{r_0}\right)+L)({\left|\frac{\mathrm{dH}}{\mathrm{Kdt}}\right|}_{\mathrm{end}}-1)-H_{\mathrm{end}}---\left(13\right)$

$H_{\mathrm{end}}=H_0-K(t_{\mathrm{end}}-t_0)-\mathrm{L\Δ\θ}-\frac{2R_{\mathrm{end}}^3+3{\mathrm{LR}}_{\mathrm{end}}^2-2r_0^3-3{\mathrm{Lr}}_0^2}{{3r}_1^2}\mathrm{\Δ\θ}---\left(14\right)$

$t_0=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(L\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+L\×(1-\mathrm{\Δ\θ})}\right))---\left(15\right)$

式中：H_end为入渗末期单环水位；t_end为入渗时间；R_end为入渗末期球缺扩散半径。

然后再将渗透系数K、土壤吸力C及实测时间序列{t_i}代入式(1)求出一维入渗时单环内水深；代入式(8)采用Adams-Bashforth-Moulton变阶法求解该微分方程，获得球缺扩散半径R，然后再将R代入(11)求出三维入渗时单环内水深。此时，求得的单环一维与三维入渗状态下单环内水深系列即为单环变水头入渗模型中实测时间序列{t_i}对应的单环内水深变化序列{h_i}，比较其与实测水头序列{H_i}的差异调整K重新计算{h_i}，直到{h_i}与{H_i}的误差最小，此时的K即为土壤的饱和垂向渗透系数。

所述单环插入土壤的深度2.0～5.0cm。

所述往单环内注水，注到5.0～20.0cm深。

所述自记水位计为HOBO公司生产的U20自记水位计，用于记录单环内水位变化。

所述土壤含水量测定仪为Trase System公司生产的时域反射仪(TDR)，用于测定土壤初始及饱和含水量。

本发明具有如下优点：

由于一种用于测量土壤垂向饱和渗透系数的方法不需要维持定水头，因此需水量较少；并且该发明中使用的设备简易，仅包括平衡杆、锤、单环等器材，野外携带方便，操作也非常简单，如果选用不同半径的单环，还可用于测定不同尺度的土壤渗透系数，非常便于野外土壤垂向饱和渗透系数的测定。

附图说明

图1、本发明的数据处理的流程图；

图2、本发明中将单环插入土壤的示意图；

图3、砂壤土实测与模拟单环水深结果对比图；

图4、粉壤土实测与模拟单环水深结果对比图。

具体实施方式

本发明上述的实例是对本发明的说明而不能限制本发明，在与本发明相当的含义和范围内的任何改变和调整，都应认为是在本发明的范围内。

实施例1

单环1垂直置于土壤上，将平衡杵置于单环1上端，平衡杆采用强度较高的金属材料制作，呈十字架型，用锤敲击平衡杵使单环1缓慢均匀的插入土壤一定深度。记录单环插入土壤的深度2.0～5.0cm；将土壤含水量测定仪3(TDR)插入单环内，测量并记录土壤初始含水量；将自记水位仪2(HOBO U20)放入单环1内，用于记录单环1内水深变化，即记录时间序列{t_i}及对应的单环内水深变化序列{H_i}；然后往单环1内注水，瞬时注到5.0～20.0cm的高度，然后让单环1内水流自由入渗；等单环1水流入渗结束时再用土壤含水量测定仪3(TDR)测定土壤含水量(近似为土壤饱和含水量)。整理试验数据，先计算初始球缺半径 $r_0=\sqrt{\frac{L^2}{4}+\frac{r_1^2}{2}}-\frac{L}{2}$ (其中，r₁为单环半径，L为单环插入土壤的深度)；再根据实测时间序列{t_i}及水头序列{H_i}计算出入渗末期单环水位入渗率dH/dt|_end；然后假设土壤饱和渗透系数K，再联立求解下列方程组①②③求出土壤吸力C：

$C=(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R_{\mathrm{end}}}{R_{\mathrm{end}}+L}\frac{r_0+L}{r_0}\right)+L)({\left|\frac{\mathrm{dH}}{\mathrm{Kdt}}\right|}_{\mathrm{end}}-1)-H_{\mathrm{end}}$ ①

$H_{\mathrm{end}}=H_0-K(t_{\mathrm{end}}-t_0)-\mathrm{L\Δ\θ}-\frac{2R_{\mathrm{end}}^3+3L{R}_{\mathrm{end}}^2-2r_0^3-3{\mathrm{Lr}}_0^2}{3r_1^2}\mathrm{\Δ\θ}$ ②

$t_0=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(L\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+L\×(1-\mathrm{\Δ\θ})}\right))$ ③

式中：H_end为入渗末期单环水位；t_end为入渗时间；R_end为入渗末期球缺扩散半径；Δθ为土壤缺水量等于饱和含水量减初始含水量；H₀为单环初始水深。然后将渗透系数K、土壤吸力C及{t_i}代入式④计算出一维入渗阶段单环内水深h：

$t=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(\frac{H_0-h}{\mathrm{\Δ\θ}}\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+(H_0-h)\×(1-\mathrm{\Δ\θ})/\mathrm{\Δ\θ}}\right))H_0-h\≤\mathrm{L\Δ\θ}$ ④

再采用Adams-Bashforth-Moulton变阶法求解微分方程⑤，求得土壤水扩散半径R：

$\left\{\begin{array}{c}\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{H_0+C-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{{2R}^3+3{\mathrm{LR}}^2-{2r}_0^3-3{\mathrm{Lr}}_0^2}{3r_1^2}\mathrm{\Δ\θ}}{(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)+L)\frac{2R(R+L)\mathrm{\Δ\θ}}{K{r}_1^2}}\\ {R|}_{t=t_0}=r_0\end{array}\right.$ ⑤

然后再将渗透系数K、R及{t_i}代入⑥求出三维入渗阶段单环内水深h：

$h=H_0-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{2R^3+3L{R}^2-2r_0^3-3{\mathrm{Lr}}_0^2}{{3r}_1^2}\mathrm{\Δ\θ},H_0-h>\mathrm{L\Δ\θ}$ ⑥

上述步骤计算的单环内水深变化序列{h_i}即为单环变水头入渗模型中实测时间序列{t_i}对应的单环内水深序列，比较其与单环内水深序列{H_i}的差异调整K重新计算{h_i}，直到{h_i}与{H_i}的误差最小，此时的K即为土壤的饱和垂向渗透系数。

基于数值模拟结果的MN方法验证

Carsel根据大量实测资料总结了12种典型土壤的水力参数。从中选取了三种典型土壤水力参数(表1)，假设土壤初始饱和度为30％。利用美国地质调查局(USGS)建立的三维非饱和带水流计算程序VSF-MODFLOW进行非饱和水分运动数值模拟。该软件借助VGM模型描述土壤水力特性，土壤吸力采用如下公式计算：

式中

为土壤吸力；

为相对渗透性系数；α，n，m为VGM模型参数，其中m＝1-1/n。

VSF-MODFLOW采用有限差分法求解三维理查德(Richards)方程，广泛应用于饱和-非饱和地下水数值模拟。根据土壤经验水力参数(表1)，利用VSF软件模拟单环变水头入渗时单环水深变化过程，以此作为MN方法的拟合目标，推求土壤水力参数，并与一维变水头GA方法(采用最小二乘法求解C，K)、Philip方法和Nestingen方法推求的参数进行比较。

表1土壤水力特征参数

由于单环变水头入渗具有对称性，三维VSF-MODFLOW模型可简化为水平垂向二维数值模型。根据入渗水量在土壤中影响范围大小，模拟范围选为70cm×70cm，剖分为176层×175列。假设单环半径为15cm，插入土壤5cm，初始水深15cm。采用VSF软件中的湖泊包(Lake package)模拟单环水位变化，计算间隔设置为30s，单环水深等于0cm时计算终止。

表2为以数值模拟得到的单环水位变化为目标，MN方法、一维变水头GA方法、Philip方法和Nestingen方法推求的C、K结果。与表1设定的Carsel值相比较，本发明(MN方法)计算垂向饱和渗透系数精度最高，但各种方法推求的C值误差都较大。这是由于上述各种方法都基于GA模型假设：土壤水力特征函数的渐变过程(如VGM模型)概化为脉冲函数(δ函数)，这种概化会导致土壤吸力的估算误差。

表2基于数值模拟水头变化过程参数率定结果

表中：C的单位为cm；K的单位为cm/d。

单环降水头入渗方法应用

利用河海大学水文水资源及水利工程国家重点实验室降雨大厅可变坡土槽及大型斜坡土槽，分别做两次单环变水头入渗试验。可变坡土槽长12m、宽1.5m、高1.5m，土槽内土壤于2004年填充，历经几十次降水入渗试验，土壤特性接近自然状态；大型斜坡土槽长10m、宽3.3m、高2.2m，土槽内土壤于2007年填充，经历过一次长时间降水入渗试验。两次试验采用的单环直径分别为30cm和50cm，高为20cm，插入土壤深度分别为5cm及3cm。两土槽土壤特性见表3。

单环降水头入渗试验分如下四步：(1)利用TDR观测土壤初始含水量θ₀；(2)向单环中投放HOBO U20自记水位仪，自动记录环内水深变化；(3)将单环水深瞬时加注到H₀，然后任其自由下渗；(4)单环中的水入渗完后利用TDR观测土壤含水量，此时土壤含水量即可认为是饱和含水量θ_S，土壤缺水量则为Δθ＝θ_S-θ₀。

根据实测单环水位变化数据，利用MN方法推求的土壤垂向渗透系数见表3，MN方法模拟与实测水位变化拟合结果见图3和图4。Schaap等通过统计分析大量试验结果后得出了12种典型土壤的渗透系数，其中砂壤土为1.58±0.66log(cm/d)、粉壤土为1.26±0.74log(cm/d)，这与表3中MN参数率定结果非常接近。

表3试验土槽土壤物理特性及参数推求结果

Claims (3)

1.一种用于测量土壤垂向饱和渗透系数的方法，其特征在于该方法包括以下步骤：

(1)将单环(1)缓慢均匀插入土壤，记录单环(1)插入土壤的深度L；

(2)将土壤含水量测量仪(3)插入单环(1)内的土壤中，测量记录土壤初始含水量；

(3)将自记水位计(2)放入单环(1)内，往单环(1)内注水，记录单环(1)水位变化，即时间序列{t_i}及单环内水深变化序列{H_i}；

(4)在单环(1)中水入渗结束时，用土壤含水量测量仪(3)测定单环(1)内土壤含水量；

(5)根据单环(1)变水头入渗试验中水深变化过程的计算原理，将单环(1)内水流入渗分为两个阶段，单环(1)中的一维入渗及出单环(1)后的三维入渗；单环(1)中一维入渗采用一维变水头Green-Ampt入渗模型计算：

$t=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(\frac{H_0-h}{\mathrm{\Δ\θ}}\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+(H_0-h)\×(1-\mathrm{\Δ\θ})/\mathrm{\Δ\θ}}\right)),$ 其中H₀-h≤LΔθ(1)

式中：Δθ为土壤缺水量，H₀为单环初始水深，C为土壤吸力，h为单环内水深，t为入渗时间，K为土壤垂向饱和渗透系数，L为单环插入土壤深度；

出单环(1)后的三维入渗包含两部分：重力引起的下渗土壤水以及单环水压力与土壤吸力引起的土壤水扩散，重力入渗水量等于土壤垂向饱和渗透系数乘入渗时间，土壤水扩散采用Nestingen提出的球缺扩散模型计算，此时段单环水位变化计算公式为：

$h=H_0-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{{2R}^3+3{\mathrm{LR}}^2-{2r}_0^3-{3\mathrm{Lr}}_0^3}{{3r}_1^2}\mathrm{\Δ\θ},$ 其中H₀-h＞LΔθ  (2)

式中：LΔθ为一维入渗水量，K(t-t₀)为重力引起的扩散水量，其中t₀为一维土壤水入渗结束时间，计算公式为：

$t_0=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(L\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+L\×(1-\mathrm{\Δ\θ})}\right))---\left(3\right)$

为单环水压力与土壤吸力引起的土壤水扩散，其中r₁为单环半径，R为土壤水球缺扩散半径，r₀为初始球缺半径，R可通过求解如下微分方程获得：

$\left\{\begin{array}{c}\frac{\mathrm{dR}}{\mathrm{dt}}=\frac{H_0+C-K(t-t_0)-\mathrm{L\Δ\θ}-\frac{{2R}^3+3{\mathrm{LR}}^2-{2r}_0^3-3L{r}_0^2}{{3r}_1^2}\mathrm{\Δ\θ}}{(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)+L)\frac{2R(R+L)\mathrm{\Δ\θ}}{{\mathrm{Kr}}_0^2}}\\ {R|}_{t=t_0}=r_0\end{array}\right.---\left(4\right)$

式中C为土壤吸力，计算公式为：

$C=(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R}{R+L}\frac{r_0+L}{r_0}\right)+L)\left(\right|\frac{\mathrm{dh}}{\mathrm{Kdt}}|-1)-h---\left(5\right)$

(6)整理试验数据，根据实测时间序列{t_i}及单环内水深变化序列{H_i}计算出入渗末期单环水位入渗率dH/dt|_end；然后假设土壤垂向饱和渗透系数K，再联立求解下列方程组(6)、(7)和(8)，得出土壤吸力C：

$C=(\frac{{\mathrm{\π}}^2}{8}\frac{r_0(L+r_0)}{L}\ln \left(\frac{R_{\mathrm{end}}}{R_{\mathrm{end}}+L}\frac{r_0+L}{r_0}\right)+L)({\left|\frac{\mathrm{dH}}{\mathrm{Kdt}}\right|}_{\mathrm{end}}-1)-H_{\mathrm{end}}---\left(6\right)$

$H_{\mathrm{end}}=H_0-K(t_{\mathrm{end}}-t_0)-\mathrm{L\Δ\θ}-\frac{{2R}_{\mathrm{end}}^3+{3\mathrm{LR}}_{\mathrm{end}}^2-{2r}_0^3-{3\mathrm{Lr}}_0^2}{{3r}_1^2}\mathrm{\Δ\θ}---\left(7\right)$

$t_0=\frac{\mathrm{\Δ\θ}}{K{(1-\mathrm{\Δ\θ})}^2}(L\×(1-\mathrm{\Δ\θ})+(H_0+C)\×\ln \left(\frac{H_0+C}{H_0+C+L\×(1-\mathrm{\Δ\θ})}\right))---\left(8\right)$

式中：H_end为入渗末期单环水位，t_end为入渗时间，R_end为入渗末期球缺扩散半径；

然后再将土壤垂向饱和渗透系数K、土壤吸力C及实测时间序列{t_i}代入式(1)求出一维入渗时单环内水深；再代入式(4)采用Adams-Bashforth-Moulton变阶法求解该微分方程，获得球缺扩散半径R，然后再将R代入(2)求出三维入渗时单环内水深；此时，求得的单环一维与三维入渗状态下单环内水深系列即为单环变水头入渗模型中实测时间序列{t_i}对应的单环内水深变化序列{h_i}，比较其与实测内水深序列{H_i}的差异调整K重新计算{h_i}，直到{h_i}与{H_i}的误差最小，此时的K即为土壤垂向饱和渗透系数。

2.根据权利要求1所述的一种用于测量土壤垂向饱和渗透系数的方法，其特征在于：所述单环(1)插入土壤的深度2.0～5.0cm。

3.根据权利要求1所述的一种用于测量土壤垂向饱和渗透系数的方法，其特征在于：所述往单环(1)内注水，注到5.0～20.0cm深。

Images