基于EFDC模型的洱海水温模拟

张锦鹏，吴越1,，田泽斌1,，储昭升1,*，袁静1,，侯泽英1,

1.湖泊水污染治理与生态修复技术国家工程实验室 2.中国环境科学研究院湖泊环境研究所

水温是湖库研究中重要的环境参数[1]，也是水体蒸发、水质分析和水土界面物质交换研究不可缺少的研究资料[2]。水温的变化将直接或间接地影响湖泊的营养盐循环、浮游植物群落结构变化等过程[3-4]，垂向水温剖面影响水体的分层，是影响水体垂向混合的重要因素[5]。实地监测能如实反映湖泊水温的特点，但需要耗费大量的人力物力，监测结果存在时间间隔长、采样点有限等问题，且只能反映湖体某时间段的局部水温特征。为适应当前湖泊全面、动态评价和管理的需求，湖泊水动力的模拟研究应运而生。目前国内外许多湖泊已经开展了水温模拟方面的研究[6-14]，如程伟平等[10]对高原浅水湖泊——滇池的水温进行了模拟，发现对流传热相关系数对滇池表层水温影响较显著；甘衍军等[11]对深水水库——二滩水库的水温分层规律进行了分析，并比较了建库前后下游河道的差别；段扬等[12]对丹江口水库的水温进行了模拟，发现夏季库区的水温分层现象最为明显；赵一慧等[13]对不同水文年小浪底水库的水温分层现象进行了模拟；Yang等[14]对比利奇努克湖的水温分层和密度驱动环流进行了研究，并成功模拟了湖泊的表层、内部和底部3种水流。

洱海是云南省第二大高原湖泊，也是大理市的重要水源地，具有航运、供水、发电、旅游、调节气候等功能。近几十年来，洱海水体富营养化不断加重，已经从20世纪50年代的贫营养阶段逐渐进入到目前的富营养化初期阶段[15-16]。与其他亚热带深水湖泊不同，洱海没有明显的温跃层，因此，准确地模拟洱海水温，分析水温时空分布规律及可能产生的结果，可以为洱海水质模拟和管理评估提供科学依据。

1 EFDC模型

EFDC模型[17]最初由威廉玛丽大学弗吉尼亚海洋科学研究所的John Hamrick博士等开发，该模型涵盖了水动力、水质、沉积物和有毒物质等内容，在国际上获得了广泛的应用[18-19]。近年来，该模型在国内湖泊[20-22]、河流与河口[23]等地表水体的模拟、评价和管理等方面的应用也有增加的趋势。

EFDC模型在水平方向上使用笛卡儿坐标系或正交曲线坐标系，垂向上采用σ坐标系统，可在垂直方向上提供均匀的分辨率。σ坐标与直角坐标的转换公式为：

z=(z*+h)(ζ+h)=(z*+h)H

(1)

式中：z为σ坐标，无量纲；z*为相对于参考高度的垂向直角坐标，m；h和ζ分别为基于参考高度以下的水深和相对参考高度的水面高程，m；H为总水深，m。

EFDC模型中，σ坐标下水柱中温度和热传递的基本方程为[5]：

(2)

式中：x、y为水平方向的正交-曲线坐标，m；mx、my为坐标变换系数；T为温度，℃；w为垂向速度分量，m/s；I为太阳短波辐射强度，W/m2；ST为热交换的源汇项，J/s；P和Q分别为x、y方向的质量通量分量，m2/s。

基于立面二维水动力学和水质模型(CE-QUAL2-W2)[24]，EFDC模型采用温度平衡法计算水体表面的热交换，线性化后的表面热交换率为：

Haw=-Kaw(Ts-Te)

(3)

式中：Haw为表面热交换率，W/m2，Kaw为表面热交换系数，W/(m2·℃)；Ts为水面温度，℃；Te为平衡温度，℃。太阳短波辐射进入水体后被水柱吸收，从而使水柱水温上升。水柱中的光强度可用下式表示：

(4)

式中：Ke为消光系数，m-1。这里z*特指水面以下的深度。

包含沉积床温度的模拟可使深水湖库的水温模拟更加准确，EFDC模型中沉积床和水柱底层热交换率可表示为：

Hb=-(Kb,vU+Kb,c)(Tw-Tb)

(5)

(6)

式中：Hb为沉积床-水柱热交换率，W/m2；Kb,v为对流热交换系数，W/(m2·℃)；Kb,c为传导热交换系数，W/(m2·℃)；u1和v1分别为底层水流在x、y方向上的速度分量，m/s；Tw为底层水温，℃；Tb为沉积床温度，℃。紊流过程在水体的垂直混合中具有重要影响，由紊流扩散引起的垂向运输足以使水体完全混合[5]。

EFDC中的紊流模型将垂向紊动涡黏与扩散和紊动强度、紊动长度尺度、理查森数联系起来，紊动扩散的主要表达式为:

Av=φAA0ql

(7)

(8)

Ab=φKK0ql

(9)

(10)

(11)

式中：q为紊动强度，m2/s2；l为紊动长度尺度，m；Av为垂向紊动扩散系数；φA为稳定黏性系数；A0、R1、R2、R3、K0为常数；φK为稳定扩散系数；Rq为理查森数。

2 洱海三维水温模型构建

2.1 研究区域

洱海流域汇水面积为2 565 km2，湖泊平均水深为10.8 m，最大水深为20.8 m。当水位为1 965.8 m(八五高程，全文同)时，洱海水面面积为252 km2，库容约为27亿m3。洱海水源主要来自于北部的3条入湖河流(弥苴河、罗时江、永安江)，其次为西部的苍山十八溪，东部和南部分别有3条和2条主要入湖河流。洱海唯一的天然出湖口为洱海西南部的西洱河，另有人工开凿的引洱入宾出湖口。自1973年西洱河修筑大坝发电以来，洱海的水位转由人为控制。洱海地处高原亚热带季风区，干湿气候明显，且夏季暴雨多、冬季降水少，多年平均气温维持在15 ℃左右，洱海地区的风向以西南风为主。

2.2 地形构建与网格划分

采用1∶5 000数字化洱海水下地形图，利用邻域插值法进行插值，完成洱海水下地形构建。网格划分对于后续EFDC建模至关重要。目前常用的网格类型包括非结构化网格(三角形)、矩形网格和正交曲线网格。三角形网格和矩形网格都具有生成快速的特点，但三角形网格计算密集，需要较长的运行时间；矩形网格边界附近可能会出现虚假水流流动的现象，且不容易控制和修改网格密度。正交曲线网格通过正交化实现对不规则边界的适应性，同时由于网格和主流方向一致使得计算更容易收敛[24]，因此，采用正交曲线网格，利用Delft3D软件完成洱海水平方向上的网格划分。研究区域最终被划分为1 056个正交曲线网格，平均网格分辨率为380 m×604 m，垂向平均分为3层(图1)。

图1 洱海三维计算网格示意Fig.1 Diagram of the three-dimensional computing grids of Lake Erhai

2.3 条件设置

气象条件方面，太阳辐射和风场采用中国环境科学研究院大理洱海基地的实测数据，气压、蒸发、降雨和气温数据来自国家气象科学数据中心(http://data.cma.cn/)大理站，其中太阳辐射为6 h间隔的数据，其余为日测数据。对日测数据进行了线性插值以满足模型输入的需要。水文条件方面，采用洱海主要入湖河流(27条)的流量和水温作为入流边界，洱海南部的西洱河和引洱入宾2个出湖口作为出流边界(图2)。

各入湖河流的流量和水温均为每月1测的巡测数据，数据来源于云南省水文水资源局大理分局。模型采用热启动。首先采用洱海2016年1月实测的水温和水位数据作为初始条件冷启动运行1年，再将得到的结果作为热启动的初始条件输入到新模型中，运行时段为2017年1月1日—11月30日。将干边界设为0.05 m，湿边界设为0.10 m。

3 结果与讨论

3.1 参数率定与敏感性分析

注：括号内数字为国控和省控站点代号。图3 洱海水位和水温监测点分布Fig.3 Distributions of water level and water temperature determination stations of Lake Erhai

由于目前洱海水位由人为控制，水位空间变化不大，加之收集的数据受到限制，因此本研究仅采用西洱河和引洱入宾2个监测点(图3)的水位日测数据(来自大理市洱海保护管理局)进行率定。水温采用每月一测的11个监测点(图3)的湖体巡测数据(来自大理白族自治州环境监测站)，其垂向采样位置为水体表层0.5 m和底层0.5 m。湖体底部粗糙度是影响水位的重要参数，其经验值为0.01～0.10 m，经多次运行调试设为0.017 m。洱海水位模拟值与实测值对比如图4所示。由图4可知，洱海西洱河和引洱入宾2个水位监测点的水位模拟值与实测值接近，模拟结果的平均相对误差为0.003%，能较准确地反映洱海水位的时间变化趋势，可认为模拟结果是可信的。但2017年下半年水位模拟值和实测值出现了一定的偏差，模拟值低于实测值。水位的模拟偏差一般由地形、流量边界条件2类因素引起。洱海水位下降过程拟合得较好，说明地形模拟得较准确，可良好体现湖体库容变化；水位上升过程的偏差出现在雨季(7—8月)，这可能是由于洱海西侧受山地地形影响，主要入湖河流流程短，部分地表径流未汇入河道而通过地表漫流直接入湖，该部分流量缺乏实测数据造成的，这有待在后续模型模拟中进一步完善。

图4 2017年洱海水位模拟值与实测值对比Fig.4 Comparison between water level simulation and actual measurements of Lake Erhai in 2017

热通量参数能够调节湖体与外界的热量交换，决定湖体的整体热量水平；湍流参数则会影响热量在湖体内部的分布情况，二者都是水温模拟的重要调整项。通过文献调研和模型实际运行调整，将下列参数设为常数，即最大理查森数设为0.1，水平涡黏系数设为6，垂向分子扩散系数设为3.04×10-9，底床厚度设为2.7 m，传导热交换系数设为0.36，对流热交换系数设为0.03。拉丁超立方采样(LHS)方法是一种分层的蒙特卡洛抽样(MCS)方法，在避免MCS方法坍缩性的同时能够有效提升抽样的效率和代表性[24-25]。标准秩回归(SRR)方法通过求秩解决了输入与输出之间的非线性关系，是一种常用的敏感性分析方法[22]。通过文献调查[9-10,24]获得需要进行调整的主要水温参数(太阳短波辐射快速衰减系数、太阳短波辐射中的快速波占比、太阳短波辐射慢速衰减系数、背景垂向涡黏系数等)及其最大值和最小值(表1)，采用LHS[25-27]和SRR方法[24]，利用LHS随机抽样生成50组参数，运行50次模型并得到50组结果，以进行SRR敏感性分析。SRR的平方(SRRC2)表示每个参数的变化对整个输出结果方差的贡献率，其值越大则表明输出目标对该参数越敏感。以洱海全湖11个水温监测点模拟结果的平均纳什效率系数为目标，得到的敏感性分析结果如图5所示。

表1 模型输入参数的范围 Table 1 Range of model input parameters

图5 以全湖水温监测点平均水温纳什 效率系数为指标的敏感性分析结果Fig.5 Results of sensitivity analysis with Nash coefficient of average water temperature of the lake as index

由图5可知，FSWRATF的SRRC2最大，为18.7%，表明全湖水温对FSWRATF最敏感，不确定性贡献最大。甘衍军等[11]对二滩水库的模拟研究也表明，FSWRATF的增加会使水温温差增大。SWRATNF的SRRC2为14.7%，仅次于FSWRATF，说明全湖水温对SWRATNF较为敏感。穿过水面的太阳短波辐射能被水体吸收，这个过程可以加热一定深度的水体[5]。SWRATNS的SRRC2仅为0.4%，说明SWRATNS对全湖水温模拟结果的不确定性贡献较小，敏感性较差。最终确定的参数取值SWRATNF为0.055，FSWRATF为0.095，SWRATNS为0.055，AVO为1.55×10-6，Ws为0.99。

利用纳什效率系数和相对误差对模拟结果进行评估，洱海2017年11个水温监测点模拟误差统计见表2。由表2可知，大部分监测点的纳什效率系数在0.900左右，仅个别站点存在纳什效率系数有一定偏离的现象，平均相对误差为7.700%,总体上模拟结果较为准确，可反映洱海水温的时空变化情况。

表2 洱海11个监测点水温模拟误差统计

3.2 水温分布规律

图6 2017年洱海主要监测点水温模拟值与实测值对比Fig.6 Comparison of simulation and observation results at main stations in 2017 in Lake Erhai

图7 洱海不同季节模拟水温场Fig.7 Simulated water temperature field of Lake Erhai in different seasons

洱海夏末秋初(7—9月)为高水温期，正午全湖水温在20 ℃以上；冬春季(1—4月)为低水温期，水温常保持在10 ℃以上，无结冰现象。湖心、南湖心、北湖心3个监测点表层和底层水温的模拟值与实测值对比如图6所示。由图6可知，水温模拟结果较好地反映了洱海水温的实际变化趋势，1—2月为水温低谷期，随后水温缓慢上升，7—9月为高水温期，随后水温快速下降。

湖泊水温的空间分布特点往往与地形有关[28]。洱海不同季节水温场模拟结果如图7所示。由图7可知，洱海水温的空间分布主要为南北方向上的变化，且具有一定的季节特征。春夏季水温呈自北向南递减的特点，秋冬季呈现出南北低、中间高的特点。由于洱海中部为湖心盆地，是洱海的深水区，上述特点说明洱海水温的空间分布与湖盆地形和水深存在一定的关联性。此外，洱海水温的南北向分布特点可能与主要河流入流位置和出湖位置有关：北部的北三江是洱海主要水源补给渠道，入湖口分布在洱海北部，而主要出湖口——西洱河位于洱海的西南部，另一人工出湖口——引洱入宾出口位于洱海东南部，这导致洱海水体流动的总体趋势为自北向南，并对水温的空间分布造成一定影响。

水温分层对水体的垂向营养盐交换过程具有重要影响[8]。一般湖泊深达几米时，在夏秋季节上下层间就会出现明显的温差[5]。如二滩水库的表层、底层温差达14 ℃[11]，丹江口水库夏季的表层、底层温差达11 ℃[12]，抚仙湖表层、20 m水深处温差达8 ℃左右[29]。洱海湖心监测点(水深为20.8 m)多年实测水温剖面如图8所示。由图8可知，洱海水温的垂向温差并不明显，多在1.5 ℃以内，这与前人的研究结果基本一致[30-32]。

图8 洱海湖心监测点多年实测水温剖面Fig.8 Multi-years measured water temperature profile at lake center monitoring station of Lake Erhai

2017年洱海湖心监测点垂向水温模拟结果如图9所示。由图9结合图8可知，洱海表层、底层实测水温全年温差很小，绝大多数在1.5 ℃以内，水温垂向模拟结果很好地贴合了垂向表层、底层实测温差较小的特点。个别时段的模拟结果存在一定的误差，可能的原因包括可用的数据资料不全面、采样时的水温波动、采样误差等。

注：均为当日10:00水温数据。图9 2017年洱海湖心监测点垂向水温模拟结果Fig.9 Simulations of vertical water temperature at lake center monitoring station of Lake Erhai in 2017

3.3 水温变化的影响因素与结果

湖泊水温变化过程很大程度上由外部热通量、河流流入/流出和水动力过程决定，影响水温的因素包括气温、水深、季节变化、水平耗散、风场和潮汐产生的垂向混合、层化、入流带进的温度及人类活动的影响[5]，其中气温、太阳辐射和风场变化是影响湖泊水温的重要因素[11-12]。水温的季节变化主要受太阳辐射等热通量因素的影响[5]，夏秋季洱海地区太阳辐射强，湖体水温处于高位；冬春季太阳辐射弱，湖体水温偏低。水温的季节变化对于藻类的生长具有重要意义。许多研究[33-34]已经证明了蓝藻的高温生长优势，原因是蓝藻是原核生物，与绿藻、硅藻等其他真核藻类相比更偏好在高水温下生长。夏秋季洱海的高水温耦合暴雨径流带来的营养盐为藻华的发生创造了条件[31]。

风力扰动可能是洱海垂向水温弱分层现象的重要原因。洱海地区风大、风期长，年平均风速在2.3 m/s左右，1年中风速大于17 m/s的大风天数可超过100 d[16]，风对湖水的扰动强烈。研究表明，风场变化会加剧湖体垂向上的湍流运输，风场作用力导致的湖面波动及产生的水流剪切力足以使湖泊水体混合[12,35]。此外，洱海盆地西部和南部存在地热异常区[36-37]，可能促使垂向水温温差减小。而洱海营养盐的垂向变化并不明显[16]，这可能与洱海水温的垂向弱分层、湖区风力扰动强等相关。

4 结论

(1)根据洱海的地形、气象、水文等实测数据，基于EFDC模型建立了洱海三维水动力模型，并对2017年洱海水位、水温进行了模拟。采用拉丁超立方采样和标准秩回归方法对影响水温的重要参数进行了敏感性分析，发现湖体水温对太阳短波辐射中的快速波占比和太阳短波辐射快速衰减系数2个参数的变化较为敏感，而对太阳短波辐射慢速衰减系数变化的敏感度较弱。

(2)洱海全湖平均水温变化特点是春夏季水温升高，在8月达到峰值，秋冬季逐渐下降。洱海水温的空间分布随季节变化明显，春夏季呈自北向南递减的特点，秋冬季呈南北低、中间高的特点，这与洱海出入湖口的位置及水下地形有关。

(3)洱海水温主要受风力扰动、太阳辐射等因素的影响，水温的季节变化和垂向特征对洱海营养盐垂向变化和水华的发生具有重要意义。