劉小璐, 聶 銘*, 羅嘯宇, 楊 劍, 段忠東
(1.廣東電網(wǎng)有限責任公司電力科學研究院, 廣州 510080; 2.哈爾濱工業(yè)大學(深圳)土木與環(huán)境工程學院, 深圳 518055)
在臺風預報及其相關災害的應急、預防等方面,數(shù)值臺風模式已成為不可或缺的一部分[1-5]。得益于現(xiàn)代計算機技術的發(fā)展,多種中尺度模式[如MM5(the fifth-generation Penn State/NCAR mesoscale model)及WRF(weather research and forecasting model)]被應用在臺風相關的研究中。目前,在臺風研究中,使用最廣泛的中尺度模式為WRF[6],該模式由美國國家大氣研究中心 (National Center For Atmospheric Research, NCAR)、美國國家海洋和大氣管理局 (National Oceanic And Atmospheric Administration, NOAA) 等機構聯(lián)合多所大學共同完成。WRF模式可以使用歷史再分析資料對發(fā)生過的天氣現(xiàn)象進行重演,且其繼承和發(fā)展了MM5中的多種微物理參數(shù)化方案和多重嵌套網(wǎng)格技術等優(yōu)勢,成為眾多學者研究臺風熱動力結構的重要工具[7-9]。WRF模式采用了地形跟隨質量坐標,可以模擬復雜地形條件下的氣流運動。Xue等[10]采用WRF模式模擬1213號臺風“啟德”和1409號臺風“威馬遜”,并通過修改物理參數(shù) (如:臺風中心氣壓、地形下墊面等) 模擬了瓊州海峽百年一遇的臺風風速,結果表明WRF模式模擬的風攻角、最大風速以及風速風向時程等,可應用于工程抗風設計。目前,WRF模式對于單個臺風案例可以有較好的模擬結果;但受限于歷史再分析資料較短的觀測年限,很難實現(xiàn)WRF模式對較長時間年限 (如100年) 的眾多臺風案例模擬,因此,很難滿足臺風危險性分析的要求。
從20世紀開始,隨著人們對臺風觀測手段的不斷提升,逐步發(fā)展出了一套簡化的、可適用工程大量模擬的參數(shù)化臺風風場模型。這種參數(shù)化臺風風場模型僅需輸入臺風中心經(jīng)緯度、臺風中心氣壓等數(shù)個信息就可以大致模擬出臺風風場結構。對于重點關注極值風速、需要大量模擬臺風案例的工程抗風設計來說,這種簡化的參數(shù)化模型表現(xiàn)出了極大的優(yōu)勢。參數(shù)化臺風風場模型模擬近地面風速通常分為以下兩步:首先,通過臺風中心經(jīng)緯度、臺風中心氣壓 (臺風中心附近最大風速) 等,根據(jù)梯度風平衡條件,反演出梯度風速;然后,通過風速折減系數(shù)[10]或邊界層模式[11-12]把梯度風轉換到近地面。這種方法重點考慮了臺風的動力學特性,而忽略了不關心的熱力學效應和水汽循環(huán)等過程。近些年了,學者們重點研究了臺風風場的邊界層模式。Wu等[13]采用Arya邊界層方案建立了梯度風速和近地表風速的關系,并開發(fā)了一套簡化公式,以提高模型計算效率。Fang等[14]在Meng風場[11]的基礎之上考慮了氣壓隨高度變化,并利用觀測資料擬合出了渦旋黏滯系數(shù)K,驗證結構表明該模型能用較好的重現(xiàn)出風速分析時程及邊界層臺風風場結構。Hong等[15]通過邊界層動力學平衡方程,在柱坐標下采用隱式差分方法,建立了數(shù)值邊界層風場模型。相比于其他模型,該方法直接從控制方程出發(fā),減少了經(jīng)驗性公式的應用,從某種程度上來說更能客觀反映臺風的動力學本質。
然而,這些參數(shù)的臺風風場模型,通常都是建立在平坦下墊面的,如海洋表面。這使得當臺風登陸以后,依然進行平坦下墊面假設會變得于實際不符,特別是當臺風經(jīng)過高山時必然會對邊界層氣流產(chǎn)生阻礙效應。鄒振操等[7]使用WRF模式對比了不同精度的地形資料,表明高精度的地形數(shù)據(jù)能在一定程度上改善近地面風場的模擬效果。另外,陸地上變化著的地表覆蓋也會對臺風風場產(chǎn)生顯著影響。馮簫等[16]采用WRF模式研究了海南島三種土地覆蓋資料對模擬結果的影響。對比地面10 m 高度風速場的模擬結果,他們發(fā)現(xiàn)采用具有更精細的土地覆蓋分類數(shù)據(jù),可以增大下墊面的不均勻性,提高了地表的粗糙長度,使得模擬的地面風與實際風更貼近。
因此,在參數(shù)化臺風模型中,引入地形和土地覆蓋,將更有利于刻畫出與實際情況相符的臺風風場。參考WRF模式的垂直坐標系,將把參數(shù)化臺風模型建立在隨地形變化的坐標系下來考慮地形起伏對風場的影響。這種坐標系可將模型控制方程投影到一個隨地形起伏變化而變化的計算網(wǎng)格上,避免了氣壓或者是等溫線與地面相交而帶來的復雜邊界,同時也可以較好地描述連續(xù)場 (如空氣的流動和氣溫的改變) 的變化[17]。對于土地覆蓋的影響,較為成熟的做法是把不同土地覆蓋類型等效為地表粗糙長度。李軍等[18]將21種土地覆蓋類型對應到8種地表粗糙長度上,研究中國地區(qū)的地表粗糙長度變化。李沁怡等[19]利用歐洲空間局 (European Space Agency, ESA) 發(fā)布的高分辨率全球土地覆蓋數(shù)據(jù),編制了水平分辨率約為300 m的中國地表粗糙長度分布圖,并與實測數(shù)據(jù)進行了對比驗證。因此,將土地覆蓋變化的信息隱含到地表粗糙長度中,并在隨地形變化的坐標系下進行參數(shù)化臺風風場運算,可得到考慮地形起伏變化影響的臺風風場。
在建模過程,首先需要定義計算域。水平方向有:xmin≤x≤xmax,ymin≤y≤ymax;垂直方向有hs(x,y)≤z≤H,其中,hs(x,y)為地形海拔高度 (圖1),H為上邊界,設置為10 000 m。直接在hs(x,y)上建模會導致非正交網(wǎng)格,不方便模型運算。為此,Gal-Chen等[20]提出了一種坐標轉換方法,將hs(x,y)上的非平坦計算域轉換到隨地形起伏的坐標系下 (σ坐標系),這樣就保證了整個計算域可建立在正交網(wǎng)格上。保持水平方向不變,即xσ=x,yσ=y(xσ、yσ為σ坐標系下的兩個水平坐標);垂直方向變?yōu)?/p>
圖1 地形剖面及垂直分層示意圖
(1)
式(1)中:z為原始坐標系下的垂直坐標,zσ為σ坐標系下的垂直坐標。
由國際水道測量組織(International Hydrographic Organization, IHO)和聯(lián)合國教科文組織下屬的政府間海洋委員會(Intergovernmental Oceanographic Commission, IOC)聯(lián)合指導的大洋水深制圖項目 (General Bathymetric Chart of the Oceans, GEBCO)可以提供分辨率為15″的全球地形數(shù)據(jù)。該套數(shù)據(jù)是目前公開發(fā)布的,地形分辨率最高的一類,滿足中尺度模型的對地形分辨率的需要。使用GEBCO 2020柵格數(shù)據(jù),其中雷州半島及海南島地區(qū)的地形高程示意如圖2所示。
該圖基于國家測繪地理信息局標準地圖服務網(wǎng)站下載的審圖號為GS(2020)4634號的標準地圖制作,底圖無修改
歐洲空間局 (ESA) 氣候變化協(xié)議署 (Climate Change Institute, CCI) 最新公布了2018年全球土地覆蓋圖,其中共包含了22種土地覆蓋分類,水平分辨率最高可達到300 m,是目前公開發(fā)布的,地表覆蓋分辨率最高的一類,滿足中尺度模型對不同地表覆蓋的識別需要。由于李沁怡等[19]已經(jīng)參考WRF模式,建立并驗證了ESA地表覆蓋數(shù)據(jù)與地表粗糙長度的映射關系。因此,直接使用這種映射關系[19]。考慮到登陸中國的臺風多發(fā)生在夏季,選取EW中的較大值是合適的,最后可得到雷州半島及海南島地區(qū)的地表粗糙長度示意,如圖3所示。
該圖基于國家測繪地理信息局標準地圖服務網(wǎng)站下載的審圖號為GS(2020)4634號的標準地圖制作,底圖無修改
Meng等[11]提出了參數(shù)化臺風風場模型,由于其控制方程簡單、計算效率高、沿高度方向解析等優(yōu)點,在工程抗風中得到了普遍應用。在Meng模型的基礎之上,Kepert[21]模型包含了更多的物理過程,即考慮了邊界層內的水平對流。因此,基于Kepert模型對參數(shù)化臺風風場模型進行改進使之能考慮地形起伏的影響。
Meng模型[11]和Kepert模型[21]都將邊界層內的臺風風速分解為梯度風速vg和摩擦引起的風速(u′,v′)。梯度風速可根據(jù)梯度高度處氣壓梯度力、科氏力和離心力等的平衡求得[11]
(2)
式(2)中:f為科氏參數(shù);ρ為空氣密度;r為梯度風速離臺風中心的距離;p為臺風氣壓場;cλ為臺風移動速度沿λ的分量。摩擦引起的風速(u′,v′)經(jīng)簡化后,在σ坐標系(r,λ,zσ)下可表示為[21]
(3)
(4)
式中:K為渦旋黏滯系數(shù),可近似取K=50 m2/s。求解式(3)、式(4)需要先設定下邊界條件,在近地面(z→0)時[21],得
(5)
(6)
式中:κ=0.4為von K?rm?n常數(shù);z0為地表粗糙長度;cr為臺風移動速度沿r的分量。這里需要說明的是,參數(shù)化臺風風場中的地表粗糙長度z0可用圖3中的空間地表粗糙長度代替。參考Kepert[21]的求解方法,把控制式(3)、式(4)代入邊界條件式(5)、式(6)中,通過傅里葉變換和代數(shù)運算,可得臺風邊界層風速(u,v)為
u=u′+cr
(7)
v=v′+vg+cλ
(8)
為了驗證改進臺風風場模型的可靠性,將模擬對比歷史上登陸廣東省的3場臺風,即2018年的臺風山竹、2015年的臺風彩虹和2014年的臺風威馬遜。
2018年9月,臺風山竹 (Mangkhut) 在西北太平洋洋面生成;于2018年9月16日17:00在廣東臺山海宴鎮(zhèn)登陸,登陸時臺風中心附近最大風速達到45 m/s,臺風中心最低氣壓為955 hPa,造成廣東、廣西、海南、湖南、貴州等地區(qū)近300萬人受災。2015年10月1日,超強臺風彩虹 (Mujigae) 在西北太平洋菲律賓群島附近形成;于2015年10月4日14:10在廣東省湛江登陸,是1949年以來10月份登陸中國的最強臺風。超強臺風威馬遜 (Rammasun) 在2014年7月9日于西北太平洋海面生成,隨后其一路向西北偏西移動,穿越菲律賓中部,進入南海后,臺風強度迅速加強,并先后于海南、廣東、廣西登陸。臺風威馬遜是1973年來登陸華南地區(qū)的最強臺風,登陸時臺風中心最低氣壓將近 900 hPa,給海南、廣東、廣西等地帶來嚴重災害和巨大經(jīng)濟損失[22]。
北京時間2018年9月16日19時,臺風山竹的中心位于陽江市,其中心位置為(21.9°N,112.0°E),中心最低氣壓為960 hPa,臺風移動速度為 9.14 m/s,移動方向為北偏西71.55°,地形和地表粗糙長度同圖2、圖3一致。由改進的臺風風場模型計算得的不同高度處的臺風風場如圖4所示?!昂Q笙聣|面”是指計算域下墊面設置與海洋情況相同,即hs(x,y)=0,z0=0.000 2;“僅考慮土地覆蓋”是指計算域下墊面考慮了地表覆蓋的變化但忽略了地形的變化,即hs(x,y)=0;“僅考慮地形變化”是指計算域下墊面考慮了地形的起伏變化,而地表覆蓋類型假設為海洋,即z0=0.000 2;“考慮土地覆蓋和地形變化”是指計算域中考慮了圖2、圖3中的實際地形和實際土地覆蓋類別。
在圖4中,對比了不同下墊面條件下,水平臺風風場沿高度變化的情況。從圖4中可以明顯地看出,在10 m高度的風場中,僅考慮土地覆蓋的變化時,會使得模擬的洋面風速大于陸地風速;但隨著離地高度的增加,土地覆蓋對風場的影響逐漸減弱,直至消失 (500 m高度處“僅考慮土地覆蓋”情況和“海洋下墊面”情況的風場基本相同)。在本文模型中,當僅有地形起伏作用時,風場的最大風速位置與“海洋下墊面”模擬情況類似;但僅考慮地形時,模擬的風場風速要比“海洋下墊面”情況的風速偏小,特別是在臺風中心區(qū)域。這表明在本模型中,地形對臺風風場主要起阻擋作用。地形起伏變化,會導致上升/下降氣流的加強/減弱,進而影響水平氣流的運動,但在本文模型中,忽略了空氣垂直對流作用,這一簡化可能是導致地形對風的加速作用沒有體現(xiàn)出來的原因之一。當考慮實際地形和土地覆蓋時,模擬得的10 m高度處風場比“海洋下墊面”情況不規(guī)則很多,說明臺風風場受到了局部地形和土地覆蓋的修正,使風場更趨于“本地化”。同樣的,隨著高度的增加,下墊面對風場的影響逐步減弱,邊界層臺風風場結構逐漸向梯度風場結構靠近。
為驗證改進的臺風風場模型,將模擬的3個歷史臺風風場與位于廣東省內的氣象站觀測結果做了對比,如圖5~圖7所示。對于山竹臺風,對比了其在北京時間2018年9月16日19:00的風場 (與圖4相同),選取的測站點位于臺風中心500 km范圍內,共1 077個[圖5(a)]。對于臺風彩虹,對比了北京時間2015年10月4日15:00的風場,臺風中心位于 (21.3°N,110.3°E),臺風中心氣壓為 945 hPa,臺風移動速度為 8.45 m/s,移動方向為北偏西45°;其中位于臺風中心500 km范圍內的測站點有52個[圖6(a)]。對于臺風威馬遜,對比了北京時間2014年7月18日20:00的風場,其中心位置為 (20.3°N,110.3°E),中心氣壓為 910 hPa,移動速度為5.25 m/s,移動方向為北偏西68.2°;位于臺風中心500 km范圍內的測站點有41個[圖7(a)]。需要說明的是,這3個時刻的臺風均已登陸,并采用圖2、圖3中的地形和地表粗糙長度。
圖5 臺風山竹的風速對比
圖5(b)、圖6(b)、圖7(b)為各臺風的觀測與模擬結果的比較,其中臺風山竹的相關系數(shù)R為0.4,與實測值的相關程度較低;臺風彩虹的相關系數(shù)R為0.75,表現(xiàn)為中等程度相關;臺風威馬遜的相關程度最高為0.92。均方根誤差(RMSE)表現(xiàn)在7 m/s左右,約為最大風速的15%。另外,3個臺風時刻的散點較為均勻分布在45°線左右,沒有出現(xiàn)明顯的系統(tǒng)偏差。
圖6 臺風彩虹的風速對比
圖7 臺風威馬遜的風速對比
臺風風速和風向的時程對比,成為驗證臺風模型的重要環(huán)節(jié)。利用兩個觀測站點的風速風向數(shù)據(jù),分別對比了臺風彩虹和臺風威馬遜在登陸過程中的風速風向。圖8為這兩場臺風的移動路徑和測站點的位置。臺風的中心點位置和中心氣壓數(shù)據(jù)來源于中國氣象局發(fā)布的最佳路徑數(shù)據(jù)集[23],地形和地表粗糙長度同圖2和圖3保持一致。
圖8 臺風移動路徑和測站點位置
圖9(2015年10月4—5日)和圖10(2014年7月18—19日)分別為臺風彩虹和臺風威馬遜的兩個測站點的風速風向圖。結果顯示本文模型模擬的風速風向時程與觀測結果較為吻合,說明采用本文方法模擬臺風風場具有一定的可靠性。
圖9 臺風彩虹模擬與觀測的風速風向時程對比
圖10 臺風威馬遜模擬與觀測的風速風向時程對比
在臺風登陸過程中,地形和土地覆蓋的不均勻性會對臺風風場造成顯著影響。平坦下墊面和單一地表粗糙長度假設已經(jīng)不適用于登陸臺風的模擬。為此,將目前工程上常用的參數(shù)化臺風風場模型建立在隨地形變化的σ坐標系下,并通過土地覆蓋類型和地表粗糙長度的映射關系,考慮了臺風登陸過程中土地覆蓋變化的影響。為驗證改進后臺風風場模型的可靠性,針對影響廣東省地區(qū)的三場歷史臺風[臺風山竹 (2018)、臺風彩虹 (2015)和臺風威馬遜 (2014)]做了研究和對比。得出如下結論。
(1) 地形和土地覆蓋對近地表的臺風風場有很大影響。隨著離地高度的增加,地形和土地覆蓋影響逐漸減弱,邊界層臺風風場逐漸和梯度高度風場趨于一致。
(2) 臺風登陸以后,相同位置處模擬的風速與觀測的風速散點分布接近45°直線,說明考慮地形起伏和土地覆蓋變化后的臺風風場模型并沒有表現(xiàn)出明顯的系統(tǒng)偏差。
(3)臺風登陸過程中,4個測站點的風速風向對比結果表明,本文模型能夠較好地模擬出臺風的風速風向隨時間的變化情況。