西北干旱-半干旱地区，特别是晋陕蒙宁甘地区，富集了全国65%的煤炭资源，却仅占有全国3.9%的水资源量[1-2]。水资源匮乏，已经并将持续成为能源“金三角”地区可持续开发道路上遇到的主要限制因子。
与此同时，西北富煤地区的煤炭开采实践表明，地下赋存着丰富的地下水。采矿过程扰动下，地下水系统遭到破坏，形成大量矿井水。传统上这些矿井水作为煤矿水害外排到地表，造成水资源浪费。矿井水资源的浪费与该地区水资源短缺形成了一对矛盾。此外，水利部于2013年提出了“大型煤电基地开发规划水资源论证的意见”，明确提出“煤电基地燃煤电厂建设与煤矿开采等项目用水应统筹安排、综合利用”，使得矿井水外排成为煤矿的负担，大大限制了煤炭资源开采。煤矿地下水库建设的理念和实施使得以上问题得到高效解决。煤炭开采后形成了具有很大空隙的采空区，同时，通过地下煤柱坝体和人工混凝土坝体将矿井水安全的储存在采空区，形成煤矿地下水库。为了保证采矿安全，防止采空区持续涌水导致煤矿地下水库失稳，提出了建设多个地下水库，并通过相应通道连接，使不同储水的采空区可以形成相互之间贯通的地下储水空间。根据煤炭开采生产的安排计划，可以将矿井涌水分时分地存储在这些地下采空区的空间内，通过地面与矿井中建设的各种辅助工程措施来实现矿井涌水的过滤净化、抽提使用与回灌贮存，形成有效的地下矿井涌水的暂存空间，也将其称之为“煤矿分布式地下水库”[3-6]。
以某矿区为例，其正常日涌水量为450 m3/h左右。一旦将这些矿井涌水外排到地面，水会很快蒸发消失殆尽，造成水资源的浪费。笔者以该矿区分布式地下水库建设为例，主要研究:① 根据矿井开采区地质资料，计算不同高程区域储水能力，确定分布式地下水库的库容-水位曲线，为水库安全评价提供依据;② 建立分布式地下水库的地下水分析模型，预测矿区的涌水量，为分布式地下水库系统中每个水库的实际库容评估提供依据;③ 针对矿区的涌水量，结合地表水资源，按照该地区生产和生活需水量，进行矿区水资源的优化调度，使得矿区涌水得到合理利用。以上研究内容为煤矿地下水库的储水性能(库容评价和水量预测)及水资源高效利用提供基础理论和技术支撑。
1 矿区概况
该矿区地处陕西省榆林市神木县西北角，气候干旱缺少降水，地表水蒸发强烈，夏秋季节降雨集中。该矿区水系，西部属于乌兰木伦河水系，东部属于悖牛川水系。每年6—9月为丰水期，7，8月份雨量集中，经常出现山洪;10月至次年5月为枯水期。矿井地层共有5段延安组，共含煤20余层。2-2及5-2煤层，全层范围内均可进行煤炭开采，为主要可采煤层。
矿区东侧和北侧，发育有很多支沟。北侧基岩，呈现裸露状;该区域内潜水，形成了一个比较完整、可以局部完成水循环过程的水文地质单元。矿区地层的含水层和隔水层，基本呈交替式分布。含水层，主要有第四系河谷冲积层(Q4al)潜水、上更新统萨拉乌苏组(Q3s)潜水，下更新统三门组(Q1s)砂砾层含水层，中生界的中侏罗统直罗组(J2z)裂隙潜水，中下侏罗统延安组(J1-2y)裂隙承压水。隔水层，主要有第四系的中更新统离石黄土(Q2l)隔水层、中下侏罗统延安组(J1-2y)隔水层。如图1所示。
图1 某矿区沉积-构造剖面
Fig.1 Geology and structure profile of a mine
2 库容-水位关系曲线
煤矿地下水库，主要依靠垮落带、导水裂隙带和弯曲带的空隙、孔隙、裂隙等空间来储水，笔者将这些可以储水的空间体积称为地下水库的库容。库容的大小，直接决定了地下水库储水性能，因此，研究地下水库库容-水位之间关系，对于水库储水量规划和储水安全具有重要意义[7]。
2.1 分布式地下水库简介
目前，该矿区分布式地下水库主要建设在2-2煤层。已建设的水库有3座，依次为四盘区、老六盘区、新六盘区。本文中，将四盘区、老六盘区、新六盘区依从左到有编号为1～3号水库。5-2煤层三盘区已经完成采煤作业，正在建设4号水库(图2)。通过水库之间的管道进行水量调配，2-2煤层的3个水库与5-2煤层的1个水库实现了联网运行，这为将来矿井水在地下采空区之间的统筹优化调度奠定了重要基础。大量的采空区为水库提供了足够的储水空间，经过测算大约可储水约230万m3[8]。
图2 矿区地下水库分布
Fig.2 Distribution of the underground reservoirs in the mining area
2.2 库容计算
地下水库的储水系数，表征垮落带内的孔隙、空隙和裂隙带及弯曲下沉带内的裂隙占相应空间的比例。通过孔隙介质来进行储水的地下水库，其库容为
V=μV0
(1)
式中，V为地下水库库容(储水量)，m3;V0为开采扰动区总体积，m3;μ为煤矿地下水库储水系数。
矿区地下水库储水系数主要取决于采空区岩体空隙率，由采后时间、岩层性质、开采工艺参数、冒落岩体块度及堆积方式、矿山压力等因素决定[8-13]。根据现场的实测数据和模拟研究，矿区地下水库μ的取值范围在0.25～0.35[8]。
将开采底板作为水库底部(库容计算0高程点)，通过逐层累加的方法计算地下水库库容，确定了库容-水位曲线。矿区实测水位高程7 m时的储水能力约为230万m3，计算结果为235.18万m3。其5～10 m高程的总库容如图3所示。
图3 分布式地下水库5～10 m的库容-水位关系曲线
Fig.3 Relationship between reservoir capacity and water level for 5-10 m underground reservoir in mining area
3 三维地下水数值模拟
对煤矿地下水库区域地下水进行数值模拟是了解水库建设对地下渗流场影响的最重要方法之一[14-15]。本文基于矿区实际工程地质水文地质条件及采动影响后采空区分布，建立了三维地下水数值模型。地下水在三维各向异性空间流动的控制方程为
(2)
式中，h为地下水水头，m;kxx，kyy，kzz为地下水主渗透系数，m/s;ss为释水率，m-1;h1为水头边界条件，q为流量边界条件，h0为初始水头分布。
3.1 水文地质概念模型
矿区模拟区域范围，包含矿井边界内的所有区域(图4)。x轴为东偏北30°方向，与运输大巷一致;y轴为北偏西30°方向;z轴为竖直方向，数值与海拔高程一致。模型x向的长度为14 km，y向长度为11 km，总面积154 km2;z方向从海拔高程500 m的地层至最高约1 320 m的地表。基于各地层分布、采空区底板高程信息及水文地质条件，确定河流、水库、泉、井等补给区、排泄区的位置。图中包含的要素有:① 4座地下水库;② 泉和水井;③ 乌兰木伦河和悖牛川两条主要水系。
图4 渗流模拟补排水边界
Fig.4 Boundary of water charge and discharge
根据区域内各岩层的特性，将其划分为煤层、含水层、相对隔水层。煤层，重点放在2-2和5-2煤层，并规划各盘区作为分布式地下水库。考虑因开采扰动后渗透系数的变化，模型共分为36层。其中32～36层为拓展地层，主要目的是为减小下边界对5-2煤模拟结果的影响。
由于缺少初始渗流场资料，在模拟实际工况之前，将各个计算单元体的底部高程作为初始渗流场的水头值。在稳定流计算中，设定与计算结果相近的水头初始值用以提高计算效率。首先，将模型边界设定为隔水边界。根据区域内乌兰木伦河与悖牛川水位，提取边界网格水头赋予模型定水头边界，以减小人工设定边界对模拟区域的影响，并依据此边界条件按照实际工况进行模拟计算。
3.2 数学建模
根据以上边界构建3D计算网格，并将各地层空间坐标特征点导入计算模型，建立三维模型(图5)。
图5 矿区三维模型
Fig.5 3D groundwater simulation model for the mining area
本次模拟的水量补给来源，包括河流和降水。河流的补给，通过设定定水头边界来实现。降水的补给，通过设定日均降雨来实现。根据水文资料，该地区年平均降水量为450 mm，取渗透补给系数为0.333，则日降水入渗补给量为0.41 mm/d。蒸发是水量排泄的主要去处，根据水文资料蒸发量为降水量的4～5倍，初步拟定日蒸发量为降水量的4.5倍，即1.85 mm/d，另取蒸发深度1.5 m。
模拟过程中水量排泄区(边界)，还包括地下水库、泉、水井。地下水库边界，通过给定地下水位(高程)和水力传导系数来实现;泉边界，通过给定日均涌水量来实现;水井边界，通过给定地下水位(高程)来实现。河流是一个特殊的源汇项。河流的水头高程，比2-2煤层的高程低，又比5-2煤层的高程高。河流，既接受降水形成地表径流及较高高程处的地下水入渗补给，又会给5-2煤层的4号水库进行河流补给。
3.3 计算结果分析
初取河流的水力传导系数为2 m2/d、地下水库区域的水力传导系数为0.001 m2/d进行计算，补给量与排泄量的误差为4.69%。通过现场试验测量数据进行反演，取河流的水力传导系数为2.48 m2/d、地下水库区域的水力传导系数为0.001 05 m2/d，使误差降至3.15%。
根据区域气象站多年月平均降水量数据，输入不同降水量和相应的蒸发量，对每月的涌水量进行计算，得出模型矿区的补给量与排泄量、以及各水库涌水量，如图6所示。
对结果进行分析，可以发现:① 河流向4号水库的补给流量很稳定;② 矿区涌水量保持稳定;③ 蒸发排泄量，随着降水量的增大呈线性增大趋势，变化幅度较大;④ 降水补给量、河流排泄量，其变化幅度基本与降水量一致;⑤ 通过与实测涌水量进行比较，当降水量较小时，计算出的涌水量值偏小。
3.4 涌水量修正
GM模型(灰色系统模型)，可借助离散随机数通过生成变为随机性已经得到明显弱化并且比较有规律的一系列生成数，建立微分方程形式的数学模型，GM模型方便对事物变化的过程进行探讨[16-17]。
根据已有的2008—2012年各个月份涌水量，建立GM模型，在此基础上对2013—2020年涌水量进行预测，并与GMS模型的预测结果对比，如图7所示。
从图7可以看出:① 降水量不是特别小时，两模型的预测结果拟合度较好，但GM模型预测偏大;② 当降水过小时，GMS模型随降水减少而迅速减少，这与实际的涌水量不相符。
因为GM模型预测结果偏大，故3—10月份的涌水量选取GMS模型预测值;11月至次年2月，GMS模型预测的涌水量明显偏小，取GM的预测结果作为该时段的涌水量。上述处理方法可使数值模拟矿区涌水总量结果更加准确，为以下煤矿分布式地下水库水量预测和优化调度提供依据。
图6 模拟矿区补给量、排泄量以及各水库日均涌水量 年际变化曲线
Fig.6 Simulated annual variation curves of mine supply, excretion and daily inrush water
图7 GM模型与GMS模型的结果值对比
Fig.7 Comparation between the GM results and GMS (Groundwater Modelling System) results
4 分布式地下水库优化调度设计
矿区水资源需水量，主要包括以下几部分组成:生产用水、生活用水、工业用水、生态绿化用水，见表1。
表1 各种用水量供水优先顺序说明
Table 1 Water supply of different kinds of usage
种类生产用水生活用水工业用水绿化用水绿化生态用水———1工业复用水1—12生活水源—1——联网自来水2223
注:“—”表示相应水源不供应该用水需求;数字1，2，3表示优先顺序依次降低，只有当优先权大的供水水源没有富余时，才使用优先权小的供水水源。
4.1 矿区各用户需水优化调度
优化调度，需要有相对应的数学模型要素。其要素主要有决策的变量、由决策变量建立的目标函数、决策变量需要满足的约束条件[18-20]。根据这3者，即可列出约束方程。优化调度数学模型表达式为
(3)
式中，
为n个决策变量向量;
为目标函数向量集(对于多目标问题，m≥2;对于单目标规划问题，m=1);V-min为极小化向量目标函数;Gi(x)为约束函数。
分布式地下水库水资源的优化调度模型，首先是在保证地下水库结构安全的条件下，使得水资源供给保证率最高、供水成本最低。
(4)
(5)
式中，M为用户类型总数;N为与用户类型相对应的水资源类型总数;m为用户类型编号，m=1，2，3，4;n为水源类型编号，n=1，2，3;xmn为第m个水源向第n个用户的供水量;f1(x)为全区供水量;αn为第n种水源供单位水量的成本系数，根据供水源与需水方的距离相关;h1(x)为全区供水量费用。
主要有水量需求、各水源供水能力、决策变量特定等3方面的约束方程。在保证各类用水需求得到满足的前提下，最大供给量不超过用户需水量120%的上限范围，这样做是为了节约水资源、降低用水的成本。另外，为表示严格执行供水顺序，不能供用时的未知量设定为0，其余应不小于0。
根据各用户的用水量需求，按照约束方程和目标函数建立方程组，对方程组进行求解，得出2015—2020年的供水源向各用户供水方案，其中表2和3分别为2015和2020年的供水明细。
表2 2015年供水源向各用户供水明细
Table 2 Water supply volumes of different kinds of usage in 2015m3/d
种类生产生活工业绿化自来水2 0272 6264 102192工业复用水1 25702 883206绿化复用水0001 032
表3 2020年供水源向各用户供水明细
Table 3 Water supply volumes of different kinds of usage in 2020m3/d
种类生产生活工业绿化自来水2 4493 2824 151101工业复用水1 65403 080133绿化复用水0001 531
4.2 各水库联合调水方案
联合调度方案的目标，是使矿区各水库在丰水期时能维持地下水位尽可能上升得缓慢一些，以保证工程安全。而枯水期通过增加调用水量使地下水位下降到水文年初的水位。
各个水库均从7.0 m水位开始调试。调试目标为:经过1个水文年，水库水位又恢复到7.0 m。经过反复调试，确定1～4号水库的供水量分别为4 725.6，2 190.8，2 608.7，1 533.0 m3/d，能使水库水位基本保持稳定。
图8 联合调试的各水库月末水位
Fig.8 Water level estimate by combined scheduling method
4.3 优化结果综合分析
通过观察联合调试的各水库月末水位图(图8)，可以发现，矿区涌水高峰时水库水位变化也不明显，这说明水库的调蓄能力比较强。在降水段集中的雨季，水库可充分将多余的涌水暂存。地下水库多年日平均供水能力为11 058 m3，2015—2020年预测的生产、工业和绿化供水需求总量分别为11 698，11 987，12 249,12 632，12 854，13 097 m3。可见，仅地下水库的水量不足以维持矿区用水需求，还需通过联网自来水进行补充。
5 结 论
(1)提出了煤矿分布式地下水库的储水形式和库容-水位曲线确定方法。
(2)建立了分布式地下水库的渗流场模拟模型，分析了矿区涌水量。提出了GM(1,1)涌水量预测模型。
(3)根据矿区各种用水需求，编制程序实现了整个矿区范围内的优化调度。并根据矿区的涌水量，进行地下水库的联合调度调试，得出确保地下水库水位稳定条件下可以供应的水量。
[1] 王永光.中国煤炭市场供需平衡分析与对策研究[D].太原:山西财经大学，2007.
WANG Yongguang.The analyse and countermeasures study on the balance of coal market supply and demand in China[D].Taiyuan:Shanxi University of Finance & Economics,2007.
[2] 顾大钊.能源“金三角”煤炭开发水资源保护与利用[M].北京:科学出版社，2012.
[3] 顾大钊.煤矿分布式地下水库技术及其应用[A].国家能源局第二届能源论坛论文集[C].北京:煤炭工业出版社,2012:68-74.
GU Dazhao.The coal-mine distributed underground reservoir’s technology and its application[A].The second session of the National Energy Board Energy Forum[C].Beijing:Coal Industry Press,2012:68-74.
[4] 范立民.保水采煤面临的科学问题[J].煤炭学报,2019,44(3):667-674.
FAN Limin.Some scientific issues in water-preserved coal mining[J].Journal of China Coal Society,2019,44(3):667-674.
[5] 杜献杰,冯国瑞,戚庭野,等.结构充填“保水-储水”采煤顶板稳定性分析[J].煤炭学报,2019,44(3):820-829.
DU Xianjie,FENG Guorui,QI Tingye,et al.Roof stability analyses of“water-preserved and water-stored”coal mining with constructional backfill mining[J].Journal of China Coal Society,2019,44(3):820-829.
[6] 庞义辉，李全生，曹光明，等.煤矿地下水库储水空间构成分析及计算方法[J].煤炭学报,2019,44(2):557-566.
PANG Yihui,LI Quansheng,CAO Guangming,et al.Analysis and calculation method of underground reservoir water storage space composition[J].Journal of China Coal Society,2019,44(2):557-566.
[7] 李旺林,刘长余,汤怀义.地下水库设计理论与工程实践[M].郑州:黄河水利出版社，2012.
[8] 顾大钊.能源“金三角”煤炭现代开采水资源及地表生态保护技术[J].中国工程科学,2013，4:102-107.
GU Dazhao.Water resource and surface ecology protection technology of modern coal mining in China’s energy “Golden Triangle”[J].Chinese Engineering Science,2013,4:102-107.
[9] 黄庆享.浅埋煤层覆岩隔水性与保水开采分类[J].岩石力学与工程学报,2010,29(S2):3622-3627.
HUANG Qingxiang.Impermeablility of overburden rock shallow buriedcoal seam and classification of water conservation mining[J].Chinese Journal of Rock Mechanics and Enigeering,2010,29(S2):3622-3627.
[10] LUAN Hebing,LIU Xiaoli,WANG Enzhi,et al.The effect of discontinuity density on key strata threshold span in longwall panels[J].Journal of Geophysics and Engineering,2019,16(2):289-296.
[11] LIU Xiaoli,HAN Guofeng,WANG Enzhi,et al.Multiscale hierarchical analysis of rock mass and its mechanical and hydraulic properties prediction[J].Journal of Rock Mechanics and Rock Engineering,2018,10(4):694-702.
[12] LIANG Tao,LIU Xiaoli,WANG Sijing,et al.Study on the fractal characteristics of fracture network evolution induced by mining[J].Advances in Civil Engineering,2018:9589364.
[13] 缪协兴,王安,孙亚军.干旱半干旱矿区水资源保护性采煤基础与应用研究[J].岩石力学与工程学报,2009,28(2):217-227.
MIAO Xiexing,WANG An,SUN Yajun.Research on basic theory of mining with water resources protection and application to arid and semi-arid mining areas[J].Chinese Journal of Rock Mechanics and Enigeering,2009,28(2):217-227.
[14] 易立新,徐鹤.地下水数值模拟:GMS应用基础与实例[M].北京:化学工业出版社，2009.
[15] ALEXEY Chernyshenko,MAXIM Olshanskii.An unfitted finite element method for the Darcy problem in a fracture network[J].Journal of Computational and Applied Mathematics,2020,366:112424.
[16] 周永正,詹棠森,方成鸿，等.数学建模[M].上海:同济大学出版社，2010.
[17] LIU Xiaoli,WANG Sijing,WANG Enzhi.A study on the uplift mechanism of Tongjiezi dam using a coupled hydro-mechanical model[J].Engineering Geology,2011,117(1-2):134-150.
[18] 尚文绣，王远见，贾冬梅.基于天然水文情势的水库调度规则研究[J].人民黄河,2019,41(6):34-37.
SHANG Wenxiu,WANG Yuanjian,JIA Dongmei.Study on reservoir operation rules based on natural flow regime[J].Yellow River,2019,41(6):34-37.
[19] BARNABY Dobson,THORSTEN Wagener,FRANCESCA Pianosi.An argument-driven classification and comparison of reservoir operation optimization methods[J].Advances in Water Resources,2019,128:74-86.
[20] 艾学山，董祚，莫明珠.水库多目标调度模型及算法研究[J].水力发电学报,2017,36(12):19-27.
AI Xueshan,DONG Zuo,MO Mingzhu.Multi-objective reservoir operation model and algorithm[J].Journal of Hydroelectric Engineering,2017,36(12):19-27.
