基于大井法和地下水模型系统数值模拟方法的某矿坑涌水量预测对比分析

孙刚友, 胡清珍, 康钦容, 夏缘帝, 袁威, 张卫中*

(1.中核第四研究设计工程有限公司, 石家庄 050021; 2.武汉工程大学资源与安全工程学院, 武汉 430073)

在矿山开采的过程中,矿井突涌水这一灾害时常发生,不仅阻碍了资源的开发,而且对人员的安全也存在隐患[1]。如若在矿坑建设生产前期,通过合理科学的方法对矿坑涌水量进行预测[2-3],将极大程度地服务于矿坑后期各阶段的水文地质调查工作的开展,对矿山在实际开采时安全地开展工作,以及防止灾害事故发生,都有极大的作用[4]。

矿坑涌水量的预测方法有数值模拟法[5-7]、水文地质比拟法[8-9]、水均衡法、解析法[10-12]。近年来也研究出了很多相关分析方法,施龙青等[13]对不同时间下的涌水量进行统计,分析两者间的函数关系;李占利等[14]提出互补集合经验模态分解与门控循环单元相结合的矿井涌水量预测模型;姜小妮等[15]采用非完整大井法计算矿坑涌水量,提高了计算精度;李建林等[16]用称重标极差分析(rescaled range analysis,R/S)法分析涌水量序列的周期变化趋势并采用灰色预测模型对矿井涌水量进行了预测,提供了一种新思路。但这些方法也存在不足之处,水文地质比拟法适用于开采同一矿井或者地质环境基本相同的区域,每个矿区都有其不同的水文地质条件,所以此方法的使用具有局限性;和时间序列有关的分析法需要建立在观测数据之上,观测时间久,数据收集困难,结果易受影响;大井法[17-18]因其简捷的计算得到了广泛的应用,但矿山实际上并没有完全满足公式的前提条件,因而其预测结果将会与矿山实际涌水量偏差较大,单一方法计算涌水量的说服力不强,准确性不够。

因此现以矿坑涌水量预测精确度[19-20]为主,研究调查该矿坑丰富的水文地质数据,在采用大井法预测后又使用适应性强、预测精度高的地下水模型系统(groundwater modeling system,GMS)水文地质软件对矿坑涌水量进行计算预测。使用大井法进行预测由于实际环境条件不完全匹配公式的要求,会出现误差导致精确度不够,采用GMS数值模拟软件加以检验,两种方法的相互补充验证使该矿坑在预测涌水量的精确范围上得到提升,为矿山开采提供有力的技术支撑。

1 研究区概况

1.1 区域地质及水文地质条件

矿区位于内蒙古自治区阿拉善右旗东北部某地段,地势较为平坦,地表大多被第四系沙土和沙丘覆盖,区内水系多为季节性水流和干河床,植被稀少,断裂构造以压扭型断裂为主。

1.2 矿床水文地质特征

1.2.1 主要含、隔水层及富水性

据钻探结果显示该矿床主要的含水组分为3层:下层埋藏较深,以泥岩和粉砂岩为主,泥岩稳定性较好,具有分隔作用,而且该层位于矿床下部,因此作为隔水底板;中层含矿量较多,无明显挤压破碎带发育,岩心块度相对完整,通过岩心编录对6个钻孔的节理和裂隙进行了详细研究,发现所有水文地质孔的过滤器段裂隙发育程度低甚至不发育,开启裂隙占裂隙总数比重较低,是主要的含水层;上层由含泥砂岩组成,不易透水的特质使其作为隔水层。

1.2.2 地下水补径排条件

由于泥岩的厚度大且分布稳定,作为隔水顶板可隔绝大气降水的渗入,含矿含水层地下水的主要补给是地下水的侧向径流补给,在矿床的西南部约10 km的地方,设置了一个孔隙为800 m的地质孔,在该孔中出现了冒水及气泡溢出的情况,说明在地下水径流过程中在矿床的西南部遇到阻隔,最终确定总体径流方向为由北东向南西径流。

2 大井法预测矿坑涌水量

2.1 计算公式的选择

在实际开采时,开采场地需达到合理安全的开采条件,地下水位的降低使在矿坑和远离矿坑的地下水部位的压强不一致。因此在进行矿坑涌水量预测时,需选用地下水动力学的大井法的承压转无压公式,公式如下。

(1)

R0=R+r0

(2)

(3)

(4)

式中:Q为矿坑涌水量预测值,m3/d;Kcp为平均渗透系数,m/d;K为渗透系数,m/d;S为地下水位降深,m;M为含水层平均厚度,m;R0为引用影响半径,m;r0为大井引用半径,m;R为影响半径,m;P为圈定的大井面积,m2。

2.2 参数确定

2.2.1 抽水试验取得的资料

统计整理了每个孔承压含水层的厚度,以及作为抽水孔、观测孔以及采用水位恢复法时的渗透系数,如表1所示。

表1 抽水试验成果表

2.2.2 渗透系数的确定

在预测矿坑涌水量时,因各段含水层渗透性分布不均,渗透系数的取值由抽水孔、观测孔以及水位恢复法这3种方法每个观测点的渗透系数值与含水层厚度的加权平均值确定,其值见表2。

表2 渗透系数的加权平均值计算表

2.2.3 矿床地下水平均水位标高的确定

10个水文地质孔一个水文年的地下水位标高观测值,如图1所示,取这10个水文地质孔地下水位标高观测值的平均数,得出矿床地下水平均水位标高为1 273.66 m。

图1 地下水位标高统计图

2.2.4 地下水位降深的确定

为保证矿床正常生产,在生产时地下水位必须降至矿体底板以下。对8个主要矿体埋深的最低标高和最高标高进行了统计,如图2所示,计算得出主要矿体埋深最低标高平均值为773.49 m,最高标高平均值为793.90 m。

图2 主要矿体埋深标高统计图

矿床在生产时要保持采场作业环境,因此需要地下水位降至矿体埋深底板以下5 m位置,即水位降深应在上述标高平均值的基础上再加5 m,故所得地下水位降深为地下水位标高与矿体底板标高之差,如表3所示。

表3 常规开采时地下水位降深值

2.2.5 含水层厚度的确定

在施工过程中,对含有水文地质孔的岩心进行了标记,在进行涌水量预测时,确定了各个水文地质孔的含水层厚度,取13个水文地质孔的算术平均值M为82.45 m。

2.2.6 大井引用的确定

(1)矿体面积的圈定。根据矿体分布形态来圈定矿体面积见图3,最终测得矿床疏干面积P为1 549 726 m2。

图3 矿床疏干范围图

(2)大井引用半径的确定。根据前述大井法引用半径的计算公式,求得大井引用半径r0为703 m。

2.2.7 影响半径及引用影响半径的确定

结合确定的渗透系数及降深,根据前面的公式分别确定最高矿体标高及最低矿体标高的影响半径及引用影响半径数值如表4所示。

表4 影响半径及引用影响半径表

2.2.8 矿坑涌水量计算

根据公式和确定的参数,对矿坑涌水量进行计算,计算结果如表5所示。

表5 矿坑涌水量预测表

从表5可以看出,3种方法求得的渗透系数计算的矿坑涌水量差别并不大,在同一降深深度内,差值分别为2 904 m3/d和2 802 m3/d,误差分别为5.9%和5.8%。最终由2个降深、3个渗透系数计算的最大降深平均涌水量为47 948 m3/d,最小降深平均涌水量为46 501 m3/d,从而得知由于承压水头高,在抽水孔进行抽水试验时,水头损失并不明显,同时也说明该矿床矿坑涌水量较大。

3 地下水数值模拟预测矿坑涌水量

抽水试验得到了很多的渗透系数值,在以上工作的基础上结合前人工作资料,对该矿床的水文地质条件进行分析,初步分析了地下水的补给、径流、排泄,对水文地质单元进行了详细划分,建立三维可视化地下水模型。

3.1 水文地质概念模型

3.1.1 研究区范围

以某矿区所在地为凹陷盆地作为一个完整的水文地质单元,以第四系松散沉积物与白垩系基岩为分界线,模型范围确定了一个水文地质单元,矿床位于模型范围的中部,避免边界条件在计算中对结果的影响。

3.1.2 水文地质条件概化

将研究区水文地质条件进行逐一分析,对水文地质条件中的含水层结构、地下水径流及动态特征、模型边界条件等进行确定。

(1)含水层结构。模拟的深度以研究区所得的钻孔深度为参考,结合DEM提取标高确定,矿区的开采深度约在地下500 m。含水层水平方向不再分层,含水层结构垂向上分为上下两层,第一层含水少,砂体含量少;第二层为主要含矿含水层,从模型角度分析主要含水层的渗透系数大于第一层,地下水垂向流动速度相比于第二层中的地下水流动速度差2个数量级。研究区含水层受构造影响,导水构造对含水层的影响很大,故在模型中将构造影响加入计算中。含水层结构在三维地下水模型中概化为两层,见图4,地下水的垂向运动速度比水平运动速度低2个数量级。

图4 含水层结构图

(2)地下水径流及动态特征。地下水径流方向为北东至南西方向流动,受水文地质单元地理条件影响,有一定的侧向补给来源。北东方向为补给来源,南西为排泄方向,由于构造影响,上下两层在构造发育位置存在一定的水力联系。

(3)模型边界条件。在侧向上,北部是宗乃山-沙拉扎山隆起区,为地下水的补给区,因此北部边界可作为流入边界;南部为巴音诺尔公隆起,东部由那仁哈拉凸起和毕级耳台凸起相隔成半封闭状,作为隔水边界或者弱透水边界;西部为巴丹吉林沙漠,可作为隔水边界,如图5所示。

图5 模型边界划分三维图

3.2 地下水流数值模拟

3.2.1 地下水的数学模型

承压含水层在进行矿井排水时,地下水会由承压水变为潜水,首先承压水可用式(5)控制方程描述为

(5)

式(5)中:Ss为贮水率,1/m;H为承压含水层的水位标高,m;Kx、Ky、Kz分别为水平和垂向渗透系数,m/d;ε为含水层的源汇项,1/d。

潜水可用式(6)控制方程描述为

(6)

式(6)中:S为贮水系数;h为承压含水层的水位标高,m。

3.2.2 模拟计算软件

1)模型网格剖分

本次研究区的范围为矿床勘探区域,根据GMS数值模拟软件要求对研究区进行矩形网格剖分,初步确定在平面上进行正方形网格剖分,矿区作为典型重点研究地段,可采用局部加密技术对模型进行局部加密,见图6。网格剖分以200×200尺度进行,剖分后对某矿床研究区进行了适当的网格加密,即对网格进行细化。

图6 三维水文地质模型网格剖分示意图

2)水文地质参数

用于该模型的水文地质参数主要是承压含水层的渗透系数(包括水平渗透系数和垂向渗透系数)和储水系数。水文地质参数的表示包括空间位置和参数值,空间位置可用Shape格式中的线或面给出,水平渗透系数和储水系数参照前人的工作成果以及本次工作钻孔的抽水试验计算得到初值以及矿区环境条件。建模时将得到的初值根据矿坑环境条件进行参数分区,然后通过水位拟合进行参数识别,参数调整时考虑到了构造影响,在构造发育位置对透水构造进行了详细刻画,见图7,最后对各分区的识别参数值进行确定。

图7 水文地质渗透系数分区图

3)源汇项确定

根据地下水补径排条件,确定地下水均衡要素,并采用相应计算公式,研究区承压含水层的补给量为侧向流入量,根据各层地下水位等值线,用达西断面法求各边界流入量初值,并通过模型进行校正。研究区承压含水层的地下水排泄量主要如下。①断层排泄:先根据研究区的水文地质条件估测一个初值,然后通过代入模型进行识别验证;②矿井抽水排泄:通过识别验证后的模型进行逆运算。

4)模型检测

调参过程中对于研究区的水文地质条件进行不断修正,达到水位线与观测值基本吻合,经过多次检测对比以后确定,通过图8最终可以看出模型经过10阶段调参后,模拟的东部等水位线水力梯度与观测值基本吻合,地下水模拟流场的西部位置流场与观测值亦基本吻合,地下水的运动趋势和观测值基本一致,误差较小。

图8 第10阶段调参后等水位线图

5)涌水量预测

对调参后基本吻合地下水观测值的水文地质模型进行降深运算,逐步得到涌水量数值,涌水量计算在模型中进行,抽水量通过逐步分析得到。涌水量计算数值通过数次计算得到,仅展示计算变化比较大的图件,图件数与实际计算次数比例约为1∶20。在模型调参过程中导出矿体模型,将矿体模型输入了GMS软件中,对矿体和水文地质模型关系进行分析,从而更形象地得到矿体与含水层的关系,如图9所示。

图9 矿体与水文地质三维模型关系图

如图9所示的矿体位置可以看出矿体的一部分处于含水层中,另一部分不在含水层中,属于非渗透性矿体,理清了矿体和含水层的关系后,进一步对涌水量数值进行确定。

经过反复调参,模型运行多次后,对模型运行和计算结果分析后得到了矿床开采时的涌水量数值,调参数过程中确定了涌水量数值范围在26 000～43 000 m3/d。在第11次进行涌水量数值32 000 m3/d输入后,含矿含水层的等水位线疏干形态几乎覆盖了整个模型的东部区域,此时在三维图(图10)中可以看出在矿床位置第一层含水层和第二层含矿含水层的疏干范围都覆盖了矿床,第二层满足矿床的范围处于完全疏干的状态,由于构造影响,第二层构造位置也有一定的疏干,符合地下水与构造有水力联系的实际情况,确定矿床开采时的涌水量为32 000 m3/d。

图10 涌水量计算第11次输入后含矿含水层底板等水位线三维图

4 结果对比评价

通过大井法预测的矿坑最大降深平均涌水量为47 948 m3/d,最小降深平均涌水量为46 501 m3/d;GMS软件预测的平均矿坑涌水量为32 000 m3/d,最大涌水量为38 000 m3/d。两种方法预测的矿坑涌水量有所差异,大井法预测公式是在含水层水平、等厚、均质等严格条件下推导出的,矿山实际上很难完全满足公式的前提条件,因而其预测结果将会与偏大,但也提高了一定范围的参考。随着计算机的发展,数值模拟方法适应性强、预测精度高,已经被大量应用于各个方面,而且此次对该矿坑水文地质数据的收集较为全面,因此由数值模拟得出的参考范围将会更精确。两种预测方法相结合,互相补充验证,为实际的矿坑涌水量提供了有力的参考。

5 结论

(1)矿坑涌水量随矿坑开拓面积的扩大而增加,单一预测方法所建立的模型存在预测精度不高、误差大的问题,且其适应范围窄, 大井法预测公式是在严格条件下推导出来的,其预测结果偏大;数值模拟方法预测精度高,由该方法得出的数据将会比大井法得出的数据范围小且更为精确。

(2)在计算矿坑涌水量的过程中,大井法和GMS数值模拟法两种预测方法相互补充验证,得出该矿坑涌水量的大范围和小范围预测值,计算得到的矿坑涌水量更加可靠,可为实际开采矿坑时涌水量数值提供精确范围的参考。通过双重预测方法的结果表明该矿床矿坑涌水量较大,研究成果将为开采矿山提供技术支撑。