UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analytic modelling of the deformation due to plate coupling at subduction zones 1995

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


ubc_1996-0046.pdf [ 5.41MB ]
JSON: 1.0052925.json
JSON-LD: 1.0052925+ld.json
RDF/XML (Pretty): 1.0052925.xml
RDF/JSON: 1.0052925+rdf.json
Turtle: 1.0052925+rdf-turtle.txt
N-Triples: 1.0052925+rdf-ntriples.txt

Full Text

A N A L Y T I C M O D E L L I N G O F T H E D E F O R M A T I O N D U E T O P L A T E C O U P L I N G A T S U B D U C T I O N ZONES by JULIAN DOUGLASS B.Sc.E, Queen's University, 1992 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH C O L U M B I A December 1995 © Julian J . Douglass, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ABSTRACT A majority of the largest earthquakes result from the build-up of coupling stresses between tectonic plates at subduction zones. These earthquakes arise, periodically, when the accumulating stress reaches a level which exceeds the rupture strength of the interface. The stresses along the interface may be inferred from geodetic observations ob- tained at the Earth's surface. The process of strain accumulation along subduction boundaries is of particular interest due to the hazard associated with great earth- quakes. Hence, a significant set of geodetic data has been obtained above subduction zones. However, successful interpretation of this data for the magnitude and distribu- tion of stress along the interface requires that a suitable physical model of the strain accumulation process be available. Here, two approaches to this modelling problem are considered. The first is a kinematic approach to the problem in which dislocation theory is used to represent the stick-slip cycle. According to this representation, the deformation due to locking is constrained by the kinematics of plate convergence, with the most important parameter being the plate convergence rate. Hence, an assumption implicit to models that are based on this representation is that the dynamics of subduction can be physically represented by a model constrained only by kinematic parameters. To evaluate this assumption, model solutions are calculated for the case of an infinite half-space Earth. New analytic expressions are derived for this case which allow for calculation of the stress distribution over the complete half-space. The model solutions reveal a subsurface stress field which is difficult to justify on physical grounds. These difficulties raise questions regarding previous interpretations of geodetic data which were based on this model. By exclusively using this model, other possible dynamic scenarios, other than that which is prescribed by the kinematic approach, are excluded from consideration by the interpreter. ii Given the limitations of the dislocation represention, a second, alternative ap- proach to the problem is proposed. This approach is based on a dynamic representation of the process of strain build-up. The overthrust wedge is treated as a free body with tractions applied along its surface. The deformation is then calculated for a prescribed distribution of force along the plate interface. A uniform elastic rheology is assumed. Thus, the fundamental elasticity problem to be solved is that of an infinite elastic wedge subjected to traction boundary conditions. Solutions are obtained by two different, but complementary techniques. In the first, the wedge is conformally mapped to an infinite strip and solutions obtained by Fourier transform methods. Equivalent solutions are also obtained by applying the Mellin transform directly to the biharmonic equation for the Airy stress function. The results of these calculations show that entirely different distributions of basal tractions can lead to very similar deformation signals along the top surface of the wedge. Thus, an ambiguity in the interpretation of the geodetic data is confirmed which is not accounted for in deformation predictions based on the dislocation approach. iii T A B L E O F C O N T E N T S Abstract i i Table of Contents iv List of Tables v i i List of Figures vi i i Chapter 1. I N T R O D U C T I O N 1 1.1 A silence disturbed by the occasional loud bang? 1 1.2 Geophysical data 2 1.2.1 Seismic constraints 2 1.2.2 Thermal constraints 3 1.2.3 Geodetic constraints 4 1.3 Thesis scope 4 Chapter 2. M O D E L L I N G I S S U E S 7 2.1 Introduction 7 2.2 Stick-slip subduction 7 2.3 The geodetic data set 9 2.3.1 Uplift measurements 10 2.3.2 Strain measurements 11 2.4 Model geometry and rheology 11 2.4.1 Boundary constraints 13 2.5 Mathematical approach 13 Chapter 3. K I N E M A T I C A P P R O A C H : T H E D I S L O C A T I O N M O D E L 15 3.1 Introduction 15 3.2 A dislocation model of the earthquake cycle 15 3.2.1 Edge dislocation in the elastic half-plane 19 3.2.2 Calculation of stress and displacement 20 iv 3.3 The stress field 22 3.4 Discussion 23 Chapter 4. T H E D Y N A M I C A P P R O A C H 29 4.1 Introduction 29 4.2 Conceptual representation 30 4.3 Traction at the base of the overthrust 32 4.4 The stress field within the elastic wedge 33 4.4.1 Green's functions 34 4.4.2 Variation of surface strain with depth of applied force 35 4.4.3 Dynamic implications of the dislocation 38 4.5 Discussion 38 Chapter 5. C O M P L E X V A R I A B L E A P P R O A C H T O 2D E L A S T I C I T Y 4 1 5.1 Introduction 41 5.2 Plane strain 41 5.3 Airy stress function 43 5.4 General solution to the biharmonic equation 45 5.5 Stresses and displacements in terms of the analytic functions 47 5.6 Force on a contour 49 5.7 The general problem 51 5.8 Application of the method 53 Chapter 6. D I S L O C A T I O N I N T H E E L A S T I C H A L F P L A N E 56 6.1 Dislocations in continuous media 56 6.2 Dislocations and branch cuts 57 6.3 Edge dislocation in the infinite plane 58 6.4 Edge dislocation in the infinite half-plane 61 Chapter 7. F I R S T F U N D A M E N T A L P R O B L E M F O R T H E W E D G E 64 7.1 Introduction 64 v 7.2 Statement of the problem 65 7.3 Complex variable approach 67 7.3.1 Solution by Fourier transform 68 7.3.2 Mapping and transformation of the boundary conditions . . . . 68 7.3.3 Symmetric normal stresses with zero shear on the boundary . . 71 7.4 Solution by Mellin transform 74 7.4.1 The biharmonic equation in the Mellin domain 75 7.4.2 Symmetric normal stresses with zero shear II 78 7.5 The inverse transform 80 7.6 A note on the wedge paradox 83 Chapter 8. SUMMARY 85 8.1 Modelling of deformation above a subduction zone 85 8.2 Theory of two-dimensional elasticity 87 REFERENCES . . 88 vi LIST OF TABLES 7.1. Boundary conditions for given symmetry cases 71 7.2. Function definitions . . 73 vii L I S T O F F I G U R E S 2.1. Forces acting on lithospheric plates 8 2.2. Simplified subduction geometry 12 3.1. The dislocation representation 16 3.2. Solution domain for the dislocation model 18 3.3. Distributions of slip deficit 18 3.4. Subsurface stress field 21 3.5. Shear traction down the plate interface 23 3.6. Subsurface stress field for a model with transition zone 26 4.1. Formulation of the dynamic model 31 4.2. Point forces applied to the base of the wedge . . 34 4.3. Strain profiles for forces at different depths 36 4.4. Tractions on the plate interface for the dislocation model 37 4.6. Surface strain for opposing shear tractions 39 5.1. Contour within an elastic domain 50 5.2. A multi-connected domain 53 6.1. An elastic domain with a cut 57 6.2. Edge dislocation in the z-plane and the transform plane 59 7.1. The oo elastic wedge 66 7.2. The wedge transformed to the strip 69 7.3. Poles in the plane of the transform variable 80 viii Chapter 1 INTRODUCTION 1 .1 A silence disturbed by the occasional loud bang On the 22nd of May, 1960, a huge earthquake occurred off the coast of Chile. The main event had a moment magnitude over 9, and was followed by numerous large aftershocks. The resulting ground motion produced severe shaking in cities such as Concepcion and Valdivia, and left many low-lying areas inundated as a result of ground subsidence. The earthquake also caused ground liquefaction in some areas, produced seiches in numerous lakes, and sent a severe Tsunami wave off to wreak havoc in coastal Chile, Hawaii and even Japan [Eiby, 1980, pl72-174]. This event is an example of a megathrust earthquake. Such events periodically rock plate margins where one plate is being subducted underneath another. In doing so, they act as dramatic indicators as to the degree of coupling which must take place between convergent plates. Megathrust events are thought to arise as a result of stick-slip sliding along the plate interface. According to the stick-slip model, strain accumulates around some portion of the plate boundary which remains locked until the stress supported at the boundary reaches some critical value. At this point seismic rupture occurs. The cycle then begins anew. The period required for sufficient strain to accumulate in order to precipitate a great earthquake is believed to be between one-hundred and one-thousand years. Due to their immense magnitude, megathrust earthquakes pose a significant seis- mic hazard to nearby (and far-off, when tsunamis are generated) human populations. As a result, society is keenly interested in answering two questions: When will the big one occur? and How big will it be? These are non-trivial questions. The answers 1 Chapter 1. INTRODUCTION 2 depend on a large number of physical variables, many of which are poorly understood. However, the geophysical information available to workers studying such questions is expanding in both breadth and depth, providing workers with better means to constrain both the extent of the coupled portion of a subduction plate boundary and the degree of coupling. This is allowing for ever more precise estimates of potential megathrust magnitude and repeat time. Geodetic measurements of crustal deformation obtained above subduction zones are one geophysical data set that is employed in these efforts. However, in order for this, or any other, data set to be successfully interpreted, a reasonable physical model must be available with which to relate the data to the physical parameters of interest. In the case of subduction zones, the parameters of interest are the extent of the coupled portion of the plate boundary and the rate at which strain accumulates over this region. In this thesis, consideration is given to the evaluation and development of physical models applicable to the interpretation of deformation measurements acquired above subduction zones. 1.2 Geophysical data A variety of geophysical information is employed in the study of subduction zones. Of particular importance is the information obtained by seismic methods, from heat flow measurements, and, as mentioned above, from geodetic observations. 1.2.1 Seismic constraints Seismic data provide essential information on subduction zone geometry and rhe- ology. Various types of seismic information are employed in subduction zone studies. First, locations of earthquake hypocentres can aid in delineating the geometry of the subducting slab, since the slab and slab/continent boundary regions act as the source of many small earthquakes [Hasegawa, 1994]. In addition, slab mapping can be performed Chapter 1. INTRODUCTION 3 by inversion of delays in teleseismic arrival times. Finally, invaluable information is pro- vided by controlled source surveys which cross subduction margins. Such experiments have been performed off the coasts of British Columbia [Drew and Clowes, 1990], and Japan[Moore et. al., 1990]. Seismic data provide invaluable constraints on many model parameters, including elastic constants and the dip angle of the plate interface. This information can then be used as a priori input in formulating dynamic models of sub- duction behaviour. Without this information, models would be too poorly constrained by the available thermal and/or deformation data. 1.2.2 Thermal constraints Thermal information from heat flow measurements may be used to constrain the maximum and minimum depths of the locked zone. As rheological properties change with depth, the degree of coupling which can be supported at the plate interface must vary downdip. At very shallow depth, interplate coupling is thought to be minimal since earthquake rupture zones do not generally extend all the way to the trench axis [Byrne et al., 1988]. This lack of coupling is attributed to the presence of unconsolidated sediments and is partially thermally controlled [Hyndman and Wang, 1993], although the exact physical processes in the shallow region are still not well understood. At greater depth, the maximum nucleation depth of interplate earthquakes is limited by temperature induced rheological changes. Events arising due to stick-slip behaviour are unlikely to initiate below depths at which the temperature exceeds 300°C [Hyndman and Wang, 1993]. Combining the heat flow information with an appropriate thermal model, the extent of the locked zone at a particular subduction zone may be estimated [Hyndman and Wang, 1993, Hyndman et al., 1995]. / Chapter 1. INTRODUCTION 4 1.2.3 Geodetic constraints The geophysical data which provide the most direct indication of plate coupling behaviour are the set of geodetic observations of strain rate and displacement taken above subduction zones. Zones at which such survey campaigns have been undertaken include the Cascadia, Aleutian and Nankai margins (see section 2.3 for references). The data employed include levelling records, tide gauge data, gravity data, and position- ing information obtained using space geodetic techniques such as very long baseline interferometry (VLBI) and the global positioning system [Hager et al., 1991]. 1.3 Thesis scope The interpretation of geodetic observations requires a physical model with which to relate the data to the coupling characteristics of the plate boundary. Two criteria must be fulfilled by the model. First, the model must realistically simulate the dynamics of the subduction process. In other words, a model must predict coupling forces along the plate interface which are physically meaningful. Second, the model parametrization must be simple, since the the geodetic data is sparse, and cannot be expected to constrain large numbers of unknowns; i.e. the number of free parameters should be minimized as much as possible. Unfortunately, these two criteria are, to a certain extent, mutually exclusive, as, the more assumptions are made to simplify a model, the farther the model is moved away from the reality of the physical situation. The objective of this thesis is to formulate a model of the strain accumulation between great earthquakes at a subduction boundary which satisfies, as far as possible, both the above criteria. The thesis is divided into two main parts which are thematically linked, but quite different in content. In the first part, which includes Chapters 2, 3, and 4, the geophys- ical aspects of the modelling study are reported. Included in this part are solutions Chapter 1. INTRODUCTION 5 calculated for two different deformation models, along with discussions of the physical implications of the results. The second part of the thesis, which includes chapters 5, 6, and 7, comprises mathematical solutions to a number of interesting problems in two-dimensional elas- ticity. These were obtained in the course of the geophysical study. However, in order to maintain continuity in the discussion of the geophysics, they are treated separately. Before embarking on any modelling study, one should define the constraints under which the model must operate. This is done in chapter 2, which addresses issues relevant to the modelling problem. In particular, the physics of the stick-slip process is described as precisely as possible given the current level of understanding of subduction dynamics. In addition, consideration is given to the extent of the data set available for interpretation. Finally, general assumptions are introduced which are adopted in the models discussed in chapters 3 and 4. Dislocation models are the current state of the art in models of interseismic defor- mation at stick-slip subduction zones. These models are based on a kinematic represen- tation of the subduction process. In this representation, a dislocation (or discontinuity in displacement) is located along the seismogenic portion of the plate boundary. Strain is then built up by increasing the magnitude of the dislocation at the plate convergence rate. An assumption which is implicit to the representation is that the dynamics of sub- duction may be physically represented by a model whose primary constraint (i.e. the plate convergence rate) is purely kinematic. Many previous studies have included in- terpretations of geodetic observations that are partially based on comparisons of data with results from dislocation models. The prevalence of this model in the literature has motivated a quantitative analysis of a dislocation model in an elastic half space that is presented in Chapter 3. New analytic expressions for this case are derived which allow for calculation of the stress distribution over the complete half-space. Model solutions expose a subsurface stress field that is physically unrealistic. Furthermore, by physical Chapter 1. INTRODUCTION 6 arguments, it is demonstrated that the observed flaws in the stress field cannot be attributed to geometric and rheological assumptions of this particular model. Instead they are features inherent in the dislocation approach. Given the doubts raised as to the validity of the assumptions of the dislocation approach, attention is turned to the development of an alternative approach in chapter 4. A representation of the stick-slip subduction is proposed based on dynamic, as op- posed to kinematic, constraints. As in the previous chapter, a uniform elastic rheology is assumed. The boundary conditions employed are based on an assumed distribution of force along the plate interface. Solutions calculated for this model show that very similar profiles of strain along the surface may be obtained given very different distri- butions of traction along the base of the wedge. Thus, an ambiguity inherent to the interpretation of surface deformation data is demonstrated. The zone of deformation at a subduction zone is generally much longer, along strike, than it is wide. Hence, generally the deformation can be assumed to be invariant along strike. Over relatively short interseismic intervals, an elastic rheology is a suitable assumption. Thus, the theory of two-dimensional elasticity can be employed to obtain model solutions. This is the case for the models discussed in chapters 3 and 4. The methods used to obtain these solutions are the subject of the second part. Chapter 5 contains an introduction to the essentials of the complex variable ap- proach to 2D elasticity. In particular, it is shown that the stresses and displacements can be expressed in terms of two complex functions that are analytic in the elastic domain. These expressions are then applied, in chapters 6 and 7, to two problems relevant to the subduction models discussed in the first part. In chapter 6, the equa- tions of elasticity are solved for an edge dislocation in a half-space, while chapter 7 is devoted to solving for the deformation within an elastic wedge of infinite dimension that is subjected to traction boundary conditions on its surface. Major results of the thesis are summarized in chapter 8. Chapter 2 M O D E L L I N G I S S U E S 2.1 Introduction Quantitative models of physical phenomena accomplish two main objectives. They provide insight into the physics of the phenomenon of interest, and they act as a math- ematical link between experimental data and model parameters. The most appropriate modelling approach is best decided upon after considering a number of issues. First, a conceptual physical model of the process to be modelled must be established. Subse- quently, the quality of the available data must be assessed, focusing particular attention on the data's spatial and temporal resolution. Once these tasks have been accomplished, appropriate simplifying assumptions may be implemented, and a mathematical model may be formulated. In this chapter, these issues are considered as they pertain to the problem of interest in this thesis: the formulation of a model simulating interseismic deformation at a subduction zone. 2.2 Stick-slip subduction A conceptual model of the dynamics of plate coupling must be established in order to model the strain accumulation at a stick-slip boundary. The convergent plates are driven towards each other by some combination of remote tectonic forces, that are counterbalanced by a combination of reaction forces at the plate boundary. A number of plate driving forces have been postulated [Forsyth and Uyeda, 1975], including slab pull, ridge push, trench suction, and various drag forces (figure 2.1). The relative importance of these forces is still a subject of debate [Stefanick and Jurdy, 1992; Shimenda, 1993]. Fortunately, however, for the purpose of simulating interseismic deformation near the 7 (b) Figure 2 .1 Forces acting on lithospheric plates: (a) ridge push, (b) slab pull, (c) mantle drag, (d) trench suction, (e) friction re- sistance at transform faults [Forsyth and Uyeda, 1975]. subduction zone, the exact nature of the plate driving forces is unimportant; these forces can be considered invariant over the small spatial and time scales associated with the earthquake cycle. A conceptual model of the distribution of reaction forces along the plate inter- face is then required to model interseismic strain accumulation. Assuming that the net tectonic force is constant, it follows that the net reaction force must also remain constant. The net reaction force comprises two components: i) friction forces which are associated with plate locking and ii) velocity dependent viscous drag forces. The balance between these two components varies over the course of an earthquake cycle. It is this variation that produces the observable deformations. At the beginning of the seismic cycle (directly following the previous megathrust event), the dynamically supported viscous stress is at its maximum. Over the course of the interseismic period, this stress relaxes and the elastic stress increases. Thus, the effect of the locked zone is to produce a shear traction (or tug) along the locked zone in the direction of motion of the other plate. Seismic rupture occurs when the magnitude of this shear traction reaches the breaking strength of the fault. Chapter 2. MODELLING ISSUES 9 Hence, the deformation observed above the subduction zone is caused by the com- bined effects of two processes: i) the build-up of elastic strain due to locking, and ii) the relaxation of viscous stress caused by the previous megathrust earthquake. The surface deformation signature associated with elastic shear is likely to be the more important of the two processes. This is a consequence of the depth dependence of rock rheology. At shallow depth, pressure and temperature are low and crustal rocks behave in a brittle manner. Hence, the rheology can be assumed to be perfectly elastic. At greater depths, and as the pressure and temperature increase, the rocks begin to evidence more viscous behaviour. Hence, elastic coupling between convergent plates occurs near the subduction trench, whereas coupling at depth is predominantly viscous. As a consequence, by a simple application of Saint Venant's principle, and its analog for viscous flow, the surface deformation signature associated with viscous relaxation should be significantly smaller then that assosciated with strain accumulation along a locked zone. Based on previously reported laboratory data and field estimates for crustal rocks, Hyndman and Wang [1993] estimate a maximum temperature of 300-350°C above which a subduction plate boundary can no longer remain locked. Using a thermal model constrained by heat flow measurements, they then conclude that the maximum depth of the locked zone is approximately twelve kilometres at the Cascadia margin [Hyndman and Wang, 1993], and is approximately twenty kilometres at the slightly cooler Nankai margin [Hyndman et al., 1995]. Such estimates further suggest that, to first order, the build up of stress along the locked zone is the most important influence on surface deformation. 2.3 The geodetic data set Geodetic data are available from a number of coastal regions adjacent to subduc- tion fronts. These areas include the south coast of Alaska, southwest Japan and points Chapter 2. MODELLING ISSUES 10 along the west coast of North America near the Cascadia subduction zone. The data comes from repeat levelling surveys, tide gauge records, gravity surveys and horizontal strain networks. Two quantities are inferred from the data: rate of uplift and rate of accumulation of horizontal strain. 2.3.1 Uplift measurements The most detailed data set has been accumulated from the southwest coast of Japan. Some tide gauge records extend back to the late 1800s, and uplift records have been determined for at least fifteen stations for the interseismic period from 1951-1990 [Savage and Thatcher, 1992, Savage, 1995]. In addition, levelling lines exist which have been resurveyed up to five times since the 1890's [Thatcher, 1984]. The data are sufficient to infer mean uplift rates over four time intervals and four separate routes running roughly perpendicular to the trench. During one of these time intervals, two megathrust events occurred; one in 1944 and one in 1946 [Ando, 1975]. Hence, co- seismic as well as interseismic deformation data are available. A significant set of geodetic data also exists for the Cascadia margin [Dragert et al., 1994]. Uplift data are available from three levelling lines running roughly perpendicular to the trench, with each line having been surveyed at least twice over the last thirty years. Hence, mean uplift rates for that period may be calculated. In addition, tide gauge records spanning periods of up to 80 years (in the case of Victoria, B.C.) provide additional uplift data. Finally, the data set is complemented by repeated, high-precision gravimeter surveys which were conducted along one of the levelling lines. The gravity data provide a further means to infer uplift along the survey line. Uplift data has also been analyzed from Alaska. However, the data set is much sparser than that available from the Nankai and Cascadia margins. Savage and Plafker [1991] calculated uplift rates from thirteen tide gauge stations scattered along the south Chapter 2. MODELLING ISSUES 11 coast of the state from 135°W to 160°W. Beyond these records, little uplift information has been reported. 2.S.2 Strain measurements Few observations of horizontal strain rates have been reported. The largest data set comes, once again, from Japan where numerous strain networks have been surveyed near the Nankai subduction zone [Savage and Thatcher, 1992]. Four strain measurements have been recorded along the Cascadia margin [Dragert et al., 1994], and one in the vicinity of the Shumagin Islands of Alaska [Larson and Lisowski, 1994, Lisowski et al., 1988]. However, the density of available strain data is likely to improve substantially in the near future. At the time of writing, Japan had over 200 continuously tracking GPS antennas scattered over the country [Geophysical Survey Institute of Japan, World Wide Web page, 1995]. Initiatives had also been started in Chile to measure strain along the west coast of South America [Angermann et al., 1994], and extensive campaigns to measure strain along segments of the south coast of Alaska were underway [Sauber and Cohen, personal communication, 1995]. Satellite geodesy is rendering feasible relative positioning surveys that would have been impossible given the financial and logistical constraints of traditional survey techniques. 2.4 Model geometry and rheology The above data set is sparse even along the most robustly surveyed margins. Hence, a number of general assumptions are warranted. These assumptions are applicable to each of the models described in the succeeding chapters. First, and perhaps most im- portant, the assumption is made that surface deformation due to relaxation of viscous stress at depth is negligible. Hence, a purely elastic rheology is implemented. Further- more, any variation in the elastic constants, shear and bulk modulus, over the model Chapter 2. MODELLING ISSUES 12 domain is unlikely to be large enough to substantially alter predicted deformation pat- terns. Therefore, they are assumed constant throughout the model domain. Second, a simplified subduction geometry is assumed. The most important geo- metric parameter is the dip of the locked zone. The dip angle is generally small at stick-slip margins, often being less than 10°. The dip also varies with depth, increas- ing gradually until a depth of approximately 30 km, and then exhibiting a kink which increases the dip angle substantially, often to greater than 30°. Here, the dip angle is assumed constant, and any kink in the plate boundary is presumed to occur at depths well below the loci of elastic stress build-up along the locked zone. Considering the locked zone depth extents estimated by Hyndman et al. [1995], this assumption is well within reason. Upon implementing these assumptions, a simplified model geometry is obtained (figure 2.2). The free parameters are those associated with the extent of the locked zone and the degree of coupling between the plates. These parameter magnitudes are the unknowns to be determined by comparison of modelling results with the available geodetic data. However, the exact nature of these parameters depends on both the nature of the boundary constraints imposed on the model, and on the mathematical approach selected. Figure 2.2 Geometry of simplified subduction model with (a) and (b) indicating the overriding plate and subducting slab re- spectively, (c) is the dip angle of the interface. The hatching indicates the locked zone. Chapter 2. MODELLING ISSUES 13 2.4-1 Boundary constraints As the driving mechanism behind the subduction process is not well understood, physical assumptions must be imposed in order to establish far-field boundary condi- tions for any subduction model. Two approaches are usually taken. One is a kinematic approach in which the rate of plate convergence, a known quantity, is specified. A representation of the earthquake cycle is then formulated that maintains this rate of convergence. The alternative is to approach the problem from a dynamic perspective. This approach is based on the assumption that some unknown force is acting in the far field. This force is balanced by reaction or collision forces acting at the plate boundary. The distribution of these reaction forces is assumed to vary during the earthquake cy- cle, thus producing the observed surface deformations. Hence, the far-field, or driving forces, do not need to be specified. Models based on each of the above approaches are considered in Chapters 3 and 4 of this thesis, respectively. 2.5 Mathematical approach An important aspect of any modelling problem is the selection of a mathematical approach. The available options range from pure analysis to pure numerics. The final choice of approach depends on two main factors: the level of model complexity required (this is highly dependent on the size of the data set), and the level of sophistication reached in previous modelling studies. In this case, an analytic approach is justified. This approach has the disadvantage of permitting only very simple model parametriza- tions. However, the significance of this difficulty is tempered by the limited size of the data set. In addition, few previous models have been developed to investigate subduction dynamics. Moreover, analytic models are particularly useful for analyzing relatively new types of geophysical data, such as the unprecedented geodetic data sets presently becoming available. Hence, for the subduction problem, significant benefit Chapter 2. MODELLING ISSUES may be anticipated from the additional solutions. 14 insight that is gained by pursuing analytic Chapter 3 K I N E M A T I C A P P R O A C H : D I S L O C A T I O N M O D E L S 3.1 Introduction f Significant effort has been devoted to the development of kinematic models of sub- duction deformation. Such models employ a dislocation representation of the earth- quake cycle [Savage, 1983, Maatsu'ura and Sato, 1989], and have been used extensively in efforts to interpret various geodetic data sets [Cohen, 1994; Cohen et al., 1995; Dragert et al, 1994; Hyndman and Wang, 1993; Hyndman et al., 1995; Lisowski et al., 1988; Savage, 1983,1995; Savage and Blajker, 1991; Savage and Thatcher, 1992; Savage et al, 1991; Thatcher and Rundle, 1984]. However, kinematic (or dislocation) models do not take into consideration the dynamics of the plate interaction, despite being de- signed to simulate deformation which arises as a direct result of this interaction. In effect, the dynamics of the plate interaction is imposed by the kinematic representation. In this chapter, the dynamics implied by this representation are investigated; motivated by the importance of dislocation models in the literature on subduction deformation. 3.2 A dislocation model of the earthquake cycle The investigation begins with a quantitative analysis of a simple dislocation model. The model makes use of the dislocation representation of the earthquake cycle proposed by Savage [1983]. According to this representation, strain accumulation at a subduction zone is equivalent to the linear superposition of solutions to two subsidiary problems (figure 3.1). The first problem is that of steady reverse slip along the plate bound- ary. This is the steady state situation whose solution is equivalent to the long term \ A shorter version of this chapter was published by Douglass and Buffett [1995] 15 Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 16 Figure 3.1 According to the dislocation representation, locking of the seismogenic zone (as depicted in (c)) is achieved by the superposition of a dislocation oriented in the back slip direction (b), on a state of steady slip (a). deformation signature. For the second problem, an edge dislocation, oriented in the direction opposite to that of the relative motion of the two plates, is imposed along and parallel to the locked portion of the interface. The magnitude of this dislocation is often termed the slip deficit. A no slip condition can be maintained by increasing the magnitude of the dislocation at the plate convergence rate. Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 17 Generally the rate of strain accumulation associated with steady state subduction is assumed negligible in comparison to that associated with the accumulation of slip deficit. Strain accumulation during the steady state process is not released by the model's strain release mechanism (i.e. by seismic rupture). Hence, this strain accu- mulates over a subduction zone's entire history. Thus, it is argued that, if this steady accumulation were significant, huge deformations would be observed within the over- riding plate [Savage, 1983]. Accepting the above assumption, the strain accumulation problem is reduced to a single boundary-value problem. The geometry of the problem is illustrated in figure 3.2. The domain consists of the entire two dimensional half-space, excluding an infinitely thin strip coincident with the dislocation. A solution is obtained by satisfying the equations of elasticity over this domain given the boundary conditions: a free surface along the top of the half space, a function specifying the displacement discontinuity as a function of depth along the locked zone, and conditions which maintain the finiteness of the stresses and displacements at infinity. The parameters of the problem are the slab dip angle, the position and length of the locked zone, the distribution of slip deficit along this zone, and any geometric and rheological parameters associated with the chosen Earth model. This Earth model is obtained upon implementing the geometric and rheological as- sumptions described in section 2.4. The mathematical problem is to calculate the strain build-up due to a distribution of edge dislocations within an elastic half space. The problem is most effectively tackled using a Green's function approach. First, a solution is sought to the problem of a point dislocation within an elastic half space. Solutions for more complicated distributions are then obtained by evaluating the integral, il(x,y)= f to(x,y;s)b'(s)ds (3.1) Jo Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 18 \ \ ^0 ^ ^ ^ ^ s ^ i i ^ ^ ^ ^ • y Figure 3.2 Domain of solution for a dislocation model of subduc- tion deformation. The thick dot-dashed lines delineate boundaries of the solution domain. The thick solid line indicates the location of the displacement discontinuity. where fl(x,y) is the stress, strain or displacement component as a function of posi- tion, S7(x, y;s) is the corresponding solution for an edge dislocation of unit magnitude, extending from the origin to a distance, s, down the plate interface, and b'(s) is the derivative of the dislocation magnitude as a function of distance from the surface down the plate interface. (a) o CD "D 1 0 CL (b) 0 1 Figure 3.3 Magnitude of slip deficit as a function of distance, s, down the plate interface for the two models considered in this study. Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 19 3.2.1 Edge dislocation in the elastic half-space The problem of a uniform edge dislocation, extending from the surface to some finite depth in an elastic half spaces, is discussed at length in chapter 6. Only the solution is given here. The solution is expressed in terms of two complex potentials, <f> and ip, each of which is a function of the complex variable, z — x + iy. The coordinate system is defined with the x and y axes extending along the surface and vertically downward respectively. The origin is located at the point of intersection of the plate interface and the half-space boundary. According to the general formulae of the complex variable method of two-dimensional elasticity (q.v. Chapter 5), <TXX+<Tyy = 2[<f>'(z)+<p'(z)], (3.2) <Tyy + axx + 2i<rxy = 2[z<j>"(z) + i>\z)], (3.3) and 2u.(u + iv) = K<f>(z) - z<t>'{z) - i/>(z), (3.4) where the crap are components of the stress tensor, u and v are the horizontal and vertical displacements respectively, K — Z — AU where v is Poisson's ratio, and fi is the bulk modulus. The potential functions which apply in the case of an edge dislocation are given by, 'w-^+fH? 1 ' _ (3-5) *•« = -(TO)+ (^k + t ^ ( B ( z " 1 3 + S ( z + *> ] - ( 3 ' 6 ) where B is the magnitude of the displacement discontinuity and za refers to the location of the downdip end of this discontinuity. Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 20 3.2.2 Calculation of stress and displacement Given the expressions (3.1) to (3.6), the mathematical problem is essentially solved. However, for practical applications these expressions must be reduced to their real and imaginary parts in order to extract the physical quantities of interest: the stress and displacement components. In general, this is a non-trivial task. The operation entails the reduction of lengthy rational expressions of a complex variable. The algebra is tractable for the case, y = 0. Hence, expressions for the strain and displacement components at the free surface are available [Freund and Barnett, 1976, Savage, 1983]. However, the exclusive use of these expressions does not allow a proper assessment of the dynamic implications of the dislocation representation of the subduction process. Hence, for this study, the expressions, (3.2-6), are also reduced for the general case of arbitrary x and y, thus allowing for the evaluation of the stress components over the full two-dimensional domain. The reduction is rendered feasible by the availability of modern analytic software tools. In this case, a routine was written in the Mathematica language [Wolfram Research, 1992] to reduce rational, complex expressions. The code was then implemented, and expressions were obtained for the three components of the stress tensor. The length of these expressions, being summations of roughly sixty rational expressions, precludes their inclusion in this thesis. However, the character of the solutions is illustrated by the examples discussed below. Figure 3.4 (following page) Subsurface stress field produced by a point dislocation at s=l, i.e. that associated with the distribution of slip deficit illustrated in figure 3.3a. The three plots, (a), (b) and (c), correspond to the three components of the stress tensor, crxt + o-yy, axx — ayy and axy. The colour scale is linear with deepening red and blue indicating positive and negative values respectively. Note that the scale is chosen so as to exaggerate the stress distribution away from the location of the stress singularity, and that the strain energy is largely confined to the dark blue and red areas. Chapter S. KINEMATIC APPROACH: DISLOCATION MODELS 21 1 2 3 4 5 Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 22 3.3 The stress field Dislocation models that are typically used to simulate the strain accumulation at a subduction zone are typically assigned a fairly simple distribution of slip deficit along the locked zone. Here, solutions are obtained for two such distributions. The various coordinates employed in the analyses are illustrated in figure 3.2. In the first example, the deformation is calculated for a locked zone that extends from the surface to some finite depth (see figure 3.3a). The slip deficit is assumed uniform along the locked zone. Hence, b'(s) is simply a Dirac-delta function located at the downdip end of the locked zone. It follows from (3.1) that the deformation occurs as a direct result of this jump in the slip deficit. In figure 3.4, the components of the stress tensor are plotted as a function of position within the half space. The most notable feature of these plots is the distinct lobed pattern which emanates from the stress singularity at the downdip end of the locked zone. The lobes indicate stresses which alternate in sign about the slip-deficit jump at the base of the locked zone. The magnitude of the shear traction along the base of the overriding plate is shown in figure 3.5. As anticipated based on the stress field plots, this plot shows change of sign in the shear traction at the base of the locked zone. It is also notable that the shear is actually oriented updip below this zone. Such an orientation, which is opposite to the direction of relative motion of the convergent plates, is very difficult to justify on physical grounds. In addition, the net shear supported by the locked zone is of the same magnitude as that supported outside this zone. This violates a starting assumption of the conceptual model of stick-slip subduction described in the previous chapter; i.e. that relatively little shear is supported elastically below the locked zone. Evidently, the singular nature of the stress field associated with the above solu- tion plays an important role in determining the deformation pattern. However, this singularity is a mathematical artifact of the chosen slip-distribution, resulting from the abrupt change in slip deficit at the base of the locked zone. Because the presence of Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 23 the singularity is a consequence of the simplified model parametrization, it can not be considered a major flaw in the model. Indeed, many authors have argued that there must be a transition zone between the fully locked and freely slipping portion of the interface, and have employed models incorporating such zones in order to fit the model data to surface observations [Lisowski et al., 1988, Hyndamn and Wang, 1993]. I l l Figure 3.5 Plot of the shear traction along the plate interface. The asymptote in the <rr$ profile corresponds to the location of the displacement singularity. The second example incorporates such a transition zone into the model (figure 3.3b). The amount of slip deficit is gradually reduced from its maximal value to zero. However, the solution for this case, shown in figure 3.6, exhibits the same lobed stress pattern as is observed in figure 3.4. Significant interplate traction is evidenced below the locked zone, and more significantly, the orientation of the shear traction continues to be updip along a portion of the plate boundary. 3.4 Discussion The examples presented above provide quantitative evidence of the problems with dislocation models of subduction behaviour. The reasons for these difficulties become Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 24 clear upon examination of the dynamic constraints imposed by the dislocation rep- resentation of subduction behaviour. According to this representation, deformation is produced by prescribing kinematic conditions along the plate interface. No force is applied anywhere along the boundary of the solution domain: the inner and outer contours of figure 3.2. Hence, the net traction must be zero along any internal surface which divides the domain into two parts. The plate interface is such an internal surface. As a result, any interplate traction that occurs along a segment of the plate boundary must be balanced by some opposing traction acting outside this zone. This has two important consequences. First, the tractions must have opposing orientations along different portions of the interface, and, second, any net traction along the locked zone must be counterbalanced by equal and opposite tractions acting outside this zone. The latter consequence is a violation of the starting assumption that coupling be confined to the locked zone, whereas the former implies that interplate shear must act in the direction opposite to that of the relative motion of the two plates along some portion of the plate boundary. Both features are observed in figure 3.5. The above physical arguments also serve to explain the lobed stress patterns observed in figures 3.4 and 3.6. Further shortcomings of the dislocation formulation are revealed in the limit of a very strong locked zone (or, equivalently, a very long interseismic period). As discussed in section 2.2, tectonic plates converge under the influence of some combination of remote tectonic forces. These forces are balanced by reaction forces acting along the plate interface, which can be divided into net shear and normal components. The shear component is supported by a combination of static coupling along the locked zone and dynamic or viscous forces acting downdip of this zone. Hence, if the locked zone were infinitely strong, the slab would eventually, with the build-up of elastic strain, be forced to either come to a standstill (assuming elastic rheology) or break. Thus, the dynamic forces, which are velocity dependent, would be reduced to zero. This would then require Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 25 that the entire shear component of the plate interaction force be supported along the locked zone. Dislocation models do not allow for any such 'slowing' of the downgoing slab with time. The magnitude of the slip deficit is permitted to grow at the plate convergence rate indefinitely, with no bound placed on the magnitude of the elastic strain energy that can accumulate. Some workers have championed the use of the dislocation model by noting that models based on the same theory have been successfully applied to the problem of predicting the stress drop due to earthquake rupture. According to this argument, the dislocation model must be viable because it predicts a distribution of stress equivalent to that released by earthquake rupture. However, this argument is not self-consistent. Seismic events occur as a result of catastrophic failure along a fault. When this occurs, a large amount of stress must be released in order to reach a force equilibrium along the fault in which shear is no longer supported by the locked zone. The path of least resistance for such a lowering of stored strain energy involves slip along the fault. However, during the interseismic period, there is no physical reason to anticipate that the pattern of strain build-up should take a shape exactly akin to that released during the previous (or subsequent) earthquake. The only requirement is that dynamic balance be maintained at the plate interface. This requires that the total shear supported by static and viscous effects remain constant. In addition, if one employs the dislocation formulation, one permits the possibility that a great earthquake could release all the built-up strain in the system. However, for dynamic balance to be maintained, the shear previously supported by the locked zone must be transferred to other regions of the interface. Hence, regardless of the magnitude of the earthquake slip, a non-zero stress field must persist following an earthquake. Figure 3.6 (following page) Analagous to figure 3.4, plots of the stress field induced by the distribution of slip deficit depicted in figure 3.3b. The three plots, (a), (b) and (c), correspond to the three components of the stress tensor, <rxx + <ryy, o~xx — <ryy and  Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 27 Soon after the introduction of the dislocation formulation by Savage [1983], it was realized that the elastic dislocation model could not explain all the deformation features observed in the geodetic data. Hence, Thatcher and Rundle [1984] proposed two alter- native models based on the same theory. The first model permits the occurrence of the postseismic creep downdip of the locked zone, while the other allows for the presence of viscoelastic rheology in the lithosphere. These models were introduced in an attempt to explain the more rapid deformation which is observed immediately following a great earthquake. Matsu'ura and Sato [1989] modified the dislocation formulation by explic- itly accounting for the long term deformations that accumulates due to persistant slip along the entire plate boundary. The inclusion of gravity in the analysis prevents the long term deformation from reaching large magnitudes. However, all of these models invoke the same basic assumption that subduction dynamics can be modelled using only kinematic constraints, and none have provided a completely successful fit to the available data [Savage, 1995]. No completely satisfactory explanation for the observed problems with the dislo- cation approach has been put forward. One suggestion is that three-dimensional effects may play a role, thus negating the validity of the two-dimensional assumption prevalent in present models. Three-dimensional effects can certainly not be ruled out from con- sideration. Earthquake evidence indicates that seismic ruptures sometimes take place over disconnected portions of a plate boundary fault. Given that geodetic data profiles often include information acquired over a wide segment of the overthrusting plate, non- uniformities in the stress field across this segment could result in deformation signals that are misleading in a two-dimensional context. The alternative explanation for the difficulties encountered with dislocation solu- tions, and the one preferred by this author, lies in the inherently kinematic nature of the dislocation representation. Dislocation models do not reproduce a dynamic picture Chapter 3. KINEMATIC APPROACH: DISLOCATION MODELS 28 consistent with that anticipated based on reasonable physical models of subduction be- haviour. In fact, the exclusive use of the dislocation model for interpretation actually excludes from consideration a number of possible dynamic scenarios. Aside from any drawbacks based on physical arguments, the use of the dislocation model also places restrictions on what can be learned from geodetic data. Since the model is purely kinematic, the sole parameters which are constrained by comparisons of data to model are the size of the locked zone and the slip deficit along this zone. Often this information is obtained from other data sources: analyses of previous great earthquakes, crustal seismic surveys, or thermal models, for example. Comparisons of geodetic observations to dislocation model solutions are then used to help confirm these results. However, the potential exists to extract much more information from geodetic observations. Surface deformations are a product of the complicated dynamic interaction between convergent plates, the nature of which is still not well understood. A simple model of subduction behaviour based on dynamic, rather than kinematic, considerations, would provide much more insight into the nature of the subduction process. It is this question of formulating a dynamic approach to the subduction problem which is taken up in the following chapter. Chapter 4 W E D G E M O D E L 4.1 Introduction Precise measurements of crustal deformation provide a window through which to examine the dynamics of tectonic processes within the subsurface. Stick-slip subduction is one such process. It is of particular interest due to the earthquake risk inherent to stick-slip behaviour. As a consequence, many of the areas with the longest histories of geodetic observation are located near subduction plate margins. However, in the quest to explain the observed deformation, little consideration has been given to the dynamics of the plate interaction. Modelling studies have been primarily focused on development of dislocation models based on a kinematic representation of the stick-slip cycle [Savage, 1983; Thatcher and Rundle, 1984; Matsu'ura and Sato, 1989; Cohen, 1995]. In chapter 3, it was noted that dislocation models do not adequately simulate the physics of the earthquake cycle. Hence, there is legitimate reason to question conclusions based on comparisons of data to output from such models. In this chapter, the modelling problem is reconsidered from the standpoint of subduction dynamics. An alternative representation of the plate interaction which governs stick-slip behaviour is proposed based on dynamic (as opposed to kinematic) constraints. A preliminary mathematical model is then developed based on this representation. The model is employed to explore the nature of deformation within the overthrust wedge; i.e. that section of the overriding plate located above the subducting slab. Particular attention is given to the shape of the strain pattern predicted at the surface of the wedge. 29 Chapter 4- WEDGE MODEL 30 4.1 Conceptual representation The dynamics of subduction play an important role in plate tectonics. Of the eight plate driving forces postulated by Forsyth and Uyeda [1975], four are directly linked to subduction. Hence, the question of subduction dynamics is intimately connected to the nature of the plate driving forces. Previous studies of subduction dynamics have, for the most part, adopted one of two approaches to the problem. In the first instance, an attempt is made to ascertain the relative magnitudes and, therefore, importance of the plate driving forces through parametric investigations. Such studies range from estimations of force magnitudes based on simple thermomechanical models [Davies, 1980; 1983] to finite element studies of entire plates [Stefanick and Jurdy, 1992]. The second approach to the problem is to employ numerical [Whittaker et al., 1992] or experimental [Shimenda, 1993] models to follow the evolution of subduction boundaries. Results are used to help pinpoint the mechanisms behind the less well understood driving forces, and to explain some of the past and present geotectonic features observed at convergent margins. A common feature of many studies of subduction dynamics is the large scale on which the process is modelled. Here the interest is in obtaining a better understanding of the process of stick-slip subduction, a phenomenon which is local to the shallow portion of the plate boundary. The available data for this type of study are deforma- tion observations. Clearly these data, acquired over a small region near a subduction boundary, cannot be expected to constrain the many parameters which must be consid- ered in large scale studies of subduction dynamics. Fortunately the majority of these parameters can be eliminated based on scale considerations. The driving forces of plate convergence act over broad spatial and temporal scales - orders of magnitude larger than those spanned by any set of geodetic observations. Accordingly, these forces can be assumed invariant from the point of view of modelling deformations due to stick-slip Chapter 4- WEDGE MODEL 31 behaviour. This permits formulation of a much simplified representation of the dynam- ics of plate interaction. First, the area of interest is isolated from the overall tectonic picture (figure 4.1a). Zooming in on this zone then permits the grouping of the large scale tectonic forces into a net 'tectonic' force that acts in the far field of the area of interest. This force is balanced by reaction forces acting at the plate boundary (figure 4.1c). Figure 4.1 Causative dynamics of deformation at a subduction boundary, (a) The region of interest located in the tectonic con- text, (b) The overriding and subducting plates converging under the influence of driving forces acting in the far field, (c) Isolation of the two plates as free bodies, with far-field forces balanced by tractions associated with plate interaction. Chapter 4. WEDGE MODEL 32 Finally, it is noted that the observations are located almost exclusively at the surface of the overriding plate. The subducting plate generally sits underneath the ocean. As a result, acquisition of deformation data on this plate is logistically difficult. This state of affairs suggests a model formulation whereby attention is given to the deformation arising in the overthrust wedge. This is achieved by treating the overthrust wedge as a free body, separate from the subducting slab. The influence of the slab is then accomodated by setting appropriate boundary conditions along the base of the wedge. Thus the model space is restricted to the area within the boundaries of the wedge. As discussed in section 2.4, deformation within the wedge then arises due to variations in the distribution of boundary forces along the wedge base. 4.3 Traction at the base of the overthrust For the purpose of formulating a mathematical model, a number of assumptions are made (many of which were introduced in section 2.4). Specifically, the overriding plate is assumed to take the shape of a wedge with opening angle of less than « 20°. Also, a perfectly elastic rheology is assumed. Thus, the problem is reduced to one of two-dimensional elasticity. Prior to calculation, boundary conditions must be chosen. The top surface of the wedge corresponds to the Earth's surface. Hence, free surface conditions can be assumed to hold there. At the base of the wedge, however, the situation is compli- cated by the presence of the subducting slab. The slab is deformable, having elastic properties on the same order of magnitude as those of the overthrust. Therefore, in order to precisely model this boundary, contact interface conditions would have to be incorporated into the model. However, boundary conditions of this type would greatly increase the complexity of the analysis and are not warranted at this nascent stage of model development. Instead solutions are calculated assuming prescribed stresses with no conditions placed on the displacements. These conditions correspond to an isolated Chapter 4- WEDGE MODEL 33 wedge and represents the first of two possible end member approximations. The second of these end members corresponds to a rigid subducting slab for which normal displace- ments along the base of the wedge must vanish. The actual deformation should fall between that calculated for these two end member states. Ideally, it would be useful to obtain solutions for both end member cases. However, the rigid boundary case, with its mixed (displacement and stress) boundary conditions also entails a very involved analysis. Thus, this preliminary study is focused on the characterization of solutions obtained with free surface conditions applied at the base of the wedge. 4.4 The stress field within the elastic wedge The elasticity problem to be solved is that for the deformation of an infinite wedge with stresses prescribed along the boundaries. The problem can be solved by two different transform methods. However, the analysis is fairly involved. Hence, discussion thereof is deferred to Chapter 7. Results of the analysis only are presented here. As was done for the dislocation model described in the previous chapter, the properties of the stress field are examined using a Green's function approach. The fundamental solutions are those obtained for point shear and normal forces applied to the base of the wedge (figures la and c). Solutions for more arbitrary force distributions are obtained by linear superposition of the fundamental results. In general, the net force applied to the wedge results in a non-zero torque about a moment arm which could extend to infinity within the solution domain. Thus, any calculation of the displacements has no physical meaning. For this reason, the solutions presented below axe aimed at characterizing the influence of various force distributions on the strain profile that would be observed at the surface. Both uplift and horizontal strain accumulation are measured by geodetic methods, but here only strain measure- ments may be interpreted because of the infinite nature of the calculated displacements. Chapter 4. WEDGE MODEL 34 0 2 4 x x x Figure 4.2 Point forces applied to the base of the wedge, along with profiles of the horizontal strain along the wedge surface. Plots (a) and (b) illustrate the point shear case, (c) and (d) illustrate the case of a point normal traction, while (e) and (f) illustrate the horizontal force case. Nevertheless, the strain calculations provide valuable insight into the dynamics of the build-up phase of the earthquake cycle. 4-4•! Green's functions The Green's functions are not only used to construct solutions for more general boundary conditions, but they also provide important insight into the dynamics of the Chapter 4- WEDGE MODEL 35 subduction process. The solutions highlight the differing characteristics of the defor- mations associated with shear and normal tractions applied to the base of the wedge. As one might expect, the downwards torque produced by the shear component of the coupling force causes extension at the wedge surface. This is indicated by the positive strain signal shown in figure 4.2b. However, the presence of this extensional regime runs contrary to expectations based on physical intuition. In a collisional environment, one would anticipate that compression would dominate, an expectation that is borne out by the few measurements of strain rate which have been collected in the vicinity of subduction margins [Savage et a/., 1991, Dragert, 1994]. This compression cannot arise as a result of the build-up of shear due to locking along the seismogenic zone. Instead it must arise as a result of accumulating normal tractions acting at the base of the wedge. The Green's function for the normal force case indicates a strain pro- file that is uniformly compressional (figure 4.2d). When the component solutions are superposed to obtain a net horizontal force (figure 4.2e), the calculated deformation is observed to be extensional over a narrow region where the strain signal is dominated by the influence of the shear traction. However, at greater distances from the point of application of the force, the strain is compressional (figure 4.2f). 4-4• % Variation of surface strain with depth of applied force Variation in the strain profile is observed as the depth to the point of application of the force increases. Strain profiles arising due to point forces applied at a number of depths are illustrated in figure 4.3. It is notable that the maximum magnitude of the strain does not vary as the applied force is moved down the bottom boundary of the wedge. However, the analysis does show that this magnitude is directly proportional to the strength of the applied force. Therefore, one has two properties of the deformation signal that may provide a useful 'first cut' method of interpreting a strain profile. The Chapter 4- WEDGE MODEL 36 Figure 4.3 Strain profiles corresponding to point forces applied at different distances from the wedge tip. The profiles are calculated for point (a) shear, (b) normal and (c) horizontal forces. The forces are applied at one (long dashes), three (solid), and five (dotted curve) units from the tip of the wedge. breadth of the strain signal provides an indication of the depth to the force centroid, whereas the signal strength may be used to infer the strength of the coupling force. Chapter 4. WEDGE MODEL 37 Figure 4.4 Plots of the shear and normal tractions acting on the s-axis, for a uniform dislocation. The location of the dislocation on the s-axis is depicted in (a). The shear and normal tractions acting on the s-axis are plotted in (b) and (c) respectively. The associated profile of the surface strain is given in (d). Chapter 4. WEDGE MODEL 38 4-4-3 Dynamic implications of the dislocation As discussed in the previous chapter, dislocation models are often used in efforts to predict deformation above subduction zones. Solutions obtained using the wedge model may aid in characterizing the dynamic implications of these models. Consider the case of a uniform edge dislocation extending from the surface to the downdip end of a locked zone. The tractions, which act on the s-axis (which extends down the plate interface), for this case are calculated using the methods discussed in chapter 3 (figure 4.4a-c). Numerical integration confirms that both the normal and shear tractions have a net magnitude of zero along the plate interface. The resulting profile of surface strain is highly oscillatory in the vicinity of the dislocation. However, in the far field, the strain pattern does predict the compression which is observed landwards of subduction zones (figure 4.4d). Using the wedge model, it is demonstrated that a number of possible force distributions may be proposed that result in the far field compression observed above. For example, the simple horizontal force depicted in figure 4.2e gives a compression signal landward of the point of application of the force. In addition, in figure 4.6a, a force distribution is proposed which approximates the traction profiles calculated using the dislocation theory (in particular the oppositely oriented asymptotic shears shown in figure 4.4b). As anticipated, the resulting strain profile, shown in figure 4.6b, also shows compression in the far field. 4.5 Discussion Above, a new dynamic model of deformation within an overthrust is presented. The model permits any distribution of force to be postulated along the base of the over- thrust, and is therefore more general than the dislocation model. In establishing the dynamic model, a number of important assumptions are made, notably that the bottom boundary is treated as a free surface except where stresses are prescribed. However, despite these assumptions, a number of important issues may be addressed based on the Chapter 4- WEDGE MODEL 39 qualitative characterisitics of the deformation solutions. The most important of these issues is the non-uniqueness which must be accounted for when interpreting observa- tions. Vastly different force distributions at the base of the overthrust can result in very similar strain profiles. For example, a pattern of strain such as that pictured in figure 4.2f could be attributed to a horizontal force or to the combined influence of a pair of shear tractions acting in opposing directions acting along the wedge base (figure 4.6). Results such as these serve as an important reminder that, although models may aid in developing an understanding of the nature of the plate coupling behaviour, a degree of physical intuition will always be necessary when making physical interpretations. X Figure 4.6 The surface strain profile, (b), induced by oppositely oriented point shear tractions at s = 0.8 and s = 1.2, (a). Chapter 4. WEDGE MODEL 40 An interesting question is inspired by the strain profile for a horizontal force applied to the base of the wedge. This profile implies that at significant distances inland of the locked zone the strain is compressional. However, the curve implies that a zone of extension may exist above the locus of the shear force which acts to resist the downward slip of the subducting slab. The question is whether or not the normal tractions due to collision resistance dominate this extension signal. At the time of writing of this thesis, the available geodetic data could not provide a definitive answer to this question. The strain measurements available along the Cascadia boundary all indicate compression, but were obtained well inland of the seismogenic zone [Dragert, 1994]. Chapter 5 C O M P L E X V A R I A B L E A P P R O A C H T O 2D E L A S T I C I T Y 5.1 Introduction The complex variable approach to two-dimensional elasticity was developed in the early part of this century [Muskhelishvili, 1963, Sokolnikoff, 1956]. In this chapter, the approach is introduced, starting from the static equations of equilibrium in an elastic, isotropic body. First, the equations are simplified given the assumptions of plane strain. The Airy stress function is then introduced, and the problem is reduced to one of solving the biharmonic equation. Finally, the problem is transformed into a simpler form by writing the Airy stress function in terms of two analytic functions of a single complex variable. The form of the analytic functions is discussed for the general case of a multiply connected domain. Finally, the utility of the method is demonstrated through the analysis of the simple, yet instructive, problem of a point force applied to the whole plane. 5.2 Plane strain The basic equations for the static equilibrium of an elastic, isotropic body consist of three equations: the equilibrium relation which relates the stress components to any applied forcing function, the stress-strain relationship, or Hooke's law, and the strain- displacement relation. Here the complete system is given [see, for example, Fung, 1965] <Tij,j + Xi = 0, (5.1) (5.2) and 41 Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 42 where i, j G {1,2,3} and summation over repeated indices is assumed (see the Standard Notation of Elasticity for definitions of other standard notation). This system can be significantly reduced for problems cast in two dimensions. Here, the case of plane strain is considered. The domain of solution is assumed to extend from —oo to +oo in the third coordinate direction, and no external forces are applied in that direction. Thus, all derivatives with respect to that coordinate vanish. Assuming the direction in question corresponds to the 3-axis of a Cartesian system, these conditions can be written mathematically as, (») u 3 = 0 (ii) d3q = 0 (iii) <Ti3 = C23 = 0 where q is any physical parameter characterizing the domain of solution. From the first two of the above conditions and (5.3) it is easily deduced that ezz = 0. However, the normal traction, perpendicular to the 12 plane, is non-zero. From (5.2) it is given by 0"33 = V(Taa. (5.4) Given the above simplifications, the equations of elastic equilibrium reduce to <rap,a + Xa = 0, (5.6) Co/3 = \(U<*,P + U 0 ,c«)> ( 5 - 7 ) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 43 where a, 3 £ {1,2}. Thus, under the assumptions of plane strain, the system to be solved is reduced from a set of fifteen equations to a set of just eight. 5.3 Airy stress function The system of equations for the plane strain case can be further reduced by elim- inating the displacements and strains from (5.5-7). First, consider the expression for the off-diagonal component of the strain tensor in terms of the two displacement com- ponents, «12 = ^ ( u i ) 2 - r u 2 ) i ) . (5.8) Taking the derivative of the above with respect to the x\ and Xi coordinates, this becomes C12.12 = ^( u l ,122 + U 2 ,112), (5.9) or, from (5.7), 2ei2,12 = Cll,22 + ^22,11- (5.10) This last expression is often referred to as the compatibility equation, and can also be derived from the condition that the strains correspond to a single-valued displace- ment field. Upon substitution of the stress-strain relation, (5.5), the equation expands to 0"l2,12 + 021,21 = <Tn(22 + 0-22,11 ~ ^(^11,22 + 022,22 + 0"ll , l l + 0-22.ll)- (5-H) Then, from the equilibrium relation, (5.6), it is found that 0-11,1 + 0-21,2 + X\ = 0 <72i,2i = — 0-11,11 — -X"i,i, (5.12) 0-12,1 + 022,2 + %2 — 0 012,12 = — 022,22 — %2,2- (5.13) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 44 Hence (5.11) can be written as (1 - v)(o-iltu + erU)22 + o-22,ii + 0-22,22) + Xiti + X2,2 = 0, (5.14) or (1 - vWcjift + Xa,a = 0. (5.15) 'Now assume that the body force can be expressed as the gradient of a harmonic po- tential. Then, Xa,a = <P,aa = 0, where <p is the potential. Given this assumption, the system of equations, (5.5-7), reduces to a system of three equations for the stress components, Vaotfip = 0, (5.16) <rap,a + X/3 = 0. (5-17) The equations, (5.17), are satisfied by writing the stress components in terms of a scalar function A(x,y) as d2A d2A d2A . l o . ^=QyT + ^ " w = w + t p a x y = -dx-dy- ( 5 - 1 8 ) or eP&^M 4- <p8ap (5.19) where e a0 7 is the standard permutation tensor. Thus the only equation that remains to be satisified is (5.16). Upon introducing the expressions, (5.18), and noting that <ptXX +<P,yy = 0, (5.16) reduces to A,1111 + 2A,ii22 + A ,2222 = 0 or V 4 A = 0, (5.20) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 45 which is the biharmonic equation. The function, A, is termed the Airy stress function. 5.4 General solution to the biharmonic equation At the turn of the century theorists derived general solutions to the biharmonic equation in terms of the complex variable, z — x + iy, and demonstrated the utility of these solutions in solving particular problems in two-dimensional elasticity [see, for example, Muskhelishvili, 1963; Sokolnikoff, 1956; see also Fung, 1965, p.250-9]. The derivation of the general solution to the biharmonic equation makes use of a function, g(z), which is analytic throughout the domain of solution, S. By definition, the deriva- tive of such a function exists and is unique throughout S [Churchill et al., 1974, p.46]. Thus, letting u and v be the real and imaginary parts of g, respectively, the derivative can be written as do du .dv .du dv /-«-.\ -f = — +%— = -i— + — . (5.21 dz ox ox dy dy Hence, du _ dv du _ dv . . dx dy dy dx' Upon differentiating (5.22) with respect to x and then with respect to y, it can be readily shown that u and v are harmonic, i.e. that V2u = 0 and V 2 v = 0. Now consider a function, / , which satisfies the biharmonic equation, V 4 / = V2(V2/) = o. It follows that V 2 / is harmonic. Letting P — V 2 / , by (5.22) a conjugate harmonic function Q can be determined such that P + iQ is analytic. The integral of P + iQ is also analytic. Let (f>(z) = j f (P + iQ)dz = R + il (5.23) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 46 where R(z) and I(z) are the real and imaginary parts of <p{z). Now consider Therefore, if / is defined by f = xR + yl + q (5.24) where q is an arbitrary harmonic function, then it follows that / is harmonic. Based on this result, / can be written as where the overbar denotes the complex conjugate, IRQ designates the real part, and %(z) is an analytic function with real part q. The result, (5.25), was originally obtained by Goursat. However, the derivation given here follows that of Muskhelishvili, which is slightly more general than the original [Muskhelishvili, 1963, p.110-1]. The Airy stress function, A, which satisfies the biharmonic equation, can be sub- stituted directly for the function, / , in the above solution. Hence, / = *[Z#z) + X(z)] (5.25) A = + x(z)) (5.26) or, upon expanding the right hand side, 2A = z<f>(z) + z4>{z) + X(z) + X (z) . (5.27) Thus the solution of the Airy stress function reduces to the problem of determin- ing two analytic functions, 4>(z) and x(z), which satisfy conditions placed along the boundary of the elastic region. Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 47 5.5 Stresses and displacements in terms of the analytic functions Expressions for the stresses in terms of the analytic functions can be obtained by direct substitution of (5.27) into (5.18). The partial derivatives are easily evaluated upon noting that d_ dx and d_ dy dz d dz d dx dz dx dz d_ d_ dz dz (5.28) dz d dz d dy dz dy dz _ = .(d__d_\ dz \dz dz) (5.29) A lengthy, but straightforward, calculation yields the following useful formulae for the stress components: <rxx +(Tyy=2 [<j>'(z) + ~&{zy) (5.30) o-yy - axx + 2iaxy = 2 [z<f>"(z) + $'(zj\, (5.31) where, as is customary, the function %'(z) has been replaced by rj>(z). In addition, the Laplacian of the Airy stress function can be expressed in terms of <f>(z) using (5.18) and (5.30). Thence, V 2 A = 0~xx + (Tyy = m{(t>\z)}. (5.32) where 5ft[] denotes the real part. The displacements can also be specified in terms of <f>(z) and i>(z). However, the derivation is a little more arduous. First the strains are written in terms of the Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 48 Airy stress function. The displacements are then calculated by integrating the strain- displacement relations. Expressions for the strains in terms of the Airy stress function are obtained by substituting (5.18) into (5.5), which yields £ , , = i [ ( l - , ) V * A - ^ ] , (5.33) «n = ^[( l-")^A-vg], (5.34) i B2K 2/x dxdy (5.35) The z-component of the displacement, u, is given by du dx Substitution of (5.32) and (5.33) into (5.36) yields = i - [ ( l - , W ' ( z ) ] - 0 ] , du dx 2/x (5.36) (5.37) which, upon integration, gives an expression for the the displacement, u, u 2/x ( i - *)4St[<K*)] - + f(y) (5.38) where / is an arbitrary function of y. The equation, (5.34) can be subjected to a similar analysis to obtain an expression for the y-component of displacement, v. The result is v = j - [(1 - i/)49Mz)] - ^] + g(x) (5.39) 2/xLv" dy. where denotes imaginary part and g is an arbitrary function of x. The final stress-strain relation, (5.35), limits the arbitrariness of / and g. Substi- tuting (5.38) and (5.39) into this relation and simplifying, yields the result, f'(y)-+g'(*) = o. Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 49 Therefore f(y) = ay + 7 and g(x) = —ax + 7' where a, 7 and 7' are arbitrary constants. The forms of f(y) and <?(a:) indicate that these functions represent rigid body displacement and rotation. Thus, they can be ignored in deformation calculations. The displacements can be expressed in a useful, complex form by summing expres- sions for u and iv obtained using (5.38) and (5.39). The result is 2ti(u + iv) = A(l-u)<f>(z)-(^+i^). Finally, substituting for A from (5.27), and recalling (5.28) and (5.29), the displacement expression simplifies to 2u.(u + iv) = K<f>(z) - z<p'(z)-ip(z), (5.40) where K — 3 — Av. 5.6 Force on a contour The force on an arbitrary contour is now evaluated in terms of the complex po- tentials <f> and ip. Let the contour run from A to B such that it does not cross the boundary of an elastic domain, S (see figure 5.1). A unit normal vector, n, is defined such that it points to the right when travelling in the direction of the tangent vector, t. The net force, F, on the contour is obtained by summing over the surface tractions which act along its length. Hence, (5.41) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 50 Figure 5.1 Contour running from point A to point B within an elastic domain, S. Introducing, from (5.19), the Airy stress function, this becomes, Fa = / e Q 73 epfr A yl£npds. JA (5.42) However, is equivalent to a 90° rotation tensor in two dimensions. Hence e ^ n ^ = te, and Fa = J eaj3An(t(ds, (5.43) but t{ = dx^/ds. Thus, Fa = J e ^ A ^ ^ d s (5.44) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 51 Given (5.44), it is a simple exercise in differentiation to write the components of F in terms of the analytic functions, (f> and ij>. In its most useful form, the end result of the calculation is, Fi + iF2 = -i [<f>(z) + zfiiz) + ^ ) ] BA (5.45) 5.7 The general problem The general form of the analytic functions, (j> and must allow for the possibility that the elastic body is multiply connected. Let the region occupied by the body be S, which is bounded by m + 1 contours, Lk, the outermost of which, I m +i , may or may not extend to infinite radius. In illustrative terms, such a domain can be described as a plate with holes (figure 5.2). By definition, <f> and ij; are analytic throughout S. A necessary consequence of this fact is that these functions must be single-valued on any simply connected part of S. However, on the whole of S, or on any multiply connected subregion thereof, <f> and ij> may be multi-valued. That is, upon describing a contour within S encircling any multiply connected part of S (such as L* in figure 5.2) the value of an analytic function may increase by some increment. In the following, it will be necessary to differentiate between single and multi-valued analytic functions. Hence, the following definition is introduced. Definition 5.1 A function, / , is termed holomorphic over a region, S, of the complex plane if / is single-valued over the entirety of S. \ The multi-valuedness of <j> and ij) is constrained by the requirement that the stresses be single-valued everywhere within S. The equation, (5.53), o-xx + o-yy = 4.U[<f>'(z)] t Definition 5.1 follows the usage of Muskhelishvili (or, at least, that of his translator) [1968, p.120-1]. On simply connected regions analytic functions are necessarily holomorphic. Hence, the distinction introduced in definition 5.1 is not employed in many standard texts on complex analysis. In fact, the terms analytic and holomorphic are often defined synonymously [in particular, see Marsden, 1973, p.44; Carrier et. al., 1983, p.26]. Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 52 shows that the real part of <j>'(z) must be single-valued. However, as no restriction is placed on the imaginary part, 4>'(z) may increase by some increment, say 2mAk, where Ak is a real constant, upon travelling a closed contour about one of the Lk- One multi-valued function which exhibits this property is log z. Consequently, <f>'(z) can be written as m 4>\z) = $0(z) + Y, A* l oS(z - ( 5- 4 6) fc=l where $o(z) is an analytic and holomorphic function in 5, and the z* are points located outside S. Upon integrating (5.46), one finds that $o(z)dz + ] T ] Au[(z - zfe)log(z - zfc) - (z - Zk)]'+ const. (5.47) -o Jfe=l where zo is some point in S. Since <f>(z) is analytic, so must be the first term on the right hand side of (5.47). However, unlike its derivative, $o(z), this term is not necessarily single-valued in S. Hence, the most general expression for 4>(z) is m m 4>(z) = z^Aklog(z - zfc) + ^ fklog(z - zk) + <f>o(z) (5.48) k=l k=l where <j>o(z) is holomorphic in S, and the 7*. are complex constants. A general expression for i{>(z) can be obtained by considering the other stress relation, (5.31), which, along with the requirement that the stresses be single-valued, necessitates that tp'(z) be holomorphic. Thus, ij>(z) can be written as TO Tj>(z) = Y,T'kl<>s(z-Zk) + M') (5-49) Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 53 i ' x 2 • JHT L m + i ^^^^^BH^^^^^^^^B s ^ ^ 2 ^ ~ \ t * *1 Figure 5.2 A contour, X*, encircling a multi-connected part of the domain, S. where i>o(z) is holomorphic in S, and the are complex constants. The equations (5.48) and (5.49) combine to form a complete general solution which can be applied to any problem in two-dimensional elasticity. 5.8 Application of the method The general solutions for <f> and ip are particularly useful when the problem under consideration includes singularities either on the boundary or embedded within the solution domain. Often, particular solutions to such problems can be simply obtained by colocating the logarithmic singularities with the singular features of the problem. Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 54 As an illustration of the physical significance of the singular, or logarithmic, terms in the general expressions for <j> and ip, consideration is given to the problem of an infinite plane acted on by a point force. Assume the force is acting at the origin. This point must first be excluded from the solution domain. Hence, a contour is defined of infinitesimal radius about the origin. Then, as the origin is the only point excluded from the elastic domain, the general solutions (5.48) and (5.49) must simplify to <f>(z) = zAlogz +7logz + d>0(z), (5.50) >(z) =7'logz + ^o(«). (5-51) Now, consider an arbitrary closed contour which runs in the anti-clockwise direction about the origin. The net force on such a contour can be calculated using the result, (5.45). Substitution of (5.50) and (5.51) into this expression yields an expression for the net force in terms of the arbitrary constants in the general solutions. Hence, Fi +iF2 = - 2 7 r [ ( z - r l ) A + 7fc -7*] (5.52) where the Fi are the components of the applied force. The net change in displacement upon completing a circuit of the contour is calculated using (5.40). Upon substituting the expressions, (5.50) and (5.51), one finds that 2u.[u + iv]L = 2m [(« + l)zA + KJ + K-y']. (5.53) Here, the displacements are required to be single valued. Thus, u + iv = 0, and the left hand side of (5.53) has no z dependence. Therefore, A can also be set to zero. This reduces the problem to a set of two equations for the unknowns, 7 and 7'. Not only that, but, once (5.52) and (5.53) have been solved, all the boundary conditions are satisfied. Hence, the functions <f>o(z) and ipo(z) are redundant and can be set equal Chapter 5. COMPLEX VARIABLE APPROACH TO 2D ELASTICITY 55 to zero. The final solution, which is obtained upon solving (5.52) and (5.53) for 7 and 7' is ™ = ( 5 - 5 4 J ) As a final note, it is interesting that the displacement undergone after a circuit of a closed contour need not necessarily be single-valued. The expression, (5.53) allows for the possibility of a multi-valued displacement field. As a result, the complex variable method also has important application to problems of dislocation theory. One such problem is discussed in the following chapter. Chapter 6 D I S L O C A T I O N I N T H E E L A S T I C H A L F P L A N E 6.1 Dislocations in continuous media The theory of dislocations in continuous media has played an important role in geophysics since being introduced to the field by Steketee [1958b; also see Rybicki, 1986 for a review]. Particular attention has been given to the problem of determining the deformation which arises due to dislocation sources located within an elastic half-space. For example, a comprehensive set of analytic solutions has been obtained for the case of an arbitrary point or finite rectangular dislocation located within an elastic half-space [Okada, 1985; 1992]. In the context of modelling the deformation above a subduction zone, the dislocation is assumed to be of infinite length along the strike of the fault, with a slip vector that points downdip. Thus, the problem of interest is that of an infinite edge dislocation with slip vector located in the plane perpendicular to the half-space boundary. This problem is somewhat different from those solved by Okada and his predecessors [e.g. Steketee, 1958a; Chinnery, 1963; Press, 1965] in that the dislocation is of infinite length. Hence, the deformation is invariant along the direction of strike. This has the advantage of permitting the problem to be formulated in two dimensions, which greatly simplifies the analysis. The complex methods described in the previous chapter may then be employed to obtain a solution. In this chapter, the problem of an edge dislocation within an elastic half-plane is solved in terms of the complex variable. The method of solution follows that described by Freund and Barnett [1976]. 56 Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 57 A*2 Figure 6.1 An elastic domain with a branch cut along the seg- ment AB. 6.2 Dislocations and branch cuts A dislocation is a discontinuity in displacement. In order for a dislocation to persist in an elastic domain, a non-zero stress field must be maintained. This stress field can be determined by a number of methods. The complex variable approach is particularly well suited when the problem can be cast in two-dimensions. Consider, for example, a two-dimensional domain, 5, with a dislocation defined along the line, AB (figure 6.1). Let the dislocation run from a point, zo, to a point on the boundary of S. Clearly, the stresses will be singular at ZQ. Hence, ZQ is excluded from the solution domain. This is done by introducing an infinitesimal contour about this point. If a complete circuit is made around this new contour from a point immediately above the cut to a point immediately below, then, based on (5.40), (5.48) and (5.49), the observed change in displacements is given by Au + iAv = — [(K + l)AiZ + «7i - 7̂ ] Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 58 where, in deriving this expression, it was noted that the only possible singular point located within the travelled contour is located at zo- A branch cut must be included in 5, coincident with the line of dislocation, in order to preclude the possibility of multi-valued displacements occurring anywhere else in the solution domain. 6.3 Edge dislocation in the infinite plane Here the problem of interest is that of an edge dislocation in an elastic half-space. However, before proceeding with the solution to that problem, it is useful to consider the problem of an edge dislocation within a whole space. The geometry of the problem is depicted in figure 6.2a. A simply connected domain is defined with an inner contour about the one point, zo, A branch cut is defined such that it follows a path coincident with the dislocation. The magnitude of the displacement discontinuity is prescribed as [uo + ivo] = beia, where a is the polar angle between the z-axis and the cut in the solution domain. Let 4>i(z) and i>i(z) be the complex potentials that determine the stress and dis- placement fields within the plane. A coordinate transformation can be performed which shifts the end-point of the dislocation from z = ZQ to the origin. The transformation is defined by The complex potentials, <j>(u>) and ij>(u}), m the u; plane are related to rpi(z) and tpi(z) by [Muskhelishvili, 1968, p.135] U) = z — z 0. (6.1) <£i(z) = <f>(u) = <f>(z - z0) (6.2) V>i(z) = ij>(u) - z04>\u}) = V>(z - z„) - z0<£'(z - z0) (6.3) The general form for the complex potentials which ensures single-valued stresses can be introduced directly based on the results of Chapter 5. Based on (5.46) and Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 59 (5.47), and upon noting that the point at u> = 0 is not in the solution domain, the potentials are seen to be given by 4>{w) = u>Alog(u>) + 7log(u;) + <f>*(u) V>(u>) = 7' log(u>) + ij>*(u), (6.4) (6.5) where <f>*(u) and •0*(w) are functions holomorphic in the u; plane, Ak is a real constant, and 7k and 7k are complex constants. (a) t % 1 k (b) branch cut z - plane co - plane Figure 6.2 The dislocation located in, (a), the z-plane and, (b), the u>-plane. Note the branch cut oriented along the direction of the Burger's vector. The constants in equation (6.4) and (6.5) are a function of the boundary con- straints. As no net force is applied along the boundary of the domain, it follows from (5.45) that 0 = i[<f>{u) + u<j>'(u) + (̂w)Uk (6-6) Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 60 where Lk is a contour running from a point immediately above the dislocation to a point immediately below. Substitution of the potentials, (6.5a-b), into (6.6) yields 0 = -27r[(w + l )A+7-r7 T ]. (6.7) A second constraint on the constants in equations (6.4) and (6.5) is obtained upon integrating the displacement equation, (5.40), along a similar contour. The result is the expression, 2/i[u + i v ] L k = [k#w) - w^'(w) - V>MJLk (6-8) where the left hand side can be equated to the dislocation magnitude. Substitution of this magnitude, along with (6.4) and (6.5), into (6.8) gives 2uheia = 27ri(kAu> + K~f + wA + T 7 ) . (6.9) In neither (6.7), nor (6.9), is there any dependance on the coordinate, u>, on the left hand side of the equations (i.e. no discontinuity in the rotation is prescribed in the problem). Thus, A = 0, and (6.8) and (6.9) simplify to 0 = 7 - 7 , (6.10a) 2ubeia = 2iri[Kj + T 7 ] . (6.106) This set of equations is easily solved for the two constants 7 and 7'. Hence, — ubeia 7 = 7* = S R i ^ ( 6 ' U ) where it has been recalled that K = 3 — 4i/. The complex potentials are obtained immediately upon substitution of (6.11) into (6.5). In addition, all the boundary conditions are satisfied upon determining the Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 61 constants 7 and 7'. Thus the functions <f>*(u) and i>*(w) are redundant and can be set identically to zero. Thus, 4>(w) = -iBlogw, (6.12a) V>(u;) = iBlogu;, (6.126) where B = ^ 4TT(1 - v) The above functions can be transformed to the original z coordinate system using the relations (6.2) and (6.3). Thus, <£i(z) = <j)(z - z„) = -zBlog(z - z c), (6.13a) and zBz - V»i(z) = >̂(z — z0) + z^<f>'(z — z0) = iBlog(z - z0) H —. (6.136) Z ZQ For the purpose of calculating the stresses, it is often more convenient to work with the first derivatives of the complex potentials. Upon performing the required operation, one finds that —iB = ( 6 ' 1 4 a ) z — za , 1 / tB zz0B ., , , . = ̂  - ( i ^ F - ( 6 - 1 4 J ) Based on these results, the stress and displacement fields can be calculated using the appropriate relations derived in Chapter 5. Hence, the problem is solved. 6.4 Point dislocation in the elastic half-plane The solution to the problem of an edge dislocation in a half plane is readily obtained by adapting the solution to the corresponding problem in the whole plane. At the Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 62 surface of the half-plane, free surface boundary conditions must be satisfied. Hence, by (5.31) and (5.32), 0 = T y y + irxy = 4>'{z) + fiiz) + z(j>"(z) + i>\z), (6.14) for z located along the half-plane boundary. Solutions may be sought in the form 4>{z) = <f>i{z) 4- <£0(z), (6.15a) i>{z) =V1(z) + V'o(z). (6.156) where <j>\ and are the potentials, calculated in the previous section, for the edge dislocation in a whole plane, and <f>0 and if>0 are functions holomorphic in the half plane which are necessitated by the presence of the traction free boundary. Without loss of generality, this boundary can be located at y = 0, and the domain of solution taken to be the upper half plane. Then, upon substituting (6.15a-b) into (6.14), and evaluating for points along the free surface, one finds that 0 = <f>[(x) + ^) + x4>'[{x) + +[ (x) + 4>'0{x) + #^) + xftix) + 1>'0(x). (6.16) By definition, d>0(z) and tp0(z) must be holomorphic in the upper half-plane. How- ever, these functions may possess singularities which fall in the lower half-plane, since such singularities will not affect the holomorphic character of the functions within the solution domain. Upon noting that, along the domain boundary, z = z = x, it follows that the condition, (6.14), can be satisfied by decomposing (6.16) as follows, tiXz) + 0i(z) + z<j>'{{z) + ^i(z) = 0, (6.17a) <f>'o(z)+z^(z) + ^o(z)+^) = 0. (6.176) This decomposition ensures that <j>0 and ij>0 are functions of z (as opposed to z), and have no singularities in the half plane comprising the domain of solution [Freund and Chapter 6. DISLOCATION IN THE ELASTIC HALF PLANE 63 Barnett, 1976]. In fact, as will be observed below, these functions are only singular at z = z~o, which is a point in the lower half-plane. Using the expressions for <f>(z) and ip(z), obtained in the previous section, the equations (6.17) are readily solved for <f>0 and if>0. Finally, upon substituting the results of this calculation into (6.15), and taking the derivative with respect to z, one finds that AII \ *B , *B , *B(z7-z 0) . . 4>{z) = + = + , _ s 2 , (6.18a) z — zQ z — z„ (z - zoy iB iz0B iB (z - z0) (z - z0)2 (z - z0) Given (6.18a-b) the problem is solved. Physical quantities such as the stresses, strains and displacements are obtained by substitution of these potential solutions into the formulae derived in chapter 5. However, such calculations are beyond the scope of this chapter. A number of examples are discussed in chapter 3 in the context of modelling deformation at stick-slip subduction zones. Chapter 7 T H E F I R S T F U N D A M E N T A L P R O B L E M F O R T H E W E D G E 7.1 Introduction The problem of determining the deformation within an infinite elastic wedge (figure 7.1) subject to traction boundary conditions has been studied for over a century. The first analysis of the problem is attributed to Levy who, at the end of the last century, wrote a solution for a wedge loaded by uniform normal pressure on one face [Levy, 1898]. Carothers [1912] later provided a solution to the problem of a concentrated couple applied to the vertex of the wedge. These early results were obtained by assuming a particular form of the Airy stress function, which is separable in polar coordinates. However, the assumed form of the radial dependance of the stress function leads to solutions which break down for wedge angles greater then 180° [Sternberg and Koiter, 1958, Dempsey, 1981]. In particular, for loading, antisymmetric with respect to the wedge bisector, the classic solution breaks down at a critical angle, 2a « 257°. With symmetric loading, the solution breaks down when 2a > 180°. This break down of the classic solution at large wedge angles is known as the wedge paradox. The wedge problem is also effectively tackled using transform methods. The form of the biharmonic equation in polar coordinates lends itself to analysis based on the Mellin transform [Brahtz, 1933; Shepherd, 1935; Tranter, 1948]. Alternatively, one may begin with the biharmonic solution expressed in terms of analytic functions of a complex variable. In this case, the problem lends itself to an analysis based on the Fourier transform [Green, 1968, p.304-9]. Much of the work on the wedge problem which has appeared since the early studies cited above, has been focused on resolving the wedge paradox. In this context, a 64 Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 65 concerted effort has been made to determine the order of the singularity as r —> 0, and with some success. Dempsey [1981] found admissable solutions for large wedge angles by extending the class of Airy stress functions which are assumed in the traditional separation of variables approach to the biharmonic equation. [Sternberg and Koiter, 1958] found leading order expressions for the stress field arising due to an antisymmetric loading using the Mellin transform approach. Bogy [1968] has employed the same method to examine the similar problem of determining the stress field within two dissimilar wedges bonded along an edge. Here, the wedge problem is considered in the context of the geophysical appli- cation described in Chapter 4. The objective is to calculate the deformation pattern that would be observed at the surface of a continental overthrust wedge due to some distribution of tractions along its base. Hence, a complete solution for general loading is required. In the first part of the chapter, the elasticity problem is solved by both the complex Fourier and Mellin transform approaches. Both approaches lead to solutions which are expressed as integrals with respect to the transform variables. In the latter part of the chapter, the equivalence of the solutions obtained by the two methods is demonstrated. Techniques of evaluating the inverse transform integrals are discussed. A final paragraph is included in which the expression of the wedge paradox in the transform domain is described. 7.2 Statement of the problem The objective is to determine and characterize the stress and displacement distri- butions within an infinite two-dimensional wedge arising as a result of some arbitrary forcing along a finite stretch of the boundary (figure 7.1). The wedge is assumed ho- mogeneous, elastic and isotropic. The equations to be solved are those of plane strain, Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 66 Figure 7.1 Solution domain and coordinates for the wedge prob- lem. which reduce to a single biharmonic equation for the Airy stress function (q.v. section 5.3), V 4 A = 0. (7.1) The boundary conditions for the problem can be expressed as <h{r) e = a (9l(r) 9 = a aee — \ °~rQ = \ • (7-2) (f2(r) 6 =-a (g2(r) 9 = -a Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 67 where f\, f2, g\ and g2 are prescribed functions. Finally, the solution must satisfy a regularity requirement that the integrals, I (Tee dr and / crr0 dr (-a < 0 < a), (7.3) Jo Jo be convergent. This requirement ensures that the net force acting along any radial line remains finite. 7.3 Complex variable approach In chapter 5, a general solution to the biharmonic equation was obtained in terms of two analytic functions, <f> and ip, of a complex variable, and the components of the stress tensor expressed in terms of these two functions. Transformed to polar coordinates the resulting expressions, (5.30) and (5.31), become + 0-00 = 4[<P'(z) + V(z)\, (7.4) arr - O-00 + 2i<rr0 = -A[z<P"(z) + -ip'{z)\. (7.5) z Subtraction of (7.4) from (7.5) results in the useful relation CT08 - io-ro = 2 [<p'(z) + #{z) + z<P"(z) + - ^ J ] , (7.6) z which, along with its complex conjugate and the boundary conditions, (7.2), is sufficient to allow determination of the stress field. Introducing the conditions that the stresses vanish at infinity, that there be no rotation at infinity and that the displacements be single valued, the form of of the solutions for <p(z) and ip(z) can be restricted to those functions of the form (q.v. 5.48 and 5.49), Tl <P(z) = A^log(* - zk) + p{z) (7.7a) k=l n tf(z) = £ £ l o g ( * - z f c ) + ^ ( * ) (7.76) fc=i Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 68 where <f>*(z) and V,*(z) a r e functions holomorphic in the region under consideration, and the zj. must correspond to points located outside this region. 7.3.1 Solution by Fourier transform One method of solving for <f>(z) and i>(z) is to conformally map the wedge into an infinite strip, to which Fourier methods, generalized to the complex plane, are applicable. The solution methodology was first presented by Green [1968, p.304-9]. The Fourier theorem, generalized to complex functions, states that, if we let /(z) be a function of the complex variable, z = x + iy, holomorphic in the region a < y < b where a and b are real constants, and defined such that is itself holomorphic in the region, A < y* < //, where s = x* + iy*. Also, within any infinite strip interior to this region X — • C O X — > — C O (7.8) for e > 0, then the function x* — > C O C O (7.9) and c+ioo (7.10) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 69 i C, p l a n e z p l a n e Figure 7.2 Conformal transformation of the infinite wedge to the infinite strip. where c is real and fi < c < A [Green, 1968, p.50]. 7.8.2 Mapping and transformation of the boundary conditions The conformal transformation used in the Fourier method maps the wedge geom- etry onto an infinite strip. This transformation is given by z = o>(0 = ae< (7.11) where £ — £ + it]. In the z plane, the axes are defined such that the x-axis bisects the wedge and the origin is colocated with its corner. Hence, (7.10) maps the boundary of the wedge in the z plane; (i.e. the lines 6 = ±a), to the lines rj = ±a in the ( plane. The constant, a, defines the distance from the vertex of the wedge to the two boundary points which are mapped onto the rj axis of the £ plane (figure 7.2). Multiplication of (7.5) by ae^+t^s and integration with respect to £ from —oo to +co, yields fOO {aee-i<Tre)e^8di = / oo »oo »oo — a<j>'(z)et+it'dt+ / a{(j>'(z) + z<f>"(z))e(+i(8d£+ / a - ^ ' ( z ) e * + * a d £ , (7-12) -oo J—oo J—oo Z Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 70 which suggests the introduction of the transform functions, *(«) = f <j>{zy<*di, (7.i3) J—oo -oo coo V(s) = f i){z)eiisdi (7.14) J—oo where the above functions are defined over the strip r) = ± a , in the C plane. The asymptotic behaviour of <j> and ip is determined from (7.4a-b). As \z\ —• oo, or as | oo, <P(z) == 0(log z) = 0(C) = O(log z) = 0(C) (7.15) while as |z| —> 0, or as £ —> — oo, <P(z) = 0(1) = 0 ( e ° ) TP(Z) = 0(1) = 0 ( e ° ) . (7.16) Hence, by (7.8), the integrals, (7.13) and (7.14) converge as long as 0 < < 1. (7.17) The expressions for the inverse formulae follow directly from (7.10), (•c+ioo -I /.c-l-too <K*) = 5- / (7.18) 2 7 r . /c- ioo V»(z) = r - / V(s)e-*'ds. (7.19) By integrating by parts the terms on the right hand side of (7.10) and introducing the transform relations (7.13) and (7.14), one may show that, - - / {(TOO - icrrg)e^3d( = isQWe*'-** + « a $(-«)e-" ( ' - i ) + wtf (*)e-"('-i). 2 J-OO (7.20) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 71 The above equation, upon setting rj = a and —a, yields the following expressions s$(s) - is2$(-s)e-2a(s-V + sV(s)e-2aa = £ [F^s) - iG^s)] e"̂ '-*) (7.21) s$(s) - is2^{-s)e2a^ + sV{s)e2at = ̂  [F2(s) - iG2(s)} e*-* (7.22) where the functions, Fy and G 7 , 7 € {1>2}, are functions, denned below, which depend on the boundary conditions. Noting, based on (7.11), that r = ae*, it follows from (7.2) that o~eo Therefore ( fi(aet) rj = a ( g^ae*) rj — a < <Tr0 = < . (7.23) [h(aet) TJ = -a { 92(aei) 77 = -a Fy(s) = f°° fy(aet)eX1+i*Ut and G7(s) = f°° 0.y(ae«)e«1+">d£. (7.24) J—00 J—oo 7.3.3 Symmetric normal stresses with zero shear on the boundary Equations (7.22) and (7.23) along with their complex conjugates are sufficient to permit determination of $(«) and \P(s). For the four symmetry cases listed in Table o'ee o~ro Case 1 fi(r) = • /a(r) = /(r) 9i(r) = 92(r) = 0 Case 2 fi(r) = -Mr) = /(r) 9i(r) = 92{r) = 0 Case 3 Mr) = / 2(r) = 0 9i(r) = : 92(r) = g(r) Case 4 fi(r) = /a(r) = 0 9i(?) = ~92(r) = g(r) Table 7.1: Boundary tractions for the four symmetry cases. Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 72 7.1, these functions can be obtained by inspection. Here the analysis is continued for one of these cases, that of symmetric normal stress along the boundary, with zero shear. Two solutions are given for this case, that for a general force distribution, along with that for the particular case of a point, or Dirac-delta, forcing function. The latter solution represents the Green's function for the problem. The analysis follows a number of steps. First the equations (7.21) and (7.22) are solved for the given boundary conditions. These expressions are then substituted into the stress relations, (7.4) and (7.5), thus reducing the problem of calculating the stress components to one of quadrature. Finally, the integral expressions are simplified for the particular case of a Dirac-delta distribution of traction along the two boundaries. Solutions for the other cases may be obtained by a similar procedure. The analytic functions in the Fourier domain For a symmetrical system with vanishing shear stress along the boundary, F2(s) = F^s) = F(s), G2(s) = G^s) = 0. (7.25) Substitution of the above into (7.22) and (7.23) and elimination of ^(s), yields the expression, 4s($(s) sinh 2aa + s${-s) sin 2a) = ia(eia+a" - e-ia-a*)f(s). (7.26) Taking the complex conjugate of the above and replacing s with —s, one obtains 4s($(-s) sinh 2sa + s$(s) sin 2a) = ia(eia+as - e-ia-as)f(s) (7.27) from which it follows immediately that, \ s) ~ \ s ) - ^(gink 2sa + s s i n 2a) \ • ) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 73 The expression for ̂ (s) is obtained by back substitution of (7.28) into (7.22): v ' 4s(sinh 2sa + s sin 2a) ' v ' ' Integral expressions for the stress components Employing the inverse Fourier transform expressions, (7.18) and (7.19), the stress components can be written in terms of integrals over the transform variables. Hence, F O O + 1 C r( • • • - - ' — i -oo+tc r(crrr — (Tee + 2iar0) = (7.31) + C (5(1 + w)$(«)¥(-*)e<- ,"+* ,'> - s^(s)e^-,",-i^)e-ii8ds -oo+tc Substitution of (7.28) and (7.29) into the above yields expressions for the stress com- i<rrr + 0-00) = - - f°°+tC (5$(s)e("'-*') + *$(-s)e(-"«+i"))e-i<'<fc (7.30) _ — r ponents as a function of the distribution of traction along the wedge boundary. Hence, »oo+tc 7rrl "(Orr + CT00) = " ̂ W*) + iPxi*)^*'d* (7.32) ./-oo+tc E>l(S) -̂ 4(1 + w)(Pi(«) - Ps{s))e-*'ds (7.33) -oo+tc -C'lC5) -oo+tc where the functions Pp(s), 3 € {1,2,3} are hyperbolic functions defined in table 7.2. Table 7.2: Functions of a, 9, and s employed in the text. Pi (a) = sin(a — 9) cosh(a + 9)s + sin(a + 9) cosh(a — 9)s P2(s) = cos(a — 9) sinh(a + 9)s 4- cos(a + 9) sinh(a — 9)s P3(̂ ) = sin(a — 9) sinh(a + 0)s — sin(a + 9) sinh(a — 9)s Pi(s) = cos(a — 9) cosh(a + 0)s — cos(a + 9) cosh(a — 9)s D\ (s) = sinh 2as + s sin 2a D2(s) = sinh2a5 — s sin 2a Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 74 Green's function solution Finally, the Green's function is obtained for the above case by locating symmetric point forces, directed normal to the wedge surface, at r = b. Hence, /(r) = 6(r — b), and 7.4 Solution by Mellin transform The biharmonic equation, (7.1), when expressed in polar coordinates, also lends itself to analysis based on the Mellin transform. Here, the transform is employed to obtain solutions to the wedge problem as stated in section 7.2. While there is some redundancy between the method of solution employed here and the complex variable ap- proach described in the previous section, there are also different advantages to be gained from each approach. The complex variable method provides continuity with chapters 5 and 6, further demonstrating the utility of complex variable methods. However, in the complex variable solution, calculation of the Airy stress function is bypassed. This (7.34) where it has been recalled that r = ae*. Letting b/a = e6 and substituting (7.34) into (7.32) and (7.33) gives Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 75 poses practical problems if one is interested in calculating the displacements as well as the stresses from the solution. In contrast, the Airy stress function is obtained rather naturally in the Mellin transform domain, as are the displacement components. Finally, by showing the equivalence of the two solutions, weight is lent to both analyses. 7.4.1 The biharmonic equation in the Mellin domain The Mellin transform, f(s) of a function, /(r) is given by f(s) = M[f(r)] = r fry'-1 dr. (7.37) Jo The inverse formula is given by (7.38) where c is a real constant. The following transform functions are denned: A(5) = M [A] (7.39a) aij(s) = M[r2aij] (7.396) ui(s) = M[rui] (7.39c) /*(*) = M[fi] (7.39d) &(5) = M[9i] (7.39c) where i,j € {a, 8}. The Mellin transform is directly related to the two sided Laplace transform, and, as with that transform, usually exists only for those values of 5 falling within a strip of regularity parallel to the imaginary axis [see, for example, Bracewell, 1986, p.221-223]. Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 76 The upper bound of this strip follows from the regularity requirement, (7.3). Given this condition, (Tie = 0 ( r - 1 + £ ) as r —> oo. Hence, by (7.39b) along with (7.37), > f°°0{r-1+€)r~s+1t Jv <TiO £ I u\r - • - jr" 1 dr 'p = f 0(rs+t)dr (7.40) Jp where p is some positive real number, and e > 0. It follows that existence of &ie requires that 3fi(5) < —1. The lower bound on $l(s) follows upon consideration of the relations given by (7.7), which, along with (7.4) and (7.5), imply that as r —• 0, <Tij — O(r 0). Hence, on > fP0(ri+1)dr. (7.41) Existence oi&ij requires that 3?(5) > —2. Thus, the strip of regularity for the transform functions listed as (7.39a-e) includes values of 5 such that —2 < 9ft(S) < —1. Expanding (7.1) in polar coordinates, and multiplying by r 4 one obtains d dr2 + T d r + de2 0. (7.42) Application of the Mellin transform to this result leads to the fourth order ordinary differential equation (ODE), A = 0 (7.43) In addition, transformation of (5.18) to polar coordinates, and application of the Mellin transform to the resulting expressions yields, by (7.39b) and (7.37), d2A <ree dr2 _ 1 d 2 A l_dA a r 0 ~ r drdO + r 2 86 1 8K 1 d2k 0~rr = ~— h (TOO = s(s + 1)A dk dr ee2 cr^e = (5 + 1) dd2 dB — s (7.44a) (7.446) (7.44c) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 77 The ODE, (7.43) possesses the general solution, A(s, 9) = Ax (s) sin 59 + A2 (s) cos s9 + Az(s) sm(5 + 2)6 + AA(s) cos(5 + 2)9 (7.45) The functions, Ak(5), k € {1,2,3,4}, are set by the conditions placed on the normal and shear tractions acting along the boundary. Substitution of (7.45) into (7.44) to evaluate the stresses on 0 = ± a yields a set of equations for the Ak, sin sa cos sa sin(5 + 2)a cos(s + 2)a "^i(')" — sin sa cos sa — sin(5 + 2)a cos(s + 2)a A 2 ( i ) cos sa — sin sa C O S(5 + 2)a i^±2l S i n ( j + 2)a As (5) cos sa sin sa ( 5 t 2 ) cos(S + 2)a +̂21 sin(5 + 2)a .^4(5). A(») 5(3+1) A(5) 5(5+1) 5(3+1) . ga(*) Lsrl+fj which have the solution (gi + 52)5sin(5 + 2)a - (fi - f2)(s + 2) cos(s + 2)a MS) M^ Az(s) M*) 2s(s + l)d2(s) (gi - g2)scos(s + 2)a + (/i + + 2)sin(5 + 2)a 2s(5 + l)di(5) ~(gi +52) sin sa + (fi - /2)cos sa 2(s + l)d2(s) ~(gi - 52) cos sa - (fi + /2)sin sa (7.46) 2(5 + 1)^(5) (7.47a) (7.476) (7.47c) (7.47d) where di(s) = (5 + 1) sin 2a + sin[2a(s + 1)] d2(S) = (s + 1) sin 2a - sin[2a(s + 1)]. (7.48a) (7.486) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 78 7-4-2 Symmetric normal stress with zero shear II For deformation symmetric about the wedge bisector, and excluding the possibility of any shear acting on the boundary, the boundary conditions in the Mellin domain are = A(i) = /(*), 9i(S) = h{») = 0. The transformed Airy stress function for this case is A f(~\ (̂  2)sin(a + 2)ct cos sO — S sin Sa cos(if + 2)0 . . A = f { 3 ) s(s +1)^(3) ' ( 7 > 4 9 ) Upon letting 5 = —1 + is, and following some simplification, (7.49) becomes, k = fts)*}s) + '™ (7.50) where f*(x) = f(-l+ix). It follows, based on (7.50), (7.39b) and (7.44) that Di[s) Di(s) where, in obtaining the above result, it has been noted that, dP2 d0 dPi de and dP3 Ps + 3 P 4 = - P 4 + sP3 ie + Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 79 Thus, *rr + ho = ^ 4 [P2(s) + iP^s)] (7.52a) oIt + &M + 2i&T0 = 2 ^ | ( 1 + is) [P^s) - P3(s)]. (7.526) Di(s) By combining (7.38) and (7.39b), an expression is obtained for the inverse Mellin transform of the non-deviatoric stress, pc*+too Ji. r ^ «c*-t-too (<jvr + aee) = TT1 I [*» + oree]r~'ds Z m Jc*-ioo or roo-t(c* + l) f* r(<rrr+vee) = - / j-^[P2(s) + iP1{s)]r-"ds. (7.53) 7T ^ . ^ . ^ c . + i) 1̂ 1 IS). For the case of a Dirac-delta type forcing function, f(r) = S(r — b), the transformed forcing is /.OO f(s)= 6(r-b)rs+1dr = bs+1 =6". (7.54) •/o Thence, substitution of f(s) into (7.53) yields ds (7.55) 1 P2(s) + iPi(s){ry" r(arr + aee) = - / prr^ 7 < where c* + 1 has been replaced by c. The regularity condition, (7.3), requires that —2 < c* < —1, or —2 < —c — 1 < —1. Hence, the condition placed on c is that 0 < c < 1. The above result is exactly equivalent to that obtained by the complex variable method and given by (7.35). Hence, the equivalence of the solutions obtained by the two transform methods has been demonstrated. Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 80 7.5 The inverse transform From a mathematical standpoint, the first fundamental problem for the infinite wedge is solved once results of the form of (7.35) and (7.36) have been obtained. How- ever, for application purposes, it is necessary to evaluate the inverse transforms nu- merically. The inverse Fourier transform is carried out in the complex .s-plane, while the inverse Mellin transform involves the variable s. But, as was demonstrated in the previous, section, the inverse Mellin integrals can be readily mapped onto the s-plane by the transformation, S —» —1 + is. Consequently, the quadrature problems for each solution method are equivalent. X X 1 # X -1 X X X Figure 7.3 A typical distribution of poles in the s plane for the case of deformation antisymmetric with respect to the wedge bidector. The crosses are the pole locations. The dashed fine in- dicates a possible contour of integration for the inverse transform. Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 81 The expressions for the Green's funtions for the cases listed in Table 7.1 have the general form r(crrr + aee) = - —^e^-^ds (0 < c < 1) (7.56) where N(s) is a hyperbolic function of s, 7 6 {1,2}, and the D-f(s) are listed in Table 7.2. The integrand of the above expression exhibits poles for values of s which satisfy, Dy(s) = 0, or { sinh 2as + s sin 2a = 0 7 = 1 (7.57) sinh 2as — s sin 2a = 0 7 = 2 The pole locations can be categorized according to their dependence on the wedge half- angle, a. Those whose positions in the s-plane are not dependent on a are termed a-independent, whereas those whose positions vary with a are termed a-dependent. The condition (7.57) admits a-independent poles at s € {0} 7 = 1 (7.58) s € {—i,0,i} 7 = 1 In addition, an infinite number of a-dependent poles exist for each value of a. In general, these poles have to be located by numerical means. For small wedge angles (2a << 7r ) the a-dependent zeros fall outside the strip 0 < 3?(s) < 1, so there are no poles along the path of integration. Also, each pole is mirrored across the two coordinate axes and about the origin by conjugate poles. A typical distribution of poles in the s-plane, along with a possible integration contour, is illustrated in figure 7.4. Integrals of the form, (7.56), can be simplified by dropping the integration contour to the base of the strip 0 < $s(s) < 1, and integrating along the real axis. A small perturbation in the contour is required in order to evade the pole at s = 0, which, Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 82 obeying the rules of integration in the complex plane, is accounted for by subtracting 7rz times the residue at s = 0. Letting s = p + iq, (7.56) becomes ( 7 - 5 9 ) For example, substitution of the N(p) and D^p) applicable to the case of symmetric normal loadings with zero shear, yields nr(<rrr + (T6e) = r P 2 ( P ) n + / f l ( P ) [ «>s(£ - S)p - i sin(<f - S)p] dp (7.60a) i(l + ip) 1 [ c o s ( £ -S)p-isin(£ - S)p]dp. (7.606) The residues at s = 0 for each of the above integrals are readily obtained by L'Hopital's rule. By way of illustration, the residue of the integrand of (7.60a) is calculated below: p(P2(p) + iPi(p))e-^-S^ reso = lim r-— ;— «->o sinh 2ap + p sin 2a r (P2(p) + iPi(p))e-'(^-g)P + sIRjp) = lim —— ;— *-*o 2a cosh 2ap + sin 2a 2i sin a cos 6 2a -f- sin 2a where IR(p) is a term which becomes irrelevant when the limit is evaluated. Noting that Di(p), P2(p) and Pj(p) are odd functions of p, and that Pi(p) is even, (7.60) can be simplified to obtain, 27rsinacos77 f°° P3(s)cos(£ - S)s + Pi{s)sin(f - 6)s 7rr((7rr + tree) = -z—, . 0 +2 : =—-r ds (7.61aJ 2a + sin 2a JQ T>i(s) 2TT sin a cos 77 [°° P i (s ) (cos(£ - S)s - ssin(£ - S)s) ̂  / 7 f i l M 7rr(<7 r r - aee) = ——, . 0 2 1 TTT~\ D S ( ' - o l o ) 2a-fsm2a J0 E>i(s) Jo Dl\s) Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 83 A computer code was written to evaluate integral expressions of the form those present in equations (7.61a-c). The code had to be sufficiently robust to handle inte- grands which were often highly oscillatory. Hence, a fourth order Runge-Kutta method was employed [Press et al., 1988]. The integrals, which are of the form were cast as differential equations of the form, dy _ I(p) dp (i-py with the boundary condition 1(0) = 0. The Runge-Kutte algorithm was then used to step the value of p from 0 to 1; with the value of y(p) at p — 1 being equivalent to the numerical value of the original integral. The above method was used to calculate the stress field corresponding to the Green's function for case (i). Within the accuracy limitations imposed by numerical dispersion, the boundary conditions are well reproduced by the integration method. Other solution examples are given in Chapter 4. 7.5 A note on the wedge paradox In the introduction to the chapter, mention was made of the wedge paradox, which refers to the breakdown of the classic solution to the first fundamental problem for wedge opening angles greater than 180° (or « 257° with anti-symmetric loading). The paradox can be observed in the solutions obtained in the Mellin and Fourier transform planes. Here consideration is given to the solution expressions obtained in the Fourier plane, since any discussion of the expressions in the Mellin domain would be very similar. As is observed in figure 7.3, the strip of regularity in the Fourier domain, 0 < $s(s) < 1, is free of poles. Therefore, there is no ambiguity regarding the location of the line of integration for the inverse transforms. However, at angles greater than the Chapter 7. THE FIRST FUNDAMENTAL PROBLEM FOR THE WEDGE 84 critical angle, an a-dependent pole is observed on the imaginary axis and within the strip of regularity. This then raises a question as to where the line of integration should run, as lines running both above and below this pole satisfy the regularity requirement. Here, this issue is not resolved, since, for the intended geophysical application of the above solutions, results for small wedge angles only are required. Hence, the question of resolution of the wedge paradox in the domain of the Fourier or Mellin transform variable remains open. The most promising method of attacking this problem, as envisioned by this author, will entail a careful analysis of the residue expansion of the transform integrals. Chapter 8 S U M M A R Y 8.1 Model l ing of deformation above a subduction zone In this thesis, two physical models are developed to quantify the deformation which occurs due to coupling at subduction zones. The specific purpose of the models is to provide a means for interpreting geodetic observations of deformation in terms of physical parameters related to the extent and degree of coupling along the plate interface. The effects of this coupling are of particular interest due to the earthquake hazard associated with the accumulation of stress along the interface. Eventually, the magnitude of this stress must reach a level above that which can be supported along the interface, and seismic rupture occurs. The largest earthquakes ever recorded have arisen as a result of failure along subduction faults. A number of general assumptions are instituted prior to discussing model specifics. In chapter two, a conceptual formulation of the stick-slip process is described. Issues such as those associated with the change of rheology with depth and the distribution of traction along the plate interface are considered. In addition, the set of geodetic data is characterized. The end result of this discussion is a simplified, general representation of a subduction zone. The representation is configured to provide a reasonable physical picture of a subduction zone while at the same time keeping the number of physical parameters to a minimum as necessitated by the sparse distribution of the geodetic observations. The first model considered is based on a kinematic representation of the stick- slip cycle which was introduced by Savage [1983]. According to this representation, the deformation due to coupling can be assumed equivalent to that associated with 85 Chapter 8. SUMMARY 86 a dislocation co-located with the coupled portion of the interface. The slip vector of the dislocation is oriented in the direction opposite to that of the relative motion of the two plates, while its magnitude is allowed to grow at the plate convergence rate. Thus, at least kinematically, locking together of the two plates is simulated. From a dynamic point of view, however, it is by no means clear that this formulation adequately simulates the behaviour of the plate interface. In chapter three, the dynamic implications of the dislocation formulation are investigated. Quantitative results are obtained for the case of an elastic, half-space Earth. The model results and the ensuing discussion highlight a number of problems that are inherent to the dislocation approach. The most significant drawback of the dislocation approach is the limitation which this approach places on the possible distributions of stress along the plate interface. By force balance arguments, it is demonstrated that the net shear stress along the plate boundary must total to zero. Therefore, any shear supported along the coupled portion of the interface must be balanced by shear tractions of equal and opposite orientations acting outside this boundary. This implicit character of the stress distribution limits the freedom of the interpreter to investigate different dynamic scenarios. In addition, the updip tractions which must occur downdip of the locked zone appear aphysical given the direction of the relative motion between the two plates. Given the disadvantages of the dislocation approach, an alternative approach is considered in chapter four. The new appoach is based on a dynamic formulation of the subduction process whereby the overthrust wedge portion of the overriding plate is treated as a free body. The deformation within this body is constrained by traction conditions which are introduced along the base of the overthrust. The effects of different force distributions on the surface deformation are examined. Green's functions are plotted point forces with varying orientations and positions at the base of the wedge. It is noted that, because shear forces produce extension at the wedge surface, plate coupling might produce a positive or at least neutral horizontal strain signal directly Chapter 8. SUMMARY 87 above any locked portion of the interface. The most important observation of the chapter pertains to the non-uniqueness inherent in the dynamic interpretation. Similar profiles of surface strain that are fundamentally different in character can be obtained for entirely different distributions of traction along the base of the wedge. Thus, in modelling deformation above subduction zones, a degree of physical insight will always be necessary in order to specify the dynamic character of the stress distribution along the plate interface. 8.2 Theory of two-dimensional elasticity In order to obtain the quantitative results described in the geophysical study de- scribed above, solutions to a pair of problems of two-dimensional elasticity are required. First, the solution to the problem of a uniform edge dislocation within a half-space is required as the fundamental building block of the Green's function approach to the dislocation model within a uniformly elastic earth. Second, the dynamic approach necessitates that a solution be sought to the problem of determining the deformation within a wedge subjected to traction boundary conditions. Both of the above problems are studied using the complex variable approach to two-dimensional elasticity. Hence, the theoretical basis for this approach is outlined in some detail in chapter five. Starting from the basic equations of elasticity, expressions are obtained for the components of stress and displacement in terms of two analytic functions of the complex variable, z = x+iy, where x and y are the coordinate axes. The utility of these expressions in analyses of problems with singular boundary conditions is then demonstrated by solving the problem of a point force applied to an infinite elastic solid. The solution to the problem of an edge dislocation within an elastic half space is described in chapter six. First, a solution is obtained for the case of an edge dislocation within the infnite plane. This solution is extracted easily from the general expressions Chapter 8. SUMMARY 88 for the analytic functions (which were derived in chapter five) because the geometry of the solution domain and boundary conditions play no role in determining the stress field. The half-space solution is obtained by adding a term, holomorphic in the solution domain, to the whole-space solution. The new term satisfies the free surface conditions which now exist along the half-space boundary. The wedge problem is analyzed by two transform methods. The first transform method employs the complex expressions for the stress and displacement components derived in chapter five. The wedge is conformally mapped into the interior of the infinite strip. This domain is then transformed using a complex version of the Fourier theorem. The second method employs the Mellin transform. This transform is applied to the biharmonic equation as expressed in polar coordinates. This equation is then solved in the domain of the transform variable. For both methods, the boundary conditions are satisfied in the transform domain and inverse transform integrals are obtained for the physical parameters. The equivalence of the solution integrals obtained by the two methods is demonstrated by a simple change of variable. Numerical integrations based on the Runga-Kutta method are then performed to obtain quantitative results. REFERENCES Ando, M. 1975. Source mechanisms and tectonic significance of historical earthquakes along the Nankai trough, Japan. Tectonophysics, 27, 119-140. Angermann, J. K., C. Reigber, and J. Reinking 1994. SAGA South American geody- namic activities: results of first GPS campaigns (abstract). EOS supplement, 75, 163. Bogy, D. B. 1968. Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J. appl. Mech., 35, 460-466. Bracewell, R. N. 1986. The Fourier Transform and its Applications, 2nd edition, McGraw-Hill. Brahtz, J. H. A. 1933. Stress distribution in wedges with arbitrary boundary forces. Physics, 4, 56-65. Byrne, D. E., D. M. Davis, and L. R. Sykes 1988. Loci and maximum size of thrust earthquakes and the mechanics of the shallow region of subduction zones. Tecton- ics, 7, 833-857. Carothers, S. D. 1912. Plane strain in a wedge. Proc. Roy. Soc. Edin., 23, 292. Chinnery M. A. 1963. The stress changes that accompany strike-slip faulting. Bull. Seism. Soc. Am., 53, 921-932. Churchill, R. V., J. W. Brown, and R. F. Verhey 1976. Complex Variables and Appli- cations, McGraw-Hill. Cohen, S. C. 1994. Evaluation of the importance of model features for cyclic deformation due to dip-slip faulting. J. Geophys. Res., 1 0 0 , 2031-2038. Cohen, S. C , S. Holdahl, D. Caprette, S. Hilla, R. Safford, and D. Schultz 1995. Uplift of the Kenai peninsula, Alaska, since the 1964 Prince Wiliam Sound Earthquake. Geophys. J. Int., 1 1 9 , 831-841. Davies, G.F. 1980. Mechanics of subducted lithosphere. J. Geophys. Res., 85, 6,304- 6,318. Davies, G.F. 1983. Subduction zone stresses: constraints from mechanics and from topographic and geoid anomalies. Tectonophysics, 9 9 , 85-98. Dempsey J. P. 1981. The wedge subjected to tractions: a paradox resolved. Journal of Elasticity, 1 1 , 1-10. Douglass, J. J., and B. A. Buffett 1995. The stress state implied by dislocation models of subduction deformation. Geophys. Res. Lett., in press. 89 REFERENCES . 90 Dragert, H., R. D. Hyndman, G. C. Rogers, and K. Wang 1994. Current deformation and the width of the seismogenic zone of the northern Cascaadia subduction thrust. J. Geophys. Res., 99, 653-668. Drew, J. D., and R. M. Clowes 199. A re-interpretation of the seismic structure across the active subduction zone of western Canada, in Studies of Laterally Heteroge- neous Structures Using Seismic Refraction and Reflection Data, ed. A. G. Green; Geological Survey of Canada, Paper 89-13, p.115-132. Eiby, G. A. 1980. Earthquakes, Van Nostrand Reinhold. Forsyth D., S. Uyeda 1975. On the relative importance of the driving forces of plate motion. Geophys. J. R. astr. Soc, 43, 163-200. Freund, L. B., D. M. Barnett 1976. A two-dimensional analysis of surface deformation due to dip-slip faulting. Bull. Seismol. Soc. Am., 66, 667-675. Fung, Y. C. 1965. Foundations of Solid Mechanics, Prentice-Hall, Englewood Cliffs, New Jersey. Green, A. E. , and W. Zerna 1968. Theoretical Elasticity, Clarendon Press, Oxford. Hager, B. H., R. W. King, and M. H. Murray 1991. Measurement of crustal deformation using the global positioning system. Annu. Rev. Earth Planet. Sci., 19, 351-82. Hasegawa, A., S. Horiuchi, and N. Umino 1994. Seismic structure of the northeastern Japan convergent margin: a synthesis. / . Geophys. Res., 99, 22,295-22,311. Hyndman, R. D., and K. Wang 1993. Thermal constraints on the zone of major thrust earthquake failure: the cascadia subduction zone. J. Geophys. Res., 98, 2,039- 2,060. Hyndman, R. D., K. Wang, and M. Yamano 1995. Thermal constraints on the seismo- genic portion of the southwestern Japan subduction thrust. J. Geophys. Res., 100, 15,373-15,392. Larson, K. M. , and M. Lisowski 1994. Strain accumulation measurements in the Shu- magin Islands: results of initital GPS measurements. Geophys. Res. Lett., 21(6), 489-492. Levy, M. 1898. Sur la legitime de la regie dite du trapeze dans l'etude de la resistance des barrages en magnnerie. Compte. Rend., 126, 1235-1240. Lisowski, M. , J. C. Savage, W. H. Prescott, and W. K. Gross 1988. Absence of strian accumulation in the Shumagin seismic gap, Alaska, 1980-87. J. Geophys. Res., 93, 7,909-7,922. Marsden, J. E . 1973. Basic Complex Analysis, W. H. Freeman and Company, San Francisco. REFERENCES . 91 Matsu'ura, M. , and T. Sato 1989. A dislocation model for the earthquake cycle at convergent plate boundaries. Geophys. J., 96 , 23-32. Moore, G.F., and 7 others 1990. Structure of the Nankai trough accretionary zone from multichannel seismic reflection data. J. Geophys. Res., 95 , 8,753-8,765. Muskhelishvili, N. I. 1963. Some Basic Problems of the Mathematical Theory of Elas- ticity, P. Noordhoff, Groningen, The Netherlands. Okada, Y. 1985. Surface deformation due to shear and tensile faults in a half space. Bull. Seism. Soc. Am., 75, 1,135-1,154. Okada, Y. 1992. Internal deformation due to shear and tensile faults in a half space. Bull. Seism. Soc. Am., 82 , 1,018-1,040. Press, F. 1965. Displacements, strains, and tilts at teleseismic distances. J. Geophys. Res., 70 , 2,395-2,412. Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vettering 1988. Cambridge University Press. Rybicki, R. 1986. Dislocations and their geophysical application, in Continuum Theo- ries in Solid Earth Geophysics, ed. R. Teisseyre, Elsevier, p.18-186. Savage, J. C. 1983. A dislocation model of strain accumulation at a subduction zone. J. Geophys. Res., 88 , 4,984-4,996. Savage, J. C. 1995. Interseismic uplift at the Nankai subduction zone, southwest Japan, 1951-1990. J. Geophys. Res., 100, 6,339-6,350. Savage, J. C , M. Lisowski, and W. H. Prescott 1991. Strain accumulation in Western Washington. J. Geophys. Res., 96 , 14,493-14,507. Savage, J. C , and G. Plafker 1991. Tide gueage measurements of uplift along the south coast of Alaska. J. Geophys. Res., 96 , 4325-4335. Savage, J. C , and W. Thatcher 1992. Interseismic deformation at the Nankai trough, Japan, subduction zone. J. Geophys. Res., 97 , 11,117-11,135. Shepherd W. M. 1935. Stress systems in an infinite sector. Proc. Roy. Soc. Lon. series A, 148, 284. Shimenda A. I. 1993. Subduction of the lithosphere and back arc dynamics: insights from physical modelling. J. Geophys. Res., 98 , 16,617-16,185. Sokolnikoff, I. S. 1956. Mathematical Theory of Elasticity, 2nd edition, McGraw-Hill. Stefanick, M. , and D. M. Jurdy 1992. Stress observations and driving force models for the South American plate. 97 , 11,905-11,913. Steketee, J. A. 1958a. On Volterra's dislocation in a semi-infinite elastic medium. Can. J. Phys., 36 , 192-205. REFERENCES . 92 Steketee, J. A. 1958b. Some geophysical applications of the elasticity theory of dislo- cations. Can. J. Phys., 36, 1,168-1,198. Sternberg E. , and W. T. Koiter 1958. The wedge under a concentrated couple: a paradox in the two-dimensional theory of elasticity. J. appl. Mech., 25, 575-581. Thatcher, W., and J. B. Rundle 1984. A viscoelastic coupling model for the cyclic defor- mation due to periodically repeated earthquakes at subduction zones. J. Geophys. Res., 89, 7,631-7,640. Tranter, C. J. 1948. The use of the Mellin transform in finding the stress distribution in an infinite wedge. Q. J. Mech. appl. Math., 1, 125-130. Whittaker, A., M.H.P. Bott, and G.D. Waghorn 1992. Stresses and plate boundary faults associated with subduction plate margins. / . Geophys. Res., 97, 11,933- 11,944. Wolfram Research, Inc. 1992. Mathematica. version 2.1.


Citation Scheme:


Usage Statistics

Country Views Downloads
Japan 4 0
China 4 8
United States 2 1
City Views Downloads
Tokyo 4 0
Beijing 4 0
Sunnyvale 1 0
Unknown 1 3

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}


Share to:


Related Items