川西螺髻山清水沟倒数第二次冰期以来的冰川规模与古气候重建

唐倩玉， 张 威， 刘 亮， 柴 乐， 李亚鹏，张廉卿， 孙 波， 张宏杰

（1.辽宁师范大学地理科学学院，辽宁大连 116029； 2.东华理工大学地球科学学院，江西南昌 330013；3.中国科学院青藏高原研究所青藏高原环境变化与地表过程重点实验室，北京 100101）

0 引言

第四纪全球气候的显著特征之一是冰期和间冰期交替，丰富的第四纪冰川作用遗迹为探讨第四纪期间的气候环境变化特征与形成机制提供了良好的载体［1-2］，第四纪也是人类出现和迅速发展的时代，研究这一时期的气候变化对探讨自然环境变化对人类的影响有重要意义。随着地貌与第四纪研究的不断深入，采用相应的数学模型，并辅助地理信息系统（geographic information system，简称GIS），第四纪冰川研究已经进入了定量化研究阶段［3-10］。尤其是根据古今冰川平衡线高度变化，定量化推算不同冰期时段的气温与降水，明确过去和现在气候变化的差异，可为未来气候变化趋势预测提供参考［11-13］。

螺髻山位于青藏高原东缘、横断山的中部，气候受东南季风和西南季风的共同影响，加上高大山体对气流的抬升作用使得这里降水充沛，11 月至3月该地区受高空西风环流南支急流控制，以晴朗天气为主，降水少，4 月东南季风和西南季风北上，雨季开始，9—10 月印度洋暖湿气流更强烈的侵入成为年降水的第二个高峰，夏半年的降水量占全年的75%～90%［14］。该区域是青藏高原东缘有明确古冰川地形而无现代冰川发育的地区之一［2］，在第四纪期间经历了多次冰川作用，海拔2 500 m 以上保留有序列齐全、对气候变化敏感的海洋性冰川遗迹。近年来，一些学者对螺髻山第四纪冰川地貌进行了深入考察与研究，进行了第四纪冰期系列划分并建立了绝对年代学框架［15-17］（图1）。依据不同期次的冰川遗迹，运用冰川纵剖面模型、积累区面积比率法（accumulation area ratio，简称AAR）和面积-高程平衡率法（area altitude balance ratio，简称AABR），以及气温与降水关系模型（简称P-T 模型）、气温递减率模型（lapse-rate model，简称LR）等，可以获得不同阶段古冰川规模及冰川物质平衡线高度（equilibrium-line altitude，简称ELA），进而得到不同冰川作用期的气候环境演变特征。

图1 螺髻山清水沟地理位置与冰川地貌Fig.1 Geographical location and glacial landforms of the Qingshui Valley，Luoji Mountain

1 研究区概况

清水沟谷地的第四纪冰川沉积物保存较好，从源区冰斗到沟谷谷口依次保存有三套冰碛物（图1），张威等［15］对这三套冰碛物的测年结果显示：第一套冰碛物上采集两个电子自旋共振样品（electron spin resonance，简称ESR）测年结果为（13±1）ka 和（23±2）ka，以及两个OSL（光释光）测年结果为（15.74±1.15）ka 和（23.24±1.68）ka，因此根据相对地貌和测年结果可确定第一套冰碛物形成于末次冰期晚期，对应于海洋氧同位素阶段（marine oxygen isotope stage，简称MIS）2。第二套冰碛物分布于海拔3 400～3 550 m大海子冰坎下部，终碛与侧碛相连，部分终碛被后期流水冲开，第二套冰碛上的两个ESR 测年结果为（72±30）ka 和（81±11）ka，结合地貌可判断其形成于末次冰期早期（MIS 4）。第三套冰碛物的侧碛垄从海拔3 400 m延伸至海拔2 300 m，终碛堤位于清水沟西北2 200 m 的山包上，横剖面呈月牙状堆积体，相对高度50 m［18］。在冰碛垄的两个样品ESR 年龄为（234±45）ka 和（191±157）ka，结合地貌判断其形成于倒数第二次冰期（MIS 6）。倒数第二次冰期冰川从大海子冰斗（下层冰斗）向下游溢出形成了上槽谷，从海拔3 740 m 延伸到2 500 m 左右，末次冰期早期形成了清水沟3 800 m左右的冰斗群以及下槽谷［19］。

2 研究方法

2.1 冰川规模定量恢复方法

使用Pellitero 等［20］开发的冰川规模重建工具箱Glacier Reconstruction Tools（简称GLARE），恢复冰川表面高程时需要古冰川底部的地形数据和古冰川分布范围。倒数第二次冰期的冰川范围以下层冰斗后壁、上槽谷、侧碛垄以及谷口终碛为参考，末次冰期早期的冰川范围以上层冰斗后壁、下槽谷、侧碛垄为参考，末次冰期晚期的冰川范围以上冰斗附近的侧碛和终碛为参考。先沿着冰川中流线计算冰川厚度，之后插值至整个冰川分布范围，最终得到整条冰川的厚度。冰川纵剖面模型是根据冰川动力学方程建立的［21］。

式中：τ为剪切应力；ρ为冰川冰的密度；g为重力加速度；h为冰川的厚度；α为冰川表面的坡度。

具体步骤为：①使用Flowline Creation 工具，生成古冰川中流线；②使用Shear Stress 工具中的Construct Interval Nodes 得到中流线上的节点；③使用Shear Stress 工具中的Define Shear Stress 对节点的剪切应力进行赋值；④使用Glacier Reconstruction中的Flowline ice thickness tool 得到流线上的冰川海拔高度；⑤考虑槽谷形态对冰川的影响，使用F因子中的Automatic ice thickness 工具对流线上的冰川海拔高度进行校正；⑥使用Glacier Reconstruction中的Glacier surface interpolation 输入中流线和冰川作用范围，选择克里金插值方法，得到冰面高程；⑦冰面高程与冰川底部地形相减可以得到冰川厚度；⑧通过得到的冰面高程可以在ArcGIS 上计算冰川表面高程栅格数据的面积；⑨在ArcGIS上计算冰川厚度乘以冰川表面插值网格面积得出冰川体积。

2.2 冰川平衡线估算方法

ELA 是指冰川上积累量与消融量相等的点连线的高度［22-24］。本文采用目前较为常用的AAR 和AABR法来恢复清水沟古冰川物质平衡线高度。当冰川处于稳定状态时，冰川积累区面积与整个冰川面积的比值即为AAR 值［25］，古冰川ELA 的计算参考Kern 等［26］提出的AAR 值与冰川面积间的对数关系式。

AABR 法是由Furbish 等［27］提出的根据冰川消融区和积累区的物质平衡梯度来计算冰川物质平衡线高度的方法，Rea［28］在计算全球65 条冰川的物质平衡时发现BR 的平均值为1.75±0.71，因此本文参考此平均值计算清水沟的BR 值。本文基于恢复后的古冰面高程计算不同时期的古ELA，可以使计算结果更为精确。

2.3 基于冰川平衡线变化的气候重建模型

（1）P-T模型

气温和降水是影响冰川ELA 变化的重要因素［29-30］，因而现代冰川ELA 处气温和降水的关系也被广泛地应用于气候重建［31-32］。古冰川运用P-T 模型时要结合ELA 的变化量及气温递减率，从而得到因为ELA 位置的变化而产生的气温变化量。在假定古冰川前进时ELA 和现代冰川ELA 处气温和降水关系一致的基础上，结合气候代用指标得到降水的变化量，带入到P-T 模型中得到气温变化量。然后将位置变化产生的气温变化量和降水波动产生的气温变化量相加，就是气温波动值［33］。本文在恢复不同阶段古气候时参考了施雅风等［29］的研究，建立了中国西部冰川年降水量P（mm）和年均温T（℃）之间的关系式。

（2）LR模型在使用LR 模型时要考虑温度递减率、积累梯度和积累量的变化［34］。该模型应用ELA 处冰川积累量与温度的转换系数f，将由于ELA 升降导致的降水变化量转化为相应的气温变化量［33］。

式中：ELA1和ELA2分别为现代和冰期时的物质平衡线高度；T1为现代冰川ELA2处的气温（℃）；T2为冰期时古冰川ELA2处的气温（℃）；c1为现代冰川在ELA2处的积累量（mm）；c2为冰期时古冰川在ELA2处的积累量（mm）；f为物质平衡线处冰川积累量与温度的转换系数（℃·mm-1）。

3 结果与讨论

3.1 冰川规模估算

利用GLARE 工具恢复螺髻山清水沟倒数第二次冰期的冰川山谷指数为0.45～0.87［图2（a）］，平均值为0.62；剪切应力范围为40～140 kPa［图2（b）］，平均值为102.10 kPa；冰川表面高程为2 125～3 826 m［图3（a）］，平均值为3 029 m，最大冰川厚度为127 m［图3（b）］，平均值为55 m，冰川表面积为3.44 km2，体积为0.19 km3。恢复末次冰期早期的冰川山谷指数为0.53～0.96［图2（c）］，平均值为0.77；剪切应力范围为50～70 kPa［图2（d）］，平均值为64.65 kPa；冰川表面高程为3 380～4 055 m［图3（c）］，平均值为3 054 m，最大冰川厚度为131 m［图3（d）］，平均值为53 m，冰川表面积为2.22 km2，体积为0.12 km3。末次冰期晚期的冰川山谷指数为0.49～0.95［图2（e）］，平均值为0.79；剪切应力范围为40～70 kPa［图2（f）］，平均值为57.71 kPa；冰川表面高程为3 765～4 043 m［图3（e）］，平均值为3 016 m，最大冰川厚度为121 m［图3（f）］，平均值为57 m，冰川表面积为1.20 km2，体积为0.07 km3。

图2 清水沟不同冰进时期的山谷指数和剪切应力Fig.2 Shape index（a，c，e）and shear stress（b，d，f）of Penultimate Glaciation，early stage of Last Glaciation，and late stage of Last Glaciation in the Qingshui Valley

图3 清水沟不同冰进时期的冰面高程和冰川厚度分布Fig.3 Distribution of ice surface elevations（a，c，e）and glacier thickness（b，d，f）of Penultimate Glaciation，early stage of Last Glaciation，and late stage of Last Glaciation in the Qingshui Valley

测量的山谷指数在末次冰期晚期最大，末次冰期早期居中，倒数第二次冰期最小，这主要是因为靠近冰川源头受到的冰川侵蚀更强烈、持续时间也更长，使槽谷底部更趋近于“U”形［35］。对于山谷冰川来说，剪切应力在50～150 kPa 的范围内［36］，恢复的各阶段冰川底部剪切应力平均值均在此范围内。研究区倒数第二次冰期以来冰川的面积和体积都不断减小，这种第四纪以来冰期规模逐渐减小的特点在中国山地冰川区具有普遍性［37］。

3.2 不同时期冰川平衡线高度

将恢复出的冰川面积带入式（2），计算出倒数第二次冰期AAR 值为0.56，末次冰期早期AAR 值为0.53，末次冰期晚期AAR 值为0.49，然后使用Pellitero 等［20］开发的ELA calculation 工具输入AAR值和BR 值可自动计算对应的ELA 高度。从表1 看出，运用AABR 法、AAR 法，结合前人使用的冰斗底部高程法（CF）、末端至冰斗后壁比率法（MELM）、侧碛垄最大高度法（THAR）和TSAM 法［18］计算出清水沟倒数第二次冰期ELA 平均值为3 132 m，末次冰期早期ELA 平均值为3 776 m，末次冰期晚期ELA 平均值为3 927 m。根据高晓昕［18］计算的结果，螺髻山现代冰川平衡线高度为4 848 m，倒数第二次冰期相对于现代ELA 下降1 716 m，末次冰期早期下降1 071 m，末次冰期晚期下降920 m。

表1 研究区MIS 6以来各次冰进ELA重建Table 1 ELA reconstruction in each glacial advance since MIS 6 in the study area

3.3 古气候环境重建

收集了研究区附近的西昌、昭觉、木里、盐源、九龙和稻城这六个气象站1990—2019 年每日的气温和降水数据（表2），所用的气象数据来自国家气象科学数据中心（http://data.cma.cn/）的中国地面气候资料日值数据集（V3.0）。

表2 研究区附近气象站气候数据概要Table 2 Summary of present climate datasets from meteorological stations used near the study area

对研究区附近气象站的多年平均气温T（℃）和海拔H（m），以及多年平均降水量P（mm）和海拔H（m）进行拟合。

使用P-T 模型对气温变化进行恢复的方法如下［33］：通过式（5）得到气温递减率为0.54 ℃·（100m）-1，结合倒数第二次冰期ELA下降高度，先得出因为位置变化导致的气温变化为9.27 ℃；然后根据式（6）计算出现代ELA 处的降水为512.38 mm，根据式（3）计算出现代ELA 处的气温为-11.10 ℃，位于青藏高原东北部若尔盖盆地RM 钻孔的孢粉分析结果［38］显示，MIS 6 年降水量为现在的60%～80%，因此在MIS 6 时现代ELA 处的降水为307.43～409.90 mm，使用式（3）得到ELA 处的年均温为-13.70～-12.23 ℃，所以降水波动产生的气温变化量为1.13～2.59 ℃；位置和降水波动产生的气温变化量相加就是气温的波动值，因此得出在倒数第二次冰期时气温比现在低10～12 ℃。末次冰期早期和晚期的降水量分别为现在的80% 左右和60%～80%［38］，按照上述方法得到的计算结果如表3 所示。图4（a）模拟了P-T模型不同冰进阶段年降水量占现代年降水量的比例由10% 变化到200%，变化的步长为10%时的气候状况。

表3 研究区MIS 6以来各次冰进的古气候状况Table 3 Palaeoclimatic condition in each glacial advance since MIS 6 in the study area

使用LR 模型对气温变化进行恢复的方法见文献［39］。LR 模型与P-T 模型相比，还需要冰川积累梯度和不同阶段ELA的乘积，以及ELA处冰川积累量与温度的转换系数（f）。由式（6）可以得到冰川积累梯度和不同阶段ELA 的乘积。根据Ohmura 等［30］提出的f值，本文采用0.0029 ℃·mm-1。降水比例与P-T 模型设置为一致，将上述数据带入LR 模型计算，结果如图4（b）和表3所示。

图4 模型模拟结果Fig.4 Simulation results of P-T（a）and LR（b）models

根据P-T 模型和LR 模型的结果（表3），可以看出：倒数第二次冰期气温降低值为8～12 ℃，这与根据冰楔假型恢复出的鄂尔多斯倒数第二次冰期时冻土边界附近的年均温比现在低10～15 ℃相符［40］，该阶段降水仅为现在的60%～80%，所以气温大幅度的下降是导致该阶段冰川规模最大的原因；末次冰期早期气温降低值为6～7 ℃，与若尔盖盆地孢粉记录显示MIS 4 阶段年平均气温比现在低4 ℃相差不大［38］，该阶段降水为现在的80% 左右，但是气温下降幅度小于倒数第二次冰期，因而冰川规模小于倒数第二次冰期；末次冰期晚期气温降低值为4～8 ℃，与若尔盖盆地孢粉记录显示MIS 2阶段的年平均气温比现在低5～6 ℃的情况一致［38］，末次冰期晚期降水仅为现在的60%～80%，降温幅度也不大，因此该阶段冰川规模最小。

4 结论

螺髻山清水沟第四纪冰川规模与古环境重建的主要结论如下：

（1）基于冰碛物的分布范围，利用冰川纵剖面模型恢复了古冰川的规模，清水沟倒数第二次冰期的冰川面积为3.44 km2，最大冰川厚度为127 m，冰川体积为0.19 km3；末次冰期早期的冰川面积为2.22 km2，最大冰川厚度为131 m，冰川体积为0.12 km3；末次冰期晚期的冰川面积为1.20 km2，最大冰川厚度为121 m，冰川体积为0.07 km3。

（2）运用AABR 法和AAR 法，并结合冰斗底部高程法（CF）、末端至冰斗后壁比率法（MELM）、侧碛垄最大高度法（THAR）和TSAM 法恢复了古冰川平衡线高度，结果显示清水沟倒数第二次冰期ELA高度为3 132 m，末次冰期早期ELA 高度为3 776 m，末次冰期晚期ELA高度为3 927 m。

（3）运用P-T 模型和LR 模型，结合气候代用指标，可以得到研究区倒数第二次冰期气温降低值为8～12 ℃，末次冰期早期气温降低值为6～7 ℃，末次冰期晚期气温降低值为4～8 ℃。