典型高温地热系统水热循环及锂同位素分馏过程模拟研究

Numerical simulation of hydrothermal cycling process and lithium isotope fractionation in a typical high-temperature geothermal system

摘要

水热系统中的多场耦合相互作用对水热循环过程和地热流体化学成分具有显著的影响。本文利用COMSOL-Multiphysics建立水—热—化学(同位素)多场耦合数值模拟模型,通过简化模型对锂同位素分馏模拟方法进行了验证。在此基础上,基于对羊八井地热田水热循环过程的认识,建立了羊八井典型剖面水热循环多场耦合模型;再现了羊八井地热系统的水热循环过程和水岩反应过程中的锂同位素分馏过程,并且讨论了主要的参数对热能聚敛效果的影响。结果表明:较高的断裂带渗透率将加快深部断裂带附近围岩温度的下降,而较低的渗透率则无法在近地表形成一定规模的水热活动;通过地表排泄量对断裂带渗透率进行约束后,发现当断裂带与深部熔融体直接沟通时,才可在近地表形成长期存在(近150 ka)的高温地热显示;在断裂带沟通了深部熔融体且熔融体热源温度不变的前提下,熔融体的具体埋深对水热活动强烈程度的影响不大;长期的水岩活动会使断裂带内锂元素大量消耗,只有当深部熔融体为断裂系统提供持续不断的物质来源时,才能保证地热流体中长期保持较高的锂元素浓度水平;基于锂同位素分馏过程估算出参与水岩反应的岩体中锂元素质量分数为25~35 mg/kg,δ⁷Li_rock为-2.0‰~0.5‰。研究成果有助于进一步认识典型高温地热系统成因机制模式。

关键词

典型高温地热系统 / 水热循环 / 水岩反应 / 多场耦合模型 / 锂同位素

Key words

typical high-temperature geothermal system / hydrothermal cycle / water-rock reaction / multifield coupling model / lithium isotope

0 引言

水热型地热系统是以水或蒸汽等形式进行热能储存、运移和转换的系统。地下水循环深度大、深部裂隙发育、水岩换热充分且揭露流体浅的构造区是优质的地热开发靶区^[1]。水热系统中发生的水热循环过程对地热资源的形成和演化具有重要意义,与传统的地下水系统不同,构成水热型地热系统的地质环境、水文环境、热量传递机制甚至热储介质都更加复杂^[2-3]。1963年Toth通过数学解析方法绘制了二维均值各向同性盆地的地下水流场,并基于此建立了区域地下水流系统理论^[4],此后该理论广泛应用在盆地水文地质研究、地热系统研究、石油勘探、核废料处置等众多方面。在地下水流系统中,重力通常作为驱动地下水流最重要的驱动力,地形或潜水面起伏决定地下水流系统的发育特征;然而在大量的研究中人们也发现地下水流系统还受到温度和化学组分等的影响;不同温度和盐度条件下水的密度和黏度会发生变化,因此温度差和盐度差可以为地下水流动提供额外的驱动力^[5-6]。在地热资源丰富的地区,地下水的温度和盐度对水流的影响不容忽视,肖巍等^[5]将物理试验和数值模拟相结合,分析了温度场和地下水流系统之间的相互影响;结果表明,较大的温差会使地下水流系统的循环深度增大,改变水流系统模式。郝奇琛等^[7]的研究表明,地下水流系统排泄区受蒸发浓缩作用的影响强烈,地下水密度的增大将削弱水流运移的驱动力。Abarca等^[8]探究了海水入侵背景下水动力弥散作用对密度驱动流的影响。

地热体系中的同位素特征通常可以记录水岩反应过程中的重要信息。锂元素作为最轻的金属元素,具有较大的同位素质量差异,其重同位素⁷Li与轻同位素⁶Li之间的相对质量差达到17%左右。水岩反应中锂元素倾向于进入流体相,同时锂同位素的分馏过程受温度影响显著,温度较低时,锂同位素分馏程度高,随着温度的升高锂同位素分馏程度逐渐降低^[9]。近些年锂同位素地球化学越来越多地被应用到了地热体系的研究中^[10]。地热流体中往往具有较高的锂元素浓度,据魏帅超等^[11]的统计,我国每年约有3 233 t锂金属随地热水排泄。田世洪等^[12]采集于拉萨地块的钾质岩样品中锂元素质量分数范围为20~35 mg/kg,δ⁷Li的范围为-3.50‰~0.5‰。于沨等^[13]基于锂同位素质量平衡拟合估算出与川西茶洛热泉水进行水岩反应的岩体的锂同位素组成。

羊八井地热田是我国典型的高温地热系统,根据热储层的不同埋深和流体性状的差异将羊八井热田的热储分为浅部热储和深部热储两部分,两者间关系密切,属同一个水热系统的两个不同部位。以往研究表明,羊八井地热流体主要来源于念青唐古拉山的冰雪融水下渗,在北区热沟附近地下深部沿断裂破碎带向上运移至浅部形成北区深部热储和浅层热储,热流体再沿多组南北向的断裂向南区运移,形成南区的浅部热储^[14-15]。本文考虑水动力场、温度场和水化学场之间的相互作用,建立羊八井典型剖面水—热—化多场耦合数值模拟模型,再现了羊八井地热系统的水热循环过程和水岩反应过程中的锂同位素分馏过程,并且讨论了主要的参数对热能聚敛效果的影响。

1 研究区概况

羊八井地热田位于地中海—喜马拉雅地热带东段(图1),属于陆—陆碰撞板缘非火山型高温地热田,其产生与欧亚板块和印度板块的陆—陆碰撞隆升过程密切相关。羊八井盆地的形成与念青唐古拉山脉隆升和花岗岩侵位之间存在动力学成因联系。自中新世开始,羊八井及邻区由区域挤压构造环境转变为区域伸展构造环境,伸展走滑构造运动渐居主导地位,逐步形成典型盆—山构造体系。受断裂带构造控制,地热流体由念青唐古拉山大量冰雪融水沿断裂带下渗并经深部热源加热后流向排泄区,形成了羊八井深部热储和浅部热储。其中浅部热储为第四系孔隙型热储,温度介于130~170 ℃;深部热储为裂隙型热储,围岩岩性为糜棱岩化花岗岩、花岗质糜棱岩和碎裂花岗岩,作为深部热储主要岩性的花岗岩可能是羊八井地热系统水岩反应过程中锂元素的主要来源,羊八井深部热储的最高温度超过了300 ℃^[15]。

目前认为最可能的热源是壳内尚未冷却的部分熔融层,岩石圈地球物理探测资料揭示了藏南地壳中部分熔融层的存在^[16]。王刚等^[17]通过岩石圈电性结构的研究指出在亚东谷露裂谷下方0~20 km深度存在一个向东倾斜的低阻通道。薛国强等^[18]的大地电磁探测结果表明,羊八井深部局部熔融体埋藏深度为6.2~14 km,这些证据表明控制羊八井水热循环的断裂带系统可能沟通了深部的局部熔融体热源(图2^[18-19])。

羊八井地热流体中含有较高浓度的锂元素,郑绵平和刘文高^[20]曾对羊八井地热田中的锂资源量进行估算,结果显示羊八井地热田锂资源量达到39万t。魏帅超等^[11]通过对全球大陆地热系统的锂同位素数据进行总结,发现地热水中锂质量浓度一般在0.5~190 000 μg/L范围内,均值为10 280 μg/L,具有岩浆型热源的地热系统往往具有更低的δ⁷Li值。Wang等^[21]的测试结果显示深部热储ZK4001井中的锂元素质量浓度高达24 000 μg/L。本研究于2022年在羊八井采集了9组地热流体进行锂元素及锂同位素测试,结果显示羊八井浅部热储地热流体中锂元素平均质量浓度为7 009 μg/L,锂同位素δ⁷Li的平均值为1.23‰。

2 地热系统水—热—化多场耦合模型

2.1 模拟工具及控制方程

利用COMSOL建立水—热—化多场耦合模型,该程序曾被广泛应用于地热等场景的多物理场耦合模拟研究中且被证明是有效的^[22⇓-24]。地热系统的水热循环过程涉及地下水流动过程和热量的传递过程,以及地热流体中化学物质的传递过程。将地下水流动过程概化为等效多孔介质中的流动,符合达西定律^[25]:

(1) ∇(ρu)=Q_m

式中:ρ为流体密度,kg/m³;u为达西流速,m/s;Q_m为流体质量源,kg/(m³·s)。

达西流速可根据下式计算:

(2)u=- \frac{k}{\mu }(∇p+ρg)

式中:μ为流体动力黏度,Pa·s;k为储层渗透率,m²;p为流体压力,Pa;g为重力加速度,m/s²。

多孔介质中的热量传递过程通过能量守恒方程描述^[26]:

(3)(ρC_p)_eff \frac{\partial T}{\partial t}+ρ_fC_P,fu∇T=∇(K_eff∇T)+Q

式中:(ρC_p)_eff为储层等效体积热容,kJ/(m³·K);T为温度,K;C_p,f为流体恒压热容,J/(kg·K);K_eff为热储层有效导热系数,W/(m·K);Q为热源,W/m³。

多孔介质中稀物质传递场中溶质的运动可以通过对流弥散方程来描述^[27]:

(4) \frac{\partial c_i}{\partial t}=D∇c_i-u∇c_i

式中:c_i为流体中的溶质浓度,mol/m³; D为流体扩散系数,m²/s。

模型中考虑水动力场、温度场和化学传递场之间的相互耦合,根据温度和盐度对流体密度的影响^[28],Kuder^[29]曾总结出流体密度随流体温度、压力和盐度变化的公式:

(5) ρ=ρ_T+Δp_p+ρ_TDS

式中:ρ_T表示特定温度下水的密度,kg/m³;Δp_p为压力校正项,相较于温度对流体密度的影响来说,压力对流体密度的影响很小,本文暂不考虑压力校正项;ρ_TDS表示流体的盐度,g/L。

温度对流体密度的影响则采用NIST官网数据进行拟合得出的经验公式:

(6) ρ_T=1 023.8-1.083 5T+0.008 4T²-
4×10^-5T³+3×10^-8T⁴

Batzle在相关研究中定义了水的黏度随温度和盐度的变化关系^[30]:

(7) μ=(0.1+0.333×ρ_TDS/10⁶)+
[1.65+91.9×(ρ_TDS/10⁶)³]×exp(-α)

(8) α={0.42×[(ρ_TDS/10⁶)^0.8-0.17]²+0.045}×T^0.8

式中μ表示特定温度和盐度下水的流体黏度,单位为mPa·s。

以往研究表明,在水岩反应过程中,锂元素优先进入流体相,不同的同位素进入流体的速率不同导致锂同位素分馏。模型中考虑了水岩反应中锂同位素的分馏过程,将锂的重同位素(⁷Li)与轻同位素(⁶Li)视为两个体系,假设水岩反应为一级反应,化学反应速率与物质浓度呈正比^[31],则有:

((9)) R_r^i=k_fC_i

式中: R_r^i表示物质i参与水岩反应的速率,mol/(m³·s);k_f表示反应速率常数,s^-1;C_i表示物质i的浓度,mol/m³。反应速率常数k_f可由Arrhenius公式求得:

(10)k_f=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}

式中:A_f表示Arrhenius指前因子,s^-1;E_a表示化学反应活化能,J/mol;R表示摩尔气体常数,J/(mol·K);T表示温度,K。

根据同位素分馏系数的定义,锂同位素分馏系数可以表示为:α=(⁷L {\mathrm{i}}_{\mathrm{a}\mathrm{l}\mathrm{t}}^t/⁶L {\mathrm{i}}_{\mathrm{a}\mathrm{l}\mathrm{t}}^t)/(⁷L {\mathrm{i}}_{\mathrm{w}}^t/⁶L {\mathrm{i}}_{\mathrm{w}}^t),若假设参与水岩反应的岩石中锂元素浓度保持不变,且水岩反应达到平衡的条件下,存在平衡常数 K_{\mathrm{p}}^6 {=}^6L {\mathrm{i}}_{\mathrm{a}\mathrm{l}\mathrm{t}}^t/⁶L {\mathrm{i}}_{\mathrm{w}}^t; K_{\mathrm{p}}^7=⁷L {\mathrm{i}}_{\mathrm{a}\mathrm{l}\mathrm{t}}^t/⁷L {{\mathrm{i}}_{\mathrm{w}}^{\mathrm{t}}}^{\mathrm{}}^[13],结合同位素分馏系数的定义有

(11) K_p^7=α· K_p^6

地热流体中同位素的分馏过程表现为水岩反应中不同核素进入流体的速率不同,对于⁶Li体系来说,分别有溶解速率 R_{\mathrm{d}}^6和沉淀速率 R_{\mathrm{p}}^6:

(12) R_{\mathrm{d}}^6=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}·L {\mathrm{i}}_{\mathrm{r}}^6/ K_{\mathrm{p}}^6

(13) R_p^6=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}·L {\mathrm{i}}_{\mathrm{w}}^6

式中:L {\mathrm{i}}_{\mathrm{r}}^6表示岩石中⁶Li的浓度;L {\mathrm{i}}_{\mathrm{w}}^6表示流体中⁶Li的浓度。

结合式(12)和(13),有净反应速率:

(14) R_w^6=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}(L {\mathrm{i}}_{\mathrm{r}}^6/ K_p^6-L {\mathrm{i}}_{\mathrm{w}}^6)

同理,对于⁷Li体系来说,存在净反应速率:

(15) R_w^7=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}(L {\mathrm{i}}_{\mathrm{r}}^7/ K_p^7-L {\mathrm{i}}_{\mathrm{w}}^7)

与式(11)结合后得

(16) R_w^7=A_f {\mathrm{e}}^{-\frac{E_a}{RT}}[L {\mathrm{i}}_{\mathrm{r}}^7/(α· K_p^6)-L {\mathrm{i}}_{\mathrm{w}}^7]

根据质量守恒定律,在水岩反应中,流体中获得的锂元素质量就等于岩石中损失的锂元素质量,于是有质量平衡方程:

(17)ρ_wθ \frac{d{C}_w}{\mathrm{d}t}=ρ_r(1-θ) \frac{d{C}_r}{\mathrm{d}t}

式中:ρ_w表示水的密度;θ表示孔隙率; \frac{\mathrm{d}{C}_w}{\mathrm{d}t}表示水中锂元素浓度的变化率,可以用反应速率R_w来代替;ρ_r表示岩石的密度; \frac{\mathrm{d}{C}_r}{\mathrm{d}t}表示岩石中锂元素浓度的变化率,可以表示为岩石中的反应速率R_r。

根据上述质量平衡方程,可以将岩石中锂元素浓度的变化率(反应速率)表示为

(18)R_r= \frac{{\rho }_w\theta }{{\rho }_r(1-\theta )}R_w

以往的试验研究证明水岩作用过程中锂同位素分馏值(ΔLi_{solution-solid})对温度具强烈依赖性。已有的锂同位素分馏试验多是利用玄武岩进行的,玄武岩中主要的含锂矿物为锂辉石,而研究区酸性岩储层中主要的含锂矿物为云母。相关研究表明,锂云母和锂辉石都具有六面体配位,同位素分馏系数随温度的变化趋势较为相似。高温条件下的玄武岩水岩反应试验结果显示:当T=2 ℃时,Δ=+19‰;当T=160 ℃时,Δ=+9‰;当T=350 ℃时,Δ=+4‰^[9,32⇓-34]。Wunder等^[35-36]通过试验得到了高温(500~900 ℃)条件下锂同位素分馏系数,结果表明,即使在非常高的温度(900 ℃)条件下,依然存在约1‰的同位素分馏。

2.2 模型验证

COMSOL中多孔介质的渗流和传热过程模型已经在大量的研究中被验证^[22,25-26],本节建立了一个简化的水—热—化(同位素分馏)多场耦合模型,验证本文提出的锂同位素分馏模型。模型几何面为长方形区域(图3),长度为10 km,宽度为1 km。左端为流体入口,右端为流体出口,温度设置为固定值(100、200、300、400和500 ℃),以验证不同温度条件下同位素分馏系数计算结果;入口锂元素组成设为初始流体的锂同位素组成(锂元素质量浓度185 μg/L,δ⁷Li为32‰);岩石中的锂元素质量分数设置为15 mg/kg,锂同位素组成δ⁷Li设置为0‰。水岩反应方程中的E_a设置为80 kJ/mol,R设置为8.314 56 kJ/(mol·K)^[37];调整A_f来控制水岩反应的速率,并使同位素分馏体系尽快达到平衡状态。

图4显示,水岩反应过程中岩石中的锂元素转移到反应流体中被消耗带走。A_f取值越大,水岩反应速率越快,流体中锂元素浓度的过渡区域也越小,这意味着水岩反应在较短时间内完成。当A_f增大到0.01 s^-1时,锂浓度的过渡区域不足2 km,这意味着水岩反应在不到2 km的流程内完成。

水岩反应中锂元素在岩石和流体中的分配比例由平衡系数 K_p^6来控制,分别计算了不同温度条件下出口处岩石和流体中的锂元素浓度。结果(图5)显示,不同温度条件下岩石与流体中锂元素含量之比与给定的分配系数K_p相符合。

图6显示:在水岩反应锋面未到达出口附近时,岩石中锂同位素δ⁷Li值未发生变化,流体中锂同位素δ⁷Li值较好地满足同位素分馏系数α的要求;当水岩反应锋面到达出口后,流体变为初始流体的锂同位素组成。

通过以上对简化模型中水岩反应过程的模拟,可以认为本文对水岩反应中锂同位素分馏过程的概化达到了预期的效果:水岩反应可以在较短时间内达到平衡,且水岩平衡中锂元素在岩石与流体中的分布满足平衡系数K_p的要求,水岩反应过程中的锂同位素分馏程度满足同位素分馏系数α的要求,因此可将这种方法应用于地热流体循环模拟研究中。

3 羊八井地热系统典型剖面模拟

3.1 典型剖面模型概化

为了定量刻画羊八井地热田的水热循环模式,基于念青唐古拉山伸展构造图(图2)和羊八井邻区地表DEM高程数据,建立了羊八井地热田典型剖面水热循环多场耦合数值模拟模型(图7)。模型宽度约40 km,顶部边界为地表高程,模型高度10~12 km。在模型中概化出深部的断裂带系统和深部的熔融体热源。采用自由三角形网格剖分方式将模拟区域共剖分为50 936个单元。

在渗流场中,将底部和侧边界设置为无流动边界,顶部边界设置为给定压力边界(大气压力);在传热场中,将底部边界设置为给定热流边界,侧边界设置为热绝缘边界,顶部边界温度设置为当地年平均气温;在化学物质传递场中,将补给区大气降水/冰雪融水中锂同位素组成设置为水岩反应初始流体的值,并假定水岩反应只在断裂带范围内进行。

地热活动区的泉华记录了水热活动的时间和期次^[38]。李振清^[39]根据西藏南部地区泉华年龄特征构建了水热活动时空格架,将西藏南部的水热活动分为4个时期,其中最后一个时期始于15万年以前,同时该期水热活动相较于前两期更强烈。此外,李家振等^[40]通过铀系法测定的羊应乡水热活动同位素年代为14.6万年左右。因此参照前人数据,将羊八井地热系统典型剖面模型模拟时间跨度设定为15万年。

3.2 模型参数

研究水热系统的基础是背景地温特征,大地热流是分析地热系统深部温度场非常重要的参数^{[41⇓⇓-44]}。以往研究表明,大地热流的主要贡献包括来自地幔的深部热流和来自地壳岩石的放射性生热^[45⇓-47]。我国青藏高原巨厚的增厚地壳使得该地区壳内热流值份额占比较高,此外青藏高原地壳广泛分布的“壳内低速层”也是非常重要的恒常热源^{[48⇓⇓⇓-52]}。羊八井盆地受热对流效应影响显著,因此热流数值变动幅度较大(表1);沈显杰等^[53]1989年在羊八井进行的钻孔热流测量结果显示,羊八井深部提供了量级为83~108 mW/m²的传导型热流。基于本模型的假设进行分析,这部分热流的来源可能包括下部边界向上传导的深部热流、熔融体热源提供的热流和岩石放射性元素衰变产生的热量。

以往研究表明上地壳岩石放射性生热量差别不大,平均为1.58~1.67 μW/m^{3 [54-55]},模型中取值1.6 μW/m³。羊八井地区上地壳主要发育的岩石岩性为花岗岩,热导率取值范围一般在2.2~3.2 W/(m·K),模型取中间值2.5 W/(m·K),大量研究成果表明,羊八井深部存在局部熔融体热源^[56-57]。花岗岩的深熔温度为600~700 ℃^[15]。根据薛国强等^[18]2023年最新的大地电磁探测成果,羊八井深部熔融体的深度范围在6.2~14 km,平均深度为10 km。Su和Tan^[58]对亚东谷露裂谷的数值研究也表明,深度超过10 km的岩浆对地热异常的影响有限,因此本文讨论熔融体埋深6~10 km的情况。

对于对流型地热系统来说,渗透率是决定流体循环快慢并影响地热发育特征十分关键的参数。以往研究中对羊八井浅部沉积型热储层渗透率的估计值在(1~10)×10^-15 m²范围内;断裂带的渗透率根据经验值在模型中被设定为10^-15~10^-13 m^{2 [58]};深部热储层的渗透率被视为与断裂带渗透率相同;对于结晶基底来说,地层渗透率一般随深度的增加而降低。本模型中基底渗透率根据Stober等^[59]利用深井实验数据拟合出的渗透率与深度之间的关系给定:

(19) lg k=-1.38×lg z-15.4

式中:k表示渗透率, m²;z表示深度, km。

模型采用的其他参数依照经验值给定,主要参数见表2。

4 结果与讨论

4.1 背景地温场模拟

基于前文对模型和参数进行的概化,通过稳态模拟刻画了不考虑对流的背景地温场(图8)。考虑局部熔融体和岩石放射性生热对热流的贡献且模型底部边界采用西藏地热地质大队在拉萨北郊实测的大地热流值(66 mW/m²)时,利用稳态模型计算出的羊八井地区传导型热流范围为82~130 mW/m²,平均值约为109 mW/m²,与沈显杰等^[53]的实测数据接近。将此地温场作为初始条件,模拟讨论断裂对流作用下典型高温地热田的热能聚敛机制。

4.2 断裂带渗透率对地热系统的影响

首先讨论断裂渗透率对水热活动的影响。设定的断裂沟通深度约为6 km,熔融体深度约为10 km。

图9显示了不同渗透率条件下断裂带热对流的温度场分布及其随时间的变化。随着念青唐古拉山的大气降水和冰雪融水沿着断裂带下渗,断裂带附近温度降低,下渗流体在断裂带深部被加热并经排泄区断裂向上运移至羊八井深部热储层,并与来自山区的冷地下水进一步混合后向浅部热储层排泄。随着水热循环,断裂带深部及其周围岩体的温度出现明显下降,导致深循环地热流体所能带到地表的热量减少,表现为高温异常的强度减弱和范围缩小。

图10显示,断裂带渗透率越大,近地表高温异常形成越快,同时地热显示可持续的时间也越短。当断裂带渗透率为1×10^-13 m²时,150 ka后浅部热储平均温度降低到约40 ℃,深部热储中的最高温度降低到约76 ℃。

较大的断裂带渗透率加快了深部水热循环,无法维持近地表长期的高温异常;而较小的渗透率会使水热循环过慢,无法在近地表形成一定规模的水热活动。羊八井曾是我国地表水热活动最强烈的地区之一,根据1977年的实测数据,羊八井热水湖与热泉等地表排泄量约为107.5 kg/s^[60]。

根据羊八井地热田水热活动范围,垂直于模型剖面方向上的延伸长度为2~3 km,因此剖面模型的排泄量应在0.036~0.054 kg/(m·s)之间。结合改变断裂带渗透率的模拟结果(图11),可以得出在基于本文假设的断裂发育情况条件下,断裂带渗透率在(1~5)×10^-14 m²范围内时,模拟出的地表排泄量与羊八井实际的水热活动规模相符。

4.3 断裂带沟通深度对地热系统的影响

采用固定的断裂带渗透率(2×10^-14 m²)和固定的熔融体深度(10 km),改变断裂带沟通深度(6、8和10 km)模拟水热循环过程(图12)。随着断裂带沟通深度的增加,流体接触的岩体温度越来越高,热储温度也越高。当断裂带未与深部熔融体热源直接接触时,其间的热量传递主要通过热传导的方式进行,长时间的水热循环后传导热流不足以支撑断裂对流形成的地热异常,表现为地表高温异常持续时间较短;而当断裂带与深部熔融体直接接触时,便可维持更长时间的高温异常。

图13显示,断裂带沟通深度越大,热储温度越高。在断裂带沟通深度为6和8 km的情况下,模型运行150 ka后深部热储中的最高温度都不超过250 ℃;而断裂带深度在10 km的情况下,由于与深部熔融体热源直接接触,150 ka后深部热储最高温度约为357 ℃。羊八井ZK4002测井成果显示,该孔在深1 850 m处的井温高达300 ℃以上^[15]。由此可推断羊八井深部断裂带应该是直接沟通了深部熔融体热源或距离熔融体热源非常近,才可以在近地表形成如此高温且持续的地热异常。

4.4 熔融体热源深度对地热系统的影响

保持断裂带与深部熔融体之间的接触关系,改变熔融体热源的深度(6、8和10 km)进行水热循环的模拟。模拟结果(图14)显示,不同熔融体埋深情况下的水热循环过程表现出相似的温度场变化。

图15显示,3种熔融体深度方案下模型运行150 ka后浅部热储中的平均温度在140 ℃左右,同时表现出非常缓慢的温度下降趋势,熔融体埋深从6 km增加到10 km,150 ka后浅部热储中平均温度从140.90 ℃变化为139.36 ℃,平均温度之差约1.54 ℃;深部热储中最高温度的变化趋势也十分相似,熔融体埋深从6 km增加到10 km,模型运行150 ka后深部热储中的最高温度由350.97 ℃变为337.87 ℃,最高温度之差约13 ℃,相较于高于300 ℃的温度异常来说差别很小。与传导型热流差异情况相比,这部分差异对温度异常的影响十分有限。也就是说,在断裂带与深部熔融体直接沟通的情况下,近地表可形成长期存在的高温地热显示,熔融体温度决定地表热异常的强烈程度,而具体熔融体埋深对地表水热活动强烈程度的影响可能不大。

4.5 水岩反应中的锂同位素分馏模拟

为了讨论锂同位素组成与深部水热循环和水岩反应之间的相互影响,参照前文计算结果将熔融体埋深固定为10 km,断裂带渗透率固定为2×10^-14 m²。根据田世洪等^[12]对拉萨地块大陆岩石圈锂同位素的研究,采集于羊应地区的钾质岩样品中锂元素质量分数范围为20~35 mg/kg,δ⁷Li范围为-3.50‰~0.5‰,在模型中先取中间值锂元素质量分数27.5 mg/kg,δ⁷Li为-1.5‰。参照于沨等^[13]的研究,将稀释海水作为水岩反应初始流体(锂元素质量浓度185 μg/L,δ⁷Li=32‰)来模拟水热循环过程中的锂同位素分馏过程。

断裂带的孔隙率决定了水岩反应中参与的流体和岩体的比例。不考虑深部熔融体为断裂带提供物质来源的模拟结果(图16a)显示,随着水热循环及水岩反应的进行,排泄区地热流体中锂元素浓度经历了短暂的上升(受温度变化影响)后出现快速下降,这是由于岩石中的锂元素被水岩反应消耗所致。孔隙率越小,断裂带中岩石所占比例越大,相应的水岩反应过程中锂元素的消耗越慢。当孔隙率为10% 时,60 ka后排泄区锂元素质量浓度下降至不足1 000 μg/L。然而大规模断裂带的孔隙率通常大于10%,具有较好的导水性能^[61],这意味着在深部熔融体未向断裂带提供锂元素的条件下,长期的水岩活动会大量消耗断裂带内的锂元素,无法满足现今地表观测到的锂元素浓度情况。

考虑熔融体为断裂带提供物质来源(假设断裂带与熔融体接触边界处岩石中的锂元素浓度和锂同位素组成不变)的计算结果(图16b)显示,由于深部熔融体为水岩反应提供持续的物质来源,经历长期的水热循环/水岩反应后排泄区流体中的锂元素浓度依然处在较高的水平。此外,不同孔隙率条件下150 ka后排泄区流体中的锂元素浓度和锂同位素组成非常接近。

调整参与水岩反应岩石锂元素质量分数和锂同位素组成并与实测数据进行拟合。模拟结果(图17)显示,岩石中的锂元素质量分数及锂同位素组成对排泄区地热流体的锂同位素组成具有明显的控制作用。当岩石锂元素质量分数在25~35 mg/kg范围时,模型稳定后地表排泄区地热流体的平均锂元素质量浓度保持在6 120~8 517 μg/L,与实测数数值相近(实测浅部热储锂元素质量浓度平均值约为7 009 μg/L);当参与水岩反应岩石中锂同位素组成δ⁷Li_rock在-2.0‰~0.5‰范围时,排泄区地热流体的平均锂同位素组成δ⁷Li在0.50‰~3.50‰之间,与实测数值相近(实测锂同位素组成δ⁷Li平均值约为1.23‰)。据此推测锂元素的物质来源应该是青藏高原富含锂的新生代长英质岩浆岩,田世洪等^[12]采集于羊应乡的钾质岩样品中锂元素质量分数范围在20~35 mg/kg,δ⁷Li_rock的范围在-3.50‰~0.5‰,这一组成类似于印度洋下地壳的Li同位素组成,指示热量来源可能与俯冲过程有关。

羊八井水热系统的地热流体经历了念青唐古拉山大量冰雪融水沿断裂带下渗并经深部热源加热后流向排泄区的过程,并形成羊八井深部热储和浅部热储^[15]。泉华和水热活动期次的研究表明,包括羊八井地热田在内的青藏高原中部水热活动至少持续了几十万年,其中最后一个期次的水热活动大约从150 ka前持续至今^[39]。沟通深部的断裂系统是控制水热系统发育规模的关键因素,造山带中的水热系统通常表现为通过深大断裂的水热循环过程^[62]。驱动水流运动的水力梯度来自大气降水补给区高程、地势差和导水断层的共同作用;深大断裂的存在将地下水补给区与深部的热源连接起来,导致地下水可能循环到的深度增大。以往研究表明,断裂带型地热系统的循环深度可能超过10 km^[63]。根据薛国强等^[18]最新的大地电磁探测结果,羊八井地热田深部熔融体埋深为6.2~14 km,因此羊八井断裂带直接沟通深部熔融体热源是很有可能的。本文中多场耦合数值模型的计算结果也从热量来源和锂元素的物质来源方面为断裂带直接沟通深部热源提供了证据。模拟结果表明,当断裂带未达到熔融体热源深度时,地表难以维持如此长时间的高温地热活动,同时地热流体中也无法长期保持较高的锂元素浓度水平。

5 结论

本文通过建立水—热—化多场耦合模型再现了羊八井深部水热循环过程,分析了不同影响因素对水热循环过程的影响,认识到不同假设条件下热储温度的变化和水岩反应中锂同位素分馏过程的模拟结果共同约束了羊八井深部水热循环过程。具体结论如下。

(1)沟通深部的断裂系统是控制水热系统发育规模的关键因素,过小的断裂带渗透率无法满足近地表水热活动排泄量的约束;断裂带沟通深度不够则无法提供长期高温水热活动所需的热量来源;当断裂带沟通到深部的高温熔融体热源且断裂带具有足够的渗透性能时,可维持近地表长期高温水热活动。

(2)对水岩反应中锂同位素分馏过程的模拟,进一步表明深部熔融体热源与断裂系统之间存在接触关系。水岩反应中锂元素倾向于进入流体相,只有深部熔融体为断裂系统提供持续不断的物质来源才能保证地热流体中长期保持较高的锂元素浓度水平;进一步地,基于锂同位素分馏过程估算出参与水岩反应的岩体中锂元素质量分数为25~35 mg/kg,δ⁷Li_rock为-2.0‰~0.5‰。

(3)本文对于锂同位素分馏过程的模拟主要考虑到温度对同位素分馏系数的影响。实际上,在水岩相互作用过程中锂元素向流体迁移的水化学行为还受到吸附等作用的共同影响,未来针对这一复杂的水化学行为还需要更加深入的研究。

参考文献

[1] 多吉, 王贵玲. 加大深部热能探采技术攻关持续推进地热资源规模化开发[J]. 科技导报, 2022, 40(20): 1.

[2] 王贵玲, 蔺文静. 我国主要水热型地热系统形成机制与成因模式[J]. 地质学报, 2020, 94(7): 1923-1937.

[3] WANG G L, GAN H N, LIN W J, et al. Hydrothermal systems characterized by crustal thermally-dominated structures of southeastern China[J]. Acta Geologica Sinica (English Edition), 2023, 97(4): 1003-1013.

[4] TÓTH J. A theoretical analysis of groundwater flow in small drainage basins[J]. Journal of Geophysical Research, 1963, 68(16): 4795-4812.

[5] 肖巍, 孙蓉琳, 陈明霞, 不同地下水流系统模式渗流场和温度场的互相影响[J]. 地质科技通报, 2022(1): 251-259.

[6] BAHLALI M L, SALINAS P, JACKSON M D. Efficient numerical simulation of density-driven flows: application to the 2- and 3-D elder problem[J]. Water Resources Research, 2022, 58(8): e2022W-e32307W.

[7] 郝奇琛, 崔伟哲, 黄林显. 盆地地下水密度变化对水流驱动力的影响[J]. 济南大学学报(自然科学版), 2020, 34(6): 595-602.

[8] ABARCA E, CARRERA J, SÁNCHEZ-VILA X, et al. Anisotropic dispersive Henry problem[J]. Advances in Water Resources, 2007, 30(4): 913-926.

[9] MILLOT R, SCAILLET B, SANJUAN B. Lithium isotopes in island arc geothermal systems: Guadeloupe, Martinique (French West Indies) and experimental approach[J]. Geochimica et Cosmochimica Acta, 2010, 74(6): 1852-1871.

[10] 余小灿, 刘成林, 王春连. 锂同位素地球化学在大陆地热体系研究中的应用[J]. 地球科学进展, 2020, 35(3): 246-258.

[11] 魏帅超, 张薇, 付勇, 我国地热水中锂元素分布特征及资源开发利用[J]. 中国地质, 2024, 51(5): 1527-1553.

[12] 田世洪, 路娜, 侯增谦, MC-ICP-MS锂同位素溶液分析技术与应用[J]. 地质论评, 2021, 67(5): 1441-1464.

[13] 于沨, 于扬, 王登红, 锂同位素地球化学在地热流体水岩反应中的应用: 以川西现代富锂热泉研究为例[J]. 岩石学报, 2022, 38(2): 472-482.

[14] 陈卫营, 薛国强, 赵平, 西藏羊八井地热田SOTEM探测及热储结构分析[J]. 地球物理学报, 2023, 66(11): 4805-4816.

[15] 多吉. 典型高温地热系统: 羊八井热田基本特征[J]. 中国工程科学, 2003, 5(1): 42-47.

[16] UNSWORTH M J, TEAM T I M, JONES A G, et al. Crustal rheology of the Himalaya and Southern Tibet inferred from magnetotelluric data[J]. Nature, 2005, 438(7064): 78-81.

[17] 王刚, 魏文博, 金胜, 冈底斯成矿带东段的电性结构特征研究[J]. 地球物理学报, 2017, 60(8): 2993-3003.

[18] 薛国强, 陈卫营, 赵平, 西藏羊八井地热田三维电性结构模型: 来自大地电磁的证据[J]. 中国科学: 地球科学, 2023, 53(8): 1859-1871.

[19] 吴珍汉. 青藏高原腹地的地壳变形与构造地貌形成演化过程[M]. 北京: 地质出版社, 2003.

[20] 郑绵平, 刘文高. 西藏发现富锂镁硼酸盐矿床[J]. 地质论评, 1982, 28(3): 263-266.

[21] WANG X, WANG G L, LU C, et al. Evolution of deep parent fluids of geothermal fields in the Nimu-Nagchu geothermal belt, Tibet, China[J]. Geothermics, 2018, 71: 118-131.

[22] XU H R, LIU G H, ZHAO Z H, et al. Coupled THMC modeling on chemical stimulation in fractured geothermal reservoirs[J]. Geothermics, 2023, 116: 102854.

[23] SHI H L, WANG G L, LU C. Numerical investigation on delaying thermal breakthrough by regulating reinjection fluid path in multi-aquifer geothermal system[J]. Applied Thermal Engineering, 2023, 221: 119692.

[24] SAEED M, MRITYUNJAY S, AYSEGUL T, et al. Simulations and global sensitivity analysis of the thermo-hydraulic-mechanical processes in a fractured geothermal reservoir[J]. Energy, 2022, 247: 123511.

[25] WANG G L, LIU G H, ZHAO Z H, et al. A robust numerical method for modeling multiple wells in city-scale geothermal field based on simplified one-dimensional well model[J]. Renewable Energy, 2019, 139: 873-894.

[26] MA F, LIU G H, ZHAO Z H, et al. Coupled thermo-hydro-mechanical modeling on the Rongcheng geothermal field, China[J]. Rock Mechanics and Rock Engineering, 2022, 55(8): 5209-5233.

[27] ZHANG D X, KANG Q J. Pore scale simulation of solute transport in fractured porousmedia[J]. Geophysical Research Letters, 2004, 31(12): 289-302.

[28] RANJRAM M, GLEESON T, LUIJENDIJK E. Is the permeability of crystalline rock in the shallow crust related to depth, lithology, or tectonic setting?[J]. Geofluids, 2005, 15: 106-119.

[29] KUDER J. Methoden zur berechnung von fluidparametern[R]. Hannover: Methoden zur Berechnung von Fluidparametern, 2011.

[30] BATZLE M, WANG Z J. Seismic properties of pore fluids[J]. Geophysics, 1992, 57(11): 1396-1408.

[31] 徐浩然, 程镜如, 赵志宏. 华北地区碳酸盐岩热储层酸化压裂模拟方法与应用[J]. 地质学报, 2020, 94(7): 2157-2165.

[32] LUI-HEUNG C, GIESKES J M, YOU C F, et al. Lithium isotope geochemistry of sediments and hydrothermal fluids of the Guaymas Basin, Gulf of California[J]. Geochimica et Cosmochimica Acta, 1994, 58(20): 4443-4454.

[33] JAMES R H, RUDNICKI M D, PALMER M R. The alkali element and boron geochemistry of the Escanaba Trough sediment-hosted hydrothermal system[J]. Earth and Planetary Science Letters, 1999, 171(1): 157-169.

[34] LUI-HEUNG C, EDMOND J M. Variation of lithium isotope composition in the marine environment: a preliminary report[J]. Geochimica et Cosmochimica Acta, 1988, 52(6): 1711-1717.

[35] WUNDER B, MEIXNER A, ROMER R L, et al. Temperature-dependent isotopic fractionation of lithium between clinopyroxene and high-pressure hydrous fluids[J]. Contributions to Mineralogy and Petrology, 2006, 151(1): 112-120.

[36] WUNDER B, MEIXNER A, ROMER R L, et al. Li-isotope fractionation between silicates and fluids: pressure dependence and influence of the bonding environment[J]. European Journal of Mineralogy, 2011, 23(3): 333-342.

[37] JENKIN G R T, LINKLATER C, FALLICK A E. Modeling of mineral δ¹⁸O values in an igneous aureole: closed-system model predicts apparent open-system δ¹⁸O values[J]. Geology, 1991, 19(12): 1185-1188.

[38] 王潇. 西藏羊易地热田泉华地球化学特征及其指示意义[D]. 北京: 中国地质科学院, 2018.

[39] 李振清. 青藏高原碰撞造山过程中的现代热水活动[D]. 北京: 中国地质科学院, 2002.

[40] 李家振, 孙善平, 张有瑜, 西藏羊应乡地热田形成特点及评价探讨[J]. 现代地质, 1994, 8(1): 49-56.

[41] 宋嘉佳, 王贵玲, 邢林啸, 岩石热导率校正对大地热流计算值的影响: 以渤海湾盆地冀中坳陷为例[J]. 地质论评, 2023, 69(4): 1349-1364.

[42] FURLONG K P, CHAPMAN D S. Heat flow, heat generation, and the thermal state of the lithosphere[J]. Annual Review of Earth and Planetary Sciences, 2013, 41: 385-410.

[43] 王贵玲, 蔺文静, 刘峰, 地热系统深部热能聚敛理论及勘查实践[J]. 地质学报, 2023, 97(3): 639-660.

[44] 王贵玲, 马峰, 侯贺晟, 松辽盆地坳陷层控地热系统研究[J]. 地球学报, 2023, 44(1): 21-32.

[45] LIN W J, WANG G L, GAN H N, et al. Heat source model for Enhanced Geothermal Systems (EGS) under different geological conditions in China[J]. Gondwana Research, 2023, 122: 243-259.

[46] 蔺文静, 王贵玲, 甘浩男. 华南陆缘火成岩区差异性地壳热结构及地热意义[J]. 地质学报, 2024, 98(2): 544-557.

[47] 王贵玲, 刘峰, 蔺文静, 我国陆区地壳生热率分布与壳幔热流特征研究[J]. 地球物理学报, 2023, 66(12): 5041-5056.

[48] 王强, 苟国宁, 张修政, 青藏高原中北部地壳流动与高原扩展: 来自火山岩的证据[J]. 中国科学基金, 2017, 31(2): 121-127.

[49] WANG Q, CHUNG S L, LI X H, et al. Crustal melting and flow beneath northern Tibet: evidence from mid-Miocene to quaternary strongly peraluminous rhyolites in the southern Kunlun range[J]. Journal of Petrology, 2012, 53(12): 2523-2566.

[50] 龙登红, 周小龙, 杨坤光, 青藏高原东北缘深部地质构造与地热资源分布关系研究[J]. 中国地质, 2021, 48(3): 721-731.

[51] 邱楠生, 胡圣标, 何丽娟. 沉积盆地地热学[M]. 青岛: 中国石油大学出版社, 2019.

[52] 沈显杰, 张文仁, 杨淑贞, 青藏高原南北地体壳幔热结构差异的大地热流证据[J]. 中国地质科学院院报, 1990, 11(2): 203-214.

[53] 沈显杰, 张文仁, 杨淑贞, 西藏中部地热区的钻孔热流测量[J]. 地质科学, 1989, 24(4): 376-384.

[54] JAUPART C, MARESCHAL J C, IAROTSKY L. Radiogenic heat production in the continental crust[J]. Lithos, 2016, 262: 398-427.

[55] 何丽娟, 胡圣标, 汪集旸. 中国东部大陆地区岩石圈热结构特征[J]. 自然科学进展, 2001, 11(9): 72-75.

[56] WANG G, WEI W B, YE G F, et al. 3-D electrical structure across theYadong-Gulu Rift revealed by magnetotelluric data: new insights on the extension of the upper crust and the geometry of the underthrusting Indian lithospheric slab in southern Tibet[J]. Earth and Planetary Science Letters, 2017, 474: 172-179.

[57] CHEN L S, BOOKER J R, JONES A G, et al. Electrically conductive crust in southern Tibet from INDEPTH magnetotelluric surveying[J]. Science, 1996, 274(5293): 1694-1696.

[58] SU J B, TAN H B. The genesis of rare-alkali metal enrichment in the geothermal anomalies controlled by faults and magma along the northernYadong-Gulu Rift[J]. Ore Geology Reviews, 2022, 147: 104987.

[59] STOBER I, BUCHER K. Hydraulic conductivity of fractured upper crust: insights from hydraulic tests in boreholes and fluid-rock interaction in crystalline basement rocks[J]. Geofluids, 2015, 15(1/2): 161-178.

[60] 姚足金, 张钖根, 安可士, 西藏羊八井地热资源评价[R]. 石家庄: 水文地质工程地质研究所, 1984.

[61] 许增光, 曹成, 柴军瑞, 断层带破碎岩体非达西渗流特性及模型研究[J]. 岩石力学与工程学报, 2023, 42(增刊2): 4099-4108.

[62] MENZIES C D, TEAGLE D A H, CRAW D, et al. Incursion of meteoric waters into the ductile regime in an activeorogen[J]. Earth and Planetary Science Letters, 2014, 399: 1-13.

[63] DIAMOND L W, WANNER C, WABER H N. Penetration depth of meteoric water inorogenic geothermal systems[J]. Geology, 2018, 46(12): 1063-1066.

基金资助

国家重点研发计划项目(2021YFB1507300)

国家重点研发计划项目(2021YFB1507302)

中国地质调查局地质调查项目(DD20230019)

中国地质调查局地质调查项目(DD20221676)