长江干流铜陵段突发水污染情景模拟分析

齐于顺，刘仁志*，张启月，崔保山，郭忠，王飞，方丽，方陵

1.北京师范大学环境学院

2.铜陵市生态环境局

突发环境事件是指由于污染物排放或自然灾害、生产安全事故等因素，导致有毒有害物质大量进入环境介质，突然造成或可能造成环境质量下降，危及公众身体健康和财产安全，需要采取紧急措施予以应对的事件[1]。与其他普通污染事件不同的是，突发环境事件具有随机性、复杂性、高强度、高风险的特点[2]，会在短时间内造成周边环境污染物超标，甚至会造成长期影响[3]。其中，突发水环境事件最易发生，李旭等[4]对2011—2017年突发环境污染事件的受体类型进行统计发现，水是主要污染受体，占比高达85.59%。长江沿岸城市“化工围江”、上下游排污口与取水口交错分布等布局性风险突出，导致长江水环境风险形势不容乐观[5-6]。此外，由于长江水体流动性强的特点，突发水环境事件的影响范围更大，响应时间更短，如果不能及时妥善处置，短时间内可能导致饮用水水源地污染、水域生态环境破坏，甚至造成跨界污染，产生极大的不良社会影响[7-8]。目前突发水环境事件仍是影响供水和人民安全的重要威胁因素[9-10]。因此，有必要通过多种方式提前做好突发水环境事件的预测预警与应急准备，有效防范突发水污染事故，降低突发水环境风险[11]。

水环境模型可对水环境系统及其内部发生的复杂过程进行定量化描述[12]，模拟污染物在一定时间和空间范围内的迁移转化过程，计算污染物到达下游环境敏感目标的时间、污染团的历经时间和污染物的浓度分布等[13]，有助于帮助应急管理人员了解事故可能的发展趋势并做出合理的处置。常用的河流水质模型主要有EFDC、WASP、QUAL2E、QUAL2K、MIKE、BASINS、SWAT 等[13-17]。其中，丹麦水力研究所开发的平面二维模型MIKE21，因用户界面友好，具有多种计算网格、模块供用户选择，且有强大的前、后处理功能，被广泛用于河流[18-19]、湖泊[20-21]、水库[22]、海湾[23]等水质模拟。如袁玥等[18]应用MIKE21模型对长江蕲春段非正常排污情况进行了模拟，舒长莉等[24]应用MIKE21模型对赣江南昌段突发污染事故进行了模拟，均实现了对污染影响后果的量化表达。为贯彻落实习近平总书记“不搞大开发，共抓大保护”的重要指示精神，着力解决长江水环境、水生态、水风险突出问题[25]，为长江经济带高质量发展提供保证[26]，笔者以长江生态环境保护修复驻点跟踪研究工作为基础，利用MIKE21构建长江干流铜陵段突发水污染事故水动力-水质耦合模型，考虑了多种可能的泄漏情景，并对不同水文期进行模拟分析，实现对突发水污染事故的动态模拟，以期为长江干流铜陵段水域突发水污染事故风险防范、预测预警和应急处置提供参考。

1 研究区概况

铜陵市位于安徽省中南部、长江下游（117°04’E～118°09’E，30°38’N～31°09’N），是长江经济带重要节点城市。长江干流铜陵段全长142.6 km，过境的长江水资源丰富，是铜陵市主要供水来源，最小日平均流量为10 700 m3/s，最大日平均流量为68 400 m3/s，年平均流量为29 250 m3/s（2019年），水量年内变化显著。研究河段内环境敏感目标包括铜陵市淡水豚国家级自然保护区和3处集中式生活饮用水水源地（铜陵市水厂、新三水厂和第五水厂）。淡水豚国家级自然保护区范围上始铜陵市枞阳县老洲镇，下至铜陵市义安区金牛渡，全长58 km，其中核心区、缓冲区和实验区面积分别为9 534、6 360、15 624 hm2[27]。横港扫把沟地区紧临长江，是铜陵市最早的工业聚集区之一，由于历史原因，该区域功能布局不合理，工业用地与居住用地混杂，多家企业距长江直线距离不足1 km。目前该区域仍存留多家风险企业，也没有成立工业园区管委会，这给下游集中式饮用水水源地、淡水豚国家级自然保护区以及长江水体水质安全带来重大威胁。研究区概况如图1所示。

图1 研究区概况Fig.1 Schematic map of the study area

2 数据来源与研究方法

2.1 数据来源

研究区水文数据采用2019年大通水文站实测逐日流量数据和横港、荻港2个水位站实测逐日水位数据；地形文件为1∶25 000实测水下地形数据；风险企业、饮用水水源保护区和淡水豚保护区等相关资料由铜陵市生态环境局提供。

2.2 MIKE21模型基本原理

长江铜陵段江面宽度远大于深度，且不存在明显分层现象，因此采用二维水动力模型可以满足研究需要[28]。选取MIKE21 Flow Model FM模型进行模拟，二维水动力模型基于三向不可压缩和Reynolds值均布的Navier-Stokes方程，并服从于Boussinesq假定和静水压力的假定，水动力控制方程组（浅水方程）为[24,29-31]：

式中：t为时间；x、y为笛卡尔坐标； η为水位，m；d为静止水深，m；h为 总水深，h=d+η，m；u、v分别为x、y方向上的速度分量；f为哥氏力系数，f=2ωsinφ，其中 ω 为地球自转角速度，φ为当地纬度；g为重力加速度，取9.81 m /s2；ρ和 ρ0分别为水体和参考水的密度， k g/m3；Pa为大气压力，Pa；sxx、sxy、syx、syy分 别为x轴的法线方向、y轴的法线方向、xy表面上的切线方向和yx表面上的切线方向的辐射应力分量；Ax、Ay为应力分项；Txx、Txy、Tyx、Tyy为黏滞切应力分量；S为源汇项；us、vs为源项水流流速，m/s；u¯、v¯为沿水深平均的流速；τsx、τsy为水面风应力在x、y方向上的分量;τbx、τby为河床 ( 海底) 底部应力的分量。

u¯、v¯计算公式如下：

水质模块建立于二维空间下的输移扩散方程为：

式中：C为各典型污染物浓度；Dx、Dy为各典型污染物在x、y上的扩散系数；P为各典型污染物降解项。

MIKE21模型采用有限体积法对计算区域进行空间离散，将该连续统一体细分为若干个不重叠的三角形或四边形单元。模型计算的时间和精度取决于计算数值方法所使用的求解格式精度，浅水方程的时间积分和输移扩散方程基于半隐格式求解，相应平流项采用显式格式求解，垂直对流项采用全隐格式求解。受显式格式稳定性的限制，为保持模型计算的稳定性，模型中时间步长的设定必须保证CFL（Courant-Friedrich Levy)数小于1。浅水方程和输移扩散方程在笛卡尔坐标上的CFL分别定义为：

式中： Δx和 Δy为x、y方向的特征长度； Δt为时间间距。 Δx和 Δy近似于三角形网格的最小边长，水深和流速值为三角形网格中心的取值。

2.3 模型建立与验证

利用MIKE21建立水动力-水质耦合模型，模型构建流程如图2所示。基于铜陵市水利工程图和实测水下地形数据，利用ArcMap与AutoCAD软件提取生成长江铜陵段水陆边界线（Land.xyz文件），并读取相应水深散点数据（Water.xyz文件）。在Mesh Generator 中导入水陆边界线（Land.xyz文件）并进行光滑处理，采用非结构三角网格对河道地形进行处理，然后生成网格；插入地形水深散点数据（Water.xyz文件）并进行地形插值后导出mesh文件（图3），计算区域共有10 896个三角网格，6 886个网格节点。模型上边界采用2019年大通水文站逐日流量时间序列，下边界采用2019年荻港水位站逐日水位时间序列。选择Flow Model FM建立水动力-水质耦合模型，经过多次调试后选择合适参数；将模拟时长设为24～30 h，时间步长选择60 s；为保证模型稳定运行，设置CFL为0.8，最大时间步长为60 s；干水深(hdry)为0.005m，淹没水深（hflood）为0.05m，湿水深(hwet)为0.1m；底床摩擦力采用曼宁系数，取值32 m1/3/s；初始水位根据不同水文期选择模拟对应开始日期水位值；其他参数采用推荐或默认值。在对流扩散模块，添加需要模拟的风险物质硫酸，从更安全的角度考虑预测结果，将硫酸降解系数设为0[32]；并在源汇项中添加风险物质初始浓度和泄漏流量等相关数据。

图2 模型构建流程Fig.2 Model construction flowchart

图3 水动力模型网格地形高程Fig.3 Grid terrain elevation of hydrodynamic model

因长江铜陵段水文情势变化显著，有明显的丰、平、枯水文情势节律性变化，且枯水期部分时间段水位站处于不工作状态。故对丰、平、枯水文时期部分时间段进行模拟，并选取2019年横港水位站逐日水位数据进行验证，将实测水位数据与模拟水位数据进行对比。经对照发现（图4），模拟水位与实际水位数值相差不大，平均相对误差为2.4%，曲线趋势基本相似，波峰与波峰相对，波谷与波谷相对。模拟值与实测值吻合性良好，可作为水质模块的水动力基础。

图4 横港水位站2019年模拟水位与实测水位Fig.4 Simulated and measured water level of Henggang water-level station in 2019

2.4 泄漏情景设定与水质模拟

横港扫把沟地区多家风险企业的储罐内存在大量风险物质，给下游集中式饮用水水源地、淡水豚国家级自然保护区以及长江水体水质带来巨大安全隐患。针对某化工企业硫酸储罐，设定3种泄漏情景，分别为罐体20%管径破裂、罐体100%管径破裂和储罐完全泄漏，泄漏时间均设为10 min。在丰水期、平水期和枯水期对3种泄漏情景分别进行模拟，具体泄漏时间分别设定为7月10日、4月16日和10月12日的06:00。所设定的泄漏情景如表1所示。

表1 3种设定泄漏情景Table 1 Three supposed leakage situations

硫酸罐体20%管径和100%管径破裂时泄漏速率用伯努利方程[33]计算（液体在泄漏口没有急骤蒸发），公式如下：

式中：QL为液体泄漏速率， k g/s ；P1为容器内介质压力，Pa；P0为环境压力，Pa；ρ为泄漏液体密度， k g/m3；h为 裂口之上液位高度， m ；Cd为液体泄漏系数，取0.62；A为裂口面积， m2。

根据实地调研，判断突发泄漏情况下的可能入江位置如图1所示，将不同泄漏情景下计算得到的泄漏流量-时间序列等相关数据输入构建的水动力-水质耦合模型，可得到不同水文期的风险物质迁移扩散分布结果。

3 结果与讨论

采用MIKE21水动力-水质耦合模型对3种泄漏情景在丰水期、平水期和枯水期的风险物质迁移扩散情况分别进行了模拟，定量化表达出风险物质到达下游环境敏感目标的时间、最大污染峰团浓度和污染持续时间等。泄漏硫酸对长江水质的影响可通过对水体pH的影响来反映[32]，GB 3838—2002《地表水环境质量标准》规定地表水Ⅰ类～Ⅴ类水域水体pH为6～9，当pH为6时，通过计算得到对应硫酸的浓度为0.049 mg/L，即硫酸浓度大于0.049 mg/L时会对下游各风险受体和长江水质造成影响。通过水质模型模拟结果可得到储罐泄漏后各时间节点的浓度，进而预测和判断影响程度。下面针对同一水文期不同泄漏情景和同一泄漏情景在不同水文期的水质模拟结果进行说明。

3.1 不同泄漏情景下水质模拟结果

以3种不同泄漏情景在丰水期模拟结果为例进行说明，模拟结果如图5～图7所示。3种泄漏情景下最大污染峰团到达下游敏感受体所需时间相同，到达三水厂取水口、市水厂取水口、淡水豚核心区、五水厂取水口和市出境断面的时间分别为89、93、332、677、751 min，这是由于同一时间段的水动力结果相同。但在储罐完全泄漏情景下，风险物质到达下游敏感受体时间稍早于其他2种泄漏情景，风险物质完全流过敏感受体时间稍晚于其他2种泄漏情景，即影响时间更长。此外，该情景下最大污染峰团浓度也明显高于其他2种泄漏情景下2～3个数量级。在三水厂取水口处，20%管径破裂时最大污染峰团浓度为1.891 mg/L，100%管径破裂时最大污染峰团浓度为47.273 mg/L，储罐完全泄漏时最大污染峰团浓度为1 381.39 mg/L，这是由于不同泄漏情景下风险物质泄漏总量不同。随着时间的延长，风险物质被稀释，到达下游各敏感受体时的最大污染峰团浓度会逐渐降低。在3种不同泄漏情景下，风险物质到达市出境断面时的最大污染峰团浓度分别为0.090 5、2.262 4和66.123 4 mg/L，均超过0.049 mg/L的标准要求，会造成不同程度的跨界污染。淡水豚核心保护区范围最大，其受影响持续时间也最长，分别为 251、313、343 min。

图5 丰水期罐体20%管径破裂时水质模拟结果Fig.5 Simulation results of water quality when 20% pipe diameter of tank breaks in wet period

图6 丰水期罐体100%管径破裂时水质模拟结果Fig.6 Simulation results of water quality when 100% pipe diameter of tank breaks in wet period

3.2 不同水文期泄漏水质模拟结果

以不同水文期储罐完全泄漏情景的模拟结果为例进行说明，模拟结果如图7～图9所示。根据模拟结果可知，丰水期风险物质到达下游敏感受体时间最短，其中风险物质到达三水厂取水口仅用75 min，影响持续时间为43 min。相较于丰水期，平水期和枯水期的流速更小，风险物质需更长时间到达下游敏感受体，如平水期、枯水期风险物质到达三水厂取水口的时间分别为103和111 min。但在平水期和枯水期，风险物质对敏感受体的影响持续时间会更长，其对三水厂取水口的影响持续时间分别为65和72 min。同时因流量更小，风险物质在市水厂取水口、三水厂取水口和淡水豚核心区的最大污染峰团浓度会更高，枯水期的浓度分别为1 883、1 736和426 mg/L。而到达下游更远的五水厂取水口时，丰水期、平水期、枯水期浓度依次降低，其主要原因可能是停留时间较长，风险物质被稀释[24]。平水期和枯水期，在五水厂取水口、市出境断面处会先后出现2次污染峰团，这是由于江心洲两侧水体流速相差较大，部分风险物质随右侧支流先行到达敏感受体，随后左侧支流中的风险物质再次到达敏感受体。

图7 丰水期储罐完全泄露时水质模拟结果Fig.7 Water quality simulation results when the tank breaks completely in wet period

图9 枯水期储罐完全泄露时水质模拟结果Fig.9 Water quality simulation results when the tank breaks completely in dry period

图8 平水期储罐完全泄露时水质模拟结果Fig.8 Water quality simulation results when the tank breaks completely in normal period

综上，风险物质对下游风险受体的影响程度与泄漏总量和水文期密切相关，所得结论与长江其他段类似研究[29,31]相同。同一水文期，风险物质泄漏总量越大，对下游敏感受体的影响时间越长，污染团浓度也越高。泄漏风险物质在平水期和枯水期时，对受体的影响时间更长；在丰水期时，到达受体的时间更短，即应急反应时间更少。

4 结论

（1）利用MIKE21水动力-水质耦合模型，对长江铜陵段某化工企业储罐内风险物质在不同水文期、3种不同泄漏情景下的突发水污染事故进行了模拟，所建立的水动力模型模拟水位结果与实测结果吻合性良好，曲线趋势基本相似，表明模型选取参数较为合理，可作为水质模块的水动力基础。

（2）在设定的3种不同泄漏情景下，储罐完全泄漏时因其泄漏总量大，导致风险物质浓度最高，影响时间最长。在不同水文期的模拟结果显示，丰水期风险物质到达下游敏感受体的时间最短，留给应急人员处置时间也最少；平水期和枯水期风险物质的最大污染峰团浓度会更高，影响时间会更长，且在五水厂取水口和市出境断面处先后有2次污染峰团到达。

（3）本研究基于MIKE21构建了适用于长江干流铜陵段的水动力-水质模型，通过对风险企业内多种风险物质在不同水文情势和泄漏情景下的模拟分析，可建立当地的模拟预测预警资料库，为突发水污染事故的预测预警、应急处置提供决策参考。结合本次模拟结果，因风险物质在丰水期到达受体时间更短，平水期和枯水期对受体影响时间更长，当接到突发水环境事件报告时，在丰水期应立即采取措施保护风险受体，在平水期和枯水期则首先对泄漏风险源进行控制，从而减少影响时长。此外，有必要通过加强日常环境安全检查等措施，及时发现并整改问题，从而有效降低对水环境的影响程度。后续研究还需采集更为丰富的水文数据、水质监测数据等，不断进行完善以进一步提高模型模拟的准确性。