尾矿坝加高对重金属Cu2+迁移影响的数值模拟研究

2022-04-22陈红丹

能源环境保护 2022年2期

孙杰，郝喆，陈红丹

(辽宁大学环境学院，辽宁沈阳 110036)

0 引言

尾矿库用于堆放矿山开采后产生的尾矿，是矿山开采中的重要设施[1]。它也是矿山开采中产生的人造污染源[2]，若管理不当，可能会造成尾矿砂外泄、污染物扩散等后果，严重威胁周围环境、资源及下游人民的生命财产安全[3]。尾矿坝是尾矿库中最主要的构筑物，由初期坝和后期尾矿材料形成的坝体共同构成[4]。尾矿废水中含有重金属污染物及其他有毒有害物质，其经过迁移作用会对周围的水和土壤造成污染[5-7]，因此，尾矿库引发的重金属污染问题愈加受到关注。

目前，许多学者在尾矿库污染物运移的数值模拟方面已有一定的研究进展。张建军[1]等应用GEO-SLOPE软件中的SEEP/W和CTRAN/W模块分析了尾矿砂各层渗透系数比、干滩长度等对渗流场的影响，并且模拟了考虑和不考虑吸附两种情况下100天、1 000天坝体内铜离子的迁移过程；陈红丹[8]等模拟了降雨、蒸发与地下水耦合作用下的大型尾矿库重金属Cu2+的迁移过程；周泽超[9]等利用TOUGH2软件根据建立的地下水流动模型与地下水污染物迁移模型，模拟计算出了各时刻铜离子在坝体中的分布情况；Yu Liren等[10]模拟研究了阿纳马市附近的索利莫斯河的流动与污染物的运移；B.O.Tsydenov[11]通过数值模拟分析了某海域秋季热坝发展过程中，污染物在不同风况下的空间迁移规律及分布特点；Ben Abdelghani Farouk[12]等运用Hydro-Geosphere软件模拟了污染物在露天矿山中随时间推移的迁移情况。上述研究多基于现有尾矿库研究重金属污染物的对流-弥散作用，未见在尾矿库堆积过程中，尾矿坝加高对于重金属迁移范围的影响研究。本文以本钢集团歪头山尾矿库为例，利用GeoStudio软件中的SEEP/W和CTRAN/W模块建立尾矿库数值模型，研究尾矿坝加高对于重金属Cu2+在尾矿库中迁移范围的影响，为大型尾矿库环境地质评价和污染预测防治提供依据。

1 研究区域及方法

1.1 研究区域

研究区域位于辽宁省本溪市溪湖区，该尾矿库属于二等尾矿库。库区占地面积1.62 km2，现库容量约为1亿m3。尾矿库初期坝由透水堆石坝构成，初期坝坝顶标高为191.9 m，坝底标高180.0 m，相对高度11.9 m，顶宽4.0 m，坝长410.0 m，内坡比1∶2，外坡比1∶2.5。尾矿库的主坝尾矿堆积，坝顶标高为280.0 m，坝库坝底标高180 m，主坝高度差约为100 m，坝顶坝长约为1 600 m。

1.2 研究方法

Geostudio系统软件是由加拿大岩土软件GEO-SLOPE公司开发的面向岩土、采矿、交通、水利、地质、环境工程等领域的仿真分析软件，是全球最知名的岩土工程分析软件之一。它包含以下八种专业分析模块：SLOPE/W(边坡稳定性分析)、SEEP/W(地下水渗流分析)、SIGMA/W(岩土应力变形分析)、QUAKE/W(地震动力响应分析)、TEMP/W(地热分析)、CTRAN/W(地下水污染物传输分析)、VADOSE/W(综合渗流蒸发区和土壤表层分析)。当用户需要进行多模块工作时，只需将模型的边界条件和材料参数等定义在一个模块中，使用其他模块时不需重复输入相关参数，因此该软件建模高效、操作简单[13]。本研究中所使用的Geostudio版本为2012版。

SEEP/W模块主要用于分析多孔渗水材料，它可以分析饱和稳态及饱和-不饱和瞬态问题，几乎可以对所有地下水问题进行建模分析，其计算原理是通过内插函数来表证该单元上每一节点处的水头值，单元内的水头分布用线性差值函数来近似。CTRAN/W模块主要用于分析污染物通过土层和岩石等渗流介质外界环境中的水体，经过地面和地下水的非饱和区域进入地下水体的问题，其计算原理是伽辽金有限元法，该方法从剩余加权法出发对对流-弥散方程进行数值离散，是求解微分方程近似解的一种方法。CTRAN/W和SEEP/W相结合，可以用来分析污染物的运移。SEEP/W模块计算水的流速、水体积含量和水流量，CTRAN/W模块则用这些参数计算污染物的迁移。SEEP/W和CTRAN/W模块均属于有限元分析软件。

运用GeoStudio中的SEEP/W和CTRAN/W两个模块开展尾矿坝加高对重金属污染物Cu2+在尾矿库中迁移规律影响的数值模拟研究。在SEEP/W模块中设置尾矿库模型材料参数、边界条件等参数，通过稳态渗流计算得出尾矿库孔隙水压力等势图、体积含水率等结果；利用CTRAN/W模块建立对流-弥散和粒子示踪模型，通过计算重金属Cu2+的对流-弥散过程，分析对流-弥散过程计算出的重金属浓度等势图、粒子路径迁移图等结果，得到Cu2+在尾矿库中迁移的情况，最终通过分析得出尾矿坝加高过程中Cu2+的迁移规律。

2 基础原理及模型建立

2.1 Darcy定律

Darcy定律是渗流理论中的基本定律[14]，是形容饱和土中水的渗流速度与水力坡降之间的联系的。其最初是从饱和土中得到的，同样也适用于非饱和土，不同之处在于渗流系数是否为常量[15]。表达式为：

Q=ki

(1)

式(1)中：Q——流速，m/d；k——渗流系数，m/d；i——渗透坡降，无量纲。

2.2 渗流基本方程

渗流计算方程[16]：

(2)

式(2)中：H——总水头，m；kx——x方向的渗流系数，m/d；ky——y方向的渗流系数，m/d；Q——边界流量，m3/d；θ——体积水容量，L/m3；t——时间，d。

2.3 对流-弥散方程

污染物的迁移机理受到对流-弥散、吸附-沉淀、解吸、络合、生物降解等作用。根据尾矿库各参数性质，本文仅考虑重金属的对流-弥散作用。对流-弥散方程为[17]：

(3)

式(3)中：Dx——横向弥散系数，m2/d；Dy——纵向弥散系数，m2/d；C——污染物浓度，g/m3；ux——横向渗流速度，m/d；uy——纵向渗流速度，m/d。

2.4 模型建立

本文以歪头山尾矿库为例，模拟在尾矿坝加高过程中重金属Cu2+的迁移过程，通过对比在不同尾矿坝坝高下重金属的迁移范围及其浓度分布情况，研究尾矿坝加高对重金属污染物迁移的影响规律。如图2所示，本次研究共选取三种不同坝高，工况1、工况2、工况3的顶部高程分别为236、263、280 m，三个概化计算模型剖面的地面高程均为180 m，计算区域的水平距离均为595 m。模型网格剖分为三角形或四边形，在模拟中已做网格无关化检验，网格尺寸大小的划分对重金属迁移的影响可忽略不计，本文尾矿库模型全局单元尺寸取5 m，工况1、工况2、工况3剖分单元分别为1 212个、1 450个、1 545个，节点分别为1 320个、1 551个、1 644个。

图2 尾矿库概化计算模型Fig.2 Generalized calculation model of the tailings pond

2.5 模型边界条件

在进行稳态渗流、对流-弥散模拟时，设置排水边界、零流量边界、上游水头、下游水头等水力边界条件和污染物边界条件。如图2所示，在稳态渗流模拟中，AB为零流量边界，BC为上游水头，DE为下游水头，EF为排水边界。在进行对流-弥散模拟时，BC为重金属Cu2+污染物边界，初始浓度为50 g/m3。

2.6 材料参数

在概化计算模型中，包含尾细砂、尾粉质粘土、基岩等三种不同土层，其中前两层为饱和/不饱和模型，基岩为饱和模型，各介质水文地质参数如表1所示，土水特征曲线及各土层渗透系数-基质吸力曲线详见图3和图4。

表1 各介质水文地质参数

图3 各介质土水特征曲线Fig.3 Soil and water characteristic curves of each medium

图4 各介质渗流系数-基质吸力特征曲线Fig.4 Seepage coefficient-matric suction characteristic curve of each medium

3 数值模拟计算结果分析

3.1 渗流分析

图5为不同坝高下的尾矿库孔隙水压力等势图，图中所示剖面1和剖面2交点高程坐标分别为(211.93 m，240.31 m)、(146.35 m，223.01 m)。图6为两个剖面在工况1高度、工况2高度和工况3高度下的垂直方向孔隙水压力等势分布情况。

根据图5和图6可知，在同一尾矿库模型中，垂直方向上孔隙水压力随着尾矿库深度的增加而增大，尾矿库底部的孔隙水压力最大，尾矿库顶部的孔隙水压力最小。且孔隙水压力随高程变化的曲线大致呈一条直线，所以孔隙水压力随高程变化的关系为一次函数关系，经过数据分析可得函数斜率为-9.77，即高程每增加1 m，孔隙水压力减少9.77 kPa。

图5 稳态渗流作用下尾矿库孔隙水压力等势分布图Fig.5 Equipotential distribution diagrams of pore water pressure of the tailing pond under steady seepage

图6 垂直剖面孔隙水压力分布Fig.6 Pressure distribution of pore water in vertical section

通过图6可知，在剖面1处，当高程为199.27 m 时，工况1、工况2、工况3的孔隙水压力分别为29.43、180.69、244.08 kPa；当高程为240.20 m 时，工况1、工况2、工况3的孔隙水压力分别为-378.70、-217.80、-148.49 kPa。在剖面2处，当高程为169.97 m时，工况1、工况2、工况3的孔隙水压力分别为251.81、362.86、408.23 kPa；当高程为223.00 m时，工况1、工况2、工况3的孔隙水压力分别为-272.86、-148.59、-96.36 kPa。由此可知，孔隙水压力会受到尾矿坝高度的影响，当尾矿坝高度增加时，相同位置的孔隙水压力增大。

3.2 重金属Cu2+迁移规律分析

根据实际水质检测结果，在不考虑吸附和离子间相互作用的情况下，研究区域的重金属Cu2+初始浓度为50 g/m3。通过对尾矿库进行对流-弥散和粒子示踪模型计算，分析研究尾矿坝加高对于重金属污染物的迁移规律。根据GeoStudio数值模拟2 920 d计算结果，得出重金属Cu2+在不同坝高下的浓度分布情况和粒子迁移路径，如图7和8所示。图7中粉色阴影表示Cu2+浓度高于40 g/m3，蓝色阴影表示Cu2+浓度低于10 g/m3。

通过观察图7可以看出，在同一尾矿库模型中，污染物浓度沿初始浓度边界向初期坝方向逐渐降低。通过对比发现，不同坝高的尾矿库重金属污染物Cu2+迁移的质量浓度等势分布不同，尾矿坝加高会对污染物迁移产生影响。在剖面1和剖面2处可以看出，当尾矿库深度相同时，尾矿坝越高，重金属Cu2+浓度越大。由图8粒子迁移路径可知，重金属Cu2+是随地下水沿尾矿坝方向迁移的，且随着尾矿坝的加高，污染物的迁移距离越大。

图7 重金属Cu2+迁移2 920 d浓度等势图Fig.7 Concentration equipotential diagrams of heavy metal Cu2+ after migration for 2 920 days

根据图9可知，在同一尾矿库模型中，Cu2+的质量浓度随着高程的增加而减小，尾矿库底部浓度最高，尾矿库表层浓度最低且接近于0。在剖面1处，当高程为169.97 m时，工况1、工况2、工况3的Cu2+质量浓度分别为6.59、28.65、35.14 g/m3；当高程为199.25 m时，工况1、工况2、工况3的Cu2+质量浓度分别为3.63、26.27、33.97 g/m3。在剖面2处，当高程为169.97 m时，工况1、工况2、工况3的Cu2+质量浓度分别为2.72、20.72、27.83 g/m3；当高程为200.90 m时，工况1、工况2、工况3的Cu2+质量浓度分别为0.59、16.20、24.08 g/m3。由此可见，在同一深度，经过相同的时间，尾矿坝加高会使重金属污染物的质量浓度增大。通过图10可知，随着时间的推移，在同一高程位置处，污染物Cu2+浓度呈上升趋势。

图8 重金属Cu2+迁移2 920 d路径分布Fig.8 Path distribution of heavy metal Cu2+ after migration for 2 920 days

图9 垂直剖面重金属Cu2+浓度分布(2 920 d)Fig.9 Concentration distribution of heavy metal Cu2+ in vertical profile(2 920 days)

图10 不同迁移时间下Cu2+浓度分布(剖面1)Fig.10 Concentration distribution of Cu2+ after migration for different time(profile1)

图11 各粒子移距变化(2 920 d)Fig.11 Displacement variation of each particle(2 920 days)

由图11和图12可以看出，尾矿坝加高会使污染物颗粒的迁移距离和平均迁移速度增大，且尾矿坝越高，颗粒的迁移距离和平均迁移速度越大。当尾矿坝加高到工况2时，1～6号颗粒在2 920 d 内的移动距离和平均迁移速度比工况1时分别增加了190.35%、84.14%、73.29%、74.95%、75.64%、76.21%；当尾矿坝加高到工况3时，1～6号颗粒在2 920 d内的移动距离和平均迁移速度比工况2时分别增加了57.92%、21.85%、23.57%、24.10%、24.64%、25.20%。