鲁东招远地热田地热通量及地热成因研究

钟振楠，康凤新，宋明忠，郎旭娟，柳禄湧，傅朋远，李志杰

1) 山东省地质矿产勘查开发局第六地质大队，山东威海，264209；2) 山东省地质矿产勘查开发局，济南，250013；3) 河北地质大学水资源与环境学院，石家庄，050031；4) 山东省深部金矿探测大数据应用开发工程实验室，山东威海，264209；5) 河北省水资源可持续利用与产业结构优化协同创新中心，石家庄，050031；6) 河北省水资源可持续利用与开发重点实验室，石家庄，050031；7) 河北省水文学及水资源重点学科，石家庄，050031

内容提要: 鲁东地热区在地质构造单元上位于沂沭断裂带昌邑—大店断裂以东，地热资源丰富。本文采集了鲁东地热区招远地热田内一眼2000 m深地热井(DRZK01)中的40块岩芯样品进行岩石热导率、岩石生热率测试及分析，结合测温资料及收集资料对该区地热通量构成及分层热流进行了分析研究；采集区内典型地热流体样品进行水化学分析并采用合适的地热温标估算了该区热储温度；综合研究成果建立了该区地热成因概念模型。研究结果显示，该区岩石热导率数值较高，测量值在2.8～5.7 W/(m·K)之间，普遍高于上地壳平均热导率，地温梯度较高，为31.8℃/km；利用热导率和地温梯度计算的地热通量102 mW/m2中热传导分量为(73.2±6.18) mW/m2，对流分量为(28.76±6.18) mW/m2，其中热传导分量中地壳热流为22.5 mW/m2，地幔热流为(50.74±6.18) mW/m2，二者比值为1∶1.98～1∶2.52，为“热幔冷壳”型热结构。石英温标计算热储平均温度为128.6 ℃，热循环深度约3634 m。研究结果丰富了该区地热系统理论并为该区地热资源开发利用提供一定的理论支撑。

在全球气候变化的大背景下，地热资源由于具备清洁、方便、蕴藏量大且可再生等优点而日益受到世界各国重视(王贵玲等，2000；庞忠和等，2012,2014；袁利娟等，2020; 李泓泉等，2020； 原若溪等,2021)。招远地热田位于鲁东地热区，鲁东地热区是山地陆块地壳演化最为复杂的地区，尤其是中生代以来，地壳活动尤为频繁剧烈，地热资源丰富(徐希强等，2015)。为了科学高效利用该区地热资源，查明该区“源、通、盖、储”因素，开展地热资源赋存成因机理研究是不可缺少的关键环节(金秉福等，2000；邱楠生，2001；张涛，2011；史猛等，2019；赵辉等，2019)。大地热流、地温梯度、岩石热导率、岩石生热率作为基本的地热学热参数，为地球动力学演化以及油气生成研究提供地热背景和基础参数(邱楠生，1998；叶正仁等，2001；蔺文静等，2016；孙旭东，2020)；李学伦等(1997)利用地下水中SiO2含量计算山东半岛硅热流值，认为胶北隆起在大地构造上属于较高热异常区。Jiang Guangzheng等(2016)依托一口位于胶东半岛莱州湾附近的金矿勘探井的相关资料，计算得到一高质量大地热流数据。赵辉等(2019)对胶东地热田地热流体的补径排特征及地热田的地热资源进行了评价研究。前人研究为本文研究提供了丰富的资料支撑但是同时仍然存在很多问题，如该区岩石圈热结构方面的研究基础比较薄弱，地热热量构成及地热成因等方面尚存在很多需要深入研究的问题。因此本文在前人研究的基础上，以招远地热田内一眼2000 m深地热井为研究对象，采集岩芯样品、典型地热流体样品、收集相关数据，分析岩层热导率、生热率及地温分布特征，计算该区地热通量；估算热储温度及热循环深度，并建立了招远地热田地热成因概念模型，丰富了该区地热系统的相关理论，为该区地热资源可持续、高效开发利用提供理论支撑。

1 研究区地质背景

鲁东地热区以牟平—即墨断裂为界分为东西两大构造单元，招远地热田位于鲁东地热区西部华北板块胶辽隆起区的胶北隆起北部。鲁东地热区主要出露太古代、元古代、中生代、新生代的地层(图1)(崔煜烽等，2018)。该区中生代印支期由于地壳发生强烈的脆性张裂变形作用，形成了一系列近SN向张性断裂以及NE向、NW向和近EW向扭性断裂，并导致下部地壳部分岩石熔融，从而产生大规模的花岗质岩浆，沿断裂入侵形成花岗岩体，有的沿NE向断裂入侵形成一些浅成脉岩，由于上述地质演变，区域上具有良好的地热形成条件(田禹，2015)。

2 样品采集与测试

本次研究中，岩芯样品采集工作在招远地热田一眼深2000 m，地理坐标为120°24′41″，37°21′38″的地热钻孔(DRZK01)中进行(地热井具体位置见图1)，50 m取芯一次，共采集了40块岩芯样品，所取岩芯样品岩性为中生代玲珑序列大庄子单元的二长花岗岩(岩性柱状图见图2)。岩芯热导率数据是由中国科学院地质与地球物理研究所岩石热物性与地热测量实验室检测，使用德国生产的TCS(Thermal Conductivity Scanning)热导率自动扫描仪，其测量范围为0.2～25 W/(m·K)，测量精度为±3%。样品规格依照实验室要求进行样品预处理，样品长度大于3 cm，至少保留了一个平整面，平整面上没有裂隙、空洞，并且未遭受蚀变。岩芯放射性元素(铀、钍、钾含量)数据是由澳实分析检测(广州)有限公司检测，其中铀、钍元素含量采用ICP-MS(电感耦合等离子体质谱仪)分析，分析仪器是由美国生产的Agilent(型号7700X)，测试范围0.05×10-6～1000×10-6，精度控制相对偏差<10%，准确度控制相对误差<10%。钾含量采用X射线荧光光谱仪分析，分析仪器是荷兰生产的PANalytical(型号PW2424)，测试范围0.01～20%，精度控制相对偏差<5%，准确度控制相对误差<2%。放射性含量样品制样方式相同，样品烘干后，破碎过筛，保留直径1～2 mm样品300 g做正样，其余样品为副样，副样留存备用。将正样用无污染钵振动碾磨至200目；取30 g正样加入偏硼酸锂—四硼酸锂熔剂，混合均匀，于熔炉中在1025 ℃熔化。熔液冷却后，用包含了硝酸、盐酸和氢氟酸的混酸消解定容，再用等离子体质谱仪定量分析。另取30 g正样，均分为两份：一份试样中加入含硝酸锂的硼酸锂—硝酸锂熔融助熔剂，充分混合后，高温熔融。熔融物倒入铂金模子形成扁平玻璃片后，再用X射线荧光光谱仪分析。同时另一份试样放入马弗炉中，于1000 ℃灼烧。冷却后称重。样品加热前后的重量差即是烧失量。本方法所测的钾元素含量为烧失量的结果和XRF测得的元素氧化物结果相加。检测数据真实可靠。

水质全分析样品用500 mL塑料瓶采集并由中国地质科学院水文地质环境地质研究所国土资源部地下水矿泉水及环境监测中心进行检测，水质全分析是按照中华人民共和国国家标准引用天然矿泉水检验方法(GB/T8538-2008)检测。

3 钻井温度测量

本次研究中井温测量工作使用Server6000便携式数控测井系统完成，测量范围0～230℃，终孔后静井18 d 测温，可保证井温恢复至稳态或似稳态。并且所测温度及对应深度都经过Matrix Logging System处理后获取，温度测量精确度为±0.1 ℃。野外测温时，探头下降速率约为6 m/min，探头每隔0.1 m记录一个温度数据，抽析间隔2 m。

钻孔测温曲线及地温梯度曲线见图2，由图2可以看出，测温曲线整体为“上凸”型，说明该钻孔温度受到区域上升流的影响，地下水作为载体将热流携带至地表浅部(王一波等，2019)，因此直接利用钻孔中岩石热导率数据和地温梯度数据通过公式

图2 招远地热田钻孔柱状图及测温曲线图Fig. 2 Drilling column diagram and temperature measurement curve of geothermal well in the Zhaoyuan geothermal field

q=-K(dT/dH)

的计算值不能称为大地热流，可称为地热通量。由图2看出，该区地温增长分为5个层段：① 0～500 m层段地温增长迅速，平均地温梯度为159.3 ℃/km，地温梯度过大，这可能是由地表温度影响强烈原因造成的，因此进行温度场分析时，为了减小误差地层表层的温度数据一般会被舍弃(胡圣标等，1994；Huang Shaopeng et al.，2000；Pollack et al.，2003)；② 500～946 m层段温度随深度平稳增加，表明热量传递主要以热传导方式进行，地温梯度在该层段较为稳定，平均值为54.7 ℃/km；③ 1300～1800 m层段由于受到上升流影响，地温梯度基本不变；④ 946～1171 m层段和1300～1809 m层段地温梯度出现负增长，可能是由于局部破碎导致冷水混入造成的；⑤ 1171～1300 m层段温度比较稳定，当温度曲线出现分段现象时，地温梯度可采用上下两段的温度曲线来进行校正(徐明等，2011)。综合上文对于地温梯度的分段分析，该区地温梯度由500～1300 m层段的测温数据计算获取为31.8 ℃/km。

4 岩石热导率特征

4.1 岩石热导率测试数据分布特征

岩石热导率是表征岩石传热特性的物理量，它是研究岩石圈热结构和地球深部热状态的重要参数之一(倪守斌等，1999；栾锡武等，2003；章邦桐等，2010)。本次研究所取岩芯热导率测量值在2.8～5.7 W/(m·K)之间，平均值为4.0 W/(m·K)(图3)。

图3 招远地热田钻孔热导率值统计直方图Fig. 3 Statistical histogram of thermal conductivity in the Zhaoyuan geothermal field

4.2 岩石热导率测试数据的校正

岩石本身特性(岩石矿物成分、颗粒粒度、孔隙度、渗透率等)、温度、压力等因素会影响岩石热导率数值，对于花岗岩等致密岩石，影响其热导率的因素主要是温度(Brigaud et al.，1989)。由于本次研究中取到的岩石样品均为致密的二长花岗岩，岩石本身特性较一致，因此仅对岩石热导率测试数据进行了温度校正，校正过程是根据Sass等通过实验得出的经验公式进行的(Sass et al.，1992)：

(1)

(2)

式中，λ(0)，λ(25)分别为0℃和25℃时岩石的热导率，T为岩石原位形成温度，λ(T)为其原位温度对应下的岩石热导率，即为校正后岩石热导率。按照上述方法对研究区测量的40个岩石热导率数据进行了校正，校正后的岩石热导率数据见图4。由图4b可以看出，本次研究获得的岩芯热导率值可以分为三部分(图4b中分别由三个椭圆圈出)，每一部分热导率值基本都呈现出随着深度增加而增大趋势；而且相同岩性的岩芯热导率数值存在一定变化范围，这可能是由相同岩性之间物质成分存在差异造成的，也可能是受物质成分分布样式等因素影响(崔景伟等，2019)。校正后热导率数值分布在2.5～5.3 W/(m·K)之间，均值为3.5 W/(m·K)，所取岩石样品热导率测量值都高于上地壳平均热导率2.5 W/(m·K)，反映出该区具有良好的地热地质背景。

图4 招远地热田校正后钻孔岩芯热导率统计直方图(a)及热导率与深度关系图(b)Fig. 4 Statistical histogram of thermal conductivity after correctionand relationship between thermal conductivity and depth of the Zhaoyuan geothermal field

5 岩石放射性生热率特征

岩石的放射性生热率(A)是单位体积岩石中所含放射性元素在单位时间由衰变所释放的能量，单位为μW/m3，它表征着岩石自身生热能力的高低。岩石放射性生热率是描述地球内部热状态的一个非常重要的热物性参数，可以为研究地球深部热分布状态、热流的构成、探讨地球动力学过程提供重要的数据支撑(胡圣标等，1994，2001；何丽娟等，2006；Tang Boning et al.，2019)。

5.1 研究区地层浅部岩石放射性生热率特征

岩石中所含的放射性元素很多，但是对生热产生贡献的元素主要是U、Th、K 三种元素，岩石放射性生热率(A)就是通过测量计算岩石中U、Th、K 这3种元素的含量获取(李文庆，2015；郎旭娟，2016)。本文岩石放射性生热率数据是按照Rybach于1976年提出的计算公式获取(Rybach，1976；邱楠生，2002)：

A=0.01ρ(9.52CU+2.56CTh+3.48CK)

(3)

其中，A为岩石放射性生热率(μW/m3)，ρ为岩石密度(g/cm3)，CU、CTh、CK分别为岩石中铀(μg/g)、钍(μg/g)、钾的含量(%)。

本次研究中岩石放射性生热率计算结果见图5。岩石生热率数值统计直方图(图5a)表明，该区岩石生热率值分布在0.4～2.2μW/m3之间，计算数据多数集中在1.2～1.6μW/m3之间，平均值为1.3μW/m3。图5b显示，生热率数值和深度之间存在弱线性关系，生热率计算结果总体上随着深度增加而减小。

图5 招远地热田钻孔生热率统计直方图(a)及钻孔生热率随深度变化图(b)Fig. 5 Statistical histogram of heat production and relationship between heat production and depth of the Zhaoyuan geothermal field

5.2 研究区地层深部岩石放射性生热率特征

地层深部的放射性生热率的分布是解决该区地壳上地幔热流配分及壳内热结构、热状态等问题的关键(赵平，1995)。由于实际情况限制，本次研究中只取到了该区0～2000m层段的岩石样品，对于2000m以下地层深部岩石放射性生热率数据利用Rybach通过实验得出的VP—A关系公式获取(Rybach et al.，1984；赵平，1995；刘峰等，2020)：

50 MPa压力条件下：

(4)

100 MPa压力条件下：

(5)

200 MPa压力条件下：

(6)

其中，A为岩石放射性生热率，VP为地震波速。

利用上述公式，代入收集的研究区地震波速数据(图6)，得出区内5000～36000 m(36000 m深度约为该区莫霍面深度(方宝明，2006))层段的岩石生热率数据，详细数据见表1。

图6 华北克拉通超长观测距探测剖面不同块体二维速度结构图(据王帅军等修改，2014) Fig. 6 Two-dimensional velocity structure diagram of different blocks in the North China Craton ultra-long observation range detection profile(after Wang Shuaijun et al, 2014&)

由表1可以看出，利用VP—A关系公式计算的岩石生热率值随着深度的增加而递减，对比0～2 km层段实测生热率曲线(表1)，可以看出公式5和公式6计算的结果偏大，利用公式4计算的生热率结果最接近实际测量值。因此将公式4(50 MPa压力条件下)的计算结果作为该区深部岩层生热率结果。

表1 招远地热田 5.0～36.0 km(莫霍面)生热率计算结果及0～2 km实测值

6 地热田热储温度及热循环深度估算

6.1 地热田热储温度估算

地热流体中，某些决定热水中溶解组分比例的水岩反应依赖于温度，因此可用这些组分的比值来估算热储温度，这种热储温度估算方法称为地热温标法，常用的地热温标有二氧化硅地热温标计、Na—K地热温标计、Na—K—Ca地热温标计等(Fournier，1977；Arnorsson，1985；刘昭等，2014)。根据研究区地热流体的水化学特征及地热地质情况，结合各地热温标适用条件，选取了合适的地热温标对该区6口地热井进行热储温度估算，估算结果见表2。

表2 招远地热田深部热储温度估算结果Table 2 Estimated results of reservoir temperature in the Zhaoyuan thermal field

由计算结果可以看出，在多种地热温标中，玉髓地热温标计算结果偏低，因为钻孔DRZK01井深1000 m以上时，实测温度达120 ℃左右，故玉髓温标计算结果不是热储温度的真实反映；Na—K及Na—K—Ca地热温标计算出的结果偏高，石英温标和Na—K—Ca温标计算结果相差较大，且Na—K—Ca温标计算结果偏高，远高于石英温标计算结果，而石英温标计算结果高于K—Mg温标计算结果，故最终结果选用石英温标计算的热储温度(102.88℃～148.84℃)为该区热储温度的真实反映，均值为128.63℃。

6.2 地热田热循环深度

地热田内地热流体补给来源主要为大气降水(Sanliyuksel et al.，2011；Dotsika，2012)，大气降水入渗后经过深循环上升至地表形成地热水，循环深度越大，地热水温度越高。地下热水循环深度可由下列公式推算出来：

(7)

其中，H为地下热水循环深度，m；t1为地下热水深部最高温度，℃，取上文计算的热储温度均值128.63 ℃，t2为该区恒温带温度，取当地常年平均气温14 ℃；h为恒温度厚度，m，取值为30 m；I为地温梯度，取上文分析的0.0318 ℃/m。将以上数值代入公式(7)中得出该区地下热水循环深度约为3634.7m。

7 招远地热田地热通量构成及地热成因分析

7.1 地热田地热通量及构成分析

7.1.1地热田地热通量中对流分量

该区地热通量利用公式

来计算，q为地表热流(mW/m2)，k为岩石热导率[W/(m·K)]，dT/dH为地温梯度(℃/km)，负号表示热流传导方向。式中地温梯度取值为上文利用钻井(DRZK01)500～1300 m 层段测温曲线分析的31.8 ℃/km，岩石热导率数据采用钻井(DRZK01)500～1300 m层段岩芯热导率检测平均值3.2 W/(m·K)，最终利用公式q=-k(dT/dH)计算该区地热通量约为 102 mW/m2，具有较高的地热背景。由于该区内未开展过系统的大地热流测量工作，该区大地热流值可取利用距离该地热田约40 km的莱州三山岛黄金科钻的数据计算的热流值 73.24±6.18 mW/m2(Jiang Guangzheng et al.，2016)，由此可推算出本次研究中计算的地热通量102 mW/m2中对流分量为(28.76±8.76) mW/m2。

7.1.2地热田地热通量中传导分量

7.1.1 节中采用的大地热流值73.24±6.18 mW/m2即为该区地热通量中的传导分量，即在地表观测到的大地热流，主要由地壳中放射性元素衰变产生热能(地壳热流qc)和地幔深部热能(地幔热流qm)两部分组成(Jaupart et al.，2007；Czechowski et al.，2012；姜光政等，2016；刘绍文等，2017)。地壳中放射性元素衰变产生的热能

qc=Σqa

式中：qa=AD， 为地壳中各结构层由放射性元素衰变产生的热流，A为生热率，D为相应层段厚度。地幔热流由公式

qm=q-qc

获取(Birch et al.，1968；Roy et al.，1968；徐青等，1992；Seipold et al.，1998 ；Ranalli，2005；Duchkov et al.，2009)。进而可获得一个地区地壳热流与地幔热流配分比，分析其岩石圈热结构类型(Huang Shaopeng et al.，1992)。

本次研究中，0～2 km层段岩石生热率数值采取本文实际测算的岩石生热率均值1.3 μW/m3；2 km以下每个层段的生热率数值用表1中利用公式4(50 MPa压力条件下)计算的结果，取值为各层段上下层面生热率的均值。利用公式分别计算每一层段的地壳热流值、地幔热流值和热流值，具体计算结果见表3。

表3 招远地壳热结构计算结果表Table 3 The results of thermal structure of the lithosphere in Zhaoyuan

由计算结果可以看出，招远地区放射性元素衰变产生的热能和为22.5 mW/m2，即地壳热流为22.5 mW/m2，区内大地热流采用73.24±6.18 mW/m2，就可估算出区内地幔热流为(50.74±6.18) mW/m2，进一步估算出区内壳、幔热流比为1∶1.98～1∶2.52，岩石圈热结构属于“热幔冷壳”型，与汪集旸等分析的中国东部地区热结构都为为“热幔冷壳”型研究结果一致(汪集旸等，1986；唐显春等，2020)。

由上述分析可知，该区计算的地热通量102 mW/m2由三部分热量构成，对流分量为28.76±8.76 mW/m2; 传导分量中地壳放射性元素衰变产生热量为 22.5 mW/m2; 深部地幔热量为 50.74±6.18 mW/m2，三者比例约为1∶0.8∶1.76，地热田内热量约有二分之一是来自深部地幔。

7.2 地热田地热成因概念模型

综合上文研究，研究区地热成因概念模型如图7所示：在较高的区域热背景下(73.24 mW/m2)，地热田周围山区的降水入渗后通过断层或构造有利地段进入热田地下水流动系统，然后沿着热田东西两侧2 km处出露的深大断裂—招平断裂和玲珑断裂向深部运移，在运移过程中，不断接受来自地幔深部的热量(50.74±6.18 mW/m2)、岩石放射性元素衰变产生热量(22.5 mW/m2)和对流热量(28.76±8.76 mW/m2)，温度不断升高(被加热至约120～148 ℃)，在约3.6～3.8 km深的部位两条断裂相交，从而使得地热水沿着两组深大断裂交汇处形成的良好通道向上运移近地表处形成地热异常区(地热田)。两条深大断裂——招平断裂和玲珑断裂对该区地下热水的形成起到了重要作用：① 两条断裂是大气降水(冷水)径流到深部从而进行加温的主要通道；② 两条断裂在约3.8km深处交汇，从而形成深部热流上涌的有利通道。两处断裂多期活动，尤其是玲珑断裂在第四系以来仍有活动的痕迹，致使处于两条断裂中间位置的地热田内构造发育，形成以NNE向为主的构造群，对其深入调查发现，NNE向F1断裂(图7)为主要构造，破碎带宽15～20 m，发育深度大于2 km，断裂北部有砾石充填，南段张性较好。从构造行迹及附近节理互切关系分析，F1断裂为玲珑断裂的伴生断裂，至少经历了先压后张两次地壳运动，是一条长期活动至第四系且仍有活动的断裂构造。依据水文地质资料显示，地热田内F1断裂处富水性最强，大于10 m3/h，远离F1断裂，富水性逐渐减弱，进一步推断F1断裂为地热流体上涌的主要通道。

图7 招远地热田地热成因概念模型Fig. 7 The conceptual model of geothermal genesis in Zhaoyuan thermal field

8 结论

(1)通过测试分析显示，该区岩石热导率数值呈现出随着深度加深而增大的规律，数值分布在 2.5～5.3 W/(m·K)之间，平均值为4.0 W/(m·K)，高于上地壳平均热导率2.5 W/(m·K)；依据钻孔测温资料计算，该区地温梯度值为31.8 ℃/km，结合热导率测试数据，估算该钻孔地热通量为102 mW/m2，反映出该区良好的地热地质背景。

(2)在多种地温温标中，石英温标最适合用来计算该区热储温度，估算热储温度约为128.6 ℃；地下热水循环深度约3634 m。

(3)利用地震波推算出的该区深部岩石生热率分布规律与实测情况基本一致，并以此为依据推算出该区各层段的热流分布规律。通过分析，该区地热通量中既有来自传导热量(73.24 mW/m2)的贡献也有对流热量(28.76±8.76 mW/m2)的贡献，其中传导热量中地壳热流与地幔热流配分比约为1∶1.98～1∶2.52，岩石圈热结构属于“热幔冷壳”型。

(4)该区地热系统属“断控型”热储，地下热水补给通道及上涌通道主要是招平断裂和玲珑断裂两条深大断裂，并且玲珑断裂至今仍在活动；物探资料显示两条深大断裂在地下约3.8 km处交汇，这也正是本文计算的该区地下热水的热循环深度范围，验证了本文研究的正确度及可信度。