淺層包氣帶土壤水晝夜分布變化規(guī)律模擬研究

馬稚桐，王文科，張在勇，趙 明，陳 立，宋 浩，虞佩媛

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

1 研究背景

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

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

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

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

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

2.1 試驗區(qū)概括

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

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

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

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

表1 土柱填充風積沙物理特性

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

2.3 水分運移模型建立

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

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

(1)

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

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

(2)

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

水蒸汽流定義為：

(3)

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

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

(4)

(5)

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

土壤水分特征曲線是模擬包氣帶土壤水分運移的基礎，包含了反映土壤水力特征的6個參數(shù)(θr，θs，a，n，l，Ks)；試驗采集原狀土樣獲得實測含水率和基質(zhì)勢，通過MATLAB軟件包(curve fitting工具箱)對實測值進行擬合，得到6個模型參數(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軟件建立一維水汽熱傳輸模型，并通過實測土壤含水率和溫度進行識別驗證。模型采用等間隔方法，模型域剖面長為190 cm，將剖面按1 cm等間距剖分為190個網(wǎng)格(圖2)，模擬時長為360 h。初始條件設為8月1日實測剖面土壤含水率和溫度。

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

圖2 原位試驗場土柱剖面示意圖(單位：cm)

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

(6)

(7)

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

3 結(jié)果與討論

3.1 模擬結(jié)果

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

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

3.2 土壤含水率的晝夜分布

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

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

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

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

3.4 包氣帶水分布特征討論

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

(8)

(9)

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

圖3 試驗期土壤剖面不同埋深的含水率實測值與模擬值對比

圖4 試驗期土壤剖面不同埋深的溫度實測值與模擬值對比

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

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

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

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

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

圖7 土壤表層含水率與空氣溫度的相關關系

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

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

4 結(jié) 論

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

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

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

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