Plateau Meteorology

Plateau Meteorology >

2020 , Vol. 39 >Issue 4: 706 - 718

DOI: https://doi.org/10.7522/j.issn.1000-0534.2020.00022

Temporal and Spatial Characteristics of Ground Surface Soil Heat Flux over the Qinghai-Tibetan Plateau

  • Cheng YANG , 1, 2 ,
  • Tonghua WU , 1 ,
  • Jimin YAO 1 ,
  • Ren LI 1 ,
  • Changwei XIE 1 ,
  • Guojie HU 1 ,
  • Xiaofan ZHU 1 ,
  • Junming HAO 1, 2 ,
  • Jie NI 1, 2 ,
  • Xiangfei LI 1, 2 ,
  • Wensi MA 1, 2 ,
  • Amin WEN 1, 2 ,
  • Chengpeng SHANG 1, 2
Expand
  • 1. Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, State Key Laboratory of Cryospheric Sciences, Lanzhou 730000, Gansu, China
  • 2. University of Chinese Academy of Science, Beijing 100049, China

Received date: 2019-11-18

  Revised date: 2020-04-24

  Online published: 2020-08-28

Fold

Highlights

We have collected the measured data of 9 observation stations on the Qinghai-Tibetan Plateau.The measured data of these sites on the Qinghai-Tibetan Plateau were used to analyze the seasonal and daily variation characteristics of ground surface soil heat flux, and we discussed the reasons for the differences in the temporal characteristics of ground surface soil heat flux between stations.Then MODIS data (including the Terra Moderate Resolution Imaging Spectroradiometer vegetation indices version 6 data MOD13Q1 and the surface spectral reflectance of Terra Moderate Resolution Imaging Spectroradiometer bands 1 through 7 product MOD09CMG), daily land surface temperature dataset in Western China under all-sky condition with 1 km spatial resolution, and assimilation data (ITPCAS-SRad and ITPCAS-LRad) were applied to simulate the spatial distribution of ground surface soil heat flux over the Qinghai-Tibetan Plateau by the scheme Ma on 12 July and 16 October 2014, 1 January and 7 April 2015.The results shown that: The amplitude of ground surface soil heat flux varies with the season, the summer is larger, the winter is the smallest, and the amplitude between the stations is different from the underlying surface.The higher the coverage of the underlying surface, the smaller the amplitude; ground surface soil heat flux is positive in spring and summer and throughout the year, while ground surface soil heat flux in autumn and winter is negative.Ground surface soil heat flux is greater than 0, which means that the underground soil absorbs heat from the surface, otherwise, the underground soil releases heat.Ground surface soil heat flux on the the Qinghai-Tibetan Plateau shows an obvious daily variation curve of inverted U-shape, and its change of the night is relatively flat compared with the daytime.There is a significant seasonal difference in the duration of the daily variation curve of ground surface soil heat flux.The order of the four seasons is: summer > spring > autumn >winter.There is a good positive correlation between the spatial distribution characteristics of ground surface soil heat flux and that of land surface temperature on the the Qinghai-Tibetan Plateau.Station data shows that for every 1 °C increase in land surface temperature, ground surface soil heat flux will increase by 2 ~ 5 W·m-2.This study further deepens the understanding of ground surface soil heat flux over the Qinghai-Tibetan Plateau, in order to provide basis and help for the prediction of frozen soil evolution, accurate estimation of evapotranspiration and geogas interaction.

Cite this article

Cheng YANG , Tonghua WU , Jimin YAO , Ren LI , Changwei XIE , Guojie HU , Xiaofan ZHU , Junming HAO , Jie NI , Xiangfei LI , Wensi MA , Amin WEN , Chengpeng SHANG . Temporal and Spatial Characteristics of Ground Surface Soil Heat Flux over the Qinghai-Tibetan Plateau[J]. Plateau Meteorology, 2020 , 39(4) : 706 -718 . DOI: 10.7522/j.issn.1000-0534.2020.00022

1 引言

青藏高原（简称高原）约占中国陆地面积的四分之一, 平均海拔超过4000 m(马伟强等, 2005; 肖瑶等, 2011; 李黎等, 2019)。因其海拔高而形成独特的气候环境, 高原热力、 动力作用以及地-气间的物质能量交换过程对我国、 亚洲乃至全球的气候变化均有重大影响, 被称为地球第三极（叶笃正等, 1979; 吴国雄, 2004）。仅高原发育的多年冻土就有约1.06×106 km2, 占整个高原面积的40.2%（Zou et al, 2017）, 随着全球变暖, 人们非常关注气候变暖对冻土变化的影响（Yang et al, 2019）。表层土壤热通量（G 0）是阐明其影响的切入点之一。此外, G 0作为估算蒸散发的输入参数, 对精确估算蒸散发也至关重要。
20世纪70年代以来, 中外科学家多次进行关于高原的气象科学试验, 如两次青藏高原大气科学实验（QXPMEX, TIPEX）（章基嘉等, 1988; 陶诗言等, 1998）、 全球能量水分循环亚洲季风之青藏高原试验研究（GEWEX/GAME-Tibetan）、 “全球协调加强观测计划之亚澳季风青藏高原试验”（CEOP/CAMP-Tibetan）和黑河地区地一气相互作用野外观测实验（HEIFE）, 获得了大量的宝贵数据资料（仲雷等, 2006; 肖瑶等, 2011）, 极大地提高了人们对高原G 0的认识（胡隐樵等, 1994; Tsukamoto et al, 1995; 季国良等, 1995, 1999; 沈志宝等, 2002; 马伟强等, 2005）。李明财等（2008）发现土壤热通量的变化与净辐射直接相关, 但由于植被类型的不同, 两者同步性存在明显差异; 刘宏谊等（2009）指出土壤热通量与地表反照率反相关; 肖瑶等（2011）在高原腹地的研究表明G 0的变化主要受地气温差的影响, 同时下垫面状况也起到重要作用; Ma（2006）借助Landsat-5（7）TM（ETM）和NOAA/AVHRR数据反演了黑河实验区、 敦煌试验区、 GAME-Tibetan和CAMP-Tibetan试验区的G 0估算值与实测值的误差率不超过10%。仲雷等（2011）利用AVHRR和MODIS两种卫星资料发现西藏中部“一江两河”地区土壤热通量的空间分布与地表特征参数（植被指数、 植被覆盖度和地表反照率）的空间分布具有很好的对应关系; Zhong et al（2019）通过静止和极轨卫星数据估算高原小时级别的G 0, 均方根误差约为37 W·m-2。然而由于高原特殊的地理条件和恶劣的气候环境, 目前很多研究仅集中在城市和道路沿线的区域, 尽管也有研究给出高原面上的定量结果, 但对整个高原G 0时空分布特征的了解仍具有一定局限性。
本研究利用高原9个站点较新的实测数据, 结合遥感和同化数据, 分析了高原G 0的时空分布特征及其影响因素, 有助于进一步加深对高原G 0的理解, 以期对高原冻土演变预测和蒸散发的精确估算, 甚至对大、 中尺度陆面过程在高原G 0参数化方案挑选与改进及地气相互作用的深入研究有所帮助。

2 资料选取与方法介绍

2.1　研究区概况

利用的9个自动气象站在高原的空间分布如图1所示, 图1背景图为青藏高原冻土空间分布图（Zou et al, 2017）。安多（Amdo）观测场地势较为平整, 视野开阔, 地面被高原稀疏草甸覆盖; 狮泉河（Shiquanhe）位于青藏高原西部, 气候属高原寒区和寒区季风干旱性气候, 太阳辐射强, 地面风速大, 空气干燥, 雨量少, 昼夜温差大(段丽君等, 2017), 下垫面稀疏而短, 一般只有几厘米高; 唐古拉（TGLMS）位于连续多年冻土区, 被高寒草甸覆盖; 西大滩（XDTMS）位于高原上的岛状冻土区, 靠近多年冻土区北部下限; 班戈站（BGGRASS）周围观测区域较为平坦, 被高山草原覆盖, 草高2~4 cm; 比如（Biru）主要下垫面类型为高山草原。嘉黎（Jiali）地处高原亚寒带半湿润季风气候区, 日照充足, 辐射强。那曲（Naqu）地处那曲河谷的高寒草原, 处于高原亚寒带的半湿润气候区(朱志鹍等, 2015); 聂荣站（Nierong）受地形影响较大, 气候寒冷干燥, 属高原亚干旱、 半干旱季风气候区, 观测场周围的地形较宽且相对平坦, 下垫面类型为高寒草甸(杨丽薇等, 2017)。
图1 研究区和站点分布（来自Zou et al, 2017）

Fig.1 Range of the Qinghai-Tibetan Plateau (QTP) and the spatial distribution of observed sites (From Zou et al, 2017)

2.2　资料选取

所用数据有观测站点提供的实测数据和反演青藏高原表层土壤热通量所需的区域驱动数据。站点提供降水, 土壤温度和湿度剖面, 地表温度和土壤热通量（表1）, 数据采集间隔为30 min。降水数据用于剔除特殊天气下的数据; 用土壤热通量, 土壤温度和土壤水分计算表层土壤热通量; 地表温度用于后期分析。TGLMS和XDTMS两个站点的观测数据由中国科学院青藏高原冰冻圈研究站提供, 其他7个站点的实测数据可从国家气象信息中心（http: //tipex.data.cma.cn/tipex）下载。此处需要说明的是部分站点没有地表温度的实际观测, 但可以根据史蒂芬-玻尔兹曼公式计算获取, 公式的具体形式和各站点的地表发射率在Yang et al（2020）中有详细介绍。
表1 观测站点的基本情况、 数据有效时段, 观测项目

Table 1 Basic information of observation sites, data duration and observation items

观测站点 位置（经度/ 纬度/海拔） 下垫面 类型 数据有效时段 观测项目 观测深度或高度 时间 分辨率
TGLMS

91°56′E

33°04′N

5100 m

 高寒草甸 2015年1月1日至12月31日 土壤热通量 5, 10, 20 cm 30 min
土壤含水量 5, 10, 20, 40 cm
土壤温度 5, 10, 20, 40 cm
地表温度 1.5 m
降水量 1.5 m
XDTMS

94°08′E

35°43′N

4538 m

 高寒草原 2015年1月1日至12月31日 土壤热通量 5, 10, 20 cm 30 min
土壤含水量 5, 10, 20, 40 cm
土壤温度 5, 10, 20, 40 cm
地表温度 1.5 m
降水量 1.5 m
Amdo

91°37.2'E

32°14.4'N

4695 m

 高寒草原 2014年7月1日至2015年9月30日 土壤含水量 5, 10, 20, 40, 80, 160 cm 30 min
土壤温度 5, 10, 20, 40, 80, 160 cm
土壤热通量 5, 10 cm
降水量 1.5 m
BGGRASS

90°1.68'E

31°25.07'N

4700 m

 高寒草原 2014年7月15日至2015年9月23日 土壤含水量 5, 10, 20, 40, 100 cm 30 min
土壤温度 5, 10, 20, 40, 100 cm
土壤热通量 5, 10 cm
降水量 1.5 m
Biru

93°9.5'E

31°39.8'N

4408 m

 高寒草原 2014年9月14日至2015年10月2日 土壤含水量 5, 10, 20, 40, 100 cm 30 min
土壤温度 5, 10, 20, 40, 100 cm
土壤热通量 5, 10 cm
降水量 1.5 m
Jiali

93°13.95'E

30°38.456'N

4509 m

 高寒草甸

2014年9月17日至2015年6月28日,

2015年9月6日至10月11日

 土壤含水量 5, 10, 20, 40, 100 cm 30 min
土壤温度 5, 10, 20, 40, 100 cm
土壤热通量 5, 10 cm
降水量 1.5 m
Naqu

91°54'E

31°22.2'N

4509 m

 高寒草原

2014年7月1日至9月10日,

2014年9月16日至12月7日,

2015年1月28日至8月1日,

2015年9月1日至9月30日

 土壤含水量 5, 10, 20, 40, 80, 160 cm 30 min
土壤温度 5, 10, 20, 40, 80, 160 cm
土壤热通量 5, 10 cm
降水量 1.5 m
Nierong

92°18.27'E

32°7.33'N

4607 m

 高寒草甸 2014年7月11日至2015年10月21日 土壤含水量 5, 10, 20, 50, 100 cm 30 min
土壤温度 5, 10, 20, 50, 100 cm
土壤热通量 5, 10 cm
降水量 1 m
Shiquanhe

80°6'E

32°29.4'N

4278.6 m

 高寒荒漠 2014年9月16日至2015年9月30日 土壤含水量 5, 10, 20, 40, 80 cm 30 min
土壤温度 5, 10, 20, 40, 80 cm
土壤热通量 5, 10, 20, 40, 80 cm
地表温度 1.5 m
降水量 1.5 m

观测深度/高度属性一栏中, 土壤热通量、 土壤温度和土壤含水量全部是深度, 其余地表参数是高度

所需区域驱动数据包括中国西部逐日1 km空间分辨率全天候地表温度数据集, MODIS数据（MOD09CMG和MOD13Q1）和中国区域地面气象要素驱动数据集。中国西部逐日1 km空间分辨率全天候地表温度数据集（2003 -2018年）V1（http: //data.tpdc.ac.cn/zh-hans/data/d09ebfe3-a8c9-49b5-aa22-1fd484751d70/）提供青藏高原区域的地表温度数据。该数据是一种基于地表温度时间成分分解模型的热红外与被动微波地表温度集成的全天候、 具有较好的精度和图像质量的地表温度产品, 时间跨度为2003 -2018年, 时间分辨率为逐日2次, 空间分辨率为1 km (Zhou et al, 2017; Zhang et al, 2019)。MOD09CMG产品用于反演区域地表反照率, 该产品提供了MODIS波段1~7的地表光谱反射率的估计值, 重采样到5600 m的像素分辨率, 并根据大气条件(如气体、 气溶胶和瑞利散射)进行校正。MOD13Q1提供研究所需的归一化植被指数（NDVI）和红光、 近红外波段的地表反射率, NDVI用于反演区域地表发射率, 红光、 近红外波段的地表反射率反演区域改良的土壤调整植被指数（MSAVI）。上述遥感反演算法在章节2.4中有详细介绍。本研究还使用中国区域地面气象要素驱动数据集(CMFD, http: //tipex.data.cma.cn/tipex)中向下短波辐射（CMFD-SRad）和向下长波辐射（CMFD-LRad）作为模型部分参数的输入。该数据集融合了来自中国气象局(CMA)740个业务站点的观测数据, 以及相应的普林斯顿全球气象强迫数据集(何杰, 2010; Yang et al, 2010; Zhong et al, 2019), 数据空间分辨率为10 km, 时间分辨率为3 h(表2)。本研究中MODIS数据提供或反演获取的参数都被重采样到10 km, 一天之内的地表发射率, 植被指数和地表反照率可以认为变化不大（Zhong et al, 2019）, 选取中国西部逐日1 km空间分辨率全天候地表温度数据集白天的数据, 在青藏高原对应的北京时间范围是14:30 -16:30, 本研究选取14:00和17:00两个时刻强迫数据（CMFD）均值作为向下短波、 长波辐射的数值。
表2 遥感和同化数据的基本信息

Table 2 Basic information of remote sensing and assimilation data

产品名称 包含内容 数据有效时段 时间分辨率 空间分辨率
MOD13Q1 归一化植被指数（NDVI）红光波段和近红外波段的地表反射率

2014年7月12日

2014年10月16日, 2015年1月1日, 2015年4月7日

 16-day 250 m
中国西部逐日1km空间分辨率全天候地表温度数据集 白天/夜间地表温度 2014年7月12日, 2014年10月16日, 2015年1月1日, 2015年4月7日 Daily 1000 m
MOD09CMG 波段1~7的地表反射率 Daily 5600 m
ITPCAS-SRad 向下短波辐射 3-hour 10000 m
ITPCAS-LRad 向下长波辐射 3-hour 10000 m

MOD是Terra卫星的MODIS数据

2.3　站点 G 0的计算

土壤的一维热传导(或扩散)方程(Thermal Diffusion Equation, TDE)为（Yang et al, 2008）:
∂ ρ s c s T ∂ t = ∂ G ∂ z
G = λ s ∂ T ∂ z
式中: t为时间(单位: s); z为土壤深度(向下为正)(单位: m); T为土壤温度(单位: K); c s是土壤热容量(单位: J·kg-1·K-1); λs为土壤热传导系数(单位: W·K-1·m-1); G为土壤热通量(向下为正)(单位: W·m-2)。
对公式(1)积分, 得到表层土壤热通量G 0:
G 0 = G ( z r e f ) + ∫ Z r e f 0 ∂ ρ s c s T ( z ) ∂ t d z
式中: z ref为参考深度(在本研究中参考深度是10 cm); G(z ref)为参考深度处的土壤热通量。
由于青藏高原的近地表土壤存在冻融过程, 所以高原上观测站点G 0的计算也应考虑冰的热容量(Yang et al, 2008; Yao et al, 2011), 故公式(3)的离散形式可以表示为:
G 0 = G 0.10 m + ( ρ d r y c d r y + ρ w c w θ 0.05 m + ρ i c i ϑ i , 0.05 m ) ∂ T 0.05 m ∂ t × 0.10
式中: θ 0.05m为深度0.05 m处的未冻水体积含水量(单位: m3∙m-3); ρ dry c dry为干土的热容量, 约取0.90×106 J∙m-3∙K-1 (Tanaka et al, 2003); ρ w c w=4.2×106 J∙m-3∙K-1为液态水的热容量; ρici=1.89×106 J∙m-3∙K-1为冰的热容量。由于没有蒸散发数据, 本研究计算0.05 m深度土壤含冰量时忽略了蒸散发对5 cm土壤水分的影响。下面是冻融过程0.05 m深度土壤含冰量的计算(Yao et al, 2011):
(1) 当0.05 m深度以上的土壤层处于完全融化阶段时, 此土壤层是没有冰的存在的, 也就是土壤中含冰量为0( ϑ i , 0.05 m = 0);
(2) 当0.05 m深度以上的土壤层处于日冻融循环阶段时, 土壤中的含冰量可以根据对应深度的土壤中相邻时刻之间未冻水含量的差值来求取, 具体公式如下所示(Tom, 1987):
∂ ρ i ϑ i , 0.05 m ¯ ∂ T = ∂ ρ w θ 0.05 m ¯ ∂ T
(3) 当土层10 cm以上完全冻结, ϑi , 0.05m等于日冻融循环阶段最后一天的土壤中含冰量。需要注意的是, 即使土壤冻结, 土壤中的含冰量也是变化的, 但是在冻结阶段, 冰的含量变化很小, 因此在计算中, 土壤中含冰量可以被视为是一个常数(Yao et al, 2011)。
日冻融循环时段可以根据深度为0.05 m处的土壤温度T 0.05m和土壤湿度θ 0.05m来确定(表3)。
表3 站点日冻融循环时段起止对应的时间

Table 3 The start and end date of daily freezing-thawing cycle at sites

站点 年份
2014年 2015年
XDT - 3月25日至5月1日, 10月14日至11月6日
TGL - 3月17日至5月20日, 10月9日至11月8日
Amdo 9月28日至12月9日 4月12日至6月16日
BGGRASS 10月18日至11月20日 2月24日至4月16日
Biru 10月19日至11月23日 3月10日至4月20日
Jiali 11月6日至12月17日 2月23日至4月15日
Naqu 10月20日至12月8日 2月3日至3月29日
Nierong 10月24日至11月19日 2月22日至4月10日
Shiquanhe 10月13日至12月7日 2月5日至4月21日

- 表示无数据

2.4　区域 G 0的计算

当前G 0遥感估算模型主要基于G 0与净辐射R n之间的比值, 譬如利用植被指数来估算这一比值的遥感模型Choudhury（Choudhury et al, 1987）, Clawson（Jackson et al, 1987; Moran et al, 1989）和SEBS（Su, 2002）, 以及除了植被指数还考虑地表温度、 地表反照率的遥感模型SEBAL（Bastiaanssen, 2000）和Ma（Ma et al, 2003）, 各模型的具体介绍见Yang et al(2019)。而且Yang et al(2019)还证明Ma方案在高原的适用性最佳, 故本研究直接采用Ma方案来估算高原的表层土壤热通量G 0。Ma方案的具体形式如下所示:
G 0 = T s α ⋅ ( 0.0087 α ¯ 2 + 0.00454 α ¯ + 0.00029 ) ⋅ ( 1 - 0.964 M S A V I 4 ) ⋅ R n
式中: T s为地表温度（单位: ℃）; α和 α ¯分别是地表反照率的瞬时值和日均值; MSAVI（Modified soil adjusted vegetation index）为改良的土壤调整植被指数。地表温度T s可以由产品中国西部逐日1 km空间分辨率全天候地表温度数据集提供, 但是真正的地表温度值必须通过遥感影像值乘以缩放因子0.01获取(单位: K）。下面详细介绍各地表参数的遥感反演算法。
净辐射R n可以通过公式（7）计算获取（Amaty et al, 2015）:
R n = D S R - U S R + D L R - U L R U S R = α ⋅ D S R U L R = σ ⋅ ε s ⋅ T s 4 + ( 1 - ε s ) ⋅ D L R R n = ( 1 - α ) ⋅ D S R + ε s ⋅ D L R - ε s ⋅ σ ⋅ T s 4
式中: USR, DSR, DLR和ULR分别是向上短波辐射, 向下短波辐射, 向下长波辐射和向上长波辐射(单位: W∙m-2)。
在一天10:00 -16:00时间范围内, 地表反照率α的瞬时值大约等于同一天的地表反照率日均值 α ¯（Yang et al, 2019）。地表瞬时反照率α可以从MODIS产品MOD09CMG中的波段反射率数据和经验公式得出（Liang, 2001）:
α=0.160R 1+0.291R 2+0.243R 3+0.166R 4
+0.112R 5+0.018R 7-0.0015 (8)
式中: R 1、 R 2、 R 3、 R 4、 R 5和R 7为MOD09CMG的第1、 2、 3、 4、 5和7波段的地表反射率。
ε s是地表发射率, 根据以下公式计算得出（覃志豪等, 2004）:
ε s =           0.973 , N D V I < 0.05 0.986 + 0.004 ⋅ f c , 0.05 ≤ N D V I ≤ 0.7         0.99 , N D V I > 0.7
式中: NDVI为归一化植被指数, 可从MODIS植被指数产品MOD13Q1获得, 其时间间隔为16天, 空间分辨率为250 m。植被盖度f c可以根据NDVI估算（Carlson et al, 1997）:
f c = N D V I - N D V I m i n N D V I m a x - N D V I m i n 2
式中: NDVI max和NDVI min分别是完全植被覆盖和裸露土壤的NDVI值。
MSAVI也是由MODIS植被指数产品MOD13Q1通过以下公式计算得出（Qi et al, 1994）:
M S A V I = 2 N I R + 1 - 2 N I R + 1 2 - 8 N I R - R E D 2
式中: RED和NIR是MOD13Q1的RED波段和NIR波段的地表反射率值。

3 结果与讨论

3.1　站点 G 0与 G 5的比较

由于表层土壤热通量无法实际观测, 本研究以TGLMS站为例, 利用最浅的土壤热通量观测值深度5 cm的土壤热通量G 5来检验一下表层土壤热通量G 0的单点计算结果。从图2(a)中可以看出, G 0与G 5的无截距线性关系式中的斜率略大于1, 表明G 0是略大于G 5。2015年表层土壤热通量G 0年平均日变化曲线的振幅要高于实测5 cm土壤热通量G 5的振幅, G 0年平均日变化曲线要先于G 5达到峰值[图2(b）]。这些G 0与G 5比较得出的结论与李明财等（2008）研究研究结果果相一致, 说明计算的表层土壤热通量是符合实际情况的。
图2 TGLMS站2015年（除去降水）表层土壤热通量G 0与5 cm土壤热通量G 5的散点分布和年平均日变化曲线比较

Fig.2 The comparison of scatter diagram and the annual average daily change curves of the surface soil heat fluxes G 0 and 5 cm depth soil heat fluxes G 5 with an interval of 30 min at TGLMS in 2015 (excluding precipitation data)

3.2　 G 0的季节变化特征

选取地处多年冻土区的TGLMS和季节冻土区的Amdo两个站点春季（3 -5月）、 夏季（6 -8月）、 秋季（9 -11月）和冬季（12月至次年2月）(Wu et al, 2013)研究表层土壤热通量的季节变化特征（图3）。
图3 TGLMS和Amdo站点表层土壤热通量G 0的季节变化（时间分辨率是30 min）

Fig.3 Seasonal variation of ground surface soil heat flux G 0 at TGLMS and Amdo sites (Time resolution is 30 min)

表层土壤通量振幅随季节变化（图3）, 夏季较大, 冬季最小, 这是由于表层土壤热通量变化受净辐射的影响。净辐射主要分为潜热、 感热和表层土壤热通量, 而且表层土壤热通量与净辐射极显著正相关（郝雅婕等, 2019）, 影响净辐射季节变化最显著的是太阳高度角(肖瑶等, 2011)。由于高原大气的透明度很高, 削光作用较弱, 到达地面的总辐射（太阳短波辐射）要大于同纬度其他地区（季国良, 1985）, 总辐射随太阳高度角的变化而呈现出特别显著的季节波动, 从而直接影响净辐射季节变化（肖瑶等, 2011）。此外, 地表反照率也是影响净辐射季节变化的一个重要因素, 夏季植被茂盛, 反照率小, 植被秋季变黄, 地表反照率变大, 春季开始返青, 反照率也比较大, 冬季地表植被枯死, 地表反照率最大, 如果有雪, 地表反照率急剧增大, 可以达到0.9（杨成等, 2016）。Amdo表层土壤热通量的振幅明显高于TGLMS(图3), 这可能是由于两站的下垫面不同所致（肖瑶等, 2011）, TGLMS和Amdo分别位于高寒草甸区和高寒草原区, TGLMS的植被覆盖度高于Amdo, TGLMS到达土壤表层的净辐射峰值相对较小, 随之G 0峰值也会变小。
由表4可知, 不同季节土壤与地上气体热量交换模式有明显的差异, 表层土壤热通量G 0在春季和夏季整体为正（G 0的传播方向是向下为正）, 表明能量由地上部分向土壤传递, 土壤为热汇; 而秋季和冬季G 0则为负, 能量由土壤向地上部分传递, 土壤为热源（李明财等, 2008）; 相对全年而言, 绝大部分站点（Jiali除外）的G 0是在由近地表大气向土壤传递热量, 年均值量级不超过5 W·m-2, 单方面来讲, 这不利于高原上大面积多年冻土的存在与发育。多年冻土对全球气候变暖来说像是一个缓冲区, 因为多年冻土活动层的融化过程需要相当的热量, 所以如果冻土退化或消失, 可能会加速气候变暖（Eugster et al, 2000; 谷良雷, 2011）。
表4 各站点地表土壤热通量 G 0的季节和年平均值

Table 4 Seasonal and annual average of ground surface soil heat flux G 0 at each site

站点 G 0
春季 夏季 秋季 冬季 年
XDTMS 9.14 8.62 -3.26 -7.78 4.03
TGLMS 5.12 9.67 -3.34 -7.83 0.97
Amdo 5.35 7.53 -3.61 -6.34 0.78
BGGRASS 6.28 6.11 -2.40 -4.97 1.30
Biru 8.22 7.62 -5.97 -11.45 0.35
Jiali 4.63 - -6.88 -5.74 -
Naqu 8.54 10.65 -6.03 -7.91 2.31
Nierong 3.16 5.51 -1.09 -2.43 1.22
Shiquanhe 3.89 5.24 -4.78 -3.62 0.37

- 表示无数据

3.3　 G 0的日变化特征

为了更清楚地了解表层土壤热通量在不同季节的日变化特征, 分别选取TGLMS和Amdo两个站点每个季节每0.5 h平均的资料进行分析（图4）。结果表明, 四季均呈现明显的倒立“U”字型的日变化曲线, G 0随着日出后太阳高度角的增大而逐渐增大, 一般在12:00 -15:00达到最大值, 然后随着太阳高度角的减小, G 0又趋于减小。G 0在一天之内两次通过零点, 日落后逐渐转为负值。G 0变化曲线在夜间比较平缓, 表明表层土壤热通量在夜晚的变化相对白天而言比较平缓（仲雷等, 2006）。而且G 0日变化曲线为正值的时长存在明显季节差异, 四个季节的顺序是夏季>春季>秋季>冬季（图4）, 最大9 h, 最小6 h。
图4 TGLMS和Amdo两个站点表层土壤热通量G 0的季节平均的日变化

Fig.4 Diurnal variation of seasonal average G 0 at TGLMS and Amdo sites

表层土壤热通量日变化的季节区别也较为明显, 各季节G 0振幅不同, 冬季振幅最小, 夏季振幅最大, 但是Amdo春季振幅最大, 这可能与Amdo夏季净辐射有关。天空的总云量、 云状等气象要素是引起日平均净辐射高频波动的主要因子（肖瑶等, 2011）, 高原夏季季风期, 降水和云含量会增多, 气象要素变化极快, 导致总辐射日变化剧烈且地面接受的总辐射也会降低, 到达地面的净辐射会有所降低, 从而影响净辐射使日平均净辐射也呈现出高频波动且降低的特征（李宏毅等, 2018）。G 0又与净辐射正相关（郝雅婕等, 2019）, 因此Amdo夏季表层土壤热通量振幅会小于春季。
四季中G 0日变化振幅Amdo都大于TGLMS, 而且Amdo的峰值也早于TGLMS, 前面季节变化特征部分已经讨论过, 这可能是与两站下垫面不同有关。由于植被的遮挡, 高寒草甸的TGLMS相对于高寒草原的Amdo, 到达地面的净辐射会有所减少, G 0变化对净辐射的反馈延时会有所增加（李明财等, 2008）。

3.4　 G 0的空间分布特征

考虑实测数据的有效时段和高原G 0的季节特征, 选取了2014年7月12日和10月16日、 2015年1月1日和4月7日来反映高原G 0的空间分布特征。利用估算表层土壤热通量G 0的方法, 将区域参数作为输入驱动, 进而模拟获取所选4天（2014年7月12日和10月16日、 2015年1月1日和4月7日）高原表层土壤热通量的空间分布（图6）。通过G 0的遥感估算值和实测值的比较（图5, 其中表格是圆圈对应的站点和日期）发现, 两者线性关系的斜率是0.89, 均方根误差RMSE是18.78 W·m-2, 估算值是比较接近对应的实测值, MBE（-7.19 W·m-2）小于0, 说明方案Ma存在一定的低估现象, 这与Yang et al（2020）在高原得出的结论一致。方案Ma模拟结果出现低估可能与G 0与净辐射Rn日变化之间存在相位偏移（Zhong et al, 2019）和方案输入的遥感数据精度（Yang et al, 2020）有关。
图5 Ma方案估算表层土壤热通量值G 0_Ma与对应表层土壤热通量实测值G 0_obs的比较

Fig.5 Compare the estimated G 0_Ma by the scheme Ma with the measured G 0_obs

图6 2014年7月12日和10月16日, 2015年1月1日和4月7日表层土壤热通量G 0的空间分布特征

Fig.6 Spatial distribution of ground surface soil heat flux G 0 on 12 July, 16 October in 2014 and on 1 January, 7 April in 2015

所选4天的区域模拟结果（图6）可以反映出青藏高原表层土壤热通量空间分布特征的季节变化, 春夏较大, 冬季最小, 这与在单点得出G 0的季节变化规律一致。而且发现模拟的高原G 0（图6）与对应时刻高原地表温度T s的空间分布（图没有显示）的相关矩阵中的相关系数都大于0.5（表5）, 这说明G 0空间分布受T s的影响显著。由各观测站点的“实测值”表层土壤热通量与地表温度的线性关系（表6）可知, 表层土壤热通量G 0与地表温度T s呈正相关（P<0.01）, 即T s每增加1 ℃, G 0便随之增大2~5 W·m-2（表6中线性关系的斜率）, Yang et al（2019）在高原也得出表层土壤热通量与地表温度的空间分布具有较高的一致性。
表5 四期（2014年7月12日和10月16日、 2015年1月1日和4月7日）区域 G 0_Ma与地表温度 T s的空间相关矩阵

Table 5 Correlation matrix of regional simulated ground surface soil heat flux G 0_Ma by the scheme Ma and land surface temperature T s at four times (on 12 July and 16 October in 2014, on 1 January and 7 April 2015)

区域表层土壤热通量G 0_Ma
2014-07-12 2014-10-16 2015-01-01 2015-04-07
区域地表温度T s 2014-07-12 0.67 - - -
2014-10-16 - 0.70 - -
2015-01-01 - - 0.88 -
2015-04-07 - - - 0.55

- 表示无数据

表6 9个观测站点“实测”表层土壤热通量 G 0与对应地表温度 T s的关系

Table 6 The relationship between the "measured" ground surface soil heat flux G 0 and the corresponding land surface temperature T s at 9 observation stations

站点 线性关系 相关系数R 显著性检验 N
TGLMS G 0=2.26T s+9.69 0.71 P<0.01 11290
XDTMS G 0=2.94T s+11.29 0.64 P<0.01 15873
Amdo G 0=3.28T s-1.46 0.57 P<0.01 12698
BGGRASS G 0=2.90T s-8.64 0.68 P<0.01 20315
Biru G 0=2.75T s-0.32 0.59 P<0.01 17513
Jiali G 0=1.88T s-3.12 0.60 P<0.01 14800
Naqu G 0=5.04T s-18.46 0.73 P<0.01 9442
Nierong G 0=2.13T s-1.15 0.55 P<0.01 21451
Shiquanhe G 0=2.16T s-25.24 0.60 P<0.01 12048

各个站点数据时段如表1, 时间分辨率是30 min, N是样本数

4 结论

通过青藏高原9个站点的实测数据, 并结合遥感（MODIS和中国西部逐日1 km空间分辨率全天候地表温度数据集）和同化数据（ITPCAS-SRad和ITPCAS-LRad）, 对青藏高原表层土壤热通量的时空分布特征进行了分析, 得出以下主要结论:
(1） 表层土壤通量G 0振幅随季节变化, 夏季较大, 冬季最小, 站点之间振幅不同与下垫面有关, 下垫面植被覆盖度越高, 振幅越小; G 0在春季和夏季以及全年整体为正, 而秋季和冬季G 0则为负, G 0的年均值不超过5 W·m-2。
(2） 高原表层土壤通量呈现明显倒立的“U”型的日变化曲线, 夜晚的变化相对白天而言比较平缓; G 0日变化曲线为正值的时长存在明显季节差异, 四个季节的顺序是: 夏季>春季>秋季>冬季。
(3） 高原G 0的空间分布特征与高原地表温度的空间分布有较好的正相关, 站点数据显示地表温度每增加1 ℃, G 0便随之增大2~5 W·m-2。

非常感谢中国科学院青藏高原冰冻圈观测研究站, 寒区旱区科学数据中心（http: //westdc.westgis.ac.cn）和国家气象信息中心（http: //www.nmic.cn）提供的站点实测数据, 感谢美国宇航局（https: //modis.gsfc.nasa.gov）提供的MODIS数据和国家青藏高原科学数据中心（http: //data.tpdc.ac.cn/）提供的中国西部逐日1 km空间分辨率全天候地表温度数据集和同化数据。

References

Publishing order | Descend order by publishing year | Descend order by cited within
Amaty P M, Ma Y M, Cunbo H, al et, 2015.Recent trends (2003 -2013) of land surface heat fluxes on the southern side of the central Himalayas, Nepal[J].Journal of Geophysical Research: Atmospheres, 120, 11957-11970.DOI: 10.1002/2015jd023510.

Bastiaanssen W G M, 2000.SEBAL-based sensible and latent heat fluxes in the irrigated Gediz Basin, Turkey[J].Journal of Hydrology, 229, 87-100.DOI: 10.1016/s0022-1694(99)00202-4.

Carlson T N, Ripley D A, 1997.On the relation between NDVI, fractional vegetation cover, and leaf area index[J].Remote Sensing of Environment, 62, 241-252.DOI: 10.1016/s0034-4257(97)00104-1.

Choudhury B J, Idso S B, Reginato R J, 1987.Analysis of an empirical model for soil heat flux under a growing wheat crop for estimating evaporation by an infrared-temperature based energy balance equation[J].Agricultural and Forest Meteorology, 39, 283-297.DOI: 10.1016/0168-1923(87)90021-9.

Eugster W, Rouse W, Pielkesr R A, al et, 2000.Land-atmosphere energy exchange in arctic tundra and boreal forest: Available data and feedbacks to climate[J].Global Change Biology, 6(1): 84-115.DOI: 10.1046/j.1365-2486.2000.06015.x.

Jackson R D, Moran M S, Gay L W, al et, 1987.Evaluating Evaporation from field crops using airborne radiometry and ground-based meteorological data[J].Irrigation Science, 8 (2): 81-90.DOI: 10.1007/BF00259473.

Liang S, 2001.Narrowband to broadband conversions of land surface albedo I Algorithms[J].Remote Sensing of Environment, 76 (2), 213-238.DOI: 10.1016/s0034-4257(00)00205-4.

Moran M S, Jackson R D, Raymond L H, al et, 1989.Mapping surface energy balance components by combining Landsat thermatic mapper and ground-based meteorological data[J].Remote Sensing of Environment, 30, 77-87.DOI: 10.1016/0034-4257(89)90049-7.

Ma Y M, 2006.Determination of regional surface heat fluxes over heterogeneous landscapes by integrating satellite remote sensing with boundary layer observations[D].The Netherlands: Wageningen University, 1-181.

Ma Y M, Ishikawa H, Tsukamoto O, al et, 2003.Regionalization of surface fluxes over heterogeneous landscape of the Tibetanan Plateau by using satellite remote sensing data[J].Journal of the Meteorological Society of Japan, 81(2): 277-293.DOI: 10.2151/jmsj.81.277.

Osterkamp T, 1987.Freezing and thawing of soils and permafrost containing unfrozen water or brine[J].Water Resources Research, 23(12): 2279-2285.DOI: 10.1029/wr023i012p02279.

Qi J G, Chehbouni A, Huete A, al et, 1994.A modified soil adjusted vegetation index[J].Remote Sensing of Environment, 48 (2), 119-126.DOI: 10.1016/0034-4257(94)90134-1.

Su Z B, 2002.The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes[J].Hydrology and Earth System Sciences, 6(1), 85-100.

Wu T H, Lin Z, Ren L, al et, 2013.Recent ground surface warming and its effects on permafrost on the central Qinghai-Tibetan Plateau[J].International Journal of Climatology, 33: 920-930.DOI: 10.1002/joc.3479.

Tanaka K, Tamagawa I, Ishikawa H, al et, 2003.Surface energy budget and closure of the eastern Tibetanan Plateau during the GAME-Tibetan IOP 1998[J].Journal of Hydrology, 283: 169-183.

Tsukamoto O, Sahashi K, Wang J M, 1995.Heat budget and evaptranspiration at an oasis surface surrounded by desert[J].Journal of the Meteorological Society of Japan, 73(5): 925-935.

Yang C, Wu T H, Wang J M, al et, 2019.Estimating surface soil heat flux in permafrost regions using remote sensing-based models on the Northern Qinghai-Tibetanan Plateau under clear-sky conditions[J].Remote Sensing, 11 (4), 416.DOI: 10.3390/rs11040416.

Yang C, Wu T H, Yao J M, al et, 2020.An assessment of using remote sensing-based models to estimate ground surface soil heat flux on the Tibetanan Plateau during the freeze-thaw process[J].Remote Sensing, 12, 501.

Yang K, Wang J M, 2008.A temperature prediction-correction method for estimating surface soil heat flux from soil temperature and moisture data[J].Science in China (Series D: Earth Sciences), 05: 99-107.

Yang K, He J, Tang W J, al et, 2010.On downward shortwave and longwave radiations over high altitude regions: Observation and modeling in the Tibetanan Plateau[J].Agricultural and Forest Meteorology, 150(1), 38-46.DOI: 10.1016/j.agrformet.2009. 08.004.

Yao J M, Zhao L, Gu L L, al et, 2011.The surface energy budget in the permafrost region of the Tibetanan Plateau[J].Atmospheric Research, 102(4): 394-407.DOI: 10.1016/j.atmosres.2011. 09.001.

Zhang X D, Zhou J, Frank-Michael G, al et, 2019.A method based on temporal component decomposition for estimating 1-km all-weather land surface temperature by merging satellite thermal infrared and passive microwave observations[J].IEEE Transactions on Geoscience and Remote Sensing, 57: 4670-4691.

Zhong L, Ma Y M, Hu Z Y, al et, 2019.Estimation of hourly land surface heat fluxes over the Tibetanan Plateau[J].Atmospheric Chemistry and Physics, 19: 5529-5541.

Zhou J, Zhang X D, Zhan W F, al et, 2017.A thermal sampling depth correction method for land surface temperature estimation from satellite passive microwave observation over barren land[J].IEEE Transactions on Geoscience and Remote Sensing, 55: 4743-4756.

程国栋, 赵林, 2000.青藏高原开发中的冻土问题[J].第四纪研究, 20(6): 521-531．

段丽君, 段安民, 胡文婷, 等, 2017.2014 年夏季青藏高原狮泉河与林芝降水低频振荡及陆—气过程日变化特征[J].大气科学, 41 (4): 767-783.

谷良雷, 2011.青藏高原季节冻土区和多年冻土区地表能量收支和蒸散发对比研究[D].兰州: 中国科学院寒区旱区环境与工程研究所, 1-136.

郝雅婕, 邓巧玲, 王艳霞, 等, 2019.元江干热河谷稀树灌丛土壤热通量特征[J].西北林学院学报, 34(5): 23-28.DOI: 10.3969/j.issn.1001-7431.2019.05.04.

何杰, 2010.中国区域高时空分辨率地面气象要素数据集的建立[D].北京: 中国科学院青藏高原研究所, 1-78.

胡隐樵, 高由禧, 王介民, 等, 1994.黑河实验(HEIFE)的一些研究成果[J].高原气象, 13(3): 225-236.

季国良, 1999.青藏高原能量收支观测实验的新进展[J].高原气象, 18(3): 333-340．

季国良, 1985.青藏高原主体的大气透明度特征[J].高原气象, 4(2): 122-129.

季国良, 吕兰芝, 1995.藏北高原太阳辐射能收支的季节变化[J].太阳能学报, 16(4): 340-346.

李黎, 吕世华, 范广洲, 2019.夏季青藏高原地表能量变化对高原低涡生成的影响分析[J].高原气象, 38(6): 1172-1180.DOI: 10.7.522 /j.issn.1000-0534.2018.00154.

李明财, 罗天祥, 郭军, 等, 2008.藏东南高山林线冷杉原始林土壤热通量[J].山地学报, 26(4): 490-495.DOI: 10.3969/j.issn.1008-2786.2008.04.017.

李宏毅, 肖子牛, 朱玉祥, 2018.藏东南地区草地下垫面湍流通量和辐射平衡各分量的变化特征[J].高原气象, 37(4): 923-935.DOI: 10.7252/j.issn.1000-0534.2017.2007.

刘宏谊, 杨兴国, 张强, 等, 2009.敦煌戈壁冬夏季地表辐射与能量平衡特征对比研究[J].中国沙漠, 29(3): 558-565.

马伟强, 马耀明, 胡泽勇, 等, 2005a.藏北高原地区辐射收支和季节变化与卫星遥感的对比分析[J].干旱区资源与环境, 1: 110-116.DOI: CNKI: SUN: GHZH.0.2005-01-021.

马伟强, 马耀明, 李茂善, 等, 2005b.藏北高原地区地表辐射出支和能量平衡的季节变化[J].冰川冻土, 27(5): 673-679.http: //

沈志宝, 成天涛, 王可丽, 2002.青藏高原地面-对流层系统的能量收支[J].高原气象, 21(6): 20-25.

覃志豪, 李文娟, 徐斌, 等, 2004.陆地卫星TM6波段范围内地表比辐射率的估计[J].国土资源遥感, 15(3): 28-32.

陶诗言, 陈联寿, 徐祥德, 等, 1998.第二次靑藏高原大气科学试验理论硏究进展(一) [M].北京: 气象出版社, 1-348.

吴国雄, 2004.我国青藏高原气候动力学研究的近期进展[J].第四纪研究, 24(1): 1-9.

翁笃鸣, 高歌, 2001.利用卫星资料试作青藏高原地表净辐射场的气候反演[J].气象科学, 21 (2): 162-168.DOI: 10.3969/j.issn.1009-0827.2001.02.005.

肖瑶, 赵林, 李韧, 等, 2011.青藏高原腹地高原多年冻土区能量收支各分量的季节变化特征[J].冰川冻土, 33(5): 79-85.

杨成, 姚济敏, 赵林, 等, 2016.藏北高原多年冻土区地表反照率时空变化特征[J].冰川冻土, 38(6): 1518-1528.

杨丽薇, 高晓清, 惠小英, 等, 2017．青藏高原中部聂荣亚寒带半干旱草地近地层湍流特征研究[J].高原气象, 36(4): 875-885.DOI: 10.7522 /j.issn.1000-0534.2016.00089.

叶笃正, 高由禧, 等, 1979.青藏高原气象学[M].北京: 科学出版社, 30-55.

章基嘉, 朱抱真, 朱福康, 等, 1988.青藏高原气象学进展[M].北京: 科学出版社, 92-123.

仲雷, 马耀明, 马伟强, 等, 2006.珠峰北坡地区近地层大气湍流与地气能量交换特征[J].地球科学进展, 12: 84-94.

仲雷, 马耀明, 马伟强, 等, 2011.西藏中部“一江两河”地区地表通量的卫星遥感估算[J].冰川冻土, 33(2): 309-317.

朱志鹍, 马耀明, 胡泽勇, 等, 2015.青藏高原那曲高寒草甸生态系统CO2净交换及其影响因子[J].高原气象, 34(5): 1217-1223.DOI: 10.7522/j.issn.1000-0534.2014.00135.

Options
Outlines

/

〈 〉

The system is designed and developed by Beijing magtec Technology Development Co., Ltd. with technical support: support@magtech.com.cn