学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

jh130024-NA16

ポストペタスケールシステムを目指した
二酸化炭素地中貯留シミュレーション技術の研究開発
山本肇（大成建設株式会社）概要 地球温暖化対策の二酸化炭素の地中貯留技術では、地中圧入後の CO2 の挙動を正確にシミュレートする技術が重要になる。本共同研究ではポストペタ・エクサスケール
システム上で大規模シミュレーションを実施するために改良版「TOUGH2-MP」の高度化（OpenMP/MPI ハイブリッド並列化，並列メッシュ生成・領域分割手法，最適化された並列反復法による疎行列ソルバー，並列リアルタイム可視化手法）を実施する。本
年度は現在, 最適化された並列前処理付き反復法による疎行列ソルバーを構築中である。
今後は解析コード TOUGH2MP の全体にわたって OpenMP/MPI ハイブリッド並列化を施すとともに、試解析を実施して並列リアルタイム可視化手法を適用する予定である。

1. 研究の目的と意義ットとする計算機は Fujitsu

1-1.

大）である。

研究の目的

PRIMEHPC FX10（東地球温暖化対策の一つとして，火力発電所などから排出される二酸化炭素ガスを，地下深部の地

1-2.

研究の意義

層中に圧入するという二酸化炭素の地中貯留が国二酸化炭素地中貯留は地球温暖化対策の一つと内外で認知されつつある。その実用化にあたってして有効な手段であるが，その実施にあたっては，は，社会的受容の観点からも，圧入後の二酸化炭

その 100 年～1000 年以上における長期的な安定性

素の地中挙動を正確にシミュレートする技術の確や副作用として生じる環境影響についても十分に立が重要な課題となる。

吟味されなければならない．

著者等は米国ローレンスバークレイ国立研究所

二酸化炭素の地中貯留シミュレーションでは，で開発された有限体積法に基づくオープンソース

CO2 が地中にトラップされるメカニズムを考慮する

多 相 流 体 シ ミ ュ レ ー タ TOUGH2-MP(Zhang et

ために，温度や圧力条件で変化する超臨界 CO2 の流

al.2008)を元に物理モデルの追加，疎行列ソルバ

体物性（密度，粘性，地下水への溶解度）や，地ーの改良を実施し，地球シミュレータ，T2K オープ

下水と CO2 の相互作用（水相－CO2 相間の成分移動

ンスパコン（東大）を使用して最大１千万メッシ

や二相流特性：相対浸透率，毛細管圧力の効果等）

ュ規模の大規模非定常シミュレーションを実施しに伴う強い非線形性を考慮する必要がある．加えてきた（例えば、Yamamoto et al. 2013）。改良版て，地中での長期的な鉱物変質や水質変化を予測「TOUGH2-MP」は Flat MPI による並列化を採用しするための地球化学計算との連成解析や，大量のており，領域分割には METIS を使用しているため，

流体圧入に伴う地盤中の応力や変形，破壊を評価

問題規模，計算機資源の規模に制約があった。本する地盤力学との連成解析が必要になっている．

共同研究ではポストペタスケール，エクサスケー

このようなシミュレーションは計算負荷が高いた

ルシステム上でより詳細な大規模シミュレーショ

め，工学的利用はもとより，科学研究においてもンを実施するために改良版「TOUGH2-MP」の高度化

相当のモデルの単純化が求められるのが現状であ（OpenMP/MPI ハイブリッド並列化，並列メッシュ

る。

生成・領域分割手法，最適化された並列反復法に

本共同研究を通じて大規模二酸化炭素地中貯留

よる疎行列ソルバー，並列リアルタイム可視化手シミュレーション技術を確立することによって，法）を実施し，大規模二酸化炭素地中貯留シミュ

ポストペタスケール，エクサスケールシステム上

レーション技術の確立を図るものである。ターゲでより詳細な大規模シミュレーションを実施する
1

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

ことが可能となり，二酸化炭素地中貯留の有効性

京大学情報基盤センターの共同研究として実施しや安全性を詳細に検討するための強力な高速計算ている。役割分担は以下の通り。

ツールが提供できる。例えば，地層の高度な不均 山本肇（大成建設）：課題代表者、TOUGH2-MP

質性や重力不安定性によるフィンガリングや対流

Hybrid 並列化，CO2 地中挙動シミュレーション実現象が CO2 の長期安定性に大きな影響を与えるこ施とが指摘されているが(例えば、Ennis-King and 中島研吾（東京大学）：課題副代表者、並列反復

Paterson, 2005)，非常に微細な格子解像度を要求

法による疎行列ソルバー最適化、並列メッシュ

されるため，必要とされる計算時間の問題から十生成・領域分割手法, 並列リアルタイム可視化

分な検討がなされていない。また，二酸化炭素の手法

地中貯留では大量の流体を地下に圧入すると，広 片桐 孝洋（東京大学）：並列反復法による疎行

域の地下水圧を上昇させる可能性が指摘されてお列ソルバー最適化

り，そのような周辺環境影響の解析を行う上でも 七井 慎一（大成建設）： OpenMP/MPI ハイブリッ

大規模なシミュレーションモデルを解く必要があ

ド並列化る。このように，高解像度・大規模モデルでの CO2

 増岡 健太郎（大成建設）： CO2 地中挙動シミュレの地中挙動解析は，これまでの計算技術では未解ーション、OpenMP/MPI ハイブリッド並列化

明な要因の科学的理解の促進や広域の環境影響評

 平塚 裕介（大成建設）： CO2 地中挙動シミュレー

価上で重要であり，二酸化炭素回収・貯留で生じション、並列リアルタイム可視化手法

る地中での現象の科学的理解や社会的受容性の確保の面など，二酸化炭素地中貯留の実現に向けた(2) 共同研究分野

大きなインパクトがある。

超大規模数値計算系応用分野

本共同研究で開発する高度化技術のうち，(3) 当公募型共同研究ならではという事項など

OpenMP/MPI ハイブリッド並列化，並列メッシュ生
成・領域分割手法，最適化された並列反復法によ本研究は CO2 の地中貯留のシミュレーションをる疎行列ソルバーによってよりポストペタスケー

行う地下水工学の技術者と計算機科学者の共同研

ル，エクサスケールシステム上での大規模シミュ究である。地下水工学分野での超並列スパコンのレーションを高速に実施することが可能となる。

潜在ニーズは、CO2 地中貯留だけでなく放射性廃棄

また並列リアルタイム可視化手法によって，大規物の地層処分などで見いだされる。これまで国内

模シミュレーション結果をディスク上に保存するで上記分野での超並列スパコンの利用事例が少な

ことなく効率的なポスト処理が可能となる。これ

い理由として、モデリングを担う地下水技術者とらの機能はオープンソースのライブラリとして公計算機科学者との技術交流が十分でない点が指摘

開することによって，二酸化炭素地中貯留の他にされる。本共同研究は、第一線の計算機科学者がも様々な大規模シミュレーションに適用すること

民間の地下水技術者と直接連携する事例形成にもが可能であり，ポストペタスケール，エクサスケ

なっており、超並列スパコンの産業利用拡大にもールシステム上での大規模シミュレーションの普貢献するものと期待できる。この点が当公募型共及に貢献するものである。

同研究ならではという特徴である。

2. 当拠点公募型共同研究として実施した意義

3. 研究成果の詳細と当初計画の達成状況(1) 共同研究を実施した拠点名および役割分担

3-1．研究成果の詳細について

本研究は、大成建設株式会社技術センターと東(1) 解析コード2

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

本研究では、地中での超臨界 CO2 の挙動を数値F    u   k

解析するため、多成分・多相系の地下流体シミュ

k r  

レータ TOUGH2（Pruess et al., 1999）の並列化

P    g (5)

バージョンである TOUGH2-MP（Zhang et al., 2008）

ここに、μ:相の粘性係数(kg/ms)、u：流速、k：絶を用いている。二酸化炭素は地下深部約 1km に圧対浸透率(m2)、kr：相の相対浸透率、P：相の圧力入するのが安全面、効率面において最適と考えら(MPa)である。熱フラックスは、熱伝導と熱対流をれているが、その温度・圧力条件下において二酸考えて、F h  T   h  F化炭素は超臨界状態となる。超臨界状態の二酸化(6)

炭素の密度や粘性は、温度・圧力に対して敏感に変化する。また、水と二酸化炭素の相互作用（溶ここに、：熱伝導率(W/m℃)、T：温度(℃)、h：

解、混相流特性）も考慮する必要がある。今回のβ相の比エンタルピー（J/kg）である。

計算では、ECO2N という EOS モジュールを用い、(1)式を、積分差分法により時間と空間について

非等温系での３成分（H2O, NaCl, CO2）２相（水相、離散化すると、次のような残差形式の非線形代数

CO2 相）系流体を取り扱う。

方程式が得られる。

閉境界n に囲まれた体積 Vn に対して積分されR n (x t 1 )  M n (x t 1 )  M n (x t )

n



Vn
d
κ
q

Γ
d

n

Vn

  F n
Γ



κ

Vn
d
κ
M
Vn
t
d d

たエネルギー・質量保存則は次式で表される。

(1)

ここに、M：成分の単位体積当たりの質量(kg/m3)
またはエネルギー (J/m )、F ：成分のフラックス

圧力、飽和度など）、Rn：グリッドブロック n に(kg/m2s または J/m2s)、q：成分のシンク／ソース

おける成分に関する残差、Vn：ブロック n の体積、項（kg/m3s または J/m3s）、n は部分面要素 dn の法

qn：質量またはエネルギーのシンク／ソース、Δ

線ベクトルである。左辺の Mk は次のように書ける。
質量： M    S    X

熱： M  (1   )  R C R T  h(7)

ここに、ベクトル xt：時間 t における解（温度、



3

t 
t 1

 , t 1 
 A nm Fnm (x )  Vn q n   0


Vn m

t：時間ステップ、Anm：ブロック n と m が接する

(2)

面積、Fnm：ブロック n と m 間の平均フラックス、である。そして、(7)式をニュートン・ラプソン法

 S    U  (3)により解く。

ここに、：空隙率、S：相の飽和度、：相の

密度（kg/m ）、X：相における成分の質量分率、3i

R n , t 1
x i

( x i ,p 1  x i ,p )  R n , t 1 ( x i ,p )

(8)

p

3

ρR：岩石粒子の密度（kg/m ）、T：温度（℃）、CR：

ここに、xi，p+1 は i 番目の変数の p 回目の反復計算

岩石の比熱（J/kg℃）、U：相の内部エネルギー

ステップでの値である。

(J/kg)、上付き添字の h は熱を表す。
移流による質量フラックスは、全ての相のフラ(2) 計算性能測定

ックスを合計して求められる。
Fadv

  X  F

現時点までに得られた FX10 上での TOUGH2-MP の(4)

計算性能を示す。Yamamoto et al.（2013）は，図

-1 に示すような節点数や自由度の異なる４つのモ

β相の質量フラックスは、次式の多成分に拡張

デルについて計算性能を測定している。そのうち、したダルシー則により計算する。

(1)のモデルについて、ES2 と FX10（最適化前）での計算時間を比較した結果を図-2 に示す。これを3

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

見ると、FX10 では 2048 コアまで良好なスケーラ
ビリティが得られている。しかし，現状では ES2(3) 重力対流を考慮した CO2 長期的安定性のシミの 16 コアは FX10 の 512 コアよりも高速である。

ュレーション

ES2 と FX10 では１コア当たりのピーク性能が ES2

地下に圧入された超臨界状態の CO2 は、密度がは FX10 の約 7 倍高速（ES2:102.4GFlops、FX10：

水よりもやや小さい（比重 0.5～0.8）ため、浮力14.78Flops）である点に加え、１コア当たりの実によって地表へ上昇するポテンシャルを有してい効メモリバンク幅も ES: 約 200GB/s/コアに対しる。そこで、CCS では、貯留層を覆うキャップロて、FX10：約 5GB/s/コアと約 40 倍違う点に注意

ックと呼ばれる低浸透性の地層により、CO2 の上昇が必要である。

を抑え込むことが必要である。一方、CO2 は周辺の

今後は、OpenMP/MPI ハイブリッド並列化など現地下水に溶け込むが、CO2 が溶け込んだ水の比重は

在実施中の最適化を施した後、再度の計算性能測若干増加するため（地下水の塩分濃度によるが、定を行う予定である。

0.1%～1%程度）、負の浮力を有する安定性の高い貯(1) Tokyo Bay Model–Impact on regional groundwater flow –留状態になると考えられている。この点に関し、(2) CLASTIC Model–Geothermal Production–

上側の流体が下側の流体より重いときに生じる重Injection
Well

力不安定性（レイリー・テーラー不安定性,
Rayleigh- Taylor instability）が CO2 の地下水への溶解を促進することが指摘されている

Audigane et al. (2011)
Yamamoto et al. (2009)

10km(Ennis-King and Paterson, 2005)。

25.6 million dof

30 million dof (10 M grids, 3components)(6.4 M grids, 3components+thermal)(3) DDC (Diffusion‐Dissolution‐Convection)–Enhancement of CO2 Dissolution –(4) SPE 10 Model–Highly Heterogeneous Reservoir –Producer

Injection Well

Local‐scale
Model

Original
Reservior Model
Christie and Blunt (2001)
Qi et al. (2009)

CO2 Stored (tons)

Reservoir

Heterogeneous (original)
Homogeneous (averaged)

6E+6

Injector

Native Groundwater (Brine)度の格子を用いる必要があることが知られている。

7E+6

Caprock (Low permeable seal)
Supercritical
CO2

重力不安定性を正確に数値計算するには高解像

5E+6

例えば、線形安定性解析の理論によれば、比較的

4E+6
3E+6
2E+6
1E+6

浸透率の高い貯留層（絶対浸透率 1darcy）では 1cm

0E+0
0

5

10

15

20

Time (years)

Injector

Injector

オーダー以上の格子間隔が必要とされる Pruess

Calc. Timeand Zhang(2008）は、数ｍ程度の局所モデルで密X 300

度対流による CO2 長期安定性への効果を検討してHeterogeneous(No upscaling)

Pruess and Zhang (2008)

Homogeneous

12 million dof (4.1 M grids, 3components)

3.3 million dof (1.1 M grids, 3components)

いる。しかし、超並列計算機を用いることで、局

図-1 計算性能の測定に用いたモデル一覧

所モデルをスケールアップし、重力対流が CO2 の（Yamamoto et al. 2013）

長期安定性に及ぼす効果を貯留層全体のスケール

1.E+5で計算できる可能性がある。そこで、実際の CO2

Calculation Time (sec)

貯留層をモデル化した３次元のフィールドモデル
を用いて、CO2 を溶解した地下水の重力対流をシミ
1.E+4

FX10
(Scalar Machine)ュレートすることを試みた（山本ほか. 2013）。モ
デル化範囲は CO2 圧入井から半径約 1.2km 範囲の厚さ約 80m の CO2 貯留層であり、CO2 が到達する圧

1.E+3

入井近傍の領域の格子サイズを約 50cm に詳細化
The Earth Simulator2(Vector Machine)

1.E+2

した。CO2 の圧入量 10 万 t/年、圧入期間 1 年とし、貯留層の有効孔隙率と絶対浸透率は各々20%、

10

100

1000

10000

100mD（鉛直、水平方向とも）とし、相対浸透率曲

Number of Processors

図-2 計算性能の測定結果

線には Corey の式を用い、毛管力の効果は無視し4

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

た。

③最適化された並列前処理付き反復法による疎行計算結果の一例として、図-3(a)に、圧入から

列ソルバー

200 年後における CO2 貯留層内での CO2 飽和度なら④並列リアルタイム可視化手法

びに地下水中の CO2 濃度を示す。同図中(b)に、地

改良版「TOUGH2-MP」は既に MPI によって並列化

下水中の CO2 濃度の等値面を示す。貯留層上部にされているため，その並列データ構造を元に，①超臨界状態の CO2 が滞留した状態から、地下水中

では疎行列ソルバー以外の要素積分，係数行列生に CO2 が拡散しながら溶解範囲を広げ、拡散域が成，保存量計算，非線形処理などの諸機能を OpenMP

ある厚さに達すると重力不安定性による地下水対によって並列化するとともに Fujitsu FX10 向け最流が発生しフィンガリングが成長する様子がシミ適化を実施する。②では、並列メッシュ生成にあュレートされ、貯留層上部に貯留された超臨界 CO2

たっては分散メッシュを細分化する機能を開発し，の地下水への溶解が促進される様子が確認された。 これに ParMETIS，PT-SCOTCH 等の既存の並列グラ
今後は、さらに格子を詳細するとともに、実デー

フ分割ライブラリを組み合わせて領域分割機能をタに基づく地層モデルのキャリブレーションを行

開発，実装する。③では OpenMP/MPI ハイブリッド

う予定である。

並列化が施された ILU 前処理付き反復法による疎行列ソルバーを開発する。TOUGH2-MP は有限体積法によって離散化されており，疎行列を係数とす

Injection Well

る大規模連立一次方程式の求解が計算コストを要するプロセスである。TOUGH2-MP は多相流体シミ
ュレータであるため，疎行列ソルバーに対しては多自由度であることの特性を生かしてブロック化

Injection Wellを適用し Fujitsu FX10 向け最適化を実施する。
④についてはボクセル格子型背景格子を利用し，ハイブリッド並列化が施された並列可視化手法を(a) CO2 貯留層内での CO2 飽和度（上図）ならびに地下水中の溶存 CO2 の質量分率（下図）（圧入井を通る鉛直断面）実装する。
2013 年度は，②を除く各項目について研究開発を実施する計画であった。これまでに、③最適化
された並列前処理付き反復法による疎行列ソルバ
ーの構築を進めているところであるが，当初計画
より時間を要している。最適化前の計算性能の測定結果からは，計算速度向上の余地が大きいと考えられ，引き続き作業を継続する予定である。

(b) 地下水中の CO2 質量分率の等値面(緑色：0.02）
図-3 貯留層スケールの重力対流のシミュレーシ
ョン結果（約 200 年後）

4. 今後の展望
今後は、現在作成中の③最適化された並列前処
理付き反復法による疎行列ソルバーと並行して，
解 析 コ ー ド TOUGH2MP の 全 体 に わ た っ て ①

3-2. 当初計画の達成状況についてOpenMP/MPI ハイブリッド並列化を施し、試解析を本共同研究の実施項目は以下の４項目である。

実施して並列リアルタイム可視化手法を適用する①OpenMP/MPI ハイブリッド並列化予定である。また、②を開発，実装することによ②並列メッシュ生成・領域分割手法
5

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月って更なる大規模問題へ適用を図る。

廣川雄一, 緒方隆盛, 中島研吾，二酸化炭素地

東大情報基盤センターで実施している公募型研究

中貯留に関する大規模シミュレーション技術「大規模 HPC チャレンジ」に応募し，Fujitsu FX10の開発, 平成 24 年度地球シミュレータ産業戦

4,000 ノード以上を使用して諸機能の検証を実施

略利用プログラム成果報告書, 2013する。

5. 研究成果リスト(1) 学術論文
Yamamoto, H., S. Nanai, K. Zhang, P. Audigane, C.

参考文献

Chiaberge, R. Ogata, N. Nishikawa, Y. Hirokawa, S.

1) Ennis-King JP, Paterson L. Role of convective
mixing in the long-term storage of carbon dioxide

Shingu, and K. Nakajima, 2013, Numericalin deep saline formations. SPE J, 10(3):349–56.

Simulation of Long-term Fate of CO2 Stored in

SPE-84344-PA, 2005.

Deep Reservoir Rocks on Massively Parallel Vector
Supercomputer,2) Pruess, K., C. Oldenburg, and G. Moridis,
Lawrence

Berkeley

Notes

in

Computer

Science Vol. 7851, pp.80-92, 2013.

TOUGH2 User’s Guide, Version 2.0, Rep.
LBNL-43134,

LectureMukhopadhyay, S., C. Doughty, D. Bacon, J. Li, L.

National

Wei, H. Yamamoto, S. Gasda, S.A. Hosseini, JP

Laboratory, Berkeley, CA, 1999.
3) Pruess, K., K. Zhang, Numerical modeling studies

Nicot, J.T. Birkholzer, Comparison of Selectedof the dissolution–diffusion–convection process

Flow Models of the S-3 Site in the Sim-SEQ

during CO2 storage in saline aquifers. Technical

Project, Transport inReport LBNL-1243E, Lawrence Berkeley

Porous Media (submitted)

Yamamoto, H., S. Nanai, K. Zhang, P. Audigane, C.
Chiaberge, R. Ogata, N. Nishikawa, Y. Hirokawa, S.

National Laboratory, California, 2008.

Shingu, and K. Nakajima, High-performance4) Yamamoto, H., Zhang, K., Karasaki, K., Marui, A.,
Numerical

supercomputing as a risk evaluation tool forinvestigation for the impact of CO2 geologic

geologic carbon dioxide storages, Energy Procedia

sequestration on regional groundwater flow,

37, pp. 3997-4005, 2013.

Uehara,H.,

Nishikawa,

N.,

International Journal of Greenhouse Gas Control,

Mukhopadhyay, S., Z. Hou, L. Gosink, D. Bacon, C.

Vol.3, pp.586-599, 2009.

Doughty, J.J. Li, L. Wei, S. Gasda, G. Bacci, R.

5) Yamamoto, H., S. Nanai, K. Zhang, P. Audigane,

Govindan, J.-Q. Shi, H. Yamamoto, R. Ramanathan,C. Chiaberge, R. Ogata, N. Nishikawa, Y.

JP Nicot, S.A. Hosseini J.T. Birkholzer, A.

Hirokawa,

Nakajima,

Bonneville, Model Comparison and Uncertainty

High-performance supercomputing as a risk

Quantification for Geologic Carbon Storage: The

evaluation tool for geologic carbon dioxide

Sim-SEQ Initiative, Energy Procedia 37, pp.

storages, Energy Procedia 37, pp. 3997-4005,

3867-3874,2013.

S.

Shingu,

and

K.

2013.

(2) 国際会議プロシーディングス

6) Zhang, K., Y.S. Wu, and K. Pruess, User’s Guide

Yamamoto,H., K. Nakajima, S. Nanai, K. Masuoka, Y.

for TOUGH2-MP -. A Massively Parallel Version

Hiratsuka, T. Katagiri, K. Zhang, Post-peta scale highof the TOUGH2 Code, LBNL-315E, Lawrence

performance computing for coupled multi-physics

Berkeley National Laboratory Report, 2008.

modeling

7) 山本 肇，七井慎一,

of

CO2

geologic

sequestration,

12th

International Conference on Greenhouse Gas Control

Keni Zhang, 西川憲明,
6

学際大規模情報基盤共同利用・共同研究拠点 平成 25 年度共同研究 最終報告書 2014 年 5 月

Technologies (accepted)(3) 国際会議発表
山本 肇，CCS (carbon dioxide capture & storage) -A
viable option for mitigating the contribution offossil fuel emissions to global warming -),上海同済大学

機 械 ・ エ ネ ル ギ ー 工 学 科 、 ME

(Mechanical Engineering) Seminar, 2013 年 6 月
17 日
山本 肇，High-performance computing for predicting
the long-term safety and permanence of CO2 stored
in deep reservoir rocks，北京師範大学 水科学学
院 セミナー，2013 年 6 月 21 日
(4) 国内会議発表
山本 肇, Mukhopadhyay, S.，二酸化炭素地中貯留モ
デリングの国際比較研究 Sim-SEQ―事前予測
シミュレーション―, 第 68 回土木学会年次学
術講演会，土木学会，2013.
山本 肇, 中島研吾，超並列コンピュータ上における地質媒体中の多成分・多相流体解析の計算性
能－二酸化炭素地中貯留を対象とした事例－,
第 69 回土木学会年次学術講演会，土木学会，(投稿中).
(5) その他（特許，プレス発表，著書等）
山本 肇，七井慎一,

Keni Zhang, 西川憲明, 廣川雄一, 緒方隆盛, 中島研吾，二酸化炭素地中貯留に関する大規模シミュレーション技術の開発, 平成 24 年度地球シミュレータ産業戦略利用プログラム成果報告書, 2013.

7