基于等效渗流通道模型的地热尾水回灌理论模型

赵志宏，刘桂宏,2，谭现锋，张平平

(1.清华大学土木工程系，北京 100084；2.中国矿业大学力学与建筑工程学院，江苏 徐州 221116；3.中国地质大学(武汉)工程学院，湖北 武汉 430074；4.山东省鲁北地质工程勘察院，山东 德州 253072)

基于等效渗流通道模型的地热尾水回灌理论模型

赵志宏1，刘桂宏1,2，谭现锋3,4，张平平4

(1.清华大学土木工程系，北京 100084；2.中国矿业大学力学与建筑工程学院，江苏 徐州 221116；3.中国地质大学(武汉)工程学院，湖北 武汉 430074；4.山东省鲁北地质工程勘察院，山东 德州 253072)

深层地热能是一种清洁环保的可再生能源，尾水回灌是深层地热可持续开发利用的重要保证，但目前国内地热尾水回灌还处在起步阶段，地热尾水回灌方案优化设计的理论模型尚不成熟。文章结合示踪试验，提出了基于等效渗流通道模型的热储参数反演与开采井热突破预测的完整理论框架。应用多孔介质溶质运移理论，修正了等效渗流通道中的溶质浓度解析解，并联合移动渐近线法提出了热储参数反演理论模型；推导了等效渗流通道中的对流传热解析模型，可对开采井中的水温变化进行预测。依托山东省德州市平原县魏庄社区地热对井示踪试验，应用该理论框架进行了示踪试验反演及开采井热突破预测，并对影响开采井热突破的主要因素进行了敏感性分析。

尾水回灌；砂岩热储；示踪试验；参数反分析；热突破

深层地热资源是一种清洁环保的可再生能源，而地热尾水回灌是深层地热可持续开发利用的关键环节。目前，能准确指导地热尾水回灌方案优化设计的理论模型尚不成熟，相关工作多依据经验进行，故长期尾水回灌造成开采井热突破的风险较高。准确掌握回灌井与开采井之间的水力连通情况是地热尾水回灌科学管理的前提，但由于埋深和现场测量困难等限制，直接探明地热井间的水力连通情况还比较困难。基于示踪试验的水文地质参数反分析是解决上述问题的有效手段之一，并已在地热田管理中得到了应用[1]。

Axelsson等[2]利用碘化钾和荧光素钠在冰岛的Laugaland地热田进行了三次示踪试验，采用等效渗流通道模型分析了采灌井之间的连通程度，并对开采井的热突破以及开采能量进行了预测。曾梅香等[3]利用碘化钾和放射性同位素35S在天津市王兰庄地热田进行了示踪试验，通过分析示踪剂回收曲线的峰值浓度，计算了回灌流体的“概化”速度，初步判断了采灌井之间热储层的连通性。庞菊梅等[4]利用荧光素钠在雄县地热田进行了示踪试验，分析了回灌井与多口开采井之间存在多渗流通道的情况，并对通道距离与渗流速度进行了反分析。李元杰[5]利用钼酸铵和碘化钾在北京城区地热田进行了示踪试验，结合抽水试验数据，构造了高斯-塞德尔迭代公式进行参数反演计算，并对热储温度场进行了模拟研究。王树芳等[6]利用荧光素钠在雄县地热田进行了示踪试验，计算了示踪剂在通道中的最大流速以及地热水的流向，并运用水平裂隙型介质模型分析了不同对井间距对热突破的影响。

以上研究在地热田示踪试验方面积累了丰富的经验，但对于示踪剂在热储渗流通道中的运移规律及开采井热突破预测等理论模型方面的研究尚不完善。本文结合现场示踪试验，提出了基于等效渗流通道模型的热储参数反演与开采井热突破预测的完整理论框架。依托山东鲁北地质工程勘察院在德州市平原县魏庄社区的地热对井系统，应用该理论框架对砂岩热储地热尾水回灌进行了示踪试验反演及开采井热突破预测，并对影响开采井热突破的主要因素进行了敏感性分析。

1 理论模型

1.1 等效渗流通道模型

等效渗流通道模型假设地热对井之间存在若干优势渗透通道，其主导了热储中的地热水渗流(图1)。优势渗流通道可能是近乎垂直的断裂带的一部分，或是水平岩层的一部分，并且优势流动通道内的渗流可近似等效为一维达西流[1]。

图1 等效渗流通道模型示意图Fig.1 Sketch of the equivalent flow channel model

1.2 溶质运移理论

示踪剂在一维等效渗流通道中的主要运移机理为对流和水动力弥散。若忽略通道与其它热储岩石间的扩散过程及通道内的横向水动力弥散作用，等效渗流通道内的示踪剂运移过程可由经典的多孔介质中对流-扩散方程来描述[7]：

(1)

式中：C——等效渗流通道中的示踪剂浓度/(kg·m-3)；

t——时间/s；

x——通道轴向坐标/m；

D——通道轴向弥散系数/(m2·s-1)；

u——通道中的平均流速/(m·s-1)。

在t=0时刻从回灌井瞬时投入质量为M的示踪剂，其他初始条件和边界条件可以表示为：

C(x,y,t=0)=C0;C(x=0,y,t)=0;

C(x,y=∞,t)=0

(2)

式中：y——通道纵坐标/m。

结合初始条件和边界条件，偏微分方程(1)的解为：

(3)

其中：

(4)

D=αLu

(5)

式中：Mc——进入渗流通道中的示踪剂质量/kg；A——渗流通道的截面积/m2；φ——孔隙度；qc——渗流通道中的回灌率/(kg·s-1)；ρ——水的密度/(kg·m-3)；αL——纵向弥散度/m。

考虑到开采井中示踪剂的质量守恒，在开采率为Q(kg·s-1)的情况下，有下列关系成立：

C(t)·qc=c(t)·Q

(6)

开采井中示踪剂的浓度表达式为：

(7)

式中：c——开采井中示踪剂的浓度/(kg·m-3)。

Axelsson[1,8]曾提出过以下相似的表达式：

(8)

但上式左右两端的量纲不统一。庞菊梅等[7]给出的开采井中示踪剂浓度表达式为：

(9)

其将式(7)中x/t一项替换成了u，虽然对量纲没有影响，但式(8)中的x与t均为变量，与通道中的平均流速u不能相互替换。关于式(7)与式(9)在应用方面的具体差别，将在2.2节中详述。

如果地热对井之间存在有n条流动通道，开采井中示踪剂的浓度表达式为：

(10)

其中：

(11)

(12)

Di=αLiui

(13)

(14)

式中：xi——第i条渗流通道的轴向坐标/m；qi——第i条渗流通道中的回灌率/(kg·s-1)；qin——回灌率/(kg·s-1)；ui——第i条渗流通道中的流速/(m·s-1)；αLi——第i条渗流通道中的纵向弥散度/m；Ai——第i条渗流通道的截面积/m2；Mi——第i条渗流通道中示踪剂的质量/kg。

1.3 反演模型

式(10)中的未知参数(如通道长度xi、平均流速ui、弥散度αi以及通道截面积Ai等)无法直接测定，需通过示踪试验反演求得。本文引入移动近似法(The Method of Moving Asymptotes，MMA)建立热储参数反分析理论框架。目标函数即原始优化问题P可定义为：

Subjectto:h1(x)≥0

(15)

拟合优度R2是衡量拟合程度的重要指标，其定义如下[9]：

(16)

其中

(17)

R2的取值范围为[0,1]，R2的值越接近1，说明拟合程度越好；反之，R2的值越接近0，说明拟合程度越差。

MMA算法由Svanberg提出[10]，特别适合应用于多变量的优化问题。MMA算法需要先指定优化变量的上下界，之后按照如下步骤进行优化计算：

①选择初始计算点x(1)。

②计算在该迭代点处的目标函数和约束函数各自的函数值fi(x(k))和梯度值▽fi(x(k))。

④求解子问题P(k)，并令所得到的最优解成为下一个迭代点x(k+1)，然后重复步骤②～④。

⑤满足收敛判据时，迭代停止。计算所得子问题P(k)的最优解就是原始问题P的唯一最优解。

Svanberg[11]给出了如下收敛判据的形式：

(18)

其中，ε对于所有的j=1,…,n是一个极小值。

Axelsson[1]曾开发了基于非线性最小二乘法的示踪试验分析软件TRINV来进行示踪试验反分析，但对复杂问题非线性最小二乘法易出现迭代不收敛的情况，而MMA算法的优势在于可同时优化多个变量，对优化变量的上下界不敏感，且极少出现迭代不收敛的情况。

1.4 对流传热理论

假设水-岩界面温度相等，等效渗流通道中的对流换热控制微分方程为：

(19)

其中：

(ρc)f=ρwcwφ+ρrcr(1-φ)

(20)

式中：T——温度/℃； (ρc)f——渗流通道中材料的体积热容/(J·m-3·℃-1)；

kr——热储岩石的导热系数/(W·m-1·℃-1);

b——渗流通道的宽度/m；

ρw——水的密度/(kg·m-3)；

ρr——热储岩石密度/(kg·m-3)；

cw——水的质量热容/(J·kg-1·℃-1)；

cr——热储岩石的质量热容/(J·kg-1·℃-1)。

式(19)可进一步化简为如下形式：

(21)

令

(22)

式中：h——渗流通道的高度/m。

则式(21)可化简为：

(23)

当h>>b时，初始条件和边界条件为：

T(x,y,t=0)=T0;T(x=0,y,t)=Tin;

T(x,y=∞,t)=T0

(24)

所以，渗流通道中水的温度表达式为：

T(x,t)=T0+(Tin-T0)·

(25)

其中，erfc()为误差补函数，它与误差函数erf()的关系为：

(26)

所以，式(25)可以写为

(27)

Axelsson[1,8]等人曾提出过以下相似的表达式：

(28)

其中，α为热扩散系数，可具体表示为：

(29)

所以，式(28)最后可以化为：

(30)

笔者认为，(ρc)f表示流动通道中材料的体积热容，这其中的材料不仅有水，还有多孔的热储岩石，因此，kr/(ρc)f这一项其实是考虑了流动通道中材料的热扩散系数。而kr/ρrcr只是考虑了热储岩石的热扩散系数，整个热储层是由热储岩石和地下水共同组成的，因此，在理论上，考虑整个热储层的热扩散系数要更严谨一些。关于式(27)和式(30)的具体区别，将在2.4节详述。

考虑到开采率为Q，则在不同时刻开采井中的水温Tp为：

(31)

如果两井之间有i条流动通道，则开采井中的水温Tp可以通过以下表达式计算：

(32)

2 工程实例

2.1 示踪试验

试验场地位于山东省德州市平原县魏庄社区，取水层为新近系馆陶组热储，开采井深1 450 m，回灌井深1 400 m，两井直线距离231 m。2015年1月20日至3月16日，采用钼酸铵为示踪剂在该地热对井系统开展了示踪试验。2015年1月20日在回灌井一次性投入钼酸铵50 kg，从2015年1月26日起，每天取两件水样；随着时间的推移，取样频率逐渐加密，至2月17日，每两小时取样一件；后期随着示踪剂含量的降低，取样频率减小；至2015年3月16日试验结束，共收集水样273件。测试数据整理后形成了地热水中钼含量与时间关系散点分布图(图2)。钼元素的浓度峰值大概在750 h左右到达开采井，由于本次试验是生产性回灌，供暖需于2015年3月15日24点停暖，不再进行地热回灌试验，示踪试验也到此终止。但至供暖结束时钼含量有升高迹象，表明回灌水中钼酸铵第二次到达了开采井中。

2.2 参数反演结果与分析

由钼元素的两个浓度峰值推测该地热对井之间存

图2 示踪试验数据和拟合曲线对比图Fig.2 Profiles showing model data(open circles)for tracer testing and model fits for Eqs. (7) and (9)

在2条优势渗流路径，参数反演中设定的参数取值范围如表1所示。分别采用式(7)与式(9)对示踪试验结果进行拟合。综合表1与图2可以看出：采用式(7)得到的拟合曲线的拟合优度R2=0.613>0.6，拟合程度较好。两条渗流通道的长度约为241 m和304 m，说明该地热对井之间存在一条较为直接的渗流路径和另外一条较为曲折的渗流路径。长度较大的渗流通道内水流速度和弥散度都较小，但通道截面积较大。采用式(9)得到的拟合曲线与采用式(7)得到的拟合曲线相近似，参数优化结果与式(7)相比，最大差别不超过4%，但计算所得式(9)曲线的拟合优度R2=0.577<0.613，这说明式(7)不仅在理论上比式(9)更严谨，而且也能更好地用于示踪试验反演。

表1 参数取值范围及反演结果Table 1 Best-fit parameters for the equivalent flow channels between geothermal wells

2.3 热突破预测结果与分析

回灌的“冷水”通过渗流通道被加热，这取决于渗流通道的表面积而不是它们的体积，因此，必须对流动通道的几何特征做出一定的假设。考虑到初始条件和边界条件存在的前提是h>>b，结合示踪试验反分析出的参数取值，对开采井中水的温度变化进行预测，相关参数取值如表2所示。

通过分析示踪试验数据可以得到示踪剂的回收率Mi/M，计算所得示踪剂回收率为0.4323%，计算所得开采水温在100 a内的变化情况如图3所示。本文采用的热突破定义为：回灌引起开采井水温降低[8]。示踪剂回收率较低表明采灌井之间的水力联系较差，开采井热突破时间约为20 a。基于目前的研究预测，该地热对井系统在长期回灌条件下开采水温降低并不明显，100 a之后仅降低了0.5 ℃。但是，需指出以上热突破预测的正确性还需要结果开采井水温的监测数据进行校正，由于该地热对井运行时间只有2 a，相关数据很多，这方面的研究留待今后进行。

表2 参数取值列表Table 2 List of parameter values

图3 开采井温度变化曲线Fig.3 Temperature change in the production well

2.4 参数敏感性分析

为了对比式(27)和式(30)的区别，同时讨论渗流通道宽度b、回灌量q和对井间距x的大小对热突破的影响程度。下面将用通道1的相关数据对以上参数进行敏感性分析，为进一步优化地热回灌提供理论依据。

(1)计算当b取0.1 m、0.3 m、0.5 m，对应的h取6.4 m、2.13 m、1.28 m时开采井中温度的变化情况，结果见图4。运用式(27)和式(30)所计算的开采井水温变化曲线的趋势基本一致，式(27)所计算的温度变化幅度略小于式(30)的结果，最大差值在0.15 ℃左右，但从理论上而言，式(27)比式(30)要更严谨一些。此外，还可以看出，随着b的不断增大，开采温度下降幅度越大，这是由于该模型只考虑上下平板的热传导作用，当b越大，上下平板的间距就越大，在相同时间内通过渗流通道的水被加热得就越不充分，开采水温下降的幅度就越大。

图4 开采井中温度随时间变化规律Fig.4 Prediction of temperature in the production well

(2)开采量保持65 m3/h(18.06 kg/s)不变，回灌量分别取值20 m3/h(5.56 kg/s)、40 m3/h(11.12 kg/s)、60 m3/h(16.68 kg/s)时开采井温度在10 a内的变化曲线如图5所示。随着回灌量的不断增大，开采井温度下降幅度也在不断增大(表3)。这是由于采灌井之间流动通道中的回灌水流速大小主要由回灌量大小决定，由式(22)、式(27)可知，在开采量保持不变的情况下，回灌量q越大，则β越大，因此，t就越小，回灌水到达开采井的时间就越短，换热时间越短，故温度下降幅度也就越大。

图5 不同回灌量下开采井温度变化曲线Fig.5 Changes in temperature in the production well under different reinjection rates

q/(m3·h-1)204060ΔT/℃46110175

(3)开采量保持不变65 m3/h(18.06 kg/s)，回灌量保持20 m3/h(5.56 kg/s)不变，对井间距x分别取值300 m、600 m、900 m时开采井温度在10 a内的变化曲线见图6。随着对井间距的不断增大，开采温度变化的幅度在不断减小。根据式(27)得出，在β保持不变的前提下，x值越大，则t越大，因此，回灌水到达开采井的时间就越慢，在流动通道中被加热的时间越长，最后开采温度变化幅度也就越小(表4)。

图6 不同对井间距下开采温度变化曲线Fig.6 Changes in temperature in the production well considering different well distances

x/m300600900ΔT/℃422311

3 结论

(1)示踪试验反演分析中得到的拟合优度R2=0.613>0.6，说明优化结果的拟合程度较好；运用修正前的公式进行相同的试验数据进行反分析，计算所得的R2=0.577<0.6，说明本文提出的修正模型可更好地运用于示踪试验反分析。

(2) 示踪剂回收率较低表明采灌井之间的水力联系较差，本文初步计算开采井热突破时间约为20 a，但该地热对井系统在长期回灌条件下开采水温降低并不明显。

(3)流通道宽度b越大，开采井温度下降的幅度就越大，故在热突破预测中b的取值要合理；回灌量q越大，回灌地热尾水到达开采井的时间就越短，开采井地热水温度下降的幅度就越大，故在实际工程中砂岩热储回灌量不宜过大，可采用“一采两灌”或“两采多灌”的砂岩热储地热尾水回灌模式；对井间距x越大，回灌水到达开采井的时间就越长，开采温度下降的幅度就越小，故在实际工程中应尽量设计较大的采灌井间距。

[1] Axelsson G, Björnsson G, Montalvo F. Quantitative interpretation of tracer test data[C]//Proceedings World Geothermal Congress. Antalya, Turkey, 2005:24-29.

[2] Axelsson G, Flovenz O G, Hauksdottir S.etal. Analysis of tracer test data, and injection-induced cooling, in the Laugaland geothermal field, N-Iceland[J]. Geothermics,2001, 30(6): 697-725.

[3] 曾梅香,阮传侠,田光辉. 井间示踪技术在岩溶裂隙热储层采灌系统中的应用——以天津市王兰庄地热田回灌井HX-25B示踪测试为例[J].地质与勘探,2008, 44(2):105-109.[ZENG M X，RUAN C X，TIAN G H,etal. Connect test between karst cranny reservoir pumping well and injection well[J]. Geology and Prospecting, 2008,44(2):105-109.(in Chinese)]

[4] 庞菊梅,庞忠和,孔彦龙,等. 岩溶热储井间连通性的示踪研究[J]. 地质科学,2014,49(3):915-923.[PANG J M,PANG Z H, KONG Y L,etal. Interwell connectivity in a karstic geothermal reservoir through tracer tests[J]. Chinese Journal of Geology,2014,49(3):915-923.(in Chinese)]

[5] 李元杰. 地热回灌示踪技术及热储模拟实验研究[D]. 北京:中国地质科学院,2010.[LI Y J. Study on geothermal reinjection of tracer technology and simulation of thermal energy storage[D]. Beijing: Chinese Academy of Geological Sciences, 2010.(in Chinese)]

[6] 王树芳,刘久荣,林沛,等. 岩溶热储回灌实验与示踪试验研究[J]. 水文地质工程地质,2013,40(6):129-133. [WANG S F,LIU J R, LIN P,etal. A study of reinjection experiment and tracer test in a karst geothermal reservoir[J]. Hydrogeology & Engineering Geology, 2013,40(6):129-133.(in Chinese)]

[7] Bodin J, Delay F, Marsily G. Solute transport in a single fracture with negligible matrix permeability: 2. mathematical formalism[J]. Hydrogeology Journal, 2003,11(4):434-454.

[8] Axelsson G, Björnsson G, Flovenz O G,etal. Injection experiments in low-temperature geothermal areas in Iceland[C]//Proceedings of the World Geothermal Congress. Florence, Italy,1995:1991-1996.

[9] Zhao Z H, Skelton A. Simultaneous calculation of metamorphic fluid fluxes, reaction rates and fluid-rock interaction timescales using a novel inverse modeling framework[J]. Earth and Planetary Science Letters, 2013, 373:217-227.

[10] Svanberg K. The method of moving asymptotes-a new method for structural optimization[J]. International Journal for Numerical Methods in Engineering, 1987,24(2):359-373.

[11] Svanberg K. A class of globally convergent optimization methods based on conservative convex separable approximations[J]. SIAM Journal on Optimization, 2002, 12(2):555-573.

责任编辑：汪美华

Theoretical model of geothermal tail water reinjection based onthe equivalent flow channel model

ZHAO Zhihong1，LIU Guihong1,2，TAN Xianfeng3,4，ZHANG Pingping4

(1.TsinghuaUniversity，DepartmentofCivilEngineering，Beijing100084，China；2.ChinaUniversityofMiningandTechnology，SchoolofMechanicsandCivilEngineering，Xuzhou,Jiangsu221116，China；3.ChinaUniversityofGeosciences(Wuhan)，FacultyofEngineering，Wuhan，Hubei430074，China；4.ShandongProvincialLubeiGeo-engineeringExplorationInstitute,Dezhou,Shandong253072，China)

Deep geothermal resources are clean, environmentally friendly and renewable, and geothermal tail water reinjection plays an important role in the sustainable management of deep geothermal reservoirs. However, the progress of geothermal tail water reinjection is slow in China, and theoretical models guiding geothermal tail water reinjection is still under development. Based on tracer tests and equivalent flow channel models, this study proposes a theoretical modeling framework in order to back-calculate the unknown parameters of geothermal reservoirs and predict the thermal breakthrough in the production well. The inverse modeling module combines the analytical solution of the tracer transport equation with the method of moving asymptotes, and an analytical model considering heat advection and transfer can be used to predict the thermal breakthrough in the production well. Based on the results of tracer tests in the double geothermal wells near Dezhou of Shandong, two equivalent flow channels were determined and their average flow velocity, dispersion coefficient, section area and channel length are estimated by using this theoretical framework. The possible thermal breakthrough in the production well and its main influencing factors are also discussed.

tail water reinjection; sandstone reservoir; tracer test; inverse modelling; thermal breakthrough

2016-05-26;

2016-07-28

国家自然科学基金项目(51509138)；北京市自然科学基金项目(8152020)；北京市优秀人才培养资助(201400002012G115)

赵志宏(1983-)，男，博士，讲师，主要从事岩石力学与地下工程方面的研究与教学工作。E-mail：zhzhao@tsinghua.edu.cn

谭现锋(1977-)，男，研究员，在职博士研究生，主要从事地热、岩盐、页岩气、干热岩的勘探开发及相关研究工作。E-mail：geotan1977@126.com

10.16030/j.cnki.issn.1000-3665.2017.03.23

P314.1

A

1000-3665(2017)03-0158-07