地球物理学报  2019, Vol. 62 Issue (1): 63-77   PDF    
引用本文
史亚男, 魏东平, 皇甫鹏鹏, 等. 2019. 海洋板块俯冲作用下上覆大陆岩石层减薄机制的动力学模拟. 地球物理学报, 62(1): 63-77, doi: 10.6038/cjg2019L0785.  
Shi Y N, Wei D P, Huangfu P P, et al. 2019. Dynamics of thinning of overriding continental lithosphere induced by oceanic plate subduction: Numerical modeling. Chinese J. Geophys. (in Chinese), 62(1): 63-77, doi: 10.6038/cjg2019L0785.  
海洋板块俯冲作用下上覆大陆岩石层减薄机制的动力学模拟
史亚男, 魏东平, 皇甫鹏鹏, 李忠海, 刘梦雪     
中国科学院计算地球动力学重点实验室; 中国科学院大学地球与行星科学学院, 北京 100049
收稿日期 2017-12-27, 2018-01-30 收修定稿
基金项目: 国家自然科学基金项目（41474086，41622404，41688103）和中国科学院与国家外国专家局创新团队国际合作伙伴计划项目（KZZD-EW-TZ-19）资助
第一作者简介: 史亚男, 女, 1990年生, 博士研究生, 主要从事俯冲带数值模拟、地球动力学等方面的研究工作.E-mail:shiyanan13@mails.ucas.ac.cn
通讯作者: 魏东平, 男, 1966年生, 教授, 主要从事东亚地区微板块运动、板块动力学驱动机制、俯冲带与地震各向异性等方面的研究工作.E-mail:dongping@ucas.ac.cn
摘要：几乎所有大陆岩石层的减薄现象，可能都与海洋板块的俯冲作用相关，但是两者之间的内在联系迄今仍不十分明确，为此，我们设计了一系列包含洋-陆俯冲系统的二维数值模型，来探讨海洋板块的俯冲作用对上覆大陆岩石层变形行为的影响，尤其对大陆岩石层减薄效应的制约.模型结果表明，海洋板块俯冲过程中的地幔楔熔体对大陆岩石层地幔的热侵蚀以及由熔体上升所诱发的地幔局部对流的强烈扰动会导致上覆大陆岩石层的减薄效应.这种效应不仅表现在横向上的向陆内蔓延，还表现在垂向上的向浅部发展.且多类动力学参数都能制约大陆岩石层的减薄效应.具体地，随着汇聚速率和洋壳厚度的增加，上覆大陆岩石层在横向上的减薄范围越大，在垂向上的减薄程度也越深；而随着俯冲海洋板块年龄的增加，上覆大陆岩石层在横向上的减薄范围增大，但在垂向上的减薄程度会减小；随着上覆大陆岩石层厚度的增加，其横向减薄范围会减小，但在垂向上的减薄程度会加深.本文研究成果能为揭示华北克拉通减薄/破坏的动力学过程提供一定的理论参考依据.
关键词: 海洋板块俯冲      上覆大陆岩石层      减薄/破坏      数值模拟     
Dynamics of thinning of overriding continental lithosphere induced by oceanic plate subduction: Numerical modeling
SHI YaNan, WEI DongPing, HUANGFU PengPeng, LI ZhongHai, LIU MengXue     
Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences; College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Much evidence from geophysical, geochemical and geological studies have suggested that the thinning of continental lithosphere is generally related to oceanic subduction. While the dynamic process of such thinning remains in debate. Using a 2-D thermo-mechanical model, we systematically investigated the mechanisms of continental lithospheric thinning under oceanic subduction and explored their key constraints. The numerical models indicate that the erosion by partial melting in the mantle wedge and by the small-scale mantle convection could induce lithospheric thinning both horizontally and vertically. Moreover, the convergence velocity, thickness of oceanic crust, the age of subducting oceanic plate and the thickness of overriding continental lithosphere play important roles in the dynamics of lithospheric thinning. Specifically, with the increase of convergence velocity and the thickness of oceanic crust, the range of lithospheric thinning increases both horizontally and vertically. Besides, the range of lithospheric thinning increases horizontally and decreases vertically with the growing age of the subducting oceanic plate. Nevertheless, with increasing the thickness of overriding continental lithosphere, the range of lithospheric thinning decreases horizontally and increases vertically. Our model results could provide some insights into the dynamic processes of the thinning and destruction of the North China Craton during geologic epochs.
Keywords: Oceanic plate subduction    Overriding lithosphere    Thinning/destruction    Numerical modeling    
0 引言

全球板块运动模型中的海洋板块俯冲带，对应着地幔对流的下降流, 它是地球内部最重要的物理化学系统，这种海洋板块俯冲不仅是板块运动和洋脊扩张的最主要动力, 也是地球循环系统的重要构成因素之一, 同时提供了地球各圈层之间物质和化学元素的迁移通道(Stern, 2002).此外, 海洋板块的俯冲作用同样影响上覆大陆岩石层和地幔楔的热-动力学行为.一方面, 俯冲海洋板块的脱水作用可以降低上覆地幔楔岩体的熔点, 产生大量的幔源岩浆, 进而广泛发育活动大陆边缘的岩浆作用；另一方面, 板块富水流体或含水熔体的加入以及由此引发的地幔岩体的部分熔融，也会显著降低地幔楔岩石的流变强度和黏滞度, 从而强烈扰动地幔楔角流场, 最终发生固态上地幔岩石层的热侵蚀作用以及多种构造变形.

洋-陆俯冲作用会引发上覆大陆岩石层发生一系列复杂的变形行为，主要包括：应力应变状态变化(Stern, 2002；张克亮和魏东平, 2011)、上覆大陆岩石层的拉张和挤压(Sobolev and Babeyko, 2005；Li et al., 2013；Han et al., 2016)、全球范围内的地形地貌变化(Stern, 2002；Shi et al., 2018；Zhang and Wei, 2012)以及弧后盆地的发育和演化(Karig, 1971；Uyeda and Kanamori, 1979；Dvorkin et al., 1993).前人采用数值模拟方法对上覆大陆岩石层的动力学变形进行了广泛的研究(Van Hunen et al., 2000, 2004；Clark et al., 2008；Gerya and Meilick, 2011；Sobolev and Babeyko, 2005；Zhang and Wei, 2011；Li et al., 2015).在俯冲带，大量观测资料显示弧后地区岩石层厚度较小且为高温热状态，Arcay等(2005)采用二维数值模型研究了俯冲板块脱水对上覆大陆岩石层的减薄作用，认为这可能是由于水化地幔楔的低黏滞度的自由热对流引起的.在此基础之上，Arcay等(2008)进一步使用热-化学侵蚀模型研究了俯冲带的弧后应变特征.Gerya和Meilick(2011)基于二维岩石-热-力学耦合模型研究了洋-陆俯冲作用对上覆大陆岩石层动力学行为的影响.模型结果显示海沟后撤会导致上覆大陆岩石层流变强度的弱化，进而发生减薄作用.Gorczyk等(2007)基于二维岩石-热-力学耦合模型研究了海洋板块俯冲作用下板块中水的迁移速度对弧后扩张过程的影响.研究结果认为板块的耦合作用强烈依赖于板块汇聚速率和水的迁移速率之比.Clark等(2008)采用三维数值模型分析了海洋板块俯冲作用下弧后区域的阶段性，将其分为三种类型：无阶段性、类阶段性和超阶段性.Sobolev和Babeyko(2005)利用热-力学耦合模型研究了上覆南美大陆挤压缩短的动力学制约因素，认为俯冲界面上的剪切耦合作用是控制其缩短的主要因素.

关于海洋板块俯冲作用对上覆大陆岩石层的影响，其中比较特殊的是对克拉通(岩石层厚度~200 km)减薄/破坏的影响.二维数值模拟结果显示，地壳榴辉岩的重力失稳可能是引起岩石层大规模减薄的原因(乔彦超等，2012)；更进一步，具有不同流变性质的克拉通重力失稳过程的数值模拟结果发现，这种岩石层减薄/破坏具有分期并且多阶段的特征(Wang et al., 2015).而对下地壳榴辉岩的黏性、密度以及岩石层地幔黏性等参数的敏感性试验结果则表明，岩石层地幔的黏性会显著影响拆沉时间, 影响范围从几个到十几个百万年不等(程华冬等，2017).另一方面，乔彦超等(2013)通过数值模型研究了岩石层的热侵蚀过程, 发现小尺度的地幔对流现象可以使岩石层发生减薄，其减薄量可以达到100 km.程华冬等(2016)利用二维数值模型研究了克拉通岩石层的对流减薄过程, 讨论了岩石层宽度、厚度、黏性差异以及密度差异等参数的动力学制约，模型结果显示黏性差异对大陆岩石层减薄的时间尺度有很大影响.He(2014)采用二维数值模型研究了俯冲板块与大地幔楔之间的关系，发现板块俯冲过程可以促进低黏性大地幔楔的形成.Wang等(2016)采用数值模拟的方法研究了水化作用在大陆岩石层减薄中的作用，模型结果显示水化作用与克拉通地区大规模的岩石层减薄过程是密切相关的.

随着研究的日益深入，已有越来越多的学者认识到，海洋板块的俯冲作用可能是诱发上覆大陆岩石层减薄/破坏的主要控制因素，并用各有特点的地质模型进行数值模拟.但海洋板块俯冲如何导致上覆大陆岩石层发生减薄，甚至大规模破坏，其一般性的作用机理是什么，仍存在较大争议.为此，本文试图通过建立二维洋-陆俯冲数值模型, 来探讨海洋板块俯冲过程中，其地幔楔内热-熔体扰动引发的不稳定地幔角流体系对上覆大陆岩石层的影响, 并在此基础上进一步对上覆大陆岩石层减薄的动力学制约因素进行讨论.

1 数值模拟方法介绍 1.1 控制方程

对于板块俯冲碰撞相关的动力学数值模型，一般对三组控制方程进行求解，包括斯托克斯流体动力学方程、物质守恒方程以及热量守恒方程.本文中，将采用I2VIS程序, 对这些方程在不规则的欧拉网格节点上采用有限差分算法进行离散化和近似求解，具体的算法可见Gerya和Yuen(2003a)或Li(2014).

1.1.1 斯托克斯方程

(1)

其中, g是重力加速度, 密度ρ依赖于温度T、压力P、岩石类型C和部分熔融比例M, σ′xx, σ′yy和σ′xy是偏应力张量.

1.1.2 非可压缩流体物质守恒方程

(2)

其中, vx和vy分别表示水平速度和垂向速度.

1.1.3 热量守恒方程

(3)

其中, Cp是等压热容, DT/Dt表示温度对时间的物质导数, qx, qy代表热流值, k是热传导率，依赖于温度T, 压强P和岩石类型C.

Hr表示放射性生热，假设为依赖于岩石类型的常量：

(4)

Ha表示绝热变压生热，即伴随着压力的变化(增压或减压)而导致的热量生成或消失：

(5)

HS表示摩擦生热，代表了黏性形变中机械能转化为热能：

(6)

其中，σ′xx, σ′yy, σ′xy是偏应力张量，是应变率张量.

1.2 黏-塑性流变学性质 1.2.1 黏性形变

数值模型中的流变学关系采用整合的黏-塑性本构关系.其中，韧性流变的黏滞系数是一个与温度、压力、物质成分、应变率和熔融程度相关, 可表示为

(7)

式中, 是应变率张量的二阶不变量；AD、E、V和n是实验岩石学所得到的流变参数，分别表示物质常数、活化能、活化体积和应力指数；F是根据实验类型确定的无量纲系数；R是气体常数.

1.2.2 塑性形变

上述的黏性流变需要与塑性流变相结合，从而形成实际的黏-塑性流变特征，这里采用改进的Drucker-Prager屈服准则(如Ranalli, 1995)：

(8)

该方程中, σyield表示屈服应力；表示应变率张量二阶不变量；P是总压力；C0是P=0条件下的岩石剩余强度；φ是内摩擦角；λ是孔隙流体系数(Brace and Kohlstedt, 1980).

1.2.3 形变机制的整合

基于ηductile和ηplastic, 对于模型区域中的某一点来讲，其黏-塑性流变关系的最终等效黏滞系数可以定义为它们之中的最小者(Ranalli，1995)：

(9)

1.3 部分熔融

基于实验岩石学的约束条件, 本算法考虑了地壳岩石的部分熔融行为(Gerya and Yuen, 2003b；Burg and Gerya, 2005).该算法中近似认为部分熔融体积比例M与温度间存在如下线性关系：

(10)

式中, Tsolidus和Tliquidus分别代表特定岩性的湿固相线温度和干液相线温度.

部分熔融岩石的有效密度取决于熔融比例：

(11)

式中, ρsolid和ρmolten分别代表固相和熔融岩石的密度, 是关于温度和压力的函数：

(12)

式中, ρ0代表岩石在P0=0.1 MPa和T0=298 K温压条件下的标准密度；α和β分别代表热膨胀系数和可压缩系数.

2 初始模型和边界条件

基于要研究的海洋板块俯冲作用下上覆大陆岩石层减薄机制的动力学背景，我们设计了多组海洋板块向上覆大陆板块俯冲的数值模型，模型空间宽度为4000 km，深度为670 km.参考模型中洋壳厚度为8 km，大陆地壳厚度为35 km，由20 km的上地壳和15 km的下地壳组成.俯冲海洋板块的年龄为60 Ma，上覆大陆岩石层的厚度为120 km.在地壳表面之上，与自由滑动的模型顶界面之间，设计有一层相对高黏滞度的伪空气层，其与上地壳的接触面用以模拟地形起伏面，该地形起伏包含了近似的剥蚀和沉积作用(Li and Gerya, 2009).数值模型中所采用的沉积和剥蚀速率均设置为0.3 mm·a-1(Gerya and Yuen, 2003a；Burg and Gerya, 2005).数值模型的初始物质场和温度场设计如图 1所示, 模型采用的黏滞性流变参数见表 1，模型中的主要材料参数见表 2.

表 1 模型采用的黏滞性流变参数 Table 1 Viscous flow parameters used in numerical experiments
表 2 模型中的主要材料参数 Table 2 Material properties used in numerical experiments

模型中除了底部边界，其他边界均采用自由滑动的速度边界条件.底部边界是渗透性边界，采用近无限深度的外部自由滑动边界条件，这意味着在计算模型区域的下面满足自由滑动边界条件(如Li et al., 2016).与普通的自由滑动边界条件相同，外部自由滑动条件也将满足计算区域内的物质守恒.

图 1 初始数值模型和边界条件 模型空间尺度为4000 km×670 km, 数值模型的网格节点采用不均一的空间分辨率, 俯冲核心区域的分辨率为2 km×2 km, 并向两端逐渐减小至30 km×30 km.白色线条是等温线, 单位是℃.模型中颜色代表不同的岩石组成, 岩性色标如下：1空气；2水；3大陆上地壳；4大陆下地壳；5洋壳；6岩石层地幔；7软流层地幔；8初始薄弱带；9, 10沉积物；11沉积物部分熔融；12大陆上地壳部分熔融；13, 大陆下地壳部分熔融；14, 洋壳部分熔融. Fig. 1 Initial numerical model and boundary conditions Model size is 4000 km×670 km. Grid cell is 2 km×2 km in the central subduction zone and 30 km×30 km elsewhere. White lines are isotherms in ℃. Colors indicate medium types. 1 Air; 2 Water; 3 Upper continental crust; 4 Lower continental crust; 5 Oceanic crust; 6 Lithosphere mantle; 7 Asthenosphere mantle; 8 Initial weak zone; 9 and 10 Sediment; 11 Partial molten sediment; 12 Partial molten upper continental crust; 13 Partial molten lower continental crust; 14 Partial molten oceanic crust.

对于热量守恒方程的边界条件，模型顶部为固定温度(0 ℃)，大陆岩石层底边界温度为1300 ℃(Turcotte and Schubert, 2002)；两侧边界的水平方向温度梯度为零(即零热流).底部边界采用的是外部边界固定温度条件，关于外部边界的理论和应用可见Li等(2010, 2016).模型自初始条件开始，俯冲岩石圈受到一个向右的推力作用，该推力主要通过在俯冲海洋板块的远端施加一个固定的速率(Vx)而实现(图 1).

3 模型结果 3.1 参考模型

参考模型中, 俯冲海洋板块的汇聚速率为6 cm·a-1, 洋壳厚度为8 km, 海洋板块的年龄为60 Ma, 上覆大陆岩石层厚度为120 km, 该模型演化结果如图 2所示.海洋板块俯冲过程中, 洋壳和大洋沉积物质不断进入俯冲带, 受上覆地幔楔热传导加热的影响发生部分熔融(近似认为是实际地质过程中交代地幔岩的部分熔融), 并在自身正浮力作用下逐渐脱离俯冲板块进入地幔楔(图 2a), 最终到达大陆岩石层底部.当俯冲时间到达16.2 Myr时, 在快速上升熔体的冲击、由此诱发的地幔楔内局部对流及熔体和软流层物质对岩石层地幔底部热侵蚀的共同作用下, 上覆岩石层底部受到一定程度的减薄/破坏(图 2b).随后，熔体的持续上升以及由此引发的地幔楔局部对流作用的增强导致上覆岩石层底部减薄/破坏逐渐向浅部扩张, 减薄程度逐渐加深；同时, 上覆岩石层地幔减薄也逐渐向陆内方向扩展.模型运行至24.0 Myr时, 上覆大陆岩石层在垂向上减薄约30 km, 而横向上的减薄范围可达到约300 km, 如图 2c-2d所示.对应的黏滞度场清晰显示, 板块熔体的加入以及由此诱发的局部对流作用共同导致上覆地幔楔黏滞度降低约1/1000~1/100.随着岩石层地幔减薄向浅部和向陆内双重发展, 具有低黏滞度的软流层物质也相应的向浅部和陆内方向蔓延, 如图 2a′-2d′所示.

图 2 参考模型中的物质场(a—d)和等效黏滞系数场演化(a′—d′) t表示板块俯冲时间, ηeff表示等效黏滞系数.模型参数：汇聚速率为6 cm·a-1, 洋壳厚度为8 km, 俯冲海洋板块年龄为60 Ma, 上覆大陆岩石层厚度为120 km. Fig. 2 Evolution of the reference model. The color shows the evolution of material field (a—d) and the viscosity structure (a′—d′), respectively ηeff is the effective viscosity. Time (Myr) of subduction is given in each panel. Convergence velocity is 6 cm·a-1. Thickness of oceanic crust is 8 km. Age of subducting oceanic plate is 60 Ma. Thickness of overriding continental plate is 120 km.
3.2 汇聚速率对上覆大陆岩石层减薄的影响

在参考模型的基础上, 我们分别设计了一个慢速汇聚模型(汇聚速率Vx=2 cm·a-1)和一个快速汇聚模型(汇聚速率Vx = 10 cm·a-1), 来探讨汇聚速率对上覆大陆岩石层减薄的影响.模型结果如图 3所示, 具体而言, 俯冲海洋板块的汇聚速率分别为2 cm·a-1、6 cm·a-1和10 cm·a-1的三个模型在24.0 Myr时，上覆大陆岩石层在横向上的减薄范围分别达到约50 km、300 km和400 km, 在垂向上的减薄深度依次为10 km、20 km和30 km(图 3a-3c).此外，当三个模型岩石层的俯冲汇聚量相同，即汇聚速率分别为2 cm·a-1、6 cm·a-1和10 cm·a-1的三个模型分别在72.0 Myr，24.0 Myr和14.4 Myr时，上覆大陆岩石层在横向上的减薄范围分别达到约200 km、300 km和400 km, 在垂向上的减薄深度依次为18 km、20 km和30 km(图 3a′-3c′).模型结果显示随着海洋板块汇聚速率的增加, 上覆大陆岩石层的减薄范围在横向上显著增大, 且垂向上的减薄深度随之加深.可以看出，在汇聚速率不同的情况下，无论俯冲汇聚量是否一致，模型结果的变化趋势均保持一致.所以相对俯冲汇聚量而言，汇聚速率对上覆大陆岩石层减薄的影响占据主要地位.

图 3 (a—c)俯冲海洋板块的汇聚速率依次为2 cm·a-1, 6 cm·a-1, 10 cm·a-1的三个模型在汇聚时间t=24.0 Myr时的物质场演化结果；(a′—c′)三个模型在俯冲汇聚量相同时的物质场演化结果，汇聚时间分别为72.0 Myr、24.0 Myr和14.4 Myr.其他参数设置与参考模型中一致 Fig. 3 (a—c) Evolution of models of various convergence velocities when the time (Myr) of subduction is 24.0 Myr in each model. (a′—c′) Evolution of models of various convergence velocities when the length of subduction is same in each model. The time (Myr) of subduction is 72.0 Myr, 24.0 Myr and 14.4 Myr respectively. Other parameters are identical to the reference model

造成这种趋势的原因是, 在海洋板块汇聚速率较低的情况下, 俯冲洋壳和沉积物在俯冲带浅部就能被强烈增温弱化而发生部分熔融, 熔体随即进入地幔楔并聚集在大陆岩石层地幔底部, 这样就导致进入软流层深部的洋壳和沉积物的减少, 造成深部板块熔体的减少以及对地幔楔角流场扰动的减弱, 并最终体现在上覆大陆岩石层地幔减薄的水平范围的减小.这与冷的海洋板块俯冲带相对于热俯冲带更能发育弧岩浆作用有相似机制.在冷俯冲带, 俯冲洋壳在弧前地幔浅部不能发生显著脱水, 造成俯冲进入弧下地幔深度的洋壳仍然能析出大量流体交代上覆地幔楔, 随之导致大量弧岩浆的发育(Peacock and Wang, 1999;Syracuse et al., 2010)；相反地, 在热俯冲带, 洋壳在弧前地幔浅部就能发生显著脱水, 致使俯冲进入弧下地幔深度的洋壳很难析出大量流体交代地幔楔, 进而导致弧岩浆作用的相对缺乏(郑永飞等, 2015).

3.3 俯冲洋壳厚度对上覆大陆岩石层减薄的影响

在参考模型的基础上, 我们将俯冲洋壳的厚度分别设置为4 km、8 km、12 km和16 km, 来探讨俯冲洋壳厚度对上覆大陆岩石层减薄的影响.模型结果如图 4所示, 随着俯冲洋壳厚度的升高, 上覆大陆岩石层在横向上的减薄范围显著增加, 同时岩石层地幔在垂向上的减薄程度也相对加深.具体而言, 在模型运行至24.0 Myr时, 洋壳厚度分别为4 km、8 km、12 km和16 km的四个模型所导致的上覆大陆岩石层地幔在横向上的减薄范围分别达到约280 km、300 km、320 km和350 km, 垂向上的减薄程度依次约为10 km、30 km、40 km和50 km(图 4a-4d).相应的黏滞度场也清晰显示, 随着俯冲洋壳厚度的增加, 上覆地幔楔黏滞度降低约1/10000~1/1000, 具有低黏滞度的软流层物质也相应的向浅部和陆内方向扩张, 如图 4a′ -4d′所示.此外, 增加俯冲洋壳的厚度会对地幔熔体的分布区域和体积含量产生一定影响.随着洋壳厚度的增加, 会直接导致熔体体量的增加, 并且相应造成聚集在大陆岩石层地幔底部的熔体体积的增加(图 4c, 4d).重要的是，板块物质的拆离、部分熔融及向上迁移程度的增加，伴随着由此诱发的地幔楔内局部对流作用的增强, 从而造成上覆大陆岩石层在横向上减薄范围的扩大以及垂向上减薄深度的加深.

图 4 洋壳厚度依次为(a) 4 km, (b) 8 km, (c) 12 km, (d) 16 km的四个模型在汇聚时间t= 24.0 Myr时的物质场(a—d)和等效黏滞系数场(a′—d′)演化结果.其他参数设置与参考模型中一致 Fig. 4 Evolution of models of various thicknesses of oceanic crust. (a) 4 km, (b) 8 km, (c) 12 km, (d) 16 km. The color shows the evolution of material field (a—d) and the viscosity structure (a′—d′), respectively. Time (Myr) of subduction is 24.0 Myr in each model. Other parameters are identical to the reference model
3.4 俯冲海洋板块年龄对上覆大陆岩石层减薄的影响

在参考模型的基础上, 我们将俯冲海洋板块年龄分别设置为20 Ma、40 Ma、60 Ma、80 Ma和100 Ma, 来探讨其对上覆大陆岩石层减薄的影响.在我们的数值模型中，海洋板块年龄的改变会直接影响海洋岩石层的厚度，同时初始温度结构也会发生变化.模型结果如图 5所示, 随着俯冲海洋板块年龄的升高, 上覆大陆岩石层在横向上的减薄范围显著增加, 但在垂向上的减薄深度却有所减弱.具体而言, 在模型运行至24.0 Myr时, 俯冲海洋板块的年龄分别为20 Ma、40 Ma、60 Ma、80 Ma和100 Ma的五个模型中的上覆大陆岩石层地幔在横向上的减薄范围分别达到约210 km、260 km、300 km、330 km和360 km, 同时垂向上的减薄深度依次约为50 km、40 km、30 km、20 km和10 km(图 5a-5e).相应的软流层黏滞度场也同样体现出岩石层减薄的横向变化和垂向变化(图 5a′-5e′).产生这种变化的原因是, 当俯冲海洋板块年龄较小时, 俯冲板块比较热, 俯冲洋壳和沉积物在俯冲带浅部就能被强烈增温弱化、脱水，致使地幔楔内部分熔融的发生, 随后楔入大陆岩石层地幔之中, 从而导致上覆大陆岩石层地幔减薄的水平范围相对较小, 但在垂向上的减薄深度较深.相反地, 当俯冲海洋板块的年龄比较大时, 俯冲板块比较冷, 俯冲洋壳在弧前地幔浅部没能显著增温, 进而俯冲至弧下地幔深度才逐渐增温、脱水，发生部分熔融，进而引发更大规模地幔对流，最终导致上覆大陆岩石层地幔减薄的水平范围相对较大, 但垂向上的减薄深度相对较浅.

图 5 俯冲海洋板块年龄 依次为(a) 20 Ma; (b) 40 Ma; (c) 60 Ma; (d) 80 Ma; (e) 100 Ma的五个模型在汇聚时间t=24.0 Myr时的物质场(a—e)和等效黏滞系数场(a′—e′)演化结果.其他参数设置与参考模型中一致. Fig. 5 Evolution of models of various ages of subducting oceanic plate (a) 20 Ma; (b) 40 Ma; (c) 60 Ma; (d) 80 Ma; (e) 100 Ma. The color shows the evolution of material field (a—e) and the viscosity structure (a′—e′), respectively. Time (Myr) of subduction is 24.0 Myr in each model. Other parameters are identical to the reference model.
3.5 上覆大陆岩石层的厚度对其减薄的影响

在参考模型的基础上, 我们将上覆大陆岩石层的厚度分别设置为100 km、120 km、140 km、160 km、180 km和200 km, 来探讨上覆大陆岩石层厚度对其自身变形行为的影响.模型结果如图 6所示, 随着上覆大陆岩石层的厚度的增加, 上覆大陆岩石层在横向上的减薄范围会有所减小, 但是在垂向上的减薄深度却会显著加深.具体而言, 模型运行至24.0 Myr时, 厚度分别为100 km、120 km、140 km、160 km、180 km和200 km的上覆大陆岩石层在横向上的减薄范围分别达到约320 km、300 km、280 km、260 km、230 km和200 km, 在垂向上的减薄深度依次约为20 km、30 km、40 km、50 km、60 km和70 km(图 6a-6f), 对应的软流层黏滞系数场也会有相关的反映(图 6a′-6f′).增加上覆大陆岩石层的厚度会对板块熔体的分布区域和熔体的含量以及由此诱发的地幔局部对流程度产生一定影响.在俯冲海洋板块的年龄一定的情况下, 当上覆大陆岩石层相对较薄时, 拆离的洋壳和沉积物可聚集至大陆岩石层的底部, 致使大陆岩石层底部横向减薄范围相对较大, 但是垂向上的减薄深度相对较浅；而当上覆大陆岩石层比较厚时, 洋壳和沉积物拆离和部分熔融在大陆岩石层的中部即开始发生, 从而导致大陆岩石层横向减薄范围相对较小, 但在垂向上的减薄深度相对较深.

图 6 上覆大陆岩石层厚度 依次为(a)100 km; (b) 120 km; (c) 140 km; (d) 160 km; (e) 180 km; (f) 200 km的六个模型在汇聚时间t=24.0 Myr时的物质场(a—f)和等效黏滞系数场(a′—f′)演化结果.其他参数设置与参考模型中一致. Fig. 6 Evolution of models of various thicknesses of overriding continental plate (a) 100 km; (b) 120 km; (c) 140 km; (d) 160 km; (e) 180 km; (f) 200 km. The color shows the evolution of material field (a—f) and the viscosity structure (a′—f′), respectively. Time (Myr) of subduction is 24.0 Myr in each model. Other parameters are identical to the reference model.
4 讨论 4.1 海洋板块俯冲为何能造成上覆大陆岩石层底部减薄？

数值模拟试验结果显示，海洋板块的俯冲作用通常会造成上覆大陆岩石层一定程度的减薄.造成这一现象的原因主要包括两方面：(1)被俯冲板块带至上地幔深度的洋壳和沉积物在热传导作用下迅速增温、弱化、脱水，导致部分熔融的发生，向上运移，并最终聚集在大陆岩石层地幔底部；而此过程中，上升的熔体对大陆岩石层底部的硬“冲击”和热侵蚀作用可以促使岩石层地幔底部减薄现象的发生；(2)更重要的是，在熔体上升过程中所诱发的地幔楔内部局部对流现象的产生，能加速大陆岩石层地幔底部物质的迁移、扩散，从而造成大陆岩石层进一步的减薄，甚至破坏.

另外，从海洋板块俯冲的地质演化来看，海洋板块俯冲过程中洋壳和大量沉积物被拖曳进入俯冲隧道，并在地幔深度处发生显著增温弱化和大量脱水，引起上覆地幔楔物质的水化蚀变, 并最终导致蚀变地幔的部分熔融.地幔熔体物质随后向上运移，到达大陆岩石层地幔底部时，在热侵蚀作用下造成其一定程度减薄.同时，地幔楔内部的局部对流作用又进一步加剧大陆岩石层减薄.并且这种减薄效应一方面表现在横向上空间减薄范围的扩张, 另一方面表现为垂向上向浅部发展(图 2d).

4.2 海洋板块俯冲对大陆岩石层减薄效应的范围

持续的海洋板块俯冲对上覆大陆岩石层的减薄效应可以表现为横向上的减薄范围和垂向上的减薄程度，并且这种减薄效应会随着动力学制约因素的改变而发生变化.

汇聚速率变化能显著影响上覆大陆岩石层的减薄效应.模型结果显示，随着汇聚速率从2 cm·a-1增加到10 cm·a-1，上覆大陆岩石层在横向上的减薄范围可以从50 km扩展到400 km, 在垂向上的减薄程度可以从10 km增加到30 km(图 3a-3c).在快速俯冲的条件下, 上覆大陆岩石层无论是在横向还是垂向上都会发生显著的减薄现象.

洋壳厚度变化同样制约上覆大陆岩石层的减薄效应.随着洋壳厚度从4 km到16 km的增加, 上覆大陆岩石层在横向上的减薄范围可以从280 km增加到350 km, 在垂向上的减薄程度可以从10 km增加到50 km, 聚集在大陆岩石层地幔底部的熔体含量也会随之增加(图 4a-4d).

随着俯冲海洋板块年龄的增加，上覆大陆岩石层减薄的横向范围和垂向程度表现出截然相反的变化趋势.模型结果显示，当海洋板块年龄从20 Ma增加到100 Ma时, 上覆大陆岩石层在横向上的减薄范围随之增加, 从210 km逐渐扩张到360 km, 但是, 对应的上覆大陆岩石层在垂向上的减薄程度则随之减小, 从50 km逐渐减小到10 km(图 5a-5e).造成这一现象的原因主要是，随着海洋板块年龄的增加，洋壳脱水和由此诱发的地幔楔部分熔融的深度也会逐渐增加，进而制约着由熔体上升迁移而引发的地幔楔内部局部对流的强度，并进一步表现在对上覆大陆岩石层的减薄效应上.

上覆大陆岩石层自身的厚度也会影响其在海洋板块俯冲作用下自身的变形行为和减薄效应.模型试验表明，随着上覆大陆岩石层厚度从100 km增加到200 km，由海洋板块俯冲作用所导致的上覆大陆岩石层减薄效应横向上的减薄范围会从320 km减小到200 km, 但是在垂向上的减薄程度却可以从20 km增加到70 km(图 6a-6f).

4.3 洋-陆俯冲系统中大陆减薄的地貌响应

由海洋板块俯冲作用而导致的上覆大陆岩石层减薄效应在大陆地貌演化上也会表现出一定响应(如Li et al., 2015).我们以参考模型结果为例，来简要阐述洋-陆俯冲系统中大陆减薄的地貌变化.图 7是参考模型地形演化的结果.从初始俯冲至20.0 Myr，海洋板块俯冲作用会在大陆板块边缘形成一狭长造山带.但随着海洋板块俯冲的持续进行，造山带会逐渐向陆内扩展.在俯冲时间到达25.0 Myr时, 距离海沟约200 km内的大陆显著抬升，平均达4 km；而当模型运行至40.0 Myr时, 距离海沟500 km内大陆范围内都发生显著的地形抬升.结合图 2所示的参考模型的岩石层尺度演化结果可知，在海洋板块俯冲作用下的大陆边缘造山带向陆内扩展是大陆岩石层减薄效应的浅表响应.

图 7 参考模型下的地形演化 颜色代表地形起伏, 如右侧色标.横轴是截取的模型宽度：2500≤x≤3700 km.纵轴是时间轴：0≤t≤40.0 Myr.模型参数(与图 2相同)：汇聚速率为6 cm·a-1, 洋壳厚度为8 km, 俯冲海洋板块年龄为60 Ma, 上覆大陆岩石层厚度为120 km. Fig. 7 Evolution of the topography of the reference model The axes: horizontally 2500≤x≤3700 km and vertically 0≤t≤40.0 Myr. Model parameters are same as Fig. 2. Convergence velocity is 6 cm·a-1. Thickness of oceanic crust is 8 km. Age of subducting oceanic plate is 60 Ma. Thickness of overriding continental plate is 120 km.

以华北东部古地貌为例(吴智平等, 2007；朱日祥和郑天愉, 2009；陈斌等, 2013)，中生代以来经历了印支运动、燕山运动和喜马拉雅运动三次主要的构造运动.起始时期，华北地区经历了海西运动，接受了广泛的陆表海沉积.早-中三叠世华北东部地区基本与晚海西期以来的构造格局保持一致，继续接受沉积，形成了大型的华北内陆盆地.中三叠世末期的印支运动使得华北地区东部整体抬升，晚三叠世的华北盆地退缩到太行山以西，东部开始遭受剥蚀.晚侏罗世-早白垩世时期，太平洋板块的活动对华北地区的构造演化占据控制地位，造山带增高，形成燕山期高原，同时华北东部盆地的发育主要以隆起在高原之上的山间盆地为主.到晚白垩世时期，华北地区巨大的岩石层遭受各种侵蚀作用发生了拆沉作用，华北地区整体抬升，全区遭受剥蚀，东部造山带坍塌，前期东高西低的古地形不复存在，地形上整体夷平.中新世以来，东部地区的断陷作用减弱，进入断坳转换期.可以看出，在太平洋板块俯冲之后的后撤作用下，华北东部岩石层受到强烈的拉张作用，从而发生了大规模的减薄.因此，现今华北东部的地形是经过了强烈的伸展改造后的结果，故与本文地形模拟结果所得到的初始高原和造山带(图 7)不一致.可以推测，到达海洋板块俯冲后期，板块后撤，上覆大陆岩石层受到强烈的拉张作用，随之可以发生沉降.

通过系列模型试验结果，我们知道多类动力学参数均能制约由持续海洋板块俯冲而诱发的大陆岩石层减薄效应，且这些参数值的变化影响着大陆岩石层减薄的横向范围和垂向程度，进而影响着造山带演化.以汇聚速率的变化为例，对于海洋板块的汇聚速率分别为2 cm·a-1、6 cm·a-1和10 cm·a-1的三个模型，在模型运行至24.0 Myr时, 上覆大陆岩石层的地貌抬升高度最高可以达到4 km、6 km和7 km(图 8).并且随着汇聚速率的增大, 海沟的位置会逐渐向陆内迁移, 同时上覆大陆岩石层发生变形的区间也会随之向陆内更远处扩展.地形演化的结果也进一步表明, 随着汇聚速率的升高, 上覆大陆岩石层在垂向上的减薄深度会加深, 横向上的空间减薄范围也会发生显著扩张.

图 8 海洋板块汇聚速率分别为2 cm·a-1、6 cm·a-1、10 cm·a-1的三个模型在汇聚时间t=24.0 Myr时的地形演化结果.其他参数设置与参考模型中一致 Fig. 8 Evolution of the topography of models with various convergence velocities: 2 cm·a-1, 6 cm·a-1, and 10 cm·a-1. Time (Myr) of subduction is 24.0 Myr in each model. Other parameters are same as the reference model
4.4 对华北克拉通减薄/破坏的启示

最近二十年来, 多个研究小组都对华北克拉通减薄/破坏的动力学机制及时空演化过程进行了研究和争论，主要包括：拆沉作用(Gao et al., 1998, 2002)；热侵蚀作用(Menzies and Xu, 1998；Xu, 2001, 2007；Xu et al., 2008)；橄榄岩-熔体相互作用(Zhang et al., 2002)；机械拉张作用(吴福元等, 2008)；岩浆提取作用(Chen et al., 2004)以及岩石层地幔水化作用(Niu, 2005)等.通过对中国东部岩石层减薄的时间、范围以及机制等问题进行讨论, 吴福元等(2003)认为岩石层减薄发生在晚中生代, 主要与古太平洋板块的俯冲有关，并且这种板块俯冲机制一定程度上可以推广到全球范围的克拉通减薄/破坏(吴福元等，2014).局部区域的华北克拉通曾经受到古特提斯洋和古亚洲洋俯冲的影响(Windley et al., 2010；陈斌等, 2013)，大范围的华北克拉通减薄/破坏则主要与太平洋的西向俯冲有关(徐义刚等, 2009；郑永飞和吴福元, 2009；朱日祥和郑天愉, 2009；朱日祥等, 2012).此外, Dickinson(2004)认为怀俄明克拉通及其北侧的Wopmay和南侧的Yavapai-Mazatzal古元古代造山带的减薄/破坏，主要受到中生代太平洋的东向俯冲所控制，这与华北克拉通的减薄/破坏具有一定程度的相似性(吴福元等, 2014).

本文的数值模型并非主要针对华北克拉通而展开，但有的模型可以作为古太平洋俯冲对华北克拉通减薄/破坏效应的简化，我们在上文将上覆大陆岩石层的厚度最高设置到200 km(图 6f).模型试验结果表明, 在海洋板块俯冲作用下, 板块熔体(可近似为交代地幔岩部分熔融产生的熔体)以及由熔体上升而诱发的地幔楔角流场的强烈扰动可以对上覆大陆岩石层在横向空间范围和垂向范围产生一定程度的减薄/破坏作用.当俯冲时间达到40.0 Myr时, 上覆大陆岩石层在横向空间上的减薄范围可以达到500 km, 并且对应着地形的局部抬升(图 9)；模型在垂向上的减薄量可达到70 km(图 6f).然而，从中生代至今，华北克拉通减薄量普遍超过100 km(徐义刚等, 2009), 这说明古太平洋板块的长时间持续俯冲可能会对华北克拉通减薄/破坏产生一定影响(徐义刚等, 2009；朱日祥和郑天愉, 2009；朱日祥等, 2012), 但并非唯一因素.华北克拉通发生减薄/破坏的另一种主要机制可能是太平洋板块俯冲之后的后撤作用.海洋板块的俯冲后撤不仅会引起大规模的软流层物质上涌，同时岩石层也会受到强烈拉张作用，进而导致华北克拉通的减薄/破坏.总之，导致华北克拉通岩石层减薄的因素是多元化的, 比如晚石炭纪古亚洲洋的南向俯冲、蒙古-华北碰撞、古生代特提斯洋的北向俯冲、三叠纪时期华南-华北的陆陆碰撞也都会对华北克拉通周边岩石层的物理化学性质及其稳定性产生一定的影响(翟明国等, 2008).

图 9 上覆大陆岩石层厚度为200 km的模型下的地形演化 颜色代表地形起伏, 如右侧色标.横轴是截取的模型宽度：2500≤x≤3700 km.纵轴是时间轴：0≤t≤40.0 Myr.模型参数：汇聚速率为6 cm·a-1, 洋壳厚度为8 km, 俯冲海洋板块年龄为60 Ma, 上覆大陆岩石层厚度为200 km. Fig. 9 Evolution of the topography of models with 200-km overriding continental lithosphere The axes: horizontally 2500≤x≤3700 km and vertically 0≤t≤40.0 Myr. Convergence velocity is 6 cm·a-1. Thickness of oceanic crust is 8 km. Age of subducting oceanic plate is 60 Ma. Thickness of overriding continental plate is 200 km.

由于模型算法本身的局限性, 本文没有考虑流体活动作用，所以海洋板块在俯冲过程中，没有真实反映板块流体交代上覆地幔楔岩石形成交代地幔岩的部分熔融过程, 以及由此造成的对上覆岩石层地幔的热侵蚀破坏过程(Li et al., 2015；Liu et al., 2017；Yang et al., 2018), 这会对实验结果产生一定影响.此外, 本文重点探讨的是海洋板块俯冲对上覆大陆岩石层的影响及过程，即海洋板块向大陆板块俯冲过程中，海洋板块的洋壳及沉积物发生熔融上升，进而引起地幔热物质上涌，从而导致上覆大陆岩石层减薄.虽然本文没有针对克拉通的真实地质条件来建立模型进行讨论，但模拟结果能为克拉通的减薄/破坏机制提供一些启示.在下一步工作中将考虑俯冲后撤的影响，结合华北克拉通比较真实的地质条件来建立数值模型，从而进一步来探讨华北克拉通减薄/破坏的可能机制.

5 结论

通过建立二维洋-陆俯冲的数值模型, 本文探讨了海洋板块俯冲过程中的地幔楔内的热-熔体扰动所引发的不稳定地幔角流体系对上覆大陆岩石层的影响, 并对上覆大陆岩石层减薄的动力学制约进行了讨论, 得到以下结论：

(1) 海洋板块俯冲过程中的地幔楔熔体对大陆岩石层地幔的热侵蚀效应，以及由熔体上升所诱发的角流区的强烈扰动，都会造成上覆大陆岩石层的减薄效应，这种效应不仅表现在横向上的减薄范围，还表现在垂向上的减薄程度.

(2) 洋-陆俯冲过程中, 随着汇聚速率和洋壳厚度的增加, 上覆大陆岩石层在横向上的减薄范围越大, 在垂向上的减薄程度也越深；随着俯冲海洋板块年龄的升高, 上覆大陆岩石层在横向上的减薄范围会增大, 但在垂向上的减薄程度会减小；而随着上覆大陆岩石层厚度的升高, 其横向减薄范围会减小, 但在垂向上的减薄程度会加深.

致谢  感谢审稿专家对本文提出的宝贵建议，同时感谢刘明启、王振山、王少坡给予本文的建议和帮助.
References
Arcay D, Tric E, Doin M P. 2005. Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics. Physics of the Earth and Planetary Interiors, 149(1-2): 133-153. DOI:10.1016/j.pepi.2004.08.020
Arcay D, Lallemand S, Doin M P. 2008. Back-arc strain in subduction zones: Statistical observations versus numerical modeling. Geochemistry, Geophysics, Geosystems, 9(5): Q05015. DOI:10.1029/2007GC001875
Bittner D, Schmeling H. 1995. Numerical modelling of melting processes and induced diapirism in the lower crust. Geophysical Journal International, 123(1): 59-70. DOI:10.1111/gji.1995.123.issue-1
Brace W F, Kohlstedt D L. 1980. Limits on lithospheric stress imposed by laboratory experiments. Journal of Geophysical Research: Solid Earth, 85(B11): 6248-6252. DOI:10.1029/JB085iB11p06248
Burg J P, Gerya T V. 2005. The role of viscous heating in Barrovian metamorphism of collisional orogens: Thermomechanical models and application to the Lepontine Dome in the Central Alps. Journal of Metamorphic Geology, 23(2): 75-95. DOI:10.1111/jmg.2005.23.issue-2
Chen B, Jahn B M, Arakawa Y, et al. 2004. Petrogenesis of the Mesozoic intrusive complexes from the southern Taihang orogen, north China Craton: Elemental and Sr-Nd-Pb isotopic constraints. Contributions to Mineralogy and Petrology, 148(4): 489-501. DOI:10.1007/s00410-004-0620-0
Chen B, Niu X L, Wang Z Q, et al. 2013. Geochronology, petrology, and geochemistry of the Yaojiazhuang ultramafic-syenitic complex from the North China Craton. Science China Earth Sciences, 56(8): 1294-1307. DOI:10.1007/s11430-013-4603-8
Cheng H D, Huang J S, Fu R S. 2016. Numerical simulation on convective thinning of the cratonic lithosphere. Chinese Journal of Geophysics (in Chinese), 59(4): 1293-1308. DOI:10.6038/cjg20160412
Cheng H D, Huang J S. 2017. Numerical simulation on the detachment of eclogitic lower crust. Scientia Sinica Terrae (in Chinese), 47(1): 82-94. DOI:10.1360/N072016-00019
Clark S R, Stegman D, Müller R D. 2008. Episodicity in back-arc tectonic regimes. Physics of the Earth and Planetary Interiors, 171(1-4): 265-279. DOI:10.1016/j.pepi.2008.04.012
Clauser C, Huenges E. 1995. Thermal conductivity of rocks and minerals. //Rock Physics & Phase Relations: A Handbook of Physical. AGU, 3: 105-126.
Dickinson W R. 2004. Evolution of the North American cordillera. Annual Review of Earth and Planetary Sciences, 32: 13-45. DOI:10.1146/annurev.earth.32.101802.120257
Dvorkin J, Nur A, Mavko G, et al. 1993. Narrow subducting slabs and the origin of backarc basins. Tectonophysics, 227(1-4): 63-79. DOI:10.1016/0040-1951(93)90087-Z
Gao S, Luo T C, Zhang B R, et al. 1998. Chemical composition of the continental crust as revealed by studies in East China. Geochimica et Cosmochimica Acta, 62(11): 1959-1975. DOI:10.1016/S0016-7037(98)00121-5
Gao S, Rudnick R L, Carlson R W, et al. 2002. Re-Os evidence for replacement of ancient mantle lithosphere beneath the North China Craton. Earth and Planetary Science Letters, 198(3-4): 307-322. DOI:10.1016/S0012-821X(02)00489-2
Gerya T V, Yuen D A. 2003a. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. Physics of the Earth and Planetary Interiors, 140(4): 293-318. DOI:10.1016/j.pepi.2003.09.006
Gerya T V, Yuen D A. 2003b. Rayleigh-Taylor instabilities from hydration and melting propel 'cold plumes' at subduction zones. Earth and Planetary Science Letters, 212(1-2): 47-62. DOI:10.1016/S0012-821X(03)00265-6
Gerya T V, Meilick F I. 2011. Geodynamic regimes of subduction under an active margin: Effects of rheological weakening by fluids and melts. Journal of Metamorphic Geology, 29(1): 7-31. DOI:10.1111/jmg.2010.29.issue-1
Gorczyk W, Guillot S, Gerya T V, et al. 2007. Asthenospheric upwelling, oceanic slab retreat, and exhumation of UHP mantle rocks: Insights from Greater Antilles. Geophysical Research Letters, 34(21): L21309. DOI:10.1029/2007GL031059
Green D H, Ringwood A E. 1967. An experimental investigation of the gabbro to eclogite transformation and its petrological applications. Geochimica et Cosmochimica Acta, 31(5): 767-833. DOI:10.1016/S0016-7037(67)80031-0
Han P, Wei D P, Zhang K L, et al. 2016. Lattice-Preferred orientations of olivine in subducting oceanic lithosphere derived from the observed seismic anisotropies in double seismic zones. Earthquake Science, 29(4): 243-258. DOI:10.1007/s11589-016-0160-5
He L J. 2014. Numerical modeling of convective erosion and peridotite-melt interaction in big mantle wedge: Implications for the destruction of the North China Craton. Journal of Geophysical Research: Solid Earth, 119(4): 3662-3677. DOI:10.1002/2013JB010657
Hermann J. 2002. Experimental constraints on phase relations in subducted continental crust. Contributions to Mineralogy and Petrology, 143(2): 219-235. DOI:10.1007/s00410-001-0336-3
Ji S C, Zhao P L. 1993. Flow laws of multiphase rocks calculated from experimental data on the constituent phases. Earth and Planetary Science Letters, 117(1-2): 181-187. DOI:10.1016/0012-821X(93)90125-S
Karig D E. 1971. Origin and development of marginal basins in the Western Pacific. Journal of Geophysical Research, 76(11): 2542-2561. DOI:10.1029/JB076i011p02542
Kirby S H, Kronenberg A K. 1987. Rheology of the lithosphere: selected topics. Reviews of Geophysics, 25(6): 1219-1244. DOI:10.1029/RG025i006p01219
Li Z H, Gerya T V. 2009. Polyphase formation and exhumation of high- to ultrahigh-pressure rocks in continental subduction zone: Numerical modeling and application to the Sulu ultrahigh-pressure terrane in eastern China. Journal of Geophysical Research: Solid Earth, 114(B9): B09406. DOI:10.1029/2008JB005935
Li Z H, Gerya T V, Burg J P. 2010. Influence of tectonic overpressure on P-T paths of HP-UHP rocks in continental collision zones: Thermomechanical modelling. Journal of Metamorphic Geology, 28(3): 227-247. DOI:10.1111/jmg.2010.28.issue-3
Li Z H, Xu Z Q, Gerya T V, et al. 2013. Collision of continental corner from 3-D numerical modeling. Earth and Planetary Science Letters, 380: 98-111. DOI:10.1016/j.epsl.2013.08.034
Li Z H. 2014. A review on the numerical geodynamic modeling of continental subduction, collision and exhumation. Science China Earth Sciences, 57(1): 47-69. DOI:10.1007/s11430-013-4696-0
Li Z H, Liu M Q, Gerya T V. 2015. Material transportation and fluid-melt activity in the subduction channel: Numerical modeling. Science China Earth Sciences, 58(8): 1251-1268. DOI:10.1007/s11430-015-5123-5
Li Z H, Liu M, Gerya T V. 2016. Lithosphere delamination in continental collisional orogens: A systematic numerical study. Journal of Geophysical Research: Solid Earth, 121(7): 5186-5211. DOI:10.1002/2016JB013106
Liu M Q, Li Z H, Yang S H. 2017. Diapir versus along-channel ascent of crustal material during plate convergence: Constrained by the thermal structure of subduction zones. Journal of Asian Earth Sciences, 145: 16-36. DOI:10.1016/j.jseaes.2017.02.036
Menzies M A, Xu Y G. 1998. Geodynamics of the North China craton. //Mantle Dynamics and Plate Interactions in East Asia. Washington, DC: American Geophysical Union, 155-165.
Niu Y L. 2005. Generation and evolution of basaltic magmas: Some basic concepts and a hypothesis for the origin of the Mesozoic-Cenozoic volcanism in eastern China. Geological Journal of China Universities, 11(1): 9-46.
Peacock S M, Wang K L. 1999. Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan. Science, 286(5441): 937-939. DOI:10.1126/science.286.5441.937
Qiao Y C, Guo Z Q, Shi Y L. 2012. Numerical simulation of a possible thinning mechanism of the North China Craton-Gravitational instability delamination induced by lower crust eclogite. Chinese Journal of Geophysics (in Chinese), 55(12): 4249-4256. DOI:10.6038/j.issn.0001-5733.2012.12.036
Qiao Y C, Guo Z Q, Shi Y L. 2013. Thermal convection thinning of the North China Craton: Numerical simulation. Science China Earth Sciences, 56(5): 773-782.
Ranalli G, Murphy D C. 1987. Rheological stratification of the lithosphere. Tectonophysics, 132(4): 281-295.
Ranalli G. 1995. Rheology of the Earth, Deformation and Flow Process in Geophysics and Geodynamics. 2nd ed. London: Chapman & Hall.
Shi Y N, Wei D P, Li Z H, et al. 2018. Subduction mode selection during slab and mantle transition zone interaction: Numerical modeling. Pure and Applied Geophysics, 175(2): 529-548. DOI:10.1007/s00024-017-1762-0
Sobolev S V, Babeyko A Y. 2005. What drives orogeny in the Andes?. Geology, 33(8): 617-620. DOI:10.1130/G21557.1
Stern R J. 2002. Subduction zones. Review of Geophysics, 40(4): 3-1.
Syracuse E M, Van Keken P E, Abers G A. 2010. The global range of subduction zone thermal models. Physics of the Earth and Planetary Interiors, 183(1-2): 73-90. DOI:10.1016/j.pepi.2010.02.004
Turcotte D L, Schubert G. 2002. Geodynamics. Cambridge: Cambridge University Press.
Uyeda S, Kanamori H. 1979. Back-arc opening and the mode of subduction. Journal of Geophysical Research: Solid Earth, 84(B3): 1049-1061. DOI:10.1029/JB084iB03p01049
Van Hunen J, Van Den Berg A P, Vlaar N J. 2000. A thermo-mechanical model of horizontal subduction below an overriding plate. Earth and Planetary Science Letters, 182(2): 157-169. DOI:10.1016/S0012-821X(00)00240-5
Van Hunen J, Van Den Berg A P, Vlaar N J. 2004. Various mechanisms to induce present-day shallow flat subduction and implications for the younger Earth: A numerical parameter study. Physics of the Earth and Planetary Interiors, 146(1-2): 179-194. DOI:10.1016/j.pepi.2003.07.027
Wang Y M, Huang J S, Zhong S J. 2015. Episodic and multistaged gravitational instability of cratonic lithosphere and its implications for reactivation of the North China Craton. Geochemistry, Geophysics, Geosystems, 16(3): 815-833. DOI:10.1002/2014GC005681
Wang Y M, Huang J S, Zhong S J, et al. 2016. Heat flux and topography constraints on thermochemical structure below North China Craton regions and implications for evolution of cratonic lithosphere. Journal of Geophysical Research: Solid Earth, 121(4): 3081-3098. DOI:10.1002/jgrb.v121.4
Windley B F, Maruyama S, Xiao W J. 2010. Delamination/thinning of sub-continental lithosphere mantle under Eastern China: The role of water and multiple subduction. American Journal of Science, 310(10): 1250-1293. DOI:10.2475/10.2010.03
Wu F Y, Ge W C, Sun D Y, et al. 2003. Discussions on the lithospheric thinning in Eastern China. Earth Science Frontiers (in Chinese), 10(3): 51-60.
Wu F Y, Xu Y G, Gao S, et al. 2008. Lithospheric thinning and destruction of the North China Craton. Acta Petrologica Sinica (in Chinese), 24(6): 1145-1174.
Wu F Y, Xu Y G, Zhu R X, et al. 2014. Thinning and destruction of the cratonic lithosphere: A global perspective. Science China Earth Sciences, 57(12): 2878-2890. DOI:10.1007/s11430-014-4995-0
Wu Z P, Hou X B, Li W. 2007. Discussion on Mesozoic basin patterns and evolution in the Eastern North China block. Geotectonica et Metallogenia (in Chinese), 31(4): 385-399.
Xu Y G. 2001. Thermo-tectonic destruction of the Archaean lithospheric keel beneath the Sino-Korean craton in China: Evidence, timing and mechanism. Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy, 26(9-10): 747-757. DOI:10.1016/S1464-1895(01)00124-7
Xu Y G. 2007. Diachronous lithospheric thinning of the North China Craton and formation of the Daxin′anling-Taihangshan gravity lineament. Lithos, 96(1-2): 281-298. DOI:10.1016/j.lithos.2006.09.013
Xu Y G, Blusztajn J, Ma J L, et al. 2008. Late Archean to Early Proterozoic lithospheric mantle beneath the western North China Craton: Sr-Nd-Os isotopes of peridotite xenoliths from Yangyuan and Fansi. Lithos, 102(1-2): 25-42. DOI:10.1016/j.lithos.2007.04.005
Xu Y G, Li H Y, Pang C J, et al. 2009. On the timing and duration of the destruction of the North China Craton. Chinese Science Bulletin, 54: 3379.
Yang S H, Li Z H, Gerya T V, et al. 2018. Dynamics of terrane accretion during seaward continental drifting and oceanic subduction: Numerical modeling and implications for the Jurassic crustal growth of the Lhasa Terrane, Tibet. Tectonophysics, 746: 212-228. DOI:10.1016/j.tecto.2017.07.018
Zhai M G. 2008. Lower crust and lithospheric mantle beneath the North China Craton before the Mesozoic lithospheric disruption. Acta Petrologica Sinica (in Chinese), 24(10): 2185-2204.
Zhang H F, Sun M, Zhou X H, et al. 2002. Mesozoic lithosphere destruction beneath the North China Craton: Evidence from major-, trace-element and Sr-Nd-Pb isotope studies of Fangcheng basalts. Contributions to Mineralogy and Petrology, 144(2): 241-253. DOI:10.1007/s00410-002-0395-0
Zhang K L, Wei D P. 2011. A kinematic thermal model for descending slabs with velocity boundary layers: A case study for the Tonga subducting slab. Acta Geologica Sinica, 85(1): 211-222.
Zhang K L, Wei D P. 2011. On the influence factors of double seismic zones. Chinese Journal of Geophysics (in Chinese), 54(11): 2838-2850. DOI:10.3969/j.issn.0001-5733.2011.11.014
Zhang K L, Wei D P. 2012. Correlation between plate age and layer separation of double seismic zones. Earthquake Science, 25(1): 95-101.
Zheng Y F, Wu F Y. 2009. Growth and reworking of cratonic lithosphere. Chinese Science Bulletin, 54(19): 3347-3353.
Zheng Y F, Chen Y X, Dai L Q, et al. 2015. Developing plate tectonics theory from oceanic subduction zones to collisional orogens. Science China Earth Sciences, 58(7): 1045-1069.
Zhu R X, Zheng T Y. 2009. Destruction geodynamics of the North China Craton and its Paleoproterozoic plate tectonics. Chinese Science Bulletin, 54: 3354.
Zhu R X, Xu Y G, Zhu G, et al. 2012. Destruction of the North China Craton. Science China Earth Sciences, 55(10): 1565-1587. DOI:10.1007/s11430-012-4516-y
陈斌, 牛晓露, 王志强, 等. 2013. 华北克拉通北缘姚家庄过钾质超镁铁岩-正长岩杂岩体的锆石U-Pb年代学、岩石学和地球化学特征. 中国科学:地球科学, 43(7): 1073-1087.
程华冬, 黄金水, 傅容珊. 2016. 克拉通岩石圈对流减薄的数值模拟. 地球物理学报, 59(4): 1293-1308. DOI:10.6038/cjg20160412
程华冬, 黄金水. 2017. 榴辉岩下地壳拆沉的数值模拟. 中国科学:地球科学, 47(1): 82-94.
乔彦超, 郭子祺, 石耀霖. 2012. 数值模拟华北克拉通岩石圈减薄的一种可能机制--下地壳榴辉岩重力失稳引起的拆沉. 地球物理学报, 55(12): 4249-4256. DOI:10.6038/j.issn.0001-5733.2012.12.036
乔彦超, 郭子祺, 石耀霖. 2013. 数值模拟华北克拉通岩石圈热对流侵蚀减薄机制. 中国科学:地球科学, 43(4): 642-652.
吴福元, 葛文春, 孙德有, 等. 2003. 中国东部岩石层减薄研究中的几个问题. 地学前缘, 10(3): 51-60. DOI:10.3321/j.issn:1005-2321.2003.03.004
吴福元, 徐义刚, 高山, 等. 2008. 华北岩石圈减薄与克拉通破坏研究的主要学术争论. 岩石学报, 24(6): 1145-1174.
吴福元, 徐义刚, 朱日祥, 等. 2014. 克拉通岩石圈减薄与破坏. 中国科学:地球科学, 44(11): 2358-2372.
吴智平, 侯旭波, 李伟. 2007. 华北东部地区中生代盆地格局及演化过程探讨. 大地构造与成矿学, 31(4): 385-399. DOI:10.3969/j.issn.1001-1552.2007.04.001
徐义刚, 李洪颜, 庞崇进, 等. 2009. 论华北克拉通破坏的时限. 科学通报, 54(14): 1974-1989.
翟明国. 2008. 华北克拉通中生代破坏前的岩石圈地幔与下地壳. 岩石学报, 24(10): 2185-2204.
张克亮, 魏东平. 2011. 双地震带的影响因素探讨. 地球物理学报, 54(11): 2838-2850. DOI:10.3969/j.issn.0001-5733.2011.11.014
郑永飞, 吴福元. 2009. 克拉通岩石层的生长和再造. 科学通报, 54(14): 1945-1949.
郑永飞, 陈伊翔, 戴立群, 等. 2015. 发展板块构造理论:从洋壳俯冲带到碰撞造山带. 中国科学:地球科学, 45(6): 711-735.
朱日祥, 郑天愉. 2009. 华北克拉通破坏机制与古元古代板块构造体系. 科学通报, 54(14): 1950-1961.
朱日祥, 徐义刚, 朱光, 等. 2012. 华北克拉通破坏. 中国科学:地球科学, 42(8): 1135-1159.