1 引言

土壤温度是陆气相互作用中的关键物理量(陈海山等, 2010; 郭东林等, 2010; 李剑铎等, 2013), 可以通过影响地表能量和水分收支的变化来影响气候变化。浅层土壤温度的变动在还未传递到深层土壤前, 可向大气中释放能量, 深层土壤还可通过逐步释放能量影响浅层土壤。因此, 浅层土壤作为交界面, 其温度对地气间的热量输送非常关键(Zhou et al, 2004; 张伟等, 2012; 孟现勇等, 2017)。土壤温度时空变异大, 地面站点观测难以获得空间上连续的数据(Moiwo et al, 2012), 卫星遥感技术又往往因空间分辨率、卫星过境时间、反演精度和积雪覆盖等因素的影响, 有较大局限性(马思源等, 2016)。陆面模式模拟可以充分利用观测和遥感数据, 并能获取时空连续、高分辨率、不同深度的土壤温度格点化资料, 受到研究者的广泛重视。目前已有较多利用陆面模式模拟土壤温度的相关研究, 如夏坤等(2011)发现CLM3.0模式能模拟出青藏高原东北部各层土壤的冻融变化趋势; 熊建胜等(2014)发现CLM模式能较好地模拟出玛曲的土壤温湿度、湍流通量等的变化趋势及特征; 杨扬等(2016)通过比较NOAH、SHAW、CLM三个陆面过程模式在半干旱区的表现后发现CLM模式模拟的浅层土壤温度最优; Fang et al(2016)通过改进CLM4.5中土壤水力性质参数化方案研究青藏高原地区土壤冻融过程; 陈渤黎等(2017)使用青藏高原东部若尔盖站的观测数据驱动CLM3.5模式进行模拟研究, 并分析了冻融过程在土壤温度变化和能量分配中的作用。

陆面模式模拟结果的精确性与大气驱动场的质量息息相关(Wang et al, 2011)。李明星等(2011)建立了基于气象台站观测资料的大气驱动场(ObsFC), 并驱动CLM3.5模拟了中国区域土壤湿度变化; 杨晓春等(2015)利用静止气象卫星FY-Y2的反演产品, 结合NCEP/NCAR再分析数据构造了新的大气强迫数据, 并驱动CLM3.0模式模拟土壤湿度; Yin et al(2018)分析了四种大气强迫数据(GSWP3、PGF、CRU-NCEP、WFDEI)对模拟中国区土壤湿度的影响, 结果表明不同数据带来的差异较大。这些研究都表明大气强迫数据对于模式模拟效果有重要影响。此外, 地表覆盖作为陆面模式的地表参数, 决定了地表反射率、蒸发、地面吸收太阳辐射和发射长波辐射的能力(Gibbard et al, 2005; Sr et al, 2011), 其输入资料的质量影响模拟结果的可靠性。Zhao et al(2001)的研究表明中国区域大范围地表覆盖的变化对于气温变化的影响非常明显; Benítez et al(2007)选取了四种地表覆盖资料研究全球碳排放问题, 发现不同地表覆盖数据最大可导致碳排放45 %的变化幅度; Li et al(2017)利用WRF耦合Noah模式分析了不同地表覆盖数据对柏林市热岛效应模拟的影响。

青藏高原地表类型多样, 环境条件恶劣, 土壤温度观测困难, 而且地形高度和独特地理位置使其陆面过程具有特殊性(Huang et al, 2016)。因此, 利用高时空分辨率的大气强迫数据和高空间分辨率的地表植被功能型数据, 驱动模式开展土壤温度模拟研究对于提高模拟精度具有重要意义。本文利用最新研制的高时空分辨率(1 km, 1 h)中国气象局陆面数据同化系统(HRCLDAS-V1.0)大气近地面强迫数据和新的地表植被功能型融合数据(MVEG), 驱动通用陆面模式CLM, 开展了青藏高原10 cm土壤温度模拟研究, 探讨不同的大气强迫和植被功能型数据对模式模拟土壤温度的影响。

2 数据来源与方法介绍 2.1 研究区概况

世界海拔最高、中国面积最大的高原-青藏高原(图 1)被称为“世界屋脊”、“第三极”(陈宇航等, 2016), 平均海拔在4000 m以上, 地势西高东低。地域辽阔, 主体属高原气候, 季风气候性质较为显著, 干湿季分明, 降水多集中于夏季。高原东岸南部地区年平均气温可达到20 ℃, 年降水量为2000 mm; 而西北部地区年平均气温在-6 ℃以下, 年降水量很少, 只有50 mm。总体上, 气候由东南暖湿向西北干冷递变(于海英等, 2009; 解晋等, 2018)。与之相对应, 植被类型也从东南到西北, 逐步由森林向草原、荒漠过渡, 植被种类丰富。其中, 草地分布最广, 主要分布在高原中部; 灌木则主要位于南部地区(于伯华等, 2009)。

2.2 研究数据

观测数据: 2015年逐小时的青藏高原土壤温度观测站点数据由国家气象信息中心资料服务室提供, 共39个观测站点, 每个站点共测量了5, 10, 15, 20, 40, 80, 160和320 cm 8个深度的土壤温度, 其中5~20 cm深度的土壤温度为浅层土壤温度, 40~320 cm深度的土壤温度为深层土壤温度。本研究利用10 cm处的土壤温度观测数据作为验证数据。

大气强迫数据: HRCLDAS-V1.0(High Resolution CMA Land Data Assimilation System Version1.0)是国家气象信息中心最新研制的陆面数据同化系统。CLDAS系统目前发展有三个版本, 分别是CLDAS-V1.0、CLDAS-V2.0和HRCLDAS-V1.0。研究对CLDAS产品质量已经进行了较多的验证与评估, 证实该数据具有较高的质量和应用价值(Qin et al, 2017; Yang et al, 2017; 朱智等, 2014, 2017), 其中HRCLDAS-V1.0产品包括地面气压、2 m气温、2 m比湿、降水、10 m风速和太阳短波辐射等。气温、气压、湿度和风速产品主要通过引入1 km分辨率地形数据, 采用多重网格变分分析制作而成; 短波辐射首先利用风云二号气象卫星可见光信息通道、高分辨率地形、地表反照率等数据来改善地面入射太阳辐射产品质量, 然后应用多重网格变分分析技术将其与Hybrid模型模拟的地面太阳辐射进行融合得到; 降水驱动数据则采用多卫星集成降水产品EMSIP作为背景场, 与全国4万多个站点观测数据进行融合获取(韩帅等, 2018)。与之前的CLDAS系统数据相比, HRCLDAS-V1.0产品具有高时空分辨率特征, 空间分辨率由6 km提高到1 km, 为更高精度土壤温度模拟提供了基础, 本文利用了其最新生成的2015年大气强迫数据开展研究。

地表覆盖和中国植被功能型数据:本文利用的2013年MODIS地表覆盖类型产品MCD12Q1, 是根据MODIS Aqua和Terra卫星数据利用监督决策树分类方法制作的, 空间分辨率为500 m, 共包括5种分类方案, 本研究采用依据IGBP方案分类得到的地表覆盖数据集。该方案中青藏高原地区草地分布最广, 占比为55.55 %, 其次为裸土(37.72 %), 林地为3.55 %, 灌木和农作物较少。此外, 该地区海拔高, 冰雪分布比较广泛, 约为1.10 %, 水体只占0.05 %, 城市和湿地占比相同, 均为0.04 %。

2015年中国土地利用现状遥感监测数据(中国科学院资源环境科学数据中心, 数据下载网址: http://www.resdc.cn)是在2010年数据(以Landsat TM/ETM遥感影像为主要数据源人工目视解译生成)基础上, 基于Landsat8遥感影像, 通过人工目视解译生成。土地利用类型包括耕地、林地、草地、水域、居民地和未利用土地6个一级类型以及25个二级类型, 该数据与植被功能型数据分类体系不一致, 不能用于模式模拟, 本文利用其中1 km分辨率的栅格化数据对新的地表植被功能型融合数据对主要类别进行比较。

中国植被功能型图(1 km)(DOI: 10.3972/westdc.001.2013.db, 寒区旱区科学数据中心)是冉有华等(2009)采用IGBP分类系统(Ran et al, 2012)研制了中国1 km地表覆盖资料(MICLCover), 并在此基础上, 将依Chen et al(2011)制作的1981—2008年的中国大气驱动数据作为气候因子, 根据Bonan(2002)提出的地表覆盖转换为植被功能型的气候规则对MICLCover进行重新划分制作而成的。

中国生态地理分区数据:中国生态地理分区图(中国科学院资源环境科学数据中心, 数据下载网址: http://www.resdc.cn)根据大于10 ℃的积温、大于10 ℃的天数、最冷月平均气温和最暖月平均气温等, 将中国区域划分为16个温度带; 根据年干燥指数、天然植被类型、是否有部分次生盐渍化、是否可旱作以及是否需灌溉等指标将中国区域划分为潮湿、湿润、半湿润、半干旱、干旱和极干旱区6个区域。

2.3 研究方法 2.3.1 新的地表植被功能型融合数据(MVEG)的生成

MCD12Q1产品的数据格式和分类与陆面模式所需的植被功能型数据不同, 因此, 本研究以MODIS地表覆盖产品为基础、中国植被功能型图(1 km)和中国生态地理分区图进行融合, 生成一套分辨率为1 km的新的地表植被功能型融合数据MVEG。首先对MCD12Q1产品进行拼接、投影转换、裁切、格式转换等预处理, 并根据中国生态地理分区图将中国区域分为温带、寒带和热带, 然后如图 2所示, 对林地、灌木、草地等进行重新分类。其中, (1)混交林利用中国植被功能型图替代; (2)利用2013年冬季MODIS的NDVI产品将灌木分为落叶灌木和常绿灌木; (3)极地C₃草位于寒带, 非极地C₃草位于温带, C₄草位于热带(孟现勇等, 2017); (4)在每个格网上, 计算各地表类型占格网的百分比(公式2.1), 最终得到能够应用于陆面模式的高分辨率植被功能型数据。

$ {\left( {{\rm{PC}}{{\rm{T}}_ - }{\rm{PFT}}} \right)_{ij}} = {\rm area }{({\rm{PFT}})_{ij}}/A \times 100\% \;, $

(1)

式中: (PCT_PFT)_ij为第i个网格第j种PFT的百分比; area(PFT)_ij为第i个网格第j种PFT的面积; A为该网格的面积。

2.3.2 CLM模式

通用陆面模式CLM(Community Land Model)模式是目前国际上较为完善可靠的陆面模式之一(Kumar et al, 2011; Hoppe et al, 2014)。模式中用嵌套的次网格来描述空间异质性, 每个网格有3个层次, 第一层为陆地单元, 分为植被、湿地、冰川、湖泊和城市; 第二层是柱单元, 由10层土壤、地表覆盖层和最多5层积雪(积雪层数由雪深决定)表示; 第三层为植被功能类型, 按照植被的生物物理以及化学性质将植被分为不同的植被功能类型, 共定义了16种植被类型(包括裸土)(王平等, 2013)。模式中地表覆盖类型包括湿地、冰川、湖泊、城市和植被五种, 根据这五种地物在网格中所占的比例进行存储, 总和为100 %; 其中植被又被细分为15种植被功能类型, 每个网格中植被功能类型占比最大的四种PFT(plant functional type)进入模式参数的计算, 每个网格中所有PFT的总和为100 %。

CLM中土壤温度的计算(Oleson et al, 2004)采用了土壤热传导方程, 热传导第一定律如式(2)所示:

$ F = - \lambda \nabla T\;, $

(2)

式中: F为单位时间内通过单位横截面的热传导; λ为热导率; ▽T为温度空间梯度; 热传导第一定律一维形式式(3)所示:

$ Fz = - \lambda \frac{{\partial T}}{{\partial z}}\;, $

(3)

式中: z为垂直方向深度(单位: m), 以向下为正, F_z以向上为正。

为了更好的描述不稳定或瞬时状态, 能量原则在模式中采用连续性方程来描述, 如式(4)所示:

$ c\frac{{\partial T}}{{\partial t}} = - \frac{{\partial Fz}}{{\partial z}}\;, $

(4)

式中: c为等体积雪或土壤的热容量; t表示时间。

联合式(3)和式(4), 可得到热传导第二定律形式, 如式(5)所示:

$ c\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial z}}\left[ {\lambda \frac{{\partial T}}{{\partial z}}} \right]\;, $

(5)

CLM模式中计算温度剖面时起初没有考虑相变, 然后重新调整了相变。

2.3.3 土壤温度模拟设计

研究设计了三组实验(表 1), 第一组实验利用CLDAS-V2.0大气近地面强迫数据(6 km, 1 h)驱动CLM作为参照实验, 记为CLDAS/CLM; 第二组实验利用最新的高时空分辨率(1 km, 1 h)的CMA高分辨率陆面数据同化系统(HRCLDAS-V1.0)大气近地面强迫数据驱动CLM, 记为HRCLDAS/CLM; 前两组除大气强迫数据外其他均采用模式自带资料, 在实验二的基础上增加实验三, 用MVEG数据替换模式自带的植被功能型数据进行模拟, 记为HR-MVEG/CLM。三组实验均每6 h输出一次10 cm土壤温度的模拟结果, 即每日00 : 00(世界时, 下同), 06 : 00, 12 : 00和18 : 00 4个时次。

2.3.4 模拟结果评价指标

模拟结果主要利用平均偏差(Bias)、均方根误差(RMSE)和相关系数(R)进行评估(崔园园等, 2018; 朱智等, 2018), 计算公式如下:

$ Bias = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{S_i} - {G_i}} \right)} \;, $

(6)

$ RMSE = \sqrt {\frac{{\sum\limits_{i = 1}^N {{{\left( {{S_i} - {G_i}} \right)}^2}} }}{N}} \;, $

(7)

$ R = \frac{{\sum\limits_{i = 1}^N {\left( {{S_i} - \bar S} \right)} \left( {{G_i} - \bar G} \right)}}{{\sqrt {\sum\limits_{i = 1}^N {{{\left( {{S_i} - \bar S} \right)}^2}} } \sqrt {\sum\limits_{i = 1}^N {{{\left( {{G_i} - \bar G} \right)}^2}} } }}\;, $

(8)

式中: N为样本数量; S_i为模拟结果; G_i为观测站点数据; S和G分别为模拟数据和观测站点数据的平均值。

3 结果与分析 3.1 植被功能型数据比较

CLM模式自带的植被功能型数据集是由Lawrence and Chase(Lawrence et al, 2015)依据MODIS产品制作的(简称LC)。整体来看(表 2、图 3), 青藏高原地区两种植被功能型数据中各种类型的空间分布较为相似, 但是面积占比存在一定的差异, 特别是裸地和草地, 裸地主要分布在高原西北, 并向东南逐渐减少, MVEG中的裸地比LC减少29.2 %; MVEG草地占比明显增多, 且分布更广, 面积占比较LC增加32.3 %; 林地主要分布于高原的东南部, MVEG林地面积有所减少, 但东南边缘区域林地占比却更大; 高原上农作物则主要零星分布在东部部分地区。

由于没有同期的植被功能型数据, 而且地表覆盖或土地利用数据的分类方式不同, 对每种植被功能类型的验证困难, 故研究利用2015年中国科学院资源环境科学数据中心制作的, 由陆地遥感卫星目视解译的中国土地利用现状遥感监测数据, 按裸地、草地、灌木、林地和农作物等大类进行对比。MVEG数据各类别占比较与遥感监测数据较一致(表 2), 相差变化在0.2 % ~3.4 %, 灌木、林地、农作物差异较小, 相对来说, 裸地和草地差异稍大, 而LC数据与遥感监测数据差距较大, 说明MVEG数据精度相对较高, 并有一定合理性, 可以用于模式的模拟。

3.2 逐日土壤温度变化

从2015年逐日10 cm土壤温度的时间变化曲线(图 4)可以看出, 土壤温度呈现出近似正弦曲线式的变化, 且与太阳辐射的变化基本一致, 这与杨梅学等(2000)和陈海存等(2013)的研究结果一致, 这是因为浅层土壤温度的长期变化主要受太阳辐射年变化的影响。一年中的土壤温度最高可达到约23 ℃, 最低温度-5 ℃左右(表 3)。三组实验土壤温度模拟值均能反映随时间的变化规律, 而且反映出土壤温度的短期振荡, 但不同的处理实验略有差异, 其中以HR-MVEG/CLM模拟效果最佳, 趋势变化与观测最接近, 但在第1~60日、320~345日, 模拟的土壤温度略高于观测, 其余时间都与和观测值较吻合; HRCLDAS/CLM在整个研究时间段内与HR-MVEG/CLM变化趋势相似, 但在第80~100日高估了土壤温度峰值, 第130~180日内对峰值和谷值均有所低估, HR-MVEG/CLM则改善了这些时段内对极值的估计。CLDAS/CLM整体上高于观测值, 在第23~28日之间的土壤温度比其他两组最大偏低0.3 ℃, 在这之后明显偏高1 ℃左右。

3.3 模拟结果的相关系数和误差分析

从月尺度的相关系数和误差统计结果(图 5)分析, 三组实验模拟的土壤温度与观测值相关系数均大于0.84, 这说明三组实验在不同月份均有较好的模拟效果, 但相对来说, 8月和9月的相关系数较其他月份偏低。与CLDAS/CLM和HRCLDAS/ CLM相比, HR-MVEG/CLM在大多数月份的相关系数有较为明显的增大, 尤其是4—10月。模拟值与观测的偏差大致介于-1.2~1.6 ℃之间, HRCLDAS/CLM偏差最大减小1.26 ℃, 除在4、5、6月偏差较大外, 其余月份均优于CLDAS/CLM; HR-MVEG/CLM的偏差整体上最低。HR-MVEG/CLM的均方根误差最大不超过2.0 ℃, 比CLDAS/CLM明显减小, 11月甚至减小0.72 ℃, HR-MVEG/CLM还改善了HRCLDAS/CLM在5—6月的表现, 偏差和均方根误差在5月分别比HRCLDAS/CLM减小0.54和0.29 ℃, 6月则分别减小0.48和0.23 ℃。总体上模式模拟的结果在4—10月的表现略逊于其他月份, 与这一时期始末分别是土壤开始消融和冻结的时间, 土壤温度的波动较大, 气温日差相对较大有关。

从季节和年均结果分析(表 4), CLDAS/CLM与HRCLDAS/CLM相关系数大致相近, HR-MVEG/CLM则有所改进, 最大提高0.02(夏季)。HRCL- DAS/CLM和HR-MVEG/CLM弱化了CLDAS/CLM模拟值各季均高估的情况, 但春夏两季模拟值偏低, 且HRCLDAS/CLM春季偏差绝对值高于CLDAS/CLM, 其他时期均有明显减小, 尤其秋季Bias减小了1.10 ℃, HR-MVEG/CLM偏差进一步降低, 但秋季偏差却增大了0.08 ℃(仍低于CLDAS/CLM), 导致年平均偏差也略有升高。各实验都在夏季达到均方根误差的最大值, 且各时期RMSE均为CLDAS/CLM>HRCLDAS/CLM>HR-MVEG/ CLM, HRCLDAS/CLM年平均RMSE减少0.18 ℃, HR-MVEG/CLM又进一步比HRCLDAS/CLM减少了0.09 ℃。

3.4 空间分布模拟分析

从年平均土壤温度空间分布分析(图 6), 模拟值从东南向西北逐渐降低, 与高原地区年均辐射从东南向西北逐渐减少实际相符合。东南地区年均土壤温度在5~20 ℃, 西北地区在0 ℃以下。模式模拟值能够较好地反映出高原地区山脉的走向, 在山脉线上年平均土壤温度明显高于其周围地区, 这是由于高原海拔较高, 其山脉上接受的太阳辐射更多, 因此, 土壤温度较高。由于HRCLDAS/CLM和HR-MVEG/CLM高分辨率优势, 能够模拟出更多的细节分布特征, 在中部地区, HR-MVEG/CLM有明显的白斑缺值, 这是由于高分辨率植被功能型数据中存在湖泊等水体, 模式模拟结果为缺值。HRCLDAS/CLM模拟值略低于CLDAS/CLM, 特别是在东南部的横断山脉, 而西北土壤温度低值区则略高, 低于-10 ℃的区域减少。从HR-MVEG/CLM和HRCLDAS/CLM模拟值差值分布进一步对比可看出[图 6(d)], 自东南向西北, 差值由正逐渐变为负值, 东南地区土壤温度偏差能够达到0~2 ℃, 中部大部分地区差值为-2 ℃左右, 说明植被功能型数据改进对模拟结果的空间分布有较大影响。

4 结论与讨论

在MODIS地表覆盖产品(MCD12Q1)基础上制作了新的地表植被功能型融合数据(下称MVEG), 然后利用最新的高时空分辨率HRCLDAS-V1.0(1 km, 1 h)、CLDAS-V2.0(6 km, 1 h)两种大气近地面强迫数据和MVEG数据设计了三组模拟实验, 驱动CLM模式模拟了2015年青藏高原地区10 cm土壤温度, 并用观测数据进行验证, 结果表明:

(1) 从逐日土壤温度模拟结果时间曲线来看, CLDAS/CLM模拟的土壤温度整体偏高, HRCLDAS/CLM明显降低, HR-MVEG/CLM最接近观测值, 且改善了HRCLDAS/CLM对极值的估计。

(2) 月、季、年平均时间尺度上三组模拟实验结果表明, 三组实验与观测的相关系数均大于0.84。HR-MVEG/CLM误差最小, 年平均偏差为0.07 ℃, 年平均均方根误差仅为1.60 ℃。HRCLDAS/CLM和HR-MVEG/CLM改变了CLDAS/CLM模拟值各季均高估的情况, 大部分时期HRCLDAS/CLM模拟的土壤温度偏差都低于CLDAS/CLM, HR-MVEG/CLM偏差除秋季外均减小, 并且改善了HRCLDAS/CLM在5、6月模拟误差增大的情况。不同季节均方根误差均为HR-MVEG/CLM < HRCLDAS/CLM < CLDAS/CLM。

(3) 从空间分布来看, 青藏高原年平均土壤温度从东南向西北逐渐降低。东南部横断山脉HRCLDAS/CLM土壤温度略低于CLDAS/CLM, 西北区则略高; HR-MVEG/CLM土壤温度在东南部分区域高于HRCLDAS/CLM, 其他区域则略低。说明植被功能型对土壤温度的空间分布有较大影响, 植被茂密区域的浅层土壤温度相对较低。

研究利用了高时空分辨率(1 km, 1 h)的大气驱动数据开展了青藏高原的土壤温度模拟研究, 并且分析了高时空分辨率大气驱动数据和地表植被功能型数据同时改进后对模式模拟的影响。但陆面模式的模拟结果精度不仅与大气强迫和地表覆盖类型有关, 还依赖于模式物理过程参数化方案的改进。高原西部海拔高, 环境恶劣, 设置观测站点难度大, 故观测站点相对集中于东部, 给分析带来一定的不确定性, 后期有待利用更多站点和长时间序列资料进行分析, 并开展多参数对模式模拟的影响研究。