华北克拉通位于欧亚大陆东侧边缘, 包含内蒙古南部、渤海湾、黄海北部等地区, 是现今最为古老的克拉通之一, 可划分为西部陆块(Western Block)、中部造山带(Trans NCO)和东部陆块(Eastern Block)三部分(Zhao et al., 2001).许多研究表明, 相对表现仍较为稳定的华北克拉通西部陆块, 其东部陆块自中生代以来发生了明显的破坏, 表现在地表热流相对其他典型稳定克拉通明显偏高, 岩石圈厚度强烈减薄, 并伴随有大规模的岩浆活动和强烈的岩石圈伸展(Menzies et al., 1993; Xu, 2001; Ren et al., 2002; He and Wang, 2003; Wu et al., 2005; Chen, 2010; Qi and Yang, 2010; Windley et al., 2010; 朱日祥等, 2011, 2012; Li et al., 2012, 2014; He, 2015).对华北克拉通破坏的机制提出有诸多的假说, 但对东部破坏严重而西部却仍较为稳定的原因仍不十分清楚, 是目前亟需深入探讨的问题(吴福元等, 2008; 朱日祥等, 2011; Tang et al., 2013).其中作为对克拉通的稳定/破坏具有重要影响的岩石圈热结构特征普遍受到人们的关注, 为此结合不同的方法对此展开了诸多研究(何丽娟等, 2001; 臧绍先和刘永刚, 2002; An and Shi, 2006; Wei et al., 2008; Huang and Xu, 2010; 汪洋和程素华, 2011; Sun et al., 2013; 杨嵩等, 2013).但华北克拉通东部二维精细岩石圈热结构的数据目前仍相对较少, 更缺少对华北克拉通东西陆块差异性破坏在热结构特征上表现的对比分析.

本文基于长度约850 km, 东西向横穿整个华北克拉通东部区域的文登—阿拉善左旗二维地震剖面的宁津—忻州段(Jia et al., 2014; 王帅军等, 2014; 刘志等, 2015)，对华北克拉通东部的岩石圈热结构进行模拟计算, 获得了二维空间上较为精细的岩石圈热结构特征, 弥补了克拉通东部二维空间上热模拟的不足, 并与西部鄂尔多斯盆地热结构进行对比, 分析热岩石圈与地震岩石圈厚度差异变化的特征, 探讨可能存在的原因及其对华北克拉通差异性破坏的启示.

空间直角坐标系下, 二维稳态热传导方程表述为：

(1)

k为热导率(W·m^-1·K^-1)；A为热源(文中指放射性生热，μW·m^-3)；x为横向距离(km)；z为垂直深度(km)；T为温度(K).

边界条件为：

(2)

z₀为模型表面；z_m为模型底部；x₀为左边界；x₁为右边界；T₀为表面温度；Q_m为底部热流.

取常年平均表面温度值T₀=15 ℃(http://www.ncdc.noaa.gov/)作为上部温度边界条件, 模型左、右边界为绝热边界条件, 底部热流边界为Q_m.根据地震剖面, 将二维有限元模型划分为五层, 分别为沉积层、上地壳基底、中地壳、下地壳和岩石圈地幔(图 1c), 给定初始模型厚度为180 km.在模拟计算中通过不断调整迭代底部热流Q_m来拟合地表热流值Q_c最终得到岩石圈的热结构特征.

图 1

Fig. 1

图 1 (a)华北克拉通区域图, 虚线为克拉通边界, 蓝色和红色实线分别为克拉通东/西部地震剖面线(据黄方等, 2015; 刘志等, 2015); (b)华北克拉通大地热流分布特征; (c)二维岩石圈结构剖面, a中红线(据刘志等, 2015) Fig. 1 (a) Map of NCC regions. Dashed lines outline the North China Craton (NCC). Solid lines show location of two-dimensional seismic profile (from Liu et al., 2015; Hang et al., 2015). (b) Surface heat flow map in the NCC. (c) Two-dimensional profile of lithospheric structure (from Liu et al., 2015)

在模拟计算中, 我们仅使用了最新统计可信度较高的A类和B类热流数据(胡圣标等, 2001; 王良书等, 2002; 龚育龄等, 2003; 王永新等, 2003)以保证数据的可靠性.统计分析显示华北克拉通地区热流值符合正态分布.整体平均热流61.9±12.5 mW·m^-2, 西部61.3±9.4 mW·m^-2, 东部地表热流值略微偏高平均63.7±13.4 mW·m^-2, 具有自西向东增加的趋势, 与网格加权等统计方法得到结果一致(Gong et al., 2011).为了更加直观地得出热流的空间变化特征，我们使用具有相对较强保真能力的克里金插值方法对华北克拉通地区的热流进行插值计算处理, 得到了区域上的热流变化特征(图 1b).在此基础上提取地震剖面上的热流数据, 获取的地表热流值变化范围在56.1~75.0 mW·m^-2.如此长距离和大范围的热流变化很难用地壳生热率的变化来解释, 并且华北克拉通东部作为一个独立完整的陆块, 生热率结构研究至今没有发现明显的不同(Gao et al., 1998; 迟清华和鄢明才, 1998), 故而这种差异性特征应该是受到来自深部地幔热流供应差异的影响.

岩石热物性参数的研究早在20世纪就已经开始(Birch and Clark, 1940; Lachenbruch, 1970; Rybach and Buntebarth, 1984), 近年来不断有新的发现(Gibert et al., 2003; Pertermann and Hofmeister, 2006; Merriman et al., 2013; Miao et al., 2014; Wang et al., 2014).不同的学者见解不同, 采用的热物性参数也略有所差异(Artemieva and Mooney, 2001; Merriman et al., 2013; He, 2014; Miao et al., 2014), 针对华北克拉通, 本文中采用热物性参数如表 1.

表 1(Table 1)

表 1 模型热物性参数及值 Table 1 Values of physical properties used in modeling 岩石圈结构 热导率(W·m-1·K-1) 生热率(μW·m-3) 沉积盖层 2.2 1.08 上地壳基底 3.0 1.25 中地壳 2.8 0.84 下地壳 2.0 0.32 岩石圈地幔 3.2 0.03 注：热导率参考迟清华和鄢明才, 1998；生热率参考Furlong and Chapman, 2013；Vilà et al., 2010；Liu et al., 2001；迟清华和鄢明才, 1998.

表 1 模型热物性参数及值 Table 1 Values of physical properties used in modeling

地表热流依据来源主要分为两部分, 一是来自于地壳放射性元素生热, 二是来自于深部地幔的热通量供给.岩石圈的整体热状态主要由深部热流值控制, 相对地壳浅部生热元素贡献的热通量, 深部地幔的热通量供给大小更能反映出其所在地区的深部构造活动性.一般而言来自于地幔的热供给比例越大该地区构造就越活跃.拟合结果显示沿宁津—忻州剖面自西向东地幔热流值(Moho面处垂向传导热流值)逐渐增大, 由西部约24 mW·m^-2增大到东部郯庐断裂带附近约44 mW·m^-2, 远高于鄂尔多斯盆地(20.5~24.5 mW·m^-2)(黄方等, 2015).东部地区热流壳幔比值(Q_c/Q_m)由西侧逐渐向东减小(图 2b)变化范围在1.61~0.70, 明显不同于西部较高的壳幔比值(1.84~1.51)(黄方等, 2015).相比典型稳定克拉通地壳热流占主导(Q_c>Q_m)、低地幔热流值( < 25 mW·m^-2)的特征(Jaupart et al., 1998; Rudnick et al., 1998; Artemieva and Mooney, 2001; Hasterok and Chapman, 2011; Hacker et al., 2015), 华北克拉通东部来自深部地幔热流的贡献较大.

图 2

Fig. 2

图 2 (a)热流空间变化特征；(b)热流壳幔比(WNCC数据来自黄方等, 2015).剖面位置见图 1a Fig. 2 (a)Spatial variation of heat flow; (b)Ratio of heat flow contribution of crust to mantle(data of WNCC come from Huang et al., 2015), profile position as shown in Fig.1a

温度场计算结果显示，华北克拉通东部地区温度梯度远高于南非等典型稳定克拉通(对应于地表热流约40 mW·m^-2)(图 3)且横向波动明显(图 4a).华北克拉通东部Moho面温度变化范围在560~710 ℃, 同一深度温度最高值出现在郯庐断裂带附近, 靠近太行山温度逐渐降低(图 4a).对比西部鄂尔多斯温度场结构, 华北东部同一深度温度明显较高, 在100 km深度处温度相差可达300 ℃以上(图 4).拟合得到的华北克拉通东部地区地温曲线与新生代捕虏体显示的深度(压力)-温度数据吻合较好(Huang and Xu, 2010), 表明模拟计算结果的合理性(图 3).

图 3

Fig. 3

图 3 华北克拉通东部岩石圈地温曲线及捕虏体深度-温度数据(捕虏体数据参考Huang and Xu, 2010,位置见图 4a) Fig. 3 Lithospheric geothermal curve and xenolith P-T data in eastern NCC (xenolith P-T data from Huang and Xu, 2010, position as shown in Fig. 4a)

图 4

Fig. 4

图 4 (a)宁津—忻州剖面岩石圈温度场分布, 星号和虚线分别为图 3捕虏体和地温曲线位置；(b)鄂尔多斯剖面岩石圈温度场分布(黄方等, 2015, 改).剖面位置见图 1a Fig. 4 Lithospheric temperature field of (a) Ningjin-Xinzhou and (b) Ordos profiles (modified from Huang et al., 2015). Asterisk and dash line: location of xenolith and geothermal curve as shown in Fig. 3, profile position as shown in Fig. 1a

在对岩石圈的研究中, 其厚度特征一致备受关注, 基于不同的原理与方法提出了不同意义上的岩石圈厚度, 如地震、热、弹性及电性等岩石圈厚度(Artemieva, 2009).实际上定义岩石圈厚度的许多参数都很大程度上依赖于温度场变化造成的岩石物性的“突变”, 因此热岩石圈厚度是最直接也是争议最少的定义方式, 并得到了广泛的应用(Pollack and Chapman, 1977; Artemieva, 2009, 2011).地热学上常取到达特定温度等温面(~1200/1300 ℃或0.8倍地幔岩熔点温度)的深度或地温梯度线向下延伸与地幔绝热线交点的深度作为热岩石圈厚度的底界, 温度范围一般介于1200~1400 ℃之间.本文中取T=1300 ℃等温线作为热岩石圈的底界(Artemieva and Mooney, 2001), 得到华北克拉通东部热岩石圈厚度变化范围在75~139 km(图 4a), 同时以T₁=1200+0.5z地幔绝热线约束得到的热岩石圈厚度上限为71~133 km, T₂=1300+0.4z约束得到的下限为78~149 km.相比于华北克拉通西部鄂尔多斯盆地(图 4b)在128~162 km之间的热岩石圈厚度变化(黄方等, 2015), 华北克拉通东部热岩石圈厚度明显较薄.计算结果与东部新生代玄武岩地幔捕虏体制约的岩石圈厚度( < 80 km)(Xu, 2001)吻合较好.对比古生代金伯利岩中地幔包裹体和金刚石中含矿物捕虏体制约的~200 km的岩石圈厚度(Menzies et al., 1993), 在东部发生了>100 km的减薄.

岩石圈热结构主要受控于深部热通量特征.造成华北克拉通东、西部岩石圈深部热通量差异的主要因素是软流圈地幔有效黏度的差异, 软流圈有效黏度的降低促使地幔对流更为活跃, 流变边界层减薄以提供更多的热量(He, 2014).通过计算我们得到华北克拉通东部岩石圈底部垂向传导热通量为19.3~44.6 mW·m^-2, 西部约10.1~25.1 mW·m^-2, 造成两者差异需要软流圈有效黏度变化可能达2~3个数量级(He, 2014).同时流变边界层的存在也是造成地震岩石圈和热岩石圈厚度差异的主要因素, 较小的流变边界层厚度对应着较小的地震/热岩石圈厚度差异(Artemieva, 2009; Wang, 2010; He, 2014).对比华北克拉通地震岩石圈与热岩石圈厚度, 发现两者均呈现出由西向东逐渐减薄的趋势, 但之间的差值存在变化, 在鄂尔多斯盆地最大差异可能达140 km(黄方等, 2015), 较小也有80~90 km, 西部汾渭地堑差异略有降低约64 km(何丽娟, 2014), 而华北克拉通东部两者差异最大不到40 km, 且向东逐渐减小, 靠近郯庐断裂带附近两者差异只有不到10 km.利用流变边界层厚度与软流圈有效黏度对数的线性相关关系(何丽娟, 2014), 近似估算得到造成华北克拉通东、西陆块热/地震岩石圈厚度差异达数十千米的变化需要软流圈上部有效黏度差异达三个数量级，与底部热通量差异估算得到的黏度差异一致.

软流圈地幔的黏度受成分、温度、应变速率、含水量等多种因素的影响.实验证明含水量的增加能够显著地降低岩石的黏度(Karato and Wu, 1993; Karato and Jung, 2003; Korenaga and Karato, 2008; Li et al., 2008; Karato, 2010; Wang, 2010; Burov, 2011).华北克拉通长时间以来受到太平洋板块西向俯冲的影响(朱日祥等, 2012)，俯冲脱流体作用可以有效地增加软流圈地幔的含水量, 含水量的增加可能是导致华北克拉通东部被破坏的一个重要因素(Zhu and Zheng, 2009; Windley et al., 2010).若仅考虑流体作用的影响, 由经验蠕变公式：η∝C_OH^r-1, 取r=1.93(Korenaga and Karato, 2008), C_OH含量增加约11~35倍就能使软流圈有效黏度降低2~3个数量级, 加速软流圈对流, 压缩流变边界层厚度, 为克拉通东部岩石圈底侵至现今厚度提供所需的能量.然而考虑到现今华北克拉通东、西陆块存在热结构上的差异, 尤其是在西部减薄较为明显的郯庐断裂带附近, 岩石圈温度场较高软流圈上涌诱发部分熔融加上触发其他流变机制作用, 可能仅需要极少量的水(~200H/10⁶Si)软流圈顶部有效黏度就可以降低至10¹⁷ Pa·s(Wang, 2010), 足以保持现今克拉通东部陆块高地幔热通量、低地震/热岩石圈厚度差异的特征.

通过二维稳态热模拟研究获得了华北克拉通东部宁津—忻州二维空间上较为详细的岩石圈热结构特征.对比华北克拉通西部的低地幔热流值(21.2~24.5 mW·m^-2)、高壳幔热流比(1.51~1.84)以及厚的热岩石圈厚度(128~164 km), 华北克拉通东部呈现出高地幔热流值(24~44 mW·m^-2)、低壳幔热流比(1.61~0.70)和薄热岩石圈厚度的热结构特征, 热岩石圈厚度由西约139 km向东逐渐减薄到郯庐断裂带附近约75 km, 差异达64 km.

对比地震岩石圈与热岩石圈厚度差异, 由克拉通西部向东部逐渐降低, 这种差异性特征与岩石圈地幔底部垂向传导热流的变化趋势一致, 是由于软流圈黏度变化造成的, 初步估算有效黏度的差异达2~3个数量级, 若仅由水含量的不同来解释, 可能需要11~35倍的差异.考虑到现今的岩石圈热结构特征, 华北克拉通东部可能并不需要巨量的水来降低软流圈黏度以维持目前较高的地幔热流和较低的地震/热岩石圈厚度差异.降低的软流圈黏度增大了岩石圈底部热通量的供应, 为克拉通破坏提供了能量来源, 可能是华北克拉通东部遭受破坏而西部保持相对稳定的原因之一.

感谢审稿专家提出的宝贵意见.