近年来，我国的突发水污染事故进入了密集高发时段^[1-2]。长江作为黄金水道，水上运输发达，且沿江多工业企业，突发污染事故风险高。据统计，自2000年以来长江流域较大的水污染事件就发生了50多起^[3]。不同于常规污染，突发污染事故具有不确定性与随机性^[4]，污染物在瞬间或短期内大量进入水体，在破坏水域生态环境的同时也对水源地的供水安全造成了严重威胁，影响社会经济稳定。

为了降低突发水污染事故的危害，对事故的模拟研究也逐渐成为热点，目前主要以水动力、水质数学模型为基础^[5-9]，结合情景设定^[10]来研究污染物的迁移扩散过程，确定污染物的运动规律和时空分布特征，进一步为突发污染事故预警与应急提供依据^[11]。不少学者对长江突发污染事故进行了模拟研究，例如袁玥等^[12]应用MIKE21FM二维模型对长江蕲春段排污事故进行模拟，量化污染事故的影响范围与强度；张楚天等^[13]采用平面二维水质模型与GIS集成，实现对长江武汉段突发污染事故的动态可视化模拟。不同于长江中上游河段，长江下游河段多汊道地形^[14]且受到潮汐作用的影响^[15]，水动力条件复杂，目前对长江下游感潮河段这一复杂水文条件下的突发污染事故模拟较少。且由于地形资料、水文资料等的限制，突发污染事故的水动力基础多选取事故发生的特定河段^[16]，模拟范围较小，构建大范围尺度的水动力模块作为基础来模拟突发污染事故的研究鲜有报道。

笔者构建了大通到徐六泾段的二维水动力模型，在水动力模块的基础上，选取事故风险高且水文条件复杂的南京段进行突发污染事故模拟。基于长江南京段潜在的事故风险类型，同时设定工业和船舶2种突发污染事故情景，模拟污染团的迁移扩散过程以及到达各个水厂取水口的时间及浓度，研究结果可对长江水源地取水用水进行安全预警，同时为水源地应急调水提供依据。

1 数据来源与研究方法 1.1 数据来源

研究区的地形文件由南京市航道局2009年编制的长江下游航行参考图数字化生成，总长度约531 km，共用航道图20张，比例尺皆为1:40 000。流量数据采用2014年大通水文站实测日流量数据。潮位数据采用2014年南京站、镇江(二)站、江阴站、徐六泾站实测日潮位数据，每日有4个潮位数据，包括2个高潮位和2个低潮位。

1.2 研究方法 1.2.1 模型选取

长江下游干流的感潮河段范围广，且江面宽度远大于深度，因此采用二维水动力模型作为突发污染事故的研究基础。选取MIKE21 Flow Model FM模型进行模拟，该模型采用的是FM网格结构，可对网格信息进行编辑，在研究区适用性更高。

二维水动力模块建立于Navier-Stokes方程上，将三维流动控制方程简化为二维浅水流动质量和动量守恒控制方程(公式1~3)。水质模块建立于二维空间下的对流扩散方程(公式4)。

连续性方程:

$ \frac{{\partial h}}{{\partial t}} + \frac{{\partial h\bar u}}{{\partial x}} + \frac{{\partial h\bar v}}{{\partial y}} = hS。$

(1)

x方向动量方程:

$ \begin{array}{l} \frac{{\partial h\bar u}}{{\partial t}} + \frac{{\partial h{{\bar u}^2}}}{{\partial x}} + \frac{{\partial h\bar u\bar v}}{{\partial y}} = f\frac{{g{h^2}\partial \rho }}{{2{\rho _0}\partial x}} + \frac{{{\tau _{sx}}}}{{{\rho _0}}} - \\ \frac{{{\tau _{bx}}}}{{{\rho _0}}}\left( {\frac{{\partial {S_{yx}}}}{{{\partial _x}}} + \frac{{\partial {S_{yy}}}}{{{\partial _y}}}} \right) + \frac{\partial }{{{\partial _x}}}\left( {h{T_{xy}}} \right) + \frac{\partial }{{{\partial _y}}}\left( {h{T_{yy}}} \right) + h{v_s}S。\end{array} $

(2)

y方向动量方程:

$ \begin{array}{l} \frac{{\partial h\bar v}}{{\partial t}} + \frac{{\partial h\overline {uv} }}{{\partial x}} + \frac{{\partial h{{\bar v}^2}}}{{\partial y}} = f\frac{{\partial {P_{\rm{a}}}}}{{\partial y}} - \frac{{g{h^2}\partial \rho }}{{2{\rho _0}\partial y}} + \frac{{{\tau _{sy}}}}{{{\rho _0}}} - \frac{{{\tau _{by}}}}{{{\rho _0}}} - \\ \frac{1}{{{\rho _0}}}\left( {\frac{{\partial {S_{yx}}}}{{{\partial _x}}} + \frac{{\partial {S_{yy}}}}{{{\partial _y}}}} \right) + \frac{\partial }{{{\partial _x}}}\left( {h{T_{xy}}} \right) + \frac{\partial }{{{\partial _y}}}\left( {h{T_{yy}}} \right) + h{v_s}S。\end{array} $

(3)

对流扩散方程:

$ \frac{{\partial C}}{{\partial t}} + u\frac{{\partial C}}{{\partial y}} + v\frac{{\partial C}}{{\partial y}} = {D_x}\frac{{{\partial ^2}C}}{{\partial {x^2}}} + {D_y}\frac{{{\partial ^2}C}}{{\partial {y^2}}}。$

(4)

式(1)~(4)中, t为时间，s；(x，y)为直角或斜角坐标系下的坐标，m；(u，v)为(x，y)方向水流流速，m·s^-1；h为总水深，m；η为潮位，m；S为源汇项，kg·m^-3·s^-1；(u_s，v_s)为源汇项在(x，y)方向水流速度，m·s^-1；g为重力加速度，m·s^-2；ρ和ρ₀分别为水体和参考水密度，kg·m^-3；P_a为水表面大气压，P_a；τ_sx、τ_bx、τ_sy、τ_by为底部剪切应力，N；S_xx、S_xy、S_yx、S_yy为辐射应力分量，N；f为科氏力系数，m^-1；T_xx、T_xy、T_yx、T_yy为粘滞切应力分量，N；C为物质质量浓度，kg·m^-3；(D_x，D_y)为(x，y)方向上的扩散系数，m²·s^-1。字母上方有一横杠的表示平均值。

1.2.2 水动力模块构建与验证

采用非结构三角网格对河道地形进行处理，三角网格总数为79 685个，最小网格分辨率为100 m。如图 1所示，上边界设为大通水文站，下边界设为徐六泾潮位站，模型计算的时间步长设置为10 min。

长江中下游的糙率系数(n)一般取值在0.014~0.030之间^[17]。该研究的二维水动力模型中采用曼宁系数(M)，M与n之间存在倒数关系，即M=1/n。由于选取的河道较长，且上下游河道的地形差距较大，因此将河底地形考虑在内，使用水深数据对公式进行优化。改进后的曼宁系数公式为

$ M = \frac{1}{{{n_0} + {n_1}/|d|}}。$

(5)

式(5)中，M为曼宁系数；n₀和n₁为底部的粗糙程度，n₀的取值范围一般为0.01~0.02，n₁的取值范围一般为0.005~0.01；d为水深，m。在取值范围内选取不同的n₀和n₁值对水位进行模拟，经过率定得出n₀=0.014，n₁=0.01。

选取2014年2月南京站、镇江(二)站、江阴站潮位数据进行验证，将这3个站点的实测潮位数据与模型的模拟水位数据进行对比，结果如图 2所示。模拟值与实测值吻合性良好，可作为水质模块的水动力基础。

1.2.3 突发污染事故情景设计与参数设置

1989年扬子石化公司人员误将设备剩余物料倒入马汊河，造成附近水厂供水无法饮用；2014年，扬子石化炼油厂酸罐爆炸，导致5个储罐泄漏着火，污染物通过排水口进入马汊河；2013年“鑫川8号”船舶碰擦南京长江大桥后，行至八卦洲洲头沉没，导致沉船内80多t燃料重油溢出。综合来看，长江南京段水源安全威胁主要来自于2个方面:一是由于企业或其他污染源超常规排放；二是由于车辆、船舶运输事故造成的污染物泄漏。突发污染物类型集中在油类、化学品(硫酸、苯酚和苯类等)和重金属(以镉为代表)上。结合长江南京段水源地的分布，设定2种突发污染事故情景(图 3):一为马汊河工业事故，二为南京长江三桥船舶事故，同时选取对长江南京段水源威胁大的苯酚作为特征污染物。

考虑水动力条件最弱情况下污染物扩散对水源地产生的影响最大^[18-19]，选取枯水期(2014年2月)作为水动力基础，假定100 t苯酚在2 h内进入水体。考虑事故发生点到水源地的距离(表 1)，在马汊河工业事故中设定落潮和涨潮2个事故发生的时间点，分别为2月5日06:00与2月5日12:00；在三桥船舶事故中，设定事故发生时间点为落潮段2月5日06:00。

在泄露污染初期污染源的扩散情况尚不清晰，且长江水量丰沛，故仅考虑苯酚的横向扩散，以污染最大化来建立预测^[20-21]。另由于苯酚属于非常规污染物类型，无常规水质监测数据资料，故结合其他参考文献设定苯酚的白天降解系数为0.20 d^-1[22]，夜晚降解系数为0.176 d^-1[23]，扩散系数为10 m²·s^-1[24]。

2 结果与分析

根据GB 3838—2002《地表水环境质量标准》，在水质检测中，Ⅲ类水允许的苯酚极限值为0.005 mg·L^-1。因此在输出浓度结果时，将浓度值超出0.005 mg·L^-1的区域定义为受影响区域。

2.1 情景一:马汊河工业事故

如图 4所示，在落潮和涨潮的水动力条件下上坝水源地均能检测到污染物，但污染物浓度量级仅为10^-11 mg·L^-1，因此可认为在潮汐作用中污染团会随着潮汐上溯，但对上坝水源地基本无影响；龙潭水源地在落潮和涨潮下均受到影响，表现为落潮情况下污染团先到达龙潭水源地，浓度峰值低于涨潮，2种情况下污染物的整体运动特征一致。

以落潮情况为例进行分析(表 2):龙潭水源地在事故发生后24~54.5 h内受到影响，总时长约为30.5 h。在此过程中，污染物浓度峰值发生了5次连续波动。考虑到水流的往复运动，可概化成2个连续的峰值过程:事故发生后30.6~35.3 h污染物质量浓度达到首个峰值，为0.196 mg·L^-1，波峰持续时间为4.7 h；事故发生后41.5~44.3 h污染物质量浓度到达第2个较小的峰值，为0.094 mg·L^-1，波峰持续时间为2.8 h。

第1次峰值过程污染团受径流作用，从上游向下游运动〔图 5(a)〕；第2个峰值过程污染团受到涨潮流控制，从下游向上游回溯^[25] 〔图 5(b)〕。第2个峰值浓度比第1个小，且对龙潭水源地的影响时间更短，这是由于在污染团的来回运动中污染物被稀释和降解，使得回溯过程对水源地的影响较小。

由此可见，由于潮流的控制，污染团会发生往复运动，水源地多次出现污染物浓度峰值，受影响的时间也随之加长。在马汊河工业事故情景下，污染团仅对龙潭水源地产生影响，并未上溯到马汊河上游的上坝水源地和燕子矶水源地。因此在龙潭水源地受到影响的30.5 h内可就近使用上坝和燕子矶水源地作为龙潭水源地的应急水源。

2.2 情景二:南京长江三桥处船舶事故 2.2.1 江浦—浦口、燕子矶、上坝、龙潭水源地

如表 3所示，情景2南京长江三桥事故对下游5个水源地均有影响，除夹江水源地外，三桥北侧、南侧事故对其他4个水源地的影响基本相同。以北侧事故为例，江浦—浦口水源地最先受到影响，持续时间为18.5 h；燕子矶为第2受影响水源地，持续时间为24.5 h；第3受影响的是上坝水源地，持续时间为27.8 h；最后受影响的是龙潭水源地，持续时间为145.8 h。

各水厂在受影响的过程中污染物浓度均出现多次峰值。以北侧事故为例(图 6)，其中北河口、城南、江浦、浦口、上元门水厂出现2次峰值，城北、远古水厂出现3次峰值，龙潭水厂出现6次以上峰值。综合水源地来看，江浦—浦口水源地在事故发生后11.5 h出现浓度峰值，为0.667 mg·L^-1；燕子矶水源地在事故发生后16.8 h出现浓度峰值，为0.440 mg·L^-1；远古水源地在事故发生后29.3 h出现浓度峰值，为0.283 mg·L^-1；龙潭水源地在事故发生后47.8 h出现浓度峰值，为0.193 mg·L^-1。

同情景1，江浦—浦口、燕子矶、上坝水源地水厂取水口出现2~3次峰值，为径流、潮流共同作用的结果，峰值的出现是连续的。但值得注意的是，龙潭水源地峰值的出现不是完全连续的，这是由于污染团在经过八卦洲汊道时分为2团分别沿着长江干流、八卦洲左汊进行运动(图 7)。长江干流携带的污染团先到达龙潭水源地，在径流、潮流的共同作用下来回运动，使得龙潭水源地出现前5个峰值，峰值质量浓度为0.045~0.189 mg·L^-1；八卦洲左汊的污染团后到达龙潭水源地，使龙潭水源地出现后4个峰值，质量浓度低于0.019 mg·L^-1。

综上可知，水源地距离事故点的位置越远，受事故的影响越晚，但受影响的时间会加长；水源地出现的峰值个数随距离增加而增加，但峰值浓度值减小。这是由于在污染团的运动过程中一方面污染带长度随着迁移扩散距离的增加而增加，在径流、潮流共同作用下往复过程多，污染物浓度不断下降；另一方面，长江南京段汊道较多，污染团通过汊道时分散，导致位于汊道下游的水源地经过多个污染团，下游水源地受影响的时间加长、污染物浓度峰值个数增加。

2.2.2 夹江水源地

如图 8所示，夹江水源地在2次事故中受到的影响差距很大。在北侧事故中，夹江水源地受影响的时间为2月5日18:20—2月6日15:00，持续时间为20.7 h；位于水源地内的北河口水厂出现2个质量浓度峰值，最大值为0.032 mg·L^-1，城南水厂也出现2个峰值，最大值为0.006 mg·L^-1。在南侧事故中受影响的时间为2月5日14:50—2月12日07:40，持续时间为160.8 h；北河口水厂出现6个以上质量浓度峰值，最大值为0.055 mg·L^-1，城南水厂出现6个峰值，最大值为0.199 mg·L^-1。

分析二维水质模型，在三桥北侧事故中污染团少量进入夹江水源地，仅对上游城南水厂产生影响，未到达下游北河口水厂，北河口水厂的污染团来自于长江主江段潮流回溯过程；在三桥南侧事故中污染团进入夹江水源地，对城南水厂、北河口水厂均产生影响，并在夹江水源地内来回运动。

夹江水源地地形狭窄，水流量小，一旦污染团进入，会在水源地内来回运动，持续影响水源地水质。因此，对于夹江水源地而言，要着重防止污染团进入夹江段，同时要防止由于潮流作用污染团上溯进入北河口水厂。

3 结论

基于构建的大通到徐六泾二维水动力模型，结合长江南京段实际污染来源，设定马汊河工业污染和长江三桥船舶污染2种突发污染事故情景，模拟污染物苯酚的迁移扩散过程。

(1) 该研究构建的大通到徐六泾二维水动力模型中，南京站、镇江站和江阴站的模拟潮位与实测数据吻合度良好，为感潮河段下突发污染事故污染物迁移扩散的模拟奠定了可靠的基础。

(2) 在设定的2种典型突发污染事故情景中，马汊河工业事故仅对下游龙潭水源地有30.5 h的影响；三桥船舶事故中，江浦—浦口、燕子矶、上坝、龙潭水源地受影响的时间分别为18、24、28和146 h，夹江水源地在三桥北侧事故中受影响时间为18 h，在三桥南侧事故中为154 h。

(3) 长江南京段受到径流和潮汐作用的共同影响，污染物在河道内往复运动，随着污染物的迁移扩散，污染带长度增加，下游水源地受影响时间加长。长江南京段具有独特的汊道地形，污染团进入不同汊道，对下游龙潭水源地造成多次污染。

该研究构建的大通到徐六泾的二维水动力模型为长江下游突发污染事故的模拟提供了一定的研究基础。同时模拟过程中发现潮汐作用对污染物的迁移扩散产生一定影响，但均未上溯到事故发生点之上的水源地，后续可进一步进行不同潮汐类型以及涨潮落潮情况下的研究，为突发污染事故的应急、水源地选址等提供依据。