マントル対流数値シミュレーション概論

1 マントル対流を記述する方程式の導出

(February 9, 2022)

まず、3次元デカルト座標系内での流体の運動を記述する一般的な方程式の形を 示す。 その後、マントル対流の性質を利用して方程式を簡単化する。

なお、この章の記述は [6, 13]を参考にした。

1.1 基礎方程式

流体の運動を記述する基本的な方程式は、質量、運動量、エネルギーの保存則で ある。 これに加え、状態方程式、及び輸送係数(熱伝導率や粘性率)を規定する方程式を 含めて、閉じた基礎方程式系を構成する。

1.1.1 質量保存則

ρt+xi(ρvi)=0 (1.1)

ここで、ρは流体の密度、𝒗=(v1,v2,v3)は流体の速度である。

よく知られているLagrange微分 (または実体微分)

DDt=t+vixi

を使って式(1.1)を書き直すと

1ρDρDt+vixi=0 (1.2)

を得る。 この式の左辺第1項は、ある一定の質量を持つ流体素片が占める体積の変化を表 している。 特に非圧縮(ρが一定)の場合にはさらに

vixi=0 (1.3)

と簡略化される。

1.1.2 運動量の保存則

ρDviDt=-pxi+τijxj+ρgi (1.4)

この式の左辺は単位体積の流体素片の質量(ρ)と加速度の積、即ち流体素片 の持つ運動量の時間変化を表す。

右辺第1項と第2項は流体素片にはたらく表面力(面を通してはたらく力)を表す。 このうち、第1項は流体の圧力pによる等方的な応力、第2項は偏差応力テンソ ル (deviatoric stress tensor)τijで表される非等方的な応力である。 流体素片にはたらくトルクのつりあいを考えると、偏差応力テンソルは対称テン ソル

τij=τji (1.5)

であることがわかる。

右辺第3項は流体素片にはたらく体積力を表す。 ここでは体積力は重力のみを考え、その他の力(電磁気的な力やコリオリ力など) は考えない。 後に示すように、マントル対流のような遅い流れでは、コリオリ力の影響は無視 できる。

1.1.3 構成方程式とNavier-Stokes方程式

物体の力学的応答、例えば作用する応力と変形速度を結びつける式を構成方程式 (constitutive equation)という。 運動量の保存則(1.4)に構成方程式を加えることで、流体の速度に関す る閉じた方程式を構成することができる。

流体中の速度場から速度勾配テンソル(velocity gradient tensor)Lij

Lij=vixj

のように定義する。これを用いて2つのテンソル𝑬𝑾

𝑬=12(𝑳+𝑳T)oreij=12(vixj+vjxi) (1.6)
𝑾=12(𝑳-𝑳T)orwij=12(vixj-vjxi)

のように定義する。 𝑾は回転速度テンソルと呼ばれ、剛体的な回転を表わしている。 一方、𝑬は歪速度テンソル(strain rate tensor)と呼ばれ、非剛体的な変 形を表わす。 この両者のうち、𝑾は(剛体回転を表わしているから)流体中の偏差応力に は寄与しない。

多くの流体では、歪速度テンソルeijと偏差応力テンソルτijの間 には線型な関係がある。 このような流体をNewton流体と呼ぶ。 さらに流体が等方的であるとすれば、eijτij

τij=2μeij+λekkδij (1.7)

なる関係式で結ばれる。 ここでδijはKroneckerのδであり、μは粘性率、λ は第2粘性率と呼ばれる。 また偏差応力テンソルτijの対角和(trace)をとり、空間の3方向に平均 すると、

τii3=eii(λ+23μ) (1.8)

となるが、ここで現れるλ+23μを体積粘性率 (bulk viscosity) kB と呼び、kBと書く。 体積粘性率kBは流体の体積変化に伴なう粘性応力の目安である。 多くの流体ではkBは非常に小さいことが知られているので、kB=0とする (Stokes近似)。 これより、Newton流体かつStokes流体の場合には、偏差応力テンソルの成分は

τij=2μeij-23μekkδij=μ(vixj+vjxi-23vkxkδij) (1.9)

となる。

式(1.9)を式(1.4)に代入すると、よく知られたNavier-Stokes 方程式

ρDviDt=-pxi+xj[μ(vixj+vjxi-23vkxkδij)]+ρgi (1.10)

を得る。 特に非圧縮を仮定すると、式(1.10)は

ρDviDt=-pxi+xj[μ(vixj+vjxi)]+ρgi (1.11)

と簡単化される。 さらに(マントル対流らしくはなくなるが)粘性率μが一定であるとすると

ρDviDt=-pxi+μ2vixj2+ρgi (1.12)

と表わされる。

1.1.4 状態方程式

熱力学的考察によれば、物質の状態は2つの状態変数(温度T、圧力pなど)に より定められる。 即ち、状態方程式は一般的には3つの変数を含んだ式、例えば

f(p,ρ,T)=0 (1.13)

のように表わされる。

状態方程式を具体的に指定する際には、熱膨張率αや等温圧縮率χT を用いるのが便利である。 これらは以下のように定義される。

α1V(VT)p=-1ρ(ρT)p (1.14)
χT-1V(Vp)T=1ρ(ρp)T (1.15)

ここで、V=1/ρは比体積である。 上式で( )の右下つきの添字は、その変数を一定に保ったままの微分を表わす。

1.1.5 熱エネルギーの保存則

マントル対流で考慮すべきエネルギーには、熱エネルギーと運動エネルギーの2 種類があるが、ここでは熱エネルギーの保存則のみを考える。 運動エネルギーの保存則は運動量の保存則から導かれるものなので、ここでは考 えない。

熱力学の第2法則より、熱エネルギーの保存則は

ρTDsDt=Φ-qixi+ρH (1.16)

と表わされる。 ここで、Tは温度、sは比エントロピー(単位質量あたりのエントロピー)、 Hは単位質量あたりの内部加熱率である。 右辺第1項は粘性散逸による発熱で、

Φ=τijvixj (1.17)

と書ける。 また右辺第2項は伝導による熱輸送であり、熱流量𝒒はFourierの法則から

qi=-kTxi (1.18)

と表わされる。 ここでkは熱伝導率である。

粘性散逸Φをあらわに書くと、

Φ =[μ(vixj+vjxi)+(kB-23μ)vkxkδij]vixj
=kB(vixi)2+2μ[12(vixj+vjxi)-13vkxkδij]2 (1.19)

ここで右辺第1項は体積変化に起因する項、第2項は剪断変形(体積変化を伴なわ ない)に起因する項を表わす。 さらにStokes近似(kB=0)かつ非圧縮性を仮定すれば、式(1.19)は

Φ=μ2(vixj+vjxi)2 (1.20)

と簡単化される。

式(1.16)の左辺を具体的にどう書き表わすかによって、熱エネルギーの 保存則はさまざまな形で書き表わすことができる。 導出の詳細は補遺Aを参照のこと。 例えば、 温度Tと圧力pを独立変数として選ぶと、

ds=CpTdT-αρdp (1.21)

となる。 ここで、Cpは定圧比熱である。 これを式(1.16)に代入すれば

ρCpDTDt-αTDpDt=xi(kTxi)+Φ+ρH (1.22)

を得る。

1.2 マントル対流向けの近似

以下では、マントル内の熱対流を考える場合に有効な近似を導入し、前章で導出 した基礎方程式を簡略化する。

まず、状態方程式(1.13)の簡単化を考える。 マントル物質の相転移や組成の違いによる密度変化を無視すれば、マントルの密 度ρは第一義的には温度Tと圧力pの関数で与えられる。

ρ=ρ(T,p)

密度変化は主に静水圧下の圧縮によるものが大きい(マントルの上面と下 面で約65%程度)ことが知られている。 これより、式を簡単にするため、ρ=ρ¯+ρT=T¯+Tp=p¯+p の如く、基本場 (¯つき)と基本場からのずれ(つき)の2つに分ける。

ρ =ρ¯+ρ=ρ(T¯+T,p¯+p)
ρ¯(T¯,p¯)+(ρp)T¯p+(ρT)p¯T
=ρ¯(T¯,p¯)+ρ¯χ¯Tp-α¯ρ¯T (1.23)

その際、定常かつ運動のない状態を基本場として採用するのが便利である。 さらに簡単のため、基本場の空間変化は重力の方向にのみ起こると仮定する。 この仮定から、圧力の基本場が満たすべき式は

0=-p¯xi+ρ¯g¯i (1.24)

で与えられる。

以降は、式(1.23)を基にして、基礎方程式を簡略化していく。 その際、方程式のスケーリングを行ない、各項の大きさを見積っていくのが便利 である。 表1及び2にスケーリングに用いる諸量をまとめておく。

Table 1: スケーリングに用いる基準量
記号 意味 値の目安
ρr 対流している状態での特徴的な密度 4×103 kg m-3
ΔTr 対流を駆動する特徴的な温度差 103 K
χTr 等温圧縮率χTの基準値 3×10-12 Pa-1
αr 熱膨張率αの基準値 3×10-5 K-1
μr 粘性率μの基準値 1021 Pa s
kr 熱伝導率kの基準値 4 W m-1K-1
Cpr 定圧比熱Cpの基準値 103 J kg-1K-1
b 対流層の厚さ 3×106 m
Table 2: 基礎方程式の無次元化に用いるスケール
スケール 値の目安
長さ b 3×106 m
時間 b2ρrCpr/kr 9×1018 s
速度 kr/bρrCpr 3.3×10-13 m s-1
圧力 μrkr/b2ρrCpr 1.1×102 Pa

以下では、無次元化された量をで表記する。

線型化された状態方程式(1.23)をスケーリングすると、

ρρ¯ =1+χ¯TpkrμrχTrρrCprb2-α¯TαrΔTr
=1+χ¯TpM2Pr-α¯Tϵ (1.25)

ここで新たに3つの無次元量を定義した。

M2kr2χTrρrCpr2b2,PrμrCprkr,ϵαrΔTr (1.26)

MはMach数(系の特徴的な速度と音速の比)、PrはPrantdl数(運動量の拡散率 と熱拡散率の比)であり、ϵは熱的な浮力の目安を表わす。 表1の値を用いると、これらの無次元量は M210-33Pr2.5×1023ϵ3×10-2 と推定できる。 これより、式(1.25)内でM2Prϵを含む項を無視してもよ いことが理解できる。 M2Pr0 かつ ϵ0 とおく近似を非弾性近似 (anelastic approximation) という。

次に質量保存則(1.1)をスケーリングする。 スケーリングされた状態方程式(1.25)を代入して整理すると、

ρ¯(χ¯M2Prpt-α¯ϵTt)+xi[ρ¯(1+χ¯TpM2Pr-α¯Tϵ)vi]=0 (1.27)

式(1.27)に非弾性近似を適用すると

xi(ρ¯vi)=0 (1.28)

を得る。 この式はρの時間変化を無視したことになっており、これによ り基礎方程式から弾性波を除去することができる。 実際、マントル内の弾性波(要するに地震波)は数km/秒で伝播する現象であり、 マントル対流が対象とする時間スケール(数百万年以上)では当然除去すべきもの である。 さらに、ρ¯の空間変化が無視できる(Boussinesq近似)と仮 定すれば、式(1.28)は

vixi=0 (1.29)

となり、非圧縮性流体の条件と一致する。

ここで、Boussinesq近似の意味を考えてみる。 断熱的な場合におけるρ¯の空間変化は

ρ¯xi=(ρp)s¯(pxi)=ρ¯χa¯(pxi)

ここでχaは断熱圧縮率である。 圧力の空間変化が静水圧的であると仮定すると、

ρ¯xi=ρ¯χa¯ρ¯g¯i=ρ¯2χa¯g¯i (1.30)

式(1.30)から、ρ¯の変化が起こる空間スケール hd¯を見積ると、

hd¯ =[1ρ¯|ρ¯xi|]-1 (1.31)
=1ρ¯χa¯g¯=γ¯Cp¯α¯g¯ (1.32)

ここで、Grüneisen parameter γ

γ=αρCvχT=αρCpχa (1.33)

を導入した。 式(1.32)で定義されたhd¯bよりも十分大きければ、 ρ¯の空間変化を無視してもよいであろう。 即ち、

1bhd¯=bα¯g¯γ¯Cp¯

右辺をスケーリングしてやると

1α¯g¯γ¯Cp¯αrgrbγrCpr=α¯g¯γ¯Cp¯Dγr (1.34)

ここで

DαrgrbCpr (1.35)

はDissipation numberと呼ばれる量である。 Dissipation number と Grüneisen parameter の比D/γrはマントルの 厚さとマントル内の密度変化の空間スケールとの比を表わす。 式(1.34)より、Boussinesq近似とは、D/γr0とし た極限に相当する。 しかし、実際のマントルではD0.5γ1と見積られる。 Boussinesq近似は式が簡便なためよく用いられるが、実際はその程度でしか正し くないことに注意しよう。

Navier-Stokes方程式(1.10)のスケーリングを行なう。 スケーリングされた状態方程式(1.25)を代入して整理すると、

1Prρ¯(1+χ¯TpM2Pr-α¯Tϵ)DviDt =-pxi+xj[μ(vixj+vjxi-23vkxkδij)]
+ρrχTrgrbg¯iρ¯χ¯Tp-b3ρr2CprgrαrΔTrμrkrρ¯α¯g¯iT

式(1.33)と(1.35)より、右辺第3項は

DγrCprCvrg¯iρ¯χ¯Tp

となり、右辺第4項から新たな無次元量Ra

Rab3ρr2CprgrαrΔTrμrkr (1.36)

と定義する。 RaはRayleigh数と呼ばれる。 これらを用いると、スケーリングされたNavier-Stokes方程式は結局

1Prρ¯(1+χ¯TpM2Pr-α¯Tϵ)DviDt =-pxi+g¯iρ¯χ¯TpDγrCprCvr
-ρ¯α¯g¯iTRa+xj[μ(vixj+vjxi-23vkxkδij)] (1.37)

となる。 ここで 1/Pr0 であることを考えると、式(1.37)の左辺 は無視することができる。 これより、非弾性近似を適用した場合には式(1.37)は

0= -pxi+g¯iρ¯χ¯TpDγrCprCvr-ρ¯α¯g¯iTRa
+xj[μ(vixj+vjxi-23vkxkδij)] (1.38)

となる。 式(1.38)の右辺より、対流を駆動する浮力には、(i)温度変化による密 度変化(第3項)、及び (ii) 圧力変化による密度変化(第2項)、の2つの寄与があ ることが分かる。 しかし一般的には、(ii) の寄与を無視して、

0= -pxi-ρ¯α¯g¯iTRa+xj[μ(vixj+vjxi-23vkxkδij)] (1.39)

とした方程式を用いることが多い。 これを「打ち切りされた非弾性近似 (truncated anelasic liquid approximation)」と呼ぶ。 この場合、浮力に寄与するのは温度変化による密度変化のみになる。 さらに、Boussinesq近似を施すと、

0=-pxi-g¯iTRa+xj[μ(vixj+vjxi)] (1.40)

と簡単化される。 ただしここでは Boussinesq 近似を利用して、 ρ¯=α¯=const.1 かつ vkxk=0 を用いた。

最後に、熱エネルギーの保存則(1.22)式のスケーリングを行なう。 スケーリングされた状態方程式(1.25)を代入して整理すると、

ρ¯(1+χ¯TpM2Pr-α¯Tϵ)Cp¯DDt(T¯+T)-αrμrkrb2ρr2Cpr2=ϵDRaα¯(T¯+T)DDt(p¯+p)
=xi[k¯xi(T¯+T)]+μrkrb2ρr2Cpr2ΔTr=DRaΦ+ρ¯(1+χ¯TpM2Pr-α¯Tϵ)H(b2HrρrkrΔTr)

ここで、内部発熱率の基準値Hrを導入した。 さらに左辺第2項に含まれる、圧力の基本場p¯のLagrange微分 の項を変形する。 基本場は時間変化しないと仮定しているので、

Dp¯Dt=p¯t+vip¯xi=vip¯xi

一方、式(1.24)をスケーリングすると

p¯xi=b3ρr2Cprgrμrkrρ¯g¯i=Raϵρ¯g¯i

であるから結局

Dp¯Dt=Raϵρ¯g¯ivi

これを代入して整理すると、スケーリングされた熱エネルギーの保存則は結局、

ρ¯(1+χ¯TpM2Pr-α¯Tϵ)Cp¯DDt(T¯+T)-α¯(T¯+T)DpDtϵDRa-α¯(T¯+T)g¯iviρ¯D
=xi[k¯xi(T¯+T)]+ΦDRa+ρ¯(1+χ¯TpM2Pr-α¯Tϵ)H(b2HrρrkrΔTr) (1.41)

これに非弾性近似を適用すると、

ρ¯Cp¯DDt(T¯+T)-α¯(T¯+T)g¯iviρ¯D=xi[k¯xi(T¯+T)]+ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.42)

を得る。

温度の基本場T¯をしかるべく仮定すると、式(1.42) はさらに簡単化できる。 基本場は時間変化しないと仮定しているので、

ρ¯Cp¯DT¯Dt=ρ¯Cp¯viT¯xi

加えて、基本場は断熱的であると仮定すると、断熱温度勾配の式(A.19) をスケーリングすることにより、

T¯xi=αrgrbCprα¯g¯iCp¯T¯=Dα¯g¯iCp¯T¯

となる。 これより、

ρ¯Cp¯DT¯Dt=ρ¯α¯T¯Dg¯ivi (1.43)

が得られる。 式(1.43)を式(1.42)に代入してやると、非弾性近似における 熱エネルギーの保存則

ρ¯Cp¯DTDt-α¯Tg¯iviρ¯D=xi[k¯xi(T¯+T)]+ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.44)

が得られる。 左辺第2項が断熱温度勾配に沿った温度変化の効果を表わしている。 加えて熱伝導率などが一様であると仮定すれば、右辺第1項内の熱伝導項を簡単 にすることもできる。 さらにBoussinesq近似を適用すると、D0の場合には断熱温度勾配 が0になるので、

ρ¯Cp¯DTDt=xi(k¯Txi)+ρ¯H(b2HrρrkrΔTr) (1.45)

となる。

以上で、マントル対流向けの基礎方程式系の近似が完了した。

これまでに適用した非弾性近似とBoussinesq近似の他に、マントル対流シミュレー ションには「拡張されたBoussinesq近似」(extended Boussinesq approximation) なる近似が使われることもある。 この近似は、Boussinesq近似で課された D/γr0の条件を D0 ではなく、γ かつ D0 と して適用するものである。 この場合、質量保存則と Navier-Stokes方程式はBoussinesq近似によるもの(式 (1.3)及び(1.40))と同じになり、熱エネルギー保存則は非弾 性近似によるもの(式(1.44))を用いることになる。

1.2.1 てっとり早く知りたい方のための基礎方程式のまとめ

  1. 1.

    基礎方程式の一般形 (次元つき)

    • 質量保存則

      1ρDρDt+vixi=0 (1.2)
    • 運動量保存則 (Navier-Stokes 方程式)

      ρDviDt=-pxi+xj[μ(vixj+vjxi-23vkxkδij)]+ρgi (1.10)
    • 熱エネルギー保存則

      ρCpDTDt-αTDpDt=xi(kTxi)+Φ+ρH (1.22)
  2. 2.

    非弾性近似 (ϵ1, M2Pr1) かつPrandtl数無限大の 場合 (次元なし)

    • 基本場は運動なし、定常、静水圧的、かつ断熱的とした

    • 質量保存則

      xi(ρ¯vi)=0 (1.28)
    • 運動量保存則 (Navier-Stokes 方程式)

      0= -pxi+g¯iρ¯χ¯TpDγrCprCvr-ρ¯α¯g¯iTRa
      +xj[μ(vixj+vjxi-23vkxkδij)] (1.38)
    • 熱エネルギー保存則

      ρ¯Cp¯DTDt-α¯Tg¯iviρ¯D =xi[k¯xi(T¯+T)]
      +ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.44)

    ちなみに、次元つきの量に戻して書き直すと以下のようになる。

    • 質量保存則

      xi(ρ¯vi)=0 (1.46)
    • 運動量保存則 (Navier-Stokes 方程式)

      0= -pxi+ρ¯g¯i(χ¯Tp-α¯T)
      +xj[μ(vixj+vjxi-23vkxkδij)] (1.47)
    • 熱エネルギー保存則

      ρ¯Cp¯DTDt-α¯Tg¯iviρ¯ =xi[k¯xi(T¯+T)]+Φ+ρ¯H (1.48)
  3. 3.

    打ち切りされた非弾性近似 (ϵ1, M2Pr1, 圧力変化 に起因する浮力を無視) かつ Prandtl数無限大の場合 (次元なし)

    • 基本場は運動なし、定常、静水圧的、かつ断熱的とした

    • 質量保存則

      xi(ρ¯vi)=0 (1.28)
    • 運動量保存則 (Navier-Stokes 方程式)

      0= -pxi-ρ¯α¯g¯iTRa+xj[μ(vixj+vjxi-23vkxkδij)] (1.38)
    • 熱エネルギー保存則

      ρ¯Cp¯DTDt-α¯Tg¯iviρ¯D =xi[k¯xi(T¯+T)]
      +ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.44)
  4. 4.

    Boussinesq 近似 (ϵ1, M2Pr1, D1) かつ Prandtl数無限大の場合 (次元なし)

    • 基本場は運動なし、定常、 ρ¯,C¯p, α¯が 一定

    • 質量保存則

      vixi=0 (1.29)
    • 運動量保存則 (Navier-Stokes 方程式)

      0=-pxi-g¯iTRa+xj[μ(vixj+vjxi)] (1.40)
    • 熱エネルギー保存則

      ρ¯Cp¯DTDt=xi(k¯Txi)+ρ¯H(b2HrρrkrΔTr) (1.45)

    ちなみに、次元つきの量に戻して書き直すと以下のようになる。

    • 質量保存則

      vixi=0 (1.29)
    • 運動量保存則 (Navier-Stokes 方程式)

      0=-pxi-ρ¯g¯iα¯T+xj[μ(vixj+vjxi)] (1.40)
    • 熱エネルギー保存則

      ρ¯Cp¯DTDt=xi(k¯Txi)+ρ¯H (1.45)

ただし、以下の無次元パラメータを含んでいる。

Pr =Prandtlnumber μrCprkr (1.26)
ϵ αrΔTr (1.26)
M2 kr2χTrρrCpr2b2 (1.26)
D =Dissipationnumber αrgrbCpr (1.35)
Ra =Rayleighnumber b3ρr2CprgrαrΔTrμrkr (1.36)
RaH =internalheatingRayleighnumber b5ρr3CprgrαrHrμrkr2 (1.49)

1.3 マントル対流の基本場の設定

ここでは、先の導出では十分に検討しなかった、基本場の設定について論じる。 以下のいずれの場合も、運動がなく静水圧的で、かつ定常な状態を基本場とする。

基本場とは、基礎方程式の扱いを容易にするための数学的な方便の一つというべ きものであるので、必ずしも現実に存在しなくてもよい。 しかしながら、特に圧縮性を考える際には注意が必要である。 なぜなら、状態方程式を指定するのみならず、必要な熱力学変数を正しく (熱力 学に則って) 扱う必要があるからである。

1.3.1 Boussinesq 近似の場合

Boussinesq 近似の場合には、基本場の設定は極めて簡単である。 その理由の1つは、基本場での密度分布 ρ¯ を一定とみなしてよ いからである。

まず基本場での圧力分布を考えてみよう。 そのためには静水圧を記述する式(1.24)を積分すればよい。 式(1.24)の積分には重力加速度が必要だが、デカルト座標系の場合には これを一定とみなしてもよいであろうから、結局

p¯=ρ¯g¯z

と書ける。 ここで重力はz軸の正の向きに働くものとし、かつp¯(z=0)=0とし た。 一方、3次元球殻モデルを考える場合には、半径方向に重力加速度が変化するこ とを考慮する必要があるであろう。 例えば外径a、内径cの3次元球殻領域をとり、マントルの密度をρm、 コアの密度をρcとすると

g¯(r)=GM¯(r)r2 =43πρcGr(for0rc) (1.50)
=43πG[rρm+c3r2(ρc-ρm)](forcra) (1.51)

となる。 またこの場合の静水圧を記述する式は

dp¯dr=-ρ¯g¯(r) (1.52)

となるから、(1.51)と(1.52)より

p¯(r) =43πρmGc3(ρc-ρm)(1r-1a)+23πGρm2(a2-r2)(forcra)
=23πGρc2(c2-r2)+23πGρm2(a2-c2)+43πρmGc3(ρc-ρm)(1c-1a)(for0rc) (1.53)

を得る。 ただしここでもp¯(r=a)=0とした。

次に温度の基本場について考える。 Boussinesq 近似では断熱温度勾配が0であるから、T¯ を一定にと るのが最も簡単である。 しかしこれとは異なるT¯を採用することも可能である。 例えば

0=xi(k¯T¯xi)+ρ¯H (1.54)

で与えられる1次元定常熱伝導状態の解をとることもできる。 デカルト座標系のもとで具体的に書くと

0=ddz(k¯dT¯dz)+ρ¯H (1.55)

あるいは3次元球殻モデルでは

0=1r2ddr(r2k¯dT¯dr)+ρ¯H (1.56)

である。 これに加えてk¯Hの関数形を与えればT¯が求まる。 例えばk¯Hの両方とも一定であるとすれば簡単に積分できて結 局、

T(z) =-ρ¯H2k¯z2+c1z+c2 (1.57)
T(r) =-ρ¯H6k¯z2+c1r+c2 (1.58)

を得る。 ここでc1及びc2は境界条件によって決まる積分定数である。

1.3.2 非Boussinesq近似の場合

次に圧縮性流体あるいは非弾性流体の場合を考える。 この場合は、密度ρや熱膨張率αなどの物性値の深さ(圧力)変化を自 然に取り込める利点はあるものの、そのために基本場の状態が簡単に求まらない という難点もある。 また、研究者によって基本場のとり方もまちまちであったりもする。

ここでは熱力学に極力則った形式で、基本場の状態を求めていくことにする。 この場合、基本場を規定する以下の方程式系

  1. 1.

    静水圧の関係式(1.24)

    dp¯dz=ρ¯g¯ (1.59)
    ordp¯dr=-ρ¯g¯ (1.52)
  2. 2.

    断熱的な圧力変化と温度変化の関係(A.16)

    dT¯dp¯=α¯T¯ρ¯Cp¯ (1.60)
  3. 3.

    断熱圧縮率の熱力学的な定義(A.12)

    1ρ¯dρ¯dp¯=χa¯ (1.61)
  4. 4.

    Grüneisen parameter (1.33)

    γ¯=α¯ρ¯Cp¯χa¯ (1.62)

を適当な仮定のもとで積分することにより、基本場を構成することになる。 ただし、8個の未知量 (p¯ρ¯g¯T¯α¯Cp¯χa¯γ¯) を決定するには、上の4式の他に4 つの条件を指定する必要がある。

デカルト座標系の場合によく用いられるもの

1つの簡単な例は、g¯Cp¯α¯γ¯ を一定と仮定してしまうことにより得られる。 まず式(1.62)、(1.59)、(1.61)から Adams-Willamson の式

1ρ¯dρ¯dz=ρ¯g¯χa¯=g¯α¯Cp¯γ¯ (1.63)

が導かれる。 仮定により式(1.63)の右辺は定数なので、解析的に積分してやると

ρ¯(z)=ρ¯(z=0)exp[g¯α¯Cp¯γ¯z]=ρ¯(z=0)exp[zhd¯] (1.64)

を得る。 ここでhd¯は式(1.31)で定義した、密度変化のスケール ハイトである。 またzを対流層の厚さbで無次元化してやると、

ρ¯(z)=ρ¯(z=0)exp[zbbhd¯]=ρ¯(z=0)exp[zbD¯γ¯] (1.65)

となる。 ここで

D¯=α¯g¯bCp¯ (1.66)

は基本場の dissipation number (式(1.35)を見よ) である。 式(1.64)で得られたρ¯を式(1.59)に代入し て積分すると、基本場の圧力分布が

p¯(z)=ρ¯(z=0)g¯hd¯[exp(zhd¯)-1] (1.67)

と求まる。 ただしp¯(z=0)=0とした。 基本場の温度分布についても同様である。 式(1.59)と(1.60)より

1T¯dT¯dz=α¯g¯Cp¯ (1.68)

が得られ、さらに積分すると

T¯(z)=T¯(z=0)exp(α¯g¯Cp¯z)=T¯(z=0)exp(γ¯zhd¯)=T¯(z=0)exp(zD¯b) (1.69)

となる。 圧力変化と密度変化のスケールハイトはhd¯で共通だが、温度変化 のスケールハイトはhd¯/γ¯=b/D¯であ る。 この基本場の欠点は、深部で対流が起こりにくくなることである。 なぜなら、式(1.68)、(1.69)より断熱温度勾配を求めてみる と

(dT¯dz)a=α¯g¯Cp¯T¯=α¯g¯T¯(z=0)Cp¯exp(γ¯zhd¯)=(dT¯dz)a(z=0)exp(γ¯zhd¯) (1.70)

即ち、断熱温度勾配が深さとともに指数関数的に増大してしまうことになるから である。

これに対し、Tackley [9] では、上とは異なった方法が用いられ ている。 そこでは g¯Cp¯ は一定とみなすものの、 α¯γ¯ は以下のような関数形で密度に依 存して変化すると仮定する。

ααo =exp{-δToκ[1-(ρoρ)κ]} (1.71)
γγo =ρoρ (1.72)

ここで、式(1.71)は経験的に得られている Anderson-Grüneisen parameter δT

δT(lnαlnρ)T=δTo(ρoρ)κ (1.73)

が断熱線に沿っても成り立つとの仮定による。 式(1.73)に現われる δT はおよそ46κはおよそ1.4程度の値を持つことが実験的に知られている。 また式(1.72)は Mie-Grüneisen の状態方程式からの派生であ る (おそらく)。 以上の条件式を数値的に積分することにより基本場が構成できる。

1.4 相転移の取り扱い

相転移を厳密に熱力学に則って計算するのは繁雑である。 ここではマントル対流業界でよく用いられている、拡張された Boussinesq 近似 のもとでの phase function に基づく扱いについて紹介する。 以下の記述は主に [2] を参考にした。 また以下では簡単のため相転移は1つのみと仮定しているが、相転移が複数ある 場合も同様な定式化により取り扱いが可能である。

仮定している相平衡図、及び相転移の計算に用いる “excess
pressure”
Figure 1: 仮定している相平衡図、及び相転移の計算に用いる “excess pressure” π の定義の概念図。
よく用いられる phase function
Figure 2: よく用いられる phase function Γ とその微分 dΓdπのグラフ。

方程式を簡単化するため、以下では Clapeyron 勾配 γ を一定にとり、 相転移が起こる圧力ptを温度Tの関数として

pt(T)=po+γT (1.74)

と仮定する。 ここでpoは温度T=0において相転移が起こる圧力である。 この値からの圧力のずれによって相転移の起こり方が決まると仮定しよう。 そこで以下のように余剰圧力(excess pressure)とでもいうべき量πを定義す る。

π=p-pt(T)=p-po-γT (1.75)

(実際には静水圧を仮定して圧力を深さに読み換える。) さらに、2つの相の分率を表わす phase function Γ (ただし 0Γ1Γ=0 で低圧相、Γ=1 で高圧相)を π の関 数として導入する。 よく使われる形は

Γ(π)=12[1+tanh(πΔp)] (1.76)

である。 ここで、Δp は相転移にかかる圧力の幅に相当するパラメータである。 また、式(1.76)をπで微分すると、

dΓdπ=12Δp[1-tanh2(πΔp)]=2ΔpΓ(1-Γ) (1.77)

となる。 これらのグラフを図2に示す。 この関数 Γ は相転移を数値的に簡便に取り扱うために恣意的に導入した に過ぎないことに注意しよう。 しかも後に見るように、Γ の関数形そのものではなく、その微分 dΓdπ のみが意味を持っているので、式 (1.76) を盲目的に信じる必要もない。 実際に Honda et al. [4] などでは、dΓdπ を Gaussian と仮定してシミュレーションを行っている。

phase function Γ を用いて、マントル対流の基礎方程式に相転移の効果 を導入することを考えよう。 以下では簡単のため、密度と比エントロピー以外の全ての熱力学的な量は相転移 によって不変である(1次の相転移に相当)と仮定する。 まず状態方程式を修正しよう。 マントル物質の密度は温度Tと圧力pに加え、phase function Γ にも 依存する関数で与えられるとしよう。

ρ=ρ(T,p,Γ)

1.2章の取り扱いに倣い、密度分布を基本場(¯つき)と基本 場からのずれ(つき)の2つに分離すると、

ρ =ρ¯+ρ=ρ(T¯+T,p¯+p,Γ¯+Γ)
ρ¯(T¯,p¯,Γ¯)+(ρp)T,Γ¯p+(ρT)p,Γ¯T+(ρΓ)p,T¯Γ
=ρ¯(T¯,p¯,Γ¯)+ρ¯χ¯Tp-α¯ρ¯T+Δρ¯Γ (1.78)

ここで、 Δρ¯ρh¯-ρ¯ は高 圧相と低圧相の密度差である。 (今後も同様に、下つき添字で低圧相、hで高圧相を示す。) 以下では簡便のため Δρ¯ρ¯=Δρρ と書き、これを一定とみなす。 式(1.25)に倣ってスケーリングを施すと、

ρρ¯=1+χ¯TpM2Pr-α¯Tϵ+ΓΔρρ (1.79)

を得る。 同様に、エントロピーの変化量(式(1.21))は以下のように修正されるべ きであろう。

ds=CpTdT-αρdp-ΔsdΓ (1.80)

ここで、Δss-sh は低圧相と高圧相の比エントロピーの差 である。 なお、Clausius-Clapeyron の式(A.24)より、ΔsΔρの間には

γ=s-shρ-1-ρh-1=ρρhs-shρh-ρ=ρρhΔsΔρ (1.81)

なる関係がある。 これを用いれば式(1.80)は、

ds=CpTdT-αρdp-γΔρρ2dΓ (1.82)

とも表わすことができる (ここでρρhρ2と近似した)。

これらを用いて、拡張されたBoussinesq近似のもとでの基礎方程式を修正してみ よう。 この場合、質量保存則はBoussinesq近似と同じ非圧縮の仮定(式(1.29)) を用いてよい。 以下では Navier-Stokes方程式と熱エネルギーの保存則がどのように修正される べきかを論じる。 相転移による密度変化の影響を考慮すると、スケーリングされた非弾性近似のも とでの Navier-Stokes 方程式は

0=-pxi+g¯iρ¯χ¯TpDγrCprCvr-ρ¯α¯g¯iTRa+ρ¯g¯iΓRb+xj[μ(vixj+vjxi-23vkxkδij)] (1.83)

と書ける。 ここで新たな無次元量Rb

Rbb3ρr2CprgrμrkrΔρρ (1.84)

と置いた。 これは相転移による密度変化が対流を促進(または抑制)する効果の強さを表わす 量で、phase boundary Rayleigh 数と呼ばれることもある。 さらに Boussinesq 近似を施すと、

0=-pxi-ρ¯α¯g¯iTRa+ρ¯g¯iΓRb+xj[μ(vixj+vjxi-23vkxkδij)] (1.85)

を得る。 これと式(1.37)と比較すると、温度変化による密度変化に加え、相転 移による密度変化が浮力に寄与していることが分かる。

次に、熱エネルギーの保存則を修正する。 式(1.82)を式(1.16)に代入すると、

ρCpDTDt-αTDpDt-TγΔρρDΓDt=xi(kTxi)+Φ+ρH (1.86)

となる。 左辺第3項が相転移によって新たにつけ加わった項であり、これは相転移による エントロピー変化の効果を表わしている。 次にこの式のスケーリングを試みる。 左辺第1項、第2項、及び右辺は式(1.41)などと同じなので省略し、式 (1.41)の導出にならって左辺第3項のスケーリングを考えると、

-TγΔρρDΓDt -ΔTr(T¯+T)ρrgrbΔTrγΔρρkrb2ρrCprDΓDt
=-krΔTrb2grbCprΔTrΔρρ=DRbRaγ(T¯+T)DΓDt (1.87)

ここで、Clapeyron 勾配 γ に含まれる圧力は静水圧 ρrgrb に よりスケーリングした。 これより、非弾性近似のもとでのスケーリングされた熱エネルギーの保存則は、

ρ¯Cp¯DDt(T¯+T)-Dα¯(T¯+T)ρ¯g¯ivi-DRbRaγ(T¯+T)DΓDt
=xi[kxi(T¯+T)]+ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.88)

相転移の効果は左辺第3項に含まれている。 この項は Dissipation number D でスケーリングされていることに注意しよう。 これより、D0の極限では熱輸送に相転移の効果は現われないこと が分かる 11 1 この議論から、Boussinesq 近似のもとでは相転移の潜熱の効果を無 視することが多い。ただし、Ogawa and Nakamura [5] によれば、固 相-液相の相転移では Boussinesq 近似であっても潜熱の効果は必ずしも無視す べきではないと指摘されている。

これまでに議論した基礎方程式の導出は、phase function Γ の具体的な 関数形にはよらないものであった。 また式(1.88)には Γ の Lagrange 微分が含まれているの で、解くのが少々面倒である。 そこで、式(1.75)及び(1.76)を用いてこ れらの式をより簡便な形に書き直してみる。 式(1.75)内の余剰圧力を静水圧ρrgrbでスケーリン グすると、

π=1ρrgrbμrkrb2ρrCpr=ϵRa(p-pt)-γT

ここで、ϵRa なる因子は圧力のスケーリングに用いた基準値 の違いに起因している。 これを用いると、phase function Γ の微小変化は

δΓ=dΓdπδπ=dΓdπ(ϵRaδp-γδT)

と与えられる。 これをNavier-Stokes方程式(1.83)に代入、整理すると、

0=-pxi+g¯iρ¯χ¯TpDγrCprCvr-ρ¯g¯iTRa(α¯+γRbRadΓdπ)
+xj[μ(vixj+vjxi-23vkxkδij)] (1.89)

を得る。 これと式(1.38)を比較すると、相転移による密度変化の影響は、実効 的な熱膨張率を

αeffα¯+γRbRadΓdπ=α¯+PdΓdπ (1.90)

の如く定義し直したことと同じである。 ここで新たな無次元量P

PγRbRa (1.91)

と置いた。 これは時に phase buoyancy parameter とも呼ばれる。 熱輸送方程式(1.88)の左辺第3項も同様に変形しよう。

DΓDt =dΓdπDπDt=dΓdπ[ϵRaDpDt-γDTDt]=dΓdπ[ϵRaDp¯Dt=ρ¯g¯ivi+ϵRaDpDt=0-γDDt(T¯+T)]
=dΓdπ[ρ¯g¯ivi-γDDt(T¯+T)] (1.92)

これを(1.88)に代入して整理すると、

[ρ¯Cp¯+DRbRaγ2(T¯+T)dΓdπ]DDt(T¯+T)-D(T¯+T)(α¯+γRbRadΓdπ)ρ¯g¯ivi
=xi[kxi(T¯+T)]+ΦDRa+ρ¯H(b2HrρrkrΔTr) (1.93)

を得る。 これと式(1.42)を比較すると、相転移によるエントロピーの出入りは、 式(1.90)で与えられる実効的な熱膨張率に加え、実効的な熱容量を

(ρCp)effρ¯Cp¯+DRbRaγ2(T¯+T)dΓdπ (1.94)

の如く定義し直したことと同じであることが分かる。

式(1.89)及び(1.93) をまとめると、相転移の効果は以下の3つに分かれていることが理解できる。

  1. 1.

    温度の違いによる相転移境界のずれ。 これは式(1.89)の右辺第3項の、相転移によ る浮力の変化として現われる。

  2. 2.

    潜熱の出入りによる温度変化。 これは式(1.93)の左辺第2項の、相転移に よる断熱温度勾配の変化として現われる。

  3. 3.

    潜熱の出入りによる温度変化が原因となる相転移境界のずれ。 これは式(1.93)の左辺第2項の、相転移に よる熱容量の変化として現われる。

ただし Boussinesq 近似を扱う場合には潜熱の出入りは無視されるので、上記の (2)及び(3)の効果は含まれない。 また、(3)の効果を残すと式(1.93)が温度に関し て非線型になってしまうため、この効果は多くの場合は無視されている。 とはいえ無次元化された Clapeyon 勾配の値は通常小さいので、 γ2を含む項を落とすことによる悪影響は小さいと期待される。

相転移が流れに及ぼす影響の概念図。Schubert et
al. 
Figure 3: 相転移が流れに及ぼす影響の概念図。Schubert et al. [6]より改。

上記の(1)と(2)の効果が対流に及ぼす影響を図3に概念的に示す。 一例として低温の下降流の挙動にのみ注目して考える。 まず(1)温度の違いによる相転移境界のずれの効果のみを考えよう。 γ>0の場合、低温の下降流は周囲より浅い部分で高密度の相へ転移する。 即ち、より大きな負の浮力を獲得することができる。 逆にγ<0の場合には下降流は周囲より深部で高密度の相へ転移するため、 負の浮力を獲得できない。 これより(1)の効果では、γ>0の相転移はこれを横切る流れを加速し、そ の反対にγ<0の相転移は阻害する方向に働く。 一方、(2)潜熱の出入りによる温度変化は(1)と正反対の効果を持つ。 低温の下降流がγ>0の相転移を通過する際に潜熱が解放される。 即ち下降流内の温度が上がり、熱的な負の浮力が小さくなる。 逆にγ<0の相転移を通過する際には潜熱が吸収され、下降流内の温度が下 がり、熱的な負の浮力が大きくなる。 これより(2)の効果では、γ>0の相転移はこれを横切る流れを阻害し、そ の反対にγ<0の相転移は加速する方向に働く。 ちなみに(3)の効果は、γの符号によらず、相転移を横切る流れを阻害す る方向に働く。

また、さらなる簡略化として、相転移面の起こる位置の温度変化を無視するとい う近似もある。 この場合、dΓdπ の計算において温度 T を一定とみな し、これを深さのみの関数と仮定する。 さらに、この dΓdπ をもとにして得られた αeff を用いて、相転移における密度変化の影響を取り入れる。 この近似では、温度分布が変化しても相転移面の位置を不変に保つことができる ことに加え、相転移による対流の層構造の有無に関してよい近似となっている [2] ので、計算の簡略化としてよく用いられている (例えば [4, 7]など)。