基于Delft3D水动力模型对青草沙水库水温变化的模拟和预测

毕潇伟，徐 聪，何义亮

（上海交通大学环境科学与工程学院，上海 200240）

青草沙水库是上海市四大水源地之一，自2011年全面投入使用后日益突出的水质问题成为水库生态环境的重点。根据近年上海水务部门监测，库内水质已呈现营养化的趋势，张琴娟［1］对水库的水质监测发现藻类属于影响水库的重要污染物之一。所以对青草沙水库水质的研究显得尤为重要。

水库的蓄水调节过程使得库区呈现水位偏高、水流流速降低等现象，对水体的自净能力有相当大的影响［2］。水库的水动力特征，包含水位、流速等，对水库的水环境有着重要影响。流场的变化会影响水温的变化［3］，同时影响了水质点的迁移［4］。尤其是在短时间内，水体流场对水温的影响尤为明显。

通过流场的模型可以对水温进行精确模拟，水温的改变也会对水环境产生多方面的影响，比如水体中的溶解氧、水体中有机物以及水体中的各种微生物和浮游生物等［2］。浮游藻类是水环境质量的一个重要指标，在生物学水质评价中广泛运用［5］。水温对藻类的种类和数量都会产生较为明显的影响，尤其是夏季高温条件会引起藻类数量的大量增加［6］。可以通过藻类数量较为直接地反映青草沙水库的水质问题，所以本文以藻类作为典型问题来探讨温度对水质的潜在影响。因为水温和气温有着强烈的关系，大气和水面的热交换直接影响了水温的升高［7］，所以我们可以通过水动力模型来预测未来气温条件下青草沙水库的水温变化。

青草沙水库处于长江口的南北港分流口下方，由于同时受到了径流和潮流的作用，所以兼有了河、海的流场特征［8］。随着计算技术的广泛运用，数值模拟在水环境中的运用越来越成熟，在水动力数值模拟中最常用的模型有：EFDC模型、Mike系列软件、POM模型和Delft3D软件。本文选用的Delft3D采用ADI计算方法（隐式方向交替法），计算稳定性好，可以模拟复杂的边界条件。

国内外有很多关于全球未来气温的研究，张国庆等［9］基于模型预测了未来气温以每10年0.24℃的增长速度上升；IPCC第五次报告［10］预估21世纪末，全球气温升高至 2.6～4.8℃的可能性最大；郯俊岭等［11］对中国未来气温变化进行预估，提出在21世纪末上升温度很可能超过2℃。高学杰等［12］以温室效应为机理，以长江中下游地区气候的变化为研究对象，探究在温室效应加重的情况下气温的变化，发现长江中下游区域气温升高值比内陆高2.3℃ 以上。以国外IPCC的报告对未来全球气候变化的预测为权威，并基于国内外的研究结果，考虑到青草沙处在长江口这个特殊位置，受到长江径流和海水的影响，在全球气候变化下，长江口的气温也会略微高于平均。

Takemura等［13］进行了水温与藻类生长关系的试验，发现温度高于11℃、并低于一定温度下，光合作用效率随温度的变化而呈现比例关系。多项研究都表明大多数藻类最适宜生长温度为30℃［14］。徐良等［15］在Arrhenius规律基础上对Eppley探究水温对藻类生长规律的影响做了进一步地探究，修正了新的公式用作藻类增殖系数。在藻类生长过程中，当温度超过最适温度时，藻类生长符合公式Fp（T）＝ pt－T；当温度低于最适温度时，藻类生长符合公式Fp（T）＝pT－t（T 是最适温度，t为实际温度［16］）。当水温超过30℃，水温的升高反而会抑制藻类的生长［17］。

本研究利用Delf3D-Flow模块，在夏季对青草沙水库运行情况进行水动力模拟。通过流速率定和验证得到良好的水动力模型，在验证好的模型基础上对温度进行模拟。为了更好地探究水温变化，在原模型基础上改变空气温度用以模拟未来气温下的水温变化，从而更好地探究青草沙水库的水温变化和潜在的水质问题。

1 青草沙水库流场的模拟和验证

1.1 模型的建立

1.1.1 Flow 模块的基本原理

在模块的流体基本计算中，将水体假设为不可压缩流体，并忽略温度的改变对流场的基本影响。Delft3D-flow模块建立在Navier-Stokes方程基础上，符合浅水特性和Boussinesq假定。通过ADI计算方法对设定坐标系下的控制方程组进行离散求解，在忽略垂向加速度影响的前提下，推导出静水压强假定下的水流方程［18］。

1.1.2 边界条件和参数设定

由于青草沙水库位于长江口潮流区域内，受到径流、潮流的影响，水库附近长江水位呈现波动［19］。水库一般运行规律为：在水位高的时候取水，水位越高，可取流量越大。青草沙水库的引排水方式由泵、闸相结合构成，分别在咸潮初期、咸潮期、咸潮末期进行提水、抢水和补水。上下游水闸共同维护水库水体的交换和库内水体的流动，只有在非咸潮期通过下游水闸排水。因此，具体流量及取水时间主要根据附近水位具有的潮汐变化规律进行设定。通过2011年的流场模拟与实测分析，夏季入水流量约为 700～900 m3／s，出水口流量约为 100～200 m3／s，进出水流量保持总量守恒。

初始条件指定了水体的初始状态，本文模拟采用的初始条件是“冷启动”（流速为零和自由水面为静止的状态）。使用库内的平均水位作为水库初始时刻水位，将模型计算分为10层，按照随着水深不同，每层深度随之变化，但是所占百分比不同原理设置。基于涡黏性系数和涡扩散系数、速度尺以及长度尺度成正比的假设，采用k-ε紊流模式来确定这些系数。黏度系数选取经验值，因为涡黏性的各向异性，选取的水平黏度系数远大于垂直黏度系数［18，20-21］。基于前人的经验参考和多次模型的尝试率定，最后选定水平黏度系数为 2 m2／s，垂直黏度系数为 0.03 m2／s，水平扩散系数为 2 m2／s，垂直扩散系数为 0.001 m2／s，曼宁系数为0.02。在物理参数的考虑中，重力加速度选取 9.81 m／s2，水体密度选取 1 000 kg ／m3［22］。

1.1.3 区域网格划分

青草沙水库地处长江的江心口，由两个库区构成，青草沙库区以及中央沙库区。中央沙库区作为库头，处在长江口南北港分汊口的中央水区域，但是由于滩面较高，并不记入模拟计算区域。青草沙库区属于北港，和长兴岛之间有一道潮沟，所以主要的计算区域选取的是青草沙库区。由于在青草沙中存在垦区，垦区的特性近似于陆地，所以垦区本身不记入计算区域，并在垦区和水区域之间进行了边界处理。计算区域采用正交模式，最后划分的网格如图1所示。利用Regrid模块建立正四边形网络，在垂向作分层处理，层数为10。在网格绘制中，对靠近边角的网格进行光滑处理，对网格适当加密，并去除不合适网格。为了提高计算精确度，网格横纵向的夹角余弦值都控制在0～0.02。计算区域网格的划分和实际情况拟合良好，能够较好地满足流场模型的计算要求和精确度运算。

图1 青草沙水库网格划分图Fig.1 Mesh Graph of Qingcaosha Reservoir

1.2 流场模拟和验证分析

Delft-Flow模块的主要功能是模拟不同条件下的水流，通过将地形等实际情况简化为网格等计算条件，在计算中加入不同的边界设定和初始条件，来模拟风速等外界因素对水体表面以及内部的作用，这也是本文的核心内容［22］。由于青草沙水库受到海水的潮汐影响［23］，而潮汐的循环是从小潮到中潮，再到大潮，最后恢复到小潮。所以大潮时水库外的水位较高，引入的水流量在这个时候比较大，库内的流场也因此受到影响。对水库实测数据进行分析发现，小潮、中潮、大潮时库内流场的变动尤为明显，所以在模拟时考虑了潮汐的问题。

运行时间选取2011年7月2日～2011年9月30日，对这个时间段的流场进行计算。验证时间选取了2011年9月23日～2011年9月30日，采用实测数据进行对比分析。考虑到实测位置的比较，选取了具有代表性的R1、R2、R3、R4点。R1点是在库首，靠近上游水闸的位置，R2、R3点是在经过垦区分汊后的两侧位置，R4则是在靠近下游水闸的位置。分别选取了表层和底层的模拟和实测数据进行对比，但是因为青草沙水库本身属于浅水型水库，所以分层现象也不是特别明显。R1点因为靠近上游水闸，所处位置为浅水，大部分时间段处于低流速状态，分层现象并不明显，所以没有使用表层和底层的变化来分析。利用实测数据对模拟数据进行拟合对比计算，依据式（1）［24］计算拟合系数 NSE。

图2 2011年夏季R1流速模拟实测对比图Fig.2 Simulation and Actual Measurement of R1 Flow Rate in 2011 Summer

图3 2011年夏季R2～R4底层、表层流速实测和模拟对比图Fig.3 Comparison of Simulation and Actual Measurement of Bottom and Surface Layers’Velocity of R2～R4 in 2011 Summer

2 青草沙水动力模型的温度模拟

通过流场的验证，得到适用于青草沙水库的水动力模型。基于该模型，采用2011年的实测气温对青草沙水库进行水温模拟。通过2011年的模拟结果对水温的分布进行分析。

2.1 Ocean模式的添加

此次温度模拟添加了Ocean模式。该模型不仅考虑了温度的横向传递和在水—空气界面的垂向传递，还考虑了太阳辐射在水面的接受率。有效的太阳辐射依据太阳的垂直辐射强度、大气湿度和覆云率，从而更加合理地模拟出水库整体的温度变化。考虑自由曲面的热交换太阳能（短波）、大气（长波）辐射的单独效应，以及由于回辐射、蒸发、对流产生的热损失，计算如式（2）。

其中：Qsn—净太阳辐射；

Qan—净入射大气辐射；

Qbr—反向辐射；

Qev—蒸发热通量；

Qco—对流热通量。

本试验通过添加Ocean模式，模拟了青草沙水库库区上方的太阳辐射条件，水温的模拟因此更具准确性。本试验采用的Ocean模式中的气象数据全部来源于中国气象数据网。

2.2 实际气温下水温模拟

2.2.1 水温验证

选取R3和R4的表层和底层温度进行验证和分析。如图 4所示，水温的模拟值维持在 22～31℃，只有在高温时期的短时间内水温温度会有超过30℃。通过对比实测水温值和模拟水温值可以看出拟合效果良好，该水温模型符合实际情况。

图4 2011年夏季R3、R4表层和底层温度实测和模拟对比图Fig.4 Comparison Between Simulation and Actual Measurement of Surface and Bottom Layers’Temperature of R3，R4 in 2011 Summer

2.2.2 水温水平分布

从模拟结果可以直观了解水温分布情况。如图5所示，水温变化较大的时间阶段是7月初以及9月末。水体从上游开闸口进入后，温度随着水体向下游扩散而慢慢变化。上游引入清水，水在进入水库后在库区内慢慢形成扩散状态。随着水体扩散，加之水流流速的影响，从入水口流入后水温在上游就有了较快地变化，逐步扩散到垦区附近。由于空气的温度、太阳辐射等综合因素的影响，水温逐步升高。7月初的前10 d，整体水域的水温就已经全部发生了改变。在20～30 d，库区内的水温逐步趋向于升高。8月份水库的水温达到高峰值，并在8月份维持这样的高温值。进入9月后，又呈现出跟7月份相似的变化规律，随着气温的降低，水温也逐渐降低，从入水口开始温度慢慢降低，然后扩散到全部计算区域。

图5 2011年夏季模拟水温分布变化Fig.5 Simulation of Water Temperature Distribution in 2011 Summer

2.2.3 水温垂向分布

从各点的垂向分布看，R1处在上游水闸附近，靠近入水口，水流量短时间内较大，水深较浅，水温模拟在该模型模拟期间在垂向上几乎没有分层现象出现；R2处在垦区南北两侧，水深相对都较浅，分层现象在短时间内也不明显；R4较另外三个点水深较深，有一定的分层，但是最高温度差不超过1℃。如图6所示，R1几乎没有分层现象，R2、R3的分层现象不明显，R4有一定分层现象。

图6 2011年夏季模拟水温表层底层对比Fig.6 Comparison of Simulated Water Temperature Between Surface and Bottom Layers’in 2011 Summer

3 未来不同气温条件下青草沙水库水温变化及潜在的水质影响

基于国内外文献资料，未来全球气温在逐步升高。为了探究气温变化对水库水温的影响，本文通过改变模型的气温条件来预测水温的变化。从Nasa Observatory［24］提供的全球海洋表面温度变化来看，青草沙所处的长江口附近海域温度从7月～8月是一个升温的过程。在夏季的7月、8月、9月期间，可以看出细微变化，如图7所示。长江流域的温度一般来说略高于全球海洋平均温度，但是总体变化趋势是一致的。这个细微的变化和本文所采用的水动力模拟下的水温变化趋势是一致的，也进一步验证可以在全球气候变化下探讨该水动力模型。

3.1 2011年夏季水温变化下藻类问题

图7 2011年夏季全球海洋温度变化Fig.7 Gobal Ocean Temperature Change in 2011 Summer

不同温度条件下藻类的生长大不相同，而夏季温度的升高对藻类的生长提供了有利条件。以R4点为例，对R4附近实测藻类的数量变化曲线（图8）以及R4点水温模拟变化曲线对比分析发现（图9），在模拟的三个月中，藻类个数和温度的变化规律是相似的，在一定程度上说明，随着温度的变化，藻类的个数相应有一定的变化，并且在一定范围内，温度升高，藻类也相应有所增加。

图8 2011年R4附近的藻类个数变化曲线Fig.8 The Curve of Algae Change in R4 in 2011

图9 2011年夏季R4表层温度模拟和实测的变化Fig.9 Simulation and Measurement of R4 Surface Layers’Temperature in 2011 Summer

3.2 21世纪末的水温预测

本文假设21世纪末气温升高4.5℃，以此为21世纪末水温预测的基础。对21世纪末青草沙水库水温进行预测分析。在模拟结果中，选择R2点的表层水温变化进行分析。模拟结果如图10所示，从21世纪末气温上升的情况下的模拟来看，水温从7月中旬～9月下旬都处于30℃以上。从初始温度开始，整个夏季处于一直升高的状况，只在夏季末有些许下降趋势。而当夏季温度稳定后，通过对比2011年气温下的水温情况，21世纪夏季模拟水温总体偏高。对比最高温度发现，假设气温下水温的最高值相比2011年气温下水温的最高值增加了约3℃，对比R2点两次模拟的水温曲线发现，未来气温条件下的水温在进入稳定期后，要经历比2011年气温条件下更长的高水温时间段，而且进入高温时间段的时间有所提前。

图10 21世纪末气温条件和2011年气温条件下模拟R2表层水温对比图Fig.10 Comparison of R2 Surface Layers’Water Temperature Under the Temperature Condition at the End of 21st Century and in 2011

3.3 探究21世纪内气温影响下藻类的变化可能性

为了探究随全球气候变化导致的青草沙水库水温变化的具体情况，在本文模拟方案中，基于参考文献假定未来气温升高趋势前提下，分别对气温升高1、2.5、4.5℃进行水温模拟。基于藻类的适宜温度和增殖系数公式作为计算基础，来推断藻类的可能变化。

图11 21世纪末气温条件下模拟R2、R3、R4点表层水温变化情况Fig.11 Simulation of R2，R3 and R4 Surface Layers’Water Temperature at Temperature Condition of 21stcentury

模拟结果如图11所示，当气温升高 1℃时，水温在7月底达到30℃，在8月上旬都维持在高温时间段，在高温时间段的水温波动不大，最高温度不超过31℃。当气温升高2.5℃时，进入30℃的时间点提前了，高温时间段变长，一直维持到9月中旬。当气温升高4.5℃时，对比前2个气温模拟下的高温值和高温时间段又有所加长，最高水温已经达到34℃。通过对比3个气温条件下的结果，发现随着气温升高，水温的高温时间段拉长，水温上升会比之前的更早。对比3次模拟结果，发现温度升高后导致的直接结果就是水温高温时间段在夏季变长，而且水温上升时间点提前。对于藻类的影响在于，在原本适宜藻类生长的夏季温度过高，超过藻类生长最适宜温度，适宜藻类生长的时间段提前，藻类的暴发也因此可能提前。

4 结论和建议

本文采用Delft3D-Flow模块对青草沙水库流场进行了验证和分析，模拟了夏季青草沙流场动态，得到拟合良好的水动力模型。基于该水动力模型，在添加Ocean模式的基础上分别模拟不同气温条件下青草沙水库水温情况，并以藻类问题作为水质的典型案例进行分析和预测。

通过青草沙水库夏季流场模型的构建，模拟了水库库区水流流态和流速。模拟流场和实测数据的对比拟合在误差范围内，证实该水动力模型可用于青草沙的水动力模拟。通过模拟夏季期间2011年实测气温下的水体温度情况和不同气温假设条件下的水体温度，主要得到以下结论。

（1）流速在水库上中下游几乎没有分层现象。温度分层上，R1靠近入水口，流速大，水深浅，几乎没有分层现象；R2、R3水深也较浅，只有在白天时间段略有影响；R4较其他点水深，有分层现象。

（2）水温水平分布呈现一定的规律，夏季初水温从上游水闸口开始渐变，慢慢扩散到库区内部水体；夏季中期库区内全部水温达到最高值；夏季末，随着气温的降低，水温又从上游水闸口开始慢慢降低。

（3）对比夏季中2011年和21世纪末气温条件下的青草沙垦区附近点水温变化发现，21世纪末最高水温值要超过 2011年最高水温值，差值接近3℃。

（4）对比分析不同气温条件模拟下青草沙水库水温的变化发现，随着气温的升高，水温高温时间段变长，升温时间点提前。藻类生长高峰期可能不再是夏季，藻类暴发时间点会提前。

［1］张琴娟.青草沙水库的长江原水水质调查［J］.中国科技投资，2016（22）：375，388.

［2］王辉.大伙房水库流场及水温分布的数值模拟研究［D］.大连：大连理工大学，2015.

［3］陈泾.青草沙水库流场、滞留时间及盐水入侵来源［D］.上海：华东师范大学，2014.

［4］李杰，袁建忠.青草沙水库流场和水质点数学模型计算分析［C］.苏州：中国水利学会2007学术年会人类活动与河口分会场，2007.

［5］李宝林，王玉亭，张路增.以浮游植物评价达赉湖水质污染及营养水平［J］.水生生物学报，1993，39（1）：27-34.

［6］黄晓琛.藻类生长过程中的生态特性及水华预测方法研究［D］.上海：上海交通大学，2009.

［7］杨东方，吴建平，曲延峰，等.地球生态系统的气温和水温补充机制［J］.海洋科学进展，2007，25（1）：121-126.

［8］李茂学，刘小梅，付新永，等.青草沙水库取水泵闸规模论证［J］.给水排水，2009，35（2）：50-54.

［9］张国庆，李鹏飞，黄立，等.全球未来50年平均气温的时间序列分析与预测［J］.甘肃科技，2008，24（17）：72-74.

［10］Stocker T F，Qin D，Plattner G K，et al.IPCC，2013：Climate change 2013：The physical science basis.Contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change［J］.Computational Geometry，2013，18（2）：95-123.

［11］郯俊岭，江志红，马婷婷.基于贝叶斯模型的中国未来气温变化预估及不确定性分析［J］.气象学报，2016，74（4）：583-597.

［12］高学杰，林一骅，赵宗慈，等.温室效应对我国长江中下游地区气候的影响——数值模拟［J］.自然灾害学报，2004，13（1）：38-43.

［13］Takamura N，Iwakuma T，Yasuno M.Photosynthesis and primary production of Microcystis aeruginosa Kütz.in Lake Kasumigaura［J］.Journal of Plankton Research，1985，7（3）：303-312.

［14］逯南南，褚福敏.温度及光照强度对藻类生长的影响［C］.山东：2010“黄河杯”城镇饮用水安全保障技术论坛暨城市供水水质监测技术交流会，2010.

［15］徐良，冯平，孙冬梅，等.水温对藻类生长变化影响的数值模拟［J］.安全与环境学报，2013，13（5）：76-81.

［16］柴小颖.光照和温度对三峡库区典型水华藻类生长的影响研究［D］.重庆：重庆大学，2009.

［17］廖振华.青草沙水库水动力数值模拟［D］.上海：上海交通大学，2011.

［18］乔飞，郑丙辉，雷坤，等.长江下游及河口区水动力特征［J］.环境科学研究，2017（3）：389-397.

［19］Maa P Y.An efficient horizontal two-dimensional hydrodynamic model［J］.Coastal Engineering，1990，14（1）：1-18.

［20］许旭峰，刘青泉.太湖风生流特征的数值模拟研究［J］.水动力学研究与进展，2009，24（4）：512-518.

［21］张仁徽，陈宇里，耿钊.网格尺寸对Delft3D有限元水流场仿真精度影响的分析［J］.应用科技，2015，42（1）：57-61.

［22］刘小梅，王晓鹏，关许为，等.青草沙避咸蓄淡水库库容与特征水位研究［J］.水利水电技术，2009，40（7）：5-8.

［23］Moriasi D N，Arnold J G，Van Liew M W，et al.Model evaluation guidelines for systematic quantification of accuracy in watershed simulations［J］.Transactions of the Asabe，2007，50（3）：885-900.

［24］NASA Observatory ［DB／OL ］.http： ／／www.earthobservatory.nasa.gov／.