淺談MIKE 11水動(dòng)力學(xué)模型的搭建

姚 建

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

1 基本情況

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

2 模型搭建

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

(1)

式中：x、t分別為計(jì)算點(diǎn)空間和時(shí)間的坐標(biāo)；A為過水?dāng)嗝婷娣e；Q為過流流量；h為水位；q為旁側(cè)入流流量；C為謝才系數(shù)；R為水力半徑；α為動(dòng)量校正系數(shù)；g為重力加速度。

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

圖1 Abbott格式水位點(diǎn)、流量點(diǎn)交替布置圖

2.1 連續(xù)性方程求解

對(duì)每一h點(diǎn)求解連續(xù)性方程。h點(diǎn)處過流寬度bs可以描述為：

(2)

則連續(xù)方程可以寫為：

(3)

這里空間步長上，只有對(duì)Q求導(dǎo)，如圖2所示，則在時(shí)間步長n+1/2時(shí)，空間步長對(duì)Q的導(dǎo)數(shù)為：

(4)

而bs又可以寫為：

(5)

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

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

(6)

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

圖2 6點(diǎn)Abbott格式求解連續(xù)性方程

2.2 量方程的求解

對(duì)每一個(gè)q點(diǎn)求解動(dòng)量方程，如圖3所示。

圖3 6點(diǎn)Abbott格式求解動(dòng)量方程

通過數(shù)值變換，動(dòng)量方程可以寫為：

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

(7)

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

α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河流水動(dòng)力學(xué)模型的搭建主要包括4個(gè)文件：河網(wǎng)文件(River Network)、斷面文件(Cross-section)、邊界文件(Boundary Conditions)和水動(dòng)力參數(shù)文件(HD Parameters)的搭建。現(xiàn)分述如下：

1)河網(wǎng)文件：

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

圖4 模型河網(wǎng)圖

表1 河道一覽表

2)斷面文件：

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

表2 斷面位置一覽表

3)邊界文件：

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

表3 MIKE 11水動(dòng)力學(xué)模型開邊界設(shè)置情況

4)水動(dòng)力參數(shù)文件：

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

3 模型參數(shù)識(shí)別和驗(yàn)證

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

3.1 模型識(shí)別

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

表4 模型參數(shù)表

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

2)模型識(shí)別階段典型地下水站點(diǎn)實(shí)測水位與模擬水位對(duì)比關(guān)系見圖5。

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

(a) (b)

其中：ME—平均誤差，MAE—平均絕對(duì)誤差，RMSE—均方根差，STDres—?dú)埐畹臉?biāo)準(zhǔn)偏差，R—相關(guān)系數(shù)

3.2 模型驗(yàn)證

本模型以2001年1月1日地下水位作為識(shí)別的初始水位，以2005年至2010年的實(shí)測水位作為對(duì)比水位來進(jìn)行模型的驗(yàn)證。模型驗(yàn)證階段典型地下水站點(diǎn)實(shí)測水位與模擬水位對(duì)比關(guān)系見圖7。

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

4 模型率定結(jié)果分析

4.1 模型成果分析

模型率定階段共取實(shí)測觀測點(diǎn)3個(gè)，識(shí)別期和驗(yàn)證期的誤差詳見表5。

表5 模型率定階段誤差表

從表5可以看出：識(shí)別期平均誤差ME的平均值為-0.17m，平均絕對(duì)誤差MAE的平均值為0.34m，相關(guān)系數(shù)R的平均值0.79；驗(yàn)證期平均誤差ME的平均值為-0.11m，平均絕對(duì)誤差MAE的平均值為0.31m，相關(guān)系數(shù)R的平均值為0.79。模型識(shí)別和檢驗(yàn)結(jié)果證明所建立的數(shù)學(xué)模型、邊界條件及調(diào)定的水文地質(zhì)參數(shù)和源匯項(xiàng)基本都是正確可靠的。

4.2 存在問題分析

從識(shí)別期和驗(yàn)證期文屯地下水監(jiān)測井實(shí)測水位與模擬水位對(duì)比圖中不難發(fā)現(xiàn)雖然總體趨勢較好，但模擬水位過程線與實(shí)測水位過程線在年內(nèi)的擬合不是太好，模擬出來的線型總體比實(shí)際的要超前幾個(gè)月。其主要原因是：研究區(qū)范圍內(nèi)巖性具粗細(xì)相間的結(jié)構(gòu)特點(diǎn)，有多層亞黏土、砂、砂礫石互層，而我們?cè)谀Ｐ椭性O(shè)的地質(zhì)層由于資料的限制只是概化成均勻的兩層，這樣就不可避免對(duì)模擬產(chǎn)生影響，模擬出來的結(jié)果比實(shí)際的要來的敏感。

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

雖然在模擬的結(jié)果中有以上一些問題，但由于線型及趨勢總體較好，各個(gè)誤差統(tǒng)計(jì)參數(shù)比較理想，所以這些問題對(duì)水平衡及地下水總量的控制不會(huì)產(chǎn)生大的影響。