冰川荷载作用下地壳回弹问题的混合间断有限元

Mixed discontinuous Galerkin element for crustal rebound problem caused by glacier load

摘要

随着全球变暖加剧，冰川的持续退缩造成冰川覆面地壳回弹及冰缘区崩塌滑坡灾害频发。在降水作用下，这些灾害可能进一步转化为破坏性更强的地质灾害链。重调和方程常被用于建模冰川载荷作用下的地壳回弹问题，但主流的数值解法难以兼顾高精度和低复杂度。本文采用混合间断有限元方法对该问题进行数值求解，通过引入中间变量降低方程对有限元空间光滑性的要求，且间断有限元的选取比传统有限元更具灵活性，能在保持计算精度的前提下极大降低计算复杂度。本文以中国最大的海洋性冰川——恰青冰川为例研究了冰川完全消融时地壳的回弹形变规律。参数敏感性分析表明，地壳厚度对回弹形变的影响较大，而地壳岩石的杨氏模量对回弹形变的影响则相对较小。本文的结果为冰川消融引发滑坡灾害的预测与防治提供了重要依据。

Abstract

With the intensification of global warming， continuous melting of overlying glaciers leads to crustal rebound and results in the development of collapse and landslide disasters in the periglacial area.These disasters can further transform into more destructive geological hazard chains under the influence of precipitation.The biharmonic equation is often used to model the crustal rebound problem caused by glacier load.Nowadays， mainstream numerical methods are difficult to balance accuracy and complexity.To address this problem， a mixed discontinuous finite element method is adopted by introducing an intermediate variable to reduce the requirement of smoothness of finite element space for the biharmonic equation.In comparison with the traditional finite element methods， the selection of discontinuous Galerkin element has greater flexibility and can greatly simplify the complexity of numerical simulation， while maintaining high calculation accuracy.Furthermore， this study simulates the rebound of crust when the largest marine glacier in China， the Qiaqing Glacier， completely melts.Sensitivity analysis indicates that the thickness has a significant impact on the rebound deformation and the Young’s modulus of crust is relatively minor.The obtained results are expected to provide an important reference for the prediction and prevention of landslide disasters triggered by glacier melt.

Graphical abstract

关键词

冰川荷载 / 地壳回弹 / 重调和方程 / 混合间断有限元

Key words

glacier loading / crustal rebound / biharmonic equation / mixed discontinuous Galerkin element（2020 MSC 65M60）

青藏高原现今的地壳隆起状态是地球科学的一个重要关注点，其中由地表负荷改变引起的地壳形变会直接体现在青藏高原的地壳隆升观测中^［1-3］。青藏高原有着中低纬度地区最大的现代冰川分布^［4-5］，全球变暖明显造成了冰川退缩及地壳回弹^［6-7］，导致地震活动发生的可能性增加^［8-9］。地壳回弹还可能导致崩塌滑坡等地质灾害^［10］，这些灾害在降水的作用下会进一步形成泥石流灾害链。因而冰川退缩引起的地壳回弹问题一直是学界关注的焦点。

通常而言，地球可以看作是一种弹性体，地球表面负荷改变引起的地壳形变主要发生在垂直方向^［11］。根据地壳均衡理论，常将地壳看作弹性薄板，用四阶重调和方程来描述其受载形变^［12］。由于重调和方程的解析求解十分困难，差分法^［13-14］、有限体积法^［15］、半解析法^［16］、谱方法^［17］以及有限元法^［18-19］等数值解法成为了目前求解重调和方程的主流方法。其中，有限差分法具有方法简单、容易实现等优点，但对边界区域的规则性要求很高，难以处理复杂网格，并且可能存在精度不高的问题。有限体积法虽然自身包含几何信息，易于处理复杂网格，但计算复杂度高，并且不易提高精度。谱方法将解近似展开成光滑函数，求解精度较高，但存在伪震荡现象，不适用于复杂几何模型。有限元法是目前应用最多的方法，对不规则区域边界的处理较为方便，有着求解精度高等优势。

在四阶重调和方程的有限元求解中，协调元的构造需要具有C^{\mathrm{1}}-连续性，如Argyris元^［20］ 和Bogner-Fox-Schmidt元^［21］。一般而言，重调和方程协调元的构造较为困难^［22］，且实现过程中容易出现计算非常密集的问题。混合有限元方法是一种处理四阶问题的经典方法，通过引入一个中间变量将四阶重调和问题转换成二阶系统，在此基础上建立混合变分方程。作为传统有限元方法的改进，间断有限元（Discontinuous Galerkin，DG）方法允许有限元函数在单元界面上不连续，有着构造简单的优势。鉴于混合有限元方法和间断有限元方法求解重调和方程的便利性，Gudi等^［23］利用混合间断有限元方法求解了重调和方程，Feng等^［24］利用混合间断有限元方法求解了重调和方程的特征值问题。

本文参考文献［23］中的变分格式，对之进行改进，采用混合间断有限元方法对四阶重调和型方程进行数值求解，并利用数值算例对算法进行了验证。在此基础上，以恰青冰川消融为例模拟了冰川消融后的地壳回弹程度，为进一步的地质研究提供一定的参考。

1 混合间断有限元离散

地壳均衡是地球科学的一个基本概念。经典Vening Meinesz区域均衡理论^［25］将地壳视作一个漂浮在流塑地幔之上的弹性薄板，地壳的形变可以由弹性薄板的形变来表示，薄板在上层地质荷载和下层软流层的补偿作用下达到力学平衡。

1.1 静压平衡下的地壳弯曲方程

将地壳近似为一个弹性薄板，基于Kirchhoff-Love假设，可以将三维板壳模型转换为二维平面问题。由均衡理论可知，弹性薄板弯曲问题的控制方程^［26］为：

∆D\left(∆W\right)=P（1）

其中，算子∆=\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}};W为垂直位移，取正向向下；x和y为水平坐标；P为单位面积的荷载合力；D为板的抗弯刚度。

D=\frac{E{h}^{\mathrm{3}}}{\mathrm{12}\left(\mathrm{1}-{\nu }^{\mathrm{2}}\right)}（2）

式中，h为地壳厚度，E为地壳的杨氏模量，\nu为地壳的泊松比。此外，荷载合力为：

P=P_{\mathrm{0}}-\left({\rho }_m-{\rho }_c\right)gW（3）

其中，P_{\mathrm{0}}为和位移无关的地壳上方的地质荷载，g为重力加速度，{\rho }_m为地幔密度，{\rho }_c为地壳密度。若取k=\left({\rho }_m-{\rho }_c\right)g，结合方程（1）和（3）可得：

∆D\left(∆W\right)+kW=P_{\mathrm{0}}（4）

这是一个变系数的重调和型方程。为方便处理，取D(x,y)为计算区域平均的地壳厚度下的抗弯刚度，即D\left(x,y\right)=D为一个常数。此时，可将方程（4）改写成：

{∆}^{\mathrm{2}}W+\frac{k}{D}W=\frac{\mathrm{1}}{D}{P}_{\mathrm{0}}（5）

其中，k/D和\mathrm{1}/D均为常数。

1.2 预备知识

假设在区域\mathrm{\Omega }上作正则剖分T_h=\{K_i,\mathrm{1}\le i\le N\}，其中N为区域\mathrm{\Omega }的剖分网格数。定义区域边界上的单元边界为{\mathrm{\Gamma }}_D=\{e_k:e_k=\partial K_i\bigcap \partial \mathrm{\Omega }\}，内部单元边界为{\mathrm{\Gamma }}_I=\{e_k:e_k=\partial K_i\cap \partial K_j,\mathrm{1}\le i,j\le N\}。令\mathrm{\Gamma }={\mathrm{\Gamma }}_I\bigcup {\mathrm{\Gamma }}_D为所有单元边界的集合。定义H^{s_i}(K_i)为单元K_i上标准Sobolev空间，H^s\left(\mathrm{\Omega },T_h\right)=\{u\in L^{\mathrm{2}}\left(\mathrm{\Omega }\right):{\left.u\right|}_{K_i}\in H^{s_i}\left(K_i\right),\forall K_i\in T_h\}，其中s=\{s_i\ge \mathrm{0}:i=\mathrm{1},\dots ,N\}。对任意u\in H^s\left(\mathrm{\Omega },T_h\right)，分别定义单元边界 e_K 上的跳量和平均量如下：

1） 如果e_K=\partial K_i\cap \partial K_j\left(ij,\mathrm{1}\le i,j\le N\right)\in {\mathrm{\Gamma }}_I，则

\left[u\right]={\left.u\right|}_{K_i}-{\left.u\right|}_{K_j}（6）

\left\{u\right\}=\frac{{\left.u\right|}_{K_i}+{\left.u\right|}_{K_j}}{\mathrm{2}}（7）

2） 如果e_K=\partial K_i\bigcap \partial \mathrm{\Omega }\in {\mathrm{\Gamma }}_D， 则

\left[u\right]=\left\{u\right\}={\left.u\right|}_{K_i\bigcap \partial \mathrm{\Omega }}（8）

为了进行误差分析，定义DG空间中的范数和半范数如下：

{‖u‖}_{H^s\left(\mathrm{\Omega },T\right)}={\left({\displaystyle \sum _{i=\mathrm{1}}^N}{‖u‖}_{H^{S_i}\left(K_i\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}（9）

{\left|u\right|}_{H^s\left(\mathrm{\Omega },T\right)}={\left({\displaystyle \sum _{i=\mathrm{1}}^N}{\left|u\right|}_{H^{S_i}\left(K_i\right)}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}（10）

1.3 混合间断有限元格式

为了便于有限元方法实现和简化复杂工程问题的计算，本文用混合间断有限元方法求解如下带边界条件的控制方程：

\left\{\begin{array}{l}{∆}^{\mathrm{2}}W+\frac{k}{D}W=\frac{\mathrm{1}}{D}{P}_{\mathrm{0}},\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ W=g_{\mathrm{1}},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega },\\ \nabla W\cdot n=g_{\mathrm{2}},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\end{array}\right.（11）

其中，\mathrm{\Omega }\in {\mathbf{R}}^{\mathrm{2}}为一个有界的凸区域，\partial \mathrm{\Omega }为区域\mathrm{\Omega }的光滑边界，n为边界\partial \mathrm{\Omega }上的单位外法向量。 方程（11）为一个四阶问题，为达到降阶目的，引入中间变量Q=∆W，可将方程（11）改写成：

\left\{\begin{array}{l}Q=∆W,\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ ∆Q+\frac{k}{D}W=\frac{\mathrm{1}}{D}{P}_{\mathrm{0}},\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ W=g_{\mathrm{1}},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega },\\ \nabla W\cdot n=g_{\mathrm{2}},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\end{array}\right.（12）

并假设函数P_{\mathrm{0}}\mathrm{、}{g}_{\mathrm{1}}\mathrm{、}{g}_{\mathrm{2}}都足够光滑，以确保方程（12）有唯一解W\in H^{\mathrm{4}}。

取V=H^{\mathrm{2}}(\mathrm{\Omega },T_h)。对方程（12）中第一式左右两边同时乘以u\in V，并在区域\mathrm{\Omega }上积分，由格林公式知：

{\displaystyle \sum _i^N}\underset{K_i}{\int }\nabla W\cdot \nabla u\mathrm{d}x-\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\left\{\frac{\partial W}{\partial n}\right\}\left[u\right]\mathrm{d}s+\underset{\Omega }{\int }Qudx=\mathrm{0}（13）

由边界条件\nabla W\cdot n=g_{\mathrm{2}}，上式可以改写成：

{\displaystyle \sum _{i=\mathrm{1}}^N}\underset{K_i}{\int }\nabla W\cdot \nabla u\mathrm{d}x-\sum _{e_K\in {\mathrm{\Gamma }}_{\mathrm{I}}}\underset{e_K}{\int }\left\{\frac{\partial W}{\partial n}\right\}\left[u\right]\mathrm{d}s+\underset{\mathrm{\Omega }}{\int }Qu\mathrm{d}x=\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }{g}_{\mathrm{2}}u\mathrm{d}s（14）

由于在内部边界{\mathrm{\Gamma }}_{\mathrm{I}}上有\left[W\right]=\mathrm{0}，由边界条件W=g_{\mathrm{1}}，有

-\sum _{e_K\in \Gamma }\underset{e_K}{\int }\left\{\frac{\partial u}{\partial n}\right\}\left[W\right]\mathrm{d}s=-\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }\frac{\partial u}{\partial n}{g}_{\mathrm{1}}\mathrm{d}s（15）

由（14）和（15）式可知：

{\displaystyle \sum _{i=\mathrm{1}}^N}\underset{K_i}{\int }\nabla W\cdot \nabla u\mathrm{d}x-\sum _{e_K\in {\mathrm{\Gamma }}_{\mathrm{I}}}\underset{e_K}{\int }\left\{\frac{\partial W}{\partial n}\right\}\left[u\right]\mathrm{d}s-\sum _{e_K\in \Gamma }\underset{e_K}{\int }\left\{\frac{\partial u}{\partial n}\right\}\left[W\right]\mathrm{d}s+\underset{\mathrm{\Omega }}{\int }Qu\mathrm{d}x=\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }{g}_{\mathrm{2}}u\mathrm{d}s-\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }\frac{\partial u}{\partial n}{g}_{\mathrm{1}}\mathrm{d}s（16）

同理，对方程（12）中的第二式左右两边同时乘v\in V并在\mathrm{\Omega }上积分，有

{\displaystyle \sum _{i=\mathrm{1}}^N}\underset{K_i}{\int }\nabla Q\nabla v\mathrm{d}x-\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\left\{\frac{\partial Q}{\partial n}\right\}\left[v\right]\mathrm{d}s-\underset{\mathrm{\Omega }}{\int }\frac{k}{D}Wv\mathrm{d}x=-\underset{\mathrm{\Omega }}{\int }\frac{\mathrm{1}}{D}{P}_{\mathrm{0}}v\mathrm{d}x（17）

在内部边界{\mathrm{\Gamma }}_{\mathrm{I}}上有\left[Q\right]=\mathrm{0}，于是

\sum _{e_K\in {\mathrm{\Gamma }}_I}\underset{e_K}{\int }\left\{\frac{\partial v}{\partial n}\right\}\left[Q\right]\mathrm{d}s=\mathrm{0}（18）

\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\alpha \left[W\right]\left[v\right]\mathrm{d}s=\sum _{e_K\in {\mathrm{\Gamma }}_{\mathrm{D}}}\underset{e_K}{\int }\alpha g_{\mathrm{1}}v\mathrm{d}s（19）

其中，\alpha为惩罚参数，可取{\left.\alpha \right|}_{e_K}={\alpha }_K{\left|e_K\right|}^{-\mathrm{3}}{p}_K^{\mathrm{2}},{\alpha }_K为稳定参数；p_K为有限元空间多项式的次数；|e_K|为单元所在边的长度。由（17）~（19）式可知：

\begin{array}{l}{\displaystyle \sum _{i=\mathrm{1}}^N}\underset{K_i}{\int }\nabla Q\nabla v\mathrm{d}x-\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\left\{\frac{\partial Q}{\partial n}\right\}\left[v\right]\mathrm{d}s-\underset{\mathrm{\Omega }}{\int }\frac{k}{D}Wv\mathrm{d}x-\sum _{e_K\in {\mathrm{\Gamma }}_I}\underset{e_K}{\int }\left\{\frac{\partial v}{\partial n}\right\}\left[Q\right]\mathrm{d}s-\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\alpha \left[W\right]\left[v\right]\mathrm{d}s=-\underset{\mathrm{\Omega }}{\int }\frac{\mathrm{1}}{D}{P}_{\mathrm{0}}v\mathrm{d}x-\\ \sum _{e_K\in {\mathrm{\Gamma }}_{\mathrm{D}}}\underset{e_K}{\int }\alpha g_{\mathrm{1}}v\mathrm{d}s\end{array}（20）

令有限元解空间V_h=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{{\psi }_{\mathrm{1}},\dots ,{\psi }_{N_b}\}，其中N_b=\mathrm{d}\mathrm{i}\mathrm{m}(V_h)。令W_h={\displaystyle \sum _{i=\mathrm{1}}^{N_b}}{X}_i{\psi }_i为W的有限元解，Q_h={\displaystyle \sum _{i=\mathrm{1}}^{N_b}}{Y}_i{\psi }_i为Q的有限元解。又令X={\left(X_{\mathrm{1}},X_{\mathrm{2}},\dots ,X_{N_b}\right)}^{\mathrm{T}},Y={\left(Y_{\mathrm{1}},Y_{\mathrm{2}},\dots ,Y_{N_b}\right)}^{\mathrm{T}}，并给出如下记号：

A_{lm}=\sum \underset{K_i}{\int }\nabla {\psi }_l\cdot \nabla {\psi }_m\mathrm{d}x-\sum _{e_K\in {\mathrm{\Gamma }}_I}\underset{e_K}{\int }\left\{\frac{\partial {\psi }_m}{\partial n}\right\}\left[{\psi }_l\right]\mathrm{d}s-\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\left\{\frac{\partial {\psi }_l}{\partial n}\right\}\left[{\psi }_m\right]\mathrm{d}s（21）

M_{lm}=\underset{\mathrm{\Omega }}{\int }{\mathrm{\psi }}_l{\mathrm{\psi }}_m\mathrm{d}x（22）

{\tilde{M}}_{lm}=\frac{k}{D}\underset{\mathrm{\Omega }}{\int }{\mathrm{\psi }}_l{\mathrm{\psi }}_m\mathrm{d}x（23）

\begin{array}{c}{J}_{lm}=\sum _{e_K\in \mathrm{\Gamma }}\underset{e_K}{\int }\alpha \left[{\psi }_l\right]\left[{\psi }_m\right]\mathrm{d}s\end{array}（24）

N_l=\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }{g}_{\mathrm{2}}{\psi }_l\mathrm{d}s-\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }\frac{g_{\mathrm{1}}\partial {\psi }_l}{\partial n}\mathrm{d}s（25）

F_l=-\sum _{e_K\in {\mathrm{\Gamma }}_D}\underset{e_K}{\int }\alpha g_{\mathrm{1}}{\psi }_l\mathrm{d}s-\frac{\mathrm{1}}{D}\underset{\mathrm{\Omega }}{\int }{P}_{\mathrm{0}}{\psi }_l\mathrm{d}x（26）

则求解变分问题（16）~（20）等价于求解代数方程组：

\left\{\begin{array}{l}{\mathit{A}}^T\mathit{X}+\mathit{M}\mathit{Y}=\mathit{N},\\ \mathit{A}\mathit{Y}-\mathit{J}\mathit{X}-\tilde{\mathit{M}}\mathit{X}=\mathit{F}\end{array}\right.（27）

由于质量矩阵\mathit{M}一定为一个可逆的分块对角矩阵，方程组（27）等价于求解：

\left(\mathit{A}{\mathit{M}}^{-\mathrm{1}}{\mathit{A}}^T+\mathit{J}+\tilde{\mathit{M}}\right)\mathit{X}=\mathit{A}{\mathit{M}}^{-\mathrm{1}}\mathit{N}-\mathit{F}（28）

上述非齐次线性方程组可以用列主元消去法直接求解。

2 数值算例

取双线性矩形单元对方程（12）进行离散求解。在单元区域\mathrm{\Omega }=\left[\mathrm{0,1}\right]\times \left[\mathrm{0,1}\right]上考虑方程（11），并取D=k=\mathrm{1}，真解W=\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }y\right)。则

\left\{\begin{array}{l}{P}_{\mathrm{0}}=\mathrm{4}{\mathrm{\pi }}^{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }y\right)+\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }y\right),\\ g_{\mathrm{1}}=\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }y\right),\\ g_{\mathrm{2}}=\left(\mathrm{\pi }\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\pi }x\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }y\right),\mathrm{\pi }\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\pi }x\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\pi }y\right)\right)\cdot \mathit{n}\end{array}\right.（29）

其中，\mathit{n}为区域\mathrm{\Omega }边界上的单位外法向导数。

定义W和有限元解W_h间的误差为e_h=W-W_h。记e_h的L^{\mathrm{2}}范数为||e_h|{|}_{L^{\mathrm{2}}}， e_h的H^{\mathrm{1}}半范数为|e_h{|}_{H^{\mathrm{1}}}：

{‖e_h‖}_{L^{\mathrm{2}}}={\displaystyle \sum _{i=\mathrm{1}}^{N_h}}\underset{K_i}{\int }{\left(W-W_h\right)}^{\mathrm{2}}\mathrm{d}x（30）

{\left|e_h\right|}_{H^{\mathrm{1}}}={\displaystyle \sum _{i=\mathrm{1}}^{N_h}}\underset{K_i}{\int }{\left(\nabla W-\nabla W_h\right)}^{\mathrm{2}}\mathrm{d}x（31）

将区域\mathrm{\Omega }作均匀矩形网格剖分，如图1所示，网格尺寸分别取h=\frac{\mathrm{1}}{\mathrm{2}}\mathrm{、}\frac{\mathrm{1}}{\mathrm{4}}\mathrm{、}\frac{\mathrm{1}}{\mathrm{8}}\mathrm{、}\frac{\mathrm{1}}{\mathrm{16}}\mathrm{、}\frac{\mathrm{1}}{\mathrm{32}}。

构造双线性有限元空间并统一取稳定参数{\alpha }_K=\mathrm{1}。Q1元下相应误差的L^{\mathrm{2}}范数和H^{\mathrm{1}}半范数以及收敛阶如下表所示，可以看到，混合间断有限元格式在Q1元下误差的L^{\mathrm{2}}范数有2阶收敛速度，H^{\mathrm{1}}半范数有1阶收敛速度。

3 应用

藏东南地区冰川主要分布在海拔4000 m以上，其中面积大于100 \mathrm{k}{\mathrm{m}}^{\mathrm{2}}的冰川有雅弄冰川和恰青冰川，分别位于岗日噶布山和念青唐古拉山。藏东南冰川绝大多数是海洋性冰川，显著特征为高积累、高消融，对气候变化敏感。近年来，随着气候变暖，藏东南海洋性冰川出现强烈的面积萎缩和冰量损失^［27］。本节将模拟恰青冰川持续消融直至消失，即区域地壳上覆荷载消失的情况下地壳的回弹情况。恰青冰川的地形和冰川分布如图2所示。

要计算（2）式中的抗弯刚度，地壳厚度、泊松比以及杨氏模量是3个必需的计算要素。研究表明，在青藏高原地区，地壳泊松比约为0.25且不同区域变化不大^［28］。对于杨氏模量值的选取，已有文献表明\mathrm{120}\mathrm{G}\mathrm{P}\mathrm{a}是地面附近岩石的材料特性，而深处岩石的杨氏模量可能因高压而变得更大^［29］。因此，本文取泊松比为0.25，模拟杨氏模量的取值范围为\mathrm{90}~\mathrm{140}\mathrm{G}\mathrm{P}\mathrm{a}，地壳厚度为\mathrm{10}~\mathrm{35}\mathrm{k}\mathrm{m}。

取计算区域为\mathrm{\Omega }=\left[\mathrm{0,31}\mathrm{200}\right]\mathrm{m}\times \left[\mathrm{0,31}\mathrm{200}\right]\mathrm{m}的正方形区域，将此计算区域均匀划分成规格为\mathrm{600}\mathrm{m}\times \mathrm{600} m的小正方形，并取零边界条件：

\left\{\begin{array}{l}W=\mathrm{0},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega },\\ \frac{\partial W}{\partial n}=\mathrm{0},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\mathrm{。}\end{array}\right.

在以上假设下，由式（2）可知抗弯刚度：

D=\frac{E{h}^{\mathrm{3}}}{\mathrm{11.25}}（32）

取地幔密度{\rho }_m=\mathrm{3200}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，地壳密度{\rho }_c=\mathrm{2800}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}，重力加速度g=\mathrm{9.8}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}。由（4）式计算可得：

k=\mathrm{3920}（33）

则由（5）式可得方程：

{∆}^{\mathrm{2}}W+\frac{\mathrm{3920}\times \mathrm{11.25}}{E{h}^{\mathrm{3}}}W=\frac{\mathrm{11.25}}{E{h}^{\mathrm{3}}}{P}_{\mathrm{0}}（34）

其中，P_{\mathrm{0}}为冰川所造成的荷载。

首先，取定泊松比\nu =\mathrm{0.25}，杨氏模量E=\mathrm{1}\times {\mathrm{10}}^{\mathrm{11}}\mathrm{P}\mathrm{a}。对地壳厚度\mathrm{10}~\mathrm{35}\mathrm{k}\mathrm{m}时冰川完全消融时地壳的回弹形变进行模拟，结果如表2所示。

然后，取泊松比\nu=0.25，地壳厚度h=\mathrm{10}\mathrm{k}\mathrm{m}。对杨氏模量为90~140 GPa时，冰川完全消融时地壳的回弹形变情况进行模拟，结果如表3所示。

由表2和表3的结果可知，随着地壳厚度增大及杨氏模量增加，地壳的形变逐渐缩小。本文尝试用幂律函数f\left(x\right)=k{x}^n对表2和表3中的数据进行拟合，结果如图3和图4所示。

由图3和图4可以看到，地壳厚度-最大形变的幂律拟合结果为f_{\mathrm{1}}\left(x\right)=\mathrm{96.71}{x}^{-\mathrm{3.00}}，杨氏模量-最大形变的幂律拟合结果为f_{\mathrm{2}}\left(x\right)=\mathrm{9.67}{x}^{-\mathrm{1.00}}。这和（2）式中抗弯刚度关于地壳厚度和杨氏模量的定义次数吻合。总的来看，在表2和表3的数据下，对模型W=a{h}^b{E}^c进行拟合，其中a\mathrm{、}b\mathrm{、}c均为常数，拟合结果为W=\mathrm{9}\mathrm{658.03}{h}^{-\mathrm{2.99}}{E}^{-\mathrm{0.99}}。进一步，利用\mathrm{S}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{l}方法^［30］对模型参数h\mathrm{、}E进行敏感性分析，并令\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}代表\mathrm{95}\mathrm{\%}置信水平下的置信区间。总体敏感性系数（\mathrm{S}\mathrm{T}）、一阶敏感性系数（\mathrm{S}\mathrm{1}）、二阶敏感性系数（\mathrm{S}\mathrm{2}）结果参见表4，其中“-”表示值未测得。

由敏感性分析可知，地壳厚度对地壳回弹的影响较大，杨氏模量对地壳回弹的影响较小，且地壳厚度和杨氏模量之间的交互影响不大。具体而言，取杨氏模量E=\mathrm{1}\times {\mathrm{10}}^{\mathrm{11}}\mathrm{P}\mathrm{a}，泊松比\nu =\mathrm{0.25}，地壳厚度为\mathrm{10}\mathrm{k}\mathrm{m}。计算可得区域最大形变为9.68×10^-2\mathrm{m}，形变分布如下图所示。

由图5可知，当恰青冰川完全消融时，计算区域中心处地壳回弹形变较大，并由中心向四周形变逐渐减小。图6中的地壳回弹形变分布图更加清晰地显示了这一趋势。

4 结论

本文利用混合间断有限元方法，对重调和型方程进行数值求解，给出了数值算例，并验证了数值格式的收敛阶。然后本文模拟了我国最大的海洋性冰川——恰青冰川完全消融时冰川所覆区域地壳的垂直形变情况。计算结果表明，冰川消融时，地壳厚度对回弹形变的影响较大，而地壳岩石的杨氏模量对回弹形变的影响则相对较小。本文的结果为青藏高原地区的地灾预测和防控提供了一定的参考价值。

参考文献

[1] Fu Y N， Freymueller J T.Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurements ［J］.J Geophys Res Solid Earth， 2012， 117： B03407.

[2] Peltier W R， Argus D F， Drummond R.Space geodesy constrains ice age terminal deglaciation： The global ICE‐6G_C （VM5a） model ［J］.J Geophys Res Solid Earth， 2015， 120（1）： 450-487.

[3] Watson C， Tregoning P， Coleman R.Impact of solid Earth tide models on GPS coordinate and tropospheric time series ［J］.Geophys Res Lett， 2006， 33（8）： L08306.

[4] Qiu J.China： The third pole ［J］.Nature， 2008， 454： 393-396.

[5] Yao T， Thompson L G， Mosbrugger V， et al.Third Pole Environment （TPE） ［J］.Environ Dev， 2012， 3： 52-64.

[6] Rao W L， Liu B， Tang H， et al.Progress in studies on crustal uplift of the Qinghai-Xizang Plateau based on GNSS and GRACE ［J］.Reviews of Geophysics and Planetary Physics， 2025， 56（1）： 26-44.

[7] 饶维龙， 刘斌， 唐河， 等.青藏高原地壳隆升的GNSS与GRACE联合研究进展［J］.地球与行星物理论评，2025， 56（1）： 26-44.

[8] Wang L.Glacier Variational in the Yarlung Zangbo River Basin in China from 2000 to 2016 ［D］.Xi’an： Northwest University， 2021.

[9] 王璐.2000—2016 年中国境内雅鲁藏布江流域冰川变化研究［D］ 西安： 西北大学， 2021.

[10] Bungum H， Olesen O， Pascal C， et al.To what extent is the present seismicity of Norway driven by post-glacial rebound？ ［J］.J Geol Soc London， 2010， 167（2）： 373-384.

[11] Fjeldskaar W， Lindholm C， Dehls J F， et al.Postglacial uplift， neotectonics and seismicity in Fennoscandia ［J］.Quaternary Sci Rev， 2000， 19（14-15）： 1413-1422.

[12] Cossart E， Mercier D， Decaulne A， et al.Impacts of post‐glacial rebound on landslide spatial distribution at a regional scale in northern Iceland （Skagafjörður） ［J］.Earth Surf Proc Land， 2014， 39（3）： 336-350.

[13] 许厚泽 ．固体地球潮汐［Ｍ］．武汉： 湖北科学技术出版社， 2010.

[14] Shu Y H， Shi X H， Chen H L， et al.Isostasy and its application in tectonic geomorphology research ［J］.Geological Review， 2022， 68（4）： 1171-1190.

[15] 舒远海， 石许华， 陈汉林， 等.地壳均衡理论及其在构造地貌研究中的应用［J］.地质评论， 2022， 68（4）： 1171-1190.

[16] Li Y H.On the mixed generalized difference method for biharmonic equations ［J］.Acta Scientiarum Naturalium Universitatis Jilinensis， 1993， 31（3）： 19-30.

[17] 李永海.解双调和方程的一种混合广义差分法［J］.吉林大学自然科学学报， 1993， 31（3）： 19-30.

[18] Chen G， Li Z， Lin P.A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow ［J］.Adv Comput Math， 2008， 29（2）： 113-133.

[19] Wang T.A mixed finite volume element method based on rectangular mesh for biharmonic equations ［J］.J Comput Appl Math， 2004， 172（1）： 117-130.

[20] Liu J K， Lai X J.Semi-analytic solution and error estimation for the biharmonic equations ［J］.Journal of Tianjin University， 1993， 1： 92-99.

[21] 刘嘉焜， 赖学坚.重调和方程的半解析解法及误差估计［J］ .天津大学学报， 1993， 1： 92-99.

[22] Nunn J A， Sleep N H.Thermal contraction and flexure of intracratonal basins： A three-dimensional study of the Michigan basin ［J］.Geophys J Int， 1984， 76（3）： 587-635.

[23] Wang C， Wang J.An efficient numerical scheme for the biharmonic equation by weak Galerkin finite element methods on polygonal or polyhedral meshes ［J］.Comput Math Appl， 2014， 68（12）： 2314-2330.

[24] Guo H， Zhang Z， Zou Q.A C 0 - linear finite element method for biharmonic problems ［J］.J Sci Comput， 2018， 74（3）： 1397-1422.

[25] Argyris J H， Fried I， Scharpf D W.The TUBA family of plate elements for the matrix displacement method ［J］.Aeronautical J， 1968， 72（692）： 701-709.

[26] Bogner F K， Fox R L， Schmit L A.The generation of interelement compatible stiffness and mass matrices by the use of interpolation formulas ［C］// Proceedings of the Conference on Matrix Methods in Structural Mechanics.Washington： OTS， 1965： 397-443.

[27] Andreev A B， Lazarov R D， Racheva M R.Postprocessing and higher order convergence of the mixed finite element approximations of biharmonic eigenvalue problems ［J］.J Comput Appl Math， 2005， 182（2）： 333-349.

[28] Gudi T， Nataraj N， Pani A K.Mixed discontinuous Galerkin finite element method for the biharmonic equation ［J］.J Sci Comput， 2008， 37（2）： 139-161.

[29] Feng J， Wang S， Bi H， et al.An hp-mixed discontinuous Galerkin method for the biharmonic eigenvalue problem ［J］.Appl Math Comput， 2023， 450： 127969.

[30] Meinesz F A V.Une nouvelle methode pour la reduction isostatique regionale de l’intensite de la pesanteur ［J］.Bulletin Geodesique， 1931， 29（1）： 33-51.

[31] Nadai A.Theory of flow and fracture of solids （Ⅰ） ［M］.New York ： McGraw-Hill， 1950.

[32] Zhao J B.Study on characteristic of glacier change and movement in Southeast Tibet ［D］.Huainan： Anhui University of Science and Technology， 2024.

[33] 赵晋彪.藏东南地区冰川变化与运动特征研究［D］.淮南：安徽理工大学， 2024.

[34] Li C J， Wang Y， Liu L J， et al.Lithospheric deformation and corresponding deep geodynamic process of the SE Tibetan Plateau ［J］.Science China Earth Sciences， 2025， 55（5）： 1351-1376.

[35] 李长军， 王洋， 刘丽军， 等.青藏高原东南缘岩石圈变形特征及其深部动力学过程［J］.中国科学： 地球科学， 2025， 55（5）： 1351-1376.

[36] Zheng W H.Computational simulation on geodynamic problems ［D］.Wuhan： Huazhong University of Science and Technology， 2009.

[37] 郑文衡.地球动力学若干问题的计算仿真研究 ［D］.武汉：华中科技大学， 2009.

[38] Sobol I M.Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates ［J］.Math Comput Simulat， 2001， 55（1/2/3）： 271-280.

基金资助

国家自然科学基金(42471091)

西藏自治区科技重大专项(XZ202201ZD0003G)

中国科学院成都山地灾害与环境研究所自主部署项目(IMHE-ZDRW-02)

数值仿真四川省教育厅高校重点实验室开放研究项目(2025SZFZ004)