海湾由于受海洋动力(主要是潮流和波浪)和河流径流的作用, 其海底地貌呈现出复杂多样的形态, 如涨潮三角洲、落潮三角洲及浅滩和深槽。随着海湾的开发利用, 海湾地貌形态演变的预报越来越重要。对这些地貌形态的研究, 以往多采用现场观测和历史资料收集。近些年, 随着计算机技术的发展进步, 采用数值模拟来研究海湾地貌形态的演变逐渐成为一种有效和实用的方法。Lanzoni等^[1]采用宽度积分的一维数值模型对海湾地形的长期演变进行了计算模拟, 通过对比收缩型海湾和矩形海湾的水底平衡剖面, 分析了海湾的收缩形态对水底平衡剖面的影响。张心凤^[2]采用数值模拟方法以2003年的地形为初始地形, 对黄茅海20 a后的地形演变进行了预报。

泻湖型海湾是一种广泛存在的海湾地貌类型, 该类型海湾的湾口两侧常常存在沙坝或沙堤, 这使得海湾仅能通过穿越沙坝(沙堤)的潮汐通道与外海相连, 从而导致湾内涨潮三角洲, 潮流深槽和湾外落潮三角洲的形成。为了研究这些海湾地貌形态的形成机理和演化过程, 本研究以水东湾地貌形态为模拟对象, 对湾中涨潮三角洲(大洲岛)和落潮三角洲的形成演化过程进行了计算模拟, 所得结果与实际地形符合, 这为进一步研究泻湖型海湾地貌的动力过程提供了参考。

1 水东湾地貌特征

水东湾位于广东省西部电白县境内, 东侧距离珠江口约300 km, 西侧距离湛江约90 km。海岸地貌属于沙坝-泻湖地貌, 如图 1所示。由图可见, 水东湾地貌体系发展全面, 包括了涨潮三角洲、潮汐通道和落潮三角洲等地貌单元^[3]。其海湾湾口较小和沿岸输沙较强的特征使得湾口东西两侧存在由沿岸流导致的延伸至湾口的沙坝, 可在低潮时露出水面。西沙坝从湾口以西至晏境岭, 总长约9 km, 宽2~3 km; 东沙坝从湾口以东至博贺, 总长约13 km, 宽2.4~4.4 km^[4]。海湾纳潮面积20~34 km²(面积随大小潮而变化)^[5], 因而有利于湾中岛等涨潮三角洲的形成。海湾通过潮汐通道与外海相通, 潮汐通道长7 km、宽500~800 m、深5~15 m。通道的两端发育有涨潮三角洲和落潮三角洲堆积体^[6], 前者即为大洲岛, 岛的南侧和东侧各有一条宽度为50~160 m、深度为2~5 m的潮流深槽, 深槽环绕湾中岛呈树状分布。落潮三角洲为湾口外形成的沿潮流方向的扇形泥沙堆积体^[4], 它的主体是一宽阔平坦的细砂堆积的平台, 面积约22 km², 半径可达5 km, 外缘水深小于2 m。在波浪和潮流的作用下, 落潮三角洲上存在沙嘴和沙波等次一级地貌。

以上水东湾地貌体系的形成涉及潮流和波浪的混合作用, 整体上以潮流作用为主。波浪作用表现为在落潮三角洲区域发生破碎而引起泥沙悬移并通过伴随波浪破碎而产生的沿岸流形成沿岸输沙。海湾湾内地形主要由潮汐动力控制, 而湾口两侧沙坝及湾外落潮三角洲除受潮流作用外还受波浪作用影响。后者的存在可以从晏境岭附近岸段的沿岸总输沙量达28×10⁴ m³看出。这些沿岸输沙成为了潮汐通道在没有来自泻湖的泥沙时的泥沙供给源, 所供给的泥沙涨潮时被带入泻湖内, 落潮时被带入落潮三角洲区域, 对三角洲泥沙沉积产生贡献。

2 计算模型及验证 2.1 计算模型

所采用的计算模型包括三个部分:水动力方程、输沙模式和地形演变方程, 三者是相互耦合的。

由于所计算的水东湾区域水深较浅, 速度沿水深方向的变化并不大, 所以可采用二维水深平均浅水方程来计算潮流运动, 方程的表达式为:

连续方程:

$\frac{{\partial \eta }}{{\partial t}} + \frac{{\partial (Hu)}}{{\partial x}} + \frac{{\partial (Hv)}}{{\partial y}} = 0$ .

(1)

x方向动量方程:

$\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}} = fv - {\rm{g}}\frac{{\partial \eta }}{{\partial x}} - {\rm{g}}\frac{{u\sqrt {{u^2} + {v^2}} }}{{C_z^2H}}$ ,

(2)

y方向动量方程:

$\frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}} = - fu - {\rm{g}}\frac{{\partial \eta }}{{\partial y}} - {\rm{g}}\frac{{v\sqrt {{u^2} + {v^2}} }}{{C_z^2H}}$ ,

(3)

式中, η为水面升高; u、v分别为流速在x, y方向的分量; H为总水深, $H = \eta + h - {z_b}$, h为静水深, z_b为水底床面高程; f为科氏力系数, 由于本文所研究区域较小, 故科氏力系数不予考虑; g为重力加速度; C_z为谢才系数, ${C_z}$$ = {H^{1/6}}/{n_r}$, n_r为曼宁系数, 曼宁系数n_r的一般取值范围为0.005~0.02, 这里的取其上限, 即n_r=0.02。

地形演变方程由水底变形质量守恒得到, 其表达式为

$(1 - {{\rm{n}}_p})\frac{{\partial {z_b}}}{{\partial t}} + \frac{{\partial {q_x}}}{{\partial x}} + \frac{{\partial {q_y}}}{{\partial y}} = 0$

(4)

其中, $\boldsymbol{q}=\left(q_{x}, q_{y}\right)$为泥沙体积输沙率; ${{\rm{n}}_p} = 0.4$为泥沙孔隙率, 沙质泥沙的孔隙率一般在0.36~0.4左右, 对于一般的天然沙质海床, 可假定泥沙孔隙率为0.4^[7]。

对输沙率q的计算, 计算模拟分别采用了推移质输沙和全沙输沙两种输沙模式, 以考虑不同输沙模式对计算结果的影响。推移质输沙率的表达式取为基于能量原理的常用的表达形式^[8]:

$\boldsymbol{q}={\rm{m}}|\boldsymbol{u}|^{2} \boldsymbol{u}$

(5)

式中, $\boldsymbol{u}=(u, v)$为速度矢量; m为系数, 取0.001 s²/m, 其包含泥沙颗粒大小等的影响。这一取值来自于文献[9], 该值的合理性可通过Meyer-Peter和Muller公式^[10]得到了解, 即在该式中忽略起动流速项(该项数值一般很小), 则该输沙率公式可以写为公式(5)的形式(水底摩擦系数用谢才系数表达, 曼宁系数n_r为0.02)。

全沙输沙率的表达式取为^[11]

$\boldsymbol{q}_{s}=\frac{0.05 \alpha|\boldsymbol{u}|^{4} \boldsymbol{u}}{\mathit{\boldsymbol{g}}^{0.5} C_{z}^{3} \Delta^{2} \mathit{\boldsymbol{d}}_{50}}$

(6)

其中, $\boldsymbol{q}=\left(q_{s x}, q_{s y}\right)$为全沙体积输沙率; α为常数, 取0.25;泥沙相对密度$\Delta = ({{\rm{ \mathsf{ ρ} }}_s} - {\rho _w})/{\rho _w}$; ρ_s=2 650 kg/m³为泥沙的密度; ρ_w为水的密度; d₅₀为泥沙中值粒径, 因水东湾潮汐通道处泥沙的中值粒径为0.10~0.15 mm, 湾内泥沙的中值粒径为0.004~0.029 mm^[12], 综合考虑取d₅₀=0.12 mm。

以上方程的数值求解的计算区域的具体大小和边界如图 1所示, 其长为19 km、宽为15 km。水动力方程的离散采用矩形交错网格, 并利用ADI格式(交替方向隐式差分格式)^[13]求解, 空间步长为Δx= Δy=100 m, 时间步长为Δt=4 s。地形演变方程采用WENO格式^[14]进行求解, 其可以考虑地形变形中可能出现的的较陡的地形变化, 求解中空间步长取为Δx=Δy=100 m, 时间步长取为Δt=8 s。

2.2 边界条件

计算的初始地形根据茂名港及其附近海区的海图确定(1︰4 000, 船讯网电子海图软件), 但湾内地形不包括湾中岛(该岛将在数值模拟过程中自动生成), 见图 2。计算的实际水深值为图 2中海图水深再加上当地海平面水深(1.94 m)。水平二维坐标系的原点取为计算域的左下角, 其地理坐标为21°25′N、110°59′E, x轴正向为正东方向, y轴为正北方向。

计算区域边界的边界条件包括:

潮流边界条件:计算区域左右两侧和下侧边界(东西南三面)设置为潮流入射边界, 即给定潮流的潮位:

$\eta = {\rm{A}}\sin \omega t$

(7)

式中, A和$\omega = 2{\rm{ \mathsf{ π} }}/T$分别为潮流速度幅值和频率, T为潮汐周期, 取为半潮周期, 因当地潮汐为不规则半日潮。数值模拟中由于缺乏各分潮潮位数据, 因此采用具有代表性的M₂分潮进行计算, 潮流幅值采用文献[6]中所采用的平均潮幅: A=0.88 m。

径流及来沙边界条件:径流和来沙条件取为零, 因为海湾内陆并没有汇入水东湾的大型河流, 仅在湾顶有四条长不足10 km的小溪流(旦场河、寨头河、水东河及南海河)汇入, 但小溪流的流量均很小, 除偶遇暴雨时有明显径流之外, 径流对水东湾的影响可以忽略^[15]。又由于四条小溪流所包含的泥沙量很少, 故水东湾内陆来沙也很少。在海湾的外海水域, 水东湾海域平时海水很清澈, 一般情况下平均含沙量为0.008 4 kg/m³, 这导致进入海湾所在水域的泥沙也较少, 不考虑海域供沙对计算结果影响不大, 所以在外海边界处不考虑来沙的输入输出。

固壁边界:海湾四周边界和潮汐通道边界以及湾口两侧海岸的边界均设为固壁边界。

干湿边界:湾中岛在形成的过程中边界是不断变化的, 所以其边界处理为干湿边界, 这样可以考虑在湾中岛形成过程中泥沙输移导致的湾中岛四周边界不断淤积、水深不断变小、水底出水和出现干湿边界转化的情况。干湿边界条件可以自动识别变化的水域边界, 从而模拟对应的边界变化。

波浪条件:本研究没有考虑波浪对区域的影响, 仅考虑了潮流对水东湾地形演化过程的作用。这是考虑到水东湾口门区域较为狭窄, 潮汐通道内侧的泻湖湾主要受到潮汐动力作用的影响。波浪影响主要体现在湾口门外侧的落潮三角洲:湾口两侧波浪破碎产生的沿岸流及沿岸输沙会使泥沙在湾口处聚集, 从而使落潮三角洲的淤积得到加强。由于本研究的重点是湾内落潮三角洲的形成和演化, 所以不考虑波浪影响不会对湾内地形演变产生显著影响, 而仅会对湾外落潮三角洲的结果有一定影响。

2.3 模型验证

为了验证潮流计算模型, 将图 2中水东湾潮流通道入口处H点潮位和流速的模拟结果与H点实际潮位和流速值^[16]进行了对比, 如图 3所示。由于目前还不存在湾内的潮流实测结果, 所以验证仅针对上述潮汐通道处的潮流。这一验证事实上也可以保证湾内的潮流流场具有可靠的精度, 因为在控制方程是适用的条件下, 湾内计算结果的精度是依赖于边界入流(潮汐通道内潮流)的精度的。图 3中实测数据取自秦福寿等^[16]于2012年9月现场测量的实测数据, 由图可见, 潮位和流速计算值与实测值基本吻合, 所以可以认为本文模型具有可靠的精度。

为了考察整个海湾内的潮流结果, 图 4给出了水东湾涨、落急流速分布图。由图可见, 湾内潮流运动呈现典型的往复流态, 潮汐通道口受两岸束流的影响, 通道及泻湖湾内沟槽处流速较高, 且以落潮流占优势。这与实际水东湾潮流流态是一致的。根据2008年实测数据显示, 潮汐通道内的涨潮、落潮的最大流速分别为0.6 m/s和0.9 m/s^[6], 图 4中计算的最大涨潮流速为0.62 m/s, 最大落潮流速为0.87 m/s, 与上述实际的潮流流速最大值和最小值一致。

3 水东湾计算模拟结果

由上述图 4中结果也可见, 潮流进入湾内后涨潮流随着湾内水域的扩大而流速逐渐变小, 在湾内出现分叉, 使得一部分潮流沿潮汐通道由西北向西偏转, 然后流向湾顶; 另一部分潮流沿潮汐通道绕开大洲岛后再向西北偏转流向湾顶, 这一流动特征导致了分叉之处泥沙随着水流速度的减小而产生淤积, 导致湾中岛(大洲岛)的形成。本节给出采用方程(1)~(6)所得到的该湾中岛形成和演化的数值模拟结果。为了考虑不同输沙成分(推移质和悬移质)的影响, 模拟分别采用了全沙输沙(6)和推移质输沙(5)两种不同的输沙模式, 通过分析这两种输沙模式对应的地形演变可以看出不同输沙成分的贡献。另外, 计算的初始地形取为将湾中岛去掉后的水东湾地形, 以便研究湾中岛由无到有的成长过程, 见图 5a。

3.1 全沙输沙模式的结果

图 5给出了采用全沙输沙公式(6)所得到的5 a、15 a和20 a演变时间的水东湾地形结果, 对应的湾中岛的面积和增长率由表 1列出。由图可见, 5 a的模拟结果对应湾中岛初始的形态, 该形状为两个独立的小岛。5~15 a为湾中岛成长发展阶段, 两个小岛的面积逐渐相连而形成一个较大面积的湾中岛, 由表 1可见其增长率最大可达0.066 km²/a; 15~20 a是岛的发展成熟阶段, 增长率变小(0.054 km²/a, 见表 1), 到20 a时刻已经达到了稳定期(增长率接近于零), 最终形成了一个长1.7 km、宽0.5 km的条形湾中岛。

为了将以上结果与实际的湾内地形进行对比, 图 6给出了水东湾实际地形图。将两者对比可见, 图 5中20 a湾中岛的位置与实际的大洲岛的位置基本一致, 长度也与实际湾中岛的基本一致(前者为1.7 km, 后者为1.4 km), 但宽度比实际的要小(前者为0.5 km, 后者为1 km), 导致整个岛的面积比实际的要小。

3.2 推移质输沙模式的结果

图 7给出了采用推移质输沙模式(5)所得到的5 a、15 a和25 a演变时间的水东湾地形结果, 对应的湾中岛的面积和增长率由表 1列出。由图可见, 该输沙模式所得到的湾中岛的总的发展进程与上述全沙模式的类似, 但在具体细节上存在差别。为了说明这一差别, 图 8给出了岛面积的局部图, 该图为两种输沙模式所对应的不同年份的湾中岛的0 m等深线图。同时由表 1给出了对应的岛的形状、位置、大小及增长率随模拟时间的变化。由该图和表可以看出: ①形状方面。全沙模式得到的岛的形状为西北向东南伸展的条形(10 a后), 而推移质模式得到的形状为近似的三角形或矩形(10 a后); ②位置方面。全沙模式得到的湾中岛的中心位于x=5.45 km, 而推移质模式得到的湾中岛中心位于x=5.92 km, 二者相差约0.47 km; ③大小方面。推移质输沙模式得到的岛的面积比全沙输沙模式得到的要小, 5 a、10 a、15 a、20 a、25 a岛的面积分别为后者的60.0%、44.2%、42.9%、54.2%、89.4%; ④增长率方面。除25 a外, 推移质输沙模式的岛的增长率要比全沙模式的小: 10 a、15 a、20 a岛的增长率分别为后者的39.4%、38.5%、77.8%。这导致该输沙模式的湾中岛达到稳定形态所需要的时间比全沙输沙模式的要长(前者为25 a, 后者为20 a)。造成以上差别的原因是全沙模式同时考虑了推移质和悬移质, 而推移质模式仅仅考虑了推移质输沙。由于悬移质输沙的加入使得全沙模式对应的湾中岛的发展速度变快, 所以全沙模式前期所得到的岛的面积及增长率要大, 且能更快的达到稳定时刻, 另外也使得输沙的空间分布出现差别, 即导致了全沙输沙模式所得到的岛的形状呈现出沿潮流方向的伸展(可由图 8中潮流场速度矢量分布看出), 呈现长条形状; 而没有考虑悬移质输沙的推移质输沙模式所导致的泥沙淤积更为集中于围绕岛中心点的区域, 岛的发展围绕岛中心点向四周扩张, 而不是沿潮流方向伸展, 最终近似于矩形。将图 8中计算所得湾中岛形状与图 1中实际湾中岛的形状对比可见, 推移质输沙率对应的矩形湾中岛形状与实际岛的形状更为接近, 只是面积要比实际的岛的面积要小(实际约1.25 km², 而图 8中推移质输沙模式25 a的面积0.76 km²)。注意到, 采用全沙输沙模式所得到的岛的面积(0.85 km²)也小于实际值, 可知产生这一差别的原因是由于本模拟中的形成湾中岛的泥沙来源完全来自于湾内地形冲刷, 而实际湾中岛的形成可能还有来自于岛外泥沙的贡献, 如风暴潮导致的泥沙在岛内的泥沙淤积。由此可见, 采用全沙输沙模式和采用推移质输沙模式所得到的岛的模拟结果各有其优缺点, 前者在面积大小上与实际岛的接近, 而后者在形状上与实际岛的接近。产生这些差别的原因是由于目前的输沙率公式还属于半经验公式, 其还不能够精确的表达实际海洋中复杂的输沙情况。

4 讨论

上一节计算模拟给出了两种不同模式下水东湾的地形演变特征, 即涨潮三角洲的形成演变趋势, 这里将对两种模拟过程做进一步的分析。

4.1 涨潮三角洲模拟分析

对比两种计算模式不难发现, 全沙输沙模式和推移质输沙模式状态下, 水东湾泻湖内均形成了一个独立的湾中岛。注意到两种模式均包含了推移质的影响, 这表明水东湾泻湖内涨潮三角洲形成的主要贡献是来自于推移质泥沙, 本计算也曾考虑过单独由悬沙所引起的湾内泥沙淤积, 即通过求解悬沙输移方程得到泥沙含沙量, 然后通过计算该含沙量所产生的悬移质输沙率来计算湾内地形的冲刷和淤积, 但结果表明仅考虑悬移质并不能够导致湾中岛的形成, 其原因是悬移质泥沙虽然被带入湾内, 但因为落潮流速大于涨潮流速, 其会随着涨落潮再次输出到湾外。这也是湾口外落潮三角洲形成的泥沙来源。以上结果从另一方面反映了推移质是形成湾中岛的主要贡献而不是悬移质。

悬移质贡献小的另一解释可以归因为湾内悬移质含沙量很小, 不足以形成导致泥沙大量沉积的环境条件, 其原因之一是由于水东湾外海海域水体平时比较清澈, 含沙量较低, 因而通常海况情况下潮流携带入湾内的泥沙量很少。

由于计算中没有考虑径流, 所以不存在径流的贡献。这一处理也是与实际情况符合的, 因为该湾内无大河注入。

4.2 落潮三角洲模拟分析

由图 5和图 7中的落潮三角洲的模拟结果可见, 落潮三角洲面积是不断向离岸方向伸展的, 外侧轮廓线呈不断扩张形态, 与之相伴随的是三角洲的高度也在不断减小(由于三角洲总体体积增长不大), 逐渐形成了扁平状的扇叶形三角洲。产生这一演化形态的原因是由于三角洲所在海域的落潮潮流流速大于涨潮潮流流速(图 4), 使得口门处的泥沙在落潮潮流携带下, 不断往离岸方向发展, 导致三角洲的面积不断扩大。本模拟仅考虑潮流作用, 没有考虑波浪的影响, 而后者会产生沿岸流和沿岸输沙, 从而为口门处淤积提供沙源。所以, 本文以上落潮三角洲演化模拟结果比实际的落潮三角洲高程要低。另一个原因是, 本模拟计算对落潮三角洲的初始地形取为了当前的落潮三角洲的地形, 这是因为没能找到在此之前的落潮三角洲的地形来作为本模拟的初始地形, 这样所给出的落潮三角洲地形演化实际上是今后25 a内的地形变化趋势, 即给出的是对今后发展态势的预报。

5 结论

本文数值模拟了水东湾涨潮三角洲和落潮三角洲地貌形态的产生和演变过程, 模拟采用潮流和水底地形演变耦合计算模型, 特别是再现了湾中岛(大洲岛)的形成过程, 所得研究结论如下:

模拟中湾中岛是由于涨潮流携带泥沙在湾内不断淤积而形成的, 这与该岛实际形成的物理过程符合。模拟所得到的湾中岛的面积要比实际的小, 原因可能是实际湾中岛的形成包含更多的泥沙来源, 这些来源在数值模拟中还没有全面考虑到。计算中考虑了全沙、单独推移质的两种输沙模型, 结果表明全沙模型能够得到较大面积的湾中岛, 这是由于全沙模型中包含了悬沙输沙的贡献, 说明悬移质在湾中岛沉积过程中也具有一定的贡献。但推移质对湾中岛形成的贡献是主要的, 这从单纯推移质输沙所产生的湾中岛面积也具有较大值看出, 特别是其形成的岛的形状是与实际岛的形状接近(比全沙的长条形岛更符合实际)。

以上水东湾的数值模拟结果表明, 适当设计计算模拟模型, 可以实现对海湾地貌形态演化的成功计算及数值模拟。但本文模拟仍然还有一些没有考虑的因素, 包括波浪对水东湾地形演变的贡献以及人类活动的影响(包括围垦填海), 这需要在今后的研究中进一步改善。