下边界条件对多年冻土温度场变化数值模拟的影响

doi:10.7522/j.issn.1000-0240.2021.0048

下边界条件对多年冻土温度场变化数值模拟的影响

孙哲^,¹^,², 赵林^,³, 胡国杰¹, 乔永平¹, 杜二计¹, 邹德富¹, 谢昌卫¹

1.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室/藏北高原冰冻圈特殊环境与灾害国家野外科学观测研究站，甘肃兰州 730000

2.中国科学院大学，北京 100049

3.南京信息工程大学地理科学学院，江苏南京 210044

Influence of lower boundary conditions on the numerical simulation of permafrost temperature field changes

SUN Zhe^,¹^,², ZHAO Lin^,³, HU Guojie¹, QIAO Yongping¹, DU Erji¹, ZOU Defu¹, XIE Changwei¹

1.State Key Laboratory of Cryospheric Science / Cryosphere Research Station on the Qinghai-Tibet Plateau，Northwest Institute of Eco-Environment and Resources，Chinese Academy of Sciences，Lanzhou 730000，China

2.University of Chinese Academy of Sciences，Beijing 100049，China

3.School of Geographical Sciences，Nanjing University of Information Science and Technology，Nanjing 210044，China

通讯作者: 赵林，教授，主要从事多年冻土与环境研究. E-mail： lzhao@nuist.edu.cn

编委: 武俊杰

收稿日期: 2021-01-05 修回日期: 2021-03-03

基金资助:

第二次青藏高原综合科学考察研究项目.  2019QZKK0201
国家自然科学基金项目.  41931180.  42071094
中国科学院战略性先导科技专项.  XDA2002010202

Received: 2021-01-05 Revised: 2021-03-03

作者简介 About authors

孙哲，博士研究生，主要从事多年冻土模拟研究.E-mail：sunzhe@lzb.ac.cn , E-mail：sunzhe@lzb.ac.cn

摘要

在气候变暖背景下，北半球多年冻土呈现不同程度的退化趋势，冻土升温、活动层增厚、地下冰消融改变了区域工程地质条件、地形地貌，不仅对寒区环境和工程稳定性造成潜在的威胁，还影响着这些地区的气候、水文和生态过程。因此，准确评估和预估多年冻土热状况的变化具有重要科学和实践意义。现有用于模拟多年冻土热状况的各类模式重点考虑了近地表温度场变化对多年冻土的影响，主要集中于对气温和浅表层物理过程和参数化方案等改进和优化，而对于下边界条件设置对多年冻土热状况模拟的影响少有讨论。基于一维热传导冻土模型，以五道梁地区的多年冻土为研究对象，通过设置不同的下边界方案进行模拟实验，定量评估百年尺度气候变化下不同下边界条件对多年冻土温度场变化数值模拟的影响。结果表明：近地表层（<3 m）的温度场完全由年际气候变化决定，浅层（3~15 m）及中层（15~30 m）的多年冻土温度场受下边界条件的影响逐渐显著，深层（>30 m）的地温对百年尺度气候变化的响应不仅与气候变化的幅度有关，还与多年冻土相变热的多少有直接的关系。下边界条件不恰当的设置方式会对大尺度的气候变化下多年冻土消融程度的计算造成较大的影响，进而可能对深层地温乃至多年冻土区面积变化造成严重的误判。因此，开展百年尺度多年冻土温度场变化模拟时，应采用深层或多年冻土底板以下融土层的稳定地热流作为下边界条件。

关键词： 多年冻土 ; 地温场 ; 数值模拟 ; 下边界条件

Abstract

Permafrost in the Northern Hemisphere has degraded in different degree under the climate warming. Permafrost degradation may have changed the regional engineering geological conditions， topography and geomorphology， which may have a great impact on regional climate， hydrological and ecological processes， and engineering infrastructures. It is essential to assess and predict the permafrost thermal state accurately by models. Most present models used to simulate permafrost thermal state focus on the effect of ground surface temperature on permafrost changes， so their research direction is the improvement and optimization of the parameterization for air temperature and the physical process in shallow soil layer. However， the influence of the temperature field in the lower part of soil layer， which is strongly affected by historical climate and geothermal energy， on the permafrost simulation is rarely discussed. A comprehensive observation station has been established in Wudaoliang region on the Qinghai-Tibet Plateau by the Cryosphere Research Station on the Qinghai-Tibet Plateau， Chinese Academy of Sciences. Long-term continuous monitoring data， such as meteorological factors， temperature and moisture content in the active layer， ground temperature of permafrost in different depths， were recorded in the comprehensive observation station since 2004. It provides a basic-data support for the establishment of permafrost model and the sensitivity analysis of each parameter for the simulation of permafrost temperature field. In this paper， we established an one-dimensional heat conduction permafrost model to assess the influence of different lower boundary conditions on the simulation of the permafrost temperature field in Wudaoliang. In our model， not only it takes the phase change， differences in thermophysical properties between frost and thaw soil， distribution of underground ice in consideration， but also it is flexible enough to set different boundary conditions. Nine different lower boundary condition schemes were designed to evaluate the influence of different lower boundary conditions on the numerical simulation of permafrost temperature field changes from bottom to top quantitatively. The simulated permafrost temperature fields in different schemes were validated by the measured ground temperature series from 2005 to 2015. Moreover， a reconstructed historical ground surface temperature series reaching back to 1960s and a future ground surface temperature series reaching to 2100 were used to force the model to evaluate the influence of climate change on a centennial scale on the simulation of permafrost temperature field from top to bottom quantitatively. The result shows that the temperature field in the near ground surface layer （<3 m） is completely determined by the interannual climate change， while the temperature field of the permafrost in shallow layer （3~15 m） and transitional layer （15~30 m） is affected by lower boundary conditions significantly. The response of temperature field in the deep layer （>30 m） to the climate change on a centennial scale is not only related to the magnitude of climate change， but also directly related to the phase change heat. Improper lower boundary condition may have a great impact on the calculation of the permafrost degradation degree under large-scale climate change， which may lead to serious misjudgement of the temperature field in the deep layer and permafrost area. Therefore， we give some suggestions about the setting of lower boundary condition according to permafrost simulation in different scenarios as follow： if simulation the thickness or of the active layer or the interannual temperature change in the active layer， the lower boundary condition could ignore geothermal energy （set to zero geothermal flux） and be set in transitional layer or below the depth of zero annual amplitude； if simulating the change of mean annual ground temperature （MAGT） of permafrost for ten years or more， the lower boundary condition should be set to constant geothermal flux and be set below the transitional layer； if simulating the permafrost changes on centennial scale， the lower boundary condition should be set to constant geothermal flux and be set in the deep layer or below the permafrost base. This study would provide theoretical information for the selection of the lower boundary conditions in permafrost simulation researches.

Keywords： permafrost ; ground temperature field ; numerical simulation ; lower boundary condition

PDF (4753KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

孙哲, 赵林, 胡国杰, 乔永平, 杜二计, 邹德富, 谢昌卫. 下边界条件对多年冻土温度场变化数值模拟的影响[J]. 冰川冻土, 2021, 43(2): 357-369 doi:10.7522/j.issn.1000-0240.2021.0048

SUN Zhe, ZHAO Lin, HU Guojie, QIAO Yongping, DU Erji, ZOU Defu, XIE Changwei. Influence of lower boundary conditions on the numerical simulation of permafrost temperature field changes[J]. Journal of Glaciology and Geocryology, 2021, 43(2): 357-369 doi:10.7522/j.issn.1000-0240.2021.0048

0 引言

多年冻土约占北半球陆地面积的25%^［1］，其对气候变化尤为敏感^［2］。青藏高原是全球海拔最高的地区，分布面积高达约1.06×10⁶ km²的多年冻土^［3］，并以高温高含冰量为特征^［4-5］。全球气候变暖的背景下，多年冻土呈现退化趋势，表现为多年冻土升温、埋藏深度增大和地下冰融化等，对寒区的工程建筑的稳定性、气候以及水文生态过程都有着重大影响。多年冻土退化导致地基的力学稳定性降低以及引发的地表沉降是影响寒区工程稳定性和安全性的主要因素^［6］；多年冻土与寒区生态系统之间也存在着相辅相生的共生关系，青藏高原多年冻土退化已经引起高寒草甸生态系统发生变化，甚至出现荒漠化、沙漠化^［7-8］；随着多年冻土消融，冻土中的碳储量可能被释放，从而放大高原气候变暖效应^［9］；大量赋存于多年冻土中的地下冰转化为液态水，部分参与到活动层的冻融过程中，而其余部分将参与到区域乃至全球的水循环过程中^［10］；此外，多年冻土退化导致高原下垫面的能量平衡发生变化，最终可能触发整个北半球气候系统的变化^［11-12］。因此，准确评估和预估全球气候变暖背景下大时间尺度多年冻土热状况的变化已经成为冰冻圈科学研究的核心环节之一^［13］。

瞬态数值冻土模式是研究多年冻土温度场变化的一个有效手段，因为其能灵活描述多年冻土剖面的异质性、几何形状和边界条件^［14］。已有多种冻土模式应用于模拟多年冻土的分布及变化^［15-19］。在现有的冻土模式中，考虑相变的热传导方程是计算多年冻土温度场的核心。一般的模拟步骤为首先采用初始恒定的上边界条件进行反复迭代数百次后，直至多年冻土垂直剖面温度场达到稳态作为初始条件，然后赋予相应的上边界条件变化模拟不同情景下冻土温度场的变化^［20］。现有用于模拟多年冻土的研究主要关注近地表温度场变化的影响，主要集中于对气温和浅表层物理过程和参数化方案等改进和优化^{［18，20-21］}。而对于下边界条件设置对模拟多年冻土热状况的影响尚缺乏深入评估。例如大气环流模式中的陆面模式一般采用较浅的下边界，通常为3~10 m，采用零地热通量或恒定温度作为下边界条件^［21-22］，而一些专门针对多年冻土计算的数值模式则采用较深的下边界，通常大于50 m，采用恒定热通量或恒定温度作为下边界条件^{［15，19，23-24］}。对深层多年冻土物理过程不恰当的表现，是当前冻土模式模拟结果主要的不确定性来源之一^{［22，25-28］}。多年冻土是历史气候及地质过程的综合产物，由于土壤是不良的热导体，热传导沿岩土层深度具有滞后性及衰减性。浅层冻土的温度变化主要受控于年际气候变化，而随着岩土层加深，年际气候波动信号会被衰减过滤，深层多年冻土的温度主要取决于长期历史气候以及岩土层下部地热能量共同作用^［29］。因此，如果研究只关注活动层或较浅的多年冻土层年际地温场的变化，也许可以忽视岩土层下部地热能量的影响，相反，如果关注点在于在气候变暖的背景下，模拟和预测百年或千年尺度下多年冻土退化的时机和程度，其取决于深层多年冻土地温场变化，则应必须重视下边界条件设置的影响^［26］。定量界定“浅层”和“深层”多年冻土，同样需要评估在所研究的时间尺度内，不同的下边界条件设置对不同深度多年冻土温度计算的影响程度。

因此，本文拟利用藏北高原冰冻圈站在青藏高原五道梁地区架设的综合观测站的长期地温监测资料，基于一维热传导冻土模型，以重建过去50年的地表温度序列以及未来两种升温情景作为上边界驱动条件，在控制其他变量一致的基础上，通过设计不同的下边界方案进行模拟实验，以定量评估下边界条件对多年冻土地温场数值模拟的影响。分析结果可为冻土模拟研究中下边界条件的选取提供一定参考。

1 研究区与数据

青藏高原五道梁地区（35°13′N、93°05′E，海拔4 783 m）位于昆仑山和唐古拉山脉之间的可可西里地区的东边（图1）。该地区位于高原北部连续多年冻土区，年均气温-5.1 ℃，年降水量为284.8 mm，下垫面为盖度稀疏的高寒荒漠草原。其气候特征与多年冻土特征有较好的代表性^［30-31］。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 研究区地理位置及五道梁综合观测站照片

Fig.1 Location of the study area and the photo showing the comprehensive observation station in Wudaoliang

藏北高原冰冻圈站在五道梁架设了综合观测站，对五道梁气候变化、地气物质能量交换以及多年冻土的动态变化进行了长期的监测及研究^［32-41］。冻土监测包括活动层水热监测和钻孔地温监测。活动层水热监测采用美国Campbell公司的CR1000数采仪监测5~60 cm共12层的土壤温度及水分，监测时段2009—2015年。钻孔深40 m，并未打穿多年冻土，地温探头由冻土工程国家重点实验室加工，精度为±0.05 ℃，测量共26层（0~40 m）的土壤温度，监测时段为2005—2015年。根据2015年的监测数据，该年的地表平均温度为-1.77 ℃，活动层为2.4 m，地温年变化深度（ZAA）大约为7 m，年平均地温（MAGT）为-1.5 ℃。

2 多年冻土数值模型

2.1　考虑相变的一维热传导控制方程

假设土壤中的热量传输完全由傅里叶定律所决定，忽视热对流及热辐射的影响。模型的热传导控制方程表示为

c_{e f f} \frac{\partial T}{\partial t} - \frac{\partial}{\partial z} (λ \frac{\partial T}{\partial z})

=0（1）

式中：c_eff为有效体积热容量（kJ·m^-3·℃^-1）；λ为导热系数（W·m^-1·℃^-1）；T为温度（℃）；t为时间（s）；z为深度（m）。

通过显热容法考虑考虑土壤中的相变过程以及冻土融土热参数的差异。c_eff与λ为温度T的分段函数为

c_{e f f} (T) = \{\begin{array}{l} c_{f}, & T \leq T_{1} \\ \frac{L ∆ w}{T_{2} - T_{1}} + \frac{c_{f} + c_{u}}{2}, & T_{1} < T < T_{2} \\ c_{u}, & T \geq T_{2} \end{array}

（2）

λ (T) = \{\begin{array}{l} λ_{f}, & T \leq T_{1} \\ λ_{f} + \frac{λ_{u} - λ_{f}}{T_{2} - T_{1}} (T - T_{1}), & T_{1} < T < T_{2} \\ λ_{u}, & T \geq T_{2} \end{array}

（3）

式中：下标f与u分别表示土壤的冻结与融化状态；L为单位体积冰水相变潜热（334.54 MJ·m^-3）；Δw为单位体积相变水量；［T₁，T₂］为相变温度区间。

从五道梁观测站现场采集的岩芯土壤样品进行实验分析，主要包括土壤质地、干容重和含水量等。土壤样品通过直径5 cm、高5 cm的采样器从探坑或钻孔岩芯中采集。土壤质地由德国UGT科学仪器公司生产的SEDIMAT4-12粒度分析仪进行分析，其测量原理为湿筛法结合移液管法。干容重和含水/冰量采用烘干法（60 ℃，24 h）测定。模型中所有的土壤热参数基于以上测定的土壤性质资料，参考《冻土地区建筑地基基础设计规范》选取。模型中土壤参数取值见表1。

表1 研究区各岩土层的物理参数

Table 1 Physical parameters of soil layers in the study area

深度/m	容重/（kg·m^-3）	λ_f/（W·m^-1·℃^-1）	λ_u/（W·m^-1·℃^-1）	c_f/（kJ·m^-3·℃^-1）	c_u/（kJ·m^-3·℃^-1）	体积含水量/%	土质
0~0.2	1 650	1.40	1.12	1 979	2 457	3	砂土
0.2~2.0	1 600	2.40	1.90	1 900	2 200	21	粉质砂土
2.0~2.4	1 200	2.22	2.10	2 100	2 240	45	粉质黏土
2.4~3.0	1 550	2.40	1.90	1 950	2 530	30	砾砂
>3.0	2 210	1.45	1.15	1 913	2 227	15	碎石砂土

注：λ_f与λ_u为岩土层冻结和融化时的导热系数；c_f与c_u为岩土层冻结和融化时的体积热容量。

新窗口打开| 下载CSV

2.2　上边界条件

为充分考虑历史气候对多年冻土模拟的影响，本文提取位于五道梁的国家气象站（观测数据（http：//data.cma.cn/）中相应时间段的日气温数据，与综合观测场的日地表温度观测数据进行线性拟合，关系式为

GST=2.274+0.854X， R²=0.89（4）

式中：GST为五道梁综合观测场的实测地表温度；X为五道梁国家气象台站的历史气温序列。然后，根据式（4）与五道梁的国家气象站自1965年的历史日气温数据序列重建过去50年（1965—2015年）的日地表温度序列（图2）。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 重建1965—2015年五道梁日地表温度序列（a）以及地表温度实测值与计算值的对比图（b）

Fig.2 Reconstructed ground surface temperature series from 1965 to 2015 （a） and scatter plots showing measured vs. calculated ground surface temperature during the observation period （b）

上边界的地表温度变化T用下面正弦方程进行拟合。

T = T_{0} + \frac{k}{365} t + A s i n (\frac{2 π}{365} t + φ)

（5）

式中：T₀为初始地表年均温度（℃）；k为升温速率（℃·a^-1）；t为时间（d）；A为地表温度年较差的一半（℃）； $φ$ 为初始相位。根据重建1965—2015年的日地表温度序列的拟合结果，T₀取-2.7 ℃，k取0.025 ℃·a^-1，A取11.7 ℃， $φ$ 取-1.94。

此外，根据第5次耦合模式比较计划（CMIP5）11个气候模式的预测结果，2011—2100年中国气温升温率为0.006~0.063 ℃·a^-1［42］。在RCP8.5情景下，青藏高原地表气温最大升温率最大能达到0.057 ℃·a^-1［43］。在本研究中，为探究百年尺度气候变化对模拟多年冻土温度场变化的影响，假设两个未来（2015—2100年）升温情景：一是假设重建的历史地温升温率（k=0.025 ℃·a^-1）在未来一致持续，二是在RCP8.5极端气候升温情景下假设地表温度与气温以相同速率上升（k=0.057 ℃·a^-1），驱动采用不同下边界方案的模型模拟2015—2100年的未来多年冻土的变化（图3）。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 上边界地表温度的升温趋势

Fig.3 Rise trend of the surface temperature as the upper boundary ground

2.3　数值求解

空间步长设置为多年冻土垂直剖面5 m以上取0.05 m，5 m以下取0.5 m；时间步长设置为1 d。热传导偏微分方程通过MATLAB中ode23t函数（自由内插法的梯形法则）求解。

3 实验设计

3.1　9种不同下边界方案

现有冻土数值模拟计算中，下边界条件的设置差异主要包括设置深度不同和设置方式不同。综合现有主要的冻土模式对于下边界条件的设置，本文设计3个下边界深度：10 m（大气环流模式中的冻土层一般设置的范围^［22］，且位于研究区地温年变化深度以下）、40 m（研究区钻孔底部）和100 m（一些专门针对百年尺度多年冻土变化数值模拟研究所采用下边界深度^{［17， 20］}）；3种下边界条件的设置方式：零地热通量（不考虑地热影响，第二类边界条件）、恒定地热流（第二类边界条件）和相应岩土层的恒定温度（第一类边界条件）。3个深度和3种设置方式进行相互组合，共设计9种下边界条件方案（表2）。根据五道梁地区的深孔（大于60 m）的相关文献^［44-45］，深层地热梯度约为0.03 ℃·m^-1，换算成地热流（方案4~6）。

表2 各方案下边界条件的设置

Table 2 Lower boundary conditions of the nine projects

不考虑地热影响

考虑地热，第二类下边界条件

考虑地热，第一类下边界条件

方案1

z_e=10 m

G=0

方案4

z_e=10 m

G=0.03 ℃·m^-1

方案7

z_e=10 m

T_e=-1.4 ℃

方案2

z_e=40 m

G=0

方案5

z_e=40 m

G=0.03 ℃·m^-1

方案8

z_e=40 m

T_e=-0.8 ℃

方案3

z_e=100 m

G=0

方案6

z_e=100 m

G=0.03 ℃·m^-1

方案9

z_e=10 m

T_e=1 ℃

注：z_e为下边界设置的深度；G为下边界恒定地热梯度；T_e为下边界恒定温度。

新窗口打开| 下载CSV

3.2　实验步骤

第一步：根据9种不同下边界方案，设置9种多年冻土几何垂直剖面。

第二步：在上边界用重建历史地表温度序列的第一年地表温度正弦曲线反复迭代600~1 000次（年），使9种不同方案的多年冻土垂直剖面温度场达到稳态（各层温度的年际变化小于0.001 ℃）。

第三步：以第二步的各多年冻土垂直剖面的稳态温度场作为初始条件，施加重建历史地表温度序列作为上边界驱动，模拟不同方案的多年冻土过去50年的温度场变化。

第四步：在第三步模拟结果的基础上（所为新的初始条件），分别进一步施加k=0.025 ℃·a^-1和k=0.057 ℃·a^-1两种未来地表升温情景，模拟不同方案的多年冻土2015—2100年温度场变化。

第五步：各方案第三、四步的1965—2100年多年冻土温度场模拟结果结合实测值进行相互比较分析。

由于本文主要关注较深的多年冻土地温场的变化，其年内温度振幅较小或为零。故本文使用平均绝对误差（MAE）来判断模拟值与基准值之间的偏差。

M A E = \frac{\sum_{i = 0}^{n} | B_{i} - S_{i} |}{n}

（6）

式中：B_i 为基准值；S_i 为模拟值。

4 结果与分析

4.1　各方案模拟结果及与实测值的对比

图4展示各方案在实验第二步预热结束后多年冻土垂直剖面达到稳态的地温曲线-深度（1965年）、第三步模拟至2015年与实测的2015年地温-深度曲线，以及第四步在两种未来升温情景下模拟至2100年（末年）地温-深度曲线。图5进一步详细展现各方案模拟值与实测值的对比。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 各方案模拟的1965年（预热完成）的地温-深度曲线（a1， b1， c1），模拟和实测的2015年地温-深度曲线（a2， b2， c2），在k=0.025 ℃·a^-1 （a3， b3， c3）和k=0.057 ℃·a^-1 （a4， b4， c4）未来升温情景下模拟的2100年（末年）地温-深度曲线

Fig.4 Comparison of simulated 1965 （initialization complete） temperature-depth curve（a1， b1， c1）， simulated and measured 2015 temperature-depth （a2， b2， c2）， and simulated 2100 （last year） temperature-depth curve as k=0.025 ℃·a^-1 （a3， b3， c3） and k=0.057 ℃·a^-1（a4， b4， c4）

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 各方案3 m （a1， b1， c1）、5 m （a2， b2， c2）、10 m （a3， b3， c3）、16 m （a4， b4， c4）和30 m （a5， b5，c5）地温模拟结果与实测对比（图中MEA值计算以实测值为基准）

Fig.5 Comparison of observation and simulation of the ground temperature in 3 m （a1， b1， c1）， 5 m （a2， b2， c2）， 10 m （a3， b3， c3）， 16 m （a4， b4， c4） and 30 m （a5， b5， c5） from each project （the MAE in the figure calculated based on the observation value）

方案1~3不考虑来自岩土层下部地热能量影响，采用零地热通量作为下边界条件。在预热结束后方案1~3计算的1965年地度-深度曲线基本重合［图4（a1）］，而随着地表升温，至2015年，方案1~3的多年冻土地温模拟结果相对于实测值明显偏小，且下边界设置越深，偏小越显著［图4（a2），图5（a1~a5）］。原因为实验第二步的预热计算的初始稳态地温场的物理意义为在之前百年千年长时间尺度下处于恒定的地气热交换所达到稳定状态的过程，实验第三步反映了气候变暖改变了上边界条件导致多年冻土温度场的变化。因此方案1~3在实验第二步和第三步中忽视了来自岩土层下部地热能量对多年冻土地温场的加热作用，模拟的地温场会显著偏小。而下边界设得越深，更深层的岩土层将吸收更多来自地表升温的能量，故模拟深层岩土层的温度就越偏小^［22］。

从相变热传导理论的角度，采用恒定地热流作为第二类下边界条件的方案4~6相对更符合实际的物理过程 ^{［14，26］}。在预热结束后方案4~6计算的1965年地温-深度曲线也基本重合［图4（b1）］，而随着地表升温，方案4~6的多年冻土2015年地温模拟结果能较好的符合实测值［图4（b2）］，也能较好的再现各岩土层实测近十年的升温情况［图5（b1~b5）］。就方案4~6相互比较而言，下边界深度设置越浅，多年冻土地温模拟结果就会越偏高。这是因为在实验第三步地温自上而下逐步升高时，上层的地温梯度会减小，同时又受到岩土层下部地热流作用，使得整个地温曲线随时间呈现“弓形”^［46］。实测的上层多年冻土垂直剖面的地温梯度小于深层的垂直剖面地温梯度。若下边界深度设得浅，又采用深层的地温梯度换算的恒定地热流作为下边界条件，则会高估来自岩土层下部地热流输入的能量，进而模拟出偏高的多年冻土温度。

采用2015年恒定地温作为第一类下边界条件的方案7~9的模拟结果能很好的符合2015年实测的地温场［图4（c1）］。但由于固定了下边界的温度，其无法反映下边界附近岩土层的升温情况［图5（c1~c5）］。方案8把下边界设在40 m，方案9设为在100 m，与各岩土层实测的近十年升温速率相比，方案8中模拟8 m以下岩土层的升温速率相对于方案9与实测值开始偏小，越往深的岩土层，偏小幅度越大；而将下边界设置与10 m的方案7则无法模拟5 m以下实测的升温变化。特别是上边界有明显的升温或降温的情景时，定温边界会对模拟结果有很大影响。如图4（c3）和图4（c4）所示，当模拟未来百年尺度多年冻土地温场变化时，无论是k=0.025 ℃·a^-1还是k=0.057 ℃·a^-1升温情景，方案9显示10 m 和40 m 的岩土层温度会发生较大的变化，故采用较浅深度的定温边界条件的方案7和8所预估的结果是不具备参考价值的。

4.2　以方案6为基准各方案的相互对比

基于研究目的，本文近地表层的各种陆面过程考虑得比较简单，使得各方案的近地表层温度模拟结果以实测值作为基准值时MAE较大。方案6在100 m深度以恒定地热流作为下边界条件，相比其他方案更接近实际的岩土层下部能量传输物理过程，其模拟结果也能更好的重现研究区实测的多年冻土温度场。为从理论上定量评估不同下边界设置对计算多年冻土温度场的影响程度，以方案6的模拟结果作为基准值（理论基准值），定量比较其他方案不同深度岩土层温度模拟结果相对于方案6的MAE（表3），同时也参考各方案模拟结果以实测为基准值时的MAE（图5），从热传导的角度对多年冻土的垂直剖面进行分层。初步设定3个MAE的阈值：0.1 ℃、0.5 ℃以及1 ℃。根据目前使用的温度探头的精度为±0.05 ℃，以及多年冻土的工程分类：极不稳定（>-0.5 ℃），不稳定（-1～-0.5 ℃），稳定过渡（-3～-1 ℃）和稳定（<-3 ℃）；将0.1 ℃的MAE作为在仪器误差范围内判定模拟值与基准值是否有偏差的阈值，将0.5 ℃的MAE作为判定模拟值与基准值是否有显著偏差的阈值（偏差大于0.5 ℃就可能误判多年冻土分类），将1 ℃作为模拟值与基准值是否有极显著偏差的阈值。

表3 以方案6为基准值，计算的各方案历史模拟结果的平均绝对误差

Table 3 The MAE of historic simulation from each project calculated based on the sixth project as the truth value

深度/m	MAE/℃
深度/m	方案1	方案2	方案3	方案4	方案5	方案7	方案8	方案9
1	0.01	0.02	0.02	0.03	0	0.03	<0.01	<0.01
2	0.04	0.06	0.06	0.07	0.01	0.09	0.01	<0.01
3	0.06	0.09	0.09	0.10	0.01	0.13	0.01	<0.01
4	0.08	0.12	0.12	0.14	0.01	0.17	0.02	<0.01
5	0.10	0.15	0.16	0.18	0.01	0.21	0.02	<0.01
6	0.12	0.18	0.19	0.22	0.01	0.25	0.03	<0.01
8	0.16	0.25	0.25	0.30	0.02	0.33	0.04	<0.01
10	0.20	0.31	0.32	0.40	0.02	0.39	0.05	<0.01
12		0.37	0.38		0.03		0.06	<0.01
14		0.43	0.44		0.03		0.07	<0.01
16		0.49	0.51		0.04		0.08	<0.01
18		0.55	0.57		0.05		0.09	<0.01
20		0.61	0.63		0.05		0.11	<0.01
25		0.76	0.79		0.07		0.16	<0.01
30		0.91	0.95		0.09		0.21	0.01
35		1.05	1.11		0.13		0.26	0.01
40		1.20	1.27		0.20		0.31	0.01

新窗口打开| 下载CSV

如表3所示，以方案6为基准值时，其他方案的3 m以上岩土层温度场模拟结果的MAE均小于0.1 ℃，由于所有方案的上边界条件设置相同，说明3 m以上岩土层温度变化基本完全由上边界年际气候变化所决定，从热传导的角度可把3 m以上的岩土层划为“近地表层”，其不受岩土层下部的地热能量影响。对于不考虑地热影响的方案1~3，其3~15 m岩土层温度模拟结果的MAE在0.1~0.5 ℃之间（若以实测为基准值，MAE在0.27~0.44 ℃之间），说明该深度范围的岩土层温度变化虽然还是主要受控于年际气候变化，但随着深度增加，下边界设置的影响开始逐渐显现，该深度范围的岩土层可划为“浅层”；15~30 m岩土层温度模拟结果的MAE在0.5~1 ℃之间（若以实测为基准值，MAE在0.44~0.97 ℃之间），说明该深度范围的岩土层温度变化已经开始受到岩土层下部地热能量的显著影响。青藏高原的多年冻土地温年变化深度通常在12~20 m之间^［47］，计算15~30 m的多年冻土温度时，如果不考虑地热影响，可能会导致较大的误差进而误判多年冻土的类型，该深度范围可划为“中层”；30 m以下岩土层的温度模拟结果的MAE已经大于1 ℃（若以实测为基准值，方案3的MAE同样大于1 ℃），说明使用不考虑地热影响的方案计算30 m以下岩土层的温度应当是不可取的。因此，30 m以下岩的深度范围可以划为“深层”。

对于同样以恒定地热流作为第二类下边界条件的方案4和5，下边界设置10 m深度的方案4的浅层地温模拟结果的MAE在0.1~0.5 ℃之间；下边界设置40 m深度的方案5的深层地温模拟结果的MAE开始大于0.1 ℃（若以实测为基准值，则20 m以下岩土层的MAE开始大于0.1 ℃）。

对于以恒定温度作为第一类下边界条件的方案7~9，同样将下边界设置在100 m处的方案9与作为基准的方案6的模拟结果几乎完全一致（MAE<0.01 ℃）；将下边界设置在40 m处的方案8中，20 m以下岩土层温度模拟结果开始显现较小的偏差（0.1 ℃<MAE<0.31 ℃）；而将下边界设置在10 m处的方案7中，则3 m以下（浅层范围内）地温模拟结果开始出现偏差（MAE>0.1 ℃），7 m以下地温模拟结果的MAE大于0.3 ℃。

4.3　百年尺度气候变化对多年冻土垂直剖面温度变化的影响

下边界深度应设置在地热流或温度比较稳定的岩土层，即在所研究的时间尺度内，上边界气候变化难以影响到的岩土层。目前对于气候变暖的背景下青藏高原多年冻土变化的模拟及预估的研究，主要还是集中于近百年尺度内^［48］。因此，下面探讨在两种未来升温情景下，不同下边界条件的设置对百年尺度多年冻土温度场模拟的影响。

表4和表5分别展示各方案模拟在k=0.025 ℃·a^-1和k=0.057 ℃·a^-1两种未来升温情景下不同深度的岩土层未来近100年来总升温量。可以明显看出恒定温度（第一类边界条件）作为下边界条件的方案7、8的对靠近下边界岩土层的升温变化模拟的局限性。若以方案6为基准，在两种未来升温情景下，70 m以下岩土层的升温总量小于0.01 ℃，这也解释了采用100 m恒定温度作为下边界条件的方案9与方案6模拟结果基本一致的原因。但是方案3的模拟结果显示，70 m以下岩土层的升温总量均大于1 ℃，百年尺度内的气候变化依旧会影响这些深层岩土层的温度。在k=0.057 ℃·a^-1升温情景下，方案6和9的模拟结果显示升温的近100年间多年冻土消融量（多年冻土厚度减小）为20.15 m，而方案3显示升温的近100年间多年冻土消融量为15.20 m。这说明百年尺度上边界变化对深层多年冻土温度的影响不仅与气候变化的幅度有关，而且与多年冻土消融产生的相变热的多少有也直接的关系，气候变化信号导致传到深层导致相应位置的负热流（热流方向以向上为正）发生了改变，进而引起岩土层的温度发生变化，若多年冻土相变耗热多，则传到此处的负热流就较小，所以深层地温变化小，如果相变热小，深层地温可能就会有显著的变化。

表4 各方案模拟在k=0.025 ℃·a^-1升温情景下多年冻土垂直剖面2015—2100年总升温量

Table 4 Simulated warming of the geothermal field from the nine projects during 2015—2100 as k=0.025 ℃·a^-1

深度/m	总升温量/℃
深度/m	方案1	方案2	方案3	方案4	方案5	方案6	方案7	方案8	方案9
1	1.89	1.89	1.89	1.87	1.87	1.87	1.81	1.87	1.87
3	1.62	1.64	1.63	1.54	1.55	1.56	1.32	1.54	1.56
5	1.41	1.45	1.44	1.27	1.30	1.31	0.93	1.28	1.31
10	1.31	1.40	1.39	1.10	1.09	1.13	0	1.06	1.12
20		1.48	1.37		0.85	0.95		0.72	0.95
30		1.53	1.32		0.63	0.76		0.37	0.75
40		1.55	1.26		0.53	0.54		0	0.53
50			1.20			0.32			0.31
60			1.14			0.16			0.16
70			1.09			0.10			0.12
80			1.05			0.09			0.08
90			1.02			0.09			0.04
100			1.01			0.09			0

新窗口打开| 下载CSV

表5 各方案模拟在k=0.057 ℃·a^-1升温情景下多年冻土垂直剖面2015—2100年总升温量

Table 5 Simulated warming of the geothermal field from the nine projects during 2015—2100 as k=0.057 ℃·a^-1

深度/m	总升温量/℃
深度/m	方案1	方案2	方案3	方案4	方案5	方案6	方案7	方案8	方案9
1	5.07	4.96	4.96	5.07	4.94	4.94	4.74	4.93	4.95
3	4.84	4.43	4.40	4.87	4.35	4.36	3.55	4.31	4.36
5	4.65	3.96	3.91	4.70	3.83	3.84	2.47	3.76	3.83
10	4.45	3.02	2.92	4.52	2.74	2.77	0	2.63	2.77
20		1.60	1.61		1.01	1.09		0.95	1.08
30		1.54	1.54		0.65	0.76		0.47	0.75
40		1.56	1.48		0.56	0.55		0	0.53
50			1.41			0.33			0.32
60			1.35			0.17			0.17
70			1.30			0.11			0.13
80			1.26			0.10			0.08
90			1.23			0.10			0.04
100			1.22			0.09			0

新窗口打开| 下载CSV

各方案预估结果的比较显示，当上边界有长时间尺度显著升温的情景时，不同的下边界设置会对多年冻土消融及其产生的相变热的计算会有显著影响，进而模拟结果中百年尺度气候变化对深层多年冻土地温的影响程度也会有很大的差异。如图4（a4）和图4（b4）所示的结果，模拟RCP8.5升温情景（k=0.057 ℃·a^-1）下多年冻土百年尺度下未来变化时，若采用第二类边界条件且下边界设在浅层的方案1和方案4的模拟结果显示到2100年五道梁地区的多年冻土以完全消融，变为非多年冻土地区；然而若下边界位置设置在深层以下的方案2~5的模拟结果显示到2100年五道梁地区依旧存在着残余多年冻土，不能划为非多年冻土区。目前气候模式中计算的多年冻土层设置很薄的情况，无疑会对气候变暖的背景下全球多年冻土面积变化的预估结果造成影响并引发较大的争议^［26］。

5 讨论与结论

本文基于一维热传导冻土温度数值模型，以青藏高原五道梁地区多年冻土为背景，设计9种不同下边界条件设置方案进行控制模拟实验，定量评估在百年尺度气候变化下不同下边界条件设置对多年冻土温度场模拟的影响。结果显示，五道梁地区近表层（<3 m）的地温完全受控于年际气候变化，浅层（3~15 m）和中层（15~30 m）的多年冻土温度场受下边界条件的影响逐渐显著，下边界条件的设置对计算百年尺度气候变化下多年冻土消融相变热以及深层（>30 m）的温度场具有重大影响。不同多年冻土地区由于下垫面、地层性质和地表温度振幅等差异，可能导致具体的分层情况有所差异，限于篇幅原因，本文未探讨这些因素对下边界条件设置的影响，这也是下一步研究的关注点。

最后根据不同的条件与目的，对于采用基于一维热传导方程为核心的冻土模式对多年冻土温度场进行模拟的研究，本文给出如下的下边界条件设置建议：

（1）在条件允许的情况，建议采用实测的深层恒定地热流作为下边界条件且下边界设置在深层，使得多年冻土温度场模拟结果既能反映岩土层下部的地热能量的影响又能反映在气候影响下的变化的趋势。

（2）近表层的地温完全受控于年际气候变化，而近表层和浅层的热流为负值，且随深度递减、周期波动，在这个深度范围内采用单一热流边界并不合适。如果研究关注当前研究区活动层厚度或年际活动层温度场变化，考虑计算效率，可以忽略岩土层下部地热能量影响，建议下边界采用零地热通量，设置在中层或地温年变化深度以下。

（3）如果研究关注十年或十年以上多年冻土年平均温度（MAGT）的变化，应当不可忽视来自岩土层下部地热能量，建议下边界采用恒定地热流，设置在中层或深层。恒定温度（第一类边界条件）作为下边界条件时，无法反映下边界位置附近的多年冻土温度的变化。因此，仅仅对于地温年际没有变化或变化很小的情景，恒温边界才是可行的。对于处于零地热梯度相变过渡型多年冻土^［49］，则建议下边界设置在多年冻土底板以下的融土层。

（4）如果研究关注模拟及预估世纪尺度或更长时间尺度，多年冻土随时间退化甚至完全消失过程，则应慎重考虑下边界条件设置对多年冻土消融相变热以及深层多年冻土温度场计算的影响，建议下边界设置在深层或多年冻土底板以下的融土层，采用恒定地热流作为下边界条件。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Zhang

Tingjun

， Barry

R G

， Knowles

， et al.

Statistics and characteristics of permafrost and ground-ice distribution in the Northern Hemisphere

［J］. Polar Geography， 1999， 23（2）： 147-149.