Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
The impact of faulting and folding on earthquake cycles in collisional orogens
(USC Thesis Other)
The impact of faulting and folding on earthquake cycles in collisional orogens
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
THE IMPACT OF FAULTING AND FOLDING ON EARTHQUAKE CYCLES IN COLLISIONAL OROGENS by Sharadha Sathiakumar A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (GEOLOGICAL SCIENCES) May 2023 Copyright 2023 Sharadha Sathiakumar Dedication To the love of my life, Hari. ii Acknowledgements “Justlikeariverthatabsorbsnutrientsfromthesoilinwhichitflows,humanbeingsassimilate wisdom and intelligence from their savant companions.” This saying, loosely translated from the Tamizh-language literary workThirukurral, accurately articulates the spirit of collaboration and friendship I have had the pleasure of enjoying during this PhD journey. The strong force of advisers, collaborators, colleagues, friends, and family have played a vital role in my pursuit of knowledge. I’d like to start by saying that no amount of “thank you” would convey the depths of gratitude I have for my advisor, Dr Sylvain Barbot. As I hail from a computer engineering background with no prior edu- cation or experience in Earth Science, I owe him for inspiring and mentoring me all through my PhD. He recognized my strengths and passion and helped me channel them to a fruitful outcome: this dissertation. I am, and will always be, grateful for his words of wisdom, encouragement, technical guidance, and un- derstanding, especially during the pandemic, which proved to be very tough for graduate students. He is a friend and a confidante, guiding and supporting me every step of the way. Thank you, for all the lessons about this Blue Planet, and French wine and cheese! I would like to thank Dr Judith Hubbard, who is a great source of inspiration. I learned to appreciate the beauty of (structural) geology from her, and I am thankful to her for instilling in me the important quality of scientific rigor. She taught me to investigate and interrogate aspects of science and life through the lens of critical reasoning. iii The following list of people, my friends, and family, are my go-to people and I thank them all for always being there for me and for all the amazing conversations about movies, music, science, and every other trivial and non-trivial aspect of life that has added value beyond measure. Thank you, Manish and Ishani, for being my pillars of strength through good and bad times. Thank you, Basu, Maggie, and Siddhu, for always being there by my side, I cannot imagine going through a tough moment without your listening ears. Thank you, Vinay, Ravi, Mrinal, Nandhu, Archana, Raja, Madhura, Abhishek, Kimaya, Laxmi, Naren, Rumi, Sufi, Maya, Mirai, and Kavi for your love, friendship, incredible laughs and mind-blowing conver- sations. Thank you, Mathilde, for being my person and, Jaclyn, for your thoughtful and caring gestures that are always accompanied by bright smiles. Thank you, Naomi, Priya, Rishav, Sagar, Bruce, Judith, and Chris for your cheerful company, friendship, and love. Thank you, Binhao, Thom, and Shiying for being understanding and supportive office mates. I am eternally thankful to my family for their everlasting and unconditional love. My parents, Sathi- akumar and Mythili, truly believe in my “superhero” capabilities (I fondly rememberamma always saying “sahasampurivaay”, which means “You will achieve great feats”) and this faith will always keep the flame of desire to dream and succeed, and do great things with my life. I am indebted to my sister and brother-in- law, Madhu and Vijay, for their unwavering support, love and affection. Thank you, my beautiful nieces, Mithra and Neha, whose smiles and love melt my heart. Thank you to the most understanding and warm people, my in-laws’ family, Shantha, Vishnu, Karthi, Pramod, and Adhithyaveer. I sincerely thank the University of Southern California (USC), and especially the Department of Earth Sciences at USC, for accepting me into their family. Insightful discussions with the faculty and colleagues of USC during coursework and department seminars and at the coffee pantry are invaluable moments of my grad life. I greatly appreciate the kind support extended by the administrative staff at USC: Cynthia Waite, Karen Young, John Yu, and Darlene Garza. In the same vein, I’d also like to thank the people at the Earth Observatory of Singapore, where I started my Earth Science career. iv Just as we all reserve the decadent dessert for the end of the meal, I end this acknowledgement with a thank you to the one person I love and adore the most in this world, Hari. Truly, words fail me, but you know that you are my role model and I always look up to you. We have faced a lot during this PhD journey, living apart, and making a lot of sacrifices, and I wouldn’t be here without you and your support. Love, forever. v TableofContents Dedication .................................................................................. ii Acknowledgements ......................................................................... iii ListofTables ............................................................................... viii ListofFigures .............................................................................. x Abstract ....................................................................................xviii Chapter1: Introduction ................................................................... 1 1.1 Motivation for this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Fault-bend folds in the shallow crustal depths of collisional zones . . . . . . . . . . . . . . 6 1.3 A new geophysical modeling framework for earthquake sequence simulations with permanent deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Layout of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5 Dissertation Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Chapter2: Earthquakecyclesinfault-bendfolds .......................................... 14 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Kinematics of fault-bend folds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.1 Fault-bend folds and axial surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.2 Short-term deformation in fault-bend folds . . . . . . . . . . . . . . . . . . . . . . 24 2.3.3 Unified representation of plasticity for faults and folds . . . . . . . . . . . . . . . . 29 2.4 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.1 Constitutive behavior for fault slip . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.2 Constitutive behavior for axial surfaces . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.3 Controls on fault dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5 Dynamics of fault-bend folds with active axial surfaces . . . . . . . . . . . . . . . . . . . . 37 2.5.1 Role of friction on the seismic cycle of non-planar faults . . . . . . . . . . . . . . . 38 2.5.2 Role of the fault bend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.5.3 Role of the sediment cut-off angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.6.1 The effect of friction and fault bend magnitude on the seismic cycle . . . . . . . . . 48 2.6.2 The effect of variable loading rate on the seismic cycle . . . . . . . . . . . . . . . . 52 2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 vi 2.8 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 3: The stop-start control of seismicity by fault bends along the Main Himalayan Thrust .......................................................................... 61 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Seismic cycle modeling of the Gorkha earthquake . . . . . . . . . . . . . . . . . . . 65 3.3.2 Impact of geometry on seismic super-cycle . . . . . . . . . . . . . . . . . . . . . . 70 3.3.3 Aftershocks at fault bends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Chapter 4: Reconciling geologic and geodetic shortening rates with active folding in the Himalayanorogenicwedge ..................................................... 80 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.1 Architecture and structural development of the Himalayan orogenic wedge . . . . 83 4.3.2 Earthquake cycles with spatially variable slip rate and active folding . . . . . . . . 88 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Chapter5: ConcludingRemarks ........................................................... 97 Bibliography ............................................................................... 102 APPENDICES ............................................................................... 124 AppendixA................................................................................. 124 Supplementary Information for chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 A.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 AppendixB ................................................................................. 128 Supplementary Information for chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 B.1 Methods: Exploring the frictional regime with rate- and state-dependent friction . . . . . . 129 B.1.1 Constitutive framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 B.1.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 B.1.3 Controls on fault dynamics from dimensional analysis . . . . . . . . . . . . . . . . 130 AppendixC................................................................................. 138 Supplementary Information for chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 vii ListofTables 2.1 Frictional and geometric parameters characterizing a synclinal fault-bend fold. The static friction coefficient, effective normal stress, reference velocity, a and b frictional parameters, Poisson’s ratio, and the down-dip width of the velocity weakening region are the same for all experiments. The magnitude of fault bend, the sediment cut-off angle and the characteristic weakening distance are varied in some experiments. . . . . . . . . . . . 39 3.1 Summaryoffrictionalparametersforthepreferredmodels: Frictional parameters and its values using model E [Elliott et al., 2016] and model H [Hubbard et al., 2016] sub-surface fault geometry in the simulations presented in Figures 2-5. The static friction coefficient, reference velocity, shear wave speed, Poisson’s ratio, effective normal stress, frictional parametersa andb, and characteristic weakening distanceL are the same for both geometry. The down-dip width of the velocity weakening region is minimally varied, such that the depth of the transition zone is∼15 km for both geometry models. . . . . . . 69 B.1 Summary of frictional parameters that explore the non-dimensional parameter R u at fixed R b = 0.2. Frictional parameters and its values using model E and model H sub-surface fault geometry in the simulations presented in Supplementary Figure B.2. The static friction coefficient, reference velocity, shear wave speed, and Poisson’s ratio are the same for all simulations. We vary the characteristic weakening distance L for simulations that explore a range of non-dimensional parametersR u . The down-dip width of the velocity-weakening region for geometry E and H is fixed for all simulations such that the depth of the transition zone is at∼ 15 km. . . . . . . . . . . . . . . . . . . . . . . 131 B.2 Summary of frictional parameters that explore the non-dimensional parameter R b at fixed R u . Frictional parameters and its values using model E and model H sub-surface fault geometry in the simulations presented in Supplementary Figure B.3. The static friction coefficient, reference velocity, shear wave speed, and Poisson’s ratio are the same for all simulations. We vary the effective normal stress, frictional parameters a andb, and characteristic weakening distanceL for simulations that explore a range of non-dimensional parametersR u andR b . The down-dip width of the velocity weakening region for models E and H is fixed for all simulations, such that the depth of the transition zone is at∼ 15 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 viii C.1 Summary of frictional parameters for Central Nepal. Frictional parameters and its values used in simulations CN1-CN5. The axial surface patches are strongly rate- strengthening simulating a visco-elasto-plastic behavior. Fault patches are velocity weakening until a depth of∼15 km, and are velocity-strengthening below that depth. . . . 139 C.2 SummaryoffrictionalparametersforEastNepal. Frictional parameters and its values used in simulations EN1-EN3. The axial surface patches are strongly rate-strengthening simulating a visco-elasto-plastic behavior. Fault patches are velocity weakening until a depth of∼15 km, and are velocity-strengthening below that depth. . . . . . . . . . . . . . 139 ix ListofFigures 1.1 Folds are ubiquitous on Earth and represent the geologic expressions of plate tectonic activity. Pictures of fold scarps in California (A, C, E), South Africa (B), and Malaysia (D and F) showing curved or bent surfaces that were previously planar to accommodate tectonic forces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Broad classification of active fault-related folding. A) Undeformed sedimentary layers of the upper crust, featuring alternating weak and competent layers. B) As slip progresses along a weak detachment fault, folding of the layers above it produces detachment folds. C) When the fault propagates upwards, slip magnitude along the fault varies, eventually dying out at the fault tip. This forces the material above the propagating tip to buckle and fold, creating fault-propagation folds. D) Slip along a flat-ramp-flat geometry produces a fold with an antiformal geometry. This type of fault-related folding is termed fault-bend folds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Mechanismofflexuralslip. Photograph from Los Angeles, California, at the Highway 14 road cut near the San Andreas Fault (Location 34.56195° N, 118.13339° W.) (A) and Nepal (C), showing flexural slip. These photographs illustrate the flexural slip mechanism with bedding-parallel slip (offsets are shown with white arrows). Photograph C is published in [Sathiakumar et al., 2020]. Photo credits for (A): Sylvain Barbot. B) Illustration of bedding-parallel slip during the development of flexural slip folds. Slip occurs at frictional bedding interfaces or contacts between layers, with no slip at the points of the smallest radius of curvature (also known as the hinge or fold axis) . . . . . . . . . . . . . . . . . . 7 1.4 Basicconceptoffault-bendfolding. When hanging wall rocks encounter fault bends, rigid block translation along fault segments creates physically impossible situations. However, folding of the hanging wall rocks accommodates fault slip, creating fault-bend folds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 x 2.1 (a) Fold scarp at the Main Boundary Thrust near Islamabad, Pakistan. The folded sediment forms a cut-off angle with the ramp (white arrows). Location: 73 ◦ 8’7.44"E, 33 ◦ 45’54.3”N. (b) Nearby fold scarp. The axial surface (dashed green line) bisects the fold limbs. Location: 73 ◦ 8’4.22"E, 33 ◦ 45’39.432”N. Photos taken by S. Barbot on 8/25/2019. (c) A seismic profile from the Hikurangi subduction zone, taken from line PEG09-25, showing a thrust ramp rising from a décollement, with the development of a synclinal fault-bend fold in the hanging wall. The axial surface bisects the fold. Red line shows décollement, white arrows show thrust ramp (fault plane reflectors), blue line shows schematic folded stratigraphy. Scale is approximately 1:1. The line was collected as part of the New Zealand Government’s PEG09 2D reconnaissance seismic survey in the Pegasus Basin. (d) Photograph from Nepal, showing flexural slip. The photo shows sub-vertical bedding (seen as color banding) with a clastic dike (thin blocky gray line) cutting across the bedding. The dike is offset by 10s of cm of slip along the bedding planes. White arrows show slip, black arrows show terminations of the dike at the slip plane. This slip occurred to accommodate the folding of this stratigraphic section from its original orientation (horizontal) to its current orientation (vertical). The scale is shown with a geological hammer. Photo taken by J. Hubbard on 04/29/2017. Location: 24.74615 ◦ N, 82.83798 ◦ E. . . . 18 2.2 Synclinal (a) and anticlinal (b) fault-bend folds used in our model. The axial surface (black dotted line) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angleγ andϕ represents the change in dip between the two fault segments. Here, the uniform color in the hanging wall represents uniform local slip rate (v 1 =v 2 ) on the fault segments represented by the red solid line. Synclinal (c) and anticlinal (d) fault-bend folds with angular fault cut-off angle θ . The axial surface (black dotted line) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angleγ . Here, the color of the hanging wall changes across the axial surface, showing the increase (v 1 <v 2 in the case of c) and decrease (v 1 >v 2 in the case of d) of local slip rate on the ramp compared to the décollement as a result of folding through the fault bend. (e, f) The three velocities are plotted together, forming a closed loop. . . . . . . . . . . . . . . . . . . 22 2.3 Velocity vector diagrams of a syncline with zero cut-off angle to depict the relative displacements (velocities) across the fault and axial surface, assuming it is a fixed block. v A represents the velocity of block A with respect to block B,v B represents the velocity of block B with respect to block C, andv C represents the velocity of block C with respect to block A. If we move in clockwise direction instead,v 1 represents the velocity of block 1 with respect to block 2,v 3 represents the velocity of block 2 with respect to block 3, and v 2 represents the velocity of block 3 with respect to block 1. . . . . . . . . . . . . . . . . 26 2.4 Synclinal (a) and anticlinal (b) fault-bend folds in three-dimensions. The axial surface (green surface) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angle γ , and the fault segments are represented by red solid line. (c) The 3-D perspective view of an axial surface in the reference coordinate system (s,d,n) that forms an orthonormal right-handed system of coordinates, wheres,d, andn are the unit vectors in the strike, dip, and normal directions of the axial surface, andb is the Burger’s vector showing slip direction. (d) Planar view of the Burger’s vectorb that describes the direction of motion on the axial surface. The take-off angle β is defined between the dip vector and normal vector. (e) The traction vector (red arrow) on the axial surface and its decomposition into normal (orange arrow) and shear (blue arrow) components. . . . . . . 28 xi 2.5 Simulation results for a synclinal fault-bend fold, with a ramp dip of 40 ◦ and lower segment dip of 30 ◦ (change in dipϕ = 10), for varyingR u values. An illustration of the geometry and fault friction properties is shown in the top left panel. Velocity-weakening regions are shown in red and velocity-strengthening regions in black. The dotted black line indicates the axial surface, which is tied to the bend in the fault. (a)-(o) Seismic cycles: The color indicates the fault slip rate; red reflects seismic velocities and dark blue fault locking. Each panel shows about 2000 years of simulation for this geometry as the slip weakening distance L is varied. The x-axis depicts time steps that are non-uniformly spaced, not time itself, with fast events having a finer resolution than during interseismic periods. As R u increases, we find bimodal seismicity, with partial and through-going ruptures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6 Simulation results for the four geometries used in our study (ϕ = 10 ◦ ,ϕ = 20 ◦ ,ϕ = 30 ◦ , ϕ = 40 ◦ ) for various R u values. Geometry and fault friction properties used in the model are shown in the first row. Velocity-weakening regions are shown in red and velocity-strengthening regions in black. The dotted black line indicates the axial surface, which is tied to the fault bend. (a)-(l) Seismic cycles: The color indicates the fault slip rate; red reflects seismic velocities and dark blue, fault locking. Each panel represents about 2,000 years of simulation for a given geometry. The slip weakening distanceL is varied, showing changes in earthquake slip behavior as the magnitude of the fault bend increases. 43 2.7 Simulation results for the model (ϕ = 40 ◦ , and θ = 10 ◦ , θ = 20 ◦ , θ = 30 ◦ , θ = 40 ◦ ) for various R u values. The geometry and fault friction properties are shown in the first row. VW regions are shown in red and VS regions in black. The dotted black line indicates the axial surface, also denoting the position of the kink. Panels A to O) Seismic cycles: The color indicates the fault slip rate: red reflects seismic velocities and dark blue fault, locking. Each panel represents 2,000 years of simulation for a given geometry. The slip weakening distanceL is varied, showing changes in earthquake slip behavior as the cut-off angle increases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.8 The effect of R u and change in fault dip (ϕ) on slip behavior of non-planar faults with one fault bend and planar faults with no fault bends (ϕ=0 ◦ ). Asϕ increases, the rupture complexity increases. Rupture segmentation depends onϕ forR u values greater than 10. For smaller values of R u , multi-pulse ruptures are more common for large fault bends than for gentle fault bends. In the case of planar faults, there are no multi-pulse or partial ruptures in the shallow parts of the fault. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.9 The effect of change in fault dip (ϕ) and cut-off angle (θ ) on the slip behavior of non-planar faults for uniform (θ = 0 ◦ ) and variable (θ ̸= 0 ◦ ) local slip rate on the two segments of the fault for a small value ofR u = 6.6 (left) and a large value ofR u = 44.4 (right). As the dip of the incoming thrust sheet increases (asθ increases), more multi-pulse ruptures are observed forR u = 6.6. Similarly, steep incoming thrust sheets allow shallow partial ruptures in contrast to gentler incoming sheets that show smaller or no ruptures in the shallow region. Therefore, as θ increases, ramp ruptures for all values of ϕ become prevalent. Large cut-off angles produce surface-breaking shallow ramp ruptures. . . . . . . 54 xii 2.10 Slip accumulation on the fault and axial surface plotted against the down-dip distance for a model withϕ = 40 ◦ ,θ = 40 ◦ , andR v = 1.6. Slip is much higher on the ramp than the décollement due to the high cut-off angle and large fault bend. The geometry of the fault and the simulation results can be found in Figure 2.7o. . . . . . . . . . . . . . . . . . . . . 56 2.11 (a). The 3-D view of a décollement-ramp structure showing the hanging wall deformation due to folding for the model withϕ = 40 ◦ ,θ = 40 ◦ ,R v = 1.6. Slip on the fault is shown by the color: aseismic slip in gray, yellow lines indicate through-going ruptures, green lines are partial ruptures on the décollement, blue lines indicate partial ruptures on the ramp that do not break the surface, and pink lines are surface-breaking partial ruptures on the ramp. (b) Simulation results computed for 2000 years for the above geometry, showing through-going surface ruptures and partial ruptures with varying magnitudes (cumulative force in 2-D). Repeating rupture patterns analogous to supercycles are highlighted with green dashed boxes. (c, d). Maximum slip velocity plotted against time on the ramp (c) and décollement (d). The geometry of the fault and the simulation results can be found in Figure 2.7o. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.1 Seismo-tectonic setting and fault geometry of Nepal’s Kathmandu Valley. a Distribution of coseismic slip [Qiu et al., 2016a] in the source region of 2015 M w 7.8 Gorkha, Nepal earthquake (background colors). GNSS measurements are represented by circles for vertical displacement and arrows for horizontal displacements. Near-field GNSS station names are highlighted in yellow and the far-field stations in gray. The black dots represent the aftershock sequence recorded by the Nepal Array Measuring Aftershock Seismicity Trailing Earthquake (NAMASTE) network [Mendoza et al., 2019]. The focal mechanisms are for the immediate aftershock sequence [Wang et al., 2017a]. b Sub-surface geological cross-section along profile A-A’ constrained using surface geology [Hubbard et al., 2016]. c Geological cross-section constrained by coseismic geodetic data [Elliott et al., 2016]. The velocity-weakening regions are indicated in red and velocity-strengthening regions in black. The focal mechanisms (colored by focal depth and scaled to magnitude) are superimposed on the cross-section. These cross-sections show the major thrust faults in the region, including the surface-emergent MFT (Main Frontal Thrust), Main Boundary thrust (MBT), and Main Central Thrust (MCT). . . . . . . 64 3.2 Comparison between geodetic data and simulated surface slip for the 2015 Gorkhaearthquake. a,b The two tracks of Advanced Land Observing Satellite (ALOS) line-of-sight displacements that recorded the coseismic displacements. c, d Simulated surface displacements using model E [Elliott et al., 2016] and model H [Hubbard et al., 2016] compared with ALOS2-157 and ALOS2-048 InSAR scenes. e, f, g Simulated horizontal (north and east components) and vertical displacements indicated using black solid (model H) and dashed (model E) profiles compared against GNSS horizontal and vertical data. The blue data points correspond to near-field stations indicated by yellow labels in Figure 3.1 and red data points correspond to far-field stations indicated by gray labels. h Simulated fault slip using our numerical method, indicated in black solid (model H) and dashed (model E) lines. Fault slip estimated by inverting geodetic data with geometry H (orange solid) and geometry E (orange dashed) lines. The error bars indicate an uncertainty of 5 percent added to the data. . . . . . . . . . . . . . . . . . . . . . . . . . 71 xiii 3.3 Time history of the megathrust behavior, showing earthquake cycles using two end-membergeometrymodels. a Spatio-temporal evolution of slip on the Himalayan megathrust using geometry model E [Elliott et al., 2016]. b Similar simulations using geometry model H [Hubbard et al., 2016]. Earthquakes of different sizes are differentiated using colors; red indicates giant surface-breaking ruptures, purple represents partial ruptures of varying sizes, blue represents Gorkha-like events, and black represents aftershocks. Stars indicate the hypocenter locations. The initial and final slip velocities for seismic events are 0.1 and 0.001 meters per second, respectively. . . . . . . . . . . . . . 72 3.4 Faultslipandshearstressillustratingfaultbendsattractingpersistentseismicity. Slip accumulated on the fault plotted against the down-dip distance using geometry E (a) and geometry H (c). The color indicates the fault slip rate, with red indicating seismic slip and blue represents interseismic creep. The black contours are cumulative coseismic slip values plotted every 5 seconds. The white contours are cumulative aseismic slip values plotted every 50 years. The large slip deficit in the shallow regions of the fault is recovered during through-going ruptures. The stars represent hypocenter locations. (b, d) Shear stress accumulation on the fault; red regions indicate areas of high stress accumulation and blue indicates areas of low stress. Previous ruptures and geometric complexities causes stress concentrations where earthquakes often nucleate. Hypocenter locations are represented by stars, yellow stars represent Gorkha-like events. . . . . . . . . . . . . . . . 74 3.5 Aftershocks at fault bends following the Gorkha earthquake. a Geological cross- section along profile A-A’ that describes the MHT geometry using model H [Hubbard et al., 2016]. b The 2015 Gorkha event and the seven early aftershocks during the 12-hour time window relocated using advanced waveform techniques [Wang et al., 2017a]. We monitor the slip velocity at various observation points on the fault surface. The middle ramp consists of three observation points: Points A and C represent the position of the kink formed at the two fault bends of the middle ramp, and B is located at its center. Points D and E are on the décollement that ruptured during the 2015 event. Point E experiences the coseismic slip of the simulated partial event H65, while D represents the termination point that is completely locked. c The evolution of peak slip velocity representing the maximum slip speed attained anywhere in the model, showing the Gorkha mainshock and the subsequent updip aftershock that happened after 14 hours. d,e The evolution of slip velocity is plotted with the colors indicating the position of the observation points. . . 76 4.1 Perspective view of the Himalayan subsurface and observational constraints on the shortening-rate in Nepal. The three-dimensional fault geometry and the cross-sectional surface across Central Nepal are based on Hubbard et al. [2016]. The fault geometry at two representative cross-sections in Central and East Nepal (CN and EN, respectively) is highlighted, with the frictionally unstable seismogenic zone in red and the creeping region in black. The contours for the last known earthquakes in these segments, theM w 8.2 1934 Nepal-Bihar earthquake (EN) and theM w 7.8 2015 Gorkha earthquake (CN), are represented in yellow. Central Nepal features a middle ramp that is absent in East Nepal. Fault system evolution allows break-forward episodes of one or more ramps, allowing accretion of Indian lithosphere onto the upper plate. The local long-term slip on individual fault segments, which may vary down-dip due to internal deformation of the orogenic wedge, is represented by double arrows. . . . . . . . . . . . . . . . . . . . . . . . 82 xiv 4.2 StructuralevolutionoftheMHTinCentralandEastNepal. a) The locations of the two-dimensional profiles in Central (CC’) and East (EE’) Nepal are shown on a geological map of Nepal. b) Illustration of fault-bend fold theory concepts applied to reconstruct the structural development of the Himalayan megathrust. The long-term slip rate of fault segments varies depending on the fold shape and the history of deformation. The Main Frontal Thrust (MFT) can slip at a rate faster (pink) or slower (yellow) than the long-term loading rate due to active folding and shortening. Reconstruction of hanging wall evolution along profile EE’ in East Nepal (c-e) and profile CC’ in Central Nepal (f-k) . 86 4.3 Seismotectonic implications of down-dip long-term slip-rate variations. a and b) Illustrations of two-dimensional fault geometry models for Central and East Nepal with fault segments and sample active axial surfaces. Unstable weakening is denoted in red, velocity strengthening in black. R u and R b are non-dimensional numbers that control the dynamics of fault behavior in the rate-and-state friction framework. R u is the ratio of seismogenic zone width to characteristic nucleation size, and R b is the ratio of static to dynamic stress drop. c and d) Down-dip rupture width illustration for megathrust-spanning full ruptures (great ruptures), multiple-segment large ruptures (major ruptures), and single-segment ruptures (strong ruptures); stars indicate hypocenter locations. e) Colored histogram indicating frequency of ruptures of a given width for each of the stages in our simulations (CN1-CN5, EN1-EN3). Darker colors indicate more frequent ruptures. The lower bar indicates the slip rate of the MFT for each stage; colors of the lower bar are taken from Figure 4.2f. Recurrence times of great ruptures for all models (black: Central Nepal stages CN1-CN5; blue: East Nepal stages EN1-EN3). . . . . . 89 4.4 Variationsoflong-termslip-ratesontheMFTovermulti-millionyeartimescales due to active folding and duplexing compared to the geodetic loading rate. a and b) Fault-bend folding predicts different frontal fault slip rates depending on the history of kinematic development (black and blue horizontal lines representing Central and East Nepal, respectively). In Central Nepal (a) the geologic shortening rate at the MFT (estimated from fluvial terraces, indicated by the shaded brown region) is 20-30% higher than the geodetic loading rate for the MHT (shaded yellow region). An increase in frontal long-term slip rate is possible due to the uplift of inclined rock layers in Central Nepal in stages CN4-CN6. In East Nepal, the geologic rates are highly uncertain, and so it is not currently possible to distinguish between different stages of evolution based on rates alone. c) Evolution of long-term slip rate on the shallow MHT and MFT in Central Nepal (black) and East Nepal (blue). The modified slip rates due to structural evolution are sustained for many seismic super-cycles, spanning hundreds of thousands to a few million years. d) Geological cross-section of Central Nepal Hubbard et al. [2016], showing inclined beds climbing the middle-ramp, causing a long-term slip rate increase in the up-dip segments, similar to Stage CN4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 A.1 Convergence test performed on a synclinal fault-bend fold, with a ramp dip of 40 ◦ and lower segment dip of 0 ◦ (change in dipϕ = 40),θ = 40 ◦ , andR v = 1.6 for varyingR u values and patch sizes. Large patch sizes do not resolve the model well and numerical convergence is achieved for smaller patch sizes. . . . . . . . . . . . . . . . . . . . . . . . . 127 xv B.1 Two-dimensional fault geometry models E and H. The fault geometry along a representative cross-section perpendicular to the plate boundary that crosses the Gorkha- Pokhara Anticlinorium in the middle of the Gorkha coseismic rupture. The dip angles of the individual fault segments are marked for model E (blue) and model H (red). The transition zone of both models is located between 80 and 100 km from the surface trace of MFT (Main Frontal Thrust). The depth of the transition zone is about 15 km in both models. Both the shallow bend and the mid-crustal ramps have similar dip angles and are present at comparable depths. The main difference in fault structure between the two models is in the middle of the seismogenic zone: model E has a single décollement that connects the shallow bend to the mid-crustal ramp in comparison to the more complex fault structure with a middle ramp in model H. . . . . . . . . . . . . . . . . . . . . . . . . . 132 B.2 Earthquake cycle dynamics exploring a large range ofR u values. a-n) Simulation results for model H (left column) and model E (right column). The dotted white lines indicate the position of the fault bends and the transition zone that represent areas of high stress concentrations. The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The slip-weakening distanceL ranges between 0.5 m to 0.01 m which results inR u values to range between 1 and 60. The x axis in the plots show time steps rather than time itself, which is non-uniformly spaced such that fast events have a finer resolution than during interseismic periods. As R u increases, slip complexity increases for both geometry models. The frictional parameters are listed in Supplementary Table B.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 B.3 EarthquakecycledynamicsexploringR b values. a-j) Simulation results for model H (left column) and model E (right column). The dotted white lines indicate the position of the fault bends and the transition zone that represent areas of high stress concentrations. The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The friction parametersa,b,L are varied to explore theR u − R b phase space. The x-axis in the plots show time steps rather than time itself, which is non- uniformly spaced such that fast events have a finer resolution than during interseismic periods. As R b increases, pulse-like ruptures and fault bend seismicity increases. The frictional parameters are listed in Supplementary Table B.2. . . . . . . . . . . . . . . . . . . 134 B.4 Pulse-like vs crack-like rupture propagation. With higher R b = 0.9 values, the rupture transitions from crack-like to pulse like with both model E and model H (a and c). In contrast, with lower values ofR b = 0.2, the rupture is fully crack-like (b and d). The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The x axis in the plots show time steps rather than time itself, which is non-uniformly spaced such that fast events have a finer resolution than during interseismic periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 xvi B.5 Statistical properties of simulated events using the two end-member geometry models. Simulation results computed for 5,000 years for a) model E and b) model H, showing giant surface ruptures (red stars), partial ruptures of varying sizes (purple stars), Gorkha-like events (blue stars), and aftershocks (black stars) with varying magnitudes (cumulative force in 2-D). c) Cumulative histogram of all events for a catalog spanning 45,000 years for both geometries. Model H approximately follows a power-law distribution indicated by the dashed black line. The Pareto distribution (dashed line) is for a power exponentα =0.3 for moment per unit length above5× 10 7 N. . . . . . . . . . . . . . . . 136 C.1 Two-dimensionalfaultgeometrymodelsalongCentralandEastNepal. The fault geometry with the dip angles of each fault segment marked in Central (CN) and East Nepal (EN). At every stage of the structural fault evolution, the fault geometry and frictional properties of Central and East Nepal remain constant, but the orientation of the axial surface tied to each fault bend changes according to the orientation of the thrusted sheets. Velocity weakening (unstable fault slip) regions are marked in red, velocity strengthening regions are marked in black. The transition zone where the behavior changes from velocity strengthening to velocity weakening is at∼15 kms depth. . . . . . . . . . . . . . . 140 C.2 EarthquakecycledynamicswithslipratesegmentationinCentralNepal. a). The fault geometry of Central Nepal featuring the middle ramp in the seismogenic zone, based on Hubbard et al. [2016] fault model. The fault bend locations and the transition between velocity-strengthening and velocity-weakening zones are highlighted. b-f). Seismic cycle simulations of the various stages in Central Nepal (CN1-CN5) to illustrate the evolution of seismic behavior as the hanging wall deforms. The dashed white lines correspond to the locations of fault bends, and solid line represents the transition zone. The x-axis represents time index, not time itself, that is non-uniformly spaced, so that seismic events have finer resolution compared to aseismic periods and the y axis represents the down-dip fault distance in km. The color indicates the fault slip velocity, red indicates fast seismic events and dark shades of blue indicates fault locking. The frictional parameters are listed in Supplementary Table C.1. Variations in slip rate impact earthquake cycle characteristics like rupture width and magnitude, hypocenter location, and the recurrence times of great and partial ruptures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 C.3 Earthquake cycle dynamics with slip rate segmentation in East Nepal. a). The fault geometry of East Nepal without the middle ramp in the seismogenic zone, based on Hubbard et al. [2016] fault model. The fault bend locations and the transition between velocity-strengthening and velocity-weakening zones are highlighted. b-d). Seismic cycle simulations of the various stages in East Nepal (EN1-EN3) to illustrate the evolution of seismic behavior as the hanging wall deforms. The dashed white lines correspond to the locations of fault bends, and solid line represents the transition zone. The x-axis represents time index, not time itself, that is non-uniformly spaced, so that seismic events have finer resolution compared to aseismic periods and the y axis represents the down-dip fault distance in km. The color indicates the fault slip velocity, red indicates fast seismic events and dark shades of blue indicates fault locking. The frictional parameters are listed in Supplementary Table C.2. Variations in slip rate impact earthquake cycle characteristics like rupture width and magnitude, hypocenter location, and the recurrence times of great and partial ruptures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xvii Abstract Our planet Earth writhes with continuous tectonic motion; mountains rise in collisional margins, oceans cease to exist in subduction zones, and new oceans are born in divergent margins. Most of the tectonic activity primarily occurs on a network of interconnected crustal faults that together accommodate the relative motion between two or more tectonic plates. Fault motion generates natural catastrophes like earthquakes, landslides, and tsunamis, as well as builds mountain chains that fabricate the Earth’s aerial landscapes. In this thesis, I study the subsurface fault structure complexity and its impact on seismotec- tonics by integrating geological and geophysical tools to better understand crustal deformation processes. Convergent tectonic margins like fold-and-thrust belts and subduction zones host some of the largest and most disastrous seismic hazards in the world. In these regions, geological and geophysical models and observations of fault systems reveal several geometric features such as fault bends, jogs, branches and splays, but the role of these subsurface structural features on earthquake behavior has often been overlooked. Fault bends and other structural complexities not only affect seismic cycle characteristics, but are also foundational to shaping the Earth’s subterranean and aerial surfaces. Therefore, as a first step, I develop physics-based numerical models of seismic cycles that incorporate fault structure complexity. To achieve this, I use fault-bend folding theory, based on observations of structural geology, to define a new geophysical modeling framework that accounts for changes in fault orientations and the accompany- ing folding processes in the crust. The new framework allows us to explore how changes in fault geometry can impact the earthquake cycle. These more realistic rupture models with a single fault-bend fold show xviii that fault geometry and fault friction properties interact to shape the seismic cycle. Earthquake behavior is sensitive to fault friction parameters, the magnitude of fault bending, and the angularity of the hanging wall sedimentary packages. To understand the role of fault bends on the short-term seismic cycle processes, I study the Main Himalayan Thrust (MHT), a highly corrugated active fault system characterized by multiple fault bends, that separates the Indian subcontinent from the Eurasian Plate. Fault bends are usually thought to be segment barriers, affecting earthquake segmentation processes and coseismic rupture extent. The 2015 M w 7.8 Gorkha, Nepal earthquake was one such event where the rupture extent was hindered by fault bends. I use physics-based numerical models of earthquake sequences to explore a wide range of frictional and structural complexities to explain the Himalayan megathrust behavior in the Kathmandu section, the source region of the 2015 event. I show that fault bends are the loci of persistent seismicity, affecting inter-, co-, and post-seismic processes. These models tuned to geological and geophysical data and observations reveal that a ramp in the middle of the seismogenic zone explains the abrupt termination of the Gorkha earthquake and the source mechanism of up-dip aftershocks. This study helps better constrain the MHT geometry and illustrates the importance of high-resolution fault geometry to build more reliable models of seismic hazard assessment at collisional margins. Next, I investigate the Himalayan megathrust further to understand how the long-term mountain- building processes affect short-term seismic cycle processes. Measurements of shortening from geologic and geodetic methods provide important constraints on fault behavior and seismic hazard. However, we often encounter discrepancies in these estimates of fault motion. For instance, in Central Nepal, geological datasets provide fault slip rates that are 20-30% higher than recent geodetic fault slip rates. In this study, we build comprehensive structural models to reconcile these contrasting observations by accounting for folding in the orogenic wedge during the structural evolution of this range-bounding megathrust. We show that the mountain-building processes due to active shortening and folding along the MHT can alter xix the local slip rate of individual fault segments. During periods spanning several hundred thousand to millions of years, the frontal regions of the MHT within the seismogenic zone can slip faster or slower than the long-term shortening rate, depending on the internal deformational phase of the fault. Therefore, geodesy provides reliable estimates of long-term convergence between the Indian and Eurasian plates, while geologic datasets that are most sensitive to the shallow deformation, provide important clues about the wedge deformation. The spatial and temporal variations in the local loading rate of fault segments can induce different patterns of earthquake segmentation, with implications for seismic hazards. I conclude my thesis by summarizing the findings from the above research topics. As a young and highly challenging scientific field of inquiry, where useful and reliable datasets for testing our hypothesis are hard to come by, we stand to gain better fundamental insights about fault behavior by integrating existing geophysical and geological models and datasets. A multidisciplinary approach to understanding crustal deformational processes would help build better hazard and risk assessment models that have great humanitarian relevance and value. xx Chapter1 Introduction 1.1 Motivationforthiswork The surface of the Earth exhibits spectacular landforms as a product of active plate tectonics [Molnar, 1988, Davies, 1992, Korenaga, 2013]. Without the continuous slow movements and the occasional joggle of these tectonic plates, mountains would eventually disappear [Suppe et al., 1981, Jordan et al., 1983, Avouac, 2003, Dal Zilio et al., 2021]. As the Earth’s crust warps itself into folded structures to build complex terrain, environmental forces like the wind and water slowly act on the surface. The combined effect of tectonics and erosion built the landscapes around us and shaped an environment that favored human evolution [King and Bailey, 2006, Trauth et al., 2010, Bailey and King, 2011]. Faults and folds are geologic expressions of plate tectonic activity [Rich, 1934, Rettger, 1935]. They ac- commodate the tectonic forces and contribute to the growth of topography on the Earth’s surface [Stearns and Matthews, 1978, Golombek et al., 1991]. Geologists and geophysicists study the structure and behavior of faults and folds because they play an important role in trapping oil and gas [Hopkins, 1917, Mitra, 1990, Ingram et al., 2004], creating minerals of economic value [Cox et al., 1991, Kisters, 2005], and generating natural hazards like earthquakes [Allmendinger and Shaw, 2000, Chen et al., 2001], landslides [Keefer, 1 1984, Basharat et al., 2021], and tsunamis [Qiu and Barbot, 2022]. An in-depth understanding of fault- ing and folding can help decipher the complex crustal deformation processes that build and maintain this ever-evolving system. Folds in nature can range in size from microscopic wrinkles to Everest-sized ones (Figure 1.1). Folds in the lower crust are formed under high pressure and temperature conditions to accommodate ductile deformation [Hobbs et al., 2007, Platt et al., 2015, 2018]. However, the upper crust is dominated by brittle deformation (fractures, joints, and faults) and folding is a consequence of intense shear of soft and weak sedimentary materials [Brandes and Tanner, 2014, Nabavi and Fossen, 2021]. Folding is also common in the upper crust due to faulting, and active fault-related folding is a dominant deformation mechanism [Yeats, 1986, Van der Pluijm and Marshak, 2005, Fossen, 2016] occurring in various tectonic regimes, from fold- and-thrust belts [Suppe, 1981, Tapponnier et al., 1990, Cattin and Avouac, 2000, Avouac, 2003, Lacassin et al., 2004] and accretionary wedges in subduction zones [Biju-Duval et al., 1982, Gulick et al., 2004, Bradley et al., 2019], to strike-slip [Wilcox et al., 1973, Sylvester, 1988, Tindall and Davis, 1999], and ex- tensional regimes [Fossen and Holst, 1995, Howard and John, 1997, Sharp et al., 2000, Harris et al., 2002, Ferrill et al., 2012]. 2 Figure 1.1: Folds are ubiquitous on Earth and represent the geologic expressions of plate tectonic activity. Pictures of fold scarps in California (A, C, E), South Africa (B), and Malaysia (D and F) showing curved or bent surfaces that were previously planar to accommodate tectonic forces. 3 Many kinds of fault-related folds exist [Jamison, 1987, Brandes and Tanner, 2014, Nabavi and Fos- sen, 2021], namely detachment folds, fault-propagation folds, and fault-bend folds (Figure 1.2). Tectonic shortening along fault segments with bends produces fault and fold scarps on the surface, and these mor- phological features help us identify subsurface fault shape and the kinematics and history of deforma- tion. Detachment folds most commonly occur above a low-angle fault (often referred to as detachment) or bedding-parallel faults (referred to as décollement) [Poblet and Stuart, 1995, Poblet et al., 1997, Rowan, 1997, Mitra, 2003, Scharer et al., 2004, Li and Mitra, 2020]. Here, a ramp does not cut across the layers, and the fault runs along the near-horizontal stratigraphic layers. The competent rock layers above the weak detachment or décollement fold or buckle in response to shortening. Usually, the weak layers are over-pressurized salt or shale, while the competent layers are made of sandstone or limestone (Figure 1.2b). When faults start to propagate upwards from the deeper to the shallower regions of the crust, the slip along the fault reduces and eventually becomes zero at the fault tip. The absence of slip at the propagating tip is compensated by the folding of layers above this tip. This type of folding above a fault tip is termed as fault-propagation folds and is commonly associated with low-angle thrust faults [Suppe and Medwedeff, 1990, Mitra, 1990, Shaw et al., 2005, Hughes et al., 2014, Hughes and Shaw, 2015]. These folds are char- acterized by an asymmetric geometry that features a gently-dipping backlimb and a steep or overturned forelimb (Figure 1.2c). 4 A B C D undeformed sedimentary layers detachment fault detachment fold fault-bend fold change in fault orientation décollement décollement ramp fault tip buckling above fault tip fault-propogation fold Figure 1.2: Broad classification of active fault-related folding. A) Undeformed sedimentary layers of the upper crust, featuring alternating weak and competent layers. B) As slip progresses along a weak detachment fault, folding of the layers above it produces detachment folds. C) When the fault propagates upwards, slip magnitude along the fault varies, eventually dying out at the fault tip. This forces the material above the propagating tip to buckle and fold, creating fault-propagation folds. D) Slip along a flat-ramp- flat geometry produces a fold with an antiformal geometry. This type of fault-related folding is termed fault-bend folds. 5 This thesis focuses on fault-bend folds commonly associated with collisional zones [Suppe, 1983]. Fault-bend folds in convergent margins are formed when rocks climb non-planar fault surfaces (Fig- ure 1.2d). They record the long-term tectonic deformation over millions of years and represent the complex interaction of fault processes in geological timescales [Suppe, 1983, Shaw et al., 2005]. These morpholog- ical complexities affect both the long-term mountain-building processes and the short-term seismic cycle processes. In this thesis, I investigate the role of fault-bend folds in earthquake behavior during the seismic cycle and its role in orogenic processes in the long term. I combine geophysical and geological observa- tions and models, and numerical modeling to study these geologic structures that have implications for seismic hazards. 1.2 Fault-bendfoldsintheshallowcrustaldepthsofcollisionalzones In active tectonic margins, fault segments often encounter abrupt changes in orientation. When tectonic blocks slip along such bends, map-scale folding occurs. This process creates permanent deformation of the crust and contributes to building topography. A fundamental mechanism in the formation of folds at shallow crustal depths is flexural slip [Donath and Parker, 1964, Tanner, 1989, Couples et al., 1998, Chester et al., 1991, Damasceno et al., 2017, Johnson, 2018]. Flexural slip is the layer-parallel slip occurring along weak bedding interfaces, as a consequence of material anisotropy, creating flexural-slip folds. This process happens because the hanging wall rocks are made of sedimentary thrust stacks with weak mechanical interfaces between well-lithified rock layers. Evidence of folds and flexural slip abounds in outcrops and subsurface seismic reflection profiles (Figure 1.3). 6 Cut-off angle ~1 ft ~1 ft B bedding-parallel slip bedding-parallel slip Bedding C A hinge no slip Figure 1.3: Mechanism of flexural slip. Photograph from Los Angeles, California, at the Highway 14 road cut near the San Andreas Fault (Location 34.56195° N, 118.13339° W.) (A) and Nepal (C), showing flexural slip. These photographs illustrate the flexural slip mechanism with bedding-parallel slip (offsets are shown with white arrows). Photograph C is published in [Sathiakumar et al., 2020]. Photo credits for (A): Sylvain Barbot. B) Illustration of bedding-parallel slip during the development of flexural slip folds. Slip occurs at frictional bedding interfaces or contacts between layers, with no slip at the points of the smallest radius of curvature (also known as the hinge or fold axis) Consider a simple case of two fault segments in cross-section in the hanging wall with different orienta- tions, connected together (Figure 1.4a). When the two blocks slide past one another, rigid block translation along the upper fault segment produces a void between the two blocks (Figure 1.4b). However, rocks in nature are not strong enough to sustain such a void in the crust. If the rigid translation occurs along the lower segment, it produces an overlap, which is also not supported by the crustal rocks (Figure 1.4c). In reality, folding of the hanging wall rock layers accommodated by layer-parallel flexural slip, through the development of a kink band, allows for permanent deformation in the crust (Figure 1.4d). This process of 7 folding due to bends in the underlying fault is described as fault-bend folding [Suppe, 1983, Shaw et al., 2005]. At the fault bend, rock layers are folded and then transported along the upper segment. The fold axis along which folding occurs is called the active axial surface. The inactive axial surface keeps track of the rock particles that were folded along the active axial surface and eventually moves away as fault slip progresses. Therefore, the width of the kink band allows us to probe into the amount of fault slip. slip void slip overlap upper segment lower segment change in dip inactive axial surface active axial surface kink band slip Figure 1.4: Basicconceptoffault-bendfolding. When hanging wall rocks encounter fault bends, rigid block translation along fault segments creates physically impossible situations. However, folding of the hanging wall rocks accommodates fault slip, creating fault-bend folds. Fault-bend fold theory describes the kinematic deformation of hanging wall rocks as they ride over non-planar fault segments using simple geometric relationships. The theory makes useful assumptions about how mass is conserved. The assumed incompressibility of rocks results in the conservation of layer thickness, conservation of bedding length, and conservation of cross-section area of the overall upper plate over long time scales [Suppe, 1983]. Following these assumptions, the fold axis or the axial surface bisects the angle between the two limbs of the fold. Axial surfaces are the loci of intense deformation, strain accumulation, and off-fault deformation processes. In continental collisional regions and subduction zones, the décollement along which shortening oc- curs, transfers displacement to a higher stratigraphic level by thrusting along high-angle ramps. This ramp-décollement structure encourages the formation of fault-bend folds in these regions and accom- modates shortening via folding. Nepal [Hubbard et al., 2016], Taiwan [Suppe and Namson, 1979, Chen et al., 2007], and the outer wedge of subduction zone [Bradley et al., 2019, Watt and Brothers, 2021, Qiu and Barbot, 2022], are some prominent examples of convergent tectonic boundaries with an abundance 8 of fault-bend folds. Fault-bend folds are more commonly used to build regional and local balanced cross- sections when there is negligible internal deformation in the thrust stacks and the fold limbs are relatively symmetric [Endignoux and Mugnier, 1990, McDougall and Hussain, 1991, Mukhopadhyay and Mishra, 2004, Almeida et al., 2018]. The simple geometric relationships between the fault and the fold allow for calibrating slip and erosion rates [Lavé and Avouac, 2000, Avouac, 2003], and for calculating uplift rates using geomorphological indicators [Bollinger et al., 2013]. It is also widely used to simulate the structural progression of these ramp-décollement structures [Hubbard et al., 2016]. However, the earthquake behavior of fault-bend folds and the dynamics of off-fault deformation along axial surfaces have not been investigated thus far. Therefore, there is a need to develop new and im- proved geophysical modeling frameworks that take into account the fault structural complexities and the permanent deformation associated with folding processes. 1.3 A new geophysical modeling framework for earthquake sequence simulationswithpermanentdeformation The evolution of fault slip from the time the fault plane accumulates stress and strain until its culmination as an earthquake to release this accumulated stress is defined as the seismic cycle [Avouac, 2015, Barbot et al., 2012a, Perfettini and Avouac, 2014, Van Dinther et al., 2013, Dal Zilio et al., 2019, Shi et al., 2020a]. The series of events in this cycle is termed interseismic (the quiescent period of strain accumulation), co-seismic (the rapid release of stress during a seismic instability or earthquake), and post-seismic (the accelerated response of the crust and the mantle to the earthquake). These processes embody physical mechanisms of different time and length scales and generate distinct geodetic and seismological signals. The seismic cycle is touted as the short-term deformational process spanning time scales of a few seconds to a few hundred years or millennia. 9 To simulate the dynamics of fault behavior during the seismic cycle and reproduce these cyclic events during the different phases, geophysical numerical models based on laboratory-based empirical relation- ships are often employed. The rate- and state-dependent friction framework allows us to simulate the evolution of fault strength and fault slip, and study the inherent stress build-up and release mechanisms in the crust [Dieterich, 1979, Ruina, 1983, Rice and Ruina, 1983]. These models have helped us understand a plethora of geophysical concepts, from the slow slip and slow earthquake phenomenon [Li and Liu, 2016, Liu and Rice, 2005, 2007, Rubin, 2008, Matsuzawa et al., 2010, Segall et al., 2010, Wei et al., 2013, 2015, Li and Liu, 2017, Goswami and Barbot, 2018, Wei et al., 2018, Romanet et al., 2018, Barbot, 2019a, Nie and Barbot, 2021] to the dynamic rupture process [Okubo, 1989, Lapusta et al., 2000a, Bizzarri and Cocco, 2003, Tinti et al., 2005, Gabriel et al., 2012a]. However, typical seismic cycle models based on rate- and state-dependent friction use simplified, smooth fault planes and often omit the structural and morphological variability that impacts several aspects of the seismic cycle [Rice and Tse, 1986, Lapusta and Rice, 2003, Hori et al., 2004, Barbot et al., 2012a, Mele Veedu and Barbot, 2016, Jiang and Lapusta, 2016]. Even though fault-bend folding and other kinematic models are widely used to describe several long-term tectonic models, they are not in- corporated into models that investigate short-term fault dynamics. To fill this gap, we develop geophysical models that can describe on-fault elastic strain and off-fault plastic strain using a unified framework. In this thesis, I combine geometry-based fault-bend folding theory with geophysical modeling frame- works using rate-and-state dependent friction. Fault structural features like fault bends and kinks are responsible for building topography in the long term and impacting the spatial and temporal characteris- tics of earthquake cycles in the short term. This sort of unified framework to describe the elastic faulting processes and inelastic folding processes in the crust can help us better assess the impact of faulting and folding on earthquake behavior and study the link between long-term geological processes and short-term earthquake cycle processes. It would also help us to better estimate seismic hazards because these more realistic models of earthquake sequences provide useful statistical data on earthquake behavior. 10 1.4 Layoutofthethesis In this thesis, I compile my research papers on the impact of fault structural complexities in the short and long term, that are either published or in review. In Chapter 2, my co-authors and I developed unified models of plasticity for faults and folds and provide a mathematical framework to deal with the mechanical interactions between faulting and folding. We develop two-dimensional, quasi-dynamic numerical simu- lations of slip evolution under the rate- and state-dependent friction framework. We find that earthquake behavior is sensitive to fault friction parameters, the magnitude of fault bending, and the cut-off angle between incoming sediments and the fault. Fault bends can play an important role in earthquake seg- mentation as a result of nonlinear fault dynamics, affecting the initiation and termination of earthquakes and the details of interseismic behavior. Shallow earthquakes that initiate, propagate, and terminate near the surface are facilitated when the stratigraphy within incoming thrust sheets is not parallel to the un- derlying fault, as this can change the loading rate across the fault bend. This study provides a first step towards creating more realistic rupture models that account for realistic geological structures and off-fault deformation. In Chapter 3, we delve into the Himalayan fault system to probe its structural complexities that played an important role in the 2015M w 7.8 Gorkha earthquake. This earthquake was an unfortunate reminder of the imminent seismic hazard hovering over the densely populated regions of Nepal due to the abutting of Indian and Eurasian tectonic plates along this collisional margin. The geometry of the subsurface fault structure in the source region of this earthquake remains debated. The fault system features multiple ramp- décollement pairs, but the exact details like the number and length of each of these segments are not well- known. Therefore, we use two end-member geometry models that are supported by various geophysical and geological datasets and observations, with and without a ramp in the middle of the seismogenic zone, and run a suite of physics-based numerical models to consistently explain the Gorkha earthquake sequence characteristics. Our models show that the ramp in the middle of the seismogenic zone is necessary to 11 elucidate the coseismic rupture extent and the postseismic behavior. This study helps us better understand the role of fault bends on the seismic cycle and helps better constrain Himalayan fault geometry at the source location of the Gorkha earthquake. The Himalayan megathrust forms and evolves over millions of years and undergoes a structural and behavioral metamorphosis as it continuously builds this behemoth mountain chain. These geological timescale processes interact with fault behavior, but unraveling the complex interplay between the long- term and short-term processes is challenging. In Chapter 3, I construct an elaborate two-dimensional structural evolution along two representative profiles in Central and East Nepal that embody different structural complexities using the fault-bend folding theory. These kinematic structural models help al- leviate the apparent discrepancy in fault motion rates in Central Nepal. Here, geologic investigations of abandoned Holocene fluvial terraces yield a convergence rate of 21± 1.5mm/year across the MHT, in con- trast with the much lower convergence rate of 15.5± 2mm/year from geodetic data. Our study shows that such discrepancies can be reconciled if we include the internal dynamics of the wedge that enable frontal thrusts to slip at different rates depending on the stage of structural evolution. These long-term processes have a profound impact on the recurrence times of earthquakes and other short-term seismic cycle characteristics. I conclude my thesis in Chapter 5 with a summary of important lessons that I amassed about the earthquake behavior of single fault-bend folds and large collisional zones with fault bends. I illustrate the impact of structural complexities on seismic cycle processes, and the role they play in mountain building. This dissertation contributes to an interdisciplinary approach, combining several fields of tectonics and earthquake physics, to develop advanced tools and methods, thereby, improving our ability to make robust seismic hazard and risk assessments. 12 1.5 DissertationPublications This dissertation contains the following published and in-review manuscripts: 1. Sathiakumar,S., Barbot, S. and Hubbard, J., 2020. Earthquake cycles in fault-bend folds. Journalof Geophysical Research: Solid Earth, 125(8), p.e2019JB018557. (Chapter2) 2. Sathiakumar,S. and Barbot, S., 2021. The stop-start control of seismicity by fault bends along the Main Himalayan Thrust.Communications Earth & Environment, 2(1), pp.1-11. (Chapter3) 3. Sathiakumar, S., Barbot, S. and Hubbard, J., 2022. Reconciling geologic and geodetic shortening rates with active folding in the Himalayan orogenic wedge, Science advances, submitted (Chapter 4) 13 Chapter2 Earthquakecyclesinfault-bendfolds This chapter has been published as Sathiakumar, S., Barbot, S. and Hubbard, J., 2020. Earthquake cycles in fault-bend folds. Journal of Geophysical Research: Solid Earth, 125(8), p.e2019JB018557. 2.1 Abstract In fold-and-thrust belts and accretionary prisms, fault bends induce folding in the hanging wall that can alter the long-term loading rate on the megathrust and profoundly influence earthquake-related processes. To understand the impact of non-planar faults and off-fault deformation on the seismic cycle, we incorpo- rate fault-bend fold theory into fault dynamics and develop two-dimensional numerical simulations of slip evolution under a physics-based rate- and state-dependent friction law. Fault bends can play an impor- tant role in earthquake segmentation as a result of nonlinear fault dynamics, affecting the initiation and termination of earthquakes and the details of long-term interseismic behavior. Shallow earthquakes that initiate, propagate, and terminate near the surface are facilitated when the stratigraphy within incoming thrust sheets is not parallel to the underlying fault, as this can change the loading rate across the fault bend. 14 2.2 Introduction Faulting and folding of the crust drive the growth of topography on Earth, including the high mountain ranges associated with fold-and-thrust belts in convergent margins [Suppe, 1981, Tapponnier et al., 1990, Cattin and Avouac, 2000, Avouac, 2003, Lacassin et al., 2004]. Folding mechanisms are also reported in ex- tensional [Howard and John, 1997, Sharp et al., 2000, Ferrill et al., 2012] and strike-slip environments [Tin- dall and Davis, 1999]. There is abundant evidence from geomorphology [Avouac et al., 1993, Lavé and Avouac, 2000, Shaw et al., 2004, Lai et al., 2006, Hossler et al., 2016] and seismic reflection profiles [Baker et al., 1988, Moore et al., 1990, Lai et al., 2006, Charreau et al., 2008, Singh et al., 2008] that fault-related folding is an important deformation mechanism at active plate boundaries, which has led to numerous kinematic analyses of folds and their structural evolution [Suppe, 1983, Sanderson, 1982, Brandes and Tan- ner, 2014, Betka et al., 2015, 2018] (Figure 2.1). The relationships between faults and folds have been used to infer sedimentation and uplift rates [Hurtrez et al., 1999, Lavé and Avouac, 2000, Hubbard and Shaw, 2009] and fault slip rates to understand long-term deformation [Suppe et al., 1992, Chen et al., 2007, Le Béon et al., 2014, Copley, 2014]. A variety of kinematic models of fault-related folds have been developed to explain the long-term crustal deformation inferred from seismic reflection profiles, topographic features, radiometric dating, well data, and seismicity [Epard and Groshong Jr, 1995, Poblet and McClay, 1996, Mueller and Suppe, 1997, Souter and Hager, 1997, van der Pluijm et al., 2001, Mitra, 2003, Miller and Slingerland, 2006, Bernard et al., 2007, Okamura et al., 2007, Wang et al., 2013, Forte et al., 2013]. Fault-bend fold theory [Suppe, 1981, Medwedeff and Suppe, 1997] provides a kinematic modeling framework to describe the deformation of rocks as they ride over non-planar fault surfaces [e.g., Rich, 1934, Rodgers, 1950, Avouac et al., 1992]. According to this theory, slip on a non-planar fault leads to hanging wall deformation, and the shape of the fold can be related to the shape of the fault by assuming conservation of cross-sectional area and bed length. These assumptions reflect observations from field relationships, as well as seismic reflection and 15 well data, and form the basis for a geometrical understanding of fault-related folding [Shaw et al., 2005]. One important take-away from these studies is that the long-term slip-rate on a fault can be variable in the down-dip direction when the fault bend cuts across the bedding in the folded sheet. This implies that the local slip-rates on a fault may differ from the average shortening rate due to folding in the upper plate, depending on the morphology of the fold. Although kinematic models of faults and folds trace the evolution of deformation over geological time scales, how the structural and geometric changes influence earthquake cycles is still unclear. Fault bends and their associated off-fault deformation may play a role in the segmentation of megath- rust seismicity and the so-called supercycles [e.g., Sieh et al., 2008, Philibosian et al., 2017]. The fact that some events break a single along-strike fault segment or several at once has been attributed to fric- tional heterogeneity [Ariyoshia et al., 2009, Kaneko et al., 2010, Dublanchet et al., 2013, Noda and Lapusta, 2013], morphological gradients [Duan and Oglesby, 2005, Lozos et al., 2011, Qiu et al., 2016b], fault step- overs [Klinger et al., 2006], nonlinear dynamics [Michel et al., 2017, Wu and Chen, 2014, Noda and Hori, 2014, Barbot, 2019b], or a combination of these factors [Duan and Oglesby, 2006, Mitsui, 2018]. The me- chanical conditions that lead to a wide range of earthquake sizes along a plate boundary and the role of fault geometry in earthquake segmentation remain poorly understood. Fault bends are relevant to earthquake initiation and termination processes [King and Nábělek, 1985, King, 1986], and studies of dynamic earth- quake ruptures with geometric complexities have shown significant off-fault deformation [Aochi et al., 2002, Kase and Day, 2006, Duan and Day, 2008, Zhang et al., 2014, Preuss et al., 2019, van Zelst et al., 2019]. Previous work on the dynamics of ruptures at fault bends is often limited to single-cycle ruptures [Poliakov et al., 2002a, Lozos et al., 2011, Aochi et al., 2000, Rousseau and Rosakis, 2009, O’Reilly et al., 2015], but a multi-cycle approach is required to discuss rupture dynamics and supercycles [e.g., Herrendörfer et al., 2015, Ong et al., 2019, Dal Zilio et al., 2019]. 16 Active plate boundaries involve a network of interconnected faults that together accommodate the far-field relative motion between two or more plates. Secondary faults may develop off the main strand because of the possible misalignment between the main fault and the relative plate motion [Maerten et al., 2002, Mitra, 2002, Daout et al., 2016a]. These complexities can even include the formation and activation of multiple fault types in the same region. For instance, the Los Angeles and Ventura basins in southern California host both blind and surface-emergent active thrust faults, although the regional tectonics are dominated by strike-slip behavior [Davis and Namson, 1994, Leon et al., 2007, Pratt et al., 2002, Shaw and Suppe, 1996, Hubbard et al., 2014]. These examples demonstrate that secondary faults or fault systems can develop to accommodate kinematic boundary conditions from large-scale plate tectonics and, more fundamentally, to conserve mass within the system. Incorporating these non-trivial kinematic boundary conditions into models of fault dynamics may help us better understand earthquake behavior in accre- tionary prisms and in fold-and-thrust belts. 17 1 km c d Fold limb Axial surface Axial surface a b ~1 ft ~1 ft Fold limb Fold limb Cut-off angle Figure 2.1: (a) Fold scarp at the Main Boundary Thrust near Islamabad, Pakistan. The folded sediment forms a cut-off angle with the ramp (white arrows). Location: 73 ◦ 8’7.44"E, 33 ◦ 45’54.3”N. (b) Nearby fold scarp. The axial surface (dashed green line) bisects the fold limbs. Location: 73 ◦ 8’4.22"E, 33 ◦ 45’39.432”N. Photos taken by S. Barbot on 8/25/2019. (c) A seismic profile from the Hikurangi subduction zone, taken from line PEG09-25, showing a thrust ramp rising from a décollement, with the development of a synclinal fault-bend fold in the hanging wall. The axial surface bisects the fold. Red line shows décollement, white arrows show thrust ramp (fault plane reflectors), blue line shows schematic folded stratigraphy. Scale is approximately 1:1. The line was collected as part of the New Zealand Government’s PEG09 2D reconnais- sance seismic survey in the Pegasus Basin. (d) Photograph from Nepal, showing flexural slip. The photo shows sub-vertical bedding (seen as color banding) with a clastic dike (thin blocky gray line) cutting across the bedding. The dike is offset by 10s of cm of slip along the bedding planes. White arrows show slip, black arrows show terminations of the dike at the slip plane. This slip occurred to accommodate the folding of this stratigraphic section from its original orientation (horizontal) to its current orientation (vertical). The scale is shown with a geological hammer. Photo taken by J. Hubbard on 04/29/2017. Location: 24.74615 ◦ N, 82.83798 ◦ E. 18 Complex fault geometries and dynamics may also help to explain the mechanism behind tsunami earthquakes, those shallow events in submarine accretionary prisms that generate low frequency seis- mic waves and larger seafloor uplift and tsunami excitation than expected for their magnitude [Kanamori, 1972, Kanamori and Kikuchi, 1993]. The root cause for shallow frictional instabilities for these shallow ruptures have been attributed to frictional contrasts at the smectite-illite transitions [e.g., Bilek and Lay, 2002, Saffer and Marone, 2003, Ikari et al., 2013, Ikari and Kopf, 2017]. However, efficient seafloor uplift could be the result of slip on ramps and other fault bends [Behrmann and Kopf, 2001, Kopp, 2013, Hubbard et al., 2015, Bradley et al., 2019, Geersen, 2019]. The nature of down-dip segmentation between shallow earthquakes and seismogenic zone earthquakes and the type of earthquake patterns it produces is still poorly understood [Lay et al., 2013, Sladen and Trevisan, 2018, Martin et al., Scholz, 2014, Ammon et al., 2006a, Kozdon and Dunham, 2013]. While some ruptures are confined to the lower edge of the seismo- genic zone, the largest earthquakes break all the way to the surface [Lay et al., 2005, Ammon et al., 2005, Feldl and Bilham, 2006, Simons et al., 2012, Mugnier et al., 2013]. Other types of events, slow-slip and slow earthquakes, also take place in the shallow section [Vallee et al., 2013, Wallace et al., 2017, Toh et al., 2018, Nakano et al., 2018]. These observations not only challenge the simplistic view that shallow regions like accretionary prisms in subduction zones inhibit unstable slip [Hubbard et al., 2015], but also suggest that other factors like geometric complexities may play an important role. The aim of this study is to investi- gate the impact of fault bends and associated folding on earthquake cycles, and to perform a sensitivity test on geometric elements that influence different earthquake processes. The geometric elements stud- ied here include the change in fault dip and angularity of the fault cut-off. We use a two-dimensional, quasi-dynamic numerical model under rate-and-state friction that includes structural changes to capture the ways in which fault geometry and fault friction properties interact to shape the seismic cycle. One distinguishing feature of this study is that we explore and analyze a comprehensive model space of both friction and geometry. We study the effect of frictional properties on fault dynamics for a range of slip 19 behaviors controlled by the ratio of asperity size to nucleation size. In addition, we incorporate the axial surface associated with fault bending into our model to simulate earthquake cycles with off-fault deforma- tion around fault bends [Suppe, 1983, Shaw et al., 2005]. In the next section, we describe the kinematics of axial surfaces and how they are incorporated into our models. Subsequently, we simulate seismic cycles with non-planar faults based on the principles of fault-bend fold theory [Suppe, 1983, Shaw et al., 2005] to disentangle the effects of friction, fault bends, and loading rates. 20 2.3 Kinematicsoffault-bendfolds 2.3.1 Fault-bendfoldsandaxialsurfaces The morphology of overthrust terranes is often complicated by a long history of shortening that involves fault growth, multiple bends, and inactive features from previous stages of deformation [Medwedeff and Suppe, 1997, Hubbard et al., 2016]. For the sake of clarity, we consider a single fault bend consisting of a flat décollement and a shallow ramp rising to the surface (Figure 2.2). Without off-fault deformation, the slip around fault bends would create a space problem over long time scales because mass is not conserved in the system or the deformation is unrealistic [Suppe, 1983]. To avoid such outcomes, anelastic deformation must accumulate in the upper plate around the bend. This can happen by faulting, for example with a backthrust or a splay fault [Park et al., 2002], or by folding [Lavé and Avouac, 2000, Shaw et al., 2005]. Faulting in the upper plate creates a triple junction in a vertical cross-section, a geometric setup that is not stable over long time scales because the secondary fault that cuts the upper plate is advected along the ramp on the main fault [Cox and Hart, 1986]. Indeed, faults are Lagrangian objects, and they move along with the host rocks. Folding associated with fault bends is commonly observed at shallow confining depths, such as in the sedimentary muds that accumulate at accretionary prisms [Davis et al., 1983, Biju- Duval et al., 1982, Kopp and Kukowski, 2003, Saffer and Bekins, 2002] or fold-and-thrust belts [Suppe, 1981, Tapponnier et al., 1990, Cattin and Avouac, 2000, Avouac, 2003, Shaw et al., 2004, Lai et al., 2006]. The active axial surface that is tied to the fault bend is an Eulerian object that is not advected along with the thrust sheet, leading to a stable geometric configuration for long-term deformation. The folding of the upper plate in response to a step in décollement proceeds in distinct stages. These in- volve simultaneous crestal narrowing and uplift, followed by crestal broadening over tectonic timescales [Suppe, 1983, Shaw et al., 2005], leaving behind recognizable features in the landscape. In the case of a sin- gle décollement-ramp structure, the orientation of the long-term velocity across the ramp varies from 21 décollement-parallel to ramp-parallel, generating net horizontal shortening near the bend. This short- ening is accommodated by folding of the upper plate, where mass is conserved through the growth of topography (Figure 2.2). Closed odograph v B v c β d n v 1 v 2 v 3 Closed odograph v B v c d n v 1 v 2 v 3 Domain 1 v 1 Domain 2 v 2 θ Hanging wall Footwall θ: Cut-off angle SYNCLINE ANTICLINE Domain 1 v 1 Domain 2 v 2 γ γ Footwall Hanging wall γ γ Domain 2 v 2 Domain 1 v 1 Hanging wall Footwall Domain 1 v 1 Domain 2 v 2 θ Footwall Hanging wall θ: Cut-off angle γ γ γ γ v A v B v c v B v c v A v A v A β Θ v 1 = v 2 v 1 = v 2 v 1 < v 2 v 1 > v 2 a b c d e f φ φ φ φ φ Θ φ Figure 2.2: Synclinal (a) and anticlinal (b) fault-bend folds used in our model. The axial surface (black dotted line) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angleγ and ϕ represents the change in dip between the two fault segments. Here, the uniform color in the hanging wall represents uniform local slip rate (v 1 = v 2 ) on the fault segments represented by the red solid line. Synclinal (c) and anticlinal (d) fault-bend folds with angular fault cut-off angle θ . The axial surface (black dotted line) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angle γ . Here, the color of the hanging wall changes across the axial surface, showing the increase (v 1 <v 2 in the case of c) and decrease (v 1 >v 2 in the case of d) of local slip rate on the ramp compared to the décollement as a result of folding through the fault bend. (e, f) The three velocities are plotted together, forming a closed loop. 22 Fault-bend fold theory describes the deformation of hanging wall rocks as they ride over fault bends by making specific assumptions regarding how mass is conserved in cross-sectional two-dimensional models. Specifically, the model assumes incompressibility of rocks and hence, conservation of layer thickness, conservation of bedding length, and conservation of cross-section area of the overall upper plate over long time scales [Suppe, 1983]. As the hanging wall rocks move through a fault bend, an active axial surface forms, tied to the fault bend. Because layer thickness is preserved, this axis bisects the angle between the fold limbs. This folding is accommodated by flexural slip, i.e., small amounts of slip parallel to bedding planes, concentrated along the active axial surface [Suppe, 1983]. As slip progresses, the inactive axial surface moves away from the active axial surface. Fault bends can alter not just the orientation of the slip vector, but also its magnitude. When there is no angular fault cut-off and the beds are parallel to the underlying fault, there is no change in the magnitude of fault slip. However, in many cases of structural fault evolution, the bedding cuts the fault at an angle. The conservation of layer thickness assumption of fault-bend fold theory implies that the active axial surface forms at the bisector of the two limbs of the bedding plane. The theory predicts different long-term slip- rates on each fault segment, based on the cut-off angle θ of the approaching sediments, defined as the angle between the hanging wall stratigraphy and the fault segment before the fault bend, and the angleϕ formed by the fault bend (Figure 2.2). This effect can cause the long-term slip rate on the ramp to increase, in the case where the fault steepens like in Figure 2.2c, or decrease, in the case where the fault flattens like in Figure 2.2d, compared to the décollement. Long-term variations in slip rate may be important for the seismic cycle because the recurrence time of any slip instability is controlled to the first order by the long-term local loading rate [e.g., Barbot et al., 2012b, Lapusta and Barbot, 2012, Cattania and Segall, 2018, Barbot, 2019b]. Therefore, the ratio of the amplitude of long-term slip rates before and after the bend, R v = v 1 v 2 , (2.1) 23 wherev 1 andv 2 are the magnitudes of the long-term slip velocity vectors below and above the fault bend, respectively, will play an important role. Fault-bend fold theory provides a simple way to relate the fault shape to the fold shape, such that the slip ratio can be defined as a function of the overall geometry [Suppe, 1983] R v = sin(γ − θ ) sin(ϕ +γ − θ ) , (2.2) whereγ is the angle between the two fold limbs,θ is the cut-off angle between the fault and the beds below the fault bend, andϕ is the fault bend angle (Figure 2.2). Because of the above considerations, fault bends may trigger significant departures from the typical fault dynamics inferred for planar faults, because of the effects of the geometry on stress transfer, off-fault deformation, including along the active axial surface, and variations in long-term loading rates. Our goal is to incorporate the kinematic constraints from fault-bend fold theory into dynamic models of fault slip during the seismic cycle. The deformation and stress on the surrounding rocks caused by folding along the axial surface can be described if the anelastic strain-rate that accompanies folding can be identified [Barbot and Fialko, 2010a, Barbot et al., 2017, Barbot, 2018]. A remaining issue is how to ascribe a frictional or rheological relationship to capture the dynamics of flexural slip along the axial surface. In the following presentation, we first describe the kinematics of deformation, which will lead to an anelastic strain-rate representation of folding. The dynamic and rheological aspects are discussed in the following sections. 2.3.2 Short-termdeformationinfault-bendfolds To connect the long-term deformation predicted by fault-bend fold theory with the short-term deformation that occurs during the seismic cycle, we equate the long-term deformation to the so-called eigenstrain that describes the anelastic deformation in a visco-elasto-plastic material [Andrews, 1978, Nemat-Nasser and Hori, 1999, Nemat-Nasser, 2004, Barbot and Fialko, 2010b, Barbot et al., 2017, Barbot, 2018]. In addition, 24 we use the long-term deformation as the loading rate during the seismic cycle in a manner similar to the back-slip model. We consider a simple décollement-ramp structure in two dimensions, as shown in Figure 2.2. More realistic cases involving more fault bends and axial surfaces can be treated by linear superposition. Syn- clinal fault-bend folds are structures with concave-upward fault bend, and anticlinal fault-bend folds have concave-downward fault bends. In either case, the active axial surface bisects the two limbs around the bend in order to preserve layer thickness. Consider the velocitiesv 1 andv 2 before and after the bend in domains1 and2, respectively. The long-term, anelastic velocity field in the upper plate can be expressed as v(x)=H(n· x)(v 1 − v 2 )+v 2 , (2.3) where the vectorsv 1 andv 2 are assumed to be spatially uniform, the spatial variations being captured by the Heaviside functionH(x), n is the vector normal to the axial plane, andx is a coordinate in the hanging wall in an Eulerian reference frame. The active axial surface is an Eulerian object, i.e., it is not advected with the thrust sheet. In contrast, inactive axial surfaces are advected, but they do not generate internal deformation. Therefore, an Eulerian coordinate system is most applicable to represent the long- term velocity field. The anelastic strain-rate that takes place along the active axial surface results from the spatial velocity gradient tensor L(x)=∇⊗ v, (2.4) where∇ is the gradient operator and⊗ is the dyadic product. Combining (2.3) and (2.4), we can find a closed-form expression for the spatial velocity gradient L(x)=[n⊗ (v 1 − v 2 )]δ (n· x), (2.5) 25 whereδ (x) is the Dirac delta function, the derivative of the Heaviside function, indicating localization of strain along the axial surface, i.e., infinitely large across an infinitesimally narrow thickness, in a manner similar to faulting [Steketee, 1958a,b, Barbot and Fialko, 2010a,b]. The symmetric part of the spatial velocity gradient gives the anelastic strain-rate tensor across the axial surface [Barbot and Fialko, 2010b] ˙ ϵ = 1 2 L+L T , (2.6) where the superscript T means transposition. At this point, we do not need to consider infinitesimal or finite deformation separately because the spatial velocity gradient is always infinitesimal. Closed odograph v B v c d n v 3 v A v B v c v A φ Θ Closed odograph d n v 1 v 2 φ Θ v 1 v 3 v 2 Block A Block B Block C Block 1 Block 2 Block 3 Figure 2.3: Velocity vector diagrams of a syncline with zero cut-off angle to depict the relative displace- ments (velocities) across the fault and axial surface, assuming it is a fixed block. v A represents the velocity of block A with respect to block B, v B represents the velocity of block B with respect to block C, and v C represents the velocity of block C with respect to block A. If we move in clockwise direction instead, v 1 represents the velocity of block 1 with respect to block 2,v 3 represents the velocity of block 2 with respect to block 3, andv 2 represents the velocity of block 3 with respect to block 1. We now parameterize the deformation in a manner consistent with dislocation theory, using the usual conventions that describe faults in field and theoretical studies [e.g., Okada, 1992, Wang et al., 2003, Smith and Sandwell, 2004]. The change in slip direction across the fault bend creates a resultant, or residual, velocityv 3 to conserve mass in the system, assuming incompressible rocks. For any choice ofv 1 andv 2 , we require v 3 =v 1 − v 2 (2.7) 26 to close the kinematics and avoid unrealistic voids and overlaps over long time-scales. The odograph that arranges the three velocity vectors sequentially in a closed loop [e.g., Cox and Hart, 1986, Daout et al., 2016b,a] predicts the direction and amplitude of long-term motion of the particles crossing the axial surface as a function of the velocity on the ramp and the décollement (Figure 2.2, 2.3). We callb the Burger’s vector that represents the instantaneous velocity across the axial surface. For a long-term average, we have b=v 3 . We usev 3 as the instantaneous loading rate for an axial surface. However, the Burger’s vector will be allowed to vary dynamically during the seismic cycle, such thatv 3 represents the asymptotic average velocity. 27 n d b β SYNCLINE ANTICLINE a b c e d Décollement Ramp Décollement Ramp axial surface n n axial surface γ γ γ γ e 3 e 2 e 1 0 s d b n β t n t d e 2 e 3 e 1 0 s d n t Figure 2.4: Synclinal (a) and anticlinal (b) fault-bend folds in three-dimensions. The axial surface (green surface) bisects the two limbs of the bedding plane in the hanging wall with the inter-limb angle γ , and the fault segments are represented by red solid line. (c) The 3-D perspective view of an axial surface in the reference coordinate system (s,d,n) that forms an orthonormal right-handed system of coordinates, wheres,d, andn are the unit vectors in the strike, dip, and normal directions of the axial surface, andb is the Burger’s vector showing slip direction. (d) Planar view of the Burger’s vectorb that describes the direction of motion on the axial surface. The take-off angle β is defined between the dip vector and normal vector. (e) The traction vector (red arrow) on the axial surface and its decomposition into normal (orange arrow) and shear (blue arrow) components. Consider the two-dimensional reference coordinate system(d,n) tied to the fault, whered andn are the unit vectors in the dip and normal direction of the axial surface, respectively (Figure 2.4). The Burger’s vector can be written as b=v 3 (cosβ d+sinβ n) , (2.8) 28 whereβ is the take-off angle between the axial surface and the Burger’s vector and v 3 is the norm of the velocity vectorv 3 (Figure 2.2). In turn, the unit vectors can be described in the global coordinate system (e 2 ,e 3 ) as follows: n=+e 2 sinΘ − e 3 cosΘ , d=−e 2 cosΘ − e 3 sinΘ , (2.9) whereΘ is the dip angle of the axial surface,e 2 is horizontal, ande 3 is vertical pointing down. Combining (2.4), (2.6), (2.7), and (2.8), we obtain the anelastic strain-rate tensor in closed form: ˙ ϵ =bR, (2.10) whereb is the norm of the Burger’s vectorb, and R= 1 2 cosβ (n⊗ d+d⊗ n)+sinβ n⊗ n (2.11) is a unit eigenstrain tensor [Barbot and Fialko, 2010b]. When a particle crosses the axial surface along the directionn, it experiences an offset at the rate and direction b (Figure 2.4). Depending on the geometry of the fault bend, the componentn⊗ n of the eigenstrain can contribute to tensile or compressive stresses in the case of synclines and anticlines, respectively (Figure 2.2). In the next section, we use the eigenstrain formulation for axial surfaces to represent the dynamics of folding at fault bends, and numerically simulate earthquake cycles on faults with bends to investigate how this particular type of off-fault deformation impacts the seismic cycle. 2.3.3 Unifiedrepresentationofplasticityforfaultsandfolds Although active axial surfaces are planar deformation features that might resemble faults, there are fun- damental differences in their kinematics. Kink faults that form at fault bends and any other fault are 29 Lagrangian objects that are advected along with the surrounding rocks. In addition, fault slip is always parallel to the fault surface. Small compaction and dilatancy may occur, but there is no wholesale de- formation in the fault-normal direction. This important property allows virtually unbounded offsets to accumulate along the fault. In contrast, active axial surfaces are Eulerian objects. They are not advected with the displacement of the surrounding rocks. Rather, the rock moves through the axial surface and is folded in the process. New particles constantly cross the surface as the hanging wall shortens, typically with a non-trivial component of surface-perpendicular uniaxial anelastic strain. This allows an unbounded accumulation of anelastic strain along the axial surface in an Eulerian reference frame. However, in a La- grangian reference frame, the strain accumulation for a labeled particle is bounded. These differences may seem trivial from a mathematical point of view, but they manifest themselves in profoundly different ex- pressions in the field [Lavé and Avouac, 2000, Chen et al., 2007, Shaw and Suppe, 1994, Rowan and Linares, 2000]. Despite these differences, the deformation caused by faults and axial surfaces can be written in a similar form in the context of dislocation theory. Indeed, an axial surface exhibits simultaneous simple shear, reflected by the symmetric components n⊗ d and d⊗ n in the strain-rate tensor, and uniaxial strain, captured by the componentn⊗ n in equation (2.9). These two components develop at once, with a fixed ratio controlled by the take-off angle β . In addition, the eigenstrain is localized along the deformation surface for both axial surfaces and faults. A unified representation of slip along a fault and plastic strain along an axial surface can be derived using dislocation theory by considering a surface element with a linear combination of tensile and shear components. For faults, only shear is allowed, implyingβ = 0. In this case, the amplitude of eigenstrain describes the cumulative slip across the fault, andb represents the instantaneous slip velocity. For axial surfaces, the take-off angle β is controlled by fault-bend fold theory and the amplitude of the Burger’s vectorb represents the velocity of particles crossing the axial surface. The displacements and stress in 30 the surrounding rocks caused by these elementary surface elements can be calculated using the Green’s function for fault slip [Chinnery, 1961, Savage and Hastie, 1966, Sato and Matsu’ura, 1974, Okada, 1985, 1992, Nikkhoo and Walter, 2015] or volume elements with a small aspect ratio [Nozaki and Taya, 2001, Kuvshinov, 2008, Barbot et al., 2017, Barbot, 2018]. This generalization provides a mathematical framework to deal with the mechanical interactions among faults and axial surfaces. 2.4 Modeldescription We now explore the dynamics of fault-bend folds during the seismic cycle in order to understand the effects of the fault bend, folding of the thrust sheet in the hanging wall, and variable loading rates along the fault. We restrict our analysis to a single fault bend with a décollement-ramp structure, where the ramp subsequently intersects the free surface. We also limit our study to two-dimensional cross-sections to keep the calculation tractable and to facilitate the interpretation. 2.4.1 Constitutivebehaviorforfaultslip To capture fault dynamics over multiple seismic cycles, we adopt the rate-and-state friction framework [Di- eterich, 1979, Ruina, 1983, Rice and Ruina, 1983], a set of constitutive relationships that explain an abun- dance of laboratory observations [Dieterich, 1978, Marone et al., 1990, Marone and Kilgore, 1993, Blan- pied et al., 1995, Dieterich and Kilgore, 1994, 1996, Nakatani, 2001, Carpenter et al., 2016, Scuderi et al., 2017]. The rate-and-state constitutive framework is supported by various micro-physical models of fault strength [Chester, 1989, 1994, Sleep, 1995, 2006, Nakatani, 2001, Barbot, 2019c] and has been successful at explaining fault slip during all stages of the seismic cycle [Rice and Tse, 1986, Lapusta and Rice, 2003, Hori et al., 2004, Barbot et al., 2012b, Mele Veedu and Barbot, 2016, Qiu et al., 2016b, Jiang and Lapusta, 2016], including slow slip and slow earthquakes [Li and Liu, 2016, Liu and Rice, 2005, 2007, Rubin, 2008, Matsuzawa et al., 2010, Segall et al., 2010, Wei et al., 2013, 2015, Li and Liu, 2017, Goswami and Barbot, 31 2018, Wei et al., 2018], and quasi-static, dominantly aseismic processes [Perfettini and Avouac, 2004, Bar- bot et al., 2009, Bruhat et al., 2011, Rousset et al., 2012, Rollins et al., 2015, Feng et al., 2016, Masuti et al., 2016, Lui and Lapusta, 2016, Salman et al., 2017, Muto et al., 2019, Agata et al., 2019]. Since its inception, several formulations of rate-and-state friction have been proposed [Ruina, 1983, Linker and Dieterich, 1992, Rice and Ben-Zion, 1996, Kato and Tullis, 2001, Bizzarri, 2011, Nagata et al., 2012], based on assumptions about the number of state variables, the effect of temperature, and the be- havior at stationary contact. Here, we follow the model derived by [Barbot, 2019c], where fault strength depends explicitly on the area of contact A= µ 0 ¯σ χ ψV 0 L b µ 0 , (2.12) whereµ 0 is the coefficient of friction, ¯σ is the effective normal stress, χ is the indentation hardness,V 0 is a reference velocity,ψ is a state variable associated with contact aging,L is the characteristic weakening distance, andb is a power exponent. The real area of contactA represents a density per unit nominal fault surface area. The yield surface and the resulting sliding velocity depend on the area of contact, following V =V 0 τ Aχ µ 0 a exp − Q R 1 T − 1 T 0 , (2.13) wherea≪ 1 is a power exponent,Q is an activation energy,R is the universal gas constant, andT is the absolute temperature. Combining equations (2.12) and (2.13), we obtain the following constitutive form for rate-and-state friction V =V 0 τ µ 0 ¯σ µ 0 a ψV 0 L − b a exp − Q R 1 T − 1 T 0 , (2.14) 32 or equivalently, in terms of shear traction, τ =µ 0 ¯σ V V 0 a µ 0 ψV 0 L b µ 0 exp a µ 0 Q R 1 T − 1 T 0 . (2.15) To capture the evolution of the real area of contact with sliding, the state variable follows the evolution law ˙ ψ =exp − H R 1 T − 1 T 0 − ψV L , (2.16) where the first and second terms on the right-hand side represent healing and weakening, respectively. At steady-state, the friction reduces to µ 0 at velocity V 0 and temperature T 0 . In isothermal conditions, equation (2.16) becomes the aging law [Ruina, 1983]. 2.4.2 Constitutivebehaviorforaxialsurfaces To simulate the dynamics of axial surfaces, that is, the evolution of the corresponding anelastic strain as a function of space and time, we need to prescribe the rate of anelastic deformation as a function of stress, and potentially other state variables. However, the frictional behavior or rheological properties of axial surfaces are largely unknown. We therefore proceed conservatively by assuming an overall velocity-strengthening behavior, which prevents the development of instabilities along the axial surface. Following previous work [e.g., Hudleston and Lan, 1995, Lan and Hudleston, 1995, Fuchs et al., 2015], we assume that the deformation of the axial surface is visco-elasto-plastic, represented by a non-linear, power-law relationship between stress and strain-rate, as is common for rock deformation by dislocation creep, diffusion creep, pressure-solution creep, and other micromechanical processes relevant to crustal and mantle deformation. To simplify the numerical treatment we adopt a non-linear rate-strengthening rheology of the form ˙ ϵ =Aτ n R, (2.17) 33 where˙ ϵ is the strain-rate tensor, i.e., the instantaneous Burger’s vector isb=Aτ n ,A is a material property, n is a stress power exponent, and the unit tensorR is given by R= 1 2 cosβ (n⊗ d+d⊗ n)+sinβ n⊗ n, (2.18) where d and n are the unit vectors in the dip and normal directions of the axial surface, and τ is the effective stress along the axial surface. The effective stress simplifies to τ = R:σ =cosβ n· σ · d+sinβ n· σ · n, (2.19) whereσ is the Cauchy stress tensor and the symbols· and: represent the vector and tensor dot products, respectively. Equation (2.19) represents the projection of the stress tensor in the directions specified by the orientation of the axial surface and the instantaneous particle velocity. 2.4.3 Controlsonfaultdynamics The dynamics of fault slip during the seismic cycle can exhibit tremendous complexity, with spontaneously emerging slow-slip events, slow earthquakes, period-two cycles, and deterministic chaos, depending on the distribution of physical properties [Kato, 2003, Lapusta and Rice, 2003, Liu and Rice, 2005, Chen and Lapusta, 2009, Mitsui and Hirahara, 2011, Kato, 2014, Wu and Chen, 2014, Kato, 2016, Cattania and Segall, 2018]. For a simple fault system with a single velocity-weakening asperity, the dynamics are chiefly con- trolled by two non-dimensional parameters [Barbot, 2019c]. The Dieterich-Ruina-Rice number, R u = (b− a)¯σ G W L , (2.20) 34 where G is the rigidity and W is the size of the velocity-weakening region, controls the importance of non-local stress transfer and depends on the spatial distribution of the frictional parameters along the fault. NegativeR u values correspond to velocity-strengthening behavior; vanishing positive values correspond to stable weakening; and large positiveR u values to unstable weakening. TheR b number R b = b− a b , (2.21) controls the importance of the transient weakening and strengthening effects. As velocity-neutral behavior corresponds toR b = 0 and velocity-weakening to0 < R b < 1, theR b number quantifies the proximity to velocity-neutral. In controlled velocity step laboratory experiments, the difference (a− b) controls the static stress drop andb controls the dynamic stress drop that occurs after the direct strengthening effect at increased slip speed. Therefore, theR b number also controls the ratio of dynamic to static stress drops during rupture in some conditions, thereby affecting the style of the rupture. Using typical values ofa = 0.01 andb = 0.014 for granite, we expect the value ofR b to be 0.286 in the seismogenic zone [Blanpied et al., 1995]. At greater depths,R b approaches zero and becomes negative across the transition from velocity-weakening to velocity-strengthening, associated with increased tem- perature [Blanpied et al., 1995, Scholz, 1998]. However, the parameters that constitute theR u number are poorly known in nature. Indeed, there can be considerable variations of effective normal stress on a fault, depending on pore pressure and confining depth. The characteristic weakening distance may also vary by orders of magnitude, depending on the thickness of the fault core and the gouge layer. Therefore, theR u number represents an important degree of freedom to characterize fault behavior, as different values can be associated with different styles of coseismic rupture and different behaviors during the inter-seismic period. 35 The characteristics of the seismic cycle of a single asperity on a planar fault as a function ofR u and R b have been studied extensively [e.g., Kato, 2003, Lapusta and Rice, 2003, Mitsui and Hirahara, 2011, Wu and Chen, 2014, Barbot, 2019b, and references therein]. For1≤ R u ≤ 2, the cycle is periodic and includes short-term or long-term slow-slip events [e.g., Liu and Rice, 2005]. For R u > 17, the cycles include full and partial ruptures of the velocity-weakening region [e.g., Michel et al., 2017]. Other work has focused on seismic cycles on faults with interacting asperities [Dublanchet et al., 2013, Kato, 2016, Lui and Lapusta, 2016, Mitsui, 2018]. Here, we consider non-planar faults and explore the effects of friction, fault geometry, and spatially variable loading rates for a fault-bend-fold system. Fault-bend fold theory assumes sharp fault bends and angular kink folds. In most cases, this simpli- fication yields simple and useful geometric representations of natural faults. However, faults and folds frequently present a larger range of geometries that do not conform to the idealized model [Butler et al., 2018]. We do not either explore the effect of fault curvature, although it can play an important role in coseismic processes [Poliakov et al., 2002b, Duan and Day, 2008, Templeton et al., 2009]. The effect of curvature on fault dynamics depends on both the degree of curvature and the frictional regime [Ong et al., 2019]. Another limitation of this model is unbounded stress concentrations at the kink with accumulating slip. This is especially problematic when tracking dynamic changes of normal stresses in the system. This unphysical behavior occurs in models of non-planar faults with and without axial surfaces. However, in natural fault systems, some physical mechanisms may occur along the fault, including fluid migration, anelastic compaction, and mass transport by pressure-solution creep, to diffuse any large stress concen- tration. In our numerical simulations, we assume that these effects take place efficiently, even though they are not resolved explicitly. Accordingly, we bound the variations of normal stress near the fault bend. Another limitation of this model is the assumption of conservation of layer thickness and bedding length. This assumption does not accurately describe folding processes in nature, as small-scale distortions of the beds happen when they are folded. But the assumption is useful, especially in the context of large 36 fault systems considered in this thesis. The internal strain (strain that is responsible for altering layer thickness) is much smaller than the strain accumulated at the fold axis. The assumption tries to make use of the observation that changes in layer thickness are orders of magnitude smaller than the scale of folds considered here. Regional and local balanced cross-sections constructed based on fault-bend fold theory have successfully used these assumptions. However, alternative geometrical formulations that include changes in bedding thickness and length should be considered for relevant scenarios. 2.5 Dynamicsoffault-bendfoldswithactiveaxialsurfaces We now explore the dynamics that emerge in a simple fault-bend fold system. The numerical simulations are based on the boundary integral method (see Appendix A). We use a single fault with two segments, composed of a lower segment rising onto a surface-emergent ramp. Slip on this structure generates a synclinal fault-bend fold (Figure 2.2). In general, décollements are weak structures that generally dip less than 10 ◦ . We use larger dip angles for the lower segment to explore a wide range of fault bends without affecting the interaction between the shallow ramp and the free surface, which would otherwise make the interpretation of the controlling factors more difficult. Admittedly, real faults are much more complex and may contain imbricate structures with more than one fault bend [e.g., Medwedeff and Suppe, 1997, Hubbard et al., 2016, and references therein]. We first consider the case where the loading rate along the fault is uniform, corresponding to a zero cut-off angle between the lower segment and the incoming sediments. This is the most common situation for this geometry, as décollements form parallel to stratigraphy, taking advantage of pre-existing weak bedding planes. This situation allows us to investigate the individual roles of friction and geometry. We devise several sets of experiments that map out the influence of friction and geometry on the styles of rupture that develop in this system. We first vary the frictional parameters of the fault, particularly the characteristic weakening distance L and hence, the controlling parameter R u for a fixed geometry. We 37 then vary the shape of the fold, i.e., the angleϕ between the two fault segments. Finally, we explore the effect of the initial cut-off angle θ between the sediment bedding and the lower segment of the fault. This situation can arise when the fault has additional bends down-dip, and there has been sufficient slip on the décollement to carry non-parallel beds up to and through the shallower fault bend being considered here. Alternatively, this geometry is observed when a fault exhibits a change in orientation, but neither fault segment is a décollement. Non-zero initial cut-off angles produce contrasting loading rates on the two fault segments. 2.5.1 Roleoffrictionontheseismiccycleofnon-planarfaults We investigate the combined effects of fault friction properties and the geometry of the fault bend on the dynamics of the seismic cycle for a simple synclinal fault-bend fold structure. We consider a fault-bend fold system with a lower and an upper fault segment. The fault is continuously velocity-weakening from the surface to a down-dip distanceW , and velocity-strengthening beyond that point. The velocity-weakening segment includes the ramp and a section of the lower segment. We include an active axial surface that extends from the surface to the fault bend. In all the following experiments, the geometry, the characteristic weakening distance, and the loading rate on the ramp are changed, but many parameters otherwise stay the same, as summarized in Table 2.1. First, we document how various frictional parameters shape the seismic cycle for a modest fault bend of ϕ = 10 ◦ . We explore the effect of different values of the R u number by changing the characteristic slip weakening distance L (Figure 2.5). For R u ≤ 2 which corresponds to large values of L, sliding is single-periodic with a peak velocity well below seismic speeds (Figure 2.5a and 2.5b ). Slow-slip events nucleate from the free surface and creep systematically penetrates far into the lower segment during the inter-event period. No seismic slip takes place for these lowR u numbers because they are associated with large nucleation sizes relative to the dimensions of the velocity-weakening region, which damp or even 38 Parameter Symbol Value Static friction coefficient µ 0 0.6 Effective normal stress ¯σ 100 MPa Reference velocity V 0 10 −6 m/s Steady-state rate dependence in velocity-strengthening region a− b +4× 10 −3 Steady-state rate dependence in velocity-weakening region a− b −4 × 10 −3 Down-dip width of the velocity-weakening region W 8000 m Poisson’s ratio ν 0.25 Simulations in Fig. 2.5 Fault bend magnitude ϕ 10 ◦ Sediment cut-off angle θ 0 ◦ Characteristic weakening distance L 0.006− 0.3 m Simulations in Fig. 2.6 Fault bend magnitude ϕ 10− 40 ◦ Sediment cut-off angle θ 0 ◦ Characteristic weakening distance L 0.006− 0.3 m Simulations in Fig. 2.7 Fault bend magnitude ϕ 40 ◦ Sediment cut-off angle θ 0− 40 ◦ Characteristic weakening distance L 0.006− 0.3 m Table 2.1: Frictional and geometric parameters characterizing a synclinal fault-bend fold. The static friction coefficient, effective normal stress, reference velocity, a and b frictional parameters, Poisson’s ratio, and the down-dip width of the velocity weakening region are the same for all experiments. The magnitude of fault bend, the sediment cut-off angle and the characteristic weakening distance are varied in some experiments. prevent any frictional instability. ForR u = 2.67,L = 10 cm, the main ruptures are firmly in the seismic range (Figure 2.5c). These ruptures initiate near the free surface and propagate unilaterally down the ramp and onto the lower segment. During the inter-seismic period, creep grows at the bottom of the ramp and also penetrates into the lower segment from the velocity-strengthening region. This creates a situation of non-stationary and laterally varying partial coupling. Despite these complex features, the cycle appears largely periodic, with both events and inter-events repeating from one cycle to the next. For5.33≤ R u ≤ 14.81, or18 ≤ L≤ 50 mm, through-going ruptures tend to nucleate near the bend and propagate bilaterally (Figure 2.5d-2.5i). For R u ≥ 6.67, the interseismic period features a complex sequence of slow-slip events. Deep creep propagates into the lower segment through a sequence of pulses. Slow-slip events also take place on the ramp, initiating near the bend. 39 Figure 2.5: Simulation results for a synclinal fault-bend fold, with a ramp dip of 40 ◦ and lower segment dip of 30 ◦ (change in dip ϕ = 10), for varying R u values. An illustration of the geometry and fault friction properties is shown in the top left panel. Velocity-weakening regions are shown in red and velocity-strengthening regions in black. The dotted black line indicates the axial surface, which is tied to the bend in the fault. (a)-(o) Seismic cycles: The color indicates the fault slip rate; red reflects seismic velocities and dark blue fault locking. Each panel shows about 2000 years of simulation for this geometry as the slip weakening distanceL is varied. The x-axis depicts time steps that are non-uniformly spaced, not time itself, with fast events having a finer resolution than during interseismic periods. As R u increases, we find bimodal seismicity, with partial and through-going ruptures. 40 As the value of R u increases, the system exhibits a range of different types of seismic cycles, from single-event periodic to more complex sequences of through-going and partial ruptures. ForR u ≤ 6.67 (Fig- ure 2.5e), simple through-going ruptures break the entire fault in a single pulse (see the cases in Fig- ure 2.5c-2.5e). ForR u =8.89 (Figure 2.5f), more complex through-going ruptures break the entire velocity- weakening region with two distinct pulses, with the first one limited to the lower segment region, whereas the second pulse breaks the ramp and the lower segment together. ForR u = 10.66 (Figure 2.5g, the two pulses form a closely spaced foreshock/mainshock sequence. Once again, the foreshock breaks only the lower segment, while the main shock initiates at the bend and breaks both the ramp and the lower segment. For increasing R u ≥ 14.81, the partial rupture becomes increasingly separated from the through-going rupture, and the cycle is better characterized by an alternating sequence of full and partial ruptures. How- ever, for sufficiently large R u ≥ 33.33, the through-going ruptures propagate through multiple pulses with short periods of acceleration and deceleration. The seismic cycle is not entirely periodic, with variations of hypocenter location from one through-going rupture to the next. In addition to seismic cycle complexity, earthquake segmentation is also controlled by theR u number. For R u ≤ 6.67, all events rupture the entire velocity-weakening region. Rupture segmentation into full and partial ruptures occurs forR u ≥ 14.81. For14.81≤ R u ≤ 26.66 (see Figures 2.5i-2.5o), fast partial ruptures are limited to the lower segment. However, for R u ≥ 33.33, partial ruptures are followed by aftershocks that rupture a small area near the bend. As R u increases, the number of events per cycle increases, i.e., the number of partial events increases. We do not explore the fault dynamics for larger valuesR u ≥ 45 due to the large computational burden. However, given insights from previous studies [Lapusta and Rice, 2003, Wu and Chen, 2014, Barbot, 2019c], it is reasonable to expect that the seismic cycle will become increasingly complex, with the emergence of a wider spectrum of rupture sizes. 41 2.5.2 Roleofthefaultbend We now evaluate the effect of the magnitude of the fault bend on the seismic cycle (Figure 2.6). We fix the dip angle of the ramp to40 ◦ . Through different independent numerical simulations, we explore different dip angles for the lower segment fromϕ =0 ◦ toϕ =30 ◦ in increments of 10 ◦ . In each case, the dip of the axial surface is adjusted to bisect the two limbs of the fold. In this case, the fold limbs have the same dip as the fault segments, as the cut-off angle is zero. The resulting four geometries have varying magnitudes of fault bends, representing geometries that are commonly observed within the accretionary prisms of subduction zones and continental fold-and-thrust belts [Suppe, 1983, Lavé and Avouac, 2000, McQuarrie, 2004, Moore et al., 1990, Grando and McClay, 2007]. With these geometries, we evaluate fault behavior as the angle between the two fault segments varies fromϕ =0 toϕ =40 ◦ . For the whole range of modeled fault bends, we find that down-dip earthquakes are segmented for higher values ofR u , whereas through-going ruptures occur for lower values ofR u (Figure 2.6). However, for the same value ofR u , rupture styles change depending on the magnitude of the fault bend. ForR u = 6.67, faults with gentle bends, for exampleϕ = 10 ◦ andϕ = 20 ◦ (Figures 2.6a and 2.6b), exhibit single- pulse through-going ruptures, whereas faults with larger bends (ϕ ≥ 30 ◦ ) produce multi-pulse through- going ruptures (Figures 2.6c and 2.6d). The model with the largest bend (ϕ =40) has a distinct foreshock- like pulse before the through-going rupture. We do not find this behavior for other geometries with the sameR u number (Figure 2.6d). For R u = 17.7, a range of slip behavior is observed depending on the fault bend. For ϕ = 10 ◦ , fault segmentation is characterized by a sequence of small partial ruptures that break the lower segment, followed by a through-going rupture, forming a periodic two-event cycle. All ruptures initiate at the bottom of the velocity-weakening region. Forϕ = 20 ◦ and30 ◦ , the partial ruptures represent foreshocks of the through-going rupture. The through-going rupture initiates at the bend where the foreshock rupture 42 Figure 2.6: Simulation results for the four geometries used in our study (ϕ = 10 ◦ ,ϕ = 20 ◦ ,ϕ = 30 ◦ ,ϕ = 40 ◦ ) for variousR u values. Geometry and fault friction properties used in the model are shown in the first row. Velocity-weakening regions are shown in red and velocity-strengthening regions in black. The dotted black line indicates the axial surface, which is tied to the fault bend. (a)-(l) Seismic cycles: The color indicates the fault slip rate; red reflects seismic velocities and dark blue, fault locking. Each panel represents about 2,000 years of simulation for a given geometry. The slip weakening distanceL is varied, showing changes in earthquake slip behavior as the magnitude of the fault bend increases. 43 just terminated (Figures 2.6f and 2.6g). For the largest bend of ϕ = 40 ◦ , the cycle is single-periodic and through-going ruptures exhibit a complex pattern of acceleration and deceleration (Figure 2.6h). ForR u = 44.44, the cycles include a wide range of rupture sizes. However, the pattern also evolves with the magnitude of the fault bend. Forϕ = 10 ◦ , the through-going ruptures are separated by partial ruptures confined to the lower segment (Figure 2.6i). For ϕ = 20 ◦ , in addition to the partial ruptures in the lower segment, there are aftershock-like small events immediately following the partial rupture that nucleate at the fault bend, propagating up-dip through the ramp (Figure 2.6j). Faults with larger bends (ϕ ≥ 30 ◦ ) produce partial ruptures on both the ramp and lower segment. Depending on the angle of the fault bend, one, two, or three partial ruptures occur between two through-going ruptures. In these cases, the hypocenters of small and large ruptures are closely associated with either the down-dip limit of the velocity-weakening region or the fault bend. Forϕ = 40 ◦ (Figure 2.6l), partial ruptures on the ramp are more frequent and significantly larger than in the case of ϕ = 30 ◦ (Figure 2.6k). However, these ramp ruptures never break the surface. In contrast to the seismic processes, the pattern of interseismic creep seems relatively unaffected by the magnitude of the fault bend. ForR u = 6.67, creep propagates into the velocity-weakening section of the lower segment during the seismic cycle, and another smaller stripe of creep concentrates near the fault bend, a pattern observed for all modeled geometries. ForR u =17.7, two creep regions also spread laterally during the interseismic period, but a through-going rupture typically takes place before creep covers the entire lower segment. ForR u =44.44, the down-dip boundary of high coupling closely coincides with the down-dip boundary of the velocity-weakening region. However, creep still takes place around the fault bend. Interseismic creep seems to always be associated with the fault bend. 44 Figure 2.7: Simulation results for the model (ϕ = 40 ◦ , andθ = 10 ◦ ,θ = 20 ◦ ,θ = 30 ◦ ,θ = 40 ◦ ) for variousR u values. The geometry and fault friction properties are shown in the first row. VW regions are shown in red and VS regions in black. The dotted black line indicates the axial surface, also denoting the position of the kink. Panels A to O) Seismic cycles: The color indicates the fault slip rate: red reflects seismic velocities and dark blue fault, locking. Each panel represents 2,000 years of simulation for a given geometry. The slip weakening distanceL is varied, showing changes in earthquake slip behavior as the cut-off angle increases. 45 2.5.3 Roleofthesedimentcut-offangle In this section, we investigate the changes in cut-off angle θ and the associated changes in fault loading rates on the seismic cycle. As slip progresses along multiple fault bends, non-parallel beds can be brought up to shallower regions or faults break at an angle to stratigraphy. In both cases, the local loading rate on the fault changes depending on the cut-off angle of the approaching sediments; this is a kinematic effect of area conservation (Figure 2.2). This situation is common in fold and thrust belts [Suppe, 1983, Shaw et al., 2005, Hubbard et al., 2016], and we suggest that it may be accompanied by systematic changes in local fault dynamics. We first discuss the simulation results for a large fault-bend fold (ϕ = 40 ◦ ) to understand the range of fault behavior as θ varies (Figure 2.7). For this geometry, the through-going ruptures that occur for small values ofR u break the entire fault when the incoming sediments are gently sloping (see results from Sections 2.5.1 and 2.5.2). As the incoming limbs steepen, partial ruptures alternate breaking the ramp and décollement (i.e., the lower segment), but no single event ruptures the entire fault. This is a significant deviation from the rupture behavior observed for a fault with identical frictional properties and geometry but a cut-off angle of zero (R u =5.3 andϕ =40 ◦ , butθ =0 ◦ ; see Figure 2.7). ForR u = 22.2, we observe a range of behaviors that are distinct for each value ofθ . Whenθ = 0 ◦ , we see the multi-pulse through-going ruptures that are characteristic of large bends (Figure 2.7f). For low cut-off angles (θ =10 ◦ orθ =20 ◦ ), a small foreshock on the lower segment precedes a large rupture that breaks the entire velocity-weakening region (Figures 2.7g and 2.7h). Forθ =30 ◦ , through-going ruptures break the fault in multiple pulses and decelerate on the ramp segment (Figure 2.7i). Forθ = 40 ◦ , we see alternating partial ruptures (similar to the previous case ofR u =5.33 andθ =40 ◦ , Figure 2.7e), but these ruptures are faster than forR u =5.3 (brighter red in Figure 2.7j). For higher values of R u = 44.4, a wide range of ruptures break different regions of the fault, as ex- pected. Earthquake segmentation is observed for all cases ofθ > 0 ◦ , consistent with previous simulations 46 of planar and non-planar faults with θ = 0 ◦ . As the recurrence time is controlled by the local slip rate, which itself is controlled by the geometry of the fault and cut-off angle, significant changes in the earth- quake cycles take place as the cut-off angle increases. We observe a complex earthquake sequence with a range of partial events that break several parts of the fault before culminating in a through-going rupture. The simulation results for other fault bends (not shown) reveal similar trends. For all values of ϕ, the complexity of earthquake behavior increases as θ increases. Therefore, θ and the associated ratio of loading velocity R v may have a profound effect on the down-dip segmentation of a megathrust and the rupture style during the seismic cycle. Also, as the magnitude of the bend increases, the effects of θ and R v become more prominent. 47 2.6 Discussion Our simulations reveal several key aspects related to the dynamics of fault-bend folds. The model predicts a wide range of slip behaviors based on both frictional properties and geometric complexities. The map- ping between frictional and geometric variations and the corresponding fault behavior is useful to better understand earthquake source mechanisms. 2.6.1 Theeffectoffrictionandfaultbendmagnitudeontheseismiccycle The relation between frictional properties, captured by the Dieterich-Ruina-Rice numberR u , and the mag- nitude of fault bendingϕ can be visualized in theR u − ϕ parameter space to understand the dynamics of fault-bend folds (Figure 2.8). TheR u numbers can be associated with specific styles of rupture, producing a range of slip behaviors, from short and long-term slow events (forR u <2), to periodic through-going rup- tures (for2<R u <17), and complex earthquake cycles with through-going and partial events (R u >17). This mapping betweenR u and fault behavior is mostly based on simulations of planar faults. Our study reveals that even thoughR u is one of the dominant controls on long-term fault behavior, significant varia- tions in fault dynamics arise from geometric complexities due to fault bends. We attribute these variations in fault dynamics to stress heterogeneities due to fault-bend folding. First, we find that gentle fault bends produce more characteristic through-going ruptures in R u space than larger bends, which exhibit more complicated ruptures for the same value ofR u (see Figure 2.8 for R u = 5.33,6.67). While the bend slows down the initial rupture, the additional stress arising from the bend and the termination point of the initial slip phase triggers a larger through-going event that eventually ruptures the entire fault. With gentle fault bends, stress at the kink does not act as a barrier, and ruptures can extend through the entire velocity-weakening region. Therefore, fault bends not only play a role in coseismic processes, like initiation and termination of earthquakes [King and Nábělek, 1985, King, 1986], they also shape the details of the long-term fault behavior. As a result, for the same frictional properties 48 andR u number, the rupture style can vary depending on the magnitude of the fault bend. The persistence of seismicity at fault bends and other structural gradients may be explained by these effects [Lavé and Avouac, 2000, Wang et al., 2017b, 2018]. Shear stress variations on a fault can lead to earthquake segmentation processes and a diversity of earthquake sizes [Wu and Chen, 2014, Herrendörfer et al., 2015, Cattania and Segall, 2018, Dal Zilio et al., 2019, Barbot, 2019b]. For planar faults, partial ruptures may occur due to regions of low stress that can slow or stop rupture propagation. Low-stress barriers may form when stress accumulation is non-uniform, for example, when the velocity-weakening region is loaded from distant, creeping fault sections. Relatively small nucleation size compared to the velocity-weakening region, or, equivalently, large R u numbers, facilitates non-uniform stress accumulation, leading to earthquake segmentation. The recurring pattern of bimodal seismicity with full and partial ruptures emerges from the cumulative stress transfers of previous small events that eventually allow the velocity-weakening region to fail completely, producing a through- going rupture. Our findings qualitatively agree with frictional properties significantly affecting earthquake segmenta- tion processes on non-planar faults, as we observe segmentation only forR u > 10 in this case. However, the R u value at which faults start to show partial ruptures depend on the magnitude of the fault bend. Gentle fault bends show segmentation at lower values (R u = 14 for ϕ = 10 ◦ ) than large fault bends, which show segmentation at much higher values (R u =26 forϕ =40 ◦ ). This is because large fault bends have an overall higher shear stress on the fault and can create through-going ruptures, in contrast to gentle fault bends that have low stress barriers for the same value ofR u (see the cases in Figures 2.6e and 2.6h). Planar faults do not exhibit significant changes in earthquake segmentation with respect to R u when the dip of the fault is varied. This shows the importance of geometric complexities in models of the seismic cycles. 49 Rv=1.0 (Bedding parallel) Large bends No earthquake segmentation for both gentle and sharp bends Multi-pulse rupture Characteristic ruptures Through-going giant ruptures Through-going giant ruptures 10 20 30 40 0 Change in dip ( φ) Partial rupture on décollement Small events on ramp Partial rupture on ramp Partial rupture at the base of VW region No shallow ruptures on planar faults Earthquake segmentation for both gentle and sharp bends 2 5 6 8 10 13 14 17 22 26 33 38 44 Dieterich-Ruina-Rice number (Ru) Characteristic ruptures Gentle bends Planar fault (dip=40) Earthquake segmentation Fault Foreshocks/Aftershocks Approximate location of earthquake initiation Characteristic Multi-pulse Through-going giant ruptures Cycles with both characteristic and multi-pulse Ramp ruptures Surface-breaking ruptures Partial ruptures Décollement ruptures Figure 2.8: The effect of R u and change in fault dip (ϕ) on slip behavior of non-planar faults with one fault bend and planar faults with no fault bends (ϕ=0 ◦ ). Asϕ increases, the rupture complexity increases. Rupture segmentation depends onϕ forR u values greater than 10. For smaller values ofR u , multi-pulse ruptures are more common for large fault bends than for gentle fault bends. In the case of planar faults, there are no multi-pulse or partial ruptures in the shallow parts of the fault. 50 Finally, we find that all events on a planar fault originate from the deeper sections of the velocity- weakening region, compatible with previous studies. This is because the shear stress accumulation is highest near the boundary of the velocity-weakening region [Liu and Rice, 2005, 2007, Cattania and Segall, 2018]. In contrast, for non-planar faults, ruptures can initiate at various locations along the fault, depending on the state of shear stress on the fault, which depends largely on the magnitude of the fault bend ϕ. Our simulations also show that as ϕ increases, the size and frequency of ramp ruptures also increase, presumably due to gradients of stress accumulation at the bend being controlled by the magnitude of the fault bend. Therefore, larger fault bends are more likely to produce shallow events than gentle bends. This is in contrast to planar faults or faults with gentle bends, which always show partial ruptures near the transition boundaries between the velocity-weakening and velocity-strengthening regions. These ruptures never break the ramp segment nor reach the surface (Figures 2.5, 2.6, and 2.8). As seen in Figure 2.8, large bends can produce shallow ruptures due to the additional stress heterogeneity near the bend. This may be important when studying megathrust earthquakes, both on-land and offshore. Fault structures with significant fault bends may be more susceptible to shallow ruptures. These findings provide a theoretical context to understand the emergence of shallow ruptures at continental fold-and-thrust belts like Nepal and Taiwan, and at subduction zones around the world [Lavé and Avouac, 2000, Yue et al., 2005, Chen et al., 2007, Seeber et al., 2007, Hubbard et al., 2015, 2016, Bradley et al., 2019]. In summary, Figure 2.8 highlights the different dynamics of non-planar and planar faults, demonstrat- ing the importance of understanding the geometric complexities of natural faults. TheR u number for a given fault segment is poorly constrained in nature, even though observations of fault behavior during all phases of the seismic cycle may eventually help place useful bounds. However, our study also reveals that a singleR u number may not be sufficient to characterize the entire fault behavior when structural complex- ities are involved. Proper quantification of the geometry and loading rates are required to better predict the fault behavior. While real-world megathrusts may not be directly mapped to this simplified phase 51 space, our models emphasize that frictional and structural complexities both strongly control earthquake behavior. 2.6.2 Theeffectofvariableloadingrateontheseismiccycle When the angleθ formed between the stratigraphy and the fault below the fault bend is non-zero, there are significant changes in fault dynamics. Figure 2.9 summarizes the effects of change in dip ϕ and the hanging wall cut-off angle θ on slip behavior for a single fault-bend fold system. We will discuss two cases ofR u corresponding to low and high values (R u =6.6 andR u =44.4). We first summarize the behavior for R u =6.6. Forθ =0 ◦ , we usually observe through-going ruptures that are either characteristic or complicated, as seen in Figures 2.6 and 2.8 depending on the value of ϕ. Small foreshocks on the décollement preceding a larger through-going rupture are seen for ϕ = 40 ◦ . However, there is no earthquake segmentation for either gentle and large fault bends. As the value of θ increases, the ramp accrues slip faster than the décollement according to fault-bend fold theory. The ratio between slip on the ramp and the décollement, given byR v , is constrained by the geometry of the system [Suppe, 1983]. Gentle fault bends do not show significant changes in rupture dynamics with respect toθ . This is mainly because large values ofθ for gentle fault bends with low values ofϕ correspond to low values ofR v (refer equations 2.1 and 2.2), making the effect of θ minimal. Ruptures initiate at shallower depths, with an approximate 10% increase in slip on the ramp for a fault with ϕ = 10. Shallow rupture initiation can be attributed to heterogeneous stress accumulation on the ramp, which is slipping at a higher rate than the décollement. A 20% increase in slip on the ramp leads to more complicated through-going rupture for ϕ = 20 ◦ , in contrast to the characteristic through-going ruptures observed in the case of bedding-parallel thrust sheets with a gentle fault bend. The effect of θ increases when the value of ϕ increases, as seen in the case of ϕ = 30 ◦ and ϕ = 40 ◦ . For ϕ = 40 ◦ and θ = 40 ◦ , which correspond to R v = 1.6, there is earthquake segmentation, with partial ruptures that break the two segments of 52 the fault individually, but there are no through-going earthquakes. This behavior is unique among all the simulations considered, and shows that down-dip earthquake segmentation can be affected by the properties of the hanging wall. More specifically, these properties refer to the steeply dipping incoming thrust sheets that cause the ramp to slip at a higher rate than the décollement in areas where active folding of thrust sheets is common [Lavé and Avouac, 2000, Tapponnier et al., 1990, Medwedeff and Suppe, 1997, Wang et al., 2013]. In the case of R u = 44.4, complex earthquake cycles with full and partial ruptures are common in all cases considered. From the phase diagram, we identify two kinds of shallow ruptures: partial ramp ruptures that do not reach the surface, and “surface-breaking” ruptures. In nature, both surface-breaking and blind ruptures are capable of creating vertical seafloor displacements that can lead to the generation of large tsunamis [Nanayama et al., 2003, Fujiwara et al., 2011]). Fault bends withϕ = 10 ◦ show no shallow ruptures even with very steep incoming thrust sheets of θ = 30 ◦ or θ = 40 ◦ . In the case of ϕ = 20 ◦ andθ =30 ◦ , small foreshock/aftershock-like events as well as individual ramp ruptures can be observed. These shallow ramp ruptures did not occur for the bedding-parallel case, indicating that higher slip rates on ramps due to steep incoming thrust sheets can promote shallow slip. This effect is most pronounced when the fault bend is large, withϕ =40 ◦ . When the hanging wall cut-off angle is greater than zero for a large fault bend ( ϕ = 40 ◦ ), the ramp, which is accumulating slip more rapidly, hosts more frequent events (Figure 2.10). Non-zero angular fault cut-offs are possible in areas of significant shortening and multiple bends. In these areas, even gently dipping incoming thrust sheets with large fault bends can result in surface-breaking shallow ramp rup- tures. If these shallow ruptures produce significant seafloor displacements, they could lead to destructive tsunamis [Kanamori, 1972, Kanamori and Kikuchi, 1993, Polet and Kanamori, 2000, Ammon et al., 2006b]. Even though the characteristic properties and dynamics of accretionary prisms in subduction zones are 53 Fault Foreshocks/Aftershocks Approximate location of earthquake initiation Characteristic Multi-pulse Through-going giant ruptures Cycles with both characteristic and multi-pulse Ramp ruptures Surface-breaking ruptures Partial ruptures Décollement ruptures 10 20 30 40 0 Cut-off angle ( θ) Ru = 6.6 Gentle bends Large bends Change in dip ( φ) 10 20 30 40 Rv=1.0 (Bedding parallel) Earthquake segmentation Earthquake doublets Surface breaking shallow rupture Shallow rupture initiation Earthquake segmentation Characteristic through-going giant ruptures Multi-pulse through-going giant ruptures Small foreshocks in the shallow region No earthquake segmentation for both gentle and sharp bends Multi-pulse ruptures Multi-pulse ruptures Small foreshocks in the deeper region Gentle incoming thrust sheets Steep incoming thrust sheets Ru = 44.4 Gentle bends Large bends Ramp ruptures More ramp ruptures Trench-breaking ruptures No surface breaking shallow ruptures in gentle bends No shallow ruptures Change in dip ( φ) 10 20 30 40 Surface breaking shallow rupture No shallow ruptures Shallow rupture More shallow ruptures Shallow rupture Rv = 1.1 Rv = 1.05 Rv = 1.2 Rv = 1.3 Rv = 1.4 Rv = 1.5 Rv = 1.6 Rv = 1.7 Rv = 1.1 Rv = 1.05 Rv = 1.2 Rv = 1.3 Rv = 1.4 Rv = 1.5 Rv = 1.6 Rv = 1.7 Rv=1.0 (Bedding parallel) Alternating rupture initiation Gentle incoming thrust sheets Steep incoming thrust sheets Figure 2.9: The effect of change in fault dip (ϕ) and cut-off angle (θ ) on the slip behavior of non-planar faults for uniform (θ =0 ◦ ) and variable (θ ̸=0 ◦ ) local slip rate on the two segments of the fault for a small value ofR u = 6.6 (left) and a large value ofR u = 44.4 (right). As the dip of the incoming thrust sheet increases (asθ increases), more multi-pulse ruptures are observed forR u =6.6. Similarly, steep incoming thrust sheets allow shallow partial ruptures in contrast to gentler incoming sheets that show smaller or no ruptures in the shallow region. Therefore, as θ increases, ramp ruptures for all values of ϕ become prevalent. Large cut-off angles produce surface-breaking shallow ramp ruptures. 54 largely unknown, studies of shallow earthquakes that break to or near the surface often suggest that fric- tional heterogeneties are the main underlying cause [Bilek and Lay, 2002, Saffer and Marone, 2003, Ikari et al., 2013, Ikari and Kopf, 2017]. Here, we show that deformation in the upper plate with active folding can lead to the initiation and propagation of shallow earthquakes that reach the trench. We highlight the fact that the structural evolution of faults during subduction can enhance near-surface slip and therefore seafloor uplift, and possibly contribute to tsunami generation. The case of R u = 44.4 may also be associated with the seismogenic regions of the megathrust that host a wide variety of ruptures that are spatio-temporally separated. Paleoseismic records have suggested that some faults produce earthquake supercycles, where the fault breaks in a piecewise fashion separated by spatio-temporal scales [Sieh et al., 2008, Philibosian et al., 2017]. Simple numerical models tend to produce simple earthquake cycles or supercycles. In Figure 2.11, we observe earthquakes of various sizes that break different regions of the fault, followed by a giant rupture that breaks the entire fault. However, the dynamics of non-planar faults are more complex, causing ruptures that do not fit into this supercycle model. Shear stress variations that arise due to both changes in dip and the resulting variations in long- term slip can lead to chaotic cycles that do not repeat periodically. Even faults with the same geometry can produce different rupture styles as they evolve, due to the changes in hanging wall deformation. For instance, faults that broke in through-going ruptures in the past may start to produce partial ruptures as the hanging wall cut-off angle changes (Figure 2.7). If this pattern holds for more realistic, 3D fault surfaces, it will be important to understand not just the geometry of a fault in three dimensions but also the hanging wall structure in order to forecast earthquake patterns. Although our model based on a single fault bend is not representative of a full megathrust, insights from this model can be applied to segments of a megathrust to understand fault dynamics. The knowledge of subsurface geometric features and the history of structural evolution can help better understand fault behavior as the fault evolves using numerical models. 55 Through-going rupture Down-dip distance (km) 2 4 6 Ramp partial rupture Surface-breaking partial rupture Ramp: more slip 100 110 120 130 Down-dip distance (km) 10 26 18 2 0 4 70 80 65 75 70 80 65 75 85 Décollement: less slip Axial surface: Off-fault deformation Aseismic slip seismic slip Through-going rupture Décollement partial rupture Aseismic slip Cumulative slip (m) Cumulative slip (m) a b c Ramp partial rupture Rupture nucleation φ = 40 θ = 40 Rv=1.6 Slip velocity (m/s) 10 -6 10 0 10 -12 Figure 2.10: Slip accumulation on the fault and axial surface plotted against the down-dip distance for a model withϕ = 40 ◦ ,θ = 40 ◦ , andR v = 1.6. Slip is much higher on the ramp than the décollement due to the high cut-off angle and large fault bend. The geometry of the fault and the simulation results can be found in Figure 2.7o. 56 1 16 31 33 36 10 Supercycle 0 2400 2600 2800 3000 3200 3400 3600 3800 4000 2 3 4 5 6 7 8 9 10 11 12 13 14 15 18 17 19 21 22 23 24 25 26 27 28 20 29 30 32 34 35 37 38 39 41 40 V threshold Velocity (m/s) Velocity (m/s) Cumulative force (N) 10 -8 10 4 10 -12 0 10 -8 10 4 10 -12 7 8 9 Full rupture Full rupture Ramp rupture Full rupture Full rupture Ramp rupture 1 3 4 6 8 10 12 2 5 7 9 11 13 15 16 18 17 19 21 22 23 24 25 26 27 28 29 30 16 14 20 Full rupture Full rupture 31 31 32 33 34 35 36 37 38 39 41 40 33 36 Décollement rupture 1 V threshold Time (year) Décollement rupture 30 2500 3000 3500 4000 20 10 Distance (km) a b c d Panel D Supercycle Panel C Time (year) local slip rate local slip rate Supercycle ? Aseismic slip Full rupture Décollement rupture Ramp rupture 20 mm/year 32 mm/year Surface-breaking Supercycle Kink Surface-breaking Surface-breaking 1 16 3133 36 Figure 2.11: (a). The 3-D view of a décollement-ramp structure showing the hanging wall deformation due to folding for the model with ϕ = 40 ◦ , θ = 40 ◦ , R v = 1.6. Slip on the fault is shown by the color: aseismic slip in gray, yellow lines indicate through-going ruptures, green lines are partial ruptures on the décollement, blue lines indicate partial ruptures on the ramp that do not break the surface, and pink lines are surface-breaking partial ruptures on the ramp. (b) Simulation results computed for 2000 years for the above geometry, showing through-going surface ruptures and partial ruptures with varying magnitudes (cumulative force in 2-D). Repeating rupture patterns analogous to supercycles are highlighted with green dashed boxes. (c, d). Maximum slip velocity plotted against time on the ramp (c) and décollement (d). The geometry of the fault and the simulation results can be found in Figure 2.7o. 57 2.7 Conclusions Subduction zones and fold-and-thrust belts host some of the biggest earthquakes, as well as a multitude of events that break different segments of the megathrust. Upper plate deformation in these margins is often characterized by active fault-related folding processes. Geological studies of faults and folds combined with numerical modeling of multiple seismic cycles can help improve our understanding of fault dynam- ics, including earthquake segmentation, off-fault plasticity, and the propagation of shallow ruptures. We perform a series of numerical experiments to study the behavior of non-planar faults with fault bends using the rate-and-state friction law. We systematically analyze fault slip behavior with fault-bend folds and investigate the effects of frictional and structural complexities on seismic cycles. We numerically in- corporate axial surfaces in our models to simulate the kinematics of off-fault deformation associated with folding. The non-dimensional rate-and-state frictional parameterR u plays an important role in controlling the style of rupture segmentation, and fault bends add complexity to the rupture styles. Our two-dimensional quasi-dynamic model compares the behavior of four fault geometries with varying degrees of bends to un- derstand how the change in dip affects earthquake rupture styles. We find that even though the underlying earthquake slip regime depends primarily on the frictional parameterR u , the presence of a bend changes the rupture behavior and introduces new features depending on the morphology of the fold. Gentler bends produce seismic cycles more similar to that of planar faults. SmallR u values produce through-going rup- tures for gentle fault bends and multi-pulse through-going ruptures for larger bends. The value ofR u at which a fault starts to show bimodal seismicity depends on the change in fault dip. Large bends show segmentation at higher values ofR u than for gentle bends. Therefore, the details of rupture segmentation depend on the change in dip between the two fault segments described by ϕ. Planar faults, i.e., when ϕ =0 ◦ , exhibit no shallow ruptures, in comparison to non-planar faults, which show infrequent and small shallow ramp ruptures in the case of gentle fault bends. As ϕ increases, the size and frequency of ramp 58 ruptures increase. However, there are no surface-breaking shallow events when the bedding is parallel to the fault for either small and large bends. If there is a non-zero cut-off angle between the hanging wall stratigraphy and the fault below the fault bend, the resulting change in the slip rate on the ramp can produce significant changes in rupture styles. In the case of synclinal fault-bend folds, the slip accumulated on the ramp will be higher than the décollement, resulting in heterogeneous loading rates. Irrespective of the frictional parameters, the stress heterogeneity caused by variable slip rates can alter the earthquake behavior. The effect of the cut-off angle is much more significant for large fault bends than for gentle fault bends. In most cases, the bend plays a dominant role in rupture initiation and termination processes. The shape of the fold, controlled byϕ, and the cutoff angle, θ , together influence the rupture style. We find that folding in the upper plate can initiate and propagate shallow ramp ruptures that break the surface. This behavior could enhance seafloor uplift for submarine faults and potentially produce more energetic tsunamis. The current study combines knowledge from geological observations and geophysical modeling to understand how fault behavior is affected by fault geometry and hanging wall deformation. In many cases, the structural evolution of faults is known from geological observations. We find that earthquake behavior is sensitive to fault friction parameters, the magnitude of fault bending, and the hanging wall cut-off angle with incoming sediments. This study provides a first step towards creating more realistic rupture models that account for geological structure and off-fault deformation. 2.8 Acknowledgments We thank Associate Editor Alice Gabriel, Brendan Meade, and an anonymous reviewer at Journal of Geophysical Research: Solid Earth for their constructive feedback that helped improve the quality of the manuscript. The study benefited from funding from the National Science Foundation, under award 59 number EAR-1848192. This research was also supported by the National Research Foundation (NRF) Sin- gapore under its NRF Fellowship scheme (award no. NRF-NRFF2013-06) and by the Earth Observatory of Singapore (EOS) and the Singapore Ministry of Education under the Research Centers of Excellence initiative. The software used in this study is hosted at https://bitbucket.org/sbarbot/unicycle, changeset 891:79134975108c. 60 Chapter3 Thestop-startcontrolofseismicitybyfaultbendsalongtheMain HimalayanThrust This chapter has been published as Sathiakumar, S. and Barbot, S., 2021. The stop-start control of seismicity by fault bends along the Main Himalayan Thrust. Communications Earth & Environment, 2(1), pp.1-11. 3.1 Abstract The Himalayan megathrust accommodates most of the relative convergence between the Indian and Eurasian plates, producing cycles of blind and surface-breaking ruptures. Elucidating the mechanics of down-dip segmentation of the seismogenic zone is key to better determine seismic hazards in the region. However, the geometry of the Himalayan megathrust and its impact on seismicity remains controversial. Here, we develop seismic cycle simulations tuned to the seismo-geodetic data of the 2015 M w 7.8 Gorkha, Nepal earthquake to better constrain the megathrust geometry and its role on the demarcation of partial rup- tures. We show that a ramp in the middle of the seismogenic zone is required to explain the termination of the coseismic rupture and the source mechanism of up-dip aftershocks consistently. Alternative models with a wide décollement can only explain the mainshock. Fault structural complexities likely play an im- portant role in modulating the seismic cycle, in particular, the distribution of rupture sizes. Fault bends are 61 capable of both obstructing rupture propagation as well as behave as a source of seismicity and rupture initiation. 3.2 Introduction At subduction zones and continental collisional margins, the convergence between tectonic plates is ac- commodated by shortening along megathrusts characterized by extensive seismogenic zones that produce some of Earth’s largest earthquakes. Megathrusts oftentimes undergo seismic super-cycles, alternating frequent moderate-size events that break the seismogenic zone in a piece-wise fashion with occasional gi- ant earthquakes that break the entire seismogenic width [Satake, 2015, Philibosian et al., 2017, Philibosian and Meltzner, 2020]. Understanding what determines the spatial extent of seismic ruptures is important to produce realistic bounds on seismic hazards, as larger earthquakes produce greater levels of shaking hazard over wider areas with different frequency contents. In particular, events that rupture to the sur- face, for example, the 2005 M w 7.6 Kashmir earthquake [Avouac et al., 2006], exhibit greater hazards to life and property compared to blind ruptures like the 2015 M w 7.8 Gorkha earthquake that do not reach the surface. However, the underlying mechanics of full and partial ruptures of the seismogenic zone is poorly constrained. The development of seismic super-cycles may be controlled by several factors, such as particular frictional properties [Kaneko et al., 2010, Noda and Lapusta, 2013], geometric complexity [Qiu et al., 2016a], off-fault deformation in the hanging wall [Sathiakumar et al., 2020], and residual stress from past earthquakes [Herrendörfer et al., 2015]. The remaining challenge is to assess the relative contribution of these factors at specific plate boundaries [Shi et al., 2020b, Barbot, 2020]. Most of the relative convergence between the Indian and Eurasian plates is currently accommodated along the Himalayan megathrust, which follows the surface-emergent Main Frontal Thrust (MFT) that soles into a gently dipping décollement called the Main Himalayan Thrust (MHT). The seismic cycle at the MHT embodies both along-strike and down-dip segmentation [Pandey et al., 1995], with surface-breaking 62 ruptures, like the 1934 Nepal-Bihar event [Chen and Molnar, 1977, Sapkota et al., 2013], that unzip the entire seismogenic width, but also blind ruptures, like the 1833 Kathmandu event [Bilham, 1995, Bilham et al., 2004], that break only a small section of the megathrust. The 2015 M w 7.8 Gorkha earthquake represents a partial rupture [Avouac et al., 2015a, Galetzka et al., 2015, Lindsey et al., 2015, Wang and Fialko, 2015, Grandin et al., 2015], but the factors responsible for the rupture termination in the middle of the seismogenic zone are still controversial [Ghoshal et al., 2020a]. Some argue for the predominant role of morphological gradients, [Hubbard et al., 2016, Mugnier et al., 2017, Dal Zilio et al., 2019] while others invoke a frictional control [Michel et al., 2017]. Fault bends along the MHT control many aspects of the Himalayan orogeny at geological time scales. For example, the ramp-décollement structure extending below central Nepal is linked to the Gorkha- Pokhara Anticlinorium [Elliott et al., 2016] and is responsible for the high Himalayan topography [Avouac, 2003]. Along-strike morphological gradients may constitute important boundaries for interseismic strain accumulation in Nepal [Dal Zilio et al., 2020], but the impact of fault bends on down-dip segmentation of the MHT is lesser known, in part because the detailed structural layout of the MHT is controversial. The MHT has been described as a ramp-flat-ramp structure [Duputel et al., 2016, Elliott et al., 2016], a duplex system [Mendoza et al., 2019], or with an additional ramp in the middle of the seismogenic zone [Hubbard et al., 2016]. Since structural complexity may play a significant role in megathrust behavior, it is paramount that we reduce uncertainties in fault geometry. We will show that partial ruptures of the seismogenic zone similar to the 2015 Gorkha earthquake can be obtained either on a wide frictionally unstable décollement or on one with additional fault bends, as long as the characteristic size for earthquake nucleation is much smaller than the seismogenic zone width. However, a sequence with a mainshock followed by shallow dipping aftershocks consistent with observa- tions can only be obtained with a middle ramp. To support these claims, we present numerical simulations of the seismic cycle in Nepal that explain the slip distribution and surface deformation associated with the 63 Figure 3.1: Seismo-tectonic setting and fault geometry of Nepal’s Kathmandu Valley. a Distribu- tion of coseismic slip [Qiu et al., 2016a] in the source region of 2015 M w 7.8 Gorkha, Nepal earthquake (background colors). GNSS measurements are represented by circles for vertical displacement and arrows for horizontal displacements. Near-field GNSS station names are highlighted in yellow and the far-field stations in gray. The black dots represent the aftershock sequence recorded by the Nepal Array Measur- ing Aftershock Seismicity Trailing Earthquake (NAMASTE) network [Mendoza et al., 2019]. The focal mechanisms are for the immediate aftershock sequence [Wang et al., 2017a]. b Sub-surface geological cross-section along profile A-A’ constrained using surface geology [Hubbard et al., 2016]. c Geological cross-section constrained by coseismic geodetic data [Elliott et al., 2016]. The velocity-weakening regions are indicated in red and velocity-strengthening regions in black. The focal mechanisms (colored by focal depth and scaled to magnitude) are superimposed on the cross-section. These cross-sections show the ma- jor thrust faults in the region, including the surface-emergent MFT (Main Frontal Thrust), Main Boundary thrust (MBT), and Main Central Thrust (MCT). 64 2015 Gorkha earthquake using two end-member fault geometry models with and without a middle ramp. The simulations reveal the respective roles of fault bends and stress residuals from antecedent seismicity on down-dip segmentation and the spatial distribution of aftershocks on the Kathmandu section of the Himalayan megathrust. 3.3 Results 3.3.1 SeismiccyclemodelingoftheGorkhaearthquake We investigate the various controls on seismicity using quasi-dynamic numerical models of the seismic cycle based on a physics-based rate- and state-dependent friction law [Barbot, 2019c,b] with the boundary integral method (see Appendix B). We consider two-dimensional models along a representative cross- section perpendicular to the plate boundary that crosses the Gorkha-Pokhara Anticlinorium in the middle of the Gorkha coseismic rupture (Figure 3.1). The distribution of frictional properties and the detailed geometry of the megathrust are poorly known. Therefore, we consider end-member scenarios that allow us to study the impact of fault bends on the seismic cycle. We base the analysis on megathrust geometry models that embody different assumptions about fault complexity. Model E [Elliott et al., 2016] is based on the Gorkha coseismic surface displacements in conjunction with other geological and geophysical data and describes the MHT with a large-scale ramp-flat-ramp geometry. Model H [Hubbard et al., 2016] as- similates regional geological constraints in a balanced cross-section and features a ramp in the middle of the seismogenic zone. Both fault geometry models have a shallow fault bend where the MFT soles into the MHT and a mid-crustal ramp that approximately coincides with the bottom of the seismogenic zone (Supplementary Figure B.1). The spatial distribution of the frictional parameters controls many aspects of fault behavior [Dublanchet et al., 2013, Ray and Viesca, 2017, Ong et al., 2019]. As the spatial variations of frictional properties on the 65 MHT are poorly constrained, we consider simple cases with piece-wise homogeneous frictional properties in the seismogenic zone and the aseismic region underneath. In this context, the style of fault dynamics can be described and predicted predominantly by two non-dimensional parameters [Barbot, 2019b], i.e., the ratio of the seismogenic width with a characteristic nucleation size [Cattania, 2019] and the ratio of frictional parameters that control the static and dynamic stress drops during seismic ruptures [Gabriel et al., 2012b] (equations (B.6) and (B.7) in Appendix B). Here, the characteristic nucleation size results from a dimensional analysis of the governing equation and represents an intrinsic length scale of the mechanical system. It is uniquely defined for the entire rupture sequence up to a constant multiplier. The actual nu- cleation size of simulated events may vary throughout the cycle due to geometric complexity and variable pre-stress at the time of rupture initiation. Numerical simulations with different physical parameters lead- ing to the same non-dimensional parameters exhibit similar dynamics [Barbot, 2019b]. For example, cycles of partial and full ruptures can be obtained on planar fault segments without structural complexity as long as the characteristic size for earthquake nucleation is much smaller than the seismogenic width [Kato, 2003, Lapusta and Rice, 2003, Wu and Chen, 2014, Barbot, 2019c, Cattania, 2019]. Following these the- oretical insights, we explore a space of constitutive parameters that produce a wide range of associated non-dimensional parameters for each of the geometry models. We simulate the seismic cycle for a period of 10,000 years for the models discussed in this section, and for a period of 45,000 years for the two models described in the next section with a background loading rate of 20 mm per year. We discard the first 5,000 years to mitigate the effect of initial conditions. First, we explore the ratio of the seismogenic width with a critical nucleation size (Supplementary Table B.1). Self-emergent partial ruptures can be obtained if this ratio is greater than 20 for model E and greater than 6 for model H (Supplementary Figure B.2). Additional fault bends in the middle of the seis- mogenic zone are responsible for fault segmentation to occur at lower ratios in model H. However, the 66 presence of a middle ramp is not a sufficient condition for partial ruptures, and specific frictional prop- erties in the presence of fault bends must be considered to explain the MHT behavior. Next, we explore the parameters that control the static and dynamic stress drops (Supplementary Table B.2, Supplementary Figure B.3). Cycles associated with a large ratio of seismogenic width to nucleation size and a large ratio of static to dynamic stress drops contain ruptures that start as crack-like and transition to pulse-like to- wards the end (Supplementary Figure B.4). The pulse-like ruptures happen without additional coseismic weakening processes and are specifically a consequence of the high ratio between the static and dynamic stress drops in our quasi-dynamic simulations (see Appendix B). Models with lower ratios do not exhibit this feature, illustrating the strong frictional controls on the details of a seismic rupture. Fault bends in both geometry models attract persistent seismicity under these conditions (Supplementary Figure B.3). For example, giant ruptures sometimes initiate where the MFT soles into the décollement with model E and many ruptures initiate around the middle ramp in model H under some conditions. Up-dip aftershocks are often associated with fault bends, while down-dip aftershocks are mostly associated with the transition zone where the fault behavior changes from velocity-weakening to velocity-strengthening. The transition zone marks the base of the seismogenic zone, but also the top of a large mid-crustal ramp in the Himalayan region, so the variations of friction and geometry are combined at this location. Among the numerical simulations that produce cycles of full and partial ruptures of the MHT, we re- tain those that best explain the geodetic measurements associated with the Gorkha earthquake (Figure 3.1), in particular, the Global Navigation Satellite System (GNSS) and Synthetic Aperture Radar interferometry (InSAR) [Avouac et al., 2015a, Lindsey et al., 2015] data. Even though the GNSS network offers limited resolution for fault imaging [Sathiakumar et al., 2017], four near-field stations show clear coseismic de- formation. The large swath of Advanced Land Observing Satellite (ALOS) line-of-sight displacements constrain coseismic fault slip [Qiu et al., 2016a, Elliott et al., 2016]. With both geometry models, the main 67 features of the coseismic geodetic data and inferred slip distributions can be explained well. Some discrep- ancies between the geodetic data and our forward models could be reduced by adapting the position of the mid-crustal ramp or other adjustments, but the two-dimensional approximation will presumably remain a limiting factor. Even with these constraints, the choices of non-dimensional parameters that can explain the mainshock slip distribution and fault slip within the larger context of super-cycles with partial and full ruptures for both model E and model H is non-unique. Therefore, we narrow our search to simulations that consistently explain both the mainshock and the early aftershocks relocated using advanced broadband waveform modeling techniques [Wang et al., 2017a]. Our exploration of a wide range of models shows that no aftershock up-dip of a coseismic rupture similar to that of the Gorkha sequence can be obtained on a flat décollement using model E. We cannot exclude that shallow aftershocks would occur because of down-dip variations of frictional properties, but we cannot justify the presence of frictional contrasts above the coseismic rupture unless they are accom- panied by morphological gradients. In contrast, shallow aftershocks can be obtained in the presence of a middle ramp, as in model H, for a particular range of frictional parameters that lead to large ratios of static to dynamic stress drops. We therefore focus the remaining discussion on simulation results based on specific frictional parameters (Table 3.1) that are similar for model E and model H, i.e., with a ratio of seismogenic width to characteristic nucleation size of 60 and a ratio of parameters controlling static and dynamic stress drops of 0.9 (Supplementary Figure B.3 d,i). This simulation explains the mainshock and shallow aftershock sequence consistently on model H but only the mainshock with model E. This allows us to investigate the specific differences introduced by morphological gradients within simulations that can explain the cycle of full and partial ruptures of the MHT and the deformation associated with the Gorkha earthquake. 68 Parameter Symbol Value Static friction coefficient µ 0 0.6 Reference velocity V 0 10 −6 m/s Rigidity G 30 GPa Poisson’s ratio ν 0.25 Shear wave speed V s 3× 10 3 m/s Loading rate V L 20 mm/yr Steady-state rate dependence in velocity-strengthening region a− b +1× 10 −2 Steady-state rate dependence in velocity-weakening region a− b −1.98 × 10 −2 Effective normal stress ¯σ 50 MPa Characteristic weakening distance L 5 cm Model H Non-dimensional parameter R u 60.7 Non-dimensional parameter R b 0.9 Down-dip width of the velocity-weakening region W 92 km Model E Non-dimensional parameter R u 59.1 Non-dimensional parameter R b 0.9 Down-dip width of the velocity-weakening region W 89.58 km Table 3.1: Summaryoffrictionalparametersforthepreferredmodels: Frictional parameters and its values using model E [Elliott et al., 2016] and model H [Hubbard et al., 2016] sub-surface fault geometry in the simulations presented in Figures 2-5. The static friction coefficient, reference velocity, shear wave speed, Poisson’s ratio, effective normal stress, frictional parameters a andb, and characteristic weakening distance L are the same for both geometry. The down-dip width of the velocity weakening region is minimally varied, such that the depth of the transition zone is∼15 km for both geometry models. 69 3.3.2 Impactofgeometryonseismicsuper-cycle We start by describing our simulations of seismic super-cycles without a fault bend in the middle of the seismogenic zone. We define a super-cycle as the period that separates two consecutive through-going ruptures and encompasses one or several partial ruptures. Our preferred model using fault geometry E (Table 3.1) produces surface displacements that show a similar signature to the GNSS and InSAR geodetic data of the Gorkha earthquake along profile B-B’ (Figure 3.2). The recurrence times of giant ruptures (e.g., events E17 and E21) that break the entire seismogenic width range between 1,300 and 1,700 years (Figure 3.3a). All the intervening partial ruptures (events E15, E19, and E23) are of similar size and extent, but event E23 matches the Gorkha earthquake data more closely (Figure 3.2). The seismogenic zone is locked during the interseismic period and remains so following the partial ruptures. However, partial coupling and interseismic creep are predominant at the base of the MFT (Sup- plementary Figure B.3), as is often found at fault bends [Sathiakumar et al., 2020, Shi et al., 2020b, Barbot, 2020]. The transition zone between velocity-strengthening and velocity-weakening friction is character- ized by high stresses (Figure 3.4b). Previous ruptures leave behind notable stress residuals that explain subsequent down-dip rupture segmentation due to non-uniform pre-stress along the fault at the time of rupture initiation. All the ruptures initiate at depths between 12 and 15 km because of the presence of the transition zone and the thrust ramp at this depth, which are two separate causes for rapid inter- seismic stress accumulation (Figures 3.3 and 3.4). Model E showcases aftershocks E18 and E22 following the through-going events E17 and E21, respectively. These aftershocks occur near the transition zone, where afterslip immediately below rapidly decays in the months and years that follow through-going rup- tures. However, no aftershocks follow the Gorkha-like event, whether up-dip or down-dip. We now compare these results with our seismic cycle simulation that includes fault bends in the middle of the seismogenic zone. The presence of the middle ramp impacts several interconnected aspects of the seismic cycle, such as rupture size and extent, and recurrence patterns. Model H produces more partial 70 Figure 3.2: Comparison between geodetic data and simulated surface slip for the 2015 Gorkha earthquake.a,b The two tracks of Advanced Land Observing Satellite (ALOS) line-of-sight displacements that recorded the coseismic displacements. c, d Simulated surface displacements using model E [Elliott et al., 2016] and model H [Hubbard et al., 2016] compared with ALOS2-157 and ALOS2-048 InSAR scenes. e,f,g Simulated horizontal (north and east components) and vertical displacements indicated using black solid (model H) and dashed (model E) profiles compared against GNSS horizontal and vertical data. The blue data points correspond to near-field stations indicated by yellow labels in Figure 3.1 and red data points correspond to far-field stations indicated by gray labels. h Simulated fault slip using our numerical method, indicated in black solid (model H) and dashed (model E) lines. Fault slip estimated by inverting geodetic data with geometry H (orange solid) and geometry E (orange dashed) lines. The error bars indicate an uncertainty of 5 percent added to the data. 71 Figure 3.3: Time history of the megathrust behavior, showing earthquake cycles using two end- member geometry models. a Spatio-temporal evolution of slip on the Himalayan megathrust using geometry model E [Elliott et al., 2016]. b Similar simulations using geometry model H [Hubbard et al., 2016]. Earthquakes of different sizes are differentiated using colors; red indicates giant surface-breaking ruptures, purple represents partial ruptures of varying sizes, blue represents Gorkha-like events, and black represents aftershocks. Stars indicate the hypocenter locations. The initial and final slip velocities for seismic events are 0.1 and 0.001 meters per second, respectively. 72 ruptures that break different segments of the MHT with varying sizes. Among these, events H43 and H65 closely match the coseismic signature of the Gorkha earthquake (Figure 3.2). Eventually, the fault completely fails in a single through-going rupture with recurrence times varying between 900 and 2,000 years (Figure 3.3), i.e., almost thrice the variability found in model E (Figure 3.3). Super-cycles are also more aperiodic, containing Gorkha-like ruptures (e.g., events H43 and H65), but also larger events (e.g., H50 and H57), with a wider variability in the number of partial ruptures. High stress concentrations at the fault bends attract persistent seismicity at the middle ramp (Fig- ure 3.4c and d), including ramp-spanning isolated ruptures (e.g., events H55 and H59) and aftershocks (e.g., events H44 and H66). The middle ramp is also capable of terminating partial ruptures (e.g., event H43, H65), thereby acting as a geometric barrier. However, partial ruptures may propagate up-dip through the ramp, creating events of larger size and extent than the Gorkha earthquake (e.g., events H50 and H57). Whether or not ruptures pass the ramp in our simulations is controlled by the spatial distribution of pre- stress, which depends on the history of previous seismicity. The complexity of the rupture sequence can be qualified by the distribution of rupture sizes (Supple- mentary Figure B.5). Chaotic seismic sequences are expected to follow a power-law statistical distribution, an example of which in nature is the Gutenberg-Richter distribution. The super-cycles with model H ex- hibit a power-law distribution of rupture size spanning four orders of magnitude, possibly truncated for small magnitude events due to numerical limitations that impose a lower bound on rupture size. In con- trast, simulations with model E produce events with only three orders of magnitude with a more uniform distribution of rupture sizes. Within the assumptions of the model, fault bends represent zones of rapid stress accumulation in the interseismic period, making these regions prone to earthquake nucleation. As a result, all rupture initiations in our models are in close proximity to a bend. Fault bends therefore play an important role during all stages of the earthquake cycle by attracting seismicity and altering the background stress level. 73 Figure 3.4: Faultslipandshearstressillustratingfaultbendsattractingpersistentseismicity. Slip accumulated on the fault plotted against the down-dip distance using geometry E (a) and geometry H (c). The color indicates the fault slip rate, with red indicating seismic slip and blue represents interseismic creep. The black contours are cumulative coseismic slip values plotted every 5 seconds. The white contours are cumulative aseismic slip values plotted every 50 years. The large slip deficit in the shallow regions of the fault is recovered during through-going ruptures. The stars represent hypocenter locations. (b,d) Shear stress accumulation on the fault; red regions indicate areas of high stress accumulation and blue indicates areas of low stress. Previous ruptures and geometric complexities causes stress concentrations where earthquakes often nucleate. Hypocenter locations are represented by stars, yellow stars represent Gorkha-like events. 74 3.3.3 Aftershocksatfaultbends The hypocenters and focal mechanisms of several relocated aftershocks of the 2015 Gorkha earthquake for the first 12-hour window following the mainshock [Wang et al., 2017a] indicate significant deviations from a planar megathrust geometry (Figure 3.5). Our simulations with model H explain some of the frictional and geometric controls of the up-dip aftershocks of the Gorkha seismic sequence. Events H44 and H66 represent small-magnitude seismic events rupturing the top of the middle ramp in the early postseismic phase of the Gorkha-like events H43 and H65, respectively. For example, the up-dip aftershock H66 happens approximately 14 hours following the mainshock at the shallow fault bend joining the upper décollement and middle ramp (Figure 3.5). Triggering of aftershocks at the boundary of coseismic ruptures can be caused by rapid coseismic stress changes, but this type of events does not occur with model E. As fault bends are associated with rapid interseismic stress accumulation, they constitute prime regions for aftershock triggering because the stress level is already elevated at the start of the mainshock. In fact, in other circumstances, the same regions become the nucleation point of large earthquakes in our simulations. The focal mechanisms of the up-dip aftershocks following the Gorkha earthquake may be explained by the fault bends at this location, confirming the ∼ 30 ◦ dip and the small width of a few kilometers of the middle ramp [Hubbard et al., 2016, Wang et al., 2017a]. The spatio-temporal evolution of the mainshock- aftershock sequence closely matches our simulations with model H. However, the numerical model only explains a single aftershock, presumably due to the computational constraints, e.g., the two-dimensional approximation or the assumption of piece-wise planar fault geometry, that limit the size of events and the geometric complexity that can be resolved numerically. 75 Figure 3.5: AftershocksatfaultbendsfollowingtheGorkhaearthquake. a Geological cross-section along profile A-A’ that describes the MHT geometry using model H [Hubbard et al., 2016]. b The 2015 Gorkha event and the seven early aftershocks during the 12-hour time window relocated using advanced waveform techniques [Wang et al., 2017a]. We monitor the slip velocity at various observation points on the fault surface. The middle ramp consists of three observation points: Points A and C represent the position of the kink formed at the two fault bends of the middle ramp, and B is located at its center. Points D and E are on the décollement that ruptured during the 2015 event. Point E experiences the coseismic slip of the simulated partial event H65, while D represents the termination point that is completely locked. c The evolution of peak slip velocity representing the maximum slip speed attained anywhere in the model, showing the Gorkha mainshock and the subsequent updip aftershock that happened after 14 hours. d, e The evolution of slip velocity is plotted with the colors indicating the position of the observation points. 76 3.4 DiscussionandConclusions The devastating M w 7.8 Gorkha, Nepal earthquake resulted in about 9,000 fatalities and injured million others [Avouac et al., 2015a]. Although the quake took place in a previously known seismic gap [Ader et al., 2012, Stevens and Avouac, 2015], the seismic rupture only relaxed the lower half of the locked MHT region, raising questions about the timing, location, and size of impending earthquakes. As the shallow section of the fault remains locked with no evidence of widespread shallow afterslip [Mencin et al., 2016, Zhao et al., 2017], we contemplate why the rupture stopped in the middle of the seismogenic zone. It is unclear whether the next event will unzip the remaining locked section, repeat the 2015 scenario, or break the entire megathrust in a single giant rupture. Our numerical simulations provide new insights into these fundamental questions. The partial rupture of the seismic gap can be understood as the natural behavior of a chaotic seismic cycle characterized by non-periodic sequences of full and partial ruptures that naturally emerge — with or without a middle ramp — because of the particularly small characteristic nucleation size relative to the seismogenic width, compatible with prior assessments [Michel et al., 2017]. The partial rupture of the lower segment of the MHT followed by full locking of the upper section therefore provides strong con- straints on the frictional parameters of the megathrust. The MHT is capable of seismic ruptures from the deep ramp to the MFT following partial ruptures similar to the 2015 event, and the physical properties that govern the frictional behavior on the fault must combine to a ratio of seismogenic width to characteris- tic nucleation size of at least 60, according to our simulations. However, additional consideration of the Gorkha aftershock sequence indicates that these particular frictional properties are also accompanied by non-planar fault geometry in the middle of the seismogenic zone. Gorkha-like ruptures followed by up-dip aftershocks are systematically followed by giant ruptures with a delay ranging from about 250 to 300 years, based on our preferred model with geometry H. Small ruptures breaking only the upper section of the MHT are therefore unlikely, based on the assumptions of the model. 77 However, more imminent seismic hazard in Kathmandu may originate from seismic ruptures that break different sections of the Himalayan megathrust, initiating outside the Kathmandu Klippe [Qiu et al., 2016a]. Different scenarios of future seismicity are also possible due to the many remaining unknown mechanical properties of the plate boundary, including deformation of the upper plate by faulting and folding [Whipple et al., 2016, Sathiakumar et al., 2020], more complex frictional behavior of the megathrust [Wang and Barbot, 2020], and interactions with other segments of the Himalayan megathrust [Li et al., 2017]. The study illustrates the complex footprint of fault bends and ramps on megathrust seismic cycles. Fault bends are often simply portrayed as geometric barriers that can abruptly terminate rupture prop- agation [King and Nábělek, 1985], but this characterization is incomplete when multiple seismic cycles are considered [Duan and Oglesby, 2007]. In addition, the impact of fault geometry is multi-faceted and cannot be deconvolved from the frictional control (Supplementary Figures B.2 and B.3). Fault bends in the seismogenic zone are always associated with accruing stress during the interseismic period in our simu- lations. This stress can be manifested by creep, for example at the bottom of the MFT, by the propagation of large earthquakes across a ramp, by arrest of ruptures at the middle ramp, by the enhancement of af- tershock activity at the ramp, or by promoting larger ruptures that propagate across the deep décollement and slightly beyond the ramp. In addition, the background stress may be sufficiently elevated across the domain to facilitate through-going ruptures that break the entire seismogenic zone. Fault ramps incite more complexity than isolated velocity-strengthening barriers [Kaneko et al., 2010, Kato, 2016], which only imply a binary cycle of ruptures that either do or do not break the obstacle. The complexity of stress interactions due to the presence of morphological gradients [Romanet et al., 2020] on the megathrust translates to a wider range of rupture sizes and recurrence patterns, manifested by a statistical distribution of earthquake sizes that approaches a power-law and a strongly aperiodic earth- quake recurrence pattern. Power-law statistics of earthquake size distribution may also be obtained on 78 a planar fault under specific physical properties, leading to a ratio of seismogenic zone width to char- acteristic nucleation size greater than about 400 [Cattania, 2019]. Under these conditions, however, the ruptures are systematically starting from the bottom of the seismogenic zone, extending up-dip to varying distances. Our study shows that with non-planar faults, the power-law distribution of earthquake size can be approached for a wider range of frictional parameters than for a planar fault, also producing more widely distributed nucleation points. Variability in hypocenter location may also arise from a damage zone surrounding the plate interface [Thakur et al., 2020], which constitutes another kind of structural complexity. Several models describe the Himalayan megathrust geometry, but the exact geometric features and their role in Himalayan tectonics remain controversial. Our work illustrates the impact of faults bends on seismic cycles and their importance to consistently explain the coseismic and postseismic behavior of the Kathmandu section of the Himalayan megathrust. Our simulations reveal that seismicity at the MHT can be explained by the complex dynamics of frictional sliding with morphological gradients. The frictional properties allow seismic super-cycles of full and partial ruptures, but the rapid inter-seismic stress accumulation at geometric gradients modulates the earthquake size and recurrence patterns. Our results help better constrain the geometry of the MHT and the impact of fault bends on the seismic cycle. Incorporation of high-resolution models of fault geometry [Kyriakopoulos et al., 2015, Handy et al., 2019] into rupture dynamics models may be key to better predict future seismicity at subduction and collision margins. 3.5 Acknowledgments We thank the reviewers for their constructive feedback that helped improve the quality of the manuscript. We thank Xin Wang and Shengji Wei for technical help with seismological data. The study benefited from funding from the National Science Foundation, under award number EAR-1848192. 79 Chapter4 Reconcilinggeologicandgeodeticshorteningrateswithactivefolding intheHimalayanorogenicwedge This chapter has been submitted for peer review as Sathiakumar, S., Barbot, S. and Hubbard, J., 2022. Reconciling geologic and geodetic shortening rates with active folding in the Himalayan orogenic wedge. Science advances, submitted 4.1 Abstract Geologic and geodetic estimates of fault slip motion provide key observational constraints on the me- chanics of crustal deformation, but are often incompatible. This conundrum is well illustrated along the collision zone between India and Eurasia: the geodetic slip rates are∼ 20-30% smaller than geologic slip rates in Central Nepal. Here, we show that these different observations can be reconciled by accounting for active folding and duplexing, key mountain building processes. Our kinematic structural models show that frontal thrusts can accumulate slip at a rate slower or faster than the long-term convergence rate during different phases of structural evolution, depending on the architectural layout of the thrust sheet. Along- strike segments at different phases of structural evolution are likely, and the corresponding geomorphic rates from frontal thrusts at these locations may reflect the internal dynamics of the orogenic wedge. The along-arc geodetic rates, therefore, constitute a better estimate of the deep loading rate between the Indian 80 plate and the Tibetan plateau. The spatial variations in thrust sheet evolution causes spatial and temporal variability of the long-term slip rate, both along-strike and down-dip. This variability can induce earth- quake segmentation patterns that are expected to shift every∼0.3-1.3 Myrs as the hanging wall evolves, with implications for seismic hazards. 4.2 Introduction Over the past 60 million years, the northward movement of the Indian plate finalized the closure of the Tethys Sea, initiating the gradual collision of the Indian and Eurasian plates. This collision has led to the uplift of the vast Tibetan Plateau, with the Himalayan mountains marking the highest global elevations along its southern border [Powell and Conaghan, 1973, Avouac, 2003]. The ongoing collision is accommo- dated across the Himalayan megathrust, one of the largest fault systems in the world. Over geological time scales spanning millions of years, the structure of the Himalayan megathrust has matured with tectonic shortening, and the resulting uplift has produced continuously evolving landscapes [Molnar and Tappon- nier, 1975, Wesnousky et al., 1999, Yin and Harrison, 2000, Lavé and Avouac, 2000, 2001, Hubbard et al., 2016]. At shorter timescales, the system features cycles of earthquakes of varying magnitudes repeating over centuries and millennia [Feldl and Bilham, 2006, Bilham, 2019, Mugnier et al., 2013, Gupta and Gaha- laut, 2014, Wesnousky, 2020]. The relationship between mountain building processes and seismic cycles is critical to understanding this fault system, as well as others around the world. Here, we seek to understand (1) how strain is partitioned between slip on the Himalayan megathrust and internal deformation of the hanging wall, leading to the creation of permanent topography [Dal Zilio et al., 2021, Lamb, 2021], and (2) how the structural complexity of the mechanical system affects short-term seismic cycles in the Nepalese Himalayas. 81 Figure 4.1: Perspective view of the Himalayan subsurface and observational constraints on the shortening-rateinNepal. The three-dimensional fault geometry and the cross-sectional surface across Central Nepal are based on Hubbard et al. [2016]. The fault geometry at two representative cross-sections in Central and East Nepal (CN and EN, respectively) is highlighted, with the frictionally unstable seismo- genic zone in red and the creeping region in black. The contours for the last known earthquakes in these segments, the M w 8.2 1934 Nepal-Bihar earthquake (EN) and the M w 7.8 2015 Gorkha earthquake (CN), are represented in yellow. Central Nepal features a middle ramp that is absent in East Nepal. Fault system evolution allows break-forward episodes of one or more ramps, allowing accretion of Indian lithosphere onto the upper plate. The local long-term slip on individual fault segments, which may vary down-dip due to internal deformation of the orogenic wedge, is represented by double arrows. The interlinks between mountain building, plate tectonics, and seismic cycles can be difficult to study due to the disparate time and length scales of deformation processes. However, the apparent contrast between geodetic and geologic slip rates in the Himalayan mountain belt allows us to probe the relationship between short-term and long-term processes (Figure 4.1). Geologic fault slip rates at the frontal section of the orogenic wedge can be estimated from abandoned fluvial terraces using geometric relationships between faulting and folding. In Central Nepal, such geological observations provide a local shortening 82 rate of21± 1.5 mm/year across the Main Frontal Thrust (MFT) [Lavé and Avouac, 2000, 2001]. In contrast, recent geodetic data suggest an average total shortening rate of just15.5± 2 mm/year [Lindsey et al., 2018, Li et al., 2018]. The shortening rate in the shallow part of East Nepal is highly uncertain, and it is unclear if geodetic and geologic rates agree or disagree. Additionally, there may be longitudinal variations in strain accumulation and partitioning [Bollinger et al., 2014, Hamahashi et al., 2022]. Resolving this conundrum is key to better understanding seismo-tectonic deformation at this collision zone. In this study, we show that the geologic and geodetic slip rates are fundamentally compatible, rep- resenting possible proxy data for folding processes taking place in the collision zone. We can reconcile these observations by accounting for the internal deformation of the orogenic wedge by active folding in the hanging wall and underplating in the footwall of the Main Himalayan Thrust (MHT). We study the structural evolution of the megathrust along two representative cross-sections in Central Nepal (CN) and East Nepal (EN) that feature distinct degrees of structural complexity. We show that the rate of strain accommodation in the hanging wall changes at timescales of millions of years during stages of uplift and duplexing. This in turn causes dynamic variability of the long-term slip rate, both down-dip and along the megathrust. Using numerical simulations of seismic cycles with active folding, we show how the in- ternal deformation of the hanging wall may lead to different spatio-temporal patterns of slip budgets and earthquake statistics. 4.3 Results 4.3.1 ArchitectureandstructuraldevelopmentoftheHimalayanorogenicwedge The Himalayan orogenic wedge consists of a forward-propagating sequence of thrust faults; from north to south, these are the Main Central Thrust (MCT), the Main Boundary Thrust (MBT), and the Main Frontal Thrust (MFT, the most active splay in the system), respectively. These splay faults sole at depth into the 83 Main Himalayan Thrust fault (MHT), which accommodates crustal shortening and related deformation (Figure 4.1). The base of the MHT consists of a mid-crustal ramp, a steeply dipping thrust fault segment that connects two levels of décollement (gently dipping thrust faults). Slip along this ramp influences the long-term uplift of Himalayan rocks, as well as short-term seismicity due to persistent strain accumula- tion [Avouac, 2003, Hubbard et al., 2016, Dal Zilio et al., 2021]. In Central Nepal, magnetotelluric data illuminate the mid-crustal ramp as a conductive zone at∼15 kilometers depth [Lemonnier et al., 1999]. The position of the mid-crustal ramp in East Nepal can be traced using the surface expression of Gorkha-Pokhara anticlinorium (GPA) as a proxy. The structural heterogeneity in the fault architecture is illustrated by the curvature of GPA along strike. Another struc- tural complexity in the MHT is the presence of a middle ramp in Central Nepal within the seismogenic zone [Hubbard et al., 2016]. The middle décollement ruptured during the recent 2015M w 7.8 Gorkha earth- quake and the 1833 rupture [Avouac et al., 2015b, Bilham, 1995]. Unlike the areas to the east and west, Central Nepal has not experienced a megathrust-spanning full rupture (from the base of the megathrust to the surface) in the last 1,000 years or more, indicating that strain remains to be released and suggesting an impending seismic hazard. This middle décollement, however, vanishes to the east, as the ramps below and above it merge into a single large mid-crustal ramp [Hubbard et al., 2016]. The last known earthquake that ruptured this section of Nepal is the 1934 Nepal-Bihar earthquake [Pandey and Molnar, 1988]. To explore the evolution of the fault system and its impact on the earthquake cycle, we use repre- sentative two-dimensional cross-sections along Central and East Nepal with variations in down-dip fault geometry. We apply quantitative models of fault-bend folding to study the geometric and kinematic de- velopment of folds, and create a sequence of cross-sections that illustrate the structural evolution of the hanging wall as shortening progresses, paying particular attention to the slip rates on different fault seg- ments. Following standard fault-bend folding methods, we hold the thickness of stratigraphy across each 84 fold and the cross-sectional area constant [Suppe, 1983]. We model axial surfaces or fold axes as infinites- imally narrow, approximating the small thickness of fault and stratigraphic bends relative to the 200-km length scale of the plate boundary [Sathiakumar et al., 2020]. Large anticlines, like the GPA, are created due to slip across multiple fault bends with a flat-ramp-flat geometry. When stratigraphic layers traverse fault bends, they are folded and uplifted, building topography. Notably, when these layers are at an an- gle to the fault, the slip magnitude on either side of the fault bend will be different, due to the geometric conservation of total cross-sectional area. In our progressive cross-sections, we track these changes in slip rate across the multiple fault bends in the MHT as shortening progresses (Figure 4.2). The ratio of slip below and above the bend depends on the shape of the fault and the associated fold [Suppe, 1983, Shaw et al., 2005]. 85 ? crestal uplift uplift velocity consumed Transition from uplift to broadening uplift crestal broadening velocity created A 1 A’ B’ B 1 C 1 C’ A 1 B 2 A” B’ C 1 A 1 B 2 A” C 2 C” Quarternary alluvium Siwaliks Lesser Himalayas Greater Himalayas Tibetan - Tethys Tertiary Plutonic Surface geology MFT MBT 5 km 15 km 25 km C C’ E E’ Topography lines Active axial surface Inactive axial surface Anticlinal fault-bend fold Synclinal fault-bend fold v 1 > v 2 v 1 = v 2 γ γ φ v 1 v 2 uniform slip rate v 1 < v 2 v 1 v 2 θ φ velocity created γ γ v 2 v 1 φ v 1 = v 2 uniform slip rate v 1 v 2 θ γ γ φ velocity consumed faster slip rate slower slip rate stage EN 1 stage EN 2 stage EN 3 E E’ 0 100 km A 1 A’ B’ C 1 C’ B 1 crestal uplift crestal uplift Transition from uplift to broadening crestal broadening uplift crestal uplift uplift Transition from uplift to broadening crestal broadening break-forward thrust D 1 D’ E 1 E’ A 1 A’ B’ B 1 C 1 C” D’ E 1 E’ D 2 A 1 A” B’ B 2 C 1 C” D’ E 1 D 2 A 1 C” D’ E 1 D 2 B 2 A” C 2 C 2 ’ A 1 C” E 2 D 2 B 2 A” C 2 C 2 ’ uplift velocity created velocity consumed velocity consumed crestal broadening velocity created stage CN 1 stage CN 2 stage CN 3 C C’ Cross-sections 40 km Fault slip rate lower higher 100 75 125 (% deep loading rate) a c d e f g b h i j k stage CN 4 stage CN 5 stage CN 6?? Figure 4.2: StructuralevolutionoftheMHTinCentralandEastNepal. a) The locations of the two- dimensional profiles in Central (CC’) and East (EE’) Nepal are shown on a geological map of Nepal. b) Illustration of fault-bend fold theory concepts applied to reconstruct the structural development of the Himalayan megathrust. The long-term slip rate of fault segments varies depending on the fold shape and the history of deformation. The Main Frontal Thrust (MFT) can slip at a rate faster (pink) or slower (yellow) than the long-term loading rate due to active folding and shortening. Reconstruction of hanging wall evolution along profile EE’ in East Nepal (c-e) and profile CC’ in Central Nepal (f-k) Our two-dimensional cross-sections of Central and East Nepal (CN and EN, respectively) are based on the three-dimensional fault geometry described by [Hubbard et al., 2016]. In the early developmental stages (stages CN1, CN2, EN1), the stratigraphic layers are parallel to the underlying fault everywhere except above the ramps, which cut through existing layers (Figure 4.2c,f, g). As the hanging walls of the ramps are pushed up, they travel over fault bends onto the shallower flat fault segments, causing folding and a change in slip rate. These types of fault bend, called anticlinal fault bends, are associated with a decrease in slip rate; in our cross-sections, this implies lower long-term fault slip rates along up-dip fault segments 86 compared to the deep loading rate (Figure 4.2). With typical Nepal fault structure parameters, we obtain a decrease of∼12-15% (stages EN1, CN2) or ∼25% (stage CN1) in frontal fault slip rate, reflecting one or two anticlinal fault bends, respectively. During these stages, the structures are in a predominantly “crestal uplift” phase, as the crests of the anticlines continue to rise, contributing to the growth of topography (Figure 4.2). As shortening continues, the kink bands widen, and the inactive axial surface associated with the base of the ramp reaches the top of the ramp. At this point, the hanging wall rocks above the ramp become parallel to the ramp, and the slip rate across the upper bend becomes constant (stages CN3, EN2). This brings the local long-term slip rate on all the fault segments back to the deep loading rate (Figure 4.2d,h). The structural signature of this process now transitions to “crestal broadening”, with passive transport of the folded kink bands along the upper detachments, and the crests of the associated anticlines grow in width but not height. In the later stages (stages CN4, CN5, EN3), inclined stratigraphic layers which already went through one folding process start to climb the up-dip ramps. The resulting folds, called “synclinal fault-bend folds”, create an increase in slip rate on the up-dip part of the fault (Figure 4.2e,i, j). In our cross-sections, this creates a∼13% (stages EN3, CN4) or ∼26% (stage CN5) increase in long-term slip rate at the frontal thrust compared to the deep loading rate, reflecting one or two synclinal fault bends, respectively. Thus, the frontal part of the megathrust can slip slower or faster than the deep loading rate, averaged across the seismic cycle. These slip rates can last for hundreds of thousands to a few million years before transitioning to a new, stable slip rate. The process described here reflects the evolution of the hanging wall over a static fault geometry. The next stage of the fault system involves breaking forward of one or several of the ramps, allowing accretion of Indian lithosphere onto the upper plate. Break forward of the mid-crustal ramp will start a new cycle of fault-bend folding (e.g., stage CN6, Figure 4.2k). 87 In the Himalayas, these long-term folding mechanisms must impose heterogeneous long-term loading rates along different segments of the plate interface that can differ measurably from the deep loading rate. As the local loading rate affects the recurrence time and other characteristics of earthquakes [Barbot, 2019b, Sathiakumar et al., 2020], the structural development of the orogenic wedge will modulate the short-term seismo-tectonic behavior. In the next section, we explore numerical models of the seismic cycle that incorporate the boundary conditions and internal deformation representative of the different stages of growth of the Himalayan orogenic wedge. We use our results to discuss the impact of off-fault deformation on seismic hazards along the Central and East Nepal representative cross-sections. 4.3.2 Earthquakecycleswithspatiallyvariablesliprateandactivefolding To study the impact of folding in the orogenic wedge on the seismic cycle, we employ physics-based nu- merical simulations of earthquake sequences using rate- and state-dependent friction coupled to internal deformation of the hanging wall across active axial surfaces [Sathiakumar et al., 2020]. In order to iso- late the effect of structure on the two-dimensional Central and Eastern cross-sections (see Supplementary Tables C.1 and C.2), we use the same effective physical properties for both models, based on earthquake sequence simulations of the 2015M w 7.8 Gorkha earthquake [Sathiakumar and Barbot, 2021]. The tran- sition from deep steady-state velocity-strengthening behavior to shallower unstable weakening is fixed at∼15 km. The deep loading rate in both models is set to 17 mm per year, which is an intermediate rate between the best-fit geodetic and geologic rates, and within the uncertainty range of geodetic models. We incorporate axial surface dynamics to account for off-fault deformation associated with fault bends. To bridge long-term (> 100,000 to a few millions of years) and shorter-term or seismic cycle timescales (seconds to thousands of years), we implement our kinematic forward model of the structural evolution of the hanging wall: we generate one model for each stage of development in Central and Eastern Nepal (stages CN 1-5, EN 1-3, Figure 4.2, Supplementary Figure C.1). For each of these models, we simulate 88 seismic cycles for 20,000 years and compare the results to study the evolution of seismic behavior as the hanging wall deforms. For each region, the underlying fault geometry and frictional properties remain the same, while the slip rates of the fault segments and the orientations of the axial surfaces within the hanging wall stratigraphy vary (Supplementary Figure C.1). This model allows us to study how long-term geological evolution impacts short-term seismic processes. Great rupture unstable weakening velocity strengthening Major rupture 1 1 1 1 1 1 1 1 2 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 5 5 1 1 1 1 1 2 2 2 3 3 2 1 1 1 2 2 2 3 3 3 1 1 1 1 2 3 Strong rupture 1 1 1 3 1 1 1 1 2 2 3 3 4 4 5 1 1 2 2 4 4 4 4 5 3 3 Central Nepal Ru ~ 60 Rb = 0.286 unstable weakening velocity strengthening 3 3 3 3 3 3 1 2 1 3 3 3 1 1 3 3 3 2 2 1 1 East Nepal EN1 EN2 EN3 CN1 CN2 CN3 CN4 CN5 20 100 60 40 80 Down-dip rupture width (km) larger SZ with middle ramp smaller SZ without middle ramp East Nepal Central Nepal MFT slip rates 500 1000 1500 2000 2500 3000 Recurrence time (years) EN1 CN1 CN2 CN3 CN4 CN5 EN2 EN3 Central Nepal East Nepal a b c d e f Ru ~ 60 Rb = 0.286 3 3 Great rupture Major rupture Strong rupture 0 10 20 30 Frequency (integer) 1 1 2 75 88 100 113 126 85 100 113 % deep loading rate (17 mm/year) Missing ruptures sizes strong: moderate sized partial ruptures great ruptures major ruptures Missing ruptures sizes great ruptures major ruptures strong: moderate sized partial ruptures Figure 4.3: Seismotectonic implications of down-dip long-term slip-rate variations. a and b) Il- lustrations of two-dimensional fault geometry models for Central and East Nepal with fault segments and sample active axial surfaces. Unstable weakening is denoted in red, velocity strengthening in black.R u and R b are non-dimensional numbers that control the dynamics of fault behavior in the rate-and-state friction framework.R u is the ratio of seismogenic zone width to characteristic nucleation size, andR b is the ratio of static to dynamic stress drop. c and d) Down-dip rupture width illustration for megathrust-spanning full ruptures (great ruptures), multiple-segment large ruptures (major ruptures), and single-segment rup- tures (strong ruptures); stars indicate hypocenter locations. e) Colored histogram indicating frequency of ruptures of a given width for each of the stages in our simulations (CN1-CN5, EN1-EN3). Darker colors indicate more frequent ruptures. The lower bar indicates the slip rate of the MFT for each stage; colors of the lower bar are taken from Figure 4.2f. Recurrence times of great ruptures for all models (black: Central Nepal stages CN1-CN5; blue: East Nepal stages EN1-EN3). 89 We simulate super-cycles of earthquakes that express megathrust-spanning full ruptures (great rup- tures), multiple-segment large ruptures that may or may not break the surface (major ruptures), and smaller ruptures that break individual fault segments (strong ruptures). Here, we observe that each seismic super- cycle starts and ends with a through-going “great” rupture. Our models show that variations in slip rate impact several earthquake cycle characteristics, including down-dip rupture width (and therefore magnitude), hypocenter location, and the recurrence times of great and partial ruptures (Supplementary Figures C.2 and C.3). First, frontal fault slip rates impact the recur- rence times of great ruptures (Figure 4.3). When the frontal slip rates are higher than the deep loading rate, our models produce more frequent great events with less slip per event, compared to the more infrequent, larger events associated with slower frontal slip (mean recurrence time of∼1,000 years vs. ∼1,800 years, respectively). The recurrence time of great ruptures is also more periodic for faster frontal slip compared to more chaotic for slower frontal slip. In addition to impacting the largest ruptures, long-term slip-rate variability also modifies the number of partial ruptures between two great ruptures (i.e., within a complete super-cycle). Models in which the frontal slip rate is lower or the same as the loading rate produce twice as many incomplete ruptures per super-cycle compared to those with faster frontal slip. These partial ruptures are of variable size, ranging from multi-segment major ruptures to single-segment strong ruptures. Models with faster frontal slip also do not encompass the full spectrum of rupture sizes observed for lower frontal slip rates and exhibit less scattering of hypocenter locations, which cluster around fault bends in all cases (Figure 4.3 and Supplementary Figure C.2). The recurrence time and behavior of single-segment strong ruptures bounded by ramps, a behavior seen in the 2015 Gorkha rupture in Central Nepal, are uniform throughout the earthquake cycle in all stages 90 of structural evolution, and nucleate at the transition between the velocity-strengthening and velocity- weakening parts of the fault (Supplementary Figure C.2). These partial ruptures of the MHT can precede or succeed great ones, with inter-event times between 100 and 500 years. The structural differences between Central and East Nepal manifest as distinct characteristics of the earthquake cycle. The seismogenic width of East Nepal is smaller, reducing the size of down-dip ruptures. It also does not contain a middle ramp; this reduces the opportunity for geometric or slip rate segmentation, and we observe a more limited range of rupture sizes as a result. Nevertheless, both the fault structure and the down-dip slip-rate segmentation determine the seismo-tectonic attributes of the system (Figure 4.3 and Supplementary Figure C.3). 4.4 Discussion Understanding the deformation of the Himalayan orogenic wedge requires integrating datasets that cap- ture widely variable timescales. Along the MFT, rock uplift data, trenching, river terraces, and subsurface geology are some of the geomorphological and paleoseismological data that have provided important con- straints on earthquake characteristics and seismic hazard in these thrust-faulted terrains [Burbank et al., 1996, Lavé and Avouac, 2001, Godard et al., 2014, Bollinger et al., 2014, Mishra et al., 2016, Rajendran and Ra- jendran, 2021, Ghoshal et al., 2020b]. A survey of the Holocene fluvial terraces perched above the Bagmati river in Central Nepal across a 20 km wide deformation zone suggests a fault slip rate of21± 1.5mm/year in the Himalayan foreland [Lavé and Avouac, 2000, 2001]. At this longitude, the most recent estimates of convergence from geodesy are∼ 20-30% smaller (15± 2.4mm/year) [Lindsey et al., 2018, Li et al., 2018]. As we show, fault-bend folding makes it possible for the frontal rate to exceed the deep loading rate, but only under specific conditions. In order for frontal fault slip rates to be higher than the deep loading rate in Central Nepal according to our model, the structural evolution must be in a stage similar to stage 91 CN4, CN5, or CN6. In these, stages, inclined beds are being refolded at one (CN4 and CN6) or two (CN5) additional updip ramps (Figures 4.2 and 4.4). A deep loading rate of 13-17 mm/year can lead to a long-term slip rate at the frontal thrust of 14.5-19.5 mm/year during stage CN4 (∼ 13% increase) and 16.5-21.5 mm/year during stage CN5 (∼ 27% increase). Structural cross-sections based on surface geology at this location [Hubbard et al., 2016] suggest that the mid-crustal ramp has moved from crestal uplift to crestal broadening, and that inclined beds are climbing the up-dip middle ramp. This increases slip updip similar to processes described in both stages CN4 and CN5 in our models (Figure 4.4a). Our simplified models capture the first-order approximations of fault evolution with a fixed geometry, and more accurate models can capture other complexities like break-forward processes of the frontal fault system. For instance, the middle ramp originally formed as the lower part of a surface-breaching Main Boundary Thrust fault (MBT) and eventually its hanging wall rocks were uplifted and eroded away. In our model, these hanging wall rocks are instead transported up the fault, and in stage CN5 cause an additional increase in slip rate updip as they reach the MFT. The slip rate increase in stage CN5 is therefore not supported at this location, although intermediate ramps like this one may be active in other locations. Thus, although the difference in slip rate from geological and geodetic estimates in Central Nepal is best fit by stage CN5, stage CN4 is more representative. 92 1.0 2.0 3.0 Time progression during structural evolution (million years) 0 4.0 5.0 Slip rate (mm/year) ~11,000 years East Nepal Central Nepal a stage EN2 stage EN1 stage EN3 stage CN2 stage CN6 Toward present time transition to crestal broadening MCT MCT MBT MFT Kathmandu Klippe velocity created stage CN4 d 22 20 18 16 14 12 24 10 stage CN1 stage CN3 stage CN4 stage CN5 ???? start of stage CN/EN Holocene tectonic geomorphology GNSS velocities Central Nepal East Nepal c 30 km geologic estimate geodetic estimate shortening rate (model) b CN2 CN1 CN3 CN4 CN5 EN2 EN1 EN3 Inclined bed -25% of long-term rate +26% of long-term rate -15% of long-term rate +13% of long-term rate Figure 4.4: Variationsoflong-termslip-ratesontheMFTovermulti-millionyeartimescalesdue toactivefoldingandduplexingcomparedtothegeodeticloadingrate. a and b) Fault-bend folding predicts different frontal fault slip rates depending on the history of kinematic development (black and blue horizontal lines representing Central and East Nepal, respectively). In Central Nepal (a) the geologic shortening rate at the MFT (estimated from fluvial terraces, indicated by the shaded brown region) is 20- 30% higher than the geodetic loading rate for the MHT (shaded yellow region). An increase in frontal long-term slip rate is possible due to the uplift of inclined rock layers in Central Nepal in stages CN4-CN6. In East Nepal, the geologic rates are highly uncertain, and so it is not currently possible to distinguish between different stages of evolution based on rates alone. c) Evolution of long-term slip rate on the shallow MHT and MFT in Central Nepal (black) and East Nepal (blue). The modified slip rates due to structural evolution are sustained for many seismic super-cycles, spanning hundreds of thousands to a few million years. d) Geological cross-section of Central Nepal Hubbard et al. [2016], showing inclined beds climbing the middle-ramp, causing a long-term slip rate increase in the up-dip segments, similar to Stage CN4. 93 In East Nepal, slip rates remain uncertain. Data from sediments in boreholes in the Himalayan foothills reveal the deformation history and MFT fault slip rates on the Bardibas fault, the southernmost segment of the frontal thrust [Hamahashi et al., 2022]. At the Bhabsi river, shortening rates are∼ 3-12 mm/year, accommodating 20-80% of the geodetic shortening rate. [Hamahashi et al., 2022] suggest that this lower slip rate is due to partitioning of shortening onto other hanging wall faults, including the Patu Thrust in the north, as this is a step-over region between the two faults; incorporating this additional shortening could plausibly make the total shortening comparable to geodetic rates. Notably, river terraces from the Ratu river, just east of the Bhabsi River suggest uplift rates of 10-12 mm/year on the Bardibas Thrust [Bollinger et al., 2014], possibly suggesting 20-35 mm/year of shortening assuming a 20-30 degrees dipping fault. These widely variable rates, from far below to far above the geodetic loading rate, highlight the difficulty of measuring geologic rates. In addition to structural uncertainties, non-tectonic climate signals can affect the uplift of river terraces [Hamahashi et al., 2022]; geologic studies have therefore often relied on the overall shortening budget, inferred from geodesy, to estimate the loading rate. Here, we demonstrate that this budget can also vary depending on the hanging wall structure. In the case of Eastern Nepal, the total MFT shortening budget could be consistent with the lower slip rate estimate from the Bhabsi River (in stage EN1) or the higher slip rate estimate from the Ratu River (in stage EN3). Given the uncertainties in estimates and the likelihood of additional active faults and sedimentary processes, all stages remain possible; further work on the large-scale structure of the orogenic wedge may be more effective than assessment of the geological rates in this case. As mountain-building processes continue in the Himalayas, ramps break forward to create duplexes, and the system transitions to a new cycle of structural evolution. A hypothetical stage CN6 shows how southward migration of the mid-crustal ramps could shift the fault structure in Central Nepal to one more like East Nepal. In all scenarios, geodetic data can estimate the long-term deep loading rate representative 94 of the India-Eurasia convergence across the MHT, while the shallow geomorphic rates can help constrain the dynamics of the hanging wall in the toe of the orogenic wedge. 4.5 Conclusions The discrepancy between geodetic and geologic estimates is illustrated by temporal variations in long- term frontal fault slip rates. In fold-and-thrust belts, we may be able to reconcile some of these disparate observations by accounting for the internal deformation of the orogenic wedge, which accommodates a variable strain rate over millions of years, leading to frontal fault slip rates ranging from 70% to 130% of the deep convergence rate. The variable rates can be stable for 100,000 to millions of years as the orogenic wedge evolves. Measurements of shortening from geodesy, which describe interseismic deformation over large spatial regions, are representative of the long-term convergence rate between tectonic plates. This deformation is accommodated elastically around the deep regions of faults, and drives intermittent slip up-dip. In contrast, measurements with sensitivity to shallow deformation, such as geomorphic signals, incorporate the kinematics of the orogenic wedge. Here, we show that these kinematics can lead to local shortening rates that are lower, comparable, or higher, than the loading rate depending on the history of internal deformation and the morphology of the underlying megathrust. These key pieces of information can be combined to build a more comprehensive picture of the dynamics of a plate boundary, from the foreland to the hinterland. In Nepal, the deep convergence rate may also change with the life cycle of faults in the Tibetan Plateau and other changes in distant boundary conditions [Tapponnier et al., 2001]. These effects will be modulated by the internal deformation of the orogenic wedge, which impacts the local shortening most effectively near the frontal section. Our numerical models of earthquake sequences show that spatial variations in slip rate have a major impact on earthquake recurrence patterns. For the fault system in Central Nepal, periods of time when the 95 frontal slip rate is lower than the loading rate are likely to feature more partial ruptures of the megathrust, leaving behind concentrated pockets of high stresses that affect the timing and slip of megathrust-spanning ruptures. In contrast, periods when the frontal slip rate is higher than the deep loading rate are likely to feature more frequent great ruptures that extend through the entire seismogenic zone. An important step for seismic hazard assessment would be to incorporate the structural layout of the orogenic wedge, in order to better understand the expected patterns of seismicity. Validating these models will require integrating a wide range of geologic and geophysical datasets with complementary spatial and temporal resolution. Models spanning the relevant time scales of seismic cycles and long-term geodynamics will be key to better understanding the coupling between earthquakes and mountain building at convergent boundaries. 4.6 Acknowledgments This study is supported in part by the National Science Foundation under award number EAR-1848192. This study is also supported by the Earth Observatory of Singapore via its funding from the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centers of Excellence initiative. The software code is publicly available at https://zenodo.org/record/4679513. 96 Chapter5 ConcludingRemarks The revolutionary theory of plate tectonics gained momentum merely sixty years ago because the notion that rigid rock blocks float on molten rocks, allowing continents to slowly drift, was unimaginable. How- ever, observations of matching rock formations and fossils from far-apart continents, and magnetic data from the ocean, led to the birth of the geological sciences. Today, our technological progress allows us to measure centimeter and millimeter-scale deformation patterns on the Earth’s surface. We have developed several state-of-the-art techniques to measure rates of crustal deformation, determine the age of landscape features, and monitor earthquakes, landslides, tsunamis, sea-level rise, and volcanic activity. Despite these giant strides and scientific revelations, we still lack an in-depth understanding of several complex Earth system interactions. One of the greatest challenges we face today is the lack of datasets and models that span disparate time and length scales. We need such datasets to understand the complex mechanisms that create evolving landscapes as a result of competing subterranean and surface processes. One way to make progress despite inadequate data is to integrate our knowledge, understanding, and datasets from multidisciplinary scientific fields. This thesis undertakes collaborative research in this spirit to attempt to elevate our current understanding of faulting and folding and its impact on seismotectonics in the short-term and mountain-building in the long-term. In this chapter, I briefly summarize my goals, research questions, and the results of my investigations. 97 Fault-bend folds were first observed nearly 100 years ago, and the kinematics of these map-scale folds due to bends in the underlying fault were described more than 40 years ago. However, the vast majority of dynamic rupture models typically assume simplified, smooth fault planes and do not consider the effect of off-fault deformation due to these ubiquitous geologic structures. In Chapter 2, I use a simple case of a single fault-bend fold to investigate how structural complexity impacts the earthquake cycle. In particular, I investigate the role of fault bends on fault dynamics and earthquake segmentation. Since the parameter space is flooded with an astounding range of frictional and geometric parameters, this study establishes key factors that control earthquake cycle properties: fault friction, the magnitude of fault bending, and the geometry of hanging wall bedding planes. Seismic cycle complexity increases with increase in the ratio between seismogenic width and nucleation size, a friction parameter that controls the minimum size of the earthquake. Fault bend magnitude also plays an important role: as the magnitude of fault bending increases, the complexity of the earthquake cycle increases. Another interesting takeaway from this study is the impact of sediment cut-off angle or the hanging wall geometry on fault behavior. We show that when the incoming sediments are not parallel to the underlying fault, a case that often happens when inclined sedimentary packages are brought to an upper décollement level during shortening, folding can cause changes in slip rate across the fault bend. This process has a profound impact on both the long- term and short-term processes. In this case, earthquakes that initiate, propagate and terminate near the surface become more recurrent. If these shallow earthquakes occur in the accretionary zones of subduction margins, they are capable of large seafloor uplifts that may lead to devastating tsunamis. Therefore, the interplay between the shape of the fault and fold, and the fault friction properties plays a significant role in earthquake cycle dynamics. This study constitutes an important first step towards linking the long-term deformation predicted by fault-bend fold theory with the short-term deformation that occurs during the seismic cycle. 98 Next, I study real-world fault systems like the Himalayan megathrust that is overrun with fault bends. In Chapter 3, I investigate the details of the Main Himalayan Thrust fault (MHT) geometry using physics- based numerical simulations and by comparing them against the geodetic data of the 2015M w 7.8 Gorkha earthquake. This earthquake created an opportunity to probe important and relevant questions about the megathrust behavior and helps us answer the following questions: 1. What is the role of fault bends in the Nepal seismic cycle? 2. What controls earthquake segmentation here? Why did the Gorkha earthquake stop in the middle of the seismogenic zone, even though geodetic evidence points to a fully locked fault zone? 3. What are the future seismic hazards in the region? Our study, which takes into account both frictional and geometric complexities, reveals that the im- pact of fault geometry is multi-faceted and cannot be deconvolved from the frictional control. Fault bends impose a complex footprint on the megathrust seismic cycle. A middle ramp in the seismogenic zone is essential to consistently explain both the coseismic rupture extent and the early up-dip aftershocks. The complexity of stress interactions due to morphological gradients translates to a wider range of rupture sizes and extent. We find that models with or without a middle ramp in the seismogenic zone accumulate large slip deficits in the shallow regions of the fault, but are recovered during the complete ruptures only. Isolated up-dip ruptures, therefore, seem unlikely based on the assumptions of this model. This could be mainly because when stress is brought to the front, the entire fault is ready to break all in one go, so the chances of a partial rupture of just the up-dip locked portion seem unlikely. Several missing pieces of the Himalayan megathrust puzzle still remain, as we do not fully understand the impact of off-fault deforma- tion and folding and how these long-term mountain-building processes affect the short-term seismic cycle processes. A future goal would be to study the interactions with other fault segments of the Himalayan megathrust and how the along-strike geometric and frictional gradients affect the seismic cycle. 99 In Chapter 5, I focus on the long-term structural fault evolution of the MHT to better understand mountain-building processes. Mountains are a focal point of evolutionary processes for both lifeforms and landforms, and rock records of uplift help us understand the kinematic and structural aspects of deforma- tion over long timescales. In central Nepal, geological data of abandoned fluvial terraces from the frontal sections of the megathrust provide a time-averaged estimate of fault motion, which are 20-30% higher than current geodetic rates. It is vital to understand why these rates are different or appear different because it could allow us to probe into additional physical mechanisms contributing to the overall deformation in the Himalayas. In this study, I show that the apparent discrepancy can be indicative of a deforming crust that can cause non-uniform slip rates on the basal fault (MHT) and the frontal ramps (Main Frontal Thrust or MFT). This is because the Himalayan wedge is not an elastic body, but exhibits permanent deforma- tion due to bulk folding in the hanging wall. However, different sections of the megathrust can exhibit different structural development, with other processes like duplexing and underplating also contributing to the deformation process. These temporal variations in fault slip rates can profoundly impact seismic cycle processes, altering fault behavior and seismic cycle characteristics like recurrence times and rupture extent of earthquakes. Advances in seismic hazard assessment in Nepal and other convergent margins around the world will require integrating a wide range of geologic and geophysical datasets and models with complementary spatial and temporal resolution. These improved models stand to gain better insight by incorporating the structural layout of the orogenic wedge. This allows us to better understand the cou- pling between earthquakes and mountain building. Our study is one of the first steps towards building these models that span the relevant time scales of seismic cycles and long-term geodynamics. I conclude my thesis by noting that these new and improved numerical models based on geophysical and geological datasets and kinematic models provide a multipronged approach to understanding crustal dynamics. To build better frameworks for seismic hazard analysis and to improve our understanding of Earth’s geological systems, we have to better understand the physical parameters that control the dynamics 100 of these systems. However, the arena for physical parameters is large, and observational datasets are hardly sufficient. In the future, an upsurge in our understanding and further strides in scientific discovery are possible through an interdisciplinary and iterative approach, that combines datasets and models from varied scientific fields of inquiry that study different interconnected aspects of this complex system. 101 Bibliography JR Elliott, R Jolivet, PJ González, J-P Avouac, J Hollingsworth, MP Searle, and VL Stevens. Himalayan megathrust geometry and relation to topography revealed by the gorkha earthquake. NatureGeoscience, 9(2):174, 2016. Judith Hubbard, Rafael Almeida, Anna Foster, Soma Nath Sapkota, Paula Bürgi, and Paul Tapponnier. Structural segmentation controlled the 2015 mw 7.8 gorkha earthquake rupture in nepal. Geology, 44 (8):639–642, 2016. Sharadha Sathiakumar, Sylvain Barbot, and Judith Hubbard. Earthquake cycles in fault-bend folds. Journal of Geophysical Research: Solid Earth, 125(8):e2019JB018557, 2020. Qiang Qiu, Emma M Hill, Sylvain Barbot, Judith Hubbard, Wanpeng Feng, Eric O Lindsey, Lujia Feng, Keren Dai, Sergey V Samsonov, and Paul Tapponnier. The mechanism of partial rupture of a locked megathrust: The role of fault morphology. Geology, 44(10):875–878, 2016a. MM Mendoza, A Ghosh, MS Karplus, SL Klemperer, SN Sapkota, LB Adhikari, and A Velasco. Duplex in the Main Himalayan Thrust illuminated by aftershocks of the 2015 Mw 7.8 Gorkha earthquake. Nature Geoscience, pages 1–5, 2019. Xin Wang, Shengji Wei, and Wenbo Wu. Double-ramp on the Main Himalayan Thrust revealed by broad- band waveform modeling of the 2015 Gorkha earthquake sequence. Earth Planet. Sci. Lett., 473:83–93, 2017a. Peter Molnar. Continental tectonics in the aftermath of plate tectonics. Nature, 335(6186):131–137, 1988. Geoffrey F Davies. On the emergence of plate tectonics. Geology, 20(11):963–966, 1992. Jun Korenaga. Initiation and evolution of plate tectonics on earth: theories and observations. Annual review of earth and planetary sciences, 41(1):117–151, 2013. John Suppe et al. Mechanics of mountain building and metamorphism in taiwan. Mem.Geol.Soc.China, 4 (6):67–89, 1981. Teresa E Jordan, BL Isacks, Victor A Ramos, and Richard W Allmendinger. Mountain building in the central andes. Episodes Journal of International Geoscience, 6(3):20–26, 1983. Jean-Philippe Avouac. Mountain building, erosion, and the seismic cycle in the nepal himalaya. Advances in geophysics, 46:1–80, 2003. Luca Dal Zilio, György Hetényi, Judith Hubbard, and Laurent Bollinger. Building the himalaya from tec- tonic to earthquake scales. Nature Reviews Earth & Environment, 2(4):251–268, 2021. Geoffrey King and Geoff Bailey. Tectonics and human evolution. antiquity, 80(308):265–286, 2006. 102 Martin H Trauth, Mark A Maslin, Alan L Deino, Annett Junginger, Moses Lesoloyia, Eric O Odada, Daniel O Olago, Lydia A Olaka, Manfred R Strecker, and Ralph Tiedemann. Human evolution in a variable envi- ronment: the amplifier lakes of eastern africa. Quaternary Science Reviews, 29(23-24):2981–2988, 2010. Geoffrey N Bailey and Geoffrey CP King. Dynamic landscapes and human dispersal patterns: tectonics, coastlines, and the reconstruction of human habitats. Quaternary Science Reviews, 30(11-12):1533–1553, 2011. John L Rich. Mechanics of low-angle overthrust faulting as illustrated by cumberland thrust block, virginia, kentucky, and tennessee. AAPG Bulletin, 18(12):1584–1596, 1934. Robert Ernest Rettger. Experiments on soft-rock deformation. AAPG Bulletin, 19(2):271–292, 1935. David W Stearns and V Matthews. Faulting and forced folding in the rocky mountains foreland. Geological Society of America Memoirs, 151:1–38, 1978. MP Golombek, JB Plescia, and BJ Franklin. Faulting and folding in the formation of planetary wrinkle ridges. In Lunar and Planetary Science Conference Proceedings, volume 21, pages 679–693, 1991. Oliver B Hopkins. OilandGasPossibilitiesoftheHatchetigbeeAnticlineAlabama. US Government Printing Office, 1917. Shankar Mitra. Fault-propagation folds: geometry, kinematic evolution, and hydrocarbon traps. AAPG bulletin, 74(6):921–945, 1990. GM Ingram, TJ Chisholm, CJ Grant, CA Hedlund, P Stuart-Smith, and J Teasdale. Deepwater north west borneo: hydrocarbon accumulation in an active fold and thrust belt. Marine and Petroleum Geology, 21 (7):879–887, 2004. SF Cox, VJ Wall, MA Etheridge, and TF Potter. Deformational and metamorphic processes in the forma- tion of mesothermal vein-hosted gold deposits—examples from the lachlan fold belt in central victoria, australia. Ore geology reviews, 6(5):391–423, 1991. Alexander FM Kisters. Controls of gold-quartz vein formation during regional folding in amphibolite- facies, marble-dominated metasediments of the navachab gold mine in the pan-african damara belt, namibia. South African Journal of Geology, 108(3):365–380, 2005. Richard W Allmendinger and John H Shaw. Estimation of fault propagation distance from fold shape: Implications for earthquake hazard assessment. Geology, 28(12):1099–1102, 2000. Wen-Shan Chen, Bor-Shouh Huang, Yue-Gau Chen, Yuan-Hsi Lee, Chao-Nan Yang, Ching-Hua Lo, Hui- Cheng Chang, Quo-Cheng Sung, Neng-Wei Huang, Chin-Cheng Lin, et al. 1999 chi-chi earthquake: a case study on the role of thrust-ramp structures for generating earthquakes. Bull.Seism.Soc.Am., 91(5): 986–994, 2001. David K Keefer. Landslides caused by earthquakes. Geological Society of America Bulletin, 95(4):406–421, 1984. Muhammad Basharat, Muhammad Tayyib Riaz, M Qasim Jan, Chong Xu, and Saima Riaz. A review of landslides related to the 2005 kashmir earthquake: implication and future challenges. Natural Hazards, 108(1):1–30, 2021. Qiang Qiu and Sylvain Barbot. Tsunami excitation in the outer wedge of global subduction zones. Earth- Science Reviews, page 104054, 2022. 103 Bruce Hobbs, Klaus Regenauer-Lieb, and Alison Ord. Thermodynamics of folding in the middle to lower crust. Geology, 35(2):175–178, 2007. John P Platt, Whitney M Behr, and Frances J Cooper. Metamorphic core complexes: Windows into the mechanics and rheology of the crust. Journal of the Geological Society, 172(1):9–27, 2015. John P Platt, Haoran Xia, and William Lamborn Schmidt. Rheology and stress in subduction zones around the aseismic/seismic transition. Progress in Earth and Planetary Science, 5(1):1–12, 2018. Christian Brandes and David C Tanner. Fault-related folding: A review of kinematic models and their application. Earth-Science Reviews, 138:352–370, 2014. Seyed Tohid Nabavi and Haakon Fossen. Fold geometry and folding–a review. Earth-ScienceReviews, 222: 103812, 2021. Robert S Yeats. Active faults related to folding. Active tectonics, pages 63–79, 1986. Ben A Van der Pluijm and Stephen Marshak. Earth Structure, volume 149. Norton ISBN: 0-393-92476-X 4, 2005. Haakon Fossen. Structural geology. Cambridge university press, 2016. J. Suppe. Mechanics of mountain building and metamorphism in Taiwan. Mem.Geol.Soc.China, (4):67–89, 1981. P Tapponnier, B Meyer, J Ph Avouac, G Peltzer, Y Gaudemer, Guo Shunmin, Xiang Hongfa, Yin Kelun, Chen Zhitai, Cai Shuahua, et al. Active thrusting and folding in the qilian shan, and decoupling between upper crust and mantle in northeastern tibet. Earth Planet. Sci. Lett., 97(3-4):382–403, 1990. R Cattin and JP Avouac. Modeling mountain building and the seismic cycle in the himalaya of nepal. J. Geophys. Res., 105(B6):13389–13407, 2000. R. Lacassin, F. Valli, N. Arnaud, P. H. Leloup, J. L. Paquette, L. Haibing, P. Tapponnier, M.-L. Chevalier, S. Guillot, G. Maheo, and X. Zhiqin. Large-scale geometry, offset and kinematic evolution of the Karako- rum fault, Tibet. Earth Planet. Sci. Lett., 219(3-4):255–269, 2004. B Biju-Duval, P Le Quellec, A Mascle, V Renard, and P Valery. Multibeam bathymetric survey and high resolution seismic investigations on the barbados ridge complex (eastern caribbean): A key to the knowl- edge and interpretation of an accretionary wedge. Tectonophysics, 86(1-3):275–304, 1982. SPS Gulick, NLB Bangs, TH Shipley, Y Nakamura, G Moore, and S Kuramoto. Three-dimensional architec- ture of the nankai accretionary prism’s imbricate thrust zone off cape muroto, japan: Prism reconstruc- tion via en echelon thrust propagation. Journal of Geophysical Research: Solid Earth, 109(B2), 2004. Kyle Bradley, Yanfang Qin, Hélène Carton, Nugroho Hananto, Fernando Villanueva-Robles, Fréderique Leclerc, Wei Shengji, Paul Tapponier, Kerry Sieh, and Satish Singh. Stratigraphic control of frontal dé- collement level and structural vergence, and implications for tsunamigenic earthquake hazard in suma- tra, indonesia. Geochemistry, Geophysics, Geosystems, 2019. Ronald E Wilcox, TP t Harding, and DR Seely. Basic wrench tectonics. Aapg Bulletin, 57(1):74–96, 1973. Arthur G Sylvester. Strike-slip faults. Geological Society of America Bulletin, 100(11):1666–1703, 1988. 104 Sarah E Tindall and GH Davis. Monocline development by oblique-slip fault-propagation folding: the east kaibab monocline, colorado plateau, utah. Journal of Structural Geology, 21(10):1303–1320, 1999. Haakon Fossen and Timothy B Holst. Northwest-verging folds and the northwestward movement of the caledonian jotun nappe, norway. Journal of Structural Geology, 17(1):3–15, 1995. Keith A Howard and Barbara E John. Fault-related folding during extension: Plunging basement-cored folds in the basin and range. Geology, 25(3):223–226, 1997. Ian R Sharp, Rob L Gawthorpe, John R Underhill, and Sanjeev Gupta. Fault-propagation folding in exten- sional settings: Examples of structural style and synrift sedimentary response from the suez rift, sinai, egypt. Geological Society of America Bulletin, 112(12):1877–1899, 2000. Lyal B Harris, Hemin A Koyi, and Haakon Fossen. Mechanisms for folding of high-grade rocks in exten- sional tectonic settings. Earth-Science Reviews, 59(1-4):163–210, 2002. David A Ferrill, Alan P Morris, and Ronald N McGinnis. Extensional fault-propagation folding in me- chanically layered rocks: The case against the frictional drag mechanism. Tectonophysics, 576:78–85, 2012. William R Jamison. Geometric analysis of fold development in overthrust terranes. Journal of structural Geology, 9(2):207–219, 1987. Josep Poblet and Hardy Stuart. Reverse modelling of detachment folds; application to the pico del aguila anticline in the south central pyrenees (spain). Journal of Structural Geology, 17(12):1707–1724, 1995. Josep Poblet, Ken McClay, Fabrizio Storti, and Josep Anton Muñoz. Geometries of syntectonic sediments associated with single-layer detachment folds. Journal of Structural Geology, 19(3-4):369–381, 1997. Mark G Rowan. Three-dimensional geometry and evolution of a segmented detachment fold, mississippi fan foldbelt, gulf of mexico. Journal of Structural Geology, 19(3-4):463–480, 1997. Shankar Mitra. A unified kinematic model for the evolution of detachment folds. Journal of Structural Geology, 25(10):1659–1673, 2003. KM Scharer, DW Burbank, J Chen, RJ Weldon, Charles Rubin, R Zhao, and J Shen. Detachment folding in the southwestern tian shan–tarim foreland, china: shortening estimates and rates. Journal of Structural Geology, 26(11):2119–2137, 2004. Jianjun Li and Shankar Mitra. Seismic models of detachment and faulted detachment folds. Marine and Petroleum Geology, 117:104385, 2020. John Suppe and Donald A Medwedeff. Geometry and kinematics of fault-propagation folding. Eclogae Geologicae Helvetiae, 83(3):409–454, 1990. John H Shaw, Christopher David Connors, John Suppe, et al. Seismic interpretation of contractional fault- related folds: An AAPG seismic atlas, volume 53. American Association of Petroleum Geologists, 2005. Amanda N Hughes, Nathan P Benesh, and John H Shaw. Factors that control the development of fault-bend versus fault-propagation folds: Insights from mechanical models based on the discrete element method (dem). Journal of Structural Geology, 68:121–141, 2014. Amanda N Hughes and John H Shaw. Insights into the mechanics of fault-propagation folding styles. Bulletin, 127(11-12):1752–1765, 2015. 105 John Suppe. Geometry and kinematics of fault-bend folding. Am. J. science, 283(7):684–721, 1983. Fred A Donath and Ronald B Parker. Folds and folding. GeologicalSocietyofAmericaBulletin, 75(1):45–62, 1964. PW Geoff Tanner. The flexural-slip mechanism. Journal of Structural Geology, 11(6):635–655, 1989. Gary D Couples, Helen Lewis, and PW Geoff Tanner. Strain partitioning during flexural-slip folding. Geological Society, London, Special Publications, 127(1):149–165, 1998. Judith S Chester, John M Logan, and John H Spang. Influence of layering and boundary conditions on fault-bend and fault-propagation folding. GeologicalSocietyofAmericaBulletin, 103(8):1059–1072, 1991. Davi R Damasceno, Andreas Eckert, and Xiaolong Liu. Flexural-slip during visco-elastic buckle folding. Journal of Structural Geology, 100:62–76, 2017. Kaj M Johnson. Growth of fault-cored anticlines by flexural slip folding: Analysis by boundary element modeling. Journal of Geophysical Research: Solid Earth, 123(3):2426–2447, 2018. JOHN Suppe and J Namson. Fault-bend origin of frontal folds of the western taiwan fold-and-thrust belt. Petroleum Geology of Taiwan, 16:1–18, 1979. Yue-Gau Chen, Kuang-Yin Lai, Yuan-Hsi Lee, John Suppe, Wen-Shan Chen, Yu-Nong N Lin, Yu Wang, Jih-Hao Hung, and Yu-Ting Kuo. Coseismic fold scarps and their kinematic behavior in the 1999 chi-chi earthquake taiwan. J. Geophys. Res., 112(B3), 2007. Janet T Watt and Daniel S Brothers. Systematic characterization of morphotectonic variability along the cascadia convergent margin: Implications for shallow megathrust behavior and tsunami hazards. Geo- sphere, 17(1):95–117, 2021. Lionel Endignoux and Jean-Louis Mugnier. The use of a forward kinematic model in the construction of balanced cross sections. Tectonics, 9(5):1249–1262, 1990. James W McDougall and Ahmad Hussain. Fold and thrust propagation in the western himalaya based on a balanced cross section of the surghar range and kohat plateau, pakistan. AAPG bulletin, 75(3):463–478, 1991. DILIP K Mukhopadhyay and PREMANAND Mishra. The main frontal thrust (mft), northwestern hi- malayas: Thrust trajectory and hanging wall fold geometry from balanced cross sections. J. Geol. Soc. India, 64(6):739–746, 2004. Rafael V Almeida, Judith Hubbard, Lee Liberty, Anna Foster, and Soma Nath Sapkota. Seismic imaging of the main frontal thrust in nepal reveals a shallow décollement and blind thrusting. Earth and Planetary Science Letters, 494:216–225, 2018. J Lavé and J-Ph Avouac. Active folding of fluvial terraces across the siwaliks hills, himalayas of central nepal. J. Geophys. Res., 105(B3):5735–5770, 2000. L. Bollinger, Y. Klinger, P. Tapponnier, Y. Gaudemer, and D. Tiwari. Estimating the return times of great Himalayan earthquakes in eastern Nepal: Evidence from the Patu and Bardibas strands of the Main Frontal Thrust. Nature Geosci., 6:71–76, 2013. Jean-Philippe Avouac. From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. Annual Review of Earth and Planetary Sciences, 43:233–271, 2015. 106 Sylvain Barbot, Nadia Lapusta, and Jean-Philippe Avouac. Under the hood of the earthquake machine: Toward predictive modeling of the seismic cycle. Science, 336(6082):707–710, 2012a. Hugo Perfettini and JP Avouac. The seismic cycle in the area of the 2011 mw9. 0 tohoku-oki earthquake. Journal of Geophysical Research: Solid Earth, 119(5):4469–4515, 2014. Y Van Dinther, TV Gerya, LA Dalguer, Paul Martin Mai, G Morra, and D Giardini. The seismic cycle at subduction thrusts: Insights from seismo-thermo-mechanical models. Journal of Geophysical Research: Solid Earth, 118(12):6183–6202, 2013. Luca Dal Zilio, Ylona van Dinther, Taras Gerya, and Jean-Philippe Avouac. Bimodal seismicity in the himalaya controlled by fault friction and geometry. Nature communications, 10(1):48, 2019. Qibin Shi, Sylvain Barbot, Shengji Wei, Paul Tapponnier, Takanori Matsuzawa, and Bunichiro Shibazaki. Structural control and system-level behavior of the seismic cycle at the nankai trough. Earth, Planets and Space, 72(1):1–31, 2020a. J. H. Dieterich. Modeling of rock friction 1. experimental results and constitutive equations. J. Geophys. Res., 84(B5):2161–2168, 1979. A. Ruina. Slip instability and state variable friction laws. J. Geophys. Res., 88:10,359–10,370, 1983. J. R. Rice and A. L. Ruina. Stability of steady frictional slipping. J. Appl. Mech., 50:343–349, 1983. Duo Li and Yajing Liu. Spatiotemporal evolution of slow slip events in a nonplanar fault model for northern cascadia subduction zone. J. Geophys. Res., 121(9):6828–6845, 2016. Y. Liu and J. R. Rice. Aseismic slip transients emerge spontaneously in three-dimensional rate and state modeling of subduction earthquake sequences. J. Geophys. Res., 110(B08307), 2005. doi: 10.1029/ 2004JB003424. Y. Liu and J. R. Rice. Spontaneous and triggered aseismic deformation transients in a subduction fault model. J. Geophys. Res., 112(B09404), 2007. doi: 10.1029/2007JB004930. A. M. Rubin. Episodic slow slip events and rate-and-state friction. J.Geophys.Res., 113(B11414), 2008. doi: 10.1029/2008JB005642. T. Matsuzawa, H. Hirose, B. Shibazaki, and K. Obara. Modeling short- and long-term slow slip events in the seismic cycles of large subduction earthquakes. J. Geophys. Res., 115(B12301), 2010. doi: 10.1029/ 2010JB007566. P. Segall, A. M. Rubin, A. M. Bradley, and J. R. Rice. Dilatant strengthening as a mechanism for slow slip events. J. Geophys. Res., 115(B12305), 2010. M. Wei, Y. Kaneko, Y. Liu, and J. J McGuire. Episodic fault creep events in california controlled by shallow frictional heterogeneity. Nature geoscience, 6(7):566, 2013. Meng Wei, Yajing Liu, Yoshihiro Kaneko, Jeffrey J McGuire, and Roger Bilham. Dynamic triggering of creep events in the salton trough, southern california by regional m≥5.4 earthquakes constrained by geodetic observations and numerical simulations. Earth Planet. Sci. Lett., 427:1–10, 2015. Duo Li and Yajing Liu. Modeling slow-slip segmentation in cascadia subduction zone constrained by tremor locations and gravity anomalies. J. Geophys. Res., 122(4):3138–3157, 2017. 107 A. Goswami and S. Barbot. Slow-slip events in semi-brittle serpentinite fault zones. Scientificreports, 8(1): 6181, 2018. Meng Wei, Yoshihiro Kaneko, Pengcheng Shi, and Yajing Liu. Numerical modeling of dynamically triggered shallow slow slip events in new zealand by the 2016 mw 7.8 kaikoura earthquake. Geophys. Res. Lett., 45(10):4764–4772, 2018. Pierre Romanet, Harsha S Bhat, Romain Jolivet, and Raúl Madariaga. Fast and slow slip events emerge due to fault geometrical complexity. Geophysical Research Letters, 45(10):4809–4819, 2018. Sylvain Barbot. Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. Tectonophysics, 768:228171, 2019a. Shiying Nie and Sylvain Barbot. Seismogenic and tremorgenic slow slip near the stability transition of frictional sliding. Earth and Planetary Science Letters, 569:117037, 2021. Paul G Okubo. Dynamic rupture modeling with laboratory-derived constitutive relations. Journal of Geophysical Research: Solid Earth, 94(B9):12321–12335, 1989. Nadia Lapusta, James R Rice, Yehuda Ben-Zion, and Gutuan Zheng. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction. Journal of Geophysical Research: Solid Earth, 105(B10):23765–23789, 2000a. Andrea Bizzarri and Massimo Cocco. Slip-weakening behavior during the propagation of dynamic ruptures obeying rate-and state-dependent friction laws. Journal of Geophysical Research: Solid Earth, 108(B8), 2003. Elisa Tinti, Andrea Bizzarri, and Massimo Cocco. Modeling the dynamic rupture propagation on hetero- geneous faults with rate-and state-dependent friction. Annals of Geophysics, 2005. A-A Gabriel, J-P Ampuero, Luis A Dalguer, and Paul Martin Mai. The transition of dynamic rupture styles in elastic media under velocity-weakening friction. J. Geophys. Res., 117(B9), 2012a. James R Rice and Simon T Tse. Dynamic motion of a single degree of freedom system following a rate and state dependent friction law. J. Geophys. Res., 91(B1):521–530, 1986. N. Lapusta and J. R. Rice. Nucleation and early seismic propagation of small and large events in a crustal earthquake model. J. Geophs. Res., 108(B4, 2205), 2003. Takane Hori, Naoyuki Kato, Kazuro Hirahara, Toshitaka Baba, and Yoshiyuki Kaneda. A numerical sim- ulation of earthquake cycles along the nankai trough in southwest japan: lateral variation in frictional property due to the slab geometry controls the nucleation position. Earth Planet. Sci. Lett., 228(3-4): 215–226, 2004. M. Mele Veedu and S. Barbot. The Parkfield tremors reveal slow and fast ruptures on the same asperity. Nature, 532(7599):361–365, 2016. doi: 10.1038/nature17190. Junle Jiang and Nadia Lapusta. Deeper penetration of large earthquakes on seismically quiescent faults. Science, 352(6291):1293–1297, 2016. JEAN-P Avouac, P Tapponnier, M Bai, H You, and G Wang. Active thrusting and folding along the northern tien shan and late cenozoic rotation of the tarim relative to dzungaria and kazakhstan. J. Geophys. Res., 98(B4):6755–6804, 1993. 108 John H Shaw, Enrique Novoa, and Christopher D Connors. Structural controls on growth stratigraphy in contractional fault-related folds. 2004. Kuang-Yin Lai, Yue-Gau Chen, Jih-Hao Hung, John Suppe, Li-Fan Yue, and Ya-Wen Chen. Surface defor- mation related to kink-folding above an active fault: Evidence from geomorphic features and co-seismic slips. Quaternary International, 147(1):44–54, 2006. T Hossler, L Bollinger, SN Sapkota, J Lavé, RM Gupta, and TP Kandel. Surface ruptures of large himalayan earthquakes in western nepal: Evidence along a reactivated strand of the main boundary thrust. Earth Planet. Sci. Lett., 434:187–196, 2016. Dan M Baker, Robert J Lillie, Robert S Yeats, Gary D Johnson, Mohammad Yousuf, and Agha Sher Hamid Zamin. Development of the himalayan frontal thrust zone: Salt range, pakistan. Geology, 16(1):3–7, 1988. G Fz Moore, TH Shipley, PL Stoffa, DE Karig, A Taira, S Kuramoto, H Tokuyama, and K Suyehiro. Structure of the nankai trough accretionary zone from multichannel seismic reflection data. J. Geophys. Res., 95 (B6):8753–8765, 1990. Julien Charreau, Jean-Philippe Avouac, Yan Chen, Stéphane Dominguez, and Stuart Gilder. Miocene to present kinematics of fault-bend folding across the huerguosi anticline, northern tianshan (china), de- rived from structural, seismic, and magnetostratigraphic data. Geology, 36(11):871–874, 2008. Satish C Singh, Hélène Carton, Paul Tapponnier, Nugroho D Hananto, Ajay PS Chauhan, Djoko Hartoyo, Martin Bayly, Soelistijani Moeljopranoto, Tim Bunting, Phil Christie, et al. Seismic evidence for broken oceanic crust in the 2004 sumatra earthquake epicentral region. Nature Geoscience, 1(11):777, 2008. David J Sanderson. Models of strain variation in nappes and thrust sheets: a review. Tectonophysics, 88 (3-4):201–233, 1982. Paul Betka, Keith Klepeis, and Sharon Mosher. Along-strike variation in crustal shortening and kinematic evolution of the base of a retroarc fold-and-thrust belt: Magallanes, chile 53° s–54° s. Bulletin, 127(7-8): 1108–1134, 2015. Paul M Betka, Leonardo Seeber, Stuart N Thomson, Michael S Steckler, Ryan Sincavage, and C Zoramthara. Slip-partitioning above a shallow, weak décollement beneath the indo-burman accretionary prism. Earth and Planetary Science Letters, 503:17–28, 2018. J-E Hurtrez, F Lucazeau, J Lavé, and J-P Avouac. Investigation of the relationships between basin mor- phology, tectonic uplift, and denudation from the study of an active fold belt in the siwalik hills, central nepal. J. Geophys. Res., 104(B6):12779–12796, 1999. Judith Hubbard and John H Shaw. Uplift of the longmen shan and tibetan plateau, and the 2008 wenchuan (m= 7.9) earthquake. Nature, 458(7235):194, 2009. John Suppe, George T. Chou, and Stephen C. Hook. Rates of folding and faulting determined from growth strata, pages 105–121. Springer Netherlands, 1992. Maryline Le Béon, John Suppe, Manoj K Jaiswal, Yue-Gau Chen, and Michaela E Ustaszewski. Deciphering cumulative fault slip vectors from fold scarps: Relationships between long-term and coseismic deforma- tions in central western taiwan. J. Geophys. Res., 119(7):5943–5978, 2014. 109 Alex Copley. Postseismic afterslip 30 years after the 1978 tabas-e-golshan (iran) earthquake: observations and implications for the geological evolution of thrust belts. Geophys. J. Int., 197(2):665–679, 2014. J-L Epard and RH Groshong Jr. Kinematic model of detachment folding including limb rotation, fixed hinges and layer-parallel strain. Tectonophysics, 247(1-4):85–103, 1995. Josep Poblet and Ken McClay. Geometry and kinematics of single-layer detachment folds. AAPG bulletin, 80(7):1085–1109, 1996. Karl Mueller and John Suppe. Growth of wheeler ridge anticline, california: geomorphic evidence for fault-bend folding behaviour during earthquakes. Journal of Structural Geology, 19(3-4):383–396, 1997. BJ Souter and BH Hager. Fault propagation fold growth during the 1994 northridge, california, earthquake? Journal of Geophysical Research: Solid Earth, 102(B6):11931–11942, 1997. Ben A van der Pluijm, Chris M Hall, Peter J Vrolijk, David R Pevear, and Michael C Covey. The dating of shallow faults in the earth’s crust. Nature, 412(6843):172, 2001. Scott R Miller and Rudy L Slingerland. Topographic advection on fault-bend folds: Inheritance of valley positions and the formation of wind gaps. Geology, 34(9):769–772, 2006. Sylvain Bernard, Jean-Philippe Avouac, Stéphane Dominguez, and Martine Simoes. Kinematics of fault- related folding derived from a sandbox experiment. J. Geophys. Res., 112(B3), 2007. Yukinobu Okamura, Tatsuya Ishiyama, and Yukio Yanagisawa. Fault-related folds above the source fault of the 2004 mid-niigata prefecture earthquake, in a fold-and-thrust belt caused by basin inversion along the eastern margin of the japan sea. J. Geophys. Res., 112(B3), 2007. Maomao Wang, Dong Jia, John H Shaw, Judith Hubbard, Aiming Lin, Yiquan Li, and Li Shen. Active fault- related folding beneath an alluvial terrace in the southern longmen shan range front, sichuan basin, china: Implications for seismic hazard. Bull. Seism. Soc. Am., 103(4):2369–2385, 2013. Adam M Forte, Eric Cowgill, Ibrahim Murtuzayev, Talat Kangarli, and Marius Stoica. Structural geome- tries and magnitude of shortening in the eastern kura fold-thrust belt, azerbaijan: Implications for the development of the greater caucasus mountains. Tectonics, 32(3):688–717, 2013. Donald A Medwedeff and John Suppe. Multibend fault-bend folding. JournalofStructuralGeology, 19(3-4): 279–292, 1997. John Rodgers. Mechanics of appalachian folding as illustrated by sequatchie anticline, tennessee and alabama. AAPG Bulletin, 34(4):672–681, 1950. J Ph Avouac, B Meyer, and P Tapponnier. On the growth of normal faults and the existence of flats and ramps along the el asnam active fold and thrust system. Tectonics, 11(1):1–11, 1992. K. Sieh, D. H. Natawidjaja, A. J. Meltzner, C.-C. Shen, H. Cheng, K.-S. Li, B. W. Suwargadi, J. Galetzka, B. Philibosian, and R. L. Edwards. Earthquake supercycles inferred from sea-level changes recorded in the corals of West Sumatra. Science, 322:1674–1678, 2008. Belle Philibosian, Kerry Sieh, Jean-Philippe Avouac, Danny H Natawidjaja, Hong-Wei Chiang, Chung-Che Wu, Chuan-Chou Shen, Mudrik R Daryono, Hugo Perfettini, Bambang W Suwargadi, et al. Earthquake supercycles on the mentawai segment of the sunda megathrust in the seventeenth century and earlier. J. Geophys. Res., 122(1):642–676, 2017. 110 K. Ariyoshia, T. Hori, J.-P. Ampuero, Y. Kaneda, T. Matsuzawa, R. Hino, and A. Hasegawa. Influence of interaction between small asperities on various types of slow earthquakes in a 3-d simulation for a subduction plate boundary. Gonwana Res., 16(3-4):534–544, 2009. Y. Kaneko, J.-P. Avouac, and N. Lapusta. Towards inferring earthquake patterns from geodetic observations of interseismic coupling. Nature Geoscience, 3:363–369, 2010. Pierre Dublanchet, Pascal Bernard, and Pascal Favreau. Interactions and triggering in a 3-d rate-and-state asperity model. J. Geophys. Res., 118(5):2225–2245, 2013. Hiroyuki Noda and Nadia Lapusta. Stable creeping fault segments can become destructive as a result of dynamic weakening. Nature, 493(7433):518–521, 2013. Benchun Duan and David D Oglesby. Multicycle dynamics of nonplanar strike-slip faults. J.Geophys.Res., 110(B3), 2005. J. C. Lozos, D. D. Oglesby, B. Duan, and S. G. Wesnousky. The effects of double fault bends on rupture propagation: A geometrical parameter study. Bull. Seism. Soc. Am., 101(1):385–398, 2011. Q. Qiu, E. M Hill, S. Barbot, J. Hubbard, W. Feng, E. O Lindsey, L. Feng, K. Dai, S. V Samsonov, and P. Tappon- nier. The mechanism of partial rupture of a locked megathrust: The role of fault morphology. Geology, 44(10):875–878, 2016b. Yann Klinger, R Michel, and GCP King. Evidence for an earthquake barrier model from mw 7.8 kokoxili (tibet) earthquake slip-distribution. Earth Planet. Sci. Lett., 242(3-4):354–364, 2006. Sylvain Michel, Jean-Philippe Avouac, Nadia Lapusta, and Junle Jiang. Pulse-like partial ruptures and high- frequency radiation at creeping-locked transition during megathrust earthquakes. Geophys. Res. Lett., 44(16):8345–8351, 2017. Yan Wu and Xiaofei Chen. The scale-dependent slip pattern for a uniform fault model obeying the rate-and state-dependent friction law. J. Geophys. Res., 119(6):4890–4906, 2014. H. Noda and T. Hori. Under what circumstances does a seismogenic patch produce aseismic transients in the later interseismic period? Geophys. Res. Lett., 2014. Sylvain Barbot. Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. Tectonophysics, 2019b. Benchun Duan and David D Oglesby. Heterogeneous fault stresses from previous earthquakes and the effect on dynamics of parallel strike-slip faults. J. Geophys. Res., 111(B5), 2006. Yuta Mitsui. Elastic interaction of parallel rate-and-state-dependent frictional faults with aging and slip laws: slow-slip faults can sometimes host fast events. Earth, Planets and Space, 70(1):136, 2018. Geoffrey King and John Nábělek. Role of fault bends in the initiation and termination of earthquake rupture. Science, 228(4702):984–987, 1985. GCP King. Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure. Pure Appl. Geophys., 124(3):567–585, 1986. Hideo Aochi, Raúl Madariaga, and Eiichi Fukuyama. Effect of normal stress during rupture propagation along nonplanar faults. Journal of Geophysical Research: Solid Earth, 107(B2):ESE–5, 2002. 111 Y Kase and SM Day. Spontaneous rupture processes on a bending fault. Geophys. Res. Lett., 33(10), 2006. Benchun Duan and Steven M Day. Inelastic strain distribution and seismic radiation from rupture of a fault kink. J. Geophys. Res., 113(B12), 2008. Zhenguo Zhang, Wei Zhang, and Xiaofei Chen. Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics. Geophysical Journal International, 199(2):860–879, 2014. Simon Preuss, Robert Herrendörfer, Taras Gerya, Jean Paul Ampuero, and Ylona van Dinther. Seismic and aseismic fault growth lead to different fault orientations. Journal of Geophysical Research: Solid Earth, 2019. Iris van Zelst, Stephanie Wollherr, Alice-Agnes Gabriel, Elizabeth H. Madden, and Ylona van Dinther. Modeling Megathrust Earthquakes Across Scales: One-way Coupling From Geodynamics and Seismic Cycles to Dynamic Rupture. Journal of Geophysical Research, (124), 2019. doi: 10.1029/2019JB017539. URL https://doi.org/10.1029/2019JB017539. Alexei NB Poliakov, Renata Dmowska, and James R Rice. Dynamic shear rupture interactions with fault bends and off-axis secondary faulting. Journal of Geophysical Research: Solid Earth, 107(B11):ESE–6, 2002a. Hideo Aochi, Eiichi Fukuyama, and Mitsuhiro Matsu’ura. Spontaneous rupture propagation on a non- planar fault in 3-d elastic medium. pure and applied geophysics, 157(11-12):2003–2027, 2000. Carl-Ernst Rousseau and Ares J Rosakis. Dynamic path selection along branched faults: Experiments involving sub-rayleigh and supershear ruptures. Journal of Geophysical Research: Solid Earth, 114(B8), 2009. Ossian O’Reilly, Jan Nordström, Jeremy E Kozdon, and Eric M Dunham. Simulation of earthquake rupture dynamics in complex geometries using coupled finite difference and finite volume methods. Communi- cations in Computational Physics, 17(2):337–370, 2015. R. Herrendörfer, Y. Van Dinther, T. Gerya, and L. A. Dalguer. Earthquake supercycle in subduction zones controlled by the width of the seismogenic zone. Nature Geoscience, 8(6):471–474, 2015. Miranda Ong, Su Qing, Sylvain Barbot, and Judith Hubbard. Physics-based scenario of earthquake cycles on the Ventura Thrust system, California: the effect of variable friction and fault geometry. Pure Appl. Geophys., 2019. doi: 10.1007/s00024-019-02111-9. Laurent Maerten, Paul Gillespie, and David D Pollard. Effects of local stress perturbation on secondary fault development. Journal of Structural Geology, 24(1):145–153, 2002. Shankar Mitra. Fold-accommodation faults. AAPG bulletin, 86(4):671–693, 2002. S Daout, S Barbot, G Peltzer, M-P Doin, Z Liu, and R Jolivet. Constraining the kinematics of metropolitan los angeles faults with a slip-partitioning model. Geophys. Res. Lett., 43(21), 2016a. Thomas L Davis and Jay S Namson. A balanced cross-section of the 1994 northridge earthquake, southern california. Nature, 372(6502):167, 1994. Lorraine A Leon, Shari A Christofferson, James F Dolan, John H Shaw, and Thomas L Pratt. Earthquake- by-earthquake fold growth above the puente hills blind thrust fault, los angeles, california: Implications for fold kinematics and seismic hazard. J. Geophys. Res., 112(B3), 2007. 112 Thomas L Pratt, John H Shaw, James F Dolan, Shari A Christofferson, Robert A Williams, Jack K Odum, and Andreas Plesch. Shallow seismic imaging of folds above the puente hills blind-thrust fault, los angeles, california. Geophys. Res. Lett., 29(9):18–1, 2002. John H Shaw and John Suppe. Earthquake hazards of active blind-thrust faults under the central los angeles basin, california. J. Geophys. Res., 101(B4):8623–8642, 1996. Judith Hubbard, John H Shaw, James Dolan, Thomas L Pratt, Lee McAuliffe, and Thomas K Rockwell. Structure and seismic hazard of the ventura avenue anticline and ventura fault, california: Prospect for large, multisegment ruptures in the western transverse ranges. Bull. Seism. Soc. Am., 104(3):1070–1087, 2014. Hiroo Kanamori. Mechanism of tsunami earthquakes. Phys. Earth Planet. Inter., 6(5):346–359, 1972. Hiroo Kanamori and Masayuki Kikuchi. The 1992 nicaragua earthquake: a slow tsunami earthquake as- sociated with subducted sediments. Nature, 361(6414):714, 1993. S. L. Bilek and T. Lay. Tsunami earthquakes possibly widespread manifestations of frictional conditional stability. Geophys. Res. Lett., 29(14, 1673):4, 2002. doi: 10.1029/2002GL015215. Demian M Saffer and Chris Marone. Comparison of smectite-and illite-rich gouge frictional properties: application to the updip limit of the seismogenic zone along subduction megathrusts. Earth Planet. Sci. Lett., 215(1-2):219–235, 2003. Matt J Ikari, André R Niemeijer, Christopher J Spiers, Achim J Kopf, and Demian M Saffer. Experimental evidence linking slip instability with seafloor lithology and topography at the costa rica convergent margin. Geology, 41(8):891–894, 2013. Matt J Ikari and Achim J Kopf. Seismic potential of weak, near-surface faults revealed at plate tectonic slip rates. Science advances, 3(11):e1701269, 2017. Jan H Behrmann and Achim Kopf. Balance of tectonically accreted and subducted sediment at the chile triple junction. International Journal of Earth Sciences, 90(4):753–768, 2001. Heidrun Kopp. Invited review paper: The control of subduction zone structural complexity and geometry on margin segmentation and seismicity. Tectonophysics, 589:1–16, 2013. Judith Hubbard, Sylvain Barbot, Emma M Hill, and Paul Tapponnier. Coseismic slip on shallow décollement megathrusts: implications for seismic and tsunami hazard. Earth-Science Reviews, 141:45–55, 2015. Jacob Geersen. Sediment-starved trenches and rough subducting plates are conducive to tsunami earth- quakes. Tectonophysics, 762:28–44, 2019. Thorne Lay, Lingling Ye, Hiroo Kanamori, Yoshiki Yamazaki, Kwok Fai Cheung, Kevin Kwong, and Keith D Koper. The october 28, 2012 mw 7.8 haida gwaii underthrusting earthquake and tsunami: Slip partition- ing along the queen charlotte fault transpressional plate boundary. Earth Planet. Sci. Lett., 375:57–70, 2013. Anthony Sladen and Jenny Trevisan. Shallow megathrust earthquake ruptures betrayed by their outer- trench aftershocks signature. Earth Planet. Sci. Lett., 483:105–113, 2018. Stacey Servito Martin, Linlin Li, Emile A Okal, Julie Morin, Alexander EG Tetteroo, Adam D Switzer, and Kerry E Sieh. Reassessment of the 1907 sumatra “tsunami earthquake” based on macroseismic, seismological, and tsunami observations, and modeling. Pure Appl. Geophys., pages 1–38. 113 Christopher H Scholz. The rupture mode of the shallow large-slip surge of the tohoku-oki earthquake. Bull. Seism. Soc. Am., 104(5):2627–2631, 2014. Charles J Ammon, Hiroo Kanamori, Thorne Lay, and Aaron A Velasco. The 17 july 2006 java tsunami earthquake. Geophys. Res. Lett., 33(24), 2006a. Jeremy E Kozdon and Eric M Dunham. Rupture to the trench: Dynamic rupture simulations of the 11 march 2011 tohoku earthquakerupture to the trench: Dynamic rupture simulations of the 11 march 2011 tohoku earthquake. Bull. Seism. Soc. Am., 103(2B):1275–1289, 2013. Thorne Lay, Hiroo Kanamori, Charles J Ammon, Meredith Nettles, Steven N Ward, Richard C Aster, Susan L Beck, Susan L Bilek, Michael R Brudzinski, Rhett Butler, et al. The great sumatra-andaman earthquake of 26 december 2004. Science, 308(5725):1127–1133, 2005. Charles J Ammon, Chen Ji, Hong-Kie Thio, David Robinson, Sidao Ni, Vala Hjorleifsdottir, Hiroo Kanamori, Thorne Lay, Shamita Das, Don Helmberger, et al. Rupture process of the 2004 sumatra-andaman earth- quake. Science, 308(5725):1133–1139, 2005. Nicole Feldl and Roger Bilham. Great himalayan earthquakes and the tibetan plateau. Nature, 444(7116): 165, 2006. M. Simons, S. Minson, A. Sladen, F. Ortega, J. Jiang, S. E. Owen, L. Meng, J.-P. Ampuero, S. Wei, R. Chu, D. V. Helmberger, H. Kanamori, E. Hetland, A. W. Moore, and F. H. Webb. The 2011 Magnitude 9.0 Tohoku-Oki Earthquake: Mosaicking the Megathrust from Seconds to Centuries. Science, 332(6036): 1421–1425, 2012. J-L Mugnier, A Gajurel, P Huyghe, R Jayangondaperumal, F Jouanne, and B Upreti. Structural interpreta- tion of the great earthquakes of the last millennium in the central himalaya. Earth-ScienceReviews, 127: 30–47, 2013. Martin Vallee, Jean-Mathieu Nocquet, Jean Battaglia, Yvonne Font, Monica Segovia, Marc Regnier, Patricia Mothes, Paul Jarrin, David Cisneros, Sandro Vaca, et al. Intense interface seismicity triggered by a shallow slow slip event in the central ecuador subduction zone. J.Geophys.Res., 118(6):2965–2981, 2013. Laura M Wallace, Yoshihiro Kaneko, Sigrún Hreinsdóttir, Ian Hamling, Zhigang Peng, Noel Bartlow, Elis- abetta D’Anastasio, and Bill Fry. Large-scale dynamic triggering of shallow slow slip enhanced by over- lying sedimentary wedge. Nature Geoscience, 2017. Akiko Toh, Koichiro Obana, and Eiichiro Araki. Distribution of very low frequency earthquakes in the nankai accretionary prism influenced by a subducting-ridge. Earth Planet. Sci. Lett., 482:342–356, 2018. Masaru Nakano, Takane Hori, Eiichiro Araki, Shuichi Kodaira, and Satoshi Ide. Shallow very-low- frequency earthquakes accompany slow slip events in the nankai subduction zone. Nature commu- nications, 9(1):984, 2018. Jin-Oh Park, Tetsuro Tsuru, Shuichi Kodaira, Phil R Cummins, and Yoshiyuki Kaneda. Splay fault branching along the nankai subduction zone. Science, 297(5584):1157–1160, 2002. A Cox and RB Hart. Plate tectonics: How it works: Blackwell scientific publications. Inc. Boston, 1986. Dan Davis, John Suppe, and FA Dahlen. Mechanics of fold-and-thrust belts and accretionary wedges. J. Geophys. Res., 88(B2):1153–1172, 1983. 114 Heidrun Kopp and Nina Kukowski. Backstop geometry and accretionary mechanics of the sunda margin. Tectonics, 22(6), 2003. Demian M Saffer and Barbara A Bekins. Hydrologic controls on the morphology and mechanics of accre- tionary wedges. Geology, 30(3):271–274, 2002. S. Barbot, N. Lapusta, and J. P. Avouac. Under the hood of the earthquake machine: Towards predictive modeling of the seismic cycle. Science, 336(6082):707–710, 2012b. N. Lapusta and S. Barbot. Models of earthquakes and aseismic slip based on laboratory-derived rate and state friction laws. In A. Bizzarri and H. S. Bhat, editors, The mechanics of Faulting: From Laboratory to Real Earthquakes, pages 153–207. Research Signpost, Trivandrum, Kerala, India, 2012. Camilla Cattania and Paul Segall. Crack models of repeating earthquakes predict observed moment- recurrence scaling. J. Geophys. Res., 2018. S. Barbot and Y. Fialko. Fourier-domain Green’s function for an elastic semi-infinite solid under gravity, with applications to earthquake and volcano deformation. Geophys. J. Int., 182(2):568–582, 2010a. doi: 10.1111/j.1365-246X.2010.04655.x. S. Barbot, J. DP Moore, and V. Lambert. Displacement and stress associated with distributed anelastic deformation in a half-space. Bull. Seism. Soc. Am., 107(2):821–855, 2017. doi: 10.1785/0120160237. Sylvain Barbot. Deformation of a half-space from anelastic strain confined in a tetrahedral volume. Bull. Seism. Soc. Am., 108(5A):2687, 2018. doi: 10.1785/0120180058. URL http://dx.doi.org/10.1785/ 0120180058. DJ Andrews. Coupling of energy between tectonic processes and earthquakes. J. Geophys. Res., 83(B5): 2259–2264, 1978. S. Nemat-Nasser and M. Hori. Micromechanics: overall properties of heterogeneous materials. Elsevier, 2 edition, 1999. S. Nemat-Nasser.Plasticity.ATreatiseonFiniteDeformationofHeterogeneousInelasticMaterials. Cambridge Monographs on Mechanics. Cambridge University Press, 2004. S. Barbot and Y. Fialko. A unified continuum representation of postseismic relaxation mechanisms: semi- analytic models of afterslip, poroelastic rebound and viscoelastic flow. Geophys.J.Int., 182(3):1124–1140, 2010b. doi: 10.1111/j.1365-246X.2010.04678.x. J. A. Steketee. On Volterra’s Dislocations in a Semi-Infinite Elastic Medium. Can. J. Phys., 36:192–205, 1958a. J. A. Steketee. Some geophysical applications of the elasticity theory of dislocations. Can. J. Phys., 36: 1168–1198, 1958b. Y. Okada. Internal deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am., 82: 1018–1040, April 1992. R.J. Wang, F.L. Martin, and F. Roth. Computation of deformation induced by earthquakes in a multi-layered elastic crust - FORTRAN programs EDGRN/EDCMP. Comp. Geosci., 29:195–207, 2003. B Smith and D. Sandwell. A three-dimensional semianalytic viscoelastic model for time-dependent analy- ses of the earthquake cycle. J. Geophys. Res., 109, 2004. 115 S Daout, R Jolivet, C Lasserre, M-P Doin, S Barbot, P Tapponnier, G Peltzer, A Socquet, and J Sun. Along- strike variations of the partitioning of convergence across the haiyuan fault system detected by insar. Geophys. J. Int., 205(1):536–547, 2016b. John H Shaw and John Suppe. Active faulting and growth folding in the eastern santa barbara channel, california. Geological Society of America Bulletin, 106(5):607–626, 1994. Mark G Rowan and Roberto Linares. Fold-evolution matrices and axial-surface analysis of fault-bend folds: application to the medina anticline, eastern cordillera, colombia. AAPG bulletin, 84(6):741–764, 2000. MA Chinnery. The deformation of the ground around surface faults. Bull. Seism. Soc. Am., 51(3):355–372, 1961. J. C. Savage and L. M. Hastie. Surface deformation associated with dip-slip faulting. J. Geophys. Res., 71 (20):4897–4904, 1966. R. Sato and M. Matsu’ura. Strains and tilts on the surface of a semi-infinite medium. J. Phys. Earth, 22(2): 213–221, 1974. Y. Okada. Surface deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am., 75(4): 1135–1154, August 1985. Mehdi Nikkhoo and Thomas R Walter. Triangular dislocation: an analytical, artefact-free solution. Geo- phys. J. Int., 201(2):1117–1139, 2015. H Nozaki and M Taya. Elastic fields in a polyhedral inclusion with uniform eigenstrains and related problems. Journal of Applied mechanics, 68(3):441–452, 2001. Boris N Kuvshinov. Elastic and piezoelectric fields due to polyhedral inclusions. International Journal of Solids and Structures, 45(5):1352–1384, 2008. J. H. Dieterich. Time-dependent friction and the mechanics of stick-slip. Pure Appl. Geophys., 116(4-5): 790–806, 1978. C. Marone, B. C. Raleigh, and C. H. Scholz. Frictional behavior and constitutive modeling of simulated fault gouge. J. Geophys. Res., 95(B5):7007–7026, May 1990. C. Marone and B. Kilgore. Scaling of the critical slip distance for seismic faulting with shear strain in fault zones. Nature, 362:618–620, 1993. M. L. Blanpied, D. A. Lockner, and J. D. Byerlee. Frictional slip of granite at hydrothermal conditions. J. Geophys. Res., 100(B7):13045–13064, 1995. J. H. Dieterich and B. D. Kilgore. Direct observation of frictional contacts: New insights for sliding memory effects. Pure Appl. Geophys., 143:283–302, 1994. James H Dieterich and Brian D Kilgore. Imaging surface contacts: power law contact distributions and contact stresses in quartz, calcite, glass and acrylic plastic. Tectonophysics, 256(1-4):219–239, 1996. M. Nakatani. Conceptual and physical clarification of rate and state friction: Frictional sliding as a ther- mally activated rheology. J. Geophys. Res., 106(B7):13347–13380, 2001. BM Carpenter, MJ Ikari, and C Marone. Laboratory observations of time-dependent frictional strength- ening and stress relaxation in natural and synthetic fault gouges. J. Geophys. Res., 121(2):1183–1201, 2016. 116 Marco M Scuderi, Cristiano Collettini, C Viti, Elisa Tinti, and Chris Marone. Evolution of shear fabric in granular fault gouge from stable sliding to stick slip and implications for fault slip mode. Geology, 45(8): 731–734, 2017. Frederick M Chester. Dynamic recrystallization in semi-brittle faults. Journal of structural geology, 11(7): 847–858, 1989. F. M. Chester. Effects of temperature on friction: Constitutive equations and experiments with fault gouge. J. Geophys. Res., 99(B4):7247–7261, 1994. N. H. Sleep. Ductile creep, compaction, and rate and state dependent friction within major fault zones. J. Geophys. Res., 100(B7):13065–13080, 1995. Norman H Sleep. Real contacts and evolution laws for rate and state friction. Geochemistry, Geophysics, Geosystems, 7(8), 2006. S Barbot. Modulation of fault strength during the seismic cycle by grain-size evolution around contact junctions. Tectonophysics, 765:129–145, 2019c. H. Perfettini and J.-P. Avouac. Postseismic relaxation driven by brittle creep: a possible mechanism to rec- oncile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan. J. Geophys. Res., 109(B02304), 2004. S. Barbot, Y. Fialko, and Y. Bock. Postseismic Deformation due to the Mw 6.0 2004 Parkfield Earthquake: Stress-Driven Creep on a Fault with Spatially Variable Rate-and-State Friction Parameters. J. Geophys. Res., 114(B07405), 2009. doi: 10.1029/2008JB005748. L. Bruhat, S. Barbot, and J. P. Avouac. Evidence for postseismic deformation of the lower crust following the 2004 Mw6.0 Parkfield earthquake. J. Geophys. Res., 116:10, 2011. B. Rousset, S. Barbot, J. P. Avouac, and Y.-J. Hsu. Postseismic Deformation Following the 1999 Chi-Chi Earthquake, Taiwan: Implication for Lower-Crust Rheology. J. Geophys. Res., 117(B12405):16, 2012. J. C. Rollins, S. Barbot, and J.-P. Avouac. Mechanisms of postseismic deformation following the 2010 el mayor-cucapah earthquake. Pure App. Geophys., page 54, 2015. Lujia Feng, S. Barbot, E. M Hill, I. Hermawan, P. Banerjee, and D. H Natawidjaja. Footprints of past earthquakes revealed in the afterslip of the 2010 mw 7.8 mentawai tsunami earthquake. Geophys. Res. Lett., 43(18):9518–9526, 2016. S. Masuti, S. Barbot, S. Karato, L. Feng, and P. Banerjee. Upper mantle water stratification inferred from the 2012 Mw 8.6 Indian Ocean earthquake. Nature, 538:373–377, 2016. doi: 10.1038/nature19783. Semechah KY Lui and Nadia Lapusta. Repeating microearthquake sequences interact predominantly through postseismic slip. Nature communications, 7:13020, 2016. Rino Salman, Emma M Hill, Lujia Feng, Eric O Lindsey, Deepa Mele Veedu, Sylvain Barbot, Paramesh Banerjee, Iwan Hermawan, and Danny H Natawidjaja. Piecemeal rupture of the mentawai patch, suma- tra: The 2008 mw 7.2 north pagai earthquake sequence. J. Geophys. Res., 2017. J. Muto, J. D. P. Moore, S. Barbot, T. Iinuma, Y. Ohta, S. Horiuchi, and H. Iwamori. Persistent deep afterslip driven by nonlinear transient mantle flow after the 2011 Tohoku earthquake. Science Advances, 5(9): eaaw1164, 2019. 117 Ryoichiro Agata, Sylvain D Barbot, Kohei Fujita, Mamoru Hyodo, Takeshi Iinuma, Ryoko Nakata, Tsuyoshi Ichimura, and Takane Hori. Rapid mantle flow with power-law creep explains deformation after the 2011 tohoku mega-quake. Nature communications, 10(1):1385, 2019. M. H. Linker and J. H. Dieterich. Effects of variable normal stress on rock friction: Observations and constitutive relations. J. Geophys. Res., 97:4923–4940, 1992. J. R. Rice and Y. Ben-Zion. Slip complexity in earthquake fault models. Proc. Nat. Natl. Acad. Sci., 93(9): 3811–3818, 1996. N. Kato and T. E. Tullis. A composite rate- and state-dependent law for rock friction. Geophys. Res. Lett., 28:1103–1106, 2001. A. Bizzarri. On the deterministic description of earthquakes. Rev. Geophys., 49(RG3002): 10.1029/2011RG000356, 2011. K Nagata, M Nakatani, and S Yoshida. A revised rate-and state-dependent friction law obtained by con- straining constitutive and evolution laws separately with laboratory data. J. Geophys. Res., 117(B2), 2012. Peter J Hudleston and Labao Lan. Rheological information from geological structures. In Mechanics Prob- lems in Geodynamics Part I, pages 605–620. Springer, 1995. Labao Lan and Peter J Hudleston. The effects of rheology on the strain distribution in single layer buckle folds. Journal of Structural Geology, 17(5):727–738, 1995. Lukas Fuchs, Hemin Koyi, and Harro Schmeling. Numerical modeling of the effect of composite rheology on internal deformation in down-built diapirs. Tectonophysics, 646:79–95, 2015. Naoyuki Kato. Repeating slip events at a circular asperity: Numerical simulation with a rate-and state- dependent friction law. Bull. Earthq. Res. Inst. Univ. Tokyo, 78:151–166, 2003. T. Chen and N. Lapusta. Scaling of small repeating earthquakes explained by interaction of seismic and aseismic slip in a rate and state fault model. J. Geophys. Res., 114(B01311):12 PP., 2009. Yuta Mitsui and Kazuro Hirahara. Fault instability on a finite and planar fault related to early phase of nucleation. J. Geophys. Res., 116(B6), 2011. Naoyuki Kato. Deterministic chaos in a simulated sequence of slip events on a single isolated asperity. Geophys. J. Int., 198(2):727–736, 2014. Naoyuki Kato. Earthquake cycles in a model of interacting fault patches: Complex behavior at transition from seismic to aseismic slip. Bull. Seism. Soc. Am., 106(4):1772–1787, 2016. C. H. Scholz. Earthquakes and friction laws. Nature, 391:37–42, January 1998. Robert WH Butler, Clare E Bond, Mark A Cooper, and Hannah Watkins. Interpreting structural geometry in fold-thrust belts: Why style matters. Journal of Structural Geology, 114:251–273, 2018. Alexei NB Poliakov, Renata Dmowska, and James R Rice. Dynamic shear rupture interactions with fault bends and off-axis secondary faulting. Journal of Geophysical Research: Solid Earth, 107(B11):ESE–6, 2002b. 118 Elizabeth L Templeton, Aurélie Baudet, Harsha S Bhat, Renata Dmowska, James R Rice, Ares J Rosakis, and Carl-Ernst Rousseau. Finite element simulations of dynamic shear rupture experiments and dynamic path selection along kinked and branched faults. Journal of Geophysical Research: Solid Earth, 114(B8), 2009. Nadine McQuarrie. Crustal scale geometry of the zagros fold–thrust belt, iran. Journal of Structural Geology, 26(3):519–535, 2004. Gianluca Grando and Ken McClay. Morphotectonics domains and structural styles in the makran accre- tionary prism, offshore iran. Sedimentary Geology, 196(1-4):157–179, 2007. Xin Wang, Shengji Wei, and Wenbo Wu. Double-ramp on the main himalayan thrust revealed by broadband waveform modeling of the 2015 gorkha earthquake sequence. Earth Planet. Sci. Lett., 473:83–93, 2017b. Xin Wang, Kyle Edward Bradley, Shengji Wei, and Wenbo Wu. Active backstop faults in the mentawai re- gion of sumatra, indonesia, revealed by teleseismic broadband waveform modeling. EarthandPlanetary Science Letters, 483:29–38, 2018. L.F. Yue, J. Suppe, and J.-H. Hung. Structural geology of a classic thrust belt earthquake: the 1999 Chi-Chi earthquake Taiwan (M w 7.6). J. Struct. Geol., 27:2058–2083, 2005. Leonardo Seeber, Christian Mueller, Toshiya Fujiwara, Kohsaku Arai, Wonn Soh, Yusuf S Djajadihardja, and Marie-Helene Cormier. Accretion, mass wasting, and partitioned strain over the 26 dec 2004 mw9. 2 rupture offshore aceh, northern sumatra. Earth and Planetary Science Letters, 263(1-2):16–31, 2007. Futoshi Nanayama, Kenji Satake, Ryuta Furukawa, Koichi Shimokawa, Brian F Atwater, Kiyoyuki Shigeno, and Shigeru Yamaki. Unusually large earthquakes inferred from tsunami deposits along the kuril trench. Nature, 424(6949):660, 2003. Toshiya Fujiwara, Shuichi Kodaira, Tetsuo No, Yuka Kaiho, Narumi Takahashi, and Yoshiyuki Kaneda. The 2011 tohoku-oki earthquake: Displacement reaching the trench axis. Science, 334(6060):1240–1240, 2011. Jascha Polet and Hiroo Kanamori. Shallow subduction zone earthquakes and their tsunamigenic potential. Geophysical Journal International, 142(3):684–702, 2000. Charles J Ammon, Hiroo Kanamori, Thorne Lay, and Aaron A Velasco. The 17 july 2006 java tsunami earthquake. Geophys. Res. Lett., 33(24), 2006b. Kenji Satake. Geological and historical evidence of irregular recurrent earthquakes in Japan. Philosophical TransactionsoftheRoyalSocietyA:Mathematical,PhysicalandEngineeringSciences, 373(2053):20140375, 2015. doi: 10.1098/rsta.2014.0375. Belle Philibosian and Aron J Meltzner. Segmentation and supercycles: A catalog of earthquake rupture patterns from the sumatran sunda megathrust and other well-studied faults worldwide. Quaternary Science Reviews, 241:106390, 2020. doi: 10.1016/j.quascirev.2020.106390. J.-P. Avouac, F. Ayoub, S. Leprince, O. Konca, and D. V. Helmberger. The 2005, Mw 7.6 Kashmir earthquake: Sub-pixel correlation of ASTER images and seismic waveforms analysis. Earth. Plan. Sc. Let., 249(3-4): 514?528, 2006. Q. Shi, S. Barbot, B. Shibazaki, T. Matsuzawa, S. Wei, and P. Tapponnier. Structural control and system- level behavior of the seismic cycle at the Nankai trough. Earth Planets Space, 72(1):1–31, 2020b. doi: 10.1186/s40623-020-1145-0. 119 S. Barbot. Frictional and structural controls of seismic super-cycles at the Japan trench. Earth Planets Space, 72(63), 2020. doi: 10.1186/s40623-020-01185-3. MR Pandey, RP Tandukar, JP Avouac, J Lave, and JP Massot. Interseismic strain accumulation on the himalayan crustal ramp (nepal). Geophysical Research Letters, 22(7):751–754, 1995. Wang-Ping Chen and Peter Molnar. Seismic moments of major earthquakes and the average rate of slip in central asia. J. Geophys. Res., 82(20):2945–2969, 1977. SN Sapkota, L Bollinger, Y Klinger, P Tapponnier, Y Gaudemer, and D Tiwari. Primary surface ruptures of the great himalayan earthquakes in 1934 and 1255. Nature Geoscience, 6(1):71–76, 2013. Roger Bilham. Location and magnitude of the 1833 nepal earthquake and its relation to the rupture zones of contiguous great himalayan earthquakes. Current Science, 69(2):101–128, 1995. Roger Bilham et al. Earthquakes in India and the Himalaya: tectonics, geodesy and history. Annals of Geophysics, 2004. Jean-Philippe Avouac, Lingsen Meng, Shengji Wei, Teng Wang, and Jean-Paul Ampuero. Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake. Nature Geosci., 2015a. John Galetzka, Diego Melgar, Joachim F Genrich, Jianghui Geng, Susan Owen, Eric O Lindsey, Xiaohua Xu, Yehuda Bock, J-P Avouac, Lok Bijaya Adhikari, et al. Slip pulse and resonance of the Kathmandu basin during the 2015 Gorkha earthquake, Nepal. Science, 349(6252):1091–1095, 2015. Eric O Lindsey, Ryo Natsuaki, Xiaohua Xu, Masanobu Shimada, Manabu Hashimoto, Diego Melgar, and David T Sandwell. Line-of-sight displacement from ALOS-2 interferometry: Mw 7.8 Gorkha Earthquake and Mw 7.3 aftershock. Geophys. Res. Lett., 2015. Kang Wang and Yuri Fialko. Slip model of the 2015 Mw 7.8 Gorkha (Nepal) earthquake from inversions of ALOS-2 and GPS data. Geophys. Res. Lett., 2015. Raphaël Grandin, Martin Vallée, Claudio Satriano, Robin Lacassin, Yann Klinger, Martine Simoes, and Lau- rent Bollinger. Rupture process of the mw= 7.9 2015 gorkha earthquake (nepal): Insights into himalayan megathrust segmentation. Geophys. Res. Lett., 42(20):8373–8382, 2015. Suryodoy Ghoshal, Nadine McQuarrie, Delores M Robinson, DP Adhikari, Leah E Morgan, and Todd A Ehlers. Constraining central Himalayan (Nepal) fault geometry through integrated thermochronology and thermokinematic modeling. Tectonics, page e2020TC006399, 2020a. doi: 10.1029/2020TC006399. Jean-Louis Mugnier, François Jouanne, Roshan Bhattarai, Joaquim Cortes-Aranda, Ananta Gajurel, Pascale Leturmy, Xavier Robert, Bishal Upreti, and Riccardo Vassallo. Segmentation of the himalayan megath- rust around the gorkha earthquake (25 april 2015) in nepal. JournalofAsianEarthSciences, 141:236–252, 2017. Luca Dal Zilio, Romain Jolivet, and Ylona van Dinther. Segmentation of the Main Himalayan Thrust illuminated by Bayesian inference of interseismic coupling. Geophys. Res. Lett., 47(4):e2019GL086424, 2020. doi: 10.1029/2019GL086424. Zacharie Duputel, Jérôme Vergne, Luis Rivera, Gérard Wittlinger, Véronique Farra, and György Hetényi. The 2015 Gorkha earthquake: a large event illuminating the Main Himalayan Thrust fault. Geophys. Res. Lett., 43(6):2517–2525, 2016. 120 Sohom Ray and Robert C Viesca. Earthquake nucleation on faults with heterogeneous frictional properties, normal stress. J. Geophys. Res., 122(10):8214–8240, 2017. C Cattania. Complex earthquake sequences on simple faults. Geophys. Res. Lett., 46(17-18):10384–10393, 2019. A-A Gabriel, J-P Ampuero, Luis A Dalguer, and Paul Martin Mai. The transition of dynamic rupture styles in elastic media under velocity-weakening friction. J. Geophys. Res., 117(B9), 2012b. Sharadha Sathiakumar, Sylvain Denis Barbot, and Piyush Agram. Extending resolution of fault slip with geodetic networks through optimal network design. J. Geophys. Res., 122(12):10–538, 2017. T. Ader, J.-P. Avouac, J. Liu-Zeng, H. Lyon-Caen, L. Bollinger, J. Galetzka, J. Genrich, M. Thomas, K. Cha- nard, S. N. Sapkota, S. Rajaure, P. Shrestha, L. Ding, and M. Flouzat. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard. J. Geophys. Res., 117(B04403):16, 2012. doi: 10.1029/2011JB009071. VL Stevens and JP Avouac. Interseismic coupling on the main Himalayan thrust. Geophys. Res. Lett., 42 (14):5828–5837, 2015. David Mencin, Rebecca Bendick, Bishal Nath Upreti, Danda Pani Adhikari, Ananta Prasad Gajurel, Roshan Raj Bhattarai, Hari Ram Shrestha, Tara Nidhi Bhattarai, Niraj Manandhar, John Galetzka, et al. Himalayan strain reservoir inferred from limited afterslip following the Gorkha earthquake. Nature Geoscience, 9(7):533–537, 2016. Bin Zhao, Roland Bürgmann, Dongzhen Wang, Kai Tan, Ruilin Du, and Rui Zhang. Dominant controls of downdip afterslip and viscous relaxation on the postseismic displacements following the mw7. 9 gorkha, nepal, earthquake. J. Geophys. Res., 122(10):8376–8401, 2017. Kelin X Whipple, Manoochehr Shirzaei, Kip V Hodges, and J Ramon Arrowsmith. Active shortening within the himalayan orogenic wedge implied by the 2015 gorkha earthquake. NatureGeoscience, 9(9):711, 2016. L Wang and S Barbot. Excitation of San Andreas tremors by thermal instabilities below the seismogenic zone. Science Advances, (eabb2057), 2020. doi: 10.1126/sciadv.abb2057. Lu Li, Dongdong Yao, Xiaofeng Meng, Zhigang Peng, and Baoshan Wang. Increasing seismicity in Southern Tibet following the 2015 Mw 7.8 Gorkha, Nepal earthquake. Tectonophysics, 714:62–70, 2017. doi: 10. 1016/j.tecto.2016.08.008. Benchun Duan and David D Oglesby. Nonuniform prestress from prior earthquakes and the effect on dynamics of branched fault systems. J. Geophys. Res., 112(B5), 2007. doi: 10.1029/2006JB004443. Pierre Romanet, Dye SK Sato, and Ryosuke Ando. Curvature, a mechanical link between the geometrical complexities of a fault: application to bends, kinks and rough faults. Geophys. J. Int., 2020. P Thakur, Y Huang, and Y Kaneko. Effects of low-velocity fault damage zones on long-term earth- quake behaviors on mature strike-slip faults. J. Geophys. Res., page e2020JB019587, 2020. doi: 10.1029/2020JB019587. C Kyriakopoulos, AV Newman, AM Thomas, M Moore-Driskell, and GT Farmer. A new seismically con- strained subduction interface model for Central America. J. Geophys. Res., 120(8):5535–5548, 2015. doi: 10.1002/2014JB011859. 121 Mark R Handy, Joerg Giese, Stefan M Schmid, Jan Pleuger, Wim Spakman, Kujtim Onuzi, and Kamil Us- taszewski. Coupled crust-mantle response to slab tearing, bending, and rollback along the Dinaride- Hellenide orogen. Tectonics, 38(8):2803–2828, 2019. doi: 10.1029/2019TC005524. C McA Powell and PJ Conaghan. Plate tectonics and the himalayas. Earth and Planetary Science Letters, 20(1):1–12, 1973. Peter Molnar and Paul Tapponnier. Cenozoic tectonics of asia: Effects of a continental collision: Features of recent continental tectonics in asia can be interpreted as results of the india-eurasia collision. Science, 189(4201):419–426, 1975. Steven G Wesnousky, Senthil Kumar, R Mohindra, and VC Thakur. Uplift and convergence along the himalayan frontal thrust of india. Tectonics, 18(6):967–976, 1999. An Yin and T Mark Harrison. Geologic evolution of the himalayan-tibetan orogen. Annualreviewofearth and planetary sciences, 28(1):211–280, 2000. Jérôme Lavé and Jean-Philippe Avouac. Fluvial incision and tectonic uplift across the himalayas of central nepal. Journal of Geophysical Research: Solid Earth, 106(B11):26561–26591, 2001. Roger Bilham. Himalayan earthquakes: a review of historical seismicity and early 21st century slip poten- tial. Geological Society, London, Special Publications, 483(1):423–482, 2019. Harsh Gupta and VK Gahalaut. Seismotectonics and large earthquake generation in the himalayan region. Gondwana research, 25(1):204–213, 2014. Steven G Wesnousky. Great pending himalaya earthquakes. Seismological Research Letters, 91(6):3334– 3342, 2020. Simon Lamb. The relation between short-and long-term deformation in actively deforming plate boundary zones. Philosophical Transactions of the Royal Society A, 379(2193):20190414, 2021. Eric O Lindsey, Rafael Almeida, Rishav Mallick, Judith Hubbard, Kyle Bradley, Louisa LH Tsang, Yixiang Liu, Roland Burgmann, and Emma M Hill. Structural control on downdip locking extent of the himalayan megathrust. Journal of Geophysical Research: Solid Earth, 123(6):5265–5278, 2018. Shuiping Li, Qi Wang, Shaomin Yang, Xuejun Qiao, Zhaosheng Nie, Rong Zou, Kaihua Ding, Ping He, and Gang Chen. Geodetic imaging mega-thrust coupling beneath the himalaya. Tectonophysics, 747: 225–238, 2018. Laurent Bollinger, Soma Nath Sapkota, P Tapponnier, Yann Klinger, Magali Rizza, Jerome Van der Woerd, DR Tiwari, R Pandey, Adnand Bitri, and S Bes de Berc. Estimating the return times of great himalayan earthquakes in eastern nepal: Evidence from the patu and bardibas strands of the main frontal thrust. Journal of Geophysical Research: Solid Earth, 119(9):7123–7163, 2014. Mari Hamahashi, Judith A Hubbard, Rafael V Almeida, Samuel H Haines, Lewis A Owen, Sanjita Mishra, and Soma Nath Sapkota. Fluvial sedimentary response to late quaternary climate and tectonics at the himalayan frontal thrust, central nepal. Geochemistry, Geophysics, Geosystems, 23(9):e2022GC010366, 2022. Carole Lemonnier, Guy Marquis, Frédéric Perrier, Jean-Philippe Avouac, Gyani Chitrakar, Basantha Kafle, Som Sapkota, Umesh Gautam, Dilliram Tiwari, and Maksim Bano. Electrical structure of the himalaya of central nepal: High conductivity around the mid-crustal ramp along the mht. Geophysical Research Letters, 26(21):3261–3264, 1999. 122 Jean-Philippe Avouac, Lingsen Meng, Shengji Wei, Teng Wang, and Jean-Paul Ampuero. Lower edge of locked main himalayan thrust unzipped by the 2015 gorkha earthquake. NatureGeoscience, 8(9):708–711, 2015b. MR Pandey and Peter Molnar. The distribution of intensity of the bihar-nepal earthquake of 15 january 1934 and bounds on the extent of the rupture zone. JournalofNepalGeologicalSociety, 5(1):22–44, 1988. Sharadha Sathiakumar and Sylvain Barbot. The stop-start control of seismicity by fault bends along the main himalayan thrust. Communications Earth & Environment, 2(1):1–11, 2021. Douglas W Burbank, John Leland, Eric Fielding, Robert S Anderson, Nicholas Brozovic, Mary R Reid, and Christopher Duncan. Bedrock incision, rock uplift and threshold hillslopes in the northwestern himalayas. nature, 379(6565):505–510, 1996. Vincent Godard, Didier L Bourlès, Françoise Spinabella, Douglas W Burbank, Bodo Bookhagen, G Burch Fisher, Adrien Moulin, and Laëtitia Léanni. Dominance of tectonics over climate in himalayan denuda- tion. Geology, 42(3):243–246, 2014. Rajeeb Lochan Mishra, I Singh, A Pandey, PS Rao, HK Sahoo, and R Jayangondaperumal. Paleoseismic evidence of a giant medieval earthquake in the eastern himalaya. Geophysical Research Letters, 43(11): 5707–5715, 2016. CP Rajendran and Kusala Rajendran. Paleoseismic context of the 1950 earthquake: Implications for seismic gaps in the eastern himalaya. Physics and Chemistry of the Earth, Parts A/B/C, 124:103055, 2021. Suryodoy Ghoshal, Nadine McQuarrie, Delores M Robinson, DP Adhikari, Leah E Morgan, and Todd A Ehlers. Constraining central himalayan (nepal) fault geometry through integrated thermochronology and thermokinematic modeling. Tectonics, 39(9):e2020TC006399, 2020b. Paul Tapponnier, Xu Zhiqin, Francoise Roger, Bertrand Meyer, Nicolas Arnaud, Gérard Wittlinger, and Yang Jingsui. Oblique stepwise rise and growth of the tibet plateau. science, 294(5547):1671–1677, 2001. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. NumericalRecipesinC:TheArtofScientific Computing. 2nd ed., 994 pp., Cambridge Univ. Press, New York, 1992. J. R. Rice, N. Lapusta, and K. Ranjith. Rate and state dependent friction and the stability of sliding between elastically deformable solids. J. Mech. Phys. Solids, 49:1865–1898, 2001. N. Lapusta, J. R. Rice, Y. BenZion, and G. Zheng. Elastodynamics analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction. J. Geophys. Res., 105 (B10):23765–23789, October 2000b. A. M. Rubin and J.-P. Ampuero. Earthquake nucleation on (aging) rate and state faults. J. Geophys. Res., 110(B11312):24 PP., 2005. Robert C Viesca. Stable and unstable development of an interfacial sliding instability. Physical Review E, 93(6):060202, 2016. Brittany A Erickson, Junle Jiang, Michael Barall, Nadia Lapusta, Eric M Dunham, Ruth Harris, Lauren S Abrahams, Kali L Allison, Jean-Paul Ampuero, Sylvain Barbot, et al. The community code verification exercise for simulating sequences of earthquakes and aseismic slip (SEAS).SeismologicalResearchLetters, 91(2A):874–890, 2020. 123 AppendixA SupplementaryInformationforchapter2 ∗ ∗ Earthquake cycles in fault-bend folds 124 A.1 Numericalmethod The main constituents of the model to simulate fault dynamics and the particle flow across the axial surface are discussed in the main text. Here, we discuss the numerical method in more detail. The governing equations described below are solved numerically using the Runge-Kutta solver with adaptive time steps and four/fifth order accuracy [Press et al., 1992]. The multiplicative form of rate-and-state friction of equation (2.15) allows truly stationary contact, which causes numerical issues when the dynamic variable is the logarithm of the velocity. Indeed, for this choice of parameterization, the dynamic range of logarithm of velocity is unbounded as velocity vanishes. To avoid this issue, we employ the additive form of rate- and-state friction [Dieterich, 1979, Ruina, 1983]. The additive and multiplicative forms of rate-and-state friction produce similar results on strong faults, i.e., when the stress drops are much smaller than the long- term fault strength [Barbot, 2019c,b]. We also ignore the evolution of temperature with shear heating and therefore consider the constitutive laws in isothermal conditions τ = µ 0 +alog V V 0 +blog ΨV 0 L ¯σ, (A.1) with the state evolution law simplifying to the aging law ˙ Ψ=1 − VΨ L , (A.2) where the parameters are described in the main text. The rate of change of shear traction on the fault is due to the rate-and-state dependence of friction and dynamic changes in normal stress due to elastic interactions on a non-planar fault ˙ τ = ˙ µ ¯σ +µ ˙ ¯σ, (A.3) where µ is the friction coefficient and ¯σ is the effective normal stress, accounting for a constant pore pressure. The change of traction due to fault slip and flow across the axial surface are obtained using the integral method ˙ τ (x)= Z ∂Ω K(x;y) (V(y)− V L (y)) dy− G 2V s ˙ V(x), (A.4) where K a traction kernel that includes fault-to-fault, fault-to-axial surface, axial surface-to-fault, and axial surface-to-axial surface stress interactions. The traction kernel corresponds to the solution of the distribution of stress in the surrounding medium due to a point double-couple source. Solutions to these problems are provided by various authors and summarized by Okada (1992). The last term in equation (A.4) denotes the radiation damping that approximates the inertial effects at seismic velocities [Rice et al., 2001], with the shear wave velocity V s . The change of normal stress due to fault slip and flow across the axial surface are obtained using another integral ˙ ¯σ = Z ∂Ω J(x;y) (V(y)− V L (y)) dy, (A.5) where J is a traction kernel that includes fault-to-fault and axial surface-to-fault interactions. With our assumptions, the rheology of axial surfaces does not depend on pressure or normal stress. Writing ˙ q = Z ∂Ω J(x;y) (V(y)− V L (y)) dy (A.6) 125 for convenience, and combining equations (A.1), (A.3), (A.5), and (A.6), the acceleration of faults slip be- comes ˙ V V = ˙ τ − b¯σ ˙ Ψ/Ψ − µ ˙ ¯σ a¯σ +GV/(2V s ) . (A.7) For the axial surface, the particle velocity is directly given by equation 2.17. We calculate the evolution of frictional parameters, stress, velocity on the fault segments, and velocity across the axial surface numerically using the boundary integral method. For convenience, we treat faults and axial surfaces similarly in numerical calculations, i.e., by varying the take-off angle (faults have zero take-off angle, β =0 ◦ , whereas axial surface have non-zero take-off angle β ̸=0 ◦ ). This allows us to treat the kinematics, stress interactions, and rheological aspects of both types of structures identically in the numerical model. We approximate the integrals in equations (A.4) and (A.5) by finite sums on piecewise constant fault elements and axial surface elements. The fault and axial surface segments are discretized into small 20 m patches in all our simulations to attain numerical convergence and avoid numerical noise due to unresolved models. The fault discretization has to be sufficiently small in order to resolve both the nucleation size h ∗ and the cohesive or breakdown zone L c [Lapusta et al., 2000b, Rubin and Ampuero, 2005, Rubin, 2008, Viesca, 2016], given by h ∗ ≈ GL (b− a)¯σ , (A.8) and L c ≈ R u h ∗ . (A.9) We choose 20 m so that our patch size is less than the minimum value ofh ∗ /20 andL c /5. We verify that this discretization converges by running simulations with varying patch sizes, both smaller and larger than 20 m. We find that patch sizes larger than 100 m do not converge (Figure A.1). Numerical convergence is obtained when the fault behavior is consistent across the seismic cycle for decreasing sampling size. Our investigation of the effect of sampling reveals small differences among simulations due to round-off errors. These differences are minor and do not affect our conclusion. At each time step, we calculate the rate of change of shear traction ˙ τ , rate of change of state variable ˙ Ψ, and slip acceleration ˙ V . Our simulations are quasi-dynamic for we do not resolve the wave-mediated stress transfer. 126 Figure A.1: Convergence test performed on a synclinal fault-bend fold, with a ramp dip of 40 ◦ and lower segment dip of 0 ◦ (change in dipϕ =40), θ = 40 ◦ , and R v = 1.6 for varying R u values and patch sizes. Large patch sizes do not resolve the model well and numerical convergence is achieved for smaller patch sizes. 127 AppendixB SupplementaryInformationforchapter3 ∗ ∗ The stop-start control of seismicity by fault bends along the Main Himalayan Thrust 128 B.1 Methods:Exploringthefrictionalregimewithrate-andstate-dependent friction B.1.1 Constitutiveframework We use the rate- and state-dependent friction framework, first introduced by Dieterich, Ruina, and Rice [Di- eterich, 1979, Ruina, 1983, Rice and Ruina, 1983], supported by extensive experimental data, to model fault dynamics. Seismic cycle simulations using this framework is capable of explaining a gamut of rupture modes and styles observed in nature, which includes fault slip occurring at all stages of a seismic cycle like slow earthquakes, slow-slip events, and quasi-static aseismic processes. We follow the physics-based formulation of [Barbot, 2019c], whereby the fault strength depends on the area of contacts that bear the normal and shear loads. The real area of contact per unit fault area is given by A= µ 0 ¯σ χ ψV 0 L b µ 0 , (B.1) whereµ 0 is the coefficient of static friction, ¯σ is the effective normal stress, χ is the indentation hardness,V 0 is a reference velocity,ψ is a state variable associated with contact aging,L is the characteristic weakening distance, andb≪ 1 is a power exponent. The constitutive law for the sliding velocity is given by V =V 0 τ Aχ µ 0 a exp − Q R 1 T − 1 T 0 , (B.2) where a ≪ 1 is a power exponent, Q is an activation energy, R is the universal gas constant, and T is the absolute temperature. We obtain the constitutive equations for rate- and state-dependent friction by combining (B.1) and (B.2) as V =V 0 τ µ 0 ¯σ µ 0 a ψV 0 L − b a exp − Q R 1 T − 1 T 0 , (B.3) The state variable captures the evolution of the real area of contact with sliding history and follows the evolution law ˙ ψ =exp − H R 1 T − 1 T 0 − ψV L , (B.4) where the first and second terms on the right-hand side correspond to healing and weakening of the fault surface, respectively. In this study, we ignore changes of temperature induced by shear heating or other reactions and limit ourselves to isothermal conditions withT =T 0 . The model is more adequate for seismic cycles simulations than previous formulations [Ruina, 1983] because it avoids issues with vanishing velocities. We also ignore the changes in normal stress along the fault even though changes in normal stress at fault bends can make the contact friction-less or in unbounded compression after only a few seismic cycles, making the numerical simulation intractable. In reality, these stresses are relaxed by folding and other visco-elasto-plastic off-fault processes, but these are challenging to model because experimental data on low temperature plastic deformation mechanisms are not available. In our simulations, we assume that the stress perturbations at fault bends are effectively counteracted by other processes that relieve the bends of these perturbations instantaneously, maintaining the effective normal stress constant. 129 B.1.2 Governingequations We simulate fault slip and rupture dynamics assuming isothermal conditions in a quasi-dynamic set-up. Our model does not include the wave-mediated stress transfers and we employ the Runge-Kutta solver with adaptive time steps and four/fifth order accuracy to simulate earthquake cycles numerically. We track the changes in shear traction due to fault slip using the integral method ˙ τ (x)= Z ∂Ω K(x;y) (V(y)− V L (y)) dy− G 2V s ˙ V(x), (B.5) whereK is the traction kernel or Green’s functions tensor that describes the fault stress interactions,V L is the loading rate on the fault due to the far-field tectonic loading, V denotes the instantaneous velocity vector, and V s is the shear wave velocity. The traction kernel in equation B.5 describes the distribution of stress in the surrounding medium due to a point double-couple source. Radiation damping that is used to approximate the inertial effects due to wave-mediated stress transfers is denoted by the last term on the right hand side of equation B.5. We numerically solve the kinematics, stress interactions, velocity, and constitutive behavior using the boundary integral method, which was partially benchmarked against other methods [Erickson et al., 2020]. B.1.3 Controlsonfaultdynamicsfromdimensionalanalysis Simulations with rate- and state-dependent friction often display complexity depending on the distribution of physical properties [Wu and Chen, 2014, Cattania and Segall, 2018]. Different physical parameters con- trol various different aspects of the seismic cycle, for example, stress drop, recurrence times, and whether the rupture is slow or fast. Dimensional analysis of the physical parameters that control fault dynamics and rupture styles is useful to reduce the vast parameter space involved and identify the combinations of parameters that control rupture dynamics and rupture style. We use the two non-dimensional parameters that exhibit the most variability in nature [Barbot, 2019c] to explore the frictional regime that controls rupture styles. The Dieterich-Ruina-Rice numberR u is de- fined for a fault system with a single velocity-weakening asperity as R u = (b− a)¯σ G W L , (B.6) where G is the rigidity, and W is the size of the velocity-weakening region, i.e., the width of the seis- mogenic zone. TheR u number represents a ratio of the seismogenic zone size to a characteristic rupture nucleation dimension. Stable fault slip in velocity-weakening regions are associated withR u < 1, owing to its large nucleation size compared to the fault dimension. The value of R u increases for increasing weakening behavior with unstable slip represented byR u ≫ 1 values. Physical properties of the velocity- weakening region that fall in increasing values ofR u exhibit cycles of characteristic and periodic ruptures, transitioning to more chaotic sequences [Barbot, 2019b, Cattania, 2019]. For example, higher values ofR u give rise to chaotic cycles with full ruptures followed by partial ruptures of varying sizes. The second non-dimensional number considered is R b = b− a b , (B.7) which describes the ratio of frictional properties that control the static and dynamic stress drops. The ratio of these fundamental source properties control the style of rupture and the relative fraction of the rupture energy consumed in fracture. Large R b values favor ruptures with strong-weakening behavior 130 Parameter Symbol Value Effective normal stress ¯σ 50 MPa Reference velocity V 0 10 −6 m/s Rigidity G 30 GPa Poisson’s ratio ν 0.25 Shear wave speed V s 3× 10 3 m/s Loading rate V L 20 mm/yr Steady-state rate dependence in velocity-strengthening region a− b +1× 10 −2 Steady-state rate dependence in velocity-weakening region a− b −4 × 10 −3 Characteristic weakening distance L 0.01− 2 m Down-dip width of the velocity-weakening region for model H W 92 km Down-dip width of the velocity-weakening region for model E W 89.58 km Model H Dieterich-Ruina-Rice number R u 0.3− 61.3 Non-dimensional parameter R b 0.2 Model E Dieterich-Ruina-Rice number R u 0.3− 59.7 Non-dimensional parameter R b 0.2 Table B.1: Summaryoffrictionalparametersthatexplorethenon-dimensionalparameterR u at fixed R b = 0.2. Frictional parameters and its values using model E and model H sub-surface fault geom- etry in the simulations presented in Supplementary Figure B.2. The static friction coefficient, reference velocity, shear wave speed, and Poisson’s ratio are the same for all simulations. We vary the character- istic weakening distanceL for simulations that explore a range of non-dimensional parametersR u . The down-dip width of the velocity-weakening region for geometry E and H is fixed for all simulations such that the depth of the transition zone is at∼ 15 km. with proportionally limited fracture energy. LowR b values represent conditions close to velocity neutral representative of the transition zone between velocity-weakening and velocity-strengthening friction. The R u -R b phase space can fill up the spectrum of fault behavior ranging from simple and periodic to deterministic chaos. Large R u values result in complex seismic cycles with full and partial ruptures compared to low values of R u that result in slow-slip and slow earthquakes. Although it is difficult to constrain these values in natural settings, an exploration of this phase space, would help us understand rupture dynamics and behavior. Also, complex fault geometry and shear stress variations due to fault non- planarity can affect rupture styles in more complicated ways than observed for planar faults [Sathiakumar et al., 2020]. We design a careful study that takes into account both fault geometry and the frictional regime in theR u -R b phase space considering that we are interested in the complex rupture behavior of the MHT. 131 0 20 40 60 80 100 120 140 160 35 30 25 20 15 10 5 0 Distance from MFT (km) model E VW VS model H VW VS Transition zone 21 o ~23 o 3 o 7 o 26 o 7 o ~7 o ~16 o ~3 o Depth (km) 26 o Middle ramp Shallow bend Mid-crustal ramp Figure B.1: Two-dimensionalfaultgeometrymodelsEandH. The fault geometry along a represen- tative cross-section perpendicular to the plate boundary that crosses the Gorkha-Pokhara Anticlinorium in the middle of the Gorkha coseismic rupture. The dip angles of the individual fault segments are marked for model E (blue) and model H (red). The transition zone of both models is located between 80 and 100 km from the surface trace of MFT (Main Frontal Thrust). The depth of the transition zone is about 15 km in both models. Both the shallow bend and the mid-crustal ramps have similar dip angles and are present at comparable depths. The main difference in fault structure between the two models is in the middle of the seismogenic zone: model E has a single décollement that connects the shallow bend to the mid-crustal ramp in comparison to the more complex fault structure with a middle ramp in model H. 132 Figure B.2: Earthquake cycle dynamics exploring a large range of R u values. a-n) Simulation re- sults for model H (left column) and model E (right column). The dotted white lines indicate the position of the fault bends and the transition zone that represent areas of high stress concentrations. The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The slip-weakening distanceL ranges between 0.5 m to 0.01 m which results inR u values to range between 1 and 60. The x axis in the plots show time steps rather than time itself, which is non-uniformly spaced such that fast events have a finer resolution than during interseismic periods. As R u increases, slip complexity increases for both geometry models. The frictional parameters are listed in Supplementary Table B.1. 133 Figure B.3: EarthquakecycledynamicsexploringR b values. a-j) Simulation results for model H (left column) and model E (right column). The dotted white lines indicate the position of the fault bends and the transition zone that represent areas of high stress concentrations. The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The friction parametersa,b, L are varied to explore theR u − R b phase space. The x-axis in the plots show time steps rather than time itself, which is non-uniformly spaced such that fast events have a finer resolution than during interseismic periods. AsR b increases, pulse-like ruptures and fault bend seismicity increases. The frictional parameters are listed in Supplementary Table B.2. 134 Figure B.4: Pulse-like vs crack-like rupture propagation. With higherR b = 0.9 values, the rupture transitions from crack-like to pulse like with both model E and model H (a and c). In contrast, with lower values of R b = 0.2, the rupture is fully crack-like (b and d). The color indicates the fault slip rate; red represents seismic velocities and dark blue represents fault locking. The x axis in the plots show time steps rather than time itself, which is non-uniformly spaced such that fast events have a finer resolution than during interseismic periods. 135 5000 6000 7000 8000 9000 10 7 Moment per unit length (N) model H model E b a 10000 100 1000 10 Moment per unit length (N) Occurrence (count) 10 8 10 9 10 10 10 11 10 12 16 17 15 18 19 20 21 22 23 24 25 26 10 7 10 8 10 9 10 10 10 11 10 12 Moment per unit length (N) 62 68 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 63 64 65 66 67 69 70 71 71 72 73 10 7 10 8 10 9 10 10 10 11 10 12 Time (years) c model E model H Power-law distribution Figure B.5: Statisticalpropertiesofsimulatedeventsusingthetwoend-membergeometrymod- els. Simulation results computed for 5,000 years for a) model E and b) model H, showing giant surface ruptures (red stars), partial ruptures of varying sizes (purple stars), Gorkha-like events (blue stars), and aftershocks (black stars) with varying magnitudes (cumulative force in 2-D). c) Cumulative histogram of all events for a catalog spanning 45,000 years for both geometries. Model H approximately follows a power- law distribution indicated by the dashed black line. The Pareto distribution (dashed line) is for a power exponentα =0.3 for moment per unit length above5× 10 7 N. 136 Parameter Symbol Value Reference velocity V 0 10 −6 m/s Rigidity G 30 GPa Poisson’s ratio ν 0.25 Shear wave speed V s 3× 10 3 m/s Loading rate V L 20 mm/yr Steady-state rate dependence in velocity-strengthening region a− b +1× 10 −2 Down-dip width of the velocity-weakening region for model H W 92 km Down-dip width of the velocity-weakening region for model E W 89.58 km Simulations A and F Effective normal stress ¯σ 50 MPa Characteristic weakening distance L 0.01 m Steady-state rate dependence in velocity-weakening region forR b =0.2 a− b −4 × 10 −3 Simulations B and G Effective normal stress ¯σ 50 MPa Characteristic weakening distance L 0.04 m Steady-state rate dependence in velocity-weakening region forR b =0.8 a− b −1.76 × 10 −2 Simulations C and H Effective normal stress ¯σ 100 MPa Characteristic weakening distance L 0.04 m Steady-state rate dependence in velocity-weakening region forR b =0.8 a− b −1.76 × 10 −2 Simulations D and I Effective normal stress ¯σ 50 MPa Characteristic weakening distance L 0.05 m Steady-state rate dependence in velocity-weakening region forR b =0.9 a− b −1.98 × 10 −2 Simulations E and J Effective normal stress ¯σ 100 MPa Characteristic weakening distance L 0.05 m Steady-state rate dependence in velocity-weakening region forR b =0.9 a− b −1.98 × 10 −2 Table B.2: Summaryoffrictionalparametersthatexplorethenon-dimensionalparameterR b at fixed R u . Frictional parameters and its values using model E and model H sub-surface fault geometry in the simulations presented in Supplementary Figure B.3. The static friction coefficient, reference velocity, shear wave speed, and Poisson’s ratio are the same for all simulations. We vary the effective normal stress, frictional parametersa andb, and characteristic weakening distanceL for simulations that explore a range of non-dimensional parametersR u andR b . The down-dip width of the velocity weakening region for models E and H is fixed for all simulations, such that the depth of the transition zone is at ∼ 15 km. 137 AppendixC SupplementaryInformationforchapter4 ∗ ∗ Reconciling geologic and geodetic shortening rates with active folding in the Himalayan orogenic wedge 138 Parameter Symbol Value Effective normal stress ¯σ 125 MPa Reference velocity V 0 10 −6 m/s Rigidity G 30 GPa Poisson’s ratio ν 0.25 Shear wave speed V s 3× 10 3 m/s Loading rate V L 17 mm/yr Steady-state rate dependence in velocity-strengthening region a− b +1× 10 −2 Steady-state rate dependence in velocity-weakening region a− b −4 × 10 −3 Steady-state rate dependence in axial surfaces a− b +5× 10 −2 Characteristic weakening distance L 2.5 cm Down-dip width of the velocity-weakening region for model H W 92 km Dieterich-Ruina-Rice number R u 61.33 Non-dimensional parameter R b 0.2857 Table C.1: Summary of frictional parameters for Central Nepal. Frictional parameters and its val- ues used in simulations CN1-CN5. The axial surface patches are strongly rate-strengthening simulating a visco-elasto-plastic behavior. Fault patches are velocity weakening until a depth of∼15 km, and are velocity-strengthening below that depth. Parameter Symbol Value Effective normal stress ¯σ 125 MPa Reference velocity V 0 10 −6 m/s Rigidity G 30 GPa Poisson’s ratio ν 0.25 Shear wave speed V s 3× 10 3 m/s Loading rate V L 17 mm/yr Steady-state rate dependence in velocity-strengthening region a− b +1× 10 −2 Steady-state rate dependence in velocity-weakening region a− b −4 × 10 −3 Steady-state rate dependence in axial surfaces a− b +5× 10 −2 Characteristic weakening distance L 1.87 cm Down-dip width of the velocity-weakening region for model H W 68.9 km Dieterich-Ruina-Rice number R u 61.4 Non-dimensional parameter R b 0.2857 Table C.2:SummaryoffrictionalparametersforEastNepal. Frictional parameters and its values used in simulations EN1-EN3. The axial surface patches are strongly rate-strengthening simulating a visco- elasto-plastic behavior. Fault patches are velocity weakening until a depth of∼15 km, and are velocity- strengthening below that depth. 139 0 20 40 60 80 100 120 140 160 40 30 20 10 0 0 20 40 60 80 100 120 140 160 40 30 20 10 0 0 20 40 60 80 100 120 140 160 40 30 20 10 0 0 20 40 60 80 100 120 140 40 30 20 10 0 0 20 40 60 80 100 120 140 160 40 30 20 10 0 0 20 40 60 80 100 120 140 160 40 30 20 10 0 0 20 40 60 80 100 120 140 40 30 20 10 0 0 20 40 60 80 100 120 140 40 30 20 10 0 Distance from the MFT (km) Distance from the MFT (km) Depth (km) Depth (km) Depth (km) Depth (km) Stage CN1 Stage CN2 Stage CN3 Stage CN4 Stage CN5 Stage EN1 Stage EN2 Stage EN3 Velocity-weakening Velocity-strengthening 7 26 7 26 3 21 7 26 3 21 7 26 3 21 7 26 3 21 7 26 7 26 3 21 7 26 7 26 3 21 7 26 7 26 3 21 7 26 7 26 3 21 a c d e b f g h Figure C.1: Two-dimensional fault geometry models along Central and East Nepal. The fault ge- ometry with the dip angles of each fault segment marked in Central (CN) and East Nepal (EN). At every stage of the structural fault evolution, the fault geometry and frictional properties of Central and East Nepal remain constant, but the orientation of the axial surface tied to each fault bend changes according to the orientation of the thrusted sheets. Velocity weakening (unstable fault slip) regions are marked in red, velocity strengthening regions are marked in black. The transition zone where the behavior changes from velocity strengthening to velocity weakening is at∼15 kms depth. 140 Figure C.2: Earthquake cycle dynamics with slip rate segmentation in Central Nepal. a). The fault geometry of Central Nepal featuring the middle ramp in the seismogenic zone, based on Hubbard et al. [2016] fault model. The fault bend locations and the transition between velocity-strengthening and velocity-weakening zones are highlighted. b-f). Seismic cycle simulations of the various stages in Central Nepal (CN1-CN5) to illustrate the evolution of seismic behavior as the hanging wall deforms. The dashed white lines correspond to the locations of fault bends, and solid line represents the transition zone. The x-axis represents time index, not time itself, that is non-uniformly spaced, so that seismic events have finer resolution compared to aseismic periods and the y axis represents the down-dip fault distance in km. The color indicates the fault slip velocity, red indicates fast seismic events and dark shades of blue indicates fault locking. The frictional parameters are listed in Supplementary Table C.1. Variations in slip rate impact earthquake cycle characteristics like rupture width and magnitude, hypocenter location, and the recurrence times of great and partial ruptures. 141 Figure C.3: Earthquake cycle dynamics with slip rate segmentation in East Nepal. a). The fault geometry of East Nepal without the middle ramp in the seismogenic zone, based on Hubbard et al. [2016] fault model. The fault bend locations and the transition between velocity-strengthening and velocity- weakening zones are highlighted. b-d). Seismic cycle simulations of the various stages in East Nepal (EN1-EN3) to illustrate the evolution of seismic behavior as the hanging wall deforms. The dashed white lines correspond to the locations of fault bends, and solid line represents the transition zone. The x-axis represents time index, not time itself, that is non-uniformly spaced, so that seismic events have finer resolution compared to aseismic periods and the y axis represents the down-dip fault distance in km. The color indicates the fault slip velocity, red indicates fast seismic events and dark shades of blue indicates fault locking. The frictional parameters are listed in Supplementary Table C.2. Variations in slip rate impact earthquake cycle characteristics like rupture width and magnitude, hypocenter location, and the recurrence times of great and partial ruptures. 142
Abstract (if available)
Abstract
Faults and folds are geologic expressions of plate tectonic activity. They accommodate relative motion between two or more tectonic plates, are responsible for catastrophic hazards like earthquakes and tsunamis, and build mountain chains that fabricate Earth’s aerial landscapes. In this thesis, I study the subsurface fault structure complexity and its impact on seismotectonics by integrating geological and geophysical tools. First, I develop physics-based numerical models of seismic cycles that incorporate fault structure complexity. To achieve this, I use fault-bend folding theory, based on observations of structural geology, to define a new geophysical modeling framework that accounts for changes in fault orientations and the accompanying folding processes in the crust. These more realistic rupture models with a single fault-bend fold show that fault geometry and fault friction properties interact to shape the seismic cycle. Next, I study the Main Himalayan Thrust (MHT) fault geometry to understand the role of fault bends on the short-term seismic cycle processes. I use numerical models tuned to geological and geophysical data and observations to explain the MHT fault behavior in the Kathmandu section, the source region of the 2015 Mw 7.8 Gorkha earthquake. Next, I build comprehensive structural models to reconcile contrasting observations of fault slip rates in the Nepal Himalayas by accounting for folding in the orogenic wedge during its structural evolution. Advances in seismic hazard assessment in convergent margins around the world will require integrating a wide range of geologic and geophysical datasets and models with complementary spatial and temporal resolution.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Paleoseismology of blind-thrust faults beneath Los Angeles, California: implications for the potential of system-wide earthquakes to occur in an active fold-and-thrust belt
PDF
Heterogeneity of earthquake stress drops, focal mechanisms and active fault zones
PDF
Towards an understanding of fault-system mechanics: from single earthquakes on isolated faults to millenial-scale collective plate-boundary fault-system behavior
PDF
From laboratory friction to numerical models of fault dynamics
PDF
Volumetric interactions between major ruptures and fault zones illuminated by small earthquake properties
PDF
Microseismicity, fault structure, & the seismic cycle: insights from laboratory stick-slip experiments
PDF
Stress-strain characterization of complex seismic sources
PDF
Analysis of waveform and catalog data of aftershocks for properties of earthquakes and faults
PDF
Dynamic rupture processes and seismic radiation in models of earthquake faults separating similar and dissimilar solids
PDF
Quantifying ground deformation of large magnitude earthquakes using optical imgaging systems
PDF
Multi-scale imaging and monitoring of crustal and fault zone structures in southern California
PDF
Paleoseismologic and slip rate studies of three major faults in southern California: understanding the complex behavior of plate boundary fault systems over millenial timescales
PDF
Integration and validation of deterministic earthquake simulations in probabilistic seismic hazard analysis
PDF
Multi-scale imaging of major fault zones in Southern California
PDF
High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars
PDF
Detection and modeling of slow slip events as creep instabilities beneath major fault zones
PDF
Forecasting directivity in large earthquakes in terms of the conditional hypocenter distribution
PDF
Multi-scale imaging of the fault zone velocity structure: double-difference tomography, inversion of fault zone headwaves, and fault zone sensitivity kernels
PDF
Determining fault zone structure and examining earthquake early warning signals using large datasets of seismograms
PDF
Crustal-scale ductile shear zones in four dimensions
Asset Metadata
Creator
Sathiakumar, Sharadha
(author)
Core Title
The impact of faulting and folding on earthquake cycles in collisional orogens
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Geological Sciences
Degree Conferral Date
2023-05
Publication Date
02/03/2023
Defense Date
01/04/2023
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
earthquake cycle,fault-bend folds,faulting,folding,Geology,geophysics,Himalayas,main Himalayan thrust fault,mountain-building,Nepal,OAI-PMH Harvest,seismic hazards
Format
theses
(aat)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Barbot, Sylvain (
committee chair
), Dolan, James (
committee member
), Lynett, Patrick (
committee member
), Platt, John (
committee member
)
Creator Email
sathiaku@usc.edu,sharadhask03@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC112724697
Unique identifier
UC112724697
Identifier
etd-Sathiakuma-11464.pdf (filename)
Legacy Identifier
etd-Sathiakuma-11464
Document Type
Dissertation
Format
theses (aat)
Rights
Sathiakumar, Sharadha
Internet Media Type
application/pdf
Type
texts
Source
20230206-usctheses-batch-1006
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
earthquake cycle
fault-bend folds
faulting
folding
geophysics
main Himalayan thrust fault
mountain-building
seismic hazards