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

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

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

0 引 言

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

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

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

1 材料與方法

1.1 試驗(yàn)材料

土樣取自三峽大學(xué)校園，選取多個(gè)地形點(diǎn)進(jìn)行采樣，清除采樣點(diǎn)處表面覆蓋物并用環(huán)刀[12-17]取深度為5～15 cm 間的土樣，將所有土樣放在潔凈的塑料板上均勻混合，經(jīng)過風(fēng)干、磨碎并過2 mm 標(biāo)準(zhǔn)篩網(wǎng)后密封保存[18]。本試驗(yàn)采用篩分法與密度計(jì)法相結(jié)合分析土壤顆粒，確認(rèn)試驗(yàn)土樣為砂壤土，采用環(huán)刀法測(cè)得土樣的容重為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 土柱試驗(yàn)方法

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

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

1.3 數(shù)值模擬方法

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

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

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

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

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

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

（3）能量平衡方程。溫度在土壤中的傳輸?shù)男问街饕▊鲗?dǎo)、輻射和對(duì)流，非飽和土壤能量傳輸包括土壤固相、液相、干空氣和水汽的熱運(yùn)移及濕潤(rùn)熱[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時(shí)的蒸發(fā)潛熱，J/kg；W為局部濕潤(rùn)熱，J/kg，λeff為包含固、液、氣的有效熱傳導(dǎo)系數(shù)，W/(m·℃)。

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

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

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

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

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

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

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

1.3.3 初始條件和邊界條件

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

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

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

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

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

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

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

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

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

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

1.3.4 時(shí)空離散

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

1.3.5 模型基本參數(shù)

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

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

1.3.6 模型求解

1.3.7 模型精度驗(yàn)證

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

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

2 結(jié) 果

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

（1）降雨階段，土壤溫度受大氣溫度的影響較大。本實(shí)驗(yàn)中，大氣溫度高于土壤溫度，雨滴進(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ì)受環(huán)境溫度顯著影響所致。

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