川滇地区位于青藏高原东南缘,作为中国大陆强震活动最显著的区域之一,区域内复杂的地质构造背景、物质组成、流变学特征和独特的深部物理状态,使其成为研究大陆动力学的最佳实验窗口之一^[1-3]。随着中国大陆构造环境监测网络的建立和发展,川滇地区的全球导航卫星定位系统（global navigation satellite system,GNSS）站点密度也愈加丰富,为研究区域现今地壳变形模式提供了新的契机^[4-6]。从现有GNSS速度场数据结果来看,在川滇地区地表水平速度场存在着北强南弱、西强东弱的模式,整体上呈现为围绕喜马拉雅东构造结做顺时针旋转的特征^[7-9],其速度场观测结果所揭示的研究区域内现今水平运动及形变与其所处的板块运动环境及其内部地质构造等动力学环境紧密相关^[10-12],将地表的运动趋势和岩石圈深部的构造特征联系起来,加上岩石介质的约束,有助于人们深入研究地球内部动力学过程^[13]。由于川滇地区构造的复杂性,虽然目前研究区域已经拥有了较高密度的地表变形观测数据,但对该区域深部的动力学驱动机制在认识上依然存在争议。

关于川滇地区形变的动力学模式主要存在两类假说：第一类为块体挤出模型^[14-16],文献[17]提出了在欧亚板块碰撞下,分析了青藏高原侧向挤出的力学机制,为高原东南缘的川滇地区地壳运动研究提供了基础的理论框架,从大量的地质学观测中也逐渐认识到来自于印度板块强烈的汇聚作用主要被大型走滑断层所分解和吸收。GNSS数据结果也观测到在断裂带处的大规模滑动位移,同时地壳的挤压缩短导致青藏高原地壳由南向北逐渐增厚和隆升^[18-20],该理论很好地解释了现今研究地区观测到的走滑、旋转、拉伸和挤压的地壳变形过程以及沿块体边界的强震活动。但川滇地区内的部分次级块体则表现为不同于块体挤出模式的运动趋势,例如在菱形块体则表现为围绕着喜马拉雅东构造结做顺时针旋转的运动趋势,而在其西南方的红河断裂带及附近呈现出连续、弥散的变形特征。因此研究发现除了浅部地壳作用外,深部过程对研究区域的影响也不可忽略,便提出了第二类连续变形模型^[21-24],这种特征力主要来源于印度板块向北推挤,强烈的挤压作用对青藏高原冷的岩石圈地幔形成对流剥离和拆沉,青藏高原岩石圈由于负浮力而抬升,导致地壳层增厚,随后由于地壳增厚的不均匀性而形成重力势能差高原逐步朝外扩展；同时随着深部成像技术的发展^[25-26]（地震成像、电磁测深、深反射剖面等）,一些研究工作也强调了川滇地区的中下地壳具有较低粘度的物质或者存在熔融状态的低速体^[27]。在板块碰撞作用力的驱动下,川滇地区的较软的中下地壳部分物质以通道流的方式向东南方向流出^[28-30],遇到稳定的华南块体后导致了川滇地区地表的抬升和地壳增厚,从而引起形变^[31]。另外还有更深部的地幔对流模式也被提出来解释川滇地区地壳运动模式^[32],在印度板块和欧亚板块的碰撞下,印度岩石圈俯冲而带动了地幔软流圈层的运动,川滇地区的地幔软流层对地壳岩石所施加拖曳力的作用对地壳形变的影响也逐渐受到重视,但是以上几种作用在深部的动力学模式假说的驱动机制仍需更详细的数据资料和数值模拟实验来确认。

为了验证以上各类动力学模式假说,讨论不同力源对现今地壳变形影响,面对川滇地区复杂的三维岩石圈形变运动,传统的解析解和半解析解方法存在着一定的局限性,而随着计算机硬件技术的发展,尤其是高性能数值模拟计算平台的建立,近年来利用有限元数值模拟的方法来解决地壳形变问题一直是备受关注的重要科学热点。由于具备在解决高度非线性或者求解域几何形状比较复杂等问题方面的优越性^[29],相比于实验里进行的物理模拟,数值模拟分析方法具有便于多次实验、重复验证的特点,已成为解决各类地学问题的强大的手段之一^[33]。随着空间大地测量技术、地震层析成像技术、岩石力学物性实验分析等手段在川滇地区的深入研究获得了更成熟和丰富的地学数据,使得数值模拟方法的优势逐渐突显出来,即通过模拟实验具体分析川滇地区由于地质构造、地形载荷、重力加载及岩石流变性质等空间分布的不均匀性对该区域表层物质运动以及深部构造应力场格局的影响；厘清川滇地区构造演化的可能过程和驱动机制,深入认识川滇地区构造运动所反映的普遍性规律。

总得来说,关于川滇地区的构造变形特征和动力学机制的争议,可以贯穿于地质学、大地测量学、地震学、大地电磁学及数值模拟等大部分地学领域,厘清以上多种动力学机制对地壳形变的贡献对于研究地壳变形特征、地震活动性分析等理解都有着重要的指示意义。以上各类川滇地区动力学模式假说是否有着合理的观测资料论或实验数据支撑?是否有必要综合考虑多种动力学模式才能真实反映川滇地区现今的地壳形变特征?这些依然是需要讨论的问题。本文收集整理了已有的研究成果,旨在总结前人采用有限元数值模拟方法对川滇地区地壳形变动力学模式所展开的讨论,并指出了其中存在的问题及对未来研究进行了展望。

1 构造背景及岩石圈结构研究现状

过去的约50 Ma里,印度板块和欧亚板块之间的碰撞和持续挤压导致了青藏高原地区强烈的地壳形变、快速抬升,高原物质向东南方挤出,在川滇地区及邻区内发育了众多大规模的活动断层^{[3, 14, 34]}。不同走向的众多断层将川滇地区划分成了各级活动块体,而地壳变形则主要受控于这些块体边界的断裂带^[20],以此控制着研究区域应力应变^[35]。在川滇地区最早认为呈现出左旋走滑的鲜水河-小江断裂系是菱形块体的东边界,而表现为右旋走滑性质的金沙江-红河断裂系统则为菱形块体的西边界,在菱形块体内部北东向的丽江-小金河断裂将其分割成南、北两部分^[36]。文献[37]基于GPS约束反演得到鲜水河断裂是中国大陆内部活动最剧烈的断裂带之一,其平均滑动速率高达约12.0 mm/a,速率从西北部到东南部段逐渐下降。中部的安宁河断裂、南边相连的泽木河断裂以及新发育的大凉山断裂均有着相似的滑移速率约为4~5 mm/a,小江断裂界定了川滇菱形块体东南部边界,GPS数据显示其滑动速率由北向南逐渐减小,平均滑移速率约为8~12 mm/a。同时随着研究的深入,目前认为菱形块体的西边界并不清晰,主要因为红河断裂带表现出显著的分段特征,对其滑动速率的研究仍存在争议,且强震分布均明显弱于鲜水河-小江断裂。此外,北东走向的龙门山断裂是分隔松潘-甘孜块体和四川盆地的边界断裂。从历史地震分布看,7 级以上地震主要分布在菱形块体的边界断裂以及龙门山断裂等主要的大型断裂带上,这一块体划分模型较好地解释了强震活动的分布^[23]。但对于6 级左右地震在空间上则更为分散,使用仅有的几条主干断裂很难概括所有强震的分布区域,也有必要对川滇地区的断裂构造进一步梳理。

除了川滇地区的地震活动频度高、地形起伏大,地壳结构研究也表明区域地壳厚度和岩石介质变化也十分剧烈。接收函数结果^[25]表明,青藏高原东南缘地区的地壳厚度差异较大：松潘-甘孜块体的地壳厚度可以高达60~70 km,而云南南部的地壳厚度则只有30~40 km；龙门山断层以西地壳厚度约55 km,而以东的区域地壳厚度则迅速减小到约40 km；北川滇菱形块体下方地壳厚度在60 km以上,而南川滇菱形块体地壳厚度则变为40 km左右。同时,大量面波和体波层析成像结果表明^[38-40],在川滇菱形块体的北部以及南部的滇中地区均存在着低速层,这些低速层软弱物质的囤积导致了地壳增厚,在中下地壳所形成的软弱物质的通道也在不同程度上影响着地壳的形变。此外,体波成像结果^[41-43]表明,在扬子克拉通的西南部65~250 km深度范围内存在大范围的低速区域,这可能表明青藏高原的软流圈物质在南北向的挤压下也发生了增厚,并向东南方向的川滇地区逃逸,形成了与中下地壳流类似的地幔流,并且该地幔流的东南向运移可能导致了稳定的扬子克拉通的西缘在上地幔范围内出现了破坏,从而使得软流圈物质能顺利地向南运移。基于多数据联合反演方法,中国科技大学姚华建团队^[44]相继推出了区域高精度三维公共速度模型1.0和2.0版本（SWChina CVM-1.0和SWChina CVM-2.0）,随着地壳深部数据的累积,模型也将越来越精确,地壳模型的横向分辨率达到了0.2°。但中下地壳和上地幔部分的低速体空间分布形态仍然十分复杂,其深度和速度值也呈区域性变化,并且由于数据（台站数量、密集程度）和数据质量等以及方法不同,使得模型之间存在差异且缺乏系统对比和统一验证。

2 动力学模式分析

2.1   块体挤出模式

块体挤出模式的主要内容是：在印度板块向欧亚板块碰撞、楔入的过程中,其向北挤出的作用力导致川滇菱形块体物质向外挤出,受到稳定的华南块体阻挡,在发育形成了一系列大型走滑断裂带的作用下,各次级块体表现出了顺时针转动的运动特征。川滇地区的形变特征主要是以大型断层局部剪切为主导的板块形变来调节川滇地区内部的构造变形过程,最早有文献[3]利用塑性介质对印度板块和欧亚大陆的碰撞变形过程进行模拟还原了块体挤出过程,川滇地区的各级块体在这种碰撞作用下沿着区域内断裂带被挤出。

针对川滇地区众多复杂的断裂带,众多学者采用有限元数值模拟方法建立多断层模型来研究其相互作用影响,其优势在于可以根据断层真实具体的几何走向或形状灵活进行模拟,缺点在于存在较多的小规模、较活跃的断裂带,且对其几何走向和活动性质尚未达成共识。文献[45]建立了川西及邻区的三维有限元模型,以GPS观测数据为约束,主要针对鲜水河-小江断裂带建立了多组模型,从断层滑动速率的角度来讨论安宁河、则木河、大凉山断裂带在川西地区应变分配和地壳形变中分别起到的作用,认为安宁河-则木河断裂带与大凉山断裂带对鲜水河-小江断裂系统整体的应变分配等起到了一定程度上的影响。文献[46]建立了包含川滇块体、松潘甘孜块体和华南块体的三维有限元模型,着重针对鲜水河-小江断裂系进行模拟,认为大凉山断裂的形成和发育使鲜水河-小江断裂系形成平滑的几何形状,增加了鲜水河-小江断裂系大部分断层的滑移率,但减少了安宁河-泽木河断裂的滑移率。在断层系统中,一个断层的滑移率变化可以重新分配其他断层的滑移率,在运动学上连接并耦合的。断层相互作用显著影响藏东南缘的地壳挤压速率,川滇地块的挤压主要由其东部边界的鲜水河-小江断裂系和内部的丽江-小金河断裂带控制。文献[47]建立青藏高原东南缘的三维有限元模型进行了模拟,结果认为红河断裂较低的滑移速率并不能有效地调节区域形变,而是由实皆断裂的快速右旋走滑与鲜水河-安宁河-则木河-小江断裂系统共同形成一个大型走滑剪切带,该剪切带控制着区域浅部地壳层的变形过程,可以很好地解释现今川滇地区大地测量手段所观测的走滑、旋转、拉伸和挤压的地壳变形过程。文献[39]建立川滇地区的三维层状弹性有限元模型模拟了断层滑动速率和应变速率场,结果表明藏东南缘以剪切应变为主,在小江断裂带西侧应变速率最集中,与该区域的构造挤压变形模型一致,结合刚性块体模型和连续块体变形的构造模型可以最好地描述川滇地区的地壳形变。

在有限元模型的建立时,首先需要涉及块体划分和断层选取的问题,地块边界带由一系列活动构造断裂带所构成,总体形成该地区重要的活动断裂体系,均需纳入有限元模型框架中,但川滇地区发育有众多复杂的断裂带,前人建立的模型中大多只模拟了青藏高原东南缘一系列大尺度的走滑断层,并未能反映局部的应变细节,许多次级断层不包括在内。川滇菱形块体作为藏东南缘块体划分的核心区域,其东边界由甘孜玉树断裂、鲜水河断裂、安宁河-则木河断裂及小江断裂组成,这一观点目前已基本达成共识,但对于西部边界仍不明确或存在争议：有学者认为菱形块体的西北部边界并不存在,即该块体的西北边是不封闭的,与往东南方向运动的西藏块体相连；也有学者通过野外地质调查结果认为金沙江断裂或者澜沧江断裂是菱形块体的西北部边界,但目前尚未有定论。对于位于川滇地区西南位置的红河断裂带在很长时间以来被众多学者认定为川滇菱形块体的西南边界^[37,48-49],但红河断裂带中部和南部均没有历史强震的记录,并且在其中部出现了显著的地震空区^[50],在红河断裂带东西两侧还存在一系列近似平行的断层,如无量山断裂、楚雄-建水断裂及强震活动活跃的龙陵-澜沧断裂,究竟该如何进行块体划分对有限元模型建立至关重要,影响着区域构造运动加载和应力分配。因此需要结合多学科研究成果酌情选择合适的断层进行块体划分,分析断层间相互作用对区域构造形变的影响。

2.2   中下地壳流模式

中下地壳流模式的主要内容是：在青藏高原的挤压作用及地形的差异下,低粘度的中下地壳物质从高原向着东/东南方向迅速流出,而较坚硬的上地壳和地表则伸展缓慢,即大陆岩石圈的强度呈现为上地壳强、中下地壳弱和上地幔较强的“三明治”结构^[51]。这一模式在众多学者的数值模拟研究过程中对于地表GNSS速度场拟合效果较好,也能较合理地解释现今地表的主要变形特征和形成机制。文献[29]最早采用了数值模拟方法和下地壳流模型解释了高原的形变特征,由于较软下地壳导致上地壳的运动变形和地慢的运动解耦。早期对于中下地壳流模式的数值模拟大多通过测试下地壳流的粘滞系数或者调整中下地壳部分侧面边界条件和底部拖曳作用来实现,通过和地表GNSS观测的水平速度场的对比判断下地壳流模式是否合理。文献[52]运用遗传算法反演了川滇地区三维有限元模型的边界条件,发现在模型下地壳底部施加拖拽作用后能有效改善对GPS实测观测结果的拟合。文献[53]在考虑地形和莫霍面界面形态的情况下,利用三维粘弹性模型测试了数值模型的动力学边界,对模型底部施加拖曳作用后能更好地拟合观测结果；文献[54]利用粘弹性模型模拟了青藏高原的变形过程,模拟结果显示,青藏高原在印度板块的推挤下,软的下地壳使得高原整体隆升,柔软的下地壳和软流层物质在这一驱动力作用下向东南流动,并拖曳川滇地区上地壳运动。文献[55]针对川西地区建立了三维粘弹性模型,探讨了川西地区下地壳流动与地表变形之间的力学关系,对下地壳边界条件进行多次实验结果表明,在下地壳比上地壳流动快约11 mm/a时,模拟结果与GPS观测数据吻合较好,可用于了解川滇地区地表变形的可能机制；文献[56]在数值模拟川滇地区形变场时,考虑下地壳拖曳作用,并得到拖曳力主要集中在川滇菱形块体内,在块体北部拖曳力以南向为主,在块体南部拖曳力逐渐减小并且方向上逐渐向东西两侧分流。

从数值模拟角度讲,“中下地壳流模型”实际是地壳介质流变性差异的具体表现,中下地壳层被多数学者认定为一个低粘度的均匀层状结构,而一些研究工作也提出了川滇地区下地壳低粘度物质是以通道流的形式运动。文献[28]在中下地壳流动模型的基础上提出管道流模型用于解释青藏高原周缘地形梯度的差异。文献[57]构建了三维有限元模型模拟了3年汶川地震周期不同阶段的地壳变形,为了研究中下地壳流的形式和大小,在川滇地区西侧的中下地壳设定了一个圆形区域,并在地表以下40 km的节点上输入了恒定速度的水平物质流动,模拟更快的物质流动速度会导致区域更快的隆起速率,随后对下地壳流模式进行模拟,并搜索得到最佳的流域半径和物质流速组合分别为200 km和45 mm/a,认为仅靠刚性块体挤出模式不足以解释青藏高原东部快速隆升率,而当中下地壳物质流动速度从30 mm/a增加到50 mm/a时,地表的平均垂直隆起速率从约1 mm/a增加到约3 mm/a,地壳缩短模型结合中下地壳流动模型可以很好地解释青藏高原东部的快速隆起。相比于规则的“三明治结构”,随着地震层析成像的结果丰富,众多学者认为较弱的中下地壳层以部分不规则的熔融状态存在,文献[58]认为弱地壳层的存在只能归因于额外局部的熔融或流体,即中下部地壳中的部分熔融形成低速体,这可能是由于高温和含水矿物的脱水造成的,地壳流处于中地壳浅得多的深度（约10~40 km）,并且没有达到莫霍面。文献[59]根据川滇地区地震层析成像的研究成果,发现川滇地区下地壳部分存在两个局部的低速体,即存在部分岩石熔融和粘度降低的情况,采用调整下地壳岩石介质参数中S波速度参数的方法轮廓绘制了低速体的形状加入到三维有限元模拟中,模拟结果表明,当模型同时加入断层和低速体模型进行模拟时,可以在大多数位置更好地拟合地壳形变运动趋势,认为块体挤出模型和下地壳流模型并不排斥,在浅部上地壳的运动主要由断层控制,而下地壳的变形则是由局部剪切带和下地壳流因素共同驱动,都是解释川滇地区构造模式的必要条件。

在中下地壳流理论提出后,采用有限元数值模拟方法一般将地壳的流变性质简化为有效黏度,并大多采用层状结构来划分地壳各层,目前众多学者建立的地球动力学模型中选取的低粘滞系数存在显著差异。文献[28]采用管道流理论模拟的下地壳流粘滞系数约为1×10¹⁸ Pa·s；文献[53]将青藏高原下地壳作为整体来模拟下地壳粘滞系数和流速,根据对比试验下地壳的流动速度比地表运动速度快约10 mm/a,其粘滞系数为1×（10¹⁸~10¹⁹） Pa·s；文献[59]建立川滇地区三维有限元模拟,多次测试得到下地壳流部分粘滞系数约为（1~3）×10²⁰ Pa·s。不同模型对于低粘度的中下地壳的范围设定的不同可能是造成粘滞系数差异的主要原因。而根据接收函数、大地电磁测深、地震层析成像等多种方法对该地区岩石圈结构的分析,川滇地区地壳内存在几何形态较复杂的大范围低速层。文献[60]认为在川滇地区中下地壳内存在两个低速体通道,一个低速体从松潘-甘孜地带穿过丽江和大理盆地,再到云南西部,另一个低速体沿小江断裂带延伸,穿过红河断裂带到思茅地块,然而这些低速体是否相互连接、来自何处,受限于观测数据和手段,目前尚无法对这一科学问题给出明确的结论,因此对中下地壳流变性质的规模和尺度有待于深部地球物理结果的发展和丰富,还需得到更可靠的中下地壳流范围为有限元模型的数值模拟过程提供约束。

2.3   重力扩散模式

重力扩散模式的主要内容是：青藏高原及周边地区的地壳变形受到高原的过度重力势能驱动,青藏高原在受到印度板块向北强烈挤压的作用下,岩石圈地幔形成对流剥离和拆沉后由于负浮力而抬升,由于地壳增厚的不均匀性而形成重力势能向高原逐步朝外扩散。文献[61]对青藏高原现今构造应力状态及构造运动的三维模型的模拟结果显示,青藏高原巨大的重力势能是造成川滇物质往南东运动的重要动力。川滇地区地表高程在西北侧靠近高原的地区最高达到4 km以上,而向东南方向迅速降低,直至四川盆地地表高程约为500 m,在这样陡峭的地形海拔差异下,重力是川滇地区数值模拟过程中有必要考虑的一个因素。文献[62]建立了青藏高原地区的三维粘弹性模型研究了导致高原扩散的主要因素,研究在当前地形和GNSS速度场运动趋势的条件下,模型预测了高原地壳内的拉伸应力状态,这是由高原重力与构造压应力之间的平衡造成的,经过多次实验模拟,当青藏高原地形海拔降低到目前海拔的约50%时,走滑断层和逆冲断层在整个高原的形变趋势上占主导地位,只有当青藏高原海拔达到与目前海拔高度的约75%时,才会发生显著的地壳扩张。重力扩散模式也会在岩石流变性质的影响下共同对运动变形产生影响。文献[30]利用三维球壳模型模拟了印度和缅甸俯冲板块与地形引起的重力势能对青藏高原地区的影响,首先设定印度俯冲板块的粘度约为1×10²² Pa·s,下地壳的粘度约为1×10²⁰ Pa·s,除了川滇地区之外,它整体比较能拟合地表速度数据,随后减低了喜马拉雅东造山结附近的上地壳粘度,由于重力势能引起的地形坍塌,加大了地表速度的顺时针旋转。

除了上述对大尺度青藏高原的大陆俯冲和碰撞进行模拟外,也有聚焦于川滇地区局部的重力扩散模式的模拟。文献[63-64]利用热力学数值模拟建立了川滇地区模型对其东南向挤出的地表运动特征进行了动力学解释,该区域变形机制更合适用重力驱动来解释。文献[65]利用三维接触有限元模型模拟了青藏高原东南缘所受的动力学环境,研究结果表明川滇地区在下地壳的拖曳和重力共同作用下,高原物质向东南部溢出。文献[66]建立了川滇地区的三维粘弹塑性有限元模型,并引入了重力因素,探讨了印度板块碰撞下侧向驱动力在青藏高原东南缘地壳形变起到的作用,首先在不考虑边界构造驱动力时,对模型施加重力作用,模型整体已能够粗略模拟出现今川滇地区的地壳形变趋势,认为观测到的当今地壳运动在很大程度上可以用重力扩散来解释；然后在模型域的北部和西北边缘施加20 mm/a的速度来模拟印亚碰撞的横向驱动力,此时并未施加重力影响,模拟结果仅在川滇地区西北侧能粗略模拟出现今的运动趋势,而其他地区则并未有明显运动,因为在构造挤压模型中剪切应变沿着主要的走滑断层扩散,而运动形变在南部已基本消失,并未能反映出顺时针旋转的趋势；最后将碰撞横向驱动力和重力因素共同施加到模型上,模拟结果基本上拟合了现今GNSS速度场的主要特征,认为青藏高原东南部的当今地壳变形在很大程度上可以用重力扩散来解释,仅靠构造挤压模型中描述的边缘推力无法解释青藏高原东南部观察到的应变率划分和应变传递。

随着近年来的数值模拟研究工作中,地形载荷和重力加载因素逐渐凸显其显著的影响,尤其在藏东南缘地区,侧向压力梯度导致西藏东南部地壳明显向东南方塌陷,再加上中下地壳流效应引起局部物质快速流动,从高地向低地流动。但目前已有的数值模拟研究对模拟结果大多以GNSS地表水平速度场的拟合程度作为评判标准,均认为重力驱动的地壳变形主要依赖于地壳黏度,尤其是上地壳的粘度较低,引起的水平速度较高；而降低中下地壳的粘度引起的水平形变差异并不大。造成这种现象的主要原因可能是中下地壳低速体通道中的物质流动的水平变化很难反映到地表的水平形变上,因此这一重力驱动过程是否有中下地壳流参与,目前已有的研究中并不确定。文献[67]基于川滇地区高精度三维大地测量数据特别是垂直方向速度,通过三维有限元建模定量研究重力、构造挤压和中下地壳流对藏东南地区当前地壳变形的影响,结果表明上地壳黏度对垂直变形的影响较小,地表垂直速度主要取决于中下地壳的流变性,深层地壳中的快速流动促进了更大的地表隆起（约5 mm/a）。因此要判断重力驱动过程中影响范围可能还需要更多三维大地测量数据及丰富的深部地球物理数据,综合水平形变、垂向形变、地壳各向异性数据结果等来提供进一步约束,才能有助于了解藏东南地区岩石圈变形特征。

2.4   地幔对流模式

不同于以上3种构造变形模式基本聚焦于青藏高原东南缘地壳部分进行数值模拟,前人的工作中往往忽略了深部的软流圈运动对地表的影响,地幔物质具有高粘性,在岩石圈重力作用的影响下,深处热的上地幔岩石相较于地壳部分冷的岩石便可能发生地幔对流。从大尺度的青藏高原受到的板块俯冲作用来看,理论上青藏高原和印度板块长期向北运动趋势所造成的软流圈剧烈运动有可能对喜马拉雅东构造结和青藏高原向东南缘挤出的川滇地区地壳形变产生影响^[67]。因此地幔对流模式的主要内容是：川滇地区的地壳构造运动可能会受到软流圈的对流的拖曳作用。总的来说,地幔对流模式是对以上3种浅部模式的一种补充和完善,更适用于较大范围区域的青藏高原整体动力学模型进行模拟讨论。理论上若清晰地了解了地壳上地幔岩石介质的密度、粘滞度和泊松比等物理参数及岩石圈内部各种物理化学反应过程,通过求解复杂的控制方程组便可以获得地幔运动的清晰图像,因此当前利用数值模拟方法研究地幔对流模式对地壳形变的作用渐渐得到重视,基本是通过模拟求解上地幔对地壳部分的拖曳力或运动的形式来实现。

文献[69-70]利用三维热力学数值技术模拟了大陆碰撞-海洋俯冲和板片撕裂过程,这些模型表明,由于海洋板块俯冲后退和板片撕裂引起的地幔流可以拖拽上覆板运动完成旋转变形。文献[71]针对青藏高原地区进行了三维地幔动力学模拟,结果表明印度板块之下软流圈地幔对岩石圈向北的拖曳力可能是造成高原隆起的驱动力的主要来源,在碰撞边界下,印度岩石圈向北俯冲带动了地幔软流层运动,这种软流圈流动是川滇地区的构造挤压和喜马拉雅东构造结地形抬升并从西藏东南侧向溢出的作用力之一。文献[72]对青藏高原地区地壳运动的模拟实验结果显示,同时考虑板块推挤和地幔对流施加于岩石层底部的拖曳力这两种载荷条件时,模拟结果对地表GNSS观测结果拟合效果更佳。地幔对流拖曳力对青藏高原内部的地壳运动方向有着显著的影响,并且由于下地壳和上地幔中低粘度的岩石区域分布并不均匀,导致不同区域的壳幔耦合程度强弱存在显著差异,在局部的熔融或流体区域可能存在解耦现象,而对于没有低速体的地方,地幔对流模式的拖曳作用会对地壳的构造变形模式产生一定的影响。目前对川滇地区局部精细的地幔对流模式研究还存在空缺,虽然对于青藏高原大尺度的碰撞-俯冲模型已可以较好模拟整体地幔对流趋势,但缺少对局部区域细节的考虑,很难反映川滇地区精细的拖曳力分布细节。

3 现今动力学模拟存在的问题与挑战

综上所述,现有研究表明川滇地区地壳运动的动力学模式的驱动因素包括：在构造加载下断层的调节作用、中下地壳流变性质和岩石介质力学强度差异、地形差异引起的重力扩散和地幔对流在地壳底部施加的拖曳力作用^[47]。早期针对川滇地区的动力学数值模拟研究大多针对印度板块与欧亚板块相互作用和较软的中下地壳流变性质展开讨论,其中板块间力的相互作用在模型里体现为GNSS速度场对研究区域的外部约束,在模型内部主要断层和块体划分基本由现有大地测量学对地表形变观测反演得到,而对于深部的中下地壳流模式,通常采用规则的平板模型来模拟较软的中下地壳部分,从而导致中下地壳的厚度以及其粘滞系数的取值出现较大差异,随着地球物理探测资料的丰富^[25],能够明显地看出地壳深部介质属性的强弱,使得准确获取中下地壳内低速体形态成为可能,但其粘滞系数还需要根据拟合效果分析获得。而川滇地区的重力因素和地幔拖曳作用近年来越来越受到动力学研究的重视^[66],在模拟工作中忽略地形起伏因素引起的重力扩散作用或者忽略地幔对流作用对地壳形变的影响可能导致块体挤压模式和中下地壳流模式的贡献被高估,因而可能得到错误的次级块体划分结论或者获得存在较大误差的低速体粘滞系数^[69-70]。因此要透彻理解川滇地区复杂的地壳动力学机制,有必要综合考虑上述多种驱动机制,才能厘清各类外部和内部动力学模式对地壳形变的贡献。

目前,学者在川滇地区数值模拟研究中大多从某一种或者两种模式展开针对性讨论,这主要是因为前文中介绍的4种动力学模式各自在作用的空间尺度和时间尺度上存在差异。如对于块体挤出模式而言,其作用范围大多集中于表现为弹性的上地壳部分,但是在更长的地质时间尺度上,即便是上地壳也呈现出粘性流体的运动特征,对于中下地壳流模式和重力扩散模式,不仅仅空间尺度上作用范围扩展到整个地壳范围,在时间尺度上也需要更长的模拟时间（数万年至数十万年）才能使得模型完成演化达到稳态^[73-74]。此外,大部分现有的数值模型以观测到的GNSS地壳水平运动为约束,由于大地测量学观测到的现今地壳变形观测只有短短几十年,现今的形变并不代表过去几百万年能一直保持这样的运动趋势,块体挤出模式也许并不能用现今短短几十年的大地测量观测来确定,而地质学得到的数据如断层滑动速率等则是更长期的结果。对于一种正确变形机制,理论上讲它应该能够用所有的观测来证明,而不仅仅是地表现今的水平变形,而对于地壳垂直变形、地形、断层地质滑动速率、长期重力变化、SKS、SWS等地震学观测则考虑较少。因此需要充分考虑模拟结果和计算效率的平衡,逐步优化有限元模型的构建和网格划分细节,综合多种动力学模式来设置模拟方案,有必要结合多学科、不同时空尺度的观测资料进行约束,以验证模型的有效性。

此外,好的模型是数值模拟成功的基础,在实际建模过程中,大多数模型为了针对某一待解决的具体科学问题,往往对模型的几何特征和边界条件进行了较大程度的抽象,弱化了与研究兴趣无关的部分物理参数的变化对模拟的影响,比如采用简单的断层几何模型、固定的运动速度或应变速率作为边界条件或者简化内部物质物理化学属性等。这使得模型能够有针对性地对感兴趣的参数影响或者某些约束条件下模型变形的整体趋势的宏观变化规律进行分析,但这类模型不可避免地会导致模型对实际过程失真。对于更加复杂的川滇地区动力学模拟过程,模型中各类参数或影响因素的关联性更加紧密。由于川滇地区内除龙门山断裂带外,大多数断层均表现为走滑性质,断层倾角较大,在模型中便普遍将这些断裂带倾角简化为垂直,但其几何形态的细微变化往往能很大程度地影响模拟结果,目前已有由中国地震局地质研究所等单位给出了川滇地区（中国地震科学实验场）主要活动断裂的三维模型V2.0+版本等断层数据,使建立高精度三维断层模型成为可能；同时也已有诸多文献给出了全球或区域性的岩石圈速度模型,包括川滇地区精细三维速度结构数据及地壳各界面的埋深数据,目前姚华建团队^[44]给出的川滇地区公共速度模型的分辨率已达到0.2°×0.2°,使得模型介质更接近真实的岩石物性。但过于复杂的模型和参数往往使得计算量和模拟时间极大地增加,对计算效率和数值模拟算法等提出了更高的考验,如何权衡取舍需要在实践中进行检验。

4 总结与展望

川滇地区一直都是地球科学领域探讨的一个热点地区,其复杂的动力学模式仍是一个尚未解决的科学问题。在过去的几十年里,针对川滇地区构造变形动力学模式的研究取得了丰富的成果,本文聚焦于在该区域采用有限元数值模拟方法所获得的研究结论,在川滇地区从地壳到上地幔可能正在同时经历4种不同模式的驱动力：在高度差异的地表海拔起伏下,重力因素驱动着浅层地形流动；在上地壳,刚性块体沿着主要的走滑断层运移,呈现刚体挤出模式；在中下地壳存在大范围的壳内低粘性物质,在区域应力场的控制和影响下流动,在局部区域符合中下地壳流模式；到了上地幔,南部的印支块体发生了显著的岩石圈拆沉减薄以及软流圈物质上涌,并逐渐影响到扬子克拉通的西缘,地幔在岩石圈底边界施加的拖曳力驱动着整体形变。在前人的有限元数值模拟工作中,以上4种动力学假说都得到了一定程度上的验证,各自都有其特定适用的空间尺度和时间尺度,但同时每种动力学模式目前在认识和理解上仍存在一定的局限性。已有的数值模拟研究中表明以上多种动力学模式并不互斥,而是共同驱动造成了现今川滇地区的形变特征,综合考虑川滇地区的多种动力学模式是今后数值模拟研究的趋势。

总得来说,受当前的观测数据和数值模拟技术所限,对未来的发展趋势将聚焦于以下3个方面：

1）多物理场耦合模拟动力学过程。目前大多数川滇地区有限元模型基本是通过简单的物理-力学模型来进行运动及形变分析,模拟地壳和地幔在地质时间尺度上的变形,通过将岩石圈简化为基于弹性-粘弹性模型来展开动力学数值模拟讨论。但对于在地壳深部的效应是一个包含化学反应、热力学效应、流体流动和应力变形等多物理场耦合的复杂动力学系统,这些变形的过程中通常伴随着热传输、地球深部的相变、复杂的流变学行为、熔融和熔体迁移、化学反应以及固体运动等现象^[75]。川滇地区中下地壳流的存在表明该地区地壳深部岩石处于熔融或部分熔融的高温状态,中下地壳流的高温物质也为浅部脆性地壳发生局部熔融提供了热量,因此在地壳可能存在着一种动态的热力学-流体流变学的变化,简单用通道流或者单一的粘滞系数来表征低速体的流变性质可能会存在一定的争议。因此,有必要在传统的固体地球动力学研究的基础上,结合地壳深部的化学反应和流体流动等研究多物理场耦合的综合动力学过程,对川滇地区动力学驱动效应进行多时空尺度定量化表征和分析,这将有助于更全面地理解川滇地区动力学模式的复杂性,揭示其内在机理和演变规律。

2）多学科资料综合应用。现有的地表大地测量学观测结果已基本揭示了川滇地区地表运动特征,而对整体岩石圈动力学模式的数值模拟讨论更着眼于深部数据的丰富,数值模拟的模型建构、运动学分析和动力学模式的讨论等过程都需要充分考虑运动学模型、板块构造模型、活动地块模型、应力-应变场、岩石流变模型、地壳介质模型和震源模型等数据资料。随着地球物理学、地震学、构造地质学、重力学和地热学等多学科等对地壳深部观测的理解和认识逐步提高,多学科交叉的实测资料会为数值模型的建立提供更为丰富的数据,将数值模拟结果与实测资料进行各方面的对比,使数值模拟更具实际的科学意义。完善并融合川滇地区的多学科资料数据,共同丰富动力学数值模型的构建过程和模拟结构后处理验证的可靠性,深化川滇地区岩石圈内部驱动过程的研究和理解,是未来真正解决川滇地区动力学模式这一重要科学问题的希望所在。

3）高性能计算技术及算法效率和稳定性的提升。有限元数值模拟因其高度的灵活性和适应性,在处理复杂地质结构和边界条件方面表现出色,尤其适用于大尺度和长时间尺度的模拟。随着川滇地区各学科资料的丰富,其所建立的三维有限元模型的几何形态更为复杂、模型分辨率逐渐提高,其自适应网格划分也更加密集,这便导致了有限元模型单元数量激增,模拟计算量大大增加；同时随着各类几何和物理（结构和材料）非线性问题的出现,计算求解的收敛难度也是当前所面临的一个重大考验。随着现代计算机能力的逐渐提高,动力学过程的数值模拟也由简单变得复杂,使模拟过程更真实地反映深部地壳的物理过程。在此背景下,为了适应更加复杂的模拟需求,模拟算法必须具备更高的精度和稳定性。在计算机技术高速发展的今天,随着当前各类硬件和软件的开发和完善,有限元数值模拟方法在地球科学领域发挥着日益重要的作用,未来需要再进一步优化程序内部架构,实现程序的大规模并行计算,提高数值模拟的效率和精度,最终有效实现动力学数值模拟过程。