中国大陆是多块体拼合的产物，显生宙以来，分别受古亚洲洋、特提斯洋和环太平洋三大构造域所控制，形成横跨东西和贯穿南北的不同造山带。中国东部和南部分别受到太平洋板块和印度板块俯冲—挤压作用，所处的特殊地理位置和构造背景赋予造山带不同的构造特点，是揭示和研究全球造山带动力机制和演化史的理想场所之一^[1~3]。由于地表和岩石圈作用过程的特殊组合，不同造山带和造山高原的演化常存在不同的机制^[4~6]，同时野外获取和观察到的地质数据和现象，可能仅反映了其演化进程中的某一个阶段^[7]。随着科学技术的提高，物理模拟和计算机模拟技术逐渐为地球动力学的成因量化提供了解释的可能性，但由于物理模拟实验存在时空尺度的局限，还难以有效还原地质演化过程。计算机模拟具有不受时空限制的优势，各种参数包括受力大小、规模、材料属性等可直接来源于研究对象的测量和实验，有利于保存研究对象模型的相关信息^[8~9]，正是因为这些特点，它已成为当前国内外解决地球动力问题的强有力工具。尽管如此，数值模拟在应用于解决中国造山带复杂的动力演化过程中仍存在一些亟待解决的关键问题，如多期造山的边界条件设置，复杂造山带的多场耦合处理，高精度模拟过程中的模拟效率与可靠性等问题。造山带动力演化是一个极其复杂的过程，不同研究人员从不同侧重点来研究中国造山带动力机制，很难判断模型的优劣。因此，文章在梳理中国造山带的数值模拟技术方法基础上，针对典型的动力学模拟成果和存在的应用问题，进一步探讨数值模拟在造山带动力学模拟的研究及发展趋势，以期为推动相关领域的深入研究提供参考和借鉴。

1 造山带动力学数值模拟的主要方法

造山带动力学数值模拟基于连续介质方法，须遵循能量、质量、动能等定律方程^[7]，通过计算机数值计算方法，如有限元、有限差分、边界元等(表 1)，对模型域内设定的控制方程进行求解，从而得到求解结果。而造山带动力数值模拟中，数值计算方法常见的是有限元(FEM)和有限差分(FDM)，两种计算方法将求解域分成有限的离散单元体进行近似求解，其差异在于求解的过程。有限元为隐性求解的过程，有利于解决线性问题；而有限差分为显式求解，有利于解决非线性和大形变的问题，研究人员根据所解决的实际问题选择的不同的计算方法^[10]。

1.1 基本原理

造山带动力模拟遵循着自然科学法则，即不同的物理介质遵循不同的能量、质量和动能方程，如固体中的以X_i表示空间位置的平衡方程表示为：

$ \frac{{\partial {\delta _{ij}}}}{{\partial {X_j}}} + pg = 0 $

(1)

其中δ_ij为应力张量，pg为体力向量的分量；在岩石圈低速流动特征的平衡方程为斯托克斯公式(Stokes)：

$ \frac{{\partial {\tau _{ij}}}}{{\partial {X_i}}} - \frac{{\partial P}}{{\partial {X_i}}} + pg = 0 $

(2)

τ_ij为剪切力；X_i为空间位置；P为流体压力；pg为体力向量的分量。物理介质对外界或内部作用的响应遵循着能量、质量、和动能守恒定律，而物理介质的形变是通过本构方程来确定。在造山带动力模拟中，本构关系常见的有弹性体、弹塑性体、黏性体本构关系。其中弹性体的本构方程为：

$ {\delta _{ij}} = \omega \theta + 2\mu {\varepsilon _{ij}} $

(3)

δ_ij为应力张量，ε_ij为应变空间张量，ω和μ为拉梅系数，θ为体积模量；弹塑体的本构方程为

$ {\delta _{ij}} = {D_{ijkl}} \cdot {\varepsilon _{kl}} - d\omega \frac{{\partial F}}{{\partial {\varepsilon _{ij}}}} $

(4)

D_ijkl为弹性系数，ε_kl为应变张量，dω为条件确定的尺度因子，F为屈服函数，ε_ij为应变空间变量。黏性体的本构方程为：

$ \dot \varepsilon = {\rm{B}}\delta _{ij}^n $

(5)

$\dot \varepsilon $为应变率，B为常数，δ_ij为应力张量，n为应力指数；在深部岩石圈形变以韧性变形为主，可以引入温度和活化能参数方程，即幂律流变体方程来描述：

$ \dot \varepsilon = {\rm{A}}\delta _{ij}^n\exp \left( { - \frac{Q}{{{\rm{RT}}}}} \right) $

(6)

A为比例常数，δ_ij为应力张量，n为应力指数，Q为活化能，R为气体常数，T为绝对温度；黏弹体方程为：

$ {\sigma _{ij}} = E{\varepsilon _{ij}} + \eta \dot \varepsilon $

(7)

其中E为弹性模量，η为黏滞系数，$\dot \varepsilon $为应变率。不同本构关系存在不同的特点(表 2)。

1.2 几何非线性的计算

造山带动力模拟涉及到复杂的介质划分和大形变等非线性问题，尤其是对大形变中的网格形变和接触非线性的求解，而在模拟中求解方法有拉格朗日法(包括完全拉格朗日和修正的拉格朗日法)和拉格朗日-欧拉法(Arbitrary Lagrangian-Eulerian，ALE)。拉格朗日法以t_o初始构型作为参考构型，网格和物质节点随形变而运动，物质节点的位移矢量为X=(X₁, X₂, X₃)，X_i为物质节点的空间坐标，其数学表达式为x=x(X, t)^[19]；而接触面两侧的节点要求是保持重合一致的，但拉格朗日方法在接触的求解区域是未知量，在大位移或大形变中极容易发生网格和节点不一致或形变过大导致网格畸变，使收敛性差甚至不收敛^[20]。为了解决因大形变而导致网格畸变的问题，研究人员开发了ALE法，该种方法结合了拉格朗日和欧拉方法的特征，引进了相对固定不变的参考构型，参考构型由位移矢量ξ确定，其表达式为ξ=ξ(X, t), ξ在空间的运动规律为x=x(ξ，t)，最大的特点是网格可以独立于物质节点而移动，能将网格形变控制在最小的幅度^[19]。现有的商业模拟软件如Comsol、LS-DYNA等，能提供ALE算法解决接触、自由表面等非线性形变问题。

在岩土工程领域中，研究人员也开发了耦合的欧拉-拉格朗日方法(Coupled Eulerian-Lagrangian, CEL)，该法结合了欧拉和拉格朗日方法的优点，将形变区域自动划分为欧拉区域和拉格朗日区域进行有限元计算，采用显式动力差分方法求解。与ALE方法相比，CEL的网格是固定不动的，物质节点可在网格中自由形变，从而进一步降低网格畸变的概率^[21]。在不同介质界面的求解中(如流体和固体)，CEL对欧拉域(流体)与拉格朗日域(固体)之间的接触界面进行离散化，通过罚接触法的通用接触算法，计算作用于拉格朗日域和欧拉域的节点的接触力即罚力F_P：

$ {F_P} = {K_P} \times {d_P} $

(8)

K_P为罚刚度, d_P为罚位移；罚接触法允许欧拉域的节点进入拉格朗日域，避免了在接触边界形变过大的问题^[22]。由于CEL能够有效地解决大变形的非线性问题，在岩土工程领域得到了广泛地应用^[21~22]。

1.3 边界条件的处理

边界条件指约束模拟求解结果的条件，是判断模型是否具有真实性的重要参数之一。造山带动力模拟中，对大尺度复杂的边界条件处理，如各单元复杂的物理参数、地质构造等，研究人员一般对其进行简化和扩大模型规模或固定未知边界条件。随着模拟精度的提高，实际观测的资料为模拟提供了重要的边界条件，如GPS观测的地表位移和地震波速资料等^{[12, 23~24]}。

同时，地质演化过程中受剥蚀作用影响，为还原其初始边界条件，要对现今剥蚀地形进行恢复处理。研究人员在模型浅表引入剥蚀修正边界条件，如Kaus等^[25]在探讨新生代台湾山脉隆升机制中，采用了剥蚀与地形局部坡度成正比的扩散模型：

$ \frac{{\partial h}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{{\rm{K}}_{\rm{e}}} \times \frac{{\partial h}}{{\partial x}}} \right) $

(9)

h代表地表高程，x代表水平位移量，t表示时间尺度，K_e表示剥蚀常数。模拟结果认为新生代菲律宾板块在俯冲过程，弱岩石圈和强剥蚀率使台湾山脉的隆升具有不稳定性。尽管剥蚀与地形隆升的表达式还没能达成一致的认识，但大部分研究人员采用地形隆升与剥蚀率成正比的模型，如傅容珊等^[26]在研究青藏高原的隆升模型。

另外，当造山带剥蚀作用较小时，研究人员通过将地形起伏的模拟结果与实际地形进行修正，如Lu等^[27]在模拟武陵—雪峰造山带挤压形变机制中，通过计算了模拟与实际地形之间的关系数以及差值方差来调整造山带的隆升率，表达式如下：

$ {\rm{corrEff}} = \frac{{\sum\nolimits_{i = 1}^N {\left( {{Z_1}\left( i \right) - {{\bar Z}_1}} \right) * \left( {{Z_2}\left( i \right) - {{\bar Z}_2}} \right)} }}{{\sqrt {\sum\nolimits_{i = 1}^N {{{\left( {{Z_1}\left( i \right) - {{\bar Z}_1}} \right)}^2}} * \sum\limits_{i = 1}^N {{{\left( {{Z_2}\left( i \right) - {{\bar Z}_2}} \right)}^2}} } }} $

(10)

$ {\rm{variance}} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left( {{Z_{12}}\left( i \right) - {{\bar Z}_{12}}} \right)}^2}} $

(11)

Z₁表示模拟的地形，Z₂表示实际观测地形，Z₁和Z₂表示它们对应地形的平均值，Z₁₂和Z₁₂分代表模拟地形与观测地形之间的差值和平均差值，corrEff表示两种地形的敏感性，而variance值表示对两种地形高程差异的敏感性，当corrEff值越大同时variance值越小，说明模型与观测地形的吻合较好。修正后的地形与实际地势有较好的吻合度，通过模拟结果与实际对比，他们认为武陵—雪峰造山带下地幔岩石圈的黏度可能低于克拉通岩心下地幔岩石圈黏度。

1.4 网格划分

网格划分就是按一定的规律(三角网、四面体、六面体等)，将模型对象离散为待求解的单元体。网格的质量对模拟的精度和求解效率有直接的影响，过大的形变可能会引起网格扭曲畸变，从而影响模拟的求解效率和精度。在模拟中，网格重构是对网格质量进行调整的常见方法，即对小形变后的网格形变X进行重新划分并插入应力变量^[28]，其表达式为：

$ X = F\left( {{X_O},P,W} \right) $

(12)

X_O为形变结果，P和W分别为控制单元位置和插值函数；通过形变后的网格进行重划和新插值函数，在很大程度保持了计算网格的质量。同时，当前可用的ALE算法能降低网格畸变，但在不同介质的界面形变中，网格重划和ALE法也很难保证网格畸变的问题^[29]；为了解决大形变中网格畸变问题，研究人员把FDE和FED与MIC(Marker-in-cell，网格—粒子)结合，而MIC是一种基于欧拉方法描述物质节点形变的方法，通过FDE或FED线性插值求解网格内部物质节点的粒子未知量(速度、时间步长、位置)，来描述大形变中网格质点运动^{[12, 29~30]}, 其表达式为：

$ \begin{array}{*{20}{c}} {{B_m} = \left( {1 - \frac{{\Delta {x_m}}}{{{d_x}}}} \right)\left( {1 - \frac{{\Delta {z_m}}}{{{d_z}}}} \right){B_i} + }\\ {\frac{{\Delta {x_m}}}{{{d_x}}}\left( {1 - \frac{{\Delta {z_m}}}{{{d_z}}}} \right){B_j} + \frac{{\Delta {x_m}}}{{{d_x}}}\frac{{\Delta {z_m}}}{{{d_z}}}{B_{kL}} + }\\ {\left( {1 - \frac{{\Delta {x_m}}}{{{d_x}}}} \right)\frac{{\Delta {z_m}}}{{{d_z}}}B} \end{array} $

(13)

B_m为粒子参数；B_P(P=i; j; k; L)为描述粒子所在的四个节点，d_x和d_z为单元体在水平和垂向的位移描述，Δz_m和Δx_m为在垂向和水平上，粒子到左下角处第i个节点的所需的距离。

另外，无网格法在地学领域研究逐渐兴起，在运用中只需获得介质节点信息和边界条件, 在获取节点下通过节点与权函数近似，从而消除对网格单元质量的限制，避免了对网格的依赖性^[31]。

1.5 物理场的耦合处理

造山带动力模拟中所涉及的基本物理场分析，主要包括固体静力场(应力、应变和位移等)、流体场、温度场等，以及它们之间的相互作用(耦合关系)在模拟求解中，如何考虑和处理不同介质物理场的耦合，是一个复杂的问题。就方法而言，在地球动力学模拟中，对于不同介质的耦合求解(如流体和固体)，研究人员多通过建立由不同方程控制的模型域(如固体域和流体域)，将一个域的求解结果通过介质接触面作为另一个域输入的条件，或通过不同的控制方程进行联合求解(强耦合求解)，如引入热力方程(公式5~6)，从而达到得到共同求解的结果^[32]。为了提高求解的收敛性，有的研究人员在求解中引入并行计算方法，如Morra和RegenauerLieb^[33]使用边界元和有限元相结合的方法，对流—固接触进行求解；Liu和Currie^[34]把上覆岩石圈弹性形变与软流圈黏性形变相结合，通过ALE法进行求解。

2 中国造山带动力学数值模拟的应用现状

造山带按动力机制可划分为洋—陆俯冲造山带和陆—陆碰撞造山带。典型的洋—陆俯冲造山带有太平洋板块东缘的科迪勒拉造山带和安第斯造山带，陆—陆碰撞造山带有喜马拉雅造山带、美国Laramide造山带、南非纳米比亚Damara造山带、澳大利亚中部Alice Spring造山带等^[3]。

Tapponnier和Molnar^[35]通过数值模拟研究了印度—亚欧板块碰撞下东亚大陆的挤出模式，相对于传统地质的研究方法，其模拟结果更具直观性，此后研究人员开始运用数值模拟方法研究地学领域的重大科学问题。Willett等^[36]提出了“双向挤压造山”(doubly vergent compressional orogens)动力模型，其模型描述了刚性块体向相对静止的楔状块体俯冲的动力过程，这为洋—陆俯冲、陆—陆碰撞模拟研究提供了重要的动力参考模型。其他研究人员在此基础上开发了不同的力学模型，如安第斯造山带是Nazca板块向南美大陆俯冲所形成的的典型双向挤压造山带^[37]；Sobolev等^[38]强调了南美板块向西漂移对新生代安第斯造山带的动力作用；Martinod等^[39]认为新生代安第斯造山带的隆升与平俯冲板块的几何形态改变有关；而Faccenna等^[37]的模拟结果突出了俯冲板块与上覆板块之间的剪切对安第斯造山带动力的作用。尽管研究人员从不同的侧重点对模拟结果进行阐述，但数值模拟的直观结果无疑带动研究人员们的思考。

数值模拟技术在地球动力学的运用与研究中不断地被开发和完善，从数理模型可以划分为2D平面应力模型、2D热力模型、3D热力模型等(表 3)；从模拟的内容可以分为板块汇聚和碰撞模拟、岩石圈流变模拟、地幔对流模拟、地震机制模拟、陆内深俯冲模拟、构造应力场模拟、构造形变模拟等。针对中国造山带典型的动力机制问题，研究人员开展了众多的定量、定性数值模拟研究，并取得了很多重要成果和认识。

2.1 造山带的隆升模拟—青藏高原

青藏高原是中国进行动力数值模拟研究的天然场所，也是数值模拟在中国造山带动力学运用研究的起点。

2.1.1 板块汇聚和碰撞模型

Tapponnier和Molnar^[35]把大陆岩石圈视为刚塑性，模拟了印度—亚欧板块碰撞下东亚大陆的挤出模式，结果认为包括青藏高原的东亚大陆在作用力下呈近似滑移线形变，并认为隆升中的刚性印度地壳没有被吸收。而Zhao等^[47]认为印度地壳是可以刚性缩进的，因此在模拟中考虑了岩石圈物质纵向的差异，把印度地壳视为刚性体，西藏地壳为低黏性质体，模拟刚性印度地壳注入西藏地壳的动力过程，结果显示青藏高原的隆起可以是印度地壳注入西藏地壳而形成“双地壳”的结果。

England和Houseman^[48]应用薄粘性层模型探讨了板块碰撞下青藏地区及周边地区的变形过程与机制，论证了青藏高原新生代岩石圈变形以增厚作用为主的动力模式，模拟结果显示在印度板块作用下的高原，由南向北呈均匀隆升。Houseman和England^[49]进一步利用薄板粘性模型，在模型中考虑板块的上浮力条件，模拟限定了青藏高原在陆陆碰撞过程中的地壳增厚与侧向挤出变形量的比例为3:1至4:1，模拟结果认为青藏高原的隆升在时空上是均匀的过程。傅容珊等^[26]采用幂律流变薄层模型，模拟了印度板块向北恒速推进对高原隆升过程的影响，模型引入了剥蚀修正参数，通过模拟结果认为高原隆升因受到多种因素制约，在空间和时间上应该是不均匀和不稳定的过程。

以上的学者们在模拟中，把印度板块向北推进视为高原隆升的主要或唯一机制，主要是通过模拟大空间尺度的印度—亚洲大陆碰撞作用来探讨高原形变过程及其机制问题。

2.1.2 岩石圈流变机制模型

随着模拟研究的扩展，学者们开始关注岩石圈流变性质对高原的隆升作用机制。最初是由Royden等^[40]通过模拟认为青藏高原东缘的下地壳可能是由相对较弱的物质组成，陆陆碰撞挤压作用可引起下地壳软弱的物质发生流动，从而提出了下地壳流模型。孙玉军等^[41]认为横向地壳流变强度的差异对高原形变起重要作用，并建立了二维热力模型，在考虑岩石圈的热状态下，在非线性介质形变中引入ALE和MIC相结合的方法进行求解，模拟结果显示在高原东北部周缘的岩石圈的扩张和增厚与其岩石圈流变性质存在极大关联，地壳低黏、高速流变性质控制高原整体隆升速度，而地幔流对高原隆升作用不大。

2.1.3 重力均衡机制模型

罗焕炎等^[50]认为重力均衡在高原隆升中具有重要的作用，因此在模型中仅考虑重力作用和在印度板块的挤推对高原的影响，模型中也考虑断裂带的存在，模拟结果认为印度板的作用力是高原隆升的主要驱动力，但重力均衡调整所形成的水平压应力才是高原隆升的主要机制，同时重力均衡使在印度板块和欧亚板块过渡带处的莫霍面形成了40 km高差。Liu和Yang^[51]在模拟青藏高原的伸展塌陷机制中，考虑了岩石圈流变的横向和纵向变化以及相关的边界条件，使用GPS速度场为约束条件等，认为当高原低于其目前海拔的50%时，整个高原以走滑和逆断层作用占主导，未出现东西向地壳扩张，而当高原达到目前海拔的75%时便会发生显著的地壳伸展变形。

2.1.4 构造应力场模拟

构造应力场模拟是揭示地壳或岩石圈构造应力场成因与形成机制的重要手段。汪素云和陈培善^[52]通过二维有限元数值模拟计算了东亚地区应力场分布，并探讨了其成因机制。但其模拟未考虑东亚地形差异的影响。据此，张东宁和高龙生^[53]通过网格划分法区分了东亚的地表特征和横向不均匀性，其模拟结果显示东亚地区水平向最大主压应力轴走向的分布是以青藏高原为中心的辐射状格局，并且水平向最大主压应力值往东向逐渐减少。之后张东宁和许忠淮^[54]进一步利用三维粘弹性有限单元模型，模拟了青藏高原应力状态及构造运动。其模拟综合考虑了水平向和垂直向荷载力作用，结果显示在印度洋板块的水平挤压和重力势能的共同作用下，青藏高原物质存在向东挤出的水平运动，同时认为太平洋板块向西对亚欧板块的俯冲可能影响到青藏东北缘。

2.1.5 三维构造热力模型

Lechmann等^[42]采用平行三维有限元岩石圈和地幔演化模型(LaMEM)，研究影响印度—亚洲大陆碰撞模型参数，本构关系为随温度变化的黏性流变体，并考虑了岩层性质和印度下地壳强度与上地壳间的差异，活动断层通过GPS施加速度约束。结果显示由流变场和温度场控制的三维模型能有效地模拟真实的地表速度，重力势能的横向变化对高原压强变化具有重要作用，垂向的粘度和密度差对地表速度场影响小。同时认为仅靠高原地表速度场不足以使高原深部岩石圈发生解耦。Pusok等^[6]也采用了热力LaMEM模型研究印度—亚洲碰撞带的动力机制，模型考虑了地壳流和重力坍塌作用，模拟分析了海洋俯冲面(海沟退缩)和大陆碰撞带(海沟前进、板快脱离、地形抬升、侧向挤压)的动力学特征。模拟结果认为青藏高原的隆升不仅与来自印度板块的作用相关，地壳横向介质的差异(刚性的塔里木盆地)对隆升具有重要的作用。Chen等^[43]通过三维热力模型来研究从汇聚边缘到高原内部流变特征，为了模拟更真实的岩石圈层流变性质，利用了室内实验结果和在真实高原介质对比下，用辉绿岩代表高硬度的岩石圈层，用基性麻粒岩代表较为软弱的岩石圈层，俯冲板块的下地壳流变强度均使用斜长石来描述。模拟结果显示，在碰撞作用下，在较弱的亚洲地壳模型中，张力出现在汇聚边界的以北地区，并可形成高原；相反，在较强的亚洲地壳模型中，强硬地壳会抑制高原的形成。

从后来学者们所采用的力学模型和模拟方向可以看出，青藏高原的动力模拟研究不再局限单一物理场，物理参数贴近实际和量化过程更细化。

2.2 地幔对流和构造形变数值模拟—天山造山带

新生代天山的陆内再造山，普遍认为与印度—欧亚大陆碰撞所产生的远程效应有关，但对于天山“复活”机制存在不同看法，有的认为直接受印度—亚洲大陆碰撞影响，地壳受到挤压发生缩短从而使天山隆升^[49]；也有的认为在天山地区上存在地幔对流，天山隆升的主导因素是由于地幔对流上涌所导致，而不是挤压缩短致使天山地壳增厚^[55]。

2.2.1 构造形变模型

Neil和Houseman^[17]采用薄黏层描述陆壳形变特征，在忽视应力垂直变化下，模拟了在印度—亚洲大陆碰撞下塔里木盆地和天山的应变场，结果显示尽管在碰撞作用下塔里木盆地的形变不明显，天山地区也可以形成显著的地壳增厚。但因模型局限于平面分析，未能全面反映出青藏高原、塔里木盆地和天山岩石圈间形变能量传递关系，因此文章在采用Neil和Houseman^[17]模型的初始本构关系下，利用Ansys软件设计了经过青藏高原、塔里木盆地和天山在内的小尺度简易剖面模型，模型中把塔里木盆地视为刚性体，结果显示在挤压作用下形变能量从青藏高原越过刚性盆地到达天山，从而使天山隆升(图 1)。

Lei等^[23]在三维模型中考虑了区域活动断层的作用，利用GPS速度场作为约束条件，模拟了天山地壳的应变场。结果显示，天山地壳运动方向总体为北北东，区域构造应力场为近南北向压缩，所得的结果支持天山地壳形变是由印度—亚洲大陆碰撞所产生的远程效应的观点。但因该模型考虑的是短期动力效应，并没考虑岩石圈温度和韧性形变特征。因此，Lei等^[44]进一步改进了模型，引入岩石圈流变因素，并考虑了天山地区与南北向岩石圈物质差异和在重力的作用下引起的蠕变对地壳隆升的作用，建立了三维热力模型，结果发现帕米尔和塔里木地块的共同挤压作用是天山隆起的主要原因。卢双疆和何建坤^[56]鉴于此模型主要考虑了板块碰撞产生的远程效应力或天山域内岩石圈流变性质，而对区域物质差异对天山隆升影响的认识仍不足，因此通过在模拟中考虑块体间的地壳性质差异，得到了新的模拟结果。该结果虽与Neil和Houseman^[17]有很多相似性，但发现不同块体的地壳流变性质差异对天山隆起具有明显影响，印—藏汇聚的应变经青藏高原吸收后可以跃过刚性塔里木导致较软的天山地区发生强烈变形，进而提出天山地区新生代以来强烈地壳变形可能与青藏高原—塔里木—天山一带岩石圈流变学存在横向不均一性有关。

2.2.2 地幔对流模型

地幔对流涉及热流固耦合模拟，一般的方法是将地幔视为不可压缩的黏性流体，由纵向顶底部位存在温差所引起或由于板块俯冲引起所形成的地幔扰动(图 2)。

Liu等^[12]针对天山地区是否存在地幔对流和对天山形变的程度进行了模拟研究，模型中通过地震波波速约束材料属性，设置上下层间存在温差引发地幔扰动，其中瑞雷系数可判断是否形成对流，结果显示北天山—准噶尔盆地下方的上地幔中存在一个相对较强的逆时针对流, 但上地幔小尺度对流并不是天山复活隆升的主要原因，塔里木盆地向北挤压才是关键因素。

2.3 岩石圈拆沉数值模拟—昆仑造山带

自二叠纪晚期东昆仑进入陆内造山阶段，而三叠纪东昆仑岩浆活动却十分活跃，大量发育晚三叠花岗岩^[57]；刘成东^[58]认为东昆仑晚三叠强烈的岩浆活动与幔源底侵有关，同时底侵的幔源岩浆物的冷却又可能导致拆沉的发生。

詹华明等^[13]将拆沉作用与幔源底侵作用联系起来，通过数值模拟探讨了东昆仑三叠纪以来，幔源底侵的岩浆物冷却是否能导致拆沉的发生等问题，模型考虑了岩石圈层流变性质，同时模拟中运用库仑-莫尔强度线的破坏准则和垂向位移量来判断拆沉是否发生。模拟结果显示，东昆仑造山带在幔源底侵以后可能发生过拆沉作用，而拆沉量在地区间存在差异，同时认为底侵作用是导致拆沉作用的重要前提，底侵幔源的岩浆在冷凝过程中导致岩石密度增加，进而会引发东昆仑岩石圈出现重力不稳从而触发拆沉作用(图 3)，而拆沉作用是否会进一步触发东昆仑大规模岩浆活动，则与拆沉量有关。

2.4 地震机制模拟—龙门造山带

龙门造山带作为阻挡巴颜喀拉块体向东挤出的东界，发生了强烈的走滑—逆冲构造形变。数值模拟方法通过模拟构造应力场的分布与演化规律，是揭示地震机制的重要的工具。

杨辉等^[15]采用二维黏弹模型，考虑了块体间岩石圈流变差异和重力作用下，模拟结果认为重力均衡调整了龙门山深部静力结构，在印度板块碰撞挤压作用与青藏高原东缘重力势能差的共同作用下，高原东缘的深部物质在龙门山深处可潜入地幔，使龙门山上、中地壳层出现应力集中；而青藏高原和四川盆地间的介质属性、构造格局和流变属性的差异等对汶川大地震的孕育和发生也起到了重要作用；但该模型没有反映出局部断裂活动特征，也回避了地形起伏对应力场分布的影响。权凯^[59]选取龙门断裂中南段的剖面为模型对象，考虑地形起伏所引起的重力势能差以及岩性差异对应力分布的影响，模拟结果显示该区的现今构造应力场发育主要受到龙门山断裂带的活动特性、区域应力场分布状态、隆升变形、构造特征及区内岩石物理性质等因素控制，在前山断裂及后山断裂带的深部转折部位均有水平应力集中现象的出现。现根据文献[59]的模拟条件重新进行简单模拟后，得到与其基本一致的结果(图 4)，并显示龙门断裂总体位移向东运动，而深部地壳运动与浅层地壳位移可能存在差异，深部地壳位移不一定引起浅层地壳发生较大运动，只有当深部地壳位移量达到极限时，才会通过突发地震的方式释放应力。

考虑到在水平面上，虽然二维或平面模型能模拟高角度的走滑型地震，也可以在垂直剖面上模拟倾滑地震的发生, 但却无法模拟以走滑与倾滑分量兼有的真实三维地震位错^{[16, 60]}。柳畅等^[16]通过三维黏弹性模型研究了龙门山断裂带的应力积累和大地震复发周期，其中模型考虑了区域流变差异性，模拟结果认为龙门断裂两侧岩石圈层介质差异以及莫霍面深度在龙门断裂突变，是使龙门断裂出现应力集中的主要因素。上地壳深部的应力积累率大于浅部弹性地壳，区域压应力由南西向北东方向逐渐减小，而剪应力从南西向北东部不断的增大，并通过应力积累和地震应力降的计算认为8.0级大地震的复发周期为5400年。但模拟中没有考虑重力作用，而尹力和罗纲^[61]认为岩石圈在重力作用下和长期流变特征下趋于静岩压力(无偏张应力)，因而把模型的初始应力状态设置为静岩压力状态，模型中还考虑了区域黏滞大小、结构和断层几何形态对地表形变的影响，最后将模拟结果与短期地壳形变和长期构造形变相联系，认为震间加载的应力会在大地震中释放，震后岩石圈形成松弛效应，青藏东部在震后持续向东位移和抬升，在多个地震循环短期形变后，龙门地区整体抬升，此外认为断层的几何形状对龙门山形变具有重要的影响。

2.5 陆内深俯冲模拟—大别造山带

上世纪八十年代后期，大别造山带地区发现超高压变质矿物，尤其是含柯石英榴辉岩的发现，指示其形成深度被认为可能达到上百公里^[62]。

王飞等^[14]通过模拟结果认为超高压变质岩的折返至少经历了两个阶段，初始阶段洋壳俯冲把陆壳物质带到地幔深处，经过超高压变质作用后，在地幔和地壳物质之间密度所产生的浮力作用下折返到地壳浅层。范桃园等^[18]则提出了陆壳俯冲增生、角流俯冲和折返、浮力抬升折返和构造伸展抬升折返四阶段模型，各个阶段设定不同的边界条件，把增生楔及地幔物质以牛顿粘性流体来描述，通过数值热模拟得到了P-T-t轨迹，轨迹呈现顺时针演化，认为板块在俯冲的过程中，所携带陆源物质在初期形成增生楔，随后以角流为主，当洋壳物质拆离后，在浮力作用下使超高压变质岩快速地从深部折返到浅部，该模拟结果也认同浮力作用是超高压变质岩折返的重要因素。

后造山期延伸是每个威尔逊循环的最后阶段，但它们的动力学控制因素以及造山后延伸过程的作用仍然存在争议，而Dai等^[45]基于西大别造山带双缝合构造带(同向双俯冲带)，采用二维热力模型与MIC相结合的方法，以黏塑体描述岩石圈流变特征，研究了后造山期延伸的动力学机制，其中模型的顶部表面由12公里厚的“黏性空气”层组成，采用简化的总尺度侵蚀沉积法计算表面过程，最终模拟结果显示西大别造山带的后造山阶段可分为两个阶段，一阶段碰撞阶段包括超高压变质岩出露、俯冲块体断裂和残余岩石圈的折返过程，第二阶段代表了在受到残余岩石圈折返影响后，不同汇聚速率和地壳密度差异的再俯冲作用对后造山期延伸过程影响阶段。同时通过与西大别造山带实测资料的对比分析，认为西大别造山的后伸展动力学受华北地块再俯冲和高压/超高压碰撞后发生折返作用的控制。

2.6 弧形构造数值模拟—大巴造山带

大巴山位于秦岭造山带南部，平面上表现为一大尺度向南向突出的弧形构造带，关于大巴山前陆弧形构造动力演化机制存在不同观点。张国伟等^[63]认为大巴弧形构造与其前缘东西两端砥柱隆起有关，张岳桥等^[64]认为弧形构造的形成与造山前的伸展作用有关。武红岭等^[46]认为前陆砥柱在弧形构造形成机制中起到重要作用，因此模型把大巴山东西两侧的两个砥柱设置为刚性体，把大巴岩石圈描述为幂指数蠕变材料，模拟结果显示由于东西向两个砥柱的存在，在力作用下前陆形变朝着两砥柱中间较软地区发育，形成弧形构造(图 5)。但王瑞瑞等^[11]认为，仅靠两侧隆起阻挡模式无法完全解释大巴构造带构造特征，城口—房县弧形断裂在大巴弧形构造形成过程中扮演了重要作用。据此，其进一步建立了四个不同的参数模型进行了模拟对比分析，四个模型分别考虑了砥柱、弧形边界、砥柱叠加弧形边界以及综合砥柱叠加弧形边界与伸展断裂作用的共同作用，结果显示大巴山地区印支—燕山主造山期的最大水平主应力轨迹受到城口—房县弧形断裂和前缘两侧隆起的共同控制，进而提出大巴弧形构造应该是在早期伸展背景下形成的弧形断裂边界、前缘两侧隆起砥柱作用和底部滑脱作用为共同控制因素下形成的。总体来看，虽然众多学者在大巴弧形构造形成机制的量化模拟方面做了大量研究工作，也取得了许多重要认识，但已有的模型仍局限于二维平面上。

3 造山带动力模拟的趋势

(1) 模型由以二维向三维甚至“真三维”发展。在当前关于三维建模技术已经取得一定突破，如通过数字地质填图获取PRB(Point；Routing；Boundary)数据并投影于三维地质模型中^[65]。又如G-Plates开源软件与地理空间信息系统(GIS)相结合，通过向量几何和栅格数据投影至地球表面的构造板块上，从而形成反应真实地形等“真三维”模型^[66-67]；同时该软件设置了一个全球板块框架，模型中的边界约束条件可通过地震等资料来控制，模拟中岩层的厚度变化可以通过分布在变形网格中的点来跟踪；最近Dietmar Müller等^[67]就通过最新的G-Plates2.0软件研究了中—新生代澳大利亚和新西兰构造演化史。虽然由大量数据源支撑的“真三维模型”生成的浅表点线面精度较高，但目前只局限于浅表地壳的建模，而涉及深部构造投影或者通过浅部地表推测深部建模技术，需要进一步的试验和开发。

(2) 高精度(高分辨率)的模型使用。主要体现在模型的几何构造区分和介质辨识度细化；目前在国际上出现了一些高精度数值模拟技术，主要是与通过网格划分方法相关，如Koketsu等^[68]通过密集的监测系统构建了高精度的日本岛屿及其周边三维地壳结构模型，该模型的特点是通过优化网格划分提高了模拟精度，使得不同性质地层表面和界面的分辨率可达1 km左右。Ichimura等^[69]开发出无需计算机辅助或手动调整的三维高精度有限元模型的自动网格生成方法，该方法基于几何多重网格，变量预处理和多精度算法的迭代求解器加速计算来实现，并且被成功应用于大地震后的地壳形变模拟研究中。

(3) 模型边界条件的复杂化。随着数值分析方法完善和计算机运算能力的提高，造山带动力模型的边界条件从简化不断向复杂方面发展。如造山带隆升与剥蚀的修正和造山带热力条件的约束。在热力俯冲模型中包含了对造山带岩石圈的温度约束，但板块汇聚边界处的放射性热元素的富集是影响大陆俯冲—碰撞动力学过程和大陆碰撞区域热结构演化的重要因素^[70], 如宋菁菁等^[71]通过热力模型探讨了上、下地壳放射性生热率的敏感性对大陆俯冲的动力演化的影响。结果显示与下地壳相比，上地壳生热率的变化对大陆俯冲—碰撞动力学演化过程的影响较为显著，主要包括进入俯冲通道内的上地壳体积大小、碰撞区域内地壳熔融范围、俯冲下地壳物质折返的规模和两大陆的耦程度等四个方面。模型边界条件的复杂化往往意味着对造山带的描述更加详细、真实。

(4) 高级求解技术与传统FEM或FDM等方法相结合。造山带复杂、三维、高精度模型的使用，意味着对非线性的材料和接触问题求解的难度增强。为了提高计算效率，MIC、CLE等技术与FEM或FDM相结合，在未来造山带动力模拟中是一大趋势^{[12, 15, 29~30, 61]}；同时，利用超算硬件设施对地学领域里重大科学问题进行模拟，也是未来的趋势之一^[72]。

4 存在的问题

数值模拟技术在中国造山带动力学的运用，已经从静力结构物理场分析向考虑深部韧性形变方面等多场耦合发展。但从全局来看，在研究运用中仍存在一些问题。

(1) 复杂造山带仍然缺少整体有力的数值模型。尽管中国境内的一些造山带有很多的模拟研究，如青藏高原的动力机制等，但这些模拟研究都具有局部性。缺少全局的整体的动力模拟模型，更有些局部区域缺少相关的模拟研究，如位于中国华南大陆多期变形的数值模拟相对较少。最近Axen等^[73]通过二维热力模型模拟了晚白垩世Farallon板块平俯冲到北美大陆之下的动力学演化过程，结果显示平俯冲能将上覆板片岩石圈地幔较轻的物质刮掉并填充于软流楔，被刮下的岩石圈地幔物质抑制了岛弧岩浆作用。该模型为北美晚白垩Laramide造山运动提供了来自数值模拟的结果，而太平洋板块的俯冲如何作用于中国大陆的数值模拟仍需要进一步的研究。

(2) 数值模拟结果的可靠性。造山带动力学是一门多学科交叉的研究，动力模型主要是基于一种理想地球物理状态下的概念性的模型^[9]；在简化后的模型中，浅表的边界条件可通过大地测量技术(如GPS等)得到较好的验证；而涉及深部韧性形变、复杂或未知边界条件的处理时，还不能得到较全面的模拟效果和地质验证。通过约束等方法来解决未知边界效应，或者基于实验采用某种物理本构方程控制形变特征，但这容易造成夸大某一方面而忽视另一方条件的作用，从而影响模拟结果的可靠性。

(3) 计算能力满足不了复杂、高精度模型的需求。在深部韧性岩石圈模拟中，涉及复杂的非线性、多物理场耦合的问题，高级的求解方法与传统算法相结合虽能够有效地解决此类问题，但随着数据源的补充，如何在对多个非线性问题的求解中保证运算的稳定性和求解精度，使造山带数值模拟成为大陆动力学的研究难点之一。

5 结论

中国分布着不同构造类型的造山带，研究人员根据造山带特定的动力学问题进行了针对性模拟研究，如板块汇聚和碰撞、岩石圈流变机制、重力均衡、构造应力场、构造形变、拆沉、地震机制、地幔对流等。这些数值模拟有不同的控制方程、几何非线性计算、边界条件、网格划分和耦合等处理方法。数值模拟技术在中国造山带动力学的运用，已经从单一的静力结构物理场分析向力场、热场、流场等多物理场耦合方面发展，但在研究运用中受限于局部，缺少全局整体数值模拟，同时模拟结果的可靠性以及三维、高精度和非线性模型的运算效率等都亟待提高。这要求造山带地球物理化学等资料的不断完善以及更高级的数值计算方法。随着计算机的发展和学科之间的交叉合作，数值模拟运用于造山带动力机制的研究程度和理论认识必然会越来越深入。

致谢: 感谢匿名审稿专家的修改意见及编辑部的耐心编校修改。