摘要:本文描述了夏玉米生长期间麦秸条带覆盖(行间覆盖)条件下土壤水热动态的田间试验.在此基础上,采用交替方向有限差分法建立了模拟苗期条带覆盖下土壤二维水热迁移的数值模型.其中覆盖部分土壤表面的边界条件由能量平衡方程确定.模型的输入只需常规气象资料、土壤水热迁移参数及覆盖参数.数值计算结果表明,该模型具有一定的实用价值.
关键词:条带覆盖 水热动态 模型 模拟
蒸发条件下土壤表层的水分运动受温度和水势梯度的控制,但传统的描述非饱和土壤中水分运移规律的模型均假定土壤为等温系统,如:Gardner(1959)[1]、Hanks等(1969)[2]、Staple(1971)[3]、任理(1991)[4]. Philip(1957)从理论上得出,除非在非常低的含水率条件下,否则由温度引起的水分流动相对来说并不重要[5,6].然而,Jackson等(1974)的计算表明,在田间条件下,由温度引起的水汽流动在中等土壤含水量时就能显著地表现出来[5,7].裸土或作物苗期土壤水分的损失,大部分是穿过表层10cm左右的干土层发生的.日温度大的变化发生在上层15cm的土壤中,它可以影响穿过干燥层的水汽流动(Papendick,1973)[8],因此,模拟和分析蒸发条件下的土壤水分动态应同时考虑水和热的传输.
Philip与de Vries(1957,1958)提出了描述土壤水热耦合迁移的理论[9,10],近二十年来,国内外学者对蒸发条件下土壤水热迁移的耦合计算进行了广泛的研究[11,12,13,14,15,16,17,18]. 在二维土壤水热迁移问题的研究方面,Jury和Bellantuoni(1976)发展了一个反映表面铺盖矩形岩块的均匀田间土壤在温度梯度下热流和水汽运动的二维数学模型,结果发现,只有考虑包括温度与热传导关系时,计算值才与实测值有很好的一致性[19,20]. Chung和Horton(1987)对地面采用部分覆盖下的土壤水热流动进行了数值试验,但没有田间实测资料的检验[21].杨邦杰(1989)对土壤不均匀、地表平坦或起伏不平时的二维土壤蒸发过程的数值模型进行了研究[22]. Sui Hongjian和Zeng Dechao等(1992)用数值模型对不同覆盖下土壤温度和水分动态进行了模拟[23].
为了探讨行间条带覆盖对夏玉米生长条件下的土壤水热动态的影响,作者在北京通县永乐店试验站进行了田间试验,并本着简捷实用的原则,依据Philip和de Vries(1957,1958)提出的土壤水热流动理论和已有的研究成果,以夏玉米生长前期麦秸条带覆盖下的田间试验为背景,建立了土壤二维水热迁移的数值模型.
2 田间试验
2 .1 试验布置 田间玉米行间裸地的麦秸覆盖宽度约30cm(玉米行距为60cm).覆盖量相当于400kg/亩.在试验小区内,沿覆盖层中线、边缘及无覆盖的裸地设3个土壤温度剖面,这3个剖面水平相距分别为15cm和10cm.剖面上测点埋深为5cm、10cm、20cm、30cm、50cm.在覆盖层与土壤交界面处用曲管地温计量测界面处的地表温度,在对照区地表和覆盖层表面采用直管温度计测定温度.用于测量土温的铂热电阻安装前均进行了率定.观测时使用万用表测定铂热电阻值,然后依据分度表及田间校正值拟合的标准曲线换算出相应的土壤温度.中子管埋设在麦秸覆盖层中线,水分动态由标定后的中子仪测量.
2.2 试验结果分析 图1反映了麦秸覆盖层中线下土壤温度随时间的变化过程.图2、
图1 覆盖层中线下土壤剖面实测温度(1993.7.3-7.4)
图2 覆盖层边缘下土壤剖面实测温度(1993.7.3-7.4)
图3 距覆盖层边缘10cm处裸地土壤剖面实测温度(1993.7.3-7.4)
图3分别为覆盖层边缘下及距离覆盖层边缘10cm处裸地土壤剖面的温度动态.此时夏玉米为苗期,其遮荫作用很微弱,这样只有覆盖层对太阳辐射具有“屏蔽”作用.由图3可见,在距覆盖层边缘10cm处的玉米幼苗附近,裸地温度随时间的变化幅度明显大于覆盖层中线以下地温的变化幅度(图1).因为裸地土壤较干燥,表面温度可达到42℃以上,而在覆盖层内的土壤表面,最高温度约为32℃左右.从图2可见,覆盖层边缘下土壤表层的温度变化幅度明显小于裸地(图3)而大于覆盖层中线下的温度变幅(图1).此外,地温动态的观测表明,随着深度增加,土壤温度变幅减小,增加了相位滞后,这是土壤一个周期温度波的典型传播.
图4为条带覆盖、全覆盖与无覆盖土壤表面的温度变化过程图,图示表明,条带覆盖条件下土表温度介于全覆盖和无覆盖之间,它与无覆盖相比,可起到降低表土水分蒸发的作用,但同时又较全覆盖情况下的表土温度高,有利于玉米出苗、生长.
图5为条带覆盖、全覆盖与无覆盖条件下玉米最终产量比较图,图示明显可见,条带覆盖的玉米产量最高,说明虽然与全覆盖的覆盖量(400kg/亩)相同,条带覆盖对节水、保墒,促进农业增产更加有效。麦秸覆盖对节墒是有效措施,这一点早已被证实,但由于麦秸覆盖会降低土壤温度,对夏玉米前期生长是不利的。条带覆盖仅铺设在作物行间,一方面可以减少行间土面的无效蒸发;另一
图4 不同处理土壤表面温度
图5 不同覆盖处理产量
方面,植株部分可以充分接受太阳辐射.在夏玉米生长后期,由于覆盖层的压实,对土壤通气和热状况均有不良影响,而条带覆盖却可免除,也许这就是条带覆盖产量较高的原因.所以,对于条播作物,这种覆盖形式显然是值得推荐的.
3.1 控制方程及数值格式 夏玉米生长前期作物的根系吸水可以忽略,因此所研究的系统仅考虑土壤、覆盖和大气因素,由于田间麦秸覆盖条带是平行和空间上等距的,基于对称性,只分析流动区域的一半即可[21].
Philip和de Vries(1957)提出了非稳定耦合的土壤水热流动方程如下[21]:
这里C是土壤体积热容量(J/m·℃),T是土壤温度(℃),t是时间(s),λ是热传导度(W/m·℃),L是汽化的体积潜热(J/m),θ是体积含水量(m/m),Dθv是等温水汽扩散度(m2/s),K是水力传导度(m/s),h为负压(m),Z为垂直距离,向下为正(m),为梯度算子.
本文只在土壤表面考虑水汽对热和水分传输的影响,不考虑地下水汽流动[21],这样方程(1)可写成:
方程(2)又可写为:
Milly(1984)指出,在大多数土壤含水量情况下,土壤热液体流动并不重要[13],故(4)式可简化为:
采用交替方向隐式(ADI)有限差分法离散方程(3)和(5),则将二维问题降为一维问题来处理,ADI方程如下:
X方向隐式,Z方向显式:
Z方向隐式,X方向显式:
(9)
式中上标代表时间,下标代表空间,i为行标记,j为列标记,F为容水度(m-1).
因为方程(6)到(9)中的系数依赖于变量本身,所以方程为非线性的.本文采用显式线性化,即以前一时间步的值来近似方程(6)到(9)中的系数.经整理,方程(6)至(9)可写成:
式中Ebs和Ems分别为裸土和有覆盖的土壤表面的蒸发通量(m/s),Ho为土壤表面空气的绝对湿度(kg/m),Ha为土壤表面之上空气的绝对湿度(kg/m),ra是土壤表面和其上空气之间的空气动力学边界层阻力(s/m).rm是覆盖层的水分扩散阻力(s/m).
Ho和ra的计算公式为[21]
Ho=H*oexp〔h1/46.97(Ts+273.16)〕, (16)
ra=〔ln(2.0/Zo)〕2/(0.16Ws) , (17)
这里H*o是在土壤表面温度时的饱和温度(kg/m),h1是土壤表面的负压(m),Zo是粗糙度长度(m),Ws是风速(m/s).
空气的绝对湿度Ha和在土壤表面温度时的饱和湿度H*o由下式计算[21]:
Ha=1.323exp〔17.27Td/(Td+237.3)〕/(Ta+273.16),(18)
H*o=1.323exp〔17.27Ts/(Ts+237.3)〕/(Ts+273.16), (19)
式中Td,Ta,Ts分别是露点温度(℃)、空气温度(℃)、地表温度(℃).
为简化计算,本文把能量平衡方程仅用于覆盖层和土壤层的界面上.在此我们假设条带麦秸覆盖层为不透明覆盖层,这样辐射便不能穿透到覆盖表面之下.于是,对于覆盖层与土壤的交界面,能量平衡方程为[21]:
Ms-LEms-G=0 , (20)
这里Ms为覆盖热通量(w/m2),向下为正,LEms为潜热通量(向上为正),L、Ems意义同前,G为土壤热通量(向下为正).Ms、L和G的表达式如下[21]
Ms=λm(Tm-Ts)/THK, (21)
L=2.4946×109-2.247×106+6Ts , (22)
G=λ(Ts-T2)/(ΔZ)+ρsCps(Ts-T0s)(ΔZ)/(2Δt),(23)
式中λm是覆盖层的热传导度(W/m℃),Tm是覆盖层表面的温度(℃),THK是覆盖层厚度(m),后两个参数均由田间实测.T2是前一时间步在土壤表面以下ΔZ处结点的温度(℃),T0s是前一时间步的Ts(℃),ρs为土壤密度(kg/m),其它符号意义同前.Cps是常压下土壤的比热(J/kg ℃),其计算公式为[24]:
Cps=1000(0.2+θo/1.36)/〔0.238846(1+θo/1.36)〕, (24)
式中θo是地表含水量(m/m).
裸土表面的温度,根据气象观测数据由如下正弦函数确定:
Ts=s+Assin(2πt/86400+1.5π), (25)
这里s为模拟期间裸土表面温度的平均值(℃),As为Ts的变幅,分别为28.2℃和11℃.
条带覆盖与土壤交界面的温度采用如下步骤确定,首先由实测的麦秸覆盖层表面温度和覆盖层厚度确定覆盖层的热通量,然后将式(22)、式(23)、式(21)和式(15)代入式(20),使用二分法得到覆盖层与土壤交界面的温度Ts.
在求得裸土表面温度及覆盖与土壤交界面的温度后,由式(14)、(15)可分别得到裸土部分和覆盖部分土壤表面的蒸发通量.
其它特征量包括:XL(计算域宽度),ZL(计算域深度),Δx、ΔZ和Δt(空间和时间步长),Zo(粗糙度长度),TL(模拟总时间).
空气温度和露点温度变化用如下正弦函数来确定[16]:
这里a和d分别为日平均气温和日平均露点温度(℃),Aa和Ad分别代表各自的变化幅度(℃),t是从午夜开始一天的时间(s).
土壤热力传导度由以下经验方程计算:[21]:
这里λ是热传导度(W/m℃),θ是体积含水量(m/m),b1\,b2,b3为回归参数.
根据de Vries(1963)[25]、吴擎龙(1993)[26]土壤体积热容量的计算公式可简化为:
式中θs为土壤饱和含水量(m/m).
土壤水分特征曲线、水力传导度和容水度由 van Genuchten(1980)提出的经验方程来描述[27]:(以下依次为(30),(31),(32)):
(30)(31)(32)
这里θs和θr是饱和及残余含水量(m/m), Ks是参考温度时的饱和水力传导度(m/s), α、n、m是描述土壤水分特征曲线形状的非线性回归参数.考虑到温度,水力传导度应校正为[21]:
式中μ为粘度,T0为参考温度.
覆盖层的热传导度、水分扩散阻力及粗糙度长度的数值选自有关文献.
4 模型的验证
对于整个二维水热迁移模型,不存在解析解.本文首先只对ADI数值模型中的热流方程进行验证[21],其次运行整个模型与田间实测资料进行对比.
考虑到田间热传输问题的边界条件为Dirichlet条件和Neumann条件,所以取两个热传导算例检验之.算例1[28]的问题是方形板(边长2l为5)的热流传输,其初始条件为Ti=1,边界条件为Tb=0. Kt/l2=0.08,这里K是物质的温度计传导度,取K=0.005,求t=100时板的温度分布.算例2[29]为一个长钢棒的热传导问题,由于传导热流是对称的,所以只分析钢棒横截面的1/4区域(0.5m×0.25m),数值模拟使用的时间步长Δt=5sec,空间步长Δx=0.01m、ΔZ=0.01m.此钢棒的热力参数为:λ=20W/m·℃,ρ=3000kg/m, C=1000J/kg·℃.边界条件包括绝热边界(Neumann条件)和对流热传输边界(Cauchy条件).对流热传输系数h=10W/m·℃.钢棒的初始温度是300℃.环绕钢棒的空气流温度保持在20℃.模拟t=3600sec时的温度分布.下面给出两算例的解析解与数值解(图6、图7),可见两者吻合很好.
根据试验资料,确定数值模拟的定解条件和参数.具体地,以麦秸覆盖第二天上午8时的土壤水分剖面(假设x方向均匀,Z方向变化)为数模的土壤水分初始条件.
图6 方形板的温度分布
图7 钢棒中的热流分布
田间土壤的水热参数见表1:
参数* | Ks/(m/s) | θs/(m/m) | θr/(m/m) | a/(m) | n | m | b1 | b2 | b3 |
粉砂土 | 0.00001 | 0.48 | 0.12 | 0.6892283 | 2.170972 | 0.5393769 | 0.243 | 0.393 | 1.534 |
* Ks、θs、θr值均为田间实测,a、n、m是van Genuchten方程的参数,拟合得到,b1、b2、b3是热传导度公式中的系数,引自文献[21].
模型中输入的有关参数和数据分别列于表2和表3.
表2 模型输入参数
符 号 | 数 值 | 备 注 | |
DX | x坐标空间步长 | 0.05m | |
DZ | z坐标空间步长 | 0.05m | |
XL | x坐标长度 | 0.25m | |
ZL | z坐标长度 | 0.90m | |
DT | 时间步长 | 1.0s | |
TL | 模拟时间 | 172800s | |
To | 参考温度 | 20℃ | 引自[21] |
ρs | 土壤密度 | 1360kg/m | 实测 |
ρa | 空气密度 | 1292.8kg/m | 引自[30] |
Cpa | 空气的定压比热 | 1006.09J/kg℃ | 引自[30] |
ML | 覆盖层宽度 | 0.30m | 实测 |
THK | 覆盖层厚度 | 0.10m | 实测 |
λm | 覆盖层的热传导度 | 0.126W/m·℃ | 引自[21] |
rm | 覆盖层的水分扩散阻力 | 4800s/m | 据[21],假定 |
Zo | 土壤表面的粗糙度长度 | 0.01m | 引自[21 |
日期 | a/(℃) | Aa/(℃) | d/(℃) | Ad/(℃) | Ws/(m/s) | Tm/(℃) |
6.26 | 26.25 | 8.25 | 15.05 | 2.55 | 1.3 | 43.5 |
6.27 | 27.25 | 7.75 | 13.45 | 2.35 | 1.5 | 41.0 |
模拟时段内(6月25日至6月27日)的表土含水量用取土称重法加以校正.
模拟结果如下图所示.由图8可见,模拟的表层埋深10cm处的土壤水分横向分布值与实测值趋势有较好的一致性.图9所示为表层不同深度土壤温度的分布,计算与实测值吻合良好.图10为无覆盖处(距条带覆盖中线25cm)土壤表层温度分布,结果很好.由此可见,条带覆盖部分土壤含水率高于无覆盖区,地温则低于未覆盖部分,地表温度变幅较大,越向下层温度变幅越小,说明对条带覆盖只有用二维模型才能较真实地刻划土壤水热运动规律,特别是表层土壤的水热动态.
图8 表层土壤水分分布
1 Gardner W R. Solution of the flow equation for the drying of soils and other porous media. Soil Sci. Soc. Am. Proc., 1959,23,183-187.
2 Hanks R J, Klute, A and Bresler E. A numeric method for estimating infiltration, redistribution drainage, and evaporation of water from soil. Water Resour. Res., 1969.5,1065~1069
3 Staple W J. Boundary conditions and conductivities used in the isothermal model of evaporation from soil. Soil Sci. Soc. Am. Proc., 1971.35,853~855.
4 任 理. 野外非饱和土壤水动态的数值模拟. 武汉电力学院学报,1991,24(3):354-360.
5 Hammel J E, Papendick,R I and Campbell. G S Fallow tillage effects on evaporation and seedzone water content in a dry summer climate. Soil Sci.Soc.Am.J., 1981,45,1016-1022
6 Philip J R. Evaporation, and moisture and heat fields in the soil. J.Meteorol., 1957,14(4):354~366.
7 Jackson R D, Reginato, R J Kimball, B A and Nakayama. F S Diurnal soil-water evaporation comparison of measured and calculated soi-|water fluxes. Soil Sci.Soc.Am. Proc., 1974,38,863~866.
8 Papendick, R I, Lindstrom,M J and Cochran.V L Soil mulch effects on seedbed temperature and water during fallow in eastern Washington. Soil Soc. Am. Proc.,1973,37,307~314.
9 Philip J R, and de Vries.D A Moisture movement in porous materials under temperature gradients. Eos Trans. AGU, 1957,38(2):222~232.
10 de Vries D A. Simultaneous transfer of heat and moisture in porous media. Eos Trans. AGU, 1958,39(5):909~916.
11 Sophocleous M. Analysis of water and heat flow in unsaturated-saturated porous media. Water Rescur. Res., 1979,15(5):1195~1206.
12 Milly P C D. Moisture and heat transport in hysteretic., inhomogeneous porous media, Water Resour. Res., 1982,18(3):489~498.
13 Milly P C D. A simulation analysis of thermal effects on evaporation from soil. Water Resour,Res., 1984,20(8):1087~1098.
14 Passerat de Silans et al., Numerical modeling of coupled heat and water flows during drying in a stratified bare soil-|Comparison with field observations. J.Hydrol., 1989,105,109~138.
15 孙菽芬. 土壤内水分流动及温度分布计算——耦合型模型.力学学报,1987,19(4):374-380.
16 杨金忠,蔡树英. 土壤中水、汽、热运动的耦合模型和蒸发模拟.武汉电力学院学报,1989,22(4):35-44.
17 张瑜芳,蔡树英等.表土低含水率条件下土壤非稳定蒸发研究.武汉电力学院学报,1991,24(2):157-164.
18 蔡树英,张瑜芳,温度影响下土壤水分蒸发的数值分析.学报,1991,(11),1-8.
19 Jury W A, and Bellantuoni.B Heat and water movement under surface rocks in a field soil: Ⅰ.Thermal effects. Soil Sci.Soc.Am.J., 1976,40,505~509.
20 Jury W A, and Bellantuoni. B Heat and water movement under surface rocks in a field soil:Ⅱ.moisture effects. Soil Sci.Soc.Am.J., 1976,40,509~513.
21 Chang S O, and Horton. R Soil heat and water flow with a partial surface mulch. Water Resour. Res. 1987,23(12):2175~2186.
22 杨邦杰. 土壤蒸发过程的数值模型及其应用.北京:学术书刊出版社,1989年.
23 Sui H, Zeng, D and Chen F. A numerical model for simulating the temperature and moisture regimes of soil under various mulches. Agric. For. Meteorol.,1992,61,281~289.
24 汉克斯 R J,阿希克洛夫编著,G L,杨诗秀,陆锦文,等译,华孟校,应用土壤物理.北京:电力出版社,1984.
25 de Vries D A. Thermal properties of soils, in Physics of Plant Environment, edited by R.W. van Wijk, pp. 210-235, North-Holland, Amsterdam, 1963.
26 吴擎龙. 田间腾发条件下水热迁移数值模拟的研究.[学位],清华大学,1993.
27 van Genuchten, M Th. A closed\|form equation for predicting the hydraulic conductivity of unsaturated soil, Soil Sci.Soc.Am.J., 1980,44,892~898.
28 Carslaw H S and Jaeger. J C Conduction of Heat in Solids. Oxford University Press, pp.174,Second. edition, 1959.
29 Benjamin J G, Ghaffarzadeh M R and Cruse. R M Coupled water and heat transport in ridged soils. Soil Sci.Soc.Am.J., 1990,54,963~969.
30 朱炳海,王鹏飞,等主编. 气象学词典.上海:上海辞书出版社,1985.
31 任 理. 地下水溶质运移计算方法及土壤水热动态数值模拟的研究.[学位].武汉电力大学,1994.
Field experiments and numerical simulation of soil water and heat
regimesunder the condition of summer corn partially covered by mulch strips
Abstract The field experiments of soil moisture and temperature regimes under strip wheat straw mulch are introduced in this paper. A numerical model using the alternating direction implicit(ADI) finite difference method to study two-dimensional soil heat and water flow with a partial surface mulch cover is established. In the model, a soil surface energy balance equation is used to determine soil surface boundary condition of mulch data. The inputs required for the numerical simulations are common weather data, soil thermal and hydraulic parameters as well as mulch data. The heat flow model is verified with analytical solutions. The proposed model successfully simulates the field moisture content and temperature regimes under the condition of strip wheat straw mulch at summer seeding stage.
Key words strip mulch, soil moisture and temperature regimes, model, simulation