冰期河湖與大氣的熱交換模型

楊開林

（流域水循環模擬與調控國家重點實驗室，中國水利水電科學研究院，北京 100038）

1 研究背景

河湖與大氣的熱交換模型是計算分析冰凌形成、發展、消融時空變化規律的基礎，也是太陽能利用、氣候變化、水文、水資源研究的重要課題。河湖熱交換是指敞露水面、冰蓋和雪蓋與太陽輻射、長波輻射、表面蒸發和對流、降雪或降雨之間的熱交換。

Ashton［1］歸納總結了當時河湖與大氣的熱交換模型，建議采用Angstr?m日照百分率模型和Wun?derlich 云量百分率模型計算太陽輻射，并在分析冰情的過程中采用線性化模型描述冰期河湖與大氣的熱交換，即取河湖與大氣的凈熱通量?n=-hsa(Ts-Ta)，式中：hsa為熱交換系數，Ts為河湖表面溫度，Ts為氣溫。不過，他沒有給出如何選取hsa的方法。

Angstr?m 和Wunderlich 模型的特點是采用晴天直射為云天太陽輻射計算的起始值，沒有考慮晴天散射輻射的影響。王炳忠等［2-3］以理想大氣晴天日射為云天太陽輻射計算的起始值，考慮了晴天散射輻射的影響，成果被《凌汛計算規范（SL428—2008）》采用，不過計算的過程比較復雜，需要查專用圖表。和清華等［4］利用中國54 個站1961—2000年的逐日太陽總輻射和日照百分率資料，以線性日照百分率模型為基礎，分別以天文輻射、晴天日射和理想大氣日射作為起始值建立了各站的太陽輻射回歸方程，推薦采用以天文輻射為起始值計算云天輻射。這就產生了一個問題，在河湖冰工程中采用那一個太陽輻射模型較好。

Ashton 提出的線性化模型形式已經在冰工程中得到廣泛應用［5-8］。沈洪道［6］針對北美地區，采用hsa=20.0W/（m2·℃）。蔡琳等［7］和冀鴻蘭等［8］針對黃河，當水面敞露時，取hsa=10.0W/（m2·℃）；當冰封時，取hsa=1.5 W/（m2·℃）。顯然，從上述中會發現一個問題，hsa的取值相差相當懸殊，解決這個問題的關鍵有兩個，一是河湖與大氣熱交換模型的選取；二是熱交換模型的更加合理的線性化方法。

本文的主要目的是，基于太陽能利用、氣候學、水文和水資源等學科的研究成果，建立適用于冰封河湖與大氣的非線性熱交換模型，然后利用歷史天氣資料，以及筆者團隊進行的太陽輻射、水溫、冰溫、雪溫、氣溫和風速的原型觀測，探索對長波輻射、蒸發和對流熱交換模型的更加合理的線性化方法，以揭示和發現一些重要的普遍性規律。

2 太陽輻射模型

氣候學計算云天太陽輻射的日照百分率模型是由Angstr?m（1924）提出的［1］，數學描述是：

式中：?sc為云天的日射熱通量，W/m2；?si=Pm?s0為晴天直射熱通量，W/m2；?s0為天文輻射熱通量，W/m2；P 為所有波長范圍內的大氣平均透明系數；m 為光學大氣質量；S 為日照百分率，表示實際日照時間Treal與理論日照時間Tsun的比值，即，有0 ≤S ≤1.0。

天文輻射?s0是指完全由地球天文位置決定的到達地球大氣頂界的太陽輻射，與日地距離的平方成反比，隨太陽高度角的增加而增加［9］。

和清華等［4］推薦采用以天文輻射為起始值計算云天輻射：

采用上式全國54 個站的相對誤差變化在3.33%～18.75%之間，平均為8.39%。

2.1 云量百分率模型在我國設有日照觀測站的地區很少，用云量百分率代替日照百分率已成為計算云天太陽輻射的首選因子，因為云量是氣象站的常規預報項目。

基于Angstr?m日照百分率模型，Wunderlich 提出的云量百分率模型［1］是：

式中：C 為云量，%，0 ≤C ≤100%。

基于和清華等的日照百分率模型，當取S=1-C2，由式（2）得下述云量百分率模型：

由于Wunderlich 云量百分率和Angstr?m日照百分率模型一樣，都是以晴天直射作為云天輻射的起始值，沒有考慮晴天散射輻射的影響。

當考慮晴天散射的影響時，晴天散射可按Berlage 公式（1928）計算［2］：

2.2 大氣平均透明系數P 和光學大氣質量m 的計算王炳忠等［3］研究了我國大氣平均透明系數P 和光學大氣質量m 的計算，不過過程比較復雜，Pm的計算需要查專用圖表。鑒于此，本文將采用Glover和McCulloch 的方法計算Pm［1］：

需要指出的是式（7）只有在Pm>0 才有意義。根據Klein 的研究，m 可用下式計算：

式中：Pa為觀測地的大氣壓，hPa；P0為海平面大氣壓，hPa，一般取1013.25 hPa。

2.3 模型的比較在2016年1—4月期間，筆者團隊在黑龍江省漠河市北極村黑龍江邊開展了太陽輻射的實測，下面將以此為依據比較上面的3 個云量百分率模型：以晴天直射為起始值的Wunderlich 式（3）；以天文輻射為起始值的和清華和謝云式（4）；本文晴天日射為起始值的式（6）。

圖1（a）示出了2016年4月白天平均云量隨日期的變化（資料來源https：//tianqi.911cha.com /mohe/2016 -4.html，以下歷史氣象資料來源相同），其中：n=92 表示2016年4月1日；云量C=1.0 是天降大雪的完全陰天；云量C<0.6 可視為晴天。圖1（b）（c）（d）分別示出了式（3）（4）（6）計算的2016年4月份太陽輻射的輻照度與實測值的比較，其中天文輻射?s0的計算可參考文獻［9］。

根據圖1，可得重要結論如下：

圖1 黑龍江省漠河市北極村黑龍江邊實測太陽輻射與理論計算比較

（1）式（6）計算太陽輻射結果與實測值吻合最好，式（3）和式（4）次之，例如，當n=115 時，C=0 為完全晴天，實測最大?sc=1042.1 W/m2，而式（3）（4）（6）計算的最大?sc值分別為814.6、811.5、932.1 W/m2。這表明所提出的式（6）用于計算晴天日射的模型及大氣平均透明系數和光學大氣質量的方法是比較符合實際的。

（2）當白天平均云量C≥0.6，太陽輻射理論計算值與實測值均偏差較大，特別是連續2 天C=1.0的降雪陰天。原因可能是：①白晝日平均云量資料來源于氣象站觀測資料，與觀測點存在一定偏差；②白晝不同時間的云量相差較大。改進的途徑是用白天分時段平均云量代替白天平均云量。

3 長波輻射模型

河湖長波輻射包括水體或者冰蓋和雪蓋表面的長波輻射與大氣長波逆輻射。

3.1 河湖的長波輻射河湖作為一個近乎黑體，長波輻射的熱通量可由斯蒂弗-博爾茨曼四次方定律［1］計算：

式中：?b為河湖長波輻射的熱通量，W/m2；σ=5.67×10-8W/（m2·K4）為斯特凡·玻爾茲曼常數；εw為修正系數，對水、冰和雪可取εw=0.97；Ts為河湖表面溫度，℃，對敞露水面指水面溫度，冰封時指冰蓋表面溫度，有雪蓋時指雪面溫度。在水面溫度Ts=0℃時，?b=306 W/m2。

3.2 大氣長波逆輻射大氣長波逆輻射也可按斯蒂弗-博爾茨曼（Stefan-Boltzmann）定律［1］計算：

式中：?a為長波大氣逆輻射到河面的熱通量，W/m2；Ta為河面上1.5 m 高測得的氣溫，℃；εa為晴天的大氣發射率；γa為河湖表面對大氣長波逆輻射的反射率，γa≈0.03；K 為系數。Wunderlich 根據美國Tennessee 地區觀測資料，取K=0.17［1］。

對于晴天的大氣發射率εa，我國常常采用Idso-Jackson公式，εa只是氣溫Ta的函數［10］，而Ashton［1］推薦采用Brant 公式，其εa只是水汽壓ez的函數。黃妙芬等［11］利用實測資料比較研究了國際上常用的10 個計算εa的經驗公式，包括Angstr?m 公式、Brant 公式、Idso-Jackson 公式和Iziomon 公式等，結果表明：Iziomon 公式具有最小的平均絕對誤差和均方差，具有最大的一致性指數IA 和線性相關系數LCC；Brant 公式次之；Idso-Jackson 公式具有最大的平均絕對誤差和均方差，以及最小的一致性指數IA 和線性相關系數LCC，因此，推薦采用Iziomon 公式：

式中：ez=Rhes為離河面高度z=1.5 m 處的空氣中水汽壓，hPa；es為空氣的飽和水汽壓，hPa；Rh為空氣的相對濕度。當空氣濕度飽和時，Rh=1.0。

由式（11）可知，εa不僅與氣溫Ta有關，而且與水汽壓εz或相對濕度Rh和飽和水汽壓es有關。本文在下面的分析中，也采用Iziomon 公式計算晴天的大氣發射率εa。

飽和水汽壓es是隨表面溫度Ts的升高而增加，隨Ts的降低而減小，世界氣象組織（WMO）1996年推薦采用馬格納斯-蒂托斯（Magnus-Tetons）公式計算［12］：

式中：Ts為河湖表面溫度，℃；es0≈6.11 hPa 為Ts=0 ℃時的飽和水汽壓；a 和b 為常數，對水面：Ts＞0 ，a=17.62，b=35.86；對冰面：Ts≤0，a=21.88，b=7.66。

4 蒸發和對流模型

水面蒸發量模型的研究已有200 多年的歷史。Dalton 根據水面蒸發形成原理和維持機理，綜合考慮風速、氣溫、濕度對蒸發量的影響，提出了Dalton 模型，該模型對近代蒸發理論的創立起到了決定性的作用。我國自1950年代后期也開展了水面蒸發的研究，目前已經提出了很多水面蒸發的經驗公式［13］。

Ashton［1］總結前人研究，認為只有兩個蒸發公式可以應用于負氣溫條件下，一是Rimsha 和Don?chenko 提出的俄羅斯冬季公式，另一個是由Ryan-Harleman 提出的公式，并指出只有俄羅斯冬季公式應用于實際，在北美和冰島的長期應用中，效果很好。我國東北地區也有學者采用俄羅斯冬季公式［14］。陳惠泉等［15］理論分析和實驗驗證表明，Ryan-Harleman 公式比較適用的條件是沒有風速和虛溫度>5℃（近似為），但是，對于冰封河湖，常常小于5℃。后來，也有一些專家學者提出了新的或改進的水面蒸發公式，例如陳惠泉［15］、李萬義［16］等。不過，陳惠泉和毛世民公式主要針對發電廠冷卻水系統，李萬義公式沒有考慮水面和大氣溫差的影響。 現有的研究表明溫差對蒸發的作用很大，不容忽略［15］。

基于俄羅斯冬季公式，蒸發的熱通量可描述為［1］：

式中：?e為蒸發的熱通量，W/m2；Vz為距離河湖表面1.5 m 處的風速，m/s。

根據Bowen 對流與蒸發的關系，河湖表面對流的熱通量是：

式中：?h為對流熱通量，W/m2；Pa為當地大氣壓，hPa，隨海拔高程增加而減小。

由式（13）和式（14）可見，Vz在蒸發和對流熱交換過程中的影響很大，因此，氣象站風級的觀測資料能不能應用于河湖冰情計算分析是一個引人關注的問題。為此，在2016年1—4月在漠河北極村黑龍江岸邊進行了風速的實測。如圖10所示，漠河氣象站2016年1—4月風級在2 級和10 級之間，理論平均風速與風級的計算公式是：

式中N 為風級，所以漠河理論風速Vm在2.4～26.5 m/s 之間，但是，由于漠河屬于山區，實測黑龍江邊瞬時最大風速Vz≈5.3 m/s，日平均風速Vz在0.13～1.9 m/s 之間，1月至4月的平均風速Vzm=0.77 m/s，遠遠小于理論風速。由此，可得重要結論：氣象站風級或者風速資料一般不能正確估計地勢低洼河面1.5 m 的風速，主要原因是氣象站一般建立在地勢較高、通風好的地方，而河道地勢低，通風較差。

河道冰隨時間的發展過程，除冰蓋形成初期和開河外，是一個比較平穩的過程，短暫大風速對長期冰情發展影響有限。當以冰期平均風速評估風對蒸發和對流的影響，由于漠河冰期平均風速Vz≈0.77 m/s，所以風速影響占比為：。

圖2 2016年1月—4月漠河氣象站風級與北極村黑龍江邊實測瞬時風速和日平均風速

5 河湖與大氣的熱交換模型

河湖與大氣的凈熱交換，包括太陽輻射、長波輻射及蒸發和對流是：

式中：?n為河湖表面日平均凈熱通量，W/m2；為太陽輻射凈熱通量，W/m2，asm為太陽輻射的日平均反照率。當遇到降雪天，上式中還需要增加降雪熱通量?P［1］。我國東北、西北、內蒙古和青藏高原等地區，冬季，甚至春季，降雪頻繁且降雪量較大，冰封河湖雪蓋反照率asm對氣候變化和冰厚變化具有很大影響［19-20］。新雪覆蓋的冰的表面的反射率asm約為0.9，由雪冰組成的冰蓋的反射率為0.6～0.8，與此形成對比的是，由清澈的柱狀冰（黑冰）組成的冰蓋的反射率可能較低，只有0.2，可以采用楊開林［9］提出的水冰雪太陽輻射反照率的通用模型估計。

當天氣資料齊全，包括表面溫度Ts、氣溫Ta、云量C、相對濕度Rh、風速Vz、當地大氣壓Pa等，云天太陽輻射?sc、大氣長波逆輻射?a、河湖長波輻射?b、蒸發?e和對流?h分別采用式（6）（9）（10）（13）（14）計算。

6 河湖與大氣熱交換模型的線性化

雖然我國一些地方正在建立冰情觀測站，例如南水北調工程和黃河冰情觀測站，但尚無測站能夠提供全面定時觀測Ts、Ta、C、Rh、Vz、Pa隨時間變化的資料，因此，如何把河湖與大氣的凈熱交換多參數非線性模型式（16）轉化為只是兩個自變量Ts和Ta的線性模型，且具有較高的精度，不僅具有重要的理論意義，而且具有廣泛的實用價值。

在冰期，河湖表面存在下述情況：

（1）存在敞露水面，例如南水北調中線工程漕河渡槽，全年不封凍，在氣溫Ta<0℃的情況下，水溫Ts→0℃，這時，敞露水面與大氣的溫差(Ts-Ta)較大，特別是寒潮天氣，這種情況水面蒸發和對流模型不適合采用線性化模型，否則，可能產生較大計算誤差；

（2）河湖表面冰封，氣溫Ta>0℃，這時雪面和冰面將融化，冰蓋上會形成薄薄的水膜，甚至星羅棋布小水洼，稱為水覆冰現象，這時，冰面上溫度接近融點溫度，但冰上水體溫度Ts將趨近Ta；

（3）河湖表面冰封，氣溫Ta<0℃，甚至存在雪蓋，雪蓋或裸冰表面溫度Ts與氣溫Ta的關系。

以下僅研究情況（2）和（3）。

肖建民等［14］在黑龍江省勝利水庫測得冰面溫度Tis=0.55Ta，白乙拉等［17］在黑龍江省紅旗泡水庫測得冰面溫度Tis=0.26Ta-4.6。

在2016年1—4月上旬，筆者團隊在漠河北極村黑龍江段實測了冰情，包括冰蓋下水溫、冰蓋和雪蓋內部和表面溫度的垂向分布以及氣溫，如圖3所示，其中：Tw為水溫或冰蓋、雪蓋、氣溫；hi和hs分別為冰厚和雪厚；L 為太原理工大學研制的R-T-O 冰雪情傳感器刻度。R-T-O 冰雪情傳感器長200 cm，垂直江面穿透雪蓋和冰蓋放置，下端（圖3（a））測量水溫，中段測量冰溫和雪溫，上端（圖3（b））測量氣溫Ta和雪面或冰面溫度Ts。

根據圖3，可得下述結論：

（1）當Ta<0℃時，冰蓋和雪蓋內部溫度近似線性分布；雪蓋表面溫度Ts隨著氣溫Ta的變化而變化，前者與后者的變化存在一定的時間滯后，Ts接近于Ta，相差約為1.0～3.0℃，如圖3 中A 區和B區所示；

（2）當Ta>0℃時，雪蓋逐漸消失，冰蓋內部溫度成非線性分布，如圖3 中2016/4/11 曲線所示。

圖3 實測漠河北極村黑龍江段冰情

（3）當Ta<0℃時，雪蓋下冰面溫度Tis（如點C所示處）遠遠大于氣溫Ta，根據3 個觀測點的數據，線性回歸得Tis≈0.3Ta-2.30 ，這表明雪蓋具有很好的保溫作用，雪蓋下冰蓋表面溫度遠遠大于氣溫。

（4）11月至4月初，冰蓋上都有較厚積雪，且冰蓋下水溫約等于0℃。

由于勝利水庫、紅旗泡水庫、漠河均屬于黑龍江省，冰封期降雪頻繁，冰面上厚雪覆蓋，由上面結論（3）可以判斷出肖建民等［14］和白乙拉等［17］測得的冰面溫度與氣溫關系僅適用于雪蓋下冰面溫度。

綜上所述，由于冰封河湖表面溫度Ts接近Ta，包括裸露冰蓋表面，這就產生一個思路，在Ts=Ta點將式（16）線性化，即：

當假設Ts=Ta是河湖表面的平衡溫度，即，則式（17）可簡化為Ashton［1］形式的線性化模型：

對長波輻射式（9）和式（10）在Ts=Ta求偏導數，得：

由式（11）和式（12）分別得：

對蒸發和對流式（13）和式（14）在Ts=Ta點求偏導數：

由于Ta、C 和Rh是相互獨立且隨時空變化，鑒于此，下面將研究利用氣象站歷史日平均天氣資料獲得熱交換系數hsa的方法。需要說明的是，氣象站歷史天氣資料Ta、C 和Rh與河湖實際情況是存在一定差異［18］，但其長期統計特性是比較接近的。

令式（21）中風速Vz=0，可得圖4 和圖5所示北京地區2015.12.1—2016.2.29 與2016.12.1—2017.2.29兩個冰期hsa的線性回歸經驗式：

相關系數R2=0.21～0.33，線性相關性一般。

式（22）—式（23）得：

圖4 北京2015.12.1—2016.2.29日平均天氣資料及?sa0 和hsa 的線性回歸

圖5 北京2016.12.1—2017.2.28日平均天氣資料及?sa 和Ta的線性回歸

采用與上類似的方法，可得一些其他地區hsa與Ta的線性關系。表1 中列出了一些典型地區，包括漠河、沈陽、包頭、北京、保定、拉薩的天氣特征及不同風速Vz時hsa的線性公式，其中風速假設為冰期日平均風速。黑龍江漠河是我國緯度最大、氣候最寒冷、濕度大、云量大、冰期時間長的地區，hsa的平均值最小；包頭屬于內蒙古高原，氣候干燥，少云，冰期較長；西藏拉薩緯度小，空氣干燥、少云、冰期時間長，但大氣壓較低，hsa的平均值較小；沈陽、北京、保定屬于平原地區，常年大氣壓較高，位于我國的中部和北部，空氣干燥、少云，冰期時間較短，但hsa的平均值較大。

根據表1 可得下述結論：（1）河湖與大氣的熱交換系數hsa與日平均氣溫Ta成正比；（2）當已知Vz=0.0 m/s 時的熱交換系數hsa0，則不同風速的熱交換系數可描述為：

表1 典型地區氣象特征參數及河湖與大氣熱交換模型的線性化結果

由于北京、保定、沈陽、包頭冰期日平均氣溫Ta=-17～10℃之間、云量C=0.20～0.22、Rh=0.37～0.40，hsa的值相差不大，可以采用統一的參數，當Vz=0.0 時，hsa0約為10.0 W/m2·℃。漠河和拉薩的hsa略小于北京等地區的，前者是因為冰期相對濕度較大，后者是因為海拔高程較高。

7 結束語

基于原型冰情觀測和現有研究成果，建立了適用于冰期河湖與大氣的熱交換數學模型。太陽輻射現場觀測資料證明本文提出的計算晴天日射的模型及大氣平均透明系數和光學大氣質量的方法是實用的。實測風速與氣象站風級或者風速資料的對比表明兩者差別很大。當氣溫Ta<0.0℃時，冰蓋和雪蓋內部溫度近似線性分布，且雪蓋下冰面溫度Tis≈0.34Ta-2.30 ，但是雪蓋或裸露冰蓋表面溫度Ts接近于Ta，因此在Ts=Ta點將河湖與大氣的熱交換模型線性化是最好的近似。對于一個地區，采用典型年歷史日平均天氣資料得到的熱交換系數hsa0可用于預測其他年的熱交換。 hsa與Ta成正比，當已知風速Vz=0.0 時的熱交換系數 hsa0，則不同風速的熱交換系數。北京、保定、沈陽、包頭的平均hsa0≈10.0 W/（m2·℃），漠河和拉薩的hsa略小于北京的。