1 引 言
黄河源区位于青藏高原腹地， 是黄河流域主要产流区、 水源涵养区和草原区， 区域内包括高山、 盆地、 草原、 沙地、 湖泊、 湿地、 冰川和冻土， 是重要生态屏障， 对生态环境有深远影响和作用（赵成章和贾亮红， 2009； 唐恬等， 2013）。植被多为高寒草甸、 高寒草原和高寒灌丛， 生态环境极为脆弱， 对全球气候变化较为敏感， 不利于生态系统的稳定（魏卫东等， 2019）。黄河源区为多年冻土区， 近几十年来， 由于气候变暖以及各种因素， 导致大面积冻土退化、 草场退化及逐步沙漠化（黄荣辉等， 2012； 孟宪红等， 2020）， 严重地冲击着黄河源区的水源的涵养功能（曹文炳等， 2006）。多年冻土变化表现为由片状逐渐变为岛状、 斑状分布， 面积缩小且厚度变薄， 或者退化为季节性冻土（金会军等， 2010； 李静等， 2016）， 而季节性冻土也在逐年退化（Qin et al， 2017）。冻土与气候变化相互影响， 一方面， 气候变化影响冻土， 另一方面， 冻土通过冻融过程影响着气候系统（Chen et al， 2014； 徐洪亮等， 2021）。近地层土壤冻融过程是高原地表基本陆面特征之一， 可以作为判断冻土发育、 存在、 变化及反映气候变化的指标（陈瑞等， 2020）。因此， 研究黄河源区土壤冻融过程对地表水热交换的影响， 可为进一步利用观测和模拟方法研究黄河源区域土壤冻融过程的陆面特征及对气候的影响提供参考。
前人已在土壤冻融过程对气候变化的响应方面和对气候系统的影响展开了很多研究。Yang et al（2007）等研究发现土壤冻融状态的变化会影响地气之间的水热交换过程。土壤冻融过程会改变地表与大气之间的辐射、 能量传输特征的改变， 例如感热、 潜热、 动量通量变化及长波辐射等（陈渤黎等， 2017）， 导致地气间水热平衡的变化， 从而影响地表水热过程、 生态环境、 天气及气候系统。土壤冻融过程的存在通过土壤中水分的相变产生储水保温的作用， 影响地表能量交换（Luo et al， 2016； Yang and He， 2019）。陈渤黎（2013）通过对青藏高原有无土壤冻融过程的敏感性实验， 验证了土壤冻融过程会影响地表能量的分配， 加强地气之间的能量传输。冻结（消融）过程会释放（吸收）大量水分相变的能量， 减缓大气对土壤的冷却（加热）作用， 加强地表土壤热通量的传输。Guo et al（2011）对那曲BJ站土壤冻融过程的研究显示， 土壤中的日冻融循环会影响水分相变， 从而显著影响地表热通量， 日冻融循环也会影响地表能量收支平衡。因此随冻土退化， 作为控制地气水热交换的重要过程， 土壤冻融过程必然会受到影响。土壤冻融过程通过影响土壤中温度及湿度的时空变化引起水分相变和影响地表能量平衡， 从而影响表层大气， 因此研究土壤冻融过程对地表水热交换的影响有助于进一步揭示陆面过程对气候的影响机制。由于冻土水热过程的复杂性、 土壤冻融过程的地域空间异质性以及长期连续高密度观测资料的缺乏， 对黄河源区土壤冻融过程的水热性质、 地表能量通量的变化特征以及土壤冻融过程影响地表能量收支行为的认识还不是十分清楚。CLM陆面模式作为目前较为先进和完善的陆面模式， 其包含的生物地球物理和化学过程等， 有利于模拟研究复杂的地表水热交换过程。且CLM模式在青藏高原有较好的适用性， 可以模拟出地表水热参量的变化特征。赖欣等（2014）及袁源等（2019）验证了CLM4.0及CLM4.5模式在青藏高原对土壤湿度的模拟性能， 谢志鹏等（2017）验证了CLM4.5对青藏高原高寒草甸地表能量有较好的模拟性能。CLM5.0是最新版的陆面模式， 相比于CLM4.5， CLM5.0对青藏高原土壤温度及湿度的模拟值更接近观测值（Deng et al， 2020）。因此， 本文首先验证CLM5.0模式在玛曲的适用性， 再结合观测及CLM5.0模拟资料， 对黄河源区玛曲2017年8月至2018年7月土壤冻融过程中土壤温度和湿度、 辐射通量和能量通量的日变化及年变化特征进行分析， 探讨了土壤冻融过程中不同阶段内地表水热交换的变化特征。
2 研究区域概况、 模式介绍和数据来源
2.1 研究区域概况
黄河源区位于青藏高原东北部， 为多年冻土区边缘， 地势西高东低， 海拔为2743~5542 m。气候属于青藏高原亚寒带半干旱半湿润型（文军等， 2011； 胡光印等， 2014）。玛曲位于黄河源区域东南部， 气候寒冷， 冻土类型主要为季节性冻土（郝晓华等， 2009； 郑景云等， 2013）。冬季干冷夏季湿暖， 地表有繁盛的植被覆盖， 下垫面为高寒草甸， 水资源丰富。受大气环流和地形影响， 主要分为冷季和暖季， 冷季时间较长， 年平均气温为1.8 ℃， 年降水量约为630 mm（王少影等， 2012； Li et al， 2021）。
文中涉及的地图是基于中华人民共和国自然资源部地图技术审查中心标准地图服务系统下载的审图号为GS（2020）4621号的中国地图制作， 底图无修改。
2.2 CLM5.0模式介绍
CLM5.0是由公用地球系统模式（Community Earth System Model， CESM）陆面模型工作组（Land Model Working Group， LMWG）开发并由美国国家大气研究中心（National Center for Atmospheric Research， NCAR）维护的一系列全球陆面模型中的最新版本。CLM5.0在CLM4.5版本的基础上对多数主要部分进行更新， 特别是对土壤和植物水文、 雪密度、 河流建模、 作物建模、 碳和氮循环及耦合进行了显著更改（Lawrence et al， 2019）。Pelletier et al（2016）利用空间土壤厚度数据产品使土壤厚度由原来的空间均匀改为随空间变化。CLM5.0用零通量边界条件及饱和区和非饱和区建模。默认模型土壤层分辨率增加， 尤其是在顶部3 m内， 更明确地表示永久冻土带内的活动层厚度。引入了理查德方程的自适应时间步长解， 提高了数值土壤水的精度和稳定性（Jackson et al， 1996）。本研究主要使用陆面模式中的“I”组件进行离线单点模拟， spin-up 3个月后进行一年的模拟。大气强迫数据主要使用站点的涡动数据中的气温、 湿度、 气压、 风速、 降水、 入射太阳辐射和向下长波辐射。地表数据主要由三层次网格（陆地单元、 各单元中要素的变化及植被功能类型）嵌套而成， 用于表征地表异质性、 土壤垂直结构及组成等， 包括土地覆盖、 植被类型、 土壤分层及组成成分等信息。其中土壤质地包括沙土、 黏土和泥沙含量（单位： %）对模拟土壤温度及湿度影响较大， 因此用中国土壤特性数据集中的土壤质地替换。地表数据中土壤质地的土壤层分为10层（0.01~1.36 m）， 更深层土壤质地数据在模式运行时使用第十层数据， 中国土壤特性数据集中土壤各层厚度与地表数据中土壤层厚度不相同， 用0~0.045 m替换0.01 m及0.04 m、 0.045~0.091 m替换0.09 m、 0.091~0.166 m替换0.16 m、 0.166~0.289 m替换0.26 m、 0.289~0.493 m替换0.4 m、 0.493~0.829 m替换0.58 m和0.8 m、 0.829~1.383 m替换1.06 m和1.36 m的土壤质地数据。
2.3 数据来源
所用的观测资料来自于“黄河源区玛曲陆面过程与遥感试验平台”， 该观测平台位于甘肃省甘南藏族自治州玛曲县境内（图1 ）。数据主要包括以下部分： （1）涡动相关系统观测系统， 使用EC150开路式分析仪（Campbell scientific， USA）、 CSAT3三维超声风速仪（Campbell scientific， USA）和HMP115A相对湿度\温度传感器（Campbell scientific， USA）。用涡动相关数据处理软件对通量进行了误差修正。其中感热通量、 潜热通量、 净辐射通量， 测量高度为2.85 m。输入驱动CLM5.0大气强迫数据的长、 短波辐射、 气压、 相对湿度和风速， 辐射使用四分量净辐射计（CNR-1， Kipp and Zonen， Delft， Netherlands）， 测量高度为1.5 m， 采样时间为30 min， 在2017年8月1日至2018年7月31日， 辐射通量缺测率为6.18%， 感热、 潜热、 净辐射通量缺测率为10.02%， 2017年11月17日至12月17日连续缺测， 其他时段部分时刻缺测。（2）土壤温度/湿度观测系统观测资料。使用的是ECH2O 5TM传感器（METER Group， Inc.USA）测量土壤温度及体积含水量， 精确度3%。土壤湿度数据主要测量的是土壤中的液态水含量， 故土壤水分含量指液态水含量。采样时间为15 min， 考虑到数据的连续性， 主要选取0.05 m和0.1 m层数据。（3）土壤热通量和积雪深度观测资料。土壤热通量使用热通量板（HPF01， Wohlwend Engineering， Sennwald， Switzerland）数据， 精确度为-15%~5%， 时间分辨率为15 min， 包括5 cm、 10 cm、 20 cm和40 cm深度数据。积雪深度为日资料。考虑到研究目的及数据的连续性， 统一选取2017年8月至2018年7月的数据。由于缺乏降水资料， 且入射短波辐射、 向下长波辐射、 气压、 相对湿度和风速在2017年11月17日至12月17日连续缺测， 在制作驱动CLM5.0模式的大气强迫资料时， 检验数据连续性后， 使用国家青藏高原科学数据中心的中国区域地面气象要素驱动数据集（1979 -2018）（Yang and He， 2019）中的降水数据代替2017年8月至2018年7月的观测降水数据， 用入射短波辐射、 向下长波辐射、 气压、 相对湿度和风速替换缺测时间段数据， 空间分辨率为0.1°×0.1°， 时间分辨率为3 h。使用中国土壤特性数据集（Shangguan et al， 2013）中的土壤质地（沙土、 黏土和泥沙）数据替换地表数据中的土壤质地。
3 数据处理及分析方法
3.1 土壤热通量的计算
由于地表热通量难以直接测量， 根据Tanaka et al（2003）和Evett et al（2012）使用0.05 m的土壤热通量以及结合0~0.05 m层的热储量对地表热通量进行校正， 其中0~0.05 m的土壤温度由地表及0.05 m土壤温度的均值表示。且假定土壤中的水分变化仅由水的相变导致（陈星等， 2014； 武月月等， 2020）。
地表温度利用地表长波辐射进行计算：
式中： R UL为向上长波辐射； R DL为向下长波辐射； 地表比辐射率 为0.96； 斯蒂芬玻尔兹曼常数σ为5.67×10-8 W·m-3·K-4。
地表热传导方程表示为蓄热的垂直积分：
最右侧一项为冻结潜热项， 冻结-潜热常数Lf 为3.337×105 J·kg-1； 为土壤含冰量（单位： m3·m-3）； 为冰的质量密度（单位： kg·m-3）。在土壤冻融过程中， 水分相变， 热量发生变化， 所以当土壤温度小于0 ℃时， 加入未冻水含量变化项， 含冰量由相邻时刻的含水量之差表示（李光伟等， 2019）， 土壤热容量为：
式中： 为未冻水含量（单位： m3·m-3）， 和 分别为干土、 液态水和冰的质量密度（单位： kg·m-3）， 、 和 分别为干土、 液态水和冰的比热容（单位： J·kg-1·K-1）。干土热容量 约为0.09×106 J·m-3·K-1（Lawrence et al， 2019）， 液态水热容量 约为4.02×106 J·m-3·K-1， 冰的热容量 约为1.89×106 J·m-3·K-1。
参考深度为0.05 m， 式（2） 的离散形式可以表示为：
3.2 土壤冻融阶段的划分及典型晴天的定义
土壤冻融阶段定义为四个阶段， 分别为冻结过程、 完全冻结、 消融过程和完全消融阶段（Guo et al， 2011）。各冻融阶段的划分主要依据0.05 m土壤温度的日变化， 在不考虑盐度对冻融过程的影响的情况下， 当0.05 m土壤温度日最高温度小于0 ℃时为完全冻结阶段， 当土壤日最低温度大于0 ℃时为完全消融阶段， 当日最高温度大于0 ℃， 且最低温度小于0 ℃， 完全冻结前为冻结过程， 之后为消融过程。为了回避随机天气过程的可能影响， 当连续三天满足冻融状态判定条件时， 取第一天作为下一冻融状态的开始日。对2017年8月1日至2018年7月31日冻融过程阶段具体划分情况如表1 所示， 冻结过程阶段为2017年11月26日至12月1日、 完全冻结阶段为2017年12月2日至2018年2月27日、 消融过程阶段为2018年2月28日至3月3日、 完全消融阶段为2017年8月1日至11月23日以及2018年3月4日至7月31日。
表1 土壤冻融过程阶段划分及典型晴天日Table 1 Date of soil freezing and thawing process stage process and typical clear day |
| 土壤状态 | 起止时间 | 日数/天* | 典型晴天日 | |
|---|---|---|---|---|
| 开始时间 | 结束时间 | |||
| 完全消融 | 2017-08-01 | 2017-11-23 | 115 | 2017-09-11、 2017-10-06、 2017-11-01、 2017-11-13 |
| 冻结过程 | 2017-11-24 | 2017-12-02 | 9 | - |
| 完全冻结 | 2017-12-03 | 2018-02-26 | 86 | 2017-12-10、 2017-12-12、 2017-12-20、 2017-12-24、 2017-12-27、 2018-01-5、 2018-01-12、 2018-01-22~2018-01-24 |
| 消融过程 | 2018-02-27 | 2018-03-15 | 17 | 2018-03-05 |
| 完全消融 | 2018-03-16 | 2018-07-31 | 137 | 2018-04-08、 2018-04-16、 2018-05-12、 2018-05-13 |
*日数： 仅对2017年8月1日至2018年7月31日中各土壤冻融过程阶段日数进行统计 |
平均日变化中包含了不同天气过程的影响， 因此用不同冻融阶段典型晴天的日变化来进行对比分析。当日降水量为0 mm， 且当天太阳辐射呈现光滑的倒“U”型， 确定为典型晴天日。完全冻结阶段的典型晴天日为2017年12月10日、 12日、 20日、 24日、 27日， 2018年1月5日、 12日、 22 -24日， 共10天。消融过程阶段典型晴天日为2018年3月5日， 完全消融阶段典型晴天日为2017年9月11日、 10月6日、 11月1日、 13日， 2018年4月8日、 16日、 2018年5月12日、 13日， 共8天。由于太阳辐射资料在冻结过程阶段缺测， 故缺少此阶段典型晴天日。玛曲土壤冻融过程日数共112天， 约占全年的1/3， 其中完全冻结阶段占冻融期的3/4， 冻结过程和消融过程阶段所占日数较少。
4 土壤冻融过程中地表水热交换的观测及模拟研究
4.1 CLM5.0模拟结果验证
由于高原特殊的地理位置， 仪器维护难度大， 仪器在运行过程中易受人为因素等不可控因素的影响导致获取的观测数据存在部分甚至连续多日的缺测， 因此对相应气候要素进行模拟及验证后用模拟结果补充缺测的观测数据。土壤冻融过程的实质是由于土壤温度变化所引起的土壤水分发生相变的过程， 因此选择土壤温度和土壤湿度来验证模式在玛曲地区的适用性（陈瑞等， 2020）。图2 为2017年8月至2018年7月0.1 m土壤温度和土壤湿度模拟和实测数据的对比图。如图2 （a）所示， 模式可以模拟出土壤温度的动态变化， 尤其在非冻融期， 模式不仅显示了土壤温度的变化趋势， 且模拟土壤温度在数值上与观测值更相近。如图2 （b）所示， 模式模拟的土壤湿度与观测土壤湿度变化趋势一致， 当土壤温度大于0 ℃， 模式可以模拟出土壤湿度的动态变化， 冻融期土壤湿度更接近观测值。土壤温度及湿度模拟值偏小。综上所述， CLM5.0模式在玛曲站有一定的适用性， 可以再现玛曲土壤温度和土壤湿度等物理量的变化特征。
定量分析也显示模式在玛曲的适用性。表2 为玛曲入射短波辐射（DS）、 反射短波辐射（US）、 向下长波辐射（DL）、 向上长波辐射（UL）、 净辐射通量（Rn）、 感热通量（SH）、 潜热通量（LE）、 波文比（β）和土壤热通量（G0）在不同冻融阶段的阶段平均值。根据表2 ， .1 m土壤温度模拟和观测数据的相关系数为0.991， 通过95%的显著性检验。土壤温度的模拟值较观测值的均方根误差为1.78 ℃。0.1 m土壤湿度模拟和观测数据的相关系数为0.962， 土壤湿度的模拟值较观测值的均方根误差为0.07 m3·m-3。土壤温度、 湿度的平均偏差为-1.31 ℃和-0.06 m3·m-3， 模拟比观测数据偏小。多种因素会造成模拟值与观测值差异， 在非冻融期主要受输入模式的大气强迫资料中气温和降水与实际气温和降水的差异的影响， 在冻融期可能受冻融期土壤导热率与未冻水参数化方案的影响， 以及模式的植被类型、 土壤属性和土地利用分类等数据与玛曲站点实际情况有差异等因素的影响（李时越等， 2018； 李光伟等， 2019； 袁源等， 2019； 陈瑞等， 2020； Deng et al， 2021）。
表2 土壤温度、 湿度、 辐射通量及能量通量日均值误差统计Table 2 Error statistics of daily mean soil temperature， soil moisture， radiant flux and energy flux for Maqu site between the simulation data and observations |
| 气象要素 | R | RMSE | Bias |
|---|---|---|---|
| 土壤温度 | 0.991* | 1.78 ℃ | -1.31 ℃ |
| 土壤湿度 | 0.962* | 0.07 m3·m-3 | -0.06 m3·m-3 |
| 入射短波辐射 | 0.966* | 19.3 W·m-2 | 1.2 W·m-2 |
| 反射短波辐射 | 0.737* | 21.4 W·m-2 | -10.6 W·m-2 |
| 向下长波辐射 | 0.998* | 4.5 W·m-2 | 0.5 W·m-2 |
| 向上长波辐射 | 0.979* | 8.6 W·m-2 | 3.0 W·m-2 |
| 感热通量 | 0.553* | 29.1 W·m-2 | 16.2 W·m-2 |
| 潜热通量 | 0.901* | 11.1 W·m-2 | -0.8 W·m-2 |
| 净辐射通量 | 0.945* | 21.6 W·m-2 | 9.0 W·m-2 |
*表示通过95%的显著性水平检验 |
图3 2017年8月至2018年7月玛曲CLM5.0模式和观测数据的入射短波辐射（a）、 反射短波辐射（b）、 向下长波辐射（c）、 向上长波辐射（d）、 感热通量（e）、 潜热通量（f）和净辐射通量（g）日平均变化Fig.3 Changes of daily mean downward shortwave radiation （a）， upward shortwave radiation （b）， downward longwave radiation （c）， upward longwave radiation （d）， sensible heat flux （e）， latent heat flux （f） and net radiation flux （g） from CLM5.0 model and observations at Maqu from August 2017 to July 2018 |
4.2 土壤温度和湿度日变化特征
为了探究土壤冻融过程对地气系统中水热变化的特征， 将对包含一个完整冻融期的土壤温度、 湿度、 地表能量通量的日变化及日均值年变化特征进行分析。冻融过程各个阶段有明显的日变化特征， 选取四个冻融阶段土壤温度及湿度阶段平均日变化以及各冻融阶段的典型晴天日变化进行分析。由图4 可知， 浅层各土壤冻融阶段土壤温度有明显的日变化， 日最高温度出现在17:00（北京时， 下同）左右， 日最低温度出现在09:00左右。冻结过程阶段及消融过程阶段中， 日最高土壤温度大于0 ℃， 日最低土壤温度小于0 ℃， 此时土壤湿度随温度上升（下降）而增加（减少）， 土壤湿度日变化明显， 当温度降至0 ℃以下时土壤中存在液态水冻结， 当温度随之上升时土壤中冰融化， 土壤中存在明显的日冻融循环。水分相变中吸收或释放的能量又会影响土壤温度的变化， 使其日变幅减小。完全冻结阶段土壤中水分总量较少以及完全消融阶段存在大量水分但不涉及水分相变， 相较于冻结过程和消融过程阶段， 土壤温度日变幅较大， 土壤湿度日变化不明显。消融过程阶段土壤温度和土壤湿度的变幅比冻结过程阶段大。与阶段平均日变化相比， 典型晴天日变化趋势与阶段平均日变化相同， 但典型晴天日变幅更大， 且受当天天气的影响， 曲线有明显的波折。
4.3 辐射通量日变化特征
到达地表的太阳辐射相当于能量的收入部分， 对地表辐射收支、 能量交换以及天气气候形成有决定性的意义。入射短波辐射、 反射短波辐射及向上长波辐射在四个冻融阶段有明显相同的日变化特征， 随着太阳辐射呈单峰型变化， 在12:00 -14:00出现最大值。入射及反射短波辐射在夜间随太阳高度角为0°， 短波辐射值为0 W·m-2。长波辐射易受大气透明度、 气溶胶、 水汽、 气温和二氧化碳等的影响。向下长波辐射日变幅相对较小， 没有明显的单峰型日变化， 最大值出现在16:00 -18:00， 变化趋势与短波辐射等有所差异。辐射四分量在典型晴天日变化与各土壤冻融阶段的平均日变化特征相同。除向下长波辐射通量外， 典型晴天日的辐射通量变幅大于各阶段的平均日变化， 且波动更加剧烈。当土壤处于土壤冻融过程时， 土壤中存在冰会加大地表反照率。反照率在冻融过程阶段大于完全消融阶段（表3 ）。通过对玛曲2017年8月至2018年7月各阶段积雪日数统计（表3 ）， 冻结过程阶段没有积雪， 积雪主要分布2018年3月18日至5月7日， 即完全消融阶段。完全融化阶段的积雪日数更多， 地表有积雪时会影响地表反照率， 从而使反射短波辐射增加。当剔除积雪日后， 反照率在四个阶段分别为0.25、 0.29、 0.30和0.22， 相对于表3 中的没有剔除积雪日的反照率， 完全消融阶段反照率明显减小。而其他三个阶段反照率变化较小， 冻结过程和完全冻结阶段的入射短波辐射差距较小， 反射短波辐射在完全冻结阶段明显大于冻结过程阶段， 可能的原因是完全冻结过程阶段的土壤冰含量更多， 反射短波辐射增加地表反射率增加。冻融过程影响地表辐射收支。
表3 2017年8月至2018年7月玛曲各冻融过程阶段的辐射通量、 能量通量、 反照率、 波文比平均值及降雪日数统计Table 3 Average radiation flux， energy flux， albedo and Bowen ratio and statistics of snowfall days in different freeze-thaw process stages at site Maqu from August 2017 to July 2018 |
| 冻融阶段 | DS /（W·m-2） | US /（W·m-2） | Albedo | DL /（W·m-2） | UL /（W·m-2） | Rn /（W·m-2） | SH /（W·m-2） | LE /（W·m-2） | β | G 0 /（W·m-2） | 积雪日数 /天 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 冻结过程 | 141.3 | 34.4 | 0.25 | 203.4 | 293.9 | 16.4 | 39.3 | 12.0 | 3.3 | -9.1 | 0 |
| 完全冻结 | 146.8 | 43.7 | 0.29 | 197.8 | 285.2 | 15.9 | 30.1 | 5.0 | 6.1 | -5.6 | 6 |
| 消融过程 | 201.2 | 60.7 | 0.30 | 225.8 | 315.4 | 50.1 | 38.8 | 6.6 | 5.9 | 3.4 | 1 |
| 完全消融 | 216.9 | 52.0 | 0.24 | 295.9 | 360.9 | 117.9 | 30.4 | 89.5 | 0.3 | 1.2 | 18 |
除土壤热通量及降雪日数外， 冻结过程阶段缺测的观测数据由模拟数据代替 |
图5 各冻融阶段玛曲辐射能量通量阶段平均 （a~d）及典型晴天（e~h）日变化其中冻结过程阶段缺失的观测数据用模拟结果代替 Fig.5 Diurnal variations of radiant energy flux at Maqu averaged （a~d） and typical clear day （e~h） in different freeze-thaw process stages.The missing observation data in the freezing process stage is replaced by the simulation results |
4.4 能量通量日变化特征
通过四个冻融阶段能量通量平均日变化和典型晴天日变化， 探讨各个阶段能量通量特征及土壤冻融过程对能量通量的影响。如图6 ， 感热、 潜热、 净辐射和土壤热通量有相同的日变化， 能量通量阶段平均和典型晴天的日变化呈现单峰形， 最大值出现在14:00， 能量通量在典型晴天峰值更大。感热通量在消融过程和完全冻结过程阶段较大， 在完全消融和冻结过程较小［图6 （a）］。潜热通量在完全消融阶段变化明显， 土壤中水分充足， 水分频繁相变使潜热通量变化幅度明显大于其他阶段［图6 （b）］。受土壤中较小的土壤湿度限制， 在完全冻结阶段潜热通量及变幅最小。白天净辐射通量大于0 W·m-2， 受太阳辐射的影响净辐射通量在完全消融阶段最大， 在完全冻结阶段的最小， 冻结过程阶段大于消融过程阶段［图6 （c）］。土壤热通量的日变幅小于其他能量通量。白天土壤热通量大于夜间， 白天土壤热通量大于0 W·m-2获得热量， 夜间土壤释放热量［图6 （d）］。完全消融阶段土壤热通量最大， 冻结过程阶段土壤热通量最小。
图6 各土壤冻融阶段玛曲感热通量、 潜热通量、 净辐射通量和土壤热通量阶段平均（a ~d）及典型晴天（e ~h）日变化除土壤热通量外， 在冻结过程阶段缺失的观测数据由模拟结果代替 Fig.6 Diurnal variations of sensible heat flux， latent heat flux， net radiation flux and soil surface heat flux at Maqu averaged（a ~d） and typical clear day （e ~h） in different freeze-thaw process stages.The missing observation data in the freezing process stage is replaced by the simulation results， except soil surface heat flux |
4.5 土壤温度、 湿度及能量通量的年变化特征
土壤湿度和温度的变化反映了土壤水分和能量状况。土壤的热容量和含水量均大于大气， 所以土壤能量和湿度的变化是陆-气系统能量和水的循环过程中的一个重要方面（Yang et al， 2003）。如图7 （a）所示， 2017年8月到2018年7月， 土壤温度和湿度变化整体呈现“U”型。从2017年11月至2018年3月， 随着土壤温度逐渐低于0 ℃， 土壤经历了从开始冻结到消融的一个完整的冻融期。完全消融阶段， 地表植被丰富， 土壤中水分充足， 受充足的入射太阳辐射以及较低的地表反照率影响， 地表吸收更多的辐射能量， 地表土壤水分蒸发以及通过植被的蒸腾作用进入大气。地表蒸发依赖于太阳辐射的变化（Guo et al， 2011）。当处于冻结过程及消融过程阶段时， 土壤中的液态水昼融夜冻， 白天土壤湿度较高， 蒸发潜热较大， 夜间土壤湿度较低， 蒸发潜热减小。完全冻结阶段土壤中水分最少， 土壤中水分蒸发较少， 土壤中的水分受太阳辐射的影响存在相变。
图7 2017年8月至2018年7月0.05 m层土壤温度和湿度变化（a）及能量通量变化（b）F为冻结过程阶段, CF为完全冻结阶段, T为冻结过程阶段, CT为完全消融阶段; 冻结过程阶段除土壤热通量外用模拟结果代替缺失的观测数据 Fig.7 Variation of soil temperature and soil moisture at 0.05 m （a） and variation of energy flux （b） at site Maqu from August 2017 to July 2018.F refers to freezing stage， CF refers to completely frozen stage， T refers to thawing stage， CT refers to completely thawed stage.The missing observation data of energy flux in the freezing process stage is replaced by the simulation results， except soil surface heat flux |
土壤中水分相变会影响近地表的能量通量分配， 日冻融循环通过影响水分相变影响热量交换强度（李述训等， 2002）。净辐射为地表辐射能量收支的差额， 净辐射与地表热通量的差值是直接驱动地气相互作用的可利用能量（Ma et al， 2005； 王少影等， 2012）， 相较于净辐射通量， 土壤热通量数值上较小， 地表可利用能量变化特征主要受净辐射通量的影响。净辐射受太阳辐射的影响在夏季最大， 冬季最小。2017年8月至2018年7月呈现潜热通量大致呈现“U”型变化， 而感热通量与之变化相反。潜热通量变化与0.05 m土壤湿度的变化趋势相似， 在土壤湿度维持较大值， 此时下垫面植被较多， 水分通过蒸散发带走热量， 此阶段潜热通量较大且有明显的日变化， 潜热通量与净辐射通量日变化趋势一致， 地表可利用能量主要用于潜热通量。当土壤中水分较少时， 抑制水分相变， 潜热通量较小， 地表能量易受土壤温度变化的影响， 尤其是完全冻结阶段地表可利用能量主要用于感热通量。在冻结过程阶段， 土壤水分减少、 净辐射通量、 潜热通量和感热通量日均值均减小， 潜热通量随土壤湿度减小而减小， 与净辐射通量变化一致。相比于冻结过程阶段， 消融过程阶段潜热通量随土壤湿度增加而增加， 但此阶段净辐射通量增加趋势更大， 通过影响地表空气温度影响感热通量， 感热通量变幅较大。以上说明， 日冻融循环引起的土壤湿度的变化会影响潜热通量的变化， 影响地表能量分配， 这与Guo et al（2011）的结果是一致的。
地表热通量变化幅度最小， 地表热通量正负仅代表能量传输方向， 当地表土壤热通量为正值时， 土壤表面得到热量并向下传输， 当地表热通量为负值时， 土壤表面失去热量， 能量向上传输。土壤热通量的变幅最小， 在冻结过程阶段和完全冻结阶段土壤热通量小于0 W·m-2， 在完全消融阶段和消融过程阶段较小。从表3 中的土壤热通量在冻结过程阶段和完全冻结阶段平均值为-9.1 W·m-2和-5.6 W·m-2均小于0， 在消融过程阶段和完全消融阶段平均值为3.4 W·m-2和1.2 W·m-2。冻融过程阶段G 0（绝对值）大于完全消融阶段， 冻融过程阶段会加强地表热量交换。
4.6 土壤冻融过程中地表能量收支平衡的变化特征
地表能量闭合率为潜热和感热通量之和与净辐射与土壤热通量之差的比值。基于能量平衡基本原理， 近地面可利用能量应该与可分配能量相等， 能量闭合有利于模式及许多应用性研究的发展（王介民等， 2009）。然而在Foken（2008）研究中， 近地表能量是不闭合的， 不闭合率为10%~30%， 主要与仪器精度、 土壤热通量计算及时空尺度有关（Stannard et al， 1994）。如图8 （a）所示， 除冻结过程阶段能量通量均由模拟结果代替， 其余阶段闭合率在1左右， 本套观测资料质量通过评估， 观测资料质量较好。冻结过程阶段闭合率较大可能受以下因素影响： 模式在模拟过程中假设能量闭合率为1， 感热通量和净辐射通量模拟结果均大于观测结果， 且感热模拟的平均偏差更大（表2 ）。根据表2 计算各阶段闭合率， 冻结过程阶段、 完全冻结阶段、 消融过程阶段和完全消融阶段平均闭合率分别为1.44、 1.56、 0.99和0.81。一般情况下当处于冻结过程阶段时， 土壤水结冰释放热量， 从而高估土壤热通量进而高估能量闭合率。当处于消融过程阶段时， 土壤中的冰吸收热量转化为液态水， 低估了土壤热通量进而低估能量闭合率（Guo et al， 2011）。夜间土壤温度较低， 土壤湿度较小， 感热和潜热约为0 W·m-2， 地表失去热量， 净辐射通量小于0 W·m-2， 白天闭合率较大， 夜间闭合率较小［图8 （b）］（李泉等， 2008）。因此， 当仅考虑白天能量闭合率时， 完全冻结阶段的能量闭合率更接近1， 能量闭合更好。白天土壤温度上升， 土壤中的冰吸收热量转化为液态水， 尤其是12:00 -18:00冻结过程和消融过程阶段的土壤热通量绝对值较小进而低估能量闭合率［图6 （d）］， 图8 （b）中除了冻结过程阶段外， 夜间能量闭合率约为0， 而冻结过程夜间能量闭合率远小于0， 可能的原因是模式高估冻结过程阶段的感热通量， 从而低估了能量闭合率。
图8 地表能量闭合率日均值年变化（a）及日变化（b）F为冻结过程阶段， CF为完全冻结阶段， T为冻结过程阶段， CT为完全消融阶段； 冻结过程阶段闭合率由模拟的感热通量、 潜热通量、 净辐射通量及计算的土壤热通量计算 Fig.8 Annual variation of daily average CR （a）， diurnal variations of closure ratio （CR） averaged in different freeze/thaw stages （b）.F refers to freezing stage， CF refers to completely frozen stage， T refers to thawing stage， CT refers to completely thawed stage.The closure rate of freezing process is calculated from the simulated sensible heat flux， latent heat flux， net radiation flux and calculated soil heat flux |
5 结论
本文使用CLM5.0模式对玛曲站进行了单点模拟实验。通过验证， 该模式可以较好地模拟玛曲土壤温湿度及能量通量。利用观测资料及模拟结果统计分析了黄河源区玛曲站2017年8月至2018年7月的土壤温湿度、 辐射通量及能量通量在土壤冻融过程中的时间演变特征， 对玛曲站土壤冻融过程中地表水热交换的变化有了一定的认识。以下为主要结论：
（1） 通过对土壤温度及湿度模拟， 验证了CLM5.0模式在玛曲的适用性。模式可以模拟土壤温度和土壤湿度的变化趋势， 土壤温度在非冻融期模拟结果更好， 土壤湿度在冻融期数值上更接近观测值。模式可以模拟辐射通量及能量通量的变化趋势及状态。模式低估土壤温度、 湿度、 反射短波辐射及潜热通量。模式高估入射短波辐射、 向下长波辐射、 向上长波辐射、 净辐射通量及感热通量。
（2） 土壤中的水分受温度的影响发生相变， 对地表潜热通量的影响较大， 从而影响地气之间的能量分配。0.05 m土壤温度和湿度有明显的日变化， 在冻结和消融过程阶段中有日冻融循环特征。完全消融阶段， 土壤水分充足， 水分吸收热量蒸发， 能量通量以潜热通量为主。在完全冻结阶段， 受土壤湿度限制潜热通量最小， 地表能量通量以感热通量为主。在冻结过程， 净辐射通量较少， 地表能量主要还是受潜热影响， 而消融过程阶段， 受较大的净辐射通量的影响， 地表能量以感热通量为主。
（3） 土壤中存在日冻融循环， 水分相变导致潜热通量随土壤湿度变化， 当浅层土壤湿度增加时表层蒸发潜热增加， 土壤中水冻结减少导致蒸发潜热通量减小， 日冻融循环会通过影响水相变导致频繁的潜热通量释放影响地表热通量。冻融过程土壤热通量绝对值大于非冻融过程阶段， 即日冻融循环时表层能量交换更强。各冻融阶段典型晴天变化趋势与阶段平均值一致， 但幅度要大于阶段平均值。
（4） 冻结过程、 完全冻结、 消融过程和完全消融阶段平均能量闭合率分别为1.44、 1.56、 0.99和0.81， 消融过程和完全消融过程能量闭合率较好。日冻融循环会影响地表能量平衡， 由于土壤中的水相变释放（吸收）热量， 导致高估（低估）土壤热通量从而高估（低估）能量闭合率。
由于冻融过程易受冻土类型、 径流、 积雪、 植被类型和土壤有机质等局地因素的影响， 对于黄河源区复杂的下垫面， 在未来的工作中一方面计划结合更多的具有代表性的站点观测资料探究黄河源区土壤冻融过程及其对地表水热交换的影响。这也有利于进一步改进土壤冻融过程参数化方案， 提升模式在区域尺度上的适用性。