鸭绿江公路大桥溢油漂移扩散三维数值模拟

2012-10-31陈士谦朱建荣

华东师范大学学报（自然科学版） 2012年6期

陈士谦，朱建荣

（华东师范大学河口海岸学国家重点实验室，上海 200062）

0 引言

为了增进中朝友谊，促进中朝交流，鸭绿江浪头港附近将建造连接中国和朝鲜的公路大桥（见图1）.鸭绿江公路大桥建成后，因桥墩的存在，将增大船舶航运的危险.若发生原油泄漏事故，会对周围水域生态环境产生严重的后果，特别是会污染上游丹东市取水口（四道沟取水口、丹东造纸厂取水口和燕窝取水口），影响居民生活及工业用水的安全.因此，研究鸭绿江公路大桥在发生溢油事故后油膜漂移扩散，可为溢油事故的处理提供油膜面积、到达取水口时间和油膜厚度变化等数据，对鸭绿江环境的保护和取水的安全具有实际意义.

图1 鸭绿江河口示意图Fig.1 Sketch map of the Yalu Estuary

目前预测溢油事故发生后油膜漂移扩散的主要方法是数值模拟，一般以垂向平均的二维数值模式为主［1－5］.二维数值模式不能真实的反映风生流产生的垂向上的显著变化，从而会影响模拟油膜扩散的精度.近年来，三维模式也逐渐开始运用于油膜扩散的计算中［6－10］，三维模式能提供较为真实的风生流垂向变化，能更为准确地模拟表层流场对溢油漂移扩散的作用.本文基于改进的三维数值模式 ECOM－si［11］，加入油膜计算模块［12，13］，计算鸭绿江大桥发生溢油事故后油膜的输移.

2 数值模式设置和验证

2.1 模式设置

模式的计算区域包括整个鸭绿江河口，上边界位于荒沟水文站，外海开边界西至123.6°E，东至124.6°E，南到39.4°N（见图2）.模式最高网格分辨率为35 m，网格拟合岸线且具有良好的正交性和光滑性.水深通过数字化高分辨率海图及相应海区1分精度的水深数据插值获得.垂向均匀分为6层，时间步长设置为5 s.模式采用北京－54坐标系，后处理和绘图时转为经纬度坐标系.

图2 模式计算区域、网格及浪头港局部网格Fig.2 Model calculation domain and grids with detail grid near Langtou Port

考虑潮汐、径流和风应力对油膜漂移扩散的作用.外海开边界以潮位驱动，考虑16个主要分潮 M2、S2、N2、K2、K1、O1、P1、Q1、U2、V2、T2、L2、2N2、J1、M1和 OO1，由各分潮调和常数合成给出（调和常数来自全球潮汐模式NAOTIDE，资料来源网站为http：／／www.miz.nao.ac.jp／staffs／nao99／index＿En.html）：

其中，ζ为潮位，f为节点因子，Hi为振幅，ωi为角频率，gi为迟角，Vi＋ui为订正角（它可以由具体的年、月、日求得）.本文假设夏季和冬季突发溢油事故后油膜的漂移扩散情况，因主要取水口均位于大桥的上游（见图1），考虑到对上游取水口的最坏影响，假定溢油时间均为落憩时刻.上游开边界条件以荒沟水文站实测多年月平均径流量给出，分别取8月和2月为夏季和冬季的代表月份.鸭绿江8月平均径流量为1 400 m3／s，2月平均径流量为700 m3／s，夏季径流量明显大于冬季.风场根据多年月平均风资料，夏季为南风4.6 m／s，冬季为西北风4.9 m／s.初始条件水位、流速均取为零.鸭绿江口及其上游河段存在较大滩涂（潮间带），模式采用干湿判断法处理潮滩动边界.

模式运行约7 d，水位和流场稳定后，开始计算溢油.溢油事故假设瞬时溢油100 t，油种为重质原油，密度为950 kg／m3，低于海水密度.溢油的漂移轨迹采用质点跟踪的方法计算，溢油质点的移动由水流、风应力和随机扩散决定.本文模拟溢油扩散时间尺度较短，暂未考虑挥发、乳化、光化学氧化分解、沉积和生物降解等长时间尺度的物理化学过程.每个质点的下一时刻位置按以下公式计算：式中，X0、Y0为油粒质点的初始坐标；U、V 分别为X、Y方向的流速分量；W10为海面上方10 m处的风速；A为风向；α为修正系数；r为随机扩散项，r＝RE，R为0～1之间的随机数，E为扩散系数；B为随机扩散方向，B＝2πR.只要质点数足够多，油膜漂移的模拟就越精确.本文设置了一百万个油膜质点.每个质点代表一个恒定的体积大小，每步根据网格的水平面积大小和所含质点数，求出油膜在该网格的厚度值：

2.2 模式验证

采用2000年5月6日到5月7日鸭绿江口门外实测站点流速流向资料验证数值模式，测站位置见图1.从实测资料看，鸭绿江河口落潮历时远长于涨潮历时，涨潮历时短（见图3），这与以往研究结果一致［14］.从模式验证结果看，总体上模式计算的流速和流向与实测值吻合良好（见图3）.误差主要是由地形、底摩擦、开边界条件误差和模式计算误差等因素造成.验证结果表明本文建立的数值模式能较好地模拟鸭绿江河口的水动力过程.

3 油膜漂移扩散的数值模拟

3.1 夏季

鸭绿江河道内为不规则半日潮，落潮历时远长于涨潮历时，涨潮历时短［14］.假设溢油发生在落憩转涨时刻，2 h后已基本处于涨憩状态，故考虑溢油后2 h涨潮流期间油膜的漂移扩散过程.夏季在潮流、径流1 400 m3／s和南风4.6 m／s（记为数值试验A）作用下，油膜向上游漂移扩散，在溢油后2 h时刻油膜分布见图4左图.此时溢油已经经过四道沟取水口，影响到丹东造纸厂取水口，丹东造纸厂取水口附近油膜厚度超过55μm.3个取水口油膜厚度随时间变化（实线）见图5.在涨潮过程中，油膜在溢油发生后0.9 h到达四道沟取水口，之后油膜厚度迅速上升至最大值125μm，然后迅速降低，影响持续至1.5 h.丹东造纸厂取水口约在1.5 h受到油膜影响，略超过2 h时油膜厚度达到最大值160μm，持续影响至3.2 h.溢油没有漂移扩散到燕窝取水口.统计夏季20μm和1μm厚度油膜的面积、相对于溢油发生点最远和最近距离随时间的变化，见表1.在0.5 h至2 h内，20μm、1μm油膜距离溢油点的最远距离不断增大，其中2 h时刻最远距离分别为8.93 km、10.21 km.

图3 垂向平均流速和流向模式计算值（实线）和观测值（点）的比较测站Y1（左上），测站Y2（右上），测站Y3（左下）和测站Y4（右下）Fig.3 Comparison of the vertically averaged current speed and direction between the mode1 calculated ones（solid line）and observed ones（point）in the station Y1（left upper），Y2（right upper），Y3（left lower），and Y4（right lower）

图4 落憩后2 h油膜厚度分布图Fig.4 Thickness of oil film contour at 2 hours after ebb slack

图5 落憩后6 h内取水口油膜厚度随时间变化Fig.5 Temporal variation of oil film thickness at the water intakes inner 6 hours after ebb slack.

表1 夏、冬季落憩后2h内油膜面积、漂移最远和最近距离Tab.1 The area，the longest and shortest distance of oil film inner 2 hours after ebb slack in the case of summer and winter

3.2 冬季

冬季在潮流、径流700 m3／s和西北风4.9 m／s（记为数值试验B）作用下，油膜向上游漂移扩散，在溢油后2 h时刻油膜分布见图4右图.此时溢油已经经过四道沟取水口，影响到丹东造纸厂取水口，丹东造纸厂取水口附近油膜厚度接近45μm.3个取水口油膜厚度过程线（虚线）见图5.在涨潮过程中，油膜在溢油发生后1.2 h到达四道沟取水口，之后油膜厚度迅速上升至最大值75μm，然后迅速降低，影响持续至1.8 h.丹东造纸厂取水口约在1.7 h受到油膜影响，2.2 h时油膜厚度达到最大值50μm，持续影响至4.2 h.溢油没有达到燕窝取水口.统计冬季20μm和1μm厚度油膜的面积、相对于溢油发生点最远和最近距离随时间的变化，见表1.在0.5至2 h，20μm、1μm油膜距离溢油点的最远距离不断增大，其中2 h时刻最远距离分别为7.86 km、9.19 km.

3.3 冬夏季差异

比较夏季和冬季溢油扩散的数值模拟结果，可以看出在涨潮过程中，夏季油膜到达四道沟取水口的时间比冬季早0.3 h，四道沟取水口最大油膜厚度冬夏季分别为125μm和75 μm；夏季油膜到达丹东造纸厂取水口的时间比冬季早0.2 h，丹东造纸厂取水口最大油膜厚度冬夏季分别为160μm和50μm（见图5）.图中油膜厚度在丹东造纸厂取水口达到峰值的时刻略超过2 h，这说明2 h内是一个持续的涨潮过程.

为了直观的显示冬夏季油膜向上输运的差异，将冬季涨潮2 h时刻油膜厚度减去夏季涨潮2 h时刻油膜厚度，得到油膜厚度差值图（见图6）.从图中可以看出，丹东造纸厂取水口附近及以上区域油膜厚度变化为负值，绝对值最大超过15μm，说明冬季油膜上溯比夏季更近；四道沟取水口到接近丹东造纸厂取水口的区域油膜变化为正值，说明油膜厚度在该河段冬季比夏季厚，因为冬季向上游漂移扩散距离比夏季短，而冬夏季溢油总量相同，满足质量守恒的原理.

然而夏季径流量（1 400 m3／s）远大于冬季（700 m3／s），径流越大对油膜向上输运的阻碍作用越强，理应油膜向上游漂移扩散夏季比冬季近.为什么夏季油膜的上溯距离反而更远？影响油膜漂移扩散的动力因子除径流外，还有潮汐和风应力.下面对油膜漂移扩散的动力机制作进一步的模拟和分析.

图6 落憩2 h后冬夏季油膜厚度差值图（数值试验B－A）Fig.6 Difference of oil film thicknesses at 2 hours after ebb slack（numerical experiment B minus A）

4 机制分析

4.1 潮汐的影响

本节设置一数值试验，模式动力因子为2月大潮，径流为1 400 m3／s，南风4.6 m／s（记为数值试验A1）.试验A1与试验A区别仅在于潮的不同.将两个数值试验计算的涨潮两小时后的油膜厚度相减（试验A1减去试验A），可以看到四道沟取水口到丹东造纸厂取水口之间油膜差值为正，丹东造纸厂取水口附近及上游油膜差值为负（见图7）.这说明在风、径流相同的条件下，夏季大潮溢油上溯距离比冬季大潮溢油上溯更远.

根据2007年潮汐表，丹东新港（位置见图1）8月大潮期间3日平均潮差为5.64 m，2月大潮期间3日平均潮差为4.95 m（见图8）.夏季潮差比冬季潮差大，潮汐产生的涨潮流更强，说明8月涨潮流对油膜的上溯距离更远.

图7 落憩2 h后油膜厚度差值图（数值试验A1－A）Fig.7 Difference of oil film thickness at 2 hours after ebb slack（numerical experiment A1 minus A）

图8 大潮期间丹东新港潮位及平均潮差（单位：米）Fig.8 The tidal level and mean tidal range at Dandong during spring（Unit：m）

4.2 风的影响

本节考虑风速和风向的变化对油膜上溯距离的影响.风对油膜的影响前人已做过室内试验［15－17］，说明风对油膜的漂移起着重要作用.为了突出风的影响并简化讨论，潮汐均取为8月大潮，径流量均为1 400 m3／s.风取无风、南风9.2 m／s和西北风4.6 m／s，分别命名为数值试验A0、A2和A3，计算结果与A0的计算结果比较.

4.1 风速的影响

从不同南风风速下的油膜厚度差值图可以看到（见图9），无风相对于南风4.6 m／s相同厚度的油膜漂移扩散的距离变近，南风9.2 m／s相对于南风4.6 m／s相同厚度的油膜漂移扩散的距离变远.这说明南风风速的加大有利于油膜向上游漂移扩散更远的距离.

本文使用油膜最远距离差来表示不同条件下油膜漂移扩散结果的差别.表2定量说明了不同速度的南风对20μm（可代表油膜主体）和1μm（可代表油膜边缘）的油膜漂移扩散距离的影响，可以看出南风风速对油膜主体及边缘向上游运移都有促进作用.由表2可知，溢油发生1 h后，油膜边缘的最远距离差绝对值迅速增大；溢油发生1.5 h后，油膜主体的最远距离差绝对值显著增大.这是因为溢油发生后，油膜向上游漂移扩散，上游河道变为近南北走向，使得南风对油膜向上游的促进作用加强.

表2 溢油后2h内20μm和1μm油膜的最远距离差Tab.2 The longest distance difference of oil film at 2 hours after oil spill km

图9 落憩2 h油膜厚度差值图（左：数值试验A0－A；右：A2－A）Fig.9 Difference of oil film thicknesses at 2 hours after ebb slack（left：numerical experiment A0－A；right：numerical experiment A2－A）

4.2 风向的影响

将数值试验A3计算结果与试验A计算结果作比较，得出2 h内油膜上溯的距离试验A3要小于试验A（见图10）.这表明西北风相比南风，阻碍了油膜的向上输移.表3定量说明了不同风向（西北风和南风）对20μm和1μm的油膜漂移扩散距离的影响.由表3可知，西北风与南风相比，1.5 h时刻油膜主体的最远距离差绝对值迅速增至0.94 km.由表2可知，南风和无风相比，1.5 h时刻油膜主体的最远距离差绝对值仅为0.22 km.可以看出，西北风在溢油前1.5 h对油膜向上游的阻碍作用远大于南风对油膜向上游的促进作用.这是因为溢油漂移扩散初期河段为西北走向，导致西北风对油膜向上游阻碍作用明显.

表3 溢油2h内20μm和1μm油膜的最远距离差Tab.3 Difference the longest distance of oil film inner 2 hours after oil spill

图10 落憩后2 h油膜厚度差值图（数值试验A3－A）Fig.10 Difference of oil film thicknesses at 2 hours after ebb slack（numerical experiment A3 minus A）

鸭绿江公路大桥溢油事故发生后，油膜漂移扩散主要由径流、风和潮流共同控制.风的变化对油膜上溯的作用（见图9左图）明显大于冬夏季潮的变化对油膜上溯的作用（见图7），是造成夏季油膜上溯距离远的主要原因.

5 结论

基于三维河口海岸海洋模式ECOM－si，加入油膜输运模块，计算并分析鸭绿江浪头港附近瞬时溢油100 t情况下，2 h内油膜的漂移扩散情况，结论如下.

在涨潮流的作用下，夏季溢油发生后约0.9 h油膜到达四道沟取水口，持续时间约为0.6 h，取水口油膜厚度最大值为125μm；丹东造纸厂取水口约在1.5 h受到油膜影响，略超过2 h时油膜厚度达到最大值为160μm，油膜持续影响至3.2 h.冬季溢油发生后约1.2 h油膜达到四道沟取水口，持续时间约为0.6 h，取水口油膜厚度最大值为75μm；丹东造纸厂取水口约在1.7 h时受到油膜影响，2.2 h时取水口油膜厚度达到最大值50μm，油膜持续影响至4.2 h.冬夏季燕窝取水口均未受到溢油的影响.尽管鸭绿江夏季径流量大于冬季，存在着油膜向上游漂移扩散夏季比冬季更快和更远的现象.

夏季潮汐比冬季强是上述现象的次要原因，风是首要原因.鸭绿江大桥上游河道先是西北走向，然后是近南北走向，南风促进南北走向河道溢油向上游输运，西北风阻碍西北向河道溢油的上游输运.本文研究结果对发生溢油后的应急处理工作具有指导作用和实践意义.

［1］娄安刚，王学昌，孙长青，等.胶州湾海面溢油轨迹的数值模拟［J］.黄渤海海洋，2001，19（1）：1－8.

［2］张和庆，李福娇.近海海面油类漂流扩散的研究和预测实践［J］.热带气象学报，2001，17（1）：83－89.

［3］孙长青，赵可胜，郭耀同.渤海湾海面溢油数值计算［J］.海洋科学，2003，23（11）：63－67.

［4］ PERIANEZ R，PASCUAL－GRANGED A.Modelling surface radioactive，chemical and oil spills in the Strait of Gibraltar［J］.Computers ＆ Geosciences，2007，34：163－180.

［5］ VIOLEAU D，BUVAT C，ABED－MERAIM K，et al.Numerical modelling of boom and oil spill with SPH［J］.Coastal Engineering，2007，54：895－913.

［6］ WANG S D，SHEN Y M，GUO Y K，et al.Three－dimensional numerical simulation for transport of oil spills in seas［J］.Ocean Engineering，2008，35：503－510.

［7］刘彦呈，殷佩海，林建国，等.基于GIS的海上溢油扩散和漂移的预测研究［J］.大连海事大学学报，2002，28（3）：463－465.

［8］熊德琪，杜川，赵德祥，等.大连海域溢油应急预报信息系统及其应用［J］.交通环保，2002，23（3）：5－7，24.

［9］娄厦，刘曙光.溢油模型理论及研究综述［J］.环境科学与管理，2008，33（10）：33－37.

［10］李大鸣，陈海舟，付庆军.海上溢油数学模型的研究与应用［J］.哈尔滨工程大学学报，2008，29（12）：1291－1297.

［11］朱建荣，朱首贤.ECOM 模式的改进及在长江河口、杭州湾及邻近海区的应用［J］.海洋与湖沼，2003，34（4）：245－257.

［12］陈义中，朱建荣，陈昞睿，等.渤海原油码头溢油漂移扩散的数值模拟［J］.华东师范大学学报：自然科学版，2006（6）：25－33.

［13］傅利辉，朱建荣，陈昞睿，等.崇启大桥建成后流场变化及溢油的数值模拟.华东师范大学学报：自然科学版，2008（4）：114－123.

［14］白凤龙，高建华，汪亚平，等.鸭绿江口的潮汐特征［J］.海洋通报，2008，27（3）：7－13.

［15］季荣，陈国华，胡雅蓓，等.流—风—波共存下溢油漂移的实验室模拟［J］.青岛海洋大学学报：自然科学版，1996，26（3）：353－360.

［16］郭运武，刘栋，钟宝昌，等.风对河道溢油扩展、漂移影响的实验研究［J］.水动力学研究与进展：A辑，2008，23（4）：446－452.