ＩＳＳＮ１３４６－７８４０

港湾空港技術研究所
資料
TECHNICAL NOTEOFTHE PORT AND AIRPORT RESEARCH INSTITUTE

No.1300

March 2015

海上流出油の移流及び拡散に関する数値計算法の開発

松崎

独立行政法人義孝

港湾空港技術研究所

Independent Administrative Institution,

Port and Airport Research Institute, Japan

目

次

要旨

4

1. 序論

5

1.1

背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

既往の研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6(1)

海上流出油の移流及び拡散に関する数値計算法の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6(2)

油膜自身の特性による油拡散に関する研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7(3)

油の乱流拡散に関する研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7(4)

数値計算法の有用性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

本論文の目的と構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3

2. 油拡散粒子モデルの開発及び計算方法

92.1

モデル概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.2

移流外力の計算方法・入力方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.3

油膜自身の特性による油拡散の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.4

油の乱流拡散の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.5

海岸線，構造物等境界における計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.6

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

3. 油膜自身の特性による油拡散の推定方法

12

3.1

概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

3.2

油膜自身の特性による油拡散に関する Fay のモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.3

室内実験による油拡散実験方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.4

室内実験の結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.5

油膜自身の特性による油拡散に関する数値計算モデルの開発 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.6

定常流における油拡散実験と数値計算モデルの比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.7

考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.8

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4. 海水面極近傍における水平乱流拡散係数の推定方法

24

4.1

概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.2

風洞水槽実験による疑似油の選定. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.3

実海域における疑似油を用いた水平乱流拡散実験の方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.4

実海域における疑似油を用いた水平乱流拡散実験の結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.5

水平乱流拡散係数の推定方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

4.6

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

5. 韓国泰安沖 Hebei spirit 号油流出事故の再現計算

29

5.1

事故概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

5.2

計算条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

5.3

計算結果と観測結果の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

5.4

計算に用いるパラメータの検討 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

5.5

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

1

6. 東京湾 Diamond Grace 号油流出事故の再現計算

33

6.1

事故概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

6.2

計算条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

6.3

計算結果と観測結果の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

6.4

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

7. 東京湾 Beagle3 号油流出事故の予測計算

37

7.1

事故概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

7.2

計算条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

7.3

計算結果と観測結果の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

7.4

結語 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

8. 議論

40

8.1

既往の研究からの進展 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

8.2

数値計算法の考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

8.3

即時の予報体制の構築 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

8.4

油回収船の油回収操船支援を行う漂流油捕捉システムの開発 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

9. 結論

43

謝辞

44

記号表

48

A 付録：潮汐流の予測方法

50

A1

潮汐流の数値計算モデル及び境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50A2

潮汐流計算結果からの再現方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

2

A simulation model of oil transport on the water surface

Yoshitaka MATSUZAKI*

Synopsis

A simulation model was developed to predict spilled oil transport on the water surface using particle tracking methodof Port and Airport Research Institute (PARI) for helping oil recovery by oil recovery ships, named OIL-PARI. The OILPARI formulates three major factors that aﬀect the floating oil movement, the oil advection, the mechanical oil spreading,and the oil diﬀusion terms. This paper is intended the latter two terms.
First is the mechanical oil spreading term. We conducted the mechanical oil spreading experiments and the results
were compared with Fay’s model (1969, 1971). As a result, oil slick spreading velocity was well simulated by the Fay’s
third stage model of surface tension - viscous region when the oil slick thickness is smaller than 10 −4 m. The spreading
coeﬃcient of Fay’s third stage model was also derived. Next, we developed a new simulation model to calculate the eﬀectof the net surface tension and oil buoyancy to modify the random walk technique and Fay’s model to meet the particle
expression.
Second is the oil diﬀusion term. A diﬀusion model at the sea surface was constructed to predict the turbulent diﬀusion
coeﬃcient of spilled oil, referring results from in situ drift experiments using pseudo oil made of sponge rubber. And we
discussed the adaptively about derived diﬀusion model.
And also, we discussed the validity of the simulation model. To verify the new model, hindcast simulations were applied
to the three actual oil spill incidents. First, Hebei Spirit oil spill incident, which occurred oﬀshore from the Taean coast inSouth Korea in 2007, was simulated and the result was compared with observed ones of oil gravity position and oil area.
Second, we simulated Diamond Grace oil spill incident, which occurred in Tokyo bay in 1997, we discussed oil advection
area from observed and simulated results. And third, we simulated Beagle 3 oil spill incident, which occurred at the mouthof Tokyo bay in 2014. The error of oil gravity position between observed and calculated results was about 1 - 2 km, andadvection velocities were almost same between them.
Key Words: oil spill, numerical simulation, particle tracking method, hindcast simulation* Senior Researcher, Marine Environmental Information Group, Marine Information and Tsunami Division
3-1-1 Nagase, Yokosuka 239-0826 Japan
Phone: +81-46-844-5065 Fax: +81-46-844-0575 E-mail: matsuzaki-y@pari.go.jp

3

海上流出油の移流及び拡散に関する数値計算法の開発

松崎 義孝*

要旨

本研究では，油回収船の油回収操船支援を目的とした，海上に流出した油の漂流予測を行う油拡散粒子モデ
ル，OIL-PARI を開発した．油拡散粒子モデルは海表面の油に関して，海表面の流れによる油の移流，油膜自身の特性による油拡散，油の乱流拡散を計算できる．本論文において既往の研究から進展させた内容は，油膜自身の特性による油拡散と油の乱流拡散の予測方法についてである．
油膜自身の特性による油拡散について，Fay（1969，1971）のモデルの適応を平面 2 次元油拡散の室内実験結果から議論した．実験において油膜は厚さが 10 −4 m オーダー以下で拡がっており，その際の拡がり速度は Fayの提案する第 3 段階の表面張力-粘性領域のモデルで説明できることを示した．また，従来行われていなかった，Fay の第 3 段階の表面張力-粘性領域における拡がり係数を実験結果から導出した．次に，油膜自身の特性による油拡散について，Fay のモデルに基づいた新しい数値計算モデルを開発した．開発したモデルは，表面張力によって拡がる効果を計算するランダムウォークモデルと，重力によって拡がる効果を計算する斥力モデルの 2 種類である．
次に，乱流拡散について，数値計算モデルを用いた計算で必要となる水平乱流拡散係数について，導出した方法は擬似油を用いた実海域実験を行い，流出油のように海水面極近傍に限った水平乱流拡散係数を導出する方法を開発した．また，導出した拡散係数の計算方法について，その適応範囲を議論した．
また，実際の油流出事例 3 件と数値計算結果の比較を行い，数値計算法の有効性を検証した．まず，2007 年に韓国泰安沖で発生した油流出事故の再現計算を行い，油の移流重心及び拡散面積について比較を行った．また，1997 年に東京湾で発生した油流出事故の再現計算を行い，油の移流位置と観測結果の関係について議論した．さらに，2014 年 3 月に東京湾で発生した油流出事故の予測計算を行い，移流重心が 1〜2 km の誤差に収まっており，移流速度がほぼ一致していることを確認した．

キーワード:油流出，数値計算，粒子モデル，再現計算* 海洋環境・津波研究領域 海洋環境情報研究チーム （併）新技術研究開発領域 油濁対策研究チーム （併）アジア・太平洋沿岸防災センター 主任研究官〒 239-0826 神奈川県横須賀市長瀬 3-1-1 港湾空港技術研究所
電話:046-844-5065 Fax: 046-844-0575 E-mail: matsuzaki-y@pari.go.jp

4

序論
1.1

囲で被害があった．ナホトカ号事故のときは風浪が強く，洋上での油回収作業が困難であった．また流出した油が C 重

背景油であり，流出した油がエマルション化し高粘度となったた

タンカーの海難事故や油井の事故による油流出は従来に比め，油回収作業は困難を極めた（佐尾, 1998）．同じく 1997

較して減っているけれども，依然日本及び世界各地において年の 7 月には東京湾においてダイヤモンドグレース号の油流発生している．図-1.1 は，年ごとの世界のタンカー油流出事

出事故があり，約 1,550 kL の油が流出している．

故の回数（Musk, 2012）を示す．これをみると 1970 年代かタンカー以外の油流出では，2011 年の東日本大震災の際

ら 2000 年代にかけて，大規模な油流出事故の回数は減ってに各地で大きな被害があり，発電所，工場，船舶，自動車等きている．しかしながら，大規模な油流出が完全になくなっからの油流出が多数発生している．消防庁が行ったアンケー

たわけではなく，今後も油流出が発生すると予想される．

ト結果（畑山, 2012）によると，屋外タンク貯蔵所からの油2000 年代に入ってからは，2002 年にスペイン沖で Prestige

流出量は確認されているだけで 46,000 kL を超えている．こ号の油流出事故が発生し，約 63,000 kL もの油が流出した．

れはナホトカ号重油流出事故の油流出量の約 10 倍である．

また，2007 年に韓国西岸の泰安沖で Hebei spirit 号とクレーン

災害時に油流出が発生し港湾内に拡がると，船舶の港への台船の衝突による油流出事故が発生し，約 12,000 kL もの原

入港や係留ができないため，緊急物資の運搬に大きな影響を油が流出した．流出した油は韓国西海岸の広範囲に拡散し，及ぼす．また，油が陸上にまで運ばれると，港湾施設が使用

大きな環境被害を及ぼした．

できない，悪臭，毒性による人体への影響が懸念される等とタンカー以外の油流出では，2010 年にメキシコ湾の MC252

いった問題が発生するほか，火災が発生する可能性がある．

において油井 Deep Water Horizon が爆発炎上し，水深約 1,500

さらに，津波の濁水中の土粒子に重油が付着し，海底に沈降

m の海底から油が流出した．油の流出量は約 779,100 kL とすることで，長期的な環境問題，風評被害を引き起こすと考いわれており世界最大級の油流出事故となった．

えられる．

我が国近海においてもタンカーからの油流出が発生してい今後も日本近海において大規模油流出事故が発生する可能る．1997 年 1 月に日本海でナホトカ号の油流出事故があっ性があると考えられる．図-1.2 は近年の我が国周辺におけるた．油の流出量は約 6,240 kL であり，日本海沿岸の広い範

油汚染発生件数（海上保安庁警備救難部環境防災課, 2014）である．油汚染報告件数は年間約 300 件程度と非常に多い．これらは小規模な油流出がほとんどであるけれども，小さな事

80

故に大きな事故の要因が隠れており，今後も大規模油流出の60リスクがあると考えられる．そのため，油流出に関する十分

40

な準備が必要となる．
油流出災害が発生した際は初動体制の構築，防除計画の策20定が行われ，その際に数値計算を用いた流出油の漂流シミュ
レーションが情報として用いられる．詳細な防除計画の策定

2010

2005

2000

1995

1990

1985

1980

1975

0
1970

Number of oil spill

100のためには油の漂流位置を高精度に予測できることが望ましYearい．また，地震，津波等自然災害による港湾での油流出被害

2013

2012

2011

2010

2009

2008

2007

2006

Year

2005

450
400
350
300
250
200
150
100
50
0
2003

2010

2005

2000

1995

1990

1985

1980

1975

Reported numberof oil pollution

40
35
30
25
20
15
10
5
0
1970

Number of oil spill(a) 油流出量：7〜700kL

2004

1.

Year

(b) 油流出量：700kL 以上

図- 1.2: 近年の我が国周辺海域における油汚染発生件数（海図- 1.1: 世界のタンカー油流出事故の回数（Musk, 2012）

上保安庁, 2014）

5移流(Advection, Drift)
油膜自身の特性による油拡散(Spreading)
油の乱流拡散(Diffusion)

蒸発(Evaporation)
酸化(Oxidation)

乳化(Emulsification)
溶解(Dissolution)

分散(Dispersion)

生物分解(Biodegradation)
堆積(Sedimentation)

図- 1.3: 国土交通省が所有する大型油回収船（新保, 2003）

図- 1.4: 流出油の物理的，化学的，生物的状態変化過程を予測するツールとして，高精度に予測できる数値計算の開漂着リスクの評価の指標等に用いられる．状態変化を予測す発が望まれている．

るモデルは油回収資機材の選定，環境影響評価等の指標に用国土交通省では図-1.3 のような大型の油回収船を所有し，いられる．数値計算においては移流と拡散を予測するモデル

外洋を含む日本近海において油流出が発生した際は出港後 48（例えば，AL-Rabeh, 1994； Wu and Wang, 2010），状態変化

時間以内に油回収作業を行える 3 船体制をとっている．まを予測するモデル（例えば，ADIOS2，http://archive.orr.noaa.

た，東京湾，伊勢湾，瀬戸内海といった内湾を対象とした小gov/book_shelf/538_adios.pdf ），両者を同時に計算するモデル

型の油回収船を所有しており，海上保安庁の要請により油回（例えば，Johansen, 1987； 石油連盟, 1992； Chao 他, 2008；

収作業を行う．そのため，油流出事故が発生した際に流出油

Wang 他, 2008）がある．

の拡散や移流を予測し，効率的に油回収を行うために，流出

移流と拡散を予測するモデルには平面 2 次元の挙動を対象

油の漂流シミュレーションを用いることが期待されている．

としたモデル（例えば，Johansen, 1987；AL-Rabeh, 1994； Wu

また，平時においては油回収船による防除訓練を行っており，

and Wang, 2010）と，水中の挙動を含めた 3 次元の移動を対訓練の想定として流出油の漂流シミュレーションの活用も期象としたモデル（例えば，石油連盟, 1992； Chao 他, 2003；

待されている．

Wang 他, 2008）がある．タンカー事故等のように，油の流出

そこで，著者が所属する独立行政法人港湾空港技術研究所

元が海上にある場合，油の移流及び拡散は平面 2 次元の挙動

では 2007 年度から油回収船による油回収作業の支援を目的が支配的である．

とした，漂流予測に関する数値計算法の開発に取り組んで数値計算を用いた流出油の漂流予測において，平面 2 次元

きた．

1.2

の移動を対象とした海表面での流出油の移流と水平拡散は，

3 つの効果を考慮する必要がある．，1 つ目は海流，潮汐流，既往の研究

風による表面流といった，海表面の流況による移流，2 つ目(1) 海上流出油の移流及び拡散に関する数値計算法の概要

は重力（油と海水の密度差による浮力）や表面張力による油

海上に流出した油は，流れ，風，波等の物理的な過程，気温，水温，紫外線等の化学的な過程，生物分解等の生物的な拡散（Spreading），3 つ目は風や流れの乱れによる油の乱流過程により様々に状態を変化させる．また，油の種類（軽質

拡散（Diﬀusion，渦拡散ともいう）である．
数値計算の手法は，流出油を粒子の集まりと考える Lagrange

油，重油，原油等）の違いによりその過程が異なる． 図-1.4
は流出油の物理的，化学的，生物的な状態変化過程を示す．

的計算手法（例えば，石油連盟, 1992；Chao 他, 2003； Wang

移流，拡散，蒸発，沈降，溶解，エマルション化（乳化，も他, 2008）と，計算領域を格子で区切り，油の位置を Euler 的

しくはエマルジョン化ともいう），堆積，生物分解といったに計算する手法（例えば，Tkalich 他, 2003）がある．第 2 章
で述べるとおり，移流及び拡散の効果はそれぞれ速度とし変化が生じる．

て求めて，それらをベクトル合成して油の移動速度を計算

流出油の物理的，化学的，生物的な状態変化過程の予測すする方法が合理的で，計算方法として有効であり，その場合るモデルとして，主に移流と拡散を予測するモデル（Transport

Lagrange 的計算手法が適している．

model）と，蒸発，分散，溶解，エマルション化，堆積，生流出油を粒子の集まりと考える Lagrange 的計算手法を用

物分解等の状態変化を予測するモデル（Fate model）がある．

いた数値計算法は，数多く提案されている．既往の Lagrange

拡散と移流を予測するモデルは油回収位置の選定，環境脆弱的計算手法を用いた数値計算法における問題点は，油膜自身性の高い海域でのオイルブームの展張位置の決定，季節別油6

の特性による油拡散，油の乱流拡散の取り扱い，及び数値計

油膜自身の特性による油拡散実験の既往の研究はいくつ算法の有用性の検証であると考えられ，次項以下で詳しく説か行われている（例えば，高橋他, 1979； Akiyama 他, 1983；

明する．

Buist, 1989）けれども，それらはいずれも 1 次元方向を対象としており，平面 2 次元方向の油膜自身の特性による油拡散(2) 油膜自身の特性による油拡散に関する研究に関する実験は著者らの知る限りでは，埜口他 (1979) のみ

油は流れがない静かな水面においても拡がる．静水面上のである．

油は重力（油と水の密度差による浮力）と，正味の表面張力

埜口らは，重油を用いて静水面上での 1 次元及び平面 2 次によって拡がっており，この効果を油膜自身の特性による

元方向の拡がり実験を行い，実験結果と Fay の理論解との比油拡散，もしくは油の拡がりといい，英語では Spreading や較検証を行った．平面 2 次元方向の実験は水平方向に 6×6

Mechanical Spreading という．

m，水深 50〜70 cm の平面水槽において行った．埜口らは，油膜自身の特性による油拡散の再現モデルはいくつか提案(a) 重力‐慣性領域，及び (b) 重力-粘性領域において，実験されている．まず最も有名なものと思われるのは油拡散半径

結果は Fay のモデルと良好に一致し，(c) 表面張力-粘性領域の時間変化を求める Fay（1969, 1971）のモデルであり，多において，実験結果は Fay のモデルと一致しないことを示し

くの数値計算法に採り入れられている（例えば，石油連盟,

た．また，実験において周辺部が薄い層になった際の油膜形

1992； Wang 他, 2005）．Fay は，水面の油は重力と表面張力状の観察が困難で，詳細な運動の機構は検討できなかったとによって拡がり，慣性力，粘性力によって拡がりが抑えられしている．

ると考えた．また，油流出からの時間経過に伴い，それぞれ以上のことから，油膜自身の特性による油拡散に関するの力の大きさの相対的な関係が変化すると考えた．そこで，

Fay のモデルにおける拡がり係数及び数値計算モデルについFay は油が流出し拡がる過程において，支配的な因子は時間て，十分な議論がなされていない状況にある．

経過に伴い (a) 重力‐慣性領域， (b) 重力-粘性領域，(c) 表面
張力-粘性領域の 3 つに分けられるとし，3 つの領域について(3) 油の乱流拡散に関する研究

それぞれ油流出からの経過時間と拡散距離（拡散面積）の関

海には潮汐流，海流，風等の影響による大小様々な乱れ係を導出した．

（渦）が存在し，この乱れにより海上に流出した油は拡散しまた，Fay のモデルを改良し，風の影響を考慮した油拡散モていく．これを乱流拡散といい，英語では Diﬀusion という．

デルが Lehr et al. (1984) から提案されており，こちらも数値

油の乱流拡散の計算について，既往の多くの数値計算法

計算に広く採り入れられている（例えば Chao 他, 2003； Guo

では 2.4 節で説明するランダムウォークモデルを用いている

and Wang, 2009）．

（例えば，石油連盟, 1992； Chao 他, 2003； Wang 他, 2005；

Fay に代表されるモデルは油流出からの経過時間と拡散距

Guo and Wang, 2009）．ランダムウォークモデルを用いた方離の関係を求めているため，粒子モデルを用いた数値計算に法においては水平乱流拡散係数のみが油の拡がりの大きさをおいてどのように計算すればよいか，既往の研究では未解明決定するため，拡散係数は計算結果に大きく影響を及ぼす．

のままである．そのため Fay のモデルを粒子を用いた数値

けれども，流出油の数値計算に用いる係数の設定方法につい計算法に適応するには数値計算モデルが必要となる．また，て，十分な議論がなされていないままであった．

Fay に代表されるモデルにおいて，油は瞬間的に流出すると

例えば ASCE Task Committee on Modeling of Oil Spills of the

考えている．しかしながら，タンカーやパイプラインから油

Water Resources Engineering Division (1996) は水平乱流拡散係流出事故が発生すると，油は連続的に流出するため，瞬間流数を 1 - 100 m2 /s として計算することを提案している．しか出を対象としたモデルをそのまま用いることができない．ましながら，これは設定値により拡散面積の計算結果が 100 倍

た，海上には海流，潮汐流，風の影響による表面流があるた

異なることを意味する．

め，油が流れの中で連続流出するといった複雑な現象についケーススタディを見ると，Chao et al. (2003) はシンガポーて再現できる数値計算モデルであることが求められる．

ル海峡での油流出の再現計算において，水平乱流拡散係数を10 m2 /s として油の乱流拡散を計算している．

また，Fay のモデルで経過時間から拡散距離を求める際には拡がり係数が必要である．Fay は 1 次元方向の油膜自身のしかしながら，他の既往研究には具体的な水平乱流拡散係特性による油拡散については実験結果と比較して拡がり係数数の記載が見つからなかった．よって油拡散を精度よく計を求めているけれども，軸対象（平面 2 次元）方向について算するために，拡散係数を適切に見積もる方法が求められては実験を行っておらず，理論的に拡がり係数を求めている．

いた．

そのため，平面 2 次元方向の油膜自身の特性による油拡散に一方で，染料や漂流物を用いた実海域における拡散の観測

関する実験を行う必要がある．

が多く行われており，水平乱流拡散係数の推定方法が提案さ7れている．例えば，Okubo (1971) は染料を用いた拡散実験を予測モデルを開発している．Varlamov らもまた日本海で発行い，拡散スケールと水平乱流拡散係数の関係について整理生したナホトカ号の油流出事故の再現計算を行っている．こ

している．

こでは海流と風を考慮して再現計算を行い，油の移流と拡散

柳・樋口 (1982) は海面下 1 m に抵抗板をつりさげた 60 cmの計算結果を説明しているけれども，海上での分布状況の比四方の耐水ベニヤ板を浮子とし，その位置を航空機で観測す

較は明示されていない．また，他の海域においての適応性にる方法で拡散実験を行っている．

ついては明らかではない．

道田他 (2009) は海面からそれぞれ 1 m，6 m，11 m，16 m

Wang et al. (2008) は 1990 年 6 月に渤海で発生した油流出の位置にドローグ (抵抗体) をとりつけた GPS 搭載漂流ブイ

事故との比較計算を行っている．ここでは潮汐流と風を考慮の放流実験を行っている．

して再現計算を行い，油の移流と拡散について観測結果と比

これらの実験は，いずれも水面から数 m における水平乱較している．

流拡散係数を推定している．しかし，油拡散のような海水面

Proctor et al. (1994) は 1991 年 1 月に Arabian Gulf において極近傍に限った水平乱流拡散係数の推定は従来行われてい発生した油流出に関して観測結果と数値計算の比較を行ってなかった．道田らは水平乱流拡散の水深依存性を指摘しておいる．

り，特に海面付近では風などの影響が大きいと考えられるた

Chao et al. (2003) はシンガポール海峡周辺で発生した油流

め，水平乱流拡散係数が海水中より大きくなることが予想さ出事故との比較計算を行っている．ここでも潮汐流と風を考れる．

慮して再現計算を行い，油の移流と拡散について観測結果と以上のことから，海上流出油の乱流拡散を予測するための比較している．

ランダムウォークモデルを用いた計算に必要となる，海水面

Guo and Wang (2009) は 2005 年に発生した中国大連沖の油極近傍の水平乱流拡散係数について，その導出方法が求めら流出の再現計算を行っており，観測結果と計算結果はよく一れている．

致している．
しかしながら，実際の油流出の再現計算及び精度検証を複(4) 数値計算法の有用性

数の海域や油流出事例について行っている既往研究は見当た

既に紹介しているように，海上に流出した油の拡散，移流らず，上述したように 1 つの油流出の再現計算を行っているに関する数値計算法は多くの事例がある（例えば，石油連盟,

か，再現計算を行わずモデルを提案するのみにとどまってい

1992； Proctor 他, 1994； Varlamov 他, 1999； Chao 他, 2003；

る．海上での油流出事故は日本近海において内湾，外洋のいCastanedo 他, 2006； Wang 他, 2008； Sotillo 他, 2008；Guo and

ずれにおいても発生する可能性がある．そのため，開発する

Wang, 2009）．

数値計算法は，油回収船での作業に必要となる日本近海のい既往の数値計算法の概要と，実際の油流出事故の予測計算

ずれにおいても計算ができることが求められる．すなわち，や再現計算に関して，以下の事例が挙げられる．

特定の海域，油流出事故のみでなく，複数の海域及び油流出

石油連盟 (1992) は 1991 年から流出油の漂流予測を行う

事故を同じ数値計算法を用いて，同様の計算方法で再現でき流出油の拡散・漂流予測モデル を開発しており，現在に至ることを確認することによって，有用性を示すことができるるまで改良を行っている．流出油の挙動について，移流，油と考える．

膜自身の特性による拡散，流れの乱れによる乱流拡散，蒸発，鉛直分散，乳化，海面下での挙動について予測対象とし1.3 本論文の目的と構成ている．予測モデルを初めて扱うユーザーでも操作が簡易で本論文は，以上のような背景と既往の研究を踏まえて，海あり，非常に使いやすい GUI となっている特徴がある．石

上流出油の移流及び拡散に関する数値計算法を開発すること

油連盟は本モデルを用いて 1997 年 1 月に日本海で発生したを目的とする．

ナホトカ号の油流出事故の再現計算を行っている．その際に基本方針について，本論文で開発する数値計算法は，油回海流と風を考慮して再現計算を行い，油の移流と拡散の計算

収船において流出油の効率的な回収を支援するためのものと結果を説明している．特に海上での分布状況の比較は行われする．そこで，開発する数値計算法は以下の通りとする．まず，府県ごとの油回収量を示している．また，1997 年 7 月

ず，海上流出油の移流及び拡散のみを予測対象とするため，に東京湾で発生したダイヤモンドグレース号油流出事故（第油の挙動は平面 2 次元を対象とする．予測海域は油回収船の6 章参照）の流出油漂流シミュレーションを実施しているけ

油回収作業範囲となりうる，日本近海の閉鎖性海域及び外洋

れども，計算結果と観測結果との比較は行われておらず，計とする．予測期間は油流出から 1〜2 日後について精緻に予算結果が示されているのみである．

測できることとする．また，流出油の漂流予測を行う数値計

Varlamov et al. (1999) は日本海を対象とした流出油の漂流

算は Lagrange 的に行う方がよいと考えられるため，数値計

8

算法の構築には Lagrange 的な計算手法を採用する．

計算に用いるパラメータについて検討を行う．

次に (2) 項において示した油膜自身の特性による油拡散に

第 6 章では 1997 年に東京湾で発生した Diamond Grace 号

関する既往の研究の問題点である，Fay（1969, 1971）のモデ

油流出事故を対象に再現計算を行い，観測結果と比較する．

ルのモデルの適応性について検討を行うとともに，Fay のモ

第 7 章では 2014 年に東京湾で発生した Beagle3 号油流出

デルに基づいた数値計算モデルを新たに開発する．

事故を対象に予測計算を行い，観測結果と比較する．

次に (3) 項において示した油の乱流拡散に関する既往の研第 8 章では本研究で得られた海上流出油の移流及び拡散究の問題点である，ランダムウォークモデルにおける水平乱に関する研究成果を考察し，既往の研究から進展させた点に流拡散係数の導出方法について検討する．

ついて述べる．次に，開発した数値計算法について考察を行

さらに，(4) 項において示した数値計算法の有用性についう．また，本論文で提案する油拡散粒子モデルを用いた即時て，開発した数値計算法に関して，2 つの油流出事例に関すの流出油漂流予測体制について説明する．さらに，油拡散粒る再現計算と 1 つの油流出事例に関する予測計算の結果を示子モデルを応用した油回収船における油回収の操船を支援す

し，議論する．

る，漂流油捕捉システムについて説明する．

本論文は 9 章より構成されており，以下では各章の概要を

第 9 章では本論文で得られた研究成果を総括するととも示す．

に，今後の展望について述べる．

第 2 章では海上流出油の移流及び拡散に関する数値計算法

2. 油拡散粒子モデルの開発及び計算方法

である，油拡散粒子モデルの内容と計算方法について説明する．まず海流，潮汐流，風による表面流による油の移流の計

2.1 モデル概要

算方法について説明する．次に，粒子を用いた数値計算モデ

第 1 章で述べた通り，海上に流出した油は流れ，風，波等

ルの概要を示す．次に，水平方向の油の拡がりを支配する油の物理的な過程，気温，紫外線，生物分解等の化学的，生物

膜自身の特性による油拡散と油の乱流拡散について計算方法的な過程により様々に状態を変化させる．今回開発する数値を説明する．最後に，海岸線，構造物等境界における計算方

計算は海表面の流況による移流，油膜自身の特性による油拡法について説明する．

散及び流れの乱れによる油の乱流拡散を予測の対象とする．

第 3 章では油膜自身の特性による油拡散に関する新しい数

海面上の油は流れ，風，波の影響により移流する．移流外値計算モデルを提案する．まず，Fay の油拡散の考え方を用力のうち，流出油の漂流予測に影響が大きいのは海流，潮汐

いて，粒子間の斥力を考慮した数値計算モデルを新たに導出

流及び風による表面流である．海流は定性的には予測できるする．次に，平面水槽を用いた室内実験により油膜自身の特ものの，漂流予測の観点では定量的な予測は難しい．潮汐流性による油拡散実験を実施し，Fay の式との比較を行うととは往復流で，移流のほか油の拡散範囲に大きく影響を与える．

もに，拡がり係数を求める．また，油拡散に関する支配的な海上における風速は陸上と比較して大きいため，風によって力について考察する．次に，室内実験結果を踏まえて，ラン

発生する表面流は，油の移流と拡散に大きく影響を与える．

ダムウォークの考え方に基づいた新しい数値計算モデルを構

油拡散に与える波の影響は，波の質量輸送による移流，水中築する．最後に，定常流条件下における油の連続流出拡散実分散等が考えられるが，今回のモデルでは波の影響を考慮し験を行うとともに，数値計算の結果と比較し，数値計算モデていない．

ルの妥当性について検証する．

流出した油は重力，表面張力といった油膜自身の特性によ

第 4 章では油の乱流拡散に関して，従来は検討されていなる油拡散と，流れの乱れによる油の乱流拡散との影響によりかった海水面極近傍の水平乱流拡散係数を導出する方法につ拡がる．英語ではそれぞれ Spreading（もしくは Mechanical

いて検討する．そのために，実海域において疑似油の拡散実

Spreading）及び Diﬀusion と表現される．海上での油流出の験を実施し，実験結果に基づき，数値計算に導入するための場合，流れの乱れによる油の乱流拡散のほうが支配的である

水平乱流拡散係数の導出方法について考察する．

場合が多いけれども，流れの乱れの少ない湾や大量の油流出

第 5 章から 7 章では実際の油流出事故を対象として数値計の場合は油膜自身の特性による油拡散も無視できないと考え算を実施し，観測結果と比較することで，数値計算法の妥当られる．

性について検証する．比較対象とする油流出事故としては，

図-2.1 は油拡散粒子モデルの概要を示している．流出油の比較検討できるデータが整っている事例が望ましいため，観移流及び拡散は，油を粒子の集合体に見立て，ラグランジュ

測結果がある油流出事故を対象とした．まず，第 5 章では的に計算する．すなわち，各々の粒子 i の位置 X i は，移動2007 年に韓国西岸の泰安沖で発生した Hebei Spirit 号油流出

速度 Ui に対して，事故を対象に再現計算を行い，観測結果と比較する．また，

9

表- 2.1: 風係数法を用いた流出油漂流予測モデルと風係数

油拡散粒子モデル
油粒子の流出開発者

年

風係数 (%)

Wang and Hwang (1974)

1974

3.0

Johansen (1987)

1987

3.0

Shen and Yapa (1988)

1988(1.0〜6.0 で油粒子の移動速度の計算
油拡散•油膜自身の特性による拡散•流れの乱れによる乱流拡散
海表面の流況による移流

油粒子の移動

3.0
変化する)

移流外力
海表面流況
•潮汐流
•海流
•風による
表面流

図- 2.1: 油拡散粒子モデルの概要

∂ Xi
= Ui
∂t

Reed and Gundlach (1989)

1989

3.0

AL-Rabeh (1994)

1994

3.1

Sotillo et al. (2008)

2008

2.6

Wang et al. (2008)

2008

3.0

Guo and Wang (2009)

2009

3.0

Wu and Wang (2010)

2010

3.0

海流及び潮汐流の流速は，海水面において流れに対する油(2.1)の滑りや遅れはないと考えて，移流速度としてそのまま与で表わされる．ここで，t は時間を表す．時間積分は，Euler

える．

陽解法とし，着目時間 t の次の計算ステップにおいて，油の風による移流を計算する方法は 2 通りあると考えられ，一つは風係数法による表面流の推定法，もう一つは 3 次

X it+1 = X int + Ui ∆t(2.2)

元流体計算モデル（例えば，富田・柿沼 (2005)）を用いて直とする．ここで ∆t は計算時間刻みである．この計算方法は接海表面の吹送流を計算する方法である．以下でそれぞれの

1 次精度の陽解法であるため，精度の低い計算手法である．

計算方法について説明する．

しかしながら，数値計算法の精度よりも，風や流れの移流外風係数法は表-2.1 の通り，従来から流出油の漂流予測に多力の再現精度のほうが数値計算結果に与える影響が大きいたく使われてきた手法であり，近年も用いられている．水面近め，数値計算法が 1 次精度であることについて特に議論し傍の風速 W と風による表面流 Uw の関係は，大気から見たない．

海面への接線応力と海面から見た大気への接線応力の釣り合
いから，式 (2.4) のように求められる．

油粒子の移動速度 Ui は水平 2 次元のみを考慮する．粒子の移動速度は海表面の流れに追従し，海表面の流れは各流れ

2
τ = ρa Cda W 2 = ρ w Cdw Uwのベクトル合成で求まると考え，Ui を，

Ui = Uoi + Ut i + Uwi + Usi + Udi + Ubi

ここで τ は接線応力， ρa ， ρ w はそれぞれ大気，水の密度，(2.3)

Cda ，Cdw はそれぞれ大気から海面，海面から大気の抵抗係数である．Cda と Cdw が等しいと仮定すると，と計算する．ここで，Uoi は海流による移流速度，Ut i は潮汐流による移流速度，Uwi は風による移流速度，Usi は油膜√
Uw =自身の特性による油拡散速度，Udi は油の乱流拡散速度，Ubiは境界からの距離に依存して付加した強制流速ベクトルを表ρaW ≒ 0.035W
ρw

(2.5)

となり，水面近傍の風速 W の 0.035 倍の吹送流が発生する

わす．次節以降において，それぞれの取り扱いについて説明

ことがわかる．本計算モデルでは風による油粒子の移流速度

する．

2.2

(2.4)

Uwi は式 (2.6) で計算する．
移流外力の計算方法・入力方法

Uwi = Cw W10i

油拡散粒子モデルにおいて，移流外力は海表面の流れであ(2.6)

る．海表面の流れは対馬海流，日本海流，千島海流といった

ここで，Cw は風係数， W10i は粒子 i の位置における海上

海流，潮汐によって発生する潮汐流及び風による表面流を対10 m での風速を表わす．海面近傍より海上 10 m での風速の象としている．移流外力は数値計算，現地観測等から作成しほうが大きくなることから，本計算モデルでは風係数 Cw をたデータを油拡散粒子モデルに入力する．

0.03 とする．

10

次に，3 次元流体計算モデルである富田・柿沼 (2005) のSTOC-ML（付録 A 参照）を用いた吹送流の直接推定方法について説明する．計算において水面上で式 (2.7) のように風応力 τ を与える．
2
τ = ρa Cda W10(2.7)

ここで，W10 は海上 10 m における風速である．海面抵抗係

y+2Δy
y+Δy

数 Cda は，本多・光易 (1980) の式を用いて，以下のように

Δy y

与える．

y-Δy
y-2Δy

Δx


−3


 (1.29 − 0.024|W10 |) × 10
Cda = 
 (0.581 + 0.063|W10 |) × 10 −3


W10 < 8 m/s
W10 ≥ 8 m/s

x-2Δx

x+Δx

x

x-Δx

x+2Δx

図- 2.2: ランダムウォークの導出の定義(2.8)

風係数法と直接海表面の吹送流を計算する方法を用いた流
出油の移流予測については，第 5 章において検証を行う．

2.3

油膜自身の特性による油拡散の計算方法

t
t
Cx±1,
y = C x, y ± ∆x

油膜自身の特性による油拡散は，重力（油と海水の密度差）
及び表面張力によって拡がる現象のことをいう．本論文で提ttCx,y ±1 = C x, y ± ∆y

案する新しい数値計算モデルは第 3 章で詳しく説明する．

2.4

t
∂Cx,
y

t+∆t
t
Cx,
y = C x, y + ∆t
t
∂Cx,
y

∂x
t
∂Cx,
y
∂y

(2.10)∂t∆x 2 ∂ Cx, y
2
∂x 2
2
t
2
∆y ∂ Cx, y
2 t

+
+

2

∂y 2

(2.11)(2.12)

よって，式 (2.9) は以下のように書ける．

油の乱流拡散の計算方法

海には大小様々な乱れ（渦）が存在し，この乱れにより海

t
∂Cx,
y

上に流出した油は拡散していく．数値計算で求めた流況デー∂t

タのメッシュ間の流況や，観測値の計測点間の流況は不明

 2 t∂ Cx, y= D H 

∂x 2

+

t
∂ 2 Cx,y∂y 2 

(2.13)

ここで，∆=∆x=∆y とすると， D H は以下のように書ける．

なため，線形補間するなどして対応するけれども，メッシュ
間や計測点間の再現できていない流れ（乱れ）によっても油DH =は移流，拡散する．このような，再現できているグリッドス

1 ∆2
4 ∆t

次に，濃度 C の拡散速度を VC を考える．VC は x ，y 方向

ケールの流れのみで再現できない，サブグリッドスケールので等しく，距離 ∆ を時間 ∆t で移動することから，

流れによって生じる油の移流及び拡散を乱流拡散として補間
している．

VC =

本数値計算モデルにおいて，乱流拡散を考慮する際はランダムウォークで計算する．ランダムウォークは拡散方程式を

√

きる．

VC =

2 次元空間を図-2.2 のように，∆x ，∆y の空間格子で定義(2.16)

向における油の乱流拡散速度 Udi は式 (2.16) と等しく，平面

ある時刻 t から次の時間ステップへ進むとき，

4

4D H
∆t

次に，本稿で提案する数値計算法において，平面 2 次元方

後に濃度がとなりに進む確率は 1/4 とする．

ttt
Cxt −1, y + Cx+1,
y + C x, y −1 + C x, y+1(2.15)となる．

し，ある濃度 C が拡散する場合を考える．ある時間から ∆t

(

∆
∆t

と表される．式 (2.14)，式 (2.15) より，VC は，計算していることと同義であり，それは以下の通りに説明で

1
t+∆t
Cx,
y =

(2.14)

方向にランダムに拡がると考えて，

)

√
Udi =

(2.9)

となる．ここで，以下のようにテイラー展開を行う．


4D H cos(2πRr and ) ∆t  sin(2πRr and ) 


(2.17)となる．ここで，Rr and は 0 から 1 の一様乱数である．sin，cos は平均 0，分散 0.5 の乱数であるため，式 (2.17) は式 (2.18)と書き換えることができる．

11

まず初めに，海上に流出した油の挙動を予測するためには，自然外力による移流，重力，表面張力といった油膜自身の特海岸線影響なし

海岸線

海岸線法線方向斥力性による油拡散，油の乱流拡散をそれぞれ速度ベクトルで求め，それらをベクトル合成して再現することを説明した．

Ubi

海岸漂着したとして停止

次に，それぞれの速度ベクトルの計算方法について説明し

rbi

海域

た．まず，移流外力の導入方法について述べた．移流外力は

rb1
rb2

海流，潮汐流，風による表面流を考慮する必要があることを示した．その中で風による表面流の再現方法は 2 通り考えられ，一つは風係数法による表面流の推定法，もう一つは 3 次

陸域
海岸線

元流体計算モデルを用いて直接海表面の吹送流を計算する方

図- 2.3: 海岸線，構造物境界での計算方法法を説明した．風係数法と直接海表面の吹送流を計算する方

√
Udi =

法のどちらが有効であるかの具体的な考察は，第 4 章におい
2D H  Rn1 ∆t  Rn2 

て検証を行う．

(2.18)

続いて，油膜自身の特性による油拡散の速度ベクトルを計算する新しい数値計算モデルについて示した．また，数値計

ここで，D H は海水面極近傍の水平乱流拡散係数 (m2 /s)，Rn1 ，算モデルの詳細については第 3 章で紹介することを述べた．

Rn2 はそれぞれ独立した平均 0，分散 1 の正規乱数である．
√
平均 0，分散 1 の正規乱数を用いて表現するとき， 2 で除さらに，油の乱流拡散の速度ベクトルの計算方法を示した．
また，乱流拡散速度の計算に必要となる海水面極近傍の水平すことで検算できる． D H の求め方は，第 4 章で新しく提案

乱流拡散係数の新しい導出方法については第 4 章で説明するしており，そちらで詳しく説明する．

2.5

と述べた．
最後に，海岸線，構造物等境界における計算方法について

海岸線，構造物等境界における計算方法

海岸線や構造物等を考慮した計算を行えるように，海岸位示した．海岸線，構造物等境界では，境界からの距離に応じ

置の緯度・経度を点情報として入力し，海岸位置（点とそのて強制流速（強制流速ベクトル）を与えて，油の漂着を計算

線分）を漂流境界とした．油粒子は各海岸位置との距離ベクする方法を示した．

トル r bi に応じて，以下のように処理した（図-2.3 参照）．す3. 油膜自身の特性による油拡散の推定方法

なわち，境界から r b1 離れている場合は境界による影響なし（式 (2.19)），境界から r b1 より近づくと境界からの距離に応3.1 概要じた斥力（強制流速ベクトル）を与える（式 (2.20)），境界か海上に流出した油は，海表面の流れや風により油が乱流拡ら r b2 より近づくと海岸あるいは構造物に漂着したとして粒散し，重力（油と海水の密度差による浮力）や表面張力といっ子の移動を停止する（式 (2.21)）．

た，油膜のもつ特性によって拡がる．油膜のもつ特性によって拡がる効果を油膜自身の特性による油拡散，もしくは油の

Ubi = 0Ubi = k br

r bi
|r bi | s

Ui = 0拡がりという．

r bi > r b1

(2.19)

r b1 ≥ r bi > r b2

(2.20)

いて，最も有名なモデルとして Fay（1969, 1971）のモデルがr b2 ≥ r bi

(2.21)

挙げられる．Fay は，油膜自身の特性による油拡散は重力，第 1 章で述べたとおり，油膜自身の特性による油拡散につ

慣性力，粘性力，表面張力によって拡がり速度が決まり，油
ここで，s，k br は入力パラメータである．式 (2.20) は海岸付が流出し拡がる過程で支配的な因子が時間経過に伴い (a) 重

近では波の反射等の影響で，海岸法線に向かって移流しにく

力‐慣性領域， (b) 重力-粘性領域，(c) 表面張力-粘性領域とくなる効果を想定している．また，式 (2.21) は r b2 より油粒分けられるとして，拡がりに関するモデルを導出した．

子が近付いた場合，粒子は漂着し再流出しないとしている．

2.6

Fay のモデルは拡散半径と経過時間の関係について定式化しているため，粒子を用いたモデルに Fay のモデルを当ては結語

めるには数値計算の上で工夫が必要となる．すなわち，Fay

本章は海上に流出した油について，移流及び拡散を予測すのモデルは油が瞬間流出した時からの時間の関数となっているための数値計算法である，油拡散粒子モデルの開発に向けるけれども，タンカーやパイプラインから油流出事故が発生て，その計算方法を示したものである．以下では，本章で得すると，油流出は連続的に発生するため，Fay のモデルをそられた主要な結論を要約する．

12

lのまま用いることはできない．また，油が流れの中で連続流出するといった複雑な現象について，どのように計算すれば

油の移動

油の移動

よいかよくわかっていなかった．

h

油

一方，Fay は後述する拡散半径を求める式 (3.6)，(3.7)，(3.8)の係数を理論的に導出しており（Fay, 1971），2 次元方向の油

の拡がりに関する実験を行っていない．Fay は 2 次元の係数は実海域での油流出事故における油の拡がりと比較して求めているけれども，油の乱流拡散及び移流の効果も含まれるため，再現精度に疑問が残る．

重力
表面張力

そこで本章において，以下のことについて検討を行う．ま

慣性力
粘性力

ず，Fay のモデルについて簡単に紹介する．次に，油膜自身の特性による油拡散に関する実験を行い，Fay のモデルの適応性を評価するとともに，Fay のモデルにおける拡がり係数
図- 3.1: 油膜自身の特性による拡散の定義を評価する．次に，Fay のモデルを基に，数値計算において粒子の拡がり速度をランダムウォークを用いた方法で求めるモデルと，粒子間の斥力で求めるモデルを開発する．最後

性力は経過時間を t ，加速度を l/t 2 とすると，式 (3.2) で表わに，定常流条件下における油拡散の実験結果と，開発した数

される．

値計算モデルによる計算結果を比較し，計算結果の妥当性を

ρw l
t2

検証する．

3.2

(3.2)

ここで，ρ w は水の密度である．粘性力は速度を l/t ，水の動

油膜自身の特性による油拡散に関する Fay のモデル

粘性係数を νw ，水の境界層厚さを δ = (νw t) 1/2 とすると，式

本節では，Fay（1969, 1971）のモデルについて簡単に紹介(3.3) で表わされる．

する．Fay は瞬間的に流出した水面の油膜を円柱として考え，
それが時間経過とともに拡がっていく現象を考えた．その際ρ w νwに，水面の油は重力と表面張力によって拡がり，慣性力，粘1/2 3
ρ w νwll/t=δh
Vt 3/2(3.3)

性力によって拡がりが抑えられると考え，単位体積当たりになお，式 (3.3) から明らかな通り，ここでの粘性は水の粘性かかる力を図-3.1 のように定義し，それぞれモデル化を行っ抵抗を表しており，油の粘性ではない．

た．また，油流出からの時間経過に伴い，それぞれの力の大表面張力は式 (3.4) で表わされる．

きさの相対的な関係が変化すると考えた．そこで，Fay は油

σn
σn l
=
lh
V

が流出し拡がる過程において，支配的な因子は時間経過に伴い (a) 重力‐慣性領域， (b) 重力-粘性領域，(c) 表面張力-粘性(3.4)

ここで σ n は正味の表面張力係数（英語で Net Surface Tension）

領域の 3 つに分けられるとし，3 つの領域についてそれぞれであり，式 (3.5) で表される．

油流出からの経過時間と拡散距離（拡散面積）の関係を導出した．この場合において，Fay は以下の仮定を用いている．

σ n = σ wa − (σoa + σow )すなわち，（1）油の拡がりは方向を持たず一様に水平方向に(3.5)

拡がり，油層厚の分布は考慮しない，（2）油の粘性は水の粘ここで σ wa ，σoa ，σow はそれぞれ水-大気間，油-大気間，性に比べて大きい，（3）油の状態変化（粘性の変化，蒸発，

油-水間の表面張力係数である．
流出開始から約 1 時間後までを第 1 段階と呼ぶことにし，

溶解）はない．
以下で，油に働く力について説明する．まず，重力は油と重力と慣性力が卓越する時間帯とした．重力 (式 (3.1)）と慣水の密度差による浮力のことを表わしており，圧力勾配は式

性力（式 (3.2)）の釣り合いから油層半径 l は式 (3.6) で表わ(3.1) で表わされる．

される．

∆ρgh ∆ρgV=ll3(
l = k1

(3.1)

) 1/4
∆ρ
gV
t 1/2
ρw

(3.6)

ここで ∆ρ は水と油の密度差，g は重力加速度，l は拡散半

それから約 1 週間後までを第 2 段階と呼ぶことにし，重力径， h は油層厚，V は油層の体積で，V = l 2 h で表わす．慣と粘性力が卓越する期間とした．重力と粘性力（式 (3.3)）の釣り合いから，油層半径 l は式 (3.7) で表わされる．

13

カメラ
カメラ
定量油投入器
油
2.5m

30mm

油

定量油投入器

油

油
20m

平面図

1740mm

断面図

図- 3.3: 油回収実海域再現水槽
表- 3.1: 油回収実海域再現水槽の仕様

830mm
6m
図- 3.2: 油膜自身の特性による油拡散実験の模式図（左：小型水槽，右：大型水槽）

規模 (m)

幅 13, 長さ 32, 高さ 3.5

計測水面 (m)

幅 6, 長さ 20, 水深 2.5

造波

最大波高 0.5 m， 周期 1〜5 秒

流速)∆ρ 2 −1/2 1/6 1/4
l = k2
gV νw
tρw

最大 1 m/s

水温調整

(

(3.7)

(2 knots)

5〜30 ℃調整可

使用水

海水

供試油

機械油，重油 (A〜C)

約 1 週間後以降は第 3 段階と呼ぶことにし，粘性力と表面張力（式 (3.4)）が卓越する期間とした．粘性力と表面張力の釣り合いから，油層半径 l は式 (3.8) で表わされる．



σ n2

（1）油投入量を精密に測るため，投入用の油および油投入
用容器の重量を測定する (図-3.5)．

 1/4


l = k 3  2
ρ w νw t3/4（2）内径 26 mm の筒状の油投入器を水槽中央に設置する(3.8)（図-3.6）．設置後，投入容器内に油を注入する．

k 1 ，k 2 ，k 3 は拡がり係数である．Fay は拡がり係数を理論的（3）水槽天井から鉛直下向きに設置したカメラ（Nikon D5200）に求めており，k 1 = 1.14，k 2 = 1.45，k 3 = 2.30 としている．

で動画撮影をスタートする．

3.3（4）動画撮影スタートから約 3 秒後に油投入器を水面を乱室内実験による油拡散実験方法

さないように気を付けて持ち上げる．持ち上げると油投入器

本節では，油膜自身の特性による油拡散に関する数値計算

内の油が時間とともに拡がる．

モデルを開発するための前段階として，油膜自身の特性によ（5）油膜が水槽壁に到達したら動画撮影をストップする．

る油拡散に関する実験を行い，Fay のモデルの適応性を評価（6）投入用の油および油投入用容器の重量を再度測定し，重するとともに，Fay のモデルにおける拡がり係数を評価する．

量の差から油投入量を求める．

実験は水平方向 830×1,740 mm，水深 30 mm の小型水槽

カメラは PC と USB ケーブル接続し，ソフトウェア Cameraと（図-3.2 左参照），長さ 20 m，幅 6 m，水深 2.5 m の大型

Control Pro2 （Nikon 製）で操作した．撮影頻度 30 fps で撮影

水槽（図-3.2 右参照）を用いて行った．なお，大型水槽は油した画像から水槽の外壁の影響が無い範囲のみを解析して油回収実海域再現水槽といい，英語では Simulation Tank for Oil

膜の面積を求める．実験は油種ごとに 4〜10 回行った．使用

Recovery in Marine Situations，略して STORMS という．水槽した油の種類とその物性は表-3.2 の通りである．油の密度はの写真を図-3.3 に，仕様を表-3.1 に示す．油回収実海域再現

図-3.7 の比重計，油の粘度は Brookfield 社の B 型粘度計（型

水槽は平成 16 年度から供用されている，海水仕様としたわが国初の油回収機器類実験用の大型造波回流水槽である．実

番 HBDV-3 ULTRA，図-3.8），表面張力は Data physics 社の接際の重油を散布し，波，流速，風を与えて，海水温を管理し，触角計（型番 OCA15Plus，図-3.9）を用いて計測した．油は視認性をよくするために油用染色剤（シラド化学製 V タイプ

図-3.4 のように油回収に関する機器類の性能の試験を行うこ

Liquid Neutral Red SST-D）を用いて染色したものを使用した．

とができる．

大型水槽の実験は以下の手順で行う．

小型水槽の実験は以下の手順で行う．

14

図- 3.6: 小型水槽での油投入の様子
図- 3.4: 油回収機性能試験の様子
表- 3.2: 油の種類と物性
油種

略称

密度

粘度
（温度 11 ◦ C）

ρo

µo

kg/m3

×10 −3 Pa·s

DN002

826

3.8

VG010

860

28.8

DN460

896

2200

油-大気間

油-水間

正味の

表面張力

表面張力

表面張力 ∗

出光興産
ダフニー 2

Mobil
VG10
出光興産
図- 3.5: 油投入量の測定の様子

ダフニー 460
略称（1）定量の油を内径 23 mm の 2 重円筒状の油投入器に注ぎ，水槽中央に設置する（図-3.10）．
（2）水槽天井から鉛直下向きに設置したカメラ（Nikon D5200）

σoa

σow

σn

×10 −3 N/s

×10 −3 N/s

×10 −3 N/s

で動画撮影をスタートする．

DN002

21.0

13.4

36.7

（3）動画撮影スタートから約 3 秒後に油投入器を水面を乱

VG010

22.8

15.1

33.3

DN460

29.4

さないように気を付けて水中に下げる．それと共に，油投入∗ 水-大気間表面張力 σ

器内の油が時間とともに拡がる．

27.8
=71.1
mN/m
wa

13.9

（4）油膜が水槽壁に到達したら動画撮影をストップする．
大型水槽の場合，油自体は画像解析に適するほどの画質で次に，図-3.13 は油種と拡がり速度の関係を示している．こ撮れない．そのため，40×40×3 mm 程度の発砲ポリスチレれを見ると油の量が等しくても，物性により拡散速度は異なン製のトレーサーを油膜の外縁で押されて動くように 8 個

ることがわかる．

配置して（図-3.11 参照），トレーサーの位置を解析すること

図-3.12 及び図-3.13 の結果を踏まえて，油の拡がりに支配

で油膜の面積を求めた．トレーサーの位置は画像解析ソフト

的な力について考察する．まず油投入量と拡がり速度の関係（DITECT 製，DIPPMotionPro）を用いた．実験は油種及び油がないことから重力（油と水の密度差による浮力）は支配的

量ごとに 1〜3 回行った．

ではない．また，油の粘性のオーダーが 103 倍異なる DN002

3.4

室内実験の結果と DN460 で面積の時間変化が数倍程度であることから，油のまず，大型水槽を用いた実験結果について説明する．図-

粘性も拡がり速度に大きく影響を及ぼさないと考えられる．

3.12 は油の投入量と拡がり速度の関係を示している．これをそのため，本実験のように，油膜厚さが 10 −4 m オーダー以

見ると油量によって拡がり速度は変化せず，ほぼ一様である

下においては油の拡がりについて正味の表面張力が支配的でことがわかる．この時の油投入量と拡散面積より，少なくとあると考えられる．

も油膜厚さ 10 −4 m 以下のオーダーでは油量は拡散速度に大そこで，Fay (1969) にならい，表面張力-粘性領域のモデル

きく影響を及ぼさないことが分かる．

を小型水槽の実験結果にあてはめ，拡がり係数を評価する．

15

図- 3.7: 比重計(a) 装置(b) 液膜作成用シリンジ

図- 3.8: B 型粘度計

図- 3.9: 接触角計

また，大型水槽の実験結果と比較してその汎用性を評価する．
果の平均値と理論解である式 (3.9) の比較である．実験結果

小型水槽で行った実験結果に式 (3.8) が成り立つと仮定しは油量では違いがないと考えられることから，油種で平均して，t を変数として最小二乗法で拡散面積 A の近似を行い，て比較する．この結果を見ると，VG010 の実験結果が図-3.15

拡がり係数 k 3 を求める．なお，画像解析に適する画像が実の結果と同様に，理論解よりも小さい傾向にあった．しかし

験開始直後から撮影できなかったため，式 (3.9) のように，時ながら，DN002 及び DN460 の実験結果は理論解と一致して間方向に近似式を t a 秒シフトさせ， k 3 の値のみ評価する．
− 21

A = k 32 π ρ−1w νw σ n (t + t a ) 2
3いるように考えられる．そのため，Fay のモデルが油の拡がりのモデルとして有効であり，小型水槽で得られた拡がり係(3.9)

数は妥当であると考えられる．

近似の結果の例は図-3.14 の通りであり，本実験結果は式

3.5 油膜自身の特性による油拡散に関する数値計算モデ(3.9) とよく一致するものと考えられる．すなわち，油膜厚さルの開発

が 10 −4 m オーダー以下における油の拡がりは正味の表面張

本節では Fay（1969, 1971）のモデルに基づき，油膜自身の力と経過時間で予測が可能であるといえる．

特性による油拡散に関する数値計算モデルを導出する．

図-3.15 は式 (3.9) で得られた拡がり係数 k 3 の分布を示している．拡がり係数 k 3 にはばらつきがみられるほか，VG010

まず，Fay の第 3 段階である表面張力-粘性領域について，

は係数が比較的小さく，DN002 は比較的大きいといった傾向

ランダムウォークの考え方に基づいた数値計算モデルを導出がみられた．しかしながら，ここでは油の種類が異なってもする．Fay の考え方に基づくと，油の拡がり面積は正味の表

油の拡がり速度を拡がり係数と正味の表面張力で説明するこ面張力に比例し，経過時間の 3/2 乗に比例する．また，一様とを目指すため，結果はほぼ一様であり，油の拡がり係数は乱数を用いて粒子をランダムウォークさせると，粒子の外縁

油種によって変わらないと考え，3 種類の油の実験結果の平半径は経過時間の 1/2 乗に比例して拡がる．そこで式 (3.10)

均値をとって拡がり係数 k 3 =0.852 とする．

のように，経過時間 t の 1/4 乗を拡がり移動速度 Usi にもたせることを考える．

小型水槽で得られた油の拡がりの理論解及びその拡がり係



数の汎用性を評価するため，大型水槽での拡がり実験結果と
1
2

1
4

Usi = k 3m σ n t i ∆t

比較する．図-3.16 は大型水槽で実施した油の拡がり実験結16− 12 cos(2πRr and ) sin(2πRr and )



(3.10)

①

②

③

④

(a) 全体図

図- 3.11: 大型水槽で油が拡がる様子（時間経過と共に 1 から

Spreading area (m2)

4 の順番）(b) 2 重円筒の油投入器

図- 3.10: 大型水槽で用いた油投入器

12
8
4
0
0

10

20
Time (s)ここで k 3m は 3.4 節で説明した実験結果から導出した拡がり係数を用いた Fay のモデルとの比較で求まる拡がりモデル

25ml(1)50ml(1)
100ml(1)係数，t i は粒子 i を投入してからの経過時間を表す．具体的には，油粒子 10,000 個を用いて式 (2.3) において，Usi 以外

3040

25ml(2)50ml(2)図- 3.12: 油の投入量と拡がり速度の関係（油種 DN002）の速度ベクトルをゼロとして拡がり計算を実施し，全粒子の95%が収まる半径を経過時間ごとに求める．求まった拡がり
半径 R から拡がり面積 A を求め，式 (3.9) と一致する k 3m を円形に拡がることがほとんどない．そのため，Fay のモデル

最小二乗法で求める．

を粒子モデルに適応するには工夫が必要となる．例えば Guo

図-3.17 は提案する計算モデルを用いた油の拡がりの計算

and Wang (2009) は以下のように取り扱っている．すなわち，例である．油種は DN002，経過時間は 60 秒である．図より拡がった油の重心からの平均半径が l の時の拡がり速度 Usi

提案するモデルを用いると，ランダムウォークのため油膜外を以下のように計算する．

縁はばらつくけれども，円形に拡がる様子を計算できること

Usi =

がわかる．

r oi dll dt(3.12)

図-3.18 は油の拡がりに関する理論解と数値計算モデルの比ここで，r oi は油膜の重心と着目粒子 i の距離ベクトル，dl/dt

較である．これより，式 (3.10) で油の拡がりを計算するモデは Fay のモデルから求まる半径 l における拡がり速度である．

ルは理論解を表すことができると考える．この時の式 (3.10)

著者はこの考え方の難点は以下のとおりであると考えていの拡がりモデル係数 k 3 m は式 (3.11) の通りで次元をもつ．

る．式 (3.12) は油膜の重心から遠い粒子ほど速い移動速度を持つことを意味している．しかし，これでは粒子の間隔が粗

k 3m = 0.554 [N − 2 m 2 s − 4 ]1

3

3(3.11)

であっても速い移動速度を持つことになる．粒子が密である

次に Fay (1969) のモデルに基づき，近傍粒子から受ける斥

とき，つまり油膜の厚さが厚い時，強い発散を持ち，速い粒
子移動速度を持つはずであると考えられる．

力を計算して粒子の移動を決定する斥力モデルを導出する．
前述した Fay のモデルは油膜が円形に拡がるモデルであるけ

そこで著者は粒子の集合が円形であるかに由らず，粒子のれども，実際の油拡散は乱流拡散，移流の効果もあるため，粗密により粒子の移動速度を計算する方法を導出することと17Spreading area (m2)

Spreading area (m2)

12

12
8
4

8
4
0
0

0
0

10

20
Time (s)DN002(1)
VG010(1)
DN460(1)
DN460(3)

30

10

20

30

40

Time (s)40

DN002(2)
VG010(2)
DN460(2)Experiment (DN002)

Theory (DN002)

Experiment (VG010)

Theory (VG010)

Experiment (DN460)
Theory (DN460)
図- 3.16: 大型水槽での実験結果の平均値と理論解の比較

図- 3.13: 油種と拡がり速度の関係（油量 50 ml）

3
2

0.3

1

0.25

y (m)

Spreading area (m2)

0.35

0.2
0.15

0
-1

0.1
-2

0.05
0
0

-3

0.1 0.2 0.3 0.4 0.5 0.6 0.7

-3

-2

-1

Time (s)Experiment

0

1

2

3

x (m)

Fay's model (1969)

図- 3.14: 拡がり面積の実験結果と近似式の比較の例（油種図- 3.17: モデルを用いた油の拡がりの計算例（油種 DN002，

DN002）

60 秒後の再現）

Spreading Area (m2)

Spreading coefficientk3 (-)

35

1.2
0.9
0.6
0.3
0.0
0

10

20

30

25
20
15
10
5
0

40

0

Net surface tension σn (10-3N/m)
DN002
DN460

30

20

40

60

Time (s)VG010
Average value

図- 3.15: 実験結果に基づく拡がり係数の分布

Theory(DN002)

Simulation(DN002)

Theory(VG010)

Simulation(VG010)

Theory(DN460)

Simulation(DN460)

図- 3.18: 油の拡がりに関する理論解と数値計算モデルの比較18

yu si j ∝ r i−1
j

Nの範囲(3.18)である．すなわち，油拡散初期において粒子の移動速度は粒子間距離の-1 乗に比例する．

i

usij

rij j

Fay のモデルの重力と粘性力が卓越する第 2 段階も同様に
式 (3.7) から粒子の移動速度を求めると，粒子の移動速度は
粒子間距離の-3 乗に比例することが導出される．
数値計算の際は粒子の移動速度を以下のように求める．粒
子 i の拡がり速度 Usi は周囲の粒子 j から受ける速度 u si j をx

ベクトル合成して計算する．

図- 3.19: 粒子の発散の概念図
した．

Usi =

モデルは連続の式より導出する．ある着目粒子周りの粒子

u si j =

∑

k 1m

j

密度 no は連続であることから，油を粒子の集まりと考え，

Usi =ラグランジュ的に連続式で表わすと，

Dno+ n o ∇ · Us = 0
Dt

N
∑

N
∑
j

u si j =

∑

k 2m

ri j
|r i j | 2
ri j|r i j | 4(3.19)(3.20)

ここで，k 1m ，k 2m は Fay のモデルとの比較で求まる拡がり(3.13)

係数（単位はそれぞれ m2 /s，m4 /s）で，計算方法は後述する．

となる．ここで Us は着目する油膜位置の油拡散速度を表す．

式 (3.19) は Fay の第 1 段階，式 (3.20) は第 2 段階を表わす．

式 (3.13) の左辺をそれぞれ離散化する．まず式 (3.13) 左辺第導出した油膜自身の特性による油拡散の計算モデルに関す

2 項について，発散の概念を表したのが図-3.19 である．

る検証計算を行った．油の種類は原油アラビアンライトとし，粒子の移動速度をラグランジュ的に考えると，着目粒子 i

計算は粒子数 1000 個，粒子体積 10 m3 /個，k 1m = 1.2 × 10 −2から r i j 離れた粒子 j の速度は u si j ，着目粒子の相対速度は

m2 /s， k 2m = 2.0 × 103 m4 /s で行い，着目粒子周りの油層厚ゼロとおくと，式 (3.13) 左辺第 2 項は，

h1 が 3.2 mm 以上のとき式 (3.19)，それ以外のとき式 (3.20)

n o ∇ · Us = n ou si jri jを用いた．

(3.14)

図-3.20 は原油 1×104 m3 が静水中に瞬間流出した際の粒となる．次に，式 (3.13) 左辺第 1 項の粒子密度の時間変化

子モデルの計算結果で，6，12，24，72 時間後の拡散を表わを求める．粒子密度が高いほど時間変化が大きいと考え，式している．油粒子は流出点より円形に拡がっており本計算に(3.15) を仮定する．

おいて油が拡がる様子を再現できていると考える．
図-3.21 は本モデルを用いた計算結果と Fay のモデルの比Dno= −CnαoDt較である．計算は瞬間流出（図-3.20 の結果）と，連続流出を(3.15)

対象とした．連続流出は 1×104 m3 の原油を流出総時間に応ここで Cno は係数，α は指数である．式 (3.15) を積分すると，

1n −α+1 = −Ct + K
1−α o

じて一定流量で流出させて計算した．流出時間はそれぞれ 6
時間，24 時間，72 時間とした．解析に用いる拡散半径は全(3.16)

粒子の 95%が含まれる半径と定義した．瞬間流出の計算結果である．ここで，K は積分定数である．油層の体積が一定のは Fay のモデルとよく一致しており，図-3.20，3.21 より油拡

とき，Fay のモデルの重力と慣性力が卓越する第 1 段階にお散半径を本モデルで再現できていると言える．また，瞬間流

いて，式 (3.6) から式 (3.17) が導かれる．

出と連続流出を比較すると，連続流出の場合は瞬間流出に比較して初期の油量が少ないため拡散半径が小さく，時間が経() −1/2∆ρh ∝ l −2 ∝
gV
t −1
ρw

過し油量が多くなるにつれて瞬間流出のときの拡散半径と近(3.17)

くなっていく．本モデルにおける油膜自身の特性による油拡

油を粒子の集まりと考えており，no は h に比例すること

散速度は，水面に存在する油の集まり具合に大きく影響を受から，no は t −1 に比例する．この関係を式 (3.16) に用いる

けるため，このような計算結果になると考える．

と，α =2 を得る．よって，式 (3.13) に式 (3.14) 及び式 (3.15)

図-3.22 は x 方向に 0.1 m/s の流れがある状況で油が流量を代入すると，l と r i j は比例することから，

2×102 m3 /hour で流出した際の油の移流拡散を再現している．

19

Spreading radius (m)

2500
2000
1500
1000
500

y (m)

0
0

86400
172800259200Elapsed time (s)Instant spill
Fay equation(Instant spill)
Continuous spill
Continuous spill
6hours24hours
Continuous spill
72hours

図- 3.21: 斥力モデルによる油粒子の拡散半径と Fay のモデル

x (m)

の比較

図- 3.20: 斥力モデルによる油粒子の拡散計算結果（6，12，24，

72 時間後）がりながら下流方向へ移動していく様子がわかる．
次に，計算結果について説明する．図-3.25 は定常流におFay のモデルを直接用いた式 (3.12) の計算では表現が困難でいて油を連続流出させたときの数値計算結果である．流れのあった連続流出を計算できている．

中で移流しながら拡散する様子がよく再現されていると考え

3.6

定常流における油拡散実験と数値計算モデルの比較られる．

本節では定常流条件下における油の拡がり実験と導出した

実験結果と計算結果の比較を行う．図-3.26 は実験結果と油の拡がりの計算モデルを比較し，数値計算法の評価を行う．

計算結果の比較である．実験結果は x 方向を 0.2 m で区切

実験は 3.3 節で説明した静水時の油膜自身の特性による油拡り，その範囲のトレーサーの y 方向長さの平均値をプロット

散実験を行ったのと同じ大型水槽で実施した．図-3.23 は実した．計算結果は油流出から 40 秒後の結果を用いて，x 方向

験の模式図である．

を 0.2 m で区切り，その範囲の油粒子の 95%が含まれる y 方

定常流において油の連続流出拡散実験は以下の手順で行

向長さをプロットした．実験と計算の標準偏差は 0.26 m でう．

ある．これをみると，計算結果は実験結果とよく一致してい(1) 初期画像データを読み込む．

ると考えられる．よってこれらの結果から，汎用性の高い計(2) 流速 0.135 m/s の定常流を発生させる．

算モデルを用いて，計算結果が得られたと考える．

(3) 水槽の上流側横断方向中央部分に油投入口を設け，モン

3.7 考察ローポンプで流量 1.0 L/min で油 DN002 を連続で投入する．

本節では油膜自身の特性による油拡散に関する，実験結果(4) 水槽天井から鉛直下向きに設置したカメラ（Nikon D5200）

及び数値計算モデルの開発を踏まえて考察する．まず，実験で動画撮影をスタートする．

で得られた表面張力-粘性領域の油の拡がり係数 k 3 について，(5) 油膜の外縁を解析するためのトレーサーを上流側中央部

表-3.3 は既往の研究と本実験で得られた拡がり係数を比較しから流す．

ている．これらの結果を比較すると，拡がり係数は本実験結(6)1 分程度撮影し，動画撮影をストップする．

果と大きく異なることがわかる．埜口らの実験は， 油層先

計算は 3.5 節で説明したランダムウォークモデルを用いて端部の形状の観察が困難で，詳細な運動の機構は検討でき

速度ベクトルを求める．式 (2.3) の Uoi に水槽の定常流速でなかったとあるため，油膜の外縁を観測できた本実験の拡がある 0.135 m/s を与える．油粒子は 1 秒間に 200 個の割合でり係数のほうが埜口らの結果より大きい値となることは妥当

連続的に投入した．

であると考えられる．

実験及び計算の結果について説明する．まず，実験結果に一方，本実験においても VG010 を用いた結果については，

ついて説明する．図-3.24 は定常流において油を連続流出さ

精度よく再現されているわけではない．また，C 重油を用いせたときの油膜外縁を移動するトレーサーの軌跡の重ね合わた同様の実験を行ったけれども，提案する理論解では拡がりせ画像である．油が流出し，油膜自身の特性によって油が拡をうまく説明できなかった．そのため，本実験で得られた拡

20

カメラ
定量
ポンプ

断面図

2.5m

油
油投入口

平面図

20m

定常流
定量
ポンプ

油

y (m)

6m
図- 3.23: 定常流における油拡散実験の模式図

3
流れ方向

2

y (m)

1
0
-1
-2

x (m)

図- 3.22: 斥力モデルによる流れ場における油粒子の移流拡散-3（6，12，24 時間後）

-0.5

0.5

1.5

2.5

3.5

4.5

x (m)
がり係数 k 3 の取り扱いは慎重に行う必要があり，今後実験図- 3.24: 定常流における油拡散に関する実験結果

回数を増やす等によりそれらの実験結果を考察する必要がある．

は，精製水に SC1000 という分散剤をまぜると，水-大気の表

また，Fay（1969, 1971）は油流出からの経過時間で 3 つの面張力係数は 72 mN/m から 29 mN/m に下がった．すると，拡散過程があると示した．本実験で得られた結果は表面張

油-大気及び油-水の表面張力の合力が水-大気の表面張力より力が卓越した拡がりに関するものであると考えられる．よっも大きくなり，正味の表面張力係数は負の値をとるため，表て，Fay の示す重力が卓越した拡がりについて検討するには，

面張力は油の拡散面積を縮める方向に作用した．界面活性剤

投入する油量を多くした実験を行う必要がある．特に，Fay

以外にも何らかの影響で油流出海域の海水が変質をうけてのモデルには油の粘性の効果が含まれていないけれども，重いると，水-大気の表面張力係数が下がることがある．その力が卓越する状況における油拡散速度は油の粘性の影響が大ため，油膜自身の特性による油拡散は，理論的には計算が可きいと考えられることから，実験により確かめる必要がある

能でも，実際の油流出の状況を考えると，精緻に予測すること考えられる．

とは困難であることを念頭に数値計算を実施することが望ま

ここで，実際の油流出が発生した際の注意点について述べれる．

ておく．表面張力が卓越する際の油拡散を予測する式 (3.9)

次に，開発した油膜自身の特性による油拡散に関する数値において，正味の表面張力係数が拡がり速度を決定する．正計算モデルについて考察する．ランダムウォークモデルを用

味の表面張力は水-大気，油-大気，油-水の表面張力から決定

いた数値計算モデルは，粒子間に斥力を与える数値計算モデされるけれども，水-大気の表面張力は油分散剤のような界ルと比較して，モデルが簡易であるため実装が容易であると面活性剤によって小さくなる．例えば，著者が行った試験でいう利点がある．油回収船で油回収作業を行う際に，流出油

21

3

3
流れ方向

2

2
y (m)

y (m)

1
0

1

-1

0

-2

0
0

1

2

3

4

2

x (m)
Experiment

-3
-1

1

5

3

4
Model

図- 3.26: 定常流における油拡散に関する実験と数値計算結果

x (m)

の比較

(a) 油流出から 20 秒後

3

表- 3.3: 既往の研究と本実験で得られた拡がり係数 k 3

y (m)

2
1
0
-1

拡がり係数 k 3

実験/理論

Fay (1971)

2.30

理論

Hoult (1972)

0.128

理論

埜口他 (1979)

0.48

実験

本実験

0.852

実験

-2
-3
-1

0

1

2

3

4

場合には，ランダムウォークモデルと斥力モデルを組み合わ5

せて使い分ける．Fay の第 1 段階の式 (3.19)，第 2 段階の式

x (m)(3.20) 及び第 3 段階の式 (3.10) の使い分けは以下のように判(b) 油流出から 40 秒後

図- 3.25: 定常流における油拡散に関する数値計算結果

断する．まず，第 1 段階の式 (3.19) と第 2 段階の式 (3.20) においては，以下のように計算する．各計算ステップ，着目粒の漂流予測結果は即時に出力，表示できる必要があり，船舶子ごとに，粒子 i から近傍粒子検索半径 R 以内にある粒子のに搭載するようなパソコンは特に高速な CPU を持たないも個数 No を数え，粒子体積 Vi から油層厚 h を，のを想定しているため，計算負荷が小さいランダムウォークh=

モデルは有利である．ランダムウォークモデルの問題点は油が拡がる際の油層厚分布がランダムウォークを用いているた

No Vi
πR2

(3.21)と計算し，Fay のモデルから求まる油層厚判定値 h1 と h をめ，正規分布を仮定しているけれども，実際の油拡散は正規

比較して決定する． h が h1 より大きいとき，Fay の第 1 段階分布になるか確認されていないことが挙げられる．

である式 (3.19) を用い，h が h1 より小さいとき，Fay の第 2

一方，粒子間に斥力を与えて計算する数値計算モデルは，段階である式 (3.20) を用いる． k 1m ，k 2m ，h1 は想定してい油が大量に流出した際の重力による拡がりをよく再現でき

る油の密度と粒子の体積による． k 1m ，k 2m の決定は以下のると考えられる．しかしながら，着目粒子 i の近傍粒子 j そように行う．まず，油の拡がりの第 1 段階を表す式 (3.6) とれぞれから受ける移動速度 u si j をベクトル合成して計算す

第 2 段階を表わす式 (3.7) が等しくなるときの l を求め，l かるため，近傍粒子数の分だけ計算時間がかかることが難点で

ら h1 を求める．次に h1 を用いて k 1m を変えて式 (3.19) のある．

計算を行い，結果と式 (3.6) を比較し，誤差が最も小さい時以上の数値計算モデルの特徴を踏まえて，ランダムウォーの k 1m を求める．最後に h1 を用いて k 2m を変えて式 (3.20)

クモデルと斥力モデルの使い分けは以下のように行う．まの計算を行い，結果と式 (3.7) を比較し，誤差が最も小さい

ず，即時に結果を得たい場合や，特に高速な CPU を持たな

時の k 2m を求める．

い PC で計算を実施する際はランダムウォークモデルのみを次に，第 2 段階の式 (3.20) と第 3 段階の式 (3.10) において用いて，油膜自身の特性による油拡散を計算する．

も同様に，油層圧 h を計算し，油層厚判定値 h2 と比較して次に，計算に時間がかかってもよい場合で，大量の油流出決定する． h が h2 より大きいとき，Fay の第 2 段階である式が発生し，油の拡がりに重力領域が重要であると考えられる(3.20) を用い，h が h2 より小さいとき，Fay の第 3 段階であ22

る式 (3.10) を用いる． h2 は油膜の厚さで決まるとして，3.4

3.8 結語節の結果を考慮し h2 =1.0×10 −4 m とする．

本章では油膜自身の特性による油拡散について，室内実験

近傍粒子検索半径 R の与え方に関して，特に理論的な検による油の拡がりの評価，数値計算モデルの導出，実験結果

討は行っていないけれども，試行錯誤的に検証を行った結果と数値計算モデルの比較を行った．以下では本章で得られたによると，Vi =1 から 10 m3 のとき， R=200 m 程度で計算を主要な結論を要約する．

行うとよい．

まず，油膜自身の特性による油拡散について有名なモデル

次に油膜の拡がりの終端について考察する．時間経過と共である Fay のモデルについて紹介した．Fay は油膜自身の特に油膜は拡がり，最終的に油膜の拡がりは停止する．正味の性による油拡散は重力及び表面張力によって拡がり，慣性力

表面張力は水-大気間の表面張力が油-大気間及び油-水間の表及び粘性力によって拡がりが抑えられるとし，油が流出し拡

面張力の合力より大きい場合に油膜が拡がる方向に作用すがる過程で支配的な因子が時間経過に伴い，第 1 段階；重力る．そのため，油膜の拡がりが停止するのは油膜が単分子層‐慣性領域，第 2 段階；重力-粘性領域，第 3 段階；表面張まで拡がった時と考えられる．例えば，本章で説明した小型力-粘性領域と分けられるとして，拡がりに関するモデルを水槽実験では，油膜は厚さ 10 −6 から 10 −7 m のオーダーに導出したことを説明した．

おいても拡がっている．

次に，小型水槽及び大型水槽を用いて油膜自身の特性によ

しかしながら，前述の通り，正味の表面張力は実海域におる油拡散実験を行った．その結果，実験結果の拡がり速度はいて変化する場合がある他，Hoult (1972) は，以下のことをFay のモデルの第 3 段階；表面張力-粘性領域とよく一致し指摘している．

た．拡がり係数は Fay の提案する値よりも低い値となり，埜

すなわち，原油のように様々な構成要素からなる油は，そ

口らの結果より大きい値となった．拡がり係数は油の拡がれぞれの構成要素ごとに表面張力係数をもつため，表面張りの予測に大きく影響を与えるため，その決定は慎重に行わ

力の値を正確に見積もることが困難である．また，海水の温れる必要がある．また，実験時の油膜厚さから，油膜厚さが度，塩分濃度，生物学的物性による違いが正味の表面張力の10 −4 m オーダーの場合，油の拡がり速度は正味の表面張力

大きさを変化させる要因となる．そこで，実務面の有益性かによって決定されることがわかった．

ら 2 つの考え方が提案されている．1 つは Blokker (1964) 及次に，開発した数値計算モデルについて説明した．まず，び Allen and Estes (1970) に提案される，油膜が一定の厚さ以ランダムウォークを用いた油膜自身の特性による油拡散に関下になると停止するという考え方である．この時の油膜厚さする数値計算モデルについて説明した．計算対象は第 3 段の値は経験的に 2.5×10 −5 m とされている．もう 1 つは Fay

階；表面張力-粘性領域とし，油の拡散面積が正味の表面張の提案する油流出時から時間経過に伴い界面活性物質が水中力の 1 乗，経過時間の 3/2 乗に比例する数値計算モデルを導に溶解し，表面張力が低下して，拡がりが停止するという考出した．導出したモデルは静水時の油膜自身の特性による油え方であり，最終拡がり面積は式 (3.22) のように求められる．

拡散実験結果から得られた拡がり係数を用いた理論解とよく一致した．これは油の種類が異なっても統一のモデルを用い18
σ2 V 6A ≈  2 n 3 6 
ρ w νw D o s dて計算ができるため，有効なモデルであると考える．

(3.22)

次に，斥力の考え方を用いた油膜自身の特性による油拡散

ここで，Do は界面活性物質の水中での拡散係数，s d は界面の数値計算モデルの導出を行った．計算結果は Fay のモデル

活性物質の水中への溶解度である．なお，Fay の提案するモの第 1 段階；重力‐慣性領域，第 2 段階；重力-粘性領域とデルについて，実際に水中に界面活性物質が溶解することを一致するように拡がり係数を求めた．導出した数値計算モデ

確認しておらず，モデルの物理的意味は不明瞭であると考え

ルは油拡散について計算できることを示した．また，静水時られる．

のほか，流れがある場合の油の連続流出に関する数値計算を自然な表現で行えるようになった．

海上災害防止センター (2000) によると，水面が銀色にキラキラ光って見える状態のとき，油膜厚さは 1.0×10 −7 m であ次に，大型水槽を用いた定常流条件下の油拡散実験結果とり，光線の条件が最も良い時にかろうじてキラキラ光る油膜

導出したランダムウォークモデルの計算結果を比較した．計が見える状態のとき，油膜厚さは 5.0×10 −8 m である．

算結果は実験結果をよく再現していることが示された．
最後に，2 つの数値計算モデルの使い分けと，油の拡がり

そこで，本数値計算において，油膜の拡がりの終了についての終端について考察した．

は以下のように考える．すなわち，粒子密度から式 (3.21) で求める油層厚が 1.0×10 −7 m よりも薄くなったとき，式 (3.10)の計算を停止する．

234.

海水面極近傍における水平乱流拡散係数の推定

おいて選定する．次に，実海域において海水面極近傍におけ

方法

る水平乱流拡散の測定実験を行う．また，実験結果から，水4.1

平乱流拡散係数 D H の推定方法を考察する．

概要

4.2 風洞水槽実験による疑似油の選定

海上に流出した油の漂流予測計算において，流出した油を
粒子の集まりとしてラグランジュ的に計算する場合，海表面

実海域で海水面極近傍に限った油拡散特有の拡散係数を計の乱れによる油の乱流拡散はランダムウォークで再現するこ

測するにあたり，対象となる疑似油は海水面のみを漂い，油とがほとんどである．すなわち，ある粒子 i の流れの乱れにと同じ漂流速度でなければならないと考えられる．流れによよる水平拡散速度 Udi は 2.4 節で述べたとおり，式 (2.18) でる移流速度は疑似油と実際の油とで変わりがないと考えられ計算する．

るため，油と同じ風係数を持つ疑似油を選定することが必要であると考えた．

式 (2.17) の通りランダムウォークの計算は水平乱流拡散係
数 D H を与える必要があり，従来は定数を与えていた．乱2.2 節で述べた通り，数値計算において風による油の移動

流拡散速度 Udi は D H の 1/2 乗に比例するため，拡散面積は

速度は海面高度 10 m での風速に風係数をかけて表現する．

D H の 1 乗に比例する．そのため，水平乱流拡散係数が計算

風係数は 0.03 程度であり，例えば海面高度 10 m で風速が 10

結果に与える影響は大きいといえる．

m/s である場合，海表面の油の移動速度は 10×0.03=0.30 m/sと算出される．今回行う実海域拡散実験で用いる疑似油の風

しかしながら，従来では経験的に水平乱流拡散係数 D H を

係数も 0.03 である必要がある．

決定しており，具体的にどのように与えればよいかわかっていなかった．例えば ASCE Task Committee on Modeling of Oil

そこで，疑似油の選定実験を国土交通省国土技術政策総

Spills of the Water Resources Engineering Division (1996) におい合研究所の所有する台風防災実験水路（図-4.1 参照）で行って乱流拡散係数は 1 から 100 m2 /s としている．これでは設定

た．図-4.2 は実験水槽の全体図である．水路は長さ 28.5 m，

値により拡散面積計算結果が 100 倍異なることを意味する．

幅 1.5 m，高さ 1.3 m である．実験は水深を 0.4 m として行っ

また，乱流モデル（例えば，Smagorinsky (1963)）を用いて数た．疑似油は素材としてスポンジゴム（イアノック EPDM 系

値計算で拡散係数を決定する方法も考えられるけれども，従E-4388，硬度 20±5 度，密度 150 kg/m3 ）を用いることとした．

来その適応はよく検討されていなかった．

実海域実験では直径 1 m の円形を用いるものの，その大き一方で，染料や漂流物を用いた実海域における拡散の観測

さの疑似油は水槽実験で扱えないため，直径 30 cm，厚さ 10が多く行われており，水平乱流拡散係数の推定方法が提案さ

mm の円形の模型を用いた．また，実海域実験においては，れている．例えば，Okubo (1971) は染料を用いた拡散実験を

疑似油に小型の GPS ロガーを搭載するため，防水のケース行い，拡散スケールと水平乱流拡散係数の関係について整理にいれてマットにくくる必要がある．水槽実験では全体重量している．

（ゴムマットと防水ケースを合わせた比重）を合わせるため

柳・樋口 (1982) は海面下 1m に抵抗板をつりさげた 60cmに防水ケースにいれる GPS ロガーの模型重量を調整した．

四方の耐水ベニヤ板を浮子とし，その位置を航空機で観測す実験は以下のように行った．図-4.3 は実験の模式図である．
送風速度は 3 通り（送風機回転数：100 rpm，200 rpm，300

る方法で拡散実験を行っている．
道田他 (2009) は海面からそれぞれ 1 m，6 m，11 m，16 mrpm）で実施した．送風後，水槽内の流れが定常状態になっての位置にドローグ（抵抗体）をとりつけた GPS 搭載漂流ブから，疑似油を十分な加速区間をとった風上側から投入し，イの放流実験を行っている．

送風口から 19.5 m の位置にあるカメラ（DITECT 製）で疑似

これらの実験は，いずれも水面から数 m における水平乱油の移動を撮影した（図-4.4 参照）．撮影画像は画像解析ソ

流拡散係数を推定していると考えられる．しかし，油拡散のフト（DITECT 製，DIPPMotionPro）を用いて，疑似油の中心ような海水面極近傍に限った水平乱流拡散係数の推定は従来

位置の座標を測定し，移動速度を求めた．投入回数は各ケー行われていなかった．道田他 (2009) は水平乱流拡散の水深

スで 30 回ずつ行い，その平均値を解析に用いた．

依存性を指摘している．また，海面付近のほうが水中よりもまた，カメラ撮影位置を挟む送風口から 18.0 m 及び 21.0

乱れが大きく，拡散係数が海水中より大きくなるのではない

m 地点において風速（熱式風速計：カノマックス製プローブ

か，と予想される．

0962-00）を計測した．漂流予測を行う際は海上 10 m 風速が以上の問題を解決するため，本章では実海域において海水用いられるため，今回の実験結果も水面から 10 m に換算し面極近傍の乱流拡散を測定することを試みる．そのためにまた風速で整理する必要がある．そのため，海上 10 m 風速をず，海上流出油と同等の漂流特性を持つ疑似油を風洞水槽に

以下の方法で算出し，2 地点の平均値を解析に用いる風速とした．すなわち，風速計が水面（波面）に接しない高さから，

24

送風口

21.0m
19.5m
18.0m

カメラ撮影
風速
測定
位置

風速
測定
位置

風
疑似油移動

図- 4.3: 実験の模式図
図- 4.1: 台風防災実験水路図- 4.4: 画像計測方法照）．GPS ロガーは Wintec 製 WBT-202 で，重量 55 g，平均誤差半径（CEP）は SBAS 使用時で 2.0 m，Stand Alone で 2.5m である．各実験ケースにおいて疑似油を 20 枚程度同時刻に投入し漂流させ，各疑似油の位置を GPS ロガーで記録し図- 4.2: 実験水槽全体図

た（図-4.9 参照）．同時に図-4.10 のように，GPS コンパス，
風向風速計，ADCP を載せた観測船で疑似油の数十 m 付近を

静水面より 10 cm までを 5 mm 毎の高度で各 60 秒間の計測

ゆっくり追跡しながら風向風速，流向流速を測定した．各実を行い，それぞれの高度の平均風速を求めた．平均風速の鉛

験ケースの漂流時間は三浦市周辺の実験においては，30 分

直分布に対数則が成り立つものと仮定して，最小二乗法によから 3 時間 40 分で，合計 9 ケース行い，平塚市周辺の実験

りフィッティングを行い，海上 10 m 風速を算出した．

においては，2 時間から 5 時間で，合計 7 ケース行った．

疑似油選定実験の結果は表-4.1 及び図-4.5 の通りである．

4.4 実海域における疑似油を用いた水平乱流拡散実験の次節以降で述べる実海域実験時の平均風速は最大でも 10 m/s

結果以下であるため，今回の室内実験の風速は実海域実験時の風速が含まれている．風係数は風速が大きくなると小さくなる

GPS ロガーで取得した疑似油約 20 枚の漂流緯度経度から，傾向がみられるが，ここでは風速 10 m/s よりも小さい時をRichardson (1926) の方法に基づき，式 (4.1) を用いて水平乱対象とするため，風係数は 0.03 に近いと考え，この疑似油

流拡散係数 D H を 10 分おきに計算した．

を用いて実海域実験を行うこととした．

4.3

表- 4.1: 10 m 換算風速と風係数の関係

実海域における疑似油を用いた水平乱流拡散実験の送風機

海上 10 m

疑似油

実験は図-4.6 に示す海域で 1 回ずつ行った．1 つは 2012 年

回転数

換算風速

移動速度

10 月及び 11 月に東京湾入口に位置する神奈川県三浦市沿岸

(rpm)

(m/s)

(m/s)

周辺で，もう 1 つは 2013 年 12 月及び 2014 年 1 月に神奈川

100

3.8

県平塚市沿岸周辺で行った．図-4.7 は実験の模式図である．

200

疑似油は直径 1 m の円形，厚さ 10 mm，重量 1.2 kg のスポン

300

方法

ジゴムに GPS ロガーを取り付けたものを用いた（図-4.8 参

25

風係数の
風係数

標準偏差

0.12

0.032

0.0011

6.3

0.18

0.030

0.0022

9.4

0.25

0.026

0.00211.0

Distribution

0.8
0.6
0.4
0.2
0.0
0.020

0.025

0.030

0.035

Wind coefficient
100rpm200rpm
図- 4.5: 疑似油の風係数の分布

300rpm
図- 4.8: 疑似油の例

平塚市
周辺海域
三浦市
周辺海域

図- 4.9: 疑似油漂流実験状況（平塚市周辺）
図- 4.6: 疑似油を用いた水平乱流拡散実験を実施した海域

風速計
GPSコンパス

風向風速計
GPSコン
パス

警戒船

計測用PC
ADCP

ADCP
観測船
CTDセンサー

疑似油20枚図- 4.10: 風況及び流況データ計測方法

風・潮流

図- 4.7: 疑似油を用いた水平乱流拡散実験の模式図

26

Horizontal turbulent diffusion
coefficient (m2/s)

Horizontal turbulent diffusion
coefficient (m2/s)

10
1
0.1
0.01

0.001
0

3

6

9

12

10
1
0.1
0.01
0.001
0.00

Wind velocity (m/s)
Miura
Hiratsuka
図- 4.11: 風速と水平乱流拡散係数の関係

1 ∂σ 2

(4.1)

4 ∂t

ここで，σ 2 は疑似油の位置の不偏分散，t は時間である．σ 2は式 (4.2) を用いて計算した．

σ2 =

1
N −1

N (
∑

(x in − x) 2 + (yin − y) 2

0.20

0.30

Miura
Hiratsuka
図- 4.12: 流速と水平乱流拡散係数の関係
Horizontal turbulent
diffusion coefficient (m2/s)

DH =

0.10

Current velocity (m/s)

)

(4.2)

10
1
0.1
0.01
0.001
0.0

0.1
0.2
0.3
0.4
0.5
0.6
Pseudo oil drift velocity (m/s)
Miura
Hiratsuka
図- 4.13: 疑似油の平均移動速度と水平乱流拡散係数の関係

in=1

ここで，N は計算に用いる疑似油の個数，in は着目する疑似油，x in ，yin は疑似油 in の位置，x ，y は疑似油の位置の平均を示す．
次に結果について説明する．まず，風速，流速，移流速度

田他 (2009) の結果もあわせて記載した．

と水平乱流拡散係数を比較した．図-4.11 は風速と水平乱流

本実験結果から，拡散スケールが大きくなるにつれて水拡散係数の関係を示している．図-4.12 は流速と水平乱流拡

平乱流拡散係数は大きくなる傾向がみられた．この結果は，散係数の関係を示している．また，図-4.13 は疑似油の平均移Richardson (1926) の理論と傾向が一致していると考えられる

流速度（流速 + 風速 ×0.03）と水平乱流拡散係数の関係を示

けれども，傾きは拡散スケールの 4/3 乗よりも大きい傾向にしている．これをみると，いずれのケースにおいても顕著なある．

相関は見られない．これらの結果から，今回のような気象海

そこで，実験結果に水平乱流拡散係数が拡散スケールのべ象条件及び解析条件においては，海水面の流れの速度と水平き乗に比例すると仮定して，縦軸に対数をとり最小二乗法で乱流拡散係数の相関は見られないといえる．極端な条件を考近似し，拡散スケールと水平乱流拡散係数の関係を導出した．

えると，流れと風がない条件において乱流拡散は起こらない

導出結果が図-4.14 の太線である．海水面極近傍の水平拡散と考えることが自然であり，実験結果のように風速及び流速を表す本実験結果と，水中の水平拡散を表す Okubo (1971)，

ともに拡散に影響を及ぼさないのは不自然であると考える．

柳・樋口 (1982) の結果と比較すると，本実験の水平乱流拡

よって，海水面の流れの速度と水平乱流拡散係数の相関がみ

散係数は大きい傾向にあった．道田他 (2009) と比較すると，られない理由として，風や流れに対する相関よりも計測誤差

傾きが道田らの結果が大きいため，拡散スケールが大きくな

範囲のほうが大きく，風や流れの大きさの影響は隠れてしまると水平乱流拡散係数の大小が逆転する．道田らの傾きが正う程度の違いであったと考える．以下では，水平乱流拡散係

しく，拡散スケールの 1.94 乗に比例すると仮定すると，係数は拡散スケールのみに依存すると仮定し，海水面極近傍の数は 7.20×10 −6 となり，道田らの結果よりも大きいことがわ

水平乱流拡散係数の推定方法を検討する．

かる．よって，海水面極近傍では水中に比べて流れの乱れが

4.5

大きく，水平乱流拡散係数も大きく取る必要があると考えら

水平乱流拡散係数の推定方法

れる．

図-4.14 は三浦市周辺及び平塚市周辺の疑似油漂流実験結

また，図-4.14 を見ると三浦市周辺は東京湾の湾口，平塚

果における拡散スケールと水平乱流拡散係数の関係を示し

市周辺は相模湾で外洋であり海域が異なっているけれども，ている．疑似油の拡散スケール L は Okubo (1971) にならい

拡散スケールと水平乱流拡散係数の関係が同様であるとの結

L=3σ と定義した．また，Okubo (1971)，柳・樋口 (1982)，道27

Horizontal turbulent
diffusion coefficient DH(m2/s)

100
10
1
0.1
0.01
0.001
0.0001
10

100

1000

Diffusion scale L=3σ(m)図- 4.15: 時間と空間スケールに対する海の乱流の模式的分類（Toba et al. (1984)）

Experiment (Miura)
Experiment (Hiratsuka)
Experimental Formula離であり，擬 3 次等方性領域の上限までと考えられる．Toba

DH=1.26×10-4L1.42

et al. (1984) によると図-4.15 のように 10 km が擬 3 次等方性

=2.06×10-4L1.15

領域の上限である．したがって，10 km を上限として，流体

Okubo(1971)

DH

Yanagi et al. (1982)

DH=1.23×10-4L1.22

場の再現スケールに応じて小さくすればよいと考える．

Michida et al. (2009)

DH

時の予報体制においては，風データのメッシュサイズは 5 km

8.3 節で説明する，港湾空港技術研究所で構築している即=3.79×10-6L1.94

図- 4.14: 拡散スケールと水平乱流拡散係数の関係である．よって，L ma x は 5 km から 10 km で設定すればよいと考えられる．ここでは L ma x =6 km とすることを提案する．

果が得られた．これより本実験結果を基に水平乱流拡散係数の推定を行った結果を他の海域での油流出に関する数値計算

4.6 結語に用いても適応性が高いと考える．

本章では海表面の流れの乱れによる油の乱流拡散をランダ

ここで，導出した拡散スケールと水平乱流拡散係数の関係

ムウォークで計算する際に重要となる水平乱流拡散係数につを数値計算へ導入する際の方法と適応性を考察する．流出油

いて，従来は具体的な数値が不明であったことを説明した．

の漂流予測計算におけるランダムウォークによる乱流拡散のまた，海水面極近傍の乱流拡散を計測した実験は従来行われ計算は，小さいスケールの乱れが再現できていない部分を補ておらず，水中より拡散係数が大きくなることが予想される完するために行っている．よって考慮する拡散係数は再現すことを説明した．そこで，実海域実験を行い油の拡散スケー

る流体場の格子サイズに依存すると考えられる．

ルから水平乱流拡散係数を導出する方法について検討した．
以下では，本章で得られた主要な結論を要約する．

また，乱流拡散係数が拡散スケールの 4/3 乗に比例する乱流構造は，大きな乱れが順次小さな乱れにカスケードダウン

まず実海域において乱流拡散計測実験を行うにあたり，風し，その乱れエネルギースペクトルが波数 k の-5/3 乗に比例

洞水槽実験により油の代わりとなる疑似油の選定を行い，スする慣性中間領域（擬 3 次等方性領域）に限られると考えら

ポンジゴムシートが油と同等の漂流特性を持っており，疑似れる（Toba et al. (1984)）．

油として使用できることを説明した．

そこで，数値計算において水平乱流拡散係数 D H は，擬 3次に実海域実験の方法について説明した．実海域実験では

次等方性領域の上限までにおいて，流体場の再現スケールで疑似油を 20 枚程度投入し，それらの拡がりを小型 GPS ロ

決定される距離 L ma x まで増加し，それ以上の時は一定であガーで計測した．疑似油の拡がりから求めた海水面極近傍のると考え，式 (4.3) で計算する．

水平乱流拡散係数と風速，流速との間に特別な相関関係は見られなかった．また，海水面極近傍の水平乱流拡散と係数拡


−4
1.42


 1.26 × 10 (l i )DH = 

 1.26 × 10 −4 (L ma x ) 1.42

l i < L ma x
l i ≥ L ma x

散スケールとの相関があることを説明した．既存の水中より(4.3)

も海水面極近傍の水平乱流拡散係数が大きい傾向にあることを説明した．

ここで，l i は粒子群全体の重心と着目粒子 i の距離， L ma x

また，実験結果について，拡散スケールと水平乱流拡散係

は拡散スケールと水平乱流拡散係数の関係が一様に成り立つ数の関係について最小二乗法で近似線を求めた．さらに拡散範囲である． L ma x は流体場の再現スケールで決定される距

28

Anheung

360

18

270

12

180

6

90

0

0
7

8

9

10

Wind direction (Degree)

Wind velocity (m/s)

Oil spill

24

DateWind velocity

Wind direction

図- 5.1: Hebei Spirit 号油流出事故発生場所

360

18.0

270

12.0

180

6.0

90
0

0.0
7

8

9

10

Date

スケールと水平乱流拡散係数の関係が成り立つ空間スケール

Wind velocity

Wind directionと自然外力の再現空間スケールについて考察し，水平乱流拡(b) Gadaeam (観測器海抜高 17m)

散係数を導出する方法を求めた．

図- 5.2: 事故時の風況観測値5.

韓国泰安沖 Hebei spirit 号油流出事故の再現計算
5.1

Wind direction (Degree)

Wind velocity (m/s)(a) Ando (観測器海抜高 27m)
24.0を用いた．密度場は 1,026 kg/m3 で一様に与えた．初期流速

事故概要

場は 0 m/s とした．

2007 年 12 月 7 日 6:50 頃，韓国忠清南道泰安沖約 8 km

風による油の移流は風係数法と数値計算による直接推定で（図-5.1 参照）において錨泊中であったタンカー Hebei Spirit

再現し，計算結果を比較した．吹送流の直接推定は潮汐流計号（香港船籍，146,848 GT）にクレーン台船が衝突した．こ算と同じく富田・柿沼 (2005) の STOC-ML で潮汐流と同時にの衝突によりタンカーのタンク（1 番，3 番，5 番）が破口し，計算を行った．計算条件について，鉛直方向の格子間隔は吹

積み荷の原油，12,547 kL が黄海に流出した．油の流出時間送流を再現するために細かく設定する必要があると考えられ

は 5 番タンクが一番短く 12 月 7 日の 12:00 で止まった．破

るけれども，リアルタイムで計算する時間も考慮し，表-5.2

口寸法の小さい 1 番タンクは長時間に渡り油の漏出が続き，の通り海表面近傍で ∆z=1 m に設定した．それに伴い計算時

翌日深夜 23:40 頃に流出が止まった．

間間隔 ∆t も短くした．それ以外の計算条件は表-5.1 と同じ

韓国泰安沖の沿岸は海岸線方向に潮汐流が強い海域であである．風況データは気象業務視線センターの提供する GPV

り，流出油は広い範囲に移流した．また，事故当時は図-5.2

メソ客観解析データ（MSM）の再解析データを用いた．

のように強い北西風が吹いており，流出した油は 12 月 7 日

油拡散粒子モデルの計算条件は表-5.3 の通りである．油の深夜から 12 月 8 日未明に泰安郡の海岸線に漂着した．その流出位置は藤田他 (2008)，油の流出流量は Lee (2008) を参考

ため環境及び地域経済に多大な打撃を与える，韓国史上最大に設定した．粒子数，計算時間間隔は計算機の制約，計算時の油流出事故となった．なお，事故の詳細については藤田他

間の兼ね合いを考えて決定した．粒子体積は油流出量と粒子(2008) に詳しい記載がある．

数から逆算して決定した．

5.2

計算条件

5.3 計算結果と観測結果の比較

第 2 章で説明した油拡散粒子モデルを用いて，流出油の潮汐流の計算結果について，図-5.3 は Anheung（図-5.1 参

移流，拡散の再現計算を行った．移流外力は，次のように取照）で観測された験潮記録と潮汐の数値計算結果の比較であり扱った．まず，韓国泰安沖である黄海は海流の影響は大きる．図-5.3 における観測と計算の誤差の標準偏差は 0.27 m でくないと考えられる．また，泰安付近では大河川の流入がな

ある．油流出地点周辺の潮汐流観測記録を見つけることがでい．そのため，海流と密度流の影響を無視し，潮汐流と風にきなかったため，潮汐だけの比較にとどまるものの，潮汐のよる表面流のみを再現した．

計算結果は験潮記録をよく再現していると考える．

潮汐流は富田・柿沼 (2005) の STOC-ML を用いて計算し

図-5.4 は 3 つの時刻における油流出範囲の観測結果と漂た．潮汐流の計算条件は表-5.1 の通りである．境界条件は流予測計算結果の比較である．図-5.4 の観測結果において，Matsumoto et al. (2000) のモデルで求めた水位を開境界で与え12/7，8 は Lee et al. (2009) の目視観測，12/11 は衛星 Envisat

た．水深データは GEBCO（大洋水深総図，General Bathymetric

ASAR（http://www.esa.int/）の観測による．数値計算は油の風

Chart of the Oceans）の 30 秒格子データ（http://www.gebco.net/）

による移流の再現に風係数を用いたものと，STOC-ML （富29表- 5.4: 観測結果と計算結果における油の移流の比較

表- 5.1: 潮汐流の計算条件

日時125◦ 10′ 15′′ E

計算範囲〜 126◦ 24′ 45′′ E(観測：画像重心位置)

36◦ 00′ 15′′ N

12/8

12/11

18:00

18:00

10:40

7.6

13.4

17.5

2.0

2.8

2.5

-(数値計算：油粒子平均位置)(km)〜 37◦ 19′ 45′′ N
計算時間

12/7

時間当たりの両者のずれ (cm/s)

2007 年 12 月 4 日 0:00
〜11 日 11:00

70 × 112 × 10

格子数
格子間隔

表- 5.5: 観測結果と計算結果における油の拡散の比較∆x743 m (30 秒)∆y

925 m (30 秒)∆z

5 m(1〜8 層)

油拡散面積 (km2 )10 m(9，10 層)

日時

12/7

12/8

12/11

18:00

18:00

10:40

A:観測結果

31

110

1,195

時間間隔

∆t

10 s

B:計算結果

30

110

1,168

海底粗度係数n0.0026

観測と計算

水平渦動粘性係数

νH

30 m2 /s

結果の面積比

1.04

1.00

1.02

鉛直渦動粘性係数

νV

0.10 m2 /s

A/B

田・柿沼, 2005） を用いて吹送流の直接推定を行ったものを表- 5.2: 吹送流の直接計算における計算条件
格子間隔

時間間隔

∆z

∆t

示している．

1 m (1〜10 層)
2.5 m (11, 12 層)観測結果と風係数法を用いた計算結果を比較すると，韓国

5 m (13〜19 層)

泰安沖では北東〜南西方向の潮汐流が卓越しており，12/7，8

20 m (20 層)において潮汐流により北東〜南西方向に油が移動している．
また，事故時は北西風が卓越していたことから，風により

2.0 s

南東方向に油が移動している様子が再現されている．また，

12/11 の油拡散範囲は観測結果と一致していると考えられる．
計算結果を検証するために，表-5.4 及び表-5.5 のように観表- 5.3: 油拡散粒子モデルの計算条件
油流出位置126◦ 03.01′ E，36◦ 52.01′ N

粒子数

1,200 個粒子体積10 m3 /個

時間間隔

10 s

油流出時間及び

12/7 8:00〜12/9 0:00 211 m3 /hour

流出流量12/9 0:00〜12/11 11:00 61 m3 /hour

測結果と計算結果の比較を行った．解析方法について，観測結果の拡散面積は観測位置全面に油があると仮定し，画像のpixel 数から重心及び面積を求めた．計算結果の重心は粒子位置から求めた．計算結果の面積は，対象海域をメッシュに区切り，粒子が存在するメッシュは油が一様にあると仮定し油拡散面積とした．
まず，油の移流について，表-5.4 のように観測結果は時間経過と共に誤差が大きくなるける．しかしながら，時間当たりで考えると，2〜3 cm/s 程度と，大きな移流速度を持つ海域においても十分な精度で計算ができていると考える．

Tide level (m)

3
2
1
0
-1
-2
-3
12/7

次に，油の拡散面積について，表-5.5 より提案した数値計
算法は実験結果をよく再現できていると考えられる．
次に観測結果と吹送流を直接推定した計算結果を比較す
る．図-5.4 より吹送流を直接推定した場合，南東方向へ向か

8

9
Day
Observation

う風による油の移流がほとんど再現されていないことがわか

10

る．本数値計算では水面近傍の鉛直方向計算格子を 1 m で計算している．少なくともこの程度の格子間隔では水面極図- 5.3: Anheung の験潮記録と潮汐数値計算の比較

30

N(a) 12/7 18:00 時点
N(b) 12/8 18:00 時点
N(c) 12/11 10:40 時点

図- 5.4: 流出油の観測結果（左）と計算結果（中央：風係数法，右：吹送流の直接推定）の比較
近傍の吹送流を再現できていないと考えられる．しかしなが算結果の違いを比較し，各項の影響の大きさを検証する．

ら，水面近傍の吹送流を再現するために格子間隔を細かくす

図-5.5 は油膜自身の特性による油拡散と油の乱流拡散に関ると，時間間隔も小さくする必要があるため，計算機の能力

して，数値計算における考慮の有無と計算結果を比較したも等を考慮すると実用的でない．この結果より，現状では油ののである．計算結果の時間は油流出から 7 時間後である．計風による移流の評価は，数値計算により直接吹送流を再現す

算は油膜自身の特性による油拡散及び油の乱流拡散の両方をるよりも風係数法によって再現するほうが適している．吹送

考慮した場合，どちらかのみ考慮した場合，どちらも考慮し流の直接計算により漂流予測を行う場合は，計算機の性能向ない場合で行った．

上，コーディングの改良等の対策が必要であり，今後の課題

計算結果をみると両方を考慮した場合と比較して，油膜自であると言える．

5.4

身の特性による油拡散のみを考慮した場合（図-5.5 右上）は，
帯状に移流するのみで油粒子がほとんど拡散していない．ま

計算に用いるパラメータの検討

た，乱流拡散のみを考慮した場合（図-5.5 左下）を見てわかる

ここでは本数値計算法における油膜自身の特性による油

ように，ほとんどが乱流拡散の影響で油が拡がっている．油

拡散と油の乱流拡散の計算結果に及ぼす影響について検証

膜自身の特性による油拡散，油の乱流拡散の両者とも考慮しする．油流出事故の観測結果から油膜自身の特性による油拡ない場合（図-5.5 右下）も帯状に移流するのみであった．計散，油の乱流拡散及び移流による見かけ上の油の拡散の影響

算対象海域では潮流による移流の影響が強く，このような海を分けることができないため，数値計算結果における油膜自域においては油膜自身の特性による油拡散の効果は相対的に身の特性による油拡散と油の乱流拡散の考慮の有無と数値計

大きくなく，乱流拡散の効果が大きいということがわかった．

31

渦拡散係数(m2/s)
風向

潮流の方向

Oil spill

図- 5.6: Smagorinsky モデルを用いた水平渦動粘性係数

Diffusion Coef.
DH(m2/s)νH (m2 /s) の推定（2007 年 12 月 7 日 8:00 時点）

図- 5.5: 油膜自身の特性による油拡散と油の乱流拡散の影響．
時間 12/7 15:00（油流出から 7 時間後）

12/7 18:00 12/8 18:00
30

12/11 10:40

20
10
0

右上：油膜自身の特性による油拡散有り，油の乱流拡散無し

12 24 36 48 60 72 84 96
Elapsed time from oil spill (hour)図- 5.7: 水平乱流拡散係数の推定方法で計算された各油粒子

左下：油膜自身の特性による油拡散無し，油の乱流拡散有りの拡散係数の平均値の時系列変化

左上：油膜自身の特性による油拡散有り，油の乱流拡散有り

0右下：油膜自身の特性による油拡散無し，油の乱流拡散無し
る．時間は 12/7 8:00 の結果で，潮位変化が大きく，潮汐流
次に，水平乱流拡散係数 D H の設定方法に関して検討すが早い時間帯を取り出している．水平渦動粘性係数は局所的る．第 4 章において海水面極近傍における水平乱流拡散の推には 101 m2 /s オーダーであるけれども，ほとんどが 100 m2 /s

定方法の提案を行っている．このほかに，乱流モデルを用いオーダーである．ここで，図-5.7 は第 4 章で提案する拡散係て流況を数値計算する際に計算される水平渦動粘性係数 νH数を求める手法で計算された，各油粒子の拡散係数の平均値

が，油の乱流拡散係数 D H と等しいとして，計算に用いるこの時系列変化である．これをみると水平乱流拡散係数は 20とも可能である．そこで， Smagorinsky (1963) のモデルを用〜30 m2 /s で推移しており，平均は 22 m2 /s である．また，韓いた水平渦動粘性係数 νH の空間分布を示し，著者らの水平

国泰安沖の水平渦動粘性係数は 101 m2 /s オーダーであると方向乱流拡散の推定方法との比較を行い，Smagorinsky モデの結果を得ている（松崎・藤田 (2012)）．よって，今回の計ルの活用について考察を行った．

算条件では Smagorinsky モデルを油の乱流拡散の決定に用い

本計算においては水平方向等方乱流を仮定し，∆ は計算格

ることは適さなく，海水面極近傍における水平方向乱流拡散子に依存するとして，Smagorinsky モデルに基づき，式 (5.1)の推定方法を用いるとよいと考えられる．

のように拡散係数の算出を試みた．

最後に，導入した水平乱流拡散係数の導出方法の意義について考察する．図-5.8 は水平乱流拡散係数と拡散面積計算{(
νH = Cs2 ∆x∆y

∂u
∂x

)2

(
+

∂v
∂y

)2

(

+

1 ∂u ∂v+
2 ∂y ∂x) 2 } 21

結果の関係を示している．これを見ると，提案する水平乱流(5.1)

拡散係数の導出方法を用いても，適切な定数である D H =22m2 /s を入力しても結果は大きく変わらない．しかしながら，

ここで，Cs は Smagorinsky 定数，u，v は付録 A で説明する

定数を適切に見積もれなかった場合は，結果が大きく異なる

x ， y 方向の流速である．

ことがわかる．したがって，提案する水平乱流拡散係数の導Smagorinsky モデルを用いて水平渦動粘性係数を求めた結出方法を用いるメリットとして，定数を入力することなく，

果が図-5.6 である．この時の Cs は 0.20 とし，その他の計算

自動的に計算ができるほか，定数の見積もり誤りがなく，不条件は表-5.1 と同じである．なお Cs の理論値は 0.173 であ32

Ratio of oil slick areas(observation) / (Simulation)

2.5
Model
2.0
Diffusion Coef. 1 m²/s
1.5
Diffusion Coef. 10 m²/s
1.0
東京灯標

Diffusion Coef. 22 m²/s
0.5
Diffusion Coef. 100 m²/s

0.0
12/7
18:00

12/8
18:00

12/11
10:40

図- 5.8: 水平乱流拡散係数と拡散面積計算結果の関係

本牧

慣れな計算者が計算を行っても適切に水平乱流拡散係数を設

7/2 10:50
現在位置

定できることが挙げられる．

5.5

7/2 10:05
底触位置
第二海堡

結語

本章は油拡散粒子モデルを用いて，2007 年に韓国泰安沖
図- 6.1: Diamond Grace 号油流出事故発生場所渡辺他 (1998)

で発生したタンカー事故に伴う油流出について再現計算を行い，数値計算法の妥当性を検証した．本章で得られた主要な

モンドグレース号の座礁事故が発生し，約 1,550 kL の Umm

結論を要約する．
まず，事故の概要について説明した．事故の発生した海域

Shaif 原油が流出した．事故当時の気象は快晴で，波高は 0.5では潮汐流が強く，事故時は北西風が強く吹いていたため，

m であった．図-6.2 は，事故時の現場海域の風況観測値であ

これらの自然外力を再現する必要性があることを説明した．

る．事故発生時は南南西の風が約 9 m/s 程度吹いており，その後は南-南南東の風が 3〜8 m/s 程度吹いていた．

次に，計算条件について説明した．自然外力は潮汐流と風流出油は 3 日早朝に大黒埠頭，本牧埠頭，扇島護岸，東扇による表面流を考慮し，それぞれの計算方法について説明

島防波堤の全域に漂着した．流出油のほとんどが 4 日 21 時した．

30 分までに回収されており，東京港，千葉港，木更津港，横続いて，計算結果について説明した．まず，再現した潮汐

須賀港では流出油の漂着は確認されなかった．

は観測結果とよく一致していることを説明した．次に，計算結果について観測結果とよく一致していることを説明した．

事故当初から約 170 kL の油処理剤が散布されている．閉

また，風による表面流について，現状では 3 次元流体モデル

じた系での実験から求められた流出油に対する使用基準量はで再現することは困難であり，風係数を用いて再現するとよ原油の 10〜20%であり，今回の事故で散布された油処理剤はいことを示した．

原油 850〜1,700 kL が流出した場合に散布する量に相当する．
なお，事故の詳細については鶴谷他 (1997) に詳しい記載がまた，計算に用いるパラメータの検討を行った．まず，油ある．

膜自身の特性による油拡散と流れの乱れによる油の乱流拡散について，その効果の検証を行い，今回のように流れが大き6.2 計算条件

く乱流が大きいと考えられる海域では油の漂流の中で乱流拡移流外力は次のように取り扱った．まず，油の流出位置及

散が支配的であることを示した．さらに，開発した数値計算

び移流，拡散の位置が東京湾の横須賀から富津以北のため，法について，水平乱流拡散係数は第 4 章で説明した導出方法

海流の影響は大きくないと考えられる．また，河川流入等にで計算すると，観測結果とよく一致し，乱流モデルで有名な

よる密度流に関して，鈴木 (2013) の資料によると，HF レー

Smagorinsky モデルを用いた再現では，乱流拡散を説明でき

ダーによって計測された湾央から湾口の湾軸方向（北東-南西ないことを示した．最後に，導入した水平乱流拡散係数の導方向）表層流速の潮汐残差流の 7 月平均値は南西向きに 0.05

出方法の意義について説明した．

6.

m/s である．同様に鈴木 (2013) の資料によると，最も密度流の影響が強い位置の一つと考えられる羽田における無風時，

東京湾 Diamond Grace 号油流出事故の再現計算

すなわち主に密度流による湾軸方向表層流速は 0.078 m/s と

6.1

事故概要

見積もられている．それに対して，事故時の第二海堡におけ1997 年 7 月 2 日 10 時 5 分に東京湾中ノ瀬付近（北緯 35◦

る潮流は 0.30 m/s であり，海上 10 m 風速が 8.0 m/s である

21′ ，東経 139◦ 43′ ，図-6.1 参照）において，油送船ダイヤことから，風による表面流の平均値は 0.24 m/s と見積もられ

33

360

12

270

8

180

4

90

0

0
2

3
Date
Wind velocity

4

表- 6.1: 潮汐流の計算条件

WInd direction (Degree)

Wind velocity (m/s)

16

計算時間

1997 年 6 月 28 日 10:00〜
7 月 4 日 11 日 10:00
128 × 188 × 10

格子数
格子間隔

∆x

500 m

∆y

500 m

∆z

5 m(1〜8 層)
10 m(9 層)

Wind direction

50 m(10 層)

360

12

270

8

180

4

90

0

0
2

3
Date
Wind velocity

4Wind direction (Degree)

Wind velocity (m/s)(a) 第二海堡

16

時間間隔

∆t

10 s

海底粗度係数n0.0026

水平渦動粘性係数

νH

30 m2 /s

鉛直渦動粘性係数

νV

0.01 m2 /s

表- 6.2: 計算に使用した風向風速観測点の位置と高度
緯度Wind direction

経度

T.P. からの
風速計

観測地点名

(b) 本牧図- 6.2: 事故時の風況観測値
る．これに従うと事故発生時の海表面の流況に与える密度流

度

分

度

分

高度 (m)

東京灯標

35

33.97

139

49.68

25.0

第二海堡

35

18.75

139

44.63

29.0

本牧

35

26.33

139

41.38

56.0

の影響は 8〜13%程度であり，潮汐流と風による表面流が卓越している．そこで，本計算では計算の簡便性を優先し，移流外力は潮汐流と風による表面流のみとした．
潮汐流は 5.2 節と同様に富田・柿沼 (2005) の STOC-ML を

(
W10 =

用いて計算した．潮汐流の計算条件は表-6.1 の通りである．

(6.1)
(6.2)

m=1

Rm
Rsum

(x − x m) 2 + (y − yWx,y =

提供する 500 m 格子データ（J-EGG500）を用いた．密度場

Rm = √

は 1,026 kg/m3 で一様に与えた．初期流速場は 0 m/s とした．
風による油の移流は風係数法を用いて再現した．風況デー

3
∑

Wm
1

Rsum =

タは観測値を用いて以下のように作成した．風は海上と陸上

) 17
Wz

境界条件は Matsumoto et al. (2000) のモデルで求めた水位を東京湾口で与えた．水深データは日本海洋データセンターの

10
z

3
∑

Rm

m)

2

(6.3)

(6.4)

m=1

では全く異なることから，海上にある東京灯標，第二海堡と本牧船舶通航信号所（図-6.1，表-6.2 参照）の計 3 地点におけ

ここで，m は観測地点，x ，y はそれぞれ東西，南北の座標，

る 1 時間ごとの観測値を用いた．また計測高度から式 (6.1)

x m ， ym は観測地点 m の座標，Wm は観測地点 m におけるのように 1/7 乗則を用いて，海上 10 m 風速に換算した．風

風速， Rm は観測地点 m と対象地点の距離の逆数， Rsum は

速の空間分布は式 (6.2) から式 (6.4) を用いて，距離の逆数に3 つの Rm の和である．
油拡散粒子モデルの計算条件は表-6.3 の通りである．油の比例するよう重みをつけて風速場 Wx, y を再現した．

流出位置は渡辺他 (1998) を参考に設定した．流出流量は横浜市港湾局 (1998) の事故時の記録を参考に設定した．粒子数，計算時間間隔は計算機の制約，計算時間の兼ね合いを考えて決定した．粒子体積は油流出量と粒子数から逆算して決
定した．

34

3.5

粒子数

1,500 個粒子体積

1 m 3 /個計算時間間隔10 s

油流出時間及び

7/2 10:05〜10:50 1,333 m3 /hour

流出流量7/2 10:50〜13:05 222 m3 / hour

Tide level (m)表- 6.3: 油拡散粒子モデルの計算条件

2.5
1.5
0.5
2

3
Day
Observation30

4Calculation

Simulated current
velocity (cm/s)

図- 6.4: 事故時の横浜新港の潮位観測値と計算値の比較

20あり，計算結果は観測結果とよく一致している．図-6.3 及び図-6.4 より，潮汐流と潮位の再現結果は良好であることから，

10

同じ計算条件で行った事故時の潮汐流計算も精度よく行われているものと考える．

00

次に，流出油の拡散観測結果と計算結果を比較する．ここ

10
20
30
Observed current velocity(cm/s)

で精度検証においていくつか問題がある．まず，坂本 (1998)によると，関係各機関の油濁対応が早かったことから，3 日(a) 潮流楕円の長軸長さ間の防除活動により浮遊油の大部分が回収された．また 6.1

Simulated long axis
direction (degree)0.0

節で述べた通り油処理剤による対処も行われており，原油のみの状況と異なり流出油が水中分散をしやすい状況にあったと考えられる．さらにダイヤモンドグレース号から流出した-45.0

Umm Shaif 原油は揮発性が高い．NOAA の提供する流出油の性状変化ソフトウェア ADIOS2（http://archive.orr.noaa.gov/

book_shelf/538_adios.pdf ）によると，Umm Shaif 原油の海水

-90.0
-90
-45
0
Observed long axis direction(degree)温 25 ℃，波高 0.5 m，風速 10 m/s，海域流速 0.25 m/s の条件
下における，24 時間後の油の蒸発率は約 37%である．蒸発(b) 潮流楕円の長軸角度を考慮していない本モデルでは油の量を単純に比較すること

図- 6.3: 第二海堡の観測値と計算値の各分潮の比較はできない．今回の計算結果はこれらの影響について考慮されていない．

6.3

以上のことから，観測結果と計算結果との比較は，油回収

計算結果と観測結果の比較

潮汐流計算の再現精度について検討する．事故当日の連続作業等の影響が大きいと考えられる油量や拡散面積の評価はした潮汐流の計測データは見つからなかったため，1 年間の行わず，油回収作業等の影響が比較的小さいと考えられる移東京湾の潮汐流計算を行い，第二海堡での観測値と計算値の流のみで行うこととする．
図-6.5 は 7 月 2 日 14:30 時点及び 7 月 3 日 12:00 時点の油分潮を比較することで，潮汐流計算の精度を検証する．
図-6.3 は第二海堡における観測値と計算値の各分潮の比較

流出範囲の観測結果と同時刻の計算結果の比較である．海面である．観測値は関東地方整備局東京湾口航路事務所の設置の油膜の観測結果はいくつか発表されており（例えば，坂本,

する海象計の 2009 年 1 年間で測定された測定水深 10 m の

1998；市橋・赤松, 1997），その中で海上保安庁が発表した拡流速データを用いた．計算値は海象計の設置されている位置

散結果（横浜市港湾局, 1998）は拡散範囲のほか，油膜の外で，測定水深が同じ計算メッシュのデータを用いた．観測値

見による油膜の分類も図示されているため，この資料と比較

及び計算値の各分潮は，田中 (2008) の方法を用いて主要 8 分した．図-6.5 の観測結果で示される油膜の厚さは，表-6.4 の潮（Q1，O1，P1，K1，N2，M2，S2，K2）の潮流楕円の長油膜の外見との対応を元に作成した．
図-6.5 の 7 月 2 日 14:30 時点において，事故当時は南南西軸の長さと長軸の角度の結果を示した．図より計算値は観測

風が卓越しており，潮汐が上げ潮であったことから，油層が値とよく一致していると考える．
図-6.4 は事故時の横浜新港の潮位観測値と計算値の比較

北東方向への移流がよく一致しており，拡散する様子もよくである．計算結果と観測結果の誤差の標準偏差は 0.03 m で再現されている．また，7 月 3 日 12:00 時点において，大黒

35

N

油膜の
厚さ
(μm)

B

底触位置

底触位置

(a) 7/2 14:30 時点
N

油膜の
厚さ
(μm)

Cの中にスラッジ
Eの中に一部C
E

底触位置

底触位置

(b) 7/3 12:00 時点

図- 6.5: 油流出範囲の観測結果（左）（横浜市港湾局, 1998）と同時刻の計算結果（右）の比較
埠頭，本牧埠頭，扇島護岸，東扇島防波堤といった東京湾西側沿岸に油が移流し，漂着する様子が再現されている．特に表- 6.4: 海上災害防止センター (2000) による油膜の厚さと油油膜の厚い部分が沿岸部にある様子がよく一致していると考膜の外見の関係

える．鶴谷他 (1997) や坂本 (1998) によると，蒸発せず，回収

油膜の厚さ

油膜の外見

されなかった流出油の大半が大黒埠頭，本牧埠頭，扇島護岸，

標示

東扇島防波堤といった，川崎市，横浜市の沿岸部に移流，拡

(µm)
2.0

油膜の色が黒ずんで見える状態

A

1.0

油膜がくすんだ褐色を呈する状態

B0.3

水面に明るい褐色の帯がはっきりC散したとあり，計算結果はそれをよく再現していると考える．

6.4 結語
本章は油拡散粒子モデルを用いて，1997 年に東京湾で発

見える状態

0.15

水面がほんの少し褐色に色づいて生したタンカー事故に伴う油流出について再現計算を行い，D

数値計算法の妥当性を検証した．本章で得られた主要な結論

見える状態

0.1

水面が銀色にキラキラ光ってを要約する．

E

まず，事故の概要について説明した．事故の発生した海域

見える状態

0.05

光線の条件が最も良い時にでは潮汐流が強く，事故時は南風が強く吹いていたため，こ-れらの自然外力を再現する必要性があることを説明した．

かろうじてキラキラ光る油膜が次に，計算条件について説明した．自然外力は潮汐流と風

見える状態

による表面流を考慮し，それぞれの計算方法について説明した．
続いて，計算結果について説明した．まず，再現した潮汐

36

360

18

270

12

180

6

90

0

0
0

6

12
Time
Wind velocity

18

観音埼
港湾空港
技術研究所

WInd direction (Degree)

Wind velocity (m/s)

24

(1)
剱埼
油流出位置

図- 7.2: 調査地点と油の漂着範囲．地点 (1) から (4) は油確
認，(5)，(6) は油確認されず． (地図出典：国土地理院)

Wind directonて，事故の翌日である 3 月 19 日の午後に千葉県において調(a) 観音崎（計測器海抜高 83 m)

270

12

180

6

90

0

囲を示している．図中の (1) から (4) の海岸線において油のWInd direction (Degree)

Wind velocity (m/s)

査を実施した．図-7.2 は調査地点と油の漂着を確認できた範

360

18

漂着が確認できた．また，(5)，(6) において油の漂着が確認できなかった．(1) 以南の海岸線は未調査であったけれども，(1) 地点において千葉県館山から車で北上しながら油被害状
況調査中の千葉県職員から，(1) の地点より南側では油の漂着は確認できず，(1) の地点において初めて漂着を確認したとの情報を得た．よって，3 月 19 日午後時点の油の漂着範
囲は図-7.2 の線で示す 10 km 程度と考えられる．

0
0

6

12
Time
Wind velocity

18

油が漂着した時間に関して，東京湾フェリーの職員から 1824

日の夕刻にフェリーの入出港である金谷港に油が漂着したとWind directonの情報を得た．

(b) 剱崎（計測器海抜高 38 m)

図-7.3 から図-7.6 は現地調査時の現地の様子である。18 日

図- 7.1: 事故時の風況観測値の午前中より風浪が強くなり、図-7.3 のように護岸を超えて陸上に油が拡がっているところもあった。図-7.4 は油防除中

流は観測結果とよく一致していることを説明した．最後に，の千葉県職員の様子である．なお，保険会社の雇った業者に計算結果の油の移流位置が観測結果とよく一致していること

より流出油や流出物のほとんどが回収された．図-7.5 のようを説明した．

7.

(3)
(2)
油漂着域

24

24

(6)
(4) (5)に漂着油はエマルション化しており、Beagle3 号から C 重油が流出したものと考えられる．地点（5）では油の漂着は確

東京湾 Beagle3 号油流出事故の予測計算認されなかったけれども、図-7.6 のように油に汚染された鳥

7.1

事故概要

が何羽もいた。

2014 年 3 月 18 日午前 3:20 頃、神奈川県三浦市の沖合約

なお海上保安庁の報告によると、Beagle3 号からの油の流6 km の浦賀水道において、パナマ船籍貨物船 BEAGLE3 号

出は 3 月 19 日以降も継続しており、20 日以降は今回示した（12,630 t）と韓国船籍貨物船 PEGASUS PRIME（7,406 t）が漂着範囲以外にも油の漂着が確認されている。しかしなが衝突する事故が発生した。衝突により BEAGLE3 号が水深約

ら、著者が 3 月 25 日午後に神奈川県三浦市の沿岸で行った100m の海底に沈没し、中国人乗組員 7 人が死亡、2 人が行

現地調査によると、漂着量はタールボールやごく薄い油膜が方不明となっている。また、同船より燃料油の一部が流出し点在する程度であった。そのため、図-7.2 で示した漂着範囲

た。第三管区海上保安部によると沈没時に船内に残存したと以外では、千葉県沿岸のように大量の油は漂着しなかったと思われる油の量は C 重油約 357 kL、A 重油約 42 kL であっ考えられる。

た。事故時の天候は晴れで波も穏やかであり、視程は観音埼において 10,000 m であった。図-7.1 のように事故後は現場付

7.2 計算条件

近においておよそ 5〜20 m/s 程度の強い南西風が吹いていた。

本計算は 5.2 節及び 6.2 節と異なり，8.3 節で説明するよう

著者は油拡散初期の海岸線の油漂着状況調査を目的としに，即時に予報できる方法で行った．移流外力は，海流，潮汐流及び風による表面流のみとした．

37

図- 7.3: 調査地点（1）図- 7.5: 調査地点（3）図- 7.4: 調査地点（2）図- 7.6: 調査地点（5）

計算は油拡散粒子モデルを GUI に組み込んだシミュレー
タ（名称：漂流油捕捉システム）で行われた．基本的な計算

nt∑

方法は 5.2 節及び 6.2 節と変わらない．潮汐流の再現方法のRvel = vt

み異なり，付録 A でその方法を説明する．
油の風による移流は風係数法を用いて再現した．風況デー(Uoit − Uo)(Ucit − Uc)

it=1
nt
∑

(Uoit − Uo) 2

it=1

タは気象業務支援センターの配信するメソ数値予報モデル
油拡散粒子モデルの計算条件は表-7.1 の通りである．油のRdir =

流出位置は AIS 情報（Automatic Identification System，自動船
舶識別装置，AIS 情報出典: https://www.marinetraﬃc.com/jp/）(7.1)(Ucit − Uc) 2

it=1
nt
∑

GPV(MSM) の予報値を用いた．

v
t nt
∑

|Uoit ||Ucit | cos θ it

it=1
nt
∑

(7.2)
|Uoit ||Ucit |it=1を基に決定した．なお，後日得られた情報から，AIS 情報が

ここで，nt は解析データ個数，Uoit ，Ucit はそれぞれ時間 it の正確であったことが確認できた．油の流出流量は不明であっ際の潮汐流の観測値及び計算値，Uo，Uc はそれぞれ潮汐流のたため，表-7.1 の通り粒子数のみを設定した．流出時間も不

観測値及び計算値の平均値，θ it は時間 it の観測値と計算値の明であったため，表-7.1 の通りに仮想して設定した．粒子数，ベクトルのなす角である．解析結果はそれぞれ Rvel =0.920

計算時間間隔は計算機の制約，計算時間の兼ね合いを考えて

及び Rdir =0.959 であった．これより，潮汐流の再現結果は決定した．

高いと考える．

7.3

次に，図-7.7 は油の漂流予測計算結果である．これを見る

計算結果と観測結果の比較と，漂流位置が 1 km から 2 km ほど北方にずれている程度

調和定数から再現した潮汐流の精度について検証する．東で，漂着位置の観測結果と計算結果はよく一致している．ま京湾の第二海堡の 2010 年 1 年間の 20 分間隔の流況観測結果

た，漂着時間について，夕刻であったとの目撃証言とよく一と再現結果について，流速の相関係数 Rvel を式 (7.1)，流向

致している．これより，計算結果は実際の油の移動とよく一の内積 Rdir を式 (7.2) にて評価した致していることがわかる．
また，本計算結果は 8.3 節で説明する移流外力のデータベースを用いて，事故後に即時に行ったものである．油流出事故

38

油の漂着範囲

油粒子
油流出位置(a) 06:20(3 時間後)(b) 09:20(6 時間後)(c) 12:20(9 時間後)(d) 15:20(12 時間後)(e) 18:20(15 時間後)

(f) 21:20(18 時間後)図- 7.7: 流出油の漂流予測計算結果
まず，事故の概要について説明した．著者が行った海岸線

表- 7.1: 油拡散粒子モデルの計算条件
油流出位置139.71590◦ E，35.09753◦ Nへの油漂着状況調査の結果，事故の当日に油が多く漂着した

粒子数

100 個

のは千葉県沿岸の 10 km 程度の海岸線であったことを説明し

時間間隔

60s

油流出時間

3/18 3:20 から 3 時間一定量連続流出

た．また，漂着した時間は夕刻であったことを説明した．
次に，計算条件について説明した．自然外力は海流，潮汐流及び風による表面流を考慮し，それぞれの計算方法について説明した．本計算は即時の予報を目的としているため，潮汐流の再現方法について説明した．

が発生してから計算に必要なデータを用意するのでなく，事続いて，計算結果について説明した．まず，潮汐流を精度故の連絡，情報の入手と共にすぐに数値計算を実施できる予よく再現できたことを示した．次に計算結果について油の報体制が整ったことは大きな成果であると考える．

位置が観測結果とよく一致していることを説明した．また，

7.4

結語

油漂着時間についても実際とよく一致していることを説明

本章は，油拡散粒子モデルを用いて，2014 年 3 月に東京した．

湾浦賀水道で発生した貨物船衝突沈没事故に伴う油流出につまた，本計算結果により，油流出事故が発生してからすぐいて再現計算を行い，数値計算法の妥当性を検証した．本章に計算を実施できる予報体制が整ったことが示された．

で得られた主要な結論を要約する．

39

8.

議論
8.1

本計算モデルは有用性が高いと考える．また，第 7 章の結果
から，即時の漂流予測を行える状況であることが示された．

既往の研究からの進展

またその計算精度も十分に高く，本計算モデルの有用性が示

本論文において既往の研究から進展させた内容は，油膜自されたと考える．

身の特性による油拡散と油の乱流拡散の予測方法についてで次に，海上に流出した油の風による移流に関して考察する．

ある．

海表面の油は海表面極近傍の吹送流の影響を大きく受ける．

油膜自身の特性による油拡散について，Fay（1969, 1971）

3 次元流体モデルを用いて海表面極近傍の吹送流を再現する，のモデルの適応を平面 2 次元油拡散の室内実験結果から議ためには，鉛直方向の計算格子を小さくとる必要がある．そ論した．実験において油膜は厚さが 10 −4 m オーダー以下での場合 CFL 条件により計算時間が膨大になるため，実務に拡がっており，その際の平面 2 次元方向の拡がり速度は Fay

適応できない．海表面の油の移流は海表面の極近傍の海水のの提案する第 3 段階の表面張力-粘性領域のモデルで説明で動きに大きく影響を受けるため，油の移動を予測する際は，きることを示した．また，従来行われていなかった，Fay の

表層数 mm から数 cm オーダーの表面流を考慮する必要があ

第 3 段階の表面張力-粘性領域における拡がり係数を実験結る．そこで本研究では油の風による移流の再現に風係数法を

果から導出した．

採用した．風係数法は簡易な方法であるものの実際の油の移

次に，油膜自身の特性による油拡散について，Fay のモデ

動をよく再現しており，現場における有用性が高く，有効なルに基づいた新しい数値計算モデルを開発した．開発した計算方法であると言える．

モデルは，表面張力によって拡がる効果を計算するランダム

また，本計算においては風による移流におけるコリオリ力

ウォークモデルと，重力によって拡がる効果を計算する斥力

の評価をしていない．コリオリ力を評価しなくても計算結果モデルの 2 種類である．ランダムウォークモデルについて，

は観測結果とよく一致しているため，その影響は大きくない

静水及び定常流条件下における油の拡がりの実験結果と比較ものと考えられる．しかしながら，予測精度の向上，特にコを行った．斥力モデルについて，Fay のモデルが成り立つとリオリ力の影響の大きいと思われる外洋においては，コリオ仮定し比較を行った．

リ力の評価方法を検討する必要があると考える．数値計算に

乱流拡散について，数値計算モデルを用いた計算で必要とよる予測においてコリオリ力を評価する方法として，例えば

なる水平乱流拡散係数について，導出した方法は擬似油を用

3 次元流体モデルで鉛直方向に可能な限り細かく格子をとり，いた実海域実験を行い，流出油のように海水面極近傍に限っ

吹送流場を計算した際の計算結果をコリオリ力として油拡散

た水平乱流拡散係数を導出する方法を開発した．また，導出

粒子モデルに与える方法が考えられる．どのように計算をすした拡散係数の計算方法について，その適応範囲を議論した．

べきか今後検討していく必要がある．

本論文において既往の研究と異なる点としては，実際の油次に，本論文で提案する数値計算法の適応範囲について述

流出事例 3 件と数値計算結果の比較を行い，数値計算法の有べておく．

効性を検証している点である．本論文では，2007 年に韓国

本数値計算法の検証において，風速は大きくても 15〜16

泰安沖で発生した油流出事故の再現計算を行い，油の移流重m/s 程度であり，それ以上の風速条件下における油の移流及心及び拡散面積について比較を行っている．また，1997 年

び拡散の再現性について検証を行っていない．実海域実験にに東京湾で発生した油流出事故の再現計算を行い，油の移流おける疑似油を用いた乱流拡散実験においても，平均風速は位置と観測結果の関係について議論している．さらに，2014

10 m/s 以下で行っている．そのため，本数値計算法において

年 3 月に東京湾で発生した油流出事故の予測計算を行い，移確認の取れている適応範囲は 10 m/s であり，それ以上の強流重心が 1〜2 km の誤差に収まっており，移流速度がほぼ一風条件下における適応性は検討されていない．

致していることを確認した．

8.2

本数値計算法では波の影響は考慮されていない．波による質量輸送（ストークスドリフト）を予測するには，波の周期，

数値計算法の考察

波高，波長及び波向のスペクトルと現場海域の水深の情報が

本節では油拡散粒子モデルの有用性，適切な計算方法，適必要であり，即時の予測を行うためには，データベースを構

応範囲及びモデルの効果について考察する．
第 5，6 章の結果から，必要な移流外力を再現し，油拡散築する必要がある．また，風波と風による移流の関係につい粒子モデルを用いて計算することで，海表面の流出油の移流て明らかとなっていない．波が発達するにつれて波の質量輸

及び拡散を再現できることを示した．韓国泰安沖と東京湾の送の影響が大きくなるけれども，その際に風係数は 0.03 を油流出に関して，同じモデルを用いて同様の計算条件で数値維持するのか，波の発達にエネルギーが用いられ，風係数が計算を行い，共に精度のよい再現結果が得られたことから，小さくなるか，といった考察が現状では不十分である．この40

表- 8.1: 油拡散粒子モデルの適応範囲及び考え方
自然外力

適応範囲，計算方法

考え方

風速10 m/s まで適応を確認．

15 m/s 程度以上の風速における適応は（第 4 章実海域実験）

まだ確認していない．

15 m/s 程度 まで適応を確認．
（第 5 章油流出事例）
波

即時漂流予測のため，第 5 章の結果から，波による輸送の考慮されていない．

影響は風係数に含められることとする．

河川の影響考慮していない．

海域によっては考慮する必要がある．

海岸線での海岸線に漂着した油は再海域によっては海岸線情報を整備し，油の漂着流出しないと仮定して漂着モデルを開発して，導入する

計算している．

必要がある．

点については，検討が必要であると考える．一方で，第 5 章あると考えられる．

で説明した韓国泰安沖 Hebei spirit 号油流出事故の際は，波高以上の油拡散粒子モデルの適応範囲及び考え方についてが 3〜4 m 程度であった（藤田他, 2008）．しかしながら，計は，表-8.1 の通りにまとめられる．

算結果は観測結果とよく一致している．このことから，波に

次に，油膜自身の特性による拡散と乱流拡散のオーダーによる輸送の影響は風係数に含まれていると考えられる．よっついて考察する．油膜自身の特性による油拡散については，て即時の予報を目的とする際には波の影響を無視してしまっ

第 3 章で述べた通り，油が薄くなった状態においては正味のても，大きな問題はないといえる．

表面張力によって拡がる．正味の表面張力は式 (3.5) で表さ

油拡散粒子モデルを用いた計算を行う場合，計算対象とすれる通り，水-大気，油-大気，油-水の表面張力によって決定る海域がどういった移流外力が支配的であるか，十分考慮される．油回収の主な対象となるのは，性質が安定していてする必要がある．海上では陸上と比較して風が強いため，ほ蒸発しにくい C 重油であり，著者が試験した結果によると，C

とんどのケースで風の影響を考慮する必要があると考える．

重油の油-大気，油-水の表面張力はそれぞれ 26.8，26.2 mN/m

また，今回再現計算対象とした韓国西岸や東京湾といった潮であるため，正味の表面張力は 18.1 mN/m である．乱流拡散位変化が大きいところでは，潮汐流を精度よく再現することについては，ASCE Task Committee on Modeling of Oil Spills ofが必要である．さらに，河川流も海域に応じて再現する必要

the Water Resources Engineering Division (1996) に従い，水平がある．例えば 2010 年のメキシコ湾油流出事故の際は，ミ乱流拡散係数を D=1，10，100 m2 /s として流れの乱れスケー

シシッピ川河口において，風によって岸方向へ移動する流出

ルを表現した．図-8.1 は油膜自身の特性による C 重油の拡散

油が河川流によって妨げられたという事例がある（白石他,と油の乱流拡散の時間変化を示している．これを見ると，乱2011）．第 6 章で示した東京湾のケースにおいて，密度流がれの小さい D=1 m2 /s のときの乱流拡散と油膜自身の特性に

海表面の流れ場に与える影響は約 1 割であるとして考慮しな

よる油拡散のオーダーがほぼ同程度である．そのため，外洋

かったけれども，密度流が卓越する海域ではこれを考慮するなどのように流れの乱れスケールが大きい海域においては乱必要がある．

流拡散が卓越すると考えられる．一方，港湾内のように流れ

本計算において海岸線等における境界条件は前述の通り，の乱れスケールが小さいときは油膜自身の特性による油拡散

海岸線に近づいた際は動きにくくし，漂着した際は油を再流の影響も無視できない大きさであり，その影響を考慮した予出させないようにしている．しかしながら，防波堤，砂浜，測が必要となると考えられる．

礫海岸など，海岸線の状況により漂着後の再流出の度合いが最後に，油膜自身の特性による油拡散の計算モデルの効果

異なることが考えられる．これらを考慮した各パラメータのに関して考察する．油流出事故の観測結果から，これらの効設定は現状ではできていない．本論文では油回収船による油果を切り分けることができない．また，韓国と東京湾の事例

回収支援を目的としており，特定の海域に絞らず日本全国ではいずれも潮流と風が強い条件で油が連続流出しているた計算が可能であることを目的としており，沿岸域付近の流れ

め，見かけ上の拡散の影響が大きく，油膜自身の特性によるについて個別に検討はしていない．そのため，特定の海域に

油拡散の影響が小さくなっている．そのため今回の油流出事

おいて特に精緻に計算を行いたい場合は，密度流及び海岸線

例は油膜自身の特性による油拡散のモデルの有効性を示しにの境界等の影響について，数値計算モデルに組み込む必要が41

Oil slick radius (m)

20000
15000
10000
5000
0
0

1

2
Elapsed time (day)3

Spreading (Heavy-C oil)

Diffusion (D=1m²/s)

Diffusion (D=10m²/s)

Diffusion (D=100m²/s)

図- 8.1: 油膜自身の特性による C 重油の拡散と油の乱流拡散の拡散半径の時間変化
図- 8.2: 漂流油捕捉システムの概要
くい状況であり，実海域の現象と比較した定量的な評価ができていない．

PC を操作できる者であれば，数値計算を実施できるような

今回の 3 つの油流出事例以外に，風や流れの小さい海域に状況にある．なお，数値計算の実施は LinuxOS で動作する

おける油流出事故の観測事例があれば，再現計算を行い観測コンピュータでも可能である．

結果と比較できるけれども，そのような油流出の観測事例が8.4 油回収船の油回収操船支援を行う漂流油捕捉システ

見つからず，本計算モデルのさらなる有効性を示すことがでムの開発きない状況にある．そのため今後そのような事故事例を収集

漂流油は海上において見えにくく，遠方の油を目視で確認し比較することで，モデルの妥当性を評価する必要がある．

8.3

することは困難である．油をすぐに発見できずに海上を航即時の予報体制の構築

行していると，回収作業開始までに時間がかかり，その間に海上において油流出事故が発生した際に，油拡散粒子モデ

流出油はさらに広範囲に拡がってしまう．よって，油回収船

ルを用いてすみやかに移流及び拡散の予測を行うために，港が現地付近に到着した後に，自船で海上漂流油を検出する技

湾空港技術研究所において移流外力のデータベースを構築し術が求められていた．また，流出油の検出情報を基にシミュている．表-8.2 は，2014 年 9 月時点及び 2015 年度以降の港レーションを行うことで，精度の高い予測結果が得られるこ湾空港技術研究所における，流出油の移流に影響する海流，とも期待される．

潮汐流，風のデータベース構築状況である．海流，風のデー

そこで，国土交通省北陸地方整備局の所有する大型油回収

タは気象業務支援センターの配信サービスを利用している．

船「白山」を対象に漂流油捕捉システムを開発した．漂流油海流に関して，2015 年度以降から使用するデータの変更を

捕捉システムは，図-8.2 のように船舶用レーダーを用いた流

予定しており，より詳細な海流の再現が可能となる．風に関出油検出装置と，油拡散粒子モデルを統合し，リアルタイム

して，2015 年度以降から GPV-MSM の使用に加えて，GSMに処理，表示することにより油回収時の操船支援が可能とな（日本域）の使用を予定しており，空間的には粗いけれどもるシステムである．

1 週間程度先までの長期の再現も可能となる．潮汐流は海域

漂流油検出装置は船舶レーダーの受信信号をレーダー信号の場所ごとに周期的に変化しているので，付録 A で説明す解析装置で解析することで，海表面の油膜を検出することがる通り，一度対象海域の潮汐流をシミュレーションし，位置できる．レーダー信号解析装置はノルウェーの Miros 社の製（緯度，経度），時刻及び流速の関係から，潮汐流による移流品である（Egset 他, 2007）．流出油の検出範囲はレーダーの速度が再現できるシステムを構築している．

性能（最大 4 km）とレーダーアンテナ設置高さで決まり，大流出油の漂流シミュレーションを即時に行える海域は，東型油回収船「白山」のアンテナ設置高さ（29 m）だと，半径経 120-150 度，北緯 22.4-47.6 度の範囲である．これは表-8.2

約 3.3 km の流出油の検出が可能である．

の通り，海流，潮流，風のなかでデータ配信領域が最も狭い

流出油の検出原理について簡単に説明する．この装置は風風によるものである．シミュレーションは WindowsOS で動により発生するさざ波からの反射波を解析して流出油の位置

作するコンピュータにおいて，GUI で操作及び設定を行うを判定している．波長が約 3 cm である X バンド帯レーダー

ことができるようソフトウェアを作成しており，Windows のの電波の場合，1.5 cm 程度の波長の波が発生している場合に42

表- 8.2: 移流外力のデータベース構築状況
海流 Uoi

海流 Uoi

潮汐流（東京潮汐流（左の（2014 年度まで）（2015 年度から

湾，伊勢湾，

海域以外）

風 Uwi

使用開始）

瀬戸内海）Ut i

Ut i

データ

日本近海海流海水温・海流富田・柿沼 (2005)

Matsumoto et al.

メソ数値予報モデル

名

予報格子点資料

予報格子点資料のモデル(2000) のモデル

GPV (1) MSM，(2) GSM（GSM は 2015 年度から使用開始）

データ

気象業務支援気象業務支援モデルを用いてMatsumoto et al.

気象業務支援

配信元

センター

センター

独自に計算．

(2000) のモデルをセンター用いて再現．
配信

東経 120〜160 度

東経 116.85 度

東京湾，

東経 110〜165 度(1) 東経 120〜150 度

領域

北緯 20〜50 度〜159.75 度

伊勢湾，

北緯 20〜65 度

北緯 22.4〜47.6 度

北緯 14.85 度

瀬戸内海(2) 東経 120〜150 度〜49.7 度格子

1/4 度

1/10 度

北緯 20〜50 度

1/12 度

270 m

間隔(1) 東西方向 1/16 度
南北方向 1/20 度(2) 東西方向 1/4 度
南北方向 1/5 度
配信

月に 3 回，1 か月

1 日 1 回，30 日

ソフトウェア

ソフトウェア

港湾空港技術研究所の日時

先までの 10 日毎

先までの 1 日毎に組み込んでに組み込んでサーバーで受信．

予報値を港湾空港

予報値を港湾ある．

ある．

(1)1 日に 8 回，39 時間先技術研究所のサー

空港技術研究までの 1 時間毎の予報値．

バーで受信．

所のサーバーで(2)1 日 4 回，84 時間先

受信．

までの 1 時間毎，及び

1 日 1 回 264 時間先までの 3 時間毎の予報値．

強い後方散乱波を受信することが出来る．波長が 1.5 cm 程の挙動を想定した訓練が可能となるといった利点がある．な度の波は，図-8.3(a) のように海面を吹く風により発生するさお，漂流油捕捉システムは 2014 年度に北陸地方整備局大型ざ波に多く含まれているので，さざ波の発生時に後方散乱波

油回収船「白山」に搭載される予定である．

を受信することが出来る．そのため，油膜のない海面におい9. 結論てはさざ波によりレーダーの方向に反射して明るく写る．
一方，油膜の存在する海面では図-8.3(b) のように油膜によ本研究では，油回収船の油回収操船支援を目的とした，海上

りさざ波の発生が抑制されるため海面が平らになる．そのたに流出した油の漂流予測を行う油拡散粒子モデル，OIL-PARI

め油膜のある海面においては鏡面反射によりレーダーの方向を開発した．油拡散粒子モデルは海表面の油に関して，海表に返ってくる反射波が少ないので暗く写る．これを画像処理

面の流れによる油の移流，油膜自身の特性による油拡散，油により，より鮮明にして油膜の有無を判別している．

の乱流拡散を計算できる．

以上で説明した漂流油捕捉システムを用いると，広大な洋本研究で既往の研究より進展した部分は，第 3 章で説明し上を常に移動する漂流油の位置を捕捉し，風や海流等の自然た油膜自身の特性による油拡散，第 4 章で説明した油の乱流条件を基にした漂流油の位置予測をリアルタイムで一元的拡散について検討した点である．また，本研究の特徴としてに行うことができる．そのため，船長が針路と船速を決定す

第 5 章から第 7 章において数値計算法の有用性を示している

るのに役立ち，油回収開始までの時間を短縮できる．また，点が挙げられる．

漂流油捕捉システムを利用した油回収訓練により，実際の油油膜自身の特性による油拡散については，室内実験を行い，

43

は流出油の平面 2 次元の挙動を予測できるけれども，3 次元の挙動は予測できない．そのため，沈船や海底油井からの油流出の漂流予測が困難である．また，本計算モデルでは流出
油の蒸発や，水中に分散，溶解する現象を計算できない．海外では分散剤を用いた流出油対策が主流になりつつあり，分散剤による油の水中分散やその環境影響の評価なども行う必要が出てくる可能性がある．さらに，2011 年の東日本大震
災で顕著であった，自然災害による港湾施設からの油流出に伴う被害予測について，今後取り組んでいく必要がある．今(a) 水表面に油がない状態

回の研究成果はそれら応用的な対策へ欠かせない基礎的研究とも位置づけられ，海水面の油の挙動の解明が少なからずできたのではないかと考える．

(2014 年 10 月 27 日受付)

謝辞
著者は 2008 年 3 月に岩手大学大学院工学研究科建設環境工学専攻博士課程前期を修了した後，2008 年 4 月に独立行
政法人港湾空港技術研究所に採用され，2014 年 9 月に至る(b) 水表面に油を投入した直後まで行ってきた研究開発業務について，博士の学位論文とし

図- 8.3: 水表面の油の有無と風による水表面の波立ちの関係てとりまとめています．本稿は，著者の学位論文の書式を一(風速 9.3 m/s)

部変更し，港湾空港技術研究所資料としたものです．研究を通じて数多くの方々に御指導，御協力をいただきました．こ

Fay のモデルの評価を行うとともに，Fay のモデルにおけるこに記して感謝の意を表します．

拡がり係数を提案した．また，Fay のモデル及び実験結果に

私は課程博士ではなく，論文提出による博士学位の申請を基づいた新しい数値計算モデルを提案した．さらに，静水及行いました．申請に当たり，岩手大学工学部社会環境工学科

び定常流条件下における油の拡がりの実験結果を本数値計算

小笠原敏記准教授に推薦，審査していただき，大変お世話にモデルで予測できることを確認した．

なりました．

油の乱流拡散については，実海域における乱流拡散実験を

岩手大学 船崎健一工学部長，岩手大学工学部社会環境工

行い，海水面極近傍の水平乱流拡散係数は水中のそれよりも学科 海田輝之教授には，大変お忙しい中，審査を通じて有大きいことを示した．また，水平方向乱流拡散係数を導出す

益な助言をいただきました．

る手法を開発した．

岩手大学 堺茂樹学長には，大変お忙しい中，審査を通じ油拡散粒子モデルの有用性については，開発した油拡散粒て有益な助言をいただきました．また，岩手大学及び岩手大子モデルを用いて，油流出事故の再現計算及び予測計算を行学大学院博士前期課程在学中の頃から研究室の指導教官としい，観測結果と計算結果の比較を行っている．

て，大変多くのことを指導していただきました．

最後に，開発した数値計算法について考察を加えた．また，

みずほ情報総研株式会社 世古千博氏，山崎昇氏，荒木和博

即時の予報体制の構築状況について示した．さらに，開発し氏，石井隆氏には油拡散粒子モデルのプログラム開発，潮汐た数値計算法を応用して，油回収船の油回収操船支援を行う

流の数値計算，漂流油捕捉システムの開発等に関しまして，

漂流油捕捉システムを開発したことを示した．

補助及び指導していただきました．

以上，本論文では，海上に流出した油の移流及び拡散に関株式会社エコー ，長谷川巌氏，宮里一郎氏，小玉芳秀氏にする数値計算法の開発，油膜自身の特性による油拡散の新しは油膜自身の特性による油拡散に関する水槽実験及び解析を

い数値計算モデルの構築，海水面極近傍の水平乱流拡散係数

補助していただきました．

の導出方法の検討，数値計算法の有用性の検証を行い，研究

株式会社ソニック 高島勝美氏（現在は退職されています．
），

目標が達成された．しかしながら，本数値計算法のみで流出

磯部憲雄氏には，疑似油を用いた実海域における乱流拡散実油の挙動が全て予測できるわけではない．本数値計算モデル

験を補助していただきました．

44

一般財団法人日本気象協会 宇都宮好博氏，前田正裕氏，田聴会の発表内容について貴重なご意見をいただきました．

中秀雄氏には，潮汐流計算の実施及び潮汐流データベースの国土交通省国土技術政策総合研究所 井上岳主任研究官に構築を補助していただきました．

は，国総研資料の Latex の cls ファイルを提供していただく

東洋建設株式会社 嶋本孝氏，橋本光壽氏，大塚昌光氏，佐とともに，書き方について指導していただきました．

藤眞一氏，斉脇伸英氏には，油拡散の室内実験や実海域での最後に，本当に多くの皆様に大変お世話になりました．改

レーダー油膜検出システムの検証実験の補助をしていただき

めて深く感謝申し上げます．

ました．

参考文献

本稿の研究開発成果の一部は国土交通省港湾局及び国土交
通省北陸地方整備局の受託研究で行われたものです．特に元Akiyama, T., Mitsumori, S., Terada, T., and Koide, K. (1983): A

国土交通省港湾局技術企画課技術監理室 宮田正史氏，北陸Study of the Initial Stage of Oil Spreading on Calm Water

地方整備局 柿崎慶治氏，齋藤英治氏に大変お世話になりま

through Methods of both Instantaneous and Constant Rate ofした．

Oil Release, Journal of Chemical Engineering of Japan, Vol.

港湾空港技術研究所新技術研究開発領域 藤田勇領域長（油16, No. 6, pp. 502-507.

濁対策研究チームリーダー併任）には，油濁対策研究チーム

AL-Rabeh, A. H. (1994): Estimating Surface Oil Spill Transportのチームリーダーとして，研究に関する様々な面でお世話にdue to Wind in the Arabian Gulf, Ocean Engineering, Vol.

なりました．

21, pp. 461-465.

港湾空港技術研究所 吉江宗生研究計画官は，元油濁対策研究室長（チームリーダー）であり，油濁対策に関する研究

Allen, A. A. and Estes, J. E. (1970): Detection and Measurement

や周辺情報について，多くの指導をしていただきました．まof Oil Films, Santa Barbara Oil Symposium.

た，研究初期において資金面等の研究環境を整えていただき

ASCE Task Committee on Modeling of Oil Spills of the Water

ました．
港湾空港技術研究所 高橋重雄理事長，栗山善昭特別研究

Resources Engineering Division (1996): State-of-the-Art Re-

官，勝海務特別研究官，元港湾空港技術研究所の理事で，株view of Modeling Transport and Fate of Oil Spills, Journal of

式会社エコー 永井紀彦氏，元港湾空港技術研究所の研究主

Hydraulic Engineering, Vol. 122, pp. 594-609.

監で，横浜国立大学 中村由行教授には，研究内容の取りま

Blokker, P. C. (1964): Spreading and Evaporation of Petroleumとめに関して助言をしていただくなど，様々な場面でお世話

Products on Water, Proc. 4th Int. Harbour Conf., pp. 911-919.

になりました．
元港湾空港技術研究所で東京都市大学 田中陽二講師にはBuist, I., Potter, S., Mackay, D., and Charles, M. (1989): Labora-

調和解析の方法，潮汐流の数値計算等について指導していた

tory Studies on the Behavior and Cleanup of Waxy Crude Oil

だきました．

Spills, Proceedings of International Oil Spill Conference, pp.

国土交通省国土技術政策総合研究所 熊谷兼太郎主任研究105-113.

官には，風洞水槽実験を行う際に便宜を図っていただきま

Castanedo, S., Medina, R., Losada, I. J., Vidal, C., Mendez, F. J.,

した．

Osorio, A., Juanes, J. A., and Puente, A. (2006): The Pres-

港湾空港技術研究所新技術研究開発領域 白井一洋主任研究tige Oil Spill in Cantabria (Bay of Biscay). Part I:Operational官，田中敏成主任研究官，松本さゆり主任研究官，平林丈嗣

Forecasting System for Quick Response, Risk Assessment

研究官の皆様には機械，電気電子，コンピュータ関係といっand Protection of Natural Resources, Journal of Coastal Re-

た，私がよく知らないことを教えていただきました．

search, Vol. 22, No. 6, pp. 1474-1489.

港湾空港技術研究所 鈴木高二朗チームリーダー，本多和彦主任研究官，元港湾空港技術研究所で国土交通省 辰巳大Chao, X., Shankar, N. J., and Wang, S. S. Y. (2003): Development

介課長補佐には，研究を進めるに当たり，大変貴重なご意見

and Application of Oil Spill Model for Singapore Coastal Wa-

をいただき，ご指導をしていただきました．

ters, Journal of Hydraulic Engineering, Vol. 129, pp. 495-503.

港湾空港技術研究所 井上徹教チームリーダー，細川真也

Fay, J. A. (1969): The Spread of Oil on a Calm Sea, Plenum Press,

主任研究官には，本論文を取りまとめるに当たり，大変貴重pp. 53-64.

なご意見をいただきました．また，有川 太郎上席研究官，高川智博研究官，関克己研究官，茂木博匡専任研究員には，公45

Fay, J. A. (1971): Physical Processes in the Spread of Oil on a

Richardson, L. F. (1926): Atmospheric Diﬀusion Shown on a

Water Surface, Conference on Prevention and Control of Oil

Distance-Neighbor Graph, Proceedings of the Royal Society

Spills, pp. 463-467.

of London, Vol. 110, No. 756, pp. 709-737.

Guo, W. J. and Wang, Y. X. (2009): A Numerical Oil Spill Model

Shen, H. T. and Yapa, P. D. (1988): Oil Slick Transport in Rivers,

Based on a Hybrid Method, Marine Pollution Bulletin, Vol.

Journal of Hydraulic Engineering, Vol. 114, pp. 529-543.

58, pp. 726-734.

Smagorinsky, J. (1963): General Circulation Experiments with thePrimitive Equations, Monthly Weather Review, Vol. 91, No. 3,

Hoult, D. P. (1972): Oil Spreading on the Sea, Annual Review of

pp. 99-164.

Fluid Mechanics, Vol. 4, pp. 341-368.
Johansen, O. (1987): DOOSIM - A New Simulation Model for

Sotillo, M. G., Fanjul, E. A., Castanedo, S., Abascal, A. J., Menen-

Oil Spill Management, Proceedings of International Oil Spilldez, J., Emelianov, M., Olivella, R., Garcia-Ladona, E., Ruiz-

Conference, pp. 529-532.

Villarreal, M., Conde, J., Gomez, M., Conde, P., Gutierrez,
A. D., and Medina, R. (2008): Towards an Operational Sys-

Lee, B. (2008): Hebei Spilit 号汚染事故防除現状及び今後の計

tem for Oil-Spill Forecast over Spanish Waters: Initial Devel-

画，2008 年度韓国海洋科学技術協議会共同学術大会 (韓opments and Implementation Test, Marine Pollution Bulletin,国語)，第 56 巻, pp.10-20．

Vol. 56, pp. 686-703.

Lee, J. L., Lee, D. Y., and Kim, I. H. (2009): The Influence ofTkalich, P., Huda, M. K., and Gin, K. H. (2003): A Multiphase Oil

Waves on Oil Spill Behavior, Proceedings of the International

Spill Model, Journal of Hydraulic Research, Vol. 41, No. 2,

Conference on Oﬀshore Mechanics and Arctic Engineering,pp. 115-125.

Vol. 4, pp. 683-688.

Toba, Y., Kawamura, H., Yamashita, F., and Hanawa, K. (1984):

Lehr, W. J., Fraga, R. J., Belen, M. S., and Cekirge, H. M. (1984):

Structure of Horizontal Turbulence in the Japan Sea, Ocean

A New Technique to Estimate Initial Spill Size Using a Modi-

Hydrodynamics of the Japan and East China Seas, Vol. 39,

fied Fay-Type Spreading Formula, Marine Pollution Bulletin,

pp. 317-332. Elsevier.

Vol. 15, pp. 326-329.

Varlamov, S. M., Yoon, J., Hirose, N., Kawamura, H., and Shio-

Matsumoto, K., Takanezawa, T., and Ooe, M. (2000): Ocean Tide

hara, K. (1999): Simulation of the Oil Spill Processes in the

Models Developed by Assimilating TOPEX/POSEIDON Al-

Sea of Japan with Regional Ocean Circulation Model, Journal

timeter Data into Hydrodynamical Model: A Global Modelof Marine Science and Technology, Vol. 4, No. 3, pp. 94-107.

and a Regional Model around Japan, Journal of Oceanography, Vol. 56, pp. 567-581.

Wang, S., Shen, Y., Guo, Y., and Tang, J. (2008): ThreeDimensional Numerical Simulation for Transport of Oil Spills

Musk, S. (2012): Trends in Oil Spills from Tankers and ITOPF

in Seas, Ocean Engineering, Vol. 35, pp. 503-510.

Non-tanker Attended Incidents, Paper Presented at AMOP.

Wang, S. and Hwang, L. (1974): Numerical Simulation of Oil SlickOkubo, A. (1971): Oceanic Diﬀusion Diagrams, Deep-Sea Re-

Transport in Bays, Proceedings of the International Confer-search, Vol. 18, pp. 789-802.

ence on Coastal Engineering, Vol. 14, pp. 2227-2244.

Proctor, R., Flather, R. A., and Elliott, A. J. (1994): Modelling

Wang, W. D., Shen, Y. M., and Zheng, Y. H. (2005): Two-

Tides and Surface Drift in the Arabian Gulf -Application to

Dimensional Numerical Simulation for Transport and Fate of

the Gulf Oil Spill, Continental Shelf Research, Vol. 14, No. 5,

Oil Spills in Seas, Ocean Engineering, Vol. 32, pp. 1556-

pp. 531-545.

1571.

Reed, M. and Gundlach, E. (1989): A Coastal Zone Oil Spill

Wu, Z. and Wang, D. (2010): Simulation of the Oil Slick Move-

Model Development and Sensitivity Studies, Oil and Chemi-

ment in Tidal Waterways, Journal of Hydrodynamics, Vol. 22,

cal Pollution, Vol. 5, pp. 411-449.

pp. 96-102.

46

市橋理・赤松幸生 (1997): ダイヤモンドグレース号流出油の漂畑山健 (2012): 石油タンクの津波被害について，第 15 回消防

流状況，写真測量とリモートセンシング，第 36 巻, pp.2-4．

防災研究講演会資料, pp.85-104．

Egset，C. N.・Nost，E.・原裕英 (2007): X バンドレーダーを用藤田勇・八尋明彦・桑江朝比呂・吉江宗生 (2008): 韓国泰安

いた漂流油検出，日本マリンエンジニアリング学会誌，

沖油流出事故に関する現地調査報告，港湾空港技術研究第 42 巻，第 5 号, pp.8-13．

所資料，第 1186 号, pp.24．

海上災害防止センター (2000): 流出油事故対応防除マニュ

本多忠夫・光易恒 (1980): 水面に及ぼす風の作用に関する実アル．

験的研究，海岸工学講演会論文集，第 27 巻, pp.90-93．

海上保安庁 (1992): 日本沿岸潮汐調和定数表．

松崎義孝・藤田勇 (2012): 油拡散を考慮した流出油の数値計算法の開発，港湾空港技術研究所資料，第 1255 号, pp.26．

海上保安庁警備救難部環境防災課 (2014): 海洋汚染の現状（平成 25 年 1 月〜12 月），海上防災．

松本晃治 (2004): やさしい海洋潮汐モデリング，東京大学地震研究所研究集会地球の「流れ」を見る衛星重力ミッショ

佐尾和子 (1998): 重油汚染・明日のために-「ナホトカ」は日ン, pp.41-50．

本を変えられるか-，海洋工学研究所出版部．
榊山勉・阿部宣行・鹿島遼一 (1990): ポーラスモデルよる透

道田豊・田中潔・小松輝久・石神健二・中嶋理人 (2009): 海

過性構造物周辺の非線形波動解析，海岸工学論文集，第洋の収束・発散・拡散等による表層浮遊物の移動，海岸

37 巻, pp.554-558．

海洋研究，第 46 巻，第 2 号, pp.77-83．

坂本茂宏 (1998): ダイヤモンドグレース号原油流出事故への埜口英昭・早川典生・橋本英資・山崎宗広 (1979): 海上に流出対応，海上防災，第 96 号, pp.3-14．

した油の拡がりについて-瞬間放出による実験的研究-，
第 26 回海岸工学講演会論文集, pp.628-632．

白石哲也・藤田勇・松崎義孝 (2011): メキシコ湾油流出事故現地調査報告，港湾空港技術研究所資料，第 1239 号, pp.44．

柳哲雄・樋口明生 (1982): 現場海域の拡散係数の推定法に関

新保修・土市進・野田登・伊賀浩之・羽田広希 (2003): 大型する研究，海岸工学講演会論文集，第 29 巻, pp.580-583．

油浚渫船「白山」の油回収システムについて，海洋開発

横浜市港湾局 (1998): 東京湾原油流出事故による影響調査（横論文集，第 19 巻, pp.37-42．

浜港編）．

鈴木高二朗 (2013): 東京湾の海水交換と貧酸素化に及ぼす淡渡辺正孝・天野邦彦・石川裕二・田村正行・村上正吾・木幡

水流入と風の影響について，港湾空港技術研究所資料，

邦男 (1998): 東京湾におけるタンカー事故による原油流

第 1276 号, pp.97．

出解析，海岸工学論文集，第 45 巻, pp.926-930．

石油連盟 (1992): 流出油の拡散予測モデル調査報告書．
高橋照男・北村吉朗・中田賢三 (1979): 静止水面上での油膜の拡がり-油膜内の粘性抵抗の拡がり速度に及ぼす影響-，
化学工学論文集，第 5 巻，第 5 号, pp.526-531．
田中陽二 (2008): 最尤法による潮流の調和解析と最適な分潮の組み合わせの探索方法，海岸工学論文集，第 55 巻,

pp.381-385．
鶴谷広一・細川恭史・日比野忠史・三好英一 (1997): 東京湾ダイヤモンドグレース号油流出事故の浮遊・漂着状況の把握と水域への影響調査，港湾技研資料，第 882 号, pp.21．
富田孝史・柿沼太郎 (2005): 海水流動の 3 次元性を考慮した
高潮・津波数値シミュレータ STOC の開発と津波解析への適用，港湾空港技術研究所報告，第 44 巻，第 2 号,

pp.83-98．

47

記号表

Rr and：0 から 1 の一様乱数

A

：油層面積

Rsum

： Rm の和C：濃度s：強制流速ベクトル計算パラメータ

Cno：斥力モデル導出における係数sd：界面活性物質の水中への溶解度

Cda ，Cdw：それぞれ大気から海面，海面からt：時間

大気の抵抗係数ta：油膜拡散面積の解析における時間のシフト

Cs：Smagorinsky 定数

ti：油粒子 i の投入からの経過時間

Cw

：風係数

Ubi

：海岸線，構造物等境界における強制流速

DH：水平乱流拡散係数 (m2 /s)

Do：界面活性物質の水中での拡散係数

Ucit ,g：重力加速度

Uoit

h

：油層厚

Uc,

h1 ,h2：数値計算モデルの油層厚判別値

Uo

計算値及び観測値の平均値i：着目する油粒子

Udi：油粒子の乱流拡散速度ベクトル

in：実海域実験における着目する疑似油

Ui：着目する油粒子 i の速度ベクトルj：着目する油粒子 i の近傍粒子

Uoi：海流による油粒子の移流速度ベクトル

K

：積分定数

Us

：着目する油膜位置の油拡散速度ベクトル：それぞれ Fay(1969) の第 1，第 2，第 3

Usi：油膜自身の特性による油粒子の拡散速度

k 1 ,k 2 ,k 3

ベクトル

k 3m

潮汐流の計算値及び観測値：それぞれ潮汐流精度検証における潮汐流のベクトル段階における拡がり係数

k 1m ,k 2m ,：それぞれ潮汐流精度検証における it 時のu si j：それぞれ第 1，第 2，第 3 段階の拡がり：粒子 j による粒子 i の油膜自身の特性による油拡散速度モデル係数

u si j：粒子 j による粒子 i の油膜自身の特性にk br：強制流速ベクトルの係数

L：疑似油の拡散スケール

l：油層の拡散半径

Ut i：潮汐流による油粒子の移流速度ベクトルli：粒子群全体の重心と着目粒子 i の距離

Uwi：風による油粒子の移流速度ベクトル

L ma x：拡散スケールと水平乱流拡散係数のUw：風による表面流V：油層の体積VC：濃度 C の拡散速度

Vi

：粒子体積

よる油拡散速度ベクトル

関係が一様に成り立つ範囲

N：実海域実験結果の解析における疑似油の個数n：潮汐流計算における海底粗度計数

W：水面近傍の風速No：斥力モデルにおける近傍粒子数W10：海上 10 m での風速ベクトルno：ある着目粒子周りの粒子密度

W10i：粒子 i の位置における海上 10 m でのnt

：潮汐流精度検証に用いた，潮汐解析

風速ベクトル

Xi

データ個数：着目する油粒子 i の位置：近傍粒子検索半径x, y, z：方向r bi

：油粒子 i と各海岸位置との距離ベクトルx in , yin：実海域実験結果の解析における疑似油 inr b1 ,r b2：強制流速ベクトル Ubi の判定距離

Rdir , Rvel

：それぞれ潮汐流の流速及び

R

の位置

x,
y

：実海域実験結果の解析における疑似油の
位置の平均

Rvel

流向の相関係数ri j：粒子 i ， j 間の距離

δ

：油層下面の水の境界層厚さri j：粒子 i ， j 間の距離ベクトル

∆

：空間刻み

Rm

：風況観測地点 m と対象地点の距離の∆t：計算時間刻み∆x,∆y,∆z： x,y,z 方向空間格子∆ρ：水と油の密度差θ：潮汐流精度検証における流速観測結果と逆数

Rn1 , Rn2

：それぞれ独立した平均 0，分散 1 の
正規乱数

r oi

計算結果のなす角：油膜の重心と着目粒子 i の距離
ベクトル

48

µo

：油の粘性係数

νH , νV：それぞれ水平渦動粘性係数，鉛直渦動粘性

νw：水の動粘性係数

ρa ,ρo ,ρ w

：それぞれ大気，油，水の密度σ2：疑似油の位置の不偏分散

σn：正味の表面張力係数 (N/m)σoa , σow ,：それぞれ油-大気間，油-水間，

係数

σ wa
τ

水-大気間の表面張力係数 (N/m)
：接線応力

49

A

付録：潮汐流の予測方法

表- A1: 潮汐流の計算条件
格子数

ここでは，油拡散粒子モデルにおいて用いている潮汐流の予測方法について，用いている数値計算モデルの概要及び計算の境界条件と，計算結果からの潮汐流の再現方法について

座標系

説明する．

A1 潮汐流の数値計算モデル及び境界条件
潮汐流の数値計算モデルは富田・柿沼 (2005) の高潮・津波数

格子間隔

東京湾

750 × 930 × 5

伊勢湾

600 × 540 × 5

瀬戸内海1680 × 990 × 5

東京湾

IX 系

伊勢湾

VII 系

瀬戸内海

IV 系

∆x

270m

値シミュレータ（Storm surge and Tsunami simulator in Oceans

∆y

270m

and Coastal areas，通称 STOC）を用いている．

∆z

10m(1,2 層)

STOC の基礎方程式は，porous モデル（榊山他, 1990）を適15m(3 層)用した 3 次元の連続方程式；式 (A1) と，Reynolds 方程式；式

20m(4 層)(A2) である．なお，このモデルでは，水面位置を鉛直方向に150m(5 層)

積分した連続方程式；式 (A3) から算出するため，砕波のよ

時間間隔

∆t

3s

うに水面形が鉛直方向に多価になる場合には適用できない．

海底粗度係数n0.0026水平渦

νH

Smagorinsky (1963)∂γi ui
=0
∂x i

(A1)

動粘性係数
鉛直渦

γv

∂ui ∂γ j ui u jγv ∂p++ Ci = −∂t∂x j
ρ0 ∂x i())(∂u j
ρ w − ρ0
∂
∂ui
+
+ γv
gi +
γ j νe
ρ0
∂x j
∂x j
∂x i
γz

∂
∂η
+
∂t
∂x

∫ η
−h

γ x udz +

∂
∂y

∫ η
−h

γy vdz = 0

モデル

νV

1cm2 /s

動粘性係数

(A2)

海面高度計データを潮汐解析し，流体力学数値モデルに同化して構築されている．

(A3)
A2 潮汐流計算結果からの再現方法

ここで，x i はデカルト座標系（x，y，z），ui は x i 方向の流

潮汐流は異常潮汐を除けば周期的に変化しているので，あ速（u，v，w），γi は x i 方向の有効面透過率（γ x ，γy ，γz ，る地点における潮汐流の振幅及び位相（遅角）が既知の場合，

0 ≤ γi ≤ 1），γv は有効体積多孔率（メッシュ内を占める液それらから流速の時系列を求めることができる．そこで，あ相の体積率，0 ≤ γv ≤ 1）
，t は時間，Ci は Coriolis 項，ρ0 は

らかじめ STOC-ML を用いて 1 か月間の潮汐流計算を実施

基準密度， p は圧力， ρ w は流体密度，gi は各方向の重力加し，その結果から各格子における振幅及び位相のデータベー

速度（0，0，g ），νe は実効動粘性係数（動粘性係数と渦動スを構築した．潮汐流データベースの構築範囲は，8.3 節で粘性係数の和），η は水面変動，h は基準水面からの水深であ説明した東京湾，伊勢湾，瀬戸内海である (図-A1 参照)．境る．Ci は f 0 を Coriolis のパラメータとすると，Ci は x 及び界条件は前述の Matsumoto らのモデルで求めた水位を図-A1

y 方向でそれぞれ − f 0 v 及び f 0 u となる．渦動粘性係数は本の開境界で与えた．水深データは内閣府の南海トラフの巨大編の通り，定数で与えられている．なお，乱流モデル（SGS

地震モデル検討会の提供する 270m メッシュのデータを用い

渦動粘性係数）を使用することも可能である．

た．なお，水深が深くなると波速が大きくなり，∆t を大きくこれらを基礎方程式として，Staggerd メッシュで空間を離取る必要がある（CFL 条件）．そのため水深を 200m までとし散化し，leapfrog 法により時間発展させる．圧力について，連て，それ以上の水深を打ち切った．初期密度場は 1,026kg/m3

立 1 次方程式を解いて求める方法を STOC-IC といい，静水

で一様に与えた．初期流速場は 0m/s とした．その他の計算

圧近似を適応して求める方法を STOC-ML という．油拡散粒

条件は表-A1 の通りである．

子モデルでは STOC-ML を用いて潮汐流を計算している．

次に計算結果から主要 8 分潮（Q1、O1、P1、K1、N2、M2、潮位の境界入力条件は Matsumoto et al. (2000) の NAO.99Jb

S2、K2）の調和解析を行い，潮流に関する振幅，遅角を算出

モデルを用いた．NAO.99Jb モデルは日本周辺海域の任意のした．調和解析は海上保安庁 (1992) における日本沿岸潮汐

海洋上の点において主要 16 分潮 (M2，S2，K1，O1，N2，P1，

調和定数表の潮汐計算略説に記された考え方に基づき，最小K2，Q1，M1，J1，OO1，2N2，Mu2，Nu2，L2，T2) を予測する

二乗法により解析を行った．以下に調和解析の概要を示す．

モデルである．このモデルは約 5 年分の TOPEX/POSEIDON

50

表- A2: 主要 8 分潮の定数表（松本, 2004）
分潮名

角速度 Vt k

天文引数 χ t k

波数 nt k（deg/hour）

（deg)

（無次元）

備考

Q1

13.3986609

358.039995

1

主太陰惰率潮O113.9430356

126.4841785

1

主太陰日周潮P114.9589314

350.0224311

1

主太陽日周潮K115.0410686

9.9677947

1

日月合成日周潮N228.4397295

8.0077897

2

主太陰惰率潮M2

28.9841042

136.4519732

2

主太陰半日周潮S230.0000000

359.9902258

2

主太陽半日周潮K230.0821373

199.9355894

2

日月合成半日周潮

Uot (t) = U0 +∑At k cos(Vt k t + ut k ) +tk∑

Bt k sin(Vt k t + ut k )

tk

(A5)
式 (A5) の左辺と右辺の差の二乗和が最小となる At k ，Bt k を差の総和の偏微分にて求められる．
調和解析により求めた潮流に関する振幅，遅角から潮流を求める方法は松本 (2004) の方法に基づき行った．ある地点（緯度 ϕ，経度 λ （西経は＋、東経は−））の，時刻 t の潮流(b) 伊勢湾

Ut (ϕ, λ,t) は式 (A6) で表すことができる．

(a) 東京湾

Ut (ϕ, λ,t) =8∑

Ut k (ϕ, λ) cos[Vt k t + χ t k − κ t k (ϕ, λ) − nt k λ]

i=1

(A6)
ここで，添え字 t k は分潮の番号，Ut k (ϕ, λ) ，κ t k (ϕ, λ) はそれ
ぞれ調和解析で求めた流速の振幅と位相，t は 2000 年 1 月 1
日 00 時 00 分 00 秒（UTC）から起算した時間（hour），Vt k
は分潮の角速度， χ t k は天文引数と呼ばれる初期位相，nt k(c) 瀬戸内海は分潮の波数である．主要 8 分潮の角速度，天文引数及び波図- A1: 潮汐流計算範囲（地盤高が T.P. より低いところが白，数を表-A2 に示す．

陸地が黒）観測されたある時刻の潮流 Uot (t) は，式 (A4) で表すことができる．

Uot (t) = U0 +∑f t k Ut k cos(Vt k t + ut k − κ t k )

(A4)

tk

ここで，添え字 t k は各分潮を表し，U0 は潮汐残差流， f t k ，

Vt k ，ut k は天体の運動から求まる定数，Ut k ，κ t k は分潮の振幅と位相を表す。
式 (A4) の余弦関数を余弦関数と正弦関数の和によって表
すと，式 (A5) となる．

51

港湾空港技術研究所資料

No.1300

２０１５．３

編集兼発行人

独立行政法人港湾空港技術研究所

発

行

所

独立行政法人港湾空港技術研究所
横 須 賀 市 長 瀬 ３ 丁 目 １ 番 １ 号
TEL. 046(844)5040 URL. http://www.pari.go.jp/

印

刷

所

株 式 会 社大應

C （2015）by PARI
Copyright ○All rights reserved. No part of this book must be reproduced by any means without the written
permission of the President of PARI
この資料は、港湾空港技術研究所理事長の承認を得て刊行したものである。したがって、本報告書の全部または一部の転載、複写は港湾空港技術研究所理事長の文書による承認を得ずしてこれを行ってはならない。
E

A