基于水平吸渗确定土壤水分扩散率的解析-修正法

帅嘉伟，龙 平,※，胡世丽，王观石,，杨耀杰，罗嗣海

（1. 江西理工大学土木与测绘工程学院，赣州 341000；2. 江西理工大学江西省环境岩土与工程灾害控制重点实验室，赣州 341000）

0 引 言

研究土壤水分运动对农业灌溉、水利工程和涉及水分入渗等的实际工程问题具有重要意义。土壤水分扩散率是分析土壤水分运动的关键参数之一，准确测试土壤水分扩散率是高精度模拟土壤水分运动的前提，因而，建立准确测试土壤水分扩散率的方法具有重要意义。

现有测试土壤水分扩散率的方法主要分为两种。第一种是反分析法，Simŭ̊nek等在Richards方程的基础上，提出了估算水力参数的数值反分析法，进而确定土壤水分扩散率。大量研究表明，数值反分析法估算水力参数容易出现适定性问题（由于试验误差，使得数值反分析法存在多解性的问题，即适定性），且需要事先假设水力性质的数学模型，容易引入模型误差。第二种是水平吸渗法，Bruce等结合一维水平运动的Richards方程和Boltzmann变换，提出了测试土壤水分扩散率的水平吸渗法。水平吸渗法由于试验简单、耗时短而被广泛使用。邹小阳等采用水平吸渗法分析了残膜对土壤水分扩散率的影响。姚淑霞等运用水平吸渗法研究了科尔沁沙地不同生境条件下的土壤水分扩散率。然而测试土柱瞬时剖面的含水率比较繁琐，容易出现较大误差。针对难以准确测试土柱剖面上含水率的问题，Whisler等提出在土柱上布置含水率测试仪，测试所布置位置处含水率随时间的变化，再采用水平吸渗法分析含水率随时间变化的数据，虽然能简单、准确地获得试验数据，但对于高含水率，土壤水分扩散率易受边界效应影响而出现较大误差。

为了解决改进算法易受边界效应影响的问题，本文分析一维水平入渗Richards的解析解，在此基础上，对受边界效应影响的数据进行修正，再采用水平吸渗法确定土壤水分扩散率，综合水平吸渗法的优势，提出推导土壤水分扩散率的解析-修正法，并通过数值模拟和室内试验验证方法的合理性，以期为农业灌溉等土壤水运动的模拟提供精准参数。

1 解析-修正法的提出

1.1 一维水平吸渗Richards方程

第一类边界条件下，一维无限长水平吸渗定解问题可以采用一维水平吸渗Richards方程式（1）和式（2）描述。

式中为时间，d；表示以土柱最左端为原点、水平向右为正方向的水平坐标，m；为土柱在= 0处的含水率，m/m；为初始含水率，m/m；为土壤水分扩散率，m/d；为体积含水率，m/m。

1.2 Richards方程求解

式中erfc为互补误差函数；和为待定参数，为无量纲量；为互补误差函数中自变量值，为无量纲量；为Boltzmann变换参数，=/，m/d。

采用Boltzmann变换处理式（1），可以得到式（8）。

将式（8）写成差分形式，结果如式（9）所示，结合=/，基于～或～原始数据，采用式（9）确定土壤水分扩散率的方法称为水平吸渗法。

式中为离散节点编号（= 1, 2, 3, …,-1，+1为总结点数），为无量纲量；λ为离散节点处的Boltzmann变换参数（λ=+(-1)∆，∆为的步长），m/d；D为离散节点处的土壤水分扩散率，m/d；θ为离散节点处的体积含水率，m/m；为离散节点编号（= 1, 2, 3, …,），为无量纲量；θ为离散节点处的体积含水率，m/m；λ为离散节点处的Boltzmann变换参数（λ=+(-1)∆），m/d。

根据=/，将测点处含水率随时间变化的数据（～数据）转化为～数据，～数据分为两部分，分别为受边界和不受边界影响的数据，如图1所示。计算～曲线的斜率，在未受边界影响的区域，斜率d/d随着的减小而增加，当含水率受到边界效应影响时，斜率d/d快速减小，随着的减小，d/d由增加转变为减小的转折点对应的记为，根据～曲线确定对应的含水率，即为受边界和不受边界影响的界限含水率。先采用式（9）分析未受边界影响的数据，得到[,]区间对应的土壤水分扩散率，然后采用式（7）拟合得到的～数据，确定参数和。将式（7）所示的解析解写成差分形式，结果如式（10）所示。将式（10）代入式（9）中，整理可得式（11），采用式（11）修正[,]区间对应的含水率，再采用式（9）计算修正部分含水率对应的土壤水分扩散率，这种计算土壤水分扩散率的方法称为解析-修正法，其具体流程如图2所示。

图1 受边界影响和不受边界影响数据的示意图 Fig.1 Schematic diagram of data affected and unaffected by boundaries

图2 解析-修正法流程图 Fig.2 Flow chart of analytical-modified method

1.3 Hydrus-1D模拟含水率数据

采用Hydrus-1D软件模拟水分在壤砂土（Loamy Sand，LS）、砂壤土（Sandy Loam，SL）、砂质黏壤土（Sandy Clay Loam，SCL）、壤土（Loam，LO）、粉砂黏土（Silt Loam，SIL）和粉砂（Silt，SI）六种一维水平土柱中的入渗过程。水平土柱的长度设置为0.75 m，土柱的水力性质采用van Genuchten 模型（简称VG模型）描述，VG模型所描述土壤水分扩散率如式（12）所示，六种土柱的VG模型参数的取值如表1所示。左边界（= 0）设置为定含水率边界（=，为土柱在= 0处的含水率，m/m），右边界（=，为土柱长度，m）设置为自由排水边界，初始含水率均设置为0.150 m/m。设置时间节点数量为50，记录不同时刻下，体积含水率随水平坐标的变化情况（～数据）。在= 0.30 m、= 0.45 m和= 0.60 m处布置3个体积含水率测点，监测测点处体积含水率随时间的变化（～数据）。

表1 基本参数[24] Table 1 Basic parameters

式中为土壤水分扩散率，m/d；=/[(-)]，和分别为饱和含水率和残余含水率（m/m），、和为VG模型参数（= 1-1/），为饱和渗透系数（m/d）；为标准化含水率，= (-)/(-)，m/m；为体积含水率，m/m。

2 案例分析

2.1 供试土壤

试验选择了三种质地的土壤，分别取自江西省信丰县、江西省寻乌县和福建省屏南县的离子型稀土矿的表土，深度为2.3 m，采用比重计法测得各土壤的机械组成，结果如表2所示。根据《美国制土壤质地分类标准》，信丰（Xinfeng）、寻乌（Xunwu）和屏南（Pingnan）的土壤质地分别为砂质壤土、粉砂质壤土和壤土。

表2 土壤的机械组成 Table 2 Mechanical composition of soil

2.2 室内试验

取10根内径为0.14 m、长度为0.30 m的有机玻璃管（一端为开口端、一端为闭口端）进行水分传感器标定试验，距开口端0.15 m处开设0.04 m×0.02 m的矩形孔，并贴上硅胶垫。土柱质量含水率分别设置为10%、11%、12%、13%、15%、16%、17%、18%、19%和20%，往土样中分别加入对应质量去离子水，充分搅拌均匀后装入密封箱静置48 h。设置土柱孔隙比为1.0，土壤容重为1.35×10N/m，分5层装入土样，将水分传感器垂直插入硅胶垫处，设置采样时间间隔为5 min，待传感器读数稳定后停止采集数据，取最后5个数据的平均值作为传感器最终读数，采用烘干法测试探针周围土样的含水率。将质量含水率换算为体积含水率（质量含水率与体积含水率的换算关系为=/γ，为土壤容重，N/m；γ为水容重，N/m），即可得到水分传感器读数与体积含水率的关系数据，采用合适的数学表达式拟合试验数据，确定水分传感器的标定函数。

以信丰土壤为例，阐明室内一维水平吸渗试验过程如图3所示。取内径0.14 m、长度0.80 m的有机玻璃管（一端为开口端，一端为闭口端，两端装有法兰盘），有机玻璃管闭口端部开设一定数目直径约为2 mm的小孔，在距离有机玻璃管开口端0.45 m和0.60 m处分别开设一个0.04 m×0.02 m的矩形孔，在矩形孔处的内壁贴上硅胶垫。将土壤配制成质量含水率约为7.0%的湿土壤，采用烘干法测试土壤准确的初始含水率，设置土柱柱长为0.73 m、孔隙比为1.0。首先往管内装入0.02 m厚度的粗砂，然后分5层装入土样，第一层为0.13 m，其他层为0.15 m，装入下一层前刮毛土层上表面，待填样完成后，装入0.04 m厚度的粗砂，再放置一块直径为0.13 m、厚度为0.01 m且带有一定数目小孔的有机玻璃挡板，开口端连接直角弯管。在与土柱的上端平齐的直角弯管上开设一个溢流孔，在溢流孔处装上一根溢流管。在设定位置处分别插入一个FDS-100水分传感器（邯郸开发区清易电子科技有限公司，河北邯郸，测试精度：±3%），传感器从左往右依次记为和，设置水分传感器采样时间间隔为5 min。

图3 室内一维水平吸渗试验示意图 Fig.3 Schematic diagram of indoor one-dimensional horizontal imbibition test

直角弯管内快速加水至溢流孔处，之后采用蠕动泵全速注水。待两个传感器的读数稳定后试验结束。结合水分传感器标定函数，可以得到两个测点处体积含水率随时间的变化情况。取样测定试验结束时土柱的含水率。寻乌和屏南土壤的一维水平吸渗试验与信丰的试验步骤相同。

2.3 数据分析与统计方法

3 结果与分析

3.1 基于Hydrus-1D模拟的含水率的水分扩散率计算

采用Hydrus-1D软件模拟水分在LS、SL、SCL、LO、SIL和SI六种土柱中的入渗过程，得到不同时刻的体积含水率在剖面的分布情况如图4所示。以和为待定参数，采用式（6）拟合图4所示的～数据，结果如表 3所示。对于同一种土，三个时刻的参数的拟合结果的离散程度较小，其中SI土的参数拟合结果的离散程度最大，变异系数也仅为2.31%，由此可见，对于同一种土，参数和变化较小；对于不同质地的土，参数和取值具有明显差异，参数和的取值与土壤质地相关。六种土三个时刻的决定系数均大于0.990。由此可见，式 （4）或式（6）可以较好地描述第一类边界条件下，水分在半无限长土柱中入渗的含水率随时间和坐标的变化情况，即在第一类边界条件下，式（4）或式（6）可以作为式（1）的近似解析解。

表3 参数拟合结果 Table 3 Parameter fitting results

图4 体积含水率随水平土柱上x轴坐标长度的变化 Fig.4 Variation of volumetric water content with the length ofx-axis coordinate on horizontal soil column

= 0.30 m、= 0.45 m和= 0.60 m处体积含水率随时间变化的模拟结果如图5所示。利用方程=/将图 5的～数据转化为～数据，再对～数据进行等距节点插值（节点数量为30），结果如图6所示。采用解析-修正法分析～数据，得到～数据的修正结果同样绘制于图6中，土壤水分扩散率的计算结果如图7所示。对于同一种土，半无限长土柱内的～曲线是唯一的，即相同的，对应唯一的，不随时间或坐标的影响，当时间较小时，水分还未入渗至土柱最右端，其含水率随坐标的变化情况与半无限长土柱内的相同，因此可以采用较小时～数据确定唯一的～曲线（本文取图4中第一个时刻数据，该该时刻记为）。由图6可知，修正后的～数据与处的～数据基本完全重合，说明修正方法是合理的。越靠近入水端的测点，得到的修正后的含水率范围越大（如对于LS土柱，= 0.30 m和= 0.60 m修正后的最大含水率分别为0.373和0.346 m/m），使得土壤水分扩散率计算结果的范围也更广（如图7所示）。这是由于对于同一种土，3个测点的总时长相同，越靠近入水端，得到的越小，与呈负相关关系，因而含水率范围越大。

图5 体积含水率随时间的变化情况 Fig.5 Change of volumetric water content with time

采用水平吸渗法分析图6中的数值模拟～数据，得到土壤水分扩散率的结果如图7所示。由图7可知，采用水平吸渗法分析受边界效应影响的数据，得到扩散率的计算结果明显偏离真实值（扩散率的真实值为Hydrus-1D模拟时所采用的土壤水分扩散率，即把表1的参数代入式（12）得到的结果），出现较大的误差。由此可见，采用水平吸渗法分析～数据得到的～数据具有较大的局限性，仅能较为准确地得到不受边界效应部分含水率对应的土壤水分扩散率。

图6 体积含水率随λ参数的变化情况 Fig.6 Variation of volumetric water content withλ parameter

解析-修正法的计算结果也绘制于图7中。由图7可知，相比水平吸渗法，在受边界效应影响的含水率区域，采用解析-修正法计算得到的土壤水分扩散率也基本与1∶1线重合，计算结果精度较高；整个含水率区域，采用解析-修正法得到的土壤水分扩散率与真实值的决定系数均在0.900以上，而水平吸渗法计算结果的决定系数基本在0.600以下，甚至出现负相关的情况。由此可见，解析-修正法能利用～数据准确确定土壤水分扩散率。

图7 土壤水分扩散率的解析-修正法和水平吸渗法的计算结果 Fig.7 Results of the analytical-method and horizontal imbibition method of soil water diffusivity

3.2 解析-修正法验证结果

基于室内试验，对两个水分传感器进行标定，如图8所示，线性回归方程拟合精度较高（>0.97），能够用于土壤含水率计算。基于图8中方程得到两个测点的含水率动态变化数据（～）如图9所示。利用方程=/将～数据转化为～数据（等距差值数量为15），结果如图10所示，根据～数据可以确定和，结果也列于图10中，采用式（9）分析图10中[,]区间的数据，确定相应的土壤水分扩散率D，再采用式（10）拟合未受边界效应影响的～数据，得到近似解析解参数和的值如表4所示。再采用式（11）修正图10中[,]区间对应的含水率，修正后的～曲线绘于图10中，采用式（9）分析修正后的～曲线，得到土壤水分扩散率如图11所示，水平吸渗法的结果也绘于图11中。由图11可知，当>时，水平吸渗法得到的土壤水分扩散率随着含水率的增加而减小，土壤水分扩散率随含水率变化的趋势不相符，说明不能直接采用水平吸渗法计算受边界效应影响的土壤水分扩散率。

图8 土壤水分传感器标定曲线 Fig.8 Calibration curve of soil moisture sensor

表4 近似解析解参数的拟合结果 Table 4 Fitting results of approximate analytic solution parameters

图9 室内一维水平土柱吸渗试验θ～t数据 Fig.9θ-t data of one-dimensional horizontal soil column imbibition test

图10 体积含水率θ随Boltzmann参数λ变化的试验结果与修正结果 Fig.10 Experimental and modified results of volumetric water contentθ varying with Boltzmann parameterλ

图11 解析-修正法和水平吸渗法的土壤水分扩散率计算结果 Fig.11 Results of soil water diffusivity calculated by analytic - modified method and horizontal imbibition method

对于解析-修正法的结果，当>时，土壤水分扩散率随着含水率的增加而快速增加，满足土壤水分扩散率随含水率变化的基本趋势；对于同一土柱，根据两个测点计算得到的土壤水分扩散率吻合度很好，= 0.60 m的含水率范围，土壤水分扩散率计算结果的决定系数分别为0.958、0.962和0.951，说明不同测点的计算结果具有很好的一致性。

4 结 论

1）结合线性化水平入渗Richards方程的解析解和常数变易法，推导了一维水平入渗Richards方程的近似解析解，与Hydrus-1D模拟结果的决定系数均在0.990以上，验证了近似解析解是合理的。近似解析解中的参数与土壤性质相关，对于同一种土，参数变化较小。

2）先采用近似解析解修正受边界效应影响的数据，再采用水平吸渗法确定土壤水分扩散率，建立了解析-修正法，结果发现，解析-修正法能精确地修正受边界效应影响的数据，且越靠近入水端，修正的含水率范围越大。

3）室内试验和Hydrus-1D模拟试验都表明，含水率随时间变化的数据易受边界效应的影响，直接采用水平吸渗法确定土壤水分扩散率，在高含水率区的精度较差，整体决定系数均在0.600以下，而采用解析-修正法得到的土壤水分扩散率的决定系数均在0.900以上，弥补了水平吸渗法的不足。

4）本文方法没有事先对试验数据进行处理，减少测试误差的影响，导致得到的土壤水分扩散率的计算易受测试误差的影响。进一步研究中需引入试验数据预处理算法（如移动平均法），提高测试结果的精度。