生态脆弱矿区高强度植被恢复对地下水补给的影响

王强民,董书宁,王文科,王 皓

(1.煤炭科学研究总院，北京 100013; 2.中煤科工集团西安研究院有限公司，陕西 西安 710054; 3.陕西省煤矿水害防治技术重点实验室，陕西 西安 710054; 4.长安大学 环境科学与工程学院，陕西 西安 710054; 5.旱区地下水文与生态效应教育部重点实验室，陕西 西安 710054)

榆神矿区是我国重要的煤炭生产基地和高强度开采区，同时也是最为典型的生态环境脆弱区[1]，采矿活动、植被恢复及气候变化控制着矿区的整体生态格局[2-4]。受采煤地下水位下降、采煤塌陷及矿区土地类型转换等因素影响，煤炭开采过程中出现了植被退化的现象，为避免植被的退化及死亡，矿区开展了高强度的植被恢复工程。研究表明，截至2006年底，神东矿区累计投入生态环境建设经费2.69亿元，使得由开发初期(1985年)植被覆盖率3%～11%增长至59.4%以上(2005年)[5]，矿区和井田尺度的植被覆盖度大幅提升[6-8]，植被对降水的截留量和蒸腾量也随之增大[9-10]，造成降水对地下水的补给量减小。而研究区的主要地貌类型为风积沙，透水性好，降雨入渗补给量占地下水总补给量的78.5%[11]，且区内河流40%～70%的流量来自地下水的补给[12]，水文生态环境与地下水的关系密切。因此，分析高强度植被恢复对地下水补给的影响已成为生态脆弱矿区亟待解决的关键科学问题之一。

目前，在生态脆弱矿区地下水和植被关系的研究方面，主要集中在地下水状态变化对植被的影响[13-14]。例如通过开展野外原位监测试验，研究了采煤塌陷区土壤水的变异规律及生态效应[13];结合原位试验和数值模拟，分析了植被对采煤地下水位下降的响应机制[14]。近年来，针对植被变化对地下水的影响也开展了部分研究[9,15-16]，但研究整体偏少。例如，张俊等通过建立区域尺度剖面二维模型，分析了鄂尔多斯盆地植被变化对地下水流系统的影响[15];张建军研究了黄土高原水土保持对土壤水分的影响，认为人工刺槐的耗水量大，甚至大于当地降水量，导致深层土壤干化[16]，地下水很难得到补给[9]。但是，以上研究主要集中在黄土高原区，而榆神矿区植被和地下水的相互反馈作用较黄土高原区强烈，研究结果不能很好的指导生态脆弱矿区植被恢复实践。

笔者结合野外裸土区和植被覆盖区2个原位试验点的观测数据，以2种模式下的降水—土壤水—地下水转化过程为研究对象。通过遥感解译分析矿区的生态环境变化情况;采用原位监测手段分析植被对土壤水、地下水位的影响规律;在此基础上，利用数值仿真计算分析植被覆盖度增加对降水补给地下水的影响机制。研究成果对丰富我国西部生态脆弱矿区植被恢复与地下水资源的协调发展具有重要意义。

1 研究区概况

榆神矿区位于陕西省神木县榆阳区境内，煤炭资源储量丰富、地质构造简单、开采技术条件优越，其延安组含煤地层自上而下赋存6层可采煤层。主采2-2煤层，平均埋深261.15 m，可采厚度0.80～12.49 m，平均厚度达到10.04 m[17]。受隔水层性质、煤层厚度和埋深等因素影响，采动裂隙发育至基岩含水层内部甚至近地表的松散层含水层，改变了地下水和土壤水的运移规律，表现为地下水位下降和土壤含水率亏缺，继而加剧了矿区水资源短缺和生态环境脆弱的现状[18]。近年来，为了落实绿色矿山建设要求，各大煤矿开展了广泛的矿山植被恢复工程，矿区植被覆盖度得到整体提升。为评估矿区生态恢复植被覆盖度增加对地下水补给的影响，研究团队在榆神矿区无定河流域(距榆林市40 km)搭建了原位试验平台，试验区多年平均降水量为400 mm，平均潜在蒸散发量为1 153 mm，地貌以风积沙为主，沙柳、沙蒿为主要植被类型。

2 研究方法

2.1 原位试验

本次试验在原位试验平台的植被耗水区进行，该区设置玻璃钢材质的圆柱形监测井2个，分别为裸土区和植被覆盖区2个原位试验监测井。监测井高度为4 m，直径为2 m，侧向和底部为隔水边界，上部接受降水补给和蒸散发排泄。为合理控制和观测地下水位，监测井中间设置地下水注(抽)水孔和地下水位观测孔。植被选取研究区的优势植被沙柳作为研究对象，初始地下水位埋深设置为3 m(极限蒸发深度)，试验介质采用矿区典型的风积沙。含水率测定采用Decagon公司生产的ECH2O土壤水分传感器，地下水位测定采用斯伦贝谢水务公司生产的Mini-Diver传感器测定，测定频率均为1 h/次。裸土区含水率传感器在垂直方向上的布置为:3,10,20,30,50,80,150,250,350 cm，由于填土均匀，裸土区水平方向不设置对照传感器。植被覆盖区含水率传感器在垂直方向的布置和裸土区一致，考虑植被耗水对水平方向土壤水的影响，水平方向设置对照传感器，分别距离植被中心点0,20,50,80 cm(图1)。

图1 原位试验示意(裸土区+植被覆盖区)

2.2 数值模拟

在原位试验的基础上，笔者采用美国国家盐改中心开发的Hydrus软件[19]模拟裸土和不同植被覆盖条件下的降水—土壤水—地下水相互转化过程，并求解植被耗水量和降水入渗系数。植被覆盖条件下，一维非饱和土壤水运动的控制方程为

(1)

式中，C(h)为容水度，cm-1;h为土壤水压力水头，cm;t为时间，d;z为垂向坐标，cm，向上为正;K(h)为非饱和渗透系数，cm/d;S(z,t)为植被根系吸水速率,裸土条件下，S(z,t)为0。

文中采用van-Genuchten-Mualem模型[20]描述土壤水分特征曲线和渗透系数曲线:

(2)

(3)

(4)

式中，θs为饱和含水量，cm3/cm3;θr为残余含水量，cm3/cm3;Ks为饱和渗透系数，cm/d;Se为有效饱和度；α,m,n为相关土壤参数，其中m=1-1/n，n>1;l为弯曲度参数。

笔者采用Feddes模型[21]描述植被的根系吸水过程:

S(z,t)=α(h)b(z)Tp

(5)

式中，α(h)为水分胁迫函数;b(z)为根系分布函数;Tp为植被潜在蒸腾量，cm/d,可由下式计算:

Tp=ETp(1-e-kLAI)

(6)

式中，ETp为潜在蒸散发量，cm/d;k为冠层的消光系数;LAI为叶面积指数。

模型构建过程中的初始条件(土壤水数据)、上边界条件(气象数据)、下边界条件(地下水位数据)来源于本次原位试验的实时监测，土壤水分特征参数和植被根系密度函数见文献[3]。在模拟计算之前，本文根据原位试验的观测结果，对相关参数进行了识别和验证，30,50和80 cm处的土壤水的观测值和模拟值动态变化趋势一致(图2)。利用式(7)，(8)求得相对误差和均方根误差，见表1。

(7)

(8)

式中，RE为相对误差，%;RMSE为均方根误差，cm3/cm3;n为模拟天数;Si和Oi分别为第i天的模拟值和实测值，cm3/cm3。

图2 模型识别和验证过程土壤水的观测值和模拟值(2016年)

表1 模型识别和验证过程中的误差分析

由表1可知模型识别和验证过程中土壤水模拟值与实测值的相对误差<10%、均方根误差<0.01 cm3/cm3，反映了基于野外观测数据构建的数值模型可以较好的分析研究区裸土和植被覆盖条件下的降水—土壤水—地下水相互转化过程。

3 结果分析

3.1 矿区生态环境演变特征

文中选用数据质量高、信息量大的MODIS数据产品研究榆神矿区生态环境演变特征，其时间分辨率为16 d，空间分辨率为250 m。为增强数据的可靠程度，分别选取研究区8月份植被指数(NDVI)和植被净初级生产力(NPP)两个指标联合反映矿区生态环境演变特征，其中NPP值采用光能利用率模型(CASA)估算[22]。研究结果显示，矿区尺度的NDVI值和NPP值都呈持续增加现象，NDVI值从2000年的0.247增长到2015年的0.358，NPP值从2000年的12.67 g/m2增长到2015年的20.78 g/m2(图3)，升高幅度分别为45%和64%，表现为低覆盖区连续减小，中、高(较高)覆盖区持续增大(表2)，反映矿区生态环境持续好转。

图3 榆神矿区平均NPP和NDVI变化特征

3.2 植被影响下的土壤水变化规律

根据原位试验观测结果，分别统计相同地下水埋深和气象条件下裸土区和植被覆盖区0～80 cm内的土壤水监测数据，并绘制裸土区(一维)和植被覆盖区(二维)月平均土壤含水率分布特征图(图4)。裸土区平均土壤含水率为0.101 cm3/cm3，受气象条件和地下水位的双重影响，含水率随埋深的增大呈增加趋势;受植被根系耗水影响，植被覆盖区0～80 cm内的平均土壤含水率仅为0.062 cm3/cm3，较裸土区减少38.6%。

表2 不同植被覆盖度条件下的面积分布比例

图4 裸土区和植被覆盖区平均土壤含水率分布特征

为更好的分析不同条件下土壤水分的运移规律，利用风积沙的土壤水分特征曲线，将裸土区(一维)和植被覆盖区(二维)0～80 cm内土壤含水率换算成土壤基质势，利用式(9)计算土壤总水势。

H=h(θ)+Z

(9)

式中，H为土壤总水势，cm;h(θ)为土壤基质势，cm;θ为土壤含水率，cm3/cm3;Z为重力势，cm。

裸土区(一维)和植被覆盖区(二维)土壤总水势分布特征如图5所示(箭头为水分运移方向)。裸土区的平均土壤总水势为-92.64 cm，20 cm处总水势最大(-75.54 cm)，表现为土壤水自20 cm处向两侧运移，较易接受大气降水的补给;植被覆盖区平均土壤水势为-146.88 cm，远低于裸土区的土壤水势，在根系密集区(20 cm，50 cm)出现水势低点(-170 cm)，表现为上部、下部和侧向的土壤水向植被根系处运移，不易接受大气降水的补给。

图5 裸土区和植被覆盖区平均土壤总水势分布特征

3.3 植被影响下的地下水变化规律

图6反映了裸土区和植被覆盖区的地下水位埋深动态变化特征。在无降水期，裸土区地下水位呈现出缓慢的上升趋势，表明尽管受土壤蒸发作用，深部土壤水仍维持向下运动的趋势;受植被蒸腾耗水影响，植被覆盖区的地下水位平均以1.87 cm/d的速度呈现持续下降的趋势。在累计124 mm降水量的情景下，裸土区的地下水位表现出明显的上升趋势，上升幅度高达27 cm;相反，植被覆盖区的地下水位没有回升，呈现出较低速度(0.29 cm/d)的地下水位下降现象，表明植被蒸腾耗水缩减了降水入渗补给量。

图6 裸土区和植被覆盖区地下水埋深变化

基于地下水位动态观测数据，本次试验地下水蒸发(补给)强度可由式(10)估算:

Eg(Rg)=±μΔh

(10)

式中，Eg(Rg)为地下水蒸发(补给)强度，cm/d;μ为风积沙给水度，取值为0.26[23];Δh为地下水位降(升)幅，cm。

图7反映了裸土区和植被覆盖区的蒸发和补给强度动态变化特征。通常，裸土区的地下水以低强度的入渗补给为主，平均补给强度为0.035 cm/d，降水时补给强度大幅增加;植被覆盖区地下水表现为高强度的蒸发现象，平均蒸发强度为0.48 cm/d，降水时蒸发强度大幅减小。受植被蒸腾耗水影响，地下水面水分垂向交换量由-0.035 cm/d(入渗)增加至0.480 cm/d(排泄)，水分运移方向和水分垂向交换量都发生明显改变。由此可见，植被的存在一定程度改变了地下水的补给和排泄条件，从而引发地下水位的下降。

图7 裸土区和植被覆盖区地下水蒸发(补给)强度变化

3.4 植被覆盖度增加对地下水补给的影响

高强度植被恢复的外在表现是植被覆盖度的增加，在数值模型里面概化为植被叶面积指数的升高。杜春雨等[24]运用归一化植被指数(NDVI)推求植被叶面积指数(LAI)，并构建了线性模型LAI=7.67NDVI-4.01，其相关系数为0.66。结合遥感数据，矿区最大的NDVI值为0.98，基于上述NDVI和LAI的经验公式，推测最大LAI为3.50。笔者定义裸土条件下植被的叶面积指数(LAI)为0，模拟过程中以0.5/次的幅度逐步增加植被的叶面积指数，直至叶面积指数为3.50。考虑到当地下水埋深大于3 m时，榆神矿区降水入渗系数随埋深的增大不发生变化[25]，故本文采用3 m处的土壤水流通量估算降水入渗系数，定义为

α′=Fs/P

(11)

式中，α′为降水入渗系数;Fs为平均土壤水通量，mm;P为平均降水量，mm。

图8反映了不同叶面积指数条件下的植被耗水量和降水入渗系数变化规律。随着叶面积指数的增加，植被耗水量呈增加趋势，由水均衡原理可知土壤中向下运移的水流通量逐步减小，从而导致降水入渗补给量减小，表现出随着叶面积指数的增加，降水入渗系数逐步减小。数值模拟结果表明，裸土条件下(LAI=0)风积沙的降水入渗系数为0.54，随着叶面积指数由0增大至3.5，入渗系数减小至0.198，减小幅度高达63.3%。

图8 植被耗水量和降水入渗系数随叶面积指数的变化

4 讨 论

近年来，随着生态脆弱矿区植被恢复工程的实施，区域生态环境质量得到改善，植被覆盖度持续增加，这种现象在我国西部干旱缺水矿区尤为显著(图3)。榆神矿区地表由大面积的风积沙组成，本文研究结果显示，无植被覆盖条件下风积沙的降水入渗系数高达0.54，降水入渗直接补给地下水，使得该区域地下水埋深普遍较小，表现为地下水与植被、地表水系等关系密切[26]。在矿区植被覆盖度大幅提升的背景下，植被蒸腾耗水量的增加导致土壤水分大幅度降低[27]，从而减少降水对地下水的垂向补给[28]。结合原位试验观测和数值模拟，尹立河等认为当植被覆盖率增加至100%时，降水对地下水的补给量基本为0[9];KEESE K E通过搭建HYDRUS数值模型研究了植被对地下水补给的影响，认为当干旱区植被覆盖度增加后，地下水补给量降低幅度最高可达30倍[29]。结合本文研究成果，原位试验显示植被根系区的土壤含水率出现较大范围的低值区(图4)，造成土壤总水势的分布状况和大小出现变异(图5)，使得水分运移方向和垂向交换量发生改变，表现出裸土区和植被覆盖区截然不同的地下水位动态特征(图6);数值模拟结果表明，随着植被叶面积指数增加植被耗水量呈增加趋势，降水入渗系数呈减小趋势(图8)。因此，榆神矿区在植被恢复的工程实践中，应优先考虑植被类型和植被覆盖度两个基本要素，选择耗水量较小的植被类型和适宜的植被覆盖度，避免出现由于区域生态耗水量过大引发的降水入渗补给量大幅减小的现象发生。值得注意的是，笔者在研究中没有考虑植被降水截留对地下水补给的影响，徐先英等研究发现灌木对降水的截留量可达30%以上[30]，这同样会减少降水对地下水的补给量，因此本文的数值计算结果相对保守。

5 结 论

(1)矿区尺度的NDVI和NPP都呈持续增加现象，NDVI值从2000年的0.247增长到2015年的0.358，NPP值从2000年的12.67 g/(m2·a)增长到2015年的20.78 g/(m2·a)，表现为低覆盖区连续减小，中、高覆盖区持续增大。

(2)裸土区20 cm处土壤总水势最大，表现为土壤水自20 cm处向上、下两侧运移，较易接受大气降水的补给;植被覆盖区土壤总水势远低于裸土区的土壤水势，在根系密集区出现水势低点，表现为上部、下部和侧向的土壤水向植被根系处运移，不易接受大气降水的补给。

(3)裸土区地下水位呈现出缓慢上升的趋势，以低强度的入渗补给为主;植被覆盖区地下水位呈现出明显的持续下降趋势，以高强度的蒸发排泄为主。受植被蒸腾耗水影响，地下水面水分垂向交换量由-0.035 cm/d(入渗)增加至0.480 cm/d(排泄)，水分运移方向和水分垂向交换量都发生改变。

(4)随着叶面积指数的增加植被耗水量呈增加趋势，土壤中向下运移的水流通量逐步减小，导地下水接受降水的补给量减小，表现出随着叶面积指数的增加，降水入渗系数逐步减小。数值模拟结果表明，裸土条件下风积沙的入渗系数为0.54，随着叶面积指数由0增大至3.5，入渗系数减小至0.198，减小幅度高达63.3%。

(5)榆神矿区在植被恢复的工程实践中，应优先考虑植被类型和植被覆盖度两个基本要素，选择耗水量较小的植被类型和适宜的植被覆盖度，避免出现由于区域生态耗水量过大引发的降水入渗补给量大幅减小的现象发生，以实现生态脆弱矿区植被恢复与地下水资源保护的协调发展。