海潮顶托作用下某铝厂污水处理站泄漏引发的地下水污染预测模拟*

张钰 蓝俊康

(桂林理工大学环境科学与工程学院 广西桂林 541004)

0 引言

海滨地带含水层中污染物的运移模拟一直是个很大的难题[1-3]，这主要是因为模拟时不仅要考虑海水水质的干扰，还要考虑海潮顶托作用对地下水流速的滞缓效应。广西企沙半岛拟建一个铝生产基地，在项目建设之前需要预测项目未来在“非正常状况”(即突发泄漏的情况)下对周边地下水环境的影响程度及污染物的迁移速度和方向，以便使人们能够在泄漏事件发生前采取科学的防范措施，有效阻止地下水污染所带来的危害。为了作出这些预测，本文采用了GMS软件进行数值模拟。

1 研究区的水文地质条件

1.1 自然地理状况

项目位于广西防城港市企沙镇，研究区的地貌以侵蚀堆积平原为主，海拔标高在-0.5～31 m之间，调查区南侧、东侧与西侧均为北部湾海域，仅东北侧与陆域接壤。研究区地处湿热多雨的南亚热带季风气候区，雨量充沛，多年平均降雨量为2 401.2 mm，多集中在5—9月份。

北部湾海潮汐的变化特征：研究区的海潮汐属全日潮，其显著特点是每月大潮过后约有2～4 d每日2次的高低潮，每年每日一次高低潮的天数约占全年的60%～70%。潮位以农历5月、10月为最高，即天文大潮期，农历2月、8月为最低，即小潮期，最大潮差5.56 m，平均潮差2.39 m，设计高潮位p=5%为3.43 m，p=10%为3.31 m。

1.2 地层岩性

根据区域地质资料及现场水文地质测绘查明，研究区主要有以下3个含水岩层。

(1)第四系近现代海相沉积层(Q4-2m)砂砾孔隙含水岩组。该层主要分布在海岸地段及常被海水淹没的地段，其岩性以砂质淤泥、淤泥质砂层为主，在海边过渡为细砂层。地下水主要赋存于松散的沙层中，与海水水力联系密切，弱富水性。

(2)孔隙裂隙含水岩组(J1)。其岩性为粉砂岩、泥岩夹细砂岩、泥岩夹粉砂岩、石英砂、砂砾岩等。地下水主要赋存于岩层的孔隙和构造裂隙中，弱富水性，一般泉水流量0.092～2.724 L/s，地下水枯季径流模数0.64～10.61 L/(s·km2)。

(3)裂隙含水岩组(S1ln)。含水介质为志留系连滩组的第三组(S1lnc)：上部为泥质砂岩夹粉砂岩，中部为泥质砂岩与粉砂岩互层，下部为泥质砂岩。地下水主要赋存于构造裂隙中，弱富水性，泉水流量一般为0.034 L/s，地下水枯季径流模数0.49～19.47 L/(s·km2)。

1.3 地质构造及地震

根据区域地质资料及现场水文地质测绘查明，有一条断层通过了厂区(即图1中的F1断层)。该断层的特征如下：该断裂形成于华力西期，在燕山期活动。断层走向约为52°，产状为142∠35°，从厂区的西北侧穿过并切断了光坡背斜轴部和企沙盆地。F1也被一系列北西向断层所错断。区域内出露长度约13 km，断层带宽约9 m，在地貌反映为部分形成北东向洼地。

1—孔隙潜水，水量中等；2—孔隙潜水，水量贫乏；3—孔隙咸水，大面积分布；4—碎屑岩类裂隙孔隙水；5—构造裂隙水，地下径流模数<6 L(s·km2)；6—实测及推测正断层；7—钻孔；8—厂区边界；9—模拟区分界线和地表分水岭；10—排水沟；11—地下水流向(注：红线为厂区范围，蓝线为边界)。图1 区域水文地质情况

研究区的抗震设防烈度为6度，设计基本地震加速度值为0.05g，地震分组为第一组，自1971年以来，仅发生过的最大震级为2.9级。

1.4 地下水的补径排特征

研究区内的地下水主要是接受大气降雨补给，其径流、排泄则主要受地形的控制，受局部地势的影响，地下水流向多变，但总体的径流方向是由陆域向四周的海域扩散。厂区范围内的地下水流向大部分向北流入虾萝江，另有一小部分向东南方向渗流，通过河湾流入北部湾。

海潮顶托的证据：利用2018年研究区内地下水水位观测点(98个)实测水位数据来绘制丰、平、枯期的等水位线图，在作图时发现，只有在地下水位高于2 m以上的区域才能绘制出地下水等水位线图，而低于2 m的水位值则杂乱无章，毫无规律，无法绘制出1.5、1.0和0.5 m的等水位线。这显示出地下水位低于2 m的地段均受到了海潮顶托作用的影响，使其地下水流方向不定。

1.5 水文地球化学特征

根据研究区的7个监测孔的采样分析，结果显示：该区域地下水pH值为5.58～7.57，总硬度31～252 mg/L，总溶解性固体(TDS)101～788 mg/L，氯化物8.31～208 mg/L，硫酸盐2.76～84.4 mg/L,铁(Fe) 0～3.09 mg/L。数据显示出有部分地段地下水中的Cl、SO4、TDS偏高，显然受到了海水入侵的影响。这也表明，研究区地下水的水质受到海水水质的影响，模拟时需要考虑这个因素。

2 地下水流场的数值模拟

2.1 地下水流场的数学模型

考虑地下水垂向流动，把地下水流视为三维的非稳定流，据《环境影响评价技术导则 地下水环境》(HJ 610—2016)附录D.2.2.1得：

(1)

式中，h为地下水水头，m；Kx、Ky、Kz分别为x、y、z方向的渗透系数，m/d；B1为已知水头边界(第一类边界)；B2为隔水边界；h1为河流水位，m；W为源汇项强度，d-1；Ω为渗流区域；μs为储水率，m-1，根据经验值取0.000 5。

2.2 边界条件

平面分界：模拟区的西侧以排水沟(图1中CD段)为地下水数值模拟时第一类边界，南部(DJ段)海岸线为第一类边界，北侧以虾萝江(CE段及EF、GH段)为第一类边界，FG段为地表局部分水岭，作为第二类边界(隔水边界)。东侧以清丫江的河岸线作为第一类边界。

移动边界：受海潮波动的影响，在海滩及入海河流内，海水水位和海陆分界线都是不时地发生变化，因此南部(DJ段)、北侧的虾萝江段、东侧的清丫江河岸线等处的边界位置是不固定的，作为移动边界对待。移动边界在模拟时的处理迄今有很多学者探讨，各自方法不同[4-6]。本研究的做法是把模拟时段细分，每个时段仅为6 h，不同每个时段赋予不同的海陆边界(按地面标高等值线来确定这些边界的位置)。在预测时，由于未来的海潮位是不可预测的，本研究则选一个典型日潮(该典型的日潮差等同于多年平均潮差2.39 m)，按该典型日潮的每时刻潮位变化值来确定这些移动边界的位置(潮位值如图2中A、B、C、D点所示)。

图2 北部湾某日海潮的变化情况

垂直分界：上界面以潜水面为界，为自由边界；其高度根据实测水位值及插值求得，底部以厚层状微风化基岩(泥质砂岩、泥岩)的顶面为隔水边界，即第二类边界。

2.3 模型的概化

按含水介质的风化程度把含水层分为3层(即表层松散的粉质黏土层、强风化粉砂岩层、中强风化粉砂岩层)，而把微风化岩层的顶板作为含水层的隔水底板，微风化层不纳入到模拟的单元网格中。

2.4 水文地质参数

(1) 初值的确定：①包气带的入渗系数值：本模拟区内的表土为志留系连滩群岩石风化而成的粉质黏土，参考前人的研究和测试结果，降雨入渗系数取值如下：枯水期0.10，平水期0.09，丰水期0.07；②含水层的渗透系数：按钻孔注水试验的结果赋予初值。

(2) 反求参数：将调查区5次所实测的地下水位、海水位值输入模型中，通过调整有关参数，按非稳定流拟合出各观测点的地下水位。当达到最佳拟合效果时，所得的水文地质参数及参数分区即为所求。拟合水位与实际观测水位的拟合情况见图3(由于水位观测孔多达98个，故仅显示其中的20个监测点的拟合效果)，从该拟合图可见，拟合效果较为理想。

2.5 模型的离散化

(1) 空间的离散化。模拟平面上共为100行，100列，共计10 000个矩形单元体。地形高程以散列点的形式输入到模型中，然后运用IDW插值法进行赋值，结果得到的地形图如图4所示[7]。

图4 模拟区立体剖分图

(2)时间的离散化。在进行溶质运移模拟时，需要利用水流的非稳定流模拟结果，此时，需要对时间进行离散化处理。由于海潮波动大且无规律可循，为了对应不同海潮下的边界，把模拟时段设定为6 h。

3 溶质运移模型和地下水的预测

3.1 溶质运移模拟

三维水动力弥散问题的数学模型为[8 -10]：

(2)

初始条件：

C(x,y,z)=C0(x,y)，(x,y)∈Ω,t=0

第一类边界条件(给定浓度边界)：

C(x,y,t)|Γ1=C(x,y,t)，(x,y)∈Γ1,t≥0

(3)

第三类边界条件(给定溶质通量边界)：

t≥0

(4)

根据於红等[11]在定边东北部潜水含水层野外弥散试验测定结果，细砂、粉细砂、粉土及粉质粘土的纵向弥散度为0.005 8～0.006 4 m。因此，本项目取的纵向弥散度为0.006 m，横向弥散度则按纵向弥散度的1/8取经验值[12]。

3.2 氧化铝工程污水处理设施泄漏事故结果预测

本项目废水处理站调节池体积为30 m×36 m×10 m，其浸湿面积为2 400 m2。假设渗水速率为4.8 m3/d，且渗漏事故发生30 d后才发现并处理完毕，按表土渗透系数0.2计算得废水向地下的泄漏总量为28.8 m3。据本项目可研报告所提供的污水组分及浓度值，选取预测因子为氟化物(F-)和COD，其初始质量浓度(源强)分别为：氟化物(F-)为9.46 mg/L，COD为250 mg/L。

3.2.1 COD的预测

泄漏点距监测孔GK7较近，根据GK7孔的监测结果，地下水中COD的本底值为2.6～2.8 mg/L，因此把2.8 mg/L作为污染晕的本底值，并把COD质量浓度3.0 mg/L(地下水Ⅲ类水质标准)作为其污染晕的分界线。COD可在微生物作用下发生降解，根据目前的研究发现，COD降解的速度解符合一级反应动力学规律, 即COD降解速率方程为[13]：

Ct=C0·e-kt

(5)

式中，Ct为时间t的COD质量浓度，mg/L;C0为COD渗入含水层的初始质量浓度，mg/L;t为COD在含水层中停留的时间，d;k为衰减速率系数，d-1。

根据以往实验研究发现，衰减速率系数k与水温、COD初始浓度、含硫、氨等无机还原性物质的浓度、溶解氧的含量(DO)等多种因素有关。一般情况下，地表环境中衰减速率系数k=0.023～0.086 d-1之间[14-15]。本研究区COD分解的地下水水温为18～20°C之间，所需的溶解氧主要来自大气降雨所携带，鉴于模拟的含水介质是渗透性较差的细颗粒土层，雨水深入速度缓慢，使得地下水环境中的溶解氧较少，因此取k=0.004 d-1(相当于比在大气中的分解速度慢5.75～21.5倍)来进行模拟。

利用GMS对含水层中COD污染晕的预测结果如表1所示。

3.2.2 氟化物(F-)的预测

预测结果如图5和表1所示。

(a)第100 d (b)第1 000 d (c)第10 000 d图5 泄漏发生100 d、1 000 d、10 000 d后含水层中F-污染晕的迁移

表1 COD、F-在含水层中的扩散特征

表1显示，在泄漏事故发生后，地下水环境的污染表现出与其他地方不同的特征，即污染中心点的F-的浓度会随时间的推移呈增大的趋势。对这一不寻常现象经分析后认为，这主要是海水的的顶托影响所致，因为在海潮顶托作用下，泄漏点的地下水流基本处于停滞状态，污染物难以向下游迁移。而随着时间的推移，从包气带渗入的污染物越来越多，就使得泄漏点下方含水层中的污染物越积累越多，污染晕中心点的浓度就会随时间的推移而越来越高。

4 结论

(1)泄漏的矿液或污水在渗入地下水后，将沿地下水流向北扩散到达虾萝江，在海潮出现时由海水带进北部湾。由于被污染的虾萝江段就位于厂区范围内，因此从泄漏点至虾萝江之间并无村民居住，也无抽水的民井，因此即使偶然发生泄漏事故将造成地下水COD、F-等含量的大幅度升高，但对周围居民饮水安全的影响忽略不计。

(2)虾萝江段距很短，污染物抵达虾萝江后，将立即被带入北部湾，大量的海水完全能稀释这微小的泄露量。另外，铝元素本身的毒性小，《海水水质标准》(GB 3097—1997)中对Al也没有作限量规定。

(3)从污染物的扩散特征预测表还可以看出，泄漏的矿液或污水在渗入地下水后，其迁移速度十分缓慢，这主要是因为泄漏点及下游大部分区域的含水层受海潮的顶托作用，地下水流速十分缓慢甚至停滞不前，污染物的扩散基本是以分子扩散的方式进行，含水层的净化作用基本靠下渗雨水的稀释作用来完成。