基于能量平衡的融雪期積雪深度模擬

黃東浩，王 鵬，周麗麗，孫宏超

（沈陽農(nóng)業(yè)大學水利學院/遼寧省水土流失防控與生態(tài)修復重點實驗室，沈陽 110161）

東北黑土區(qū)地處我國中高緯度地區(qū)，雪期長度約8個月[1]，降雪量占年降水量7%～25%，由此產(chǎn)生的融雪徑流量和輸沙模數(shù)分別占全年的13.3%～24.9%和5.8%～27.7%[2]，對該區(qū)域水土可持續(xù)利用具有重要影響。隨著全球氣候變暖，該區(qū)融雪期提前與土壤凍融循環(huán)次數(shù)增多，從而加劇土壤侵蝕[3-4]；同時消融早期的融雪徑流明顯增加，導致后期積雪減少，改變了流域融雪徑流年內(nèi)分配[5]，對流域內(nèi)水資源的管理、規(guī)劃及可持續(xù)利用造成困難[6]。因此，準確預報積雪消融過程對該區(qū)水資源管理及融雪侵蝕研究具有重要的指導意義。

積雪消融預報主要依靠融雪模型，由于模型輸出表達方式的不同，模型的復雜程度也有所不同[7]，如CROCUS[8]模型可以模擬雪內(nèi)部結(jié)構變化結(jié)構特征。就目前而言，國內(nèi)外積雪消融預報模型主要分為兩類：一類是基于度日因子的概念模型；另一類是基于物理過程的能量模型。度日因子模型由于輸入?yún)?shù)少、計算簡單而被廣泛應用[9]，但度日模型存在的缺點[10]為：(1)隨著時間的推移，模擬精度明顯下降，不適用于長時間序列的融雪模擬；(2)不能定量描述積雪內(nèi)部消融狀況。針對以上不足，一些學者進行了模型優(yōu)化，AKBARI 等[11]通過優(yōu)化積雪覆蓋面積和引入流量沉降系數(shù)，來提高度日模型的模擬精度；穆艾塔爾·賽地等[12]基于地形因子將流域分帶后引入度日模型來模擬積雪消融變化。但就本質(zhì)而言度日因子模型是一個數(shù)學統(tǒng)計模型，由于積雪消融過程中雪表面時空差異較大，該模型無法描述積雪內(nèi)部的消融過程，若有降雨輸入時還會降低精度[13]。而基于能量平衡的融雪模型可定量描述融雪過程，并漸成研究趨勢[14]，模型輸入?yún)?shù)要求也較高[15]，目前研究方向主要有：(1)將能量平衡融雪模型與其他模型集成，孟現(xiàn)勇等[16]將基于能量平衡的融雪模型耦合到SWAT模型融雪模塊，提升了雪深、徑流量的模擬精度；CHEN 等[17]利用SNTHERM-MEMLS 耦合模型來模擬雪深、雪溫，取得了較好的模擬結(jié)果。(2)提高模型參數(shù)的準確性，聞昕等[18]利用機器學習的方法，從實測數(shù)據(jù)自動提取融雪模型參數(shù)來提高模型精度。(3)建立分布式融雪模型[19]，ABBOTT等[20]采用分布式的能量平衡融雪模型（SHE）來反應融雪在空間上的差異性；BROWN等[21]開發(fā)了分布式融雪模型UEBGRID在喜馬拉雅地區(qū)進行模擬試驗，發(fā)現(xiàn)模型對雪深有較好的表現(xiàn)。(4)結(jié)合地形與覆蓋物的參數(shù)，開發(fā)一種基于時間和空間的融雪模型，這種模型在融化過程中考慮了降雨、輻射大氣湍流等因素[22]。

近年來，關于東北黑土區(qū)積雪消融過程的研究中取得了一定進展，為融雪徑流預測提供了科學依據(jù)。楊鑫等[23]利用黑龍江鶴山農(nóng)場10年氣象資料采用線性回歸方法分析積雪深度和融雪徑流深的關系。田琳[24]基于遺傳算法對松花江白山流域構建積雪消融變化經(jīng)驗方程，并預測徑流量。整體而言，東北黑土區(qū)也多以遙感的方法和統(tǒng)計學方法對積雪深度進行反演模擬[25]，太陽輻射對積雪消融過程的影響以及積雪融化的物理過程等方面考慮較少。

1 研究區(qū)概況

本研究在吉林省梅河口市吉興小流域的水土保持監(jiān)測站完成數(shù)據(jù)的監(jiān)測與采集，吉興小流域面積約為15 km2，地理坐標為125°28′13″～125°31′34″E、42°10′26″～42°14′32″N，海拔高度范圍為393～970 m，該區(qū)地形起伏較大，總體呈南高北低的趨勢，為東北低山丘陵區(qū)。屬于典型的溫帶大陸性季風氣候，多年平均降水量約為708.8 mm，冬季積雪深度約為300 mm，降雪期集中在12 月到次年2 月份，降雪日數(shù)在6～32 d，為典型的季節(jié)性凍融區(qū)。年平均氣溫4.6 ℃，全年≥10 ℃積溫達到2 732 ℃，多年平均日照時數(shù)2 556 h。

在研究區(qū)吉興小流域站點架設PC400 型雪與土壤多參數(shù)系統(tǒng)采集土壤溫度、雪深、土水勢數(shù)據(jù)；氣象生態(tài)環(huán)境監(jiān)測儀采集太陽輻射、風速風向、土壤濕度、空氣溫度、相對濕度數(shù)據(jù)；使用Snow Fork雪特性分析儀對雪液態(tài)含水量、雪密度進行數(shù)據(jù)采集，觀測時間間隔為60 min，試驗儀器基本情況見表1。

表1 試驗儀器情況Table 1 Experimental equipment situation

2 數(shù)據(jù)與方法

2.1 模型計算原理

模型基于能量平衡計算（圖1為雪面能量平衡示意圖），通過輸入短波輻射、空氣溫度、風速（2 m）和相對濕度等狀態(tài)變量來驅(qū)動模型運行。模型的能量過程用于模擬融雪、再凍結(jié)以及積雪的熱量變化，模擬計算后可以得到雪面凈輻射通量、顯熱通量、潛熱通量、總能量和雪深。

圖1 雪面能量平衡示意圖[26]Figure 1 Schematic diagram of the energy balance of the snow surface

2.1.1 基本能量平衡方程 能量平衡模型：

式中：Qnet為模擬時段Δt內(nèi)積雪層單位面積上總能量的凈輸入量(W·m-2)；Qr為雪面凈輻射通量(W·m-2)；Qh為顯熱通量(W·m-2)；Qle為潛熱通量(W·m-2)；Qp為融雪期降水所帶來的熱通量(W·m-2)；Qg為來自土壤的地熱通量(W·m-2)。

雪面凈輻射通量（Qr）的計算：

式中：Qs為入射的短波輻射(W·m-2)；α為雪面反射率；εs為積雪比輻射率；Ta為空氣溫度(℃)；Ts為雪面溫度(℃)；εαc為有云覆蓋下的大氣比射率；σ為STEFANBOHZMAN 常數(shù)。雪面反射率是影響雪表面能量平衡的重要參數(shù)，采用Noah陸面模式的積雪反照率方案[27]：

式中：αmax為新雪反照率；t為距離降雪的天數(shù)；A和B為常數(shù)，分別取0.94 和0.82。本研究中采取對于雪面溫度的計算，Abdelghani發(fā)現(xiàn)當溫度在零下5 ℃時，雪面溫度和空氣溫度存在線性關系[28]，因此雪面溫度采用高黎明等[29]建立的線性關系式：

顯熱通量和潛熱通量采用在雪面能量平衡的運用較好的空氣動力學梯度方法[30]，該模型計算簡單且計算表現(xiàn)較好。

顯熱通量（Qh）的計算：

式中：Qh＞0時，能量從大氣傳向雪層，雪面獲得能量；Qh＜0時，能量從雪層傳向大氣，雪面損失能量。ρ為空氣密度；Cp為空氣的比熱容；ras為積雪層表層與附近的參考高度間的空氣動力學阻力[31]：

式中：χ為參考高度；μ為參考高度處的風速；Zm為動力阻抗系數(shù)；Zh為熱量和水汽阻抗系數(shù)；k為VON KARMAN常數(shù)。

潛熱通量（Qle）的計算：

式中：Qle＞0時，能量從大氣傳向雪層；Qle＜0時，能量從雪層傳向大氣；λv為升華潛熱；Rd為干空氣常數(shù)；e(Ta)為空氣水汽壓；e(Ts)為雪面水汽壓；

e(Ta)和e(Ts)的計算：

esvp(Ta)為飽和水汽壓，采用泰登公式[32]。

esvp(Ta)的計算：

降雨攜帶輸入能量(Qp)：

式中：ρw為水的密度；Cw為水的比熱容；Pr為降雨量；Qp＜0表示從雪層吸熱。在模擬計算時段內(nèi)基本都為固態(tài)降水，因此本研究不考慮降雨帶來的熱量。

在融雪期間，雪面溫度一般不超過0 ℃，且此時土層為凍土層，雪層底部的熱傳導度比較低，因此地熱通量對積雪消融的影響基本可以忽略不計。

融雪深度h計算：

式中：h為融雪深度(cm)；S為雪的溶解熱；ρs為積雪密度。

能量閉合率[33]：

式中：CR為能量閉合率。CR值越接近1，能量閉合率越好。

模型評價：

納什系數(shù)（Nash）：

平均絕對誤差（MAE）：

式中：n為觀測次數(shù)；H0為模擬值；Hp為觀測值；Have為觀測平均值；納什系數(shù)（Nash）反映模擬值與觀測值的比值與1∶1 曲線的靠近程度，納什系數(shù)（Nash）越接近1，相對誤差的絕對值越小，模擬精度就越高；平均絕對誤差（MAE）越接近0，模擬的效果就越好。

2.2 數(shù)據(jù)來源及模型參數(shù)取值

本研究選取2017年2月23日-3月10日和2018年3月6日-3月23日的積雪消融過程為分析時段。模型建立的基本數(shù)據(jù)主要為野外實測數(shù)據(jù)、雪面反照率初始值采用中國全球陸面再分析40年產(chǎn)品（CRA/Land），模型中參數(shù)取值見表2。

表2 參數(shù)取值Table 2 Parameter values

3 結(jié)果與分析

3.1 能量輻射模擬結(jié)果分析

積雪消融是熱傳導過程，而凈輻射通量是反應地表能量、動量、水分輸送過程的重要參數(shù)[34]，常用于計算積雪消融過程。通過模擬2017-2018 年融雪期能量(圖2)發(fā)現(xiàn)，凈輻射晝夜變化很大，日間為正值，最高可達175.6 W·m-2；夜間為負值，最低可達-79.4 W·m-2。經(jīng)計算，2017 年和2018 年凈輻射通量日均值占日總能量均值的變化趨勢存在差異，即2017呈逐漸減小趨勢；2018年變化趨勢為先減小后增大。2017年凈輻射通量日平均值為117.8 W·m-2，占總能量的74.7%，最高達80%；2018年日均值為122.5 W·m-2，占總能量的67.4%，最高達到93%。可見凈輻射通量是主導能量平衡的主要部分。總能量的輸入呈周期性的變化，日間凈輸入為正，11點達到最大值，兩年消融期的凈能量峰值均在150 W·m-2以上，2018年的峰值達到200 W·m-2；夜間基本上為負值，最低可達-76.65 W·m-2。

圖2 2017年和2018年凈輻射通量、總能量模擬值Figure 2 Simulated values of net radiative flux and total energy from 2017 to 2018

然而受湍流交換差異的影響，2018 年凈輻射通量對積雪消融能量的貢獻較2017 年有所下降，降幅達7.2%。湍流交換是雪表面能量交換過程的主要形式，并受到顯熱通量和潛熱通量影響。其中，顯熱通量受溫度控制，體現(xiàn)近地面空氣與雪層表面的熱交換。潛熱通量則取決于水氣壓梯度、空氣氣壓梯度，表現(xiàn)為雪的升華。顯熱通量和潛熱通量是積雪消融的主要影響因子，受天氣類型的影響較大。風速平均值均在2 m·s-1以下，2018 年風速最大值達到6.9 m·s-1；溫度平均值基本保持在0 ℃以下，2018 年溫度最高值達到16.1 ℃。可見，2018年風速、溫度的平均值和最大值以及變化幅度皆均比2017年高（表3和圖3）。

圖3 201年和2018年溫度、風速實測值；顯熱通量、潛熱通量模擬值Figure 3 Measured values of temperature,wind speed,simulated values of sensible heat flux and latent heat flux from 2017 to 2018

表3 2017年和2018年溫度、風速的平均值與最大值Table 3 Average and maximum values of temperature and wind speed from 2017 to 2018

2017 年的風速和溫度變幅較2018 年小，因此顯熱通量和潛熱通量的也比較穩(wěn)定。2017 年和2018 年顯熱通量平均值為正值25.6 W·m-2和35.6 W·m-2，潛熱通量為負值為-8.9 W·m-2和-9.0 W·m-2，計算結(jié)果較為接近。但整體表現(xiàn)為2018年要比2017年更大。顯熱通量的值為正時，空氣向雪面輸送熱量，這是因為空氣溫度高于雪面溫度，會形成穩(wěn)定的逆溫層[35]。潛熱通量為負值時，說明此時雪表面正在以升華和蒸發(fā)的形式消融，向空氣輸送熱量。2018年3月15日至23日出現(xiàn)了潛熱通量為正值的時段，這是因為此時飽和水氣壓大于空氣水汽壓，溫差和風速也較小，湍流交換量小，不利于雪的升華和蒸發(fā)，此時水汽發(fā)生凝結(jié)，表現(xiàn)為空氣向雪面輸送熱量[36]。

在冬季，地表與大氣之間存在著復雜的物質(zhì)和能量交換過程。其中一個重要的指標是能量閉合率，它反映了地表所接收到的凈輻射與地表釋放出去的感熱、潛熱和土壤熱之間是否達到平衡。理想情況下，這兩者應該相等，但實際觀測中往往存在一定程度的不平衡。圖4為積雪消融過程中步長為1 h數(shù)據(jù)計算的能量閉合率日平均值。由圖4 可知，2017-2018 年能量閉合率隨積雪深度的減小而減小，這是因為雪的反照率比裸露地表高，當雪深減少時，裸露面積增加，反照率降低，會導致能量閉合率變低，能量閉合率最大可達0.25～0.67。2018 年能量閉合率要優(yōu)于2017 年，2017 年能量閉合率低的原因可能是消融周期時間長，積雪覆蓋的天數(shù)較多，導致積雪表面吸收的太陽輻射量減少，因為積雪會反射大量的入射太陽輻射[37]，這將導致能量閉合率值降低。此外，較長的消融周期還會導致積雪表面形成冰層，增加雪表面反射率也會致使能量閉合率較低。

圖4 2017年和2018年雪深與能量閉合率關系Figure 4 Relationship between snow depth and energy closure rate from 2017 to 2018

3.2 積雪深度模擬

由能量輻射模擬結(jié)果可知，凈輻射能量是積雪消融的主要驅(qū)動力，凈輻射通量在日間為正值，因此積雪的消融也主要發(fā)生在日間。圖5 描述了2017 年3 月7日的能量變化和積雪的消融狀況，該天能量與積雪變化明顯，能夠較好反映積雪消融過程。由圖5 可知，凈輻射通量是先增大后減小的趨勢，在12 點達到最大值152 W·m-2，夜間為負值。顯熱通量和潛熱通量11 點達到了極值，分別為25.34 W·m-2和-24.44 W·m-2，但顯熱通量和潛熱通量出現(xiàn)了時間上明顯的滯后性，滯后約1 h，這是因為能量主要來自太陽輻射，熱傳輸速度不同步。在出現(xiàn)極值后至23∶00，顯熱通量先減小后增大，潛熱通量先增大后減小。積雪從8∶00 開始消融，在12∶00 雪深急劇減小，融雪速率達到峰值，17∶00 后積雪消融開始變緩，直至18∶00完全停止。整體看來，凈輻射通量、顯熱通量、潛熱通量集中在9∶00-15∶00變化，共計7 h，與雪深變化的時間相吻合。

圖5 凈輻射通量、顯熱通量、潛熱通量、雪深的日變化值Figure 5 Diurnal variation of net radiation flux,sensible heat flux,latent heat flux and snow depth

2017-2018 年融雪深度實測值與模擬值對比，以及1∶1 線圖見圖6。由圖6 可知，2017-2018 年融雪期積雪深度的模擬值和實測值總體趨勢較為一致。2017 年實測值與模擬值整體均勻分布在1∶1 線兩側(cè)，模擬結(jié)果較好（Nash=0.97；MAE=0.36 cm）；2018 年在雪深為小于15 cm 時，模擬值較高于實測值（Nash=0.83；MAE=1.46 cm），盡管2018 年模擬結(jié)果可滿足精度要求，但模擬值出現(xiàn)偏高或者偏低的情況，其主要原因歸結(jié)為：（1）本研究采用的是單層融雪模型，將融雪層作為一維的整體考慮，當雪層能量為正且輸入達到臨界狀態(tài)時，積雪才開始消融，而實際上整個雪層能量輸入還未達到臨界值時，表層可能會先消融，此時表層液態(tài)水受重力作用向下遷移，使模擬值比實際值略微偏高[38]；（2）另外，暴露在空氣中的雪表面產(chǎn)生較強的湍流交換也是影響雪深變化的一個重要原因。

圖6 2017年和2018年實測雪深與模擬雪深對比Figure 6 Comparison of measured and simulated snow depth from 2017 to 2018

4 討論與結(jié)論

在積雪日變化方面，張娟等[39]對三江源地區(qū)研究表明，積雪消融時間發(fā)生在10∶00-19∶00，歷時9 h；桑婧等[40]在額爾古納的研究則發(fā)現(xiàn)開始融雪時間在2∶00，消融時間近17 h。然而，本研究區(qū)積雪消融時間集中在凈輻射通量、顯熱通量、潛熱通量日內(nèi)集中變化時間段，融化時間在9∶00-15∶00，歷時7 h。在積雪消融天數(shù)方面，周楊等[41]在青藏高原瑪多地區(qū)研究結(jié)果表明集中積雪消融天數(shù)為12～18 d，而本研究區(qū)融雪時間在9～15 d。綜上所述，可以看出各地區(qū)之間積雪消融時段差別明顯。

在模擬積雪深度變化方面，東北黑土區(qū)以往的研究多以空氣溫度等指標研究積雪變化規(guī)律，此方法屬于統(tǒng)計學的方法，不能描述雪面能量收支特征。而本研究通過能量平衡原理建立的積雪消融模型能夠更加準確描述雪面能量收支與積雪深度變化的特征。該模型可為季節(jié)性凍融區(qū)建立融雪模型的建立提供理論依據(jù)。本模型的不足之處在于假設積雪的密度與溫度為均勻分布，忽略了雪層內(nèi)部溫度和密度的變化，因此在特定的環(huán)境下可能會產(chǎn)生誤差，在后續(xù)的研究中可以將雪層分層模擬，考慮積雪內(nèi)部的非均勻性，進而更準確地模擬積雪融化過程。

積雪消融主要是與大氣間物質(zhì)能量交換的過程，但由于試驗地觀測資料局限和過程的復雜性，凈輻射通量、顯熱通量、潛熱通量不是直接觀測得到，而是采用了前人研究經(jīng)驗對所涉及能量進行估算，這在今后工作中應當增加必要的能量觀測，進一步分析積雪與天氣系統(tǒng)之間的變化狀況。

本研究基于能量平衡原理建立單層融雪模型對融雪期雪深進行分析與模擬，得出以下結(jié)論：凈輻射通量、顯熱通量、潛熱通量是研究區(qū)融雪能量來源，融雪期間呈周期性的變化。其中凈輻射通量是主要的能量來源，2017 年和2018 年積雪消融觀測期，凈輻射通量平均值分別為117.8 W·m-2和122.5 W·m-2，占到總能量的74.7%～67.4%，而湍流交換分別占到了18.8%～25.8%。能量閉合率隨雪深的增大而變大，最大值為0.25～0.67，積雪覆蓋會使能量閉合率偏低。使用單層融雪模型模擬的積雪深度變化，2017年和2018年納什系數(shù)（Nash）分別達到0.97和0.83；平均絕對誤差（MAE）分別為0.36 cm和1.46 cm。該模型能較好地模擬融雪深度的變化。