半封闭海湾是指海湾口门宽度与岸线之比在0.01~0.10之间的海湾^[1], 由于此类海湾的自身特殊性, 国内外很多学者对其进行了多方面研究。Lee等^[2]通过实测水文数据资料, 研究了围填海工程对韩国安海岸瑞山湾造成的影响; 李希彬等^[3]采用三维数值模型FVCOM探讨了湛江湾海域的水动力环境, 发现湛江湾内不同区域的水交换能力相差较大, 在出口处能力较强; 张雪松等^[4]基于三维水动力数值模型Delft3D, 模拟了唐岛湾纯净潮流和波流相互作用下的淤泥滩的演变, 并比较了不同工况下河道悬浮物和沉降和回淤强度, 为该海域海岸带的合理开发提供了科学的依据; 曾相明等^[5]利用POM模型研究了1963—2003年和1963—2010年围填海工程对象山港水动力环境的影响; 何杰等^[6]基于二维潮流数值模型模拟了深圳湾工程前、后潮流的运动, 发现工程的建设对深圳湾海域环境造成的影响主要集中在施工区域附近。随着沿海地区经济的飞速发展, 土地资源紧缺, 人类便通过围填海工程来满足经济和社会发展的需求^[7]。这在一定程度上缓解了土地供求紧张问题, 获得了巨大的社会和经济效益, 但也带来了生态退化、环境恶化、资源衰退等一系列问题^[8-12], 因此研究围填海工程对海域自然生态化环境造成的负面影响显得尤为重要。

本文以东山湾为例, 在实测地形和水文数据的基础上, 采用二维潮流数值模型MIKE21探讨了围填海工程对东山湾水动力环境的影响。东山湾口门宽度仅5 km, 是典型的半封闭海湾, 其地理位置特殊, 自然条件优越, 是福建重要的天然良港, 是发展海上交通、水产养殖、旅游及能源产业的优良场所^[13]。湾内海砂资源丰富, 同时古雷半岛沿岸非金属矿产资源丰富, 尤其是硅砂矿, 储量大、质地好。自20世纪90年代以来, 该海域的围塘养殖活动日益增多, 岸线也逐渐呈现向外推移的趋势, 此外还有一些其他围填海工程对该海域的自然生态环境造成了一定的影响。前人对东山湾水动力环境的研究较少^, 陈可锋等^[14]通过对实测资料的分析, 研究了八尺门海域的潮汐潮流特征, 同时根据模型结果计算了八尺门海域海堤工程对水道和整个东山湾流场的影响; 梁群峰等^[15]利用东山湾水文资料, 分析了其海底冲淤的变化特征, 表明在水动力较弱的区域海底冲淤变化较为明显; 吴煌荣^[16]从沉积动力特征入手, 分析了不同工况对八尺门水道海域环境的影响。随着该海域围填海工程的进行, 对其自然生态环境造成的影响研究显得越来越重要, 而海域水动力状况在一定程度上能反映该海域水体的自净和污染物扩散等能力, 所以对围填海工程前、后该海域水动力环境的研究就尤为重要。本文将从潮流动力、水体半交换和纳潮量三个角度来探讨围填海工程前、后东山湾水动力环境的变化, 为工程建设后该海域自然生态环境的维护和修复提供较为科学的依据, 并能够据此提出一些更有针对性的有效修复措施。

1 实测资料分析

为了解东山湾的水文条件, 分别于2013年9月1日—9月30日、2013年12月20日—2014年1月19日进行了秋、冬季两次潮位观测; 于2013年9月6日9时~9月7日11时(秋季大潮)、2013年9月9日8时—9月10日10时(秋季中潮)、2013年9月12日8时—13日10时(秋季小潮)、2014年1月2日9时— 3日11时(冬季大潮)、2013年12月30日9时—31日11时(冬季中潮)、2013年12月28日8时—29日10时(冬季小潮)共进行了6个潮次的潮流观测。潮位及潮流观测站点的位置如图 1所示。

1.1 潮汐

潮汐性质的划分主要依据(H₀₁+H_K1)/H_M2的值来判断^[17], 其中H₀₁、H_K1、H_M2分别为O₁、K₁和M₂分潮的振幅, $ \left(H_{\mathrm{O}_{1}}+H_{\mathrm{K}_{1}}\right) / H_{\mathrm{M}_{2}} \leqslant 0.5 $为正规半日潮, $ 0.5<\left(H_{\mathrm{O}_{1}}+H_{\mathrm{K}_{1}}\right) / H_{\mathrm{M}_{2}} \leqslant 2 $为不正规半日潮。根据潮汐调和分析的结果计算得知东山湾口的古雷、刺仔尾和东山湾顶的下寨三个潮位站的(H₀₁+H_K1)/H_M2值均处于0.5~2, 因此东山湾所处海域的潮汐属于不正规半日潮。

根据2013年9月和2013年12月—2014年1月两次实测潮位数据, 统计分析东山湾口的古雷、刺仔尾和东山湾顶的下寨三个潮位站的涨、落潮历时及主要潮汐特征值, 见表 1、表 2。由表中可以看出, 下寨站秋、冬季的涨落潮历时差都小于其他两站, 三个站的涨潮历时均大于落潮历时, 古雷站的历时差相对较大。在秋、冬季, 下寨站的最高潮位、最低潮位、平均高潮位、平均低潮位、最大潮差和最小潮差均大于古雷站和刺仔尾站。

1.2 潮流

和潮汐性质分类类似, 潮流性质由比值$ ({W_{{{\rm O}_1}}}{\text{ + }}{W_{{{\rm K}_1}}})/{W_{{{\rm M}_2}}} $来划分^[17], 其中$ {W_{{{\rm O}_1}}} $、$ {W_{{{\rm K}_1}}} $、$ {W_{{{\rm M}_2}}} $分别为O₁、K₁、M₂分潮流的椭圆长半轴的长度, $ ({W_{{{\rm O}_1}}}{\text{ + }}{W_{{{\rm K}_1}}})/{W_{{{\rm M}_2}}} < 0.5 $为正规半日潮流, $ 0.5 \leqslant\left(W_{\mathrm{O}_{1}}+W_{\mathrm{K}_{1}}\right) / W_{\mathrm{M}_{2}}<2 $为不正规半日潮流。潮流调和分析结果显示, 所有潮流观测站的$ ({W_{{{\rm O}_1}}}{\text{ + }}{W_{{{\rm K}_1}}})/{W_{{{\rm M}_2}}} $值均小于0.5, 说明研究海域属于正规半日潮流。

所有潮流观测站点中, 7#站实测涨、落潮流流速相差较小, 1#、3#、9#、11#站流速为最大, 且落潮最大流速明显大于涨潮最大流速。观测期间, 这4站实测最大落潮流速为159 cm/s, 出现在冬季大潮期间, 临近漳江落潮槽的8#站落潮流明显大于涨潮流。秋季期间10#站位于分流脊上, 水深较浅, 流速较小, 3个潮次的涨、落潮流最大流速相差不大, 且不超过40 cm/s。冬季该站调整至西南侧约1.6 km的水道上, 流速增大, 且落潮流最大流速明显大于涨潮流最大流速。12#站位于漳江入海口下寨码头附近, 离岸边较近, 实测涨潮流速大于落潮流流速。

1.3 余流

秋、冬季调查期间, 7#站余流流速较小, 秋季最大值为7.1 cm/s, 冬季最大值为13.4 cm/s。秋季调查期间, 湾内径流作用较冬季强, 从大潮至小潮, 各层余流流向朝逆时针方向偏转(NNW向至WSW向), 而冬季调查期间各层余流流向偏转角度较小, 主要集中在偏N向。1#、3#、9#、11#站为调查区径潮作用最明显的站位, 两者此消彼长, 径流作用强时, 潮流作用就相对减弱, 反之, 潮流作用增强时, 径流作用就相对减弱。临近漳江落潮槽的8#站, 余流流速要小于主槽内的4个站, 最大值仅为13.1 cm/s。秋季调查期间10#站水深较浅, 各层余流流速均较小, 最大仅为3.4 cm/s, 流向基本为W-N向。冬季期间该站调整至西南侧约1.6 km的主水道上, 流速明显增大, 且流向均指向落潮方向, 最大余流流速为21.7 cm/s(SSE向)。12#站各层余流流速均较小, 最大仅为8.9 cm/s。

1.4 围填海活动简述

此次围填海区域位于福建省漳州市漳浦县古雷半岛西面的近岸海域(图 2), 经纬度范围为23°50′40.58″N~ 23°48′37.17″N, 117°35′42.06″E~117°36′52.34″E。规划面积9 km², 其中填海造地6.77 km²(其余为陆地)。此次研究主要选取受围填海影响较为明显的潮流要素如潮流流场、纳潮量、水交换时间等作为指标。

2 建立数学模型

本文采用丹麦DHI公司研发的MIKE21软件, 它采用非结构三角网格剖分计算域, 网格设计灵活且可随意控制网格疏密, 可以用来模拟一维河口水动力系统以及二维的沿海区域, 并且能实现一维和二维区域之间自由水体的自由交换, 能够模拟各种作用力引起的潮位和潮流变化, 从而更准确直观的模拟复杂的近岸海域。

2.1 控制方程

质量守恒方程:

$ \frac{{\partial \zeta }}{{\partial t}} + \frac{\partial }{{\partial x}}(hu) + \frac{\partial }{{\partial y}}(hv) = 0. $

(1)

动量方程:

$ \frac{{\partial h\bar u}}{{\partial t}} + \frac{{\partial h{{\bar u}^2}}}{{\partial x}} + \frac{{\partial h\overline{v u}}}{{\partial y}} = f\bar vh - gh\frac{{\partial \eta }}{{\partial x}} - \frac{h}{{{\rho _0}}}\frac{{\partial {p_{\mathrm{a}}}}}{{\partial x}} - \frac{{g{h^2}}}{{2{\rho _0}}}\frac{{\partial \rho }}{{\partial x}} + \frac{{{\tau _{sx}}}}{{{\rho _0}}} - \frac{{{\tau _{bx}}}}{{{\rho _0}}} - \frac{1}{{{\rho _0}}}\left( {\frac{{\partial {s_{xx}}}}{{\partial x}} + \frac{{\partial {s_{xy}}}}{{\partial y}}} \right) + \frac{\partial }{{\partial x}}(h{T_{xx}}) + \frac{\partial }{{\partial x}}(h{T_{xy}}) + h{u_s}S, $

(2)

$ \frac{{\partial h\bar v}}{{\partial t}} + \frac{{\partial h{{\bar v}^2}}}{{\partial x}} + \frac{{\partial h\overline{v u}}}{{\partial y}} = f\bar uh - gh\frac{{\partial \eta }}{{\partial y}} - \frac{h}{{{\rho _0}}}\frac{{\partial {p_{\mathrm{a}}}}}{{\partial y}} - \frac{{g{h^2}}}{{2{\rho _0}}}\frac{{\partial \rho }}{{\partial y}} + \frac{{{\tau _{sy}}}}{{{\rho _0}}} - \frac{{{\tau _{by}}}}{{{\rho _0}}} - \frac{1}{{{\rho _0}}}\left( {\frac{{\partial {s_{yx}}}}{{\partial x}} + \frac{{\partial {s_{yy}}}}{{\partial y}}} \right) + \frac{\partial }{{\partial x}}(h{T_{xy}}) + \frac{\partial }{{\partial y}}(h{T_{yy}}) + h{v_s}S. $

(3)

式中: t为时间; x、y为笛卡尔坐标系空间坐标; η为水面高程, d为水深, h为总水深h=η+d; u、v为流速在x、y方向上的分量; f为科氏参数; g为重力加速度; ρ为水体密度; ρ₀为参考密度; p_a为大气压强; s_xx、s_xy、s_yx、s_yy为辐射应力分量; T_xx、T_xy、T_yy为水平黏滞应力; (τ_sx, τ_sy)和(τ_bx, τ_by)为水面和床底的切应力在x、y方向上的分量; S为源汇项流量; u_s、u_v为源汇项对应的速度分量。

2.2 模型初始及边界条件

1) 初始条件:

$ \left\{\begin{array}{l} \left.\zeta(x, y, t)\right|_{t=0}=\zeta(x, y)=\zeta_{0} \\ \left.u(x, y, t)\right|_{t=0}=\left.v(x, y, t)\right|_{t=0}=0 \end{array} .\right. $

(4)

2) 边界条件:

开边界: 开边界采用潮位预报边界条件:

$ \zeta {\text{ = }}{A_0} + \sum\limits_{i = 1}^{16} {{H_i}{F_i}\cos [{\sigma _{it}}t - {{({v_0} + u)}_i} + {g_i}]} . $

(5)

式中, A₀为平均海面, F_i、(v₀+u)_i为天文要素, H_i、g_i为调和常数。

外海开边界潮位由16个主要分潮(M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁, MU₂, NU₂, T_2, L₂, 2N, J₁, M₁和OO₁)调和常数推算得到。

闭边界:

在闭边界取流速的法向导数为0, 在潮滩区采用漫滩边界处理。

2.3 模型配置

1) 大范围模型

本文建立了大范围的福建南部海区二维潮流数学模型, 模型东边界至119.1°E, 西边界至116.0°E, 北边界到25.3°N附近, 南边界到22.8°N附近, 计算区域见图 3。

2) 小范围计算域

计算域网格采用非结构三角形网格, 通过网格生成模块, 控制网格疏密及尺度, 在拟建工程海域进行网格加密, 网格尺度最小为10 m左右, 因此可以更好地描绘工程区域的水下地形和海岸线, 同时可以保证精确的计算精度, 在外海区域, 为了提高计算效率, 网格相对稀疏, 为1 000左右。计算区域网格见图 3。

2.4 模型的验证

在本文中, 采用研究水域的潮位和潮流的实测资料来验证模型, 以检验模型的可靠性。观测资料包括2013年9月6日9时至7日11时的实测大潮资料, 以及2013年9月12日8时至13日10时的实测小潮资料。水文测站位置如图 4所示。

图 5给出了潮位的验证结果。从图中可以看出, 大潮期间最高、最低潮位误差一般在6 cm以内, 个别误差在10 cm左右, 误差相对较小; 相对而言, 小潮期间模拟潮位与实测潮位之间误差较大, 部分误差在20 cm左右, 但整体在可接受范围内; 总体来看潮位的模拟较为理想。

图 6进一步给出了潮流的验证结果。从涨落急时段的流向变化来看, 往复流特征较为明显; 各点流向的变化趋势拟合较好, 在部分转流和流速较小时刻, 流向跳动幅度较大, 此时流向验证误差较大, 各站点流向、流速的误差统计见表 3。从表 3中可以看出, 所有站点流向和流速的误差都相对较小, 大潮和小潮期间二者的误差大多在7%左右, 少部分转流和流速较小的时刻误差相对较大, 但总的来说, 结果与实测数据资料相差不大, 拟合较好。

从潮位和潮流的模拟结果来看, 模式计算结果能够较好地反映东山湾海域潮位、潮流的运动特征, 可进一步应用于东山湾海域水动力环境的研究。

3 结果 3.1 围填海工程前、后模型结果潮流变化分析

东山湾海域主要受到来自台湾海峡的前进波, 该波动在沿福建沿岸传播过程中传入东山湾并影响其水动力环境, 潮流形态以往复流为主, 如图 7所示。研究海域外侧近岸海域涨潮流路主要受福建和台湾沿海岸线的控制, 大体上呈WS-EN走向。大范围涨潮流偏EN向, 落潮流偏WS向, 涨落潮流路清晰, 往复流特征明显。

图 8给出了围填海工程前、后东山湾海域涨、落潮流的流场变化情况, 从图中可以看出: 涨潮时, 来自台湾海峡的涨潮流沿古雷半岛南段进入东山湾, 由于东山湾湾口相对狭窄, 涨潮流相对集中, 流速较大。涨潮流进入到东山湾后由于受到整体呈S-N走向的岸线控制, 东山湾内涨潮流整体上表现为S-N走向。围填海海域位于东山湾东侧, 同样受到岸线控制, 往复流特征较为明显, 涨潮流偏N向, 落潮流偏S向, 涨落潮流路清晰。同时由于东山湾湾内水深较浅, 尤其是湾顶, 大都为浅滩, 围填海区域水深也较浅, 因此围填海局部海域潮流动力较弱, 涨、落急流速在0.3 m/s左右。

图 9给出了围填海工程前、后局部区域的涨、落潮流场变化情况, 图 10给出了涨潮、落潮和全潮流速的变化值及变化率, 从图中可以看出: 从工程区域整体上来说, 围填海导致周边海域涨、落潮流速变化表现为有增有减:

(1) 围填海对工程区域南北两侧海域潮流的遮蔽作用比较明显, 涨潮动力减弱, 因此工程区域南北两侧水域均表现出不同程度的减小, 北侧水域涨潮平均流速减小范围在0.02~0.05 m/s, 涨潮流速减小率是5%~20%; 其南侧水域涨潮平均流速减小范围是0.04~0.14 m/s, 涨潮流速减小率为10%~70%。

工程区西侧水域, 由于工程在一定程度上减小了研究海域的过水断面, 尤其是西北侧堤头还对涨潮水流起到了一定的挑流作用, 流速有所增大, 平均流速增大0.02~0.06 m/s左右, 增大率为10%~ 30%。

(2) 工程区南北两侧水域落潮流速均表现出不同程度的减小, 工程区北侧水域落潮平均流速减小范围是0.03~0.06 m/s, 落潮流速减小率为5%~30%; 其南侧水域落潮平均流速减小范围是0.04~0.17 m/s, 落潮流速减小率为10%~70%; 工程区西侧水域落潮流速有所增大, 平均流速增大0.02~0.08 m/s左右, 增大率为10%~40%。

(3) 与涨落潮动力变化一致, 全潮平均流速减小的区域主要集中在工程区南北两侧与岸线所形成的潮流阴影处; 平均流速增大的水域集中于围填海西侧。其北侧水域全潮平均流速减小范围在0.03~ 0.06, 减小率在5%~30%; 其南侧水域全潮平均流速减小范围在0.04~0.17 m/s, 流速减小率在10%~70%; 西侧水域, 全潮平均潮流速增大0.02~0.08 m/s, 增大率在10%~40%。

总体来说, 围填海工程对全潮平均流速的影响范围在工程区域周边5 000 m范围以内, 变化较大的区域主要集中在工程区域附近。

3.2 纳潮量变化影响分析

由低潮到高潮港湾所能容纳海水的数量称为纳潮量, 它是海湾开发价值和生命力的一个重要参数^[11]。

本文纳潮量采用以下公式计算:

$ Q=1/2(S_{1}+S_{2})(h_{1}–h_{2}), $

(6)

式中, Q为纳潮量, S₁为高潮水域面积, S₂为低潮水域面积, h₁为高潮高, h₂为低潮高。

利用潮流数学模型, 计算在大潮汛时刻由最低潮涨到最高潮过程中通过断面的水体体积, 则为海湾的纳潮量。本次研究在东山湾湾口处布置一条能够封闭东山湾的断面, 同时在中间以及湾顶处布置类似的断面, 断面示意图如图 11所示。

计算结果如表 4所示, 可见, 本次围填海对东山湾的纳潮量影响为–2.5%左右, 但对围填海区以北纳潮量的影响较小, 纳潮量的损失主要是由于围垦造成滩涂减少引起的。

3.3 海湾水体交换影响

水体交换能力是评价海湾环境容量和环境质量的重要指标, 交换能力的强弱直接关系到海湾的水质状况。本文采用Luff等提出的半交换时间, 即某海域保守型物质浓度降为初始浓度一半所需的时间作衡量研究区域水体交换能力的指标^[18-19]。

假定计算域内外原水和新水都是均匀分布。计算域内的初始“表征浓度”为1, 计算域外新水的初始“表征浓度”为0。在交换过程中, 计算域内“表征浓度”逐渐变为0<$ \bar c $<1($ \bar c $为原水和新水混合后的平均“表征浓度”)。这里的半交换时间, 是指计算域内一个潮周期的平均“表征浓度”由$ \bar c $=1达到$ \bar c $=0.5的时间间隔, 计算网格与水动力计算网格相同, 水体交换方程为:

$ \frac{{\partial DC}}{{\partial t}} + \frac{{\partial UDC}}{{\partial x}} + \frac{{\partial VDC}}{{\partial y}} = \frac{\partial }{{\partial x}}\left( {D{k_x}\frac{{\partial C}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {D{k_y}\frac{{\partial C}}{{\partial y}}} \right). $

(7)

式中, C为示踪剂浓度, D为水深, k_x、k_y为水平紊流扩散系数。

选取东山湾为计算区域(图 11), 模拟该计算域内的水体与外部水体交换情况。假设计算域外水体为新水, 计算域内水体为原水, 浓度初始值为1, 东山湾湾外水域浓度初始值为0。涨潮时, 计算域外的水流入计算域内, 计算域内的原水和进入的新水混合; 落潮时, 计算域内原水与新水的混合水随落潮流带出计算域。浓度统计取样点如图 12中所示(白点)。

计算结果显示: 围填海实施前的水体半交换时间为220.5 h, 实施后的水体半交换时间为239.4 h, 可见围填海实施后半交换时间较围填海之前增加约18.9 h。

图 13给出了东山湾湾内表征浓度的分布图, 从图中可以看出, 围填海工程实施前东山湾中部及西北部海域整体上各浓度线都比围填海实施后靠北, 围填海工程实施后的海水交换程度没有之前海水交换充分, 即围填海活动在一定程度上会减弱东山湾和外海的海水交换能力。

4 结论

本文利用平面二维数值模型MIKE21, 对东山湾海域围填海工程前、后水动力环境变化进行了研究, 主要从潮流动力、纳潮量和半水体交换时间三个角度对其进行了分析, 结论如下:

(1) 围填海工程后局部涨、落潮流场发生了一定的变化。围填海实施后, 围填海区域南、北两侧的流矢变化较为明显。由于受到围填海区域南、北两侧东西走向岸线的影响, 这两个区域涨潮流矢由围填海前的偏N流向改变为围填海后的偏E流向, 而落潮时刻流矢由围填海前的偏S流向改变为围填海后的偏W流向。受到围填海区域岸线的遮蔽效应, 围填海区域南、北两侧水域流速也有一定的减弱。围填海海域西侧受到围填海西侧围堤影响, 涨、落潮流矢更加平直且集中, 流速有一定的增强。

(2) 此次围填海工程对东山湾的纳潮量影响为–2.5%左右, 但对围填海区域以北纳潮量的影响较小, 同时对东山湾底部红树林区域的纳潮量影响也较小。纳潮量的损失主要是由于围垦造成的滩涂减少引起的。

(3) 该围填海工程在一定程度上减弱了东山湾和外海的海水交换能力, 围填海工程实施后海水的交换程度没有围填海实施前海水交换充分。

通过对东山湾海域围填海工程前、后潮流动力、纳潮量和半水体交换的研究, 发现本次围填海工程对海湾总体的流场影响多集中在围填海工程区域附近; 工程对全潮平均流速的影响范围在围填海区域周边5 000 m范围以内; 同时, 对东山湾海域的纳潮量影响在–2.5%左右, 减弱了海湾水体的交换能力。研究表明此次围填海工程对东山湾海域水动力环境造成的负面影响主要集中在工程区域附近, 这为工程建设后该海域自然生态环境的维护和修复提供了较为科学的依据, 并且可以据此提出一些针对性的有效补救措施。本文未考虑工程的建设对周边海域冲淤环境和悬浮泥沙扩散等方面的影响, 该部分工作有待进一步研究。

致谢: 感谢张俊彪、巩明、俞亮亮等提出的宝贵意见。