围海工程前后瓯飞滩冲淤变化特征及动力机制分析

蔡家新，潘国富，陈培雄

(自然资源部 第二海洋研究所，浙江 杭州 310013)

0 引言

我国大陆岸线长约18 000 km，其中淤泥质海岸线约2 000 km，主要分布在长江三角洲、黄河三角洲、珠江三角洲及苏北平原海岸等区域[1]。这些地区多为重要的城镇、港口或渔场的所在地，为满足经济发展需要实施了大量围海工程[2-3]。围海工程在解决“人地矛盾”、拓展发展空间的同时，也会对周边海域的水动力环境及地形产生不同程度的影响，给海岸带资源与环境的可持续发展带来挑战[4-5]。

前人就围海工程对周边水动力的影响开展过许多研究，如童朝锋 等[6]通过数值研究分析了乐清湾内外围海工程对潮不对称偏度的影响；WANG et al[7]研究了江苏沿海在持续性围海作用下的潮滩水动力的响应机制，发现持续的围海使得潮滩变窄变陡，纳潮量显著减少，从而削弱潮上带和潮间带的流速，导致沉积物迅速累积；陈道信 等[8]利用二维数值模型，针对温州市在建、拟建及规划围海工程，预测了围海工程对近海及河口水动力的影响；穆锦斌 等[9]采用二维潮流模型研究了瓯飞围海工程对周边水动力的可能影响。

瓯飞滩位于温州瓯江河口和飞云江河口之间，是我国典型的淤涨型淤泥质海岸。近年来在该海域实施了多个围海工程，如瓯飞围海工程等，对区域水动力与地形产生显著影响。为正确评价围海工程对资源与环境的影响，本文根据实测地形数据对瓯飞滩围海工程实施前后冲淤变化进行分析，通过Delft3D数学模型模拟水动力环境变化，明确冲淤变化的动力机制，为海岸保护与利用提供科学依据和技术支撑。

1 区域概况

研究区域位于温州市瑞安市东部海域，介于瓯飞一期北片与瓯飞一期南片南堤之间(图1)，该区域海岸线呈NE—SW走向，近岸海底平坦，泥面标高为0.1～3.0 m，海涂坡度介于1/540～1/850之间，向东南倾斜，东部坡度逐渐变陡，滩涂淤积泥沙主要来自海域[10]。2005年以来在本区域内开展了一系列的围海工程，其中包括阁巷围海(工程Ⅵ，2005—2007年)、丁山二期(工程Ⅱ，2007—2009年)、平阳宋埠-西湾围海(工程Ⅶ，2009—2010年)、龙湾二期(工程Ⅴ，2012—2016年)、瓯飞一期北片(工程Ⅰ，2013—2017年)、丁山三期北片(工程Ⅲ，2016—2018年)、瓯飞一期南片南堤(工程Ⅳ，2017—2018年)等。其中，瓯飞一期北片、丁山三期北片和瓯飞一期南片南堤都是近几年完成的具有代表性的围海工程，是本文重点关注的对象。

图1 研究区域及测站点位Fig.1 Location of study area and observation stations

2 材料和方法

为使描述简洁明确，本文采用罗马数字指代相应的围海工程(图1)。本文重点关注的工程Ⅰ、工程Ⅲ和工程Ⅳ，均在2013年至2019年之间实施完成，总体规模较大，下文的围海工程指这3个工程。

2.1 地形资料来源与处理

本文收集到的研究区水下地形资料为2008年 10月、2010年10月、2014年10月和2019年3月的实测资料，其中2008年只测量了中高滩区域(工程Ⅱ、Ⅲ之间区域，水深小于1 m)，2010年测量了水深 2～10 m 的区域，2014、2019年均测量了水深小于 9 m 的区域。因工程Ⅲ已于2018年完成，2019年的测量没有包含该区。水下地形的测量依据《全球定位系统(GPS)测量规范》(GB/T 18314—2009)[11]和《水运工程测量规范》(JTJ203—2001)[12]。应用GPS RTK无验潮模式、DGPS有验潮模式结合测深仪测深，按断面法进行测量。

将地形资料统一到CGCS2000坐标系下，建立数字化高程数据库，采用克里金法(Kriging)进行插值，形成DEM高程模型。

计算2008—2013年、2014—2019年、2010—2019年的海床冲淤变化，并绘制冲淤图；利用各年份DEM高程模型，提取等深线并绘制等深线变化图。沿工程Ⅰ到工程Ⅳ岸段垂直方向设置4个地形横向断面，用于研究海底地形在垂向上的演变特征。其中断面1位于工程Ⅰ附近海域，断面2位于工程Ⅰ南侧海域，断面3位于工程Ⅱ附近海域，断面4位于工程Ⅳ附近海域(图1)。

DEM高程模型的建立、实测资料的处理和分析均通过ArcGIS 10.6软件进行。

2.2 潮流数值模型

基于Delft3D开源软件构建了瓯飞滩附近海域的二维潮流数值模型，控制方程、边界条件、数值离散方法等见软件技术手册[13]。

2.2.1 数值模型计算区域与网格划分

模型计算区域为温州近海，经度范围120°E—122°E，纬度范围26.5°N—28.5°N，涵盖工程Ⅰ、工程Ⅲ及工程Ⅳ。为保证核心研究区域处于计算区域的中心处，陆域边界上至玉环市，下至苍南县。模型对重点区域进行网格加密，网格大小为100 m×100 m，其他区域网格大小为500 m×500 m(图2)。河道边界采用流量边界，开边界处使用全球潮汐模型推算的M2、S2、K1、O1、N2、P1、K2、Q1八个分潮的调和常数作为驱动条件。

图2 模型计算区域网格及局部加密示意图Fig.2 Model computing area grid and local encryption grid

2.2.2 模型验证

数值模型计算时间为2019年3月10日到2019年3月25日，模型采用热启动方式运行，初始条件由预运行结果提供，时间步长为1 min。以同步的潮位和潮流实测资料作为率定、验证模型的基础数据，站位位置见图1。

2.3 潮不对称性分析

通过计算M2分潮和M4分潮的振幅比以及二者的相位差,分析潮的不对称性[14-15]：

(1)

G=2φM2-φM4

(2)

式中：A和G均为潮变形系数，A是M4分潮与M2分潮的振幅比；G是M4分潮与M2分潮的相位差。A、G潮变形系数的变化可以量化潮不对称的变化程度，A越大表明潮汐不对称性越强。当0°

3 结果与分析

3.1 冲淤变化时空分布特征

图3为2008—2019年研究区域等深线变化。2008—2019年，0 m等深线(滩涂区域)显著向海推进，在工程Ⅰ和工程Ⅱ围堤前沿，推进超过了500 m；2 m、3 m和5 m等深线呈现整体向陆推进的趋势；6 m 及以深海域，等深线持续向海推进。

图3 2008—2019年研究区等深线变化图Fig.3 Changes in isobath of the study area from 2008 to 2019

图4为2010—2019年研究区水下地形冲淤变化，近岸海域整体呈淤积态势。工程Ⅲ和工程Ⅰ前沿的两个凸区(a、g)，冲刷趋势最为明显，冲刷厚度约为0.8 m，而在工程围堤的凹区(b)和工程围堤与岸线形成的凹区(c)，水体扰动程度小，泥沙在此落淤，最大淤积厚度达到了6 m左右。

图4 2010—2019年研究区冲淤变化图Fig.4 Change of erosion and deposition in the study area from 2010 to 2019

图5为2008—2014年、2014—2019年中高滩区域(工程Ⅱ、工程Ⅲ之间，水深小于1 m)的冲淤变化。2008—2014年这一区域表现为有冲有淤，幅度在 0.1 m 以内，其中南侧区域以淤积为主(图5a)。2014—2019年，南侧区域(靠近工程Ⅱ)仍以淤积为主，中部略有冲刷(图5b)。

图5 研究区冲淤变化图(中高滩)Fig.5 Change of erosion and deposition in the study area(mid-high beach)

从4个断面的水下地形比较可以看出(图6)，冲淤变化呈现出明显的时空差异。2010—2014年(工程Ⅰ、Ⅱ、Ⅲ、Ⅳ完工前)地形变化相对较为平缓，但在2014—2019年(工程Ⅰ、Ⅱ、Ⅲ、Ⅳ完工后)，水下地形变动幅度较大。

图6 2008—2019年研究区断面地形变化图Fig.6 Topographic change map of cross sections form 2008 to 2019

断面1(工程Ⅰ附近海域)，2010年至2014年淤积平均厚度约为0.38 m，年平均淤积约0.09 m；2014年至2019年断面冲刷平均厚度约为0.14 m，年平均冲刷约为0.03 m。

断面2(工程Ⅰ南侧海域)，2010年至2014年略有淤积，年平均淤积约0.01 m，2014年至2019年断面冲刷平均厚度约0.1 m，年平均冲刷约为0.02 m。

断面3(工程Ⅱ附近海域)，2010年至2014年以淤积为主，平均淤积厚度约0.4 m。2014至2019年近岸处淤积，外海处冲刷。冲刷处冲刷厚度平均约0.3 m，淤积处淤积厚度平均约0.26 m。

断面4(工程Ⅳ附近海域)，2010年至2014年平均淤积厚度约为0.28 m，年平均淤积约为0.07 m；2014年至2019年工程Ⅳ与工程Ⅱ围堤夹角处明显淤积，中高滩最大淤积厚度约2 m。

3.2 水动力环境变化与分析

3.2.1 模型验证

图7和图8为代表性站位的潮位、潮流验证。由图可见，计算值与实测值整体吻合较好，本文构建的潮流数值模型可以较好地反映研究区附近海域潮汐和潮流特征。研究区海域的潮流流向主要为NW—SE向，往复流特征明显(图9)。

图7 瑞安、上干山站潮位验证图Fig.7 Comparison of observed and simulated tidal level variation of station Ruian and Shangganshan

图8 D1、D4和 D8测站流向、流速验证图Fig.8 Comparison of observed and simulation current speed and direction at station D1,D4 and D8

图9 研究区大潮涨、落急流场图Fig.9 Current at maximum flood and maximum ebb during the spring tide in the study area

3.2.2 潮汐特征模拟

研究区潮波来自西太平洋，以M2分潮为主，因此本文主要研究围海工程对研究区M2分潮的影响。选择2010和2019年作为代表年份，模拟了围海工程实施前、后研究区水动力场，分析工程实施后8个测站M2分潮的变化。

由表1可知，与2010年相比，2019年研究区海域M2分潮的振幅和迟角均发生了一定变化：(a)振幅增加了0.23～0.33 m，即围海工程实施后潮波上溯到同一区域时振幅增加。(b)D6测站迟角减小35.55°,其余测站迟角减小4.57°～6.60°，即潮波到达的速度减慢，同一区域达到最大潮位的时间推迟。

表1 M2分潮调和常数表Tab.1 Tidal harmonic constant of M2

表2为围海工程前后，各站潮变形系数的变化。与2010年相比，2019年，D1、D3测站振幅比(A)分别减少了0.05和0.02，；D2、D6测站未发生改变；D4、D5测站增加了0.01，D7、D8测站增加了0.02，表明工程实施后潮不对称性增强。

表2 潮变形系数变化表Tab.2 Tidal deformation coefficient

围海工程前后，研究区均表现为涨潮占优(G<180°)，但与2010年相比，2019年各测站相位差(G)增加了60°左右。

3.2.3 围海工程前后潮流特征变化

图10为围海工程前(2010年)和围海工程后(2019年)研究区涨急与落急流场变化(由于2010年工程Ⅰ尚未施工，因此2010年该位置仍有流场分布)。从图中可以看出，围海工程没有明显改变研究区整体潮流特征，但涨潮流方向整体有向逆时针方向偏转的趋势。围海工程对工程附近局部区域影响显著，围堤附近除流速的变化外，流向均不同程度转向沿岸方向。

图10 围海工程实施前后涨、落急流场变化图Fig.10 Distribution of current at maximum flood and maximum ebb before and after projects building

在工程Ⅰ围堤前沿(d区)，涨急时刻流速由 0.4 m/s 左右减小到0.2 m/s左右，流向由向岸方向转为西南偏西方向；在落急时刻流速由0.3 m/s左右增加到0.4 m/s左右，流向由离岸方向转为东向。在工程Ⅰ围堤与工程Ⅲ围堤所形成的凹区(b区)，涨急流速由0.6 m/s左右减至0.2 m/s左右，在两围堤形成的直角处形成流速为0.2 m/s的旋转流区域。结合水下地形冲淤演变特征(图4)分析，流速的大幅度减小是造成该区域淤积的主要原因。在围堤的凸区(a、g区)，涨急流速减少了0.2 m/s，落急流速增加约0.6 m/s，该区域附近出现了明显的冲刷。在工程Ⅳ飞云江河口北岸海域(e区)，受堤身影响，涨潮流流向由西北向转为西向及西南向，流速由0.6 m/s左右下降到0.2 m/s左右；落潮流流向由东南偏东转为东南向，即水流沿着工程Ⅳ堤身流动，流速由0.2 m/s左右增加到0.55 m/s左右。在工程Ⅳ围堤与丁山二期之间流速有明显下降，导致工程Ⅳ与岸线夹角区域的淤积。

在研究区中部区域(f区)，涨急流速减少了约 0.2 m/s 左右，流向由西北向转为西北偏西向；落急流速增加了约0.3 m/s左右,流向由东向转为东南向。鉴于研究区泥沙来源主要是海域来沙，且并未发生明显变化，因此落潮流流速的增加是导致该区域冲刷的主要原因。

4 结论

(1)围海工程实施前后，研究区泥沙来源未发生变化，涨潮占优略有增强，整体呈淤积态势；局部区域流速、流向明显改变，导致地形产生相应变化。

(2)在围海工程的凸区，潮流流速明显增加，区域冲刷明显，平均冲刷约0.3 m；在工程形成的凹区，流速明显降低，泥沙在此落淤。

(3)受工程Ⅳ影响，飞云江口落潮流流向由东南偏东转为东南向，导致泥沙运移路径发生改变。围海工程实施后，研究区中部海域落潮流流速增加。这两个因素导致研究区中部冲刷，冲刷幅度在0.1 m左右。