姚 建

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

1 基本情況

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

2 模型搭建

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

(1)

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

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

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

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

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

(2)

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

(3)

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

(4)

而bs又可以寫為：

(5)

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

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

(6)

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

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

2.2 量方程的求解

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

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

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

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

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

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

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

表1 河道一覽表

2)斷面文件：

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

表2 斷面位置一覽表

3)邊界文件：

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

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

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

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

3 模型參數(shù)識別和驗證

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

3.1 模型識別

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

表4 模型參數(shù)表

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

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

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

(a) (b)

其中：ME—平均誤差，MAE—平均絕對誤差，RMSE—均方根差，STDres—殘差的標準偏差，R—相關(guān)系數(shù)

3.2 模型驗證

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

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

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

4.1 模型成果分析

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

表5 模型率定階段誤差表

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

4.2 存在問題分析

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

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

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