七里海潟湖近百年水体交换能力变化特征

摘要

基于三重嵌套的非结构网格建立七里海潟湖水动力和物质输运模型，采用实测资料对模型的潮位、流速和流向进行验证，结合卫星遥感影像资料，模拟1950—2018年七里海潟湖地貌演变过程，利用水体冲刷时间尺度探究不同工况下潟湖水体交换能力的变化特征。结果表明，自然条件下，径流改变较海平面变化对水体交换能力的影响更为显著；挡潮闸的开闭是影响潟湖水体交换能力的关键因素，同时潟湖围垦降低了纳潮量也减弱了水体交换能力，但降低幅度较小；潮汐通道工程的变化是影响水体交换能力的主要原因，潮汐通道由弯曲变为顺直及其拓宽均增强了潟湖的水体交换能力。

关键词

潟湖; 水体交换能力; 挡潮闸; 数学模型

潟湖作为陆地水域和海洋水域连接的特殊地貌形态，受到自然和人类活动的共同影响，拥有丰富的物种多样性和自然资源，其生态系统宝贵且敏感^［1-2］。随着沿海经济的快速发展，养殖废水，渔港排污和入湖河流排污使潟湖水污染问题日趋严重，藻类等富营养化问题频发^［3-5］。七里海潟湖是国内仅存的现代潟湖之一，也是昌黎黄金海岸国家级自然保护区的主要保护对象之一。然而从20世纪70年代开始，由于人类高强度地进行大规模的围垦和养殖活动，使得七里海潟湖淤积严重，水域面积严重萎缩，自然生态系统遭到严重破坏。

探明潟湖水动力和水体交换的变化机制，是保证潟湖生态修复和生态保护系统可持续发展的基础^［6］。Dias等^［7］基于拉格朗日法研究了葡萄牙Riade Aveiro潟湖的水体交换能力，得出潟湖整体的水体交换能力较强，但潟湖内部河道处的水体交换能力较弱。Jeyar等^［8］在摩洛哥Nador潟湖的水体交换研究中采用对流扩散模型，结果表明大潮较小潮条件下水体交换的过程更快。Cucco等 ^［9］和Umgiesser等^［10］分别针对Venice潟湖和Curonian潟湖采用欧拉法进行水体交换能力的研究，水体交换能力季节性特征显著，风和径流是水体交换能力主要的影响因素。Kuang等^［11］基于欧拉法研究了七里海潟湖水体交换能力的季节性变化特征及其影响因素，认为径流和潮位与水体交换能力的影响呈明显的二阶指数函数分布。匡翠萍等^［6］探究了清淤疏浚工程对七里海潟湖的影响，表明清淤疏浚工程在提升潟湖纳潮量的同时有助于改善潟湖的水体交换能力。本文基于水动力和物质输运数学模型，剖析长时间尺度下海平面上升、径流变化和人类活动作用对七里海潟湖不同时期水体交换能力变化的影响程度，为七里海潟湖湿地生态保护修复技术体系的建设提供参考依据，同时也可以为世界范围内不同区域的浅水潟湖的水体交换能力研究和生态修复技术体系的完善提供一定的理论基础。

1 研究方法

本研究采用欧拉法进行水体冲刷时间的计算。通过分析研究区域内整体的初始水体质量的变化过程，得到整体的剩余质量的指数衰减函数，将整体的初始水体质量衰减至初始质量1/e（约为37%）所需的时间确定为水体冲刷时间。水体冲刷时间越短，表明区域的水体交换能力越强。故基于变化的地形，综合考虑气候变化引起的海平面上升和径流的变化，同时考虑工程的影响，利用MIKE21软件建立七里海潟湖二维水动力和物质输运模型数学模型^［12］，利用有限体积法求解二维不可压缩Reynolds平均Navier-Stokes方程，考虑了Boussinesq近似、浅水假定和静水压力的假定，并耦合保守物质的物质守恒方程建立物质输运模型，用于评价潟湖的水体交换能力。

数学模型的具体控制方程如下：

连续性方程为

$$\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{\partial hv}{\partial y}=hS$$

（1）

式中：t为时间；x和y为笛卡尔坐标；u和v分别为x和y方向垂向平均流速分量；h为总水深，$h=d+\eta$，d为静止水深，$\eta$为水位；$S$为源流量。

物质输运模型方程为

$$\frac{\partial h{c}_m}{\partial t}+\frac{\partial hu{c}_m}{\partial x}+\frac{\partial hv{c}_m}{\partial y}=$$

$$\begin{array}{l}\frac{\partial }{\partial x}\left(D_{cx}\frac{\partial h{c}_m}{\partial x}\right)+\frac{\partial }{\partial y}\left(D_{cy}\frac{\partial h{c}_m}{\partial y}\right)+\\ h{c}_{s}S-h{k}_{p}c{}_m{}^{}\end{array}$$

（2）

式中：$c_m$为物质浓度；$D_{cx}$、$D_{cy}$为物质的扩散系数；$c_s$为源的物质浓度；$k_p$为衰减系数。

1.1　计算区域与模型设置

七里海潟湖模型平面坐标系选用北京54高斯克鲁格坐标系，采用三重非结构嵌套网格（图1）。大模型为渤海模型，中模型模拟区域为秦皇岛海域，七里海潟湖模型为加密后的小模型。网格包含6 656个节点和11 818个单元，空间步长为5~50 m。模型中变化的地形来源于七里海潟湖地貌演变的模拟结果，该模拟基于海平面和径流数据及卫星遥感影像，耦合水动力模型，综合考虑多年来海平面上升过程、径流变化及人类活动对地貌的影响。潟湖岸边以岸线即固壁边界进行考虑，外海开边界位于新开口潮汐通道口门处，采用秦皇岛海域中模型计算得到的Flather边界进行控制，4条入湖河流流量为实测月均流量。

图1  研究区域、网格及验证点位置

Fig.1  Study area, computational grid， and positions of validated points

七里海潟湖由4条河流流入，流量变化采用文献［13］及全国水文统计年报中搜集的滦河下游滦县站从1950—2018年河流径流数据河流流量进行换算得到。风应力采用欧洲中期天气预报中心（ECMWF）给出的海面10 m以上、间隔为6 h的风场资料。底部摩擦力曼宁糙率系数由水深和底部中值粒径共同控制，取值范围为0.011~0.015，时间步长由模型自动调节为0.001~60 s，CFL（courant-friedrichs-lewy）数限值为0.8。水平涡黏系数采用Smagorinsky公式计算。模型采用干湿动边界处理技术，模型中干点临界水深取0.005 m，湿点临界水深取0.05 m。

1.2　数学模型验证

模型验证采用2016年9月23日—2016年10月23日七里海潟湖潮汐通道内的潮位验证点WL的实测数据以及2016年9月23日11：00—2016年9月24日11：00的潮流验证点C1~C3（测点位置见图1）的实测数据对七里海潟湖水动力模型进行验证，验证结果如图2所示，其中潮流验证点C1位于潟湖内， C2位于潮汐通道内， C3位于秦皇岛海域。

图2  潮位及潮流验证

Fig.2  Verification of tidal level and tidal current

采用Wilmott^［14］提出的效率评价系数对计算结果进行定量评价，Skill数的具体计算方程如下：

$$S=\mathrm{1}-\frac{{\displaystyle \sum _{i=\mathrm{1}}^N}{\left|M-D\right|}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^N}{\left(\left|M-\overline{D}\right|+\left|D-\overline{D}\right|\right)}^{\mathrm{2}}}$$

（3）

式中：S为Skill数；$M$为模型计算值；$D$为实测值；$\overline{D}$为实测平均值；$N$为样本数量。当Skill值为1时，代表模型计算值和实测值之间完全一致；Skill值大于0.65时，表示模型计算结果为“极好”；Skill值在0.65至0.50之间时，表示模型计算结果为“非常好”；Skill值在0.50至0.20之间时，表示模型计算结果为“好”；Skill值小于0.20时，表示模型计算结果为“差”。各测站Skill值如表1所示，可见各测站结果均为极好，表示模型能够较好地模拟实际情况。

表1  模型评价

Tab.1  Evaluation of model

类型	点位	Skill效率系数	模型评价
潮位	WL	0.90	极好
流速	C1	0.68	极好
	C2	0.76	极好
	C3	0.75	极好
流向	C1	0.97	极好
	C2	0.98	极好
	C3	0.96	极好

1.3　水体交换能力研究

采用欧拉法计算水体冲刷时间进行水体交换能力的研究。具体方法是利用数学模型中保守物质的输运模型，求解对流扩散方程，将研究的水体染成单位浓度的示踪剂，模拟区域内水体物质浓度（$c_m$）随时间的变化过程。通过$M=V{c}_m$可得到水体质量（$M$）随时间的变化过程，其中$V$为研究区域水体的体积，剩余初始水体质量随时间的变化过程可以通过最小二乘法拟合，得到剩余质量的指数衰减函数为

$$M_t/M_{\mathrm{0}}=a{\mathrm{e}}^{-kt}+b$$

（4）

式中：$M_t$和$M_{\mathrm{0}}$分别表示水体在$t$时刻的质量和初始时刻的质量；$\mathrm{e}$为自然对数的底；$a\mathrm{、}b\mathrm{、}k$分别为通过最小二乘法拟合得到的系数。

2 1950—2018年七里海潟湖水体交换能力变化特征

1950—2018年间，由于海平面、径流以及人类活动的影响，七里海潟湖地貌发生变化，对水体交换的影响显著。20世纪50年代至60年代，七里海潟湖区域主要属自然演变。20世纪70年代，七里海潟湖潮汐通道维持自然形态，对入湖河流进行了综合治理，增强了排洪和泄洪能力，修建了25 km的围堤，此时潟湖水域面积为8 km²。同期在潟湖潮汐通道内修建挡潮闸，提高七里海调蓄能力，闸底板高程-1.5 m，闸顶高程3.2 m^［15］。由于挡潮闸的修建，减弱了潟湖的水动力及水体交换能力，泥沙在潮汐通道内发生大量淤积，致使潮汐通道向南迁移。1979年，挡潮闸被停止使用，闸门不再闭合，始终保持过水的状态^［16］。从1984年的卫星图片（图3b）可以明显看出，潟湖大部分岸线修建了围堤，水域面积缩小至7 km²，挡潮闸的存在使潮汐通道口门宽度缩窄至85 m，潮汐通道的形态发生改变^［17］。潟湖排水能力的降低和沿岸流的作用，使潮汐通道向南侧移动。

图3  七里海潟湖卫星图片

Fig. 3  Satellite images of Qilihai Lagoon

1984—1991年（图3c），七里海潟湖围垦区域主要用于水产养殖，水域面积在1990年缩窄至5.2 km²。1991年，疏浚工程使潮汐通道由弯曲变为顺直，潮汐通道外修建有防波堤。1991—2000年（图3d），由于养殖围垦的进一步加剧，潟湖水域面积缩小至2.6 km²，2001—2008年（图3e），围垦工程显著减少，水域面积仅减少了0.3 km²，缩小至2.3 km²。2008—2018年（图3f），挡潮闸上半部分拆除，保留了水下结构，并在原始的位置建桥连接潮汐通道两岸，桥梁采用墩式结构，桥面高4.7 m。同期，潮汐通道东部靠近秦皇岛海域部分拓宽，平均宽度由110 m增加至150 m，外海防波堤区域建立突堤防止海滩侵蚀。

根据卫星图片及文献数据，为研究长时间尺度下七里海潟湖水体交换能力变化特征，将该时期按照工况建设情况分为5个阶段进行模拟，分析不同工况条件下潟湖水体冲刷时间的变化。不同阶段的工况条件及模拟时间见表2。

表2  地貌演变模拟时期、参考卫星图片及主要工况条件

Tab. 2  Different periods of morphological simulations,reference satellite images, and construction conditions

模拟阶段	模拟起止日期/（年‒月‒日）	卫星图片年份/年	主要工况条件
第一阶段	1950‒01‒01—1970‒01‒01	1970	自然演变
第二阶段	1970‒01‒01—1980‒01‒01	1984	挡潮闸的建设
第三阶段	1980‒01‒01—1990‒01‒01	1990	潟湖围垦、挡潮闸始终开启
第四阶段	1990‒01‒01—2008‒01‒01	2008	围垦、潮汐通道顺直
第五阶段	2008‒01‒01—2019‒01‒01	2018	建桥、潮汐通道近海侧拓宽

2.1　第一阶段水体交换能力变化特征

图4为1950—2018年七里海潟湖水体冲刷时间及所对应的海平面、径流和主要人类活动。第一阶段（1950—1969年）七里海潟湖主要为自然演变的过程，潮汐通道入湖口门宽度约300 m，对潮流的阻碍作用较弱，平均水体冲刷时间为79 h。据图4显示，该阶段水体交换能力最强的时间为1960年，对应水体冲刷时间为46 h，1962年水体交换能力最弱，对应水体冲刷时间为122 h。但1960年的年均径流最大，为13.75 m³·s^-1，1962年的年均径流最小，为3.39 m³·s^-1，故径流量的增大会提高潟湖的水体交换能力。同时，在比较1964年和1965年的水体冲刷能力时发现，二者年均径流均为5.77 m³·s^-1，但1964年的海平面-0.12 m，较1965年的-0.20 m更高，但水体冲刷时间为82 h较1965年的86 h更小，说明1964年的水体交换能力更强。因此，在径流条件相同时，海平面升高会促进潟湖水体交换。该阶段平均海平面最高的时间为1962年，但其水体交换能力最弱，故可得出径流对水体交换能力影响较海平面的升高更为显著。

图4  1950—2018年七里海潟湖水体冲刷时间及对应的海平面、径流和主要人类活动

Fig. 4  Water flushing times and corresponding sea levels, runoff, and major human activities of Qilihai Lagoon from 1950 to 2018

2.2　第二阶段水体交换能力变化特征

第二阶段（1970—1979年），在地貌演变的模拟中考虑了潮汐通道口门处挡潮闸的建设（图3），且该阶段挡潮闸处于关闭阶段，图5为典型年份七里海潟湖涨急时刻的流场图，可以看出在1973年潟湖区域的流速较第一阶段的自然状态下显著降低，说明挡潮闸的建设对潟湖水体流入和流出有明显的阻碍作用，潟湖水动力大幅减弱，水体交换能力明显减弱。第二阶段平均径流量为6.45 m³·s^-1，与第一阶段6.99 m³·s^-1的径流量近似，但平均水体冲刷时间为259 h，是第一阶段平均水体冲刷时间的3.28倍，说明挡潮闸的建设显著降低了水体交换能力。同时，根据图4可以看出，1970年径流为6.65 m³·s^-1，较1969年径流（3.68 m³·s^-1）有大幅增加，但1970年的水体冲刷时间（247 h）较1969年（104 h）却增加了138%，水体交换能力总体呈现减弱趋势，说明挡潮闸的建设对水体交换能力的影响较径流更为明显。1973年的水体冲刷时间最大，为436 h，但年径流最小，为2.75 m³·s^-11；1979年的水体冲刷时间最小，为189 h，但年径流最大，为10.17 m³·s^-11，说明径流增加会增强潟湖的水体交换能力，此结论与第一阶段一致。该阶段的平均海平面变化范围为-0.18~-0.08 m，但是海平面的变化与水体冲刷时间的变化趋势相关性较差。因此，该阶段由于挡潮闸的建设对水体流入与流出的阻碍作用较大，对水体交换能力的影响远大于径流和海平面变化的影响。

图5  1960、1973以及2010年七里海潟湖流速分布

Fig. 5  Distribution of current velocity in Qilihai Lagoon in 1960,1973, and 2010

2.3　第三阶段水体交换能力变化特征

第三阶段（1980—1989年），由于环境及泥沙淤积问题，挡潮闸停止关闭，始终保持开启状态，且潟湖周围修建围堤使得水域面积从8.0 km²缩小至5.2 km²，平均水体冲刷时间由第二阶段的259 h降低为175 h，水体交换能力增强32%。根据图4，1980年径流和1979年径流接近，但水体冲刷时间从189 h降低至93 h，水体交换能力增加超过50%，说明挡潮闸的开启能够显著改善潟湖的水体交换能力。与自然状态（第一阶段）进行对比时发现，1952年和1988年径流均为4.58 m³·s^-1，但水体冲刷时间由1952年的88 h增加至141 h，说明围垦工程使潟湖的水体交换能力较自然条件下有明显的减弱。为了研究围垦面积与水体交换能力的关系，参考陈红霞等^［18］给出的纳潮量计算公式（公式（5）），进行了该阶段单位面积纳潮量的计算。

$$W=\frac{\mathrm{1}}{\mathrm{2}}(S_{\mathrm{1}}+S_{\mathrm{2}})(h_{\mathrm{1}}-h_{\mathrm{2}})$$

（5）

式中：W为纳潮量；S₁ 、 S₂ 分别为平均高、低潮潮位的水域面积；h₁ 、 h₂分别为S₁ 、 S₂所对应的潮高。

计算得到该阶段潟湖单位面积纳潮量由0.526 m³·m^-2降低为0.519 m³·m^-2，降低幅度仅为1%。纳潮量与围垦工程所造成的水域面积减小及海平面变化有关。据图4可知，该阶段的平均海平面变化微弱，纳潮量的改变主要来源于围垦面积的变化，也就是该阶段潟湖围垦对水体交换能力的影响有限，而挡潮闸的开启是水体交换能力改善的关键。该阶段的平均水体冲刷时间的变化范围为93~228 h，其中1980年水体交换能力最强，对应的年径流最大，为9.99 m³·s^-1；1983年水体交换能力最弱，对应的年径流最小，为1.33 m³·s^-1，故在相同的工况条件下，径流越大，潟湖水体交换能力越强，此结论与第一阶段一致。

2.4　第四阶段水体交换能力变化特征

第四阶段（1990—2007年）七里海潟湖内潮汐通道由弯曲变为顺直，潟湖进一步围垦，水域面积减小了2.9 km²，该阶段平均径流为2.55 m³·s^-1，平均水体冲刷时间由第三阶段的175 h降低至152 h，水体交换能力提升13%。据图4可以看出，虽然1990年的径流为2.28 m³·s^-1， 略小于1989年的2.86 m³·s^-1，但水体冲刷时间从179 h降低至137 h，水体交换能力增加23%，说明潮汐通道的变化显著增强潟湖水体交换能力。同时，由于该阶段继续进行潟湖围垦，海平面微弱上升，根据公式（5）计算该阶段单位面积纳潮量由0.519 m³·m^-2降低至0.474 m³·m^-2，降低幅度为9%，但水体交换能力并未随之减弱，说明潮汐通道的顺直较围垦工程来说对水体交换能力的改善更为显著。

综合第三和第四两个阶段，潟湖大规模的围垦均使得潟湖单位面积纳潮量减小，减小幅度分别为1%和9%，对水体交换能力的影响不及潮汐通道工况变化的影响，因此潮汐通道工况的变化是该阶段影响水体交换能力的关键因素。

2.5　第五阶段水体交换能力变化特征

第五阶段（2008—2018年）七里海潟湖潮汐通道拓宽，平均宽度由110 m增加至150 m，并在入湖口门处建桥，其中4个直径为2 m的圆形桥墩与水体接触，由于桥墩直径较潟湖口门宽度很小，故建桥对动力的影响较弱。根据图5c可知，该阶段水体的流入与流出有明显改善，水体交换能力增强。此阶段的径流条件趋于稳定，但海平面有了明显的上升，增加约0.15 m（图4）。由于潮汐通道的拓宽以及海平面的变化，潟湖的单位面积纳潮量为0.517 m³·m^-2，较第四阶段的0.474 m³·m^-2提升了9%，平均水体冲刷时间由第四阶段的152 h减小至126 h，水体交换能力提升17%。根据第一阶段得出的结论，海平面的变化对潟湖水体交换能力的影响较弱，故该阶段冲刷时间的降低主要来自于纳潮量的变化，即潮汐通道的拓宽提升了潟湖的纳潮量，有利于潟湖发生水体交换。对比2001年和2013年的情况，径流条件相同，均为0.82 m³·s^-1，但2013年的水体冲刷时间较2001年更小，水体交换能力更强，也就是潮汐通道拓宽增强了潟湖内的水体交换。

3 结论

本研究基于1950—2018年月径流量、海平面变化与七里海潟湖地貌演变的模拟结果，利用MIKE21建立七里海潟湖水动力模型，根据水体冲刷时间指标分析长时间尺度下潟湖的水体交换能力。研究结果表明：

在自然条件下，径流增加与海平面升高均可增强潟湖的水体交换能力，但径流对水体交换能力的影响更为显著；在相同工况的条件下，径流越大，水体交换能力越强。

受挡潮闸建设与潟湖围垦工程的影响，七里海潟湖的水体交换能力发生变化。挡潮闸建设阻碍了潟湖的流入与流出，降低了潟湖水体交换能力，但当其开启时水体交换能力得到大幅提升，因而挡潮闸的启闭成为影响潟湖水体交换能力的关键因素；同时，潟湖围垦降低了潟湖的纳潮量，但由于纳潮量变化幅度较小，减少水体交换能力有限。

潮汐通道的顺直和拓宽使得潟湖水体交换能力进一步增强。其中，第四阶段由于围垦工程所带来的纳潮量变化并未对潟湖水体交换能力造成显著影响，因而潮汐通道由弯曲变为顺直和拓宽是增强水体交换能力的主要因素。

作者贡献声明

匡翠萍：方法提出，论文撰写与修改。

王 洁：数值模拟，论文撰写。

董智超：数值模拟，图像绘制。

刘会欣：现场工作。

朱 磊：现场工作。