饒志娟 朱彬 汪如良 劉熙明

(1.江西省氣象學會秘書處，江西 南昌 330096； 2.南京信息工程大學環(huán)境科學與工程學院，江蘇 南京 210044；3.南京信息工程大學大氣物理學院，江蘇 南京 210044； 4.江西省氣象服務(wù)中心，江西 南昌 330096；5.江西省氣象科學研究所，江西 南昌 330096)

引言

大氣邊界層是人類生產(chǎn)、生活的主要空間，人類活動通過改變自然地表特征影響邊界層內(nèi)大氣湍流運動，并對天氣氣候產(chǎn)生影響[1]。地表類型復雜且非均勻的下墊面不僅使得邊界層湍流運動多變，而且對數(shù)值模式的模擬能力提出了挑戰(zhàn)。各類大氣模式的模擬效果不僅取決于氣象探測數(shù)據(jù)的精度還取決于對下墊面性質(zhì)的了解和湍流特征量的表征。

空氣動力學粗糙度是大氣邊界層的一個重要概念，反映了地表狀態(tài)的物理特性?？諝鈩恿Υ植诙茸钤缬葿agnold[2]提出，他通過大量的風洞實驗研究氣流與沙粒的相互作用，并得到了著名的空氣動力粗糙度1/3定律。早期的研究基本集中在農(nóng)田、海面、森林或城郊等起伏不太大的下表面[3-4]，Baldauf等[5]在中性層結(jié)和地形平坦的條件下，對任意給定的粗糙度分布z0(x，y)，提出了一種參數(shù)化的有效粗糙度長度計算方法，他考慮了空間結(jié)構(gòu)分布，尤其是長度尺度的影響和氣流流入的方向。近年來開展的復雜地形大氣湍流結(jié)構(gòu)的研究更具有現(xiàn)實意義。Zhong等[6]在陸面模式中發(fā)展了一種計算非均勻地表的有效粗糙度和零平面位移的新方案。黃鶴等[7]研究了渤海西岸空氣動力學粗糙度特征，對比分析了風廓線法和風速標準差法的計算結(jié)果。沙敏敏等[8]利用多層風速、風向觀測資料計算了北京市空氣動力學粗糙度，發(fā)現(xiàn)地表粗糙度、零平面位移隨著城市化發(fā)展皆有明顯增加的趨勢。

鄱陽湖是我國第一大淡水湖，是國際重要濕地，湖區(qū)呈盆狀凹地，周邊地形亦復雜，具有典型的非均勻下墊面特征。研究鄱陽湖地區(qū)非均勻下墊面地表粗糙度特征，對于局地天氣預報、湖濱地區(qū)風能資源利用、空氣污染物擴散、邊界層特征數(shù)值模擬等具有重要意義[9-10]。本文利用2013年7月1日至2014年6月30日鄱陽湖東岸70 m鐵塔渦動相關(guān)觀測資料，計算零平面位移及粗糙度，并為本地化陸面模式提供參考，以期提高模式輸出效果和預報準確性。

1 資料與方法

70 m觀測塔位于鄱陽湖東岸的江西省鄱陽縣白沙洲鄉(xiāng)境內(nèi)(29.16°N，116.61°E)，距離鄱陽縣氣象基準站18.0 km，鐵塔周邊地形平坦，下墊面主要為低矮灌木。觀驗采用的EC150開路渦動相關(guān)觀測儀安裝在鐵塔44.0 m處，超聲風速儀(CAST3，Campbell)數(shù)據(jù)采集器為CR3000 Campbell，采集數(shù)據(jù)包括風速、風向、溫度、濕度以及儀器輸出的摩擦速度，采樣頻率為10 Hz，觀測時間為2013年7月1日至2014年6月30日；降水數(shù)據(jù)來源于鄱陽縣氣象基準站。

1.1 數(shù)據(jù)處理

對原始數(shù)據(jù)剔除降水日數(shù)據(jù)、去野點、去趨勢、坐標旋轉(zhuǎn)，重新計算30 min時間長度數(shù)據(jù)。

整體湍流特征(Internal Turulence Characteristics，簡稱ITC)是檢驗湍流發(fā)展的充分性的重要參數(shù)，通常認為ITC<30%時，湍流充分發(fā)展[11]。ITC的計算式為：

(1)

式(1)中，下標“ec”和“mod”分別為觀測值和莫寧—奧布霍夫相似理論得到的計算值；σw為垂直風速標準差；u*為摩擦速度。Thomas和Foken[12]給出：

(2)

根據(jù)Gokede等[13]給出的湍流數(shù)據(jù)質(zhì)量分級指標(表1)和Irwin[14]給出的L穩(wěn)定度分類標準(表2)，QC為1—3表示質(zhì)量好的數(shù)據(jù)；4—6代表質(zhì)量較好的數(shù)據(jù)；7—8代表質(zhì)量差的數(shù)據(jù)；9表示需要剔除的數(shù)據(jù)；A、B、C、D、E、F這6類分別對應(yīng)了極不穩(wěn)定、不穩(wěn)定、弱不穩(wěn)定、中性、弱穩(wěn)定、穩(wěn)定。文中對湍流觀測資料進行檢驗(圖1)，其中79.3%為質(zhì)量好的數(shù)據(jù)(QC為1—3)，97.03%為質(zhì)量較好數(shù)據(jù)(QC為1—6)，0.63%為需要剔除的數(shù)據(jù)。根據(jù)Gokede等[13]和Irwin[14]的分類并結(jié)合當?shù)氐膶嶋H情況，文中篩選出質(zhì)量等級為1—3且為中性層結(jié)(L<-180.91或L>84.25)的數(shù)據(jù)。

表1 Gokede等(2004)給出的湍流數(shù)據(jù)質(zhì)量分級指標[13]Table 1 Turbulence data quality classification index given by Gokede et al.(2004)

表2 Irwin(1979)給出的的L穩(wěn)定度分類標準[14]Table 2 Classification criteria of Lstability given by Irwin (1979)

圖1 湍流發(fā)展的充分性與穩(wěn)定度分類統(tǒng)計Fig.1 Classification statistics of the sufficiency and stability for the turbulence development

Taylor[15]提出在滿足某些條件的情況下，當湍流流經(jīng)傳感器時，可以認為湍流是被凍結(jié)的。其含義是在空間上一固定點對湍流的觀測結(jié)果統(tǒng)計上等同于同時段沿平均風方向空間各點的觀測，也稱為定型湍流假設(shè)，當然湍流并不是真的被凍結(jié)，只是假設(shè)湍渦發(fā)展的時間尺度大于它被平流攜帶經(jīng)過探頭所需的時間，泰勒假設(shè)才適用。此假設(shè)中實際隱含著平穩(wěn)氣流和均勻湍流的條件，風速也不宜過小。為此，Willis 和 Deardorrff[16]提出以風速標準差(即湍流強度)作為使用假說的依據(jù)。湍流強度是表征湍流發(fā)展強度的量，即風速的標準偏差與平均風速之比。計算式為：

(3)

式(3)中，I為湍流強度；σ為風速標準差；u為平均風速。

圖2給出了u、v、w方向的湍流強度與風速的關(guān)系。總體而言，Iu>Iv>Iw，垂直方向上湍流較弱。當風速>2 m·s-1時，湍流強度基本在0.4以下，且隨著風速的增加，湍流強度減小，滿足泰勒假設(shè)。當風速<2 m·s-1時，湍流強度變化值在0—1，呈顯著的離散分布。隨著風速的增大，湍流強度迅速減小，約在風速>2 m·s-1時趨于穩(wěn)定。由此可知，風速較大時，觀測的相對誤差較小，湍流數(shù)據(jù)的穩(wěn)定性也較強。因此，將風速<2 m·s-1的數(shù)據(jù)剔除。

圖2 u(a)、v(b)、w(c)方向的湍流強度與風速的關(guān)系Fig.2 Relationships between turbulence intensity and wind speed in u (a)、v (b) and w (c) directions

風速、溫度和濕度的歸一化標準差也稱為湍流統(tǒng)計特征[17]，表征了對所有頻率湍流信號進行統(tǒng)計分析的結(jié)果。

圖3 歸一化的垂直速度標準差(a)和溫度標準差(b)Fig.3 Standard deviations of normalized vertical velocity (a) and temperature (b)

1.2 研究方法

零平面位移作為描述下墊面空氣動力學特征的重要參數(shù)，常常用于邊界層的湍流研究中，也是邊界層參數(shù)化的重要參數(shù)之一。確定零平面位移的傳統(tǒng)方法是利用中性近地面層風廓線進行計算，如果獲得了3個或3個以上高度的風速梯度觀測就可以用曲線回歸求解。但是實際上該方法操作起來很難，因為首先中性層結(jié)資料的獲取是十分困難的；其次也比較難于得到多個觀測高度準確風速值的廓線。超聲風溫儀由于其測量精度高、性能穩(wěn)定可靠、設(shè)計靈活等特點，是測量大氣湍流數(shù)據(jù)的有效手段之一。

1.2.1 Martano方法

根據(jù)Monin-Obukhov相似性理論，風速廓線可寫為：

u=(u*/k){ln[(z-d)/z0]-ψm(ζ)}

(4)

式(4)中，u*為摩擦速度；d為零平面位移；z0為動力粗糙度；k為Karman常數(shù)；ζ為穩(wěn)定度；ψm(ζ)為集成了穩(wěn)定度修正的函數(shù)。

文中采用Businger-Dyer-Webb[18]給出的經(jīng)驗公式。

(1)不穩(wěn)定層結(jié)下，ζ=z/L<0

(5)

其中，x=(1-16ζ)1/4。

(2)穩(wěn)定層結(jié)下，ζ=z/L>0

ψm(ζ)=-5ζ

(6)

同一高度上有多組平均量時，可用最小二乘擬合計算d和z0，計算式為：

〈{ku/u*-ln[(z-d)/z0]+ψm(ζ)}2〉=min (z0，d)

(7)

這是一個二元的非線性最小二乘問題，式(7)可以寫為:

〈[S(z0，d)-p(z0，d)]2〉=min(z0，d)

(8)

其中，p(z0，d)=ln(z-d)/z0為參數(shù)。

當滿足條件

〈S(z0，d)〉-p(z0，d)=0

(9)

或

(10)

式(8)可寫為:

〈[S(z0，d)-p(z0，d)-〈S(z0，d)-p(z0，d)〉]2〉=min (z0，d)

(11)

由于p為參數(shù)，〈p〉=p，式(11)可寫為:

(12)

z0e=(z-d)exp〈-S〉?(z-d)〈exp(-S)〉=〈z0〉m

(13)

1.2.2 TVM方法

Rotach[19]提出的通過測量溫度脈動方差從而確定零平面位移和地表粗糙度的溫度方差法，稱為TVM法(Temperature Variance Method)。根據(jù)相似理論，在不穩(wěn)定層結(jié)下，有下式：

(14)

在求得零平面位移和動力粗糙度后，可結(jié)合觀測得到的u、L，利用式(4)計算摩擦速度u*與觀測u*進行比較，其中相似性函數(shù)采用Businger-Dyer-Webb[18]給出的經(jīng)驗公式。

2 結(jié)果分析

2.1 季節(jié)變化

鄱陽湖地區(qū)四季分明，水域面積隨季節(jié)變化明顯，夏季最大、冬季最小，周邊植被在春、夏季生長旺盛，秋、冬季節(jié)凋零。因此，文中將時間劃分為春夏(3—8月)和秋冬(9月至翌年2月)兩段，對零平面位移d進行擬合。

圖4a和圖4b給出了Martano方法計算的春夏和秋冬S的標準差σS隨零平面位移d的變化，最小的σS對應(yīng)的d為最佳值。分析可知，春夏季零平面位移d為2.6 m，對應(yīng)的粗糙度z0為0.0145 m；秋冬季零平面位移d為4.3 m，對應(yīng)的粗糙度z0為0.0023 m?？梢钥闯?，春夏季零平面位移要小于秋冬季，但粗糙度遠大于秋冬季，為秋冬季的6倍多。這是因為春夏季鄱陽湖地區(qū)植被生長旺盛，但植被高度要小于秋冬季?？梢姷乇泶植诙扰c下墊面的物理性質(zhì)密切相關(guān)。

圖4c和圖4d給出了TVM方法計算的春夏和秋冬S的標準差σS隨零平面位移d的變化。分析可知，春夏季和秋冬季的零平面位移d分別為8.0 m和7.5 m，對應(yīng)的粗糙度z0分別為0.0370 m和0.0166 m，零平面位移d隨季節(jié)的變化不大，但春夏季的粗糙度為秋冬季的2倍多。

圖4 Martano方法計算的春夏(a)、秋冬(b)和TVM方法計算的春夏(c)、秋冬(d)S的標準差σS隨零平面位移d變化Fig.4 The variations of the standard deviation σS of S in spring and summer (a),autumn and winter (b) calculated with Martano method and the corresponding values in spring and summer (c),autumn and winter (d) by TVM method with zero plane displacement d

2.2 湖面和陸面的變化

由于下墊面的不均勻，在觀測塔的西邊為湖面，其他方向為低矮灌木叢。文中定義270°—360°方向為來自湖面的通量，其他方向為來自陸面通量。圖5a和圖5b給出了Martano方法計算的來自湖面和陸面S的標準差σS隨零平面位移d的變化。分析可知，湖面方向的零平面位移d為4.2 m，對應(yīng)的粗糙度z0為0.0012 m；陸面方向的零平面位移d為8.6 m，對應(yīng)的粗糙度z0為0.0122 m?？梢钥闯觯懨娣较虻拇植诙仁呛娴?0倍。

圖5 Martano方法計算的來自湖面(a)、陸面(b)和TVM方法計算的來自湖面(c)、陸面(d)S的標準差σS隨零平面位移d的變化Fig.5 The variations of the standard deviation σS of S in spring and summer (a),autumn and winter (b) calculated with TVM method and the corresponding values in lake surface (c),land surface (d) by TVM method with zero plane displacement d

圖5c和圖5d給出了TVM方法計算的來自湖面和陸面S的標準差σS隨零平面位移d的變化。分析可知，湖面方向的零平面位移d為3.2 m，對應(yīng)的粗糙度z0為0.0034 m；陸面方向的零平面位移d為4.4 m，對應(yīng)的粗糙度z0為0.0132 m。

2.3 兩種方法計算結(jié)果對比分析

表3給出了Martano和TVM方法計算的零平面位移d和粗糙度z0。分析可知，用Martano方法計算的結(jié)果隨季節(jié)變化較大，春夏季的粗糙度z0是秋冬季的6.3倍。而用TVM方法計算的零平面位移d春夏和秋冬季相差不大，粗糙度z0春夏季是秋冬季的2.2倍。用Martano方法計算的陸面方向零平面位移d和粗糙度z0分別為來自湖面的2倍和10倍，而用TVM方法計算的結(jié)果分別為1.4倍和3.9倍?？梢钥闯?，Martano方法對季節(jié)和方向更敏感。

表3 Martano和TVM方法計算的零平面位移d和粗糙度z0Table 3 Zero plane displacement d and roughness z0 calculated with Martano and TVM methods m

分別用不同方法計算得到年平均d和z0，代入Monin-Obukhov相似性理論的風廓線關(guān)系式(4)計算摩擦速度。圖6給出了計算的摩擦速度和實測摩擦速度的線性回歸。分析可知，利用Martano方法計算的摩擦速度與實測值的擬合曲線為y=1.099x，對摩擦速度造成了約9.9%的高估。TVM方法計算的摩擦速度與實測值的擬合曲線為y=1.328x，對摩擦速度造成了約32.8%的高估。由此可見，用Martano方法計算零平面位移d和粗糙度z0代入Monin-Obukhov相似性關(guān)系中計算的摩擦速度和觀測值有較好的一致性。

圖6 Martano(a)和TVM(b)方法計算的摩擦速度和實測值的線性回歸Fig.6 The linear regressions between the measured friction velocity and that calculated with Martano (a)and TVM (b) methods

文中計算的陸面和湖面的d和z0與前人研究相比存在一定的差異。Sorbjan[20]給出的平坦下墊面z0量級在10-3m左右，Stull[21]給出的廣闊水域z0量級在10-3—10-4m，Verburg和Antenucci[22]計算Tanganyika湖(世界第三大湖)的z0量級在10-5m左右。其原因可能與湖體面積及深度有關(guān)，也與大氣穩(wěn)定度有關(guān)。

3 結(jié)論

(1)Martano方法的計算結(jié)果隨季節(jié)變化較大，春夏季的粗糙度z0是秋冬季的6.3倍。TVM方法的計算結(jié)果季節(jié)變化不大，粗糙度z0春夏季是秋冬季的2.2倍。

(2)Martano方法計算的陸面方向零平面位移d和粗糙度z0分別為來自湖面的2倍和10倍，而TVM方法計算的結(jié)果分別為1.4倍和3.9倍。

(3) Martano方法計算的摩擦速度與實測值的擬合曲線為y=1.099x，對摩擦速度造成了約9.9%的高估。TVM方法計算的摩擦速度與實測值的擬合曲線為y=1.328x，對摩擦速度造成了約32.8%的高估。用Martano方法計算的摩擦速度和觀測值有較好的一致性。