Silviase/jgov_v1.4_full_ann huggingface dataset:
国立環境研究所の3次元大気輸送モデルに、 $ CH_{4} $ 濃度に加えて、炭素安定同位体比( $ \delta^{13}C $ )、水素安定同位体比( $ \delta D $ )、放射性炭素同位体比( $ ^{14}C $ )に関するプロセスを組み込んだモデルフレームを開発した。既存のメタン収支シナリオ（表(2)-1）ならびに放出源の同位体比、反応時の分別係数を用いた数値実験を行った。 $ CH_{4} $ 吸収については、土壌による酸化30Tg/yr、OHラジカルとの化学反応を考慮し、動的同位体分別効果(KIE； $ k_{13}/k_{12} $ )は表(2)-2とし、 $ ^{14}C $ のKIEは $ \delta^{13}C $ の2乗と仮定した。成層圈では、OHとの化学反応の他、KIEの大きいC1、O( $ ^{1}D $ )との化学反応が重要となるが、本研究ではまず対流圈を対象とした数値実験を目的としたため、C1、O( $ ^{1}D $ )との反応の組み込みは未着手である。 $ ^{14}C $ の放出源については、加压水型(PWR)並びに沸騰水型(BWR)の軽水炉による1960年代以降の放出が無視できないため、代表的な軽水炉について公開されている総出力[W/yr]を1988年の観測値 $ ^{12)} $ にスケーリングして放出量を推定した。

表(2)-1 本研究で用いたメタン放出量とその同位体比

放出源	放出量[Tg-CH4/yr] *1	δ13C[‰] 5)6)	δD[‰] 7)	14C[pM] 8)
湿地(泥炭地、沼地、ツンドラ)	167	-60	-370	116
シロアリ	22	-57	-390	124
バイオマス燃焼	22	-25	-225	130
水田耕作	43	-64	-320	112
家畜	105	-62	-330	120
埋め立て地	69	-55	-310	120
バイオ燃料	18	-25	-225	130
化石燃料(炭坑、天然ガス探掘・輸送、その他),工業生産	113	-40	-162	0
合計放出量、または加重平均値	559	-53.4	-299	96.8
14C放出(軽水炉)	928 [Ci/yr] (2000年の値)

*1: Patra et al. (2009) $ ^{3)} $ によるインベントリを改変。Patra et al.は、人為起源放出として the Emission Database for Global Atmospheric Research (EDGAR) v3.1、自然起源放出として NASA/GISSのインベントリデータベース $ ^{4)} $ を用いている。
表(2)-2 本研究で用いたメタン消滅の動的同位体分別効果(KIE; $ k_{13}/k_{12} $ )

消减源	KIE(13C)	KIE(14C)	KIE(D)
土壤	0.978 9)	KIE(13C)^2	1.081 10)
OH	1.0039 11)	KIE(13C)^2	1.294 11)

図(2)-7に2000年1月及び7月のモデル最下層における計算結果を示す。年間を通して、$\mathrm{CH}_{4}$ 濃度は南半球より北半球が高濃度となり、これは、$\mathrm{CH}_{4}$ の放出源が主として陸上に存在するため、陸面積の多い北半球で$\mathrm{CH}_{4}$ 放出量が多いことに起因するためである。北半球ではメタン放出源の多い陸域内部で経度方向の不均一が見られ、陸域内部で海洋よりも高濃度となっている。南半球のメタン濃度はアフリカ・南米等の陸域を除いてほぼ一様である。南北濃度勾配は、1月よりも7月に大きくなる。一方、$\delta^{13}\mathrm{C}$ およびδD分布は、北半球で低く南半球で高い緯度分布を示し、濃度の緯度分布とは対照的である。これは、北半球陸上にδ¹³C値の低い微生物起源の放出源が多く存在すること、及び、北半球地表から放出された$\mathrm{CH}_{4}$ が南半球へ輸送されるに従い大気中でOHラジカルとの反応による動的同位体分別（表(2)-2）により同位体的に重くなることにより、さらに同位体的な南北差が大きくなると解釈できる。¹⁴Cは、北半球で若干高い値を示す。これは、モデル内で化石燃料起源の¹⁴Cを含まない$\mathrm{CH}_{4}$ 放出量よりも、北半球に多く位置する軽水炉からの¹⁴C放出量が上回っているためであり、観測結果¹²と整合する。地理的分布では、韓国・日本、ヨーロッパ、アメリカ等の軽水炉の位置に対応して¹⁴Cの高い領域が存在し、化石燃料起源の$\mathrm{CH}_{4}$ の低い¹⁴Cの影響は、アラビア半島、シベリア域等の軽水炉から地理的に離れた領域で認められた。シベリアでは、冬季に低い¹⁴Cをともなう$\mathrm{CH}_{4}$ 高濃度の領域が出現している。これは、冬季の大陸内部において発達するシベリア高気圧により大気の鉛直混合が妨げられ、化石燃料起源の$\mathrm{CH}_{4}$ がモデル地表面に蓄積した結果であると解釈できる。夏季のシベリアでは、$\mathrm{CH}_{4}$ は高濃度であるが、δ¹³C、δDともに低い地域が見られ、放出源シナリオの湿地からの放出を反映している。夏季は鉛直混合が盛んであるため、放出される$\mathrm{CH}_{4}$ はモデル下層部へ蓄積されず上空へ輸送されると考えられる。このように、$\mathrm{CH}_{4}$ 濃度だけでは判別できない、放出源毎の同位体的特徴がモデルにより確認された。
(a)

CH₄ [ppbv] JAN 2000

(b)

CH₄ [ppbv] JUL 2000

(c)

(d)

(e)

(f)

(g)

(h)

図(2)-7．2000年1月及び7月におけるモデル最下層での(a)(b)CH_{4}濃度[ppbv]（コンター最大値2050ppbv、最小値1750ppbv、コンター幅30ppbv）、(c)(d) $ \delta^{13}C[\%o] $ （同-47‰、-49‰、0.2‰）(e)(f) $ \delta D[\%o] $ （同-87‰、-120‰、3‰）、(g)(h) $ ^{14}C[pM] $ （同135pM、122pM、1pM）。
年平均値緯度分布について、アメリカNOAA/ESRL/GMDによるバックグラウンドを主とした観測値 $ ^{13)(14)} $ とモデル計算値を比較した結果（図(2)-8）、陸域のMt. Waliguan, China（WLG）、Tae-ahn Peninsula, Koria（TAP）では、ローカルの放出源の影響を受けるため、モデルと観測値の差が見受けられるものの、バックグラウンドサイトにおいては、モデルは観測されるCH $ _{4} $ 濃度の南北勾配を再現することができた。一方で、 $ \delta^{13}C $ に関しては、40度以北の $ \delta^{13}C $ 値をモデル計算値が過大評価する結果となり、観測される南北勾配を再現するには、CH $ _{4} $ 放出源の多い同緯度帯でのシナリオの同位体比の調整が必要であることが分かった。 $ \delta D $ に関しては同NOAAデータベースには比較しうる観測値は登録されていないが、過去の太平洋上での観測事例 $ ^{8)} $ によると、南北両半球の差は、北半球で約10%低いとされており、本モデルで得られた南北差はこれにほぼ一致する。

また、シベリア・ノボシビルスク上空で採取された試料をサブテーマ（1）が分析したδ¹³C値、δD値と比較した結果でも、モデル計算値と観測値との差違が見いだされ（図(1)-25参照）、同位体を用いることによりCH₄放出量の高精度化し得る可能性が示唆された。

図(2)-8 アメリカNOAA/ESRL/GMD観測サイトにおける観測値とモデル計算値の年平均値（SPOからの偏差）。（上）CH_{4}濃度、（中）δ^{13}C値、（下）δD値（観測値が無いため、モデル計算値のみ）。

次に、観測とモデル計算値をディジタルフィルターを用いてフィッティングすることにより抽出した周期的季節変化成分を解析した結果を図(2)-9に示す。その結果、バックグラウンドサイトにおいて観測されるCH_{4}濃度と逆位相となるδ^{13}Cの季節変化、及び1ヶ月の位相のずれを、数値モデルで概ね再現することができた。太平洋のMauna Loa (MLO)では観測とモデルの一致が最も良く、南極SPOではモデルの位相が1ヶ月ほど早くなっており、これはモデル内の南半球への輸送の誤差も考えられる。北極域(BRW)での季節変動は、観測される濃度とδ^{13}Cの振幅及び位相のずれは再現
されたものの、モデルの位相が観測値よりも全体に2ヶ月早いことが分かった。これらの差異の原因を解析することにより、シナリオの改善が期待される。濃度- $ \delta^{13}C $ のphase diagramを作成した結果、図(2)-10に示すように、南極域(SPO)の位相変動はOHラジカルとの反応における動的同位体効果でほぼ説明がつく一方、北極域(BRW)では放出源の影響が大きく、これは特に夏季に顕著であり、陸域のCH $ _{4} $ 放出源から輸送されるCH $ _{4} $ の影響が示唆された。

図(2)-9 アメリカNOAA/ESRL/GMD観測サイトBarrow, Alaska (BRW), Mauna Loa, Hawaii (MLO), South Pole, Antarctica(SPO)における、CH_{4}濃度、δ^{13}Cの周期的季節変化成分の観測値（緑線）、モデル計算値（黒線）δDはモデル計算値のみ。

図(2)-10 CH_{4}濃度変化率と $ \delta^{13}C $ 変化の位相図。（左）BRW（右）SPO。図中の数字は月を示す。鎖線の直線はOHラジカルとの反応時のKIEのみを仮定した際の勾配を示す。
図(2)-11 観測値及びモデル計算値から得られたシベリアSurgut、Novosibirsk上空(1km、3km、7km)での周期的季節変化成分。

本サブテーマでは、シベリア西部に位置するスルグート（Surgut；SUR）、ノボシビルスク（Novosibirsk；NOV）において航空機観測を行っている。モデル計算結果では、この地域は夏季に湿地(Bogタイプ)からの放出の影響が大きく、また冬季においては、湿地からの放出は無いものの、ヨーロッパを起源とする人為起源メタンが流入している計算結果が得られた。従って、観測値は夏季の湿地からの放出を推定する束縛条件として有益であることが期待される。図(2)-11に本サブテーマによるSurgut及びNovosibirsk上空でのCH_{4}濃度観測値から得られた周期的季節変化成分をモデル計算結果と伴に示す。湿地帯に位置するSUR 1kmでは、観測値は夏季と冬季に最大となる年に2回の極大を示す。これは、上述したように、冬季には、化石燃料消費等により放出されたCH_{4}が鉛直混合の妨げにより下層に蓄積すること、夏季は、OHラジカルからの反応による消滅量よりも多いCH_{4}がシベリア域から放出されているとモデルでは説明できる。実際に同緯度帯のバックグラウンドサイトでは、夏季に極小を示す観測事実、並びにモデル計算結果と対照的である。しかしながら、モデル計算では1月と8月に極大を示すものの、観測値に対してモデル計算値の振幅が過大となっている。原因としては、モデルで使用しているシナリオの湿地からの夏季の放出が過大評価されている可能性がある。モデル内の湿地起源CH_{4}フラックスは、西シベリア域で約9Tg/yrである。一方、西シベリア湿地からのCH_{4}放出量に関して、近年の現地フラックス測定を基にしたボトムアップインベントリ $ ^{15)} $ では、同地域の湿地起源CH_{4}は3.2 Tg/yrと推定されている。しかし、西シベリア域の放出量を3.2 Tg/yrとしたモデル計算では、夏季の極大値がモデルで出現せず、同域湿地起源CH_{4}放出量を6 Tg/yrとした計算では夏の極大が出現することが分かった（図な
し）。すなわち、モデルで使用する湿地起源のCH_{4}フラックスの違いによって、西シベリア上空で観測される夏季のCH_{4}濃度のピークの再現性に違いが見られることが示され、SUR上空での観測は西シベリア湿地起源CH_{4}の束縛条件となり得ることが分かった。3km、7kmの高度では、下層の空気塊が上空へ輸送されるにつれ振幅が減衰しており、これは観測値でもモデル計算値でも同じ傾向であり、特にSURの7kmにおけるモデルの季節変化は観測値と良い一致を示した。一方、西シベリア低地の南端の森林地帯に位置するNOVにおいては、振幅は小さいものの、SURと同様に2つの極大が観測された。観測されるこの夏季の極大も湿地からの影響と思われるが、モデル計算値では、極大は秋季の1つのみであり、モデルで使用している湿地からの放出量及び位相が現実と異なっていることが示唆された。

## 5 ．本研究により得られた成果

## （1） 科学的意義

全球の温室効果ガス観測網における空白域の1つであるシベリア地域において、高精度の $ CH_{4} $ 濃度と $ N_{2}O $ 濃度のデータが継続して得られている。特に世界でも希な上空データの長期レコードであることも意義が深い。

CH_{4}濃度は2006年以降全球的に増加が確認されているが、その原因は未だ明らかでない。シベリアの湿地からのCH_{4}放出も原因の1つであるとする説もあるが、ノボシビルスクの濃度増加が全球での変動に先んじていないことやスルグートの低高度での増加が明確でないことからから、シベリア湿地からのCH_{4}放出は有力ではないと言えることを初めて明らかにしたことは意義がある。

$ N_{2}O $ は大気中での寿命が長いために時間変化や空間分布を把握しにくい物質であるが、高精度の分析技術と長期的なデータの蓄積によってその変動の実態を明らかにすることができた。特に $ N_{2}O $ 濃度の年平均値に有意な鉛直勾配があり、自由対流圈に比べて境界層内の濃度が高いことを観測結果の解析によって示したことにより、シベリアの地表面が $ N_{2}O $ の放出源であることが明らかになった。この結果は全球的な $ N_{2}O $ の放出源・吸収源分布を明らかにするにあたって、貴重な情報となることが期待される。

CH_{4}のアイソトポマーを組み込んだ大気輸送モデルは、世界でも数が限られており、今後、本研究課題で得られる最新の濃度・アイソトポマー観測データおよび放出源アイソトポマー情報を取り入れて総合的に解析することにより、観測される濃度変動の解釈やそれにともなうCH_{4}収支変動を定量化できることが期待される。

## （２）環境政策への貢献

2009年4月、2010年4月、2011年4月に国立環境研究所で行われた「一般公開」ならびに2009年7月、2010年7月、2011年7月に行われた「夏の大公開」においてシベリア観測活動についての紹介を行い、温室効果ガスが人口密集地以外でも着実に増加していることを周知することに貢献した。得られたデータは国内外の研究者と共有し、 $ CH_{4} $ 放出源や $ N_{2}O $ 放出源の推定に利用されている。これらの観測結果やフラックスの推定結果は、IPCC報告書などに引用され濃度増加のメカニズム解明に貢献する予定である。
## 6 ．国際共同研究等の状況

1) 日露科学技術協定、「シベリア生態系の影響を受けた温室効果気体の観測」、B. Belan・ロシア科学アカデミー大気光学研究所・ロシア、連携

2) 日露科学技術協定、「シベリアにおけるランド・エコシステムの温室効果ガス収支」、V. F. Galchenko・ロシア科学アカデミー微生物研究所・ロシア、連携

## 7 ．研究成果の発表状況

## （1） 誌上発表

<論文（査読あり）>

1) H-S. Kim, S. Maksutov, M.V. Glagolev, T. Machida, P.K. Patra, K. Sudo and G. Inoue: Environ. Res. Lett., 6, 035201, doi:10.1088/1748-9326/6/3/035201 (2011)

“Evaluation of methane emissions from west Siberian wetlands based on inverse modeling.”

<査読付論文に準ずる成果発表>

該当せず

くその他誌上発表（査読なし）>

1）町田敏暢、笹川基樹、下山宏、M. ARSHINOV、D. DAVYDOV、A. FOFONOV、O. KRASNOV、N. FEDOSEEV、S. MITIN、須藤洋志、勝又啓一、津田憲次、中澤高清、S. MAKSYUTOV：低温科学，68，9-20（2010）「シベリアにおける温室効果ガスの時空間分布」

## （2） 口頭発表（学会等）

1) T. Machida, M. Sasakawa, K. Shimoyama, N. Tsuda, M. Arshinov, B. Belan, D. Davidov, A. Fofonov, O. Krasnov, N. Fedoseev and S. Maksyutov: XVI International Symposium “Atmospheric and Ocean Optics. Atmospheric Physics”, Tomsk, Russia (2009)

“Long-term Monitoring of Atmospheric Greenhouse Gases in Siberia”

2) T.Saeki, S. Maksyutov and T. Nakazawa: American Geophysical Union 2009 Fall Meeting, San Francisco, California, USA (2009)

“Simulation of atmospheric methane and its carbon and hydrogen isotopes to constrain its global budget”

3）佐伯田鶴、Maksytov Shamil：第16回大気化学討論会(2010)「大気中メタンの炭素・水素同位体比に関する数値実験」

4）町田敏暢：平成22年度千葉大学CEReS寄付部門公開シンポジウム(2010)「大気中CO₂濃度の変動と民間航空機を使ったグローバルな観測展開」

5) O. Uchino, I. Morino, T. Tanaka, Y. Miyamoto, Y. Yoshida, T. Yokota, T. Machida, D. Wunch, P. Wennberg, G. toon et al.: TCON Meeting, Wollongong (2010)

“Progress of validation of GOSAT standard products”

6) H-S. Kim, M.V. Glagolev, S. Maksyutov, T. Machida, P.K. Patra, K. Sudo and G. Inoue:
ENVIROMIS-2010, Tomsk (2010)

“Evaluation of methane emission over Western Siberia”

7) O. Krasnov, S. Maksyutov, K. Shimoyama, H. Suto, A. Nadeev, V. Shelevoi, M. Glagolev, N. Kosykh, T. Machida and G. Inoue: AGU 2010 Fall Meet., San Francisco (2010)

“Automatic chamber observations of methane and carbon dioxide fluxes at West Siberian wetland”

8) S. Maksyutov, M. Glagolev, I. Kleptsova, A. Sabrekov, A. Peregon and T. Machida: AGU 2010 Fall Meet., San Francisco (2010)

“Methane emissions from the West Siberian wetlands”

9) T. Machida, M. Arshinov, B. Belan, N. Fedoseev, S. Mitin and S. Maksyutov : Sixth International Symposium on Non- $ CO_{2} $ Greenhouse Gases (NCGG-6), Amsterdam (2011)

“Temporal and spatial variations of atmospheric methane over Siberia”

10) 佐伯田鶴、Shamil Maksyutov：第17回大気化学討論会（2011）

「大気中メタンの炭素・水素同位体比に関する数値実験」

## （3） 出願特許

特に記載すべき事項はない

（４）シンポジウム、セミナーの開催（主催のもの）

特に記載すべき事項はない

（５）マスコミ等への公表・報道等

特に記載すべき事項はない

## (6) その他

特に記載すべき事項はない

### 8 \. 引用文献

1) M. Rigby et al.: Geophys. Res. Lett., 35, doi:10.1029/2008GL036037 (2008) "Renewed growth of atmospheric methane"

2) E.J. Dlugokencky et al.: Geophys. Res. Lett., 36, doi:10.1029/2009GL039780 (2009)

3) P. K. Patra, M. Takigawa, K. Ishijima, B-C Choi, D. Cunnold, E. J. Dlugokencky, P. Fraser, A. J. Gomez-Pelaez, T-Y Goo, J-S Kim, P. Krummel, R. Langenfelds, F. Meinhardt, H. Mukai, S. O'doherty, R. G. Prinn, P. Simmonds, P. Steele, Y. Tohjima, K. Tsuboi, K. Uhse, R. Weiss, D. Worthy and T. Nakazawa: Journal of the Meteorological Society of Japan, 87(4), 635-663 (2009)

“Growth Rate, Seasonal, Synoptic, Diurnal Variations and Budget of Methane in the Lower Atmosphere”

4) I. Fung, J. John, J. Lerner, E. Matthews, M. Prather, L. P. Steele and P. J. Fraser: J. Geophys. Res.,
96, 13033–13065 (1991)

“Three-dimensional model synthesis of the global methane cycle”

5) K. R. Lassey, D. M. Etheridge, D. C. Lowe, A. M. Smith and D. F. Ferretti: Atmos. Chem. Phys., 7, 2119–2139 (2007)

“Centennial evolution of the atmospheric methane budget: What do the carbon isotopes tell us?”

6) S. Houweling, Dentener, F. and J. Lelieveld: J. Geophys. Res., 105, 17243-17255 (2000)

“Simulation of prein- dustrial methane to constrain the global source strength of natural wetlands”

7) M. Whiticar and H. Schaefer: Philosophical Transactions of the Royal Society, 365, 1793-1828 (2007) "Constraining past global tropospheric methane budgets with carbon and hydrogen isotope ratios in ice"

8) P. Quay, J. Stutsman, D. Wilbur, A. Snover, E. Dlugokencky and T. Brown: Global Biogeochemical Cycles, 13, 445–461 (1999)

“The isotopic composition of atmospheric methane”

$$ ``C/^{12}C $$

10) K. Snover, P. D. Quay and W. M. Hao: Global Biogeochem. Cycles. 14, 11–24 (2000) "The D/H content of methane emitted from biomass burning"

11) G. Saueressig, J. N. Crowley, P. Bergamaschi, C. Bruehl, C. A. Brenninkmeijer and H. Fischer: J. Geophys. Res., 106, 23127–23138 (2001)

“Carbon 13 and D kinetic isotope effects in the reactions of $ CH_{4} $ with $ O(^{1}D) $ and $ OH^{-} $ : new laboratory measurements and their implications for the isotopic composition of stratospheric methane”

12) M. Wahlen, N. Tanaka, R. Henry, B. Deck, J. Zeglen, J. S. Vogel, J. Southon, A. Shemesh, A. Fairbanks and W. Broecker: Science, 245, 286–90 (1989)

13) E. J. Dlugokencky, P. M. Lang and K. A. Masarie: Atmospheric Methane Dry Air Mole Fractions from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network, 1983-2009, Version: 2010-08-12, Path: ftp://ftp.cmdl.noaa.gov/ccg/ch4/flask/event/ (2010)

14) J. W. C. White and B. H. Vaughn: University of Colorado, Institute of Arctic and Alpine Research (INSTAAR), Stable Isotopic Composition of Atmospheric Methane $ \left(^{13}\mathrm{C}\right) $ from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network, 1998-2009, Version: 2010-10-13, Path: ftp://ftp.cmdl.noaa.gov/ccg/ch4c13/flask/event/ (2009)

15) M. V. Glagolev, I. E. Kleptsova, I. V. Filippov, V. S. Kazantsev, T. Machida and S. Maksyutov: Moscow Univ. Soil Sci. Bull., 65, 86–93 (2010)

“Methane emissions from subtaiga mires of Western Siberia: the Standard Model Bc5”
## (3) 理論計算によるアイソトポマ一分別係数の決定

上智大学理工学部

物質生命理工学科 理論分子設計グループ・教授 南部 伸孝

〈研究協力者〉

新潟大学 自然科学系 数理物質科学列理学部 化学科・教授 德江 郁雄

平成21～23年度累計予算額：3，930千円

（うち、平成23年度予算額：1，310千円）

予算額は、間接経費を含む。

[要旨] 本サブテーマは、量子力学に基づく第一原理計算を実施し、大気科学において重要な役割を果たす分子の光分解反応における同位体分別係数を理論的に決定することを目的とする。さらに、得られた基礎データを用い、他のサブテーマと協力し、大気化学シミュレーションを行うことにより、多種の同位体を説明可能な新しい大気モデルの提出を目指すものである。具体的には、 $ N_{2}O $ 、COS、 $ SO_{2} $ 、SOおよび $ H_{2}S $ の光解離反応による同位体分別係数の理論的決定を実施する。研究期間中、厳密に決定した大域的ポテンシャルエネルギー曲面と遷移双極子モーメントを用い、量子力学に基づく第一原理計算を実施した。三原子分子系については量子波束動力学計算を実施し、温度に依存する光吸收断面積の決定を行うとともに、二原子分子であるSOの光解離反応において、より計算精度を上げるためにR行列伝播法に基づくプログラムを開発し、反応動力学計算を実施した。特にこの方法の特徴は、原子核配置空間内の動力学的に接近できる領域を小さな領域に分割し、各々の小さな領域で局所的なR行列を解析的に決め、さらに局所的なR行列を再帰的に組み立てることにより、全相互作用領域にわたるR行列を得ることができることにある。そして、光吸收断面積は最後のR行列から直接求められる。従って、この方法を用いることにより、重なり行列の直交化や安定化の技巧に頼らずとも、散乱行列の対称性やユニタリ性が保たれる。また、複数の電子励起状態への電子遷移を一度に計算し、振動回転励起状態からの電子遷移についても計算を可能にした。

[キーワード] 一酸化硫黄、亜酸化窒素、硫化カルボニル、光分解反応、第一原理計算

### 1. はじめに

大気科学における同位体濃縮現象は様々な分子で観測される。特に有名な反応は亜酸化窒素 $ (\mathrm{N}_{2}0) $ である。この分子は、温室効果ガスの一つであり、成層圈におけるその酸化反応が触媒となりオブン層破壊の主な要因となっている。亜酸化窒素の減少に関する主な反応は式(3-1)に示す光分解反応である。特に、 $ N_{2}0 $ の消失の90%が成層圈の紫外線窓領域におけるこの光分解反応であり、残り10%が $ 0^{(1)}\mathrm{D} $ との反応とされている。

$$ \mathrm{N}_{2}\mathrm{O}(1^{1}\Sigma^{+})+h\nu(\lambda ca.205nm)\rightarrow\mathrm{N}_{2}(1^{1}\Sigma_{g}^{+})+\mathrm{O}(^{1}\mathrm{D}) $$
図(3)-1 Yung と Millerが提唱した零点振動モデル。

歴史的には、この反応による同位体濃縮現象の原因解明が1997年にYungとMillerにより行われた。つまり、零点振動エネルギー（ZPE）モデル $ ^{1)} $ （図(3)-1を参照、分子の電子基底状態の振動準位と電子励起状態を用い、鏡像を取って吸収断面積を近似する）であり、その後、大気科学のほぼすべての研究者がそれに従って解析を続けている。その方法を簡単に示すと以下の通りである。

分別係数 $ \varepsilon $ は以下の式 (3-2) に従って決定される。ここで、 $ N_{2}0 $ の同位体種を i、j、k のインデックスで表し、 $ ^{14}N^{14}N^{16}0 $ 種を分母に用いる。J は光分解速度、 $ I(\lambda) $ は紫外線量、 $ \sigma(\lambda) $ は吸收断面積である。つまり、分別係数は「光分解速度の同位体による比ニ吸収断面積の比」で表され、図 (3)-1 を用いれば、右に示す「2つの線の下に囲まれた面積の比」となる。ところが本を正すと、

$$ \varepsilon_{ijk}=\frac{J_{^{i_{N^{j}N^{k}O}}}}{J_{^{14_{N^{14}N^{16}O}}}}-1=\frac{I(\lambda)\sigma_{^{i_{N^{j}N^{k}O}}}(\lambda)}{I(\lambda)\sigma_{^{14_{N^{14}N^{16}O}}}(\lambda)}-1=\frac{\sigma_{^{i_{N^{j}N^{k}O}}}(\lambda)}{\sigma_{^{14_{N^{14}N^{16}O}}}(\lambda)}-1 $$

簡単な零点振動エネルギーの比となるのであり、実に明快な解明方法である。例として $ {}^{14}N^{15}N^{16}O $ 分子のにおける実験および理論の結果を図(3)-2に示す。青い線が Yung と Miller によるモデルの結果である。

その一方、2003年ごろから地球温暖化問題が広く世界で認識され始め、正確な大気循環モデルを構築する必要があり、それとともに正確な大気観測が行われる。さらに、実験室でも正確な観測が可能となった。ところが、この正確なデータで解析を進めると、図(3)－2から明らかのように、このZPEモデルでは明らかに再現されないことに気づくことになる。実験結果はvon Hessbergら、SelwynとJohnson、Rockmannら、Zhangら、Kaiserらが報告している。ざっと見て、3倍以上差異がある。但し、実験結果にもかなりバラつきがあり、どの結果が真実の結果か少し疑問が残る。これに対し我々は、2004年ごろ、最新の実験結果を理論的にほぼ説明することに成功した。緑色の線が我々の結果である。理由はやはり厳密な量子動力学により光吸収断面積を決定することが不可欠であり、簡単なZPEのみでは不十分となった。 $ ^{2)} $ さらに2008年ごろ、COS分子における同位体分別係数の理論的予測に成功した。 $ ^{3、4)} $
図(3)-2^{14}N^{15}N^{16}O分子の実験及び理論の同位体分別係数(T=283K)。

ただし、図から明らかなように、最新の実験結果と比較し波長に依存する分別係数の特徴が高エネルギー側にシフトしてしまっている。そしてその理由は明確であり、両計算において電子状態の計算における基底関数の精度が不十分なため、励起状態への励起エネルギーが高エネルギー側にシフトしているため、結果が少しシフトした。さらに、より高い電子励起状態への遷移によるしているが、実際には3つの電子状態のみを考慮し、図(3)－2の

結果を得ている。したがって、あと2つの励起状態3A’及び2A”

への遷移を考慮する必要がある。

そしてその延長上に、様々な研究者が現在取り組んでいる硫黄を含む大気循環モデルとそのモデルにおける同位体効果において観測される質量非依存分別（Mass Independent Fractionation、MIF）あるいは質量依存分別（Mass Dependent Fractionation、MDF）

図(3)-3 ヤコビ座標。

についての議論への応用が可能である。特に我々の最近の結果によれば、分子の問題を理論的に解くとき原子核の質量と電子の質量に大きな差異があることから原子核の動きを止めて電子状態の問題を先に解く（ボルン・オッペンハイマー近似（あるいは断熱近似））ことが一般的であるが、この近似の破れである非断熱現象がMIFを起こす原因の1つであると考えている。そこで、その効果を取り入れたプログラム開発が必要となる。

本サブテーマの研究開発方法の概要は、量子力学に基づく第一原理計算を実施し、大気科学において重要な役割を果たす分子の光分解反応における同位体分別係数を理論的に決定することである。さらに、得られた基礎データを用い、他のサブテーマと協力し、大気化学シミュレーション
ンを行うことにより、多種の同位体分子種を説明可能な新しい大気モデルの提出を目指す。具体的には、 $ N_{2}O $ 、COS、 $ SO_{2} $ 、SOおよび $ H_{2}S $ の光解離反応による同位体分別係数の理論予測を実施する。

## 2 ．研究開発目的

理論による分子の化学反応の研究を行う上で最初に導入される近似が、ボルン・オッペンハイマー近似（断熱近似）である。つまり、電子が原子に比べ極端に軽いことから分子の運動を、電子の運動の問題と核の運動の問題とに分けて扱うことが可能となる。さらに、その電子の運動について時間非依存のシュレーディンガー方程式を固有値問題等に帰着させて解くことになる。そして一般には、その計算のことを「電子状態計算あるいは量子化学計算」と呼んでいる。一方、光解離過程等の現象を扱う上で重要な要請として、電子励起状態を正確に記述しなければならないという点がある。本研究では、多配置参照型の配置間相互作用計算(MR-CI)を用いる。電子の運動の問題を解決した後、化学反応の収率やメカニズムを把握するため、さらに量子論に従った核の運動を解く必要がある。このことを「化学反応動力学あるいは量子動力学」という。一般には、時間依存シュレーディンガー方程式を解く方法と時間非依存のシュレーディンガー方程式を解く方法の2つが存在する。時間非依存の方法としては、緊密結合方程式を数値的に解くR行列伝播法が有名である。本研究では計算量を考え、二原子分子のみこの方法に基づくプログラムを開発し実行する。その一方、三原子分子系については実空間波束発展法（あるいはChebyshev次数発展法）用いる。この様な計算手順により光吸収断面積が決定され、同位体分別係数を決定することを目的とする。

### 3. 研究開発方法

## （1） 電子状態計算 $ -N_{2}0 $ 分子における計算例

電子基底状態と電子励起状態のポテンシャルエネルギー面を決定するために量子化学計算を行う。特に、励起状態においても約0.1eV以下の計算精度を要求されるため、最低でもTriple zeta(TZ)基底を必要とする。具体的には、基底関数にDunningらのaug-cc-pVTZ $ ^{5} $ を用いた。理由は励起状態の記述に優れているからである。次に制限つきHartree-Fock計算を行い、分子軌道(MO)を調べる。核配置は、分子の安定構造付近を用いる。C $ _{s} $ 対称性のもと基底状態のMOは以下のように記述されることがわかった。

$$ \begin{aligned}&(1a')^{2}(2a')^{2}(3a')^{2}(4a')^{2}(5a')^{2}(6a')^{2}(7a')^{2}(1a')^{2}(8a')^{2}(2a')^{2}\\ &(9a')^{2}(10a')^{0}(3a')^{0}(11a')^{0}(12a')^{0}\cdots\cdots(1^{1}A'),\\ \end{aligned} $$

特にその中で重要な分子軌道は以下のように特徴付けられる。

分子軌道：その特徴

7a' : NNO σ 結合、 σ_{NN} + · σ_{NO}

8a' $ \quad $ と $ \quad $ 1a" : NNO $ \pi $ ・結合、 $ \pi_{NN} + \pi_{N0} $

9a' $ \quad $ と $ \quad $ 2a" : NN $ \pi $ と NO $ \pi $ *結合、 $ \pi_{NN} $ + $ \pi $ * $ _{NO} $

10a' $ \quad $ と $ \quad $ 3a" : NN $ \pi^{*} $ と NO $ \pi^{*} $ 結合、 $ \pi^{*}_{NN} + \pi^{*}_{NO} $

11a' : NN σ と NO σ*結合、σ_{NN} + σ_{NO}

12a' : NN $ \sigma^{*} $ と NO $ \sigma^{*} $ 結合、 $ \sigma^{*}_{NN} + \sigma^{*}_{NO} $
これらの軌道のうち活性軌道として7軌道を選び、制限付きHartree-Fock計算により求めたM0を初期波動関数にして多配置参照SCF計算を行い、分子軌道を再度決め直す。特にその7軌道は以下の通りである。

$$ [(6a^{\prime})(7a^{\prime})(1a^{\prime})(8a^{\prime})(2a^{\prime})(9a^{\prime})(11a^{\prime})]^{12} $$

図(3)-4 直線形におけるポテンシャル曲線。

図(3)-5 非直線形におけるポテンシャル曲線。

つまり、7 軌道に 12 電子詰めることによって出来る全ての電子配置を基に多配置参照 SCF 計算を行い、自然軌道の決定を行う。自然軌道を使う理由は、一般には制限付き Hartree-Fock 計算で決定される分子軌道は電子基底状態を記述するのに適したものとなるが、電子励起状態をその分子軌道を使って計算すると、かなり不安定な状態を記述してしまうことが多々あるからである。そこで、後述べる電子状態平均を取った多配置参照 SCF 計算が不可欠となる。得られた分子軌道を使い、多配置参照型の配置間相互作用計算

図(3)-6 遷移双極子モーメント（非直線形）。

(MR-CI) $ ^{6,7)} $ を行い、電子励起状態のポテンシャルエネルギー及び遷移双極子モーメントを求める。特に、得られたMR-CIのエネルギーをもとにDavidson補正を行い、二電子励起以上の効果を取り込む $ ^{8-10)} $ 。一方、解離極限で $ 0(^{1}\mathrm{D})+\mathrm{N}_{2}(\tilde{\mathrm{X}}^{1}\Sigma_{g}^{+}) $ 状態へ相関する光解離過程であることから、5つの電子状態を均等に記述できるような分子軌道が必要になる。つまり、状態平均多配置参照SCFを行う必要がある。この場合の状態数は3つのA'と2つのA''状態である $ (0(^{1}\mathrm{D}) $ の電子状態を
C_{s} 対称下で考えるとこのような組み合わせがいいことに気付く）。このような方法で量子化学計算を行うが、最終的に解かなければいけない MR-CI 行列は、大規模疎行列の固有値問題となる。この系の場合、典型的な大きさは約 1、000 万次元となった。一般な数値解法には、大規模対称行列を三重対角化する手法としてアーノルド法やランチョス法がある。特に大規模疎行列を部分空間に射影し、新しい比較的小さな行列を作る。そして，できた小さな行列の固有値を求め、それが大規模疎行列の固有値の良い近似となる。ここで、この射影操作のことをクリロフ部分空間法と呼ばれ、大規模疎行列の対角化問題によく使われてきている。一方、ハウスホルダー法等が一般的には知られているが、アルゴリズムの関係上コンピュータで利用できるメモリーの上限をすぐ超えてしまい、到底解くことができない行列のサイズである。従って、具体的にはアーノルド法に基づくクリロフ部分空間法を用い、エネルギーの低い方から必要な電子状態数（固有値数）のみを近似的に求める。本研究では、九州大学情報基盤研究開発センター所有の日立 SR16000 および富士通 Primquest および Primergy からなる 256CPU コアを用い、約 1 年間計算を行うことにより、約 5、000 の核座標配置についてのポテンシャルエネルギーデータを得ることができた。さらに、他の分子についても同様に実施し、約 2 年間計算を続けた。

## （2） SO 光解離過程

一酸化硫黄(SO)の光解離反応を、反応式（3-3）に示す。

$$ \mathrm{SO}(\mathrm{X}^{3}\Sigma^{-})\stackrel{h\nu}{\rightarrow}\left\{\begin{aligned}&\mathrm{S}(^{1}\mathrm{D})+\mathrm{O}(^{3}\mathrm{P})\\ &\mathrm{S}(^{3}\mathrm{P})+\mathrm{O}(^{3}\mathrm{P})\end{aligned}\right.\quad(3-3) $$

図(3)-7 反応に関与するポテンシャル曲線。

この反応は三重項電子状態が24個関与するが、本研究では電子基底状態 $ \left(\mathrm{X}^{3}\Sigma^{-}\right) $ から許容遷移である $ A^{3}\Pi,B^{3}\Sigma^{-},C^{3}\Pi,C^{3}\Pi,3^{3}\Sigma^{-} $ 電子状態（図(3)-7を参照）を考慮する。さらにその中でも、 $ C^{3}\Pi,C^{3}\Pi $ 、 $ B^{3}\Sigma^{-},3^{3}\Sigma^{-} $ の電子励起状態間において、ポテンシャルエネルギー曲線上の6 eV～7 eV付近に、擬交差が2つ存在する（図(3)-8を参照、赤丸内に擬交差の位置を示す）。この擬交
差において、先に述べた分子の運動を理論的に解くとき原子核の質量と電子の質量に大きな差異があることから、原子核の動きを止めて電子状態の問題を先に解く、ボルン・オッペンハイマー（BO）近似（断熱近似）の破れによる非断熱遷移現象が起きることが予想される。さらに、このような一自由度系においてこの非断熱現象により生じる量子効果が、多自由度系より明確に表れることが一般に知られている。そこで、このような特徴を持つ光解離反応を理論的に解明するにあたり、非断熱遷移を考慮した厳密化学反応動力学が必要となる。

具体的には、以下に示す非断熱結合要素（NACME）を求める必要がある。

$$ \left\langle\phi_{1}\right|\frac{\partial}{\partial Q}\left|\phi_{2}\right\rangle\qquad(3-4) $$

ここで2つφが異なる電子状態の電子波動関数を示し（図(3)-8aの場合、C状態とC'状態の電子波動関数を示す）、Qは原子核の座標を示す。つまり、電子の運動が速ければ、原子核の座標で電子の波動関数を微分しても要素(3-4)は、ほぼゼロと近似できる。これが断熱近似に相当する。ところが、この要素が非ゼロになる時が現れる。その時、粗いBO近似に従い、非断熱遷移現象を考慮することになる。そのような状況になるのが、2つの状態が近接した図(3)-8のような状況である。本研究では、理論に基づく仮想実験が可能であることから、非断熱遷移を考慮した場合と考慮しない場合の結果を比較すること可能となり、非断熱遷移現象が質量非依存分離（MIF）あるいは質量依存分離（MDF）の要因となるかを判断する大変興味が持たれる系である。

図(3)-8a 非断熱遷移に関与するポテンシャル曲線。

図(3)-8b 非断熱遷移に関与するポテンシャル曲線。

## （3） $ SO_{2} $ 光解離過程

二酸化硫黄（ $ SO_{2} $ ）の紫外光による電子遷移過程（ $ C^{1}A^{\prime}\leftrightarrow X^{1}A^{\prime} $ ）に伴う光解離反応式を以下に示す。

$$ SO_{2}(1^{1}A_{1})\xrightarrow{h\nu}\begin{cases}SO(1^{3}\Sigma^{-})+O(^{3}P)\\ O_{2}(1^{3}\Sigma^{-})+S(^{3}P)\end{cases} $$
電子励起状態 C^{1}A’ 及び D^{1}A’ 間に擬交差が存在する。そのため非断熱遷移過程が不可欠である。つまり、SO の時と同様に要素 (3-4) を求め、非断熱遷移を考慮する。

一方、近年盛んに議論がなされている現代大気循環や古代大気循環における硫黄を含む分子の連鎖反応の解明において、重要な反応の1つである。

## （4） COS 光解離過程

硫化カルボニル（COS）の紫外光による光解離反応を以下に示す。

$$ OCS(1^{1}\Sigma^{+})\xrightarrow{hv}CO(1^{1}\Sigma^{+})+S(^{1}D)\qquad(3-6) $$

この系は、 $ N_{2}O $ 分子と等原子価分子であることから電子状態の計算方法は $ N_{2}O $ 分子と似たものとなる。つまり、5つの電子状態が光解離過程に関与している。また、この分子も近年、古代大気循環にて脚光を浴びている分子の1つである。

## （5） $ H_{2}S $ 光解離過程

硫化水素（ $ H_{2}S $ ）の紫外光による光解離反応を以下に示す。

$$ \mathrm{H}_{2}\mathrm{S}(1^{1}\mathrm{A}_{1})\xrightarrow{h\nu}\left\{\begin{aligned}&\mathrm{H}_{2}(1^{1}\Sigma_{\mathrm{g}}^{+})+\mathrm{S}(^{1}\mathrm{D})\\ &\mathrm{SH}(1^{2}\Pi)+\mathrm{H}(^{2}\mathrm{S})\end{aligned}\right. $$

この系は、7つの電子状態が光解離過程に関与している。また、光吸収直後のフランク-コンドン領域に擬交差があるため、光吸収スペクトルにおいて非断熱遷移を考慮した計算が不可欠である。

上記のように分子ごとの電子状態が全く異なるため、それぞれの分子に対し一つ一つ正確な計算入力を作成し実行する必要がある。さらに我々はこれまで3原子分子以下の小さな系のみを扱ってきたが、 $ SO_{3} $ などの四原子分子系においても研究を進める必要がある。それについては後で述べる。

一方、得られたポテンシャルエネルギーと遷移双極子モーメントのデータ点から大域的ポテンシャルエネルギー曲面／曲線および大域的遷移双極子モーメント曲面／曲線を決定する。本研究では、3自由度系（3原子分子系）では直接数値内挿を行うIMLS/Shepard内挿法 $ ^{11} $ を用い、1自由度系（2原子分子系）ではスプライン内挿法を用いた。すべての量子化学計算を、非経験的分子軌道計算プログラムMolpro 2006.1を用いて行った $ ^{12} $ 。

## （6） 化学反应勐力学

第一原理計算は、大きく分けると二つの計算から構成される。先に示した電子状態計算と化学反応動力学である。ここでは、開発した化学反応動力学計算のためのプログラムを簡単に紹介する。

## （7） 1次元時間非依存シュレーディンガー方程式の解法 -S0 光解離過程へ応用

先に述べた通り、SO の紫外光による 5 つの電子遷移過程に伴う光解離過程（以下の反応過程）に
よる原子核の運動を解くことにより、光吸收断面積の決定を行う。特に4つの電子励起状態間に2つの擬交差存在する。そのため非断熱遷移過程が不可欠である。

1 次元系であることから以下に示す時間非依存シュレーディンガー方程式（緊密結合方程式）を R 行列伝播法を用い解くことにより、核の波動関数を求める。具体的には、式(3-8)を解き、

$$ \left[-\sum_{\alpha}^{M}\frac{\hbar^{2}}{2M_{\alpha}}\Delta_{a}+\mathrm{E_{m}}\right]\Psi_{m}(\mathrm{R})+\sum_{n}\Lambda_{m n}\left(\mathrm{R}\right)\Psi_{n}(\mathrm{R})=\mathrm{E}\Psi_{m}(\mathrm{R}) $$

得られた波動関数より以下の式を用いアインシュタインの B 係数を決定する。M(R)は遷移双極子モーメント関数である。

$$ \sigma(\mathbf{v})=\int\Psi_{n u c l}^{\prime}\left(\mathrm{R}\right)\mathrm{M}\left(\mathrm{R}\right)\Psi_{n u c l}^{\prime \prime}\left(\mathrm{R}\right)\mathrm{d}\mathrm{R} $$

非断熱遷移を考慮するに当たり透熱表現を用いた方が容易である。そこで、混合2状態問題を解くための透熱ポテンシャル $ V_{11} $ 、 $ V_{22} $ と非断熱結合要素 $ V_{12} $ を以下に示す式を用い決定した。

$$ \begin{aligned}\mathbf{V}_{11}=&\mathbf{E}_{11}\cos^{2}\beta+\mathbf{E}_{22}\sin^{2}\beta\\\mathbf{V}_{22}=&\mathbf{E}_{11}\sin^{2}\beta+\mathbf{E}_{22}\cos^{2}\beta\\\mathbf{V}_{12}=&\left(\mathbf{E}_{22}-\mathbf{E}_{11}\right)\cos\beta\sin\beta\end{aligned} $$

但し、状態混合角 $ \beta $ は非断熱結合要素 (3-4) を下記の式 (3-11) に沿って積分することにより決定される。

$$ -\frac{\partial\beta}{\partial Q}+\left\langle\phi_{1}\right|\frac{\partial}{\partial Q}\left|\phi_{2}\right\rangle\approx0 $$

得られた透熱ポテンシャルを用い、R行列伝播法を用い時間に依存しないシュレーディンガー方程式を解く。

## （8） R行列伝播法（透過型散乱問題） $ ^{13)} $

## 1 ）基本方程式

時間非依存シュレーディンガー方程式から出発し、図(3)-9のような1次元2準位からなる系における非断熱トンネル問題を解いてみる。但し、x=0を境に正をα-セグメントとし、負をβ-セグメントとする（図(3)-9には、透熱結合項 $ V_{12} $ および $ V_{21} $ は記載していない）。

図3-9 1次元2準位透熱ポテンシャル。