Chezy型和广义Manning型摩擦关系在渤、黄、东海陆架潮汐模拟中的应用*

滕 飞 方国洪, 魏泽勋 徐晓庆 崔欣梅 吴 頔

(1. 中国海洋大学 青岛 266100; 2. 国家海洋局第一海洋研究所 青岛 266061)

为了研究底摩擦对潮汐的影响和提高潮汐数值模拟的精确度, 我们分别应用 Chezy型和广义Manning型摩擦关系模拟渤、黄、东海的潮汐。

Chezy于 1775年首先提出定常流动水渠中摩擦力与流速的关系, 后来被称为 Chezy公式(Dronkers,1964)。在潮汐动力学中可写成

其中Fb为摩擦力, ρ为流体密度, g为重力加速度, C为 Chezy系数, u为流速。如果 C为常数, 式(1)表明摩擦力与流速的平方成比例, 为方便常常用系数r代替Chezy系数, 取

r(亦即C)取常数的情形常被称为平方摩擦, 也可称为 Chezy 型摩擦(Lefèvre et al, 2000)。Taylor首次研究了爱尔兰海的潮汐摩擦问题(Taylor, 1920), 证明采用r=0.0016至0.002所得到的潮能消耗与依据潮汐和潮流观测值所算出的潮能消耗基本一致。Proudman(1953)按照Taylor的方法进一步给出r的具体数值为0.0026。此后, 大多数潮汐数值模拟中均取r=0.002至0.003的数值。

Manning于 1890年通过实验, 提出 Chezy系数并不是常数,而是与水深 h有关(Dronkers, 1964), 他给出的关系式,

称为Manning公式, 其中n称为Manning粗糙度系数,将式(3)代入式(2), 可得

式(4)表明, 水深越浅, r值越大。

对于中国近海, Fang等(1987)用三种不同的方法计算了杭州湾潮流的摩擦系数,得到的 r值介于0.0005至0.00067之间, 显著小于Proudman(1953)和后来的数模研究所采用的数值。值得注意的是,杭州湾水深很浅, 大约10m左右。周朦等(1987)估计了渤海的摩擦系数, 得到r=0.0013, 也比 Proudman(1953)的数值小很多, 而渤海的水深大约在20m左右。这些结果表明, 与其他模拟研究相比,渤、黄、东海海底摩擦系数较小, 而且与水深的关系可能与Manning公式不一致。我们称平方摩擦系数为常数的形式为Chezy型, 称平方摩擦系数与水深的幂函数成比例形式为广义Manning型。这里“广义”指幂值可以是任意实数,不局限于Manning公式中的-1/3。本研究尝试在渤、黄、东海潮汐模拟中分别应用 Chezy型以及广义Manning型摩擦关系,并通过与实测对比来确定其最优的参数值。

1 数值模式和实测数据

1.1 基本方程

本研究采用有限体积海洋模式 FVCOM, 该模式采用σ垂向坐标来模拟不规则的海底地形, 坐标变换公式如下(Chenet al, 2006)

由于计算范围大, 考虑纬度不同所引起的科氏力变化, 采用球坐标连续方程, 球坐标系下的σ坐标的连续方程和动量方程如下

式中,R是地球半径,γ代表计算点经度,Φ代表计算点纬度,D代表未受扰动水深,ζ为相对于未受扰动海面的高度,u、v、w分别代表东向、北向和垂向速度,ρ代表海水密度,f代表科氏参数,g代表地球重力加速度,Km代表垂直涡动黏性系数。本研究采用它的二维模式。

在大范围的区域中, 引潮势产生的平衡潮是不能忽略的, 因此模拟过程必须包括引潮力。对半日潮,平衡潮为

对全日潮, 平衡潮为

其中,Ai和ωi代表第i个分潮的振幅和频率,Nsemi和Ndiurnal代表模式中参与计算的半日潮和全日潮的个数,βi为地球弹性所产生的订正因子, 对各个分潮所采用的数值见表1。

表1 平衡潮订正因子βiTab.1 Correction factor for equilibrium tides βi

1.2 模式水深

在潮汐计算中, 水深是重要的因素, 且本文实验方法中广义Manning型底摩擦系数与水深关系密切。对于陆架海域, 直接从海图上读取的水深值一般要优于国外现有数据集的数值(例 如吴 頔 等, 2015), 因此为了提高计算准确度, 本文采用的水深数据来自两部分。第一部分是林美华等(1991)水深数据集, 该数据集的水深系由海图水深读取得到, 分辨率为 5′;第二部分是 Etopo5水深数据, 对于第一部分数据中没有涵盖的区域从 Etopo5水深数据中选取, 最后得到研究区域完整的水深数据, 其分布见图1。

1.3 计算海区、边界条件的选取与计算过程

本文研究海区为 24°—41°N, 116°—128°E 的渤海、黄海和东海陆架海域。模拟区域中不包括东海的深海区, 是由于深水区的潮能耗散主要来自内潮效应, 需要另作研究。本研究采用的计算网格水平分辨率在岛屿以及大陆边界附近约为 1/20°, 其余部分约为1/10°, 共包含26851个网格节点, 51486个三角单元(图2), 垂向分12层。FVCOM模型采用内外模分离的方式求解。二维外模数值格式是基于三角形网格的有限体积法, 将连续方程和动量方程在三角形单元内积分后, 通过改进的四阶龙格-库塔方法求解。三维内模的动量方程采用简单的显式和隐式相结合的差分格式求解。外模时间步长为6s, 内外模时间步长比率为10。

图1 渤、黄、东海水深分布图(单位： m)Fig.1 Bathymetry of the Bohai, Yellow, and East China Seas

图2 计算网格Fig.2 The computation grid

计算的初始条件假设海洋是静止的, 海表面的扰动水位初始值全部为0, 所有三角形网格中心点的流速u=v=0。本研究采用正压模式, 取整个海域内的温度和盐度均为常数, 温度为18°C, 盐度为33。在台湾海峡、东海陆坡、朝鲜海峡共设有 3个开边界,开边界上每个点的水位给定,由 T_tide的预报程序给出(Pawlowiczetal, 2002),其中 M2、S2、K1和 O1调和常数来自DTU10(Chengetal, 2011)数据。

模式运行30d。对后15d时间序列进行调和分析,得到网格点上的水位的调和常数、每个三角形中心点的潮流椭圆要素以及速度各分量的调和常数。

为了对模式结果和观测值进行比较并优化模式参数, 本文在研究区域内一共选取了43个实测站, 其分布图见图3。在这43个站的数据中, 26个来自国际水文组织IHO潮汐调和常数数据集, 11个来自TOPEX/Poseidon(T/P)卫星轨道交叉点调和常数数据(Wanget al, 2012), 2个站调和常数来自Jan等(2002), 4个根据TOGA计划的实测逐时水位资料经调和分析得到。

图3 实测站位置分布图Fig.3 Deployment of observational stations

1.4 评价方法

为了评估模式结果的准确性和优化模式参数,采用模式结果和实测调和常数的向量均方根(RMS)偏差F作为价格函数, 其计算公式如下

其中Fi为各个分潮的均方根偏差

上面k代表实测站位编号, i代表分潮, H代表振幅, G代表迟角, Com代表模拟结果, Mea代表实测结果。本文中对两种底摩擦类型的结果评价和参数优化都是以价格函数F值大小为标准, F值越小代表模拟结果越好, 反之越差。

2 实验方法

2.1 Chezy型摩擦

采用Chezy型摩擦时, 底摩擦系数r为常数。采用不同的r值可以得到不同的F值, 最优的r值对应最小的F值。为了找出最优的r值, 一种较先进的方法是通过伴随方程进行反演, 但是这会涉及难度较大的伴随模式的研发。对于当前的需要优化的参数很少的情况, 我们可采用多次直接模拟的方法。我们先用范围较大、分辨率较低(即间隔较大)的不同r值进行模拟。首先, 取r=(0.0005, 0.0010, 0.0015, 0.0020,0.0025, 0.0030, 0.0035, 0.0040), 结果得到其中r=0.0010时 F值最小; 然后将分辨率增加一倍, 取r=(0.00075, 0.0010, 0.00125)再进行模拟(其中对 r=0.0010不需要重新模拟, 但参与比较), 结果得到其中r=0.0010时F值最小; 最后再将分辨率增加一倍,取 r=(0.0008725, 0.0010, 0.001125)进行模拟(和前面一样, 其中对 r=0.0010不需要重新模拟, 但参与比较), 结果得到其中 r=0.001125时, F最小, 其值为11.04cm。本次实验F值分布曲线见图4 (其中黑色圆点代表模拟过程中各个r和对应的 F值)。由于最小的几个F值实际上差别很小, 故不作更高分辨率的进一步模拟, 而认为0.001125为最优的r值。

图4 采用Chezy型摩擦模拟所得价格函数F值与摩擦系数rFig.4 The dependence of cost function F on friction coefficient r in Chezy-type relationship

2.2 广义Manning型底摩擦

为了避免当h→0和h→∞时r→0,将广义Manning型摩擦系数r取作如下形式

若 h≤h1, 取 r=r1;

若 h≥h2, 取 r=r2; (14)

为满足作为h之函数的r的连续性, a和m需分别等于

本研究中取1h=10m,2h=200m。数值模拟结果与实测值的偏差与r1和r2有关。因而价格函数F可以写成

对应着最小F值的(r1, r2)组合便是最优的组合。

为了寻找最优的(r1, r2)组合, 先取范围较大、分辨率较低的(r1, r2)组合： (r1=0.0002, 0.0004, …, 0.0014;r2=0.0005, 0.0010, …, 0.0040)这 56 组(r1, r2)值进行模拟, 并计算出对应的 F值。结果得到(r1, r2)=(0.0008,0.0015)时 F值最小。然后以该点为中心, 缩小(r1, r2)的取值范围并提高一倍分辨率继续进行数值模拟, 并计算对应的F值, 结果得到(r1, r2)=(0.0009, 0.0015)时F值最小。最后再以该点为中心, 进一步缩小(r1,r2)的取值范围并再提高一倍分辨率继续进行数值模拟, 并计算对应的F值, 结果得到(r1,r2)=(0.0009, 0.001375)时F值最小, 等于 10.88cm。由于在这个组合附近F值的变化已很小, 我们认为它已经是最优的参数组合,不再作更高分辨率组合的模拟。由式(15)可得a=0.001375,m=0.14。与原始的Manning公式不同, 该公式的幂值为负数, 而本研究得到的幂值为正。这表明在渤、黄、东海陆架区, 水深越浅摩擦系数r值越小。

全部实验结果以F值等值线的形式示于图5, 其中黑色圆点为(r1,r2)第一次实验取值点, 正方形点为第二次实验取值点, 菱形点为第三次实验取值点。为了更好显示第二和第三次实验结果, 我们将图5方框区域予以放大, 示于图6。

图5 采用广义Manning型摩擦模拟所得价格函数F值与参数组合(r1, r2)的关系： 全部实验结果Fig.5 The dependence of cost function F on parameter combination (r1, r2) in Manning-type relationship： all results

图6 采用广义Manning型摩擦模拟所得价格函数F值与参数组合(r1, r2)的关系： 图5方框区域放大图Fig.6 The dependence of cost function F on parameter combination (r1, r2) in Manning-type relationship： the enlarged view of the rectangular area in Fig.5

通过比较F值可以发现, 采用广义Manning型摩擦关系得到的结果优于Chezy型摩擦关系。对比用两种摩擦关系所得的同潮图也可以看出, 采用广义Manning型摩擦关系结果更合理。例如对M2分潮, 由于Chezy型关系给出的摩擦系数在深海区偏小, 计算结果与验潮站实测值比较在海域东南部, 特别是济州岛附近和朝鲜半岛西南沿岸显著偏大; 在渤海由于Chezy型关系给出的摩擦系数偏大, 模拟得出的两个无潮点可能离岸太近甚至消失(例如见Lefèvreet al,2000; 朱学明等, 2012)。因此下文的分析均采用由最优参数[即(r1,r2)=(0.0009, 0.001375)]广义Manning型摩擦关系计算所得的结果。

3 模拟结果验证

本文选取的43个代表性实测站大体均匀地分布在研究区域内, 我们在计算网格中找到与实测站距离最近的节点, 比较计算和实测所得的潮汐调和常数, 见表2(表中H代表振幅, 单位为cm;G代表格林尼治迟角, 单位为度(°)。由表可以得出, M2、S2、K1和 O1分潮振幅/迟角偏差绝对值的平均值分别为11.4cm/6.4°, 4.3cm/6.1°, 4.3cm/6.1°和 1.1cm/3.7°。总体上模式和实测结果具有较好的一致性。

表2 四个主要分潮计算和实测结果对比表(H: 振幅, G: 迟角)Tab.2 Results of computation and observation on harmonic constants for four principal constituents (H： amplitude, G： phase-lag)

4 渤、黄、东海海域的潮波分布

4.1 潮汐分布

对数值模拟的水位场进行调和分析, 得到 M2、S2、K1和 O1四个分潮的调和常数。根据这些调和常数绘制的同潮图分别示于图7、图8、图9和图10。对于渤、黄、东海潮汐、潮流的观测和研究工作已有很多, Fang(1986)根据大量的实测资料以及数值模拟结果绘制了一份较完整的潮汐、潮流图; 王凯等(1999)曾经用三维数值模式模拟了该区域的 M2分潮的潮汐、潮流(王凯等, 1999); Fang等(2004)基于10年的卫星高度计资料给出了该区域的潮汐同潮图(Fang et al,2004); 王永刚等(2004)用验潮站资料做了渤、黄、东海的同化数值模式; 朱学明等(2012)用 FVCOM 海洋数值模式模拟了该区域的潮汐、潮流。

图7 M2分潮振幅(cm)和格林威治迟角(°)Fig.7 The amplitude (cm) and Greenwich phase-lag(°) of M2

由M2分潮的同潮图(图7)可以看到, 本文得到的黄河口外海无潮点的位置上与 Lefèvre等(2000)及朱学明等(2012)的结果有区别。后两项研究给出的无潮点已退化到了陆地上, 而本文的结果显示无潮点不但存在, 而且离岸有一定距离。这主要是由于黄河口岸形的改变引起的, 与王永刚等(2014)模拟得出的新近情况一致。其他三个无潮点(渤海西北部、山东半岛东侧海域和黄海南部)的位置与前人结果无太大差别。此外, 在台湾岛东北有一个明显的退化的无潮点,这与前人结果也一致。

图8是模拟得到的S2分潮同潮图, 与M2分潮一样, 一共有四个无潮点。与Fang等(2004)的结果相比,在本文中渤海内部的两个无潮点完全形成, 而后者基本在海岸线上, 其他地方等迟角线差别不大。关于S2分潮的等振幅线在苏北外海的形状, Fang等(2004)的结果中等振幅线有一个明显的北向凸起, 本文则没有。但是本文的结果与王永刚等(2004)和朱学明等(2012)的结果非常接近, 这可能是由于Fang等(2004)的研究中所用的卫星高度计地面轨道较稀疏, 内插引入的误差造成的。从半日潮来看, 本文的结果与Fang等(2004)实测的结果相比, 除在渤海海域差别较大, 其他地方基本一致。

图8 S2分潮振幅(cm)和格林威治迟角(°)Fig.8 The amplitude (cm) and Greenwich phase-lag (°) of S2

图9和图10给出了K1和O1分潮的振幅和格林威治迟角分布。它们与Fang等(2004)结果非常接近。总体来说, 与半日潮相比, 全日潮在该海域结构较简单, 振幅较小, K1分潮最大 40cm 左右, O1最大在30cm左右, 因此模拟结果较好。

图9 K1分潮振幅(cm)和格林威治迟角(°)分布图Fig.9 The amplitude (cm) and Greenwich phase-lag (°) of K1

图10 O1分潮振幅(cm)和格林威治迟角(°)Fig.10 The amplitude (cm) and Greenwich phase-lag (°) of O1

4.2 潮能通量和潮汐传播

一个潮周期内单位时间通过自海底至海面单位宽度断面的潮能通量叫做能通量密度, 其算式为

其中Φx和Φy分别为潮能通量密度的东向和北向分量,H和G分别为水位的调和常数, U、ξ和V、η分别为潮流东向和北向分量的调和常数。

本文根据式(18)和式(19)计算了四个分潮的潮能通量, 其中最大的半日分潮M2和最大的全日分潮K1的能通量分布见图11和12。为了显示潮能通量区域的特征, 图中矢量的长度与能通量的平方根成正比例, 其大小比例尺见图的左下方。

从图11的M2分潮来看, 潮能输送路径主要分为三支： 第一支是通过东海陆架边缘南侧进入东海后反转, 由台湾海峡输入南海; 第二支是通过东海陆架边缘北侧经东海进而到达黄海, 在苏北外海和山东半岛南部海域形成一个逆时针旋转的潮波系统进而耗散; 第三支是沿着朝鲜半岛西海岸一直向北进入渤海耗散。从图可以看出, 大部分潮能是在黄海消耗的, 进入渤海的潮能较小; 同时在渤海海峡处, 半日潮能通量均向西, 即使在海峡南部也看不到向东的反射波能通量, 这与全日潮波有明显不同(见下)。S2分潮能通量分布特征与 M2基本相同; 量值较小, 只有M2的10%—20%。

与半日潮相比, 全日潮能量较小(其中 K1分潮见图12)。全日潮的输送路径分三支： 一支由台湾岛东侧由东海反转进入台湾海峡; 另一支在黄海中部向西北方向前进, 沿山东半岛南部沿岸和苏北沿岸向南发生耗散; 还有一支在济州岛南部沿着朝鲜半岛西岸向北前进, 通过渤海海峡北部进入渤海耗散大部分, 其余部分通过渤海海峡南部回到黄海, 并在山东半岛北侧外海消耗殆尽。全日潮波南向能通量基本上可以到达长江口北岸, 而半日潮波南向能通量只能到达苏北外海中部。这说明从反射波相对于入射波的能量来看, 全日潮要比半日潮强。

图11 M2分潮能通量密度分布图Fig.11 The energy flux density vectors of M2

图12 K1分潮能通量密度分布图Fig.12 The energy flux density vectors of K1

5 结语

基于FVCOM海洋数值模式, 建立了渤、黄、东海陆架海域的潮汐数值模型, 分别采用Chezy型和广义Manning型底摩擦系数模拟了潮汐、潮流, 并通过与实测数据对比确定了其最优参数。将计算结果与实测数据进行比较, 结果良好。最后, 用广义 Manning型底摩擦系数最优结果对该海域的潮汐、潮流和潮能进行了分析。主要结论有：

(1) 采用广义Manning型比采用Chezy型底摩擦系数得到的结果更好, 因此建议在模拟渤、黄、东海潮汐时最好采用前者。

(2) 本研究得到的底摩擦系数在0.009至0.0014之间, 显著低于 Proudman(1953)得到的 0.0026, 这也表明渤、黄、东海海底摩擦系数要比大多数模拟研究中采用的数值要小。

(3) 原始的Manning公式的幂值为负数, 而本研究得到的幂值为正。这表明在渤、黄、东海陆架, 总体上水越浅摩擦系数 r值越小, 也表明原始的Manning公式必须推广才能适用于渤、黄、东海大范围潮汐模拟。

王 凯, 方国洪, 筰冯士, 1999. 渤海、黄海、东海M2潮汐潮流的三维数值模拟. 海洋学报, 21(4)： 1—13

王永刚, 方国洪, 曹德明等, 2004. 渤、黄、东海潮汐的一种验潮站资料同化数值模式. 海洋科学进展, 22(3)： 253—274

朱学明, 鲍献文, 宋德海等, 2012. 渤、黄、东海潮汐、潮流的数值模拟与研究. 海洋与湖沼, 43(6)： 1103—1113

吴 頔, 方国洪, 崔欣梅等, 2015. 泰国湾及邻近海域潮汐潮流的数值模拟. 海洋学报, 37(1)： 11—20

林美华, 方国洪, 1991. 中国海标准经纬度水深和基准面数据表. 青岛： 中国科学院海洋研究所, 1—144

周 朦, 方国洪, 1987. |U|U的Fourier展开和渤海海底拖曳系数 CD. 海洋与湖沼, 18(1)： 1—11

Chen C S, Beardsley R C, Cowles G, 2006. An Unstructured Grid,Finite-Volume Coastal Ocean Model-FVCOM User Manual,Technical Report SMAST/UMASSD-06-0602. 2nd edn. New Bedford： School for Marine Science and Technology,University of Massachusetts Dartmouth, 6—8

Cheng Y C, Andersen O B, 2011. Multimission empirical ocean tide modeling for shallow waters and polar seas. Journal of Geophysical Research Oceans, 116(C11), doi： 10.1029/2011JC007172

Dronkers J J, 1964. Tidal Computations in Rivers and Coastal Waters. Amsterdam, North-Holland： Interscience Publishers,518

Fang G H, 1986. Tide and tidal current charts for the marginal seas adjacent to China. Chinese Journal of Oceanology and Limnology, 4(1)： 1—16

Fang G H, Li S D, Cao D M, 1987. Three methods for estimating turbulent stress and drag coefficient in tidal currents of the Hangzhou Bay. Science Bulletin, 32(4)： 252—257

Fang G H, Wang Y G, Wei Z X et al, 2004. Empirical cotidal charts of the Bohai, Yellow, and East China Seas from 10 years of TOPEX/Poseidon altimetry. Journal of Geophysical Research： Oceans, 109(C11), doi： 10.1029/2004JC002484

Jan S, Chern C S, Wang J, 2002. Transition of tidal waves from the East to South China Seas over the Taiwan Strait：Influence of the abrupt step in the topography. Journal of Oceanography, 58(6)： 837—850

Lefèvre F, Le Provost C, Lyard F H, 2000. How can we improve a global ocean tide model at a regional scale? A test on the Yellow Sea and the East China Sea. Journal of Geophysical Research： Oceans, 105(C4)： 8707—8725

Pawlowicz R, Beardsley B, Lentz S, 2002. Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE. Computers & Geosciences, 28(8)： 929—937 Proudman J, 1953. Dynamical Oceanography, London： Methuen and Co, 409

Taylor G I. 1920. Tidal friction in the Irish Sea. Philosophical Transactions of the Royal Society of London. Series A, 220：1—33

Wang Y H, Fang G H, Wei Z X et al, 2012. Cotidal charts and tidal power input atlases of the global ocean from TOPEX/Poseidon and JASON-1 altimetry. Acta Oceanologica Sinica,31(4)： 11—23