土壤水-汽-熱耦合運移的數(shù)值模擬研究與驗證

喬星華，董曉華，孫 媛，劉旋旋

（1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌443002；2.水資源安全保障湖北省協(xié)同創(chuàng)新中心，武漢430072；3.廣東省水文局肇慶水文分局，廣東肇慶526060；4.廣東省水文局清遠(yuǎn)水文分局，廣東清遠(yuǎn)511500）

0 引 言

非飽和土壤水分運動一直是水文學(xué)、土壤科學(xué)和農(nóng)業(yè)領(lǐng)域關(guān)注的熱點問題，土壤中的水-汽-熱耦合運移理論已廣泛應(yīng)用于全球氣候變化和生產(chǎn)實踐當(dāng)中[1]。土壤是一個非均質(zhì)、多相、分散和多孔的系統(tǒng)，水、汽、熱在土壤中存在并且以土壤為傳導(dǎo)介質(zhì)進(jìn)行變化，其中非飽和土壤由固相、液相和氣相組成[2，3]。作為土壤水資源補給的重要途徑，土壤入滲過程近年來被廣泛應(yīng)用于土壤水資源評價、水文模型構(gòu)建等相關(guān)領(lǐng)域中[4]。

國內(nèi)外對于土壤水分運動的研究主要考慮的是溫度的影響，一些學(xué)者對土壤水-汽-熱耦合運移規(guī)律做出了研究，楊金忠[6]等利用參數(shù)的經(jīng)驗值對水、汽、熱耦合作用的理論和數(shù)值計算進(jìn)行了研究。水建高[7]等利用非等溫的水-汽-熱耦合模型模擬了松土條件下土壤水分及溫度分布，得到地表處必須考慮溫度對水汽流動的影響的結(jié)論。曾亦鍵[8]基于PdV 理論，考慮了土壤內(nèi)部空氣，建立了蒸發(fā)條件下土壤水-汽-熱-空氣耦合模型，分析了土壤水汽運移規(guī)律，從理論上完善了非飽和帶土壤水-汽-熱耦合模型研究。高萬德[9]利用Hydrus 建立水-汽-熱耦合運移模型，對農(nóng)田灌溉條件下大埋深包氣帶的土壤水-汽-熱運移規(guī)律進(jìn)行研究。李吉祥等[10]認(rèn)為包氣帶水運移涉及液態(tài)水和水汽，而土壤水分運動涉及的條件復(fù)雜，求解困難，數(shù)值方法則是求解包氣帶水分運移的主要方法。MALIK[11]提出非飽和凍土中水-汽-熱耦合模型，得出水汽的遷移依賴于土體內(nèi)溫度梯度和初始含水量的結(jié)論。

基于以上研究可以發(fā)現(xiàn)，目前很少有針對砂壤土土壤水-汽-熱耦合運移的研究。本文開展了室內(nèi)一維土柱試驗，分降雨和蒸發(fā)兩個階段，觀測了土柱內(nèi)不同深度處的含水率、溫度、氣壓的變化過程，蒸發(fā)階段的蒸發(fā)速率，以及包括空氣溫濕度、風(fēng)速、大氣壓強(qiáng)等環(huán)境氣象因子；應(yīng)用土壤水-汽-熱耦合模型，使用上述觀測數(shù)據(jù)對模型參數(shù)進(jìn)行率定，對降雨和蒸發(fā)兩種情景下的土壤水-汽-熱耦合運移進(jìn)行模擬，利用試驗得到的土壤溫度和含水率對模擬值進(jìn)行驗證。本研究所得結(jié)論可為砂壤土土壤水-汽-熱耦合運移提供理論基礎(chǔ)，同時為農(nóng)業(yè)灌溉領(lǐng)域提供理論參考依據(jù)。

1 材料與方法

1.1 試驗材料

土樣取自三峽大學(xué)校園，選取多個地形點進(jìn)行采樣，清除采樣點處表面覆蓋物并用環(huán)刀[12-17]取深度為5～15 cm 間的土樣，將所有土樣放在潔凈的塑料板上均勻混合，經(jīng)過風(fēng)干、磨碎并過2 mm 標(biāo)準(zhǔn)篩網(wǎng)后密封保存[18]。本試驗采用篩分法與密度計法相結(jié)合分析土壤顆粒，確認(rèn)試驗土樣為砂壤土，采用環(huán)刀法測得土樣的容重為1.409 g/cm3，由達(dá)西公式得到該土壤的飽和滲透系數(shù)為Ks=0.355 mm/min。供試土樣的基本物理性質(zhì)見表1。

表1 供試土樣的基本物理性質(zhì)Tab.1 Basic physical properties of the test soil sample

1.2 土柱試驗方法

土柱實驗于2019年6月在室內(nèi)進(jìn)行，將處理后的土樣放到土柱筒中，土柱筒的高度為80 cm，直徑為14 cm。采用了模擬人工降雨蒸發(fā)的方法。

試驗分為降雨和蒸發(fā)兩個階段。第一階段為降雨階段，降雨強(qiáng)度為12 mm/h，采用馬氏瓶恒定供水，首先進(jìn)行2 h 降雨，在降雨之前對降雨裝置進(jìn)行率定，降雨強(qiáng)度是通過調(diào)整馬氏瓶的高度來控制的[19]；第二階段為蒸發(fā)階段，蒸發(fā)階段的光照輻射由距離土柱頂端5～10 cm 的遠(yuǎn)紅外燈提供，為了測量土壤的蒸發(fā)速率，將土柱置于高精度電子天平之上，蒸發(fā)時間為4 h。分別在土柱距離頂端3、8、18、28、48 cm 處放置TDR 探頭、溫度探頭和氣壓探頭。降雨和蒸發(fā)階段周期均為48 h，且降雨和蒸發(fā)均在距離開始監(jiān)測各項指標(biāo)的第10 h 進(jìn)行。試驗數(shù)據(jù)的測量由TDR 傳感器（含水率、溫度）、氣壓傳感器（空氣壓強(qiáng)、土壤氣壓）、手持式氣象站（空氣溫度、空氣濕度、風(fēng)速等）、高精度電子天平測得（蒸發(fā)速率）。

1.3 數(shù)值模擬方法

由室內(nèi)土柱模型監(jiān)測并記錄得到降雨和蒸發(fā)階段的環(huán)境氣象因子、土壤含水率、土壤氣壓和土壤溫度等。STEMMUS模型[8]是曾亦鍵提出的模擬不飽和土壤中液態(tài)水、水蒸氣、干燥空氣和傳熱耦合的程序，由數(shù)值模型控制方程（包括土壤水分運動方程、干空氣平衡方程、能量平衡方程）構(gòu)成，考慮溫度影響的土壤水分特征曲線模型影響著模型的控制方程。土壤水分特征曲線考慮常溫（20 ℃）影響，修改模型中的土壤的殘余含水率θr、飽和含水率θs和經(jīng)驗擬合參數(shù)α和n的經(jīng)驗參數(shù)，利用試驗結(jié)果對STEMMUS 的參數(shù)進(jìn)行率定，得到精度較高的模擬結(jié)果。

1.3.1 數(shù)值模型控制方程

（1）土壤水分運動方程。STEMMUS 模型中對土壤水分運動公式的改進(jìn)形式為：

式中：ρL、ρV分別為液態(tài)水、水汽密度，kg/m3；θL為體積含水率，m3/m3；θV為土壤空氣體積，m3/m3；qL、qV分別為液態(tài)水、水汽的通量，kg/(m2·s)；t為時間，s；z為垂向坐標(biāo)軸，m。

（2）干空氣平衡方程。STEMMUS 模型中考慮了土壤中空氣的傳輸及影響過程，在非飽和土壤中空氣運移的驅(qū)動力為空氣密度梯度和氣壓梯度[8]。其表達(dá)式為：

式中：ρda為干空氣密度，kg/m3；Hc為Henry常數(shù)，對于空氣，在標(biāo)準(zhǔn)大氣壓、25 ℃條件下，其值為0.02；qa為干空氣通量，kg/(m2·s)；

（3）能量平衡方程。溫度在土壤中的傳輸?shù)男问街饕▊鲗?dǎo)、輻射和對流，非飽和土壤能量傳輸包括土壤固相、液相、干空氣和水汽的熱運移及濕潤熱[20]。非飽和土壤中的能量傳輸方程可表達(dá)為：

式中：λL、λs、λg為固、液和氣的熱傳導(dǎo)系數(shù)，W/(m·℃)；θs為土壤固體體積分?jǐn)?shù)；θa為土壤空氣體積分?jǐn)?shù)，θa=θV=ε-θL；cs、cL、ca、cV為固、液、水汽、氣的比熱，J/(kg·℃)；T為參考溫度，℃；ρs為土壤固體密度，kg/m3；L0為參考溫度T時的蒸發(fā)潛熱，J/kg；W為局部濕潤熱，J/kg，λeff為包含固、液、氣的有效熱傳導(dǎo)系數(shù)，W/(m·℃)。

1.3.2 溫度影響下的土水特征曲線模型

本研究通過對VG 模型的經(jīng)驗擬合參數(shù)n和α、土壤的殘余含水率θr、飽和含水率θs進(jìn)行修改使模擬結(jié)果達(dá)到最優(yōu)化。由于試驗在室內(nèi)進(jìn)行，室內(nèi)溫度在20 ℃左右，因此當(dāng)考慮溫度為20 ℃（工程上的常溫狀態(tài)）時，對影響土壤水分特征曲線的土壤殘余含水率θr、飽和含水率θs和經(jīng)驗擬合參數(shù)α和n進(jìn)行修改。

美國學(xué)者Van Genuchten 在1980年提出了描述土壤水分特征曲線的VG模型[22]，其表達(dá)式為：

式中：θr為殘余含水率；θs為飽和含水率；α、m和n為經(jīng)驗擬合參數(shù)。

考慮溫度影響時，土壤殘余含水率θr、飽和含水率θs和經(jīng)驗擬合參數(shù)α經(jīng)驗公式如公式(5)～(7)所示，其中溫度T為20 ℃（模型中用開氏溫度293.15 K計算）。

由公式(5)～(7)可以將VG模型改寫成式(8)：

式中：αT為20 ℃下的經(jīng)驗擬合參數(shù)；θrT為20 ℃下的殘余含水率；θsT為20 ℃下的飽和含水率。

1.3.3 初始條件和邊界條件

(1)初始條件。土壤初始條件是在時間t=0 時，土壤在0 -L土柱總長度為L）時的土壤測得的各項參數(shù)，包括土壤含水率、土壤溫度、土壤氣壓，大氣溫度、濕度、空氣壓強(qiáng)等。

(2)上邊界條件。本文中土柱上表面為開放型表面，選擇第二類邊界條件，其表達(dá)形式為：

式中：E為時間為t時的蒸發(fā)速率，kg/(m2·sm2·s)；P為時間為t時的降水速率，m/s；t為時間，s。

土壤蒸發(fā)量由式(10)推求，主要是根據(jù)大氣阻力和土壤對水汽的阻力計算：

式中：ρVS為土層表面的水汽密度，kg/m3；ρVa為大氣中的水汽密度，kg/m3；ra為大氣阻力，s/m；rs為土層表面對水汽的阻力，s/m；f是常數(shù)，值為0.41；U為風(fēng)速，m/s，Zm為測量的風(fēng)速所在高度，m；d為零平面位移，m，如果為裸地，其值為0；Zom為表面粗糙長度（動量通量），取值0.001 m；ψsm為大氣穩(wěn)定性修正系數(shù)（動量通量）；Zoh表示表面粗糙長度（熱通量），取值0.001 m；ψsh為大氣穩(wěn)定性修正系數(shù)（熱通量）；rsl為水分子由水層表面擴(kuò)散到空氣中所遇到的阻力，取值10 s/m；為擬合系數(shù)，取值35.63；θmin為土壤最小含水率，取值0.15 m3/m3；θsur為土壤表層含水率，m3/m3。

空氣方程邊界條件為第一類邊界條件，表示為：

式中：Pg為時間為t時的空氣壓強(qiáng)，Pa；Γ為大氣-土壤邊界。

溫度傳導(dǎo)方程的邊界條件為第一類邊界條件，表示為：

式中：Ts為時間為t時的地表溫度，℃。

(3)下邊界條件。本次試驗采用的土柱筒下設(shè)有排氣閥，當(dāng)土壤底部含水率達(dá)到飽和時，可以自由排出。對于土壤水分運動、土壤溫度傳導(dǎo)、空氣傳導(dǎo)的下邊界條件可以寫成公式(15)～(17)的形式：

1.3.4 時空離散

實驗選用均質(zhì)砂壤土，土壤深度為70 cm，將土柱離散為25 個節(jié)點，由上至下節(jié)點間距依次增大。土壤表層的分層厚度為0.25～0.5 cm，土壤中上層的分層厚度為1 cm，土壤中下層分層厚度為2 cm，土壤下層的分層厚度為5～10 cm，將時間步長設(shè)為1 s。

1.3.5 模型基本參數(shù)

模型的基本參數(shù)參考Rosetta 數(shù)據(jù)庫的砂壤土水力特征參數(shù)[23]，修改后的經(jīng)驗參數(shù)，模型基本參數(shù)如表2所示。

表2 模型基本參數(shù)Tab.2 Basic parameters of the model

1.3.6 模型求解

1.3.7 模型精度驗證

選用平均相對誤差(AVRE)、相對均方根誤差(RRMSE)作為誤差指標(biāo)對降雨和蒸發(fā)階段含水率和溫度的實測值和模擬值進(jìn)行誤差分析[24]。當(dāng)RRMSE=0時，模擬值與實測值完全吻合，RRMSE＜10%時，模擬值與實測值的一致性很好，10%～20%為比較好，20%～30%表明模擬效果一般。AVRE值越接近于1，說明模擬效果越好[25]。

式中：Ci為模擬值；Mi為實測值；Nw為總的樣本觀測數(shù)。

2 結(jié) 果

圖1為降雨和蒸發(fā)階段空氣溫度的變化情況。

分別對降雨和蒸發(fā)階段的土壤水分和土壤溫度同時進(jìn)行模擬，得到如下結(jié)果：

圖2為降雨階段土壤含水率和溫度的模擬圖，可以看出，降雨開始3～28 cm深度的土壤含水率明顯增加，降雨結(jié)束后隨著土壤水分向下層傳播，各層的土壤含水率開始變小，減小到一定程度，各層土壤含水率接近于一個恒定值，而48 cm 處的含水率上升到一定值后保持恒定不變；降雨階段大氣溫度高于土壤的溫度，在這個時間段內(nèi)各層土壤的溫度呈上升趨勢，這種土壤溫度升高隨著降雨入滲而向土柱深處傳播。

降雨階段各層土壤含水率和溫度的總體模擬效果較好。其中土壤含水率隨著深度的增加，模擬精度變高，土壤含水率的實測值在降雨的再分布階段略高于實測值。土壤溫度在3 cm 處、8 cm 處的模擬效果相對較好，其余各層的溫度模擬值與實測值均有偏差，并且模擬精度隨著深度的增加而減小。

自改革開放以來，中共中央1982年至1986年連續(xù)5年發(fā)布以“三農(nóng)”(農(nóng)業(yè)、農(nóng)村、農(nóng)民)為主題的中央“一號文件”，2004年至2018年又連續(xù)15年發(fā)布以“三農(nóng)”為主題的中央“一號文件”，不斷指導(dǎo)“三農(nóng)”改革創(chuàng)新。“一號文件”某種意義上已成為中共中央重視“三農(nóng)”工作的專有名詞。

表3為降雨階段土壤含水率和溫度的模擬結(jié)果，可以看出，各層土壤含水率和溫度的平均相對誤差都接近于1。對于含水率的模擬，18 cm 以下深度的相對均方根誤差都小于10%，模擬效果很好；8 cm 處的相對均方根誤差在10%～20%之間，模擬效果比較好；而3 cm 處的相對均方根誤差在20%～30%之間，模擬效果一般。對于溫度的模擬，各層的相對均方根誤差均小于10%，模擬效果很好。

圖3為蒸發(fā)階段土壤含水率和溫度的模擬圖，通過監(jiān)測空氣溫度的變化，空氣溫度從第5 h 開始升高，導(dǎo)致土壤上層含水率開始減小。土壤其他層含水率變化規(guī)律類似，深度越大，土壤含水率變化幅度逐漸減小，中下層土壤含水率幾乎無變化。加熱開始時，土壤表層3 cm 處溫度迅速升高并逐漸達(dá)到最大值。加熱結(jié)束時，土壤表層3 cm 處溫度迅速降低，最終溫度與室溫達(dá)到一致。土壤8 cm 處的溫度在10.5 h 左右開始上升，在第14.5 h 處升至最大值；土壤18 cm 處的溫度在第12 h開始上升，在第16 h處達(dá)到最大值，中上層各層溫度最終都減小到室溫，且溫度減小的趨勢逐漸變緩；28 cm 處的溫度有較小的變化；而土壤48 cm處的溫度幾乎沒有變化。

表3 降雨階段含水率和溫度模擬結(jié)果Tab.3 Moisture content and temperature simulation results duringrainfall

這一階段土壤各層的含水率模擬效果較好，且下層土壤含水率的模擬精度高于上層的模擬精度。3 cm 和8 cm 處溫度模擬值與實測值非常接近，模擬效果較好，但隨深度的增加，溫度的模擬值偏差也在增加，由圖3（b）、圖3（d）、圖3（f）、圖3（h）和圖3（j）可以看出28 cm 和48 cm 處在蒸發(fā)開始后的擬合效果很差，在17 h后模擬值均大于實測值。

表4為蒸發(fā)階段土壤含水率和溫度的模擬結(jié)果，可以看出，土壤含水率的模擬值和上層土壤溫度的平均相對誤差都接近于1。對于含水率的模擬，各層相對均方根誤差均小于10%，模擬值與實測值一致性很好。對于溫度的模擬，中上層的相對均方根誤差均小于10%，模擬效果很好；但28 cm 處和48 cm 處的相對均方根誤差在20%～30%之間，模擬效果一般。

3 結(jié) 論

本文開展了室內(nèi)一維土柱試驗，分降雨和蒸發(fā)兩個階段，觀測了土柱內(nèi)不同深度處的含水率、溫度、氣壓的變化過程，蒸發(fā)階段的蒸發(fā)速率，以及包括空氣溫濕度、風(fēng)速、大氣壓強(qiáng)等環(huán)境氣象因子；然后應(yīng)用STEMMUS 模型，使用上述觀測數(shù)據(jù)對模型參數(shù)進(jìn)行率定，獲得了較好的模擬精度。主要結(jié)論如下：

（1）降雨階段，土壤溫度受大氣溫度的影響較大。本實驗中，大氣溫度高于土壤溫度，雨滴進(jìn)入土壤后，導(dǎo)致土壤溫度升高，這種土壤溫度升高隨著降雨入滲而向土柱深處傳播。

（2）蒸發(fā)階段，在光照輻射的作用下，土柱內(nèi)各層土壤的含水率緩慢下降；土柱表層土壤溫度迅速上升，熱量向土柱深層傳播，使深層土壤溫度逐次上升。

（3）降雨階段含水率和溫度的模擬效果優(yōu)于蒸發(fā)階段。降雨階段含水率、溫度模擬精度都“比較好”或“很好”，只有3 cm 處的表層土壤含水率模擬精度為“一般”（RRMSE為21.99%），原因是表層土壤的質(zhì)地不均勻所致；蒸發(fā)階段含水率和溫度大部分深度模擬精度“很好”，只有深層（28 及48 cm 處） 的溫度模擬精度為“一般”（RRMSE為21.18%及22.92%），原因應(yīng)該是土柱下部會受環(huán)境溫度顯著影響所致。

表4 蒸發(fā)階段含水率和溫度模擬結(jié)果Tab.4 Moisture content and temperature simulation results in the evaporation stage