Stefan方程在土壤凍融過程模擬中的應用

劉文惠， 謝昌衛， 劉海瑞， 龐強強， 王 武，劉廣岳， 楊雨昆， 王 銘， 張 琪

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

0 引言

多年凍土是指溫度在0 ℃或低于0 ℃至少連續存在兩年的巖土層[1］。多年凍土是冰凍圈的重要組成部分，其影響能量交換、水文過程和生態環境，進而影響全球氣候系統。活動層是指覆蓋于多年凍土之上的夏季融化冬季凍結的土層[1］。活動層的水熱動態變化過程影響著凍土區水文和生態系統的生物、物理及地球化學過程。多年凍土與大氣圈之間的相互作用主要是通過活動層中的水、熱動態變化過程而實現的。多年凍土和活動層研究很早就受到了國內外學者的廣泛關注，并加大了這方面的觀測和研究力度。多年凍土與活動層被世界氣候研 究 計 劃（World Climate Research Programme，WCRP）列入“氣候與冰凍圈計劃”（Climate and Cryosphere，CliC）的主要觀測研究內容之一。由國際凍土協會和加拿大地質調查局聯合發起的全球多年凍土監測網絡（Global Terrestrial Network for Permafrost，GTN-P）用于監測全球多年凍土地溫。環北極活動層監測網絡（Circumpolar Active Layer Monitoring Network，CALM）旨在對活動層厚度和熱狀況進行監測。活動層作為多年凍土與大氣圈進行熱、質交換的最主要場所和媒介。活動層厚度又是判斷多年凍土退化最為直觀的標志，其可以更好地反映多年凍土對氣候變化的響應狀況。隨著全球變暖加劇，多年凍土退化嚴重[2-9］，尤其是活動層厚度增厚尤為明顯[4-6］。因此，氣候變化背景下的活動層凍融過程模擬、厚度制圖和變化預測是研究凍土區生態環境、水文、工程以及碳循環的基礎，也是目前凍土學領域的研究熱點。活動層厚度可以通過鉆探、坑探等方法直接測量，也可用測溫法、地球物理勘探等間接方法來估計。但是，由于監測條件的限制和監測數據的有限，實地監測很難滿足活動層厚度空間分布的模擬需要，大尺度活動層厚度空間分布以及長時間序列活動層對氣候變化的響應還是得依靠模型來解決。

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

1 Stefan方程介紹

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

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

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

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

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

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

2 Stefan方程在多年凍土中的應用

2.1 在均質土壤中的應用

2.1.1 引入N因子

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

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

表1 融化N因子（Nt）和凍結N因子（Nf）值統計結果Table 1 The thawing N-factors（Nt）and freezing N-factors（Nf）

引入N 因子后的Stefan 方程被廣泛用于更復雜地形和下墊面的大尺度活動層厚度空間分布及變化預測模擬。龐強強等[34］模擬了青藏高原多年凍土活動層厚度的空間分布狀況。徐曉明等[35］得到青藏高原多年凍土活動層厚度平均為2.39 m，活動層厚度在羌塘盆地最小，在多年凍土區邊緣、祁連山、西昆侖山、念青唐古拉山較大。氣候變化條件下，活動層厚度呈整體增大趨勢，1981—2010 年，活動層厚度的變化量為-1.54～2.24 m，變化率為-5.90～10.13 cm·a-1，平均每年變化1.29 cm。張中瓊等[36］模擬了青藏高原活動層厚度空間分布狀況并預測了未來不同氣候情景下活動層厚度的變化情況。到2050年，A1B 情景下活動層厚度最大達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］模擬得到俄羅斯北極地區鄂畢河、葉尼塞河和勒拿河流域平均活動層厚度分別為1.80 m、1.70 m和1.66 m。

引入N 因子后的Stefan 方程的模擬精度也得到大大提高。Klene 等[37］進行庫帕河流域活動層厚度空間分布制圖時用地表溫度作為溫度驅動得到的模擬誤差為14.5%，用3 年的平均N 因子和氣溫得到的誤差為17.6%，用氣溫得到的誤差為29.2%。同時，根據不同溫度驅動估算馬銜山2010年的活動層厚度，結果顯示用地表溫度作為溫度驅動得到的活動層厚度為1.13 m，相對誤差最小（4%），最接近實測值；用氣溫得到的活動層厚度為1.01 m，相對誤差最大（13.6%），明顯偏小；而用氣溫和N 因子得到的活動層厚度為1.16 m，偏大，但是誤差要遠遠小于只用氣溫驅動得到的結果。地面溫度是土壤通過地面與大氣間熱交換特征的綜合指標，比氣溫更能體現凍土的熱狀況。而且，地表溫度對外界條件反應更為迅速，能較為準確地反映凍土的上邊界熱狀況，是眾多凍土模型最佳的選擇。因此，在完成大尺度活動層厚度空間分布模擬時，當地表溫度缺失的情況下，為了盡可能減小模擬誤差，可以考慮用N因子和氣溫作為溫度驅動。

2.1.2 引入E值

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

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

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

2.1.3 加入降雨、地形因子

山地多年凍土的形成和發展與局地微氣候和微地形（坡度和坡向）有很大的密切關系。基于此，一些學者在上式中加入了降雨和地形因子。Hinkel等[40］將阿拉斯加北方森林多年凍土區的實測活動層厚度和融化指數做最小二乘回歸擬合后得到Ste?fan方程的另一個形式：

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

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

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

2.2 分層土壤中的應用

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

2.2.1 J-L算法

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

2.2.2 分層總和法

分層總和法是通過分別計算凍融時各層土所消耗的凍融指數得到各層的凍結深度，最后各層的凍融深度之和為最大凍融深度。其推導過程見參考文獻[51］。分層總和法的研究較少，應用沒有得到推廣。

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

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

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

2.2.3 X-G算法

2013 年Xie 等[52］利用類比遞推的方法，提出了基于Stefan方程計算分層土壤凍融深度的簡易算法（X-G 算法）。X-G 算法能夠計算由任意多層不同厚度的土層組成的土壤的凍融過程，是目前唯一能將斯蒂芬方程運用到非均質土壤凍融過程的算法。它基于不同土層的物理參數將斯蒂芬方程應用到估算任意多層、任意厚度的土壤凍結/融化深度，而不是將不同層土壤參數進行平均處理，這樣避免了不必要的計算誤差。X-G 算法原理簡單，既可模擬單向凍結過程也可模擬雙向凍結過程，目前在模擬國內多年凍土凍融過程方面取得了較好的進展。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 算法發表以后，在相關領域受到了廣泛的關注，Kurylyk[58］認為X-G 算法在數學上是合理的，但不能很好解釋土壤凍融過程中溫度變化的過程，并認為N-M 算法能更好地將Stefan 方程應用到分層堆積的土壤。從非線性傅里葉熱傳導方程的斯蒂芬解入手，對Stefan 方程基本原理、X-G 算法、J-L算法和N-M 算法推導過程等進行了如下詳細分析[59］：

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

系數m 是土壤物理性質參數的函數。對上式求導數，得到：

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

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

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

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

4 存在的問題

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

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

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

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

4.2 沒有考慮降水的影響

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

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

4.3 含水量的敏感性

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

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

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

5 討論與結論

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

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

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