周祖昊,劉揚李,2,李玉慶,王鵬翔,,王 康,李 佳,朱熠明,,劉佳嘉,王富強
(1. 中國水利水電科學研究院流域水循環(huán)模擬與調控國家重點實驗室,北京 100038;2. 中國電建集團貴陽勘測設計研究院有限公司,貴州貴陽 550081;3. 西藏農牧學院水利土木工程學院,西藏林芝 860000;4. 武漢大學水資源與水電工程科學國家重點實驗室,湖北武漢 430072;5. 水利部南水北調規(guī)劃設計管理局,北京 100038;6. 華北水利水電大學水利學院,河南鄭州 450046)
目前,基于寒區(qū)水熱耦合的特點,分布式水文模型已取得一定的進展[1-2]。國際上常用的分布式水文模型都已對寒區(qū)的積雪和凍土過程進行了一定程度的拓展表達,部分模型中已加入凍土、融雪模塊來研究寒區(qū)水文循環(huán),如CRHM(Cold Regions Hydrological model)、SWAT(Soil and Water Assessment Tool)、GEOtop (Geography and Topography model)、DWHC(Distributed Water-Heat Coupled model)、VIC (Variable Infiltration Capacity model)模型等[3-4]。當前較成熟的凍土水熱耦合模型是基于非飽和土體水分遷移,用土壤未凍水含量隨溫度變化的函數來表征凍結機理的水動力學模型。其中以Harlan[5]提出的水熱耦合模型應用最為廣泛;Taylor和Luthin[6]在其基礎上用未凍水含量結合土壤水分特征曲線計算土水勢,確定水分遷移的驅動力;Horiguchi[7]又加入了溫度梯度共同作用于水分遷移;美國農業(yè)部Flerchinger和Saxton[8]綜合考慮氣相、水流遷移引起的熱量變化以及邊界處的蒸發(fā)、輻射等,建立了垂直一維水熱耦合模型,較全面地反映了影響凍土水熱變化因素;雷志棟等[9]對各向同性土柱的一維凍結問題簡化了模擬過程,忽略了水汽遷移、熱量對流作用,提出了凍土非飽和土壤水分遷移方程,對地下水淺埋條件下的土壤凍融過程進行了模擬,模型模擬效果較好。水動力學模型參數簡便易得,是較為成熟的水熱耦合模型。
青藏高原是長江、黃河、瀾滄江、雅魯藏布江等多個大江大河的發(fā)源地。作為世界面積最大、海拔最高的高海拔寒區(qū),平均海拔4 376 m[10],其廣泛分布的積雪和凍土影響了該地區(qū)土壤介質的儲水、導水、產水和熱傳導過程[11]。與中國其他地區(qū)相比,青藏高原土壤層較薄,土壤層下存在較厚的砂礫石層。與土壤相比,砂礫石具有不同的水、熱性質[12],孔隙度、密度也有所不同,對土壤含水量和導水率有較大影響[13-14]。對青藏高原的水熱運移模擬,部分陸面模式雖然考慮了砂礫石影響,通過砂礫石占比對土壤熱傳導率、土壤質地和地形等敏感參數進行調整,模擬效果會有一定提升[15-17],但由于模型將土壤概化為一維均質結構,僅通過參數調整難以體現(xiàn)該地區(qū)這種上下分層明顯的地質結構對水熱運移的影響,模擬結果均有一定偏差。因此,在研究青藏高原土壤凍融與水循環(huán)過程時,不僅要考慮土壤層內的水熱遷移,對于土壤層下伏較厚的砂礫石層也應一并考慮。與此同時,積雪較好的隔熱性和對太陽短波的反射性,也會影響土壤水熱遷移。在水熱耦合研究過程中,要把土壤和積雪作為一個復雜的、相互作用的聯(lián)合體來統(tǒng)一考慮[18]。對積雪的研究,根據其模擬的復雜程度與否,大致可分為以下兩類:一是利用強迫-恢復(force-restore)法[19-20],通過相對簡單的單層積雪模型,區(qū)分積雪和土壤的熱力學性質,模擬積雪-土壤復合層的溫度變化,以度日因子法計算融雪量[21];二是以物理機制為基礎的復雜模型,此類模型考慮了積雪內部的質量守衡和能量平衡以及積雪表面與大氣環(huán)境條件的相互作用,將積雪細分為多層,對積雪內部相變、水分運移等進行了精細描述[22-23]。但后一類模型太過復雜,影響了其推廣與應用。
本文針對青藏高原的氣候和地質特點,在土壤凍結融化期開展了野外土壤和砂礫石層水熱監(jiān)測試驗,以WEP-COR模型[1]的土壤水熱耦合模塊為基礎,構建基于“積雪-土壤-砂礫石層”連續(xù)體的水熱耦合模型,并采用試驗觀測結果對模型進行驗證。
選取雅魯藏布江的一級支流尼洋河流域作為本文研究的典型區(qū)域。尼洋河位于西藏東南林芝地區(qū)雅魯藏布江中下游左岸,北緯29°28′—30°38′,東經92°10′—94°35′,源于西藏自治區(qū)米拉山西側的錯木梁拉,河源海拔約5 000 m,在林芝市巴宜區(qū)匯入雅魯藏布江,河口海拔約2 920 m,落差2 080 m。尼洋河干流長309 km,流域平均坡降為0.73%;流域面積1.75萬km2,在雅魯藏布江眾支流中排行第四。尼洋河流域地處中緯度地帶,年均降水量1 295.1 mm,年均氣溫7.6 ℃左右[24]。
1.2.1 土壤厚度和成分調查與取樣分析
根據尼洋河流域實地考察(圖1)和部分工程設計的地質勘測資料,該流域具有土壤層厚度較薄、土壤下方的砂礫石層厚度較大的特點。尼洋河流域的砂礫石層主要包括礫石層和卵石層,礫石層主要為圓狀礫石,含卵石顆粒,礫石質量百分比約50%~65%,黏粒質量百分比5%~10%,孔隙間填充中細砂;卵石層孔隙間填充圓礫和中細砂,局部含大漂石、孤石。
針對尼洋河流域這一地質特點,為研究流域內土壤層從山腳到山頂的厚度變化規(guī)律,從尼洋河流域的源頭到河口選取了海拔不同的24個取樣點。其中點1—點16沿河道,點17—點24從山腳到山頂測量的土壤厚度(圖2)。對24個取樣點的土壤厚度及組分進行了測量和分析,尼洋河流域土壤主要為砂質壤土和壤質砂土,平均砂粒質量百分比為55.89%,粉砂粒質量百分比為31.2%,黏粒質量百分比為12.91%,厚度隨著海拔升高而減小,平均厚度為53.7 cm。
圖1 土壤-砂礫石層結構Fig.1Soil-sand gravel structure
圖2 尼洋河流域土壤取樣點、試驗點及氣象站位置Fig.2Soil sampling points and meteorological stations
1.2.2 凍結融化期“積雪-土壤-砂礫石層”監(jiān)測試驗
本次研究在尼洋河流域下游色季拉山的山腰上選擇一處開展季節(jié)凍結融化期水熱耦合過程監(jiān)測試驗。試驗點(94°21′45″E,29°27′12″N)海拔為4 607 m,試驗期為2016年11月—2017年4月,試驗期間試驗點降水8.3 mm,平均氣溫-3.4 ℃,最高6.3 ℃,最低-11.7 ℃。凍結融化期,試驗點位置見圖2。試驗人員在試驗點開挖了一個深度為1.6 m的試驗坑,在地面以下垂直深度上每隔10 cm安裝自動測定液態(tài)含水率的TDR傳感器、測量溫度的PT100傳感器和測量基質勢的TensionMark傳感器,對下墊面凍結和融化過程中水、熱、勢能過程進行自動監(jiān)測,并根據PT100傳感器和TensionMark傳感器測量結果推算土壤凍結深度,同時于開始降雪后間隔1~2周根據降雪情況通過米尺實地測量下墊面覆蓋的雪層厚度。試驗前,采用核磁共振方法對高原季節(jié)性凍融條件下水、熱過程監(jiān)測的儀器設備進行率定。儀器完成安裝后進行原狀土回填,保證回填土壤與實際土壤的參數一致。
氣象條件是水熱耦合模擬的邊界條件。由于流域內氣象站點稀少,且高程差別較大,試驗點的逐日氣溫采用考慮高程修正的RDS(Reversed Distance Squared)方法由氣象站數據插值得到,降水數據采用平面插值、結合西藏自治區(qū)水資源公報[25]中降水等值線修正得到,逐日相對濕度、日照時速和風速采用RDS方法插值得到。本文用到的基本氣象數據主要為尼洋河流域境內的林芝站和周邊嘉黎站(圖2)的降水、氣溫、相對濕度、日照時數和風速。氣象數據來源于中國氣象數據網(http:∥data.cma.cn),使用數據時間序列為2016—2017年,時間尺度為日。
本文在WEP-COR模型中凍土水熱耦合計算原理基礎上,結合青藏高原特點加以改進,構建包含“積雪-土壤-砂礫石層”連續(xù)體的高寒山區(qū)水熱耦合模型。
WEP-COR模型在保留原WEP-L模型關于地表與大氣間的能量交換計算的基礎上,系統(tǒng)添加了土壤層間的水熱通量和相關水熱參數的計算過程,并且將活動土壤層細分為11層,見圖3。其中,Rs為地表產流;Ri為第i層土壤的橫向流或壤中流,Ri與坡度和土壤含水率有關;E1為土壤蒸發(fā)量;Er為植被蒸騰量;Qi為第i層土壤的重力排水;P為降水;Ta為大氣溫度;Ti為第i層土壤的溫度;Gi為第i層土壤與相鄰土壤層間由溫度差異引起的熱通量。水熱耦合模擬采用顯式差分進行數值迭代計算,模型模擬的時間步長為1 d。
圖3 WEP-COR模型分層水熱通量計算結構示意Fig.3WEP-COR model water and heat flux simulation structure
2.1.1 土壤水熱運移方程
模型假設土壤凍融時只有液態(tài)水發(fā)生運移,土壤水分的運移主要受重力勢、基質勢和溫度勢的影響,因此添加的土壤一維垂直水分流方程可以寫成[26]:
(1)
根據能量平衡原理,凍融系統(tǒng)中每一層土壤的能量變化都用于系統(tǒng)內的土壤溫變和水分相變,溫度勢是水分相變的驅動力,而大氣與表層土壤的溫差則是熱傳導的源動力。假設土壤各向均質同性,并忽略土壤中的水汽遷移,添加的一維垂直流動的熱流運動基本方程可變形為[27]
(2)
式中:θl、θi分別為土壤中液態(tài)水、冰的體積百分比,cm3/cm3;T為土壤溫度,℃;t、z分別為時間、空間坐標(垂直向下為正);D(θl)、K(θl)分別為非飽和土壤水分擴散率、導水率,cm2/s;Cv為土壤體積熱容量,J/(m3·℃);λ為土壤導熱系數,W/(m·℃);ρi、ρl分別為冰、水密度,kg/m3;Lf為融化潛熱,Lf=3.35×106J/kg。
2.1.2 土壤水熱聯(lián)系方程
土壤凍融過程中,凍土水熱運動間的聯(lián)系主要表現(xiàn)在未凍水的含水率與土壤負溫的動態(tài)平衡中[19-20]:
θl=θm(T)
(3)
式中:θm(T)為土壤負溫對應的最大未凍水含水率。
2.2.1 模型結構
積雪作為一種良好的隔熱體,具有熱容量大、導熱系數較小等特點,對大氣和土壤間的熱傳遞具有明顯的削弱作用。試驗點所處位置積雪持續(xù)時間較長、厚度穩(wěn)定,對土壤水熱的影響不容忽視。根據尼洋河流域實地考察(圖1)和基于26個點的測量結果建立的土壤厚度—高程關系圖(圖4),發(fā)現(xiàn)尼洋河流域土壤厚度整體較薄,在土壤層下方存在厚度較大、混合砂礫石和少量土壤的砂礫石層。而且土壤層厚度在河谷區(qū)基本穩(wěn)定,在山坡上隨著高程增加呈現(xiàn)變薄的趨勢,在山頂附近厚度又基本保持不變。在試驗點處土壤層為砂質壤土,砂粒質量百分比為61.92%,粉砂粒質量百分比為26.2%,黏粒質量百分比為11.88%,土壤層厚度為40 cm;土壤下砂礫石層礫石質量百分比為53.7%,孔隙間填充中細砂。土壤與砂礫石層水分特征參數見表1。
圖4 尼洋河流域土壤厚度—海拔關系Fig.4Relationship between soil thickness and elevation in Niyang River basin
表1 尼洋河流域土壤水動力參數
根據尼洋河流域氣候和地質特征,本文將水熱耦合模擬的對象定義為“積雪-土壤-砂礫石層”連續(xù)體(圖5)。模型計算結構上,在WEP-COR模型11層土壤層分層的基礎上,于土壤層上方加入一層積雪層,同時將原有的11層土壤層細分為土壤和砂礫石層2類介質,計算結構如圖6所示。第0層為積雪,根據等高帶在計算單元中所處位置設置上部i層為土壤層,下部12-i層為砂礫石層。由于距離地表越近,水、熱變化越快,故第1、2層厚度設為10 cm,第3~11層厚度設為20 cm。值得注意的是,土壤-砂礫石層各層的厚度需根據所在位置含水層的實際厚度進行修正,也就是說,當含水層總厚度小于11層土壤-砂礫石層厚度時,須根據含水層實際厚度確定土壤-砂礫石層的計算層數和最后一層的厚度。
圖5 “積雪-土壤-砂礫石層”垂向分層概化Fig.5Vertical layering of “snow-soil-sand gravel layer”
圖6 分層水熱通量計算結構示意Fig.6Layered calculation structure of water and heat flux
2.2.2 計算原理
采用度日因子法[28]計算積雪融雪過程,具體如下:
M=k(Ta-T0)
(4)
式中:M為當日融雪當量,mm;k為積雪融化的度日因子,mm/(℃·d);T0為積雪融化臨界溫度,℃,本文取-1 ℃。
對于積雪的度日因子,一般情況下在1~7 mm/(℃·d),本次研究根據下墊面情況取為4 mm/(℃·d)。由于青藏高原地區(qū)晝夜溫差較其他地區(qū)偏大,本次研究的雨雪臨界氣溫設為2 ℃。
積雪熱量平衡公式[29]:
G=Le+Rn+Q1+Q2
(5)
式中:G為積雪層的能量變化,J/(m2·d);Le為積雪的蒸發(fā)、融化和升華潛熱通量,J/(m2·d);Rn為積雪層表面的凈輻射量,J/(m2·d),由到達積雪表面的凈長波輻射與凈短波輻射計算而得;Q1為積雪與大氣間的熱量交換量,J/(m2·d),采用強迫-恢復法計算;Q2為積雪與土壤間的熱量交換量,J/(m2·d)。
積雪和土壤間的熱量交換量采用以下公式[30-31]:
(6)
式中:ZS為積雪厚度,m;ZC為第一層土壤厚度,m;λS為冰雪層導熱系數,W/(m·℃);λC為土壤層導熱系數,W/(m·℃);RC為土壤層與積雪層的接觸熱阻,(m2·℃)/W;TC為第一層土壤層溫度,℃;TS為積雪層溫度,℃。
土壤層之間的水熱運移方程和聯(lián)系方程采用式(1)和式(2),砂礫石層之間以及砂礫石層和土壤層之間的水熱通量計算采用與土壤層相同的計算原理,但砂礫石層的各項水熱參數與土壤不同。
模型的數值模擬仍采用顯式差分進行數值迭代計算,模擬步長為1 d。每個時間步長內模型先根據初始條件從積雪層到底層砂礫石層逐層計算熱傳導、溫度和水分相變,然后以該層溫度為判定條件進行迭代計算直至收斂(步驟1);熱量計算閉合收斂后,再根據水量平衡計算各層水分運移量,并修正各土壤和砂礫石層(積雪層不考慮)的含水率和含冰率(步驟2);用液態(tài)含水率判定是否收斂,若不收斂則返回步驟1進行熱量計算,直至兩項迭代計算均閉合收斂后,“積雪-土壤-砂礫石層”連續(xù)體水熱耦合計算完成。
2.2.3 模型邊界確定
模型上、下邊界的確定方式與WEP-COR模型類似?!胺e雪-土壤-砂礫石層”連續(xù)體的上邊界為大氣,但與上邊界直接接觸的不一定只是土壤。傳入連續(xù)體的熱量,在非積雪區(qū)由地表附近的大氣和表層土壤的溫度差以及表層土壤的水熱參數決定,在積雪區(qū)由地表附近的大氣和積雪的溫度差以及積雪的水熱參數決定,均采用強迫-恢復法計算。連續(xù)體的下邊界為過渡層或者地下水層,將下邊界溫度假設為一恒定的溫度作為底層邊界溫度。
2.2.4 積雪、土壤和砂礫石層水熱參數確定
(7)
(8)
Cs=2.09×103ρs
(9)
式中:ρs為積雪密度,kg/m3;λs為積雪的導熱系數,W/(m·℃);Cs為積雪的體積熱容量,J/(m3·℃)。
土壤、砂礫石層主要的水熱參數也包括體積熱容量、導熱系數和土壤導水率等,各參數計算公式如下。
體積熱容[35]:
Cv=(1-θs)Csg+θlCl+θiCi
(10)
導熱系數計算參考了IBIS 模型,具體如下[36]:
λ=λst(56θl+224θi)
λst=1.500ωrock+0.300ωsand+0.265ωsilt+0.250ωclat
(11)
導水率[29]:
(12)
式中:Csg、Cl和Ci分別為土壤(砂礫石層)、水、冰的體積熱容量,在0℃時土壤和砂礫石層分別為1.93×103J/(m3·℃)、3.1×103J/(m3·℃),水和冰分別為4.213×103J/(m3·℃)、1.94×103J/(m3·℃);λ、λst為土壤或砂礫石層實際的導熱系數、干燥狀態(tài)下的導熱系數,W/(m·℃);ωrock、ωsand、ωsilt和ωclat分別為巖石(礫石和卵石)、砂粒、粉粒和黏粒的體積比;KS為土壤或砂礫石層經溫度修正后的飽和含水率,cm/s。
參考陳仁升等的DWHC模型[37],不同溫度條件下土壤或砂礫石層KS計算方法如下:
(13)
式中:K0為常溫下的飽和導水率,cm/s;k0為凍結條件下最小的導水率,cm/s;Tf為土壤或砂礫石層最小導水率對應的臨界溫度,℃。本文考慮到土壤和砂礫石層水動力學特性的區(qū)別,對于土壤,k0按照原公式取值為0 cm/s;對于砂礫石層,由于孔隙較大,k0取值為大于0的值。
本文采用決定系數(R2)及平均根方差(ERMS)評價模型模擬效果。R2和ERMS的計算公式如下:
(14)
(15)
試驗點凍結融化期間空氣溫度、模擬雪厚、實測雪厚如圖7所示,試驗點積雪從2016年12月3日開始,2017年4月1日完全融化,最大雪厚12.4 cm,模擬雪厚與實測值比較吻合,土壤凍結融化期間分別對試驗點土壤溫度、含水率進行了模擬對比。
圖7 試驗點凍融期間溫度及積雪厚度Fig.7Temperature and snow thickness during freezing and thawing at the experimental site
3.2.1 溫度模擬
圖8是試驗點2016—2017年凍結融化期垂向不同深度的溫度模擬值和實測值的比較結果,總的來說,
圖8 不同深度土壤和砂礫石層溫度模擬與實測對比Fig.8Simulated and observed soil and loose rock strata temperatures at different depths
模擬溫度與實測溫度較為接近,但表層(10 cm、20 cm深)溫度波動較大,模擬效果稍差,R2分別為0.83、0.89,ERMS較大分別為4.07 ℃、3.39 ℃;表層以下(40 cm深及以下)溫度較為穩(wěn)定,模擬效果較好,R2在0.96以上,ERMS均小于10 cm和20 cm層,平均ERMS為2.86 ℃。表層溫度模擬效果較差,可能與氣溫波動較大有關,由于該流域氣象站點較少,溫度空間插值很難與當地實際氣溫相符。另外,80 cm深度以下的模擬溫度略小于實際溫度,可能與溫度下邊界條件的設置有關,需要在今后的研究中進一步改進。
3.2.2 含水率模擬
圖9為試驗點2016—2017年凍結融化期液態(tài)含水率模擬和實測值對比。由圖中可以看出,在凍結階段(11月—翌年1月),受溫度降低的影響,首先是上層土壤-砂礫石層的液態(tài)含水量變小,隨后下層砂礫石層的液態(tài)含水率變小,到1月下旬達到穩(wěn)定。土壤層(40 cm以上)由于其持水能力大于砂礫石層,初始的含水率更大,且受空氣溫度影響也大,在凍結過程中含水率減小明顯。2月下旬隨著氣溫回暖,先是上層土壤-砂礫石層開始融化,液態(tài)含水量增加,隨后土壤-砂礫石層底部也開始融化。在融化期,受積雪融化下滲影響,土壤-砂礫石層上部土壤層的含水量比下部砂礫石大,比凍結前也大。在快速凍結期(2016年12月18日前)R2均值為0.38,此后穩(wěn)定凍結期與融化期R2均值為0.58??偟膩碚f,模擬的分層土壤-砂礫石層的液態(tài)
圖9 不同深度土壤液態(tài)含水率模擬與實測對比Fig.9Simulated and observed soil moisture contents at different depths
含水率與實測值基本接近,可以反映凍結融化期各階段的變化趨勢。由于土壤和砂礫石層中組成成分的不確定性較大,模型進行概化時,認為砂礫石層和土壤是均質結構,無法準確反映土壤-砂礫石層持水能力的不穩(wěn)定,這也導致了模擬結果與實測值具有一定的誤差。
3.2.3 凍結深度模擬
試驗點2016—2017年凍結融化期凍結深度模擬結果見圖10。試驗點從11月下旬左右進入凍結期,3月上旬左右融通,凍結融化期總長在4個月左右,最大凍結深度約為67.2 cm。R2為0.76,ERMS為16.36 cm。根據實測凍結深度與模擬值的對比結果來看,試驗點處模擬的凍結規(guī)律與實測結果基本一致。模擬的凍結融化期起始時間與實測結果相比略微偏早,表層土壤開始融化時間比實測略晚。
圖10 土壤和砂礫石層凍結融化期凍結深度的模擬與實測對比Fig.10Simulated and observed soil and loose rock freezing depths
本文在WEP-COR模型凍土水熱耦合計算原理基礎上,結合青藏高原土壤層較薄、下伏砂礫石層較厚以及積雪期較長的特點加以改進,構建了“積雪-土壤-砂礫石層”連續(xù)體水熱耦合模型。
(1) 模型主要具有以下特點:將均質的土壤結構細化分解成土壤和砂礫石層2類介質,水熱耦合計算考慮了土壤和砂礫石層的凍融過程;在原模型土壤層的上方加入了一層積雪層,考慮了積雪對大氣與土壤間熱量傳遞的影響;將積雪、土壤和砂礫石層概化為一個緊密的一維連續(xù)結構,提出了“積雪-土壤-砂礫石層”連續(xù)體分層水熱耦合模擬方法。
(2) 基于尼洋河流域下游色季拉山試驗點的監(jiān)測結果對本文構建的模型進行驗證,各層溫度模擬R2均值為0.91,凍結融化期內含水率模擬R2均值為0.52,土壤凍結深度模擬R2值為0.76,土壤和砂礫石層不同深度的溫度、液態(tài)含水率和凍結深度的模擬結果與實測值接近。模型能較好地反映凍結融化期各階段土壤和砂礫石層的分層水熱變化趨勢。
(3) 由于高寒山區(qū)試驗條件惡劣,目前只設置了一個試驗站點,且本次取得的監(jiān)測資料時段較短。另外,因研究流域氣象站點較少,資料較為缺乏,加上模型的下邊界條件、“積雪-土壤-砂礫石層”連續(xù)體各項水熱參數均存在一定程度的概化,導致模擬結果與實測值之間仍存在一定偏差,這些都需要在未來的研究中進一步完善。
(4) 青藏高原“積雪-土壤-砂礫石層”連續(xù)體的周期性凍融過程影響著流域產流的各個環(huán)節(jié),對水文循環(huán)過程影響很大。為研究水熱耦合運移對流域水循環(huán)的影響還需將其與水文模型結合進一步分析。