順圧Sモデルを用いて1950年から50年間の時間積分を多数行い, 仮想的外力 応答によるトレンドと, それに重なる内部変動の特徴を行本(2005)の手法で解析 した.
SVD解析の左ベクトルUと右ベクトルVを外力とした実験Aと実験Bでは, 外力に加えた少しのトレンドに応答してAOIの増加傾向が見られた. そこに大き な振幅の十年スケールの変動が確認され,観測による同様の時系列と比較すると非 常に良く対応していた.全てのメンバーとアンサンブル平均の時系列においても緩 やかな上昇トレンドが再現されている. 一方で実験Cでは外力にトレンドを加え なかったためか各メンバーとアンサンブル平均ともに上昇トレンドは再現できな かった. 内部変動と外部強制応答によるEOF1の構造は,左ベクトルUを外力とし た実験Aの外部強制応答でPNAのような構造が現れた以外は全て観測に見られ るようなAOパターンとなった.
内部変動パターンが全てAOとなることはIPCCのモデル群における実験と同 様の結果である. IPCCのモデル群のときに外部強制応答においてAOパターンが 再現されたものは, 観測や内部変動に見られるようなきれいなAOの構造は現れな かったのに対し,今回順圧Sモデルによって再現された外部強制応答のAOパター ンは実験A, 実験Cともに観測をよく再現しているのがわかる.
これらの結果から,極めて単純な順圧大気大循環モデルを用いてもIPCCのAR4 モデルと同様な特徴が再現され,北極振動の数十年スケールの長期変動が,大気の 内部変動として十分に説明できることが示唆された.
5 結論
本研究ではIPCCのAR4で使用された10個の大気海洋結合モデルを用いて20 世紀, 21世紀の地球温暖化時に見られるAOの解析を行った.
SLP場においてどのモデルも例外なく20世紀, 21世紀をつなげて計算した冬平 均のEOF1には現実的なAOの構造が現れ, AOIの年々変動における時系列には上 昇トレンドを示すものがいくつか見られた. 一方でモデルは20世紀に観測された 全球平均地上気温の変化をよく再現し, 21世紀の100年間では約2℃〜4℃の上昇 傾向を示した.
以上よりモデルのAOや気候変化における信頼性を十分評価した上で, 十年ス ケールの変動についてアンサンブルメンバーが20世紀, 21世紀ともに十分揃って いる4個のモデルのデータセットを使用し, AOを外部強制に対する応答と内部変 動に分離して解析した. その結果, 内部変動に見られるEOF1の構造はどのモデル も20世紀, 21世紀ともに共通してAOパターンとなる一方で, 外部強制応答のパ ターンはモデル間,シナリオ間で大きく異なることが示された.外部強制応答のAO 的変動は20世紀末以降, 大きな上昇トレンドとなり,今後観測に見られるAOのト レンドのかなりの部分が温室効果ガスの増加という人為的な強制力によるものと 考えられる.
また, 10個のモデルにおけるアンサンブル平均とそれを差し引いたバイアスに 見られるEOF1を計算したところ, 両者は20世紀, 21世紀ともに観測に見られる ようなAOパターンを再現し, モデルのアンサンブル平均が示すAOIの時系列は 20世紀末以降上昇傾向となった.
最後に極めて簡単な順圧大気大循環モデルを用いて同様の解析を行ったところ, IPCCのAR4モデルと同じような特徴が再現され, AOの数十年スケールの長期変 動が, 大気の純粋な内部変動として十分に説明できることが示唆された.
21世紀のSRES-A1Bシナリオのもとでの外部強制に対するAOIの増加など, 地
球温暖化に伴う中高緯度の気候変化とAOの変化にはよい対応が見られたが, その 空間パターンはモデル間で大きく異なった。また20世紀にモデルが再現したAOI が完全に観測を表現できてない部分もあり、AOの将来予測には多くの解明すべき 課題が残されている。
謝辞
本研究を進めるにあたり, 指導教員である筑波大学計算科学研究センターの田中 博教授には本研究の動機となる論文の紹介, 解析手法の提案, また筆者の質問等に おいて終始丁寧な御指導, 御鞭撻を賜り, 心より感謝しております.
また同大学生命環境科学研究科の寺崎康児さん,大庭雅道さん,加藤真悟さん,近 藤圭一さん, 山崎真吾さん, 釜江陽一さん, 同大学環境科学研究科の鈴木一歩さん, 瀬田繭美さんにはセミナーおよび様々な場において多数の御助言,御意見を頂きま した. そして, 気象研究所の松枝未遠さんには研究を進めるにあたり数多くの貴重 なアドバイス, 御指導をして頂きました.
さらに同大学大学院生命環境科学研究科の木村富士男教授,林陽生教授,上野健 一准教授,植田宏昭講師および同大学計算科学研究センターの日下博幸講師や大学 院生の皆様には, 方針発表, 中間発表,最終発表,ポスター発表の場で御助言を頂き 誠に有難うございました.
最後に,共に卒論作業を進めた地球科学専攻の4年生の皆様には時折よき相談相 手となって頂きました.
本論文は以上の皆様の御協力により完成させることができました. 心より感謝い たします.
尚, 本研究で用いた主な図は, GMT (The Generic Mapping Tools; Wessel and Smith 1991)にて作成した.
Appendix EOF 解析
EOF解析とは
AOは北半球におけるSLPをEOF解析(主成分解析)したものの第一主成分に より定義される.
EOF解析は, いくつかの地点における時系列データの主要な変動パターンを抽 出して変動の特徴を把握するための統計的手法であり,北半球における大気の変動 パターンを見るためにこれを行う.
一般に何らかの相関関係があるp個の地点でのデータx1, x2,· · ·, xp(p ≥ 2)で, 時系列方向にN個のデータx1λ, x2λ,· · ·, xpλ (λ = 1,2,· · ·, N)が得られたとする.
これらのN個の時系列データは,それぞれp個の地点での気象要素が相互に関連の ある変動を示しているとみなせるので, この変動を説明する関数として, p個の変 量の一次結合で表すことができ, これを
z =l1x1+l2x2+· · ·+lpxp (71) と仮定し,l1, l2,· · ·, lpを変化させ,
∑p i=1
l2i = 1 (72)
の条件の下で, zの分散が最大になるときのzを第一主成分という. このときの係 数をl1i(i= 1,2,· · ·, p)で表すと,
z1 =l11x1+l12x2+· · ·+l1pxp (73) となる.
次にz1とは無相関なzのうちで, 式(73)を満たす最大の分散を持つz2が決定で きる. このz2を第二主成分という. この時の係数をl2i(i = 1,2,· · ·, p)とするとz2
は,
z2 =l21x1+l22x2+· · ·+l2pxp (74) と表すことができる.以下同様にして第m成分zmまでで全変動の大部分が説明で きればこれ以上を求める必要はない.ここで,zmを,
zm =lm1x1+lm2x2+· · ·+lmpxp (75)
と書く. 従って各係数は,
l2α1+l2α2+· · ·+lαp2 =
∑p k=1
l2αk = 1 (76)
を満たし 各主成分は 以下の条件を満たすように定まる.
第一主成分z1の分散が最大である
第i主成分と第j主成分(i̸=j)は直交する(無相関である)
各主成分の計算は,p×pの分散共分散行列の固有値問題に帰着する. これについて は次節に記すこととする.
全球再解析データに対してEOF解析をする際には注意が必要である. なぜなら, 全球再解析データはどの緯度帯でも,経度方向には同じ数のグリッドが存在するか らである. 地球は球形であるので, 全球再解析データのままでは, 高緯度になるほ ど各グリッドあたりの面積が小さくなる.そのため,面積荷重を考慮せずにEOF解 析を行ってしまうと,低緯度と比較して高緯度の変動が相対的に大きく見積もられ てしまい,実際に変動を特徴付けるパターンとは異なるパターンが現れる恐れがあ る.このグリッドの偏りを修正するために, 緯度方向に面積荷重√
cosθをデータに かけて偏りを解消する.
EOF解析における固有ベクトルの計算方法
一般に,p変量データ,すなわちN個の各標本について,p種類の変量x1, x2,· · ·, xp が測定されたとし, そのデータに基づき, 主成分z1, z2,· · ·, zm(m ≤ p)を求める方 法について以下で述べる.
まず, 得られたデータより各変量の分散共分散σij を求め, それを
S=
σ11 σ12 · · · σ1p
σ21 σ22 · · · σ2p ... ... . .. ... σp1 σp2 · · · σpp
(77)
と表す. ここでのσij は各変量の分散共分散の不偏推定値であり, σij = 1
N −1
{ n
∑
λ=1
(xiλ−x¯i)(xjλ−x¯j)
}
(78)
¯ xi = 1
N
∑N λ=1
xiλ (79)
とし, さらに相関行列を求め, それを
R=
r11 r12 · · · r1p r21 r22 · · · r2p ... ... . .. ...
rp1 rp2 · · · rpp
(80)
(rii = 1;i= 1,2,· · ·, p) と表す.
ここで, p個の変量x1, x2,· · ·, xpからなるベクトル変量を
X=
x1 x2 ... xp
(81)
とする. ここで,Xのとりうる量は,
x11 x21 ... xp1
,
x12 x22 ... xp2
,· · ·,
x1N x2N ... xpN
また, 各主成分を式(71)の形で求めるために各主成分の変量xiの係数をそれぞれ ベクトル
l1 =
l11 l21 ... lp1
,l2 =
l12 l22 ... lp2
,· · ·,lm =
l1m l2m ... lpm
で表し, lαの転置行列をlTα で表すと主成分は
z1 z2 ... zm
=
lT1 lT2 ... lTm
X (82)
と書くことができる.ただし, 各lTα は単位ベクトル(lTα ·lα = 1)とする.
これより,各成分zαが条件(76)を満たすように係数を求めればよい.したがって, V{z1}=v{lT1X}=lT1V{X}l1 =lT1Sl1 (83) z1の分散V{z1}が最大になるようなl1を求めるために, ラグランジュの未定乗数 法を用いる(次節参照). ここで, ラグランジュ未定乗数をλとして,
v =lT1Sl1−λ(lT1l1−1) (84) の両辺をlT1 で微分すると,微分方程式
∂v
∂lT1 = (S−λI)l1 (85)
I:p次の単位行列
が得られ, ラグランジュの未定乗数法よりこの微分方程式が0になるときにz1の 分散V{z1}は最大になる. すなわち,以下の連立方程式
(S−λI)l1 = 0 (86)
を解けばよい. ここで,式(82)のベクトルl1の要素がすべて0以外の解を持つため にはλが固有方程式
|S−λI|= 0 (87)
の解(固有値)でなければならない.また,z1の分散V{z1}を最大にするものが存在 するとすれば, 式(85)より
V{z1}=λ (88)
であり, 式(87)の固有値に等しくなるので, z1における係数ベクトルl1としては, 式(87)の最大の固有値λ1に対応する単位固有ベクトルとして第一主成分z1を決定 することができる.一方,式(87)の固有値を展開するとλについてのp次方程式に なるので, p個の固有値(重複するものを含めて数えると)が存在し,しかもSは非 負の対称行列であるから固有値はすべて非負の実数である.その中の最大のものを λ1としてl1を求める.さらに大きさの順にλ2, λ3,· · ·, λmを抽出し,λk(1≤k ≤m) に対応する単位固有ベクトルを第k主成分zkにおける係数ベクトルlkとして, 全 部でm個の主成分を求めるとこれらは条件を満たす.