馬稚桐，王文科，張?jiān)谟?，趙 明，陳 立，宋 浩，虞佩媛

(1.長安大學(xué) 環(huán)境科學(xué)與工程學(xué)院，陜西 西安 710054；2.長安大學(xué)旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054)

1 研究背景

降雨稀少的西北干旱半干旱地區(qū)，包氣帶水汽熱的遷移轉(zhuǎn)化對(duì)蒸發(fā)過程和生態(tài)環(huán)境保護(hù)起到關(guān)鍵作用，并且土壤溫度對(duì)含水率分布影響顯著[1]。因此，考慮溫度影響的水汽運(yùn)移對(duì)干旱半干旱地區(qū)能量傳輸、質(zhì)量遷移和減少蒸發(fā)損失等具有重要的實(shí)踐意義[2-3]。

Richards 方程在早期多孔介質(zhì)中僅考慮液態(tài)水流的傳輸過程，缺少水熱耦合運(yùn)移過程，后期考慮了液態(tài)水和熱的耦合運(yùn)移，但水蒸汽在很大程度上影響著水熱運(yùn)移。當(dāng)土壤水蒸發(fā)或凝結(jié)時(shí)，水蒸汽流傳輸攜帶能量[4-5]，并影響土壤含水率的變化，尤其在干燥條件下，Liu Changming[6]、Sun Shufeng等[7]發(fā)現(xiàn)非飽和水蒸汽在運(yùn)動(dòng)過程中對(duì)水通量的貢獻(xiàn)量遠(yuǎn)大于液態(tài)水；曾亦鍵等[8]通過實(shí)驗(yàn)研究發(fā)現(xiàn)淺層包氣帶含水率晝夜變化明顯，但沒有考慮造成該現(xiàn)象的土壤內(nèi)部復(fù)雜的水汽運(yùn)動(dòng)過程。Philip等[9]提出關(guān)于液態(tài)水和汽態(tài)水通量的數(shù)學(xué)表達(dá)式，明確水通量由4部分組成,包括等溫液態(tài)水通量、等溫汽態(tài)水通量、非等溫液態(tài)水通量和非等溫汽態(tài)水通量。經(jīng)過多年發(fā)展，當(dāng)前版本的Philip &de Vries模型(PDV)得到明顯改善[10]，雖然PDV模型在實(shí)際應(yīng)用過程中還存在諸多問題，但仍被許多學(xué)者作為研究水汽運(yùn)動(dòng)的重要基礎(chǔ)[11]。以上研究表明水蒸汽對(duì)水熱傳輸過程的重要性，F(xiàn)etzer等[12]再一次重申忽略蒸汽傳輸將對(duì)近地表土壤含水率和溫度的分布產(chǎn)生重要影響。

目前，研究包氣帶水汽熱運(yùn)動(dòng)的方法主要分為直接方法和間接方法。其中，直接方法指通過儀器監(jiān)測(cè)土壤含水率和溫度變化[13-15],間接方法是通過微氣象學(xué)的方法來間接解決熱傳輸和水汽運(yùn)動(dòng)問題[16-17]。Wang Wenke等[18]、Banimahd等[19]認(rèn)為連續(xù)的含水率溫度監(jiān)測(cè)結(jié)合數(shù)值模擬對(duì)認(rèn)知復(fù)雜的水文過程具有十分重要的意義。

雖然國內(nèi)外學(xué)者就以上方法已經(jīng)對(duì)含水率變化、熱傳輸和水汽運(yùn)動(dòng)等進(jìn)行了一系列的研究，但對(duì)含水率晝夜變化的內(nèi)在原因還鮮有研究。本文采用原位實(shí)驗(yàn)結(jié)合數(shù)值模擬的方法，基于實(shí)測(cè)含水率、溫度以及地下水位變化，量化了土壤水蒸汽通量的分布規(guī)律，從而揭示其對(duì)土壤含水率晝夜變化過程的影響機(jī)制。

2 數(shù)據(jù)來源和數(shù)值模型

2.1 試驗(yàn)區(qū)概括

原位試驗(yàn)場(chǎng)位于內(nèi)蒙古鄂爾多斯市烏審旗河南鄉(xiāng)國家氣象站內(nèi)，屬典型干旱半干旱氣候。分析了氣象站多年監(jiān)測(cè)的氣象要素，包括降雨量、風(fēng)速、相對(duì)濕度、凈輻射等。試驗(yàn)區(qū)多年平均氣溫為8 ℃，最小、最大月平均氣溫分別是1月的-8.6 ℃和7月的23.9 ℃，多年平均降雨量和蒸發(fā)量為320 mm和2266 mm，蒸發(fā)強(qiáng)烈，降雨主要集中在每年的7-9月，占全年降雨量的75%左右，實(shí)驗(yàn)監(jiān)測(cè)時(shí)間為2017年8月1日至8月15日，圖1為研究區(qū)監(jiān)測(cè)期內(nèi)降雨量、氣溫和地下水位變化圖。

圖1 研究區(qū)試驗(yàn)期間降雨量、氣溫和地下水位波動(dòng)

2.2 數(shù)據(jù)來源

試驗(yàn)設(shè)置高度190 cm、直徑60 cm的土柱，中間設(shè)有地下水位監(jiān)測(cè)孔，利用DI501(Diver Inc.,±0.05%)儀器對(duì)地下水位進(jìn)行監(jiān)測(cè)；土柱剖面5、10、50 cm深度處分別布設(shè)土壤體積含水率和土壤溫度傳感器(ECH20-5TM，±0.5%)，監(jiān)測(cè)變量通過數(shù)據(jù)采集器(CR-3000，Campbell)每隔5 min自動(dòng)記錄。原位試驗(yàn)場(chǎng)土柱剖面圖見圖2，土柱填充為毛烏素沙地的風(fēng)積沙，其物理特性見表1。

表1 土柱填充風(fēng)積沙物理特性

注：θs為飽和含水率；θr為殘余含水率；Ks為飽和滲透系數(shù)。

2.3 水分運(yùn)移模型建立

HYDRUS-1 D 軟件是美國鹽度實(shí)驗(yàn)室開發(fā)的，用于模擬一維變飽和介質(zhì)中液態(tài)水、熱和溶質(zhì)的運(yùn)動(dòng)，它是基于有限元法，原始水流模型被擴(kuò)展到額外模擬一維水蒸汽流和熱量以及表面能量和水平衡的計(jì)算模型。

2.3.1 基本方程 垂直一維入滲實(shí)驗(yàn)中的水汽運(yùn)動(dòng)屬于一維水流運(yùn)動(dòng)問題，擴(kuò)展后的HYDRUS-1 D模型同時(shí)考慮了水、汽、熱耦合運(yùn)移，其控制方程為：

(1)

式中：θ為體積含水率，cm3/cm3;ql和qv分別為液態(tài)水和汽態(tài)水的水流通量，cm/h；z為空間坐標(biāo)，以向上為正,cm。

其中液態(tài)水流方程為:

(2)

式中：Klh(cm2/s)和KlT(cm2/(s·K))分別為等溫液態(tài)水力滲透系數(shù)和非等溫液態(tài)水滲透系數(shù)；h為負(fù)壓水頭，cm；T為溫度，℃。

水蒸汽流定義為：

(3)

式中：Kvh(cm/s)和KvT(cm2/(s·K))分別為等溫水蒸汽滲透系數(shù)和非等溫水蒸汽滲透系數(shù)；h為負(fù)壓水頭，cm；T為溫度，℃。

2.3.2 土壤水力參數(shù) 基于Van Genuchten和Mualem方程，土壤水分特征曲線θ(h)和非飽和導(dǎo)水率K(θ)表示如下：

(4)

(5)

式中：Ks為土壤飽和滲透系數(shù),cm/s;Se為有效飽和度；m、a和n為經(jīng)驗(yàn)形狀參數(shù)；θs和θr分別為體積飽和殘余含水率,cm3/cm3;l為孔隙連通性參數(shù)，試驗(yàn)介質(zhì)假設(shè)為均勻介質(zhì)。

土壤水分特征曲線是模擬包氣帶土壤水分運(yùn)移的基礎(chǔ)，包含了反映土壤水力特征的6個(gè)參數(shù)(θr，θs，a，n，l，Ks)；試驗(yàn)采集原狀土樣獲得實(shí)測(cè)含水率和基質(zhì)勢(shì)，通過MATLAB軟件包(curve fitting工具箱)對(duì)實(shí)測(cè)值進(jìn)行擬合，得到6個(gè)模型參數(shù)：θr=0.01 cm3/cm3，θs=0.31 cm3/cm3，a=0.045，n=1.6，Ks=8.1 cm/h，l=0.5。

2.3.3 初始條件和邊界條件 本文基于Hydrus-1 D軟件建立一維水汽熱傳輸模型，并通過實(shí)測(cè)土壤含水率和溫度進(jìn)行識(shí)別驗(yàn)證。模型采用等間隔方法，模型域剖面長為190 cm，將剖面按1 cm等間距剖分為190個(gè)網(wǎng)格(圖2)，模擬時(shí)長為360 h。初始條件設(shè)為8月1日實(shí)測(cè)剖面土壤含水率和溫度。

土柱的上邊界由微氣象數(shù)據(jù)和地表能量平衡方程組成。對(duì)于水分傳輸，上邊界為開放邊界，接受降水補(bǔ)給和土壤蒸發(fā)，因此上邊界設(shè)置為大氣邊界條件；由于地下水位常年波動(dòng)，將下邊界設(shè)置為變壓力水頭邊界。對(duì)于熱傳輸，上邊界條件定義為第一類邊界類型(Dirichlet type)，其中0表示土壤表層的位置，下邊界條件定義為第二類邊界類型(Neumann type)。

圖2 原位試驗(yàn)場(chǎng)土柱剖面示意圖(單位：cm)

2.3.4 誤差分析 在模擬過程中，為了評(píng)估模擬結(jié)果的準(zhǔn)確性，通過均方根誤差RMSE和相關(guān)系數(shù)R2來對(duì)比實(shí)測(cè)值與計(jì)算值之間的擬合程度。

(6)

(7)

式中:xi和yi分別為第i個(gè)實(shí)測(cè)值和估計(jì)值；RMSE為均方根誤差，反映模擬數(shù)據(jù)的精確程度；R2為相關(guān)系數(shù)。

3 結(jié)果與討論

3.1 模擬結(jié)果

圖3和4分別為地下5、10 及50 cm 處含水率和溫度的實(shí)測(cè)值與模擬值的對(duì)比圖，模擬值與實(shí)測(cè)值吻合較好。不同深度處土壤含水率的RMSE(土壤埋深5，10，50 cm)依次為0.0134、0.0135和0.0076，R2依次為0.976，0.977和0.947。分析結(jié)果顯示，越接近地下水位，土壤含水率越高，土壤淺層5 cm處的含水率日變化特征明顯且變幅較深層更加劇烈。當(dāng)降雨發(fā)生時(shí)(8月12日)，淺層土壤含水量的波動(dòng)明顯比埋深為50 cm處的大，深層含水率對(duì)降雨響應(yīng)有滯后效應(yīng)，這與Yeh 等[20]提出的觀點(diǎn)相似，歸其原因?yàn)橥寥缹訉?duì)包氣帶水分運(yùn)動(dòng)存在濾波作用。

圖4顯示土壤溫度在相應(yīng)位置處呈明顯的日尺度正弦分布特征。表層土壤受大氣太陽輻射和空氣溫度的直接影響，淺層5 cm處土壤溫度的變化范圍為16～40℃，地下10 cm處的變化范圍為17～35℃，而深層50 cm處溫度變幅較小，溫度變幅隨深度的增加而減小。由于淺層土壤包氣帶受外界等復(fù)雜因素的影響，所以越接近地表，模擬的溫度誤差越大(RMSE<2.243℃，R2>0.937)，但Timlin等[21]認(rèn)為模擬溫度的可允許誤差在2℃～3℃之內(nèi)，所以溫度模擬值在可接受范圍內(nèi)。

3.2 土壤含水率的晝夜分布

為分析影響淺層含水率晝夜變化的內(nèi)在原因，通過原位試驗(yàn)與數(shù)值模擬相結(jié)合的方法，本文選取模擬期蒸發(fā)較強(qiáng)的時(shí)段(8月10日-11日)進(jìn)行研究。圖5為8月10日土壤近地表處含水率模擬值隨深度變化的動(dòng)態(tài)分布，由圖5可以發(fā)現(xiàn)，近地表土壤含水率在白天存在潤濕現(xiàn)象，Goss等[22]、Zhang等[23]同樣發(fā)現(xiàn)當(dāng)大氣溫度較高時(shí)，近地表含水率變大。選取8月10-11日近地表5和10 cm處含水率、凈輻射和氣溫的日動(dòng)態(tài)分布,見圖6，圖6中顯示凈輻射在6:00-18:00之間為正值，其他時(shí)段轉(zhuǎn)為負(fù)值，其最大值出現(xiàn)在正午12:00左右；大氣溫度最大值和最小值分別出現(xiàn)在下午15:00-16:00與早晨5:00-6:00；地表5 cm處含水率的變化范圍在0.085～0.10 m3/m3，其日最大值出現(xiàn)在中午10:00-12:00，日最小值出現(xiàn)在下午18:00-20:00；地下10 cm處含水率的變化范圍在0.124～0.132 m3/m3，其最大值和日最小值分別出現(xiàn)在下午15:00-17:00之間和早晨6:00-8:00之間；含水率極值隨深度的增加發(fā)生滯后，相比5 cm處，地下10 cm的含水率最大值滯后5 h，最小值滯后12 h。圖7為土壤表層含水率與大氣溫度的相關(guān)關(guān)系，當(dāng)大氣溫度高于25℃左右時(shí)，表層土壤含水率與大氣溫度呈明顯的負(fù)相關(guān)關(guān)系，且含水率呈快速減小的趨勢(shì)，這與Wang Wenke等[18]、Chen Li等[24]提出的觀點(diǎn)相似，即25℃是影響土壤水分重分布和水汽通量的一個(gè)閾值，所以識(shí)別土壤溫度的分布特征能更好地認(rèn)知土壤水的日變化規(guī)律。

3.3 土壤剖面溫度的變化特征

太陽輻射作為地球能量的主要來源，不斷對(duì)地氣系統(tǒng)進(jìn)行加熱，土面熱狀態(tài)的改變體現(xiàn)在土壤剖面溫度的變化，本次實(shí)驗(yàn)使用凈輻射表示太陽輻射的變化(圖6)。在白天，隨著光照強(qiáng)度的增大，太陽輻射到達(dá)地面后很快被吸收，地表吸收的熱量通過多孔介質(zhì)孔隙向下傳導(dǎo)。圖8為日尺度模擬的土壤剖面溫度，由圖8可看出，10 cm以上土壤溫度梯度絕對(duì)值遠(yuǎn)大于10 cm以下，地下70 cm處的溫度基本保持不變。綜合分析圖6和8可以得出：在0:00-6:00之間，凈輻射為負(fù)，土壤表面溫度不斷下降，淺層(0～20 cm)溫度梯度為正值；在6:00-12:00之間，隨著凈輻射的逐漸增大，地面吸熱導(dǎo)致溫度逐漸增大，淺層土壤剖面溫度梯度變?yōu)樨?fù)值；12:00-18:00之間，凈輻射逐漸減小，但土壤表面吸收的熱量仍多于放出的熱量，因此到下午15:00左右，地面溫度才出現(xiàn)最大值，之后隨凈輻射的變小而逐漸減小；在18:00-24:00之間,地表放熱多于吸熱，地表溫度變小，淺層溫度梯度變?yōu)檎怠ssouline等[25]發(fā)現(xiàn)忽略溫度梯度后，隨深度變化的土壤含水率動(dòng)態(tài)分布特征不會(huì)呈現(xiàn)日尺度的潤濕現(xiàn)象，F(xiàn)etzer等[12]通過水汽熱耦合研究發(fā)現(xiàn)土壤溫度變化與水蒸汽運(yùn)動(dòng)密切相關(guān)。

下面將通過溫度梯度影響下的剖面水蒸汽通量分布規(guī)律對(duì)包氣帶水分分布進(jìn)行探討。

3.4 包氣帶水分布特征討論

在土壤溫度梯度的作用下，水蒸汽流的剖面分布決定著淺層包氣帶的水分分布特征，土壤水蒸汽流與含水率的關(guān)系表示如下：

(8)

(9)

公式(8)表明水蒸汽壓強(qiáng)大小受溫度影響顯著，忽略溫度梯度將對(duì)估算蒸汽通量產(chǎn)生較大誤差，靠近土壤表層尤為明顯[12]，因此飽和蒸汽壓與剖面溫度密切相關(guān)。結(jié)合上節(jié)對(duì)土壤剖面溫度分布的描述可知，10 cm以上的土壤溫度隨時(shí)間呈明顯的晝夜變化，當(dāng)埋深為0～10 cm時(shí)土壤溫度劇烈變化時(shí)，水蒸汽也會(huì)產(chǎn)生明顯的運(yùn)動(dòng)，導(dǎo)致土壤水進(jìn)行重新分布。圖9顯示淺層包氣帶水蒸汽運(yùn)動(dòng)過程在日尺度內(nèi)可分為3個(gè)階段：階段1(1:00-7:00)，階段2(8:00-16:00)及階段3(16:00-23:00)，這與Zeng等[15]研究淺層包氣帶的土壤水的動(dòng)力學(xué)過程所得出的通量變化模式相似。

圖3 試驗(yàn)期土壤剖面不同埋深的含水率實(shí)測(cè)值與模擬值對(duì)比

圖4 試驗(yàn)期土壤剖面不同埋深的溫度實(shí)測(cè)值與模擬值對(duì)比

圖5 模擬晝夜變化含水率沿埋深的動(dòng)態(tài)分布(8月10日)

由圖9可看出，階段1中，整個(gè)剖面均產(chǎn)生向上運(yùn)動(dòng)的水蒸汽，并且水蒸汽通量隨著時(shí)間的推移逐漸收斂，這是因?yàn)闇囟忍荻仍谠摃r(shí)段內(nèi)為正值。隨著日照的出現(xiàn)，土壤表層溫度逐漸變大，溫度梯度趨勢(shì)變反(方向向下)，階段2主要產(chǎn)生向下運(yùn)動(dòng)的水蒸汽，該時(shí)段水蒸汽保持在土壤中不被蒸發(fā)，造成日尺度土壤含水率在該時(shí)段較其他時(shí)段更高。隨著溫度梯度的變化，階段3產(chǎn)生向上運(yùn)動(dòng)的水蒸汽并被逐漸蒸發(fā)，進(jìn)而造成含水率的逐漸減小，水蒸汽在剖面運(yùn)動(dòng)所產(chǎn)生的零通量面隨時(shí)間逐漸下移，然而在17:00左右，地下10 cm附近產(chǎn)生了向下的水蒸汽運(yùn)動(dòng)，所以10 cm以下的土壤水分保持在土壤中不被蒸發(fā)，導(dǎo)致地下10 cm處土壤的最大含水率滯后于5 cm處。

分析可知，淺層包氣帶中，當(dāng)近地表處土壤溫度梯度為負(fù)時(shí)，產(chǎn)生向下的汽態(tài)水通量，造成土壤含水率變大。曾亦鍵等[8]采用原位測(cè)試方法進(jìn)行觀測(cè)，通過熱傳導(dǎo)方程得出相似的含水率晝夜變化規(guī)律，但忽略了造成該現(xiàn)象的內(nèi)在物理機(jī)制。液態(tài)水保持著水力連接性，很多研究者對(duì)溫度動(dòng)力學(xué)影響含水率變化的機(jī)理過程僅做了淺顯的解釋[11]。Assouline等[26]同樣發(fā)現(xiàn)溫度梯度的變化會(huì)造成這種潤濕現(xiàn)象；Shahraeeni等[27]利用孔隙尺度下的干燥模擬，發(fā)現(xiàn)溫度梯度對(duì)水分傳輸產(chǎn)生相似的影響。另一方面，Goss等[22]發(fā)現(xiàn)當(dāng)空氣溫度最高時(shí)，近表層土壤含水率增加，認(rèn)為向上的土壤水通量造成了地下土壤含水率的積累；Zhang等[23]認(rèn)為土壤內(nèi)部汽化帶的存在可能是造成含水率日變化升高的原因之一。本研究考慮了溫度梯度影響下的土壤水蒸汽運(yùn)動(dòng)的變化模式，闡述了水蒸汽運(yùn)動(dòng)對(duì)土壤含水率變化的影響過程，揭示了影響土壤水分晝夜分布的內(nèi)在物理過程。

圖6 土壤含水率、凈輻射和大氣溫度的日變化規(guī)律

上述結(jié)果與分析表明，淺層包氣帶的土壤水分變化特征受剖面溫度和水蒸汽運(yùn)移方向的影響較大，當(dāng)溫度梯度向上時(shí)，產(chǎn)生向上運(yùn)動(dòng)的水蒸汽通量，含水率變小，反之，含水率變大。鄂爾多斯盆地風(fēng)沙灘區(qū)位于西北地區(qū)東部，蒸發(fā)強(qiáng)烈，生態(tài)環(huán)境脆弱，植被蒸騰主要消耗包氣帶水分，而且土壤蒸發(fā)是一種在近地表附近由液態(tài)水轉(zhuǎn)為汽態(tài)水的過程，與包氣帶水的分布規(guī)律存在一定的互饋?zhàn)饔?，因此，研究淺層包氣帶的水分分布特征對(duì)維持該地區(qū)植被生態(tài)環(huán)境保護(hù)和調(diào)節(jié)土壤含水量等具有重要作用。

圖7 土壤表層含水率與空氣溫度的相關(guān)關(guān)系

圖8 剖面土壤溫度模擬值晝夜不同時(shí)刻分布規(guī)律

圖9 淺層包氣帶土壤剖面水蒸汽通量的晝夜不同時(shí)刻分布規(guī)律

4 結(jié) 論

(1)本研究在鄂爾多斯風(fēng)沙灘地區(qū)，采用原位試驗(yàn)與數(shù)值模擬相結(jié)合的方法，分析了土壤水分布及晝夜含水率變化機(jī)理。結(jié)果表明，模型能夠很好地模擬土壤內(nèi)部水蒸汽通量隨時(shí)間的變化趨勢(shì)，模型具有較好的仿真性，可用于分析土壤水晝夜分布規(guī)律。

(2)通過構(gòu)建水汽熱耦合模型得出日尺度剖面溫度的變化規(guī)律，發(fā)現(xiàn)土壤溫度受太陽輻射的直接影響，當(dāng)凈輻射為負(fù)時(shí)，土壤表層放熱導(dǎo)致溫度變??；當(dāng)凈輻射為正且逐漸變大時(shí)，土壤表層不斷吸熱導(dǎo)致地溫逐漸變大，表層土壤溫度的最大值出現(xiàn)在15:00左右。

(3)淺層包氣帶水蒸汽運(yùn)動(dòng)過程在日尺度內(nèi)分為3個(gè)階段，蒸汽通量的大小受土壤溫度的影響。當(dāng)溫度梯度為正時(shí)(階段1 1:00-7:00)，水蒸汽向上運(yùn)動(dòng)，使土壤含水率從極小值逐漸增大；當(dāng)溫度梯度為負(fù)時(shí)(階段2 8：00-16：00)，水蒸汽向下運(yùn)動(dòng)，使土壤含水率從極大值逐漸減??；受溫度梯度影響，階段3(16：00-23：00)初期地下10 cm處產(chǎn)生向下運(yùn)動(dòng)水蒸汽，造成地下10 cm處的極值含水率滯后于3 cm處的極值含水率。

(4)從生態(tài)環(huán)境保護(hù)和水能平衡等要素衡量，包氣帶土壤水分的合理分配可對(duì)進(jìn)一步維持該地區(qū)水資源可持續(xù)利用提供技術(shù)支撐。