基于Visual MODFLOW 的红层地区某钛白粉厂地下水污染物运移模拟研究

吕堂安 LV Tang-an；伍静 WU Jing；黄伟 HUANG Wei；黄华娟 HUANG Hua-juan；朱敏杰 ZHU Min-jie

（①广西大学资源环境与材料学院，南宁 530004；②中国有色桂林矿产地质研究院有限公司，桂林 541004；③广西中科地质勘查有限责任公司，南宁 530022）

0 引言

数值模拟技术能直观了解污染物泄漏后运移影响范围及浓度变化，对地下水污染的防治有着极其重要的意义，其应用的关键在于对模型的概化、边界和参数的准确选取[1]。Visual MODFLOW 由于简便的操作性和强大的可视化能力，已成为国际上最流行的溶质运移模拟软件之一[2]。国内外学者已用其对地下水污染运移开展了广泛的研究。国外有：Morsy Samah Mahmoud[3-5]等人对埃及西米尼亚地区、蒂鲁普和卡拉季地区地下水中污染物进行预测，为开发区建设提供防治依据；国内有：刘莉[6]针对淄博某化工项目，预测不同工况下的地下水中污染物溶质运移变化规律；马国明[7]运用Visual MODFLOW 对宿州市某制药厂地下水中氨氮污染物进行了模拟预测，并提出了防治措施。

某改、扩建钛白粉厂地处红层区，周边人口众多，地下水是人们生产生活的重要水源，对该厂地下水的防治显得极为敏感和重要。因此，本文以该钛白粉厂为例，在详细了解地下水水文地质特征基础上，概化了水文地质模型，通过识别与检验，利用Visual MODFLOW 软件建立了二维数学模型，开展特征污染物数值模拟，了解了该区地下水中污染因子随时间的衰减及运移范围，对指导该地区地下水的防治和保护有着重要意义。本次研究也为同类型地区如何获取准确的水文地质参数，建立符合实际地下水流场的数值模型提供参考。

1 研究区域概况

1.1 地质条件 研究区位于广西博白旺茂镇，为红层丘陵盆地地貌，依据地形及岩层特性，将其划分为两个次级水文单元（图1A）。由水文地质剖面图（图1B）可知，研究区内地层为下白垩统大坡组（K1d），上部为黄红色粉砂土及部分出露含砾砂岩；下部为红色砂岩夹砾岩、紫红色砂岩含砾砂岩，中砂－细砂，中－粗砾，磨圆度较好，分选性差，成分各异，以砂岩、石英砂岩为主，少量灰岩砾石。研究区内无大的断裂构造，白垩系砂砾岩层产状整体上较平缓。

图1 研究区水文地质简图（A）；A-A＇ 水文地质剖面图（B）

1.2 水文地质条件 研究区内含水层主要为白垩系红层，表层风化带厚度一般为10～30m，红层风化裂隙较为发育，局部形成风化层裂隙水，向下为沿构造裂隙发育的中、弱风化带，其脆性的砂、砾岩在向斜轴部及构造变形强烈部位，形成较多的节理裂隙，存在构造裂隙水。鉴于红层表层风化裂隙水含水层厚度不大，与下伏构造裂隙水之间无隔水层，存在直接水力联系，故本次研究将含水层均概化为单层潜水含水层，地下水类型主要为碎屑岩类构造裂隙水，其厚度为50～150m。浅层构造裂隙水主要接受大气降水的补给；排泄主要为蒸发和侧向径流补给。底部脆性构造不发育的致密砂砾岩层胶结完好，无裂隙发育地段，视为隔水层。

该区地下水枯水期埋深一般在3m～8.65m，丰水期埋深一般在1.52m～5.74m，年水位变幅0.85～3.48m，水质类型以HCO3—Ca 型为主。根据厂区机井J 的抽水试验结果，该区含水层渗透系数为1.663×10-4cm/s，透水性中等，给水度μ 取0.085。鉴于本次研究区的红层构造裂隙含水介质，其持水度较小，因此，有效孔隙度n 近似等于给水度，n=μ=0.085，水力坡度值I 为20.8‰，实际流速u 为0.035m/d。

2 地下水数值模型建立

2.1 水文地质概念模型 根据地质与水文地质条件调查，将其概化为一个潜水含水层，浅层地下水运动以水平方向为主；地下水流向基本与下隔水层面平行，垂向分量忽略不计，符合二维流特征，因此概化为平面二维流[8]。本次模拟将区域地下水含水系统概化为非均质、各向同性二维无压非稳定流。研究区含水层厚度为50～150m，西侧以老屋水库为地下水排泄边界，视为定水头边界；东侧及西北侧以地下水分水岭为边界，视为二类边界；东北侧以金甲河为北侧排泄边界，视为定水头边界；上边界为降水补给、蒸发排泄边界，下边界以基岩作为相对隔水边界。

2.2 地下水流数学模型 在实际情况中，由于缺少三维实测数据，在必要时须做一些假设，这可能造成三维模型的结果失真[9]。在实际工作时，当含水层与弱透水层或隔水层的贮水系数比值大于10，而弱透水层或隔水层越流系数小于1·10-6m/d 时可用二维模型[10]。本次研究建立各向同性二维非稳定流数学模型，即式（1）[11]：

式中：μd为重力给水度，l/m；Kx、Ky分别为x、y 方向上的渗透系数，m/d；t 为时间，d；h 为潜水面高度，m；z 为含水层隔水底板高度，m；ε 为含水层汇源项，m3/d；Ω—渗流区域；H0为渗流区Ω 内各点的初始水头，m；Γ1为渗流区域Ω 的第一类边界；H1为第一类边界Γ1上的已知水头函数，m；Γ2为Ω 的二类边界；K 为边界Γ2法线方向的渗透系数，m/d；n 为边界Γ2上某点处的外法线方向；q 为Γ2上的单位面积流量函数，m3/d。

本次模拟网格剖分精度为20m×20m，在分水岭及污染物渗漏处加密网格为10m×10m，模拟面积约为2.17km2。研究区主要补给来源为大气降雨，老屋水库及金甲河为主要排泄边界。降雨蒸发参数根据《玉林市地下水利用与保护规划（2016-2030）》，多年平均逐月降雨量为1583mm，多年平均蒸发量为1506mm，蒸发临界深度取5m。

2.3 数值模型动态识别与检验 本次模拟以枯水期水位作为初始水位，按照差分法得到含水层初始流场，枯水期至丰水期为模型识别期。根据水文地质试验设置含水层渗透系数Kx和Ky初始值为1.663×10-4cm/s，采用试估—校正法对模型进行检验，获得与实际流场相符程度高的模拟流场。将枯、丰水期实测点位数据与模拟水位进行检验，置信区间相关系数分别为0.995 和0.972，表明模拟流场精度较高，适用于本次模拟研究。最终确定研究区水文地质单元Ⅰ和Ⅱ的溶质模拟渗透系数Kx、Ky取值为1.039×10-4cm/s。

2.4 溶质运移数学模型 本文从保守角度考虑，本次溶质模拟计算的模型识别和计算没有考虑污染物在含水层中的吸附、挥发、生物化学反应，考虑最大值[12]。地下水以水平运动为主，垂向弥散作用忽略不计。在此前提下，本次调查研究建立地下水二维溶质运移模型，二维水动力弥散方程的数学模型见式（2）[11]：

式中：C 为模拟污染物质的浓度，mg/L；θ 为地层介质的孔隙度，无量纲；t 为时间，d；xi为沿直角坐标系轴向的距离，m；Dij为水动力弥散系数，m2/d；vi为孔隙水平均实际流速，m/d；qs为单位体积含水层流量，代表源和汇，m3/d；Cs为源或汇水流中组分的浓度，mg/L。

弥散系数结合水文地质条件和参考《水文地质手册》[13]，最终选确定为8m2/d。

3 模拟预测

根据地质和水文地质调查，渗漏点位于钛白粉厂厂区西部，属于水文地质单元Ⅰ，地下水流向为北东向西径流，污染物渗漏后将向西侧老屋水库一带运移，因此本次模拟重点关注水文地质单元Ⅰ。

3.1 污染源设定 由工艺分析可知：钛白粉厂产生的酸性废水中主要污染物为COD、氨氮和氟化物，模拟渗漏区为钛白粉污水处理车间和硫铁矿制酸污水处理站。根据《环境影响评价技术导则地下水环境》（HJ 160），当系统老化等原因导致调节池渗漏≥2L/m3·d 时视为非正常状况，结合工程验收监测数据，计算得出钛白粉污水处理车间污水渗漏量为864L/d；铁矿制酸污水处理站污水渗漏量为96L/d，预测污染物持续渗漏60 天可发现。

3.2 污染物迁移模拟预测结果 ①COD 预测评价。渗漏液初始浓度为417mg/L，持续渗漏60 天时调节池底部地下水COD 浓度达到最大值4.80mg/L；第100d 时COD中心浓度最大为3.9mg/L，渗漏点下游0～12m 范围内COD浓度高于地下水Ⅲ类标准，溶质运移范围0～25m，未出厂界（见图2A）；渗漏第1000d 时，污染物中心浓度最大为1.00mg/L，污染物运移至下游70m（见图2B）；在渗漏第7300d 时污染物COD 已运移至厂界外，到达西侧老屋水库一带，中心浓度最大为0.05mg/L，浓度低于地下水Ⅲ类标准（见图2C），扩散距离为75～269m，位于下游的敏感目标地古岭民井距研究区420m，不受到研究区渗漏影响。渗漏后COD 平均运移速度为13.45m/a。

图2 COD 运移路径及浓度变化图；100d（A）；1000d（B）；7300d（C）

②氨氮预测评价。渗漏液初始浓度为33.75mg/L，持续渗漏60 天时调节池底部地下水氨氮浓度约0.38mg/L；渗漏第100d 时氨氮中心浓度最大为0.30mg/L，污染物运移至下游0～6m 范围内，未出厂界（图3A）；渗漏第365d 时调节池底部污染物氨氮浓度为0.21mg/L，污染物运移至下游0～41m 范围内，未出厂界（图3B）；渗漏第730d 时调节池底部污染物氨氮浓度为0.05mg/L，污染物运移至下游0～58m 范围内（图3C），远低于地下水Ⅲ类标准（0.5mg/L）。渗漏后氨氮平均运移速度为2.9m/a。

图3 氨氮运移途径及浓度变化图；100d（A）；365d（B）；730d（C）

③氟化物预测评价。渗漏液初始浓度为195.25mg/L，持续渗漏60 天时调节池底部地下水氟化物浓度约2.4mg/L；渗漏第100d 氟化物中心浓度最大为1.90mg/L，污染物运移至下游0～30m 范围内，未出厂界（图4A）；渗漏第365d 污染物中心浓度为1.10mg/L，污染物运移至下游0～41m 范围内，未出厂界（图4B），渗漏第1000d 污染物中心浓度为0.50mg/L，污染物运移至下游0～56m 范围内（图4C），此时浓度已经低于地下水Ⅲ类标准（1.00mg/L）。渗漏后氟化物平均运移速度为2.8m/a。

图4 氟化物运移途径及浓度变化图；100d（A）；365d（B）；1000d（C）

4 结论

本次研究表明在符合一定水文地质条件下的红层地区，采用Visual MODFLOW 能建立与实际流场相符程度较高的二维流模型。在此基础上进行模拟预测，表明该区地下水在钛白粉厂发生渗漏后对下游红层地下水有一定范围影响，但影响未扩散至地谷岭村屯。企业应重点对污水站下游和厂区、地下水下游边界监测水点进行监测，随时掌握不同时期地下水的水质水位动态资料。