蒸發(fā)條件下非飽和土水分遷移的數(shù)值模擬

熊 義 高 抗

（1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌 443002；2.三峽大學(xué)三峽庫區(qū)地質(zhì)災(zāi)害教育部重點實驗室，湖北宜昌 443002）

大氣和土表的水量交換是水循環(huán)當(dāng)中的重要組成部分.滲透作用使大氣中的水分進入土體，蒸發(fā)蒸騰作用又使水分離開土體進入大氣，由此而引起的土體含水率的變化直接影響著土體的穩(wěn)定性.例如，膨脹土地基上土體遇水膨脹、失水收縮，可能引起上部建筑物的隆起、傾斜或開裂［1］.黃土的濕陷特性導(dǎo)致其在含水率增高時強度降低，在有荷載作用的情況下出現(xiàn)地基不均勻沉降的危害［2］.因此，研究自然條件下非飽和土的水分遷移規(guī)律對把握土體變化及工程應(yīng)用有十分重要的意義.

長期以來，人們通常采用經(jīng)典的Richards等溫流方程描述土壤中液態(tài)水的運動［3］.例如，Thomas等［4］采用該方程的一維形式模擬了在氣候變化下土體脫濕變干行為.吳宏偉等［5］采用該方程的二維形式分析了雨水入滲對非飽和土坡穩(wěn)定性影響.然而，土壤蒸發(fā)是土壤中的水由液態(tài)水轉(zhuǎn)化為氣態(tài)水的過程，是一個相變的過程.土壤水分在液相和氣相之間相互轉(zhuǎn)換時，必然伴隨著不同形式的水分（液態(tài)和氣態(tài)）及熱量的流動.因此，在土壤－大氣連續(xù)體中分析非飽和土中水分的運動時，液態(tài)水的流動，氣態(tài)水的流動和熱量的流動需要同時考慮.本文以Penman－Wilson蒸發(fā)模型預(yù)測的實際蒸發(fā)量為水流量邊界條件，利用較易獲取的土體基本物性指標(biāo)（顆粒分布、干密度、比重）預(yù)測出相應(yīng)水力參數(shù)和熱特性參數(shù)，同時結(jié)合基于Darcy定律和Fick定律而建立的非飽和土水－熱遷移模型，分析了蒸發(fā)條件下非飽和土壤水分的運動規(guī)律，為自然環(huán)境下非飽和土水分遷移問題的探究提供參考.

1 數(shù)值模型的建立

1.1 土壤水分運動控制方程

為了建立描述土體中液態(tài)水、氣態(tài)水和熱量的運動方程，需引入如下假定［6］：土體的力學(xué)行為可以用其自身土顆粒、空氣和水形成的連續(xù)單元體所描述，忽略由于滲透壓梯度引起的液態(tài)水流動，液態(tài)水和空氣的滲透系數(shù)僅僅和含水量或飽和度相關(guān)，不考慮滲透系數(shù)和吸力關(guān)系間的滯回，單元體的液態(tài)水和水蒸氣局部熱力學(xué)平衡始終成立，土體中任一點的溫度在水的融點之上和沸點之下，水中氣體的溶解不加以考慮.基于以上假設(shè)，建立了單元體中水的質(zhì)量和能量守恒的方程，其一維形式表述如下：

方程（1）中有總水頭和水蒸氣分壓兩個變量，但這兩個變量并不是獨立的，他們之間的關(guān)系可以通過Edlefsen和Anderson［7］提出的以下關(guān)系式來描述：

式中，Pvs為純水飽和蒸氣壓力（kPa）；hrair為空氣相對濕度；w為水蒸汽克分子量（0.018kg／mol）；R為通用氣體常數(shù)（8.314J／mol·K）；T為絕對溫度（K）.

1.2 邊界條件

在利用上述方程進行飽和－非飽和滲流計算時，必須給定計算模型的流量邊界條件和溫度邊界條件.

根據(jù)數(shù)值模擬的特點，可選用土表實際蒸發(fā)量作為計算模型的流量邊界條件，其值可用Penman－Wilson公式［8］計算，該公式是Wilson通過考慮土表相對濕度變化，同時引入總吸力，把土壤蒸發(fā)作為土壤濕度的連續(xù)函數(shù)而建立的：

式中，E為蒸發(fā)量（mm／day）；Γ為飽和蒸氣壓和溫度關(guān)系曲線在溫度為Ta時的斜率；Q為土表面等效凈輻射量（mm／d）；v為濕度常數(shù)，一般取0.66hPa／℃，Ea＝f（u）ea（B－A）；f（u）是風(fēng)函數(shù)，f（u）＝0.35×（1＋0.15u），u為風(fēng)速（m／s）；ea為蒸發(fā)面上方空氣中的水蒸氣分壓（kPa）；B為空氣相對濕度的倒數(shù)，即es／ea；es為空氣飽和蒸氣壓（kPa）；A為土表面相對濕度的倒數(shù)，即1／hr；當(dāng)土表飽和時，即土表相對濕度為100%時，該公式演變?yōu)?/p>

該式即為Penman推出的計算飽和蒸發(fā)量的公式，同時可以用于計算水面的潛在蒸發(fā)量.

土體和大氣間的熱量交換發(fā)生在土壤表面，因此地表土溫度可以作為計算模型溫度的邊界條件，其值可用Wilson提出的下式［6］換算：

式中，Ts為地表土的溫度（K）；Ta為土表空氣溫度（K）.

1.3 參數(shù)確定

為了求解上述模型，必需給定材料參數(shù)及蒸發(fā)邊界條件，如土水特征曲線、滲透性函數(shù)、水蒸氣擴散系數(shù)、熱特性參數(shù)和土表實際蒸發(fā)量.為避免傳統(tǒng)試驗測試手段周期長、費用昂貴的不足，本文根據(jù)相關(guān)氣象條件（溫度、相對濕度、風(fēng)速）和Penman－Wilson公式預(yù)測土表實際蒸發(fā)量，利用土壤基本物性指標(biāo)（顆粒分布、干密度、比重）預(yù)測上述參數(shù)，具體方法簡述如下：利用顆粒分布曲線和相關(guān)物理性質(zhì)指標(biāo)，根據(jù)Arya－Paris［9］模型預(yù)測出不同含水量所對應(yīng)的基質(zhì)吸力值.在此基礎(chǔ)上，利用Fredlund ＆Xing［10］模型對含水量－基質(zhì)吸力的離散點進行擬合便得到完整土水特征曲線.

土體在非飽和情況下，其滲透系數(shù)隨著其自身含水量的變化而變化.在已知土體顆粒分布、比重及干密度的基礎(chǔ)上，利用Kozeny－Carman［11］模型預(yù)測出土樣飽和滲透系數(shù)，同時結(jié)合Fredlund ＆Xing模型預(yù)測所得的土水特征曲線相關(guān)參數(shù)，便可得到完整的滲透性函數(shù)曲線.

根據(jù)Philip和deVries［12］模型的描述，在土體孔隙率（可由比重和干密度求得）已知的前提下，水蒸氣擴散系數(shù)可由其與土體含水量和環(huán)境溫度的函數(shù)關(guān)系求得.

土體的熱特性參數(shù)主要指土壤的體積比熱和熱傳導(dǎo)率，其值取決于土體中固相物質(zhì)（土顆粒）、液相物質(zhì)（孔隙水）和氣相物質(zhì)（孔隙氣）所占的體積比，以及他們各自的熱容和熱傳導(dǎo)率.水蒸氣的潛化熱是指單位質(zhì)量的液體轉(zhuǎn)變?yōu)橄嗤瑴囟鹊乃魵鈺r吸收的熱量，其值可以利用一個和溫度呈線性的方程描述［6］：Lv＝4.186×103×（607－0.7T）.

2 模型驗證

2.1 幾何模型求解及材料參數(shù)預(yù)測

為了驗證上述方程及參數(shù)預(yù)測方法的可靠性，將該理論用于模擬Wilson土柱蒸發(fā)試驗實測情況.采用有限單元法來求解偏微分方程中含水量及溫度變化率，利用Vadose／w建模并實現(xiàn)數(shù)值模擬.

試驗用土為均質(zhì)砂土，其顆粒分布曲線如圖1所示，比重為2.67，土柱模型干容重為1.64g／cm3，高度為30cm，直徑為10cm.總共劃分300個四邊形有限元單元，沿高度劃分為30個等份，沿直徑劃分10等份.

圖1 顆粒分布曲線

由土壤基本物性指標(biāo)（顆粒分布、干密度、比重），結(jié)合上述相關(guān)參數(shù)的預(yù)測方法可獲取模型計算所需參數(shù).其完整的土水特征曲線如圖2所示，滲透性函數(shù)相關(guān)參數(shù)見表1，土樣熱傳導(dǎo)率和體積比熱容與含水率關(guān)系曲線如圖3所示.

圖2 土水特征曲線

圖3 熱傳導(dǎo)率函數(shù)曲線和體積比熱容函數(shù)曲線

表1 Fredlund ＆Xing模型預(yù)測土水特征曲線和滲透性函數(shù)相關(guān)參數(shù)

影響蒸發(fā)的主要外部氣象因素主要有風(fēng)速、溫度、濕度.因此，為了求解瞬時水流量邊界和溫度邊界，需要給定實時地氣象參數(shù).Wilson土柱蒸發(fā)試驗條件為：平均環(huán)境溫度為38℃，相對濕度接近10%，風(fēng)速為0m／s.在此環(huán)境條件下土柱模型從飽和含水量開始持續(xù)蒸發(fā)40d，這一過程當(dāng)中，土柱沒有任何水分的補充，土柱水分的蒸發(fā)量即為土柱質(zhì)量的減少量，其值采用電子天平即時獲取.

2.2 計算結(jié)果分析

土柱蒸發(fā)試驗蒸發(fā)量實測值與計算值對比結(jié)果如圖4所示.從圖4中可以看出，在前4d，土柱蒸發(fā)量實測值與計算值均穩(wěn)定在7mm／d的潛在蒸發(fā)量附近.但從第4天末開始，土柱蒸發(fā)量迅速降低，在不到兩周的時間里，其值從7mm／d降低到1mm／d.隨后，土柱的蒸發(fā)量緩慢降低并趨于平穩(wěn)，在超過兩周的時間里，其值從1mm／d降低到0.5mm／d.通過對比不難發(fā)現(xiàn)，在為期40d的蒸發(fā)試驗中，蒸發(fā)量實測值和計算值的變化趨勢基本一致.同時，試驗及數(shù)值模擬結(jié)果也驗證了非飽和土蒸發(fā)呈現(xiàn)3個階段：穩(wěn)定階段、減速階段和殘余階段.

圖4 蒸發(fā)量實測值與計算值比較

圖5為蒸發(fā)試驗的不同時刻土柱縱剖面體積含水量分布情況.從圖中可以看出，在蒸發(fā)試驗進行到第4天末，土柱不同深度處體積含水量均低于初始飽和含水量，但此時不同深度的含水量降低幅度較一致，因此在縱剖面上并沒有明顯的體積含水量梯度出現(xiàn).在這一時刻，剖面含水量計算值和實測值吻合較好.但隨著蒸發(fā)試驗的持續(xù)進行，剖面含水量繼續(xù)降低，在第10天末，土柱表層體積含水量接近于0，而表層以下10cm處的體積含水量為19%，這標(biāo)志著在土柱表面有含水量極低的干燥土層的形成.此后，干燥面不斷的向下推移，在第40天末，干燥面推移至土表以下10cm處，此時的干燥土層厚度達到10cm.通過對比可以發(fā)現(xiàn)，在蒸發(fā)進行到第10天末及第40天末，土柱縱剖面含水量的計算值和實測值出現(xiàn)了明顯的差異，特別是在干燥土層以下相對濕潤的置，出現(xiàn)這一現(xiàn)象的原因有可能是隨著含水量的降低，預(yù)測的土體滲透系數(shù)減小得過快，因而在某一含水量時，輸水能力的預(yù)測值小于其真實值，在模擬計算時，同一單元節(jié)點處液態(tài)水流量小于實際情況，隨著深度的增加，這一現(xiàn)象更為突出，這也從側(cè)面說明了土體滲透系數(shù)函數(shù)對干燥土層以下土壤水分的運動起決定性的作用.而在干燥土層內(nèi)部，土體滲透系數(shù)因含水量的急劇降低而趨近于0.

圖5 不同深度處體積含水量實測值與計算值比較

圖6～7分別為1cm深度處單元節(jié)點和8cm深度處節(jié)點水流量情況.從圖中可以看出，在蒸發(fā)前期各點處水分的遷移以液態(tài)形式流動為主.但隨著含水量的降低，單元點處的滲透系數(shù)不斷減小，從而導(dǎo)致該點處液態(tài)水通量迅速降低.從某一時刻起，水蒸氣擴散成為該點處水分遷移的主要形式.以1cm深度處節(jié)點為例，在前3天，液態(tài)水流量維持在7×10－8m／s附近，水蒸汽流量接近0（此時土體飽和度高，水蒸氣擴散系數(shù)?。?從第4天起，液態(tài)水流量開始快速下降，而水蒸氣流量隨著水蒸氣擴散系數(shù)的增大而增加.在第6天，液態(tài)水流量降低至5.5×10－9m／s，而此時水蒸氣流量為2.0×10－8m／s，水蒸氣的流動占據(jù)主導(dǎo)地位.

圖6 1cm深度處節(jié)點液態(tài)水和水蒸氣隨時間變化情況

圖7 8cm深度處節(jié)點液態(tài)水和水蒸氣通量隨時間變化情況

圖8為不同時刻土體縱剖面溫度分布情況.從圖8可以看出在蒸發(fā)初期，土體溫度迅速降低，在第1天末，距離土表1cm深度處土層溫度為30.6℃，距離土表25cm深度處土層的溫度為35.7℃，可見靠近土表層位置處溫度降低的幅度大于底層的.但隨著時間的推移，土體不同深度處的溫度均開始回升，在第42天，土體縱剖面溫度重新回到了初始值附近.出現(xiàn)上述現(xiàn)象的原因是液態(tài)水轉(zhuǎn)變成水蒸氣需要消耗熱量.在試驗初期，土表保持7mm／d潛在蒸發(fā)強度，土表蒸發(fā)需要消耗大量的熱量，此時汽化潛熱成為單元體熱量變化的主要形式.而在第4天后，隨著蒸發(fā)量的迅速降低，所需消耗的熱量減小，熱傳導(dǎo)在單元體熱量變化中的作用逐漸凸顯出來.因此，在該作用下土體溫度逐漸變得均勻并再次接近于環(huán)境溫度.同時還不難發(fā)現(xiàn)，土柱縱剖面上最低點的溫度和干土層的向下推移的位置較一致，這也從另一個角度說明了蒸發(fā)面的位置隨著蒸發(fā)的進行而不斷向土體內(nèi)部移動.

圖8 不同深度處溫度實測值與計算值比較

3 結(jié) 語

利用土體基本物性指標(biāo)預(yù)測了非飽和土水－熱遷移模型計算所需的水力參數(shù)、熱特性參數(shù)，同時引入氣象環(huán)境參數(shù)預(yù)測了土表實際蒸發(fā)量和溫度值，模擬了蒸發(fā)條件下非飽和土中水分及熱量的遷移規(guī)律.

通過該理論預(yù)測得到的非穩(wěn)定蒸發(fā)量、縱剖面含水量分布、縱剖面溫度分布與Wilson室內(nèi)蒸發(fā)試驗實測吻合較好.這說明采用上述模型及相關(guān)參數(shù)預(yù)測方法來模擬蒸發(fā)條件下土壤水分的遷移是可行的.

模擬結(jié)果表明滲透系數(shù)的急劇降低導(dǎo)致干土層的出現(xiàn)，從而導(dǎo)致土壤水分的運動形式發(fā)生了根本的轉(zhuǎn)變，即在干土層內(nèi)部水蒸氣擴散成為土壤水分向上運動的主要形式.而此時在水力梯度下的滲透作用幾乎可以忽略，土壤的實際蒸發(fā)面也不再是土壤表面.

［1］ 柯尊敬，黃紹鏗.膨脹土工程性質(zhì)的研究總報告［R］.南寧：廣西大學(xué)土木系，1995.

［2］ 王鐵行，李 寧，謝定義.土體水熱力耦合問題研究意義、現(xiàn)狀及建議［J］.巖土力學(xué)，2005，26（3）：487－492.

［3］ 陳建斌，孔令偉，趙艷林，等.非飽和土的蒸發(fā)效應(yīng)與影響因素分析［J］.巖土力學(xué)，2007，28（1）：36－40.

［4］ Thomas H R，Rees S W.The Numerical Simulation of Seasonal Soil Drying in an Unsaturated Clay Soil［J］.International Journal for Numerical＆Analytical Methods Geomechanics，1993，17（1）：119－132.

［5］ 吳宏偉，陳守義，龐宇威.雨水入滲對非飽和土坡穩(wěn)定性影響的參數(shù)研究［J］.巖土力學(xué)，1999，20（1）：1－14.

［6］ Wilson G W.Soil Evaporation Fluxes for Geotechnical Engineering Problems［D］.Saskatoon：University of Saskatoon，1990.

［7］ Edlefsen N E，Anderson A B C.Thermodynamics of Soil Moisture［J］.Hilgardia，1943，15：31－39.

［8］ Wilson G W，F(xiàn)redlund D G，Barbour S L.Coupled Soilatmosphere Modeling for Soil Evaporation［J］.Canadian Geotechnical Journal，1994，31：151－161.

［9］ Arya L M，Paris J FA.Physicoempirical Model to Predict the Soil Moisture Characteristic from Particle－Size Distribution and Bulk density Data［J］.Soil Science Society of America Journal，1981，45：1023－1030.

［10］Fredlund D G，Xing A.Equations for the Soil－Water Characteristic Curve［J］.Canadian Geotechnical Journal，1994，31：521－532.

［11］Carman P C.Permeability of Saturated Sands，Soils and Clays［J］.Journal of Agricultural Science，1939，29：26.

［12］De Vries D A.Thermal Properties of Soils［M］.Amsterdam：North－Holland Publishing Company，1963.