土壤水热动态分析论文

蒸发条件下土壤表层的水分运动受温度和水势梯度的控制，但传统的描述非饱和土壤中水分运移规律的模型均假定土壤为等温系统，如：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与deVries(1957,1958)提出了描述土壤水热耦合迁移的理论［9，10］，近二十年来，国内外学者对蒸发条件下土壤水热迁移的耦合计算进行了广泛的研究［11，12，13，14，15，16，17，18］.在二维土壤水热迁移问题的研究方面，Jury和Bellantuoni(1976)发展了一个反映表面铺盖矩形岩块的均匀田间土壤在温度梯度下热流和水汽运动的二维数学模型，结果发现，只有考虑包括温度与热传导关系时，计算值才与实测值有很好的一致性［19，20］.Chung和Horton(1987)对地面采用部分覆盖下的土壤水热流动进行了数值试验，但没有田间实测资料的检验［21］.杨邦杰(1989)对土壤不均匀、地表平坦或起伏不平时的二维土壤蒸发过程的数值模型进行了研究［22］.SuiHongjian和ZengDechao等(1992)用数值模型对不同覆盖下土壤温度和水分动态进行了模拟［23］.

为了探讨行间条带覆盖对夏玉米生长条件下的土壤水热动态的影响，作者在北京通县永乐店试验站进行了田间试验，并本着简捷实用的原则，依据Philip和deVries(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数值模型的建立

3.1控制方程及数值格式夏玉米生长前期作物的根系吸水可以忽略，因此所研究的系统仅考虑土壤、覆盖和大气因素，由于田间麦秸覆盖条带是平行和空间上等距的，基于对称性，只分析流动区域的一半即可［21］.

Philip和deVries(1957)提出了非稳定耦合的土壤水热流动方程如下［21］：

C(T)/(t)=·(λT)-L·(Dθvθ)，(1)

(θ)/(t)=·(Kh)-(K)/(Z)，(2)

这里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)可写成：

C(T)/(t)=·(λT)，(3)

方程(2)又可写为：

(4)

Milly(1984)指出，在大多数土壤含水量情况下，土壤热液体流动并不重要［13］，故(4)式可简化为：

F(h)/(t)=·(Kh)-(K)/(Z).(5)

采用交替方向隐式(ADI)有限差分法离散方程(3)和(5)，则将二维问题降为一维问题来处理，ADI方程如下：

X方向隐式，Z方向显式：

(6)(7)

Z方向隐式，X方向显式：

(8)

(9)

式中上标代表时间，下标代表空间，i为行标记，j为列标记，F为容水度(m-1).

因为方程(6)到(9)中的系数依赖于变量本身，所以方程为非线性的.本文采用显式线性化，即以前一时间步的值来近似方程(6)到(9)中的系数.经整理，方程(6)至(9)可写成：

式中：

.式(10)至(13)均为三对角方程，结合边界条件，用追赶法求解.内部结点的系数由相邻结点的算术平均值来确定.

3.2上边界条件的确定在有限差分法中有效地处理通量边界条件是最困难的部分［21］.在本文中，热流问题的顶部和底部边界为Dirichlet条件.热流和水流的左、右边界使用Neumann条件，亦即没有流动的边界条件.对于水流问题，其顶部边界使用非零通量的Neumann条件，底部为Dirichlet条件.

在未覆盖的裸土表面和覆盖层与土壤层的界面上，水流问题的Neumann条件由以下公式确定［21］：

Ebs=(Ho-Ha)/(1000ra),(14)

Ebs=(Ho-Ha)/〔(1000(ra+rm)〕,(15)

式中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)可分别得到裸土部分和覆盖部分土壤表面的蒸发通量.

3.3参数的选取本文数值模型的运行只需一般的气象观测资料及覆盖和土壤参数.气象资料是日最高和最低气温、日最大和最小露点温度、日最高和最低地表温度及日平均风速.覆盖参数为覆盖宽度、厚度，覆盖层的热传导度、水分扩散阻力，覆盖表面的温度.土壤参数为土壤热力传导度、土壤体积热容量、土壤水力传导度和容水度及土壤温度和含水量的初始分布，土壤剖面下边界处的温度和含水量.

其它特征量包括：XL(计算域宽度)，ZL(计算域深度)，Δx、ΔZ和Δt(空间和时间步长)，Zo(粗糙度长度)，TL(模拟总时间).

空气温度和露点温度变化用如下正弦函数来确定［16］：

Ta=a+Aasin(2πt/86400+1.5π),(26)

Td=d+Aasin(2πt/86400+1.5π),(27)

这里a和d分别为日平均气温和日平均露点温度(℃)，Aa和Ad分别代表各自的变化幅度(℃)，t是从午夜开始一天的时间(s).

土壤热力传导度由以下经验方程计算：［21］：

λ(θ)=b1+b2θ+b3θ0.5(28)

这里λ是热传导度(W/m℃),θ是体积含水量(m/m),b1/,b2,b3为回归参数.

根据deVries(1963)［25］、吴擎龙(1993)［26］土壤体积热容量的计算公式可简化为：

C=1.925×106(1-θs)+4.184×106θ,(29)

式中θs为土壤饱和含水量(m/m).

土壤水分特征曲线、水力传导度和容水度由vanGenuchten(1980)提出的经验方程来描述［27］：（以下依次为（30），（31），（32））：

(30)（31)(32)

这里θs和θr是饱和及残余含水量(m/m),Ks是参考温度时的饱和水力传导度(m/s),α、n、m是描述土壤水分特征曲线形状的非线性回归参数.考虑到温度，水力传导度应校正为［21］：

K(h,T)=K(h)(μ(T°))/(μ(T))=K(h)(1+0.0384T+0.000211T2)/(1+0.0384T0+0.000211T20),(33)

式中μ为粘度，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：

表1田间土壤的水热参数

参数*Ks/(m/s)θs/(m/m)θr/(m/m)a/(m)nmb1b2b3

粉砂土0.000010.480.120.68922832.1709720.53937690.2430.3931.534

*Ks、θs、θr值均为田间实测，a、n、m是vanGenuchten方程的参数，拟合得到，b1、b2、b3是热传导度公式中的系数，引自文献［21］.

模型中输入的有关参数和数据分别列于表2和表3.

表2模型输入参数

符号参数定义

数值备注

DXx坐标空间步长0.05m

DZz坐标空间步长0.05m

XLx坐标长度0.25m

ZLz坐标长度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

表3模型输入的数据

日期a/(℃)Aa/(℃)d/(℃)Ad/(℃)Ws/(m/s)Tm/(℃)

6.2626.258.2515.052.551.343.5

6.2727.257.7513.452.351.541.0

模拟时段内(6月25日至6月27日)的表土含水量用取土称重法加以校正.

模拟结果如下图所示.由图8可见，模拟的表层埋深10cm处的土壤水分横向分布值与实测值趋势有较好的一致性.图9所示为表层不同深度土壤温度的分布，计算与实测值吻合良好.图10为无覆盖处(距条带覆盖中线25cm)土壤表层温度分布，结果很好.由此可见，条带覆盖部分土壤含水率高于无覆盖区，地温则低于未覆盖部分，地表温度变幅较大，越向下层温度变幅越小，说明对条带覆盖只有用二维模型才能较真实地刻划土壤水热运动规律，特别是表层土壤的水热动态.

图8表层土壤水分分布

图10裸地(x=25cm)土壤温度剖面

图9表层土壤温度分布

5小结

在夏玉米生长前期的六月份，北方降雨量往往偏少，干旱威胁玉米壮苗.覆盖不仅阻碍了土壤水分的蒸发，且由于适当降低地温也减少了水分蒸发.本文所建立的土壤二维水热迁移模型，输入参数少，相对简单，却能较好地模拟出麦秸覆盖所产生的保墒效应，因而具有一定的实用价值.

致谢本文得到张蔚榛教授的指教，田间试验承北京水利科学研究所永乐店试验站同志们的协助.

参考文献

1GardnerWR.Solutionoftheflowequationforthedryingofsoilsandotherporousmedia.SoilSci.Soc.Am.Proc.,1959，23，183-187.

2HanksRJ,Klute,AandBreslerE.Anumericmethodforestimatinginfiltration,redistributiondrainage,andevaporationofwaterfromsoil.WaterResour.Res.,1969.5，1065～1069

3StapleWJ.Boundaryconditionsandconductivitiesusedintheisothermalmodelofevaporationfromsoil.SoilSci.Soc.Am.Proc.,1971.35，853～855.

4任理.野外非饱和土壤水动态的数值模拟.武汉水利电力学院学报，1991,24(3)：354-360.

5HammelJE,Papendick，RIandCampbell.GSFallowtillageeffectsonevaporationandseedzonewatercontentinadrysummerclimate.SoilSci.Soc.Am.J.,1981,45，1016-1022

6PhilipJR.Evaporation,andmoistureandheatfieldsinthesoil.J.Meteorol.,1957，14(4)：354～366.

7JacksonRD,Reginato,RJKimball,BAandNakayama.FSDiurnalsoil－waterevaporationcomparisonofmeasuredandcalculatedsoi-|waterfluxes.SoilSci.Soc.Am.Proc.,1974，38，863～866.

8Papendick,RI,Lindstrom,MJandCochran.VLSoilmulcheffectsonseedbedtemperatureandwaterduringfallowineasternWashington.SoilSoc.Am.Proc.,1973，37，307～314.

9PhilipJR,anddeVries.DAMoisturemovementinporousmaterialsundertemperaturegradients.EosTrans.AGU,1957，38(2)：222～232.

10deVriesDA.Simultaneoustransferofheatandmoistureinporousmedia.EosTrans.AGU,1958，39(5)：909～916.

11SophocleousM.Analysisofwaterandheatflowinunsaturated-saturatedporousmedia.WaterRescur.Res.,1979，15(5)：1195～1206.

12MillyPCD.Moistureandheattransportinhysteretic.,inhomogeneousporousmedia,WaterResour.Res.,1982，18(3)：489～498.

13MillyPCD.Asimulationanalysisofthermaleffectsonevaporationfromsoil.WaterResour,Res.,1984，20(8)：1087～1098.

14PasseratdeSilansetal.,Numericalmodelingofcoupledheatandwaterflowsduringdryinginastratifiedbaresoil-|Comparisonwithfieldobservations.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.

19JuryWA,andBellantuoni.BHeatandwatermovementundersurfacerocksinafieldsoil:Ⅰ.Thermaleffects.SoilSci.Soc.Am.J.,1976，40，505～509.

20JuryWA,andBellantuoni.BHeatandwatermovementundersurfacerocksinafieldsoil:Ⅱ.moistureeffects.SoilSci.Soc.Am.J.,1976，40，509～513.

21ChangSO,andHorton.RSoilheatandwaterflowwithapartialsurfacemulch.WaterResour.Res.1987，23(12)：2175～2186.

22杨邦杰.土壤蒸发过程的数值模型及其应用.北京：学术书刊出版社，1989年.

23SuiH,Zeng,DandChenF.Anumericalmodelforsimulatingthetemperatureandmoistureregimesofsoilundervariousmulches.Agric.For.Meteorol.,1992，61，281～289.

24汉克斯RJ，阿希克洛夫编著，GL，杨诗秀，陆锦文，等译，华孟校，应用土壤物理.北京：水利电力出版社，1984.

25deVriesDA.Thermalpropertiesofsoils,inPhysicsofPlantEnvironment,editedbyR.W.vanWijk,pp.210-235,North-Holland,Amsterdam,1963.

26吴擎龙.田间腾发条件下水热迁移数值模拟的研究.［学位论文］，清华大学，1993.

27vanGenuchten,MTh.Aclosed/|formequationforpredictingthehydraulicconductivityofunsaturatedsoil,SoilSci.Soc.Am.J.,1980，44，892～898.

28CarslawHSandJaeger.JCConductionofHeatinSolids.OxfordUniversityPress,pp.174,Second.edition,1959.

29BenjaminJG,GhaffarzadehMRandCruse.RMCoupledwaterandheattransportinridgedsoils.SoilSci.Soc.Am.J.,1990，54，963～969.

30朱炳海,王鹏飞，等主编.气象学词典.上海：上海辞书出版社，1985.

31任理.地下水溶质运移计算方法及土壤水热动态数值模拟的研究.［学位论文］.武汉水利电力大学，1994.