The plate convergent system is the main driving force for plate tectonics. Some numerical studies have been conducted to reveal how the convergent system was formed. However,the formation of faults that cut through the entire lithosphere has not been adequately described in previous studies. In some studies a weak zone or a fault was set numerically without considering fault initiation (Toth and Gurnis, 1998; Doin and Henry, 2001;Gurnis and Hall, 2004). In this dissertation I solve the problem of fault formation and subduction initiation in the lithosphere using DEM (Discrete Element Method, Distinct Element Method) numerical techniques. DEM numerical methods have great advantages when simulating processes in discontinuous phenomena such as deformation, destruction and faulting. In this method the matrix is composed of many discrete elements, and the deformation of the assemblage is computed by solving the motion equations for each element.The conventional DEM is extended to construct numerical models of a plate convergent system that are composed of layers with different characteristics. Two new models, the "strain model" and the "balloon model" are developed to simulate brittle and ductile layers of the lithosphere. The strain model uses the strain of elements and the Young's modu-lus instead of the compressional distance and the elastic coecient. The balloon model treats each element as a balloon flled by gas, and the assemblage of these elements exhibit a ductile-like character: the yield strength of the assemblage is near constant at various pressures. The characteristics of the assemblages modeled using these extended DEM are tested by numerical geotechnical tests.

The numerical simulations of the development of convergent systems are conducted under four different conditions: a continental plate adjacent to an oceanic plate under a compres-sional stress regime (Experiment 1), a uniform oceanic plate in compressional stress regime (Experiment 2), two adjacent oceanic plates with different thicknesses under a compressional stress regime (Experiment 3), and adjacent continental and oceanic plates under variable stress regimes from extensional to compressional (Experiment 4). The continental plate is modeled as a 3-layered assemblage composed of the following DEM elements: a brittle upper crustal layer, a ductile lower crustal layer and a brittle mantle layer. The oceanic plate is modeled as a 2-layered assemblage: a brittle crustal layer and a brittle mantle layer. The low-strength layers such as the lower continental crust and the asthenosphere are comprised of assemblages of elements simulated using the balloon model, and other brittle layers are simulated using the strain model. Lateral compression and extension in the numerically modeled convergent systems are controlled by lateral motions of the compressional boundary. Plate subduction occurs in Experiment 1, Experiment 3 and Experiment 4. The homogeneous oceanic plate in Experiment 2 shows only thickening and folding, and plate subduction does not occur in this experiment. The forces exerted near the compressional boundary are small in Experiment 4 compared to other experiments.This might be because the extensional stage that occurred before compression in Experiment 4 results in gravitational instability of the slab, and subduction occurs more easily in this experiment.

The plate subduction process is classifed in three stages: the plate boundary is rotated by the density contrast or lateral compression between the modeled plates (Stage 1), the slab is deformed along the inclined plate boundary and the denser plate begins to subduct,intitiating the formation of a mantle wedge (Stage 2), and the denser plate subducts under the other plate continuously, taking crustal materials into the deeper part of the Earth (Stage 3). The results of deformation show that the plate can start to subduct during the initial stages of formation and deformation along the plate boundary, and no other weak zones are needed to form a plate convergent system. The stages of plate subduction might be constrained by the shortening ratio of the convergent system. The internal dynamics and deformation of the slab are analyzed by detecting the stress chain network and the slip events. The stress chain network is shown by the extent of more compressed elements and slip events, and represents the locations of apparent excess shear force compared to frictional force. During plate subduction, the orientation of stress chain network tends to be normal to the slab, and many slip events are detected beneath the plate boundary and upper half of the subducting slab. The slab-normal compressional stress regime of the modeled slabs shows that the internal stress regimes of descending young slabs are under down-dip extension. Such stress regimes in short slabs were shown by Isacks and Molnar (1971): the stress regimes in short slabs are observed predominantly as down-dip extension. The distribution of slip events in the slabs modeled in this study might reveal the characteristics of a Wadati-Benioff zone of a young slab.

本論文は6章から構成されている。第1章はプレート沈み込み初期段階に関する既存研究のレビューである。本章では、1)野外調査では難しい長時間スケールにおけるプレート収束帯の発達過程、2)アナログモデリングでは難しいスケーリングによる影響、3)既存の数値モデルでは難しいすべり面の自発的発達の過程の様相を解明するためにDEM法を導入したと述べられている。

第2章は5項目から構成され、DEM法の原理、本論文で初めて提案された拡張DEM法と計算プログラムさらに、拡張DEM法で採用した要素集合の物性を検証するための仮想一面剪断数値試験とその結果が述べられている。拡張DEM法は、地球内部のbrittleな挙動を反映させるため、要素間の弾性を扱う際に通常用いられるelastic coefficientの代わりにYoung's modulusを計算に用いた要素(strainモデル)と、地球内部のductileな挙動を反映させるため、反発力の増分が要素間距離によって変化する性質を持つ要素(balloonモデル)からなる。本論文では垂直応力に比例して最大剪断応力値も増加することをbrittleな挙動、垂直応力によらず最大剪断応力値は一定であることをductileな挙動、内部摩擦角が小さいほど、要素集合はductileである、と解釈している。

第3章は7項目から構成され、プレート沈み込みシミュレーションのモデルと、その結果について述べられている。大陸プレートはbrittleな上部大陸地殻、ductileな下部大陸地殻、brittleな大陸性マントルから構成され、海洋プレートはbrittleな軽い海洋地殻と冷えた重い海洋マントルから構成されている。大陸プレートと海洋プレートの下にはductileなアセノスフェア(ductileで高密度)が存在する。各層の物性パラメータは仮想一面剪断試験の結果により調整されている。大陸プレートおよび海洋プレートを模した要素集合を並べて収束境界モデルが構築され、これを水平方向に圧縮することによる変形様式がシミュレーションされている。シミュレーションは以下の4つの場合が行われた。1)大陸プレートと海洋プレートが接する場合、2)海洋プレートと海洋プレートが接する場合、3)厚さ(時代)の異なる海洋プレートが接する場合、4)大陸プレートと海洋プレートが一旦伸張しその後圧縮される場合。シミュレーションの結果1)と4)では海洋プレートが、3)では古いプレートの沈み込みが開始されることが判明した。

第4章はシミュレーション結果について議論されている。まず、沈み込みの様式は形態変化から3段階に区分している。第1段階では2つのプレート境界(初期条件では垂直)が沈み込むプレート側に回転する。第2段階では境界で滑りが開始され、沈み込むプレートの先端部に団子状の頭が形成される。第3段階ではプレートの先端が成長するに伴い密度差が増大し沈み込みが促進されていき、上盤側の境界近傍にマントルウェッジが形成される。また、下盤の地殻にdecollementが出現し付加プリズムが発達する。また、1)の場合では海溝地形が発達する。

次にスラブ内の応力分布をstress chain networkで解析している。その結果から、スラブ内の応力はdown-dip extensionであり、スラブ内の断層はプレート沈み込みの方向から共役45°方向に発達すると解釈している。

最後にslip event解析から要素間のすべりが頻発する領域を検出し、slip eventはスラブの沈み込み開始に伴って収束境界下部付近に集中し、その後スラブの伸長に伴い上半分にslip eventを引き込みつつ発達すると解釈している。もし、スリップイベントを地震に相当するものと考えれば、スラブ上半分に現れたスリップイベントはWadati-Benioff帯を表現していると解釈している。

本論文では新たに開発した拡張DEM法を用いてプレート沈み込みの初期過程をシミュレーションで再現している。その結果から密度差を持つプレートが接した境界では、プレート全体を横切るような断層面を前もって仮定しないでも、圧縮に伴うプレート境界面とスラブの変形によりプレートの沈み込みが開始する可能性が示唆されている。

以上のように本論文では将来、様々な展開が期待される手法を開発し、沈み込み開始の興味深い結果を得ていることは極めて高く評価される。

なお、本論文第2章は、山田泰広・徳山英一との共同研究であるが、論文提出者が主体となって理論構築・プログラミングを行ったもので、論文提出者の寄与が十分であると判断する。

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