引言

总体输送系数包括地表动量拖曳系数(C_D)、热量输送系数(C_H)和水汽输送系数(C_E)，当风速 < 10 m·s^-1时，通常认为C_H≈C_E。在陆气相互作用和大气数值模拟研究中，总体输送系数是计算不同下垫面地表与大气之间物质和能量交换的关键参数，在大气环流和气候学研究中，也是计算地表热源强度最重要的参数之一，因此得到准确的总体输送系数是陆面过程参数化研究的关键(张强和卫国安，2004；王慧等，2008)。地表粗糙度亦即空气动力学粗糙度长度，其定义为风速为0 m·s^-1的高度，与下垫面粗糙元的形态学特征和空间分布密切相关，它不仅是描绘下垫面空气动力学特征的重要物理量，而且是研究陆地与大气之间物质和能量交换过程的重要参数之一，准确获得地表粗糙度是改善陆面模式参数化方案、提升模式模拟效果的迫切需要(尚伦宇等，2010)。粗糙度随下垫面性质变化明显，与稳定度呈正相关，与风速呈负相关，摩擦速度随粗糙度增大而减小(何清等，2008)。非中性大气层结条件下，由地表粗糙度不均匀性所致的平均风速、位温梯度以及近地层大气稳定度的次网格分布都对感热通量计算产生影响(陈斌等，2010)。陆面变量(参数)扰动首先改变地表的潜热通量和感热通量，而地表通量的改变会通过陆气相互作用对局地大气的温、压、湿、风产生较大影响(王洋等，2014)。

目前大气科学界对于总体输送系数和地表粗糙度的研究已经取得了一些成果。李国平等(2002b)利用西藏的4个自动气象站的近地层梯度资料，用最小二乘法确定了各站各季节的地表粗糙度，应用廓线-通量法计算了总体输送系数并分析了其随时间的变化特征，发现青藏高原动量输送系数的多年平均值为3.53×10^-3~4.99×10^-3，热量输送系数的多年平均值为4.67×10^-3~6.73×10^-3，还讨论了总体输送系数与近地层大气层结稳定度、地表粗糙度以及地面风速等因子的关系，初步建立了可用常规气象站地面观测资料计算青藏高原总体输送系数的拟合公式。杨兴国等(2010)利用在陇中黄土高原观测资料，采用空气动力学法计算了动量和感热总体输送系数，发现陇中黄土高原半干旱区动量和感热总体输送系数受下垫面植被的影响，在一年中呈现出双峰型特征，当大气处于不稳定状态时，总体输送系数随着风速的增大而减小，当大气处于稳定状态时，随着风速的增大而增大。岳平等(2015)利用SACOL站夏季晴天近地层湍流观测资料确定了大气动力学和热力学粗糙度长度，发现总体输送系数随稳定度的增大而减小。孙俊等(2012)利用廓线法计算了黑河地区的总体输送系数和地表粗糙度，发现地表粗糙度与植被覆盖度和高度以及下垫面的性质有关，下垫面状况影响动量总体输送系数对稳定度的依赖程度。李锁锁等(2010)利用黄河源区湍流观测资料结合单层超声观测资料计算了黄河上游玛曲地区草原下垫面空气动力学粗糙度和零平面位移并应用于陆面过程模式CoLM中，改进陆面参数后的模式对感热通量和潜热通量的模拟均有明显改善。陈世强和吕世华(2013)计算了金塔试验区内戈壁和沙漠的动力学和热力学粗糙度长度，代入Noah陆面模式，模拟的戈壁、沙漠上的地表温度和感热通量同观测值较为一致，提高了该模式在沙漠、戈壁特殊区域的模拟能力，有利于将耦合了Noah模式的中尺度模式更好地应用到绿洲系统的研究中。张果等(2016)针对Noah和Noah MP两套陆面物理过程参数化方案进行了评估，认为Noah MP方案提高了土壤水分和土壤温度，在东亚区域的整体模拟效果。杨耀先等(2014)利用那曲高寒气候与环境观测站的资料，应用一种独立的确定地表动力学粗糙度的方法及两种热力学粗糙度的参数化方案，得出了动力学粗糙度、热力学粗糙度以及附加阻尼的变化规律，发现动力学粗糙度在一定时间尺度上存在着波动，热力学粗糙度在高原季风前、盛行期、衰退期有不同的日变化和季节变化特征。

青藏高原地区冻土分布广泛，冻土中冰的存在极大地改变了土壤的热力性质。土壤冻结和融化时会释放或吸收大量的热量，从而影响能量在土壤层中的分配和地表能量平衡，冻土过程的模拟对于地气相互作用、区域气候模拟和全球气候变化极为重要(李震坤等，2011)。积雪覆盖地表会阻碍地气之间的能量交换，积雪通过表面不同的反照率和不同的湍流通量形成了陆面与大气间独特的能量交换，影响近地层气象要素特征，反过来其对湍流和能量交换又有重要影响(李丹华等，2017)。由于青藏高原腹地人迹罕至，交通不便，观测资料匮乏，高原地区土壤冻结、积雪覆盖对陆气相互作用的影响研究较少。本文利用青海省气象科学研究所玉树隆宝野外观测站的微气象及涡动相关系统观测数据，通过分析探讨了未冻结、冻结和冻结有积雪覆盖三种情况下动量通量和感热通量的日变化情况，计算了三种情况下动量总体输送系数、感热总体输送系数、动力学粗糙度和热力学粗糙度，分析了附加阻尼和粗糙度雷诺数的关系，并将三种附加阻尼的参数化方案进行了比较，为全面认识青藏高原地区陆气相互作用特征提供科学支持。

1 计算方法与观测资料

直接求解总体输送系数的方法主要有三种，分别是涡动相关法、廓线-通量法和经验函数法，其中涡动相关法最为准确(高世仰等，2017)。动量总体输送系数C_D和感热总体输送系数C_H可分别通过式(1)和式(2)计算而得(岳平等，2013)：

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

(1)

$ C_{\mathrm{H}}=\frac{\overline{w^{\prime} T^{\prime} }}{u\left(T_{\mathrm{g}}-T_{\mathrm{a}}\right)} $

(2)

式中，u_*为摩擦风速(单位：m·s^-1)，u为参考高度的水平风速(单位：m·s^-1)，w′和T′分别为垂直风速脉动量(单位：m·s^-1)和温度脉动量(单位：℃)，T_g-T_a为地气温差(单位：℃)。

动力学粗糙度是近地面本身的一种特性，是指地面上方风速为0 m·s^-1的高度，当流体流经地表面时，不同地表粗糙度对流体的影响程度也不同，常常用来度量地面对气流的粗糙程度(陈金雷等，2017)。长期以来国内外学者对于动力学粗糙度从计算方法和下垫面类型等方面做了一系列工作，主要方法有：风廓线拟合法、牛顿迭代法、TVM (temperature variance method)法、Martano法、无因次化风速法等。本文采用Yang et al(2008)对数风廓线方法，根据Monin-Obukhov相似理论，含有层结稳定度订正函数的近地层风速廓线方程为：

$ \ln \left(\frac{z}{z_{0 \mathrm{m}}}\right)=\frac{k u}{u_{*}}+\psi_{\mathrm{m}}(\zeta) $

(3)

$ {\psi _m}(\zeta) = \left\{ {\begin{array}{*{20}{l}} {2\ln \left({\frac{{1 + x}}{2}} \right) + \ln \left({\frac{{1 + {x^2}}}{2}} \right) - }\\ {2{{\tan }^{ - 1}}x + \frac{\pi }{2}, \;\;\;\;\zeta < 0}&{}\\ { - 5\zeta, \;\;\;\;\;\;\zeta > 0} \end{array}} \right. $

(4)

$ x=(1-16 \zeta)^{1 / 4} $

(5)

式中，z_0m为动力学粗糙度，k为冯卡曼常数，取值0.4，ζ=z/L为大气稳定度参数，ζ>0.01时为稳定层结，ζ < -0.01时为不稳定层结，L为Monin-Obukhov长度。由式(3)~式(5)可以推出动力学粗糙度对数：

$ \ln z_{0 m}=\ln z-\frac{k}{\sqrt{C_{D}}}-\psi_{m}(\zeta) $

(6)

进而可计算出动力学粗糙度。

热力学粗糙度是指大气近地层满足Monin-Obukhov相似理论时温度廓线外延到空气温度等于地表温度时的高度(陈金雷等，2017)。根据Monin-Obukhov相似理论的方程：

$ C_{\mathrm{D}}=\frac{k^{2}}{\ln ^{2}\left[\frac{z}{z_{0 \mathrm{m}}}-\psi_{\mathrm{m}}(\zeta)\right]} $

(7)

$C_{\mathrm{H}}=\frac{k^{2}}{\ln \left[\frac{z}{z_{0 \mathrm{m}}}-\psi_{\mathrm{m}}(\zeta)\right] \cdot \ln \left[\frac{z}{z_{0 \mathrm{h}}}-\psi_{\mathrm{h}}(\zeta)\right]} $

(8)

$\psi_{h}(\zeta)=\left\{\begin{array}{l}{2 \ln \left(\frac{1+y}{2}\right), \quad \zeta<0} \\ {-7.8 \zeta, \quad \zeta>0}\end{array}\right. $

(9)

$y=0.95(1-11.6 \zeta)^{1 / 2} $

(10)

由式(7)~式(10)可以得到热力学粗糙度对数，

$\ln z_{0 \mathrm{h}}=\ln z-\frac{k \sqrt{C_{\mathrm{D}}}}{C_{\mathrm{H}}}-\psi_{\mathrm{h}}(\zeta) $

(11)

进而可计算出热力学粗糙度。

青海省气象科学研究所玉树隆宝野外观测站位于玉树州隆宝镇(33°10′N、96°34′E)，海拔高度为4212 m，下垫面为沼泽性草甸覆盖的高寒湿地区(张海宏等，2017)。图 1给出了玉树隆宝观测站的位置及照片，建站时间为2011年10月，2014年增添了涡动相关和雪深等观测系统。观测的物理量包括空气温度、空气湿度、风速、风向、大气压、短波辐射、长波辐射、三维超声风、超声虚温、土壤温度、土壤湿度、土壤热通量、积雪深度、水汽和二氧化碳通量、甲烷浓度等。观测数据由数据采集器CR5000处理并存储，所有仪器由3块35 W的太阳能板和2个120 Ah的电瓶供电，除仪器拆装和天气原因造成供电短暂中断外，一直连续进行观测。本文所选资料时间段为2014年11月至2015年1月，其中未冻结状态所用资料为2014年11月8—15日，期间5~40 cm土壤温度>0℃，冻结状态所用资料为2014年11月22—29日，期间5~40 cm土壤温度 < 0℃，冻结有积雪覆盖状态所用资料为2015年1月6—13日，期间5~40 cm土壤温度 < 0℃且SR-50A雪深仪监测结果显示积雪深度在4~6 cm。

2 结果与分析 2.1 动量通量和感热通量的日变化

陆面过程研究的核心问题是下垫面与大气之间的能量传输和物质交换。陆面与大气之间的动量、能量和物质交换通过陆气通量反映出来，在大气动力学方程中可被描述为一些与下垫面有关的源、汇项，如决定风速变化的摩擦力项，决定大气温湿变化的感热项和潜热项等(张瑛等，2011)。土壤的热容量远大于空气，土壤的热状况及其变化将会对大气的陆面下边界起重要的作用，在环境相同的条件下，雪地表面的感热通量比裸地表面的感热通量小很多(单机坤等，2013)。图 2为玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖情况下动量通量(τ)和感热通量(H)的日变化情况，所用观测资料为典型晴天条件下的30分钟数据的平均值。从动量通量的日变化情况来看，在未冻结状态下动量通量的日变化幅度为0.27 kg·m^-1·s^-2，冻结状态下的日变化幅度达到0.48 kg·m^-1·s^-2，冻结有积雪覆盖时的日变化幅度最小，仅有0.05 kg·m^-1·s^-2。从感热通量的日变化情况来看，未冻结状态下的日变化幅度为90 W·m^-2，冻结状态下的日变化幅度达到180 W·m^-2，冻结且有积雪覆盖时的日变化幅度为80 W·m^-2。冻结状态下动量通量和感热通量的日变化幅度最大，冻结有积雪覆盖时动量通量和感热通量的日变化幅度最小，未冻结、冻结和冻结有积雪覆盖三种状态下，动量通量和感热通量白天的值差异较为明显，而夜间差异较小。

2.2 总体输送系数

动量总体输送系数和热量总体输送系数分别表征了湍流动力作用和湍流热力作用，是衡量湍流强弱程度的物理量(周明煜等，2000)。通常利用总体输送系数的参数化公式确定陆气之间能量和物质的交换。总体输送系数不仅与大气动力状态存在联系, 而且与大气热力状态密切相关(岳平等，2013)。总体输送系数对层结稳定度的变化较为敏感(高世仰等，2017)。当地表有植被覆盖时，会导致动量总体输送系数增大，而感热总体输送系数减小(王澄海等，2007)。图 3为玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖状态时不稳定层结和稳定层结条件下摩擦速度平方(u_*²)与水平风速平方(u²)的关系，回归直线的斜率即代表动量输送系数C_D。图 4为玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖状态时不稳定层结和稳定层结条件下w′T′与u(T_g-T_a)的关系，回归直线的斜率即代表感热输送系数C_H。

对于青藏高原地区总体输送系数的计算，前人已做过很多研究。李国平等(2002a；2002b；2003)利用通量廓线法计算得到的青藏高原那曲、改则、狮泉河地区C_D为4.3×10^-3~4.8×10^-3，C_H为5.7×10^-3~6.6×10^-3，周明煜等(2000)利用涡动相关法计算出当雄、昌都地区C_D为1.8×10^-3~4.4×10^-3，C_H为1.5×10^-3~4.7×10^-3，高志球等(2000)利用涡动相关法计算出那曲地区C_D为3.7×10^-3，钱泽雨等(2005)利用空气动力学方法计算得到的北麓河地区的C_D为1.74×10^-3，C_H为1.37×10^-3~5.93×10^-3，前人的研究多集中于总体输送系数在不同季节的变化，对于冻结前后和积雪覆盖状态下总体输送系数的差异性研究并不多见。表 1给出了玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖状态下C_D和C_H的值。不稳定层结和稳定层结条件下，C_D和C_H的值均表现为冻结状态下最高，冻结有积雪覆盖状态下最低，这反映了高寒湿地下垫面的湍流作用在土壤冻结之后增强，有积雪覆盖时减弱。

2.3 动力学粗糙度和热力学粗糙度

图 5为玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖时动力学粗糙度的对数lnz_0m和热力学粗糙度的对数lnz_0h的频率分布及五点平滑曲线，根据滑动平均曲线最高点确定的动力学粗糙度z_0m和热力学粗糙度z_0h的值见表 2，z_0m和z_0h在冻结状况下最小，冻结有积雪覆盖时最大。z_0m在冻结和积雪覆盖前后变化差异较z_0h明显。尚伦宇等(2010)通过计算玛曲地区冻融过程的地表粗糙度并分析变化情况发现，在从未冻结、冻结至融化后的过程中，地表粗糙度呈逐渐减小的趋势。王玉玉等(2014)研究结果表明，当地表有积雪覆盖时，在风的吹拂下雪粒发生跃移会造成地表粗糙度增大，当地面积雪融化或者风速不足以吹拂雪粒运动时，粗糙度相对较小。陈金雷等(2017)计算了黄河源区曲麻莱地区夏季的z_0m和z_0h，发现z_0h小于z_0m。本文得出的z_0m和z_0h的结果与前人研究结论较为一致。

2.4 附加阻尼的参数化

附加阻尼(kB^-1)是研究地表与大气之间物质和能量交换过程的基本参数，也是陆面过程模式与地表通量遥感估算模型的重要变量之一，附加阻尼定义式为kB^-1=ln(z_0m/z_0h)，影响kB^-1的因子较多，如气象条件、植被结构以及下垫面状况等(鞠英芹等，2014)。Brutsaert(1975；1982；1998)认为kB^-1依赖于表面性质，与粗糙度雷诺数有关。Trou-fleau et al(1997)研究表明，kB^-1和许多因素有关系，包括结构参数和气象条件等。Lhomme et al(1997)利用Shuttlewoth-Wallace两层模型，结合对地表辐射温度的线性假设，获得kB^-1的解析解后得出，kB^-1对于某一类冠层不是一个常数，随叶面积指数、株高、植被覆盖率、水分胁迫以及气象条件的变化而变化。由于kB^-1不能从观测中直接获得，需要多个要素的观测通过一系列计算得到，这些输入变量在观测中的任何误差都会对kB^-1产生影响，因此kB^-1是一个难以确定的量。过去几十年，不同研究者已发展了多个kB^-1的参数化方案。周德刚(2016)选取了一些常用的热力参数化方案，通过敦煌站夏季估算的感热通量与野外观测的比较，评价了这些参数化方案在西北干旱区的适用性，结果表明Y08方案(Yang et al，2008)估算的感热通量相对比较合理，可以用来研究西北干旱区的夏季地表感热输送特征。

在青藏高原地区，由于下垫面类型较为复杂，何种参数化方案对于高寒湿地下垫面较为合适尚无研究结论。本文选取最新的三种kB^-1参数化方案见表 3。粗糙度雷诺数Re_*=u_*z_0m/ν，ν为运动学黏性系数，本文取值为1.48×10^-5 m²·s^-1。图 6给出了玉树隆宝湿地未冻结、冻结和冻结有积雪覆盖情况下kB^-1与Re_*的关系，经优选拟合得出二者的关系为幂函数型。

将三种参数化方案计算得到的kB^-1值与利用观测资料计算得出的kB^-1值进行误差分析(表 4)，未冻结、冻结和冻结有积雪覆盖状态下参数化方案1计算得到的kB^-1值与观测值最为接近，这表明青藏高原湿地下垫面在未冻结、冻结和冻结有积雪覆盖状态下kB^-1与Re_*之间满足幂函数型关系，第一种参数化方案较为合适。

3 结论

本文利用玉树隆宝湿地观测站2014年12月至2015年1月微气象及涡动相关系统的观测资料，分析了未冻结、冻结和冻结有积雪覆盖三种情况下动量通量和感热通量的日变化情况，计算了三种情况下动量总体输送系数、感热总体输送系数、动力学粗糙度和热力学粗糙度，分析了附加阻尼和粗糙度雷诺数的关系，并将三种附加阻尼的参数化方案进行了比较，主要结论有：

(1) 冻结状态下动量通量和感热通量的日变化幅度最大，冻结有积雪覆盖时，动量通量和感热通量的日变化幅度较小。

(2) 未冻结、冻结和冻结有积雪覆盖三种情况下动量总体输送系数C_D和热量总体输送系数C_H的值在冻结时最大，冻结有积雪覆盖时最小。

(3) 未冻结、冻结和冻结有积雪覆盖三种情况下动力学粗糙度z_0m和热力学粗糙度z_0h在冻结状况下最小，冻结有积雪覆盖时最大。

(4) 本文选取的三种附加阻尼kB^-1参数化方案中，幂函数型方案对高寒湿地下垫面较为合适，得出的kB^-1值与观测值最为接近。