蒸發條件下非飽和土水分遷移的數值模擬

熊 義 高 抗

（1.三峽大學水利與環境學院，湖北宜昌 443002；2.三峽大學三峽庫區地質災害教育部重點實驗室，湖北宜昌 443002）

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

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

1 數值模型的建立

1.1 土壤水分運動控制方程

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

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

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

1.2 邊界條件

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

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

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

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

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

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

1.3 參數確定

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

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

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

土體的熱特性參數主要指土壤的體積比熱和熱傳導率，其值取決于土體中固相物質（土顆粒）、液相物質（孔隙水）和氣相物質（孔隙氣）所占的體積比，以及他們各自的熱容和熱傳導率.水蒸氣的潛化熱是指單位質量的液體轉變為相同溫度的水蒸氣時吸收的熱量，其值可以利用一個和溫度呈線性的方程描述［6］：Lv＝4.186×103×（607－0.7T）.

2 模型驗證

2.1 幾何模型求解及材料參數預測

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

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

圖1 顆粒分布曲線

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

圖2 土水特征曲線

圖3 熱傳導率函數曲線和體積比熱容函數曲線

表1 Fredlund ＆Xing模型預測土水特征曲線和滲透性函數相關參數

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

2.2 計算結果分析

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

圖4 蒸發量實測值與計算值比較

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

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

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

圖6 1cm深度處節點液態水和水蒸氣隨時間變化情況

圖7 8cm深度處節點液態水和水蒸氣通量隨時間變化情況

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

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

3 結 語

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

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

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

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

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

［3］ 陳建斌，孔令偉，趙艷林，等.非飽和土的蒸發效應與影響因素分析［J］.巖土力學，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］ 吳宏偉，陳守義，龐宇威.雨水入滲對非飽和土坡穩定性影響的參數研究［J］.巖土力學，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，Fredlund 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.