引言

大气边界层是指距离地面1~2 km的大气层最底下的一个薄层，它是大气与下垫面直接发生相互作用的层次，它与天气、气候以及大气环境研究有非常密切的关系。地面和大气间进行着动量、水汽、热量和物质交换与湍流输送，其过程直接影响和决定大气边界层的形成和发展，进一步影响气候变化和大气环境等。过去的一个世纪，地球气候经历了以变暖为主要特征的显著变化，由此而导致的海平面升高，冰川退缩，湖泊水位下降以及越来越频繁的严重自然灾害给人类生存和社会安全造成了极大的威胁(郭建侠，2006)。森林作为陆地生态系统的主体，是地球生物圈的重要组成部分。森林植被改变了下垫面的物理特征，直接影响辐射平衡、热量平衡、水量平衡和湍流通量等，形成特有的小气候特征，并与大气以湍流的形式不断进行动量、热量和水分交换(张一平等，2006；Mengesha et al，2015；谢馨瑶等，2018)，其不仅具有改善和维护区域生态环境的功能，而且在调节全球气候、维持全球碳平衡等方面具有不可替代的作用(蔡慧颖，2017)。

近些年，国内外众多学者对森林生态系统的湍流特征进行了大量研究(Pita et al，2013；Richardson et al，2010)。刘和平等(1997)通过对长白山森林生态定位研究站内气象塔上湍流及风廓线资料计算了长白山森林空气动力学粗糙度长度(z₀)，零平面位移(d)及曳力系数(C_D)、热力输送系数(C_H)，并对误差进行了分析。刘树华等(2003)利用观测的湍流资料，在对森林冠层上下的速度分量和温度谱、局地各向同性、湍流能量和热量归一化耗散率特征研究的基础上，首次给出了森林冠层上的湍涡特征长度尺度、湍流动量和热量耗散率、湍流动量和热量结构参数的特征。李萍阳等(2002)通过对美国威斯康星州北部的Cheqamegon国家森林的两个气象塔塔层湍流测量值和特征分布，计算分析得出森林的拖曳系数较以往的平坦地形大得多；与平坦地形相比，森林及林木湿地下垫面的湍流强度较大，且湍流强度随高度增加而减弱；森林的无量纲风速标准差小于林木湿地，无量纲风速标准差随高度增加而增大。刘罡等(2006)研究表明水汽和二氧化碳浓度的湍流标准差相似函数比温度要分散，并且独立于稳定度，说明在森林冠层上空的大气湍流运动中，热量、水汽和二氧化碳的输送是互不相似的，并对这种不相似现象的可能原因进行了初步的研究和探索。王春林等(2007)应用最小二乘法和能量平衡比率两种方法，系统分析了广东省鼎湖山针阔叶混交林生态系统能量平衡特点，并分析各种涡度通量修正方法对能量平衡的影响。赵晓松等(2004)利用廓线法中的牛顿迭代法计算了森林下垫面的零平面位移(d)和空气动力学粗糙度长度(z₀)，结果显示d在生长季较大、z₀较小，而在非生长季刚好相反。

小兴安岭地处北温带大陆性季风气候区，是以阔叶红松林为代表的针阔叶混交林区，是国家重点林区之一，也是中国天然林生态系统的核心区域之一，植被固碳潜力巨大，在中国森林碳汇中占有重要地位(赵溪竹，2010)。但到目前为止，国内对于小兴安岭的湍流特征尚未进行深入研究。因此，对小兴安岭地区的湍流特征进行深入、细致的研究十分必要，无论是为数值模式的发展提供基础，还是对政府决策部门提供相应依据都有重要意义。

本文利用中国气象局五营森林生态监测站内梯度观测塔上的2016年12月至2017年11月一年的涡度相关资料，研究小兴安岭森林下垫面湍流的基本特征，为提高对小兴安岭森林下垫面湍流输送特征及陆气物质与能量交换的规律和机理特征的认识提供基础。

1 观测场地、仪器简介

五营镇位于汤旺河上游地段，小兴安岭南坡腹部，全境地貌除汤旺河沿岸有少量冲积开阔地外，绝大部分属小兴安岭山地。整个地势呈东南略高、西北略低、山地居多、平原较少、地势平缓，属于低山丘陵地带。五营镇属北温带大陆性季风气候，年平均气温为0.6℃，年降水量为610.7 mm，年日照时数为2 196.0 h，无霜期为117 d，大于或等于10℃的积温为2 141.8℃，全年以西南风和偏南风为主导风向。五营原始阔叶红松林面积达3.2 hm²，有“红松故乡”之称。

观测场地为中国气象局五营红松林野外生态监测站(海拔高度为342 m)，距离五营镇西北方向10 km。站内建有面积为50 m×50 m的梯度观测场，位于次生针阔叶混交林中，平均树高约为23~25 m(图 1)。观测场附近地势略有起伏，整体表现为北高南低(图 1a)。

观测场内建有一个70 m高的梯度观测塔，50 m高度处设涡度测量系统，由数据采集器(Li-7550, Li-cor, USA)、三维超声风速仪(WindMaster, Gill, USA)、CO₂/H₂O分析仪(Li-7550, Li-cor, USA)组成，用于地气间CO₂、H₂O、动量、显热、潜热和三维风速等要素的监测及森林生态系统碳循环的研究。

2 观测资料及数据质量控制方法

以ξ=z/L表示大气稳定度，z/L>0为稳定大气，z/L < 0为不稳定大气。z为仪器测量高度，对于森林下垫面，需要引入零平面位移(d)，以高度z-d置换高度z(Kustas and Brutsaert, 1986；张宏昇，2014)；L为莫宁-奥布霍夫长度，计算公式如下：

$ L=-u_*{}^3 / k \frac{g}{\theta} \overline{w^{\prime} \theta^{\prime}} $

(1)

$ u_*=\left[\left(\overline{u^{\prime} w^{\prime}}{}^2\right)+\left(\overline{v^{\prime} w^{\prime}}{}^2\right)\right]^{1 / 4} $

(2)

式中：θ为虚位温，u_*为摩擦速度，k为卡曼常数(取0.4)，g为重力加速度，u, v, w为坐标旋转后的顺风、侧风、垂直向风速(Stull，1991)。

非中性条件下的近地层风速廓线公式为(Stull，1991)：

$ \frac{\bar{u}}{u_*}=\frac{1}{k}\left[\ln \left(\frac{z-d}{z_0}\right)+\varPsi_{\mathrm{M}}\left(\frac{z-d}{L}\right)\right] $

(3)

式中: z₀为空气动力粗糙度长度，定义为风速为零的高度。Ψ_M为稳定度函数，在稳定和不稳定条件下可分别表示为：

$ \varPsi_{\mathrm{M}}\left(\frac{z-d}{L}\right)=\frac{4.7(z-d)}{L} z / L>0 $

(4)

$ \begin{gathered} \varPsi_{\mathrm{M}}\left(\frac{z-d}{L}\right)=-2 \ln \left(\frac{1+x}{2}\right)- \\ \ln \left(\frac{1+x^2}{2}\right)+2 \tan ^{-1}\left(x-\frac{\pi}{2}\right) z / L<0 \end{gathered} $

(5)

式中: x=[1-(15z/L)]^3/4。

根据莫宁-奥布霍夫相似理论，任何一个量的标准差被特征尺度参数标准化后就被成为z/L的普适函数(Andreas et al，1998)。根据经验拟合与局地自由对流假设(张宏昇，2014)：

$ \sigma_X / X_x=\phi_X(z / L)=C_{X 1}\left(1-C_{X 2}|z / L|\right)^{\pm 1 / 3} $

(6)

式中：σ_X为X的标准差，X指u、v、w、温度T、比湿q，X_*是X的特征尺度参数(对于速度分量，X_*=w′X′/u_*；对于其他变量，X_*=-w′X′/u_*)，ϕ_X为普适函数，C_X1、C_X2为经验常数(Stull，1991)。

当-z/L→∞，边界层处于自由对流状态时，对于其他变量，式(6)可简化为：

$ \sigma_X / X_*=\phi_X(z / L)=C_X(-z / L)^{\pm 1 / 3} $

(7)

引入温度尺度(T_*=-w′T′/u_*)和湿度尺度(q_*=-w′q′/u_*)，根据式(6)和式(7)，对于无因次化温度标准差σ_T/T_*和无因次化湿度标准差σ_q/q_*，在对流状态下可以表达为：

$ \sigma_T /\left|T_*\right|=C_T(-z / L)^{-1 / 3} $

(8)

$ \sigma_q /|q *|=C_q(-z / L)^{-1 / 3} $

(9)

式中: C_T和C_q为拟合确定的参数(Stull，1991)。

观测资料选取2016年12月至2017年11月一年的数据(因仪器故障原因，9月5—9日、10月11—18日、10月30日至11月1日数据缺测)，并假定50 m位于常通量层。使用Eddypro软件对10 Hz的原始数据进行处理，具体处理流程参照王寅钧等(2013)、龚玺等(2018)、孙鹏飞等(2021)。

对Eddypro软件输出的30 min数据进行以下筛选：①剔除仪器出错时的数据；②剔除夜间弱湍流的观测数据(u_* < 0.2 m·s^-1)；③计算z₀、d、C_D时剔除弱风数据(u≤1 m·s^1-)，筛选大气稳定度范围为|z/L|≤10。

3 结果分析 3.1 风向、风速特征

研究湍流的一般方法，是把温度和风等变量分解为平均和扰动两部分。平均部分表示平均温度和平均风速的影响，扰动部分则表示波的影响或叠加在平均风速上的湍流影响。因此，了解平均场的特征是研究湍流运动及其物质和能量交换的前提。研究区域春、秋和冬季(图 2a，2c，2d)以西南风为主，其次为偏西风和偏南风；夏季以西南风为主，其次为东北风(图 2b)；总体上全年基本为偏西到偏南风，西北方向和东南方向的风很少。从图 2中的风速分布可以看到，整个观测期间以0~4 m·s^-1的风速最多，其次是4~6 m·s^-1，6~8 m·s^-1较少，8 m·s^-1以上风速接近于无，表明观测期间主要以小风环境为主，强风环境很少。

3.2 湍流强度

湍流强度(I_α=σ_α/u)指被水平风速标准化后的标准差(α=u，v，w)，σ_α为相应风速分量的标准差，u为平均风速(Stull，1991)。观测区域不同月份u, v, w方向湍流强度特征如下(表略)：10月至次年4月u, v方向湍强度较弱，在0.32~0.40，5—9月u, v方向湍流强度较强，在0.40~0.50；垂直方向各月差异和水平方向较为一致，10月至次年4月w方向湍流强度在0.14~0.16，5—9月在0.17~0.20。水平风速的年平均湍流强度为0.4，垂直风速的年平均湍流强度为0.16。

3.3 零平面位移和粗糙度

根据式(3)~式(5)进行迭代拟合，得到10 d平均的标准化零平面位移(d/h)和粗糙度长度(z₀/h)的时间序列(图 3)(其中h为林木平均高度，取值为23 m)。d/h和z₀/h均有明显的季节变化，且一致表现为生长季高，非生长季低的趋势，这与赵晓松等(2004)认为的d和z₀呈相反趋势有很大不同。具体表现为：①2016年12月至2017年4月，为林木非生长季，叶片脱落，d/h和z₀/h较小，平均值分别为0.67和0.049；②2017年5月至9月下旬，为林木生长旺季，叶片密度增加，使得d和z₀迅速增大，平均值分别为0.966和0.056；③2016年10—11月，进入非生长季，d和z₀平均值分别降到0.749和0.051。d和z₀平均值分别为18.56 m和1.21 m，与刘和平等(1997)(d=19.5±1.38 m，z₀=1.60±0.25 m)和赵晓松等(2004)(d/h值为0.642~0.861，z₀/h值为0.0466~0.0736)在长白山森林下垫面(平均树高为26 m)的计算结果较为接近。

3.4 无量纲风速标准差随稳定度变化特征

图 4显示，无论夏季还是冬季，观测区风速标准差在z/L < 0不稳定条件下均满足1/3方相似规律，且经验常数C_u1≈C_v1＞C_w1，在接近中性条件时趋于常数，在稳定和不稳定条件情况时，随着稳定度增大，湍流离散度增大。观测区在不稳定条件时无量纲风速标准差表现为(表略)：经验常数C_u1在春季最大、夏季最小，经验常数C_v1在春季最大、秋季最小，经验常数C_w1在春季最大、夏季最小。从表 1中可以看到高原、草原、沙漠相对于其他下垫面C_u1和C_v1较大，而C_w1相差不大，说明地形作用对垂直方向的湍流运动影响较小。另外，相对于其他地区的森林下垫面，本文的C_u1和C_v1更大一些，C_u1和C_v1大小主要反映局地观测点周围地形的影响(张艳武等，2009)，和前文描述的观测区附近的复杂的地形相符合。

3.5 无量纲温度、湿度标准差随稳定度变化特征

在不稳定条件下(图 5a，5c)，σ_T/|T_*|随z/L变化基本满足-1/3次方相似规律，在强不稳定层结自由对流条件时，两者的拟合关系更为明显。越接近中性条件，点的离散度越大，主要是因为此时热通量接近于零，难以被测量且测量误差较大(张艳武等，2009)。春、秋、冬季拟合参数C_T的差异不大(表略)，分别为2.06、1.94、1.96，夏季拟合参数C_T相对较大，为2.45。在稳定条件下(图 5b，5d)，σ_T/|T_*|随z/L不满足-1/3次方相似规律，点较为离散。强稳定和强不稳定层结时，σ_T/|T_*|的值均较小，随着稳定度向中性条件靠近，σ_T/|T_*|的值逐渐增大。

在不稳定条件下(图 6a，6c)，σ_q/|q_*|随z/L变化也基本满足-1/3次方相似规律，但湿度标准差离散程度明显大于温度标准差。春、夏、秋、冬季拟合参数C_q的差异较明显(表略)，分别为2.67、2.18、2.85、3.00。在稳定条件下(图 6b，6d)，σ_q/|q_*|随z/L也不满足-1/3次方相似规律，点相对温度标准差更为离散，与张艳武等(2009)的结果较为一致。

从表 2不同地区不稳定条件下拟合的系数C_T和C_q可见，不同下垫面C_T和C_q相差较大，其中本站C_T小于美国威斯康星州北部国家森林的森林下垫面(李萍阳等，2002)和河北固城平原下垫面(郭建侠，2006)，大于其他下垫面，C_q没有明显的比较规律。以上特征不仅与下垫面有关，还和不同观测高度以及不同观测季节和气候条件都有很大关系。

3.6 湍流输送系数

总体输送系数包括动量拖曳系数系数(C_D)、热量输送系数(C_H)和水汽输送系数(C_E)，它们是陆面过程中计算不同下垫面地-气间物质和能量交换的重要参数，因此，获得准确的总体输送系数是陆面过程参数化研究的关键(王慧等，2008；张海宏等，2019)。本文采用涡动相关法对C_D进行计算：

$ C_{\mathrm{D}}=\frac{u_*{}^2}{U^2} $

(10)

式中U为通量塔50 m高度处的平均风速。

如图 7所示，与其他下垫面研究中(Wang et al，2016)认为的C_D随着不稳定性的增加而单调增加不同，本文中C_D先从强不稳定到弱不稳定逐渐升高，直到达到弱不稳定条件的峰值(z/L≈-0.47)后，再随稳定度的增加而减小，这与Srivastava and Sharan(2015)的研究相似。经计算，当大气处于不稳定时(z/L<0)，C_D平均值为13.6×10^-3，近中性时(-0.1<z/L < 0.1)，C_D平均值为13.6×10^-3，稳定时(z/L>0)，C_D平均值为7.3×10^-3，整体C_D平均值为9.8×10^-3。本文计算出的C_D值相比其他下垫面较大(王慧等，2008)，可能是研究区地形起伏，森林冠层的粗糙程度在各个季节变化较大，从而对各个方向的湍流风速造成较大影响(张功等，2018)。

图 8给出了C_D与风速之间的关系。在不稳定情况下，当风速较低时，大气湍流主要以热力湍流为主，因此冠层以上大气处于局地自由对流状态，此时会导致C_D较大；而随着风速增加，湍流的剪切作用增强，风切变的增加会产生更多的机械湍流(Wang et al，2015)，使C_D逐渐减小，当风速超过6 m·s^-1时，C_D随风速的变化变得很小。在稳定情况下，在低风速时，C_D随风速而增加，随着风速增大，大气层结逐渐变为中性，C_D变得稳定。z₀是影响C_D的重要因子(王慧等，2008)，图 9中显示C_D随z₀的增加而增加，且近似呈线性变化。

4 结论

利用中国气象局五营森林生态监测站内梯度观测塔上的2016年12月至2017年11月的涡度相关资料，分析小兴安岭森林下垫面的风向、风速特征，主要研究近地层湍流动力特征(湍流强度、湍流方差等)，并对空气动力学参数(d、z₀和C_D)进行了计算。主要结论如下：

(1) 研究区整个观测期间以0~4 m·s^-1的风速最多，夏季以西南风为主，其次为东北风；春、秋、冬季风向以西南风为主，其次为偏西风和偏东风；全年基本为偏西到偏南风，西北方向和东南方向的风很少。u，v，w方向的湍流强度均在5—9月较强，在10月至次年4月较弱。水平风速的平均湍流强度为0.4，垂直风速的平均湍流强度为0.16。

(2) 研究区d和z₀均有明显的季节变化，且一致表现为生长季高，非生长季低的趋势，d和z₀平均值分别为18.56 m和1.21 m，与长白山森林下垫面的结果接近。

(3) 在不稳定条件，各个季节的风速标准差均符合1/3次方相似规律，且C_u1≈C_v1＞C_w1，夏季近中性条件，无量纲三维风速标准差σ_i/u_*(i=u，v，w)为2.62、2.53、1.10；温度和湿度的标准差在不稳定条件符合-1/3次方相似规律，夏季拟合的系数分别为2.45和2.18，在稳定条件不符合相似规律，且湿度相对温度较为离散。在近中性条件时，σ_u/u_*和σ_v/u_*值相对高原、草原、沙漠和极地下垫面较小，比其他下垫面和其他地区的森林下垫面大，而σ_w/u_*相差不大。

(4) 本文中C_D先从强不稳定到弱不稳定逐渐升高，直到达到弱不稳定条件的峰值(z/L≈-0.47) 后，再随着稳定度的增加而减小，当大气处于不稳定时(z/L < 0)，C_D平均值为13.6×10^-3，近中性时(-0.1 < z/L < 0.1)C_D平均值为13.6×10^-3，稳定时(z/L>0)平均值为7.3×10^-3，整体平均值为9.8×10^-3。在不稳定情况下，当风速较低时，C_D较大，而随着风速增加，C_D逐渐减小，当风速超过6 m·s^-1时，C_D随风速的变化变得很小。在稳定情况下，在低风速时，C_D随风速而增加，随着风速增大，大气层结逐渐变为中性，C_D变得稳定。C_D随粗糙度近似呈线性增长。

本文采用梯度观测塔50 m处的涡度观测资料进行分析，由于观测塔只有一层涡度观测系统，且各层梯度数据误差较大，无法验证50 m高度处是否存在常通量层，文中假定50 m位于常通量层，后续会对此进行深入分析。文中采用50 m处的涡度观测资料对d和z₀进行了计算，得出的结果和其他研究结果较为相似，但是在个别月份d数值偏大，可能与计算方法有关，也可能和观测场地附近的地形有关，需要进行更加深入的分析。

致谢：感谢中国气象局沈阳大气环境研究所提供了研究所需的数据，数据来源于“东北地区生态与农业气象野外科研试验基地”。