COD污染物在包气带中迁移的HYDRUS-1D模型预测——以某碳酸锂项目为例

任华丽

（永清环保股份有限公司，北京 100035）

包气带指的是地表到潜水位之间的非饱和带，它既是污染物自地表进入含水层的重要媒介，也是地下水免遭污染的良好保护层[1]。地下水环境评价中，对于范围较广或巨厚包气带地区，大量现场取样非常困难甚至难以实现，故需要在参照经验参数的基础上借助计算机进行模拟计算。目前包气带中污染物运移研究的普遍方法是采用Hydrus-1D软件建立水文地质模型，模拟污染物在时间和空间上的动态变化。

Hydrus-1D是由美国农业部、农业研究会、美国盐土改良中心（US Salinity laboratory）等机构在SUMATRA、WORM及SWMI等模型的基础上创建发展而来的，模拟宏观及微观尺度上饱和及非饱和介质中一维水流、溶质、热和二氧化碳运移和反应的软件，用它可以解算在不同边界条件制约下的数学模型。在农业、水利、环境学等领域得到了广泛应用。近年来，许多研究者用Hydrus-1D对污染物在包气带运移进行了研究，并取得一定成果。

本文应用Hydrus-1D软件，主要模拟污染物随污水进入包气带后的迁移转化过程，即水流运动及溶质在运移过程中发生的反应。对研究区包气带物质组成、边界条件及初始条件等做简化，在简化的水文地质模型的基础上开展包气带污染物运移的研究。本研究的特点是，通过数值模拟的方法研究了污染物在包气带中沿垂向向下运移范围及不同位置污染物浓度随时间变化特征。意义在于了解污染物在包气带中的运移规律，为污染物对包气带环境影响的预防、防治及监测提供依据。

1 包气带特征及模型概化

项目所在区域基础下第一岩层为厚度很大的砾砂层，透水但不具备储水条件，为透水不含水的包气带层。根据厂区地层特征，包气带厚度在110m左右，污染物必须入渗穿过很厚的岩层才有可能到达下部含水层中。根据钻孔资料，将包气带概化为一维非均质各向同性的系统。由于废水池为地埋式，深3m，渗漏处距离含水层107m，因此将项目区包气带概化为1层，模拟厚度107m，岩层为砾砂，渗透系数为9.42m/d。

2 建立数学模型

2.1 数学模型的建立与求解

（1）水流模型。包气带水分运移采用Richards方程的修改形式表示。土壤水分运动的数学模型如下：

式中：h—压力水头，m；θ—土壤体积含水率，m3/m3；t—水分运移时间，d；α—水流方向与纵轴夹角，本次为水流一维连续垂向入渗，故α=0；

K（h）—土壤水的非饱和渗透系数函数，可由方程K（h，x）=Ks（x）Kr（x）计算得出，其中，Ks为饱和渗透系数，Kr为相对渗透系数。

初始条件：

上边界条件：

下边界条件：

Hydrus-1D软件中对土壤水力特性的描述提供了5中土壤水力模型，本次评价选用目前使用最广泛的Van Genuchten-Mualem模型计算土壤水利特性参数，且不考虑水流运动的滞后现象。上边界处理为定水头边界；下边界为自由排水边界。

（2）溶质运移方程。溶质运移方程：

式中：c—污染物介质中的浓度，mg/L；D—弥散系数，m2/d；q—渗流速度，m/d；z—沿z轴的距离，m；t—时间变量，d；θ—土壤含水率，%。

初始条件：

边界条件：

a第一类Uirichlet边界条件：

连续点源：

非连续点源：

b第二类 Neumann 零梯度边界条件：

将废水池看做注入的点源，上边界为变通量边界，下边界为零通量梯度边界。

2.2 模型离散

在Hydrus-1D的SoilProfile-GraphiclEditor模块中剖分包气带结构，根据模拟区水文地质钻孔揭露的地层岩性，厂区在垂向上将模拟区剖分为1层，整个模拟层厚度共计约107m，在垂向上对模拟区进行剖分，将包气带按照0.5m一层进行剖分，剖分网格数为214，单个网格厚度约为0.5m。

2.3 水文地质参数

根据地质、水文地质条件的分析，结合拟建厂址区岩土工程勘察、水文地质勘探资料、软件中给出参数组合并结合经验法，获得的水文地质参数作为初始值。

2.3.1 水分特征曲线参数

Hydrus-1D水流模块中的SoilCatalog项包含砂土、粉土、黏土等12种典型土壤介质及其土壤水分特征曲线相关参数，软件还提供神经网络算法预测的方法，输入土壤中砂土、粉土及黏土的百分比估算出土壤层的相关水分特征曲线参数。综合已有参数、预测参数和实测参数，获得模拟区岩层水分特征曲线参数，具体参数详见表1。

表1 土壤水分特征曲线参数值

2.3.2 溶质反应参数

污染物在土壤中的运移受吸附/解析作用影响大部分污染物的吸附符合Freundlish、Langmuir等温吸附规律，且以各种形式存在于土壤环境中的污染物物质，会发生一系列的物理、化学和生化作用，本次评价只考虑溶质在固液相间的线性平衡等温吸附作用，忽略物理、化学和生化作用。根据查阅文献资料给出CODMn污染物的溶质运移参数值分配系数Kd、Langmuir等温吸附系数Nu、Freundlish吸附系数Beta、液相与气象分配系数Henry、降解速率等，具体参数详见表2。

表2 土壤水分特征曲线参数值

3 观测点及时间设置

模型构建时根据相邻最近钻孔的地层资料进行概化，在模型不同深度分布设置深度不同的7个观测点：N1：3m、N2：10m、N3：20m、N4：40m、N5：60m、N6：80m、N7：107m，来分析不同污染深度污染物浓度随时间变化的情况。模型输出时间，分布结算模型运行时间为30天、90天、180天、365天（1年）、1825天（5年）、3650天（10年）、5475天（15年）、7300天（20年）。

4 污染情景设定

4.1 正常状况

正常状况下，为有效防止废水泄露对区域土壤产生影响，进行分区防腐、防渗处理，同时加强对废水输送管道的维护和管理，防止废水的跑、冒、滴、漏和非正常排水。在不受外界因素作用下，按照设备操作规程进行合理、规范操作，废水池不会发生泄漏事故。因此，正常状况下建设项目不会发生渗漏污染土壤的情景。

4.2 非正常状况

由于废水池为地下式水池，废水池发生泄漏，管理人员不能及时发现，致使少量污水通过泄漏部位渗入包气带，最终进入含水层。

废水池为钢筋混凝土构筑物，根据《给水排水构筑物工程施工及验收规范》（GB 50141-2008）中规定钢筋混凝土结构水池渗水量不得超过2L/（m2·d），非正常状况下废水池因老化或者腐蚀产生的渗漏量按正常工况下最大允许渗漏量的10倍考虑，废水池池体底面为水平矩形，容积为216m3（9m×6m×4m），有效贮水深度3m，有效容积为120m3，按涉池壁和池底的浸湿面积计算渗漏量，则渗漏量为2.88m3/d。将发现污染物泄漏并采取措施停止泄漏的时间确定为3个月。根据废水池废水水质监测，预测因子COD浓度取78mg/L。由于地下水评价工作中通常采用耗氧量（CODMn）作为评价指标，根据COD和CODMn的经验关系，认为COD浓度与4倍的耗氧量等效。预测因子源强如下：

CODMn：19.5mg/L×2.88m3/d×90d=5.05kg

表3 污染源强列表

5 污染预测结果

废水池破裂，废水中CODMn持续渗入土壤并逐渐向下运移。初始浓度为19.5mg/L，在土壤剖面不同观测点CODMn浓度随时间变化模拟结果如图1所示，在不同水平年CODMn沿土壤迁移模拟结果如图2。

图1 土壤剖面观测点土壤水中CODMn浓度变化曲线图

（T0：0d，T1：30d，T2：90d，T3：180d，T4：365d，T5：1825d，T6：3650d，T7：5475d，T8：7300d）

由图1～图2土壤模拟结果可知，CODMn在土壤中随时间不断向下迁移，随着入渗水量增加，污染物浓度逐渐增加，当停止泄漏后，流入地下水的水量减少，污染物浓度逐渐降低。渗漏30d后，在包气带深度27.5m处CODMn浓度达到3mg/L（土壤含量0.2mg/kg）；渗漏90d后，在包气带深度69.5m处CODMn浓度达到3mg/L，潜水面CODMn浓度达到检出限值0.05mg/L，潜水面CODMn浓度逐渐增大；第142d，潜水面CODMn浓度达到3mg/L；第225d，潜水面CODMn浓度达到最大值12.7mg/L（土壤含量1.662mg/kg），然后逐渐降低；第476d，潜水面CODMn浓度降低到检出限值0.05mg/L，不再渗出。

6 结论

本文采用Hydrus-1D软件对非正常工况下废水池发生泄漏后COD污染物在包气带中的运移情况进行预测，浓度随时间先增加后降低；距离入渗点越近的位置，污染物浓度越高；距离入渗点越远的位置，污染物浓度越低。在土壤纵剖面上，随着污染物入渗和停止，土壤中污染物的含量增大后减小。在潜水面上，污染物浓度先增加后减小[11]。

本文采用Hydrus-1D软件建立水文地质模型，模拟污染物在时间和空间上的动态变化，具有极强的便利性，为污染物对包气带环境影响的预防、防治及监测提供依据。