沿海地区是人类活动的枢纽^[1], 半封闭海湾作为其中一部分, 由于优越的地理位置, 对社会经济发展有更为重要的意义^[2]。20世纪50—70年代, 为促进沿海农渔业发展, 扩大社会经济效益, 实施了大量围填海工程^[3], 导致大量海湾水体恶化, 而半封闭海湾自净能力较弱, 湾内富营养化程度加剧^[4]。海湾环境主要通过涨落潮引起的潮汐交换进行控制, 通过交换进程稀释湾内富营养化程度^[5], 故水体交换强弱是评价海湾水环境优劣的重要标准之一。而潮汐通道(tidal inlet)作为海湾与外海之间的重要连接通道, 使海湾与外海之间的水体能够实现高效交换, 对湾内水环境至关重要^[6]。以往学者为提升半封闭海湾水体交换能力, 大多采用拆除现有构筑物的方式^[7], 通过潮汐通道改善自净能力的研究并不多。而可控潮汐通道技术是指可以改变潮汐通道潮动力交换水平的一系列研究手段, 包括闸控、疏浚、开挖涵洞等。因此, 深入探究可控潮汐通道技术在海湾水体交换中的作用, 对更好地理解和管理海湾生态系统具有重要意义。

评价不同修复方案对海湾内水交换能力影响时, 野外监测和数值模拟是最常用的方法。然而, 由于对长期监测数据的巨大需求以及修复工程效果的不确定性, 实施前后水质之间的比较仅限于有限范围内的区域和监测期^[8]。因此, 就综合扩散模型和时间计算而言, 数值模拟可能是评估半封闭海湾恢复工程水动力学和水质响应的有效替代方法。目前研究海湾水体交换能力的数值方法主要为欧拉法^[9]和拉格朗日法^[10]。欧拉法是基于固定的坐标系描述物体运动的方法, 而拉格朗日方法则是基于跟随物体移动的坐标系描述物体运动的方法。众多学者基于这两种方法, 从滞留时间、纳潮量、半交换周期、水龄、人类活动^[11–13]等多因素分析半封闭海湾内水体交换情况; 针对潮汐通道, 研究更多集中于演变规律、地貌影响、潮波运动规律及流场特征变化等方面, 对水体交换的影响却鲜有研究。

柘林湾是典型的多汊型潮汐汊道^[14], 自20世纪70年代起, 湾内陆续围垦, 湾内水域面积减少近一半, 围内大量排放污水, 导致水生生态系统恶化^[4]。目前国内对海湾的研究仍比较有限^{[15, 16]}, 学者们在柘林湾的研究大都集中于养殖水污染评估^[17]、富营养化特征分析^[18]等方面, 而以何种技术理念解决湾内现存水质问题暂未见过相关研究报道。故本文以柘林湾为例, 基于不规则网格TELEMAC模型, 模拟潮流动力, 并应用WAQTEL模块, 通过示踪剂模拟污染物的扩散过程。提出了一种基于潮汐通道改善湾内水体交换速率的方法, 为日后封闭、半封闭海湾生态环境恢复工程的规划和设计提供一定的理论依据。

1 研究区域 1.1 地理位置

柘林湾地理位置为23°32′~23°37′N, 116°57′~117°05′E, 位于广东省潮州市饶平县海域, 如图 1。柘林湾是典型的多汊型潮汐通道, 潮流属典型的不规则半日潮。自20世纪80年代以来, 随着海水养殖规模不断扩大, 柘林湾水质污染问题日益严重, 造成这一现象的主要原因是20世纪70年代末建设的“澄饶联围”工程。沟通联围内外水域的三百门水闸和高沙水闸“只排不纳”, 使联围内部的劣质水体长期排入柘林湾内, 严重影响湾内水质^[19]。

1.2 水体环境

根据实测, 分别以二到四类的标准对目标海域的超标因子、超标样品数等指标进行评价。表 1表明湾内最主要的超标因子是无机氮和活性磷酸盐, 以第二类标准评价时, 最大超标倍数分别达到0.99、4.10。此外, 以一定的降解系数求得总磷(TP)、总氮(TN)以及化学需氧量(COD)含量。图 2为TP、TN、COD的分布示意, 可以看出, 湾内TN、TP含量为湾外开放水域的3~5倍, 并于三百门水闸附近富集。

2 研究方法

目前, 众多学者面对日益严峻的海湾环境, 提供了一些可行性的措施。其中, Xiong等^[20]研究了切萨皮克湾与邻近海岸陆架之间、海湾内不同区域之间的水交换, 以及它们与径流、风和持续时间的关系, 发现海湾内水交换时空变化主要与淡水流量和风高度相关, 并且可以通过Michaelis-Menten/Monod方程和线性递减函数的组合来描述; 王平等^[7]针对普兰店湾内的水动力特征进行讨论, 研究不同修复方案对水交换的影响, 表明通过拆除近岸养殖池, 可以增加湾内10%的纳潮量, 湾内半交换周期降至50 d; 佘小建等^[21]改善马銮湾地区水体环境时, 研究了海湾口门闸门开闸孔数、开闸位置对湾内水动力和水体交换的关系, 通过闸门调度使湾内形成大范围的环流提高水体的交换效率等。

可以发现, 目前对近岸半封闭海域水体交换的研究依然很少。即使已经采用了一部分修复措施, 但方法仍局限于拆除人工构筑物与增加径流量来增强海湾水动力, 这些修复方案在一些大型人工海水养殖区会受到明显限制。此外, 修复措施对半封闭海湾的影响尚未形成完善的量化体系, 尤其是环境改善的有效性评估存在严重不足。本文针对研究现状, 通过增设闸门、疏浚水道改善潮汐通道, 以一种对现存生境和人工构筑物影响较低的方式, 促进湾内水体交换, 并量化交换效果, 合理评估其对湾内环境的影响。

2.1 模型简介与设置

TELEMAC (http://www.opentelemac.org/)是一款研究水动力学和水文学领域的高性能数值仿真开源软件, 本文基于TELEMAC建立二维水动力与示踪剂输运数学模型, 数值计算方法采用有限元法或者有限体积法, 是一套适用于模拟河流、河口和海岸二、三维水动力、泥沙、水质和生态等问题的模型系统^[22]。本文采用有限元法, 与有限体积法相比, 有限元对非结构化网格和复杂的边界问题效果更好, 能够提供较高的数值精度, 在实际应用中, 有限元法在TELEMAC系统中的广泛使用得到了验证, 并在涉及复杂水体流动和海洋建模的问题上表现良好。TELEMAC-2d求解浅水方程, 计算网格每个节点的主要结果是水深和深度平均速度分量, 控制方程如下:

$ \partial h/\partial t + \nabla \cdot \left( {h\boldsymbol{u}} \right) = 0, $

(1)

$ \frac{{\partial u}}{{\partial t}} + \boldsymbol{u} \cdot \nabla \left( u \right) = - g\frac{{\partial \eta }}{{\partial x}} + \frac{1}{h}\nabla \cdot \left( {h{v_t}\nabla u} \right) - \frac{{{\tau _{\text{b}x}}}}{{\rho h}}, $

(2)

$ \frac{{\partial \boldsymbol{v}}}{{\partial t}} + \boldsymbol{u} \cdot \nabla \left( v \right) = - g\frac{{\partial \eta }}{{\partial y}} + \frac{1}{h}\nabla \cdot \left( {h{v_t}\nabla v} \right) - \frac{{{\tau _{\text{b}y}}}}{{\rho h}}, $

(3)

其中, t是时间, h是水深(水面到海底的距离), (u, v)是速度矢量u的分量, $\nabla $是笛卡尔坐标系中的空间微分算子, g是重力加速度, η是水面高程, ${v_t}$是扩散系数, ${\tau _{\text{b}x}}$和${\tau _{\text{b}y}}$是BSS(床层剪切应力)矢量${\tau _{\text{b}}}$的分量, ρ是水的密度。

示踪剂和污染物的输运可以通过WAQTEL建模

$ \frac{{\partial hC}}{{\partial t}} + \nabla \cdot \left( {h\boldsymbol{u}C} \right) = \nabla \cdot \left( {h{v_c}\nabla C} \right) + {S_{\text{tr}}}, $

(4)

其中, C是示踪剂的浓度; ${v_c}$是示踪剂的扩散系数, 本研究中初始值为1.0, $ {S_{\text{tr}}} $是示踪剂浓度的总来源, 由外部来源(工业、污水处理厂等)以及物质降解或与其他物质的化学反应引起的浓度变化等共同组成。

建立大范围数值模型时, 开边界潮位条件采用“TPXO 9.0”潮汐模型^[23], Telemac-2d可以直接读取TPXO二进制数据库并生成实时潮汐层边界条件。CFL数(库朗数)是数值计算中用来衡量数值稳定性的重要参数, 描述了数值解法中波传播速度和网格离散步长之间的关系, 模型采用半隐式时间迭代, 在满足CFL < 1.0的条件下, 计算时间步长为10.0 s。使用非结构三角形网格剖分计算域, 覆盖范围从塘边湾至南门湾, 计算使用的投影坐标系为UTM50 (EPSG: 32650)。此外, 考虑黄冈河对湾内影响, 根据《广东省水文图集》推算黄冈流域多年平均流量为39.26 m³/s。

2.2 水体交换

为量化研究方案实施后的交换效果, 合理评估其对湾内环境的影响, 通过示踪剂浓度的时空分布来反映湾内、外水体的交换情况。给定湾内示踪剂的初始浓度假定为${C_0}$, 某一时刻变成了$ {C_1} $, 此时湾内水被外海水置换的比率为$R = \left( {{C_1} - {C_0}} \right)/{C_0}$, 相应余留在原位置没有被置换的水体比率为$L = 1 - R = {C_1}/{C_0}$。

此外, 利用Luff等引入的半交换周期概念^[24], 表征水体与外界的平均交换能力。湾内水被外海水置换比率R达到0.5时, 所花费的时间称为半交换周期。为了便于分析柘林湾内不同区域之间水交换差异, 将研究区域分为6个部分, 记为S1—S6。示踪剂初始浓度在图 8所示的蓝色框中设置为1。

$ \bar C = \frac{{\mathop \sum \nolimits_{i = 1}^n \;{V_i}\left( t \right) \times {C_i}\left( t \right)}}{{\mathop \sum \nolimits_{i = 1}^n \;{V_i}\left( t \right)}},t = \left[ {0,T} \right] , $

(5)

$ \bar R = \frac{{\mathop \sum \nolimits_{i = 1}^n \;{V_i}\left( t \right) - \mathop \sum \nolimits_{i = 1}^n \;{V_i}\left( t \right) \times {C_i}\left( t \right)}}{{\mathop \sum \nolimits_{i = 1}^n \;{V_i}\left( t \right)}},t = \left[ {0,T} \right], $

(6)

$ {V_i}\left( t \right) = {S_i}\left( t \right) \times h\left( t \right), $

(7)

其中, $ \bar C $是不同区域内示踪剂的平均浓度, ${S_i}$、$ {V_i} $、h分别为每个网格计算单元的面积、体积与当地水深, $\bar R$用来表征区域范围内示踪剂的流动速率。

3 模型验证

为验证模型在柘林湾海域的准确度, 选取3个潮位测站点(T1、T2、T3)和3个潮流测站点(L1、L2、L3), 依据2016年1—2月与6—7月的水文测量资料对模型结果进行率定, 站点位置如图 3。测点的潮位和潮流的选取位置, 综合考虑了潮流由外海向近岸的传播过程。本研究高程数据均基于1985年国家基准面。

测点的潮位和潮流验证结果如图 4、图 5所示。黑色实线为模型潮位, 红色空心圆为实际潮位。潮位为不正规半日潮, 潮差近3 m, 潮位零点为平均海平面。从潮位的验证结果来看, 模型计算的各站位的潮位值与实测值吻合程度良好。

4 结果与讨论 4.1 流速流场

大、小潮涨、落急流场如图 6所示。模拟结果表明, 涨潮时, 海山岛西侧的涨潮流北向运动, 流至澄饶联围附近潮流分为东西两股, 西股向义丰溪河口汇入, 东股流入笠港水道。海山岛东侧涨潮流沿后江水道向东北方向运动, 到达柘林湾口门位置转为北向, 主要由大金门水道和小金门水道进入湾内, 之后水深变浅, 潮流方向较为散乱。落潮时情况完全相反。

流速分布由海向岸逐渐减小。海域的流速最大值位于后江水道, 达到0.96 m/s; 流速较大处还包括大、小金门水道和笠港潮汐通道西南侧海域, 小门水道内流速稍小。大小潮落急时的流速显著大于涨急时的流速, 其中小金门水道的差异最显著, 大潮涨急流速比大潮落急低0.3~0.4 m/s, 而小潮落急时的流速普遍大于小潮落急时的流速0.2 m/s。另外, 三百门水闸东北侧海域受到掩蔽, 流速随时间变化波动小, 保持在0~0.2 m/s。

4.2 余流场

余流场采用欧拉余流计算方法, 对空间固定点上的潮流流速矢量和在潮周期取平均得到, 反映通过某一空间点的所有物质所表现出来的整体净输移趋势, 如图 7。

$ {u_{\text{Euler}}} = \frac{1}{T}\int {_Tu\text{d}t} , $

(8)

式中, $ {u_{\text{Euler}}} $为欧拉余流的速度, T为计算潮周期, u为某一时刻的潮流流速。

大金门水道余流方向与涨潮流向基本一致; 小金门水道内余流场较为紊乱, 与南澳岛与汛洲岛之间的复杂岸线有关, 但南侧可见明显的落潮方向流场; 小门水道余流整体表现为落潮方向; 笠港水道包括其南部湾内区域的余流也表现为落潮方向。

整体而言, 柘林湾内的复杂潮汐通道引发湾内形成明显的环流结构, 潮流主要从大金门水道向湾内输入, 进入湾内后分为东西两侧的回旋流, 东侧的环流集中于西澳岛北侧, 通过小门水道落潮向余流向湾外输出。而柘林湾西侧余流表现为三百门水道的归槽流, 引起柘林湾西南角的环流, 湾内物质通过环流从小金门水道与笠港水道向外排出, 由此可见笠港水道在一定程度上控制着柘林湾西部湾内水体与物质的净输出。这种影响从余流场角度看应该集中于三百门水道西南侧, 湾内其余部分水流的交换由小金门水道与小门水道主控。不过以笠港水道目前的通量, 三百门水道北侧水体无法与湾外进行有效交换, 工程前污染物质将随环流在此聚集, 水质较差。

4.3 水交换影响 4.3.1 闸控

“闸控”技术拟在联围两端改建闸门, 当外海潮位高于围内水位时, 开启三百门水闸引水, 利用围内水道蓄水, 当外海潮位低于围内水位时, 开启高沙水闸排水, 形成稳定单向流, 实现“纳排联动”模式。围内地势较低, 只有–0.11 m, 水位最高0.39 m, 而堤外平均高潮位1.314 m, 为了使湾内半交换周期缩小至10 d, 设计纳潮量为2.80×10⁶ m³/d, 为简化工况, 采取平均径流量60 m³/s, 经计算, 日纳潮量最大值可以达到3.00×10⁶ m³/d, 最小纳潮量也能满足2.78×10⁶ m³/d, 满足预期目标, 从而恢复潮汐通道, 改善湾内水环境。图 9为工程前湾内污染物浓度变化情况, 表明柘林湾内水体交换能力在空间上分布不均匀, 大致以口门为中心, 交换能力与口门距离呈负相关关系, 并以扇形趋势向湾内蔓延。其中, S5、S6区域水交换速度最快, 10 d便能达到75%的交换率。随后是S1与S3, 交换率最差的是S2、S4区域, 10 d仅能交换30%左右的水体。5~10 d, 存在径流以及潮汐通道影响的海域(S1, S3), 与封闭海域(S2, S4)相比, 污染物浓度明显降低, 水体交换百分比能够提高20%。湾内半交换周期(污染度含量下降至50%的时间)在没有径流、潮汐影响的封闭区域, 主要集中在10~30 d; 而存在潮汐交换的区域, 半交换周期显著减小, 集中在2~10 d。模型结果证明, 改善潮汐通道对半封闭海湾的水交换提升作用显著, 笠港潮汐通道与三百门、高沙水闸的疏通与改建具有理论支撑。

图 10为工程后湾内污染物浓度变化状况。其中S1、S6区域浓度变化最小, 与工程前相比(图 11), 各时刻浓度差值均位于5%以下; 区域S2水交换提升效果最为显著, 浓度下降10%~15%; S3、S4及S5区域变化幅度相差不大, 保持在5%左右, 一定程度上减少污染物浓度的富集, 加强水体的交换进程。

闸门改建之后(图 12与图 13), 湾内整体水域半交换周期从10~30 d减少至5~15 d, 闸门附近的半交换周期减幅最大, 时间减少15 d, 交换速率为工程前的2倍。湾内的其他区域半交换周期减幅略低, 集中在5 d附近。值得注意的是, 黄冈河口门附近的半交换周期略微增加, 可能由于三百门开闸影响此处流态, 加剧回旋流。不过该处水交换速率快, 污染物仍能在10 d内完成半交换, 影响忽略不计。

为合理评估潮汐通道改善后各时刻的提升效果, 对柘林湾各区域不同时期水交换程度进行比较, 结果如图 14。研究表明S5和S6由于初始的几日内便能交换70%, 导致从模拟的5~60 d, 整个水交换进程都没有太大增幅, 每个阶段仅提高5%左右; S1、S3、S4区域在5~10 d区间内水交换能力提升最快, 25%以上; 而S2区由于处于湾内最深度, 略微滞后, 10~20 d内水交换能力显著提高, 增加32%。工程之后各个区域能在20 d之内完成60%的水体交换, 有效改善S2区域污染物滞留的状态。

4.3.2 联控

“联控”指闸门改建与笠港潮汐通道疏浚, 二者联合控制潮汐通道交换水平。闸控模拟结果表明, 闸门“纳排联动”(三百门水闸纳、高沙水闸排)模式, 恢复了潮汐通道的作用, 极大改善初始“只排不纳”导致湾内水环境恶化的状况。图 10表明, 同一时间段笠港水道交换程度最弱, 可能由于笠港水道(S4区域)淤积严重, 年均淤泥强度超过0.25 m/a, 导致通道内潮汐交换速率始终滞后。

借鉴闸门重建增设潮汐通道改善湾内水环境的思路, 对笠港水道整体疏浚2 m。并将两种方案实施后的模拟结果与基准(联控前)进行对比, 进一步探讨潮汐通道对半封闭海湾交换能力的影响, 为后续提升封闭、半封闭海湾水交换水平提供参考。

模拟表明, 闸门与疏浚协同作用下水交换提升显著。图 15与只进行闸门工程的结果相比, 发现笠港疏浚后S1—S4交换水平均有所提升。其中S1、S2提升1%~2%, 幅度较低; S3和S4交换效率得到显著提高, 尤其是S4区域, 涨幅最大能够达到10%。研究结果表明: 增强潮汐通道作用的方案能够极大降低半封闭海湾内污染物浓度, 仅进行简单疏浚便能提升湾内交换水平达15%。

值得注意的是, S5在5 d左右, 水交换差值为负, 可能是由于潮汐通道排出水体时, 夹杂了大量污染物, 这一部分污染物存留S5内无法及时排出。随着时间的增长, 10 d左右该区域便可以恢复正常交换。此外, 工程前后各区域的水交换差值在30 d附近达到峰值, 说明工程实施后一个月左右效果最佳。20 d左右, 6个区域的水体均能够实现至少70%的水体交换, 与柘林湾一个完整大小潮周期相近, 能够有效降低湾内污染物浓度。

针对交换速率涨幅最大阶段(20 d), 比较不同方案下湾内交换能力, 如图 16所示。在未增设任何手段的基准状况下, S2与S4区交换始终处于偏低水平; 改建闸门后S2区得到改善, 提高13%; 笠港潮汐通道疏浚后S4区域交换能力提高10%左右; 而两个方案协同作用下, 整体水质进一步提升, 20 d左右便交换70%的水体, 有效控制外源对围内污染, 减少围内污染压力, 同时促使湾内潮水流动, 改善湾内海水生态环境。

5 结论

本文在总结现有研究成果的基础上, 为改善湾内水环境提出新的思路。运用TELEMAC建立平面二维数学模型, 依据实测资料对柘林湾的潮位和潮流进行了验证, 分别计算余流、半交换期、水体交换率, 结论如下。

(1) 柘林湾是典型的多汊型潮汐通道, 湾内大面积浅滩, 涨潮流基本由槽向滩方向运动, 沿西北向汇入; 落潮流由滩向槽(潮汐通道)运动, 方向相反。涨潮初期水体从南澳岛东侧通过3个汊道进入柘林湾内, 汊道处流速增大, 涨潮时流速0.35 m/s, 落潮达到0.7 m/s, 落潮流速大于涨潮流速。同时, 湾内的复杂潮汐通道引发湾内形成明显的环流结构, 导致湾内水体交换时空上分布不均匀。

(2) 闸门改建前后相比, S1、S6区域变化最小, 各时刻水交换涨幅均位于5%以下; S2区域水交换提升效果最为显著, 提升10%~15%; S3—S5区域变化幅度相差较小, 保持在5%左右。三百门水道附近半交换周期减少近20 d, 交换速率为工程前的两倍。各区域水体能在20 d之内完成60%的交换, 一定程度上减少污染物浓度的富集, 加强交换进程。笠港水道与闸门纳排协同作用下, 湾内水环境得到进一步改善, 20 d水体完成至少70%的交换。

(3) 总结而言, 增设闸门潮汐通道一定程度上增强了湾内水动力, 完整大小潮期间, 能够提升水体交换速率15%左右, 有效改善半封闭海湾的水环境。不过闸门纳潮量的增加导致围内水动力提高, 对联围内外沿岸自然景观和生态环境的影响尚不明晰。此外, 其对当地生态环境与沙滩资源的保护等存在的影响, 是否会干扰本地物种的成功定植^[25]以及闸门纳排联动导致大量污染物流出口门, 这些物质对南澳岛生境的影响等问题都有待进一步研究。

勘误

《海洋科学》2022年第46卷第5期，第122-130页李叶, 杜以帅, 徐建平等人发表的论文《清洗频率对固定床生物滤器水处理性能的影响》一文中第124-126页图 2、3、4，因作图时误差棒有误，应作者要求进行更正。为此，特将正确的图 2—图 4刊登如下。有关网站的电子版本已及时进行替换。