日本冲绳海槽和欧亚板块附加体二维构造应力场的有限元模拟及其构造意义初探

日本冲绳海槽和欧亚板块附加体二维构造应力场的有限元模拟
及其构造意义初探
卢汉堤;黄俭惠;林大五郎;木村政昭
【摘要】应用有限元方法模拟了日本冲绳海槽和欧亚大陆板块附加体二维构造应力场,并根据库仑-摩尔准则讨论了有限元模型中断裂构造发育与变形,定量化地研究了大陆岩石圈的扩张作用以及菲律宾海板块的俯冲作用对大陆板块边缘构造应力场及其变形作用的影响.数学模拟的结果表明,大陆地壳底部的扩张作用使得冲绳海槽内发育正断层,从而导致冲绳海槽这一弧后盆地的形成.菲律宾海板块的俯冲作用引起欧亚板块附加体中形成逆掩断层.有限单元模拟实验结果与研究区地震震源机制的研究结果是一致的.
【期刊名称】《桂林理工大学学报》
【年(卷),期】2008(028)001
【总页数】7页(P1-7)
【关键词】冲绳海槽;欧亚板块附加体;构造应力场;有限元方法;扩张作用;俯冲作用【作者】卢汉堤;黄俭惠;林大五郎;木村政昭
【作者单位】桂林矿产地质研究院,广西,桂林,541004;桂林矿产地质研究院,广西,桂林,541004;日本琉球大学,理学部物质地球科学系,日本,冲绳,903-0213;日本琉球大学,理学部物质地球科学系,日本,冲绳,903-0213
【正文语种】中文
【中图分类】P542.4
日本冲绳海槽的地理位置处于东中国海与琉球群岛之间，它是琉球弧和琉球海沟后侧由于大陆岩石圈的扩张作用所形成的一个弧后扩张带(图1)。

琉球弧和冲绳海槽是典型的岛弧和弧后盆地，它们共同构成大陆板块边缘［1，2］。

许多研究成果都涉及到冲绳海槽的成因问题［3-7］。

基于欧亚板块和菲律宾海板块运动的动力学机制［8］，一些地质学家提出具有不同流体学特征和边界条件的二维 (2D)有限元模型:Viallon等(1986)提出一个弹-塑性的2D有限元模型，进一步说明位于琉球海沟后侧冲绳海槽的裂陷可以用一个由于台湾岛的碰撞而产生的锚固定作用的海沟后移模型来解释［9］;Hu等 (1996)根据弹-弹 (塑)性流体力学原理，应用在平面应力条件下的2D有限元模型，分析了琉球岛弧和吕宋岛俯冲带之间发生的现代台湾岛碰撞过程中板块聚合和变形等动力学作用与应力分布状态的成因联系［10］。

此外，Hassani等 (1997)采用3种有限元模型 (线性弹性、Maxwell粘弹性和弹-塑性)模拟了板块俯冲过程中板块变形作用及其应力分布特征［11］;Berloluzza and Perotti(1997) 应用一系列有限元弹性模型模拟了与两条右行平移断层的解压超覆作用有关的拉张断陷盆地的形成过程［12］; Munoz Martin等 (1998)应用有限元方法解析了西班牙Tajo盆地东部依据断裂位移数据演绎出的第三系的三阶段古构造应力场［13］;Wang and Wang(1999)以牛顿流体介质模型研究了北美大陆板块的动力学特征［14］。

本文以有限元数学方法为手段，模拟了冲绳海槽和欧亚板块附加体二维构造应力场特征及其断裂构造的形成机制。

主要步骤为:(1)以研究区前期地球物理探测成果为基础［7，15，16］，综合整理得到冲绳海槽的P波波速结构剖面图;(2)对波速结构剖面进行线性分割，形成有限元模型的网格，沿大陆地壳底部加载扩张位移，沿菲律宾海板块俯冲带加载俯冲位移;(3)计算模型的2D和似-3D构造应力场，根据库仑-莫尔准则，分析冲绳海槽以及欧亚板块附加体中断裂构造的形成状况;(4)将数学模拟结果与研究区地震震源机制解和构造数据的综合结果进行对比。

1.1 冲绳海槽、琉球弧和琉球海沟
冲绳海槽从日本九州岛西南部一直延伸到台湾岛，该海槽在南部宽仅为60～100 km，水深约2 300 m，而在北部宽达230 km，水深逐渐降至200 m。

对冲绳海槽来说，其裂陷作用的确切时间仍是一个争论的焦点。

虽然一致认为从2 Ma至今海槽发生过两期扩张作用，但是对于早期裂陷的时代仍存在很大分歧［17］。

地
震剖面资料和地质观察结果表明冲绳海槽的裂陷作用始于中新世晚期［3-5，18］。

木村 (1985) 提出，在冲绳海槽南部的中心位置，扩张作用始于更新世早期 (约1.9 Ma)，海槽的平均扩张速度为2 cm/a(单边)［1］，冲绳海槽北部仍处于裂陷开始期［1，2］。

Miki (1995)以古地磁资料和地质年代数据为基础，针对冲绳海槽的成因提出一个两阶段的扩张模型，并认为第一阶段扩张作用发生在10～6 Ma，第二阶段扩张作用发生在约1 Ma［6］。

琉球弧从台湾岛延伸至日本九州南部，它是一个抬升的海岭，并由100多个岛屿
组成2条相互平行的岛屿链，内岛屿链广泛分布第四纪火山岩，而外岛屿链主要
出露非火山岩［19］。

整个琉球弧被Tokara Channel和Kerama Gap这两条左
侧平移断层分划成3段［20，21］。

琉球海沟通常被认为是欧亚板块与菲律宾海板块之间的边界。

除冲绳岛南部附近水深超过7 000 m外，海沟的水深一般浅于6 500 m。

海沟从南往北逐渐变得越来
越宽阔，水深越来越浅，在北部水深约为5 500 m。

1.2 地震震源机制
Fournier等(2001)以地震震源机制解和构造数据为基础，比较系统地研究了冲绳
海槽构造应力场的演化和海槽的扩张机制［22］。

其地震震源机制解和构造数据
的综合结果如图2所示，随着菲律宾海板块相对于欧亚板块进行俯冲，冲绳海槽
和琉球弧的下部地壳受到两个近于直交方向的张性变形作用，与岛弧延伸方向近于垂直的区域扩张作用和与岛弧延伸方向近于平行的局部扩张作用。

2.1 弹性模型
有限元方法是应用数学方法和计算机技术相结合的产物，也是一种广泛用来模拟大地构造变形过程的数学模型方法［23］。

在开展有限元模拟实验时，地质体可以
假设为完全弹性体［11，24，25］。

本文以平面应变条件下弹性变形的有限元模
型为基础，模拟冲绳海槽和欧亚板块附加体二维构造应力场特征及其断裂构造的形成机制。

模拟过程中构造应力场的计算由程序ELAS.F完成 (林大五郎编制)，其中考虑了大陆板块的重力效应。

2.2 模型几何特征
如图3所示，剖面S是以研究区前期地球物理探测的研究成果为基础［7，15，16］，综合整理得到的一个地震波波速结构剖面。

该剖面长约320 km，深约18 km。

其中大陆地壳的楔形沉积物由3个构造层组成，上、中和下层的P波波速分别为1.8～1.9、2.0～3.6和4.5 km/s，上层和中层被认为是砂岩，下层被认为是
灰岩。

P波波速为6.2～6.4 km/s和7.0～7.2 km/s的构造层被假定为花岗岩［26］。

在冲绳海槽南部Moho面位于海平面以下18 km深处。

菲律宾海板块相对于欧亚板块俯冲的位移速度为6.3 cm/a［8］。

2.3 有限单元网格与物理参数
模型S的有限元网格由624个节点和1 041个三角形单元组成，模型S1、S2和
S3的有限元网格与模型S相同(图4)。

模型中细、粗和中线条部分分别表示砂岩、灰岩和花岗岩，不同岩层的物理参数如表1所示［26］。

2.4 模型边界条件
有限元模型中采用2种边界条件。

扩张位移加载于冲绳海槽内大陆地壳的底部，
它代表冲绳海槽的扩张作用;由于菲律宾海板块相对欧亚板块进行俯冲，俯冲位移
加载于菲律宾海板块俯冲带。

模型S的边界条件详细叙述如下 (图4):大陆地壳底
部BG上的节点不能在垂直方向上移动，但可以在水平方向上自由移动，海槽内地
壳底部的节点加载正、负值的扩张位移(±80 m);模型左边线段AB上的节点不能在水平方向上移动，但可以在垂直方向上自由移动;节点B是固定的。

俯冲带HG 上的节点加载俯冲位移504～168 m。

冲绳海槽的扩张作用被简化为沿大陆地壳底部加载的扩张位移边界条件。

当菲律宾海板块相对于欧亚板块进行俯冲时，俯冲作用可能引起地幔对流［24］。

有限元模型中扩张位移的正、负值代表了地幔对流的2个不同方向。

这也正是模型中扩张位移的意义所在。

在有限元模型中沿菲律宾海板块俯冲带加载了俯冲位移，俯冲位移沿着整个俯冲带是变化的。

Hassani等 (1997)指出，两板块边界之间的摩擦系数是决定俯冲过程中上覆板块应力轨迹特征的一个重要因素［11］。

如果该摩擦系数沿俯冲带深度增加而增大，对抗板块俯冲的阻力也会随着深度的增加而急剧增大，从而菲律宾海板块相对于欧亚板块的俯冲作用在俯冲带深部会大大减缓，所以整个俯冲带上不同位置的位移速率就会随深度增加而减少，从而可以认为俯冲位移沿着整个俯冲带是逐渐减小的。

2.5 平面应变条件
在平面应变条件下，依据程序ELAS.F计算得到二维应力场，可以再计算得到垂直于构造剖面的第三主应力σ*:
其中ν是岩石的泊松比［27，28］。

由于有限单元的σ1、σ2和σ*已计算得到，从而可以确定最大、中间和最小主应力。

根据重新得到的3个主应力，就可以确定有限元模型的似-3D构造应力场。

2.6 断裂破裂准则
根据库仑－莫尔准则可以判定有限元模型中每个单元的断裂发育情况(图5)。

该准则可以表示为剪切应力τ和正应力σn的线性关系
其中c为岩石的聚合强度，φ为内摩擦角［29］。

当库仑－莫尔破裂线与莫尔应
力圆相切时，岩石发生断裂作用。

此时莫尔应力圆的半径(σ1－σ3)/2等于(或大于)半圆心到破裂线的垂直距离。

有限单元产生断裂的可能性参数PF可用下式计算
根据参数PF的大小可以判定有限单元中是否发生断裂作用，如果PF的值小于1.0，那将不发生断裂作用;如果PF的值达到或超过1.0，那么该单元中将产生断裂。

由于有限元模型讨论的是静力学问题，因此可以将PF值超过1.0的区域作为破裂区。

3.1 构造应力场
模型S、S1、S2和S3的似-3D构造应力场如图6所示。

扩张位移±80 m和俯冲位移504～168 m加载于模型S作为边界条件。

模型S的构造应力场表现为:冲绳海槽内σ1为近于垂直方向的压应力，σ3为近于水平方向的张应力。

在琉球弧浅
部岩层(深0～12 km)中，σ1为近于水平方向的压应力，而在深部岩层(深12～18 km)中，σ1和σ3为倾斜方向的压应力。

在欧亚板块附加体中，σ1为近于水平方
向的压应力，σ3为近于垂直方向的压应力。

模型S的构造应力场的主要特征:冲绳海槽内σ3表现为近于水平方向的张应力，欧亚板块附加体中σ1为近于水平方向
的压应力。

扩张位移±80 m和俯冲位移0 m加载于模型S1作为边界条件。

在冲绳海槽中，
模型S1的应力分布特征与模型S相似，而在欧亚板块附加体中模型S1的应力分
布特征与模型S存在很大差别，σ1在附加体中为近于垂直方向的压应力。

扩张位移0 m和俯冲位移504～168 m加载于模型S2作为边界条件，在冲绳海槽中模
型S2的应力分布特征不同于模型S，σ3在冲绳海槽表现为压应力。

而在附加体中模型S2的应力分布特征与模型S相似。

扩张位移0 m和俯冲位移0 m加载于模
型S3作为边界条件，模型S3的应力分布特征明显不同于模型S，σ3在冲绳海槽内表现为近于垂直方向的压应力，σ1在附加体中为近于垂直方向的压应力。

3.2 边界条件对构造应力场的影响
模型S、S1、S2和S3均表现出不同的构造应力场特征(图6)。

模型S和S1中加
载扩张位移±80 m作为边界条件时，σ3在冲绳海槽内表现出张应力特征;而模型S2和S3中扩张位移0 m作为边界条件时，σ3在冲绳海槽内为压应力。

虽然模型S2和S3的俯冲位移边界条件各不相同，但是在冲绳海槽它们的应力场特征表现得非常相似。

尽管模型S和S1的俯冲位移边界条件各不相同，但是在冲绳海槽它们的应力场特征也表现得非常相似。

因此，扩张位移边界条件引起冲绳海槽处于拉张应力状态，俯冲位移边界条件几乎不影响冲绳海槽内的构造应力场分布特征。

3.3 边界条件对断裂发育状况的影响
模型S、S1、S2和S3均表现出不同的断裂发育特征(图7)。

模型S和S1中加载扩张位移±80 m作为边界条件时，冲绳海槽内发育张性断裂;而模型S2和S3中扩张位移0 m作为边界条件时，冲绳海槽内不出现张性断裂。

模型S和S2中加载俯冲位移504～168 m作为边界条件时，在欧亚板块附加体中发育逆掩断层，而模型S1和S3中俯冲位移0 m作为边界条件时，欧亚板块附加体中不出现逆掩断层。

所以说，扩张位移边界条件导致冲绳海槽内产生张性断裂，而不影响附加体内逆掩断层的发育状况。

俯冲位移边界条件导致欧亚板块附加体中产生逆掩断层，也不影响冲绳海槽内断层的发育状况。

3.4 与地震震源机制解对比
有限元模拟的结果与研究区地震震源机制解得出的结论是非常相似的(图2、6)。

在模型S中，加载扩张位移±80 m和俯冲位移504～168 m作为边界条件，其构造应力场特征:在冲绳海槽内σ3表现为张应力，在欧亚板块附加体中σ1为近于水平方向的压应力。

模型S的构造应力场与研究区地震震源机制解和构造数据的综合结果是一致的［22］。

3.5 数值模拟结果
在此应用有限元方法模拟了欧亚大陆板块上部的二维构造应力场，探讨了冲绳海槽内张性断裂以及欧亚板块附加体中逆掩断层的形成机制。

根据数值模拟的结果，得
出如下结论:
冲绳海槽的扩张位移是有限元模型的一个重要边界条件。

当扩张位移边界条件加载于有限元模型时，冲绳海槽为拉张应力状态，张性断裂发育于海槽内。

扩张位移的强度控制了冲绳海槽内张性应力和张性断层的分布范围。

区域性的张性断裂往往在重力滑脱作用的影响下演化成正断层，而由于拉张地块中发育的区域性正断层往往导致地堑构造的形成［29］，冲绳海槽是由于大陆岩石圈的扩张作用所形成的。

菲律宾海板块的俯冲位移是数学模型的另一个重要边界条件。

当俯冲位移边界条件加载于有限元模型时，欧亚板块附加体中发育逆掩断层。

俯冲位移强度控制了附加体内逆掩断层的分布范围。

欧亚板块附加体内逆掩断层的形成是菲律宾海板块相对于欧亚板块进行俯冲作用的结果。

【相关文献】
［1］Kimura M.Back-arc rifting in the Okinawa Trough［J］.Mar. Pet.Geol.，1985，2:221-240.
［2］Letouzey J，Kimura M.Okinawa Trough genesis:structure and evolution of a back-arc basin developed in a continent［J］.Mar. Pet.Geol.，1985，2:111-130.
［3］Herman B M，Anderson R N，Truchan M.Extensional tectonics in the Okinawa Trough［C］//Watkins J S，Montadert L，Dickinson P W.Geological and Geophysical Investigations of Continental Margins.American Association of Petroleum Geologists Memoir 29，Tulsa，Oklahoma，1979:199-208.
［4］Lee C S，Shor G G，Bibee L D，et al.Okinawa Trough:origin of a back-arc basin ［J］.Mar.Geol.，1980，35:219-241.
［5］Sibuet J C，Letouzey J，Barbier F，et al.Back arc extension in the Okinawa Trough ［J］.J.Geophys.Res.，1987，92(B13): 14041-14063.
［6］Miki M.Two-phase opening model for the Okinawa Trough inferred from paleomagnetic study of the Ryukyu arc［J］.J.Geophys.Res.，1995，100(B5):8169-8184. ［7］Park J O，Tokuyama H，Shinohara M，et al.Seismic record of tectonic evolution and backarc rifting in the southern Ryukyu island arc system［J］.Tectonophysics，1998，294:21-42.
［8］Seno T.The instantaneous rotation vector of the Philippine seaplate relative to the
eurasian plate［J］.Tectonophysics，1977，42:209-226.
［9］Viallon C，Huchon P，Barrier E.Opening of the Okinawa basin and collision in Taiwan:a retreating trench model with lateral anchoring［J］.Earth Planet.Sci.Lett.，1986，80:145-155.
［10］Hu J C，Angelier J，Lee J C，et al.Kinematics of convergence，deformation and stress distribution in the Taiwan collision area:2-D finite-element numerical modelling ［J］.Tectonophysics，1996，255:243-268.
［11］Hassani R，Jongmans D，Chery J.Study of plate deformation and stress in subduction processes using two-dimensional numerical models［J］.Journal of Geophysical Research，1997，102(B8):17951-17965.
［12］Bertoluzza L，Perotti C R.A finite element model of the stress field in strike-slip basins:implications for the Permian tectonics of the Southern Alps(Italy)
［J］.Tectonophysics，1997，280:185-197.
［13］Munoz Martin A，Cloetingh S，De Vicente G，et al.Finite-element modelling of Tertiary paleostress fields in the eastern part of the Tajo Basin(central Spain)
［J］.Tectonophysics，1998，300:47-62.
［14］Wang S，Wang R.Joint inversion of the kinematics and dynamics of the North American plate［J］.Tectonophysics，1999，302:173-201.
［15］Kimura M.Formation of the Okinawa Trough graben(in Japanese with English abstract)［J］.The memoirs of the Geological Society of Japan，1983，22:141-157. ［16］Hirata N，Kinoshita H，Katao H，et al.Report on DELP 1988 cruises in the Okinawa Trough:part 3.crustal structure of the Southern Okinawa Trough［J］.Bull.Earthquake Res.Inst.，Univ.Tokyo，1991，66:37-70.
［17］Sibuet J C，Deffontaines B，Hsu S K，et al.Okinawa Trough backarc basin:Early tectonic and magmatic evolution［J］.Journal of Geophysical Research，1998，
103(12):30245-30267.
［18］Letouzey J，Kimura M.The Okinawa Trough:Genesis of a back-arc basin developing along a continental margin［J］. Tectonophysics，1986，125:209-230.
［19］Iwasaki T，Hirata N，Kanazawa T，et al.Crustal and upper mantle structure in the Ryukyu Island Arc deduced from deep seismic sounding［J］.Geophysical Journal International，1990，102:631-651.
［20］Wageman J M，Hilde T W C，Emery K O.Structural framework of East China Sea and Yellow Sea［J］.The American Association of Petroleum Geologists Bulletin，1970，54:1611-1643.
［21］Kobayashi K.Sea of Japan and Okinawa Trough［C］//Nairn，A E M，Stehli F G，Uyeda S.The Ocean Basins and Margins，Vol.7A.Plenum，New York，1985:419-458. ［22］Fournier M，Fabbri O，Angelier J，et al.Regional seismicity and on-land
deformation in the Ryukyu arc:Implications for the kinematics of opening of the Okinawa Trough［J］.Journal of Geophysical Research，2001，106(B7):13951-13768.
［23］Ramsay J G，Lisle R J.The Techniques of Modern Structural Geology，
Vol.3:Applications of Continuum Mechanics in Structural Geology［M］.London:Academic Press，2000: 1061.
［24］Turcotte D L，Schubert G.Geodynamics:Applications of Continuum Physics to Geological Problems［M］.John Wiley＆Sons，1982:450.
［25］Bott M H P.Stress distribution and plate boundary force associated with collision mountain ranges［J］.Tectonophysics，1990，182:193-209.
［26］Clark S P.Handbook of Physical Constants［M］.Reversed ed. The Geological Society of America，INC.，1966:587.
［27］Timoshenko S P，Goodier J N.Theory of Elasticity［M］.3rd ed.London:McGraw-Hill，1970:567.
［28］Hayashi D，Kizaki K.Numerical analysis on migmatite dome with special reference
to finite element method［J］.Journal of the Geological Society of Japan，1972，
78(12):677-686.
［29］Melosh H J，Williams C A.Mechanics of graben formation in crustal rocks:A finite element analysis［J］.J.Geophys.Res.，1989，94(B10):13961-13973.。