多年凍土區高寒草甸蒸散發模型適用性評價
——以青藏高原風火山地區為例

肖錦旺， 王根緒， 胡兆永， 郭林茂， 李金龍

（四川大學 水利水電學院 山區河流保護與治理全國重點實驗室，四川 成都 610065）

0 引言

青藏高原被稱為地球的“第三極”［1］，擁有世界上海拔最高，分布最廣的凍土層，是對氣候變化最敏感的地區之一［2］。在全球變暖背景下，青藏高原氣溫顯著升高，其速率約為北半球同緯度地區的兩倍［3］。氣溫升高將改變地表的潛熱和感熱通量，并在一定程度上對區域的水文循環產生影響。蒸散發（ET）是水文循環的重要過程之一，是理解土壤-植被-氣候相互作用的主要變量［4］。活動層厚度的增加、土壤溫度的升高以及植被覆蓋率［5］的增加將增大土壤蒸散發，進而影響區域水熱平衡［6］。由于蒸散發量大小取決于氣象因素（太陽輻射、氣溫以及濕度等）、下墊面植被以及土壤類型［7］，直接測得蒸散發十分困難。目前對于青藏高原蒸散發監測主要集中在季節凍土區［8］，然而對于廣泛分布的多年凍土區蒸散發監測極其有限［9］，Ge 等［10］指出與季節凍土相比，多年凍土活動層的凍融循環對地表能量分布影響更大。因此，開展對青藏高原典型多年凍土區不同凍融階段蒸散發模擬研究十分有必要。

蒸散發模擬研究方法主要分為以下3 類：輻射法（Priestley and Taylor［11］，PT 模型等）、綜合法（Penman-Monteit［12］，PM 模型等）、互補理論法（Advection-aridity［13］，AA 模型、Generalized Nonlinear Complementary Principle［14］，CR 模型等）。上述三種方法在結構、復雜程度上以及所需的數據各不相同，對于不同下墊面的適應性有所差異。綜合法所需的參數最多，包括空氣動力學阻力、冠層邊界層阻抗、葉面積指數和土壤信息［12］，而輻射法和互補理論法只需常規氣象數據。綜合法相比于其他方法，其對不同的下墊面適應性較強，但數據的獲取難度大。聯合國糧食和農業組織已將綜合法FAO-Penman-Monteit（FAO-PM）模型作為計算蒸散發的標準模型［15-16］。輻射法PT 模型最開始應用于在最小平流的條件計算濕潤表面的蒸散發［17-19］，是PM 模型的簡化形式。αPT為Priestley-Taylor 校準系數，該系數表示實際蒸散發速率與平衡蒸散發速率的比值［17］，其值的大小會對模型性能產生直接影響，目前已有部分研究根據實際觀測數據對青藏高原αPT值重新標定［20-21］。互補理論法包括線性互補理論法和非線性互補理論法，1979 年Bouchet 和Stricker［13］首次提出平流-干旱互補理論AA 模型（線性互補理論法），但是，Crago 等［22］和Qualls 等［23］指出在極端情況下AA 模型模擬的結果會發生偏差，極端濕潤情況下會高估蒸散發，極端干燥情況下會低估蒸散發。2015 年Brutsaert 提出了廣義非線性互補理論模型CR 模型，Wang 等［20］在原始CR 模型的基礎上，將土壤凍融過程中冰水相變所消耗的潛熱考慮在內，提出CRice，更適用于凍土地區的實際蒸散發模擬。已有研究使用上述3類模型分別對青藏高原典型多年凍土區蒸散發進行研究［24-26］，但由于不同方法的物理基礎以及需要的數據支撐條件不同，其模擬計算結果也不盡相同，存在一定的源于方法本身的不確定性，因此，迫切需要進行不同方法的綜合比較的研究，以明確不同方法的適用性和實際蒸散發的模擬能力。

為此，本文以青藏高原多年凍土區高寒草甸為研究對象，基于渦度相關法獲得實際蒸散發數據，評價FAO-PM 模型、PT 模型、AA 模型以及CRice模型在該地區不同時間尺度蒸散發模擬的適用性，明確青藏高原多年凍土區高寒草甸蒸散發模擬的最優方案，為青藏高原多年凍土區水資源開發與利用提供一定理論基礎。

1 材料和方法

1.1 研究區域概況

研究區位于青藏高原腹地北麓河一級支流左冒西孔曲風火山流域內（93°02′～92°50′ E、34°40′～34°47′ N），該區域是典型的多年凍土區（圖1）。2014 年在該區域的高寒草甸安裝一套渦度通量系統（主要包括紅外氣體分析儀（Li-7500A，LiCor，美國）和三維超聲風速儀（Windmaster Pro，Gill，英國）及氣象觀測系統。研究區域海拔4 603～5 403 m，年平均氣溫為-5.2 ℃，年平均降水為290.9 mm，但年平均蒸發量高達1 316.9 mm，年平均相對濕度為57%，多年凍土厚度為50～120 m，活動層厚度變化范圍為0.8～2.5 m［27］。研究區域植被以高寒草甸為主，高寒草甸植被群落以高山嵩草、矮嵩草和線葉草為主［27］。

圖1 研究區概況Fig. 1 Overview of the study area： the location of fenghuoshan （a）； the temperature and humidity probe and the flux tower （b）；the location of the temperature and humidity observation point and the flux tower （c）

1.2 數據來源與處理

研究區架有自動氣象站（92° 53.34′ E、34°43.07′ N，海拔4 809 m），其主要監測氣溫（℃）、相對濕度（%）、2 m 處風速（m·s-1）、凈輻射（W·m-2）與土壤熱通量（W·m-2）。在氣象站旁安裝高2 m 的渦度協方差系統監測潛熱（W·m-2）和感染通量（W·m-2），該系統由三維超聲波風速儀和紅外氣體分析組成。土壤溫度由溫濕度探頭（92°53.78′ E、34°42.87′ N，海拔4 822 m）監測獲得，其中土壤溫度測量精度±1 ℃，分辨率為0.1 ℃，測量范圍-40～60 ℃，適用于多年凍土區溫度觀測；土壤水分測量精度為±0.03 m3·m-3，分辨率為0.000 8 m3·m-3。

1.2.1 土壤溫度數據處理

分別在土壤5 cm、20 cm、50 cm、160 cm 深度布設溫濕度探頭，溫濕度探頭每4小時監測一次數據，剔除較大離群值后將4小時土壤溫度數據轉換成日平均數據。本文根據2019—2020 年5～120 cm 土壤溫度日變化（圖2），將凍融過程劃分為4 個時期：融化期（4 月20 日—5 月27 日，Tsoi5＞0.1 ℃）、完全融化期（5 月28 日—10 月17 日，Tsoil120＞0.1 ℃）、凍結期（10 月18 日—12 月8 日，Tsoil5＜-0.1 ℃）以及完全凍結期（12月8日—4月19日，Tsoil120＜-0.1 ℃）。

圖2 2019—2020年5～120 cm土壤溫度等值線Fig. 2 The 5～120 cm soil temperature contour from 2019 to 2020

1.2.2 渦度相關通量數據處理

使用EddyPro 軟件（LI-COR， USA）對渦度相關數據進行預處理，包括尖波計數和去除［28］，H2O/CO2相對于垂直風分量的滯后修正、平面擬合坐標旋轉［29］和WPL 密度波動修正［30］。過濾夜間負湍流通量以及較大的離群數據后（剔除的數據占總數據的17%），將半小時的渦度數據轉換成日平均數據。能量閉合度是評價渦度數據的重要指標，通過對研究時段內風火山觀測點能量通量數據進行閉合度分析，能量閉合度為0.75，Zhu 等［31］對中國通量網8 個站點的研究發現能量閉合度在0.58 至1.00 之間變化，通常情況下能量閉合度的可接受范圍在0.7～0.9，因此說明本研究中使用的觀測數據是可靠的。

1.2.3 氣象數據處理

氣象數據（氣溫、風速、濕度）每30 min 記錄一次，剔除較大離群值后將半小時數據轉換成日平均數據，同時篩選出每日最高氣溫和最低氣溫并將其單獨列出。

1.3 蒸散發模型

本研究中選擇基于綜合法的FAO-PM 模型、基于輻射法的PT模型，以及基于互補理論法的線性模型AA 和改進的非線性關系模型CRice，模擬青藏高原典型多年凍土區風火山2016年4月—2017年4月以及2019 年4 月—2020 年4 月蒸散發。4 種模型的輸入參數見表1。

表1 模型驅動變量Table 1 Model driving variables

1.3.1 FAO-Penman-Monteith模型

FAO-Penman-Monteith模型是1998年由聯合國糧農組織（FAO）提出的計算潛在蒸散發的標準模型［15］，考慮了空氣動力學阻力、冠層邊界層阻抗、葉面積指數和土壤水分狀況等因素對蒸散發的影響。FAO-Penman-Monteith方法的公式如下：

式中：ETFPM為潛在蒸散發量（mm）；Rn為地面凈輻射（MJ·m-3·d-1）；G為土壤熱通量（MJ·m-3·d-1）；γ為濕度計常數（kPa·℃-1）；U2為2 m高處的風速（m·s-1）；es為飽和水汽壓（kPa）；ea為實際水汽壓（kPa）；Tmean為2 m高度的平均氣溫（℃）；Δ為飽和水汽壓溫度曲線上的斜率（kPa·℃-1）。

1.3.2 Priestley and Taylor模型

Priestley and Taylor 模型［11］是Penman-Monteith模型的簡化方法，只需要常規的氣象數據就能計算蒸散發，計算公式如下：

式中：αPT為校準值常數，使用2015 年監測數據進行αPT標定，標定結果為1.18；Qne是以蒸發單位表示的可用能量（mm），Qne= (Rn-G)λ，λ為汽化潛熱（kJ·kg-1）。

1.3.3 AA模型

基于實際蒸散發與潛在蒸散發之間的反饋，Bouchet 和Stricker 提出平流-干旱模型［13］（Advection-aridity，AA），該模型認為處于潮濕環境下（水分充足），蒸散發接近潛在值；干燥環境下，用于蒸發的能量變成產生感熱通量。AA模型的公式如下：

式中：ρ為水的密度（kg·m-3），ra為空氣動力學阻力（s·m-1）。ra可由式（4）求得：

式中：κ為von karman 常數，取值0.41；z1為風速觀測高度（m）；z2為氣溫和濕度的觀測高度（m）；本文z1和z2均取值2 m；d0為零平面位移可近似為2h3，h為植被高度。zom為動力學粗糙度（m），可近似等于h8；zov為水汽傳輸粗糙度（m），zov可近似等于0.1zom。

1.3.4 CRice模型

Brutsaert在嚴格的物理邊界條件下重新定義了互補理論中的蒸散發相關變量的概念，提出了廣義非線性互補理論模型（Generalized Complementary Relationship，CR［14］）。其公式如下：

式中：ETpo為潛在蒸散發（mm），ETpo可由Priestley and Taylor 等式求得；ETpa為表觀潛在蒸散發，可用如下Penman-Monteith［12］公式：

式中：Qne=（Rn-G）/λ，f(ur)為風函數，風函數的計算選擇roma公式［12］，其公式如下：

原始的CR 模型并未考慮多年凍土區凍結期間冰相變（Qice）所消耗的能量。Wang等［20］認為在凍結期土壤熱通量G不能作為單一的能量消耗項，Qice所消耗的能量不能忽略，優化了非線性互補關系在凍土區的應用，提出了含冰能耗項的CRice模型，具體改進方法見Wang［20］等。

1.3 模型評價方法

模型評價采用確定系數（coefficient of determination，R2）、均方根誤差（root mean squared error，RMSE）、平均誤差（mean error， ME）以及納什效率系數（Nash-Sutcliffe efficiency coefficient， NSE）。R2反映模型模擬值與實測值的擬合程度，ME 反映估算值和觀測值之間的偏差，RMSE 可反映模型的絕對無偏性和極值效應，NSE 一般用來平均水文模型的模擬結果的好壞。通常情況下模型模擬結果的R2和NSE 越接近于1，ME 和RMSE 越接近0，模型的適用性越好。評價指標的計算方法見表2。

表2 模型評價指標Table 2 Model-evaluation index

2 結果與分析

2.1 風火山地表能量和氣候特征

本文選用2019年4月—2020年4月監測數據對風火山地表能量和氣候特征進行分析。圖3為風火山監測站的氣溫、降水量、風速、相對濕度及地表能量通量（凈輻射、感熱、潛熱和地表土壤熱通量）的逐日變化特征。研究區日均最高氣溫為9.8 ℃，日均最低氣溫為-22.8 ℃。年最大氣溫差達到32.6 ℃；年降雨量為385 mm，降雨與氣溫具有明顯的季節變化，呈現雨熱同期，5—10 月氣溫高降雨量大，其余時間段氣溫低，降雨十分稀少。近地層風速冬季，最大日均風速達10 m·s-1，夏季日均風速一般低于4.5 m·s-1。受冬季低溫影響，飽和水汽壓較小，相對濕度較大，多在80%以上，夏季則相對較低。地表能量分配具有明顯的季節性，冬季感熱通量大于潛熱通量，潛熱通量變化幅度趨于平緩；夏季時，隨著降雨的增加以及凍土的融化，地表含水量增加，潛熱通量上升并超過感熱通量，成為地氣間能量交換的主導能量，同時波動程度大于冬季，地表熱通量年內變化幅度較小。

圖3 2019年4月20日—2020年4月20日風火山主要氣象要素和地表能量日變化Fig. 3 Daily variations of meteorological and surface energy flux at Fenghuoshan during 2019-04-20—2020-04-20

年內蒸散發變化的趨勢總體呈現先增后減，融化期實際蒸散發波動性要大，凍結期的實際蒸散發較為平緩（圖4）。隨著氣溫升高，土壤表層以及地表積雪融化，融化期蒸散發量迅速增加，日平均蒸散發量達2.5 mm·d-1；在完全融化期，表層土壤逐漸向下融化，土壤水含量增加；同時伴隨著雨季的到來，凈輻射的增加，土壤蒸發以及植被蒸騰達到全年最大值，這一階段日均蒸散發量約為2.5 mm·d-1，蒸散發總量約為386.4 mm，約占全年的69.8%；在凍結期，土壤開始雙向凍結，土壤中液態水含量迅速減少，同時雨季結束，降雨量減少，蒸散發量迅速減少，日平均蒸散發量約為0.7 mm·d-1。完全凍結期土壤完全凍結，降雨幾乎為0，地表蒸散發量處于全年最低值，日平均蒸散發量約為0.6 mm·d-1。

圖4 2019年4月20日—2020年4月20日實際觀測蒸散發Fig. 4 Observated evaporation during 2019-04-20—2020-04-20

2.2 年尺度的模型評價

由于2017 年5 月—2018 年8 月缺測，本文選用2016年4月—2017年4月以及2019年4月—2020年4 月數據進行模型評價。由圖5 可知，AA 模擬的結果高估了實際蒸散發，FAO-PM 模擬結果低估了實際蒸散量。由4 個性能模型評價指標可知（表3），FAO-PM 線性回歸率最小為0.73，AA、PT、CRice接近于1（0.85～1.04）。AA、PT、CRice的相關性指標R2表現得較好，其值接近于0.90。納什效率系數NSE中PT、CRice比AA、FAO-PM 表現得好，其值分別為0.86、0.90。AA 高估了實際蒸散發量（ME=-0.27 mm·d-1，RMSE=0.53 mm·d-1），FAO-PM 低估了實際蒸散發量（ME=0.44，RMSE=0.74 mm·d-1），CRice、PT 模擬結果與實測值較接近。由結果分析得出，年尺度蒸散發模型對高寒草甸適應性大小排序為：CRice＞PT＞AA＞FAO-PM。

表3 年尺度蒸發模型的性能Table 3 Performance of evaporation model at the annual scale

2.3 不同凍融階段模型評價

年尺度的模擬效果要優于基于不同凍融階段（表4）。4 個模型以年尺度為時間尺度的NSE 值均大于0.5，而不同的凍融階段的NSE 值存在小于0.5。4 個模型與實測值的10 日滑動平均線比較如圖6所示，4個模型的模擬結果與實測值變化趨勢以及峰值出現時間點基本一致。

表4 不同凍融階段蒸發模型性能Table 4 Performance of evaporation model in different freeze-thawing stages

圖6 2016年、2019年模型模擬值與觀測值10日滑動平均線比較Fig. 6 Comparison of simulated evaporation and observated evaporation （10-day average） in 2016 and 2019

融化期CRice的平均值以及中位數與實測值差距較小，總體偏差優于其余3 個模型，模擬效果最優；PT、FAO-PM 不同程度低估了實際蒸散發，AA高估了實際蒸散發，其中FAO-PM 與實際蒸散發偏差最大（圖7）。CRice與PT 在2016 年融化期模擬結果在前半段高估實際蒸散發在后半段低估實際蒸散發，2019 年與之相反。該時期模型模型性能排序：CRice＞PT＞AA＞FAO-PM。

圖7 模型模擬值與觀測值對比的箱型圖Fig. 7 Comparison of simulated evaporation and observated evaporation

完全融化期波動程度大，4 個模型與實測值變化趨勢基本一致，CRice、PT 模擬結果與實際蒸散發誤差較小，FAO-PM 明顯低估實際蒸散發，AA 明顯高估實際蒸散發；該階段CRice模擬效果最優，FAOPM 總體偏差最大，AA 的中位數和均值均高于實測值，FAO-PM 均低于實測值。4 個模型性能排序：CRice＞AA＞PT＞FAO-PM。

凍結期4個模型的模擬結果與實測值均呈現下降趨勢，AA 模擬結果高估了實際蒸散發，PT、AA、FAO-PM 均低估了實際蒸散發，該階段CRice模擬效果最優。4 個模型整體偏差在不同的凍融階段較小，4個模型性能排序：CRice＞PT＞FAO-PM＞AA。

完全凍結期波動程度小，4 個模型與實測值中位數均大于均值，AA、CRice、PT 均高估實際蒸散發，FAO-PM 低估實際蒸散發。2019 年12 月—2020 年4月4個模型與實測值存在峰值不同步現象，其可能原因是受完全凍結期間高寒惡劣環境影響，實際蒸散發存在觀測誤差。該階段，CRice模擬效果最優，AA 模擬結果與實測值偏差較大；PT、FAO-PM 以及CRice與實測值總體偏差較小。4 個模型性能排序：CRice＞PT＞FAO-PM＞AA。

3 討論

本研究表明基于非線性互補理論CRice最佳。原始CR 模型在嚴格的物理邊界條件下修正干燥和濕潤條件下互補理論模型計算蒸散發所產生的誤差［14］，削減土壤凍融循環含水量變化對模型模擬性能的影響。原始CR、PT、AA 以及FAO-PM 模型都未考慮青藏高原凍融循環冰水相變中能量消耗的問題，冰從固態到氣態經歷了兩種相變，融化和蒸發，這都需要消耗能量。在計算蒸散發的理想條件下，地下冰的廣泛分布及其復雜的相變過程，需要考慮冰相變所消耗的能量，Wang等［20］在推導CRice模型時考慮了凍融循環中時期冰水相變所消耗的能量，使其模擬結果能夠更加真實地反映實際蒸散發的變化趨勢。隨著溫度升高，活動層深度的增加［5］，冰水相變可能將更加劇烈，未來青藏高原典型多年凍土蒸散發模擬研究中，考慮冰水相變所帶來的影響，對進一步減少蒸散發模擬的不確定性具有重要作用。

盡管FAO-PM 模型對不同的下墊面以及氣候條件的適應性較強，但在本文研究中發現FAO-PM在高寒草甸蒸散發模擬的應用中的與實測值存在較大偏差。在融化期、完全融化期，FAO-PM 模型模擬結果普遍低于實測值（ME 均大于0.6 mm·d-1）。姚天次等［26］在使用FAO-PM 模型對青藏高原蒸散發進行研究時，發現FAO-PM 模型所計算的蒸散發與日照時長成正相關，隨著青藏高原實際蒸散發的增加，云量增多［6］，日照時長減少，而融化期和完全融化期的蒸散發量占全年的主導地位，云量處于最大階段，因此導致FAO-PM 模型模擬結果偏小。此外，Cai 等［32］研究指出FAO-PM 應用于干旱區蒸散發模擬時，模擬結果普遍低于實測值，其原因是干旱區全年大部分時間最低氣溫Tmin＞露點溫度Tdew，低估水汽壓虧損（VPD），使得蒸散發模擬結果偏低。

AA 模型應用于高寒草甸蒸散發模擬時，結果表明融化期以及完全融化期模擬值偏大。融化期以及完全融化期時，降雨量增大，土壤及地表水分增加，此時高寒草甸所處的環境十分濕潤。在濕潤條件下AA 模型模擬的結果明顯高估了實際蒸散發，這與Crago 和Brutsaert 的研究結果相似［22］。在完全凍結期，降雨量減少，環境較為干燥，此時AA模型模擬結果偏大，然而，研究表明在干燥環境下，AA 模型模擬結果偏小［13］，產生此現象的原因可能是風函數的選取以及粗糙度長度估計的不確定性［33］，此外，受完全凍結期高寒惡劣氣候影響，實際蒸散發監測可能產生誤差。

基于輻射法的PT 模型對高寒草甸蒸散發的模擬結果次佳。PT 模型不需要空氣動力和表面阻力作為輸入參數，只需要常規的氣象數據就能驅動，粗糙長度估計的不確定性對蒸散發模擬沒有影響，同時作為PT 模型的主驅動參數凈輻射和溫度具有較低的觀測不確定性［33］，這使得PT模型模擬結果相對良好。但是PT模型忽略了陸面植被，土壤類型等對蒸散發的作用［34］。植被性質，土壤性質，會影響徑流系數，徑流系數反過來又會影響蒸散發，這會對模型的模擬性能產生一定的影響［35］。

本文對于不同方法的蒸散發模型模擬青藏高原典型多年凍土區高寒草甸實際蒸散發的效果進行了評價，結果可為該地區高寒草甸蒸散發模擬應用研究提供一定的參考。然而受限于監測地區的數量，文中僅對典型多年凍土區某一處高寒草甸進行蒸散發模擬研究，缺乏對不同植被類型的研究，但青藏高原典型多年凍土區植被主要為高寒草甸以及高寒沼澤草甸，且本文研究區域風火山具有高寒草甸、多年凍土等下墊面特征的良好代表性，因此可以認為文本研究結果在多年凍土區高寒草甸植被條件下具有較好的適用性。

4 結論

本文以青藏高原典型多年凍土區風火山小流域為研究區，使用2015 年渦度觀測數據對αPT進行標定，基于研究區2016—2017 以及2019—2020 渦度觀測數據，評價了AA、PT、FA0-PM、CRice共4 種不同蒸散發模型在高寒草甸地區，全年以及不同凍融階段蒸散發模擬過程中的模擬性能，主要結論如下：

（1）CRice模型模擬精度最高，PT 模型次之，AA、FAO-PM 模型模擬精度較差，RMSE 值分別為0.45 mm·d-1、0.49 mm·d-1、0.53 mm·d-1、0.74 mm·d-1；NSE 值分別為0.90、0.86、0.84、0.69。AA 模型模擬的蒸散發高于實測值，FAO-PM、PT 以及CRice低于實測值。

（2）在融化、完全融化、凍結以及完全凍結四個凍融階段，CRice模型模擬結果在4 個蒸散發模型中與實測值誤差較小，并且能夠反映實際蒸散發的變化趨勢，表明模型能夠較好地體現凍融過程土壤冰水相變對高寒草甸蒸散發的影響

（3）盡管本文僅采用風火山一個站點的數據對比分析了不同模型，但由于該臺站具有多年凍土、高寒草甸等下墊面特征的良好代表性，可以認為青藏高原以高寒草地植被覆蓋的多年凍土區均具有適用性。