長江口區(qū)域溢油事故對青草沙水源地影響的數(shù)值模擬

殷 健，陳志錚，梁珊珊

(1. 上海市水文總站，上海 200232;2. 上海市排水管理處，上海 200001)

據(jù)統(tǒng)計，每年全球平均發(fā)生千噸級以上的溢油事故4.4 起，百噸級溢油事故近百起［1］，我國年均發(fā)生500 起海洋溢油事故。導致海洋溢油事故的最主要原因是海上運輸，每年因航運產(chǎn)生的石油泄露達160 ～200 萬t［2］。近年來部分沿海地區(qū)的海水含油量已超標2 ～8 倍［3］，海洋石油污染形勢十分嚴峻。

溢油的風險概率與水域狀況、天氣條件、船只密度、航道管理等因素相關(guān)。溢油事故發(fā)生后，需要迅速、準確地預測溢油的變化趨勢，為科學有效的決策提供依據(jù)。

溢油的數(shù)值模擬指基于先期計算的溢油水域流場，結(jié)合實地實時的溢油品種、溢油量、風、浪、鹽度、水溫、氣溫等因素，預測溢油在遷移與風化過程中的狀態(tài)、組分、性質(zhì)的改變以及最終歸宿，為應急決策的制定、清除手段的選擇和溢油損害的評估提供科學依據(jù)。

1 長江口區(qū)域的溢油風險分析

長江口地處長江東西運輸通道與海上南北運輸通道的交匯點，是上海港及長江干線的通海咽喉，在內(nèi)陸航運及海路航運中占據(jù)重要地位，是我國通航密度最高的水域之一，年貨運吞吐量超過7 億t，日均船舶量達2 300 多艘，高峰時期更可達6 000 多艘［2］，并以大型船舶為主，發(fā)生溢油事故的潛在風險較大。

同時長江口也是上海最重要的水源地青草沙水庫［4］所在地，區(qū)域內(nèi)有崇明東灘鳥類自然保護區(qū)、中華鱘自然保護區(qū)、九段沙濕地自然保護區(qū)等全國重要的生態(tài)保護區(qū)，研究本區(qū)域船舶溢油對溢油應急處置具有重要的現(xiàn)實意義。

2 數(shù)學模型的動力

2.1 長江口水動力數(shù)學模型

建立高效穩(wěn)定、精確可靠的水動力數(shù)學模型是進行海洋溢油行為與歸宿數(shù)值模擬的重要前提。MIKE21HD 包括了廣泛的水力現(xiàn)象，可有效模擬由于各種作用力而產(chǎn)生的水位及水流變化，進而為溢油模擬模型MIKE21SA 提供可靠的水動力計算基礎。模型計算區(qū)域及水深分布如圖1，模型區(qū)域包括長江口、杭州灣以及鄰近海區(qū);采用非結(jié)構(gòu)網(wǎng)格;模型計算時間步長為30 s，網(wǎng)格節(jié)點為26 121 個，單元數(shù)為47 174 個，計算時采用高階求解格式。

圖1 模型計算區(qū)域及水深分布Fig.1 Model's Region and Depth Distribution

模型邊界條件分別對應徑流、潮汐、海表面風應力及床底摩擦應力。模型主要參數(shù)包括床面阻力系數(shù)、紊動粘性系數(shù)、引潮勢等，其中床面阻力系數(shù)反映水流與床面相互作用過程中床面邊界的形態(tài)及粗糙程度對水流阻力的綜合影響［5］，采用實測資料反推獲得;紊動粘性系數(shù)的選擇對岸線變化復雜且有回流產(chǎn)生的岸段十分關(guān)鍵，采用經(jīng)驗公式獲取。

于模擬區(qū)域內(nèi)選取不同站點進行水動力率定，如圖2 所示。模型能較為準確地模擬區(qū)域內(nèi)潮汐潮流運動，如圖3 ～圖6 所示。

圖2 模型水動力率定站點分布Fig.2 Model's Hydrodynamic Calibration Station

圖3 高橋站潮位過程Fig.3 Gaoqiao Station's Tide Process

圖4 大戢山站潮位過程Fig.4 Dajishan Station's Tide Process

圖5 中浚大潮潮流流向過程Fig.5 Zhongjun Tide's Flow Process

圖6 北港大潮潮流流速過程Fig.6 Beigang Tide's Velocity Process

2.2 長江口溢油數(shù)學模型

2.2.1 溢油過程的數(shù)值模擬

溢油在海洋上的擴散受其物理、化學、生物過程的影響，而這些過程又與溢油的性質(zhì)、水動力環(huán)境、氣象環(huán)境等因素密切相關(guān)［6］。海面溢油發(fā)生后，首先受引力、表面張力、慣性等影響，以薄膜形式散開;隨后在海流、波浪、風場的共同作用下，產(chǎn)生漂移，對油膜漂移的精確模擬是迅速有效地清除溢油，將污染降至最低程度的重要前提;石油烴的蒸發(fā)是溢油質(zhì)量傳輸?shù)闹饕^程，增加油的比重、傾點、粘性及表面張力;油膜在遷移與風化的同時還發(fā)生溶解過程和光氧化過程，光氧化產(chǎn)物雖濃度低但毒性強，長期積累危害較大;由于溢油對水面的覆蓋，減少了溢油層以下水體與空氣的交換和循環(huán)，易導致大面積水體缺氧，進而對水體環(huán)境及水生生物正常生長產(chǎn)生不利影響。海洋溢油行為與歸宿，如圖7 所示。海洋溢油過程時間尺度，如圖8 所示。

MIKE21SA 溢油模型基于歐拉-拉格朗日理論體系，通過對溢油在水體中的擴展、漂移、蒸發(fā)、分散、沉降、乳化、溶解和光氧化等過程的模擬，獲取油膜的漂移過程及厚度增減，同時計算溢油的比重、傾點、粘性等性質(zhì)的變化。

溢油模型基于油粒子模式進行數(shù)值模擬，油粒子模式是將溢油視為大量質(zhì)量不等的粒子所組成的集合，并以粒子的宏觀運動來表達溢油在水環(huán)境中行為過程的一種模擬方法［7］。油粒子模式將溢油離散為大量油粒子，每個油粒子代表一定油量，借助油粒子的隨機運動反映油膜的擴展與漂移，通過油粒子的質(zhì)量改變體現(xiàn)油膜的蒸發(fā)、乳化、分散和沉降等過程，并通過計算設定域內(nèi)油粒子的數(shù)量、體積獲取油膜的厚度分布。油粒子模式可更有效地響應海洋水動力條件，較為精確地模擬油膜運動過程中順風拉長、迎風壓縮以及扭曲斷裂等變化。

圖7 海洋溢油行為與歸宿Fig.7 Behavior and Fate of Marine Oil Spill

圖8 海洋溢油過程時間尺度Fig.8 Marine Oil Spill Process

溢油模型首先計算各油粒子的組分和位置，分析含水率與體積，隨后將運算區(qū)域劃分網(wǎng)格并于各時間步長末統(tǒng)計域內(nèi)油粒子的數(shù)量和分布，依據(jù)油粒子體積及網(wǎng)格面積計算油膜的厚度與濃度，在模擬溢油漂移與風化的同時，借助熱量遷移計算油膜的熱容量、表面張力、粘度等物理化學性質(zhì)的變化。

2.2.2 模型構(gòu)建與參數(shù)設置

溢油模型的運算基于水動力模型，兩者的計算區(qū)域相同。溢油模型網(wǎng)格分辨率取100 m，模型模擬的時間步長小于網(wǎng)格間距與計算域最大流速的比值，模擬時長為30 h。模型考慮的條件場及參數(shù)因子包括溢油特征(源通量、油粒子數(shù)量、油組分特征、油粘度);水體特征(水溫場、鹽度場);大氣特征(氣溫場、陰暗度);擴散(橫向水平擴散、縱向水平擴散、垂向擴散);對流(風場、風因子、風偏轉(zhuǎn)角、風深度);渦流(對數(shù)速度分布);蒸發(fā)(蒸發(fā)率);溶解夾帶(傳質(zhì)系數(shù)、油包水界面張力);乳化(油最大含水率、吸收系數(shù)、釋放系數(shù));熱量遷移(油初始溫度、油輻射率、水輻射率、大氣輻射率、漫射率)。

2.3 長江口船舶溢油模擬

2.3.1 5·18 溢油事故模擬

2012 年5 月18 日夜間，通銀6 加油船裝載重油于吳淞口東長江口6#錨地進水沉沒，事發(fā)位置121°36'10.46″E、31°23'31.75″N，距離上游青草沙水庫取水泵閘約17 km，存在對青草沙水源地造成影響的可能。本次事故為瞬時溢油，溢出約25 m3重油，通過模擬溢油發(fā)生后30 h 內(nèi)油膜位置、面積、形狀、厚度等變化和溢油蒸發(fā)、溶解、乳化、沉降等趨勢(如圖9 ～圖12)，與溢油事故監(jiān)測結(jié)果進行對比，初步探討溢油的遷移與風化規(guī)律。

圖9 模擬30 h 后油膜分布Fig.9 Film Distribution after 30 Hours'Simulation

圖10 各時段油膜遷移范圍Fig.10 Scope of Film with Time

圖11 油膜平均厚度變化曲線Fig.11 Thickness of Film with Time

圖12 油膜10 μm 厚度包絡面積變化曲線Fig.12 Area of 10 μm Film Cover with Time

以下3 個模擬成果表明所構(gòu)建的油粒子模式溢油模型能較好地反映溢油的遷移與風化規(guī)律，適用于長江口區(qū)域的溢油模擬。

(1)油膜擴展主要受溢油自身物理化學性質(zhì)、潮流場、風場等因素的影響。在溢油發(fā)生后初期，油膜在重力及慣性力作用下，覆蓋面積逐步增大，平均厚度迅速減小，這一階段油膜呈橢球形均勻擴展，隨著時間推移，潮流逐漸成為支配油膜擴展的主導因素，油膜在海面被潮流拉成彗星狀油膜帶，并隨潮汐來回漂移。

(2)油粒子模式中，油膜在海面的漂移受表層流場及表面風場驅(qū)動，考慮到二維水動力模型計算的潮流為垂向平均值，模型建立對數(shù)函數(shù)關(guān)系通過垂向平均流速推算表層流場，并借助數(shù)值積分獲取油粒子團的運動軌跡。潮流與風應力的矢量疊加決定了油膜的漂移方向、速度、距離，油膜的漂移隨表層流速的增大而加速，隨表層流速的減小而減速或轉(zhuǎn)向，在漲、落急以及漲、落憩等潮流特征時段表現(xiàn)得尤為顯著，與此同時表面風場的存在有效促使了油膜向下風方向漂移。溢油遷移時段與范圍表明，溢油發(fā)生后30 h 內(nèi)油膜未對青草沙水庫取水泵閘產(chǎn)生影響，與溢油事故監(jiān)測結(jié)果一致。

(3)溢油在進行擴展、漂移的同時發(fā)生蒸發(fā)與乳化。前期蒸發(fā)隨油膜面積的增大而增強，在易蒸發(fā)的輕組分快速蒸發(fā)后，溢油的蒸發(fā)過程逐漸減弱，蒸發(fā)速率逐步減緩趨于穩(wěn)定。在波動產(chǎn)生的混合能量作用下，溢油不斷乳化，含水率增加，有效體積及覆蓋面積增大，約16 h 后含水率的變化趨緩，并逐漸穩(wěn)定在80%左右。溢油的蒸發(fā)及乳化現(xiàn)象導致溢油的密度變化，油膜密度因蒸發(fā)過程造成的輕組分損失以及乳化過程引起的含水率增加而呈上升趨勢，并隨時間推移與周圍水體的密度逐漸接近。

2.3.2 溢油模擬的影響因子

環(huán)境是一個開放性系統(tǒng)，自然因素及人類活動都會對環(huán)境產(chǎn)生各種復雜的影響，這些影響往往難以精確定量，數(shù)值模擬中模型結(jié)構(gòu)、輸入條件、參數(shù)取值等也存在一定的不確定性［8］。借助影響因子分析建立低靈敏度系統(tǒng)，可有效提高模型運行的準確性，通過比對不同輸入條件下的模擬結(jié)果，掌握各種狀況下油膜輸移運動的規(guī)律，從而可在事故發(fā)生后有針對性地收集模擬預測所需的關(guān)鍵信息。

溢油數(shù)值模型主要用于突發(fā)性溢油事故的模擬預測，目的是為事故的應急響應提供支持，油膜的運動軌跡、溢油影響范圍、油膜到達敏感區(qū)域的時間等是模型模擬預測的重點，這些指標主要受溢油位置、溢油類型、溢油方式、溢油量、潮流場、潮汐形態(tài)、風場、溫度場、鹽度場等因素的影響［9］，本文選取溢油量與風場因素進行影響因子分析，通過比對長江口歷年重、特大溢油事故案例(如表1)，溢油量因子選擇50、100 m3;考慮到長江口最頻繁風向為NW-N 及ESE-SSE，頻率分別達24%和23%，年平均風速3.7 m/s，風場因子選取長江口最頻繁風向以及年平均風速，在此基礎上于事故位置建立工況1 ～工況6(如表2)進行溢油模擬運算(如圖13 ～圖26)，討論不同工況溢油事故對青草沙水源地的潛在影響。

表1 長江口溢油事故案例Tab.1 Oil Spill Case of Yangtze Estuary

表2 溢油模擬工況Tab.2 Simulated Conditions of Oil Spill

圖13 工況1 下30 h 后油膜分布Fig.13 Film Distribution under Condition 1 after 30 Hours'Simulation

圖14 工況1 下各時段油膜遷移范圍Fig.14 Scope of Film under Condition 1 with Time

圖15 工況2 下30 h 后油膜分布Fig.15 Film Distribution under Condition 2 after 30 Hours'Simulation

圖16 工況2 下各時段油膜遷移范圍Fig.16 Scope of Film under Condition 2 with Time

圖17 工況3 下30 h 后油膜分布ig.17 Film Distribution under Condition 3 after 30 Hours'Simulation

圖18 工況3 下各時段油膜遷移范圍Fig.18 Scope of Film under Condition 3 with Time

圖20 工況4 下各時段油膜遷移范圍Fig.20 Scope of Film under Condition 4 with Time

圖21 工況5 下30 h 后油膜分布ig.21 Film Distribution under Condition 5 after 30 Hours'Simulation

圖22 工況5 下各時段油膜遷移范圍Fig.22 Scope of Film under Condition 5 with Time

圖23 工況6 下30 h 后油膜分布ig.23 Film Distribution under Condition 6 after 30 Hours'Simulation

圖24 工況6 下各時段油膜遷移范圍Fig.24 Scope of Film under Condition 6 with Time

圖25 工況1 ～6 下油膜平均厚度變化曲線Fig.25 Thickness of Film under Condition 1 ～6 with Time

通過對比不同工況的模擬成果，分析溢油量及風場因素對油膜遷移與風化過程的影響，具體分析如下。

圖26 工況1 ～6 下油膜10 μm 包絡面積變化曲線ig.26 Area of 10 μm Film Cover under Condition 1 ～6 with Time

(1)溢油量的增加提高了油膜厚度與覆蓋范圍，而風場則在溢油的遷移、風化過程中起顯著作用。靜風環(huán)境下，油膜大體隨潮汐漲落方向漂移，風場的出現(xiàn)促進油膜向下風方向移動，且風速與流速的大小不同、風向與流向的夾角不同，對油膜漂移的影響不同［10］。工況3、4 中漲潮階段風向與流向反向，減緩了油膜向上游漂移的速度，落潮階段風向與流向同向，油膜加速向下游漂移;工況5、6 中漲潮階段風向與流向同向，加快了油膜向上游漂移的速度，落潮階段風向與流向反向，油膜減速向下游漂移，相比靜風條件，各時段油膜位置均偏向西北。

(2)油膜厚度變化歷經(jīng)兩個階段，第一階段油膜厚度迅速減小，持續(xù)約3 h。第二階段油膜厚度隨時間推移緩慢減小，與之相應油膜面積則逐步擴大。相比靜風環(huán)境，風場引起的水流湍動及水面波浪影響驅(qū)動油膜的蒸發(fā)和乳化過程，導致油膜厚度與面積的變化梯度時緩時陡。在潮流場與風場的共同作用下，在工況3、4 的轉(zhuǎn)潮時段，油膜各部分受到方向不同的潮流作用，發(fā)生扭曲及斷裂，并進一步相互分離形成油斑，工況5、6 中油膜抵達岸邊后受局部地形影響顯著變形，模型采用的油粒子模式借助大量獨立運動的油粒子，較好地模擬了油膜的扭曲變形及油斑形成。

(3)溢油事故發(fā)生后，短時間內(nèi)油膜在水中的主要行為為擴展及漂移，因此，建議在突發(fā)性溢油事故發(fā)生后，如若溢油源距離青草沙水庫較近，應重點考慮油膜在水中擴展漂移到達取水泵閘的時間［11］。溢油遷移時段與范圍表明，溢油發(fā)生后30 h 內(nèi)工況1、2、5 和6 油膜未對青草沙水庫取水泵閘產(chǎn)生影響，溢油發(fā)生后21 ～24 h 內(nèi)工況3、4 油膜影響至取水泵閘區(qū)域，期間油膜厚度維持在1 ～10 μm。

3 結(jié)論

運用MIKE21SA 建立長江口杭州灣二維溢油數(shù)值模型，針對5·18 溢油事故進行后評估，與事故監(jiān)測結(jié)果的對比表明，所構(gòu)建的油粒子模式溢油模型能較好反映溢油的遷移與風化規(guī)律，可適用于長江口區(qū)域的溢油模擬。事故位置不同工況的模擬結(jié)果表明溢油量的增加提高了油膜厚度與覆蓋范圍，而風場則在溢油的遷移、風化過程中起顯著作用;油膜厚度變化歷經(jīng)兩個階段，第一階段油膜厚度迅速減小，持續(xù)約3 h，第二階段油膜厚度隨時間推移緩慢減小，與之相應油膜面積則逐步擴大;如若溢油發(fā)生在青草沙水源地附近，應重點以油膜的擴展漂移作用評估其擴散至水源地的時間。

［1］王長海.溢油漂移擴散計算模式初步研究［J］.交通環(huán)保，2000，21(2):7-9.

［2］陳莎.長江口船舶溢油模擬及其對環(huán)境的風險研究［D］.上海:上海海洋大學，2008.

［3］岳文潔，劉書俊.淺析海洋產(chǎn)業(yè)的環(huán)境損耗核算［J］. 環(huán)境與可持續(xù)發(fā)展，2008，(1):44-46.

［4］上海青草沙水源地工程簡介［J］.凈水技術(shù)，2009，28(3):78.

［5］許婷.丹麥MIKE21 模型概述及應用實例［J］. 水利科技與經(jīng)濟，2010，16(8):867-869.

［6］劉偉峰，孫英蘭. 海上溢油運動數(shù)值模擬方法的探討與改進［J］.華東師范大學學報(自然科學版)，2009，(3):90-97.

［7］劉生寶.基于POM 的長江感潮河段溢油數(shù)值模擬研究［D］.南京:河海大學，2007.

［8］姜衛(wèi)星.黃浦江溢油事故的數(shù)值模擬研究［D］. 上海:同濟大學，2007.

［9］劉彥呈，殷佩海，林建國，等. 基于GIS 的海上溢油擴散和漂移的預測研究［J］.大連海事大學學報，2002，28(3):41-44.

［10］郭運武，劉棟，鐘寶昌，等.風對河道溢油擴展、漂移影響的實驗研究［J］.水動力學研究與進展A 輯，2008，23(4):446-452.

［11］龍紹橋，婁安剛，譚海濤，等. 海上溢油粒子追蹤預測模型中的兩種數(shù)值方法比較［J］. 中國海洋大學學報(自然科學版)，2006，36(s1):157-162.