耦合干空氣機制的黃土水熱運移過程模擬
——以山地蘋果園為例*

楊 博，趙西寧，高曉東，潘岱立，霍高鵬，葉苗泰

（1.西北農(nóng)林科技大學水利與建筑工程學院，陜西楊陵 712100；2.西北農(nóng)林科技大學水土保持研究所，陜西楊陵 712100）

黃土高原是全球最大的黃土堆積區(qū)[1]，黃土剖面深厚，可達200 m。該區(qū)大部分屬于半干旱半濕潤區(qū)域，降雨量在400～600 mm 之間，光熱資源豐富，是全球最大的優(yōu)質(zhì)蘋果主產(chǎn)區(qū)，種植面積和產(chǎn)量均超過全球的25%[2-3]，但是由于干旱缺水加之果樹耗水強烈，水分供需矛盾十分突出[4]。干旱缺水是制約該區(qū)農(nóng)業(yè)生產(chǎn)與生態(tài)建設的主要瓶頸之一[5]，其中土壤水熱條件是最根本因素[6-7]。因此，了解土壤水熱動態(tài)過程可為該區(qū)果園生產(chǎn)力預報奠定基礎。

數(shù)值模擬是再現(xiàn)土壤水熱動態(tài)過程的一種可行且較為經(jīng)濟的方法[8-9]。與基于數(shù)據(jù)驅(qū)動的統(tǒng)計模型相比，基于物理過程的土壤水熱耦合模型具有較強的物理意義，經(jīng)得起理論推敲與試驗驗證，可復制性可移植性強[10-12]。目前，雖然已有研究驗證了土壤水熱耦合模型的可行性[13-15]，但是大部分模型假定土壤空氣壓強與大氣壓強始終保持平衡，空氣自由出入地表界面，并認為土壤中只存在水相運移，忽略土壤空氣流動帶來的內(nèi)部壓力變化[16-17]。實際上，土壤水熱運移是水、氣兩相流體在土壤孔隙中相互驅(qū)動的復雜過程，具有二相流特質(zhì)，土壤干空氣顯著影響著土壤水分運移[9]。二相流質(zhì)熱傳輸模型STEMMUS 首次將土壤干空氣方程耦合進土壤水汽熱過程模擬中，綜合考慮了土壤水、汽、熱、氣之間的耦合過程，并已得到初步應用[9，17]。黃土是一種點楞支架式多孔疏松的土壤，氣相所占比例較高，總孔隙度可達30%以上[18]。該區(qū)在7—9 月份集中降雨期經(jīng)常發(fā)生短時強降雨現(xiàn)象，極易形成超滲產(chǎn)流，產(chǎn)生有壓入滲，使得土壤空氣顯著影響水分和熱量運移過程[19-20]。但是，目前考慮土壤干空氣機制的黃土土壤水熱運移過程模擬較為薄弱。

本文將以黃土高原山地果園為對象，在果園土壤水熱監(jiān)測數(shù)據(jù)基礎上，采用STEMMUS 模型模擬果園剖面土壤水熱運移過程，探討土壤干空氣機制對土壤水分運移過程模擬的影響，以為黃土區(qū)土壤水分管理提供科學依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

試驗在陜西省子洲縣清水溝村山地有機蘋果種植基地（37°26′N，110°02′E）進行，海拔高度1 051 m，為典型的黃土丘陵溝壑地形。該區(qū)降水年際變化較大，年均降雨量為450 mm，多集中在7—9 月，年均氣溫為9.1℃，多年平均無霜期為145 d，屬于溫帶大陸性季風氣候。園區(qū)土壤主要為黃綿土，屬于沙壤土，田間持水量約為22%（體積含水率），平均土壤容重為1.42 g·cm–3。

2018 年5 月在試驗區(qū)果園內(nèi)隨機選擇長勢相當?shù)? 株果樹（平均種植密度為2 m×3.5 m）作為3 個重復，供試蘋果樹為紅富士品種，樹齡為10 a（至2019 年），試驗觀測時間為2018 年6 月1 日至2018 年9 月28 日。果樹具體測定生長指標見表1。

表1 果樹生長指標Table1 Apple tree growth index

1.2 觀測項目及方法

土壤含水率及土壤溫度：采用EC-5 型傳感器（Decagon Devices Inc.，USA）測量果樹根區(qū)土壤含水率與土壤溫度。距樹干40 cm 處預挖剖面，將傳感器垂直插入土壤剖面10 cm、20 cm、40 cm、60 cm、80 cm 處。同時，將EC-5 所測土壤水分和烘干法所測體積含水量（體積含水量由質(zhì)量含水量與環(huán)刀所測土壤容重相乘得到）進行擬合，校正EC-5。用RR-1008 數(shù)據(jù)記錄儀（北京雨根科技有限公司，中國）自動記錄傳感器測量的數(shù)據(jù)，監(jiān)測頻率為每10 min 記錄一次。模型輸入采用1 小時內(nèi)土壤水熱數(shù)據(jù)的平均值。

氣象數(shù)據(jù)：在距離果園50 m 的空曠處安置小型氣象站（Decagon Devices Inc.，USA），氣象數(shù)據(jù)監(jiān)測項包括逐日大氣溫度（T）、空氣濕度（RH）、降雨量（P）、太陽輻射（PAR）及距地面2 米處風速（D）。系統(tǒng)每30 min 記錄一次數(shù)據(jù)。

果樹生長指標：采用LAI-2200C 冠層分析儀（Li-Cor，USA）對果樹葉面積指數(shù)（LAI）進行測量，測量間隔為7～10 d 一次。采用LI-6400XT 光合儀（Li-Cor，USA）測定果樹葉片氣孔導度。測定時間為白天10 點～12 點之間，此時果樹氣孔導度達到頂峰并保持穩(wěn)定，測量間隔為1 個月。

果樹根系密度分布：采用根鉆法對生長期末果樹根系進行實際取樣，取樣深度為2 m，將所取樣品過40 目篩網(wǎng)，經(jīng)人工揀根，利用掃描儀掃描根系圖像（300 dpi），Delta-t scan 軟件（Delta-T Devices Company，UK）分析根系圖像，由圖像分析結(jié)果剔除直徑大于2 mm 的粗根，獲取細根根長。將獲得的各土層總細根根長除以對應取樣土體體積得出果樹各土層細根根長密度。

1.3 STEMMUS 模型

在非飽和土壤水、汽、熱耦合運移理論[21-22]基礎上，STEMMUS 模型充分考慮了土壤干空氣對土壤水運動的影響，將土壤空氣壓強視為一個獨立變量，耦合了土壤空氣運動方程[23]，使之能夠成功描述非飽和土壤水、汽、干空氣和熱耦合運移過程[9]，并運用隱式有限差分法，通過Matlab 軟件編寫程序?qū)崿F(xiàn)對控制方程的數(shù)值求解。同時，在該模型基礎上，研發(fā)人員又編入Feddes 根系吸水模型[24]及兩種不同計算方法的蒸散發(fā)模塊[25]，成為能模擬土壤含水率、土壤溫度、植物蒸騰及土壤蒸發(fā)的SPAC模型[17]。

（1）邊界條件。本研究土壤水分運動上邊界選用大氣邊界描述，包括自然降雨與土壤蒸發(fā)。試驗樣地位于水平塬面上，地勢平坦，地下水埋深較大。下邊界條件假定為自由出流邊界。熱傳輸上邊界條件采用實測土壤表層溫度，熱傳輸下邊界條件采用實測底部溫度。土壤空氣壓強梯度設為零。

（2）模型參數(shù)輸入。通過環(huán)刀獲取樣地 0～20 cm 土樣，采用CR22GⅢ型離心機（Hitachi，Japan）測定土壤水分特征曲線，并用 RETC（RETention Curve）軟件擬合VG（van Genuchten）模型的土壤水分特征曲線參數(shù)。用環(huán)刀法測定土壤飽和導水率。模型率定參數(shù)是以實測參數(shù)為基礎，以實現(xiàn)土壤水熱過程模擬效果最優(yōu)為目標，采用試錯法對實測的土壤水分特征參數(shù)θs、θr、α、n 進行參數(shù)調(diào)優(yōu)，不斷優(yōu)化率定期的土壤水熱模擬結(jié)果，確定最優(yōu)參數(shù)。土壤水力特征參數(shù)經(jīng)模型率定驗證后結(jié)果如表2 所示。土壤熱特征參數(shù)參考Yu 等[17]測定值，取值見表2。

表2 試驗地土壤水力特征參數(shù)及土壤熱特征參數(shù)Table2 Soil hydraulic and thermal characteristics of the experimental orchard

為在模型中反應模擬期內(nèi)蘋果樹葉面積指數(shù)（LAI）的連續(xù)動態(tài)變化，將冠層分析儀獲取的實測數(shù)據(jù)擬合成二次函數(shù)（R2=0.79）進行線性描述，擬合結(jié)果如圖1 所示。

模型中Feddes 根系吸水模型的原模型采用標準化根系密度分布函數(shù)對植物根系密度分布的描述，這種函數(shù)只適用于描述一年生作物的根系生長，而果樹作為多年生植物，根系分布不規(guī)則，因此本研究采用相對深度和相對根長密度來描述果樹固定根長密度分布函數(shù)b(x)[26]。

式中， br( X) 為相對根長密度； bmax最大根長密度，試驗所測果樹 bmax均值為0.193 3 cm·cm–3；X 為相對深度； xm為最大根系深度，設為3 m。 br( X) 采用指數(shù)函數(shù)進行擬合（R2=0.76）。

合并式（6）和式（7）得到果樹根長密度分布函數(shù)：

圖1 模擬期內(nèi)果樹葉面積指數(shù)和根細密度分布Fig. 1 Leaf area index and root fine density distribution during the simulation period

（3）空間與時間離散。采用Galerkin 有限元法對控制方程進行空間離散。設定土壤剖面深度為3 m，將剖面設置38 個節(jié)點，劃分37 個單元，5 個觀測點，表層空間較下層空間離散性好。采用隱式后項差分法對控制方程進行時間離散。最大迭代次數(shù)為30，時間步長設為3 600 s，模擬時段為2018 年6 月1 日至2018 年9 月28 日，共2 880 h，120 d。其中，模型率定期從2018 年6 月1 日開始到2018年7 月29 日結(jié)束，模擬時間共1 416 h。驗證期從2018 年7 月30 日開始到2018 年9 月28 日結(jié)束，模擬時間共1 464 h。

（4）模型驗證與評價。采用歸一化均方根誤差（normalized root mean square error，NRMSE）、決定系數(shù)（R2）和一致性指數(shù)（agreement of index，d）評價模型土壤含水率、土壤溫度模擬效果。式中，分別為iP、iO 模擬值和實測值；mO 為實測值均值；i 為次數(shù)。NRMSE 越小，表明模型模擬偏差越小；越接近1，表明模擬結(jié)果吻合度越高。

2 結(jié) 果

2.1 土壤水分的率定與驗證

基于STEMMUS 模型兩種蒸散計算方法，將率定期與驗證期分別模擬的土壤含水率與果樹根區(qū)不同土層深度觀測值進行比較。由圖2 和表3 可以看出，受降雨影響，模擬期內(nèi)0～10 cm 淺層土壤含水率波動頻繁，隨著深度增加，波動幅度降低。在率定期內(nèi)，ETind法模擬0～10 cm 土壤含水率與實測結(jié)果吻合程度較差，雖然在幾次強降雨中，模型對降雨的響應比較敏感，能迅速模擬出水分動態(tài)變化，但模擬值較實測值存在顯著高估，NRMSE 為30.0%，d 為0.81。隨著土層深度的增加，模型模擬的效果變優(yōu)，在土層深度為20～80 cm 時，NRMSE范圍為4.7%～9.5%，d 范圍為0.84～0.96。ETdir法模擬結(jié)果與ETind法模擬結(jié)果近似，除土層深度10～40 cm 內(nèi)模擬結(jié)果較ETind法模擬精度低外，其余土層均優(yōu)于ETind法，NRMSE 范圍為10.0%～14.9%，d 均為0.88。在驗證期內(nèi)，兩種方法模擬結(jié)果與率定期結(jié)果相近，除0～10 cm 土層模擬效果明顯優(yōu)于率定期外，其他土層模擬效果與率定期相差不大，均能較好地模擬土壤水分動態(tài)變化過程。

2.2 根區(qū)土壤貯水量變化

利用土壤貯水量計算方法對模型兩種蒸散計算方法模擬的根區(qū)土壤貯水量動態(tài)與實測值進行對比（圖3）。由圖所示，率定期內(nèi)ETind法模擬效果較ETdir法更接近實測值，NRMSE 為5.6%，d 為0.98。ETdir法模擬值整體較實測值偏低，土壤水分對降雨響應估計不足，NRMSE 為7.2%，d 為0.95。驗證期內(nèi)，ETind法模擬精度較率定期有所降低，NRMSE為5.0%，d 為0.96。ETdir法模擬精度優(yōu)于ETind法，NRMSE 為5.0%，d 為0.96。與實測值相比，兩種方法模擬值整體變化幅度較率定期內(nèi)有所提高，但均存在土壤水分對降雨響應估計不足，以及持續(xù)干旱時期模擬值高于實測值的現(xiàn)象。

圖2 各土層土壤含水率的模擬值與實測值Fig. 2 Simulated values and measured values of soil moisture content relative to soil layer

表3 土壤含水率模擬精度對比Table3 Comparison of accuracy in simulation of soil moisture content

圖3 土壤貯水量的模擬值與實測值Fig. 3 Simulated and measured values of soil water storage

2.3 土壤溫度的率定與驗證

圖4 為STEMMUS 模型兩種蒸騰計算方法在不同土層深度下土壤溫度模擬值與實測值對比結(jié)果。結(jié)合表4 可以發(fā)現(xiàn)，率定期內(nèi)兩種方法模擬的土壤溫度變化一致，但整體高于實測值，尤其在強降雨時期，模型模擬值明顯存在高估現(xiàn)象，NRMSE 范圍為 5.6%～10.0%，d 范圍為 0.76～0.87。而在驗證期內(nèi)，兩種方法模擬值與實測值較為一致，除0～10 cm 土層外，其余土層NRMSE均為0.1%，d 為0.99。

2.4 土壤干空氣對土壤水分動態(tài)模擬的影響

選擇模型驗證期內(nèi)ETdir法模擬考慮干空氣機制的土壤水分動態(tài)變化（圖5）。如圖所示，土壤干空氣對土壤水分動態(tài)模擬產(chǎn)生一定影響，考慮干空氣影響的二相耦合模型模擬值均略低于不考慮空氣影響的單相模型模擬值。尤其在降雨期間，兩種模型模擬值差異較其他時段顯著，且隨著土層深度的增加，兩種模型模擬值之間的差異逐漸減小。結(jié)合表5 可以看出，在10～20 cm 和20～40 cm 土層，耦合模型模擬精度較單相模型模擬精度略低，其余土層耦合模型模擬結(jié)果較單相模型精度更高，NRMSE 范圍為4.6%～13.2%，d 范圍為0.83～0.94。

圖4 各土層土壤溫度的模擬值與實測值Fig. 4 Simulated and measured values of soil temperature relative to soil layer

表4 土壤溫度模擬精度對比Table4 Comparison of accuracy in simulation of soil temperature

圖5 考慮干空氣機制的ETdir 法土壤含水率模擬值與實測值Fig. 5 Measured values of soil moisture content and values simulated with the ETdir method taking soil dry air mechanism into account

表5 考慮干空氣機制的ETdir 法土壤含水率模擬精度對比Table5 Comparison of accuracy in simulation of soil moisture content between the ETdir method taking soil dry air

2.5 單次降雨土壤含水率垂直分布模擬

氣機制的二相耦合模型土壤含水率模擬結(jié)果，圖6b為不考慮干空氣機制的單相模型土壤含水率模擬結(jié)果。選取模擬期內(nèi)單日最大降雨量的雨天作為模擬對象，降雨量為55.4 mm。由此可以看出，表層土壤含水率均在降雨后1 d達到頂峰，在降雨后第2 天降低，并對下層土壤含水率進行補給。對于此次降雨，在降雨2 d 后，耦合模型顯示土層深度0～70 cm內(nèi)土壤含水率變化明顯，說明耦合模型模擬入滲深度達70 cm，而單相模型入滲可影響到80 cm 深土壤含水率，表明單相模型模擬土壤水分入滲速率較耦合模型快。相比單相模型，耦合模型降雨后土壤含水率增加量明顯較小。對于耦合模型，表層土壤含水率由降雨前1 d 的0.124 cm3·cm–3變化為降雨后1 d 的0.221 cm3·cm–3。對于單相模型，表層土壤含水率由降雨前1 d 的0.124 cm3·cm–3變化為降雨后1 d 的0.252 cm3·cm–3。表明降雨后單相模型較耦合模型模擬結(jié)果偏大，土壤入滲速率偏快。結(jié)合表6可以看出，在降雨當日與降雨后1d 內(nèi)，耦合模型模擬精度高于單相模型，說明耦合模型模擬降雨入滲過程更為合理。

圖6 單次降雨前后土壤含水率垂直分布Fig. 6 Vertical distribution of soil water content before and after each rainfall event

表6 單次降雨前后土壤含水率垂直分布模擬精度對比Table6 Comparison of accuracy in simulation of vertical distribution of soil water content before and after each rainfall

3 討 論

3.1 土壤水分的率定與驗證

STEMMUS 模型兩種蒸散計算方法在率定期與驗證期內(nèi)模擬0～10 cm 淺層土壤含水率均較差（圖2 和表3），根據(jù)野外觀察，冠層截留與人為因素可能是產(chǎn)生這種現(xiàn)象的主要原因。本文試驗中土壤水分傳感器安置在果樹冠層下方，果樹冠層對降雨產(chǎn)生截留使?jié)B入土壤中的降雨量較實際觀測值少，而模型在輸入降雨參數(shù)時未考慮到這種情況，所以導致表層土壤水分模擬結(jié)果偏高。表層土壤受除草、踩踏等人為因素影響嚴重，模型難以考慮到相關情況。將STEMMUS 模型率定期與驗證期分別模擬的不同土層水分結(jié)果進行對比，反應出的差異各不相同，沒有明顯優(yōu)劣趨勢，產(chǎn)生這種情況原因之一可能是模型對土壤質(zhì)地均一的假設不合適[14，27]。土壤存在空間異質(zhì)性，且果樹根系生長會影響傳感器安裝位置的土壤性質(zhì)，而本文只采用20 cm 土層處土壤水力學參數(shù)作為模型單一平均的水力特征輸入條件，難以符合實際情況。再者，本文根系吸水模塊選取果樹生長期末固定的根系分布，忽略根系動態(tài)生長過程，這對模型計算果樹根系吸水會產(chǎn)生偏差，從而影響土壤水分模擬效果[28-29]。前人已有研究采用模型模擬黃土高原蘋果園土壤水分，如李冰冰等[29]基于Hydrus-1D 模型對渭北旱塬蘋果園土壤水分進行模擬，RMSE 為0.1 cm3·cm–3，決定系數(shù)R2介于0.74～0.85；童永平等[30]利用Hydrus-1D 模型模擬黃土關鍵帶蘋果園土壤水分動態(tài)變化過程，均方根誤差RMSE 介于0.0149 cm3·cm–3～0.0168 cm3·cm–3，決定系數(shù)R2介于0.59～0.76。本文中STEMMUS模型模擬土壤水分動態(tài)結(jié)果經(jīng)計算，RMSE 介于0.006 cm3·cm–3～0.045 cm3·cm–3，決定系數(shù)R2介于0.81～0.96，評價指標中除部分土層深度RMSE 超出上述研究范圍，其余土層指標均涵蓋或優(yōu)于上述范圍。原因是上述研究中樣本量過少，且均為模擬果園深層水分變化，忽略淺層土壤水分劇烈波動對模擬結(jié)果的影響，引起評價指標RMSE 偏低。整體而言，STEMMUS 模型對模擬黃土丘陵區(qū)果園水分動態(tài)情況具有良好的適用性。

STEMMUS 模型采用兩種蒸散計算方法對不同土層土壤水分進行模擬的結(jié)果不同，根據(jù)以往研究，產(chǎn)生這種情況是因為不同蒸散方法計算得出的蒸散量不同，從而顯著影響SPAC 模型的輸出結(jié)果[14，31-32]。如Anothai 等[31]發(fā)現(xiàn)基于CSMCERES-Maize 中FAO-56 PM 蒸散方法模擬的半干旱條件不同灌溉制度下玉米田間土壤水分效果最優(yōu)；Kingston 等[32]研究發(fā)現(xiàn)不同潛在蒸散發(fā)計算方法會顯著影響區(qū)域水資源的預測結(jié)果；呂厚荃和于貴瑞[33]根據(jù)不同氣候類型的土壤水分數(shù)據(jù)對基于不同ET 計算方法的農(nóng)田土壤水分模擬結(jié)果進行分析，發(fā)現(xiàn)Priestley-Taylor 蒸散模型模擬的土壤水分動態(tài)最優(yōu)。以上研究均表明，不同蒸散計算方法對不同模型模擬土壤水分動態(tài)的影響存在差異，而本文中STEMMUS 模型兩種蒸散計算方法模擬的土壤水分效果相似，均可以較好地滿足模擬需要。

3.2 干空氣機制對土壤水分模擬的影響

在降雨期間，STEMMUS 模型中考慮干空氣機制的二相耦合模型模擬值較單相模型模擬值更接近實測值，且模擬雨水入滲速率偏小，入滲土壤水分增量偏低，影響了模擬得到的入滲深度（圖5 和圖6）。一般而言，模型誤差主要來源于三方面：模型結(jié)構誤差、輸入?yún)?shù)誤差、輸入變量的測量誤差。與原有模型相比，考慮干空氣機制的耦合模型在結(jié)構誤差會降低，但是輸入?yún)?shù)誤差和變量的測量誤差仍然存在，這些誤差的存在造成了模擬值與實測值之間的差異。在實際過程中，短時強降雨導致黃土表面易形成超滲產(chǎn)流，土壤空氣被禁錮，隨著水分不斷入滲，土壤孔隙內(nèi)空氣不斷積壓，形成對雨水入滲的阻礙作用，使土壤含水率上升出現(xiàn)遲滯現(xiàn)象。若不考慮土壤干空氣對土壤水分運移的影響，則會出現(xiàn)土壤含水率上升較快的現(xiàn)象[34]。同時，傳統(tǒng)土壤水分運動方程中液態(tài)水與汽態(tài)水通量只考慮土壤水勢、重力勢與土壤溫度梯度驅(qū)動力的影響[15]，而STEMMUS 模型假設土壤干空氣作為一個獨立狀態(tài)變量存在于土壤中，使得土壤水分通量需要考慮土壤空氣流動，不僅在汽態(tài)水通量中考慮水汽的擴散機制，增加了水汽的對流與彌散機制，還在液態(tài)水通量中加入土壤空氣壓強梯度驅(qū)動[9，17]。因此，在雨水入滲過程中土壤空氣壓強梯度增大，水汽的擴散、對流、彌散運移加強，從而導致土壤空氣流動劇烈，形成對入滲過程的阻礙作用。所以，從研究區(qū)域?qū)嶋H環(huán)境與模型機理兩方面不難看出，考慮土壤干空氣機制模擬土壤水熱動態(tài)過程是很有必要的。

4 結(jié) 論

本文采用考慮土壤水-汽-熱-氣耦合的STEMMUS 模型，以黃土高原山地果園為例，模擬了黃土剖面土壤水熱動態(tài)過程，并用田間同步實測數(shù)據(jù)進行率定與驗證，得出結(jié)論如下：（1）在率定期與驗證期，模型模擬值與實測值吻合較好，STEMMUS 模型可以較好地用于模擬黃土丘陵區(qū)山地蘋果園土壤水分運動、土壤溫度變化情況，為果園土壤水熱動態(tài)實時監(jiān)測與田間水熱環(huán)境管理提供了可靠的工具。（2）STEMMUS 模型中考慮干空氣機制的二相耦合模型模擬值較單相模型模擬值更接近實測值，模擬雨水入滲速率偏小，主要原因是模型考慮土壤空氣流動影響，考慮了土壤干空氣機制，對傳統(tǒng)水熱耦合運移理論進行完善。綜上，STEMMUS 模型無論在理論基礎還是田間數(shù)據(jù)驗證上，都表明可為黃土區(qū)田間土壤水熱實時監(jiān)測與管理提供技術理論支撐。