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

殷 健，陳志錚，梁珊珊

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

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

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

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

1 長(zhǎng)江口區(qū)域的溢油風(fēng)險(xiǎn)分析

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

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

2 數(shù)學(xué)模型的動(dòng)力

2.1 長(zhǎng)江口水動(dòng)力數(shù)學(xué)模型

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

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

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

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

圖2 模型水動(dòng)力率定站點(diǎn)分布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 長(zhǎng)江口溢油數(shù)學(xué)模型

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

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

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

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

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

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

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

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

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

2.3 長(zhǎng)江口船舶溢油模擬

2.3.1 5·18 溢油事故模擬

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

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

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

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

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

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

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

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

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

2.3.2 溢油模擬的影響因子

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

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

表1 長(zhǎng)江口溢油事故案例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 下各時(shí)段油膜遷移范圍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 下各時(shí)段油膜遷移范圍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 下各時(shí)段油膜遷移范圍Fig.18 Scope of Film under Condition 3 with Time

圖20 工況4 下各時(shí)段油膜遷移范圍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 下各時(shí)段油膜遷移范圍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 下各時(shí)段油膜遷移范圍Fig.24 Scope of Film under Condition 6 with Time

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

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

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

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

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

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

3 結(jié)論

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

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

［2］陳莎.長(zhǎng)江口船舶溢油模擬及其對(duì)環(huán)境的風(fēng)險(xiǎn)研究［D］.上海:上海海洋大學(xué)，2008.

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

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

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

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

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

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

［9］劉彥呈，殷佩海，林建國(guó)，等. 基于GIS 的海上溢油擴(kuò)散和漂移的預(yù)測(cè)研究［J］.大連海事大學(xué)學(xué)報(bào)，2002，28(3):41-44.

［10］郭運(yùn)武，劉棟，鐘寶昌，等.風(fēng)對(duì)河道溢油擴(kuò)展、漂移影響的實(shí)驗(yàn)研究［J］.水動(dòng)力學(xué)研究與進(jìn)展A 輯，2008，23(4):446-452.

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