解鸣晓，李 姗，张 弛，李 鑫，阳志文

（1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，南京210029；2.交通运输部天津水运工程科学研究所工程泥沙交通行业重点实验室，天津300456；3.交通运输部天津水运工程科学研究所港口水工建筑技术国家工程实验室，天津300456；4.中交天津港湾工程研究院有限公司中国交建海岸水动力重点实验室，天津300222；5.河海大学港口海岸与近海工程学院，南京210098）

沙质海岸破波带内底部离岸流及沙坝迁移数值模拟研究

解鸣晓1，2，3，李 姗4，张 弛5，李 鑫1，2，阳志文2，3

（1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，南京210029；2.交通运输部天津水运工程科学研究所工程泥沙交通行业重点实验室，天津300456；3.交通运输部天津水运工程科学研究所港口水工建筑技术国家工程实验室，天津300456；4.中交天津港湾工程研究院有限公司中国交建海岸水动力重点实验室，天津300222；5.河海大学港口海岸与近海工程学院，南京210098）

建立了破波带内三维水动力泥沙数学模型体系，模型中同时考虑了波浪、波生底部离岸流、非粘性泥沙运动及地貌演变等动力机制。通过波浪水槽实验实测水、沙及地形数据验证，所建模式可有效模拟破波带内底部离岸流、泥沙运动过程及离岸沙坝的随时间迁移规律。模拟成果表明波浪破碎后引起的波生流系在破波带内形成一个时均垂向环流结构，其中底部时均流向指向外海，表层则指向近岸，靠近破波点处，底部离岸流的流速逐渐降低。在时均垂向环流体系的输送下，破波带内被波浪掀起的床面泥沙在底部离岸流作用下持续向外海输送，并在破波点附近堆积形成沙坝。沙坝的形成改变了局部波流场结构，在波浪的持续作用下逐渐向外海运移直至剖面达到基本平衡。

破波带；三维水沙数学模型；底部离岸流；泥沙运动；沙坝迁移

近岸沙滩位于波浪破碎带内，泥沙的沿岸及横向运动直接控制了沙滩地貌格局和岸线的稳定性。近10年来，我国风暴潮灾害频率增加，70％以上的沙质岸线遭遇明显侵蚀。目前，数值模拟技术已广泛应用于海岸工程设计中。为进一步实现合理的海岸防护，发展合适的数学模型体系预测岸滩剖面演变十分重要，是近年来国内外较为热点的研究课题。沙滩地貌的主要塑造动力为近岸波浪，床面泥沙被波浪掀起后，可在沿岸流、底部离岸流、裂流等波生流系的作用下运移。对沙滩剖面演变和横向输沙而言，破波带内的垂向波生流结构，即Stokes漂流、波生时均垂向回流和破波水滚均可起重要作用，动力机制复杂。大量现场经验和水槽实验数据表明，在破波带内流系的输送下，破波点附近发育平行岸线的沙坝，并存在一定的迁移性。沙坝的迁移改变了剖面地貌，进而反馈至波浪传播和波生流系的运移，是一个动态平衡过程。

在对破波带内水沙动力机制的数值模式开发中，历史研究多数采用平面二维模型［1-3］，亦有学者发展了三维模型［4-5］。以往研究在波生流的模拟中多数采用辐射应力理论，然而辐射应力为平面二维概念，不能刻画时均剩余动量的垂向分布，以致无法模拟底部离岸流现象。近年来，不同学者亦基于不同理论推导了辐射应力的垂向分布公式［6-10］。基于最新研究进展，Warner等人［11］，Wu和Zhang［12］等人分别将不同表达式嵌入三维水动力模型中，但未嵌入泥沙运动模块。对于沙坝运动的模拟，尹晶［13］和张弛等人［14］分别基于高阶Boussinesq方程和边界层理论建立了一维模型体系。

笔者［15-16］建立了一个基于时均动力过程的波生流运动三维数值模式，模型中考虑了辐射应力的垂向分布、破波水滚和波浪紊动效应。本文中基于笔者所建的水动力数值模式，通过增加非粘性泥沙运动和地貌演变模块，将其拓展为一个三维水沙动力地貌数值模式，并采用实验室内沙坝变迁实测数据对模型的有效性进行检验，并进行相关讨论。

1三维水沙数值模拟建立

1.1水动力模式

水动力方程的基本形式为雷诺方程，见式（1）～（4）。在对波生流的模拟中，引入波浪时均剩余动量的垂向分布、破波水滚和波浪紊动作为驱动力。模型网格剖分中采用水平向矩形、垂向贴体坐标的形式。

式中：t为时间；x和y为水平向坐标；η为水位；U和V分别为x和y方向的水平流速分量；ω为σ坐标系下垂向流速；D为水深；M为时均波生剩余动量，可由Lin⁃Zhang公式求解［7］，表达式见式（5），R为破波水滚动量；KMc和AMc分别为波流共同作用下的垂向和水平紊动系数；ρw为海水密度。

式中：E为波能；n为波能传递率；k为波数；T为波周期；δ为Kronecker张量标记；i和j分别代表x和y方向。波流共同作用下的底部剪切力τcw采用Soulsby等人［17］的理论。

在波流共同作用的紊动描述中，分别单独求解水流与波浪引起的紊动系数AM与KM，并将其叠加［15-16］，水流引起的水平紊动系数AM（σ）采用Smagorinsky方程求解，形式见式（6），其中Δx和Δy为x和y方向的空间步长；Cs为经验系数，取值在0.1～0.2左右。波浪引起的水平紊动系数采用笔者［15-16］推导的公式，见式（7）中所示，其中λ为无因次系数，取值在0.2～0.5左右［15］。水流引起的垂向紊动系数KM（σ）采用Mellor⁃Yamada紊流闭合方程求解，波浪引起的垂向紊动系数表达式见式（8），其中b为无因次系数，取值在0.001左右［15］。

笔者［15］推导了波浪破碎引起的水滚能量传输方程，见式（9），式中α为水滚能量传递率，取值在0.0～1.0之间；AR为水滚面积；ρR为水滚密度；C为波速；KR＝0.375（0.3+2.4s），s为底坡坡度。当波浪参数确定后，取破波点外AR＝0，然后自破波点向岸迭代求解。水滚能量垂向分布采用Haas和Warner［18］建议形式。

方程求解技术采取内外模式分裂法，平面求解采用显式差分，垂向求解采用隐式差分。网格配置采用Arakawa-C网格，时均剩余动量和水滚动量项布置在网格中心。动边界处理采用Oey［19］提出的OGCM法。

1.2波浪模型

波浪模型采用折绕射联合REF/DIF数值模式，模型基于抛物缓坡方程理论，可模拟波浪的浅水变形、折射、绕射和破碎等现象，同时可考虑水流对波浪的反馈影响。

1.3非粘性泥沙运动模型

悬沙运动采用三维对流扩散方程，见式（10），其中C为含沙量；DH、DV分别为泥沙水平及垂向紊动扩散系数；ω为泥沙沉速，由张瑞瑾公式计算，见式（11）；D50为泥沙中值粒径；ν为水体运动粘性系数；ρs为泥沙颗粒密度。泥沙悬浮在水体中，将会改变海水密度，从而引起密度斜压效应，见式（12），其中Cvol为体积含沙量，模拟时根据每一时步所求解的水体含沙量数值修正含沙水体密度，并代回水动力方程中求解。悬浮泥沙与床面的交换形式采用“参考高度”概念，参考点含沙量采用式（13）～（15）求解［4］，其中ks为Nikuradse糙率；在悬沙与床面的交换项处理中，采用式（16）～（17），其中kmx代表垂向网格中高于参考高度a的第一层网格中心位置。在低于参考高度至床面的区间，认为泥沙以推移质形式运动，输沙率表达式见式（18）～（22），其中q=qbc+qbw为总输沙率；u*′为摩阻流速；ub为底层网格中心处的流速数值，vr为垂线平均流速；Uw为波浪底部最大质点流速。床面变形方程见式（23），其中zbed为床面高程。在式（17）中的波浪输沙率计算中，同时考虑了非线性波浪水质点运动不对称的影响［4］。

2沙坝变迁模拟及验证

为验证模型的有效性，采用Arcilla等人［20］于Delft水工所开展的LIP11D大型波浪水槽实验作为算例。实验所用水槽的长度、宽度和深度分别为233 m、5 m和7 m。在实验1 b中，坡脚前波高达1.5 m，波周期5.0 m、床面泥沙中值粒径0.22 mm，接近天然沙滩前缘的波浪环境，实验中持续17 h的窄谱不规则波作用于原始沙坝地形上，并同步测量了多个站点不同时刻的波高、增减水、时均流速、含沙量和地形变化等关键参数。图1中示意了实验地形和观测站位置信息，表1中给出了数学模型中所采用的主要参数情况。

为量化比较模拟值与实测值间的误差情况，采用2个误差指标进行分析，分别为均方差RMS和相关系数COR，其中前者反映模拟值和实测值间的偏差，后者则反映两套数值的相关度。

图2和图3中分别给出了T=1 h，T=9 h和T=17 h时刻的模拟与实测波高、波浪增减水的对比情况。基于对比结果，模型可较好地反映出不同时刻波浪沿水槽的传播规律和增减水沿程分布情况，模拟误差均很小。基于模拟结果，波浪在实验的沙坝地形上存在3次明显的破碎过程，分别发生在（1）x≈50 m处，此处波浪初传播至陡坡上方，在浅水变形作用下发生破碎；（2）x≈140 m附近，此处位于外侧沙坝顶端；（3）x≈160 m附近，此处位于内侧沙坝顶端。在内、外两个初始沙坝间的沟槽内（140 m＜x＜160 m），波高略有增大。至于增减水的沿程分布，则呈现出破波区内增水、浅水变形区内减水的整体趋势。

图1 LIP11D 1b实验地形及测站位置示意Fig.1BathymetryandobservationsectionlocationsoftheLIP11D1bcase

表1 LIP11D 1b实验数值模拟采用主要参数情况Tab.1 Input parameters in the numerical model for LIP11D experiment case 1b

图2 LIP11D 1b实验不同时刻实测与模拟波高分布对比Fig.2 Comparisons between the simulated and observed wave height for LIP11D 1b case

图3 LIP11D 1b实验不同时刻实测与模拟波浪增减水分布对比Fig.3 Comparisons between the simulated and observed wave setup for LIP11D 1b case

图4中分别给出了不同观测时刻、不同断面位置处的模拟与实测时均流速垂向分布的对比情况。经对比，在外侧沙坝顶端附近（130 m＜x＜150 m），模型计算所得时均流速数值较实测值略低。经分析认为，这一误差的原因在于所模拟床面为动床状态，特别是在波浪和时均流速的作用下外侧沙坝有明显迁移，从而导致模拟所得地形和实测地形存在一定差异，进而影响了时均流速的计算精度。尽管如此，总体来说所建模型可较为有效的反映出底部离岸流的垂向分布趋势，且随时间的变化规律亦与实测值较为接近。

从时均流速的垂向分布趋势来看，最大离岸流速出现在近床面附近，并向表层逐渐衰减；从沿程分布来看，最大底部离岸流出现在沙坝附近区域。至表层水体，离岸流逐渐反向，流向转为指向岸线。这一现象反映出在破波带内存在一个时均垂向回流结构，其表层向岸、底层离岸，是沙坝向海侧运移的主要驱动力。

图5中分别给出了不同观测时刻、不同断面位置处的模拟与实测时均含沙量垂向分布的对比情况。与图4中所反映出的规律类似，模拟结果表明在地形变化较为平缓的离岸区内（60 m＜X＜130 m），模拟值与实测值的吻合度较佳，最大误差出现在外侧沙坝附近（130 m＜X＜140）。

经分析，这一误差的成因主要有三：（1）根据图4中结果所示，在该区域的模拟所得底部离岸流速小于实测值，导致掀沙能力较实际略有不足；（2）实验中的实际地形形态存在沙纹现象，即床面并非严格平整，而这些沙纹将引起床面局部的小尺度涡旋，但由于本文中所采用的模拟理论为基于时均模型，难以反映一个亚波周期内的小尺度紊动结构；（3）在模型中采用的底部摩阻高度取为常值，而在实际水槽中的沙纹结构并非沿程均匀，导致了一定的不匹配。实际上，由于模型计算中将动力过程和地貌过程耦合，从而地形模拟的精度误差亦可反馈至水沙运动信息，亦将引起一定的模拟误差。

图4 LIP11D 1b实验不同时刻实测与模拟时均流速分布对比Fig.4 Comparisons between the simulated and observed undertow velocity for LIP11D 1b case

图5 LIP11D 1b实验不同时刻实测与模拟时均含沙量分布对比Fig.5 Comparisons between the simulated and observed suspended SSC for LIP11D 1b case

图6中给出了T=17 h时刻的模拟与实测输沙率的对比情况。总体来说，模型可较好地刻画输沙率的沿程分布规律和量级。在沙坝的外海侧，泥沙输移方向为指向近岸，而在沙坝的向岸侧，泥沙输移方向则为指向外海，最大输沙率出现在沙坝顶端的波浪强破碎区。图7中给出了T=17 h时刻的模拟与实测沙滩剖面地形的对比情况。经统计，模拟所得的沙坝高度略低于实测高度，且沙坝间的沟槽深度略大于实测深度，误差的成因与模拟所得底部离岸流强度和含沙量的误差所致。

尽管如此，综合以上对不同参数的误差分析，从整体上来看，模拟所得的底部离岸流流速、含沙量、输沙率和剖面地形与实测值均达到了较佳的相关度，在绝大多数断面处，相关系数COR均超过了80％，且均方差RMS不大。因此，认为所建模型可有效地复演破波带内沙滩剖面在强波浪作用下的水动力特征、泥沙输移规律和沙坝迁移、演变特征，模型精度较高。

需指出的是，在沙坝的横向迁移中，同时可受指向岸侧的Stokes漂流和指向海侧的底部离岸流双重影响。在本文采用的算例中，重点讨论了沙坝在强底部离岸流条件下的向海迁移。然而，在弱浪环境下，由于波浪破碎强度低，底部离岸流速较弱，而Stokes漂流占优，可能导致沙坝的向岸侧迁移，因此在后续工作中，应进一步针对这一现象进行相关研究工作。

图6 LIP11D1b实验T=17h时刻实测与模拟输沙率沿程分布对比Fig.6Comparisonsbetweenthesimulatedandobservedsediment transportrateforLIP11D1bcase(T=17h)

图7 LIP11D 1b实验T=17 h时刻实测与模拟地形沿程分布对比Fig.7 Comparisons between the simulated and observed bed elevation for LIP11D 1b case(T=17 h)

3结论

本文建立了一个基于动力过程的破波带内水动力及沙滩剖面地貌演变的三维数值模式，模型中同时考虑了波生流的垂向分布、悬沙和底沙运动以及剖面地形的耦合反馈。采用基于大型波浪水槽的实验数据对模型进行了详细验证，并评价了模型表现力。研究结果表明所建数值模式可较为有效地复演破波带内的水沙动力机制和沙坝发育和迁移规律，模拟精度较高。主要结论如下：

（1）波浪破碎后，在Stokes漂流、波生时均剩余动量和破波水滚动量传输的综合作用下，在破波点岸侧形成一个垂向时均环流结构，即底部时均流速指向外海，而表层时均流速则指向岸侧，最大底部离岸流速出现在破波点内侧。

（2）在破波带内底部离岸流的输送下，波浪起动的沙滩泥沙可以悬沙和底沙的形式持续向海侧输送。由于底部离岸流速向海侧逐渐降低，从而泥沙在破波点附近堆积形成沙坝。沙坝的发育改变了沙滩剖面地形，进一步改变了波浪传播和底部离岸流的分布格局，进而使得沙坝有向海迁移的趋势，直至动力环境和地貌形态达到基本稳定状态。

［1］Nam P T,Larson M,Hanson H,et al.A numerical model of beach morphological evolution due to waves and currents in the vicini⁃ty of coastal structures［J］.Coastal Engineering,2011,58:863-876.

［2］Nam P T,Larson M.Model of nearshore waves and wave⁃Induced currents around a detached breakwater［J］.Journal of Waterway, Port,Coastal and Ocean Engineering,2010,6:156-176.

［3］Ding Y,Wang S S Y.Development and application of coastal and estuarine morphological process modeling system［J］.Journal of Coastal Research,Special Issue,2008,52:127-140.

［4］Lesser G R,Roelvink J A,van Kester J A T M,et al.Development and validation of a three⁃dimensional morphological model［J］. Coastal Engineering,2004,51:883-915.

［5］Roelvink J A.Coastal morphodynamic evolution techniques［J］.Coastal Engineering,2006,53:277-287.

［6］Ardhuin F,Rascle N,Belibassakis K A.Explicit wave⁃averaged primitive equations using a generalized Lagrangian mean［J］. Ocean Modeling,2008,20:35-60.

［7］Lin P Z,Zhang D.The depth⁃dependent radiation stresses and their effect on coastal currents［C］//ICHD.Proceedings of the 6th In⁃ternational Conference of Hydrodynamics:Hydrodynamics VI Theory and Applications,2004:247-253.

［8］Mellor G L.Some consequences of the three dimensional current and surface wave equations［J］.Journal of Physical Oceanogra⁃phy,2005,33:1 978-1 989.

［9］McWilliams J C,Restrepo J M,Lane E M.An asymptotic theory for the interaction of waves and currents in coastal waters［J］.Jour⁃nal of Fluid Mechanics,2004,511:135-178.

［10］Xia H Y,Xia Z W,Zhu L S.Vertical variation in radiation stress and wave⁃induced current［J］.Coastal Engineering,2004,51: 309-321.

［11］Warner J C,Sherwood C R,Signell R P,et al.Development of a three⁃dimensional,regional,coupled wave,current,and sediment⁃transport model［J］.Computers&Geosciences,2008,34:1284-1306.

［12］Wu X Z,Zhang Q H.A three⁃dimensional nearshore hydrodynamic model with depth⁃dependent radiation stresses［J］.China Ocean Engineering,2009,23:291-302.

［13］尹晶.海岸沙坝运动的实验与数值模拟研究［D］.大连:大连理工大学,2012.

［14］张弛,郑金海,王义刚.波浪作用下沙坝剖面形成过程的数值模拟［J］.水科学进展,2012,23（1）:104-109. ZHANG C,ZHENG J H,WANG Y G.Numerical simulation of wave⁃induced sandbar formation［J］.Advances in Water Science, 2012,23（1）:104-109.

［15］Xie M X.Establishment,validation and discussions of a three dimensional wave⁃induced current model［J］.Ocean Modelling, 2011,38:230-243.

［16］Xie M X.Three⁃dimensional numerical modeling of the wave⁃induced rip currents under irregular bathymetry［J］.Journal of Hy⁃drodynamics,2012,6:864-872.

［17］Soulsby R L,Hamm L,Klopman G,et al.Wave⁃current interaction within and outside the bottom boundary layer［J］.Coastal En⁃gineering,1993,21:41-69.

［18］Haas K A,Warner J C.Comparing a quasi⁃3D to a full 3D nearshore circulation model:SHORECIRC and ROMS［J］.Ocean Mod⁃eling,2009,26:91-103.

［19］Oey L Y.An OCGM with movable land⁃sea boundaries［J］.Ocean Modelling,2006,13:176-195.

［20］Arcilla A S,Roevink J A,O'Connor B A.The Delta Flume'93 experiments［C］//ASCE.Proceedings of Coastal Dynamics,2010: 488-502.

Numerical modeling of the undertow and sandbar migration process in the surfzone

XIE Ming⁃xiao1,2,3,LI Shan4,ZHANG Chi5,LI Xin1,2,YANG Zhi⁃wen2,3
(1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Nanjing Hydraulic Research Institute,Nanjing 210029,China;2.Key Laboratory of Engineering Sediment,Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China;3.National Engineering Laboratory for Port Hydraulic Construction Technology,Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China; 4.Key Laboratory of Coastal Engineering Hydrodynamics of CCCC,Tianjin Port Engineering Institute,CCCC, Tianjin 300222,China;5.College of Harbor,Coastal and Offshore Engineering,Hohai University,Nanjing 210098, China)

A process⁃based three⁃dimensional numerical model for surfzone hydrodynamics and beach profile evolution was developed.The model incorporates the coupling process of waves,wave⁃induced undertow,sediment transport and morphology response.Using the established model,the wave⁃induced undertow structure and the sand⁃bar migration characteristics were modeled,and a series experimental datasets were applied to validate the model. The results show that the model can satisfactorily describe the surfzone hydrodynamic and the sandbar migration characteristics in the surfzone with acceptable precision.Once the wave breaks,the vertical distribution of the un⁃dertow indicates a phase⁃averaged circulation structure.Outside the breaking zone,the undertow rapidly diminish⁃es,and the undertow direction points to the offshore near the water surface.For the sediment transport characteris⁃tics along the beach profile,the sediments entrained by waves continuously move to the breaking point,and deposit near the breaking point to form a sandbar where the undertow velocity decreases.

surfzone;process⁃based three⁃dimensional numerical model;wave⁃induced undertow;sediment transport;sandbar evolution

TV 142；O 242.1

A

1005-8443（2016）04-0349-07

2016-01-27；

2016-02-25

国家自然科学基金青年科学基金（41306033）；交通运输部应用基础研究项目（2014329224330）；水文水资源与水利工程科学国家重点实验室开放项目（2014492211）；天津市自然科学基金青年项目（16JCQNJC06900）

解鸣晓（1982-），男，山东省青岛人，副研究员，博士，主要从事海岸动力学研究。

Biography：XIE Ming⁃xiao(1982-),male,associate professor.