基于质点运移模型的辽河流域区域地下水循环规律研究

李 昕

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

辽河流域是我国东北地区农业发展以及粮食生产的重要基地，其地表水资源由于水体污染严重、水资源供需关系紧张以及水资源浪费等多种因素的综合作用而日趋严峻，地下水作为该区域主要供水水源之一对于缓解供水压力、促进区域发展发挥重要作用。对辽河流域地下水循环规律进行研究对于促进区域水资源可持续发展以及科学规划具有重要意义。近年来，国内外学者针对水资源循环规律分别采用多个理论和方法从多个角度开展了大量研究，如韩冬梅等[1]采用数值模拟技术和地下水流模型对忻州第四系孔隙地下水流平均滞留时间、分布特征进行了分析模拟；孙跃等[2]依据黄河分布区域水文地质对河南省供水安全问题采用数值模拟技术进行研究，并运用局部地下水流预测模型对10个水源地的地下水开采量进行评价分析；王丹等[3]利用地下水动力学法对草甸区剖面地下水进行二维数值模拟，对地下水运动规律及影响因素进行揭示和表征；蒋小伟等[4]对地下水模型通过地下水年龄实测法进行矫正，研究表明地下水年龄从补给区到排放区整体呈现变老趋势；地下水年龄在盆地下游垂直方向可产生突变，并以此可对不同级次流动系统进行评判。对辽河流域地下水循环规律研究相对较多且主要是采用同位素追踪法和水化学法进行探讨，而对水动力场的相关研究主要侧重于利用定性分析，如刘玉等[5]对陇东盆地泾河流域地下水循环规律采用水文地球化学理论和同位素基础进行研究；徐中华等[6]以陇东盆地白垩系中层地下水为研究对象，采用水文地球化学理论对其水化学演化规律和发展进行了深入的研究分析；张彦林等[7]对陇东盆地地下水运行状况采用同位素法和水化学法进行探讨并对其运移规律以及影响因素进行了探讨。

文章采用数值模拟技术对辽河流域流速分布特征、水循环宽度、深度以及地下水运行更新时间进行定量的分析和计算，以期为该流域水资源的科学规划管理以及合理开发利用提供一定的决策依据和理论支持。

1 区域概况

辽河流域位于我国东北地区西南部，属于温带半湿润半干旱季风气候，冬季寒冷干燥漫长、夏季炎热多雨。全年降雨量和径流量时空分布极不均匀，降雨量由东南向西北方向整体呈降低趋势，多年平均降雨量为320～860mm，其中62%以上集中在7—8月份，并且多以暴雨或强降雨的形式出现；气温在平原地区较高，山地较低，由南向北整体呈降低趋势，多年平均气温为4～9℃；由南向北方向年蒸发量依次递减，并处于982～1650mm范围。流域内地形支离破碎、沟壑纵横且具有强烈的切割特征，因此为典型的低山丘陵地形，区域内海拔高度由西北向中南部逐渐降低，最高和最低处约为1650m和860m。流域内各支流纵横交错、蜿蜒曲折，主要支流有老哈河、浑河、太子河、柴河等，研究区域地质结构为白垩系，虽然在流域内白垩系分布较为广泛，但在河谷与沟谷底部或河谷底两侧多以树枝状显现，受天环向斜控制，盆地结构多以雁行排列的向斜组成。可开采利用的主要含水层位于白垩系含水岩组，而太子河底部含水层岩组和白垩系顶部的砂岩为风成砂岩，是主要含水层[8]。浑河底部多以泥岩为主，是该流域的隔水层或弱透水层。地下水多以微淡水、淡水为主，并且受补给条件、含水层厚度以及区域分布特征作用，该区域地下水可开采潜力整体处于较低水平。

2 研究方法及内容

对辽河流域剖面地下水利用Visual Modflow软件进行数值模拟，并且利用向前追踪例子算法和Modpath模块生成剖面地下水流线分布，进而建立辽河流域地下水动力场。对辽河流域剖面水头特征利用地下水动力场进行确定，对其地下水流系统采用Toth地下水流系统理论进行划分并确定各系统的循环特征和尺度特征，然后对地下水的平均滞留时间利用Modpath模块进行计算。

2.1 质点运移数学模型

Modpath是一种可对任意时刻水质点的位置和流线分布进行简便、快捷的计算和确定的模块，并可对质点向前、向后示踪流线进行准确的模拟预测。Madpath的质量平衡方程可对地下水二维稳定流动的渗流流速和有效孔隙率进行求解，计算公式如下：

(1)

式中，Vx，Vz—分别为流动流速矢量在坐标轴方向的分量，LT-1；w—含水层内单位体积水源汇流产生的水量，T-1；n—含水层有效孔隙率，%。

2.2 计算流速

采用单元格的形式对地下水流速度进行计算，即每个单元格对应一个水流速度，计算公式如下：

(2)

式中，Vx，Vz—分别为水平和垂直方向的水流流速，LT-1；n—含水层有效孔隙率，%；Qx，Qz—分别为通过水平和垂直方向的水量，L3T-1；Δx，Δy，Δz—分别为单元格在坐标轴方向的长度，L。

2.3 径流划分标准

首先采用Visual Modflow软件并结合每个单元格的径流速度和剖面地下水流数值模型提出径流划分相关标准，然后对计算结果利用surfer10进行差分计算，并完成地下水流速矢量特征的输出，将其采用色阶的模式进行表征。考虑到辽河流域水流速度基本特征和实际状况，因此将其径流划分为弱径流、中等径流和强径流带[9]，所对应的径流速度区间分别为小于0.006m/d、0.06～0.12m/d、大于0.12m/d。

2.4 平均滞留时间

在不同含水层路径上地下水沿流线方向的平均流速的大小可采用地下水的平均滞留时间进行表征，并采用Modpath放置粒子法即通过对各含水层岩组运动路径上粒子所消耗的时间进行确定，进而可对地下水的平均滞留时间进行更加直观、科学的反映。

3 数值模型

3.1 剖面概况及概念模型

剖面贯穿辽河流域的东西方向并且在剖面上共布设18处钻孔，其中用于环河含水岩组和含水层研究分析的钻孔分别为14个和4个。区域等水位线分布方向整体上与剖面走向呈垂直分布特征，并且白垩系含水岩组厚度在垂直方向上存在显著差异，并且含水岩组整体上在西部较厚而东部较薄。剖面西部区域为降雨入渗主要补给区域，而剖面东西两侧为主要径流补给区域，河道或河道交汇处为水流系统的排泄基准面并且是地下水排泄的主要排泄区域[10]。

罗汉洞含水岩组顶板为辽河流域剖面的上边界，在浑河与太子河交汇处尖灭，并且在上边界处出现降雨入渗补给，因此该处为流量边界。侧向径流补给处往往为剖面两侧，因此可称为定流量边界条件，洛河组含水岩组为流域剖面底部，因此，可称为隔水边界。侧向径流补给和大气降雨入渗补给为剖面补给的主要来源，地下水向河流的排泄为排泄相的主要来源，采用定水头的方式对辽河与其他支流交汇处以及沿河道进行处理。含水层的分布特征采用辽河流域剖面水文地质进行表征，依据相关文献和研究成果对水文地质剖面含水层组的岩性可概化为3层，即洛河、环河以及罗汉洞含水岩组，可对罗汉洞组、环河以及洛河含水岩组承压水分布区概化为3个、4个和6个渗透系数分区，各渗透系数K1～K13取值为0.005～1.2。通过构建数学模型进而使得水文地质趋于客观系统的为概念模型的近似处理本质，对自然条件下剖面上地下水循环基本特征进行客观、准确的解释为建立模型的首要目标，所以，研究采用非均质各向异性稳定流对研究区域模型进行概化，以此提高模型与实际状况的贴近程度[11]。

3.2 地下水流数学模型及求解

可采用下述微分方程对水文地质概念模型的定量问题进行求解：

(3)

式中，D，h—分别为剖面渗流区域和含水层水位标高；Kx，Kz—分别为水平方向和垂直方向的渗透系数，m/d；h1，q(x，z，h)—分别为水头边界上水头分布和流量边界单宽流量；Γ1，Γ2—分别为水头边界和流量边界。

利用Visula Modflow软件对模型进行求解，其基本过程为在时间和空间上将研究区域进行离散化处理，并对每个网格点构建均衡方程并联立形成一个大型的线性方程组，然后对每个网格的水头值进行迭代求解。文章将辽河流域剖面划分为25层、180列和1行的等面积单元格并以此更加准确、客观地揭示该区域的地下水流特征。采用与条件共轭梯度法对水流数学模型进行求解，其中内部和外部迭代次数分别为50次，水位变化和残差收敛标准均为0.01m。

依据水化学场特征对模型的水动力场运行结果进行验证分析，其中剖面上的水化学场采用剖面附近12处水井并进行水化学数据分析，并按取样点位置将所取水样的水化学类型进行放置。结果显示，地下水阴离子在太子河东部区域类型多变，而且阳离子在排泄区与过渡区变化不大，并以钠离子为主，水化学特征在整个研究区域与地下水运移规律基本保持一致。

4 结果与分析

4.1 水头分布特征

剖面模型水头分布计算结果显示辽河流域南部平泉县附近的水头为1200～1280m，水头整体上由南向北方向逐渐降低。该区域的水头分布特征为水流由南向北如平泉县、宁城县、开鲁县和通辽市等逐渐呈降低趋势，该区域属于辽河的西部发源地，共有两处水源即老哈河和西拉木伦河水源。

4.2 系统划分

结合相关文献资料将辽河流域剖面依据地下水流的源、汇关系以及实际状况可划分为区域循环、中间循环和局部循环三级循环系统，其中整个剖面内的地下水为区域循环系统，而在河道、湖泊附近多以中间和局部循环系统存在。

4.3 尺度特征

由于河流支流众多，在辽河流域剖面上河流与地下水存在水流交换并形成了多个地下水局部循环系统，此类循环系统延伸区间通常位于5～40km，其发育深度为20～180m。在老哈河、太子河、西辽河等区域的局部地下水循环系统其宽度通常低于10km，其发育深度约为500～800m；其他主要支流局部循环系统宽度为10～80km，发育深度相对较浅，为60～500m。辽河流域北部的碧流河、少冷河、查干木伦河区间发育循环系统，该区域循环系统宽度最大为100km，其发育深度可达到洛河岩组底部最大深度约1200m，各循环系统的发育宽度范围见表1。

表1 各级水流系统发育范围

4.4 水循环特征

在辽河流域剖面上，侧向径流补给为承压水的主要补给来源，而大气降雨入渗补给为流域南部河段的补给源，而北沙河、兰河、细河作为地下水排泄的主要区域。分别选取2个过水断面进行研究分析，结果表明地下水流线基本与过水断面相垂直，而断面水平径流在垂直方向流动不显著。所以，可按照剖面循环系统对两个过水断面经过的模型层进行分析，对各单元格的水平流量进行计算和求解，结果见表2。

表2 辽河剖面三级循环系统流量分析结果

研究表明，区域和中间循环系统为地下水循环的主要集中区域，并且在Z2过水断面上局部和中间循环系统的流量占比处于较低水平，而绝大多数系统流量为区域循环系统，其原因主要为洛河组岩层可对Z2断面层的东侧进行补给，而砂岩的环河组岩层的补给量相对较低，基本可忽略；结合区域实际状况，环河组区域主要为中间循环系统，而洛河组岩层位于区域循环系统。因此，造成区域循环系统流量相对中间和局部循环系统处于较高水平。

4.5 水流分布特征

辽河剖面地下水流平均速度为0.0021m/d，其变化区间为0.00006～0.012m/d。研究表明：辽河流域的西部区域属于弱径流带，其水流速度相对较低；浑河区域地下水流速有一定程度的提升并在洛河组岩层区达到最大，该区域属于强径流和中等径流交界处，流速在东辽河与西辽河区域逐渐降低为弱径流带。整体上看，上部环河组和下部洛河组分别位于弱径流带和中等径流带。

4.6 平均滞留时间

研究采用放置粒子法并结合流速分布特征一级含水层岩组分布状况对水流运移路径进行模拟，其中地下水贯穿了3个不同的含水岩层。研究表明，强径流带、中等径流带和弱径流带的循环更新时间分别为小于4218a，4218～25162a，大于25162a。

5 结论

剖面模型水头分布计算结果显示，辽河流域南部平泉县附近的水头为1200～1280m，水头整体上由南向北方向逐渐降低。按照区域特征可将辽河流域划分为区域循环、中间循环和局部循环三级循环系统，其中整个剖面内的地下水为区域循环系统，而在河道、湖泊附近多以中间和局部循环系统存在。

在老哈河、太子河、西辽河等区域的局部地下水循环系统其宽度通常低于10km，其发育深度约为500～800m；其他主要支流局部循环系统宽度为10～80km，发育深度相对较浅，为60～500m。区域和中间循环系统为地下水循环的主要集中区域，并且在Z2过水断面上局部和中间循环系统的流量占比处于较低水平，而绝大多数系统流量为区域循环系统，其原因主要为洛河组岩层可对Z2断面层的东侧进行补给，而砂岩的环河组岩层的补给量相对较低，基本可忽略。

辽河流域北部的碧流河、少冷河、查干木伦河区间发育循环系统，该区域循环系统宽度最大为100km，其发育深度可达到洛河岩组底部，最大深度约1200m；辽河剖面地下水流平均速度为0.0021m/d，其变化区间为0.00006～0.012m/d；强径流带、中等径流带和弱径流带的循环更新时间分别为小于4218a，4218～25162a，大于25162a。为提高地下水资源的开发利用效率和科学性，可对中等径流带的地下水进行合理开采利用，而对弱径流带禁止开采或少量开采。