冰川冻土, 2020, 42(1): 254-264 doi: DOI:10.7522/j.issn.1000-0240.2019.0092

冰冻圈技术

新一代冰流模式乌阿及其在南极埃默里冰架的应用

李腾,1,5, 陈卓奇2,3,5, 李慧林4, 程晓,2,3,5, 韦屹1,5, 刘岩1,3,5

1.北京师范大学 全球变化与地球系统科学研究院 遥感科学国家重点试验室，北京 100875

2.中山大学 测绘科学与技术学院，广东 珠海 519000

3.南方海洋科学与工程广东省试验室（珠海），广东 珠海 519000

4.中国科学院 西北生态环境 资源研究院 天山冰川观测试验站，甘肃 兰州 730000

5.中国高校极地联合研究中心，北京 100875

A new-generation ice flow model Úa and its application in Amery Ice Shelf， Antarctica

LI Teng,1,5, CHEN Zhuoqi2,3,5, LI Huilin4, CHENG Xiao,2,3,5, WEI Yi1,5, LIU Yan1,3,5

1.State Key Laboratory of Remote Sensing Science，College of Global Change and Earth System Science （GCESS），Beijing Normal University，Beijing 100875，China

2.School of Geospatial Engineering and Science，Sun Yat-Sen University，Zhuhai 519000，Guangdong，China

3.Southern Marine Science and Engineering Guangdong Laboratory，Zhuhai 519000，Guangdong，China

4.Tianshan Glaciological Station，Northwest Institute of Eco-Environment and Resources，Chinese Academy of Sciences，Lanzhou 730000，China

5.University Corporation for Polar Research （UCPR），Beijing 100875，China

通讯作者: 程晓， 教授， 从事极地遥感监测理论、 方法与应用研究. E-mail： chengxiao9@mail.sysu.edu.cn.

编委: 庞瑜

收稿日期: 2020-03-17   修回日期: 2020-05-13  

基金资助: 国家自然科学基金.  41830536.  41771081
国家重点研发计划项目.  2018YFC1406100
中国科学院前沿科学研究重点计划项目.  QYZDB-SSW-SYS024

Received: 2020-03-17   Revised: 2020-05-13  

作者简介 About authors

李腾（1992-），男，江苏丰县人，2015年在武汉大学获学士学位，现为北京师范大学在读博士研究生，从事南极冰盖模拟研究.E-mail：litengbnu@foxmail.com , E-mail：litengbnu@foxmail.com

摘要

南极冰盖不仅是全球环境变化的指示器， 其消融所产生成的淡水输入也是未来海平面上升的主要不确定性来源。数值模式是诊断冰流动力机制、 评估冰盖物质损耗的重要手段。本文首先介绍了乌阿（冰岛语Úa或英语Ua）冰流模式的基本原理， 并利用该模式模拟东南极埃默里冰架的动态变化。乌阿冰流模式基于质量和动量守恒方程的垂直积分， 在自适应不规则三角网格上求解微分方程， 仅用少数参数规则即可构造适应冰流动力特征的网格结构， 有效缩减运算时间。采用当前主流的模式边界数据集， 针对埃默里冰架设计了两个试验。试验一为反演试验， 试验中模式的代价函数在100次迭代后下降三个数量级， 表明模拟的流速与遥感观测吻合（RMSE = 13.35 m·a-1）， 但高频细节仍有待提高； 试验二为预测试验， 测试了模拟冰厚变化率的不确定性， 以自由漂移量接近零为标准选出一组最优模型参数， 最后假设埃默里冰架解体情景开展模拟， 结果表明冰架解体会导致海平面上升（45.36 ± 0.08） mm。随着资料更新迭代， 基于最新发布的南极底部地形数据模拟效果是否提升还有待未来检验。

关键词： 冰川学 ; 埃默里冰架 ; 南极冰盖 ; 数值模拟 ; 参数反演

Abstract

Antarctic Ice Sheet is not only the sensitive indicator of global environmental changes， but its melting is also the largest uncertainty of the future sea-level rise. The numerical model is an important methodology for diagnostic ice flow mechanism and assessing ice mass loss. This paper introduces the basic principles of a new generation ice flow dynamic model， Úa， based on which the two experiments are implemented to simulate the dynamics of Amery Ice Shelf （AIS） in East Antarctic. Leveraging the adaptive triangular finite-element meshes， such model solves vertically integrated formulations of the mass and momentum conservation equations. The generated mesh accurately captures the glacial dynamics from only a few parametric criteria， which reduces the computation time. On the basis of mainstream datasets， two experiments are designed in the AIS area. In the first part of the inverse run， the value of cost function in the model reduces for three orders of magnitude. Although the fine-scale structures remain to be improved， the output flow pattern agrees well with the satellite observation （RMSE = 13.35 m·a-1）. In the second part of the prognostic run， we test the sensitivity of ice thickness change rate， from which the optimal one with the minimum bias is selected to simulate the shelf’s collapse scenario. The results demonstrate the upstream ice system would contribute about （45.36 ± 0.08） mm to sea-level rise. Since Antarctic topography data were updated， whether our model can be improved from the new data deserves further research.

Keywords： glaciology ; Amery Ice Shelf ; Antarctica Ice Sheet ; numerical model ; parameter inverse

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

本文引用格式

李腾, 陈卓奇, 李慧林, 程晓, 韦屹, 刘岩. 新一代冰流模式乌阿及其在南极埃默里冰架的应用[J]. 冰川冻土, 2020, 42(1): 254-264 doi:DOI:10.7522/j.issn.1000-0240.2019.0092

LI Teng, CHEN Zhuoqi, LI Huilin, CHENG Xiao, WEI Yi, LIU Yan. A new-generation ice flow model Úa and its application in Amery Ice Shelf， Antarctica[J]. Journal of Glaciology and Geocryology, 2020, 42(1): 254-264 doi:DOI:10.7522/j.issn.1000-0240.2019.0092

①尽管这是通行和权威的译法， 作者认为“薄冰层近似”不仅贴切原意且更加符合中文的表达习惯。

0 引言

南极冰盖是地球上最大的淡水存储单元， 若南极冰盖全部融化将导致全球海平面上升57 m［1］。南极冰盖在重力的作用下由内陆向四周缓慢流动， 并在沿海受到地形的约束而形成漂浮在海洋上的冰架［2］。冰架表面吸收太阳辐射而产生融化， 底部因与海水接触而受到大洋环流的直接影响。当表面或底部的裂隙贯穿全部冰架时， 则会崩解产生冰山， 同时对上游冰盖的支撑作用（buttressing effect）减弱， 这是南极冰架调节冰盖稳定性的重要机制［3］。

当前全球气候变化已成为人类社会可持续发展面临的严峻挑战， 南极冰盖不仅与全球气候系统密切相关， 其物质损失是海平面上升预测当中最大的不确定性来源［4］。冰流模式是我们理解和刻画冰川运动机制及其对气候系统反馈效应的重要工具。从最简单的零阶“浅冰近似”（Shallow Ice Approximation， SIA）①到全应力分量的斯托克斯方程组（Full-Stokes Equations， FS）， 目前国际冰川学界已经开发出了数十种复杂程度各异的冰流物理模式。科研人员基于上述模式组织了多次比较计划， 用于验证数值方法的有效性。其中影响力比较大的模式比较计划包括海平面对冰盖演化的响应（Sea-level Response to Ice Sheet Evolution， SeaRISE）［5］， 三维海洋性冰盖模拟比较计划（3D Marine Ice Sheet Model Intercomparison Project， MISMIP3d）［6］等。

我国的冰流模式研究工作起步于山地冰川研究， 李慧林等［7］较早探讨了动力学模式对中国西部冰川预测的适用性， 并以天山乌鲁木齐河源1号冰川为例进行数值试验［8］。近年来， 国内一些学者也陆续将国际上先进的极地冰盖模式引入。唐学远等［9］较为系统地总结了南极冰盖数值模拟研究的前沿进展， 并指出了触地线等可能引起数值不稳定的问题， 随后将二维简化的 GLIMMER模式应用于南极的“冰期-间冰期”循环， 结果表明冰盖总物质量的波动符合古气候记录， 该模型能够刻画长时间尺度上的气候变化反馈机制［10］。季青原等［11］将并行冰盖模式（Parallel Ice Sheet Model， PISM）引入南极冰盖， 他们在埃默里冰架的试验结果表明模拟速度场能够重现冰架的流动过程， 只有在冰架前缘出现微弱差异。张良甫等［12］配置了更复杂的三维热动力耦合的全斯托克斯模式Elmer/Ice， 并尝试对兰伯特冰川上游至冰穹A的冰盖断面建模。针对冰架-海洋相互作用， Williams等［13］将三维热动力耦合模型应用于埃默里冰架地区， 发现冰腔内部由密度梯度驱动的局部环流显著影响冰海界面的温盐交换过程。Gong等［14］基于局部大气模式输出的结果作为强迫场， 模拟了气候变化情境下埃默里冰架动力学响应， 结果表明海水升温引发的底部融化可能驱动触地线退缩40 km。

乌阿冰流模式是近年来国际上新发展的一种有限元数值模式， 作为新一代冰流模式的优秀代表， 它以架构稳定、 配置简单、 计算效率高、 网格设置灵活、 支持并行运算等优势受到科研人员的青睐［15］。以最近完成的第二次线性响应模式比较计划（Linear Response Model Intercomparison Project， LARMIP-2）为例［16］， 相对于传统的粗分辨率SIA冰盖模式， 新一代冰盖模型普遍综合考虑了纵向薄膜应力和底部剪切应力， 能够更准确地抓住快速冰流的动力特征和适应触地线附件的应力突变； 在参与该计划的16组模型中， 仅有其中三组（含乌阿）在触地线附近实现了千米级尺度的精细模拟， 而以PISM和GRISLI为代表的固定网格模型典型尺寸都超过了10 km， 充分展现出自适应网格剖分方案的优势。王澄海等［17］系统总结了冰川数值建模的研究进展并将冰川模式发展归纳为四个阶段， 我国目前总体上属于第二阶段。以乌阿模式为代表的国际先进水平大多处于第三或第四阶段， 故将其介绍给国内科研人员将对促进我国冰川模式的发展有借鉴意义。本文以东南极兰伯特冰川-埃默里冰架系统为例， 设计了反演对比和前向预测两个试验， 定量评估未来100年埃默里冰架解体对海平面上升的贡献。

1 研究区介绍

埃默里冰架（Amery Ice Shelf， AIS）是东南极最大的冰流系统［图1（a）］， 面积约60 000 km2， 上游补给区分别是兰伯特冰川、 米勒冰川和费舍尔冰川， 它们由内陆的冰穹A绵延约1 500 km注入南大洋普里兹湾海域［18］。该流域内地形复杂， 气候特征差异大， 埃默里冰架的触地线最深处超过2 000 m， 冰架底部与海洋相互作用活跃［19］。由于澳大利亚和中国的考察站［见图1（a）中泰山站和中山站的位置］都位于该冰架附近， 两国都开展了系统的实地调查， 积累了宝贵的研究资料［20］。该冰架前端发育的三条裂缝吸引了众多冰川学家的注意［21］， 2019年9月底崩解的D28平顶冰山面积1 600 km2， 约占其总体冰量的1.1%［22］。根据Zhou等［23］最新估计， 目前埃默里冰架的物质状态整体处于微弱正平衡［（3.1 ± 9.4） Gt·a-1］当中。

图1

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

图1   埃默里冰架

黄框中为模拟计算的网格范围， 蓝线为模拟触地线， 红色箭头表示狄利克雷边界条件

Fig.1   The Amery Ice Shelf （AIS）

Landsat Image Mosaic of Antarctica （LIMA） map of the research area， and the location of the AIS in Antarctica （a）. The yellow polygon represents the computational domain； the unstructured triangular finite element mesh adapted in this paper （b）.The blue line represents the modeled grounding line and the red arrows represent the Dirichlet boundary conditions


2 数据与方法

2.1　模型原理

乌阿冰流模式由Hilmar Gudmundsson教授在英国南极调查局和诺森比亚大学的资助下开发维护， 最初用来研究触地线的稳定性［15， 24］。2018年公开共享全部源码， 全世界的冰川学者均可免费地下载、 修改和使用（源码仓储：https：//github.com/GHilmarG/UaSource）。该模式基于MATLAB计算软件（2017a或更高版本）， 相比于PISM和Elmer/Ice， 乌阿冰流模式的安装部署和编码计算都更为简单， 现已经广泛应用于多种类型自然冰流的模拟。在冰流控制方程方面， 乌阿冰流模式对冰流的三维斯托克斯方程全应力分量简化为基于垂向积分的二维动量方程， 简化给出了两种近似方案， 分别为“浅冰近似”和“浅冰流近似”（Shallow Ice-Stream Approximation， SSA）［25］， 见式（1）：

∇xy∙hD- tb=ρgh∇xys+12gh2∇xyρ

式中：∇xy为空间梯度算子； h为冰厚， 单位m； s为表面高程， 单位m； D为二阶应力张量， 定义为二维偏应力τij的矩阵， 单位kPa：

2τxx+τyyτxyτxx2τyy+τxx

tb为底部阻力， 单位kPa； ρ为冰的密度， 单位kg·m-3； g为重力加速度， 单位m·s-2。

在冰流的本构关系及与底部基岩的摩擦关系方面， 乌阿冰流模式引入了冰川学界广泛使用的Glen流动定律（式2）和Weertman滑动定律（式3）， 用户可以自行调节非线性方程的系数：

ϵ˙ij=Aτn-1τij

式中： ϵ˙ij为应变率， 单位a-1； τ为偏应力张量的二阶不变量， 单位kPa， 定义为τ =∑i, jτijτij/2； A是与冰温有关的各项同性黏度因子， 也是模型中的待定参数； n为自定义的非线性系数， 一般选3， 此时A的单位为a-1·kPa-3。

vb=𝒢- mCtbm

式中： vb 为底部滑动速度， 单位m·a-1， 对“浅冰流近似”而言底部滑动速度等于表面速度； C为需要反演的底部滑动系数； 𝒢为无量纲接地掩膜， 接地冰盖设为1， 漂浮冰架设为0， 亦即漂浮冰架的底部不存在剪切力也无须反演滑动系数； m为自定义的非线性系数， 一般选3， 此时C的单位为m·a-1·kPa-3。

在数值求解上， 乌阿冰流模式基于自适应不规则有限元三角网格［26］， 网格可以结合多种规则灵活调节尺寸， 比如表面坡度、 冰厚梯度、 触地线的距离等， 力求能够在模拟精度和运算效率两方面达到平衡。在其提供的多种网格控制选项中， 最常用的是依据应变率自适应缩放网格， 以及指定触地线空间范围与其对应的尺寸［图1（b）］。对于模拟山地冰川等可变计算区域的情景， 相应的网格可自动激活或失效， 计算边界亦可随之迁移。

在正演情景下， 模式联立式（1）、 （2）和（3）， 并基于全隐式或半隐式的Newton-Raphson中心差分算法迭代求解水平二维瞬时速度（U， V）和垂向冰厚变化率（h˙）。在满足Courant-Friedrichs-Lewy收敛条件的情况下， 为了加速运算， 模式可以自适应调节前向积分的时间步长， 这对冰厚演化的模拟至关重要。在反演情景下， 模型假设已知冰流初始速度， 使用基于贝叶斯理论的切线性伴随同化算法， 在计算空间内拟合最优模型参数黏度因子A和滑动系数C。反演中用于极小化计算的目标代价函数（objective cost function）形式为：

J = I + R

式中： I为失配函数（misfit function）， 用于量化模拟与观测之间的误差， 随着迭代而逐渐缩小； R为正则函数（regularization function）， 用于量化参数与先验值之间的差异， 有利于反演场抵抗高频噪声， 在二维网格面实现一定程度的空间平滑。失配函数的具体形式为：

    I=12B∫U-UmeasUerrordB+12B∫(V-VmeasVerror)dB+12B(h˙-h˙meash˙error)dB

式中：U、V、 h˙分别为水平速度在X方向的分量、 水平速度在Y方向的分量和冰厚变化率， 单位均为m·a-1， dB表示在整个计算区域上的面微分算子。

正则函数的具体形式为：

R=12B∫γs2∇p-p˜2+γa2p-p˜2dB

式中：γs 和γa 分别是调节参数变率和幅度的无量纲系数； p为控制变量集合， 在此处具体为需要反演的A和C； ∇为散度算子， dB的含义同上。因为控制变量的变化范围可能很大， 为了收缩其参数空间加速运算， 一种常见的操作是取对数， 也就是对log10(A)和log10(C)求解而不是直接针对A和C求解。

乌阿冰流模式最早成功应用于西南极海洋性冰盖的触地线稳定性研究， 近年来应用领域已逐步拓展， 比如冰架流动格局的时空演变［27］， 冰架减薄对触地线上游冰川的动力响应［28］， 裂隙的应力分布与传播机制［29］， 格陵兰峡湾冰流的表面融化［30］等。另外一方面该模式也陆续参加了一些系列冰川模式比较计划［16］， 均取得良好的效果。

2.2　数据介绍

乌阿冰流模式的输入数据主要包括计算区域边界、 冰厚、 基岩地形、 表面高程、 积累率、 融化率、 初始冰流速、 冰厚变化率等， 详见表1。为了能够代表冰架当前的动力状态， 本文尽量搜集并使用了最新的公开数据。由于内陆地区的补给冰流速度缓慢， 对冰架动力的影响较小， 本研究的计算域并没有拓展到最上游的分冰岭， 但在冰架前端根据最近的遥感影像做了适当的延伸。值得指出的是冰厚和基岩地形使用了2013年发布的BedMap2数据集， 而表面高程使用了2019年发布的REMA数据集， 在两者间不匹配的区域， 我们保持冰厚不变而调整了底部地形用以适应表面高程。同时， 第二版MEaSUREs全南极冰流速度数据集用于提供反演试验的初始速度， 并用来验证最终模拟的流速。

表1   本文所用到的数据列表

Table 1  The datasets used in this paper

变量数据集时间分辨率*
前端边界Landsat-8/Sentinel-12018/2019年15/40 m
流域范围Zwally’s Basin［31］2012年矢量线段
冰厚h、 基岩BBedMap2［32］多年合成， 2013年发布1 000 m
表面高程sREMA［33］多年合成， 2019年发布1 000 m
表面积累率asRACMO2.3p2［34］1976—2016年， 气候平均态27 km
底部融化率ab温家洪等估算结果［35］多年合成， 2010年发布20 km
初始流速（U， V）MEaSUREs-V2［36］多年合成， 2017年发布450 m
冰厚变化（h˙）Matt King等估算结果［37］多年合成， 2009年发布高度计脚印点
触地线GLMOA产品［38］2014年矢量线段

注：* 计算范围内所有变量均重新双线性插值到自适应不规则网格

新窗口打开| 下载CSV


由于东南极纬度高终年极寒， 表面融化和升华造成的物质损耗极小， 表面净物质平衡约等于表面积累。考虑到现有的表面积累率和底部融化率结果较为粗糙， 本文做了双线性插值。依据King等［37］基于长时间序列雷达高度计的研究结果， 最近40年来埃默里冰架总体处于稳定状态， 所以我们假定冰厚不随时间变化（h˙=0）。有关模式输入的几何变量， 详见图2示意， 以上输入数据最终都插值到图1（b）中的不规则三角网的节点上面。另外在全部模拟试验中始终保持海平面为0 m， 亦即忽略了潮汐波动对冰架动力的影响。

图2

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

图2   冰川几何的模式输入变量示意图

s为冰川表面高程， b为冰川底面高程， S为海洋表面高程， B为基岩地形高程， 冰川厚度为h = s - b， 海洋深度为H = S - B，吃水深度为d = S - b， 出水高度为f = s - S

Fig.2   The diagram showing the variables as model inputs for ice flow

Including the ice surface （s）， ice bottom （b）， sea level （S）， beckrock （B）， ice thickness （h = s - b）， ocean bathymetry （H = S - B）， ice draft （d = S - b）， ice freeboard （f = s - S）


2.3　试验设计

首先我们选择初始流速计算的应变率场、 漂浮与否、 到触地线距离三个参数自适应优化网格。具体规则如下： 内陆最粗的网格为25 km， 触地线附近逐步加密到3 km， 上游冰流随应变率自适应到5 km， 下游冰架在满足应变率规则的情况下设置小于8 km。图1（b）中边缘的红色箭头代表冰流微分方程组的狄利克雷边界条件， 即设置所有239个接地的边界体元速度为0 m·a-1； 根据海平面高程始终为0 m， 包括反演试验和预测实验， 自由漂浮的冰架前端应用诺伊曼边界条件， 即

σ⋅n̂=- ρwgS-b⋅n̂

式中： σ为崩解前端的柯西应力张量； n̂为指向海洋的边界单位法向矢量； ρw为海水密度； S-b为海平面与冰架底层之差， 即冰架的吃水深度（图2）。

本文的试验设计如图3所示， 分为反演试验与预测试验两个部分。对于反演试验， 使用的MEaSUREs-V2流速数据集中自带了误差分析（方程5中的Umeas与Vmeas）。试验中不确定可能来源于假设的冰架厚度稳定（h˙=0）。因此我们同时设置了四种冰厚变化率用于敏感性测试， 分别为不加垂向约束（即h˙error→∞）和h˙error取1、 0.1、 0.01 m·a-1。在反演试验中， 以模拟流速与观测流速差值最小为标准获得黏度因子A和底部滑动系数C， 它们也将被用于预测试验中。预测试验在表面积累率与底部融化率不变的条件下， 分为两种子试验。子试验一继续敏感性试验， 测试不同误差约束下对冰架自由漂移的影响； 子试验二改变原始计算域的几何边界， 模拟埃默里解体情境下对上游冰流的动力响应。两个预测子试验前向积分时间都是100年。

图3

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

图3   本文试验设计流程图

先进行4次反演冰川黏度和底部滑动系数， 并以遥感观测的MEaSUREs-V2作为基准， 评价模拟的冰流速度， 最后比较4套反演系数的自由漂移量， 以及冰架解体100后的海平面上升量

Fig.3   The flowchart of experiments in this paper

Firstly， the viscosity rate factor （A） and the basal slipperiness （C） are retrieved from 4 sets of inverse runs. Then we evaluate the modeled ice velocity against the remote sensing measurements of MEaSUREs-V2 datasets. In the end， the optimal inverse parameters is selected based on the least drifting bias before calculating the sea level contribution in the ice shelf disintegration scenario


3 结果与讨论

本节我们对反演和预测这两部分试验的计算结果进行比较分析。

3.1　网格分析

网格系统是所有有限元系统求解微分方程的基石。对于固定的计算域空间， 一套疏密得当的网格能够在抓住冰流动力的主要特征的同时减少运算时间。本研究使用的网格共有26 538个元素， 13 423个节点。从图1（b）来看， 本研究的加密规则较为合理地反映出埃默里冰架上下游的冰流动力状况， 例如在靠近上游边界的内陆冰盖， 网格尺度在20 km以上， 三条补给冰川的空间特征都能反映在网格分布上， 尤其是通量最大的兰伯特冰川， 一直向上游延伸到边界。网格最密集处为三条冰川汇流的触地线， 此处山脉密集， 特殊的地形使得此处的应变率加大， 为了适应其距离的变化梯度， 网格尺寸在5 km 左右。随着流向下游冰架逐步展宽， 网格经历了由粗糙到精细的变化， 从图1（a）中遥感影像显示的表面裂隙可以推断在冰架前缘的应变率再次增大， 冰流模式也在此处设置精细网格开展模拟。对触地线而言， 图1（a）中黑线为观测触地线， 图1（b）中蓝线为模拟触地线， 两者基本吻合， 差异主要集中在冰架中部， 例如冰架东侧MOA触地线中并未包括吉洛克岛但模拟的触地线则向下游延伸。

3.2　反演评价

反演是一个不断迭代极小化代价函数J的过程， 图4结果表明经过100次迭代， 乌阿冰流模式中目标函数下降的总体趋势。由于初始的冰川黏度A与底部滑动系数C未必符合实际情况， 目标代价函数的值高达3 616（无量纲）， 在最初的10次迭代过程中函数值呈现对数级缩减到24.47， 并在15次左右到达“拐点”， 随后的函数逐步达到收敛。从第20到第100个循环， 代价函数仅从14.92下降到9.27， 对应的控制变量A和C也接近收敛状态这时的下降特征为线性， 图4中的子图放大展示了第70 ~ 80个循环的变化。

图4

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

图4   目标代价函数随反演迭代的下降趋势

主图坐标轴为对数尺度， 放大子图则为线性尺度

Fig.4   The decreasing pattern of the objective cost function in iterations of the inverse run

The main plot employs logarithmic scale and the inset employs the linear scale


评价冰流模式的反演结果的主要手段是利用观测流速进行验证［11］。图5（a）和5（b）分别展示了利用遥感影像合成的MEaSUREs-V2与乌阿模式的模拟流速， 从空间上来看模拟流速较为准确地反映了实际冰流的分布特征， 尤其抓住上游三支冰川分支的汇流区域流速较快的特点， 但与观测值相比， 仍然出现了不同程度的低估。尤其是刚刚进入触地线区域， 低估达到了-70 m·a-1； 另一方面， 模式模拟流速在埃默里冰架东西两侧的缓慢流速区出现高估。此外， 观测流速数据中［图5（a）］中出现很多由于崎岖的基岩导致的孤立“黏滞点”（sticky spot）， 尤其是在莫森陡崖和查尔斯王子山脉的上游区域。模拟流速数据中高频特征被不同程度地抑制。图5（c）展示了模拟与观测流速之间的差异， 如上所述误差主要集中在北侧和西侧的触地线上游（图中两个红圈处）， 其中兰伯特冰川与米勒冰川的误差更显著， 而费舍尔冰川由于其本身流速缓慢而更贴近观测。总体来说， 基于乌阿冰流模型的流速场反演能够刻画冰流分布特征， 但趋势更为平滑。

图5

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

图5   观测冰流速与模拟冰流速对比图

Fig.5   The ice flow comparison between the MEaSUREs-V2 surface observation compiled with satellite remote sensing （a）；the depth-averaged velocity simulated by Úa model （b）； the difference between the modeled and measured flow speed （c）


在每个网格节点分别获取观测和模拟流速的水平（U）和竖直（V）分量并作差， 以此来评估模拟流速场的质量：理论上来讲， 如果反演结果足够精确两者之差应当接近0， 也就是失配函数［式（5）］的前两项接近0。从图6（a）和6（b）来看， 绝大多数的模拟结果都与观测流速吻合， 误差直方图特征接近以0为均值的正态分布（标准差分别为13.14和12.44 m·a-1）， 这表明模拟的流速误差主要以随机误差为主， 系统性的偏差很小（均值分别为-0.28和-1.14 m·a-1）。两者速度对比的散点图也能体现上述结论， R2高达0.99， 均方根误差为13.55 m·a-1， 在13 423个格点中仅有81个离群点的流速离差超过了50 m·a-1， 占比0.6%。

图6

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

图6   冰流速度对比图（红色虚线为1∶1比例）

（b） the scatter plot of modeled flow speed against the observation， and the red dotted line represent 1∶1 line； （c） same as （a） but for V component］

Fig.6   The ice velocity comparison［（a） the difference of U component between modeled and observed velocities；


图7（a）、 7（b）分别是反演过程迭代100次所求解的冰流黏度因子、 基岩滑动因子， 两者最突出的表现是靠近触地线处值突然增大， 这与观测事实基本一致。因为该区域是从接地冰盖向漂浮冰架的过渡地带， 相当于突然失去底部摩擦力而导致冰流加速， 由式（2）和式（3）可知为了适应加速及其引发的应变率突增， 对应的A和C也要相应增大。

图7

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

图7   反演得到的空间场黏度因子log10（A）（a）与底部滑动系数log10（C）（b）

Fig.7   The inversed spatial field of the viscosity rate factor log10（A）（a）； and the basal slipperiness factor log10（C）（b）


值得指出的是， 由于BedMap2数据集基于直接钻孔、 探冰雷达、 重力反演等多种资料合成［32］， 插值过程中不可避免地产生较大误差。在精度较差的区域， 作为底部边界条件的冰下地形同样造成上游冰川物理场偏离目前的表面遥感观测数据。在前人模拟工作中， 有时会根据实际需要手动修正地形， 比如Pittard等直接用RTOPO数据替代了触地线附近的BedMap2数据［39］。此外， 底部地形也随着观测数据的不断积累而得到更加严格约束， 最新发布的BedMachine数据集在本研究区内精度是否有所提高还有待未来检验［40］。另一方面， BedMap2的迭代产品BedMap3也即将完成， 为将来南极冰盖模拟工作提供了更多数据选择。

3.3　百年预测

在反演得到的A和C空间场的基础上， 本研究开展预测试验。子试验一采用4组冰厚变化率用来测试参数h˙的不确定性对自由漂移模拟的影响， 子试验二用于模拟冰架解体后的动力响应。在预测实验中， 选取冰架露出海面的体积（volume-above-floatation）作为刻画上游冰川与埃默里冰架系统物质总量的指标， 因为这一部分冰的体积（乘以平均密度即为物质量）能够直接影响到海平面的上升或下降。

如果将严格符合遥感观测的冰流模型看做一个完整自洽的系统， 在初始化模型时由于输入数据带入的误差而造成系统存在残差无法完全闭合， 并随着前向积分累积增大； 如上节所述， 尤其是当底部地形不准确时， 初始化完成后的系统刚进入预测状态时需要大幅调节冰厚。以上现象是由数值性的误差引起， 而非现实中存在的物理过程， 文献中通常把这种现象称为非物理性畸变（non-physical）［41］。子试验一对比4组自由漂移试验。有别于以往仅考虑流速的二维冰流模式， 新一代乌阿冰流模式增加了反演失配函数中关于冰厚变化的控制［式（5）， 相当于2.5维］， 这个特点使得数值系统能够在一定程度上减少由于输入的底部地形或冰厚而造成的非物理性畸变的影响。本研究假定目前冰架系统的厚度趋于稳定（即方程5中的h˙meas=0）但我们并对此设定的不确定性未知（h˙error）， 于是以传统不加约束为参考， 分别将其设置为1/0.1/0.01 m·a-1测试参数敏感性。理想情况下， 冰流系统在漂移试验中应当保持不变， 物质量为一条平直线； 但由于反演过程无法达到完美收敛， 一方面方程系统从反演迭代回到前向积分过程中， 地形或冰厚的随机误差会产生非物理性的畸变； 另外一方面反演造成的方向性偏差又往往随着时间积累增大。图8（a）上半部分显示了4组物质量的演化路径。其中， 前20年约为非物理畸变过程， 前5年减少约1 950 ~ 2 341 Gt（约合5.39 ~ 6.47 mm等效海平面变化）， 可以作为冰流反演的随机误差； 100年后4条曲线偏差积累不同， 对1、 0.1、 0.01 m·a-1三组对应的海平面变化分别约为+3.25 mm、 -0.08 mm和-0.98 mm。模拟结果表明冰厚变化约束越强（h˙error→0）冰流系统能更快地从非物理性畸变中恢复自然漂移， 但偏差积累也更大， 甚至超过原本随机误差。此外， 当h˙error→∞与h˙error=1 m·a-1演化曲线几乎重合， 这与当代卫星测高的表面高程误差量级一致（底部基岩不变， 冰厚变化误差直接反应在表面高程上面）。由于在100年尺度上h˙error对应的总误差最小， 本研究以这一组参数作为自由漂移的结果绘制应变率场［图8（b）］， 结果显示应变率较大的地方主要集中在冰流交汇处， 触地线上游和崩解前端； 尤其对崩解前端只存在外向的拉伸张量（蓝色）， 这与遥感观测到的裂隙发育亦可互相印证。

图8

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

图8   模式预测结果［分图（b）中， 蓝色箭头代表拉伸， 红色箭头代表压缩］

Fig.8   The model prediction result ［（a） the ice mass evolution trajectory in different scenarios； （b） the strain rate field after 100 years’ drifting， with the blue arrows for extension and red ones for compression； （c） the ice flow acceleration after disintegration］


子试验二模拟了埃默里冰架解体的理想极端情景， 图8（c）展示了相对于初始流速冰架解体100年后的冰流加速情况。上游补给冰流均出现不同程度的加速， 其中触地线南端达到600 ~ 700 m·a-1， 流速几乎翻了一倍， 兰伯特冰川响应的加速信号向上游延伸约150 km， 米勒冰川与其表现类似。与反演速度误差的表现较为类似， 费舍尔冰川的加速趋势更弱， 考虑到费舍尔冰川上游出露的裸岩面积更大， 三条补给冰川之间这种固有差异很可能受到底部地形内在控制。相对于漂移试验， 触地线后退约120 km。从总物质量的角度， 由冰架解体造成的上游物质损耗16 423 Gt， 将会引起海平面上升（45.36 ± 0.08） mm。

4 结论与展望

数值模式是刻画山地冰川和极地冰盖动力过程和预测未来冰流演化的重要工具。本文向读者介绍了新一代冰流模式乌阿的基本原理， 并基于最新的数据集， 针对东南极兰伯特冰川-埃默里冰架系统开展反演和预测两部分试验。研究表明乌阿模式中的网格设置规则能够刻画冰流动力特征的空间分布：在流速梯度大的区域自适应精细化， 而在内陆流速慢的冰盖区域能够放大以节约计算资源。模式进行100次迭代后， 反演得到冰流速场与遥感卫星的观测结果总体吻合， 速度差异普遍低于50 m·a-1， 相关系数达到0.99， 但失去了部分高频的黏滞特征。在设置4组不同冰厚变化率不确定性的基础上， 分别进行100年前向积分漂移模拟， 发现较强的冰厚约束能够使系统更快地从反演误差造成的非物理畸变中恢复， 但偏差的积累也使得结果偏离0值。模拟结果显示， 在缺少观测的情况下需要根据积分时间选择适当的误差量级。最后本文还研究模拟冰架解体后的冰流状况， 当冰架解体100年后上游冰流普遍加速100 ~ 700 m·a-1， 排泄通量增大冰盖物质损耗能够造成海平面上升（45.36 ± 0.08） mm， 同时伴随着触地线相对退缩。

本研究的试验设计还比较简单， 比如在预测试验中我们保持表面积累率和底部融化率不变， 没有考虑气候变化情境下可能导致降雪和海水温度的同时增加。冰盖模式与气候、 海洋模式的耦合是增进对地球系统理解和海平面预测精度的重点发展方向。未来我们计划基于中等复杂度的羽流参数化方案， 采用最新的CMIP6气候模式结果作为外界强迫因子， 更准确地刻画埃默里冰架流域的演化状态。

参考文献

The IMBIE team.

Mass balance of the Antarctic Ice Sheet from 1992 to 2017

［J］. Nature， 2018， 558（7709）： 219 - 222.

[本文引用: 1]

Li T， Liu Y， Li T， et al.

Antarctic surface ice velocity retrieval from MODIS-Based mosaic of Antarctica （MOA）

［J］. Remote Sensing， 2018， 10（7）： 1045 - 1063.

[本文引用: 1]

Liu Y， Moore J C， Cheng X， et al.

Ocean-driven thinning enhances iceberg calving and retreat of Antarctic ice shelves

［J］. Proceedings of the National Academy of Sciences， 2015， 112（11）： 3263 - 3268.

[本文引用: 1]

Li Teng， Cheng Xiao， Liu Yan， et al.

An overview of the Ice Sheet Model Intercomparison Project （ISMIP） in CMIP6

［J］. Climate Change Research， 2020， 16（2）： 255 - 262.

[本文引用: 1]

李腾， 程晓， 刘岩， 等.

CMIP6冰盖模式比较计划（ISMIP）概况与评述

［J］. 气候变化研究进展， 2020， 16（2）： 255 - 262.

[本文引用: 1]

Bindschadler R A， Nowicki S， Abe-Ouchi A， et al.

Ice-sheet model sensitivities to environmental forcing and their use in projecting future sea level （the SeaRISE project）

［J］. Journal of Glaciology， 2013， 59（214）： 195 - 224.

[本文引用: 1]

Pattyn F， Perichon L， Durand G， et al.

Grounding-line migration in plan-view marine ice-sheet models： results of the ice2sea MISMIP3d intercomparison

［J］. Journal of Glaciology， 2013， 59（215）： 410 - 422.

[本文引用: 1]

Li Huilin， Li Zhongqin， Shen Yongping， et al.

Glacier dynamic models and their applicability for the glaciers in China

［J］. Journal of Glaciology and Geocryology， 2007， 29（2）： 201 - 208.

[本文引用: 1]

李慧林， 李忠勤， 沈永平， 等.

冰川动力学模式及其对中国冰川变化预测的适应性

［J］. 冰川冻土， 2007， 29（2）： 201 - 208.

[本文引用: 1]

Zhou Zaiming， Li Zhognqin， Li Huilin， et al.

The flow velocity features and dynamic simulation of the Glacier No.1 at the Headwaters of Urumqi River， Tianshan Mountains

［J］. Journal of Glaciology and Geocryology， 2009， 31（1）： 55 - 61.

[本文引用: 1]

周在明， 李忠勤， 李慧林， 等.

天山乌鲁木齐河源区1号冰川运动速度特征及其动力学模拟

［J］. 冰川冻土， 2009， 31（1）： 55 - 61.

[本文引用: 1]

Tang Xueyuan， Zhang Zhanhai， Sun Bo.

Progress and prospect in numerical modelling of the Antarctic Ice-Sheet

［J］. Chinese Journal of Polar Research， 2006， 18（4）： 290 - 300.

[本文引用: 1]

唐学远， 张占海， 孙波.

南极冰盖数值模拟研究进展与展望

［J］. 极地研究， 2006， 18（4）： 290 - 300.

[本文引用: 1]

Tang Xueyuan， Sun Bo， Zhang Zhanhai， et al.

GLIMMER Antarctic Ice Sheet Model， an experimental research of moving boundary condition

［J］. Journal of Glaciology and Geocryology， 2007， 29（6）： 905 - 913.

[本文引用: 1]

唐学远， 孙波， 张占海， 等.

南极冰盖 GLIMMER 模式移动边界试验研究

［J］. 冰川冻土， 2007， 29（6）： 905 - 913.

[本文引用: 1]

Ji Qingyuan， Wang Bangbing， Sun Bo.

Applicability PISM for velocity analysis of the Amery Ice Shelf， East Antarctica

［J］. Chinese Journal of Polar Research， 2015， 27（3）： 229 - 236.

[本文引用: 2]

季青原， 王帮兵， 孙波.

PISM 冰盖模式对Amery冰架流速场模拟的适用性

［J］. 极地研究， 2015， 27（3）： 229 - 236.

[本文引用: 2]

Zhang Liangfu， Tang Xueyuan， Yang Shuhu.

Numerical simulations of East Antarctic Ice Sheet based on the Elmer/Ice Model

［J］. Chinese Journal of Polar Research， 2017， 29（3）： 390 - 398.

[本文引用: 1]

张良甫， 唐学远， 杨树瑚， 等.

基于Elmer/Ice冰盖模型的南极Gamburtsev山脉Lambert冰流区域的数值模拟研究

［J］. 极地研究， 2017， 29（3）： 390 - 398.

[本文引用: 1]

Williams M J， Grosfeld K， Warner R C， et al.

Ocean circulation and ice-ocean interaction beneath the Amery Ice Shelf， Antarctica

［J］. Journal of Geophysical Research： Oceans， 2001， 106（C10）： 22383 - 22399.

[本文引用: 1]

Gong Y， Cornford S L， Payne A J.

Modelling the response of the Lambert Glacier-Amery Ice Shelf system， East Antarctica， to uncertain climate forcing over the 21st and 22nd centuries

［J］. The Cryosphere， 2014， 8（3）： 1057 - 68.

[本文引用: 1]

Gudmundsson G H， Krug J， Durand G， et al.

The stability of grounding lines on retrograde slopes

［J］. The Cryosphere， 2012， 6（6）： 1497 - 1505.

[本文引用: 2]

Levermann A， Winkelmann R， Albrecht T， et al.

Projecting Antarctic’s contribution to future sea level rise from basal ice shelf melt using linear response functions of 16 ice sheet models （LARMIP-2）

［J］. Earth System Dynamics， 2020， 11（1）： 35 - 76.

[本文引用: 2]

Wang Chenghai， Cheng Rong， Zhao Wen， et al.

Research progress on the glacial dynamics models

［J］. Journal of Glaciology and Geocryology， 2019， 41（4）： 1 - 9.

[本文引用: 1]

王澄海， 程蓉， 赵文， 等.

冰川动力学模式模型进展及研究

［J］. 冰川冻土， 2019， 41（4）： 1 - 9.

[本文引用: 1]

Yu J， Liu H， Jezek K C， et al.

Analysis of velocity field， mass balance， and basal melt of the Lambert Glacier-Amery Ice Shelf system by incorporating Radarsat SAR interferometry and ICESat laser altimetry measurements

［J］. Journal of Geophysical Research： Solid Earth， 2010， 115（B11）： 1 - 16.

[本文引用: 1]

Shi Jiuxin.

A review of ice shelf-ocean interaction in Antarctica

［J］. Chinese Journal of Polar Research， 2018， 30（3）： 287 - 302.

[本文引用: 1]

史久新.

南极冰架-海洋相互作用研究综述

［J］. 极地研究， 2018， 30（3）： 287 - 302.

[本文引用: 1]

Tang Chengjia， Li Yuansheng， Chen Zhenlou， et al.

A review on studies of Antarctic Shelves and advances in Chinese research on Amery Ice Shelf

［J］. Chinese Journal of Polar Research， 2008， 20（3）： 265 - 274.

[本文引用: 1]

唐承佳， 李院生， 陈振楼， 等.

南极冰架研究现状与埃默里冰架研究展望

［J］. 极地研究， 2008， 20（3）： 265 - 274.

[本文引用: 1]

Zhao C， Cheng X， Hui F， et al.

Monitoring the Amery Ice Shelf front during 2004-2012 using ENVISAT ASAR data

［J］. Adv Polar Sci， 2013， 24（2）： 133 - 137.

[本文引用: 1]

Teng L， Yan L， Xiao C.

Recent and imminent calving events do little to impair Amery ice shelf’s stability

［J］. Acta Oceanologica Sinica， 2020， 42（5）： 168 - 170.

[本文引用: 1]

Zhou C， Liang Q， Chen Y， et al.

Mass balance assessment of the Amery Ice Shelf Basin， East Antarctica

［J］. Earth and Space Science， 2019， 6（10）： 1987 - 1999.

[本文引用: 1]

Gudmundsson G.

Ice-shelf buttressing and the stability of marine ice sheets

［J］. The Cryosphere， 2013， 7（2）： 647 - 655.

[本文引用: 1]

Bueler E， Brown J.

Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model

［J］. Journal of Geophysical Research： Earth Surface， 2009， 114（F3）： 1 - 21.

[本文引用: 1]

Geuzaine C， Remacle J.

Gmsh： A 3-D finite element mesh generator with built-in pre-and postprocessing facilities

［J］. International Journal for Numerical Methods in Engineering， 2009， 79（11）： 1309 - 1331.

[本文引用: 1]

Gudmundsson G H， De Rydt J， Nagler T.

Five decades of strong temporal variability in the flow of Brunt Ice Shelf， Antarctica

［J］. Journal of Glaciology， 2017， 63（237）： 164 - 175.

[本文引用: 1]

Reese R， Gudmundsson G H， Levermann A， et al.

The far reach of ice-shelf thinning in Antarctica

［J］. Nature Climate Change， 2018， 8（1）： 53 - 59.

[本文引用: 1]

De Rydt J， Gudmundsson H， Nagler T， et al.

Recent rift formation and impact on the structural integrity of the Brunt Ice Shelf， East Antarctica

［J］. The Cryosphere， 2018， 12（2）： 505 - 520.

[本文引用: 1]

Rathmann N， Hvidberg C， Solgaard A， et al.

Highly temporally resolved response to seasonal surface melt of the Zachariae and 79N outlet glaciers in northeast Greenland

［J］. Geophysical Research Letters， 2017， 44（19）： 9805 - 9814.

[本文引用: 1]

Zwally H， Jay M B G， Beckley M A， et al.

Antarctic and Greenland drainage systems

［EB/OL］. （2020-02-20） ［2020-03-12］. .

[本文引用: 1]

Fretwell P， Pritchard H D， Vaughan D G， et al.

Bedmap2： improved ice bed， surface and thickness datasets for Antarctica

［J］. The Cryosphere， 2013， 7（1）： 375-393.

[本文引用: 2]

Howat I M， Porter C， Smith B E， et al.

The reference elevation model of Antarctica

［J］. Cryosphere， 2019， 13（2）： 665 - 674.

[本文引用: 1]

Van Wessem J M， Van de Berg W J， Noël B P Y， et al.

Modelling the climate and surface mass balance of polar ice sheets using RACMO2， part 2： Antarctica （1979-2016）

［J］. The Cryosphere， 2017， 12（4）： 1479 - 1498.

[本文引用: 1]

Wen J， Wang Y， Wang W， et al.

Basal melting and freezing under the Amery Ice Shelf， East Antarctica

［J］. Journal of Glaciology， 2010， 56（195）： 81 - 90.

[本文引用: 1]

Mouginot J， Rignot E， Scheuchl B， et al.

Comprehensive annual ice sheet velocity mapping using Landsat-8， Sentinel-1， and RADARSAT-2 Data

［J］. Remote Sensing， 2017， 9（4）： 364 - 373.

[本文引用: 1]

King M A， Coleman R， Freemantle A J， et al.

A 4-decade record of elevation change of the Amery Ice Shelf， East Antarctica

［J］. Journal of Geophysical Research， 2009， 114（F1）： 1 - 13.

[本文引用: 2]

Scambos T， Haran T， Fahnestock M， et al.

MODIS-based Mosaic of Antarctica （MOA） data sets： Continent-wide surface morphology and snow grain size

［J］. Remote Sensing of Environment， 2007， 111（2-3）： 242 - 257.

[本文引用: 1]

Pittard M L， Galton-Fenzi B K， Roberts J L， et al.

Organization of ice flow by localized regions of elevated geothermal heat flux

［J］. Geophysical Research Letters， 2016， 43（7）： 3342 - 3350.

[本文引用: 1]

Morlighem M， Rignot E， Binder T， et al.

Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet

［J］. Nature Geoscience， 2020， 13（2）： 132 - 137.

[本文引用: 1]

Goelzer H， Nowicki S， Edwards T， et al.

Design and results of the ice sheet model initialisation initMIP-Greenland： An ISMIP6 intercomparison

［J］. The Cryosphere， 2018， 12（4）： 1433 - 1460.

[本文引用: 1]

/

〈 〉

陇ICP备2021001824号-10
　 甘公网安备 62010202000676号
网站版权 © 《冰川冻土》编辑部
地址：兰州市天水中路8号　
邮编：730000　电话：0931-8260767　
E-mail：edjgg@lzb.ac.cn
总访问量 　 今日访问 　 在线人数
技术支持：北京玛格泰克科技发展有限公司