基于FVCOM 多种定解条件的数值分析与评估

陈昌刚，丁 扬，姚志刚，张 丛

（1. 物理海洋教育部重点实验室，山东 青岛 266100；2. 中国海洋大学 海洋与大气学院，山东 青岛 266100；3. 山东省海洋科学研究院，山东 青岛 266104）

作为海洋中的主要动力过程，环流直接影响温度、盐度、密度等要素分布。研究海洋环流有助于加深对海洋变化过程的了解。随着科学技术的发展，海洋数值模型成为一种新兴手段，被广泛应用于海洋环流的研究领域[1-4]。

中国近海环流的数值研究可追溯到20 世纪70 年代。从热盐环流的“风旋度-热盐梯度”二维简化模型方程[5]，到基于正压模式的浅水方程[6]，再到三维斜压模式[4,7-12]，均取得了丰硕成果，对认识与理解海洋起到至关重要的作用。例如，Huang 等[7]应用汉堡海洋数值模型HAMSOM 研究渤海斜压环流并指出，渤海环流是由潮、风应力以及斜压作用共同决定的。Xia 等[11]基于波-潮-环流耦合模型MASNUM，充分考虑在表面混合效应下研究黄海环流结构，提出黄海夏季水平环流具有3 层结构。Gan 等[1]基于区域海洋数值模式ROMS（Regional Ocean Model System）研究，验证中国海域的环流具有显著的分层结构，且常年存在，强度随季节变化。

Uda[13]基于观测数据首先提出黄海暖流，并认为黄海暖流是对马暖流的北向分支。乐肯堂和毛汉礼[14]分析中美联合调查资料，提出黄海暖流分为西北支与北支两个分支，称为“双峰结构”。之后，多位学者对该结果进行了深入研究，例如，赵胜等[15]利用卫星反演温度手段统计1981 年至2010 年黄海暖流分支出现次数，提出分叉结构具有普遍性。熊学军等[16]基于南黄海西部海区水体环境综合调查资料提出黄海暖流具有“后部显著、中部分叉、颈部收缩、顶部膨胀”的形态结构。

众所周知，不同定解条件对应微分方程不同的解，不同温盐初始场或开边界强迫场也会带来不同质量的数值结果[17-18]。例如，Peng 等[19]基于普林斯顿海洋模型的敏感性实验指出，相应地调整初始条件可以显著改善风暴潮的模拟结果，Marchesiello 等[20]认为选择适当开边界的重要性可比拟是否选择开边界，Sun 等[21]指出北赤道流分叉点位置模拟不佳可以归因于开边界的选择，并认为该结论可推广至低纬度西边界流的模拟。

基于上述认识，本文使用不同海洋数据作为温盐初始场和开边界强迫场，相互搭配驱动数值模型，对比结果与实测数据，研究数值模拟渤、黄、东海环流的最优定解条件，分析包括路径与来源在内的三海环流的一些典型特征，以及冬季黄海暖流分叉现象的季节内变率，以期为深入探究东中国海的环流结构提供基础。

1 数值方法与数据

1.1 FVCOM 模型简介与参数配置

本文基于有限体积海岸模型（Finite-Volume Coastal Ocean Model，FVCOM）对渤、黄、东海环流进行模拟。该模型由美国马萨诸塞大学和伍兹霍尔海洋研究所联合开发，采用有限体积法对Navior-Stokes（NS）原始积分方程进行离散，在水平方向采用无结构三角形网格，能够更好地贴合浅海及河口区域的复杂岸线，垂向上使用σ 随底坐标刻画海底地形起伏，此外，FVCOM 模型支持干-湿网格技术解决漫滩问题，以及内外模分离技术加速计算[22]。近年来，FVCOM 模型被广泛地应用于海洋学的研究[23-25]。

计算域北至渤海辽东湾，南至台湾东南海域，西侧为辽宁至福建岸线，东至138°E，包括朝鲜半岛与九州在内的日本部分岸线。最高空间分辨率约为200 m，在离岸方向逐渐变粗。3 条开边界从北至南依次位于日本海南部与日本南部外海、琉球群岛东南海域、台湾海峡中部与台湾东南海域（图1）。设定模型自2014 年冷启动，使用不同全球海洋数据产品作为定解条件（详见1.2 节），基于全球数据产品逐时给定开边界处温度、盐度、流速与水位，形成全球-东中国海的两层嵌套结构，嵌套的开边界节点分布详见图1。风场与降水等大气强迫数据源于欧洲中期预报中心的第五代再分析数据集（The fifth generation of ECMWF reanalysis，ERA5）[26]，采用同化数据设定渤、黄、东海的底摩擦系数[27]，其余配置参数如表1 所示。运行3 a 时间足够稳定之后输出2017 年1 月至4 月的逐小时数据。

表1 FVCOM 模型参数配置Table 1 FVCOM model configuration

图1 模型计算网格及观测站位分布Fig. 1 Model domain and the location of observation stations

1.2 数据与方法

本文目标之一是评估常用海洋数据产品在渤、黄、东海数值模拟中的适用性，包括混合坐标海洋模型（Hybrid Coordinate Ocean Model, HYCOM）再分析数据[28]、涡解高精度海洋模型（OGCM for the Earth Simulator, OFES）再分析数据[29]、简单海洋数据同化资料（Simple Ocean Data Assimilation, SODA）[30]、日本沿海海洋再分析数据集（Japan Coastal Ocean Predictability Experiments Reanalysis, JCOPE2）[31]、高分辨率模型（Estimating Circulation and Climate of the Ocean, ECCO2）再分析数据[32]、世界大洋数据集2018版（World Ocean Atlas 2018, WOA18）、数字环境模型数据（Generalized Digital Environmental Model,GDEM）[33]共计7 套不同的数据集。评估步骤主要包括模型对初始场的敏感性评估和模型对开边界的敏感性评估：①模型对初始场的敏感性评估。开边界采用0.08° HYCOM 再分析资料，温、盐初始场分别采用7 套不同数据集（ECCO2、GDEM、HYCOM、JCOPE2、OFES、SODA、WOA18），驱动FVCOM 模型，结合观测资料，对温、盐初始场的敏感性进行评估，确定最优的温盐初始场方案。②模型对开边界的敏感性评估。类似对温、盐初始条件的敏感性评估，基于上一步筛选的最优初始场，分别采用5 套不同数据产品（ECCO2、HYCOM、OFES、SODA、SODA_h）作为模型的开边界输入，评估渤、黄、东海FVCOM 模型对开边界数据的敏感性，并确定最优的配置方案（表2）。

表2 海洋数据产品Table 2 Oceanographic data

本研究温盐实测数据来自2016 年至2017 年国家基金委冬季渤、黄海共享航次，包含2017 年1月共114 个大面站位的温盐深剖面仪（Conductivity-Temperature-Depth profiler, CTD）数据和4 个海床基声学多普勒流速剖面仪（Acoustic Doppler Current Profiler，ADCP）流速数据。CTD 的型号为Seabird 911 plus，ADCP 为Teledyne Sentenl 系列，工作频率为300 kHz，采样间隔为10 min，垂向分辨率为2 m，各站位置见图1。

利用偏差（Error，σ）、均方根误差（Root Mean Square Error，ERMS）与方差（Variance，S）分析对比结果，它们的计算公式为：

式中： yobs为 实测值； ymodel为模型结果。

2 结果分析

2.1 最优温盐初始场

在第一阶段试验中，采用2017 年1 月CTD 温、盐数据，对不同初始场下的数值结果进行验证（使用与站位位置、观测时间最接近的网格节点数据），结果表明，温度场的水平分布在渤海及北黄海部分均呈现自SE 向NW 的降温趋势，与冬季黄海暖水流经渤海海峡北上的热输运路径吻合（图2a）。不同初始场驱动下的数值结果（图2b～图2h）与观测结果整体相差较小，再现了黄海到渤海的降温趋势，但在南黄海及东海存在显著差异，为此，本文选择了2 个典型区域，即海州湾外海A（图2 中黑色虚线方框）和青岛外海B（图2 中黑色虚线椭圆框），分析温度数值模拟效果。分析表明，在A 区域存在NW 向延伸的暖舌，以11 ℃等温线作为暖舌边界，其主体可越过34°N。在各初始场方案中，A 区域的暖舌范围及走向与观测存在不同程度的差异，仅SODA 中暖舌伸展至34°N，其他初始场方案的暖舌范围均相对观测偏南，且存在向北的偏差，并在35°N 出现相似的低温偏差。在B 区域，山东半岛南部环绕分布着西黄海沿岸流的低温水体，外海被高温的暖流水占据。各初始场方案在山东半岛南部的B 区域存在不同程度的低温偏差，其中SODA 与WOA18 的偏差最大，HYCOM 与OFES 的偏差最小。SODA 方案的数值结果在A 区域较好，但在B 区域大幅偏离观测结果，与之相比，HYCOM 与OFES 对A、B 区域的模拟较为合理，与观测结果吻合最好。

图2 渤、黄海深度平均温度与温度偏差分布（初始场）Fig. 2 Depth-averaged temperature and temperature deviation (initial fields) in the Bohai Sea and the Yellow Sea

与温度相比，盐度的数值模拟结果与观测存在更加明显的差异（图3b～图3h）。观测（图3a）表明大部分区域盐度高于31.5，仅海州湾附近的盐度较低。各初始场方案较好地重现了渤、黄海的盐度分布特征，但低盐区（盐度低于30.5）量值及范围存在不同程度的偏差，其中SODA 方案的低盐区从莱州湾延伸至长江口，范围最大，OFES 和GDEM 次之，ECCO2 的低盐区范围最小。北黄海的盐度观测值为32.2 且分布较均匀，SODA 方案得到偏淡海水（偏差为—0.7），其他初始场方案出现高盐偏差，在ECCO2 试验中平均偏差达到0.6，HYCOM、JCOPE2 与OFES 方案同样出现高盐区域，但整体与观测值相差不大。

分析温度与盐度大面平均结果与观测值偏差（图4）可知：在50 m 水层以上，数值温度与观测温度之间偏差较小；在50 m 水层以下，SODA 方案与观测之间仍匹配较好，其余方案的偏差分布趋势相近，但呈现出低温态势。在75 m 水层时，温度的数值偏差达到2 °C。且温度偏差的方差随深度出现“低高低”变化，最高值出现在40～60 m 水层，这表明站位分布对温度模拟效果有很大的影响，在60～80 m 水层中，温度偏差较大，方差却较小，说明底层水体温度模拟精度问题是普遍存在的，相反，50 m 以上的水层不仅温度偏差较小，而且离散程度也较小，是温度模拟最优的水层。除SODA 和ECCO2 方案，其他初始场方案中盐度结果与观测值之间的偏差较小，SODA 方案的盐度整体偏低，ECCO2 方案与之相反，但两者对盐度的数值结果均不理想。SODA 方案盐度偏差的方差分布出现与温度偏差相反的“高低高”变化，极小值同样出现在40～60 m 水层，这表明该水层温度的数值结果对位置不敏感。其余初始场方案的盐度偏差方差随水深变化较小，这表明不同位置的盐度偏差随深度变化的趋势相似。

图4 大面平均温度与盐度偏差（初始场）随水深变化Fig. 4 Deviation (initial fields) of the average temperature and salinity versus depth

温度、盐度偏差分布（图2、图3）表明，HYCOM 与OFES 方案温度模拟效果较优，各方案在海州湾附近均存在低盐现象问题，其中OFES 与ECCO2 方案的偏差较小，但在北黄海区域ECCO2 方案存在不可忽略的高盐偏差，OFES 的高盐偏差相对较小。SODA 温度偏差随水深变化与观测值吻合最好，但其平均盐度明显偏小。除了ECCO2 在不同水深的平均盐度偏大、SODA 偏小外，OFES等其他初始场方案的温度误差在上层水深（水深小于50 m）较小，在中下层（水深大于50 m）较大，平均盐度与观测值吻合程度高，偏差较小。从各方面比较来看，OFES 方案表现出优势，认为在温盐初始场敏感试验中，OFES 方案表现最佳。

2.2 最优开边界强迫

在初始场敏感性试验中，OFES 方案对FVCOM 的驱动效果最佳。开边界敏感性试验在OFES 初始场的基础上建立，使用5 套开边界强迫场（ECCO2、HYCOM、OFES、SODA、SODA_h）驱动模型，并输出结果进行观测数据比较，从而得出最佳的开边界方案。

从深度平均温度（图5）和盐度（图6）数值结果的大面分布可以看出，与初始场敏感性实验结果类似，各开边界方案的温度结果从北黄海到渤海呈带状分布并逐渐降低，表现出黄海暖流及其余脉入侵带来的温度影响，与观测值一致。除了SODA 与SODA_h 方案在A 区域的暖舌与观测值吻合较好外，其他3 个开边界方案的暖舌延伸幅度均小于观测值。在开边界敏感性试验中，同样在B 区域出现低温，其中OFES 低温偏差最大，HYCOM、ECCO2 次之，SODA 与SODA_h 方案的误差最小。

图5 温度偏差分布（开边界）Fig. 5 Temperature deviation of numerical results (open boundary conditions)

图6 盐度偏差分布（开边界）Fig. 6 Salinity deviation of numerical results (open boundary conditions)

OFES 方案在黄海西部等区域存在显著的高盐误差，与观测值对比，其盐度数值结果较差。其他开边界方案的盐度在海州湾出现低盐中心，在黄海中部出现高盐误差，其中各方案的低盐偏差相差不大，高盐误差最小的是SODA 方案。

从大面平均温、盐数据的垂向变化（图7）可以看出，在25 m 水层以上，温度偏差随水深变化的趋势与观测值呈现出良好的一致性且量值差异较小，在25 m 水层以下，下降趋势的变化与观测值相对一致，但其量值上的偏差逐渐增加，其中SODA_h 方案偏差最小，最大温差仅为0.5 °C，OFES与SODA 次之，差异最大的是HYCOM，最大温差可达2.3 °C。在不同水深，除OFES 外，基于各开边界方案的盐度与观测值偏差较小，这一小偏差从表层保持到底层，且最大盐差不超过0.5。与OFES 提供的温盐初始场一样，OFES 开边界方案的盐度结果存在不可忽视的高盐误差，数值结果不可靠。

图7 大面平均盐度与温度偏差（开边界）随水深变化Fig. 7 Deviation (open boundary conditions) of the average temperature and salinity versus depth

4 个站点的流速数值结果与ADCP 观测值的合计均方根误差（ERMS）结果（表3）显示，流速的ERMS较小，在不同开边界条件下的差异也较小，最大值仅为0.2 m/s。这表明数值模拟结果能够较好地重现真实的流场大小，并且流速大小对开边界条件不敏感。但是，流向与实测流向的差异较大，ERMS约为80°，并且在不同方案中均会出现，其中，SODA 的流向ERMS最小，为75.54°。

表3 不同开边界方案流速结果与观测值的均方根误差Table 3 Root mean square error (ERMS) of the simulated currents and against the mooring ADCP for open boundary condition experiments

SODA 与SODA_h 的深度平均温度分布相似，与实测值对比效果良好。OFES 在海州湾的盐度数值效果最好，低盐偏差较小，但在黄海中部的高盐误差不可忽略。相比之下，SODA 在黄海中部的高盐偏差最小，结果令人满意，同时，在海州湾的盐度数值效果仅次于OFES。在不同水深的对比中，SODA_h 方案的温度数值效果最好，SODA 在中上水层（水深小于 50 m）与SODA_h 相似，在下部水层（水深大于50 m）的模拟存在低温偏差，但在可接受范围之内，除ECCO2 外，其余开边界方案的盐度差异较小。流速大小的数值结果令人满意，但流向与实测值的出入较大，其中SODA 流向的均方根误差最小。在不同方面的对比中，SODA 开边界方案综合表现最好，所以，我们认为SODA 提供的开边界为最优开边界强迫场。

2.3 初始场再评估

在初始场敏感性试验中，采用的HYCOM 再分析数据作为开边界强迫，评估结果显示OFES 为最优温盐初始场。而基于OFES 的开边界强迫场敏感性试验对比得到的最优开边界为SODA。为了消除初始场敏感性试验中开边界选取带来的偶然性误差，使用SODA 数据作为开边界，开展了温盐初始场的再次评估，验证OFES 作为最优初始场的普遍性。

从各数值结果总体偏差的占比分布结果（图8）可以看出，不同初始场的温度偏差对敏感性不高，各数值结果表现相差不大。盐度偏差的差异明显，这与第一次评估时所得结果类似，其中SODA 的盐度偏差占比大于其余6 个盐度偏差占比之和，因此，在采用SODA 作为初始场时会引入不可忽略的误差。最小盐度偏差出现在OFES 作为初始场时，与次小的WOA18 相差0.41%，总偏差占比中，OEFS 偏差占比最小为19.91%，说明OFES 作为最优温盐初始场不具有偶然性，并且更换开边界后得到了与第一次温盐初始场评估一致的结论，这可佐证初始场敏感性试验结果。

图8 温、盐总偏差占比（初始场再评估）Fig. 8 Percentage of total temperature and salinity bias (reevaluation of initial fields)

2.4 流场分析

根据前文分析可知，当使用OFES 温盐初始场和SODA 开边界强迫时，FVCOM 在渤海、黄海和东海的数值结果最为可靠。在此基础上，进一步讨论东中国海环流分布。

冬季渤海、黄海和东海深度平均流矢量分布（图9）显示：黑潮属于北赤道流在吕宋岛以东的北上分支，而且呈现流速强、流幅宽、流量大的显著特征。黑潮从台湾东北部流入东海大陆架，遇到陆坡阻碍后，流矢量发生反气旋式偏转，在124°～125°E 转为NE 向，并沿大陆架边缘继续前进，主轴逐渐延伸，至（128°E, 31°N）附近时分为2 支，其中一支北上进入对马海峡，另一支分别在九州西南、南部发生弯曲，之后再次折回NE 向，沿日本东岸向北太平洋发展。在黑潮N 向、NE 向流动过程中，可以在其东侧观察到与主轴流动方向相反的黑潮逆流，并伴随有涡旋产生。模拟结果与表层黑潮实测轨迹高度吻合[34]。PN 断面（图10b）冬季平均净流量Q=21.481 Sv，基于SODA3.4.2 资料计算得到Q=22.192 Sv，两者相差较小，这一结果与许达等[35]基于CTD 计算的PN 断面流量（21.5 Sv）相差不大，并与魏艳州等[36]基于长期水文资料预测的PN 断面流量变化趋势相符，进一步验证了模拟结果的准确性。黑潮主轴流速（流速大于0.2 m/s）呈单核带状分布，整个流层有400～500 m 的厚度，表层流速达1.2 m/s，黑潮主要流核位于大陆坡折带附近，并随着水深的增加，有向东倾斜的趋势。

图9 冬季的深度平均流场Fig. 9 Depth-averaged current field in winter

图10 PN 断面流速分布与流量大小Fig. 10 Vertical distribution of current velocity and current vectors along the PN section

对马暖流发生在东海NE 部，流经九州以西海域向北发展，之后转向NE，并流经朝鲜海峡，进入日本海。从图8 可知，构成对马暖流的水体主要来自九州西南海域的黑潮分支、两支北上的台湾暖流分支和南下的西朝鲜沿岸流，其中，黑潮分支的贡献最大。

台湾暖流是发生在闽浙外海、长江口以南海域，自SW 向NE 的一支暖流。从图9 可以看出，台湾暖流在（123°E, 27°N）附近分为2 支。近岸分支沿123°E 经线北上，流向舟山群岛，至30°30′N的长江口外，发生连续弯曲，转向NE，随后汇入对马暖流。外海侧分支在28°N 附近发生气旋式转向，后紧贴黑潮主轴左侧向九州岛方向延伸，并在九州以西海域汇入对马暖流。

黄海暖流是一支自济州岛南部流入、沿黄海槽向NW 延伸的高温高盐海流[13]，基于水文观测数据，前人发现黄海暖流在34°～35°N 的主段部分存在2 个分支：西北支和北支[13-14,16,37-38]，该分叉现象在海表温度数据与卫星红外影像的分析支持[15,39-40]。从图9 可以看出，冬季平均黄海暖流在34°30′N 附近出现明显的分叉。其中，北支规模较大，沿黄海槽向北发展，后进入北黄海，西北支主体可达山东半岛南岸，在向岸发展过程中，不断有偏北的小分支脱离，这些小分支经过反气旋式的弯曲后，再次折回北支，随之向北发展。

由于表现出显著的射流性质，前人研究多采用高温或高盐水舌间接表征黄海暖流路径与变化特征[15-16]，但近期有研究指出黄海暖流路径与暖舌并不一致[41]，本文基于FVCOM 的流速结果，使用流速数值结果提出一个角度指标，直观分析黄海暖流分叉强度及其冬季季节内变率。

通过统计逐小时的流场分布，确定黄海暖流分叉位置基本位于在矩形框C（图11）内，采用（122°25′48″E, 34°43′12″N）～（124°E, 33°30′N）的延长线为分界[16]划分西北分支与北分支，以0.04 m/s 为特征流速，对大于该值的数据进行线性拟合，并定义2 个分支拟合线的夹角为黄海暖流的分叉角度（图11），随后使用36 h 滑动平均的流场计算分叉角度。

图11 黄海暖流分叉角度示意图Fig. 11 Schematic diagram of the bifurcation angle of the Yellow Sea warm current

由36 h 低通滤波之后分叉角度的小波分析结果（图12）可知，分叉角度表现出“中间高、两边低”的周期能量分布，冬季开始阶段，分叉角度多为2～4 d 的弱周期变化，且周期能量散乱分布，1 月过后，逐渐出现3 d、5～8 d 的强周期变化，从2 月中旬持续到3 月中旬，到3 月末，周期能量逐渐消减。黄海暖流具有显著的季节变化，冬季最强，春季开始消退，夏秋季节消失，从冬季出现到春季逐渐消减，可大致分为4 个时期：初冬的生成期、仲冬的强盛期、冬春交际的消退期、春季的消散期，分叉结构在黄海暖流强盛期时比较稳定，所以分叉角度在1 月的周期能量较低，到了冬春交际时期，黄海环流逐渐消退形成不稳定分叉结构，分叉角度出现强周期变化，到了春季，流动消散无法维持分叉结构，可认为分叉角度基本描述出黄海暖流在冬春交际前后的消退变化。

图12 小波实部系数分布与功率谱分析Fig. 12 Wavelet real part coefficient distribution and power spectrum

从分叉角度、风向和风速（矩形框C 内的平均风场）的功率谱结果（图12b～图12d）可以看到，分叉角度的前3 个主要周期（3 d、5 d 和8 d）与风场周期具有良好的对应关系，尤其是与风向匹配最好，可以认为，风场是影响黄海暖流分叉强度的主要因素之一，分叉是补偿性黄海暖流响应冬季风季节变化的结果[16]。同时注意到，分叉角度存在20 d 以上接近月变化的周期，但由于时间序列较短，无法准确分辨。

3 结 论

本文使用有限体积海岸数值模型（FVCOM），对多套海洋数据产品进行评估。通过与观测数据进行多角度的对比验证，筛选出最适合东中国海数值模型的温盐初始场和开边界方案，并具体分析了东中国海的典型环流特征，同时，本文也深入研究了冬季黄海分叉现象的变化周期，并得出如下结论。

1）利用直观的控制变量法，发现温盐初始场和开边界方案对数值结果的质量有重要影响。经过不同角度的对比分析，我们认为在东中国海区域，使用涡解高精度海洋模型再分析数据（OFES）作为温盐初始场、简单海洋数据同化资料（SODA）作为开边界是FVCOM 最优搭配方案。

2）黑潮的主轴流速具有400～500 m 的单核垂直结构，主要位于大陆坡折带，随深度增加而向东倾斜；黑潮北上分支是对马暖流中最大的水体贡献者，台湾暖流的2 个分支最终都汇合成为对马暖流的一部分，也是对马暖流的主要水源之一。

3）黄海暖流主流存在2 个分支，根据分叉角度指标的衡量发现，风场变化是造成流动分叉的主要因素。黄海暖流的冬季分叉现象在2 月出现3 d、5～8 d 的强周期变化，而在1 月与3 月时，周期性变化现象不明显。

本文采用数值模式对东中国海的典型流动进行数据评估和深入分析。最优定解方案的结论基于2017 年冬季观测数据，对FVCOM 在该时间段内选择温盐初始场和开边界具有指导意义，但对于其他时间段仅具有一定的参考价值。研究流动和黄海暖流分叉部分的结果有助于加深对东中国海环流的认识，对该海域的后续研究具有借鉴和参考意义。