渤海湾周缘高温地热异常特征及深部热结构分析

张 健, 何雨蓓, 姜程浩, 褚 伟, 方 桂

中国科学院大学地球与行星科学学院, 北京 100049

以大地热流为基础的地热学研究表明(姜光政等, 2016; Jiang et al., 2019), 渤海湾周缘地热资源十分丰富, 且高温地热资源潜力随深度的增大而增大。在渤海湾周缘开展高温地热异常特征与深部热结构形成机制研究, 对解决我国东部高温地热勘查开发瓶颈、优化能源结构体系具有重要的现实意义。地热能资源是我国“十四五”规划的未来战略性接替能源, 高温地热资源蕴含丰富的地热能, 是国家寻求非碳基能源的主要目标, 在国家加快构建现代能源体系、努力实现“双碳”目标任务进程中, 清洁高效、低碳排放的高温地热资源日益重要, 备受关注。高温地热异常特征与深部热结构的地热学研究, 除了需要大地热流资料外, 还需要震、磁、重、电等地球物理资料, 分析与地热能赋存密切相关的居里点深度、Vp相关的浅层生热率、Vs相关的壳幔低速层、泊松比等地球物理学指标。最近, 中国科学技术大学地震与地球内部物理实验室张海江教授课题组利用中国大陆地震台网体波走时数据、全国面波频散曲线数据集, 应用体波-面波联合成像算法, 获得了新的中国大陆岩石圈速度模型(Unified Seismic Tomography models for continental China lithosphere)USTClitho2.0(Han et al., 2022), 为地热学研究提供了区域动力学和深部结构参考模型。本文基于USTClitho2.0地震学模型(图1b), 以及航磁、区域重力资料, 利用地热学方法对渤海湾周缘高温地热异常特征及深部热结构开展分析。

图1 渤海湾周缘热流点分布及地震Vs波速模型Fig. 1 Heat flow points and Vs velocity model of Bohai Bay

1 构造背景与地热地质条件

渤海湾位于华北克拉通东部, 新生代以来经历了古近纪裂谷期(66—23 Ma)、新近纪—第四纪后裂谷期(23 Ma到现在)两个构造演化阶段(Ye et al.,1985; Li et al., 2012)。裂谷期, 由于构造沉降(Qi and Yang, 2010; Li et al., 2012), 发育一系列地堑和半地堑, 沉积层序包括孔店组、沙河街组和东营组(Yang and Xu, 2004; Hao et al., 2011; Zhu et al., 2014); 后裂谷期, 由于热沉降(Qi and Yang, 2010; Huang et al., 2014), 形成以渤中坳陷为沉降沉积中心、众多沉积坳陷和凸起相间的大型盆地, 沉积层序包括新近系馆陶组、明化镇组和第四系平原组。郯庐断裂带渤海段在新生代的断层活动控制了两侧坳陷的构造-热演化, 远离郯庐断裂的坳陷, 如, 黄骅坳陷(李明刚等, 2009), 所受影响较小, 昌潍坳陷、渤中坳陷处于断裂带上, 所受影响较大。镜质组反射率和磷灰石裂变径迹计算的热演化史表明(邱楠生等,2007), 昌潍坳陷地温梯度在孔店组沉积时期为58～57 ℃/km, 沙河街末期为 47 ℃/km, 东营组至馆陶组沉积末期为43 ℃/km, 现今为42 ℃/km。渤中坳陷地温梯度在孔店组沉积时期为56～50 ℃/km,沙河街组沉积时期为50～40 ℃/km, 东营组沉积末期为37.5 ℃/km, 馆陶组沉积末期为 34 ℃/km, 现今为33 ℃/km。总体上, 在从古近纪到新近纪的断陷向坳陷构造演化过程中, 地温梯度逐渐降低。

渤海湾发育有沉积盆地型、断裂裂隙型、岩浆活动型等多种地热资源类型, 是华北重要的地热资源区。2014年, 在渤海湾南岸的陈庄凸起, 利津干热岩SDGRY1井在2492 m测温104.2 ℃(谭现锋和王浩, 2015), 终孔孔深2 500.58 m, 测温曲线上段(0—1200 m)地温梯度 47.7 ℃/km, 下段(1200—2500 m)地温梯度29.0 ℃/km, 上、下两段均为线性传导型地温曲线。推测在5120 m深度以下, 陈庄凸起温度将大于 180 ℃(谭现锋和王浩, 2015)。2019年, 在渤海湾北岸的马头营凸起, 马头营干热岩M1井在3965 m钻遇151.25 ℃高温岩体(齐晓飞等, 2020), 钻孔地层: 0—250 m为第四系松散沉积物, 250—1150 m 为明化镇组砂岩、泥岩, 1150—1350 m为馆陶组砂岩、泥岩, 1350—4000 m为太古宙变质花岗岩类或变质岩系。测温曲线表明: 新生界地温梯度为50～70 ℃/km; 新生界底至2500 m太古宙地层, 地温梯度为 12 ℃/km; 太古宙地层内2500—3500 m 深度, 地温梯度为 17 ℃/km; 深度＞3500 m, 地温梯度为27 ℃/km。预测在5074 m深度以下, 马头营凸起温度将大于 180 ℃(蔺文静等,2021)。

大地热流点计算得到的热岩石圈厚度和莫霍面温度是研究深部高温热结构的重要依据。本研究区内目前共有实测大地热流数据91个(图1a), 热流值在 40.9～113.9 mW/m2之间, 平均值为 66.88 mW/m2。高热异常主要出现在冀中、渤中等坳陷区, 其中,冀中坳陷岩石生热率在1.05～1.61 μW/m3之间, 热导率在 1.88～3.5 W/(m·K)之间(姜光政等, 2016; Jiang et al., 2019; Wang et al., 2019)。前人(何丽娟等, 2001;臧绍先等, 2002; 彭波和邹华耀, 2013; Qiu et al.,2015)依据一维稳态热传导方程, 利用钻井或插值点的地表热流, 向下外推各深度的温度, 计算热岩石圈厚度、莫霍面温度, 认为渤海湾热岩石圈厚度在61～102 km之间, 莫霍面温度在450～850 ℃之间。由于钻井数据向深部外推会带来一定误差, 比如,地壳厚度3 km的误差就将会导致11%～14%的热岩石圈厚度误差、15%～20%的 Moho面温度误差(Wang and Cheng, 2012), 且一维计算中, 地壳结构横向起伏较大时, 小距离横向热传导也会对计算结果造成较大误差。因此, 人们更倾向于由航磁异常反演居里点深度约束地壳温度(张健等, 2018)、由地震波速模型计算上地幔温度并以 1300 ℃等温线确定热岩石圈厚度(安美建和石耀霖, 2007)。在居里温度统一取475 ℃条件下, 由航磁异常反演得到的渤海湾周缘地热异常区居里面深度在18.4～26.5 km之间, 平均22.1 km, Moho面温度在600～800 ℃之间(张健等, 2022)。依据华北地区S波速结构模型, 计算的渤海湾盆地平均岩石圈厚度为 100～110 km(杨嵩等, 2013)。

上述计算中, 存在2个问题: 1)居里面反演是针对整个中国东部地区, 居里温度的设置缺乏约束;2)岩石圈厚度计算采用的波速模型分辨率仅为3°～6°, 难以分辨小构造单元间的热结构差异。本文研究将针对这 2个问题做出改进: 1)将依据居里等温面深度变化范围, 利用大地热流值与地温梯度之间的关系, 约束居里温度; 2)将依据最新的0.5°×0.5°分辨率地震学模型 USTClitho2.0(Han et al.,2022)(图1b), 计算岩石圈上地幔温度和岩石圈底界面1300 ℃等温面深度。

2 计算方法

2.1 浅层(≤5 km)地温计算方法

式(1)中, 如果已知地壳结构和第i深度层热导率Ki、生热率Ai, 则由地表温度Ts、热流Qs可以计算不同深度Zi处的地温Tz、热流Qz。

2.2 居里点深度计算方法

居里点深度 CPD(Curie Point Depth)是表征铁磁性矿物消磁温度(居里温度)所对应的深度。实验室中, 各类铁磁矿物的消磁温度大致为: 磁黄铁矿300～350 ℃, 磁铁矿 575～585 ℃, 镍铁矿 760～800 ℃。当含有这些磁性矿物地层的温度接近消磁温度时, 磁性特征逐渐消失。由居里点深度勾绘的居里等温面(Curie Point Isotherm Surface)是控制地壳热结构的重要温度界面。

计算居里面深度的常用方法是磁异常谱分析方法, 设: 磁异常ΔT的功率谱径向平均为θΔT。

式(2)中,A、B、C为可选常数,kλ为波数,Zt、Zb分别是磁性体顶、底界面,Z0是磁性层的中间深度。在短波谱段(波长小于两倍磁性层厚度), 由磁异常功率谱的斜率可以估算Zt。在长波谱段,Z0可以根据拟合曲线的斜率求出。Zb为居里点深度 CPD(磁性体底界面深度)。

计算时, 磁异常的径向功率谱的“滑动窗口”为 0.5°×0.5°, 为保证重叠 50%, 滑动距离为窗口宽度的1/2。

2.3 中-深层(5 km≤Depth≤CPD)地温计算方法

稳态热传导条件下, 大地热流Q在数值上等于地温梯度dT/dZ与岩石热导率K的乘积。居里等温面的温度TC等于磁性层厚度DC与地温梯度 dT/dZ的乘积。因此, 利用式(3)计算的居里面深度DC和统计得到的居里温度TC, 可以计算地温梯度dT/dZ,结合式(1)计算的 5 km深度的温度T5km, 可以插值得到5 km ～ CPD居里点之间任意深度的温度。

2.4 地幔(50 km≤Depth≤250 km)地温计算方法

在 50～250 km 深度范围, 受温度影响的岩石非弹性是控制地震Vs波速的主要因素(Sobolev et al., 1996; Goes et al., 2000)。高温条件下, 利用品质因子Q的非弹性校正, 得到非弹性校正后温度相关的Vs波速计算公式为:

式(4)中,A、a为非弹性常数,ω为非弹性影响频率,H为活化能,V为活化体积,R为普适气体常数。采用非弹性模型a=0.15,A=0.148,H=500 kJ/mol,V=20 cm3/mol(Goes et al., 2000), 在给定矿物组成成分和温度压力条件下, 根据矿物弹性常数随温压变化关系、以及高温条件下的非弹性影响, 由式(4)正演计算剪切波速Vs。如果已知上地幔各深度剪切波速结构, 则在给定初始条件下, 通过反演迭代,计算波速与观测波速的差值ΔVs, 不断修正初始温度模型, 拟合ΔVs可得地幔三维温度场分布(Goes et al., 2000; 安美建和石耀霖, 2007)。

3 计算结果

在反演上地幔温度的式(4)计算中, 我们采用最新的地震学模型USTClitho2.0(Han et al., 2022)的Vs作为观测波速, 约束和拟合地幔温度模型。该模型给出了横向 0.5°×0.5°、纵向 0、5、10、15、20、30、60、80、100、120、150 km 的中国大陆 72°—136°E、18°—54°N 的三维地震纵波速度(Vp)和横波速度(Vs)结构, 考虑了地形起伏对浅部速度结构的影响, 反演的 1300 ℃绝热等温深度与地震学低速带吻合较好, 对于小构造单元间的热岩石圈厚度差异具有较好的分辨率。

此外, 由于马头营凸起 M1井为渤海湾北岸深层高温地热资源勘查提供了重要线索(张保建等,2020), 陈庄凸起利津 SDGRY1井为渤海湾南岸深层高温地热资源勘查提供了重要借鉴(蔺文静等,2021), 本文在开展深部高温热结构计算研究中, 均以马头营M-1井、利津SDGRY1井为参照。

3.1 浅层高温特征

由式(1)计算浅层地温、热流过程中, 渤海湾周缘地层厚度、热导率、生热率资料均取自前人文献(胡圣标等, 1999; Hu et al., 2001; 左银辉等, 2013;常健等, 2016; 刘琼颖和何丽娟, 2019; 王朱亭等,2019; Wang et al., 2019; 陈超强等, 2022), 如表1。

表1 渤海湾周缘地层热物性参数表Table 1 Thermophysical parameters of strata around Bohai Bay

渤海湾5 km深度地温分布如图5a所示。图5a中, 渤海湾周缘5 km深度的地温大于140 ℃, 其中,北岸的马头营M1井、南岸的利津SDGRY1井所在的构造凸起部位, 地温大于 160 ℃, 且马头营 M1井所在的凸起构造部位在5 km深度达到了前人预测的180 ℃高温(蔺文静等, 2021), 而利津SDGRY1井所在的构造凸起部位没有达到前人预测的180 ℃高温(谭现锋和王浩, 2015)。图5b、5c是计算地温曲线与实测地温曲线的对比, 其中, 图5b是研究区内 14个高热流点(≥80 mW/m2)的计算结果, 图中温度-深度曲线序号与图5a中括号内高热流计算点序号对应; 图5c是马头营 M1井、利津 SDGRY1井实测地温曲线和地温梯度(蔺文静等, 2021)。对比图5b、c可以看出, 在0—1500 m、0—100 ℃区域,图5b中10、11号高热流点的计算地温曲线与图5c中利津SDGRY1井实测地温曲线具有较好的对应。由于利津SDGRY1井测温曲线线性较好, 是未受地层水热活动扰动的传导型地温曲线, 直至终孔, 其地温曲线均与图5b中10、11号高热流点的计算地温曲线吻合。马头营M1井在1500 m以深, 受到冷水降温扰动, 地温梯度降低, 地温曲线向低温偏移,2 km以深, 其地温曲线与图5b中13号高热流点的计算地温曲线不吻合。

图2 渤海湾周缘地热异常区浅层地温特征Fig. 2 Shallow geothermal temperature of Bohai Bay area

3.2 中-深层热结构

3.2.1 上地壳生热率

地壳浅层热量的主要传输方式是传导, 按傅立叶热传导理论, 在热量传导过程中, 生热率是影响沉积盖层、地壳上层热结构的重要参数。浅表层的生热率可以通过钻孔岩心测试获得, 但沉积盖层底界以下的中-深层生热率常常需要外推。表1中, 沉积层生热率 1.61～2.82 μW/m3是实测结果, 上地壳生热率 1.05～1.26 μW/m3是外推参考值。一般地, 在地壳花岗岩浆垂向分异过程中, 放射性生热元素向上迁移、在浅层富集, 因此地表生热率较高, 深部生热率较低(Hasterok and Webb, 2017; Hasterok et al., 2018)。通过回归分析, 在50 MPa条件下, 生热率 A和地震纵波Vp遵循指数规律(Rybach and Buntebarth, 1982): LnA=16.5-2.74Vp(相关系数0.9),单位:A为μW/m3,Vp为km/s。利用USTClitho2.0速度模型, 提取上地壳Vp数据, 按式LnA=16.5-2.74Vp+Ln0.2 (Ln0.2代替回归公式的系统误差)计算渤海湾周缘 5～15 km深度的生热率结构如图3所示。

图3 渤海湾周缘地热异常区中-深层生热率分布特征Fig. 3 Heat generation of upper crust in the geothermal anomaly area around Bohai Bay

图3a、b、c分别是5 km、10 km、15 km深度生热率等值线平面图。5 km深度, 生热率在0.16～2.67 μW/m3之间, 平均值为 1.04 μW/m3, 低于表1中沉积盖层1.85 μW/m3平均生热率值。10 km深度, 生热率在 0.097～1.53 μW/m3之间, 平均值为0.55 μW/m3; 15 km 深度, 生热率在 0.077～0.44 μW/m3之间, 平均值为0.21 μW/m3。计算得到的中-深层生热率符合放射性生热元素向上迁移、在浅层富集的规律。据此, 图3d给出了渤海湾周缘地热异常区上地壳生热率“自下而上增强”的生热率模式。与表 1给出的生热率外推参考值相比, 由USTClitho2.0得到的上地壳生热率分布具有更丰富的地热信息。图3是马头营M1井、利津SDGRY1井所在位置之下的生热率变化, 可以看出, 二者总体变化一致, 生热率都随深度变小, 但马头营M1井上地壳生热率小于利津SDGRY1井上地壳生热率。

3.2.2 居里温度

居里温度Tc(Curie temperature)自1903年被皮埃尔·居里发现后, 被应用于许多领域。在地热学领域, 将航磁数据反演计算的地壳磁性层底面深度MLBD(Magnetic Layer Bottom Depth)作为居里点深度CPD(Curie Point Depth), 进而得到居里等温面深度DCT(Depth to Curie Temperature isotherm)。假设:稳定陆壳地温梯度 25～30 ℃/km、地层内磁铁矿居里温度 580 ℃, 二者相乘, 则可得陆壳居里点深度19.3～23.2 km。居里面深度DCT与地幔热源密切相关, 地幔热源浅, 则居里面浅; 地幔热源深, 则居里面深。居里点深度CPD可以依据航磁异常由式(2)计算得到, 计算结果如图4所示。

图4 居里深度与居里温度分析图Fig. 4 Curie depth and Curie temperature analysis diagram

图4a是用于计算地壳磁性层底面深度MLBD的航磁异常图。受沧东断裂带、郯庐断裂带影响,研究区磁异常总体走向为 NE向, 表现为正负相间、串珠状变化的磁异常区, ΔT异常在-225.4～373.3 nT之间, 平均值为5.7 nT。其中, 利津SDGRY1井处于高值异常区, 马头营M1井处于低值异常区。航磁ΔT异常是研究深部热结构的重要资料, 通过磁异常获取居里面深度, 提供地壳磁性特征消失的温度范围, 可以约束中、下地壳温度。图4b为居里等温面深度Dc计算结果, 其深度范围在16.45～23.56 km之间, 平均值为20.83 km。居里面深度随地幔热隆升或热沉降而起伏, 利津SDGRY1井的居里面深度较浅, 马头营M1井居里面深度较深。

居里等温面深度是认识中-深层地壳热结构的重要指标, 如果已知居里等温面深度Dc, 则由居里温度Tc可以推出居里面之上陆壳的地温梯度。虽然,与居里温度相关的居里深度可以由航磁异常计算得到, 但由于地壳不同矿物的磁性消磁温度不集中,居里温度无法严格限定。实际地壳内居里温度随岩石中磁性矿物成分、含量变化, 并不是确定值, 而是统计平均值, 我们通过地表大地热流, 利用选择和统计方法约束居里温度的不确定性。图4c是给定热导率(K=3.3 W/(m·K), 参照表1, 上地壳热导率K=3.5 W/(m·K))条件下, 由式(2)计算不同居里面深度Dc、不同地表热流Q条件下的Tc曲线, 在300～500 ℃温度曲线中, 最终选择Tc=420 ℃。图4d是利用式(1)计算的研究区91个热流点在图4b上各自居里深度的温度, 在 250～650 ℃之间统计, 得到Tc=420 ℃。通过这两种方法, 我们认为渤海湾周缘居里温度平均值为420 ℃。

3.2.3 岩石圈底界深度

岩石圈底界深度随上地幔热隆升或热沉降而起伏, 是衡量深部高温热结构的重要依据。最新的地热学研究得到渤海湾盆地的热岩石圈厚度约80 km(陈超强等, 2022), 面波成像得出的渤海湾盆地的地震岩石圈厚度在 60～70 km 之间(李孟奎等,2018), 表明渤海湾周缘热岩石圈与地震岩石圈厚度十分接近, 或者说, 该区岩石圈底界流变边界层(何丽娟, 2014)厚度很小, 因此, 我们可以利用0.5°×0.5°高分辨率 USTClitho2.0地震学模型,由式(4)计算岩石圈上地幔温度, 然后利用 1300 ℃等温面确定岩石圈厚度。计算结果如图5所示。

图5 渤海湾周缘岩石圈底界深度与温度分析图Fig. 5 Lithospheric bottom boundary depth and temperature around Bohai Bay

图5a是60 km深度的地温等值线平面图, 温度范围在1100～1250 ℃之间, 平均值为1184 ℃, 其中, 马头营M1井、利津SDGRY1井在同一个高温带, 马头营M1井在60 km深度的温度约为1220 ℃,利津 SDGRY1井在 60 km深度的温度约为1200 ℃。图5b是80 km深度的地温等值线平面图,温度范围在 1230～1450 ℃之间, 平均值为1329 ℃, 其中, 马头营M1井北侧是一个1400 ℃的高温区, 利津SDGRY1井处在一个NE-SW向的高温带的中间。马头营M1井在80 km深度的温度约为1330 ℃, 利津SDGRY1井在80 km深度的温度约为1360 ℃。图5c是马头营M1井、利津SDGRY1井所在位置上地幔垂向地温曲线, 可以看出, 75 km以浅, 马头营M1井上地幔温度大于利津SDGRY1井上地幔温度, 75 km以深, 马头营M1井上地幔温度小于利津SDGRY1井上地幔温度。利用图5a、b温度结果, 可以求取1300 ℃对应的深度图, 结果如图5d所示。图5d中, 1300 ℃所对应深度变化范围在66.3～97.5 km之间, 平均值为76.8 km, 此即为由USTClitho2.0地震学模型求得的热岩石圈底界厚度。图5d表明, 渤海湾周缘岩石圈西薄东厚, 岩石圈最厚处在横穿渤海湾的郯庐断裂带之东北侧, 表明渤海湾盆地岩石圈曾在华北克拉通破坏过程中发生过的自西向东的减薄作用仍然在进行中。马头营M1井和利津SDGRY1井之间, 也是一个岩石圈相对较厚的区域, 厚度大于80 km。马头营M1井、利津 SDGRY1井各自在不同的岩石圈减薄区, 其中,马头营M1井位于一个大的北东走向减薄区的边缘,利津 SDGRY1井位于一个小的北东走向减薄区中央, 利津SDGRY1井岩石圈底界深度略浅于马头营M1井岩石圈底界深度, 与图5c的结果一致。

4 讨论与分析

4.1 Moho面与壳幔动力学

莫霍面深度与温度对岩石圈结构、壳幔构造演化、深部动力学具有重要意义。Moho面是壳、幔不同物质组分和成分结构的分界面, 对应地震波速不连续面, 通常以波速梯度 dV/dZ或密度差Δρ来定义或标志莫霍界面, 是认识深部结构与构造特征的重要界面。前人研究表明(张成科等, 2002; Zheng et al., 2005; 郝天珧等, 2014), 渤海湾周缘莫霍面埋深介于28～37 km, 地壳厚度自西北部36～37 km至东部渤海海域 28 km, 呈现逐渐减薄的趋势。渤海湾盆地冀中坳陷的Moho面深度32～37 km, 黄骅坳陷为29～33 km, 至海域中的渤中坳陷Moho面深度仅为28 km。前人(左银辉等, 2013; 常健等, 2016; 刘琼颖和何丽娟, 2019; 王朱亭等, 2019; Wang et al.,2019; 陈超强等, 2022)计算的渤海湾盆地的莫霍面温度在482～732 ℃之间, 其中东部的渤中坳陷莫霍面温度大于660 ℃, 中、西部的黄骅坳陷莫霍面温度低于 580 ℃。本文利用大地热流值及USTClitho2.0地震学模型计算得到Moho面深度与温度结果以及动力学剖面(位置如图1a所示)特征如图6所示。

图6 渤海湾AB剖面Moho深度及动力学特征分析图(剖面位置见图1)Fig. 6 Moho depth and dynamic characteristics of section AB in Bohai Bay (see Fig. 1 for section location)

图6a是由USTClitho2.0模型的Vs、Vp计算的泊松比的剖面特征。泊松比σ也称为横向变形系数,σ越大, 构造区横向变形越大, 反之亦然。图中, 泊松比σ在0.2～0.3之间, Moho面处于高、低泊松比σ分界的梯度区, Moho面之下是一个以σ= 0.28为界的高泊松比凸起区, 表明深部动力学过程较强烈。浅部5～10 km深度也存在一个以σ= 0.27为界的相对高泊松比区, 表明在马头营M1井、利津SDGRY1井浅部构造活动也较活跃。图6b是AB剖面温度结构, 图中, 居里等温面Curie的温度为420 ℃, 岩石圈底界面 Lab的温度为 1300 ℃, Moho温度在600～800 ℃之间。Moho面不是温度界面, 没有明确的地热学意义, 但可以利用 Moho面温度判断大陆岩石圈挤压收缩的变形程度, 莫霍面温度Tm≥650 ℃, 热结构是岩石圈流变强度和变形方式的主控因素; 莫霍面温度Tm≤550 ℃, 岩性结构是控制岩石圈流变强度和变形方式的主导因素。图中可以看出, 利津SDGRY1井之下以及马头营M1井以北区域, 莫霍面温度较高(Tm＞750 ℃), 1300 ℃岩石圈底界面上隆, 表明岩石圈拉张减薄幅度较大,是渤海湾盆地新生代发生拉张中心自西向东、自南往北、向海域方向渐进式迁移的重要原因。由于渤海湾周缘在华北克拉通破坏过程中发生过大规模的岩石圈减薄, 尤其在新生代发生了明显的拉张中心迁移, Moho面与420 ℃居里等温面、1300 ℃岩石圈底界面之间的位置关系, 对于岩石圈收缩变形过程、拉张程度等具有明确的指示作用。

4.2 壳-幔热流比

壳-幔热流比是分析岩石圈热结构的重要参数。地表热流可以通过地温测井, 由稳态温度曲线和相应深度的岩石热导率计算得到, 地壳热流和地幔热流则可以由3种不同方法得到: 1)由地震Vp波速与生热率关系(Rybach and Buntebarth, 1982)计算得到;2)由地表热流和岩石生热率模型线性插值得到(臧绍先等, 2002); 3)由地壳分层模型“回剥法”计算得到(邱楠生等, 2017)。本文采用第三种方法计算渤海湾壳幔热流比, 计算时, 参考前人资料(胡圣标等,1999; Hu et al., 2001; 左银辉等, 2013; 常健等,2016; 邱楠生等, 2017; 刘琼颖和何丽娟, 2019; 王朱亭等, 2019; Wang et al., 2019; 陈超强等, 2022),在USTClitho2.0模型基础上, 设定地壳分层结构参数如表2所示。

表2 渤海湾地壳分层结构参数表Table 2 Crustal structure parameters of Bohai Bay

由表2参数, 利用式(1)计算的壳幔热流变化特征如图7(剖面位置如图1a所示)。图7a为地表热流插值结果, 图7b为地震Vs波速结构, 图7c是岩石圈热流结构。

图7 渤海湾壳幔热流结构分析图(AB剖面位置见图1)Fig. 7 Analysis diagram of lithospheric heat flow structure in Bohai Bay (see Fig. 1 for section AB)

图7a热流剖面上, 莫霍面热流QMoho、5 km深度热流Q5km、地表热流Q0km均按相同模式起伏, 其中, 利津SDGRY1井位于两个高热流峰值之间的低谷点, 马头营 M1井位于高热流峰值区内, 这个特点在图1a中的地表热流分布中也可以看到。图7b地震Vs波速结构剖面上, 岩石圈底界面Lab对应上地幔低速带的顶界面, 反映了岩石圈上地幔的活动状态。Moho面与 Lab之间是Vs高速区,高速区波速在4.4～4.5 km/s之间, 该高速区厚度变化与渤海湾周缘岩石圈减薄、地壳减薄的差异性、同步性、协调性关系密切, 是岩石圈热-流变学研究(何丽娟, 2014)的重点层位。图7c给出了渤海湾周缘目前 91个实测大地热流数据(图1a)壳-幔热流比例, 其中,Qc/Q=56.4%、Qm/Q=43.6%,Qc/Q＞Qm/Q,近似“热壳”结构, 其较“热”的上地壳层热流在地表总热流中占比 23.5%、在地壳中占比41.7%(=23.5/56.4), 为该区寻找和发现中-深层高温地热资源提供了重要依据。

图7c的结果与前人的结果不同。前人认为(邱楠生等, 2017)渤海湾地表热流中, 来自地幔的热流占 49%～62%, 属于“热幔”结构, 但前人也指出(Wang and Cheng, 2012)地壳分层结构和Moho面深度, 对岩石圈热结构计算结果有很大影响。在地表热流、地壳热导率和生热率已知条件下, “热壳”还是“热幔”取决于地壳结构模型精度。我们在最新的高分辨率地震模型USTClitho2.0基础上得到渤海湾周缘“热壳”结论, 与高分辨率的USTClitho2.0地壳结构的可靠性密切相关。

4.3 中-上地壳地温场特征

居里面是地壳内重要的温度控制界面, 中、上地壳地温场界面起伏形态与 Curie界面密切相关。居里面之上的地温场, 受中-浅层热导率、生热率控制, 一定热导率条件下, 地层生热率高, 则地表地热异常背景高。但一定生热率条件下, 地层热导率大, 该地层的地温反而低。图8是基于USTClitho2.0地壳结构模型, 由密度与Vp、生热率与Vp统计关系(Rybach and Buntebarth, 1982)计算的居里面之上温度剖面。

图8 渤海湾中-上地壳温度特征分析图(AB剖面位置见图1)Fig. 8 Analysis of temperature characteristics of the middle upper crust in Bohai Bay (see Fig. 1 for section AB)

由于渤海湾盆地在中生代时期便已发生了大规模的岩石圈减薄, 至新生代初期进一步减薄, 中-上地壳拉张中心以及浅部沉积盆地的沉积、沉降中心一致向东部海域迁移。这种迁移过程造成居里等温面起伏, 并影响中-上地壳地温场的整体形态, 地温横向起伏较明显。其中, 利津SDGRY1井位之下是中、上地壳地温起伏的低值点, 马头营 M1井位之下是中、上地壳地温起伏的高值点(图8a, b)。马头营M1井位之下的地温一直高于利津SDGRY1井,直到大约75 km深度或Lab之下, 利津SDGRY1井下的地温才逐渐高于马头营M1井(图2b, c, 图5c)。图8a是依据USTClitho2.0模型Vp速度计算的密度(Rybach and Buntebarth, 1982)和温度剖面, 该剖面上, 密度值ρ在 2.25～2.90 g/cm3之间, 总体形态如同盆地, 中间低、两侧高。图8b是依据Vp速度计算的生热率和温度剖面, 剖面上, 生热率A在0.10～1.55 μW/m3之间, 总体形态如同“牛舌”, 中间高, 两侧及下部低。该生热率分布形态, 在中、上地壳热传导过程中, 增高了利津SDGRY1井至马头营 M1井位之间的上地壳浅层地热背景, 使剖面中部浅表层地温高于两侧地壳地温。但由于中、上地壳热导率较高(表2), 地壳浅层高生热率的增温效果并不十分明显。

5 主要结论

基于新的中国大陆岩石圈速度模型USTClitho2.0, 本文对渤海湾周缘开展了地热学计算和综合分析, 主要结论如下:

(1)USTClitho2.0得到的上地壳生热率分布比地热学研究通常的生热率外推结果, 具有更丰富的信息。Vp-A关系表明, 马头营M1井、利津SDGRY1井的生热率随深度变小, 但马头营 M1井上地壳生热率小于利津SDGRY1井上地壳生热率。

(2)航磁异常反演表明, 渤海湾周缘居里等温面深度Dc在16.5～23.6 km之间, 平均值为20.8 km。利津SDGRY1井的居里面深度较浅, 马头营M1井居里面深度较深。91个热流点的居里点温度计算与统计结果表明, 渤海湾周缘居里温度在250～650 ℃之间, 平均值Tc=420 ℃。中、上地壳地温场界面起伏形态与 Curie界面密切相关, 居里等温面起伏导致中-上地壳地温场横向变化明显, 其中, 利津SDGRY1井位之下是中、上地壳地温起伏的低值点,马头营M1井位之下是中、上地壳地温起伏的高点。

(3)USTClitho2.0地震学模型及大地热流计算得到 Moho面深度与温度结果表明: 渤海湾周缘莫霍面埋深介于 28～37 km, 温度在600～800 ℃之间。岩石圈1300 ℃等温界面深度在66.3～97.5 km之间,平均值为76.8 km。参考USTClitho2.0模型设定的地壳分层结构, 计算的壳幔热流比表明: 渤海湾周缘壳、幔热流在大地热流值中占比Qc/Q=56.4%、Qm/Q=43.6%,Qc/Q＞Qm/Q, 近似“热壳”结构, 较“热”的上地壳层热流在地表总热流中占比23.5%,是寻找和发现中深层高温地热资源的重要依据。

致谢:衷心感谢审稿专家对本文提出的宝贵意见和编辑部主审及各位编辑的辛勤指导。

Acknowledgements:

This study was supported by National Key Research & Development Program of China (No.2021YFA0716002), National Natural Science Foundation of China (No. 42176052), and Chinese Academy of Sciences (No. XDB42020104).