楊開林

（流域水循環(huán)模擬與調(diào)控國家重點實驗室，中國水利水電科學(xué)研究院，北京 100038）

1 研究背景

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

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

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

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

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

2 太陽輻射模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 模型的比較在2016年1—4月期間，筆者團隊在黑龍江省漠河市北極村黑龍江邊開展了太陽輻射的實測，下面將以此為依據(jù)比較上面的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］。

根據(jù)圖1，可得重要結(jié)論如下：

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

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

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

3 長波輻射模型

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

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

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

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

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

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

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

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

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

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

4 蒸發(fā)和對流模型

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

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

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

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

根據(jù)Bowen 對流與蒸發(fā)的關(guān)系，河湖表面對流的熱通量是：

式中：?h為對流熱通量，W/m2；Pa為當(dāng)?shù)卮髿鈮?，hPa，隨海拔高程增加而減小。

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

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

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

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

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

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

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

當(dāng)天氣資料齊全，包括表面溫度Ts、氣溫Ta、云量C、相對濕度Rh、風(fēng)速Vz、當(dāng)?shù)卮髿鈮篜a等，云天太陽輻射?sc、大氣長波逆輻射?a、河湖長波輻射?b、蒸發(fā)?e和對流?h分別采用式（6）（9）（10）（13）（14）計算。

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

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

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

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

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

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

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

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

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

根據(jù)圖3，可得下述結(jié)論：

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

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

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

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

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

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

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

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

對長波輻射式（9）和式（10）在Ts=Ta求偏導(dǎo)數(shù)，得：

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

對蒸發(fā)和對流式（13）和式（14）在Ts=Ta點求偏導(dǎo)數(shù)：

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

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

相關(guān)系數(shù)R2=0.21～0.33，線性相關(guān)性一般。

式（22）—式（23）得：

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

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

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

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

表1 典型地區(qū)氣象特征參數(shù)及河湖與大氣熱交換模型的線性化結(jié)果

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

7 結(jié)束語

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