一种板块俯冲热结构参数的数值模拟分析方法

摘要

本发明公开了一种板块俯冲热结构参数的数值模拟分析方法，基于纳维斯托克斯的质量、动量和能量守恒方程，板块从海沟开始俯冲至目前形状和深度并达到稳定状态后，调整修改模型参数以使地表热流量计算值达到观测值，据此得出最优模型的结果。该方法在已知板块深度和俯冲速度数据的基础上，可以计算和表现俯冲带的三维温度场、板块含水率、脱水率和岩石相变分布。该计算结果可以进一步用于地球深部温度场、岩石相变折返过程、火山岩浆的地下流体场以及变质岩矿产勘查和开发的研究。

说明书

技术领域

本发明涉及板块俯冲地球动力学热数值模拟技术领域，特别涉及一种板块俯冲热结构参数的数值模拟分析方法。

背景技术

板块俯冲带是地球表面和深部进行物质和能量循环的重要环节，也是地震能量释放的主要地带。板块俯冲是地幔对流在地表的表现形式，俯冲的陆地或海洋板块中的低密度陆壳物质经过折返成超高压变质带，或者在上覆地幔的高温下发育成熔流体到达陆壳基底形成岩浆作用，形成了各具特色的深部构造和物质交代(许志琴等，2018)。

目前国际上板块俯冲热数值模型主要为二维结构，三维结构较少；板块温度构造和岩石相变计算结果的三维化仿真更为欠缺。核心瓶颈是板块三维形状与三维俯冲速度的代码化和模型化。

如何更准确地计算板块俯冲热结构，同时还实现基于三维模型结果的仿真表达，成为同行从业人员亟待解决的问题。

发明内容

鉴于上述问题，提出了本发明以便提供一种克服上述问题或者至少部分地解决上述问题的板块俯冲热结构参数的数值模拟分析方法。

本发明实施例提供一种板块俯冲热结构参数的数值模拟分析方法，包括：

根据俯冲年龄和板块自由冷却模型计算板块边界温度初始条件；

根据板块深度和俯冲速度，设定板块内部物质的区域化和速度场，并给定模型各区域的相关参数，所述相关参数包括：密度、粘度、热传导率、比热容、热膨胀率和热辐射率；

根据所述边界温度初始条件和所述相关参数，利用质量、动量和能量守恒的纳维斯托克斯方程，迭代演算每个俯冲步骤的质量、温度、速度和能量场；

当俯冲过程完成且达到温度稳定状态后，得出最终质量场、温度场、速度场和能量场；

对比地表热流观测值，调整所述给定模型各区域的相关参数，使温度场符合观测值：

根据所述温度场和速度场，计算板块含水率、脱水率和岩石相变分布结果；

根据计算出的所述含水率、脱水率、岩石相变分布结果以及所述温度场，进行三维绘图。

进一步地，根据所述边界温度初始条件和所述相关参数，利用质量、动量和能量守恒的纳维斯托克斯方程，迭代演算每个俯冲步骤的质量、温度、速度和能量场；包括：

构建三维数值模型：

质量守恒：

动量守恒：

能量守恒：

温度初始条件根据自由冷却半空间下上覆岩石圈和俯冲岩石圈的温度计算得出；

上覆岩石圈的初始温度：

俯冲岩石圈的初始温度：

数值模型的边界条件为顶部T＝0和前左右三个侧面的上覆板块和软流圈为自由滑动边界，底部边界为渗透性边界，后侧边界为曲线会聚边界；

自板块边界温度初始条件开始，后侧边界的会聚速率为实际俯冲速度；温度边界条件，模型顶部为固定温度0℃，四周的侧面边界条件为水平方向温度梯度为零；底部边界采用的为外部边界固定温度条件1600K；

实现迭代演算每个俯冲步骤的质量、温度、速度和能量场。

进一步地，所述温度稳定状态为温度误差低于预设阈值；所述预设阈值为0.1％～10％。

本发明实施例提供的上述技术方案的有益效果至少包括：

本发明实施例提供的一种板块俯冲热结构参数的数值模拟分析方法，基于纳维斯托克斯的质量、动量和能量守恒方程，板块从海沟开始俯冲至目前形状和深度并达到稳定状态后，调整修改模型参数以使地表热流量计算值达到观测值，据此得出最优模型的结果。该方法在已知板块深度和俯冲速度数据的基础上，可以计算和表现俯冲带的三维温度场、板块含水率、脱水率和岩石相变分布。该计算结果可以进一步用于地球深部温度场、岩石相变折返过程、火山岩浆的地下流体场以及变质岩矿产勘查和开发的研究。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1为本发明实施例提供的板块俯冲热结构参数的数值模拟分析方法流程图；

图2为本发明实施例提供的实施例的原理图；

图3为菲律宾海板块在日本西南的旋转示意图。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

参照图1所示，一种板块俯冲热结构参数的数值模拟分析方法，包括：

S10、根据俯冲年龄和板块自由冷却模型计算板块边界温度初始条件；

S20、根据板块深度和俯冲速度，设定板块内部物质的区域化和速度场，并给定模型各区域的相关参数，所述相关参数包括：密度、粘度、热传导率、比热容、热膨胀率和热辐射率；

S30、根据所述边界温度初始条件和所述相关参数，利用质量、动量和能量守恒的纳维斯托克斯方程，迭代演算每个俯冲步骤的质量、温度、速度和能量场；

S40、当俯冲过程完成且达到温度稳定状态后，得出最终质量场、温度场、速度场和能量场；温度稳定状态为温度误差低于预设阈值；该预设阈值为0.1％～10％，可根据具体板块的温度场数据选择合适的阈值。

S50、对比地表热流观测值，调整所述给定模型各区域的相关参数，使温度场符合观测值：

S60、根据所述温度场和速度场，计算板块含水率、脱水率和岩石相变分布结果；

S70、根据计算出的所述含水率、脱水率、岩石相变分布结果以及所述温度场，进行三维绘图。

本发明实施例，基于地震波层析成像数据和板块三维俯冲速度的三维数值模型已经解决了这一技术难题，既能更准确地计算板块俯冲热结构，同时还实现了模型结果的仿真表达。在已知板块深度(经度、纬度和深度)和俯冲速度(经度、纬度和速度)数据的基础上，可以计算和表现俯冲带的三维温度场、板块含水率、脱水率和岩石相变分布。该计算结果可以进一步用于地球深部温度场、岩石相变折返过程、火山岩浆的地下流体场以及变质岩矿产勘查和开发的研究。

该方法基于纳维斯托克斯的质量、动量和能量守恒方程，板块从海沟开始俯冲至目前形状和深度并达到稳定状态后，调整修改模型参数以使地表热流量计算值达到观测值，据此得出最优模型的结果。

下面为了便于说明本发明的技术方案，以青藏高原底部印度板块俯冲为例，进行详细的阐述。

构建青藏高原底部印度板块俯冲的三维数值模型并计算俯冲板块的温度结构、高温变相岩分布特征和地幔流场的动力学特点，参照图2所示。

A.模型代码研发

迄今为止热结构模型较少考虑板块形状在海沟方向上的不规则变化，以及忽略俯冲角度复杂时存在的计算误差。如：当俯冲速度方向与二维计算断面不平行时，速度的第三维分量会带来俯冲增量或减量却未予以考虑；而当板块形状呈现出第三维上的弯曲时，计算断面的相邻物理场会对断面内产生影响也被忽略不计。Wada et al.(2015)和Rosas etal.(2015)的热结构模型也未反映俯冲速度的空间三维变化。本发明实施例通过之前对数值模型代码Stag3d(Tackley and Xie,2003)的延伸开发研究，导入俯冲太平洋板块三维形状和速度如图3所示。

B.板块形状和速度数据收集和导入

导入青藏高原印度板块俯冲界面的板块年龄(全球海底年龄数据库，Muller etal.2008)、板块形态(Slab2，Hayes et al.,2018，Science)、板块相对运动速度(Morvel，DeMets et al.,2010)等构建俯冲中的印度板块岩石圈域，给与合理密度和粘度设定。每个计算步时板块长度按俯冲速度不断延伸，直到穿越大陆板块侧的模型边界，并使温度场达到稳定状态。

C.模型参数选择和配置

结合提斯洋闭合和洋陆转换的俯冲过程，模型求解时间可设置为50-65个百万年或者其他前人研究成果。遵循不可压缩物质质量守恒定律、斯托克斯方程和动力学能量守恒定律等(e.g.,Gerya,2010；Yoshioka andMurakami,2007；Li,2014；Li et al.,2012,2013,2014,2015,2016,2017,2019；李忠海，2014；李忠海等，2015；Chen,2014；Chen etal.,2013,2014,2016,2017)以及黏塑性流变学性质本构关系和粘滞系数(Ranalli,1995等)，本发明可以构建三维数值模型：

质量守恒：

动量守恒：

能量守恒：

温度初始条件根据自由冷却半空间下上覆岩石圈和俯冲岩石圈的温度计算得出；

上覆岩石圈的初始温度：

俯冲岩石圈的初始温度：

数值模型的边界条件为顶部T＝0和前左右三个侧面的上覆板块和软流圈为自由滑动边界，底部边界为渗透性边界，后侧边界为曲线会聚边界；

自板块边界温度初始条件开始，后侧边界的会聚速率为实际俯冲速度；温度边界条件，模型顶部为固定温度0℃，四周的侧面边界条件为水平方向温度梯度为零(零热流量)；底部边界采用的为外部边界固定温度条件1600K；可据模型演化而动态变化。大陆板块和海洋板块的差异性如冷却年龄、厚度、岩石成分和流变性等将在本研究中作为一个重要方面考虑并探讨合理模型参数。

D.模型约束条件

作为热结构模型的基本约束条件，通过BSR等手段得到的全球地表热流量观测值数据库和居里点深度推导的全球值。此外，层析成像结果对于俯冲板片脱水位置和可能上升岩浆路径等有参考作用(Nabelek et al.,2009；Zhang et al.,2018等)；大地测量则提供了地表应变场并据以推测地壳变形和壳幔对流耦合等(Klemperer,S.L.，2006等)，地磁学结果可以推测俯冲板片的空间年龄分布，地质学研究可提供俯冲岩石相变后折返上升到地表的证据，从而估算榴辉岩和橄榄岩等变质岩相分布。青藏高原后碰撞岩浆作用形成钾质、超钾质、火山岩、钾质埃达克岩、淡色花岗岩与钾质钙碱性花岗岩为特征，集中发育于冈底斯构造岩浆带和西藏南部地区。

E.模型运行计算产出

在以上地球物理学参数选择和约束条件下，经过反复调试和最优化配置，在高性能平行设备上进行数值模拟计算，可以得到印度板块三维温度场以及高温变相岩在青藏高原底部的分布规律。产出包括：印度板块俯冲三维温度场、垂直剖面俯冲温度结构、水平剖面俯冲温度结构、印度板块三维温度场数据库、印度板块俯冲高温高压变质岩相三维分布、垂直剖面高温高压变质岩相分布、水平剖面高温高压变质岩相分布、印度板块高温高压变质岩相三维分布数据库以及地幔流三维分布等。

由于前人的二维俯冲热数值模型的精度低且不同模型间结果差别大，温度误差多达200-600摄氏度，与地表热流观测值出入较多，使得俯冲带深部热结构研究陷入停滞。其中未能导入三维板块形状和速度场同时缺乏三维化表达是温度误差的主要原因之一。本发明实施例提供的板块俯冲热结构参数的分析方法，可以克服俯冲数值模型不能导入三维板块形状和三维俯冲速度的瓶颈，使得客观正演板块三维俯冲热结构成为可能，其计算的三维板块温度场、板块含水率、脱水率和岩石相变分布比二维模型精度明显提高，温度误差降低到100摄氏度以下。以三维温度结构为基础计算的地下流体和岩石相变分布具有更高精度和可靠性，能够为变质岩矿产勘查和开发提供参考。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。