磷石膏氟化物在包气带–含水层系统中的迁移特征研究

Migration Characteristics of Phosphogypsum Fluoride in Vadose Zone-Aquifer

摘要: 为研究污染物在包气带–含水层中的迁移特征，以磷石膏渗滤液中氟化物为例，通过建立包气带–含水层污染物迁移耦合模型，对磷石膏堆存不同时期污染物的迁移规律进行模拟分析。结果表明：包气带对污染物的垂向迁移产生不同的阻滞效应，有效降低了初期进入含水层的污染物浓度；雨季含水层中污染物水平迁移方向变化复杂，污染物迁移速率更快，污染范围更大，污染管控治理难度高；旱季污染物迁移方向相对单一，表现出相对慢的迁移速率和小的扩散范围，但含水层中污染物累积程度更高，导致污染物浓度也更高，不利于污染物的扩散，污染管控治理难度相对小。研究结果对于包气带–含水层系统水污染治理和风险管控提供一定依据。

Abstract: To investigate the migration characteristics of pollutants in the vadose zone-aquifer system, this paper takes fluoride in phosphogypsum leachate as the research object. A coupled numerical model for pollutant migration in the vadose zone-aquifer system was established to simulate and analyze the migration laws of pollutants during different stacking periods of phosphogypsum. The results show that: the vadose zone exerts varying retardation effects on the vertical migration of pollutants, which effectively reduces the concentration of pollutants entering the aquifer in the initial stage; in the rainy season, the horizontal migration direction of pollutants in the aquifer is complex, accompanied by faster migration rate and larger pollution scope, leading to high difficulty in pollution control and remediation; in the dry season, the migration direction of pollutants is relatively single, showing a slower migration rate and smaller diffusion range, whereas the cumulative degree of pollutants in the aquifer is higher, resulting in higher pollutant concentration, which is not conducive to pollutant diffusion, but the difficulty of pollution control and remediation is relatively lower. The research findings provide a scientific basis for water pollution control and risk management in the vadose zone-aquifer system.

1. 引言

水资源是人类社会赖以生存和发展的基础之一[1]。近年来，我国环境问题日趋严重，水污染问题受到国内外学者的广泛关注[2] [3]，地下水资源是当前人类生活和工业生产的重要资源，很多区域地下水资源的开发利用途径都是饮用水源，污染源和地下水源之间的主要传输介质是包气带及含水层，研究包气带–含水层中污染物迁移成为了地下水污染风险评估的关键[4]。地下水和地表水之间存在频繁的水量交换，而地下水和地表水中污染物的迁移在水量交换过程中存在密切关联[5]，尤其是旱季时地表水反向补给地下水，对含水层中污染物的迁移规律产生显著影响。

云南省磷化工产业规模、磷肥产量均居全国前列。根据相关数据，每生产1吨磷酸会产生4~5吨的磷石膏[6]。磷石膏中含有氟化物、磷酸盐等污染物，磷石膏大量堆存会对周边产生很大的环境污染隐患，磷石膏淋滤液入渗进入包气带–地下水含水层后，在包气带–含水层水量交换的同时，以水作为载体的氟化物、磷酸盐等污染物的迁移规律与特征也发生显著变化，同时污染物在通过非饱和带向地下水迁移渗透过程中会受到河流水位季节性变化的显著影响[7]。

近年来，国内外学者针对包气带防污性能、污染物迁移规律开展了研究[8]-[10]，同时也有学者针对地下水–地表水系统对于农业生产、抗旱、气候等方面开展了大量研究[11]-[13]，但均未充分考虑包气带–含水层的耦合过程。基于Hydrus-1D建立垂向一维污染物迁移模型可有效模拟污染物在包气带中的迁移过程[14]，而饱和带含水层污染迁移规律数值模拟应用最为广泛的数值模拟方法主要为MT3DMS [15]。将Hydrus-1D包气带污染物迁移模型与MT3DMS含水层污染物迁移模型相耦合，以磷石膏淋滤液氟化物作为研究对象，将包气带污染物迁移模型的输出值作为MT3DMS含水层污染物迁移模型的输入值，能够有效获取包气带–含水层中污染物迁移规律，可为后续包气带–含水层系统水污染治理和风险管控提供一定依据。

2. 研究对象与方法

2.1. 研究区概况

2.1.1. 地形地貌

研究区位于滇中高原东缘象山山脉中段，呈小范围缓丘–中、低山地貌、山麓和冲洪积台地–坡谷。地形大体东南高、北西低。盆地标高1819~1900 m，东南部山区高程在2000~2166 m，盆地和山地的相对高差在300 m以上，地形坡度一般5˚~35˚。研究区最低为沙龙河北侧河段，沙龙河为区域最低侵蚀基准面。研究区磷石膏渣场位于东南部山区，利用地形依山而建，渣场总体地势南东–北西走向，东南高西北低，南侧和北侧分别发育一条季节性冲沟，旱季干涸，雨季形成临时地表径流汇入西侧沙龙河。研究区地形地貌图见图1。

Figure 1. Topographic map of the study area

图1. 研究区地形地貌图

2.1.2. 水文条件概况

研究区年平均降水量920 mm，沙龙河自西南侧边界流经至东北侧边界，5~10月份为丰水期，流量占年平均流量的85%，11月至次年4月为枯水期，流量仅占年平均流量的15%。枯季地表水对地下水的补给特征明显。

研究区包气带主要土壤类型自上至下依次为粉质粘土、粘土，粉质粘土厚4~5 cm，粘土层厚12~15 cm。地下水含水层主要为震旦系下统澄江组(Z₁c)砂岩、粉砂质泥岩的强风化层，岩石风化裂隙发育，包气带中无稳定地下水位。研究区东南侧为前震旦系美党组(Ptm)变质砂岩、板岩的强风化层。研究区包气带主要接受大气降水的直接补给，而潜水含水层则接受上覆包气带入渗和雨季东南侧美党组裂隙含水层的侧向补给，以地下岩层裂隙和层面流形式沿地势径流向西侧和北侧沙龙河排泄。

2.2. 研究方法

2.2.1. 模型简介

包气带中污染物迁移利用Hydrus-1D建立一维污染物迁移模型，Hydrus-1D软件使用修改后Richards方程分析水在包气带的变化，只考虑垂向上的一维流动，忽略水平流动和侧向流动。模型中采用经典对流–弥散方程描述一维溶质运移[16]。

$\frac{\partial \left(\theta c\right)}{\partial t}+\rho \frac{\partial s}{\partial t}=\frac{\partial }{\partial x}\left(\theta D\frac{\partial c}{\partial z}\right)-\frac{\partial \left(qc\right)}{\partial x}-\phi$ (1)

式中：c为溶质浓度；D为饱和–非饱和水动力弥散系数；s为吸附在土壤颗粒上的固态溶质浓度；q为体积流动通量密度；$\phi$ 为源汇项。

含水层平面二维模型的建立主要利用MODFLOW和MT3DMS计算机程序，MODFLOW用来模拟地下水在多孔介质中的流动[17]，求解多孔介质地下水流动问题采用单元格中心有限差分算法，这种方法是基于共轭梯度法预处理提出的，其实质是将多孔介质划分一系列网格，利用有限差分方法近似估算单元格中心的水头值，求出椭圆型偏微分方程(PDE)的近似解，得到单元格有限差分矩阵，偏微分方程为：

$-\nabla \cdot K\nabla p=W$ (2)

式中：K为渗透系数张量，单位(L/T)，其中L为长度，T为时间；p为压力水头，单位(L)；W表示源汇项，单位(T⁻¹)，即流入或流出单位体积含水层的流量。

MT3DMS是一个污染物迁移模型程序，它可以模拟多种污染物在地下水系统中对流、弥散、化学反应等作用下的迁移[18]。其建立污染物迁移模型的控制方程为：

$\frac{\partial \left(\theta C^k\right)}{\partial t}=\frac{\partial }{\partial x_i}\left(\theta D_{ij}\frac{\partial C^k}{\partial x_i}\right)-\frac{\partial }{\partial x_i}\left(\theta v_i{C}^k\right)+q_s{C}_s^k+{\displaystyle \sum R_n}$ (3)

式中：θ为含水层有效孔隙度，无量纲；C_k为预测污染物浓度(ML⁻³)；t为污染物迁移时间(T)；x_i,j为沿着i，j方向的坐标(L)；D_ij为水动力弥散系数张量(L²T⁻¹)；v_i为水流速度(L^T−1)；q_s为含水层单位体积流量(T⁻¹)，正值代表源项，负值代表汇项；C_sk为污染源污染物浓度(ML⁻³)；${\displaystyle \sum R_n}$ 为化学反应项(ML⁻³ T⁻¹)。

2.2.2. 构建包气带–含水层耦合模型

渣库磷石膏堆存历史超过20年，根据不同土体厚度，假定包气带水平方向分布均匀连续，由于渣库稳定地下水位埋深较深，雨季时稳定地下水位不会上升至包气带，地下水位波动对包气带模拟厚度影响不大，因此利用Hydrus-1D建立垂向一维污染物迁移模型，粉质粘土和粘土层厚度分别概化为5 cm和15 cm。相关参数选取模型内置参数，分别在深度3 cm、5 cm、10 cm、15 cm处设置观测点。

含水层边界条件和参数选取主要是基于已有调查数据，将研究区边界结合调查结果，圈定模型建立范围。研究区西侧和北侧为沙龙河，旱、雨季概化为定水头边界；而东侧为山脊形成的地表水分水岭，旱、雨季概化为隔水边界；东南侧美党组(Ptm)裂隙含水层旱季对研究区地下水不产生补给，概化为零流量边界，雨季对研究区地下水产生补给概化为定流量边界；另外研究区渣场北侧和南侧季节性冲沟雨季概化为排水沟，旱季不定义边界条件。在MODFLOW程序中定义水文地质参数和边界条件等，分别模拟生成研究区含水层雨季和旱季渗流场，获取研究区含水层地下水渗流特征，为污染物迁移模型的建立提供渗流环境，建立的研究区含水层平面二维模型如图2。将包气带垂向一维模型与含水层平面二维模型耦合后最终建立的包气带–含水层耦合模型如图3。

Figure 2. Two-dimensional horizontal aquifer model of the study area

图2. 研究区含水层水平二维模型

参数数据均通过现场试验、现场观测综合相关文献资料获取，澄江组(Z₁c)裂隙含水层渗透系数通过6个水文试验孔抽水试验获取，每组抽水试验进行3个不同降深的抽水，得出平均值后即为试验孔含水层渗透系数值，6个水文试验孔渗透系数为0.16~0.97 m/d，本次研究按最不利情况取最大值。污染物氟化物浓度通过对渣库中堆渣采集样品进行水平振荡试验，将堆渣样品置于水中水平震荡后取上清液，检测溶于水中氟化物的浓度，7个样品氟化物浓度范围为2.6~55 mg/L，包气带垂向一维模型氟化物源强遵循最不利原则选取最大值，含水层平面二维模型采用包气带污染物迁移模型动态输出值，具体的参数取值如下表1。

2.2.3. 模型验证

在模型建立过程中，不断进行模型参数的调整，使研究区内观测水位和模拟水位动态一致，模拟的地下水渗流场方能反映出研究区真实的地下水流场空间分布特征。通过调查研究区16口地下水井，分别获取雨季和旱季的地下水水位观测值，与数值模型模拟水位值进行拟合，拟合结果如图4，雨季和旱季的

Figure 3. Vadose zone-aquifer coupling model

图3. 包气带–含水层耦合模型

Table 1. Model parameter value table

表1. 模型参数取值表

模型	参数类型		数值	数据来源
包气带污染迁移模型	污染物浓度	氟化物/(mg·L^−^1)	55	实验
含水层地下水渗流模型	裂隙含水层	渗透系数/(m·d^−^1)	0.97	实验
		有效孔隙度/(无量纲)	0.1	砂岩经验值
		侧向补给流量/(m^3·d^−^1)	50	历史数据
	沙龙河	上游初始水头/m	1865	现场观测
		下游初始水头/m	1820	现场观测
	排水沟	南排水沟流量/(m^2·d^−^1·m^−^1)	0.5	现场观测
		北排水沟流量/(m^2·d^−^1·m^−^1)	1.5	现场观测
	年降雨量/(m·d^−^1)		0.0025	历史气象资料
含水层污染迁移模型	污染物浓度	氟化物/(mg·L^−^1)	包气带污染物迁移模型动态输出值	模拟值
	纵向弥散系数/(m^2·d^−^1)		10	千米尺度经验值

Figure 4. Water level fitting curve in the study area

图4. 研究区水位拟合曲线

水位拟合方差分别为0.9837和0.9648，模拟值和观测值具有较高的一致性，拟合精度较高，因此建立的模型可以有效模拟研究区地下水渗流场情况。

3. 结果与讨论

3.1. 含水层渗流特征

图5为数值模拟形成的研究区含水层雨季和旱季的渗流模型。从雨季含水层渗流模型(a)可以看出，雨季研究区含水层地下水水位范围为1820~1883.55 m，受大气降雨垂直入渗和东南侧相邻含水层的补给

Figure 5. Aquifer seepage model in the study area

图5. 研究区含水层渗流模型

作用，水位最高值位于研究区东南侧含水层界线处，水位最低值位于研究区东北侧沙龙河沿岸。地下水等水位线基本与沙龙河平行分布，研究区地下水整体上自东南向西北径流，对沙龙河地表水形成大量的补给。另外受渣场南北两条季节性河流(排水沟)的影响，沿排水沟两侧地下水水位出现明显滞后效应，两侧一定范围含水层地下水向排水沟补给形成临时地表水体。

从旱季含水层渗流模型(b)可以看出，旱季研究区含水层地下水水位值范围为1820~1864.98 m，地下水含水层由于缺乏相邻含水层和大气降雨垂直入渗补给，水位最高值主要分布在研究区西南侧沙龙河沿岸区域，水位最低值仍位于研究区东北侧沙龙河沿岸。此时研究区地下水等水位线基本与西侧沙龙河相交，研究区地下水整体上向北和东北方向径流，仅对沙龙河北侧河段产生一定补给，而沙龙河西南段所在区域水位较高，对研究区地下水含水层形成了明显的反向补给。

通过对比雨季和旱季含水层渗流模型，研究区地下水水位在雨季和旱季时表现出明显差异，旱季时地下水水位明显低于雨季，且水位分布区域产生明显差异。另外，由于旱季地表水对地下水含水层的反向补给，研究区旱季地下水流向也产生了较为明显的改变。

为进一步对比分析研究区水位季节性变化特征，分别定义A-B剖面为旱季渣场下游地下水径流主方向，C-D剖面为雨季地下水径流主方向，分别获取了两个径流方向的水位分布情况，如图6所示。从图中可以看出，雨季时渣场所在区域地下水水位沿A-B和C-D剖面方向均呈逐步下降趋势，越靠近上游，地下水位的差异越大，越靠近沙龙河，地下水位差异越小。在南侧1.4 km排水沟处水位出现明显波动，下降趋势突然变缓，区域地下水虽然整体向西北方向径流，受南侧排水沟影响，部分地下水呈现向西南径流的趋势，且南侧排水沟明显改变了所在区域的地下水位分布，而北侧排水沟上游对地下水分布影响较小，剖面位置并未出现明显的地下水位波动。旱季时渣场所在区域地下水水位沿A-B剖面方向呈逐步下降趋势，沿C-D剖面方向则呈逐步上升趋势，说明区域地下水水位分布和地下水流向受季节性影响明显，区域地下水流场出现了明显的改变。

Figure 6. Water level change map of aquifer profile in the study area

图6. 研究区含水层剖面水位变化图

3.2. 迁移特征及影响分析

3.2.1. 包气带垂向迁移特征及影响因素

包气带垂向迁移模拟时间为10年，时间步长为1年，共选取了4个不同埋深获取了氟化物的浓度穿透曲线。从图7可以看出，随着埋深的逐渐增大，污染物的穿透迁移时间分别约为3.4年、4.8年、8.2年、9.4年，穿透迁移速率逐渐降低，主要是因为随着深度的增加污染物的穿透路径增大。0~5 cm包气带土层性质为粉质粘土，在埋深3 cm和5 cm处，氟化物浓度迅速升高于3年左右时间即接近入渗源强；5~20 cm包气带土层性质为粘土，10 cm和20 cm处氟化物浓度增长速度明显变缓，至8年左右接近入渗源强，主要是因为包气带土体性质发生变化，粘土可以有效阻滞污染物的迁移过程。

Figure 7. Penetration curve of fluoride concentration at different burial depths in vadose zone

图7. 包气带不同埋深氟化物浓度穿透曲线

不同时间节点氟化物在包气带垂向浓度分布情况如图8。从图中可以看出，在0~5 cm范围内，氟化物浓度呈近似直线分布，表现出氟化物在粉质粘土层中迅速迁移穿透的过程。随着时间的推移，5~20 cm范围内氟化物的浓度分布特征逐步由曲线趋于直线，表现出氟化物在粘土层中的缓慢迁移逐渐累积的过程。尤其是在包气带底部，模拟期氟化物的浓度为1.44~54.62 mg/L，第1年浓度增量较小，约为1.44 mg/L，第2年浓度增量迅速增大，约为13.95 mg/L，至第3年浓度增量趋于稳定，约为15.68 mg/L，反映出污染物质首先缓慢迁移穿透，进而快速迁移的过程；而3年以后浓度增长速度逐渐变慢，则反映了污染物在迁移穿透过程中在粘土层中不断累积影响迁移速度。当底部接近源强浓度后，粘土层的阻滞作用基本消失。

Figure 8. Vertical distribution curve of fluoride concentration in vadose zone at different time nodes

图8. 不同时间节点包气带垂向氟化物浓度分布曲线

3.2.2. 包气带–含水层迁移特征及影响因素

利用MT3DMS程序将含水层平面二维模型与包气带一维垂向迁移模型相耦合，将包气带污染物动态输出结果作为含水层污染物输入源强，建立研究区雨季和旱季包气带–含水层污染物迁移模型，分别得到了研究区污染物在雨季和旱季不同时段的迁移结果如图9和图10所示，进而可以分析污染物的迁移规律、特征以及相关影响因素。

从污染物的浓度分布来看，在包气带的阻滞作用下，至第8年含水层中污染物浓度仍未达到淋滤液氟化物浓度源强，堆场渗漏的污染物在水平方向上虽然出现不同程度的迁移扩散，包气带对于污染物在含水层的迁移方向基本不会造成影响，但对污染物浓度影响较大。另外，磷石膏主要成分为硫酸钙，其渗滤液中含有大量钙离子，旱季地下水交换速率慢，渗滤液中钙离子与氟离子进入地下水后易发生沉淀反应生成氟化钙，使可溶态氟化物的浓度降低，造成预测浓度偏高；雨季地下水交换速率快，渗滤液中钙离子与氟离子进入地下水后易发生溶解反应，释放氟离子，使可溶态氟化物的浓度升高，造成预测浓度偏低。但沉淀溶解反应仅对污染物浓度造成一定影响，基本不会影响污染物的迁移速率和迁移距离。

从污染物的迁移方向来看，受不同时期补给条件的影响，雨季和旱季氟化物的迁移主方向明显不同。雨季时，污染物主要向西南–西北方向迁移，而且显现出初期主要向西南往南侧排水沟方向(C-D剖面方向)迁移，后期转为快速向西北往沙龙河方向迁移的趋势；旱季污染物的迁移主方向则变为北北西方向(A-B剖面方向)，迁移主方向较雨季出现明显变化。两个时期地下水和地表水之间补排关系的变化造成了地下水流场的分布差异，是直接导致污染物迁移主方向发生变化的主要因素。

从污染羽的形状和大小来看，在初期污染羽形状均呈现近矩形，雨季受南侧排水沟影响西南方向略突出，雨季第8年后污染羽顺地下水流方向逐渐变宽，前缘收窄，西南方向略突出，但至排水沟为界未继续向外扩散；旱季第8年后污染羽形状则呈现典型的羽毛状趋势，沿地下水流主方向向前延伸，延伸宽度逐渐降低。雨季含水层受大气降雨补给影响，地下水在水平方向径流更分散，污染物呈现出明显比旱季扩散范围更大的趋势，而排水沟则阻断了污染羽的进一步扩散。旱季时，由于含水层中地下水的

Figure 9. Pollutant migration map in rainy season

图9. 雨季污染物迁移图

Figure 10. Pollutant migration map in dry season

图10. 旱季污染物迁移图

径流主要受含水层渗透系数等水文参数影响，影响因素更单一，污染物的迁移主方向也基本不变。整体来说，雨季污染物迁移形成的污染羽范围更大，形状变化差异更大；旱季污染羽形状变化较小，但范围也更小。

从污染物的迁移速率来看，雨季约1600 d左右污染物迁移至南侧排水沟边1#观测点，最大迁移距离约为440 m；3200 d左右污染物迁移至西北侧沙龙河2#观测点，最大迁移距离约为1120 m，污染物的平均迁移速率分别为0.28 m/d和0.35 m/d。旱季2000 d左右污染物迁移至下游3#观测点，最大迁移距离约为240 m，5000 d左右污染物迁移至下游4#观测点，最大迁移距离约为650 m，污染物的平均迁移速率分别为0.12 m/d和0.13 m/d，污染物未迁移至地表水体。从以上数据可以看出，雨季污染物的迁移速率逐渐加快，主要是由于污染物的迁移主方向先向西南排水沟方向，再转为变成向西沙龙河方向，排水沟在一定程度上导致附近地下水位的波动，小范围改变了地下水流场，与西北方向的主流场互相作用，短暂地迟滞了污染物的迁移。旱季污染物基本上只沿一个主方向迁移，迁移速率基本相同，主要是因为旱季地下水流场稳定，对污染物迁移速率的影响较小。整体来说，雨季污染物迁移速率明显大于旱季，大气降雨的垂直入渗补给和相邻含水层的侧向补给明显增大了污染物向地表水体的迁移速率。

为进一步分析污染物的迁移规律和特征，分别获取了4个观测点污染物浓度穿透曲线，如图11所示，图中C1、C2、C3、C4分别代表1~4#观测点处污染物浓度。从雨季氟化物浓度穿透曲线可以看出，渣场污染物泄漏后约2050 d和4000 d后南侧排水沟和西北侧沙龙河地表水中氟化物浓度超过《地表水环境质量标准》(GB3838-2002)中Ⅲ类水质标准。1#观测点中污染物浓度变化速率明显大于2#观测点。从旱季氟化物浓度穿透曲线可以看出，渣场污染物泄漏后约2300 d和6200 d后3#观测点和4#观测点地下水中氟化物浓度超过《地下水质量标准》(GB/T14848-2017)中Ⅲ类水质标准，3#观测点中污染物浓度变化速率大于4观测点。通过对比2个时期氟化物浓度穿透曲线，雨季污染物渗漏以后快速迁移扩散，迁移速率分别为0.28 m/d和0.35 m/d，明显大于旱季(0.12 m/d和0.13 m/d)；而在模拟期内雨季污染物浓度的增长幅度远小于旱季。

Figure 11. Pollutant concentration breakthrough curve at observation point

图11. 观测点污染物浓度穿透曲线

4. 结论

(1) 污染物在包气带迁移穿透过程中在粘土层中不断累积影响迁移速度，不同的土体性质对污染物迁移产生不同的阻滞效应，当浓度接近源强后，包气带的阻滞作用基本丧失。污染物首先在包气带中垂向迁移，然后进入含水层主要进行水平迁移，包气带不会影响污染物在含水层中的迁移方向，但会显著降低初期含水层中污染物浓度。

(2) 不同时期含水层中污染物迁移方向差异较大，雨季污染物的迁移主方向为先西南后西北，污染范围更大，受排水沟影响，污染物平均迁移速率先慢后快；旱季污染物的迁移主方向为北西方向，污染范围相对小，污染物迁移速率先快后慢。雨季污染物迁移速率明显大于旱季，但污染物浓度的增长幅度远小于旱季。

(3) 整体来说，雨季有利于含水层污染物的扩散，迁移方向更为复杂，表现出更快的迁移速率、更大的扩散范围，水污染范围大，管控难度高。旱季污染物迁移方向相对单一，表现出相对慢的迁移速率和小的扩散范围，但含水层中污染物累积程度更高，污染物浓度更高，不利于污染物的扩散，水污染管控难度相对小，但对于地下水污染的治理较为不利。

基金项目

国家自然科学基金项目(42267063)。

NOTES

^*通讯作者。