基于保水剂应用的土壤水分运移模型构建方法及源汇项求解

廖人宽，杨培岭(.中国水利水电科学研究院，北京 00048；. 中国农业大学，北京 00083 )

0 引 言

通过建立模型来描述根系对于土壤水分的吸收动态以及土壤水分的运移规律是一种较为常用的研究方法[1-3]。目前，绝大多数水分运移模型都是基于Richards方程建立的，而考虑根系吸水的土壤水分运移模型通常是在Richards方程中加上一个考虑根系吸水的源汇项(根系吸水速率)来构建的[4-5]。由于根系吸水的不可实测性，Richards方程根系吸水源汇项反求估算值作为比较可靠的根系吸水“实测值”在根系吸水、吸肥规律研究方面得到了很大的认可。前人所提出的基于两个连续含水率剖面的根系吸水反求解法在小麦根系吸水规律研究上得到了很好的应用[6,7]。

在土壤中应用保水剂来保蓄土壤水分是旱地农业节水技术中一项重要的非工程技术方法。目前，该项技术已经在苹果、葡萄、小麦、玉米等作物上进行了广泛应用[8-10]，并取得了显著的节水和增产效果。在相当一部分旱地农业区，由于水源条件所限，其作物生长完全依赖于降雨的水分补给。然而，旱地农业区降雨时空分布不均和土壤保墒能力差导致雨水的利用效率很低，直接影响了农作物的生长[11,12]。研究表明，在土壤中施入保水剂能够有效提高旱地农业区雨水利用效率，其可以反复吸持进入土壤中的降雨，减少水分深层淋溶渗漏量，所吸持的水分随着土壤干旱程度加剧及作物根系根压增大会缓慢释放以供作物生长所需，从而保证了作物生育期内的持续供水[13,14]。

然而，大量的研究还显示[15-17]，保水剂的应用不仅仅起到吸水保水持续供给作物生长所需水分的作用，其自身独特的反复吸水-释水功能对土壤理化性质、微生物群落及作物根系生长等也都会产生一定的影响，进而对水分在土壤中的运移规律产生周期性和非稳定性影响，这也让保水剂应用条件下的土壤水分运移模型构建变得复杂和困难。韩玉国等[18,19]将水分运移时间t和保水剂施入土壤中的时间T视为两个独立的时间，建立了保水剂应用条件下的一维土壤水分运移模型并获取了解析解，但该模型只针对裸土条件下土壤水分运移的求解，并没有考虑根系吸水源汇项求解的问题，不利于深入开展保水剂应用条件下水分管理制度的研究。

基于此，本文以Richards水分运移方程为基础，将水分运移时间和保水剂施入土壤中的时间视为同一时间t来建立考虑根系吸水源汇项的水分运移模型，并通过反求迭代法来求解根系吸水速率的数值解，为进一步探索和建立保水剂应用条件下的根系吸水模型提供指导。

1 保水剂应用条件下的土壤水分运移模型

在保水剂应用条件下，因为存在保水剂保水作用效果随时间变化的因素，土壤水分非饱和扩散率D(θ)和非饱和导水率K(θ)不再只是一个随含水率θ的改变而变化的函数，而是既随含水率θ变化，又随时间t变化，含有时间项t的时变关系函数式D(θ,t)和K(θ,t)。因此，在保水剂应用条件下，考虑根系吸水及土面蒸发的一维垂直非饱和土壤水分运动的定解问题应描述为：

(1)

式中：D(θ,t)为非饱和扩散率时变函数；K(θ,t)为非饱和导水率时变函数；S(z,t)为根系吸水源汇项；E(t)为土面蒸发强度。

2 保水剂应用条件下的土壤水分运移模型的求解

对于(式1)，由于其为非线性方程，一般在进行解析求解时都是先用一个平均扩散度D来近似代替Richards方程中的非饱和扩散率D(θ,t)，再通过拉氏变换和分部积分法进行求解。然而，由于水力参数的动态变化规律是体现保水剂作用效应的重要特征，若对其进行平均简化处理将使得所建立的保水剂应用条件下的水分运移模型的意义大打折扣。因此，考虑到模型正确物理意义的体现和目前对于求解水力参数为时变函数的Richards方程解析解的巨大难度，本文采用可以兼顾体现模型正确物理意义且容易求解的隐式差分方法[20]进行数值求解。相较于显式差分和中心差分格式，隐式差分法在求解时是无条件收敛的。而且，本文数值试验中涉及对时间步长和距离步长取值等问题的探讨，采用显式差分或中心差分格式的话需要时间步长和距离步长满足一定的要求才能确保方程解的收敛性，会限制后文数值试验的探讨范围。

此外，在建立保水剂应用条件下的一维垂向非饱和土壤水分运移模型时，假设施入保水剂的土壤混合层为一种性质随时间发生动态变化的刚性土壤，忽略体积改变带来的影响；同时，忽略土壤水力参数测定过程的时间对参数变化的影响。

3 保水剂应用条件下的水分运移模型根系吸水源汇项求解

根系吸水源汇项求解借鉴Zuo Q等[6,7]提出的迭代反求法，即已实测得两个连续的含水率剖面分别为θ(z, 0)和θ(z,T)，假定期间由于根系吸水土壤各深度含水率的减少量为Δθ(z, 0～T)，T表示连续两次实测值的时间间隔，则0～T时段内的平均根系吸水速率可概化表示为S(z,0～T)=Δθ(z,0～T)/T。

首先令S=0，根据初始的含水率实测值(t=0时刻)，采用隐式差分法求解(式1)，算出t=T时刻第1次迭代的θ(1)(zi,T)值(第1次迭代时未考虑S的影响，即设S=0)，把此时计算出的Δθ(1)(zi,0～T)=θ(1)(zi,T)-θ(T)(zi,T)除以时间T后近似为平均根系吸水量S(1)=Δθ(1)/T并带入原方程继续进行求解，求解出第2次迭代出的θ(2)(zi,T)值(此后的迭代均考虑S的影响，即设S≠0)，此时进行判断Δθ(2)(zi,0～T)=θ(2)(zi,T)-θ(T)(zi,T)是否满足<误差控制标准(相对误差平方和的平均值)，若满足则结束迭代，若不满足则有S(2)=S(1)+Δθ(2)/T，将其继续带入原方程进行迭代求解直至满足迭代控制标准为止。

4 反求估算方法稳定性分析

4.1 数值试验理论根系吸水模型设定

(1)理论根系吸水模型。理论根系吸水模型采用邵爱军等[21,22]提出的模型，其模型主要为基于蒸腾作用与作物根系吸水关系而建立。拟采用的根系吸水模型表达式为：

S=ETA[e-B(lnZ-C)2]/Z

(2)

式中：ET为腾发量，mm/d；Z为相对深度，Z=z/Lr(t)；A，B，C为经验系数，A量纲为1/mm，B，C无量纲；Lr(t)为作物吸水层深度，cm，其表达式为：

Lr(t)=150(1-e-1.380 813 T)

(3)

式中：T为相对时间，T=t/M；M为玉米生育期天数，d；t为播种时间，d。

(2)数值试验参数选择。土壤物理参数及根系参数也选自邵爱军等的文献中。已知试验区田间土壤分为两层，其中0～70 cm为沙壤土，70 cm以下为粉砂土，0～120cm区域为土壤非饱和带。沙壤土非饱和导水率为D=0.001 48 e15.99θ[图1(a)]，水分特征曲线为θ=0.48 e-0.001 43|h|(0≤|h|≤129.2)，θ=3.915|h|-0.475(|h|>129.2)；粉砂土非饱和导水率为D=54.9θ3.614[图2(a)]，水分特征曲线为θ=0.50e-0.001 46|h|(0≤|h|≤<168.64),θ=13.893|h|-0.719(|h|>168.64)。根据文献中已知的土壤水分特征曲线和θ=θr+(θs-θr)/[1+(αh)n]m分别获取了两种质地土壤的非饱和导水率公式[图1(b)和图2(b)]。

理论根系吸水模型中的参数A=8.600 827×10-4-1.189 26×10-3(T-0.591)2，B=1.662 603，C=-1.308 06；M=88 d，ET本文设置为3 mm/d，迭代控制误差为1×10-5。

图1 沙壤土水力特征参数曲线图

图2 粉砂土水力特征参数曲线图

4.2 数值试验方案及验证步骤

本文根系吸水反求估算方法稳定性验证主要考虑的影响因素有：①时间间隔和时间步长；②土壤水力特征参数；③土壤含水率测量误差和仪器精度；④土壤层状性及边界条件。数值试验参数设置见表1。数值试验流程如下，其验证的步骤为：

(1)首先输入水分运移参数及理论根系吸水模型参数，并确定时间及空间步长、仪器精度、边界条件等反求方法的控制条件；

(2)输入初始含水率剖面θ(z, 0)和理论根系吸水值S(z, 0)，采用隐式差分法求解水分运移方程，得到时刻t土壤剖面的水分分布理论值θ(z,t)，由于土壤含水率实际测试过程中往往存在误差，因此本文采用一个服从正态分布的误差向量(vector error-VE)来对其进行误差扰动(程序中通过per值的大小来体现)，使其更加接近实际的情况，即θ*(z,t)=θ(z,t)+误差向量；

(3)采用代数多项式方程(式4)对第(2)步得到的θ*(z,t)进行拟合，从而获得连续而平滑的θ**(z, t)分布曲线;

θ**(z,t)=r1+r2(z-Z)+r3(z-Z)2+…+rn(z-Z)n-1

(4)

式中：r1，r2，…，rn为拟合参数；Z为z1，z2，…，zn累加值的均值，拟合从n=3开始，直至计算值和理论值之间的最大绝对误差小于误差控制要求w(w为仪器测试精度)为止；

(4)采用提出的反求解法，结合第(3)步得到的θ*(z,t)分布曲线计算不同时刻t1(初时)和t2(末时)间的平均根系吸水模拟值；

(5)分析平均根系吸水理论值同模拟值之间的误差。

4.3 数值试验结果

数值试验测试结果发现，吸水间隔、水力参数、测试误差和仪器精度对该反求方法稳定性的影响较大，时间步长、层状土壤及边界条件对模型稳定性的影响较小，从稳定性、计算量、精确度等方面综合考虑，建议本文所提出的反求方法适用条件为时间间隔5～17 d，步长为1 000～10 000，误差控制在0.9以内，仪器精度控制在0.03以内，土面蒸发强度控制在0.6 mm/d以内。

表1 各数值试验参数选取

5 结 论

本文将保水剂应用条件下的土壤水分非饱和扩散率及非饱和导水率分别表示为D(θ,t)和K(θ,t)，且本文将保水剂施入土壤中的时间和试验进行时间视为同一时间t，解决了之前提出的模型中将保水剂施入土壤中的时间和试验进行时间视为相对独立而阻碍模型进一步发展的问题，本文提出的模型构建思路和数值求解方法为建立保水剂应用条件下考虑作物根系吸水的水分运移模型提供了基础条件；提出了保水剂应用条件下土壤水分运移模型根系吸水源汇项的反求估算方法，并设置了数值试验对该反求方法进行了检验，结果发现吸水间隔、水力参数、测试误差和仪器精度对该反求方法稳定性的影响较大，时间步长、层状土壤及边界条件对模型稳定性的影响较小。

□

[1] Lekakis E H, Antonopoulos V Z. Modeling the effects of different irrigation water salinity on soil water movement, uptake and multicomponent solute transport[J]. Journal of Hydrology, 2016,530:431-446.

[2] Williamson T N, Lee B D, Schoeneberger P J, et al. Simulating soil-water movement through loess-veneered landscapes using nonconsilient saturated hydraulic conductivity measurements[J]. Soil Science Society of America Journal, 2014,78(4):1 320-1 331.

[3] Draye X, Yangmin K, Lobet G, et al. Model-assisted integration of physiological and environmental constraints affecting the dynamic and spatial patterns of root water uptake from soils[J]. Journal of Experimental Botany, 2010,61(8):2 145-2 155.

[4] Javaux M, Schröder T, Vanderborght J, et al. Use of a three-dimensional detailed modeling approach for predicting root water uptake[J]. Vadose Zone Journal, 2008,7(3):1 079-1 088.

[5] Li H, Yi J, Zhang J, et al. Modeling of soil Water and salt dynamics and its effects on root water uptake in Heihe arid wetland, Gansu, China[J]. Water, 2015,7(5):2 382-2 401.

[6] Zuo Q, Zhang R. Estimating root-water-uptake using an inverse method[J]. Soil Science, 2002,167(9):561-571.

[7] Zuo Q, Meng L, Zhang R. Simulating soil water flow with root-water-uptake applying an inverse method[J]. Soil Science, 2004,169(1):13-24.

[8] Eneji A E, Islam R, An P, et al. Nitrate retention and physiological adjustment of maize to soil amendment with superabsorbent polymers[J]. Journal of Cleaner Production, 2013,52(4):474-480.

[9] Namvar H, Asgharzade A, Babaeian M, et al. The effects of superabsorbent polymer on yield and yield component of two grape varieties.[J]. International Journal of Biosciences, 2014.

[10] 庄文化, 吴普特, 冯 浩,等. 土壤中施用聚丙烯酸钠保水剂对冬小麦生长及产量影响[J]. 农业工程学报, 2008,24(5):37-41.

[11] Turner N C. Agronomic options for improving rainfall-use efficiency of crops in dryland farming systems.[J]. Journal of Experimental Botany, 2004,55(407):2 413-2 425.

[12] McCutcheon, M. C, Farahani, et al. Effect of soil water on apparent soil electrical conductivity and texture relationships in a dryland field[J]. Biosystems Engineering, 2006,94(1):19-32.

[13] 廖人宽, 杨培岭, 任树梅. 高吸水树脂保水剂提高肥效及减少农业面源污染[J]. 农业工程学报, 2012,28(17):1-10.

[14] Liu F, Ma H, Xing S, et al. Effects of super-absorbent polymer on dry matter accumulation and nutrient uptake of Pinuspinaster container seedlings[J]. Journal of Forest Research, 2013,18(3):220-227.

[15] Han Y G, Yang P L, Luo Y P, et al. Porosity change model for watered super absorbent polymer-treated soil[J]. Environmental Earth Sciences, 2010,61(6):1 197-1 205.

[16] Guilherme M R, Aouada F A, Fajardo A R, et al. Superabsorbent hydrogels based on polysaccharides for application in agriculture as soil conditioner and nutrient carrier: A review[J]. European Polymer Journal, 2015.

[17] Bai W, Zhang H, Liu B, et al. Effects of super-absorbent polymers on the physical and chemical properties of soil following different wetting and drying cycles[J]. Soil Use & Management, 2010,26(3):253-260.

[18] Han Y G, Yu XX, Yang P L, et al. Dynamic study on water diffusivity of soil with super-absorbent polymer application[J]. Environmental Earth Sciences, 2012, 69(1): 289-296.

[19] 韩玉国,武亨飞,杨培岭,等. 保水剂对土壤的物理性质与水分入渗的动态影响[J]. 干旱地区农业研究, 2013,31(5):161-167.

[20] 雷志栋, 杨诗秀, 谢森传. 土壤水动力学[M]. 北京：清华大学出版社, 1988.

[21] 邵爱军. 土壤水分运动的数值模拟——以作物根系吸水项为例[J]. 水文地质工程地质, 1996,(2):5-8.

[22] 邵爱军, 李会昌. 野外条件下作物根系吸水模型的建立[J]. 水利学报, 1997,(2):68-72.