徐孟飘，东培华，马骏，罗锋，张长宽，范代读，周曾,,5*

(1.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；2.华设设计集团股份有限公司，江苏 南京 210014；3.杭州市港航管理服务中心航道管理处，浙江 杭州 310016；4.同济大学 海洋与地球科学学院，上海 200092；5.河海大学 江苏省海岸海洋资源开发与环境安全重点试验室，江苏 南京 210098)

1 引言

潮滩是细颗粒沉积物在潮流为主导的水动力条件下形成的宽阔而平坦的浅滩[1]。潮滩上的水动力过程及沉积物输移过程呈双向反馈，集中表现为潮滩剖面形态的演变及潮滩上多组沉积物的分选现象。不同类型、属性的沉积物颗粒在潮滩上呈现不同的空间分布特征，该现象通常被称为沉积物分选。沉积物分选现象主要包括水平向分布和垂向分层两个方面，常见的水平向分布现象例如潮滩沉积物粒径的陆向细化，垂向分层现象一般通过潮滩底床沉积物层理来描述，沉积物在水平向和垂向的分选过程一般是紧密联系的。

层理是描述沉积特征常用的术语[2]，潮汐沉积物中的潮汐层理展现出的周期韵律性可用于估算风暴潮事件、推算古地月轨道参数，有助于推算地史时期地−月系的演化历史[3]，同时有助于分析潮滩沉积历史及提高对海平面变化规律的认知等[3]，因此开展粉砂淤泥质潮滩上多组分泥沙的垂向分选机制的研究具有重要的实际意义。潮滩上潮汐的周期性作用（涨落潮周期、大小潮周期、季节性及更长时间尺度的周期）在沉积地貌上表现出的冲刷与淤积的交替促使了互层层理（砂质层和泥质层构成）的形成[4]。对于互层层理的形成机制及内在机理的认知仍未形成统一认识，许多学者从现场实际观测入手分析其机制并取得一些进展。Deloffre等[5]在对法国、英国等河口潮间带潮滩的研究中发现沉积物表现出大小潮周期性特征，王建等[6]对江苏中部淤泥质潮滩进行实地观测认为毫米级薄水平互层层理由半日潮产生，而厘米级厚砂泥互层层理则为半月天文潮（大小潮）的产物范代读等[7]和Fan等[8]对长江口潮滩进行观测，发现潮汐层偶的厚度与大小潮周期潮差旋回性有关，层偶厚度是潮差的函数，随潮差由小潮−大潮−小潮而发生由薄−厚−薄的规律性变化。相较于现场实际观测方面，数值模拟层面的研究相对较少，是有待进一步深入的方向。本文以潮流为主导作用的潮滩为例，建立了潮滩一维概化动力地貌模型，分析沉积物垂向分布规律，探究周期性潮汐条件及不同边界含沙量对沉积韵律层结构（层理组成、层理厚度等）的影响。

2 研究方法

本文使用Delft3D模型系统中的水动力模块、沉积物输移模块以及动力地貌演变模块来进行一维水动力及地貌演变的数值模拟计算。Delft3D模型系统是一套可进行水流、波浪、水质、生态、泥沙输移、床底地貌及各个过程之间相互作用数值模拟的强大开源软件包，在内河、河口、海岸、沉积物、地貌演变水质等各方面被广泛应用。

2.1 模型介绍

2.1.1 水动力模块

假设潮滩沿岸均匀，可使用一维模式对潮滩剖面进行模拟，在静压假定和Boussinesq近似两个假定的基础上，对不可压缩流体的Navier-Stokes方程进行求解。一维浅水方程包括连续性方程和动量守恒方程。

（1）基本方程

① 连续性方程：

② 动量守恒方程：

式中，x为空间水平坐标；u为垂线平均流速在x方向的分量（单位：m/s）；h为全水深，即水面到水底的距离（单位：m）；η为相对于某标准（例如：平均海平面，MSL）的水位值（单位：m）；v为平面紊动黏性系数（单位：m2/s）；g为重力加速度（单位：m/s2）；C为谢才系数（单位：m1/2/s）。

模型中设置某一时间步长中当计算单元的水深小于 0.1 m（临界值）时，处于“干”状态，该计算单元不参与计算。当该计算单元水深大于0.2 m（临界值）时，处于“湿”状态，该计算单元重新被激活参与模型计算。

（2）由潮流引起的底部切应力

计算公式为： τc= ρgu2/C2，式中，ρ为水的密度（单位：kg/m3）。

（3）波浪模块

Delft3D默认的波浪模块为SWAN（Simulating Waves Nearshore）波浪模型，由于长周期模拟需要的计算时长较长因此本文未使用SWAN模块，而是采用了基于 Young 和 Verhagen[9]以及 van Rijn[10]提出的公式开发的简单的风浪模型，该公式于2009年被证实适用于模拟浅滩上风浪产生的作用[11]，该公式将无量纲波能量（ε =g2E/UW4）及无量纲谱峰波频率（ς=与无量纲风区及无量纲水深联系起来，对已知风场（例如：风区长度、风速及水深）的波高及波周期进行估算。

式中，A1=0.493δ0.75，B1=3.13×10−3χ0.57，A2=0.331δ1.01，B2=5.215×10−4χ0.73；UW为海拔10 m处的风速；F为风区长度；TP为峰波周期；E为波能，其公式为E=ρgHS2/8（HS为有效波高）。当风条件给定时，解以上方程可估算有效波高HS和峰波周期TP。由于波浪在传播过程中波陡增大，波浪在浅滩上会破碎，衡量波浪破碎的典型破碎指标 γ为有效波高和水深的比值（γ =HS/h），在本文中，参考Roberts等[12]参数设置，设定 γ =0.5。在波浪破碎后，波高变为受水深限制，即HS=γh。当波浪破碎时，波浪引起的底部切应力达到最大值且随着水深的向陆减小而减小[13]。根据线性波理论，底部最大轨道速度为

且由波浪引起的底床切应力为

式中，k为波数；fw为由Soulsby[14]定义的波浪摩擦因子；为非黏性沉积物的中值粒径。综合底床最大切应力是引起沉积物输移的主要因素，由潮流和风浪引起的底床切应力的叠加通常大于两者的线性总和，在本模型中参考文献[14]，综合底床最大切应力的公式为

2.1.2 沉积物输运模块

模块中的沉积物分为黏性和非黏性两种，在计算沉积物的输运时，黏性部分与非黏性部分使用不同方法进行计算，即该模块忽略了沉积物混合物的临界再悬浮条件与单个种类沉积物的临界条件的差别。黏性沙和非黏性沙的物理过程（例如侵蚀、淤积和输移）均分为不同部分考虑。

黏性沉积物的输运使用带源汇项的平流方程来描述黏性沙的侵蚀与淤积过程

式中，Qmud,e与Qmud,d分别代表侵蚀量与淤积量，该冲淤量使用被广泛运用的经典的Patheniades-Krone公式来计算，公式如下：

式中，Me为冲刷系数（单位：kg/（m2·s））；ωs为沉降速度（单位：m/s）；c是垂向平均浓度（单位：kg/m3）；τcr,e和τcr,d（ 单位：Pa）分别表示黏性沉积物的临界起动切应力及临界沉降切应力。该模型中使用Winterwerp[15]提出的“连续沉积”的概念，因此黏性沉积物的临界沉降切应力 τcr,d取 较大值（默认为1 000 Pa），则上公式可近似等于Qmud,d=ωsc。

对于非黏性沉积物，使用Sousby-van Rijn公式来计算潮流与波浪共同作用下的输移量。总沉积物输移量分为推移质和悬移质两类来计算

其中，

式中，Acal是 校准系数；CD是仅有潮流产生的阻力系数；Urms为 波质点轨道流速的均方根ucr为泥沙颗粒起动的临界流速；s为沉积物的相对密度 ；D∗无量纲粒径大小

2.1.3 地貌演变模块

床面变形方程为

式中，D是沉积物的沉积通量；E是沉积物的侵蚀通量；ρs为沉积物的密度。地貌演变的计算是基于黏性沙与非黏性沙这两种沉积物在每个水动力计算步长上 的侵蚀或淤积量的更新。

2.1.4 底床分层模块

模型中将初始底床分为若干薄层，分为3大类：冲淤层、交换层以及基准层。冲淤层的厚度 δa是一个用户根据模拟需要选取的参数，在模型计算时冲淤层中的沉积物发生侵蚀或者淤积与水体中的沉积物进行交换，但厚度始终保持不变来保证模型的稳定性，交换层在计算过程中厚度可变化，但存在一个最大厚度 δu。

模型开始计算之前，不同种类沉积物（黏性沙及非黏性沙）由初始设定的比例进行均匀混合且被分为若干层，当模型开始计算时，最顶层的冲淤层会发生侵蚀或淤积，且冲淤层下的交换层数量及厚度也会随之发生变化。底床分层模型及模型中的侵蚀或淤积的计算过程见图1。当床面发生侵蚀时，冲淤层会发生冲刷且部分沉积物会被带走，当侵蚀量大于 δa厚的冲淤层时，其下的交换层会向冲淤层进行沉积物补充使侵蚀过程继续执行，在下一个计算步长开始之前交换层会对被侵蚀过的冲淤层执行补充过程，使冲淤层厚度再次回到定值 δa，此时交换层厚度减小，以此反映床面的侵蚀过程。当床面发生淤积时，冲淤层会变厚，在下一个计算步长开始之前会执行补充过程，冲淤层厚度将保持定值 δa将多余的沉积物向下传递并入交换层，若交换层达到所设定的厚度值 δu，则底床中 将会产生一层新的交换层，以此反映床面的淤积过程。

图1 底床分层模型及其计算过程Fig.1 Bed stratigraphy modeling and schematic diagram of calculation process

2.2 模型建立

本文基于我国沿海典型粉砂淤泥质潮滩概化模型对潮滩沉积物垂向分层机制开展研究，研究使用的底床分层模型由于模型设置垂向层数较多，长期沉积过程计算时间很长，且本研究主要关注的是垂向的沉积物分布特征，因此选择一维模型更为合适，在不影响研究需求的前提下大大提高了模型计算效率。

粉砂淤泥质潮滩坡度较缓，约为1‰～3‰[16]，本文基于Zhou等[17]的研究进行模型设定，采用一维模型，初始床面为高程从–14 m 至–1 m（平均海平面以下）的缓坡，且垂向分层层厚1 mm。模型使用50 m的等距矩形网格，共计280个网格，时间步长Δt=0.25 min，满足模型计算稳定性及精度要求。图2为初始地形剖面示意图，假设初始底床由50%的黏土、25%粉砂和25%的细砂组成，且这3种沉积物沿深度方向均匀混合。

图2 初始床面组成Fig.2 Initial bed composition

模型的东侧为开边界其余均为固边界，开边界水动力条件主要使用M2及S2两个分潮，M2为太阴主要半日分潮，S2为太阳主要半日分潮。为研究垂向分层结构对不同潮流条件的响应，设置了6组工况进行计算，工况设置如表1所示。

表1 水动力条件工况设置Table 1 Hydrodynamic condition setting

由于粉砂淤泥质潮滩主要由黏土、粉砂及细砂构成，因此模型同时考虑黏性沙与非黏性沙，设定3种沉积物：黏土、粉砂和细砂。粉砂介于黏性和非黏性之间，处理时有时被当做黏性砂，但也有一些学者发现由于粉砂的成分与砂相似，透水性强、黏性差，一些国内外同行在研究中也有将之处理为非黏性砂。本研究中的粉砂粒径设置为30 μm，在粉砂界定的3.9～62.5 μm范围内，颗粒相对较大，将其考虑为非黏性沙。其中，黏土设置为黏性沙且以悬移质的方式输移，而粉砂与细砂均设置为非黏性沙，既可以悬移质的方式输移也以推移质的方式输运。初始底床由2 m的黏土、1 m的粉砂和 1 m的细砂组成，3种沉积物均匀混合。经过调试，黏土沉速取为0.5 mm/s，临界起动切应力为 0.2 N/m2，冲刷率系数取为 5×10−4kg/(m2·s)，为了使黏性沙始终处于沉降状态，因此黏性沙的沉降临界切应力取较大值设定为 1 000 N/m2（默认值），黏性沉积物干密度为500 kg/m3。属于非黏性沙的细砂的中值粒径D50取为90 μm，粉砂的中值粒径D50取为30 μm，由于Delft3D软件内默认的非黏性沙计算公式不能设置粒径小于 63 μm的颗粒，因此本文调用Van Rijn[10]公式进行非黏性沙的计算，非黏性沉积物干密度为1 600 kg/m3。在东侧开边界给予泥沙含沙量边界条件：黏土为 0.006 kg/m3，粉砂为 0.004 kg/m3，细砂为0.001 kg/m3。由于本研究是概化模型，侧重于揭示一些机制，若边界含沙量设置过大，易形成边界的剧烈淤积造成模型无法长期运行，因此选取的工况边界含沙量相对较小，但不会影响所得结论。

3 结果与分析

3.1 垂向分层结构对潮流的响应

为研究潮流对潮滩底床沉积物垂向结构的影响，模型采用了单M2分潮作用及M2、S2共同作用两种工况，其模拟的沉积物分选、分层特征结果如图3所示，计算结果是10 a接近平衡态。由于边界含沙量相对较低，潮间带发育相对较慢，且表现出不完全，但这对本文侧重说明的泥沙垂向层理现象影响较小。

图3 单M2分潮作用下（a）以及大小潮作用下（b）潮滩沉积物分选、分层特征，代表点A、B、C将用作后续分析Fig.3 Results of sediment sorting and layering under only M2tide (a) and spring/neap tidal cycles (b),representative points A,B and C are used for subsequent analyses

在仅M2分潮作用下，选取潮间带上的3点A、B、C点来进行对比分析（A点距离岸线0.5 km，B点距离岸线 1 km，C 点距离岸线 1.5 km），图 4为 3个代表点仅在M2分潮作用下垂向沉积物占比图。

图4 仅 M2分潮作用下 A、B、C 3 点处垂向沉积物占比Fig.4 Percentage of vertical sediments at points A,B and C under the action of M2 tidal constituent

由图4可见，潮间带上A、B、C 3点，由于底床高程的不断抬高，同一位置的水动力条件发生改变，各类沉积物垂向含沙量也发生改变。A、B、C 3点细砂垂向含沙量占比几乎为0。A点位于高程最大处，潮间带上部黏性沉积物占主导，因此垂向沉积物含沙量中黏土占比均大于粉砂占比。B点位于黏土与粉砂分布的过渡区域，垂向上在高程较小处粉砂占比较多，但随着沉积物淤积底床抬升，动力条件减弱，在高程较小处沉积物的垂向含沙量发生变化（黏土占比大于粉砂占比），而C点处于粉砂主导区域，因此粉砂占比整体大于黏土占比。尽管3个代表点在垂向占比上有所不同，但在垂向上含沙量曲线较为平滑没有明显波动，且均未见到明显的层理现象，可能原因包括：（1）模型中底床分层层厚设置较大（1 mm）限制了单个潮周期内潮汐层偶现象的形成；（2）模型中将粉砂当作一个组分，未考虑粗、细粉砂的不同行为过程；（3）模型未考虑絮凝沉积作用和细粒沉积物在滩上沉积后暴露空气发生脱水固结作用对再侵蚀的影响等。

图5 大小潮作用下 A、B、C 3 点处垂向沉积物占比Fig.5 Percentage of vertical sediments at points A,B and C under the action of the superposition of M2and S2constituents

图5为M2和S2分潮叠加的大小潮作用下潮间带上离岸 0.5 km、1 km、1.5 km 的 A、B、C 3点处的垂向沉积物占比分布，由于潮间带上的沉积物主要由粉砂与黏土组成，因此图中只对比了粉砂与黏土的垂向占比分布情况。由图可见，相较于图4，垂向上黏土与粉砂的含量分布线呈锯齿状表现出明显的波动且具有一定的规律性，这种周期性的规律主要是由大小潮的周期旋回造成的。A点处高程较小处的锯齿状较高程大处更明显，是由于随着沉积物的不断淤积地貌不断变化使得水深在减小，潮动力相对的越来越弱，相对较粗颗粒的沉积物（粉砂）更难以被搬运至近岸处。从水平向沉积物的分布规律来看，B点位于黏土区域与粉砂区域的过渡区域，在此区域中黏土与粉砂对于潮汐能量的变化响应较好沉积物占比变化剧烈，因此锯齿状明显。而C点位于粉砂与细砂的过渡带，该点垂向上高程较小处细砂含量相对较大，因此由黏土与细砂形成的层理结构表现较差。

为研究沉积物垂向分布与潮汐周期的对应关系，选取离岸0.6 km、高程约为–1.4 m处的P1点以及离岸0.2 km、高程约为–1.1 m处的P2点进行分析（即模型的一个网格50 m），图6为黏土与粉砂的垂向分布情况，图中的沉积物于29 d内形成，包含两个大小潮周期，该段时间内的潮位过程、P2点累计层数及含沙量变化如图7所示。

图6a代表黏土的占比分布，颜色越红代表黏土占比越多，越蓝代表占比越少，图6b代表粉砂的占比分布，同理越红表示粉砂越多，越蓝代表粉砂越少，由图可见，P1点与P2点均包含完整的潮汐韵律层，韵律层由位于上部的砂质层及位于下部的泥质层构成，且韵律层的个数与大小潮周期数量相一致。P2点垂向层理结构中黏土占比均大于P1点，这是由于P2点较P1点高程更大，水动力条件相对弱，粗颗粒沉积物更难输移至较高处。由图7可见，该潮位过程呈现“大潮–小潮–大潮–小潮”的变化趋势，与潮位过程相对应，垂向分层结构中沉积物颗粒均呈现“较粗–较细–较粗–较细”的变化规律。同时潮汐大小潮周期性的变化规律也体现于层理厚度的旋回性变化之中，由于大潮期间潮汐能量大，沉积物在较强的动力条件下更易起动并发生输运，潮差大时层理较厚而潮差小时能量衰减，因此层理较薄。

图6 P1 点（a，b）和 P2 点（c，d）两个大小潮周期内垂向沉积物占比Fig.6 Percentage of vertical sediments in two spring-neap tidal cycles at poiots P1 and P2

图7 潮位过程（a），P2 点累计层数（b）和含沙量变化（c）Fig.7 Tidal level (a)， cumulative layers (b) and changes of percentage of vertical sediments (c)

3.2 边界含沙量对垂向分层结构的影响

现场观测表明，水体中悬浮沉积物的浓度也会影响大小潮周期层偶结构的组成。因此，为研究水体含沙量对于沉积层偶结构特征的影响，在相同水动力条件的基础上设置3种边界含沙量工况，如表2所示。3种工况中，工况三的含沙量量级与通常认为的相对接近，工况一及工况二的设置实际上是对层理结构中“粗–细”形成机制的对比试验及研究，工况二及工况三的设置实际上是对层理结构中“厚–薄”形成机制的对比试验及研究。此外，由于本研究是概化模型侧重于揭示一些机制，在边界含沙量设置过大时，容易形成边界的剧烈淤积造成模型无法长期运行，因此选取的工况边界含沙量相对较小，但不会影响所得结论。

表2 3种工况的边界含沙量设置Table 2 Three settings of sediment concentration

图8为不同工况下潮间带区域（0～4 km）黏土与粉砂占比分布对比，黏土占比图中颜色越红代表黏土占比越多，颜色越蓝代表黏土占比越少，粉砂占比图中颜色含义同理。在工况一情况下，黏土边界含沙量略大于粉砂含沙量，潮汐韵律层结构于离岸0～2 km处较为明显，主要是由黏土与粉砂组成，而2～4 km处主要是由粉砂与细砂组成的层理。在粉砂边界含沙量大于黏土边界含沙量的工况二的情况中，可见在0.5～4 km范围内粉砂占比明显增大，在离岸0～2 km处虽然仍是黏土与粉砂组成的层理结构，但不论是大潮期形成的砂质层还是小潮期形成的泥质层层理中的沉积物颗粒一定程度上均发生粗化，而在离岸2～4 km范围内由粉砂与细砂构成的韵律层中的沉积物颗粒在一定程度上均发生了细化。工况三中由于黏土与细砂的含沙量是工况二中的10倍，由图8明显可见潮间带区域内整体淤积且最高处达到1 m左右，且离岸2～4 km范围内黏土与粉砂占比均有所增大，其中粉砂占比增大更明显，该范围内变为黏土、粉砂与细砂3种沉积物构成的层理结构，韵律层内的沉积物颗粒发生了明显的细化。并且，当边界含沙量浓度增大时，整个潮间带区域内大小潮对应的层理更加完整。

图8 不同工况下潮间带区域黏土与粉砂占比分布Fig.8 Distribution of clay and silt fractions in intertidal flat under different sediment concentration conditions

图9为代表点（离岸1 km处）不同工况下两个大小潮周期内黏土与粉砂的垂向占比分布，由图可见，当粉砂边界含沙量大于黏土边界含沙量时，潮汐韵律层中的沉积物均发生粗化。而当黏土与粉砂边界含沙量整体增大时（工况三）该代表点处两个潮周期内产生的潮汐韵律层明显增厚，但整体沉积物颗粒粒径变小，这主要是由于粉砂的边界含沙量虽然较大，但水体的挟沙能力有限，很大一部分粉砂较难随潮流输移至高滩处而是停留在高程较小处，而较大部分黏土却能输移至高滩处并沉积，因此粉砂占比反而减小。

图9 代表点（离岸1 km处）不同工况下两个大小潮周期内沉积物垂向分布Fig.9 Vertical distribution of sediments in two tidalcycles under different conditions at representative points

4 结论

潮滩垂向沉积韵律层的形成是多种动力因子共同作用的结果，其中大小潮周期性变化是韵律层形成的重要原因之一，本文着重从数值模拟角度探讨了大小潮周期的影响，分析了潮滩垂向沉积物的分选与分层结构形成的机理。具体结论如下：

（1）大小潮作用下，潮间带区域可见明显的分层结构（主要由黏土与粉砂形成），黏土与粉砂含沙量曲线沿垂向表现出明显的波动呈锯齿状。

（2）结果表明一个层理结构由泥质层与砂质层构成，其中泥质层在小潮期间形成，而砂质层于大潮期间形成，同时层理的厚度也与大小潮周期相关，大潮时层理较厚而小潮时层理较薄。

（3）除了周期性的潮汐条件这个主要因素，悬浮泥沙浓度也是影响大小潮周期性层理结构的重要因子。边界含沙量中粉砂占比增大会使得大小潮周期性层理结构中泥质层与砂质层均粗化且砂质层厚度增大，当边界含沙量整体显著增大时，潮滩上的垂向潮汐韵律层会更加完整且厚度明显增大。

在本文的研究中，在计算沉积物的输运时，黏性与非黏性沙使用不同方法进行计算，忽略了泥沙混合物由于黏性沙的存在呈现出的一定的黏性，可能会使得泥沙输运通量略有增大。同时，在潮滩底床垂向沉积物观测中，除了大小潮表现出的韵律性层理结构，涨落潮期间形成的潮汐层偶（基本沉积单元）以及波浪、风暴（风暴来临的水位、风暴历时）等动力因子作用下的层理重塑或破坏也常被观测到的现象，是未来数值模拟研究的重点，进而可以更加深入地掌握潮汐层理的形成与演化机制。