汤加区域地属大洋洲，位于太平洋西南部赤道附近.本区的构造特征是一个夹在太平洋板块与印度—澳洲板块之间的板块，它位于汤加海沟以西，中心大约是南纬19°-东经173°，大地构造位置特殊.本区内的俯冲带属于扩展的澳大利亚—太平洋板块边界的一部分，它反映了与太平洋板块汇聚有关的多级构造历史.前人对本地区进行了多项研究，发现汤加地区与俯冲相关的深震活动频繁(Billington, 1980；Chase, 1971；Cherkis, 1980；Dziewonski and Woodhouse, 1983；Dziewonski et al., 1983, 1988; Eguchi, 1984；Kroenke et al., 1983；Pelletier and Louat, 1989；刘超等，2009；张勇等，2009)，为本地区的构造特征的标定提供了翔实的实际参考数据.在汤加地区内的斐济群岛和斐济盆地的扩展带上具有高地震活动，其所示的地震中心在一个特定的区域内分布.Sykes(1969)发现从汤加海沟到200 km外有一个活动的地震中心.由于火山事件通常发生在地震之前(Huppert and Frohlich, 1981)，这些地震的根源可归因于岩浆通过时导致的岩石断裂.剪切波接收函数成像则表明该区域内小于200 km深度处存在板块熔融(Chen and Forsyth, 1978；Lawver et al., 1976；Malahoff et al., 1982；Mele et al., 1997; Tibi and Wiens, 2005).在汤加(Tibi and Wiens, 2005)下大约150 km深度处存在着倾斜平台，它的倾斜度为35°.该区域地震观测台站记录的数据表明，在下沉的海洋的上部或下部岩石圈可能发生部分熔融.汤加地区的地球动力学系统中观察到的许多地球动力学过程有助于研究当前本区内的复杂的构造模式(Ruellan and Lagabrielle, 2005).汤加俯冲的北部边缘受到板块活动与深部地幔上涌等动力学机制的控制，在局部区域影响板内的应力状态.海沟北部边缘的特点是太平洋底的板块的撕裂.在地震分布和震源机制测定的基础上，Millen和Hamburger(1998)认为，太平洋板块在海沟沟槽的北部的整个深度范围内俯冲，角度下降为倾角18°左右，最大可达88 km深度处，沟槽向西延伸为具有右倾走滑运动的大变形断层，其与向西太平洋板块的运动适应.基于浅层地震活动的分布(图 1)，汤加俯冲带经常被划分为对应于主要构造单元的三个区域：即板间区域、弧后区域和西北区域(Isacks et al., 1969; Kelleher and McCann, 1976; Pelletier and Louat, 1989; Sykes，1969).图 1可见来自Engdahl等(1998)给出的EHB目录(地震深度<50 km)(图 1a)(Engdahl et al., 1998)的浅地震分布以及应力状态机制解(CMT)(http://www.globalcmt.org/CMTsearch.html)(图 1b)(Dziewonski and Chou, 1981；Oliver et al., 1969).由于当前不同的构造状态即后弧分裂阶段和扩散阶段被分别研究，Kermadec-Havre(b区)和Tonga-Lau(c区)结构域也被分离开来.

图 1

Fig. 1

图 1 汤加俯冲带的分区以及地震分布 汤加俯冲带被划分为对应于主要构造单元的三个区域，即板间区域、弧后区域和西北区域(Isacks et al., 1969; Pelletier and Louat, 1989; Sykes，1969).(a) Engdahl目录(地震深度＜50 km)(Engdahl et al., 1998)的浅地震活动；(b)应力状态机制解(CMT)(Dziewonski and Chou, 1981). Fig. 1 Partition and earthquake distribution of Tonga subduction zone Tonga subduction zone is divided into three areas corresponding to the main tectonic units, i.e. interplate area, back-arc area and the northwest area (Isacks et al., 1969; Pelletier and Louat, 1989; Sykes, 1969). (a) Shallow earthquakes from Engdahl catalog (< 50 km seismic depth) (Engdahl et al., 1998); (b) Stress state from focal mechanic solutions (CMT) (Dziewonski and Chou, 1981).

汤加地区的板块结构的复杂性可以由地震波速结构指示出来(Giardini, 1984；Gill, 1976；Johnson and Molnar, 1972；Karig, 1970).根据1964年到1999年发生的5级以上的地震数据，Chen和Brudzinski(2001)对于汤加地区的板块特征进行了三维再现.Bonnardot等(2007)发现地震活动模式与应力场存在两个不同的力学单元，南部段(22°S至28°S)的特征为从60 km到700 km深的均匀的下垂压缩.这种均匀应力状态反映了中深部板块的强耦合.作为引导条件，其与应力所导致的倾斜地震带呈平行关系，并将应力从深部传递到地震带表面(Hekinian et al., 2008).这个连续的板块由地震层析成像结果给出(van der Hilst，1995).相反，在北段(从15°S至21°S)的中部深度的板块(60~300 km)处于不同的应力状态.该结果与Isacks提出的模型(Isacks et al., 1969)一致.Molnar(1971)解释了板块内应力状态变化，该部分受压，对流体物质的渗透性的阻碍作用得到加强，这种状态在深的独立板块中持续.Tibi和Wiens(2005)发现了汤加地区410 km和660 km间断面的存在以及660 km间断面的双层结构.周元泽和臧绍先(2003)也对汤加地区复杂地幔间断面结构进行了系统的研究.鲁人齐等(2013)建立了三维波速结构，观测到了本地区存在着形态不规则的板片断块.黄骥超等(2016)研究了非均匀应力场，发现其与波速结构的相关性较好，并且有效地勾勒出了复杂的构造.Okal和Kirby(1998)则认为地震活动位于太平洋板片的拆离断块内，拆离板片碎块中的相变引起了这些地震.

汤加地区的板块汇聚速度较快，深震多发，是研究地球动力学的指示性区域，对于理解俯冲动力学和地幔过程是研究的理想测试区，但是它的多级构造历史的来源不是很清晰，俯冲板片的形态的力学机制有待进一步探索.影响板块形态的力学机制是复杂的，在本研究中我们试图对一些相关的科学问题进行挖掘和探讨.对于板块的俯冲消减的角度问题，为什么沿着汤加俯冲带不是以某种固定的角度俯冲，而是在不同区域以不同的倾斜角度斜插下去，它们与本区的板块俯冲所诱导的地幔流动作用和对其上的板块产生的浮力有着什么样的相关性，从而使得这样一种相关性导致板块不垂直俯冲或呈固定的角度俯冲.另外，本地区的海沟是太平洋板块与澳洲大陆板块汇聚处，海沟位置在构造历史时期的移动与大陆板块前进的速度相当，海沟下面的大洋板块既要绕海沟转动，又要随着自身向下运动，那么这样一种关系对下伏地幔流体到底有着什么样的作用，有待进一步的分析.

本地区是研究地球动力学的指示性区域，对于理解俯冲动力学和地幔过程是研究的理想测试区.通过对汤加地区的地幔结构和运动学模式的研究，我们需要弄清研究区的复杂的俯冲模式，通过量化评估本地区的地形、板块的运动速度、岩石结构、相变和浮力比等因素对于板块的运动学模式的影响，进一步弄清俯冲板块倾角随着沟弧边缘发生变化的特征，分类探讨力学机制的影响因素，包括浮力异常、海沟等构造特征，与俯冲板块匹配的上覆板块流体的吸力等等.另外，我们试图在汤加北部重建一个二维剖面的板块的演化历史，弄清该地区的相应板块的运动方向，探索拆离块体的来源以及对后期的板块俯冲的形态有着怎样的作用，为后人对本地区进一步研究提供一定的数值模拟的参考.

由于研究地区的复杂性，要求计算模型必须具有较好的精度.因此我们先建立了全球区域的标准地幔对流模型进行测试，再对汤加地区建立量化动力学模型来模拟本区的俯冲历史.本文以地幔流体力学基本理论作为基础，通过实际数据预处理、数值分析核心算法和后处理可视化VTK的编写来建立不同边界条件下的模型，并通过与实际观测数据的对比来进行综合研究.

我们用三维的球坐标地幔对流模型有限元开源代码CitcomS (McNamara and Zhong, 2004；Tan et al., 2006; Zhong et al., 2000)来建立地幔对流和板块俯冲的模型.首先建立了全球标准地幔对流模型.对于热化学地幔对流模型的解，考察了质量、动量和能量守恒方程.这里假设地幔流体是不可压缩的，黏度结构是可变的，并且满足Boussinesq近似(除温度造成的密度变化(浮力)作为驱动力之外，其余地方的密度视为常数)，此时地幔流体力学的控制方程如下：

(1)

(2)

(3)

(4)

其中u是速度，P是动压力，T是温度，ρ_m是地幔密度，Δρ_C是组构的密度异常，ΔT是温度异常，η是黏度，α是热膨胀系数，k是热扩散系数，g是重力加速度，C是成分.相关参数值见表 1(Turcotte and Schubert, 2002；孙荀英，2013).值得注意的是，在方程(2)里，热和组构浮力都要考虑，方程(4)描述了化学颗粒的水平对流，用拉格朗日比值法求解(Tackley and King, 2003).

表 1 Table 1

表 1 模型相关参数值(Turcotte and Schubert, 2002；孙荀英，2013) Table 1 Parameter values of model (Turcotte and Schubert, 2002; Sun, 2013) 相关 参数 密度 热传导率 热扩散系数 热膨胀系数 数值 3300 kg·m-3 4 W·m-1·K-1 10-6m2·s-1 2×10-5K-1

表 1 模型相关参数值(Turcotte and Schubert, 2002；孙荀英，2013) Table 1 Parameter values of model (Turcotte and Schubert, 2002; Sun, 2013)

对于全球区域的标准地幔对流模型，我们用地幔背景黏度作为初始条件加入模型中.初始条件还包括与温度和速度相关的扰动数和瑞利数，在这里不加任何板块速度边界条件和黏度构造单元(包括海底岩石年龄等).对于网格剖分, 我们将全球分成12个块体，每个块体上单独剖分次级网格，分辨率km取决于网格数(图 2).当考虑全球尺度的模型时，网格数极大，需要并行条件，因此采用北京大学地球物理系的368核的ITAG集群进行计算(http://162.106.73.181)，模型用12个CPU并行.对于模型输出结果的原始数据，用Python编制后处理程序转换输出格式为可视化的vtk.

图 2

Fig. 2

图 2 全球区域的标准地幔对流模型的网格剖分 Fig. 2 Gridding of standard mantle convective model in the global region

研究地幔对流时，我们用温度异常分布图来反映流体的演化(图 3).由图 3可以看出，在对流开始前的稳定状态时，温度分布为同心球内的平行的层状特征(图 3a).随着时间的增加，核幔边界处的高温地幔物质上涌，并穿过周围较冷的地幔，形成地幔柱(图 3c).在演化后期，地幔柱的顶部形成帽状结构，并在接近岩石圈的底部向两侧散开(图 3e)，这也是板块运动的力学机制的来源，为下一步研究提供了模型框架.

图 3

Fig. 3

图 3 全球区域的标准地幔对流的温度场演化 Fig. 3 Temperature field evolution of standard mantle convection in the global region

根据俯冲板块的位置，我们确定了研究区的范围，并对地形数据(GTOPO30)进行处理(图 4).模型区域从140°E到150°W，60°S到10°N，深度上包括整个地幔(0~2890 km).网格剖分特征是将本地区离散化成257×257×65个节点(纬度×经度×深度)，用一个空间变量的有限元网格，以达到27 km×27 km×8 km的分辨率.模型所加的边界条件主要考虑以下因素，一方面为板块运动学和海底年龄演化，依据来自于Müller等(2008)的古板块构造重建的研究，另一方面为模型表面的速度边界条件，它随着时间而更新.除此之外，所有其他的边界条件是自由的(图 5).核幔边界的温度设置为固定的地幔温度背景值.用开源的古地理软件Gplates(http://www.gplates.org/)将板块表面速度和海底年龄按8Ma的间断来输出重构的原始数据，并用GMT进行后处理.我们也考虑汤加俯冲带的变形历史, 根据海沟剖面的形态演化来进行重构.

图 4

Fig. 4

图 4 数值模拟区域范围与地形 Fig. 4 Range and topography of the modeling area

图 5

Fig. 5

图 5 40 Ma至现在的表面速度边界条件和海底年龄分布 白色箭头显示表面速度的方向和大小.背景颜色代表海底年龄. Fig. 5 Surface velocity boundary conditions and seafloor ages from 40 Ma to present White arrows show the direction and magnitude of surface velocities. Background color represents the seafloor ages.

对于地幔黏度，在模型中将其设置为依赖于温度和组构.地幔黏度背景值设为4层(图 6)，分为岩石圈(0~100 km)，软流圈(100~300 km)，转换带(300~660 km)和下地幔(660~2890 km)，各自的黏度分别是10²³, 3×10¹⁹, 10²¹和3×10²²Pa·s，图中深度用归一化方法表示.这样的黏度剖面与如今的地幔地震结构和岩石学实验相匹配(Liu and Stegman, 2011; Hu et al., 2016).为了促进非对称俯冲，模型设定了一个沿着俯冲带边界的板间接口和一个在俯冲的板块之上的弱的地幔楔的块体，类似的设置参考以前的模型(Liu and Stegman, 2011).另外，为了模拟来自于本区内的澳洲等大陆内部克拉通的动态吸力，设置大陆单元是高黏度体.这样一种高黏度的特征调节了上地幔流体的性质，产生了改变俯冲板块形态和俯冲角度的吸力.

图 6

Fig. 6

图 6 地幔背景黏度，黏度用对数方式表示，深度用标准化表示(Liu and Stegman, 2011；Hu et al., 2016) Fig. 6 Mantle background viscosity expressed logarithmically and normalized depths (Liu and Stegman, 2011; Hu et al., 2016)

计算采用128核并行，对区域进行约束，后处理编制程序的框架与前文的全球标准模型类似，另行添加了黏度可视化的模块等.计算从0步运算到4000步，每一步代表 0.01 Ma.

图 7为计算得到的黏度演化图，地幔背景黏度值来源于图 6.在不同的位置作纵剖面(图 7b)，可以得到俯冲板块随深度变化的黏度分布情况.这里选择与汤加俯冲带垂直的中部地区的25°S处作剖面，可以看出，在板块俯冲之前的40 Ma时，模型的地幔黏度是近乎水平分布，与地幔背景黏度类似，其中大陆为高黏度体.随着时间的演化，在32 Ma时，一个黏度弱化的解耦带在表面靠西的位置出现，它指示了板块俯冲的后缘位置.演化到16 Ma时，黏度弱化带的前端出现了一个红色的三角形区域，这在构造上代表了地幔楔，它和弱化带控制了板块俯冲的位置.在演化的后期到现在，高黏度的物质穿透地幔转换区，俯冲形成规模并继续进行.同时海沟的位置在黏度演化的过程中也发生着变化，由图可以看到汤加海沟是由西向东发展的.

图 7

Fig. 7

图 7 黏度随时间的演化 (a)模型表面的黏度；(b)垂直于汤加俯冲带的25°S的黏度剖面. Fig. 7 Viscosity evolution over time (a) Viscosity of model surface; (b) Viscosity profile perpendicular to the Tonga subduction zone at 25°S.

对于俯冲结构的判断，我们需要用温度异常剖面进行研究.为了探索模型与实际观测资料的相关性，采用当前(0 Ma)的温度异常剖面与实际观测资料进行比较.图 8显示了2个不同纬度(25°S和35°S)的深度剖面以及来源于实际观测数据GAP_P4(Fukao and Obayashi, 2013)的波速结构.由图可见，模型的温度异常剖面与实际波速结构具有相关性.因此我们可以用这样的模型来探索0 Ma之前的俯冲带演化.图 9为所取俯冲剖面地区的实际震源位置分布图(Brudzinski and Chen, 2003；Chen and Brudzinski, 2001).通过不同深度的震源位置与模型俯冲板块的位置对比可以看出，模型的俯冲板块与实际的贝尼奥夫带的分布有相关性.

图 8

Fig. 8

图 8 模型温度异常与实际波速结构 (a) 25°S深度剖面；(b) 35°S深度剖面.空心圆表示地震分布，数据来源EHB(http://www.isc.ac.uk/isc-ehb/). Fig. 8 Model temperature anomalies and real velocity structure (a) Cross section along 25°S; (b) Cross section along 35°S. Hollow circles indicate earthquakes (from EHB, http://www.isc.ac.uk/isc-ehb/).

图 9

Fig. 9

图 9 汤加北部地区不同深度的震源分布(Chen and Brudzinski, 2001) Fig. 9 Seismic source distributions at different depths in the northern Tonga region (Chen and Brudzinski, 2001)

为了探索复杂结构的演化，图 10针对汤加俯冲带北部的局部区域进行研究.剖面所取的俯冲结构位于15°S附近，为两个向东西的离散的双重板片.温度剖面的演化结果显示板块俯冲发生在24 Ma之前(图 10c)，并在16—8 Ma(图 9d, 9e)之间产生分离，最终演化成现在的结构(图 10f).

图 10

Fig. 10

图 10 汤加北部地区复杂俯冲结构的演化与实际波速结构(15°S附近) 空心圆表示地震分布, 数据来源EHB(http://www.isc.ac.uk/isc-ehb/). Fig. 10 Evolution of complex subduction structure and real velocity structure near 15°S in the northern Tonga region Hollow circles indicate earthquakes (from EHB, http://www.isc.ac.uk/isc-ehb/).

我们的工作讨论了俯冲板块的形态与海底的岩石年龄、上覆俯冲板块产生的流体吸力和浮力比相关性，发现在三个不同的剖面处的海底岩石年龄、上覆板块的形态和其所产生的流体吸力都是不同的，同时还分析了板块汇聚速度对俯冲带结构形态的影响.

首先，不同的海底岩石年龄和全球俯冲板块的倾斜角度有一定的相关性(Jarrard, 1986; Lallemand et al., 2005).较老的海底岩石年龄容易引起更平的俯冲(Billen and Greg, 2007)，此时俯冲板块的强度主导了控制俯冲板块倾角的负浮力.其次，插入地幔的俯冲块体所导致的流体吸力也控制了沿着海沟的俯冲板块的倾角的变化.图 8与图 9的俯冲板块插入地幔的深度不同，35°S处的板块俯冲时间不长，俯冲板片所到达的地幔深度较浅，此时角度并没有发生太多的变化，而25°S和15°S的板块演化的时间较长，插入深度较深，较大地扰动了地幔流体，流固体产生的耦合作用导致俯冲角度发生了不同于刚开始俯冲时的变化.在遇到地幔转换区间断面时，由于有较大的相变情况发生，俯冲角度进一步的受到了相变界面阻挡的影响，俯冲板片的形态由于所处区域的不同发生了不同的变化.相应的，Stevenson和Turner(1977)，Taramón等(2015)的研究指出，板块变平至一定的临界角时会引发较平的板片和增加的吸力之间的正反馈，进一步减少俯冲板块的倾角.最后，对于浮力比特征是否会导致俯冲板块变的平整仍然存在争论.虽然早期的二维数值模式(Van Hunen et al., 2000)支持这个假设，但是最近的一项研究(Gerya et al., 2009)表明，作为对俯冲有影响作用的浮力比因素，其有助于局部水平板片的形成的条件是要同时考虑其他物理机制.我们的研究表明吸力等其他因素的存在也是必要的.

另外，对俯冲角度的其他影响因素包括上覆板块的水平运动速度和海沟后退速度，前人的研究提供了一定的佐证，例如Kameyama和Nishioka(2012)的研究指出对于地幔对流导致的俯冲的标准模型，相对俯冲速度与俯冲角度变平缓呈正相关.当俯冲速度足够大的时候，例如V_t≥3.15 cm·a^-1时，板片的水平趋势非常明显，当V_t接近于0 cm·a^-1，水平板片的趋势仅小幅出现，不同的俯冲速度与俯冲板片角度的相关性的比较指示了俯冲板块的角度和海沟的收缩后退的速度有关.对于本文的研究，由于模型用了统一的板块运动学边界条件，上覆板块的速度为Gplates给定的固定数值，没有独立的在某一固定剖面测试这些速度对板块俯冲作用的影响.可以由图 5a的平面图看到，在35°S—15°S三个剖面处的海沟的后退速度由慢到快，俯冲角度可能也受到了一些影响，对于这个现象的一个很好的理解将要求更好的实际边界条件的数据来进行约束.

我们模型的计算结果与地震波速异常剖面在一些局部的小区域也存在着一些与波速结构不太相关的地方(图 8, 9)，这可能一方面是由于地震波速结构的精确度决定，因为在汤加板块之下，地震发生的深度不是均一的，数据采集的分辨率有一定的差别.另一方面，模型本身的边界条件需要进一步加密，只有当表面板块的速度和岩石的年龄等越精确的时候，俯冲的形态才能越接近真实，而且模型的初始条件如物质构成、上覆构造单元的黏度分布和浮力比等需要进一步的精确的实际数据.

通过对汤加地区地幔结构的地球动力学模拟，我们探索了研究区局部区域的剖面上的复杂板片形态的演化，得出结论如下：

(1) 汤加地区附近的板块汇聚速度、关键构造单元的黏度结构和浮力比对板块形态和俯冲角起着较大的影响作用，俯冲板片的形态在沿海沟的方向上不同.

(2) 在近于垂直于汤加海沟的25°S的黏度剖面上可以看到，汤加地区的板块俯冲是在俯冲带西部先形成弱化带，随后形成地幔楔，并导致高黏度的俯冲板块穿透地幔转换区进入下地幔中.

(3) 在汤加地区的北部存在着一个双重结构的板块，板块的俯冲发生在24 Ma之前，板片的撕裂发生在18 Ma到6 Ma之间.

感谢北京大学的张南研究员和UIUC(University of Illinois at Urbana-Champaign，伊利诺伊大学厄巴纳-香槟分校)的LiJun Liu教授在地球动力学模型建立上的指导，感谢审稿老师和编辑的详细意见与指导.