火星における対流励起重力波と熱圏への影響
渡邊歩佳

東京大学理学系研究科
地球惑星科学専攻
宇宙惑星科学講座
平成 26 年 2 月提出

要旨
惑星の熱圏の熱構造は大気散逸や大気化学に影響しており、これがどのようなエネルギー収支で決まるのかは重要課題である。一般的には、太陽紫外線による加熱と熱伝導による冷却のバランスが基本であり、そこに大気大循環による熱輸送が影響すると考えられている。地球ではさらに、対流圏で生じた重力波が鉛直伝搬した
のち散逸することで熱圏下部を加熱していることが指摘されており、そのような内部波による大気加熱も注目されている。他の惑星において波動による大気加熱の研究は少ないが、火星の熱圏下部においては重力波と考えられる密度変動がとらえられており、そのような波動が熱圏の加熱に寄与する可能性が指摘されている。特に短周期かつ鉛直波長の長い重力波は分子拡散による減衰を受けにくいため、熱圏まで伝搬することができ、エネルギー収支に影響することが考えられる。ここでは火星大気を対象とし、波動の励起プロセスと鉛直伝搬特性について理論的考察を行う。
まず、非静力学の 2 次元対流モデルを構築し、高高度まで伝搬するような重力波の生成を考察した。対流モデルにおける加熱源は、地表面が太陽光により暖められることを想定した一様加熱と、局所的に分布したダストが太陽光の熱を吸収し大気を暖めることを想定した局所加熱を考えた。こうして再現された重力波の熱圏下部における特性は、火星熱圏大気中の密度変動の観測結果と、水平波長・密度振幅がよく一致していた。このことから、火星の熱圏で観測された密度変動には対流によって励起された重力波による寄与が含まれる可能性が示された。
次に、対流による重力波励起実験で観察された波のパラメーターをもとに、重力波の鉛直伝搬特性について線形解を用いて調べた。線形モデルは、水平方向と時間方向に振動解を仮定したうえで鉛直構造を行列演算によって解くものである。熱圏に到達した重力波は分子粘性により減衰し、大気加熱をもたらすが、加熱率の高度分布には 2 つの極大が見られた。一つは振幅が最大となる高度付近で生じる顕熱フ
ラックスの発散によるものであり、もう一つはそれより高い高度で生じる運動エネ
ルギーの粘性散逸によるものである。これらの加熱率は紫外線によるものと比べても無視できない大きさとなった。これらの結果から、波動による熱圏へのエネルギー
輸送が熱圏のエネルギー収支を考える上で無視できない重要なものであることが示された。

Abstract
Considering that the thermal structure of a planetary thermosphere inﬂuences
atmospheric escape and atmospheric chemistry, what kind of energy balance determines the thermal structure is an important issue. In general, it is considered that abalance between heating by solar ultraviolet radiation and cooling by molecular heat
conduction is fundamental, and the energy transport by atmospheric general circulation also plays an important role. As for the Earth s atmosphere, it is also pointed
out that gravity waves which are generated in the troposphere propagate verticallyupward and dissipate in the lower thermosphere to heat this region. Though only afew studies have been done on the heating of atmospheres by internal waves in otherplanetary atmospheres, wavelike density ﬂuctuations indicative of gravity wave wereobserved in the Martian lower thermosphere, and the possibility that these waves
contribute to the heating of the thermosphere has been suggested. Especially gravity
waves with short periods and long vertical wavelengths will be able to propagate to
the thermosphere because they suﬀer less damping by various dissipation processes;
such waves can play key roles in the thermosphere energy balance. Here we focus on
the Martian atmosphere, and examine theoretically the excitation and propagationof gravity waves.
In this study, I constructed a two-dimensional convection model based on a nonhydrostatic meteorological model to examine the excitation of waves which can propagate to high altitudes. The heating source of the convection model is considered ahorizontally-uniform heating of the atmosphere near the surface by sunlight and alocalized heating caused by absorption of sunlight by a dust cloud. The horizontal
wavelengths and the amplitudes of the gravity waves reproduced in the lower thermosphere well coincide with the results of observations of density ﬂuctuations in thelower thermosphere. Based on this result, it is probable that the waves observed in
the Martian thermosphere are excited by convection.
Next, I investigated the properties of vertical propagation and dissipation using
linear wave solutions bases on the wave parameters observed in the numerical experiments of the convective generation of gravity waves. The linear model solves
the vertical structure by matrix operation under the assumption of a wave solution
with respect to the horizontal direction and the time. Gravity waves reaching thethermosphere are damped by molecular diﬀusion and heat the atmosphere. Thevertical proﬁle of the heating rate shows two maxima: the lower one is located near
the altitude where the amplitude peaks, and is generated by the sensible heat ﬂux
divergence, and the upper one is generated by the viscous dissipation of kinetic energy. These heating rates are comparable to other dominant processes such as EUV
heating. The results suggest that the energy transport by gravity waves can play
an important role in the energy balance of the Martian thermosphere.

3

目次
1

2

序論

1.1

火星大気の構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

熱圏のエネルギー収支 . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

火星大気における重力波 . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4

本研究の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 次元対流実験

11

2.1

基本方程式

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2

計算領域と加熱強制

2.3

2.4

2.5
3

4

1

. . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2.1

一様加熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2.2

局所加熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

その他の計算設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3.1

初期値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3.2

境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.3.3

時間刻みと数値粘性 . . . . . . . . . . . . . . . . . . . . . . . 18

結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.1

一様加熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.2

局所加熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

重力波の線形波動解

32

3.1

基本方程式

3.2

解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.3

計算領域と波動強制

3.4

その他の計算設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

. . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.4.1

背景場 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.4.2

分子粘性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4.3

境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5

結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.6

考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

まとめ

43

i

序論

1

惑星大気中を鉛直伝搬する波動は、運動量やエネルギーを上層まで輸送し、波動の種類によっては熱圏大気にまで影響を及ぼすことが予想される。地球では波動による運動量やエネルギーの輸送が大気大循環を考える上で重要となり、観測と数値モデリングによる研究が多くなされている。しかし、他の惑星を対象としたそのような研究はほとんどなされていない。
本研究では火星大気を対象とし、対流から生じた重力波が熱圏へエネルギーを輸送するプロセスを、対流による励起過程も含めて数値計算により明らかにする。
本節では、1.1 節で火星大気の基本的な構造を概観し、1.2 節で熱圏のエネルギー
収支について説明し、惑星の気候を考える上で熱圏温度が重要となることを示す。

1.3 節では火星大気中を伝搬する波動の観測と波動が熱圏に及ぼす影響を数値計算した先行研究を紹介する。1.4 節では本研究の目的について述べる。

1.1

火星大気の構造

火星は太陽からの距離が地球-太陽間距離の約 1.5 倍、大きさ (半径 3396 km) は地
球 (半径 6378 km) の約半分で、質量は約 10 分の 1 である。自転軸が 25 度傾いているため地球と同様に四季が存在する。
火星の大気は地球に比べて希薄であり、大気圧は 6-9 hPa である。それに対して
スケールハイトは約 11 km あり、約 6 km である地球よりも高い。これは火星の重力が地球の約 4 %しかないためである。大気組成は約 95 %が二酸化炭素、2.7 %が窒素、残りはアルゴン (1.6 %) などであり、これもまた地球とは大きく異なっている。さらに両極には二酸化炭素の氷が存在し、それが季節によって昇華と凝結を繰り返すため、大気量は大きく変化する。二酸化炭素が大気のほとんどを占めているが、希薄であるために温室効果の影響は小さく、平均的には地球よりも低温である(図 1.1)。気温は平均的には地表から高度約 100 km までは高度とともに低下し、およそ 120 km 以上の熱圏では太陽紫外線による加熱のために再び温度が上昇する (図

1.2)。

1

図 1.1: 火星着陸機が観測した火星の大気密度の高度分布 (左) と温度の高度分布(右)(Schoﬁeld et al., 1997)

図 1.2: 火星熱圏の温度分布 (Bougher et al., 1998)

1.2

熱圏のエネルギー収支

熱圏の温度構造は、大気力学のみならず大気種ごとのスケールハイト、化学反応
係数、電離気体の再結合係数などを通じて大気組成の分布にも関わっている。熱圏の大気組成は大気散逸への影響という観点から特に興味深い。
火星では現在も大気が散逸していることがわかっているが、その詳細なメカニズ

2

ムは未だ不明である。ジーンズエスケープや流体力学的散逸のような熱的散逸のみ
ならず、スパッタリングのような非熱的散逸さえも熱圏の温度構造に大きく左右される可能性がある。Leblanc and Johnson (2002) は、ピックアップイオンによる火星大気のスパッタリングが紫外線強度によってどう変わるかを調べた。ピックアッ
プイオンとは大気上層の中性の原子・分子が太陽紫外線により電離されることで生
じる O＋ イオンなどが太陽風電場で加速されたもので、これらが大気中の中性粒子に衝突することで新たに原子や分子を宇宙へと追い出すのがスパッタリングである。

Leblanc and Johnson (2002) はピックアップイオン一個あたりがはじき出す粒子の数を、Monte Carlo model を用いて紫外線強度別に算出している。紫外線強度は、現在を 1EUV として火星生成初期を 6EUV、その中間を 3EUV と想定している。紫外
線強度別に算出された火星熱圏の温度分布を図 1.3 に、ピックアップイオン 1 個がはじき出す粒子の数を表 1-1 に示す。これらによると、紫外線強度が増す、すなわ
ち熱圏の温度が上がるほど大気の散逸量は全体的に増え、中でも O 原子の散逸する割合が増す。これは温度が高いほど exobase より上に存在する O 原子の量が増えるためである。つまり、非熱的散逸においても熱圏の温度が大気の散逸量や散逸成分に大きく影響していると考えられる。このように、熱圏の温度構造を理解することは大気散逸の理解につながり、それは惑星の気候を理解する上で重要となる。
地球熱圏の温度構造を決めるエネルギー収支においては、太陽紫外線による加熱と分子熱伝導による冷却、それに大規模な大気循環による熱輸送が支配的とされてきたが、対流圏で励起されて上方伝搬する重力波による直接的な加熱や冷却も大きく寄与することが Yigit and Medvedev (2009) などによって指摘されている (図 1.4)。
火星については、Bougher et al. (2009) が火星熱圏の大循環モデルによりエネルギー収支を計算した結果として図 1.5 を示している。ここでは惑星スケールの循環が太陽紫外線による加熱や熱伝導による冷却と比べて無視できない寄与をすることが示されているが、対流圏から伝搬する波動による直接的なエネルギー供給は考慮されていない。火星でも地球と同様に重力波によるエネルギー輸送が寄与する可能性があるが、それに関する研究はほとんどなされていない。
対流圏における擾乱の熱圏への影響を示唆する例として、Mars Global Surveyor(MGS) による大気密度観測の結果を図 1.6 に示す。これは高度 130-160km での密度変化を示しており、どの高度においても Ls=224.4 で密度が倍増しているが、この直前に下層でダストストームが生じている。この結果から下層での局所的な大気現象と上層での密度変化には何らかの関係があると考えられる。下層と上層をつなげる
プロセスとしては、子午面循環の変化や重力波によるエネルギー輸送が考えられる。

3

図 1.3: 紫外線強度別に算出された火星熱圏の温度分布 (Leblanc & Johnson, 2002)

表 1-1: O+イオン1個がスパッタリングで散逸させる粒子数 (Leblanc & Johnson,2002)

図 1.4: 地球における重力波の鉛直伝搬 (Yigit and Medvedev, 2009)4図 1.5:

火星の熱圏大循環モデルによる赤道域の昼側におけるエネルギー収支(Bougher et al., 2009)

図 1.6: MGS 搭載の加速度計で計測された火星熱圏の大気密度の時間変化 (Withers

and Pratt, 2013)。 Ls = 224.4 で密度の倍増が生じている。

5

1.3

火星大気における重力波

火星においても地球と同様に重力波と思われる波動が観測されている (図 1.7, 図1.8)。図 1.7 は MGS による電波掩蔽観測の結果である。3 次関数でフィッティングされた背景場の温度分布 (左図実線) からのずれをプロットしたものが右図になる。鉛
直波長が 10 km 程度の波の形が見えており、重力波と解釈したうえで振幅の季節・緯度依存性などが調べられている。ただし電波掩蔽による観測は高度およそ 40 km
以下に限定されている。
それより高い高度での重力波の観測は限られている。図 1.8 は Mars Odyssey による Aerobraking 時の大気密度計測の結果であり、高度 100-130 km において波動と思われる密度振動が検出されている (Fritts et al., 2006)。水平スケールが 20-200 km
であることを考慮すると、この波動は重力波であると考えられる。このように水平波長が比較的短い重力波は、電波掩蔽では電波経路に沿った 200 km 程度の積分効
果のために均されてしまい、観測することができない。これらの重力波の励起源として Fritts et al. (2006) は山岳波の可能性を挙げているが、励起過程の定量的な考察はされておらず、実際にどのような波が観測結果を説明しうるのかはわかっていない。

Walterscheid et al. (2013) は火星大気中を鉛直伝搬する重力波が大気に及ぼす影
響を線形波動解により算出している。図 1.9 は前述の Aerobraking 時の大気密度計
測結果 (図 1.8) のデータを参考とし、適当な水平波長と周期と密度振幅を与えて波動による大気の加熱率を計算した結果である。高度 170 km 付近と高度 300 km 付
近の 2 カ所に極大が見られるが、前者は顕熱フラックスの発散によるもの、後者は波の力学エネルギーの粘性散逸によるものである。いずれも重力波が熱圏大気において大きな分子粘性のために減衰することでもたらされる。この結果は、太陽紫外線による加熱と比較しても無視できない大きさの加熱を熱圏にもたらしうることを示している。これが現実に起きているのであれば、熱圏のエネルギー収支に考慮されるべきである。しかし、Walterscheid et al. (2013) は、水平波長と振幅について

Aerobraking による観測結果と矛盾しないものを与えたうえで、周期については根
拠なく仮定したものである。すなわち、観測された波がどのように励起されたものであるのかも、それらの波の特性が Walterscheid らが仮定したようなものであるのかも、現時点では不明である。現実的な考察を行うには波動の励起過程から考える必要がある。
重力波の励起源としては、山岳波に加えて地表付近で生じる対流も考えられる。
数値実験によれば、火星の大気は熱容量が小さいため、太陽光による加熱により風6

速が 10 m s−1 を越える激しい対流を生じうる (Odaka et al., 1998)。そのような対流は重力波を効果的に励起することが予想される。対流の形態としては、地表面が太陽光により加熱されることによって広い範囲でほぼ一様に発達する対流が考えられる。それに加えて、火星特有の対流として、大気中に浮遊するダストが太陽光により加熱されることで生じるものが考えられる。たとえば、Spiga et al. (2013) では、Mars Climate Sounder によって観測された高度数十キロに分布する薄いダスト
層 (図 1.10) の生成メカニズムを説明するものとして ロケットダストストーム を提案している。これによれば、地表面付近に局所的に存在するダストが太陽光を吸収し、それが大気を温める。その結果対流が生じ、ダストが高度数十キロまで巻き上げられる (図 1.12, 図 1.11)。このような対流が下層で生じているならば、そこから波動が励起されると考えられる。

7

図 1.7: Mars Global Surveyor による電波掩蔽観測 (Creasey et al., 2006)。 (左図) 丸
印が観測結果, 実線は３次関数でフィッティングした背景場の温度。 (右図) 背景場
からのずれ。

図 1.8: Mars Odyssey による高度 100-130 km での Aerobraking 時の大気密度計測(Fritts et al., 2013)。 滑らかに変化する成分を差し引いて擾乱成分を取りだしたもの。準周期的な構造は水平スケール 20-200 km, 振幅 5-50 %の重力波であると考えられる。

8

図 1.9: 線形解をもとに計算された鉛直伝搬する重力波による加熱率 (Walterschiedet al., 2013)

図 1.10: Mars Climate Sounder により観測された高度数十キロに分布する薄いダス
ト層 (Heavens et al., 2011)

図 1.11: 3 次元の領域モデルによる、太陽光を吸収した局所的なダスト塊よる大気加熱で生じた対流がダストが巻き上げる様子を計算 (Spiga et al., 2013)。

9図 1.12: 局所的なダスト塊の太陽光加熱により対流を生成する数値実験における、加熱率 (K/hour)。 加熱率の大きい場所には濃いダスト雲がある。(Spiga et al., 2013)

1.4

本研究の目的

以上のように、火星の熱圏のエネルギーバランスにおいて重力波によるエネルギー輸送が無視できない寄与をする可能性があり、これを定量的に評価することは重要課題である。1.3 節で述べたように、Walterscheid et al. (2013) が重力波による加熱率を計算しているが、彼らが与えたような波の特性がどれほど妥当であるのかも、
これらの波がどのように励起されるのかも分からない。そこで、私の研究では重力波が対流によって励起される可能性に注目し、高高度まで伝搬する波動の励起を対流計算により再現し、加速度計による密度変動計測の結果と比較する。加熱源としては、地表面が太陽光で加熱されることを想定した水平一様な場合と、局在した大気ダスト雲が太陽光加熱を受けることを想定した局所的な場合の 2 ケースを想定し、結果を比較する。また、対流計算から得られた波動のパラメータを用いて線形波動解を求め、波動が地表から熱圏高度まで鉛直伝搬して散逸するプロセスを調べる。
この手法によって、火星における重力波の励起源と伝搬過程を考察するとともに、
より現実的な加熱率の値を得る。
第 2 節では構築した対流モデルの説明と計算結果を、第 3 節では線形モデルの説
明と計算結果を示し、第 4 節で本研究のまとめを述べる。

10

2 次元対流実験

2

本研究では、火星の下層大気における対流を調べるために地球のメソ気象モデル

CReSS (Cloud Resolving Storm Simulator) Version 2.3 (Tsuboki and Sakakibara,2007) を用い、高度 0-300 km、水平 600 km の領域の 2 次元対流モデルを構築した。
このモデルにより, 高高度まで伝搬する波動の励起条件やメカニズムを調べ、より現実的な波動のパラメーターを得る。本章では、2.1 節でモデルの基本方程式について、2.2 節では計算領域と加熱強制について、2.3 節でその他のモデル設定について述べる。2.4 節で計算結果を示し、2.5 節で考察を述べる。

2.1

基本方程式

モデルの独立変数は空間の水平座標 x、鉛直座標 z 、時間 t である。これらの関数として定義される従属変数は、速度の水平成分 u、と鉛直成分 w、温位 θ、気圧 p、
密度 ρ である。これらの従属変数のうち圧力と温位と密度ついては静力学平衡

∂p
= −ρg
∂z

(2.1)

を満たす基準状態と、それからの偏差の関係が、

p = p + p′
θ = θ − θ′
ρ = ρ + ρ′のように与えられる。密度は、状態方程式

ρ=

p
Rd T

(2.2)

により求められる。ここで g は重力加速度、T は温度、Rd は乾燥空気の気体定数である。密度以外の従属変数は以下に示す時間発展方程式系で表現される。
運動方程式:

(
)
∂ρu
∂u
∂u
∂p′
= −ρ u+w+ T urb.u
−
∂t
∂x
∂z
∂x
)
(
)
( ′
∂w
∂ρw
∂w
∂p′
θ
p′
= −ρ u+w−
ρg− 2 + T urb.w
∂t
∂x
∂z
∂z
ρcs
θ

(2.3)
(2.4)

気圧偏差の方程式:

(
(
)
)
∂p′
∂p′
∂p′
∂u ∂w
ρc2 dθ
2=− u+w++ ρgw − ρcs
+ s
∂t
∂x
∂z
∂x
∂z
θ dt
11

(2.5)

ここで、cs は空気中の音速であり、

√
cs =

cp
Rd T
cv

(2.6)

である。cp は定圧比熱、cv は定積比熱である。
温位偏差の方程式:

(
)
∂θ
∂θ
∂ρθ
∂θ
= −ρ u+w− ρw+ T urb.θ
∂t
∂x
∂z
∂z

(2.7)

気圧の方程式以外に現れる T urb.ϕ はサブグリッドスケールの拡散項であり、以下にその定式化について簡単に述べる。
運動方程式中の拡散項は応力テンソル τij を用いて以下のように表現される。
)(∂τ11 ∂τ12+(2.8)
T urb.u =∂x∂z()∂τ21 ∂τ22
T urb.w =
+
(2.9)
∂x
∂z
ここで、応力テンソル τij は、分子拡散によるせん断応力とサブグリッドスケールの運動によるレイノルズ応力からなる。レイノルズ応力は平均量からの変動成分からなるので、以下のように拡散係数を用いた平均量の勾配拡散の形式にモデル化することができる。

()}∂u 1 ∂u ∂w
2ρvτ
−+∂x 3 ∂x
∂z
)
(
∂u ∂u
ρvτ
+
∂z
∂z
)
(
∂u ∂u
ρvτ
+
∂z
∂z
(
)}
{
∂w 1 ∂u ∂w
2ρvτ
−
+
∂z
3 ∂x
∂z
{

τ11 =
τ12 =
τ21 =
τ22 =

(2.10)(2.11)(2.12)(2.13)

ただし、vτ は運動量に関する渦拡散係数である。分子拡散は、対流が重力波を励起する下層大気においては渦拡散に比べてはるかに小さいため、ここでは考慮しない。
熱圏高度では分子拡散が無視できない大きさになるが、そのような高度での重力波の伝搬については 3 節において線形波動解を用いて別途調べる。
温位偏差の方程式中の渦拡散係数については以下のように表現される。
{ ()} { ()}∂∂θ∂∂θ
T urb.θ =
ρvH
+
ρvH
∂x
∂x
∂z
∂z

(2.14)

ただし、vH は温位に関する渦拡散係数である。ここでも分子拡散係数は渦拡散係数に比べて非常に小さいと考えて無視する。

12

本研究では、渦拡散係数を表現する手法として、乱流運動エネルギーを用いた 1.5
次クロージャーの渦拡散モデル (Klemp and Wikhelson, 1978) を使用する。このモデルでは、乱流運動が卓越している領域で拡散が作用することを踏まえ、渦拡散係
数が乱流運動エネルギーの平方根に比例すると仮定し、乱流運動エネルギーをモデ
ル内の物理量を用いた時間発展方程式で解く。
乱流運動エネルギーは、各速度成分について、平均流からの偏差に”を付して以下のように与えられる。

E=

)
1( 2
u” + w”22(2.15)

乱流運動エネルギーの時間発展方程式は乱流運動エネルギーの移流、浮力による生成、流れの変形速度による生成、拡散、散逸の効果を含め、以下のように与えられる。
(
)
∂E
∂E
∂E
gvH ∂θ
= −
u
+w
−
∂t
∂x
∂z
θ ∂z
{( )
(
)2 }
)2
(
)
(
2
∂u
∂w
∂u ∂w
∂u ∂w
2
+ 2vE
+
+ vE
+
− E
+
∂x
∂z
∂z
∂x
3
∂x
∂z
(
)
(
) ( )
3
∂
∂E
∂
∂E
ce
+
E2
(2.16)
vE
+
vE
−
∂x
∂x
∂z
∂z
l
ここで、vE は乱流エネルギーに対する渦拡散係数である。散逸項の係数 ce は乱流エネルギーが熱エネルギーへと変わる程度を表現する量で、CReSS では以下のように与えられる。

ce = 0.93(2.17)

このようにして求められた乱流運動エネルギーを用いて、運動量に対する渦拡散係
数 vτ は以下のように与えられる。
1vτ = 0.1E 2 l
ここで、l は混合長スケールで、

 ∆s,
l=
 分 (∆s, ls ) ,(2.18)

不安定または中立の場合
安定な場合(2.19)のように与えられる。ただし、
1(2.20)∆s = (∆x∆z) 2
1

ls = 0.76E 2

13

g ∂θ
θ ∂z

− 12

(2.21)

である。ここで、∆x,∆z はそれぞれ x,z 方向の格子間隔である。また、温位に対する渦拡散係数 vH と乱流運動エネルギーに対する渦拡散係数 vE は

vτ
vH
vE

=

1

2l
1+
∆s
= 2vτ

(= P r)(2.22)(2.23)のように求める。

2.2

計算領域と加熱強制

加熱強制に関して、地表面の太陽光加熱を想定した一様加熱とダスト雲の太陽光加熱を想定した局所加熱の２ケースを考える。現実には完全に水平一様な加熱は存
在せず、またダスト雲の加熱のみに限られるわけでもない。つまりこれらの実験は両極端のケースであり、実際にはこれらの中間的な状況が一般的であると考えられる。
高度方向には 0-200 km を物理領域とするが、その上に波の反射を防ぐ厚さ 100

km のスポンジ層を設けるので、実際の計算領域は高度 0-300 km である。水平方向の領域は 600 km であり、格子間隔は x 方向に ∆x = 1 km、z 方向に ∆z = 0.5 kmとした。中立成層を高度 0-5 km に設定したが、これは図 2.1 に示した電波掩蔽による昼側の観測結果において、中立成層が地表からおおよそ 3-7 km の高度まで存在していることを参考としている。この中立成層の範囲が日中の対流境界層に対応すると考える。
二種類の加熱強制に関して以下で述べる。

2.2.1

一様加熱

太陽光で加熱された地表面から大気に熱エネルギーが渡される過程を計算する代わりに、高度 0-1 km の大気を水平方向一様に加熱する (図 2.2)。加熱率は 10 K/hour としているが、これは境界層における夜間冷却の値から見積もったものである。Haberle

et al. (1993) の境界層モデリングによれば、日中に境界層が到達する高度領域における夜間の放射冷却は 50 K/day 程度である。この冷却により夜間に高度 0-5 km の境
界層全体から失われる熱は、日中に地表から地表近くの大気に与えられると考えられるので、この熱を 0-1 km の大気層に与えることにすると、その加熱率は 50(K/day)× 5(km)/1(km) = 250(K/day)∼10(K/hour) となる。

14

2.2.2

局所加熱

局所加熱は 1.3 節で紹介したロケットダストストームを想定したものである。加熱領域は水平領域中央に高度 0-3 km、幅 60 km の範囲で与え、加熱時間は 120 分とした (図 2.3)。ただし、高度方向、水平方向、時間方向いずれも次のような cos 関数
で与えている。

|x − xcen | < xwidth /2, |z − zcen | < zwidth /2, |t − tcen | < twidth /2 のとき
q = Qρ0 cos

2π (z − zcen )2π (x − xcen )2π (t − tcen )cos
cos
2zwidth2xwidth
2twidth

それ以外のときq = 0

Q は加熱の極大値、ρ0 は高度 0 km における密度、zcen = 1.5 km, xcen ＝ 300 km は
加熱領域中央の座標、zwidth = 3 km, xwidth = 60 km は加熱領域の厚さと幅を表す。

tcen = 60 分 は加熱が最大となる時刻で、twidth = 120 分 が加熱時間を表す。加熱率の最大値としては、Spigaet al. (2013) の実験の加熱率極大値 20 K/hour を参考に、標準実験として Q=10 K/hour、感度解析のために Q=5 K/hour の 2 ケースを採用した。
加熱率の変化は、立ち上がりは日の出とともにダストが加熱され始めることを意
識し、終息はダストが上昇して地表付近からなくなることを意識している。Spiga etal. (2013) の実験では 1-2 時間の時間スケールでダストが移動していることから、上述の加熱時間を採用した。

15

図 2.1: MGS による電波掩蔽で観測された下層大気の気温と温位の高度分布 (Hinson

et al. 2008)

図 2.2: 一様加熱実験の計算領域と加熱領域

図 2.3: 局所加熱実験の計算領域と加熱領域

16

2.3

その他の計算設定

2.3.1

初期値

CReSS に与える物理パラメーターは各高度での初期温位、下端気圧、下端密度、
定圧比熱、気体定数、重力加速度である。各高度での初期温位としては、図 2.4(左)で示した火星の温度分布を以下の式で近似し、そこで得られた温度の値を温位に変換して代入している。ここで T の単位は K、z の単位は m である。近似した結果を
図 2.4(右) に示す。

1
T (z) = 100 × tanh2(z − 125000
20000

)

(

z )
+ 185 + 70 × exp −30000(2.24)

その他のパラメーターを表 2.1 に示す。下端気圧と下端密度は計算領域下端の高度０ km における値である。計算開始時においては最大値 0.01K のランダムな初期
温位擾乱を全領域に与えた。

図 2.4: (左) 火星熱圏の温度分布 (Bougher et al. 1988)。(右) 計算に用いた温度分布。

17

表 2.1: CReSS に与えた物理パラメーター
パラメータ

2.3.2

値

下端気圧

600 [Pa]

下端密度

0.013 [km m−3 ]

定圧比熱

904 [J kg−1 K−1 ]

気体定数191 [J kg−1 K−1 ]重力加速度3.72 [m s−2 ]境界条件

側面は、一様加熱では周期境界条件、局所加熱では放射境界条件とした。上下はともに固定端としているが、上には波の反射を防ぐために厚さ 100 km のスポンジ層を与えている。スポンジ層では、物理量を初期値へ緩和させることで計算領域内に生じる波を吸収することを想定し、次のように基礎方程式の付加的な項として与
えられる。

∂ρϕ
= … − γ (z) ρ (ϕ − ϕ0 )∂t

(2.25)

ここで、ϕ は任意の予報変数、ϕ0 は当該予報変数の初期値、γ (z) は計算領域の減衰係数である。上端スポンジ層の減衰定数の大きさは、スポンジ層の下端高度を zlow 、
上端高度を ztop とすると z ≥ zlow に対して以下のように与えられる。
{()}
z − zlowγ (z) = α 1 − cos πztop ―zlow(2.26)

ここで、α は減衰定数で時間の逆数の次元を持ち、本研究では 0.003 s−1 とした。

2.3.3

時間刻みと数値粘性

CReSS では格子点の構造と変数の配置は水平・鉛直ともにスタッガードを用いている。速度の水平成分 u は鉛直方向にはスカラー量 (気圧偏差、温位偏差) と同じ
高度で、水平方向には半格子ずらして配置される。鉛直成分 z は水平方向にはスカラーと同じ位置で、鉛直方向には半格子ずらして配置される。CReSS で解く圧縮方
程式系は音波を解として含むため、時間積分をすべて陽解法で行おうとすると最も位相速度の早い音波の速度が CFL 条件を規定し、時間間隔が非常に短くなってしまう。この音波による制約を回避するために、モード別時間分割法を用いた。この方18

法は、音波に関する項についてのみ短い時間間隔で、それ以外の項を長い時間間隔で時間積分を行うものである。長い時間間隔の時間積分は、リープフロッグスキー
ムに Asselin のタイムフィルターを併用することで、奇数タイムステップと偶数タイムステップの間で起こる解の分離を抑え、短い時間間隔の時間積分は陽解法の時間
前進・空間後退スキームを用い、鉛直についてのみ陰解法の Crank-Nicolson スキー
ムを用いている。短い時間間隔の陽解法においては、運動方程式に現れる圧力項に
人為的な発散項を加えることで音波を抑えている。本研究では、局所加熱計算において長い時間ステップを 2.5 s、短い時間ステップを 0.25 s とし、一様加熱計算において長い時間ステップを 0.5 s、短い時間ステップを 0.05 s として計算を行った。
また、モデル内の移流項を中心差分で表現するにあたって、計算不安定を避けるために、人工的な 4 次の数値粘性を採用する。数値粘性項は、任意の物理量を ϕ として次のように表される。

∂ 4 (ρϕ)∂ 4 (ρϕ)Dif f.ϕ = −vh
− vv
∂x4
∂z 4

(2.27)

ここで、vh ,vv はそれぞれ水平方向と鉛直方向の粘性係数を表し、次のように定義さ
れる。

∆x4
∆t
∆z 4
vv = αv∆tv h = αh(2.28)(2.29)

ただし、αh ,αv は無次元量で、値は以下のように設定される。

αh = αv = 0.001

2.4

結果

2.4.1

一様加熱

(2.30)

対流が生じる下層 (高度 0-20 km) と重力波が伝搬する上層 (高度 50-200 km) に分けて結果を示す。図 2.5 は上層での鉛直速度、図 2.6, 図 2.7 は下層での鉛直速度と温位の変動成分を示している。時間はどれも加熱開始から 140 分後であり、暖色が正の値、寒色が負の値を表している。下層では高度方向・水平方向ともにおよそ 5km のスケールを持つ対流セルが並んで生じている。対流の中では上昇流と高温域、下降流と低温域がだいたい対応しており、対流によって熱が上向きに運ばれている。
対流のすぐ上の安定層に重力波が励起されている様子が観察される。上層では多くの対流から生じた波動が重なり合っているため綺麗な波形は見えていない。しかし、

19

図 2.8 に示したホフメラー図を見ると、上昇流と下降流が交互に現れており、波の伝搬が確認できる。ここで、エネルギーは上方へ位相は下方へ伝搬している様子が見られるが、これは群速度と位相速度の鉛直成分が逆向きになるという内部重力波の特徴に一致する。したがって、ここで見えている波動は内部重力波であると推定できる。
上層における波動の振幅は時間とともに大きくなり、計算開始からおよそ 140 分
あたりで領域上部において波が崩れ、計算が破綻する。このときには波の影響で等
温位面が大きく波打って対流不安定な領域が作られており (図 2.9)、これが原因と考えられる。波の振幅が飽和によって決まる上限値に達していると考えられるため、以下ではこの時刻における結果をもとに考察を行う。
図 2.10 には、高度 130 km における密度変動の水平構造を示す。これは Fritts etal.(2013) の観測で示されている火星の熱圏大気中の密度変動 (図 1.8) と対応している。観測された密度変動は振幅が 5-50%、波長が数十 km〜100 km 程度のものが卓越していたが、計算結果はこれらとおおよそ一致している。このことから、実際に火星大気中に存在しうる波動を対流計算により再現できていると言える。
図 2.5 から、卓越する重力波の水平波長は〜150 km、図 2.8 から、鉛直波長は〜

100 km、周期は〜20 分である。

20

図 2.5: 一様加熱における上層での鉛直速度。 計算開始から 140 分後の結果。 図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

図 2.6: 一様加熱おける下層での鉛直速度。計算開始から 140 分後の結果。 図中のInterval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

図 2.7: 一様加熱おける下層での温位の変動成分。計算開始から 140 分後の結果。 図中の Interval は等高線間隔 (K) を、Min と Max は最大値と最小値 (K) を表す。

21

図 2.8: 一様加熱での x = 300 km における高度 25-200 km の鉛直速度のホフメラー
図。図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

図 2.9: 計算開始から 140 分後の等温位面。

図 2.10: 一様加熱での高度 130 km における密度変動。 計算開始 120 分後の結果。

22

2.4.2局所加熱

まずは対流生成の様子を見るために、下層 (高度 0-20 km) の結果を示す。図 2.11,
図 2.13 は、加熱率 10 K/h における温位の擾乱と鉛直速度である。計算開始からの経過時間は 30 分, 60 分, 90 分, 120 分で、暖色が正の値、寒色が負の値を表している。温位擾乱を見ると、中央で温度の上昇が確認できる。またその時の鉛直速度を見ると中央で上昇流が発生し、それに伴って下降流が両側に生じており、対流が生じているのがわかる。時間が経過すると上昇流と下降流の構造が細かくなる様子が確認できる。
また、加熱率 5 K/h での鉛直速度を図 2.14 に示す。10 K/h の場合に比べて、対流の生成が遅くなっているのがわかる。

23

図 2.11: 加熱率 10 K/hour における下層での温位擾乱。計算開始から 30 分, 60 分,

90 分, 120 分後の結果。 図中の Interval は等高線間隔 (K) を、Min と Max は最大値と最小値 (K) を表す。

図 2.12: 加熱率 10 K/hour における下層での鉛直速度。計算開始から 30 分, 60 分,

90 分, 120 分後の結果。 図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

24

図 2.13: 加熱率 5 K/hour における下層での温位擾乱。計算開始から 30 分, 60 分, 90
分, 120 分後の結果。 図中の Interval は等高線間隔 (K) を、Min と Max は最大値と最小値 (K) を表す。

図 2.14: 加熱率 5 K/hour における下層での鉛直速度。計算開始から 30 分, 60 分, 90
分, 120 分後の結果。 図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

25

次に波動伝播の様子を見るために、上層 (高度 50-200 km) の結果を示す。図 2.15,
図 2.16 は、加熱率 10 K/hour と 5 K/hour における鉛直速度で、計算開始からの経
過時間は、10 K/hour については 75 分, 150 分, 225 分、5 K/hour については 75 分,

150 分, 320 分である。いずれについても上昇流と下降流が交互に現れており、明らかに波の構造が見えている。領域中央下部の対流域から放射状に波面が伸びており、
群速度が波面に対しておよそ平行方向である重力波が励起されていることがわかる。

75 分と 150 分における結果を 2 つの加熱率の間で比較すると、同じ経過時間では 5
K/hour よりも 10 K/hour のほうが 1 桁くらい振幅が大きいことがわかる。水平波
長、鉛直波長は時間とともに短くなり、水平波長はおよそ 20-150 km、鉛直波長はおよそ 20-100 km である。振幅は時間とともに大きくなり、10 K/hour では計算開始からおよそ 225 分で領域上部において波が崩れ (図 2.17)、このときには波の影響で等温位面が大きく波打って対流不安定な領域が作られており、これが砕波の原因と考えられる。すなわち、波の振幅は飽和によって決まる上限値に達していると考えられる。この時点の振幅は 10 K/hour と 5 K/hour の間で 2 倍程度の違いしかない。

26

図 2.15: 加熱率 10 K/hour における上層での鉛直速度。計算開始から 75 分, 150 分,

225 分後の結果。 図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

図 2.16: 加熱率 5 K/hour における上層での鉛直速度。計算開始から 75 分, 150 分,

320 分後の結果。 図中の Interval は等高線間隔 (m/s) を、Min と Max は最大値と最小値 (m/s) を表す。

27

図 2.17: 加熱率 10 K/hour における計算開始から 225 分後の等温位面。

図 2.18, 図 2.19 には、x = 100 km における鉛直速度のホフメラー図を示した。一様加熱計算の結果と同様にエネルギーは上方へ、位相は下方へ移動している様子が見られることから、ここで見えている波動は内部重力波であると推定できる。また、この図から読み取れる波動の周期は 10-30 分で、加熱率による違いははっきりしない。図 2.20 には、加熱率 10 K/h での高度 130 km における密度変動の水平構造を示す。振幅はおよそ 1-15 %、水平波長はおよそ 30-60 km であり、この結果もやはり、

Fritts et al.(2013) により観測された波の特性 (図 1.8) の範囲内である。

28

図 2.18: 加熱率 10 K/hour での x = 100 km における鉛直速度のホフメラー図。

図 2.19: 加熱率 5 K/hour での x = 100 km における鉛直速度のホフメラー図。

図 2.20: 加熱率 10 K/h での高度 130 km における密度変動。 計算開始 200 分後の結
果。

29

2.5

考察

一様加熱実験では、上層ほど波の構造が大きくなっている。初期に設定した厚さ5 km の中立層の中に風速が 15 m s−1 に達する強い対流が発生した。火星大気における過去の数値実験では太陽光による加熱により生じる対流の風速は 10 m/s くらいであり (Odaka et al., 1998)、今回の結果は先行研究と整合的である。個々の対流
セルは高さ・幅ともに 5 km 程であり、重力波のスケールとは大きく異なる。図 2.6を見ると、対流セルの高さが領域内で均一ではなく数十 km〜100 km のスケールで変化しており、このような対流セルの集団の大きな構造が重力波の水平波長を決めている可能性がある。高度およそ 150 km 以上で等位相面の傾きが急になる様子が見られるが、これは砕波により激しい鉛直混合が起こることに加え、スポンジ層の影響で波が減衰がする効果により位相の高度変化が少なくなることが理由として考えられる。
局所加実験では、波長が高度や時間とともに変化している。これは、下層で生成する対流が非線形性の効果により時間とともに細かくなっていき、励起される波動の水平スケールもそれに伴って大きくなっていくことである程度説明できると考えられる。つまり、伝搬しながら波長が変化しているのではなく、高高度の波ほど初期に生成された水平スケールな大きな波が見えているので波長が長いということになる。実際に、下層で見えている波の水平波長, 鉛直波長, 周期を読み取り、ある時間後の位置を算出して、そこに見えている波の水平波長, 鉛直波長, 周期を測るとだ
いたい同じ値になる。
局所加熱実験における重力波の振幅の大きさは時間とともに増大するが、波が崩れるところで頭打ちになる。図 2.15 における 225 分、図 2.16 における 320 分はそのような状態である。下層での対流生成にかかる時間の違いと励起される波動の伝搬速度の違いのため、砕波する前には同じ時刻で振幅を比べると 10K/h のほうが 5K/hに比べて 1 桁程度大きい。しかし、砕波が起こって振幅が頭打ちになると両者の振幅に大きな違いはなくなる。よって、最終的な振幅の大きさは加熱の大きさにはあまり依らないことがわかる。
一様加熱と局所加熱は両極端の例であり、一般的には両者の中間的な状況が多い
ものと考えられる。しかし、今回の実験ではどちらのケースにおいても観測とそれほど違わない大きさの水平波長と密度振幅を持つ波が見られた。ということは、その中間的な状況においてもそのような波が励起されると考えるのは自然である。今回の結果から、観測された密度変動には対流によって励起された重力波による寄与も含まれている可能性が示唆される。

30

卓越する周期は数十分程度の比較的短いものであった。そのために鉛直波長は数十 km と長く、これは重力波が減衰せずに高い高度まで到達する可能性を示唆する。
また、ここで見られた波は水平波長が 20-200 km 程であった。このような水平波長の短い波は電波掩蔽で観測することができないために観測データが少なく、惑星大
気研究においてはほとんど研究の対象とされてこなかった。しかし計算結果は熱圏で観測された波動の水平構造をよく説明することから、このような波が上層大気に与える影響を調べる必要がある。

31

重力波の線形波動解

3

重力波の鉛直伝搬プロセスとそれが散逸することによる加熱率を調べるために、波の強制項を加えた 2 次元のプリミティブ方程式を線形化して用いた。線形化したの
ち離散化することで方程式系は代数方程式となり、水平方向と時間方向に振動解を仮定したうえで鉛直構造を行列演算によって解くことができるようになる。与えるパラメータは水平波長と周期であり、対流実験から得られた値をもとに計算を行った。
本章では、2.1 節でモデルの基本方程式について、2.2 節でその解法、2.3 節で計算
領域と波動強制について、2.4 節で分子粘性や境界条件など残りの計算設定について説明し、2.5 節,2.6 節で結果と考察を述べる。

3.1

基本方程式

ここでは full-wave モデル (Hickey et al., 2000; Schubert et al., 2003) を用いている。波の反射を考慮しない WKB モデルに対して、このモデルでは波の反射も計算に含まれているため、より現実的な結果が得られると考えられる。独立変数は水平
座標 x、鉛直座標 z 、時間 t として定義され、従属変数は速度の水平成分 u、鉛直成
分 w、圧力 p、温度 T 、密度 ρ である。ここでは以下の方程式系を考える。
運動方程式:dv+ ∇p − ρg − ∇ · σm + ρKR v = 0
dt

(3.1)

dρ
+ ρ∇ · v = 0
dt

(3.2)

d (cv T )+ p∇ · v − σm : ∇v − ∇ · (km ∇T ) + cv ρKR T = q
dt

(3.3)

ρ
連続の式:

熱力学の式:

ρ

状態方程式:

p = ρRT(3.4)

ここで、v = (v, w) は速度ベクトル、σm は分子粘性応力テンソル、KR はレイリー
摩擦係数、km は分子熱伝率、cv は定積比熱、Q は非断熱加熱である。分子粘性応力テンソルはクロネッカーのデルタ δij を用いて以下で与えられる。
(
)
∂vi
∂vj
2
σij = µm+
− δij ∇ · v∂xj ∂xi 3µm は分子粘性係数である。
32

(3.5)

3.2

解法

まず方程式を線形化する。各変数を a = a + a′ のように、基本場と擾乱成分に分けて各方程式に代入する。ただし、背景風は考えないため u = 0, w = 0 とし、擾乱
の 2 次以上の項は微小量として無視する。また水平座標と時間について振動解を仮定して以下の形で書くことができる。

ã′ = a′ (z) exp(ikx − iωt)

(3.6)

k は水平波数、ω は振動数である。これにより水平微分と時間微分は以下のように置き換えられる。

∂
→ ik
∂x

,

∂
→ − iω
∂t(3.7)

さらに、大気密度の指数関数的な高度変化による振幅の大きな高度変化が解に与える誤差を抑えるために、以下のような変数変換を行う。

u′ = e 2H ũ′zw′ = e
p

′

= e

z
2H

(3.8)

w̃′

z
− 2H

p̃

(3.9)
′

(3.10)

T′ = ez
2HT̃ ′(3.11)

q′ = ez2Hq̃ ′(3.12)

ここで、H は平均的なスケールハイトである。擾乱をを加えた方程式から基本場が満たす方程式を引くことにより、擾乱が満たす方程式が以下のように導かれる。
)
(
)
)
( ′
( ′
ũ
1 dũ′ d2 ũ′dµm ũ′
dũ′
ik
w̃
dw̃′
−µm
+
+ 2 −
+
+ µm
+
4H 2 H dz
dz{dz
2Hdz } 3
2H
dz
z
dµ
4
m ′
w̃ − ike− H p̃′ = 0 (3.13)+ iρ (ω − iKR ) + µm k 2 ũ′ + ik
3
dz
( ′
)
( ′
(
)
)
w̃
1 dw̃′ d2 w̃′4 dµm w̃′
dw̃′
ũ
dũ′
4
ik
+
+
−
+
+
− µm
+ µm
3
4H 2 H dz
dz 2
3 dz
2H
dz
3
2H
dz
( ′
)
′
′
dp̃
2ik dµm ′ − z p̃ dp − z p̃gρũ −e H+e H+
− T̃ ′ = 0 (3.14)+{iρ (ω−iKR )+µm k 2 }w̃′ −
3 dz
p dz
2H dz
T

dw̃′
d lnρ ′ iω ′ iω − z ′w̃′+− ikũ′ +w̃ − T̃ + e H p̃ = 0 (3.15)
2H
dz
dz
p
T
(
) (
(
)
)
(
)
T̃ ′
1 dT̃ ′ d2 T̃ ′w̃′ dw̃′d lnT
dkm
T̃ ′ dT̃ ′km++v 2 −p++km+++ikpũ′4H 2 H dz
dz
2H dz
dz
dz
2H dz
[
{
(
)
}]
dT
1 dT d lnT d lnkm d lnT
2−cvρ w̃+ −icv ρ(ω−iKR)+km+−−k T̃ ′ = q̃ ′ (3.16)
dz
dz
dz
dz
T dz

33

次に、鉛直 N 層の格子をとって方程式の離散化を行う。鉛直微分の項は有限差分で表し、式 (3.13) から式 (3.16) は第 k 層において
(
)
(
)
ũk
dµm ũkũk+1 − ũk−1
1 ũk+1 − ũk−1 ũk+1 − 2ũk + ũk−1
−
−µm
+
+
+
4H 2 H
2∆z (
2H
2∆z}
(∆z)2 ) { dzikw̃kw̃k+1 − w̃k−1
4+ µm++ iρ (ω − iKR ) + µm k 2 ũk
3
2H
2∆z
3
z
dµm
+ik
w̃k − ike− H p̃k = 0 (3.17)dz()
4
1 w̃k+1 − w̃k−1 w̃k+1 − 2w̃k + w̃k−1
w̃k
− µm
+
+
2
2
3
4H
H(
2∆z
) (∆z) ()
4 dµm w̃kw̃k+1 − w̃k−1
ikũkũk+1 − ũk−1
−
+
+ µm
+
3 dz
2H
2∆z
3
2H
2∆z
z p̃k dp2dµmũk − e− H+{iρ (ω − iKR ) + µm k 2 }w̃k − ik
3 dz
p )dz
(
z
p̃k
p̃k+1 − p̃k−1
gρ
−H
+e
+
− T̃k = 0 (3.18)
2H
2∆z
T

w̃k
iω z
w̃k+1 − w̃k−1
d lnρ
iω
+
− ikũk +w̃k − T̃k + e− H p̃k = 0 (3.19)
2H
dz
dz
p
T
)
(
(
)
T̃k
w̃k
1 T̃k+1 − T̃k−1 T̃k+1 − 2T̃k + T̃k−1
w̃k+1 − T̃k−1
−p
km
+
+
+
4H 2 H
2∆z
2H
2∆z
(∆z)2
)
)(
(
T̃k
T̃k+1 − T̃k−1
dkm
d lnT
dT
+ km
+
+ ikpũk − cv ρ w̃k
+
dz
dz
2H
2∆z
dz
[
(
)
}]
{
1 dT
d lnT d lnkm d lnT2T̃k = q̃k (3.20)+ −icv ρ (ω − iKR ) + km
+
−
−k
dz
dz
dz
T dz
と書き換えられる。それぞれの方程式に対し k = 1, …, N の N 個の方程式ができるので、合計 4N 個の方程式ができる。また、未知数の数も各変数に対し k = 1, …

, N の N 個で合計 4N 個となるので、この連立方程式は解くことができる。プログ
ラム上では

AX = B(3.21)′, p′1 , …, p′N , T1′ , …, TN′ ]X T = [u′1 , …, u′N , w1′ , …, wN(3.22)の行列式を解く。X はであり、A は方程式の係数を成分とした 4N ×4N の行列、B は強制項

B T = [0, …, 0, 0, …, 0, 0, …, 0, q1 , …, qN ](3.23)である。
また、波動による大気の加熱率 Q[K/day] は以下の式で与えられる。

ρcp Q = ⟨σm ′ : ∇ v ′ ⟩ −ddp ⟨w′ ρ′ ⟩{cp ρ⟨w′ T ′ ⟩} + ⟨v ′ · ∇p′ ⟩ −
dz
dz ρ
34

(3.24)

第一項は運動エネルギーの粘性散逸 (viscos heating), 第二項は顕熱フラックスの発
散 (sensible heat ﬂux divergence), 第三項は圧力勾配がなす仕事（work done by thewave-induced pressure gradients), 第四項は質量の鉛直輸送がなす仕事 (work doneby the second-order wave-induced Eulerian drift) である。

3.3

計算領域と波動強制

高度 0-500 km を物理領域とするが、下端に波の反射を防ぐためのスポンジ層を厚さ 200 km で入れるため、計算領域は-200 km から 500 km となる。格子間隔は 1 kmとする。スポンジ層はレイリー摩擦係数として運動方程式と熱力学の式に含まれており、以下の式で与えられる。

KR = ω sin

2πz
4Hl

(3.25)

ここで、Hl はスポンジ層の厚さ (=200 km) である。上端に関しては分子拡散により波が十分に減衰するため、スポンジ層を設けていない。
また、波動は物理領域の下端近くにおいて加熱強制で励起させており、その振幅の高度分布は以下の式で与えている。

q ′ = cv ρq̂ exp−(z − zc )2
δz 2(3.26)

ここで、q̂ は加熱率の極大値、zc は加熱の振幅が最大となる高度 (=5 km)、δz は加熱を与える高度の広がり (=2 km) を表す。加熱率 q̂ の大きさは、高度 120 km での鉛直風振幅が対流実験の結果と合うように調整して与えられる。ただし、対流実験
においてこの高度で砕波が起こっている場合、それより高い高度でこの振幅を越えることはありえないため、線形解においては高度 120 km 以上での振幅の最大値が飽和振幅とおおよそ一致するように振幅を与える。

3.4

その他の計算設定

3.4.1

背景場

大気組成は CO2 と O を想定し、高度とともに O 原子の占める割合が増える様子を再現している (図 3.1 左)。背景場の温度分布は対流実験と同様に図 2.4(左) を参考にして与えている (図 3.1 右)。気圧は高度 0 km で 600 Pa とし、静水圧平衡しているとして高度とともに減少する。また、背景風は無いものとする。

35

図 3.1: 線形モデルで与えた平均分子質量分布 (左) と温度分布 (右)

3.4.2

分子粘性

大気中を伝搬する波は粘性によって散逸される。粘性には、大気の乱れによる渦粘性と分子の運動による分子粘性があり、分子の平均自由行程の長い高高度ほど分子粘性が強く効いてくる。火星大気中の渦拡散係数はよくわかっていないため、仮に地球における値を参考に渦拡散係数を 10 m2 /s として (Fukao et al., 1994) 鉛直波
長 60 km の波について渦拡散の時定数を計算するとおよそ 2.8 日であり、これは本研究で対象としている波動の周期 (数十分) に比べて十分長い。そのため渦拡散の影響はほとんど受けずに伝搬していくと考えられるので、このモデルでは渦粘性は考えず分子粘性の効果のみ考慮する。また、放射減衰も波を散逸させる要因の一つだが、
図 3.2 に示すように、この時定数も高度 120 km 付近で局所的に 1.2 時間程度である
がそれ以外の高度では数時間以上と周期に比べて長い (Eckermann et al., 2011)。完
全に無視できる大きさではないが、ほとんど効かないとしてこのモデルでは考慮しない。
分子粘性係数は、分子拡散係数 Dm を用いて

µm = ρDm(3.27)と書ける。分子拡散係数は、Chapman and Cowling (1970) から引用して以下のよ
うに与えた。
−5 1

Dm = 1.38×10

36

ρ

(

T
273

)0.933(3.28)

この高度分布を図 3.3 に示す。また、分子熱伝導率も分子粘性係数を用いて表され、

km = 1.64cv µm(3.29)と与えた。

図 3.2: 火星における放射減衰率 (day−1 ) の高度と鉛直波長への依存性 (Eckermann

et al., 2011)

図 3.3: 線形モデルで与えた分子拡散係数 (Chapman and Cowling, 1970)37

3.4.3境界条件

境界条件は上端と下端で振幅ゼロとした。ただし、下端には厚さ 200 km のスポンジ層を入れることで波の反射を防いでいる。上端では本来、放射境界条件が適切であるが、分子粘性の影響により波が十分に減衰するため上端の境界条件も振幅ゼ
ロとし、スポンジ層も設けていない。

3.5

結果

与える波のパラメータとしては、対流実験における高度 120 km 付近での典型的な水平波長・鉛直波長・周期・鉛直速度振幅を読み取って用いた。高度 120 km 以
下では分子粘性の影響が小さいので、波動励起実験の結果からこの高度において波のパラメータを読み取ることは妥当である。ただし振幅については、この高度で砕波が起こって振幅が上限値に達する場合、それより高い高度でこの振幅を越えることはありえないため、線形解においては高度 120 km 以上での振幅の最大値がこの上限値になるように振幅を与える。また、局所加熱実験で見られた水平波長・鉛直
波長・周期の組み合わせは線形解とは完全には一致しない。たとえば、線形解において水平波長と周期を局所加熱実験結果に合わせると、計算される鉛直波長は局所
加熱計算結果に比べて半分程度になる。この理由としては、局所加熱実験では波源から波面が放射状に広がっており、平面波の解とは異なる状況であることが考えられる。熱圏での伝搬特性においては鉛直伝搬速度が重要であるため、ここでは鉛直波長と周期を対流実験結果に合わせた。結果として局所加熱実験を想定した計算における水平波長は対流計算結果に比べて 2 倍程度となっている。パラメータの値を
表 3.1 に示す。Walterscheid et al. (2013) による火星重力波の線形解 (図 1.9) に比べ、ここで採用した波長と周期はかなり短い。
一様加熱実験においては、上空での波の振幅が飽和振幅に達して砕波が生じる直前の 140 分後の結果を用いた (Case 1)。局所加熱実験においては、砕波が始まるよ
り前の 150 分後の結果 (Case 2-1) と波が飽和して振幅がほぼ定常に達した 225 分後の結果 (Case 2-2) を用いた。各ケースで計算された鉛直風振幅と位相、波が背景大気にもたらす加熱率の高度分布を図 3.4 - 3.6 に示す。火星の大気では高度 250 km あたりが exobase となるので、300 km より上の計算結果は考えないことにして、0-300

km の結果を提示する。また、加熱率のグラフでは４成分の和を赤線で示しており、
最小値, 最大値とその時の高度を表 3.2 にまとめた。

Case 1 では高度約 150 km で振幅が最大となり、位相は高度約 150 km より上で変化が緩やかになっている。つまり鉛直方向に波として伝搬しにくくなっている。この38

ことは、これより上で振幅が減少していくことと整合的である。また振幅が最大となる高度のすぐ下とすぐ上で加熱率はそれぞれ極大, 極小となり、さらに上空での振幅の減衰過程において二つ目の極大が存在している。加熱率の各成分を見ると、一つ目 (高度約 140 km) は顕熱フラックスの発散によるもの、二つ目 (高度約 180 km)
は運動エネルギーの粘性散逸によるものであることがわかる。このような特徴は、Walterscheid et al. (2013) がやはり火星大気中の重力波の線形解で示したもの（図1.8）と定性的には同じである。
Case 2-1 では高度約 150 km で振幅が最大となり、高度約 150 km より上ではほぼ波として伝搬できていない。また Case 1 と同様に、振幅が最大となる前後で加熱率は極大, 極小となり、極小では負の値を取り冷却を引き起こしている。このような冷却は、重力波の減衰過程で鉛直方向に熱輸送が生じる結果としてもたらされる(Yigit and Medvedev, 2009)。振幅の減衰過程における二つ目の極大では、一つ目の極小と同じくらいの加熱をもたらしている。

Case 2-2 では高度約 120 km で振幅が最大となり、高度約 150 km より上では波として伝搬しにくくなっている。このケースでもやはり振幅が最大となる前後で加熱率は極大, 極小となっているが、極大は一つしか存在せず、Case 1, Case 2-1 に比べて大きな冷却をもたらしている。冷却部分では顕熱フラックスの発散が大きく寄与している。

表 3.1: 線形モデルで与えた波（鉛直波長は水平波長と周期を与えた結果として計算された値）
水平波長 鉛直波長

周期

水平位相速度

鉛直位相速度

鉛直風振幅

Case 1∗1

95 km

60 km

15 min105.6 m/s

66.7 m/s

50 m/sCase 2-1∗2

110 km

45 km

20 min91.7 m/s

37.5 m/s

4 m/s

Case 2-2∗3

54 km

20 km

20 min45.0 m/s

16.9 m/s

20 m/s*1: 一様加熱, 140 分後、*2: 局所加熱, 150 分後、*3: 局所加熱, 225 分後39

図 3.4: Case 1 の計算結果。左から鉛直風振幅、位相、加熱率. 加熱率の図では 4 成
分の和を赤線で示している。

図 3.5: Case 2-1 の計算結果。左から鉛直風振幅、位相、加熱率. 加熱率の図では 4成分の和を赤線で示している。

40

図 3.6: Case 2-2 の計算結果。左から鉛直風振幅、位相、加熱率. 加熱率の図では 4成分の和を赤線で示している。

表 3.2: 加熱率の最大値・最小値

3.6

最大値 , 高度

最小値 , 高度

Case 1

1943.1 K/day , 139 km

467.1 K/day , 161 km

Case 2-1

16.4 K/day , 132 km-1.0 K/day , 160 km

Case 2-2

554.9 K/day , 120 km-393.0 K/day , 131 km

考察

いずれのケースにおいても分子粘性の効果により振幅が減衰する高度領域において大気加熱が見られるが、その高度分布の形には違いがある。Csse 1 と Case 2-1 では、高度 130-140 km 付近に見られる顕熱フラックスの発散による加熱率の極大に加えて、さらに上空の高度 200-250 km 付近で運動エネルギーの粘性散逸による加熱が生じている。それに対して Case 2-2 では、高度 120-130 km において顕熱フラッ

41

クスの発散による加熱と冷却が見られるものの、その上空に運動エネルギーの粘性散逸による極大は生じない。Hickey et al. (2011) は地球の熱圏を伝搬する重力波の理論的考察により、水平位相速度が速いほど運動エネルギーの粘性散逸による加熱が顕著になると結論しているが、この結果はそれと整合的である。これは、近似的な分散関係 cz = ωcx /N （cz は鉛直群速度、ω は周波数、cx は水平位相速度、N は

Brunt-vaisalla 周波数）から理解される。周波数が今回のように同程度の場合には、
水平位相速度が大きいほど鉛直群速度が大きくなるため分子粘性の影響が小さくな
り、より高い高度まで到達したのち大きな分子粘性により散逸するからと考えられる。実際、Case 2-2 では高度 200 km 以上で振幅が非常に小さくなっているが、Case

1 と Case 2-1 では高度 300 km 付近まで振幅が顕著な大きさを保っている。
Case 2-1 と Case 2-2 を比較すると、Case 2-1 が Case 2-2 に比べてより高い高度で大気加熱をもたらす一方で、加熱率の大きさは Case 2-2 の方が一桁大きい。このよ
うな違いは、局所加熱実験の結果においてパラメータの値を読み取った時刻が異なることに起因する。つまり、対流発生から時間とともに波動が大気に影響を与える
高度とその加熱率が変わっていくことを示している。対流生成初期に励起された波
動は鉛直波長が長く、鉛直群速度が大きいためより高い高度まで伝搬できるが、振幅が小さいため加熱率は小さい。時間の経過とともに励起される波動の鉛直波長は短くなり、鉛直群速度が小さくなる一方で、飽和するまで振幅は増大し続けるので波動の到達高度は低く、加熱率は大きくなると考えられる。このような鉛直波長の時間変化は、水平波長が時間とともに小さくなることと対応している。このように、対流生成初期に励起された波ほど高い高度において運動エネルギーの粘性散逸による大気加熱を引き起こしやすいという結果が得られた。
加熱率の大きさは、どのケースにおいても紫外線による加熱など熱圏のエネルギー
収支において従来考慮されてきたプロセスと比べて無視できない大きさである。特に一様加熱実験を想定した Case 1 では 2000 K/day に達しており、現実的であるとは考えにくい大きさになっている。ただし、今回は対流計算も線形波動解も二次元の計算であり、三次元的に波が広がっていく場合には振幅の高度分布が大きく変わってくる可能性もある。大きな影響は無いとして考慮しなかった放射減衰の影響も再考の余地がある。また、今回考えなかった背景風を考慮すると、高高度まで伝搬できる波が制限されうるので、これも加熱率の分布に影響を与えると考えられる。

42

4

まとめ
本研究では、火星大気中の重力波が熱圏に与える影響を調べるために、2 次元対流実験により高高度まで伝搬する重力波の生成を考察し、そこで観察された重力波の鉛直伝搬特性について線形モデルを用いて解析した。
まず非静力学気象モデル CReSS を用いて 2 次元対流モデルを構築し、地表面が太陽光により暖められることを想定した一様加熱と、局所的に分布したダストが太陽光の熱を吸収し大気を暖めることを想定した局所加熱を考えた。一様加熱の加熱率
は Haberle et al. (1993) により見積もられた境界層での夜間放射冷却の値を参考として与えて、局所加熱の加熱率は Spiga et al. (2013) においてダストストームをシ
ミュレーションする際に算出された値を参考とした。いずれの実験においても高高度まで伝搬できる内部重力波が見られ、高度 100 km 程度より上では砕波が生じて振幅が飽和する様子が見られた。卓越する水平波長は 100 km 前後、鉛直波長は数十 km、周期は 20 分前後であった。このような波は速い鉛直伝搬速度を持つため、あまり減衰せず熱圏高度にまで達すると考えられる。熱圏下部における密度変動をFritts et al. (2013) の観測で示されている火星熱圏大気中の密度変動と比較したとこ
ろ、水平波長・密度振幅共によく一致していた。このことから、火星の熱圏で観測された密度変動には対流によって励起された重力波による寄与が含まれている可能性が高いと言える。
次に、対流実験の結果から熱圏下部における典型的な波のパラメーターを読み取
り、線形波動解により波の鉛直伝搬特性を調べた。線形解としては full-wave モデル(Hickey et al., 2000; Schubert et al., 2003) を用い、水平方向と時間方向に振動解を仮定したうえで鉛直構造を行列演算によって解いた。熱圏に到達した重力波は分子粘性により減衰し、大気加熱をもたらすが、加熱率の高度分布には 2 つの極大が見られた。一つは振幅が最大となる高度付近で生じる顕熱フラックスの発散によるものであり、もう一つはそれより高い高度で生じる運動エネルギーの粘性散逸によるものである。鉛直伝搬速度の速い波は振幅の減衰が緩やかであるために、顕熱フラックスの発散による加熱に加えて高い高度での運動エネルギーの粘性散逸による加熱が目立ち、一方で鉛直伝搬速度の小さい波は顕熱フラックスの発散による加熱
しか見えないという結果となった。局所加熱実験においては、対流発達の初期に励起される波は振幅が小さく鉛直伝搬速度は大きいという特徴を持つため熱圏の高高度まで到達して比較的小さな加熱率をもたらし、一方で対流発達の後期の波は振幅が大きく鉛直伝搬速度は小さいという特徴を持つため比較的低高度で大きな加熱率をもたらすという結果が得られた。

43

本研究では火星の熱圏まで伝搬する重力波の対流励起を初めて再現し、その結果
から波のパラメータを決定することによって加熱率の計算を行った。このようなアプローチは、Walterscheid et al. (2013) が恣意的に波動のパラメータを仮定して加熱率を評価したものとは全く異なるものである。結果として得られた加熱率は太陽紫外線によるものと比べても無視できない大きさとなった。今回の結果を踏まえると、今後は火星熱圏のエネルギー収支において重力波によるエネルギー輸送も考慮する必要があると言える。また、今回重要性が示されたような鉛直伝搬速度の速い
波は他の惑星においても励起されている可能性が高く、広く惑星熱圏の構造を理解するうえで重要な物理過程である可能性がある。

44

謝辞
本研究の遂行にあたり、熱心かつ丁寧にご指導いただきました指導教官である今村剛先生に心より御礼申し上げます。また、数値モデル CReSS の使用方法や研究方
針に関する助言を下さいました前島康光先生に御礼申し上げます。中村正人先生、
佐藤毅彦先生、阿部琢美先生にはセミナーを通じて研究内容への指摘や助言いただ
きました。心より感謝申し上げます。日々の研究において多くの助言、励ましを頂きました安藤紘基様を始めとする宇宙科学研究所・太陽系科学研究系・惑星大気・
電離グループの皆様に心より感謝申し上げます。STP セミナーを通じて有益な助言を下さいました宇宙プラズマグループの皆様、並びに東京大学大学院地球惑星科学
専攻、宇宙惑星科学講座太陽地球系物理コロキウムの皆様に感謝申し上げます。また、様々な面で研究生活を支えて下さいました宇宙科学研究所秘書である石川美穂
様に心より御礼申し上げます。

45

参考文献
[1] Bougher, S.W., McDunn T.M., Zoldak K.A. and Forbes J.M. : Solar cycle
variability of Mars dayside exospheric temperatures: model evaluation of underlying thermal balances, Geophys. Res. Lett., 36, 2009[2] Bougher, S.W., Dickinson R.E. : Mars, mesosphere and thermosphere 1. Global
mean heat buget and thermal structure, J. Geophys. Res., 93,7325-7337, 1988
[3] Chapman, S. and T. G. Cowling : The Mathematical Theory of Non-Uniform
Gasses, Cambridge Univ. Press, New York, 1970.
[4] Creasey, J.E. , J. M. Forbes, and D. P. Hinson : Global and seasonal distributionof gravity wave activity in Mars’ lower atmosphere derived from MGS radio
occultation data, Geophys. Res. Lett., 33, 2006
[5] Eckermann, S. D. , J. Mab, X. Zhu: Scale-dependent infrared radiative damping
rates on Mars and their role in the deposition of gravity-wave momentum ﬂux,
Icarus, 211 429- 442,2011[6] Fritts, D. C., L. Wang, and R. H. Tolson : Mean and gravity wave structures andvariability in the Mars upper atmosphere inferred from Mars Global Surveyor
and Mars Odyssey aerobraking densities, J. Geophys. Res., 111, 2006
[7] Heavens, N. G., M.I.Richardson , A.Kleinbhl ,D. M. Kass ,D. J. McCleese , W.
Abdou,J. L. Benson ,J. T. Schoﬁeld ,J. H. Shirley , and P. M.Wolkenberg : The
vertical distribution of dust in the Martian atmosphereduring northern springand summer: Observations by the Mars Climate Sounder and analysis of zonal
average vertical dust proﬁles, J. Geophys. Res.,116,2011[8] Hinson, D.P. , M.Ptzold , S.Tellmann ,B. Huslerd ,G.L.Tyler : The depth ofthe convective boundary layer on Mars Icarus, 198, 57-66, 2008
[9] Leblanc, F. and R. E.Johnson : Role of molecular species in pickup ion sputtering of the Martian atmosphere, J. Geophys. Res., 107(E2), 2002
[10] Schoﬁeld, J.T., J.R. Barnes, D. Crisp , R.M.Haberle, S.Larsen , J. A. Magalhaes , J. R. Murphy , A.Seiﬀ , G. Wilson : The Mars Pathﬁnder Atmospheric
Structure Investigation/Meteorolog, ISSN 0036-8075, 199746

[11] Spiga Aymeric , Julien Faure, Jean-Baptiste Madeleine, Anni Maattanen, andFrancois Forget : Rocket dust storms and detached dust layers in the Martian
atmosphere,Planets,2013[12] Odaka Masatsugu, Kensuke Nakajima, Shin-ichi Takehiro, Masaki Ishiwatari,
and Yoshi-Yuki Hayashi : A numerical study of the Martian atmospheric convection with a two-dimensional anelastic model, Earth Planets Space, 50, 431―437, 1998

[13] Yigit Erdal and Medvedev S. Alexander : Heating and cooling of the thermosphere by internal gravity waves, Geophys. Res. Lett., 36, 2009[14] Walterscheid, R. L., G. Schubert, and D. G. Brinkman : Acoustic waves in theupper mesosphere and lower thermosphere generated by deep tropical convection, J. Geophys. Res., 108(A11), 1392, 2013[15] Withers Paul , R.Pratt : An observational study of the response of the upper
atmosphere of Mars to lower atmospheric dust storms, Icarus, 225, 378-389,201347