陳 軍，劉 意，李 婷

(成都信息工程大學(xué)資源與環(huán)境學(xué)院,成都 610225)

0 引 言

近年來，因人類活動導(dǎo)致氣候異常，暴雨洪澇災(zāi)害頻發(fā)。為減少暴雨洪澇導(dǎo)致的財(cái)產(chǎn)和人員損失，亟須利用徑流匯流模型對流域進(jìn)行暴雨洪澇災(zāi)害預(yù)測。基于物理機(jī)制的徑流匯流模型一般以網(wǎng)格為單元進(jìn)行匯流模擬，可獲取全流域的水流場空間和時(shí)間變化，逐漸成為暴雨洪澇災(zāi)害預(yù)測的重要手段[1-6]。

然而，流域匯流是一個(gè)非常復(fù)雜的過程，除降水之外，還存在上游河道來水、冰川融水、地下水補(bǔ)給等。由于上述因素難以精確地在流域上按網(wǎng)格單元量化，現(xiàn)有徑流匯流模型采用近似的方法進(jìn)行模擬。以SWAT模型為例，它將多年的降雨數(shù)據(jù)代入到無水流域中反復(fù)迭代，將模型的最終輸出作為流域前期水情。這種方式不僅運(yùn)算量大，還存在模擬誤差的迭代問題，因而所得到的前期水情與實(shí)際存在較大差距[7-11]；作為當(dāng)前著名的水動力模型之一，F(xiàn)loodArea由于能提供流域任一柵格單元的水深模擬結(jié)果，被廣泛應(yīng)用于暴雨洪澇模擬、山洪風(fēng)險(xiǎn)分析[12-14]，但該模型在徑流模擬時(shí)，前期水情僅作為輸入?yún)?shù)由用戶提供。

本文以納維-斯托克斯方程(N-S方程)為指導(dǎo)，提出一種基于穩(wěn)定流場的流域徑流匯流模擬方法，該方法以穩(wěn)定流場作為前期水情生成和徑流匯流模擬的基本依據(jù)。最后，以岷江上域鎮(zhèn)江關(guān)-校場壩段為例，分析方法的有效性。

1 研究區(qū)概況

岷江上游位于四川省北部，流域主體主要位于松潘縣和茂縣，地形以山地地貌為主，地勢北高南低。在西風(fēng)南支急流、暖濕東南季風(fēng)以及地形條件的共同影響下，其氣候特征主要表現(xiàn)為冬干夏濕，干濕季節(jié)分明，屬于非地帶性的干旱河谷氣候。流域年降水量800～1 200 mm，夏季降水高度集中，降水曲線呈單峰變化，鎮(zhèn)江關(guān)以上春季有極少量的雪水補(bǔ)入，年內(nèi)徑流季節(jié)變化與降雨季節(jié)相應(yīng)。岷江上游多年平均年徑流量松潘縣鎮(zhèn)江水文站約為18 億m3，匯口下不遠(yuǎn)的五通橋水文站多年平均年徑流量約780 億m3。區(qū)域大面積暴雨發(fā)生或隨洪峰傳播，干流峰高量大，嚴(yán)重威脅中下游安全[15]。為降低人類活動對徑流匯流模型精度驗(yàn)證的影響，選取岷江上游的鎮(zhèn)江關(guān)和校場壩兩個(gè)水文站，將其所在河道位置分別作為流域出口，計(jì)算流域重疊范圍作為研究流域，如圖1所示。流域上游入水口為鎮(zhèn)江關(guān)水文站，下游出口為校場壩水文站，干流河道總長約35 km，匯流面積約為580 km2。

圖1 岷江上游研究區(qū)域Fig.1 The study area of the upper reaches of Minjiang River

2 基于納維-斯托克斯方程的徑流匯流模型

納維-斯托克斯方程(Navier-Stokes equations)是描述黏性不可壓縮流體動量守恒的經(jīng)典運(yùn)動方程，反映了黏性流體(又稱真實(shí)流體)流動的基本力學(xué)規(guī)律，在流體力學(xué)中有十分重要的意義[16]。但它是一個(gè)非線性偏微分方程，只有在極少數(shù)簡單流動的情況下才有解析解[17-20]。本文以N-S為指導(dǎo)，探索柵格模式下的徑流匯流模擬。

2.1 水流運(yùn)動的N-S方程

N-S方程是一個(gè)用來描述流體運(yùn)動特征的最基本力學(xué)方程組。由于水流運(yùn)動是同性不可壓的流體運(yùn)動，可以視為是無黏性的流體，N-S公式表達(dá)為方程組：

(1)

▽μ=0

(2)

式中：u為速度場；t為時(shí)間；ρ為流體密度；p為壓強(qiáng)；F是作用于流體上的質(zhì)量力(外力)。

在徑流匯流模擬過程中，流體為實(shí)際的水體，采用水深和水速描述水體的屬性。N-S方程中通過外力項(xiàng)和壓力項(xiàng)求解以模擬水體在重力、摩擦力、壓力梯度等作用下的水速大小及變化；通過平流項(xiàng)求解以模擬水流因流動在流域各位置上水深大小及變化。

2.2 柵格模式下的N-S方程分解

柵格模式的徑流匯流模擬將流域從空間上和時(shí)間上均劃分為相等的片段。在空間上，按一定分辨率將流域劃分為大小相等的矩形網(wǎng)格，將坡面匯流和河道匯流統(tǒng)一描述為網(wǎng)格的水流運(yùn)動；在時(shí)間上，徑流匯流過程分解為相等間隔的時(shí)間片。在每一個(gè)時(shí)間片上，流域徑流匯流過程分解為柵格單元為基本單位的模擬計(jì)算，上一個(gè)時(shí)間片柵格單元模擬的水深和水速，作為下一個(gè)時(shí)間片初始水深和水速。通過時(shí)間片的迭代，實(shí)現(xiàn)流域的徑流匯流過程模擬。

設(shè)流域中任一柵格單元(中心柵格)為c，時(shí)間片初始時(shí)間為tn，終止時(shí)間為tn+1，時(shí)間片的時(shí)間間隔為Δt，則：

Δt=tn+1-tn

(3)

中心柵格c的徑流匯流模擬需要依據(jù)tn時(shí)刻中心柵格及8鄰域柵格單元b的水深和水速，模擬得到tn+1時(shí)刻中心柵格的水深和水速。

(4)

(5)

式中：Δx、Δy為鄰域柵格坐標(biāo)相對中心柵格坐標(biāo)的偏移，單位為像素；vmax表示時(shí)間片上允許的最大速度標(biāo)量。

中心柵格剩余水量比例Sc為：

(6)

表1給出了柵格單元不同的水速方向上，流出水量和剩余水量與初始水量之間的關(guān)系。從表1可見，以式(5)和式(6)為基礎(chǔ)的多流向法平流項(xiàng)模擬算法滿足質(zhì)量守恒定律。

表1 不同水速方向上流出水量和剩余水量關(guān)系表Tab.1 The relation table of outflow and remaining water amount in different water velocity directions

依據(jù)式(5)和式(6)，平流計(jì)算后中心柵格水深為：

(7)

(8)

依據(jù)式(5)和式(6)，當(dāng)模擬網(wǎng)格水速大小等于vmax時(shí)，時(shí)間片上的網(wǎng)格水流全部流出。因此，vmax決定了水質(zhì)模擬的時(shí)間片間隔，如下式：

(9)

式中：C為柵格分辨率。

從式(9)可見，vmax越大，時(shí)間片時(shí)間間隔越小，時(shí)間精度越高，但1 h需要迭代的次數(shù)越多；另外，為保證水量平衡，當(dāng)柵格單元c的模擬速度大于vmax時(shí)，設(shè)置為vmax。

(3)中心柵格c的水動力模擬。水動力模擬對應(yīng)于N-S方程的壓力項(xiàng)和外力項(xiàng)求解。設(shè)中心柵格高程為Hc，鄰域柵格高程為Hb，則水體高度差ΔHb為：

(10)

(11)

式中：α為與水體密度、重力加速度、摩擦力等相關(guān)的正數(shù)常量，通過流域?qū)嶋H水文數(shù)據(jù)進(jìn)行率定。

地表水流由于受到地表摩擦力的影響，使水流速度發(fā)生垂直變化。在徑流匯流模擬過程中，當(dāng)水深低于給定閾值時(shí)，將水速乘以衰減系數(shù)以減小水流速度。衰減系數(shù)ε近似定義為：

(12)

式中：dmin和dmax分別為水流摩擦力的下界深度和上界深度，m；σ為摩擦比例系數(shù)，通過流域?qū)嶋H水文數(shù)據(jù)進(jìn)行率定。

3 基于穩(wěn)定流場的徑流匯流模擬

3.1 穩(wěn)定流場的基本概念

假定流域上游持續(xù)不間斷地提供穩(wěn)定水源，通過徑流匯流模型不斷迭代，流域出口處的水深和水速最終趨于穩(wěn)定。此時(shí)的流域水流場稱為穩(wěn)定水流場。

顯然，現(xiàn)實(shí)中不存在穩(wěn)定水源，穩(wěn)定流場也不可能存在。但在流域上游設(shè)定穩(wěn)定供應(yīng)水源，通過模型迭代，如果能在流域出口處最終模擬出穩(wěn)定的水深和水速，則徑流匯流模型具有一定的合理性。

設(shè)上游來水點(diǎn)穩(wěn)定供水量為ξ(單位為m3/h)，則考慮上游來水的柵格單元水量收支由式(4)修改為式(13)：

(13)

式中：ξc為時(shí)間片上柵格單元分?jǐn)偟纳嫌蝸硭俊?/p>

根據(jù)上游來水的定義，在上游來水點(diǎn)附近，ξc>0，其余位置上ξc=0，并且ξc滿足：

ξ=∑ξcN

(14)

式中：N為小時(shí)需迭代的時(shí)間片數(shù)，由式(9)換算得到。

合理的徑流匯流模型是穩(wěn)定水流場建立的關(guān)鍵。穩(wěn)定水流場創(chuàng)建過程中，采用式(15)定量描述穩(wěn)定流場是否創(chuàng)建完成。

(15)

式中：W為流域下游控制點(diǎn)水深；n為當(dāng)前已模擬的小時(shí)數(shù)；h為前推小時(shí)數(shù)；A為累積誤差閾值。

3.2 徑流匯流模型完善與穩(wěn)定流場生成

直接按N-S方程分解得到的徑流匯流模型存在“匯流水流方向異?！眴栴}。在流水不暢的河道區(qū)域，如果模擬水速過大，上游柵格水流快速向中心柵格聚集，中心柵格水流又不能及時(shí)向下游流動，導(dǎo)致短時(shí)間內(nèi)產(chǎn)生較大的水體高差，進(jìn)而對周圍8鄰域柵格產(chǎn)生更大速度增量。隨著迭代次數(shù)的增加，這些位置上很容易產(chǎn)生水流方向異常，甚至出現(xiàn)相鄰上下游柵格流速方向相反并頻繁交替的現(xiàn)象，導(dǎo)致水深和水速變化劇烈。表2展示了模型完善前某柵格單元典型時(shí)間片上的柵格水速和水深變化。從表2可見，現(xiàn)有模型違背了現(xiàn)實(shí)徑流匯流規(guī)律，即使在上游穩(wěn)定供水情況下也難以生成相對穩(wěn)定的水流場(如圖2所示的模型完善前流域出口水深和水速曲線)，從而不能精確模擬流域的徑流匯流過程。

為避免上述問題，在徑流匯流模型中添加水流混合模擬。水流混合模擬由柵格單元水速增量修正和鄰域柵格單元水流混合模擬兩步組成。

(1)柵格單元水速增量修正。假定水深大于1 m的柵格水流與周圍水團(tuán)混合后水速增量減小，并與水深呈反比關(guān)系。其修正公式為式(16)：

(16)

(2)鄰域柵格單元水流混合模擬。假定每一個(gè)時(shí)間片上，各柵格單元均與8鄰域柵格水流混合，混合過程中遵循動量守恒定律。首先計(jì)算中心柵格和8鄰域柵格的總動量和總水量，兩者相除得到中心柵格水流混合后的水速。

表2 模型完善前后典型柵格水速和水深變化Tab.2 The water speed and depth change of typical grid before and after the completion of runoff model

從表2可見，模型完善后，典型柵格的水流不斷累積，水深穩(wěn)定增加，水速逐漸減小，水流方向保持穩(wěn)定，符合徑流匯流的實(shí)際規(guī)律。

以岷江上游2018年5月20日12時(shí)鎮(zhèn)江關(guān)水文站瞬時(shí)流量作為岷江干流入流量(ξc=86.4)，在h=10、A=0.05的條件下,迭代15 h后，穩(wěn)定水流場生成，如圖3所示的模型完善后流域出口的水深和水速曲線。相對于模型完善前流域出口在上游穩(wěn)定供水條件下仍然有約0.08 m幅度的水位變化(圖2)，完善后的模型水深和水速均趨于穩(wěn)定，為基于穩(wěn)定流場的前期水情生成和徑流匯流模擬奠定了基礎(chǔ)。

圖2 模型完善前流域出口水深與水速變化Fig.2 The change of water depth and velocity at pour point before the model was modified

圖3 模型完善后流域出口水深與水速變化Fig.3 The change of water depth and velocity at pour point after the model was modified

3.3 基于穩(wěn)定流場的徑流匯流模擬

為客觀模擬流域前期水情，依據(jù)流域匯流時(shí)間，在上游河道輸入模擬時(shí)段前數(shù)小時(shí)的流量序列，通過模型迭代建立水流場作為流域前期水情。在模擬流域前期水情過程中，首先將前期流量序列的第一個(gè)小時(shí)流量代入模型作為上游穩(wěn)定供水，創(chuàng)建流域的穩(wěn)定流場；在此基礎(chǔ)上，將前期上游流量序列和面雨量序列逐小時(shí)迭代，模擬流域的前期水情。

以流域前期水情為基礎(chǔ)，將模擬時(shí)段逐小時(shí)上游流量序列和面雨量序列代入模型中迭代，模擬流域的徑流匯流過程。

為提高徑流匯流模擬精度，利用校場壩水文站水位數(shù)據(jù)對模型參數(shù)分別率定。在率定某一個(gè)參數(shù)時(shí)，其余參數(shù)設(shè)為固定值，取與實(shí)際水位誤差最小的參數(shù)值作為率定結(jié)果。率定結(jié)果如表3所示。

表3 模型參數(shù)率定結(jié)果Tab.3 The calibration results of model parameters

4 徑流匯流模擬精度評價(jià)

為客觀評價(jià)徑流匯流模擬的精度，將徑流匯流模擬按無穩(wěn)定流場無上游來水(A)、無穩(wěn)定流場有上游來水(B)、穩(wěn)定流場(C)3種方式，以2018年5月21日12時(shí)至2018年6月15日11時(shí)的長時(shí)間序列(600 h)為例，分別做出3種模擬方法的流域出口校場壩的水位曲線，并與實(shí)際水位曲線進(jìn)行對比，制作小時(shí)水位曲線(圖4)和誤差統(tǒng)計(jì)表(表4)。其中，A和B均未在模型中模擬水流混合；B和C均采用24 h前期上游流量和面雨量模擬前期水情。

從圖4可見，由于未考慮水流混合模擬，A和B模擬出的流域出口水位存在劇烈波動，與實(shí)際水位變化情況并不吻合，這與直接分解N-S方程產(chǎn)生的“匯流水流方向異常”有關(guān)；作為本文的徑流匯流模擬方法，C模擬的結(jié)果無水位劇烈波動，與實(shí)際水位曲線最接近。逐小時(shí)統(tǒng)計(jì)模擬水位與實(shí)際水位誤差，求取3種方法的平均水位誤差及與實(shí)際水位曲線的相關(guān)系數(shù)。從表4可見，基于穩(wěn)定流場的徑流匯流模擬的平均水位誤差最小，與實(shí)際水位曲線的相關(guān)系數(shù)達(dá)0.86，為高度正相關(guān)，其模擬結(jié)果具有一定的可信度，能在一定程度上指導(dǎo)暴雨洪澇災(zāi)害預(yù)測工作。

圖4 典型時(shí)段校場壩模擬水位與實(shí)際水位對比Fig.4 The comparation between simulated water level and actual water level of Jiaochangba

表4 3種徑流匯流方法的水位誤差統(tǒng)計(jì)
Tab.4 The water level error statistics of threerunoff simulation methods

模擬方法平均水位誤差/m相關(guān)系數(shù)A0.310.53B0.230.72C0.180.86

5 結(jié) 語

針對前期水情誤差影響徑流匯流模擬精度的問題，提出了一種基于穩(wěn)定流場的徑流匯流模型。首先探討了柵格模式下N-S方程的平流項(xiàng)、壓力項(xiàng)和外立項(xiàng)的具體實(shí)現(xiàn)，建立了初步的徑流匯流模型。實(shí)驗(yàn)發(fā)現(xiàn)，直接分解后得到的模型存在“匯流水流方向異?！眴栴}。通過在徑流匯流模擬中添加水流混合環(huán)節(jié)，不僅解決了該問題，也為穩(wěn)定流場創(chuàng)建奠定了基礎(chǔ)。然后，將徑流匯流模擬前的上游流量序列作為上游來水量，通過模型迭代模擬流域前期水情；最后，通過岷江上游鎮(zhèn)江關(guān)至校場壩流域徑流匯流模擬驗(yàn)證，發(fā)現(xiàn)本文方法更真實(shí)地反映了流域徑流匯流過程，模擬精度得到了提升，在流域暴雨洪澇預(yù)警預(yù)測中具有一定的應(yīng)用價(jià)值。

然而，模型只是原型的模擬，不可能反應(yīng)原型的一切特征，尤其像流域徑流匯流這樣復(fù)雜的現(xiàn)象。本文將上游入流作為流域前期的唯一水源，并未考慮流域冰川融水、地下水等產(chǎn)匯流因素，僅能得到近似的前期水情。另外，本文選取的實(shí)驗(yàn)流域無明顯的閘壩調(diào)蓄，人類活動影響較小。要進(jìn)一步提高流域徑流匯流模擬精度，需要在模型中進(jìn)一步考慮冰川融水、人類活動等因素。

□