某制药厂废水排放对地下水污染的数值模拟预测研究

杨 鑫,宁立波

(1.中国地质大学(武汉)地球科学学院，湖北 武汉 430074;2.中国地质大学(武汉)环境学院，湖北 武汉 430074)

制药工业关系着国计民生，是经济发展的支柱行业之一。但是，制药行业属于精细化工，制药废水的组成十分复杂，具有有机污染物种类多、浓度高，COD和BOD5值高，NH3-N浓度和固体悬浮物SS浓度高、色度深、毒性大等特征[1-3]。据统计，我国制药行业废水的排放量较大，约占全国工业企业污水排放量的2%[2]。制药废水在排放过程中污水管道发生渗漏可能造成的地下水污染问题一直备受关注，也是目前地下水污染中较难解决的问题之一，其中的关键问题在于地下水中污染物的运移特征及其范围较难确定。数值模拟方法是解决该类问题的途径之一，前人已做了大量的研究工作[4-6]，主要运用的数值模拟软件包括Modflow、GMS、Tough2等[6-9]，但针对制药废水泄漏而产生的地下水污染范围进行数值模拟的研究较少。新郑市属于郑州市规划的航空港区，其功能庞大，对环境质量的要求较高，但尚未进行地下水污染方面的环境影响评价工作。基于此，本文以拟进驻该区的某制药厂为研究对象，运用GMS软件对该制药厂废水排放过程中污水管道发生渗漏可能造成的地下水污染范围进行了数值模拟研究，以为该地区地下水资源的合理开发与保护提供科学依据。

1 研究区概况

该制药厂选址于河南省新郑市规划新区(见图1)，该地区属暖温带大陆性季风气候，多年平均蒸发量为1 859.7 mm，多年平均降雨量为683.64 mm。研究区地处黄河和双洎河冲洪积平原区，地形平坦，地表径流迟缓，地下水埋深较浅，且包气带岩性大部分为粉土，局部为粉砂，结构松散，极有利于大气降水的入渗补给。该地区地层以第四系为主，主要含水岩组为松散岩类孔隙含水岩组，其补给来源主要为大气降水的入渗补给，地下水的水化学类型较简单，以HCO3-Ca、HCO3-Ca·Mg型水为主。

图1 研究区位置平面图Fig.1 Location of the study area

2 研究区数值模型的建立

2.1 研究区水文地质条件与模型概化

水文地质概念模型是将含水层实际的边界性质、内部结构、渗透性质、水力特征和补给排泄等条件进行概化，便于进行数学与物理模拟[4-5]。水文地质概念模型的核心要素是边界条件、内部结构和地下水流态，通过对调查区的岩性构造、水动力场、水化学场的分析，可以确定水文地质概念模型的要素[6]。

2.1.1 研究区水文地质条件

研究区地下水属松散岩类孔隙水，含水层由第四系上更新统和全新统的细砂、中砂、砂砾石层构成。研究区双洎河河谷区的含水层岩性为砂、砂砾石，厚度为15～20 m，地下水水位埋深浅。双洎河河谷区的外围地带和双洎河北部支流的河谷地带，含水层岩性主要为粉土，局部夹少量薄层粉细砂和卵砾石，厚度为20～30 m，地下水水位埋深为25 m左右；研究区内的其他地区地下水水位埋深大，含水层岩性为黄土，夹少量的砂层，富水性差。研究区的浅层含水岩组主要由黄河多次迁徒、改道、泛滥冲积而成，在含水层的颗粒粗细、埋藏分布、渗透性能等方面具有明显的分带性；含水层大体由西南向东北方向延伸，呈单层出现，局部地段夹有粉土或粉质黏土薄层透镜体，厚度为20～40 m，个别地段厚度小于20 m；含水层顶板埋藏深度为5～50 m不等,岩性以粉土、粉砂为主，从上游至下游颗粒稍有变细，整体厚度约为60 m，属潜水含水层，下伏粉质黏土、黏土隔水层。

研究区浅层地下水主要补给来源是大气降水的入渗补给，此外双洎河的侧向渗漏也直接补给该区浅层地下水；浅层地下水的径流方向大致与地形变化趋势一致，总体方向由西南、西北山前向东、东南方向径流，浅层地下水的排泄方式主要有人工开采、越流排泄和径流排泄等，地下水埋深基本都大于4 m，蒸发排泄量很小，可以忽略。

2.1.2 研究区水文地质概念模型的建立

根据研究区的水文地质条件，构建了研究区的水文地质概念模型，见图2。

图2 研究区水文地质概念模型图Fig.2 Conceptual model map of hydrogeology in the study area

(1) 模拟范围：研究区模拟范围为40 km2，由于该地区浅层地下水含水层的埋深约60 m，结合钻孔资料，模拟垂直深度为60 m。

(2) 边界条件：①水平边界，即模拟区西与南左侧边界为定流量补给边界，东、北与南右侧边界为定流量排泄边界；②垂直边界，模拟区垂向地下水补给包括大气降水入渗补给、河流侧向入渗补给、灌溉入渗补给；地下水排泄主要为侧向径流排泄与人工开采。

(3) 含水层结构特征：含水层大体由西南向东北方向或由西向东延伸，本次模型含水层为层1(粉细砂)，隔水层为层2(粉质黏土、黏土)。

研究区水文地质概念模型的水力特征参数见图3，污染源为某药业污水处理厂管道破裂渗漏形成，为一线形连续性污染源。

图3 研究区水文地质概念模型的水力特征参数图 (负值表示排泄)Fig.3 Hydraulic characteristic parameters of the study area hydrogeological conceptual model (negative values represent excretion)

2.2 数学模型的构建

有机污染物在地下水中的运移非常复杂，影响因素除对流、弥散作用以外，还存在物理、化学、微生物等作用，这些作用常常会使地下水中污染物总量减少，运移扩散速度减慢。本次模拟中假设污染物质在地下水运移过程中不与含水层介质发生反应，可以被认为是保守型污染物质。所以，本文选择地下水流模型和溶质运移模型两个模型，联合求解水流方程和溶质运移方程，即可得到污染物质的空间分布。

2.2.1 地下水流模型

三维、非均质、各向异性的层流、非稳定潜层地下水流数学模型可表示为

式中：Kx、Ky、Kz分别为含水层在x、y、z方向的渗透系数(m/d)；Kxx为Kx在x方向上的投影；Kyy为Ky在y方向上的投影；Kzz为Kz在z方向上的投影;x，y，z方向为各向异性介质渗透系数的主方向;h为地下水水位(m);h0为地下水初始水位(m)；Kn为边界上垂直于等水位线方向的渗透系数(m/d);Ω为研究区范围；Γ2为边界；w为降雨入渗补给量(m3/a)。

2.2.2 溶质运移模型

不考虑污染物在含水层中的吸附、交换、挥发、生物化学反应，地下水中溶质运移的数学模型可表示为

式中：α为含水层的弥散度(m)；vi，vj分别为i和j方向上的速度分量(m/d)；|v|为速度模；C为模拟污染物的浓度(mg/L)；ne为有效孔隙度；i,j为1，2，3，x1=x，x2=y，x3=z；t为模拟时间(d)；C′为源汇的污染物浓度(mg/L)；W为源汇单位面积上的通量；Dij为弥散系数(m2/d)。

2.3 模型识别与参数确定

2.3.1 模拟流场及初始条件

本次模拟以2016年5月的地下水流场作为初始流场(见图4)，以2017年5月统测的地下水流场作为模拟流场。

图4 模型初始地下水流场Fig.4 Initial groundwater flow field of the model

2.3.2 模拟区网格剖分

模拟区网格剖分单元格为50 m×50 m，其网格剖分见图5。

图5 模拟区网格剖分图Fig.5 Mesh subdivision map of simulation area

2.3.3 模型识别与参数确定

模型的识别与验证是整个模拟中极为重要的工作，通常需要在反复地调整参数和某些源汇项的基础上才能达到较为理想的拟合结果[10]。本次模型识别与验证过程采用试估-校正法，该法属于反求参数的间接方法之一[11]。本次模拟时间为2018年1月到2022年12月，每个时间段内包括若干时间步长，时间步长由模型自动控制，严格控制每次迭代的误差。模型最终识别的研究区含水层水文地质参数分区见图6，模型识别的降雨入渗系数、含水层渗透系数和给水度见表1，灌溉入渗系数为0.10，该区地下水水位埋深普遍大于5 m，潜水蒸发排泄小可忽略不计。

图6 研究区含水层水文地质参数分区图Fig.6 Zoning maps of hydrogeological parameters of aquifers in the study area

编号降雨入渗系数水平方向渗透系数/(m·d-1)给水度层1层2层1层210.150.80.008 60.0470.03420.202.10.001 00.0470.03430.151.00.001 00.0470.03440.153.50.001 00.0340.03450.101.00.001 00.0340.03460.121.00.001 00.0400.03470.101.50.008 60.0470.034

研究区不同岩性介质弥散度的初值依据前人研究成果及经验值确定，最终确定研究区不同岩性介质的平均纵、横向弥散度，见表2。

表2 研究区不同岩性介质弥散度经验参数表

3 模拟结果与分析

污染源位置为遂成药业污水处理厂，假定污水管道出现长10 m、宽2 cm的裂缝，管道天然基础层的渗透系数取值为0.047 m/d，则渗漏量约为10×0.02×0.047×1 000=9.4 kg/d，管道在无防渗设置的情况下发生渗漏，即为连续的线性污染源。根据污水中主要污染物监测指标，考虑污染物与岩土之间的反应和自身降解性，选取CODCr为地下水污染预测因子，根据遂成药业污水处理厂进水水质状况，进水中CODCr浓度为419 mg/L。

本文运用GMS软件对拟建的遂成药业污水处理厂废水排放过程中污水管道发生渗漏可能造成的地下水污染范围进行了数值模拟预测研究，模型运行时间起始于2018年1月1日，模拟预测了100 d、730 d 和1 825 d(2022年12月31日)三个时间段内研究区地下水中污染物CODCr的最大运移距离和污染晕的运移分布情况，其预测结果见表3和图7。

污染组分的扩散具有弥散特征，初期是以污染源为中心向周围弥散，之后随地下水的流向污染晕的主轴方向与其趋向一致(向南运移)，所以，污染组分CODCr的检出范围以整个扩散范围进行统计，最大运移距离以污染源(见图7中颜色最深的中心点)为起点向东南方向的运移距离为准。预测结果表明：该制药厂废水排放过程中管道渗漏发生100 d后，地下水中污染物CODCr的检出范围为148 520 m2，超标范围为0 m2，地下水中污染物CODCr的最大运移距离为238 m；渗漏发生730 d后，地下水中污染物CODCr的检出范围为165 500 m2,超标范围为0 m2，地下水中污染物CODCr的最大运移距离为250 m；渗漏发生1 825 d后，地下水中污染物CODCr的检出范围224 000 m2，超标范围为0 m2，地下水中污染物CODCr的最大运移距离为362 m(见表3)。受该地区西侧与西南侧两个地下水降落漏斗的影响，该污染物CODcr渗漏后倾向于西南侧运移且运移速度缓慢，至规划年2022年末距东侧地表水体距离较远，对周边水体的影响较小。

表3 研究区地下水中污染物CODCr的预测结果

图7 研究区地下水中污染物CODCr的污染晕运移分布 预测图Fig.7 Prediction map of CODCr pollution halo transport in groundwater in the study area

由图7可见，研究区地下水中污染物CODCr的污染晕随地下水的流向而扩散，其扩散方向与地下水的流向一致；地下水中污染物CODCr的运移距离随时间的增加而增加，但整体上看在预测期内地下水中污染物CODCr的运移距离仍处在该制药厂的建设区范围内，不会对区域地下水造成明显的污染。

4 结论与建议

(1) 制药厂建设所排放的废水含有复杂且大量的有机污染物，可能会对厂区范围内的地下水造成污染，所以必须在制药厂建设前对废水排放可能造成的地下水污染范围进行预测研究。

(2) 本文运用GMS软件对拟建的遂成药业污水处理厂废水排放过程中污水管道发生渗漏可能造成的地下水污染范围进行了数值模拟研究，模型运行时间起始于2018年1月1日，模拟预测了100 d、730 d和1 825 d(2022年12月31日)三个时间段内研究区地下水中污染物CODCr的最大运移距离。模拟结果表明：渗漏发生100 d后，地下水中污染物CODCr的最大运移距离为238 m；渗漏发生730 d后，地下水中污染物CODCr的最大运移距离为250 m；渗漏发生1 825 d后，地下水中污染物CODCr的最大运移距离为362 m。

(3) 研究区地下水中污染物CODCr的污染晕随地下水的流向而扩散，其扩散方向与地下水的流向一致；地下水中污染物CODCr的运移距离随时间的增加而增加，但整体上看，在预测期内不会对区域地下水造成明显的污染。但是，鉴于该制药厂污染物的极大危害性，建议在厂区范围内及其周边设置监测井，定期检测地下水中污染组分的变化，一旦发现地下水受到污染，应立即采取相应的措施进行治理。