基于Visual MODFLOW的黑龙江省巴彦县地下水流数值模拟

娄安宇，陈 末，臧珊珊，吴志刚

(黑龙江大学 a.水利电力学院; b.寒区地下水研究所, 哈尔滨 150080)

0 引 言

随着社会的发展，巴彦县的城镇化建设也迅速发展，地表水资源已经供不应求，地下水资源的利用在该地区已达到不可替代的地位。不过，地下水不是无穷无尽的，所以要对该地区的地下水进行分析研究。地下空间错综复杂，地下水含水层中有着各种形状复杂的空隙[1]。因此，分析研究区地下水的补给来源、排泄途径，进行边界条件的概化，构建研究区地下水流概念模型，从而得到相关的水文地质参数是一个极为简单的方法[2]。

数值模拟方法在地下水研究方面的效率、灵活性和低成本管理等方面已成为一种重要的方法，并日益引起注意和广泛的应用[3]。数值模拟方法是一种数值模拟程序的微分方程，根据边界和初始基本条件，来描述地下水水位随着时间的推移而产生连续变动得数学模型,并且通过电脑软件中的离散化方法来解决这个数学问题[4]。

本文应用数值模拟方法，通过使用Visual MODFLOW 软件建立巴彦县地下水流数值模型，得到该地区相应的水文地质参数，并对这一地区地下水水位进行预测。

1 研究区概况

巴彦县位于黑龙江省西部偏南，隶属于哈尔滨市，巴彦县境南北长85 km、东西宽72.7 km，幅员3 137.7 km2。境内主要有“一江三河”，其中“一江”是松花江，“三河”分别为泥河、少陵河和五岳河。见图1。

图1 研究区示意图

巴彦县地区属于中温带大陆性季风气候，每个季节的温差很大，温度最高可以达到36.5℃，最低可以到-40.2℃，年平均气温约为2.9℃[5]。地下水开采量达12 734.19 m3，主要为农业灌溉用水，部分用于乡镇生活用水、工业生产和生态环境改善等方面。全县地下水资源开发利用量约为7 400×104m3[6]。

地下水贮存条件及分布规律为，一级阶地区含水层分布均匀，厚度稳定，堆积物松散，水量较丰富，地下水类型为阶地松散层孔隙承压水。漫滩区孔隙潜水在2月下旬至3月初为枯水期，进入雨季后水位逐渐上升，至8-9月初上升到最高值，水位年变幅1～4.5 m。阶地区孔隙承压水在4月下旬至5月上旬为枯水期，9月末至10月初出现最高值，地下水动态类型为降水入渗-径流型[7]。

综合全区地貌成因和形态，将巴彦县划分为3种成因类型，3个形态单元见表1。

表1 地貌单元分区表

2 地下水流数值模型

2.1 水文地质概念模型

以相对完整的水文地质单元为数值模拟计算区，尽量将计算区边界设置在自然边界处，或者设置在容易确定流量或地下水位的人为边界处。本次研究区位于哈尔滨市巴彦县，北部边界为泥河，西部边界为漂河和少陵河，南至松花江干流。位于松嫩平原东部，松花江中游左岸，依据地形地貌特征、地层结构及水文地质特征，将区内地下水划分为河漫滩松散层孔隙潜水、一级阶地松散层孔隙承压水和高平原(松花江二级阶地)松散层孔隙承压水。随着系统的变化，三维空间参数向量的一般地下水流动速度反映了输入和输出系统的异质性，输入和输出系统随时间和空间的变化而变化，来显示地下水的流动情况。

2.2 含水层系统的概化

该研究区的地下水主要是以水平运动为水流的运动方向，而垂直方向的运动则很少，水文地质参数伴随空间分布的变化而改变，这反映了地下水系统的均质性。地下水系统输入和输出随时间和空间的变化而变化，地下水流是不稳定的[8]。综上所述，将该模拟区的含水系统概化为均质的各向同性的三维非稳定流。

2.3 模型边界条件的概化

地下水流向主要从泥河到松花江，泥河、大荒沟、少陵河以及漂河的流向切割了该区域的含水层并且地下水与它们有补排关系。当河流或湖泊切割含水层，两者有直接水力联系时，可作为第一类边界条件处理[9]，因此将泥河、大荒沟、少陵河以及漂河作为定水头边界。由于其他边界和地下水流方向平行，所以把其定为地质零通量边界，没有侧向流量的补给，所以作为隔水边界。

2.4 数学模型

数学模型描述了一个研究领域即地下流动的水流，并确定数学方程的解决方案，这些条件构成一个真正的数学结构问题。含水系统概化为均质的各向同性的三维非稳定流，其相应的偏微分方程以及定解条件如下：

式中：h为水位，m；t为时间，d；W为源汇项，1/d；K为渗透系数，m/d；Ω为模型模拟区；Γ1为第一类边界；Γ2为第二类边界；n为边界面的外法线方向；H0为地下水初始水位函数；q(x,y,z,t)为第二类边界上已知流量函数。

3 模型的识别与验证

3.1 研究区的剖分

对于研究区的剖分采用有限差分方法求解，并且采用强隐式法联立迭代求解代数方程组[10]，在水平方向上对潜水含水层用相互垂直的平行线对研究区进行网格剖分，对研究区剖分成100×100的矩形单元格，其中有效单元计算点共3 164个。

3.2 计算时限与步长

识别期为2015年1月1日至2016年12月31日，验正期是2017年1月1日至2017年12月31日。由于地下水位受地下水开采与利用的影响，所以在对模型进行输入参数时将每个季度作为一个地下水开采期。

3.3 参数的选取

3.3.1 渗透系数K

根据巴彦县当地的水文地质条件以及含水层岩性，将其划分为3个不同渗透系数值的区域，见图2。通过查阅资料以及查看巴彦县的水文地质图可以发现，巴彦县的含水层岩性基本为亚黏土、粉土质砂及粉砂等地质；根据水文地质手册可以得到，所有分区的渗透系数经验值的区间为0.5～2 m/d。不同的分区输入相应的渗透系数初始值，由于此次模型属于各向同性，所以水平渗透系数KX、KY与垂向渗透系数KZ相同。

图2 研究区渗透系数分区图

3.3.2 降雨入渗系数

通过查找资料可知，该地区是含水层研究的重要组成部分，横向渗透、边境供应和其他形式的补给作用相对较弱，且观测井的埋深都在1～5 m之间，所以判断巴彦县的地下水补给主要为大气降水补给。根据巴彦县水文地质条件、流域水系、行政分区现状、地貌和岩性等特征，把研究区的大气降水渗透性分为6个区。见图3。

图3 研究区大气降水渗透性分区图

3.3.3 初始水头

此次模型中初始水头为各个观测井中2015年1月1日的水头值。观测井位置见图4。

图4 观测井位置图

3.4 模型的识别

模型采用间接求参法达到识别目的, 通过不断调整参数误差的方法对测井计算的初值，使其计算水头和测量水头符合拟合要求, 即计算水头和实测水头之间相差很小。若计算水头和实测水头相差较大，说明设定参数或边界条件不合适，需调整，只能通过模型识别、调整、再识别、再调整的多次反复过程，一直到计算水位与实测水位误差满足要求为止，此时所选用的参数即为校正后的参数。观测井校准水位图见图5，误差分析见表2，观测井水头时间图见图6，单个观测井水位拟合图见图7(以5号观测井为例)。

图5 观测井水位校准图

图6 观测井水头时间曲线

图7 5号观测井水位拟合曲线

表2 误差分析表

通过以上图表可以发现，各个观测井的计算水位与观测水位变化已基本一致，除个别点处的水位变化步调有偏差外，大部分点的变化基调基本上是同步的。通过该图可知，该模型具有良好的拟合效果，说明建立的数学模型、各类边界条件概化、水文地质参数的选择、源汇条件的处理是合理的，基本反映了研究区内地下水真实的变化特征[11]。

3.5 模型的验证

为使建立的数值模型更为准确地描述出地下水流系统的真实情况，需要利用现有的另一个时段的地下水动态观测资料对已经建立的模型进行检验，以此来进一步证明所建立的数值模型和所得模型参数的可靠性，本次验证的初始时间为2017年1月1日。观测井校准水位图见图8，误差分析见表3，观测井水头时间图见图9，单个观测井水位拟合图见图10(以5号观测井为例)。

图8 观测井水位校准图

表3 误差分析表

图9 观测井水头时间曲线

图10 5号观测井水位拟合曲线

从模型识别结果和拟合效果可以看出，地下水位拟合程度非常高，建立的模型可以真实地反映该地区实际的地下水情况，同时也论证其数学模型、边界条件、水文地质参数和源汇条件建立的正确性。

3.6 水文地质参数

经过模型的识别与验证可以得到，研究区各个分区的渗透系数和降雨入渗系数值见表4、表5。

表4 分区渗透系数表

表5 分区降雨入渗系数表

4 地下水位预测

本文建立的数值模型通过不同监测井的结果经过分析与计算得到含水层的含水量变化情况，并根据线性回归方法对研究区地下水位进行预测。大多数情况下，仅有流经方向和降雨量等一些自然环境影响地下水位。将模型初始时间设定为2015年1月1日，不改变已经建立的模型的任何条件与参数，对未来10年的地下水位进行预测。研究区未来10年地下水位预计变化情况见表6,地下水位折线图见图11，单个观测井地下水位折线图见图12-图14。

表6 地下水位预计变化情况

图11 地下水位预测折线图

图12 4号观测井水位预计变化图

图13 5号观测井水位预计变化图

图14 3号观测井水位预计变化图

根据表6、图11总体来看，巴彦县地下水位将会逐年下降，造成这一现象的根本原因是地下水开采量大于补给量。巴彦县经济主要以农业生产为主，是黑龙江省的粮食主产区，由于地表水的供水能力有限，需要不断开采地下水来维持农业的发展。地下水开采主要用途是农业灌溉，部分用于乡镇生活用水、工业生产和生态环境改善等方面。

根据图12-图14可以看出，巴彦县不同地区地下水位的变化情况并不相同，大部分地区地下水位有下降趋势。4号井位于洼兴镇附近，地下水不仅用于农业灌溉，附近的城镇居民生活用水以及工业用水也需要通过地下水进行补给，所以位于城镇附近的地下水位下降很明显。5号井位于村庄附近，地下水基本只用于农业灌溉，地下水位只是略有下降。3号井位于少陵水库和漂河附近，能够得到足够的地表水补给，这一地区的地下水位会有些上升。

针对地下水过度开采的情况，应该严格控制地下水的开采量，可以通过增加工厂企业的中水回用，取缔新水，逐渐降低对地下水的依赖。还可以考虑通过建设各种大中小型蓄水工程，在汛期时，对过境洪水进行拦蓄使用。同时，相关部门应加强监测站网的建设，进一步积累和完善监测资料，提高地下水资源评价成果的可靠性和精度。

5 结 论

本文以巴彦县为研究区域，根据前期水文地质资料的搜集、水文地质勘察与试验结果，建立了水文地质概念模型， 并运用Visual MODFLOW软件进行三维地下水流模拟，经过模型识别与校正，建立了符合目标的地下水流模型。得到以下结论：

1) 农业用水主要来源是地下水，同时随着城市规模的不断扩大，人口增长和经济建设的迅猛发展使得工业和生活用水量日益增大，区域地下水资源供需矛盾日益突出。

2) 通过建立巴彦县地下水流数值模型，经过模型的识别与验证，可以确定巴彦县地区的渗透系数为：河漫滩地区0.8 m/d；山区1.4 m/d；山前地区1.7 m/d。