鄂尔多斯地块位于华北、华南和青藏三大亚板块的交接部位(图 1a)，地块内部断裂不发育，构造运动不活跃，史上未见有Ms＞6 的地震发生，是中国陆内著名的稳定块体，被认为是中朝地台活化后残存的古克拉通^[3~6].鄂尔多斯地块周缘为断陷盆地所围绕:西为银川盆地，东为山西盆地，北为河套盆地，南为渭河盆地.这些盆地以深大断裂为界，构造活跃，是强震活动区，历史上曾多次发生8级以上强震(图 1b).

图 1

Fig. 1

图 1 (a)中国大陆活动块体分布[1]及相对于欧亚板块稳定点的GPS速度场[2]，虚线矩形框表示数值模拟范围，实线矩形框表示鄂尔多斯地块及其邻近地区；(b)鄂尔多斯地块及其邻近地区的活动构造及地震分布情况 Fig. 1 (a) Active blocks[1] and GPS velocity tield of continental China[2] ，the dash rectangle is the whole area ofmodeling, and the solid rectangle is Ordos block and its adjacent area； (b) The earthquake distribution aroundOrdos blockand its adjacent active areas

长期以来，鄂尔多斯地块受到了广泛的研究关注.有些研究认为鄂尔多斯地块存在厚厚的岩石圈根，相当稳定^{[3， 4]}；有些研究认为鄂尔多斯地块存在较强的整体转动，争论较为激烈^[6~15].大部分研究认为鄂尔多斯地块的现今运动状态主要体现为上升运动^[6~8]、逆时针转动^{[9， 10]}、或者顺时针运动背景之上的复合叠加逆时针运动^{[11， 12]}等；少部分研究认为鄂尔多斯地块存在顺时针旋转^{[13， 14]}；另有学者提出鄂尔多斯地块相对于周围地块在作交替运动:时而顺时针转动，时而逆时针转动^[15].尽管侧重点有所不同，但所有这些研究中，鄂尔多斯地块被视为整体性很好的岩石圈块体.

多数学者认为青藏高原北东向的挤压对鄂尔多斯地块的现今构造变形状态有着重要的作用^{[7， 8， 12， 15]}.然而，鄂尔多斯地块在青藏高原的强烈推挤下为什么能够保持稳定?而稳定块体内部的运动状态又是怎样的?这些问题目前尚未得到确切的结论，值得进一步深入研究.

依据地面的地质和地球物理资料只能推测地表附近的构造活动状况和岩石圈现今结构分布，难以对三维空间的运动过程进行确切的讨论.已有很多文献报道了数值模拟技术在岩石圈动力学、运动学、地震模拟等各方面的应用^[16~21].早期的研究主要基于弹性体模型，随着流变学研究的发展，黏弹性体模型被越来越多地用于模拟构造的演化及其运动过程.然而，迄今为止的绝大多数研究都是基于概念模型，在模型几何结构、物理属性以及边界约束条件等方面都作了很大的简化，与实际情况相差较大.近几年对鄂尔多斯地块及其周边地区进行了大量的地球物理观测，得到了它们的深部结构和物质性质^[22~28].在综合整理这些资料的基础上，本文利用三维数值模拟技术来研究鄂尔多斯地块的运动和形变状态，揭示其不同圈层间的相互作用和运动特征.

三维黏弹性体模型的水平范围在100°E至124°E、30°N 至43°N 之间，东西跨度2311.1km, 南北向跨度1445.5km(见图 2).模型纵向上分为五层，分别为上、中、下地壳、地幔岩石圈和软流圈地幔，总深度300km.上、中、下地壳分界面、莫霍面、岩石圈底界面等重要间断面采用臧绍先等^[23]的研究结果.横向上依据张培震等的结果^[1]并参考GPS 速度场计算的应变率^[29]共划分出29 个块体，断裂带深度止于莫霍面.弹性参数依据PREM 模型，黏度参数依据臧绍先等^[24]给出的华北地区三维流变结构确定.

图 2

Fig. 2

图 2 三维黏弹性体模型(a)及其对应的有限元网格(b) 有限元网格模型中，纵向上共划分了14层:上地壳3层，中地壳2 层，下地壳3层，岩石圈地幔3层，软流圈3层，共计40992个单元，44880个节点，ANSYS数值模拟中选用Visco88三维黏弹性实体单元. Fig. 2 A 3D viscoelastic model (a) andcorresponding FE mesh (b) 14 layers in vertical in the FE mesh: 3 layers for upper crust, 2 layers for middle crust, 3 layers for lower crust, 3 layers for lithosphere, 3 layers for asthenosphere.There are total 40992 elements and 44880 nodes.A Visco88 solid element type was used in ANSYS modeling.

模型边界条件为:上边界应力自由；底边界垂向固定，水平方向自由；四个侧边界垂向固定，水平方向速度为全国GPS速度场^{[2， 29]}内插获得:岩石圈部分的水平速度保持不变，岩石圈底界以下，水平速度自上而下作线性衰减，至500km 处为0，这种方式下，模型底部(300km)仍有一定的水平速度.本研究只考虑水平运动，没有考虑重力的影响.

采用ANSYS软件进行三维数值计算.总时间30万年，在稳态计算的最后5 万年内，时间子步进行了加密，最后共分29个时间子步.经过多次试算，最终获得了三维速度分布，并由速度结果计算了各块体运动的欧拉矢量，以便比较其水平运动随深度的变化.

本文主要分析鄂尔多斯地块的数值模拟情况.图 3给出了数值模拟得到的地表速度与GPS 观测值的比较，可见模拟结果与观测值基本一致.表 1为数值模拟与实测的GPS速度的拟合误差分析.表 2为依据数值模拟地表速度计算的欧拉矢量与依据GPS计算的欧拉矢量的对比.在利用GPS资料反演欧拉极时，对资料进行挑选，如某台站的速度拟合残差大于所有台站拟合残差均方差的2 倍，即自动将其剔除.鄂尔多斯地块中共有GPS台站27个，经过欧拉极反演自动剔除后，最后剩下19个台站参与计算了欧拉极参数.从图 3可以看出，去掉的8个台站基本都靠近边界，形变大.对于数值模拟结果，由于介质相对均匀，欧拉极反演总体拟合残差均方差很小，块体内所有251个节点均参与了计算.

图 3

Fig. 3

图 3 稳态时数值模拟结果与GPS实测值的对比 鄂尔多斯地块内部所有观测台站比较(a)、鄂尔多斯地块内部参与计算欧拉极的观测台站的拟合对比情况(b).黑色箭头为GPS观测值，灰色箭头为计算结果.A、B、C、D 为鄂尔多斯地块内4点，在计算速度随深度的变化时需用到. Fig. 3 Comparisons of modeling results and observed GPS data Comparison for all the GPS stations inside Ordos block (a) ，comparison for the GPS stations used in calculating the Euler Poles inside Ordos block (b).The black arrows are observed GPS velocities» and the grey arrows are modeling results.A，B，C and D are 4 sites inside Ordos block.They wil be used to calculate the horizontal velocities along depth.

表 1(Table 1)

表 1 GPS台站处数值模拟结果与观测值的拟合误差分析 Table 1 Errors analysis of modeling resultsand observed data for the GPS stations

表 1 GPS台站处数值模拟结果与观测值的拟合误差分析 Table 1 Errors analysis of modeling resultsand observed data for the GPS stations

由图 3和表 1可知，数值模拟结果与GPS结果具有良好的一致性.鄂尔多斯地块内部所有GPS观测点处的平均拟合中误差与GPS 平均中误差基本一致，表明数值模拟总体拟合程度已落入观测误差以内.对于地块内部参与计算欧拉极的19 个台站，无论是东分量还是北分量，其总体拟合程度均已落入观测误差以内，这表明数值模拟结果很好地模拟了鄂尔多斯地块的刚性运动，同时也表明是局部形变或其他偶然因素导致北分量拟合稍差.数值模拟与观测结果之间的拟合差小于观测中误差，表明模拟结果已满足地球物理反演拟合的要求，这在三维动力学数值模拟中是很难达到的.此外，自动删除拟合残差较大的台站后，拟合差得到了显著的减少，说明在求解欧拉极过程中我们所采取的措施产生了显著的效果.

表 2表明数值模拟获得的鄂尔多斯地块欧拉参数与GPS资料的计算结果吻合得很好.二者的欧拉参数非常接近，数值模拟欧拉参数也完全落入了GPS欧拉参数的误差范围以内.从各自拟合的情况看，GPS观测值速度幅度的相对标准偏差为5.7%，速度方向的标准偏差为3.7°，而数值模拟结果的误差更小.无论是观测抑还是计算结果，欧拉极反演中速度幅度的拟合残差均小于20%，表明在鄂尔多斯块体内部，至少80%以上的速度值是由刚性运动形成的，亦即鄂尔多斯地块相对于欧亚稳定板块的运动中主要是刚性运动成分，故利用这些数据反演得到的欧拉极参数是可信的.

表 2(Table 2)

表 2 数值模拟结果计算的鄂尔多斯地块欧拉极与GPS速度场结果计算的对比 Table 2 Comparison of Euler poles calculated by modeling results and observed GPS data

表 2 数值模拟结果计算的鄂尔多斯地块欧拉极与GPS速度场结果计算的对比 Table 2 Comparison of Euler poles calculated by modeling results and observed GPS data

为了研究鄂尔多斯地块的水平运动随深度的变化，利用数值模拟结果，我们计算了鄂尔多斯地块不同深度的欧拉极参数.图 4 显示鄂尔多斯地块的欧拉运动在下地壳附近出现了快速的变化，而在中、上地壳和岩石圈地幔内部，欧拉极参数保持较好的稳定性.壳幔间较大的速度梯度带则表明它们之间存在较大的水平运动差异.为了进一步验证这一结果的稳定性，我们在稳态计算的基础上，在模型的西部速度边界进行了4种改变:在速度方向保持不变的情况下，速度幅度分别增大0.5cm/a和减小0.5cm/a；在速度大小保持不变的情况下，速度方向分别逆时针旋转20°和顺时针旋转20°.从图 4可以看出，在这些不同的边界条件下，鄂尔多斯岩石圈内部纵向水平运动差异均完好存在，表现出非常好的稳定性和系统性.因此，可以认为在选定的边界条件范围内，数值模拟获得的鄂尔多斯地块壳幔间的纵向水平运动差异主要由内部的物质性质所决定，并描述了其层间相互作用的主要特征.

图 4

Fig. 4

图 4 鄂尔多斯地块的欧拉极参数随深度的变化 (a)(b)(c)三个图分别表示欧拉极坐标经度、纬度以及鄂尔多斯地块绕该极点的转动角速度.图中曲线表示稳态和西边界4种突变后5万年的结果，曲线上圆圈自浅至深表示地表、上地壳底界、中地壳底界、下地壳底界(莫霍面)、岩石圈底界面. Fig. 4 Variation of Ordos block^s Euler pole along depth (a)，(b) and (c) are longitude, latitude and angle rate of Euler Pole respectively.The curves are modeling results for stable status and50000 a after the western boundary velocity disturbed for 4 modes respectively, and circles in these curves are surface and the bottom of uppercrust, middle crust, lower crust (Moho) and the bottom of lithosphere from shallow to deep of the earth.

图 5进一步表明，在鄂尔多斯地块南部，速度随深度的变化在下地壳处存在较大的梯度带，岩石圈地幔较中、上地壳运动快，速度方位角逐渐增大，向东南偏转；在鄂尔多斯地块北部，岩石圈地幔较中、上地壳的运动速度略有减小，但变化不大，下地壳处不存在明显的速度梯度带，速度方位角逐渐减小，向西北偏转.上述特征表明鄂尔多斯地块岩石圈地幔相对于中上地壳存在逆时针方向的扭动.

图 5

Fig. 5

图 5 鄂尔多斯地块内A、B、C、D 四点稳态时速度幅度(|V|)和方位(Φ)随深度的变化 图中曲线上圆圈自浅至深表示地表、上地壳底界、中地壳底界、下地壳底界(莫霍面)和岩石圈底界. Fig. 5 Variation of velocity amplitude( |V|) and azimuth (Φ)for A，B，C and D sites inside Ordos block modeling in stable status Circles in these curves are earth surface and the bottom of upper crust ，middle crust ，lower crust (Moho) and the bottom of lithosphere from shallow to deep of the earth.

本文所构建的三维黏弹性模型直接基于观测资料的计算结果.如模型横向分块直接采用张培震等的活动块体划分数据^[1]，并依据GPS计算的应变率分布情况和数值模拟结果对其做了更细致的划分；纵向分层界面则依据观测反演的三维速度结构^{[25， 26]}和热结构^[23]确定，均为空间三维曲面；三维流变结构是依据华北地区大地热流数据及流变学实验结果计算获得^[24]，计算中考虑了破裂机制的作用，得到的流变强度更为可靠，模型岩石圈部分的黏度值直接来自于该三维流变结构的插值，未作其他任何简化；岩石圈弹性参数选自公认的PREM 模型.此外，本文直接利用GPS数据内插结果作为边界条件，同时以GPS观测数据作为内部约束条件.这些措施在最大程度上减少了数值模拟过程中主观因素的影响.

在计算华北地区三维流变结构时，摩擦滑动机制假设岩石圈中断裂均匀分布，而实际情况下，华北地区存在构造活动强烈的大的活动断裂，这些断裂带上的变形要比别的地方更为强烈.数值模拟中对这些断裂带做了细致的划分，以GPS 速度场为约束，通过多次试算，对断裂带和青藏东北缘块体的黏度值进行了调整，在原黏度值的基础上乘以一个系数，使得数值模拟结果与GPS观测速度值之间得到了更好的拟合.

上述这些措施使得本文的数值模拟过程较之于其他研究更接近于实际情况，增加了数值模拟结果的可靠性，本文数值模拟得到的速度与GPS观测结果良好的一致性也支持这一结论.通过这一研究工作，本文结果印证了横向上鄂尔多斯地块内部形变小、块体保持着整体运动状态、以刚性运动为主、形变主要集中在其周缘的断陷盆地内等观点，同时还发现在纵向上，鄂尔多斯地块岩石圈中，中、上地壳与下地壳以下的岩石圈地幔之间存在较大的水平运动差异，岩石圈地幔相对于中上地壳存在着逆时针扭动趋势.由于鄂尔多斯地块内部构造简单，中上地壳横向整体性好，刚性强度大，在外力作用下，不易运动和变形，而同时周缘软弱的断裂带为其构成了外力缓冲区，在青藏高原北东向的推挤下，鄂尔多斯地块的中上地壳受到的力较岩石圈地幔弱，下地壳韧性层的存在导致中上地壳与岩石圈地幔的整体运动出现了较大的差异，从而使得岩石圈地幔相对于上地壳产生逆时针扭动.因此，我们认为稳定的鄂尔多斯岩石圈块体的水平运动在横向上具有良好的整体性，而纵向上壳幔之间存在一定的水平运动差异，由此推测岩石圈地幔与中、上地壳之间存在着较强的剪切作用力.在这个剪切力长期作用下，由于摩擦生热的原因，可能在鄂尔多斯地块下地壳处形成韧性剪切带.这个热的韧性剪切带的存在可能是鄂尔多斯地块较强的大地热流^[27]且下地壳存在显著高导层^[28]的原因.

中国地震局地质研究所王敏博士为本文提供了中国地壳运动观测网络的GPS数据，在此表示衷心感谢.同时感谢两位匿名评审人提出的宝贵意见!