沙质海岸港口工程上下游岸滩演变模拟敏感性参数研究

王宁舸，宋丽佳

（1.南京水利科学研究院，江苏 南京 210029；2.中国船舶重工集团公司第七二四研究所，江苏 南京 211153）

0 引言

在沙质海岸建设港口工程时，挡沙堤工程会引起上游岸滩淤积，需准确判断上游淤积防护年限；下游岸滩冲刷后退，需采取防护措施。自20世纪70—80年代服务于毛里塔尼亚友谊港以来，我国已积累了丰富的物理模型试验研究成果和试验经验[1-2]。近20年来，随着泥沙运动基本理论的深入研究和计算机运算能力的不断提高，数值模拟成为相对更常规的手段。基于过程的岸滩长历时演变数值模拟难度相对较大，计算耗时较长，数模研究大多从一线理论或岸滩短期演变角度探讨泥沙运动的基本规律[3]。解鸣晓等[4]采用一线模型对西非塞拉利昂码头工程所在海岸的沿岸输沙能力进行了模拟。张海文等[5]对友谊港北侧沉船附近岸滩模拟了60 d的变化。Tang等[6]研究了离岸堤与直型丁坝对岸滩短历时冲淤的影响。也有少数沙质海岸岸滩长历时演变模拟的成果，如Stive等[7]介绍了人工养滩新方式“补沙引擎”，预测了沙源投放后20 a的岸滩演变趋势。

沙质海岸岸滩长历时演变模拟仍有许多技术环节尚未明确。本项研究通过数值模拟，以期探讨部分技术环节的处理方式，为沙质海岸港口工程相关数值模拟研究工作提供一定的参考。

1 研究区域概况

本项研究以毛里塔尼亚某港口工程（以下简称“A港”）为例，工程位于友谊港以南175 km的沙质海岸，岸线大致呈12°N—192°N走向，港口工程地理位置示意见图1。A港口工程建设挖入式港池，采用防浪挡沙堤掩护进出航道。防浪挡沙堤堤长初步设计为630 m，对上游来沙进行挡沙防护。下游采用150 m直型丁坝，对下游岸滩进行冲刷防护。A港初步设计方案见图2。

图1 毛里塔尼亚A港工程位置Fig.1 Location of Port A in Mauritania

A港海域面向开敞大西洋，沿岸波浪较强，外海大浪方向一般出现在WNW—NNW向范围内，这3个方向的波浪出现频率合计约为87%。港区海域波高Hs在1.0～2.0 m之间的频率为74%，2.0～3.0 m的波高出现频率为11%。经外海波浪推算以及波能加权平均统计方法计算等，判断A港代表波方向与岸线法线方向夹角为19°，H1/10代表波高为1.26 m，代表波周期为10.88 s。A港海域潮汐潮流弱，多年平均潮差不足1.0 m，潮流平均流速不足0.1 m/s[8]。

A港海域近岸泥沙中值粒径为0.23 mm，泥沙运动以波浪作用下的沿岸输沙为主要特征，沿岸输沙率约100万m3/a，方向总体为自北向南。图2还给出了工程前水深等深线图，岸线及水下地形较平直，近岸地形较陡，平均坡度约1/50，泥沙运动上界高度和下界水深分别约为+3 m和-12 m[8]。

图2 A港工程平面布置与工程前水下地形等深线示意图Fig.2 Layout of Port A and isobath of underwater terrain before the project

2 沙质海岸岸滩演变模拟敏感性分析

2.1 模型概况

本项研究采用Delft3D数值模拟软件[9]建立了平面二维沿岸输沙与岸滩演变数学模型。水流模块基于沿水深平均的浅水方程。波浪模块采用SWAN模式，基于动谱平衡方程。泥沙计算采用无黏性沙模块，基于Van Rijn（1993）公式。

数学模型采用矩形网格，模型范围覆盖全部A港港区、港区上下游各10 km及离岸6 km。离岸边界附近网格步长为50 m，近港区逐步加密至10 m。模型地形采用现场实测地形。

2.2 岸滩演变模拟敏感性参数分析

本项数学模型采用2016年现场实测潮位和流速资料进行了验证，并耦合了代表波的传播变形计算，可模拟近岸破波沿岸流。模型中泥沙中值粒径设置为0.23 mm，通过对各项输沙量计算系数的率定，最终控制模型上游来沙大致稳定在100万m3/a。

2.2.1 地貌加速因子敏感性分析

地貌加速因子（Morphological scale factor）是一个无量纲的系数，通过对每个时间步长床面的冲淤量进行该系数的倍乘，起到加速地形变化、减少模型计算时间的作用，在较长时间序列岸滩演变模拟中尤为重要。地貌加速因子取值太小会影响模型计算效率，但取值过大会导致岸滩变化对水动力的急剧影响，从而影响模型计算精度和稳定性。本项研究地貌加速因子考虑4、7和10三种情况。

不同加速因子条件下，地形等深线3 a计算结果对比见图3。图中给出了工程前岸线（+2 m）位置，挡沙堤上游岸线和水下等深线前进，表明岸滩淤积，下游反之表明岸滩冲刷。这3种地貌加速因子条件下的岸滩冲淤演变趋势一致，冲淤幅度和速率比较接近，大致在3 a后挡沙堤堤头开始出现绕沙。具体差别来看，对于上游淤积，加速因子采用10时岸滩淤积略快，采用7时岸滩淤积速率居中，采用4时岸滩淤积较前两者略慢，而下游岸滩冲刷结果与上游淤积呈相反趋势。因此，随着地貌加速因子的增大，岸滩演变计算结果呈上游淤积越快、下游冲刷越慢的趋势。

图3 不同加速因子地形等深线计算结果对比Fig.3 Comparison of terrain isobath with different acceleration factors

为了说明计算结果的合理性，采用JTS 145—2015《港口与航道水文规范》推荐的沙质海岸突堤式建筑物上游岸线演变预报计算方法，估算A港挡沙堤完全有效拦沙时间。挡沙堤堤长630 m，垂直于岸线方向的有效长度约570 m，减去上游岸滩淤积后的剖面宽度约375 m（变形高度15×淤积剖面坡比25），则挡沙堤有效拦沙长度约195 m。再根据波浪入射角19°和岸滩淤积角约13°，得到A港挡沙堤完全有效拦沙时间约2.9 a，数模计算结果（接近3 a）与公式预报结果接近。

从计算稳定性和精度看，这3种因子取值均合适，并未出现因取值变化而导致计算结果较大差异。从工程安全和实际应用角度，如何选取地貌加速因子值得进一步讨论。对于下游岸滩演变，因其还涉及到波浪绕射、下游反向输沙等复杂水沙问题，数学模型重在反映趋势，工程实际应用还需要结合物理模型进行综合分析。而对于上游岸滩演变，当发生较大淤积后，堤头绕沙将逐步发生，并产生港池航道回淤问题。因此，从工程安全角度，数模宜重点考虑上游淤积问题，在保证计算精度的情况下，加速因子可取值偏大一些。本项数模研究也对加速因子20进行了计算，但计算过程中出现发散而报错，因此地貌加速因子选择10对于工程实际应用而言是相对合适的。

综合上述分析，以A港工程为例，沙质海岸岸滩演变模拟中地貌加速因子可初步选择10。

2.2.2 模型边界条件敏感性分析

为了防止边界条件对模型计算结果的影响，需根据工况选择合适的边界条件或边界远离工程区域。这涉及到两个问题，其一是应该选择什么样的边界条件，以保证边界的稳定并尽可能符合现场实际情况，这是本小节讨论的重点；其二是边界应该远离工程多远，既保证工程影响范围不会触及边界，又使得模型范围不会太大以致降低模型计算效率，这将在下一小节重点讨论。

采用地貌加速因子10，探讨不同边界条件设置对沙质海岸岸滩演变模拟的影响。数学模型常用边界条件如潮位-流速组合边界条件，除此之外，Delft3D数值模拟软件还提供了纽曼边界（Neumann boundary conditions），可与潮位边界组合。

纽曼边界设置示意见图4中的AA′和BB′，位于垂直于岸线的两侧边界，设置为水位梯度条件，在风场和波浪场较为稳定、沿岸潮波变化不大的情况下，水位梯度可概化为0。离岸方向开边界设置潮位，见图4中的AB边界。纽曼边界可用以专门解决相对平直海岸的横向流速、潮位分布问题，防止边界扰动的发展。沙质海岸多呈顺直或顺直微弯走向，根据上述介绍，纽曼边界类型可适用。

图4 纽曼边界设置示意Fig.4 Neumann boundary setting

分别采用潮位-流速组合和潮位-纽曼组合边界条件，进行A港工程实施后沿岸输沙模拟，并统计沿岸输沙率分布，结果见表1、表2。表明两种边界条件下，距离工程较远海域的沿岸输沙率均能达到100万m3/a，并呈现上游沿岸输沙受阻挡而逐步减小、下游沿岸输沙逐渐恢复的趋势。

表1 防浪挡沙堤上游不同边界条件模型沿岸输沙率分布Table 1 Longshore sediment transport rate in the upstream of breakwater under different boundary condition models 万m3/a

表2 防护丁坝下游不同边界条件模型沿岸输沙率分布Table 2 Longshore sediment transport rate in the downstream of groyne under different boundary condition models 万m3/a

两种边界条件的沿岸输沙率分布差异主要集中在模型上下游边界处，或者说是100万m3/a沿岸输沙率的实现方式不同。潮位-流速组合边界条件下，模型上下游边界的沿岸输沙梯度很大，上游100万m3/a输沙率是通过边界附近泥沙沿程迅速补充得到的，下游则表现为向边界附近沿岸输沙急剧衰减，这显然并不合理。尽管在工程主要影响范围内，潮位-流速组合边界条件的沿岸输沙分布尚合理，但在模型边界附近体现不出现场沿岸输沙平衡状态。相比之下，潮位-纽曼组合边界条件的处理结果明显改善，边界附近沿岸输沙大致平衡，能够反映上游平衡来沙和下游输沙平衡过境的状态。

两种边界条件对边界附近沿岸输沙梯度的影响最终体现在地形冲淤上，计算3 a后水下地形等深线对比见图5。潮位-流速组合边界条件下，由于边界沿岸输沙不平衡，导致上、下游边界水下地形等深线分别呈冲刷后退和淤积前进的趋势，岸滩不稳定。相比之下，潮位-纽曼组合边界条件下边界附近等深线总体冲淤稳定并保持平顺，基本体现了自然岸滩的稳定状态，更加符合实际。

图5 不同边界条件地形计算结果对比Fig.5 Comparison of terrain isobath with different boundary conditions

据此可进一步推测，当模拟更长历时的岸滩冲淤变化时，潮位-流速边界条件引起的边界不稳定会影响模型计算的准确性，如上游沿岸输沙需通过边界岸滩的持续冲刷进行泥沙供给，但这种供给量是有限的，随着时间的推移，当边界岸滩冲刷殆尽导致泥沙供给不足后，沙源需求及相应的冲刷范围将会向工程区延伸，最终影响工程区附近的计算精度。

综合上述分析，对于自然状态下沿岸输沙基本平衡、岸滩基本冲淤稳定的沙质海岸，其岸滩演变模拟需合理选择边界条件。采用Delft3D数值模拟软件计算时，宜采用潮位-纽曼组合边界条件。

2.2.3 模型计算范围敏感性分析

港口工程对沙质海岸上下游沿岸输沙和岸滩冲淤具有一定的影响范围，因此岸滩演变模拟需要考虑合适的计算范围。根据以上敏感性分析结果，模型地貌加速因子采用10，边界采用潮位-纽曼组合边界条件，分别计算模型上、下游各5 km和10 km范围对沙质海岸岸滩演变模拟的影响。

两种模型计算范围下，港口工程上下游沿岸输沙分布见表3、表4。从中可见，在港口工程上下游各10 km计算范围下，在上下游各6 km之外，沿岸输沙大致呈不受工程影响的自然平衡状态。自6 km位置起向港口方向，沿岸输沙开始逐步减弱，表明港口工程影响范围大致达到上下游各6 km，模型10 km计算范围是合适的。相比之下，在港口工程上下游各5 km计算范围下，尽管模型边界能够大致达到100万m3/a的沿岸输沙强度，但由于工程实际影响范围已经达到上下游各6 km，因此模型范围采用上下游各5 km是偏小的，边界已经受到工程影响。

表3 防浪挡沙堤上游不同计算范围模型沿岸输沙率分布Table 3 Longshore sediment transport rate in the upstream of breakwater under different calculation range models 万m3/a

表4 防护丁坝下游不同计算范围模型沿岸输沙率分布Table 4 Longshore sediment transport rate in the downstream of groyne under different calculation range models 万m3/a

两种模型范围下岸滩冲淤演变6 a的计算结果对比见图6。从中可以看出，由于港口工程对沿岸输沙影响范围大致在上下游各6 km，因此在边界处输沙量大致均在100万m3/a的情况下，模型上下游各5 km计算范围导致泥沙淤积和冲刷的空间范围由实际的6 km被强行压缩至5 km，从而导致上游岸滩淤积前进、下游岸滩冲刷后退的速率偏快。如图6显示，模型上下游各5 km范围条件下，上游水下地形等深线淤积前进幅度略大一些，并且上游泥沙更早越过堤头并淤积在航道内，下游等深线冲刷后退幅度也相对大一些。

图6 不同计算范围地形计算结果对比Fig.6 Comparison of terrain isobath with different model ranges

综合上述分析，沙质海岸岸滩演变模拟需考虑工程对沿岸输沙和岸滩冲淤的影响范围。根据A港海域计算分析，计算范围可初步采用上下游各10 km。当然，范围超过10 km也是可行的，但考虑到模型的计算效率，10 km已基本满足需要。

3 友谊港上下游岸滩演变验证

A港海域实测地形资料相对缺乏，为了进一步证实以上敏感性分析结论的可靠性，建立友谊港海域岸滩演变数学模型，采用友谊港实测地形对模型进行验证。

3.1 友谊港概况

友谊港海域的泥沙及水动力环境与A港相似，港区近岸泥沙中值粒径为0.25 mm，岸滩坡度为1/30，近岸自北向南沿岸输沙量约100万m3/a。友谊港-10 m水深的H1/10代表波高为1.10 m，波周期为10.8 s，波浪方向与岸线法线方向的夹角为26°。

3.2 岸滩演变验证

3.2.1 模型设置

在A港海域模型的基础上，将港口工程、水深地形、波浪要素、泥沙粒径等参数替换为友谊港海域的相关参数，地貌加速因子采用10，模型边界采用潮位-纽曼组合边界条件，模型计算范围为港区上下游各10 km。

3.2.2 岸滩冲淤验证

模型以1986年12月友谊港上下游实测地形为初始地形，采用1990年7月实测地形作为验证地形。1990年7月友谊港上下游水下地形验证结果见图7，上游最大淤积点R1和下游最大冲刷点L2的岸滩剖面计算结果与实测对比见图8。结果表明，模型上下游水下地形等深线分布的计算结果与实测接近，上游最大淤积点和下游最大冲刷点的剖面变化也与实测吻合良好。

图7 1990年7月友谊港上下游岸滩地形验证结果Fig.7 Validation of upstream and downstream morphology of Friendship Port in July 1990

图8 岸滩剖面计算与实测对比Fig.8 Comparison between calculation and field measurement beach profile

综上所述，友谊港海域模型基本复演了1986—1990年友谊港上下游岸滩冲淤变化，采用经敏感性分析确定的各项设置条件是大体合适的，进行沙质海岸岸滩演变模拟时可初步考虑这种设置方式。

4 结语

通过敏感性分析，探讨了沙质海岸岸滩演变模拟中几个重要的数学模型设置条件。研究认为，综合考虑模型计算精度与计算效率，地貌加速因子可以取10。采用Delft3D数值模拟软件进行计算时，宜采用潮位-纽曼组合边界条件。综合考虑港口工程对沿岸输沙、岸滩冲淤的影响范围以及模型计算效率，计算范围可以采用港口工程上下游各10 km。

沙质海岸岸滩演变模拟的模型影响参数较多，未来还可探讨泥沙扩散系数、泥沙输运公式等设置条件，以及考虑垂向分层的三维模型与平面二维数学模型计算结果的差异。从动力角度，还可对潮位变化、波要素分级设置等进行研究。通过以上研究工作，以期进一步增强对沙质海岸岸滩演变模拟的认识，提高精细化模拟水平。