雅鲁藏布江中游风沙区典型下垫面空气动力学参数研究

张正偲， 张 焱,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）植被区与流沙区的空气动力学粗糙度对平均风速的响应方式不同，与下垫面粗糙元性质以及是否形成风沙流密切相关。