基于水熱耦合的青藏高原分布式水文模型
——I. “積雪-土壤-砂礫石層”連續(xù)體水熱耦合模擬

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

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

目前,基于寒區(qū)水熱耦合的特點,分布式水文模型已取得一定的進(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é)合土壤水分特征曲線計算土水勢,確定水分遷移的驅(qū)動力；Horiguchi[7]又加入了溫度梯度共同作用于水分遷移；美國農(nóng)業(yè)部Flerchinger和Saxton[8]綜合考慮氣相、水流遷移引起的熱量變化以及邊界處的蒸發(fā)、輻射等,建立了垂直一維水熱耦合模型,較全面地反映了影響凍土水熱變化因素；雷志棟等[9]對各向同性土柱的一維凍結(jié)問題簡化了模擬過程,忽略了水汽遷移、熱量對流作用,提出了凍土非飽和土壤水分遷移方程,對地下水淺埋條件下的土壤凍融過程進(jìn)行了模擬,模型模擬效果較好。水動力學(xué)模型參數(shù)簡便易得,是較為成熟的水熱耦合模型。

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

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

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

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

針對尼洋河流域這一地質(zhì)特點,為研究流域內(nèi)土壤層從山腳到山頂?shù)暮穸茸兓?guī)律,從尼洋河流域的源頭到河口選取了海拔不同的24個取樣點。其中點1—點16沿河道,點17—點24從山腳到山頂測量的土壤厚度(圖2)。對24個取樣點的土壤厚度及組分進(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 尼洋河流域土壤取樣點、試驗點及氣象站位置Fig.2Soil sampling points and meteorological stations

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

本次研究在尼洋河流域下游色季拉山的山腰上選擇一處開展季節(jié)凍結(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 ℃。凍結(jié)融化期,試驗點位置見圖2。試驗人員在試驗點開挖了一個深度為1.6 m的試驗坑,在地面以下垂直深度上每隔10 cm安裝自動測定液態(tài)含水率的TDR傳感器、測量溫度的PT100傳感器和測量基質(zhì)勢的TensionMark傳感器,對下墊面凍結(jié)和融化過程中水、熱、勢能過程進(jìn)行自動監(jiān)測,并根據(jù)PT100傳感器和TensionMark傳感器測量結(jié)果推算土壤凍結(jié)深度,同時于開始降雪后間隔1～2周根據(jù)降雪情況通過米尺實地測量下墊面覆蓋的雪層厚度。試驗前,采用核磁共振方法對高原季節(jié)性凍融條件下水、熱過程監(jiān)測的儀器設(shè)備進(jìn)行率定。儀器完成安裝后進(jìn)行原狀土回填,保證回填土壤與實際土壤的參數(shù)一致。

1.3 氣象資料

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

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

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

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

WEP-COR模型在保留原WEP-L模型關(guān)于地表與大氣間的能量交換計算的基礎(chǔ)上,系統(tǒng)添加了土壤層間的水熱通量和相關(guān)水熱參數(shù)的計算過程,并且將活動土壤層細(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ù)值迭代計算,模型模擬的時間步長為1 d。

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

2.1.1 土壤水熱運移方程

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

(1)

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

(2)

式中:θl、θi分別為土壤中液態(tài)水、冰的體積百分比,cm3/cm3；T為土壤溫度,℃；t、z分別為時間、空間坐標(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)系方程

土壤凍融過程中,凍土水熱運動間的聯(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ù)較小等特點,對大氣和土壤間的熱傳遞具有明顯的削弱作用。試驗點所處位置積雪持續(xù)時間較長、厚度穩(wěn)定,對土壤水熱的影響不容忽視。根據(jù)尼洋河流域?qū)嵉乜疾?圖1)和基于26個點的測量結(jié)果建立的土壤厚度—高程關(guān)系圖(圖4),發(fā)現(xiàn)尼洋河流域土壤厚度整體較薄,在土壤層下方存在厚度較大、混合砂礫石和少量土壤的砂礫石層。而且土壤層厚度在河谷區(qū)基本穩(wěn)定,在山坡上隨著高程增加呈現(xià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)。模型計算結(jié)構(gòu)上,在WEP-COR模型11層土壤層分層的基礎(chǔ)上,于土壤層上方加入一層積雪層,同時將原有的11層土壤層細(xì)分為土壤和砂礫石層2類介質(zhì),計算結(jié)構(gòu)如圖6所示。第0層為積雪,根據(jù)等高帶在計算單元中所處位置設(shè)置上部i層為土壤層,下部12-i層為砂礫石層。由于距離地表越近,水、熱變化越快,故第1、2層厚度設(shè)為10 cm,第3～11層厚度設(shè)為20 cm。值得注意的是,土壤-砂礫石層各層的厚度需根據(jù)所在位置含水層的實際厚度進(jìn)行修正,也就是說,當(dāng)含水層總厚度小于11層土壤-砂礫石層厚度時,須根據(jù)含水層實際厚度確定土壤-砂礫石層的計算層數(shù)和最后一層的厚度。

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

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

2.2.2 計算原理

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

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á)積雪表面的凈長波輻射與凈短波輻射計算而得；Q1為積雪與大氣間的熱量交換量,J/(m2·d),采用強(qiáng)迫-恢復(fù)法計算；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為積雪層溫度,℃。

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

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

2.2.3 模型邊界確定

模型上、下邊界的確定方式與WEP-COR模型類似。“積雪-土壤-砂礫石層”連續(xù)體的上邊界為大氣,但與上邊界直接接觸的不一定只是土壤。傳入連續(xù)體的熱量,在非積雪區(qū)由地表附近的大氣和表層土壤的溫度差以及表層土壤的水熱參數(shù)決定,在積雪區(qū)由地表附近的大氣和積雪的溫度差以及積雪的水熱參數(shù)決定,均采用強(qiáng)迫-恢復(fù)法計算。連續(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ù)計算公式如下。

體積熱容[35]:

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

(10)

導(dǎo)熱系數(shù)計算參考了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℃時土壤和砂礫石層分別為1.93×103J/(m3·℃)、3.1×103J/(m3·℃),水和冰分別為4.213×103J/(m3·℃)、1.94×103J/(m3·℃)；λ、λst為土壤或砂礫石層實際的導(dǎo)熱系數(shù)、干燥狀態(tài)下的導(dǎo)熱系數(shù),W/(m·℃)；ωrock、ωsand、ωsilt和ωclat分別為巖石(礫石和卵石)、砂粒、粉粒和黏粒的體積比；KS為土壤或砂礫石層經(jīng)溫度修正后的飽和含水率,cm/s。

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

(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 模型驗證與結(jié)果分析

3.1 評價方法

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

(14)

(15)

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

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

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

3.2.1 溫度模擬

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

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

3.2.2 含水率模擬

圖9為試驗點2016—2017年凍結(jié)融化期液態(tài)含水率模擬和實測值對比。由圖中可以看出,在凍結(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?？偟膩碚f,模擬的分層土壤-砂礫石層的液態(tài)

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

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

3.2.3 凍結(jié)深度模擬

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

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

4 結(jié) 論

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

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

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

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

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