摩擦熱を考慮した剪断変形のシミュレーション
亀山 真典，堀 高峰，Phil R. Cummins，平野 聡，馬場 俊孝，金田義行(海洋科学技術センター 海底下深部構造フロンティア ∗ )Numerical experiments of shear deformation with frictional
heating
Masanori Kameyama, Takane Hori, Phil R. Cummins, Satoshi Hirano, Toshitaka
Baba and Yoshiyuki Kaneda(Frontier Research Program for Subduction Dynamics, Japan Marine Science andTechnology Center)∗〒 237-0061 神奈川県横須賀市夏島町 2-15

submitted to Bull. Earthq. Res. Inst., Univ. Tokyo
初稿: 平成 12 年 3 月 25 日
第 2 稿: 平成 12 年 5 月 29 日
受理: 平成 12 年 5 月 30 日

欄外表題: 摩擦熱を考慮した剪断変形のシミュレーション

final version of Kameyama et al.

1

Abstract
We develop a thermal-mechanical model for describing the formation of shear zones. We consider shear
deformation of a two-dimensional rectangular region composed of a viscous fluid under a constant velocity
at the boundary. Viscosity of the material is assumed to depend only on temperature for simplicity. In
order to enhance the shear deformation, we included a small inclusion whose viscosity differs from thatof the surrounding material. We carried out time-marching simulations and monitored the evolution oftemperature and strain around the inclusion. Our results show that the deformation localizes in a narrow
region when sufficient heat is generated by viscous dissipation. In the zone of localized deformation,
temperature increases by several hundred degrees owing to strong dissipative heating. We found that thezone of localized deformation develops in the region where the greatest heat is generated at the initialstage of evolution. The zone of localized deformation pierces the inclusion when the inclusion is weaker
than the surrounding material, while it develops away from the inclusion when the inclusion is stiffer
than the surrounding material. These findings, though qualitatively, suggest that heating by viscous
dissipation may play an important role in the formation of decollements around subducted seamounts.

Key words: frictional heating, shear deformation, decollements, earthquakes, two-dimensional finite difference methodfinal version of Kameyama et al.

1

2

はじめに
岩石の変形が狭い領域に集中して起こるという現象は，多くの地学的現象でみられる．断層運動や地震は，脆性領域で起こっている変形の集中の顕著な例であるが，変形の集中は脆性領域だけでなく塑性領域でも起こることが知られている [Sibson, 1977]．塑性領域での変形の集中は，プレート境界のモデル [例えばYuen et al., 1978]，プレート沈み込み開始の原因 [例えば Tackley, 1998]，あるいは深発地震の原因 [Ogawa,
1987; Hobbs and Ord, 1988] として重要な意味を持つと考えられている．ここでは，塑性変形領域での変形の集中がどのようにして起こるかを数値シミュレーションにより検討する．
変形の集中が起こるためには，変形に伴なう強度の低下が鍵となる．そうした強度の低下を塑性領域で起こす有力な機構の 1 つとして，熱による変形の集中 [例えば Poirier, 1980] が挙げられている．それによれば，変形に際して発生する摩擦熱によって，変形した岩石は温められ，温度が上昇するであろう．それと同時に，岩石の強度は温度に強く依存しているので，温度が上昇すると強度は低下するであろう．この両者の間のカップリングによって，発熱による温度上昇がきっかけとなって力学的な変形が促進され，変形が集中して起こる可能性がある．過去に行なわれた 1 次元数値シミュレーション [例えば Ogawa, 1987;

Kameyama et al., 1997, 1999] によれば，変形と温度上昇のカップリングによって変形が促進された結果，温度上昇と変形速度が集中する狭い領域 (シアゾーン) が形成しうることが示された．しかし，1 次元モデルを用いたこれらのシミュレーションの結果が，有限の大きさをもつプレートの変形に実際に適用できるかどうかは自明ではない．
この研究では，上のモデルを発展させたモデルで剪断変形シミュレーションを行ない，変形の集中がどのように起こるかを検討する．モデルの対象として，海洋プレートの厚さ程度の空間スケールをもつ系を，
プレート運動程度の速度で剪断変形させる状況を考える．ただし，摩擦熱を考慮した 3 次元粘弾性体の剪
断変形シミュレーションは数値計算上非常に困難であるので，ここではより現実的なシミュレーションを将来行なうための準備段階として，モデルを 2 次元粘性流体に限定して計算を行なう．

2

数値シミュレーションモデル
計算に用いたシミュレーションモデルの概略を Fig. 1 に示す．粘性流体で満たされた，厚さ 2d=100km，Fig. 1

幅 λd=200km の 2 次元長方形領域の剪断変形を考える．流体の Prandtl 数 (動粘性率と熱拡散率の比) は無限大と仮定し，系の時間変化は温度変化によって規定されるとする．剪断方向に x 軸をとり，それと直交する方向に z 軸をとる．領域の外側境界 z = ±d を一定の剪断速度 ±U で動かし，内部の流れと温度分布の変化を追跡する．変形を促進するため (詳細は第 3.1 章参照)，z = 0 沿いに小領域 δs (Fig.1 の陰ツキ領域) を置き，計算の初期段階では周囲と粘性率が異なるような領域を設けておく．なお簡単のため，領域の中心 z = 0 での対称性を仮定し，0 ≤ z ≤ d の領域のみを考えることにする．この章に現われるパラメータの値は Table 1 に示してある．

Table 1

final version of Kameyama et al.
2.13

基礎方程式系

流体の運動方程式は

−

· µ
¶¸
∂p
∂
∂vj
∂vi
+
η
+
=0
∂xi
∂xj
∂xi
∂xj
∂vk
=0
∂xk

(1)
(2)

と書ける．ここで，p は圧力，vi は i (i = x, z) 方向の速度，η は粘性率である．Prandtl 数無限大を仮定し∂viiているために，運動方程式から慣性項 (ρ[ ∂v∂t + vj ∂xj ]) が消え，結果的に力のつりあいの式と同等になっていることに注意されたい．流体の粘性率 η は温度 T に依存するとし，

η = ηo exp[−E(T − Tb )](3)のように置く．ここで ηo ，E ，Tb は定数である．定数は dry olivine [Karato and Wu, 1993] に近い粘性率を与えるように選んである．
温度の時間変化は以下の式で与えられる．
µ 2
¶
∂T
∂ T
∂2T
∂
∂
Φ
=κ
+
−
(vx T ) −(vz T ) +
∂t
∂x2
∂z 2
∂x
∂z
ρCp

(4)

ここで κ は熱拡散率，ρ は密度，Cp は比熱である．式 (4) のうち，右辺第 1 項は熱拡散による冷却，第 2
項，第 3 項は移流による熱輸送，第 4 項は粘性散逸 (摩擦) による発熱の寄与を表わす．摩擦による単位体
積あたりの発熱量 Φ はΦ =σII ε̇II= η(ε̇II )2(5)で与えられる．ここで，σII ，ε̇II はそれぞれ応力テンソル，歪速度テンソルの第 2 不変量で，各々の最大剪
断成分に相当する．
式 (4) の境界条件は，(1) 外側境界 z = d での温度を T = Tw に固定，(2) 中心 z = 0 では熱流束 0(z 方向の温度勾配が 0) とする．また，x = 0 と x = L の面で周期境界条件を課す．初期条件として，温度分布はT = Tw で一様とおく．ただし，小領域 δs に対応する (x − λd/2)2 + z 2 ≤ r2 の範囲では，T = Tw + δT のように温度を周囲より δT だけ変化させておく．

2.2

シミュレーションの方法

基礎方程式は差分法により離散化した．x 方向に 160 個，z 方向に 48 個の長方形のメッシュに分割して計算を行なった．メッシュは δs を含む 83 λd ≤ x ≤ 58 λd，0 ≤ z ≤ 21 d で細かく，それ以外では 2 倍の粗さにとった．計算に用いた最小のメッシュ幅は d/64 である．計算結果の妥当性を調べるためにメッシュ幅を
上記の半分にとって Case A1 (後述) の計算を行なったところ，変形の集中の起こり方には変化がないことを確認した．実際の計算は無次元量で行なっているが，計算の結果は次元量で表示した．無次元量から次元量への変換は，Table 2 に示した諸量で行なった．

Table 2

final version of Kameyama et al.

4

結果

3

本研究では剪断速度 U ，外側境界の温度 Tw ，小領域 δs に与えた温度コントラスト δT の 3 つの量を変化させて計算を行なった．以下で示す計算で用いた値は Table 3 にまとめてある．

3.1

Table 3

対照実験

最初に，この系の基本的性質を理解するため，小領域 δs に温度コントラスト δT を与えない場合 (Case

A0) と正の δT を与えた場合 (Case A1) での系の時間発展の違いをみてみる．(δT が負の場合については第3.3 章で論じる．) 計算に用いた U ，Tw の値はそれぞれ 12.6cm/yr，1000K で共通である．
Fig. 2 及び Fig. 3 に Case A0 の結果を示す．Fig. 2 は各時刻における，x(横) 方向に平均した温度 hT i のFig. 2z 方向のプロファイルを，Fig. 3 は領域内の (a) 最高温度 Tmax と平均温度 Tav ，(b) 最大歪速度 ε̇max の時

Fig. 36間変化を示す．剪断変形を開始させた時刻 t = 0 から 7.34 × 10 年までの間，摩擦熱によって領域全体で温度が上昇している (Fig. 3a)．z = 0 付近では摩擦により発生した熱が逃げにくいため，温度が大きく上昇しているが，外側境界 (z = d) 付近では境界からの冷却のため，温度の上昇が抑えられている．冷却の効果
は時間の経過とともに領域内部まで浸透していき，温度上昇の抑えられた領域は z = d 付近から徐々に広がっていく．z = 0 付近と z = d 付近の間の温度差は徐々に大きくなり，その結果領域内の粘性率コントラ
ストも大きくなって，最大歪速度 ε̇max も時間とともに増大していく．時刻 t = 7.34 × 106 年には，z = 0付近の温度は約 1360K まで上昇し，最大歪速度 ε̇max も平均歪速度 U/d の約 5 倍にまで大きくなる．しかし，この時刻以降は z = 0 での温度上昇は起こらなくなる．これは熱伝導による冷却の効果が z = 0 付近に到達し，摩擦熱による温度上昇と，熱伝導による冷却とがつりあってしまうからである．その後は系は徐々に定常状態に向かっていく．そして t = 5.34 × 107 年には，摩擦による加熱と，熱伝導による冷却が系内のあらゆる部分で等しくなった熱平衡状態に至る．
次に，δT =140K なる温度コントラストを与えた Case A1 の結果を示す．この δT により初期状態での

δs 内部と周囲との粘性率コントラストを e−4 (約 50 分の 1) にとっている．Fig. 4 は領域内の (a) 最高温度
Tmax と平均温度 Tav ，(b) 最大歪速度 ε̇max の時間変化を表わす．Fig. 4 より，Case A1 での系の進化は，
長い時間スケールの変化と短い時間スケールの変化が重なって起こっていることがわかる．初期の短い (約105 年) 期間では系内の最大歪速度 ε̇max と最高温度 Tmax ，平均温度 Tav が急速に増加するが，その後の長い (約 5 × 107 年) 期間では，ε̇max ，Tmax ，Tav は徐々に一定の値へと変化している．しかも，定常状態での Tmax ，Tav ，ε̇max の値は，Case A0 での定常状態での値と同じになっている．

Case A1 での系の時間発展を詳細に見るために，Fig. 5 に各時刻における (a) 温度 T ，歪速度 ε̇II の分布，
(b)x 方向に平均した x 方向の速度 hvx i の z 方向のプロファイルを示す．Fig. 5a では色で温度，コンターで歪速度の分布を表わしている．図より，初期の段階の時間発展とは z = 0 付近の薄い層内でのみ温度と歪速度が大きくなる段階で，その後の段階とは熱平衡状態に至るまでの段階であることがわかる．初期の段階
では，剪断変形を開始させた時刻 t = 0 より，摩擦熱によって領域全体で温度の上昇が起こる．面 z = 0 付

Fig. 4

final version of Kameyama et al.

5

近に注目すると，初期状態で温度コントラストを与えて粘性を低くしたことにより，δs に向かって流れが集中する (Fig. 10a も参照)．これによって δs の周囲での流れが加速され，面 z = 0 沿いに大量の摩擦熱が発生する．大量の摩擦熱の発生により，時刻 t = 2.20 × 104 年では，z = 0 沿いに温度の高い部分が δs を中心に横に広がるように発達する．この領域では歪速度が大きいため，周囲よりも大量の摩擦熱が発生し，温度も周囲よりも大きく上昇する．温度が上昇することにより粘性率が一層低下し，更に歪速度が大きくなる．このような，摩擦熱による温度上昇と，温度上昇による粘性率低下との間に正のフィードバックがはたらき，一旦大きく変形した領域がさらに変形しやすくなり，ε̇max ，Tmax が加速度的に増大する．ε̇max ，Tmax が極大となった 8.43 × 104 年では，z = 0 付近のごく薄い層の中に温度上昇と変形のほとんどが集中しており，温度は周囲より 200K 以上も高く，歪速度も周囲より 2 桁以上大きくなっている．また Fig. 5bから，変形の集中した層の中でのみ大きな速度勾配が存在し，それ以外の領域ではほとんど剛体的に動いている，すなわち「すべり面」に似た構造が発達していることがわかる．以下，このような高温かつ歪速度の大きい層ができることを，「変形の集中が起きた」ということにする．しかし，時刻 t = 106 年以降では

Tmax はほぼ一定の値となり，ε̇max は減少を始める．これは，z = 0 付近での粘性率の低下のために摩擦熱の発生量が小さくなり，摩擦による加熱が熱伝導による冷却にうち消されてしまうからである．z = 0 以外の領域では依然として摩擦熱によって温度が上昇し，z = 0 付近との温度差が小さくなる．そのため，系内の粘性率のコントラストも小さくなり，z = 0 付近への歪速度の集中も弱まっていく．そして t = 5.08 × 107
年後には，Case A0 で得られた熱平衡状態 (Figs. 2 及び 9 参照) と非常に近い状態になっている．
正の温度コントラスト δT を与える領域 δs の大きさを変えても，系の時間発展の特徴は大局的には変わらないことを確認した．δs の半径 r を 0.2d，0.02d に変えた計算 (Cases A1b, A1t，Table 3 参照) を行なったところ，いずれの場合も，熱平衡状態に至るまでに変形の集中を経由していった．ただし，変形の集中が起こるまでの時間は，r が大きいほど短くなった．これは δs が小さいほど z = 0 面の温度を上げるのに長い時間が必要となるからであると理解できる．同様に系の縦横比 λ を 2 にとった計算 (Case A1n) を行な
い，系の横幅を変えても変形の集中を経由していくことは変わらないことも確認した．
以上をまとめると，正の δT を与えることによって，系の時間変化は，熱平衡に至るまでの長い時間ス
ケールの変化に加えて，短い時間スケールでの変化が重なったものとなる．短い時間スケールの変化では，ある条件下 (第 3.2 章参照) では正の δT を与えた部分に沿って急激な温度上昇と歪速度の集中が起こる．また長い時間スケールをかけて達成される熱平衡状態は，境界条件 (Tw と U ) で決まり，δT の有無とは無関係であることもわかる．

3.2

変形の集中が起こる条件

正の δT を与えたときに見られた変形の集中がどのような条件下で起こるか調べるために，外側境界の温度 Tw と剪断速度 U をさまざまに変えて計算を行なった．この節では δT =140K を固定してある．
まず，Tw を一定に保った上で剪断速度 U のみを変化させて行なった計算の結果を Fig. 6 に示す．Fig.

6a(左) には系内の最大歪速度 ε̇max が極大となった時刻における，温度 T と歪速度 ε̇II の分布を表している．

Fig. 6

final version of Kameyama et al.

6

面 z = 0 沿いにある太線で描かれた半円は，初期条件で温度コントラスト δT を与えた領域 δs を表わす．

Fig. 6b(右) は，hvx i の z 方向のプロファイルを表わす．図より，U がある閾値 Uth (この場合は 6.31cm/yr)
よりも大きい場合 (Cases A2, A1, A3) には変形の集中が起きるが，U が Uth より小さい場合 (Case A4) では起こらないことがわかる．また，剪断速度が大きいほど，変形の集中した層が薄くなり，また温度の上昇
量も大きくなっている．
剪断速度 U が閾値 Uth より小さいときに変形の集中が起こらない原因をみるために，Fig. 7 に Case A4

Fig. 7での系の時間発展を示す．Case A4 でも計算開始直後から摩擦熱により Tav は徐々に上昇していくが，Tmax
は逆に徐々に低下している．これは，Case A4 では歪速度が小さく，摩擦熱の発生も少ないため，摩擦熱による温度上昇よりも熱伝導による冷却の効果が卓越するためである．そのため，領域 δs に与えた温度コン
トラストは熱伝導によりうち消されてしまう．温度コントラストが消えてしまったために，この場合の系の進化は δT を 0 とした場合と同様，変形の集中を経由せずに熱平衡状態に至ることになる．言い換えると，

U < Uth の場合は発生する摩擦熱の量が小さいために変形の集中が起こらないことがわかる．なお，境界条件が Case A0 とは異なる (Case A0 では U =12.6cm/yr) ために，最終的な熱平衡状態は Case A0 のそれとは異なったものとなっている．
変形の集中に必要な剪断速度の閾値 Uth が何によって決まるかを調べるために，Tw を 900K から 1100Kの範囲で変化させて同様の計算を行なった．その計算結果を Fig. 8 にまとめた．図の▲は変形の集中が起
こった Tw ，U の値を示し，●は起こらなかった値を示す．この図から，Uth は温度 Tw が高くなるほど大きくなることがわかる．このことは発生する摩擦熱の量の違いによって説明できる．Tw が高く粘性率が低いと，発生する摩擦熱は少なくなる (式 (5) 参照) ので，温度が高いほど変形の集中は起こりにくくなり，言い換えれば，変形の集中に必要な U も大きくなる．また，Fig. 8 の実線は，時刻 t = 0 で発生する単位体積あたりの摩擦熱 Φo が

Φo = σII (t = 0)ε̇av = η(T = Tw )

µ

U
d

¶2

= 4 × 10−4 W/m3

(6)となる Tw ，U の関係を表わす．この実線が変形の集中する場合としない場合を分けていることから，変形の集中が起きるためにはおよそ 10−4 W/m3 オーダーの摩擦熱が発生することが必要であり，この摩擦熱の大きさ Φo は Tw にはほとんどよらないこともわかる．
系の厚さ d を 50km から 1km に変えて同様の計算を行なったところ，変形の集中に必要な摩擦熱の大きさは d の 2 乗に反比例することを確認した．これは線型安定性解析 [Gruntfest, 1963; Ogawa, 1987] 及び

1 次元モデルでのシミュレーション [Yuen et al., 1978; Kameyama et al., 1999] から予想される結果と調和的である．この結果は摩擦熱による温度上昇と熱伝導による冷却の競合から理解できる．線型安定性解
析によれば，熱伝導による冷却の時間スケールと摩擦熱による温度上昇の時間スケールの比で定義される

Gruntfest 数 G がある値 Gth よりも大きくなれば変形の集中が起こる．熱伝導による冷却の時定数は d2 /κ
と表わされ [例えば Turcotte and Schubert, 1982]，また摩擦熱による温度上昇の時間スケールは Φ−1o に比例することから，G > Gth となるためには Φo ∝ d−2 でなければならないことが示される．

Fig. 8

final version of Kameyama et al.

7

ここまでの結果をまとめると，大量の摩擦熱が発生する場合，摩擦熱による温度上昇と，温度上昇による粘性率低下の間の正のフィードバックがはたらき，温度上昇と歪速度が狭い領域に集中することが示さ
れた．

3.3

温度コントラストを変化させた場合前小節までに示した計算では，小領域 δs に与えた周囲との温度コントラスト δT を非負の一定の値に固定してきた．ここでは，δT を変化させたときに系の時間変化がどう変わるかを調べる．

Fig. 9 に Tw ，U を一定にした上で，温度コントラスト δT の大きさを変化させて行なった計算の結果をFig. 9示す．図には，系内の最大歪速度 ε̇max が極大となった時刻における系の状態を示す．初期状態での δs とそれ以外の領域との粘性率コントラストは Case A1 で e−4 (約 50 分の 1)，Case D1 で e−1 (約 3 分の 1)，Case

A0 で 0，Case D2 で e，Case D5 で e4 となっている．すなわち，上 2 つは面 z = 0 付近に周囲よりも粘性の低い領域を挿入した場合で，下 2 つが周囲よりも粘性の高い領域を挿入した場合を考えていることに対応する．図より，δT の符号によって，温度上昇と歪速度の増大の起こる場所が変わること，すなわち δT が正の場合 (Cases A1，D1) では柔らかい領域 δs を貫くように z = 0 沿いに起こり，δT が負の場合 (Cases D2，

D5) では固い領域 δs を避けるように面 z = 0 から離れたところで起こること，また δT の絶対値によって，温度上昇と歪速度の増大の起きる場所は変わらないこと，がわかる．

Fig. 10 に Case A1 と D5 における，計算初期の段階での流れ場と，摩擦による発熱量の分布を示す．図から，δT の違いによって変形の集中する場所が異なるのは，最大の摩擦熱が発生する場所の違いによることがわかる．正の δT を用いた Case A1 では，粘性率の低い領域 δs に向かって流れが曲げられている．これにより，流れの集中する面 z = 0 沿いに大量の摩擦熱が発生し，そこで温度が大きく上昇することにな
る．その結果 Fig. 9 のように，z = 0 沿いに変形が集中した層が発達するようになる．一方，負の δT を用いた Case D5 では，領域 δs の粘性率が高いため，面 z = 0 付近での横方向の流れが妨げられ，流れは δsを迂回するように生じる．このため，摩擦による発熱は，流れが最も集中する，δs の直上で最大になり，そこでの温度上昇が大きくなる．温度の上昇によって，先に示した強度低下との正のフィードバックにより，時間の経過とともに，突起物の直上で温度と歪速度の上昇とが起こり，1.47 × 106 年後 (Fig. 9 参照) には
温度と歪速度の大きい領域で幅の広い帯状に発達するようになる．(ただし，この場合でも最終的な熱平衡状態は Case A0 のそれと同じになることに注意されたい．)

δs の半径 r を 0.2d，0.02d に変えて同様の計算を行なったところ，変形の集中の起こる場所と δs との関係は変わらないことを確認した．δT が正の場合 (Cases A1b, A1t)，変形の集中は δs を貫くように起こり，負の場合 (Cases D5b, D5t) は δs を迂回するように起こる．また δs の場所を x = λd/2, z = d/2 に変えて計算を行なった場合 (Cases A1s, D5s) でもその結果が変わらないことも確認した．

Fig. 9 はまた，同じ Tw ，U を与えたとき，δT が負の場合は，高温かつ歪速度の大きい領域が δT が正の場合のそれよりも非常に厚くなっていることも示している．この理由は，負の δT を与えた場合は正のδT を与えた場合と比べて変形の集中が起こりにくくなっているからである．負の δT を与えた場合，領域

Fig. 10

final version of Kameyama et al.

8δs の直上に最大の摩擦熱が生じる領域ができた (Fig. 10b)．この領域の z 座標を zf とすると，この状況はz = zf に正の温度擾乱を置き，みかけの領域の厚さが d − zf と小さくなったことに対応する．領域の厚さが小さくなれば変形の集中に必要な摩擦熱は大きくなる (第 3.2 章の議論を参照) ので，同じ Tw ，U の値で比較すれば δT が負の場合は正の場合と比べて変形の集中は起こりにくくなる．そのため，Fig. 9 に示した
条件下では，負の δT を与えた場合，高温の領域が薄く発達するかわりに，熱平衡状態と同様，幅の広い領域での温度の上昇が起こる．ただしこの場合も，発生する摩擦熱が十分大きくなれば，変形の集中が起こる．実際，剪断速度 U を 10 倍に大きくして同様の計算を行なうと，負の δT を与えた場合でも高温の領域の幅は非常に薄くなることを確認した．
この小節の結果をまとめると，変形の集中は，は初期の段階で最大の摩擦熱が発生するところに起こることがわかる．粘性率が周囲よりも低い領域があれば，そこに流れが集中して最大の摩擦熱が発生し，その結果変形の集中した層はその領域を貫くように発達する．逆に粘性率が周囲よりも高い領域があれば，そこを迂回するように流れが集中し，その結果変形の集中した層はその領域から離れた外側のところで発達する．ここでは δs に与える粘性率コントラストの原因として温度のみを考えていたが，温度以外の要因で δsに粘性率コントラストがつけられた場合でも，結果は同様である．δs に温度コントラストを与えずに，δsの領域のみ粘性率を強制的に変えた計算を行ない，変形の集中の起こる場所と δs との関係は変わらないことを確認した．

4

議論 · 考察
本研究での数値シミュレーションの結果，(1) 大量の摩擦熱が発生する場合には，摩擦熱による温度上昇と，温度上昇による強度低下との間の正のフィードバックにより，層状の薄い領域の中に温度上昇と変形が集中すること，また，(2) 初期に与える温度コントラストの符号によって，変形の集中の起こり方が異なる
こと，すなわち，周囲よりやわらかい物質があるときは，そこを貫くように変形の集中が起こり，周囲よりも固い物質があるときは，そこから離れたところで変形の集中が起こること，がわかった．変形の集中が起こるために必要な摩擦熱の量は，dry olivine の粘性パラメータを用いた場合はおよそ 10−4 W/m3 程度と見積もられた．この値は同様の粘性パラメータを用いて行なった 1 次元シミュレーション [Kameyama et al.,1999] から見積もられた値 (O(10−5 )W/m3 ) よりも大きい．この原因として，(a)1 次元モデルと 2 次元モデルの違い，(b) 粘弾性レオロジーと粘性レオロジーの違い，が考えられる．1 次元モデルでは予め面的に広がった温度コントラストが与えられているのに対し，ここで用いた 2 次元モデルではパッチ状の温度コン
トラストを与えている．パッチ状の温度コントラストから面的な高温の領域を発生させるために，より大量の摩擦熱が必要になったと考えられる．また，粘弾性レオロジーを用いた場合は，予め蓄えられた弾性エネ
ルギーが熱エネルギーに変換され，摩擦以外の熱源も系の温度上昇に寄与する [Ogawa, 1987; Kameyama

et al., 1999]．しかし純粘性レオロジーを用いた場合は，弾性エネルギーの寄与がないため，それに相当す
る量のエネルギーを摩擦によって供給する必要があるからと考えられる．
粘性流体によるシアゾーン形成モデルは，塑性変形が卓越する地球深部での変形の集中 (例えば深発地震

final version of Kameyama et al.

9の原因 [Ogawa, 1987]) に適用されてきた．ここでは，本研究の結果を地球表層での変形に適用することを考える．本研究で用いた数値シミュレーションモデルでは変形メカニズムとして粘性流動しか考慮していな
いため，弾性変形が卓越すると思われる低温の地球表層での現象と定量的に比較することはできない．しかし，摩擦熱の発生を考える上で不可欠なのは粘性変形 (不可逆変化) であるから，弾性変形 (可逆変化) を無視したことによる本質的な影響は少なく，定性的には比較が可能になると考えられる．そこで，付加体の内部で生じる変形への適用を考える．付加体中の堆積物の固結の度合が小さければ，付加体の変形は主に流動によって起こり，弾性的な変形を無視したことによる影響は小さいと期待できるからである．
温度コントラスト δT を負にとって，周囲よりも固い領域を挿入した計算 (Cases D4，D5) は，柔らかい付加体の中を海洋性地殻からなる固い海山が沈み込んでいる状況と比較できる．この計算の結果から，表面に海山を持つ海洋性地殻が沈み込むことによって付加体が剪断を受けると，剪断方向に延びた帯状の高温の領域が海山の上部に発達することが期待される．こうして発達した高温の帯は，海山の周辺に発達するデコルマ面に相当すると考えられる．近年の反射地震探査 [Park et al., 1999] により，室戸沖南海トラフの付加体内を沈み込む海山と，海山の上部のデコルマ面の存在が明らかになっている．また高温の帯の中では，含水鉱物の脱水が盛んに起こると期待される．含水鉱物の脱水により生じた水はデコルマ面の強度を低下させ，大きなすべり量を生み出すことが可能になると考えられる [Hubbert and Rubey, 1959]．加えて，デコルマ面に大量の水があれば，それによって高電気伝導度異常が観測される可能性もある．ただし，堆積物の粘性率およびその温度依存性の推定が困難であるので，摩擦熱による変形の集中が実際に付加体の中で起こるかどうかは明らかではない．また，実際の付加体と海洋プレートがプレート境界面に対して対称であるとは考えられないので，対称性を仮定したモデルの結果がどの程度適用できるかも明らかではない．ここでは摩擦熱による変形の集中とデコルマ面形成の因果関係を定性的に議論したが，今後より厳密な定量的議論をするためには，付加体の変形に特化したモデルで検証する必要がある．
この研究では，摩擦熱と変形の相互作用モデルにより，断層でのすべりに似た速度分布 (Fig. 6b) を再現することができた．一方で，地震がきっかけとなって温度変化が引き起こされる可能性が種々の観測事実や理論的考察から示唆されている．McKenzie and Brune [1972] によれば，地震時の断層のすべりによって，局所的，瞬間的ではあれ鉱物の融点を超えるほどの高温を発生させるのに十分な量の熱が摩擦によって発生しうる．また，地震断層の岩石サンプルの分析 [例えば大槻, 1998] や深発地震の震源パラメータ [Kanamori

et al., 1998] といった観測からも，地震のすべりによって鉱物の融解が起きた可能性が示唆されている．では逆に，地震のような断層のすべりが，ここで示した摩擦熱と変形の相互作用モデルにより励起されうるだ
ろうか? ここで考えている数値モデルでは，変形は粘性流動でのみ起こると仮定している．粘性流動のみを考えたモデルで得られたすべりは，地震で発生するすべりの特徴を満たしていない．粘性流体では応力が下がると歪速度も低下してしまい，応力の低下がすべりの促進に結びつかない．また，外側境界に与えた剪断速度より大きいすべり速度は得られない．従って，このモデルを地震発生メカニズムのモデルに発展させるためには，弾性変形の効果をとり入れることが不可欠である．モデル物質を粘弾性体に変えることにfinal version of Kameyama et al.

10より，主に弾性変形していると考えられているリソスフェアでの変形を扱うことが可能になる．また，弾性変形の効果を含めることにより，プレート速度より大きいすべり速度が弾性反発によって発生すると期待さ
れ [Ogawa, 1987; Kameyama et al., 1999]，高速なすべりによってより大量の摩擦熱が発生すれば，鉱物の融点を超えるところまで温度が上昇する可能性も考えられる [McKenzie and Brune, 1972]．今後は地震発
生メカニズムへの適用も視野に入れ，粘弾性流体を用いて同様のシミュレーションを行なう予定である．

謝辞
この計算のアイデアは，カリフォルニア工科大学の金森博雄先生との議論の中で生まれた．また，計算結
果とデコルマ面が対応する可能性は，海洋科学技術センターの鶴哲郎博士，朴進午博士に指摘していただ
いた．第 3.3 章の結果の解釈には，ミネソタ大学の Arkady Ten 博士との議論が非常に役にたった．本稿の査読にあたって下さった山野誠助教授，および匿名の査読者からは，本稿を改善するために多くの有益なコメントを頂いた．記して感謝の意を表わします．

final version of Kameyama et al.

11

参考文献
大槻憲四郎, 1998, 断層岩が示す地震性摩擦すべりの物理過程, 月刊地球, 号外 20 号, 213-218.

Gruntfest, I. J., 1963, Thermal feedback in liquid flow: Plane shear at constant stress, Trans. Soc.
Rheol., 7, 195-207.
Hobbs, B. E., and A. Ord, 1988, Plastic instabilities: Implications for the origin of intermediate anddeep focus earthquakes, J. Geophys. Res., 93, 10521-10540.
Hubbert, M. K., and W. W. Rubey, 1959, Role of fluid pressure in mechanics of overthrust faulting I:
mechanics of fluid-filled porous solids and its application to overthrust faulting, Geol. Soc. Amer.
Bull., 70, 115-166.
Kameyama, M., D. A. Yuen, and H. Fujimoto, 1997, The interaction of viscous heating with grain-size
dependent rheology in the formation of localized slip zones, Geophys. Res. Lett., 24, 2523-2526.
Kameyama, M., D. A. Yuen, and S.-I. Karato, 1999, Thermal-mechanical effects of low-temperatureplasticity (the Peierls mechanism) on the deformation of a viscoelastic shear zone, Earth Planet.
Sci. Lett., 168, 159-172.
Kanamori, H., D. L. Anderson, and T. H. Heaton, 1998, Frictional melting during the rupture of the1994 Bolivian earthquake, Science, 279, 839-842.
Karato, S.-I., and P. Wu, 1993, Rheology of the upper mantle: a synthesis, Science, 260, 771-778.
McKenzie, D. P., and J. N. Brune, 1972, Melting on fault planes during large earthquakes, Geophys. J.
Roy. Astr. Soc., 29, 65-78.
Ogawa, M., 1987, Shear instability in a viscoelastic material as the cause of deep focus earthquakes, J.
Geophys. Res., 92, 13801-13810.
Park, J.-O., T. Tsuru, Y. Kaneda, Y. Kono, S. Kodaira, N. Takahashi, and H. Kinoshita, 1999, A
subducting seamount beneath the Nankai accretionary prism off Shikoku, southwestern Japan,
Geophys. Res. Lett., 26, 931-934.
Poirier, J.-P., 1980, Shear localization and shear instability in materials in the ductile field, J. Struct.
Geol., 2, 135-142.
Sibson, R. H., 1977, Fault rocks and fault mechanisms, J. Geol. Soc. London, 133, 191-213.
Tackley, P. J., 1998, Self-consistent generation of tectonic plates in three-dimensional mantle convection,
Earth Planet. Sci. Lett., 157, 9-22.

final version of Kameyama et al.

12

Turcotte, D. L., and G. Schubert, 1982, Geodynamics: Applications of continuum physics to geologicalproblems, 450 pp., John Wiley, New York.
Yuen, D. A., L. Fleitout, G. Schubert, and C. Froidevaux, 1978, Shear deformation zones along majortransform faults and subducting slabs, Geophys. J. Roy. Astr. Soc., 54, 93-119.

final version of Kameyama et al.

13

Table 1: The adopted values of physical parameters.
symboldescription

value

Tw
δT

temperature at outer boundary
initial temperature-perturbation

900∼1100 K-140∼140 KrηoE

radius of region where initial temperature is Tw + δT
reference viscosity
temperature-dependence of viscosityTbρ

reference temperaturedensity

0.05d1022 Pa s
2.843 × 10−2 K−1

Cpκ

specific heat
thermal diffusivity

1127 K
3000 kg/m3
800 J/kg K
10−6 m2 /s

Table 2: The values of scaling parameters.
quantity

unit

lengthdtime
velocity

2

d /κ
κ/d

temperaturestrain rate
viscosity

1/E
κ/d2
ηo

stressηo κ/d2

values
5 × 104m
77.93 × 10
6.31 × 10−2

yr
cm/yr

35.17
4 × 10−16
1022

K
s−1
Pa·s

4 × 106

Pa

final version of Kameyama et al.

14

Table 3: The adopted values of parameters U , Tw and δT employed in numerical experiments presented
in this paper.
case

U [cm/yr]

Tw [K]

δT [K]

A0

12.6

1000

0

A2
A1

31.5
12.6

1000
1000

140
140

A3
A4
A5

6.31
3.15
1.26

1000
1000
1000

140
140
140

B1
B2

3.15
1.26

900
900

140
140

B3
B4
C1

0.631
0.315
31.5

900
900
1100

140
140
140

C2
D1

12.6
12.6

1100
1000

140
35

D2
D3
D4

12.6
12.6
12.6

1000
1000
1000

-35
-70
-105

D5
E1

12.6
12.6

1000
1100

-140
-140

A1b
A1t
A1n

12.6
12.6
12.6

1000
1000
1000

140
140
140

A1s
D5b
D5t

12.6
12.6
12.6

1000
1000
1000

140
-140
-140

D5s

12.6

1000

-140

final version of Kameyama et al.

15zd =50kmU (fixed)

T=Tw (fixed)
periodic

periodicη=η (T)δs0-dλd =200km

-U

T=Tw

Figure 1: Schematic of numerical model in this study.

xfinal version of Kameyama et al.

d

16

5.59x104
1.61x106
7.34x106
5.34x107

z

0
1000

temperature [K]1400

Figure 2: Plots of distributions of hT i, temperature averaged in x-direction, against z for Case A0 where

no temperature-perturbation is given. The elapsed time is shown by the numbers in the figure. Parametervalues are Tw =1000K and U =12.61cm/yr.

final version of Kameyama et al.

17

1400

[K]

T max
T av
(a) temperature1000
10 -11[s -1](b) strain rate

10 -14

.
ε max

U/d
0

time [x10 4 yrs]

5000

Figure 3: Temporal evolution of (a) maximum temperature Tmax (solid line) and spatially averaged
temperature Tav (dashed line), and (b) maximum strain rate ε̇max (solid line), for Case A0. In (b), theaverage shear strain rate ε̇av ≡ U/d is shown by a dashed line for comparison.

final version of Kameyama et al.

18

1800
1400

[K]

1000

0

T av

10
T max

(a) temperature

1000
10 -8
10 -11[s-1 ]

10 -14

0

10
.
ε max

U/d
(b) strain rate

10 -140time [x10 4 yrs]

5000

Figure 4: The same as Fig. 3 but for Case A1 where temperature-perturbation δT =140K is given.

final version of Kameyama et al.

19(b) hvx i/U(a) T and ε̇IIt = 0 yr

d
z

0

[ND]

1

t = 2.20 × 104 yr

0

d
z

-13

0

[ND]

1

t = 8.43 × 104 yr

d

-14.5

z

4

-1

-14
-13.5
-13
-12.5
-12

0

[ND]

1

t = 5.08 × 107 yr

0

d

-16.5
-16
-15.5
-15
-14.5
-14
-13.5
-13
-12.5

z

0
1000

0

1100

1200

temperature [K]

1300

1400

[ND]

1

0

horizontal
velocity [ND]

Figure 5: (a) Snapshots of distributions of temperature T (indicated by shadings) and strain rate ε̇II(indicated by solid lines), and (b) the profile of hvx i, velocity in x-direction averaged over 0 ≤ x ≤ λd,against z, for Case A1. The elapsed time is shown by numbers attached on each row. In (a), contours
indicate the logarithm of ε̇II in the unit of s−1 . The meaning of the gray-scale is indicated by the scale-bar
at the bottom of the figure. In (b), hvx i is normalized with U , velocity at the outer boundary z = d.
Parameter values are Tw =1000K, U =12.61cm/yr and δT =140K.

final version of Kameyama et al.

20(b) hvx i/U(a) T and ε̇II
4Case A2 (t = 5.25 × 10 yr)

-14.5

d

U=31.54cm/yr

z
-14
-13.5
-13
-12.5
-12
-11.5

-14.5

0

[ND]

1

Case A1 (t = 8.43 × 104 yr)
-14.5

0

d

U=12.61cm/yr

z

4

-1

-14
-13.5
-13
-12.5
-12

0

[ND]

1

Case A3 (t = 3.23 × 105 yr)
-15
-14.5

0

d

U=6.31cm/yr
-1

4

z

-14
-13.5
-13
-12.5

0

[ND]

1

Case A4 (t = 3.54 × 106 yr)
-16
-15.5
-15

d

U=3.15cm/yr

-14.5

z

-14
-13.5

-13

1000

0

0
1100

1200
1300
temperature [K]

1400

[ND]

1

0

horizontal
velocity [ND]

Figure 6: (a) Snapshots of distributions of temperature T (indicated by shadings) and strain rate ε̇II(indicated by solid lines), and (b) the profile of hvx i against z, for the cases indicated by the run numbersin the figure. The elapsed time is shown by the numbers in the figure.In (a), the meanings of gray-scaleand contours are the same as in Fig. 5a. The region outlined by solid semicircle indicates δs, where
temperature-perturbation is given by the initial condition. In (b), hvx i is normalized with U , velocity at
the outer boundary z = d. The values of U are shown by the numbers in the figure. Parameter values
are Tw =1000K and δT =140K.

final version of Kameyama et al.

21

1800
1400

[K]

1000

0

10
T max
T av
(a) temperature

1000
10 -9
10 -12[s-1 ]

10 -150

U/d
10 -15
0

.
ε max

10

(b) strain rate

time [x10 4 yrs]

Figure 7: The same as Fig. 4 but for Case A4 where U =3.15cm/yr.

5000

shearing velocity U [cm/yr]

final version of Kameyama et al.

22

102

101

100

10-1

d=50kmlocalized
not localized
900
1000
1100
temperature Tw [K]Figure 8: Diagram for the occurrence of the localization of deformation in terms of the velocity at theouter boundary U versus the temperature at the outer boundary Tw . Solid circles indicate that the
localization does not occur, and the solid triangles indicate that the localization occurs. The solid line
indicates the values of U and Tw which give the initial frictional heating rate Φo = 4 × 10−4 W/m3 .

final version of Kameyama et al.

23(b) hvx i/U(a) T and ε̇II
Case A1 (t = 8.43 × 10 yr)
4

d

-14.5

z

4

-1

-14
-13.5
-13
-12.5
-12

0

[ND]

1

Case D1 (t = 3.44 × 105 yr)
-15.5
-15
-14.5

0

d

U=12.61cm/yr

z
-14
-13.5
-13
-12.5
-12

0

[ND]

1

Case A0 (t = 5.34 × 107 yr)
-16.5
-16
-15.5
-15
-14.5
-14
-13.5
-13
-12.5

0

d

U=12.61cm/yr

z

0

[ND]

1

Case D2 (t = 1.49 × 106 yr)
-15.5
-15
-14.5
-14
-13.5
-13

0

d

U=12.61cm/yr

z

-13
-13.5

0

[ND]

1

Case D5 (t = 1.47 × 106 yr)
-15.5
-15
-14.5
-14
-13.5
-13

d

U=12.61cm/yr

z

-13
-13.5

0
1000

0

1100

1200

temperature [K]

1300

1400

[ND]

1

0

horizontal
velocity [ND]

Figure 9: (a) Snapshots of distributions of temperature T (indicated by shadings) and strain rate ε̇II(indicated by solid lines), and (b) the profile of hvx i against z, for the cases indicated by the run numbersin the figure. In (a), the meanings of gray-scale and contours are the same as in Fig. 5a. The region
outlined by solid semicircle indicates δs, where temperature-perturbation is given by the initial condition.
In (b), hvx i is normalized with U , velocity at the outer boundary z = d. Parameter values are Tw =1000K

and U =12.6cm/yr.

final version of Kameyama et al.

24(a) Case A1 (t = 5.33 × 102 yr)(b) Case D5 (t = 4.85 × 102 yr)

-4

-3.5

-3

3

-2.5

-2

log Φ[W/m ]Figure 10: Snapshots of distributions of frictional heating Φ (indicated by shadings) and streamlines(indicated by solid lines) for the cases indicated by the run numbers in the figure. The logarithm of Φ inthe unit of W/m3 is shown. The meaning of the gray-scale is indicated by the scale-bar at the bottomof the figure. The region outlined by solid semicircle indicates δs. Note that the distributions in therange 23 λd ≤ x ≤ 52 λd and 0 ≤ z ≤ d/2 are shown in the figure. Parameter values are Tw =1000K and
U =12.6cm/yr.