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

doi:10.7522/j.issn.1000-0240.2019.0042

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

王澄海^,¹^,², 程蓉¹^,², 赵文¹^,², 孙超¹^,²

1.兰州大学大气科学学院甘肃省干旱气候变化与减灾重点实验室, 甘肃兰州 730000

2.兰州大学地球系统模式研发中心, 甘肃兰州 730000

Research progress on the glacial dynamics models

WANG Chenghai^,¹^,², CHENG Rong¹^,², ZHAO Wen¹^,², SUN Chao¹^,²

1.Research and Development Center of Earth System Model （RDCM），Lanzhou University，Lanzhou 730000，China

2.Key Laboratory of Arid Climate Change and Disaster Reduction of Gansu Province，College of Atmospheric Sciences，Lanzhou University，Lanzhou 730000，China

编委: 周成林

收稿日期: 2018-11-20 修回日期: 2019-04-20

基金资助:

国家自然科学基金. 91837205. 41661144017. 41805032

Received: 2018-11-20 Revised: 2019-04-20

作者简介 About authors

王澄海（1961-），男，甘肃天水人，教授，2002年在中国科学院寒区旱区环境与工程研究所获博士学位，从事陆气相互作用、寒区气候学、数值模式及模拟研究.E-mail：wch@lzu.edu.cn. , E-mail：wch@lzu.edu.cn

摘要

冰冻圈是气候系统中的一个重要圈层，其中冰川又是冰冻圈的重要组成部分，冰川、尤其是山地冰川的本构方程和建模一直是冰川动力学的核心任务。首先，简要回顾冰川模型的研究和发展，简要介绍了基于Navier-Stokes方程耦合温度场的三维冰川模型。然后，介绍了冰川建模过程中的常用的静水压力近似、一阶近似、浅冰近似等的基本概念，总结了冰川的动力数值模式建立的主要方法，对于常用的GLIMMER冰盖模式的物理框架及其应用进行了介绍。最后，针对目前的简化模型难以准确地描述山地冰川的物理过程及其变化的问题，提出了一个基于全Navier-Stokes方程的山地冰川模型及其动力框架、边界条件处理的设想。本文可为建立、发展冰川及冰架模型，尤其建立和发展山地冰川模型提供基础知识和参考。

关键词： 山地冰川 ; 冰川模型 ; 发展 ; 展望

Abstract

The cryosphere is an important component of the climate system， and the glaciers are the crucial part of the cryosphere. The constitutive equations and modeling of glaciers， especially mountain glaciers， have always been the core task of glacial dynamics. This study briefly reviewed the research and development on glacier model which includes the three-dimensional glacier model used the Navier-Stokes equation coupled with the temperature field and some basic concepts of model simplification， such as the hydrostatic pressure approximation， the first order approximation， and the shallow ice approximation. Based on the summary on the primary method of establishing the dynamic numerical model of glaciers this study also introduced a typical example of ice sheet model GLIMMER and its application in research area. Finally， considering the poor representation of simplified numerical model of mountain glaciers， especially in the descriptions of their physical processes and spatiotemporal variations， a dynamical framework of mountain glacial model based on the complete Navier-Stokes equation and the corresponding treatment of boundary condition are proposed in this study. This study could provide basic knowledge and reference for the establishment and development of the glacier and ice shelf models， especially for the mountain glacier models.

Keywords： mountain glaciers ; glacier model ; research progress ; future development

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

本文引用格式

王澄海, 程蓉, 赵文, 孙超. 冰川动力学模式模型进展及研究[J]. 冰川冻土, 2020, 42(1): 43-52 doi:10.7522/j.issn.1000-0240.2019.0042

WANG Chenghai, CHENG Rong, ZHAO Wen, SUN Chao. Research progress on the glacial dynamics models[J]. Journal of Glaciology and Geocryology, 2020, 42(1): 43-52 doi:10.7522/j.issn.1000-0240.2019.0042

0 引言

冰冻圈是气候系统中重要的圈层之一。冰冻圈中以冰川和大气之间的作用最为直接，受气候变化影响也最为显著。近年来，观测到的两极冰川（架、盖）消融引起的海平面上升，洋流发生变化引起了人们的广泛关注。在气候系统中冰冻圈的变化估算和模拟是最具不确定性的，也是气候系统中亟需改进模拟的圈层。图1可以看出，冰冻圈过程的时空尺度分布差异巨大，从几天到千年，几米到数千公里，其物理本质和变化具有及其复杂的过程。

图1

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

图1 冰冻圈的时空分布图（根据IPCC， 2007修改）

Fig.1 Spatial and temporal distributions of the cryosphere （changed according to IPCC， 2007）

冰冻圈和气候之间相互作用，相互影响。当气候发生变化时，冰川会发生一系列响应；同时，冰川变化又会影响全球的总反照率、地形的结构和其上的大气层结，进而对区域大气环流状况产生影响。同样，冰冻圈和气候之间的相互作用的时间尺度可达几十年甚至几千年。冰冻圈和气候之间的这种关系我们简单的归纳为图2。由于冰冻圈各分量，尤其是冰川，位于高纬度和高海拔地区，观测资料缺少，研究和认识冰冻圈变化及其与气候之间相互作用的有效途径是数值模式。

图2

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

图2 气候与冰冻圈相互关系示意图

Fig.2 Relationship between climate and cryosphere

1 冰川模式的发展

冰川观测研究起源于18世纪，早期的研究以观测为主，主要是收集关于冰川演变发展的数据。Saussure（1740—1899）就曾经对冰川的滑动进行了观测。由于冰川的演变非常缓慢，现有的观测只能提供过去一段时间内的有限了解，要定量确定冰川未来的变化，只有通过建立冰川的数学物理模型，进行长期数值积分的模拟研究。

因此， 20世纪以来，人们进行了大量探索，建立、发展了冰川数值模式用来模拟冰川和冰盖在各种条件下的物理特征及其运动和变化规律。Kuhn（1762—1825）试图从理论的角度分析冰川的运动，他应该是第一个将热力学理论和力学理论结合起来分析冰川运动的学者。此后， Scheuchzer（1672—1733）， Charpentier（1786—1855）和Agassiz（1807—1873）提出了膨胀理论，从理论上解释冰川的运动。这一理论认为冰川运动的驱动力是缝隙中水的冻结膨胀效应，以及冰川中的毛细作用。早期的工作对冰川的变化并没有细化，另外，冰盖（包括冰架）是独立的固体单元，便于建模试验，因此主要关注的是冰盖的变化。此类模型经历了以下几个发展和研究阶段：

第一个阶段，开始于20世纪60—70年代，主要是对浅冰近似的模型的应用，包括Oerlemans在70—80年代用一维的冰川流线模型模拟冰川末端的变化；

第二个阶段，二维的冰川模型，主要用来得到观察场的诊断结果（Waddington等^［1］、 Johannesson等^［2］、 Abe-Ouchi^［3］）；

第三个阶段，增加了垂直积分的平面模型，用来耦合气候模型（Fastook等^［4］、 Verbitsky等^［5］、 MacAyeal^［6］、 Marsiat^［7］）；

第四阶段，三维模型的建立，其包括了冰川流场和温度场的耦合（Huybrechts^［8］、 Calov^［9］、 Fabre等^［10］）。

尽管这些模型使用同样的守恒定律，但使用了不同的关系式和简化方法，最主要的区别是选择全Navier-Stokes方程进行计算^［11］，但还是忽略冰川的 $z$ 平面的剪切应力 $τ_{x z}$ 和 $τ_{y z}$ （远小于 $τ_{z z}$ ）进行计算以及是否需要显式计算流线进行分类等^［12］。1997年提出的EISMINT（European Ice Sheet Modeling Initiative）科学计划，其主旨就是比较、评估当时的冰盖模型，主要比较的参数包括冰层厚度、速度和温度。

伴随着计算机计算量的提高和有限元技术的发展，为了掌握山地冰川在不同气候情景下的进退消融，国际上提出了许多用有限元方案求解全Navier-Stokes方程的模式。目前，使用有限元方法求解冰川问题主要是基于成熟的有限元商业软件包，通用的商业软件包含三个步骤：前处理、方程求解以及后处理。前处理可以方便的建立所求解问题的几何模型，例如可建立复杂的冰川几何模型，或者可以通过与专业的建模软件结合（如Ug， Pro/E等），建立复杂的冰川模型。方程的求解是问题的核心，冰川问题需要求解全Navier-Stokes方程并考虑Glen定律，这就限制了一些商业软件包的使用，因为有些有限元软件不能求解此类问题；后处理则是能将结果很好得呈现出来，这部分在现有的商业软件包中都做得比较好。多物理场耦合有限元软件COMSOL在冰川问题的求解中使用得较多，但是COMSOL的缺点是不能将气候变化对冰川的影响考虑进去。此外， Jarosch^［13］提出的Icetool模型，以及由Zwinger等^［14］在有限元软件ELMER（Open Source Finite Element Software for Multiphysical Problems）的基础上建立的有限元冰川模型。该类模型考虑了冰川中存在的所有应力，使得该类冰川模型可以应用在山地冰川的模拟上。为了比较不同模型的模拟结果， Pattyn等^［11］提出了高阶模型比较计划ISMIP-HOM，它针对已经提出的全Navier-Stokes方程以及高阶Navier-Stokes方程的冰川模型进行了模拟比较，为山地冰川在未来气候变化情景下消融进退的情况判断提供了定量客观的依据，并且有助于反演古气候变迁下冰川发生改变的情况。

我国的冰川研究还处在第二个阶段，即根据质量守恒方程，将观测数据进行诊断计算的物质平衡方法。这类方法定量化程度不高，预测未来冰川的变化需要气候模式的结果，不利于系统研究冰川和气候系统之间关系，也不能有效地和气候模式进行耦合。目前，中国尚未参与ISMIP-HOM计划。不过近年来，中国学者在冰川研究中取得了较好的进展^［15］，例如Zhang等^［16］用二维高阶热力流动带模型研究了珠穆朗玛峰的东绒布冰川的温度和速度分布； Wang等^［17］用原位测量与数值模拟的方法研究了祁连山老虎沟12号冰川的热力特征。这是可喜的进展，但仍然没有解决冰川剖分的问题。

我国是一个冰川资源丰富的大国。尤其青藏高原地区的冰川是我国西部地区和周边各国的主要水资源。而青藏高原地处中纬度高海拔地区，下垫面情况复杂，冰川高差较大，冰川的动力和热力过程极为复杂。为了提升我国气候预测能力，定量预估冰川水资源的变化和减轻自然灾害的能力，需要预估冰冻圈对未来气候变化的响应程度，从而迫切需要建立在观测基础上的冰川物理模型，进而为发展耦合冰冻圈各分量的气候系统模式提供必要的分量模式。

2 主要的冰川动力数值模式的建立及模拟

2.1　模型的一般形式

建立冰川模型时，如果忽略地球的曲率效应，引入正交直角坐标系，对于垂直部分，由于冰川的厚度远小于地球的半径，所以曲率效应可以忽略。在水平方向，冰川的扩展可以通过选择合适投影方法和垂直坐标表达（图3）。

图3

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

图3 冰川模型示意图（引自Dynamics of ice sheets andglaciers^［18］， 2009）

Fig.3 Schematic diagram of glacier model （cited from Dynamics of ice sheets and glaciers^［18］， 2009）

假设冰川为不可压缩，所以质量守恒方程为：

\nabla ∙ ν = 0

（1）

Glen流动定律则可以写做：

t^{D} = 2 η (T^{'}, d_{e}) D

（2）

式中： $t^{D}$ 为（Traceless）应力偏量； $η$ 为粘滞系数； $D$ 为应变率张量。

考虑Navier-Stokes方程：

ρ \frac{d ν}{d t} = - g r a d p + η Δ ν + (g r a d ν + {(g r a d ν)}^{T}) g r a d η + f

（3）

式中： $f$ 为地球引力和科氏力的合力：

f = ρ g - 2 ρ Ω \times ν

（4）

使用冰川的特征值 $[W] / [L]$ ，进行量纲分析，比较加速度项和气压梯度力项，得到加速度项可以忽略。适用同样的处理方法，分析冰川的Rossby数，科氏力项可以忽略。所以Navier-Stokes方程变为：

- g r a d p + η ∆ ν + (g r a d ν + {(g r a d ν)}^{T}) ∙ g r a d η + ρ g = 0

（5）

由于Navier-Stokes方程是一个速度场的微分方程，所以更方便使用依赖于速度梯度的切变黏度的公式。切边黏度的计算公式采用Glen型流定律，其指数n=3：

η (T^{'}, d_{e}) = \frac{1}{2} B (T^{'}) d_{e}^{- (1 - \frac{1}{n})}

（6）

$B (T')$ 是流动因子：

B (T') = A {(T')}^{- 1 / n}

（7）

$A (T')$ 是根据Hooke^［19］的研究结果用Arrhenius定律计算得到：

A (T') = A_{0} e^{- Q / R / (T + β p)}

（8）

式中： $β = 9.8 \times 10^{- 8} k P a^{- 1};$ $Q$ 为冰蠕动活化能； $R$ 为气体普适常数； $A_{0}$ 为常数。

由于黏性是依赖温度的，所以存在一个温度耦合的问题，其完全的方程组需要一个描述温度场的方程组：

\frac{d u}{d t} = c \frac{d T}{d t}, d i v q = - d i v (k g r a d T)

（9）

式中： $u$ 为内能； $q$ 为热通量，考虑到除了最上面的几厘米的冰层，辐射项可以忽略，所以可有演变方程：

ρ c \frac{d T}{d t} = d i v (k g r a d T) + 4 η d_{e}^{2}

（10）

假定冰川没有融化，所以有约束条件 $T \leq T_{m}$ 。将连续方程式（1），斯托克斯方程式（3），黏性计算公式（6）和温度演变方程式（9）综合起来，我们就得到了一个含有6个方程来求解6个未知数 $v_{x}, v_{y}, v_{z},$ $η, P, T$ 的封闭的热力斯托克斯方程组。

2.2　模型的简化

2.2.1　静水压力近似

由于冰川的 $z$ 平面的剪切应力 $τ_{x y}$ 和 $τ_{y z}$ 远小于 $τ_{z z}$ ，所以再次忽略切向应力，将式（5）垂直方向的方程简化成：

\frac{\partial t_{z z}}{\partial z} = ρ g

（11）

动量守恒方程变为：

$\frac{\partial t_{x x}}{\partial x} + \frac{\partial t_{x y}}{\partial y} + \frac{\partial t_{x z}}{\partial z} = 0$ ；

$\frac{\partial t_{y x}}{\partial x} + \frac{\partial t_{y y}}{\partial y} + \frac{\partial t_{y z}}{\partial z} = 0$ ；

\frac{\partial t_{z z}}{\partial z} = ρ g

（12）

经整理基本动量方程式（5）可以写作：

$2 \frac{\partial t_{x x}^{D}}{\partial x} + \frac{\partial t_{y y}^{D}}{\partial x} + \frac{\partial t_{x y}^{D}}{\partial y} + \frac{\partial t_{x z}^{D}}{\partial z} = ρ g \frac{\partial h}{\partial x}$ ；

2 \frac{\partial t_{y y}^{D}}{\partial y} + \frac{\partial t_{x x}^{D}}{\partial x} + \frac{\partial t_{x y}^{D}}{\partial x} + \frac{\partial t_{y z}^{D}}{\partial z} = ρ g \frac{\partial h}{\partial y}

（13）

Glen流公式（2）则简化成了：

$t_{x x}^{D} = 2 η \frac{\partial v_{x}}{\partial x}$ ； $t_{y y}^{D} = 2 η \frac{\partial v_{y}}{\partial y}$ ； $t_{z z}^{D} = 2 η \frac{\partial v_{z}}{\partial z}$ ；

$t_{x z} = η (\frac{\partial v_{x}}{\partial z} + \frac{\partial v_{z}}{\partial x})$ ； $t_{y z} = η (\frac{\partial v_{y}}{\partial z} + \frac{\partial v_{z}}{\partial y})$ ； $t_{x y} = η (\frac{\partial v_{x}}{\partial y} + \frac{\partial v_{y}}{\partial x})$

（14）

在这样的情况下的数值模型是静力近似下的数值模型，虽然这样的方程依然很复杂，但是相对于原来的Navier-Stokes方程，已经做了很大的近似。并且这样的近似对于大多数空间尺度的冰川是合理的，所以这样的近似也就成了其他冰川模型的一条基本假设（对于山地冰川而言，由于有槽谷的影响，在靠近冰川的边缘，这样的假设不一定成立）。但是由于这样的简化并没有有效地简化方程的复杂性，所以这样的模式并不实用，需要进一步简化。

2.2.2　一阶近似

在浅冰近似的基础上，由于垂直速度的水平方向的变化量（ $d w / d x ≪ d u / d z$ ）远小于水平速度在垂直方向上的偏导，所以对于黏性方程，可以做进一步的简化。对于Glen流公式有：

\begin{array}{l} t_{x x}^{D} = 2 η \frac{\partial v_{x}}{\partial x}; t_{y y}^{D} = 2 η \frac{\partial v_{y}}{\partial y}; t_{z z}^{D} = 2 η \frac{\partial v_{z}}{\partial z}; \\ t_{x z} = η \frac{\partial v_{x}}{\partial z}; t_{y z} = η \frac{\partial v_{y}}{\partial z}; t_{x y} = η (\frac{\partial v_{x}}{\partial y} + \frac{\partial v_{y}}{\partial x}) \end{array}

（15）

这样的近似方案可以对方程进一步简化，减少计算量，同时这样的尺度分析对于各种尺度的冰川都是有效的，又由于基本上保持了Navier-Stokes方程的特征，所以这样的近似结果可以作为别的数值模式模拟结果的检验标准。

2.2.3　浅冰近似

目前流行的模型多是基于浅冰近似模型。浅冰模型（图4）是现代冰川模型的基础，其完整的数学物理基础模型是由Hutter^［20］于1983年提出的。由于求解全部的方程组需要大量的数值计算，为了节省计算，提出了浅冰模型。早期的浅冰模型只考虑了各向同性的情况，后来增加了各向异性的情况的分析^［21］，并且耦合了温度场。

图4

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

图4 浅冰近似示意图（引自Dynamics of ice sheets andglaciers^［18］， 2009）

Fig.4 Shallow ice approximation （cited from Dynamics of ice sheets and glaciers^［18］， 2009）

当考虑到水平尺度较小的冰川如山地冰川的时候，浅冰模型的尺度分析基础不再适合于这样的情形，于是人们开始对完全的Navier-Stokes方程构造的数值模式模拟的结果和浅冰模型的模拟结果进行比较^［22］。结果表明，对于浅冰模型影响较大的是底面的坡度、冰川底部的槽谷。当底面的坡度较大的时候，浅冰模型模拟的结果也较差。同时也发现，降雪等气候原因造成的冰川厚度的增加对于最后模拟结果的影响不大。

对于山地冰川，气候变动可以导致他们形态和动力结果的改变，目前有研究表明，山地冰川可以很好地反映出气候的变化^［23］。由于冰川观测的困难性，数值模式模拟被认为是分析冰川对气候复杂响应的最主要、甚至是唯一的方式和手段。

对于大尺度的冰盖，除了在冰穹附近的区域（大约在水平10 km范围之内）以及冰的边缘区域，冰的流场是近似平行于底面的。同时，自由面和冰的底面间的坡度也很小。在这种情况下，通过尺度分析，法向应力的偏移量和垂直平面的切向应力可以忽略，所有法向应力则等于负的压力^［20］：

t_{x x} = t_{y y} = t_{z z} = - p

（16）

所以垂直方向的动量平衡就简化成：

\frac{\partial p}{\partial z} = - ρ g

（17）

对它积分可以得到静止压力的积分：

p = p_{h y d} = ρ g (h - z)

（18）

水平方向的动量平衡就简化成：

$\frac{\partial t_{x z}}{\partial z} = \frac{\partial p}{\partial x} = ρ g \frac{\partial h}{\partial x}$ ；

\frac{\partial t_{y z}}{\partial z} = \frac{\partial p}{\partial y} = ρ g \frac{\partial h}{\partial y}

（19）

考虑到地面的坡度非常小，可以忽略，所以表面的垂直法向量就可以简化成只有垂直方向的。自由界面的条件就可以简化成：

{p|}_{z = h} = 0; {t_{x z}|}_{z = h} = 0; {t_{y z}|}_{z = h} = 0

（20）

根据以上条件，积分式（19）就可以得到：

t_{x z} = - ρ g (h - z) \frac{\partial h}{\partial x}; t_{y z} = - ρ g (h - z) \frac{\partial h}{\partial y}

（21）

式（19）和（21）表明在浅冰近似里，应力场完全依赖于冰川的几何形态。已知有效应力则等于：

σ_{e} = \sqrt[]{t_{x z}^{2} + t_{y z}^{2}} = - ρ g (h - z) |g r a d h|

（22）

将以上结果代入Glen公式：

D = A (T') σ_{e}^{n - 1} t_{D}

（23）

则有：

\frac{1}{2} (\frac{\partial v_{x}}{\partial z} + \frac{\partial v_{z}}{\partial x}) = A (T^{'}) σ_{e}^{n - 1} t_{x z} =

$- A (T^{'}) [ρ g {(h - z)}^{n}] {|g r a d h|}^{n - 1} \frac{\partial h}{\partial x}$ ；

\frac{1}{2} (\frac{\partial v_{y}}{\partial z} + \frac{\partial v_{z}}{\partial y}) = A (T^{'}) σ_{e}^{n - 1} t_{y z} =

- A (T^{'}) [ρ g {(h - z)}^{n}] {|g r a d h|}^{n - 1} \frac{\partial h}{\partial y}

（24）

根据量纲分析的结果，垂直速度的水平梯度被忽略，这样就有：

$\frac{1}{2} \frac{\partial v_{x}}{\partial z} = A (T^{'}) σ_{e}^{n - 1} t_{x z} = - A (T^{'}) [ρ g {(h - z)}^{n}] {|g r a d h|}^{n - 1} \frac{\partial h}{\partial x}$ ；

\frac{1}{2} \frac{\partial v_{y}}{\partial z} = A (T^{'}) σ_{e}^{n - 1} t_{y z} = - A (T^{'}) [ρ g {(h - z)}^{n}] {|g r a d h|}^{n - 1} \frac{\partial h}{\partial y}

（25）

对质量守恒方程式（1）积分，可以得到冰川厚度方程：

\frac{\partial}{\partial x} \int_{b}^{h} v_{x} d z + \frac{\partial}{\partial y} \int_{b}^{h} v_{y} d z + \frac{\partial h}{\partial t} - N_{s} a_{\frac{1}{s}} - \frac{\partial b}{\partial t} + N_{b} a_{\frac{1}{b}} = 0

（26）

对温度场进行类似的量纲分析，可以得到：

ρ c \frac{\partial T}{\partial t} = \frac{\partial}{\partial z} (k \frac{\partial T}{\partial z}) + 4 η A^{2} (T^{'}) σ_{e}^{2 n}

（27）

其中 $σ_{e}$ 由式（22）得到。

式（21）代表了浅冰近似中平行于底面的应力，计算这个公式在底部的结果可以得到向量：

τ_{d} = - ρ g H (\binom{\frac{\partial h}{\partial x}}{\frac{\partial h}{\partial y}})

（28）

在浅冰近似里，拖曳应力等于负的底部拖拽，也即假设驱动力和阻抗力是平衡的。

目前有许多冰川（盖）模式模拟冰川（盖）的变化，其中有代表性的冰盖模式是GLIMMER数值模式， GLIMMER模式是浅冰理论模型的应用。

3 GLIMMER模式

由图5可知， GLIMMER 模式是一个三维有限差分冰盖模式，它源于Payne^［24］。该模式在研究GENIE地球系统模式时，对南极陆冰部分模拟工作。GLIMMER模式的目标是要建立一个可供其他地球系统模拟模式调用的标准模式。GLIMMER模式通过Fortran95程序库实现，可以被其他赋有边界条件的模式调用，它还包含从简单的EISMINT类型到耦合GENIEESM的复杂驱动器的驱动项。输入与输出数据通过netCDF格式编写。在2006年，该模式已经通过在EISMINT21和EISMINT22基准下进行了严格测试^［25］。

图5

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

图5 GLIMMER 模型示意图（引自GLIMMER 1.5.1 Documentation， 2010）

Fig.5 Schematic diagrams of the GLIMMER Model： Cartesian coordinates（a）； Sigma coordinate system（b）（cited from GLIMMER 1.5.1 Documentation， 2010）

该模式是基于浅冰模型的理论，它包括几个基本部分：（1）GLIDE（General Land Ice Dynamic Elements）部分是整个模型的核心。这个部分包括一个冰盖模型，它负责计算冰的速度，内部冰的温度分布，各向异性的调整，以及融化水的产生。这一模型需要气候模式的驱动才能运行。（2）SIMPLE：是一个简单的气候驱动的界面程序，以适应EISMINIT计划，并提供了理想的地理条件。（3）GLINT：是GLIMMER的界面程序，原来为了GENIE3地球系统所设计。（4）EIS：这是Edinburgh冰盖气候驱动项，它基于对平衡线高程、海平面表面温度和海平面变化作参数化处理。上述后面二者是气候驱动项，另外模式还包括一个可被其他部分调用的多功能模块GLUM和一些通用的可视化程序。同时GLIMMER还包括一个补充的质量平衡参数化方案。

显然，该模型虽然是冰盖模型，但可以经过改进应用到研究单条冰川模型的变化，其结果直观，操作简单。然而，该模型不能用于冰盖（川）对大气的反馈研究。

4 山地冰川模型

针对山地冰川，传统的浅冰近似已近无法做出合理的模拟描述，需要高阶的模型或者是全Navier-Stokes方程来模拟。

4.1　模型的动力框架

所谓冰盖或山地冰川高阶模型是指考虑除水平方向以外的其他法向应力梯度效应的模型^［26］。其模型依据质量守恒公式（1）和动量守恒公式（29），如下：

ρ \frac{d v}{d t} = \nabla ∙ T + ρ g

（29）

式中： $ρ$ 代表冰的密度； $g$ 代表重力加速度； $v$ 为速度矢量； $T$ 为应力张量，相应的常数见表1。

表1 冰川常用物理参数

Table 1 Common physical parameters of glaciers

参数	量级	单位
$A$ （冰流常数）	10^-16	Pa^-n·a^-1
$ρ$ （冰密度）	910	kg·m^-3
$g$ （重力加速度）	9.81	m·s^-2
$N$ （Glen流定律常数）	3

新窗口打开| 下载CSV

通过使用冰川的特征值，进行量纲分析，比较加速度项和气压梯度力项，得到加速度项可以忽略。采用同样的处理方法，分析冰川的Rossby数，科氏力项可以忽略。

由于是不可压缩的流体，所以应力张量可以分解成偏应力部分和各向同性的压力部分：

T = T^{'} - P I

（30）

其中的对于冰雪的黏性偏应力和应变率的本构方程：

T^{'} = 2 η \dot{e}

（31）

式中： $T^{'}$ 为偏应力； $\dot{e}$ 为应变率； $η$ 为有效黏性。线性和非线性的冰川流体都被考虑了进来。现在普遍使用的是Glen流定律：

η = \frac{1}{2} A^{- 1 / n} {\dot{ε}}_{e}^{(1 - n) / n}

（32）

式中： ${\dot{ε}}_{e}$ 为应变率的第二不变量，对于线性的情形（即牛顿流体），方程取 $n$ =1， $A$ 变成了整个冰盖区域都为常数的情况，忽略加速度项，则整个动量方程可以写作：

d i v T + ρ_{i} g = d i v T^{'} - g r a d P + ρ g = 0

（33）

由于只有重力加速度项这一项外力作用，所以有：

$\frac{\partial t_{y x}}{\partial x} + \frac{\partial t_{y y}}{\partial y} + \frac{\partial t_{y z}}{\partial z} = 0$ ；

$\frac{\partial t_{z z}}{\partial z} = ρ g$ ；

\frac{\partial t_{x x}}{\partial x} + \frac{\partial t_{x y}}{\partial y} + \frac{\partial t_{x z}}{\partial z} = 0

（34）

4.2　上边界条件

类似于所有边界条件，冰川的自由边界可以认为是一个奇异边界，如果我们将冰川表面的梯度向量代表成：

n = \frac{1}{N_{s}} (\begin{matrix} \binom{- \frac{\partial h}{\partial x}}{- \frac{\partial h}{\partial y}} \\ - 1 \end{matrix})

（35）

其中指向大气为正方向。如果我们用隐式表达奇异面：

F_{s} (x, t) = z - h (x, y, t) = 0

（36）

由于其局地导数为0，所以有：

\frac{\partial h}{\partial t} + v_{x} \frac{\partial h}{\partial x} + v_{y} \frac{\partial h}{\partial y} - v_{z} = N_{s} a_{\frac{1}{s}}

（37）

其中 $a_{\frac{1}{s}}$ 是垂直方向的由于降水导致的冰川的累积。

在动力方面，由于自由面不受力，所以有零应力条件：

t ∙ n = 0

（38）

在热力方面，需要和自由大气一致，即冰雪的温度和大气的温度是一样的：

T = T_{s}

（39）

通过数值模式实验也证明了，这样的近似对于0 ℃以下的大气温度是一个很好的近似。

4.3　下边界条件

类似于上边界的奇异面，在下边界也有相同的边界条件：

\frac{\partial h}{\partial t} + v_{x} \frac{\partial h}{\partial x} + v_{y} \frac{\partial h}{\partial y} - v_{z} = N_{s} a_{\frac{1}{b}}

（40）

其中 $a_{\frac{1}{b}}$ 是地热等作用导致冰川底部的消融，其计算式为：

a_{\frac{1}{b}} = \frac{q^{\frac{1}{g e o}} - k (g r a d T ∙ n) - v_{b} ∙ t ∙ n}{ρ L}

（41）

对于冰川底部的动力条件是：

t ∙ n = t_{l i t h} ∙ n

（42）

这说明冰川底部的应力条件是连续的，但很显然，冰川底部应力很难测得（需要打探钻孔）。所以又有了相应的经验的参数化方案来解决冰川底部应力的给定^［18］，在已经求得冰川底部温度的情况下可以使用如下的经验公式来给出冰川的边界条件：

v_{b} = \{\begin{array}{l} 0, & i f T_{b} < T_{m} \\ x - C_{b} \frac{τ_{b}}{N_{b}^{q}} e_{t}, & i f T_{b} = T_{m} \end{array}

（43）

此为数学上的简化，假设底部的温度条件是已知常数，其中e_t 表示在切向于冰川底部平面基底的剪切应力的方向（e_t ⊥n）。

5 山地冰川模拟存在的问题及一个新模型的设想

目前，我国在山地冰川模型的应用和建模方面开展了初步的探索工作^［15-16］。如前所述，总体上还在使用公用计算工具（如COMSOL）对冰川体进行剖分，然后再计算其固体本构方程。这种方法有如下限制： ①该途径对单个冰川有效； ②这种途径不能和其他模式进行耦合，其剖分结果也取决于边界条件。因此，尚不能成为与气候模式实现耦合的“模式”，也就不能实现和大气之间的相互作用的研究。因此，我们试图通过以下思路，建立一个基于有限体积法剖分求解Navier-Stokes方程、可以任意和其他模式灵活耦合、改变边界条件的建模思路。

第一，在空间尺度上，山地冰川远小于冰盖，这导致在物理本质上浅冰模型并不适合于山地冰川。现在也有一些工作在研究浅冰模型可以在不同程度上对山地冰川能进行近似模拟，但整体上，浅冰模型并不适宜使用在山地冰川上。因此，全Navier-Stokes方程的求解是山地冰川模型建立的首选。

第二，由于山地冰川的底部状况不同于冰盖模型，可以认为是一个坡度可忽略的平面，其底部的地形结构对于山地冰川的模拟有很大的影响。另一方面，由于野外工作的艰巨性，仅有部分是通过雷达声纳技术得到的山地冰川底部的地形资料，这就严重的影响了山地冰川数值模式的下边界条件的给定。现有的方式是将下边界条件的求取当做一个反问题，先假设一个下边界条件，通过将模拟结果和实测数据的比较，再逐步订正下边界条件^［18］。但由于冰川是一个很好的高频滤波器，底部的许多短波起伏数据是无法用这一方法来得到的。因此，我们拟用临近边界值和冰川中轴线的线性变化反演底部边界条件。

第三，山地冰川内部的速度流场很难给定，常规的打孔测量数据对于为数众多的山地冰川并不合适，实际上只有极少数的山地冰川有速度流场资料而且精度较低，这给山地冰川数值模式的初始化提出了问题，同时模拟结果的校验也比较困难。因此，在一级近似的基础上，针对山地冰川的特点建立一个计算量适中的有限体积上的数值模式，设计出一整套合理的求取下边界条件这一反问题的数值模拟方案，合并到整个模式之中。

第四，除了动力方面，在热力方面，山地冰川和冰盖也有很大的区别，大多数冰盖中的冰是冷冰，即其温度在局地压力下的熔点以下，而温度达到熔点的暖冰仅存在于冰川的底部和较薄的冰层附近，或者存在于积累区或者消融区的表层。对于较高纬度的山地冰川，这样的问题同样存在。然而，对于许多低纬度的山地冰川，它们相当一部分是由暖冰组成。暖冰含有一小部分的液态水，因此热量部分的变化将会导致水含量的变化，相对于冷冰而言，热量部分的变化仅仅会导致温度的变化，而暖冰的热力动力学性质都会随着含水量的多少而改变，为全面地模拟出这种流场的特性，有必要对于水含量的空间分布有所掌握。在此基础上，已有一些研究提出了针对冰盖的数值模式以模拟含有暖冰的冰川^［27］。通过冰川底部基岩的均一性，建立底部边界条件中热通量的边界值。

第五，冰川模型的理论已基本完备，其中方程的计算格式也有比较成熟的稳定和守恒格式，如何将一个不规则形、热固耦合的物理过程计算出来，核心问题之一是对冰川的剖分。对Navier-Stokes方程在山地冰川的模拟，我们初步提出了一个有限体积法进行剖分的思路，比较了六面体和四面体的剖分结果，后者似乎比前者具有优势，通过平面高度恰当地找到剖分点单元和单元数据^［28］。

根据上述思路，我们建立了一个山地冰川的模型（图6）并开展了初步试验。

图6

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

图6 山地冰川模型示意图

Fig.6 Schematic diagrams of mountain glaciers， projected to： xz plane（a）； yz plane（b）； xy plane（c）

我国是中、低纬度地区冰冻圈最发达的国家，冰川面积达59 425 km²，占全球中、低纬度冰川面积的50%以上；冰川冰储量5 600 km³。冰川、积雪年融水量达1 360×10⁸ m³。近百年来，我国冰冻圈显著萎缩，已对区域气候、水资源、生态与环境产生了重大影响，在未来全球气候变暖背景下，随着冰冻圈萎缩加剧，对气候和环境的影响也将更为突出。冰冻圈变化对我国的水资源安全，维系我国西部高寒和干旱区生态系统稳定的基本保障，对我国西部生态安全的威胁日益凸显。我国作为中低纬度冰冻圈最发达的国家，其变化对我国及周边地区气候有重要影响。因此，在对冰川进行加密观测和监测的同时，研发、建立具有自主知识产权的冰川模式具有紧迫性。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Waddington

E D

， Clarke

G K C

Stable-isotope pattern predicted in surge-type glaciers

［J］. Canadian Journal of Earth Sciences， 1988， 25（5）： 657 - 668.