Stefan方程在土壤凍融過程模擬中的應(yīng)用

劉文惠， 謝昌衛(wèi)， 劉海瑞， 龐強(qiáng)強(qiáng)， 王 武，劉廣岳， 楊雨昆， 王 銘， 張 琪

（1.青海大學(xué)地質(zhì)工程系，青海西寧 810016； 2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國家重點實驗室藏北高原冰凍圈特殊環(huán)境與災(zāi)害國家野外科學(xué)觀測研究站，甘肅蘭州 730000； 3.青海大學(xué)生態(tài)環(huán)境工程學(xué)院，青海西寧 810016）

0 引言

多年凍土是指溫度在0 ℃或低于0 ℃至少連續(xù)存在兩年的巖土層[1］。多年凍土是冰凍圈的重要組成部分，其影響能量交換、水文過程和生態(tài)環(huán)境，進(jìn)而影響全球氣候系統(tǒng)。活動層是指覆蓋于多年凍土之上的夏季融化冬季凍結(jié)的土層[1］?；顒訉拥乃疅釀討B(tài)變化過程影響著凍土區(qū)水文和生態(tài)系統(tǒng)的生物、物理及地球化學(xué)過程。多年凍土與大氣圈之間的相互作用主要是通過活動層中的水、熱動態(tài)變化過程而實現(xiàn)的。多年凍土和活動層研究很早就受到了國內(nèi)外學(xué)者的廣泛關(guān)注，并加大了這方面的觀測和研究力度。多年凍土與活動層被世界氣候研 究 計 劃（World Climate Research Programme，WCRP）列入“氣候與冰凍圈計劃”（Climate and Cryosphere，CliC）的主要觀測研究內(nèi)容之一。由國際凍土協(xié)會和加拿大地質(zhì)調(diào)查局聯(lián)合發(fā)起的全球多年凍土監(jiān)測網(wǎng)絡(luò)（Global Terrestrial Network for Permafrost，GTN-P）用于監(jiān)測全球多年凍土地溫。環(huán)北極活動層監(jiān)測網(wǎng)絡(luò)（Circumpolar Active Layer Monitoring Network，CALM）旨在對活動層厚度和熱狀況進(jìn)行監(jiān)測?；顒訉幼鳛槎嗄陜鐾僚c大氣圈進(jìn)行熱、質(zhì)交換的最主要場所和媒介?；顒訉雍穸扔质桥袛喽嗄陜鐾镣嘶顬橹庇^的標(biāo)志，其可以更好地反映多年凍土對氣候變化的響應(yīng)狀況。隨著全球變暖加劇，多年凍土退化嚴(yán)重[2-9］，尤其是活動層厚度增厚尤為明顯[4-6］。因此，氣候變化背景下的活動層凍融過程模擬、厚度制圖和變化預(yù)測是研究凍土區(qū)生態(tài)環(huán)境、水文、工程以及碳循環(huán)的基礎(chǔ)，也是目前凍土學(xué)領(lǐng)域的研究熱點?；顒訉雍穸瓤梢酝ㄟ^鉆探、坑探等方法直接測量，也可用測溫法、地球物理勘探等間接方法來估計。但是，由于監(jiān)測條件的限制和監(jiān)測數(shù)據(jù)的有限，實地監(jiān)測很難滿足活動層厚度空間分布的模擬需要，大尺度活動層厚度空間分布以及長時間序列活動層對氣候變化的響應(yīng)還是得依靠模型來解決。

用于模擬活動層凍融過程和厚度的模型較多，主要分為兩大類：經(jīng)驗半經(jīng)驗公式和以求解熱傳導(dǎo)方程為基礎(chǔ)的數(shù)學(xué)物理方法[10-11］。由于數(shù)學(xué)物理方法的初始邊界條件在時空上變化較大，因此，常用于計算單點的凍結(jié)融化深度。對于空間上高度非均勻的凍土空間的變化，經(jīng)驗和半經(jīng)驗型公式是一個較為實際的選擇，尤其是半經(jīng)驗半理論的計算方案，既充分考慮了凍土的物理特性和過程，又可以用大尺度參量獲取凍土在水平和垂直空間的變化，更適宜于凍土的實際研究[11］。目前，Stefan 方程是國內(nèi)外用于計算多年凍土凍結(jié)和融化深度最常用的經(jīng)驗公式。它充分考慮了氣候條件、土壤熱屬性和水分條件，形式簡單，驅(qū)動參數(shù)少，模擬效果較好，既可以用于模擬單點的凍結(jié)融化深度，也可以較方便地模擬大尺度的活動層厚度空間分布。本文就不同修正形式的Stefan方程在國內(nèi)外多年凍土活動層凍融過程模擬中的研究進(jìn)展和應(yīng)用中存在的問題進(jìn)行了探討，旨在為今后的研究提供參考。

1 Stefan方程介紹

斯蒂芬方程是奧地利科學(xué)家Josef Stefan 在研究北極地區(qū)湖冰形成過程時提出的[12］。假設(shè)冰體內(nèi)熱量傳導(dǎo)非常迅速，并且冰體內(nèi)的溫度變化是線性的。當(dāng)冰的表面溫度低于相變溫度時，冰體下部與湖水接觸面處的溫度則等于相變溫度。在一個給定的時間內(nèi)，冰體從下面湖水中得到的熱量和其從表面排出的熱量是相等的，這一關(guān)系可以表達(dá)為：

由此得出冰體厚度隨時間的公式：

事實上，式（1）僅可以描述冰體內(nèi)簡化的熱傳導(dǎo)過程，對于常規(guī)的熱傳導(dǎo)過程，通常用一維傅力葉熱傳導(dǎo)方程描述。在給定的邊界條件和假設(shè)下，Nemman 和Stefan 先后給出了這一非齊次二階微分方程的解。當(dāng)假設(shè)凍結(jié)體內(nèi)溫度梯度為線性的條件下，Stefan 解仍然可以簡化成式（2）的形式。因此，一般認(rèn)為Stefan 方程是一維熱傳導(dǎo)方程Stefan解簡化后的公式表達(dá)。

Stefan 方程提出后的數(shù)十年里面，一般用于計算湖冰、海冰等冰體的厚度。1943年，Berggren[13］對式（2）中純冰體的熱容用土壤中冰體的熱容代替，將Stefan方程應(yīng)用到淺層土層凍結(jié)和融化過程的計算中。在理論上，冰-水相變伴隨的潛熱釋放或者吸收要遠(yuǎn)大于干土本身熱容的變化，因此，當(dāng)土層內(nèi)含水量較大時，將Stefan 方程應(yīng)用到土壤中與應(yīng)用到冰體內(nèi)造成的差別并不大。由此，式（1）轉(zhuǎn)化為式（3），即為廣泛應(yīng)用于凍土學(xué)中的Stefan方程的通用形式[14］。

式中：Z 為凍結(jié)/融化深度（m）；K 為導(dǎo)熱系數(shù)（W·m-1·℃-1）；QL為土壤水分相變引起的潛熱變化，QL=Lρ(ω ?ωu)，L 為冰的融化潛熱（3.3×105J·kg-1）；ρ為土壤的干容重（kg·m-3）；ω 為總的含水量（%）；ωu為未凍水含量（%）；DDF 和DDT分別為凍結(jié)指數(shù)和融化指數(shù)（℃·d），指地表溫度日均值連續(xù)低于0 ℃的溫度累計之和地表溫度日均值連續(xù)高于0 ℃的溫度累計之和；Tt為日平均地表溫度（℃）；T0= 0。

Stefan 方程首次將地表（或者大氣）溫度的變化與冰層（或者土層）的凍結(jié)融化過程以簡單公式的形式聯(lián)系起來，極大地簡化了土壤凍結(jié)融化深度的分析計算過程。因此，Stefan 方程在推出后被廣泛運用到寒冷地區(qū)的工程建設(shè)中，并逐漸被應(yīng)用到冰川消融、凍土變化研究中[15-18］。目前，在多年凍土活動層凍融過程的研究中，Stefan 方程仍然是應(yīng)用最廣泛的模擬計算方法。

2 Stefan方程在多年凍土中的應(yīng)用

2.1 在均質(zhì)土壤中的應(yīng)用

2.1.1 引入N因子

Stefan 方程最初被運用到土層凍融過程時，為了避免因未考慮外部熱量交換因素引起的誤差，建議采用地表以下5~10 cm 的日平均溫度作為溫度驅(qū)動[19］。由于日尺度該深度的日平均溫度獲取困難，一般用日平均地表溫度計算凍結(jié)-融化指數(shù)。相比地表溫度或者一定深度內(nèi)土層溫度，氣溫更容易觀測，實際應(yīng)用中則大多采用氣溫來計算計算凍結(jié)-融化指數(shù)。這在一定程度上造成了額外的計算誤差。針對這種狀況，Carlson[20］最早提出N 因子，對Stefan方程做了進(jìn)一步的修正：

式中：Nt為融化N 因子，是地面融化指數(shù)與氣溫融化指數(shù)之比；Nf為凍結(jié)N 因子，是地面凍結(jié)指數(shù)與氣溫凍結(jié)指數(shù)之比。從概念上講，N 因子用一個簡單的數(shù)值代表了特定下墊面地表能量通量的總和。一般情況下，融化N 因子大于凍結(jié)N 因子。不同下墊面類型的N因子不相同，但是，同一區(qū)域相同地表條件的N 因子相差不大，基本可以取固定值。從表1的統(tǒng)計結(jié)果可以發(fā)現(xiàn)，青藏高原的融化N 因子和凍結(jié)N因子普遍大于達(dá)拉斯加和加拿大等高緯度地區(qū)的融化N 因子和凍結(jié)N 因子；土壤顆粒越粗，融化N因子越大。裸地融化N 因子最大（＞1.28），其次是灌叢和苔原，河道和沼澤的融化N 因子最小（＜0.7）[31］。北極地區(qū)的N 因子值比較穩(wěn)定，其年變化幅度小于10%[32］。山地多年凍土區(qū)的N因子值年際變化比較大，尤其是冬季隨著積雪的年變化而變化較大[33］。

表1 融化N因子（Nt）和凍結(jié)N因子（Nf）值統(tǒng)計結(jié)果Table 1 The thawing N-factors（Nt）and freezing N-factors（Nf）

引入N 因子后的Stefan 方程被廣泛用于更復(fù)雜地形和下墊面的大尺度活動層厚度空間分布及變化預(yù)測模擬。龐強(qiáng)強(qiáng)等[34］模擬了青藏高原多年凍土活動層厚度的空間分布狀況。徐曉明等[35］得到青藏高原多年凍土活動層厚度平均為2.39 m，活動層厚度在羌塘盆地最小，在多年凍土區(qū)邊緣、祁連山、西昆侖山、念青唐古拉山較大。氣候變化條件下，活動層厚度呈整體增大趨勢，1981—2010 年，活動層厚度的變化量為-1.54～2.24 m，變化率為-5.90～10.13 cm·a-1，平均每年變化1.29 cm。張中瓊等[36］模擬了青藏高原活動層厚度空間分布狀況并預(yù)測了未來不同氣候情景下活動層厚度的變化情況。到2050年，A1B 情景下活動層厚度最大達(dá)10.20 m，增加約0.30~0.80 m；B1 情景下增加0.20~0.50 m；A2 情景下增加0.20~0.55 m。Klene等[37］完成了庫帕河流域活動層厚度的空間分布制圖，Shur 等[32］完成了俄羅斯多年凍土分布制圖，張廷軍等[38］模擬得到俄羅斯北極地區(qū)鄂畢河、葉尼塞河和勒拿河流域平均活動層厚度分別為1.80 m、1.70 m和1.66 m。

引入N 因子后的Stefan 方程的模擬精度也得到大大提高。Klene 等[37］進(jìn)行庫帕河流域活動層厚度空間分布制圖時用地表溫度作為溫度驅(qū)動得到的模擬誤差為14.5%，用3 年的平均N 因子和氣溫得到的誤差為17.6%，用氣溫得到的誤差為29.2%。同時，根據(jù)不同溫度驅(qū)動估算馬銜山2010年的活動層厚度，結(jié)果顯示用地表溫度作為溫度驅(qū)動得到的活動層厚度為1.13 m，相對誤差最小（4%），最接近實測值；用氣溫得到的活動層厚度為1.01 m，相對誤差最大（13.6%），明顯偏?。欢脷鉁睾蚇 因子得到的活動層厚度為1.16 m，偏大，但是誤差要遠(yuǎn)遠(yuǎn)小于只用氣溫驅(qū)動得到的結(jié)果。地面溫度是土壤通過地面與大氣間熱交換特征的綜合指標(biāo)，比氣溫更能體現(xiàn)凍土的熱狀況。而且，地表溫度對外界條件反應(yīng)更為迅速，能較為準(zhǔn)確地反映凍土的上邊界熱狀況，是眾多凍土模型最佳的選擇。因此，在完成大尺度活動層厚度空間分布模擬時，當(dāng)?shù)乇頊囟热笔У那闆r下，為了盡可能減小模擬誤差，可以考慮用N因子和氣溫作為溫度驅(qū)動。

2.1.2 引入E值

Harlan 和Nixon 認(rèn)為Stefan 方程可以表示為土壤特征與融化指數(shù)之間的線性函數(shù)[39］：

式中：E 為與土壤熱參數(shù)、含水量和地表覆蓋類型有關(guān)的比例參數(shù)，表示活動層的融化速率。已知任意點的實測活動層厚度和融化指數(shù)，可以得到E 值。已知E 值和其他土壤參數(shù)的情況下，便可得到其中的任一土壤熱參數(shù)。Hinkel等[40］得到阿拉斯加北部森林多年凍土區(qū)1992 年和1993 年的“E”值分別為1.22和1.21。假定土壤孔隙度分別為0.5和0.6時，得到導(dǎo)熱系數(shù)分別為0.126和0.151 W·m-1·℃-1。

已知E值和融化指數(shù)就可以得到活動層厚度的空間分布狀況。Peng 等[41］得到黑河流域的E值范圍為0.028~0.053，其中砂礫石E值最大，裸巖次之，荒漠的最小。并結(jié)合氣溫得到黑河流域2000—2008 的年平均凍結(jié)深度的空間分布范圍為1.0~3.5 m。Shiklomanov 等[42-43］成了阿拉斯加庫帕河1987—1999 年的年平均活動層厚度高精度制圖（50 m 分辨率），并發(fā)現(xiàn)在年變化尺度上，活動層厚度與融化指數(shù)具有高度一致性。1989 年和1998 年的融化指數(shù)最大，對應(yīng)最大的活動層厚度，1991 年和1996 年的融化指數(shù)最小，對應(yīng)最小的活動層厚度。Brown等[44］和Hinkel等[45］基于極地活動層監(jiān)測數(shù)據(jù)和氣溫融化指數(shù)做擬合分析得到，活動層厚度與融化指數(shù)具有很好的對應(yīng)關(guān)系，阿拉斯加、加拿大北部、北歐地區(qū)和俄羅斯地區(qū)的一些監(jiān)測點均顯示活動層厚度在1998年達(dá)到最大值（對應(yīng)最大融化指數(shù)），在2000年達(dá)到最小值（對應(yīng)最小融化指數(shù)）。

2.1.3 加入降雨、地形因子

山地多年凍土的形成和發(fā)展與局地微氣候和微地形（坡度和坡向）有很大的密切關(guān)系?；诖耍恍W(xué)者在上式中加入了降雨和地形因子。Hinkel等[40］將阿拉斯加北方森林多年凍土區(qū)的實測活動層厚度和融化指數(shù)做最小二乘回歸擬合后得到Ste?fan方程的另一個形式：

認(rèn)為α項表示為降雨對活動層厚度的影響。暖季一天內(nèi)不同時刻的降雨對活動層會產(chǎn)生不同的影響效果。日出前后，氣溫很低，此時形成的降水雨水溫度較低（可稱之為冷降雨事件），而形成于午后氣溫相對較高時的降水溫度較高，可稱之為暖降雨事件。其中，暖降雨雨水溫度高，降落地面后會帶入大量感熱進(jìn)入深層土壤，從而加速了活動層融化速率，使得活動層厚度增厚。但是，這個研究結(jié)果只限于阿拉斯加森林環(huán)境下泥炭層較厚的多年凍土區(qū)，其他地區(qū)沒有嘗試過。Shiklomanov 等[42］完成庫帕河流域活動層厚度高精度制圖時用降水代替含水量（該區(qū)域蒸發(fā)很小，假定降水全部用于下滲）發(fā)現(xiàn)活動層融化速率與前一年融化期末的降水有很大的相關(guān)性，1996 年融化期末降水最少，對應(yīng)1997年的活動層融化速率最大。

Nelson 等[46］在上式中加入了地形因素，完成了阿拉斯加庫帕河流域更詳盡的活動層厚度分布制圖（1 km 分辨率），得到北坡（陰面）的活動層厚度在減?。ㄘ?fù)值），南坡（陽面）的活動層厚度在增大。

式中：r 為地形因子引起的輻射因素，是個無量綱參數(shù)；Rs為斜坡上的潛在太陽輻射，與緯度、坡度和坡向相關(guān)；Rh為水平面上的潛在太陽輻射。

2.2 分層土壤中的應(yīng)用

土壤不同于冰體那樣是由均質(zhì)成分構(gòu)成，在不同深度巖土成分、結(jié)構(gòu)以及水分條件等通常有較大的差異。將本來用于計算均質(zhì)冰體凍融深度的Ste?fan方程應(yīng)用到非均質(zhì)的巖土中，必然會引起一定的計算誤差。因此，不同時期的不同學(xué)者提出了一些將Stefan方程用于計算分層堆積土壤凍融深度的算法，如被廣泛應(yīng)用的J-L算法[47-49］、N-M算法[50］、分層總和法[51］和X-G算法[52］。

2.2.1 J-L算法

J-L 算法（也被稱為St Paul 方程）由Jumikis 和Lunardini自1950年代提出，在工程計算中得到了廣泛的應(yīng)用，并被應(yīng)用到許多大型的數(shù)學(xué)模型中。原理是通過計算凍融到第n層所需要的凍融指數(shù)來得到最大凍融深度，其推導(dǎo)過程見參考文獻(xiàn)[47-48］。Woo 等[53-54］利用J-L 算法模擬了加拿大西北部馬更些河流域的活動層厚度并分析了有機(jī)質(zhì)層的影響。當(dāng)泥炭層厚度為0.2 m 時，A2 情景下，苔原活動層厚度由0.68 m 增大到2100 年的0.96 m；森林活動層由1.35 m 增大到2100 年的1.92 m。當(dāng)泥炭層厚度增厚到1 m 時，A2 情境下苔原的活動層厚度由0.39 m 增大到2100 年的0.50 m；森林活動層由0.65 m 增大到2100 年的0.75 m。同時，分析了土壤質(zhì)地、含水量和溫度等對該區(qū)域不連續(xù)多年凍土活動層厚度的影響[55］。Woo等[56］于2004將J-L算法做了改進(jìn)后用于模擬雙向的凍結(jié)融化過程，得到北美自北向南多年凍土區(qū)到季節(jié)凍土區(qū)6種下墊面的凍結(jié)融化過程。對比單向的凍結(jié)融化過程，雙向的凍結(jié)融化過程大大提高了模擬精度，誤差范圍為0.16~0.58，而單向模擬的結(jié)果誤差高達(dá)0.53~1.34。

2.2.2 分層總和法

分層總和法是通過分別計算凍融時各層土所消耗的凍融指數(shù)得到各層的凍結(jié)深度，最后各層的凍融深度之和為最大凍融深度。其推導(dǎo)過程見參考文獻(xiàn)[51］。分層總和法的研究較少，應(yīng)用沒有得到推廣。

Stefan 方程通用形式[式（3）］可以改寫成如下形式：

從式（8）中可以看出，凍融指數(shù)可以看成是凍融深度的冪函數(shù)，即凍融指數(shù)與凍融深度的平方成正比，對于厚度為z = z1+ z2的均質(zhì)土壤，如下關(guān)系式成立：

從式（9）可以看出，即使對均質(zhì)土壤，要凍結(jié)/融化厚度為z 的土壤，其所需要消耗的凍結(jié)/融化指數(shù)要大于分別凍結(jié)/融化厚度為z1和z2的土壤所需要的凍融指數(shù)之和，因此，簡單地利用分層總和法估算非均質(zhì)土壤凍結(jié)/融化深度是錯誤的。

2.2.3 X-G算法

2013 年Xie 等[52］利用類比遞推的方法，提出了基于Stefan方程計算分層土壤凍融深度的簡易算法（X-G 算法）。X-G 算法能夠計算由任意多層不同厚度的土層組成的土壤的凍融過程，是目前唯一能將斯蒂芬方程運用到非均質(zhì)土壤凍融過程的算法。它基于不同土層的物理參數(shù)將斯蒂芬方程應(yīng)用到估算任意多層、任意厚度的土壤凍結(jié)/融化深度，而不是將不同層土壤參數(shù)進(jìn)行平均處理，這樣避免了不必要的計算誤差。X-G 算法原理簡單，既可模擬單向凍結(jié)過程也可模擬雙向凍結(jié)過程，目前在模擬國內(nèi)多年凍土凍融過程方面取得了較好的進(jìn)展。Xie 等[52］模擬得到馬銜山、兩道河、昆龍山埡口和西大灘的活動層厚度分別為1.12 m、1.23 m、1.72 m和1.86 m，與實測值的相對誤差范圍為4.2%~17%。劉文惠等[57］利用該算法模擬了馬銜山2010—2013 年的活動層厚度，模擬值分別為1.06 m、99 m、1.16 m 和1.04 m，與實測值很接近，相對誤差范圍為1%~9%。

X-G 算法發(fā)表以后，在相關(guān)領(lǐng)域受到了廣泛的關(guān)注，Kurylyk[58］認(rèn)為X-G 算法在數(shù)學(xué)上是合理的，但不能很好解釋土壤凍融過程中溫度變化的過程，并認(rèn)為N-M 算法能更好地將Stefan 方程應(yīng)用到分層堆積的土壤。從非線性傅里葉熱傳導(dǎo)方程的斯蒂芬解入手，對Stefan 方程基本原理、X-G 算法、J-L算法和N-M 算法推導(dǎo)過程等進(jìn)行了如下詳細(xì)分析[59］：

分層堆積的土壤中凍結(jié)融化深度作為時間的分段函數(shù)，其函數(shù)形式為：

系數(shù)m 是土壤物理性質(zhì)參數(shù)的函數(shù)。對上式求導(dǎo)數(shù)，得到：

可見其導(dǎo)數(shù)形式也是不連續(xù)的。分層堆積土壤中的凍融過程是時間的分段函數(shù)，并不是連續(xù)函數(shù)，從而最終證明基于連續(xù)函數(shù)推導(dǎo)得到，即建立在連續(xù)微分形式上的J-L 算法和N-M 算法是錯誤的，X-G 算法是目前唯一能將斯蒂芬方程運用到分層堆積土壤的算法。

3 與陸面和水文模型的耦合

目前的陸面過程模型和水文模型的研究對均勻、覆蓋稠密的下墊面條件研究比較成熟，且較多，但忽略了凍土和積雪等復(fù)雜因素或僅做了十分簡單的概化。考慮各圈層相互作用，發(fā)展多圈層綜合陸面過程面模型和水文模型成為一種必要。多年凍土對地氣之間能量交換、水循環(huán)有很重要的影響，活動層凍融過程的準(zhǔn)確模擬可以很好地研究多年凍土區(qū)地-氣-能-水的交換過程。因此，在陸面過程面模型和水文模型中耦合凍土模型，反映氣候變化背景下土壤凍融過程中水熱過程遷移等對陸地-大氣水熱交換過程和寒區(qū)環(huán)境水文過程的影響是目前大氣和水文學(xué)者的關(guān)注重點。許多研究者也在陸面過程模型框架下發(fā)展了凍土參數(shù)化方案，發(fā)現(xiàn)在陸面水文過程模型中考慮凍土作用能顯著地增強(qiáng)模型模擬能力[60-63］。不少研究人員嘗試在分布式水文模型中添加凍土過程，以適應(yīng)寒區(qū)環(huán)境水文過程模擬。研究表明，包含凍土凍融模塊的模型能夠成功模擬由凍土融化引起的春季洪水過程，而運用不具備凍土模擬功能的分布水文模型在凍土明顯的寒區(qū)流域進(jìn)行模擬時，不能準(zhǔn)確捕捉融雪徑流過程，由于缺失凍土模塊而嚴(yán)重低估融雪期間的徑流洪峰，但暖期徑流量又明顯偏高[64-67］。

Stefan 方程由于其原理簡單、計算簡便、所需要參數(shù)易于獲得，在一些陸面過程和水文模型中得到應(yīng)用。多年凍土與氣候模式的耦合有兩種方法，一種是將凍結(jié)融化過程融合到氣候模式中；一種是基于氣候模式的溫度數(shù)據(jù)作為凍結(jié)融化過程的溫度驅(qū)動，這種方法目前應(yīng)用較多，比較成熟。Stendel等[68］基 于OAGCMs（ocean-atmosphere general cir?culation models）和Stefan 方程得到當(dāng)前氣候條件下（1961—1990 年）北半球多年凍土活動層厚度范圍為0~1.2 m；IPCCA2情景模式下2071—2100年活動層厚度將增加1.1~1.6 倍，其中俄羅斯北極地區(qū)活動層厚度將增加30%~40%，增加幅度最大的是西伯利亞東北部和加拿大西部，其次是中國和蒙古北部。Li等[69］將Stefan 方程的常用表達(dá)式耦合到陸面模式SiB2 中發(fā)現(xiàn)，加入凍土參數(shù)化后的模擬精度大大提高，模擬結(jié)果誤差不到9%。其中，模擬得到的表層和底部土壤濕度的絕對誤差分別為0.020和0.013，遠(yuǎn)遠(yuǎn)小于沒有加入凍土參數(shù)化的模擬結(jié)果。Yi等[70］考慮泥炭層、未凍水和冠層儲熱后將雙向的J-L算法耦合到陸面模式CLM3中，對比沒有加入J-L算法的模擬結(jié)果，改進(jìn)后的CLM3 消除了凍結(jié)融化鋒面在3、4月份的較大波動以及在1月和3月出現(xiàn)稽延期后較大的跳躍下降現(xiàn)象；并且在模擬土壤溫度、土壤含水量和雪深方面的精度大大提高。Fox[71］將J-L算法融合到水文模型中去，基于土壤水熱相互作用，分析了凍結(jié)融化過程對土壤水量平衡要素的影響，并通過敏感性實驗分析得出凍融過程對土壤徑流變化和預(yù)測極端徑流事件具有重要的潛在影響。與最新版本的大尺度寒區(qū)水文模型（cold region hydro?logical model，CRHM）凍土模塊相比，盡管耦合Ste?fan 方程后的水文模型在描述凍土活動層凍融過程對水文過程的影響方面有較大優(yōu)勢，但不能準(zhǔn)確描述凍土融水過程、水分遷移入滲過程、凍土坡面匯流以及產(chǎn)匯流過程等，更不能定量化地表徑流和壤中流的量及凍土融水對徑流的貢獻(xiàn)值。

4 存在的問題

4.1 模擬誤差分析及可能的原因

土壤中的凍融過程是非常復(fù)雜的，許多因素起著關(guān)鍵的影響作用。Stefan 方程以一維熱傳導(dǎo)方程為理論基礎(chǔ)，假設(shè)地面吸收的熱量全部用于多年凍土的融化且溫度變化是線性得出的，這顯然與實際凍融過程不一致。因此，利用Stefan 方程模擬得到的凍結(jié)融化深度與實際凍融深度之間通常有或大或小的誤差。一般認(rèn)為，Stefan 方程忽略了感熱變化，沒有考慮外部熱量交換和凍結(jié)巖層與下覆融土層的熱量交換，最終導(dǎo)致計算結(jié)果偏大。但是，在實際的應(yīng)用過程中卻出現(xiàn)了模擬值偏小的情況，如蒙古北部地區(qū)和青藏高原的個別點（表2）。針對這種情況，可能的解釋是同一土壤剖面內(nèi)不同深度的土壤含水量、導(dǎo)熱率和容重等參數(shù)存在較大的差異，將這些參數(shù)取平均值或者依照一層同性選取參數(shù)有可能拉低了整個參數(shù)取值，造成計算結(jié)果有一些偏小。

表2 基于Stefan方程的活動層厚度估算值小于實測值的結(jié)果統(tǒng)計Table 2 The relative errors between calculated active layer thickness and measured active layer thickness

凍融指數(shù)和導(dǎo)熱系數(shù)作為Stefan 方程兩個重要輸入?yún)?shù)，其計算的準(zhǔn)確性和獲取方法差異影響模擬結(jié)果，進(jìn)而引起一定的模擬誤差。凍融指數(shù)在一定程度上可指示凍融作用的深度、強(qiáng)度及持續(xù)時間，其變化深刻影響凍融作用下形成的冰緣地貌和寒區(qū)地質(zhì)環(huán)境[72］。凍融指數(shù)計算方法包括經(jīng)典日計算方法、月平均方法和年振幅方法[23］。凍融指數(shù)計算時針對溫度數(shù)據(jù)缺測時間長短采用不同的插補(bǔ)方式，缺測1天，選擇其前后各一天數(shù)值取平均值插補(bǔ)；缺測2天，缺測第一天選擇該日前兩天的數(shù)值取平均，缺測第二天選擇該日后兩天的數(shù)值取平均；缺測超過2天但不超過一個月，選擇前后各一年該月數(shù)據(jù)進(jìn)行插值補(bǔ)充[73］。此外，不同學(xué)者針對不同研究區(qū)域的凍結(jié)期和融化期開始結(jié)束時間的界定方法不同。凍融指數(shù)不同計算方法和溫度數(shù)據(jù)缺測不同插補(bǔ)方式得到的凍融指數(shù)有一定的差異，進(jìn)而影響凍融指數(shù)的計算準(zhǔn)確性。凍土導(dǎo)熱系數(shù)反映了地層凍結(jié)過程中包括土中相變潛熱影響的綜合導(dǎo)熱能力，直接影響凍結(jié)冷量在地層中的傳遞過程。多年凍土導(dǎo)熱系數(shù)基于野外測試方法和計算模型得到?，F(xiàn)存的適用于凍土區(qū)的導(dǎo)熱系數(shù)計算模型多以一種或幾種土壤條件為前提，或者多考慮局地因素影響。同時，多年凍土區(qū)土壤受凍融循環(huán)影響較大，多年凍土內(nèi)部水熱傳輸過程復(fù)雜，模型沒有考慮未凍水含量、土骨架組成及凍土結(jié)構(gòu)等對凍土導(dǎo)熱系數(shù)影響，使得模型得到的導(dǎo)熱系數(shù)精度不高、適用性具有局限性[74-75］。

4.2 沒有考慮降水的影響

一般來講，降雨增加，活動層含水量增加，土壤凍結(jié)需要釋放更大的潛熱，使得土壤凍結(jié)過程中的溫度降低大大減小；同時，降雨減少，秋季凍結(jié)的冰含量減少，來年用于融化冰消耗的熱量減少，更多的熱量用于加熱活動層，使得活動層融化速率增大。然而，Subin 等[76］卻認(rèn)為降雨會增大地表感熱傳遞和土體融化潛熱量，使得活動層和多年凍土溫度升高。

降雨對活動層水熱的影響比較復(fù)雜。不同區(qū)域的影響不同；同一區(qū)域不同降雨強(qiáng)度和降雨時長、一年中暖季降雨和冷季降雨以及一天中不同時刻的暖降雨和冷降雨均對活動層水熱有影響且影響機(jī)理不同。國外有關(guān)這方面的研究僅在阿拉斯加北方森林有過。青藏高原降水較多，主要集中在5—9 月，東部降水多，西部降水少[77］。降水的這種時空分布不均勻性對活動層水熱影響差異較大。目前，除了北麓河地區(qū)外，整個青藏高原上有關(guān)這方面的研究至今是空白。張明禮等[78-79］在北麓河的研究認(rèn)為夏季短期、高頻次降雨有減少地表凈輻射、增加地表蒸發(fā)潛熱、降低土壤表層溫度的作用。李德生等[80］在北麓河地區(qū)發(fā)現(xiàn)暖季的高頻率、小雨量降雨對活動層具有顯著的冷卻效果，且凌晨2 點左右的降雨產(chǎn)生能量變量很小，而14點左右的降雨產(chǎn)生的能量變量最大。Wen 等[81］通過凍土監(jiān)測發(fā)現(xiàn)，北麓河夏季降雨入滲對流作用降低了地表溫度梯度、減小土壤熱通量，降雨增加減緩凍土的退化。因此，可以嘗試在Stefan 方程中加入降雨后定量分析降雨尤其是極端降雨事件對活動層融化速率的影響。

4.3 含水量的敏感性

土壤含水量和未凍水含量作為Stefan 方程的重要輸入項，考慮了冰水二相轉(zhuǎn)換釋放的潛熱對融化深度的影響。當(dāng)其他條件恒定，土壤含水量發(fā)生變化時，季節(jié)凍結(jié)和季節(jié)融化深度相應(yīng)變化。Woo等[56］也發(fā)現(xiàn)相比干容重、有機(jī)質(zhì)和礦物質(zhì)含量，J-L算法對地表溫度和土壤含水量更敏感，得到當(dāng)?shù)乇頊囟壬? ℃時，最大凍結(jié)深度減小了0.22 m，當(dāng)土壤重量含水量分別增加50%和150%時最大凍結(jié)深度減小了0.06 m和0.24 m。但是，在少數(shù)研究中往往沒有考慮土壤含水量和未凍水含量的影響。Ro?manovsky 等[18］忽略含水量后模擬阿拉斯加北極海岸平原自北向南的West Dock、Franklin 和Franklin Bluffs 三個點1987—1992 年的活動層厚度，得到的最大誤差分別高達(dá)26%、3%和71%，存在明顯高估的情況。Stefan 方程對含水量具有高度依賴性。在水分較充足的區(qū)域模擬得到的誤差較小，在相對干旱的區(qū)域得到的誤差較大。正如圖1 所示，模擬值與實測值的相對誤差隨含水量的增加而降低。這符合Stefan 方程理論上的定義，Stefan 方程只適用于計算湖水（冰）的凍結(jié)（融化）厚度，如果用以計算土層的凍結(jié)或融化深度，則其計算結(jié)果的誤差隨土層含水（冰）量的減少而增大。

圖1 青藏高原活動層厚度模擬值與實測值的相對誤差與土壤含水量的關(guān)系Fig.1 Relationship between relative error and soil water content of permafrost on the Qinghai-Tibet Plateau

在多年凍土區(qū)活動層底部及上限附近最冷時期存在著較高的未凍水含量，這部分未凍水并沒有參與實際的凍結(jié)融化過程。在計算中應(yīng)該除掉這部分未凍水含量。Wilhelm 等[82］估算了南極半島西部阿姆斯勒島活動層厚度，由于該區(qū)域含水量比較少（重量含水量0~8.4%），忽略了含水量后得到的活動層厚度明顯偏大。Uxa[83］對Wilhelm 等的模擬結(jié)果做了糾正，發(fā)現(xiàn)該區(qū)域多年凍土中還存在2%的未凍水，假定這部分未凍水參與凍融過程，得到的活動層厚度更接近實測值?；赬-G 算法不考慮未凍水得到馬銜山2010—2013 四年的活動層厚度均小于實測厚度，但實際上馬銜山多年凍土區(qū)活動層下部最冷時期始終保持0.1 m3·m-3的體積未凍水。將這部分體積未凍水轉(zhuǎn)換成重量未凍水并假定其不參與融化過程，得到的活動層厚度為118 cm；而當(dāng)以總的含水量為輸入項開展模擬時，由于高估了凍融過程中冰水二相轉(zhuǎn)化消耗的潛熱，導(dǎo)致模擬深度減小，模擬的活動層厚度為104 cm，明顯小于未凍水不參與的情況。不同土壤質(zhì)地的未凍水含量是不相同的，同一土壤類型的未凍水含量隨溫度而變化。一般情況下，由于缺乏土壤水分的觀測資料，尤其是未凍水含量的監(jiān)測更少，很難滿足凍土空間分布上的確定，在計算過程中未凍水含量常取固定值，必定使得模擬結(jié)果存在較大的誤差和不準(zhǔn)確性。

5 討論與結(jié)論

針對Stefan 方程進(jìn)行的多種改進(jìn)措施雖然增強(qiáng)了方程對影響凍融過程的因素的包容性，但改進(jìn)后的方程往往限制更多，不同地區(qū)應(yīng)用時要引用更多的經(jīng)驗性因子。事實上，由于影響凍融過程因素很多，復(fù)雜的方程并不一定會取得更好的模擬效果。如俄羅斯凍土學(xué)家Kudryavtsev 于1974 年提出的Kudryavtsev 方程[84］是在凍土學(xué)中與Stefan 方程并列的計算活動層凍融深度的重要方法。該方程既可以模擬活動層厚度，也可以模擬多年凍土頂板溫度。計算中不僅考慮了潛熱，而且考慮了積雪、土壤的熱傳導(dǎo)和熱容量效應(yīng)，輸入?yún)?shù)較多，更接近于真實情況。過去20 多年來，Romanovsky 方程被廣泛用于北極圈活動層厚度和多年凍土頂板溫度的模擬、不同氣候情景模式下活動層厚度的預(yù)測[85-86］、活動層加厚引發(fā)的潛在危險評估[87-88］，以及全球變暖背景下俄羅斯北極濕地溫室氣體的排放評估[89］。在青藏高原地區(qū)也有多次應(yīng)用，如Pang等[90］模擬了青藏高原活動層厚度的空間分布；Luo等[91］得到黃河源區(qū)多年凍土頂板溫度和活動層厚度的空間分布。從現(xiàn)有的文獻(xiàn)來看，一般認(rèn)為Ro?manovsky 方程的模擬精度要好于Stefan 方程。如Romanovsky 等[18］、Shiklomanov 等[92］在模擬阿拉斯加北坡多年凍土活動層厚度時得到Kudryavtsev 模型模擬精度遠(yuǎn)遠(yuǎn)大于Stefan方程。但也有的研究認(rèn)為Kudryavtsev 方程所需要參數(shù)太多從而導(dǎo)致了模擬誤差。如Yin 等[93］發(fā)現(xiàn)Steafn 方程在模擬五道梁多年凍土活動層厚度時的誤差（＜7%）要小于Kudry?avtsev 模型（2.8%~27.4%）。從理論上講，Kudry?avtsev 方程可以看作是Stefan 方程的一種全面改進(jìn)，Romanovsky 等[18］對Kudryavtsev 方程與Stefan方程之間的轉(zhuǎn)換方法進(jìn)行了較為詳細(xì)的介紹，本文不再贅述。

實際上，多年凍土的凍融結(jié)過程是土壤水分相變的過程，當(dāng)水分由冰相轉(zhuǎn)化為水相時，相對應(yīng)的土壤熱容量變小，熱傳導(dǎo)加快，從而使得凍結(jié)鋒面的位置發(fā)生變化。在較小的日時間尺度上，含水量的影響也許可以忽略，但在時間較長的年尺度上具有至關(guān)重要的作用。因此，在一些含水量較高的多年凍土區(qū)，不僅要考慮總含水量，還要考慮到未凍水含量的年變化。未凍水不僅直接影響Stefan方程的模擬結(jié)果，而且通過影響導(dǎo)熱系數(shù)進(jìn)而再次影響模擬結(jié)果。因此，Stefan 方程在凍土模擬中未來考慮和改進(jìn)的關(guān)鍵是未凍水，如何通過大尺度參量來確定土壤未凍水含量的變化是一個基本而又很重要的問題，這也是目前Stefan 方程最急需考慮和解決的重中之重。同時，多年凍土作為冰凍圈很重要的主體之一，其與其他圈層的相互作用越來越受到重視。一些學(xué)者試圖在寒區(qū)陸面模式和水文模型中考慮加入多年凍土的影響，但是僅僅將多年凍土的影響作為其中一個單獨的子模塊。如何將多年凍土的諸多參數(shù)真正融入陸面模型和水文模型中是目前亟待解決的問題。

Stefan 方程參數(shù)少、形式簡單、模擬效果可靠，是活動層厚度計算中運用最廣泛、使用最方便的公式之一。但在國內(nèi)的應(yīng)用相對較少，現(xiàn)有的研究大多只是將此公式簡單應(yīng)用于模擬多年凍土活動層厚度的空間分布狀況。隨著青藏高原地區(qū)多年凍土變化研究的深入，Stefan 方程的應(yīng)用必將日趨廣泛。本文簡要介紹了Stefan方程的推導(dǎo)背景和在凍土研究中的一些改進(jìn)，希望起到拋磚引玉的作用，未來相關(guān)學(xué)者可在此基礎(chǔ)上，對這一歷史悠久的凍融過程模擬計算方法開展更深入地研究，使其在青藏高原多年凍土研究中得到更廣泛地應(yīng)用。