中国极地冰盖数值模拟研究进展

赵励耘，唐学远，张通，冷伟，艾松涛，刘成彦，王玉哲，刘岩，岳超

（1.北京师范大学全球变化与地球系统科学研究院，北京 100875；2.自然资源部极地科学重点实验室中国极地研究中心，上海 200136；3.北京师范大学地表过程与资源生态国家重点实验室，北京 100875；4.中国科学院数学与系统科学研究院，北京 100190；5.武汉大学中国南极测绘研究中心，湖北 武汉 430079；6.南方海洋科学与工程广东省实验室（珠海），广东 珠海 519080；7.山东师范大学地理与环境学院，山东 济南 250014）

0 引言

极地冰盖（图1）对全球热量平衡和海平面上升具有重要的影响。IPCC第六次评估报告指出，气候变化评估中最大的不确定性之一来自对极地冰盖对未来海平面上升贡献量的估计，格陵兰冰盖物质损失对20世纪全球海平面变化贡献了29%，南极冰盖物质损失区域差异性大，其中西南极冰盖底部由于受到海洋环流影响融化速率相对较快［1-2］。此外，东南极尤其是威尔克斯盆地（Wilkes Basin）也发现潜在的不稳定性［3］。

图1 南极冰盖（a）；北极的格陵兰冰盖、斯瓦尔巴群岛和冰岛（b）［在图中标注了我国建立的台站（长城站、中山站、昆仑站、泰山站、罗斯海新站、黄河站）以及本文提到的一些研究区域或冰川位置（威尔克斯盆地、伊丽莎白公主地、埃默里冰架、普利兹湾、Totten冰川、雅各布港冰川、瓦特纳冰帽）。Dome A与昆仑站非常近，Austre Lovénbreen冰川和Pedersenbreen冰川与黄河站距离非常近，不再标注］Fig.1 Antarctic Ice Sheet（a）；Greenland Ice Sheet，Svalbard Islands and Iceland in the Arctic（b）［The stations established by China（Great Wall Station，Zhongshan Station，Kunlun Station，Taishan Station，Ross Sea New Station and Yellow River Station）and some study regions or glacier locations mentioned in this paper（Wilkes Basin，Princess Elizabeth Land，Amery ice shelf，Prydz Bay，Totten Glacier，Jakobshavn Glacier，Vatnajökull ice cap）are marked in the figure.Dome A，very close to Kunlun Station，Austre Lovénbreen Glacier and Pedersenbreen Glacier almost next to Yellow River Station，are not marked］

由于极地冰盖气候恶劣、难以接近，评估当前极地冰盖变化的主要手段是基于遥感监测，在部分区域结合野外实地测量。遥感技术可以对冰盖实现全覆盖、连续性监测。地面实地测量难以大范围实施，但是可以更精确地探测冰盖内部构造和底部地形等关键物理信息。遥感技术和地面实地测量等多源观测为研究极地冰盖变化提供了必备的数据资料（表1）。

表1 冰盖关键动力过程、模型中的相关变量和主要观测手段Table 1 Key dynamic processes of ice sheet，relevant variables in the model and main observational approaches

基于多源观测数据作为模式边界条件和外部驱动，可以通过数值求解冰流动力学模型，详细描述冰流的运动或预测冰盖的未来演化。冰盖数值模拟通常包括两种：一种是诊断模拟（diagnostic modeling），用以诊断冰盖在稳态假设和特定边界条件及参数下的动态平衡状态；也可以用来测试冰盖对某些参数或边界条件的敏感性以及微小变化对模型行为的影响，从而确定可能反馈机制的性质。诊断模型可用于提高对控制特定冰流行为的参数的理解，或研究一般冰盖中一个或多个物理参数与边界条件的重要性。另一种是预估模拟（prognostic modeling），用以描述冰盖运动的过程，通过调整模型并使其包括不同的物理过程，以便更好地将模拟预估与实际情况相匹配，从而更好地预测冰盖的未来演变。预估模型既可直接用于预估（或预测）冰盖的未来变化；也可用于重建冰盖的历史演化过程，从而有助于准确理解冰盖未来可能发生的变化。

如果要考虑到作用在冰上的所有应力分量，full-Stokes模型可能是最合适的工具。在某些情形下忽略部分应力分量，可以得到简化近似方程。基于full-Stokes模型或者简化近似的冰流模型，各国研究团队陆续研发了一些用于模拟冰盖状态或者预估其未来变化的数值求解方案，即冰流模式［4］。虽然我国在冰盖数值模拟方向的研究起步较晚，但是目前也形成了三个自主研发的冰流模式——多温型陆地冰模型PoLIM［5］、三维并行冰盖演化模型［6］和高精度并行有限元冰盖模型［7-9］。PoLIM为二维Blatter-Pattyn高阶近似热-力耦合冰流模型，以焓方法计算冰川热力场，包含冰下水文模型，采用有限差分方法求解计算。张怀等开发了三维并行冰盖演化模型，进行了格陵兰冰盖的模拟演化［6］。冷伟等研发的三维full-Stokes和高阶近似冰盖模式求解器，使用了适用于大规模数值模拟的快速求解方法［7-8］，并采用了保持局部质量守恒的数值离散方法以提高计算精度［9］，具有热-力耦合和处理触地线进退过程的功能［10-11］。

冰盖与大气、基岩、海洋接触，在边界处存在复杂的过程，包括冰盖表面物质平衡、冰盖水文过程、冰下水热环境、冰架崩解、冰架底部消融和再冻结等［12-13］（表1、图2）。这些过程都会对冰流系统的物质平衡和稳定性产生重要的影响，是最近冰盖动力学研究领域的热点和难点。目前在理论和观测上对上述关键动力过程和机制的认识都存在很大不足，严重制约了模式的发展和模拟结果的置信度。尽管国际冰盖模式发展时间较长，但是其对上述过程多采取简单的参数化处理或者尚未考虑，其模拟结果仍存在很大的不确定性［4］。特别的，近年来开展的冰盖-海洋耦合模式比较计划，不同耦合模式对海洋强迫的模拟结果差异很大。因此，模式的改进需要以冰盖典型区域长期持续的高精度精细测量作为基础支撑以提升对冰盖关键过程和机制的理解。

图2 冰盖关键动力过程示意图（本示意图的地形更适合西南极冰盖，而不适用于东南极和格陵兰冰盖大部分地区）Fig.2 Schematic diagram of key dynamic process of ice sheet（The topography in this diagram is more suitable for West Antarctic Ice Sheet than most region in East Antarctic and Greenland Ice Sheet）

前些年已有数篇针对极地冰盖数值模拟的评述文章［14-15］。近十余年来，我国冰川学家在极地冰盖数值模拟方面产出了不少新的研究成果。本文简要介绍冰盖数值模拟方法，进而梳理了我国在极地冰盖数值模拟方面的新进展，并对未来的研究前景进行了展望。

1 冰盖数值模拟方法简述

开发冰盖动力学数值或解析模型的主要用途是为了更好地了解冰盖的流动行为以及它们如何对外部环境的强迫做出响应。这些模型通常基于一些描述冰流的基本定律或假设，经过简化后可以得到解析或数值解。冰盖动力学模型包括动力学方程和初始条件以及边界条件。动力学方程主要包括冰盖的质量、动量、能量守恒定律和本构方程Glen流动定律。目前最为完整的动力学方程是热-力耦合的三维full-Stokes方程，它将描述质量、动量守恒的Stokes方程和描述能量守恒的热传导方程进行耦合，适用于各种情形的冰盖模拟。在Stokes方程的基础上进行不同程度的简化近似，可以得到静水压力近似（hydrostatic approximation）、一阶近似（first order approximation）、浅冰盖近似（shallow ice approximation）、浅冰架近似（shallow shelf approximation）等［16-17］。在full-Stokes方程中忽略垂直切应力τzx和τzy沿着运动方向的导数，可以得到静水压力近似。在静水压力近似基础上，忽略速度竖直分量的水平导数，即假设速度竖直分量vz仅依赖于竖直坐标z和时间t，vz=vz（z，t），得到一阶近似。在一阶近似基础上进一步简化，把冰盖看作是平行底面的剪切流，水平切应力τxz和τyz对水平方向的动量平衡起到重要作用，忽略法向应力偏量和竖直切应力τxy，可以得到浅冰盖近似模型；而如果应用在冰架上，假设冰架底面切应力为零，流速的水平分量不依赖于深度，即vx=vx（x，y，t），vy=vy（x，y，t），则得到浅冰架近似模型。浅冰盖近似和浅冰架近似可以结合成混合模型（hybrid model），应用到整个冰盖的模拟［18-20］。

黏性系数是动力学方程的关键参数，直接影响着冰流的速度大小，并且敏感地依赖于冰温和冰晶组构等因素。因此，热-力耦合的动力学方程的模拟结果相对更接近实际情形。full-Stokes方程及其近似方程是建立在冰各向同性的假设之上，由于实际冰盖冰晶组构各向异性等原因，冰流动力学方程中的重要参数（例如冰的黏性增强因子、冰底摩擦系数）可以通过求解反问题来反演，使得表面冰流速的模拟值与观测值误差最小［21-22］。原则上，对于full-Stokes冰盖模型而言，随着计算机技术的发展使得完整应力分布的充分计算已成为可能并得到逐步实现。但对于百年以上的积分，某些大陆尺度的冰盖演化模拟在数值计算上仍不可行。

技术上，冰盖动力学模型求解的空间域是基于地形数据建立的研究对象的数字高程模型。基于冰厚、海平面高程和冰底高程数据，可以判断出冰盖触地部分和冰架部分。在研究区域的边界（冰面、冰底、冰崖、侧面等）要给出合理的边界条件，主要包括：（1）冰面温度、表面物质平衡；（2）冰底触地部分的热力平衡方程（包括地热通量等）、运动状态等；（3）冰盖底部的压融点温度和物质平衡等；（4）冰崖处的海水压强和崩解过程等。其中，表面物质平衡和冰架底面物质平衡体现的是大气和海洋对冰盖作用的强迫场。冰盖表面物质平衡的获取途径较多，包括实地测量、再分析资料和区域气候模式的模拟结果等。冰架底面物质平衡的直接监测条件受限（主要是回声测深仪［23］和钻孔），缺乏基础数据。以往和目前的研究多采用参数化方案或者海洋模式的模拟结果进行估算。如何给出较为精确的冰架底面物质平衡是目前研究的主要难点之一。另外，冰盖底部的运动状态很难监测，而且对冰下水热环境变化非常敏感等，具有较大的不确定性。以往研究中的常用方法是在假定的滑动定律形式下通过数据同化反演底部滑动系数［24-25］。冰崖处的崩解过程对南极冰盖的稳定性非常重要，也是南极冰盖物质损失很大的不确定性来源。但是崩解过程如何在冰盖模型中精确地刻画是迄今尚未较好解决的难题。在以往的冰盖模拟研究中常见的处理方式有：冰崖的位置保持固定，或者利用冰裂隙深度、冰面径流、应变率、应力、冰厚等对崩解速率进行参数化［4，26-28］。

冰盖模型的初始条件是确定起始时间冰盖的状态。如何得到符合冰盖实际情况的初始状态，即模式的初始化，是一个研究难题。前两年国际专门针对模式初始化开展了模式比较计划initMIP［29-30］。常见的方法有长期的预热模拟（spin-up）和基于数据同化方法的稳态模拟，即融合观测数据和数值模拟得到同化的模式参数和较为接近实际情形的初始状态。同化的参数在预估冰盖未来演化状态中是否仍然适用，是有待商榷的问题。

根据基本物理定律而构造的用来模拟冰盖状态或者预估其未来变化的一套流体力学和热力学偏微分方程组及其求解方案，称之为冰流模式。针对冰盖研发的冰流模式也称为冰盖模式。冰盖模式通过数值求解带有初始条件和边界条件的物理模型，模拟冰盖的速度场、温度场、应力场等主要物理量的时空分布以及冰盖形态和触地线位置的变化，从而有效地理解冰盖动力演化过程。冰盖模式计算的大致流程见图3。

图3 冰盖动力学数值模式计算流程图Fig.3 Flow chart of numerical simulation of ice sheet model

2 中国在南极冰盖数值模拟领域的主要进展

我国早期开展的南极冰盖数值模拟研究主要是一些理想化的冰盖诊断模拟试验，探索冰盖运动变化特征。唐学远等［31］利用基于浅冰近似模型的三维冰盖模式GLIMMER，计算了方形冰盖积分区域的流场特征量，研究了冰盖在大尺度长时间序列条件下对气候变化的反馈，考察了稳定态下冰盖演化的周期性行为。唐学远等［32］采用三维冰盖模式GLIMMER，对南极冰盖开展了理想化冰盖数值试验，给出了冰盖处于稳定状态下的各种演化特征曲线，并将该模式从三维情形简化为二维情形，在耦合温度场条件下模拟了二维冰流的演化，探讨了简化后的二维模型在模拟Dome A地区年代-深度-古积累率关系的应用前景。Zhang等［11］比较了三维full-Stokes冰盖模式求解器FELIX-S［7-9］和冰流模式Elmer/Ice两个不同Stokes冰盖模式对简化的海洋性冰盖的模拟结果，研究表明：对于诊断模拟，使用相同的几何形状和网格时，两个模型给出了相似的结果（差异＜2%），原因是两个模型的结构是相似的；对于预测模拟，发现FELIX-S（Elmer/Ice）接地线相对更后退（前进），结果与在诊断试验结果中观察到的微小差异一致，这是由于两个模型中底部边界条件的不同选择造成的［11］。

在过去四十年，我国南极科学考察以中山站、昆仑站、泰山站等为支撑平台，在东南极伊丽莎白公主地和Dome A地区逐步实施了冰芯钻探和强化观测［33-34］。沿中山站—Dome A断面，已逐步建设了星-空-地一体化的业务监测空间网络。基于这些极地考察区域的观测优势，我国在南极内陆冰盖典型区域开展数值模拟研究，取得了一些国际水平的研究成果。

Dome A在昆仑站附近，位于东南极冰盖分冰岭的中心，被认为是获取古老冰芯记录的理想地点［35］。Sun等［36］利用冰流模式Elmer/Ice，基于各向同性和简单的各向异性的冰晶组构的假设条件，对Dome A深冰芯所在区域冰温和冰龄的空间分布进行诊断模拟，模拟结果表明昆仑站可能提供过去60～70万年垂直钻探的高分辨率记录。Zhao等［37］采用冰流模式Elmer/Ice，将Wang等［38］发现的各向异性冰晶组构随深度交替现象进行参数化，并利用雷达等时层的年代信息作为约束条件，重新模拟了Dome A地区的冰温和冰龄空间分布，估算了昆仑站冰芯位置的底部融化速率，指出超过一百万年的古老冰有望在距离昆仑站5～6 km以内、冰底以上200 m深度处找到。Wang等［39］利用内部等时层和Vostok冰芯年代资料，得到Dome A六条等时层位上的深度-年代关系，并将其输入到Dansgaard-Johnsen冰流模式，揭示了Dome A地区和昆仑站深冰钻探点在16万年以来6个不同历史时期的古积累率的时空变化。张良甫等［40］采用基于full-Stokes方程的三维热力-动力耦合冰流模式Elmer/Ice，模拟了甘布尔采夫山脉的速度场和温度场。泰山站是中国于2014年新建立的科考站。Tang等［41］对泰山站的冰川与气象条件开展了综合的分析，并采用一维垂向模型模拟了冰龄和冰底温度。

在过去十几年，我国对东南极的埃默里冰架（Amery Ice Shelf）和普里兹湾（Prydz Bay）开展了长期的观测，对以这两个区域为代表的南极冰盖边缘开展了数值模拟研究。基于高分辨率的海洋-海冰-冰架耦合数值模式，Liu等［42］首次命名了普里兹湾东岸流（PBECC），继而诊断了PBECC的海洋热输送能力及其对Amery冰架底部质量平衡的影响［43］。Cheng等［44］利用自主研发的冰架海洋边界层模式对Amery冰架底部的边界层冰晶羽流过程开展了深入的研究和数值模拟，提高了我国自主研发冰架海洋边界层模式的技术水平。季青原等［45］利用PISM冰盖模式模拟了Amery冰架的速度场，模拟结果总体上与观测结果较为接近。Li等［46］利用Úa冰流模式模拟和分析了Amery冰架在三次崩解事件前后的变化，发现崩解事件对Amery冰架的稳定性影响很小；也对Amery冰架开展了反演和预测试验，模拟了冰架解体极端情景下的冰流加速状况［47］。Sun等［48］基于冰盖模式BISICLES，采取多种冰架底面融化率参数化方案，对南极包含Totten冰川的奥罗拉（Aurora）流域开展动力学模拟和预估研究。

3 中国在北极的冰盖数值模拟领域主要进展

我国在北极地区的冰盖模拟研究主要集中在格陵兰冰盖和斯瓦尔巴群岛、冰岛等地区的冰川。在针对整个格陵兰冰盖的数值模拟方面，Yan等［49］利用冰盖模式SICOPOLIS，探讨了不同现代气候强迫和初始化方案对于现代格陵兰冰盖模拟结果的影响，预估了21世纪格陵兰冰盖物质平衡的变化以及冰盖融化对全球海平面上升的贡献。

表面融水和溢出冰川动力过程的物质输送是造成格陵兰冰盖物质损失的两大主要因素。格陵兰冰盖表面消融产生大量融水，形成冰面水系并输送到冰裂隙、竖井、冰面湖，下渗到冰盖内部和底部，通过促进底部滑动和提升冰温，对冰盖运动速率与稳定性产生影响。格陵兰冰盖溢出冰川末端大多延伸到狭长的峡湾，通过峡湾连通到开阔的海域。海洋环流的变化带动峡湾水体的变化，对冰川的物质平衡和动力过程起到重要的影响。运用数值模拟预估格陵兰冰盖表面径流和溢出冰川的演化是理解和评估格陵兰冰盖未来物质平衡变化的重要方面。

在针对格陵兰冰盖表面径流的模拟研究中，我国在不同的时空尺度上开展了研究。Smith和Yang等［50］利用卫星遥感与无人机影像，提取了一个典型集水区的冰面水系，结合实地测量的冰面河出口位置径流量，提出了一种测算冰面径流量的新方法，研究结果表明，目前常用的区域气候模式过高地估算了冰面实际径流量，冰面融水输送的方式和效率对于冰盖物质平衡以及动力过程具有重要影响。Moore等［51］利用冰盖表面能量-物质平衡模型，研究平流层注入硫酸盐气溶胶的地球工程试验G4对格陵兰冰盖表面径流的影响，模拟结果显示在G4情景下冰盖表面径流比RCP4.5情景下减少了20%。Yue等［52］利用度日模型，开展了21世纪和22世纪不同气候情景下格陵兰冰盖表面径流的模拟和预估研究。

在溢出冰川动力过程的数值模拟方面，我国主要对格陵兰冰盖最快的溢出冰川——雅各布港冰川进行了数值模拟。Guo和Zhao等［53］采用冰盖模式BISICLES和同化的参数方案，重现了该冰川在2004—2013年间由于海水增暖造成的冰流速度加速等行为过程，在季节尺度上模拟的冰流速度、冰架末端位置变化、触地线位置变化与观测值具有很好的一致性，并且在RCP4.5情景下预估了冰川的动力学演化。

除了格陵兰冰盖，我国科研工作者依托黄河站为支撑平台，对北极斯瓦尔巴群岛的冰川等开展了实地观测和数值模拟有效结合的研究。Ai等［54］利用Elmer/Ice冰流模式对北极黄河站附近的Austre Lovénbreen冰川开展了数值模拟，发现该冰川除了以前实测发现的快速冰流区，还有另一个速度更快的区域。基于该模拟发现，2016年起开展了高精度GPS加密测量，证实了该发现，新发现的快速冰流区比以前测得的快速流速区流速要快8%左右。此外，以历史观测记录和Elmer/Ice冰流模式为参照，艾松涛等［55］在忽略动力学因素的条件下，利用Arc-GIS软件在冰川轮廓内的每个像素处计算冰面高程变化，进而模拟山地冰川面积、冰厚和体积的长期变化，并探讨了该方法的可行性和适用条件。在2015/2016年的北极异常升温背景下，黄河站Austre Lovénbreen冰川和Pedersenbreen冰川的流速都有较大增加。Wang等［56］通过模拟发现，底部滑动是冰流增速的主要影响因素，直接原因是冰下水量（subglacial water）的增加。随着冰川融水渗透到冰川内部，还可能出现冬季缓慢释放至冰川末端形成冰湖的现象［57］。通过冰川生命周期的模拟，斯瓦尔巴群岛上类似Austre Lovénbreen冰川的小冰川很可能在一百余年后彻底消失［58］。另外，Yue等［59］应用冰盖模式PISM，模拟和预估了冰岛冰帽在21世纪物质平衡和动力过程的演化。

4 总结与展望

回顾过去十几年，我国在极地冰盖数值模拟领域的研究取得了明显的进步。我国已经自主研发了三个冰盖模式，并在持续改进。同时，在极地冰盖及其典型区域上应用冰盖模式模拟，得到了一些重要的科学认识，取得了具备国际水平的研究成果。经过数十年的极地科考活动，我国在极地冰盖的监测方面形成了强化监测优势区域，其中包括南极的Dome A、Amery冰架、伊丽莎白公主地和北极的黄河站等，积累了大量宝贵的数据资源。近年来在格陵兰冰盖典型区域也开展了空-天-地一体化监测［60］。这些都为研究冰盖变化关键过程提供了良好的数据平台。

与此同时，我国在极地冰盖数值模拟领域也面临着一些问题，制约着当前和未来的工作，其中包括：研究人员和人才很少，没有形成实力较强的团队；拥有自主知识产权的冰盖模式在研发上处于起步阶段，需要持续投入；我国在国际冰盖模式比较计划中的参与度不足；研究成果的数量和质量尽管有很大进步，但是涉及入海冰川或冰架崩解和消融的模拟研究较少；国际上已有研究团队探索了大气环流模式与冰盖模式的耦合［61-64］以及海洋模式与冰盖模式的耦合［65-67］，提高模式模拟的精度，而我国目前鲜有海-冰-气耦合模式方面的研究。

针对这些问题，需要提高对冰盖模式模拟研究领域的重视程度，持续加大投入力度，加强国内研究院所的交流合作，有效结合实地观测、遥感监测和数值模拟，形成研究合力；加强与大气、海洋等其他圈层领域互动，推进海-冰-气多圈层耦合和相互作用研究，提升极地冰盖模拟在气候变化科学领域的重要性；在我国传统观测优势区域，通过观测、模拟等综合研究手段开展联合攻关，有望在冰盖关键动力过程和机制的科学认识上有所突破，提升冰盖模式对物理过程的刻画能力，产出具有中国特色的世界先进水平的研究成果。展望未来，可能需要在以下几个方面加强研究力量：

（1）充分利用多源数据和区域气候模式，构建长时间高分辨率时空连续的冰盖表面物质平衡，准确刻画区域尺度空间变化特征；开展格陵兰冰盖表面融水和冰面水系的观测和模拟，关联冰面-冰下融水，探索冰面融水对冰盖运动的影响。

（2）开展冰盖典型流域的动力学模拟，研究海洋强迫对入海冰川末端或者冰架的作用机理，改进冰架消融和崩解在冰盖模式中的实现，探索入海冰川或者冰架与海洋之间的相互作用和耦合模拟，加深对关键动力学机制的理解。

（3）研究冰盖冰下水系分布、冰下水文过程及其对冰盖不稳定性的影响，开展冰下水文过程模拟，在冰盖模式中考虑冰下水文过程的影响。