一種基于穩(wěn)定流場的流域徑流匯流模擬方法

陳 軍，劉 意，李 婷

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

0 引 言

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

然而，流域匯流是一個非常復(fù)雜的過程，除降水之外，還存在上游河道來水、冰川融水、地下水補給等。由于上述因素難以精確地在流域上按網(wǎng)格單元量化，現(xiàn)有徑流匯流模型采用近似的方法進行模擬。以SWAT模型為例，它將多年的降雨數(shù)據(jù)代入到無水流域中反復(fù)迭代，將模型的最終輸出作為流域前期水情。這種方式不僅運算量大，還存在模擬誤差的迭代問題，因而所得到的前期水情與實際存在較大差距[7-11]；作為當前著名的水動力模型之一，F(xiàn)loodArea由于能提供流域任一柵格單元的水深模擬結(jié)果，被廣泛應(yīng)用于暴雨洪澇模擬、山洪風(fēng)險分析[12-14]，但該模型在徑流模擬時，前期水情僅作為輸入?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)以上春季有極少量的雪水補入，年內(nèi)徑流季節(jié)變化與降雨季節(jié)相應(yīng)。岷江上游多年平均年徑流量松潘縣鎮(zhèn)江水文站約為18 億m3，匯口下不遠的五通橋水文站多年平均年徑流量約780 億m3。區(qū)域大面積暴雨發(fā)生或隨洪峰傳播，干流峰高量大，嚴重威脅中下游安全[15]。為降低人類活動對徑流匯流模型精度驗證的影響，選取岷江上游的鎮(zhèn)江關(guān)和校場壩兩個水文站，將其所在河道位置分別作為流域出口，計算流域重疊范圍作為研究流域，如圖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)典運動方程，反映了黏性流體(又稱真實流體)流動的基本力學(xué)規(guī)律，在流體力學(xué)中有十分重要的意義[16]。但它是一個非線性偏微分方程，只有在極少數(shù)簡單流動的情況下才有解析解[17-20]。本文以N-S為指導(dǎo)，探索柵格模式下的徑流匯流模擬。

2.1 水流運動的N-S方程

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

(1)

▽μ=0

(2)

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

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

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

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

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

Δt=tn+1-tn

(3)

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

(4)

(5)

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

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

(6)

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

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

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

(7)

(8)

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

(9)

式中：C為柵格分辨率。

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

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

(10)

(11)

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

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

(12)

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

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

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

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

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

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

(13)

式中：ξc為時間片上柵格單元分攤的上游來水量。

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

ξ=∑ξcN

(14)

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

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

(15)

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

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

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

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

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

(16)

(2)鄰域柵格單元水流混合模擬。假定每一個時間片上，各柵格單元均與8鄰域柵格水流混合，混合過程中遵循動量守恒定律。首先計算中心柵格和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)定，符合徑流匯流的實際規(guī)律。

以岷江上游2018年5月20日12時鎮(zhèn)江關(guān)水文站瞬時流量作為岷江干流入流量(ξ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ù)小時的流量序列，通過模型迭代建立水流場作為流域前期水情。在模擬流域前期水情過程中，首先將前期流量序列的第一個小時流量代入模型作為上游穩(wěn)定供水，創(chuàng)建流域的穩(wěn)定流場；在此基礎(chǔ)上，將前期上游流量序列和面雨量序列逐小時迭代，模擬流域的前期水情。

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

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

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

4 徑流匯流模擬精度評價

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

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

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

表4 3種徑流匯流方法的水位誤差統(tǒng)計
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àn)，建立了初步的徑流匯流模型。實驗發(fā)現(xiàn)，直接分解后得到的模型存在“匯流水流方向異?！眴栴}。通過在徑流匯流模擬中添加水流混合環(huán)節(jié)，不僅解決了該問題，也為穩(wěn)定流場創(chuàng)建奠定了基礎(chǔ)。然后，將徑流匯流模擬前的上游流量序列作為上游來水量，通過模型迭代模擬流域前期水情；最后，通過岷江上游鎮(zhèn)江關(guān)至校場壩流域徑流匯流模擬驗證，發(fā)現(xiàn)本文方法更真實地反映了流域徑流匯流過程，模擬精度得到了提升，在流域暴雨洪澇預(yù)警預(yù)測中具有一定的應(yīng)用價值。

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

□