雅魯藏布江中游風沙區典型下墊面空氣動力學參數研究

張正偲， 張 焱,2， 馬鵬飛， 潘凱佳,2， 扎 多，益西拉姆， 仁青桑布

（1.中國科學院西北生態環境資源研究院沙漠與沙漠化重點實驗室，甘肅 蘭州 730000；2.中國科學院大學，北京 100049；3.西藏自治區氣候中心，西藏 拉薩 850000；4. 山南市氣象局，西藏 山南 856000）

風沙活動是影響西北干旱、半干旱區生態環境與社會可持續發展的主要因子之一［1］。青藏高原地區同樣存在風沙活動，因高海拔、低氣壓的大氣環境導致其風沙活動過程與西北低海拔地區明顯不同［2］。

空氣動力學參數摩阻風速（u*）和空氣動力學粗糙度（z0）是風沙活動研究的重要內容［3-4］。正確計算u*和z0，有利于預測和量化輸沙量和粉塵通量［5-6］、土壤風蝕潛力［7］，是評價防沙治沙工程的重要指標［8］。u*是氣流對地表剪切力的反映，當摩阻風速超過臨界摩阻風速時，地表微粒脫離地面進入空中，發生風蝕［5］。z0是地表風速減小到零時的某一高度，是大氣邊界層湍流屬性通量參數化方案中常用的參數之一［7］。z0對地表粗糙元變化的響應敏感，反映地表粗糙元對風速的減弱作用［7］。u*、z0受下墊面性質影響顯著。茅宇豪等［9］對不同干擾程度的草原以及沙丘等11種下墊面進行空氣動力學參數研究，發現沙丘地表u*和z0最小，且u*和z0與下墊面的生物量以及起伏程度密切相關。何玉斐等［10］認為，戈壁地表粗糙度為10-4～10-3m。周杰等［11］研究了策勒縣高風速事件中沙漠-綠洲過渡帶4種下墊面的u*、z0，結果表明u*和z0流動沙丘＞半流動沙丘＞固定沙丘＞綠洲。上述研究表明，不同下墊面空氣動力學性質差異顯著。因此，對風沙區不同下墊面近地層氣流動力學參數的觀測與計算有利于認識不同粗糙元與風動力之間的能量轉換關系，進而評估區域風沙活動強度、估算沙塵釋放量等。然而，上述研究主要是針對低海拔的風沙區，高海拔地區空氣密度（0.84 kg·m-3）小于低海拔地區（1.23 kg·m-3）；高海拔地區風沙沉積物富含粉沙和黏土（5.8%±7.2%）［12］，遠高于低海拔地區（＜1%）。根據Shao 等［13］的起動風速計算公式，高海拔地區沙粒的起動過程不同于低海拔地區。但目前關于高海拔地區的野外實測研究很少，限制了我們對高海拔地區風沙災害形成機理的認知。

雅江中游風沙活動強烈［14］，是除中國北方荒漠區外的重要粉塵源區之一［15］。其風旱同期的氣候［14］和廣布的松散沉積物［16-17］以及人為破壞地表植被與土層結構等［18］，造成該區域風沙災害嚴重，但對該區域風沙問題的研究相對薄弱［16］。目前，雅江風沙研究主要集中在以下幾方面：（1）沉積物物源。研究認為河谷底部（河心洲、河漫灘）是發生風沙活動的物源［2,12,17-19］。（2）風沙地貌。風沙地貌類型復雜多樣［2,14］，但最為典型的是爬坡沙丘。（3）沙漠化。該區域土地沙漠化是氣候變化與人為活動共同作用的結果［12,19］。（4）風沙運動特征。不同地表風沙運動特征不同［20］。但上述研究不能很好闡明不同地表的沙塵運動機理，主要原因是對不同下墊面的空氣動力學性質的研究較少。因此，本文通過野外實測的風廓線數據計算了植被區、人為干擾區和流沙區3 種下墊面的空氣動力學參數，為區域輸沙量計算、粉塵通量及區域風沙災害評估提供參考。

1 研究區概況

研究區位于雅魯藏布江（雅江）中游山南段，西起貢嘎縣，東至乃東區（90°30′～92°40′E、29°10′～29°30′N），海拔高于3550 m，且自西向東逐漸降低（圖1），北側為念青唐古拉山，南側為喜馬拉雅山。研究區年平均氣溫為8.7 ℃，年平均降水量為378 mm，年內變化曲線呈單峰型，旱（冬、春季）、雨（夏）季分明；年平均相對濕度為42%，夏季較濕、冬季較干；年平均風速為2.6 m·s-1，冬春季大，夏秋季??；風向復雜，受地形影響顯著，主風向為西風［16］。

雅江中游水系辮狀發育，受徑流季節性影響，冬、春季水位低，河心洲及河漫灘大面積裸露［16-17］，表層干燥松散的沉積物在強氣流作用下向北面山麓甚至山坡輸運，風沙活動頻繁［16］。有關研究表明，該地區的風沙災害可能始于以青稞農耕為代表的歷史時期，并隨著人類開發強度的加大，風沙災害日趨嚴重；近年來，風沙活動頻次有增多的趨勢［19］。

2 研究方法

2.1 野外測量

研究區冬、春季風沙活動頻繁，尤以正午之后更為嚴重，故選擇2021 年1 月26 日14:35—18:04 和1月27日12:35—15:34分別位于山南市貢嘎縣森布日（圖1b）和昌果附近的河漫灘處（圖1c），利用二維超聲波風速儀（風速量程分別為0～30 m·s-1，分辨率為0.01 m·s-1；風向量程為0～359°；分辨率為1°；數據記錄頻率為1 min）對不同高度風速進行觀測。觀測期間，兩地均形成明顯的風沙流。

圖1 研究區地理位置及各站點地表狀況Fig.1 Location of study area and sites and land surface properties in Semburi and Changguo,respectively

森布日地表有植被覆蓋，以蘆葦（Phragmites australias）為主，植被高度為20～40 cm，植被覆蓋度為10%～30%（圖1b）。昌果地表為松散的風沙沉積物，無植被覆蓋（圖1c）。2 個測點地表起伏程度均較小。在森布日選擇植被區（S1，植被蓋度＞20%）和人為干擾區（S2，植被蓋度＜20%）2種下墊面進行測量；設有5 層高度，風速傳感器安裝高度分別為8 cm、20 cm、40 cm、90 cm和180 cm。雖然，昌果下墊面性質相近，但為了對比，同樣測量了2 組風廓線（C1，C2；均無植被），且風速傳感器安裝高度分別為5 cm、10 cm、30 cm、80 cm、170 cm。兩地參考風速的風速傳感器架均設在距地面200 cm 高度處。森布日和昌果分別記錄210條和180條數據。

2.2 數據分析

高風速條件下大氣層結呈中性［5］，中性層結條件下近地層風廓線方程可以用下式［3］表示：

式中：uz表示高度z處的風速（m·s-1）；u*表示地表摩阻風速（m·s-1）；z表示高度（m）；z0表示空氣動力學粗糙長度（m）；k是馮·卡曼常數（k=0.4）。

風速與高度的關系可以按照下式［7］利用最小二乘法進行回歸擬合。同時，本文引入判定系數（R2）和均方根誤差（RMSE）以評價擬合優度。

式中：a，b為擬合系數。

u*和z0可分別用下式計算：

在用風廓線方程計算參數時，為了控制數據質量，選擇2 m 高度風速≥4 m·s-1且風向為200°～290°范圍內的數據（圖2 中陰影區域）。其中，森布日有158 條≥4 m·s-1的風速記錄（平均風速為7.20±2.17 m·s-1，主風向為248°±21°），昌果有173 條≥4 m·s-1的風速記錄（平均風速為7.98±1.56 m·s-1，主風向為240°±14°）。

圖2 觀測期間風速和風向Fig.2 Wind speed and direction during field experiments

為了對比不同下墊面的風速廓線，將不同高度的風速和最高處的風速按照下式進行歸一化處理：

式中：i表示從低到高第i層（i=1，2，3，4，5）；ui是第i層風速（m·s-1）；ui′是第i層風速歸一化之后的無量綱風速；uf是參考風速（m·s-1）。

3 結果與分析

3.1 風速廓線特征

植被區和流沙區風速隨高度降低而降低，但風速降幅差異顯著（表1）。各站點風速降幅為：植被區＞人為干擾區＞流沙區（C1、C2），分別為：52%±21%、40%±19%、18%±12%、18%±13%。

表1 森布日和昌果不同高度的平均風速及風速降幅Tab.1 Averaged five-layer wind speeds and its reduction amplitude at Semburi and Changguo

2種下墊面的風速廓線均符合對數線性函數關系（uz=bln(z) +a，R2＞0.93，RMSE 為0.41～0.52，圖3a）。下墊面性質不同擬合函數的參數明顯不同。風速廓線的斜率越大，風速隨高度降低的幅度越大。植被區相對人為干擾區植被高度高、覆蓋度大，故S1點斜率（1.83）大于S2點斜率（1.47）。昌果地表均一，均為流沙，故昌果兩點斜率大小相近（0.91、0.87）。歸一化的風速廓線同樣表現出相似的規律（圖3b）：植被高度越高、覆蓋度越大，則斜率越大，S1 點和S2 點斜率分別為0.28、0.23。C1 點和C2 點斜率均為0.11，說明本研究的結果具有可重復性，風速廓線系數的差異是由下墊面屬性不同造成的。

圖3 風速廓線Fig.3 Wind profile

3.2 空氣動力學參數特征

下墊面性質不同，u*差異顯著。u*在植被區＞人為干擾區＞流沙區（圖4a，圖4b）。植被區u*約為人為干擾區的1.2 倍，約是流沙區的2.1 倍。植被區u*的分布范圍廣，呈單峰分布，存在粗尾（圖4a）。S1點的u*最小值、平均值和最大值分別是0.20 m·s-1、0.73 m·s-1和1.76 m·s-1（圖4b）。S2 點的u*最小值、平均值和最大值分別是0.19 m·s-1、0.59 m·s-1和1.54 m·s-1。流沙區u*的分布范圍比植被區的窄，u*分布集中。C1點的u*最小值、平均值和最大值分別是0.19 m·s-1、0.35 m·s-1和0.69 m·s-1。C2 點的u*最小值、平均值和最大值分別是0.19 m·s-1、0.36 m·s-1和0.67 m·s-1。

下墊面性質不同，z0差異亦顯著（圖4c，圖4d）。z0在植被區＞人為干擾區＞流沙區，且植被區z0約是人為干擾區的2倍，約是流沙區的150倍。植被區z0的變化區間大，概率密度曲線存在細尾（圖4c）。S1的z0變化區間為：0.5×10-2～13.1×10-2m（圖4d），平均值為7.1×10-2m；人為干擾區（S2）z0變化區間為：0.1×10-2～10.4×10-2m，平均值為3.7×10-2m。流沙區z0的變化區間比植被區小，z0分布較植被區集中，近似單峰分布。C1的z0變化區間為：1.1×10-4～39.3×10-4m，平均值為3.2×10-4m；C2 的z0變化區間為：1.2×10-4～47.8×10-4m，平均值為3.5×10-4m。

圖4 摩阻風速（u*）和空氣動力學粗糙度（z0）概率密度Fig.4 Probability density curves of u*and z0

4 討論

4.1 空氣動力學參數與下墊面性質的關系

空氣動力學參數受地表屬性和氣流狀態共同影響，其變化是粗糙元與氣流相互作用的結果［7］。為了闡明下墊面對z0和u*的影響，本文計算參數時不僅嚴格限制了風速和風向（圖2），還對風速進行了歸一化處理（圖3b）。結果表明，植被區和流沙區的空氣動力學參數差異顯著（P＜0.05），反映了下墊面性質對氣流狀態影響顯著。植被區和流沙區地表主要差異在于植被蓋度和植被高度，植被區植被蓋度約為10%～30%，植被高度約為20～40 cm，流沙區植被蓋度為0（圖1c）；其實質反映了地表粗糙元的不同。Lettau等［21］認為影響地表空氣動力學參數大小的原因很復雜，但最關鍵的決定因素是下墊面的粗糙元。Li等［22］得出相似結論，認為地表粗糙元是影響空氣動力學參數的關鍵因子，地表粗糙元越大，u*和z0愈大。植被的存在增加了地表粗糙元，增強了近地層氣流的雷諾剪切引力且造成零平面高度上升［23-24］，故植被區的u*和z0大于流沙區。如S1植被蓋度和高度明顯大于S2，導致S1 的z0是S2 的1.9 倍。C1 和C2 為無植被，所以z0平均值變化不明顯，分別為3.2×10-4m 和3.5×10-4m。植被區的頂層和底層風速降幅可達89%，而流沙區僅為40%。這種差異導致植被區的u*明顯大于流沙區（S1和C1的u*平均值分別為0.73 m·s-1和0.35 m·s-1）。張強等［25］認為造成這種差異的主要原因是柔性植被可以通過改變自身的結構和形態更多地消耗近地表氣流的能量，減少了直接到達地表的能量，從而增大了近地表風速的降幅，使得近地表風速快速降至0，導致u*和z0增大。

4.2 空氣動力學參數與平均風速的關系

4.2.1 摩阻風速與平均風速的關系 植被區和流沙區的u*與2 m高度平均風速呈顯著的線性正相關關系（P＜0.05），這與前人的研究結果相似［9,26］，但2 種下墊面的u*隨風速增大程度不同（圖5）。植被區u*隨風速增大程度是流沙區的2倍（S2和C2斜率分別為0.12和0.06）。此外，S1增加幅度略大于S2（斜率分別為0.13 和0.12），C1 與C2 地表均是由平坦的流沙組成，下墊面相似，兩地u*隨風速增大的程度相近。據此判斷，u*與風速的增大程度與地表類型密切相關，植被覆蓋的地表u*隨風速增大程度比流沙

圖5 摩阻風速（u*）與平均風速（u）的關系Fig.5 Relationship between threshold friction velocity(u*)and averaged wind speed(u)

地表大；近地層氣流湍流強度增強是地表摩阻風速增大的主要原因。在相同地表，近地層湍流強度隨著風速增大而增大，導致地表摩阻風速增大，表現出摩阻風速隨風速增大而增大。除此之外，筆者還觀察到植被區摩阻風速隨風速的增大程度比流沙區大；其原因是粗糙度大的植被區增強了近地層氣流的湍流強度，且隨風速增大而增強。

4.2.2 空氣動力學粗糙度與平均風速的關系 植被區和流沙區z0與2 m高度平均風速的關系明顯不同（圖6）。植被區的z0隨風速增大而緩慢增大（圖6a）。流沙區的z0隨風速增大先緩慢減小，當風速增大到閾值（6.5 m·s-1）后，z0迅速增大（圖6b）。張強等［25］通過研究不同植被下墊面的空力動力學參數時發現，近地層風速會明顯改變柔性粗糙元的結構和形態，從而導致z0在低風速時較大，高風速時較小。楊興華等［27］在塔克拉瑪干沙漠過渡帶研究了無風沙運動的床面表明，大氣層結呈中性時，z0與2 m 高度處風速呈顯著的負指數關系。Zhang 等［28］在騰格里沙漠東南緣觀測風沙活動條件下，草方格覆蓋的地表和沙質地表的空氣動力學參數時發現，低風速（＜6 m·s-1）時草方格覆蓋的地表z0隨風速迅速增大，高風速（＞6 m·s-1）時緩慢增大；但沙質地表z0在低風速（＜8 m·s-1）時隨風速增大而緩慢增大，在高風速（＞8 m·s-1）時z0迅速增大。高詠晴等［29］通過風洞實驗研究了風沙流和凈風場中空氣動力學參數發現，風沙流中z0大于凈風場，且z0在凈風場中隨著風速增大而減小，但在風沙流中z0隨著風速的增大而增大。王翠等［30］的研究也證明了這一點。綜上所述，筆者認為：（1）植被區z0是風沙流和植被綜合作用的結果。（2）流沙區有無風沙流對計算的粗糙度值影響顯著。無風沙流時計算的z0可能會偏大，形成風沙流時，z0隨風速呈指數增加。而臨界值（6.5 m·s-1）指示了起動風速，這與計算的臨界起沙風風速相似（流沙地表平均粒徑為292.81 μm。沙粒密度為2650 kg·m-3。觀測期間平均氣壓為655.7 hPa，平均氣溫為10.1 ℃。根據氣體狀態方程計算出空氣密度約為0.81 kg·m-3。據此，利用起沙風公式［3］計算流沙地表2 m 高度處臨界起沙風風速為6.6 m·s-1）。

圖6 平均風速（u）與空氣動力學粗糙度（z0）的關系Fig.6 Relationship between aerodynamic roughness length(z0)and averaged wind speed(u)

4.3 海拔對風速廓線的影響

為了闡明海拔與近地層風速的關系，對不同區域流沙地表的野外實測資料按照公式（5）做歸一化處理（圖7）。研究區流沙地表（昌果）與塔中［31］、敦煌［32］、庫布齊風電場［33］、巴丹吉林拐子湖［34］、中衛風沙觀測場［35］的歸一化風速廓線均符合對數線形規律（u′=cln(z)+d；c，d為標準化風速廓線擬合系數，R2＞0.98，RMSE 為0.03～0.06）（圖7a、表2）；參數c隨海拔增加呈指數遞減，但d呈指數遞增，揭示了近地層風速隨海拔高度變化的差異（圖7b、圖7c）；c值越大（d值越?。?，風速減小越快，反之則亦然。這表明相對于低海拔地區，高海拔地區空氣密度小，黏度小，氣流內摩擦力小，氣流能量向下傳遞時被消耗的少，根據雷諾數公式判斷高海拔地區更易于發生湍流。說明高海拔地區氣流隨高度增加（減?。┏潭纫鹊秃０蔚貐^慢，能量消耗較少，從而導致高海拔地區風沙沉積物比低海拔地區運動的更高、更遠［20］。

表2 歸一化風速廓線參數（u′=cln(z)+d）Tab.2 Normalized wind speed profile parameters for each station

圖7 歸一化風速廓線，風速廓線參數（c,d）與海拔（h）的關系Fig.7 Normalized wind speed profiles of naturally flat bare sand surfaces at high and low altitudes areas and coefficient c and d versus altitude

綜上所述，森布日和昌果的u*和z0差異顯著，與平均風速的關系明顯不同。盡管在高風速條件下，森布日和昌果地表都有可能發生風蝕，但通過對森布日和昌果u*和z0的分析，昌果地表發生風蝕的潛力更大。結合研究區的地表類型時空變化特征，冬、春季雅江中游類似昌果地表環境的區域，河漫灘、河心洲大面積裸露［17］，建議風沙災害防治重點放在河心洲和河漫灘。

5 結論

（1）植被區近地層風速的降低幅度大于流沙區，體現了植被對近地層風速的減弱作用，從而減弱土壤風蝕。

（2）植被區的u*、z0大于流沙區，植被區的u*、z0分別約是流沙區的2倍、150倍。植被區和流沙區的u*平均值分別為0.73 m·s-1和0.36 m·s-1，對應z0平均值分別在10-2、10-4量級。據此，風沙災害防治重點應該在河心洲和河漫灘。

（3）植被區與流沙區的摩阻風速均與平均風速呈顯著的線形關系，隨平均風速增大而增大，但植被區增大的程度大于流沙區。

（4）植被區與流沙區的空氣動力學粗糙度對平均風速的響應方式不同，與下墊面粗糙元性質以及是否形成風沙流密切相關。