基于Visual Modflow的太子河流域地下水流数值模拟分析

徐世俊

(辽宁省辽阳水文局，辽宁 辽阳 111000)

水是人们赖以生存的必要条件，日常生活用水以及农业灌溉用水更为重要。地下水在生活用水、工业用水以及农田灌溉用水等方面发挥着至关重要的作用[1]。为了预防水资源污染、过度开采水资源引起的严重缺水并合理化配置水资源，水流动态模拟分析一直被国内外学者认为是一种非常重要的数学模拟方法。该数值模拟软件是由美国著名学者采用三维有限差分法开发研制而成的，其适用范围为孔隙介质地下水。Visual Modflow模型软件是一种地下水模拟可视软件，由多种地下水模型软件的功能属性如Modflow、Modpath、MT3D、RT3D、WIN PEST等共同形成[2]。

模拟分析范围的确定受多个水文地质因素影响，如地形地貌、地质结构、岩层分布、天气环境、农业实用、地下水流情况等，其中地质条件属性受均质性、各向性、孔隙裂隙、流体密度等因素影响，地下水流情况又可分为稳定流态和非稳定流态水流、饱和流态和非饱和流态水流。软件模型由一维的点式模型逐渐发展成二维平面模型，随着科学计算机的进步和发展逐渐形成了现在适用性最好的三维立体模型。边界范围和初始条件也是软件进行模拟分析范围确定的主要参数。同时为了提高软件模拟的精确性和可适用性，还应考虑野外现场测试以及试验测量数据如弥散系数、渗透系数、分配系数、反应速率常数等[3]。

文章以太子河流域地下水为例，考虑以上对模拟范围确定的影响因素，分析该流域的水文地质构造条件、地下水流特性、农业用水以及生产生活用水等情况。通过分析太子河流域水文地质、岩层结构、水流动态方向、地理位置等条件，制定各参数取值范围，建立了Visual Modflow软件模型，然后对模型的计算结果与实际检测结果进行对比验证，并结合当地水土资源监管部门对地下水位的控制要求，为太子河流域地下水资源的开采和利用提出了合理化建议[4]。

1 太子河流域地下水状况

1.1 水文地质条件

处于我国辽宁省东部的太子河流域其大气环境属于温暖带湿润半湿润气候区，气候适宜，平均温度为5～10℃，全年降水量充足，其中71.2%降雨量集中在6～8月。太子河全长413km，流经本溪、辽阳以及鞍山等市，流域面积约为13880km2。太子河流域上游是以落叶阔叶林为主的山地森林富水区，植被所占面积超过50%，为山地森林多水区；低山丘陵区即为太子河中游地区，所占面积较少，为6.1%；其下游区为农作物平原区，也称为少水区，比较适合农作物生长及城镇居住，该区域面积适中，占流域总面积的24.9%。

在分析研究区域含水介质以及水流特性的基础上，将太子河流域地下水层合理地划分为三类即第四系松散岩类孔隙潜水、白坚系碎屑岩裂隙孔隙潜水及承压水。含水层受地质条件和沉积岩体变化影响，其分布规律和形成条件相差较大。综合考虑，太子河上游区为第四系松散岩类孔隙水，该层地质基岩以砂岩居多，其中含有少量泥质或钙质胶结的泥岩和碎岩，厚度约为200～350m。地表具有很厚的松散堆积物，为地下水的形成和存储提供了必要的先决条件。该区域地势平坦，降雨量较多，有利于对降雨进行存储补给，是该流域富水地区，有地下水库之称。中游区为白坚系碎屑岩裂隙，岩层厚度较小，约10～30m，最大不超过50m，上层含有少量泥沙，透水性较差。该区域分布狭窄，宽度较小，一般在800mm以内，低山丘陵地面高低起伏，土质结构疏松，虽有一定降水量但受自然地貌环境因素制约，对水分的吸收和储存能力较差，富水性差，不利于大气环境的水气循环。下游区域土地平缓，为承压水层，水层厚度20～80m，水位埋深5～20m，适于农作物生长及人们生活居住。该层地下水主要被用于植被生长灌溉、供应生活饮用水、工厂生产化工用水。据统计，该流域各方面总供水量高达53.382×108m3/a，其中农作物灌溉用水47.993×108m3/a，占用水量的90%，生活用水为4.382×108m3/a，工业用水2.074×108m3/a，林牧业用水5.558×108m3/a。

1.2 孔隙潜水运动特征及化学特性

由于含水层之间存在相对不稳定的含水层，故第三层含水层与太子河潜水层在相同部位的结构动态特征大致一样，且其地下水运动规律同样受天气气候影响。中游山地丘陵地区，其降水期主要集中在3～5月，地下水位在5月份开始逐渐上升至7月份达到极值点，然后开始下降，到第二年3月份下降到最低。在少水区即平原地区，地下水不仅与降水量有关，而且与农业灌溉以及生产生活用水量有关，同样也受气象气候影响[5]。降雨补给、河流运动以及日照蒸发是研究地区地表水汇入和流出的重要途径，占总径流量的92%以上。根据大气环境和年降水量，每年地表水流经径流量变化幅度不大，属于均衡分配。结合每年水流量变化情况，当年8月份至来年5月份是流量最大时期，流量在6月份会急剧下降，其原因可能是农牧业使用或灌溉用水。6月中旬至7月中旬流量将至最低，8月份开始逐渐上升至来年5月份。白坚系碎屑岩裂隙潜水存于低山丘陵区，该地区以风积碎石、胶泥砂以及吸粉砂为主，地表结构质地松散，对于降水和地下流经方面产生不利影响，不利于地下水的交流循环以及盐含量较低。因此在个别山丘低洼地区，地下水位较高，阳光蒸发扩散作用强烈，风化作用严重。本地区矿物化学成分主要为HCO3- SO4- Na- Ca，矿化度大于1200mg/L，溶解性固体271.25～393.84mg/L，小于600mg/L，水质情况一般。总硬度在172.35～198.86mg/L范围之内。pH值在7.25～8.18区间，属于弱碱性微硬型淡水。

2 建立Visual Modflow模型

2.1 概念模型与数学模型

2.1.1 概念模型

文章通过以上分析并结合实际情况，建立Visual Modflow模型，模型模拟的范围为太子河流域下游区域，区域属于含水层平原断裂性水文地质，其边界上没有侧向地下水流入，仅有侧向流出。含水层埋藏于盆地及坚硬的沉积岩下面，形成了第一承压含水层和第二含水层，地下水位因含水层边界变化所引起的流量变化量相对较小。下游区域地势平坦，从而形成了可约束的地质边界条件[6]。研究区域包括本溪、辽阳等所有的灌溉区域，其中包括少量个别的鞍山灌溉渠，覆盖面积近8600km2。下游地区为居民生活区，其地下水由北向东南方向运动，东南部的本溪地下水自东南向西北方向运动，最终均汇入太子河。太子河地下水源汇流包括：气候降雨渗流、生活用水排放、农业灌溉汇入、冰雪融化流入、地下深水供给、日照蒸发、植被吸收、人工开采利用等，其中气候降雨、人工开采利用以及农业灌溉是地下水汇流的主要影响因素。降水变化引用系数及降水汇流系数计算结果见表1。

表1 降水汇流系数计算表

通过大量的查阅资料并对水文地质进行现场勘探观测，将水文地质边界条件概化为以南部、西部以及北部作为流入边界，因为此边界水流汇入相对较小，结构层相对致密[7]。东部为地下水流流出边界，东北部可视为零流量边界。不同含水层之间存在的水头差相互作用，产生不同的影响是将其作为一个水系统的关键因素。研究区内存在断裂层，断裂层与含水层之间相互贯通、相互作用，因此将两者概化为一个含水层进行模拟分析。由于Visual Modflow无法在一个界面上做一对多层编辑，故将南区含水层划分为两层，孔隙潜水层和承压含水层之间是存在密切关系的含水层。

2.1.2 数学模型

文章在充分考虑了研究区域边界概念模型后，将太子河流域地下水流按照一定的顺序逐级划分为潜水层、弱透水层以及深水层。为了简化计算，减少计算工作量，将水文地质条件进行了简化和概化，并对不同含水层的参数进行了分区。根据软件的计算原理和方法，并采用数学计算模型，则简化后的地下水流数学方程的建立和求解过程如式(1)所示。

(1)

式中，h—水头；K—渗透系数；*μ—含水层的单位储水系数；G—辅助计算项；W—地下水流排泄项；h0—初始水头；h1(x，y，z，t)—第一类边界上的水头；f1—第二类边界上的水流通量；n—边界外法线方向；Γ1，Γ2—第一、第二类边界。

(2)

式中，n—有效孔隙度；C(x，y，z，t)—污染浓度；Ce—污染源的浓度；D—弥散系数张量。

上述公式可用如下矩阵表示：

(3)

Visual Modflow软件是通过将含水层模拟成像形成一个三维的网格图形，将研究区域的地下含水层按照一定的参数设置并结合地质条件划分为若干个小的含水层，每个小的含水层又根据坐标和数值大小划分为行和列。因此含水层通过参数设置计算被划分为多个小长方体，其中心位置被称为格点，格点的多少代表被划分的多少。然后将划分后的格点通过有限元差值分析法求解地下水流方程，求解方程并建立该区域的地下水模型。为更好地标记和计算地下水流方程，用字母i代表格点所在位置的行号，字母j代表所在位置的列号，长方体的的厚度用k表示。研究规定自上而下第一层为模型的顶层即k=1，随深度的降低k值逐渐增大，并用x和y方向代表行和列的方向[8]。概念模型如图1所示。

图1 Visual Modflow含水层三维单元格划分

2.2 模型的识别与验证

2.2.1 研究区域的剖分

根据太子河流域的真实情况，并考虑了该地区的复杂性，文章建立的太子河流域下游地下水水文地质概念模型将地下水流视为三维非稳定流态，并将其地下含水层概化为单层非均质各向异性潜水含水层[9]。采用有限差分求解，并采用强隐式法联立迭代求解代数方程组，在垂直空间上，采用X、Y正交网格对模拟区域进行平面上的等间距划分，△X=△Y=1km。将太子河流域下游平原区划分为210×180的矩形网格单元，因此位于计算区内的单元有3614个在垂直方向上，考虑含水层结构特性和概化措施，将含水层划分为2层，每层有效计算单元为3076个，所以该研究区域内有效单元计算格点总共为7228个。

2.2.2 计算的时限与步长

受研究区域内水文地质历史资料完整性、人工开采数据真实性以及统计科学的发展局限性，并考虑本模型研究的使用性范围和精确性，模型模拟期选择2000年1月～2005年12月，共72个月，采用简介识别法验证模型中渗透参数和弹性储水率。为降低模型模拟与实测水头之间的误差，文章选择地下水流变化和降水较小的时间段作为模型识别期和验证期，故将模型识别期定为2005～2007年，模型验证区间为2008年。因地下水开采利用对水位影响较大，故考虑了地下水开采对水位变化影响，并以每季度视为一个地下水开采期进行模型参数输入。

2.2.3 初始压力水头的确定

对格点初始压力水头的确定可通过计算含水层平均水头减去相应格点的上下标高求得，不同格点相应的平均总水头值是通过模拟含水层与隔水层在格点处厚度的方法求得。计算含水层地下水2005年的平均水头可参考不同含水层所在水位观测孔的数值[10]。

2.2.4 Visual Modflow模型识别结果

为了达到模型的识别目的，文章采用间接法求各参数值。通过试错法给参数设定变化范围，并不断调整，将拟合计算的水头和实测水头进行对比分析，最终以拟合结果误差最小值作为参数初始值。分区所计算的第Ⅰ、Ⅱ承压水层共含有25个参数区间，弱水层含有10个参数区间。通过识别、验证以及计算，本次模拟潜水层共有7个参数，其中水平方向的渗透系数区间为5～135m/d，垂直方向渗透系数区间为0.03～15m/d，模拟计算弹性储水系数区间为0.21～0.28。弱水层横向渗透系数区间为0.03～1.8m/d，弱水层竖向渗透系数区间为0.01～0.12m/d，弹性区间系数约为0.18。Visual Modflow模型对实测点进行模拟并拟合，模拟计算水位与实测水位拟合曲线如图2、图3所示。

图2 模型第Ⅰ承压水层观测孔水位拟合图

图3 模型第Ⅱ承压水层观测孔水位拟合图

为更好地分析拟合情况和实际情况之间的误差，先将统误差统计结果进行汇总分析，见表2。

表2 模拟水头与实测水头误差绝对值统计结果

3 地下水资源预测

采取不同抽水量测试模拟方法是对太子河流域地下水进行数值模拟时减少模拟结果偏差的重要方法。因此文章建立的Visual Modflow模型是通过对不同抽水量的不同测试孔结果进行分析计算该含水层的含水量变化，以避免过度抽取地下水造成抽水量超出实际资源量。同时采用线性回归方程分解法对太子河流域边界进行预测，考虑到该河流特征基本保持不变，所以对其他参数不进行设置，仅考虑河流水位变化情况。一般情况下，地下水位仅受降雨量、流经方向等自然环境影响，故不做更加深入复杂的参数设定。综合考虑以上因素，结合当前对太子河开采利用现状对研究区域地下水位进行预报，其各承压水层地下水位变化情况见表3。

表3 模型预测各承压含水层地下水位变化情况 单位：m

由表3可知，在目前地下水开采利用条件下，各承压水层最低水位逐渐降低，地下储水量逐渐减少，如果不能进一步加强用水控制管理并采取合理规划措施，未来太子河流域地下水将面临严重缺水现象。

4 结论

文章通过对Visual Modflow软件的研究应用，建立了适用于太子河流域的地下水流数值预测三维模型，结果表明该模型能够准确反映太子河流域的水文地质条件和特征；能够客观反映本地区实际情况，可以对该流域地下水流变化情况进行有效模拟。该软件操作简便、计算简单、功能性强、立体感好、适用性强，可对其他流域进行推广应用，能够及时准确地对地下水流进行预测。

文章通过Visual Modflow软件对太子河流域地下水流的发展变化情况进行预测，结果显示在各承压含水层水位逐渐降低，提出应及时采取有效措施加强对地下水的开采和管控。

[1] 涂亮, 宋汉周. 基于Visual Modflow的黑河中游地下水流数值模拟[J] . 勘察科学技术, 2009(02): 19- 23.

[2] 党素珍, 刘昌明, 王中根, 等. 近10年黑河流域上游积雪时空分布特征及变化趋势[J]. 资源科学, 2012(08): 1574- 1580.

[3] 邵景力, 崔雅莉, 赵云章, 等. 黄河下游影响带(河南段)三维地下水流数值模拟模型及其应用[J]. 吉林大学学报(地球科学版), 2003, 33(01): 51- 55.

[4] 王巍巍. 锦州市地下水水质现状与变化趋势研究[J]. 水利技术监督, 2016(05): 59- 62.

[5] 王志民, 于福兰, 刘志芬, 等. 阜新市地表水资源质量评价[J]. 东北水利水电, 2003(09): 38- 39.

[6] 高素丽. 辽阳市水资源开发利用和管理保护对策[J]. 水土保持应用技术, 2011(04): 45- 47.

[7] 杨锦华. 佛寺水库水源饮水安全与水环境保护[J]. 水利技术监督, 2006(03): 50- 52.

[8] 温树影. 大伙房水库水体富营养化现状分析及对策研究[J]. 水土保持应用技术, 2015(02): 31- 32.

[9] 马宇, 王淑伟. 辽宁省水资源现状分析及保护措施研究[J]. 水利规划与设计, 2015(11): 42- 44.

[10] 徐斌, 刘克岩, 时晓飞, 等. 地表水资源水质水量结合评价方法研究[J]. 水利规划与设计, 2003(02): 25- 30.