降雨入渗补给地下水数值模拟研究

大致在 0.18～0.27，分析认为影响降雨入渗强度的因素有降雨量、潜水埋深、包气带岩性等。其中在研究区广泛分
布的风积沙对地下水资源起到了一定的保护作用。
关键词:包气带；数值模拟；降雨入渗；土壤水力学参数；神东矿区
中图分类号：P641.131；O241
文献标识码：A
Study on Rainfall Infiltration Recharge Numerical Simulation
7期
何 渊：降雨入渗补给地下水数值模拟研究
37
横交错，沟谷陡峻狭窄，植被稀少，水土流失严重。 项目组于 2013 年 8 月在研究区由南至北垂直乌
兰木伦河按地形地貌特征布置了三条剖面线进行包 气带岩性调查，同时取样。取样范围基本控制研究 区域，每条剖面布置 4 个取样点，沿剖面布置的每个 取样点取包气带埋深 0～3m 颗粒分析样品，同时在 示范矿井（补连塔与大柳塔煤矿）范围及包气带岩性 （风积沙和黄土状图）典型区域取 0～3m 包气带原状 样，每个深度原状样均取双样(图 1)。
气象要素 平均数
标准差 变异系数 变异性评价
降雨量 366.0mm
88.5ห้องสมุดไป่ตู้
0.3
弱
蒸发量 2221.3mm
178.8
0.1
弱
无降水日
342d
3.9
0.04
弱
通过对研究区气象站的降水、蒸发和无降水日 的变异性分析可见，区内各气象要素均为弱变异。 由于区内各气象要素均为弱变异，且变异系数较小， 因此该频率年的各气象要素具有代表性。 3.2 水文地质概念模型
2012 年）的气象资料，气象要素包括蒸发量、降雨
量、无降水日，各要素均精确到月。
∑ 对研究区气象数据取平均，即
-x
=
1 n
n i=1
xi
，公式
中 x 表示某一个气象要素的测量值，数目为 n 个，-x
表示某一个气象要素的 n 个测量值的平均值。记气
象 要 素 的 标 准 差 为 σ，该 值 计 算 采 用 公 式 ：
①原状土样土壤水分特征曲线测定。非饱和土 壤中的基质吸力或压力水头和含水率之间有一定的 对应关系，这种关系称为水分特征曲线 。 [5] 非饱和 渗透系数可根据水分特征曲线进行推定。土壤水分 特征曲线成为非饱和渗流数值模拟的重要基础。
此次土壤水分特征曲线测定试验采用日本 hi⁃ mac 公司生产的 CR21G 离心机进行测定，测定不同 压力下土壤体积含水量，绘制成曲线如图 4 所示：
0 引言
随着神东矿区煤炭资源开发规模逐渐增大，水 资源短缺问题日趋严重，实现地下水资源可持续开 发利用，为能源基地建设与开发提供水资源保障已 迫在眉睫。而实现此目标的前提是准确计算和评价 现有水资源。研究区内地下水循环方式主要以垂向 交换为主。作为地下水资源评价中补给量的主要构 成部分，以往的评价大多采用经验估算方法,这样概 化存在较大误差。此次研究是在广泛调查神东矿区 包气带岩性结构并结合野外取样、室内参数测定的 基础上，采用数值模拟方法建立包气带水分运移数 值模拟模型，通过模拟计算得出针对神东矿区内不 同地段降雨入渗强度及入渗系数，并对影响降雨入 渗强度的主要因素进行探讨。研究中获得的土壤水 力学参数可用于同类型区域的水文地质计算。
②原状土样土壤饱和导水率测定。土壤饱和导 水率测定采用环刀法，用环刀取原状土样，浸水后， 在单位水压梯度下，根据达西定律，求得通过垂直于 水流方向的单位土壤截面积的水流速度，称为土壤 的饱和导水率或渗透系数。测得结果见表 2：
③颗粒分析样品土壤水力学参数获取。颗粒分 析采用世界通用的马尔文激光颗粒分析仪测定(型
何渊
（中煤科工集团西安研究院有限公司，陕西 西安 710054）
摘 要: 针对在地下水资源评价中入渗补给量采用经验估算存在较大误差的问题，以神东矿区为例，在广泛调查矿
区包气带岩性结构并结合野外取样、室内参数测定的基础上，采用数值模拟的方法建立包气带水分运移数值模拟
模型，定量模拟矿区内不同地段降雨入渗强度，探讨影响降雨入渗强度的主要因素，计算得出研究区降雨入渗系数
基金项目：国家科技支撑项目（2012BAC10B00）。 作者简介：何渊（1980—），女，山西太原人，硕士，工程师，主要从事
水资源环境及矿井水害防治方面的研究工作。 收稿日期：2014-05-17 责任编辑：樊小舟
1 研究区概况
神东矿区地处毛乌素沙漠的东南缘、陕北黄土 高原北部。地形总的趋势为中部高而东西低，中西 部地形平坦。地貌以风沙滩地及半固定沙丘为主， 该地貌类型在神东矿区主力矿井中以补连塔煤矿最 为典型。矿区东部地形复杂，为黄土高原丘陵沟壑 区, 梁峁起伏, 沟道密集，地形破碎，此种地貌类型仅 在矿区的东部局部地区分布，以大柳塔煤矿最为典 型。神东矿区所在区域属干旱半干旱温带大陆性气 候, 降雨稀少，蒸发强烈，生态环境极其脆弱。
边 界 条 件 ：模 型 的 上 边 界 条 件 由 气 象 条 件 决 定，HYDRUS-1D 软件代码根据连续的气象数据采 用能量平衡方程确定地表的能量通量。同时，计算
得到的地表蒸发作为水分运移的上边界条件。模型 下边界考虑为自由排泄边界(图 2)。 3.3.3 数学模型的求解
上述数学模型求解采用美国盐土实验室开发的 hydrus-1D 软件。数值离散采用有限元法，由程序 自动完成。长度、时间单位分别选用 cm、d。初始时 间步长为 0.01d，计算过程中，根据计算收敛的情况， 时 间 步 长 由 程 序 自 动 调 整 ，最 小 时 间 步 长 为 0.000 01d，最大时间步长为 1d[7]。 3.4 参数获取 3.4.1 土样参数的实验室测定
为：C≤0.1 属弱变异；0.1＜C≤1 属中等变异；1＜C 属 强变异[1]。
按上述方法计算得到研究区气象站 95%降水频 率年，各气象要素的变异程度载于表 1 中。
表 1 研究气象站区 95%频率年气象要素变异性分析
Table 1 Variability analysis of 95% frequency yearly meteorological factors of weather station spot in study area
图 1 包气带调查实际材料图
Figure 1 Map of original investigation data in aeration zone
3 包气带水分运移数值模拟
3.1 研究区气象数据的统计分析
为了获得研究区域的气象要素特征，项目组收
集 了 研 究 区 气 象 站（伊 旗 气 象 站）近 20a（1993-
根据野外典型点取样的包气带岩性结构概化出 的水文地质概念模型如图 2 所示。由于所计算区域 包气带中水分运移主要是垂向运动，因此，将其概化 为垂向一维流。试验介质概化为 2 层，z 轴坐标原点 选在地表，向上为正。已有研究认为[4]，在包气带厚 度小于潜水蒸发极限深度的条件下,随着包气带厚 度的增大,降雨入渗强度逐渐减小。大于 3m 后,年降 雨入渗补给强度将趋于稳定，与埋深无明显相关关 系。根据本次研究收集区域内地质、水文地质资料 分析，神东矿区内地下水位埋深除沟谷地段外基本 均大于 3m，故本次数值模拟的包气带计算剖面厚度 概化为 3m。
∑ σ =
1 n
n
(xi
i=1
-
-x)2
。记变异性系数为 C，计算公式为
C = σ/-x ，按水资源中的常规处理对变异性系数分级
图 2 水文地质概念模型示意图
Figure 2 A schematic diagram of hydrogeological conceptual model
3.3 数学模型 3.3.1 控制方程
图 4 原状样土壤水分特征曲线
Figure 4 Undisturbed sample soil moisture characteristic curves
7期
何 渊：降雨入渗补给地下水数值模拟研究
39
表 2 原状样土壤饱和导水率测定结果
Table 2 Undisturbed sample soil saturated hydraulic conductivity tested results
第 26 卷 7 期 2014 年 7 月 doi:10.3969/j.issn.1674-1803.2014.07.09 文章编号:1674-1803(2014)07-0036-03
中国煤炭地质
COAL GEOLOGY OF CHINA
Vol.26 No.07 Jul. 2014
降雨入渗补给地下水数值模拟研究
2 包气带结构特征
根据野外调查及结合前人已做过的地质地貌 图，遥感解译图等资料分析[1-3]，研究区在地貌单元 上可划分为风沙地貌和黄土地貌两大类。以风沙地 貌（土壤质地分类为砂土、壤质砂土）为主，覆盖区内 90％以上面积，其地表形态表现为风沙滩地、沙丘及 沙梁，在风沙滩地区，湖泊海子星罗棋布；黄土地貌 （土壤质地分类为砂质壤土、壤土为主）仅在区内东 南部沿沟壑地带分布，地形切割较深，冲沟发育，纵
He Yuan
(Xian Research Institute, China Coal Technology & Engineering Group Corp, Xian, Shaanxi 710054)
Abstract: In allusion to the inaccuracy existed in empirical estimation of the recharge from rainfall infiltration, A research had been car⁃ ried out in Shendong mining area. Based on broad inquiry of aeration zone lithologic structure in the mining area, combined with field sampling, laboratory parameter measurement carried out aeration zone moisture migration numerical simulation modeling. Quantitative⁃ ly simulated rainfall infiltration intensity, discussed main factors impacting the intensity, estimated rainfall infiltration coefficient ap⁃ proximately at 0.18~0.27 in the study area. Analysis has considered that factors impacting rainfall infiltration intensity have rainfall rate, phreatic water buried depth, aeration zone lithology etc. Furthermore, the eolian sand is widely distributed in the study area has played a certain protective effect for groundwater resources. Keywords: aeration zone; numerical simulation; rainfall infiltration; hydraulic parameters of soil; Shendong mining area
为原点，向上为正；
3.3.2 定解条件
初始条件：模型的初始条件设置为计算剖面的
负压值，假定地表为负压值为 0m，计算剖面底部
（3m）处负压值为-3m，线性分布(图 3)。
图 3 计算剖面初始负压分布
Figure 3 Initial negative pressure distribution on rated section
包气带水分运移基本方程[5-6]：
38
中国煤炭地质
第 26 卷
∂θ ∂t
=
∂ ∂z
éêK ë
(
h,
z
)æè
∂h ∂z
+
1öøùûú
,
式中：θ 为土壤体积含水率，cm3 ·（cm）-3；K 为
非饱和渗透系数，cm·h-1；t 为时间变量，h；z 为空间
变量，cm；h 为压力水头，在包气带中称为负压；地表