Topography in plate subduction zones is formed as the result of interaction between the internal crustal processes such as igneous activity and tectonic deformation due to plate subduction and the external surface processes such as erosion and sedimentation controlled by climate and environment. The internal processes cause surface uplift and subsidence, while the external processes even the Earth's surface through mass transfer from highlands to lowlands. Therefore, we can regard the internal and external processes as active and passive processes, respectively. In the present paper we focused on the internal crustal processes in plate subduction zones.

In numerical simulation for the development of mechanical and thermal structure in plate subduction zones, the following points are essential: 1) rational representation of plate-to-plate interaction, 2) realistic modeling of rheological structure, and 3) interaction between mechanical and thermal processes. For the first point, we rationally represented the plate-to-plate interaction by the increase of tangential displacement discontinuity (fault slip) at plate interfaces. For the second point, we modeled the realistic rheological structure of subduction zones with a finite element method (FEM). As to the third point, we constructed a mechanical-thermal interaction model by combining a thermal FEM model for temperature changes due to thermal diffusion and advection with a mechanical FEM model for internal velocity fields due to plate subduction through a temperature-dependent constitutive equation prescribing rheological properties of materials.

In Chapter 2, first, we developed a mechanical FEM model to compute the long-term internal velocity fields for a simple subduction zone, which is modeled by an elastic continental plate and an elastic oceanic plate descending into an underlying viscoelastic half-space. The steady plate subduction is represented by the steady increase of tangential displacement discontinuity at the interface between the continental and oceanic plates (Matsu'ura & Sato, 1989). The most direct way to obtain the internal velocity fields is to solve a boundary value problem for the elastic-viscoelastic composite medium in the time domain, but it is very difficult even in the simple case. In order to avoid the difficulty in computation, we applied the corresponding principle of linear viscoelasticity (Lee, 1955; Radok, 1957) and the equivalent theorem (Fukahata & Matsu'ura, 2006) to the elastic-viscoelastic problem, and reduced it to a simple elastic problem, which is directly obtained from the original elastic-viscoelastic problem by regarding the underlying viscoelastic half-space as the elastic half-space with a normal bulk modulus and a very small rigidity. Then, we can use a standard elastic FEM to construct a mechanical model for computing long-term internal velocity fields due to steady plate subduction in a realistic situation. In the construction of the mechanical FEM model, we represented the tangential displacement discontinuity at the plate interface with the split node technique (Melosh & Raefsky, 1981). The negative and positive buoyancy effects related to surface uplift and subsidence were incorporated into the stress-free conditions at the Earth's surface (Williams & Richardson, 1991). When we decrease the rigidity of the underlying elastic half-space to zero, the computational instability called "volumetric locking" arises. We resolved this computational instability by using the selective-reduced integration scheme (Hughes, 1980).

Next, with the mechanical FEM model we computed the internal velocity fields due to steady plate subduction. We took a 2000 km (horizontal) x 500 km(vertical) fine-mesh model region composed of 800x500 elements as shown in Fig. 1, where the blue thick solid line indicates the split nodes to represent tangential displacement discontinuity (fault slip) at the plate interface. In order to suppress boundary echoes we added a broad coarse-mesh region surrounding the fine-mesh model region. The right side of the total model region is set to be horizontally fixed but vertically free, the base to be vertically fixed but horizontally free, and the left side to be free both for horizontal and vertical displacements.Therefore, the computed results represent the velocity fields relative to a remote reference point fixed to the continental plate. The traction free conditions at the upper surface were modified to incorporate the negative and positive buoyancy effects related to surface uplift and subsidence. From numerical computations we found the following. In the case of a horizontally layered model, the long-term internal velocity fields calculated from the FEM model almost completely agree with those from analytical expressions by Fukahata & Matsu'ura (2006). The pattern of velocity fields in the case of a descending slab model are similar to that in the case of the horizontally layered model, if the thickness of plates is the same in both models, as shown in Fig. 2. A significant difference between these two models is in uplift rates on the overriding plate: the maximum uplift rate for the descending slab model is much faster than that for the horizontally layered model. Then, we examined the effects of the thickness of descending slab on surface velocity patterns. The results are shown in Fig. 3. If the descending slab is thicker, the uplift zone in the overriding plate and the subsidence zone around the trench become broader, and the uplift and subsidence rates become faster. If the overriding plate is thicker, the subsidence pattern around the trench does not change significantly, but the uplift pattern in the overriding plate changes notably. These results are consistent with geophysical observations in subduction zones.

In Chapter 3, first, we developed a thermal FEM model to compute temperature changes due to thermal diffusion and advection. In the development of the thermal FEM model, we considered three sources of advection;the subduction of the cold oceanic slab, the mantle flow induced by slab subduction , and the verticaldeformation of the overriding plate. With the thermal FEM model we c omputed temporalchange i n thermal structure for given internal velocity fields. The computed results in Fig. 4 show that the third source is less effective in heat transfer, but crucial in geomorphic development in plate subduction zones. Next,combining the thermal FEM model with the mechanical FEM model, we constructed a mechanical-thermal interaction model. In this model the mechanical-thermal interaction is described through a temperature-dependent constitutive equation prescribing rheological properties of materials; that is, the temperature-dependence of the mantle viscosity (Courtney & Beaumont, 1983). With the mechanical-thermal interaction model, we numerically simulated the development of mechanical and thermal structure in a plate subduction zone. In this simulation, at a certain time step, we compute internal velocity fields due to plate subduction. Then, by using the computed internal velocity fields, we evaluate temperature changes due to thermal diffusion and advection, and update the boundaries between the lithosphere and the asthenosphere to compute internal velocity fields at the next time step. Through the numerical simulation with such a computation algorithm, we revealed the evolution process of mechanical and thermal structure in the plate subduction zone over a span of 5 Myr. In each diagram of Fig. 5 we show the surface uplift rate (top) and the mechanical and thermal structure development (bottom) together with internal velocity fields. In the early stages of plate subduction (0-2 Myr),the cooling of the mantle wedge l e ads to the thi ckening of the continental lithosphere near the plate boundary and increases the upli ft rates of the lithosphere beneath the island arc. Then, as time goes on, the thinning of the lithosphere beneath the island arc proceeds, and the uplift rates further increase there. This simulation result indicates that the mechanical-thermal interaction is crucial to understand the geomorphic evolution in plate subduction zones.

Fig. 1 The FEM model used for numerical computation.

Fig. 2 Internal velocity fields for a descending slab model.

Fig. 3 The effects of the thickness of descending slab.

Fig. 4 Thermal structure change due to steady plate subduction.

Fig. 5 The development of mechanical and thermal structure.

本論文は全五章からなる。第一章では、先行研究のレヴューおよび本研究での戦略、基本的考え方が述べられている。沈み込み帯では、数千年から数百万年スケールにおいて、沈み込みに伴う力学的作用により流動・変形場が形成され、温度構造が改変される。さらに、温度構造の改変により流動学的な構造が改変され、変形場を改変する。この相互作用が進行する沈み込み帯の変形・温度構造の発達を数値シミュレーションする際には、1)プレート間相互作用の合理的表現、2)物性構造の適切なモデル化、3)力学的過程と熱的過程の相互作用、の3点が重要となる。第一章では、本研究でこれら3点がどのように扱われたかが概観されている:第一の点をプレート境界における変位不連続により合理的に表現し、第二の点を複雑な物性構造を扱える有限要素法 (FEM) によりモデル化した。第三の点については、力学FEMモデルと温度変化を扱う熱的FEMモデルとを、流動学的特性を規定する温度依存構成式により結びつけ、力学的過程と熱的過程の相互作用を扱うカップルモデルを構築した。

第二章では、まず定式化および長期の内部変形速度場を計算する力学FEMモデルが示されている。沈み込みプレート(弾性)が上盤側プレート(弾性)の下に沈みこみ、粘弾性領域であるアセノスフェアへと降下していると仮定する。プレートの定常沈み込みは、プレート境界における変位不連続により表現する。この設定で、変形速度場は、時間領域において弾性・粘弾性の混合媒体で境界値問題を解くことにより求められるが、この解を得るのは一般に非常に難しい。そこで、弾性・粘弾性問題の等価定理を適用し、解きやすい形式に変換することで、プレート沈み込みに伴う内部変形速度場を求めた。その上でFEMを用い、プレート境界における変位不連続をSplit Node 法により表現した。第二章の後半では、この手法が沈み込み帯における長期の変形評価に応用されている。計算条件は、2000 km × 500 km の2次元領域を計算領域として密に分割し、境界の影響を減衰させるため、密分割領域の周囲に粗い分割域を加えた。水平成層構造を仮定した場合、今回構築したFEMモデルにより計算した長期の速度場はほぼFukahata & Matsu'ura (2006) による解析解による計算結果と一致した。弾性領域であるスラブの、上盤側プレート下への沈み込みを考慮したモデルでは、水平成層構造モデルに比べ、上盤側の隆起速度は顕著に速くなった。次に、沈み込みプレートの厚さを変えて計算したところ、プレートが厚くなるにつれ、上盤側の隆起域と海溝付近の沈降域は広がり、それぞれの隆起速度と沈降速度は速くなった。上盤側プレートの厚さを変えた場合は、プレートが厚くなっても沈み込みパターンは変わらないが、上盤側プレートにおける隆起速度は有意に変化した。これらの結果は、実際の地球において、リソスフェアの厚さに対応する海洋底年代と、隆起沈降速度に対応する重力異常とが、対応関係にある可能性を示唆する。

第三章では、前章で述べられた力学的FEMモデルに加えて、熱拡散・移流による温度変化を計算する熱的FEMモデルが構築され、さらに沈み込み帯における力学的-熱的構造の連性発達モデルが構築される。先ず、ある時間におけるプレート沈み込みによる内部変形速度場を計算する。次に、得られた速度場から熱の拡散・移流による温度変化を計算する。さらに、得られた温度場からリソスフェア-アセノスフェア境界を更新し、次の時間ステップにおける内部変形速度場の計算に用いる。このような計算アルゴリズムに基づき、500万年間のプレート沈み込み帯における力学的-熱的構造の発達過程を明らかにした。その結果、沈み込みの開始から時間が経過するにつれ、上盤側隆起域のリソスフェアが薄くなり、この地域の隆起速度が速くなった。このシミュレーション結果は、力学問題と熱問題の相互作用がプレート沈み込み帯における地形発達の重要な支配要因であることを示している。

第四章では、モデルの適用性、限界を踏まえた上で、計算結果と観測(地形、隆起・沈降、重力異常、熱流量等)との対比が行われている。特に、沈み込むプレートの厚さと重力異常の強度・分布の間に明瞭な対応関係が存在し、南米、アリューシャン、日本列島において、古く厚いプレートの沈み込みが、より強い重力異常に対応していることが確かめられた。第五章では、本論文全体のまとめがなされている。

以上のように、これまで取り扱いの難しかった粘弾性体としての沈み込み帯の変形を適切にモデル化し、さらにそれらを熱構造とカップリングさせて解くことに初めて成功した。なお、本論文第二章と第三章は、松浦充宏氏、深畑幸俊氏、中島研吾氏との共同研究であるが、論文提出者が主体となって定式化、モデル構築、解析を行ったもので、論文提出者の寄与が十分であると判断する。

したがって、博士(理学)の学位を授与できると認める。