榆神矿区地下水埋深上限阈值

马雄德，祁 浩，郭亮亮，迟宝锁，王宏科，朱占荣，曹虎生

(1.长安大学 水利与环境学院，陕西 西安 710054; 2.长安大学 旱区地下水文与生态效应教育部重点实验室，陕西 西安 710054; 3.自然资源部煤炭资源勘查与综合利用重点实验室，陕西 西安 710016; 4.陕西省一八五煤田地质有限公司，陕西 榆林 719000; 5.陕西陕煤陕北矿业有限公司,陕西 榆林 719000)

分布在黄河流域中上游的晋、陕、蒙、宁、甘等5省区的煤炭基地2019年产量27.49亿t，占全国原煤产量的73.4%(2020年中国原煤行业分析报告)，有力支撑了区域经济和国家能源需求。但该区地处我国西部干旱-半干旱区，降水稀少、水资源匮乏，主采煤层上覆隔水层厚度偏薄，煤层开采极易破坏上覆含水层和生态环境。据测算，在晋陕蒙地区，平均每生产1 t原煤直接损失2 t地下水资源[1]，已直接或间接对黄土高原丘陵沟壑水土保持生态功能区的水源涵养、防风固沙等功能产生了重要胁迫[2]。由于干旱少雨，补给量有限，水资源先天不足已经成为制约黄河中上游地区煤炭工业发展的资源瓶颈[3]。因此，提高对西部矿区高强度煤炭开采过程中水资源及生态破坏机理的认识，事关黄河流域高质量发展与生态保护之大计。

煤层采出以后，采空区地表会产生一定程度的下沉[4]，从而形成采矿沉陷盆地。依据隔水岩组厚度与导水裂隙带高度的关系，采煤对地下水位的扰动主要有2种表现形式。一种是煤层开采后地下水位埋深相对变浅，有时伴随着一定范围的地表积水。当导水裂隙带高度小于隔水岩组厚度，且地下潜水位较高时，地下水很容易在沉陷盆地中出露，形成人工地表水体[5]。这在我国东部地区比较常见，如安徽省两淮、山东省兖滕、江苏省徐州等矿区，地下水水位埋深较浅，而地表下沉系数大，每生产10 000 t煤，就会产生133～200 m2沉陷区[6]，近一半的沉陷区都会出现地表积水，从而会形成湿地、永久或季节性池塘、湖泊等人工地表水体[7]。由于沉陷区人工地表水体破坏了优质耕地[8]，改变了生态系统构成[9]，高潜水沉陷区大量研究主要围绕土地复垦[10]和生态恢复而展开[11]。

另一种形式是煤层开采后地下水位埋深增加。当煤层上覆隔水岩组厚度较小，大采高形成的导水裂隙带破坏潜水含水层隔水底板时，导致隔水岩组渗透率急剧增大，最多可达初始值的数十倍乃至百倍[12]，地下水迅速沿导水裂隙泻入矿井，使采空区周边一定范围地下水位降低[13]。这类问题在我国西部陕西、山西、内蒙古、宁夏等地区大型矿井开采过程中较为突出。范立民等[14]对比了1994年和2015年榆神府矿区地下水流场，指出该区内煤层开采已经造成地下水位大范围下降，水位降深大于8 m的区域面积达到758.9 km2。这进一步说明浅埋煤层开采引起的地下水位下降现象普遍存在，而且这种变化逐渐从孤立的煤矿区向整个水文单元甚至流域扩展[15]，正在剧烈改变着区域水循环模式。由于地处西部干旱-半干旱区，地下水埋深较浅时(小于4 m)，一些植被依赖地下水而生存[16]，地下水位下降后植被失去了稳定的供水水源，由于缺水而呈现生理性退化，因而地下水位下降问题被广泛关注。钱鸣高[17]、范立民[18]等呼吁在西部地区煤炭开采过程中要保护地下水，20 a以来，“保水采煤”已成了西部矿区最重要的命题之一。

最近，在榆神矿区西部地下水监测中发现:由于煤层开采地面沉降导致第四系地下水在地面成片出露，并能长期存在。类似地，杨志[19]发现金鸡滩煤矿6个监测孔水位在工作面经过以后都呈现不同幅度的上升(即水位埋深变浅)，最大上升了0.57 m。马立强等[20]也指出，补连塔煤矿开采后2个监测孔水位分别上升0.66 m和0.87 m。地下水位埋深变浅带来一个重要的问题，即地下水蒸发速率增加。地下水埋深越小蒸发速率越大[21]，地下水的无效蒸发量就越大。李莹[22]在计算榆神矿区地下水资源量时，为不同潜水埋深给定不同的蒸发系数，即潜水埋深<1 m，1～2 m和2～3 m，潜水蒸发系数分别为0.17，0.12和0.06。由此估算，如地下水位埋深从2～3 m上升到0～1 m，单位面积上蒸发量将增加183.3%，这非常不利于干旱地区浅层地下水资源管理。因此，控制合理的水位埋深上限已成为榆神矿区西部煤层开采中面临的新的科学问题。

为了减少蒸发，防治土壤盐渍化，我国许多地区开展了地下水位埋深上限研究。如，WANG等[23]认为有植物存在时三江平原地下水埋深上限为0.5～2.8 m。李明等[24]认为新疆南部阿克兰干地区土壤为粉砂夹黏土时地下水上限埋深为1.93 m，而土壤为粉砂时地下水埋深上限2.4 m。程双虎等[25]认为河北省平原区水位上限埋深为2 m。鉴于以往临界水位埋深的确定都是通过实测毛细上升高度或野外调查统计后获取的，结论往往难以统一。笔者通过潜水蒸发实验和理论推导:① 揭示风积沙区潜水蒸发过程;② 提出地下水埋深上限阈值计算公式;③ 确定榆神矿区地下水位上限阈值。

1 榆神矿区水位变化分区

根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》，隔水层厚度与导水裂隙带高度之差是否满足安全防水煤岩柱厚度，是判断含水层结构是否被破坏的标准，即

hsh≥hli+hb

(1)

式中,hsh为隔水岩组厚度;hli为导水裂隙带高度;hb为防水安全煤岩柱高度；隔水层为黏土时可取3倍采高，无黏土层时取5倍采高[26]。

式(1)两边同除以M，依此来表征采后隔水岩组隔水性稳定程度k，则有

(2)

式中，c为常数，取3或5。

范立民等[27]分析了40多个钻孔实测导水裂隙带高度，指出榆神府矿区采高3～6 m时，导水裂隙带高度一般为采厚的22～28倍，则有:

(1)k≥33时(隔水层为黏土层时为31)，煤层开采后隔水岩组隔水稳定性保持不变，含水层结构完整。但由于采后上覆岩层发生沉陷，地表高程降低，使潜水埋深相对变浅。在地势低洼处，地下水将会出露地表后形成地表水体。如小保当一号井开采后，地下水埋深变浅，甚至出露[28]。

(2)k≤22时，导水裂隙将完全贯通隔水岩组，含水层结构完全破坏，地下水位下降。如锦界煤矿SMJ107 号监测井水位最大降幅达12.67 m[14]。

(3)当22

据此，以采高5.0 m为例，隔水岩组厚度以22倍采高(110 m)和33倍采高(165 m)(有黏土层时为155 m)为标准，可将榆神矿区煤层开采后地下水位变化分为3个区，即水位埋深变浅区、过渡区和水位埋深增加区，如图1所示。

图1 榆神矿区水位变化预测分区Fig.1 Prediction zone of water level change in Yushen mining area

由图1可知，水位埋深变浅主要分区在榆神三期和四期规划区。该区首采煤层埋深300 m以上，上覆基岩隔水层厚度大，开采后导水裂隙带高度小于隔水岩组厚度，煤层开采后对潜水含水层结构影响小，因而开采沉陷区地下水位埋深一般会相对变浅。

图2为2019年6月统测的榆神矿区地下水位埋深图。由图2可知，榆神矿区西部风沙滩地区地下水位埋深普遍较小，其中水位埋深<4 m的面积为3 109.1 km2，占矿区面积的59.1%。据宋世杰等[29]的估算，榆神矿区埋深300 m，采厚5 m条件下，地表下沉系数平均为0.69，最大可达0.783，地表最大下沉值近3.9 m。因此，榆神三、四期地下水位小于4 m的区域都存在大范围地下水出露或浅埋的可能。

图2 榆神矿区地下水位埋深分区Fig.2 Zoning for depth to groundwater table

结合野外调查成果，水位埋深变浅区地表一般为厚度在10～146 m的风积沙，其粒径主要分布在0.074～0.250 mm。以下针对风积沙开展潜水蒸发试验，以获取风积沙区的地下水埋深上限阈值。

2 潜水蒸发规律

2.1 潜水蒸发试验

实验设备如图3所示。蒸渗柱采用直径为11 cm的有机玻璃制成。蒸渗柱与马氏瓶连接，用来记录蒸渗柱内的蒸发量。蒸渗柱顶部40 cm处放置红外灯以增加蒸发量，在等温条件下，用蒸馏水开展了一组初始饱和裸土蒸发实验(不定水位)和7组定水位蒸发实验(1组水面蒸发，6组潜水蒸发，水位埋深分别为0.1,0.3，0.4，0.5，0.6，0.7 m)。该实验是在一个小的、隔热良好的房间进行的，温度大致恒定(28± 2.7 ℃)，由湿度计测量的相对湿度在0.15～0.30。

图3 蒸渗仪组成Fig.3 Composition of lysimeter

实验用沙采自榆林市小壕兔乡(39.00°N，108.75°E)。参照国内外常用做法[30]，在研究区开挖1 m深的剖面，采集沙样运回实验室。沙样装填完成以后，为蒸渗柱连续充水和排水2～3次，以保证风积沙自然密实，使干容重与野外条件的实测值相差小于5%。

2.2 潜水蒸发过程

初始饱和风积沙的蒸发过程如图4所示。由图4可知，初始饱和裸土的蒸发过程中，蒸发速率存在明显的转折，大致可以分为2个阶段[31]。第1阶段(稳定蒸发阶段)蒸发速率较高，而第2阶段(水汽扩散阶段)蒸发速率十分缓慢，2者之间存在一个转化阶段。第1阶段的蒸发速率较高，SHOKRI[32]认为这是由于表土层中细小孔隙中的液体弯月面与大气保持联系，这些具有水力联系的小孔隙产生的毛细水力梯度将土壤内部的水分输送到地表以满足蒸发需求。LEHMANN等[33]认为当地下水降低到一定深度时，向下的重力和黏滞力将克服向上的毛细作用力，进而阻断了地下水与土表面之间的水力联系。这时，表层土的水分供应受限，因而蒸发面将由地表进入土壤内部，并在土壤表层形成一定厚度的干层[34]。此时土壤剖面上，水分以液态形式到达干层底部，以气态形式穿过干层进入大气，蒸发进入第2阶段，蒸发速率变得十分缓慢(图4)。

图4 毛乌素沙的蒸发过程Fig.4 Evaporation process of the Mu Us sand

2.3 不同埋深条件下潜水蒸发

从马氏瓶获得的读数换算到蒸渗仪的蒸发强度，图5为不同地下水位埋深条件下潜水蒸发损失情况。为了直观对比，将每个埋深的潜水累计蒸发量(∑m)用水面累计蒸发量(∑m0)进行了标准化。由图5可知，地下水位埋深从10 cm降低到30 cm时，总蒸发水损失减少了近20%，对总蒸发水损失的影响较小。然而，将地下水位埋深从30 cm降低到50 cm，总蒸发水损失减少了近90%，导致总蒸发水损失显著减少。这说明对于试验沙而言，地下水位埋深0.5 m左右时，蒸发过程中的水分传输机制发生了转变，地表开始形成干层，蒸发进入第2阶段。

图5 地下水位与潜水蒸发关系Fig.5 Relationship between groundwater level and evaporation rate

3 矿区地下水埋深上限阈值的确定方法

3.1 矿区地下水埋深上限阈值的概念

在常温常压下，毛细上升高度与毛细管直径关系[35]为

(3)

式中,hc为毛细上升高度;d为毛细管直径。

可以看出，毛细管直径是决定毛细上升高度的主要因素，毛细管直径越小，毛细上升高度越大。从土壤孔隙失水的角度来看[36]，蒸发过程中毛细上升高度范围内的大孔隙先于小孔隙失水，因此在蒸发第1阶段，一定深度范围的土壤剖面上大孔隙内的水分完全耗散后，在毛细力的作用下小孔隙之间仍然存在持续的水力联系，只要地下水位保持不变，地下水与地表之间的毛管通量就不受限，因而蒸发速率较高(图6(a))。随着地下水位埋深下降，小孔隙中的水分也完全散失，此时毛管水力联系断裂，毛管流停止，蒸发进入第2阶段，潜水的蒸发速率急剧减小。因此，可以认为在土壤蒸发过程中存在一个临界值Dcrit(图6(b))，当地下水位埋深大于Dcrit时，表层土与地下水之间的毛细水力联系断裂，蒸发进入第2阶段，蒸发速率十分缓慢。

图6 不同蒸发阶段土壤剖面上水分分布示意Fig.6 Water distribution in the soil profile at different evaporation stages

干旱矿区煤层开采导致地下水位埋深d相对变浅时，如果d

3.2 地下水埋深上限阈值的确定

当均质土处于稳定蒸发时，以潜水面为原点，向上为正，其水分运移方程可表示为

(4)

式中，E为蒸发强度;K(h)为非饱和导水率;h为基质势;z为位置变量。

对方程(4)两边积分，得

(5)

要求得式(5)的解，就必须给定非饱和渗透系数的表达式。由于地下水面以上毛管边缘带土壤接近饱和状态，因此将土壤剖面以进气压力值划分为非饱和带和饱和带[37-38]。当0≤h≤hb，用Ks(饱和含水量对应的导水率)代替非饱和导水率K(h)。于是有

(6)

式中，h0为毛细管上升的最终湿润锋位置处的土壤吸力水头。

假定当地下水位埋深为Dcrit时表层土与地下水之间最小孔隙毛管水力联系中断，此时表层土的基质势为hmax，由图7可知，hmax=hb+hc,则有

(7)

Gardner指数水力传导模型[39]因其指数性质而在文献中被广泛使用，是因为这便于求得解析解。Gardner指数水力传导模型为

K(h)=Ksexp[-αG(h-hb)]

(8)

式中，αG为拟合参数，代表介质的孔径分布。

将式(8)代入式(7)，整理得

(9)

据前文，蒸发第2阶段，土壤剖面上最小孔隙中的毛管水断裂后，潜水蒸发强度E非常小,远低于饱和渗透系数Ks，即E/Ks≈0，代入式(9)有

Dcrit≈hb+hc

(10)

式(10)即为毛细水力联系断裂时临界地下水位埋深，即地下水埋深上限阈值的计算公式。可以看出，地下水埋深上限阈值与进气压力值和毛细上升高度有关。在蒸发第2阶段，土壤水分运动速率十分缓慢，因而可以忽略黏滞力[32]，毛细水力联系断裂时临界地下水位埋深在数值上等于表层土与地下水之间毛管水力联系中断时表层土的基质势为hmax(以长度单位表示)，如图7所示。

图7 典型的水分特征曲线概念模型Fig.7 Conceptual model of water characteristic curve

根据MCQUEEN和MILLER[40]提出的土壤水分特征曲线形状的概念模型，土壤基质吸力与含水率在半对数坐标系中由3段直线构成，其中进气压力值与残余含水率对应基质吸力之间的线段为毛细作用区，因而可以通过水分特征曲线来确定毛细水力联系断裂时临界基质吸力hmax，而量化进气值和残余含水率可以通过在水分特征曲线的拐点上构建双切线来完成[41]。

定义一个无量纲的含水率变量Θ，由饱和含水率和残余含水率定义成标准化体积含水率的形式为

(11)

式中，θ为体积含水率;θr和θs分别为残余含水率和饱和含水率。

可知，随着体积含水率逐渐趋近于残余含水率θr，标准化含水率Θ趋于0;若体积含水率趋于θs，则标准化含水率Θ趋近1。通过土壤水分特征曲线在θr处的转折点做切线，切线与饱和度Θ=0(对应残余含水率)的交点等于hmax[33]，如图7所示。

以van Genuchten方程[42]来描述土壤水分特征曲线，其表达式为

(12)

式中，α，n为拟合参数。

为了求得土壤水分特征曲线的拐点及切线的斜率，需要对式(12)求导，即求一阶导数得到切线的斜率，求二阶导数得到拐点坐标。由此可以写出过拐点的切线表达式:

(13)

于是，求式(13)与Θ=1的交点，得

(14)

式(14)给出了以van Genuchten方程为基础的矿区地下水埋深上限阈值通用计算方法，对干旱半干旱地区需要减少地下水蒸发，降低土壤盐渍化风险的区域都适用。其中，参数α,n为van Genuchten方程中的拟合参数，α与土的进气值有关，n与土的孔径分布有关，不能直接测得，一般通过实测土壤一系列含水率θ及其对应的基质势h后，利用方程(12)拟合得到。

4 榆神矿区地下水埋深上限阈值

在榆神煤矿三期规划区随机选择5个采样点，地表下0～80 cm采集土壤样品10 cm3;将每个点采集的土壤样品搅拌均匀，采用砂性漏斗法测量各样品在脱湿过程中的水头差和排水体积，见表1。实验结束后采用烘干法测定土壤含水率，再根据表1中排水量计算出每次降低出口时(对应土壤吸力)的含水率。

表1 出水口位置下降高度与排水量统计Table 1 Statistical table of outlet falling height and discharge

采用最小二乘法将换算的土壤吸力与含水率数据与van Genuchten方程进行拟合，得到拟合参数α=0.028 5，n=5.93。代入式(14)，得hmax=49.3 cm。

不同地下水位埋深条件下潜水蒸发速率如图8所示。图8中的阴影部分是利用式(14)计算的地下水埋深上限，即49.3 cm。在该深度范围内，地下水与表土之间有水力联系。当地下水位埋深大于该深度以后，地下水与表土之间的水力联系断裂，蒸发速率明显减弱。

蒸发实验结果表明(图8)，蒸发速率随地下水位埋深增加呈递减规律，地下水位埋深0.1 m时地下水蒸发速率约为2.35 mm/h，而当地下水位埋深大于0.5 m时，蒸发速率接近0.01 mm/h。说明对于毛乌素风积沙而言，地下水位埋深0.5 m时，地下水与表土之间的水力联系断裂，蒸发进入第2阶段。

图8 地下水位与蒸发速率关系Fig.8 Relationship between groundwater level and evaporation rate

杨泽元等[43]认为毛乌素风积沙干表层形成的临界地下水位埋深为50 cm。何渊[44]在乌审旗河南乡气象站开展了潜水蒸发原位实验，并指出，对于风积沙而言，潜水埋深40 cm是蒸发速率的转折点，潜水埋深50 cm被当作地下水极限蒸发深度。上述研究成果与本文的结论基本一致。据此，笔者建议将地下水位埋深50 cm作为榆神矿区煤层开采过程中地下水埋深上限阈值。

5 结 论

(1)以隔水岩组厚度22倍采高(110 m)和33倍采高(165 m)(有黏土层时为155 m)为标准，可将榆神矿区煤层开采后地下水位变化分为3个区，即水位埋深变浅区、过渡区和水位埋深增加区。

(2)榆神矿区潜水位埋深变浅区主要分布在榆神三、四期规划区，地表组成物质为风积沙。风积沙的蒸发过程分为2个阶段，即“稳定蒸发”阶段和“水汽扩散”阶段。水位埋深50 cm左右时，蒸发进入“水汽扩散”阶段，蒸发速率趋近于0。

(3)通过在水分特征曲线的转折点处构建双切线，推导了求取地下水埋深上限阈值的解析公式。地下水埋深上限阈值与进气压力值和毛细上升高度有关。

(4)利用毛乌素风积沙进行了脱水试验，根据实测含水率-负压数据进行拟合，求取了相关参数，确定了榆神矿区煤层开采过程中地下埋深上限阈值为50 cm。