大型浅水湖泊三维水环境数值模拟及关键参数研究

王海燕， 陈凯麒， 祁昌军，张贝贝， 黄 茹， 温静雅

(1:环境保护部环境工程评估中心，北京 100012;2．国家环境保护环境影响评价数值模拟重点实验室，北京 100012;3.水电环境研究院，北京 100012; 4.河海大学水文水资源学院,南京 210098)

1 引 言

当前，数值模拟已经成为研究水环境问题的一项重要手段，以往针对浅水湖泊的研究，通常采用二维数值模拟以减少计算量，同时因实测数据难以获取而无法进行更深层次的机理性研究。随着计算机水平发展和数据获取手段的丰富，水环境数值模拟的精度亟需进一步提升。对于大型浅水湖泊来说，研究域水平空间尺度远远大于垂向尺度，如太湖南北长68.5km，东西平均宽34km，而湖泊平均水深只有1.9m，平面和垂向尺度相差万倍。根据2003年7月秦伯强[1]等在太湖乌龟山附近开展的湖流垂直结构观察发现，即使是太湖这样最大水深不足3m的浅水湖泊，也存在垂向湖流分层的现象，且溶解氧、叶绿素及营养盐浓度等要素也与水深有着密切的响应关系[2]，为了进一步研究浅水湖泊流场内部结构，如补偿流和上、下流矢的切变等，更精确的描述富营养化的形成及转化过程，势必要更多的采用三维数值模拟，建立浅水湖泊生态水动力及水质模型，模拟水动力及水质要素迁移转化过程。

本文从水动力以及与垂向流场分布相关的水环境要素角度出发，总结、概括并分析了不同水质要素垂向分布规律及与水深响应关系，归纳了建模过程中关键参数的取值，并提出下一步浅水湖泊水环境模拟值得进一步深入研究的方向。

2 浅水湖泊三维水环境数值模拟研究

水环境数值模拟研究始于1925年，至今已有92年发展历程，美国、英国、韩国、澳大利亚、丹麦等国家曾先后发布了相应的环境模型开发及应用指南，并推荐了多种不同类型的数值模拟软件。我国大陆地区的数值模拟研究始于20世纪70年代，在水环境数值模拟研究方面也已开展了大量工作[3]，其中针对湖泊的数值模拟已取得很多成果[4-5]。早期二维水环境数值模拟被广泛用于浅水湖泊研究，但研究与水深有相关关系的要素时，则更多地采用三维数值模拟。

2.1 浅水湖泊三维风生流结构

湖泊流场主要有三种：风生流、密度流和吞吐流，其中浅水湖泊中风生流占主导地位。在定常风作用下，浅水湖泊表层为风生流，流向与风向平行，中间为过渡层，底层流向相反，为补偿流[6]。在大小相同，方向相反的定场风作用下，湖流在水平方向上呈现流型相似，但环流方向相反的特征[7]，且风速主要影响环流流速，风向影响环流结构[6, 8～10]，在垂向上，增水区域下沉运动强烈，减水区域存在较强的上升补偿流，形成垂向封闭环流[11]。此外，当风生流主导湖泊流场时，存在临界风速值，超过临界值的风速会使流场结构发生显著改变[12]。关于盐度对艾比湖[13]、美国切萨皮克湾[14]三维流场影响的研究结果表明，咸水湖与淡水湖水动力空间分布相似，水体盐度增加对湖泊环流结构基本无影响。黄平[15]建立了基于交错网格差分方法的湖泊风生流三维数值模型，成功运用于武汉市墨水湖；曹慧江等[9]用三维数值模型计算了夏季滴水湖吞吐流、密度流和风生流三者共同作用下的流场；汤露露等[16-17]建立了基于ECOMSED模式的太湖风生流三维数值模型，其计算的太湖三维水流结构与现有成果基本吻合[8, 18]，田勇[19-20]提出了能应对实际工程计算需要的湖泊三维水动力计算模型，该模型能够有效模拟湖泊三维流场结构以及风生环流、补偿流、吞吐流等复杂流速，揭示不同分层流场的流速差异；张利民、濮培民等建立了日本琵琶湖[21～23]水动力模型，基本掌握了其水动力过程及机制；Schmalz[24]建立了Okeechobee湖水位模型。

2.2 沉积物及水体污染物和蓝藻空间分布

风场可以影响水动力过程，使湖泊局部水位产生变化，继而改变水体运动速度和方向，影响各种物质在湖泊内的输移扩散，对湖泊的水质及生态系统的结构和功能产生重要影响[12]。水动力引起的沉积物再悬浮是一次污染物大量释放的过程，其污染物释放量远远大于静态条件下由孔隙水和上覆水体的浓度差而引起的释放量[25]。研究表明太湖底泥悬浮的临界风速在5～6.5m/s之间，一旦超过临界风速，底部的泥沙会发生悬浮，引起内源释放[26-27]。三维数值模拟可以更精确的模拟水土界面的物质交换过程及污染物到上覆水中的时空分布规律，Kopasakis K等[28]建立了浅水湖泊生态水动力模拟系统，包括水动力、沉积物和生态动力学模块；陈彧等[29](2010)运用生态动力学模型CAEDYM与三维水动力模型ELCOM的耦合模型模拟了太湖三维流场下的水质分布，并获得一套适用于太湖的水质模型参数；冯峰等[30]以浅水湖泊武汉东湖为研究对象，分析并比较了沉积物中碳、氮、磷等生源要素主要形态的垂向分布，结果表明pH值和氨氮都随深度增加而有所增加最终趋于稳定。陆健刚[31]等研究表明中水动力和强水动力下，水体中溶解态Cu、Zn、Pb含量从底层水体至表层水体呈对数增长。

底泥悬浮引起的泥孔隙水中的营养物质释放到上覆水中，水体营养物质浓度大大增加，为藻类提供营养元素，进而与湖泊富营养化产生密切联系。吴挺峰等[32]2008年9月在梅梁湾及贡湖湾水文水质及气象的观测结果表明，观测期间水动力强度对太湖北部湖区叶绿素a浓度垂向分层及蓝藻水华水平漂移均具有重要影响。王超等[33]建立了基于水动力空间分布的浅水湖泊模型，研究了大型浅水湖泊中微囊藻素的运输过程，结果表明，微囊藻素垂向分布与风力条件和对应的水动力条件有明显的响应关系，微风条件(<3m/s)下，蓝藻产生的区域从表层变为水下0.3m[34]。在水动力滞缓水域，蓝藻水华易在水表发生漂移堆积，在水动力强度较大水域，强烈的垂向混合作用能使蓝藻沿水深方向混合均匀，降低水华暴发风险。

风的持续时间及其动力性也是影响太湖蓝藻水华时空间变化的关键因子之一[35]。Wu等[36]利用卫星图像结合实地抽样调查研究发现，蓝藻水华面积与风速成负相关，较小风速有利于蓝藻水华形成，但短时间强风(30分钟平均风速超过6m/s)引起的强烈混合作用则会增加大面积蓝藻水华形成的机会，并且太湖地区4月和10月大风频繁可能正是蓝藻水华大面积发生的重要原因。李未等[37](2016)建立了太湖三维水动力—水质耦合模型，对短期内太湖湖泛易发水域及发生面积进行预测，准确率在80%以上。欧阳潇然等[38]利用FVCOM对太湖梅梁湾三维水温及水体中溶解氧的时空分布进行了模拟，相关研究标明，水下光场对于混合藻类初级生产力具有明显的影响，由于水体对于光照的衰减作用，太湖梅梁湾口出现光抑制的深度主要分布在水下0.1～0.3m区域[39]。

3 浅水湖泊数模关键参数研究

数值模拟过程中参数的率定至关重要，直接影响到模型的稳定性及计算结果的精度。目前对数值模拟过程中的参数研究很多，但在规范化和标准化上仍有欠缺，亟需形成一套完整的适用于浅水湖泊水动力模型的参数方案，本节内容对前人研究过程中的相关参数取值进行梳理和对比，以期为今后研究中的建模过程提供一定参考。

3.1 湖泊水动力模型关键参数

水动力过程主要涉及连续方程和动量方程，连续方程如下：

∂t(mζ)+∂x(myHu)+∂y(mxHv)+∂z(mw)=0

(1)

(2)

其中H是水深，u、v分别是流速分量；w是垂向速度。mx和my是拉梅系数。动量方程如下：

∂t(mHu)+∂x(myHuu)+∂y(mxHvu)+∂z(mwu)-(mf+v∂xmy-u∂ymx)Hv=-myH∂x(gζ+p)-my(∂xh-z∂xH)∂xp+∂z(mH-1Av∂zu)+Qu

(3)

∂t(mHv)+∂x(myHuv)+∂y(mxHvv)+∂z(mwv)+(mf+v∂xmy-u∂ymx)Hu=-mxH∂y(gζ+p)-mx(∂yh-z∂yH)∂zp+∂z(mH-1Av∂zv)+Qv

(4)

其中f是柯氏力系数，Av是垂直涡动粘滞系数，Qu和Qv是动量源汇项，g为重力加速度。

有关切应力项的处理如下，水面风切应力：

(5)

(6)

式中wx、wy为10m高空处风矢量在x、y 方向上的分量；Cf为风的拖曳系数，由下式计算：

(7)

其中，ρa与ρw分别是空气的密度与湖水的密度。

床面切应力：

(8)

(9)

式中，cb为湖底床面的拖曳系数，由下式计算：

(10)

其中，Δbl是无量纲底部层的厚度，σ0是无量纲的底部粗糙度。水动力数值模拟计算过程涉及多个参数输入，选取其中对数值模拟计算结果影响较大的5个参数进行分析。

(1)风应力拖曳系数

风是大型浅水湖泊水流运动最主要的驱动力之一，在风生流的形成中风应力起决定性作用。风应力拖曳系数决定了大气与湖泊间的动量传输率。在选择大型浅水湖泊水动力模型参数时，要着重率定风场参数，其中风应力拖曳系数对计算结果不确定性影响最高[40]。

在浅水湖泊水动力模拟中，风拖曳系数取值从1×10-6～3×10-6[41～43]各不相同，也有把拖曳系数作为随风速变化的变量，周婕等[44]指出，风应力拖曳系数为常数时，风应力和风速平方呈线性关系，风应力拖曳系数取为和风速有关的表达式时，两者为非线性关系，模拟结果显示，非线性关系的考虑相对更为合理。

(2)底部拖曳系数

底部拖曳系数Cb是空间变量，反映了水底糙度对能量损失的影响，文献中的底部拖曳系数通常采用公式C=h1/6/n计算[45]，h可近似取计算区域内的平均水深，n为曼宁系数(即底部糙率)，其取值直接影响底部拖曳系数(式10)，底部糙率是影响模型稳定的重要参数，浅水湖泊底部通常比较平坦，故文献中多采用经验常数[46]，也有根据不同水生植物分布及下垫面不同采取分区取值方案[47]。

(3)涡粘系数

涡粘系数分为水平涡粘系数和垂向涡粘系数，部分研究者将涡粘系数取为常数，另有部分研究者将垂向涡粘系数取为与风速有关的值，Koutitas[48]采用如下公式计算垂向涡粘系数

(11)

王惠中等[49]采用垂向变化的涡黏系数取值对太湖进行数值模拟，结果表明，水平流速的梯度变化在变涡黏系数时要小于常涡黏系数，底部切应力在常涡黏系数时要大于变涡黏系数，说明取常涡黏系数实际上夸大了湖水流层间的切变作用和作用于湖底泥层的作用力，认为在进行太湖三维风生流场模拟时不宜简单将垂向涡黏系数取常数值；Yongqi Wang[50]建立了Constance湖三维模型，采用Smagorinsky和2阶Mellor-Yamada模型计算水平和垂向涡粘系数，计算结果表明，水平和垂向湍流涡动粘滞系数对环流和温度的垂向分布起着重要的作用，而不能取为常数。

(4)扩散系数

扩散系数主要是对水质扩散而言，其对水动力模型计算结果影响较小，当河流/水库断面形态单一，岸坡倾角大，水深增加快，扩散系数呈现常数扩散系数特征，可简化取为常数[7, 17, 40, 51]，对滩槽形复式断面，深槽远离岸边，岸坡倾角小，滩地水深浅，水深增加缓慢，扩散系数呈现变扩散系数特征[52]。

3.2 参数经验值及经验公式

本文查阅了前人在浅水湖泊数值模拟研究中的参数取值及经验公式，现汇总供参考(见下表)。

表 浅水湖泊水动力模型参数取值表Tab. Values of hydrodynamic model parameters in shallow lakes

4 展 望

随着现场观测手段的发展和计算机水平的提高，浅水湖泊数值模拟研究受到的限制越来越小，富营养化是湖泊特别是浅水湖泊面临的主要环境问题，治理难度很大，其相关要素又与水深及水力条件空间分布有密切的响应关系，对浅水湖泊进行精细化的三维数值模拟变得越来越有必要。针对现有研究已取得的成果和存在的局限性，提出以下展望。

4.1 复杂边界条件下的生态动力学模型构建

由于湖泊现场观测难度大，大型浅水湖泊全湖区气象及水文观测数据极难获得，通常采用假设的定常风作为气象边界条件，未来可结合水利、气象等管理部门，开发数据共享平台，利用好现有的浅水湖泊现场观测数据，开发基于大气-水-底泥多边界全过程生态模型，研究水-气交界面动力条件及碳循环过程，探索变化边界条件下的气候与富营养化及沉积物再悬浮的循环和响应机制，建立更符合天然情况的全过程生态动力学数值模型。

4.2 多条件耦合的富营养化微观响应机制

富营养化产生需要在合适的光照、温度、营养盐浓度及微生物生长所需水力条件，水体中的氮磷营养盐浓度对微生物的降解功能有显著影响，营养盐分布也与微生物分布有密切联系，水中植物会影响水流结构分布，需要进一步开展上述要素之间的协同作用和迁移转化规律，研究多条件耦合下的富营养化产生微观机制。

4.3 数值模拟精细化及参数规范化研究

现有数学模型均基于一定假设条件下，且模拟过程中各项参数取值差异较大，尚未针对浅水湖泊形成一套规范化的参数取值方案，可结合前人研究，通过物理水槽试验，得出参数初始值或公式，并通过湖泊原型观测数据进行验证，最终得到适用于浅水湖泊的规范化参数取值方案。