渤海湾盆地临清坳陷的数值模拟与异常沉降

王子扬邱亮王昆王彦举苏子俊向禹艮

渤海湾盆地临清坳陷的数值模拟与异常沉降

王子扬1,邱亮2*,王昆1,王彦举1,苏子俊1,向禹艮1

1. 山东科技大学地球科学与工程学院,山东 青岛 266590 2. 中国地质大学（北京）地球科学与资源学院, 北京 100083

为探索渤海湾盆地南部的临清坳陷机理，本文利用回剥法和改进的有限元应变速率反演技术，基于3个剖面14口实际井和合成井的数据，获取了拉伸因子、应变速率、计算沉降量等数据，揭示了基于观测沉降量和计算沉降量井差值的异常沉降现象：中新世时期负差值的减速沉降，和上新世至第四纪时期正差值的加速沉降。平衡剖面恢复和区域热流值计算则指示了中新世减速沉降的机制主要为热源的持续流入，抑制了岩石圈的热收缩；中新世以来的加速沉降则与后裂谷期的正断层活动密切相关。

临清坳陷; 地面沉降; 数值模拟

裂后沉降是指发生在断裂构造控制的同裂谷阶段之后的热收缩沉降。前人研究认为，岩石圈伸展主导的裂谷沉降后，热收缩作用主导了裂后沉降阶段[1,2]。然而，最近的研究表明，野外观测获取的实际裂后沉降与麦肯齐热收缩沉降模型的预测结果存在差异，这种差异被称为加速沉降[3]。加速沉降在被动大陆边缘盆地和大型裂谷盆地均有发现，如欧洲的罗克尔盆地和法罗设得兰盆地，南海的马来盆地、珠江口盆地，及东亚的松辽盆地和渤海湾盆地等陆相裂谷盆地中也有发现[4-6]。渤海湾盆地的研究表明，加速沉降现象集中于后期，并没有贯穿整个裂后沉降阶段[7]。那么，裂后阶段是否存在减速沉降的现象，加速沉降和减速沉降及其转变的的地球动力学机制是什么。

临清坳陷位于太平洋板块在东亚下方地幔过渡区的滞留板片西缘与重力异常交界处，是研究深部地质过程与上层地壳盆地演化之间相互作用的理想区域。本文选取位于渤海湾盆地西南部的临清坳陷为研究对象，重点研究了后裂谷阶段加速沉降和减速沉降的产生及其机制，利用包含14口真实或合成井的3个地震截面进行建模，从中获得拉伸因子、应变速率、预测沉降和古热流等数据，揭示临清坳陷的后裂谷期构造演化。

1 研究概况

临清坳陷是一个中、新生代复合型断陷沉降区，在构造区划上属于渤海湾盆地二级负向构造单元，位于渤海湾盆地的西南收敛端，呈北北东向展布，北接黄骅坳陷、沧县隆起，南至内黄隆起，西邻太行山隆起，东至鲁西南隆起[7,8]。临清坳陷构造演化经历了古近纪始新世～渐新世快速沉降与堆积的断陷阶段以及新近纪中新世一第四纪连续沉积的区域性坳陷阶段[9]。区内主要发育中-新生代盖层，新生代地层是中国重要的油气富集层位，主要由古近纪沙河街组、东营组，新近纪馆陶组和明化镇组及第四纪地层组成[10]。

2 材料与方法

2.1 数据来源及参数

本次研究中用于沉降分析的数据来自临清坳陷的3个地震剖面（从南到北依次为1、2和3）中的3口实井（标记为s）和11口合成井（标记为h）。基于地震剖面和实际井的分析，对模拟井的深度和岩性进行了推算。通过对收集到的地层岩性、厚度和孔隙度等基础数据进行分解和回剥，重建观测到的构造沉降厚度。地壳和岩石圈厚度及其他参数均来自前人研究（表1）；为了简化计算，我们将古水深定为零。

表 1 参数及取值[8-10]

2.2 实际构造沉降量

用于计算构造沉降的软件是斯伦贝谢公司开发的基于回剥法的PetroMod 2012.1。回剥法是一种基于Airy均衡模型的反演方法，包括以下几种步骤：分解、古水深校正、沉积界面温度校正、剥蚀校正、海面起伏和构造沉降提取[11]。在回剥过程中还应考虑跨越780万年，从距今23.8百万年至距今16.0百万年（23.8～16.0 Ma）的剥蚀作用（即东营运动）。为了便于比较不同剖面的应变速率和断层活跃率，假设东营运动期间沉积和剥蚀速率为零。

表 2 回剥法中使用的岩石力学参数[11]

2.3 模拟构造沉降量

确定预测沉降需要两个步骤。第一步，通过输入背剥曲线进行改良的一维应变速率反演，获得应变总量和应变速率。第二步，采用改良的有限可拓模型，以为输入参数，计算预测沉降。

2.3.1 一维应变速率反演前人研究提出了一种基于构造沉降史计算应变总量和应变速率的新算法（式1）[12,13]。

式中，()是应变速率对时间的函数，由于方程1难以得到解析解，可采用鲍威尔算法进行数值迭代计算。然而，鲍威尔算法作为一种共轭收敛加速方法，需要大量数据，而沉降曲线无法提供这些数据。此外，沉降曲线的线性依赖性可能导致搜索退化，导致不收敛问题[14]。本文对传统方法进行了改进，采用牛顿二阶插值迭代法求得数值解；将积分转化为和，就可以求得的值（式2）。

式中，是中点回归参数，G是中点应变速率，,和是鲍威尔和谐矩阵的初始取值，T是伸展起始时间，其它字符见表1。

2.3.2 改进的有限扩展模型有限扩展模型是Jarvis等提出的经典均匀伸展模型的改进版本。有限扩展模型不仅考虑了均匀伸展模型中包含的热平流项，还在热流方程中加入了热扩散函数，以反映岩石圈结构中温度随时间变化的情况[15]。但是，来自地壳的放射性生热在盆地发育中具有重要作用，在模拟盆地演化过程时，必须考虑放射性生热。通过确定大陆地温并使其等同于基础热流，可以用以深度为变量的负指数模型来模拟地壳放射性生热。

式中，等式右侧第一项代表水平热流值，第二项代表热扩散函数，第三项代表了放射性热流值。

由于应变总量是沉降量的函数，因此对上式进行变形，得到的表达式：

式中，y是利用回剥法计算得出的实际构造沉降量，其它参数含义见表1。

3 构造沉降

3.1 凹陷南部剖面A

剖面A横切整个坳陷盆地，西部受分别西倾和东倾的断层共同控制，形成不对称地堑，东部受控盆断层控制，形成地堑构造，中央隆起充当两侧断层的共同下盘。剖面西缘断层活动不占优势，古生代地层与断层带重叠，形成西部斜坡带（图1a）。

图 1 横切坳陷南部的构造剖面A

（a）根据地震剖面解释的构造剖面，井sa1～sa6的位置；（b～g）sa1～sa6井的伸展总量、应变速率和观测沉降（蓝线）与理论沉降（橙色线）曲线。

剖面A的裂谷期地层为古近系沙河街组和东营组组成，平均应变总量约为1.75。裂谷期的应变速率变化指示了两幕次的裂谷活动时期，对应形成了沙河街组三段～二段（Es3-Es2）和东营组（Ed）地层。Es3-2活动时期以Sa6井为中心，最大应变速率达～2×10-15s-1，值最大为1.80。Ed时期以Ma3井为中心，最大应变率是～1.3×10-15s-1，值最大为1.70。根据不同的应变率可以观察到不同井构造沉降的变化。Es3-2应变率的递增顺序是Ma6、Ma5、Ma3、Ma1、Ma2、Sa4。Ed应变率的递增顺序为Ma3、Ma2、Ma1、Ma6、Sa4、Ma5（图1b-g），两幕活动的应变速率表现此消彼长的数值关系。

东营组和馆陶组间的平行不整合分隔了裂谷期与后裂谷期。中新世的后裂谷期应变率为零，值保持不变，与McKenzie的理论模型一致。然而，上新世和第四纪的应变率和值明显增加，表明该时期出现快速沉降（图1b-g）。根据沉降的观测值（蓝线）与预测值（粉线）的区别，可以确定中新世以来的理论预测值与实际观测值之间存在差异（图1b-g）。沉降差在中新世期为-96～328 m，在上新世期间为+237～407 m。中新世观测沉降值的倾角平缓，而上新世-第四纪则较为陡峭。

3.2 凹陷中部剖面B

凹陷中部是同裂谷期断层活动最强烈的区域。中部岩陷的结构是典型的半地堑，以西倾的铲状兰聊断层为主。同裂谷期的两幕活动模式（Es3-2和Ed）由应变速率图清晰可见（图2b-f），平均应变率和值分别为2.2×10-15s-1和1.78，Es3-2时期断层活动较强烈；Ed期以弱断层为特征，应变率只有0.5×10-15s-1。Es3-2和Ed的应变速率呈现正相关，变化模式与剖面A相反，整个凹陷中部的沉降受兰聊断层的完全控制。

图 2 横切坳陷中部的构造剖面B

后裂谷期的应变速率在中新世均为零，而在上新世-第四纪期间明显增加，平均值为～0.4×10-15s-1；值在中新世保持不变，在上新世-第四纪明显增加，平均值为1.82。中新世观测沉降直线的倾角平缓，而上新世-第四纪更为陡峭，蓝、粉线之间的区域代表观测沉降量和计算沉降量之间的差值，中新世两者的沉降差为-266～343 m，上新世则为+123～247 m。

3.3 凹陷北部剖面C

北部凹陷是由兰聊断层控制的复合半地堑，结合上盘发育的对立的次级断层。这些次级结构以旋转的平面正断层为标志，形成多米诺骨牌断层系统。单峰应变率模式表明，该地区出现了独特的Es3-2断层期，与南部和中部的裂谷发育模式明显不同。Ed时期多为弱断层，特点是应变率低，约为0.13×10-15s-1。

图 3 横切坳陷北部的构造剖面C

后裂谷阶段的应变速率模式和应变总量模式与南部和中部凹陷类似。中新世时期应变速率为零应变总量保持不变，在上新世-第四纪期间增至平均值约为0.4×10-15s-1，应变总量持续增长至第四纪。中新世两线之差为-185～218 m，在上新世-第四纪为+266～291 m（图3b-d）。通过北部凹陷地层的厚度和应变率可知，在新近纪期间，沉积中心迁移到了西部洼陷。

4 讨论

4.1 加速沉降与减速沉降

根据均匀伸展模型，岩石圈在同裂谷末期断层活动逐渐停止，后裂谷初期热收缩沉降启动。应变速率ε和应变总量是评价同裂谷期断层伸展强度的主要定量指标。当应变速率为零，应变总量保持不变，即表示正断层活动的终止。由于总热流随时间呈现负指数型变化，热沉降曲线因此表现为负指数模型[16-18]。

尽管研究区的中新世后裂谷期应变速率为零、应变总量值保持稳定，符合均匀伸展模型的预测，但观测沉降量与计算沉降量存在较大差异，观测沉降量明显小于计算沉降量，该负差值范围为-343～96 m，这指示了一种非常规的减速沉降模式。由于在23.8～16 Ma阶段整个渤海湾盆地进入了一个剥蚀夷平期，该减速沉降起始于16 Ma。新近纪-第四纪阶段，应变速率>0、应变总量不断增长的情况下，观测沉降量与计算沉降量之间的正差值范围是123～407 m，这表示上新世-第四纪期间出现了加速沉降。渤海湾盆地的济阳凹陷和渤中凹陷均已观测到加速沉降现象，减速沉降则是首次被发现。

在C剖面的所有井与剖面A的Ma1到Ma4井中，当前（0 Ma时）的观测沉降量与计算沉降量之间均为正差值，其余剖面的井均为负差值。该现象具有两种不同的分布类型，正差异的特点是蓝线（观测沉降值）与粉线（预测沉降值）出现交叉点；负差异则以不相交为特征。蓝、粉线相交的地质意义是，后期的加速沉降补充了早期减速沉降造成的厚度不足，剩余量为正差值表明后裂谷期总沉降高于均匀伸展模型预测的热沉降量。两线不相交则表明，后期加速沉降量不足以弥补早期减速沉降的缺失，产生了50～150 m的负差异。但按照当前蓝线、粉线的延长趋势，未来加速沉降可能能够补偿减速沉降造成的沉降量缺失。

4.2 裂后期正断层活动与加速沉降

根据均匀伸展模型，正断层活动在后裂谷期停止，但渤海湾盆地的渤中凹陷和冀中凹陷均发现了裂后期的正断层活动[12,13]。因此可以证明，正断层的重启造成了渤海湾盆地的裂后期加速沉降。

在上新世-第四纪期间，临清坳陷也出现了裂后期的正断层活动。如图4所示，基底断层快速复活的特点是上盘的沉积层比下盘更厚。断层的活动速率表明断层活动相对强烈，这一时期的活动速率甚至大于Es4、Es2等裂谷期断层活动相对平静阶段。因此，后裂谷期的正断层活动与加速沉降密切相关。后裂谷期正断层活动的机制较为复杂，人们提出了几种构造模型来解释这种现象，如郯庐断层的右旋运动，新裂谷活动期和华北克拉通破坏的重新启动等，对该问题的成因还有争议。

平衡剖面和断层活动率表明，中新世期缺少正断层活动。同时，正断层活动的结果只能导致加速沉降发生，与中新世的减速沉降因果关系较小。因此，中新世的减速沉降应归因于其他构造因素。

图 4 后裂谷期平衡地质剖面

4.3 古热流与减速沉降

岩石圈减薄后的热收缩是造成裂后沉降的主要机制，由于热流减少，热收缩沉降的速度随时间呈指数下降。因此，热流是裂后沉降的直接控制因素，加速和减速沉降的出现可能会反映热流的转变。为了探究沉降变化与古热流之间的关系，本文设计了基于热流模型预测值和观测值的综合分析。以平均应变总量=1.75的剖面A为例，使用改进的模拟热流方程对剖面A的预测平均热流和沉降进行估算。

由于后裂谷期沉降是同裂谷期沉降的函数，可以据此调整模型获得热流的观测值：给从应变率反演中获得的值增加一个Δ值（如0.05），计算出的沉降量也相应的增加。当结果与观测的裂后沉降相吻合时，可认为此时的热流与沉降的观测值相对应。该值以式5的左边表示。

其中式5等号左侧是与观测的沉降相对应的热流。六口井的数据处理完成后，可采用线性指数拟合的方法得到与沉降观测值相对应的热流曲线（图5a）。临清坳陷的古热流模型得出的数值为61～94 mW/m2，与本研究估计的热流值一致。

a 比例地质年代尺度为20 Ma，比例热流值以裂谷发育前的初始热流为1

b 示模型计算的理论热流值，红线表示根据实际沉降量反演的实际热流值

根据剖面A的模拟结果，在中新世期间，热流呈线性下降，随后在上新世-第四纪期间表现为曲线下降的趋势。在整个裂后阶段，热流的观测值比预测值要大（图5b）。在7.7 Ma前后出现了两条曲线的最大差异（约13.67）。在加速沉降期（18.5～5.3 Ma），观测和预测的热流差在3.01～13.67 m之间，而该值在减速沉降期约为10.44～12.67 mm2。

在中新世，热收缩损失的热量将由其他热源继续提供，尽管不足以阻止热流的减少，但可以减缓热流减少的速度，抑制岩石圈的热收缩。因此，热流减速过慢导致热收缩不充分，沉降量并不明显。在上新世早期，热流冷却模型从线性过渡到曲线，其他热源消失，热收缩以McKenzie的负指数模型重新启动。热流值约为78.5 m2，比预测值高16%（图5b），其产生的热收缩能量更强，导致明化镇组出现了超充分沉降。因此，加速、减速沉降的出现与热流变化有明显的相关性。

5 结论

本文针对渤海湾盆地临清坳陷的后裂谷阶段，利用包含14口真实或合成井的3个地震截面，采用回剥法和改进的有限元应变速率反演技术进行数值模拟，获取了拉伸因子、应变速率、计算沉降量和古热流等的数据。后裂谷期观测沉降量与计算沉降量的对比揭示了中新世的减速沉降，沉降差值为−96～343 m，和上新世-第四纪的加速沉降，沉降差值为+123～407 m。通过平衡剖面恢复和区域后裂谷期热流值计算，认为中新世后裂谷期减速沉降主要归因于热源的持续流入，抑制了岩石圈的热收缩；中新世加速沉降则与后裂谷期的正断层活动密切相关。

[1] McKenzie D. Some remarks on the development of sedimentary basins [J]. Earth Planetary Science Letters, 1978,40(1):25-32

[2] McKenzie D, Jackson J, Priestley K. Thermal structure of oceanic and continental lithosphere [J]. Earth Planetary Science Letters, 2005,233(3-4):337-349

[3] McKenzie D, Watts A, Parsons B,. Planform of mantle convection beneath the Pacific Ocean [J]. Nature, 1980,288:442-446

[4] Madon M, Watts A. Gravity anomalies, subsidence history and the tectonic evolution of the Malay and Penyu Basins ( offshore Peninsular Malaysia ) [J]. Basin Research, 1998,10:375-392

[5] 漆家福.渤海湾新生代盆地的两种构造系统及其成因解释[J].中国地质,2004,31(1):15-22

[6] Qi J, Yang Q. Cenozoic structural deformation and dynamic processes of the Bohai Bay Basin province, China [J]. Marine and Petroleum Geology, 2010,27(4):757-771

[7] Qiu NS, Xu W, Zuo YH,. Meso-Cenozoic thermal regime in the Bohai Bay Basin, eastern North China Craton [J]. International Geology Review, 2015,57(3):271-289

[8] Qiu NS, Zuo YH, Chang J,. Geothermal evidence of Meso-Cenozoic lithosphere thinning in the Jiyang sub-basin, Bohai Bay Basin, eastern North China Craton [J]. Gondwana Research, 2014,26:1079-1092

[9] 任凤楼,柳忠泉,邱连贵,等.渤海湾盆地新生代各坳陷沉降的时空差异性[J].地质科学,2008,43(3):546-557,575

[10] Ren J, Tamaki K, Li S,. Late Mesozoic and Cenozoic rifting and its dynamic setting in Eastern China and adjacent areas [J]. Tectonophysics, 2002,344(3-4):175-205

[11] 史卜庆,吴智平,王纪祥,等.渤海湾盆地东营运动的特征及成因分析[J].石油实验地质,1999,21(3):196-200

[12] White N. Recovery of strain rate variation from inversion of subsidence data [J]. Nature, 1993,366(6454):499-452

[13] White N. An inverse method for determining lithospheric strain rate variation on geological timescales [J]. Earth Planetary Science Letters, 1994,122(3-4):351-371

[14] Silva J. Influence diagnostics and estimation algorithms for powell's SCLS [J]. Journal Bus. Economic Statistics, 2001,19(1):55-62

[15] 肖国强,杨吉龙,赵长荣,等.天津滨海地区G2孔磁性地层年代及其构造指示[J].地质通报,2014(10):1642-1650

[16] Song B, Chen L, Zhang J,A matlab program for 1d strain rate inversion [J]. Computer Geoscience, 2010,36(1):16-23

[17] 万天丰.论构造事件的节律性[J].地学前缘,1997,4(3-4):257-263

[18] Wheeler P, White N. Quest for dynamic topography: observations from Southeast Asia [J]. Geology, 2000,28(11):963-966

Numerical Simulation and Abnormal Sedimentation of Linqing Depression on Bohai Bay Basin

WANG Zi-yang1, QIU Liang2*, WANG Kun1, WANG Yan-ju1, SU Zi-jun1, XIANG Yu-gen1

1.266590,2.100083,

In order to explore the mechanism of the Linqing depression in the southern Bohai Bay Basin, this paper uses the back-stripping method and the improved finite element strain rate inversion technique. Based on the data of 14 actual and synthetic Wells in 3 profiles, the data of tensile factor, strain rate and calculated settlement are obtained, and the abnormal settlement phenomenon based on the well difference between observed settlement and calculated settlement is revealed: The deceleration of negative difference in Miocene period, and the acceleration of positive difference in Pliocene to Quaternary period. The restoration of the equilibrium profile and the calculation of regional heat flow indicate that the mechanism of the slow subsidence in Miocene is mainly the continuous inflow of heat source, which inhibits the thermal contraction of the lithosphere. The accelerated subsidence since Miocene is closely related to the normal fault activity in the post-rift period.

Lingqing depression; ground subsidence; numerical simulation

P642.26

A

1000-2324(2022)06-0956-07

2022-05-02

2022-07-21

国家自然科学基金委青年基金(41702207)

王子扬(2002-),男,本科生,主要从事工程地质及地面沉降机理探究. E-mail:2652382461@qq.com

Author for correpondence. E-mail:qiuliang2011@gmail.com

10.3969/j.issn.1000-2324.2022.06.023