红树林海湾潮汐潮流特征数值模拟研究

2021-11-24李谊纯刘世通

人民珠江 2021年11期

李谊纯，刘世通

(北部湾大学建筑工程学院，广西钦州 535011)

红树林海岸是生物海岸的一种，在中国广泛分布于低纬度的广西、广东、海南及福建等沿海地区。红树植物密集生长于潮间带上，形成的红树林对于防护海岸、保养潮滩乃至水产经济具有重要的意义。鉴于红树林系统在近岸区域的重要角色，相关研究颇为丰富。然而相对而言，红树林水动力学的研究相对于红树林研究领域中的其他分支，进展较为缓慢。究其原因，主要有：红树林海岸中的水文泥沙等相关实测资料准确获取的难度大；相对而言，红树林是空间尺度较小的研究对象，地形、地貌、水文、生物等要素的空间差异较大，不利于形成较为普适的统一结论；部分关键物理量难以参数化，尤其难以在水动力计算中予以概化。Wolanski是红树林水动力学研究的先驱，也较早地应用数学模型开展了天然红树林水动力过程的模拟研究[1-2]。此后，Kjerfve、Golbuu、Furukawa等学者在此领域做了一些探索性研究[3-6]。其中，Mazda等[7]建立了理想红树林海湾的数值模拟研究，阐释了潮流过程的“双峰”现象。Wolanski等[8]利用二维数学模型模拟了Hinchinbrook湾的水动力过程，然而其网格尺度较大(344 m)，导致在数值模拟中，岸边的红树林仅对水体的连续性起作用，而对红树林水动力的特征刻画不显著。近年来，亦有利用三维数学模型的类似应用，但结论较为初步[9]。另一方面，河口海岸学中的潮汐/潮流不对称理论可为红树林海湾水动力研究提供重要参考，如：在宽而浅、潮滩广阔的河口，潮汐可呈现落潮主导型；在窄而深的河口，潮汐多呈现涨潮主导型。为了探讨小尺度红树林海湾的水动力特征，建立了二维水动力数学模型，在理想地形上进行了数值模拟研究。就分布位置而言，红树林可分为凹地型、河岸边缘型、海岸边缘型[10]。在此，本文仅对河岸边缘型红树林展开研究，以期对河口防护及环境安全提供参考。

1 数学模型简介

1.1 模型方程

红树林海湾的水深一般很浅，平面尺度亦不大。因此，选择平面二维水动力学模型对红树林海湾进行模拟，暂不考虑近岸盐度变化导致的斜压效应以及水动力要素的垂向变化特征。平面二维水动力学控制方程可写为：

(1)

(2)

(3)

式中η、u、v、H——水位、x和y方向的流速、总水深；f——柯氏力系数；n——曼宁系数。

在数值模拟中，x和y方向分别设置为海湾的纵向和横向，正方向分别指向东向和北向。海湾不考虑横向变化，因此在数值模拟中柯氏力系数f取零。

1.2 方程的离散

水动力控制方程的数值离散格式关乎求解的精度。在此，动量方程的对流项采用Van Leer 限制器处理以使其数值精度在理论上达到二阶，连续方程采用显格式求解。红树林水动力的数值模拟中，底摩阻是关键的作用力。红树林的曼宁系数约为0.1～0.7，一般约为0.3～0.4，这高于正常底床糙率一个量级[8]。因此，离散方程中底摩阻项的处理尤为关键。在此，动量方程的时间积分采用时间分裂两步格式处理：第一步仅隐式求解底摩阻项，可保证当底摩阻很大时，流速只相应变小而不产生虚假流动；第二步则对动量方程中的其他诸项采用显格式求解。离散方程可表示如下，第一步求解：

(4)

(5)

第二步求解：

(6)

(7)

式中 Δt——时间步长；m——时间层。

为与式(2)(3)保持一致，式(6)(7)中保留了柯氏力项，具体计算中f取零。

1.3 动边界处理

红树林生长于潮间带之上，加之潮滩坡度较缓，因此岸边界随着潮位涨落变动显著。动边界的处理方法有很多种，如虚拟水深法、冻结法、水位平铺法、渗透介质法等[11]。综合而言，基本可以分为两类：干湿网格法和窄缝法。干湿网格法是在计算过程中判断计算点(网格)的总水深，进而确定该点(网格)是否参与计算以实现对岸边界变动的处理。干湿网格法在计算中对水位点或流速点进行判断，有多种实现形式，因此可认为干湿网格法是一类方法。这类方法具有实现灵活，移植性强的特点。然而，由于具体判断方法和判断准则的原因，干湿网格法对岸线移动速度的控制可能存在较为明显差异[12]。干湿网格类方法本质上是“间断”的方法，岸边界的移动依赖于计算点的干湿判断。窄缝法是在岸滩上人为设置窄缝，将水引入窄缝模拟岸边线变动的一种方法。窄缝法通过“化引水深”实现方程在固定范围内模拟计算并实现岸边界变动，本质上是一种“连续”的方法。陶建华等[13]最早将其应用于波浪爬坡的计算，何少苓等[14]将其应用于近岸潮流模拟研究，孙琪等[15]将窄缝法应用于近岸泥沙输运的研究中，均取得了较好的效果。Shi等[16-17]应用窄缝法于近岸的波浪计算取得了一定的效果，但发现窄缝的选取是需要注意的问题。图1为窄缝法的示意，式(8)为窄缝宽度系数的计算公式。“化引水深”是考虑了窄缝后的等效水深，具体计算方法可参考文献[13-16]。

(8)

式中f(z)——窄缝宽度系数；z、zb——水面高度和河床高度(图1)；ε——窄缝参数，代表窄缝最低处的宽度；α——常数，控制窄缝缩窄的速度。

图1 窄缝法示意(z0为人为设置的窄缝最低点，应低于最低水位)

窄缝法允许水在窄缝内流动，利用“化引水深”代替总水深，实现岸边界的自动判别。Shi等[16-17]认为窄缝参数确定不当可能导致计算失稳，此外也可能导致“化引水深”与总水深出现较大差异。在此，通过调整ε和α的值将最小“化引水深”控制为约0.01 m。此外，由于红树林根系及生物活动(如蟹孔)的存在，潮滩在涨落潮过程中较一般的潮滩存在更多的水的渗入和排出。窄缝法中“化引水深”的引入，在一定程度上模拟了红树林潮滩上水的渗入和排出。

2 计算区域及计算方案

2.1 计算区域

图2为计算区域的示意，红树林海湾长度约5 km，中间深槽深度为5 m，深槽两侧红树林潮滩宽度约为600～800 m，海湾顶部红树林宽度150 m。湾口外水域约为27 km×19 km，开边界仅在东侧设置，给定潮位过程。图2中S1—S3为采样断面，分别代表湾口、中部及湾顶。图3为红树林海湾的横断面示意。图中h为深槽水深，b为深槽宽度，α为潮滩自滩肩向两侧的坡度。L1—L3为图2中所示断面上的采样点位置，分别位于深槽中部、滩肩及潮滩内部。假设海湾深槽及潮滩的水深自湾口向湾顶不存在纵向变化，湾口外水深自湾口向外海线性变化(坡降为0.000 5)。

图2 计算区域及水深示意(m)

图3 海湾横断面示意

2.2 计算条件及计算方案

为研究不同因素对红树林海湾水动力特征的影响，将考虑深槽宽度、潮滩坡度以及曼宁系数的作用。深槽宽度反映槽滩体积比，是潮汐不对称的重要指示性参数；潮滩坡度对海湾槽滩的涨落潮过程具有重要影响；曼宁系数是量化红树林及其密度的直接参数。为此，根据广西钦州湾红树林潮汐水道情况及参考前人文献[7],设置红树林湾内的深槽宽度、潮滩坡度和曼宁系数的取值见表1，共计180组模拟计算。计算网格全域为矩形网格。在红树林海湾内深槽横向设置10个网格，相应的空间步长最小为5 m、最大为40 m。海湾内潮滩上的空间步长为50 m，海湾外水域空间步长逐渐增大，外海开边界处为500 m。开边界只在东侧开边界取潮位过程控制，不考虑上游径流注入。开边界潮差取1.6 m，中潮位位于滩肩下0.4 m。计算基面置于滩肩所处平面。为数据处理的简便，潮周期取24 h。

表1 计算条件及计算方案

3 结果与讨论

3.1 潮流基本特征及红树林的影响

红树林的存在对潮滩上及深槽的水流过程会导致明显的改变。一般的，对于潮滩宽阔的河口海岸，由于潮滩在水位较高时能容纳较多的水体，可能导致潮汐呈落潮主导型，即落潮历时小于涨潮历时、落潮最大流速大于涨潮最大流速。在深且窄、潮滩不发育的河口海岸，潮汐一般呈现涨潮主导型。对于小尺度的红树林海湾，潮滩分布宽阔(图3)，具备产生落潮不对称型潮汐的潜在条件。在此给出并对比了红树林小尺度海湾的潮汐潮流过程的基本特征及其变化。

图4、5分别给出了潮滩坡度为0.001 5情况的采样点上(图1、2)的潮位和流速过程。由于深槽处横向流速很小，图中未给出。对于潮位过程，可以看出红树林无明显影响。另一方面，对于流速，红树林产生了非常明显的影响。首先，对于滩肩(L2)和滩面上(L1)的采样点，在没有红树林的情况下，纵向流速u可超过0.1 m/s，最大可近 0.5 m/s (图4)；而在红树林的影响下，纵向流速减弱至小于0.05 m/s(图5)。对于横向流速v，红树林的存在与否并无明显量值的变化，一般小于0.05 m/s。由此可以认为，红树林对潮滩上流速的影响主要表现为：纵向流速明显减弱，横向流速变化很小；潮滩上水体由纵向输运为主变为横向和纵向输运并重。对于深槽中(L1)的流速，红树林的存在导致纵向流速有所增加。需要注意的是，在湾口断面(S1)的深槽处，流速虽在量值上相差不大，但却呈现不同的涨落潮不对称性质。在无红树林情况下(图4)，最大涨落潮流速分别为-0.64、0.74 m/s，最大落潮流速大于最大涨潮流速，呈现落潮流主导的不对称型式。而在有红树林存在的情况下(图5)，最大涨落潮流速分别为-0.91、0.87 m/s，最大涨潮流速大于最大落潮流速，呈现弱的涨潮流主导的不对称型式。潮流不对称型式的变化对于水体及物质输运具有重要的指示意义，对于红树林尤其重要。在此处的模拟中，外海边界为对称型的潮位过程，在潮波向近岸传播的过程中，由于水深变浅、底摩阻等因素的影响，会逐渐向涨潮主导型发展。在海湾口门处应表现为涨潮主导型。有红树林存在的情况下的结果符合这一特征。在无红树林的情况下，湾口处的潮流为落潮主导型。如前所述，宽阔的潮滩可促使落潮主导型产生，湾口处的流速表现为此特征。但在海湾内部的潮流仍呈涨潮主导型，这种不对称型式的变化可能是由于纵向变化隐含着更为复杂的机制。

图4 潮滩无红树林时(n=0.025)采样点的潮位和流速过程

图5 潮滩有红树林时(n=0.3)采样点的潮位和流速过程

图6、7分别给出了有无红树林情况下海湾水量输运情况。对比可以发现在无红树林的情况下，平均水体输运在潮滩上有明显的向上游的净输运，在深槽则为向下游的净输运。这一特征符合河口海湾正压流的基本特征[18]。而在有红树林的情况下，此特征显著减弱。这种现象应是前述深槽处潮流不对称型式变化的直接原因。此外，由图6、7还可看出，在涨潮时，无红树林潮滩存在较强的纵向输运，横向输运相对较弱。红树林潮滩则基本表现为相反的特征。落潮时，无红树林的情况下，纵向输运则主要集中于深槽处，而滩面上的横向输运则在有红树林时更为突出，且分布面积更大。

a)纵向输运 b)横向输运

图8、9分别为潮滩是否存在红树林情况下的潮位-流速，其中流速仅给出纵向流速。从图中可以得出潮流在涨落潮过程中的变化。可以看出，红树林的存在导致深槽纵向流速的增大及潮滩流速的大幅减小。此外，最大流速出现于中潮位和高潮位之间，表明潮波既非驻波亦非行进波，而是处于二者之间的混合型。由于潮滩的存在，最大涨落潮流速出现在较高水位，这符合一般情况下的河口近岸水域的潮波特征[19]。图8、9还显示涨潮最大流速一般出现在潮位为0.5～1.0 m，落潮最大流速则多出现在潮位0.0 ～0.3 m。

图8 无红树林情况下(n=0.025)的潮位-流速

图9 有红树林情况下(n=0.3)的潮位-流速

为进一步研究最大流速与潮位的关系，从水量平衡出发建立一个简单的理论模型，研究了最大流速出现时的潮位情况。对于此类地形，湾口的断面平均流速可表示为[19]：

(9)

式中η——潮位；ua——断面平均流速；As、Ac——海湾的表面积和横断面面积。

对于图3所示的断面可有:

(10)

其中，T=η-h，L为河(湾)长。符号见图3。上式右端对时间求导并令其为零，可得最大流速出现的时刻以及最大流速的值，但式(10)较难求得此解析解。在此，考虑到河口潮波多为半驻波半行进波的混合型式，最大流速一般出现在中潮位前后的时刻，所以假设在此时段内η随时间t线性变化。因此，求解：

(11)

(12)

(13)

(14)

因为在一般情况下有：tanβ>>tanα

(15)

(16)

由此可以看出，最大流速出现的潮位与深槽形态(深宽比tanβ)、潮滩坡度(tanα)以及深槽尺度(h或b)相关。此外，式(11)—(16)中假设了潮位是随时间线性变化的，然而对于不同的涨落潮过程，此变化的速率存在差异，如涨落潮不对称导致的水位变化的快慢、半日潮海区和全日潮海区的差异等。同时，式(9)假设了整个河口(海湾)的水位是同步涨落且没有水面坡度。因此，式(16)仅是从水体连续性出发得出的最大流速出现时刻的一个简单估计，更细节性的关系尚需进一步研究。对上述所数值模拟的情况：tanβ=0.1，tanα=0.001 5，b=50 m，可得T≈0.57 m。此数值与图8或图9所示情况呈现较好的一致性。但是式(16)未考虑涨落潮不对称的变化，所以其估计未能表现图中所示的涨落潮过程中最大流速出现的水位存在差异的现象。

为进一步探讨流速和潮位的横向变化，图10、11给出了S2 断面的流速情况，因为采用了对称地形，所以图中仅给出了自深槽中心向北的半个断面的流速情况。首先，定义潮流不对称指标为：

(17)

式中Uem、Ufm——最大落潮流速和最大涨潮流速；U——既可以代表纵向流速也可代表横向流速；TA——潮流不对称指标，当其大于1时，代表落潮流主导，反之则为涨潮流主导。

S2断面位于红树林海湾中部，图10、11均显示在此断面，纵向涨落潮流速均自滩肩向潮滩内部逐渐减小。潮流不对称均表现为涨潮不对称型，这有利于物质向上游输运。对于横向流速，在无红树林时，横向流速在断面上无明显变化趋势，但红树林会导致横向流速自滩肩向潮滩内部缓慢减小。潮流不对称则表现为在邻近深槽部分为落潮主导型，而在潮滩的上部则为涨潮主导型。这表明在潮滩的中上部有利于物质的向岸输运，但在滩肩及其邻近部位则可能向深槽净输运。红树林的存在使此转折点明显向深槽移动，表明红树林的存在有利于潮滩的发育。但是，需要注意的是，横向流速一般很小，不足以起动潮滩滩面上的泥沙，因此，上述的潮流不对称应主要作用于水体中已存在的悬沙及其他物质。

图10 无红树林情况下(n=0.025)的断面最大流速分布

图11 有红树林情况下(n=0.3)的断面最大流速分布

在红树林潮滩上，水流缓慢，可以预计压强梯度力是流速的决定性动力之一。图12给出了S2断面上L2点的涨落潮过程中的水面比降的变化。可以看出，在无红树林时，涨潮过程中纵比降一般大于横比降。在有红树林时，纵比降无明显变化，但横比降显著增大，超过纵比降。这与前述的无红树林时纵向流速主导输运一致。总体而言，红树林的存在增强了横向流动的作用。

图12 水面比降的时间变化

3.2 潮滩坡度和红树林密度的影响

图13给出了S2断面L1 点的潮位和纵向流速随潮滩坡度的变化。可以看出，随着坡度的增加，潮差逐渐增大，高潮位的增幅小于低潮位的降幅。在坡度为0.000 5时涨潮历时较小，其他情况变化不大。但是落潮过程中具有明显的差异。坡度越大，潮位下降越快，落潮历时越短。落潮最大流速小幅增大，但涨潮最大流速增加明显。落潮过程中落潮流速随坡度的减小而显著增大。这应源于坡度小的情况下，红树林海湾内潮棱体的体积显著地增大。相应的涨潮过程中的流速亦表现为随坡度减小而增大。但涨潮过程中的变化较落潮过程小。因此，可以认为，对于此类小尺度的海湾，涨潮过程主要受海湾湾口处潮位控制，而落潮过程则主要是湾内的水体体积的影响。湾口处的过水面积是另一个关键的要素。图13中，并未表现出来一般意义上的落潮主导型。原因可能由于湾口处入射潮波为涨潮主导型的影响以外，湾口处过水断面的面积应为一个控制因素。有红树林存在情况下，潮滩上的纵向水体输运大幅减小；虽然随着水位的增高，湾口断面宽度增大，但主要的过水断面仍为深槽的宽度。所以当湾内的体积增大时，湾口处的涨潮流速亦增大。而落潮过程与之类似，湾口过水断面的有效宽度基本不受潮滩坡度的影响，而湾内水体体积的增大导致落潮历时的延长和落潮流的增大。图14为滩肩和滩面上的横向流速变化，均表现为随坡度减小而逐渐增大的趋势。其原因亦可归于湾内水体体积的变化。

此外，据前人研究，红树林密度不同，其导致的阻力亦不同，若以曼宁系数表达，则约为0.1～0.7。但是，对于S2断面上的采样点的比对显示，在n=0.1～0.5的区间内，潮位和流速过程变化很小(图未列出)。结合前述，红树林的有无对水流影响显著，但红树林的密度影响相对较小。此外，红树林对水流实际的阻力作用于整个水体(水柱)，而曼宁系数仅是将其概化为底摩阻。虽然所采用的模型为垂向平均的二维模型，但是这种阻力的不同作用方式亦可能存在潜在的差异，尚需进一步研究。

图13 S2断面L1 点的潮位和纵向流速随潮滩坡度的变化

图14 S2断面L2和L3 点横向流速随潮滩坡度的变化

3.3 深槽宽度的影响

图15 给出了S2断面L1 点的潮位和纵向流速随深槽宽度的变化。可以看出，随着宽度的增大，潮差逐渐增大，但主要表现为高潮位的抬高。从涨落潮历时来看，随着宽度增加，涨潮历时略有减小，潮汐过程向涨潮主导型发展。这时由于深槽宽度增大，减小了滩槽体积比。流速则表现为显著的减小，从趋势上亦表现出涨潮最大流速逐渐大于落潮流速的变化。图16给出的滩肩上(L2)和潮滩上(L3)的横向流速过程随深槽宽度的变化。滩肩流速总体而言变幅不大且无明显规律，潮滩滩面上流速随着宽度的增大而逐渐增大，这可能与深槽中的潮差增大导致的涨落潮过程中横向比降增大有关。此外，在潮滩上，水位受滩面地形的影响，其变幅不如深槽中显著。