青藏高原的热力和动力作用对大气环流和天气气候的影响一直是气象科学的关键问题, 也是世界气象界关注的重点和难点(段安民等, 2003; 马伟强等, 2004; 马耀明等, 2006; 阳坤等, 2010)。地面感热通量在高原大气热量源汇和地面热量平衡中起着至关重要的作用(李栋梁等, 2003), 在地表感热通量的估算中, 一个关键参数是地表热力输送系数(C_H), 它在陆气相互作用和大气数值模拟研究中也起着重要作用。

近年来, 对于各种主要类型下垫面热力输送系数的估算及其区域分布和时空变化特征的分析等方面都取得了一定的研究成果(王慧等, 2008)。研究表明, 植被属性如植株类型、覆盖度等会影响粗糙度, 从而影响热力输送系数的大小。例如, Kondo and Kawanaka(1986)研究了植被类型和热力输送系数的关系, 结果表明热力输送系数随着植被密度而变化。李国平等(2002)得到拉萨站的热力输送系数季节差异较大, 并分析这与地表植被的变化状况有关。除了植被属性外, 微尺度地形和土壤类型也会对热力输送系数产生影响。李国平等(2002)计算了拉萨(河谷), 日喀则(盆地), 那曲(高原)等不同下垫面的热力输送系数, 结果差异较大。直接求解热力输送系数的方法主要有三种, 分别是湍流脉动相关法(涡动相关法)、廓线-通量法(总体输送法)和常规资料参数化法(经验函数法)。其中湍流脉动相关法最为准确, 但对观测仪器的精度和性能要求较高, 只能适用于个别有条件的台站。例如, 岳平等(2013)根据空气动力学法和涡动相关法计算了锡林郭勒草原的热力输送系数, 并建立了热力输送系数与大气稳定度和近地层气温的关系。冯健武等(2012)采用涡动相关法计算了通榆地区退化草地和农田下垫面的地表空气动力学粗糙度(z_0m)、热力学粗糙度(z_0h)、地表拖曳系数(C_D)和热力输送系数(C_H), 并发现热传输附加阻尼(kB^-1)和总体输送系数都存在明显的季节和年际变化。常规资料参数化法是利用风速、气压或地形高度等参数的拟合函数, 精度不高, 而且不同拟合函数的结果差别较大。例如, 王慧和李栋梁(2010)拟合得到了针对不同下垫面类型的区域地表C_H与归一化植被指数(NDVI)之间的经验关系式。目前应用比较广泛的是廓线-通量法, 例如, 李国平(2000)利用廓线-通量方法计算了改则和狮泉河地区的热力输送系数, 并根据总体方程估算了地表动量通量、感热和潜热。

由于青藏高原海拔高、范围广、地形复杂, 高原上不同区域的地形及地表状况的不同造成气候条件的较大差异, 加之缺少对高原的详细、持续、大范围的观测, 并且不同作者使用的研究方法各不相同, 对于高原热力输送系数的估算结果之间差别较大, 有时偏离实际值较多。而且目前多使用台站资料进行研究, 局地观测确定出的总体输送系数无法有效推广至区域尺度和非均匀下垫面(李国平等, 2002; 岳平等, 2013; 冯健武等, 2012)。因此, 建立一种适用于整个青藏高原地区非均匀下垫面的热力输送系数的估算方法, 对准确估算青藏高原地表感热通量以及改进全球气候模式和数值预报模式中高原区域的参数化方案具有重要的实际意义。由于卫星遥感资料具有探测范围大、获取资料快、受地面条件限制小、获取信息量大的优点, 本文结合多种卫星遥感资料, 对廓线-通量法加以改进, 通过建立更准确的动力粗糙度和热传输附加阻尼的估算方案, 以求得适用青藏高原非均匀下垫面的热力粗糙度, 从而进一步估算热力输送系数。

在2005-2009年期间, 中日科学家联合开展了JICA(Japan international cooperation agency)“中日气象灾害合作研究中心”试验, 建立了青藏高原及其东部大气综合监测网。本文选取JICA项目数据集中2008年大理(DL)、林芝(LZ)、温江(WJ)三个站的边界层通量塔观测系统数据, 包括边界层观测系统向上、下长、短波辐射(长、短波辐射计: CG4、CM21, Kipp & Zonen), 超声风、温(10次·s^-1)(三维超声仪: CSAT3, Campbell), 铁塔梯度风(MetOne034B)、温、湿(HMP45) 度数据、土壤热通量(HFP01SC-L50, Hukseflux)及土壤温度(STP01-L50, Hukseflux)、湿度(CS616-L, Campbell)等综合观测数据。大理边界层通量观测系统位于(25.42°N, 100.11°E), 海拔1900.5 m, 观测场主要为农田生态。林芝边界层通量观测系统位于(29.46°N, 94.44°E), 海拔3327 m, 观测场下垫面较为平坦, 四周为林地, 地表为高原草甸。温江边界层通量观测站位于(30.42°N, 103.5°E), 海拔约539 m, 下垫面地势平坦, 以农田生态为主(蔡雯悦等, 2012)。本文利用这三个试验站的涡动相关观测资料对kB^-1的计算方案进行对比与改进。

此外, 选取中日合作项目“全球能量水分循环亚洲季风青藏高原试验”(GAME-Tibet, 1996-2000)中2008年的BJ(那曲布交)、MS3478/NPAM、D105、Amdo(安多)4个站的梯度观测数据进行验证。BJ站位于那曲河谷高寒草原, 观测点经纬度为31.37°N, 91.90°E, 海拔为4509 m。观测站场地平坦, 四周开阔, 地表主要为砂土, 有稀疏的细石块分布, 不均匀的生长着高度10~20 cm的高寒草甸。D105站位于(32.06°N, 91.94°E), 海拔约5039 m, 试验场地平坦, 视野开阔, 东边约20 km处有山丘, 地面为高原稀疏草甸。MS3478/NPAM站位于(31.56°N, 91.43°E), 海拔约4620 m, 试验场地面平坦、四周开阔, 东边约5 km、西边约30 km、南北约10 km处分布着相对高度约100~200 m的山丘, 地面被高约15 cm高原草甸覆盖。安多(Admo)站位于唐古拉山南坡, 观测点经纬度为(32.14°N, 91.37°E), 海拔高度为4696 m。地处安多草原, 地表均匀覆盖约2~5 cm的高原草甸。试验场地的地形和下垫面特征基本代表了藏北高原大面积特征, 即高原典型的高寒草甸下垫面。来自四个试验站的数据集提供了梯度数据和其他物理量, 如风速、风向、湿度、相关高度上的气温、地面气压和土壤温度等, 数据间隔为30 min。这四个站的梯度观测数据作为非均匀地形地表的代表, 用来检验本文的估算方案。另外, BJ站同时具有大气边界层塔, 提供了涡动相关观测数据集, 包括摩擦速度(u_*)、温度尺度(T_*)、感热通量及其他通量、超声波风速、温度和湿度等等。本文采用该站的涡动相关观测资料来对比和检验梯度观测数据得到的结果。

由于青藏高原地区气象台站的稀缺, 站点资料难以覆盖整个高原, 因此, 采用遥感资料来估算整个高原的热力输送系数, 以分析其分布特征。MODIS植被指数合成产品中, 最为常用的有两种, 分别是传统NDVI(归一化植被指数)和通过对大气与土壤校准因素自我纠正改善后的EVI(增强型植被指数)(Leeuwen et al, 1999)。EVI是对NDVI的发展和延续, 从植被指数计算公式和合成方法两方面做了改进, 避免了NDVI在植被高覆盖区易饱和的问题, 考虑了土壤背景对植被指数的影响, 对气溶胶等残留做了进一步校正。EVI时间序列相较于NDVI时间序列季节性更明显, 能够更好地反映高植被覆盖区的季节性变化特征, 并且很少有突降现象, 时间序列曲线较平滑(李文梅等, 2010)。因此本文选取Aqua/MODIS 1km分辨率的MYD13A2 EVI 16天最大值合成数据集, 作为植被覆盖度数据; 从MODIS产品MCD12Q1获得的国际地圈-生物圈计划(IGBP)约1 km分辨率的8天合成植被类型数据, 用于估计植被高度; 结合MOD15A2的1 km分辨率的叶面积指数(LAI)8天合成产品计算粗糙度; 陆面温度采用分裂窗算法的MYD11A1的8天1 km分辨率的数据; 发射率和反照率产品分别采用MYD11A2和MCD43B3的8天1 km分辨率的数据集; 欧洲中期天气预报中心(ECMWF)ERA-Interim的08:00(世界时, 下同)地面温度和2 m气温数据用来对比和检验06:00估算的地-气温差, 10 m风速数据用来估算相关子参数和粗糙度, 分辨率为0.75 km; 地形资料采用GLCF/SRTM3(Shuttle Radar Topography Mission)即“航天飞机雷达地形测绘使命数据”数字高程产品, 选取30 m分辨率DEM资料, 提取地形特征。本文结合上述多种资料, 全部升尺度到1 km分辨率的尺度上, 建立基于地形特征因子和EVI的粗糙度估算模型; 结合LAI和植被类型信息, 分析其对表层土壤温度及地-气温差的影响, 用LAI修订地表土壤温度估算方法和地-气温差(Zhang et al, 2015), 最后研究非均匀地形和植被背景下的空气动力学粗糙度、热力输送系数。

理论上, 空气动力学粗糙度z_0m是平均风速随高度减小到零时的高度。z_0m是陆气相互作用中一个重要的地表参数, 它不仅能反映下垫面的状况, 而且表征陆气间物质、能量以及水分交换过程的强弱、大小(杨耀先等, 2014)。研究表明, 地表粗糙度的大小并不与地形高度完全成正比, 那些用地形高度推算地表粗糙度的经验函数关系只适用于估计较大区域的平均地表粗糙度(李国平等, 2002)。空气动力学粗糙度对植被区域植株的密度、高度、叶面积指数、覆盖度、风速等的响应都很敏感。粗糙单元分布密度对流场影响较大, 随粗糙单元分布密度的增大, 空气动力学粗糙度将逐渐增大(张杰等, 2010)。由于季节而导致地物状况的变换, 空气粗糙度就具有季节性的变换(张强等, 2012)。因此如何连续地确定出能反映非均匀下垫面逐年各季节甚至各月变化特征的地表粗糙度, 关系到进一步提高热力输送系数计算的客观性和准确性。计算空气动力学粗糙度的方法有多种, 最直接的方法是拟合风廓线法。

地气间的能量传输的研究通常建立在莫宁一奥布霍夫相似(MOS)理论基础上, 近地层相似理论的成就在于理论上阐明了无量纲风、温梯度是M-O参数ζ的函数, 非中性大气层结下的风温廓线公式通常写成(赵鸣, 2006):

式中: k(≈0.4) 为Karman常数, u为平均风速, u_*为地表摩擦风速, θ为高度z处位温, θ_s为地表位温, T_*为温度尺度, d为零平面位移高度, z_0m为空气动力学粗糙度, z_0h为热力学粗糙度, ${\zeta = \frac{{z - d}}{L}}$, 其中L为奥布霍夫长度。φ_m(ζ)、φ_h(ζ)分别为动量、热量稳定度修正函数, 在中性条件下为零。Su et al(2001)根据Massman(1997)的方法推导得到u_*与u的关系式:

式中: C₁, C₂, C₃分别为0.32, 0.264和15.1, C_d=0.2为植被的拖曳系数。中性大气层结时容易得到空气动力学粗糙度如下:

根据Massman(1997)的植被动量传输模型, 得到${\frac{d}{z} = 1 - \frac{1}{{2n\left( {1 - {e^{ - 2n}}} \right)}}}$, 其中, n为冠层内风速剖面消光系数, 定义$n = \frac{{{C_{\rm{d}}}{\rm{LAI}}}}{{2u_*^2/{u^2}}}$。

本文使用增强型植被指数(EVI)作为植被生长状态及植被覆盖度的指示因子, 表征植被的季节变化, 给定一个比率f_a=0.7, 给出不受大尺度地形影响的水平均一下垫面动力学粗糙度如下:

式中:右端前项为裸地微尺度地形作用项, 后项为植被作用项, a=0.104, b=0.31, h_v为由IGBP的植被类型数据估算的植被高度(表 1)。

表 1(Table 1)

表 1 IGBP植被覆盖类型与植被高度对照 Table 1 The type of vegetation cover from IGBP and corresponding height 编号 植被类型 植被高度/m 1 常绿针叶林 5 2 常绿阔叶林 5 3 落叶针叶林 5 4 落叶阔叶林 5 5 混交林 5 6 封闭灌木 0.3 7 开阔灌木 0.3 8 多树草原 0.2 9 稀树草原 0.2 10 草原 0.2 11 永久性湿地 0.4 12 农田 0.9 13 城市建筑 10 14 庄稼/自然栽培植被 0.9 15 水或冰或雪 0.001 16 贫瘠或稀疏植被 0.15

表 1 IGBP植被覆盖类型与植被高度对照 Table 1 The type of vegetation cover from IGBP and corresponding height

热传输附加阻尼kB^-1是对地表和大气垂直热量交换研究中的重要参数之一, 也是陆面过程模式与地表通量遥感估算模型的重要变量之一。一般可表示为${\rm{k}}{{\rm{B}}^{ - 1}} = \ln \frac{{{z_{0{\rm{m}}}}}}{{{z_{0{\rm{h}}}}}}$, 因此通过kB^-1和z_0m就可以得到z_0h。kB^-1在不同的下垫面的值有着很大的差异, 在一个均匀的植被地表下, 其是一个定值, 而在其它下垫面, kB^-1都是一个难以确定的变化量, 需要建立一个适合的kB^-1参数化模型。

多年来, 有许多学者提出过各种关于kB^-1的计算方案, 并广泛应用于陆面模式中(表 2)。

表 1(Table 1)

表 2 不同的kB-1参数化方案 Table 2 Different parameterization schemes of kB-1 参数化方案 参考文献 缩写 备注 kB-1=ln(Pr×Re*) Sheppard, 1958 S58 Pr=0.7, Re*=u*z0m/ν kB-1=kα(8 Re*) 0.45Pr0.5 Owenand Thomson, 1963 O63 k=0.4, α=0.52, Pr=0.7 kB-1=2.46Re*0.25-2 Brutsaert, 1982 B82 Re*=u*z0m/ν kB-1=Au(Ts-Ta) Kustas et al, 1989 K89 A=0.17 kB-1=kαRe*0.45 Zengand Dickinson, 1998 Z98 k=0.4, α=0.52 ${\rm{k}}{{\rm{B}}^{ - 1}} = \frac{{k{C_d}}}{{4{C_t}\frac{{{u_*}}}{u}\left( {1 - {e^{ - n/2}}} \right)}}f_c^2 + \frac{{k\frac{{{u_*}}}{u}\frac{{{z_{0m}}}}{h}}}{{C_t^*}}f_c^2f_s^2 + {\rm{kB}}_s^{ - 1}f_s^2$ Su et al, 2001 S01 $n = \frac{{{C_d}{\rm{LAI}}}}{{2u_*^2/{u^2}}},{C_d} = 0.2,{C_t} = 0.01,{f_c} = \frac{{{\rm{EVI}} - {\rm{EV}}{{\rm{I}}_{\min }}}}{{{\rm{EV}}{{\rm{I}}_{\max }} - {\rm{EV}}{{\rm{I}}_{\min }}}},{f_s} = 1 - {f_c}$ kB-1=1.29Re*0.25-2 Kanda, 2007 K07 Re*=u*z0m/ν z0h=(70ν/u*)×exp (-βu*0.5|T*|0.25) Yang et al, 2007 Y07 β=7.2

表 2 不同的kB-1参数化方案 Table 2 Different parameterization schemes of kB-1

早期的参数化方案一般是根据实验结果与粗糙雷诺数Re_*建立简单的关系, 如S58、O63、B82、Z98、K07, 其中最常被研究者使用的是B82方案。后来有学者将kB^-1与其他常规气象参数建立关系, 如风速、温度、摩擦速度u_*等, 如K89、Y07, 因Y07在稀疏植被区能够取得较好的结果, 近几年常被使用。

其中S01方案是Su et al(2001)根据Massman(1999)方法建立的一个物理模型, 公式如表 2中所示, 其中f_c为植被覆盖度, 通常利用植被指数近似估算植被覆盖度, 本文使用EVI得到f_c。C_t为叶片热传输系数, 取0.006n。C_t^*为土壤热传输系数, 这是一个关于普朗特数和雷诺数的函数。S01公式中第一项采用的是Choudhury and Monteith(1988)的植被冠层模型; 第二项是Massman(1999)方法中用以表征植被和土壤混合作用的项; 第三项采用的是Brutsaert(1982)的裸地模型, 即kB_s^-1=2.46Re_*^0.25-2。本文将第三项kB_s^-1采用Yang(2008)的方法, 即表 2中的Y07方案, 根据:

得到

式中: ν是空气运动粘性系数, 取1.328×10^-5$\left( {\frac{{{P_0}}}{P}} \right){\left( {\frac{T}{{{T_0}}}} \right)^{1.754}}$m·s^-1, P₀=1.013×10⁵ Pa, T₀=273.15 K, P(Pa)为地表气压, T(K)温度; β=7.2 s^1/2·m^-1/2· K^-1/4。

另外, 令$B_H^{ - 1} = \frac{1}{k}\ln \frac{{{z_0}}}{{{z_{\rm{h}}}}}$, 由(2) 式可推导出下式(赵鸣, 2006):

令地-气温差dT_sa=θ_s-θ_a, 从(8) 式第一项中可以看出dT_sa/T_*对于kB^-1的作用, 即热力因子的贡献, 第二项为动力因子的贡献。此外, 本文利用JICA试验三个站点(DL、LZ、WJ)的观测资料探究了dT_sa/T_*与dT_sa的关系。图 1是dT_sa/T_*与dT_sa的拟合曲线, 发现二者存在非线性关系, 由于温度容易由卫星资料得到, 为表征热力学因素对kB^-1的作用, 将(8) 式修正为:

图 1(Figure 1)

图 1 地-气温差和温度尺度的比与地-气温差的散点关系 Figure 1 Scatter relationships between temperature scale for heat and temperature difference between surface and air temperature

该式即为本文最终使用的kB^-1计算方案, 修正项$\frac{1}{{{{(d{T_{{\rm{sa}}}})}^{0.{\rm{ }}1}}}}$是通过敏感性试验模拟的结果。

图 2分别给出了使用涡动相关法感热观测值倒推(图 2a)、本文方案(图 2b)、Y07方案(图 2c)、S01(图 2d)方案计算的林芝站的kB^-1。涡动相关法与本文方案得到的值均约为2~4, Y07方案得到的值约为4~7, B82方案的值约为13~17。可以看出本文方案在时间变化趋势和数值大小方面都与实际的kB^-1值有比较好的对应关系, 相对于其他方案而言, 本方案的估算效果最佳。本文的结果更接近实际值的原因主要在于使用地-气温差进行了修正, 结果也证明这种改进是合理的。

图 2(Figure 2)

图 2 涡动相关法(a)与修正方案(b)、Y07 (c)、S01 (d)参数化方案计算的林芝站2008年的热传输附加阻尼kB-1 Figure 2 The kB-1 of LZ station in 2008 calculated by eddy correlation techniques (a) and several parameterization schemes of correction scheme (b), Y07 (c) and S01 (d)

给定大气状态参数、地表动力粗糙度(z_0m)和热力粗糙度(z_0h)等, 热力输送系数(C_H)可用廓线-通量法表示为:

式中: P_r是普朗特数, ${ {\frac{z}{L}} }$≥0时, P_r=1; ${ {\frac{z}{L}} }$ < 0时, P_r=0.95。z_m和z_h分别为风速和气温的测量高度; φ_m和φ_h是对风速和温度廓线的稳定度校正项。式中z_0m由式(5) 得到, z_0h由z_0m和由式(9) 计算的kB^-1而得。

此外, 湍流数据处理采用涡动相关法可以得到式(11)(岳平等, 2013), 本文用该方法来进行合理性检验。

式中: θ_a、θ_s分别是空气和地表位温, ω′是垂直脉动风速, θ′是脉动温度。

利用本文的kB^-1方案估算了2008年四个试验站(BJ、MS3478、D105、Amdo)的z_0m、kB^-1和z_0h(图 3~5), 并观察其季节演变规律。由2008年典型站点的动力粗糙度z_0m日平均估算值及趋势线(图 3)可以看出, 动力粗糙度在夏季较大, 冬季较小。z_0m的值随着植被生长而发生变化, 夏季植被生长旺盛, 植被密度和高度增大, 因此粗糙度较大。冬季基本没有植被影响, 只体现出微尺度地形的作用。Amdo站的植被高度在四个站中最稳定, 粗糙度的季节变化不明显。

图 3(Figure 3)

图 3 2008年典型站点的动力粗糙度z0m日平均估算值及趋势线 (a) BJ, (b) MS3478, (c) D105, (d) Amdo Figure 3 Daily average estimation and trend line of the roughness length for momentum of typical sites in 2008

图 4(Figure 4)

图 4 2008年典型站点的热传输附加阻尼kB-1日平均估算值及趋势线 (a) BJ, (b) MS3478, (c) D105, (d) Amdo Figure 4 Daily average estimation and trend line of the additional damping of heat transfer of typical sites in 2008

图 5(Figure 5)

图 5 2008年典型站点的热力粗糙度z0h日平均估算值及趋势线 (a) BJ, (b) MS3478, (c) D105, (d) Amdo, 图 5a中“+”和虚线分别为涡动相关法计算值和趋势线 Figure 5 Daily average estimation and trend line of the roughness length for heat of typical sites in 2008. In Fig. 5a, the "+" and dashed lines are values and trend line calculated by the eddy correlation method

从2008年典型站点的热传输附加阻尼kB^-1日平均估算值及趋势线(图 4)可以看到, kB^-1的值在夏季较小, 冬季较大, 与粗糙度相反。kB^-1受很多因素影响, 包括微尺度地形和土壤特征、植被、风速以及温度因子。夏季偏小的原因可能是植被覆盖度和密度的增大造成风速的减小, 从而导致摩擦速度u_*和雷诺数Re减小, 另外, 夏季地-气温差大也是原因之一。可以看到冬季Amdo站的值最大, BJ和MS3478站的值也较大, 而D105最小, 并且季节变化也最小。由于冬季可以不考虑植被作用, 因此原因可能是受到微尺度地形、风速和地-气温差不同程度的影响。基于z_0m和kB^-1的估算值, 得到了四个试验站的热力学粗糙度z_0h(图 5)。z_0h的值基本在夏季较大, 冬季较小, 其中Amdo站较为特殊, 可能受唐古拉山地形和高植被的影响, 中性大气层结条件较弱, 因此点的分布不够集中。而且Amdo站的值大于其他三个站, 反映了该站植被和地形更为复杂的特点。为检验本文kB^-1方案估算的z_0h的可靠性, 给出BJ站的涡动相关观测结果(图 5(a))。通过对比发现, 本文估算结果接近涡动相关观测结果, 表明本文对于kB^-1方案的改进是合理的, 该方案可以有效应用于梯度观测资料中, 其中误差主要来自植被类型和高度的估算等。

通过估算的2008年7月青藏高原z_0m、kB^-1和z_0h以及对应的EVI(图 6)可以看出, 夏季在高原东南部z_0m的值明显大于西北部, 东南部边缘可达0.06 m左右, 东南大部分地区为0.03 m左右, 西北部小于0.02 m(图 6b)。对照EVI(图 6a)来看可知主要原因很可能与高原植被分布状况有关。由于青藏高原四周存在着差异极端显著、对比十分强烈的气候带和植被地理区域, 因此该区域的植被区系较为复杂。资料显示(吴征镒, 1995), 高原东南边缘(海拔3000~4000 m)发育着以森林为代表的山地垂直带植被; 向西、向北, 地势升高, 海拔一般为4000~4500 m, 代之而起的是高寒灌丛和高寒草甸; 继续向北、向西, 平均海拔达4500~5000 m, 大面积分布着高寒草原和高寒荒漠植被; 在海拔较低的藏南谷地, 分布着温性草原和温性干旱落叶灌丛植被; 在青藏高原的最西北部, 即喀喇昆仑山与昆仑山之间的山原和湖盆区, 平均海拔在5000 m以上, 有大面积的冻土分布, 发育着高寒荒漠植被。这样复杂的植被和地形特征也造成了粗糙度分布的不均匀, 整体来看南部大于北部, 东部大于西部。

图 6(Figure 6)

图 6 2008年7月青藏高原EVI(a)、z0m (b)、kB-1(c)和z0h(d)分布 Figure 6 Distribution of EVI (a), z0m (b), kB-1(c) and z0h(d) over Qinghai-Tibetan Plateau in July 2008

7月kB^-1(图 6c)的值基本上都在3~7, 西北部基本大于5, 东南部基本小于5。kB^-1空间分布基本与EVI呈反相关, 原因可能是受风速的影响。kB^-1与风速呈正相关关系, 本研究应用的ECMWF风速在夏季呈西北高、东南低的分布, 导致kB^-1西北高东南低的分布形势。高原上z_0h(图 6d)的分布型与z_0m类似, 也具有南部大于北部的特征, 但比z_0m大约小一个量级, 数值基本上小于0.001 m。

本文计算了青藏高原2000-2014年的月平均热力输送系数, 并选取四个月(4月、7月、10月、12月)为春夏秋冬各季节代表月份, 用其15年平均值来分析青藏高原的热力输送系数的时空分布特征(图 7)。在4月(图 7a)即春季高原整体C_H较小, 除南部边缘地区外基本上都小于0.002; 7月(图 7b)即夏季整个高原的C_H显著增大, 基本在0.002~0.004, 期间C_H值的变化与z_0m的变化趋势密切相关; 10月(图 7c)即秋季随着高原地面温度的下降, 整个高原的C_H值也随之减小, 但南部依然大于北部地区; 12月(图 7d)即冬季的C_H值明显减小, 并且小于春季的值。总体而言, 青藏高原的热力输送系数的季节变化特征是夏季>秋季>春季>冬季。

图 7(Figure 7)

图 7 20002014年4月(a)、7月(b)、10月(c)、12月(d)青藏高原月平均热力输送系数的分布 Figure 7 The distribution of themonthly mean heat transfer coefficient of January (a), April (b), July (c) and October (d) over Qinghai-Tibetan Plateau from 2000 to 2014

生态地理区域划分同时涵盖了气候要素特征和植被信息, 因此研究不同生态地理区域的气候因素的差异是有必要的。在这里, 基于我国生态地理区域系统的划分(张戈丽等, 2010), 在青藏高原上选取三个代表性区域: Ⅰ区(30°N-32°N, 97°E-100°E)位于高原温带湿润半湿润地区, 属于川西藏东高山深谷针叶林区; Ⅱ区(30°N-32°N, 87°E-90°E)位于高原亚寒带半干旱地区, 属于羌塘高原湖盆高寒草原区; Ⅲ区(36°N-38°N, 94°E-97°E)位于高原亚寒带干旱地区, 属于昆仑高山高原高寒荒漠区。为进一步研究青藏高原不同地区的差异, 图 8给出了2000-2014年多年夏季平均的EVI分布, 并标注了所选取的三个特征区域的位置。因为植被指数能够表征区域植被覆盖的整体状况, 从图上也可以看到不同生态地理区域的EVI值具有显著差异, 从Ⅰ区到Ⅱ区到Ⅲ区EVI值逐渐递减, 即植被覆盖状况随着水热条件变差而变差。

图 8(Figure 8)

图 8 2000-2014年青藏高原夏季平均EVI分布及选取特征生态区域 Figure 8 The distribution of average EVI in summer of the Qinghai-Tibetan Plateau during 2000-2014 and the selected characteristic ecological regions

仍选取4月(图 9a)、7月(图 9b)、10月(图 9c)、12月(图 9d)作为春夏秋冬各季节代表月份, 计算了2000-2014年各月的热力输送系数平均值, 来研究其在不同地理生态区域不同季节的时间变化特征。从不同分区的C_H年际变化来看, 在各个季节都是Ⅰ区即湿润半湿润地区的C_H最大, 而对于Ⅱ区和Ⅲ区即半干旱地区和干旱地区而言, 在春季二者的值大小近似, Ⅲ区的C_H稍大; 夏季和秋季都是Ⅱ区的C_H大, 并且在夏季Ⅱ区的值甚至与Ⅰ区接近, 说明该区域夏季植被生长状况较好; 而冬季是Ⅲ区即干旱地区的C_H比Ⅱ区的值大, 这可能与其地形及土壤状况有关, 造成粗糙度较大, 并且干旱地区的地-气温差较大也是原因之一。同时三个区域的值都是夏季较大, 冬季较小。从年际变化来看, 没有明显的变化趋势。不同生态地理区域的热力输送系数的差异, 原因来自于植被高度、密度等状况以及地表的热力特征的差异。植被高度越高, 密度越大, 地表的热力特征越强, 造成热力粗糙度越大, 从而热力输送系数越大。

图 9(Figure 9)

图 9 2000-2014不同生态地理区域4月(a)、7月(b)、10月(c)、12月(d)月均热力输送系数变化 Figure 9 Variation of monthly mean heat transfer coefficient in different eco-geographical regions of April (a), July (b), October (c) and December (d) from 2000 to 2014

选取BJ站的湍流观测资料对本文的估算方案进行合理性检验。因为涡动相关法被认为是最准确的方法, 因此利用涡动相关法得到的C_H值对本文的估算方案得到的结果进行检验(图 10)。从图 10中可以看出, 二者具有较好的相关性, 因此本文的结果较为准确。表 3给出了近年来对青藏高原不同地区的C_H在中性大气条件下的估算值, 改则、当雄、昌都、北麓河四个站分别位于青藏高原的西、南、东、北部, 下垫面属性也各不相同, 因此它们的值具有一定的代表性。其中前三个站是周明煜等(2000)通过涡动相关法的得到的, 北麓河站是钱泽雨等(2005)通过廓线-通量法得到的。将本文的计算值与之对比, 当雄和北麓河站的结果非常接近, 而改则和昌都站偏差较大, 其原因可能和该研究所采用的ECWMF再分析数据的不确定性有关, 再分析风速资料在高原地区可能存在偏差, 导致部分地区的计算值偏离实际值, 同时其他人的结果可能也存在不确定性。

图 10(Figure 10)

图 10 涡动相关法与本文计算的BJ站2008年热力输送系数的月平均值拟合关系 Figure 10 The fitting relation between the heat transfer coefficient calculated by the eddy covariance technique and method in this paper of BJ station in 2008

表 1(Table 1)

表 3 青藏高原部分地区热力输送系数中性计算值 Table 3 The heat transfer coefficient in parts of Tibet Plateau in neutral atmospheric stratification 站点 位置 下垫面属性 文献 本文 改则 32.18°N, 84.03°E 荒漠沙石 2.15×10-3 1.2×10-3 当雄 30.48°N, 91.1°E 开阔干河谷 1.5×10-3 1.5×10-3 昌都 31.15°N, 97.17°E 河谷草地 4.7×10-3 1.6×10-3 北麓河 34.85°N, 92.94°E 稀疏短草沙壤土 1.37×10-3 1.4×10-3   注:改则, 当雄, 昌都周明煜等(2000), 北麓河钱泽雨等(2005)

表 3 青藏高原部分地区热力输送系数中性计算值 Table 3 The heat transfer coefficient in parts of Tibet Plateau in neutral atmospheric stratification

在使用廓线-通量法计算热力输送系数的过程中, 空气动力学粗糙度和热力学粗糙度的计算是否准确直接影响了热力输送系数的准确性, 而热力学粗糙度又与kB^-1有很大关系。本文基于地形特征因子和植被特征因子, 结合地理高程、植被叶面积指数、植被类型与植被覆盖度等多种资料建立了动力学粗糙度、热传输附加阻尼的估算模型, 得到了青藏高原非均匀地形和植被背景下的空气动力粗糙度和热力输送系数。

(1) 在对粗糙度的估算中, 对原始廓线公式进行修改, 增加了EVI和植被高度表征的植被作用项, 更能体现植被覆盖度及类型变化导致的粗糙度的变化。在对热传输附加阻尼的估算中, 对前人的方案进行改进, 使用地-气温差修正, 得到了更接近实际值的kB^-1方案。

(2) 通过分析青藏高原相关参数的空间分布特征, 发现动力粗糙度、热传输附加阻尼、热力粗糙度的空间分布与EVI密切相关, 夏季动力粗糙度和热力粗糙度的分布呈东南大, 西北小的特征, 而热传输附加阻尼基本呈相反的空间型。各参数在高原不同地区差异较大。

(3) 通过对四个典型试验站点相关参数的季节演变规律分析, 动力粗糙度和热力粗糙度都是夏季大冬季小的趋势, 而热传输附加阻尼相反。受到微尺度地形、植被、风速以及地-气温差等多种因素影响, 不同地点的数值大小及趋势程度都存在差异, 体现出高原参数估算中的非均匀性和复杂性。

(4) 青藏高原热力输送系数的空间分布特征主要有两个特点, 一是南北差异, 南部大于北部, 二是东西差异, 东部大于西部。季节演变特征主大致为夏季>秋季>春季>冬季。对于不同的生态地理区域, 热力输送系数的年际变化差异较大, 主要体现为:湿润半湿润地区最大, 半干旱地区与干旱地区较小。导致空间差异较大的根本原因是植被高度、密度和地表的热力特征, 高的植被高度、密度和地表的热力特征对应高的热力粗糙度和热力输送系数。

本文仅讨论了中性大气下的热力输送系数, 但热力输送系数对层结稳定度的变化较为敏感, 需要接触稳定度修订因子进行参数化改进。然而在高原区域尺度上, 再分析资料的不确定性较大, 不能获取大气层结稳定性的空间分布结果, 因此, 非中性大气层结下参数的计算有待完善。但总体来讲, 本文的方法可以获得中性层结下非均匀网格点上连续的热力输送系数, 这对于提高高原区域尺度的模式模拟的精度有一定的意义, 也为其他条件下热力输送系数的估算以及感热通量的研究提供参考依据。