初道忠，赵忠琦,2

(1.山东理工大学 资源与环境工程学院，山东 淄博 255049；2.山东黄金集团玲珑金矿，山东 烟台 265400)

随着信息技术的飞速发展，利用数值法计算矿井涌水量的方法得到了广泛的应用。这类方法能在考虑含水层的各向异性的同时，处理水文地质边界条件较复杂的矿区，模拟矿山开采过程中矿井疏排水情况及涌水量变化特征。

地下水的数值模拟计算方法主要有有限单元法、有限差分法、边界元法以及有限分析法等。本文利用有限差分法进行数值模拟计算，其基本思想为把连续定解区域用有限个离散点构成网格代替，把连续定解区域上的连续变量函数用离散变量函数近似表示。本文利用GMS软件中的BOREHOLE模块及MODFLOW模块对大庄子矿区的矿井涌水量进行预测，并利用模拟结果对矿山开采进行指导。

1 矿区基本水文地质条件简介

矿区分布于华北板块(Ⅰ级)、胶辽隆起区(Ⅱ级)、胶北隆起(Ⅲ级)、胶北断隆(Ⅳ级)、明村-但山凸起(Ⅴ级)处。该区域属河流冲洪积堆积地貌类型，大多为平原，地势起伏情况较缓，总体呈现东南高西北低的地势特点，最低侵蚀基准面标高位于矿体以上，为+10 m左右。

矿区属于北温带半湿润季风气候，根据该区域气象局1992～2014年气象资料，该区域多年平均气温12～12.2 ℃，1998年8月13日达到高温最值，为38.6 ℃，1992年1月19日达到低温最值，为-18.3 ℃。区域降水情况如图1所示，年平均降水量680 mm，年平均蒸发量1 711.1 mm，蒸发量大于降水量,年平均相对湿度69%。矿区地表排水设施完善，可轻松排走由于降水所存积的积水。该区域距离大规模水体较远，区域西部有一条流入莱州湾的胶莱河，区域南部有一条汇入胶莱河的淄阳河，区域北部有一条流入莱州湾的沙河，三者都为间歇性河流，旱季大都干涸，仅雨季具有较大流量，不会影响到目标区域矿床的开采。

图1 1981—2010年区域降水量柱形图Fig.1 Column chart of regional precipitation from 1981 to 2010

依据地下水的赋存条件、水理性质、水力特征等，将矿区地下水大致划分为三大类型：松散岩类孔隙水、碎屑岩类孔隙裂隙水和基岩裂隙水。各层的补给条件主要为大气降水、基岩风化裂隙水以及相邻的裂隙含水层的地下水静储量等。通过对钻孔终孔地下水位埋深的研究发现，在枯水期时各钻孔水位埋深大多低于第四系底板，而丰水期时水位埋深则超过第四系底板。故目标区域地下水不承压，为潜水。

2 建立矿区水文地质模型

2.1 含水层概化

根据矿区含水层的主要分布特点，对复杂含水层分布情况进行概化[1-2]，在不失真的前提下，将矿区含水层概化为四层，由上到下分别为第四系松散孔隙含水层、林家庄组砂岩含水层及荆山群基岩裂隙含水层联合含水层、云山单元二长花岗岩隔水层以及荆山群基岩裂隙弱含水层。具体概化分层方法如图2所示。

图2 含水层概念模型Fig 2 Aquifer conceptual model

在四层概化地层之下为另一隔水层，由于矿体赋存于荆山群基岩裂隙弱含水层内，故分析四层概化地层之下的隔水层无意义。其中第四系松散孔隙含水层、林家庄组砂岩含水层受季节影响较大，降水补给量占比较多，地下水为潜水；由于有云山单元二长花岗岩隔水层遮挡，荆山群基岩裂隙含水层地下水具有微承压性，但由于其基岩裂隙不发育，含水量较少，故承压性质不明显[3-4]。

2.2 确定边界条件

矿区属鲁东低山丘陵水文地质区胶莱拗陷水文地质亚区胶莱盆地孔隙裂隙弱—强富水地段，地处滨海平原区，地势比较平缓，最大比高10 m左右。矿区内最低侵蚀基准面标高为+10 m。南2 km有淄阳河自东向西流入泽河，为间歇性河流，大部分河段常年干涸，故地表水与矿区地下水水力联系不密切。

根据矿区的水位等值线图(图3)及矿区所在区域方位，矿区西部边缘主要为裸露的云山单元二长花岗岩，因二长花岗岩裂隙不发育，隔水性良好，故以其为矿区西部第二类隔水边界；矿区南北部边缘以其等水位线为边界，由于水头已知，故为第一类边界；由于矿区东部为另一矿区的勘探区域，地质资料不详，故以矿区东部边界为界，根据原始流场模拟情况确定为第一类边界[5]。

图3 矿区水位等值线Fig.3 Water level contour of mining area

2.3 建立渗流数学模型

根据矿区的含水层情况，结合地下水流补给及排泄情况，评估区视为三维非稳定流系统[6]，利用三维非稳定流微分方程建立数学模型如下所示：

(1)

H(x,y,z,t)E,S,N=φ(x,y,z,t)

(2)

(3)

H(x,y,z,t)|t=0=H0(x,y,z)

(4)

其中:K为含水层渗透系数；H为评估区水头分布情况；Ss为含水层贮水率；H0为t=0时的水头分布情况；W为单位时间单位体积含水层补给或排泄水量；(x,y,z)为空间坐标，E,W,S,N代表评估区东西南北边界；t代表时间。

渗流数学模型是地下水模拟的理论指导，是模型的概念基础。该数学模型中，公式(1)表示评估区域各含水层采用潜水三维均质非稳定流微分方程计算方法进行模拟；公式(2)、(3)表示评估区域各个方位的边界条件情况；公式(4)表示评估区域各位置的初始水头情况。

2.4 确定各含水层渗透系数

根据已经建立的渗流模型，对各含水层赋予初始渗透系数，按照抽水试验中抽水井及观测井的位置，在渗流模型中同样布置观测井，不断修改各含水层的渗透系数，对渗流模型进行反演迭代计算，使模拟结果与抽水试验实测结果相一致，进而确定各含水层渗透系数。

利用观测孔JG07、JG11、JG24在2016年2月19日～2016年2月27日时间段内的抽水试验水位数据，对各个测量水位的时间节点进行编号，结合渗流模型模拟结果，其拟合曲线如图4所示。

(a)JG07观测孔水位拟合曲线

求得含水层1、含水层2、含水层3、含水层4的渗透系数分别为：3.54 m·d-1、0.81 m·d-1、0.000 52 m·d-1、0.13 m·d-1。

3 抽水数值模拟

3.1 原始流场分析

以地下水模拟软件GMS为平台，利用已知的钻孔数据结合borehole模块[7]生成钻孔模型以及概化地层模型如图5、6所示。

利用GMS的MODFLOW模块[8-9]对含水层模型赋以边界条件及降水条件见表1。

根据原始流场情况可知，流场与边界水位反映情况一致，但由于降水缘故，流场中潜水水位面高于边界水位面，说明降水对于含水层1、2层影响较大，而含水层3、4受到降水影响较小，主要受边界径流以及泄流影响。

模拟非抽水条件下原始流场如图7所示。

3.2 抽水数值模拟

在原始流场的基础上，设立假想抽水井，不断调试假想抽水井的日抽水量，观察流场的变化情况，进而得到最优日抽水量。假想抽水井的最优日抽水量见表2。根据评估区域的气象资料可知，评估区域降水量最大的时间段为六到八月份，故以该时间段评估矿井的最大涌水量。

表2 日抽水量汇总表Tab.2 Daily pumping summary table

根据最优日抽水量模拟流场情况如图8、9所示。

图8 抽水数值模拟结果Fig.8 Pumping numerical simulation results

根据抽水模拟结果，红色区域代表地下水疏干区域，为四个抽水井形成的水位降落漏斗范围，周围彩色单元表示由于边界水头条件及降水条件仍然存留的地下水含水单元。由于矿体赋存于荆山群野头组基岩裂隙弱含水层的围岩蚀变带中，故以第四含水层的疏水情况作为排水标准。

(a) 含水层1模拟结果

由以上红色区域分布情况可知，红色区域分布于四个含水层中的每一层，且至少占每层面积的80%以上。由此可知地下水日抽水量为1 715～2 498m3/d，疏水时间为六月到九月。

对比抽水数值模拟结果与该矿区地质详查报告中利用大井法[10]预测抽水结构可知，大井法预测最大抽水量为1 532 m3/d，与数值模拟结果基本吻合。两种方法存在偏差主要是由于大井法计算过程较为简单，故精确度较低造成的。

4 结束语

通过对大庄子区域矿井涌水量的数值模拟过程可得出以下结论：(1)该矿区的地下水日疏水量为1 715～2 498 m3/d，矿井最大涌水量为2 498 m3/d；(2)利用GMS对地下水进行数值模拟的方法需要考虑地层结构、地下水的径流影响、降水情况、边界条件以及渗流模型的选取等问题，适用于地质资料较为详细时使用，故其计算结果更为可信；(3)数值模拟过程中，如果地层情况较为复杂，概化地层难以表达，可将复杂地层详细划分后分别进行概化。