珠江口盆地白云凹陷新生代构造-热演化模拟

胡杰， 龙祖烈， 黄玉平， 张林友， 胡迪， 王一波， 胡圣标*

1 中国科学院地质与地球物理研究所， 岩石圈演化国家重点实验室， 北京 100029 2 中国科学院地球科学研究院， 北京 100029 3 中国科学院大学， 北京 100049 4 中海石油深海分公司， 广东 深圳 518000 5 中国地质调查局水文地质环境地质调查中心， 河北 保定 071051 6 油气资源与勘探技术教育部重点实验室， 长江大学地球科学学院， 武汉 430100

0 引言

珠江口盆地是南海北部最大的中新生代沉积盆地，因为良好油气勘探前景，其新生代以来的构造-热演化过程一直备受中外学者关注.Nissen等(1995)利用纯剪—单剪模型，讨论了南海北部大陆边缘的拉张机制，提出纯剪切模型在解释南海大陆边缘大尺度张裂过程上优于单剪模型，而这两种端元的联合使用则可能会更好地解释南海大陆边缘的拉张变形(Nissen et al., 1995).Clift和Lin(2001)通过分析珠江口盆地与北部湾盆地30口钻井的构造沉降数据，指出南海北部岩石圈是不均一拉张的(地壳拉张强于地幔拉张)(Clift and Lin, 2001).Shi等(2005)和赵中贤等(2010)研究了珠江口盆地的异常沉降现象，认为异常沉降与17 Ma左右的短期高速沉积速率有关(Shi et al., 2005; 赵中贤等, 2010).张云帆等(2014)通过计算白云凹陷的拉张因子，推测在23.8和13.8 Ma两个时期，下地壳伸展及流动引起的盆地沉积沉降作用是控制陆架坡折发育位置的一个重要因素(张云帆等, 2014).He等(2017)通过对珠江口盆地的构造沉降研究识别出两期快速的沉降，分别对应南海洋脊跃迁和台湾的弧陆碰撞(He et al., 2017).但目前针对白云凹陷的构造-热演化模拟并不多，且前人的研究主要集中在构造沉降，对其基底热流史涉及较少，而基底热流对探讨盆地构造背景及指导油气勘探有重要意义.

构造-热演化模拟是一种岩石圈尺度的研究盆地热史的方法，不仅可以恢复盆地基底热流史，还可得到不同时期的应变速率及拉张系数，以探讨盆地演化的地球动力学背景.构造-热演化模拟能摆脱样品数量限制，对研究区不同区域进行对比.白云凹陷地层完整、研究程度高，大量钻井、地震剖面和热流数据为构造-热演化模拟提供了充足的基础数据(Dong et al., 2009; Shi et al., 2017; Tang et al., 2014) .

本文收集了白云凹陷最新的三维地震剖面、钻井岩性数据，选取6条地震解释剖面通过多期非瞬时伸展模型进行构造-热演化模拟，获得了白云凹陷新生代各裂陷幕的构造沉降，进而预测了盆地基底热流的演化史，在此基础上对比了白云凹陷各区域拉张强度和构造-热演化的时空差异性.研究成果可为探讨盆地构造演化及深部动力过程提供启示，并可指导油气勘探.

1 地质背景

珠江口盆地位于南海北部被动大陆边缘，其构造格局具有明显“东西分块，南北分带”特征，呈现典型“三隆两坳”地形(吴景富等, 2013; 张健和宋海斌, 2001).珠江口盆地的构造演化与南海扩张密切相关，总体上，经历了裂前期、裂陷期、断拗转换期、拗陷期4个发育阶段，各个阶段发生多幕构造运动，其中包括神狐运动、珠琼运动一幕、珠琼运动二幕、南海运动、白云运动和东沙运动，最终盆地由早期分散分布的地堑、半地堑相互连通、扩展形成现今地貌(Hou et al., 2019;吕宝凤 等, 2012; 廖杰等, 2011; 张功成, 2010; 赵阳慧, 2016).珠江口盆地现今大地热流平均 71.5±13.4 mW·m-2，地温梯度平均37.87±7.35 ℃·km-1，属于典型“热盆”，其中白云凹陷大地热流平均值为 76 mW·m-2(Jiang et al., 2016，2019; 唐晓音等, 2014).珠江口盆地的热流分布格局具有北低南高的特点，整体受地壳减薄控制，局部受到热流体及岩浆作用影响(唐晓音等, 2018; 饶春涛和李平鲁, 1991) .

图1 珠江口盆地构造图及相关钻井Fig.1 Geological map of PRMB and the boreholes

图2 珠江口盆地新生代地层柱状图，CR，MI和KO分别代表Cramer等(2009), Miller等(2011) and Kominz 等(2008)Fig.2 Generalized Cenozoic stratigraphic column. CR, MI and KO denote Cramer et al. (2009), Miller et al. (2011) and Kominz et al. (2008), respectively

白云凹陷位于珠江口盆地南部珠Ⅱ坳陷，水深大部分为500～1500 m，地壳厚度为18～28 km，从北到南厚度减薄(卫小冬等, 2011).凹陷西接云开低隆，东邻东沙隆起，北接番禺低隆，南连荔湾凹陷，整体为东西走向(唐晓音等, 2016; 朱伟林, 2010).白云凹陷中心新生代沉积厚度可达13 km, 是盆地的沉降中心及沉积中心.其中，古新世一渐新世早期为湖泊充填期；渐新世晚期(珠海组沉积时)的海侵事件沉积了海相砂泥岩；早中新世珠江组下部为海相沉积；珠江组沉积晚期的再次海侵形成了开阔浅海，此后以泥质沉积为主(能源等, 2013; 张功成, 2010; 赵阳慧, 2016).

2 模拟方法及基础数据

2.1 构造沉降恢复

本文使用回剥法(Backstripping)来恢复构造沉降.回剥技术是根据沉积压实原理法则，在保持地层骨架厚度不变的条件下，以盆地内地层分层为基础，按地质年龄从新到老把地层逐层剥去，从而恢复每个时代末所有沉积地层的形态及古厚度(Clift and Lin, 2001; Sclater and Christie, 1980).该方法需要进行四种校正：(1)脱压实校正；(2)沉积物负荷校正；(3)古水深、古海平面校正(Allen and Allen, 2013).

(1) 脱压实校正

在正常压实情况下，孔隙度和深度关系服从指数分布：

φ=φ0e-cz，

(1)

脱压实校正中假设岩石的体积分为两部分：一是骨架颗粒Vs，认为在压实过程中保持不变；二是孔隙体积Vh，随着沉积压实厚度呈指数降低，方程(2).

(2)

单位体积骨架占据的厚度表示为

(3)

方程式中，y1、y2：所回剥地层的顶界与底界；Vs：岩石骨架厚度；φ：地层的孔隙度.

按照假设，当岩石回剥到新的深度D′(y′2,y′1为上下层深度)时，其骨架厚度不变，仅孔隙度发生了膨胀，进一步展开求积分，即有

(4)

对上式采用数值求解的方式不断迭代，即可得到去压实的每一沉积时期内的地层沉积厚度和深度信息.

(2) 沉积物负载校正

根据均衡原理，沉积物的负载沉降量可以表述为

(5)

式中，S*为岩层去压实后的地层厚度，ρs为盆地沉积层的平均密度，ρm为地幔密度，ρw为水的密度.

(3) 古水深、古海平面校正

古水深类似于沉积物负载，同样会导致附加的沉降值，其大小表示为

(6)

其中，ΔSL：相对现今海平面，古海平面的高程.

经校正后的构造沉降量可表示：

(7)

其中，Y：构造沉降量；S′：压实校正后的沉积物厚度；ρm：上地幔密度；ρs：沉积物垂向平均密度；ρw：水的密度；Wd伸展时古水深.

本次研究采用中海油深圳提供的六口实钻井砂泥比值，同时参考珠江口盆地沉积相及基底岩性资料进行赋值(表1).

表1 白云凹陷地层岩性综合统计表Table 1 Statistical table of lithology in the Baiyun Sag

孔隙度随深度变化关系采用指数形式.不同岩性的孔隙度、压实系数使用He等(He et al., 2017).

表2 白云凹陷砂、泥岩初始孔隙度及压实系数Table 2 Initial porosities and compaction factors for sandstone and mudstone in Baiyun Sag

对于深水盆地而言，古水深估算是构造沉降中的一个决定性的因素.此次研究，我们先在白云凹陷各区域选取资料丰富的井作为控制井，确定其各时期古水深，再通过插值获得整个区域古水深.中海油深圳分公司根据测井资料和地震层序研究得到的沉积相系统图和钻井样品的古生物资料，是本次研究古水深的主要依据.一般来说，冲积-河流相沉积可忽略水深，滨-浅海水深20～200 m，上陆坡 200～500 m左右，陆坡至深海平原的浊积相沉积>500～2000 m(Allen and Allen, 2013; Xie et al., 2014; 龚再升等, 2004).裂陷期主要为陆相及洋陆过渡相，古生物丰富、古水深小且容易确定，对构造沉降影响也小.所以本次研究主要难点在裂后古水深.ODP 1148现今水深3294 m，Zhao(2005)指出古水深从32.8 Ma开始增长至今，中中新世(14 Ma)之后水深达到2500到3500 m, 超过现今水深的76%(Zhao, 2005).另外根据中海油深海分公司内部的最新研究成果，认为13.8 Ma以后陆架坡折一直维持在白云凹陷北侧，与今天的位置基本一致，白云凹陷处于陆坡深水环境(Zhang et al., 2013; 柳保军等, 2019).因此，我们断定白云凹陷古水深陡增发生在14 Ma之前，并结合沉积学资料和海底地形，确定了白云凹陷的水深演化.

2.2 伸展模型

本文基于多期次非瞬时伸展模型(Chen, 2014; He et al., 2001)，在模型中我们假设岩石圈水平方向上的伸展导致垂向上的减薄，热的软流圈物质上涌补偿缺失的岩石圈，温度场随之升高，伸展结束后由于热衰减导致盆地热沉降(Jackson et al., 2008; McKenzie, 1978).伸展过程中岩石圈的温度场主要受控于平流传热，即软流圈上涌，而在伸展结束后岩石圈温度场的演化依赖于传导传热作用，即热沉降(He et al., 2001; 张林友等, 2016; 袁玉松等, 2009).热传导方程采用质点网格有限差分法求解，相对以往的模型有了更高的精度，同时可以实时追踪模型中任何质点的温度、压力以及速度演化信息(Duretz and Gerya, 2013; Gerya, 2009).

岩石圈温度场通过求解瞬时热传导方程来获得，平流项由质点移动来实现(Gerya, 2009)，方程(8)可简化为方程(9)：

(8)

(9)

其中，v为速度，k为热导率，κ为热扩散系数，ρ为密度，Cp为比热，A为生热率.

在得到温度场信息后，盆地伸展导致的理论构造沉降量ys由艾瑞均衡获得：

(10)

(11)

在做反演计算时，通过不断地调整伸展速率建立不同的沉降曲线(Hmod)，并与回剥得到的构造沉降曲线(Hbackstripping)对比，寻找使两者最接近的伸展速率.本次研究使用misfit(所有对比年龄点Hmod与Hbackstripping差的绝对值)来描述其相似度，那么上述问题转化为对misfit进行最优化，求取其最小值对应的伸展速率，并以此计算热史.本文使用最小二乘序列二次规划SLSQP(Sequential Least Squares Programming)(Kraft, 1988)作为最优化算法来实现上述过程，其目标函数为

(12)

其中Hbackstripping为观察到的构造沉降史，Hmod为模拟得到的构造沉降史，i为对比年龄点，n为对比年龄点数目.

在此次研究中模型上边界为恒温(0 ℃)，模型下边界固定在90 km, 温度保持在1300 ℃，初始温度场根据稳态热传导方程计算得到.由于南海北部大陆边缘构造演化复杂，通过珠江口盆地钻遇的基底花岗岩样品的地球化学分析，揭示它们形成于一个安第斯型的大陆弧构造背景(Holloway, 1982).因此岩石圈拉张开始前的岩石圈厚度并不一定为前人普遍使用的125 km(McKenzie, 1978)，岩石圈极有可能在拉张开始前就经历过减薄.因为构造沉降对初始岩石圈厚度、初始地壳厚度相对敏感，在使用多组初始地壳及岩石圈厚度值进行试算之后，将沿各剖面的现今热流及地壳厚度 (OBS1993剖面，(Yan et al., 2001))作为约束，最终将初始地壳厚度设置为28 km，其中上地壳17 km，下地壳11 km，初始岩石圈厚度90 km.地壳深部热物性参数未有实测数据，岩石热物性(生热率和热导率)10%的偏差将造成总构造沉降及热流1%和3%的不确定(刘琼颖和何丽娟, 2019)，因此本次模拟热物性采用普适性参数.通过在热传导方程中提高热导率以等效模拟软流圈对流的热效应(刘琼颖和何丽娟, 2015)，通过多次试算，认为白云凹陷等效热导率设置为6 W·(m·k)-1比较合理，模型中使用的参数如表3.

表3 地壳、岩石圈和软流圈的初始厚度及模型中使用的热物性参数Table 3 Initial thickness of crust, lithosphere and asthenosphere and thermal property values and used in the modelling

2.3 模拟虚拟井选点

本次研究在白云凹陷共收集选取6条侧线(图3)，并识别出T30(5 Ma)，T32(10 Ma)，T35(13.82 Ma)，T40(15.97 Ma)，T50(19.1 Ma)，T60(23.03 Ma)，T70(33.9 Ma)，T80(38 Ma)，WCSB4(41 Ma)，Tg(47.8 Ma)地震反射界面.本次选点为达到更好的拟合效果，仅选择构造沉降较大人工井.同时为简化模拟计算过程，本次模拟未考虑剥蚀.主要选择根据中海油深海分公司提供的古水深以及岩性数据，共建立人工井 84口，并计算其构造沉降及热流史.

图3 白云凹陷剖面位置图，红点为人工井位置，等值线图代表白云凹陷基底深度(m)Fig.3 Location of modeling profile in Baiyun Sag, the red points represent the artificial wells, the contour map indicates the depth of Tg in (m)

3 模拟结果

3.1 构造沉降特征

通过对白云凹陷84口虚拟井的模拟，得出白云凹陷始新世以来主要存在两期快速沉降，并且构造沉降具有明显的东、西差异.两期快速沉降分别发生在始新世(47.8～33.9 Ma)与中新世(23.03～13.82 Ma).并且不同区域拉张因子各不相同，总体呈现北低南高，中间高、边缘低.白云西洼第一期拉张强，拉张因子约1.4～1.6；第二期拉张因子约1.2～1.3，总体拉张因子1.8～2.0.白云东洼第一期拉张弱，拉张因子1.1～1.2；第二期1.3～1.4，总体1.6左右.白云主洼中心拉张因子最高，最大可达到3.5以上，第一期1.6～1.7，第二期1.9～2.1，总体3左右.白云南洼第一期拉张弱，拉张因子1.1～1.3，第二期强1.6～1.7，总体2左右.

3.2 基底热流史及其特征

白云凹陷存在两期拉张过程，第一期拉张自始新世 (47.8～33.9 Ma)，强烈裂谷作用，其中白云西洼和主洼拉张强，凹陷中心古基底热流快速上升至～82 mW·m-2.第二期拉张中新世(23～13.8 Ma)，此次拉张在白云凹陷南部更强烈，白云凹陷主体在13.8 Ma达到最高古基底热流，其中白云主洼中心最高古基底热流～100 mW·m-2，此后进入冷却阶段，古基底热流值缓慢下降、逐渐接近稳态.白云凹陷都在13.8 Ma达到最高古热流，只是不同区域因为拉张强度不同导致热流值存在差异.白云凹陷南部现今超过80 mW·m-2的基底热流也一定程度反映第二期的强烈的拉张作用.

4 讨论

4.1 深部变形的机制

在伸展过程中，地壳变形受控于伸展强度、地壳温度场和地壳厚度(Choi et al., 2013; Lavier and Buck, 2002; McIntosh et al., 2014).Clift等(2002)认为白云凹陷为非均一拉张，地壳拉张强度大于地幔，且可能存在下地壳流(Clift et al., 2002).张云帆等(2014)认为白云凹陷的初始流变结果为减薄型地壳，下地壳的韧性伸展变形占总伸展变形更高的比例(张云帆等, 2014).但是前人在讨论白云凹陷伸展时常忽略高速层(HVL)的影响.Zhao等(2010)通过对比Vp/Vs，认为珠江口盆地下地壳的成分为镁铁质，高速层为辉长质(Gabbro)(Zhao et al., 2010).Hu等(2020)以此计算流变学强度，得出现今珠江口盆地下地壳主体处于脆性域，不太可能发生韧性变形，而高速层强度低且处于韧性域(Hu et al., 2020).如此看来，前人研究的下地壳韧性层有可能为高速层.Wan等(2017)研究表明洋陆转换带(OCT)的HVL是与伸展相关的新生代岩浆活动造成的(Wan et al., 2017).考虑岩浆作用的时间，现今的测年显示玄武岩火山活动幕在约56～38 Ma，早于南海海底扩张(～33 Ma)(Wang and Shu, 2012; 邹和平, 1993).由此第一期拉张过程极有可能与高速层形成同期.在伸展过程中，形成高速层的岩浆作用通过平流作用将深部热带入地壳，极大地改变了地壳的温度场.Moho面温度升高，下地壳可能落入韧性域，流动性增强，脆性上地壳的高角度正断层数量也减少.同时，上下地壳较大的流动性差异可能造成在拉张过程中上下地壳的解耦，既不均一拉张.而在第二期(～23 Ma)拉张开始前，白云凹陷下已经存在3～5 km厚度高速层.根据Hu等 (2020)的研究高速层现今强度弱易发生韧性流动，那么在第二期拉张时期，因为更大的深度，理论上处于比现今更高的温度，应该具有与现今相当甚至更高的流动性 (Hu et al., 2020; 施小斌等, 2000).那么“软弱”的高速层将极大地增加地壳中韧性厚度的比例，在裂后期当沉积速率很大时, 由于深部对地表过程即沉积物的加载的响应，韧性层在容纳形变的过程中会发生横向流动, 从而造成裂谷的进一步沉降(Li et al., 2019; Tetreault and Buiter, 2018).理论上，HVL的横向流动遵循体积守恒的地壳均衡，也会造成软流圈上拱，在我们的模型中，对计算结果影响主要来自于其物性参数与下地壳的差异，但如前所述10%的热物性差异对构造沉降和热流影响为1%和3%(刘琼颖和何丽娟, 2019)，因此其对基底热流的绝对值影响并不会很大，相对值和趋势是可靠的.

4.2 裂后异常沉降

白云凹陷新生代沉降特征表现为“北低南高，中间高、边缘低”的特征.“中间高，边缘低” 主要由始新世强烈裂谷作用造成的，裂谷作用从中央到边缘逐渐减弱.而“北低南高”的格局为裂后异常沉降导致，在南部深水区，存在更大的裂后沉降.白云凹陷中新世(约23～13.8 Ma)可以发现明显的快速沉降.此次沉降可以大致与前人的研究对应，包括Clift和Lin(2001)关于珠江口盆地浅水区研究中的12～14 Ma的event 3(Clift and Lin, 2001)；Shi等(2005)研究中18.5～15.1 Ma的短期快速沉降(Shi et al., 2005)；Liao等(2011)白云凹陷研究中17.5～13.8 Ma的快速沉积事件(Liao et al., 2011)；Chen(2014)白云凹陷中17.5～10.5 Ma的快速裂后沉积(Chen, 2014)；He等(2017)中的第一期裂后快速沉积(He et al., 2017).此次沉降由北到南逐渐增强，但并没有接受巨厚的沉积而是造成古水深的陡增.构造沉降误差主要来源于初始孔隙度、压实系数、古海平面和古水深矫正.其中古水深的不确定是造成白云凹陷构造沉降不确定的主要因素，特别是在白云凹陷南部(现今水深超过2000 m)(Higgs et al., 1999).因此不同学者对白云凹陷裂后沉降时间的差异极有可能来源于古水深的差异.关于裂后异常沉降的原因有多种，Clift和Lin(2001)和Chen(2014)认为是重启的轻度的伸展或裂谷作用(Chen, 2014; Clift and Lin, 2001)；Shi等(2005)使用岩浆侵位来解释其17 Ma的快速沉降(Shi et al., 2005)；Liao等(2011)推荐沉积负载造成的下地壳流来解释裂后期的异常沉降，因为岩浆侵位不足以造成如此规模的沉降(Liao et al., 2011)；He等(2017)将异常沉降与白云运动联系，提出坡折带的跃迁是造成白云凹陷异常沉降的原因，但并没有解释其动力学机制(He et al., 2017).Xie等(2006)提出地幔动力地形来解释此异常沉降，得出40 Ma以来造成约300 m的异常沉降，但也无法完全解释深水区的巨大沉降(Xie et al., 2006).地幔动力地形是由板片俯冲扰乱地幔流造成的(Bertelloni and Gurnis, 1997; Burgess et al., 1997)，白云凹陷中生代处于俯冲带边缘，所以其在新生代继承了中生代形成地幔动力地形(Zhou et al., 2006).而裂后异常沉降主要发生于～23 Ma到13.8 Ma，且未找到明显证据支持裂前形成的地幔动力地形在裂后短暂的时间窗口快速作用.因此，地幔动力地形不太可能是裂后异常沉降的主要原因.此次沉降起始时间～23 Ma与南海洋脊跃迁(24～26 Ma)相近(Briais et al., 1993)，终止时间(～13.8 Ma)与南海停止扩张时间(～15 Ma)接近(Taylor and Hayes, 1983).我们推测南海扩张极有可能是裂后异常沉降的驱动力，白云凹陷此时处于伸展背景，伸展作用是造成异常沉降的主要原因.同时结合上小结论述，裂陷期形成的高速层可产生“下地壳流”，其也可产生相当的异常沉降.综上，裂后的异常沉降可能是重启的伸展作用与“下地壳流”共同造成的.

图4 白云凹陷各剖面人工井的构造沉降曲线Fig.4 Tectonic subsidence evolution curves of artificial wells for the profiles in Baiyun Sag

图5 白云凹陷各剖面人工井的基底热流史,色标与图4对应Fig.5 Basal heat flow history of artificial wells for the profiles in Baiyun Sag, same color in Fig.4

5 结论

1)新生代期间，白云凹陷的总拉张因子为1.4～3.5.第一期拉张被巨厚的裂陷期沉积物记录，第二期拉张造成古水深陡增，达到1000～2000 m；

2)白云凹陷基底热流经历两次增高事件，新生代初期(47.8 Ma)初始热流为55 mW·m-2；第一期(47.8～33.9 Ma)增温后热流处于60～81 mW·m-2，平均～68 mW·m-2；第二期(23.8～13.82 Ma)达到最高古热流，热流区间73～100 mW·m-2，平均～80 mW·m-2；

3)白云凹陷拉张中心在始新世位于白云西洼和白云主洼，在中新世位于白云主洼和白云南洼；

4)裂后异常沉降可能由重启的伸展作用和高速层引起的“下地壳流”共同控制的.

致谢感谢三位审稿专家提出的建设性修改意见,研究过程中得到中国地质科学院地质力学所唐晓音老师、中海油深圳分公司杨学奇工程师意见及帮助，在此表示感谢！