针对城市污水管网渗漏预测的数字模拟研究

宋雷震

(淮南联合大学 智能制造学院，安徽 淮南 232038)

城市化进程给人们生活带来便利的同时，也引发城市污水管网外渗污染问题，同时建筑工程开发过程中的基坑降水也导致了地面沉降和地下水降落漏斗[1]等问题。地下水模拟系统(Groundwater Modeling System，GMS)能够对地下水位进行数值模拟[2]。GMS软件是在现有地下水模型的基础上建立的全面综合性的图形界面软件，是目前已经被广泛应用于世界范围内的地下水模拟方法[3]，它具备强大的数据分析处理功能、方便快捷的操作界面、三维可视化的效果。本文利用GMS软件对城市污水管网外渗污染地下水进行地下水流和水质的数值模拟分析，并运用溶质运移模型预测排管泄漏的情况。

1 管网外渗污染地下水的数值模型

1.1 水文地质条件和地下水环境

本研究选取某市区域对城市污水管网外渗污染地下水情况进行数值模拟和预测。该地区气候为温带湿润季风气候，海拔高度约为60 m，相对高度差为8，全年平均气温为6.3 ℃，最低气温和最高气温分别为-31 ℃和35 ℃，无霜期和封冻期均约为150 d，全年平均日照时间约为2 700 h[4]。采集气象站该区域2005—2018年的年降水量进行处理分析，分别得到偏态系数、变异系数、降水量均值3个统计参数，数值分别为0.30，0.75，656.68 mm。2008—2018年系列降水频率曲线如图1所示。从图1可以看出，降水频率范围为3%～96%且集中分布在15%～85%。当降水频率为75%时，相应年平均降水量为524.7 mm。年降水量和频率可以拟合为一条光滑的曲线。

图1 2005—2018年系列降水频率曲线

研究区域的地下水广泛存在于圆砾土层、粉细砂、中粗砂、粉细砂中。自然状态下水位年幅度变化为3 m，同时基坑降水产生降落漏斗，该漏斗以城区为中心，面积约为16 km2。整个研究时间段内，地下水埋深为1.49～10.25 m，温度约为8～12 ℃，水质主要含有重碳酸钙钠和重碳酸钙，矿化度为0.5 g/L，铁离子含量丰富，质量浓度通常为2～11 mg/L，部分地段浓度更高。研究区域处于建设阶段，排泄方式为基坑降水、外围机井灌溉、植物蒸腾作用。地下水通过径流直接流入研究区域，补给源为侧向径流、水田区灌溉回渗水、降雨[5]。城区内地下水类型为人工开采型，城区外为气候型。利用渗水试验得到该区域包气带的渗透系数范围为0.08～1.44 m/d，渗透深度范围为40～140 m。含水层水文地质参数值如表1所示，展示9种钻孔编号在不同降深下的不水文参数值。通过附加人工流场下稳定流抽水的野外弥散试验可知，纵向和横向弥散系数分别为0.5 m2/d和0.015 m2/d，纵向弥散度和横向弥散度分别为0.212 8 m和0.006 4 m。

表1 水文地质参数情况

依据GB/T 14848—2017《地下水质量标准》对研究区域的地下水水质进行评价，通过指标值的限值范围得到地下水质量情况，检测得到地表水质量情况如表2所示。因此，研究区域菌群总数、总大肠菌群、铝、锰、铁等指标均超出集中式饮用水源限定标准，该区域的地下水质情况被人类生产生活影响。通过地下水水位等值线图，分析得到该区域形成以2个施工工地为中心的沉降漏斗，中心水位埋深为9.14～10.20 m。但区域历史资料显示漏斗中心下降的深度约为6 m，地下水流场从天然状态变成从4个方向流向新城。

表2 地表水质情况

1.2 地下水流和水质数学模型

GMS模拟系统包括UTCHEM、RT3D、SEEP2D、MT3DMS、MODFLOW等计算模块和Solid、MAP、PEST等辅助模块。研究依据区域含水层情况，将其设置为松散岩类孔隙微承压和潜水层含水层。设置地下水水力坡度在2/1 000以内，水流运动方式为平面二维流，水力特征遵从达西定律，同时当前和预期内水位均为非稳定流运动。研究依据建设区阶段和建设完成后等水位情况，均将研究区域边界设置为二类边界[6-7]。整体区域面积23.6 km2，人为干扰下和天然状态下边界条件划分情况如图2所示。根据图2中线条和箭头所表示的河流边界及水流方向，将区域划分为I1和I22个部分。

(a)人为干扰下边界条件

结合前述分析，研究建立非匀质各向同性松散岩类孔隙微承压水和潜水的地下水流模拟数值模型。孔隙潜水和孔隙微承压水的数学模型计算公式分别如式(1)和式(2)所示。

(1)

(2)

式中：K、T分别表示潜水含水层的渗透系数和承压含水层导水系数；μ*和μ表示承压含水层弹性释水系数和潜水含水层给水度；H、H0、h分别表示承压水水头、潜水水位、地下水位；B表示含水层底板标高，Qd和Qr表示排泄强度和入渗补给强度；Qi表示大井开采量；h0和h1表示初始水位和一类边界点的水位；q表示二类边界单宽流量；x，y为坐标；D表示计算区范围；Γ1和Γ2分别表示一类边界和二类边界；t表示时间；n表示开采井总数。

通过一类、二类边界条件、初始条件、偏微分方程即可求解问题。具体求解过程为：首先在区域范围进行线性插值和三角形剖分，利用伽辽金有限单元法离散数学模型为单元方程组，即可通过计算机程序求解问题。全部等值线图利用Surfer软件绘制，经分析检测没有问题的情况下，通过GMS软件实现数据自动采集和单元自动剖分。区域剖分使用Modflow模块，含水层分为2层，每层含有30 000单元，单元平均面积为1 405 m2，活动单元格16 807个，该种求解方法在保证计算准确率的前提下，也能够提升运行速度[8-9]。鉴于研究区域水位信息统计资料只有2019年1月的数据，因此不能实现模型识别和检测。研究利用区域平水年的源汇项进行模型识别，假如地下水位通过平水年的源汇项运行呈现稳定的变化状态，则模型参数设置合理，具有很好的可行性。在此基础上，初始流场设置为地下水流场，接着将2008年作为模型的源汇项，连续运行10 a。1个应力期为1个月。降水入渗补给量的计算公式为：

Qpr=0.1×α×P×F。

(3)

式中：α和P分别表示降水量和降水入渗补给系数；F表示一分区面积减去水面的面积。侧向径流补给量为：

Qlr=36.5×K×I×M×L。

(4)

式中：K和I表示含水层渗透系数和垂直于剖面方向的水力坡度；M和L分别表示含水层厚度和剖面长度。研究分析水位低于4 m的蒸发排泄情况，潜水蒸发排泄量Qe为：

Qe=c·F·E601。

(5)

式中：c表示潜水蒸发系数；F表示计算面积；E601表示水面蒸发量。在此基础上，利用MT3DMS模块构建地下水溶质运移模型，完成对污染的预测，选用的模拟因子为化学耗氧(O2计)[10]。二维非稳定流溶质运移方式下弥散数学模型为：

(6)

求解问题采用3rd-orderTVD方法和有限差分析方法。I1区域预测5 a后完成建设，因此5 a内地下水流场仍然会受到人工因素影响，研究选取2种特征污染物完成管网外渗过程，设置模型运行10 a。

2 外渗污染地下水数值模拟分析

2.1 水质模型识别结果

研究首先沿验证所提出地下水流数值模型的可行性，2层区域的初始设置值如下：第1层分区I1区域和I2区域面积依次为12.25 m2和11.35 m2，渗透系数依次为14.43 m/d和24.48 m/d，给水度分别为0.10和0.12。第2层分区I1区域和I2区域面积与第1层分区一致，渗透系数依次为80 m/d和100 m/d，给水度分别为0.18和0.20，弹性释水系数均为0.000 1。经侧向径流补给量计算可知，渗透系数为14.43 m/d，水力坡度为2%，含水层厚度为30 m，侧向补给值为688.12万m3/a。农业灌溉回渗量计算可知，农业灌溉量和回渗量为134.46万m3/a和20.17万m3/a。人工开采包括基坑降水和农业灌溉，总值为836.2万m3/a。

I1区域的降水入渗系数弱于I2区域，这是由于I1区域比I2区域早建设完成，大面积硬化路面覆盖导致降水入渗情况较差，而I2区域由于还存在大规模的耕地，因此降水入渗情况更好。2个区域的降水量补给计算结果如图3(a)所示。I1区域和I2区域降水入渗补给量分别为38.7万m3/a和107.58万m3/a。潜水蒸发排泄量情况如图3(b)所示，I1和I2的面积分别为1.8 km2和0.11 km2，地下水蒸发量分别为21.78 万m3/a和1.33 万m3/a，水面蒸发量均为1 210 mm，蒸发系数均为0.1。最终地下水补给资源总量和地下水总排泄量达到动态平衡，补给误差仅为4.74 万m3/a。

(a)降水入渗补给资源量

漏斗中心和QT05标记处的地下水位变化情况分别如图4(a)和4(b)所示。从2008—2018年整体来看，2个地段的地下水位表现出动态平衡的现象，同时实际动态情况也和变化幅度接近，因此模型的参数设计较为合理，所建立的数学模型和水文地质模型可靠性较高。

(a)漏斗中心

2.2 溶质运移模型结果

通过GMS软件中的MT3DMS模块完成三维地下水流数值模拟。结合尺度效应影响，潜水层纵向和横向弥散系数初始值设置为3.98 m2/h和0.313 m2/h，横纵弥散度比值为0.079。微承压层的层纵向和横向弥散系数初始值设置为27.34 m2/h和2.875 m2/h，横纵弥散度比值为0.105。区域内潜水含水层渗透系数为9.8 m/d，有效孔隙度为0.3。在人工加速流场下，纵向和横向弥散系数初始值设置为0.004 3 m2/h和0.001 5 m2/h，纵向和横向弥散度分别为0.004 6 m和0.001 6 m。区域内承压，含水层渗透系数为47 m/d，有效孔隙度为0.192。在人工加速流场下，纵向和横向弥散系数初始值设置为0.5 m2/h和0.015 m2/h，纵向和横向弥散度分别为0.212 8 m和0.006 4 m。

前述可知11 a降水系列数据可以计算得知降水均值、变态系数、偏态系数。实验选取丰平枯指标评判降水情况，分别将降水频率为25%、50%、75%的年份设置为枯水年、平水年、丰水年，计算结果如表3所示。枯水年、平水年、丰水年对应的降水量分别为506.8，621.3，738.5 mm。因此研究选取2008年作为污染物质的预测年份，通过GMS软件中的MT3DMS模块实现对城市污水管网外渗情况下点源连续渗漏对地下水污染情况预测，检测出研究区域有32处污水管道渗漏，模型预测污染物在5 a和10 a内的迁移过程。人为干扰下，10 a内化学耗氧量最大污染范围为0.29 km2，最大运移距离是1.88 km；天然状态下，10 a内化学耗氧量最大污染范围为0.14 km2，最大运移距离是0.99 km，如图5所示。

表3 丰水年、平水年、枯水年结果

管道泄露5 a后 管道泄露10 a后

3 结论

此次研究以城市污水管网外渗污染地下水源为例，利用渗水试验、水质分析等步骤完成区域水文水质分析，建立水流和水质数学模型并设置参数进行水源污染模型识别和预测。源汇项处理得到补给资源总量中大气降水入渗、侧向流入、灌溉回渗量分别为146.28万，688.12万，20.17万m3/a，排泄项中基坑降水、蒸发量、农业灌溉分别为701.74万，23.11万，134.46万m3/a，两者基本达到平衡状态，水文地质模型具有较高的识别效果。枯水年、平水年、丰水年对应的年份分别为2015年、2008年、2012年，相应的降水量分别为506.8，621.3，738.5 mm。由于所获取的数据有限，未对水流数学模型进行参数调优，这在后续研究需要完善。