浅谈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年的大断面资料加上这两个位置的高程常数差值作为此处的断面。出流断面设在河流入海口上游一段距离附近，断面也同样借用沙河水文站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 存在问题分析

从识别期和验证期文屯地下水监测井实测水位与模拟水位对比图中不难发现虽然总体趋势较好，但模拟水位过程线与实测水位过程线在年内的拟合不是太好，模拟出来的线型总体比实际的要超前几个月。其主要原因是：研究区范围内岩性具粗细相间的结构特点，有多层亚黏土、砂、砂砾石互层，而我们在模型中设的地质层由于资料的限制只是概化成均匀的两层，这样就不可避免对模拟产生影响，模拟出来的结果比实际的要来的敏感。

而大房身监测井实测水位与模拟水位对比图是这三个地下水位监测井中最差的一个，究其原因除了有刚才文屯的因素，另外一个重要原因是大房身监测井不是专用监测井，它平时有抽用水情况。特别是在天旱少雨的年份，它会被大量的抽水浇灌葡萄等作物，所以测出的地下水位受人为影响较大。虽然我们在设置葡萄灌溉区域及灌溉量时加以了考虑，但由于实际情况随机而复杂，所以模拟出来的单井效果并没有很大改善。

虽然在模拟的结果中有以上一些问题，但由于线型及趋势总体较好，各个误差统计参数比较理想，所以这些问题对水平衡及地下水总量的控制不会产生大的影响。