近50 Ma以来，印度—欧亚板块持续碰撞不仅塑造了现今高海拔的青藏高原(Molnar and Tapponnier, 1975; England and Houseman, 1986)，同时驱使高原物质向东扩展.受刚性四川盆地的阻挡，在青藏高原东缘形成大规模北东-南西向的逆冲推覆构造(邓起东等, 1994; Tapponnier et al., 2001).松潘—甘孜地块和四川盆地之间地形最为陡峻的龙门山断裂带(王二七和孟庆任, 2008)，是青藏高原东缘典型的构造活跃区，2008年M_S8.0汶川地震和2013年M_S7.0芦山地震都发生在龙门山断裂带上.青藏高原东缘作为高原物质向东扩展的前缘过渡带，研究其构造特征和变形机制有助于更好地认识青藏高原侧向扩展模式及高原地块边缘的发震机制.

晚新生代以来，龙门山快速隆升形成显著的陡变地形起伏，主要表现为松潘—甘孜地块至四川盆地近30 km范围内高达4 km的地形差异.而晚新生代以来龙门山地区的地壳缩短量不足以产生高原边缘如此大的地表抬升和地壳增厚(Burchfiel et al., 1995; 王二七和孟庆任, 2008).针对青藏高原东缘的陡变地貌特征和隆升机制，主要存在两种典型动力学模式：一种为逆冲断层滑动模式，即龙门山地区将高原东缘的扩展挤压作用以高倾角逆冲断裂活动形式转化为较大的垂向抬升驱使龙门山迅速隆升形成高陡的地形(Tapponnier et al., 2001; Hubbard and Shaw, 2009; Wang et al., 2011; Zhang, 2013); 另一种是下地壳流模式，即青藏高原柔性中-下地壳(部分熔融)向东流动受刚性四川盆地的阻挡驱使地壳物质抬升，形成高陡地形(Bird, 1991; Royden et al., 1997, 2008).随着中国大陆地表形变观测的不断发展，大量高精度GPS观测和水准观测提供了青藏高原东缘现今丰富的三维形变数据.区域GPS水平速度场显示松潘—甘孜地块至四川盆地现今水平缩短变形较弱(小于3 mm·a^-1)，GPS垂向速度观测及区域精密水准观测数据显示龙门山断裂带西北侧的松潘—甘孜地块地表隆升速率为2~4 mm·a^-1，而东南侧的四川盆地呈现0~2 mm·a^-1的相对下沉趋势(Shen et al., 2005; Hao et al., 2014).显然，相对于龙门山两侧高达6 mm·a^-1的地表垂向变形差异，区域现今水平缩短速率同样不足以产生龙门山断裂两侧如此显著的地表垂向变形量(Shen et al., 2005; Liang et al., 2013; Hao et al., 2014).

青藏高原东部岩石圈深部探测结果显示地壳结构在龙门山地区存在明显过渡：松潘—甘孜地块地壳厚度约为60 km，向四川盆地逐渐减小至40 km (楼海等, 2010; 陈石等, 2013; Zhang et al., 2010).龙门山断裂带两侧地壳地震波速结构也存在较大差异，相对于刚性的四川盆地，松潘—甘孜地块具有较弱的力学性质，中地壳低速层和高导层等地球物理观测证据指示区域可能存在地壳流(嘉世旭等, 2014; Bai et al., 2010; Wang et al., 2007, 2008; Liu et al., 2014).Wang等(2018)利用P波接收函数方法分析认为松潘—甘孜块体地壳结构在空间上具有很强的横向不均匀性，块体内部的垂向增厚和块体间强的相互作用共同主导了地壳的变形.另外，青藏高原东缘地壳密度结构的相关研究表明松潘—甘孜地块地壳密度整体低于四川盆地，其中松潘—甘孜地块上地壳存在约5 km厚的低密度层(Zhang et al., 2014)，可能对应部分熔融的地壳流.布格重力数据显示松潘—甘孜地块处于负异常的重力非均衡状态(Zhang et al., 2014).

青藏高原东缘松潘—甘孜地块至四川盆地地壳结构和力学性质的空间差异可能对区域变形具有重要影响.据此，本文结合区域岩石圈结构、密度分布、地形起伏和地表三维形变等资料，建立横跨龙门山断裂带的二维有限元岩石圈剖面模型，分析地壳密度结构、陡变地形和岩石圈流变性质对区域变形的影响，进而探讨青藏高原东缘现今垂向变形机制.

为计算青藏高原东缘陡变地形和地壳密度结构空间差异控制下的区域岩石圈变形，本文采用已广泛应用于模拟地质构造演化和地壳流等问题的黏性流体模型(Bird et al., 1991; Clark and Royden, 2000; Bendick and Flesch, 2007; Copley, 2008)，利用有限元开源程序Gale (https://geodynamics.org/cig/software/gale/)来计算分析青藏高原东缘岩石圈变形问题.该软件通过引入单元内追踪点技术实现模型密度和黏滞系数等力学参数随时间和位置的变化，可有效模拟岩石圈的变形演化等问题.由于温度场和岩石圈流变性质在短时间尺度内变化较小，本文采用简化的牛顿黏性流体控制方程，对应岩石圈各分层流变结构采用等效黏滞系数.

主要控制方程如下：

(1)

(2)

(3)

(4)

其中，σ_ij为应力张量的分量，P为压强，τ_ij为偏应力张量的分量，δ为克罗内克函数，ν_{i, j}为速度梯度，η为等效黏滞系数，f_i为体力项.

研究区由西向东依次包括松潘—甘孜地块、龙门山断裂带和四川盆地.参考龙门山周缘水准观测到的地表垂向变形结果，选取接近水准观测点沿线且横跨龙门山断裂带的一条剖面(如图 1所示)，建立了二维有限元模型.模型(北西-南东向)总长390 km，由左向右分别为松潘—甘孜地块、龙门山和四川盆地，深度上由地表延伸至岩石圈地幔110 km处.模型考虑地壳各分层界面起伏、地形起伏及地壳力学结构的横向差异等因素.地震学研究结果显示青藏高原地壳温度较高，而四川盆地地壳温度较低(Liang et al., 2012; Shen et al., 2015; 沈旭章, 2013)，因此模型中松潘—甘孜地块相对于四川盆地具有较低的岩石圈黏滞系数(见图 2b).模型密度参考Zhang等(2014)给出的青藏高原东缘岩石圈剖面精细密度结构，其中，参考松潘—甘孜地块中地壳低密度层的分布，对模型中松潘—甘孜地块地壳20km深处设置了约5km厚的低密度层，并赋予较低的等效黏滞系数，如图 2b所示.

图 1

Figure 1

图 1 青藏高原东缘构造背景和地表形变特征 白色箭头为水准观测获取的垂向速度(Hao et al., 2014)，黑色箭头为GPS垂向速度，灰色箭头为GPS水平速度分布(Liang et al., 2013)，插图(b)中黑色线框为研究区位置. Figure 1 The geological settings and surface deformation in Eastern Tibet The white arrows are vertical velocities from precise leveling observations (Hao et al., 2014), the black arrows are vertical velocities from GPS observations, the gray arrows are horizontal velocities from GPS observations (Liang et al., 2013), the black rectangular in the inset (b) is the location of the research domain.

为对比分析地壳密度的横向差异对青藏高原东缘岩石圈垂向变形的影响，本文设定了三组对比模型，分别给定三类驱动力：构造驱动、重力驱动和重力及构造共同驱动的数值模型，如图 2(c—e)所示.其中，(1)构造驱动的模型(Model Rt，图 2c)，不考虑重力，参考龙门山两侧GPS速度场资料，对模型左侧边界施加8.3 mm·a^-1的水平速度，右侧边界施加7.0 mm·a^-1的水平速度，水平速度不随深度变化，两侧边界切向自由；(2)重力驱动的模型(Model R, 图 2d)，考虑重力，模型的变形驱动力为地壳密度横向差异造成的区域压力梯度；模型两侧边界水平速度固定为0，切向自由；(3)构造和重力共同作用的模型(Model Rtg, 图 2e)，考虑重力，模型变形驱动力为构造加载和重力调整共同作用；参考区域GPS速度场，对模型两侧边界施加水平构造约束(同Model Rt模型)，切向自由.三组对比模型的上边界地表均自由；底边界垂向速度约束为0，切向自由.此外，为分析重力驱动作用下岩石圈流变性质和地形起伏对模型变形的影响，在Model R模型的基础上，本文通过改变松潘—甘孜地块中地壳和模型岩石圈地幔黏滞系数及不考虑地形起伏设置了多组对比模型，各模型描述见表 1.

表 1 Table 1

表 1 计算模型参数设置 Table 1 Model parameters 计算模型 地形 起伏 Moho起伏 构造 加载 重力 效应 中地壳黏滞系数 (1020Pa·s) 岩石圈地幔黏滞系数 (1021Pa·s) Model Rt √ √ √ × 1.0 1.0 Model R √ √ × √ 1.0 1.0 Model Rtg √ √ √ √ 1.0 1.0 Model 1 √ √ × √ 0.01 1.0 Model 2 √ √ × √ 0.1 1.0 Model 3 √ √ × √ 10 1.0 Model 4 √ √ × √ 1.0 5.0 Model 5 √ √ × √ 1.0 10.0 Model T0 √ √ × √ 1.0 1.0

表 1 计算模型参数设置 Table 1 Model parameters

图 2

Figure 2

图 2 二维岩石圈剖面模型 (a)横跨龙门山断裂的AA′剖面(见图 1)地形和Moho起伏；(b)模型(Model R)岩石圈分层黏滞系数和密度分布，密度数据参考(Zhang et al., 2014)；(c)—(e)分别为构造加载(Model Rt)、重力效应(Model R)、构造加载及重力效应之和(Model Rtg)对应的模型变形驱动力示意图. Figure 2 Model setup of the 2D lithospheric profile (a) The topography and Moho of profile AA′ (in Fig. 1) across the Longmenshan fault; (b) The distribution of viscosities and densities in the layered lithosphere, with the crustal densities acquired from (Zhang et al., 2014); (c)—(e)The diagrams of driving conditions from tectonic loading (Model Rt), gravity adjustment (Model R), and the combination of tectonic loading and gravity adjustment (Model Rgt), respectively.

稳定的GPS速度场可有效反映区域构造变形，为数值模型提供可靠的构造边界约束.构造加载作用下的模型(Model Rt)变形结果显示，在水平挤压缩短控制下近地表呈现微弱的抬升现象，如图 3所示.水平速度场由松潘—甘孜地块向四川盆地方向呈递减趋势，由于模型力学结构的非均匀性，模型水平速度场的递减梯度也呈现空间上的非均匀分布；其中，松潘—甘孜地块中地壳(20 km深处)黏滞系数较低的柔性层，相对于其上覆及下部地壳和岩石圈地幔，水平运动速度较快.垂向运动表现为近地表抬升速率由松潘—甘孜地块向四川盆地递减，且由浅至深抬升速率不断减小；其中，松潘—甘孜地块地表垂向隆升速率约为0.65 mm·a^-1，而四川盆地地表抬升速率约为0.2 mm·a^-1.由此可知，构造加载作用下的区域抬升量微弱，地表抬升速率明显小于青藏高原东缘现今观测到的垂向变形速率，同时水平构造驱动产生的垂向变形的空间差异也明显小于实际观测值.

图 3

Figure 3

图 3 水平构造驱动作用下模型变形分布 (a)模型水平速度分布；(b)模型垂向速度分布；箭头代表速度方向. Figure 3 Model deformation driven by horizontal tectonic loading (a) The horizontal velocities; (b) The vertical velocities of the model, the arrows indicate the directions of model velocities.

图 4给出了仅考虑重力作用的模型(Model R)变形结果，该模型驱动力来源于地壳非均匀密度结构和陡变地形起伏造成的横向压力梯度.岩石圈横向压力梯度驱使模型深部岩石圈地幔由四川盆地向松潘—甘孜地块呈现近水平方向的快速流动，水平运动速率高达5 mm·a^-1.受陡变地形起伏的影响，松潘—甘孜地块东边缘至龙门山地区低密度力学性质较弱的中地壳向四川盆地方向流动，水平运动速率高达2.0 mm·a^-1，可认为是陡变地形驱动产生的局部地壳流动.而模型内松潘—甘孜地块西侧和四川盆地地壳范围及近地表的水平运动速度较小.在重力效应驱动下，龙门山至四川盆地呈地壳规模的下沉趋势，近地表最大下沉集中在四川盆地东南侧，速率高达1.0 mm·a^-1；而龙门山西北侧的松潘—甘孜地块呈现地壳规模的整体抬升现象，近地表最大上升速率约2.0 mm·a^-1.重力效应造成龙门山两侧地块地表近3 mm·a^-1的垂向运动速度差异，而水平构造加载造成的模型两侧垂向运动速率差异小于0.5 mm·a^-1，由此可知，青藏高原东缘地壳密度结构差异对龙门山两侧垂向变形影响显著.

图 4

Figure 4

图 4 重力驱动作用下模型变形分布 (a)模型水平速度分布；(b)模型垂向速度分布，箭头代表模型速度方向. Figure 4 Model deformation driven by gravity (a) The horizontal velocities; (b) The vertical velocities of the model, the arrows indicate the directions of model velocities.

区域变形往往是重力和水平构造驱动等因素共同作用的结果，因此我们计算了以上两个因素共同控制下的模型变形(Model Rtg模型)，如图 5所示.在重力和水平构造挤压共同作用下，模型地壳水平运动速率由松潘—甘孜地块向四川盆地逐渐递减，地壳整体以北西-南东方向的水平挤压缩短变形为主.在地壳密度差异造成的压力梯度的调整作用下，松潘—甘孜地块至龙门山深部岩石圈地幔的水平运动速率约为3~4 mm·a^-1，明显低于上覆地壳的水平运动速率(大于6 mm·a^-1).松潘—甘孜地块与龙门山过渡区的中地壳呈现局部的水平高速流动(约10 mm·a^-1)，其主要受控于龙门山地区的陡变地形起伏.由于水平构造驱动产生的垂向变形较微弱，模型垂向变形分布与单纯重力驱动的模型垂向变形结果近似，松潘—甘孜地块近地表呈高达2.5 mm·a^-1的整体抬升，而龙门山构造带以东至四川盆地呈现区域性下沉，其中四川盆地最东边缘的下沉速率高达0.8 mm·a^-1.该模型计算得到的垂向变形与龙门山构造带两侧地表水准观测的垂向变形速率在空间分布上近似.

图 5

Figure 5

图 5 重力调整和构造驱动共同作用下的模型变形分布 (a)模型水平速度分布；(b)模型垂向速度分布；其中箭头代表模型速度方向. Figure 5 Model deformation driven by horizontal tectonic loading and gravity (a) The horizontal velocities; (b) The vertical velocities of the model, the arrows indicate the directions of model velocities.

地壳和岩石圈地幔流变性质是岩石圈变形和地质构造形成的重要控制参数(Bürgmann and Dresen, 2008; Pang et al., 2018)，因此在上述Model R模型的基础上，我们建立了2组数值对比试验来分析松潘—甘孜地块中地壳软弱层和模型岩石圈地幔流变性质对龙门山两侧垂向变形的影响，由于重力对模型垂向变形的影响较为突出，在对比模型中变形驱动力仅考虑重力效应，暂时不考虑构造加载.

为分析区域陡变地形和地壳密度结构横向差异等因素控制下松潘—甘孜地块中地壳流变性质对龙门山地区及两侧地块变形的影响，本文建立了松潘—甘孜地块中地壳不同黏滞系数的多个数值对比模型，计算结果见图 6.当松潘—甘孜地块中地壳黏滞系数较低时(10¹⁸Pa·s，Model 1，图 6a—6b)，局部陡变地形起伏驱使黏滞系数较低的中地壳形成显著的局部地壳流动，水平流速可达36 mm·a^-1；深部岩石圈地幔由四川盆地向松潘—甘孜地块呈现数mm·a^-1量级的水平流动.低黏滞系数的中地壳使松潘—甘孜地块上地壳与中、下地壳垂向运动解耦，主要表现为中、下地壳明显抬升，而上地壳垂向运动相对微弱.龙门山地区深部地壳局部下沉，稳定的四川盆地地壳相对整体下沉.随着松潘—甘孜地块中地壳黏滞系数的增加(10¹⁹Pa·s，Model 2，图 6c—6d)，地壳流空间规模及水平流速均减小，深部岩石圈地幔运动模式受中地壳流变性质的影响较小.由于中地壳黏滞系数增加，松潘—甘孜地块上、下地壳变形解耦程度减弱，岩石圈压力梯度驱使模型左侧松潘—甘孜地块地壳整体抬升，上地壳抬升速率明显小于中-下地壳.当松潘—甘孜中地壳黏滞系数增加至10²¹Pa·s (Model 3，图 6e—6f)，模型中地壳集中流动现象消失，陡变地形起伏驱使松潘—甘孜地块中、下地壳以小于2 mm·a^-1的速度向四川盆地方向运动，地壳各分层垂向变形强烈耦合呈整体抬升；龙门山至四川盆地地壳整体下沉.

图 6

Figure 6

图 6 松潘—甘孜地块中地壳流变性质对模型变形的影响 (a—b), (c—d)和(e—f)分别为重力调整作用下中地壳低密度层黏滞系数为1018, 1019, 1021 Pa·s对应模型(Model 1-3)的水平和垂向速度分布；(g—h)分别为Model 1—3和Model R地表水平和垂向速度对比. Figure 6 The influence of rheology of the middle crust in Songpan-Garzê block on model deformation (a—b), (c—d) and (e—f) are the horizontal and vertical velocities of model 1-3 with viscosities of 1018, 1019, 1021 Pa·s in the middle crust of Songpan-Garzê block; (g—h) the comparison of horizontal and vertical velocities on the surface of Model 1-3 and ModelR, respectively.

松潘—甘孜地块中地壳不同流变性质控制下区域地表变形也有所差异，如图 6(g, h)所示.在重力驱动作用下，模型地表表现为由松潘—甘孜地块向四川盆地的水平运动，垂向变形以松潘—甘孜地块的抬升和四川盆地的下沉为主.随着中地壳黏滞系数的增加，松潘—甘孜地块地表向右侧四川盆地方向的水平运动速度不断增强，垂向抬升速率也逐渐增加；中地壳黏滞系数为10¹⁹~10²⁰Pa·s的两个模型(Model 2和Model R)对应的模型地表水平和垂向运动速率相近.而中地壳黏滞系数较低的(10¹⁸Pa·s，Model 1)模型，由于松潘—甘孜地块上地壳与深部地壳的变形解耦，模型地表水平运动和垂向抬升微弱，仅在龙门山地区存在局部显著的水平运动和垂向抬升.黏滞系数较高的中地壳(10²¹Pa·s，Model 3)模型中，模型地表水平运动显著但垂向变形减弱.

地壳密度结构差异造成的压力梯度驱使岩石圈变形，不同流变性质的岩石圈地幔对压力梯度的变形响应也有所差异，图 7为岩石圈地幔不同黏滞系数对应模型的变形结果.当岩石圈地幔黏滞系数增加至5.0×10²¹Pa·s (Model 4)，深部岩石圈地幔水平运动速率约为1.0 mm·a^-1，对应松潘—甘孜地块近地表抬升速率约为0.6 mm·a^-1，四川盆地近地表下沉速率约为0.5 mm·a^-1(图 7(a—b)).当岩石圈地幔黏滞系数继续增加至1.0×10²²Pa·s(Model 5)，深部岩石圈地幔水平运动速率降低至约0.7 mm·a^-1，松潘—甘孜地块近地表抬升速率降至小于0.2 mm·a^-1，四川盆地下沉速率下降至0.2 mm·a^-1，局部地形起伏驱动的中地壳水平流速同样呈小幅降低趋势，见图 7(c—d).因此，随着岩石圈地幔黏滞系数的增加，岩石圈地幔深部由四川盆地向松潘—甘孜地块方向的水平运动速度逐渐降低，相应的松潘—甘孜地块地表抬升速率和四川盆地及龙门山地区地壳下沉速率也逐渐减小.而岩石圈地幔流变性质对局部地形起伏造成的中地壳流动的影响相对较小.

图 7

Figure 7

图 7 岩石圈地幔流变性质对模型变形的影响 (a—b)和(c—d)分别为岩石圈地幔黏滞系数为5.0×1021Pa·s (Model 4)和1.0×1022Pa·s (Model 5)对应模型水平和垂向速度分布. Figure 7 The influence of rheology of the lithospheric mantle on model deformation (a—b) and (c—d) are the horizontal and vertical velocities of models with viscosities of5.0×1021 Pa·s (Model 4) and 1.0×1022 Pa·s (Model 5) in lithospheric mantle.

陡变地形起伏造成地壳局部压力梯度，进而驱使柔性地壳变形.为定量分析地形起伏对模型变形的影响，本文建立了无地形起伏的对比模型.在不考虑地形起伏的模型中(Model T0)，模型深部岩石圈地幔同样由四川盆地向松潘—甘孜地块水平流动，流速高达4.6 mm·a^-1；以水平运动为主的深部岩石圈地幔拖曳其上覆地壳向松潘—甘孜地块方向运动，水平运动速率由深至浅逐渐降低，近地表与深部岩石圈地幔水平运动方向一致(图 8a).垂向运动上，松潘—甘孜地块地壳呈整体抬升，速率高达3 mm·a^-1；四川盆地地壳整体以高达1.0 mm·a^-1的速率下沉，垂向运动趋势呈现各地块整体的一致性(图 8b).在此模型中，松潘—甘孜地块至龙门山地区黏滞系数较低的中地壳不存在局部地壳流动.对比有、无地形起伏的两个模型地表变形分布可知(图 8(c—d))：在局部地形起伏作用下，龙门山两侧地表呈现由松潘—甘孜地块向四川盆地方向的水平运动，同时松潘—甘孜地块的抬升速率明显小于无地形起伏的模型，体现了高地形的高原在地壳范围内的“垮塌”现象；四川盆地和龙门山均呈下沉趋势，但未考虑地形起伏的模型对应的下沉速率较为显著.由此可知，龙门山地区陡变地形起伏造成地壳范围的高原“垮塌”对区域地壳局部变形存在重要的调整作用，但陡变地形的“垮塌”效应对青藏高原东缘垂向变形的影响小于区域重力的调整作用.

图 8

Figure 8

图 8 无陡变地形起伏的模型变形分布 (a—b)不考虑地形起伏的模型水平和垂向速度分布；(c—d)地形起伏模型(Model R)和无地形起伏模型(Model T0)地表水平和垂向运动速度对比. Figure 8 The deformation of model without topographic relief (a—b) The horizontal and vertical velocities of the model without topographic relief; (c—d) The comparison of horizontal and vertical velocities on the surface of models with (Model R) and without (Model T0) topographic relief.

青藏高原周缘重力观测结果显示四川盆地存在约-110~-80 mGal的布格异常，西部的松潘—甘孜地块存在约-450 mGal的布格异常(Zhang et al., 2014; Fu et al., 2014)，青藏高原东部构造运动引起的重力变化约为0.34 mGal·a^-1(Yi et al., 2016).相对接近重力均衡的四川盆地，松潘—甘孜地块具有较轻的地壳结构，存在明显的负重力异常.岩石圈结构、区域地球物理场和构造变形息息相关(宋成科等，2017)，区域地壳密度结构的横向差异通过改变静岩压力梯度进而可能驱动岩石圈变形.参考青藏高原东缘精细地壳密度结构(Zhang et al., 2014)，本文的数值模拟结果显示在地壳精细密度结构控制下龙门山两侧垂向变形显著，其中，地壳密度结构差异造成的岩石圈压力梯度引起深部岩石圈地幔由四川盆地向松潘—甘孜地块高达约4 mm·a^-1的水平对流运动，其水平运动规模大于地壳内的变形；同时，松潘—甘孜地块地壳呈现整体隆升，而四川盆地和龙门山地区整体下沉.

本文一系列数值对比模型结果显示岩石圈地幔和中地壳流变性质会影响区域重力调整和局部地形起伏控制下的模型变形.其中，在重力驱动作用下，深部岩石圈地幔黏滞系数越小，龙门山两侧地壳垂向变形差异越显著.松潘—甘孜地块中地壳流变性质通过调整上、下地壳变形耦合程度影响松潘—甘孜地块上地壳垂向运动速率，松潘—甘孜地块中地壳黏滞系数越低，地形起伏驱动形成的局部地壳流动越明显，同时也会引起松潘—甘孜地块上、下地壳显著的变形解耦造成微弱的上地壳抬升；而在中地壳相对高黏滞系数的控制下，松潘—甘孜地块上、下地壳呈现耦合变形模式，深部岩石圈地幔对流驱动地壳尺度内的整体抬升.区域陡变地形起伏驱动产生松潘—甘孜地块中地壳的局部流动及地壳尺度的局部变形，并驱使模型近地表呈现与深部岩石圈地幔方向相反的水平运动.岩石圈流变性质、陡变地形起伏和地壳密度结构差异是控制龙门山两侧地块垂向变形的重要因素.

从区域质量守恒角度分析，地表形变观测资料揭示的地壳缩短不足以产生龙门山地区如此显著的地表隆升(如，Royden et al., 1997, 2008)，因此龙门山的抬升被认为与地壳深部结构和动力学过程存在密切关联，如地壳和岩石圈地幔尺度的黏性流动驱动地形抬升(Bird, 1991; Royden et al., 1997; Burchfiel, 2004)，或延伸至地壳底部的逆冲推覆构造的逆冲运动过程促使龙门山的快速隆升(Tapponnier et al., 2001; Hubbard and Shaw, 2009).而针对青藏高原东缘现今变形模式，本文模拟计算得到龙门山两侧垂向变形趋势主要受控于地壳密度结构差异造成的区域岩石圈压力梯度，该压力梯度驱使松潘—甘孜地块地壳呈高达2 mm·a^-1的整体抬升，而四川盆地地壳呈现高达1 mm·a^-1的下沉现象，如图 9所示.此外，龙门山地区陡变地形驱使松潘—甘孜地块中地壳流动，对上覆龙门山地表垂向变形具有一定的调整作用.

图 9

Figure 9

图 9 龙门山周缘垂向运动速度观测值与计算值对比 其中，白色箭头为水准观测值，黑色为GPS观测值，灰色箭头为计算值. Figure 9 The comparison of the observed and calculated vertical velocities The white arrows are vertical velocities from precise leveling observations, the black arrows are vertical velocities from GPS observations, and gray arrows are the calculated vertical velocities.

对比龙门山两侧垂向变形观测值和计算值(见图 9)，尽管两者在变形布局上近似，但四川盆地地表垂向变形量值存在较大差异.水准观测数据显示，

龙门山断裂带以东的四川盆地地表呈现突变式的下沉，与其西侧的松潘—甘孜地块形成显著的非连续变形.而数值计算结果则表现为连续变形，四川盆地内下沉速率由西北向东南侧逐渐增加.数值计算结果与实际观测量值的差异，可能归因于数值模型中相对简化的地壳结构和连续变形的本构关系进而忽略了龙门山断裂逆冲非连续的变形滑动等构造因素.但总体上，本文根据青藏高原东缘现有的地壳密度结构模型，模拟发现地壳密度结构和地形起伏等因素对区域垂向变形的重要影响.

本文通过建立横跨龙门山断裂带的二维有限元岩石圈模型，利用牛顿黏性流体控制方程计算分析了青藏高原东缘在重力效应控制下的岩石圈变形特征，主要结论如下：

(1) 水平构造挤压驱动下，青藏高原东缘地壳呈整体微弱抬升，说明区域水平构造挤压不足以解释龙门山断裂两侧显著的垂向变形差异.

(2) 区域地壳密度结构横向差异造成的压力梯度驱使松潘—甘孜地块地壳整体抬升，速率达2.0 mm·a^-1；而四川盆地整体下沉，速率达1 mm·a^-1；因此地壳密度结构是青藏高原东缘垂向变形差异的重要驱动因素.

(3) 松潘—甘孜地块较低的中地壳黏滞系数会引起地壳内的变形解耦，导致地块上地壳抬升微弱；区域岩石圈地幔黏滞系数越低，重力调整下松潘—甘孜地块和四川盆地垂向变形差异越大.松潘—甘孜地块中地壳和区域岩石圈地幔流变性质影响重力调整作用下的模型变形量值和岩石圈变形耦合程度，是区域垂向变形的重要控制参数.