地热井井间干扰是指某一地热井抽(放)水时，对相邻地热井的水量、水位及水温等产生影响的现象。重庆作为著名的温泉之都，其地热资源丰富，品质优良，地热水资源开发带来的经济效益也令人瞩目。然而，在地热水资源勘查开发中，井间干扰现象普遍存在，其干扰程度如何？合理的井间距是多少？目前鲜有研究，这对合理布井和开发也造成了一定的困难。由于地热水主要赋存于深部岩溶含水介质中，具有高度的非均质性和各向异性，定量刻画其地下水水流特征，便成为多年来水文地质工作的热点与难点，因此，开展地热水水流特征模拟研究，定量评估地热井间的干扰程度，具有重要的科学与现实意义。

Visual Modflow软件是由加拿大Waterloo Hydrogeologic Inc在美国地质调查局Modflow软件(1984年)的基础上应用可视化技术开发研制的，其基本原理就是应用有限差分法对渗流场进行离散求解，得到离散点上的近似值，它已广泛应用于模拟二维、三维饱和流状态下的水流和污染物运移，如刘淑琴等^[1]以挥发酚为例，建立了地下水流运动和污染物迁移数学模型，利用Visual Modflow中的Modflow和Mt3dms模块求解，对气化炉闭炉后的污染物迁移规律进行了数值模拟，结果拟合性良好；朱常春等^[2]基于Visual Modflow软件分析了相邻矿山干扰开采对地下水流场的影响；刘永良等^[3]基于Visual Modflow软件准确预测鹤壁矿区寺湾矿矿井涌水量和水位降深。

目前，针对孔隙介质，地下水数值模拟较为成熟，而针对管道-裂隙-空隙三重介质并存的岩溶含水系统，其应用相对较少。不少学者在系统概化和模块改进等方面进行了尝试与探索，并取得了一定的成效：1998年，成建梅等^[4]以广西环江北山岩溶含水系统为例，建立岩溶管道-裂隙-孔隙三重介质地下水模型；2013年，Gallegos等^[5]以室内水箱模型为基础，应用Modflow-CFP探索岩溶管道水流模拟方法；2014年，肖斌等^[6]基于Modflow对岩溶管道概化进行探讨，从理论上阐述不同模块概化岩溶管道的优缺点；赵良杰等^[7]以Modflow为基础，以广西寨底岩溶地下河系统为实例，对比分析不同方法概化岩溶管道结构，认为相对Drain模块而言，应用River模块概化模拟岩溶管道更加精确。目前国内对深层地热水数值模拟研究较少，为此，笔者选用地热勘查程度较高的重庆观音峡背斜北段东翼静观片区作为研究对象，尝试利用Visual Modflow软件建立该片区的地下水流数值模拟模型，预测分析不同开采方案下的井间干扰程度，以期为地热资源勘查合理布井提供依据。

1 研究区概况

研究区位于重庆市北碚区东北部，嘉陵江东岸，面积45 km²。地势总体趋势北东高，向南西延伸逐渐降低，属亚热带季风气候，极端最高气温42 ℃左右，极端最低气温-2 ℃左右，年均降雨量1 000 mm以上。本区属于长江水系，嘉陵江为长江流域主要支流，切割深度400~500 m，为本区最低侵蚀基准面。本区处于观音峡背斜北段东翼，主要发育第四系(Q)、侏罗系(J)、三叠系(T)地层(图 1)。区内地热水热储层为三叠系下统嘉陵江组(T₁j)和中统雷口坡组(T₂l)，热储盖层为三叠系上统须家河组(T₃xj)和侏罗系(J)，热储下部隔水层为三叠系下统飞仙关组(T₁f)，其中主要热储层为三叠系下统嘉陵江组第二、三、四段(T₁j²⁺³⁺⁴)，次要热储层为三叠系下统嘉陵江组一段(T₁j¹)及中统雷口坡组(T₂l)。研究区已成功施工地热井3口，分别为ZK1、ZK1-1、ZK3(表 1)，其钻井热储层的柱状对比图见图 2，各井的主要出水段分别为T₁j²⁺³⁺⁴、T₁j²、T₁j²⁺³，其岩性主要为灰岩、白云质灰岩及白云岩，据钻井过程中水文地质观测，由于灰岩、白云岩中岩溶管道发育具有极大的空间异质性，致使相邻区域统一地层的钻井涌水量差异较大。完井后各井分别进行了抽(放)水试验，其降深及对应的出水量差异极大(表 1)，以此为基础计算出的渗透系数差别也较大。

2 概念模型的建立 2.1 水文地质概念模型

水文地质概念模型是对水文地质条件的简化，是对地下水系统的科学概化^[8]，研究区地热水热储层为嘉陵江组(T₁j)和雷口坡组(T₂l)，也是本次模拟研究的目的层，整体上该含水层可概化成非均质-各向异性的潜水-承压水含水层。

研究一个目标系统对地下水的影响时，当周边没有第一类边界条件或者距离较远时，可以以目标系统的边界为基础往四周延伸3~5 km作为模拟区的边界，概化为二类流量边界^[9-10]。因此，模型的南、北两侧作为二类流量边界进行处理。含水层上部J+T₃xj地层和下部的T₁f地层，作为隔水层处理。研究区边界概化如图 3所示。

2.2 数学模型的建立

根据水文地质概念模型及含水岩组的水力性质，可将模拟区地下水流概化成非均质各向同性、非稳定地下水流系统，并建立相应的数学模型：

$ \left\{ \begin{array}{l} \frac{\partial }{{\partial x}}\left( {{k_x}\frac{{\partial h}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{k_y}\frac{{\partial h}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {{k_z}\frac{{\partial h}}{{\partial z}}} \right) + W = {u_{\rm{s}}}\frac{{\partial h}}{{\partial t}}, \\ h\left( {x, y, t} \right){|_{t = 0}} = {h_0}\left( {x, y} \right), \left( {x, y} \right) \in D, \\ h\left( {x, y, z, t} \right)|{\mathit{\Gamma} _1} = h\left( {x, y, z, t} \right), \left( {x, y, z} \right) \in {\mathit{\Gamma} _1}, t \ge 0, \\ k\frac{{\partial h}}{{\partial n}}|{\mathit{\Gamma} _2} = q\left( {x, y, z, t} \right), \left( {x, y, z} \right) \in {\mathit{\Gamma} _2}, t > 0。\end{array} \right. $

式中：k_x、k_y、k_z分别为渗透系数在x、y、z方向的分量，m/d；h为地下水水位，m；W为单位体积流量，用以表示流进源或流出汇的水量；u_s为含水岩组的储水率，L/m；h₀(x，y，z)为已知水位分布，m; t为时间，h；D为模拟区范围；Γ₁为一类边界；Γ₂为二类边界；n为边界上的外法线方向；q为二类边界上的已知流量分布。

2.3 网格的剖分

根据Visual Modflow地下水流数值模型要求，将模型3层采用网格剖分，并对ZK1、ZK1-1、ZK3地热井所在位置进行网格细分，平面上将计算区域剖分成为75行29列，共2 175个网格。

2.4 参数的初值

利用ZK1井2010年1月19日~22日、ZK1-1井2012年11月25日~28日、ZK3井2016年12月17日~20日抽(放)水试验获取水文地质参数初值(表 1)。3口井渗透系数差异很大，分析其直接原因是深部白云岩、灰岩中岩溶管道、岩溶裂隙等(岩溶发育程度)发育极不均一，即具有极大的空间异质性。但受目前勘查技术手段及深井钻探成本的限制，无法详细了解整个研究区(3个钻井以外)的深部岩溶发育情况，故采用插值法确定区域渗透系数。插值后生成渗透系数等值线作为参数初值，目标含水层的渗透系数为0.01~0.82 m/d(图 4)。

2.5 模型的识别与验证

模型识别与检验的目的是检验所建数学模型在识别与验证期受各种因素激励后地下水动力场与实测值时空分布的一致性，分析数学模型对地下水系统的行为和功能的适应性，深化对水文地质条件的认识，确定水文地质参数^[11-12]。

利用ZK1井2010年12月、ZK1-1井2013年7月、ZK3井2017年5月抽(放)水试验资料对模型进行识别和验证。通过调整水文地质参数和边界条件，使计算的水头值与实测的水头值之差最小，以取得最佳的拟合效果。ZK1、ZK1-1水头观测值与计算值的偏离误差见图 5~6及表 2~3。

通过数据分析，ZK1观测值与计算值的残差均值为1.08 m，ZK1-1观测值与计算值的残差均值为2.31 m，说明模型的拟合程度较好。

通过反复调整参数，识别了研究区的水文地质条件，确定了模型结构、参数，在此基础上，对研究区含水层的模拟水头值与初始水头值进行对比分析(图 7)。通过模型识别，模型的水头动态曲线与实测的水头曲线达到了较好的拟合，两者的动态变化过程也比较吻合，模拟的流场基本反映了地下水流的现状特征，流场拟合情况较好。其成果可以用于在不同干扰条件下地下水流场预测。

2.6 识别的结果

通过模型的识别和验证，获得了含水层的水文地质参数，水平渗透系数K_xx=K_yy=0.011 0 m/d，垂直渗透系数K_zz=0.001 1 m/d，每年降雨补给量500 mm。

3 模型预测

地热井井间距和开采量是影响井间干扰的重要因素，笔者拟通过改变地热井的井间距和开采量，预测不同条件下地热水流场的变化。本次假定在区内增加一个钻孔，其开孔层位均为侏罗系中统上沙溪庙组(J₂s²)，终孔层位为三叠系下统嘉陵江组第一段(T₁j¹)，当该井开采时，其余现有钻井均不开采。拟定了3种开采方案，对假定新增的钻井设置不同的井间距与开采量(表 4)，预测不同方案下区内地下水流场变化，从而判定新增井是否对现有ZK1井产生干扰。

按照拟定的开采方案，运用识别和验证的地下水流数值模型进行模拟，预测了拟定开采井对已有井ZK1的影响(图 8)。

由图 8(a)(b)可知，当距离ZK1井1 km的KC1井开采量为500 m³/d时，约8.0 h后以KC1井为中心形成降落漏斗，水头约下降76 m，ZK1井水头未发生变化，始终维持在+453 m左右，对ZK1附近流场无影响；假定开采量为1 000 m³/d，约14.0 h后以KC1井为中心会形成降落漏斗，水头约下降154 m，ZK1井地下水头也发生变化，即以ZK1井为中心形成降落漏斗，约1.6 h后水头约下降19 m，对ZK1井产生干扰。

由图 8(c)(d)可知，当距离ZK1井2 km的KC2井假定开采量为1 000 m³/d时，约10.0 h后以KC1井为中心形成降落漏斗，水头约下降115 m，ZK1井水头未发生变化，始终维持在+453 m左右，对ZK1附近流场无影响；假定开采量为1 500 m³/d时，约20.0 h后以KC2井为中心会形成降落漏斗，水头约下降154 m，ZK1号井水头也发生变化，即以ZK1井为中心形成降落漏斗，约1.6 h后水头约下降19 m，对ZK1井产生干扰。

由图 8(e)~(h)可知，当距离ZK1井5 km的KC3井假定开采量为1 500，2 000，2 500 m³/d时，分别在约25.0，28.0，30.0 h后以KC3井为中心形成降落漏斗，水头分别会下降58，111，165 m，ZK1井水头均未发生变化，对ZK1附近流场均无影响；假定开采量为3 000 m³/d时，以KC3井为中心会形成降落漏斗，约31.0 h后水头约下降111 m，ZK1井水头也发生变化，即以ZK1井为中心形成降落漏斗，约3.5 h后水头约下降8 m，对ZK1井产生干扰。

4 结论与展望 4.1 结论

应用Visual Modflow软件建立了观音峡背斜静观片区的地热水水流数值模型，利用后期抽水试验资料对模型进行识别和验证，确定了静观片区地热水模型结构和水文地质参数，在此基础上预测了不同条件下的井间干扰程度：

1) 在距已有地热井ZK1井1 km处布置钻孔，按500 m³/d开采，对ZK1附近流场无影响；按1 000 m³/d开采，以ZK1井为中心会形成降落漏斗，约1.6 h后水头约下降19 m。

2) 在距已有地热井ZK1井2 km处布置钻孔，按1 000 m³/d开采，对ZK1附近流场无影响；按1 500 m³/d开采，以ZK1井为中心会形成降落漏斗，约2.5 h后水头约下降13 m。

3) 在距已有地热井ZK1井5 km处布置钻孔，按1 500，2 000，2 500 m³/d开采，对ZK1附近流场无影响；按3 000 m³/d开采，以ZK1井为中心会形成降落漏斗，约3.5 h后水头约下降8 m。

模拟结果表明，地热井井间距越小、开采量越大，对相邻井产生干扰的可能性越大，反之越小，这与水文地质规律也是相符的，说明利用Visual Modflow模型模拟地热水流场特征，预测地热井干扰程度这一方法是可行的。

4.2 展望

目前国内对深层地热水数值模拟研究较少，笔者尝试利用Visual Modflow软件建立该片区的地下水流数值模拟模型，预测分析不同开采方案下的井间干扰程度，得到了预期的结论，为地热资源勘查合理布井提供了方法与理论参考。

地热井井间干扰也受到岩溶含水介质特征(岩溶管道的连通性)、补径排条件、抽水时间等因素的影响，笔者未考虑这些因素，这也是下一步研究的目标。另外，由于岩溶含水介质的极大空间异质性，在其他地区建模过程中，需尽可能收集地热钻探编录、物探、水文地质观测等资料，获取尽可能多的水文地质参数，提高模型预测的准确性。