基于四维变分同化方法的南海中尺度涡后报实验*

赵 军 高 山,4① 王 凡,4

(1. 中国科学院海洋研究所 青岛 266071; 2. 中国科学院海洋环流与波动重点实验室 青岛 266071; 3. 中国科学院大学北京 100049; 4. 青岛海洋科学与技术试点国家实验室 海洋动力过程与气候功能实验室 青岛 266237)

南海是中国三大边缘海中最大最深的一个, 同时也是西北太平洋最大的半封闭式边缘海。南海作为中尺度涡旋活跃地, 已经引起越来越多海洋学者的关注。目前南海涡旋的形成机制主要包括有: 黑潮环流脱落(Zhanget al,2017)、风应力旋度(Wanget al,2005)、地形风急流(Wanget al,2008)、黑潮入侵引起的锋面不稳定(Wuet al, 2007)、黑潮流径变化(Nanet al, 2011)等等。此外, 中尺度涡在南海北部大尺度环流的调整中也起着重要作用。Xiu 等(2010)使用ROMS (regional ocean modeling system)模式结合AVISO (archiving validation and interpolation of satellite oceanographic data)资料提供的融合海面高度异常数据, 统计分析了1993—2007 年南海中尺度涡的特征。模式模拟的涡旋数量与卫星观测数相差不大,涡旋的生命周期和涡旋大小以及涡旋的垂直范围和涡旋大小都存在线性关系。Lin 等(2015)采用ROMS模式对南海中尺度涡的三维特征进行了统计分析,并将涡旋的形成机制分为3 种类型, 表面风应力输入、海流与海底地形相互作用以及黑潮入侵。虽然数值模拟可以很好地再现南海中尺度涡特征, 但是由于各种不确定性, 诸如垂直混合参数化、大气强迫、分辨率等等都会对模拟结果产生影响, 导致模拟的中尺度涡旋与实际观测很难在时间和空间上较好地对应。研究表明, 只有在模式中加入数据同化才能模拟到与实际观测相近的中尺度涡(Okeet 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 是一个三维的、自由海面、静水压原始方程模型, 水平方向为正交曲线网格, 具有基于地形跟随的坐标系统(Shchepetkinet al, 2005)。ROMS 包含多种参数化垂直混合方案和平流方案。在本次模拟中,水平(垂直)方向的动量、温盐平流方案采用的是三阶迎风(四阶中心)平流方案, 垂向涡黏性参数化方案采用非线性KPP 混合方案(K- profile parameterization)(Largeet 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 滤波器(Penvenet 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 °C之间, 从北边的陆架向南温度逐渐增加至26—28 °C 之间。而在西北太平洋附近, 等温线的分布基本上与纬度平行且由北向南逐渐递增。南海的气候态SSS 则呈现自东北吕宋海峡附近向西南递减的团带状分布模式, 而且南海SSS 自东北向西南的梯度差别非常小, 基本维持在33.5—34.5。此外, 图1c、1d 中均可以看到西北太平洋在8°—16°N 附近有一个明显的低盐度区, 这主要是由于西向的北赤道流导致的。

图1 长期平均的海表面温度(sea surface temperature, SST)和海表面盐度(sea surface salinity, SSS)Fig.1 Comparison of long-term mean SST (°C) from ROMS (a) and WOA13 (b) and SSS from ROMS (c) and WOA13 (d)

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

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

图3 夏季长期平均的表面流场比较Fig.3 Comparison of long-term mean surface velocity from ROMS (a) and OFES (b) in summer

图4 冬季长期平均的表面流场比较Fig.4 Comparison of long-term mean surface velocity from ROMS (a) and OFES (b) in winter

1.3 涡旋的识别追踪方法

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

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

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

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

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

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

目前, ROMS 模式主要有3 种与其匹配的四维变分(four dimensional variational, 4D-Var)同化模块(Mooreet 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 的增量进行准确的描述。假如目标是在同化窗口[tj,tj+1]内通过循环估计找到J最小值, 通常采用分段滚动式同化方案(Broquetet al, 2009), 具体的循环步骤如下:

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 数据同化(I4D-Var)结果、AVISO 实际观测和控制实验(CR)中涡旋个数统计Tab.1 Statistics of the number of vortices in I4D-Var, AVISO, and CR

由表1 的对比可以看到, 在AVISO 观测的涡旋中冷暖涡的个数相当, 约1︰1, 该结果与全球范围内冷暖涡的比例接近(Cheltonet 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), 上述结果充分说明经过数据同化的模拟场可以再现真实海洋中的中尺度涡旋, 而且与观测涡旋的位置、振幅、半径大小也都十分吻合。

图5 数据同化方法(I4D-Var)和AVISO 对应涡旋的涡心差距(a)、半径相对误差(b)和振幅相对误差(c)Fig.5 Eddy center distance (a), radius relative error (b), and amplitude relative error (c) of eddies corresponding to I4D-Var and AVISO

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

图6 I4D-Var (a)、AVISO(b)和控制实验(CR)(c)涡旋出现频率Fig.6 Frequency of eddies appearance in I4D-Var (a), AVISO (b), and CR (c)

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

图7 I4D-Var (a)、AVISO(b)和CR(c)涡旋半径比例图Fig.7 Scale diagram of eddies radius in I4D-Var (a) , AVISO (b) , and CR (c)

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

图8 I4D-Var (a)、AVISO(b)和CR(c)涡旋半径空间分布图Fig.8 Spatial distribution of eddies radius in I4D-Var (a), AVISO (b), and CR (c)

下面可通过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) (Chenet al,2010)。Chen 等(2010)研究表明越南涡旋对通常发生在夏秋季节, 局地风应力是该涡旋对年际性出现的主要原因。另外一个典型的涡旋就是南海东北部发展起来一个比较强的暖涡, 该暖涡位于台湾西南部, 所以也常被成为台湾暖涡(Taiwan anticyclonic eddy,TAE)。TAE 经常发生在冬季, 有关TAE 的生成机制有许多的研究, 其中Nan 等(2011)提出了涡旋与黑潮相互作用模型, 认为黑潮路径的变化是该地区涡旋形成的主要诱因。冬季由于受东北季风影响, 黑潮环流较为频繁的出现在台湾东南部, 极易发生流套脱离现象, 从而发展为暖涡。

图9 AVISO 海面高度异常SSHA(单位: cm)以及地转流(单位: cm/s)演变Fig.9 Evolution in SSHA (sea surface height anomalies) (unit: cm) and geostrophic current (unit: cm/s) of AVISO

图10 是图9 在相同时间段由I4D-Var 同化得到的SSHA 和地转流分布图, 从图中可以看出, I4D-Var的SSHA 和流场的空间分布非常接近AVISO 的观测结果。图中各个涡旋的位置可以互相对应, 除了上述VEP 和TAE 都得到很好地模拟, 其余涡旋也均被很好地模拟, 如11 月13—22 日南海东部吕宋海峡口附近吕宋冷涡(Luzon Cyclonic Eddy, LCE, Wanget 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”型涡列呈周期性出现。

图10 I4D-Var 海海面高度异常(单位: cm)以及地转流(单位: cm/s)演变Fig.10 Evolution in SSHA (unit: cm) and geostrophic current (unit: cm/s) of I4D-Var

为了进一步检验模式加入同化后对温、盐场的调整, 选取了某时刻沿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°C 等温线最深可达90 m 的位置; 而在未同化的CR 实验中(图11b)两个阴影位置的等温线都比较平缓, 没有出现明显的弯曲情况。同样地, 在盐度剖面图中,I4D-Var 中冷涡内(图11c 蓝色阴影)的盐度变大, 下层的高盐水有强烈的上升趋势; 暖涡内(图 11c 红色阴影)的盐度等值线明显向下弯曲, 但是没有冷涡弯曲的强烈; 反观CR 的盐度剖面图(图11d)两个阴影内的等值线都较为平缓, 说明这里没有涡旋存在。经过上述同化结果与AVISO 观测的一系列比较表明, 相比单纯模式模拟的结果, 加入四维变分同化后的模式可以极大地提高对南海中尺度涡旋的模拟能力。在一个涡旋剖面个例的验证中发现同化后的模式相比不同化对SSHA 模拟的优化可达8%, 而且模式加入同化后不仅可以改变海面高度, 其最终还会导致整个温度、盐度和密度场的改变。由此可以证明模式加入数据同化之后可以得到一个接近真实的海洋动力学状态, 在该基础上可以展开中尺度涡的后报模拟及时效性检验工作。

图11 沿115°E 的温盐剖面图Fig.11 Profiles of temperature and salinity along 115°E

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), 相应的计算公式如下:

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

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

图14 涡旋动力学参数平均相关系数的时间变化Fig.14 Time series of mean correlation coefficient of eddy kinematic parameters

通过图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 模拟的涡旋半径和振幅又恢复到和观测接近的状态。上述结果表明后预报场可以预报模拟个别强壮涡旋整个生命周期的移动路径, 并且在短时间内高度还原真实涡旋的大小和振幅变化, 其他运动学参数的预报模拟结果还需要进一步的验证。

图12 AVISO 和HC(hindcast control)的SSHA 对比Fig.12 Comparison in SSHA between AVISO and HC(hindcast control)

图13 HC 和AVISO 对应涡旋的SSHA 相关系数(a)、涡心差距(b)和半径绝对误差(c)的时间变化Fig.13 Time series of SSHA correlation coefficient (a), eddy center distance (b), and radius absolute errors (c) of eddies corresponding to HC and AVISO

图15 台湾暖涡运动轨迹图比较Fig.15 Comparison of TAE trajectory

图16 涡旋半径和振幅相对误差的时间变化Fig.16 Time series of relative error of eddy radius and amplitude

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

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

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

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

3 结论

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

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

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

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

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