地下水中污染物溶质运移模拟

喻 佳 魏善彪 麻胜坤

（1．江西核工业环境保护中心，江西 南昌 330002；2．江西省地质局实验测试大队，江西 南昌 330002；3．万载红狮环保科技有限公司，江西 宜春 336100）

0 引言

在建设项目运行过程中出现非正常工况，导致废液或废水泄漏后，进入土壤及地下水环境中，会对土壤及地下水环境造成污染。《环境影响评价技术导则 地下水环境》（HJ 610—2016）的评价对象为地面以下饱和含水层中的重力水，未考虑非饱和带对废液或废水的吸附、弥散、降解等作用。该文以某火电厂脱硫废水收集池为例，利用GMS中FEMWATER模块将研究区饱和及非饱和带作为一个整体进行模拟计算，评价废水泄漏后对地下水环境的影响。

1 FEMWATER模拟求解原理

FEMWATER为基于三维有限元的计算模式。有限元法为采用“分片逼近”方法求解水流、溶质运移耦合的偏微分方程。该求解方法的原理如下：假设模拟范围是由很多很小的可是彼此联系的三维的多面体单元组成的，各个单元的顶点称为节点，不同单元之间则通过节点进行联系；先采用变分法或加权剩余法来建立各节点单元的系数矩阵，对各单元进行离散求解，然后将各单元矩阵集合，合并成描述整个求解范围的方程组，同时把边界条件归并到总矩阵中，进而建立整个求解范围的系数矩阵，用FEMWATER中的求解器对方程组进行求解，最后得出结果。

因此，FEMWATER数值模拟的原理为在一定的初始、边界条件下，用三维有限元法对多孔介质中的饱和-非饱带的水流及溶质运移方程进行求解。下面介绍水流和溶质控制方程。

1.1 地下水水流运动控制方程

地下水水流运动控制方程如公式（1）所示[1]。

式中：μs为贮水率，1/m；h为水位，m；Kx、Ky、Kz分别为x、y、z方向上的渗透系数，m/d；t为时间，d；W为源汇项，m3/d。

1.1.1 初始条件

式中：h0（x，y，z）为已知水位分布；Ω为模型模拟区域。

1.1.2 边界条件

第一类边界如公式（3）所示。

式中：Γ1为一类边界；h（x，y，z，t）为一类边界上的已知水位函数。

第二类边界如公式（4）所示。

式中：Γ2为二类边界；k为三维空间上的渗透系数张量；为边界Γ2的外法线方向；q（x，y，z，t）为二类边界上的已知流量函数。

第三类边界如公式（5）所示。

式中：α为已知函数；Γ3为三类边界；k为三维空间上的渗透系数张量；为边界Γ3的外法线方向；q（x，y，z）为三类边界上的已知流量函数。

1.2 溶质运移控制方程

溶质运移控制方程如公式（6）所示[2]。

式中：R为迟滞系数，量纲为1。为介质密度，kg/（dm）3；θ为介质孔隙度，无量纲；c为组分的浓度，g/L；为介质骨架吸附的溶质浓度，g/kg；t为时间，d；x、y、z为空间位置坐标，m；Dij为水动力弥散系数张量，m2/d；vi为地下水渗流速度张量，m/d；W为水流的源和汇，1/d；Cs为组分的浓度，g/L；λ1为溶解相一级反应速率，1/d；λ2为吸附相反应速率，1/d。

1.2.1 初始条件

式中：C0（x，y，z）为已知质量浓度分布；Ω为模型模拟区域。

1.2.2 边界条件

第一类边界为Dirichlet条件，即指定边界浓度；第二类边界为Neumann条件，即指定边界上的浓度梯度；第三类边界为Cauchy条件，即同时指定边界浓度以及边界上的浓度梯度，也就是第一类边界条件与第二类边界条件的结合。

基于FEMWATER可运用三维有限元法对饱和及非饱和带多孔介质中的水流及溶质运移进行求解，国内不少学者利用其进行研究[3-4]。

2 研究区概况

该文以某火力发电厂脱硫废水收集池为例进行研究。根据研究区工程勘察报告，其地层自上而下为杂填土，层厚1.00m～2.80m；粉质黏土，层厚约1.20 m～4.10 m，局部夹薄层灰黄色粉土、粉砂；淤泥质黏土，层厚0.90 m～3.50 m；粉细砂，层厚约35.0 m；中砂，层厚约7.10 m～13.50 m；泥质砂岩，此层未揭穿。地下水类型主要包括松散岩类孔隙水及红层孔隙裂隙水。松散岩类孔隙水由第四系冲积成因的粉质黏土、粉细砂、砂砾石及卵砾石组成，由于表层发育一定厚度的淤泥质和粉砂质黏土层，使其潜水含水层具有微承压性，含水层厚度一般为30m～40m；红层孔隙裂隙水由于岩石泥质含量较高，抗风化能力差，裂隙多被充填，因此地下水富水性为贫乏。

3 数值模型

3.1 水文地质概念模型

水文地质概念模型为将研究区域地下含水系统的边界、内部的含水层结构、地下水水动力和化学特征、渗透系数等水文地质参数的空间分布以及补径排条件等概化为便于转换为数值模型的基本模式。

该文根据研究区地下含水系统特征建立水文地质概念模型。首先，确定研究区模拟边界。边界条件概化为研究区位于平原区，水文地质条件较简单，取西北侧等水位线为西北边界，该边界为透水边界，反映上游地下水对本研究区域的补给；东南方向的长江划定为排泄边界，作为第一类边界；两侧边界垂直于等水位线，划定为零通量边界；取模型表层为上边界，接受大气降水的入渗补给，取为第一类边界；模型底部为泥质砂岩，该边界划为零通量边界。源汇项：研究区内地下水主要是接受大气降水的入渗补给，主要排泄方式为流入长江，伴有少部分的潜水以蒸发方式排泄。研究区域水文地质较简单，在水平方向上不进行参数分区；根据水文地质剖面图及其显示的地下水含水层岩性分布情况，模型概化为潜水-微承压含水层，垂向上分为3层，自上而下0 m～5 m为杂填土粉质黏土、淤泥质黏土，5 m～45 m为粉细砂、中砂，45 m～53 m为泥质砂岩。

3.2 用TIN生成三维网格

GMS为FEMWATER提供了用TIN生成三维网格、由实体（Solid）模型产生三维网格、由节点产生三维网格、直接由Modflow有限差分网格转为三维有限元网格，共4种可生成三维有限元网格的方式。该研究采用TIN生成三维网格方法。

TIN生成三维网格利用二维网格（2-Dmesh）模块及三角形不规则网格剖分（TIN）模块一起生成三维网格（3-Dmesh），该方法可相对比较方便及快速自动生成大范围的三维网格，局部网格还可以根据需要使用GMS提供的工具进行加密等处理。由于二维网格（2-Dmesh）是三维网格（3-Dmesh）在x-y平面上的投影，因此，需要先生成二维网格（2-Dmesh）。首先，在确定好的多边形（研究区）上建立二维网格（2-Dmesh），然后再将二维网格（2-Dmesh）转换为一对TIN，这一对TIN在x-y平面上必须保持一致，否则生成三维网格（3-Dmesh）会出错，无法生成三维网格。可以通过数据点的高程差值或者输入一个确定的高程值来确定TIN的高程值，这一对TIN中间的区域就是在这一层上所要建立的三维网格（3-Dmesh）区域。然后选中这一对TIN，使用TINS菜单中的FILL BetweenTINS-＞3-Dmesh就可以生成所要的三维网格[5]。

研究区三维网格生成情况如图1所示。

图1 研究区平面、垂向网格剖分图

由图1可知，创建网格时对局部研究区域进行加密，可使预测结果更精确。

3.3 水流数值模型识别与校正

概化的水文地质概念模型，通过FEMWATER界面转化的水流数值模型，须反映研究区域实际流场情况。因此，进行水质模拟预测前，必须对水流模型进行识别、校正，主要是对其参数以及边界条件等进行校正，源汇项包括降水及蒸发等。本次模拟使用的是多孔介质模型，将研究区概化为非均质、各项同性、空间三维结构、非稳定地下水流。根据水文地质资料，初步选定的参数较客观地反映研究区实际的流场情况，加之细致地调参进行模拟后，校正后的模型取得了符合要求的效果。研究区参数见表1。

表1 模型计算主要参数一览表

3.4 水质模型条件的概化

水流模型经识别校正满足要求后，进行溶质模拟前，需要输入水质模型相关参数。通常，以评估潜在的污染源为目的的许多溶质迁移问题都将地下水中污染物浓度设为0；因此，本次水质模型将地下水中氟化物初始溶度设为0。模型将污染源以面源形式设定浓度边界，污染源位置按实际设计概化，并在污染源位置处对网格进行加密处理。在模拟污染物扩散时，不考虑吸附作用、化学反应等因素，重点考虑了对流、弥散作用。

3.5 预测情景

模拟情景设置为在有防渗条件下，脱硫废水收集池防渗破损5%发生泄漏情景下污染物运移模拟。

3.6 预测源强

由于建立的模型没有考虑吸附作用及化学反应等因素，因此在其他条件（包括水动力条件、弥散及泄漏量等）相同的情况下，污染物在地下水中的扩散主要由污染物初始泄漏浓度决定。污染物浓度的高低是相对的，该文以超标倍数这个参数选择预测因子。根据脱硫废水收集池中污染物超标倍数情况，该文选取氟化物为预测因子。在防渗破损情况下，脱硫废水收集池氟化物浓度取50 mg/L。

3.7 泄漏时间确定

脱硫废水收集池下游约30 m处设置了观测井，假定废水持续泄漏，观测井中监测到地下水中氟化物浓度约在273天出现超标，结合监测井每逢单月监测一次的监测频率，确定泄漏时间为330天。

3.8 预测时间

根据电厂运行年限，确定预测时间为7300天。预测污染物泄漏后第100天、第860天、1000天、7300天氟化物在地下水中运移情况。

4 预测结果

将校正后的水流模型带入水质模型进行模拟计算，得到氟化物在地下水中运移的计算结果（图2～图3），以下各图说明了废水泄漏第100天、860天后氟化物在地下水中水平及垂向上的运移情况。氟化物污染物模拟期内运移距离及浓度随时间变化见表2。

图2 废水发生泄漏后，第100天地下水中氟化物浓度分布

图3 废水发生泄漏后，第860天地下水中氟化物浓度分布

表2 不同时间地下水中氟化物运移距离及浓度情况

由图2可知，氟化物发生泄漏后，由于废水持续泄漏，氟化物在地下水对流及弥散作用下，向各个方向扩散，泄漏面源中心点处氟化物浓度最大，为43.91 mg/L，中心点水平运移距离为0 m。

由图3可知，废水泄漏330天后不再发生泄漏，运移至地下水中氟化物为仅有的污染源在地下水对流及弥散作用下，向各个方向扩散，由于不再有泄漏源强进入地下水中，污染羽中心点在地下水对流的作用下向下游迁移，迁移距离为6.38 m，中心点处氟化物浓度最大，为1.0 mg/L，将至标准限值；垂直方向上，氟化物对承压层地下水影响相对较小。由于氟化物浓度将至限值1.0 mg/L，对环境影响很小，因此，未对第1000天、7300天进行评价。根据模型计算结果可知，自860天后，地下水中氟化物浓度逐渐降低，直至降为0。

分析不同时间氟化物在地下水中运移范围可知，在平面上氟化物主要向地下水下游扩散，范围相对较小，在水动力条件及弥散作用下，废水泄漏后地下水中氟化物的迁移对潜水、微承压含水层地的影响在厂区内。

5 结论

该文使用FEMWATER模块对研究区域进行三维有限元的数值模拟计算，利用GMS软件的图像处理功能，将模拟结果可视化输出，操作直观简便。

该文利用FEMWATER对研究区饱和-非饱和带作为一个整体进行模拟计算，通过概化的模型，基于水流控制方程、溶质控制方程的耦合，可视化的输出结果较为完整地呈现了污染物泄漏后，在饱和-非饱和带中的迁移情况。

建议将土壤环境与地下水环境作为一个整体，结合地下水监控井确定废液或废水泄漏时间，评价污染物对土壤及地下水环境的影响。