高温地热生产井碳酸钙结垢定量评价：水文地球化学

Quantitative Assessment of Calcite Scaling of a High Temperature Geothermal Production Well: Hydrogeochemistry

摘要

高温地热井井筒结垢是地热开发遇到的突出问题之一，其中涉及到复杂的水文地球化学过程.本文首先联合井筒两相流动、水‒气‒结垢物水文地球化学和井壁粘附模型，建立了井筒碳酸钙结垢定量评价的耦合模型.针对西藏羊八井地热田典型井开展了井口取样和分析，结果显示井筒结垢物为碳酸钙，地热流体中碳酸盐矿物过饱和，非凝结气体主要成分是CO₂.利用建立的耦合模型及实测流体数据，定量评价了井筒碳酸钙的结垢位置和速率，并给出了其在井筒中的结垢形状.CO₂分压对碳酸钙的析出有控制性的作用，最大结垢厚度位置发生在闪蒸面以上10~20 m位置，持续1年的开采最大结垢厚度在考虑井筒粘附动力学时为14~25 mm，假设完全粘附条件下为200 mm.流体中CO₂含量越高，结垢厚度也越大.

关键词

高温地热 / 井筒结垢 / 碳酸钙 / 水文地球化学 / 定量评价 / 水文地质学

Key words

high-temperature geothermal / wellbore scaling / calcite / hydrogeochemistry / quantitative assessment / hydrogeology

0 前言

地热井井筒结垢是高温地热开发普遍遇到的问题之一.比如，国内羊八井地热田开发遇到了较为严重的井筒碳酸钙结垢问题，需要定期采用机械空心通井器除垢（周大吉，2003）；甘孜康定地热井在放喷阶段，工作人员在管线和阀门处发现了大量的碳酸钙结垢，其影响了流体的流动（王延欣等，2015）.国外，美国（McLin et al.， 2012）、土耳其（Demir et al.， 2014）、肯尼亚（Jones and Renaut， 1998）和日本（Fukuyama and Chen， 2021）等高温地热开发也都存在严重的井筒和管线结垢问题.

高温地热生产井一般为碳酸钙结垢，其与温度、压力、流体闪蒸、CO₂析出和流体的化学组分都有非常直接的关系，其涉及到两相流动和水文地球化学过程.Abouie et al.（2017）耦合井筒模拟器UTWELL（Coelho et al.， 2021）和水文地球化学计算模块IPhreeqc（Charlton and Parkhurst， 2011）用于评价油井和地热井筒的结垢问题.Zolfagharroshan and Khamehchi（2020）基于井筒流动模型和地球化学计算程序WATCH（Iceland Water Chemistry Group， 2010）建立了地热井结垢模型，并评价了伊朗西北部地热井的结垢问题.李义曼和庞忠和（2018）、Li et al. （2020）联合井筒流动模拟器WELLSIM（Gunn and Freeston， 1991）、HOLA（Björnsson， 1987）、地球化学模拟器WATCH（Iceland Water Chemistry Group， 2010）和PHREEQC（Parkhurst and Appelo， 2013）定量评价了康定地热井筒碳酸钙结垢.基于井口水气取样测试结果和井筒温压测试，Akın and Kargı（2019）采用PHREEQC评价了地热井筒流动过程中的地球化学演化和可能的结垢趋势.Wanner et al.（2017） 根据井口水样、井中结垢物的化学和矿物特征以及结垢物中的气体包裹体，结合模拟计算评价了德国南部高温地热井筒钙垢的形成原因.

国内，《地热资源地质勘查规范》（GB/T 11615-2010）中提到的拉申指数和雷兹诺指数（以下简称“规范法”），以及通过水文地球化学计算获得的矿物饱和指数法常用来简单评价地热流体的结垢趋势.比如，张恒等（2016）通过地热流体矿物饱和指数法计算评价了康定地热井可能的结垢类型.韦梅华等（2012）采用规范法和热储矿物平衡分析评价了四川康定地热区地热结垢的可能性.于湲等（2007）基于规范法和WATCH的计算结果对北京城区地热田10个代表性地热井进行了结垢趋势判定.

从已有的研究工作可以看到大多数有关高温地热井井筒结垢都只是根据温度、压力和组分评价矿物的饱和状态，缺乏对两相流动、气‒水‒结垢物相互化学作用和颗粒粘附，以及它们之间耦合的定量研究.因此，本文在两相流动的基础上，从水文地球化学的角度分析了井筒结垢过程.首先考虑多个过程建立井筒碳酸钙结垢的定量评价模型，并给出相关的计算流程；然后基于羊八井地热田典型井测试和分析获取评价的关键参数；再基于两相流动的计算结果，结合水文地球化学评价井筒结垢位置和厚度，揭示井筒结垢机理和规律.

1 井筒碳酸钙结垢评价模型

高温地热井井筒碳酸钙结垢从物理过程上涉及到井筒水气两相流动、结垢物过饱和析出的水文地球化学反应和固体在井筒壁上的粘附等多个过程.首先，含有CO₂和多种组分的高温地热流体（单相液体）从储层进入到井筒，再由井底往上流动过程中，压力降低；当压力降低到CO₂‒水体系对应的相变压力时，液相开始闪蒸形成气相，同时伴随温度的降低，这种闪蒸过程一直会进行到井口或是所有液相转化为气相；其次，碳酸钙矿物由于温度和CO₂分压的变化过饱和析出，并逐渐结晶，一部分随着流体往上迁移，一部分粘附到井筒壁上形成结垢物（图1）.井筒壁上的结垢物逐渐堆积会减小井筒直径，严重时会堵塞井筒，对产能产生重大影响.

1.1　井筒流动模型

高温地热生产井流体大部分为H₂O，并含有一定比例的非凝析气体.对碳酸钙结垢有重要影响的为H₂O和CO₂，因此井筒结垢流动过程中至少需要考虑这两个组分.根据质量守恒、能量守恒和动量守恒方程，可以建立CO₂-H₂O体系一维垂直井筒相变流动控制方程.考虑稳定状态，则：

总的质量守恒： \frac{\mathrm{d}}{\mathrm{d}z} \left[{\rho }_{\mathrm{l}}{u}_{\mathrm{l}}\left(\mathrm{1}-\alpha \right)+{\rho }_{\mathrm{g}}{u}_{\mathrm{g}}\alpha \right] = \mathrm{0}，(1)

CO₂质量守恒： \frac{\mathrm{d}}{\mathrm{d}z} \left[{\rho }_{\mathrm{l}}{u}_{\mathrm{l}}\left(\mathrm{1}-\alpha \right)x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}+{\rho }_{\mathrm{g}}{u}_{\mathrm{g}}\alpha y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}\right] = \mathrm{0}，(2)

动量守恒： \frac{\mathrm{d}}{\mathrm{d}z} \left({\dot{m}}_{\mathrm{l}}{u}_{\mathrm{l}}+{\dot{m}}_{\mathrm{g}}{u}_{\mathrm{g}}\right) + A \frac{\mathrm{d}P}{\mathrm{d}z} + {\rho }_{\mathrm{m}} g A + \frac{{\rho }_{\mathrm{m}}fA}{\mathrm{2}d} u_{\mathrm{m}}^{\mathrm{2}} = \mathrm{0}，(3)

能量守恒： \frac{\mathrm{d}}{\mathrm{d}z} [ {\dot{m}}_{\mathrm{l}} \left(\frac{u_{\mathrm{l}}^{\mathrm{2}}}{\mathrm{2}}+gz+h_{\mathrm{l}}\right) + {\dot{m}}_{\mathrm{g}} \left(\frac{u_{\mathrm{g}}^{\mathrm{2}}}{\mathrm{2}}+gz+h_{\mathrm{g}}\right) ] + q = \mathrm{0}，(4)

其中： {\rho }_{\mathrm{l}}和 {\rho }_{\mathrm{g}}分别为液相和气相密度，单位kg/m³； u_{\mathrm{l}}和 u_{\mathrm{g}}分别为液相和气相速度，单位m/s； \alpha为气相体积分数； x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}和 y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}分别为液相和气相CO₂的质量分数； {\dot{m}}_{\mathrm{l}}和 {\dot{m}}_{\mathrm{g}}为液相和气相质量速度，单位kg/s； A为面积，m²； P = P_{\mathrm{H}\mathrm{2}\mathrm{O}} + P_{\mathrm{C}\mathrm{O}\mathrm{2}}为总压力， P_{\mathrm{H}\mathrm{2}\mathrm{O}}和 P_{\mathrm{C}\mathrm{O}\mathrm{2}}分别为H₂O和CO₂的分压，单位Pa； f为井壁摩擦系数，当雷洛数 ， f = \frac{\mathrm{16}}{Re}，否则 \frac{\mathrm{1}}{\sqrt[]{f}} = - \mathrm{4} \mathrm{l} \mathrm{g} [ \frac{\mathrm{2}\epsilon /d}{\mathrm{3.7}} - \frac{\mathrm{5.02}}{Re} \mathrm{l} \mathrm{g} ( \frac{\mathrm{2}\epsilon /d}{\mathrm{3.7}} + \frac{\mathrm{13}}{R\mathrm{e}} ) ]， {\mu }_m为混合物粘度，单位Pa⋅s， \epsilon为井壁粗糙度，单位m； d为井筒直径，单位m； h_{\mathrm{l}}和 h_{\mathrm{g}}分别为液相和气相热焓，单位J/kg； {\rho }_{\mathrm{m}}为混合密度， {\rho }_{\mathrm{m}} = \left(\mathrm{1}-\alpha \right) {\rho }_{\mathrm{l}} + \alpha {\rho }_{\mathrm{g}}，单位kg/m³； u_{\mathrm{m}}为混合物的平均速度， u_{\mathrm{m}} = \frac{\left(\mathrm{1}-\alpha \right){\rho }_{\mathrm{l}}{u}_{\mathrm{l}}+\alpha {\rho }_{\mathrm{g}}{u}_{\mathrm{g}}}{{\rho }_{\mathrm{m}}}，单位m/s； q为井筒与围岩单位长度的热交换量，单位J/（s⋅m）， q = U ( T_{\mathrm{f}} - T_{\mathrm{w}\mathrm{b}} ) , U为集中参数，受井径、以及围岩、套管和固井水泥等热传导系数和比热容综合影响，一般需要采用实际数据校正，单位为W/（ m⋅ ℃）， T_{\mathrm{f}}和 T_{\mathrm{w}\mathrm{b}}分别为井筒流体和井壁的温度，单位为 ℃； z为井筒高程，单位为m； g = \mathrm{9.8}m²/s为重力加速度.有关井筒流动的具体细节请参考雷宏武等（2023）.

1.2　井筒中水-气-结垢物地球化学模型

含有CO₂的流体由井底向井口流动过程中，伴随着压力和温度的降低，CO₂逐渐脱气由液相进入或形成气相，液相中pH升高，碳酸钙矿物过饱和析出形成固体.这个过程涉及到气、液和固体（碳酸钙）之间的水文地球化学平衡，气相中考虑 \mathrm{C} {\mathrm{O}}_{\mathrm{2}(\mathrm{g})}、液相中考虑与碳酸钙结垢相关的离子及络合物，包括： {\mathrm{H}}^{+}、 \mathrm{C} {\mathrm{a}}^{\mathrm{2}+}、 \mathrm{C} \mathrm{a} \mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{+}、 \mathrm{O} {\mathrm{H}}^{-}、 \mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{-}、 \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{\mathrm{2}-}、 \mathrm{C} {\mathrm{O}}_{\mathrm{2}(\mathrm{a}\mathrm{q})}、 \mathrm{C} \mathrm{a} \mathrm{C} {\mathrm{O}}_{\mathrm{3}(\mathrm{a}\mathrm{q})}^{\mathrm{0}}、 \mathrm{C} \mathrm{a} {(\mathrm{O}\mathrm{H})}^{+}，固相中只考虑 \mathrm{C} \mathrm{a} \mathrm{C} {\mathrm{O}}_{\mathrm{3}(\mathrm{s})}，则主要水文地球化学过程为：

气-液平衡： \mathrm{C} {\mathrm{O}}_{\mathrm{2}(\mathrm{a}\mathrm{q})} \leftrightarrow \mathrm{C} {\mathrm{O}}_{\mathrm{2}(\mathrm{g})}，(5)

水溶液组分平衡：

{\mathrm{H}}_{\mathrm{2}} \mathrm{O} + \mathrm{C} {\mathrm{O}}_{\mathrm{2}(\mathrm{a}\mathrm{q})} \leftrightarrow \mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{-} + {\mathrm{H}}^{+},(6)

\mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{-} \leftrightarrow \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{\mathrm{2}-} + {\mathrm{H}}^{+},(7)

{\mathrm{H}}_{\mathrm{2}} \mathrm{O} \leftrightarrow {\mathrm{H}}^{+} + \mathrm{O} {\mathrm{H}}^{-},(8)

\mathrm{C} \mathrm{a} \mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{+} \leftrightarrow \mathrm{C} {\mathrm{a}}^{\mathrm{2}+} + \mathrm{H} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{-}，(9)

\mathrm{C} \mathrm{a} \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{\mathrm{0}} \leftrightarrow \mathrm{C} {\mathrm{a}}^{\mathrm{2}+} + \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{\mathrm{2}-},(10)

\mathrm{C} \mathrm{a} {(\mathrm{O}\mathrm{H})}^{+} \leftrightarrow \mathrm{C} {\mathrm{a}}^{\mathrm{2}+} + \mathrm{O} {\mathrm{H}}^{-},(11)

固-液平衡： \mathrm{C} \mathrm{a} \mathrm{C} {\mathrm{O}}_{\mathrm{3}(\mathrm{s})} \leftrightarrow \mathrm{C} {\mathrm{a}}^{\mathrm{2}+} + \mathrm{C} {\mathrm{O}}_{\mathrm{3}}^{\mathrm{2}-},(12)

需要说明的是虽然其他矿物（比如白云石）也可能过饱和，但是在水文地球化学结垢模型中没有考虑，主要是因为白云石析出成晶过程主要受动力学控制，在井筒这么短的流动时间内来不及析出沉淀.另外，其他离子和相关的络合物也未在这里详细列出，但是在模型中仍然考虑，因为其存在会影响离子活度，从而影响结垢过程.这些化学反应的平衡常数主要来源于TOUGHREACT（Xu et al.， 2006， 2011）自带的数据库，其是温度的函数.

1.3　结垢物井壁粘附模型

结垢物首先在水溶液中过饱和析出，然后粘附到井壁形成结垢，其中主要有3个重要的机制过程（Jamialahmadi et al.， 2009）：（1）结垢物在运移过程中横向运移粘附到井筒的转化过程，该过程与流态、固体析出物颗粒大小、流体的粘度等相关；（2）结垢物粘附到井壁的随机过程（以概率量化）；（3）流体中和井壁面结垢物的浓度梯度.具体粘附模型（Alhosani and Daraboina， 2020）为：

J_{\mathrm{d}} = K_{\mathrm{t}} \left(SP\right) \left(C_{\mathrm{b}}-C_{\mathrm{f}}\right),(13)

式中： J_{\mathrm{d}}为井壁粘附通量，单位为kg/（m²⋅s）； K_{\mathrm{t}}为运移系数，单位为m/s； S P为粘附概率，无单位； C_{\mathrm{b}}和 C_{\mathrm{f}}分别为井壁和流体中结垢物的浓度，单位为kg/m³，一般井壁表面的浓度假设为0.

根据Cleaver and Yates（1975）的模型，运移系数可以表达为：

K_{\mathrm{t}} = \frac{\mathrm{0.084}{v}_{\mathrm{a}\mathrm{v}\mathrm{g}}\sqrt[]{f/\mathrm{8}}}{S{c}^{\mathrm{2}/\mathrm{3}}}，(14)

式中： v_{\mathrm{a}\mathrm{v}\mathrm{g}}为平均速度，单位m/s； f为范宁摩擦因子，无单位，可表达为 f = \frac{\mathrm{0.316}\mathrm{4}}{R{e}^{\mathrm{0.25}}}， R e为雷洛数； S c为施密特数，无单位，可表达为 S c = \frac{\mu }{\rho D_{\mathrm{B}}}， \mu为粘度， \rho为流体密度，布朗扩散率 D_{\mathrm{B}} = \frac{k_{\mathrm{B}}T}{\mathrm{3}\pi \mu d_{\mathrm{p}}}，其中玻尔兹曼常数 k_{\mathrm{B}} = \mathrm{1.38} \times {\mathrm{10}}^{-\mathrm{23}}J/K， T为温度，单位K， d_{\mathrm{p}}为结垢物粒子直径，单位为m，本文取0.5 µm.

根据Watkinson（1970）的定义，粘附概率可表达为：

S P = k_{\mathrm{d}} \frac{e^{-E_{\mathrm{a}}/R{T}_{\mathrm{s}}}}{v_{\mathrm{a}\mathrm{v}\mathrm{g}}^{\mathrm{2}}},(15)

式中： k_{\mathrm{d}}为常数，单位为m²/s²； E_{\mathrm{a}}为粘附活化能，单位为kJ/mol； R为气体常数，为0.008 314 kJ（mol⋅K）； T_{\mathrm{s}}为井壁温度，单位为K. k_{\mathrm{d}}和 E_{\mathrm{a}}一般可以根据实验获取，根据Jamialahmadi et al.（2009）实验结果本文取9.76✕10⁸和65.3.

1.4　组分的运移

流体中的各种物质（包括Na、Ca等元素，以及结垢物）在井筒中会随着流体同时运移，根据物质守恒，针对一维井筒可以得到以下方程：

\frac{\partial [\left(\mathrm{1}-\alpha \right){\rho }_{\mathrm{l}}{C}_k]}{\partial t} = - \frac{\partial \left[C_k{\rho }_{\mathrm{l}}{u}_{\mathrm{l}}\left(\mathrm{1}-\alpha \right)\right]}{\partial z} + q，(16)

式中： C_k为组分 k的浓度，单位kg/kg， q为源汇项，刻画结垢物析出和粘附过程，单位为kg/（m³⋅s）.由于井筒中的流速快，上述方程仅考虑对流过程，不考虑由于浓度梯度引起的扩散过程.在达到稳定状态后，方程（16）可以写为：

对于不参与结垢的元素（比如Na）：

C_k^i = C_k^{i-\mathrm{1}} \frac{{\rho }_{\mathrm{l}}^{i-\mathrm{1}}{u}_{\mathrm{l}}^{i-\mathrm{1}}(\mathrm{1}-{\alpha }_{\mathrm{l}}^{i-\mathrm{1}})}{{\rho }_{\mathrm{l}}^i{u}_{\mathrm{l}}^i(\mathrm{1}-{\alpha }_{\mathrm{l}}^i)}，(17)

对于结垢物(CaCO₃)：C_k^i = [ C_k^{i-\mathrm{1}} {\rho }_{\mathrm{l}}^{i-\mathrm{1}} u_{\mathrm{l}}^{i-\mathrm{1}} \left(\mathrm{1}-{\alpha }_{\mathrm{l}}^{i-\mathrm{1}}\right) - \frac{\mathrm{4}{K}_{\mathrm{t}}^{i-\mathrm{1}}{(SP)}^{i-\mathrm{1}}{C}_k^{i-\mathrm{1}}∆Z}{d} + C_{\mathrm{C}\mathrm{a}}^{i-\mathrm{1}} {\rho }_{\mathrm{l}}^{i-\mathrm{1}} u_{\mathrm{l}}^{i-\mathrm{1}} \left(\mathrm{1}-{\alpha }_{\mathrm{l}}^{i-\mathrm{1}}\right) - C_{\mathrm{C}\mathrm{a}}^i {\rho }_{\mathrm{l}}^i u_{\mathrm{l}}^i \left(\mathrm{1}-{\alpha }_{\mathrm{l}}^i\right) ) ] / [ {\rho }_{\mathrm{l}}^i u_{\mathrm{l}}^i ( \mathrm{1} - {\alpha }_{\mathrm{l}}^i ) ]，(18)

式中： i表示计算节点，由井底向井口编号； ∆ Z表示两个节点之间的间距.对于每个节点的钙浓度，在假设水‒气‒结垢物平衡条件下由化学程序计算.由上式可以看到，对于不参与结垢的元素浓度主要是由浓缩过程控制，对于参与结垢的元素主要由水‒结垢物‒气相互作用控制，而对于结垢物的浓度，则由浓缩、水‒结垢 物‒气相互作用和井壁粘附3个过程控制.在获得结垢物的浓度分布后，由式（13）可以获得每一个节点的结垢速率，从而计算出结垢厚度.

1.5　计算流程

井筒结垢涉及到前面提到的3个过程的耦合，首先由井底到井口流动过程的计算，获得温度、压力、CO₂分压、气相质量分数、气相和液相的速度.然后，进行水文地球化学计算.在未知井底化学组分的情况下，需要首先进行井底流体化学组分的计算.井口流体是在井底流体进行脱气、碳酸钙析出和浓缩3种主要作用下形成的，因此，除了Ca、C和H外，其他元素可以表达为：

C_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}} = \frac{(\mathrm{1}-x_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}})}{(\mathrm{1}-x_{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{d}})} C_{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{d}},(19)

式中： C_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}和 C_{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{d}}分别为井底和井口元素浓度，单位为mg/kg； x_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}和 x_{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{d}}分别为井底和井口气相质量分数.在假设井底化学组分与储层岩石满足化学平衡态时，可以通过在已知CO₂分压和储层温度条件下地热流体与过量的碳酸钙进行平衡而获得井底流体的C和Ca.在获得井底流体化学组分后，便可以由下至上计算流体化学变化和井筒粘附.水文地球化学的计算采用TOUGHREACT（Xu et al.， 2006， 2011）软件进行，该软件已被广泛应用到地热水‒岩‒气相互作用评价中（Dobson et al.， 2004； Xu et al.， 2004， 2014）.井筒碳酸钙结垢计算流程见图2.

2 羊八井ZK304井井筒碳酸钙结垢评价

我国高温地热田主要集中在川藏地区（汪新伟等， 2022），其中的羊八井地热田是我国最早开发的高温地热田，自1977年9月羊八井第一台 1 MW试验机组发电成功，迄今最大装机容量已扩大到25.18 MW.浅部地热田可分为南、北两区：南区是第四系孔隙型热储，北区由第四系孔隙型和喜山期花岗岩裂隙型两类热储组成.浅部热储埋深在地表以下180~280 m，岩性由第四系冲洪积砂砾石层、冰碛砂砾石、基岩顶部花岗岩风化壳组成.热储顶部由厚度不等的泥砾层或粉砂质粘土层构成盖层.热储底部基岩为喜马拉雅山早期花岗岩和凝灰岩，热田北部局部见有糜棱岩化花岗岩.浅层热储流体温度一般在140~160 ℃，流体主要为液态，气体含量较少，水质以Cl·HCO₃-Na型为主，属深部流体与浅表冷水混合产物.

地热田在开发早期井筒就遇到了较为严重的碳酸钙结垢问题.孙本达和袁义方（1987）在羊八井10号井开展了防垢和除垢的试验研究，评价了不同化学阻垢剂的阻垢效果和酸洗除垢的可行性.本文以北区典型井ZK304为对象，通过放喷过程中两相流动、井口取样测试分析和模拟计算，定量评价了井筒结垢过程.

2.1　井筒两相流动特征

根据两相流动和历史拟合的结果（雷宏武等，2023）得到关键认识：（1）ZK304井压力由井底的7.3 bar逐渐下降到井口的2.5 bar，压力的下降主要由3部分贡献：克服重力，克服井壁摩擦和动量的变化，其中克服重力占到了绝对主导的地位；（2）温度由井底的152 ℃下降到井口的128 ℃，在不考虑围岩换热的情况下，温度的下降主要由相变热、CO₂析出热和动能变化引起，其中液体水蒸发为气态的相变吸热占绝对主导作用；（3）闪蒸点在不含CO₂的情况下为埋深97.7 m，在CO₂质量分数0.04%时为107.2 m.在闪蒸点以上，CO₂分压快速下降.具体细节可参考文献雷宏武等（2023），在此不再累述.

2.2　热储流体组分

2.2.1　井口取样和测试

为了获取地热水和气的化学组分，2021年1月在典型地热井的井口进行了采样和测试.采样井ZK324位于地热田南区，ZK304和ZK303位于地热田北区.水样的温度、pH和电导采用Hach仪器在现场直接测量；Ca²⁺、HCO₃ ^-和CO₃ ^2-离子浓度在现场通过滴定法测试.取样、保存和测试方法与郭清海和杨晨（2021）文中描述一致.

从表1可以看到，井口液体的温度约为85 ℃，与当地大气压对应的相变温度一致，说明井口流体处于两相状态.Ca²⁺现场滴定的结果与实验室离子色谱仪测量值基本一致，介于1~2 mg/L.由于压力和温度的变化导致CO₂的解析，HCO₃ ^-和CO₃ ^2-的测量有很大的不确定性.为了分析这种不确定性，本文通过理论计算与测量值进行比较.根据碳酸平衡理论，无论是开放系统还是封闭系统，与碳相关的3种主要组分CO₂（aq）、HCO₃ ^-和CO₃ ^2-含量的比例仅取决于溶液的pH.因此，根据获得的pH比例，结合溶液阴阳离子电荷平衡计算得到溶液中HCO₃ ^-和CO₃ ^2-的含量.从表1可看到，滴定的CO₃ ^2-浓度要明显大于计算值，但碳总量除了ZK304相差较大外（误差25%），另外两口井ZK303和ZK324差异不大，误差分别为4%和0.4%.阳离子中，Na⁺占主导；阴离子中，Cl^-和HCO₃ ^-占主导.三口井流体的化学组分几乎一致，反映了羊八井浅层热储化学的一致性.基于化学组分计算井口典型矿物的饱和指数，可以看到碳酸钙和白云石均过饱和，这也是井筒结垢的必要条件.非凝结气体中，CO₂为主要组分，三口井的含量均超过90%.

在现场井口（图3a）和分离器（图3b）中观察到明显的结垢物.通过XRD分析，碳酸钙的体积占比超过99%，可见井筒结垢主要为钙垢.

2.2.2　井下取样测试

20世纪80年代，在联合国的支援下，中国和意大利技术人员在羊八井地热田开展了高温地热资源的勘探、开发和利用工作，并获取了大量关键参数，包括对井下样品的测试和分析（表2）.测试结果显示现场热储层的平均温度为152.9 ℃，pH为6.13，Ca的浓度为14.49 mg/L.

2.2.3　模拟评价获得储层流体组分

由于未获取井底流体样品，储层流体组分需要根据井口的数据进行反演.由于取样时井口汽化程度较低，不参与结垢的组分（比如Na）可以认为与井口和储层中的一致，而参与结垢的Ca和C浓度需要根据水‒岩‒气平衡获得.根据两相流动计算评价，CO₂的井底含量有一定的不确定性.为了分析这种不确定性，与两相流动一致，设置了两种CO₂含量，最大井底CO₂含量（0.08%）和其一半含量（0.04%），分别对应的井底CO₂分压为2.0和1.0 bar.从图4可以看到，井底Ca的浓度随着CO₂分压升高而增加，由1 bar下的5.56 mg/kg增加到2 bar下的10.4 mg/kg.pH由1 bar下的6.6下降到2 bar下的6.3.与以前的测试结果（表2）相比，计算的Ca偏低，pH偏高，主要原因可能有：（1）在当时测量时，地层中原本的CO₂含量就高于现在，在多年的开采条件下，CO₂含量降低；（2）当时测量时可能有误差，井底样品拿到地面冷凝后，部分悬浮的碳酸钙与流体发生了再平衡.

计算的井底流体化学组分见表3，可以看到，Na、K、Li、Cl和Si分别主要以Na⁺、K⁺、Li⁺、Cl^-和SiO₂（aq）形式存在，Ca元素中以Ca²⁺和CaHCO₃ ⁺为主，Mg元素以MgSO₄（aq）和Mg²⁺为主，C元素主要以CO₂（aq）和HCO₃ ^-为主，S元素以SO₄ ^2-和NaSO₄ ^-为主.CO₂分压的增高，对C元素的组分比例有较大影响，CO₂（aq）的比例由54.6%提高到69.9%，HCO₃ ^-的比例由45.1%降低到29.8%.

2.3　井筒流体地球化学特征

由于井筒碳酸钙结垢受pH影响大，对Ca的浓度有决定性的影响，其他离子只是简单的浓缩，因此，为了节约篇幅，本文仅分析pH和Ca的分布情况.图5显示了两种CO₂含量下的Ca浓度分布，可以看到Ca的分布与井筒CO₂分压（雷宏武等，2023）有极高的相似性，说明水中Ca的浓度受CO₂分压控制.虽然温度的降低有利于碳酸钙的溶解，但这里这种影响基本可以忽略.在井底CO₂分压为1 bar和2 bar的情况下，井口计算的Ca浓度分别为 0.43 mg/kg和0.55 mg/kg，小于井口测量的约 1.5 mg/kg，主要原因是测量条件是在井筒自流条件下，不是在放喷条件下，自流时流速低，围岩换热大大降低了流体的温度，低的温度碳酸钙结垢程度低，水中的Ca浓度也就高.图6显示了溶液的pH分布，随着CO₂的析出和分压的降低，溶液中pH逐渐升高，井口计算值在井底CO₂分压为2.0和1.0 bar情况下分别为8.14和8.39，接近于测量的8.55.

2.4　结垢过程定量评价

在考虑井筒吸附过程条件下，CO₂含量为0.04%时评价的结垢厚度如图7a所示，整体结垢呈现出上下厚度小且不对称，中间厚度大，这与现场观测结果，与现场观测结果（Benoit，1989）具有很高的一致性.最大结垢速率为0.039 7 mm/d，发生在埋深96.9 m，比闪蒸深度107.2 m要浅10 m左右.持续开采1年，最大结垢厚度为14.5 mm.为了评价可能的最大结垢厚度，假设析出的碳酸钙瞬时全部粘附到井筒壁上，计算结果显示（图7b）最大结垢位置与闪蒸位置一致，结垢速率为0.58 mm/d，持续开采1年结垢厚度达到211 mm.

当井底CO₂含量为0.08%时，结果整体趋势与CO₂为0.04%的类似.最大结垢速率为0.067 mm/d，发生在埋深97.5 m.持续开采1年，最大结垢厚度为24.5 mm（图8a），可见结垢厚度随着CO₂的增加而增加.如果析出的碳酸钙全部吸附在井筒壁上，结垢速率为0.56 mm/d，持续开采1年结垢厚度达到205 mm（图8b），与CO₂为0.04%情况非常接近.

3 结论

高温地热井井筒结垢是地热田开发遇到的突出问题之一，严重时会影响地热发电性能.井筒结垢涉及到复杂的井筒多相流动和水文地球化学过程，本文从水文地球化学角度，结合多相流动和结垢粘附，建立了井筒结垢定量评价模型，并结合西藏羊八井地热田开展了典型井的碳酸钙结垢评价.

（1）根据井筒两相流动、水‒气‒结垢物水文地球化学和井壁粘附模型，建立了井筒碳酸钙结垢评价的耦合模型，并开发了相关的计算程序.

（2）羊八井3口地热典型井的流体和结垢物取样和测试分析显示该地区井口流体化学具有较高的相似性，都以Cl·HCO₃-Na水类型为主，典型碳酸盐（碳酸钙和白云石）矿物过饱和，非凝析气体CO₂占绝对主导地位，井筒结垢物为碳酸钙.

（3）基于耦合模型评价了羊八井典型地热井筒碳酸钙结垢，认识到CO₂分压对碳酸钙的析出有控制性的作用，最大结垢厚度位置并不是在闪蒸位置，而是发生在闪蒸面以上10~20 m位置，持续1年的开采最大结垢厚度在考虑井筒粘附动力学时为14~25 mm，假设完全粘附条件下为200 mm.流体中CO₂含量越高，结垢厚度也越大.

本文虽然成功评价了羊八井典型井筒结垢问题，但是仍然存在较多的不确定性：①流体中CO₂含量不确定，这个会大大影响整个的评价结果，需要在下一阶段通过井下保真取样和测试确定；②井筒结垢物的粘附模型中仍然有较多参数无法确定，下一步可以通过相关机理实验获得；③水文地球化学模型中碳酸钙假设为平衡控制的，这种假设是否合理需要通过实验验证；④结垢后井径减小，其对流动的影响需要进一步评价.

参考文献

[1] Abouie, A., Korrani, A. K., Shirdel, M., et al., 2017. Comprehensive Modeling of Scale Deposition by Use of a Coupled Geochemical and Compositional Wellbore Simulator. SPE Journal, 22(4): 1225-1241. https://doi.org/10.2118/185942-pa

[2] Akın, T., Kargı, H., 2019. Modeling the Geochemical Evolution of Fluids in Geothermal Wells and Its Implication for Sustainable Energy Production. Geothermics, 77: 115-129. https://doi.org/10.1016/j.geothermics.2018.09.003

[3] Alhosani, A., Daraboina, N., 2020. Unified Model to Predict Asphaltene Deposition in Production Pipelines. Energy & Fuels, 34(2): 1720-1727. https://doi.org/10.1021/acs.energyfuels.9b04287

[4] Benoit, W. R., 1989. Carbonate Scaling Characteristics in Dixie Valley, Nevada Geothermal Wellbores. Geothermics, 18(1-2): 41-48. https://doi.org/10.1016/0375-6505(89)90008-4

[5] Björnsson, G., 1987. A Multi-Feedzone Geothermal Wellbore Simulator (Dissertation). Lawrence Berkeley Laboratory, Berkeley.

[6] Charlton, S. R., Parkhurst, D. L., 2011. Modules Based on the Geochemical Model PHREEQC for Use in Scripting and Programming Languages. Computers & Geosciences, 37(10): 1653-1663. https://doi.org/10.1016/j.cageo.2011.02.005

[7] Cleaver, J. W., Yates, B., 1975. A Sub Layer Model for the Deposition of Particles from a Turbulent Flow. Chemical Engineering Science, 30(8): 983-992. https://doi.org/10.1016/0009-2509(75)80065-0

[8] Coelho, F. M. C., Sepehrnoori, K., Ezekoye, O. A., 2021. Coupled Geochemical and Compositional Wellbore Simulators: A Case Study on Scaling Tendencies under Water Evaporation and CO₂ Dissolution. Journal of Petroleum Science and Engineering, 202: 108569. https://doi.org/10.1016/j.petrol.2021.108569

[9] Demir, M. M., Baba, A., Atilla, V., et al., 2014. Types of the Scaling in Hyper Saline Geothermal System in Northwest Turkey. Geothermics, 50: 1-9. https://doi.org/10.1016/j.geothermics.2013.08.003

[10] Dobson, P. F., Salah, S., Spycher, N., et al., 2004. Simulation of Water-Rock Interaction in the Yellowstone Geothermal System Using TOUGHREACT. Geothermics, 33(4): 493-502. https://doi.org/10.1016/j.geothermics.2003.10.002

[11] Fukuyama, M., Chen, F. Y., 2021. Geochemical Characteristics of Silica Scales Precipitated from the Geothermal Fluid at the Onuma Geothermal Power Plant in Japan. Journal of Mineralogical and Petrological Sciences, 116(3): 159-169. https://doi.org/10.2465/jmps.201130b

[12] Gunn, C., Freeston, D., 1991. An Integrated Steady-State Wellbore Simulation and Analysis Package. The 13th New Zealand Geothermal Workshop, Auckland.

[13] Guo, Q. H., Yang, C., 2021. Tungsten Anomaly of the High-Temperature Hot Springs in the Daggyai Hydrothermal Area, Tibet, China. Earth Science, 46(7): 2544-2554 (in Chinese with English abstract).

[14] Iceland Water Chemistry Group, 2010. The Chemical Speciation Program WATCH, Version 2.4. website: ÍSOR - Iceland GeoSurvey, Reykjavik. www.geothermal.is/ software.

[15] Jamialahmadi, M., Soltani, B., Müller-Steinhagen, H., et al., 2009. Measurement and Prediction of the Rate of Deposition of Flocculated Asphaltene Particles from Oil. International Journal of Heat and Mass Transfer, 52(19-20): 4624-4634. https://doi.org/10.1016/j.ijheatmasstransfer.2009.01.049

[16] Jones, B., Renaut, R. W., 1998. Origin of Platy Calcite Crystals in Hot-Spring Deposits in the Kenya Rift Valley. Journal of Sedimentary Research, 68(5): 913-927. https://doi.org/10.2110/jsr.68.913

[17] Lei, H. W., Bai, B., Cui, Y. X., et al., 2023. Quantitative Assessment of Calcite Scaling of a High Temperature Geothermal Production Well: Two-Phase Flow—Application to the Yangbajing Geothermal Fields, Tibet. Earth Science, 48(3): 923-934 (in Chinese with English abstract).

[18] Li, Y. M., Pang, Z. H., 2018. Carbonate Calcium Scale Formation and Quantitative Assessment in Geothermal System. Advances in New and Renewable Energy, 6(4): 274-281 (in Chinese with English abstract).

[19] Li, Y. M., Pang, Z. H., Galeczka, I. M., 2020. Quantitative Assessment of Calcite Scaling of a High Temperature Geothermal Well in the Kangding Geothermal Field of Eastern Himalayan Syntax. Geothermics, 87: 101844. https://doi.org/10.1016/j.geothermics.2020.101844

[20] McLin, K. S., Moore, J. N., Bowman, J. R., et al., 2012. Mineralogy and Fluid Inclusion Gas Chemistry of Production Well Mineral Scale Deposits at the Dixie Valley Geothermal Field, USA. Geofluids, 12(3): 216-227. https://doi.org/10.1111/j.1468-8123.2012.00363.x

[21] Parkhurst, D. L., Appelo, C. A. J., 2013. Description of Input and Examples for PHREEQC Version 3: A Computer Program for Speciation, Batch-Reaction, One-dimensional Transport, and Inverse Geochemical Calculations. U.S. Geological Survey, Denver.

[22] Sun, B. D., Yuan, Y. F., 1987. Study on Preventing Scaling and Descaling of Geothermal Fluid. Thermal Power Generation, 16(4): 15-19 (in Chinese).

[23] Wang, X. W., Wang, T. H., Gao, N. A., et al., 2022. Formation Mechanism and Development Potential of Geothermal Resource along the Sichuan-Tibet Railway. Earth Science, 47(3): 995-1011 (in Chinese with English abstract).

[24] Wang, Y. X., Liu, S. L., Bian, Q. Y., et al., 2015. Scaling Analysis of Geothermal Well from Ganzi and Countermeasures for Anti-Scale. Advances in New and Renewable Energy, 3(3): 202-206 (in Chinese with English abstract).

[25] Wanner, C., Eichinger, F., Jahrfeld, T., et al., 2017. Causes of Abundant Calcite Scaling in Geothermal Wells in the Bavarian Molasse Basin, Southern Germany. Geothermics, 70: 324-338. https://doi.org/10.1016/j.geothermics.2017.05.001

[26] Watkinson, A. P., 1970. Particulate Fouling of Sensible Heat Exchangers. University of British Columbia, Vancouver.

[27] Wei, M. H., Tian, T. S., Sun, Y. D., 2012. A Study of the Scaling Trend of Thermal Groundwater in Kangding County of Sichuan. Hydrogeology ＆ Engineering Geology, 39(5): 132-138 (in Chinese with English abstract).

[28] Xu, T. F., Feng, G. H., Shi, Y., 2014. On Fluid-Rock Chemical Interaction in CO₂-Based Geothermal Systems. Journal of Geochemical Exploration, 144: 179-193. https://doi.org/10.1016/j.gexplo.2014.02.002

[29] Xu, T. F., Ontoy, Y., Molling, P., et al., 2004. Reactive Transport Modeling of Injection Well Scaling and Acidizing at Tiwi Field, Philippines. Geothermics, 33(4): 477-491. https://doi.org/10.1016/j.geothermics.2003.09.012

[30] Xu, T. F., Sonnenthal, E., Spycher, N., et al., 2006. TOUGHREACT-A Simulation Program for Non-Isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media: Applications to Geothermal Injectivity and CO₂ Geological Sequestration. Computers & Geosciences, 32(2): 145-165. https://doi.org/10.1016/j.cageo.2005.06.014

[31] Xu, T. F., Spycher, N., Sonnenthal, E., et al., 2011. TOUGHREACT Version 2.0: A Simulator for Subsurface Reactive Transport under Non-Isothermal Multiphase Flow Conditions. Computers & Geosciences, 37(6): 763-774. https://doi.org/10.1016/j.cageo.2010.10.007

[32] Yu, Y., Zhou, X., Fang, B., 2007. Judgement and Analysis of the Scaling Trend of Thermal Groundwater in Beijing’s Urban Geothermal Fields. City Geology, 2(2): 14-18 (in Chinese with English abstract).

[33] Zhang, H., Hu, Y. Z., Yun, Z. H., et al., 2016. Applying Hydro-Geochemistry Simulating Technology to Study Scaling of the High-Temperature Geothermal Well in Kangding County. Advances in New and Renewable Energy, 4(2): 111-117 (in Chinese with English abstract).

[34] Zhou, D. J., 2003. Operation, Problems and Countermeasures of Yangbajing Geothermal Power Station in Tibet. Electric Power Construction, 24(10): 1-3, 9 (in Chinese with English abstract).

[35] Zolfagharroshan, M., Khamehchi, E., 2020. A Rigorous Approach to Scale Formation and Deposition Modelling in Geothermal Wellbores. Geothermics, 87: 101841. https://doi.org/10.1016/j.geothermics.2020.101841

[36] 郭清海, 杨晨, 2021. 西藏搭格架高温热泉中钨的水文地球化学异常. 地球科学, 46(7): 2544-2554.

[37] 雷宏武, 白冰, 崔银祥, 等, 2023. 高温地热生产井碳酸钙结垢定量评价：两相流动——以西藏羊八井为例. 地球科学, 48(3): 923-934.

[38] 李义曼, 庞忠和, 2018. 地热系统碳酸钙垢形成原因及定量化评价. 新能源进展, 6(4): 274-281.

[39] 孙本达, 袁义方, 1987. 防止地热流体结垢和除垢的研究. 热力发电, 16(4): 15-19.

[40] 汪新伟, 王婷灏, 高楠安, 等, 2022. 川藏铁路沿线地热资源形成机理与开发潜力, 地球科学, 47(3): 995-1011.

[41] 王延欣, 刘世良, 边庆玉, 等, 2015. 甘孜地热井结垢分析及防垢对策. 新能源进展, 3(3): 202-206.

[42] 韦梅华, 田廷山, 孙燕冬, 等, 2012. 四川省康定地区地热水结垢趋势分析. 水文地质工程地质, 39(5): 132-138.

[43] 于湲, 周训, 方斌, 2007. 北京城区地下热水结垢趋势的判断和分析. 城市地质, 2(2): 14-18.

[44] 张恒, 胡亚召, 云智汉, 等, 2016. 水文地球化学模拟技术在康定某高温地热井结垢研究中的应用. 新能源进展, 4(2): 111-117.

[45] 周大吉, 2003. 西藏羊八井地热发电站的运行、问题及对策. 电力建设, 24(10): 1-3, 9.

基金资助

国家自然科学基金项目(41972316)

中国核工业集团有限公司集中研发项目“地热勘查开发利用关键技术研究”

四川省科技创新人才项目(2022JDRC0027)