2018, Vol. 34
2. 空军研究院战场环境研究所，北京 100085;
3. 95788部队气象台，四川 成都 610043;
4. 95968部队气象台，北京 100097;
5. 空军气象中心，北京 100085;
6. 32022部队气象台，湖北 武汉 430074
2. Institute of Battlefield Environment, Air Force Academy, Beijing 100085, China;
3. Meteorological Observatory, Army Unit 95788, Chengdu 610043, China;
4. Meteorological Observatory, Army Unit 95968, Beijing 100097, China;
5. Meteorological Center of Air Force, Beijing 100085, China;
6. Meteorological Observatory, Army Unit 32022, Wuhan 430074, China
台风是天气尺度上最显著的海气相互作用过程，通常伴随着大风、巨浪、海表面温度下降等现象，因此对台风的预报必须充分考虑海气间双向响应和反馈的影响[1]。目前，受限于海上观测资料的缺乏，中尺度海气浪耦合模式是研究台风的主要手段。国外在这方面的工作开展较早，美国海军于1997年发展了一个双向耦合模式COAMPS[2]，经过多年的发展完善，初步在西太平洋和地中海区域取得了良好的模拟效果。进入21世纪后，耦合模式在耦合形式、耦合技术和应用范围上得到了更全面的发展，由NCEP开发完成的大气、海洋、飓风预报系统HWRF[3]，极大地提高了飓风的强度和结构的预报水平。我国气象工作者也积极开展了海气耦合模式的研究。蒋小平等[4]将大气模式MM5和海洋模式POM进行耦合，对台风Krovanh进行了模拟；关皓等[5]建立了MM5-POM-WW3海气浪耦合模式，对南海典型台风个例进行了模拟研究；刘磊等[6]构建了基于MCT耦合器的WRF-POM并行区域海气耦合模式对台风“格美”进行了数值模拟试验。
海洋上层的海表混合层与上覆大气直接相互作用，其中包括了海风搅拌驱动的混合、不稳定浮力强迫、流切变不稳定、湍流平流、非局地混合等过程，而垂直混合是海洋内部混合的重要部分[7]，它决定了海洋的层结和海水上翻速率，集中反映了大尺度运动和小尺度过程之间的关系，由于模式分辨率的制约，需要对其进行参数化。根据混合系数的选取方式不同，垂直混合参数化方案可分为：常混合系数、依赖于Richardson数的简单混合方案和湍流闭合模型。湍流闭合模型的求解基于湍流特征长度闭合假设，因此混合系数受长度尺度变量的影响很大。
不同的垂直混合参数化方案直接影响海洋模式的预报结果并使得海洋对大气强迫的响应产生影响，因此，垂直混合参数化方案的研究一直是国内外众多学者关注的重点。Bryan[8]通过研究原始方程大洋一般环流模型的数值解对垂直扩散量值的敏感性问题，发现无论是经圈方向还是纬圈方向的积分翻转环流都对垂直扩散敏感；金向泽等[9]在以前20层模型的基础上设计了一个新的30层大洋模型，并进行了垂直混合过程参数化的敏感性试验，研究表明水平黏性的减弱在赤道大西洋会减小北大西洋深水的返回率；李阳春等[10]发现垂直混合方案GM90主要受到等密度面扩散系数和等密度面厚度扩散系数的影响。
Mellor[11]提出了多阶预报方程求解湍流场，通过计算出湍流混合长度尺度和速度尺度。Rodi[12]引入了湍能耗散率ε方程来确定湍流长度尺度。Umlauf[13]发现湍流长度尺度的选择存在约束条件，并基于此提出了通用尺度模型(Generic Length Scale, GLS)，可以推导出Mellor和Yamada的k-kl方案以及Rodi的k-ε方案。中尺度海气浪耦合模式提供了6种垂直混合参数化方案，结合前人经验以及预报海域地形，本文拟采用Mellor-Yamada2.5及GLS-ε两种垂直混合参数化方案，对西北太平洋海域进行了预报试验，然后比较两者的应用效果，对耦合模式的垂直混合参数化方案进行分析和研究。
2 模式和方案 2.1 模式介绍区域中尺度海气浪耦合模式以中尺度耦合器MCT为主导，以区域气象模式WRF、海洋环流模式ROMS和海浪预报模式SWAN作为分量构成。海洋模式ROMS[14]是三维非线性、自由表面和基于地形跟随坐标的斜压原始模式。作为新兴的海洋模式，ROMS可模拟不同尺度的运动，在国外已被广泛地应用于海洋生物、海洋地质以及海冰研究的诸多领域。模式选取能够更加适当地描述流场受地形影响的S(Stretched)坐标系，在垂直方向上采用非等比例分层的方式，相对于采用水深等比例分层的经典模式POM，ROMS能够使温跃层和底边界层具有更高的解析度[15]。此外，ROMS使用新的高阶水平压力梯度算法，相比POM可有效减少模式在地形变化大的区域的水平压力梯度计算误差累积。大气模式WRF[16]则是一个完全可压的非静力模式，控制方程组都写为通量形式。模式采用Arakawa C网格，这种网格有利于在具有高分辨率的模拟中提高准确性。目前模式已经实现了地形跟随高度坐标(欧拉高度坐标)方案和地形跟随静力气压坐标(欧拉质量坐标)方案两种方案以供选择，时间积分采用三阶Runge-Kutta方案。而海浪模式SWAN[17]是由荷兰Delft大学开发的第三代近岸海浪模式。通过给定流场、水深和边界条件，可用于陆架浅水区、海岸、湖泊、河口水域的风浪、涌浪及混合浪的计算，能得出模拟海域内的各种波要素，具有稳定性好、计算精度较高等特点。此外，SWAN模式相对于WW3模式受水深地形影响较小，同时其对浪高的模拟能力更好[18]。
构建中采用中尺度多模式耦合器MCT，通过“即插即拔”方式方便地连接各分量模式，同时通过对耦合模式各分量的积分控制函数中加入耦合控制函数并予以调用，实现了WRF、ROMS和SWAN三者间双向参数传递和物理量交换，以及对耦合模式积分运行的控制。耦合模式的组成和变量传递如图 1所示。
|
图 1 耦合模式机制示意图 |
本研究模式区域范围基本覆盖西北太平洋地区海域，为避免各子模式间变量传递的插值误差，三个子模式均采用分辨率为6 km的水平网格，其中WRF子模式垂直方向分50层，低层加密，ROMS子模式在垂直方向分为40层，上层海洋加密，并将最大水深截断为4 500 m，同时，取北、东、南三个边界为开边界，以确保海上预报区域可更好地考虑潮流、潮汐、大气等的强迫作用，而将西边界选择在东亚大陆，设为固边界，避免由于开边界的不确定性对该海域海洋特征预报效果的影响，海洋子模式采用的垂直混合参数化方案见表 1。耦合模式运行中各分量模式每隔1 h进行变量交换。
|
表 1 两种不同垂直混合参数化方案 |
本文针对表 1的两个方案分别利用耦合模式进行了相同时段(2014年9月7—10日和17—20日)为期8天的回报试验，起报时间均为12 UTC，预报时长为72 h。其中，2014年9月7—10日无台风过程，17—20日有一次台风过程，这两个时段分别可代表不同天气形势下的耦合模式预报效果。WRF子模式采用T799全球预报场驱动[19]，T799预报场的水平分辨率为0.5 °×0.5 °，垂直层次为27层，时间分辨率为6 h。ROMS子模式的驱动背景场来自总站西北太平洋海洋预报场[20]，其水平分辨率为0.25 °×0.25 °，垂直层次为35层，时间分辨率为12 h。SWAN子模式采用提前24 h热启动方式，所需的驱动风场提取自T799预报场。用于检验大气预报效果的FNL再分析资料[21]融合了大量的观测资料以及卫星反演资料，其水平分辨率为0.25 °×0.25 °，已被广泛应用于数值模式的诊断分析中。此外，台风系统路径和强度的观测资料来自中央气象台台风网[22]。
3 结果分析 3.1 海表气象要素为了更全面评估不同参数化方案对模式预报效果的影响，文中选用评价模式比较常见的均方根误差(RMSE)和皮尔逊相关系数(R)作为效果评估参数，针对9月7日12 UTC—10日12 UTC和17日12 UTC—20日12 UTC共8天的2 m温度、2 m湿度和10 m风u、v分量的12~72 h预报效果进行了逐12 h分析。
参数计算公式：
|
其中，n为观测次数，xi为观测值，
均方根误差给出统计区域内误差幅度的平均情况，通过图 2可看出，除9月9日P1发生溢出外，P1和P2的各项均方根误差在各预报时段的变化较一致。7—10日的2 m温度均方根误差相对较低，17日之后开始逐渐增大，同时可看出存在一个较明显的昼夜变化；2 m相对湿度的均方根误差维持在9%左右，在18日和19日达到了最大；10 m风u分量和v分量分别在19日和10日的60 h预报和72 h预报存在一个明显的异常增大，考虑到实际天气情况，可能是台风“凤凰”的影响所导致的。
|
图 2 大气要素预报结果均方根误差分布对比 a. 2 m温度；b. 2 m相对湿度；c. 10 m风u分量；d. 10 m风v分量。 |
为了进一步明确两者之间的差异，将P1的均方根误差减去P2得到了两者的差异分布(图 3)。在7—10日，P1的2 m温度和2 m相对湿度均方根误差基本上小于P2，并且偏差值随预报时效的增加而增大，最大达到了0.1 K和0.45%，而在18日和19日，P1则明显大于P2，最大偏差值分别高达0.2 K和0.55%。
|
图 3 P1与P2均方根误差差异分布 a. 2 m温度；b. 2 m相对湿度；c. 10 m风u分量；d. 10 m风v分量。 |
相关系数反映的是预报场和客观分析场变化趋势的相似程度，由图 4可看出，P1与P2的2 m温度相关系数的贴合程度最高，均保持在0.8以上，10 m风u分量次之，2 m相对湿度较差，而10 m风v分量存在明显的随预报时效快速下降。
|
图 4 大气要素预报结果相关系数分布对比 a. 2 m温度；b. 2 m相对湿度；c. 10 m风u分量；d. 10 m风v分量。 |
同样地，将P1的相关系数减去P2得到了两者的差异分布(图 5)。在7—10日，P1的2 m温度相关系数和2 m相对湿度相关系数基本上大于P2，并且偏差值随着预报时效的增加而增大，最大达到了0.011和0.027，而在18日和19日，P1则明显小于P2，最大偏差值分别为0.017和0.080；P1的10 m风u分量和v分量相关系数基本上要比P2大，在19日偏差达到最大值0.056和0.055。
|
图 5 P1与P2相关系数差异分布 a. 2 m温度；b. 2 m相对湿度；c. 10 m风u分量；d. 10 m风v分量。 |
综合比较P1和P2均方根误差和相关系数的差异，可初步得出结论，即在无台风天气的情况下，GLS-ε方案对2 m温度和2 m相对湿度的预报效果略好于MY2.5方案，而在台风天气的影响下，MY2.5方案略优于GLS-ε方案。
3.2 台风路径及强度分析耦合模式对台风的预报情况需要从模式输出结果中确定台风的路径和强度。对台风中心位置的判定往往是以海平面气压场最小值为依据，但海平面气压场最小值并不能完整地反映台风的结构，因此本文同时参考850 hPa高度上相对涡度分布和10 m高度最大风速中心，从而更合理地确定出台风中心位置。9月17—20日耦合模式72 h台风预报路径与观测路径的对比分布如图 6所示。红色点代表了台风位置的观测值，而黑色和蓝色分别代表了P1和P2的预报台风位置。台风“凤凰”9月17日12 UTC于菲律宾马尼拉以东约900 km洋面生成，随后以25~30 km/h的速度向西北方向移动，强度逐渐加强，于9月19日12 UTC擦过吕宋岛东北部沿海一带，并于9月21日02 UTC登陆台湾南部，强度随之降低，然后转向北偏西方向移动，最终于9月22日09 UTC二次登陆浙江东部沿海地区。由图 6a、6b、6c可看出，耦合试验初始场台风中心的定位与观测值偏差较大，导致预报路径偏差很大，而P1与P2之间路径的偏差较小，即对台风路径的预报受初始场的大尺度大气环流影响较大，而受海气相互作用的影响相对较小。由图 6d可知，9月20日12 UTC—23日12 UTC的路径误差最小，前36 h基本在30 km以下，而在后36 h则保持在80 km以内。
|
图 6 台风预报路径与观测路径对比 a. 9月17日12 UTC—20日12 UTC；b. 9月18日12 UTC—21日12 UTC；c. 9月19日12 UTC—22日12 UTC；d. 9月20日12 UTC—23日12 UTC。 |
台风预报强度与观测值的对比如图 7所示。耦合模式初始场的台风强度较观测值均偏弱，9月17日12 UTC—20日12 UTC和9月18日12 UTC—21日12 UTC两个预报时段内，前60 h的强度偏差与初始时刻基本保持一致，随后达到观测强度，预报强度变化趋势与观测值吻合较好；9月19日12 UTC—22日12 UTC的台风强度的发展变化与观测值最接近；9月20日12 UTC—23日12 UTC的台风强度的发展落后于观测值。四个预报时段内，P1的强度普遍要大于P2，考虑到P1与P2引起的海表面温度(SST)水平分布不同，从而对台风的发展造成影响，接下来将对台风经过前后的SST变化进行检验。
|
图 7 同图 6，但为台风中心最低气压(单位：hPa)对比 |
考虑到台风预报强度与观测值的吻合程度，以9月19日12 UTC—22日12 UTC这一时段为例，对台风经过前后的SST变化与NOAA提供的日平均最优插值(OI) SST[23]进行了比较，其空间分布如图 8所示(并给出了台风观测路径)。对台风强度造成影响的因素有很多，包括台风下垫面SST、风切变、眼壁置换周期以及天气尺度的大气强迫等[24]。根据试验的目的，两组预报试验的区别在于海洋垂直混合参数化方案的不同。在模式其他设置相同的情况下，两组试验预报的SST存在差异。从而造成台风强度不一致。
|
图 8 OI SST(a~c)、P1(d~f)和P2(g~i)试验的SST及各自72 h变化 a、d、g对应起报时刻，即9月19日12 UTC，b、e、h对应72 h，即9月22日12 UTC，c、f、i对应72 h变化。 |
从起报时刻(9月19日12 UTC)至24 h(9月20日12 UTC)，沿台风路径周围的SST较一致，台风强度变化不大，从24 h至48 h(9月21日12 UTC)，台风先后经过台湾西南侧和东侧较温暖的海面，台风得到了显著增强，随后直至72 h(9月22日12 UTC)，台风北上经过较冷海面，强度持续下降。
将模式所用初始场SST与OI SST进行比较，前者较后者在台湾岛东侧和吕宋岛南侧海域要高出0.8~1.0 ℃，这可能是19—22日预报台风较观测值更强的因素。从台风经过前后海表面温度的变化来看，P1(图 8f)、P2(图 8i)与OI SST(图 8c)的空间分布较一致，即沿台风路径的海水变冷，但P1和P2在台风路径的右侧造成的海水降温范围和程度更大，这是在其他时段内耦合模式对台风强度的预报偏弱的主要原因。P2较P1的降温更明显，对应48~72 h的台风强度，发现P2较P1的强度偏弱，这是海洋模式不同预报结果对大气的反馈所造成的。这在一定程度上说明MY2.5方案对海洋上层的混合较GLS-ε方案更强。
3.4 混合层温度混合层对海气相互作用的影响极大，海洋与大气的热力及动力交换过程主要发生在混合层。为了研究不同垂直混合方案对混合层的影响，本文选取了台风路径附近三个不同的点位S1(121.94 °E，25.54 °N)、S2(123.25 °E，26.13 °N)和S3(123.47 °E，27.05 °N)，并在图 8c中进行了标明。以P1为例，混合层温度随时间变化如图 9所示。台风经过之后，混合层温度下降明显，对应混合层深度分别由15 m、18 m和20 m迅速增大至80 m、60 m和50 m，即台风引起的强烈垂直混合使得混合层深度增大。
|
图 9 P1不同点位混合层温度随时间变化图 |
为了进一步比较P1与P2两组试验在混合层的影响差异，将P2减去P1得到的混合层温度随时间变化如图 10所示。从图中可看出，P1与P2预报温度的差异在混合层存在一个明显的分层，差异层深度与混合层深度基本一致。在差异层以上P1略高于P2，说明GLS-ε方案引起的垂直混合较弱，使得混合层上层海温较高，结合图 8中SST的变化，可解释P1的台风强度较P2更大，即GLS-ε方案可为台风的发展提供更多的热量；在差异层内P2明显大于P1，说明受MY2.5方案影响而产生的温度下降在该位置减弱；在差异层以下两者的差别较小，即在混合层以下这两个方案引起的垂直混合程度相近。另外从图 10中可明显看出，上述现象在台风经过后得以加强。综上所述，混合层内MY2.5方案引起的垂直混合较GLS-ε方案更强烈。
|
图 10 不同点位混合层温度差异随时间变化图 |
本文以区域中尺度海气浪耦合模式为基础，通过2组数值预报试验初步探讨了不同垂直混合参数化方案对大气要素场预报效果的影响，并且讨论了台风过程中的海气相互作用，得到以下几点结论。
(1) 在无台风天气的情况下，GLS-ε方案对2 m温度和2 m相对湿度的预报效果要好于MY2.5方案，而在台风天气的影响下，则是MY2.5方案优于GLS-ε方案。
(2) MY2.5方案和GLS-ε方案对台风路径的预报的影响差距极小，台风路径的预报效果主要取决于初始场台风中心的定位及台风强度，下一步将考虑采用Bogus方案对台风路径的预报加以改进。
(3) 混合层热含量对台风的发展至关重要，在海洋对台风的响应过程中，MY2.5方案所造成的垂直混合更强烈，导致混合层温度偏低，从而使得台风预报强度偏弱，因此，从灾害预警的角度考虑，采用GLS-ε方案能使得耦合模式预报强度偏弱的情况得以改善。
| [1] | 张进峰, 黄立文, 文元桥, 等. 区域大气-海浪-海洋耦合模式的建立[J]. 武汉理工大学学报(交通科学与工程版), 2011, 35(3): 587-590. DOI:10.3963/j.issn.1006-2823.2011.03.037 |
| [2] | GOLAZ J C, WANG S, DOYLE J D, et al. Coamps ®-les: model evaluation and analysis of second-and third-moment vertical velocity budgets[J]. Boundary-Layer Meteorology, 2005, 116(3): 487-517. DOI:10.1007/s10546-004-7300-5 |
| [3] | CAVALLO S M, TORN R D, SNYDER C, et al. Evaluation of the advanced hurricane wrf data assimilation system for the 2009 Atlantic hurricane season[J]. Mon Wea Rev, 2013, 141(2): 523-541. DOI:10.1175/MWR-D-12-00139.1 |
| [4] | 蒋小平, 刘春霞, 齐义泉. 利用一个海气耦合模式对台风Krovanh的模拟[J]. 大气科学, 2009, 33(1): 99-108. DOI:10.3878/j.issn.1006-9895.2009.01.09 |
| [5] | 关皓, 周林, 薛彦广, 等. 南海中尺度大气-海流-海浪耦合模式的建立及应用[J]. 热带气象学报, 2012, 28(2): 211-218. DOI:10.3969/j.issn.1004-4965.2012.02.008 |
| [6] | 刘磊, 费建芳, 林霄沛, 等. 海气相互作用对"格美"台风发展的影响研究[J]. 大气科学, 2011, 35(3): 444-456. DOI:10.3878/j.issn.1006-9895.2011.03.06 |
| [7] | 汪雷, 王彰贵, 凌铁军, 等. 海洋模式中垂直混合参数化方案介绍[J]. 海洋预报, 2014, 31(5): 93-104. |
| [8] | BRYAN F. Parameter sensitivity of primitive equation ocean general circulation models[J]. J Phy Ocean, 1987, 17(7): 970-985. DOI:10.1175/1520-0485(1987)017<0970:PSOPEO>2.0.CO;2 |
| [9] | 金向泽, 张学洪, 周天军. Fundamental framework and experiments of the third generation of IAP/LASG world oceangeneral circulation model[J]. 大气科学进展, 1999, 16(2): 197-215. |
| [10] | 李阳春, 徐永福, 赵亮, 等. 全球海洋模式中CFC-11吸收对次网格尺度混合参数化的敏感性[J]. 海洋学报, 2007, 29(3): 31-38. DOI:10.3321/j.issn:0253-4193.2007.03.004 |
| [11] | MELLOR G L. A Hierarchy of Turbulence Closure Models for Planetary Boundary Layers[J]. J Atmos Sci, 1974, 31(7): 1791-1806. DOI:10.1175/1520-0469(1974)031<1791:AHOTCM>2.0.CO;2 |
| [12] | RODI W. Examples of calculation methods for flow and mixing in stratified fluids[J]. J Geophy Res Ocean, 1987, 92(C5): 5305-5328. DOI:10.1029/JC092iC05p05305 |
| [13] | UMLAUF L. A generic length-scale equation for geophysical turbulence models[J]. J Mar Res, 2003, 61(2): 235-265. DOI:10.1357/002224003322005087 |
| [14] | 孙一妹, 费建芳, 程小平, 等. WRF_ROMS-1.2中尺度海气耦合模式简介[J]. 海洋预报, 2010, 27(2): 82-88. DOI:10.3969/j.issn.1003-0239.2010.02.014 |
| [15] | 林夏艳, 管玉平, 刘宇. 2000年秋季东沙冷涡的三维结构及其演化过程[J]. 热带海洋学报, 2013, 32(2): 55-65. |
| [16] | MICHALAKES J, CHEN S, DUDHIA J, et al. Development of a next-generation regional weather research and forecast model[M]// Developments In Teracomputing. 2001: 269-276. |
| [17] | SWAN T. SWAN Cycle Ⅲ version 41.01AB User Manual[z]. Delft University of Technology, 2015. |
| [18] | 梅婵娟, 赵栋梁, 史剑. 两种海浪模式对中国黄海海域浪高模拟能力的比较[J]. 海洋预报, 2008, 25(2): 92-98. DOI:10.3969/j.issn.1003-0239.2008.02.013 |
| [19] | YAMASHITA C, LIU H, CHU X. Gravity wave variations during the 2009 stratospheric sudden warming as revealed by ECMWF-T799 and observations[J]. Geophy Res Lett, 2010, 37(22): 333-345. |
| [20] | CHASSIGNET E P, HURLBURT H E, SMEDSTAD O M, et al. US GODAE: Global ocean prediction with the hybrid Coordinate Ocean Model (HYCOM)[J]. Oceanography, 2009, 22(2): 64-75. DOI:10.5670/oceanog |
| [21] | 李长青, 闰之辉, 王瀛, 等. NCEP-FNL与T213L31同化资料对比分析中的地形差异影响[J]. 气象, 2007, 33(1): 62-69. |
| [22] | http://typhoon.weather.com.cn/ |
| [23] | REYNOLDS R W, SMITH T M, LIU C, et al. Daily high-resolution-blended analyses for sea surface temperature[J]. J Clim, 2007, 20(22): 5473-5496. DOI:10.1175/2007JCLI1824.1 |
| [24] | ZAMBON J B, HE R, WARNER J C. Investigation of hurricane Ivan using the coupled ocean-atmosphere-wave-sediment transport (COAWST) model[J]. Ocean Dynamics, 2014, 64(11): 1535-1554. DOI:10.1007/s10236-014-0777-7 |