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

第6期侯恩光等:沂河流域山东段径流演变特征及影响因素分析·955·
可分别得到因降水量和人类活动影响而各自产生的径流影响量。

突变年之前沂河流域的实测平均径流量为21.787×108m3，理论平均径流量为22.5×108m3；突变年之后沂河流域的实测平均径流量为
25.134×108m3，理论平均径流量为23.854×108m3。

具体如表5所示。

表5降水量和人类活动对沂河流域径流量影响的统计值
Table5Contribution rate of precipitation and human activities to runoff
年限Duration
平均降水量/mm
Average precipitation
实测平均径流量/×108m3
Live average runoff
理论平均径流量/×108m3
Theoretical average runoff
降水影响量/×108m3
Effects of precipitation
人类活动影响量/×108m3
Effects of humen
1965-1988739.93421.78722.5 2.067 1.28
1989-2016789.88225.13423.854
由表5可知，降水量对沂河流域径流量增大的影响量为2.067×108m3，贡献率为61.76%，而人类活动的影响量为1.28×108m3，贡献率为38.24%，故由此可以量化得出人类活动对沂河流域径流量增大造成了较大的影响。

4结论
（1）沂河流域1956-2016年径流量总体呈现下降趋势；1965年为径流量的最重要突变点，其次为1989年；23年、11年和6年分别对应径流量变化的第1、第2和第3主周期；
（2）沂河流域1956-2016年降水量总体呈现下降趋势，1964年为降水量的最重要突变点，其次为1989年；在1956-1969年、1970-1989年和1990-2016年分别以10-12年、14-18年和20-30年的时间尺度呈周期性变化，23年、11年和7年分别对应降水量变化的第1、第2和第3主周期；
（3）在各个5年时段沂河流域径流量的极值比和变差系数均比降水量的大，径流量与降水量的变化趋势并非完全一致，径流量的变化幅度更大且变化效果更为显著；
（4）沂河流域径流量和降水量在1965-2016年的最重要突变点为1989年；降水量和人类活动对该流域径流量增大的贡献率分别为61.76%和38.24%，由此可量化表明人类活动对沂河流域径流量增大造成了较大的影响。

参考文献
[1]Zhang HX,Zhang L,Zhu RR,et al.Responses of streamflow to climate and land surface change in the headwaters of
the Yellow River Basin[J].Water Resources Research,2009,45(7):641-648
[2]孙琪,徐长春,任正良,等.塔里木河流域产水量时空分布及驱动因素分析[J].灌溉排水学报,2021,40(8):114-122
[3]Wang SJ,Yan M,Yan YX,et al.Contributions of climate change and human activities to the changes in runoff
increment in different sections of the Yellow River[J].Quaternary International,2012,282:66-77
[4]孙洋洋,穆兴民,高鹏,等.黄河上中游区径流量格局长期变化过程[J].水土保持研究,2017,24(3):59-64
[5]张海荣,周建中,曾小凡,等.金沙江流域降水和径流时空演变的非一致性分析[J].水文,2015,35(6):90-96
[6]潘扎荣,阮晓红,朱愿福,等.近50年来淮河干流径流演变规律分析[J].水土保持学报,2013,27(1):51-55,59
[7]Tao H,Gemmer M,Bai Y,et al.Trends of stream-flow in the tarim river basin during the past50years:human impact
or climate change[J].Journal of Hydrology,2011,400(1/2):1-9
[8]宋小燕,穆兴民,高鹏,等.松花江哈尔滨站近100年来径流量变化趋势[J].自然资源学报,2009,24(10):1803-1809
[9]郭爱军,畅建霞,王义民,等.近50年泾河流域降雨-径流关系变化及驱动因素定量分析[J].农业工程学报,2015,
31(14):165-171
[10]刘昌明,王恺文,王冠,等.1956-2016年黄河流域径流及其影响因素的变化分析[J].人民黄河,2022,44(9):1-5,16
[11]姬世保,杜军凯,仇亚琴,等.基于二元水循环模拟的密云水库上游流域河川径流演变归因分析[J].水利水电技术(中
英文),2022,53(4):116-127
[12]阳扬,翟禄新,贾艳红,等.不同土地利用方式下漓江流域降水-径流特征变化分析研究[J].中国农村水利水电,2022(8):97-105
[13]吴子怡,谢平,桑燕芳,等.水文序列跳跃变异点的滑动相关系数识别方法[J].水利学报,2017,48(12):1473-1481,1489
[14]刘茜,王延贵.江河水沙变化突变性与周期性分析方法及比较[J].水利水电科技进展,2015,35(2):17-23
[15]胡彩虹,王艺璇,管新建,等.基于双累积曲线法的径流变化成因分析[J].水资源研究,2012,1(4):204-210
[16]穆兴民,张秀勤,高鹏,等.双累积曲线方法理论及在水文气象领域应用中应注意的问题[J].水文,2010,30(4):47-51
[17]张琳,马敬旭,张倩,等.近60多年洞庭湖水沙演变特征及其与人类活动的关系[J].长江科学院院报,2021,38(9):14-20
[18]杨思雨,姜仁贵,解建仓,等.渭河上游径流变化特征及归因研究[J].水资源与水工程学报,2019,30(2):37-42
Copyright©博看网. All Rights Reserved.
山东农业大学学报(自然科学版),2022,53(6):956-962VOL.53NO.62022 Journal of Shandong Agricultural University(Natural Science Edition)doi:10.3969/j.issn.1000-2324.2022.06.023
渤海湾盆地临清坳陷的数值模拟与异常沉降
王子扬1,邱亮2*,王昆1,王彦举1,苏子俊1,向禹艮1
1.山东科技大学地球科学与工程学院,山东青岛266590
2.中国地质大学（北京）地球科学与资源学院,北京100083
摘要:为探索渤海湾盆地南部的临清坳陷机理，本文利用回剥法和改进的有限元应变速率反演技术，基于3个剖面14口实际井和合成井的数据，获取了拉伸因子、应变速率、计算沉降量等数据，揭示了基于观测沉降量和计算沉降量井差值的异常沉降现象：中新世时期负差值的减速沉降，和上新世至第四纪时期正差值的加速沉降。

平衡剖面恢复和区域热流值计算则指示了中新世减速沉降的机制主要为热源的持续流入，抑制了岩石圈的热收缩；中新世以来的加速沉降则与后裂谷期的正断层活动密切相关。

关键词:临清坳陷;地面沉降;数值模拟
中图法分类号:P642.26文献标识码:A文章编号:1000-2324(2022)06-0956-07 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.College of Earth Sciences and Engineering/Shandong University of Science and Technology,Qingdao266590,China
2.School of Earth Sciences and Resources/China University of Geosciences(Beijing),Beijing100083,China Abstract: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 of14actual and synthetic Wells in3profiles,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. Keywords:Lingqing depression;ground subsidence;numerical simulation
裂后沉降是指发生在断裂构造控制的同裂谷阶段之后的热收缩沉降。

前人研究认为，岩石圈伸展主导的裂谷沉降后，热收缩作用主导了裂后沉降阶段[1,2]。

然而，最近的研究表明，野外观测获取的实际裂后沉降与麦肯齐热收缩沉降模型的预测结果存在差异，这种差异被称为加速沉降[3]。

加速沉降在被动大陆边缘盆地和大型裂谷盆地均有发现，如欧洲的罗克尔盆地和法罗设得兰盆地，南海的马来盆地、珠江口盆地，及东亚的松辽盆地和渤海湾盆地等陆相裂谷盆地中也有发现[4-6]。

渤海湾盆地的研究表明，加速沉降现象集中于后期，并没有贯穿整个裂后沉降阶段[7]。

那么，裂后阶段是否存在减速沉降的现象，加速沉降和减速沉降及其转变的的地球动力学机制是什么。

临清坳陷位于太平洋板块在东亚下方地幔过渡区的滞留板片西缘与重力异常交界处，是研究深部地质过程与上层地壳盆地演化之间相互作用的理想区域。

本文选取位于渤海湾盆地西南部的临清坳陷为研究对象，重点研究了后裂谷阶段加速沉降和减速沉降的产生及其机制，利用包含14口真实或合成井的3个地震截面进行建模，从中获得拉伸因子、应变速率、预测沉降和古热流等数据，揭示临清坳陷的后裂谷期构造演化。

收稿日期:2022-05-02修回日期:2022-07-21
基金项目:国家自然科学基金委青年基金(41702207)
第1作者简介:王子扬(2002-),男,本科生,主要从事工程地质及地面沉降机理探究.E-mail:*****************
*通讯作者:Author for correpondence.E-mail:**********************
第6期王子扬等:渤海湾盆地临清坳陷的数值模拟与异常沉降·957·
1研究概况
临清坳陷是一个中、新生代复合型断陷沉降区，在构造区划上属于渤海湾盆地二级负向构造单元，位于渤海湾盆地的西南收敛端，呈北北东向展布，北接黄骅坳陷、沧县隆起，南至内黄隆起，西邻太行山隆起，东至鲁西南隆起[7,8]。

临清坳陷构造演化经历了古近纪始新世~渐新世快速沉降与堆积的断陷阶段以及新近纪中新世一第四纪连续沉积的区域性坳陷阶段[9]。

区内主要发育中-新生代盖层，新生代地层是中国重要的油气富集层位，主要由古近纪沙河街组、东营组，新近纪馆陶组和明化镇组及第四纪地层组成[10]。

2材料与方法
2.1数据来源及参数
本次研究中用于沉降分析的数据来自临清坳陷的3个地震剖面（从南到北依次为1、2和3）中的3口实井（标记为s）和11口合成井（标记为h）。

基于地震剖面和实际井的分析，对模拟井的深度和岩性进行了推算。

通过对收集到的地层岩性、厚度和孔隙度等基础数据进行分解和回剥，重建观测到的构造沉降厚度。

地壳和岩石圈厚度及其他参数均来自前人研究（表1）；为了简化计算，我们将古水深定为零。

表1参数及取值[8-10]
Table1Parameters and values in the test
符号Symbol
含义
Implication
值
Value
y l岩石圈厚度80km
y c地壳厚度35km
ρc*地壳密度 2.8g·cm-3
T0地表温度0°C
T a软流圈温度1333°C
ρa软流圈密度 3.2g·cm-3
Ρm*地幔密度 3.33g·cm-3
αv热扩散系数 3.28×10-5°C-1
κ热传导系数10-6m2·s-1
2.2实际构造沉降量
用于计算构造沉降的软件是斯伦贝谢公司开发的基于回剥法的PetroMod2012.1。

回剥法是一种基于Airy均衡模型的反演方法，包括以下几种步骤：分解、古水深校正、沉积界面温度校正、剥蚀校正、海面起伏和构造沉降提取[11]。

在回剥过程中还应考虑跨越780万年，从距今23.8百万年至距今16.0百万年（23.8~16.0Ma）的剥蚀作用（即东营运动）。

为了便于比较不同剖面的应变速率和断层活跃率，假设东营运动期间沉积和剥蚀速率为零。

表2回剥法中使用的岩石力学参数[11]
Table2Rock mechanical parameters used in Back-stripping
岩性Lithology
压实系数
Compacting factor
初始孔隙度
Initial porosity
骨架密度
Framework density
砂岩0.270.49 2.65
页岩0.510.63 2.72
砂砾岩0.330.32 2.70
玄武岩0.250.10 3.00
2.3模拟构造沉降量
确定预测沉降需要两个步骤。

第一步，通过输入背剥曲线进行改良的一维应变速率反演，获得应变总量β和应变速率ε。

第二步，采用改良的有限可拓模型，以β为输入参数，计算预测沉降。

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

·958·山东农业大学学报(自然科学版)第53卷
⎰∆=t
dt t G 0))((exp β(1)
式中，G (t )是应变速率对时间的函数，由于方程1难以得到解析解，可采用鲍威尔算法进行数值迭代计算。

然而，鲍威尔算法作为一种共轭收敛加速方法，需要大量数据，而沉降曲线无法提供这些数据。

此外，沉降曲线的线性依赖性可能导致搜索退化，导致不收敛问题[14]。

本文对传统方法进行了改进，采用牛顿二阶插值迭代法求得数值解；将积分转化为和，就可以求得β的值（式2）。

j m Z l nn in j A nn jj t T hh jj y G rr lmd A ⨯--+⨯+⨯-⨯⨯+=-+1
12()1()1((2)式中，lmd 是中点回归参数，G in 是中点应变速率，rr ,nn 和jj 是鲍威尔和谐矩阵的初始取值，T z 是伸展起始时间，其它字符见表1。

2.3.2改进的有限扩展模型有限扩展模型是Jarvis 等提出的经典均匀伸展模型的改进版本。

有限扩展模型不仅考虑了均匀伸展模型中包含的热平流项，还在热流方程中加入了热扩散函数，以反映岩石圈结构中温度随时间变化的情况[15]。

但是，来自地壳的放射性生热在盆地发育中具有重要作用，在模拟盆地演化过程时，必须考虑放射性生热。

通过确定大陆地温并使其等同于基础热流，可以用以深度为变量的负指数模型来模拟地壳放射性生热。

⎪⎪⎭⎫ ⎝⎛--+∂∂+∂∂-=∂∂r p a a y c A z T k z T t T )(exp 022ρυ(3)
式中，等式右侧第一项代表水平热流值，第二项代表热扩散函数，第三项代表了放射性热流值。

由于应变总量β是沉降量的函数，因此对上式进行变形，得到β的表达式：
⎥⎦
⎤⎢⎣⎡-⎪⎭⎫ ⎝⎛----=****221,),])1([-11m m m l c c m s m m s T a T a y y y T a y l ρρρρρβυυυ（(4)式中，y s 是利用回剥法计算得出的实际构造沉降量，其它参数含义见表1。

3
构造沉降3.1凹陷南部剖面A
剖面A 横切整个坳陷盆地，西部受分别西倾和东倾的断层共同控制，形成不对称地堑，东部受控盆断层控制，形成地堑构造，中央隆起充当两侧断层的共同下盘。

剖面西缘断层活动不占优势，古生代地层与断层带重叠，形成西部斜坡带（图1a
）。

图1横切坳陷南部的构造剖面A
Fig.1The structural section A on south transverse depression
（a ）根据地震剖面解释的构造剖面，井sa1~sa6的位置；（b~g ）sa1~sa6井的伸展总量、应变速率和观测沉降（蓝线）与理论沉降（橙色线）曲线。

第6期王子扬等:渤海湾盆地临清坳陷的数值模拟与异常沉降·959·
剖面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~328m，在上新世期间为+237~407m。

中新世观测沉降值的倾角平缓，而上新世-第四纪则较为陡峭。

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
Fig.2The structural section B on middle transverse depression
后裂谷期的应变速率在中新世均为零，而在上新世-第四纪期间明显增加，平均值为~0.4×10-15 s-1；β值在中新世保持不变，在上新世-第四纪明显增加，平均值为1.82。

中新世观测沉降直线的倾角平缓，而上新世-第四纪更为陡峭，蓝、粉线之间的区域代表观测沉降量和计算沉降量之间的差值，中新世两者的沉降差为-266~343m，上新世则为+123~247m。

3.3凹陷北部剖面C
北部凹陷是由兰聊断层控制的复合半地堑，结合上盘发育的对立的次级断层。

这些次级结构以旋转的平面正断层为标志，形成多米诺骨牌断层系统。

单峰应变率模式表明，该地区出现了独特的
·960·山东农业大学学报(自然科学版)第53卷Es3-2断层期，与南部和中部的裂谷发育模式明显不同。

Ed时期多为弱断层，特点是应变率低，约为0.13×10-15s-1。

图3横切坳陷北部的构造剖面C
Fig.3The structural section C on north transverse depression
后裂谷阶段的应变速率模式和应变总量模式与南部和中部凹陷类似。

中新世时期应变速率为零应变总量保持不变，在上新世-第四纪期间增至平均值约为0.4×10-15s-1，应变总量持续增长至第四纪。

中新世两线之差为-185~218m，在上新世-第四纪为+266~291m（图3b-d）。

通过北部凹陷地层的厚度和应变率可知，在新近纪期间，沉积中心迁移到了西部洼陷。

4讨论
4.1加速沉降与减速沉降
根据均匀伸展模型，岩石圈在同裂谷末期断层活动逐渐停止，后裂谷初期热收缩沉降启动。

应变速率ε和应变总量β是评价同裂谷期断层伸展强度的主要定量指标。

当应变速率为零，应变总量保持不变，即表示正断层活动的终止。

由于总热流随时间呈现负指数型变化，热沉降曲线因此表现为负指数模型[16-18]。

尽管研究区的中新世后裂谷期应变速率为零、应变总量值保持稳定，符合均匀伸展模型的预测，但观测沉降量与计算沉降量存在较大差异，观测沉降量明显小于计算沉降量，该负差值范围为-343~96m，这指示了一种非常规的减速沉降模式。

由于在23.8~16Ma阶段整个渤海湾盆地进入了一个剥蚀夷平期，该减速沉降起始于16Ma。

新近纪-第四纪阶段，应变速率>0、应变总量不断增长的情况下，观测沉降量与计算沉降量之间的正差值范围是123~407m，这表示上新世-第四纪期间出现了加速沉降。

渤海湾盆地的济阳凹陷和渤中凹陷均已观测到加速沉降现象，减速沉降则是首次被发现。

在C剖面的所有井与剖面A的Ma1到Ma4井中，当前（0Ma时）的观测沉降量与计算沉降量之间均为正差值，其余剖面的井均为负差值。

该现象具有两种不同的分布类型，正差异的特点是蓝线（观测沉降值）与粉线（预测沉降值）出现交叉点；负差异则以不相交为特征。

蓝、粉线相交的地质意义是，后期的加速沉降补充了早期减速沉降造成的厚度不足，剩余量为正差值表明后裂谷期总沉降高于均匀伸展模型预测的热沉降量。

两线不相交则表明，后期加速沉降量不足以弥补早期减速沉降的缺失，产生了50~150m的负差异。

但按照当前蓝线、粉线的延长趋势，未来加速沉降可能能够补偿减速沉降造成的沉降量缺失。

4.2裂后期正断层活动与加速沉降
根据均匀伸展模型，正断层活动在后裂谷期停止，但渤海湾盆地的渤中凹陷和冀中凹陷均发现
第6期王子扬等:渤海湾盆地临清坳陷的数值模拟与异常沉降·961·了裂后期的正断层活动[12,13]。

因此可以证明，正断层的重启造成了渤海湾盆地的裂后期加速沉降。

在上新世-第四纪期间，临清坳陷也出现了裂后期的正断层活动。

如图4所示，基底断层快速复活的特点是上盘的沉积层比下盘更厚。

断层的活动速率表明断层活动相对强烈，这一时期的活动速率甚至大于Es 4、Es 2等裂谷期断层活动相对平静阶段。

因此，后裂谷期的正断层活动与加速沉降密切相关。

后裂谷期正断层活动的机制较为复杂，人们提出了几种构造模型来解释这种现象，如郯庐断层的右旋运动，新裂谷活动期和华北克拉通破坏的重新启动等，对该问题的成因还有争议。

平衡剖面和断层活动率表明，中新世期缺少正断层活动。

同时，正断层活动的结果只能导致加速沉降发生，与中新世的减速沉降因果关系较小。

因此，中新世的减速沉降应归因于其他构造因素。

图4后裂谷期平衡地质剖面
Fig.4The balance geological section at post-rift stages
4.3古热流与减速沉降
岩石圈减薄后的热收缩是造成裂后沉降的主要机制，由于热流减少，热收缩沉降的速度随时间呈指数下降。

因此，热流是裂后沉降的直接控制因素，加速和减速沉降的出现可能会反映热流的转变。

为了探究沉降变化与古热流之间的关系，本文设计了基于热流模型预测值和观测值的综合分析。

以平均应变总量β=1.75的剖面A 为例，使用改进的模拟热流方程对剖面A 的预测平均热流和沉降进行估算。

由于后裂谷期沉降是同裂谷期沉降的函数，可以据此调整模型获得热流的观测值：给从应变率反演中获得的β值增加一个Δ值（如0.05），计算出的沉降量也相应的增加。

当结果与观测的裂后沉降相吻合时，可认为此时的热流与沉降的观测值相对应。

该值以式5的左边表示。

⎥⎦⎤⎢⎣⎡--+∂∂+∂∂-=∆+∂∂r p a a z c A z T k z T H t T )(exp 022ρυ(5)
其中式5等号左侧是与观测的沉降相对应的热流。

六口井的数据处理完成后，可采用线性指数拟合的方法得到与沉降观测值相对应的热流曲线（图5a ）。

临清坳陷的古热流模型得出的数值为61~94mW/m 2，与本研究估计的热流值一致。

图5计算热流曲线Fig.5The heat flux curves
a 比例地质年代尺度为20Ma ，比例热流值以裂谷发育前的初始热流为1
b 示模型计算的理论热流值，红线表示根据实际沉降量反演的实际热流值。