少資料條件河渠冬季冰情過程模擬研究

王 濤，郭新蕾，陳玉壯，劉吉峰，佘云童，潘佳佳

（1.中國水利水電科學(xué)研究院 流域水循環(huán)與調(diào)控國家重點實驗室 北京 100038；2.阿爾伯塔大學(xué)，阿爾伯塔省 加拿大； 3.黃河水利委員會水文局，河南 鄭州 450004）

1 研究背景

北半球較高緯度的地區(qū)約有60%以上河流在冬季會經(jīng)歷冰凌過程，我國北緯30°以北占四分之三的國土面積上都有冰凌現(xiàn)象發(fā)生。河冰發(fā)展演變包括流凌、冰蓋形成、冰蓋熱增長和封河、冰蓋的熱衰減和消融開河等過程，是復(fù)雜的水文、氣象、水力學(xué)、熱力學(xué)和動力學(xué)交互作用的結(jié)果。為了準(zhǔn)確模擬和預(yù)報河渠冰凌發(fā)展全過程，國內(nèi)外學(xué)者在該領(lǐng)域做了大量工作。較為系統(tǒng)的模型可從Shen的RICE模型算起，提出的模擬河冰過程的雙層解析框架，考慮了水溫分布和冰的濃度分布以及冰蓋熱力增長和消退。Beltaos［1］研制了RIVJAM模型計算寬河型冰塞所引起的水位升高，該模型可以較好的模擬非平衡的冰塞或接地冰塞。類似的模擬模型還包括ICEJAM［2］，HEC-RAS，ICEPRO，ICESIM［3］，River 1D［4］和RIVICE［5］等。Shen等［6-8］建立了具有世界水準(zhǔn)的一維、二維DynRICE模型（CRISSP）來模擬冰凌發(fā)展過程，并考慮了河床變化和泥沙運(yùn)動，該系列模型已成功應(yīng)用于世界多條河流的冰情研究中。最近Wazney等［9］在上述模擬冰蓋形成演化時提出了新的公式以嘗試建立熱力學(xué)和動力學(xué)之間的聯(lián)系并得到驗證。國內(nèi)楊開林［10］、王軍［11］、茅澤育［12］、穆祥鵬［13］、李潤玲［14］等也建立了冰塞形成及演變發(fā)展動態(tài)數(shù)學(xué)模型、河道冰塞堆積厚度的數(shù)值模型及其冰厚增長度日法等，計算了典型河道中的冰塞堆積厚度，郭新蕾等［15］針對長距離明渠-閘門-泵站系統(tǒng)冬季正、反向輸水可能出現(xiàn)的冰問題，開發(fā)了調(diào)水工程冬季輸水冰情過程模擬平臺并在南水北調(diào)、密云水庫調(diào)蓄工程中取得較好模擬效果。此外，Morse［16］、郭新蕾［17］、王濤［18］、Sun等［19-20］先后將統(tǒng)計模型、人工神經(jīng)網(wǎng)絡(luò)、模糊數(shù)學(xué)方法等與冰情預(yù)測預(yù)報結(jié)合，較為成功的預(yù)測了封開河時間、冰厚、冰壩發(fā)生日期等冰情參數(shù)，也應(yīng)用到黃河、黑龍江等河流，取得了較好效果。

上述進(jìn)展表明冬季河冰過程的模擬研究已取得較大進(jìn)展，但上述數(shù)學(xué)模型在一些天然河道冰情模擬應(yīng)用時常常會遇到一個難題，即缺乏細(xì)致系統(tǒng)的河道邊界條件、水文、氣象和水力學(xué)基礎(chǔ)資料。具體表現(xiàn)一是國內(nèi)水文站主要設(shè)立在較大流域上，而且站與站之間斷面間距長，用于計算的水文站數(shù)據(jù)空間步長大；二是測量的水文和氣象資料通常為日均值或者每日最大、最小值，在時間分布上現(xiàn)有數(shù)據(jù)的時間步長太大；三是天然江河支流縱橫，但支流通常未設(shè)水文觀測站，支流流入和流出的水量未知；四是水文站觀測的水溫數(shù)據(jù)通常精度不夠或數(shù)據(jù)不完整；五是太陽輻射和氣溫等相關(guān)氣象資料采用氣象局發(fā)布的信息，氣象站同水文觀測站的距離較遠(yuǎn)，數(shù)據(jù)同步性差。比如我國黑龍江上游900 km也只有8個基本水位站，主要觀測水位信息［21］。黃河是國內(nèi)水文站系統(tǒng)較為健全、水文觀測數(shù)據(jù)相對系統(tǒng)的河流，但黃河內(nèi)蒙古河段823 km也只有4個水文站［22］，2014年以后也只增加了包頭水文站。1990年代，黃河水利委員會先后同芬蘭和美國合作開發(fā)了黃河下游冰情預(yù)報數(shù)學(xué)模型，但黃河上的觀測資料不能滿足模型計算中對實測資料精度和種類的需求，再加之黃河局部河段為游蕩性河道，河床地形變化頻繁，導(dǎo)致模型建成后很快也不能適應(yīng)實際運(yùn)行應(yīng)用的需要［23-25］。目前，北方天然河道冰情模擬中普遍存在著實測資料不全或缺失、河道斷面資料難以測量等條件制約。

鑒于此，本研究提出少資料條件河渠冬季冰情發(fā)展過程模擬的實用方法和處理技術(shù)，重點研究了支流流量未知時的流量分配方案、每時氣溫和太陽輻射過程的計算方法等，并將上述模型和方法應(yīng)用到黃河內(nèi)蒙河段冰情模擬計算中，目的是為這種少資料條件河渠的冬季冰情過程模擬提供方法支撐。

2 少資料條件下支流流量分配

天然江河進(jìn)行水動力學(xué)計算時如果沿途支流流量信息缺乏，有必要根據(jù)實測流量和水位進(jìn)行模擬反演，通過辨識沿途支流流量分配比例以達(dá)到計算中區(qū)間流量的平衡。

明渠非恒定流的連續(xù)性方程和運(yùn)動方程［26］為

式中：x為距離，m；t為時間，s；z為水位，m；Q為流量，m3/s；A為過流面積；B為水面寬；ql為單位流程上的側(cè)向出流量，負(fù)值表示流入。

在進(jìn)行數(shù)值計算時，需要給出研究河段的初始條件、上下游邊界、支流的流量邊界條件。邊界條件通常有3種類型，分別為水位邊界條件Z=Z（t），流量邊界條件Q=Q（t），水位流量邊界條件Q=Q（Z）。計算中通常采用上游邊界和支流邊界給定初始流量、下游邊界給定初始水位的方法，前后計算斷面應(yīng)滿足流量平衡和能量守恒。

計算河段支流邊界示意圖如圖1所示，其中邊界條件包括6個：上游邊界①，下游邊界②，支流邊界③-⑥，在冰情模擬時常會遇到只有上游邊界和下游邊界的流量信息。因支流③-⑥流量未監(jiān)測，計算中為了滿足流量平衡，需要對區(qū)間匯入和流出的流量進(jìn)行合理分配。如果僅僅簡單采用上、下游邊界流量差進(jìn)行比例分配的話，忽略了實際流量自上游向下游傳遞的時間差很可能導(dǎo)致分配方案的不合理。本研究提出下述流量分配方法進(jìn)行支流流量分配，可確保計算河段流量的平衡。

圖1 計算河系支流邊界示意圖

在計算中只輸入上游邊界條件①的流量過程和下游邊界②的水位條件，上游邊界①流量傳播到各支流邊界③-④-⑤-⑥-下游邊界②的模擬值為Qi（i=1，2，3，4），以黃河內(nèi)蒙古河段巴彥高勒至三湖河口4個分退水口作為支流邊界計算的流量變化曲線如圖2所示，其中上游邊界①為巴彥高勒站，下游邊界為②三湖河口水文站，圖2中支流邊界③-⑥模擬值和下游邊界2的模擬值分別為支流未分配流量條件下計算結(jié)果。圖2中可以看出，下游邊界②模擬流量和實測流量存在明顯差異，因為模擬中區(qū)間支流流量未計入，導(dǎo)致流量存在不平衡。為此，需給支流邊界③-⑥分配適當(dāng)流量。由于支流邊界③-⑥并無流量觀測數(shù)據(jù)，采用實測下游邊界②的流量Qdown分別減去流量傳遞到③-⑥支流斷面處的模擬流量，差值按照比例ai（i=1，2，3，4）分配到各支流邊界③-⑥作為計算中支流分配的流量qi（i=1，2，3，4），確保計算河段流量平衡，各支流分配的流量qi為：

流量分配系數(shù)αi的確定可通過調(diào)整各支流邊界αi值以計算下游邊界②流量模擬值，再與實測值對比，以使關(guān)鍵時間點上的流量值最接近率定得到，或可以嘗試?yán)米钚《朔▽?yōu)達(dá)到相對準(zhǔn)確辨識區(qū)間支流流量分配比例的目的，確保計算河段流量的平衡。通過支流邊界③-⑥流量的合理分配，計算得到支流分配流量后下游邊界②的模擬值如圖2所示，由圖2可知，下游邊界②三湖河口水文站模擬值和實測值吻合良好。

圖2 未分配支流流量傳遞過程模擬曲線（實測值暫按±5%的誤差考慮）

3 少資料條件下每時太陽輻射的計算

太陽輻射資料來自氣象局，是以日為單位的凈太陽輻射，而模型計算需要每小時的太陽輻射值，本研究采用計算和參數(shù)率定相結(jié)合解決缺少每時太陽輻射資料的問題。

太陽輻射計算，參考文獻(xiàn)［27］無云狀態(tài)下的短波太陽輻射φcl，可以由下式計算：

式中：φcl為短波凈太陽輻射，W/m2；φso為每單位面積總的外來太陽輻射，J/m2；Iso為太陽輻射常數(shù)，冬季約為1380 W/m2；ω是時角，午時為0°，每小時變化15°，早上為正值，下午為負(fù)值；δ為太陽傾角，弧度；φ為緯度，度，北緯為正，南緯為負(fù)；dn為一年中的天數(shù)，從1月1日算起，如1月20日=20 d；m為當(dāng)?shù)貧鈮?pa時海拔z米的光學(xué)氣團(tuán)，m； p0為水平面氣壓，Pa；m0為水平面處的光學(xué)氣團(tuán)；太陽緯度α=90-θz；θz=arccos( )sinδsinφ+cosδcosφcosωi，弧度； ωi為前一個小時與當(dāng)前小時時角的均值；E0為地球軌道的偏心校正系數(shù)。

φcl為無云狀態(tài)下的太陽輻射，在有云條件下，太陽輻射將減小，用下式估算［28］：

式中：φri為有云狀態(tài)下的太陽輻射，W/m2；C為云的覆蓋度，分為10個等級，C=0表示晴空無云，C=10表示烏云全部覆蓋。如果云的狀態(tài)沒有數(shù)據(jù)，通常采用估計數(shù)據(jù)。

當(dāng)太陽輻射到達(dá)水面或冰面，部分輻射將反射回大氣，太陽輻射公式再次修正為φR：

式中Rt為反射率，同緯度和照射物體表面性質(zhì)有關(guān)。

通常在冰情演變過程計算中，需將河段分為盡可能多的區(qū)域，確保計算斷面使用的氣象數(shù)據(jù)的準(zhǔn)確度，每個區(qū)域計算一個站點的太陽輻射和氣溫。通過式（4）—式（10）可計算出晴空無云狀態(tài)下冰情太陽輻射，式（11）計算實際有云時狀態(tài)太陽輻射，但是云的覆蓋程度沒有資料記錄，本研究采用每天日照數(shù)同統(tǒng)計時段2012年11月—2013年3月內(nèi)最大日照數(shù)比值衡量云層覆蓋狀態(tài)。于是，式（11）可修正為：

式中：SSH為計算日日照數(shù)；SSHmax為統(tǒng)計時段內(nèi)最大日照數(shù)。

冰情演變過程中太陽對冰面和水面的反射率不僅跟冰的表面狀況有關(guān)，還和河流中水的含沙量、水的渾濁度等相關(guān)。實際計算中不同的河流不同冰狀態(tài)反射率都不相同，一維模型計算中考慮太陽輻射時通常忽略太陽反射率的影響，導(dǎo)致熱擴(kuò)散計算中存在誤差，本研究針對不同時段、不同河道冰情狀況，采用不同的系數(shù)反應(yīng)太陽反射率的變化，通過模型率定反演得到表達(dá)太陽反射率的系數(shù)。在冰情模擬中該系數(shù)分為流凌前Rb1、流凌期Rb2、封河期Rb3、開河前期Rb4、開河期Rb5、冰蓋表面有積雪覆蓋Rb66種反射率，不同時段冰水情狀態(tài)如表1所示，則式（12）可寫為：

表1 不同冰情狀態(tài)表達(dá)太陽反射率的系數(shù)

式中Rbi為表達(dá)太陽輻射變化的系數(shù)，i=1，2，3，4，5，6。通過冰情實際情況率定得到。Rbi同冰面或水面太陽反射率Rti、河道情況Rii、河道水質(zhì)Wai、冰情情況Ici和積雪覆蓋情況Sni等有關(guān)。

黃河內(nèi)蒙河段有磴口、臨河、五原、烏拉特前旗、包頭、土右旗氣象局、托克托縣7個典型氣象站（如圖3），黃河內(nèi)蒙古河段冰情計算中，分別計算上述7個氣象站太陽輻射值，7個氣象站太陽輻射值分布規(guī)律和走勢一致，在這里只列舉出臨河的太陽輻射計算結(jié)果，臨河氣象站位于東經(jīng)40.73°，北緯107.37°，高程1041.1 m，圖4為2012年11月1日—2013年3月31日臨河每時無云狀態(tài)下凈太陽輻射φcl（W/m2）計算值，圖5列出了2012年11月臨河每時無云狀態(tài)下凈太陽輻射。為了更清楚了解太陽輻射局部變化規(guī)律，根據(jù)計算太陽輻射值對應(yīng)的年份和日期，可率定出不同結(jié)冰時段和冰蓋狀態(tài)下太陽反射率的影響，并根據(jù)日照數(shù)估算出云層遮擋度的影響，代入式（14）即可求出符合實際情況的太陽輻射值。

圖3 氣象站分布和氣象數(shù)據(jù)分區(qū)示意圖

圖4 2012年11月—2013年3月臨河每時太陽輻射計算值

圖5 2012年11月臨河每時太陽輻射計算值

4 少資料條件下每時氣溫的計算

水文站和氣象站記錄的氣溫值為日最高值、日最低值或者日均值，且站與站之間距離遠(yuǎn)，冰情模擬需要每小時的氣溫值，現(xiàn)有資料不能滿足模擬模型的需求。本研究中利用Parton and Logan（1981）［29-30］的氣溫計算公式進(jìn)行推演。

式中：Ti為白天或者夜晚的第i個小時的溫度；Y為白天時長，h；Z為夜晚時長，h；Tmax和Tmin為日最高和最低氣溫；TS為日落時的氣溫；m為最低溫度出現(xiàn)后到日落的小時數(shù)；n為日落后到最低溫度時間的小時數(shù)；a為最高氣溫的滯后系數(shù)；b為夜間氣溫系數(shù)；c為從日出時開始的最低溫度的滯后系數(shù)；HR為計算溫度對應(yīng)的小時數(shù)，取值為1～24 h。系數(shù)a、b、c是由土壤和氣溫決定的，根據(jù)當(dāng)?shù)貙崪y氣溫進(jìn)行率定得到，本研究率定得到的a=0.15，b=2.28，c=-0.2。

仍以圖3所示氣象站資料來說明，根據(jù)氣象站提供的日最大氣溫和最小氣溫值率定所需站點的時均氣溫，圖6為2012年11月1日—2013年3月31日臨河每時氣溫計算值，圖7將2012年11月每時氣溫單獨列出，便于更清楚觀察氣溫的局部變化規(guī)律。

圖6 2012年11月—2013年3月臨河每時氣溫計算值

圖7 2012年11月臨河每時氣溫計算值

5 黃河內(nèi)蒙古河段冰情模擬的應(yīng)用

將上述少資料條件河渠冬季冰情發(fā)展過程模擬的實用方法和處理技術(shù)應(yīng)用到黃河內(nèi)蒙古河段的冰情模擬中。

該段最新河道斷面數(shù)據(jù)為2012年10月的大斷面資料，從巴彥高勒到頭道拐共測量斷面166個，因此本研究冰情計算、驗證、校核及冰情模擬以2012—2013凌汛期作為研究對象。冰情演變過程模擬中必需的物理量氣溫、太陽輻射和支流流量分配數(shù)據(jù)均不能滿足計算需要，需要采用上述方法對所需資料進(jìn)行推演和計算。

5.1 少資料條件支流流量分配 模擬河段匯入支流眾多，沿途分布多處用于工業(yè)、農(nóng)業(yè)、畜牧業(yè)和居民飲水的取水和退水口，支流匯入沒有水文觀測資料，沿途取水資料模糊，因此需要通過實測流量和水位模擬，對沿途流量進(jìn)行分配。內(nèi)蒙古河段只有巴彥高勒、三湖河口和頭道拐3個水文站流量資料。為了確保巴彥高勒和三湖河口、三湖河口和頭道拐之間流量平衡，需按“少資料條件支流流量分配方法”將流量進(jìn)行分配。計算中采用的支流邊界如表2所示，巴彥高勒-三湖河口之間分配支流邊界4個，分別為支流邊界1、支流邊界2、支流邊界3和支流邊界4，三湖河口和頭道拐之間分配支流邊界6個，分別為支流邊界5、支流邊界6、支流邊界7、支流邊界8、支流邊界9、支流邊界10和頭道拐，根據(jù)研究中區(qū)間流量分配方法，經(jīng)多次反算率定，建議巴彥高勒至三湖河口區(qū)間支流流量分配系數(shù)比例αi分別為0.1∶0.1∶0.3∶0.5，三湖河口至頭道拐區(qū)間支流流量分配系數(shù)比例αi分別為0.1∶0.1∶0.1∶0.1∶0.3∶0.3。

5.2 少資料條件太陽輻射和氣溫的模擬 由于研究河段較長，流經(jīng)區(qū)域氣候條件有明顯差異，因此將計算流域劃分為8個氣象數(shù)據(jù)計算區(qū)（如圖3所示），區(qū)域1—區(qū)域8所采用的氣溫和太陽輻射資料分別對應(yīng)氣象站為：磴口、臨河、五原、烏拉特前旗、包頭、包頭、土右旗氣象局、托克托縣。根據(jù)有限的氣象資料，采用本研究少資料條件每時太陽輻射和氣溫的計算方法，可計算出上述8個水文站每時氣溫和每時太陽輻射作為8各區(qū)域計算輸入值。其中臨河氣溫和凈太陽輻射計算如圖4—圖7所示。

5.3 黃河內(nèi)蒙河段冰情發(fā)展的模擬 冰情過程模擬采用河冰過程冰水基本數(shù)學(xué)模型，包括：明渠非恒定流模型、水流的熱擴(kuò)散模型、冰花擴(kuò)散模型、冰蓋下的水流輸冰能力模型、水面浮冰的輸運(yùn)模型、冰蓋和冰塊厚度發(fā)展模型、冰塞下冰花含量和冰塞厚度模型等。通過相關(guān)參數(shù)率定選擇出適合黃河冰情演變規(guī)律的較為優(yōu)選的方案，計算中取水和空氣熱交換系數(shù)hwa=16（W/m2/°C），平鋪上溯模式發(fā)展的最大費(fèi)勞德數(shù)Fr_jux=0.06和水力加厚模式發(fā)展的最大費(fèi)勞德數(shù)，該值分段給出，即：1～50 km河段Fr_max=0.11，50～174 km河段Fr_max=0.13，174 km至下游段河段Fr_max=0.097。

巴彥高勒封河時間為2012年12月23日。圖8為巴彥高勒固定冰蓋厚度發(fā)展過程模擬值和實測值的比較，表2為冰蓋厚度模擬值和實測值的誤差。2012—2013年度冰情信息顯示，截止到2012年12月27日研究河段全線封凍，根據(jù)冰厚測量數(shù)據(jù)顯示，到2013年2月1日，冰厚已充分發(fā)展達(dá)到最大值。巴彥高勒整個冬季冰蓋模擬值和實測值均方根誤差為0.065 m，冰蓋增長期（封河—2013年2月1日）均方根誤差為0.056 m；巴彥高勒整個冬季實測量值和模擬值絕對誤差均值是0.059 m，冰蓋增長期絕對誤差均值分別為0.041 m。模擬結(jié)果顯示冰蓋增長期模擬誤差小于冰蓋消融期，原因是：在冰厚增長期，熱力作用為冰蓋厚度增加的主要因素，而在冰蓋消融期，冰蓋并不是按照靜止不動通過熱力作用就地融化消失，通常冰蓋消融期是由熱力和動力共同作用，但在一維冰情模擬中，冰蓋發(fā)展的數(shù)值模擬尚未能考慮到動力作用，所以冰蓋模擬中消融期計算誤差比增長期大。

圖8 巴彥高勒冰厚模擬值和實測值比較

表2 冰蓋厚度模擬值和實測值的誤差

內(nèi)蒙古河段最先封凍時間為2012年11月30日，位于昭君墳上游附近（109°54′12″E和40°31′53″N），該位置距離上游巴彥高勒水文站286.7 km，從該位置冰蓋向上游發(fā)展過程如圖9所示，2012年12月23日封凍上首發(fā)展到上游巴彥高勒。從圖9上看出冬季冰蓋前沿發(fā)展過程模擬值和實測值吻合較好，整個冬季模擬值和實測值的均方根誤差為18.76 km，絕對誤差均值6.43 km；2013年3月5日前模擬值和實測值的均方根誤差僅為6.27 km，絕對誤差均值2.50 km，模擬的冰蓋前沿發(fā)展過程與觀測值基本對應(yīng)，但仍然存在開河期模擬值和實測值誤差大于冰蓋增長期和穩(wěn)定期的情況。由上面的對比來看，提出的時均氣溫和時均太陽輻射計算方法、支流流量的動態(tài)分配方法能夠較好的解決冰情計算中相關(guān)少資料的難題。

圖9 冰蓋前沿的發(fā)展過程

6 結(jié)論

針對計算河段區(qū)間退水和取水資料缺少問題，本文提出了流量動態(tài)分配方法，并應(yīng)用到黃河巴彥高勒到頭道拐區(qū)間10個分水口支流流量分配中，通過合理的比例分配支流流量，確保了計算區(qū)間流量的平衡。為了得到計算中時均太陽輻射資料，在計算凈太陽輻射基礎(chǔ)上，提出了云層遮擋計算方法和不同結(jié)冰期太陽反射率影響的推演方法，通過計算水溫和實測水溫的比較，率定出冬季不同冰情時段時均太陽輻射值。時均氣溫采用Parton和Logan的氣溫計算公式通過實測日最低、最高氣溫率定得到。將上述方法得到的支流流量、氣象站的時均氣溫和太陽輻射應(yīng)用到黃河內(nèi)蒙古河段冰情模擬中，通過對比巴彥高勒冰厚和冰蓋前沿發(fā)展過程模擬值和實測值的均方根誤差和絕對誤差均值對比，表明總體上模擬值和實測值吻合較好，證實所提出的時均氣溫和時均太陽輻射計算方法、支流流量的動態(tài)分配方法能夠解決冰情計算中相關(guān)資料缺少的難題。

本文提出的不同河段支流流量分配比例的優(yōu)化方法尚待進(jìn)一步研究，不同冰蓋結(jié)構(gòu)和封河期太陽輻射率的變化將通過原型觀測和理論研究相結(jié)合進(jìn)一步驗證和率定，并探索其變化規(guī)律。隨著水文和氣象觀測儀器的改進(jìn)，觀測數(shù)據(jù)更加詳實，支流流量數(shù)據(jù)、時均氣溫和太陽輻射數(shù)據(jù)的逐步完善，冰情模擬的精度會進(jìn)一步得到提高。