川西北高寒區凍融交替作用后土壤水熱運移模擬研究

宋 潔, 李曉寧, 趙 麗, 樊 敏

(1.西南科技大學 環境與資源學院, 四川 綿陽 621010; 2.河北邯鄲市瀚林環境評估有限公司, 河北 邯鄲 056000)

川西北高寒區處于長江和黃河的發源地，能夠在黃河和長江上游的水源涵養與補給、生態平衡中發揮著重要作用[1]。高寒地區土壤長期經歷晝夜及季節性凍融作用的影響，土壤內部結構易于發生變化：孔隙度增大，容重減小，破壞團聚體水穩性，并減小其抗剪強度，增大土壤可蝕性，造成土壤退化[2-4]。同時，土壤周期性凍融作用對動植物生存存在影響：土壤熱量傳輸，水分相變和鹽分堆積，改變了土壤的水分傳輸能力，從而影響化學物質的運移[5]，導致植被退化，野生動物的棲息地質量和生物多樣性下降。因此，土壤水熱運移規律的研究對于綜合評價地表、地下水資源，有效地利用土壤水、熱資源，合理解決高寒區資源的開發、保護生物多樣性和土壤沙化防治與植被恢復等實際問題都具有重要意義。尚倫宇等[6]揭示了土壤水熱變化對青藏高原地表能量的影響。于煒[7]分析了科爾沁地區沙坨地和草甸地兩種典型地貌條件下的土壤凍融期內水熱動態變化過程及聯系。隨著科學技術發展，許多學者在研究土壤水熱運移時將模型模擬方法運用其中，Harlan[8]最早提出了基于非飽和土遷移機制的土壤凍結過程中水熱耦合運移模型，鄭秀清等[9]采用包括水遷移和熱對流遷移的水熱耦合數值模擬模型，模擬天然條件下土壤的季節性凍融過程以及其中的水熱遷移規律。李瑞平等[10]利用SHAW模型研究了內蒙古河套灌區鹽漬化土壤凍融期水熱鹽的動態變化；王宇祥等[11]利用HYDRUS軟件模擬科爾沁沙地沙丘—草甸相間區土壤水分動態變化指出：流動沙丘和草甸地降雨與表層土壤水分呈極顯著相關。目前，國內外關于高寒區凍融土壤研究多在融雪入滲水熱運移規律[12]、凍融土壤物理結構分析[13]、氣候對土壤凍融影響等問題的試驗與模擬研究，且國內研究多在藏北、東北地區進行，對于川西北高寒區凍融土壤水熱運移對土壤沙化進程的影響研究較少。本文以川西北高寒區經過室內凍融交替作用后沙化和天然草地土壤為研究對象，利用HYDRUS軟件模擬土壤剖面水熱變化，明確川西北高寒區土壤水熱運移特性，探明凍融交替條件作用后沙化和天然草地土壤體積含水率和溫度隨土層深度變化的規律，為該地區土壤沙化治理與預防，維護生態系統平衡提供科學依據和理論支持。

1 材料與方法

1.1 研究區概況

研究區位于阿壩藏族羌族自治州紅原縣(31°50′—33°22′N，101°51′—103°23′E)，位于青藏高原東部，地勢西北高、東南低，海拔在3 000～4 000 m范圍內。研究區內季節性凍土分布較廣，氣候屬大陸性高原寒溫帶季風氣候，四季變化劃分不明顯，冬長夏短，寒冷氣候占據全年大部分時間；但日照長，太陽幅射強烈，早晚溫差極大，極端最高氣溫24.6℃，極端最低氣溫達-22.8℃，年平均氣溫2.9℃[4]。阿壩藏族羌族自治州紅原縣干雨季節分明，雨熱同季，降水量較為豐沛，年降水量可達860.8 mm，年均積雪期可達76 d。研究區主要的土壤類型為沼澤土、亞高山草甸土和風沙土。按照土壤質地的國際分類制，研究區沙化草地土壤屬于砂質土；天然草地土壤屬于砂壤質土[4]。沼澤、草甸、灌叢和森林為該區主要植被類型[4]。

1.2 供試材料

1.2.1 土壤采集及物理性質測定 研究采用的土壤取自阿壩藏族羌族自治州沙化最嚴重的地區之一紅原縣瓦切鄉。依據《天然草地退化、沙化鹽漬化的分級指標》(GB19377—2003)選取沙化草地和植被覆蓋良好的天然草地2類。沙化草地的平均蓋度13%，物種豐富度低，極多裸露砂粒，干燥、人為擾動大；天然草地土層較厚，平均蓋度95%以上，植物物種豐富。植被類型主要有四川蒿草(Kobresiasetchwanensis)、細葉亞(Ajaniatenuifolia)、小柴胡群叢(Bupleurumtenue)等。

土壤采集原狀土，每15 cm為一層，共取4層土壤樣品，用于裝填實驗室水熱運移模擬土柱。每層土壤分別用環刀重復取樣3次，用于測定土壤容重和含水率。采用環刀法測定土壤含水率和容重，土壤物理性質見表1。

表1 供試土壤物理性質

1.3 研究方法

1.3.1 試驗裝置 本研究采用3個自制土柱土壤積水入滲裝置進行入滲試驗(圖1)，將原狀采集土壤各層分別風干后過2 mm篩，沙化草地和天然草地土壤各層含水率和容重按照野外測定數據設置(表1)，采用干堆法，根據公式(1)在每次裝入土壤后，壓實使其達到規定的土柱高度15 cm，保證試驗土柱各層土壤干容重與天然土壤干容重相近。

W=V·γ(1+S)

(1)

式中：W為土壤質量(g)；V為每次裝入的土體體積(cm3);γ為天然土體干容重(g/cm3);S為室內土壤含水量(%),室內土壤含水量則根據采樣時的實測值,采用稱重法,保持土壤重量含水率與實測值一致。沙化草地和天然草地土柱裝填高度均為 60 cm,每隔12 cm開設孔洞(12 cm,24 cm,36 cm,48 cm和60 cm),用于插入土壤水分傳感器探頭(型號ECH2O),水分傳感器探頭孔洞對側設置排氣孔(φ=0.2 cm)。土柱底部設置排氣孔板和排水室(h=0.2 cm,φ=10 cm),出水由燒杯直接截取,利用馬氏瓶穩壓供水,供水水溫控制為室溫25℃。

1.3.2 凍融條件設定 根據川西北高寒區實際凍融現狀，設置連續反復凍融(模擬晝夜)模式，依據凍結溫度設置為-11℃，每天凍結12 h，在空調室內恒溫25℃進行融化，消融12 h的凍融條件，連續凍融6 d，第6天融化后進行水熱運移試驗。

圖1 土柱試驗裝置示意圖

1.3.3 試驗步驟 經6 d連續反復凍融作用后，對土柱進行恒定水頭入滲試驗。土柱分別在12 cm，24 cm，36 cm，48 cm和60 cm處插入ECH2O監測探頭，數據采集儀自動記錄剖面土壤含水率和溫度變化實時監測值，讀取頻率1次/ min。控制馬氏瓶內水面高度，確保土柱表面水分入滲恒定水頭2 cm。當排水室有水流溢出時，試驗結束。由于試驗持續時間較短，試驗分析中不計蒸發對入滲過程的影響。

2 HYDRUS模型構建

2.1 基本方程

(2)

式中:θ為土壤體積含水量(cm3/cm3);t為時間(min);z為一維垂直向坐標(cm);K(h)為土壤的非飽和含水率(cm/s);h為壓力水頭(cm);z為土柱縱剖面空間坐標(cm),原點在圖層的上邊界,向下為正。本試驗不考慮作物根系對土壤水的吸收作用,故S為0。

土壤水分特征曲線和非飽和導水率用Van Genuchten方程擬合：

(3)

(4)

(5)

式中:θ為土壤體積含水率；θs為土壤飽和含水率(cm3/cm3)；θr為土壤殘余含水率(cm3/cm3)；K為土壤的非飽和導水率(cm/s)；Ks為飽和導水率(cm/s)；Se為無量綱的相對含水量；α，m(m=1-1/n)，n,l(一般取0.5)均為擬合參數。

2.1.2 熱運移方程 土體熱量運移采用的基本方程(僅考慮液態水運動)為：

C(θ)=Cnβn+Cwβw+Caav

(6)

式中:Cp(θ)為多孔介質比熱容〔J/(cm3·℃)〕;Cn為固相比熱容〔J/(cm3·℃)〕;Ca為氣體比熱容〔J/(cm3·℃)〕;Cw為液體比熱容〔J/(cm3·℃)〕;λ(θ)為土壤導熱率(cm2/s);q為水分通量(g/s)；T為土壤溫度(℃)。

2.2 初始條件與邊界條件

2.2.1 初始條件 在研究土壤入滲過程中的水熱運移規律時，初始條件為各層土壤初始含水率，不考慮蒸發、降雨和地下水對試驗的影響情況下，t=0時，土體剖面水流模型初始條件：

h(z,t)=h0(z)t=0, 0≤z≤Z,Z=66 cm

T(z,t)=t0(z)

(7)

式中:h0(z)為土體初始負壓水頭(cm);t0(z)為土壤初始溫度(℃)。

2.2.2 邊界條件 水流上邊界為恒定2 cm水頭邊界。下邊界為自由排水邊界：

h(z,t)=h0(t)z= 0

(8)

式中:h0(t)為上邊界邊界壓力水頭值(cm)；k為土壤非飽和導水率(cm/min)；q0(t)為土壤下邊界水通量(cm/s)。

T(z,t)=T0(t)z= 0

(9)

式中:T0(t)為進水流溫度(℃);T和T0分別為土壤和下邊界溫度;λ為為土壤導熱率〔J/(cm·K)〕；Cw為液體比熱容〔J/(cm·K)〕；q0為下邊界土壤水通量(cm/s)。

2.2.3 土體剖面信息 模擬土體剖面空間步長設置為1cm，60 cm土柱則剖分為60層，共60個節點，距地表不同距離處設定土壤水分、溫度及水勢變化監測點，含水量和溫度數據獲取點為4個，分別在12 cm，24 cm，36 cm和48 cm處。設置初始時間步長0.001 d，可變步長設置為0.5～0.001 d。

2.3 模型參數

2.3.1 土體水力學參數 在HYDRUS-1 D模型的水分運動模型中，輸入不同土層深度的粒級分析(砂粒、粉砂粒和黏粒的百分含量)和容重值，采用RETC軟件對已有試驗數據進行擬合，并利用軟件獲取Van Genuchten模型中的土壤含水率參數(θr，θs，α，n和l)以及飽和水力傳導系數(Ks)等相關參數，土體相關水力參數見表2。

2.3.2 土體熱力學參數 熱力學參數利用HYDURS軟件的內部程序Rosetta Lite.v.1.1獲得，土壤固相占總體積的比率為0.57；土壤縱向熱擴散率DL(cm2/s)為5.00；土壤橫向熱擴散率Dr(cm2/s)為1.00；熱導率函數系數(b1)為1.13×108；熱導率函數系數(b2)為1.53×108；熱導率函數系數(b3)為1.16×108；土壤固相熱容Cn(J/(g·℃)為6.91×108；土壤有機質熱容Co(J/g·℃)為1.04×108；土壤液相熱容Cw(J/g·℃)為1.54×108。

2.4 模型驗證與評價

2.4.1 模型驗證 對實測的土體數據進行擬合，分析土柱體積含水率和溫度模擬數據的模擬效果。本文采用均方根誤差(RMSE，公式10)及決定系數(R2，公式11)2個指標來評價模型的模擬效果。

(10)

(11)

2.4.2 土壤體積含水率和溫度動態變化的驗證結果為了驗證該模型模擬的效果，利用HYDRUS-1D模擬經過連續凍融作用后沙化草地和天然草地土壤含水率和溫度隨時間變化過程，結合實測數據，繪制土壤體積含水率和溫度模擬值與實測值對比曲線圖(圖2和圖3)。為了觀察模擬數據的模擬效果，對實測的土體數據進行擬合，采用均方根誤差及決定系數2個評價指標對模擬值和實測值之間的偏差進行評價。

表2 率定后土壤水分特征參數

從圖2和圖3可見，沙化草地和天然草地各土層體積含水率和溫度模擬值與測定值接近，且變化趨勢基本一致。天然草地15 cm處土壤體積含水率模擬值略低于測定值，這種差異的存在說明該層土體接近地表，容易受到外界環境因素的影響。60 cm處土壤體積含水率模擬中期數據高于測定數據，說明土體初始含水量越高，土柱在融化時間內未能完全融化，導致實測時體積含水率偏低，隨著深度的增加模擬精度有所提高，這可能與接近地表的熱通量受到外界的影響較多。

圖2 不同土層土壤水分隨時間變化的模擬值與測定值對比

2.4.3 模型模擬結果評價 經過參數的率定和驗證，沙化草地和天然草地土壤體積含水率和溫度在剖面各節點上的相對差值都較小，展現了比較好的擬合效果和擬合精度。模擬值與實測值的均方根誤差RMSE在合理范圍內，決定系數R2大于0.9。HYDRUS-1D模型能夠用于模擬土壤含水率和溫度的變化，具有一定的可靠性與穩定性(表3)。

3 模型應用與分析

3.1 HYDRUS-1D模擬預測土壤水分變化過程

選取2010年實測的輻射、氣溫、濕度等氣象數據輸入模型，模擬并分析土壤經過凍融循環作用后，土壤體積含水率變化特性。

圖3 不同土層土壤溫度隨時間變化的模擬值與測定值對比

表3 模擬結果的評價

圖4中顯示凍融循環作用前后各層土壤體積含水率模擬變化情況。凍融前(圖4A和圖4B)，沙化草地0—15 cm土壤體積含水率最低，45—60 cm土壤體積含水率最高，15—30 cm和30—45 cm土壤體積含水率次之，且差異較小；天然草地各層土壤呈現隨著剖面加深，土壤體積含水率逐漸降低的趨勢；沙化草地表層土壤體積含水率與天然草地表層土壤體積含水率相差約10%。反復凍融循環后(圖4C和圖4D)，兩種土壤含水率隨季節變化均呈現上下波動的趨勢，但天然草地各土層含水率范圍在0.08%～0.25%，總體均高于沙化草地含水率的變化范圍為0.06%～0.16%；天然草地各層含水率有明顯差異，0—15 cm土壤含水率最低，15—30 cm土壤含水率最高；沙化草地各層土壤呈現隨著剖面加深，且各土層含水率差異較小。

3.2 模擬預測凍融循環條件下土壤溫度變化過程

土壤凍融狀況受土壤溫度影響，在溫度梯度作用下，土壤溫度會隨著太陽輻射變化和季節更替而出現晝夜變化和季節波動。

從圖5中可看出凍融循環前后，土壤溫度模擬值呈現先波動上升后波動下降的趨勢，60～300 d兩種不同程度植被覆蓋下各層模擬土壤溫度變化趨勢基本一致，均呈現隨著剖面加深，土壤溫度逐漸降低的趨勢，且沙化草地土壤各層最高溫度平均低于天然草地各層土壤溫度1.8℃。60 d前和300 d后兩種不同程度植被覆蓋下各層土壤溫度呈現與60～300 d相反的變化趨勢，即底層溫度最高，表層溫度最低。沙化草地和天然草地表層土壤溫度最大值均出現在6月份下旬和7月份，氣溫最高分別為12.4℃和15.6℃，土壤溫度最低為1月份下旬。

4 討 論

4.1 凍融作用下土壤水分運移特性模擬

通過對沙化草地和天然草地凍融前后土壤中水分運移模擬表明，凍融方式和植被覆蓋度對水分運移的影響極大。郭志強等[14]研究表明，夏季晝夜溫差大，土壤體積含水率在強烈蒸發、強輻射作用下出現顯著變化，導致土壤晝夜凍融反復循環的發生；冉洪伍等[12]對藏北高寒草地土壤凍融過程水分變化的研究結果顯示，由于積雪覆蓋于土壤表層，土壤處于長期凍結狀態，各類草地土壤溫度結果均表現為無明顯變化；另有研究表明草甸植被覆蓋度的下降可以導致土壤含水量下降[15]。本研究綜合考慮以上因素對土壤水分運移的影響，結合HYDRUS-1D模型分析水分運移規律及其在草地退化及沙化過程中的作用。

HYDRUS模擬預測凍融前后沙化草地和天然草地的土壤水分運移特性中，不同深度土壤的含水量分別在75 d和50 d之后才明顯增大，分析原因在于一維模型忽略了土壤中空氣壓力的變化[16]和水流的滯后作用。沙化草地土壤體積含水量在120 d前后的變化趨勢相反，考慮到因前120 d為該區的1—4月份，氣溫未回升，無水分入滲和降水補給，且蒸發不強烈；120～300 d為該區的5—10月份，氣溫回升，積雪融化，雨季來臨降水量逐漸增加，氣溫逐漸達到一年中的最高溫(15.6℃)，蒸發作用增強，導致表層的土壤體積含水率低于深層的土壤體積含水率。天然草地在100～300 d內土壤體積含水率沿著剖面自上而下呈降低趨勢，分析原因為此時間段積雪融化水分入滲、降水量逐漸增加，紅原縣牧草生長的茂盛，土壤表面植被覆蓋率高，有效減緩了土壤水分的蒸發，加之植物根系發揮了涵養水源的作用。

圖4 凍融循環后各土層土壤體積含水率模擬

圖5 不同凍融循環后各土層土壤溫度模擬

反復凍融循環后，由于高寒區季節變化、降水差異、蒸發強度和土壤理化性質的差異，土壤水分的垂直分布和時間變化具有明顯的規律。3—5月份為凍土消融前期，氣溫回升雪水融化入滲，各層土壤體積含水率逐漸增加，表層土壤中變化更為明顯。此時，沙化草地和天然草地各層土壤體積含水率沿著土壤剖面向下呈現增加趨勢。6—10月份氣溫升高加速深層凍土消融，夏季降水量增加，沙化草地和天然草地土壤體積含水量均呈不同程度的上升趨勢；7月份左右氣溫達峰值，土壤體積含水率在強烈蒸發作用下出現顯著變化，加之研究區海拔高、太陽輻射強，在凍融期間表層土壤和大氣之間時刻發生著水分和能量的流動，導致了凍融期間晝夜凍融循環的發生，土壤體積含水率波動性明顯。同時土壤質地對于凍土水力傳導度的影響更為顯著，植被覆蓋率較高的天然草地土壤質地黏性大，凍結后各土層水力傳導度的變化幅度要顯著小于砂性較重的沙化草地土壤。10月份至次年3月份為土壤凍結期，積雪覆蓋于土壤表層，隔絕大氣與土層的直接接觸，土壤處于凍結狀態，水分無法入滲，沙化草地和天然草地土壤水分模擬結果均表現為無明顯變化。

研究表明，植被覆蓋率與土壤有機質含量、容重和孔隙度均有顯著相關關系，連續反復凍融過程是土體的密實過程，隨著土壤密實度增加，干容重減小，土體的孔隙率變大，在補水條件下，連續反復凍融后，冰變成水體積減小承載力下降，土體內部發生沉降。連續反復凍融下對天然草地表層土壤容重和孔隙較其他土層影響大，0—10 cm土層含水率明顯低于其他土層，主要是因為凍融循環使表層土壤的毛管孔隙度下降，持水力降低。植物根系及其生物量的大小在蓄水保肥、防風固沙方面至關重要，根系生物量及其分布特征與容重、有機質的物化性質變化相關性較高，兩者的變化均易于引起土壤持水性的變化，草甸植被覆蓋度的下降可以導致土壤含水量下降。

4.2凍融作用下土壤溫度變化特性模擬

沙化草地和天然草地土壤表層和深層溫度在時間序列上呈現截然相反的變化趨勢，這可能與接近地表的熱通量受到外界的影響較多，另外說明土壤熱參數的取值有一定的誤差[17]。沙化草地和天然草地土壤溫度在季節上呈現一定的周期性，3月份以前和9月份以后，土壤溫度呈現上層低于下層的趨勢，3—9月份土壤溫度從表層到底層逐漸降低，距表層越近溫度變化幅度越大[4]，加之熱量損耗導致土體內部隨深度增加，波動幅度呈現減小的趨勢。在整個模擬周期內，土壤各層溫度呈現不斷波動的現象，由于凍融土壤受降水、蒸發、太陽輻射影響較大，溫度變化劇烈[17]。土壤凍融循環過程中，溫度隨之上下波動，土壤中的顆粒、水分和氣體的占比及分布隨之發生變化，當土壤內的水發生固液相變化時，水或冰會對周圍的土體產生不同的擠壓力，這一作用力會破壞土壤顆粒之間的膠結作用，引起土壤的大孔隙、中孔隙、小孔隙和毛細孔隙位置及形態發生變化，從而也改變了土壤孔隙形態結構。太陽照射能夠引起土壤的溫度升高，在試驗過程中土壤物理性質的變化沒有引起土壤溫度的明顯變化。

5 結 論

基于土壤水分運動的動力學方程和土壤熱流基本方程，建立水熱運移數學模型，利用HYDRUS-1D軟件對模型進行求解，結合實測數據對模型模擬進行評價，主要結論如下：(1) 入滲水頭為2 cm時，經凍融交替作用后60cm土柱剖面各土層土壤溫度和水分的測定值和模擬值吻合較好，率定后的土壤水熱運移數學模型能夠較好地反映出土壤水熱的變化及分布狀態。(2) 利用率定后的飽和含水率(θs)、殘留含水率(θr)等設定HYDRUS-1 D軟件中的模型運行參數，模擬預測長時間(365 d)的土壤水熱運移動態。此結果有助于了解川西北高寒區凍融交替作用下土壤垂直剖面水分溫度的動態變化特征，對進一步研究該地區土壤結構，沙化治理方法和維護生態平衡具有重要的理論意義和應用價值。