基于四維變分同化方法的南海中尺度渦后報實驗*

趙 軍 高 山,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、水下滑翔機、潛標等現場觀測的溫鹽剖面數據, 進行南海中尺度渦的實際預報實驗。