随着数值模拟技术的发展，一维[1-3]及二维[4-6]数学模型被广泛应用于河流、海岸、湖泊等的水动力、水质问题研究中。一维模型具有计算效率高、能够构建复杂河网和便于设置水工建筑物调度规则的优势。二维模型具有计算网格精细，能够精确模拟局部流态和污染物质运移规律的优势。近年来将一、二维模型耦合起来发挥其各自优势，形成一二维耦合模型成为一种必然的趋势[7]。例如崔冬等[8]将一维节制闸和闸上下游二维河道耦合建模，研究了节制闸不同调度规则下，水利枢纽通航水流条件；顾杰等[9]利用一二维耦合模型研究了秦皇岛河流和海岸的水动力、水质变化规律；郭聪等[10]建立了茅洲河一维模型与口外海域二维模型相耦合的水动力、水质模型，在茅洲河水环境治理效果预测中发挥了重要作用。
本文以汾河中游二坝蓄水闸控河段和上游两个污水处理厂为研究对象，利用MIKE软件建立一二维耦合的水动力与水质数学模型；在闸前水位不变的前提下，分析不同水文条件下和污水厂发生非正常排污时，河道的水流特性和污染物的运移规律。
1 模型原理
1.1 一维模型
1.1.1水动力模型
水动力模块以圣维南方程组为基础建立，其包括了连续性方程和动量方程，分别如下：
(1)
(2)
式中：Q为流量；q为侧向入流；A为过水断面面积；h为水位；R为水力半径；C为谢才系数；α为动量修正系数。
1.1.2水质模型
对流扩散模型考虑水环境中污染物的对流扩散过程及污染物线性消减，其方程如下：
(3)
式中：C为污染物浓度；D为污染物弥散系数；A为断面过水面积；Q为流量；K为降解系数；C2污染物的电源浓度；q为污染物的点源流量。
1.2 二维模型
1.2.1水动力模型
水动力模型原理基于沿水深积分的二维不可压缩雷诺平均Navier-Stokes方程，服从Boussinesq假定和静水压力假定，其方程如下：
(4)
(5)
(6)
式中：η和d分别表示河底高程和静水深；
分别为x、y方向上的平均流速；g为重力加速度；ρ为水的密度；ρ0为水的相对密度；Pa为当地大气压强；Sij为辐射应力分量；τbx、τby为河床底部应力张量；S为点源流量；us、vs为点源在x、y方向上的水流流速；Tij为水平黏滞应力项，根据沿水深平均流速的梯度计算。
1.2.2水质模型
水质模型与水动力模型是动态连接的，可以模拟污染物随水流的运移规律，其方程如下：
(7)
式中：t为时间；C为浓度；D为水平扩散系数；u、v为流速。
1.3 一二维耦合模型
利用MIKE FLOOD的一二维标准链接模式，将一维和二维模型进行耦合建模。水动力耦合采用了双扫描计算方法[11]，在一、二维模型耦合线上，一维模型从二维模型中获得第n时步水深hn，计算出第n时步的流量Qn，一维模型再利用Qn及hn，应用式(7)预测出第n+1/2时步的流量Qn+1/2，并作为源项输入到二维模型，二维模型再利用Qn+1/2计算出第n+1时步的总水深hn+1，以此计算所有时间步的值。
(8)
水质耦合计算取决于一、二维模型间的水流方向，当水流由一维(1D)流向二维(2D)时，一维水质浓度计算结果作为源项输入到二维模型：
(9)
当水流由二维(2D)流向一维(1D)时，二维模型中扩散方程修正为：
(10)
2 模型应用
2.1 模拟区域概况
模拟河段为汾河太原城区段三期工程2#闸至介休市义棠水文站之间的河道，全长约110 km，平均河宽约300 m，同时考虑汾河的一条支流潇河，其上边界取潇河北格镇水闸控制断面。汾河太原城区段三期工程2#闸下游约12.5 km处为汾河二坝水闸控制工程，共设12孔弧形闸门，闸孔净宽10 m.水闸上游回水影响范围内有晋阳污水厂和汾东污水处理厂(以下简称A厂和B厂)，其中A厂位于汾河西岸，在二坝水闸上游约2 km处，B厂位于潇河北岸，距汾河东岸约4.5 km，模拟区域和工厂的相对位置见图1。
图1 模拟区域概况
Fig.1 Simulation area overview
2.2 模型概化
汾河二坝以上蓄水运行下回水区是水动力特性和A、B厂污染物运移规律研究的重点区域，且具有二维属性，故将汾河二坝水闸以上5.5 km河段用二维模型概化，图1所示的蓝色河段最大河宽600 m，最小140 m.采用三角形非结构化网格划分二维区域，并对污水厂排水口和水闸出口处网格局部加密，网格最小长度5 m，模型共划分网格21 425个，具体网格见图2.其余河段利用一维模型概化，共设置206个断面(见图1)，在51号断面下游汾河一维河段设置二坝水闸，方便调度控制。
图2 水闸回水区二维网格
Fig.2 Two-dimensional grid of backwater of sluice
在一、二维模型耦合建模过程中，需将二维模拟区域的上下边界层网格，与相邻一维河段边界断面相连，分别为23、50、51号断面，这样水流和浓度可在两模型之间传递。
2.3 边界条件
汾河与潇河上游边界给定流量和浓度变化过程；汾河下游水边界根据义堂水文站水位-流量关系给定，浓度设定为自由出流边界。污水厂排水口给定其排水流量和COD浓度值。
2.4 参数设置
水动力参数主要涉及曼宁系数和涡黏性系数，曼宁系数通过率定取32 m1/3/s；涡黏性系数采用Smagorinsky公式，取0.28.
水质参数主要涉及扩散系数和降解系数，一维模型扩散系数取10 m2/s，二维扩散系数与Smagorinsky系数相当[12]，取为0.28 m2/s；COD降解系数的取值采用类比法，借鉴马耀光等[13]对泾河COD降解系数研究成果，认为汾河和泾河河流大小和水环境状况具有相似性，取0.091/d.模型计算过程中时间步长取1 s.
2.5 模型验证
2.5.1模拟洪水演进过程的验证分析
本文选取了该河段2016年7月19日-2016年7月22日的水文数据对模型进行了验证计算，对比分析二坝水文站的实测值与模拟值如图3，4所示，其流量和水位的最大相对偏差分别为13.5%和12.7%，平均值分别为6.1%和5.3%.
图3 流量验证过程
Fig.3 Flow verification process
图4 水位验证过程
Fig.4 Water level verification process
2.5.2模拟恒定流动河道水面线的验证分析
针对计算情景一(见表1)的恒定水流情况，将耦合模型计算所得沿主河道的水面线与传统方法[14]计算结果进行对比分析，如图5所示。两者偏差在8%以内。
图5 水位沿程变化
Fig.5 Water level varies along the course
由于数值模拟、实测资料和传统方法均有不同程度的误差，综合分析认为耦合模型与其比较出现的偏差在合理范围内，模拟结果可信。
3 计算结果及分析
3.1 计算情景设置
3.1.1天然来水水量和水质条件
选取研究域生态基流的水量和典型洪水过程两种典型来水条件。生态基流汾河取为3.75 m3/s，潇河为0.22 m3/s.典型洪水过程选取2016年7月22日14点-2016年7月26日7点发生的洪水过程，具体见图6.汾河二坝水闸所在河段水功能区满足地表Ⅳ类水质要求，故设置其COD基础质量浓度为30 mg/L.
图6 典型洪水过程
Fig.6 Typical flood process
3.1.2污水厂水量和水质条件
根据污水厂处理规模，A厂的排水流量为5.56 m3/s，B厂排水流量为6.62 m3/s，经过处理后均排入汾河。正常排污质量浓度为30 mg/L，考虑到两个污水厂都存在非正常排污可能，设非正常排污浓度为200 mg/L，同时设置相应时间3 h.
3.1.3闸坝控制条件
为保证灌溉需要的水量，设置水闸调度原则为保持正常蓄水水位764.4 m不变，生态基流条件下开中部两孔闸和发生典型洪水条件下开中部六孔闸放水调节水位。
综合以上条件，共设置4种计算情景，具体描述见表1.
表1 模型计算情景
Table 1 Model calculation scenario
情景设置水流条件(流量)/(m3·s-1)A厂B厂汾河潇河水质条件非正常污水厂及排污时间/hA厂质量浓度/(mg·L-1)B厂质量浓度/(mg·L-1)水闸调度开闸孔数原则情景一5.566.623.750.22A厂/3200302情景二5.566.623.750.22B厂/3302002情景三5.566.62洪峰97洪峰10A厂/3200306情景四5.566.62洪峰97洪峰10B厂/3302006调节闸门开度,维持正常蓄水水位764.4m
3.2 计算结果与分析
3.2.1水动力计算结果与分析
情景一和情景二的水流及水闸调度条件相同，为生态基流下，闸门开2孔，二坝水闸回水区及闸附近流场见图7.由图可知：河道整体流速偏小，主流区平均在0.019 m/s左右，当水流通过闸门时流速最大，为0.21 m/s；在水闸左岸库区形成一个较大的死水区；A厂排水流速为0.18 m/s，其垂直河岸射流对河道局部流场产生了较大的影响，形成了一个顺时针的环流，其影响范围约为长300 m，宽270 m.
图7 情景一、二下水闸回水区及闸附近流场
Fig.7 Scenarios 1 and 2 lower sluice backwater area and flow field near sluice
情景三和情景四的水流及水闸调度条件相同，为典型洪水下，闸门开6孔，洪峰通过水闸典型时刻(t=14 h)的回水区及闸附近流场见图8.由图可知：水闸回水区的上游流速较大，主流区平均流速为0.47 m/s左右；中游流速较小，闸前1 000 m处主流区流速为0.15 m/s；水流通过闸门时流速最大，为0.59 m/s；水闸左岸库区流速最小，为0.005 m/s左右，形成较为封闭的水域；A厂排水对河道整体流场影响较小，这是由于上游洪峰来流量是污水厂排水流量的近20倍，使得A厂排水产生的垂直于河岸水流动量快速衰减。
图8 情景三、四下水闸回水区及闸附近流场(t=14 h)
Fig.8 Scenarios 3 and 4 backwater area of sluice and flow field around sluice (t=14 h)
3.2.2污染物出闸规律分析
情景一和情景二为河道生态基流下，A、B厂分别出现非正常排污情况，其闸后COD浓度随时间变化规律见图9.由图可知：A厂排污24.33 h后，污染物开始排出水闸，在第31.83 h达到峰值浓度50.65 mg/L，在52.00 h后污染物完全排出水闸。B厂排污44.50 h后，污染物开始排出水闸，在第54.33 h达到峰值浓度39.92 mg/L，在72.33 h后完全排出水闸。B厂比A厂距二坝远5.5 km，B厂污染物延后20.17 h开始排出水闸，峰值浓度低10.73 mg/L.
图9 情景一、二下闸后COD浓度随时间变化
Fig.9 COD concentration changes with time after the sluice of the first and second scenario
情景三和情景四为河道典型洪水下，A、B厂分别出现非正常排污情况，其闸后COD浓度随时间变化规律见图10.由图可知：A厂排污4.50 h后，污染物开始排出水闸，在第8.67 h达到峰值浓度36.46 mg/L，在12.17 h后污染物完全排出水闸。B厂污染物9.34 h后排出水闸，在13.50 h后达到峰值浓度36.60 mg/L，在20.50 h后完全排出水闸。B厂污染物比A厂污染物延后4.84 h开始排出水闸，峰值浓度高0.14 mg/L.
图10 情景三、四下闸后COD浓度随时间变化
Fig.10 COD concentration changes with time after the sluice of the third and fourth scenario
与情景一、二生态基流情况相比，情景三、四典型洪水情况下，A、B厂污染物开始排出水闸时间分别提前了19.83 h、35.16 h，其峰值浓度分别下降了14.19 mg/L、3.32 mg/L，说明河道来流量大更利于污染物快速排出水闸和消减其浓度。
3.2.3回水区内污染物运移规律分析
对于四种情景，分别取前期、中期和末期三个典型时刻，分析浓度在汾河二坝回水影响区内的迁移变化规律。各计算情景动态模拟结果见图11-图14.
图11 典型时刻二坝闸前COD浓度分布图(情景一)
Fig.11 COD concentration distribution diagram in front of the second dam at a typical time (scenario 1)
图12 典型时刻二坝闸前COD浓度分布图(情景二)
Fig.12 COD concentration distribution diagram in front of the second dam at a typical time (scenario 2)
图13 典型时刻二坝闸前COD浓度分布图(情景三)
Fig.13 COD concentration distribution diagram in front of the second dam at a typical time (scenario 3)
图14 典型时刻二坝闸前COD浓度分布图(情景四)
Fig.14 COD concentration distribution diagram in front of the second dam at a typical time (scenario 4)
情景一：A厂污染物前期呈团状，在向下游移动过程中，前部污染团浓度高，后部形成浓度较低的污染带。前、中、末期受污染水域面积(质量浓度大于30 mg/L的范围，下同)分别为：0.16、0.29、0.36 km2.
情景二：B厂污染物主要呈团状向下游运移，在向下游运移过程中污染团面积逐步变大，高浓度污染物主要集中在污染团中部。前、中、末期受污染水域面积分别为：0.17、0.24、0.28 km2.
情景三：A厂污染物主要呈带状向下游运移，高浓度污染物贴近右岸。前、中、末期受污染水域面积分别为：0.28、0.31、0.26 km2.
情景四：B厂污染物主要呈带状向下游运移，并逐步扩散至整个河道断面，高浓度污染物贴近河道右岸。前、中、末期受污染水域面积分别为0.36、0.65、0.70 km2.
生态基流情况下，污染物以高浓度污染点为中心呈辐射状扩散，主要以团状缓慢下移，这是由于扩散起了主导作用；A厂排污比B厂排污的三个时期平均污染水域面积大20.83%.典型洪水情况下，高浓度污染物主要贴近河岸分布，浓度向河道中部方向逐步递减，污染物以带状快速下移，这是由于对流起了主导作用；A厂排污比B厂排污的三个时期平均污染水域面积小64.05%.对比各情景下的受污染水域面积发现，情景四末期受污染水域面积最大。
4 结论
1) 本文采用MIKE软件，建立了汾河二坝水闸回水区二维模型与回水区上下游和潇河河段一维模型相耦合的水动力和水质数学模型。并利用实测水文数据对水动力学模型进行了验证计算，验证结果良好。
2) 分析各情景回水区流场，水流通过闸门时流速最大；A厂排水对二坝回水区的流场影响较大，特别是上游来水较小时。情景一、二水流通过闸门流速为0.21 m/s，A厂排水使河道形成了一个影响范围在270～300 m的顺时针环流；情景三、四洪峰通过时段闸孔出流流速为0.59 m/s，此刻，A厂排水对河道整体流场几乎无影响。
3) 河道来流量大有利于污染物快速排出二坝回水区域，消减污染物浓度。与情景一、二生态基流来流情况相比，情景三、四典型洪水情况下，A、B厂污染物开始排出闸孔时间分别提前了19.83 h和35.16 h，其峰值浓度分别下降了14.19 mg/L和3.32 mg/L.
4) 相同水流条件下，距二坝较远的B厂污染物在二坝回水区的滞留时间较长。生态基流情况下，B厂污染物比A厂污染物开始排出闸孔时间延后20.13 h；典型洪水情况下，B厂污染物比A厂污染物开始排出水闸时间延后4.84 h.
5) 河道来流量大小对污染物的对流、扩散有较大影响，生态基流来流情况下，污染物以团状缓慢下移，扩散起了主导作用；典型洪水来流情况下，污染物主要贴近河岸以带状快速下移，对流起了主导作用。
汾河由于受人类活动的影响，已经由天然河流演变为人工改造型河流，为了使河流水环境模型能够模拟受闸控河段的水动力特性及污染物运移规律，本文引入了一二维耦合的水动力水质模型进行研究。这样可以充分利用一、二维模型各自具有的优势，使所研究问题更能反映实际，同时研究成果对汾河二坝回水区污染物预测和控制有一定的参考价值。
[1] 马小雪,杨军,毛媛媛,等.平原河网区突发性水污染事件应急调水数值模拟分析[J].中国农村水利水电,2015(4)：47-50,59.
MA X X,YANG J,MAO Y Y,et al.A numerical simulation and analysis of emergency water transfer for sudden sater pollution accident in the Plain River[J].China Rural Water and Hydropower,2015(4)：47-50,59.
[2] 沈红丽,高成.基于MIKE11和AHP法的感潮河段防洪潮方案研究[J].三峡大学学报(自然科学版),2020,42(1)：36-41.
SHEN H L,GAO C.Research on flood tide control planning of tidal river based on MIKE11 Model and AHP[J].Journal of China Three Gorges University(Natural Sciences),2020,42(1)：36-41.
[3] 熊鸿斌,陈雪,张斯思.基于MIKE11模型提高污染河流水质改善效果的方法[J].环境科学,2017,38(12)：5063-5073.
XIONG H B,CHEN X,ZHANG S S.Method of improving the water quality of polluted rivers based on the MIKE11 Model[J].Environmental Science,2017,38(12)：5063-5073.
[4] 贾瑞鹏,窦明,米庆彬,等.基于MIKE21的万宝湖二维水环境数值模拟[J].环境污染与防治,2018,40(5)：527-532.
JIA R P,DOU M,MI Q B,et al.Two-dimensional simulation on water environment of Wanbao Lake based on MIKE21 software[J].Environmental Pollution and Control,2018,40(5)：527-532.
[5] RU ZHONG,SHIGEKI MASUNAGA,HONG TIANQIU,et al.Fuzzy model for two-dimensional river water quality simulation under sudden pollutants discharged [J].Journal of Hydrodynamics,2007,19(4)：434-441.
[6] ZHU C J,LIANG Q,YAN F,et al.Reduction of waste water in erhai lake based on MIKE21 hydrodynamic and water quality model[J].The Scientific World Journal,2013(7/8)：958506.
[7] 马天海,孙娟,颜剑波.基于MIKEFLOOD阳澄湖一二维水动力耦合模型研究[J].水科学与工程技术,2018(1)：25-30.
MA T H,SUN J,YAN J B.Research on 1D and 2D hydrodynamic coupling model of Yangcheng Lake based on MIKE FLOOD model[J].Water Sciences and Engineering Technology,2018(1)：25-30.
[8] 崔冬,刘新成,潘丽红.一、二维耦合数学模型在水利枢纽通航水流条件研究中的应用[J].水利水电科技进展,2009,29(1)：35-39.
CUI D,LIU X C,PAN L H.Application of 1D and 2D coupled mathematical model in navigation flow condition of hydro-junction[J].Advances in Science and Technology of Water Resources,2009,29(1)：35-39.
[9] 顾杰,胡成飞,李正尧,等.秦皇岛河流—海岸水动力和水质耦合模拟分析[J].海洋科学,2017,41(2)：1-11.
GU J,HU C F,LI Z Y,et al.Coupling simulation and analysis of hydrodynamics and water quality in Qinhuangdao rivers and coastal waters[J].Marine Sciences,2017,41(2)：1-11.
[10] 郭聪,施家月,张亚力,等.一、二维耦合模型在茅洲河水环境治理中的应用[J].环境影响评价,2019,41(4)：59-62.
GUO C,SHI J Y,ZHANG Y L,et al.Application of 1D-2D coupling model in the water environment management for maozhou River[J].Environmental Impact Assessment,2019,41(4)：59-62.
[11] DHI.MIKE FLOOD user manual[M].Copenhagen,Denmark：MIKE BY DHI,2014.
[12] 徐聪.引调水工程水环境风险与水质安全保障研究[D].西安：西安理工大学,2016：36.
[13] 马耀光,李书琴,杨宣政,等.泾河下游可降解污染物的自净特征[J].干旱地区农业研究,2003,21(1)：126-128.
MA Y G,LI S Q,YANG X Z,et al.The degradation characteristic of organic pollutants in the Jinghe River[J].Agricultural Research in the Arid Areas,2003,21(1)：126-128.
[14] 张志昌,魏炳乾,郝瑞霞.水力学(下册)[M].北京：中国水利水电出版社,2011：21-25.