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

孙 哲， 赵 林， 胡国杰， 乔永平， 杜二计， 邹德富， 谢昌卫

（1.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室/藏北高原冰冻圈特殊环境与灾害国家野外科学观测研究站，甘肃 兰州 730000； 2.中国科学院大学，北京 100049； 3.南京信息工程大学地理科学学院，江苏 南京 210044）

0 引言

多年冻土约占北半球陆地面积的25%［1］，其对气候变化尤为敏感［2］。青藏高原是全球海拔最高的地区，分布面积高达约1.06×106km2的多年冻土［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 研究区地理位置及五道梁综合观测站照片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 考虑相变的一维热传导控制方程

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

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

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

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

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

表1 研究区各岩土层的物理参数Table 1 Physical parameters of soil layers in the study area

2.2 上边界条件

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

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

图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用下面正弦方程进行拟合。

式中：T0为初始地表年均温度（℃）；k为升温速率（℃·a-1）；t为时间（d）；A为地表温度年较差的一半（℃）；φ为初始相位。根据重建1965—2015 年的日地表温度序列的拟合结果，T0取-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 上边界地表温度的升温趋势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

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）来判断模拟值与基准值之间的偏差。

式中：Bi为基准值；Si为模拟值。

4 结果与分析

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

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

图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 各方案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

如表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

表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

各方案预估结果的比较显示，当上边界有长时间尺度显著升温的情景时，不同的下边界设置会对多年冻土消融及其产生的相变热的计算会有显著影响，进而模拟结果中百年尺度气候变化对深层多年冻土地温的影响程度也会有很大的差异。如图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）如果研究关注模拟及预估世纪尺度或更长时间尺度，多年冻土随时间退化甚至完全消失过程，则应慎重考虑下边界条件设置对多年冻土消融相变热以及深层多年冻土温度场计算的影响，建议下边界设置在深层或多年冻土底板以下的融土层，采用恒定地热流作为下边界条件。