多工程和波流耦合作用下养护海滩的地貌演变特征

匡翠萍，韩雪健，宫立新，江林锋

（1.同济大学土木工程学院，上海200092；2.河北省地矿局第八地质大队，河北秦皇岛066001）

海岸带是极其重要的地理空间，大约20%的世界人口居住在距海岸线25km以内的地区［1］。因为沿海地区有利于诸多人类活动，如捕鱼、工业、旅游和运输等，因此人口聚集效应非常明显［2］，且沿海地区的人口增长趋势在未来的一段时间仍然持续。然而，海岸侵蚀影响人类生活、生物多样性、自然资源等［3–4］，是全世界沿海地区面临的一个严重问题。砂质海岸作为海岸带重要的组成部分拥有丰富的资源和秀美的风光，是人类休闲娱乐的聚集地［5］。砂砾质海岸占我国18 000km大陆海岸线的30%以上，其中七成以上的沙砾质海岸存在侵蚀现象，其强度和范围仍不断增大，采取相应的防护措施刻不容缓。

多种防护措施已经被应用到防止海岸侵蚀。这些措施经历了一个从硬到软、单一到复杂的过程。起初，为了防止海岸线的侵蚀，人们采用防波堤、丁坝、人工岬角等硬防护措施。硬防护在某些特定地区是可行的，但在其他地区可能失效并可能发生二次侵蚀的问题［6］。此外，硬防护也有建设成本高、工期长等缺点［7-10］。由于硬防护的这些问题，海岸侵蚀防护逐渐转向软防护，如海滩养护。海滩养护通常将海中或陆上的沙源用以填补侵蚀的海滩。与硬防护相比，软防护施工周期短，与当地环境兼容性强，对海岸环境变化的适应性强。然而，在近岸水动力的作用下，填充的泥沙会逐渐流失，因此海滩需要长期的定期维护。基于硬防护和软防护的利弊，最好的选择是充分发挥其各自优势，将两种措施有效地结合起来［11］。美国最早进行了海滩养护工程的实践，其数量和规模均居世界首位，累计的工程总数达2 000余项［12］。欧洲各国开始应用海滩养护技术的时间不一，其中葡萄牙于1950年首先付诸实践［13］。澳大利亚于1974年开始实施海滩养护工程，与欧美国家相比，其工程规模都较小但频率较高［14］。海滩养护在中国最早始于1992年的香港浅水湾，此后在天津、上海、大连、秦皇岛、青岛、平潭、泉州、厦门、珠海、北海、钦州、海口等35个城市均陆续实施，总数已达100余项［15］。秦皇岛从2000年开始实践人工养滩技术以来，陆续实施了十余次海滩养护工程，除滩肩补沙和水下沙坝补沙外，采取多种辅助设施，如鹅卵石护坡、预制沙袋潜堤、抛石潜堤、人工鱼礁、离岸岬头、沙质突堤、生态护岸等［16-20］。在多年的实践经验基础上，逐渐形成了“人工鱼礁—人工沙坝—滩肩补沙”人工养滩生态修复模式，既保护了海滩又使岬湾的水体交换问题得到改善［21］。

海滩养护工程在完成后受到水上、水下各种外部动力作用，其形态和物质组成的分布因此会发生变化，如海岸沙丘在风作用下发生运移及形态变化，水下海床在波流作用下发生冲淤变化从而改变水深。波流的作用相互影响：一方面，波浪通过辐射应力［22–23］的作用产生沿岸波生流和裂流等，同时波浪的存在会增大底部摩阻进而阻碍水流流动［24］、增大海表面粗糙度从而增加海气间的动量转化率等［25–26］；另一方面，共线的水流会使反方向传播波浪的波长减小、波高增大，使同方向传播的波浪波长增大、波高减小，不共线的水流会使波浪发生折射［27］，此外方向周期性变化的水流会使波浪的波高、波周期等波要素产生周期性变化的特征［28–29］。养护后的海滩泥沙组成丰富，且受到波流的共同作用，运动比较复杂，存在悬移质和推移质并存的运移形式。国内外学者对波流作用下海滩养护后工程区域的泥沙输运和分选过程［30–31］、地形演变［32–33］等进行了广泛的研究。

针对秦皇岛市戴河至洋河口岸滩侵蚀不断加剧，浴场海滩质量下降的问题，开展了基于海滩养护技术的海滩整治修复工程。工程（包含滩肩补沙、水下沙坝及人工岬头）对工程临近区域的地形地貌产生了一定的影响，同时离岸的人工岛也增加了该区域水动力环境的复杂性。本文基于验证良好的波流耦合水动力模型［34］、泥沙输运及海床演变模型，探究该区域在多工程组合影响和波流耦合作用下的地貌演变特征。

1 区域概况

秦皇岛市位于渤海湾与辽东湾交汇处，是著名的沿海城市。东北与辽宁接壤，西与唐山相邻，北连承德地区，南临渤海，万里长城横亘东西。西距北京270km，西南距天津220km，是北京的东大门，为东北与华北两大经济区的咽喉要道。随着区域经济的快速发展以及自然环境的变化，秦皇岛海域的岸线、湿地和生态环境受到不同程度的污染和损害，局部海域典型生态灾害频发，这大大影响了秦皇岛社会经济的可持续发展。秦皇岛市政府于2016年开始实施“蓝色海湾整治行动”，旨在使重点侵蚀岸滩得到有效整治修复、湿地生态环境得到明显改善。

戴河口至洋河口岸线整治修复工程作为“蓝色海湾整治行动”的重要组成部分之一，修复岸线长度约3.1km，以戴河口至仙螺岛约1.4km和洋河口至海上自助餐厅约1km两岸段为重点修复岸段，兼顾中间约700m岸段。工程内容包括滩肩补沙、水下沙坝吹填和人工岬头，如图1所示。该工程于2017年4月开始实施，至2017年6月完成全部工程结构单元，施工依照先人工岬头构筑，再滩肩补沙，最后人工沙坝吹填的顺序进行。

图1 戴河口至洋河口岸线整治修复工程平面布置Fig.1 Beach renovation and restoration project layout from Daihe Estuary to Yanghe Estuary

仙螺岛是由南戴河旅游开发区管委于1999年开发建成的一座供观光游览的人工岛，位于距岸边1km的浅海海域。仙螺岛平面形状为矩形，长100，宽70m，面积7 000m2［35］。

滩肩补沙后高程达到1.6m以上（以85高程为基准），两端共长约2.4km的重点养护岸段沙滩宽度拓宽40～60m，中间约700m岸段平均增宽10m左右。滩肩前缘向陆为坡度小于1：100的缓坡，向海坡度为1：10。滩肩补沙总方量约30×104m3。3座人工沙坝在离岸约200m处吹填而成，其中东北段布设2座，起到补充该岸段海滩泥沙以及在沿岸流作用下作为沙源对相邻岸段海滩进行补沙的作用，西南段布设1座沙坝，主要起到对此岸段海滩泥沙进行补充的作用。每座沙坝长约450m，坝顶宽约50m，底部宽约80～100m，坝顶高程为-0.9m左右。沙坝总计吹填方量约12×104m3。在戴河口西侧导流堤以南约100m外、洋河口东侧导流堤以东约100m外分别构建1座潜礁式人工岬头，形成人工岬湾，达到增加遮蔽作用，总体上又不影响观感和水体交换的效果。其中戴河口附近的人工岬头长约300m，洋河口附近人工岬头长约280m，两座人工岬头宽度均约12m，其顶部高程均为-1.2m左右。

2 数学模型建立

MIKE21软件由丹麦水环境研究所（DHI）研发，采用控制体积法对微分方程组离散，包含结构和非结构网格，主要用于海洋河流区域水动力环境的综合模拟，可模拟河流、海岸等二维水动力环境［36］。本文使用MIKE21软件建立基于水动力模型基础上的泥沙输运和海床演变模型，其中水动力模型验证良好并且应用模型得到了工程影响下的波流动力响应特征［34］，水动力模型相关的控制方程、参数设置及验证可参见［34］，这里不再赘述，仅对泥沙输运和海床演变模型进行详细介绍。

2.1 模型方程

（1）泥沙输运方程

研究海域泥沙运动以悬移质输运为主时，泥沙输运公式如下：

泥沙淤积和侵蚀公式如下：

淤积：

侵蚀：

式中：u、v分别为x、y向的水深平均流速；cˉ为水深平均泥沙浓度；εx、εy分别为x、y向的泥沙扩散系数；h为总水深；QL为源强，即单位面积的源流量；CL为源强的泥沙浓度；Ss为淤积/冲刷项；SD为淤积项，SE为冲刷项；ωs为泥沙沉降速度；Cb为近底床泥沙浓度；pd为淤积概率；τb为底床切应力；τcd为淤积临界切应力；E为侵蚀系数；a为软底床侵蚀指数；τce为侵蚀临界切应力；n为硬底床侵蚀指数。

（2）海床演变方程

海床演变方程的形式为：

式中：为泥沙干密度；ηs为床面冲淤厚度（正为淤，负为冲）。

海床演变模型是根据泥沙净沉积速率在每个时间步长内对海底地形网格数据做出更新，反馈到水动力模型，可确保水动力模型的精确性。

式中：Bn为n时刻的地形数据；Bn+1为n+1时刻的地形数据；Vn为n时刻海床冲淤变化。

海床厚度与地形一样，在同一时间进行更新。净沉积率表示为：

式中：DN为净沉积率；D为单位时间淤积厚度；E为单位时间冲刷厚度；Δt为时间步长。如果发生了净冲刷（DN＜0）同时M+DN＜0（M为河床可冲刷厚度），系统将自动调整到冲淤平衡，即既不冲刷也不淤积，也就是M+DN=0。

2.2 模型参数

泥沙输运和海床演变模型需要与水动力模型耦合计算，且共用一套网格。采用大、小双重嵌套网格进行计算。大模型网格为整个渤海海域，以大连老虎滩至山东烟台两个潮位站连线为开边界；小模型网格南起葡萄岛以南4 km处，北至金山嘴以北5 km处，其边界条件由渤海大模型提供。剖分后的网格如图2所示，大模型节点数有14 183个，三角网格数有23 419个；小模型节点数为10 988个，三角网格数为21 067个，坐标为北京54坐标系（中央子午线经度为120°E）。为了兼顾计算精度和计算效率，小模型网格由工程区向外逐渐稀疏，外海边界附近网格长度为800 m，工程区附近网格长度为10 m。在工程区域，沙坝、人工岬头和人工岛周围网格进行局部加密，以变水深刻画沙坝地形，将人工岬头设置为可变高度的堰模块，人工岛设置为陆地区域。泥沙输运和海床演变模型中沉降系数、临界起动切应力、临界沉降切应力均采用空间变化场。经过率定，沉降系数为5～10 m∙s-1；临界起动切应力根据唐存本公式计算得到，为0.06～0.4 N∙m-2；临界沉降切应力约为临界起动切应力的4/9，为0.036～0.178 N∙m-2。底床可冲厚度设置为2 m，泥沙初始浓度设定为0.02 kg∙m-3。三个开边界的泥沙浓度设定为0。

图2 模型计算网格及水深图Fig.2 The computational mesh and the bathymetry

2.3 模型验证

采用河北省地矿局第八地质大队实测泥沙资料对模型进行验证，实测泥沙资料采用2011年5—6月秦皇岛4个泥沙测站S1、S2、S3、S4的垂向平均泥沙浓度，测站位置见图3，验证结果见图4。

图3 泥沙测站位置示意图Fig.3 Locations of sediment observation stations

图4 2011年5－6月秦皇岛S1～S4站泥沙浓度验证及2016年10月戴河口站泥沙浓度验证Fig.4 Validation of sediment concentration at S1～S4 in May and June 2011 and at Daihekou Station in October 2016

采用Skill模型［37］对泥沙输运模型效率进行了评估，Skill值的范围在0～1之间。Skill值为1时，代表模型计算值和实测值之间完全相符；Skill值大于0.65时，表示模型计算结果为极好；Skill值在0.65～0.5之间时，表示模型计算结果为非常好；Skill值在0.5～0.2之间时，表示模型计算结果为好；Skill值小于0.2时，表示模型计算结果为差；Skill值为0时，代表模型计算值和实测值之间完全不相符。模型效率评估结果见表1，可以看出2011年5—6月及2016年10月5个测站的泥沙浓度在0.04 kg∙m-3左右，泥沙模型的Skill效率评价都在非常好以上。故可将其应用于工程建设数值模拟中。

表1 泥沙输运模型的Skill效率系数Tab.1 Skill efficiency coefficient of sediment transport model

2.4 波流耦合

地貌演变的数值模拟过程中，为了提高计算效率，通常对研究区域的水动力实测数据进行分析，统计分析出代表动力以表示该区域长时间的动力平均状况。代表潮主要有单一代表潮和多个代表潮两种选择方法。本文考虑到戴河口—洋河口岸线海域6月份平均潮位较高，潮差较大，海洋动力作用向岸延伸至滩肩边缘，对岸滩的塑造作用最强，故最终选择6月的自然潮型作为代表潮进行计算。波浪所引起的底部切应力对泥沙起动具有直接影响；波浪的波生时均剩余动量产生波生流，对泥沙输运具有间接影响。代表浪的选择需要同时考虑波浪在这两个方面的动力作用。波浪这两方面的动力贡献均与波高H的平方（反映波能）成比例，因此可以基于H2得到代表浪［38］。本文基于2016年实测波浪数据（期间发生了1次风暴潮），对波高的平方进行加权分别得到了代表浪的有效波高、周期和波向（见表2），其中常浪表示剔除风暴潮期间波浪序列的代表浪，强浪表示风暴潮期间波浪序列的代表浪。

表2 代表浪参数Fig.2 Parameters of representative wave

为了探究多工程组合影响和波流耦合作用下养护海滩的地貌演变特征，选取代表潮和代表浪耦合模式下对工程海域进行地貌演变数值模拟。一般来说海滩整治修复工程第1年的冲淤变化幅度最大，随后逐年减少［39］，因此选择1年为地貌演变的模拟时间。本文中由于选用参数化的代表浪来体现波浪作用，波流耦合采用波浪模型向潮流模型耦合的方式。也就是波浪模型计算所得到的波浪辐射应力加入到潮流模型计算中，从而影响潮流模型中的流场和潮位，再影响泥沙输运和海床演变。

3 结果分析

3.1 泥沙分布特征

根据验证良好的水动力模型可知［34］，工程区域的潮流为往复流，涨潮流为SW向，落潮流为NE向，流向基本与岸线平行，量级小属于弱潮区，相较而言波浪作用更为明显。涨急时刻的北突堤和落急时刻的南突堤背流侧（岬湾内侧）均出现流速极小的区域。常浪向为近SSE向，强浪向为近ESE向，且常、强波浪均在人工岛处发生绕射并在背面产生波影区。

选择具有代表性的夏季大潮分别与常浪和强浪耦合对工程后戴河口—洋河口海域泥沙浓度分布进行模拟分析。图5为夏季大潮与常浪耦合下工程后戴河口—洋河口海域涨落急时刻的泥沙浓度场。从整体分布上看，戴河口—洋河口海域的泥沙浓度较低。工程后沙坝靠岸侧泥沙浓度场在涨落急时刻存在明显的高浓度带，涨急时刻由于洋河口导堤的阻挡作用，高浓度带与工程前类似，均位于导堤东北侧，但更靠近岸边，且浓度可达0.07 kg∙m-3。落急时刻高浓度带依然纵向贯穿整个戴河口—洋河口近岸海域，但与工程前相比也更靠近岸线，且都位于水下沙坝与岸线之间海域，浓度约为0.08～0.1 kg∙m-3，落急时刻高浓度带平均泥沙浓度高于涨急时刻。水下沙坝处泥沙浓度较高，涨急时刻约为0.05 kg∙m-3，落急时刻洋河口处的沙坝可达0.08 kg∙m-3，落急时刻水下沙坝平均泥沙浓度高于涨急时刻。沙坝处及沙坝靠岸侧区域泥沙浓度较高是由于沙坝位于岸滩前，阻挡了波浪的直接侵袭，波浪先作用在沙坝上，掀起了较多泥沙，并向位于其后的岸滩输运。波浪是该区域泥沙悬扬的主要动力，在近SSE向波浪的作用下，涨急时人工岛后北半岸线水域的波影区波浪动力较弱，加上SW向的水流搬运作用，使得该区域没有出现悬沙高浓度带。而在落急时，波浪辐射应力驱动的NE向沿岸波生流和同方向的落潮流将南半岸线附近高浓度悬沙向NE方向搬运，使得悬沙高浓度带纵向贯穿该岸线。

图5 夏季大潮与常浪耦合下工程后涨急和落急时刻的泥沙浓度场Fig.5 Sediment concentration fields at the maximum flood and ebb after the project under the coupled summer spring tide and the normal waves

夏季大潮与常浪耦合情景下，将工程后泥沙浓度场减去工程前的泥沙浓度场，得到图6。可以看出工程后工程海域的泥沙浓度有所增加，特别是水下沙坝所在位置，涨急时刻最大增值可达0.04 kg∙m-3，落急时刻最大增值可达0.07 kg∙m-3，落急时刻大于涨急时刻。水下沙坝的构筑使得波浪在沙坝处破碎，破波紊动引发大量的泥沙悬扬，破波产生的波生流和潮流叠加产生向岸水流，引起沙坝之后的区域泥沙浓度增加。

图6 夏季大潮与常浪耦合下工程引起的涨急和落急时刻泥沙浓度变化场Fig.6 Sediment concentration difference fields at the maximum flood and ebb due to the project under the coupled summer spring tide and the normal waves

图7是夏季大潮与强浪耦合下工程后戴河口—洋河口海域涨落急时刻泥沙浓度场。泥沙浓度分布工程后与工程前相比很相近，区别在于工程后强浪下的泥沙浓度值整体增大，悬沙高浓度带向外海移动至离岸200 m的位置。涨急时刻最高浓度位于洋河口导堤堤头处，约为0.15 kg∙m-3，落急时刻最高浓度位于洋河口侧沙坝前，约为0.17 kg∙m-3，落急时刻高浓度带平均泥沙浓度高于涨急时刻。与常浪情况下不同，强浪下的悬沙高浓度带处于沙坝向海侧，这是由于更大的波高使得波浪在到达沙坝前已经破碎，破波紊动引发大量泥沙悬扬。涨急时刻在接近ESE向的强浪作用下，人工岛的西侧产生波影区，使得悬沙高浓度带在该处中断。在落急时刻，由于落潮流的搬运作用，中断位置向北移动，同时由于强浪下泥沙不易沉降，南突堤南部的高浓度悬沙向北运移，受到人工岛的阻碍而收缩，形成南突堤和人工岛之间的高浓度带。

图7 夏季大潮与强浪耦合下工程后涨急和落急时刻泥沙浓度场Fig.7 Sediment concentration fields at the maximum flood and ebb after the project under the coupled summer spring tide and the strong waves

图8为夏季大潮与强浪耦合下工程引起的涨落急时刻戴河口—洋河口海域泥沙浓度变化场。强浪作用下工程对泥沙浓度的影响与常浪下类似。工程后，工程海域的泥沙浓度有所上升，增加最大处位于水下沙坝附近，涨急时刻增大值可达0.043 kg∙m-3，位于戴河口西南侧水下沙坝后。落急时刻增大值可达0.065 kg∙m-3，位于洋河口水下沙坝处。落急时刻的泥沙浓度增加值高于涨急时刻。

图8 夏季大潮与强浪耦合下工程引起的涨急和落急时刻泥沙浓度变化场Fig.8 Sediment concentration difference fields at the maximum flood and ebb due to the project under the coupled summer spring tides and the strong waves

3.2 地貌演变特征

利用海床冲淤演变模型对整治修复工程后的戴河口—洋河口海域进行1年的海床演变模拟。图9为潮流和常浪耦合下工程前后研究海域1年的海床冲淤演变，相当于全年不发生风暴潮的情况。工程前情景下，南半岸线附近水域发生比较严重的侵蚀，最大侵蚀深度达到0.3 m。北半部分岸线附近水域由于原始地形具有沿岸多条沙坝的特征，因此出现与岸线平行的条状侵蚀、淤积交错的现象，在坝顶侵蚀，最大侵蚀深度达0.3 m，在槽底淤积，最大淤积厚度达0.15 m。由侵蚀范围在岸线中段的收缩现象来看，人工岛在一定程度上缓解了其后侧的侵蚀，但是对离岸200 m范围内的缓解作用不明显，滩肩区域仍然出现比较明显的侵蚀。工程后情景下，三座人工沙坝处发生明显侵蚀，沙坝1、沙坝2和沙坝3处侵蚀深度最大分别达到0.53 m、0.59 m和0.63 m。沙坝向岸侧的区域除靠近水边线的少部分滩肩区域外出现明显淤积现象，平均淤积厚度为0.1 m，最大淤积厚度达0.18 m，由此可见沙坝有效地保护了其后侧的掩护区域，且沙坝在波浪作用下侵蚀的泥沙向后侧输运喂养了掩护区域。南部人工岬头1在接近SSE向常浪作用下其后侧工程前情景下的侵蚀区域消失，北部人工岬头2与沙坝间的区域侵蚀范围减小、量级降低，由此可见人工岬头对于海床侵蚀亦有一定的抑制作用。

在常浪中耦合1 d的强浪，计算出常-强浪与潮流耦合情景（相当于1年发生1次风暴潮）下1年后戴河口—洋河口海域工程后的冲淤演变，见图10。工程前情境下，与图9相比外海出现大面积的微侵蚀区域，沙坝周围区域侵蚀范围变大，程度加深。在工程后情景下，与图9相比三座人工沙坝处侵蚀程度加深，沙坝1、沙坝2和沙坝3处侵蚀深度最大分别达到0.73 m、0.8 m和0.81 m。沙坝后侧掩护区域的侵蚀面积变大，淤积面积变小（其中滩肩处的侵蚀面积变大，程度加深），但靠近沙坝的后侧淤积厚度变大，沙坝1、沙坝2和沙坝3后侧淤积厚度分别达到0.2 m、0.39 m、0.27 m。可见，三座沙坝处均发生侵蚀，其后侧淤积，沙坝出现向岸的运移现象。也就是在风暴潮存在的情况下虽然滩肩侵蚀现象变严重，但沙坝及人工岬头掩护区域依旧存在大面积的淤积区域，因此工程仍具有一定的保护作用。

图10 潮流和常－强浪耦合作用下1年的海床冲淤演变Fig.10 Annual seabed evolution in 1 year under the coupled tides and normal-strong waves

4 结论

本文基于验证良好的波流耦合水动力模型建立泥沙输运及海床演变模型，并结合现场实测资料验证了模型的可靠性，进而探究该海域在多工程组合影响和波流耦合作用下的地貌演变特征，得到以下结论：

（1）戴河口—洋河口海域泥沙浓度整体较低，潮流和常浪耦合下的泥沙分布与潮流和和强浪耦合下的泥沙分布形式一致；平均泥沙浓度及工程前后泥沙浓度增加值落急时刻大于于涨急时刻。

（2）工程后沙坝周围出现悬沙高浓度带。常浪时由于沙坝位于岸滩前，阻挡了波浪的直接侵袭，波浪先作用在沙坝上，掀起了沙坝泥沙，向位于其后的岸滩输运，高浓度带位于沙坝向岸侧；强浪时由于更大的波高使得波浪在还未到达沙坝处已经提前破碎，破波紊动引发大量泥沙悬扬，高浓度带位于沙坝向海侧。

（3）工程后沙坝处发生侵蚀，沙坝向岸侧发生淤积：常浪与潮流耦合作用下（代表1年内不发生风暴潮）沙坝处侵蚀深度最大达0.63 m，沙坝向岸侧最大淤积厚度达0.18 m；常、强浪与潮流耦合作用下（代表1年内发生1次风暴潮），沙坝处侵蚀深度最大达0.81 m，沙坝向岸侧最大淤积厚度达0.39 m。可见，工程在有无风暴潮的情况下均能起到保护海滩的作用，达到修复海岸的目的。

作者贡献说明：

匡翠萍：项目及基金负责人，研究方法，模拟方案，数据分析，论文修改；

韩雪健：数值模拟，数据分析，论文撰写；

宫立新：现场实测，论文撰写；

江林锋：数值模拟，数据分析，论文撰写。