趙 軍 高 山,4① 王 凡,4

(1. 中國科學(xué)院海洋研究所 青島 266071; 2. 中國科學(xué)院海洋環(huán)流與波動重點(diǎn)實(shí)驗(yàn)室 青島 266071; 3. 中國科學(xué)院大學(xué)北京 100049; 4. 青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室 海洋動力過程與氣候功能實(shí)驗(yàn)室 青島 266237)

南海是中國三大邊緣海中最大最深的一個, 同時也是西北太平洋最大的半封閉式邊緣海。南海作為中尺度渦旋活躍地, 已經(jīng)引起越來越多海洋學(xué)者的關(guān)注。目前南海渦旋的形成機(jī)制主要包括有: 黑潮環(huán)流脫落(Zhanget al,2017)、風(fēng)應(yīng)力旋度(Wanget al,2005)、地形風(fēng)急流(Wanget al,2008)、黑潮入侵引起的鋒面不穩(wěn)定(Wuet al, 2007)、黑潮流徑變化(Nanet al, 2011)等等。此外, 中尺度渦在南海北部大尺度環(huán)流的調(diào)整中也起著重要作用。Xiu 等(2010)使用ROMS (regional ocean modeling system)模式結(jié)合AVISO (archiving validation and interpolation of satellite oceanographic data)資料提供的融合海面高度異常數(shù)據(jù), 統(tǒng)計分析了1993—2007 年南海中尺度渦的特征。模式模擬的渦旋數(shù)量與衛(wèi)星觀測數(shù)相差不大,渦旋的生命周期和渦旋大小以及渦旋的垂直范圍和渦旋大小都存在線性關(guān)系。Lin 等(2015)采用ROMS模式對南海中尺度渦的三維特征進(jìn)行了統(tǒng)計分析,并將渦旋的形成機(jī)制分為3 種類型, 表面風(fēng)應(yīng)力輸入、海流與海底地形相互作用以及黑潮入侵。雖然數(shù)值模擬可以很好地再現(xiàn)南海中尺度渦特征, 但是由于各種不確定性, 諸如垂直混合參數(shù)化、大氣強(qiáng)迫、分辨率等等都會對模擬結(jié)果產(chǎn)生影響, 導(dǎo)致模擬的中尺度渦旋與實(shí)際觀測很難在時間和空間上較好地對應(yīng)。研究表明, 只有在模式中加入數(shù)據(jù)同化才能模擬到與實(shí)際觀測相近的中尺度渦(Okeet al, 2013), 進(jìn)而為中尺度渦預(yù)報提供必要的初始條件。

現(xiàn)階段隨著海洋觀測技術(shù)的不斷進(jìn)步, 越來越多的觀測調(diào)查在南海展開, 尤其是衛(wèi)星海洋高度計數(shù)據(jù)提供了長期的中尺度渦觀測資料。若能將這些觀測資料和動力模型結(jié)合, 可以極大地提高模擬系統(tǒng)的準(zhǔn)確性, 同時也可以為中尺度渦的預(yù)報提供更為精準(zhǔn)的初始條件。高山等(2007)使用大氣物理研究所的OVALS (ocean variational analysis system)三維變分同化系統(tǒng), 結(jié)合POM(princeton ocean model)模式對中尺度渦進(jìn)行了同化模擬實(shí)驗(yàn)。該系統(tǒng)將衛(wèi)星高度計資料同化反演成為溫、鹽偽觀測數(shù)據(jù), 然后再次進(jìn)行常規(guī)溫鹽同化, 得到溫、鹽分析場。將實(shí)驗(yàn)結(jié)果與觀測結(jié)果比較表明, 加入高度計資料同化的模擬結(jié)果遠(yuǎn)遠(yuǎn)好于未使用同化的模擬結(jié)果, 說明高度計反演溫鹽場的三維變分同化方法對于模擬中尺度渦現(xiàn)象是非常有效的。Neveu 等(2016)基于 0.1°分辨率的ROMS 模式, 利用四維變分同化方法將兩組歷史觀測序列同化到加利福尼亞流系統(tǒng)(California current system, CCS), 發(fā)現(xiàn)兩組實(shí)驗(yàn)中的代價函數(shù)都顯著降低, 此外通過定量評估發(fā)現(xiàn)數(shù)據(jù)同化還可以改變CCS 的環(huán)流系統(tǒng)。趙福等(2017)基于1/30°分辨率的ROMS 模式, 利用集合最優(yōu)插值(ensemble optimal interpolation, EnOI)同化方法對南海北部中尺度渦進(jìn)行了同化模擬研究。模擬結(jié)果再現(xiàn)了2013 年冬季發(fā)生在臺灣島西南海域的一對冷、暖中尺度渦的生成及傳播過程。結(jié)果同樣表明模式加入數(shù)據(jù)同化能夠模擬出與觀測數(shù)據(jù)相近的中尺度渦。Xie 等(2020)將1993—2011 年的南海海表高度異常資料通過EnOI 方法同化到模式中, 發(fā)現(xiàn)可以使模擬的氣旋渦(反氣旋渦)的數(shù)量增加10.3%(13.6%), 同時還發(fā)現(xiàn)經(jīng)過數(shù)據(jù)同化后, 中尺度渦旋的發(fā)生有非常明顯的季節(jié)性變化。

盡管已經(jīng)有了如此多的研究基礎(chǔ), 但是我國海洋觀測技術(shù)與海洋再分析預(yù)報技術(shù)起步較晚, 中尺度渦旋的數(shù)值模擬與同化工作發(fā)展較慢。以往對中尺度渦的研究工作大部分屬于方法技術(shù)的研究范疇,使用四維變分方法的中尺度渦同化預(yù)報研究比較缺乏, 沒有形成實(shí)用性的海洋中尺度模擬預(yù)報系統(tǒng)產(chǎn)品, 與國外相比差距較大, 因此, 目前迫切需要研制實(shí)用性的中尺度渦數(shù)值預(yù)報應(yīng)用產(chǎn)品。本文的目的是通過四維變分?jǐn)?shù)據(jù)同化得到最優(yōu)化初始場, 從而為中尺度渦的預(yù)報提供一定的理論基礎(chǔ)和可行性方案。

本文使用ROMS 模式提供的四維變分同化模塊,針對我國南海重點(diǎn)海域的中尺度渦進(jìn)行了同化模擬實(shí)驗(yàn)。在此基礎(chǔ)上, 使用同化結(jié)果作為初始場, 對南海中尺度渦進(jìn)行了系統(tǒng)的后報實(shí)驗(yàn)。本文第1 部分主要介紹ROMS 模式、同化方案與渦旋探測識別方法,第2 部分介紹同化結(jié)果的驗(yàn)證和中尺度渦旋的后報模擬實(shí)驗(yàn), 第3 部分為結(jié)果與討論。

1 模式與方法

1.1 ROMS 模式設(shè)置

ROMS 是一個三維的、自由海面、靜水壓原始方程模型, 水平方向?yàn)檎磺€網(wǎng)格, 具有基于地形跟隨的坐標(biāo)系統(tǒng)(Shchepetkinet al, 2005)。ROMS 包含多種參數(shù)化垂直混合方案和平流方案。在本次模擬中,水平(垂直)方向的動量、溫鹽平流方案采用的是三階迎風(fēng)(四階中心)平流方案, 垂向渦黏性參數(shù)化方案采用非線性KPP 混合方案(K- profile parameterization)(Largeet al, 1994)。近年來ROMS 作為新興的海洋模式被廣泛應(yīng)用于海洋生態(tài)系統(tǒng)、沉積物、海氣相互作用等領(lǐng)域。

本文選取的模擬區(qū)域?yàn)?05.1°—147.2°E, 2.3°S—27.2°N, 水平網(wǎng)格分辨率為0.1°, 垂直方向分為50層, 表層和底層的拉伸系數(shù)分別為6.0 和0.1。模式的水深數(shù)據(jù)為美國國家地球物理數(shù)據(jù)中心提供的ETOPO2 地形資料(https://rda.ucar.edu/datasets)。為了減小隨底坐標(biāo)系中陡峭地形引起的壓力梯度力的計算誤差, 采用Shapiro 濾波器(Penvenet al, 2008)對地形進(jìn)行平滑處理。

為了得到匹配模式自身的平衡初始場, 模式用氣候態(tài)數(shù)據(jù)進(jìn)行了60 a 的起轉(zhuǎn)(spin-up)過程。其中氣候態(tài)的側(cè)邊界采用SODA (simple ocean data assimilation)再分析數(shù)據(jù)的氣候態(tài)月平均數(shù)據(jù)。起轉(zhuǎn)過程的初始場選取SODA 資料一月份的月平均數(shù)據(jù), 海表?xiàng)l件均來自氣候態(tài)NCEP(National Centers for Environmental Prediction)數(shù)據(jù)。此外, 模式在海表強(qiáng)迫中加入了海表溫度(sea surface temperature, SST)和海表鹽度(sea surface salinity, SSS)以校正海表的溫鹽通量狀況, 這里SST 和SSS 采用美國海洋資料中心(U.S. National Oceanographic Data Center)提供的分辨率為1°的海洋氣候態(tài)平均資料WOA09(World Ocean Atlas 2009)。

經(jīng)過長期的運(yùn)行調(diào)整, 其動力學(xué)熱力學(xué)過程已經(jīng)與模式相適應(yīng), 選取第60 年1 月份的結(jié)果作為模式模擬的初始場, 為了保證側(cè)向開邊界、表面強(qiáng)迫場之間數(shù)據(jù)的匹配, 模式中的大氣動力、熱通量和淡水通量強(qiáng)迫場以及側(cè)邊界的海表高度、溫度、鹽度、斜壓水平流速場均來自NCEP 數(shù)據(jù), 模式模擬的時間范圍為2003 年1 月1 日—2013 年1 月30 日。

1.2 ROMS 模式驗(yàn)證

基于以上模式的設(shè)置, 我們首先對所研究區(qū)域的表層溫鹽場和流場情況進(jìn)行了驗(yàn)證。如圖1 所示,將模式模擬的氣候態(tài)(2008 年1 月1 日—2012 年12月31 日輸出的平均結(jié)果)SST 和SSS 與WOA13 數(shù)據(jù)進(jìn)行對比。可以看到南海SST 從北向南逐漸增加(圖1a), 而且其等溫線的分布基本上和南海等深線的走向是一致的, 呈東北西南向分布。其中北部陸架區(qū)溫度較低, SST 基本在22—24 °C之間, 從北邊的陸架向南溫度逐漸增加至26—28 °C 之間。而在西北太平洋附近, 等溫線的分布基本上與緯度平行且由北向南逐漸遞增。南海的氣候態(tài)SSS 則呈現(xiàn)自東北呂宋海峽附近向西南遞減的團(tuán)帶狀分布模式, 而且南海SSS 自東北向西南的梯度差別非常小, 基本維持在33.5—34.5。此外, 圖1c、1d 中均可以看到西北太平洋在8°—16°N 附近有一個明顯的低鹽度區(qū), 這主要是由于西向的北赤道流導(dǎo)致的。

圖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)

另外西太平洋水途徑呂宋海峽進(jìn)入南海, 導(dǎo)致在呂宋海峽有一個明顯的高鹽水舌入侵現(xiàn)象。Qu 等(2000)認(rèn)為雖然存在混合效應(yīng), 但是南海的水團(tuán)分布中, 高鹽水的北太平洋熱帶水(North Pacific tropical water, NPTW)在25 等密度面(位勢密度為1 025 kg/m3的面)上依然能被追蹤到。圖2 分別比較了ROMS 和WOA13 在25 等密面上的鹽度分布圖, 發(fā)現(xiàn)25 等密度面上的鹽度分布自呂宋海峽開始向西南擴(kuò)散至加里曼丹島, 同時中國南部和越南東部沿海的鹽度又比較小, 所以呈現(xiàn)出典型的舌狀擴(kuò)散, 與已有的研究吻合(Quet al, 2000)。

圖3—4 是ROMS 氣候態(tài)的表層流場在夏季(冬季)與OFES (OGCM for the earth simulator)數(shù)據(jù)的比較結(jié)果, 結(jié)果表明, ROMS 模擬的南海表層環(huán)流結(jié)構(gòu)與OFES 數(shù)據(jù)的結(jié)果大致吻合。其中夏季受西南季風(fēng)影響, 南海南部上層環(huán)流有一個明顯的反氣旋結(jié)構(gòu)(圖3), 而冬季整個南海上層環(huán)流都是一個氣旋式結(jié)構(gòu)(圖4)。但是無論是夏季還是冬季, 在南海西邊界沿著陸架的位置總是有一條南海西邊界流存在, 其中冬季的強(qiáng)度更為強(qiáng)勁。而在西北太平洋也可以看到一條西向流動的北赤道流, 北赤道流到達(dá)西邊界后在大約14°N 的位置分叉, 然后向北的一支形成黑潮,向南的一支形成棉蘭老流, 與已有的研究類似。此外圖中還可以看到黑潮流經(jīng)呂宋海峽的時候有明顯的入侵現(xiàn)象, 實(shí)際上黑潮入侵南海不僅對整個南海環(huán)流以及南海水體性質(zhì)有影響, 對南海北部中尺度渦的產(chǎ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)年提出的流場幾何特征的基礎(chǔ), 主要由下面三個步驟完成:

(1) 渦旋中心位置速度最小, 記為渦心, 渦心東西(南北)兩側(cè)速度v(u)的數(shù)值符號相反, 大小隨著距中心點(diǎn)的距離線性增加, 在渦心附近, 速度矢量的旋轉(zhuǎn)方向必須一致。

(2) 渦旋的邊界定義為圍繞渦旋中心的最大閉合流函數(shù)。

(3) 渦旋的運(yùn)動軌跡需要在確定好渦心之后, 通過比較連續(xù)時間的渦心位置進(jìn)行算法追蹤。

簡單來說, 以某時刻檢測出來的中尺度渦心周圍150 km 范圍內(nèi), 尋找下一時刻距離最近的一個相同極性的渦旋。若能找到這樣一個渦, 則視為前一個渦旋的后繼渦旋, 且該渦旋不會成為其他渦的后繼渦旋, 然后以此繼續(xù)追蹤下去; 若未找到合乎條件的渦, 則該渦的演變終止。該方法應(yīng)用廣泛(Liuet al,2012; Linet al, 2015; 趙軍等, 2018; Sunet al, 2018),相比Okubo-Weiss 方法和Winding Angle 方法, 擁有很高的探測成功率和較低的誤報率(Nencioliet al,2010)。

1.4 同化數(shù)據(jù)與同化方案

同化數(shù)據(jù)采用了來自法國國家空間局衛(wèi)星海洋學(xué)存檔數(shù)據(jù)中心(Archiving Validation and Interpolation of Satellite Oceanographic Data, AVISO,http://www.aviso.altimetry.fr/duacs/)提供的海面高度異常資料, 該資料是0.25°分辨率的規(guī)則網(wǎng)格點(diǎn), 并且融合了Jason-1、Jason-2、Cryo Sat-2 衛(wèi)星的沿軌海面高度異常資料, 是目前研究海表高度異常使用的最廣泛的資料。本文采用衛(wèi)星高度計資料進(jìn)行同化的時間范圍是2012 年7 月1 日—2013 年1 月30 日, 本次同化中取的海表高度數(shù)據(jù)誤差為2 cm, 觀測誤差協(xié)方差為0.000 4 m2。

數(shù)據(jù)同化是將確定性模式結(jié)果和觀測數(shù)據(jù)融合在一起, 本文的后報模擬實(shí)驗(yàn)是在上述基礎(chǔ)上獲得一個對海洋狀態(tài)重新估計的初始場, 此時的初始場應(yīng)該比單純的模式結(jié)果更加接近真實(shí)情況?，F(xiàn)階段衛(wèi)星高度計資料以及現(xiàn)場觀測的溫鹽剖面數(shù)據(jù)為我們提供了很多的中尺度渦觀測資料, 如果能將其同化入海洋模式中, 將成為提高中尺度渦數(shù)值模擬精度的有效途徑。

目前, ROMS 模式主要有3 種與其匹配的四維變分(four dimensional variational, 4D-Var)同化模塊(Mooreet al, 2011): 原始方程的增量強(qiáng)約束方法(I4D-Var), 對偶方程的基于物理空間統(tǒng)計分析系統(tǒng)(physical-space statistical analysis system, 4D-PSAS)以及對偶方程的基于代表的方法(dual formu- lation representer-based variant of 4D-Var, R4D-Var)。I4D-Var是在模型控制變量的全部空間中進(jìn)行最佳循環(huán)估計的搜索, 其穩(wěn)定性比較好, 缺點(diǎn)是需要的計算量很大。而4D-PSAS 和R4D-Var 的計算量小, 但只是在觀測的模型狀態(tài)向量空間才被搜索。本文采用的是I4D-Var 四維變分資料同化方法, 該方法首先由Courtier 等(1994)提出, 包含有三個重要的模塊: 非線性模塊(nonlinear forward model, NLM)、切線性模塊(tangent linear model, TLM)和伴隨模塊。ROMS 求代價函數(shù)最小化所采用的是基于Lanczos 方法的共軛梯度下降算法(Fisher, 1998)。其中代價函數(shù)公式為

式中,J為代價函數(shù), 下標(biāo)b 和o 分別代表背景場和觀測場;δ x表示x的增量;R是觀測場誤差協(xié)方差矩陣;B際上是由初始場、表面強(qiáng)迫場、背景場和模式的誤差協(xié)方差組成的對角矩陣;G是在觀測位置采樣的切線性模型;d是第一個猜測值與觀察值之間的差異。

數(shù)據(jù)同化通常將整個周期分為多段子周期進(jìn)行分析, 以便于TLM 由NLM 的增量進(jìn)行準(zhǔn)確的描述。假如目標(biāo)是在同化窗口[tj,tj+1]內(nèi)通過循環(huán)估計找到J最小值, 通常采用分段滾動式同化方案(Broquetet al, 2009), 具體的循環(huán)步驟如下:

2 結(jié)果與分析

利用1.4 節(jié)中介紹的I4D-Var 同化方法, 采用5 d的時間窗進(jìn)行分段滾動式同化, 外循環(huán)設(shè)置為1 次, 內(nèi)循環(huán)設(shè)置為30 次, 對南海中尺度渦旋進(jìn)行了系統(tǒng)的數(shù)據(jù)同化實(shí)驗(yàn)。同化的時間范圍為2012 年7 月1 日—2013年1 月30 日, 并主要對包括南海在內(nèi)的105°—123°E,9°—26.24°N 的中尺度渦旋進(jìn)行了統(tǒng)計分析。

2.1 同化結(jié)果的驗(yàn)證

將相同時間范圍與空間范圍內(nèi)數(shù)據(jù)同化(I4D-Var)結(jié)果、AVISO 實(shí)際觀測和控制實(shí)驗(yàn)(control run, CR)中的中尺度渦旋個數(shù)進(jìn)行統(tǒng)計, 其結(jié)果如表1 所示。由于AVISO 數(shù)據(jù)的分辨率為0.25°, 因此表1 中只統(tǒng)計渦旋半徑大于60 km 的渦旋個數(shù)。

表1 數(shù)據(jù)同化(I4D-Var)結(jié)果、AVISO 實(shí)際觀測和控制實(shí)驗(yàn)(CR)中渦旋個數(shù)統(tǒng)計Tab.1 Statistics of the number of vortices in I4D-Var, AVISO, and CR

由表1 的對比可以看到, 在AVISO 觀測的渦旋中冷暖渦的個數(shù)相當(dāng), 約1︰1, 該結(jié)果與全球范圍內(nèi)冷暖渦的比例接近(Cheltonet al, 2011), 但是冷渦個數(shù)稍大于暖渦, 該現(xiàn)象在I4D-Var 和CR 中也都得到了較好地模擬。I4D-Var 統(tǒng)計到的總渦旋個數(shù)與實(shí)際觀測中的個數(shù)更為接近, 相差僅8%, 冷渦和暖渦分別比觀測多了13.6%和12%。I4D-Var 結(jié)果相比CR, 其渦旋個數(shù)更接近為真值。而未同化的結(jié)果則差別很大,CR 中的渦旋個數(shù)偏多, 相差27%。找出整個同化時間段內(nèi)每天AVISO 和I4D-Var 一一對應(yīng)的渦旋, 對于半徑在60 km 以上的渦旋, I4D-Var 的成功識別率為90%。然后分別計算了兩者的渦心距離差距、半徑相對誤差和振幅相對誤差, 如圖 5 所示。可以看到I4D-Var 與AVISO 一一對應(yīng)的渦旋中, 其平均渦心距離有4.3 海里的差距(圖5a), 平均振幅相對誤差為0.16 (圖5c), 而平均半徑相對誤差很小, 僅有0.07(圖5b), 上述結(jié)果充分說明經(jīng)過數(shù)據(jù)同化的模擬場可以再現(xiàn)真實(shí)海洋中的中尺度渦旋, 而且與觀測渦旋的位置、振幅、半徑大小也都十分吻合。

圖5 數(shù)據(jù)同化方法(I4D-Var)和AVISO 對應(yīng)渦旋的渦心差距(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 三者在同化時間段內(nèi)渦旋出現(xiàn)的位置分布進(jìn)行對比。由圖6b 的衛(wèi)星觀測數(shù)據(jù)可以看到, 本次同化時間范圍內(nèi)中尺度渦旋主要出現(xiàn)在南海南部, 尤其是在越南沿岸東側(cè)渦旋的出現(xiàn)頻率最高, 其原因主要與季風(fēng)驅(qū)動和地形作用有關(guān)(Chenet al, 2009)。其余沿著118°E 也有兩個出現(xiàn)頻率相對集中的位置分別位于13°N 和15°N附近。在該時間段內(nèi), 雖然南海北部的渦旋出現(xiàn)頻率很低, 但是在呂宋海峽出入口的渦旋出現(xiàn)頻率普遍較高, 其中臺灣西南部是南海北部渦旋出現(xiàn)頻率的高值區(qū)域。有關(guān)臺灣島西南部的暖渦生成機(jī)制仍然是許多學(xué)者關(guān)注的熱點(diǎn)問題。而I4D-Var 的渦旋出現(xiàn)頻率無論是在強(qiáng)度還是空間分布與AVISO 都有很好的對應(yīng)關(guān)系, 臺灣島西南部渦旋的出現(xiàn)頻率要略高于觀測結(jié)果。而CR 中渦旋的出現(xiàn)位置與觀測相比相差比較大。由此可見, 數(shù)值模式經(jīng)過數(shù)據(jù)同化后能夠很好地產(chǎn)生與AVISO 對應(yīng)的中尺度渦, 而為加入數(shù)據(jù)同化的模擬中渦旋出現(xiàn)位置的分布則具有非常強(qiáng)的隨機(jī)性。

圖6 I4D-Var (a)、AVISO(b)和控制實(shí)驗(yàn)(CR)(c)渦旋出現(xiàn)頻率Fig.6 Frequency of eddies appearance in I4D-Var (a), AVISO (b), and CR (c)

接著將同化時間段內(nèi)所有的渦旋半徑進(jìn)行了比較分析, 結(jié)果如圖7 所示。結(jié)果表明, I4D-Var 和CR中半徑在60 km 以下的渦旋分別占了總渦旋數(shù)量的34%和28%, 而AVISO 數(shù)據(jù)識別的渦旋大部分都在60 km 以上(92%), 造成該差異最主要因素可能是分辨率的不同, 如果只統(tǒng)計半徑大于60 km 的渦旋比例分布, 則I4D-Var 與AVISO 觀測結(jié)果非常吻合, CR 結(jié)果則相差很多(表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)

為進(jìn)一步證明同化結(jié)果的有效性, 還將I4D-Var、AVISO 及CR 中渦旋的半徑空間分布進(jìn)行對比。從圖8 中可以看到I4D-Var 的渦旋半徑分布(圖8a)與實(shí)際觀測(圖 8b)非常一致, 大半徑的渦旋主要分布于南海南部 115°E 以西的地方, 且半徑普遍要大于100 km。南海北部大半徑渦旋較少, 出現(xiàn)頻率較多的臺灣島西南部渦旋半徑主要在80 km 左右。而CR(圖8c)中無論是渦旋半徑空間位置分布還是大小分布都與觀測相差甚遠(yuǎn)。上述的統(tǒng)計結(jié)果間接的表明了I4D-Var 相對于CR 具有明顯的優(yōu)越性, 證明了在模式中加入四維變分同化對中尺度渦的模擬是非常有效的。

圖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 日時間段的個例分析來直觀的看一下加入同化的數(shù)值模擬效果。圖9 是AVISO 的海表高度異常(sea surface height anomalies, SSHA)和地轉(zhuǎn)流分布圖, 在11 月13—28 日期間南海海盆西部有一個越南冷渦(Vietnam cyclonic eddy, VCE)和一個越南暖渦(Vietnam anticyclonic eddy, VAE)組成的渦旋對, 該渦旋對通常在越南東岸周期性出現(xiàn), 所以也將其統(tǒng)稱為越南渦旋對(Vietnam eddy pair, VEP) (Chenet al,2010)。Chen 等(2010)研究表明越南渦旋對通常發(fā)生在夏秋季節(jié), 局地風(fēng)應(yīng)力是該渦旋對年際性出現(xiàn)的主要原因。另外一個典型的渦旋就是南海東北部發(fā)展起來一個比較強(qiáng)的暖渦, 該暖渦位于臺灣西南部, 所以也常被成為臺灣暖渦(Taiwan anticyclonic eddy,TAE)。TAE 經(jīng)常發(fā)生在冬季, 有關(guān)TAE 的生成機(jī)制有許多的研究, 其中Nan 等(2011)提出了渦旋與黑潮相互作用模型, 認(rèn)為黑潮路徑的變化是該地區(qū)渦旋形成的主要誘因。冬季由于受東北季風(fēng)影響, 黑潮環(huán)流較為頻繁的出現(xiàn)在臺灣東南部, 極易發(fā)生流套脫離現(xiàn)象, 從而發(fā)展為暖渦。

圖9 AVISO 海面高度異常SSHA(單位: cm)以及地轉(zhuǎn)流(單位: 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 和地轉(zhuǎn)流分布圖, 從圖中可以看出, I4D-Var的SSHA 和流場的空間分布非常接近AVISO 的觀測結(jié)果。圖中各個渦旋的位置可以互相對應(yīng), 除了上述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 結(jié)果中均可以看到一個由LCE、CCE 和VCE 貫穿南海中部的冷渦序列結(jié)構(gòu),Chu 等(2020)基于衛(wèi)星觀測和現(xiàn)場觀測數(shù)據(jù)的統(tǒng)計結(jié)果指出南海冬季沿17°N 和西邊界形成的反“L”型渦列呈周期性出現(xiàn)。

圖10 I4D-Var 海海面高度異常(單位: cm)以及地轉(zhuǎn)流(單位: cm/s)演變Fig.10 Evolution in SSHA (unit: cm) and geostrophic current (unit: cm/s) of I4D-Var

為了進(jìn)一步檢驗(yàn)?zāi)Ｊ郊尤胪髮?、鹽場的調(diào)整, 選取了某時刻沿115°E (10°—18°N, 水深<300 m)的溫鹽剖面圖, 該經(jīng)度恰好穿過一個冷渦和一個暖渦。圖11a、11b 上方是該剖面上AVISO(黑線)、I4D-Var(紅線)和CR(藍(lán)線)的SSHA 比較結(jié)果, 通過AVISO 結(jié)果可以看到, 該剖面上SSHA 分別在13°N 和14.7°N,出現(xiàn)極小值(0.02 m)和極大值(0.27 m), 對應(yīng)著一個冷渦和一個暖渦。其中, I4D-Var 的SSHA 極小值點(diǎn)與AVISO 基本一致, 極大值點(diǎn)位于15.2°N 附近, 相比AVISO 稍微向北偏移0.5°。而CR 的SSHA 普遍小于AVISO 和I4D-Var, 在16.6°N 有一個極大值點(diǎn)(0.12 m), 但是距離AVISO 和I4D-Var 的極大值距離比較遠(yuǎn),其結(jié)果不確切。經(jīng)計算, CR 的SSHA 均方根誤差(RMSE)為14%, 而I4D-Var 的RMSE 為5%, 相比CR得到了8%的優(yōu)化。比較I4D-Var 和CR 的溫鹽剖面可以看到, 圖11a 冷渦附近(藍(lán)色陰影)的等溫線有明顯的向上彎曲, 渦旋內(nèi)的水溫比周圍水體要小, 且該冷渦的溫度中軸線向南偏移; 暖渦(紅色陰影)內(nèi)的等溫線則向下彎曲, 表層的海水幅聚下沉, 其18°C 等溫線最深可達(dá)90 m 的位置; 而在未同化的CR 實(shí)驗(yàn)中(圖11b)兩個陰影位置的等溫線都比較平緩, 沒有出現(xiàn)明顯的彎曲情況。同樣地, 在鹽度剖面圖中,I4D-Var 中冷渦內(nèi)(圖11c 藍(lán)色陰影)的鹽度變大, 下層的高鹽水有強(qiáng)烈的上升趨勢; 暖渦內(nèi)(圖 11c 紅色陰影)的鹽度等值線明顯向下彎曲, 但是沒有冷渦彎曲的強(qiáng)烈; 反觀CR 的鹽度剖面圖(圖11d)兩個陰影內(nèi)的等值線都較為平緩, 說明這里沒有渦旋存在。經(jīng)過上述同化結(jié)果與AVISO 觀測的一系列比較表明, 相比單純模式模擬的結(jié)果, 加入四維變分同化后的模式可以極大地提高對南海中尺度渦旋的模擬能力。在一個渦旋剖面?zhèn)€例的驗(yàn)證中發(fā)現(xiàn)同化后的模式相比不同化對SSHA 模擬的優(yōu)化可達(dá)8%, 而且模式加入同化后不僅可以改變海面高度, 其最終還會導(dǎo)致整個溫度、鹽度和密度場的改變。由此可以證明模式加入數(shù)據(jù)同化之后可以得到一個接近真實(shí)的海洋動力學(xué)狀態(tài), 在該基礎(chǔ)上可以展開中尺度渦的后報模擬及時效性檢驗(yàn)工作。

圖11 沿115°E 的溫鹽剖面圖Fig.11 Profiles of temperature and salinity along 115°E

2.2 南海中尺度渦旋的后報模擬實(shí)驗(yàn)及其時效性檢驗(yàn)

通過上面的分析可以看出, I4D-Var 能夠很好地模擬出與觀測結(jié)果相符的中尺度渦, 但是該結(jié)果是在有觀測資料同化的前提下。那么, 如果只有較為準(zhǔn)確的初始場, 而沒有同化的幫助, 僅僅靠海洋模式是否能夠?qū)χ谐叨葴u旋的進(jìn)一步運(yùn)動進(jìn)行準(zhǔn)確地模擬呢？眾所周知, 中尺度渦雖然稱之為中尺度, 但實(shí)際上已屬于大尺度海洋現(xiàn)象的范疇, 它的運(yùn)動受到準(zhǔn)地轉(zhuǎn)平衡關(guān)系的約束。因此, 從理論上講, 像風(fēng)場、水熱通量等外在的強(qiáng)迫場雖然對中尺度渦有長期的影響, 但是短期內(nèi)對中尺度渦的運(yùn)動和發(fā)展影響不大(高山等, 2007)。因此, 可以想象, 借助中尺度渦的這種地轉(zhuǎn)平衡約束的慣性, 使用模式模擬在短期內(nèi)是應(yīng)當(dāng)可以保證渦旋模擬的大致準(zhǔn)確性。為了驗(yàn)證該假設(shè)的可行性, 我們決定使用同化模擬得到的結(jié)果作為初始場, 然后僅依靠ROMS 海洋模式, 對中尺度渦進(jìn)行后報實(shí)驗(yàn), 以初步檢驗(yàn)上述猜測的合理性, 進(jìn)而確定中尺度渦的預(yù)報時效。具體方案如下:

(1) 同化進(jìn)行到2012 年12 月8 日時停止, 得到該時間段的同化場結(jié)果;

(2) 將(1) 中的結(jié)果作為初始場, 保持海表強(qiáng)迫場和側(cè)邊界場不變, 僅靠模式本身運(yùn)行, 對中尺度渦進(jìn)行后報模擬, 運(yùn)行時間為3 個月(2012 年12 月8 日—2013 年3 月8 日);

(3) 將后報模擬(hindcast control, HC)結(jié)果與相同時間段的AVISO 高度計資料進(jìn)行對比。

需要指出的是本文中所說的預(yù)報, 實(shí)際上是在同化結(jié)果提供初始場的基礎(chǔ)上進(jìn)行的后報模擬, 主要用于檢測同化模擬的有效性。

分別選取HC 在第5、10、15 和20 天的SSHA結(jié)果與觀測結(jié)果進(jìn)行對比, 結(jié)果如圖12 所示。首先HC 模擬的SSHA 分布與AVISO 十分吻合, 該時間段內(nèi)HC 模擬的南海中尺度渦與AVISO 也有很好的對應(yīng)關(guān)系, 其中TAE、HCE、VCE 和CAE 在后預(yù)報的20 d 內(nèi)均持續(xù)存在, VAE 也明顯地存在于前10 d, 之后逐漸消亡。分別計算后預(yù)報20 d 內(nèi)HC 和AVISO之間對應(yīng)渦旋的SSHA 相關(guān)系數(shù)、渦心差距和半徑絕對誤差, 如圖13a 所示, HCE、VCE 和CAE 在前20 d內(nèi)的SSHA 相關(guān)系數(shù)均在0.7 以上。TAE 在前13 dSSHA 相關(guān)系數(shù)均在0.85 以上, 但之后有減小趨勢,其可能原因是AVISO 中TAE 逐漸向西南方向移動(圖12e 和12g), 而HC 中雖有移動, 但很緩慢。VAE 的SSHA 相關(guān)系數(shù)在前10 d 皆在0.6 以上, 10 d 以后由于VAE 開始消亡, 相關(guān)系數(shù)迅速減小。另外, 圖13b可以看出, 除了VCE 的渦心差距在大部分時間內(nèi)保持一致外, 其余渦旋的渦心差距雖然均有增大的趨勢, 但是皆在30 海里范圍內(nèi), 平均渦心差距僅為15海里。VCE、CAE 和VAE 由于渦旋形態(tài)的高時變性,導(dǎo)致三者的半徑絕對誤差隨時間略有增加, 但是基本維持在30 海里以內(nèi), 而TAE 和HCE 較為穩(wěn)定, 其平均半徑絕對誤差分別為3.5 和6.6 海里。由此可見,上述后報結(jié)果驗(yàn)證了我們的猜測, 即由同化結(jié)果而來的初始場具備了中尺度渦巨大的準(zhǔn)地轉(zhuǎn)運(yùn)動慣性,僅使用ROMS 海洋模式對南海中尺度渦旋進(jìn)行預(yù)報模擬是可行的, 對于本次南海中尺度渦旋的后報模擬, 預(yù)報時效可達(dá)10 d 以上。

為了更清晰地量化中尺度渦的預(yù)報模擬效果,本研究進(jìn)一步計算了所有中尺度渦(圖12 中標(biāo)注渦旋)運(yùn)動學(xué)參數(shù)的平均相關(guān)系數(shù), 并以此作為預(yù)報的評估指標(biāo)。所選取的運(yùn)動學(xué)參數(shù)主要包括: SSHA、相對渦度(relative vorticity, RV)、渦動能(eddy kinetic energy, EKE)和水平散度(divergence, Div), 相應(yīng)的計算公式如下:

式中,VR表示相對渦度;EEK表示渦動能;D表示水平散度; ～u和v～分別代表地轉(zhuǎn)速度異常的緯向分量和經(jīng)向分量:

4 個運(yùn)動學(xué)參數(shù)的平均相關(guān)系數(shù)隨時間變化如圖14 所示, 可以看到在后預(yù)報的20 d 中, 上述4 個參數(shù)的平均相關(guān)系數(shù)均在0.67 以上, 其中EKE 和RV 的相關(guān)性甚至皆在0.78 以上, 該結(jié)果進(jìn)一步證明四維變分同化技術(shù)可以提供合理的預(yù)報初始場, 能極大地提高中尺度渦的短期預(yù)報準(zhǔn)確率。

圖14 渦旋動力學(xué)參數(shù)平均相關(guān)系數(shù)的時間變化Fig.14 Time series of mean correlation coefficient of eddy kinematic parameters

通過圖12 和圖13 可以看到TAE 相比模擬區(qū)域中的其他渦旋都要強(qiáng)大, 半徑大小與觀測高度吻合,渦旋形態(tài)的時變性小, 并且持續(xù)地向西南方向移動,下面以TAE 作為研究個例, 繼續(xù)檢驗(yàn)經(jīng)四維變分同化優(yōu)化過的初始場對于中尺度渦旋的預(yù)報效果。繼續(xù)探測并追蹤后報模擬中渦旋的變化, 發(fā)現(xiàn) TAE 自2012 年12 月8 日開始繼續(xù)運(yùn)行了86 d, 一直到2013年3 月2 日終止。而在AVISO 觀測中該渦旋生命周期有89 d (2012 年11 月28 日—2013 年2 月24 日)與預(yù)報結(jié)果比較一致。兩者運(yùn)行路徑的比較如圖15所示, 兩者皆產(chǎn)生于南海東北部約118.5°E, 21°N 附近, 然后向西南方向運(yùn)動, 并于113.5°E, 17.5°N 附近消亡。圖15 可以看出HC 模擬的移動路徑和AVISO大致吻合, 說明同化提供的初始場對于預(yù)報中尺度渦位置的演變發(fā)展有很好的幫助。此外還比較了渦旋在HC 模擬和AVISO 觀測兩種情況下半徑和振幅相對誤差的時間變化, 如圖16 所示, TAE 振幅和半徑的變化趨勢幾乎是同步的, 其中, 前30 d 內(nèi)TAE 平均半徑相對誤差和平均振幅相對誤差分別為-0.01 和-0.13, 證明前30 d, TAE 的后報模擬半徑和振幅與觀測高度吻合。30—65 d, HC 模擬的渦旋經(jīng)歷了一個半徑和振幅減小又增大的過程, 直至最后消亡的階段,HC 模擬的渦旋半徑和振幅又恢復(fù)到和觀測接近的狀態(tài)。上述結(jié)果表明后預(yù)報場可以預(yù)報模擬個別強(qiáng)壯渦旋整個生命周期的移動路徑, 并且在短時間內(nèi)高度還原真實(shí)渦旋的大小和振幅變化, 其他運(yùn)動學(xué)參數(shù)的預(yù)報模擬結(jié)果還需要進(jìn)一步的驗(yàn)證。

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

圖13 HC 和AVISO 對應(yīng)渦旋的SSHA 相關(guān)系數(shù)(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 臺灣暖渦運(yùn)動軌跡圖比較Fig.15 Comparison of TAE trajectory

圖16 渦旋半徑和振幅相對誤差的時間變化Fig.16 Time series of relative error of eddy radius and amplitude

綜合上述對比結(jié)果進(jìn)一步表明, 模式的后報結(jié)果的確在至少10 d 內(nèi)依然與AVISO 觀測分析結(jié)果非常吻合。個別成熟且穩(wěn)定發(fā)展的渦旋, 比如TAE 甚至能夠預(yù)報整個生命周期的移動路徑。后報模擬因?yàn)橛辛怂木S變分同化給出的合理初始場, 可以極大地提高中尺度渦旋模擬的準(zhǔn)確性, 究其原因有如下兩點(diǎn):

(1) 得益于ROMS 成熟的四維變分同化模塊給后報模擬提供了一個接近真實(shí)海洋狀態(tài)的初始場,其海面高度和流場與實(shí)際觀測結(jié)果非常接近, 而且其對應(yīng)的溫度、鹽度和密度場均得到調(diào)整。

(2) 海洋中地轉(zhuǎn)平衡的調(diào)整是一個緩慢的過程,對于已經(jīng)形成的中尺度渦來說, 其運(yùn)動和發(fā)展主要受到地轉(zhuǎn)平衡的約束, 而快速變化的強(qiáng)迫場在短時間內(nèi)對其影響不大。

因此, 中尺度渦一旦形成, 則認(rèn)為模式達(dá)到一個接近真實(shí)海洋地轉(zhuǎn)平衡的狀態(tài), 在這種地轉(zhuǎn)平衡的“慣性”驅(qū)動下, 模式可以很好地預(yù)測中尺度渦在短期內(nèi)的運(yùn)動狀態(tài)。

3 結(jié)論

本文首先利用四維變分?jǐn)?shù)據(jù)同化方法將AVISO衛(wèi)星高度計資料同化到ROMS 模式中, 并對同化結(jié)果進(jìn)行了統(tǒng)計分析; 然后使用同化模擬得到的結(jié)果作為初始場, 僅依靠ROMS 海洋模式, 對中尺度渦進(jìn)行后預(yù)報實(shí)驗(yàn), 探討了經(jīng)數(shù)據(jù)同化優(yōu)化的初始場對于中尺度渦旋后預(yù)報模擬的可行性和有效性。結(jié)論如下:

加入四維變分?jǐn)?shù)據(jù)同化之后的ROMS 模式能夠極大的提高南海中尺度渦的模擬精度。對比同化期間內(nèi)衛(wèi)星觀測和同化結(jié)果的渦旋位置, 發(fā)現(xiàn)同化可以一一再現(xiàn)90%的南海中尺度渦旋, 其中平均渦心距離為4.3 海里, 平均振幅相對誤差為0.16, 平均半徑相對誤差僅為0.07。此外, 經(jīng)同化再現(xiàn)的中尺度渦旋的數(shù)量、出現(xiàn)頻率和半徑的空間分布與觀測都非常吻合, 而單獨(dú)模式模擬的渦旋產(chǎn)生位置和大小則具有很大的隨機(jī)性。

同化方法是通過對整個模擬場的海水狀態(tài)進(jìn)行調(diào)整產(chǎn)生的中尺度渦旋, 并使得渦旋的運(yùn)動得以維持。通過對一個同時穿越冷渦和暖渦的剖面對比發(fā)現(xiàn),相比未加入同化的結(jié)果, 加入同化后SSHA 的均方根誤差得到了8%的優(yōu)化, 而且同化模擬中冷、暖渦對應(yīng)的溫度和鹽度場均產(chǎn)生相應(yīng)的變化。

同化模擬可以幫助模式得到一個接近真實(shí)的中尺度渦模擬場, 以這個穩(wěn)定的平衡場作為初始場, 進(jìn)行正常的南海中尺度渦旋的后報模擬發(fā)現(xiàn), 模式的后預(yù)報能力可達(dá)10 d 以上。將個別穩(wěn)定并強(qiáng)壯發(fā)展的渦旋的運(yùn)動軌跡、半徑及振幅信息與實(shí)際觀測進(jìn)行了比較, 結(jié)果表明其生命周期和運(yùn)動軌跡幾乎與觀測吻合, 而且短時間內(nèi)高度還原真實(shí)渦旋的大小和振幅變化。

值得注意的是同化得到的預(yù)報初始場實(shí)際上是一個中尺度渦旋的地轉(zhuǎn)平衡狀態(tài), 所謂的預(yù)報, 則是中尺度渦在地轉(zhuǎn)平衡“慣性”驅(qū)動下的結(jié)果。而諸如風(fēng)暴潮、臺風(fēng)等極端天氣都有可能打破上述地轉(zhuǎn)平衡,下一步可以對這些極端天氣進(jìn)行單獨(dú)的討論研究。本文中一個隱含的假設(shè)是模式、表面邊界條件和側(cè)向開邊界條件是準(zhǔn)確的, 但實(shí)際上, 四維變分同化中初始場、表面強(qiáng)迫場、背景場、模式和觀測的誤差協(xié)方差矩陣都會影響到代價函數(shù)的計算。因此, 下一步將采用國內(nèi)外認(rèn)可的實(shí)時強(qiáng)迫場數(shù)據(jù)繼續(xù)調(diào)試優(yōu)化模式本身的模擬能力, 并在同化中加入Argo、水下滑翔機(jī)、潛標(biāo)等現(xiàn)場觀測的溫鹽剖面數(shù)據(jù), 進(jìn)行南海中尺度渦的實(shí)際預(yù)報實(shí)驗(yàn)。