构造掀斜主导的断陷湖盆缓坡层序“源—汇”正演模拟定量研究

2023-01-13谈明轩朱筱敏张自力孙浩南

沉积学报 2022年6期

关键词：层序源区样式

谈明轩，朱筱敏，张自力，孙浩南

1.河海大学海洋学院，南京 210098

2.中国石油大学（北京）地球科学学院，北京 102249

3.中国石油天然气公司西南油气田分公司勘探开发研究院，成都 610051

0 引言

构造掀斜作用是断陷盆地构造沉降和隆升的重要机制之一[1]。掀斜作用通常包括正向掀斜（向陆一侧相对抬升、向盆一侧相对沉降）和反向掀斜（向陆一侧相对沉降、向盆一侧相对抬升）两类[1￣2]。其中，正向掀斜作用主要发育于断陷盆地缓坡带上盘，也是本文的主要研究对象。历史地震监测数据证实现代活动断陷盆地存在短期时间尺度的幕式掀斜作用[3￣4]。构造掀斜作用不仅控制了断陷盆地的水系样式和剥蚀作用[5￣6]，而且对湖平面变化具有极为重要的影响[1,7￣8]。前人理论研究表明，在构造主控的断陷湖盆中，气候变化对湖平面影响相对微弱，湖平面变化与构造掀斜枢纽点的相对位置密切相关[9]。枢纽点位于湖平面之上造成缓坡带相对湖平面上升；枢纽点位于湖平面之下时造成缓坡带相对湖平面下降（图1）。然而目前大多数研究停留于几何学模型和地层定性分析层面[7,10]，缺乏相关实例的验证，其理论可靠性仍有待进一步检验。

图1 构造掀斜作用下恒定绝对湖平面的缓坡带滨线迁移模式（据文献[9]修改）（a）枢纽点位于湖平面之上；（b）枢纽点位于湖平面之下Fig.1 Schematic model showing shoreline migration influenced by tilting with a fixed absolute lake level in the hanging￣wall dipslope (modified from reference [9])(a) pivotal point above lake level; (b) pivotal point below lake level

可容空间和沉积物供给是断陷盆地充填特征及其地层建造样式的重要控制因素[5,11]。由于构造掀斜作用能够直接影响源区剥蚀量和盆地沉降量，因此是造成掀斜型缓坡层序多样性的根本原因之一。然而断陷盆地缓坡带在掀斜作用下遭受广泛的暴露剥蚀，地层发育通常并不完整，为缓坡带层序地层学研究带来多解性。地层数值模拟的优势在于其能够通过改变边界条件建立不同类型层序地层样式，能够正演地层沉积充填过程，是定量层序地层学分析的一种有效手段。地貌演化与“源—汇”系统近年来是地质学热点议题，是构造地貌学、盆地动力学与沉积学交叉研究内容[12￣13]。本次数值模拟研究主要基于“源—汇”沉积通量守恒原理的三维地层正演模拟方法[14￣16]，分析构造掀斜作用相关的剥蚀与沉积作用及其对断陷湖盆缓坡带湖平面变化、层序地层充填的直接影响，以“源—汇”系统研究手段探讨掀斜型缓坡层序样式的主要驱动机制。

1 研究方法

本次模拟所应用的Badlands软件是Python架构下以河流作用和坡面过程模拟侵蚀作用、泥沙运移和沉积过程的地貌演化开源程序[17￣20]。与其他基于扩散模型、几何模型和水力学模型的地层正演模拟软件相比[21￣22]，该地层正演模拟代码充分考虑了泥沙沉积作用与流域侵蚀作用的内在关系[18]。因此，该代码所模拟的层序地层发育样式和沉积演化特征不仅能够直接响应于构造掀斜作用，而且能将沉积物通量和流域地貌演化直接纳入“源—汇”系统研究。

1.1 剥蚀、搬运及沉积作用理论模型

该正演地层模型主要基于受“分离”限制（“D-L”模型）的水力侵蚀模型定量表征河流下切作用。岩石剥蚀速率E是源区面积（A）、净降雨量（P）、坡度（S）的函数，即：

式中：剥蚀系数ε（yr-1）是气候、底床糙度及岩石强度的函数。无量纲参数m和n默认值分别为0.5和0.1[17￣18]。由此可见，沉积物供给主要受控于源区的汇水面积、地形高差及山体陡度等地貌参数。

此外，该模型通过利用线性或非线性扩散模拟坡面过程[19]，即：

式中：z为地形高差，扩散系数κ（m2/yr）与源区岩性及降雨量有关[23￣24]。地表、水下及河流的扩散系数均有所差异。其中，冲积平原区临界沉积坡度为0.001 m/m。当沉积物进入水下，则沿多方向水流进行扩散搬运。

此外，Badlands软件后处理工具包能够自动提取正演模型中不同时期滨线位置的新增可容空间深度δA（m）和新增沉积地层厚度δS（m）[17]，为本次定量模拟地层学研究提供了新的方法。

1.2 初始盆地模型的参数设置

本研究将断陷盆地缓坡带简化为平缓的单斜地形，其初始地形设置为100 km（长度）×100 km（宽度）、地表网格间距设置为1 km。该单斜初始地形最大相对高差为200 m，最大水深为-400 m，平均坡降为0.006 m/m（图2）。为更直观地揭示地层垂向叠置特征，该正演数值模型的运行时间设置为2 Myr（三级层序时间尺度），单层沉积时间为10 kyr（准层序时间尺度），累计沉积了200层地层单元。

图2 简化的缓坡带初始三维古地形Fig.2 Simplified initial three￣dimensional paleo￣topography of the hanging￣wall dipslope

1.3 控制变量的参数设置

为突出构造掀斜对层序构型的控制作用，本次模拟研究所设置的模型均简化为恒定的绝对湖平面及模拟研究区均一的降雨量（即2 000 mm/yr），以达到控制单一变量的目的。在湖盆中，绝对湖平面是由相对与一个固定基准面（例如地心）所定义的，相对湖平面则是湖盆底面距湖面的高度[25]。绝对湖平面在自然界实例较少，但也并不是完全不可能实现的。闭流湖盆可通过降雨量和蒸发量的平衡所实现，而敞流湖盆则可通过构造稳定的出水口保持恒定的绝对湖平面。在本研究简化的构造模型中，构造掀斜的枢纽位置固定不变，位移量与距离呈线性关系。基于枢纽点和绝对湖平面的相对位置，本研究共建立两个构造模型：构造模型A（最大隆升量和沉降量分别为100 m和-400 m）（图3a）；构造模型B（最大隆升量和沉降量分别为300 m和-300 m）（图3b）。在模拟运行过程中，研究区的抬升量和沉降量呈线性的变化特征。

图3 两种掀斜构造模型的剖面示意及三维立体特征（a）枢纽点位于湖平面之上；（b）枢纽点位于湖平面之下Fig.3 Schematic profiles of two tectonic tilting models and three￣dimensional characteristics(a) pivotal point above lake level; (b) pivotal point below lake level

为系统探讨层序地层样式多样性可能，本次模拟研究还额外考虑沉积物供给对掀斜缓坡带地层的影响，因此共建立了4种地层正演模型（表1）。其中，模型1、2和模型3、4分别为高、低沉积物供给背景，主要通过源区岩性差异所造成的高、低风化强度进行调节。本研究主要参照希腊Corinth湾裂谷盆地母岩区多类型基岩剥蚀能力的相关研究[6]，将剥蚀系数ε值设置为1×10-6yr-1代表源区主体以硬岩性组成，岩石扛风化能力较强，同等气候条件下沉积物供给能力较弱；ε值设置为1×10-5yr-1代表源区为混合岩性组成，岩石扛风化能力较弱，同等气候条件下沉积物供给能力较强[26￣27]。其他一系列（地表、水下和河流）扩散系数的取值主要参考了前人三角洲数值模拟研究中所选取的数量级，并在各模型设置中保持不变（表1），便于在数值模拟中通过控制相关变量进行对照分析。通过公式（1）进行运算，模型3、4所产生的总沉积物通量仅为模型1、2计算结果的10%左右。由于该正演模拟程序的水下扩散模型较为简单，并没有分别考虑砂泥扩散特征，因此多种沉积相类型主要依据水深进行划分[17￣18,28]。

表1 地层正演模拟关键参数基本设置Table 1 Key parameters of stratigraphic forward model runs

2 研究结果

2.1 模型1结果

模型1中，湖岸线向盆地边缘后退了12 km，表明掀斜作用造成缓坡带的相对湖平面上升（图4a）。湖岸线在初始0.75 Myr向盆缘的平均迁移速率为10.7 km/Myr，而后其平均迁移速率降至4.8 km/Myr（图4b）。该模型中源区地貌最终呈夷平特征（图4c），这表明虽然掀斜作用可造成源区不断抬升，但强烈剥蚀的能力同样可导致源区快速夷平。尽管相对湖平面在不断上升，但该模型中的地层发育样式类似于一个完整的被动陆缘I型三级层序构型特征（图4b，c）。虽然湖平面并没有经历经典层序的海平面变化特征，该层序早期和晚期发育的进积型三角洲与其层序模式中的低位前积楔和高位期进积三角洲相似程度较高（图4b）。

图4 模型1地层剖面、平面模型及相关沉积物通量变化特征（a）剖面A—A’地层叠置样式；（b）Wheeler域剖面；（c）模拟区剥蚀/沉积厚度平面展布；（d）模拟区沉积物通量直方图（5万年为间隔）Fig.4 Stratigraphic profile and planar model of model run 1 and variation of sediment flux(a) stratigraphic stacking pattern of profile A￣A’; (b) Wheeler diagram; (c) distribution of denudation and sedimentation; (d) histogram of sediment flux per 50 kyr

本层序模型中类“低位体系域”（LST）持续时间约为0.3 Myr（图4c），主要发育进积型地层叠置样式。此时期研究区总沉积物通量达到峰值（约为600 km3/Myr）（图4d），三角洲进积型叠置样式与滨线迁移方向相反，整体符合低位楔的发育特征。类“湖侵体系域”（TST）持续时间达0.45 Myr（图4c），然而此时期研究区总沉积物通量骤减为370 km3/Myr（图4d）。类“高位体系域”（HST）则主要表现为加积型地层叠置样式，该时期研究区总沉积物通量降至约200 km3/Myr（图4d）。

2.2 模型2结果

模型2中，湖岸线向盆地方向推进了18 km，表明掀斜作用造成缓坡带相对湖平面下降（图5a）。湖岸线向盆地方向的平均迁移速率高达45 km/Myr（发育于初始0.4 Myr），同时伴随着广泛的陆上暴露剥蚀（图5b），随后1.6 Myr湖岸线迁移基本停止。该模型中源区地貌最终呈夷平的准平原化特征（图5c），表明源区基本已剥蚀殆尽。由于研究区发育较为广泛的湖退作用，因此该层序中地层样式主要发育类“下降期体系域”（FSST）和类“低位体系域”（LST）（图5b）。

本层序模型中类“下降期体系域”（FSST）主要发育于初始0.4 Myr，强制海退现象明显，造成了显著的陆上暴露剥蚀，甚至形成类似于层序界面的特征（图5a）。该时期滨线迁移轨迹为下降型，同时研究区总沉积物通量在短时间内达到峰值（不低于470.4 km3/Myr）。类“低位体系域”（LST）上超于类“下降期体系域”（FSST）之上，事实上已经形成新的层序界面（图5a）。由此表明，掀斜作用是构造主控型断陷湖盆层序界面成因之一[29]。由于相对湖平面持续下降，滨岸平原也造成持续遭受剥蚀（图5b）。该体系域主要表现为降积、加积—进积型的地层叠置样式，同期该模拟研究区的总沉积物通量进一步降低（图5c，d）。

图5 模型2地层剖面、平面模型及相关沉积物通量变化特征（a）剖面B—B’地层叠置样式；（b）Wheeler域剖面；（c）模拟区剥蚀/沉积厚度平面展布；（d）模拟区沉积物通量直方图（5万年为间隔）Fig.5 Stratigraphic profile and planar model of model run 2 and variation of sediment flux(a) stratigraphic stacking pattern of profile B￣B’; (b) Wheeler diagram; (c) distribution of denudation and sedimentation; (d) histogram of sediment flux per 50 kyr

2.3 模型3结果

尽管模型3中的相对湖平面不断上升，缓坡带三角洲却整体保持相对稳定（图6a，b）。在该模型中，湖岸线向盆地边缘后退了11 km，其平均迁移速率达5.5 km/Myr（图6b）。该模型中源区地貌以线性和树枝型沟谷发育为特征（图6c），表明源区仍处于强烈剥蚀阶段。尽管使用了相同的构造模型（表1），低沉积物供给的模型3与高沉积物供给下的模型1在地层样式上已存在较大的差异（图4a、图6a）。该地层正演模型并不具有经典层序三分构型特征，垂向上主要表现为退积—加积型地层样式，仅与高位体系域（HST）更为类似。与高沉积供给条件的模型1对比，该模型中研究区的总沉积物通量整体不超过60 km3/Myr（仅为模型1总沉积通量的近十分之一），其变化趋势却全然不同（图6d）。模型3中模拟研究区的沉积物通量变化在初始阶段达到峰值56.7 km3/Myr，而后缓慢降至55.6 km3/Myr，并最终缓慢上升至56.6 km3/Myr（图6d）。

图6 模型3地层剖面、平面模型及相关沉积物通量变化特征（a）剖面C—C’地层叠置样式；（b）Wheeler域剖面；（c）模拟区剥蚀/沉积厚度平面展布；（d）模拟区沉积物通量直方图（5万年为间隔）；分水岭的形成实为模拟的边界效应所致Fig.6 Stratigraphic profile and planar model of model run 3 and variation of sediment flux(a) stratigraphic stacking pattern of profile C￣C’; (b) Wheeler diagram; (c) distribution of denudation and sedimentation; (d) histogram of sediment flux per 50 kyr. The formation of drainage divide is caused by the modeling boundary effect

2.4 模型4结果

模型4中源区地貌线性和树枝型沟谷发育（图7c），表明源区仍处于强烈剥蚀阶段。该模型中湖岸线最终向盆地方向推进了模型2中同样的距离，表明相对湖平面下降的幅度是一致的（图5b、图7b）。然而其湖岸线迁移特征却迥然相异，整体具有长期、缓慢的特征。模型4在1.6 Myr中平均迁移速率为11.25 km/Myr，并在最后0.4 Myr迁移速率逐渐近乎停滞（图7b）。

该模型的地层样式整体上与模型2类似，主要包括类“下降期体系域”（FSST）和类“低位体系域”（LST）（图7a，b），但其总沉积物通量变化特征具有较大的差异（图8c）。类“下降期体系域”（FSST）发育时间为初始的1.6 Myr，主要表现为下降型迁移轨迹和显著陆上剥蚀现象（图7b）。该时期研究区沉积物通量先快速降至52 km3/Myr，并缓慢上升至72 km3/Myr（图7d）。类“低位体系域”（LST）仅在最后0.4 Myr发育，主要发育进积型地层叠置样式。同时期模拟研究区的总沉积物通量逐渐上升，并最终达到峰值75 km3/Myr（图7d）。

图7 模型4地层剖面、平面模型及相关沉积物通量变化特征（a）剖面D—D’地层叠置样式；（b）Wheeler域剖面；（c）模拟区剥蚀/沉积厚度平面展布；（d）模拟区沉积物通量直方图（5万年为间隔）；分水岭的形成实为模拟的边界效应所致Fig.7 Stratigraphic profile and planar model of model run 4 and variation of sediment flux(a) stratigraphic stacking pattern of profile D￣D’; (b) Wheeler diagram; (c) distribution of denudation and sedimentation; (d) histogram of sediment flux per 50 kyr. The formation of drainage divide is caused by the modeling boundary effect

3 讨论

3.1 滨线迁移的指示意义

虽然掀斜作用对断陷盆地缓坡带湖平面变化的影响在先前相关研究中均有所涉及[8￣9]，但是本次基于简化模型的模拟研究表明其变化特征比概念模式更为复杂。即使在相同构造模型边界条件下，仍然可能具有差异性的滨线迁移特征。模式1和模型3中相对湖平面呈上升趋势，δA为正值。正常情况下，δA与δS差值也为正值（图8a，c），表现为整体湖侵的特征（图9a）。然而实际上δA与δS差值却并不总是正值（图8a，c），表明部分时期沉积滨线处新增沉积物沉积厚度大于新增可容空间的深度。在这种情况下（即δA-δS<0），即使相对湖平面不断上升，但由于三角洲沉积速率相对较快，总可容空间也会不断降低[22]，也会产生正常湖退的特征（图9b）。

模型2、4整体呈相对湖平面下降的趋势，因而造成强制湖退现象（图9c），因此δA为负值。δA与δS差值基本均为负值，且δA与δS差值逐渐趋近于零（图8b，d）。这与两个模型滨线迁移速率降低密切相关，湖岸线开始趋于稳定，δA逐渐趋近于零值。由此可见，在基于构造模型B建立的模型2、4中，其滨线迁移轨迹主要受控于可容空间变化速率。

图9 构造掀斜作用下缓坡带湖平面变化新模式（a，b）枢纽点位于湖平面之上；（c）枢纽点位于湖平面之下Fig.9 New schematic model of lake level fluctuation of the hanging￣wall dipslope driven by tectonic tilting(a, b) pivotal point above lake level; (c) pivotal point below the lake level

3.2 沉积物供给对地层建造的作用

与经典被动陆缘层序相比，掀斜型缓坡层序样式的复杂性在于其滨线迁移轨迹与地层叠置样式的不一致性。模型1、3的滨线迁移轨迹属于向陆迁移型，但是地层地层叠置样式并不以退积型为主导；模型2、4的滨线迁移轨迹属于下降型，但其地层叠置样式亦可发育加积型特征。由此表明，迁移轨迹分析方法对于其层序构型的分析效果并不理想。这种滨线迁移和地层叠置样式不一致性现象的存在与沉积物供给的主导作用密切相关。前人认为掀斜作用会使缓坡带坡度增大，增强了沉积物供给能力，并最终造成了整体湖侵趋势下三角洲进积现象[30]。断块翘倾旋转作用使得苏伊士湾海相断陷盆地中新世碳酸盐沉积向混积岩转化，证实这种构造作用具有增强陆源碎屑沉积物供给能力的作用[31]。

构造掀斜作用并不一定造成沉积物供给能力的增强。本模型所设置的初始地貌为低地丘陵（最大地形高差仅为200 m），在高沉积物供给条件下，其地形高差和陡坡会迅速减小，源区逐渐呈夷平化地貌特征。该背景下掀斜作用所能增强的沉积物通量是极为有限的。模型1分析表明，在整体湖侵的背景下，源区的汇水面积不断缩小，沉积物供给能力并不一定提高。与前人认识不同的是，沉积物供给速率在初始0.4 Myr就可能达到峰值，有利于三角洲向盆地方向进积，因而在初始阶段就存在这种不一致性（图4a，b）。随着沉积物供给速率的降低，这种层序样式与迁移轨迹的不一致性逐渐消失。在整体湖退趋势下，模型2沉积物速率逐渐减少，使得后期三角洲沉积作用减弱，因而主要表现为加积样式，而非进积样式（图5a，b）。

在低沉积物供给条件下，源区风化剥蚀速率较低，模型3、4的源区水系在2 Myr之后仍处于壮年期（图6c、图7c）。该背景下构造掀斜作用则会进一步增强研究区的总沉积物通量（图6d、图7d）。模型3的汇水面积不断缩小，但沉积供给能力却略微上升，表明掀斜作用贡献程度相对更大。其地层叠置样式却并不整体表现为进积特征，也证实沉积物供给速率增幅较为有限。在模型4中，在旋转掀斜和强制湖退的共同作用下，沉积物供给能力快速增强，使得三角洲进积速率增大，因而在研究区与沉积滨线的迁移轨迹保持一致性。

综上所述，在两种不同沉积物供给条件下，沉积物通量大小和变化特征均有很大差异，是造成地层叠置样式和相关“体系域”类型差异的一个关键因素。

3.3 数值模拟结果相关启示及实例应用

本次正演模拟研究表明断陷盆地层序具有复杂性和多样性，并不是被动陆缘层序的简单套用。虽然通常认为构造掀斜作用直接影响可容空间的变化，但实质上是其枢纽点与湖平面相对位置进而决定了可容空间的增减。响应于源区地貌的沉积物通量变化对地层垂向分布特征具有极为重要的控制作用，尤其是在沉积物供给主导情况下形成一系列与湖平面变化相悖的地层样式（如相对湖平面上升背景下发育类“体系域”三分特征，见图4；相对湖平面下降背景下出现地层超覆特征，见图5）。因此针对这种掀斜型缓坡带层序，简单的体系域划分实际上是较为困难的，需要结合地层叠置样式和滨线迁移轨迹共同系统解释。由于本模拟研究设置的初始地形仅为相对单一的低地丘陵，源区总剥蚀产物是相对有限的。因此，本模拟研究仅代表其中一种或几种缓坡带层序可能性的发育模式，对断陷盆地内部掀斜断块形成的局部物源体系（如凹陷内的低凸起型物源体系）更有参考价值。通过多种地貌初始条件，进一步细化构造掀斜模型（甚至考虑反向掀斜、枢纽点迁移以及幕式构造运动等多因素）及气候模型（如季风影响、雨影效应、地形效应等），能够建立更完善、更具多样性的缓坡带层序，以发展和丰富断陷湖盆层序地层学理论。

由于本研究所设定的绝对湖平面边界条件十分理想，现实中很难找到实例进行佐证[1]。地处高纬度的贝加尔湖在晚更新世以来可能符合类似的条件（图10a，b）[32]，但其发育的高频层序（相当于六级层序）却远低于本次正演模拟的研究精度。然而如果将构造和沉积间隔细化到千年尺度，取得的精细模拟结果至少能够反映枢纽点位置与相对湖平面变化的关系。前人研究表明，贝加尔湖自海洋同位素阶段3期（MIS-3期）以来基本保持绝对湖平面长期稳定的状态（图10c）[33￣34]。由于湖平面下降超过2 m就使其成为闭流湖盆，因此该时期气候驱动的湖平面下降持续时间并不长[35]。其缓坡带Selenga三角洲发育区地震剖面所表现的地层叠置样式及局部相对湖平面变化曲线表明该区域相对湖平面变化频繁（图10d）[36]。基于前人研究成果，这种整体后退型的滨线迁移轨迹表明研究区的相对湖平面自MIS-3期以来整体呈逐渐上升的趋势（图10e）。虽然气候型层序观点（气候主导的相对湖平面变化）能够解释这种滨线迁移轨迹发育特征[36]，然而这种高纬度湖泊对于气候变化的敏感度并不高，本质上属于构造主控型湖盆[32]。在绝对湖平面恒定的条件下，贝加尔湖缓坡带Selenga三角洲自MIS-3期以来发育的地震地层单元IIa、Ib及Ia的垂向叠置样式可用掀斜型缓坡层序模式进行成因解释（图10d）。前人地震地层学研究亦表明，Selenga三角洲发育区具有典型的构造掀斜特征，并伴随显著的湖侵作用[32]。基于本次正演模拟的相关研究成果，其构造掀斜枢纽点推测位于湖平面之上。地震地层单元Ib降积样式的发育是海洋同位素阶段2/3期（MIS-2/3）界限处全球气候恶化所致[37]。该时期气候变化造成贝加尔湖的湖平面在短期内可能迅速下降数十米，并形成闭流湖盆[37￣38]。在构造掀斜作用下，该缓坡带的相对湖平面在新仙女木事件之后表现为持续上升的趋势（图10e）。

图10 （a，b）贝加尔湖Selenga三角洲地理位置；（c）中—晚更新世以来贝加尔湖（相对及绝对）湖平面变化特征；（d）基于地震解释成果的地层剖面E—E’；（e）贝加尔湖Selenga三角洲发育区NW区块及SE区块相对湖平面变化特征（修改自文献[33,36]）R-P型：退积—进积型叠置样式；D-A型：降积—加积型叠置样式Fig.10 (a, b) Location of Selenga delta in Lake Baikal; (c) (relative and absolute) lake level fluctuation of Lake Baikal since the Middle to Upper Pleistocene; (d) stratigraphic profile E￣E’ based on seismic interpretation; and (e) relative lake levels of NW and SE blocks of the Selenga delta in Lake Baikal (modified from references [33,36])R￣P = retrogradation-progradation stacking pattern; D￣A = degradation￣aggradation stacking pattern