Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Postseismic deformation following the 1991 Racha, Georgia earthquake Podgorski, Joel Edwin 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2006-0619.pdf [ 12.52MB ]
Metadata
JSON: 831-1.0052519.json
JSON-LD: 831-1.0052519-ld.json
RDF/XML (Pretty): 831-1.0052519-rdf.xml
RDF/JSON: 831-1.0052519-rdf.json
Turtle: 831-1.0052519-turtle.txt
N-Triples: 831-1.0052519-rdf-ntriples.txt
Original Record: 831-1.0052519-source.json
Full Text
831-1.0052519-fulltext.txt
Citation
831-1.0052519.ris

Full Text

POSTSEISMIC DEFORMATION FOLLOWING THE 1991 RACHA, GEORGIA EARTHQUAKE by Joel Edwin Podgorski B.S., University of California at Santa Cruz, 1997 A THESIS SUBMITTED IN PARITAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Geophysics) THE UNIVERSITY OF BRITISH COLUMBIA July 2006 ©Joel Podgorski, 2006 Abstract The April 29, 1991 Racha, Georgia M w=6.9 earthquake was the largest earthquake recorded in the western section of the Greater Caucasus Mountains. Lithospheric de-formation following this earthquake was recorded by a network of eight GPS stations. These data were modeled for postseismic deformation with the mechanisms of afterslip and viscoelastic relaxation. Prior to modeling, the GPS data were corrected for secular tectonic motions by two separate methods. The first was to use the 1996-2000 portion of the GPS time-series as an estimate of preseismic deformation rates. The second method involved using a regional tectonic block model to produce velocity estimates for the GPS stations. Shallow afterslip was found to best explain the data for the 1991-1994 time period with a moment of 7 x 10 1 8 N m. This is equivalent to about 20% of the coseismic moment and represents 300 times the total moment released from aftershocks during the same time period. Most of the postseismic deformation was completed by 1994. While viscoelastic relaxation was not found to be responsible for postseismic de-formation, it was determined that viscosities less than 10 1 8 Pa sec can be ruled out for possible a possible viscoelastic layer in the lower crust. ii Contents Abstract ii Contents iii List of Tables vi List of Figures viii 1 Introduction 1 2 Background 3 2.1 Tectonic Setting 3 2.2 Geology 6 2.3 Lithospheric Structure 7 2.4 Seismicity 9 3 Racha Event 11 4 Postseismic Deformation Processes 14 iii 4.1 Fault Afterslip 15 4.2 Viscoelastic Relaxation 19 4.3 Other Postseismic Processes 21 5 Other Postseismic Deformation Studies 24 6 GPS Data 28 6.1 Stations and Time Points 29 6.2 Acquisition 30 6.3 Analysis 31 6.4 Position Measurements 31 6.5 Velocity Decay in GPS Time-Series 33 6.6 Secular Correction 37 6.7 Error Combination . . . > 42 6.8 Data Used in Modeling 42 7 Methods 44 7.1 Afterslip Inversion Modeling 44 7.2 Viscoelastic Forward Modeling 48 8 Results 53 8.1 Afterslip 53 8.2 Viscoelastic Relaxation Modeling 70 9 Discussion 75 iv 9.1 Secular Correction of GPS Data 75 9.2 Afterslip and Viscoelastic Relaxation Models 77 10 Conclusions 83 Bibliography 85 Appendix A GPS Position Measurements 95 Appendix B Conversion Velocities for ITRF2000 to Eurasian Refer-ence Frame 98 Appendix C GPS Displacements 99 Appendix D Earth Model for Viscoelastic Modeling 102 V List of Tables 2.1 Crustal P and S wave velocities from Triep et al (1995) compared with typical P and S velocities from the PREM (Dziewonski and Anderson, 1981) 8 3.1 Racha source parameters according to different sources. Moment is listed with total coseismic moment above and mainshock-only mo-ment below in parentheses 12 6.1 Study duration, number of GPS stations, and average number of ob-servations of other postseismic deformation studies 29 6.2 Residuals (WSTD) and linear rates from linear and logarithmic func-tions fit to data. Amplitude of logarithmic function is also included. All values are in mm or mm/yr 36 6.3 Block model velocities (mm/yr) used in correcting GPS time-series for secular tectonic displacements 42 vi 7.1 Coseismic slip on 4 km x 4 km fault patches from model of Tan (2005) 52 8.1 Sensitivity of afterslip results to choice of fault plane as measured by changes in WRSS reduction explained by inversion models. In-versions are based on TSC data and a B (smoothing) of 10 km2/m. Fault geometries begin with the parameters of Tan (2005) and are adapted incrementally to the models of (a) Fuenzalida et al (1997), (b) Triep et al (1995), and (c) Harvard CMT. 60 8.2 Comparison of WRSS reduction and moment resulting from models with different B values. The first two B values for each data set are those determined by the methods of analyzing the trade-off between misfit and roughness and of cross validation sum of squares (CVSS). The third B value, a more conservative 10 km2/m, is that which was used in many of the tests presented in the results. This value had been chosen due to it allowing for a more reasonable-appearing slip distribution 65 8.3 Test for the sensitivity to Poisson's ratio, v, by increasing and de-creasing v from the value of 0.25 otherwise used. Differences in WRSS reduction and moment are indicated 68 vii List of Figures 1.1 Generalized tectonic setting of the Caucasus Mountain region. The Anatolia and Azerbaijan blocks displace to the west and east re-spectively while the Arabian plate proceeds northward into Eurasia (Lesser and Greater Caucasus) 2 2.1 Selected large-scale tectonic and geologic features of the Caucasus Mountains. MCT is the Main Caucasus Thrust. GDZ is the Gagra-Dzhava Zone of Jurassic sedimentary rocks just north of the Racha epicenter which is indicated by a star. DM is the crystalline Dzir-ula Massif on the south side of the Racha epicenter. BKSZ is the Borzhomi-Kazbeg Shear Zone 5 viii 4.1 Rate and state friction. Friction, /i, is represented in the vertical di-rection and displacement in the horizontal direction, a is the increase in friction due to a sudden acceleration, and b is the decrease from this peak by the end of the period of greater speed, a - b determines if a medium experiences (a) velocity weakening (negative) or (b) velocity strengthening (positive) 17 6.1 Location of GPS sites (triangles). Solid rectangle indicates the mod-eled coseismic slip plane and star the epicenter of Tan (2005), and dashed plane indicates rectangle of afterslip plane. 30 6.2 Map of part of the tectonic block model used in producing secular velocity estimates (Reilinger et al, 2006). The area of this study as well as the straight line estimate of the Main Caucasus Thrust plate boundary are indicated 40 7.1 Profile perpendicular to strike of afterslip fault plane with geometry of Tan (2005) and two months of aftershocks. 47 7.2 Coseismic fault model, from Tan (2005). Displacement is in meters. Each tile is 4 km square. 47 ix 7.3 Fault plane resolution tests using alternating 1 and 0 meter slip patches on a 120 km x 40 km plane, (a) 60 km x 20 km and (b) 20 km x 20 km sized slip patches were forward modeled (top) and the result-ing surface displacement vectors were inverted for slip (bottom) to determine which areas of the modeled fault plane the GPS network can resolveslip on. Inversions were performed using a smoothing, B, of 10 km2/m except where indicated 49 8.1 1991-1994 afterslip inversion results using (a) TSC data and (b) BMC data, (c) Coseismic slip model of Tan (2005). Star on each slip plane represents the Racha hypocenter, and each tile is 4 km square. Comparing either afterslip model (a) or (b) with the coseis-mic slip model (c) shows that the bulk of modeled afterslip occurs away and up-dip from the area of maximum coseismic slip 54 8.2 Comparison of 1991-1994 data vectors with 1-cr error ellipses (red) with model vectors (blue). Model vectors are produced by forward modeling the afterslip results, (a) is from using TSC data and (b) is from using BMC data. The star represents the Racha epicenter. . . . 55 8.3 1991-1994 slip inversion using data with no secular correction. WRSS reduction is 51.7435% and the moment is 8.987 xlO 1 8 N m 56 8.4 1991-1994 afterslip model and data vectors for all available sites for (a) TSC data (b) BMC data and (c) non-corrected (NC) data. . . . . 57 x 8.5 1991-1994 slip inversions using the same GPS sites for all data sets (without site NICH). (a) TSC model (b) BMC model (c) NC model. 58 8.6 Surface outline of the fault plane using the fault geometries of Tan, Fuenzalida, Triep, and Harvard CMT. Stars are the Racha epicenter according to the different groups 61 8.7 Patterns of modeled afterslip using different published fault geome-tries. Star indicates the Racha hypocenter. (a) Tan (2005), (b) Fuen-zalida et al (1997), (c) Triep et al (1995), (d) Harvard CMT. 62 8.8 Determination of smoothing parameter, B, for TSC data (a,b) and BMC data (c,d). The first method used, a "Tikhonov plot", de-termines where the change between the weighted residual sum of squares (WRSS) and roughness is greatest (a,c) and then use the associated 8. The second method is cross validation sum of squares (CVSS) where 8 is chosen based on which value produces the small-est sum of residuals. Between the two methods, the optimal value of B ranges from 2 to 4 km2/ni for the TSC data and 3 to 7.5 km2/m for the BMC data 66 8.9 Afterslip inversions resulting from use of different smoothing val-ues, B, for TSC and BMC data. The first two smoothing values dis-played for each data type are those corresponding to the methods of WRSS versus roughness and of CVSS (see Figure 8.8) 67 8.10 1991-1994 afterslip inversions with dip-slip and strike-slip compo-nents using (a) TSC data and (b) BMC data 69 xi 8.11 1994-1996 afterslip inversion results (a) using data TSC data and (b) using BMC data. Star on each slip plane represents the Racha hypocenter, and each tile is 4 km square 70 8.12 Sensitivity of WRSS to lower-crust viscoelastic layer thickness and viscosity. The best model with the TSC data (a) has a viscoelastic lower crust of 10 km thickness and a viscosity of about 1.0 x 1017 Pa sec. The best model (star) with BMC data (b) has a viscoelastic lower-crustal thickness of 10-18 km and a viscosity of about 2.5 -6x 1017 Pa sec 71 8.13 Displacement vectors from best fitting ViscolD models and 1991-1994 GPS site displacements, (a) TSC data with a viscoelastic lower crust at 35-45 km with a viscosity of 1 xlO 1 7 Pa sec. (b) BMC data with a viscoelastic lower crust at 25-45 km with a viscosity of 9 x 10 Pa sec. Star is the Racha epicenter. 72 8.14 WRSS versus logarithm of viscosity for forward-modeled viscoelas-tic mantle half-space models compared with (a) TSC data and (b) BMC data 73 8.15 Comparison between BMC data and models of a 20 km thick vis-coelastic lower crust (25-45 km depth), (a) WRSS versus logarithm of viscosity, (b) Corresponding WRSS reduction versus logarithm of viscosity 74 xii (a) Aftershocks over the two months following the Racha event (May 10, 1991 to July 9, 1991) (Fuenzalida et al, 1997). Afterslip fault plane is outlined. Star represents Racha hypocenter and can be used with fault plane outline to compare aftershocks with 1991-1994 (b) TSC afterslip and (c) BMC afterslip. The afterslip models are from data covering the period July 25, 1991 to October 10, 1994 xiii Chapter 1 Introduction The April 29,1991 Racha, Georgia Mw=6.9 earthquake was the largest recorded earthquake in this western section of the Greater Caucasus Mountains. This region comprises a section of the complex Arabian-African-Eurasian plate col-lision zone. Soon after the earthquake eight Global Positioning System (GPS) stations were deployed in and around the earthquake rupture area. The focus of this study is the use of measurements from these GPS stations to constrain mod-els that describe accelerated crustal deformation resulting from this earthquake. A reason for doing this is to place constraints on what deformation mechanisms and rheologies may exist in this part of the continental collision zone. This knowledge will help create a clearer picture of the mechanisms at work in this plate boundary zone in particular and in continental collision zones in general. 1 40°N 45°N Figure 1.1: Generalized tectonic setting of the Caucasus Mountain region. The Anatolia and Azerbaijan blocks displace to the west and east respectively while the Arabian plate proceeds northward into Eurasia (Lesser and Greater Cauca-sus). 2 Chapter 2 Background 2.1 Tectonic Setting The Caucasus Mountains are contained within the northern part of the middle segment of the Alpine-Himalayan fold belt and are created by the collision of Arabia into Eurasia (e.g., McKenzie, 1970) (Figure 1.1). This results in both the Greater and Lesser Caucasus being under a compressional stress regime (Philip et al, 1989). The Caucasus Mountains are bounded on the west and east by oceanic crust of the Black Sea and South Caspian basins respectively and are formed by many nappes thrust generally from north to south. Exceptions to this include northward (in addition to southward) thrusting in the eastern Greater Caucasus and in the far western portion of the western Greater Caucasus. The region contains a string of recent volcanic features that are aligned roughly N-S from the Arabian plate to the Greater Caucasus, suggesting E-W extension 3 (Fuenzalida et al, 1997). In particular there is a cluster of weak seismicity at 41.5°N, 44°E in the Lesser Caucasus that corresponds to N-S trending tensile volcanic fractures (Philip et al, 1989). The Eurasian and Arabian plates converge at a rate of 25 ± 3 mm/yr along N20°W near the Caucasus (at 38°N, 40°E) (DeMets et al, 1990; DeMets, 1994; Smith et al, 1994; Reilinger et al, 1997; Reilinger et al, 2006). Subduction of continental crust is believed to be avoided by lateral extrusion of the Anatolian block, Azerbaijan block, and northwest Iranian block {McKenzie, 1972; Jackson, 1992) and by underthrusting of accreted terrains (Triep et al, 1995). The oblique convergence relative to the Caucasus can be explained by two models: 1- indenter tectonics requiring a left-lateral strike-slip zone along the western boundary of the collisional region at the Borzhomi-Kazbeg shear zone (BKSZ) {Philip etal, 1989). (See Figure 2.1) 2- oblique convergence partitioned between Caucasus-perpendicular collision along this range and right-lateral mo-tion along ESE-WNW faults further south {Jackson, 1992). The BKSZ may be followed southwest into the Northeast Anatolian Fault. In the Greater Caucasus the shear zone splits into several segments, most notably along the Kazbeg volcano and also along the western end of the Alazani Basin {Fuenzalida et al, 1997). There are many differences between the western and eastern Greater Cau-casus, divided by the BKSZ at around 44°E. Philip et al (1989) point out the following. The western Greater Caucasus has concentrated deformation along its southern slope but simply a smooth monocline on its northern slope while 4 Figure 2.1: Selected large-scale tectonic and geologic features of the Caucasus Mountains. M C T is the M a i n Caucasus Thrust. G D Z is the Gagra-Dzhava Zone of Jurassic sedimentary rocks just north of the Racha epicenter which is indicated by a star. D M is the crystalline Dzirula Massif on the south side of the Racha epicenter. B K S Z is the Borzhomi-Kazbeg Shear Zone. the eastern Greater Caucasus has concentrated deformation on both sides. Pl io-Quaternary volcanism occurs only in the western Greater Caucasus (from E l -brus to Kazbek). Thrust faults on the northern boundary of the western Greater Caucasus dip to the south while those on the northern boundary of the eastern Greater Caucasus dip to the north. There is much more seismicity in the east-ern Greater Caucasus, and the Racha event represents the largest known event 5 in the west. Intermediate depth earthquakes are found in the eastern Caucasus (possibly related to the remnant of a subducted slab (Triep et al, 1995)) but not in the western Caucasus (where there is Quaternary volcanism that could also indicate recent subduction). Bouguer gravity anomalies and topography suggest a northward displacement of the eastern Greater Caucasus relative to the western part (Ruppel and McNutt, 1990). Finally, unroofed basement is exposed in the west but not in the east (Triep et al, 1995). 2.2 Geology The highest peaks in the Greater Caucasus, Elbrus and Kazbek, are located in the western part of the range and are Pliocene-Pleistocene volcanoes of calc-alkaline composition that have been active from late Miocene until historical times (McKenzie, 1972). The core of the Greater Caucasus's Main Range con-sists of highly metamorphosed Precambrian rocks (790 Ma) that are considered the remnants of a microcontinent. The northern Peredovoy Range contains Pa-leozoic ophiolites and island-arc complexes that represent a suture zone (Zonen-shain et al, 1990). The majority of the Greater Caucasus is underlain by Jurassic and Cretaceous rocks. The lower and middle Jurassic sequences are mostly shales and slates with the shale sequences containing lava layers up to 3 km thick. The oldest of these lavas is of calc-alkaline composition and is a remnant of a Great Caucasian island arc. Abundant basalts from the early and middle Jurassic occur in the 6 middle Greater Caucasus and are associated with rhyolite. This basalt-rhyolite series is interpreted to show the formation of the Great Caucasian sedimentary basin by extension (Zonenshain et al, 1990). On the southern side of the Greater Caucasus, continuous sedimentation is found from the middle Jurassic through the Cretaceous and Paleogene with no uplift and occurred in a marginal sea north of the Lesser Caucasus (Kopp and Shcherba, 1985). This shows that uplift of the Greater Caucasus began in the middle to late Miocene, and it has continued until the present. The rocks on the north side of the Racha rupture area are the Jurassic Gagra-Dzhava thrust sheet (Triep et al, 1995) (See Figure 2.1). To the south of the Jurassic through Paleogene sequences is the crystalline Dzirula Massif that shows similarities to the late Paleozoic Hercynian structures in Europe (Kopp and Shcherba, 1995; Dotduyev, 1986) and is different from the Gondwanan microcontinents to the south. It is overthrust on the south by the Lesser Caucasus (Burtman, 1989) and on the north by the Greater Caucasus. Pa-leomagnetic data show an intermediate position between Eurasia and Gondwana (Asanidze et al, 1980). 2.3 Lithospheric Structure It is uncertain exactly what the depth of the Moho is beneath the Greater Cauca-sus. Bouguer gravity studies indicate a Moho depth of up to 60 km (Ruppel and McNutt, 1990) with a depth of about 45 km in the vicinity of the Racha rupture 7 area (Balavadze et al, 1979). That is, once the effects of elevation and the mass of the mountains were accounted for, anomalously low gravity measurements implied a thicker than normal crust. As well, in their Racha source study Triep et al (1995) used a model borrowed from the Georgian Academy of Sciences with a Moho depth of 45 km (which they adapted to a 44 km depth), though it is unclear on what basis this model was formed. Triep et al (1995) calculated crustal P and S wave velocities by inversion of P and S arrival times from the Racha mainshock and aftershocks. These are shown in Table 2.1 and are compared with crustal values of the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson, 1981). Triep et al's velocities are generally less than those of the PREM, suggesting a lower crustal density due to either less dense material, higher temperatures, or both. Crustal P and S Velocities Triep et al PREM Depth P Velocity 5 Velocity P Velocity S Velocity (km) (km/sec) (km/sec) (km/sec) (km/sec) 0-3 4.83 2.61 5.80 3.20 3-7 5.44 3.19 5.80 3.20 7-19 5.98 3.51 6.13 3.43 19-45 6.33 3.65 6.80 3.90 Table 2.1: Crustal P and S wave velocities from Triep et al (1995) compared with typical P and S velocities from the PREM (Dziewonski and Anderson, 1981). An abrupt boundary in seismic attenuation of P and S waves exists between the stable Russian platform to the north and the Greater Caucasus (Sarker and 8 Abers, 1998). According to laboratory experiments the detected amount of at-tenuation corresponds to a temperature increase along the base of the crust of 70-400°C relative to the platform (Kampfinann and Berckhemer, 1985; J. Jack-son et al, 1992). Sarker and Abers (1998) point out that if the attenuation is due to increased temperatures that lower densities would exist (for a given chem-ical composition) thereby rendering Bouguer gravity estimates of Moho depth inconclusive. There is also tomographic attenuation evidence from Pn waves of a partially eroded to completely absent lithospheric mantle beneath the western Greater Caucasus (Al-Lazki et al, 2003). This suggests a thin lithosphere with little or no mantle component underlain by an upper mantle asthenosphere. 2.4 Seismicity There have been several large earthquakes in or near the Caucasus Mountains during the latter part of the 20th century. The 1970 Dagestan thrust event (Ms=6.5) took place on the northern slope of the eastern Greater Caucasus. In 1988 the Spitak earthquake (Ms=6.9) occurred along the southern thrust of the Lesser Caucasus. The 1990 Manjil (Iran) earthquake (Ms=7.3) took place by the south-western coast of the Caspian Sea and showed a left-lateral mechanism con-sistent with the expulsion of the northwest Iranian block. Another strike slip event (Ms=6.8) occurred in 1992 in Erzincan (Turkey) along the North Anato-lian Fault. In general, in the Greater Caucasus seismicity is concentrated along 9 its southern boundary where surface elevations are less and deformation is able to more easily occur (Allen et al, 2004). 10 Chapter 3 Racha Event The M w 6.9 Racha thrust earthquake (Ms=7.0) occurred on April 29, 1991 at 09:12:48.1 GMT along the southern slope of the eastern end of the western Greater Caucasus Mountains where the Greater Caucasus overthrust the Dzirula Massif (see Figure 2.1). The event can be divided into four distinct subevents with a total duration of 22 seconds (Fuenzalida et al, 1997). The rupture area is located in the Gagra-Dzhava thrust sheet of Jurassic sedimentary rocks which is located about 40 km southwest of the main thrust of the Greater Caucasus. This is about 120 km north of Tbilisi in the Republic of Georgia and partly in the district of Southern Ossetia. No surface rupture was observed, though there were large landslides, new groundwater sources, and cracks on top of anticlines (Fuenzalida et al, 1997). Triep et al (1995) proposed two possible fault surfaces: 1) a contact between crystalline basement of the Dzirula Massif below and the Jurassic Gagra-Dzhava thrust sheet complex above or 2) a fault contained en-11 tirely within the Dzirula Massif. Different research groups have reported different earthquake parameters. Those of Fuenzalida et al (1997), Tan (2005), Triep et al (1995), the Harvard Centroid Moment Tensor (CMT) catalog, and the National Earthquake Information Center (NEIC) are summarized in Table 3.1. Total coseismic moment, mainshock-only moment, and errors are indicated where available. The moment for NEIC was calculated from their reported magnitude, Ms = 7.30, using the relation: log M 0 = l.5Ms + 16.1 (3.1) RACHA SOURCE PARAMETERS Group Strike Dip Rake Depth (km) Moment (N m) Lat, Lon Fuenzalida 286.7 29 92.4 4.81 3.146 x 1019 42.42°N etal ±1.3° ±0.6° ±1.2° ±0.2 (2.27 x 1019) 43.69°E Tan 287° 30° 90° 3.195 x 1019 42.42°N o (1.952 x 1019) 43.67°E Triep 292.9 24 97.7 4.42 n/a 42.4238°N etal ±6.2° ±2.9° ±7.4° ±2.1 (2.25 x 1019) 43.6643°E Harvard CMT 288° 39° 106° 22.3 3.34 x 1019 (n/a) 42.60°N 43.61°E NEIC n/a n/a n/a 17 1.12 x 1020 (n/a) 42.453°N 43.673°E Table 3.1: Racha source parameters according to different sources. Moment is listed with total coseismic moment above and mainshock-only moment below in parentheses. Among Fuenzalida et al (1997), Tan (2005), and Triep et al (1995) the total 12 coseismic moments reported are similar as are the mainshock-only moments. All three groups of researchers obtained their results through inversion of teleseismic body waves. The Harvard C M T and, to a much greater extent, the NEIC reported magnitudes are larger than those of the aforementioned studies. The hypocentral depths of these latter two sources are also significantly deeper. These differences may arise due to the Harvard C M T method being being based on longer period waves (free oscillations) and the NEIC result having larger uncertainties due to it being a preliminary estimate that has not necessarily been thoroughly analyzed. 13 Chapter 4 Postseismic Deformation Processes Postseismic deformation of the earth's surface shows the earth's time-dependent response to a step stress change induced by an earthquake. Modeling postseismic deformation can help place constraints on the rheology and structure of the crust and upper mantle in the vicinity of the earthquake (e.g. Freed et al, 2006; Hearn et al, 2002). Postseismic deformation resulting from a shallow crustal earthquake can be observed on a time scale of months to decades and often occurs in a different sense to that of the background motions driving the long-term deformation. This postseismic deformation is typically described by one of two main processes: fault afterslip or viscoelastic relaxation. It is possible for both processes to be responsible for observed deformation, though afterslip usually occurs in the im-mediate months to years following an earthquake while the effects of viscoelastic relaxation can last for decades (Pollitz, 1997). 14 4.1 Fault Afterslip Fault afterslip is the accelerated motion that is often observed on the surface of the earth immediately following a large earthquake. This is detected through some means of geodetic measurement, almost exclusively nowadays by the Global Positioning System (GPS). Afterslip can occur up to several years following an earthquake and decays logarithmically with time for one-dimensional models (Marone et al, 1991). It is usually observed to occur in areas adjacent to the coseismic rupture. The generally agreed upon explanation of fault afterslip centers on the prin-ciple of rate and state variable friction (Dietrich, 1979; Ruina, 1983) which ex-plains the varying friction on a fault associated with coseismic and postseismic slip. This posits that friction along an interface is dependent on the speed at which the surfaces slide past each other as well as the state of the interface, which can be considered to be the average age of contacts on the surface (Diet-rich, 1979): H=li(V,e)=no + aln(V/Vo) + b\n(V08/Dc) (4.1) where Vo is a reference velocity and Dc is the critical distance over which the surface must move to establish a new set of uncorrelated contacts. In a tectonic setting Vo is the preseismic creep rate and is less than the long-term plate velocity on a fault that regularly experiences earthquakes. The state variable 6 is defined 15 as e = DC/V (4.2) If 6 is set to Dc/Vo in Equation 4.1, then // is equal to simply // 0 , that is, the friction at some reference velocity, Vo. Otherwise, plugging Equation 4.2 for 0 into Equation 4.1 yields an equation for friction at a steady-state velocity, V (different from the reference velocity Vo): Here the sign of (a - b) determines if the friction becomes less than the back-ground friction, / i 0 , (velocity weakening) or greater than ju 0 (velocity strengthen-ing). The quantities a and b are determined experimentally with a representing the observed increase in friction due to a sudden increase in sliding velocity, known as the direct effect, and b representing the decrease in friction from this peak to a stable level (see Figure 4.1). This means, for example, that a fault patch in a velocity weakening zone will experience a decrease in friction during the slip of an earthquake and will then itself be more likely to fail. For this reason earthquakes can nucleate only in zones of velocity weakening. Velocity strengthening, on the other hand, means that friction will rise during increased slip so that the rupture of an earthquake dies off in a velocity strengthening zone as the increased friction makes this zone only less likely to slip. The means, however, that a velocity strengthening HM =Ho + (a-b)\n(V/V0) (4.3) \ 16 Figure 4.1: Rate and state friction. Friction, p, is represented in the vertical direction and displacement in the horizontal direction, a is the increase in friction due to a sudden acceleration, and b is the decrease from this peak by the end of the period of greater speed, a-b determines if a medium experiences (a) velocity weakening (negative) or (b) velocity strengthening (positive). 17 zone absorbs some of these earthquake stresses which must be then dissipated by enhanced creep (as opposed to sudden failure) along the fault. This is afterslip. The uppermost layer of the crust is a velocity strengthening zone. Its depth can extend down to several kilometers, and it is usually associated with an area of fault gouge and/or unconsolidated sediments (Marone et al, 1991). The bound-ary between velocity strengthening and velocity weakening can be inferred from the cutoff depth of aftershocks. Above this cutoff depth, postseismic stresses are relieved not by aftershocks but rather by afterslip. Beneath the seismogenic zone (generally 3-12 km depth) is also an area of velocity strengthening friction. Tse and Rice (1986) showed that the onset of this lower zone can be explained by the geothermal gradient. That is, at temperatures at and above 300°C the frictional regime of an average granite sample changes from velocity weakening to velocity strengthening. As well, modeling confirms that postseismic motion at this depth can be regarded as deep afterslip or viscous shear zone creep along a discrete surface as it is not necessary to model a broad region of deformation to achieve a valid result. Blanpied et al (1991) built on these results by including the effect of fluid pressure and concluded that velocity weakening exists at a temperature range of 100°-350°C for a similar granite sample. Even though much of the stress drop during an earthquake occurs in the velocity-weakening zone, the overlying velocity-strengthening zone will also ex-perience displacement as the two are coupled together. The amount of displace-ment experienced at the surface (velocity-strengthening zone) can be likened 18 to a spring whose stiffness controls how much coseismic slip from the velocity-weakening zone is absorbed and how much is transmitted to the surface (Marone et al, 1991). If the spring is very stiff, the coseismic slip at the surface will be almost equal to that at the base of the velocity-strengthening zone. Conversely, in a very weak spring almost no coseismic slip will be experienced at the surface. In either case afterslip can be expected with sufficient time to make up the differ-ence between coseismic slip at the surface and coseismic slip at the base of the velocity-strengthening zone. This is, of course, assuming no plastic deformation of rock surrounding the fault. In their afterslip model, Marone et al (1991) calculated the stiffness of the velocity-strengthening zone as the ratio of the shear modulus to this layer's thick-ness. The means then that a thinner velocity-strengthening layer translates to a higher stiffness that will mean that most of the surface displacement will be the result of coseismic slip rather than afterslip. 4.2 Viscoelastic Relaxation In the context of this study viscoelastic relaxation refers to the time-dependent response of a viscoelastic layer in the lower crust or upper mantle to the stresses induced in it by a shallow earthquake in an overriding elastic layer. A viscoelas-tic material responds elastically on short time scales and viscously on long time scales. In a perfectly elastic earth, any elastic response on the surface would occur almost immediately, while in an earth with a viscoelastic layer there is an 19 initial elastic response followed by transient deformation generally governed by the viscosity, thickness, and depth of the viscoelastic layer (Elasser, 1969; Nur andMavko, 1974). The basic idea of delayed deformation following an increase in stress is similar to that with fault afterslip. There are a few different types of linear viscoelastic rheologies including Maxwell, standard linear solid, and Burgers Body. Maxwell rheology followed by the standard linear solid rheology are the most commonly used in crustal deformation modeling (Hetland and Hager, 2005). A Maxwell rheology is what is used in this study, and its fundamental rheological law relating strain rate (de/dt), stress (cr), and rate of change of stress (dcr/dt) in a one-dimensional form is (Turcotte and Schubert, 2002): — = 1. (A A\ dt ~ 2na + E dt where rj is viscosity and E is Young's modulus. The other linear rheologies differ in that they have either a yield strength (standard linear solid) or two characteristic viscosities (Burgers Body). In non-linear rheologies, which are not considered here, viscosity depends on the dif-ferential stress. The viscoelastic layer, or that which is most capable of ductile deformation, is typically taken to be the asthenosphere (Elasser, 1969), and early models fo-cused only on an elastic layer (upper crust) overlying a viscoelastic half-space 20 (e.g. Nur and Mavko, 1974). Subsequent models have included a viscoelastic layer positioned between two other elastic layers, for example, a viscoelastic lower crust situated between an elastic upper crust and an elastic upper man-tle. Rather than having a homogenous viscosity some studies have modeled a depth-varying viscosity (e.g. Pollitz et al, 2000; Freed et al, 2006). The general method of solution of the viscoelastic response to the stresses in-duced by an earthquake are to first find the static displacements and. stresses in an entirely elastic half-space. The correspondence principle (Fung, 1965) and the Laplace transform are then used to find the temporal response to these stresses (Nur and Mavko, 191 A). Use of the correspondence principle allows results from linear elasticity to be transformed to equivalent results for a Maxwell material through use of the properties of the Fourier transform. 4.3 Other Postseismic Processes Other possible postseismic deformation mechanisms include shallow crustal inelas-tic deformation (e.g. folding of weak sediments) and poroelastic rebound which refers to a reorganization of crustal fluids due to a stress regime change caused by an earthquake. Anelastic deformation can be modeled as afterslip along a plane with afterslip acting as a proxy for such deformation (Donnellan and Lyzenga, 1998). It has been shown that poroelastic deformation, where it exists, tends to counteract the vertical postseismic effects of afterslip by inducing a component of subsidence in regions of positive coseismic stress change (compression) and 21 uplift in regions of negative stress change (extension) (Peltier et al, 1996; Peltzer et al, 1998) where the opposite directions of motion would be expected in the near field from afterslip. A much smaller amount of horizontal motion can be expected to result and is in the direction from tensional to compressional coseis-mic stress change. Peltzer et al (1996, 1998) were able to simulate this effect by lowering Poisson's ratio by 0.03 to represent the change between undrained (coseismic) and drained (postseismic) crustal material. Poisson's ratio, v, represents the compressibility of a material with the (at least theoretical) range of values being 0-0.5 with increasing v corresponding to decreasing compressibility. For a uniaxial stress in the 1-direction, Poisson's ratio is defined as: 62 = 6 3 = -vei (4.5) with the e terms being normal strain. To see why the range of Poisson's ratio is 0-0.5, consider the dilation, A, of a volume: A = e 1 + 6 2 + e3 = ei(l-2v) (4.6) For a totally compressible material, v must equal zero so that dilation is rep-resented entirely by the strain in the 1-direction. Alternatively, v must equal 0.5 22 for the dilational to equal 0, an incompressible material. 23 Chapter 5 Other Postseismic Deformation Studies There are a few other examples of GPS postseismic deformation studies of con-tinental thrust earthquakes (similar to the Racha earthquake). Consideration of these results helps put those of this study into perspective. Postseismic deforma-tion studies of other types of earthquakes are also mentioned below in order to demonstrate the types of results typically found. Continental Thrust Earthquakes Chi-Chi, Taiwan An event with very similar parameters to the Racha earthquake (except magni-tude) is the September 21,1999 Chi-Chi, Taiwan earthquake (Mw=7.6). It had a hypocenter depth of 8 km, dip of 30°, and rake of S5°(Chang et al, 2000). The 24 primary process responsible for the first three months of postseismic deforma-tion was found to be fault afterslip with the majority of afterslip occurring where coseismic slip was minimal (Hsu et al, 2002). The afterslip during the first 3 months after the earthquake is equivalent to about 7% of the coseismic moment (Hsu et al, 2002), and the first 15 months constitute about 16% with this equaling 2.3 times the aftershock moment (after-shock moment = 43% of afterslip moment) (Yu et al, 2003). That is, the first 3 months make up about 44% of the afterslip of the first 15 months. This highlights that afterslip starts soon after the earthquake and decays with time. It has also been shown that a combination of afterslip and viscoelastic mod-els can explain GPS measurements in the first few months after the earthquake with afterslip being dominant around the rupture zone and viscoelastic relax-ation away from this zone (Sheu and Shieh, 2004). A lower-crustal viscosity of 5 x 1017 Pa s (determined from laboratory measurements) and the long-term strain rate were found to best fit GPS data from far-field sites. Preseismic GPS data in Taiwan have also been modeled showing that a viscoelastic asthenosphere below 27-39 km depth with a viscosity of 0.5 - 4 x 1019 Pa s best fits the data (Johnson et al, 2005). Northridge, California The January 17,1994 Northridge earthquake (Mw=6.7) had a hypocentral depth of 20 km. In their two-year postseismic deformation study, Donnellan and Lyzenga (1998) found that most postseismic deformation occurred within the 25 first year following the earthquake and that fault afterslip was the dominant pro-cess, at least for the first few years. Their modeled two-year afterslip is equal to 22% of the coseismic moment, and of this afterslip only 10% could be accounted for with aftershocks. Only sites closer to the rupture area were found to show a clear non-linear postseismic trend (Donnellan and Lyzenga, 1998). Loma Prieta, California The October 17, 1989 Mw=6.9 Loma Prieta earthquake occurred on a buried 70°-dipping fault (Dietz and Ellsworth, 1990). Shallow afterslip on the coseis-mic thrust as well as on an adjacent thrust was found to be the dominant post-seismic process and was active for about 5 years following the quake (Segall et al, 2000). The moment released by this afterslip was equivalent to about 10% of the coseismic moment. Burgmann et al (1997) also found that afterslip and not viscoelastic relaxation was responsible for the first 5 years of postseismic geodetic observations and accounted for several orders of magnitude more mo-ment release than that of aftershocks. GPS site velocities were accelerated over background levels only within 20 km of the earthquake rupture. Non-Thrust Continental Earthquakes The strike-slip August 17, 1999 Mw=7.5 Izmit, Turkey earthquake exhibited afterslip in the immediate 87-day postseismic time period equal to about 20% of the coseismic moment (Burgmann et al, 2002). Of this, only 5% of the moment could be attributed to aftershocks (Reilinger et al, 2000). 26 The first two years of postseismic deformation following the 2002 strike-slip Mw=7.9 Denali, Alaska earthquake have been modeled simultaneously with a combination of deformation mechanisms using 41 GPS sites (Freed et al, 2006). Far-field GPS sites were best explained with a viscoelastic upper mantle (3x 1018 Pa s). Shallow and deep afterslip were better fit to nearer GPS sites. Finally, poroelastic rebound was seen to be active in the immediate area around the fault rupture. Early postseismic deformation following the 1992 strike-slip Mw=7.3 Lan-ders, California earthquake had been niodeled with viscoelastic relaxation of the upper mantle (8 ± 4 x 1018 Pa s) (Pollitz et al, 2000), a viscoelastic upper mantle of 1017 Pa s (Yu et al, 1996), and a viscoelastic lower crust of 1018 Pa s (Deng et al, 1998). It has since been shown that these low viscosities in the lower-crust or upper-mantle do not fit later GPS measurements well as they predict defor-mation rates that are much too large. Deep afterslip (Savage and Svarc, 1997) a combination of deep afterslip and poroelastic rebound (Peltzer et al, 1998; Fialko, 2004), or a non-linear (power-law) viscoelastic lower or upper mantle (Freed and Burgmann, 2004) are now considered to be mechanisms more likely responsible for the postseismic deformation of Landers. Five months of postseismic deformation of the 1999 strike-slip Mw=7.1 Hec-tor Mine, California earthquake have been modeled with deep afterslip using an elastic half-space (Owen et al, 2002). The first three years of postseismic defor-mation were found to be best modeled with a non-linear (power-law) viscoelastic upper mantle. 27 Chapter 6 GPS Data The data used to constrain the models in this study are Global Positioning Sys-tem (GPS) position measurements. The technique of using GPS satellites for millimeter-scale geodetic measurements had its beginnings around the start of the 1980s (Dixon, 1991). Triangulation with at least three GPS satellites can de-termine position of a receiver on the ground with an accuracy to about one meter. This is achieved my measuring the time it takes for a signal from the satellite to reach the receiver and by taking into account satellite and receiver clock errors and well as atmospheric effects. Finer scale precision is obtained using relative positioning to the global network of fixed stations of the International GNSS Service (IGS) that are spread about a few hundred kilometers apart from each other. (GNSS stands for Global Navigation Satellite System.) 28 6.1 Stations and Time Points There were a total of eight GPS sites in and around the Racha rupture area (see Figure 6.1) that were measured at some or all of the following time points: 1991.56301, 1994.76575, 1996.73907, 1998.71096, and 2000.75819. The first time point, 1991.56301, is July 23, 1991, 2 months and 24 days after the Racha earthquake. This is a relatively small amount of data when compared with other postseismic deformation studies. Table 6.1 highlights this by showing the study duration, number of GPS stations, and data points per station for several other studies. STUDY STUDY NUM. GPS AVG. OBS. DURATION STATIONS PER STATION Chi-Chi (Yu et al, 2003) 15 months 80 7 Denali (Freed et al, 2006) 2 years 41 (daily) Hector Mine (Owen et al, 2002) 5 months 32 20 Izmit (Burgmann et al, 2002) 87 days 47 27 Loma Prieta (Burgmann et al, 1997) 5 years 50 14 Loma Prieta (Segalletal, 2000) 8.5 years 62 18 Northridge (Donnellan and Lyzenga, 1998) 2 years 8 5 Table 6.1: Study duration, number of GPS stations, and aver-age number of observations of other postseismic deformation studies. 29 4 2 E 4 3 E 4 4 E 4 5 E Figure 6.1: Location of GPS sites (triangles). Solid rectangle indicates the mod-eled coseismic slip plane and star the epicenter of Tan (2005), and dashed plane indicates rectangle of afterslip plane. 6.2 Acquisition The GPS measurements were obtained during surveys conducted by Rob Reilinger of MIT along with associates from MIT, Indiana University, and in Russia and Armenia (Reilinger et al, 1997) which lasted for about two weeks around the time of each reported measurement. The stations were observed for about 8-12 hours per day for the 1991 measurement and for 10-24 hours per day for subse-quent measurements. 30 6.3 Analysis The GPS data were analyzed by Simon McClusky of MIT using the GAMIT/ GLOBK software (King and Bock, 2004) in two steps. The first was to estimate station coordinates, satellite orbital parameters, atmospheric delay corrections, and earth orientation parameters for each observation session using GPS pseu-doranges and phases. All parameter estimates were loosely constrained. The second step was to combine parameter estimates and covariances from all ses-sions and apply International Terrestrial Reference Frame (ITRF) position and velocity constraints using 14 global IGS core sites (Reilinger et al, 1997). 6.4 Position Measurements Appendix A shows the GPS site position measurements in the ITRF2000 ref-erence frame. This reference frame is defined by a no net rotation condition applied to approximately 200 base station GPS sites on the earth's surface that allow this reference frame to be free from any plate tectonic motion model (Al-tamimi et al, 2002). This model is further constrained by Very Long Baseline Interferometry (VLBI) and Satellite Laser Ranging (SLR) instrumentation. Corrections to Sites SACC and ZELB/ZECK In an effort to make the most of a limited data set, a position value was interpo-lated for that missing at SACC in 1998. Such an interpolation was considered 31 reasonable as no significant deformation signal was expected at this stage in the postseismic period (see Other Postseismic Deformation Studies). This was done by a simple linear interpolation between the position data for 1996 and 2000. The error given to the 1998 data point was an average of the error from the other four time points. This is somewhat conservative as this results in larger errors for the north and up components than if the average of only 1996 and 2000 were considered, while the error of the east component remains about the same. As it turns out, this interpolation became redundant due to how the data were ulti-mately used in modeling, explained below. Stations ZELB and ZECK refer to the same location but with position mea-surements taken over different time periods (1991-1996 and 1998-2000, respec-tively). Their reported values reflect different baselines, that is, to different fixed stations in the IGS network. In order to produce one station (called ZELB), the station with later measurements, ZECK, was combined with ZELB to produce one continuous time-series. This was done by first extrapolating the position of ZELB to 1998 by using the 1991-1996 velocity. The error for this time point was assigned that of ZECK in 1998. The measurement assigned to ZELB at 2000 was then found using the displacement of ZECK between 1998 and 2000. The error from ZECK at 2000 was likewise assigned to that of ZELB at 2000. Conversion from ITRF2000 to Eurasian Reference Frame It was necessary to apply vectors to each station to convert from the ITRF2000 reference frame to a Eurasian reference frame to thereby allow for a more in-32 tuitive consideration of individual vectors. This means that displacement and velocity vectors could now be considered relative to a fixed and nondeforming Eurasian continent (north side of Racha rupture area). Velocities derived by Si-mon McClusky of MIT from the plate motion model of DeMets et al (1994) to perform this conversion are displayed in Appendix B. These velocities were multiplied by the time lapsed since the first position measurement (1991.56301) and subtracted from the ITRF2000 position measurements to produce new val-ues oriented to a Eurasian reference frame. The first measurement did not need to be adjusted since all subsequent use of the data involved displacements and velocities relative to this first measurement. 6.5 Velocity Decay in GPS Time-Series Once in the Eurasian reference frame the data were analyzed to see what kind of accelerated postseismic signal, if any, stood out. This would be manifested by a decaying slope on a displacement-time plot whereby postseismic velocities are initially high and then decay to a slower rate of motion later in the postseismic period. A logarithmic time-series decay can be expected to occur as a result of an afterslip process (Marone et al, 1991), and an exponential decay will result from viscoelastic relaxation (for near-field GPS sites). The priority in trying to detect accelerated postseismic motion was to see if any function with a decaying slope could better fit the GPS position data than could a straight line. In this light, exponential and logarithmic functions were 33 initially tested, but it was decided to use only a logarithmic function as either function provided a similar fit and the logarithmic function was more straight-forward to apply. It was then reasoned that if the residuals of the position data to a logarithmic function were less than those to a straight Une, that this could be in-dicative of accelerated postseismic motion, due either to afterslip or viscoelastic relaxation. The form of the logarithmic function fit to the data is: x(i) = x0 + v(t -tm) + [C + A ln(l + (t- teq/r))]®(teq) (6.1) where v is long-term velocity, C is coseismic offset, A is the amplitude of the logarithmic term, T is the time constant of the logarithmic function, Q(teq) is a Heavy-side step function at the time of the earthquake (teq), and tm is the first time of the time-series. An application developed by Thomas Herring and Simon McClusky of MIT (Tsview) facilitates the fitting of this logarithmic function. It estimates values of the amplitude, A, and long-term velocity, v, and allows the start time of the earthquake, teq, and the time constant, T , to be adjusted. teq was set to 1991.32434 (corresponding to 9:12:48.1 GMT on April 29, 1991), and T was given a typical value of 10 years (McClusky, personal comm.). To test the influence of the choice of T on the fit to data, values between 1 and 100 were used with essentially no difference in fit found. 34 Residuals to the linear and logarithmic fits were calculated using weighted standard deviation (WSTD): WSTD = N'-l (6.2) where w,- is the weight for the /th observation (in this case the inverse of the 1-cr error), N' is the number of non-zero weights, and xw is the weighted mean of the observations: The WSTD is essentially a regular standard deviation with greater empha-sis or weight placed on those observations that are better determined (that have smaller errors). The calculated residuals are shown in Table 6.2. It can be seen here that sites KHUR (N), LESO (N, E, U), SACC (E), and ZELB (E, U) show improve-ments in the fit of their position data to a logarithmic function versus a linear fit. KHUR, LESO, and SACC are located in the immediate earthquake rupture area, and therefore their results are considered a potential indication of accel-erated postseismic deformation. As ZELB is located far from the rupture area on the other side of the Greater Caucasus Mountains and had relatively small 35 displacements, it was not viewed to be as significant. Station Dir. Linear Model Logarithmic Model WSTD Linear Rate WSTD Linear Rate Amplitude N 5.16 0.22 ± 1.07 6.32 -1.08 ± 2.59 5.17 ± 8.96 KHOT E 4.29 2.57 ± 1.09 5.66 0.96 ± 4.42 8.37 ±21.66 U 16.31 2.50 ± 3.54 19.67 -1.93 ± 8.40 18.28 ± 29.86 N 4.19 -2.42 ± 1.02 2.82 -0.36 ± 1.31 -7.74 ±4.18 KHUR E 3.24 0.83 ± 0.92 3.94 -0.28 ±2.18 4.32 ± 7.27 U 14.49 2.55 ± 3.80 14.94 -3.32 ± 7.37 23.71 ± 25.24 N 4.92 -2.36 ± 1.40 0.75 -0.14 ±0.88 -7.06 ± 2.43 LESO E 4.98 2.54 ± 1.43 0.85 3.76 ± 1.01 -4.80 ± 3.45 U 22.83 4.17 ±7.21 3.74 -6.96 ± 4.75 41.59 ± 15.40 N 0.67 -0.48 ± 0.55 0.67 -1.71 ± 1.17 2.76 ± 2.32 NICH E 0.87 1.60 ±0.68 0.87 6.51 ± 1.52 -10.22 ± 2.83 U 3.18 2.13 ±3.25 3.18 -15.56 ± 6.28 53.46 ± 16.25 N 0.31 2.18 ±0.06 1.00 2.09 ± 0.90 0.23 ± 2.08 SACC E 6.63 2.84 ±2.15 1.14 2.81 ± 1.53 0.09 ± 3.53 U 12.46 5.61 ± 3.22 5.07 3.88 ± 4.84 5.84 ± 13.70 N 3.90 0.25 ± 0.77 3.88 1.68 ± 1.61 -1.61 ±8.52 SFRE E 3.62 0.41 ± 0.86 4.35 1.51 ±2.05 7.06 ± 3.75 U 20.55 8.66 ± 4.24 19.79 0.92 ± 8.27 -35.01 ± 29.85 N 1.57 -0.04 ± 0.39 0.76 -0.64 ± 0.85 1.58 ± 1.98 VANI E 3.68 3.04 ± 1.01 0.83 2.51 ±0.99 1.40 ±2 .32 U 26.48 0.33 ± 7.88 4.32 -2.21 ±5.17 8.68 ± 15.28 N 0.24 -0.45 ± 0.09 0.24 -0.52 ± 0.33 -332.69 ± 1.09 ZELB E 4.42 0.50 ± 1.66 0.35 -0.80 ± 0.49 1336 ± 1.51 U 7.42 -0.12 ±3.24 0.85 -2.55 ± 1.56 999.36 ± 4.27 Table 6.2: Residuals (WSTD) and linear rates from linear and logarithmic functions fit to data. Amplitude of logarith-mic function is also included. All values are in mm or mm/yr. 36 6.6 Secular Correction Typically when conducting a postseismic study GPS data are available from be-fore the earthquake, and these allow for an estimation of interseismic, or secular, tectonic motion. This information can be used to calculate constant preseismic velocities that are used to correct the postseismic displacement time series. In so doing the postseismic signal is isolated. Unfortunately in the case of this study, no such preseismic data are available, forcing the search for another method to estimate the secular deformation. Two such methods were explored, using data from late in the time-series and using estimates from a regional tectonic block model. Secular Correction from Time-Series One possibility is to use the latter portion of the time-series as an approxima-tion of the secular tectonic signal, that is, to assume that all or most postseismic deformation was completed by this time. This has the advantage of being able to use data from the same GPS stations, and so if the assumption of completed postseismic deformation is accurate it is potentially nearly as good as having a preseismic data set from these stations. This method of course precludes the possibility of testing for postseismic deformation late in the time-series. As ex-plained in the section on other postseismic deformation studies, afterslip was found to be almost exclusively responsible for observed postseismic deforma-tion in continental thrust earthquakes, and this afterslip was found to last for no 37 longer than five years. In considering the Racha data set for a secular deformation estimation, the two possibilities were to use the 1998-2000 measurements or those from 1996-2000. The latter was adopted as this allowed for more data to be used thereby allowing for a more reliable fit. It also made it possible to use two more sites in the analysis (LESO and VANI) by allowing a line to be fit to two of the three time points (1996.73907,1998.71096, 2000.75819) if only two data points were available. By starting at 1996 to estimate the secular tectonic signal, there is a five-year gap to the earthquake, consistent with findings at Loma Prieta (an earth-quake of the same magnitude as Racha) that all afterslip was finished after five years (Segall et al, 2000). Once the 1996-2000 velocities for each site were cal-culated, they were multiplied by either the 1991-1994 or 1994-1996 time period and then subtracted from either set of respective displacements. This method is hereafter referred to as Time-Series Correction (TSC). Secular Correction from Tectonic Block Model An alternative method to removing the secular tectonic signal from the time-series is to use horizontal velocities (no vertical velocities are available) pre-dicted by a kinematic block model of the Africa-Arabia-Eurasia region (Reilinger et al, 2006). (A kinematic model is one whose motions are valid spatially with-out regard for the processes responsible.) This model is comprised of large blocks representing coherent, undeforming tectonic units based on a best fit to GPS velocities from 1988 to 2005. In this model the blocks are separated by 38 faults slipping at a prescribed rate below a "locking depth" of 15 km. A locking depth refers to the thickness at the top of the crust which is mod-eled as not slipping interseismically. In at least a strike-slip fault situation the breadth of deformation perpendicular to the fault is proportional to the arctan of the locking depth (Savage andBurford, 1973). That is, a thicker locked surface layer means that strain will be manifested in a broader zone around the fault: ydr arctan (6.4) where y is displacement at a surface point, x is perpendicular distance from the fault, and D is locking depth. The boundary along two of the blocks in this model is a straight line passing through the Main Caucasus Thrust (MCT) representing the boundary between the Eurasian and Black Sea and Azerbaijan blocks (see Figure 6.2). As the MCT is actually north of the Racha rupture area, Philip Vernant of MIT was able to adjust this model to have the "MCT" pass through the Racha rupture area, with the site SACC on one side and the other near-field sites KHOT, KHUR, LESO, and SFRE on the other (see Figure 6.1). The temporal range of data used for this block model encapsulates and aver-ages the pre-, co-, and early postseismic deformation of the Racha earthquake. As such, it may overestimate preseismic velocities due to the expectation that coseismic and early postseismic velocities are much larger than preseismic ones. 39 Figure 6.2: Map of part of the tectonic block model used in producing secular velocity estimates (Reilinger et al, 2006). The area of this study as well as the straight line estimate of the Main Caucasus Thrust plate boundary are indicated. 40 Also, the large scale nature of this model does not necessarily lend itself well to accurate deformation predictions on the scale of the Racha GPS network. For example, the model contains one fault in the Caucasus (through the Racha rup-ture area) and does not take into account secular motion along other nearby fault structures, most notably the actual MCT to the north. Finally, there is no ver-tical component in the model, and so one-third of the GPS data are not being corrected by the use of this model, though the secular vertical signal is expected to be small. Though unable to explain all the complexity contained within, the block model may still be a potentially useful, independent approximation of the background tectonic signal that can be used in the absence of actual site-specific preseismic GPS velocities. Along the "MCT" the block model calculated a dip-slip (reverse) slip rate of 3.6 ± 0.3 mm/yr and a right lateral strike-slip rate of 0.7 ± 0.2 mm/yr. The site velocities listed in Table 6.3 were calculated using a 10 km locking depth and 30°north-dipping fault surface. Displacements calculated from the block model velocities were then subtracted from the GPS displacements (relative to the first time point, 1991.56301) to try to isolate the postseismic Racha displacements. This method is hereafter referred to as the Block Model Correction (BMC). 41 Block Model Velocities STATION North (mm/yr) East (mm/yr) KHOT 3.102 0.096 KHUR 1.927 0.114 LESO 2.333 0.059 NICH 6.981 0.650 SACC 3.603 0.398 SFRE 2.524 0.048 VANI 3.874 0.582 ZELB 0.839 0.294 Table 6.3: Block model velocities (mm/yr) used in correcting GPS time-series for secular tectonic displacements. 6.7 Error Combination Once the data had been corrected for secular deformation by the methods de-scribed above the data were differenced to produce displacement (and velocity) values. At this stage, the errors in position at each time point were combined by taking the square root of the sum of the squares of the individual errors. 6.8 Data Used in Modeling The 1991-1994 period was the primary time period on which postseismic de-formation modeling was performed. This is due to it being the earliest period following the Racha event and that most deformation observed following sim-ilar earthquakes takes place in the first few years (see Table 6.1). To test this, 42 GPS displacements from the 1994-1996 period were modeled as well, and data from both time periods are displayed in Appendix C. These data are in the form of displacements (rather than velocities) as there are only the end time points of each period available and all stations share the same exact measurement times. Both the TSC data and BMC data are shown. 43 Chapter 7 Methods 7.1 Afterslip Inversion Modeling Afterslip has been modeled kinematically by using the equations describing dis-placements on the earth's surface due to slip on a fault in an elastic half-space (Okada, 1985). The basic form of these equations comes from Steketee (1958): 1 rr dd" d i = - J ^ A U j M j k - + p 'dd{ ddf vkdL (7.1) which describes the surface displacement di {x\,X2,xf) due to a dislocation Aw/li, §2,13) on a surface, S. dj is the ith component of the displacement at (x\, X2, JC3) due to the y'th direction point force of magnitude F at (^1, § 2 , b^)- <5o is the Kronecker delta, A and p are Lame constants, is the direction cosine nor-mal to dL. This assumes uniform slip on rectangular fault patches and is based 44 on integration of the solution for a point dislocation over the uniform patch. Okada (1985) expressed Equation 7.1 in terms of strike-slip, dip-slip, and tensile dislocations on a fault surface rather than of a point force. These expres-sions are integrated to produce displacements on the earth's surface due to slip on a finite rectangular fault patch (Okada's equations 25-30, 1985). These are the equations that are typically used for elastic dislocation models and earthquake slip inversions. To model afterslip a bounded least-squares inversion approach was used (Stark and Parker, 1995). In the inversion process, Green's functions, G, are calculated to relate the slip on each fault patch, s,, to the displacement at each GPS station, dj in an elastic half-space. The commonly used value of 0.25 was used for Poisson's ratio (e.g. Hreinsdottir et al, 2003; Sheu and Shieh, 2004). Tikhonov regularization was used to enforce a smooth slip distribution. This uses a smoothing parameter,/?, to facilitate a trade-off between the best fit to data and a smoothly varying solution that is physically plausible. This is achieved by minimizing the terms of misfit and roughness: H W - k G ^ - ^ l P + ^ l M l 2 (7.2) where W r W is the data covariance matrix and L is the finite difference approx-imation of the Laplacian operator (V2). Multiplying by W _ 1 weights the misfit term. For example, if a given measurement's error (W) is small, multiplying by W"1 will produce a large value that is relatively more important in the mini-45 mization. ||x|| means (x • xy'2 and is necessary to perform in order to have both terms in Equation 7.2 be scalar quantities. Squaring both terms forces both to be positive values. Afterslip Plane It is assumed that afterslip occurred on the same plane as the coseismic slip. The coseismic fault plane of Tan (2005) was extended 60 km laterally to each side from the hypocenter, up to the surface, and down to 20 km depth along a 30°slope (120 km x 40 km). Figure 7.1 confirms that this fault plane geometry plots directly through the middle of the aftershock activity (Fuenzalida et al, 1997). This not only encloses all of the slip in Tan's 60 km x 24 km coseismic slip model (Figure 7.2) but in the course of testing also captures all slip that the inversion code produced. Tan's fault orientation (Table 3.1) was chosen on the grounds of consistency as this was the only available coseismic slip model and was therefore necessary for the viscoelastic forward-modeling portion of this project. As described subsequently, model results were tested for sensitivity to changes in fault geometry by modifying the fault parameters to match those of some of the other Racha earthquake models and doing separate inversions. Resolution Testing In order to determine which areas of the 120 km x 40 km afterslip fault plane the GPS network is able to resolve slip on, a resolution test was conducted. First, 46 SURFACE -15 -10 - 5 0 5 10 15 20 25 DISTANCE ALONG SURFACE IN DIRECTION OF DIP (KM) Figure 7.1: Profile perpendicular to strike of afterslip fault plane with geometry of Tan (2005) and two months of aftershocks. WNW 13 KM DEPTH ESE 1 KM BELOW SURFACE Figure 7.2: Coseismic fault model, from Tan (2005). Displacement is in meters. Each tile is 4 km square. 47 a "checkerboard" pattern of 20 km x 20 km tiles was constructed with alternat-ing amounts of dip-slip of one or zero meters (Figure 7.3). This was forward-modeled with the resulting GPS displacement vectors then being inverted for slip on the fault plane. If the GPS network is sufficiently dense and evenly distributed around the fault plane, one would expect most of the checkerboard patterned slip to be recovered in the inversion process. The inversion result in Figure 7.3 demonstrates that the GPS network is ca-pable of constraining slip on only the left (WNW) side of the fault plane. As indicated in the figure only the left 60% of this fault was therefore used for slip inversions. Additionally, the upper-left corner of this section of the fault appears incapable of accurately resolving slip with the given array of GPS sites. 7.2 Viscoelastic Forward Modeling Viscoelastic postseismic deformation is modeled using the ViscolD code of Pollitz (1992) which solves for displacement on the surface of the earth due to relaxation of viscoelastic layers at depth resulting from an earthquake in an upper elastic layer. The earthquake displacement field is divided into toroidal and spherical components, the equations of static equilibrium are solved for a spherically symmetric 1-D earth model, and (by using a Maxwell rheology) the correspondence principle is used to find the time-dependent solution. The dis-placement field is calculated using a summation of these normal mode solutions. This model requires that a finite source as well as elastic and viscoelastic 48 FORWARD MODELED INPUT INVERSION ci, OUTPUT SURFACE FORWARD MODELED INPUT INVERSION OUTPUT 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 meters Figure 7.3: Fault plane resolution tests using alternating 1 and 0 meter slip patches on a 120 km x 40 km plane, (a) 60 km x 20 km and (b) 20 km x 20 km sized slip patches were forward modeled (top) and the resulting surface displacement vectors were inverted for slip (bottom) to determine which areas of the modeled fault plane the GPS network can resolve slip on. Inversions were performed using a smoothing, B, of 10 km2/m except where indicated. 49 properties of the spherically symmetric layers be specified. With these parame-ters the code forward models cumulative displacements at prescribed points on the earth's surface, which can then be compared with GPS displacement data. In this model the patterns of displacement are sensitive to the earth model and the earthquake geometry. The viscosity controls the time dependence of the defor-mation (Pollitz, 1992). Three different viscoelastic lower-crust thicknesses were modeled as was a viscoelastic half-space. A viscoelastic lower-crust is often modeled for post-seismic deformation starting at the cutoff depth of seismicity (e.g. Pollitz et al, 2000; Sheu and Shieh, 2004). The May and June, 1991 aftershocks found by Fuenzalida et al (1997) do not extend beyond 15 km depth. Therefore, a 30 km thick viscoelastic lower crust (15-45 km depth) was modeled as were also 20 km and 10 km lower crusts (25-45 km depth and 35-45 km depth respectively). Justification for modeling a viscoelastic mantle half-space (45 km and deeper) is the high Pn attenuation beneath the Greater Caucasus suggesting little to no mantle lithosphere (Al-Lazki et al, 2003). Earth Model Input File The ViscolD code requires that an earth model be specified that includes the density (p), bulk modulus (K), shear modulus (p), and viscosity (77) for each spherical layer. Values of the bulk and shear moduli for the crust were obtained using P and S waves obtained by inversion of P and 5 arrival times from the Racha mainshock and aftershocks (Triep et al, 1995) (See Table 2.1). The rela-50 tions used for calculating these elastic moduli are: (7.3) (7.4) (7.5) where Vp is the P-wave velocity, is the S-wave velocity, K is the bulk mod-ulus, and A and fi and the first and second Lame constants, (JJ. is also known as the shear modulus.) Inserting Equation 7.3, Equation 7.5 becomes: Equation 7.4 is first used to find the shear modulus (JJ). Equation 7.6 is then used to find the bulk modulus (K). Values for the elastic moduli below the crust as well as density values were obtained using the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson, 1981). It was found, however, that varying values of p and K below 80 km depth had essentially no effect on the model results, and so in interest of speeding up the model's run time these moduli were assigned constant values below this depth. The values used in the earth model are listed in Appendix D. 9 4 K = pV2p--n (7.6) 51 Coseismic Slip Model The slip model input to ViscolD was that of Tan (2005) (Figure 7.2 and Table 7.1). As this slip model represents only the mainshock of the four subevents that together lasted 22 seconds, the ViscolD model displacement results were multiplied by 1.637 which is the ratio of the total moment to mainshock-only moment reported by Tan. This can be done due to the linear relation between stress and strain (Equation 4.4). COSEISMIC SLIP MODEL VALUES ALONG KM STRIKE DOWN-DIP 0-4 4-8 8-12 12-16 16-20 20-24 ESE 0.109 0.234 0.259 0.238 0.247 0.201 0.271 0.633 0.869 0.807 0.610 0.366 0.585 1.325 1.922 1.722 1.120 0.586 0.778 1.786 2.737 2.378 1.517 0.781 0.616 1.457 2.172 2.074 1.524 0.869 0.335 0.884 1.344 1.462 1.265 0.829 0.140 0.449 0.671 0.847 0.894 0.711 0.030 0.146 0.209 0.349 0.510 0.500 0.136 0.188 0.129 0.164 0.323 0.392 WNW 0.149 0.161 0.074 0.088 0.174 0.231 Table 7.1: Coseismic slip on 4 km x 4 km fault patches from model of Tan (2005). 52 Chapter 8 Results 8.1 Afterslip 1991-1994 Data The afterslip inversion results from using the 1991-1994 displacements, using the TSC and BMC data and a positivity constraint (no backward slip), are shown in Figure 8.1. Using a smoothing value of 10 km2/m, the reduction in weighted residual sum of squares (WRSS) relative to a model with zero displacements is 89.7489% for the TSC data and 65.6188% for the BMC data. The amount of afterslip moment estimated by each model is 7.585 x 1018 N m (24% of coseis-mic moment of Tan (2005)) and 6.015 x 1018 N m (19% of coseismic moment) respectively. It is notable that both secular correction methods produced similar afterslip distributions with the highest concentration of slip occurring along the surface and along strike to the WNW of the coseismic hypocenter. Both models 53 also show a lobe of afterslip occurring at about 10-16 km depth located about 32 km WNW along strike from the hypocenter. It is also seen in Figure 8.1 through comparison that most of the afterslip occurred immediately adjacent to (and not coincident with) the coseismic slip. 20 KM DEPTH WNW 13 KM DEPTH ESE ESE (C) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Slip (m) 72 km 1 KM BELOW SURFACE 0.5 1.0 1.5 2.0 2.5 Slip <m) WNW SURFACE ESE 0 0.1 0.2 0.3 0.4 0.5 Slip (m) Figure 8.1: 1991-1994 afterslip inversion results using (a) TSC data and (b) B M C data, (c) Coseismic slip model of Tan (2005). Star on each slip plane represents the Racha hypocenter, and each tile is 4 km square. Comparing either afterslip model (a) or (b) with the coseismic slip model (c) shows that the bulk of modeled afterslip occurs away and up-dip from the area of maximum coseismic slip. 54 To evaluate how well the afterslip models fit the data, modeled afterslip was forward-modeled producing displacement vectors that were compared with the data displacement vectors and 1-cr errors (Figure 8.2). Only the nearest sites to the modeled fault (KHOT, KHUR, LESO, and SFRE) were considered as these show the greatest displacements and therefore are the most affected by fault motions. Model vectors from both TSC and B M C data fit the data of these nearest sites at least just within the 1-cr error ellipses if not almost perfectly. Longitude E Longitude E Figure 8.2: Comparison of 1991-1994 data vectors with 1-cr error ellipses (red) with model vectors (blue). Model vectors are produced by forward modeling the afterslip results, (a) is from using TSC data and (b) is from using B M C data. The star represents the Racha epicenter. Comparison to Data without Secular Correction To test what effect either of the secular correction methods (TSC or BMC) have on the data set, the 1991-1994 time period was modeled without any secular correction (Figure 8.3). The basic pattern of slip as with the models of secular corrected data is seen along with some additional slip distributed along the edges. 55 The reduction in W R S S of this model, 51.7435% is less than that of either of the two models based on secular corrected data. The moment of 8.987 x l O 1 8 N m is about 30-40% larger than that of the models with secular correction. -72 km-20 KM DEPTH Displacement for 1991.563 to 1994.7658 relative to SACC <0)42.8~ 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Slip (m) 42.5 43 43.5 Longitude E Figure 8.3: 1991-1994 slip inversion using data with no secular correction. W R S S reduction is 51.7435% and the moment is 8.987 x l O 1 8 N m. To see why the data fitting (WRSS reduction) is significantly worse with the non-secular corrected (NC) data compared to the T S C and B M C data, the data and model vectors of all available sites were plotted for comparison (Figure 8.4). Here it is seen that while all data sets fit the GPS sites close to the fault plane (KHOT, K H U R , L E S O , SFRE) reasonably well , it is primarily the more distant sites ( N I C H , V A N I , Z E L B ) that are fit much less well by N C and to a lesser extent by B M C . (Site N I C H is absent from the T S C data due to a lack of data points to establish a 1996-2000 velocity for this site.) Another way of considering this is that the displacement magnitudes of these distant sites in the B M C and N C data are much larger than in the T S C data. The modeled vectors of the distant sites from all three data sets have about the same magnitude. 56 (a) 43.5 Z 43 g •o 3 J42.S 42 2EL6 ( b) ZELB 43.51 20 M M VANI5 2 43 1 rj 20MM / H D R S ^ KHOT 41.5 42 42.5 43 43.5 44 44.5 Longitude E 42 1 41,5*-VANI NICH 42 42.5 43 43.5 44 44.5 45 Longitude E (c) TzaB 43.6P 43.2 T3 1 4 2 . 8 is 42.4 42 20 M M VANI KHOT n i c h 42 42,5 43 43.5 44 44.5 45 Longitude E Figure 8.4: 1991-1994 afterslip model and data vectors for all available sites for (a) TSC data (b) B M C data and (c) non-corrected (NC) data. All Data Sets with Same GPS Sites The TSC data are lacking the relatively distant site NICH due to the lack of data at this site to determine a 1996-2000 velocity. In the previous section the TSC, B M C , and NC data sets were compared. In order to see what effect NICH has on how well the model fits the B M C and NC data sets, these other data sets were modeled without NICH and the inversion results compared (Figure 57 8.5). By removing NICH, the B M C model improved its WRSS reduction from 65.6188% to 72.3125%, and the NC model improved slightly from 51.7435% to 53.7368%. The B M C and NC slip models appear almost identical and both of these are very similar to the TSC model. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Slip (m) Slip (m) - 72 km -0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Slip (m) Figure 8.5: 1991-1994 slip inversions using the same GPS sites for all data sets (without site NICH). (a) TSC model (b) B M C model (c) NC model. Sensitivity to Fault Parameters As stated in Section 7.1, the fault geometry used for inversion modeling assumed the values of those found by Tan (2005). In order to ascertain how sensitive the 58 inversion results are to choice of fault geometry, the WRSS reduction found us-ing the TSC data for the 1991-1994 period was compared with WRSS reduction values resulting from using the fault geometries of the other models listed in Table 3.1. A B of 10 km2/m was used throughout. This was performed in a step-wise process, and Table 8.1 shows how the WRSS reduction was affected by incremental changes in fault parameters to adapt from those of Tan (2005) to either (a) Fuenzalida et al (1997), (b) Triep et al (1995), or (c) the Harvard CMT solution. The published hypocenter depths of these models were not used but rather bracketed, as shown in the last two lines of each sub-table. This was necessary to ensure that the fault plane reached the earth's surface while retain-ing the 4 km2 fault tile dimension. The last line of each sub-table (bold-faced) represents the greater diversion from the fault geometry of Tan and is therefore a conservative (over-) estimate of fault geometry difference. As seen in Table 8.1, the effects of the differences in reported fault geometries on WRSS reduction are minor. The actual difference in fault plane location and orientation is relatively small as Figure 8.6 shows. 59 From Model of Tan to Fuenzalida Change WRSS Red. (%) % Change (no change) 89.65% 0.0 lat/lon to 42.42°N/43.69°E 89.68% +0.033% and strike to 286.7° 89.69% +0.050% and dip to 29°, depth to 5.8177 km 89.65% -0.0046% and dip to 29°, depth to 3.8785 km 89.48% -0.19% (a) From Model of Tan to Triep Change WRSS Red. (%) % Change (no change) 89.65% 0.0 lat/lon to 42.4238°N/43.6643°E 89.60% -0.053% and strike to 292.9° 89.46% -0.21% and dip to 24°, depth to 4.8808 km 89.59% -0.069% and dip to 24°, depth to 3.2539 km 88.12% -1.70% (b) From Model of Tan to Harvard CMT Change WRSS Red. (%) % Change (no change) 89.65% 0.0 lat/lon to 42.60°N/43.61°E 3.25% -96.38% and strike to 288° 3.08% -96.56% and dip to 39°, depth to 22.6555 km 90.89% +1.38% (c) Table 8.1: Sensitivity of afterslip results to choice of fault plane as measured by changes in WRSS reduction explained by inversion models. Inversions are based on TSC data and a B (smoothing) of 10 km2/m. Fault geometries begin with the parameters of Tan (2005) and are adapted incrementally to the models of (a) Fuenzalida et al (1997), (b) Triep et al (1995), and (c) Harvard CMT. The reduction in WRSS obtained using the different fault geometries varies 60 / relatively little. To compare these changes with that of the relatively arbitrary task of choosing a smoothing value (fi), the original WRSS reduction of 89.6489% with 8 = 10 km2/m changes to 89.8284% with 8 = 9 km2/m or 89.4820% with 8 = 11 km2/m. These changes correspond to +0.2002% and -0.1862% from 89.6489% respectively. 42.5° N 43.0° E 43.5° E 44.0° E Figure 8.6: Surface outline of the fault plane using the fault geometries of Tan, Fuenzalida, Triep, and Harvard CMT. Stars are the Racha epicenter according to the different groups. It is not enough to compare just error reduction, rather it is also important to know how much the actual slip solution varies among the different possible fault geometries. This is shown in Figure 8.7 where it can be seen that the slip 61 -distributions resulting from the use of the geometries of Fuenzalida et al (1997) and Triep et al (1995) show the same basic slip pattern as that of Tan (2005). The main difference to be noted is that the modeled slip tends to not extend as far down the plane, especially with the geometry of Triep. The moment of the modeled afterslip is relatively similar as well with 26.0918%, 29.9303%, and 35.451% of the coseismic slip of Tan for the geometries of Fuenzalida, Triep, and Tan respectively. SURFACE SURFACE Figure 8.7: Patterns of modeled afterslip using different published fault geome-tries. Star indicates the Racha hypocenter. (a) Tan (2005), (b) Fuenzalida et al (1997), (c) Triep et al (1995), (d) Harvard CMT. 62 Determination of Best Smoothing Value p As indicated in Equation 7.2, smoothing or B acts in the inversion as a weight on how significant the roughness term should be. That is, when minimizing the misfit and roughness, the smoothing determines how important achieving (a smoothly varying result should be. As smoothing is increased, results become more evenly distributed and therefore more likely to resemble what happens in nature, though this comes at the cost of not explaining the data as well as a model with less smoothing. There is no obvious way of knowing which smoothing value gives the best trade-off, and choosing one seems to be something of an art. In all of the inversion results described above, a B of 10 km2/m was used. This was initially based primarily on trial and error with 10 km2/m seeming to be on the threshold of B values that produced slip distributions without scattered patches of unrealistic appearing slip. That is, smaller values of B gave slip dis-tributions with faint, disconnected areas of slip that appeared unrealistic in the sense that they were disjointed and erratically distributed. In order to ensure that 10 km2/m for 8 was a reasonable choice and to help determine how reliable the patch of deep slip shown in Figure 8.1 may be, two different methods of choosing an optimal smoothing were tested. The first was simply make a "Tikhonov plot" which is the model's weighted residual sum of squares (WRSS) versus roughness (m2/km4) and to find the point on the plot where the curvature between the two is the greatest and then select the corresponding B. This B value is taken as an approximate compromise (e.g. Hansen, 1998). The idea is that changing B from this point would cause either of 63 the following undesirable situations to arise: 1) allowing the misfit to steadily in-crease while producing only marginal decreases in roughness or 2) allowing the roughness to increase substantially while reaping only minimal improvements in misfit. The second method tried was the cross validation sum of squares (CVSS) (Matthews and Segall, 1993). In this method the optimal B is chosen based on that which allows the most efficient prediction of missing stations. Each station in the GPS network is removed one at a time and inversions are performed with the remaining stations for a variety of B values. The resulting inversions are then forward-modeled for the missing site to find its residuals between the model and data. Once this is done for each site all the squared residuals are summed up according to B value, and it is this sum of residuals that is used to select the optimal B. The B that produces the minimal sum of squared residuals is that considered to be the optimal smoothing value. The results of these two tests for both types of secular correction (TSC and BMC) are shown in Figure 8.8. The range of optimal B values for the TSC data is 2-4 km2/m and that for the BMC data is 3-7.5 km2/m. The differences in amount of WRSS reduction and slip moment between using these B values and 10 km2/m are shown in Table 8.2. That the deep patch of slip still exists at B = 10, a higher than "optimal" value which acts to blur the spatial distribution of slip, argues for the believability of this slip. 64 Comparison of Different Smoothing Values WRSS Moment % Coseis. (km2/m) Red. (Nm) Moment Data with Time- 2 91.08% 6.96±0.56 x 101" 21.78% Series Correction 4 90.74% 7.11±0.30x 1018 22.26% (TSC) 10 89.65% 7.58±0.13x 1018 23.74% Data with Block 3 66.85% 5.12±0.36x 101* 16.04% Model Correction 7.5 66.08% 5.74±0.16x 1018 17.97% (BMC) 10 65.62% 6.02±0.13x 1018 18.83% Table 8.2: Comparison of WRSS reduction and moment re-sulting from models with different B values. The first two 0 values for each data set are those determined by the meth-ods of analyzing the trade-off between misfit and roughness and of cross validation sum of squares (CVSS). The third 6 value, a more conservative 10 km2/m, is that which was used in many of the tests presented in the results. This value had been chosen due to it allowing for a more reasonable-appearing slip distribution. Performing the CVSS test also indicated that it is only site KHUR that is responsible for the deep patch of slip shown in Figure 8.1. KHUR happens to lie almost directly above this deeper modeled area of slip, and as seen in Figure 8.2 it is one of the better determined sites in terms of signal to noise (or displacement to error). To see what the difference in inversion results looks like from using these different B values, Figure 8.9 displays the inversions using the "optimal" B values as well as a B of 10 km2/m. 65 % 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0 Q % 2 4 6 8 10 12" 14 16 ' 18 Roughness (m?/km4) . fi (Nm2/m) Figure 8.8: Determination of smoothing parameter, B, for TSC data (a,b) and BMC data (c,d). The first method used, a "Tikhonov plot", determines where the change between the weighted residual sum of squares (WRSS) and rough-ness is greatest (a,c) and then use the associated B. The second method is cross validation sum of squares (CVSS) where B is chosen based on which value pro-duces the smallest sum of residuals. Between the two methods, the optimal value of B ranges from 2 to 4 km2/m for the TSC data and 3 to 7.5 km2/m for the BMC data. 66 T S C D A T A B M C D A T A Figure 8.9: Afterslip inversions resulting from use of different smoothing values, B, for TSC and BMC data. The first two smoothing values displayed for each data type are those corresponding to the methods of WRSS versus roughness and of CVSS (see Figure 8.8). Model Sensitivity to v To test of the sensitivity of afterslip inversion results to choice of Poisson's ratio, v, values 0.03 greater and smaller were used in inversions with TSC and BMC data. Table 8.3 shows these results. 67 Sensitivity to Poisson's Ratio Poisson's T S C Data B M C Data ~ Ratio (v) WRSS Red. Moment (Nm) WRSS Red. Moment (Nm) 028 89.707% 7.765 x 1018 65.905% 6.195 x 1018 0.25 89.649% 7.585 x 1018 65.619% 6.015 x 1018 0.22 89.584% 7.415 x 1018 65.336% 5.852 x 1018 Table 8.3: Test for the sensitivity to Poisson's ratio, v, by in-creasing and decreasing v from the value of 0.25 otherwise used. Differences in WRSS reduction and moment are indi-cated. Experiment with Slip Constraints Up to now all afterslip inversions were performed with a dip-slip-only constraint. That is, only dip-slip motion was allowed along the fault as this is consistent with Tan's (2005) coseismic slip model. As there is a small component of right-lateral strike-slip motion along the Greater Caucasus mountains (e.g. Triep et al, 1995; Fuenzalida et al, 1997), an experiment was tried allowing right-lateral as well as dip-slip motion along the fault. The result is that for the TSC and BMC data the WRSS reduction increased only modestly (89.649% to 90.750% and 65.619% to 67.130% respectively). A similar dip-slip pattern as before resulted, but all of the strike-slip component was placed in a far, deep corner of the fault plane (bottom northwest corner) (see Figure 8.10). 68 DIP-SLIP 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 (m) 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Cm) 0 0.05 0.1 0.15 0.2 0 25 0.3 0.35 0.4 0 0.02 0.04 0.06 0.08 0.1 0.12 (m) ( m ) Figure 8.10: 1991-1994 afterslip inversions with dip-slip and strike-slip compo-nents using (a) TSC data and (b) B M C data. 1994-1996 Displacements While it was not expected to be able to detect afterslip much beyond three years following the Racha event (see Other Postseismic Deformation Studies), the 1994-1996 time period was modeled for afterslip to ensure that this was a valid assumption. Figure 8.11 shows these results for the two methods of secular cor-rection. Although similar through the center of the fault plane, these afterslip models do not explain the data very well. Again using a B of 10 km 2/m the WRSS reduction for the TSC and B M C methods were 22.0902% and 24.9727% 69 respectively. The moments due to each correction type are 2.714 x 10 1 8 N m and 1.247 x 10 1 8 N m respectively. As the 1994-1996 data fit was not good, afterslip was not modeled in later time periods. 72 km 20 KM DEPTH (b)| a . D z o a WNW SURFACE 0 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 Slip (m) 0 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 Slip (m) Figure 8.11: 1994-1996 afterslip inversion results (a) using data TSC data and (b) using B M C data. Star on each slip plane represents the Racha hypocenter, and each tile is 4 km square. 8.2 Viscoelastic Relaxation Modeling 1991-1994 Displacements Viscoelastic Lower Crust The spherically symmetric viscoelastic forward-modeling code, ViscolD (Pol-litz, 1992), was used to forward model the Racha mainshock model of Tan (2005). Three different viscoelastic lower crust thicknesses were used: 10 km, 20 km, and 30 km. This thickness is measured upward from the base of the 45 km thick crust. A range of viscosities (1 x 10 1 7 to 1 x 10 2 0 Pa sec) was modeled. 70 The weighted residual sum of squares (WRSS) was calculated to assess the fit of the viscoelastic models to the data corrected with both the TSC and B M C meth-ods. Contour plots were made of the WRSS against thickness of viscoelastic lower-crust and logarithm of viscosity (Figure 8.12). Figure 8.12: Sensitivity of WRSS to lower-crust viscoelastic layer thickness and viscosity. The best model with the TSC data (a) has a viscoelastic lower crust of 10 km thickness and a viscosity of about 1.0 x 10 1 7 Pa sec. The best model (star) with B M C data (b) has a viscoelastic lower-crustal thickness of 10-18 km and a viscosity of about 2.5 - 6 x 10 1 7 Pa sec. The viscoelastic lower-crust models show the best fit to the TSC data with a layer thickness of 10 km and a viscosity of about 1 x l O 1 7 Pa sec with a WRSS reduction of about 19%. The viscoelastic lower-crust models have a minimum WRSS for the B M C data with a layer thickness of about 15-20 km and a viscos-ity of 5 x l O 1 7 to 1 x l O 1 8 Pa sec (Figure 8.12b). This corresponds to a WRSS reduction of about 21%. Neither set of GPS data is well modeled as viscoelastic relaxation of a lower-crustal layer. The vectors of these best models are compared against data in Figure 8.13. 18 18.5 19 19.5 LOG VISCOSITY (PA SEC) 18 18.5 19 19.5 LOG VISCOSITY (PA SEC) 20 71 42.4 42-6 42-8 43 43.2 43.4 43.6 43 8 44 42.4 42.6 42.8 43 43.2 43.4 43.6 43.8 44 LONGITUDE E LONGITUDE E Figure 8.13: Displacement vectors from best fitting ViscolD models and 1991-1994 GPS site displacements, (a) TSC data with a viscoelastic lower crust at 35-45 km with a viscosity of 1 x l O 1 7 Pa sec. (b) B M C data with a viscoelastic lower crust at 25-45 km with a viscosity of 9 x l O 1 7 Pa sec. Star is the Racha epicenter. Viscoelastic Mantle Half-Space The results of modeling a viscoelastic mantle half-space are indicated in Figure 8.14. The best fitting model for both TSC and B M C data sets has a viscosity of about 2 x l O 1 8 Pa sec corresponding to a WRSS reduction of about 5-8%. 1996-2000 Displacements As stated previously, the TSC data don't allow the possibility of testing for post-seismic deformation during the 1996-2000 time period. The block model cor-rected (BMC) data, however, are not subject to this restriction. Lower-crustal viscoelastic relaxation was modeled for an array of viscosity values. A layer thickness of 20 km (25-45 km depth) was used as this is the thickness that best fit the data for the 1991-1994 period. The best fitting 1991-72 LOG VISCOSITY (PA SEC) LOG VISCOSITY (PA SEC) Figure 8.14: WRSS versus logarithm of viscosity for forward-modeled vis-coelastic mantle half-space models compared with (a) TSC data and (b) BMC data. 1994 BMC model has a viscosity of 9 x 1017 Pa sec (and a 20 km thickness). The WRSS reduction with this viscosity and thickness for the 1996-2000 data is less than 0%, a very poor fit. Figure 8.15 shows that the best fit to the 1996-2000 data produces about a 25% WRSS reduction with a viscosity of about 1019 Pa sec. 73 LOG VISCOSITY (PA SEC) LOG VISCOSITY (PA SEC) Figure 8.15: Comparison between BMC data and models of a 20 km thick vis-coelastic lower crust (25-45 km depth), (a) WRSS versus logarithm of viscosity, (b) Corresponding WRSS reduction versus logarithm of viscosity. 74 Chapter 9 Discussion 9.1 Secular Correction of GPS Data Since no GPS sites were established in the region surrounding the Racha earth-quake prior to its occurrence, pre-earthquake site velocities were unavailable. This spawned the search for other correction methods, and two were used through-out this study. The first method is the time-series correction (TSC) that assumes that all or most postseismic deformation was over within the first five years and which used the 1996-2000 period of the GPS data set to estimate the secular deformation rate. The second method is the block model correction (BMC) that used the results of a tectonic block model to estimate the secular deformation. Both correction methods have their limitations. The TSC method may indeed introduce an overcompensation in that there could still be a postseismic signal contained in the 1996-2000 time-series that gets subtracted from the first period 75 of the time-series. The BMC method is based on a model that is spatially coarse relative to the physical dimensions of this study. It also contains only horizontal motions, so that the vertical component does not get corrected for secular dis-placement by this method. Finally, the BMC method is based on a model that includes the entire seismic cycle and not just the preseismic period. This is a problem because deformation rates and patterns around a fault can change be-tween earthquakes (Savage and Prescott, 1978). Due to these uncertainties in the two correction methods, both were employed throughout with the outcome that similar results were obtained. That both methods obtain similar results lends confidence to these findings that would not be achieved by either method alone. To otherwise compare how significant the choice of secular correction method is, the data were modeled for afterslip with no correction for secular deforma-tion. Although the WRSS is much less impressive (52% versus 66% or 90%), the overall pattern of afterslip is still very similar (Figure 8.5). This shows that no matter how the data are secularly corrected, there is little effect on the inferred pattern of afterslip. That the secularly corrected data are much better fit by the model indicates that the secular correction methods employed are valid, perhaps one (TSC) more so than the other (BMC). Figure 8.4 shows why the TSC correction method results in such a better WRSS reduction. While the near-field GPS sites are fit reasonably well by the two secular correction methods and by no correction, it is the far-field sites that are fit best by the TSC method. This shows again that the TSC data is probably the most appropriate or accurate to use and by extension that the assumption of 76 no postseismic deformation after five years is likely reasonable. 9.2 Afterslip and Viscoelastic Relaxation Models 1991-1994 Displacements The 1991-1994 data set was best fit by an afterslip model that explained about 90% of the data. By comparison, the best viscoelastic model produced a WRSS reduction of only 21%. To see what this difference in model fit looks like in term of displacement vectors, a comparison of Figures 8.2a and 8.13b shows that the afterslip model closely mimics the data while the viscoelastic model bears very little actual resemblance: Due to differences in the modeled afterslip moment from using different val-ues of smoothing parameter (0) and due to uncertainties in choosing an appro-priate value for 0, the moments produced by the TSC data are about 22-24% of the coseismic moment (with a 0 range of 2-10 km2/m). (0=2 km2/m results from the CVSS method, and the conservative value of 0=10 km2/m is used through-out most of the modeling and testing.) This afterslip moment release is on the high end though still similar to other studies that found afterslip in the imme-diate postseismic period with values of 10-22% of the coseismic moment being reported (e.g. Yu et al, 2003; Donnellan and Lyzenga,\99%; Segall et al, 2000). This must still be approached with some caution, however, as the start of the GPS time-series was about 2.75 months after the Racha event, and a sizable amount of afterslip may have occurred undetected during this time. For example, 77 in the first three months after the Chi-Chi earthquake about 7% of the coseismic moment was detected in afterslip (Hsu et al, 2002). So, the actual postseismic moment release following the Racha earthquake may be something larger than 22% of the coseismic moment. On the other hand, the afterslip predicted by the elastic half-space model may be an overestimate due to a rock type contrast across the fault zone. The crystalline basement Dzirula Massif that comprises the footwall may be expected to deform less than the sedimentary Gagra-Dzhava zone that makes up the hanging wall due to the rocks of the hanging wall being less dense and weaker than those of the footwall. As most of the GPS stations sit on the hanging wall the inversion is mainly trying to fit these hanging wall displacements. As the model contains no difference in rock type across the fault, the slip inversion may be prescribing too much slip as it is almost exclusively fitting displacements of the weaker hanging wall. (If there were more stations on the footwall, the inversion would be more balanced.) Considering the possibilities in choosing fault parameters based on published studies of the source geometry, the differences in data explanation with all other factors equal are generally smaller than those of the aforementioned 0 range. For example, the range of WRSS reduction from using any of the three available fault geometries based on teleseismic body waves (at 0=10 km2/m) is 88.123-89.649% while that from the a 0 range of 2-10 km2/m (with the Tan fault geom-etry) is 89.649-91.079%. When allowing strike-slip as well as dip-slip afterslip motion, the meager improvement in WRSS reduction does not justify the inclusion of the strike-78 slip component. Also, according to the kinematic model this strike-slip motion would have to occur in a completely different area of the fault plane which is the very corner (northwest) where the GPS network is not able to resolve slip (Figure 7.3). Adjusting Poisson's ratio up or down by 0.03 from 0.25 did not produce a significant effect in data explanation or afterslip moment. Therefore the choice of the typically used value of 0.25 seems to be reasonable. As expected given the relatively small size of the earthquake and local scale of deformation, viscoelastic relaxation of the mantle is not responsible for ob-served early postseismic deformation as has been suggested for other earth-quakes (e.g., Pollitz, 2000). Comparison of Afterslip with Aftershocks A comparison of afterslip with the distribution of aftershocks (Fuenzalida et al, 1997) shows that the two generally do not overlap (Figure 9.1). While these aftershock data were collected in the two months between the Racha event and the first GPS position measurement, they provide a good overall picture of the distribution of aftershocks. One can also assume that afterslip that occurred before the first GPS measurement was in a similar location. It is noteworthy that the deep patch of the afterslip happens to fall in an area with no aftershocks but which is bounded on three sides by significant aftershock activity. This seems to = indicate that this section of the fault deforms aseismically and that it loads the surrounding region. 79 (a) Figure 9.1: (a) Aftershocks over the two months following the Racha event (May 10, 1991 to July 9, 1991) (Fuenzalida et al, 1997). Afterslip fault plane is out-lined. Star represents Racha hypocenter and can be used with fault plane outline to compare aftershocks with 1991-1994 (b) TSC afterslip and (c) BMC afterslip. The afterslip models are from data covering the period July 25, 1991 to October 10, 1994. 80 The 1991-1994 afterslip moment was compared with seismic moment from aftershocks during the same time period over the area of the afterslip fault plane. The aftershocks used were from the Incorporated Research Institutions for Seis-mology (IRIS) which provide a compilation of events from different earthquake catalogs. The total moment of these events summed to 2.4461 x 1016 N m. This is far less than the moment release from afterslip (7.112 x 1018 N m) which shows that aseismic deformation is required by the GPS data. This difference in moment between afterslip and aftershocks is comparable to that found in Loma Prieta, California (Biirgmann et al, 1997). Post-1994 Displacements To test if afterslip continued beyond 1994, the 1994-1996 time period was mod-eled for afterslip with the result that no pattern of afterslip could accurately ex-plain either the TSC or BMC data. As viscoelastic relaxation can manifest itself later in the postseismic period, the 1996-2000 BMC data were modeled with the result that the best viscoelastic model could explain only 25% of the data. Very low viscosities (< 1018 Pa sec) produce models with negative WRSS reduction. That is, at these low viscosi-ties modeling produces displacement vectors that are much larger than the data. According to these results, lower-crust viscosities less than 1018 Pa sec can be ruled out. This can place constraints on theoretical models of continental colli-sion models. For example Pysklywec et al (2002) model a continental collision zone using a lower-crustal rheology of wet quartz which has a viscosity on the 81 order of 10 7 Pa sec. If the results presented here hold for other continental col-lision zones in general, such a low viscosity can be dismissed as unreasonable. Moderately low viscosities (> 1019 Pa sec) cannot necessarily be dismissed be-cause models with such viscosity values yield too little surface deformation to be resolved with campaign-model GPS. 82 Chapter 10 Conclusions It was found that postseismic deformation after the Racha earthquake was due to afterslip. Most afterslip was shallow, though a patch at 10-16 km depth was identified. Slip could be resolved only west of the hypocenter. Though the TSC data, Tan's fault geometry, only dip-slip on the coseismic rupture plane, and a B of 4 km2/m are preferred, these findings hold even when these constraints are relaxed. The final WRSS reduction and afterslip moment over the 1991-1994 time period given the above parameter choices is 91% and 7.112 ± 0.295 x 1018 N m respectively. (The uncertainty in moment results from the uncertainty in slip produced by the inversion code.) This moment corresponds to 22.2± 0.9% of the coseismic slip of Tan (2005), and only 0.34% of this afterslip can be attributed to aftershocks. Viscoelastic models did not fit the GPS data well and afterslip was found to be the dominant mechanism of postseismic deformation in the first three years 83 following the Racha earthquake. However, the existence of a viscoelastic lower crust or upper mantle in the western Greater Caucasus cannot be ruled out as it may simply be that this event was not large enough to excite a noticeable viscoelastic response. 84 Bibliography Alexidze, M.A., Gugunava, G.E., Kiria, D.K., and Chelidze, T.L., A three-dimensional stationary model of the thermal and thermoelastic fields of the Caucasus, Tectonophysics, 227, 191-203, 1993. Al-Lazki, A.L, Seber, D., Sandvol, E., Turkelli, N., Mohamad, R., and Barazangi, M., Tomographic Pn velocity and anisotropy structure beneath the Anatolian plateau (eastern Turkey) and the surrounding regions, Geophys. Res. Lett., 30, 8043-8047, 2003. Altamimi, Z., P. Sillard, and C. Boucher, ITRF2000: A new release of the International Terrestrial Reference Frame for earth science applications, /. Geophys. Res, 107, 2214-2233, 2002: Asanidze, B.Z., Pecherksy, D.M., and Adamia, Sh.A., Results of paleomagnetic studies of Paleozoic rocks in the Caucasus (in Russian), Izv. Acad. Sci. USSR, Phys. Solid Earth, Engl. Transl, 34-48, 1980. Balavadze, B.K., Shengelaya, G.Sh., and Mindel, P.Sh., A gravitational model of the earth crust for the Caucasus and Caspian Sea area, In: Gravitational Models of the Earth Crust and Upper Mantle, emphNarkova Dumka, Kiev. 85 Blanpied, M.L., D.A. Lockner, and J.D. Byerlee, Fault stability inferred from granite sliding experiments at hydrothermal conditions, Geophys. Res. Lett., 75,609-612,1991. Burgmann, R., S. Ergintav, P. Segall, E. Hearn, S. McClusky, R. Reilinger, H. Woith, and J. Zschau, Time-dependent distributed afterslip on and deep below the izmit earthquake rupture, Bull. Seism. Soc. Am., 92 126-137, 2002. Burgmann, R., Segall, P., Lisowski, M., and Svarc, J., Postseismic strain following the 1989 Loma Prieta earthquake from GPS and leveling measurements, J. Geophys. Res., 102, 4933-4955, 1997. Burtman, V.S., Kinematics of the Arabian syntaxis, Geotectonics, 23, 139-146, 1989. Chang, C.H., Y.M. Wu, T.C. Shin, and C.Y. Wang, Relocation of the 1999 Chi-Chi earthquake in Taiwan, Terr. Atmos. Oceanic Sci., 11, 581-590, 2000. Cohen, S.C., A multilayer model of time dependent deformation following an earthquake on a strike slip fault, J. Geophys. Res., 87, 5409-5421, 1982. DeMets, C , R.G. Gordon, D.F. Argus, and S. Stein, Current plate motions, Geophys. J. Int., 101, 425-478, 1990. DeMets, C , Gordon, R.G., Argus, D.F., and Stein, S., Effects of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions, Geophys. Res. Lett., 21,2191-2194, 1994. Deng, J.M., H. Kanamori, and E. Hauksson, Viscoelastic flow in the upper crust after the 1992 Landers, California earthquake, Science, 282, 1689-86 1692,1998. Dietrich, J.H., Modeling of rock friction: 1. experimental results and constitutive equations, J. Geophys. Res., 84, 2161-2168,1979.1ing of rock friction: 1. experimental results and constitutive equations,/. Geophys. Res., 84, 2161-2168, 1979. Dietz, L.D. and Ellsworth, W.L., The October 17, 1989, Loma Prieta, California, earthquake and its aftershocks: Geometry of the sequence from high-resolution locations, Geophys. Res. Lett., 17, 1417-1420, 1990. Dixon, T.H., An introduction to the global positioning system and some geological applications, Rev. of Geophys, 29, 249-276, 1991. Donnellan, A. and G.A. Lyzenga, GPS observations of fault afterslip and upper crustal deformation following the Northridge earthquake, J. Geophys. Res., 103, 21285-21297, 1998. Dotduyev, S.I., Nappe structure of the Greater Caucasus Range, Geotectonics, 20,420-430,1986. Dziewonski, A.M. and Anderson, D.L., Preliminary reference Earth model, Phys. Earth Plan. Interiors, 25, 297-356, 1981. Elasser, W.M., Convection and stress propagation in the upper mantle, in The Application of Modern Physics to the Earth and Planetary Interiors, Wiley-Interscience, 1969. Fialko, Y. Evidence of fluid-filled upper crust from observations of postseismic deformation due to the 1992 Mw 7.3 Landers earthquake, 87 J. Geophys. Res., 109, B08401, doi:10.1029/2004JB002985, 2004. Freed, A.M. and R. Burgmann, Evidence of power-law flow in the Mojave desert mantle, Nature, 430, 548-551, 2004. Freed, A.M., R. Burgmann, E. Calais, J. Freymueller, and S. Hreinsdottir, Implications of deformation following the 2002 Denali, Alaska, earthquake for postseismic relaxation processes and lithospheric rheology, /. Geophys. Res., Ill, B01401, doi:10.1029/2005JB003894, 2006. Fuenzalida, H., Rivera, L., Haessler, H., Legrand, D., Philip, H., Dorbath, L., McCormack, D., Arefiev, S., Langer, C , and Cisternas, A., Seismic source study of the Racha-Dzhava (Georgia) earthquake from aftershocks and broad-band teleseismic body-wave records: an example of active nappe tectonics, Geophys. J. Int., 130, 29-46, 1997. Fung, Y.C., Foundations of Solid Mechanics, (Prentice-Hall, Englewood Cliffs, NJ, p. 425, 1965) Hansen, R C , Rank-Deficient and Discrete Ill-Posed Problems, Society for Industrial and Applied Mathematics, 1998. Hearn, E.H., B.H. Hager, R.E. Reilinger, Viscoelastic deformation from North Anatolian Fault Zone earthquakes and the eastern Mediterranean GPS velocity field, Geophys. Res. Lett., 29, 10.1029/2002GL014889, 2002. Hetland, E.A. and B.H. Hager, Postseismic and interseismic displacements near a strike-slip fault: A two-dimensional theory for general linear viscoelastic rheologies, 7. Geophys. Res., 110,B 10401, 88 doi:10.1029/2005JB003689, 2005. Hreinsdottir, S., Freymueller, J.T., Fletcher, H.J., Larsen, C.F., and Burgmann, R., Coseismic slip distribution of the 2002 Mw7.9 Denali fault earthquake, Alaska, determined from GPS measurements, Geophys. Res. Lett., 30, 1670-1674. Hsu, Y.J., N. Bechor, P. Segall, S.B. Yu, L.C. Kao, and K.F. Ma, Rapid afterslip following the 1999 Chi-Chi, Taiwan earthquake, Geophys. Res. Lett, 29(16), 1974, doi:10.1029/2002GL014967, 2002. Jackson, I., Peterson, M.S., and FitzGerald, J.D., Seismic wave dispersion and attenuation in Aheim dunite: an experimental study, Geophys. J. Int., 108, 517-534,1992. Jackson, J., Partitioning of strike-slip and convergent motion between Eurasia and Arabia in eastern Turkey and the Caucasus, J. Geophys. Res., 97, 12471-12479, 1992. Johnson, K.M., Segall, P., and Yu, S.B., A viscoelastic earthquake cycle model for Taiwan, J. Geophys. Res., 110, 1029/2004JB003516, 2005. Kampfmann, W. and Berckhemer, H., High temperature experiments on the elastic and anelastic behavior of magmatic rocks, Phys. Earth Planet. Inter., 40, 223-247, 1985. King, R.W. and Y. Bock, Documentation of the MIT GPS analysis software: GAMIT, Mass. Inst, of TechnoL, Cambridge, 2004. Kopp, M.L. and I.G. Shcherba, History of late Alpine development of the 89 eastern Caucasus, Geotechnika, 6, 94-108,1985. Marone, C.J., Scholz, C.H., Bilham, R., On the mechanics of earthquake afterslip,/. Geophys. Res., 96, 8441-8452, 1991. Matthews, M. and R Segall, Statistical inversion of crustal deformation data and estimation of the depth distribution of slip in the 1906 earthquake, J. Geophys. Res., 98, 12153-12163, 1993. McKenzie, D.R, Plate tectonics of the Mediterranean Region, Nature, 226, 239-243, 1970. McKenzie, D.P., Active tectonics of the Mediterranean region, Geophys. J. R. Astron. Soc., 30, 109-185, 1972. McClusky, S., Personal Communication, 2006. Okada, Y., Surface deformation due to shear and tensile faults in a half-space, Bull. Seis. Soc. Amer., 75, 1135-1154, 1985. Owen, S., G. Anderson, D.C. Agnew, H. Johnson, K. Hurst, R. Reilinger, Z.-K. Shen, J. Svarc, and T. Baker, Early postseismic deformation from the 16 October 1999 Mw 7.1 Hector Mine, California earthquake as measured by survey-mode GPS, Bull. Seis. Soc. Amer., 92, 1423-1432, 2002. Peltzer G., P. Rosen, and F. Rogez, Poroelastic rebound along the Landers 1992 earthquake surface rupture, / Geophys. Res., 103,30131-30145,1998. Peltzer G., P. Rosen, F. Rogez, and K. Hudnut, Postseismic rebound in fault step-overs caused by pore fluid flow, Science, 273, 1202-1204, 1996. Philip, H., Cisternas, A., Gvishani, A., and Gorshkov, A., The Caucasus: an actual example of the initial stages of continental collision, Tectonophysics, 90 161, 1-21, 1989. Pollitz, F.F., Gravitational viscoelastic postseismic relaxation on a layered spherical Earth, J. Geophys. Res., 102, 17921-17941, 1997. Pollitz, E E , Postseismic relaxation theory on the spherical Earth, Bull. Seismol. Soc. Am., 82, 422-453, 1992. Pollitz, E E , G. Peltzer, R. Burgmann, Mobility of continental mantle: Evidence from postseismic geodetic observations following the 1992 Landers earthquake, J. Geophys. Res., 105, 8035-8054, 2000. Pysklywec, R.N., C. Beaumont, and P. Fullsack, Lithospheric deformation during the early stages of continental collision: Numerical experiments and comparison with South Island, New Zealand, J. Geophys. Res., 107, 10.1029/2001JB000252, 2002. Reilinger, R., S. Erigintav, R. Burgmann, S. McClusky, O. Lenk, A. Barka, O. Gurkan, E. Hearn, K.L. Feigl, R. Cakmak, B. Aktug, H. Ozener, and M.N. Toksoz, Coseismic and postseismic fault slip for the 17 August 1999, M = 7.5, Izmit, Turkey earthquake, Science, 289, 1519-1524, 2000. Reilinger, R., McClusky, S.C., Oral, M.B., King, R.W., Toksoez, M.N., Barka, A.A., Kinik, I., Lenk, O., and Sanli, I., Global Positioning System measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone, J. Geophys. Res., 102,9983-9999,1997. Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., Kadirov, E , Guliev, I., Stepanyan, R., Nadariya, M., Hahubia, G., Mahmoud, S., Sakr, K., ArRajehi, A., Paradissis, D., Al-Aydrus, A., 91 Prilepin, M., Guseva, T., Evren, E., Dmitrotsa, A., Filikov, S.V., Gomez, E , Al-Ghazzi, R., and Karam, G., GPS constraints on continental deformation in the Africa-Arabia, Eurasia continental collision zone and implications for the dynamics of plate interactions, in press, 2006. Ruina, A.L., Slip instability and state variable friction laws, J. Geophys. Res., 88, 10359-10370, 1983. Ruppel, C. and McNutt, M., Regional compensation of the Greater Caucasus mountains based on an analysis of Bouguer gravity data, Earth Planet. Sci. Letters, 98, 360-379, 1990. Sarker, G. and Abers, G.A., Deep structures along the boundary of a collisional belt: attenuation tomography of P and S waves in the Greater Caucasus, Geophys. J. Int., 133, 326-340,1998. Savage, J.C. and Burford, R.O., Geodetic determination of relative plate motion in central California, J. Geophys. Res., 78, 832-845, 1973. Savage, J.C and W.H. Prescott, Asthenosphere readjustment and the earthquake cycle, J. Geophys. Res., 83, 3369-3376, 1978. Savage, J.C. and J.L. Svarc, Postseismic deformation associated with the 1992 Mw = 7.3 Landers earthquake, southern California, J. Geophys. Res., 102, 7565-7577, 1997. Segall, P., Burgmann, R., and Matthews, M., time-dependent triggered afterslip following the 1989 Loma Prieta earthquake, J. Geophys. Res., 105, 5615-5634, 2000. Sheu, S.Y. and Shieh, C.F., Viscoelastic-afterslip concurrence: a possible 92 mechanism in the early post-seismic deformation of the Mw 7.6,1999 Chi-Chi (Taiwan) earthquake, Geophys. J. Int., 159, 1112-1124, 2004. Smith, D.E., Kolenkiewics, R., Robbins, J.W., Dunn, RJ., and Torrence, M.H., Horizontal crustal motion in the central and eastern Mediterranean inferred from Satellite Laser Ranging measurement, Geophys. Res. Lett., 21, 1979-1982.1994. Stark, RB. and R.L. Parker, Bounded-variable least-squares: An algorithm and application, Comput. Stat., 10, 129-141, 1995. Steketee, J.A., On Volterra's dislocation in a semi-infinite elastic medium, Can. J. Phys., 36, 192-205,1958. Tan, O., personal communication, 2005. Triep, E.G., G.A. Abers, A.L. Lerner-Lam, V. Mishatkin, N. Zakharchenko, and O. Starovoit, Active thrust front of the Greater Caucasus: The April 29, 1991 Racha earthquake sequence and its tectonic implications, J. Geophys. Res., iW, 4011-4033, 1995. Tse, S.T. and J.R. Rice, Crustal earthquake instability in relation to the depth variation of frictional slip properties, J. Geophys. Res., 91,9452-9472,1986. Turcotte, D.L. and G. Schubert, Geodynamics, Cambridge University Press, 2002. Yu, S.B., Hsu, Y.J., Kuo, L.C., and Chen, H.Y., GPS measurement of postseismic deformation following the 1999 Chi-Chi, Taiwan, earthquake, J. Geophys. Res., 108, 2520-2535, 2003. Yu, T.-T., J.B. Rundle, and J. Fernandez, Surface deformation due to a strike-93 slip fault in an elastic gravitational layer overlying a viscoelastic gravitational half-space, J. Geophys. Res., 101, 3199-3214, 1996. Zonenshain, L.P., M.I. Kuzmin, and L.M. Natapov, Alpine-Himalayan Fold-belt within the USSR, in Geology of the USSR, A Plate-Tectonic Synthesis, Geodyn. Sen, vol. 21, edited by B.M. Page, pp. 167-179, AGU, Washington, D.C., 1990. 94 Appendix A GPS Position Measurements GPS measurements in ITRF2000 reference frame: top number is the measure-ment (m) and the second number the 1-cr error. ' N D ' means that no data were no available for these stations/times. GPS POSITION DATA IN ITRF2000 SITE DIR 1991.563 1994.766 1996.739 1998.711 2000.758 N 4.84527 4.89303 4.92217 4.93942 4.96993 ± 0.00575 0.00206 0.00191 0.00169 0.00298 KHOT E 1.84394 1.94697 2.00204 2.06145 2.10662 ± 0.02016 0.00199 0.00182 0.00198 0.00520 U 1.86131 1.92139 1.90792 1.91190 1.93312 ± 0.03317 0.01131 0.00933 0.00901 0.01610 N 2.36334 2.37622 2.39601 2.41212 2.43549 ± 0.00451 0.00193 0.00109 0.00145 0.00283 KHUR E 7.61391 7.70199 7.75830 7.80795 7.85687 . ± 0.00626 0.00190 0.00112 0.00171 0.00522 U 2.39668 2.46641 2.45516 2.45531 2.47293 ± 0.02813 0.01087 0.00532 0.00817 0.01499 95 SITE DIR 1991.563 1994.766 1996.739 1998.711 2000.758 N 7.85791 7.87400 7.89575 7.91382 ND ± 0.00390 0.00164 0.00116 0.00131 ND LESO E 0.81795 0.89633 0.95746 1.00992 ND ± 0.00621 0.00160 0.00149 0.00140 ND U 0.65077 0.74676 0.73975 0.74823 ND ± 0.02743 0.00805 0.00550 0.00681 ND N 7.52641 7.57876 7.60789 ND ND ± 0.00308 0.00084 0.00121 ND ND NICH E 1.27515 1.35474 1.41534 ND ND ± 0.00347 0.00109 0.00161 ND ND U 9.79525 9.88275 9.87612 ND ND ± 0.02637 0.00376 0.00614 ND ND N 1.61635 1.66415 1.69377 ND 1.75325 ± 0.00277 0.00254 0.00130 ND 0.00284 SACC E 8.96247 9.04982 9.11276 ND 9.21484 ± 0.00519 0.00231 0.00139 ND 0.00652 U 7.01767 7.05985 7.05494 ND 7.08324 ± 0.02367 0.01290 0.00614 ND 0.01476 N 2.03446 2.06599 2.09086 2.11243 2.14293 ± 0.00449 0.00305 0.00216 0.00157 0.00262 SFRE E 4.76390 4.84061 . 4.89107 4.94660 4.99285 ± 0.00741 0.00275 0.00203 0.00202 0.00582 U 9.18377 9.27296 9.27767 9.27820 9.30947 ± 0.02764 0.01623 0.01116 0.00846 0.01395 N 9.42078 9.46442 9.49105 9.51535 ND ± 0.00267 0.00156 0.00122 0.00140 ND VANI E 2.45266 2.54597 2.60689 2.66022 ND ± 0.00322 0.00164 0.00120 0.00187 ND U 2.00303 1.99439 2.02384 1.99901 ND ± 0.02576 0.00925 0.00651 0.00769 ND 96 SITE DIR 1991.563 1994.766 1996.739 1998.711 2000.758 N ND ND ND 1.51676 1.53706 ± ND ND ND 0.00040 0.00041 ZECK E ND ND ND 1.55894 1.60892 ± ND ND ND 0.00061 0.00065 U ND ND ND 6.27418 6.27466 ± ND ND ND 0.00121 0.00138 N 1.85230 1.88473 1.90398 ND ND ± 0.00157 0.00072 0.00068 ND ND ZELB E 0.21499 0.30104 0.35676 ND ND ± 0.00227 0.00095 0.00079 ND ND U 5.21060 5.28061 5.26822 ND ND ± 0.02157 0.00343 0.00346 ND ND 97 Appendix B Conversion Velocities for ITRF2000 to Eurasian Reference Frame Velocities (mm/yr) to convert position measurements of each station from the J.TRF2000 to a Eurasian reference frame. CONVERSION V E L O C I T I E S F R O M ITRF2000 T O EURASIA STATION N O R T H (mm/yr) E A S T (mm/yr) KHOT 9.18 25.67 KHUR 9.12 25.68 LESO 9.09 25.71 NICH 8.86 25.88 SACC 9.12 25.71 SFRE 9.12 25.70 VANI 9.33 25.64 ZECK 9.53 25.35 ZELB 9.53 25.35 98 Appendix C GPS Displacements GPS displacements (in mm) used in postseismic deformation modeling, (a) 1991-1994 without secular correction (b) 1991-1994 with secular correction from time-series, (c) 1991-1994 with secular correction from block model, (d) 1994-1996 without secular correction (e) 1994-1996 with secular correction from time-series, (f) 1994-1996 with secular correction from block model. 1991-1994 NON-CORRECTED (NC) GPS DISPLACEMENTS NORTH Cr-N EAST CT-E UP CT-U KHOT 20.820 6.108 18.360 20.258 60.080 35.045 KHUR 5.830 4.906 -16.330 6.542 69.730 30.157 LESO -3.960 4.231 -13.020 6.413 95.990 28.587 NICH -3.300 3.193 23.970 3.637 87.500 26.637 SACC 5.010 3.758 18.590 5.681 42.180 26.957 SFRE -5.600 5.428 2.320 7.904 89.190 32.053 VANI 11.190 3.092 13.760 3.614 -8.640 27.370 ZELB 4.860 1.727 1.910 2.461 70.010 21.841 (a) 99 1991-1994 TIME-SERIES CORRECTED (TSC) GPS DISPLACEMENTS NORTH CT-N EAST CT-E UP cr-U KHOT 9.640 6.108 19.773 20.258 39.915 35.045 KHUR -18.613 4.906 9.545 6.542 55.484 30.157 LESO -13.259 4.231 -6.825 6.413 82.217 28.587 SACC 0.402 3.758 6.005 5.681 19.628 26.957 SFRE -10.003 5.428 -4.340 7.904 63.699 32.052 VANI 4.172 3.092 6.692 3.614 31.689 27.370 ZELB 0.544 1.727 3.292 2.461 50.633 21.841 (b) 1991-1994 BLOCK MODEL CORRECTED (BMC) GPS DISPLACEMENTS NORTH CT-N EAST cr -E UP cr-U KHOT 8.420 6.108 20.510 20.258 60.080 35.045 KHUR -22.500 4.906 5.470 6.542 69.730 30.157 LESO -20.490 4.231 -4.150 6.413 95.990 28.587 NICH 1.620 3.192 -5.380 3.637 87.500 26.637 SACC 7.050 3.758 3.730 5.681 42.180 26.957 SFRE -5.760 5.428 -5.750 7.904 89.190 32.053 VANI 1.350 3.092 9.330 3.614 -8.640 27.370 ZELB -0.780 1.727 3.920 2.461 70.010 21.841 (c) 1994-1996 NON-CORRECTED (NC) GPS DISPLACEMENTS NORTH Cr-N EAST cr -E UP Cr-U KHOT 4.41 2.809 1.102 2.697 -13.470 14.661 KHUR 5.640 2.217 1.790 2.205 -11.250 12.102 LESO 1.039 2.009 3.810 2.186 -7.010 9.749 NICH 9.530 1.473 1.165 1.944 -6.630 7.120 SACC 1.220 2.853 1.162 2.696 -4.910 14.287 SFRE -2.500 3.737 6.870 3.418 4.710 19.697 VANI 1.033 1.980 8.220 2.032 29.450 11.311 ZELB 5.700 0.990 4.400 1.236 -12.390 4.872 (d) 100 1994-1996 TIME-SERIES CORRECTED (TSC) GPS DISPLACEMENTS NORTH cr-N EAST o--E UP cr-U KHOT 5.653 2.809 3.773 2.697 -25.894 14.661 KHUR 0.386 2.217 7.922 2.206 -20.027 12.102 LESO 3.667 2.009 8.632 2.186 -15.496 9.749 SACC 0.416 2.853 12.820 2.696 -18.805 14.287 SFRE -0.720 3.737 0.522 3.418 -10.996 19.697 VANI 2.312 1.980 7.551 2.032 54.298 11.311 ZELB -0.396 0.990 4.730 1.236 -24.327 4.872 (e) 1994-1996 BLOCK MODEL CORRECTED (BMC) GPS DISPLACEMENTS NORTH cr-N EAST o--E UP cr-V KHOT 4.910 2.809 4.220 2.697 -13.470 14.661 KHUR -2.010 2.217 5.410 2.205 -11.250 12.102 LESO -0.800 2.009 0.103 2.186 -7.010 9.749 NICH -2.130 1.473 8.250 1.944 -6.630 7.120 SACC 4.510 2.853 11.420 2.696 -4.910 14.287 SFRE 1.900 3.737 -0.350 3.418 4.710 19.697 VANI 0.580 1.980 9.170 2.032 29.450 11.311 ZELB -1.210 0.990 5.120 1.236 -12.390 4.872 (f) 101 Appendix D Earth Model for Viscoelastic Modeling Earth model used in viscoelastic forward modeling with the ViscolD program. This 1-D model ranges from the Earth's surface to a depth of 1086.6 km. Earth Model for Viscoelastic Forward Modeling D E P T H 1 D E P T H 2 DENSITY B U L K M O D . S H E A R M O D . (km) (km) (g/cm3) (xlO1 0 Pa) (xl0 1 0Pa) 0 3 2.60 3.7 1.77 3 4 2.60 4.17 2.65 4 7 2.60 4.17 2.65 7 8 2.60 5.03 3.2 8 10 2.60 5.03 3.2 10 12 2.60 5.03 3.2 12 15 2.60 5.03 3.2 15 16 2.90 5.61 3.57 16 20 2.90 5.61 3.57 20 25 2.90 5.61 3.57 25 30 2.90 5.61 3.57 102 DEPTH 1 DEPTH 2 DENSITY BULK MOD. SHEAR MOD. (km) (km) (g/cm3) (X10 1 0 Pa) (xlO1 0 Pa) 30 35 2.90 5.61 3.57 35 45 2.90 5.61 3.57 45 80 3.38 13.07 6.77 80 113.5 3.37 15.00 7.00 113.5 139.1 3.37 15.00 7.00 139.1 164.7 3.38 15.00 7.00 164.7 190.4 3.39 15.00 7.00 190.4 216.0 3.39 15.00 7.00 216.0 241.6 3.40 15.00 7.00 292.9 267.2 3.42 15.00 7.00 267.2 292.9 3.44 15.00 7.00 292.9 318.5 3.47 15.00 7.00 318.5 344.1 3.51 15.00 7.00 344.1 369.8 3.55 15.00 7.00 369.8 395.4 3.60 15.00 7.00 395.4 421.0 3.66 15.00 7,00 421.0 452.3 3.71 15.00 7.00 452.3 483.5 3.76 15.00 7.00 483.5 514.7 3.81 15.00 7.00 514.7 546.0 3.85 15.00 7.00 546.0 577.2 3.90 15.00 7.00 577.2 608.5 3.96 15.00 7.00 608.5 639.7 4.03 15.00 7.00 639.7 671.0 4.11 15.00 7.00 671.0 705.6 4.21 15.00 7.00 705.6 740.3 4.32 15.00 7.00 740.3 774.9 4.40 15.00 7.00 774.9 809.5 4.47 15.00 7.00 809.5 844.2 4.51 15.00 7.00 844.2 878.8 4.54 15.00 7.00 878.8 913.4 4.54 15.00 7.00 913.4 948.1 4.54 15.00 7.00 948.1 982.7 4.54 15.00 7.00 982.7 1017.3 4.55 15.00 7.00 1017.3 1051.9 4.57 15.00 7.00 103 DEPTH 1 DEPTH 2 DENSITY BULK MOD. SHEAR MOD. (km) (km) (g/cm3) (xlQ1 0Pa) (xlQ1 0Pa) 1051.9 1086.6 4.59 15.00 7.00 104 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052519/manifest

Comment

Related Items