長江口寶鋼碼頭溢油事故油膜漂移擴散數值模擬

顧靖華, 朱建榮, 金智

河口海岸學國家重點實驗室(華東師范大學), 上海 200241

碼頭、港口、航道和橋墩等鄰近水域是溢油事故多發之處。溢油事故發生后, 油膜的漂移擴散會對周圍水域和敏感目標造成污染。因此, 研究溢油事故后的油膜漂移擴散, 可為溢油事故的處理提供油膜面積、到達敏感目標時間和油膜厚度變化等數據, 對于事故應對和環境保護具有實際意義。

研究溢油事故發生后油膜運動的主要方法有現場實驗、水槽試驗、理論分析和數值模擬。溢油事故后油膜漂浮于水面, 其運動主要是由水流的平流,以及風引起的漂移和隨機擴散決定, 還受蒸發、分散、乳化、沉降及生物降解等過程的影響。張和慶 等(2001)在珠江口通過現場溢油實驗, 分析了油類入海后漂流和擴散過程。郭運武 等(2008)通過水槽實驗, 分析了在不同風速和水流速度條件下油膜的擴展、漂移運動狀況, 結果表明風在溢油油膜擴展、漂移中起著重要的作用。劉棟 等(2006)通過潮汐水槽模擬實驗, 觀察和記錄油膜在潮汐水流中兩個半日潮潮周期內的擴展、漂移情況, 分析了溢油油膜在潮汐水流中的擴展、漂移過程。婁廈 等(2008)理論分析了油在水中的行為和歸宿, 概括了多個溢油模型的理論依據、基本假設和應用范圍。Violeau 等(2007)通過室內試驗和應用平滑粒子水動力拉格朗日方法模擬溢油的漂移擴散。在溢油運動的研究和預測中, 數值模擬是一種應用最為廣泛的研究方法。較早的一些模式應用垂向平均的二維模式較多,例如婁安剛 等(2001)建立了膠州灣內的二維預報潮流場, 再考慮風場對海面溢油的影響, 模擬了溢油事故后溢油的油膜漂移軌跡。孫長青 等(2003)建立渤海灣二維溢油飄移數值模型, 計算在不同風場和潮流場的共同作用下, 溢油飄移的輸移軌跡、可能抵岸的時間、地點及殘留量等。李大鳴 等(2008)建立了在潮汐、風浪作用下油膜運動二維數學模型,對大亞灣海區油碼頭瞬時溢油事故后的影響范圍進行了模擬。Periá?ez 等(2008)建立了二維水動力模式,模擬直布羅陀海峽潮流和余流, 并質點跟蹤技術模擬潮流、混合和風對溢油運動的影響。垂向平均的二維數值模式不能真實反映風生流和密度流等在垂向上的變化, 從而影響模擬油膜漂移擴散的精度。因此, 將三維模式運用到油膜漂移擴散的模擬中,能較為真實地模擬風生流和密度流等的垂向變化,并更為準確地模擬表層流場對溢油漂移擴散的作用。例如, Wang 等(2008)采用粒子方法在二維溢油模式的基礎上發展三維溢油模式, 應用到渤海溢油事故的數值模擬。應用三維數值模式, 陳義中 等(2006)模擬了渤海原油碼頭溢油漂移擴散, 傅利輝等(2008)模擬了長江口北支崇啟大橋建成后流場變化及溢油漂移擴散, 陳士謙 等(2012)模擬了鴨綠江公路大橋溢油的漂移擴散。

長江口為三級分汊四口入海的大型潮汐河口,環境敏感保護目標眾多(表1)。環境敏感保護目標位于溢油事故發生點的不同方向, 不同的方向對應不同的敏感目標。數值預測油膜運動需要考慮主導風和不利風, 其對油膜的運動起著十分重要的作用。對溢油事故附近氣象站實測風作統計分析, 可得出各個方向的頻率和平均風速, 第一頻率的風取為主導風。不利風向取決于溢油事故發生點和敏感目標的連線, 與指向敏感目標最接近的風向為不利風向。主導風和不利風可從實測風的統計結果給出。本文采用崇明東灘氣象站實測的風資料(測站見圖1), 給出主導風和不利風向。在此基礎上, 應用長期改進和應用的河口海岸水動力三維數值模式

表1 長江口環境敏感保護目標Tab. 1 Environmental sensitive protection targets in the Changjiang Estuary

ECOM-si(semi-implicit estuarine, costal and ocean model), 耦合溢油漂移擴散模塊, 模擬并分析上海港寶鋼原料區域碼頭船舶發生溢油事故后油膜的漂移擴散(溢油事故點位置見圖1), 以及對環境敏感保護目標(圖2)的影響。

圖1 長江口形-勢圖及上海港寶鋼原料區域碼頭溢油事故泄漏點Fig. 1 Sketch map of the Changjiang Estuary, and the site of oil leakage accident at the Baosteel raw material regional wharf of the Shanghai Port. The green dot marks the weather station at the Chongming eastern shoal; the red dot represents the oil leakage point, black dots stand for the water intakes of the reservoirs

圖2 長江口環境敏感保護目標位置Fig. 2 Positions of the environmental sensitive protection targets in the Changjiang Estuary

1 數值模式設置

采用作者課題組長期改進和應用的河口海岸水動力三維數值模式ECOM-si, 該模式在長江口水動力過程和鹽水入侵等方面研究中經過了大量的率定和驗證, 得到了廣泛的應用, 取得了大量研究成果

(Lyu et al, 2018; Chen et al, 2019; Zhu et al, 2020; Ma et al, 2022) 。在ECOM-si 的基礎上, 耦合油膜漂移擴散模塊。該模塊將油分成有限個粒子, 每個粒子的位置按以下公式計算:

式中X0、Y0為某粒子的上一時刻坐標;X、Y為某質點的下一時刻坐標;U、V分別為X、Y方向的流速分量, 由水動力數值模式計算得出;W10為海面上10m 處的風速;A為風向;α為風漂移系數。在風的影響下, 油膜漂移速度的增加量約為風速的2%～3%,漂移方向與風向成0°～40°夾角。由于油膜中心的漂移速度和方向是表面水流和風所引起的流速矢量之和, 在模擬中使用風漂移系數α為風與海流的角度修正系數, 一般在0.02～0.03 之間, 本文取0.03;r為隨機擴散項,r=RE,R為0～1 之間的隨機數,E為擴散系數, 計算公式為E=, 其中Aα為水平擴散系數, 由水動力模型提供;B為隨機擴散方向, 計算公式為B=2πR。

在油膜的漂移擴散過程中, 還需考慮油膜的蒸發。油膜蒸發受溢油種類、氣溫、溢油面積、風速、太陽輻射和油膜厚度等因素的影響。本模塊考慮溢油的蒸發公式為(Stiver et al,1984)

其中,Fv為溢油蒸發體積分數;T是環境溫度(空氣溫度), 本文長江口冬季取283K;θ為無因次時間(蒸發暴露時間);A1、B1、T0、TG為各種依賴于油品的常數,A1=6.3,B1=10.3;T0為油在Fv=0 時的初始沸點溫度(單位: K);TG為蒸餾曲線的坡度(單位: K);T0、TG根據下列經驗公式計算:

其中API 為比重指數, 對重油、重油和輕油可分別取值為16.9162、26.8458 和36.1882。蒸發公式中的θ由下式計算:

其中K是空氣端傳質系數,A2是溢油面積,t是時間,V是溢油泄漏的體積。針對公式實際使用中蒸發過快的現象, Fingas(1995)對K值進行了調整, 為了解決無風速情況下傳質系數公式無法計算的問題, 孫寶楠 等(2010)假設在傳質系數上加一個小量R, 從而傳質系數可以表示為

其中U10為海表面 10m 風速。R1的范圍在(0.1～10)×10–4m·s–1的量級(本文K值取10–4),t為溢油發生后的時間, 單位為h。每個油粒子初始時刻的體積(總溢油體積/劃分總粒子數)減去t時刻的溢油蒸發量, 即為t時刻每個油粒子的體積。

溢油事故發生點為上海港寶鋼原料區域碼頭石灰石及砂巖卸船碼頭, 坐標31°28′6″N, 121°25′33″E,泄漏量1000t, 0#輕質柴油。在徑流作用下, 油膜向海方向輸運, 落憩時刻之后為漲潮, 相對于漲憩時刻,油膜會有更長時間滯留在口門內, 為不利情況, 故假設溢油事故發生在大潮期間的落憩時刻。水動力模式先運行10d, 水位和流場準穩定后, 再釋放油粒子。只要質點數足夠多, 油膜漂移的模擬就越精確, 因此本文設置了一千萬個油膜質點。每個質點代表一個恒定的體積大小, 每步根據網格的水平面積大小和所含質點數, 可以求出油膜在該網格的厚度值: 厚度=質點數×t時刻單個質點代表的體積/網格的水平面積。

模型計算范圍包括整個長江河口、杭州灣和鄰近海域(圖3a), 東至125°E 附近, 北至接近34°N, 南至27°30′N 左右。模型水平方向上采用曲線非正交網格, 較好地擬合了長江河口的岸線, 并對南北支分汊口、深水航道等區域加密(圖3b、c)。長江口內網格分辨率最高可至約100m, 長江口外網格分辨率最高可至10km 左右。模型在垂直方向上采用σ坐標, 由表至底均勻分成10 層。計算時間步長取為60s。由于長江河口潮灘分布較廣, 模型采用干濕網格判斷實現潮灘動邊界, 最小臨界水深閾值設置為0.2m。岸線和水深數據為2017 年長江河口實測資料,水深數據以黃海85 高程為基面, 見圖1。

圖3 模型計算范圍和網格(a)以及局部放大的南北支分汊口附近網格(b)和南支口門附近網格(c)Fig. 3 Model calculation domain and grids (a), and the enlarged view of grids near the bifurcation between the South Branch and North Branch (b) and near the mouth of the South Branch. The red lines are the guide dikes and groins of the Deepwater Channel

本文模擬和分析冬季1 月徑流量和主導風、不利風作用下油膜的漂移擴散。選取大通作為上游邊界, 模式的徑流邊界條件可直接采用大通水文站的月平均徑流量資料。1950—2020 年多年月平均長江入海徑流量1 月為11700m3·s–1。外海開邊界潮位考慮由16 個分潮(M2, S2, N2, K2, K1, O1, P1, Q1, MU2,NU2, T2, L2, 2N2, J1, M1和OO1) 的調和常數合成得到, 資料來自全球潮汐模式NAOTIDE(http://www.miz.nao.ac.jp)計算結果。

對主導風和不利風向, 本文采用長江口崇明東灘氣象站(測站位置見圖1)2007 年至2021 年14 年的實測風速風向資料, 統計1 月8 個羅盤方向風的頻率和平均風速。1 月第一至第八風向頻率分別為26.9%、26.1%、16.5%、8.3%、7.7%、6.3%、5.4%、2.8%, 風向分別是N、NW、NE、E、SE、W、S、SW, 風速分別是5.6、4.8、5.0、4.1、4.1、3.6、4.0、3.4m·s–1。圖4 為1 月風玫瑰圖, 形象給出了8 個方向風的頻率和平均風速, 以及各個方向上不同風速(間隔2m·s–1)及其對應的頻率。

圖4 崇明東灘氣象站2007—2021 年統計得出的1 月風玫瑰圖Fig. 4 Wind roses in January calculated from 2007 to 2021 at weather station on Chongming eastern shoal

主導風取為第一頻率的風。1 月第一頻率(26.9%)風向為北風, 風速5.6m·s–1。不利風由溢油事故發生點和環境敏感保護目標位置確定。由于長江口環境敏感保護目標眾多(圖2), 需要給出多個不利風向。上海港寶鋼原料區域碼頭船舶發生溢油事故發生點的上游敏感目標為陳行水庫、太倉水庫、東風西沙水庫、刀鱭核心區和實驗區, 不利風向為東南風, 為第5 頻率(7.7%), 風速4.0m·s–1。對下游東南向敏感目標為青草沙水庫、刀鱭實驗區、吳淞濕地和九段沙濕地, 不利風向為西北風, 為第2 頻率(26.1%),風速4.8m·s–1。對下游東北向敏感目標為青草沙水庫、刀鱭實驗區、崇明東灘保護區和中華鱘保護區,不利風向為西南風, 為第 8 頻率(2.8%), 風速3.2m·s–1。

2 油膜漂移擴散的數值模擬

本文模擬和分析上海港寶鋼原料區域碼頭綜合改造工程碼頭溢油事故發生后油膜的輸運擴散及其對敏感目標的影響, 對模擬結果繪制油膜分布和水源地取水口油膜厚度隨時間變化過程圖進行分析。

2.1 主導風北風5.6m·s–1 風速

油膜的運動主要受水流的平流、風的漂移和隨機擴散作用。本文假設寶鋼碼頭的溢油事故發生在大潮落憩時刻, 長江河口為半日潮, 事故發生后3h內為漲潮流, 在潮流的平流作用下向上游運動(圖5)。油膜中心位于陳行水庫水源地保護區, 厚度超過2μm。之后在徑流的作用下, 油膜在隨漲落潮流上下移動過程中向下游輸運, 至12h 油膜已進入南港,中心位于外高橋附近(圖6)。至72h, 油膜繼續向下游漂移擴散, 大部分位于南港中下段和南槽中上段的南側, 小部分進入北槽南側(圖7)。在油膜的運動過程中, 絕大部分位于長江口南岸, 在主導風北風的作用下向南漂移起著重要的作用。

圖5 溢油事故后3 小時油膜分布Fig. 5 Distribution of oil film thickness after 3 hours of the oil spill accident

圖6 溢油事故后24 小時油膜分布Fig. 6 Distribution of oil film thickness after 24 hours of the oil spill accident

圖7 溢油事故后72 小時油膜分布Fig. 7 Distribution of oil film thickness after 72 hours of the oil spill accident

從長江口水庫取水口油膜厚度隨時間變化過程看(圖 8), 因受徑流和北風的作用, 溢油事故對南支上端北側的東風西沙水庫取水口無影響。太倉水庫位置南支中段南側, 離溢油事故發生點較遠, 同樣受徑流和北風的作用, 油膜最大厚度小于 1μm, 影響微小。青草沙水庫取水口位于溢油點東北側不遠處, 受北風作用, 油膜最大厚度小于 1μm, 影響微小。陳行水庫取水口位于溢油點上游不遠處, 油膜2h 隨漲潮流到達取水口, 之后隨著潮流周期性起伏, 最大厚度約為 68μm。31h后, 油膜的影響消除, 期間最長持續影響10h, 總影響時間為29h。

圖8 油膜厚度在東風西沙水庫取水口(a)、太倉水庫取水口(b)、陳行水庫取水口(c)和青草沙水庫取水口(d)隨時間變化情況黑線: 主導風(北風); 紅線: 不利風向(東南風); 綠線: 不利風向(西北風); 紫線: 不利風向(西南風)Fig. 8 Temporal variation of oil film thickness at the water intakes of the Dongfengxisha reservoir (a), Taicang reservoir (b),Chenhang reservoir (c) and Qingcaosha reservoir. Black line: prevailing wind (northerly wind); red line: unfavorable wind direction (southeasterly wind); green line: unfavorable wind direction (northwesterly wind); purple line: unfavorable wind direction (southwesterly wind)

對青草沙飲用水水源保護區, 溢油事故發生后11h 油膜到達, 36h 后離開, 持續影響25h, 期間最大厚度為1.42μm。對陳行飲用水水源保護區, 溢油事故發生后 1.5h 油膜到達, 53h 后離開, 持續影響51.5h, 期間最大厚度為151.31μm。對長江太倉瀏河飲用水水源保護區, 溢油事故發生后3.5h 油膜到達,5.5h 后離開, 持續影響 2h, 期間最大厚度為2.99μm。溢油事故對東風西沙飲用水水源保護區都無影響。

對長江刀鱭國家級水產種質資源核心保護區,溢油事故發生后2.5h 油膜到達, 29h 后離開, 持續影響26.5h, 期間最大厚度為47.84μm。對長江刀鱭國家級水產種質資源實驗保護區, 溢油事故發生后0.5h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為52.99μm。對上海吳淞炮臺灣國家濕地公園, 溢油事故發生后8h 油膜到達, 持續影響時間超過第72h, 期間最大厚度為32.70μm。對上海九段沙濕地國家級自然保護區, 溢油事故發生后44.5h 油膜到達, 影響時間超過溢油發生后的72h,期間最大厚度為1.47μm。溢油事故對上海崇明東灘鳥類國家級自然保護區、上海市長江口中華鱘省級自然保護區無影響, 緣于它們遠離溢油點, 位于東北的北港口門北側。

2.2 不利風向東南風4.0m·s–1 風速

在不利風向東南風4.0m·s–1作用下, 油膜朝西北方向漂移, 聚集于南支北岸, 同時受徑流作用向下游輸運,在溢油事故后72h 油膜位于南支的上下段的中北側和北港上段的中北側(圖9)。這與主導風北風作用下油膜的分布十分不同, 表明風在油膜漂移擴散中起著十分重要的作用, 迫使油膜克服徑流朝風向的方向漂移。

圖9 溢油事故后72h 油膜分布Fig. 9 Distribution of oil film thickness at 72 hours after the oil spill accident

從長江口水庫取水口油膜厚度隨時間變化過程看(圖8), 對東風西沙水庫取水口, 溢油事故發生后16h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為10μm。與主導風北風結果相比, 東南風使得油膜向西北漂移更遠的距離。對太倉水庫取水口, 溢油事故發生后4h 油膜到達, 8h 后離開, 持續影響4h, 期間最大厚度為1.8μm。對陳行水庫取水口, 溢油事故發生后2h 油膜到達, 12h 后離開, 持續影響10h, 期間最大厚度為7μm。對青草沙水庫取水口, 溢油事故發生后10h 油膜到達, 39h 后離開,持續影響29h, 期間最大厚度為2.3μm。

對青草沙飲用水水源保護區, 溢油事故發生后9h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為10.45μm。對陳行飲用水水源保護區,溢油事故發生后1h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為149.46μm。對長江太倉瀏河飲用水水源保護區, 溢油事故發生后3h 油膜到達, 55.5h 后離開, 持續影響52.5h, 期間最大厚度為52.27μm。對東風西沙飲用水水源保護區, 溢油事故發生后17.5h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為37.89μm。

對長江刀鱭國家級水產種質資源核心保護區,溢油事故發生后2h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為84.47μm。對長江刀鱭國家級水產種質資源實驗保護區, 溢油事故發生后0.5h 油膜到達, 影響時間超過溢油發生后的72h, 期間最大厚度為110.00μm。溢油事故對上海吳淞炮臺灣國家濕地公園、上海崇明東灘鳥類國家級自然保護區、上海九段沙濕地國家級自然保護區、上海市長江口中華鱘省級自然保護區均無影響。這與油膜的平面分布圖是一致的。

2.3 不利風向西北風4.8m·s–1 風速

溢油事故后72h 油膜分布與主導風北風作用下的分布相似, 這里不再給出。

從長江口水庫取水口油膜厚度隨時間變化過程看(圖8), 西北風阻止了油膜向上游漂移, 溢油事故對東風西沙、太倉水庫取水口無影響, 對陳行水庫和青草沙水庫影響微小, 最大油膜厚度小于1μm。與主導風北風結果相比, 西北風與北風作用下陳行水庫取水口油膜厚度差別顯著。

對青草沙飲用水水源保護區, 溢油事故發生后3h 油膜到達, 23h 后離開, 持續影響20h, 期間最大厚度為5.12μm。對陳行飲用水水源保護區, 溢油事故發生后10h 油膜到達, 12.5h 后離開, 持續影響2.5h, 期間最大厚度為4.90μm。對長江太倉瀏河飲用水水源保護區和東風西沙飲用水水源保護區無影響。

對于其他環境敏感保護目標, 西北風減輕了溢油事故點上游和北港、崇明東灘外側敏感目標的影響, 但加重了南港和南槽的環境敏感保護目標的影響。

2.4 不利風向西南風3.2m·s–1 風速

溢油事故發生后, 油膜受漲潮流的作用向上游運動, 同時受西南風的作用向東北漂移, 大部分油膜繞過青草沙水庫的西南角, 進入北港。在溢油事故后72h(圖10), 油膜受向海的徑流和西南風的作用, 大部分油膜在北港向下游漂移擴散, 位于崇明島南岸的北港中段中北側, 小部分油膜在北槽上端的中北側?？梢? 西南風使油膜漂移至北港, 并和徑流一起作用向海輸運。

圖10 溢油事故后72h 油膜分布Fig. 10 Distribution of oil film thickness at 72 hours after the oil spill accident

從長江口水庫取水口油膜厚度隨時間變化過程看(圖8), 溢油事故對東風西沙水庫和太倉水庫取水口無影響。對陳行水庫取水口, 溢油事故發生后2h油膜到達, 8h 后離開, 持續影響6h, 期間最大厚度為6.3μm。對青草沙水庫取水口, 溢油事故發生后9h 油膜到達, 54h 后離開, 持續影響45h, 期間最大厚度為8.5μm。

對青草沙飲用水水源保護區, 溢油事故發生后8h 油膜到達, 67h 后離開, 持續影響59h, 期間最大厚度為17.79μm。對陳行飲用水水源保護區, 溢油事故發生后1h 油膜到達, 41h 后離開, 持續影響40h,期間最大厚度為122.60μm。對長江太倉瀏河飲用水水源保護區, 溢油事故發生后3.5h 油膜到達, 6.5h后離開, 持續影響3h, 期間最大厚度為15.02μm。溢油事故對東風西沙飲用水水源保護區無影響。

對其他環境敏感保護目標的影響, 西南風減輕了溢油事故點上游和下游南港、南槽環境敏感目標的影響, 但加重了北港和崇明東灘外側敏感目標的影響。

3 結論

本文應用河口海岸三維水動力模式ECOM-si,耦合溢油漂移擴散模塊, 模擬和分析長江口寶鋼碼頭船舶發生溢油事故后油膜的漂移擴散, 以及對環境敏感保護目標的影響?；陂L江口崇明東灘氣象站2007—2021 共14 年的實測風速風向資料, 給出影響溢油漂移的主導風和不利風向。在冬季多年平均1月徑流量11700m3·s–1和主導風、3 個不利風作用下,數值模擬和詳細給出了上海港寶鋼原料區域碼頭溢油事故發生后油膜的平面分布、油膜到達和離開4個水庫取水口和飲用水水源保護區以及其他環境敏感保護區的時間、持續影響時間和最大油膜厚度。

長江口寶鋼碼頭溢油事故發生后, 油膜的運動主要受水流的平流、風的漂移和隨機擴散作用。在1 月主導風北風5.6m·s–1風速作用下, 油膜沿長江口南岸向下游漂移擴散, 在溢油事故后第72h 大部分位于南港中下段和南槽中上段的南側, 小部分進入北槽南側。在油膜的運動過程中, 絕大部分位于長江口南岸, 在主導風北風的作用下向南漂移起著重要的作用。在不利風向東南風4.0m·s–1風速作用下,油膜西北方向漂移, 聚集于南支北岸, 同時受徑流作用向下游輸運, 在溢油事故后72h 油膜位于南支的上下段的中北側和北港上段的中北側。在不利風向西北風4.8m·s–1風速作用下, 西北風減輕了溢油事故點上游和北港、崇明東灘外側敏感目標受到的影響, 加重了對南港和南槽的環境敏感保護目標的影響。在不利風向西南風3.2m·s–1風速作用下, 在溢油事故后72h 大部分油膜在北港向下游漂移擴散,位于崇明島南岸的北港中段中北側, 小部分油膜在北槽上端的中北側。西南風減輕了溢油事故點上游和下游南港、南槽環境敏感目標的影響, 加重了北港和崇明東灘外側敏感目標的影響。油膜隨漲落潮流作上下游來回振蕩過程中, 徑流使油膜向海輸運,風使油膜朝風向方向漂移。不同風向作用下油膜的分布和對環境敏感保護目標影響顯著不同, 表明風在油膜漂移擴散中起著十分重要的作用。