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

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

其中F_b为摩擦力，ρ为流体密度，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，该模式采用σ垂向坐标来模拟不规则的海底地形，坐标变换公式如下(Chen et al，2006)

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

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

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

对全日潮，平衡潮为

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

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。

计算的初始条件假设海洋是静止的，海表面的扰动水位初始值全部为0，所有三角形网格中心点的

流速u=v=0。本研究采用正压模式，取整个海域内的温度和盐度均为常数，温度为18°C，盐度为33。在台湾海峡、东海陆坡、朝鲜海峡共设有3个开边界，开边界上每个点的水位给定，由T_tide的预报程序给出(Pawlowicz et al，2002)，其中M₂、S₂、K₁和O₁调和常数来自DTU10(Cheng et al，2011)数据。

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

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

1.4 评价方法

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

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

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

2 实验方法 2.1 Chezy型摩擦

采用Chezy型摩擦时，底摩擦系数r为常数。采用不同的r值可以得到不同的F值，最优的r值对应最小的F值。为了找出最优的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值。

2.2 广义Manning型底摩擦

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

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

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

对应着最小值的(r₁，r₂)组合便是最优的组合。

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

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

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

3 模拟结果验证

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

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

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

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

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

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

4.2 潮能通量和潮汐传播

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

其中，和分别为流速和能通密度向量，T为潮波周期，ρ为海水密度，本文取1025kg/m³。将调和常数代入式(17)，可得

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

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

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

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

5 结语

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

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

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

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