中小河流水体重金属Zn、Pb、As沿程迁移扩散过程模拟<br/>——以沘江为例

中小河流水体重金属Zn、Pb、As沿程迁移扩散过程模拟
——以沘江为例

2022-02-21王瑞芳赵筱青李嘉旺窦小东

云南大学学报(自然科学版) 2022年1期

易琦，王瑞芳，赵筱青**，李嘉旺，窦小东

(1.云南大学地球科学学院，云南昆明 650500；2.云南省气象服务中心，云南昆明 650034)

重金属污染物毒性强、不易降解，其浓度达到一定阈值，不仅会危害水体及水生生态系统，还能经过食物链传递，危害人体健康[1-3].河流重金属元素的转化主要为水体中溶解态、悬浮态以及底泥中沉积态三者之间的相互转化[4-6].目前，对于大型河流水体中重金属迁移模拟多采用二、三维模型，中小型河流多采用一、二维模型.黄本生等[7]建立河流重金属随水−悬浮物−底泥迁移转化耦合模型.武万国[8]建立了湘江长株潭河段二维水动力−水质模型.无论哪一种模型，确定水体重金属污染物的迁移转化系数，均是确定水体重金属污染物迁移规律的关键因素.沘江为中小河流，是澜沧江上游的一级支流，也是沿岸居民生产生活的主要水源，但由于近年来重金属对河水的污染，河水使用率急剧降低，导致沿岸的兰坪、云龙两县出现不同程度的水质性缺水[9].因此，摸清沘江水体重金属污染物的迁移转化规律，对沘江水体重金属污染治理具有积极意义.本文在枯水期对沘江水体采样分析基础上，对河水中重金属污染物的综合衰减系数进行率定，从而建立沘江水体重金属污染物的迁移转化方程，研究结果可为沘江流域环境污染防治及生态恢复等提供科学依据.

1 研究区概况

沘江属澜沧江水系，流域面积约2 400 km2，干流长约150 km，总体流向由北向南，途经兰坪、云龙两县后汇入澜沧江[10-11].金顶铅锌矿位于沘江上游河段的东部，其储量居中国已探明铅锌矿储量的首位.矿区三面环山，大部分可以露天开采，雨季时非常容易汇水，导致大量的矿渣、碎石等被冲入河道，沘江水体已鱼虾绝迹，污染严重，无法饮用和灌溉，沘江水质性缺水问题严重[12-14].

2 研究方法

2.1 样品采集与分析2010 年1 月底，按照矿区对沘江影响大小，以距离矿区远近为原则，距离越远样点越稀疏，反之样点越密集，共取得沘江水样15 个（图1）.其中，矿区上游采样点3 个，基本不受矿业活动影响，取其均值以反映矿业开采前的沘江水质状况，作为水体的背景值；矿区下游采样点12 个，其中，1～11 号样点分别位于兰坪县和云龙县境内，12 号样点为澜沧江与沘江交汇处水样.水样于岸边采集并贮存于2.5 L 的高密度聚乙烯瓶中.水样的重金属质量浓度由昆明矿产资源监督检测中心检测分析，检出限分别为：ρ（Zn）<1.0 μg/L，ρ（Pb）<0.01 μg/L，ρ（As）<0.01 μg/L.

图1 研究区域及采样点分布Fig.1 Research area and distribution of sampling points

2.2 污染物迁移扩散模拟方法

2.2.1 水体重金属迁移扩散模型的选取天然中小河流多为浅宽河流，重金属污染物进入水体后很快在垂向均匀混合[6]，因此在考虑河流重金属污染物质量浓度变化时只考虑横向（y）与纵向（x）方向即二维状态即可.而对于天然宽度较窄的小河流而言，可以近似地认为重金属污染物进入水体后迅速在横向（y）与垂向（z）方向均匀混合，其水体中重金属污染物迁移转化仅考虑纵向（x）方向即一维状态即可，沘江即属此类河流.

综合考虑水体中重金属迁移转化的水动力学特性，及重金属在水体−悬浮物−底泥间的吸附、解吸、沉降、再悬浮过程[15]，同时，结合重金属自身的化学转化因素的影响，可知沘江水体重金属污染物迁移的一维模拟方程可用下式表示：

式中，D为分子扩散系数（单位：m2/s）；ρ为溶质的质量浓度（单位：mg/L）；X为河流纵向上距排污口附近距离（单位：m）；u为河水的平均流速（单位：m/s）；K为综合衰减系数（单位：s−1）；ρ0为矿区下游排污口附近重金属质量浓度值（单位：mg/L）；ρb为沘江河流水体重金属背景值（单位：mg/L）.

若不考虑水体重金属背景值的影响，式（1）的解析解为：

若考虑到背景值的影响则公式（2）、（3）可分别表示为：

本次沘江实地采样，1 号采样点是矿区下游水体中重金属与沘江水体重金属背景值的总含量，因此ρ0的值已经包含ρb的影响.因此，对于沘江水体重金属迁移转化的的解析解应采用公式（2）、（3）模拟，鉴于公式（3）中纵向扩散系数D的确定需要考虑河水深度及比降的影响，而本次采样并未测得沘江深度及比降的具体数据，同时河流的比降及深度对重金属迁移的影响也主要是通过流速实现.因此，本次沘江水体重金属迁移转化的解析解采用公式（2）进行.

2.2.2 模型参数的率定及检验模型中边界值的选取方法是取矿区下游的1 号样点的Zn、Pb、As 的质量浓度作为连续面源污染的初始值ρ0.

模型中综合衰减系数K的率定一般包括经验公式法、实验法、试错法[16].经验公式法是根据前人总结的一些经验关系式或近似计算方法确定某一参数的取值.所谓实验法是通过相关实验确定某一参数的具体值，而试错法则是通过查阅相关文献，确定该参数的范围或该参数的取值大概在哪个数值附近，进而通过尝试不同的赋值即试错，并对赋值后的结果检验，若达到预期要求或通过检验，则所赋值即为该参数值，若达不到预期要求或通不过检验，则要继续赋予其它的值，直到取得满意效果或通过检验为止.鉴于目前并无成熟的综合衰减系数K的经验公式，且由于自身实验条件的限制，本文采取试错法来确定K的取值，并运用PASW18.0软件不断对K的不同试错值所生成的各样点模拟值与实测值进行配对样本T检验，直到通过检验进而确定K值的范围，并将配对T检验效果最理想时所用的K值作为K的首选值.

配对T检验是对相同总体在两种不同的条件下进行相同研究，对比两种条件下总体均值是否存在显著差异的方法.因其计算简便、检验性能高效，因而得到广泛应用[17].配对样本T检验的零假设H0认为，两总体均值之间没有显著差异.检验思路为：首先计算每对观察值的差值，得到差值序列；再计算差值的均值；最后检验均值是否和0 存在显著差异.若存在显著差异，则认为两总体均值间有显著差异，反之，没有显著差异.

应用PASW18.0 软件分别进行Zn、Pb、As 的实测值与不同K值下其模拟值的配对样本T检验，并根据T分布表给出T值所对应的概率值P.依据为若P<显著水平α，则拒绝零假设，认为两总体均值有显著差异；若P>显著水平α，则不能拒绝零假设，认为两总体均值不存在显著差异，判定K值范围[18].

2.2.3 基于MATLAB 对沘江水体重金属迁移全程模拟 MATLAB 软件即矩阵实验室，作为科技应用软件，其集程序设计语言、图形处理以及数学于一体[19-20]，把科学计算、结果的可视化和编程集中在一个非常方便的环境中，该软件已经在许多领域得到了广泛的应用.宋新山等[21]对应用该软件构建地表水环境系统的模拟模型做了详细的实例演示，效果较好.因此，本文选用该软件对沘江水体重金属迁移进行全程模拟.

3 结果与分析

3.1 沘江水体中Zn 的模拟结果与实测数据的配对样本T 检验模型中各样点距1 号样点的距离由采样时GPS 所测得的各样点的地理坐标经Arcgis 软件求得.由于天然河流的水文条件及周边环境各不相同，因此目前并无普遍适用的K值范围，但可以确认在有连续稳定污染源排放的天然河流中K>0.本文尝试性选取步长0.001 s−1逐次递增K值，运用公式（2）计算2～12 号采样点Zn 的模拟值.PASW18.0 软件对K值取0.001～0.04 s−1时，2～12 号采样点Zn 的模拟值与实测值的配对T检验结果显示（图2（a））：当K值介于0.013～0.017 s−1时，P>α=0.05，模拟值与实测值不存在显著差异，即当沘江水体中Zn 的综合衰减系数K值介于0.013～0.017 s−1时，公式（2）可以较好地反映沘江水体中Zn 的迁移转化规律；当K=0.015 s−1时，P=0.762，公式（2）的模拟值最接近真实值.因此，沘江水体中Zn 的综合衰减系数K的最优值取0.015 s−1.

3.2 沘江水体中Pb 的模拟结果与实测数据的配对样本T 检验鉴于10～12 号采样点间Pb 的迁移转化受上游拦河坝的影响较大的事实，沘江水体中Pb 的综合衰减系数K值采取分段率定，即分别率定2～9 号采样点与10～12 号采样点的综合衰减系数K值.仍然选取步长0.001 逐次递增K值，并运用公式（2）分别计算2～9 号及10～12 号采样点Pb 的模拟值.

PASW18.0 软件对K值取0.001～0.04 s−1时，2～9 号采样点Pb 的模拟值与实测值的配对T检验结果显示（图2（b））：当K值介于0.003～0.009 s−1时，P>α=0.05，模拟值与实测值不存在显著差异，即当沘江水体2～9 号采样点间Pb 的综合衰减系数K值介于0.003～0.009 s−1时，公式（2）可以较好地反映沘江水体中2～9 号采样点间Pb 的迁移转化规律；当K=0.006 s−1时，P=0.692，公式（2）的模拟值最接近真实值.因此，沘江水体中2～9 号采样点间Pb 的综合衰减系数K的最优值取0.006 s−1.

PASW18.0 软件对K值取0.001～0.02 s−1时，10～12 号采样点Pb 的模拟值与实测值的配对T检验结果显示（图2（c））：当K值介于0.002～0.008 s−1时，P>α=0.05，模拟值与实测值不存在显著差异，即当沘江水体10～12 号采样点间Pb 的综合衰减系数K值介于0.002～0.008 s−1时，公式（2）可以较好地反映沘江水体中10～12 号采样点间Pb 的迁移转化规律；当K=0.004 s−1时，P=0.667，公式（2）的模拟值最接近真实值.因此，沘江水体中10～12 号采样点间Pb 的综合衰减系数K的最优值取0.004 s−1.

图2 不同K 值时2～12 号采样点模拟值与实测值配对T 检验结果Fig.2 PairedT test results of simulated and measured values at samples 2-12 at differentK values

3.3 沘江水体中As 的模拟结果与实测数据的配对样本T 检验10、11、12 号采样点间As 的迁移转化亦受上游拦河坝的影响较大，沘江中As 的综合衰减系数亦采取分段率定，即分别率定2～9 号采样点与10～12 号采样点的综合衰减系数K.同样选取步长0.001 逐次递增K值，并运用公式（2）分别计算2～9 号及10～12 号采样点As 的模拟值.

PASW18.0 软件对K值取0.001～0.04 s−1时，2～9 号采样点As 的模拟值与实测值的配对T检验结果显示（图2（d））：当K值介于0.013～0.027 s−1时，P>α=0.05，模拟值与实测值不存在显著差异，即当沘江水体2～9 号采样点间As 的综合衰减系数K值介于0.013～0.027 s−1时，公式（2）可以较好地反映沘江水体中2～9 号采样点间As 的迁移转化规律；当K=0.018 s−1时，P=0.918，公式（2）的模拟值最接近真实值.因此，沘江水体中2～9 号采样点间Pb 的综合衰减系数K的最优值取0.018 s−1.

PASW18.0 软件对K值取0.001～0.04 s−1时，10～12 号采样点As 的模拟值与实测值的配对T检验结果显示（图2（e））：当K值介于0.001～0.04 s−1时，P>α=0.05，模拟值与实测值不存在显著差异，实际上，当K=0.023 s−1时，10～12 号采样点As 模拟值的数量级已达10−5，已非常接近0，即K>0.023 s−1时，K值的增加对T检验结果已不会形成影响.因此，可将K值限定在0.001～0.023 s−1之间，即当沘江水体10～12 号采样点中As 的综合衰减系数介于0.001～0.023 s−1时，公式（2）可以较好地反映沘江水体中10～12 号采样点间As 的迁移转化规律；当K=0.002 s−1时，P=0.859，公式（2）的模拟值最接近真实值.因此，沘江水体中10～12 号采样点间As 的综合衰减系数K的最优值取0.002 s−1.

综上所述，矿区下游沘江水体中Zn 的综合衰减系数K值介于0.013～0.017 s−1时，其最优值为0.015 s−1.矿区下游至果朗村（9 号采样点）附近江段水体Pb 的综合衰减系数K值介于0.003～0.009 s−1时，其最优值为0.006 s−1，As 的综合衰减系数K值介于0.013～0.027 s−1时，其最优值为0.018 s−1.此外，果朗村附近至沘江与澜沧江交汇江段水体Pb的综合衰减系数K值介于0.002～0.008 s−1时，其最优值为0.004 s−1，As 的综合衰减系数K值介于0.001～0.023 s−1时，其最优值为0.002 s−1.

3.4 沘江全程水体中重金属污染物迁移的沿程模拟在MATLAB 中调用Function 函数，并对ρ(x)=各相关参数进行定义后，即可得到沘江全程水体中重金属污染物沿江迁移的模拟模型，也即矿区下游沘江水体中Zn 的综合衰减系数K=0.015 s−1；矿区下游至果朗村(9 号采样点)江段水体Pb 的K=0.006 s−1，As 的K=0.018 s−1；果朗村9 号采样点至沘江与澜沧江交汇江段Pb 的K=0.004 s−1，As 的K=0.002 s−1.依此模型及拟定的相应K值，只需在矿区下游附近河流水体中某处分别测得该点Zn、Pb、As 的质量浓度值ρ0，并在模型中输入该值，模型即可自动模拟该测点下游水体中任一点的Zn、Pb、As 质量浓度值，并绘制该点下游河流沿程水体中Zn、Pb、As 质量浓度变化曲线.

以本次采样1 号样点作为ρ0，则矿区下游ρ0点以下沿程Zn、Pb、As 的模拟值与本次采样点实测值的对比结果见图3.Zn 的沿程模拟值与本次采样点的实测值极为相近；而Pb 的沿程模拟值与本次采样点实测值的变化趋势亦较为接近.虽然个别样点的实测值与模拟值略有偏差，但下游波动幅度较大的10 号采样点与沿程模拟值（ρ0点下游88 km处）尚且较为吻合，可见Pb 的沿程模拟值基本能客观反映水体中Pb 的质量浓度变化，模拟效果较好；除个别采样点外，矿区下游ρ0点以下As 的沿程模拟值与对应的实测值也基本吻合，即便在下游出现实测值显著波动的10～12 号采样点也能较好反映其变化趋势，模拟计算能够较好地拟合水体中As的质量浓度沿程变化.由此，经本文率定综合衰减系数K值后的沘江水体重金属迁移模型，能够较好地模拟沘江枯水期水体重金属Zn、Pb、As 的质量浓度变化，该模型可用于对沘江枯水期水体重金属Zn、Pb、As 质量浓度值的模拟.

图3 沿程模拟值与金顶矿区下游采样点实测值对比Fig.3 Comparison of simulated values along the way with measured values at downstream sampling points in the Jinding Mining area

4 讨论

沘江受兰坪金顶铅锌矿（铅锌储量居世界第二位、亚洲第一位）矿区影响，水体重金属Zn、Pb、As 污染严重[10]，而沘江流域的生态健康是维系澜沧江上游区域生态安全的重要保障.因此，摸清沘江水体重金属Zn、Pb、As 沿程迁移扩散规律，建立适合沘江水体Zn、Pb、As 迁移转化的模拟模型对流域生态安全等的研究至关重要，建模后准确率定模型中各重金属的综合衰减系数，则是模型能否准确模拟水体中重金属迁移扩散的关键所在.

何绍福等[22]在2015 年研究沙溪沙县段水环境容量时，建立了该河段的一维水质模型，并采用试错法率定了CODCr的耗氧速率、氨氮的硝化速率、无机磷的沉积速率等水质因子的衰减速率，进而采用试错法确定了沙溪沙县段水体中CODCr、氨氮、总磷等的环境容量，取得了较好的效果.徐凌云等[23]在2017 年研究温岭市平原河网水环境容量时，对温岭市平原河网水量及水质进行了一维建模，并采用试错法率定了CODMn、氨氮等的衰减系数，模拟值与监测值拟合较好.本文同样采用试错法率定了沘江流域水体重金属Zn、Pb、As 沿程迁移扩散的综合衰减系数，率定结果也较为理想.因此，在进行水体重金属迁移扩散的综合衰减系数率定时，若实验条件不允许，则采用试错法确定水体重金属迁移扩散的综合衰减系数也是一个较为稳妥的方法，但也要注意试错法所采取的步长选择难免会对K值的准确性产生一定的影响.

由于沘江为天然浅窄的河流，因此本次沘江水体重金属迁移扩散建模采用的是一维模型，即只考虑纵向（x）重金属的微量扩散，并未考虑重金属在横向（y）及垂向（z）的微量扩散，这也会在一定程度上增加研究结果的不确定性.本文未考虑河流沿岸工农业生产对沘江水体造成的影响，有待于今后进一步的研究.

此外，沘江丰水期和枯水期的水体流速等水文参数存在一定的差异，重金属在水体−悬浮物−底泥间的吸附、解吸、沉降、再悬浮过程及重金属污染物的化学转化过程亦会不同.若能够在丰水期对沘江水体再次采样，并与本次采样数据对比分析，模型将会取得更好的适用性.因此，在考虑对沘江丰水期重金属Zn、Pb、As 迁移转化进行模拟时，应在丰水期再次采样，并对模型中Zn、Pb、As 的综合衰减系数K分别在其取值范围内重新率定，或取个别样点进行适用性验证后，再应用本模型进行模拟.

5 结论

根据水体中重金属沿程迁移的普适模型，结合沘江水体重金属的实测数据，本文通过试错法率定了该数学模型的综合衰减系数K的取值范围及其最优值，确定了沘江水体中重金属迁移的数学模型，并运用MATLAB 建立了适合沘江全程的水体重金属迁移转化的模拟模型，得到如下结论.

（1）枯水期沘江流域金顶矿区下游Zn 的综合衰减系数K值介于0.013～0.017 s−1时，其最优值为0.015 s−1；矿区下游至果朗村（2～9 号采样点）附近江段水体Pb 的综合衰减系数K值介于0.003～0.009 s−1时，其最优值为0.006 s−1，As 的综合衰减系数K值介于0.013～0.027 s−1时，其最优值为0.018 s−1；果朗村附近至沘江与澜沧江交汇江段（10～12 号采样点）水体Pb 的综合衰减系数K值介于0.002～0.008 s−1时，其最优值为0.004 s−1，As的综合衰减系数K值介于0.001～0.023 s−1时，其最优值为0.002 s−1.即枯水期沘江流域金顶矿区下游重金属Zn、Pb、As 的综合衰减系数K分别为Zn 全程0.015 s−1；果朗村上游Pb 为0.006 s−1，As为0.018 s−1；下游Pb 为0.004 s−1，As 为0.002 s−1.