CN117784239B - 过去甲烷水合物稳定带底界面位置的判断方法 - Google Patents

Abstract

本发明所述的过去甲烷水合物稳定带底界面位置的判断方法，引入底层水温度与海平面变化参数，在二维剖面中具体时刻甲烷水合物稳定带底界面深度的求解。包括步骤1、确定模型输入参数；步骤2、确定模拟起始时间；步骤3、获取底层水温度和相对海平面变化；步骤4、建立数值模型；步骤5、求解主控方程；步骤6、检测数值模型；步骤7、温度变化赋值并运行主控方程；步骤8、获取所需时刻t₁的温度矩阵T₁；步骤9、计算第一压力矩阵M1；步骤10、计算第二压力矩阵M2；步骤11、获得所需时刻t₁甲烷水合物稳定带位置；步骤12、结果检查和保存；步骤13、计算t₂至t_n时间范围内的水合物稳定带深度；步骤14、判断过去大范围水合物发生分解的情况。

Description

过去甲烷水合物稳定带底界面位置的判断方法

技术领域

本发明涉及一种为研究过去大范围甲烷水合物发生分解而提供有效计算工具的判断稳定带底界面位置的新方法，属于海洋地质领域。

背景技术

目前地球中存在的天然气水合物（以下简称为水合物），主要是烃类分子（通常是甲烷）和水分子在低温高压下形成的似冰状笼形化合物，其广泛地存在于陆架边缘的海洋环境中。由于外界物理条件的改变容易造成水合物失去稳定性，且会释放出来大量的甲烷气体，进入到海水中的甲烷会造成海洋缺氧和酸化。如果大量甲烷进入到大气当中，还会加剧全球温室效应。考虑到水合物的亚稳定性和巨大存量，普遍认为过去地质历史中水合物在改变全球气候变化的过程中扮演了重要的角色，故而长期受到世界各国政府和学术团体的关注。

现有技术研究是否有区域性大规模水合物发生分解和甲烷泄漏，主要是依靠判断自然界中水合物稳定带底界面深度随着时间的变化过程。天然气水合物稳定带（GasHydrate Stability Zone，简称GHSZ）是满足水合物形成和发育的物理条件的空间，其底界面（Base of Gas Hydrate Stability Zone，简称BGHSZ）可由水合物相平衡曲线和实际温度曲线在温度—深度图版中交汇而确定（如后附图1所示）。即在二维剖面中体现为一条海底下方数十米至数百米的曲线，在深水地区该底界面与海底基本平行，在浅海部分逐渐变浅直至与海底相交。水合物稳定带底界面代表着流体相态的分界，其上方是水合物+水+甲烷，其下方是水+甲烷，由此产生的波阻抗差异可被声学地震剖面捕捉到从而表征为一条高亮反射层，即水合物似海底反射层（BSR）。该底界面在过去地质历史中主要受到温度和压力的影响，在自然界中分别体现为海洋底层水温度和相对海平面高低的变化，会发生向上或向下的位置，其中，前者被认为是区域性大规模水合物发生分解的必要条件，这是因为大量水合物富集在底界面附近，一旦该界面发生向上迁移，会有大量水合物被移出稳定带从而失去稳定性并发生分解。可见，判断水合物稳定带底界面深度随时间变化是研究水合物在过去环境效应的重要步骤。

目前通过水合物稳定带底界面位置深度变化来判断过去是否有大范围甲烷水合物发生分解情况，主要有以下两种方法：

一是，依靠声学地震剖面中发育的残余似BSR特征去反推位置，该特征通常出现在现今BSR的上方或下方，是底界面曾经发育的地质记录，直接指示了过去水合物稳定带底界面的位置。但是此种方式，声学地震剖面中发育的残余似BSR特征存在的情况极为稀少，自然界中现今BSR上下方仅有不到2%的情况能够观察到该特征，多个BSR现象同时出现需要满足苛刻的地质条件才能发育，极少的发育数量远不足以满足研究的需要。特别是浅海区稳定带接近海底的部分，全球研究表明尚未在此处观察到该特征，而此处稳定带深度变化和水合物分解频繁，是重点研究区域。此外，即使观察到该现象，由于缺乏测年分析手段，导致无法确定该底界面对应的过去时间。

二是，使用数值模拟方法结合外界物理条件参数计算该底界面位置，该类方法通常由两种途径实现，分别是使用Tough+Hydrate计算机软件和稳态温度相平衡匹配计算；其中，前者是由美国加州伯克利实验室开发，使用有限元分析建模思路；后者是假设给定时刻模型热传导为零（ΔH=0），计算海底下方沉积物节点温度，之后转化为水合物平衡临界压力后确定该底界面位置深度随时间变化规律。一方面，Tough+Hydrate软件可以计算底界面深度随时间变化规律，但该软件模式对象主要为短时间维度内（小于数年）地质储层尺寸（数百米），在计算自然界中大范围真实情况，通常需要考虑长时间维度（数百至数万年）和大空间尺度（数十至数百公里），导致使用该软件的算力和时间过高（超级计算机上运行十天以上），且研究中会调整输入参数多次模拟，从而进一步加剧这一问题。而且软件设置参数复杂，次要参数无法准确约束时常会导致程序运行中止。另一方面，使用稳态温度相平衡计算是之前研究中具有一定可行性的方法，该方法由于假设了稳态，即忽略了热量传导过程，导致水合物稳定带底界面深度模拟结果可能出现不连续，无法反映该界面在自然中真实的迁移过程，有时甚至会得出和实际情况相反的结果（如后附图2所示），这是因为海平面变化带来的压力影响是瞬时的，而热量传递带来的温度影响是有时间滞后的。因此，此类方法的研究结果并不全面且准确性并不高。

有鉴于此，特提出本专利申请。

发明内容

本发明所述过去甲烷水合物稳定带底界面位置的判断方法，在于解决上述现有技术存在的问题而提出在二维数值模型中引入自然条件底层水温度与海平面变化参数，以期实现在二维剖面中具体时刻甲烷水合物稳定带底界面深度的准确求解，从而达到便于获得随时间连续变化规律与多个时刻的深度进行直观有效对比的设计目的。

为实现上述设计目的，所述过去甲烷水合物稳定带底界面位置的判断方法包括下述步骤：

步骤1、确定模型输入参数；

依据目标二维剖面情景，沿横向和纵向分别建立x轴和z轴，确定数值模型的总长度L和总高度H，获取当前海底x和z坐标位置；

步骤2、确定模拟起始时间；

确定需要模拟的起始时间t₀、结束时间t_end、以及需要计算的稳定带深度时刻，即t₁、t₂……t_n；将需要模拟的地质时间转化成从零开始的顺时时间；

步骤3、获取底层水温度和相对海平面变化；

在模拟时间范围中确定海洋底层水温度和相对海平面变化数值；

步骤4、建立数值模型；

根据模型尺寸确定x轴和z轴方向上的节点数量；

采用二维热传导稳态方程建立数值模型，

，其中，T是沉积物温度，x和z分别是沉积物的长度和距离海底的深度，ρ_b是沉积物密度，C_b是比热容，t是时间；k_x和k_z分别是沉积物在水平和垂直方向上的热导率；

步骤5、求解主控方程；

将主控方程进行离散化处理，在编程软件中完成主控方程求解，将步骤4中的数值模型设置边界条件；

步骤6、检测数值模型；

假设温度变化条件运行程序，例如海底温度随时间呈线性增加，检查海底下方沉积物多个具体节点温度随时间变化的曲线，以验证所述数值模型的可靠性；

步骤7、温度变化赋值并运行主控方程；

在确认数值模型可靠性的前提下，将底层水温度变化曲线按从t₀至t_end时间范围内、以秒为单位进行线性插值，方便微积分中迭代计算；

步骤8、获取所需时刻t₁的温度矩阵T₁；

在数值模型运行完成后，获取所需具体时刻t₁所有模型节点的温度，形成与节点位置一一对应的温度矩阵T₁，将T₁中海水节点数据中z<z_SB的部分删除，此时第一温度矩阵T₁代表对应时刻t₁海底下方沉积物节点的模拟温度；

步骤9、计算第一压力矩阵M1；

根据实验室内水合物相平衡方程，将沉积物节点温度矩阵T1中的各数值转化成压力值，按照节点位置形成二维第一压力矩阵M1，数值单位为兆帕（MPa）；

步骤10、计算第二压力矩阵M2；

按下述公式计算海底下方沉积物内各节点的静水压力，

，其中，ρ_sw是海水密度，g是重力加速度，h是现今海水深度，Δh是t = t₁时相对海平面变化；

各节点压力值根据t₁时刻相对海平面高度进行校正，计算后按照节点位置形成二维第二压力矩阵M2，结果单位为兆帕（MPa）；

步骤11、获得所需时刻t₁甲烷水合物稳定带位置；

将第一压力矩阵M1和第二压力矩阵M2相减得到第三压力矩阵M3，即第三压力矩阵M₃ = M₁-M₂；

在第三压力矩阵M₃中使用插值方法获得数值为零的节点位置，结果即为t = t₁时甲烷水合物底界面的位置，将符合条件的底界面位置连接起来为海底下方一条连续的曲线；

将符合上述条件的位置x、z坐标导出；

步骤12、结果检查和保存；

将步骤11中的计算结果和海底位置进行投图，观察二者空间位置关系，保存计算结果正确的坐标至.txt文件中；

步骤13、计算t₂至t_n时间范围内的水合物稳定带深度；

针对步骤2中的其余时刻t₂至t_n，重复执行上述步骤7至步骤11，分别计算得出各个时刻对应的水合物稳定带深度；

步骤14、判断过去大范围水合物发生分解的情况；

针对目标位置将水合物稳定带底界面深度随时刻t₁至t_n进行投图，获得水合物稳定带深度减小的时间段，即为发生大范围水合物发生分解的时间；判断流程结束。

进一步地，所述的步骤3，海域底层水温度变化数值是通过沉积物底栖有孔虫的氧同位素求得；相对海平面变化数值引用全球海平面变化曲线数据。

进一步地，所述的步骤5，在MATLAB软件中使用ODE45函数进行微积分计算以完成主控方程求解。

进一步地，所述的步骤5，将步骤4中的数值模型设置下述边界条件，

当0 ≤ x ≤ L，且z取值为数值模型的底边界位置时，z轴方向模型节点温度斜率为地温梯度G；

当0 ≤ x ≤ L，且z取值为海底位置时，T_SB海底温度应为所在时刻的温度赋值。

进一步地，所述的步骤6，假设的温度变化条件包括海底温度随时间呈线性增加。

进一步地，所述的步骤9，相平衡方程采用下述多项式拟合公式：

，

其中，P、T为水合物相平衡状态的临界压力和温度。

进一步地，所述的步骤12，以稳定带底界面随海水变浅逐渐靠近海底并最终与之相交为计算结果正确的判断条件。

综上，本申请所述过去甲烷水合物稳定带底界面位置的判断方法具有以下优点与有益效果：

本申请能够为判断过去大范围甲烷水合物发生分解提供有效直接的计算工具，而且计算和判断过程不必借助声学地震数据等地质记录证据，可以直接根据外界海洋底层水温度和相对海平面变化计算出稳定带深度，同时考虑了热量传递带来的热滞后效应，因此在各具体时刻均可获得对应的底界面深度，从而求得其随时间的连续变化规律，方便在实际研究中对多个时刻的深度进行对比；

本申请忽略了难以约束的参数条件，如水合物厚度和饱和度等，相较于考虑此类约束条件的模拟结果，本申请能够在尽可能少影响计算结果的前提下简化了模型建立时的参数赋值，因此判断操作性更强，较为适合研究初期中约束参数较难获得但却需快速判断多种地质物理变化条件下稳定带深度的情形；

本申请使用模型节点计算手段而未使用有限元方法，有限元计算方法针对的是力学-热学情形，而本申请针对的稳定带底界面深度计算不需要考虑力学因素，因此本申请能够大幅度地减少程序运行时间（至少30%）、节约算力，且取消了在超级计算机上运行程序的硬性条件，所生成的结果数据兼容性更好，在多种软件中均可实现可视化，易于推广。

附图说明

图1为温度-深度图版中确定甲烷水合物稳定带底界面深度的示意图；

图2为相同底层水和海平面变化假设条件下水合物稳定带模拟的不同结果对比图；其中左图为热量稳态条件下的模拟示意图；右图为热传导真实情景中的模拟示意图；

图3为本申请所述判断方法建立的数值模型示意图；

图4为本申请所述过去甲烷水合物稳定带底界面位置的判断方法流程图；

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

实施例1，如图3和图4所示，本申请提出的过去甲烷水合物稳定带底界面位置的判断方法，包括下述实施步骤：

步骤1、确定模型输入参数；

依据目标二维剖面情景，沿横向和纵向分别建立x轴和z轴，确定数值模型的总长度L和总高度H，获取当前海底x和z坐标位置；

步骤2、确定模拟起始时间；

确定需要模拟的起始、结束时间，即t₀和t_end、以及需要计算的稳定带深度时刻，即t₁、t₂……t_n；

将需要模拟的地质时间转化成从零开始的顺时时间，例如需要模拟过去2万年至1万年的时间，应将模拟时间设置为0万年至1万年；

随后，将时间单位转换成秒，方便在数值模型中进行阶步计算；

步骤3、获取底层水温度和相对海平面变化；

在模拟时间范围中确定海洋底层水温度和相对海平面变化数值；

具体地，海域底层水温度变化数值，是通过沉积物底栖有孔虫的氧同位素求得；相对海平面变化数值，引用全球海平面变化曲线数据，例如使用Lambeck（2004）文献中的研究结果；

步骤4、建立数值模型；

如图3所示，根据模型尺寸确定x轴和z轴方向上的节点数量；通常x轴和z轴均应有不少于数百个节点，z轴方向上节点间隔应在1至5米之间；

采用二维热传导稳态方程建立数值模型，

，其中，T是沉积物温度，x和z分别是沉积物的长度和距离海底的深度，ρ_b是沉积物密度，C_b是比热容，t是时间；k_x和k_z分别是沉积物在水平（x轴）和垂直（z轴）方向上的热导率；

步骤5、求解主控方程；

将主控方程进行离散化处理，在编程软件中完成主控方程求解，例如在MATLAB软件中使用ODE45函数进行微积分计算；

将步骤4中的数值模型设置下述边界条件，

，其中，G是地温梯度，T是数值模型节点温度，T_SB是海底温度，z_SB(x)是现今海底位置的x、z坐标；

具体地，数值模型设置的边界条件如下：

当0 ≤ x ≤ L，且z取值为数值模型的底边界位置时，z轴方向模型节点温度斜率为地温梯度G；

当0 ≤ x ≤ L，且z取值为海底位置时，T_SB海底温度应为所在时刻的温度赋值；

步骤6、检测数值模型；

假设温度变化条件运行程序，例如海底温度随时间呈线性增加，检查海底下方沉积物多个具体节点温度随时间变化的曲线，以验证所述数值模型的可靠性；

步骤7、温度变化赋值并运行主控方程；

在确认数值模型可靠性的前提下，将底层水温度变化曲线按从t₀至t_end时间范围内、以秒为单位进行线性插值，方便微积分中迭代计算；

即将t₀时刻温度赋予数值模型海水和海底节点作为数值模型的初始条件，运行数值模型；

步骤8、获取所需时刻t₁的温度矩阵T₁；

在数值模型运行完成后，获取所需具体时刻t₁所有模型节点的温度，形成与节点位置一一对应的温度矩阵T₁，将T₁中海水节点数据中z<z_SB的部分删除，此时第一温度矩阵T₁代表对应时刻t₁海底下方沉积物节点的模拟温度；

步骤9、计算第一压力矩阵M1；

根据实验室内水合物相平衡方程，将沉积物节点温度矩阵T1中的各数值转化成压力值，按照节点位置形成二维第一压力矩阵M1，数值单位为兆帕（MPa）；

上述相平衡方程可采用如Moridis（2008）在实验室求得的多项式拟合公式：

，其中，P、T为水合物相平衡状态的临界压力和温度；

步骤10、计算第二压力矩阵M2；

按下述公式计算海底下方（即沉积物内各节点）静水压力：

,

其中，ρ_sw是海水密度，g是重力加速度，h是现今海水深度，Δh是t = t₁时相对海平面变化；

各节点压力值根据t₁时刻相对海平面高度进行校正，例如现今水深1000米，在t =t₁时相对海平面为-100米，此处水深应为1000-100 = 900米；

计算后按照节点位置形成二维第二压力矩阵M2，结果单位为兆帕（MPa）；

上述第一压力矩阵M1和第二压力矩阵M2对应的节点位置应该相同、且具有相同的行数和列数；

步骤11、获得所需时刻t₁甲烷水合物稳定带位置；

将第一压力矩阵M1和第二压力矩阵M2相减得到第三压力矩阵M3，即第三压力矩阵M₃ = M₁-M₂；

在第三压力矩阵M₃中使用插值方法获得数值为零的节点位置，结果即为t = t₁时甲烷水合物底界面的位置，将符合条件的底界面位置连接起来为海底下方一条连续的曲线；

将符合上述条件的位置x、z坐标导出；

步骤12、结果检查和保存；

将步骤11中的计算结果和海底位置进行投图，观察二者空间位置关系；

以稳定带底界面随海水变浅逐渐靠近海底并最终与之相交为计算结果正确的判断条件，保存计算结果正确的坐标至.txt文件中；

步骤13、计算t₂至t_n时间范围内的水合物稳定带深度；

针对步骤2中的其余时刻t₂至t_n，重复执行上述步骤7至步骤11，分别计算得出各个时刻对应的水合物稳定带深度；

步骤14、判断过去大范围水合物发生分解的情况；

针对目标位置将水合物稳定带底界面深度随时刻t₁至t_n进行投图，获得水合物稳定带深度减小的时间段，即为发生大范围水合物发生分解的时间；判断流程结束。

如上所述，结合附图和描述给出的方案内容，可以衍生出类似的技术方案，均仍属于本发明技术方案的权利范围。

Claims (6)

1.一种过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：包括下述实施步

骤，

步骤1、确定模型输入参数；

依据目标二维剖面情景，沿横向和纵向分别建立x轴和z轴，确定数值模型的总长度L和总高度H，获取当前海底x和z坐标位置；

步骤2、确定模拟起始时间；

确定需要模拟的起始时间t₀、结束时间t_end、以及需要计算的稳定带深度时刻，即t₁、t₂……t_n；将需要模拟的地质时间转化成从零开始的顺时时间；

步骤3、获取底层水温度和相对海平面变化；

在模拟时间范围中确定海洋底层水温度和相对海平面变化数值；

步骤4、建立数值模型；

根据模型尺寸确定x轴和z轴方向上的节点数量；

采用二维热传导稳态方程建立数值模型，

;

其中，T是沉积物温度，x、z分别是x轴、z轴坐标值，ρ _b 是沉积物密度，C _b 是比热容，t是时间；k _x 和k _z 分别是沉积物在水平和垂直方向上的热导率；

步骤5、求解主控方程；

将主控方程进行离散化处理，在编程软件中完成主控方程求解，将步骤4中的数值模型设置边界条件；

步骤6、检测数值模型；

假设海底温度随时间呈线性增加，运行程序，检查海底下方沉积物多个具体节点温度随时间变化的曲线，以验证所述数值模型的可靠性；

步骤7、温度变化赋值并运行主控方程；

在确认数值模型可靠性的前提下，将底层水温度变化曲线按从t₀至t_end时间范围内、以秒为单位进行线性插值，方便微积分中迭代计算；

步骤8、获取所需时刻t₁的温度矩阵T₁；

在数值模型运行完成后，获取所需具体时刻t₁所有模型节点的温度，形成与节点位置一一对应的温度矩阵T₁，将T₁中海水节点数据中z<z_SB的部分删除，其中z_SB是现今海底位置的z轴坐标值，此时第一温度矩阵T₁代表对应时刻t₁海底下方沉积物节点的模拟温度；

步骤9、计算第一压力矩阵M1；

根据实验室内水合物相平衡方程，将沉积物节点温度矩阵T1中的各数值转化成压力值，按照节点位置形成二维第一压力矩阵M1，数值单位为兆帕（MPa）；

步骤10、计算第二压力矩阵M2；

按下述公式计算海底下方沉积物内各节点的静水压力，;

其中，ρ_sw是海水密度，g是重力加速度，h是现今海水深度，Δh是t = t₁时相对海平面变化；

各节点压力值根据t₁时刻相对海平面高度进行校正，计算后按照节点位置形成二维第二压力矩阵M2，结果单位为兆帕（MPa）；

步骤11、获得所需时刻t₁甲烷水合物稳定带位置；

将第一压力矩阵M1和第二压力矩阵M2相减得到第三压力矩阵M3，即第三压力矩阵M₃ =M₁-M₂；

在第三压力矩阵M₃中使用插值方法获得数值为零的节点位置，结果即为t = t₁时甲烷水合物底界面的位置，将符合条件的底界面位置连接起来为海底下方一条连续的曲线；

将符合上述条件的位置x、z坐标导出；

步骤12、结果检查和保存；

将步骤11中的计算结果和海底位置进行投图，观察二者空间位置关系，保存计算结果正确的坐标至.txt文件中；

步骤13、计算t₂至t_n时间范围内的水合物稳定带深度；

针对步骤2中的其余时刻t₂至t_n，重复执行上述步骤8至步骤11，分别计算得出各个时刻对应的水合物稳定带深度；

步骤14、判断过去大范围水合物发生分解的情况；

针对目标位置将水合物稳定带底界面深度随时刻t₁至t_n进行投图，获得水合物稳定带深度减小的时间段，即为发生大范围水合物发生分解的时间；判断流程结束。

2.根据权利要求1所述的过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：所述的步骤3，海域底层水温度变化数值是通过沉积物底栖有孔虫的氧同位素求得；相对海平面变化数值引用全球海平面变化曲线数据。

3.根据权利要求1所述的过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：所述的步骤5，在MATLAB软件中使用ODE45函数进行微积分计算以完成主控方程求解。

4.根据权利要求1或2所述的过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：所述的步骤5将步骤4中的数值模型设置下述边界条件，

当0 ≤x ≤ L，且z取值为数值模型的底边界位置时，z轴方向模型节点温度斜率为地温梯度G；

当0 ≤x ≤ L，且z取值为海底位置时，海底温度T_SB应为所在时刻的温度赋值。

5.根据权利要求1所述的过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：所述的步骤9，相平衡方程采用下述多项式拟合公式，

；

其中，P、Te分别为水合物相平衡状态的临界压力和温度。

6.根据权利要求1所述的过去甲烷水合物稳定带底界面位置的判断方法，其特征在于：所述的步骤12，以稳定带底界面随海水变浅逐渐靠近海底并最终与之相交为计算结果正确的判断条件。