李文钰，侯精明，刘占衍，张 松，栾广学，杜颖恩，韩占涛

（1.西安理工大学 西北旱区生态水利国家重点实验室，陕西 西安 710048；2.河北省邯郸水文勘测研究中心，河北 邯郸 056011；3.四川水利职业技术学院，四川 成都 611830；4.河南智河工程技术有限公司，河南 郑州 450003；5.生态环境部 土壤与农业农村生态环境监管技术中心，北京 100020）

随着人类社会的发展，对自然资源的需求不断提升，这一过程中矿产资源的开发利用造成大量矿渣堆积，威胁生态安全［1－3］。如富含多种矿产资源的文峪河流域，因早年废弃矿渣未得到有效处理而致使其在河道堆积，上游来水或降水对其冲刷、淋溶和浸泡后产生了大量高浓度污染物，对下游河道水环境及居民饮用水安全造成了极大威胁。为遏止矿渣堆积对生态环境的不利影响［4－8］，提出高效精准的水环境治理方案，需要定量分析矿渣堆积后流域水系污染物的分布特征及迁移总量。

近年来国内外学者常使用的水环境模型有SWMM、SWAT、MIKE 模型等，如：曾晓岚 等［9］使用SWMM 模型对滇池东岸花卉大棚种植区降雨径流水质、水量进行了模拟，并分析了氮磷的单位面积负荷率及输出特征；刘骞等［10］应用SWAT 模型构建了岷江流域分布式水文和污染负荷模型，并计算了各项减排措施及气象驱动条件对流域水环境改善的贡献率；刘晨辉等［11］采用MIKE 模型模拟了COD、BOD5、氨氮等在长江的迁移和扩散情况。以上模型对特定区域具有较好的模拟效果，但SWMM 等一维模型无法展示污染物的对流扩散过程；SWAT 等水文模型参数较多，建模困难，且根据其本底特征确定污染物输移路径后无法在模拟过程中依照水动力条件自动调整路径，适用于大尺度流域水文模拟，无法对中小尺度流域的水动力和水质变化过程进行高分辨率模拟；MIKE 等水动力水质模型虽然可以进行高分辨率水环境模拟，但其模型重点在于水动力模块，且对于百万数量级以上网格，模型运行速度较慢。基于此，本文采用GPU 加速技术构建了小文峪河流域水环境数值模型，并充分利用现有数据模拟分析小文峪河流域污染物空间分布及污染物来源，以期为小文峪河流域水环境治理提供参考。

1 材料与方法

1.1 研究区概况

文峪河位于湖北省十堰市竹山县西北部，上游主要支流有小文峪河和界岭河，下游汇入堵河。文峪河流域属于亚热带季风气候区，雨量充沛，热量充足。水量多寡取决于气候条件，呈季节性变化。小文峪河流域地表水主要污染源为矿洞和矿渣堆，矿渣堆沿自然坡面、河道沟谷及近沟坡面溜坡放置，降水时水流从矿渣堆表面和内部流过，形成矿渣堆淋溶水，严重影响小文峪河的水环境。小文峪河从上游到下游主要支流有竹园沟、清石沟、孙家沟、宝泉寺及郭家沟（见图1），各支流矿渣堆面积分别为8 402、5 997、951、16 638、20 824 m2，矿渣方量分别为50 423、19 679、3 377、96 788、84 500 m3。本研究选取超标污染物硫酸根离子和铁离子为目标污染物，对研究区污染物时空分布进行高分辨率模拟并量化分析小文峪河污染物来源。

图1 研究区矿渣堆分布及监测点布设

1.2 流量及污染物数据

本研究采用2021 年野外单场次监测数据，监测点15 个，见图1。监测指标主要有硫酸根离子、铁离子、流量等。监测时间为旱季，流量较为稳定，因此以该场次监测数据为模型率定数据。监测点主要分布在矿渣堆上下游及支流与主河道交汇处，因此以支流矿渣堆上游监测水质和流量数据为模型边界条件，认为无矿渣堆积的支流上游来水污染物浓度为0。支流上游来水对部分矿渣进行冲刷与矿洞水混合后经径流过程到达下游，发生复杂的物理变化，采用一级冲刷过程来综合概化地表水流经矿渣堆并与矿洞水混合到达下游的污染物浓度变化规律。一级支流在向主河道汇流过程中会有二级支流的汇入，由于缺乏二级支流监测数据，因此采用一级衰减过程来描述污染物浓度下降过程。小文峪河主河道上游流量和水质边界条件根据监测数据设置，即不考虑研究区域上游污染物空间分布及水量汇集特征。

1.3 数值模型及求解方法

模型的地表水动力模块控制方程为平面二维浅水方程（简称SWEs）［12－14］，主要针对具有自由液面且以平面运动为主的水流，只考虑水平方向流速，忽略垂向流动，且忽略运动黏性项、紊流黏性项、风应力和科氏力。二维浅水方程的守恒格式可用如下矢量形式来表示：

式中：t为时间；qx、qy分别为x、y方向的单宽流量；h为水深；q为变量矢量，包括水深及x、y方向的单宽流量；F、G分别为x、y方向的通量矢量；g为重力加速度；S为源项，包括降水源项、底坡源项及摩阻源项；u、v分别为x、y方向的流速；i为入渗源项；zb为河床底面高程；Cf为谢才系数，Cf＝gn2／h1／3，其中n为曼宁系数。

水质模块控制方程为对流扩散方程：

式中：C为污染物垂向平均浓度；Dx、Dy分别为x、y方向的扩散系数，因对流较强，故忽略扩散作用；qin为点源排放的流量强度；Cin为点源污染物垂向平均浓度；S为污染物浓度对时间的全导数dC／dt，即转化项，如污染物生化反应、生物作用等。

仅考虑衰减与冲刷反应时，转化项S可表示为

式中：K1为衰减系数，K2为冲刷系数。

采用Godunov 格式有限体积法求解二维浅水方程，选用二阶MUSCL 方法进行变量重构［15－16］，通过HLLC 近似黎曼求解器求解控制单元界面上的质量通量与动量通量。为提高计算效率，模型结合CUDA 架构实现GPU 加速［17］，同时保持精确度和高稳定性。

1.4 污染物通量计算

污染物通量是指水体中的污染物在单位时间内通过某一断面的总质量，是水文、地质、化学及生物等综合作用的结果［18－20］。污染物通量计算公式为

式中：L为污染物通量；Q为平均流量；m为转化系数，值为0.086 4［18］。

2 模型构建

地形数据的分辨率为2 m，共有269 437 个网格，模型计算采用开边界，仅在支流监测点位置设置入流流量和污染物初始浓度。研究区处于山区，且河道长期有水，可认为河道水流无下渗；参照文献［21－22］，根据当地实际情况，取曼宁系数为0.04。由于矿洞水仅有极少部分能汇入地表水，且流量较小，因此将其与各支流的发源地进行合并建模。研究区以矿渣堆为主的污染源释放的污染物会进入地表水和地下水，再随地表水和地下水向下游迁移，使下游地表水和地下水中污染物超过相应水质标准。由于矿渣堆积方量不同，因此采用GIS 技术确定各支流矿渣堆积的范围，分别率定水流经过各支流矿渣堆及其上下游时污染物的冲刷系数和衰减系数，当地表水以一定流速通过矿渣堆时，模型将计算出各支流矿渣堆析出的污染物量。

为分析矿渣堆中硫酸根离子和铁离子对流域水质的影响，采用基于GPU 加速技术的水动力水质耦合模型对污染物的迁移转化过程进行模拟。为验证模型模拟结果的准确性，根据郭家沟、宝泉寺、孙家沟、清石沟及竹园沟的水质监测数据，对模型进行参数率定，结果见表1。

表1 水质参数率定结果 10－3／s

3 结果分析

水质模拟结果见表2，可以看出各监测点污染物模拟结果相对误差均在10%以内，表明所构建模型能较好地模拟小文峪河流域硫酸根离子和铁离子在河网中的迁移转化过程。

表2 硫酸根离子和铁离子模拟结果和实测结果对比

3.1 污染物空间分布特征

研究区硫酸根离子和铁离子空间分布见图2、图3，可以看出矿渣堆上游污染物质量浓度均较低，经过矿渣堆途中污染物质量浓度逐渐上升，在矿渣堆下游污染物质量浓度达到峰值。竹园沟、清石沟、孙家沟、宝泉寺、郭家沟5 条支流硫酸根离子质量浓度分别从26.0、23.5、134.0、19.5、20.3 mg／L 上升至222.0、121.1、170.0、453.0、409.0 mg／L，铁离子质量浓度分别从0.13、0.07、6.94、0.06、0.15 mg／L 上升至10.60、28.90、10.50、83.10、59.60 mg／L。矿渣堆被冲刷过程中，各支流污染物质量浓度总体变化趋势大致相同，在入流点到矿渣堆上游污染物质量浓度呈下降趋势，在矿渣堆上游至矿渣堆下游污染物质量浓度呈上升趋势。其中竹园沟上游矿渣堆积范围较广，水流到达该支流与主河道交汇处时，污染物质量浓度达到最大，之后受上游水流汇入、污染物衰减过程及河网本底特征影响，污染物质量浓度呈减小趋势。清石沟矿渣堆面积较小，在水流经矿渣堆到达主河道之前，两种污染物质量浓度经稀释和底泥吸附等作用有所下降。支流孙家沟地表水的硫酸根离子和铁离子质量浓度均小于上游来水的，水流汇合后对主河道污染物具有一定稀释作用，使主河道污染物质量浓度下降；宝泉寺和郭家沟矿渣堆直达下游与主河道交汇处，故在交汇前硫酸根离子和铁离子质量浓度均较大，且宝泉寺和郭家沟支流流量较大，因此水流进入主河道后，硫酸根离子和铁离子质量浓度均发生了明显变化，致使主河道下游污染物质量浓度增大。矿渣堆下游是河网污染物集中整治的关键区域。

图2 硫酸根离子质量浓度空间分布

图3 铁离子质量浓度空间分布

3.2 水量及污染物通量分析

各支流对小文峪河干流水量贡献不一，由大到小依次为宝泉寺、郭家沟、竹园沟、清石沟和孙家沟，占比分别为21.97%、21.21%、20.46%、18.94%、17.42%。整体水量差距相对不大，但由于各支流下游污染物质量浓度不同，因此各支流硫酸根离子和铁离子负荷对主河道的贡献率差距较大。各支流水流和日污染物通量见图4。各支流输入主河道的污染物通量由大到小依次为宝泉寺、郭家沟、竹园沟、清石沟、孙家沟，其中硫酸根离子通量分别为1.135、0.989、0.518、0.157、0.079 t／d，铁离子通量分别为0.204、0.139、0.026、0.019、0.007 t／d。硫酸根离子和铁离子入河通量在不同支流变化趋势较为一致，最大值均出现在宝泉寺。受水量和矿渣堆积方量不同影响，污染物通量在各支流呈现不同变化，以硫酸根离子为例，其污染物入河通量的最大值出现在宝泉寺，最小值出现在孙家沟，最大值为最小值的14.4倍，而宝泉寺水量仅为孙家沟的1.26倍，表明污染物通量的差异主要与矿渣堆积方量有关。

图4 污染物通量及水量

3.3 空间分布特征

各支流硫酸根离子、铁离子入河通量大小顺序均为宝泉寺＞郭家沟＞竹园沟＞清石沟＞孙家沟，可见入河污染物主要来源于宝泉寺和郭家沟，二者占比合计在70%以上，其他支流的占比相对较小。硫酸根离子污染负荷主要来自郭家沟和宝泉寺，其对小文峪河主河道的贡献率分别为34.36%、39.42%；除郭家沟和宝泉寺外，竹园沟对河网中硫酸根离子的贡献率也较大，为17.99%。在各支流中，宝泉寺对小文峪河主河道铁离子的贡献率最大，为51.64%，其次为郭家沟，为35.10%。由各支流对小文峪河主河道污染物的贡献率可知，整体上，宝泉寺的污染负荷大于郭家沟的，其次为竹园沟、清石沟和孙家沟；流域内污染负荷主要来自下游支流宝泉寺和郭家沟，其产生的硫酸根离子、铁离子污染负荷分别占整个研究区的73.78% 和86.74%。

各支流污染物通量及贡献率不同的主要原因是，矿渣堆积方量不同，其中宝泉寺和郭家沟矿渣堆占地面积分别为16 638、20 824 m2，其他3 条支流矿渣堆占地面积之和为15 350 m2；宝泉寺和郭家沟流量较大，水流对矿渣堆的冲刷较为严重，导致污染物质量浓度增大，污染负荷增大。因此对小文峪河流域的水环境治理方案制定过程中，应首先考虑宝泉寺和郭家沟两条支流的整治。

4 结论

以十堰市小文峪河流域为例，基于二维水动力水质耦合模型，充分利用各支流矿渣堆属性、流域河网本底特征及地表水、矿洞水、淋溶水水质和水量监测数据等，构建了基于GPU 加速技术的矿山区河网水环境高分辨率数值模型，对小文峪河流域硫酸根离子和铁离子的空间分布进行了可视化模拟，并量化分析了小文峪河污染负荷来源。模型模拟结果与污染物实测值相比，相对误差小于10%，表明模型可以对矿山区河网水环境变化过程进行模拟。小文峪河主河道的水量贡献率由大到小分别为宝泉寺、郭家沟、竹园沟、清石沟和孙家沟，占比分别为21.97%、21.21%、20.46%、18.94%、17.42%。硫酸根离子和铁离子对小文峪河贡献率由大到小均为宝泉寺＞郭家沟＞竹园沟＞清石沟＞孙家沟，其中宝泉寺和郭家沟硫酸根离子贡献率分别为39.42%和34.36%，铁离子贡献率分别为51.64%和35.10%。

综上，基于GPU 加速技术的水动力水质耦合模型对矿山区河网水环境进行模拟效果较好，在小文峪河流域水环境治理方案制定过程中，应首先考虑宝泉寺和郭家沟两条支流的整治。