地膜覆盖条件下SPAC系统水热耦合运移模型的研究

吴从林,黄介生，沈荣开
(武汉水利电力大学水利系)

摘要：本文建立了地膜覆盖条件下SPAC系统水热耦合运移模型。该模型由大气、植物冠层、地膜覆盖层和土壤四层组成。根据能量平衡原理建立了各层的能量平衡方程。采用牛顿-莱普森方法和有限差分法求解了各能量平衡方程和土壤水分及温度剖面，并用冬小麦的实测资料对模型进行了验证，结果表明模型具有一定的精度。此外，本文对SPAC系统中的阻力项进行了探讨。

关键词：地膜覆盖；SPAC系统；腾发量；水热运移

收稿日期：2000-02-22
作者简介：吴从林(1975-)，男，研究生。

　　自地膜覆盖技术出现以来，世界各国对其增温保墒效应及其机理开展了广泛的研究。但以往的研究更多的是偏重于试验，而对其增温保墒机理的研究相对较少。1979年，以色列学者Maherer等人首次从理论上对这一问题进行了研究，他们建立了透明塑膜覆盖层的热平衡方程，从物理的角度描述了大气、覆盖层、土壤的热交换过程，为探索地膜覆盖的机理奠定了基础；1984年，Maherer等又将Philip de Vries水热耦合方程应用于透明塑膜覆盖条件下的土壤水热分布问题，将地膜覆盖机理的研究又向前推进了一步。之后，Chung和Horton于1987年又针对覆盖条件下的二维土壤水热耦合运动进行了数值计算。国内王树森和邓根云于1988年对地膜覆盖后土壤能量平衡的变化特点进行了理论分析。1992年，隋红建等将油纸，塑膜及土面增温剂等均视为薄膜类材料，建立了不同覆盖条件下土壤二维水热运动数学模型。但以往的研究都忽略了一个很重要的因素，那就是作物。显然，有作物生长的情况下，地膜与大气，地膜与土壤间的能量交换与没有作物生长的情况是不同的。从另一个角度来讲，地膜的存在改变了土壤-植物-大气间的能量交换和水分传输规律，同时也就改变了作物的耗水规律 ^[1]。研究有作物生长情况下地膜覆盖的增温保墒机理，探讨地膜覆盖条件下SPAC系统水热运移的规律，对于了解地膜覆盖条件下的作物耗水规律，科学地进行水管理具有十分重要的意义。

　　本文将地膜覆盖农田概化为由大气、植物冠层、地膜覆盖层和土壤组成的一维系统，根据能量平衡原理建立了各层的能量平衡方程，并与土壤水热运移方程相结合，建立了地膜覆盖条件下SPAC系统水热运移模型，同时，采用冬小麦地膜覆盖农田的实测资料对模型进行了验证。

1 地膜覆盖条件下SPAC系统水热运移模型的建立

　　地膜覆盖条件下SPAC系统的能量分配与转换可以由图1所示的四层结构来描述。第一层为位于参考高度处的大气层；第二层为位于动量传输汇处的植物冠层；第三层为地膜覆盖层；第四层为土壤层，计算深度取1m，其顶部为土壤表面，底部的土壤水分和温度由实测资料给出。为简化计，对冠层作如下假设：(1)冠层在水平方向上是足够均匀的；(2)空气中分子输送项远比湍流输送项小，可忽略不计；(3)植物冠层内植物体吸收的净辐射能将全部用于与周围空气进行潜热和显热交换，而植物体本身没有净能量贮存；(4)热量、水汽在垂直方向的输送能量可用梯度-扩散理论表示即

图1 地膜覆盖条件下SPAC系统的能量分配与转换

　　在以上简化和假定的前提下，可根据能量平衡原理和土壤水动力学中土壤水热运移的基本方程分别建立各层的能量平衡方程及土壤水热运移方程。

1.1 各层能量平衡方程 根据能量平衡原理，图1所示各层的能量平衡方程可分别表示如下 ^[2]：

R_S+R_l-ε_pσT⁴_p-H_pa-λE_p-G=0	(1)
(α_pg+r_mgt_pg)R_g+(α_pl+r_mlt_pl)R_L+(r_ml-2)ε_pσ⁴_p+ε_mσT⁴_m+t_mlε_sσT⁴_s-λE_p-H_lp=0	(2)
(α_mg+r_sgt_mg)t_pgR_g+α_ml(t_plR_l+ε_pσT⁴_p)-2ε_mσT⁴_m+α_mlε_sσT⁴_s+H_sm-H_mp=0	(3)
(1-r_sg)t_mgt_pgR_g+t_ml(t_plR_l+ε_pσT⁴_p)+ε_mσT⁴_m-(1-r_ml)ε_sσT⁴_s-H_sm-G=0	(4)

式中:R_g、R_l、R_s和R_n分别为到达地面的太阳辐射、长波辐射、短波辐射和净辐射(W/m²)；T_a、T_p、T_l、T_m和T_s分别为空气、冠层、叶片、地膜和地表温度(℃)；r_ml、r_mg和r_sg分别为地膜的长、短波反射率和土壤短波反射率；a_pg、a_pl、a_mg和a_ml分别为冠层和地膜对短波和长波的吸收率；t_pg、t_pl、t_mg和t_ml分别为冠层和地膜对短波和长波的透射率；ε_p、ε_m和ε_s分别为冠层、地膜和土壤的散射率；σ为Stefan Boltzman常数(5.67×10^-8w／m²·k⁴).

　　式(1)～(4)中各辐射项按下列公式计算^[2]：

R_s=(1-r_pg)R_g	(5)
R_l=σ(T_a+273.16)⁴[0.605+0.048(1370H_a)^1/2]	(6)
R_n=R_s+R_l-ε_pσT⁴_p	(7)
	(8)

以上各式中：H_a为空气绝对湿度(kg/m³)；k_h和T_a分别为土壤导热率(J/(s.k.m))和空气温度(℃).

　　潜热及显热项按下列公式计算：

H_pa=ρc_pT_p-T_a/r^a_a	(9)
H_lp=ρc_pT_l-T_p/r^c_a	(10)
λE=ρc_p/γe_p-e_a/r^a_a	(11)
H_mp=ρc_pT_m-T_p/r^s_a	(12)
H_sm=ρc_pT_s-T_m/r^m_s	(13)
λE_p=ρc_p/γe_sat(T_l)-e_p/r^p_s+r^c_a	(14)
λE=λE_p	(15)
H_pa=H_lp+H_mp	(16)

以上各式中：ρ为空气密度(kg/m³)，γ为比湿计常数(mb/k)，C_p为水蒸汽比热(j/kg·k),e_sat(T₁)为温度为T_l时的饱和水汽压(mba)，e_a和e_p分别为参考高度和冠层处空气的水汽压(mba)，λ为水蒸汽潜热(j/kg)，H和E分别表示各部分显热和潜热(w/m²);r^a_a、r^c_a、r^s_a、r^m_s和r^p_s分别为冠层到参考高度处空气动力学阻力、边界层阻力、冠层内空气动力学阻力、土壤与地膜之间的阻力和冠层气孔阻力(s/m).

1.2 土壤水热运移方程^[3，4] 土壤水热运移的控制方程已发展到比较完善的程度，有等温和非等温两种模型。关于等温和非等温模型的适用条件国内外都有研究。Philip的研究表明，等温模型在描述土壤未达到十分干燥的情况下的蒸发过程是比较有效的；水建高等对松土覆盖条件下土壤水热运移的模拟表明，当土壤含水率较高时使用等温模型精度较高，但当土壤含水率较低时误差比较明显^[5]。考虑到本文研究对象的土壤水分是从接近饱和到接近凋萎点的整个过程，这里采用非等温模型。由于非等温模型涉及的因素非常复杂，为简化起见，模型中将忽略某些影响较小的因子。根据Milly的研究，忽略温度梯度引起的液态水流动项q_wT造成的误差小于1.5%；忽略温度梯度造成的汽态水流动项q_vT将使计算的蒸发量增大1%～15%；而忽略水势梯度下的汽态水流动项q_vh将使蒸发量明显减少^[3]。因此，这里忽略q_wT和q_vh，可得如下水热耦合控制方程：

	(17)
	(18)

式(17)和(18)中水、热方程的系数分别为：

	(19)
	(20)
	(21)
	(22)

其中：C=C_sρ_s(1-n)+C_Lρ_lθ_l+C_aρ_aθ_v,G=(C_Lρ_l-C_aρ_s)T-Lρ_v

式中：C_s，C_l和C_v分别为土壤基质、土壤水及土壤空气的比热(j/kg.k);ρ_s是土壤基质的密度(kg/m³)；ρ_a是土壤空隙内的空气密度(kg/m³)，ρ_v是土壤空隙内的水汽密度(kg/m³)；L是蒸发潜热；h,T分别为土壤基质势和温度，θ_l和θ_v分别为土壤中液态水和汽态水体积含水率(m³/m³)；k是非饱和导水率(m/s)，D_hv和D_Tv分别为水蒸汽传导率和热蒸汽扩散率；z是垂直坐标，向下为正；S_r为根系吸水速度(m/s).

1.3 土壤水热运移方程的定解条件

初始条件：

h(z,t)｜_t=0=h(z,0) 0≤z≤100cm

T(z,t)｜_t=0=T(z,0) 0≤z≤100cm

下边界条件：

h(z,t)｜_z=L=h(L,t) t≥0,L=100cm

T(z,t)｜_z=L=T(L,t) t≥0，L=100cm

上边界条件：

(23)

当地表有地膜覆盖时

q_m｜_z=0=0

(24)

式中：qm为地表处水分通量。

1.4 SPAC系统能量方程中阻力项的确定

1.4.1 零平面位移d和粗糙高度z₀ 零平面位移d表示把冠层简化为一个大叶层时的有效高度，粗糙高度z₀表示冠层从气流中吸收动量的效率尺度^[3]，分别由下列公式计算：

d=0.719H	(25)
z₀=(z_m-d)exp[-(κU_H/U_*+1-1/G)]	(26)

式中：z_m=d+G(H-d),G=2.0～3.0，H为冠层顶高度(m)；U_H为冠层顶处风速(m/s)；κ为卡尔曼常数(0.41)；U为摩擦风速(m/s).

1.4.2 冠层叶片边界层阻力r^c_a 根据Jones的研究，冠层叶片边界层阻力r^c_a或传导率g_b是风速的函数^[3]，由边界层理论得：

r^c_a=1/LAI	(27)
=2a/n[U_H/w]^0.5[1-exp(-n/2)]	(28)

式中：w是叶片宽度(m)，a=0.01ms^-0.5，U(z)是冠层内高度为z处的风速。

1.4.3 土壤水汽蒸发扩散阻力r^s_s r^s_s是水汽从土壤中水汽源面运移到土壤表面时所受到的阻力。孙菽芬^[7]和Camillo^[8]的研究认为，土壤水分蒸发扩散阻力主要与土壤含水率有关。

　　本文采用孙菽芬的计算公式：

r_s=3.5[θ_s/θ]^2.3+33.5

(29)

1.4.4 气孔总阻力r^p_s 气孔总阻力r^p_s的大小，除由植物本身的特性决定外，还受环境因子的影响。卢振民提出一种非线性模型，他认为影响气孔阻力的环境因子有五个：净辐射R_n、空气饱和差D、CO₂浓度、叶水势φ_e、土壤水势φ_s^[9]。本文采用他提出的公式计算冠层气孔传导率和总阻力：

r^p_s=100.0/g_c

(30)

其中：g_c=[2.867LAI+0.0277(1-e^-KLAI)R_n/k][1-0.0254D]/[1+(Ψ_e/(-31529))^4.58]

式中:R_n的单位为W/m²，D的单位为hPa，φ_e的单位为cmH₂O.

1.4.5 边界层空气动力学阻力r^a_a和冠层内空气动力学阻力r^s_a 本文中为了简化起见，不考虑层结的稳定度订正，中性层结下，采用Shuttleworth和Wallace^[10]的公式：

	(31)
LAI>4.0	(32)

式中:r^a_a0,r^s_a0,r^a_a1,r^s_a1分别为裸地和冠层完全覆盖情况下的边界层空气动力学和冠层内空气动力学阻力(s/m)；LAI为叶面积指数。

1.4.6 覆膜条件下的阻力项r^p_m和r^m_s 覆膜下冠层内空气动力学阻力的求解同不覆膜的情况相同，可由式(31)、(32)计算。

　　土壤与膜之间的显热交换是借助空气介质，为自由对流形式。

r^m_s=ρc_pd/k_hNu	(33)
Nu=c(G_r·P_r)^m	(34)
G_r=αgd³(T_s-T_m)/v²	(35)
P_r=0.71,α=1/273,v=0.15,c=0.15,m=0.25	(36)

式中：k_h为热传导率;d=0.125m，为特征长度;Nu为Nusselt数;Gr为Grashof数;Pr为Prandtl数;v为液体粘滞系数(cm²/s).

1.5 SPAC系统中能量参数和水热参数的确定 模型中的能量参数是一个动态变化过程。本文根据在北京水科所永久店试验站的实际观测结果拟合出下列计算公式。

r_pg=-8×10^-5Id²+0.0079Id+0.0828	(37)
k=-3×10^-6Id³+0.0002Id²+0.0148Id	(38)
r_mg=-0.0003Id²+0.0175Id+0.0398	(39)

式中：r_pg，r_mg分别为冠层和地膜的短波反射率；k为消光系数；Id为日序数，从3月15日开始计算。

　　其他光学参数的取值分别为：t_pg=1-k-r_pg;t_p1=t_pg;t_mg=0.85；r_sg=0.20；r_m1=0.12；t_m1=0.78；ε_p=k；ε_m=0.12;ε_s=0.80或ε_s=0.80+0.18θ；a_p1=1-r_p1-t_pl;a_pg=k;a_mg=1-t_mg-r_mg；a_m1=1-r_m1-t_m1。

　　土壤水分特征曲线采用Van Genuchten M. Th^[11]提出的公式，式中的参数根据实测数据拟合确定：

θ-θ_r/θ_s-θ_r=[1+(αh)ⁿ]^-m	(40)
m=1-1/n	(41)

式中：θ_s、θ_r分别为土壤饱和含水率和残余含水率(m³/m³);h为土壤基质势的绝对值(或土壤水吸力s);α的单位为(1/cm);n、m为无量纲。

表1 土壤各种成分的性质分数成分


成分	I(下标)	C_I/(J/cm³·C)

液体水	1(l)	4.18
空气	2(a)	1.3×10^-3
石英	3(s)	1.92
其它矿物值	4(s)	1.92
有机质	5(s)	2.51

上表中的下标l-表示液体；a表示空气；s表示固体。

　　各层土壤非饱和水力传导度由式(42)进行拟合：

K_w(h)=p｜h｜^Q

(42)

　　土壤热传导系数采用Chung^[12]给出的经验公式：

λ=(b₁+b₂θ_L+b₃θ²_L)/100

(43)

　　土壤体积热容量C主要与土质有关，一般而言，在土壤组成中，固体部分变化不大，土壤热容量的大小主要决定于水和空气所占的比例，可采取下式计算：

(44)

各参数的取值见表1.

　　水蒸汽传导率D_hv和热蒸汽扩散率D_Tv可分别表示为：

D_hv=ρ_L^-1D_aΩθ_vρ_vg/(R·(T+273.16))	(45)
D_Tv=ρ^-1_LfD_aθ_vρ_v(R_l-hg/R)/(T+273.16)²	(46)

式中：g为重力加速度(9.81m/s²)；R为水蒸汽汽体常数(416.5m²/s²·k^-1)；D_a为水汽分子扩散度(m²/s)；Ω为充汽孔隙区域的弯曲度；f为Philip de Vres对水汽扩散引入的修正因子。

_{根系吸水项可由下式计算^[9]：

s_r(z,t)=A/v(t)Tp·e^[-B(-C)2]
(47)

式中：T_p为作物蒸腾速率(m/s)；A，B，C为经验常数，根据实测资料分析，可求出，A=1.3412，B=2.6818，C=0.2703.}

2 模型求解

　　前面已对求解模型所涉及到的条件及参数作了说明。模型的求解即是对能量平衡方程(1)～(4)和土壤水热运移方程(17)、(18)的求解。

　　能量平衡方程(1)～(4)采用Newton Raphson非线性方程组迭代法求解。求解时，时间步长取200s，空间步长取2cm，方程的控制精度取为：ε_T=0.01℃，ε_h=0.01cmH₂O.由此可解出T_p,T_l,T_m,T_s四个变量。求解出的时段末地表温度T_s即为求解土壤水热方程(17)、(18)时温度的第一类上边界条件。

　　土壤水热运移方程(17)(18)采用有限差分法进行数值求解。

3 模型验证

　　现采用本文提出的模型，对北京通县永乐店地区1998～1999年覆膜下的冬小麦生育期内的田间土壤水热运移问题进行数值模拟。模拟从1999年3月15日开始，到6月10日结束，分别对覆膜条件下的土壤体积含水率剖面、土壤温度剖面和SPAC系统中上边界温度进行模拟，选取有实测资料可比较的时段输出结果，分别见图2，图3，图4，图5a和图5b.

15t02.gif (6195 bytes)

图2 冬小麦返青～拔节期土壤体积含水率剖面分布的计算与实测值比较


图3 冬小麦返青～拔节期土壤温度剖面分布的计算与实测值比较	图4 冬小麦拔节～抽穗期土壤温度剖面分布的计算与实测值比较

图5(a) SPAC系统中Tm实测值与模拟值的比较	图5(b) SPAC系统中Ts实测值与模拟值的比较

　　从图2、图3和图4可看出土壤体积含水率剖面和土壤温度剖面的计算值与相应的实测值吻合都比较好。由图5(a)和图5(b)可以看出，上边界温度的计算尽管个别时段与实测值的偏差也比较大，但总体上来看，两者差异不大。由此可见本文所建立的模型基本上反映了地膜覆盖条件下SPAC系统水热运移这一复杂问题的物理本质，从而为研究地膜覆盖条件下作物的耗水规律及地膜覆盖农田的水管理提供了可靠的分析工具。应该指出，地膜覆盖条件下SPAC系统的水热运移规律是一个十分复杂的问题。本文所建立的模型作了许多必要的简化。模型中涉及大量的参数，有些参数本身也有待进一步研究和完善。

4 结语

　　本文将地膜覆盖条件下的SPAC系统概化为由大气、植物冠层、地膜覆盖层和土壤四层结构组成的系统，并根据能量平衡原理和土壤水热运移理论建立了地膜覆盖下SPAC系统水热运移模型。采用北京市水科所永乐店试验站1998 1999年冬小麦覆盖农田的试验资料对模型进行了验证。结果表明，该模型具有较好的仿真性，所模拟的土壤水分剖面和温度剖面以及上边界温度与实测数据吻合较好，从而为分析研究地膜覆盖条件下作物的耗水规律提供了理论依据。

参考文献：

[1] 黄介生，沈荣开。地膜覆盖技术的研究现状和展望，二十一世纪农田水利学术研讨会论文集[C]。1997，80-81.

[2] Huang Jiesheng, Shen Rongkai。Analysis of crop evapotranspiration under transparent polyethylene mulch, 17^th ICID Congress on Irrigation and Drainage[C]。1999,Vol.1A,327-338.

[3] 黄洪峰。土壤植物大气相互作用原理及模拟研究[M]。北京：气象出版社，1997.6.

[4] 张蔚榛。地下水和土壤水动力学[M]。北京：中国水利电力出版社，1999.9.

[5] 王康。塑膜覆盖条件下作物耗水规律的研究[D]。武汉：武汉水利电力大学硕士论文，1999，6.

[6] Alves I, Perrier A, Pereira L S. Aerodynamic and surface resistance complete cover crops:How good is the “Big Leaf”?. Transaction of ASAE, Vol.41(2):345-351.

[7] 孙菽芬。土壤内水分流动及温度分布计算耦合模型[J]。力学学报，1987，4.

[8] Peter J, Camilla, Robert J, Gurney.A resistance parameter for bare soil evaporation models[J]。Soil Science, 141(2).

[9] 李会昌。SPAC中水分运移与作物生长动态模拟及其在灌溉预报中的应用研究[D]。武汉：武汉水利电力大学博士论文，1997.9.

[10] James Shuttleworth W, Wallace J S.Evaporation from spare crops anenergy combination theory[M]。J. R. Met. Soc., 1985.839-855.

[11] Van Ghenuchten M. Th.A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J]。Soil Sci. Soc. Am. J., 1982,44:892-898.

[12] Chung S O, Horton R. Soil heat and water flow with a partial surface mulch[J]。Water Resour. Res., 1987,23:2175-2186.