土壤水熱動態分析論文

導語：土壤水熱動態分析論文一文來源于網友上傳，不代表本站觀點，若需要原創文章可咨詢客服老師，歡迎參考。

蒸發條件下土壤表層的水分運動受溫度和水勢梯度的控制，但傳統的描述非飽和土壤中水分運移規律的模型均假定土壤為等溫系統，如：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.