南海是中国三大边缘海中最大最深的一个, 同时也是西北太平洋最大的半封闭式边缘海。南海作为中尺度涡旋活跃地, 已经引起越来越多海洋学者的关注。目前南海涡旋的形成机制主要包括有: 黑潮环流脱落(Zhang et al, 2017)、风应力旋度(Wang et al, 2005)、地形风急流(Wang et al, 2008)、黑潮入侵引起的锋面不稳定(Wu et al, 2007)、黑潮流径变化(Nan et al, 2011)等等。此外, 中尺度涡在南海北部大尺度环流的调整中也起着重要作用。Xiu等(2010)使用ROMS (regional ocean modeling system)模式结合AVISO (archiving validation and interpolation of satellite oceanographic data)资料提供的融合海面高度异常数据, 统计分析了1993—2007年南海中尺度涡的特征。模式模拟的涡旋数量与卫星观测数相差不大, 涡旋的生命周期和涡旋大小以及涡旋的垂直范围和涡旋大小都存在线性关系。Lin等(2015)采用ROMS模式对南海中尺度涡的三维特征进行了统计分析, 并将涡旋的形成机制分为3种类型, 表面风应力输入、海流与海底地形相互作用以及黑潮入侵。虽然数值模拟可以很好地再现南海中尺度涡特征, 但是由于各种不确定性, 诸如垂直混合参数化、大气强迫、分辨率等等都会对模拟结果产生影响, 导致模拟的中尺度涡旋与实际观测很难在时间和空间上较好地对应。研究表明, 只有在模式中加入数据同化才能模拟到与实际观测相近的中尺度涡(Oke et al, 2013), 进而为中尺度涡预报提供必要的初始条件。

现阶段随着海洋观测技术的不断进步, 越来越多的观测调查在南海展开, 尤其是卫星海洋高度计数据提供了长期的中尺度涡观测资料。若能将这些观测资料和动力模型结合, 可以极大地提高模拟系统的准确性, 同时也可以为中尺度涡的预报提供更为精准的初始条件。高山等(2007)使用大气物理研究所的OVALS (ocean variational analysis system)三维变分同化系统, 结合POM(princeton ocean model)模式对中尺度涡进行了同化模拟实验。该系统将卫星高度计资料同化反演成为温、盐伪观测数据, 然后再次进行常规温盐同化, 得到温、盐分析场。将实验结果与观测结果比较表明, 加入高度计资料同化的模拟结果远远好于未使用同化的模拟结果, 说明高度计反演温盐场的三维变分同化方法对于模拟中尺度涡现象是非常有效的。Neveu等(2016)基于0.1°分辨率的ROMS模式, 利用四维变分同化方法将两组历史观测序列同化到加利福尼亚流系统(California current system, CCS), 发现两组实验中的代价函数都显著降低, 此外通过定量评估发现数据同化还可以改变CCS的环流系统。赵福等(2017)基于1/30°分辨率的ROMS模式, 利用集合最优插值(ensemble optimal interpolation, EnOI)同化方法对南海北部中尺度涡进行了同化模拟研究。模拟结果再现了2013年冬季发生在台湾岛西南海域的一对冷、暖中尺度涡的生成及传播过程。结果同样表明模式加入数据同化能够模拟出与观测数据相近的中尺度涡。Xie等(2020)将1993—2011年的南海海表高度异常资料通过EnOI方法同化到模式中, 发现可以使模拟的气旋涡(反气旋涡)的数量增加10.3%(13.6%), 同时还发现经过数据同化后, 中尺度涡旋的发生有非常明显的季节性变化。

尽管已经有了如此多的研究基础, 但是我国海洋观测技术与海洋再分析预报技术起步较晚, 中尺度涡旋的数值模拟与同化工作发展较慢。以往对中尺度涡的研究工作大部分属于方法技术的研究范畴, 使用四维变分方法的中尺度涡同化预报研究比较缺乏, 没有形成实用性的海洋中尺度模拟预报系统产品, 与国外相比差距较大, 因此, 目前迫切需要研制实用性的中尺度涡数值预报应用产品。本文的目的是通过四维变分数据同化得到最优化初始场, 从而为中尺度涡的预报提供一定的理论基础和可行性方案。

本文使用ROMS模式提供的四维变分同化模块, 针对我国南海重点海域的中尺度涡进行了同化模拟实验。在此基础上, 使用同化结果作为初始场, 对南海中尺度涡进行了系统的后报实验。本文第1部分主要介绍ROMS模式、同化方案与涡旋探测识别方法, 第2部分介绍同化结果的验证和中尺度涡旋的后报模拟实验, 第3部分为结果与讨论。

1 模式与方法 1.1 ROMS模式设置

ROMS是一个三维的、自由海面、静水压原始方程模型, 水平方向为正交曲线网格, 具有基于地形跟随的坐标系统(Shchepetkin et al, 2005)。ROMS包含多种参数化垂直混合方案和平流方案。在本次模拟中, 水平(垂直)方向的动量、温盐平流方案采用的是三阶迎风(四阶中心)平流方案, 垂向涡黏性参数化方案采用非线性KPP混合方案(K- profile parameterization) (Large et al, 1994)。近年来ROMS作为新兴的海洋模式被广泛应用于海洋生态系统、沉积物、海气相互作用等领域。

本文选取的模拟区域为105.1°—147.2°E, 2.3°S—27.2°N, 水平网格分辨率为0.1°, 垂直方向分为50层, 表层和底层的拉伸系数分别为6.0和0.1。模式的水深数据为美国国家地球物理数据中心提供的ETOPO2地形资料(https://rda.ucar.edu/datasets)。为了减小随底坐标系中陡峭地形引起的压力梯度力的计算误差, 采用Shapiro滤波器(Penven et al, 2008)对地形进行平滑处理。

为了得到匹配模式自身的平衡初始场, 模式用气候态数据进行了60 a的起转(spin-up)过程。其中气候态的侧边界采用SODA (simple ocean data assimilation)再分析数据的气候态月平均数据。起转过程的初始场选取SODA资料一月份的月平均数据, 海表条件均来自气候态NCEP(National Centers for Environmental Prediction)数据。此外, 模式在海表强迫中加入了海表温度(sea surface temperature, SST)和海表盐度(sea surface salinity, SSS)以校正海表的温盐通量状况, 这里SST和SSS采用美国海洋资料中心(U.S. National Oceanographic Data Center)提供的分辨率为1°的海洋气候态平均资料WOA09(World Ocean Atlas 2009)。

经过长期的运行调整, 其动力学热力学过程已经与模式相适应, 选取第60年1月份的结果作为模式模拟的初始场, 为了保证侧向开边界、表面强迫场之间数据的匹配, 模式中的大气动力、热通量和淡水通量强迫场以及侧边界的海表高度、温度、盐度、斜压水平流速场均来自NCEP数据, 模式模拟的时间范围为2003年1月1日—2013年1月30日。

1.2 ROMS模式验证

基于以上模式的设置, 我们首先对所研究区域的表层温盐场和流场情况进行了验证。如图 1所示, 将模式模拟的气候态(2008年1月1日—2012年12月31日输出的平均结果)SST和SSS与WOA13数据进行对比。可以看到南海SST从北向南逐渐增加(图 1a), 而且其等温线的分布基本上和南海等深线的走向是一致的, 呈东北西南向分布。其中北部陆架区温度较低, SST基本在22—24 ℃之间, 从北边的陆架向南温度逐渐增加至26—28 ℃之间。而在西北太平洋附近, 等温线的分布基本上与纬度平行且由北向南逐渐递增。南海的气候态SSS则呈现自东北吕宋海峡附近向西南递减的团带状分布模式, 而且南海SSS自东北向西南的梯度差别非常小, 基本维持在33.5—34.5。此外, 图 1c、1d中均可以看到西北太平洋在8°—16°N附近有一个明显的低盐度区, 这主要是由于西向的北赤道流导致的。

另外西太平洋水途径吕宋海峡进入南海, 导致在吕宋海峡有一个明显的高盐水舌入侵现象。Qu等(2000)认为虽然存在混合效应, 但是南海的水团分布中, 高盐水的北太平洋热带水(North Pacific tropical water, NPTW)在25等密度面(位势密度为1 025 kg/m³的面)上依然能被追踪到。图 2分别比较了ROMS和WOA13在25等密面上的盐度分布图, 发现25等密度面上的盐度分布自吕宋海峡开始向西南扩散至加里曼丹岛, 同时中国南部和越南东部沿海的盐度又比较小, 所以呈现出典型的舌状扩散, 与已有的研究吻合(Qu et al, 2000)。

图 3-4是ROMS气候态的表层流场在夏季(冬季)与OFES (OGCM for the earth simulator)数据的比较结果, 结果表明, ROMS模拟的南海表层环流结构与OFES数据的结果大致吻合。其中夏季受西南季风影响, 南海南部上层环流有一个明显的反气旋结构(图 3), 而冬季整个南海上层环流都是一个气旋式结构(图 4)。但是无论是夏季还是冬季, 在南海西边界沿着陆架的位置总是有一条南海西边界流存在, 其中冬季的强度更为强劲。而在西北太平洋也可以看到一条西向流动的北赤道流, 北赤道流到达西边界后在大约14°N的位置分叉, 然后向北的一支形成黑潮, 向南的一支形成棉兰老流, 与已有的研究类似。此外图中还可以看到黑潮流经吕宋海峡的时候有明显的入侵现象, 实际上黑潮入侵南海不仅对整个南海环流以及南海水体性质有影响, 对南海北部中尺度涡的产生也具有重大作用。

1.3 涡旋的识别追踪方法

本文用到的涡旋识别和追踪方法采用Dong等(2011)的涡旋探测方法, 该方法基于Nencioli等(2010)年提出的流场几何特征的基础, 主要由下面三个步骤完成:

(1) 涡旋中心位置速度最小, 记为涡心, 涡心东西(南北)两侧速度v(u)的数值符号相反, 大小随着距中心点的距离线性增加, 在涡心附近, 速度矢量的旋转方向必须一致。

(2) 涡旋的边界定义为围绕涡旋中心的最大闭合流函数。

(3) 涡旋的运动轨迹需要在确定好涡心之后, 通过比较连续时间的涡心位置进行算法追踪。

简单来说, 以某时刻检测出来的中尺度涡心周围150 km范围内, 寻找下一时刻距离最近的一个相同极性的涡旋。若能找到这样一个涡, 则视为前一个涡旋的后继涡旋, 且该涡旋不会成为其他涡的后继涡旋, 然后以此继续追踪下去; 若未找到合乎条件的涡, 则该涡的演变终止。该方法应用广泛(Liu et al, 2012; Lin et al, 2015; 赵军等, 2018; Sun et al, 2018), 相比Okubo-Weiss方法和Winding Angle方法, 拥有很高的探测成功率和较低的误报率(Nencioli et al, 2010)。

1.4 同化数据与同化方案

同化数据采用了来自法国国家空间局卫星海洋学存档数据中心(Archiving Validation and Interpolation of Satellite Oceanographic Data, AVISO, http://www.aviso.altimetry.fr/duacs/)提供的海面高度异常资料, 该资料是0.25°分辨率的规则网格点, 并且融合了Jason-1、Jason-2、Cryo Sat-2卫星的沿轨海面高度异常资料, 是目前研究海表高度异常使用的最广泛的资料。本文采用卫星高度计资料进行同化的时间范围是2012年7月1日—2013年1月30日, 本次同化中取的海表高度数据误差为2 cm, 观测误差协方差为0.000 4 m²。

数据同化是将确定性模式结果和观测数据融合在一起, 本文的后报模拟实验是在上述基础上获得一个对海洋状态重新估计的初始场, 此时的初始场应该比单纯的模式结果更加接近真实情况。现阶段卫星高度计资料以及现场观测的温盐剖面数据为我们提供了很多的中尺度涡观测资料, 如果能将其同化入海洋模式中, 将成为提高中尺度涡数值模拟精度的有效途径。

目前, ROMS模式主要有3种与其匹配的四维变分(four dimensional variational, 4D-Var)同化模块(Moore et al, 2011): 原始方程的增量强约束方法(I4D-Var), 对偶方程的基于物理空间统计分析系统(physical-space statistical analysis system, 4D-PSAS)以及对偶方程的基于代表的方法(dual formu- lation representer-based variant of 4D-Var, R4D-Var)。I4D-Var是在模型控制变量的全部空间中进行最佳循环估计的搜索, 其稳定性比较好, 缺点是需要的计算量很大。而4D-PSAS和R4D-Var的计算量小, 但只是在观测的模型状态向量空间才被搜索。本文采用的是I4D-Var四维变分资料同化方法, 该方法首先由Courtier等(1994)提出, 包含有三个重要的模块: 非线性模块(nonlinear forward model, NLM)、切线性模块(tangent linear model, TLM)和伴随模块。ROMS求代价函数最小化所采用的是基于Lanczos方法的共轭梯度下降算法(Fisher, 1998)。其中代价函数公式为

式中, J为代价函数, 下标b和o分别代表背景场和观测场; δx表示x的增量; R是观测场误差协方差矩阵; B际上是由初始场、表面强迫场、背景场和模式的误差协方差组成的对角矩阵; G是在观测位置采样的切线性模型; d是第一个猜测值与观察值之间的差异。

数据同化通常将整个周期分为多段子周期进行分析, 以便于TLM由NLM的增量进行准确的描述。假如目标是在同化窗口[t_j, t_j+1]内通过循环估计找到J最小值, 通常采用分段滚动式同化方案(Broquet et al, 2009), 具体的循环步骤如下:

(1) 在第j个循环中, 初始条件来自上次的同化循环, 记为x_tj^{ini, j}, 可以用来计算NLM同化窗内的后预报场从而为该循环提供背景场。

(2) I4D-Var在第j个循环, 使用外循环(N_outer)共轭梯度下降算法得到最佳增量的估计δx_tj^{N_outer, j}, 用x_tj^{ini, j} + δx_tj^{N_outer, j} = x_tj^{N_outer, j}得到最优分析初始条件的估计。

(3) 当同化窗内有可用观测值时, 用(2)中的结果进行NLM计算得到j+1时刻最优观测的模式轨迹x_tj^{N_outer, j}。

(4) 将(3)得到的观测模式最优估计作为j+1时刻的初始点: x_{tj + 1}^{ini, j + 1} = x_tj^{N_outer, j}, 循环往复进行。

2 结果与分析

利用1.4节中介绍的I4D-Var同化方法, 采用5 d的时间窗进行分段滚动式同化, 外循环设置为1次, 内循环设置为30次, 对南海中尺度涡旋进行了系统的数据同化实验。同化的时间范围为2012年7月1日—2013年1月30日, 并主要对包括南海在内的105°—123°E, 9°—26.24°N的中尺度涡旋进行了统计分析。

2.1 同化结果的验证

将相同时间范围与空间范围内数据同化(I4D-Var)结果、AVISO实际观测和控制实验(control run, CR)中的中尺度涡旋个数进行统计, 其结果如表 1所示。由于AVISO数据的分辨率为0.25°, 因此表 1中只统计涡旋半径大于60 km的涡旋个数。

由表 1的对比可以看到, 在AVISO观测的涡旋中冷暖涡的个数相当, 约1︰1, 该结果与全球范围内冷暖涡的比例接近(Chelton et al, 2011), 但是冷涡个数稍大于暖涡, 该现象在I4D-Var和CR中也都得到了较好地模拟。I4D-Var统计到的总涡旋个数与实际观测中的个数更为接近, 相差仅8%, 冷涡和暖涡分别比观测多了13.6%和12%。I4D-Var结果相比CR, 其涡旋个数更接近为真值。而未同化的结果则差别很大, CR中的涡旋个数偏多, 相差27%。找出整个同化时间段内每天AVISO和I4D-Var一一对应的涡旋, 对于半径在60 km以上的涡旋, I4D-Var的成功识别率为90%。然后分别计算了两者的涡心距离差距、半径相对误差和振幅相对误差, 如图 5所示。可以看到I4D-Var与AVISO一一对应的涡旋中, 其平均涡心距离有4.3海里的差距(图 5a), 平均振幅相对误差为0.16 (图 5c), 而平均半径相对误差很小, 仅有0.07 (图 5b), 上述结果充分说明经过数据同化的模拟场可以再现真实海洋中的中尺度涡旋, 而且与观测涡旋的位置、振幅、半径大小也都十分吻合。

随后将I4D-Var、AVISO及CR三者在同化时间段内涡旋出现的位置分布进行对比。由图 6b的卫星观测数据可以看到, 本次同化时间范围内中尺度涡旋主要出现在南海南部, 尤其是在越南沿岸东侧涡旋的出现频率最高, 其原因主要与季风驱动和地形作用有关(Chen et al, 2009)。其余沿着118°E也有两个出现频率相对集中的位置分别位于13°N和15°N附近。在该时间段内, 虽然南海北部的涡旋出现频率很低, 但是在吕宋海峡出入口的涡旋出现频率普遍较高, 其中台湾西南部是南海北部涡旋出现频率的高值区域。有关台湾岛西南部的暖涡生成机制仍然是许多学者关注的热点问题。而I4D-Var的涡旋出现频率无论是在强度还是空间分布与AVISO都有很好的对应关系, 台湾岛西南部涡旋的出现频率要略高于观测结果。而CR中涡旋的出现位置与观测相比相差比较大。由此可见, 数值模式经过数据同化后能够很好地产生与AVISO对应的中尺度涡, 而为加入数据同化的模拟中涡旋出现位置的分布则具有非常强的随机性。

接着将同化时间段内所有的涡旋半径进行了比较分析, 结果如图 7所示。结果表明, I4D-Var和CR中半径在60 km以下的涡旋分别占了总涡旋数量的34%和28%, 而AVISO数据识别的涡旋大部分都在60 km以上(92%), 造成该差异最主要因素可能是分辨率的不同, 如果只统计半径大于60 km的涡旋比例分布, 则I4D-Var与AVISO观测结果非常吻合, CR结果则相差很多(表 1)。

为进一步证明同化结果的有效性, 还将I4D-Var、AVISO及CR中涡旋的半径空间分布进行对比。从图 8中可以看到I4D-Var的涡旋半径分布(图 8a)与实际观测(图 8b)非常一致, 大半径的涡旋主要分布于南海南部115°E以西的地方, 且半径普遍要大于100 km。南海北部大半径涡旋较少, 出现频率较多的台湾岛西南部涡旋半径主要在80 km左右。而CR (图 8c)中无论是涡旋半径空间位置分布还是大小分布都与观测相差甚远。上述的统计结果间接的表明了I4D-Var相对于CR具有明显的优越性, 证明了在模式中加入四维变分同化对中尺度涡的模拟是非常有效的。

下面可通过2012年11月13日—2012年11月28日时间段的个例分析来直观的看一下加入同化的数值模拟效果。图 9是AVISO的海表高度异常(sea surface height anomalies, SSHA)和地转流分布图, 在11月13—28日期间南海海盆西部有一个越南冷涡(Vietnam cyclonic eddy, VCE)和一个越南暖涡(Vietnam anticyclonic eddy, VAE)组成的涡旋对, 该涡旋对通常在越南东岸周期性出现, 所以也将其统称为越南涡旋对(Vietnam eddy pair, VEP) (Chen et al, 2010)。Chen等(2010)研究表明越南涡旋对通常发生在夏秋季节, 局地风应力是该涡旋对年际性出现的主要原因。另外一个典型的涡旋就是南海东北部发展起来一个比较强的暖涡, 该暖涡位于台湾西南部, 所以也常被成为台湾暖涡(Taiwan anticyclonic eddy, TAE)。TAE经常发生在冬季, 有关TAE的生成机制有许多的研究, 其中Nan等(2011)提出了涡旋与黑潮相互作用模型, 认为黑潮路径的变化是该地区涡旋形成的主要诱因。冬季由于受东北季风影响, 黑潮环流较为频繁的出现在台湾东南部, 极易发生流套脱离现象, 从而发展为暖涡。

图 10是图 9在相同时间段由I4D-Var同化得到的SSHA和地转流分布图, 从图中可以看出, I4D-Var的SSHA和流场的空间分布非常接近AVISO的观测结果。图中各个涡旋的位置可以互相对应, 除了上述VEP和TAE都得到很好地模拟, 其余涡旋也均被很好地模拟, 如11月13—22日南海东部吕宋海峡口附近吕宋冷涡(Luzon Cyclonic Eddy, LCE, Wang et al, 2008)、11月13—19日南海中部113.6°E, 16°N附近的中部气旋涡(center cyclonic eddy, CCE), 11月13—28日南海东部118.2°E, 13°N黄岩岛附近有一个肉眼可识别的冷涡(Huangyan Dao cyclonic eddy, HCE)。有趣的是, 在I4D-Var和AVISO结果中均可以看到一个由LCE、CCE和VCE贯穿南海中部的冷涡序列结构, Chu等(2020)基于卫星观测和现场观测数据的统计结果指出南海冬季沿17°N和西边界形成的反“L”型涡列呈周期性出现。

为了进一步检验模式加入同化后对温、盐场的调整, 选取了某时刻沿115°E (10°—18°N, 水深 < 300 m)的温盐剖面图, 该经度恰好穿过一个冷涡和一个暖涡。图 11a、11b上方是该剖面上AVISO(黑线)、I4D-Var (红线)和CR(蓝线)的SSHA比较结果, 通过AVISO结果可以看到, 该剖面上SSHA分别在13°N和14.7°N, 出现极小值(0.02 m)和极大值(0.27 m), 对应着一个冷涡和一个暖涡。其中, I4D-Var的SSHA极小值点与AVISO基本一致, 极大值点位于15.2°N附近, 相比AVISO稍微向北偏移0.5°。而CR的SSHA普遍小于AVISO和I4D-Var, 在16.6°N有一个极大值点(0.12 m), 但是距离AVISO和I4D-Var的极大值距离比较远, 其结果不确切。经计算, CR的SSHA均方根误差(RMSE)为14%, 而I4D-Var的RMSE为5%, 相比CR得到了8%的优化。比较I4D-Var和CR的温盐剖面可以看到, 图 11a冷涡附近(蓝色阴影)的等温线有明显的向上弯曲, 涡旋内的水温比周围水体要小, 且该冷涡的温度中轴线向南偏移; 暖涡(红色阴影)内的等温线则向下弯曲, 表层的海水幅聚下沉, 其18℃等温线最深可达90 m的位置; 而在未同化的CR实验中(图 11b)两个阴影位置的等温线都比较平缓, 没有出现明显的弯曲情况。同样地, 在盐度剖面图中, I4D-Var中冷涡内(图 11c蓝色阴影)的盐度变大, 下层的高盐水有强烈的上升趋势; 暖涡内(图 11c红色阴影)的盐度等值线明显向下弯曲, 但是没有冷涡弯曲的强烈; 反观CR的盐度剖面图(图 11d)两个阴影内的等值线都较为平缓, 说明这里没有涡旋存在。

经过上述同化结果与AVISO观测的一系列比较表明, 相比单纯模式模拟的结果, 加入四维变分同化后的模式可以极大地提高对南海中尺度涡旋的模拟能力。在一个涡旋剖面个例的验证中发现同化后的模式相比不同化对SSHA模拟的优化可达8%, 而且模式加入同化后不仅可以改变海面高度, 其最终还会导致整个温度、盐度和密度场的改变。由此可以证明模式加入数据同化之后可以得到一个接近真实的海洋动力学状态, 在该基础上可以展开中尺度涡的后报模拟及时效性检验工作。

2.2 南海中尺度涡旋的后报模拟实验及其时效性检验

通过上面的分析可以看出, I4D-Var能够很好地模拟出与观测结果相符的中尺度涡, 但是该结果是在有观测资料同化的前提下。那么, 如果只有较为准确的初始场, 而没有同化的帮助, 仅仅靠海洋模式是否能够对中尺度涡旋的进一步运动进行准确地模拟呢？众所周知, 中尺度涡虽然称之为中尺度, 但实际上已属于大尺度海洋现象的范畴, 它的运动受到准地转平衡关系的约束。因此, 从理论上讲, 像风场、水热通量等外在的强迫场虽然对中尺度涡有长期的影响, 但是短期内对中尺度涡的运动和发展影响不大(高山等, 2007)。因此, 可以想象, 借助中尺度涡的这种地转平衡约束的惯性, 使用模式模拟在短期内是应当可以保证涡旋模拟的大致准确性。为了验证该假设的可行性, 我们决定使用同化模拟得到的结果作为初始场, 然后仅依靠ROMS海洋模式, 对中尺度涡进行后报实验, 以初步检验上述猜测的合理性, 进而确定中尺度涡的预报时效。具体方案如下:

(1) 同化进行到2012年12月8日时停止, 得到该时间段的同化场结果;

(2) 将(1) 中的结果作为初始场, 保持海表强迫场和侧边界场不变, 仅靠模式本身运行, 对中尺度涡进行后报模拟, 运行时间为3个月(2012年12月8日—2013年3月8日);

(3) 将后报模拟(hindcast control, HC)结果与相同时间段的AVISO高度计资料进行对比。

需要指出的是本文中所说的预报, 实际上是在同化结果提供初始场的基础上进行的后报模拟, 主要用于检测同化模拟的有效性。

分别选取HC在第5、10、15和20天的SSHA结果与观测结果进行对比, 结果如图 12所示。首先HC模拟的SSHA分布与AVISO十分吻合, 该时间段内HC模拟的南海中尺度涡与AVISO也有很好的对应关系, 其中TAE、HCE、VCE和CAE在后预报的20 d内均持续存在, VAE也明显地存在于前10 d, 之后逐渐消亡。分别计算后预报20 d内HC和AVISO之间对应涡旋的SSHA相关系数、涡心差距和半径绝对误差, 如图 13a所示, HCE、VCE和CAE在前20 d内的SSHA相关系数均在0.7以上。TAE在前13 dSSHA相关系数均在0.85以上, 但之后有减小趋势, 其可能原因是AVISO中TAE逐渐向西南方向移动(图 12e和12g), 而HC中虽有移动, 但很缓慢。VAE的SSHA相关系数在前10 d皆在0.6以上, 10 d以后由于VAE开始消亡, 相关系数迅速减小。另外, 图 13b可以看出, 除了VCE的涡心差距在大部分时间内保持一致外, 其余涡旋的涡心差距虽然均有增大的趋势, 但是皆在30海里范围内, 平均涡心差距仅为15海里。VCE、CAE和VAE由于涡旋形态的高时变性, 导致三者的半径绝对误差随时间略有增加, 但是基本维持在30海里以内, 而TAE和HCE较为稳定, 其平均半径绝对误差分别为3.5和6.6海里。由此可见, 上述后报结果验证了我们的猜测, 即由同化结果而来的初始场具备了中尺度涡巨大的准地转运动惯性, 仅使用ROMS海洋模式对南海中尺度涡旋进行预报模拟是可行的, 对于本次南海中尺度涡旋的后报模拟, 预报时效可达10 d以上。

为了更清晰地量化中尺度涡的预报模拟效果, 本研究进一步计算了所有中尺度涡(图 12中标注涡旋)运动学参数的平均相关系数, 并以此作为预报的评估指标。所选取的运动学参数主要包括: SSHA、相对涡度(relative vorticity, RV)、涡动能(eddy kinetic energy, EKE)和水平散度(divergence, Div), 相应的计算公式如下:

式中, V_R表示相对涡度; E_EK表示涡动能; D表示水平散度; 和分别代表地转速度异常的纬向分量和经向分量:

4个运动学参数的平均相关系数随时间变化如图 14所示, 可以看到在后预报的20 d中, 上述4个参数的平均相关系数均在0.67以上, 其中EKE和RV的相关性甚至皆在0.78以上, 该结果进一步证明四维变分同化技术可以提供合理的预报初始场, 能极大地提高中尺度涡的短期预报准确率。

通过图 12和图 13可以看到TAE相比模拟区域中的其他涡旋都要强大, 半径大小与观测高度吻合, 涡旋形态的时变性小, 并且持续地向西南方向移动, 下面以TAE作为研究个例, 继续检验经四维变分同化优化过的初始场对于中尺度涡旋的预报效果。继续探测并追踪后报模拟中涡旋的变化, 发现TAE自2012年12月8日开始继续运行了86 d, 一直到2013年3月2日终止。而在AVISO观测中该涡旋生命周期有89 d (2012年11月28日—2013年2月24日)与预报结果比较一致。两者运行路径的比较如图 15所示, 两者皆产生于南海东北部约118.5°E, 21°N附近, 然后向西南方向运动, 并于113.5°E, 17.5°N附近消亡。图 15可以看出HC模拟的移动路径和AVISO大致吻合, 说明同化提供的初始场对于预报中尺度涡位置的演变发展有很好的帮助。此外还比较了涡旋在HC模拟和AVISO观测两种情况下半径和振幅相对误差的时间变化, 如图 16所示, TAE振幅和半径的变化趋势几乎是同步的, 其中, 前30 d内TAE平均半径相对误差和平均振幅相对误差分别为-0.01和-0.13, 证明前30 d, TAE的后报模拟半径和振幅与观测高度吻合。30—65 d, HC模拟的涡旋经历了一个半径和振幅减小又增大的过程, 直至最后消亡的阶段, HC模拟的涡旋半径和振幅又恢复到和观测接近的状态。上述结果表明后预报场可以预报模拟个别强壮涡旋整个生命周期的移动路径, 并且在短时间内高度还原真实涡旋的大小和振幅变化, 其他运动学参数的预报模拟结果还需要进一步的验证。

综合上述对比结果进一步表明, 模式的后报结果的确在至少10 d内依然与AVISO观测分析结果非常吻合。个别成熟且稳定发展的涡旋, 比如TAE甚至能够预报整个生命周期的移动路径。后报模拟因为有了四维变分同化给出的合理初始场, 可以极大地提高中尺度涡旋模拟的准确性, 究其原因有如下两点:

(1) 得益于ROMS成熟的四维变分同化模块给后报模拟提供了一个接近真实海洋状态的初始场, 其海面高度和流场与实际观测结果非常接近, 而且其对应的温度、盐度和密度场均得到调整。

(2) 海洋中地转平衡的调整是一个缓慢的过程, 对于已经形成的中尺度涡来说, 其运动和发展主要受到地转平衡的约束, 而快速变化的强迫场在短时间内对其影响不大。

因此, 中尺度涡一旦形成, 则认为模式达到一个接近真实海洋地转平衡的状态, 在这种地转平衡的“惯性”驱动下, 模式可以很好地预测中尺度涡在短期内的运动状态。

3 结论

本文首先利用四维变分数据同化方法将AVISO卫星高度计资料同化到ROMS模式中, 并对同化结果进行了统计分析; 然后使用同化模拟得到的结果作为初始场, 仅依靠ROMS海洋模式, 对中尺度涡进行后预报实验, 探讨了经数据同化优化的初始场对于中尺度涡旋后预报模拟的可行性和有效性。结论如下:

加入四维变分数据同化之后的ROMS模式能够极大的提高南海中尺度涡的模拟精度。对比同化期间内卫星观测和同化结果的涡旋位置, 发现同化可以一一再现90%的南海中尺度涡旋, 其中平均涡心距离为4.3海里, 平均振幅相对误差为0.16, 平均半径相对误差仅为0.07。此外, 经同化再现的中尺度涡旋的数量、出现频率和半径的空间分布与观测都非常吻合, 而单独模式模拟的涡旋产生位置和大小则具有很大的随机性。

同化方法是通过对整个模拟场的海水状态进行调整产生的中尺度涡旋, 并使得涡旋的运动得以维持。通过对一个同时穿越冷涡和暖涡的剖面对比发现, 相比未加入同化的结果, 加入同化后SSHA的均方根误差得到了8%的优化, 而且同化模拟中冷、暖涡对应的温度和盐度场均产生相应的变化。

同化模拟可以帮助模式得到一个接近真实的中尺度涡模拟场, 以这个稳定的平衡场作为初始场, 进行正常的南海中尺度涡旋的后报模拟发现, 模式的后预报能力可达10 d以上。将个别稳定并强壮发展的涡旋的运动轨迹、半径及振幅信息与实际观测进行了比较, 结果表明其生命周期和运动轨迹几乎与观测吻合, 而且短时间内高度还原真实涡旋的大小和振幅变化。

值得注意的是同化得到的预报初始场实际上是一个中尺度涡旋的地转平衡状态, 所谓的预报, 则是中尺度涡在地转平衡“惯性”驱动下的结果。而诸如风暴潮、台风等极端天气都有可能打破上述地转平衡, 下一步可以对这些极端天气进行单独的讨论研究。本文中一个隐含的假设是模式、表面边界条件和侧向开边界条件是准确的, 但实际上, 四维变分同化中初始场、表面强迫场、背景场、模式和观测的误差协方差矩阵都会影响到代价函数的计算。因此, 下一步将采用国内外认可的实时强迫场数据继续调试优化模式本身的模拟能力, 并在同化中加入Argo、水下滑翔机、潜标等现场观测的温盐剖面数据, 进行南海中尺度涡的实际预报实验。