张 越, 许向科, 孙雅晴
(1.中国科学院青藏高原研究所青藏高原地球系统与资源环境全国重点实验室,北京 100101;2.中国科学院大学,北京 100049)
冰川对气候变化非常敏感[1],青藏高原是地球上除南北两极外冰川发育最广泛的区域[2]。末次冰盛期(Last Glacial Maximum,LGM)发生在距今18~24 ka[3],期间全球大范围降温,两极以及山地冰川大规模扩张。根据古冰川作用形成的地貌能够定量重建古冰川的范围、冰储量和平衡线高度(equilibrium-line altitude,ELA)等参数。将其与现代冰川的相关参数进行对比,尤其是平衡线高度的变化,能够进一步应用冰川气候模型可以直接地定量地重建冰期时的气候条件,最终可以和湖泊孢粉记录等气候代用指标[4-11]以及大气环流模型[12-15]的结果互相对比验证。这也是研究冰期时古气候的一种重要思路,对了解高原上冰川-气候相互作用机制有重要的意义。
青藏高原东南部是高原上海洋型冰川发育最广泛的地区之一[16],主要受印度季风控制[17],是了解印度季风与冰川变化关系的关键区域。Hu等[18-19]通过10Be暴露测年方法测得了青藏高原东南地区巴松措流域与南迦巴瓦峰附近的派山谷冰川作用的时间。另外,在巴松措流域,Hu等[18]还应用光释光(OSL)测年技术测定了冰碛物的年代。其中,10Be暴露测年法所测得的是冰川形成该规模的最小年代,而OSL法测得的是冰碛物的沉积年代,Hu等[18]综合判定此次冰川作用发生在LGM时期,并根据侧碛-终碛等地貌大致限定了LGM冰川作用的范围,并估算了当时的ELA。由于大多情况下冰川地貌存在缺失或者后期地质改造而退化的情况[20],而冰川平衡剖面模型仅需要部分保存下来的地貌证据即可完成对古冰川规模的重建,因此最好将这些地貌、地形证据与数值模型方法相结合用于古冰川的重建工作[21]。对于稳定态冰川的模拟,冰川纵剖面模型与高阶冰流模型运行结果极其相似[22],而且所需参数少且易于获得,运算简单快捷。因此,近年来冰川平衡纵剖面模型在对青藏高原范围内的古冰川重建方面受到了广泛的应用[23-32]。本文将以前人的10Be测年结果为年代学基础,应用冰川平衡纵剖面模型结合有关地貌证据,定量重建这两条冰川在LGM时期的范围、冰储量和ELA等参数,由此进一步恢复该地区LGM时期的气候条件。
青藏高原东南部北接念青唐古拉山脉,西接喜马拉雅山脉,东接横断山脉,最高峰为南迦巴瓦峰,海拔高度为7 782 m(图1)。通过国家气象科学数据中心(http://data.cma.cn/)获得的距离研究区最近的林芝(海拔2 993.0 m)、波密(海拔2 737.2 m)、米林(海拔2 951.2 m)3个气象站的气象观测记录,记录显示年平均气温(MAAT)分别为9.1℃、9.0℃和8.8℃,平均年降水量分别为729.7 mm、890.1 mm和708.8 mm,降水主要集中在每年的4月到9月,约占总降水量的73.4%~85.8%(图2),印度季风是该地区主要的水汽来源[19]。
图1 巴松措及派山谷地理位置[图1(b)和图1(c)中数字为测得的10Be年龄,M1、M2、M3分别代表派山谷保留的三组冰碛]Fig.1 Geographical location of the Basongcuo Catchment and Pai Valley[The numbers in Fig.1(b)and Fig.1(c)are the10Be ages,and M1,M2,and M3 represent the three groups of moraine of Pai Valley,respectively]
图2 林芝(a)、波密(b)、米林(c)三站月平均气温及月降水量Fig.2 Mean monthly air temperature and mean monthly precipitation values for Linzhi(a),Bomi(b)and Milin(c)
巴松措流域位于念青唐古拉山脉东段与喜马拉雅山脉交界处的南坡,南迦巴瓦峰西北方向约115 km处。根据GLIMS(Global Land Ice Measure⁃ments from Space)数据,该流域现代冰川覆盖范围约220.2 km2,其中冰川面积大于0.5 km2的有82条。巴松措湖沿岸有一组侧碛和终碛存在,也是该流域最外围的冰碛,代表了该流域最广泛的冰川活动。侧碛在湖的南岸保存较好,长19 km,北岸侧碛只有约2.4 km,与侧碛构成一体的终碛位于湖的出口处,由几组冰碛丘陵组成,表明冰川在该阶段到达巴松措现在的出口位置,距离主要支流(中措、白兰沟、新措)现代冰川末端分别为~33、~36和~48 km,其上分布有花岗岩漂砾。为了约束该次冰期的数值年代,Hu等[18]在最外围的丘状终碛顶部采集了5个漂砾样品,使用宇宙成因核素10Be暴露测年法测得了其暴露年龄,为了方便统一比较,本研究利用CRONUS-Earth 3.0计 算 程 序(http://hess.ess.washington.edu),Lal(1991)/Stone(2000)生产率随时间变化的模型[33-34]重新计算了该组冰碛的10Be暴露年龄:(20.1±1.4)ka、(19.1±1.4)ka、(20.1±1.4)ka、(26.4±1.7)ka和(20.5±1.4)ka(表1);另外,Hu等[18]对南北两岸的侧碛采集的23个沉积样品开展了光释光测年,结果显示冰进最大规模在16~30 ka之间,综合判定这组终碛-侧碛垄形成的年代是LGM时期。
派山谷发源于南迦巴瓦峰西坡,呈NW走向,山谷内分布有三组冰碛,最外层终碛位于该谷出口,末端海拔约2 970 m,距离现代冰川约7.1 km,表面较平坦,顶部有大量花岗岩漂砾,Hu等[19]采集了5块花岗岩漂砾样品,进行10Be暴露测年法测年,重新计算后结果分别为(33.1±2.3)ka、(23.7±1.6)ka、(22.7±1.4)ka、(22.6±1.6)ka和(22.7±1.5)ka(表1),经过Peirce检验后,33.1 ka被认为是异常值而剔除,判定此次冰进发生在距今~23 ka,即LGM时期。
表1 10Be年龄重新计算结果Table 1 The recaculation of10Be exposure dating results
本研究所使用的巴松措流域及派山谷的现代冰川范围界定主要依据全球冰川编目资料(Ran⁃dolph Glacier Inventory 6.0,RGI 6.0),可从GLIMS数据集(http://www.glims.org/)获得。冰底地形可以由现代冰川数字高程模型(Digital Elevation Mod⁃el,DEM)(http://www.gscloud.cn/)减去现代冰川厚度得到,其中现代冰川厚度数据在(https://www.research-collection.ethz.ch/)中获得。流域边界由LGM终碛垄外围以及流域的分水岭连接绘制,中流线根据冰川谷中等高线上凹凸特征及综合考虑3D地形及流域边界手工绘制。由国家气象科学数据中心(http://data.cma.cn/)下载得到1998—2019年周边8个气象站的日值气象数据,处理后计算得到年平均6—8月气温及年平均降水量数据。
本研究应用冰川的二维平衡剖面模型方法,沿中心流线重建冰厚。该模型不考虑冰川底部滑动(sliding),冰川流动靠冰川的变形(deformation)来驱动,并假设:冰川为完全塑性体,具有特定的屈服应力(τY),要使冰川发生变形运动,其受到的驱动应力(τD)应该达到屈服应力的这个门槛,冰的τD由冰的重量和表面梯度决定。如果驱动应力小于屈服应力,那么冰不会流动,而是通过表面变厚变陡,增加τD。另一方面,冰面形态会不断调整以保持τD和τY二者平衡的状态保持冰川的变形运动。这个状态可以表示为:
式中:ρ为冰密度(约为900 kg·m-3);g为重力加速度;H为冰厚;h为冰面高程;x为水平坐标(x轴与冰的流动方向平行,冰川上面的部分为正)。然后,从终碛向上沿冰川逐级迭代求解[35]:
式中:b为冰川底床高度;∆x为步长;τˉav为平均基底剪应力;F为形状因子;i为迭代次数。对于山谷冰川和其他受地形限制的冰川来说,侧向阻力也可以提供显著的流动阻力,可以将形状因子(F因子)加入公式,表示床体所支撑的驱动应力的比例:
式中:A为冰川截面积;p为截面积的周长;H为某一点的冰厚。对于冰盖或冰原等不被约束的冰体,F=1。
应用以上原理,“GlaRe”工具是在ArcGIS中使用的用Python编写的程序,根据有限的地貌信息,就可以快速实现上述方程的运算,生成古冰川3D表面[21]。模型至少需输入的参数有:冰底地形DEM、冰川流线矢量文件。具体操作步骤如下(图3):
图3 古冰川重建流程图Fig.3 Flow chart of paleo-glacier reconstruction
(1)在图新地球中绘制流线以及冰川流域边界的矢量文件,其中流线的绘制尽可能细致,有利于后期插值的可靠性。需要说明的是,“GlaRe”模块虽然可以根据终碛垄位置自动生成冰川流线,但对于形状复杂的冰川来说,自动生成的中流线容易出错[36],因此本研究选择了手动绘制的方法;
(2)使用“Construct Interval Nodes”工具,输入流线,自定义步长,输出流线的点矢量文件;
(3)通过“Define Shear Stress”定义基底剪应力,本研究将剪应力值分别设置为75 kPa、100 kPa和125 kPa,分析其对重建结果的影响;
(4)在“Flowline ice thickness tool”中输入无冰地形和流线的点矢量文件,得到流线逐点的冰厚度;
(5)在冰川流动受地形约束的地方,应用“Ffac⁃tor correction with user given cross-sections”,得到考虑F因子情况下的冰厚度点矢量文件;
(6)最后在“Glacier surface interpolation”中,结合冰川流域范围对冰流线厚度的点矢量文件应用“Topo to raster”插值方法进行插值,得到古冰川的三维表面DEM;
(7)古冰川厚度可用古冰川表面DEM减去无冰地形得到,冰储量可通过ArcGIS中的“Surface Volume”求得。
Pellitero等[21]对GlaRe模块进行了测试,结果显示对冰川储量、表面积和ELA的重建误差分别在25%、20%和75 m以内,考虑F因子后误差显著降低,分别降至10%、6%和10~15 m以内。
冰川平衡线指某一时段内冰川上物质平衡为零的所有点的连线,一般有年平衡线或多年平衡线,ELA即冰川表面平衡线所对应的海拔高度[37],ELA以上属于冰川的积累区,以下属于冰川消融区。本研究采用AAR、AABR两种方法计算ELA。AAR法的原理为:冰川处于稳定状态时,积累区面积占整个冰川面积的比值(AAR值)是固定的[38],这样就可以在GIS中通过冰川表面积计算冰川的积累区范围,从而确定ELA。全球范围来看,冰斗和山谷冰川的AAR值的范围在0.5~0.8之间[39],张鲜鹤等[40]研究发现念青唐古拉山脉东段的AAR值在0.5~0.7之间,而喜马拉雅山脉东段的AAR值都在0.5以下。AABR法利用冰川积累区面积以及消融区和积累区物质平衡梯度(冰川消融和积累随海拔高度的变化)的比值,即BR值来计算ELA。Rea[41]计算了全球具有物质平衡观测记录的冰川AABR值,认为全球平均值为1.75±0.71。由于夏季的高消融率,建议把AABR值为1.8~2.2作为中纬度海洋型冰川的代表[39,41-42]。本研究使用Pellitero等[42]编写的ArcGIS模块来实现ELA的重建计算,本文分别计算了不同AAR(0.5~0.7)和AABR(1.8~2.2)值的ELA值,通过两种方法对比得出冰川在LGM时的ELA。
鉴于冰川ELA对气候变化的敏感性,ELA的变化值(∆ELA)可以很好地反映气温和降水的变化情况[43]。Ohmura等[44]总结了平衡线处的气温降水关系,建立了平衡线处的气温、降水关系的经验回归模型(P-T模型)。施雅风等[45]对中国西部17条现代冰川的记录进行研究,得出了ELA处的夏季(6—8月)平均气温(T,单位:℃)和年降水量(P,单位:mm)的对数回归模型:
为了分析方便,另外考虑到ELA变化导致的气温变化,将P-T模型变形为ELA处气温变化与降水变化的关系式:
式中:∆T为夏季气温变化量(T1-T2);∆P为年降水量变化百分比(P1P2);T1、P1和T2、P2分别代表古今的夏季气温和年降水量。由于气温和降水是决定ELA的最重要的因素,本研究采用ELA与气候变化(气温和降水)的关系式来进一步约束古气候条件,即LR模型[46-48]:
式中:∆ELA、∆T、∆P分别为LGM以来的ELA的变化量、夏季平均气温变化量和年降水变化量,Ohmu⁃ra等[44]研究了全球70条冰川的平衡线上气温降水关系,得出范围在2.5×10-3~3.3×10-3℃·分别表示年降水量、夏季平均气温随海拔的变化情况。本研究假设这些参数在LGM时期与现在一致,上述两个经验模型,可以模拟出LGM时期的气温和降水的组合情况。
根据冰川地貌学证据,本文应用“GlaRe”模块模拟得到了LGM时期巴松措流域与南迦巴瓦峰附近的派山谷冰川规模的各项参数。Pellitero等[21]通过对现代冰川的重建进行模型验证,对比发现用Topo to Raster插值法和克里金插值法及IDW插值法对冰储量和冰川面积的模拟结果相差不大,且对于后期ELA的计算来说Topo to raster法相对更优秀;并且本研究通过实验发现,Topo to raster插值方法得到的冰川表面更加符合冰川的实际形态,因此此次重建均使用该方法进行插值。根据文献[49],大多选择100 kPa的基底剪应力,如此重建得到的巴松措流域的冰川冰储量达到了274.4 km3,面积达982.3 km2,约是现代冰川的4.5倍;派山谷在LGM时期冰储量达0.51 km3,面积达5.76 km2。而一般冰川的基底剪应力在50~150 kPa之间[49],考虑到冰底地形的起伏会导致基底剪应力变化,本研究分别选取75 kPa、100 kPa和125 kPa的基底剪应力,重建出两研究区末次冰盛期期间的冰川表面(表2~3)。在巴松措流域,增加或减少25 kPa的基底剪应力重建得到的末次冰盛期时的冰储量与100 kPa相比,得到的冰储量结果相差大概在30%以内,冰川面积的差值在15%以内。统一使用AABR法进行ELA的计算,得到75 kPa和125 kPa与100 kPa基底剪应力得到的ELA值分别为4 388~4 413 m、4 472~4 497 m和4 580~4 605 m,可以看出差异在84~108 m。派山谷冰川在末次冰盛期时,在100 kPa的基础上上下变化25 kPa,造成的重建得到的冰储量差异也在30%以内,冰川面积的差异在13%以内,AABR计算得到的三种基底剪应力对应的ELA分别为3 566~3 616 m、3 619~3 669 m、3 672~3 722 m,依次相差53 m。可以看出,对巴松措这种复杂的大流域冰川的重建,整体来说基底剪应力的差异会导致重建所得的较大的ELA值差异。由于巴松措湖(最大水深约120 m[51])和派山谷底碛的影响,重建时所输入的冰底地形在靠近冰川末端处会比LGM时期真实的冰底地形偏高,相应地会导致重建出的冰川在冰舌处厚度偏小。另外,Finlayson[52]评估了冰流模型在不同底碛及冰川消退后地形改变情形下对冰川厚度重建结果的影响,认为底碛及冰川消退后地形改变对该方法重建的古冰川表面高程整体上有不到1%的影响。将AAR和AABR法计算得到的ELA值取平均,得到LGM时期巴松措流域冰川和派山谷冰川的ELA分别为4 460~4 547 m和3 569~3 694 m。
表2 巴松措冰川流域不同基底剪应力的重建结果Table 2 Reconstruction of the LGM glaciers of BasongcuoCatchment under different basal shear stresses
图4 LGM时期巴松措流域[(a),(b)]和派山谷[(c),(d)]冰川表面和厚度重建结果Fig.4 The reconstruction of ice surface and thickness of Basongcuo Catchment[(a),(b)]and Pai Valley[(c),(d)]
表3 派山谷冰川不同基底剪应力的重建结果Table 3 Reconstruction of the LGM glaciers of Pai Valley under different basal shear stresses
通过使用AAR和AABR两种方法计算两地在末次冰盛期时的冰川平衡线高度(表4),可以看出对于两个流域来说AABR法得出的ELA值不确定范围更小,而且两种方法计算得到的∆ELA相差很小,而本研究后期对于气候的模拟主要输入的是∆ELA的值,因此本研究选择AABR法得到的∆ELA值代入模型进行温度和降水条件的计算。为了尽量减小系统误差,现代ELA也将冰面DEM输入模型进行计算,其中得到巴松措流域和派山谷的∆ELA分别为535 m和1 059 m。其中由于谷源部分地形陡峭,派山谷内已无现代冰川发育,仅有少量季节性积雪,因此其现代ELA是通过输入相邻冰川表面计算得到的。本研究重建得到的巴松措地区的∆ELA与Hu等[18]的估算结果小约155 m,由于本研究和Hu等[18]对于巴松措流域范围的划定有一定差别,可能是造成该差异的原因。而本文对派山谷地区∆ELA重建结果比Hu等[18]的估算结果大150 m左右,除了底碛的影响外,对于ELA计算方法的不同以及现代冰川ELA的确定也可能是造成该差异的原因。
表4 研究区LGM时期冰川平衡线高度重建结果(单位:m)Table 4 ELA reconstruction during LGM in the study area(unit:m)
利用研究区附近的八个气象站点的气象观测数据(表5),线性拟合得到该区域1998—2019年的平均夏季气温(T,℃)随海拔(z,m)的变化情况,即夏季气温递减率
表5 本研究所用气象站点气象数据(1998年1月—2019年12月)概要Table 5 Summary of the modern climate data(1998-01—2019-12)for meteorological stations used in the study
研究区的降水量(P,mm)随海拔(z,m)升高而降低[50],利用上述八个站点的气象数据计算得到1998—2019年平均的年降水量梯度a-1·km-1),可拟合成以下关系式:
本研究将∆P设定为现代降水量的30%~取值来自公式的值设定为(2.9±0.4)×10-3℃·mm-1,代入模型[式(6)],得到了LGM时期两地冰川平衡线处的气温和降水组合(图5)。可以看出,降水减少30%~70%的情况下,不考虑参数引起的不确定性,要达到LGM时期的冰川规模,巴松措流域和派山谷地区所需要的夏季平均气温分别比现在低2.14~2.59℃和4.08~4.78℃,可以看出气温的降低是LGM时期青藏高原东南地区发生冰进的主要驱动力。
仁错湖(30.73° N,96.68° E,海拔4 450 m,图1)的孢粉记录显示LGM时期青藏高原东南的地区的年降水量是现在的40%[54],在这样的降水条件下,LR模型模拟得到巴松措和派山谷地区LGM时期的气温分别比现在降低了2.48℃和4.61℃。LR模型的不确定性主要由参数和∆ELA决定,由模型构成可以看出,参数大小引起的模拟结果的不确定性主要与降水梯度和∆P有关,其范围为±0.05℃;由∆ELA的不确定性引起的结果的不确定性主要体现在派山谷地区,约为±0.08℃,主要与气温递减率和降水梯度有关,不随∆P变化而变化。
P-T模型模拟结果显示(图5),在LGM时期降水量是现在降水量40%的情况下,巴松措和派山谷地区的夏季气温分别较现代低4.41℃和6.51℃。P-T模型的不确定性主要由∆ELA的不确定性决定,主要体现在派山谷地区,且只与气温递减率有关,不确定性值为0.1℃。由于LGM以来研究区所在地块以至少5 mm·a-1速率遭受侵蚀[55],因此本研究使用120 m的剥蚀来修正ELA的计算,因此∆ELA结果可能被低估(不超过4%)。此外考虑到气温递减率(4.0℃·km-1)下,由于侵蚀所引起的温度变化的不确定性,最终降温大小大约被低估0.48℃。因此经过校正后,巴松措和派山谷LGM时期温度比现在分别降低了2.96~4.89℃和5.09~6.99℃。
图5 LR模型模拟结果(红色、黑色、蓝色分别代表的值为3.3×10-3℃·mm-1、2.9×10-3℃·mm-1、2.5×10-3℃·mm-1)Fig.5 Simulation results of LR model(The red lines were plotted by usingvalue of 3.3×10-3℃·mm-1,the black lines by 2.9×10-3℃·mm-1,and the blue lines by 2.5×10-3℃·mm-1)
以上通过LR和P-T模型计算得到了LGM时期两地可能的气候情景,这些气温和降水组合有助于评估其他气候代用指标的定量重建结果。由海登湖和仁错的高分辨率孢粉记录可知,青藏高原东南地区LGM期间与现在相比,7月气温降低了2~5℃,年降水量减少60%[54]。25次古气候模拟比较计划(Paleoclimate Modelling Intercomparison Project,PMIP)的平均值显示LGM时期研究区年平均气温约比现代低4.0~4.5℃,年降水量相对现代减少约15%[14],本研究重建的LGM时期的气温与孢粉记录和PMIP模拟结果相差不大。Xu等[56]通过重建念青唐古拉山脉西段东南坡LGM时期的冰川规模,得到当时比现代气温低2.9~4.6℃。Chen等[57]计算得到青藏高原东南帕隆藏布流域五条冰川LGM时期的∆ELA平均为917 m,并重建得到当时的夏季气温至少比现在低6.3℃。Zhou等[58]重建得到波堆藏布LGM时期的∆ELA约为600 m,重建得到LGM时期气温比现在低6.6℃左右,这与本研究所重建的派山谷地区的∆ELA和气温下降值很相似,但很明显,巴松措流域LGM时期的∆ELA、气温变化相比以上偏小,这可能是因为重建的两个流域的冰川发生在LGM时期内不同的阶段,但整体来看,两地LGM气温下降值的重建结果均处于以往研究所重建得到的气温下降值的范围内。
从10Be年龄来说,巴松措流域冰川约在19~20 ka达到现在湖出口的位置,而派山谷的冰碛10Be年代学结果显示冰川前进发生在末次冰盛期早期(~23 ka),恰好为古里雅冰芯连续记录显示的LGM最低温所出现的时段[53],两地区的气温和降水组合体现了青藏高原东南地区具体到LGM不同阶段的气候状况(图5~6)。另外,张廉卿等[59]对念青唐古拉山西段冰川的研究发现平衡线高度变化主要受气候因素的影响,两地均在印度季风控制区,相对来说派山谷更靠近高原边缘,相比巴松措受到印度季风的影响更大,降水更加充沛,平衡线高度上升快。本研究的结果也支持了前人的工作,即从高原边缘到内部,从东南到西北LGM期间的∆ELA逐渐减小[1]。除气候因素外,两地冰川的范围以及冰底地形的复杂性都有很大不同,尤其是研究区位于东喜马拉雅构造结,地形复杂陡峻,这些因素也会对平衡线高度变化量造成一定的影响。总的来说,气温降低是该地区LGM时期ELA下降的主要原因,即在气候干燥(冰川积累量减少)的情况下,足够的低温增加了积累区面积并减少了冰川的消融[60],使冰川净积累增加,从而导致该时期冰川发生大规模冰进。
图6 P-T模型模拟结果(图中红蓝色线为∆ELA不确定性造成的模拟结果的不确定性边界值)Fig.6 Simulation results of P-T model(The red and blue lines represent the uncertainty boundary values of the simulation results caused by the uncertainty ofΔELA,respectively)
本文根据冰川地貌证据,应用“GlaRe”模块,对LGM时期青藏高原东南的巴松措流域以及派山谷流域的冰川进行了重建,得到以下结论:
(1)100 kPa的基底剪应力下,重建得到末次冰盛期时巴松措流域的冰川面积达到982.3 km2,约是现代冰川的4.5倍,冰储量约为274.4 km3;派山谷无现代冰川分布,其LGM时期的冰川面积达5.76 km2,冰储量约为0.51 km3。由AAR法和AABR法计算得到的LGM时期两冰川的平衡线高度分别为4 460~4 547 m和3 569~3 694 m,相比现在分别降低了535 m和1 034~1 184 m。
(2)基于冰川ELA的变化,结合仁错湖花粉记录显示的降水信息,运用冰川气候模型,得到了LGM时期两流域冰川平衡线处的气温和降水组合。在降水减少60%的情况下,考虑到LGM以来的构造剥蚀对平衡线高度变化的影响,模拟得到LGM时期巴松措流域和派山谷冰川ELA处的夏季平均气温分别比现在低大约2.96~4.89℃和5.09~6.99℃。