神东矿区降雨入渗补给强度研究*

何 渊

(中煤科工集团西安研究院有限公司，陕西 西安710054)

0 引 言

神东矿区地处陕北黄土高原北部、鄂尔多斯高原毛乌素沙漠的东南缘。该区属干旱半干旱温带大陆性气候，降雨稀少，蒸发强烈，生态环境脆弱，矿区的发展带来了一系列水资源及生态环境问题，合理开发利用地下水资源需要对区域的地下水资源做出科学的评价，已有研究表明，区域地下水循环方式主要以垂向交换为主［1］。作为地下水资源评价中补给量的主要构成部分，以往的计算大多采用经验公式计算的方法简单的将降雨入渗补给地下水量概化为降雨入渗系数，误差较大。本次研究采用野外调查包气带岩性结构及取样，采用数值法模拟降雨入渗过程，模型中参数由实验室测定给出，对于土样的实验室测定主要采用离心机、环刀法及马尔文激光颗粒分析仪测定各项土壤水力学参数，在此基础上结合当地的气象资料建立研究区包气带数值模拟模型，计算出区域不同地段降雨入渗强度及入渗系数，同时得出的参数也可用于类似区域的水文地质计算［2］，为地下水资源评价提供科学依据。

1 野外调查及取样

项目组于2013 年在研究区由南至北垂直乌兰木伦河按地质地貌特征布置了3 条剖面线进行包气带岩性调查，同时取样。取样范围基本控制研究区域，每条剖面布置4 个取样点，沿剖面布置的每个取样点取包气带埋深0 ～3 m 颗分样，同时在示范矿井(补连塔与大柳塔煤矿)范围及包气带岩性(风积沙和黄土状图)典型区域取0 ～3 m 包气带原状样，每个深度原状样均取双样。

2 研究区包气带岩性分布

根据野外调查及结合前人已做过的地质地貌图，遥感解译图等资料分析，研究区在地貌单元上可划分为风沙地貌和黄土地貌2 大类。以风沙地貌(土壤质地分类为砂土、壤质砂土)为主，覆盖区内90%以上面积，其地表形态表现为风沙滩地及沙丘-沙梁，在风沙滩地区，湖泊海子星罗棋布;黄土地貌(土壤质地分类为砂质壤土、壤土为主)仅在区内东南部沿沟壑地带分布，地形切割较深，冲沟发育，纵横交错，沟谷陡峻狭窄，植被稀少，水土流失严重［3-5］。

图1 水文地质概念模型示意图Fig.1 Hydrogeological conceptual model

3 包气带水分运移数值模拟

3.1 水文地质概念模型

根据野外典型点取样的包气带岩性结构概化出的水文地质概念模型如图1 所示。由于所计算区域包气带中水分运移主要是垂向运动，因此，将其概化为垂向一维流［6-8］。试验介质概化为2 层，z 轴坐标原点选在地表，向上为正。已有研究认为，在包气带厚度小于潜水蒸发极限深度的条件下，随着包气带厚度的增大，降雨入渗强度逐渐减小，大于3 m 后，年降雨入渗补给强度将趋于稳定，与埋深无明显相关关系。根据本次研究收集区域内地质、水文地质资料分析，神东矿区内地下水位埋深除沟谷地段外基本均大于3 m，故本次数值模拟的包气带计算剖面厚度概化为3 m.

3.2 参数获取

3.2.1 土壤水分特征曲线

此次土壤水力学参数获取的试验采用日本himac 公司生产的CR21G 离心机测定土壤水分特征曲线，测定时温度恒定20 ℃.对于原状土，首先将标准环刀(面积100 cm2)取得原状土样浸水饱和，并测定其饱和含水量。然后将饱和的原状土样品放入离心机装置中，根据设计压力设定离心机转速，加压范围为10 ～1 000 kPa，分别为10，20，40，60，80，100，200，400，600，800，1 000 kPa.在每个压力达到平衡后，取出样品通过称重获得质量含水量。试验测得8 个原状样土壤水分特征曲线如图2 所示。

本次计算利用RETC 软件对8 个原状样利用实验室测得的土壤水分特征曲线进行参数反演，得出8 个原状样的土壤水力学参数。其结果见表1.

表1 土壤水力学参数测定结果表Tab.1 Soilhydraulic parameter(original state samples)

图2 原状样土壤水分特征曲线Fig.2 Soil water characteristic curve of the original state sample

3.2.2 土壤饱和导水率

土壤饱和导水率测定采用环刀法，用环刀取原状土样，浸水后，在单位水压梯度下，根据达西定律，求得通过垂直于水流方向的单位土壤截面积的水流速度，称为土壤的饱和导水率［9-10］。测得结果见表2.

表2 原状样土壤饱和导水率测定结果Tab.2 Saturated Soil hydraulic conductivity of original state sample

3.2.3 土壤颗粒分析

颗粒分析采用世界通用的马尔文激光颗粒分析仪测定(型号为MS2000)，该分析仪是采用激光衍射法来测定土壤颗粒之间的粒度大小及粒度分级。

本次计算中颗粒分析土样的土壤水力参数获取是利用RETC 软件中的土壤转换函数功能，即根据土壤的颗粒级配中砂粒、粉粒、粘粒的百分含量以及土壤容重等土壤物理性质数据，直接输出Van-Genuchten 模型中的参数。

颗粒分析土样参数反演结果见表3.

3.3 数学模型

式中 θ 为土壤体积含水率，cm3·cm-3;k 为非饱和渗透系数，cm·hr-1;t 为时间变量，hr-1;z 为空间变量，cm;h 为压力水头(在包气带中称为负压);地表为原点，向上为正。

表3 颗粒分析样土壤水力学参数测定结果表Tab.3 Soil hydraulic parameter(particle-size analysis of soil samples)

)

3.4 模型的初始条件和边界条件

模型的初始条件设置为计算剖面的负压值，假定地表为负压值为0 m，计算剖面底部(3 m)处负压值为-3 m，线性分布。

模型的上边界条件由气象条件决定(本次计算采用研究区内气象站近20 a 的气象资料)，HYDRUS-1D 软件采用能量平衡方程将连续的气象资料转换为地表的能量通量。并计算出地表蒸发作为模型的上边界条件。模型下边界条件考虑为自由排泄边界［11］。

3.5 数学模型的求解

上述数学模型求解采用美国盐土实验室开发的hydrus-1D 软件解算此模型。数值离散采用有限元法，由程序自动完成。长度单位定为cm，时间单位定为d.初始时间步长为0.01 d，计算过程中，程序会根据计算收敛的情况自定调整时间步长，设定最小时间步长为0.000 01 d，最大时间步长为1 d［12-16］。

经数值模拟计算，研究区95%降水频率年地下水入渗补给量计算的结果见表5.

设地下水获得的降雨入渗补给量为pr，大气降水量记为p，降雨入渗补给系数为，根据降雨入渗补给系数的定义，有α =pr/p. 计算得出区域各计算点地下水不同的入渗补给系数，见表4.

表4 研究区降雨入渗强度及入渗系数汇总表Tab.4 Rainfall infiltration intensity and coefficient in study area

4 结 论

1)研究区在地貌单元上可划分为风沙地貌和黄土地貌两大类。在风沙滩地区以风沙地貌(土壤质地分类为砂土、壤质砂土)为主，覆盖区内90%以上面积，其地表形态表现为风沙滩地及沙丘，沙梁;黄土地貌(土壤质地分类为砂质壤土、壤土为主)仅在区内东南部沿沟壑地带分布。

2)本次数值模拟的包气带计算剖面厚度取为3 m，结合当地的气象条件及典型点岩性结构计算研究区降雨入渗强度及入渗补给系数。从研究区降雨入渗系数计算结果可以看出，区域降雨入渗系数大致范围在18.4% ～27.7%.

3)区域降雨入渗强度主要受到当地气象条件、包气带岩性结构及地下水位埋深(即包气带厚度)的影响。由本次计算结果可以看出，在相同降雨特征和水位埋深条件下，粗颗粒岩性的降雨入渗系数大于细颗粒岩性，即砂土＞壤质砂土＞砂质壤土＞壤土＞粉壤土。所以在研究区内地表广泛分布的风积沙对降雨入渗十分有利，对地下水资源的保护起到积极的作用。

4)从计算的结果来看，在研究区除沟谷地段只要地下水埋深大于潜水蒸发的极限埋深，潜水蒸发微弱，对潜水面来讲降水补给是主要的，研究区广泛分布的风积砂对区域降雨入渗起到了积极的作用。

5)文中拟合出的研究区各岩性水文地质参数可用于该区域及类似区域的水文地质计算。

6)包气带水分运移规律受多种因素影响，本次研究仅考虑了降雨量、蒸发量、包气带岩性、厚度等因素的影响，对于气、热、植被等的影响因素未来尚需开展专门的研究工作。

References

［1］ 何 渊.鄂尔多斯盆地沙漠高原区湖泊和潜水面蒸发能力研究［D］.西安:长安大学，2006.HE Yuan.Research on evaporation intensity of lakes and groundwater in desert plateau of Ordos Basin［D］. Xi’an:Chang’an University，2006.

［2］ 黄金廷.鄂尔多斯盆地沙漠高原区降雨入渗补给地下水研究［D］.西安:长安大学，2006.HUANG Jin-ting.Research on rainfall recharge in desert plateau of Ordos Basin［D］.Xi’an:Chang’an University，2006.

［3］ 肖 晖，李建新，韩 伟，等.鄂尔多斯盆地南缘渭北隆起中新生代构造抬升及演化［J］.西安科技大学学报，2013，33(5):359 -363.XIAO Hui，LI Jian-xin，HAN Wei，et al. The tectonic uplift time and evolution characteristics of Weibei uplift in the south edge of Ordos Basin［J］.Journal of Xi’an University of Science and Technology，2013，33(5):359-363.

［4］ 崔宏俊，刘宝平，欧阳畹町.鄂尔多斯盆地延长气田石盒子组沉积特征分析［J］. 西安科技大学学报，2012，32(6):726 -731.CUI Hong-jun，LIU Bao-ping，OUYANG Wan-ding. Analysis of the sedimentary characteristics of Shihezi formation in Yanchang gasfield，Ordos Basin［J］.Journal of Xi’an University of Science and Technology，2012，32(6):726 -731.

［5］ 王 云，梁 明，汪桂生.基于ArcGIS 的流域水文特征分析［J］.西安科技大学学报，2012，32(5):581 -585.WANG Yun，LIANG Ming，WANG Gui-sheng. Analysis of hydrological features based on ArcGIS［J］.Journal of Xi’an University of Science and Technology，2012，32(5):581 -585.

［6］ 李 翔，王金安，张少杰.复杂地质体三维数值建模方法研究［J］.西安科技大学学报，2012，32(6):676-681.LI Xiang，WANG Jin-an，ZHANG Shao-jie.3D modeling method of complicated geological body［J］. Journal of Xi’an University of Science and Technology，2012，32(6):676 -681.

［7］ 计 宏，余学义.厚黄土覆盖层采动破坏数值模拟研究［J］.西安科技大学学报，2012，32(2):209 -213.JI Hong，YU Xue-yi. Numerical simulation of surface subsidence in thick loess mining seam［J］. Journal of Xi’an University of Science and Technology，2012，32(2):209 -213.

［8］ 白如红，郝君，何龙伟，等.浅埋煤层综采工作面矿压显现规律数值研究［J］.西安科技大学学报，2012，32(2):676 -681.BAI Ru-hong，HAO Jun，HE Long-wei，et al.Numerical study of mine pressure regularity in fully-mechanized face of shallow coal seam［J］.Journal of Xi’an University of Science and Technology，2012，32(2):676 -681.

［9］ 李志强，唐 旭.水射流卸压增渗及抽采瓦斯效果的渗流力学数值解［J］. 西安科技大学学报，2012，32(4):464 -469.LI Zhi-qiang，TANG Xu. Stress relief and permeability increasing by water jet and numerical method of seepage mechanics equation for gas extraction effect［J］.Journal of Xi’an University of Science and Technology，2012，32(4):464 -469.

［10］雷志栋，杨诗秀，谢森传. 土壤水动力学［M］. 北京:清华大学出版社，1988.LEI Zhi-dong，YANG Shi-xiu，XIE Sen-chuan.Soil water dynamics［M］.Beijing:Tsinghua University Press，1988.

［11］唐丽云，杨更社，赵 涛，等.重力坝坝体水平变形模型研究［J］.西安科技大学学报，2012，32(2):193 -197.TANG Li-yun，YANG Geng-she，ZHAO Tao，et al.Horizontal deformation model of gravity dam［J］. Journal of Xi’an University of Science and Technology，2012，32(2):193 -197.

［12］王松阳，李树刚，王红胜.水渗流对薄层复合顶板的弱化影响［J］.西安科技大学学报，2012，32(6):686-690.WANG Song-yang，LI Shu-gang，WANG Hong-sheng.Weaken influence of water seepage on thin-layer compound roof［J］. Journal of Xi’an University of Science and Technology，2012，32(6):686 -690.

［13］ 周 明.利用地中渗透仪观测资料进行降雨入渗补给规律分析［J］. 地质科技情报，2002，21(1):37 -40.ZHOU Ming.Analysis of precipitation recharge using observed data of lysimeter［J］. Geological Science and Technology Information，2002，21(1):37 -40.

［14］赵贵章.鄂尔多斯盆地风沙滩地区包气带水——地下水转化机理研究［D］.西安:长安大学，2011.ZHAO Gui-zhang.Study on transformation mechanism of vadose zone water-groundwater in the wind-blown sand area of the Ordos Basin［D］.Xi’an:Chang’an University，2011.

［15］冯 平，李建柱.土壤水分下渗机制及其在半干旱区产流模拟中的应用［J］.干旱区资源与环境，2008，22(6):95 -98.FENG Ping，LI Jian-zhu. Ifiltration mechanism of soil moisture and its application in run off generation modelling［J］. Journal of Arid Land Resources and Environment，2008，22(6):95 -98.

［16］杜涛涛，窦林名，蓝 航.定向水力致裂防冲原理数值模拟研究［J］. 西安科技大学学报，2012，32(4):444 -449.DU Tao-tao，DOU Lin-ming，LAN Hang. Simulation on rockburst prevention by directional hydraulic fracture［J］.Journal of Xi’an University of Science and Technology，2012，32(4):444 -449.