华北地区中生代经历了大规模的裂陷伸展、岩石圈减薄及岩浆活动，地质构造复杂，构造活动强烈，是中国大陆东部主要地震活动区.客观认识该区现今地壳运动动力学机理无疑具有重要的现实意义.事实上，随着中国大陆及华北地区地质、地球物理及形变观测等工作的展开，相关研究逐渐增多^[1~10]. Wang等^{[1, 2]}认为太平洋板块俯冲作用对于形成中国大陆东部地区构造应力场至关重要.黄立人、陈连旺等^[3~6]研究断层运动与区域构造应力场、表面形变.冯向东等^[7]认为华北地区应力分布受到上涌地幔沿盆地边缘法向引张作用影响.何建坤等^[8]通过模拟太行山厚度转变带构造应力场分布，分析流变分层对构造应力场影响.Liu等^[9]模拟计算了长时间尺度下华北地区水平运动场及应变能积累.以上研究显示，华北地区构造运动动力学过程非常复杂，受板块运动、断层活动、岩石层弹-粘性以及重力作用等诸多因素影响.然而，不同因素对现今地壳运动及形变态势、构造应力场空间分布的具体影响程度如何？尤其是，层析成像结果显示华北地区地表以下100~400 km广泛存在低速区^[11]，浅层构造运动是否受到岩石层深部物质运移的重大影响？要深人认识上述问题，仍需更为接近实际的数值模型及模拟计算，以便进一步综合分析与研究.

华北地区于20世纪90年代引人GPS观测，地壳运动网络工程项目实施后，全区GPS流动站增至200个以上，并完成若干期复测，丰富的观测资料为研究工作奠定了基础.本文利用三维有限元模型MR (实际分层模型）与MH (均匀分层模型），结合不同边界条件及参数，对华北地区晚近时期(~Ma)构造运动进行数值模拟，将模拟结果与GPS结果、断层形变观测、构造应力场结果相比较，寻求最符合实际的计算结果，分析板块运动、周边构造块体运动、断层活动、流变分层、岩石层深部运动等对表面运动及应力分布的影响，认识其动力学机理.

本文在108°E~127°E、32°N~42°N准矩形区域建立有限元模型，厚度取100 km.为方便分析，具体构建MR与MH两个模型.MR (实际分层模型）包含区内主要活动断层，并沿垂向分7层，分别为上、中、下地壳和地幔(4层），其中地表和上、中、下地壳界面依照现有研究资料确立(详细内容参看文献[8]). MH (均匀分层模型）与MR (实际分层模型）基本类似，唯沿垂向均匀分层，其中上、中、下地壳厚度各11 km，上地幔分4层，厚度由上至下依次17、17、17、16 km.

长时间载荷作用尤其在高温压条件下岩石具有明显流变特性，满足

(1)

(2)

式中为等效应变率，σ等效应力，B为蠕变系数，R为气体常数，E为活化能，T为绝对温度，A，n为与岩石有关的常数.

本文设上地壳为弹性，中、下地壳及地幔具有蠕变性.弹性参数(杨氏模量、泊松比）由波速计算得到，中、下地壳和地幔各层蠕变系数B采用温度-深度曲线^[8]、岩石常数A及活化能E(见表 1)通过公式(2)计算得到.两者构成一组基本的材料参数(见表 2).在具体模拟计算中，须在该基本参数上实施调整，分析各类因素的具体影响.这包括：适当降低断层单元的杨氏模量，分析断层运动的影响；周整上地幔顶部单元(模型第4层单元）蠕变系数，分析地壳与上地幔顶岩石强度相对大小的影响；忽略介质流变性，视所有单元为弹性体，将岩石弹性与流变性影响进行对比分析等等.

表 1(Table 1)

表 1 岩石常数A与活化能E Table 1 Constant A and active energy E

表 1 岩石常数A与活化能E Table 1 Constant A and active energy E

表 2(Table 2)

表 2 基本材料参数 Table 2 Bascc material parameters

表 2 基本材料参数 Table 2 Bascc material parameters

计算采用有限元静力分析方法，计算软件为ANSYS.采用ANSYS软件单元库中solidl85划分单元，按照表 1、表 2设定单元材料参数.单元solidl85可实现蠕变材料(creep)计算，将蠕变系数B、常数n直接赋值即可.计算中考虑大变形导致的位移-应变非线性，取重力加速度为9.8 m/s².

如何确定边界约束对于模拟计算至关重要.我们注意到，地质资料推断的中国大陆地壳运动^[12]与地震资料计算结果^{[13, 14]}、GPS观测结果^[15]均较接近，表明地质时间尺度上中国大陆运动、构造应力场总体格局应处于基本稳定状态.其次，大陆东部地震层析成像结果显示，太平洋俯冲板片并未穿过550 km间断面，而是沿间断面向西延伸至鄂尔多斯高原下部^[11].如此大规模、长距离的板片俯冲可能引发地幔小尺度对流.数值模拟也显示，印-欧板块碰撞、青藏高原强烈隆升亦可能引发其下部物质向大陆东部运移，并成为华北地区新生代以来大规模断陷、拉张运动的主要动因^[16].由此可见，华北及周边区域岩石层下部可能存在较大范围及较大规模物质运移，并造成该区域岩石层运动在不同水平分层上的差异.

基于以上认识，本文利用GPS观测结果(采用华北及周边区域1999~2004年GPS观测结果，详情参看文献[6])给定模型四侧面及底面速度约束进行计算.第一步，选择位于模型周边稳定构造块体上的GPS测点结果，采用线性插值确定模型表面东、西、南、北四条边框上的速度约束(见图 1a).第二步，假定周边构造块体深部与浅部的水平运动方向是一致的，仅速率可能有变化，通过设定速率沿深度的线性变化确定模型四侧面速度约束.这可模拟岩石层水平差异性运动对浅表的影响.第三步，确定底面约束.

图 1

Fig. 1

图 1 模型约束示意图 （a）表面约束；（b）底面约束（第2组）；（c）侧面约束（共3组）. Fig. 1 The sketch of model condition (a) Surface condition; (b) Bottom condition (the second group); (c) Lateral face condition (three groups in total).

本文设计了3组侧面约束方式(见图 1c), 第一组, 侧面速度不随深度变化, 即设定周边构造块体呈上下一致的整体性运动.第二组, 西侧面速度不随深度变化, 东、南、北侧面速度随深度变化, 即设定岩石层水平向差异运动对东部沿海地区影响较大, 对西部山区影响较小.第三组, 模型四侧面速度均随深部发生变化, 即设定水平差异运动对整个区域产生影响.本文设计了2组底面约束方式, 第一组, 垂向固定(垂向速度为零), 水平向自由滑动, 即设定岩石层底部相对自由.第二组, 垂向固定(垂向速度为零), 水平向约束则通过侧面速度约束利用二维线性插值得到(见图 1b)，以分析地幔对岩石层底部拖曳作用的影响.

经反复试算，本文给出较典型的22个计算方案(方案编号为N01~N11，E01~E11).表 3为方案N01~N11.方案E01~E11的边界条件及参数分别与方案N01~N11的一致，但将岩石层视为弹性体.统一取加载时间4Ma的计算结果，该加载时间为中、下地壳及上地幔松弛时间的3倍以上，此时运动、形变及应力分布均已稳定，详细参看文献[17].为表述方便，模拟结果编号采用模型与计算方案编号相组合的方式，如模拟结果MR_N01为实际分层模型(MR)与计算方案N01的计算结果.共实施计算MR_N01~MR_N11、MR_E01~MR_E11、MH_ N01~MH_N11、MH_ E01~MH_E11.

表 3(Table 3)

表 3 有限单元模拟计算方案 Table 3 Simulation scheme by the finite element method

表 3 有限单元模拟计算方案 Table 3 Simulation scheme by the finite element method

华北地区1999~2004年GPS观测结果显示，全区速度方向以东南向为主，速度值由北向南逐渐增大.值得注意的是，多期GPS观测均显示速度变化具有一定非均匀性.速率等值线在区域中部发生明显弯曲，表明中部运动速率略高^[8].模拟结果MR _N01~MR_N11、MR_E01~MR_E11、MH_N01~MH_N11、MH_ E01~MH_E11均显示，模拟速度场由北而南逐渐增大，整体格局与GPS结果一致，并且除个别测点以外，模拟结果与GPS结果的方向夹角一般小于12. 5°(见图 2)，不同模拟的差别主要体现在速度值上.本文利用平均离散度D (见公式(3))衡量GPS观测结果与模拟结果之差异，结果见图 3.

图 2

Fig. 2

图 2 模拟计算MR_N09结果与GPS观测结果对比 空心箭头为计算结果，实心箭头为观测结果. Fig. 2 The contrast between the results of MR_N09 simulation and the GPS data Horrow narrows denote simulation results, solid narrows denote GPS data.

图 3

Fig. 3

图 3 各模拟结果的平均离散度 Fig. 3 Average dispersion of simulations

(3)

式中，D为模拟结果与GPS结果的平均离散度；m为区域内GPS测点总数；V_hci为第i个GPS测点处的模拟速率，V_hoi为该测点的GPS观测结果.

根据图 3我们可知：①考虑岩石层流变特性的模拟结果(MR_N01~MR_N11，MH_N01~MH_ N11)与观测结果较接近(D值分别为1. 153~1. 666 mm²/a，1. 145~1. 577 mm²/a)，仅将岩石层视为弹性体的模拟结果(MR_E01~MR_E11，MH_ E01~MH_E11)则与观测结果存在较大差异(D值分别为1. 172~4. 430 mm²/a，1. 19~4. 461 mm²/a).②若计算方案相同，实际分层模型结果(MR_N01~MR_ N11，MR_E01~MR_E11) -般优于相应的均匀分层模型结果(MH_N01~MH_N11，MH_ E01~MH _E11).③考虑断层运动的模拟结果(MR_N02、ME_ E02)优于不考虑断层运动的模拟结果(MR_N01、ME _E01)，前者D值分别为1 1881 181 mm²/a，后者D值分别为1. 189、1. 182 mm²/a.④若改变蠕变系数B，使地幔顶部岩石强度等于或高于下地壳的岩石强度(计算方案N04、N05)，则模拟结果(MR_N04~MR_N05, MH_N04~MH_N05)略差于相同约束条件(计算方案N03)的模拟结果(MR_N03、MH_ N03)，前者D值范围分别为1 179~1 181、1. 180~1.182 mm²/a，后者D值分别为1.176、1.175 mm²/a.⑤对于模型相同的计算，若考虑下部较上部稍快的岩石层水平差异性运动(计算方案N03~N09，E03~E09)则模拟结果与实测结果更接近.在所有模拟计算中，MR_N09结果与实测结果最接近，D值为1.145 mm²/a.该计算采用MR模型(实际分层模型），考虑断层运动(区分断层单元、非断层单元）、地形重力以及岩石流变性，边界约束不仅考虑模型四侧面的水平向差异运动(100 km深度处的速度值为表面速度的1. 2倍），并且考虑模型底面的整体移动.

综合以上结果，研究区表面水平运动与周边构造块体及深部岩石层的相对运动方式密切相关.如果仅将周边块体的水平运动视为“上下一致”的整体性运动，所得结果与实测结果的差异往往较大，然而，若视周边块体运动为深部快、浅部慢的差异性整体运动，并将区域下部物质的整体性水平移动也纳入考虑，则可有效提高模拟结果与实测结果的吻合度.此外，断层运动、重力作用、岩石层流变分层对研究区表面运动也有一定影响，将其纳入考虑，可使模拟结果得到一定程度改善，但与前述水平差异性运动方式相比，当属次要因素.

为观察断层运动性质及水平速率，本文在模型表面设定若干条测线(见图 4a)分别横跨张家口-渤海构造带、山西断裂带、唐山-磁县断裂带、郯庐断裂带.得到测线各点沿断层方向的速度分量，设张渤带自W向E为断层正向，山西带、唐-磁带和郯庐带自S向N为断层正向.图 4(b~f)为模拟计算MR_N09的结果.

图 4

Fig. 4

图 4 断层水平运动模拟结果 (a)在模型表面设定的跨断层测线AA′、BB′、DD′、EE′、FF′; (b~f)模拟计算MR_N09各测线沿断层方向的速度分布. Fig. 4 Fault horizontal velocities resulted from finite element simulation (a) The locations of the survey lines across the faults on the model surface which are named AA′, BB′, DD′, EE′and FF′, respectively; (b~f) The distributes of the velocity component along the trend of the faults of each survey line.

采用计算方案N02~N11、E02~E11，断层单元的弹性模量为周围介质的十分之一，根据模拟结果可知：①断裂带两侧水平速度发生较大变化，且在断裂带上变化最为剧烈，断层宽度约5~6 km，揭示出断裂带鲜明的扭动特性.各模拟结果均显示张家口-渤海断裂带为左旋，山西断裂带、唐山-磁县断裂带和郯庐断裂带为右旋，与地质研究及实际观测一致^[18].②各模拟计算得到的断层水平错动速率均较低，为0. 01~0.135 mm/a.③若计算方案相同，MR模型计算(实际分层模型）显示的断层速率高于MH模型结果(均匀分层模型).如模拟计算MR_N09显示张-渤带、山西带、唐-磁带、郯庐带速度依次为0. 130、0. 135、0. 098、0. 070 mm/a，而模拟计算MH_ N09依次为0. 120、0. 121、0. 078、0. 065 mm/a.另外，若模型相同，考虑介质流变性(计算方案N02~N11)的断层速度高于完全弹性的计算结果(计算方案E02~E11).④断层速度明显受垂向约束的影响.若区域周边为上慢下快的水平向差异运动(计算方案N03~N11、H03~H11)则断层速度较快.如模拟计算MR_N02显示张-渤带、山西带、唐-磁带、郯庐带速度依次为0. 083、0. 104、0. 073、0. 080 mm/a，明显低于模拟计算MR_N09、MH_N09结果，后者考虑了岩石层差异运动.④地质资料认为华北平原各盆地断层活动的平均速率一般小于1 mm/a，有些则在0. 1 mm/a以下.根据断层形变观测^[18]，最近20年张-渤带断层水平速率为0. 13~0.14 mm/a，山西带为0.14~0.15 mm/a, 唐-磁带约为0.12 mm/a, 郯庐带为0.07~0.08 mm/a.与其他模拟结果相比，模拟计算MR_N09与断层形变观测最为接近.并且，该模拟计算显示山西带断层速率最高，郯庐带最低，亦与实际观测相符.

以上表明，研究区呈共轭分布的活动断裂带的现今活动特性，EW向断层左旋运动(张家口一渤海断裂带）、NE向断层右旋运动(山西断裂带、唐山一磁县断裂带和郯庐断裂带），明显受控于周边构造块体相对运动的大环境，并具有继承性.同时，研究区重力作用和流变分层的横向非均匀分布导致断层运动速度普遍增大.此外，研究区周边及下部的水平差异性运动同样会提高断层运动速度.

已有研究表明^[19]，华北地区主张、主压应力大致水平，方向分别为NNW-NS、NEE-EW，显示受统一的应力场作用.此外，某些局部区域主压应力接近直立，与整个区域存在差异，如山西带南部.本文结果显示，①各模拟结果得到的约为10 km深度处主压、主张应力方向均很接近，且全区一致性较好，主张应力接近水平，方向NNW-NS，主压应力在水平面的投影为NEE-EW向.以上与现有结果相吻合，表明研究区构造应力场的形成，主要为周边构造块体相对运动和相互作用的结果.②断层对应力场主方向的影响局限于断层附近及其内部.另外，若采用相同的计算方案，模型MR结果(MR_N01~MR_ N11)与模型MH结果(MH_H01~MH_H11)给出的主应力方向差别很小.表明断层分布、地形重力以及岩石分层结构的横向非均匀性对研究区构造应力场方向的影响较小.

需要指出的是，模拟计算得到的主压应力仰角在全区是变化的，北部地区仰角较高，南部地区则接近水平，该仰角变化与现有研究结果存在一定差别. 图 5a为模拟计算MR_N09在约10 km深度应力状态和震源机制解结果(其中，红白球为模拟结果，兰白球为哈佛大学震源机制解，白色区中轴为主压轴）.本文认为，其一，由于模拟计算边界约束源于GPS观测结果，可能在一定程度上对研究区水平向构造挤压作用(如太平洋板块俯冲作用)估计不够充分，从而影响了计算结果.其二，模拟计算揭示，造成山西带南部主压应力高仰角的原因是复杂的，除地形重力以及岩石层流变特征等因素以外，应该存在与局部构造运动有关的动力学过程，且后者的影响更为重要.

图 5

Fig. 5

图 5 应力场模拟结果 (a)模拟计算MR_N09~10 km处应力场与震源机制解对比(红白球为模拟结果，兰白球为震源机制解，主压轴位于白区中轴）(b)模拟计算MR_N01~10 km处最大剪应力分布(单位MPa); c)模拟计算MR_N09~10 km处最大剪应力分布(单位MPa). Fig. 5 The stress field resulted from finite element simulation (a) The stress field in~10 km depth resulted from MR_ N09calculationcompared with focal mechanism (Red-white ball denotes simulation result, blue-white ball denotes focal mechanism. The compress principal stress axis is located in the middle of white region, respectively); (b) The maximum shear stress distribution in~10 km depth resulted fromMR _ N01 calculation (unt, MPa.); (c) The maximum shear stress distribution in~10 km depth resulted from MR_N09calculation (unit, MPa).

各模拟结果显示的最大剪应力分布存在较大差异.①断层对最大剪应力分布影响非常明显.图 5(b、c)为模拟计算MR_N01、MR_N09给出的约10 km深度处最大剪应力等值线图.可以看出等值线高梯度带的走向与断层分布密切相关.实际上，若不考虑断层运动，其等值线在全区平缓变化，无法形成高梯度带.②垂向约束对最大剪应力分布的影响亦很显著.若不考虑岩石层水平向差异运动(计算方案N02, 水平速度约束不随深度变化），则高梯度带主要沿NNE向分布，且南部地区剪应力值、高梯度带的等值线密集度高于北部地区，见图 5a.若考虑差异运动(计算方案N03~N11，水平速度约束随深度逐渐增大），则剪切高梯度带沿NNE、NWW呈共轭分布.北部地区剪应力值增大，其中张家口一渤海断裂带、郯庐断裂带北段附近尤为突出.现有研究显示上述断裂带现今活动性较强，先后发生1975年M_s7. 3海城地震、1976年M_s7. 8唐山地震.由此可见，华北地区岩石层水平差异运动，也即岩石层下部对表层的拖动，对于形成NE-NNE、NWW共轭高剪切带体系、促进内部活动地块(太行山、冀鲁、豫淮次级活动地块、鲁东-黄海活动地块）间的相对运动起到不可忽视的作用.③若计算方案相同，模型MR结果(MR_N01~MR_N11)与模型MH结果(MH_ H01~MH_H11)差异不大，表明地形重力、流变分层的横向非均匀性影响是次要的.

本文在三维有限元MR (实际分层模型）、MH (均匀分层模型）框架下，采用不同侧面、底面约束组合以及岩石弹性、流变性设置对华北地区晚近时期(~Ma)构造运动进行了数值模拟，并将计算结果与GPS观测、跨断层观测以及应力场研究结果等进行比对，分析诸如块体运动、断层、地形重力、岩石流变结构等因素对该区现今运动及构造应力场的具体影响.综合模拟结果，对于该区现今地壳运动动力学机理的初步认识概括如下：

(1)华北地区现今地壳运动受控于周边大型构造块体(鄂尔多斯地块、华南活动地块、东北亚活动地块）相对运动的基本态势，这表现在以周边区域GPS观测结果作为基本约束条件所给出的表面速度、应力场模拟结果的总体特征与现有结果是比较接近的.一方面，模拟速度场与GPS结果的分布格局一致，速度值的平均离散度也处于较小范围(在1.145~4.43 mm²/a).另一方面，模拟应力场全区一致性较好，主张应力接近水平，方向NNW-NS，主压应力为NEE-EW向，与现有研究结果一致.

(2)太平洋板块俯冲、印度-欧亚板块碰撞挤压可能引发小尺度地幔对流等动力学过程，造成华北及周边地区岩石层深部、浅部水平运动具有差异性，深部物质沿SEE-SE向移动速度相对于浅部较快，并对浅表起拖动作用，从而对整个区域构造运动、应力分布以及内部次级块体相对运动产生深刻影响.本文认为，广泛存在于华北及周边地区的差异性水平运动是导致本区表层速度在空间上非均匀变化的主要因素~并对剪应力高梯度带分布产生较大影响. ①模拟计算显示，若区域周边以深部、浅部完全一致的方式进行整体性水平运动，模拟结果与实测的差异往往较大~若以深部稍快、浅部稍慢的方式运动，则模拟结果得到明显改善，尤其考虑了岩石圈底部的拖动作用后，效果会更好~最佳模拟结果为MR_ N9, 该计算设定100 km深度处的速度值为表面速度的1.2倍，并且考虑岩石层底部的整体性移动.②最大剪应力分布的模拟结果显示，深部速度稍快于浅部，可显著提高区域北部的最大剪应力值，并沿张家口-渤海断裂带、郯庐带北段形成高值区及高梯度带，从而使全区形成NE-NNE、NWW共轭高剪切带体系，对内部次级活动地块(太行山、冀鲁、豫淮次级活动地块、鲁东-黄海活动地块）间相对运动起到有效促进作用.③根据本文计算结果，岩石层水平差异运动可以提高断层运动速率，使其与断层形变观测结果更为一致.

(3)断层分布对于区域表层运动及应力格局的影响各具特点.一方面，断层分布对表面运动的影响一般集中于断层两侧较小范围内(图 4显示，部分跨断层速率变化在距离断层稍远区域很快恢复，另一部分虽未恢复，但变化量很小，即断层运动速率约0.1 mm/a), 故对表层运动总格局的影响是有限的.断层对运动场的影响主要表现在促使区域内次级活动块体周边产生微小形变，使表层速度空间分布更具非均匀性，从而更接近实际情况.考虑断层运动的模拟结果与GPS观测更相符即说明此点.但总体上，与上文第一点、第二点相比，断层分布对表层速度空间变化的影响属于较次要地位.另一方面，断层分布对于区域应力场影响显著，主要表现为对剪切应力的影响.模拟计算显示，沿断裂带容易形成剪应力高值带及高梯度带，促使应变能在局部快速积累.尤其是区内呈共轭分布的高剪切带体系的存在，对于促使区域内次级构造块体间相对运动，以及地震孕育、发生起到非常重要的作用.

(4)地形重力作用及岩石流变结构的横向非均匀性的影响主要表现于促使区域内断层运动速度有所增强.但总体而言，上述因素对于华北地区现今构造运动的影响应处于较为次要的地位.计算显示，采用相同的计算方案，与模型MH结果(均匀分层模型）相比，模型MR (实际分层模型）计算对模拟结果的改善仍是有限的.另外，本文尝试改变上地幔与下地壳的相对强度(计算方案N04、N05)，结果显示在强下地壳、弱上地幔情形下，模拟结果与现有GPS观测等结果更接近，反之稍差.表明对于华北地区而言，下地壳岩石强度有可能高于上地幔顶的岩石强度.

根据本文结果，华北地区现今表面运动及应力场除受周边构造块体(包括鄂尔多斯地块、华南地块、东北亚地块等）相对运动的重大影响外，岩石层深部运动的影响是不可忽视的.模拟结果显示，区域下部较快速的“拖动”作用可以使表面速度场得到更好模拟，并在地表以下约10 km深度处形成共轭分布的剪应力梯度带.可见太平洋板块的俯冲作用、印-欧板块的碰撞挤压作用等可能造成区域深部、浅部水平运动差异，从而对研究区现今地壳运动产生深刻影响.模拟结果同时显示，地形重力作用、断层分布及岩石层流变分层的非均匀性也具有一定影响作用，但应处于较为次要地位.

本文利用数值模拟对华北地区现今运动及形变的动力学机理进行分析与研究，在建模、设定计算方案、参数设置及与实际观测对比等环节，笔者尽量做到与现有研究成果相结合.尽管如此，仍存在需改进的地方，如模拟计算假定华北地区晚近时期周边构造块体、岩石层下部物质的运动状态保持恒定不变，忽略了该运动状态在短时间尺度上的波动，这与实际情况可能存在差别.又如，受到现有研究结果的限制，在模型构建上无法做到与现今华北地区的实际岩石层流变分层、断层分布完全一致.另外，本文没有考虑热过程的影响~这些因素均会导致模拟计算本身存在一定偏差.此外，受条件限制，现有GPS结果、跨断层形变观测及构造应力场研究尚无法非常充分地反映华北地区地壳运动、应力场全貌，使得模拟结果与实际情况的对比不能更为精细.相信随着各方面资料的积累，这项工作会取得更好的效果.