The Role of Surface Processes in Partitioning of Strain and Uplift , St . Elias Orogen , Southern Alaska

Three-dimensional thermo-mechanical numerical simulations of the ongoing Yakutat–North America collision are used to identify the role of surface processes in triggering localized rapid uplift, exhumation, and strain observed within the St. Elias orogen of southern Alaska. Thermochronological data reveal localized rapid exhumation associated with the Seward-Malaspina and Hubbard Glaciers within a tectonic corner structure where transpressional motion to the south along the Fairweather Fault system transitions to shortening to the north and west within the active fold-and-thrust belt of the St. Elias orogen. The modeled deformation patterns are characteristic of oblique convergence within a tectonic corner, recording the transition from simple shear to contractional strain within a zone spatially consistent with the highest exhumation rates suggesting the corner geometry is the primary control of strain partitioning. The relative roles of surface-related processes versus tectonics-related processes in the development of this pattern of deformation were tested with the inclusion of an erosional surface model. The presence of surface processes enhanced the uplift and development of a localized rapid exhumation. When spatially and temporally erosion models are employed, the location of maxima is shifted in response. This indicates that efficient erosion, and resultant deposition and material advection can influence the localization of strain and uplift.


Introduction
An enduring question of landscape evolution regards the relative role of surface processes (e.g.erosion, landscape development, and mass redistribution) versus tectonic (e.g.kinematics and strain partitioning) in influencing the geomorphic and geologic evolution of orogens.Both processes are involved at the orogen scale, where strongly asymmetric orographic gradients impose a major control on orogen evolution (e.g.Koons, 1990;Koons, 1995;Willett et al., 1993;Whipple, 2009), as well as at catchment scales as a local response to intense erosion (Zeitler et al., 2001;Stewart et al., 2008;Enkelmann et al., 2009;Enkelmann et al., 2011).A positive feedback between erosion and tectonics has been suggested to occur in the eastern syntaxis of the Yakutat -North American collision zone within the St. Elias Mountains in southern Alaska based upon a suite of relatively young zircon fission track ages (< 4 Ma) for sediment derived from the evolving orogen (Enkelmann et al., 2009;Enkelmann et al., 2010).The St. Elias Mountains are relatively young, with subduction-collision ongoing since <10 Ma (e.g.Plafker et al., 1994;Pavlis et al., 2004), thus the suggested feedback between erosion and tectonic processes might represent a nascent stage of the complex coupling of tectonics and erosion.
The goal of this study is to investigate the relative roles of surface and tectonic processes at the Yakutat plate corner as it collides obliquely with the North American margin.Three-dimensional (3D) dynamic thermo-mechanical numerical simulations are employed to model the mechanical and thermal behavior of an oblique collisional orogen.These simulations utilize simple fully coupled mechanical and thermal models of Earth materials to simulate the partitioning strain into fault zones, development characteristic surface uplift and rock exhumation patterns, and the evolution of the thermal structure as the model deforms.The boundary conditions applied to the top surface of the model (depth = 0 km) are varied to test the relative roles of erosion and tectonics.The boundary conditions utilized include a passively deforming boundary (no erosion), a simple erosional landscape model (diffusion-based erosion model), and a more complex, glacially focused erosional landscape model.The results of these simulations are then compared with field observations, seismic data, geodetic velocities, thermochronometry, and structural measurements.
Figure 1.This topographic map of southern Alaska shows the main tectonic and topographic features (Smith & Sandwell, 1997) Note.The map inset shows the approximate location of the model geometry with the blue dashed line.The hatched area depicts the surface area of the Yakutat block.Major faults are shown as red lines (after Plafker et al., 1994).Plate vectors for the Pacific Plate (Pac) and Yakutat Terrane (Yak) are from Fletcher & Freymueller (1999).The black triangles denoting Mt.Logan (5960m; L), Mt.St. Elias (5500m; E), and Mt.Fairweather (4600m; F; Fairweather Range), together with the trace of the coastline, are used for spatial reference for later figures.Gray shaded areas depict current glacial cover (NSDIC, 1999).Labels indicate locations of important features: B: Bagley Glacier; BR: Border Ranges Fault; CS: Chugach-St.Elias Fault; DF: Denali Fault; FQC: Fairweather-Queen Charlotte Fault; KI: Kayak Island Zone; M: Malaspina Glacier; P: Pamplona Zone; RMF: Ragged Mountain Fault; S: Seward Glacier; TF: Transition Fault; Y: Yakutat Bay; YF: Yakutat Foothills.
This study addresses the hypothesis that efficient glacial erosion and subsequent mass redistribution within the active orogen in southern Alaska has altered the spatial patterns of exhumation and strain partitioning (e.g.Pavlis et al., 2004;Berger et al., 2008a,b;Chapman et al., 2008;Enkelmann et al., 2010).This project relies on the development and testing of a series of dynamic numerical models using varied model surface boundary conditions (e.g.erosion and deposition) to determine their overall effects on the long-term partitioning of strain within the evolving orogen through time.Determining what controls exhumation is seldom a possibility in mature orogens (e.g.England & Molnar, 1990).

Recent Tectonics of Southeast Alaska
The St. Elias Mountains, southern Alaska, have been uplifted by the ongoing collision and accretion of the Yakutat terrane, a wedge shaped microplate consisting of a thickened oceanic plateau (~30 km), sliver of allochthonous North American continental basement, and associated sedimentary cover, onto North America (Figure 1; Plafker et al., 1994;Pavlis et al., 2004;Christeson et al, 2010).This orogenic event began with the initiation of the dextral Queen Charlotte-Fairweather fault system at ~30 Ma, which reorganized the western margin of North America by separating the Yakutat terrane and transferring it northward (Plafker et al, 1994).During the late Miocene (~10 Ma), the Yakutat terrane arrived at the Aleutian subduction zone and has since been colliding with and shallowly subducting beneath North America, causing the sedimentary cover of the Yakutat terrane to be reworked in an evolving foreland fold-and-thrust belt (Pavlis et al, 2004).The high topography of southeastern Alaska, including the St. Elias and Alaska Ranges, has been driven by this collision and shallow angle subduction event (Plafker et al., 1994).

Observed Strain Patterns
The spatial patterns of strain observed within southern Alaska are characteristic of those expected from an oblique collision with a syntaxial corner structure.Generally the orogen can be divided into four areas based upon the partitioning of strain and geomorphological expression within the region (please refer to Figure 1 for locations):

Eastern Transpressive Boundary
The southeastern to eastern boundary is characterized by a narrow orogenic highlands, the Fairweather Range, interpreted to be a dextral transpressive boundary related to restraining bend within the Fairweather-Queen Charlotte fault system (Plafker et al., 1978).Geodetic and paleoseismicity estimates along strike motion of 45-50 mm/yr (Fletcher & Freymueller, 2003) with some of the convergence partitioned into a narrow thrust belt in the Yakutat foothills (Chapman et al., 2008;Bruhn et al., 2004).The transpressive Fairweather Fault disappears north of Yakutat Bay within the syntaxial corner of the plate boundary (Bruhn et al., 2012) and the plate motion is inferred to be taken up on oblique convergent structures, where there is an abrupt increase in topography to form the Chugach-St.Elias and Wrangell Mountains.This transition in strain regime to oblique convergence causes a dextral shear zone that formed beneath the Bagley glacier to be reactivated as an oblique thrust fault over last 5 to 6 Ma (Bruhn et al., 2012).This newly activated thrust forms a fault-bounded sliver capable of partitioning ~50 km of dextral plate motion and several km of thrust motion over last 20 Ma (Bruhn et al., 2012).

Fold-and-Thrust Belt and Tectonic Backstop
The main portion of the orogen consists of an actively deforming highland where a thick sequence (7-10 km) of pre-to synorogenic marine to fluvial to glaciomarine sediment, underlain by shallowly subducting basement of the Yakutat Terrane, is being actively deformed with a fold-and-thrust belt that youngs to the south toward the present day deformation front in the offshore Pamplona zone (Plafker, 1987;Chapman et al., 2008;Lowe et al., 2008).The Pamplona Zone, and its onshore equivalent the Malaspina Fault, cut across the grain of the orogen and are considered active structures (Estabrook et al., 1992;Chapman et al., 2012).While the Pamplona Zone is considered the active deformation front, it has undergone only a small amount of shortening over the last half million years (e.g.Lowe et al., 2008), suggesting that most of the deformation is accommodated elsewhere in the orogen (Pavlis et al., 2004), likely the onshore fold-and-thrust belt.This has been suggested to be likely due to a structural reorganization in response to intense glacial erosion and mass redistribution (Pavlis et al., 2004;Berger et al., 2008a).Pavlis et al. (2004) interprets the total amount of shortening recorded in the fold-and-thrust belt (~150-200 km) is only ~33% of the expected total shortening since the Miocene and that erosion and along strike removal accounts for the remaining shortening.This observation is consistent with the erosive removal of material above growing duplex and can be readily seen in the thermochronological record, which shows high rates of exhumation within the orogen (e.g.Berger et al., 2008;Enkelmann et al., 2011).
Two other trends are evident from west to east across the orogen.First, there is an observed increase in elevation.This topographic gradient mirrors the angle of subduction of the Yakutat Terrane and Pacific plate, suggesting the high elevations and strain patterns to the east are related to strain transferred from the shallowly subducting Yakutat terrane.Second, there is an increase in metamorphic grade from west to east, culminating in amphibolite grade rocks in the vicinity of Mt.St. Elias (Dusel-Bacon et al., 1995), suggesting an increase in uplift and exhumation over that region.
The large tectonic escarpment formed by the Chugach-St.Elias Fault (CSEF) suggests that it forms the tectonic backstop for the orogenic wedge, however, thermochronological data suggests that the fault has not been active for the last 5 Ma (Berger et al., 2008a).The Bagley-Contact Fault (CF), originally formed within the pre-existing Aleutian subduction complex (Plafker et al., 1994), has been reactivated during the Neogene collision event and shows substantial southside up motion (Berger et al., 2008).The CF has been suggested to form the backstop to the orogen (Bruhn et al., 2012;Chapman et al., 2012).

Transition Fault
The exact nature of the Transition Fault, which places Oligocene Pacific crust against the Mesozoic to Paleocene basement of the Yakutat Terrane (Bruns, 1983), is controversial (e.g.Pavlis et al., 2004;Gulick et al., 2008).The main controversy centers around the nature of motion on the fault, which has been interpreted as either a relict transform boundary (Bruns, 1983) or an active thrust fault that forms a low-angle decollement that places Yakutat Terrane over the Pacific plate with little to no internal deformation (Plafker et al., 1994;Bruhn et al., 2004).Seismicity and bathymetric data suggests that the fault is active in the southeastern segment (e.g.Doser & Lomas, 2000); however, undeformed sedimentary rocks overlap the fault trace (Pavlis et al., 2004;Gulick et al., 2008).Gulick et al. (2008) theorized that the motion on the Transition Fault is the result of stalling of subduction of the Yakutat terrane and translation of motion from the Fairweather Fault.

Western Boundary
The western boundary of the orogen, or western syntaxis, is defined by the Ragged Mountain Fault segment of the CSEF and takes a 90º counterclockwise bend and forms a N-S structure that cuts across the general pattern of the orogen.The orogen in this area is characterized by relatively subdued topography and diffused strain partitioned onto small-to mid-scale (~10 3 m) tectonic structures (Bruhn et al., 2004;Chapman et al., 2008;Li et al., 2010).
The most prominent feature within this area is the Ragged Mountain Fault, an active thrust suggested to record ~100 km indention of the Alaskan margin by the Yakutat Terrane.The western boundary continues offshore as the Kayak Island Zone (Figure 1).

Quaternary Structural Reorganization
Chapman et al (2008) illustrated that the present deformation front within the orogen does not correspond to the finite strain front, suggesting a complex, large-scale structural reorganization of the orogen.This Quaternary structural reorganization is thought to be driven by orogenic topographic growth and mass redistribution linked to glaciation (Chapman et al., 2008;Berger et al., 2008a).The structural reorganization includes development and growth of the syntaxial corner structures, shifting of the deformation front in response to Pleistocene erosion and sedimentation, and retrothrust motion within the fold-and-thrust belt forming a doubly vergent tectonic wedge (Chapman et al., 2008).The progressive structural reorganization in response to climatic triggers and the overall growth of the orogenic wedge demonstrates that the orogen has widened through time, which is consistent with other studies (e.g.Meigs et al., 2008).

Thermochronological Record of Uplift and Exhumation
The dynamic feedback between climate and orogenesis has become a readily accepted theory (Zeitler et al., 2001;2014;Koons et al., 2002;Stewart et al., 2008;Enkelmann et al., 2011).Many studies have interpreted that the partitioning of strain and exhumation has been driven by reorganization of mass within an orogen (e.g.Zeitler et al., 2001;Berger & Spotilla, 2008).Abundant low temperature thermochronologic data has been collected from the St.
Elias orogen in an effort to determine the rates and patterns of uplift and exhumation within the evolving orogen (Berger et al., 2008a;Enkelmann et al., 2009;McAleer et al., 2009).These data indicate a strong coupling exists between glacial episodes and denudation and deformation.
2.3.1 Exhumation within the Fold-and-Thrust Belt Berger et al (2008a;2008b) collected low-temperature thermochronologic data from bedrock within the thin-skinned fold-thrust belt to the west of Mt.St. Elias and Mt.Logan.This dataset shows that there is a zone of high exhumation rates on the windward flank of the orogen associated with the inferred glacial equilibrium line altitude (ELA) where glacial discharge and abrasion should be at a maximum resulting in efficient erosion and exhumation (4-5 mm/yr exhumation; Hallet, 1979;Berger & Spotilla, 2008).The spatial coverage of ice has waxed and waned over the last 5 Ma with variations in climate, which has caused the location of the ELA to also shift through time (Pewe, 1975;Lagoe et al., 1993;Meigs & Sauber, 2000;Reece et al., 2011;Berger & Spotilla, 2008).Berger et al. (2008a) interpreted that the variation in ELA has altered the rate and location of glacial erosion and redistributed mass within the orogen triggering a structural reorganization with the main zone of deformation moving toward the indentor.

Exhumation within the Eastern Syntaxial Corner
Rapid exhumation is inferred in the area of the Seward and Malaspina Glaciers where short term (~10 4 yr) erosion rates of 10-100 mm/yr within the glacial valleys (Hallet et al., 1996), and long term (~10 6 yr) thermochronological data indicate rapid cooling from crustal levels of >5 km depths in less than 2 Myr (Enkelmann et al., 2009;2010).Enkelmann et al. (2009;2010) provided evidence of the existence of a zone of focused exhumation associated with the Seward-Malaspina glacial system that was interpreted to be the result of coupling between the active tectonic uplift and efficient glacial erosion.

Exhumation within the Western Syntaxial Corner
Arkle et al ( 2013) presented low temperature thermochronological data for the area around the western syntaxis within the Chugach terrane.They found a zone of relatively fast exhumation rates (~0.7 mm/yr) centered on the Contact fault system.The authors suggested the patterns of uplift and exhumation over the last ~5 Ma are consistent with shallow subduction of the Yakutat terrane starting in the Oligocene and underplating Generally, the observed pattern of exhumation determined through detailed study of low-temperature thermochronologic indicates there are high rates of uplift and exhumation within the syntaxial corner structures and the rates decrease into the fold-and-thrust belt.Additionally, thermochronologic studies suggest a strong link between glacial erosion and evolving strain and exhumation patterns.This has been interpreted as the result of a dynamically coupled system where efficient erosion has enhanced the general pattern of strain and uplift (Pavlis et al., 2004;Berger et al., 2008a;Koons et al., 2010).The series of numerical simulations presented here demonstrates that: the strain patterns within Southern Alaska are characteristic of oblique convergence within a syntaxial corner structure; the presence of an erosional surface effects strain patterns, differing in both strain rate and spatial partitioning; and temporally and spatially variable focused erosion can drive a structural reorganization similar to that seen in the St. Elias orogen.

Thermal-Mechanical Modeling Methods
The framework of this study is based upon a series of three-dimensional, dynamically coupled, thermo-mechanical continuum mechanics simulations focusing on the St. Elias Mountain region of southern Alaska that reproduce the temporal and spatial evolution of the entire region as the subduction and collision of the Yakutat terrane progressed (Figure 2; Koons et al., 2010).Simultaneous solution of kinematic, thermal, and mechanical equations is accomplished through mixed implicit-explicit techniques within a modified Lagrangian framework (FLAC 3D ; Cundall & Board, 1988;Koons et al., 2002).The simulations are homogeneous and isotropic, with no explicitly defined pre-existing weaknesses and all deformation within the simulation, including the formation of faults or shear zones, forms dynamically with the accumulation of strain.

Model Geometry and Constitutive Relationships
The three-dimensional solution domain is standardized in size and extends 800 km normal to the boundary (Y-axis) by 700 km parallel to the boundary (X-axis) by 50 km depth (Z-axis) with a 10 km horizontal discretization and 5 km vertical discretization.The location of the solution domain spans the width of the Yakutat terrane such that the boundaries are far from the area of interest (Figures 1 and 2).The simulations use a succession of spatially and temporally varying surface boundary conditions that explore the controls of erosion, mass redistribution, and varying topographic surfaces on the formation and evolution of deformation patterns.
The simulation domain is defined with coupled advective-conductive thermal and thermally dependent rheological constitutive models.The rheological model used is similar to the "jelly sandwich" layered model (Chen and Molnar, 1983;Burov and Watts, 2006) and is a function of temperature, strain rate, and accumulated strain.The mechanical model assigned to the cold upper crust (< 350 ºC) is the standard associated plasticity of the pressure dependent Mohr -Coulomb model, where σ 1 and σ 3 are the major and minor principal stresses derived from the model (MPa), C is the cohesion (MPa) and N = (1asin )/(1  -sin ) where is the internal angle of friction (Table 1).  A series of simulations for this model include a strain-dependent material weakening that simulates the formation of faults and shear zones (e.g.Montesi & Zuber, 2002) by changing the material cohesion and internal angle of friction.As at higher temperatures, the viscous strength of rock depends to the first order on temperature, the hotter lower crust (> 350 ºC) is assigned a Newtonian temperature-dependent viscosity based upon the Frank-Kamenetskii approximation: where η 0 is the reference viscosity (at a reference T), is an exponential factor (i.e.characteristic temperature), and temperature (T) is derived thermal solution.Values for the reference viscosity and characteristic temperature are chosen to yield an effective viscosity of 10 22 to 10 19 Pas (Table 1; Turcotte & Schubert, 2002).While this model is likely much more simple than the crustal rheology, the sharp reduction in strength due to the exponential dependence upon temperature is a robust feature of quartz and feldspar dominated crust (Sibson, 1982).The model basal and surface boundaries are assigned fixed temperatures of 0 C and 750 C, respectively.An average initial geothermal gradient of 15 C/km is applied to the simulation domain.The assigned thermal model is a conductive-advective solution that is mechanically coupled through the kinematic solution, Where is thermal conductivity, is density, t is time, C v is heat capacity, is the velocity field taken from the mechanical solution, and A is the radiogenic heat source (Table 1; Turcotte & Shubert, 2002).

Applied Model Boundary Conditions
The initial velocity conditions are imposed on the base of the model over an area corresponding to the spatial extent of the Yakutat terrane.The basal kinematic conditions are applied so that the surface velocity results are consistent with the observed motions of the Yakutat terrane (Figure 2; ~45 mm/yr; DeMets et al., 1994;Fletcher & Freymueller, 1999;Freymueller et al., 2008;Elliot et al., 2010).
Variable erosional conditions are applied to the model surface (Table 2).The erosion models are based upon 1) the spatial extent of current glaciations with the assumption that glacial erosion maintains a near constant elevation during exhumation (Local glacial; Hallet et al, 1996;Meigs and Sauber, 2000) and 2) diffusion as a function of surface slope, relief, and surface conditions (Regional diffusion; e.g.Burbank and Anderson, 2012).The regional erosion condition is applied to the entire model surface.This erosional model utilizes simple diffusion conditions: where K e is the diffusional erosion coefficient (Ahnert, 1970;Montgomery & Brandon, 2002) and ∇ 2 z is curvature of the surface topography.Finally, the employed erosional model tracks the movement of material and accounts for deposition of sediment within basins.The rate of diffusional erosion is dependent upon: 1) the curvature of surface topography, 2) erosion coefficient, 3) the elevation of the node, 4) the surface geology (e.g.sediment or bedrock), 5) the elevation of ELA, and 6) the presence of glaciers (NSDIC, 1999; Table 3).The K e values are shown in Table 3 and yield maximum erosion rates consistent with observed rates (e.g.Hallet et al., 1996).3.75 m/year 1) Example erosion rate given for 2 z = 1.
For most models with applied erosion conditions, the ELA is set to 2500 m above mean sea level.However, to demonstrate the effects of temporally variable erosion rates driven by periods of glacial advance, a model was developed that includes variable ELA (Table 4).It should be noted that the goal of this model was not to reproduce the observed glacial pattern within southern Alaska, but to apply variable erosion conditions that would provide model results that illustrate their effects on rates and patterns uplift and strain.

Comparison of Model Results for Varied Surface Processes
The series of models results shown here used variable surface boundary conditions, including: 1) free surface (Reference), 2) diffusion-based erosion (Erode), and 3) combined diffusion-based and glacial erosion (Local Erode).Generally, as seen by the vertical velocity (e.g.uplift; Figure 3) and shear strain rate (Figure 4) results at the end of the 5 Ma long model simulations, the inclusion of erosion (either diffusion-based or focused glacial erosion) caused only minor changes to the resultant uplift rate and strain partitioning (Figures 3 and 4).The overall pattern of uplift (Vz; mm/year) was similar for each model simulation irrespective of surface boundary conditions (Figure 3).The Reference Model (Figure 3A), which included a free surface without erosion, develops a zone of uplift (4-5 mm/year peak rate) extending from the Yakutat Bay (YB) area westward.The geometry of the uplifted zone changes from a narrow zone of uplift with relatively high amplitude to a wide, low amplitude zone of uplift westward.Finally, a zone of relatively low uplift extends southward from YB toward Mt.Fairweather (MF).
The strain pattern develop by this model consists of a zone of high strain along the southern transpressive boundary (Figure 4A).This high strain zone splays to a series of compressive faults north of MF.The zone of strain is bound to the south by a northward-dipping shear zone in the vicinity of the coast and to north by a southward-dipping shear zone in the vicinity of Mt.Logan (ML).These two shear zones merge at depth to form an orogenic wedge structure that bounds the area of uplift.The strain partitioning changes from the southern shear being dominant in the eastern part of the model to the northern shear zone being dominant in the west.
The inclusion of a simple diffusion-based erosion scheme (model Erode) did not alter the spatial pattern of uplift or strain (Figs.3B and 4B) in minor ways.However, in inclusion of both diffusion-based and local glacial erosion (model Local Erosion) caused an increase in the peak uplift rates to approximately 5-6 mm/yr (Figure 4C) and the zone of peak uplift occurred over a much larger area than the Reference model (Figure 3A).The spatial partitioning of strain within this model was similar to the results of the Reference model (Figure 4B).The main differences were that the Erosion model favored the accumulation of strain further to the west on the southern compressional zone and further to the east on the northern compressional zone (Figure 4C).
Overall, the patterns of uplift and strain developed in the Reference model (Figures 3A and 4A) held through all model simulations.The general pattern of uplift was altered by the inclusion of glacial erosion conditions by increasing the rate of uplift, widening the uplifted orogen, and shifting the zone of uplift toward the indentor (Figure 3).The strain pattern was characterized in all models by a lateral shear zone along the eastern boundary that splays into a northern and southern shear zone that define the orogenic wedge.That overall pattern was altered through the inclusion of erosion by transferring strain to the southern shear zone and shifting this zone of strain progressively toward the indentor (Figure 4).

General Temporal Changes in All Models
All model simulations had 5 Ma run times.Each model shows variation within rates and patterns of strain and uplift during the first 2 Ma and little variability after 2 Ma, with one notable exception.The model with spatially and temporally variable erosion shows significant variation in the location and rates of uplift and strain through the entire model run as a result of variable erosion conditions (See Section 3.2 and Table 4).The temporal changes in uplift and strain rates are given for the 5 Ma simulations for the Local Erode (Figure 5) and Variable Erode (Figure 6) models.Both models have similar applied erosion conditions, with the only difference is the Variable Erode model alters the ELA as the model simulation proceeds.

Temporal Changes in Strain Partitioning for Local Erode Model
In order to demonstrate the effects of temporally variable erosion on the partitioning of strain within the orogen, the results for uplift and shear strain (Figure 5) are given for the Local Erode model.Generally, as noted in Sections 4.1, the spatial partitioning of strain within all of the models reflects the basic pattern that is developed in the Reference model.The pattern of uplift is fully established within the first 2 Ma (Figure 5A) and the long-term rate of uplift by 3 Ma (Figure 5C).The early stages of model development, through 2 Ma, display a moderately slower rate of uplift (~3.5 mm/yr versus ~6 mm/yr; Figs.5A and 5C).Additionally, the zone of uplifted area, which defines the orogenic wedge structure, increases in width spatially from east to west and temporally as the model evolves (Figure 5A, 5C, 5E and 5G).The general pattern of strain is also developed early during the model run (Figure 5B).Strain is progressively partitioned into the northern and southern shear zones as the model progresses (Figs.5B, 5D, 5F, and 5H).Generally, the strain and uplift patterns are developed early in the model run and only the rates change as the model evolves.

Temporal Changes in Strain Partitioning for Variable Erode Model
The Variable Erode model includes a spatially and temporally variable erosion model.The coefficient of diffusion is variable based upon glacial cover, slope, elevation, sediment thickness, and ELA (Section 3.2 and Tables 3 and 4).The ELA and glacial cover vary with time as prescribed in the model code.
The pattern of uplift associated with model (Figure 6) was significantly different from the other models (Figs. 3  and 5).All models develop the same elongated zone of peak uplift, however, for this model the zone was shifted southwards towards the indentor and has a slightly higher rate developed much earlier in the model run (~7 mm/year at 2 Ma).The uplift pattern was established by 2 Ma (Figure 6A) with peak rates of 6-7 mm/year occurring the furthest south of any model.This maximum in rate was observed at 2 Ma following an extended period when the ELA dropped to sea level (Figure 6A and 6C and Table 4).At 4 Ma, moderate uplift rates (2-3 mm/yr) are observed south of the present day coastline (Figure 6E) and continue through 5 Ma (Figure 6G).This change in the uplift pattern is also consistent with a drop in the elevation of the ELA to sea level (Table 4).At 5 Ma, the maximum uplift rate shifts away from the corner structure or vicinity of ML further to the west (Figure 6G).

Temporal Changes in Strain Partitioning for Variable Erode Model
The Variable Erode model includes a spatially and temporally variable erosion model.The coefficient of diffusion is variable based upon glacial cover, slope, elevation, sediment thickness, and ELA (Section 3.2 and Tables 3 and  4).The ELA and glacial cover vary with time as prescribed in the model code.The pattern of uplift associated with model (Figure 6) was significantly different from the other models (Figs. 3  and 5).All models develop the same elongated zone of peak uplift, however, for this model the zone was shifted southwards towards the indentor and has a slightly higher rate developed much earlier in the model run (~7 mm/year at 2 Ma).The uplift pattern was established by 2 Ma (Figure 6A) with peak rates of 6-7 mm/year occurring the furthest south of any model.This maximum in rate was observed at 2 Ma following an extended period when the ELA dropped to sea level (Figure 6A and 6C and Table 4).At 4 Ma, moderate uplift rates (2-3 mm/yr) are observed south of the present day coastline (Figure 6E) and continue through 5 Ma (Figure 6G).This change in the uplift pattern is also consistent with a drop in the elevation of the ELA to sea level (Table 4).At 5 Ma, the maximum uplift rate shifts away from the corner structure or vicinity of ML further to the west (Figure 6G).The strain pattern that develops is similar to other models; however, the southern shear zone is dominant and shifted further to the south (Figure 6B, 6D, 6F, and 6H).Additionally, apart from the high strain rates along the lateral boundary, there is a zone of high strain in the location of present day YB that develops by 2 Ma (Figure 6B) and remains through the entire 5 Ma model run (Fig 6G).The strain rates are relatively constant within the major zones of strain from 3 Ma to the end of the model run (Figs.6D, 6F, and 6G).Finally, there is a significant amount of strain (with respect to the other models), albeit with low rates, found south of the main zone of strain following the ELA change from 2.5 Ma to 3.5 Ma (Figure 6F and 6G).

Strain Partitioning for Variable Erode Model
The partitioning of strain within the model can be seen within the individual components of the strain rate tensor: ἐ XY , ἐ YZ , and ἐ XZ (Figure 7).The eastern lateral boundary is well defined in the model results by a dextral transpressive zone as shown by the horizontal component ἐ XY (Figure 7A).This dextral motion splays into a series of oblique motion thrust faults to the north of present day YB (Figure 7B and 7C).The two vertical components show thrust faulting that verges eastward (negative values of ἐ XZ ), westward (positive values of ἐ XZ ), northward (positive values of ἐ YZ ) and southward (negative values of ἐ YZ ).The coastal shear zones are evident on Figure 7B.The offshore fault consists mostly motion in the YZ plane, whereas the thrust fault along the coast has a considerable component of motion along the XZ plane (Figure 7C), especially to the east.The peak rates for vertical strain rates correspond to the area around YB where the transpressional motion of the lateral shear zone is partitioned into oblique compressional strain of orogen further to the west.Finally, minor oblique motion to eastern directed thrust faulting occurs inboard (north and east) from the main zone of strain within the orogen (Figure 7C).

Discussion
The primary goal of this study was to demonstrate how coupled feedback between tectonics and surface processes affect deformation patterns in the active St. Elias Orogen.The simulation results indicate that the patterns of uplift and strain (Figs. 3,4,and 5) are controlled on the first order by the structure of the boundary and on the second order vary systematically as a function of surface boundary conditions.In all cases the eastern lateral boundary is defined as a dextral transpressive system.The main differences are within the main zone of convergence.Simulations that do not include erosion develop a dominant backthrust structure, whereas those that include erosion favor the development a dominant frontal thrust.When variable erosional conditions are included, such as those that would occur with variable glacial cover, the localization of uplift and strain changes to reflect the changing surface processes (Figs. 3 through 7).This series of simulations demonstrate the role played by erosion in how strain is partitioned, which is consistent with the interpretation presented by Chapman et al (2008).

Evolution of Orogenic Geometry
The orogenic width developed by the model has a systematic change from the narrow Fairweather Range associated with the dextral transpression, into the corner structure, and widening of the orogen into an evolving fold and thrust belt further to the west, which can be seen in the simulated uplift rates (Figs. 3,5,and 6).The natural topography (mean elevation) mirrors the structural width of the orogen as the distribution of shortening changes with structural style westward along strike (Meigs & Sauber, 2000;Bruhn et al., 2004) and strain becomes more distributed into a wider foreland fold-and-thrust belt.This is mirrored in the simulation results as the orogen widens with a decrease in magnitude of surface uplift as strain is distributed over wider area (Figs. 1,3,5,and 6).

Structural Evolution of the Orogen
The model develops a strain pattern that at least partially corresponds to large-scale mapped structures in southern Alaska (Figs. 4 through 7).The most obvious correlation is the development of the lateral dextral transpressional shear zone that corresponds to the Queen Charlotte-Fairweather Fault zone (Figs. 1, 4, and 7).As mass is carried into the corner by the oblique part of the orogen (Figs. 4 and 7), strain is transferred westward into a bivergent thrust system that is broadly consistent with observed deformation within the St. Elias Range, where the shear zones define oppositely dipping thrust faults located south and northeast of the highest peaks, Mt.
St. Elias and Mt.Logan, respectively (Figs. 4 and 7).The transition from transform to contractional strain results in a maximum in strain rate within the corner structure (Koons et al., 2010), which is also the modeled area where vertical flow causes the greatest rate of uplift (Figure 3) and exhumation rates.This can be seen within the individual components of the strain rate tensor: ἐ XY , ἐ YZ , and ἐ XZ (Figure 7).The dextral transpression continues further to the northwest (Figs. 4 and 7), corresponding to the feature that links the deformation within the St.
The models predict development of a southward vergent thrust system (Figs. 4 and 7), which coincides with the well-recognized fold-and-thrust belt (e.g.Plafker et al., 1994) and the Malaspina Fault and its offshore continuations (e.g.Worthington et al., 2009).The secondary southward vergent thrust system is located to the north of Mt.St. Elias within the vicinity of the known Contact fault and/or proposed back-thrust (Figure 7; Berger et al., 2008a,b), both located under the 15 km wide and >100 km long Bagley Ice field.Additionally, the model develops a northward vergent thrust (Figure 7) at the northern flanks of the mountain range coinciding with the Chugach terrane and the Border Range Fault system at the topographic front of the Chitina valley.There is little evidence of active deformation along the northern flanks from either neotectonics (e.g.Plafker et al., 1994) or GPS studies (e.g.Freymueller et al., 2008), but northward increasing cooling ages suggest deformation over longer time scales (Enkelmann et al., 2008), and high fracture density suggests active deformation (Champagnac et al., 2012).Spotila and Berger (2010) presented low-temperature cooling ages for the St. Elias orogen that show a symmetrical pattern of rock uplift around the eastern syntaxial corner.The authors suggest that partitioning of transcurrent deformation into this obliquely transpressive restraining bend is driving the observed deformation and exhumation within the corner.This model is consistent with the numerical results presented here where a relatively narrow zone of uplift occurs along the eastern lateral boundary of the indentor culminating in a bulls-eye pattern of high uplift within the corner.The transition from this lateral oblique transcurrent strain to contractional strain within the fold and thrust belt to the west is evident within the shear strain components of the strain rate tensor (Figure 7a and 7b).The numerical model results suggest that the backstop of deformation associated with the Yakutat convergence and collision is located just south of the Chitina River valley and thus coincides roughly with the Border Range fault zone, the suture between the Chugach terrane and the Wrangellia terrane.This interpretation is in agreement with detrital zircon fission track and K-feldspar modeling data that suggest that there has been no significant cooling due to exhumation of the Wrangellia terrane rock west of Mt.Logan since the Late Cretaceous (Enkelmann et al., 2010).

Role of Erosion
The thermal-mechanical simulations presented here have captured most of the variance in the signal of accretionary tectonics in southeast Alaska (Figs. 3 and 4).The results show the formation of strain rate maxima in the tectonic corner that spatially coincides with the observed zone of rapid exhumation in the Seward Glacier area (Figs. 3,5,and 6).Inclusion of erosion conditions alter these tectonically developed strain patterns and captures the evolution of local topography, observed fault zones, and cooling age patterns (Figs. 4 and 5).In particular, the simulations reveal focused uplift within the tectonic corner, which is roughly coincident with the present highest topography (Mt.Logan and Mt St. Elias) and an extended ice field region of the Bagley Ice Field.
The spatial coincidence of the corner structure with the transition from dextral shear strain along the Fairweather Range to the convergent strain associated with the St. Elias orogenic wedge in all models indicates deformation and uplift patterns are dominated by the tectonic signal.This pattern of focused strain demonstrates the first order control provided by the tectonic geometry on the focusing of strain and the secondary influence of topographic load and surface processes.
The simulations produce a focused zone of uplift that is consistent with the expected area of focused exhumation underneath the Seward Glacier (Figs. 3,5,and 6).The results further suggest that the region of intense uplift at the Yakutat corner is larger than just the Seward Glacier catchment and extends further northward and eastward.Analysis of deformed syncollisional deposits within the fold-and-thrust belt revealed a localized exhumation zone with cooling rates similar to the current rates existed 3 to 4 Ma (Enkelmann et al., 2008).Assuming that this sediment source was the same as today indicates that the observed zone of exhumation at the Yakutat indenter corner has existed since the Pliocene and is thus not solely a tectonic response that developed due to Quaternary climate change as suggested in the case for the evolution of the fold-and thrust belt (Berger et al., 2008b).
Including the effects of local glacial erosion, these simulations predict that a local perturbation of the thermal structure sufficient to produce a thermochronological record similar to the observed zircon fission track cooling age data can develop within 5 Myr.This zone of high exhumation is formed in all simulations, indicating that its development is related to tectonics and is not likely due to surface processes.

Limitations of the Model
The development of the model requires a balance between computational efficiency and resolution, which can potentially have strong effects on simulation results.Discretization that is too coarse may not capture the system sufficiently, and too fine may lead to premature computational failure.As such, the resolution chosen for the model is such that groups of individual structures (e.g.fold and thrust belts) are discernable as broad (10-50 km) shear zones.The sensitivity analysis completed on the variable discretization indicated that there were no induced errors associated with the range of grid spacing that were used and only the resolution of the solutions was affected.
The diffusion based erosion scheme that is used in the simulations is simplistic.The rate of erosion applied to the model is based upon local topography and is not a function of rock type.Additionally, the resolution is restricted to the 10 km horizontal grid spacing of the numerical solutions.That discretization is larger in scale than individual glaciers and therefore cannot adequately develop glacial landscapes.While the erosion model is relatively simplistic, it meets the purpose for which it was designed by displaying the general linkages between the large-scale landscape development processes and underlying tectonics.
Finally, the large-scale geological structure of the region is not included in the model.Due to numerical constraints, the structures of the Yakutat terrane and North America (e.g.thickness, rock type, pre-existing structures) are not included within the model geometry.These structural heterogeneities could provide a strong control on strain partitioning, however, the older structures within the study area appears to have little current strain accumulation and likely do not affect partitioning of strain within the orogen.

Conclusions
The framework of the three-dimensional models presented here allows for the investigation of coupling between tectonics, rheological and thermal evolution, and surface processes in the evolution of an active orogen during the early stages of convergence.The strain and uplift pattern observed within southern Alaska has been interpreted as a dynamically coupled system where efficient erosion has enhanced the general pattern of strain and uplift (Pavlis et al., 2004;Berger et al., 2008a;Koons et al., 2010).This series of numerical models demonstrate that the strain patterns within Southern Alaska are characteristic of oblique convergence within a syntaxial corner structure where dextral transpressive strain to the south along the Queen Charlotte and Fairweather fault systems is partitioned into the active fold-and-thrust belt of the St. Elias orogen.Additionally, the model results indicate that the presence of an erosional surface effects strain patterns, differing in both strain rate and spatial partitioning, and indicates that temporally and spatially variable focused erosion can drive a structural reorganization similar to that seen in the St. Elias orogen.

Figure 2 .
Figure 2. Simplified 3-D cartoon plot of the model geometry, including topography Note.The black dashed line in the bottom level indicates the region defined as the Yakutat Terrane projected onto the surface, to the southeast onto which the initial velocity condition (vectors; Yak velocity from Figure 1; Fletcher & Freymueller, 1999) are imposed at the base of the model.The model mechanics and dimensions are as noted.See text for additional details.

Figure 3 .
Figure 3. Contour plots of Vz (uplift; mm/year) taken from the a) Reference model, b) Erosion model, and c) Local Erosion model at the model surface.Note.The location of important landmarks (Mt.Logan, Mt.St. Elias, and Mt.Fairweather) and the coastline (dark line) are shown for spatial reference; please refer to Figure 1 for further details.

Figure 4 .
Figure 4. Contour plots of strain rate (nStrain/year) taken from the a) Reference model, b) Erosion model, and c) Local Erosion model at the model surface Note.The location of important landmarks (Mt.Logan, Mt.St. Elias, and Mt.Fairweather) and the coastline (dark line) are shown for spatial reference; please refer to Figure 1 for further details.

Figure 5 .
Figure 5. Contour plots of Vz (uplift; mm/year) and strain rate (nStrain/year) for the Local Erosion model taken from the a) Vz at 2 Ma, b) strain rate at 2 Ma, c) Vz at 3 Ma, d) strain rate at 3 Ma, e) Vz at 4 Ma, f) strain rate at 4 Ma, g) Vz at 5 Ma, and h) strain rate at 5 Ma, all taken at the model surface Note.The location of important landmarks (Mt.Logan, Mt.St. Elias, and Mt.Fairweather) and the coastline (dark

Figure 6 .
Figure 6.Contour plots of Vz (uplift; mm/year) and strain rate (nStrain/year) for the Variable Erosion model taken from the a) Vz at 2 Ma, b) strain rate at 2 Ma, c) Vz at 3 Ma, d) strain rate at 3 Ma, e) Vz at 4 Ma, f) strain rate at 4 Ma, g) Vz at 5 Ma, and h) strain rate at 5 Ma, all taken at the model surface Note.The location of important landmarks (Mt.Logan, Mt.St. Elias, and Mt.Fairweather) and the coastline (dark line) are shown for spatial reference; please refer to Figure 1 for further details.

Figure 7 .
Figure 7. Contour plots of the strain rate tensor components in the a) horizontal (ε XY ) plane, b) vertical North-South (ε YZ ) planes, and c) vertical East-West (ε XZ ) planes from the Variable Erosion model at the model surface.Negative values in a) track locations of dextral transcurrent shear strain, while in b) and c) positive values are northward/eastward thrusts and negative values are southward/westward thrusts Note.The location of important landmarks and the coastline (dark line) are shown for spatial reference; please refer to Figure 1 for further details.

Table 1 .
Material properties and constants assigned to the thermo-mechanical numerical models.

Table 2 .
Applied Surface Boundary Conditions

Table 3 .
Erosion model parameters used to estimate K e

Table 4 .
Variable ELA conditions for the Variable Erode model.