周祖昊,劉揚(yáng)李,2,李玉慶,王鵬翔,,王 康,李 佳,朱熠明,,劉佳嘉,王富強(qiáng)

(1. 中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室,北京 100038；2. 中國電建集團(tuán)貴陽勘測設(shè)計(jì)研究院有限公司,貴州貴陽 550081；3. 西藏農(nóng)牧學(xué)院水利土木工程學(xué)院,西藏林芝 860000；4. 武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430072；5. 水利部南水北調(diào)規(guī)劃設(shè)計(jì)管理局,北京 100038；6. 華北水利水電大學(xué)水利學(xué)院,河南鄭州 450046)

目前,基于寒區(qū)水熱耦合的特點(diǎn),分布式水文模型已取得一定的進(jìn)展[1-2]。國際上常用的分布式水文模型都已對寒區(qū)的積雪和凍土過程進(jìn)行了一定程度的拓展表達(dá),部分模型中已加入凍土、融雪模塊來研究寒區(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]。當(dāng)前較成熟的凍土水熱耦合模型是基于非飽和土體水分遷移,用土壤未凍水含量隨溫度變化的函數(shù)來表征凍結(jié)機(jī)理的水動力學(xué)模型。其中以Harlan[5]提出的水熱耦合模型應(yīng)用最為廣泛；Taylor和Luthin[6]在其基礎(chǔ)上用未凍水含量結(jié)合土壤水分特征曲線計(jì)算土水勢,確定水分遷移的驅(qū)動力；Horiguchi[7]又加入了溫度梯度共同作用于水分遷移；美國農(nóng)業(yè)部Flerchinger和Saxton[8]綜合考慮氣相、水流遷移引起的熱量變化以及邊界處的蒸發(fā)、輻射等,建立了垂直一維水熱耦合模型,較全面地反映了影響凍土水熱變化因素；雷志棟等[9]對各向同性土柱的一維凍結(jié)問題簡化了模擬過程,忽略了水汽遷移、熱量對流作用,提出了凍土非飽和土壤水分遷移方程,對地下水淺埋條件下的土壤凍融過程進(jìn)行了模擬,模型模擬效果較好。水動力學(xué)模型參數(shù)簡便易得,是較為成熟的水熱耦合模型。

青藏高原是長江、黃河、瀾滄江、雅魯藏布江等多個(gè)大江大河的發(fā)源地。作為世界面積最大、海拔最高的高海拔寒區(qū),平均海拔4 376 m[10],其廣泛分布的積雪和凍土影響了該地區(qū)土壤介質(zhì)的儲水、導(dǎo)水、產(chǎn)水和熱傳導(dǎo)過程[11]。與中國其他地區(qū)相比,青藏高原土壤層較薄,土壤層下存在較厚的砂礫石層。與土壤相比,砂礫石具有不同的水、熱性質(zhì)[12],孔隙度、密度也有所不同,對土壤含水量和導(dǎo)水率有較大影響[13-14]。對青藏高原的水熱運(yùn)移模擬,部分陸面模式雖然考慮了砂礫石影響,通過砂礫石占比對土壤熱傳導(dǎo)率、土壤質(zhì)地和地形等敏感參數(shù)進(jìn)行調(diào)整,模擬效果會有一定提升[15-17],但由于模型將土壤概化為一維均質(zhì)結(jié)構(gòu),僅通過參數(shù)調(diào)整難以體現(xiàn)該地區(qū)這種上下分層明顯的地質(zhì)結(jié)構(gòu)對水熱運(yùn)移的影響,模擬結(jié)果均有一定偏差。因此,在研究青藏高原土壤凍融與水循環(huán)過程時(shí),不僅要考慮土壤層內(nèi)的水熱遷移,對于土壤層下伏較厚的砂礫石層也應(yīng)一并考慮。與此同時(shí),積雪較好的隔熱性和對太陽短波的反射性,也會影響土壤水熱遷移。在水熱耦合研究過程中,要把土壤和積雪作為一個(gè)復(fù)雜的、相互作用的聯(lián)合體來統(tǒng)一考慮[18]。對積雪的研究,根據(jù)其模擬的復(fù)雜程度與否,大致可分為以下兩類:一是利用強(qiáng)迫-恢復(fù)(force-restore)法[19-20],通過相對簡單的單層積雪模型,區(qū)分積雪和土壤的熱力學(xué)性質(zhì),模擬積雪-土壤復(fù)合層的溫度變化,以度日因子法計(jì)算融雪量[21]；二是以物理機(jī)制為基礎(chǔ)的復(fù)雜模型,此類模型考慮了積雪內(nèi)部的質(zhì)量守衡和能量平衡以及積雪表面與大氣環(huán)境條件的相互作用,將積雪細(xì)分為多層,對積雪內(nèi)部相變、水分運(yùn)移等進(jìn)行了精細(xì)描述[22-23]。但后一類模型太過復(fù)雜,影響了其推廣與應(yīng)用。

本文針對青藏高原的氣候和地質(zhì)特點(diǎn),在土壤凍結(jié)融化期開展了野外土壤和砂礫石層水熱監(jiān)測試驗(yàn),以WEP-COR模型[1]的土壤水熱耦合模塊為基礎(chǔ),構(gòu)建基于“積雪-土壤-砂礫石層”連續(xù)體的水熱耦合模型,并采用試驗(yàn)觀測結(jié)果對模型進(jìn)行驗(yàn)證。

1 研究區(qū)與資料

1.1 研究區(qū)選取

選取雅魯藏布江的一級支流尼洋河流域作為本文研究的典型區(qū)域。尼洋河位于西藏東南林芝地區(qū)雅魯藏布江中下游左岸,北緯29°28′—30°38′,東經(jīng)92°10′—94°35′,源于西藏自治區(qū)米拉山西側(cè)的錯木梁拉,河源海拔約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 試驗(yàn)布設(shè)

1.2.1 土壤厚度和成分調(diào)查與取樣分析

根據(jù)尼洋河流域?qū)嵉乜疾?圖1)和部分工程設(shè)計(jì)的地質(zhì)勘測資料,該流域具有土壤層厚度較薄、土壤下方的砂礫石層厚度較大的特點(diǎn)。尼洋河流域的砂礫石層主要包括礫石層和卵石層,礫石層主要為圓狀礫石,含卵石顆粒,礫石質(zhì)量百分比約50%～65%,黏粒質(zhì)量百分比5%～10%,孔隙間填充中細(xì)砂；卵石層孔隙間填充圓礫和中細(xì)砂,局部含大漂石、孤石。

針對尼洋河流域這一地質(zhì)特點(diǎn),為研究流域內(nèi)土壤層從山腳到山頂?shù)暮穸茸兓?guī)律,從尼洋河流域的源頭到河口選取了海拔不同的24個(gè)取樣點(diǎn)。其中點(diǎn)1—點(diǎn)16沿河道,點(diǎn)17—點(diǎn)24從山腳到山頂測量的土壤厚度(圖2)。對24個(gè)取樣點(diǎn)的土壤厚度及組分進(jìn)行了測量和分析,尼洋河流域土壤主要為砂質(zhì)壤土和壤質(zhì)砂土,平均砂粒質(zhì)量百分比為55.89%,粉砂粒質(zhì)量百分比為31.2%,黏粒質(zhì)量百分比為12.91%,厚度隨著海拔升高而減小,平均厚度為53.7 cm。

圖1 土壤-砂礫石層結(jié)構(gòu)Fig.1Soil-sand gravel structure

圖2 尼洋河流域土壤取樣點(diǎn)、試驗(yàn)點(diǎn)及氣象站位置Fig.2Soil sampling points and meteorological stations

1.2.2 凍結(jié)融化期“積雪-土壤-砂礫石層”監(jiān)測試驗(yàn)

本次研究在尼洋河流域下游色季拉山的山腰上選擇一處開展季節(jié)凍結(jié)融化期水熱耦合過程監(jiān)測試驗(yàn)。試驗(yàn)點(diǎn)(94°21′45″E,29°27′12″N)海拔為4 607 m,試驗(yàn)期為2016年11月—2017年4月,試驗(yàn)期間試驗(yàn)點(diǎn)降水8.3 mm,平均氣溫-3.4 ℃,最高6.3 ℃,最低-11.7 ℃。凍結(jié)融化期,試驗(yàn)點(diǎn)位置見圖2。試驗(yàn)人員在試驗(yàn)點(diǎn)開挖了一個(gè)深度為1.6 m的試驗(yàn)坑,在地面以下垂直深度上每隔10 cm安裝自動測定液態(tài)含水率的TDR傳感器、測量溫度的PT100傳感器和測量基質(zhì)勢的TensionMark傳感器,對下墊面凍結(jié)和融化過程中水、熱、勢能過程進(jìn)行自動監(jiān)測,并根據(jù)PT100傳感器和TensionMark傳感器測量結(jié)果推算土壤凍結(jié)深度,同時(shí)于開始降雪后間隔1～2周根據(jù)降雪情況通過米尺實(shí)地測量下墊面覆蓋的雪層厚度。試驗(yàn)前,采用核磁共振方法對高原季節(jié)性凍融條件下水、熱過程監(jiān)測的儀器設(shè)備進(jìn)行率定。儀器完成安裝后進(jìn)行原狀土回填,保證回填土壤與實(shí)際土壤的參數(shù)一致。

1.3 氣象資料

氣象條件是水熱耦合模擬的邊界條件。由于流域內(nèi)氣象站點(diǎn)稀少,且高程差別較大,試驗(yàn)點(diǎn)的逐日氣溫采用考慮高程修正的RDS(Reversed Distance Squared)方法由氣象站數(shù)據(jù)插值得到,降水?dāng)?shù)據(jù)采用平面插值、結(jié)合西藏自治區(qū)水資源公報(bào)[25]中降水等值線修正得到,逐日相對濕度、日照時(shí)速和風(fēng)速采用RDS方法插值得到。本文用到的基本氣象數(shù)據(jù)主要為尼洋河流域境內(nèi)的林芝站和周邊嘉黎站(圖2)的降水、氣溫、相對濕度、日照時(shí)數(shù)和風(fēng)速。氣象數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn),使用數(shù)據(jù)時(shí)間序列為2016—2017年,時(shí)間尺度為日。

2 “積雪-土壤-砂礫石層”連續(xù)體水熱耦合模型構(gòu)建

本文在WEP-COR模型中凍土水熱耦合計(jì)算原理基礎(chǔ)上,結(jié)合青藏高原特點(diǎn)加以改進(jìn),構(gòu)建包含“積雪-土壤-砂礫石層”連續(xù)體的高寒山區(qū)水熱耦合模型。

2.1 WEP-COR中的凍土水熱耦合計(jì)算原理

WEP-COR模型在保留原WEP-L模型關(guān)于地表與大氣間的能量交換計(jì)算的基礎(chǔ)上,系統(tǒng)添加了土壤層間的水熱通量和相關(guān)水熱參數(shù)的計(jì)算過程,并且將活動土壤層細(xì)分為11層,見圖3。其中,Rs為地表產(chǎn)流;Ri為第i層土壤的橫向流或壤中流,Ri與坡度和土壤含水率有關(guān);E1為土壤蒸發(fā)量;Er為植被蒸騰量;Qi為第i層土壤的重力排水;P為降水;Ta為大氣溫度;Ti為第i層土壤的溫度;Gi為第i層土壤與相鄰?fù)寥缹娱g由溫度差異引起的熱通量。水熱耦合模擬采用顯式差分進(jìn)行數(shù)值迭代計(jì)算,模型模擬的時(shí)間步長為1 d。

圖3 WEP-COR模型分層水熱通量計(jì)算結(jié)構(gòu)示意Fig.3WEP-COR model water and heat flux simulation structure

2.1.1 土壤水熱運(yùn)移方程

模型假設(shè)土壤凍融時(shí)只有液態(tài)水發(fā)生運(yùn)移,土壤水分的運(yùn)移主要受重力勢、基質(zhì)勢和溫度勢的影響,因此添加的土壤一維垂直水分流方程可以寫成[26]:

(1)

根據(jù)能量平衡原理,凍融系統(tǒng)中每一層土壤的能量變化都用于系統(tǒng)內(nèi)的土壤溫變和水分相變,溫度勢是水分相變的驅(qū)動力,而大氣與表層土壤的溫差則是熱傳導(dǎo)的源動力。假設(shè)土壤各向均質(zhì)同性,并忽略土壤中的水汽遷移,添加的一維垂直流動的熱流運(yùn)動基本方程可變形為[27]

(2)

式中:θl、θi分別為土壤中液態(tài)水、冰的體積百分比,cm3/cm3；T為土壤溫度,℃；t、z分別為時(shí)間、空間坐標(biāo)(垂直向下為正)；D(θl)、K(θl)分別為非飽和土壤水分?jǐn)U散率、導(dǎo)水率,cm2/s；Cv為土壤體積熱容量,J/(m3·℃)；λ為土壤導(dǎo)熱系數(shù),W/(m·℃)；ρi、ρl分別為冰、水密度,kg/m3；Lf為融化潛熱,Lf=3.35×106J/kg。

2.1.2 土壤水熱聯(lián)系方程

土壤凍融過程中,凍土水熱運(yùn)動間的聯(lián)系主要表現(xiàn)在未凍水的含水率與土壤負(fù)溫的動態(tài)平衡中[19-20]:

θl=θm(T)

(3)

式中:θm(T)為土壤負(fù)溫對應(yīng)的最大未凍水含水率。

2.2 “積雪-土壤-砂礫石層”連續(xù)體水熱耦合模型原理

2.2.1 模型結(jié)構(gòu)

積雪作為一種良好的隔熱體,具有熱容量大、導(dǎo)熱系數(shù)較小等特點(diǎn),對大氣和土壤間的熱傳遞具有明顯的削弱作用。試驗(yàn)點(diǎn)所處位置積雪持續(xù)時(shí)間較長、厚度穩(wěn)定,對土壤水熱的影響不容忽視。根據(jù)尼洋河流域?qū)嵉乜疾?圖1)和基于26個(gè)點(diǎn)的測量結(jié)果建立的土壤厚度—高程關(guān)系圖(圖4),發(fā)現(xiàn)尼洋河流域土壤厚度整體較薄,在土壤層下方存在厚度較大、混合砂礫石和少量土壤的砂礫石層。而且土壤層厚度在河谷區(qū)基本穩(wěn)定,在山坡上隨著高程增加呈現(xiàn)變薄的趨勢,在山頂附近厚度又基本保持不變。在試驗(yàn)點(diǎn)處土壤層為砂質(zhì)壤土,砂粒質(zhì)量百分比為61.92%,粉砂粒質(zhì)量百分比為26.2%,黏粒質(zhì)量百分比為11.88%,土壤層厚度為40 cm；土壤下砂礫石層礫石質(zhì)量百分比為53.7%,孔隙間填充中細(xì)砂。土壤與砂礫石層水分特征參數(shù)見表1。

圖4 尼洋河流域土壤厚度—海拔關(guān)系Fig.4Relationship between soil thickness and elevation in Niyang River basin

表1 尼洋河流域土壤水動力參數(shù)

根據(jù)尼洋河流域氣候和地質(zhì)特征,本文將水熱耦合模擬的對象定義為“積雪-土壤-砂礫石層”連續(xù)體(圖5)。模型計(jì)算結(jié)構(gòu)上,在WEP-COR模型11層土壤層分層的基礎(chǔ)上,于土壤層上方加入一層積雪層,同時(shí)將原有的11層土壤層細(xì)分為土壤和砂礫石層2類介質(zhì),計(jì)算結(jié)構(gòu)如圖6所示。第0層為積雪,根據(jù)等高帶在計(jì)算單元中所處位置設(shè)置上部i層為土壤層,下部12-i層為砂礫石層。由于距離地表越近,水、熱變化越快,故第1、2層厚度設(shè)為10 cm,第3～11層厚度設(shè)為20 cm。值得注意的是,土壤-砂礫石層各層的厚度需根據(jù)所在位置含水層的實(shí)際厚度進(jìn)行修正,也就是說,當(dāng)含水層總厚度小于11層土壤-砂礫石層厚度時(shí),須根據(jù)含水層實(shí)際厚度確定土壤-砂礫石層的計(jì)算層數(shù)和最后一層的厚度。

圖5 “積雪-土壤-砂礫石層”垂向分層概化Fig.5Vertical layering of “snow-soil-sand gravel layer”

圖6 分層水熱通量計(jì)算結(jié)構(gòu)示意Fig.6Layered calculation structure of water and heat flux

2.2.2 計(jì)算原理

采用度日因子法[28]計(jì)算積雪融雪過程,具體如下:

M=k(Ta-T0)

(4)

式中:M為當(dāng)日融雪當(dāng)量,mm；k為積雪融化的度日因子,mm/(℃·d)；T0為積雪融化臨界溫度,℃,本文取-1 ℃。

對于積雪的度日因子,一般情況下在1～7 mm/(℃·d),本次研究根據(jù)下墊面情況取為4 mm/(℃·d)。由于青藏高原地區(qū)晝夜溫差較其他地區(qū)偏大,本次研究的雨雪臨界氣溫設(shè)為2 ℃。

積雪熱量平衡公式[29]:

G=Le+Rn+Q1+Q2

(5)

式中:G為積雪層的能量變化,J/(m2·d)；Le為積雪的蒸發(fā)、融化和升華潛熱通量,J/(m2·d)；Rn為積雪層表面的凈輻射量,J/(m2·d),由到達(dá)積雪表面的凈長波輻射與凈短波輻射計(jì)算而得；Q1為積雪與大氣間的熱量交換量,J/(m2·d),采用強(qiáng)迫-恢復(fù)法計(jì)算；Q2為積雪與土壤間的熱量交換量,J/(m2·d)。

積雪和土壤間的熱量交換量采用以下公式[30-31]:

(6)

式中:ZS為積雪厚度,m；ZC為第一層土壤厚度,m；λS為冰雪層導(dǎo)熱系數(shù),W/(m·℃)；λC為土壤層導(dǎo)熱系數(shù),W/(m·℃)；RC為土壤層與積雪層的接觸熱阻,(m2·℃)/W；TC為第一層土壤層溫度,℃；TS為積雪層溫度,℃。

土壤層之間的水熱運(yùn)移方程和聯(lián)系方程采用式(1)和式(2),砂礫石層之間以及砂礫石層和土壤層之間的水熱通量計(jì)算采用與土壤層相同的計(jì)算原理,但砂礫石層的各項(xiàng)水熱參數(shù)與土壤不同。

模型的數(shù)值模擬仍采用顯式差分進(jìn)行數(shù)值迭代計(jì)算,模擬步長為1 d。每個(gè)時(shí)間步長內(nèi)模型先根據(jù)初始條件從積雪層到底層砂礫石層逐層計(jì)算熱傳導(dǎo)、溫度和水分相變,然后以該層溫度為判定條件進(jìn)行迭代計(jì)算直至收斂(步驟1)；熱量計(jì)算閉合收斂后,再根據(jù)水量平衡計(jì)算各層水分運(yùn)移量,并修正各土壤和砂礫石層(積雪層不考慮)的含水率和含冰率(步驟2)；用液態(tài)含水率判定是否收斂,若不收斂則返回步驟1進(jìn)行熱量計(jì)算,直至兩項(xiàng)迭代計(jì)算均閉合收斂后,“積雪-土壤-砂礫石層”連續(xù)體水熱耦合計(jì)算完成。

2.2.3 模型邊界確定

模型上、下邊界的確定方式與WEP-COR模型類似?！胺e雪-土壤-砂礫石層”連續(xù)體的上邊界為大氣,但與上邊界直接接觸的不一定只是土壤。傳入連續(xù)體的熱量,在非積雪區(qū)由地表附近的大氣和表層土壤的溫度差以及表層土壤的水熱參數(shù)決定,在積雪區(qū)由地表附近的大氣和積雪的溫度差以及積雪的水熱參數(shù)決定,均采用強(qiáng)迫-恢復(fù)法計(jì)算。連續(xù)體的下邊界為過渡層或者地下水層,將下邊界溫度假設(shè)為一恒定的溫度作為底層邊界溫度。

2.2.4 積雪、土壤和砂礫石層水熱參數(shù)確定

(7)

(8)

Cs=2.09×103ρs

(9)

式中:ρs為積雪密度,kg/m3；λs為積雪的導(dǎo)熱系數(shù),W/(m·℃)；Cs為積雪的體積熱容量,J/(m3·℃)。

土壤、砂礫石層主要的水熱參數(shù)也包括體積熱容量、導(dǎo)熱系數(shù)和土壤導(dǎo)水率等,各參數(shù)計(jì)算公式如下。

體積熱容[35]:

Cv=(1-θs)Csg+θlCl+θiCi

(10)

導(dǎo)熱系數(shù)計(jì)算參考了IBIS 模型,具體如下[36]:

λ=λst(56θl+224θi)
λst=1.500ωrock+0.300ωsand+0.265ωsilt+0.250ωclat

(11)

導(dǎo)水率[29]:

(12)

式中:Csg、Cl和Ci分別為土壤(砂礫石層)、水、冰的體積熱容量,在0℃時(shí)土壤和砂礫石層分別為1.93×103J/(m3·℃)、3.1×103J/(m3·℃),水和冰分別為4.213×103J/(m3·℃)、1.94×103J/(m3·℃)；λ、λst為土壤或砂礫石層實(shí)際的導(dǎo)熱系數(shù)、干燥狀態(tài)下的導(dǎo)熱系數(shù),W/(m·℃)；ωrock、ωsand、ωsilt和ωclat分別為巖石(礫石和卵石)、砂粒、粉粒和黏粒的體積比；KS為土壤或砂礫石層經(jīng)溫度修正后的飽和含水率,cm/s。

參考陳仁升等的DWHC模型[37],不同溫度條件下土壤或砂礫石層KS計(jì)算方法如下:

(13)

式中:K0為常溫下的飽和導(dǎo)水率,cm/s；k0為凍結(jié)條件下最小的導(dǎo)水率,cm/s；Tf為土壤或砂礫石層最小導(dǎo)水率對應(yīng)的臨界溫度,℃。本文考慮到土壤和砂礫石層水動力學(xué)特性的區(qū)別,對于土壤,k0按照原公式取值為0 cm/s；對于砂礫石層,由于孔隙較大,k0取值為大于0的值。

3 模型驗(yàn)證與結(jié)果分析

3.1 評價(jià)方法

本文采用決定系數(shù)(R2)及平均根方差(ERMS)評價(jià)模型模擬效果。R2和ERMS的計(jì)算公式如下:

(14)

(15)

3.2 土壤和砂礫石層水熱耦合模擬

試驗(yàn)點(diǎn)凍結(jié)融化期間空氣溫度、模擬雪厚、實(shí)測雪厚如圖7所示,試驗(yàn)點(diǎn)積雪從2016年12月3日開始,2017年4月1日完全融化,最大雪厚12.4 cm,模擬雪厚與實(shí)測值比較吻合,土壤凍結(jié)融化期間分別對試驗(yàn)點(diǎn)土壤溫度、含水率進(jìn)行了模擬對比。

圖7 試驗(yàn)點(diǎn)凍融期間溫度及積雪厚度Fig.7Temperature and snow thickness during freezing and thawing at the experimental site

3.2.1 溫度模擬

圖8是試驗(yàn)點(diǎn)2016—2017年凍結(jié)融化期垂向不同深度的溫度模擬值和實(shí)測值的比較結(jié)果,總的來說,

圖8 不同深度土壤和砂礫石層溫度模擬與實(shí)測對比Fig.8Simulated and observed soil and loose rock strata temperatures at different depths

模擬溫度與實(shí)測溫度較為接近,但表層(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 ℃。表層溫度模擬效果較差,可能與氣溫波動較大有關(guān),由于該流域氣象站點(diǎn)較少,溫度空間插值很難與當(dāng)?shù)貙?shí)際氣溫相符。另外,80 cm深度以下的模擬溫度略小于實(shí)際溫度,可能與溫度下邊界條件的設(shè)置有關(guān),需要在今后的研究中進(jìn)一步改進(jìn)。

3.2.2 含水率模擬

圖9為試驗(yàn)點(diǎn)2016—2017年凍結(jié)融化期液態(tài)含水率模擬和實(shí)測值對比。由圖中可以看出,在凍結(jié)階段(11月—翌年1月),受溫度降低的影響,首先是上層土壤-砂礫石層的液態(tài)含水量變小,隨后下層砂礫石層的液態(tài)含水率變小,到1月下旬達(dá)到穩(wěn)定。土壤層(40 cm以上)由于其持水能力大于砂礫石層,初始的含水率更大,且受空氣溫度影響也大,在凍結(jié)過程中含水率減小明顯。2月下旬隨著氣溫回暖,先是上層土壤-砂礫石層開始融化,液態(tài)含水量增加,隨后土壤-砂礫石層底部也開始融化。在融化期,受積雪融化下滲影響,土壤-砂礫石層上部土壤層的含水量比下部砂礫石大,比凍結(jié)前也大。在快速凍結(jié)期(2016年12月18日前)R2均值為0.38,此后穩(wěn)定凍結(jié)期與融化期R2均值為0.58。總的來說,模擬的分層土壤-砂礫石層的液態(tài)

圖9 不同深度土壤液態(tài)含水率模擬與實(shí)測對比Fig.9Simulated and observed soil moisture contents at different depths

含水率與實(shí)測值基本接近,可以反映凍結(jié)融化期各階段的變化趨勢。由于土壤和砂礫石層中組成成分的不確定性較大,模型進(jìn)行概化時(shí),認(rèn)為砂礫石層和土壤是均質(zhì)結(jié)構(gòu),無法準(zhǔn)確反映土壤-砂礫石層持水能力的不穩(wěn)定,這也導(dǎo)致了模擬結(jié)果與實(shí)測值具有一定的誤差。

3.2.3 凍結(jié)深度模擬

試驗(yàn)點(diǎn)2016—2017年凍結(jié)融化期凍結(jié)深度模擬結(jié)果見圖10。試驗(yàn)點(diǎn)從11月下旬左右進(jìn)入凍結(jié)期,3月上旬左右融通,凍結(jié)融化期總長在4個(gè)月左右,最大凍結(jié)深度約為67.2 cm。R2為0.76,ERMS為16.36 cm。根據(jù)實(shí)測凍結(jié)深度與模擬值的對比結(jié)果來看,試驗(yàn)點(diǎn)處模擬的凍結(jié)規(guī)律與實(shí)測結(jié)果基本一致。模擬的凍結(jié)融化期起始時(shí)間與實(shí)測結(jié)果相比略微偏早,表層土壤開始融化時(shí)間比實(shí)測略晚。

圖10 土壤和砂礫石層凍結(jié)融化期凍結(jié)深度的模擬與實(shí)測對比Fig.10Simulated and observed soil and loose rock freezing depths

4 結(jié) 論

本文在WEP-COR模型凍土水熱耦合計(jì)算原理基礎(chǔ)上,結(jié)合青藏高原土壤層較薄、下伏砂礫石層較厚以及積雪期較長的特點(diǎn)加以改進(jìn),構(gòu)建了“積雪-土壤-砂礫石層”連續(xù)體水熱耦合模型。

(1) 模型主要具有以下特點(diǎn):將均質(zhì)的土壤結(jié)構(gòu)細(xì)化分解成土壤和砂礫石層2類介質(zhì),水熱耦合計(jì)算考慮了土壤和砂礫石層的凍融過程；在原模型土壤層的上方加入了一層積雪層,考慮了積雪對大氣與土壤間熱量傳遞的影響；將積雪、土壤和砂礫石層概化為一個(gè)緊密的一維連續(xù)結(jié)構(gòu),提出了“積雪-土壤-砂礫石層”連續(xù)體分層水熱耦合模擬方法。

(2) 基于尼洋河流域下游色季拉山試驗(yàn)點(diǎn)的監(jiān)測結(jié)果對本文構(gòu)建的模型進(jìn)行驗(yàn)證,各層溫度模擬R2均值為0.91,凍結(jié)融化期內(nèi)含水率模擬R2均值為0.52,土壤凍結(jié)深度模擬R2值為0.76,土壤和砂礫石層不同深度的溫度、液態(tài)含水率和凍結(jié)深度的模擬結(jié)果與實(shí)測值接近。模型能較好地反映凍結(jié)融化期各階段土壤和砂礫石層的分層水熱變化趨勢。

(3) 由于高寒山區(qū)試驗(yàn)條件惡劣,目前只設(shè)置了一個(gè)試驗(yàn)站點(diǎn),且本次取得的監(jiān)測資料時(shí)段較短。另外,因研究流域氣象站點(diǎn)較少,資料較為缺乏,加上模型的下邊界條件、“積雪-土壤-砂礫石層”連續(xù)體各項(xiàng)水熱參數(shù)均存在一定程度的概化,導(dǎo)致模擬結(jié)果與實(shí)測值之間仍存在一定偏差,這些都需要在未來的研究中進(jìn)一步完善。

(4) 青藏高原“積雪-土壤-砂礫石層”連續(xù)體的周期性凍融過程影響著流域產(chǎn)流的各個(gè)環(huán)節(jié),對水文循環(huán)過程影響很大。為研究水熱耦合運(yùn)移對流域水循環(huán)的影響還需將其與水文模型結(jié)合進(jìn)一步分析。