粉砂淤泥质潮滩对于波浪作用的响应

汪 垚，谢 婕

（中交上海航道勘察设计研究院有限公司，上海 200120）

引言

潮滩的演变对海岸防护有着重要的意义，其演变机制的研究已成为了多学科关注的焦点[1-2]。潮滩水动力、泥沙以及地貌演变对波浪过程较为敏感，波浪对于潮滩演变起着重要的作用，是维持潮滩稳定的重要因素，在开敞潮滩的数值模拟中应考虑波浪的作用[3]。

一般认为，在多因素的共同作用下，潮滩剖面形态发育为上凸形和上凹形两种[4]。在泥沙供应充分条件下，潮流动力强而波浪作用弱，潮滩往往发育为上凸形剖面；而当波浪作用强且潮流动力较弱时，潮滩往往发育为上凹形剖面[5]。国内外对于砂质海岸的研究较早[6-7]，研究成果也相对丰富，但对于粉砂淤泥质潮滩的研究成果仍然较少。潮滩上动力条件复杂，现场恶劣的自然环境很大程度上限制了持续的现场观测，亟需通过数值模拟的方式研究潮滩演变的动力地貌过程。目前数学模型研究中所考虑的动力因素也以潮流为主[8-9]，Prichard等[10]以欧洲西北部的潮滩为原型，建立了概化的一维横向剖面的数学模型，通过一系列的数值实验，研究了在单一横向潮流作用下，潮滩剖面形态与潮差、潮汐不对称性、大小潮以及沉积物供应的关系，并指出，潮滩的宽度与潮差无关，与沉积物的供应呈正相关关系；涨落潮不对称使得潮滩坡度变陡，落潮占优还将导致潮滩的蚀退，大小潮的变化也会减小潮滩的淤积速率。该模型虽然对实际情况作了较大的简化，但已经可以反演出潮流控制的潮滩的主要形态特征，不过，该模型的缺陷之处在于其忽略了波浪对于泥沙运动的影响，而Dyer[11]和Kirby[12]等人指出，波浪和潮流作用之间的平衡控制着实际中大多数的潮滩系统，波浪在泥沙起动中的作用尤为重要。

现有研究成果中波浪对潮滩演变的动力机制研究关注较少，Waeles[3]等人建立了考虑表面重力波和潮流的一维概化横向剖面数学模型。计算结果表明，单一潮流作用下的潮滩会持续淤长，而潮流、波浪共同作用下的潮滩可以达到稳定，为了进一步研究波浪的作用，分别在外海边界计算波高为5 cm、10 cm的两种工况，前者的剖面仍能保持上凸形剖面特性，但后者潮滩的剖面变为上凹形。这一研究表明，即使在较小的波浪作用下，潮滩的剖面仍会受到波浪的影响，同时也通过数值模拟的手段表明波浪对于潮滩塑造的作用不可忽略。Rossington[13]在前人所建立的模型的基础上，加入了沿岸潮流的影响因素，建立了一套考虑横向潮流、沿岸流以及波浪等因素的一维横向剖面数学模型，与实际情况更为接近，该模型被运用于英国塞文河口潮滩的季节性冲淤变化的计算，并模拟了该潮滩的均衡态剖面，与实测值吻合良好。但该模型的缺陷在于，沿岸向潮流仅仅增加了床面切应力，并不输运泥沙。因此，该模型在运用时仍存在一定的局限性，其仅适用于沿岸向净输沙很小的潮滩。可见，一维横向剖面模型虽然取得了较大的成果，但其局限性也十分明显，模型中无法反映出沿岸方向的水沙输运，对于沿岸方向水沙输运占主导的潮滩（例如江苏中部的大丰潮滩）并不完全适用。总的来说，目前在二维模型方面开展的研究仍然较少。

1 研究区域

大丰潮滩为典型的堆积型粉砂淤泥质潮滩，宽度约为2～6 km，坡度平缓，约为0.1 %～0.3 %，近岸潮汐类型为不规则半日潮，平均潮差3.68 m，冬季有效波高小于1 m，其他季节小于0.5 m[14]。近岸潮流为与岸平行的往复流，涨潮历时普遍小于落潮历时，平均涨落潮历时之比为0.73，涨潮流速一般大于落潮流速，两者之比约为1.4。总体来说，研究区域动力作用以潮流作用为主，波浪作用较弱[15]。潮滩在南北向的涨落潮流的影响下，泥沙输运以沿岸向的输沙为主。在充足的泥沙供应下，剖面形态呈现明显的双凸形（见图1），分别在高潮水边线和平均低潮位线附近形成凸点[16]。滩面的沉积物粒径范围为0.001～0.05 mm，以粉砂和淤泥为主，由于滩面高程、动力条件的差异，滩面沉积物呈现出明显的分带性，自岸向陆可划分为泥滩、泥沙混合滩、粉砂细砂滩[17]。

图1 观测断面和剖面各点高程

2 模型建立

本文建立了平面二维水动力、泥沙输运及潮滩中长期演变数学模型，在微时间步长内耦合计算水流、泥沙运动和地形变化，实现中长期地貌演变过程的高效模拟[18]。波浪计算采用SWAN波浪模型，控制方程为动谱平衡方程，采用二维作用量谱密度N(σ,θ)来描述随机波浪场，模型采用有限差分方法，在时间、空间和谱空间五个变元对控制方程进行数值处理。对非粘性沙组份采用van Rijn方法计算悬沙和底沙运动[19-20]；对粘性沙组份采用Partheniades-Krone公式[21]计算悬沙运动，忽略底沙运动。

本文所选择的研究对象——大丰潮滩位于西洋海域，此处的岸线基本平直，且涨落潮流方向与岸线平行，因此，根据江苏中部大丰王港河口潮滩遥感图像，将研究区域概化为沿岸向12 km、横向10 km的矩形纳潮盆地。参照2006年实测的JD33、JD34剖面（图2，采用国家85高程，以下同），假设初始横剖面高程从2 m向海线性降低至-14 m，初始坡度为1.6 ‰。

综合考虑模型稳定性（CFL＜10）和计算效率的要求，模型中水流、泥沙以及地貌的时间步长Δt=0.5 min。由于波浪计算量较大，为提高计算效率，波浪计算的时间步长为Δt=745 min，也即为与水流模块的耦合时间步长，且每次均在高潮位时计算波浪，经过敏感性测试，满足模型计算精度要求。根据大丰港潮位站2006年9月～2007年10月实测潮位资料，经调和分析计算得到分潮调和常数。

图2 2006年实测的JD33、JD34剖面

图3 观测断面及观测点位置

模型的东边界考虑了M2、M4、S2和MS4分潮，其振幅分别为1.7 m、0.2 m、0.6 m、0.1 m，相位分别为0°、-142°、-295°、-141°，并考虑了潮波在沿岸方向的相位差。在南、北侧向开边界，采用Nuemann水动力条件，即保持沿岸方向水位梯度为常数[22]。研究区域波浪作用较弱，研究区域盛行偏北风[23]。模型设置波浪方向为NNE方向，为分析不同波浪方向对于潮滩演变的影响，而研究区域南部受到辐射沙脊群的掩护，因此，模型中以E向为对比工况。波浪大小由边界上的有效波高控制，分别为0.1 m、0.15 m、0.2 m三种工况。波浪频谱周期均取为3 s。模型中忽略风应力和径流的影响。

江苏中部的大丰潮滩为典型的粉砂淤泥质潮滩，2006年Y9站涨潮时悬沙d50=0.01 mm，粒径超过0.062 mm泥沙组分不超过5 %。2008年观测的潮间带悬沙中淤泥含量超过85 %。悬沙中主要以粘性沙为主，非粘性沙主要以推移质的形式输移，因此，模型中同时考虑粘性沙和非粘性沙。初始底床由5 m的粘性沙和5 m的非粘性沙组成，不考虑床面分层情况，泥沙均匀混合，当床面发生冲刷时，粘性沙和非粘性沙单独计算冲刷量，且不考虑沙、泥间相互作用。经调试，粘性沙沉速取为0.6 mm/s，冲刷临界切应力为 0.2 N/m2，冲刷率系数取5×10-4kg/(m2·s)，淤积临界切应力取为较大值[24]，保证粘性沙颗粒始终处于沉降过程；非粘性沙中值粒径d50=90 μm。研究表明，江苏中部的外海泥沙供应充分，根据2006年现场测量结果[25]，西洋水道涨落潮垂线平均最大含沙量约1.0～1.6 kg/m3。经调试，模型东边界上粘性沙含量取为1.0 kg/m3，南北侧向开边界含沙量从0.1 kg/m3向东边界作线性插值。三个边界上非粘性沙均采用Nuemann含沙量条件，即边界处的含沙量取值为其内侧相邻网格的含沙量。

3 潮流和波浪作用下潮滩演变模拟

3.1 波浪场分析

选择波向为NNE（22.5°），波浪谱峰周期T=3s，边界波浪Hs=0.10 m、0.15 m、0.20 m三种工况，对比分析演变初期（自初始地形开始演变30天后）潮滩上的波浪变化。

图4为三种波浪工况下的波浪场，波浪在向岸传播过程中，由于床面摩擦的作用，能量耗散，有效波高不断减小，同时，波浪的方向也发生偏转，尤其是在高滩区域，波浪方向几乎与岸线垂直。

图4 计算区域波浪场（有效波高单位：m）

前人的研究表明[26]，波浪对于泥沙运动的影响主要通过增大床面切应力，增强床面扰动作用，波浪引起的切应力计算公式如下：

式中：

τw——波浪引起的床面切应力；

ρ——水体密度；

fw——底部摩擦系数；

uw,b——波浪底部质点轨道最大流速；

Hs——波浪有效波高；

Tp——波浪谱峰周期；

k——波数；

h——水深。

由公式（1）可知，波浪引起的切应力大小主要与波浪质点轨道的最大流速有关，且τw u2w,b。因此，uw,b的大小可反映出波浪对于床面扰动作用的强弱。而由式（2）可知，uw,b又与水深h有关，且在一定范围内，随着水深的减小而增大。这使得波浪在向岸传播过程中对底部的扰动作用发生变化。图5为当边界波高Hs=0.1 m时，演变30 d后断面B上uw,b横向变化（均为高潮位时）。可以看出，大小潮时的uw,b变化趋势一致，波浪产生的床面切应力向岸方向先增大，后迅速减小。在西洋深槽由于水深较大，波浪对床面几乎没有扰动作用，随着水深向岸逐渐减小，uw,b的大小向岸增加，从而导致床面切应力增加，可以看到，波浪的扰动作用几乎影响整个潮间带。而当水深减小到一定程度，波浪发生破碎，波高迅速衰减，从而导致uw,b迅速减小，波浪对于床面的扰动作用急剧减弱并消失。小潮时由于水深较浅，波浪在离岸较远处破碎，无法到达潮间带上部区域。

以大潮为例，比较不同波浪大小条件下的uw,b的横向变化分布，如图6所示。不同波浪大小条件下，uw,b的变化趋势一致。波高越大，波浪所产生的近底扰动越大，且随着波浪的增大，波浪破碎点的位置向外海移动。总体而言，本模型对于潮滩上波浪过程的模拟结果与前人的定性认识一致，可用来进一步模拟分析波浪对于潮滩水沙过程以及剖面形态的影响。

3.2 水沙过程分析

选择初始地形条件下离岸2.5 km且处于潮间带的P1点，对比不同波浪工况下模型演变初期（自初始地形开始演变30 d后）的水沙过程，包括水位、流速、床面底部切应力和含沙量过程，由于模型中波浪作用较弱，因此不同波浪工况下的水位和流速过程差异很小，图7仅列出了不同波浪作用下的底部切应力和含沙量过程。图8为仅考虑潮流作用下高潮位及低潮位下的流场图。

图5 断面B上uw,b 横向变化（Hs=0.1 m）

图6 不同波高条件下断面B上uw,b 横向变化对比

图7 演变初期P1点的水沙过程对比

图8 计算区域流场图（流速单位：m/s）

床面切应力过程随着波浪作用的增大明显增大，且在小潮时波浪的作用更为明显，这主要是因为小潮时水深较浅，波浪对于床面的扰动更强烈，对床面切应力的贡献也就越大。含沙量过程与床面切应力过程的变化趋势相一致，两者存在较好的对应关系。随着波浪作用的增强，含沙量明显增加。可以认为，由于波浪作用导致的床面切应力增加，从而导致泥沙再悬浮作用增强，悬沙中含沙量增加。同样地，小潮时含沙量的增加更为明显，原因与床面切应力的变化一致。总的来说，波浪作用对于水位和流速的影响不大，主要的作用在于增大了床面的切应力，增强了泥沙的再悬浮作用，悬沙中含沙量增加。根据计算区域流场图可以看到，西洋深槽内涨落潮主要为N-S方向的往复流，在潮间带为较弱的逆时针旋转流，与前人的现场观测结果基本一致[27]。

表1、表2分别列出了边界有效波高Hs=0和Hs=0.1 m/s两种工况下计算域北、南、东开边界断面粘性沙悬沙通量、水通量及计算域内净通量，统计时段为地貌演变30 d后的一个大潮、小潮。可以看出，悬沙通量与水通量变化保持同步；沿岸方向通量为横向通量的数倍，即沿岸向水沙输运相对于横向输沙占优[28]。涨潮优势特性使得北边界净通量为流入，而南边界净通量为流出计算域。大、小潮的水沙通量对比，由于小潮水动力弱，水通量和悬沙通量总体减小，大潮时悬沙净通量为负值，表明潮滩整体处于冲刷，小潮时悬沙净通量变为正值，表明潮滩整体上处于淤积。当考虑波浪作用后（表2），大小潮的冲淤趋势一致，但大潮时潮滩整体的冲刷量有所增加，小潮时淤积量略有减少，总体上潮滩整体的淤积量略有减小。（表1、表2中斜杠前数值代表大潮，斜杠后数值代表小潮；流入、流出计算域通量分别以正、负值表示。）

表1 粘性悬沙通量及水通量统计（Hs=0）

表2 粘性悬沙通量及水通量统计（Hs=0.1m）

3.3 波浪对地貌演变的影响

在沉积物供应充分的条件下，潮流的作用使得泥沙净向岸输运，潮滩不断淤长，并形成双凸形剖面。由水沙过程分析可知，波浪的作用增大了床面切应力，增强了泥沙再悬浮作用，增强了床面的冲刷作用。本节将重点讨论不同波浪条件下潮滩剖面形态的变化。

1）考虑波浪前后的地貌演变比较

图9为计算区域模拟2年后考虑波浪与否的地貌演变对比。总体上对比来看，无论是否考虑波浪作用，整个计算区域均处于淤积环境，但考虑波浪作用后计算区域的淤积量明显减小，尤其是在离岸约3～6 km的潮间带区域。图11为距离南边界6 km的断面B在不考虑波高以及波高为0.1 m两种工况下每1年的剖面地形演变的对比。对比剖面形态的演变过程，可以看到，在较小的波浪作用下，潮滩演变的趋势未有变化，整个潮间带依然处于淤积环境中，西洋深槽处于轻微冲刷状态，潮滩总体仍然能够保持淤长的状态。从各局部位置来看，与仅考虑潮流作用相比，除了潮间带上部的沉积量略有增大外，剖面其他位置的沉积量均明显减小，这主要是由于波浪的作用使得破波点附近及其向海一侧的冲刷作用增强，沉积量减小，而涨潮流占明显优势的潮流作用将更多由波浪悬浮起来的泥沙挟带至潮间带上部堆积，且潮间带上部的波浪已经破碎，对泥沙的再悬浮作用较弱，几乎无法使泥沙起动。正是由于波浪增强了破波点附近及其向海一侧的床面泥沙的再悬浮作用，导致了上凸点向海淤进的速率明显减小，下凸点向上抬升的速率减小。而西洋深槽的地形差异不大，始终处于轻微冲刷的状态，主要是因为西洋深槽水深较大，受波浪的影响较小。对比波浪场的分析，可以发现地貌演变的结果与波浪横向切应力变化分布规律相一致。波浪的作用增大了潮间带的冲刷，减弱了泥沙向岸输移作用，但对深水区以及潮间带上部影响不大。总的来说，与仅考虑潮流作用相比，在较小的波浪作用下，剖面形态整体上更接近于实际现场观测的结果。本研究所建立的模型能够较好地复演出江苏中部潮滩的剖面形态特征（见图1）。

2）波高对于剖面形态的影响

图11为模拟2年后三种不同边界波高的潮滩剖面形态对比。当边界处有效波高增大至0.15 m，剖面由原来的淤长型发展为侵蚀性剖面，形态由上凸形转变为凹形，离岸2.5 km以外的区域始终处于冲刷状态，并在2.5 km处形成一陡坎。边界处有效波高继续增大至0.2 m，冲刷作用增强，陡坎位置向陆方向移动。模拟的剖面形态与Roberts等人[29]一维模型的结果定性吻合。

图9 计算区域模拟2年后的地貌演变（高程单位：m）

图10 考虑波浪与否的断面B剖面地形每1年的变化对比

江苏中部海岸波浪条件的季节性变化明显，冬季波浪相对较强，其它季节波浪较弱。因此，为更接近实际情况，在模型中考虑波浪的季节性变化，分析季节性波浪对潮滩剖面的影响。模型中设置12月-次年2月的波浪为冬季浪，波高为0.15 m，其他季节波高为0.1 m，波向均为NNE向（暂未考虑波向变化，波向影响见3.3.3），分析演变2年后剖面形态的变化。图12为是否考虑季节性变化的潮滩剖面形态的对比。结果表明，在考虑季节性变化后，滩面的沉积速率介于两种波浪工况之间，与Roberts等人[29]一维模型中计算间歇性波浪的结果定性吻合。这主要是由于冬季波浪冲刷作用较强，而夏季波浪作用弱，潮滩仍以淤积为主。潮滩不再处于持续的淤积或者冲刷的单一变化，而是处于一种动态的冲淤变化中。由此看出，潮滩地貌演变对波浪的变化十分敏感。但本文中并未考虑沉积物供应以及植被等因素的季节性变化，因此对潮滩剖面的季节性冲淤变化不做具体分析。

图11 不同边界波高的潮滩剖面形态对比图

图12 考虑季节性变化的潮滩剖面形态的对比

3）波向对地貌演变的影响

本节在已建立的模型基础上对波浪方向的影响机制进行探讨。选择NNE向浪和E向浪两种工况，边界波高均取为0.1 m，且外海沉积物供应与前文所述保持一致。由不同波向下的P1点的切应力对比（图13），在大小潮时，E向浪作用下的滩面切应力均明显小于NNE向的波浪作用。这主要是因为波浪、潮流共同作用下的床面切应力大小与波浪、水流之间的夹角有关（见公式3），即波向影响着床面切应力。

其中：τm——平均床面切应力；

τw——波浪产生的床面切应力；

φ——波浪、水流之间的夹角。

为进一步分析波向对于整个计算区域的潮滩剖面形态的影响，选择距离南边界3 km、6 km、9 km的A、B、C三个断面，对比不同波浪方向下的潮滩剖面形态，见图14。可以看到，E向浪作用下，各断面的滩面淤积量均增大。由此可见，波向的改变影响着波、流共同作用下的床面切应力的大小，从而影响潮滩上的泥沙运动过程，并最终影响着潮滩的地貌演变。因此，可以认为，波向也是潮滩演变的一个重要影响因素。

图13 不同波向下的P1点的切应力对比

图14 模拟2年后潮滩剖面形态对比

4 结语

本文以江苏中部的粉砂淤泥质潮滩为原型，建立平面二维潮流、波浪、泥沙耦合地貌演变的概化数学模型，较为全面地讨论了波高、波向对于潮滩演变的影响，探讨了开敞海域粉砂淤泥质潮滩在复杂动力条件下的演变机制，重点研究了潮滩地貌演变对于波浪作用的响应。

1）潮滩水动力、泥沙以及地貌演变对波浪过程较为敏感，波浪对于潮滩演变起着重要的作用，是维持潮滩长期稳定的重要因素，在开敞潮滩的数值模拟中应考虑波浪的作用。

2）在较小的波浪作用下，潮滩仍然能够保持淤长的趋势，与不考虑波浪作用相比，整体淤积量明显减小，剖面形态仍保持双凸形；随着波浪作用的增强，床面冲刷作用增强，淤长型潮滩转变为侵蚀性潮滩，双凸形剖面逐步转变为凹形剖面，并在离岸数公里处形成陡坎，边界处有效波高增大，冲刷作用增强，陡坎位置向陆方向移动。

3）波向的变化影响着波、流共同作用下的床面切应力的大小，从而影响潮滩上的泥沙运动过程，并最终影响着潮滩的地貌演变。