模拟非饱和土壤水分运动，可通过求解符合边界条件的Richards方程实现^[1-4]，该方法对涉及水分入渗的实际工程具有重要意义。而测定土壤水力特性是求解Richards方程的关键之一^[5-8]，其中土壤水分扩散率是模拟土壤水分水平入渗运动的重要参数，并且可作为测试土水特征曲线的中间参数，用于简化土水特征曲线测试^[9]。因此，建立准确测定土壤水分扩散率的方法具有重要意义。

测试土壤水分扩散率的方法主要分为基于水流入渗过程的反演方法和基于水平入渗试验的直接测试法。Gardner^[10]提出利用压力盘出流数据来反演土壤水力特性，其中包括反演土壤水分扩散率。反演方法通过水流入渗试验中的特定数据结合水流入渗模型能有效估计相应的土壤水力特性参数^[11-14]。但是反演方法受解的非唯一性影响，需要严格的边界条件和特定的水流入渗过程数据才能减小其影响^[15-18]。Bruce和Klute^[19]提出利用一维水平试验测试土壤水分扩散率的水平入渗法，该方法操作简单，但是在收集试验数据时容易出现较大的误差^[20-23]。Evangelides等^[24]通过记录湿润锋随时间运动情况，结合累计入渗量，利用其构造的函数计算土柱剖面含水率分布，最终计算出土壤水分扩散率，该方法避免了测试土壤含水率的误差，但并没有较可靠的记录湿润锋位置变化的方法，可能会影响结果的准确性。Villarreal等^[25]在试验土柱上布置水分传感器记录特定位置处含水率随时间变化数据，以替代土柱剖面含水率分布数据，但由于存在边界效应，该方法的传感器摆放位置会对结果造成一定影响。为减小边界效应影响，帅嘉伟等^[26]提出求解一维水平入渗Richards方程在第一类边界条件下的近似解析解，基于近似解析解，修正受边界影响的数据，再采用水平入渗法确定土壤水分扩散率，虽然该方法能降低边界效应影响，但是受含水率测试误差的影响较大。

为了降低测试误差对水平入渗法测试结果的影响，本文结合Van Genuchten(VG)模型^[27]和Philip法^[28]，先预设模型参数，采用Philip法计算不同体积含水率对应的Boltzmann参数^[19]，量化计算值与测试值的误差，再通过优化算法确定模型参数，以此建立具有抵抗测试误差能力的迭代法确定土壤水分扩散率，并通过数值模拟和文献数据验证该方法的合理性，以期为涉及水分入渗的实际模拟提供准确参数。

1 研究方法与数据来源 1.1 一维水平入渗基本方程

一维水平入渗可以采用一维水平Richards方程描述，如式(1)所示；对于第一类边界条件下，无限长水平入渗的定解条件如式(2)所示^[29-30]。

$ \frac{{\partial \theta }}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {D\frac{{\partial \theta }}{{\partial x}}} \right) $

(1)

$ \left\{ \begin{array}{l} \theta \left( {x = 0, t \geqslant 0} \right) = {\theta _{\text{0}}} \hfill \\ \theta \left( {x > 0, t = 0} \right) = {\theta _{{\text{ini}}}} \hfill \\ \theta \left( {x \to \infty , t \geqslant 0} \right) = {\theta _{{\text{ini}}}} \hfill \end{array} \right. $

(2)

式中：t为时间，d；x表示水平坐标，本文以土柱最左端为原点、水平向右为正方向，m；θ为体积含水率，m³/m³；D为土壤水分扩散率，是体积含水率的函数，D=D(θ)，m²/d；θ₀和θ_ini分别为进水端边界土壤含水率和初始含水率，m³/m³。

第一类边界条件下，采用Boltzmann变换处理式(1)，可以得到式(3)^[31]：

$ \int_{{\theta _{{\text{ini}}}}}^\theta {\lambda {\text{d}}\theta } = - 2D(\theta )\frac{{{\text{d}}\theta }}{{{\text{d}}\lambda }} $

(3)

式中：λ为Boltzmann变换参数，λ = x/t^0.5，m/d^0.5。

1.2 Philip法

式(3)给出了λ随θ变化的隐式关系，当给定土壤水分扩散率D时，为了计算不同θ对应的λ，Philip提出了一种求解方法，具体计算步骤如下^[31]：

1) 将含水率区间[θ₀, θ_ini]等分为n份，分别记为θ₀、θ₁、θ₂、…、θ_i、…、θ_n，其中θ₀ = θ₀，θ_n = θ_ini，步长Δθ = (θ₀−θ_ini)/n，将各个体积含水率代入D – θ的关系式中，即可得到D_i(i = 0, 1, 2, …, n)，将D_i代入式(4)中，得到土壤水分扩散率的加权平均值$\bar D$；采用线性差分计算D_i+1/2的值(式(5)，i = 0, 1, 2, …, n−1)。

$ \bar D = \frac{2}{{{n^2}}}\sum\limits_{i = 0}^n {\left( {n - i} \right){D_i}} $

(4)

$ {D_{i + {1 / 2}}} = \frac{{{D_i} + {D_{i + 1}}}}{2} $

(5)

式中：n和i分别为含水率区间等分份数及编号；$\bar D$为土壤水分扩散率的加权平均值。

2) 采用式(6)计算λ₁，然后交替采用式(7)和式(8)计算I_i+1/2和λ_i(i = 1, 2, …, n−1)。

$ {\lambda _1} = {{2{D_{{1 / 2}}}} / {{I_{{1 / 2}}}}} $

(6)

$ {I_{i + {1 / 2}}} = {I_{\left( {i - 1} \right) + {1 / 2}}} - {\lambda _i} $

(7)

$ {\lambda _{i + 1}} = {\lambda _i} + {{2{D_{{{i + 1} / 2}}}} / {{I_{{{i + 1} / 2}}}}} $

(8)

式中：I_i+1/2为简化计算引入的迭代变量。

3) 采用式(9)计算I_i+1/2的迭代初始值I_1/2。

$ {I_{{1 / 2}}} = 2n\sqrt {{{\bar D} / {\rm{ \mathsf{ π}}}}} $

(9)

4) 采用式(10)计算参数y，采用式(11)计算参数A；再分别采用式(12)和式(13)计算$I_{n - {1 / 2}}^{\text{F}}$和$I_{n - {1 / 2}}^{\text{S}}$的值，计算$I_{n - {1/ 2}}^{\text{F}}$和$I_{n - {1 / 2}}^{\text{S}}$的相对误差(式(14))。

$ y = \frac{{{\lambda _{n - 1}}}}{{2\sqrt {{D_{n - {1 / 2}}}} }} $

(10)

$ A = \frac{{2y{{\text{e}}^{ - {y^2}}}}}{{\sqrt {\rm{ \mathsf{ π}}} {\text{erfc}}\left( y \right)}} - 2{y^2} $

(11)

$ I_{n - {1 / 2}}^{\text{F}} = {I_{\left( {n - 2} \right) + {1 / 2}}} - {\lambda _{n - 1}} $

(12)

$ I_{n - {1 / 2}}^{\text{S}} = \frac{{{\lambda _{n - 1}}}}{2} + \frac{{2{D_{n - {1 / 2}}}A}}{{{\lambda _{n - 1}}}} $

(13)

$ \Delta = \left| {I_{n - {1 / 2}}^{\text{F}} - I_{n - {1 / 2}}^{\text{S}}} \right| $

(14)

式中：Δ为迭代误差；$I_{n - {1 / 2}}^{\text{F}}$和$I_{n - {1 / 2}}^{\text{S}}$分别为迭代值和相对精确值；y和A为简化计算引入的迭代变量。

5) 给定最大允许误差值[Δ]，若Δ≤[Δ]，θ_i对应的λ_i即为所求结果；若Δ > [Δ]，采用式(15)修正I_1/2的值，重复上述步骤3)、4)和5)，直至Δ≤[Δ]。

$ {I_{{1 / 2}}} = {I_{{1 / 2}}} - \frac{{\mathit{\Delta }}}{2} $

(15)

1.3 确定土壤水分扩散率的迭代法

1) 试验数据为θ – x数据。某一时刻，得到了水平入渗土柱上不同位置的体积含水率(即θ – x数据)，采用λ = x/t^0.5，将θ – x数据转化为θ – λ数据，λ_{test, 0}、λ_{test, 1}、…、λ_{test, k}、…、λ_{test, N}对应的体积含水率分别为θ_{test, 0}、θ_{test, 1}、…、θ_{test, k}、…、θ_{test, N}，测试土柱的初始含水率θ_ini、稳定含水率θ₀和饱和含水率θ_s。采用VG模型量化土壤水分扩散率随体积含水率的变化规律(式(17))，给定参数[θ_r, m, D_s]，设定n和[Δ]的值，采用1.2节介绍的Philip法，即可得到θ₁、θ₂、…、θ_i、…、θ_n对应λ的计算值，分别记为λ_{cal, 1}、λ_{cal, 2}、…、λ_{cal, i}、…、λ_{cal, n}。同时，当θ = θ_0­时，对应λ的值为0，则λ_{cal, 0} = 0；当θ = θ_ini­时，对应λ的值为∞，为了避免Philip法的计算结果出现∞，采用试验数据中λ的最大值λ_max代替λ_{cal, n}。采用插值算法处理θ – λ数据的计算值，得到λ_{test, 0}、λ_{test, 1}、…、λ_{test, k}、…、λ_{test, N}对应体积含水率的计算值θ_{cal, 1}、θ_{cal, 2}、…、θ_{cal, k}、…、θ_{cal, N}。以体积含水率的试验值和计算值的平方和作为计算误差(式(18))，采用内点法搜索计算误差最小时对应的参数[θ_r, m, D_s]，将得到的参数代入式(17)，即可得到土壤水分扩散率D(θ)。由于采用的Philip法是通过迭代的方式计算λ，因而本文将提出的方法称为迭代法。

$ D = \frac{{{D_{\text{s}}}{{\left[ {1 - {{\left( {1 - {{\mathit{\Theta}} ^{1/m}}} \right)}^m}} \right]}^2}}}{{{{\mathit{\Theta}} ^{\left( {m + 2} \right)/\left( {2m} \right)}}{{\left( {{{\mathit{\Theta}} ^{ - 1/m}} - 1} \right)}^m}}} $

(17)

$ \xi = \sum\limits_{k = 1}^N {{{\left( {{\theta _{{\text{test}}, k}} - {\theta _{{\text{cal}}, k}}} \right)}^2}} $

(18)

式中：Θ为标准化含水率，Θ = (θ–θ_r)/(θ_s–θ_r)，θ_r和θ_s分别为残余和饱和含水率，m³/ m³；D_s = K_s/[nmα (θ_s–θ_r)]，m²/d，其中n、m和α为土柱孔隙分布相关的参数，m = 1–1/n，K_s为饱和渗透系数，m/d；θ_test和θ_cal分别体积含水率的试验值和计算值，m³/ m³；ξ为体积含水率计算误差。

2) 试验数据为θ – t数据。在土柱上布置监测点，监测该位置体积含水率随时间的变化(即θ – t数据)，采用λ = x/t^0.5，将θ – t数据转化为θ – λ数据，采用文献[26]的方法确定受边界效应影响的临界含水率θ_lim，以[θ_ini, θ_lim]范围内的θ – λ数据作为基础数据。之后，计算土壤水分扩散率的过程与试验数据为θ – x数据时相同，不再赘述。

1.4 数据来源

通过数值模拟数据及他人一维水平入渗试验数据验证本文方法，其中，数值模拟数据通过Hydrus-1D软件获取，他人试验数据引用Villarreal等^[25]、帅嘉伟等^[26]和雷志栋^[31]研究中的一维水平入渗试验数据。利用数值模拟及他人试验中的θ – x数据和θ – t数据确定一维水平入渗率，再与对应一维水平入渗率的真实值比较，验证方法的可靠性。

1.5 Hydrus-1D数值模拟试验

采用Hydrus-1D软件模拟水分在壤砂土、砂壤土、砂质黏壤土和壤土4种土柱中的一维水平入渗过程。土柱的长度设置为1 m，土壤水分扩散率与体积含水率的关系采用VG模型描述，4种土柱的VG模型参数如表 1所示^[27]。左边界(x = 0)设置为定含水率边界(θ = θ₀)，右边界(x = L，L为土柱长度，m)设置为自由排水边界，初始时刻土柱的含水率为θ_ini，4个土柱θ₀和θ_ini的取值也列于表 1。在x = 0.25 m处布置测点，监测测点处体积含水率θ随时间t的变化(θ – t数据)。

2 结果与分析 2.1 测试误差的影响

4种土体分别在t=0.6、1.6、12和6.0 d时，体积含水率随坐标变化的Hydrus-1D模拟结果如图 1所示。模拟结果通过式(19)引入随机误差，随机误差分别为0%、2%、4%、6%、8% 和10% 的结果也绘制于图 1。

$ {\theta _\varepsilon } = \left[ {1 - 2\left( {{\text{rand}} - 0.5} \right)\varepsilon } \right]\theta $

(19)

式中：rand为取[0, 1]范围内的随机数运算；ε为随机误差；θ_ε为带有随机误差的体积含水率，m³/m³。

采用上述建立的迭代法分析图 1所示的数据，得到[θ_r, m, D_s]结果如表 2所示，将参数[θ_r, m, D_s]代入式(17)中，即可得到土壤水分扩散率的计算值，以理论值为横坐标，以计算值为纵坐标，展示土壤水分扩散率计算值与理论值的差异，结果如图 2所示。由图 2可知，不同测试误差的计算结果均分布在1︰1线附近，采用决定系数(式(20))量化计算值与理论值的误差，决定系数的结果也列于表 2。计算结果表明，当测试误差在[0, 10%]范围内时，决定系数均大于0.990 0，由此可见，迭代法的计算结果受体积含水率测试误差的影响较小。

$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^N {{{\left( {{D_{{\text{theory, }}i}} - {D_{{\text{cal, }}i}}} \right)}^2}} }}{{\sum\limits_{i = 1}^N {{{\left( {{D_{{\text{theory, }}i}} - \frac{1}{N}\sum\limits_{i = 1}^N {{D_{{\text{theory, }}i}}} } \right)}^2}} }} $

(20)

式中：R²为决定系数；D_theory为土壤水分扩散率的理论值；D_cal为土壤水分扩散率的计算值；i为体积含水率的编号；N为总体积含水率数据个数。

由表 2可知，对于壤土，不同的随机误差情况下，θ_r、m和D_s的计算结果分别在[0.000, 0.062]、[0.253, 0.448]和[0.188, 0.838]范围内变化，θ_r、m和D_s的变异系数CV(标准差与均值之比)分别为70.1%、19.8%和45.3%，由此可见，m值计算结果的波动最小，θ_r计算结果的波动最大。对于m值，CV∈[16%, 35%]，说明m的计算结果具有中等变异性，对于θ_r和D_s，CV > 35%，说明θ_r和D_s的计算结果具有高度变异性。虽然3个参数计算结果的变异较大，但是土水特征曲线计算结果的决定系数大于0.990 0，说明总体计算误差较小，以下通过响应曲面法分析原因。

Hydrus-1D中，壤土VG模型参数[θ_r, m, D_s]的理论值为[0.078, 0.359, 0.352]^[27]，将[θ_r, m, D_s]的理论值代入式(17)中，即可得到土壤水分扩散率的理论值D_theory。同时改变两个参数的值，另一个参数取理论值，通过式(17)得到这种参数组合下土壤水分扩散率D，采用决定系数量化D与D_theory的差异，则可以得到θ_r–m–R²、θ_r–D_s–R²和D_s–m–R²的响应曲面，结果如图 3所示。由图 3A和图 3B可知，响应曲面的等值线几乎相互平行，表明决定系数的变化对于θ_r不敏感，因而θ_r计算结果的变异系数较大，但是m和D_s取合适值时，不同的θ_r均能达到较大的决定系数。由图 3C可知，当给定θ_r时，决定系数为0.995 0对应的[m, D_s]参数组合呈月牙状，即在该范围的[m, D_s]参数值均能得到较为准确的土壤水分扩散率，说明参数[m, D_s]取值在有较大差异时，仍可得到与真实值相比差异较小的土壤水分扩散率。

2.2 临界含水率的影响

x=0.1、0.3、0.5、0.7和0.9 m处的体积含水率随时间的变化情况如图 4所示，采用文献[26]的方法确定受边界效应影响的临界含水率θ_lim，结果也列于图 4。随着监测点越靠近出水口，临界含水率越小，受边界效应影响的体积含水率范围越大。采用上述迭代法分析[θ_ini, θ_lim]范围的θ – t数据，得到[θ_r, m, D_s]的计算值如表 3所示。仍然以壤土为例进行分析，θ_r、m和D_s的变异系数CV分别为96.7%、15.9%和40.2%，即当选择不同位置的θ – t数据分析时，θ_r和D_s具有较大的变异性。随着测点与入水口的距离增加，决定系数没有明显变化，决定系数均大于0.980 0。对于所分析的4种土体类型，当x=0.9 m时，不受边界效应影响的含水率区间分别为[0.150, 0.235]、[0.150, 0.241]、[0.150, 0.266]和[0.150, 0.281]，得到的土壤水分扩散率的决定系数分别为0.998 5、0.998 7、0.997 9和0.998 4，均大于0.990 0。由此可见，受边界效应影响的临界含水率对计算结果的影响较小。

2.3 迭代法的试验验证

为了进一步通过试验数据验证本文提出方法的合理性，利用他人研究中的一维水平入渗试验数据进行验证(图 5，4种土柱分别为CHA^[25]、DOR^[25]、XF^[26]和LEI^[31])。采用迭代法分析时所需要的参数列于表 4，[θ_r, m, D_s]的计算结果也列于表 4。对于LEI^[31]，[θ_r, m, D_s]迭代法的计算结果为[0.034, 0.222, 8.131]，将参数代入式(17)，得到土壤水分扩散率如图 6所示，与水平吸渗法的结果进行对比，决定系数达到了0.891 5，将迭代法得到的土壤水分扩散率作为已知数据，采用Philip法计算λ随θ变化，结果如图 5D所示，与试验数据的决定系数为0.975 7。对于CHA^[25]和DOR^[25]，土壤水分扩散率的计算结果分别如图 5A和图 5B所示，与试验数据进行比较，本文方法的决定系数分别为0.850 0和0.962 0，文献中经验模型的决定系数分别为0.500 3和0.713 3^[25]，本文方法优于经验模型，且对于体积含水率大于0.350 m³/m³对应土壤水分扩散率的结果，经验模型的决定系数分别仅为0.148 8和0.405 9，而本文方法分别为0.746 7和0.930 4，由此可见，本文方法可以弥补经验模型难以准确确定高体积含水率区间的土壤水分扩散率的不足。对于XF^[26]，测试数据为x = 0.45 m处的θ – t数据，临界含水率为0.321 m³/m³，[θ_r, m, D_s]迭代法的计算结果为[0.000, 0.218, 10.650]，将参数代入式(17)，得到土壤水分扩散率如图 6C所示，采用式(1)模拟一维水平入渗，得到x = 0.45 m处θ – t数据的计算值，并转化为θ–λ数据，结果如图 5C所示，其与试验值的决定系数为0.993 4。

综上所述，采用迭代法分析试验数据，土壤水分扩散率的决定系数均在0.850 0以上，具有较高的精度，与Villarreal等^[25]和帅嘉伟等^[26]提出的方法相比，迭代法可以提高高含水率区域土壤水分扩散率计算结果的精度，同时受测试误差影响较小，迭代法具有较高的合理性。

3 结论

本文结合VG模型和Philip法，先预设模型参数，采用Philip法得到θ – λ数据的计算值，再以决定系数量化计算值与测试值的误差，通过内点法优化模型参数，以此建立具有抵抗测试误差能力的迭代法，并通过数值模拟和文献数据验证了该方法的合理性，得到以下结论：

1) 采用迭代法分析测试误差0% ~ 10% 的θ – x数据，θ_r和D_s的变异系数大于35%，m的变异系数∈[16%, 35%]，说明θ_r和D_s具有高度变异性，m具有中等变异性，将[θ_r, m, D_s]的结果代入VG模型，得到土壤水分扩散率计算结果的决定系数均大于0.995 0。由此可见，土壤水分扩散率的计算结果精度较高，且受测试误差的影响较小。

2) 对于θ – t数据，随着与入水口的距离增加，临界含水率θ_lim减小，不受边界效应影响的含水率范围[θ_ini, θ_lim]越窄。采用迭代法分析不同位置处[θ_ini, θ_lim]范围的θ – t数据，确定土壤水分扩散率，结果表明，土壤水分扩散率计算结果的决定系数均大于0.990 0，说明临界含水率θ_lim对土壤水分扩散率结果的影响较小。

3) 与现有确定土壤水分扩散率的方法进行对比，当分析θ – x数据时，迭代法与现有方法的的结果相近，均能得到精度较高的结果；当分析θ – t数据时，本文方法能通过确定后的参数计算出全含水率范围内的土壤水分扩散率，且在低含水率时结果与解析–修正法相近，说明本文方法可以弥补现有方法难以准确确定高体积含水率区间的土壤水分扩散率的不足，这进一步论证了迭代法的合理性。