作为欧亚板块和印度一澳大利亚板块之间的主要边界^[1~3], 巽他弧的构造演化历史^{[4, 5]}以及与俯冲板块的俯冲形态、俯冲角度、俯冲方向密切相关的一系列构造现象, 例如緬甸微板块的拆离^[6]、安德曼海盆的打开^{[7, 8]}、平行海沟方向的断裂滑移和地震的发生等, 一直都为该区研究热点.尤其在2004年巽他弧中部苏门答腊地区发生了M_w9.3的强烈地震之后, 对于该区发震机制的解释、贝尼奥夫带构造特征的认识以及应力应变场的分布特征等得到了更加广泛的关注.例如, 朱守彪等9利用有限元方法模拟出在闭锁状态下苏门答腊地区地震孕育、发生、迁移的一种可能过程.Grevemeyer等^[10]根据巽他弧各段的布格重力异常提出俯冲带耦合面积的大小决定了地震的强弱及其规模.ChHeh^[11]利用大地测量数据和数值模拟方法描绘出苏门答腊地区俯冲带各段的耦合强度, 并以此为基础讨论了耦合强度与地震发生的空间分布关系.而Shapro^[12]则提出巽他海沟强耦合带的形成以及地震的发生都与俯冲到巽他板块之下的Wharton洋中脊^{[13, 14]}有关, 同时还指出这种强耦合带控制了北西向断层和盆地的发育.

由以上前人结果可知, 无论采用何种方法对苏门答腊及其周边地区的地震机制、构造特征进行分析时, 都不应忽略俯冲带耦合作用的影响.而这种影响是如何控制弧内的应变分离^{[15, 16]}、应变强度分段、主应变方向上的偏转以及地震破裂带的空间分布特征的呢~虽然前人有所论述, 但并无直接证据.本文将通过采用有限元数值模拟的方法来直观地验证这种耦合作用的影响.

东南亚位于欧亚大陆东南部, 主要由巽他次级板块、緬甸微板块等拼合而成, 其边界多为碰撞和俯冲带(图 1).由于受到欧亚板块、印度一澳大利亚板块以及菲律宾板块的联合挤压作用, 该区构造活动异常强烈, 多伴有强烈地震的发生.尤其在其西部边界, 地震活动更加明显, 近200年来已知7级以上的地震就达十几次.而作为地震孕育、发生的直接力源, 印度一澳大利亚板块向巽他板块之下不断地斜向俯冲^{[17, 18]}, 俯冲推移距离自南向北有所变化, 但平均约为60mm/a.

图 1

Fig. 1

图 1 东南亚地区大地构造立体模式图 图中虚线框表示模型研究区. Fig. 1 3D tectonic model of the Southeast Asia Rectangle with dash line s the range of study area.

东南亚西边界的北部为安达曼海盆(图 1).海盆底部发育有北东向的雁列式扩张脊, 该扩张脊被一系列的北北西和近南北向的转换断层所切割^[19].对于其扩张机制, 一些学者认为它是受到苏门答腊断裂和实皆断裂活动影响的拉分盆地^[20], 而另一些学者认为该海盆的形成主要受到印度板块斜向俯冲作用的影响^[21].总之, 安达曼海盆并不是一个典型的弧后盆地, 因为其扩张方向并没有垂直于海沟方向, 这明显不同于其他活动的弧后盆地^[22].安达曼海盆平均的扩张速率约为18mm/a^[23], 且该区地震类型多为拉张型和走滑型.

东南亚西边界的东南部为苏门答腊(图 1), 从北西至南东可延伸约1600 km.根据板块构造模型可知, 印度一澳大利亚板块的北东向运动在该区可被分解成两部分位移.一部分平行于海沟方向, 且能够被右行走滑的苏门答腊断裂和明打威(Mentawei)断裂调节吸收^[24].另一部分向下俯冲至巽他次级板块的地幔中.在俯冲过程中, 由于材料属性、摩擦系数以及温度上的变化能够导致上驮板块和下俯板块之间产生耦合接触, 这就为应变的积累、地震的发生提供了有利的环境.苏门答腊断裂是该区内一条非常重要的断裂, 其发育于火山弧上, 调节吸收了俯冲带内北西向的走滑运动.沿苏门答腊断裂的右旋位移量由南向北逐渐增大, 在最南端滑移率仅为6mm/a.

本文有限元模型(图 2)是基于沟-弧系统的概念理论模型, 研究区北起緬甸中部的莫则沃(约20°N), 南至东爪哇的Sumba (约10°S), 西部边界沿巽他弧向西扩展了约400 km, 而东部边界则可到达巽他板块中部地区的中途岛附近(约110°E).模型区内主要包括大陆岩石圈(巽他板块)、上浮于上地幔之上的大洋岩石圈(印度一澳大利亚板块)、俯冲至上地幔300 km处的俯冲板片、软流圈上地幔以及位于俯冲板块之上的低黏度地幔楔^[25].模型总厚度700 km, 其中大洋岩石圈厚度为80 km, 大陆岩石圈地壳厚度为35 km, 岩石圈地幔厚度为65 km, 软流圈上地幔厚度约为600 km.据现今的层析成像技术以及布格重力异常可知, 沿整个巽他弧俯冲带的俯冲迹线和俯冲角度是不断变化的, 而这种几何形态上的变化必定会影响弧内的应力-应变环境.因此, 为了能够真实地反映这种变化的影响, 本文根据前人的研究成果^{[12, 26, 27]}拟合出俯冲板片近真实的地质几何学形态.如图 2b所示, 从Sumba至苏门答腊北部的俯冲带俯冲角度变化不大, 约为30°而从苏门答腊北部至緬甸弧中部(横跨整个安达曼弧)的俯冲角度则在30°和53°之间变化.

图 2

Fig. 2

图 2 有限元模型及其边界条件 图中箭头表示运动方向；三角形表示固定方向. Fig. 2 Finite element model and is boundary conditions Arrow s motion sense; Triangle is a fixed direction.

模型中考虑了3条主要断层, 分别为俯冲带巨型逆冲断层、苏门答腊断层以及实皆断层(如图 2a所示).其中实皆断层自新生代以来滑移距离大约430 km^{[28, 29]}, 苏门答腊断层北部自中第三系以来滑移距离大约130 km, 而南部滑移距离大约35km^{[30, 31]}, 这些断层的周边地区经常伴有地震发生, 且多为走滑型地震, 断面近直立, 因此在模型中将它们的倾角设为90°以方便计算.与上述两条断层相比, 俯冲带巨型逆冲断层的形成是与俯冲板块的空间形态直接相关, 因此在模型中该条断层的倾向、倾角及走向伴随俯冲板块空间形态变化而变化(图 2b).

对于上述3条断层的接触行为的判定, 本文使用库仑摩擦模型^{[32, 33]}和耦合法来描述.库仑摩擦模型的表达式如下:

(1)

其中p表示断层接触面间的正应力, μ表示断层接触面间的摩擦系数.在模拟过程中, 断层接触面受力方式可分解为法向力和切向力, 其中法向力表现为垂直于断层面的压力, 而切向力与库仑摩擦模型(公式(1))剪应力的对比则描述断层面的接触状态.即, 当切向力小于剪应力τ时, 断面处于粘合状态; 当切向力大于剪应力τ时, 断面的接触状态由粘合转为滑动, 滑动距离的最大值设为每一迭代步中断面接触单元平均长度的1%.而对于断面接触力的计算, 本文采用了有限元接触算法中的惩罚函数^[9]法, 其法向力的表达式为

(2)

其中n表示接触面法线方向, f_n表示法向接触力, k_n表示法向接触刚度, h_n表示穿透距离.切向力(摩擦力)的表达式为

(3)

其中t表示接触面切向方向, f_t表示切向接触力, k_t表示切向接触刚度, h_t表示滑移距离.模型中耦合法用来描述俯冲带耦合区的接触行为, 在本实验中使用了全耦合法.即, 断层上盘接触面耦合区节点的自由度完全耦合于断层下盘接触面耦合区节点.因此在模拟过程中, 上覆板块耦合区节点的运动方向会伴随俯冲板块运动方向的改变而改变.

对于上述两种接触方法的使用, 本文是这样给定的.首先, 将断层上下盘接触面区分为普通接触区和全耦合区.然后, 在普通接触区将断面接触力设置为惩罚函数法, 并以库仑摩擦模型为准则判断断面接触行为, 而在耦合区设置为全耦合法.最后, 通过对比模拟结果调整普通接触区和全耦合区面积的大小.由于本模型将各板块视为以千年为尺度运动的准静态模型, 因此忽略重力、地震可能对断层接触面及其周边地区的影响.

岩石和矿物的变形实验研究和大量的野外观察证明, 中-深度地壳和上地幔的大多数岩石变形, 都是通过组成岩石的矿物的塑性变形和流动变形完成的^[34], 而在微观尺度上主要变现为位错蠕变和扩散蠕变两种过程.为了能够正确描述这种塑性流变学行为, 本文采用了稳态幂指数流变速率公式^{[35, 36]}, 其表达式如下:

(4)

式中表示应变率, A表示物质结构常数, R表示气体常数Q表示活化能, σ表示差应力, n表示应力指数.在模拟过程中, 以上参数的选择主要是基于前人的研究成果^[37~39](表 1).

表 1(Table 1)

表 1 模型参数 Table 1 Model parameters

表 1 模型参数 Table 1 Model parameters

而对于公式(4)中温度T, 在忽略重力及温度演化的情况下, 该值的高低将决定岩石圈和上地幔横向和纵向上的不均匀性, 而这种不均匀性又表现为不同区域间不同的流变学行为.为了能够描述这种温度上的变化, 本文首先通过有限元方法(模型底面温度设为1900℃和顶面温度设为0℃模拟出岩石圈和上地幔的温度场(图 3), 然后将这一结果作为初始温度代入公式(4)中求解稳态幂指数流变速率模型.初始温度场模型参数的选择如表 2所示.

图 3

Fig. 3

图 3 初始温度场 Fig. 3 The initial temperature field

表 2(Table 2)

表 2 初始温度场模型参数 Table 2 Model parameters of the initial temperature field

表 2 初始温度场模型参数 Table 2 Model parameters of the initial temperature field

对于有限元模型的边界条件本文是这样给定的:为了不让外载荷直接作用于研究区内, 本模型将巽他弧以西400 km内的大洋板块作为模型的西部边界, 并根据GPS^[40]数据的计算结果, 将运动速率设为由南部的65mm/a逐渐递减到北部的55mm/a, 由于本模型没有考虑初始应力及形变的影响, 因此将模拟时间设置为6000年, 总位移量在330~390 m之间, 计算步数为1000, 而运动方向则设为由南部的北东向逐渐偏转为北部的北北东向.模型的东部边界位于巽他板块中部地区的中途岛附近, 根据前人模拟速度场^[41]可知, 该区的运动方向向东, 几乎未受到菲律宾板块俯冲作用的影响, 因此将其设为自由边界.模型的北部边界位于20°N附近, 由于受到青藏地块和华南地块的挤压阻挡作用, 因此将该边界设为固定边界.而对于模型的南部边界, 虽然受到澳大利亚板块正北向俯冲作用的影响, 但通过GPS测量结果可以看到, 在该作用下地表却表现为向东的运动, 因此在模型中将该边界设置为具有一定量的运动速率, 模型的上层表面本文将其设为自由边界, 模型下层底面垂向本文将其设为固定, 而水平方向设为自由界面.

俯冲带的耦合作用对苏门答腊地区地震的发生、孕育以及弧内一系列构造现象的演化起着重要的控制作用, 而这种耦合作用的强弱又是与耦合面积的大小直接相关的.Grevemeyer^[10]曾指出苏门答腊俯冲带地震规模明显高于Java俯冲带的原因正是由于耦合面积大小决定的.虽然这种观点已经被广泛接受, 但在苏门答腊地区2004年M_w9.3地震上的应用却存在一些争议, Shapiro^[12]认为强震的产生源于高温、低密度的俯冲板块与上驮板块之间的充分耦合.而McCafrey^[17]和Chlieh^[11]则分别利用GPS数据得出苏门答腊南部应处于弱耦合带, 但对为什么会在弱耦合带发生如此强烈的地震没有做更多的解释.由此可见, 在苏门答腊地区对耦合面积和区域应力、应变场之间关系的认识是十分必要的.为此, 本文根据前人的研究成果测试了俯冲带五种不同的耦合面积, 如图 4所示.其中第一和第二种耦合面积分别表示在苏门答腊2004年地震主震区和整个苏门答腊俯冲带处于弱耦合状态下的耦合模型, 这样设置的目的是通过两者的对比查看俯冲带有无耦合对苏门答腊前弧的影响.第三种耦合模型的设置是基于Shapio^[12]的假设, 即苏门答腊2004年地震主震区俯冲带处于一种高强度耦合状态~第四和第五种耦合模型则是基于Chlieh^[11]和McCafrey^[17]的测量结果, 即苏门答腊俯冲带的中南部处于强耦合状态, 而2004年地震主震区则处于弱耦合或无耦合状态.

图 4

Fig. 4

图 4 俯冲带内五种不同的耦合面积模型 图中白色框表示耦合面积，黑色虚线框表示摩擦区a、b、c、d、e分别对应于五种测试模型，A、B、C、D分别代表緬甸弧、安达曼弧、苏门答腊2004年地震的主震区以及苏门答腊南部带. Fig. 4 Five kinds of models with different coupled area In picture, rectangle with white line is a coupled area. Rectangle with black dash line is a friction zone. a, b, c, d, e correspond to models (1, 2, 3, 4, 5), respectively. A, B, C, D correspond to Burma-Arc, Andaman-Arc, The northern Sumatra, The southern Sumatra, respectively.

在印度一澳大利亚板块斜向俯冲作用下, 不同耦合面积对上驮板块应变场分布特征的影响显然不同.有些区域可能应变过于集中从而产生破裂或者屈服, 而有些区域可能应变相对较弱从而保持完整.为了能够正确描绘出这一应变分布特征, 本文采用了应变强度这一概念, 其表达形式如下:

(5)

式中ε₁表示第一主应变, ε₂表示第二主应变, ε₃表示第三主应变图 5即为五种不同耦合面积对研究区应变强度影响的等值线图.

图 5

Fig. 5

图 5 五种模型的应变强度分布图 图 5a中，黑色虚线框表示表面应变场主应变方向的提取范围，黑色实线表示表面应变强度的测线位置；图 5e中，黑色虚线表示垂直位移场剖面位置. Fig. 5 The distributions of strain intensity of five models with different coupled area In picture 5a, rectangle with black dash line is the range of the direction of principal strain field, black lines is the location of measure lines of surface strain intensity; In picture 5e, black dash line is the location of profile of the vertical displacement field.

从应变强度分布图中(图 5)可以看到, 在苏门答腊南部带(介于测线aa'和dd'之间)印度一澳大利亚板块的俯冲方向与苏门答腊海沟呈45°俯冲速度约为65 mm/a, 在这种斜向俯冲作用下该区应变强度值的分段性特征非常明显.应变强度首先集中于苏门答腊南部带最南端(介于测线aa'和bb'之间), 强度值伴随耦合面积的增大不断增高(图 6a).从测线aa'(图 6a)和测线bb'(图 6b)、cc'图 6c)的对比可以看到, 伴随南部耦合面积的增加能够导致中部地区(介于测线bb'和cc'之间)靠近海沟一侧应变强度值的降低, 从而形成一个应变低值区.经过应变强度较弱的中部, 北部地区(介于测线cc'和dd'之间)的应变强度值在强耦和作用下能够再次升高, 但最大值明显小于南部.这说明南部强耦合带对应变积累的影响最远能扩展至中部地区, 而北部的应变强度只取决于该区耦合面积的大小(图 6d、图 6e), 但由于俯冲方向和俯冲速度上的变化, 该区应变强度值小于南部地区.由此我们推测苏门答腊南部带地震破裂模式以及构造演化特征很可能受到区域间耦合面积相互作用的影响, 而呈现出分段性.即, 在印度一澳大利亚板块斜向俯冲过程中, 应变可集中于俯冲带内某强耦合接触界面(对于本文模型该区位于苏门答腊南部带最南端), 该界面能够使上驮板块产生一个垂直海沟方向的挤压应变以及平行海沟方向的剪切应变^[16].该种应变分离模式控制了弧内一系列逆冲断层、走滑断层以及拉分盆地的生成^[42].当应变积累达到破裂极限时, 便可导致地震的发生, 并在震源周边地区产生一定长度的破裂带(对于本文模型该破裂带可能位于南部带中部地区).伴随该破裂带能量的释放应变可重新集中, 但集中区应该只能位于破裂带的两侧(对于本文模型新的集中区可能位于南部带北部地区).同样, 当新的应变集中区能量积累到一定程度时, 便可导致地震的发生, 从而完成一次地震破裂的迁移.

图 6

Fig. 6

图 6 不同模型间测线的表面应变强度曲线 在图(a~g)中，x轴对应图 5a中的7条测线位置, y轴表示应变强度值.每幅图中五种颜色的应变强度曲线分别代表五种模型. Fig. 6 The variation of surface strain intensity in five models with difference measure lines In pictures (a~g), x-axis correspond to the position ofmeasure line, y-axis refers to value of strain intensity.Five lines with difference colors in those picturs correspond to models (1, 2, 3, 4, 5), respectively.

在苏门答腊2004年地震主震区, 俯冲板块的俯冲方向与海沟夹角约为40°俯冲速度略减为60mm/a, 在这种条件作用下该区应变场的分布特征较为特殊.从测线dd'(图 6d)和测线ee'(图 6e)的对比以及图 7中可以看到, 伴随苏门答腊南部带耦合面积的变化, 该区应变场的主应变方向能够发生改变, 而应变强度值却变化不大.但伴随2004年地震主震区耦合面积的变化, 该区应变场的主应变方向和应变强度值都发生了明显改变.这说明2004年地震主震区的应变场分布特征主要受控于该区的耦合强度, 但同时不能忽略苏门答腊南部带对其主应变方向带来的影响.在安达曼弧地区, 俯冲板块的俯冲方向与海沟夹角从40°逐渐转变为近平行展布, 而俯冲速度也从60mm/a逐渐递减为55mm/a, 从图 6中可以看到在此种条件作用下不同耦合面积的应变强度差自南向北逐渐减小, 说明耦合面积的影响范围只能扩展至安达曼弧的中部地区.而从图 5中可以看到实皆断裂和苏门答腊中央断裂对研究区主应变方向的分布特征起到重要的控制作用.

图 7

Fig. 7

图 7 五种模型的主应变方向差异对照图 A、B分别对应图 5a中A、B黑色虚线框↔.表示应变拉伸方向，>- < 表示应变挤压方向.(a)模型3的应变方向；(b)模型3和模型1的主应变方向差异，其中红色表示模型3的主应变方向，绿色表示模型1的主应变方向；(c)模型3和模型2的主应变方向差异；(d)模型3和模型4的主应变方向差异；(e)模型3和模型5的主应变方向差异. Fig. 7 The comparison of the directions of principal strains among five models A, B correspond to rectangle (A, B) with black dash line. ↔ is extension, >- < is contraction. (a) The direction of principal strain of model 3;(b) The difference of the directions of principal strain between model 3 and model 1, red arrow is the direction of principal strain of model 3, green arrow s the direction of principal strain of model 1;(c) The difference of the directions of principal strain between model 3 and model 2;(d) The d i ference of the directions of principal strain between model 3 and model 4;(e) The difference of the directions of principal strain between model 3 and model 5.

通过上节的论述能够看出在不同的耦合关系作用下, 巽他弧应变强度整体上呈现出分段性分布特征.虽然这种分段性能够用来说明巽他弧应变积累、破裂、屈服的优势方位, 但并不能够解释为什么苏门答腊南部带前弧盆地的断裂走向为北西向, 而至2004年地震主震区部分断层又偏转为近北北西向, 甚至北东向.因此本节将通过对不同耦合面应变方向上的差异来回答这一问题, 并以此为基础提出较合理的耦合模型.

图 7a为第三种耦合关系作用下前弧盆地表面应变场的主应变方向, 提取范围如图 5a中黑色虚线框所围限的A、B区, 其中A区包括了安达曼海盆的大部分区域, B区包括了苏门答腊2004年地震主震区以及南部带的北部和中部地区, 选取这两个区域的原因主要是基于该区域伴随耦合面积的变化应变差异较明显而考虑的.从图中可以看到, 在14~19区(苏门答腊南部带)主应变方向总体上呈现出北西向拉伸(平行海沟方向)、北东向挤压, 这一结果与McCary^[17]通过GPS测量以及有限元数值模拟方法所计算结果十分相似, 而事实上在该区也确实存在一系列的北西向的走滑断层以及北东向的逆冲断层^[43], 由此可见北西向拉伸的应变方向结果应该是合理的.在11~13区(苏门答腊2004年地震主震区), 主应变方向在俯冲带强耦合的作用下总体呈现出北西向拉伸、北东向挤压, 只是在局部地区偏转为近北北西向.如果假设这一计算结果是正确的, 那么该区也应发育有大量的北西向的走滑断层.但通过海底考察^[42]发现2004年地震主震区的走滑断层多为北北西和北北东向, 而且GPS观测数据也指出该区速度场方向与海沟的夹角不超过20°16], 明显小于苏门答腊南部带30°的夹角, 这从另一方面说明该区主应变方向理应不同于苏门答腊南部带.由此可见Grevemeyer^[10]、Shapiro^[12]等提出的2004年地震主震区处于强耦合状态并不符合实际情况, 该种耦合模型也并不可取.在1~10区(安达曼海盆地区), 应变场的总体拉伸方向由靠近苏门答腊中央断裂一侧的北西向逐渐偏转为靠近海沟一侧的南北向, 这一分布特征与安达曼海盆的拉张方向十分相似, 因此可视为较合理的计算结果.

图 7b为模型1和模型3的主应变方向差异对照图.从图中可以看到模型1在11~13区(2004年地震主震区)中间段的主应变拉伸方向总体表现为近北西向, 这与该区断层走向相矛盾, 因此为不合理模型.图 7c为模型2和模型3的主应变方向差异对照图, 通过对比可以发现两种模型的主应变分布特征在11~13区差异较大, 因此模型2可视为不合理模型.图 7d为模型4和模型3主应变方向差异对照图, 从图中可以看到模型4的主应变方向在15~19区总体表现为北西向拉伸, 只是在靠近海沟一侧发生了一定量偏转(与模型3的应变方向相对照).而在11~14区拉伸方向则明显发生变化, 由14区的近东西向逐渐偏转为11区的近南北向, 这一结果与该区断层走向有较好的对应.在9~10区, 由于模型4在2004年地震主震区设为无耦合状态, 因此其主应变方向在局部区域与模型3差异较大, 但总体规律变化不明显, 都呈现出近北西向的拉伸.图 7e为模型5和模型3主应变方向差异对照图.可以看到, 模型5在11~19区靠近海沟一侧的应变拉伸方向与模型4十分相似, 都呈现出由南北向逐渐偏转为东西向的分布特征.而在11~13区中间段, 由于模型在2004年地震主震区设为弱耦合状态, 因此其主应变拉伸方向由13区的北北西向逐渐偏转为11区的北东向, 虽然这一结果与模型4的结果有所差别, 但考虑其与2004年地震主震区的断层走向具有一定的相似性, 因此也可认为是一个合理模型.

通过以上分析, 能够认识到模型4和模型5的空间耦合关系能够导致2004年地震主震区应变方向上的偏转, 而该偏转与区内一系列的北北西向和北北东向的走滑断层有较好的对应.从模拟结果上来看, 本文认为造成主应变方向发生偏转的原因主要有以下四个因素: (1)俯冲板块的斜向俯冲; (2)俯冲带空间形态的改变; (3)苏门答腊南部带处于强耦合状态.(4)苏门答腊2004年地震主震区处于弱耦合或者无耦合状态.

在4.1和4.2节本文分别论述了苏门答腊前弧盆地的表面应变强度和应变方向的分布特征, 并以此为基础得出了苏门答腊2004年地震主震区应处于弱耦合状态的结论.但为什么在这种弱耦合状态下会引发海啸呢?还尚未提及.我们知道, 2004年海啸的发生主要是由于地震发生时大陆岩石圈垂向位移的瞬态变化所引发的^[44], 因此说大陆岩石圈挠曲特征将决定地震以及海啸的激发模式.为此本文在模型4和模型5中过苏门答腊前弧提取了5条剖面(剖面位置如图 5e所示), 以查看垂直位移场的分布特征.

图 8a、8b分别为模型4和模型5的垂直位移场剖面图.从图 8中可以看到, 无论处于何种耦合状态下苏门答腊地区大陆岩石圈的挠曲量和挠曲范围自南向北都呈现出逐渐递减的分布特征.对于该特征形成的原因, 本文通过不同边界条件模拟结果对比发现, 主要取决于俯冲板块的俯冲方向、俯冲速度以及俯冲角度.即, 俯冲板块的俯冲速度越大、俯冲角度越小、俯冲方向越接近垂直于海沟方向, 其上驮板块的挠曲量和挠曲范围也就越大.由此可见, 伴随俯冲方向和俯冲速度的变化, 整个巽他弧的垂直位移场也应呈现出类似的分布特征.而在苏门答腊2004年地震主震区, 由于耦合面积的变化两种模型的垂直位移场有所差异.在模型4中最大负位移(向下挠曲)出现在前弧盆地中间段, 对应于俯冲板片的最大弯曲处, 从4.2节的模拟结果上来看, 该区所对应的应变强度值也较为集中, 说明在此模型中2004年地震主震区最大负位移的出现主要受到俯冲板块推拽作用的影响.而在模型5中最大负位移出现在前弧盆地靠近海沟一侧, 对应于模型中俯冲带的耦合区.虽然从4.2节可知这种弱耦合状态对该区应变强度影响较小, 但依然能够导致上驮板块在靠近海沟一侧发生挠曲.由此可见, 俯冲带耦合面积存在与否能够影响2004年地震主震区地震、海啸的发生模式, 据此本文提出了两种可能的模式, 如图 9所示.

图 8

Fig. 8

图 8 垂直位移场剖面图 (a)表示模型4的垂直位移场，(b)表示模型5的垂直位移场. Fig. 8 Profiles of the vertical displacement field (a) is the vertical displacement field of model 4. (b) is the vertical displacement field of model 5.

图 9

Fig. 9

图 9 地震机制图 1、2、3、4分别对应地震发生前后的4个阶段. Fig. 9 Mechanism map 1, 2, 3, 4 correspond to the four stages of pre-and post-earthquake.

图 9a为模型4的地震发生模式, 可划分为4个阶段.第1阶段对应于俯冲带的初始状态, 也可表示一次大地震发生后的平静期.在第2阶段, 伴随印度一澳大利亚板块不断俯冲拖拽, 在前弧盆地的中部地区出现了一定量的挠曲.挠曲量在俯冲作用下逐渐变大直到达到破裂极限(对应于第3阶段), 从而发生弹性回跳进入新一轮的平静期(对应于第4阶段地震、海啸的发生).图 9b为模型5的地震发生模式, 与模型4的地震发生模式相似, 也可划分为4个阶段, 但与之不同的地方在于阶段2和阶段3的挠曲方式.在模型5中由于俯冲带耦合面积的存在, 因此在前弧盆地前部靠近海沟一侧首先出现了一定量挠曲.伴随俯冲作用的不断进行, 该区岩石圈挠曲量不断增大, 当达到破裂极限时, 发生地震进入平静期.然而, 对于哪一种模式更符合实际情况, 还需对海啸被激发时的空间形态进行深入的研究.

虽然上述两种模型在2004年苏门答腊地震主震区的挠曲量和挠曲位置有所差异, 但挠曲范围却十分相似, 都远小于苏门答腊南部地区.本文认为正是该区域小范围挠曲的存在为海啸的激发提供了有利的构造环境.这一特征同样说明了苏门答腊前弧的挠曲主要取决于俯冲板块的俯冲方向、角度以及速度, 而低度耦合面积存在与否只能影响局部区域的挠曲方式.

有关巽他俯冲带动力学过程对前弧盆地应变场影响的有限单元法模拟, 前人已经做了不少工作, 如Ghose等^[5]利用三维有限元方法模拟出了巽他弧应力场的分布特征, 并给出了板块作用的边界力及结构上的横向不均匀性; McCaffrey等^[17]通过二维有限元方法模拟了在不同条件作用下苏门答腊前弧盆地的应变分布特征, 并以此为基础探讨了该区的动力学演化过程, 而且对前弧盆地应变分离现象也给出了较合理的解释; 朱守彪等^[9]通过采用断层滑移技术模拟了俯冲板块与上驮板块之间的闭锁、解锁、滑动到再闭锁这一周期性过程, 即地震孕育、发生过程, 并由此提出苏门答腊大地震在时间上具有准周期性、空间上有迁移特征、破裂由深部到浅部进行.虽然上述研究成果对各自所关心的问题都做了较完满的解答, 但他们在模拟过程中同时也都忽略了一个关键问题, 即苏门答腊俯冲带内不同区域间耦合状态的相互作用方式对前弧盆地应变场的影响.因此对该区域内与耦合作用密切相关的以下构造现象的发生, 这些研究成果也并没有做更多的解释.

(1)地震破裂带的分段性.从苏门答腊地区南部近200多年来8级以上巨大地震的分布特征上来看, 其明显具有空间上的分段性, 而该种分段性在时间域和空间域上却并没有表现出一定的规律性(顺序发生), 而呈现出跳跃式发展的分布特征.即在某区域发生一次地震后, 下一次地震位置可能位于其南部, 也可位于其北部, 而两次强震的破裂区域有可能重合, 也有可能分离, 但震中位置基本上都位于俯冲带强耦合区.这就说明虽然印度一澳大利亚板块的斜向俯冲作用是整个巽他弧构造演化的直接动力源, 但对于局部地区地震的发生, 耦合面积的大小也应起着重要的控制作用.耦合面越大, 应变强度越大, 发生地震的可能性也就越大.当然, 这一影响也并不是一成不变, 还要考虑耦合面积间相互作用方式对应变积累的影响.从本文的模拟结果中也能看出, 不同的耦合面积对应不同的应变强度.尤其在苏门答腊南部带, 其南部的高应变强度值能够导致中部的应变强度值变弱, 而北部地区的应变强度能够在强耦合作用下再次升高, 总体上呈现出应变积累的分段性和跳跃性, 这一特征与该区地震的发生模式具有较好的对应.当然, 在本文模型中只是根据前人的研究成果描述了一种静态的空间耦合状态, 而在实际演化中这种空间耦合状态应该会伴随俯冲方向、俯冲速度、俯冲角度、俯冲板块的材料属性及俯冲带温度的变化而变化.当一次地震事件发生后, 同样也会导致俯冲带耦合面积的重新分配, 从而影响下一次地震的发生.因此说对俯冲带耦合状态的空间分布规律以及影响因素的研究是十分必要的, 这也是我们未来的研究方向.

(2)苏门答腊2004年地震主震区应变方向的偏转.McCaffrey^[17]和Baroux^[16]分别根据GPS数据得出苏门答腊北部地区(本文中的苏门答腊2004年地震主震区)的运动方向与苏门答腊南部地区(本文中的苏门答腊南部带)相比存在一定的偏差, 而这种偏差能够造成应变方向上的偏转.通过海底调查发现, 在苏门答腊北部地区存在一些发生偏转的走滑断层, 如Mentawei断层, 由南部的北西向逐渐偏转为北部的北北东向.对于这种现象的发生, 前人提出了不同的结论, McCaffrey^[17]认为是俯冲带弱耦合状态导致了上驮板块运动方向的改变.Baroux^[16]认为是苏门答腊北部地区的塑性行为导致了运动方向上的差异.而从本文的研究结果上来看, 影响苏门答腊北部地区应变方向偏转的因素主要有四个.首先, 该区俯冲带不应该处于强耦合状态, 如果假设该区处于强耦合状态, 其应变方向和运动方向必然会与南部地区保持一致(模型2的应变方向结果), 这将与实际情况相冲突.其次, 苏门答腊南部地区应处于强耦合状态.通过4.2节几种应变场方向的比较能够发现, 只有在南部带处于强耦合的状态下才能导致2004年地震主震区应变方向上的偏转.再次, 俯冲板块的斜向俯冲, 因为只有在斜向俯冲作用下前弧盆地才能产生应变分离或者不完全分离现象, 从而为应变方向的偏转提供了一个有利的动力学背景.最后, 俯冲带空间形态的改变.从区域地质背景图(图 1)中可以看到, 苏门答腊海沟的空间形态至2004年地震主震区时突然发生了转折, 而这一转折可能为应变方向的偏转提供了一个有利的环境背景.本文认为正是以上这四个必要因素的联合作用才最终导致苏门答腊北部地区应变方向和运动方向上的偏转.

(3)弱耦合带强震强海啸的发生.通过本文的数值模拟发现, 伴随俯冲板块的俯冲速度、俯冲角度以及俯冲方向上的变化, 模型中俯冲带不同区域间上驮板块的挠曲量和挠曲范围各不相同.据此我们拟合出了苏门答腊及其周边地区大陆岩石圈挠曲的理想模式图(图 10).从图中可以看到, 最大挠曲量和挠曲范围出现在苏门答腊南部地区.这是因为在该区印度一澳大利亚板块的俯冲速度最快、俯冲角度最小、俯冲方向与海沟夹角最大, 从而下俯板块对上驮板块的影响最大, 该结果与Franco等^[39]通过二维有限元方法模拟结果十分相似.最小挠曲量和挠曲范围出现在安达曼海盆地区, 这同样与俯冲板块的俯冲速度、方向、角度有关.而在苏门答腊2004年地震主震区情况则较为特殊, 该区大陆岩石圈的挠曲量和挠曲范围明显小于苏门答腊南部地区而大于安达曼海盆地区.本文认为正是这一挠曲特征为海啸的激发提供了有利构造环境, 原因如下:首先, 虽然该区垂向上的最大负位移与南部地区相比相对较小, 但其与周边地区的垂向位移差则较大; 其次, 虽然该区的挠曲范围较小, 但与南部地区相比其地形起伏更为剧烈, 当俯冲带发生破裂时大陆岩石圈可能更容易发生弹性回跳现象, 从而激发海啸的产生; 最后, 该区可能会受到俯冲带耦合作用的影响而改变其挠曲行为, 从而有利于海啸的产生.

图 10

Fig. 10

图 10 苏门答腊及其周边地区大陆岩石圈挠曲模式图 Fig. 10 Model showing the deflection of continental lithosphere around Sumatra

同时, 本文还存在以下一些不足之处: (1)由于目前我们无法了解研究区的应力状况, 因此忽略了初始应力场的影响, 而初始应力场可能对苏门答腊地区的应变积累、地震的发生起着重要的作用; (2)用俯冲带俯冲速率来代替重力和洋中脊推力对俯冲板片的影响, 因此在模拟过程中忽略了这两个力, 但在实际情况下这两个力可能对俯冲板块的俯冲形态及其周边地区带来重大影响; (3)在模型中没有考虑安达曼海盆中转换断层对模拟结果的影响; (4)没有考虑沉积盆地对俯冲带的影响.

总之, 本文模拟了在印度一澳大利亚板块斜向俯冲作用下整个巽他弧的应变场分布特征, 模拟结果与观测结果有很大的相似性, 同时存在着一定的差异.在未来的工作中, 建立更为复杂的、更为符合地质实际的模型, 加之综合考虑初始应力、重力、沉积盆地以及由多条断层组合而成的断裂带等的影响, 将能够逐步重现更接近真实的苏门答腊地区的构造、运动和演化动力学过程.