淺談MIKE 11水動力學模型的搭建

姚 建

(遼寧省營口水文局，遼寧 營口 115000)

1 基本情況

文中模型以蓋州市-鲅魚圈區山前平原為研究區，西邊為渤海，面積大約為300km2。本模型研究區域包含3個地下水監測井：文屯、大房身和蘭東，兩條主要河流熊岳河、沙河。

2 模型搭建

MIKE 11水動力計算模型基于垂向積分的物質和動量守恒方程，即一維非恒定流Saint-Venant方程組來模擬河流或河口的水流狀態[1]。

(1)

式中：x、t分別為計算點空間和時間的坐標；A為過水斷面面積；Q為過流流量；h為水位；q為旁側入流流量；C為謝才系數；R為水力半徑；α為動量校正系數；g為重力加速度。

方程組利用Abbott-Ionescu六點隱式有限差分格式求解，如圖1所示。該格式在每一個網格點不是同時計算水位和流量，而是按順序交替計算水位或流量，分別稱為h點和Q點。Abbott-Ionescu格式具有穩定性好、計算精度高的特點。離散后的線形方程組用追趕法求解。

圖1 Abbott格式水位點、流量點交替布置圖

2.1 連續性方程求解

對每一h點求解連續性方程。h點處過流寬度bs可以描述為：

(2)

則連續方程可以寫為：

(3)

這里空間步長上，只有對Q求導，如圖2所示，則在時間步長n+1/2時，空間步長對Q的導數為：

(4)

而bs又可以寫為：

(5)

式中：Aoj為計算點j-1和j之間的面積；Aoj+1為計算點j和j+1之間的面積；Δ2xj為計算點j-1和j+1之間的空間步長。將以上各式代入連續性方程得出：

αjQj-1n+1+βjhjn+1+γjQj+1n+1=δj

(6)

式中α，β，γ是b和δ的函數，并隨n時刻Q和h及n+1/2時刻Q的大小而變化。

圖2 6點Abbott格式求解連續性方程

2.2 量方程的求解

對每一個q點求解動量方程，如圖3所示。

圖3 6點Abbott格式求解動量方程

通過數值變換，動量方程可以寫為：

αjhj-1n+1+βjQjn+1+γjhj+1n+1

(7)

式中：(各參數符合意義同上)

αj=f(A)

βj=f(Qjn，△t,△x,C,A,R)

γj=f(A)

δj=f(A，△x，△t，α,q,v,θ,hj-1n，Qj-1n+1/2，Qjn，hj+1n，Qj+1n+1/2)

Mike11河流水動力學模型的搭建主要包括4個文件：河網文件(River Network)、斷面文件(Cross-section)、邊界文件(Boundary Conditions)和水動力參數文件(HD Parameters)的搭建。現分述如下：

1)河網文件：

模型研究區域較小，研究區內主要河流只有熊岳河和沙河，河網清晰。這樣模型的河網文件就比較簡單明了，詳見圖4模型河網圖、表1河道一覽表。

圖4 模型河網圖

表1 河道一覽表

2)斷面文件：

在熊岳河進入模型研究區不遠處便是熊岳水文站。因此模型采用了熊岳水文站2010年的大斷面資料作為河流的入流斷面。熊岳河由東南往西北穿過模型區域注入渤海，在入海口處設出流斷面。由于在該處沒有實測斷面資料，所以在處理時我們先從地形圖上量出熊岳河入?？诘拇笾聦挾?，再利用arcgis在數字高程圖中提取熊岳河入海口處的若干高程值，以此概化出河流入海口處的斷面情況。并且由于同一條河段河底坡度變化有的急劇，有的平緩，所以在河底坡度有顯著變化的位置也需要假設斷面，才能更好地描述河底高程的變化情況，為此在河段河底比降急劇變化(36300米)處假設了一個斷面，此斷面是由模型自動根據上下兩個斷面插補而來。沙河的入流斷面設在兩個支流的交界處，由于此處沒有水文站，故借用了其下游不遠處的沙河水文站2010年的大斷面資料加上這兩個位置的高程常數差值作為此處的斷面。出流斷面設在河流入?？谏嫌我欢尉嚯x附近，斷面也同樣借用沙河水文站2010年的大斷面資料減去這兩個位置的高程常數差值作為此處的斷面。詳見圖4(圖中紅色方框即為斷面)和表2。

表2 斷面位置一覽表

3)邊界文件：

MIKE 11水動力學模型在上游入流處和下游出口處都需要設置邊界條件，一般將水文站上游實測流量和下游實測水位作為模型的開邊界。本模型中熊岳河上游入流處采用熊岳水文站2001-2010年實測流量資料作為流量開邊界。下游出口處，由于沒有水文站或實測水位資料，所以只能借用營口潮水位站的水位資料作為水位開邊界。沙河上游入流處沒有實測流量資料，但考慮到其與沙河水文站距離不遠，且中間沒有水流匯入或分出，故采用沙河水文站2001-2010年實測流量資料作為流量開邊界。下游出口處設了一個水位開邊界，其設置方法同熊岳河。

表3 MIKE 11水動力學模型開邊界設置情況

4)水動力參數文件：

水動力參數文件的設置主要包括初始條件和河床糙率的設定。給定初始條件的目的是讓模型平穩啟動，所以原則上初始水位(或水深)和流量應盡可能與模擬開始時刻的實際河網水動力條件一致。初始水位不能高于河岸或低于河床，否則可能導致模型無法順利起算。如果有計算得到的結果文件，可以將其引入成熱啟動文件作為模型計算的初始條件。本模型中，根據2001年1月1日的水文資料，熊岳河設置初始水深0.20m，初始流量0.100m3/s；沙河設置初始水深0.18m，初始流量0 m3/s(死水)。另一個需要設置的參數是河床糙率。由于缺少實測糙率資料，所以模型中糙率M統一設成經驗值30。

3 模型參數識別和驗證

由于河流入?？诟浇鼪]有水文站，缺乏實測水位和流量資料，所以沒辦法單獨對MIKE 11模型進行識別和驗證，這無疑會對模型總體的精度產生一定影響，但后面我們把MIKE 11與MIKE SHE進行了耦合，并對整個模型進行了綜合率定，效果良好。

3.1 模型識別

為使模型具有足夠的精度，必須進行模型識別。本模型以2001年1月1日地下水位作為識別的初始水位，以2001年至2005年的實測水位作為對比水位來進行模型的識別。在這個階段，要調整輸入量和邊界性質，以驗證模型的適用性和所求參數的準確性。模型最終確定參數見表4。

表4 模型參數表

1)從模型識別階段2001～2005年中選取了2005年8月1日，2004年4月1日及2002年12月1日豐、平、枯三個時段的模擬地下水流場。研究區的地下水走向大致是從山丘區流向海邊，并且地下水位在豐水期較高，在平水期一般，在枯水期較低，而且靠近山丘的區域由于含水層厚度很薄，地下水流量較小，這些與實際情況是相吻合的。

2)模型識別階段典型地下水站點實測水位與模擬水位對比關系見圖5。

從圖5可以看出：識別階段平均誤差最大值為文屯站的-0.319908，最小值為蘭東站的0.0631603；平均絕對誤差的最大值為大房身站的0.448679，最小值為蘭東站的0.178977。相關關系R普遍達到0.80以上，只有大房身站較低，為0.65。根據以上統計數據來看，模型的率定效果較好。

(a) (b)

其中：ME—平均誤差，MAE—平均絕對誤差，RMSE—均方根差，STDres—殘差的標準偏差，R—相關系數

3.2 模型驗證

本模型以2001年1月1日地下水位作為識別的初始水位，以2005年至2010年的實測水位作為對比水位來進行模型的驗證。模型驗證階段典型地下水站點實測水位與模擬水位對比關系見圖7。

從圖6可以看出：驗證階段平均誤差最大值為文屯站的-0.283844，最小值為蘭東站的0.0444185；平均絕對誤差的最大值為文屯站的0.411728，最小值為蘭東站的0.16355。相關關系R普遍達到0.8以上，只有大房身站較低。

4 模型率定結果分析

4.1 模型成果分析

模型率定階段共取實測觀測點3個，識別期和驗證期的誤差詳見表5。

表5 模型率定階段誤差表

從表5可以看出：識別期平均誤差ME的平均值為-0.17m，平均絕對誤差MAE的平均值為0.34m，相關系數R的平均值0.79；驗證期平均誤差ME的平均值為-0.11m，平均絕對誤差MAE的平均值為0.31m，相關系數R的平均值為0.79。模型識別和檢驗結果證明所建立的數學模型、邊界條件及調定的水文地質參數和源匯項基本都是正確可靠的。

4.2 存在問題分析

從識別期和驗證期文屯地下水監測井實測水位與模擬水位對比圖中不難發現雖然總體趨勢較好，但模擬水位過程線與實測水位過程線在年內的擬合不是太好，模擬出來的線型總體比實際的要超前幾個月。其主要原因是：研究區范圍內巖性具粗細相間的結構特點，有多層亞黏土、砂、砂礫石互層，而我們在模型中設的地質層由于資料的限制只是概化成均勻的兩層，這樣就不可避免對模擬產生影響，模擬出來的結果比實際的要來的敏感。

而大房身監測井實測水位與模擬水位對比圖是這三個地下水位監測井中最差的一個，究其原因除了有剛才文屯的因素，另外一個重要原因是大房身監測井不是專用監測井，它平時有抽用水情況。特別是在天旱少雨的年份，它會被大量的抽水澆灌葡萄等作物，所以測出的地下水位受人為影響較大。雖然我們在設置葡萄灌溉區域及灌溉量時加以了考慮，但由于實際情況隨機而復雜，所以模擬出來的單井效果并沒有很大改善。

雖然在模擬的結果中有以上一些問題，但由于線型及趨勢總體較好，各個誤差統計參數比較理想，所以這些問題對水平衡及地下水總量的控制不會產生大的影響。