SUBDUCTION ACROSS TIMESCALES: GEODETIC AND SEISMIC INVESTIGATIONS OF ONGOING AND TERMINAL-STAGE SUBDUCTION IN ALASKA AND NORTHWEST CANADA By Connor Drooff A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Earth and Environmental Sciences—Doctor of Philosophy 2024 ABSTRACT Tectonic activity in Southeast Alaska and Northwest Canada is a consequence of the relative motion between the North American and Pacific plates. In this dissertation we present the results of three geodetic and seismic studies that investigate earthquake processes in the northernmost section of the Cordilleran orogenic belt within the Alaskan and Northwest Canadian convergent margins. Chapter 2 focuses on geodetic coupling of the Aleutian trench outboard of the Alaska Peninsula. Geodetic coupling is a measure of how elastic strain accumulates along a subduction interface between megathrust (M > 9.0) earthquakes and is a critical tool for monitoring earthquake hazards in subduction zones. We present evidence that the subduction interface in this area is comprised of a sequence of discrete segments and show how these correspond with observed rupture characteristics of the 2021 M 7.8 Simeonof earthquake. Chapter 3 explores how the segments identified in Chapter 2 have been updated in successive work and how the coupling segments first presented in Chapter 2 are shown to closely correlate with the 2021-2023 Alaska Peninsula earthquake sequence. Chapter 4 transitions some 1500 km to the northwest into the Northern Canadian Cordillera (NCC). This complex tectonic setting is undergoing active tectonic deformation characterized by diffuse patterns of seismicity. We use broadband seismic data from Macken- zie Mountain Earthscope Project (MMEP), USArray Transportable Array (TA), and per- manent Canadian National Seismic Network stations to present a local earthquake catalogue with high sensitivity to small regional earthquakes. Deep learning techniques are adopted for both seismic phase detection and event association. Event relocations are performed to provide well constrained estimates of earthquake depth distributions. Clusters of microseis- micity spanning the upper crust are located in the central Richardson Mountains, in the vicinity of the Tintina fault, and in the northeast Selwyn Basin. These clusters suggest that the core of the Richardson Anticlinorium is tectonically active and that the Tintina fault is a locus for low levels of active deformation. We interpret seismicity in the northeast Selwyn Basin as partially occurring along the Plateau thrust at depth or in its hanging wall. We suggest that some combination of localized duplex structures and lithological strength con- trasts both within the Selwyn Basin and between the Selwyn Basin and abutting Paleozoic shelf sequences may be responsible for seismicity in the Mackenzie Mountain foreland. Chapter 5 is also situated within the NCC and investigates GPS site velocities in order to place constraints on measurable geodetic deformation patterns occurring within the moun- tain belt. First order observations of site velocities are on the order of 1-3 mm/yr. We use these to distinguish between 3D predictions of GIA models for the effects of both post-Last Glacial Maximum deglaciation on the North American continent as well as post-Little Ice Age deglaciation in Southeast Alaska. When GPS site velocities are corrected for the 3D effects of GIA a clear effect of strain accumulation at the level of 3-5 mm/yr is observed which extends at least as far inboard as the Mackenzie Mountain foreland. We use these new site velocities to constrain a regional tectonic block model and evaluate whether the Tintina fault and Southern Richardson Mountains could plausibly define present-day block boundaries. Our findings suggest that inclusion of the Tintina fault represents a significant improvement in fit to GPS site velocities and that maximum allowable strike slip and tensile deformation rates on the Tintina fault are 1 mm/yr. We also find that GPS site velocities in the Richardson Mountains favor a more southward directed block motion relative to sites further south, however the inclusion of a separate Richardson block does not represent a significant improvement in model fit based on currently available data. Copyright by CONNOR DROOFF 2024 ACKNOWLEDGEMENTS The work contained in this dissertation is the culmination of 6+ years of graduate studies at two separate institutions. The list of people that have been influential to me in both the professional and personal spheres in that time is too numerous to exhaustively fit here. My family including my mother Marilyn and father Michael first and foremost have been a bedrock of support and encouragement throughout my graduate career. It is hard to imagine that I would have been able to finish up a doctoral degree without their presence in my life. My sister Carly and Aunt Pam have also been great outlets of emotional support and lively conversation. My friend group from the time in Alaska including Julia, Ragen, Jordan, and Cole did so much to make Fairbanks feel like a home to me and although our time together was relatively brief I have found a lot of meaning in maintaining those friendships as the years have gone by. My Michigan friends Gabriel, Jackson, and Sahira among many others have made my time in Lansing so rich and it is hard to imagine completing graduate school without them in my life. My partner Julia has been a boundless source of support and encouragement and I am forever grateful to all that she has done for me during my time in graduate school. And lastly my cat Hugo has been an absolute bedrock in my life throughout the last 10 years and has been the one constant throughout 4 moves in graduate school, two of which were across the country. v 1 7 29 37 TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 2 NEW CONSTRAINTS ON SLIP DEFICIT ON THE ALEUTIAN MEGATHRUST AND INFLATION AT MT. VENIAMINOF, ALASKA FROM REPEAT GPS MEASUREMENTS . . . . . . . CHAPTER 3 UPDATES ON THE ALEUTIAN MEGATHRUST COUPLING MODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 4 CHAPTER 5 NEW INSIGHTS INTO THE ACTIVE TECTONICS OF THE NORTHERN CANADIAN CORDILLERA FROM AN ENHANCED EARTHQUAKE CATALOG . . . . . . . . . . . . . . . . . . . . . GEODETIC STRAIN TRANSFER, TECTONIC BLOCK MOTION, AND 3D GLACIAL ISOSTATIC ADJUSTMENT IN THE NORTHERN 72 CANADIAN CORDILLERA . . . . . . . . . . . . . . . . . . . . . CHAPTER 6 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 vi CHAPTER 1 INTRODUCTION 1 Tectonic activity in Southeast Alaska and Northwest Canada is a consequence of the relative motion between the North American and Pacific plates. The direction of conver- gence between the plates is oriented towards northwest in the northern Gulf of Alaska and serves as a fundamental control on deformation. In the southeast Gulf of Alaska offshore of British Columbia the tectonic regime is defined by right-lateral strike-slip motion that transitions to the northwest into classic oceanic-continental subduction in the Aleutian arc. This heterogeneity gives rise to a variety of geological processes which can shed light onto the spatial scales to which tectonic strain is able to penetrate into a continental interior. In this dissertation we present the results of three geodetic and seismic studies that investi- gate earthquake processes along the Alaskan and Northwest Canadian margin. The tectonic processes which are presented and modeled in this dissertation span the range of spatiotem- poral scales over which subduction can evolve over geologic time. Subduction zones host the vast majority of naturally-occurring earthquakes observed globally and are thus of primary importance when it comes to understanding and mitigating earthquake and tsunami-related hazards. The earthquake cycle at subduction zones is capable of producing large (Mw > 9.0) megathrust earthquakes which have recurrence intervals on the order of 10’s to 1000’s of years. With each event that ruptures the plate interface meters to 10’s meters of oceanic lithosphere is pushed slightly deeper into the underlying mantle. The continuity of this pro- cess on timescales of 10’s to 100’s of millions of years allows for entire oceanic plates to be cycled. The Pacific plate is the largest tectonic plate on the surface of the Earth at present. It underlies the vast majority of the Pacific ocean basin and spans a total area of ∼ 103 million km2 or roughly 22% of the Earth’s surface. This has not always been the case. Tectonic reconstructions of prior oceanic plate configurations suggest that the Farallon plate could have been as large or potentially larger than the modern day Pacific plate (citations). Starting during the breakup of supercontinent of Gondwana ∼ 200 Ma the Farallon plate was progressively consumed beneath the western margin of ancestral North America, giving 2 shape to many of the features of Cordilleran mountain belt as they exist today. Global seismic tomography models have produced images of the Farallon slab in the mantle beneath the North American continent, suggesting that it may extend as far down into the Earth’s interior as the core-mantle boundary (Grand, 2002; Van Der Lee and Nolet, 1997). As it converged with North America, the Farallon plate acted as a tectonic conveyor belt and was responsible for amalgamating a complex sequence of exotic terranes and island arcs to the margin. These accreted terranes, today only fractionally preserved in the exposed geologic record, hint at great complexity in the configuration of the oceanic lithosphere outboard of ancestral North America. The Mesozoic configuration of island arcs, oceanic plateaus, and microcontinents out- board of North America has been likened to the modern Western Pacific (Monger and Price, 2002). These features were a result of complex interactions between variably sized sections of oceanic lithosphere (Sigloch and Mihalynuk, 2013). Plate reconstructions suggest that by the Late Cretaceous a triple junction of oceanic spreading ridges as well as an offshore subduction zone similar the modern Sea of Japan (Monger, 1997) was located in the Eastern proto-Pacific. The triple junction separated the Farallon plate in the West from the Kula plate in the north and the Izanagi plate to the west. The Farallon and Kula plates as well as the spreading ridge that existed between them were progressively consumed beneath the ancestral North American margin throughout Mesozoic time. This process was responsible for elevating the thermal state of the mantle lithosphere beneath Western North America in both the Northern Canadian Cordillera as well as in the modern day basin and range province in Utah/ Nevada. This elevated thermal state is thought to be responsible for active deformation in the NCC. The inherited structures from Mesozoic convergence are lo- cated throughout the Mackenzie Mountains of present day Yukon and Northwest Territories and are presumably being reactivated by the inboard stress transfer caused by high heat flow. 3 As the subduction of the Farallon and Kula plates progressed the body forces which held them together became lower in magnitude and led to fragmentation and realignment of the configuration of oceanic lithosphere. This was a long-lived process which took 10’s of millions of years to play out, but eventually led to the initiation of subduction along the Aleutian arc at ∼ 60 Ma as well as the fragmentation of the Farallon plate into the modern day Juan de Fuca, Nazca, and Cocos plates. The relict triple junction which was active during Mesozoic time was never fully consumed, but has instead been transported northwards and became a coherent unit within the modern day Pacific plate (Atwater and Menard, 1970). It is presently located outboard of the Shumagin Islands and the southeast Alaska peninsula and is postulated to have a strong influence on the seismic potential of the Alaska-Aleutian subduction zone (i.e. Shillington et al., 2015). This is thought to be due to the alignment of rift-related normal faults inherited from Mesozoic spreading, which are now being reactivated by bending forces within the Pacific plate as it subducts, influencing the efficiency by which fluids are able to enter into the downgoing plate. This is shown to be a control on interseismic coupling characteristics in chapter 2 of this dissertation. Chapter 2 presents a study of geodetic coupling of the Aleutian trench outboard of the Alaska Peninsula. Geodetic coupling is a measure of how and where elastic strain is actively accumulating along a subduction interface. This strain is mostly recoverable and often referred to as ’slip deficit’, which can be thought of as energy that is available to be released in the form of future great earthquakes. The Aleutian trench outboard of the Alaska peninsula lies directly west of the 1964 M 9.2 Great Alaska earthquake rupture patch and displays a remarkable transition from full coupling behavior beneath the Kodiak archipelago to full creeping behavior beneath the Sanak islands over a distance of 500 km. This region has garnered particular scientific attention in recent years due to a complex series of earthquakes along the subduction interface that include the July 22, 2020 M 7.8 Simeonof earthquake, the July 28, 2021 M 8.2 Chignik earthquake, and most recently the 2023 M 7.2 July 16, 2023 Alaska Peninsula earthquake. The Chignik event was the largest earthquake to occur 4 in the US since 1965 and the 7th largest recorded event in US history. We present evidence that the subduction interface in this area is comprised of a sequence of discrete segments, which are likely the result of variations in hydration state at the plate interface. This is influenced by alignment of tectonic fabrics within the incoming Pacific plate which are being reactivated by bending forces as the plate is subducted. Chapter 3 further elaborates on contributions to publications that directly stemmed from the model presented in Chapter 2. We show that by adjusting the assumptions made regarding the downdip distribution of plate coupling that we were able to present a slip deficit model that aligns very closely with observed rupture behavior of the 2021 M 7.8 Simeonof earthquake. This updated slip model also closely aligns with the rupture of the 2022 M 8.2 Chignik earthquake. The breakup of the Farallon plate (∼ 60 Ma) was also coeval with the creation of an oceanic plateau called the Yakutat block outboard of modern day Vancouver island (Bruns, 1983). The Yakutat block has since been progressively transported northwards along the Queen Charlotte margin and is now accreting to and subducting beneath the southern Alaskan margin. It is presently a primary driver of active tectonics in Alaska and Northwest Canada. Chapter 4 focuses on the Northern Canadian Cordillera (NCC), located hundreds of km northeast of the Yakutat collision zone. The tectonic history of the NCC records episodes of subduction dating back well into the late Triassic, most notably the consumption of the Kula-Farallon spreading ridge and opening of a slab window at ∼ 60 Ma(Thorkelson and Taylor, 1989). This has profoundly influenced the thermal state and strength of the lithospheric mantle which is observed to be unusually hot and is modeled to be capable of transferring tectonic stress far inboard into the North American continental interior (Maz- zotti and Hyndman, 2002). We leverage the most dense deployment to date of seismic instruments throughout the region in order to investigate the spatial distribution of low magnitude earthquakes. Our results suggest that the Tintina fault, which was previously assumed to be inactive, hosts many low magnitude events at moderate crustal depths. This represents some of the first evidence of direct observations of tectonic activity along this 5 inherited structure and has implications for earthquake hazard assessments in the vicinity. We also show that low magnitude earthquakes extend northeast well past the Tintina fault into the Mackenzie Mountain foreland. While this area has been previously observed to be seismically active, previous studies have lacked the constraints to propose specific structures that are actively deforming. We show that low magnitude Chapter 5 is also situated within the NCC and investigates GPS site velocities in order to place constraints on measurable geodetic deformation patterns occurring within the moun- tain belt. First order observations of site velocities are on the order of 1-3 mm/yr. We use these to distinguish between 3D predictions of GIA models for the effects of both post-Last Glacial Maximum deglaciation on the North American continent as well as post-Little Ice Age deglaciation in Southeast Alaska. When GPS site velocities are corrected for the 3D effects of GIA a clear effect of strain accumulation at the level of 3-5 mm/yr is observed which extends at least as far inboard as the Mackenzie Mountain foreland. We use these new site velocities to constrain a regional tectonic block model and evaluate whether the Tintina fault and Southern Richardson Mountains could plausibly define present-day block boundaries. Our findings suggest that inclusion of the Tintina fault represents a significant improvement in fit to GPS site velocities and that maximum allowable strike slip and tensile deformation rates on the Tintina fault are 1 mm/yr. We also find that GPS site velocities in the Richardson Mountains favor a more southward directed block motion relative to sites further south, in agreement with predictions of geodynamic models. 6 CHAPTER 2 NEW CONSTRAINTS ON SLIP DEFICIT ON THE ALEUTIAN MEGATHRUST AND INFLATION AT MT. VENIAMINOF, ALASKA FROM REPEAT GPS MEASUREMENTS 7 ABSTRACT We employ an enhanced dataset of GPS velocities to reassess the slip deficit along the plate interface of the Alaska-Aleutian subduction zone. An examination of velocities of sites near Mt Veniaminof, located in the southwest Alaska peninsula, shows that existing models for tectonic deformation do not accurately predict the local velocity field along this section of the peninsula. In a combined model, we solve for the volcanic inflation signal on the edifice of Mt. Veniaminof and reassess the variations in trench coupling outboard of the volcano. The interseismic deformation near Veniaminof and the eastern Shumagin islands requires an additional model segment with higher slip deficit than the western Shumagin islands. This segment appears to have been the main rupture patch of the 2020 M 7.9 Simeonof earthquake. The volcano inflated over the time period of 2005-2017, although at a lower rate than in 2002-2005. Plain Language Summary Subduction zones globally display marked variations in slip behavior along their length. The subduction zone off the coast of the Alaska Peninsula is one such example, where the shallow section of the plate boundary transitions from being fully locked to fully creeping over the span of roughly 250 km. We show that surface deformation data gathered from the region near Mt. Veniaminof and in the eastern Shumagin islands indicates a particular region of ’sticking’ behavior, which appears to have been the focus of rupture in the 2020 M 7.8 Simeonof Earthquake. We find that the data from the eastern Shumagin islands are best explained by these sites lying directly on top of an abrupt along-strike boundary in the slip behavior of the subduction zone, providing an important observation of how the earthquake cycle manifests itself in regions where plate coupling is not uniform along the length of the subduction zone. 2.1 Introduction The Alaska-Aleutian subduction zone accommodates relative motion between the Pacific and North American plates along the southern margin of Alaska. This convergent bound- 8 ary transitions along strike from that of Yakutat oceanic plateau collision and accretion in the northwest Gulf of Alaska to subduction beneath mainland Alaska and the Aleutians (Eberhart-Phillips et al., 2006). The dip of the subduction interface is more shallow in the Kodiak segment, gradually increasing westward (Hayes et al., 2018). Other properties of the megathrust vary along strike, in particular the distribution of slip deficit. Slip deficit is a measure of the degree to which slip in the shallow region of the subduction interface is impeded by frictional coupling between the downgoing and overriding plates. Generally speaking, highly coupled sections of a megathrust accumulate high levels of tectonic stress, which is available to be released at a later date in the form of large magnitude earthquakes. Much of the Alaska Peninsula segment of the subduction zone ruptured in an Mw 8.3 event in 1938 (Figure 4.1), but the Shumagin segment of the plate interface has not ruptured in a great earthquake since at least the 1700s (Davies et al., 1981). The lack of seismic energy release along this section of the subduction zone has been well documented by previous workers i.e. (Nishenko and Jacob, 1990), (Abers et al., 1995) and has led to it being termed the ‘Shumagin Seismic gap’ due to an apparent absence of great (Mw > 8.0) earthquakes from the historical record. Paleoseismic studies i.e. (Witter et al., 2014) have found no evidence of uplifted marine terraces in the archipelago, suggesting that rupture of the megathrust beneath the islands does not occur in large magnitude events. Previous work has suggested that the along-strike variations in slip deficit are related to changes in hydration of the downgoing plate that are connected to the orientation of the pre- existing seafloor fabric (Shillington et al., 2015); (Li and Freymueller, 2018). The seafloor fabric of the downgoing Pacific plate is inherited from a relict triple junction of spreading ridges between the Farallon, Pacific, and Kula plates (Atwater and Menard, 1970) directly outboard of Veniaminof (Figure 4.1). Subduction of the westernmost seafloor generated at the Farallon-Pacific ridge terminates directly outboard of the Semidi islands, which coincides with the transition from full to partial trench locking. (Shillington et al., 2015) suggested that this transition in the dominant orientation of remanent rift structures exerts a control 9 Figure 2.1 Tectonic setting and velocity observations of the subduction zone outboard of the Alaska peninsula. Slip deficit along the plate interface as determined by (Li and Freymueller, 2018) is shown by gray shading, where black refers to >99% coupled, dark gray to 66-99% coupling, and light gray to 33-66% coupling. Boundaries between segments of the (Li and Freymueller, 2018) model shown as dashed orange lines. Velocities of all sites used in this study are shown relative to the Peninsula block (Li et al., 2016). Rupture zones of the 20th century great earthquakes from (Davies et al., 1981) are shown with thinner black dashed lines. The Pacific plate velocity relative to the Peninsula block shown offshore, and magnetic anomaly patterns on the Pacific plate are shown from the EMAG2 model (Maus et al., 2009). Lower right inset shows regional setting of the study area, with modeled blocks labeled such that BRNG = Bering block, SOAK = Southern Alaska block, PENN = Peninsula block, YB = Yakutat Block, NOAM = stable North American Plate, PCFC = Pacific Plate. Site velocities in the region near Mt. Veniaminof plotted in the top left inset (note different scale). Red vectors in top left inset show the predictions of the (Li and Freymueller, 2018) model. The five sites closest to the volcano have been occupied in 2002, 2005, and 2017. Site velocities include a radial outward influence caused by volcanic inflation. 10 on the hydration state of the downgoing plate, in that the reactivation of faults in crust with the favorably-oriented Kula-Pacific fabric by flexure through the outer rise creates additional pathways for fluids to enter into the downgoing slab. The previous geodetic study of (Li and Freymueller, 2018) excluded all of the data on and around Mt. Veniaminof (inset of Figure 4.1) out of concern that volcanic inflation might bias the estimated slip deficit model. However, this in effect created an along-strike data gap, and we found that their model systematically misfit the velocities of these sites by underpredicting the trench-normal component of velocity field by an average of 1.3 mm/yr for all excluded sites. This motivated a reassessment of the slip deficit and volcanic source models. We added the data from the sites on and near Mt. Veniaminof to the dataset of (Li and Freymueller, 2018), and re-assessed the segmentation of slip deficit along the Alaska Peninsula using an iterative solution to model both the volcanic and tectonic deformation together. 2.1.1 Mt. Veniaminof Eruptive and Deformation History Mt. Veniaminof lies eastward and inboard of the Shumagin segment of the Aleutian trench. It is a broad stratovolcano that measures roughly 35 km in diameter and stands over 2.5 km tall, with a 4 km-wide glacier-capped summit caldera. An eruptive cone on the west end of the caldera has produced small volume basaltic to andesitic flows in historic time (Bacon et al., 2007) and eruptive activity is typically sustained on the order of months with recurrence every few years. The volcano recently underwent an explosive eruption from June to October 2013 that produced a ∼ 5 km tall ash column and several lava flows estimated to total roughly 500 m3 of extruded volume (Dixon et al., 2015). (Fournier and Freymueller, 2008) used repeat GPS measurements at an array of stations on the volcanic edifice to show that Veniaminof exhibited an inflationary geodetic signal between 2002-2005 once the predictions of a coupling model of the Aleutian trench are removed from the data. Five sites on the volcanic edifice were re-surveyed in 2017, and the 2005-2017 velocities were substantially different than those from 2002-2005, suggesting a much smaller volcanic 11 inflation signal. Velocities for those 5 sites and others on the Pacific and Bering Sea coasts of the Peninsula near the volcano are shown in the inset of Figure 4.1. 2.1.2 Previous Modeling Work A westward transition from locked to creeping behavior of the megathrust outboard of the Alaska Peninsula has been observed since the region began to be surveyed by geodetic instruments i.e. (Freymueller and Beavan, 1999), (Fletcher et al., 2001). Further work incor- porating larger datastets has placed additional constraints on the locations of locking vari- ations along the plate interface i.e. (Fournier and Freymueller, 2007), (Li and Freymueller, 2018), which grades from fully coupled near the Kodiak archipelago to fully creeping near Sanak island. The model of (Li and Freymueller, 2018) showed that the data were fit best with discrete and relatively sharp transitions in slip deficit: one aligned with the Semidi is- lands, a second east of the Shumagin islands and another west of the Shumagins marking the transition to fully creeping (Figure 4.1). However, that study did not use data from volcanic sites on Veniaminof in order to minimize the effects of volcanic deformation on the slip deficit model, and this introduced a spatial data gap. A subtle, but persistent underprediction of the trench normal velocities at excluded sites suggests the presence of along-strike variation in the slip deficit not captured in their model. We improve on their model by considering those data along with a model for the volcanic inflation signal. GPS observations on the direct volcanic edifice of Mt Veniaminof are confined to 3 surveys conducted during the summer field seasons of 2002, 2005, and 2017. Using the 2002-2005 data from an array of 10 sites on the volcanic edifice, (Fournier and Freymueller, 2008) showed that a deep Mogi point source or a shallower sill source both fit the data well. The shallower sill source provided the lowest misfit to the data, but the relative reduction in error was small given the increased number of model parameters. It is also important to note that this survey interval represents a relatively short time frame in which the volcano produced a few smaller explosions that involved elevated seismicity and ash emissions but no observed lava flows, and also preceded a 2013 eruptive event which involved the extrusion of 500 m3 of lava 12 flows, so it is reasonable to assume that the 2002-2005 data should record a short term inter- eruptive inflationary signal. SAR interferograms over the volcano generally lack coherence due to snow and ice cover and infrequent temporal sampling. (Lu and Dzurisin, 2014) show that stacking coherent summer-to-summer interferograms from independent satellite passes could produce a reasonable estimate of the long-term deformation field and found that a deep (∼ 8km) Mogi point source provided a reasonable fit to the data. The 2005-2017 deformation signal from the volcano is clearly smaller than that from 2002-2005, so we will use only those data in the models shown in this paper. 2.2 Data We used data from campaign and continuous GPS sites along the Alaska Peninsula, over the time period 1992-2017 i.e. (Freymueller and Beavan, 1999), (Fletcher et al., 2001), (Fournier and Freymueller, 2008) to estimate our combined trench locking and volcanic source model. Repeat campaign GPS measurements at 12 sites near Mt. Veniaminof were performed in 2002 and 2005 (Fournier and Freymueller, 2008). Five of these sites on the direct volcanic edifice were re-occupied in the summer of 2017 and we incorporate these data into our study. We used the GIPSY/OASIS goa-5.0 software package developed by the Jet Propulsion Laboratory (JPL) to obtain daily coordinate and covariance estimates for continuous and campaign GPS stations using point positioning. Data processing methods are discussed in greater detail in (Fu et al., 2012). We adopted JPL nonfiducial orbits and clock products and estimated a transformation for each daily solution into the ITRF2008 reference frame (Altamimi et al., 2011) using our own solutions. We calculated velocities from the time series of each individual station using a weighted least squares approach. We fit linear velocities (Table S1) to all data within the survey period (1992–2017) except in the case of the sites on Mt. Veniaminof, where we fit separate velocities to stations between 2005-2017 and 2002- 2005 because the volcanic signal was variable in time. We use the 2005-2017 Veniaminof velocities in the combined model. The volcanic source is clearly smaller in the later time 13 period than it was in 2002-2005. Due to the temporal sparseness of surveys at some sites we assumed a white noise model for sigma variation in daily point positions in the velocity fits, and scaled the uncertainties based on the post-fit misfit so that the reduced chi-square statistic was equal to 1.0 for each fit. We also removed the geocenter translation rate error in ITRF2008 estimated by (Argus et al., 2010) and calculate velocities in a North America-fixed reference frame defined by the GEODVEL model (Argus et al., 2010). Site motions are well described by linear motion with time, except for the sites on Veniaminof volcano itself (see the supplement of (Li and Freymueller, 2018)). Figure 1 shows the velocities with the Peninsula block motion of (Li and Freymueller, 2018) subtracted to isolate the subduction and volcanic related signals. Models for Glacial Isostatic Adjustment (Hu and Freymueller, 2019) and 1964 postseismic deformation (Suito and Freymueller, 2009) were removed from all GPS velocities. 2.3 Combining Models for Slip Deficit and Volcanic Deformation Velocities of sites on and near Mt. Veniaminof fill in a critical data gap. The dominant component of the geodetic signal at these sites is due to slip deficit along the plate interface outboard of the volcano (Figure 4.1), however previous work (Fournier and Freymueller, 2007), (Lu and Dzurisin, 2014) showed that Mt. Veniaminof also exhibits a localized defor- mation field that is the result of processes operating at depth beneath the volcano. In order to gain better insight into both volcanic processes as well as slip deficit on the megathrust we adopted an iterative approach that allows us to model both processes simultaneously. The iterative approach works well because the volcanic deformation is highly localized, and affects only sites on the edifice of the volcano, while the elastic strain from the megathrust is smooth in the vicinity of the volcano as it is far from the trench. Because the volcanic signal is small and localized, uncertainties in the volcanic source model result in very small perturbations to the volcanic model velocity predictions. As slip deficit on the megathrust outboard of the volcano exerts the largest influence on the observed velocity field (Figure 4.1), we solved for slip deficit first and volcanic inflation 14 second. We subtracted out the predictions of the slip deficit model from the data, and then inverted the residuals for the volcanic source model as described below. We then iterated the solution by subtracting the volcanic inflationary signal from the original dataset and re-inverting for a revised slip deficit model. The residuals from this second slip deficit model were then inverted for an updated volcanic model. This process converged following the second inversion for slip deficit. The results for the final slip deficit model are shown in Figure 2, with residuals and the volcanic source model predictions shown in Figure 3. In supplemental table S4, we compare the volcanic source models for 2002-2005 and 2005-2017, assuming that the slip deficit model is time independent. 2.3.1 Megathrust Slip Deficit Inversion We estimated the slip deficit along the subduction interface (Figure 2.2) using the TDEFNODE software package (McCaffrey et al., 2002), following the same procedure as (Li and Freymueller, 2018). In particular, we define segment boundaries and assume that the slip deficit distribution is uniform along-strike within each segment, and apply smoothing in the downdip direction with an assumption that slip deficit is highest near the trench and cannot increase downdip. The previous study showed that this type of model with abrupt along-strike changes fit the data better than a model that was smoothed along strike without explicit segmentation. TDEFNODE models deformation within the overriding North American plate by combin- ing block rotation and the elastic effects of interseismic slip deficit. Faults are approximated as a series of planar dislocations in an elastic halfspace using the equations of (Okada, 1992). We defined block boundaries and fixed the angular velocities using the published Penin- sula (Li et al., 2016), Southern Alaska (Fletcher, 2002), and Bering (Cross and Freymueller, 2008) block boundaries and rotational poles as well as the Pacific and North American plate boundaries and rotational poles defined by (Argus et al., 2010). All sites in this study lie on the Peninsula block, and the relative motion of the Peninsula block and Bering Sea crust is small (Elliott and Freymueller, 2020a), so the only substantial elastic strains come from 15 the subduction zone. We used the Slab2.0 interface geometry (Hayes et al., 2018) for the Aleutian megathrust and created a grid of nodes at isodepths every 5 km down dip spaced 10 km apart along strike extending from the 10 km depth contour to 75 km. The final grid of model parameters has 73 nodes in the along-strike direction and 14 nodes in the downdip direction. We used the same smoothing methods as (Li and Freymueller, 2018) in order to fit a geologically reasonable locking distribution to the plate interface. In this approach, we divided the megathrust into a series of discrete segments within which slip deficit is held to be constant along strike within each segment and we allow down-dip smoothing parameters and segment boundaries to vary. Our inversion for the slip deficit parameters was performed on the horizontal components of the observed velocity field only, for the reasons discussed in the supplement of (Li and Freymueller, 2018). Vertical velocities have a small, systematic, long wavelength misfit that might be caused by inadequacies in the GIA models (collapse of the Laurentide ice sheet forebulge). 16 Figure 2.2 Best fit slip deficit model in this study. Warm colors represent a highly coupled interface, cool colors indicate creeping behavior. Aleutian trench shown as a thicker black line. Segments are numbered from right to left in order of decreasing slip deficit. Segment boundaries of (Li and Freymueller, 2018) plotted as dashed orange lines. The July 22, 2020 M 7.8 Simeonof earthquake centroid and GlobalCMT focal mechanism are shown, and aftershocks over the first 10 days are scaled by magnitude and shaded based on time after the mainshock. Location of Mt Veniaminof shown as red triangle. 0 nT magnetic anomaly contours from EMAG2 (Maus et al., 2009) plotted in the Pacific plate. Error bars give uncertainty ranges for boundary locations which are identical to (Li and Freymueller, 2018) except for the 3/4 boundary which has an uncertainty range of +30/-20 km (see figure S3), corresponding to 10% increases in overall misfit. 17 The four segment model of (Li and Freymueller, 2018) had all of the Shumagin Islands lying within the same segment, while having the Semidi segment terminate outboard of Veniaminof (Figure 4.1). Their model systematically underpredicted velocities on and around Veniaminof by up to 2-3 mm/yr (residuals to this model are shown in the inset of Figure 2.3) and had small but systematic misfits of velocities in the Shumagin Islands. We inserted a distinct segment outboard of Veniaminof and the eastern Shumagins (Segment 3 in Figure 2.2), and optimized the segment boundary locations in order to better fit these features of the data. Compared to (Li and Freymueller, 2018), we also used an updated subduction geometry (Slab 2.0), and a 10 km along-strike discretization of the interface rather than 25 km. However, it is the addition of the new data near Mt. Veniaminof that results in the largest change in the model. Aside from the additional data near Mt. Veniaminof used in this study, we used the same data as (Li and Freymueller, 2018) in the other regions, so we held the boundaries between segments 1/2 and 4/5 fixed in our inversions. The segment 1/2 boundary is well to the east of our focus area and was fixed to the same location as (Li and Freymueller, 2018). The data coverage immediately west of the Shumagin Islands is poor, so the exact location of the segment 4/5 boundary remains uncertain on the order of 10’s of km, although it must lie west of the islands. We set the boundary to lie directly in the middle of the gap in the observations. We then tested a variety of models in which the boundaries between segments 2/3 and 3/4 were shifted in 10 km increments in the along-strike direction, and selected the model with the lowest total misfit, shown in Figure 2.2. We found the location of the boundary between segments 2/3 to be weakly constrained due to a lack of near-trench data, but our best-fit location for this boundary is about 20 km east of the corresponding boundary in the model of (Li and Freymueller, 2018) (Figure 2.2). We find that this boundary must be located eastward of Veniaminof volcano, as the model for segment 2 would over-predict the observations there. The best-fit location for the boundary between segments 3/4 was found to be in the eastern Shumagin islands. Compared 18 to the model of (Li and Freymueller, 2018), our model finds lower slip deficit beneath the western Shumagins and higher slip deficit beneath the eastern Shumagins and the region directly east of the islands. Residuals for most sites are within their error ellipses (Figure 2.3). Trench-parallel residu- als on Kodiak Island may result from strike-slip faulting not considered in this model (Elliott and Freymueller, 2020a). Within our main study area, one of a cluster of sites in the area of Cold Bay has a velocity that is discrepant with the others, but there are no systematic residuals to our best-fit model. We tested the significance of the additional segment by comparing the total misfit of the 5-segment model with that of the 4-segment model of (Li and Freymueller, 2018). As the data and models are identical at the eastern and western ends, we considered the misfit of sites close enough to segments 3 and 4 of our final model to be influenced by the differences in models. We found that the misfit of the 5-segment model was 21% lower than the 4- segment model across that region, and that residuals to the 5-segment model were randomly distributed while those of the 4-segment model had systematic deviations. 19 Figure 2.3 Residual velocities for all sites used in this study after removal of predictions from the combined model. The inset box shows velocities for sites on Mt Veniaminof with the slip deficit model shown in Figure 2.2 removed as black vectors. Red inset vectors give predictions of the volcanic model in this study. Blue inset vectors give the velocity residuals between forward model predictions of (Li and Freymueller, 2018) and the observations. Note different velocity scale on inset. 2.3.2 Volcanic Source Inversion After removal of the best fit slip deficit model alone, there was a clear radial pattern in the residual velocities for the sites around Veniaminof (Figure 2.3, inset). We modeled the volcanic inflation source using a point source approximation (Mogi, 1958) in an elastic 20 halfspace. Due to sparse data we used the previously published horizontal source location (Fournier and Freymueller, 2008), which was located under the volcano’s broad summit, and solved for depth and source strength. In our inversion, depth was allowed to vary between 2 - 10 km and source strength between 1e+06 - 5e+06 m3/yr, and best fit parameters were obtained by a grid search. The forward model predictions from this model are shown in the inset of Figure 3. We repeated this process for the 2002-2005 data as well (Text S1). 2.4 Assessing Model Segment Boundaries Our preferred model contains five discrete along-strike segments, each with a unique down-dip profile of coupling (Figure 2.2). As in (Li and Freymueller, 2018), the general trend describes a plate interface that is fully locked in the east near Kodiak and transitions to fully creeping outboard of the western end of the Alaska peninsula. The width of the locked zone and amplitude of locking decrease westward along the subduction interface through a series of abrupt transitions. The main new feature of this model is the insertion of a distinct locking transition in the eastern Shumagin islands to explain the higher velocities of sites in the eastern part of the archipelago relative to the westernmost part. This change is required in order to better fit the data from the Veniaminof area and the eastern Shumagins. 2.4.1 A Slip Deficit Transition beneath the Shumagin Islands In order to assess why our 5-segment model fits the data better, we performed a posterior examination of the dataset in the Shumagin islands. We observe a distinct change in the trench normal elastic velocities of sites in the eastern Shumagin islands, which is evident when these are examined relative to continuous site AB07 located in Sand Point (Figure 2.4). The eastern Shumagin sites display trench normal velocities larger than those in the western Shumagins at the same distance from the trench (Figure 2.5). Figure 2.5 shows data and the model predictions computed at the center of each segment (away from segment boundary effects). The sites in the eastern Shumagins (magenta vectors in Figure 2.4 and magenta dots in Figure 2.5d) have trench-normal velocities that are intermediate between those predicted for the centers of segments 3 and 4. Our model fits these sites due to the 21 edge effect of the abrupt transition in slip deficit at the segment boundary. Figure 2.4 Elastic velocities of sites in the Shumagin islands relative to Sand Point corrected for far field block rotation. Heavy dashed line gives the location of the segment 3/4 locking transition shown in Figure 2.2. Magenta vectors give site locations in the Eastern Shumagins and correspond with the magenta dots in Figure 2.5. 2.4.2 Impact of Moving the Shumagin Segment Boundary Location Because the sites in the eastern Shumagin islands (and Outer Shumagins closest to the trench) are nearly atop the segment boundary and affected by the contributions of both segments, moving the segment 3/4 boundary from the best-fit location would change the slip deficit model for both segments. Here we assess the impact of those changes and the corresponding changes in model predictions. We use the test models described in the last section, and focus here on the boundary between segments 3 and 4. The best fit slip deficit 22 models are shown in panels a and b of Figure 2.5, while the corresponding velocity predictions are shown in panels c and d. Moving the segment 3/4 boundary eastward from the best-fit location results in a higher slip deficit estimate than the preferred model, for both segments 3 and 4 (red curves). When the boundary is shifted to the east, the eastern Shumagin sites become more effected by the segment 4 model, and this model must have higher slip deficit (and start to misfit the western Shumagin sites) to compensate. The segment 3 model also shifts to slightly higher slip deficit, as the impact of segment 3 on those sites is also lower. Moving the segment 3/4 boundary westward results in a lower slip deficit estimate for both segments (blue curves). In this case, the eastern Shumagins sites fall within segment 3 and the western Shumagin sites then lie atop the boundary, so the segment 3 model needs a lower slip deficit to fit these data. That results makes predictions that over-predict the eastern Shumagins and under-predict the Veniaminof data. In this case, with the edge effects from segment 4 contributing to the western Shumagin velocities, the estimated slip deficit on segment 4 also decreases. 2.5 Discussion The velocities of sites located on and around Mt Veniaminof fill in an important gap in observations, and this larger data set requires additional complexity in the megathrust slip deficit model. Specifically, our new model more effectively captures the gradient in observed velocities from the western to eastern Shumagin islands, and more accurately predicts the velocities in the region immediately east of the islands. In this section, we discuss the relationship between this revised slip deficit model and the 2020 earthquake, the causes and implications of segmentation, and for the seafloor-fabric hypothesis for the control on slip deficit. 2.5.1 The M 7.8 Simeonof Earthquake - Slip and Slip Deficit An M 7.8 thrust event ruptured the plate interface of the Alaska-Aleutian subduction zone on July 22 2020, initiating at about 20 km depth 50 km east of Simeonof island, the eastern island in the outer Shumagins. Rupture initiated close to our model’s segment 2/3 boundary, 23 Figure 2.5 Results of shifting boundary locations along strike, showing slip deficit in the form of locking fraction in the left panels (a and c), and predicted trench-normal velocity in the right panels (b and d) with the observations shown as dots. Each individual line corresponds to a model with unique locations of the boundary segments 3 and 4. Model segment boundaries were systematically shifted in increments of 10 km in the along strike direction in order to obtain best-fit locking profiles and the same set of models is shown in each panel. The best fit reference model (Figure 2.2) is denoted by thick black lines. Panels a and b are for the Veniaminof segment (segment 3) and panels c and d are for the Shumagin segment (segment 4). Lines are colored such that models with higher locking at the trench in segment 4 than the best fit model are plotted in red and lower values of locking at the trench in blue. Magenta dots denote sites in the eastern Shumagin islands. The dotted line in panel d shows the trench normal velocity predicted for the Veniaminof segment from the best fit model and is identical to the black line in b. 24 and propagated westward (Figure 2.2 and Figure S2), with the GlobalCMT centroid located just east of the islands. There are two prominent clusters of aftershocks, one east of the islands and the other west of the islands, with fewer events underneath the islands. (Crowell and Melgar, 2020) determined a slip model for the earthquake (see Supplementary Figure S2). Their primary slip patch corresponds to our model’s segment 3, and slip extends across our segment boundary at greater depth. However, there is a close correspondence between the edge of the main slip patch and our segment 3/4 boundary. The western aftershocks also appear to be at greater depth than the main rupture, and have continued with a different temporal pattern, making it possible that these represent small asperities in a region of triggered postseismic creep because they occur at the extreme depth limit of slip deficit (Figure 2.2). Our model predicts a moment deficit rate for segment 3 of 1.38e+19 Nm/yr. If we assume that the 1917 Shumagin earthquake was the last time that this portion of the plate interface ruptured, then the segment has over that time interval accumulated 1.42e+21 Nm of slip deficit, which is the equivalent of a Mw 8.1 earthquake if it were to rupture all at once. Within the rupture area of (Crowell and Melgar, 2020), the model predicts a cumulative slip deficit since 1917 equivalent to an Mw 7.65 event, slightly smaller than the Mw 7.8 GCMT determination, suggesting that the 2020 event released all available slip within it’s slip patch. The additional moment in the 2020 earthquake could reflect residual slip deficit left over in 1917, and the earthquake may have had modest additional slip outside of our segment 3 or at shallower depth. The rupture was confined mainly to the 30-50 km depth interval and did not propagate updip, where our model predicts higher slip deficit, although this is heavily influenced by our assumption that the slip deficit rate can only decrease in the downdip direction (thus the maximum must occur at the trench). In total, over the 10-30 km depth range (updip of the 2020 event) we predict enough cumulative slip deficit since 1917 to host a Mw 7.95 event, so regardless of the rupture beneath the Shumagins the updip section of the 25 megathrust may still accommodate significant levels of slip deficit, sufficient to host Mw ∼ 8.0 future events. This value is likely an upper bound, because slip deficit at the trench might be smaller and seafloor geodetic measurements would be needed to confirm our model assumption. However, in any case our model suggests that events much larger than Mw ∼ 8 are unlikely in this segment, consistent with the lack of geological evidence for such large events (Witter et al., 2014). We also note that the model of (Li and Freymueller, 2018) had a substantially lower estimates of an accumulated Mw of 7.43 in the rupture area and 7.68 in the updip region, totaling Mw 7.78 for the entire section in the time interval from 1917-2020. While such a value is close to the determined moment magnitude of the earthquake, it would imply that 100 years of accumulated slip deficit was entirely released by the earthquake and would be inconsistent with observed aftershock patterns. Further studies that incorporate multidisci- plinary datasets are necessary in order to better constrain rupture complexity arising from this event. 2.5.2 Nature of and Controls on Along-Strike Segmentation A geodetic slip deficit model provides an estimate of a kinematic property: whether or not slip is occurring on the megathrust. The plate interface may include regions that are mechanically locked (by friction), which will be surrounded by regions that are not slipping because of the stress shadowing effect of the mechanically locked regions (Hetland and Simons, 2010), (Takagi et al., 2019) Although most such studies have focused on updip and downdip stress shadows, (Herman et al., 2018) show (their Figure 3b) that a clear stress shadow is expected in the along-strike direction as well. In their model, a fully mechanically locked patch would be flanked by a gradient to 50 percent slip deficit over a distance of 20-30 km. This gradient might be more abrupt if the region is partially locked or the mechanically locked region is smaller. Our data are fit better by a model in which the along-strike boundary is even more abrupt than that. There are likely two implications of this. One is that the mechanically locked region of the segment is almost certainly smaller than the total 26 width of the segment, but the actual along-strike gradient is not imaged in the inversion because of along-strike smoothing within each segment. The second is that the along-strike width of the gradient in slip deficit between segments is likely not broader than shown in the model of Herman et al. (2018); further analysis of partially-locked segments and the interplay between smoothing and distributions of creep should be investigated. As with (Li and Freymueller, 2018), our model puts the transition between the fully- locked segment 1 and the partially-locked segment 2 at the boundary between N-S Farallon- Pacific and E-W Kula-Pacific seafloor fabric trends. The location of our segment 3 corre- sponds to the northernmost E-W anomalies of the Kula-Pacific trend. The width of segments 3 and 4 totals about 150 km, or roughly ten times the effective elastic thickness of the oceanic lithosphere inferred from the shape of the outer rise, but roughly comparable to the total width of the outer rise bulge (Garcia et al., 2019). Overall, our model shows a slightly better correspondence between the segment boundaries and the locations of the changes in seafloor fabric than did the previous model. 2.6 Conclusions We present an updated and optimized model for the along-strike variation in slip deficit on the megathrust near the Shumagin segment of the Alaska-Aleutian subduction zone. We show that a larger dataset incorporating observations from Mt. Veniaminof requires a distinct segment of higher slip deficit outboard of the volcano and the eastern Shumagin islands in order to better fit observed trends in the interseismic velocities. Velocities of sites in the eastern Shumagins are intemediate between those in the western Shumagins and those near Veniaminof, and are best explained by an abrupt along-strike boundary in the slip deficit that lies directly under the islands. This segment appears to have been the main rupture area in the July 22, 2020 M 7.8 Simeonof earthquake. We also show that longer period, decadal-scale inflationary signal at the volcano between the campaign surveys of 2005 and 2017 is driven by a deeper ∼9km point source, in agreement with InSAR results of (Lu and Dzurisin, 2014). The average volcanic deformation rate from 27 2005-2017 was substantially slower than the 2002-2005 inflation rate, although the lack of time resolution makes it difficult to constrain the pressurization history further. When taken together, our approach of modeling tectonic and volcanic deformation together through an iterative process allows us to arrive at better models for both. 28 UPDATES ON THE ALEUTIAN MEGATHRUST COUPLING MODEL CHAPTER 3 29 3.1 Introduction The Alaska-Aleutian megathrust outboard of the Alaska Peninsula has ruptured in a sequence of earthquakes that began with the Simeonof earthquake on July 22, 2020. The Simeonof earthquake was followed a year later by the July 28, 2021 M 8.2 Chignik earthquake and the 2023 M 7.2 July 16, 2023 Alaska Peninsula earthquake. At of the time of publication of Chapter 2 of this dissertation only the Simeonof event had yet occurred. The portion of the subduction interface that ruptured during the Simeonof earthquake had previously been identified as the ’Shumagin Seismic Gap’ (i.e. Davies et al., 1981), one of the few regions that had not hosted an instrumentally recorded M > 8.0 earthquake in the 20th century. Seismic gaps have been observed in subduction systems globally and are commonly thought to contain elevated risk of future great earthquakes (i.e. McCann et al., 1979). The 2021 M 8.2 Chignik earthquake, which occurred almost a year after the Simeonof earthquake, was the largest earthquake located within the Alaska-Aleutian subduction zone since 1965 and the 7th largest recorded earthquake in US history. Its USGS hypocenter was located around 50 km east of the Simeonof hypocenter, downdip of the zone of maximum coupling in segment 2 of the model presented in Chapter 2. The hypocenter for the 2023 Alaska Peninsula earthquake was located around 50 km west of the Simeonof hypocenter at the edge of the plate coupling model presented in Chapter 2. The USGS finite fault model for the 2023 Alaska Peninsula earthquake showed rupture that broadly co-located with elevated aftershock activity from the 2020 Simeonof earthquake. The close temporal relationship between these events suggests that they may be related to one another. The last time the plate interface outboard of the Alaska peninsula ruptured in a sequence of earthquakes began with the 1938 Mw 8.3 earthquake west of Kodiak. Tide gauge records suggest that the 1938 rupture could plausibly have been located within and updip of the 2021 Chignik rupture patch (Freymueller et al., 2021; Elliott et al., 2022). The 1938 event was the first in a remarkable sequence of 5 large (M > 8.0) events across the Alaska-Aleutian subduction zone that persisted until 1965 (i.e. Sykes, 1971), the largest of which was the 30 1964 M 9.2 great Alaskan earthquake. The spatial extent of aftershock activity from these 5 events alone covered almost the entire length of the subduction zone from Anchorage to Attu, a distance of over 3500 km (Tape and Lomax, 2022). Given this prior rupture history it is entirely plausible that the rupture sequence touched off by the 2020 Simeonof event is ongoing and may include future large (M > 8.0) earthquakes. The rupture patches for the 2020-2023 sequence of temporally clustered events did not overlap with each other in space. Calculations of Coulomb failure criteria along the plate interface suggest that it is likely that failure occurred in the latter 2021 Chignik and 2023 Alaska Peninsula events due to stress changes that were induced by previous earthquakes in the sequence (i.e. Xiao et al., 2021; Elliott et al., 2022). Due to both the Simeonof/ Chignik events as well as a wealth of new ocean bottom seismic data recorded during the AACSE (Alaska-Aleutian Community Seismic Experiment) deployment of ocean-bottom seismometers off the coast of the Alaska peninsula a number of additional opportunities related to exploring new variations of the trench coupling model have presented themselves in the years since. This chapter serves as a follow-on to Chapter 2 in order to elaborate on additional work on co-authored publications aimed at exploring plausible variations of trench coupling models. 3.2 Data & Methods The trench coupling models presented in this chapter were obtained using mostly the same set of GPS site velocities presented in Chapter 2. The first 5 segment model was fit to the exact same dataset, whereas the latter 6 segment model was fit to an extended dataset that includes additional sites on Kodiak that were not considered in the inversion presented in Chapter 2 because they were outside of the model space. Those site velocities were taken from (Elliott and Freymueller, 2020b) and were time windowed in order to reflect interseismic behavior. The choice around site velocities was made because of the 2020-2023 rupture sequence, which has introduced both coseismic offsets and postseismic transients to post- July 2020 GPS timeseries. Since trench coupling is fundamentally an interseismic 31 process we avoid contamination of model solutions with non-interseismic effects by excluding all GPS data from the period after July 2020 from our input dataset. It is outside the scope of our studies to fit postseismic curves to GPS observations, which would be premature given only 3-4 years of observations since the Simeonof event and 2-3 years since the Chignik event. Furthermore, it would be difficult to disentangle the two postseismic trends from each other. We use the TDEFNODE software package to model deformation within the overriding North American plate by combining block rotation and the elastic effects of interseismic slip deficit. Faults are approximated as a series of planar dislocations in an elastic halfspace using the equations of (Okada, 1992). We defined block boundaries and fixed the angular velocities using the published Peninsula (Li et al., 2016), Southern Alaska (Fletcher, 2002), and Bering (Cross and Freymueller, 2008) block boundaries and rotational poles as well as the Pacific and North American plate boundaries and rotational poles defined by (Argus et al., 2010). All sites in this study lie on the Peninsula block, and the relative motion of the Peninsula block and Bering Sea crust is small (Elliott and Freymueller, 2020a), so the only substantial elastic strains come from the subduction zone. We used the Slab2.0 interface geometry (Hayes et al., 2018) for the Aleutian megathrust and created a grid of nodes at isodepths every 5 km down dip spaced 10 km apart along strike extending from the 10 km depth contour to 75 km. The final grid of model parameters has 73 nodes in the along-strike direction and 14 nodes in the downdip direction. We use the same segment boundaries as presented in Chapter 2 of this dissertation, and discuss below choices about the down-dip smoothing parameters. A key assumption of the model detailed in Chapter 2 of this dissertation is that slip deficit is held to monotonically decrease down-dip within each segment of the subduction interface. This assumes that maximum slip deficit will be present in the shallowest portion of the megathrust closest to the trench. In practice most subduction earthquakes nucleate netween 15-30 km depth and it is relatively rare for megathrust earthquakes to rupture updip to the trench due to velocity-strengthening slip properties for the shallowest portion of the 32 subduction interface between 0-15 km depth (i.e. Lay et al., 2012; Marone and Scholz, 1988). Our modeling approach is useful in that it enables us to solve for along strike variations in trench coupling properties while providing a maximum possible estimate within each segment. In the absence of near-trench GPS velocity observations in the Shumagin and Veniaminof segments our ability to distinguish between such models is limited. The rupture of the July 22, 2020 M 7.8 Simeonof earthquake directly challenged the downdip-decreasing assumption about coupling distribution along the trench interface. Peak slip was located between 30-40 km depth, with slip extending more broadly along strike than updip (Figure 3.1). Shallow rupture from the Simeonof earthquake was highly limited (DeSanto et al., 2023), and thus a reconsideration of the assumptions made about how slip deficit is distributed downdip along the subduction interface was warranted. The goal of the first study presented in this subchapter was to fit an alternative coupling distribution following a Gaussian curve in the downdip direction. Initial results for the slip patch of the Simeonof earthquake based on InSAR interferograms, static GPS displacements, and seismic waveforms suggested that significant slip was located east of the Shumagin islands but downdip of the 20 km depth contour of the Aleutian trench (Figure 3.1). In light of this, we ran some tests for alternate coupling models that include gaussian down-dip coupling distributions within each segment. Gaussian fits to the data were determined in an iterative least squares & simulated annealing inversion and the model is shown in Figure 3.1 panel b. 3.3 Results 3.3.1 Gaussian 5 Segment Model The gaussian model more closely tracks rupture in the 2021 Simeonof earthquake when compared to the original model presented in Chapter 2 (Figure 3.1). The main rupture patch beneath the eastern Shumagin islands lies directly downdip of the peak coupling of the middle segment, and the lateral extent of rupture into the two neighboring segments generally tracks the peaks in partial coupling. The cumulative slip deficit in each segment of the gaussian 5 segment model is substantially lower compared to the original. Every 33 Figure 3.1 From (Xiao et al., 2021) slip model for the 2020 Simeonof earthquake overlain on alternative coupling distribution models. segment west of the fully coupled kodiak section has a lower peak slip deficit and a more smooth distribution about the peak. For example, the segment underlying the western Shumagin islands (segment 1 on Figure 3.1) has a peak coupling ratio of 0.6 in the original model and a peak coupling ratio of 0.4 in the gaussian model. This difference has important implications for the seismogenic potential of the subduction zone, suggesting that seismic moment sufficient to rupture in a future Mw > 8.0 earthquake was not necessarily building up in the shallowest portion of the subduction interface updip of the Simeonof rupture. Another notable feature of the Gaussian coupling model was the partially locked patch 34 Figure 3.2 Slip models for the 2020 Simeonof, 2021 Chignik, and 2023 Alaska Peninsula earthquakes overlain on the coupling model from (Xiao et al., 2021). Purple contours give 1 m slip contours for the Chignik and Simeonof events. Orange contours give 0.5 m slip contours for the 2023 Mw7.2 Alaska Peninsula earthquake. Cyan dashed line gives the rupture area of the 1964 Mw9.2 Alaska earthquake from (Davies et al., 1981). Coupling distribution is colored such that solid black represents full plate interface coupling while white represents full creeping behavior. in the segment directly east of the Simeonof rupture. This had a peak locking fraction of ∼ 0.6 and turned out to be the locus of the 2021 Mw 8.2 Chignik event (Elliott et al., 2022) (Figure 3.2). Slip modeling work suggests that this segment was almost entirely ruptured in this event. The 2023 Mw7.2 Alaska peninsula earthquake ruptured the easternmost extent of the coupling model. Whether the Alaska Peninsula earthquake represents that last in this multi-year sequence of earthquakes remains to be seen. 35 3.3.2 Gaussian 6 Segment Model The AACSE ocean-bottom seismometer array represented the most dense deployment of seismic instruments to date, offering unprecedented opportunities to model seismic ve- locity structure in the shallow portion of the Alaska-Aleutian megathrust. Work that was performed modeling Vp/Vs variations along the subduction interface by (Wang et al., 2022) suggested that a substantial heterogeneity in fluid content exists offshore of the western edge of the Kodiak archipelago. As part of this work, we tested a new coupling model that al- lowed for an extra variation of slip deficit beneath Kodiak island separate from the segment outboard of Chirikof island. This slip deficit model is shown in Figure 3.2. All other coupling boundaries to the east of the model space were held to the same as the other previous models and a gaussian downdip coupling distribution was used. We find that the geodetic data prefers a highly coupled shallow section of the subduction interface in the Chirikof segment, while the Kodiak segment does not necessitate full shallow coupling. Peak slip deficit ratio in the Kodiak segment reaches a value of ∼ 0.9 at around the 15 km depth contour and decays to 0 by about the 40 km depth contour. The coupling fraction at the trench is close to 0.7, however this is prone to a high degree of uncertainty due to the lack of observations outboard of Kodiak island. The Kodiak segment in the model almost perfectly aligns with the western extent of the 1964 M 9.2 Great Alaskan earthquake, so it is entirely plausible that this is the same asperity that ruptured in the 1964 earthquake. Interseismic coupling in Chirikof segment directly west of the Kodiak segment is a robust feature that is virtually identical between the 6 segment and 5 segment models presented in this chapter, as well as the downdip-decreasing model presented in Chapter 2 of this dissertation. 36 CHAPTER 4 NEW INSIGHTS INTO THE ACTIVE TECTONICS OF THE NORTHERN CANADIAN CORDILLERA FROM AN ENHANCED EARTHQUAKE CATALOG 37 ABSTRACT Seismic activity in the Northern Canadian Cordillera is characterized by diffuse earthquakes that extend hundreds of km northwest from the Yakutat collision zone. We use 25 months of broadband seismic data from Mackenzie Mountain Earthscope Project (MMEP), USAr- ray Transportable Array (TA), and permanent Canadian National Seismic Network stations to present a local earthquake catalog with high sensitivity to small regional events. Deep learning techniques are adopted for both seismic phase detection and association. Event relocations are performed to provide well constrained estimates of earthquake depth distri- butions. Clusters of seismicity spanning the upper crust are located in the central Richardson Mountains, along the Tintina fault, and in the northeast Selwyn Basin. These clusters sug- gest that the core of the Richardson Anticlinorium is tectonically active and that the Tintina fault is a locus for low levels of active deformation. We interpret seismicity in the northeast Selwyn Basin as primarily occurring in the hanging wall of the Plateau thrust fault and suggest that some combination of localized duplex structures and lithological strength con- trasts both within the Selwyn Basin and between abutting Paleozoic shelf sequences may be responsible for seismicity in the Mackenzie Mountain foreland. Plain Language Summary The Northern Canadian Cordillera is situated inboard of a transition from a right-lateral tectonic regime in coastal southeast Alaska to a subduction setting beneath southern Alaska. It presents a unique case study to examine the length scales over which tectonic deformation is able to permeate into a continental interior. In this study we use data from a relatively dense regional network of seismometers in order to present the most spatially complete catalog of earthquakes in the Northern Canadian Cordillera to date. We use machine learning techniques for the detection and classification of small magnitude earthquakes and use our results to analyze regional tectonics. Active faults that host earthquakes are identified and we use the earthquake depth distributions to further distinguish between mechanisms for regional stress transfer. We find evidence that the Tintina fault is presently active along a 38 confined area at a low rate. We also find that low magnitude earthquakes do not generally occur within the Mackenzie Mountain Fold-Thrust belt but are prevalent within the foreland Selwyn Basin. We also detect a pattern of relatively deep (∼ 30 km) earthquakes in the northern Richardson Mountains. 4.1 Introduction The Northern Canadian Cordillera (NCC) marks the northeast extent of the North Amer- ican Cordillera orogenic belt, a chain of mountains and active tectonics along the western margin of North America that stretches from Arctic Alaska to Panama. The fundamental driver of deformation throughout the Cordillera is the interaction of the Pacific and North American tectonic plates. In the northern Gulf of Alaska the tectonic regime is characterized by right-lateral strike-slip motion north of Cascadia that transitions to the northwest into subduction beneath Alaska and the Aleutian arc. The Northern Canadian Cordillera sits inboard and northeast of this transition and transpressive tectonic stress has been inferred to extend from the Yakutat collision zone in the northeast Gulf of Alaska for 100’s of km inboard into the North American continental interior. The highest levels of historic seismicity are observed seaward of the Denali fault, how- ever more diffuse patterns of M > 5.0 historic earthquakes extend into the Mackenzie and Richardson Mountains in Yukon and Northwest Territories (Figure 4.1). The permanent seismic network in this region has always been sparse, which has resulted in large location uncertainties and a limited ability to detect smaller events. For the region northeast of the Denali fault, previous studies have not attempted to connect active zones of seismicity to specific active faults and the sparse spatial distribution of recorded events has left it un- clear which structures remain active. Across the central part of the study region near the Tintina fault, only a few larger earthquakes were recorded previously and with variable focal mechanisms (Figure 4.1). The hot backarc mantle which underlies the Northern Canadian Cordillera makes it a unique case study into the length scales and mechanisms by which forces from plate boundary 39 Figure 4.1 Seismicity and tectonic setting of the Northern Canadian Cordillera. All earth- quake hypocenters for events M>3.5 between 1985 and 2022 plotted in blue. Focal mech- anisms shown for events from 1976-2023. Blue focal mechanisms for events M>5.0 be- tween 1976-1995 are from the GCMT catalog. Black focal mechanisms for events M>4.0 from 1995-2023 are from the Pacific Geoscience Centre (PGC) CMT catalog (i.e. Kao et al., 2012). GOA=Gulf of Alaska, YB=Yakutat Block, CDF=Cordilleran Deformation Front, RM=Richardson Mountains, WM=Wernecke Mountains, DRF=Duke River Fault, LTZ=Liard Transfer Zone, GBL=Great Bear Lake. 40 zones are able to permeate into continental interiors. Backarc settings are a fundamental feature of subduction systems globally, yet aspects of their evolution remain uncertain (i.e. Clark et al., 2008). In this study we use 25 months of broadband seismic recordings to present an earthquake catalog for the NCC with a high sensitivity to small magnitude earthquakes. This time interval was chosen because it covers the entire deployment of the Mackenzie Mountain Earthscope Project (MMEP) PASSCAL deployment (i.e. Baker et al., 2020), which coupled with the Earthscope TA deployment in northwest Canada represents the densest available set of regional broadband seismic recordings. We also perform double-difference event relocation to recover the most precise 3D distribution maps of seismicity to date. We analyze this catalog to identify the current active tectonic structures and evaluate the implications for stress transfer models. 4.1.1 Previous Work Earthquake focal mechanisms (Figure 4.1) suggest thrusting behavior in the Mackenzie Mountains, transitioning to the northwest into right-lateral shearing in the Richardson Moun- tains (i.e. Hyndman et al., 2005; Ristau et al., 2007). Estimates of shortening rates based on recorded seismicity rates are on the order of 1-10 mm/yr in the Mackenzie and Richardson Mountains and < 1 mm/yr within the Selwyn Basin and along the Tintina Fault (Leonard et al., 2008). The regional earthquake patterns were explained by (Mazzotti and Hyndman, 2002) using an orogenic float model (i.e. Oldow et al., 1990) in which tectonic strain from the Yakutat collision zone is transferred at depth along a weak lower crustal detachment into the Mackenzie and Richardson Mountains. In this scenario the North American craton acts as a rigid backstop and forces the localization of strain along inherited fold-thrust structures within the relatively weak crust in the abutting Mackenzie Mountains. Key to this model are observations of high regional heat flow, which suggests significant lithospheric strength only in the upper 15 - 20 km of Cordilleran crust (i.e. Lewis et al., 2003; Mazzotti and Hyndman, 2002). An alternative though not necessarily contradictory model proposed by (Finzel et al., 2015) suggests that significant regional stress accumulates throughout north- 41 ern Alaska and the Canadian Cordillera due to mantle traction forces acting at the base of the lithosphere. These traction forces are thought to be the result of global mantle flow that is oriented north-south through Arctic Canada and Alaska and is deflected by subduct- ing Yakutat lithosphere beneath Alaska. (McConeghy et al., 2022) presented a geodynamic model for surface displacements based on mantle tractions that showed that uplift in the Mackenzie Mountains does not require strain transfer along a lower crustal detachment, and that instead the flow of mantle material away from the Yakutat collision zone possibly in- teracting with the north-south oriented flow originating from the Arctic alone is enough to produce thrust earthquakes in the Mackenzie Mountains. Body wave tomography work by (Est`eve et al., 2020) shows two thick high-velocity anomalies inferred to be cratonic fragments in the asthenosphere to the northwest and southeast of the Mackenzie Mountain fold-thrust belt, and that the belt itself is under- lain by low-velocity anomaly in the asthenosphere. This result suggests that mantle flow may indeed be channeled along the base of the lithosphere underlying the Mackenzie Moun- tains, which is observed to be relatively warm and thin (i.e. Lewis et al., 2003; Audet et al., 2020). Shear wave splits at MMEP stations from (Bolton et al., 2022) show an allignment of asthenospheric fabric oriented southwest-northeast in the Mackenzie Mountains, subparallel to the direction of North American absolute plate motion. This is inferred to be indicative of asthenospheric flow beneath the Mackenzie Mountains that is largely coupled to North America. A rotation of splitting directions is observed towards the southwest in the vicin- ity of the Tintina fault, which does not readily align with models of regional mantle flow and may instead represent fossilized lithospheric fabrics related to displacements along the Tintina fault that may have acted as conduits for asthenospheric flow related to the opening of a slab window (i.e. Thorkelson and Taylor, 1989) in the early Miocene (∼ 20 Ma). 42 Figure 4.2 Paleozoic depositional features from (Cecile et al., 1997) plotted on top of modern terrane boundaries from (Colpron et al., 2007) NE of the Tintina Fault. Basinal strata of the Selwyn Basin and Richardson Trough colored in grey. Thick dashed line gives the axis of the Paleozoic Niddery High after (Cecile and Norford, 2000). AA = Anvil Allochthon, MCE = Misty Creek Embayment, CDF = Cordilleran Deformation Front. 4.2 Regional Geology 4.2.1 Paleozoic Evolution The Cordilleran orogen in the Canadian Yukon & Northwest Territories initially formed due to incipient seafloor spreading during the late Proterozoic breakup of the Rodinian su- percontinent (∼ 750 Ma). This process created the ancient continent of Laurentia, which was the ancestor to modern day North America. Rifting gradually progressed offshore of the margin of Laurentia throughout Paleozoic time (Bond et al., 1985). Marine platformal to basinal sediments deposited during this interval reflect a tectonically stable environment likened to modern eastern North America (i.e. Ross, 1991; Monger and Price, 2002). These sequences from the peri-Laurentian realm are today located in the most inboard regions 43 of the modern NCC, northeast of the Tintina fault. The early Paleozoic Mackenzie Plat- form (Figure 4.2) is composed of shallow marine miogeoclinal sequences containing platform carbonates and mature clastic sediments that thin northeast towards the Canadian Shield (Gabrielse and Yorath, 1989). To the southwest the Mackenzie platform grades into the Selwyn Basin (Figure 4.2), which is composed of late Neoproterozoic to late Devonian deep water argillaceous shales, cherts, and limestones (i.e. Cecile et al., 1982; Mair et al., 2006). The Paleozoic environment outboard of Laurentia in which the Mackenzie Platform and Sel- wyn Basin were deposited was a marginal to deep marine setting with a rather narrow hinge line ∼ 20 − 40 km wide between shelf and basin (Cecile et al., 1982). The Misty Creek Embayment is a 100x150 km northwest-trending depositional feature that forms the northeast extension of the Selwyn Basin (Figure 4.2). It is interpreted as a distinct rift basin that developed from late Early Cambrian until Early Silurian time and is comprised of interbedded quartzites and siltstones in the northeast near the Mackenzie Platform that grade southwest into shales and siltstones (Fritz et al., 1991). The Niddery High forms the southwest boundary of the Misty Creek Embayment and is interpreted as a local submarine high with anomalously thin basinal strata (Cecile and Norford, 2000). By the early-mid Devonian the Misty Creek Embayment had been re-established as a marine platform and mid-Devonian to Carboniferous basinal deposition was confined to the region southwest of the Niddery High (Cecile and Norford, 2000). The Meilleur River Embayment forms the southeast extension of the Selwyn basin, located in the Nahanni region of the southern Mackenzie Mountains where it forms a ∼ 200 x 100 km convex-eastward protrusion into the Mackenzie Platform (Figure 4.2). The Meilleur River Embayment underwent basinal deposition from early Ordivician to mid-Devonian time with a brief interlude of platformal sedimentation dated to the middle Ordivician (Cecile et al., 1997). In the Wernecke Moun- tains a band of shelf sequences called the Ogilvie Platform separates the Selwyn basin in the south from the Richardson Trough to the north (Figure 4.2). Stratigraphic relationships suggest that both the Richardson Trough and Ogilvie Platform were originally part of the 44 Mackenzie Platform but became separate elements by the mid-Cambrian as a result of out- board rifting processes (Fritz et al., 1991). The Richardson trough has similar lithology to the Selwyn Basin, and may have developed as an aulacogen during mid-Cambrian extension outboard of Laurentia (Gabrielse and Yorath, 1989). 4.2.2 Structural Evolution Major structures in the Mackenzie, Wernecke, and Selwyn mountains reflect Mid-Mesozoic contraction and a subsequent transition to Late Mesozoic-Cenozoic transcurrent motion. Thrust faulting in the Selwyn Basin originated in the mid-Jurassic and is largely coeval with the emplacement of the outboard Yukon-Tanana terrane (Mair et al., 2006). An inversion of regional gravity data performed by (Hayward, 2019) revealed the presence of a regional-scale decollement beneath this inboard region that is inferred to have been active from the Mid- Cretaceous to the Paleocene, largely pre-dating the maturation of the Tintina fault. This decollement is likely to have facilitated large-scale stress transfer into the Selwyn Basin, leading to the development of ramp-like thrust faults through reactivation of rift-related basin-bounding normal faults. The Plateau Thrust fault outcrops in the central Mackenzie Mountains and juxtaposes Mesoproterozoic strata in its hanging wall with middle Devonian units in its foot wall (Cecile and Cook, 1981). It is inferred to have accommodated ∼ 30 km of shortening (Cecile and Cook, 1981) and is projected to extend at a low angle from its outcrop location in the Central Mackenzie Mountains southwest into the Selwyn Basin near MacMillan Pass, where it is presumably connected to the regional decollement at depth. In the northwest Yukon Territory, the Dawson Thrust is interpreted as a reactivated Neopro- terozoic basin-bounding fault that, along with with the overlying Tombstone and Robert Service Thrusts has accommodated around 100 km of shortening since the mid-Cretaceous (Mair et al., 2006). These structures correlate with a shallowing of the decollement identified by (Hayward, 2019). Beginning in the Late Mesozoic, the direction of convergence outboard of the Cordilleran margin began to rotate from dominantly east-west to north-south (i.e. Engebretson et al., 45 1985) in orientation, leading to the development and maturation of lithospheric-scale tran- scurrent faults. These first-order tectonic structures have cumulatively produced hundreds of km of dextral (right-lateral) displacements since the Mid-Cretaceous (i.e. Gabrielse et al., 2006). The two largest such features are the Tintina and Denali faults (Figure 4.1). The Tintina fault strikes northwest-southeast and runs from the Liard zone of northwest British Columbia for over 1200 km before terminating in the Yukon Flats of central Alaska and records displacements as far back as the late Early Cretaceous (110 Ma). Late Cretaceous offsets along the Tintina fault are geometrically complex and linked to widespread granitic plutonism (Gabrielse et al., 2006). Offsets along the McQuesten plutonic suite (ca. ∼67 Ma) suggest that 400 - 430 km of displacement along the Tintina fault has occurred throughout Cenozoic time (i.e. Murphy and Mortensen, 2003), the majority of which is inferred to have occurred in the Eocene (65 - 40 Ma) (Gabrielse et al., 2006). In northwest Canada, the Tintina fault generally forms the boundary between autochthonous sequences deposited along the ancient Laurentian margin and allochthonous terranes that have been accreted to the Cordilleran margin since the early Mesozoic. The Teslin fault strikes subparallel to the Tintina fault (Figure 4.1) and is located 50-150 km outboard of it in Yukon Territory. The Teslin fault is associated with dextral displacements of up to 130 km in British Columbia (Gabrielse, 1985) and extends northwest for 500 km, terminating in the north near Carmacks in Yukon Territory. Compared to the Tintina fault, the Teslin fault has been interpreted to be a relatively thin-skinned structure confined to the upper crust (Cook et al., 2004; Snyder et al., 2005). The Denali fault trends subparallel to the Tintina fault and is located outboard of it by 250 - 300 km (Figure 4.1). The Eastern Denali fault strikes northwest from Haines, Alaska where it is presumably connected to the strike-slip tectonics of the Queen Charlotte margin via the (now inactive) Coast Shear Zone and Chatham Strait faults (Brothers et al., 2018). It continues up through Yukon Territory to a linkage with the Totschunda fault near Mentasta Lake in the Alaska Range some 50 km southwest of Tok, Alaska. Complex modern 46 deformation histories are recorded throughout this region both along the Denali fault itself and along several subsidiary structures, most notably the Totschunda and Duke River faults (Matmon et al., 2006). The Totschunda and Denali Faults are highly interconnected, as evidenced by the rupture of the 2002 Denali fault earthquake, which initiated in the central Alaska range and propagated eastward before jumping over to the Totschunda fault east of the Denali/Totschunda junction. To the southeast, the Totschunda fault is presumed to be connected to the Fairweather fault via the Connector fault, which runs through the southwest corner of Yukon Territory. Tectonic strain is accommodated along all of these structures, with the Eastern Denali fault accommodating ∼ 8 mm/yr of slip and the Totschunda and (presumably the Connector fault as well) each accommodating ∼ 6 mm/yr of slip. West of the confluence of the Totschunda and Denali faults the slip rate generally decays from ∼ 12 mm/yr west of Mentasta lake to ∼ 9 mm/yr in the central Alaska range (Matmon et al., 2006). 4.2.3 Modern Tectonic Setting The Mackenzie Mountain Fold-Thrust belt is situated mostly in Northwest Territories, and in map view forms a ∼ 300 km convex-northeast protrusion into the North American continent bounded by the Richardson and Wernecke Mountains in the northwest, the Selwyn Mountains to the west, and the Liard transfer zone in the southeast. Thrust structures found throughout the belt are the result of Mesozoic contraction, and many were likely basin-bounding normal faults that were subsequently reactivated during accretion of the Intermontane terranes (Dusel-Bacon et al., 2002). The Nahanni region in the southeast Mackenzie Mountains is a zone of anomalously high rates of seismicity compared to the rest of the belt. This includes the 1985 M 6.6 and 6.9 Nahanni earthquakes, which are the two largest earthquakes in the instrumental record located inboard of the Tintina fault. The Selwyn Basin is located southeast of the Mackenzie Mountain foreland and is bounded to the southwest by the Tintina fault (Figure 4.2). Sedimentary packages within the basin are comprised of clay-rich shales, cherts, and limestones that are inferred to have deposited 47 throughout the Paleozoic in deep water marine basins outboard of the Laurentian margin (i.e. Fritz et al., 1991). Following deposition, these sequences were strongly deformed and thrust faulted during Mesozoic contraction. Prior studies have suggested that the Selwyn Basin is seismically inactive, with virtually no M > 3.5 earthquakes being located within a zone ∼ 150 km wide extending northeast from the Tintina fault (Figure 4.1). The Richardson Mountains form a north-south trending band of topographic relief along the northern segment of the Yukon/ Northwest Territories border (Figure 4.1). Stratigraphic units in the core of the Richardson mountains comprise the Richardson Anticlinorium, a rather narrow (∼ 50 km) north-plunging structural anticline interpreted by (Hall and Cook, 1998) as a ’pop-up’ structure that may have formed over a ramp-like lower crustal detachment extending at depth west beneath Eagle Plain and has accommodated ∼ 30 km of cumulative convergence. A positive Bouguer gravity anomaly is centered in the Richardson Anticlino- rium (Pinet, 2021a) and the sequence of mostly north-south trending faults that cross-cut it is referred to as the Richardson Array. Previously reported regional seismicity catalogs (i.e. Leonard et al., 2008) record high concentrations of earthquakes in the Richardson Mountains and focal mechanisms suggest dominantly dextral strike-slip deformation (Hyndman et al., 2005b). However, (Leonard et al., 2008) did not attempt to relate these events to specific faults, because of the large location uncertainties that resulted from the sparse network. Outstanding questions remain with regards to the relationship between seismicity in the Mackenzie and Richardson Mountains. (Leonard et al., 2008) calculated a strike-slip rate of 2.1 mm/yr in the Richardson Mountains and a thrusting rate of 1.8 mm/yr in the Mackenzie Mountains, yet the two are linked by a significantly higher zone of transpressive deformation in the eastern Wernecke Mountains with an inferred shortening rate of 8.2 mm/yr. The direction of convergence in this zone, which is broadly co-located within the northern Ogilvie Platform (Figure 4.2) is intermediary between the dextral north-south oriented strike slip to the north in the Richardson Mountains and thrusting to the southeast in the Mackenzie Mountains. (Leonard et al., 2008) attribute the discrepancy in magnitude of shortening 48 to assumptions regarding earthquake return periods made in their calculations. Geodetic strain transfer into the NCC was detected by (Leonard et al., 2007) using a relatively sparse network of campaign and continuous GPS stations. Their results show strain transfer as far as Whitehorse (∼ 150 km outboard of the Tintina fault) but not further inboard into the Mackenzie and Richardson Mountains. (Fl¨uck et al., 2003) presented a regional analysis of effective elastic thickness (Te), which varies from 15-20 km in the central NCC to ∼ 30 km in the northeast Selwyn Basin and > 30 km in the Richardson Mountains. These estimates of Te were derived by computing the coherency between topography and Bouguer gravity anomalies, which can be problematic in regions with ongoing orogenic activity and thus the southeastern portion near the Denali fault was windowed out of their analysis. 4.3 Data and Methods We derived the seismicity catalog from broadband seismic recordings of the Macken- zie Mountain EarthScope Project (MMEP) temporary deployment of broadband seismome- ters as well as from the EarthScope TA and permanent CNSN (Canadian National Seismic Network) stations. A flowchart describing the analysis steps is shown in Figure 4.3. We performed seismic phase picking and association using recent software packages exploiting machine-learning techniques, and events were relocated using a 3D velocity model. We use data from the stations shown in Figure 4.4 for the time period between August 1, 2016 through to September 1, 2018, when all MMEP stations were active. The MMEP project deployed a dense line of 41 broadband seismometers (Baker et al., 2020) that extended from northwest British Columbia for over 1000 km northeast to Great Bear Lake in the Northwest Territories. These data are used alongside recordings from the Canadian National Seismic Network and the EarthScope TA footprint in northwest Canada. The combination of these seismic networks represents the most dense regional sampling of seismic instruments to date. We then relocated all events to produce more accurate relative locations. 49 Figure 4.3 Workflow for the methods used in this study. 4.3.1 Phase Detection We adopted a hybrid approach to the detection of seismic phase arrivals in order to leverage the capabilities of the EQTransformer phase detection package (Mousavi et al., 2020) while also being cognizant of its current limitations. EQTransformer is a deep neural network that reads in seismic waveforms and generates 3 different prediction curves that give the likelihood of a detection, P-, and/ or S-arrival being present within the data. The detection probability is a continuous curve while the P- and S- detections are discrete points in time that exceed a user-defined detection threshold and are assigned probabilities based on the characteristics of the waveform. The algorithm operates on each waveform trace individually in order to generate phase detections. We chose to use EQTransformer due to its low false detection rate and relative high computation speed when compared with other Deep Learning-based pickers (i.e. Garc´ıa et al., 2022; M¨unchmeyer et al., 2022). 50 In total, EQTransformer returned 1,103,599 phase detections for all stations over the 25 month interval. EQTransformer was highly effective at picking P waves, but picked relatively few S waves. From our examination of event records, we found picking problems for events as small as magnitude 3.3 at event-receiver spacings greater than 200 km. Erroneous picks were due to EQTransformer’s mischaracterization of head waves, in particular Pn phase arrivals, as direct Pg phases. As it is currently distributed, EQTransformer only outputs prediction curves for direct P (Pg) and S (Sg) phases. Since the method uses a built-in ‘attention mechanism’ that looks for S arrivals directly following a P arrival, a cascading effect of the mischaracterization of Pn phase arrivals as Pg phases was that the ‘true’ Pg phase arrival often was confused for an Sg arrival. For larger events, we found that EQTransformer also sometimes identified a part of the wavetrain between the true P and S arrivals as an additional Pg phase, often paired with a following Sg phase. These mis-identified phases resulted in inaccurate associations and spurious events with unacceptably large travel time residuals. Thus we needed a way to eliminate these erroneous picks while retaining the excellent performance of EQTransformer in locating smaller events. It was not possible automatically to window out all stations further than 200 km from event hypocenters, because we could not know the hypocenter location until after events were located. Instead, we excluded all picks generated by EQTransformer that occurred within 5 minutes after an event with magnitude greater than ML 3.3 in the Canadian National Earthquake catalog. We also address this problem through our association procedure, as detailed below. This magnitude threshold is substantially higher than the CNSN magnitude of completeness (2.9), so we assumed that all events in the range that we excluded were reliably recorded. This screening eliminated the problem of erroneous picks and associations, however it is possible that a few early aftershocks after every large event could be missed. Since we adopt all CNSN events M > 3.3 into our catalog regardless of their time relationship to other events, the number of events lost due to this time windowing is minor. For the events exceeding the M 3.3 threshold, we used a mixture of phase picks that were generated 51 by analysts at the Pacific Geoscience Center for all CNSN sites, and arrivals we picked by hand for the sites along the MMEP transect. In total, we characterized 274 events using CNSN and manual picks. 4.3.2 Phase Association and Event Relocation Phase association is the procedure by which a timeseries of earthquake phase detections across a receiver network are linked into common events. We use the PhaseLink software package (Ross et al., 2019) in order to perform this task for the EQTransformer pick database that we generated. PhaseLink is a deep learning approach for associating seismic phase arrivals. It works by generating a series of synthetic earthquakes and resulting moveout patterns for phase arrivals across a user-defined seismic network. This is performed across a regular grid that encloses the study area and is constituted of a user-defined 1-D velocity profile. A neural network is then trained to associate these picks together into their ’true’ corresponding earthquakes. Since each model trained with PhaseLink is unique to a par- ticular study region and receiver geometry that corresponds with the desired use case, it is necessary for the user to train models themselves. Part of our approach to windowing out mis-identified head wave phases involved training the associator models on synthetic data generated for event-receiver distances up to 200 km. This almost entirely eliminated the problem of mis-identified head waves while allowing for some flexibility in associating phase arrivals into events that were more than 200 km away from their source. Outside of this choice we adhere to almost all recommended values from within the PhaseLink docu- mentation for a volume corresponding to our study area. We give a complete list of the parameters used to both train and apply the associator model in supplementary table S18, with justifications for any modifications from recommended values. We chose our best fit model as the one that achieved the lowest loss values during training. After selection of our best fit associator model we then applied it towards the pick database that we generated using EQTransformer in order to assemble a catalog of potential earthquakes. We chose to define an event as containing at least 6 associated P arrivals, with 52 no constraints on the number of S arrivals present, as EQTransformer picked a relatively small number of S arrivals. After the initial run of the association scheme, we then back- projected through time for each associated event in order to incorporate potentially missed phases. We did this because our model had a higher precision than recall, meaning excel- lent rejection of false associations but more prevalent discarding of true associations. The final run of the association scheme on our pick database generated 5142 candidate earth- quake events. From there we solved for initial event locations using the hypoinverse software package (Klein, 2002) using a 1-D regional seismic velocity model (Ma and Audet, 2017). We then screened out poorly constrained events by filtering out all earthquakes with root mean square travel-time residuals greater than 1 second. Local magnitudes for events were computed in a format consistent with the Canadian National Seismic Network using the formulation of (Nuttli, 1973). P displacements for the magnitude calculations were computed for each event within a 2 second window of detected Pg phase arrivals using bandpass filtered waveforms between 1.0-12 Hz with cosine tapers between 1.0-1.5 Hz and 10-12 Hz on respective ends of the frequency band. Local magnitude uncertainties were estimated from the variance of individual station magnitudes and are on the order of 0.2 - 0.4 for most events. Since the local magnitude scale that we use is not calibrated for our specific region we do not focus our analysis on moment release rates. We performed event relocations using the TomoDD-SE package (Zhang and Thurber, 2003) using the absolute Vs model of (Est`eve et al., 2021). As no published models for regional-scale absolute P velocity currently exist, we needed a method to extrapolate this Vs model into a workable absolute Vp model using published Vp/Vs ratios. The approach that we used was to spatially average the Vp/Vs ratios published in (Fern´andez-Viejo et al., 2005) for the SNoRCLE seismic line 31, which ran from ∼ 100 km east of Whitehorse in the central Yukon Territory to ∼ 100 km southwest of Macmillan Pass at the Yukon/ Northwest Territories border. This profile was chosen as it evenly samples the Cordilleran lithosphere on either side of the Tintina fault and is taken to be the most representative published model of 53 our study region. While this assumption undoubtedly leads to an oversimplification of Vp/Vs structure, we note that lateral variations in Vp/Vs are very small within the Cordilleran crust and only show significant variation (∼ 5%) in the Mackenzie foreland belt close to the Cordilleran Deformation Front (Fern´andez-Viejo et al., 2005). Since our study is focused on local earthquake distributions and it is outside of our scope to solve for velocity structure we hold the velocity model fixed in the tomoDD relocation process. Error estimates for initial event locations are computed using hypoinverse (Klein, 2002), which approximates error via the eigenvalues of the spatial elements of the covariance matrix for each location solution. These principal values are used to define an error ellipsoid, of which the largest such value is always the depth error. Since the covariance matrix is computed from partial derivatives around the event hypocenter it is not based on discontinuities within the earth or path-related uncertainties so it is always an underestimate (Klein, 2002). The tomoDD event relocation scheme operates on events by minimizing iteratively re-weighted location errors based on user-defined parameters (Zhang and Thurber, 2003; Waldhauser, 2001). A strength of this method is that velocity model features are accounted for, however it begins its iteration with the error estimates from hypoinverse and is thus biased by the problems mentioned above. We denote catalog error statistics in Supplemental Figures S18 and S19 and show location uncertainties for key regions in Supplemental Figures S13 and S14. We report initial event location errors for events that were not relocated and relocated errors for events that were relocated. We further note that the highest 1-sigma uncertainty values in the velocity model are ∼ 0.3 km/s, so propagating this error outwards over most event-receiver distances yields a maximum 1-sigma location uncertainty of ∼ 5 km due to velocity model uncertainties alone. Adding this to the reported catalog uncertainties gives error estimates of < 7 km in both horizontal and vertical components for the vast majority of events. 54 Figure 4.4 Top: Regional seismic catalog Aug 2016 - Sep 2018. Cross section lines for A-A’ and B-B’ profiles shown. Dashed box shows inset region plotted below in the bottom left. Solid black lines are major regional faults from Figure 4.1. Stations are color coded by network and earthquake sizes are scaled by magnitude. TC = Tintina Cluster, SC = Selwyn Cluster. Bottom Left: Close-up view of events in the The Eastern Denali Fault system. Blue box gives location of bottom right inset Bottom Right: Relocated hypocenters in the May 1, 2017 doublet aftershock sequence. Focal mechanisms show GCMT solutions for the two main shocks. DRF = Duke River Fault, KRF = Kluhini River Fault. 55 4.4 Results The catalog that we present here contains the most spatially dense sampling of low- magnitude seismicity across the NCC to date, and with smaller location errors than the CNSN catalog locations. In total we detect 4107 events over the 25 month interval with 6 or more associated phase arrivals that pass our quality control criteria (Figure 4.4). Most of the earthquakes occur seaward of the Denali Fault in the southwest portion of the study region. This particular area in the eastern St. Elias orogen shown in the bottom panels of Figure 4.4 contains 3279 events, around 75% of the catalog, and is dominated by the May 1, 2017 doublet earthquakes and their associated aftershock sequence (discussed further in section 4.3). The CNSN catalog covers our entire study area, while the Alaska Earthquake Center (AEC) catalog also overlaps with the western portion. Our catalog includes 1849 events that are not included in any other catalog. Of the 4107 events that we present in this study, 2035 correspond with events in the AEC catalog and 995 correspond with events in the CNSN catalog. All of the corresponding events between our catalog and the AEC catalog are located in the St. Elias orogen, with only one event located east of the TA footprint. Of the 995 events that correspond with events in the CNSN catalog, 274 of these were adopted into our catalog, whereby we used our location/ relocation procedure to arrive at better depth constraints. We find a magnitude of completeness (Mc) of 1.8 for our catalog, substantially lower than the Mc 2.8 from the CNSN catalog over the same time period (Figures S2 & S3). Northeast of the Denali fault, earthquake activity becomes spatially diffuse and maximum event magnitudes rarely exceed M 4.0. The most obvious zone of seismicity inboard of the Denali fault is located in the Richardson and Wernecke Mountains of northeast Yukon Territory (Figure 4.4). That cluster, previously analyzed by (Leonard et al., 2008), will be discussed further in section 5.1.2. Two distinct clusters of low magnitude (M < 3.0) earthquakes are located northeast of the Denali fault, which we refer to as the Tintina and Selwyn clusters (TC and SC Figure 4.4). The Tintina cluster lies within 50 km on either 56 side of the surface trace of the Tintina fault and strikes southeast across south-central Yukon Territory for ∼ 400 km from the transfer zone with the Teslin fault to Watson Lake. The Selwyn cluster is located within the Selwyn Basin ∼ 170 km northeast of the Tintina fault. The Selwyn cluster is slightly more diffuse than the Tintina cluster and occupies an area that is approximately 150 km x 150 km in map view. Seismicity is broadly present within the Selwyn Basin, however the area occupied by the Selwyn cluster is confined to its northeast extent and is bounded approximately to the northeast by the Niddery High. 4.4.1 Earthquake Depth Distributions As a result of the depth relocation process and use of more accurate velocity models, we are the first study to be able to report reasonably well constrained depth solutions for earthquakes in the NCC. Figure 4.5 shows depth sections from lines indicated in the top panel of Figure 4.4. Profile A-A’ gives depth distributions for events along the primary transect spanning from the northeast Gulf of Alaska to Great Bear Lake. Intense seismicity including numerous M > 4.0 events are detected throughout the crust southwest of the Denali fault. The Denali fault itself is characterized by a relatively narrow (∼ 40 km wide) band of seismicity which is dominantly made up of the May 1, 2017 doublet aftershock sequence (discussed below) and extends at depth almost to the Moho. Northwest of the Denali fault, a ∼ 50 km-wide band of seismicity is observed to decay in concentration away from the surface trace of the fault. Most of these events just inboard of the Denali fault are located in the uppermost 15 km (Figure S1). The Tintina fault hosts a substantial amount of low-magnitude seismicity throughout the uppermost ∼ 25 km of crust. The width of the densest portion of this band of seismicity is around 30 km, similar to that which is associated with the Denali Fault (Figure 4.5). More diffuse earthquakes associated with the Tintina cluster extend for an additional ∼ 30 km northeast of the core of the Tintina fault. Northeast of the Tintina cluster is a relatively narrow (∼ 100 km wide) region with few earthquakes. The Selwyn cluster is located just northeast of this quiet zone. The earthquakes in the Selwyn cluster extend to ∼25 km depth 57 Figure 4.5 Cross section A-A’ and B-B’ profiles from top panel of Figure 4.4. Solid black lines gives depth to Moho from (Audet et al., 2020). Magnitude scale shown on the right. Intersections of major structural features labeled on top of each profile line. Top: A-A’ profile from the Gulf of Alaska to Great Bear Lake. Inferred geometry of the Plateau Thrust after (Cecile and Cook, 1981) shown as solid black line. CDF = Cordilleran Deformation Front. Bottom: B-B’ profile running from southeast Yukon Territory to the Richardson Mountains. Inferred position of the Mackenzie Craton from (Schaeffer and Lebedev, 2014) indicated. 58 and are more diffuse in their depth distribution when compared to the Tintina cluster. By around km 800 of the A-A’ profile line (Figure 4.5) the trend of Selwyn seismicity and indeed seismicity in general terminates some 300 km southwest of the Cordilleran Deformation Front. The MMEP line of seismometers extends past the Cordilleran Deformation Front, so this lack of earthquakes is not related to detection thresholds, and represents a true lack of earthquakes. Profile B-B’ (Figure 4.5) gives depth distributions for events along a profile extending from the Liard zone of northern British Columbia to the Richardson Mountains. Earthquakes in the southeast Mackenzie Mountains are sparse and lie within the uppermost 10 km of crust. The seismicity of the Selwyn cluster is contained between km 250 and km 400 and extends to a depth of 25 km. This cluster includes many events located near the profile, with excellent station coverage for determination of depth. A gap in seismicity at depth is present northwest along the profile transect line between km 400 and km 550. To the north of this gap is the relatively intense activity of the Wernecke Mountains, which predominantly lies within the upper 20 km of crust. Further to the north in the Richardson Mountains we detect small-magnitude, relatively deep and diffuse seismicity between 10-30 km depth. Figure 4.6 compares earthquake depth distributions across various spatial bins. Peak seismicity for the Denali fault lies between 9-17 km depth, whereas peak seismicity for the Tintina fault is between 9-15 km depth. Most seismicity around the Tintina fault is shallower than 15 km, although there is a tail to the depth distribution that extends to 21-24 km depth. The Selwyn Basin exhibits an almost linear decay in seismicity with depth beneath the seismogenic zone, which we interpret to be indicative of motion along inherited structures. We detect a bimodal earthquake depth distribution with a significant depth peak centered at 27-30 km in the Richardson Mountains, which may be indicative of regional strain transfer at depth. 59 Figure 4.6 Depth distributions of only relocated seismicity across spatial bins. Note the progressively decreasing x-axis limits between the top and lower panels. 60 4.4.2 May 1 2017 M 6.2 and 6.3 Earthquake Doublet On May 1, 2017 a pair of M 6.2/6.3 earthquakes ruptured the eastern portion of the Denali fault system at the junction of the Denali and Duke River Faults. The first was a M 6.2 thrust earthquake located between the Denali and Duke River Faults and the second was a M 6.3 strike-slip event located along the Duke River Fault (Figure 4.4). The rupture source properties were estimated by (Feng et al., 2019) and (He et al., 2018), who both found that the first earthquake was a thrust event along a steeply dipping (∼ 60◦ − 65◦) fault and that the second was a strike slip event along a near-vertical fault. Rupture directivity analysis performed by (He et al., 2018) found that the first thrust event rupture progressed updip along a southwest-dipping fault plane, and that the second strike-slip event ruptured east-southeast along a left-lateral transform fault. The relocated aftershocks that we detect delineate two distinct fault planes that trend sub-parallel to the GCMT nodal planes (Figure 4.4). The aftershocks from the first thrust event are almost exclusively located between the Denali and Duke River faults and trend sub-parallel to the strikes of both. The aftershocks from the second strike-slip event lie mostly east-southeast of the main shock, in agreement with the directivity analysis of (He et al., 2018), however several hypocenters also extend to the west and almost all lie within 0-5 km south of the main shock. We note that this trend of aftershocks strikes subparallel to the most southeastern extension of the Duke River fault. Supplemental Figure S1 shows orthogonal depth sections both across and along the Denali fault centered on the doublet region. The pattern of low-magnitude aftershocks that resulted from the doublet events continued for almost a year and extended to ∼ 25 km in depth, suggesting that almost the entire crustal width of the Denali fault was ruptured. It is difficult to discern a time migration pattern for the primary swarm of aftershocks, however a small number of events occurring around a year later located at ∼ 17 km depth extend to near the Denali/ Totschunda fault junction. Interestingly, the Denali/ Totschunda fault junction hosted 3 M<3 foreshocks prior to the doublet event as well as a M 4 aftershock that occurred around 2 months after the doublet main shocks. This suggests a high degree 61 of connectivity between the states of stress of the Denali and Totschunda faults. 4.5 Discussion The first order pattern of newly detected M < 3.0 events is generally consistent with that of previous regional catalogs, although ours has improved completeness for the 2016- 2018 time period and smaller location uncertainties. To the southeast it is clear that the Eastern Denali fault system and the Fairweather fault play primary roles in accommodating transpression caused by ongoing Yakutat collision. This transpression manifests itself both in the observable geodetic displacement field (i.e. Elliott et al., 2010; Leonard et al., 2007) as well in complex earthquake sequences such as the aforementioned May 1, 2017 doublet. A clear and rather abrupt cessation of seismicity is located in the southeast of the study region which roughly corresponds to the Liard Transfer Zone (Figure 4.1). The Liard Transfer Zone is an inherited structure dating back to Proterozoic rifting of the Laurentian margin, and has been interpreted as having separated different styles of extension (i.e. Hansen et al., 1993; Lund, 2008). A strong contrast in upper mantle anisotropy has been observed across the LTZ (Audet et al., 2016) which is likely inherited from this Proterozoic rifting episode and fundamentally limits the amount of tectonic strain that is able to permeate south. The LTZ is also roughly co-located with the extent to which tectonic stress from Yakutat collision in the far field is modeled to be essentially nonexistent (Koons et al., 2010; Marechal et al., 2015). Perhaps the most surprising absence of seismicity is in the Mackenzie Mountains northeast of the Selwyn Basin. This area has hosted moderate magnitude (M>4) earthquakes in the past(Figure 4.1), which are noted to be diffuse and sporadic in nature. We locate no events in the region from 2016-2018. We attribute the lack of observed seismicity in the northeast Mackenzie Mountains to the short temporal sampling of our catalog, but this also suggests that M>4 events occur relatively rarely. Our reported event depths are generally in good agreement with estimates of Te from topography-gravity coherence (Fl¨uck et al., 2003). The relatively narrow band of peak 62 Figure 4.7 Seismicity in the Richardson Mountains and northwest Mackenzie Mountains. Major named faults plotted in black, all mapped faults plotted in grey. KF = Knorr Fault, DeF = Deception Fault, DF = Deslauriers Fault, TF = Trevor Fault, TT = Tombstone Thrust, RST = Robert Service Thrust, SRF = Snake River Fault. Events are color coded by depth, color scale shown on bottom. Cyan Squares give station locations. Focal mechanisms give body wave (20-50s) solutions for M>4.4 events (see supplemental figures S6-S8). 63 seismicity with depth between 9-15 km along the Tintina fault matches closely with the reported Te values between 10-15 km. Further northeast in the Selwyn basin the increase in earthquake depths also corresponds to a thickening of Te to 20-25 km. In the Richardson Mountains the reported Te values are remarkably high, on the order of 40-50 km, which is much deeper than the deepest events that we detect. We note that topography-gravity coherence methods of determining Te can be problematic in regions of ongoing orogenic activity. The high seismic activity of the Richardson Mountains and its connection to active tectonics of the Beaufort margin (Esteve et al., 2022) suggest that this is indeed a region of active orogeny. We further note that the Moho in this region determined from receiver function analysis lies at ∼ 37 km depth (Audet et al., 2020), which corresponds with the maximum relocated depths for regional events that we report (Figure 4.6.) 4.5.1 Correlation With Mapped Structures Regional stress throughout the NCC is generally understood to be oriented dominantly towards the Northeast (i.e. Mazzotti and Hyndman, 2002; Leonard et al., 2008) as a result of Yakutat transpression. Recent work has highlighted the importance of the additional effect of southward-oriented mantle traction forces (i.e. Finzel et al., 2015) that extend through the Richardson mountains and are deflected by subducting Yakutat lithosphere in central Alaska. This deflection is modeled to re-orient towards the northeast in the NCC (McConeghy et al., 2022). In this section, we correlate the observed seismicity with previously mapped structures to better describe how these regional stresses relate to tectonically active faults. 4.5.1.1 Tintina Fault Lithoprobe reflection lines suggest that the damage zone that cores the Tintina fault is ∼ 30 km wide (Snyder et al., 2005), which is in agreement with the width of the core of the Tintina cluster (Figure 4.5). We find that most of the seismicity that locates closely with the surface trace of Tintina fault occurs near Ross River in the central Yukon Territory, where our detection capability is greatest. This narrow band of earthquakes extends for ∼ 50 km along strike and is co-located with the intersection of the Tintina fault with the MMEP 64 Figure 4.8 Seismicity overlain on terrane boundaries from Colpron & Nelson (2007) inboard of the Tintina Fault. NAb = North American Basinal strata, filled in grey. NAm = North American miogeoclinal strata, IM = Intermontane terrane sliver, MRE = Meilleur River Embayment. seismometer line (Figure 4.4). For greater distances from the MMEP line, we suspect that very small earthquakes could be occurring, but not detectable. Southeast of Ross River we note a parallel band of earthquakes that are located ∼ 20 km southwest and outboard of the Tintina Fault, suggesting activity along a parallel structure such as the St. Cyr Fault. Northwest of Ross River we observe a broad trend of low magnitude earthquakes that extend on either side of the Tintina Fault. The moderate- magnitude earthquakes in the transfer zone between the Teslin and Tintina faults suggests that transfer faults may play a role in accommodating regional stresses. The Tintina fault 65 is oriented almost orthogonally to the regional stress field, so the fact that it displays such an abundance of low magnitude earthquakes suggests that it represents a zone of relative weakness within the Cordillera. The Tintina cluster contains the majority of the seismicity located between the Denali Fault and the Selwyn Basin, suggesting that the Tintina fault is presently active at a low rate. (Leonard et al., 2008) estimated a deformation rate of 0.5 mm/yr for the Tintina fault in the Yukon Territory based on instrumantally recorded seismicity, which is generally consistent with our findings. River drainage development modeling in the vicinity of the Tintina fault suggests that lateral displacements along the Tintina fault postdate the Paleogene and may continue to present (Ryan et al., 2017). Late Pleistocene and Holocene offsets are also recorded along the Preacher Fault and Medicine Lake Lineament sections of the Tintina fault in eastern Alaska (Plafker et al., 1994). While these offsets are attributable to the regime of bookshelf faulting in interior Alaska inboard of the Denali fault, they provide additional evidence of a sub-optimally oriented regional stress field re-activating sections of the Tintina fault. 4.5.1.2 Northwest Yukon Territory We detect almost no seismicity in the vicinity of the Dawson, Tombstone, and Robert Service thrust faults, which appear to be relatively inactive at present (Figure 4.7). In the Richardson Mountains, we note that the northern extent of seismicity correlates with the surface trace of the Knorr Fault, and that the general core of the Richardson Anticlinorium appears to be the regional locus of deformation. Further southeast, a M 5.1 earthquake and aftershock sequence appears to be closely located with the intersection of the Trevor fault and Snake River fault in the Eastern Wernecke Mountains. The trend of moderate- magnitude earthquakes continues further south, generally occurring southwest of the Snake River fault and focal mechanisms suggest a combination of thrusting and dextral strike-slip motion. The structures in the Richardson Mountains that we identify as displaying high levels of seismicity are all aligned essentially north-south, subparallel to the orientation of 66 the inferred mantle traction field in the north of our study region (i.e. McConeghy et al., 2022). 4.5.1.3 Selwyn and Mackenzie Mountains The Selwyn cluster appears to be generally confined in map view to the northeast Selwyn Basin in the foreland of the Mackenzie Mountains (Figure 4.4). The Plateau thrust is a major thrust fault that outcrops further northeast in the Mackenzie Mountains, and was suggested by (Leonard et al., 2008) to be an active structure based on seismicity patterns. A duplex structure for the Plateau thrust beneath MacMillan Pass was suggested by (Cecile and Cook, 1981) based on the presence of an antiform northeast of MacMillan Pass that is interpreted as a ’pop-up’ structure formed over a duplex/ ramp geometry at depth. The events that comprise the Selwyn cluster are mostly located along or in the hanging wall of the Plateau thrust (Figure 4.5), which is estimated to lie at 10 km depth beneath MacMillan Pass by (Cecile and Cook, 1981) and 12 km by (Cook et al., 2004). Only a small number of earthquakes are located in the footwall of the Plateau thrust (Figure 4.5). Splay faults are commonly observed to extend into the hanging wall of regional-scale thrust faults and we interpret seismicity in the hanging wall of the Plateau Thrust as partially occurring along blind thrust faults at depth. The vast majority of seismicity inboard of the Tintina fault occurs within Paleozoic strata of the Selwyn Basin (NAb on Figure 4.8). The Selwyn Basin is composed of early to mid- Paleozoic interbedded shales, cherts, and limestones that were heavily deformed and thrust faulted during Mesozoic convergence. The juxtaposition of these rocks against the slope- to-shelf sequences of the Mackenzie Platform evidently represents a significant contrast in crustal strength. The Misty Creek Embayment of the northeast Mackenzie Mountains (Fig- ure 4.2) is entirely seismically quiet, implying that lithological controls alone are not adequate to explain the observed earthquake clustering behavior. In the region near MacMillan Pass, the Niddery High appears to roughly correlate with the northeast extent of Selwyn Basin seismicity. The Niddery High is an ancient submarine high that separated Paleozoic rift 67 basins outboard of Laurentia (Cecile and Norford, 2000). At present it forms the boundary between the Selwyn Basin and the North American miogeoclinal sequences (NAm on Figure 4.8). We detect several (M> 3.0) earthquakes just northeast of the axis of the Niddery High, which are higher in magnitude than all of the other seismicity located within the Selwyn Basin. Further to the southeast the Meilleur River Embayment does not appear to host sig- nificant levels of seismicity. This is in spite of a rich history of earthquakes in the southeast Mackenzie Mountains including the 1985 Nahanni Earthquakes. We attribute this discrep- ancy to the short temporal sampling of our catalog, but this also suggests that sequences such as the Nahanni Earthquakes are rare. 4.5.2 Implications for Asthenospheric Processes The spatial locations of the two seismicity swarms that we identify roughly co-locate with the center line of the Mackenzie Mountain fold-thrust belt. The more diffuse Selwyn cluster is in the near foreland of the belt while the Tintina cluster is located some 300 km outboard of it. Seismic tomography indicates the presence of two coherent blocks of thick cratonic lithosphere that underlie the two ends of the Mackenzie Mountain fold-thrust belt (i.e. Schaeffer and Lebedev, 2014; Est`eve et al., 2021). These are interpreted to act as a channel for mantle flow in the asthenosphere, which would be directed away from the Gulf of Alaska and towards the Mackenzie Mountain fold-thrust belt (i.e. McConeghy et al., 2022; Finzel et al., 2015). Inherited structures in the Selwyn and Mackenzie Mountains are oriented favorably for reactivation by such stresses, and the observed high heat flow and depth distribution of seismicity serve to further reinforce this mechanism of stress transfer into the Northern Cordillera. This picture is complicated by the more concentrated small earthquakes that we detect along the Tintina fault, which is oriented almost orthogonally to the modeled principal stress direction. We interpret this phenomenon to the fact that the Tintina fault represents a zone of relative weakness which may be reactivated by smaller components of non-orthogonally oriented regional stresses. A similar phenomenon is also observed along the Northern San 68 Andreas fault by (Castillo and Ellsworth, 1993), who noted a rotation in deformation along inherited structures to conform with the modern stress regime. It is also evident from depth distributions of events that the seismicity of the Selwyn cluster generally picks up some 100 km northeast of the Tintina Fault. This shadow of earthquakes may suggest that the Tintina fault indeed accommodates a small amount of regional strain despite its almost orthogonal orientation to the modeled regional stress field (Mazzotti and Hyndman, 2002). The Selwyn cluster spans the uppermost 25 km of crust and we interpret it to be a consequence of localized thrusting along inherited structures consistent with (Mazzotti and Hyndman, 2002). The overwhelming majority of these earthquakes are located some 200 - 300 km outboard of the Cordilleran Deformation Front (Figure 4.8), implying that the strength contrast between the Mackenzie Mountain fold-thrust belt and the cratonic margin is not a primary cause for localization of strain. Instead, it appears that the forces driving regional deformation extend only as far as the near foreland of Mackenzie Mountains and are largely confined to the Selwyn Basin. The instrumentally recorded earthquake record complicates this picture and instead suggests that significant moment release occurs up to 100 km inboard of the Selwyn Basin (Figure 4.1). We attribute this discrepancy to the short temporal sampling of our catalog, and would expect to see discrete patterns of seismicity occurring within the Mackenzie Mountains given a longer observational period. Despite this, we interpret the paucity of seismicity inboard of the Selwyn Basin to be indicative of a decay in the magnitude of strain transfer in the Mackenzie Mountains. We note that the Plateau Thrust extends at depth beneath the Selwyn Basin and likely acts as a vehicle for stress transfer into the Mackenzie Mountain fold-thrust belt. The nature of the linkage between seismicity in the Richardson and Mackenzie Moun- tains is shown in the bottom panel of figure 4.5. A zone of sparse seismicity separates the two regions at depth, which generally corresponds with the spatial extent of the proposed Mackenzie craton (Schaeffer and Lebedev, 2014). The Mackenzie craton is a block of cratonic lithosphere that was interpreted as having been chiseled off and displaced by the Tintina 69 fault to its present location (Est`eve et al., 2021). To the northwest of this quiet zone is the relatively intense seismic activity of the Wernecke/northwest Mackenzie Mountains, which includes a M 5.0 event and aftershock sequence. We interpret these high rates of seismic- ity to be the result of juxtaposing relatively weak basinal strata of the Richardson trough against the Ogilvie Platform (Figure 4.2). This activity grades to the northwest into deep (10-30 km), small-magnitude events in the Richardson Mountains (Figure 4.5). The linkage between these zones suggests that a southward-oriented regional force, ostensibly due to mantle tractions acting at the base of the lithosphere, may be a significant driver of tec- tonic stress in the Richardson and Wernecke/northwest Mackenzie Mountains. It is difficult to infer how far south this force extends, or if indeed it is even a factor southeast of the Mackenzie craton. 4.6 Conclusions The catalog that we present here is the most spatially dense compilation of regional earthquakes across the Northern Canadian Cordillera to date. About 1849 of the events were not in any previous catalog, 341 of which are located across the Mackenzie Mountains. The location accuracy of the catalog allows us to propose for the first time which specific faults are the most likely candidate structures that take up active tectonic motions in the region inboard (landward) of the Denali fault. Our improved depth resolution offers insight into the active seismogenic depths, and are in good agreement with independent estimates of effective elastic thickness (Te) through- out the Cordillera. We detect relatively high levels of low magnitude seismicity along a ∼ 400 km stretch of the Tintina Fault, implying that both it and nearby subsidiary struc- tures are presently active. We interpret activity along the Tintina fault as reactivation of a relatively weak structure that is sub-optimally oriented to the regional stress field. We also detect a seismicity swarm in the northeast Selwyn Basin, which is well delineated by a lithological transition from Paleozoic basinal to marginal strata in the Mackenzie Mountain foreland. This transition may be due to some combination of a lithological strength contrast 70 or localized basement ramp geometry along the northeast extent of the Selwyn basin. The connection of these basement structures to thrust faults further northeast in the Mackenzie Mountains may provide a critical link between seismicity identified in this study and prior earthquakes located further inboard. 71 CHAPTER 5 GEODETIC STRAIN TRANSFER, TECTONIC BLOCK MOTION, AND 3D GLACIAL ISOSTATIC ADJUSTMENT IN THE NORTHERN CANADIAN CORDILLERA 72 ABSTRACT Seismicity patterns and geodynamic models suggest that the Northern Canadian Cordillera is undergoing active tectonic deformation. We present the most spatially dense map of GPS observations across the region including 3 new continuous sites installed in the central Mackenzie Mountains in order to investigate the length scales over which geodetic strain permeates into the North American continental interior. We use these data along with data from GPS sites located further east in the Canadian Shield in order to evaluate models of 3D predictions for the effects of Glacial Isostatic Adjustment (GIA) resulting from both post- Last Glacial Maximum deglaciation on the North American continent as well as post-Little Ice Age deglaciation in Southeast Alaska. When GPS site velocities are corrected for the 3D effects of GIA a clear pattern of northeast-oriented convergence at a rate of 3-5 mm/yr is observed across most of the Northern Canadian Cordillera, extending at least as far inboard as the Mackenzie Mountain foreland. We use these new site velocities to constrain a re- gional tectonic block model and evaluate whether the Tintina fault and southern Richardson Mountains could plausibly define present-day block boundaries. Our findings suggest that inclusion of the Tintina fault as a block boundary represents a significant improvement in fit to GPS site velocities and that maximum allowable strike slip and tensile deformation rates on the Tintina fault are 1.5 mm/yr. GPS site velocities in the Richardson Mountains favor a more southward directed block motion relative to sites further south in the Cordillera, however this feature is not significant enough to warrant a separate block due to large gaps in station coverage. 5.1 Introduction 5.1.1 Tectonic Setting Tectonic activity in the Northern Canadian Cordillera (NCC) is commonly linked to the Yakutat collision zone in the northern Gulf of Alaska. The simultaneous subduction and accretion of the Yakutat block to the southeast Alaskan margin is a primary driver of present-day deformation throughout the region (Mazzotti and Hyndman, 2002; Elliott 73 and Freymueller, 2020b). The Yakutat collision zone also marks a transition from dextral transcurrent motion along the Queen Charlotte-Fairweather fault system in the northeast Gulf of Alaska to oceanic-continental subduction in the Aleutian trench. The unique setting of the NCC, which is inboard and northeast of this transition makes it a case study into the factors which can influence length scales over which convergent tectonic deformation is able to extend into a continental interior. Seismicity patterns in the NCC are characterized by distributed M > 5.0 events which ex- tend northeast from the Yakutat collision zone for up to 800 km inboard into the Mackenzie and Richardson Mountains (i.e. Leonard et al., 2008) (Figure 5.1). (Mazzotti and Hyndman, 2002) explained these observed earthquake patterns using an orogenic float model which posits that tectonic stress is transferred inboard into the NCC along a lower crustal detach- ment. In this model the presence of the cold, thick cratonic lithosphere at the Cordilleran Deformation Front causes the concentration of tectonic strain within the relatively weak Mackenzie Mountain fold-thrust belt. Patterns of microseismicity detected by (Drooff and Freymueller, 2023) suggest that this deformation may only be accommodated as far inboard as the Mackenzie Mountain foreland, however it is uncertain whether these events are more closely linked to the Yakutat collision zone or to north-south transpression from the Richard- son Mountains to the north. Global mantle convection models predict that present day subduction of the Pacific plate beneath the Aleutian arc results in southward-directed mantle flow beneath northern Alaska (Forte et al., 2010). This buoyancy-driven flow field is thought to give rise to mantle traction forces (Finzel et al., 2015) which operate at the base of the lithosphere. Geodynamic model- ing of the interaction between the Yakutat flat slab and this southward-directed mantle flow suggests the existence of an inflection point in mantle flow orientation just north of the cen- tral Denali fault in southeast Alaska (McConeghy et al., 2022). Northeast of this inflection point, mantle flow is oriented broadly eastward toward the Mackenzie Mountains and the traction force exerted at the base of the lithosphere produces horizontal surface deformation 74 at a rate of ∼ 3 mm/yr throughout most of the NCC (McConeghy et al., 2022). Prior studies of regional GPS site velocities have suggested that observable geodetic deformation originating from the Yakutat collision zone is accommodated at least as far inboard as Whitehorse in Yukon Territory (Leonard et al., 2007). The velocities first pre- sented in (Leonard et al., 2007) were corrected for the postseismic trend of the 2002 M 7.9 Denali Fault earthquake using the model of (Freed et al., 2006), which is a good fit to early postseismic transient deformation close to the rupture but a poor fit to sites over 200 km northeast of the rupture in the NCC. Corrections applied by (Leonard et al., 2007) for the 3D effects of Glacial Isostatic Adjustment (GIA) only included the post-Little Ice Age (LIA) component of ongoing deglaciation in Southeast Alaska from ∼ 1500 AD - present and did not consider the longer lived effects of post- Last Glacial Maximum (LGM) unloading of the North American continental ice sheets. The block model of Alaska and Northwest Canada of (Elliott and Freymueller, 2020b) was the first to present corrected site velocities in the NCC for the horizontal effects of post-LGM deglaciation, however the corrections that they applied were for the ICE-3G deglaciation model (Tushingham and Peltier, 1991) which has since updated in successive iterations. The block model of (Elliott and Freymueller, 2020b) placed most NCC sites on a North Cordillera block, which is moves east-northeast relative to stable North America and also included a separate northwest coast block which contains GPS sites in the Richardson Mountains of northern Yukon Territory. The northwest coast block rotates with more southerly orientation of motion relative to sites further south in the NCC. In this study we present the most spatially dense distribution site velocities from of all available campaign and continuous GPS sites in the NCC and the Northern Canadian Shield. This includes timeseries from 3 new continuous stations in the central Mackenzie Mountains as well as an additional 15 years of data at sites in the Yukon Territory and southeast Alaska surveyed by (Leonard et al., 2007) and (Larsen et al., 2005). The more temporally complete timeseries at campaign sites which we present are less affected by data quality variability 75 Figure 5.1 Seismicity and tectonic setting of the Northern Canadian Cordillera. All earth- quake hypocenters for events M>3.5 between 1985 and 2022 plotted in blue. Focal mech- anisms shown for events from 1976-2023. Blue focal mechanisms for events M>5.0 be- tween 1976-1995 are from the GCMT catalog. Black focal mechanisms for events M>4.0 from 1995-2023 are from the Pacific Geoscience Centre (PGC) CMT catalog (Kao et al., 2012). GOA=Gulf of Alaska, YB=Yakutat Block, CDF=Cordilleran Deformation Front, RM=Richardson Mountains, WM=Wernecke Mountains, DRF=Duke River Fault, LTZ=Liard Transfer Zone, GBL=Great Bear Lake. 76 and time-dependent postseismic deformation than those presented in previous work. 5.1.2 Glacial Isostatic Adjustment Glacial Isostatic Adjustment (GIA) is the solid earth response to the redistribution of ice sheet and ocean water mass on the earth’s surface. It manifests itself in both a rapid elastic response (i.e. Thomas et al., 2011) as well as a long-lived viscoelastic signal that can persist for thousands of years after the disappearance of continental-scale ice sheets (i.e. Peltier, 1985). The vertical ’rebound’ effect of unloading is the primary component of GIA signal and is thus most commonly studied, however horizontal effects are also long-lived and can serve as an important constraint on both Earth structure and the spatial distribution of past loads. 5.1.2.1 Loading Histories The complex late Pleistocene to present glacial history of the North American continent is defined by episodic advances and retreats of continental-scale ice sheets. These are commonly broken down in the literature into four components known as the Laurentide, Cordilleran, Innuitian, and Greenland ice domes. Most reconstructions of ice thickness maintain that these were all connected to each other at Last Glacial Maximum (LGM), also known as the Wisconsin glacial episode from approximately 25-21 ka. The Innuitian ice dome was centered on the Canadian arctic archipelago and the Greenland ice dome was centered over present day Greenland. The Laurentide ice dome was centered over Hudson Bay in Northern Canada and extended as far south as the Hudson straits in New York and as far west as Great Bear lake in the Canadian Northwest Territories. The Cordilleran ice sheet was centered over the Canadian Cordillera in present day British Columbia and Yukon Territory, extending from the Mackenzie Delta in the north to the Oregon-Washington border in the south. The Cordilleran ice sheet terminated to the northwest in central Yukon Territoy (Roy and Peltier, 2017), leaving much of interior Alaska and western Yukon Territory deglaciated at LGM. This is likely due to a rain shadowing effect from the Alaska Range and St. Elias orogen to the southeast. For simplicity we adopt the naming convention of (Lambeck et al., 2017) and 77 will refer to the conglomeration of North American continental ice sheets both during and after LGM as the Late Wisconsin (LW) ice complex. The principal datasets used to constrain the evolution of LW ice mass are geologic ob- servations of changes in relative sea level. These include both near-field observations at pro-glacial lakes as well as far-field observations of global sea level. Ice budgets were first constrained by fitting observations of far-field sea level using the sea level equations (i.e. Spada and Stocchi, 2007). Regional observations enter into the workflow in later steps in order to fine-tune the loading functions for studies on the regional scale. Present-day obser- vations of geodetic uplift rates are sometimes introduced in the latter stages of the inversion in order to further constrain the distribution of ice load on a regional scale. This is part of the workflow used to obtain the ICE-6G (Peltier et al., 2015) and ICE-7G models (Roy and Peltier, 2017). These models also have the distinction of providing 3D predictions of present-day GIA signal on a global scale. Three different loading models for Laurentide and Cordilleran ice sheet deglaciation are examined in this study: The aforementioned ICE-6G and ICE-7G models (Peltier et al., 2015; Roy and Peltier, 2018) and the ANU-ICE model (Lambeck et al., 2017). The underlying mathematical formulations of both models are similar, however there are differences in the treatment of how ice loads behave near continental shelves and shallow basins as well as how observational datasets are related to solutions for past load distribution (Johnston, 1993; Peltier, 1998; Purcell et al., 2016). The ICE-6G model was optimized over the North American continent and in combination with it’s partner VM5a mantle viscosity model was constrained in the final stages of the workflow using space-geodetic observations. ICE-7G further updated the ICE-6G loading model by shifting a component of the Laurentide ice dome at LGM away from the area west of Lake Winnipeg to instead be located over Hudson Bay. This was motivated by a need to fit far field observations of Relative Sea Level in the western Mediterranean basin (Roy and Peltier, 2018). It is important to note that the loading history of ICE-6G and ICE-7G converge over North America at ∼ 10 ka and that 78 reconstructions of ice history are virtually identical from 10 ka to present (Roy and Peltier, 2017). The ANU-ICE loading model is comprised of a series of regionally optimized loads which separately account for the evolution of all major global ice sheets as well as localized vis- coelastic structure in each case (Lambeck, 1993; Lambeck et al., 1998; Lambeck et al., 2017). The North American component of the ANU-ICE model presented in (Lambeck et al., 2017) was optimized by assuming a different mantle viscosity structure beneath North America compared to the other contemporaneous continental ice sheets. While it stands to reason that different continents will have different underlying mantle viscosity structure, this ap- proach introduces different tradeoffs between ice load history and viscosity structure than unified global models like ICE-6G/ICE-7G. Ongoing mass loss at glaciers in the St. Elias orogen and southeast Alaska is a continuous process spanning back in time to the Little Ice Age (LIA) (∼ 1500 AD). Modeled GIA signal resulting from this post-LIA deglaciation peaks in amplitude at > 3.5 cm/yr of uplift in parts of Glacier Bay (Larsen et al., 2005; Elliott and Freymueller, 2020b). (Hu and Freymueller, 2019) presented the first model of regional post-LIA GIA signal. Much of the input dataset came from a digitization of previous ice extent maps and satellite imagery compiled by (Berthier et al., 2010). With the loading history treated as a known parameter, (Hu and Freymueller, 2019) solved for solid earth structure beneath southeast Alaska in order to fit GPS observations of present-day uplift. 5.1.2.2 Earth Models Global models of Glacial Isostatic Adjustment are presently obtained by assuming a radi- ally symmetric 1-D viscoelastic structure of the solid Earth. This is due to the prohibitively large computational expense of solving for both the spatiotemporal distribution ice load- ing function as well as the underlying 3D earth viscosity structure. The models discussed above all employ an iterative workflow by which a viscoelastic structure for the solid Earth is obtained in conjunction with a loading distribution. This is often performed on a regional 79 basis in order to overcome this oversimplification of earth structure (i.e. Lambeck et al., 2017). Numerous lines of evidence from seismology and mineral physics suggest that the composition of the Earth’s lithosphere and mantle is heterogeneous. Studies have only re- cently begun to appear in the literature which investigate the effects of 3D variations in the Earth’s viscosity structure on GIA signal at both the regional (Wal et al., 2015; Marsman et al., 2021) and global scale (Austermann et al., 2021). These studies typically only report differences in observed uplift signals between models that account for 3D viscosity variations as opposed to assuming a radially symmetric viscosity structure. Lithospheric thickness and asthenospheric viscosity are the key solid earth parameters in determining GIA signal. Both are commonly inferred from seismic wavespeeds, which are highly sensitive to temperature and density but not necessarily the mechanical boundary layer represented by the Lithosphere-Asthenosphere Boundary (LAB). Thus there is often an assumption made where authors infer a particular isotherm as representing the LAB. This approach approximates the boundary layer reasonably well on a global scale but likely falls short in regional-scale applications or when dealing with exceptionally old and thick cratonic cores such as is present in the inboard regions of the NCC. No current model exists which gives 3D predictions of GIA signal over the North American continent using a realistic 3D solid Earth viscoelastic structure. The Northern Canadian Cordillera is situated in a backarc setting with hot, thin litho- sphere that overlies a low-viscosity mantle. This is juxtaposed to the east near the Cordilleran Deformation Front against the cold, thick lithosphere of the Canadian shield (Hyndman et al., 2005a; Hyndman, 2010). Estimates of lithospheric thickness based on seismic tomogra- phy vary dramatically between ∼ 50 km in the Cordillera (Audet et al., 2019) and ∼ 150 km in the Canadian Shield (Yuan and Romanowicz, 2010). The transition between thin, hot Cordilleran lithosphere and thick, cold cratonic lithosphere happens somewhere beneath the Mackenzie Mountains, where a westward-dipping low-velocity structure is observed from seismic tomography to extend from the mantle to upper crust (Audet et al., 2019; Schutt 80 et al., 2023). 5.2 Data We use the most spatially dense available set of GPS stations throughout Eastern Alaska and Northwest Canada in this study. The dataset is comprised of both campaign and continuous measurements spanning from 1991 to 2024. This includes an additional 10-15 years of repeat measurements at all of the campaign GPS sites in Yukon Territory first surveyed by (Leonard et al., 2007) as well as data from sites along the Alaska highway first presented in (Larsen et al., 2005) with occupations running through 2017 at campaign sites and continuous recordings to 2024. We present new secular velocities for 3 continuous sites (MMEP, CJOU, NWEL) that lie on a transect of the Mackenzie Mountain fold-thrust belt. These were installed as in 2017 as a part of the MMEP project, which augmented the regional footprint of the USArray Transportable Array seismic deployment in Northwest Canada. We also utilize all available data from continuous GPS sites located throughout the Northern Canadian Shield operated by the Canadian Geodetic Survey. 5.2.1 GPS Data Processing Continuous and campaign GPS data were processed at the Michigan State University Geodetic Laboratory using the GIPSY/OASIS goa-6.4 software package (i.e. Bertiger et al., 2010) in Precise Point Positioning (PPP) mode. Daily solutions for site coordinates were obtained using reprocessed JPL orbit and clock products and Vienna Mapping Functions VMF1GRID (B¨ohm et al., 2009) were applied to correct for atmospheric path delays. Cor- rection models are also applied for solid earth tides as well as ocean tidal loading. Daily solutions were then aligned into the ITRF 2020 reference frame (Altamimi et al., 2023) in order to produce daily coordinate and covariance estimates of site positions. We refer the reader to (Elliott et al., 2024) for more technical detail on the data processing procedure. 81 5.2.2 GPS Velocity Estimation We use all available campaign and continuous GPS data in order to solve for secular site velocities across the NCC. Timeseries are corrected for coseismic offsets resulting from the 2002 Mw 7.9 Denali Fault earthquake, the 2011 Mw 7.5 Craig earthquake, and the 2018 Mw 7.9 Gulf of Alaska earthquake. GPS velocities are independently solved for in each spatial component (east, north, vertical) by parameterizing each timeseries using an equation of the form: y(t) = a + bt + c sin(2πt) + d cos(2πt) + e sin(4πt) + f cos(4πt) + H(t − td)[g + h ln(1 + ( t − td τL ) + k(1 − e−(t−td)/τE )] (5.1) where y(t) gives the estimated position of a single spatial component at time t, expressed in years. The constants a, b, c, d, e, f, g, h, k are all estimated by least-squares inversion using a program called tsfit in MATLAB. Not all of the terms in Equation 5.1 are used to obtain velocity solutions at all sites. The periodic terms describe annual and semiannual variations and are only estimated for continuous stations with year-round recordings. The postseismic decay terms are given in the last item of equation 5.1 and are used to describe transient deformation resulting from the 2002 Denali fault earthquake. H(t − td) is the Heaviside function where td is the time of the earthquake and τL and τE terms are respectively the logarithmic and exponential decay time constants used to describe relaxation times for postseismic motion. These are only used at sites that contain postseismic motion in their timeseries. See Section 5.2.2.1 for information about which sites were fit using postseismic transient terms. All site velocities are estimated assuming a colored noise model consisting of white plus flicker noise for potential time dependent noise in the data (Mao et al., 1999). The b term in equation 5.1 gives the secular velocity, which should be indicative of time- independent tectonic deformation. Secular velocities are then corrected for the effects of elastic strain accumulation from the Yakutat collision along structures southwest of the study area from the predictions of (Elliott and Freymueller, 2020b) as well as for the effects 82 of GIA according to a suite of earth and load models as detailed in section 5.3.2. 5.2.2.1 Postseismic Motion from the 2002 M 7.9 Denali fault earthquake Postseismic deformation resulting from the November 2002 M 7.9 Denali Fault earth- quake is pervasive in the western portion of our study area. Sites close to the rupture of the earthquake display obvious postseismic motion, the magnitude of which decays further away from the rupture. Sites located within ∼ 150 km of the earthquake contain obvious postseismic transient deformation signals and we model those velocities by including both logarithmic and exponential decay terms in Equation 5.1. For sites greater than 150 km from the rupture the postseismic signal is less obvious. Postseismic displacements detected from InSAR in the time period directly following the earthquake (2003-2005) suggest that transient deformation extended at least 200 km away from the surface trace of the Denali fault (Biggs et al., 2009), with this 200 km distance not necessarily representing an upper bound on the lateral extent of deformation. (Johnson et al., 2009) report postseismic dis- placements in excess of 200 km from the surface trace of the Denali fault and coseismic offsets are observed well in excess of 200 km away from the rupture (Hreinsd´ottir et al., 2006). In the absence of a physics-based model of Denali postseismic signal which adequately fits the long timeseries data located (> 150 km) from the rupture we needed to a method to determine whether the inclusion of postseismic terms into the velocity fit (Equation 5.1) represents a significant improvement in the secular velocity solution. A simple windowing of stations based on distance from the rupture is insufficient due to both the complex nature of postseismic processes as well as the variable time history of observations at GPS stations, as much of our dataset is comprised of campaign-style sites established in the months and years following the November 2002 earthquake. We perform a statistical F test in order to test the significance of improvement in fit to timeseries data related to the inclusion of postseismic terms in the velocity solution and report the results in Table 5.1. We show a map of sites tested relative to the surface trace of the Denali Fault rupture in Figure 5.2. The F statistics reported in Table 5.1 are computed between velocity solutions both including and excluding 83 Figure 5.2 Site locations in eastern Alaska/ Yukon Territory tested for postseismic signal. Circles denote continuous-style GPS sites whereas squares denote campaign-style GPS sites. Epicenter of 2002 M 7.9 Denali Fault earthquake shown as yellow star. Orange line overlain on Denali Fault shows lateral extent of the rupture. Sites plotted in white are within 100 km of the rupture and display obvious postseismic signal. Sites plotted in grey were tested for the presence of postseismic signal as reported in Table 5.1. Sites plotted in black are too distant from the rupture to display postseismic signal. site ti 2002.88 DAWS TOWH 2004.71 2000.68 8130 2001.78 I177 2006.74 AB41 2002.87 AC61 Y565 2002.34 DEST 1999.38 NSLM 2000.40 MOTD 2002.34 FE 22.37 4.54 32.74 0.05 8.61 4233.66 2.21 4.38 1.32 0.02 pE 3e-3 0.08 7e-4 0.83 2e-4 1e-5 0.04 0.047 0.30 0.98 FN 0.24 1.78 0.14 0.86 42.52 478.36 2.07 8.07 1.73 4.76 pN 0.64 0.23 0.72 0.66 1e-5 1e-5 0.14 0.01 0.21 0.03 significance? yes no yes no yes yes yes yes no yes Table 5.1 F and p statistics computed between models that include or exclude postseismic terms in their respective velocity solutions. Fe denotes the F statistic computed for data in the east-west spatial component, Fn is likewise for data in the north-south component. p value significance is attributed to any p < 0.05. ti gives first observing day at site expressed in decimal year. Column on right indicates whether postseismic terms were adopted at the particular site. 84 postseismic terms. The F statistic is computed for each individual component of each site according to: F = (WRSS1 − WRSS2) /(p2 − p1) WRSS2/(n − p2) (5.2) Where WRSS signifies the Weighted Residual Sum of Squares misfit for an individual spatial component of data. WRSS1 is computed for the less complex model that does not include postseismic terms and WRSS2 is computed for the more complex model that does include postseismic terms. p1 and p2 signify the number of model parameters (the number of terms in Equation 5.1) for velocity solutions that respectively do not and do include post- seismic terms. n is the number of independent observations contained within the timeseries. Because continuous GPS timeseries have a degree of temporal correlation between daily solu- tions we determine that the time interval at which the correlation degrades to insignificance is 5 days from an examination of covariance of the colored noise model. For campaign-style sites the number of independent observations is defined as the number of site occupations because the length of each occupation is typically less than 5 days. Positive F values signify an improvement in fit to the data by including postseismic terms in the velocity solution. p values are reported in order to assess the significance of improvement in fit to the data. We interpret significant improvement in fit to the data for any site that exhibits a p statistic < 0.05 in either of the horizontal components. This occurs at all sites depicted in Figure 5.2 except for TOWH, I177, and NSLM. Since all of these sites are located closer to the rupture of the earthquake than some of thse that pass the F test we choose to adopt postseismic terms at all sites. The postseismic trend at continuous sites and those close to the rupture are fit by introducing log and exp terms to Equation 5.1. These sites are fit well by the inclusion of both of these terms. For campaign sites further away from the rupture we found a strong tradeoff to exist between the log term and the secular velocity estimation. For this reason we fit the campaign sites I177, TOWH, DAWS, and 8130 using just the exp term for the postseismic correcition. We use uniform 85 time constant values of 0.12 yrs for τL and 15 yrs for τE at all sites that were determined to contain postseismic motion. This approach allows us to fit and subsequently remove the postseismic signal reasonably well at all sites while maintaining low uncertainties in secular velocity solutions. 5.3 Methods 5.3.1 Testing Reference Frames for Stable North America The theory of plate tectonics asserts that motion at the surface of the earth can be described by a series of plate rotations and that deformation is concentrated within relatively narrow plate boundary zones. The spheroidal shape of the Earth requires that motion of each stable plate interior be described by a rotation about an Euler pole (Morgan, 1968). Space geodetic techniques including global observations of GPS site velocities have been used to derive a series of global plate motion models that give the motion of each major tectonic plate as a rigid block rotation about a pole. These include GEODVEL (Argus, 2007), MORVEL (DeMets et al., 2010), and ITRF-PM (Altamimi et al., 2017; Altamimi et al., 2023). To first order there is general agreement between these models as to the general locations of each plate’s respective Euler pole and magnitude of rotation, however for some plates the choices made surrounding input datasets in the inversion have a substantial impact on the resulting solutions. None of these solutions were obtained after applying corrections from geophysical models to input datasets. Authors are instead selective in discarding certain data from sites or regions that are suspected of being influenced by processes such as GIA or tectonic deformation. In North America, global plate motion models do not consider the vast majority of available geodetic data from Alaska and Canada due to the presumed presence of GIA and tectonic motion in site velocities. As such, these models fit site velocities observed in the southeastern United States quite well, but can exhibit systematic residuals at higher latitudes (Ding et al., 2019; Kreemer et al., 2018). Another important feature that comes out of global plate motion models is an estimation of geocenter translation rate, also known as an Origin Rate Bias (ORB). Earth-observing 86 satellites all orbit the center of mass of Earth system so any reference frame that is defined from space geodetic techniques will have this point as the origin. Since the Earth is a dynamic system undergoing time-dependent redistribution of mass the observed geocenter has an apparent time-dependent translation effect. This manifests itself as a systematic pattern of residual site velocities that is observed across multiple tectonic plates once the predictions of a global plate motion model are removed. The ORB could represent error in the geodetic reference frame, or a systematic motion between the center of plate rotation and the center of mass of Earth system. We correct initial observations of GPS site velocities across Northern Canada for three different definitions of stable North America including respective estimates of the ORB correction in order to investigate systematic residual site velocities. The frames that we test are GEODVEL (Argus et al., 2010), ITRF2020-PM (Altamimi et al., 2023), and (Ding et al., 2019). After removal of North American plate motion, residual site velocities should be indicative of GIA and unaccounted-for ORB. We find systematic southward-directed residuals at all sites across northern Canada after correcting for all of the definitions of stable North America, with the GEODVEL model correction resulting in the lowest magnitudes. We presume that these are not due to GIA as they are orthogonally orientated to the direction of advance/retreat of the Laurentide ice dome. We instead postulate that this could be due to a bias in the plate motion estimate caused by intraplate strain within the North American continent, consistent with (Ding et al., 2019; Kreemer et al., 2018). Because the primary goal of this study is to investigate tectonic deformation in the NCC we adopt the GEODVEL definition of stable North America as this reference frame correction results in the lowest magnitude of southward-directed residuals throughout the Canadian Shield. 5.3.2 Evaluating GIA Models Forward models of GIA on a global scale for a variety of different load and solid earth viscosity models are given by (Steffen, 2021). The load models which were tested are ICE- 6G/ICE-7G and ANU-ICE as described in section 5.1.2.1. The earth models test a set of 87 lithospheric thicknesses that range from 60 - 150 km in increments of 30 km as well as two values for lower mantle viscosity (2 and 20 ∗1021 Pa*s). Forward predictions are also provided for the ICE-6G model and ICE-7G models using their companion solid earth structures: VM5a and VM7, respectively. All earth models tested are radially symmetric 1D profiles with Maxwell viscoelasticity, rotational feedback, and time-dependent coastlines. The Earth’s core is assumed to be inviscid and incorporated as a lower boundary condition (Steffen, 2021). Rheological parameters including depth-dependent density, Young’s modulus, etc., are taken from PREM (Preliminary Reference Earth Model) (Dziewonski and Anderson, 1981). We estimate velocities at GPS stations spanning across the Canadian Shield from the NCC to Hudson Bay in order to evaluate the different GIA models for our study region. Residual velocities fit to GPS sites after removal of predictions for all loading models and endmember lithospheric thickness models of 60 km and 150 km are shown in (Figure 5.3). Long wavelength horizontal velocities resulting from GIA should be oriented west-southwest across Canada, away from Hudson Bay. We observe a clear systematic trend of east-oriented residuals pointed towards Hudson Bay after removal of GIA model predictions for all earth models that are not VM5a and VM7 (Figure 5.3 panels c-f). We interpret this to a short- coming in the 1D viscoelastic solid earth parameterization used by (Steffen, 2021) which leads to an overprediction of the horizontal velocity component. It is outside of the scope of this study to solve for an optimized viscoelastic solid earth structure beneath the Canadian shield. We choose to adopt the ICE-6G with VM5a model as our best approximation of the post- Laurentide (LW) component of GIA over northern Canada. This is due to the VM5a model delivering the lowest magnitude of southward-directed residuals in the Canadian Shield when compared to other solid earth models. It is important to note that the ICE-6G and ICE- 7G loading models converge by ∼ 10 ka and are identical thereafter so it is impossible to distinguish between them from present day geodetic observations alone (Roy and Peltier, 2017). Since we only have 3D predictions of GIA for the ICE-6G and ICE-7G loading 88 Figure 5.3 Horizontal residual velocities fit to GPS sites in northern Canada after removal of the predictions for GIA signal resulting from all loading models combined with endmember lithospheric thickness models. Load and earth models indicated on the bottom left of all subplots. Panel b shows our preferred combination of load and solid earth models for LW deglaciation. 89 models alongside their companion mantle viscosity models (VM5a and VM7 respectively) our only option for using the vm5a earth model is to use the ICE-6G loading model. This model combination provides the lowest eastward-directed residual velocity field in Northern Canada (Figure 5.3) while also fitting the observed uplift rates reasonably well (see section 5.3.3). 5.3.3 Constraining Lower Mantle Viscosity Using Vertical GPS Velocities We use vertical GPS site velocities in order to distinguish between predicted GIA signals for different models of solid earth structure in the NCC. To do this we sample predictions of post-Laurentide GIA signal from the suite of earth and load models detailed in section 5.3.2 along transects covered by our GPS sites (Figure 5.4). All transects begin in the northeast Gulf of Alaska and terminate in the North American continental interior, with the exception of profile A which terminates in the Canadian Arctic. Figure 5.4 shows the predictions of only the LW component of GIA for the suite of earth models with observed GPS site velocities corrected for post-LIA signal predictions of (Hu and Freymueller, 2019) overlain. The predicted vertical LW GIA signal along Profile A is very low, around 0 mm/yr of uplift for all earth models. This is because ice mass from Laurentide glaciation was largely absent from western Yukon/ interior Alaska due to rain shadowing from coastal mountain belts to the north and south. Profiles B-D in Figure 5.4 show a clear bifurcation in predicted vertical velocities which arises from the two different lower mantle viscosities tested. The differences in predicted uplift are relatively small along Profile B, which similarly to profile A was largely unglaciated at LGM. The difference between estimates of lower mantle viscosity is most pronounced along Profiles C and D, in which the vertical velocity data show a clear preference for the lower 2E21 Pa*s lower mantle viscosity value. The section of North America sampled along profiles C and D was largely covered by the Cordilleran and/or Laurentide ice sheets at LGM and the observed uplift signal is thus reflective of underlying solid earth structure. There are small variations in predicted uplift for different lithospheric thickness estimates (Figure 5.4), however these are too small and our spacing of GPS stations 90 Figure 5.4 Vertical predictions of post-Laurentide GIA signal for the ICE-6G loading model and a suite of earth viscosity models sampled across profiles in the study region. Panel (a) gives profile locations along which GIA models have been sampled. GPS station locations are shown as cyan squares. Panel (b) gives predictions for each solid earth model. Profiles are color coded such that red lines correspond to a lower mantle viscosity value of 2E22 Pa*s and blue lines correspond to a lower mantle viscosity of 2E21 Pa*s. Color weights increase with increasing lithospheric thickness values as shown in subplot legends. Black line gives predictions for the VM5a earth model. Grey dots give observed GPS station velocities corrected for post-LIA deglaciation in southeast Alaska (Hu and Freymueller, 2019). 91 is too sparse to meaningfully distinguish between them. 5.3.3.1 Accounting for Lithospheric Thickness Variations Lithospheric thickness is a key parameter in determining GIA signal. The NCC is un- derlain by thin ∼ 50 km lithosphere (Audet et al., 2019) which abuts the exceptionally thick (∼ 150 km) cratonic lithosphere of the Canadian Shield (Yuan and Romanowicz, 2010). The difference between predicted present-day horizontal GIA deformation for the ICE-6G loading model interacting with a 60 km thick lithosphere as opposed to a 150 km thick lithosphere is ∼ 2 mm/yr (Figure 5.5), which is similar in magnitude to inferred tectonic deformation in the NCC (Leonard et al., 2008). To test the effect of solid earth structure on GIA signal in the NCC we correct observed GPS site velocities for the predictions of the ICE6G loading model interacting with 3 separate models of solid Earth viscoelastic structure: (1) The vm5a solid earth model- which was determined to be the best fit for sites on the Canadian Shield. (2) 60 km thick lithosphere model- which is taken to be the most representative solid earth model for the NCC. (3) A third ’joined’ model which was constructed by joining together the forward predictions of both models in their representative spatial areas. The forward predictions for each of these models are shown in Figure 5.5. The vm5a earth model was first presented in (Peltier et al., 2015) and further constrained in the North American context in (Roy and Peltier, 2017). It consists of 5 spherically symmetric layers that span from the lower mantle to the crust. The uppermost elastic crust layer is 60 km thick and is underlain by a 40 km thick mantle lithosphere layer with a high viscosity of 1 ∗ 1022 Pa*s. The upper mantle layers from 100 - 420 and 420 - 670 km depth respectively have the same viscosity values of 5 ∗ 1020 Pa*s. The lower mantle from 670 - 1260 km depth has a viscosity of 1.57 ∗ 1021 Pa*s and from 1260 - 2885.5 km depth has a viscosity of 3.23 ∗ 1021 Pa*s. The 60 km lithosphere model consists of 3 layers: a 60 km thick elastic crust, an upper mantle from 60 - 660 km depth with a viscosity of 4 ∗ 1020 Pa*s, and a lower mantle from 660 - 2880 km depth with a viscosity of 2 ∗ 1021 Pa*s. The ’joined’ model 92 Figure 5.5 Predicted 3D velocities for GIA models with differing lithospheric thickness values. First row gives vertical predictions, second row gives east-west predictions, and bottom row gives north-south predictions. Note that vertical velocities shown in panels a-c are plotted using a separate color scale from the horizontal velocities shown in panels d-i. Panels a,d,g give 3D predictions of the ’joined’ model with GIA signal determined by lithospheric thickness. Panels b,e,h give 3D predictions for a uniform lithospheric thickness of 60 km. Panels c,f,i give 3D predictions for a uniform lithospheric thickness of 150km according to vm5a. All GIA predictions assume an ICE-6G loading function. Magenta lines bound the zone within which GIA signal was fit with a spline function in the ’joined’ model. 93 was created by masking out predictions for the vm5a model west of the Canadian Shield and masking out predictions for the 60 km lithospheric thickness model east of the NCC. These predictions were then ’joined’ together by fitting a bicubic spline function using GMT surface with a tension factor of 0.25. The Cordilleran Deformation Front (CDF) was used to mask each grid file (Figure 5.5). To define the region with thickened North American cratonic lithosphere we mask out all predictions within 3 of longitude of the CDF. To define the region with thinned Cordilleran lithosphere we mask out all predictions within 5 of longitude of the CDF. This leaves a region 8 wide between craton and Cordillera within which GIA signal is approximated by the spline function (Figure 5.5). This width was chosen because it is the same wavelength over which a discrepancy between the forward predictions of GIA models which account for either 1D or 3D viscosity variations in the solid earth was identified by (Austermann et al., 2021). It is important that the fit surface in this 8-wide region not be overly rough, as this could insert artificial discontinuities into a corrected velocity field and thus bias block model results. 5.3.4 Tectonic Block Model Setup Tectonic convergence across the NCC can be modeled using a block model (i.e. Meade and Hager, 2005) whereby observed GPS site velocities are fit by a series of block rotations. Blocks boundaries are constructed by referencing geologically mapped faults, trends in seis- micity, or compelling evidence of regional-scale trends in geodetic velocities. GPS sites are then assigned to blocks based on their location and observed site velocities are fit by angu- lar rotations in a least squares inversion. block Euler pole locations and rotation rates are then obtained from these angular velocities and fault slip rates along block boundaries are derived. We use the BLOCKS package (Meade and Loveless, 2009), which is written in MATLAB, in order to constrain a tectonic block model using GPS site velocities. GPS site velocities are corrected for (1) North American plate rotation according to the GEODVEL realization, (2) elastic effects of St. Elias collision as predicted by (Elliott and Freymueller, 2020b), and 94 Figure 5.6 Block model configurations tested in this study. Panel a shows shows the most basic configuration with no Tintina or Richardson block boundaries. Panel b shows the configuration with no Tintina fault block boundary but a separate Richardson block. Panel c shows the configuration with a Tintina fault block boundary but no Richardson block. Panel d shows the configuration which includes both Tintina and Richardson block boundaries. Blocks labeled in panel d have angular velocities and Euler pole locations that are solved for in our inversion. Unlabeled blocks in panel d are held to a priori values from (Elliott and Freymueller, 2020b). 95 (3) regional GIA signal as decribed in section 5.3.3.1. All block boundaries are defined to be vertical (90 deg dip) and are held to a locking depth of 10 km, consistent with estimates of locking depths for the Eastern Denali, Totschunda, and Duke River faults (Elliott et al., 2010). Away from the Denali fault we lack the site distribution to perform a meaningful inversion for fault locking depth and thus hold all fault locking depths to be fixed at 10 km. GPS sites are too sparsely distributed throughout the region in order to perform a meaningful inversion for internal strain accumulation within our block solutions so we a priori enforce full internal rigidity for all blocks. We construct our model by starting with the eastern Alaska/northwest Canada block boundaries from (Elliott and Freymueller, 2020b). These adopted blocks include the Fair- weather block, the Eastern Denali block, and the Fairbanks block which are a priori held to the Euler pole and rotation values reported in (Elliott and Freymueller, 2020b). Since all site velocities have been corrected for motion of the North American plate we a priori enforce zero block rotation at all sites inboard of the inferred deformation front. We then test the effects of including two additional block boundaries in the NCC along the Tintina fault and at the boundary of the Richardson and Mackenzie Mountains (Figure 5.6). The Richardson block was a feature first shown in (Elliott and Freymueller, 2020b), who called it the Northwest Coast block. We will hereafter call it the Richardson block as it is located mostly in the Richardson Mountains of northern Yukon Territory. Our more spatially dense dataset in the Canadian Cordillera is well posed to further interrogate whether such a feature is warranted in a regionally optimized model. We subdvide Yukon Territory into two distinct tectonic blocks in model geometries that consider the Tintina fault as a block boundary (Figure 5.6 panels c and d). These are the Selwyn and Intermontane blocks. The Selwyn block is located northeast of the Tintina fault and the regional geology is mostly comprised of Paleozoic marine platformal to basinal se- quences of the Selwyn Basin. We choose to limit the northeast extent of the Selwyn block to the Mackenzie Mountain foreland, consistent with observed patterns of microseismicity 96 reported in (Drooff and Freymueller, 2023). This boundary is also roughly co-located with a dramatic change in inferred lithospheric thickness from seismic tomography (Schutt et al., 2023). The regional geology of the Intermontane block is comprised of the Intermontane group of terranes, a complex amalgamation of mid-Paleozoic to early Mesozoic sedimentary and volcanic sequences which were formed outboard of the ancestral North American mar- gin and have since been either uplifted and/or accreted due to outboard tectonic activity (Colpron et al., 2022). Major dextral transcurrent offsets along the Tinina fault which sepa- rates the Selwyn and Intermontane blocks are mostly surmised to have occurred during late Mesozoic time, with the Tintina fault becoming largely inactive by the Eocene (∼ 60 Ma) (Gabrielse et al., 2006). Prior studies of earthquake statistics have inferred significant levels (∼ 4 mm/yr) of tectonic deformation in the Nahanni region of the SE Mackenzie Mountains (Leonard et al., 2008). This is mostly due to the 1985 Nahanni earthquake sequence which was constituted of two events of Mw 6.6 and Mw 6.8 that occurred on Oct 5 and Dec 23 1985 respectively (Figure 5.1). These events along with their associated aftershock sequence represent a significant percentage of the historic seismicity recorded in the Mackenzie Mountains. Recent work by (Drooff and Freymueller, 2023) which detected microseismicity in the NCC found few small magnitude earthquakes in the Nahanni region and thus inferred that active deformation must be lower than that estimated by (Leonard et al., 2008). In any case, we lack the GPS station coverage to perform a meaningful analysis regarding rates of active deformation in the Nahanni region. We thus choose to limit our Cordillera block boundary to the northeast microseismicity limit detected in (Drooff and Freymueller, 2023) which extends as far east as Tungsten near the Yukon/ Northwest Territory border. Future studies that utilize a more dense GPS station distribution in the Nahanni region will be better posed to address questions about active deformation and recurrence intervals for earthquakes like the 1985 Nahanni sequence. 97 Figure 5.7 Top: GPS site velocities and corresponding 95% confidence ellipses corrected for predictions of the ’joined’ GIA model. Velocities shown are relative to the GEODVEL definition of stable North America. Red dashed line gives location of projection profile shown on bottom. Bottom: GPS site velocities projected onto profile line (shown above). Blue dashed lines give intersection of Denali and Tintina faults with the profile line. 98 5.4 Results East-northeast convergence at a rate of 3-5 mm/yr throughout eastern Alaska and Yukon Territory is observed from the corrected GPS velocity field (Figure 5.7). Site velocities in southeast Alaska are dominantly oriented north-northeast and rotate eastward through the NCC. The Denali fault marks a clear boundary across which observed site velocities rotate from northeast to east-northeast in orientation and decay in magnitude (Figure 5.7). A clear variation in fault-normal motion is not detected across the Tintina fault (Figure 5.7 panel b). While sites outboard of the Tintina fault have northeast-directed velocities of ∼ 4 mm/yr and those inboard of the Tintina fault are closer to ∼ 3 mm/yr this change in velocity trend is within the margin of error for our velocity estimates. Northeast-directed motion across the Mackenzie Mountains decays from ∼ 3 mm/yr within 100 km of the Tintina fault to ∼ 2 mm/yr by the Yukon/Northwest Territory border and then further to less than 1 mm/yr at Norman Wells near the Cordilleran Deformation Front. Site velocities in the Richardson Mountains are moving almost purely eastward, whereas those in the Mackenzie Mountain foreland are oriented east-northeast. 5.4.1 Block Model Solutions We test the four different block model configurations (Figure 5.6) using the corrected site velocities shown in Figure 5.7. For each of the given block configurations we obtain three separate solutions, each of which corresponds to a different set of GIA corrections applied to the input GPS data. The three applied GIA models are described in Section 5.3.3.1 and forward predictions from each of these models are plotted in Figure 5.5. As a significant component of the observed eastward-oriented site velocities in the NCC are a result of the applied horizontal GIA corrections, these model runs are illustrative of uncertainties in our block model solutions that result from the GIA correction. The resulting misfit is represented by the summed Weighted Residual Sum of Squares (WRSS) for each model configuration in Table 5.2. Model configurations which include neither the Tintina nor Richardson block boundaries have a 20 − 25% greater total WRSS 99 60 km lith GIA Model Block Config WRSStot dof 3 ’joined’ 6 6 9 3 6 6 9 3 6 6 9 NT, NR NT, R T, NR T, R NT, NR NT, R T, NR T, R NT, NR NT, R T, NR T, R 604.21 546.84 492.14 484.12 603.95 546.69 491.85 483.90 575.61 529.95 484.78 478.94 vm5a F1 n/a 2.01 5.01 2.31 n/a 2.30 5.01 2.60 n/a 1.90 4.12 2.12 p1 n/a 0.12 0.003 0.019 n/a 0.12 0.002 0.019 n/a 0.20 0.009 0.07 F2 n/a n/a n/a 2.85 n/a n/a n/a 2.72 n/a n/a n/a 2.23 p2 n/a n/a n/a 0.021 n/a n/a n/a 0.026 n/a n/a n/a 0.08 F3 n/a n/a n/a 0.36 n/a n/a n/a 0.34 n/a n/a n/a 0.26 p3 n/a n/a n/a 0.93 n/a n/a n/a 0.93 n/a n/a n/a 0.99 Table 5.2 Misfit Values, F , and p statistics for all tested block model configurations. Block configurations are labeled according to whether they contain separate Tintina (T) or Richard- son (R) blocks. Conversely, No Tintina (NT) or No Richardson (NR) codes signify that the model configuration does not include separate Tinitna or Richardson blocks. dof gives de- grees of freedom for given block model, or the number of free parameters in the inversion. WRSStot gives total summed weighted residual sum of squares misfit for given model. F statistics are computed using Equation 5.2. F1 and p1 give statistical significance values for improvement of model fit for increase in model complexity relative to the base configuration with neither a Tintina block boundary nor a separate Richardson block as shown in panel a of Figure 5.6. F2 and p2 are computed relative to the model with a Richardson block but no Tintina block boundary as shown in panel b of Figure 5.6. F3 and p3 are computed relative to the model with a Tintina block boundary but no separate Richardson block as shown in panel c of Figure 5.6. relative to the models that include both. Elevated misfit values in the less complex models come from the sites that are in the key areas of our study region. Velocities fit to GPS sites on either side of the Tintina fault when considered alone display a 130% increase in misfit when the Tintina fault is not considered as a block boundary. The GPS site velocities on the Richardson block also display a 130% increase in misfit when the Richardson block is not present in the block model, however their contribution to overall WRSS is low. We calculate p and F statistics in order to assess whether the addition of the Tintina and/or Richardson block boundaries represent significant improvement in fit to the data. Summary statistics computed for model results from all block configurations and GIA cor- rections applied to input GPS velocities are reported in Table 5.2. We attribute statistical 100 significance to the improvement of fit between models for any p value that is less than 0.05. GPS site velocities which were corrected for both the ’joined’ and 60 km lithosphere GIA models range between 3-6 mm/yr, whereas GPS velocities corrected for the VM5a GIA model are ∼ 1 − 3 mm/yr. It is thus important to note when interpreting Table 5.2 that a comparison of WRSS misfit values between model solutions resulting from different GIA corrections is not one-for-one as the higher input site velocities are fit by increased angular velocities of blocks. Comparing WRSS misfit values among model runs using the same set of GIA corrections is valid. We find statistically significant p values for all model configurations that include the Tintina fault as a block boundary (Table 5.2) over the base configuration. This is due to the relatively dense set of GPS site velocity observations on either side of the Tintina fault, which helps to constrain relative block motion between the Selwyn and Intermontane blocks. The inclusion of a Richardson block by itself does not result in a statistically significant improvement in fit to the data. This is due to both the small number and sparse distribution of GPS sites in northernmost Yukon Territory. The 3 campaign-style sites that lie on the Richardson block (EPLT, EAGP, RCHR) have timeseries exceeding 15 years in length and display eastward motion at a rate of 3 − 4 mm/yr. The lack of observations in the northern Selwyn Mountains just south of the Richardson block leaves any potential relative motion between the Richardson and Selwyn blocks to be underconstrained. As a result of the sparse distribution of data, rotational parameters obtained in our inversion for the Richardson block have higher degrees of uncertainty compared to the Selwyn and Intermontane blocks. Despite this, model configurations that include both the Richardson block and the Tintina fault which were fit to site velocities corrected for the ’joined’ and 60 km lithosphere GIA models do result in significant improvements over the base configuration (p2 statistics on Table 5.2), but not over the configuration that includes a Selwyn block. 101 Figure 5.8 Block model predictions (top) and residuals (bottom) fit to GPS velocities cor- rected for the ’joined’ GIA model (left) and the vm5a GIA model (right). Residual velocities in bottom panels are shown along with corresponding 95% confidence ellipses. Blue lines give block boundaries for our block model that includes the Tintina fault but does not have a Richardson block. Block boundaries in E and SE Alaska are from (Elliott and Freymueller, 2020b), except are truncated around the perimeter of the model space in order to enforce closure. These block rotation parameters are a priori fixed. 102 Fault E Denali N of Duke River E Denali S of Duke River Tintina NE Selwyn SE (Nahanni) E Richardson S Richard- son* E son* Richard- Selwyn vm5a Strike Slip Rate (mm/yr) -2.0 ’joined’ Strike Slip Rate (mm/yr) -3.0 vm5a Tensile Slip Rate (mm/yr) 0 - -0.6 ’joined’ Tensile Slip Rate (mm/yr) -0.2 - -1.3 σ Strike Slip Rate (mm/yr) 0.1 σ Ten- sile Slip Rate (mm/yr) 0.1 -1.0 - -1.5 -1.0 - -2.0 -0.7 - -3.5 -0.7 - -2.3 0.1 -1.0 0 - 1.2 -0.3 - -0.5 -1.5 0 - 1.7 -0.6 - -0.8 -1.1 - -1.5 -1.2 - -1.4 -1.6 - -1.8 -0.8 - -1.3 -2.0 - -3.0 -3.0 - -4.0 -0.2 0.3 0 -0.5 0.1 -0.4 -1.1 -0.5 - -0.8 -2.5 -0.6 - -1.2 -1.0 -2.5 0.2 0.2 0.2 0.2 0.3 0.3 0.2 0.2 0.2 0.2 0.2 0.3 0.3 Table 5.3 Fault slip rates and associated formal uncertainties along block boundaries. Pos- itive strike slip rates signify left-lateral motion and negative strike slip rates signify right- lateral motion. Positive tensile slip rates signify extension and negative tensile slip rates signify compression/contraction. Rates and uncertainties are only reported for faults that bound the blocks tested in this study. * signifies slip rates from the model configuration that includes the Richardson block. All other fault slip rates come from the model solution that includes the Tintina fault but not the Richardson block. The Eastern Denali fault is formally defined as the Denali fault segment that lies east of the Totschunda fault junction. Our entire study area is located east of the Totschunda fault junction, so we use the term Eastern Denali fault in this table in order to be consistent with prior literature despite it lying in the western portion of our study area. 5.4.2 Fault Slip Rate Estimates The slip rates along defined block boundaries generally decay towards the north and the east. As the magnitude of block rotations are influenced by the choice of GIA correction applied to the input GPS velocities, so too are the estimated fault slip rates. We thus report slip rates for all sets of GIA corrections in Table 5.3 in order to communicate uncertainties in the solutions for slip rates resulting from the choice of applied GIA correction. In most cases the application of VM5a GIA corrections results in lower estimates of fault slip rates compared to the ’joined’ GIA corrections. The exceptions are tensile slip along both the southern East Denali fault and the Tintina fault. This is a product of the lower VM5a 103 GIA predictions for horizontal motion in Yukon Territory contrasting with convergence in southeast Alaska, leading to a higher differential in northeast-oriented site velocities across the Denali and Tintina faults (Figure 5.8). Formal uncertainties are between 0.1-0.3 mm/yr for all faults, which is significantly lower than the uncertainties from the velocity solutions (1-2 mm/yr) and much smaller than the difference in slip rates between models obtained from different GIA corrections (∼ 3 mm/yr). Estimates of both strike slip and tensile slip rates along the Eastern Denali fault vary substantially between the segments north and south of the junction with the Duke River fault. North of the junction with the Duke River fault, slip is concentrated in the strike slip component and is constant along strike at 2 - 3 mm/yr of dextral (right-lateral) motion. Contraction rates along the Eastern Denali fault north of the Duke River junction increase towards the northwest and are 0 - 0.6 or 0.2 - 1.3 mm/yr depending on the GIA corrections. Dextral strike slip motion decays along the Eastern Denali fault decay south of the junction with the Duke River fault from either 2.0 or 1.5 mm/yr depending on GIA correction to 1.0 mm/yr by Taku Inlet near Juneau. Tensile slip along the Eastern Denali fault is at a maximum of 2.3 or 3.5 mm/yr of compression just south of the Duke River junction and decays to the southeast to 0.7 mm/yr of compression near Taku Inlet. Dextral strike slip rates along the Tintina fault are consistent along strike at 1 mm/yr for the VM5a GIA modal and 1.5 mm/yr for the ’joined’ GIA model. Contraction along the Tintina fault increases from 0.8 to 1.3 mm/yr to the northwest for the ’joined’ GIA corrections and 1.1 to 1.5 mm/yr for the VM5a GIA corrections. Contraction along the northeast edge of the Selwyn block ranges from 1 to 3 mm/yr and dextral stike slip across the northeast edge of the Selwyn block ranges from 1.2 to 3.0 mm/yr. Slip is concentrated in the strike slip component where the boundary is oriented mostly east-west and transitions to mostly compression where the boundary is subparallel with the Tintina fault. There is a change in the orientation of strike slip motion associated with the concavity reversal of the block boundary between the northeast and southeast (Nahanni) edges of the Selwyn block. 104 Strike slip rates are lower in the Nahanni region at 0.5 - 1.0 mm/yr and tensile slip increases southward from 1.5 - 4.0 mm/yr. We report slip rates along the eastern edge of the Richardson block for model configu- rations that both do and do not include the Richardson block. For the configuration that does not include the Richardson block the ’joined’ GIA corrections yield dextral strike slip at a low level of 0.5 mm/yr and compression at a significantly higher at 2.5 mm/yr along the eastern extent of the boundary. For the vm5a corrections dextral strike slip along the eastern extent is very low at 0.2 mm/yr and compression is at 1.1 mm/yr. Model configurations that do include the Richardson block predict lower levels (0 - 0.4 mm/yr) of dextral strike slip along the eastern boundary and tensile slip remains relatively constant at 2.5 mm/yr. At the southern end of the block estimates of strike slip are very low at less than 0.3 mm/yr and tensile slip ranges from 0.5 - 1.2 mm/yr depending on the GIA correction. 5.5 Discussion 5.5.1 GPS Site Velocities in the NCC Site velocities that are corrected for the 3D effects of GIA strongly suggest that geodetic strain is accommodated inboard into the NCC at least as far as the Mackenzie Mountain foreland. This observation is consistent with observed patterns in seismicity (i.e. Leonard et al., 2008; Drooff and Freymueller, 2023) and represents the first clear observation of geodetic strain inboard of the Tintina fault. The Tintina fault itself does not appear to be accommodating noticeable geodetic strain as clearly as the Denali fault (Figure 5.7), however we lack the station distribution in the near field of the fault to definitively rule it out. The east-northeast convergence that we observe is in good agreement with both the orientation and magnitude of predicted surface displacements due to mantle traction forces acting at the base of the lithosphere from (McConeghy et al., 2022) and extends well into the Mackenzie Mountain fold-thrust belt. 105 5.5.2 Implications of GIA Model Tests We are unable to distinguish between the examined loading models (ICE-6G, ICE-7G, and ANU-ICE) from the presented 3D GPS site velocities and instead find that differences in observable present day deformation are driven by viscoelastic structure of the solid earth. We find that correcting for the ICE-6G/ VM5a model combination produces the lowest magnitude residual GPS site velocities oriented towards the former location of the Laurentide ice dome in Hudson Bay, an effect which is driven by the VM5a solid earth model. We lack 3D predictions for the other loading models interacting with the VM5a earth model and are thus limited in our ability to interrogate this matter further. The ICE6G/VM5a GIA model was optimized over North America and was further constrained by space geodetic measurements (Roy and Peltier, 2017) so it stands to reason that this would be the best model for LW deglaciation in northern Canada. Residual site velocities after removal of our preferred LW GIA model are dominantly oriented southwards (Figure 5.3), orthogonal to the direction of advance/retreat of LW deglaciation, so we interpret these to not be representing an unaccounted-for aspect of the GIA process. We instead suggest that this could be due to error in the model for North American plate rotation, estimation of ORB, or that it could represent internal deformation of the North American plate interior. If this long wavelength pattern in residual site ve- locities is indeed due to internal deformation then they are not exactly in agreement with those observed by (Kreemer et al., 2018), which are lower in magnitude and directed towards the southeast. We are also able to constrain the viscosity structure of the mantle beneath the western margin of the North American craton. We find that observed vertical GPS site velocities favor a lower mantle viscosity on the order of 2E21 Pa*s and we can definitively rule out a lower mantle viscosity of 2E22 Pa*s. The GIA corrections that we examine in this study have a substantial impact on resulting block model solutions. Our preferred ’joined’ GIA model, which was constructed in order to reflect the dramatic change in lithospheric thickness between the North American craton and 106 the NCC, predicts eastward-oriented velocities between 2-3 mm/yr in the region covered by our block model. The vm5a model, which represents a lower bound on the estimated post- LGM GIA signal, predicts eastward-oriented velocities at a lower ∼ 1 mm/yr, but critically also predicts northward motion at a rate of ∼ 1 mm/yr. This lower eastward-oriented velocity prediction results in lower block rotation and fault slip rates in block model solutions. For all GIA motion corrections the pattern of east-northeast convergence is observed at GPS sites, however the magnitude of this effect is impacted by the effect of lithospheric thickness on the GIA process. 5.5.3 Implications of Block Model Results Our block model results clearly favor the inclusion of the Tintina fault as a block bound- ary. Trends in microseismicity that were recently detected along the Tintina fault by (Drooff and Freymueller, 2023) suggest that the structure is presently active, albeit at a low rate. Historic seismicity rates lead (Leonard et al., 2008) to infer a ∼ 0.5 mm/yr deformation rate along the Tintina fault, which is in general agreement with the ∼ 1.0 − 1.5 mm/yr rate estimated this study. The key difference that we are able to resolve is that this slip is equal in magnitude in both the strike slip and tensile components, which implies that a significant component of northeast-directed strain is accommodated along the Tintina fault or a nearby subsidiary structure. We find that the amount of contraction between the Cordillera and the North American craton along the northeast extent of the Selwyn block is between 1 - 3 mm/yr depending on the GIA correction, which is in general agreement with the 1.8 mm/yr estimated by (Leonard et al., 2008). We find that the convergence rate increases to the south in the Nahanni region to 1.5 - 4 mm/yr. The rate estimated by (Leonard et al., 2008) is 4 mm/yr, which was influenced by the 1985 Nahanni earthquake sequence and thus their estimate likely represents an upper bound on the amount of allowable deformation, similar to the result obtained from our ’joined’ GIA model correction. Intense seismic activity in the Richardson Mountains (Figure 5.1) has been interpreted as compelling evidence for active tectonic deformation (i.e. Leonard et al., 2008; Mazzotti 107 Block Longitude (oE) Latitude (oN) Rate (o/Ma) Ω(X, Y, Z)(10−3 rad/Ma) Richardson 39.87 -42.83 0.067 0.66, 0.55, -0.80 Selwyn 64.07 76.37 0.055 0.09, 0.60, 0.22 Intermontane 81.442 20.33 0.037 0.10, 0.20, 0.93 Covariance Ω (xx,xy,xz,yy,yz,zz) (10−6 rad/Ma2) 1.24, 1.15, -3.87, 1.08, -3.61, 12.1 0.02, 0.02, -0.05, 0.02, -0.05, 0.11 0.01, 0.11, -0.28, 0.12, -0.31, 0.79 Table 5.4 Euler Pole Locations, Rotation Rates, Angular Velocities, and Angular Velocity Covariance For Preferred Block Model. and Hyndman, 2002). Geodynamic models predict a significant component of southward- directed regional strain in the Richardson Mountains that arises due to mantle traction forces acting at the base of the lithosphere (McConeghy et al., 2022). We find that GPS site velocities in the Richardson Mountains have a somewhat higher southward-directed component than those located further south on the Selwyn block. While this data is fit better by separating them onto a separate tectonic block, we do not interpret this to be a significant feature in the model. This is due to the low number of GPS site observations in the Richardson Mountains and the clustering of those sites that are available towards the southeast portion of the theoretical block. Slip rates that we obtain from block models that do not include the Richardson block suggest that slip rates are dominantly oriented in the tensile component, which is at odds with the strike-slip deformation at a rate of 2.1 mm/yr inferred from seismicity (Leonard et al., 2008) and recorded fault offsets (Pinet, 2021b). When the Richardson block is included in the block configuration, the estimate of strike slip rate along the eastern boundary decays further, suggesting that the discrepancy between our models and independent observations cannot be resolved by the addition of the extra model feature. We interpret this to be a result of the sparse regional station distribution and lack of observations directly to the east of the inferred block boundary. Future work that includes the establishment of new GPS sites in northernmost Yukon Territory will be well posed towards addressing potential tectonic block rotation in the Richardson Mountains. 108 5.6 Conclusions The 3D effects of GIA are pervasive in both the North American continental interior as well as the Northern Canadian Cordillera. We test forward predictions of all available models for post-Laurentide deglaciation and find that the ICE6G/VM5a GIA model provides the best fit to GPS site velocities in the North American continental interior. When corrected for the 3D effects of post-Laurentide deglaciation, GPS velocities indicate that measurable tectonic deformation is translated into the Northern Canadian Cordillera at least as far inboard as the Mackenzie Mountain foreland. We fit a regionally optimized tectonic block model to GPS site velocities in order to investigate the number of blocks warranted and to estimate deformation rates along structures such as the Tintina fault, Eastern Denali fault, and the inferred northeast extent of present day deformation. The deformation rates that we obtain are in general agreement with those generated from earthquake statistics and provide an upper bound on the amount of allowable tectonic deformation. The range of solutions that we obtain are illustrative of the effect that GIA corrections have on the resulting tectonic block model solutions. 109 CHAPTER 6 CONCLUSIONS 110 The models and datasets presented in this dissertation are representative of tectonic processes that occur along the same convergent margin but at substantially different stages of the subduction process. In the Aleutian trench the convergence of the Pacific plate with North America is responsible for some of the largest magnitude earthquakes ever recorded. Using site velocities data from onshore GPS sites we presented the most detailed model to date of geodetic coupling along the Alaska-Aleutian margin to date. This model was iteratively updated twice since publication in order to reflect observed rupture in the 2020 Mw 7.8 Simeonof earthquake as well as observations of Vp/Vs variations offshore of the western Kodiak archipelago. Throughout these models we explored the endmember constraints that are able to placed on interseismic coupling due to onshore geodetic data. Future observations at seafloor geodetic sites will be well posed to address questions about how far updip the interseismic coupling distribution extends in this highly segmented region of the megathrust. Longer timescale observations alongside accurate models of postseimic processes will also be well posed to address how long-lived these coupling segments are. Yakutat collision in southern Alaska is a primary driving force behind the neotectonic character of the region. This effect has been extensively studied in the near field, however this dissertation has provided valuable insight into how far northeast into the North American continental interior this process extends. We have been able to model the process extending at least as far inboard as the Mackenzie Mountain foreland, and to also suggest that the Tintina fault is actively accommodating a measurable portion of deformation. Our results are in broad agreement with previous work which has characterized deformation rates using seismicity (Leonard et al., 2008) and has also allowed us to place constraints on the 3D effects of North American plate motion and Glacial Isostatic adjustment. 111 BIBLIOGRAPHY Abers, Geoffrey A et al. (1995). “Source scaling of earthquakes in the Shumagin region, Alaska: time-domain inversions of regional waveforms”. In: Geophysical Journal Inter- national 123.1, pp. 41–58. Altamimi, Zuheir et al. (2011). “ITRF2008: an improved solution of the international terrestrial reference frame”. In: Journal of Geodesy 85.8, pp. 457–473. Altamimi, Zuheir et al. (2017). “ITRF2014 plate motion model”. In: Geophysical Journal International 209.3, pp. 1906–1912. Altamimi, Zuheir et al. (2023). “ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions”. In: Journal of Geodesy 97.5, p. 47. Argus, Donald F (2007). “Defining the translational velocity of the reference frame of Earth”. In: Geophysical Journal International 169.3, pp. 830–838. Argus, Donald F et al. (2010). “The angular velocities of the plates and the velocity of Earth’s centre from space geodesy”. In: Geophysical Journal International 180.3, pp. 913–960. Atwater, Tanya and HW Menard (1970). “Magnetic lineations in the northeast Pacific”. In: Earth and Planetary Science Letters 7.5, pp. 445–450. Audet, Pascal et al. (2016). “Control of lithospheric inheritance on neotectonic activity in northwestern Canada?” In: Geology 44.10, pp. 807–810. Audet, Pascal et al. (2019). “Seismic evidence for lithospheric thinning and heat in the northern Canadian Cordillera”. In: Geophysical Research Letters 46.8, pp. 4249–4257. Audet, Pascal et al. (2020). “Moho variations across the northern Canadian Cordillera”. In: Seismological Research Letters 91.6, pp. 3076–3085. Austermann, Jacqueline et al. (2021). “The effect of lateral variations in Earth structure on Last Interglacial sea level”. In: Geophysical Journal International 227.3, pp. 1938–1960. Bacon, Charles R et al. (2007). “Young cumulate complex beneath Veniaminof caldera, Aleutian arc, dated by zircon in erupted plutonic blocks”. In: Geology 35.6, pp. 491–494. Baker, Michael G et al. (2020). “The Mackenzie Mountains EarthScope Project: Studying active deformation in the northern North American Cordillera from margin to craton”. In: Seismological Research Letters 91.1, pp. 521–532. Berthier, Etienne et al. (2010). “Contribution of Alaskan glaciers to sea-level rise derived from satellite imagery”. In: Nature Geoscience 3.2, pp. 92–95. 112 Bertiger, Willy et al. (2010). “Single receiver phase ambiguity resolution with GPS data”. In: Journal of geodesy 84, pp. 327–337. Biggs, Juliet et al. (Feb. 2009). “The postseismic response to the 2002 M 7.9 Denali Fault earthquake: constraints from InSAR 2003–2005”. In: Geophysical Journal International 176.2, pp. 353–367. issn: 0956-540X. B¨ohm, Johannes et al. (2009). “Forecast Vienna Mapping Functions 1 for real-time analysis of space geodetic observations”. In: Journal of Geodesy 83, pp. 397–401. Bolton, Andrew R et al. (2022). “Evidence for asthenospheric flow rotation in northwest In: Geophysical Journal International insights from shear wave splitting”. Canada: 228.3, pp. 1780–1792. Bond, Gerard C et al. (1985). “An early Cambrian rift to post-rift transition in the Cordillera of western North America”. In: Nature 315.6022, pp. 742–746. Brothers, Daniel S et al. (2018). “Strain partitioning in southeastern Alaska: Is the Chatham Strait fault active?” In: Earth and Planetary Science Letters 481, pp. 362–371. Bruns, Terry R (1983). “Model for the origin of the Yakutat block, an accreting terrane in the northern Gulf of Alaska”. In: Geology 11.12, pp. 718–721. Castillo, David A and William L Ellsworth (1993). “Seismotectonics of the San Andreas fault system between Point Arena and Cape Mendocino in northern California: Implications for the development and evolution of a young transform”. In: Journal of Geophysical Research: Solid Earth 98.B4, pp. 6543–6560. Cecile, MP and DG Cook (1981). Structural Cross-section, Northern Selwyn and Mackenzie Mountains, Yukon and Northwest Territories. Cecile, MP and BS Norford (2000). Geology of the northeastern Niddery Lake map area, east-central Yukon and adjacent Northwest Territories. Vol. 553. Geological Survey of Canada. Cecile, MP et al. (1982). “The Lower Paleozoic Misty Creek Embayment, Selwyn Basin, Yukon and Northwest Territories”. In. Cecile, MP et al. (1997). “Early Paleozoic (Cambrian to early Devonian) tectonic framework, Canadian cordillera”. In: Bulletin of Canadian Petroleum Geology 45.1, pp. 54–74. Clark, Stuart R et al. (2008). “Episodicity in back-arc tectonic regimes”. In: Physics of the Earth and Planetary Interiors 171.1-4, pp. 265–279. 113 Colpron, Maurice et al. (2007). “Northern Cordilleran terranes and their interactions through time”. In: GSA today 17.4/5, p. 4. Colpron, Maurice et al. (2022). “Late Triassic to Jurassic magmatic and tectonic evolution of the Intermontane terranes in Yukon, northern Canadian Cordillera: Transition from arc to syn-collisional magmatism and post-collisional lithospheric delamination”. In: Tectonics 41.2, e2021TC007060. Cook, Frederick A et al. (2004). “Precambrian crust beneath the Mesozoic northern Cana- dian Cordillera discovered by Lithoprobe seismic reflection profiling”. In: Tectonics 23.2. Cross, Ryan S and Jeffrey T Freymueller (2008). “Evidence for and implications of a Bering plate based on geodetic measurements from the Aleutians and western Alaska”. In: Journal of Geophysical Research: Solid Earth 113.B7. Crowell, Brendan W and Diego Melgar (2020). “Slipping the Shumagin Gap: A kinematic coseismic and early afterslip model of the Mw 7.8 Simeonof Island, Alaska, earthquake”. In: Geophysical Research Letters 47.19, e2020GL090308. Davies, John et al. (1981). “Shumagin seismic gap, Alaska Peninsula: History of great earthquakes, tectonic setting, and evidence for high seismic potential”. In: Journal of Geophysical Research: Solid Earth 86.B5, pp. 3821–3855. DeMets, Charles et al. (2010). “Geologically current plate motions”. In: Geophysical journal international 181.1, pp. 1–80. DeSanto, John B et al. (2023). “Limited Shallow Slip for the 2020 Simeonof Earthquake, Alaska, Constrained by GNSS-Acoustic”. In: Geophysical Research Letters 50.16. Ding, Kaihua et al. (2019). “Glacial isostatic adjustment, intraplate strain, and relative sea level changes in the eastern United States”. In: Journal of Geophysical Research: Solid Earth 124.6, pp. 6056–6071. Dixon, James P et al. (2015). 2013 volcanic activity in Alaska: Summary of events and response of the Alaska Volcano Observatory. Tech. rep. US Geological Survey. Drooff, Connor and Jeffrey T Freymueller (2023). “New insights into the active tectonics of the Northern Canadian Cordillera from an enhanced earthquake catalog”. In: Journal of Geophysical Research: Solid Earth 128.12, e2023JB026793. Dusel-Bacon, Cynthia et al. (2002). “Mesozoic thermal history and timing of structural events for the Yukon–Tanana Upland, east-central Alaska: 40Ar/39Ar data from meta- morphic and plutonic rocks”. In: Canadian Journal of Earth Sciences 39.6, pp. 1013– 1051. 114 Dziewonski, Adam M and Don L Anderson (1981). “Preliminary reference Earth model”. In: Physics of the earth and planetary interiors 25.4, pp. 297–356. Eberhart-Phillips, Donna et al. (2006). “Imaging the transition from Aleutian subduction to Yakutat collision in central Alaska, with local earthquakes and active source data”. In: Journal of Geophysical Research: Solid Earth 111.B11. Elliott, Julie and Jeffrey T. Freymueller (2020a). “A Block Model of Present-Day Kinematics of Alaska and Western Canada”. In: Journal of Geophysical Research: Solid Earth 125.7. e2019JB018378 2019JB018378, e2019JB018378. Elliott, Julie and Jeffrey T Freymueller (2020b). “A block model of present-day kinematics of Alaska and western Canada”. In: Journal of Geophysical Research: Solid Earth 125.7, e2019JB018378. Elliott, Julie L et al. (2010). “Tectonic block motion and glacial isostatic adjustment in southeast Alaska and adjacent Canada constrained by GPS measurements”. In: Journal of Geophysical Research: Solid Earth 115.B9. Elliott, Julie L et al. (2022). “Cascading rupture of a megathrust”. In: Science advances 8.18, eabm4131. Engebretson, David C et al. (1985). “Relative motions between oceanic and continental plates in the Pacific basin”. In. Est`eve, C et al. (2020). “The upper mantle structure of northwestern Canada from teleseis- mic body wave tomography”. In: Journal of Geophysical Research: Solid Earth 125.2, e2019JB018837. Est`eve, C et al. (2021). “Surface-Wave Tomography of the Northern Canadian Cordillera Us- ing Earthquake Rayleigh Wave Group Velocities”. In: Journal of Geophysical Research: Solid Earth 126.8, e2021JB021960. Esteve, Cl´ement et al. (2022). “Seismic evidence for a weakened thick crust at the Beaufort Sea continental margin”. In: Geophysical Research Letters 49.16, e2022GL100158. Feng, Wanpeng et al. (2019). “Source parameters of the 2017 M w 6.2 Yukon earthquake doublet inferred from coseismic GPS and ALOS-2 deformation measurements”. In: Geo- physical Journal International 216.3, pp. 1517–1528. Fern´andez-Viejo, Gabriela et al. (2005). “Constraints on the composition of the crust and uppermost mantle in northwestern Canada: Vp/Vs variations along Lithoprobe’s SNor- CLE transect”. In: Canadian Journal of Earth Sciences 42.6, pp. 1205–1222. 115 Finzel, Emily S et al. (2015). “Surface motions and intraplate continental deformation in Alaska driven by mantle flow”. In: Geophysical Research Letters 42.11, pp. 4350–4358. Fletcher, Hilary J et al. (2001). “High interseismic coupling of the Alaska subduction zone In: Geophysical research letters 28.3, SW of Kodiak island inferred from GPS data”. pp. 443–446. Fletcher, Hilary Jane (2002). “Crustal deformation in Alaska measured using the Global Positioning System”. PhD thesis. Fl¨uck, P et al. (2003). “Effective elastic thickness Te of the lithosphere in western Canada”. In: Journal of Geophysical Research: Solid Earth 108.B9. Forte, AM et al. (2010). “Deep-mantle contributions to the surface dynamics of the North American continent”. In: Tectonophysics 481.1-4, pp. 3–15. Fournier, Thomas and Jeff Freymueller (2008). “Inflation detected at Mount Veniaminof, Alaska, with campaign GPS”. in: Geophysical research letters 35.20. Fournier, Thomas J and Jeffrey T Freymueller (2007). “Transition from locked to creeping subduction in the Shumagin region, Alaska”. In: Geophysical Research Letters 34.6. Freed, Andrew M et al. (2006). “Implications of deformation following the 2002 Denali, Alaska, earthquake for postseismic relaxation processes and lithospheric rheology”. In: Journal of Geophysical Research: Solid Earth 111.B1. Freymueller, Jeffrey T and John Beavan (1999). “Absence of strain accumulation in the western Shumagin segment of the Alaska subduction zone”. In: Geophysical Research Letters 26.21, pp. 3233–3236. Freymueller, Jeffrey T et al. (2021). “Constraints on the slip distribution of the 1938 MW 8.3 Alaska peninsula earthquake from tsunami modeling”. In: Geophysical Research Letters 48.9, e2021GL092812. Fritz, WH et al. (1991). “Cambrian to middle Devonian assemblages”. In: Geology of the Cordilleran Orogen in Canada. Edited by H. Gabrielse and CJ Yorath. Geological Survey of Canada, Geology of Canada 4, pp. 151–218. Fu, Yuning et al. (2012). “The effect of using inconsistent ocean tidal loading models on GPS coordinate solutions”. In: Journal of Geodesy 86.6, pp. 409–421. Gabrielse, H and CJ Yorath (1989). “DNAG# 4. The Cordilleran Orogen in Canada”. In: Geoscience Canada. 116 Gabrielse, Hubert (1985). “Major dextral transcurrent displacements along the Northern Rocky Mountain Trench and related lineaments in north-central British Columbia”. In: Geological Society of America Bulletin 96.1, pp. 1–14. Gabrielse, Hubert et al. (2006). “Cretaceous and Cenozoic dextral orogen-parallel displace- ments, magmatism, and paleogeography, north-central Canadian Cordillera”. In: Pa- leogeography of the North American Cordillera: Evidence For and Against Large-Scale Displacements: Geological Association of Canada Special Paper 46, pp. 255–276. Garcia, Emmanuel Soliman M et al. (2019). “Outer trench slope flexure and faulting at Pacific basin subduction zones”. In: Geophysical Journal International 218.1, pp. 708– 728. Garc´ıa, Jos´e Enrique et al. (2022). “Performance of deep learning pickers in routine network processing applications”. In: Seismological Society of America 93.5, pp. 2529–2542. Grand, Steven P (2002). “Mantle shear–wave tomography and the fate of subducted slabs”. In: Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 360.1800, pp. 2475–2491. Hall, KW and Frederick A Cook (1998). “Geophysical transect of the Eagle Plains foldbelt and Richardson Mountains anticlinorium, northwestern Canada”. In: Geological Society of America Bulletin 110.3, pp. 311–325. Hansen, Vicki L et al. (1993). “Asymmetric rift interpretation of the western North American margin”. In: Geology 21.12, pp. 1067–1070. Hayes, Gavin P et al. (2018). “Slab2, a comprehensive subduction zone geometry model”. In: Science 362.6410, pp. 58–61. Hayward, Nathan (2019). “The 3D geophysical investigation of a middle Cretaceous to Paleocene regional d´ecollement in the Cordillera of northern Canada and Alaska”. In: Tectonics 38.1, pp. 307–334. He, X et al. (2018). “The 1 May 2017 British Columbia-Alaska earthquake doublet and implication for complexity near southern end of Denali fault system”. In: Geophysical Research Letters 45.12, pp. 5937–5947. Herman, Matthew W et al. (2018). “The accumulation of slip deficit in subduction zones in the absence of mechanical coupling: implications for the behavior of megathrust earth- quakes”. In: Journal of Geophysical Research: Solid Earth 123.9, pp. 8260–8278. Hetland, EA and M Simons (2010). “Post-seismic and interseismic fault creep II: Tran- sient creep and interseismic stress shadows on megathrusts”. In: Geophysical Journal International 181.1, pp. 99–112. 117 Hreinsd´ottir, Sigr´un et al. (2006). “Coseismic deformation of the 2002 Denali fault earth- quake: Insights from GPS measurements”. In: Journal of Geophysical Research: Solid Earth 111.B3. Hu, Yan and Jeffrey T Freymueller (2019). “Geodetic Observations of Time-Variable Glacial Isostatic Adjustment in Southeast Alaska and Its Implications for Earth Rheology”. In: Journal of Geophysical Research: Solid Earth 124.9, pp. 9870–9889. Hyndman, RD (2010). “The consequences of Canadian Cordillera thermal regime in recent tectonics and elevation: a review”. In: Canadian Journal of Earth Sciences 47.5, pp. 621– 632. Hyndman, Roy D et al. (2005a). “Current tectonics of the northern Canadian Cordillera”. In: Canadian Journal of Earth Sciences 42.6, pp. 1117–1136. Hyndman, Roy D et al. (2005b). “Subduction zone backarcs, mobile belts, and orogenic heat”. In: Gsa Today 15.2, pp. 4–10. Johnson, Kaj M et al. (2009). “Coupled afterslip and viscoelastic flow following the 2002 Denali Fault, Alaska earthquake”. In: Geophysical Journal International 176.3, pp. 670– 682. Johnston, Paul (1993). “The effect of spatially non-uniform water loads on prediction of sea-level change”. In: Geophysical Journal International 114.3, pp. 615–634. Kao, Honn et al. (2012). “Regional centroid-moment-tensor analysis for earthquakes in In: Seismological Research Letters 83.3, Canada and adjacent regions: An update”. p. 505. Klein, Fred W (2002). User’s guide to HYPOINVERSE-2000, a Fortran program to solve for earthquake locations and magnitudes. Tech. rep. US Geological Survey. Koons, Peter O et al. (2010). “Three-dimensional mechanics of Yakutat convergence in the southern Alaskan plate corner”. In: Tectonics 29.4. Kreemer, Corn´e et al. (2018). “A robust estimation of the 3-D intraplate deformation of the North American plate from GPS”. in: Journal of Geophysical Research: Solid Earth 123.5, pp. 4388–4412. Lambeck, Kurt (1993). “Glacial rebound of the British Isles—II. A high-resolution, high- precision model”. In: Geophysical Journal International 115.3, pp. 960–990. Lambeck, Kurt et al. (1998). “Tests of glacial rebound models for Fennoscandinavia based on instrumented sea-and lake-level records”. In: Geophysical Journal International 135.2, pp. 375–387. 118 Lambeck, Kurt et al. (2017). “The North American Late Wisconsin ice sheet and mantle viscosity from glacial rebound analyses”. In: Quaternary Science Reviews 158, pp. 172– 210. Larsen, Christopher F et al. (2005). “Rapid viscoelastic uplift in southeast Alaska caused by post-Little Ice Age glacial retreat”. In: Earth and planetary Science letters 237.3-4, pp. 548–560. Lay, Thorne et al. (2012). “Depth-varying rupture properties of subduction zone megathrust faults”. In: Journal of Geophysical Research: Solid Earth 117.B4. Leonard, Lucinda J et al. (2007). “Current deformation in the northern Canadian Cordillera inferred from GPS measurements”. In: Journal of Geophysical Research: Solid Earth 112.B11. Leonard, Lucinda J et al. (2008). “Deformation rates estimated from earthquakes in the northern Cordillera of Canada and eastern Alaska”. In: Journal of Geophysical Research: Solid Earth 113.B8. Lewis, TJ et al. (2003). “Heat flow, heat generation, and crustal temperatures in the northern Canadian Cordillera: thermal control of tectonics”. In: Journal of Geophysical Research: Solid Earth 108.B6. Li, Shanshan and Jeffrey T Freymueller (2018). “Spatial Variation of Slip Behavior Beneath the Alaska Peninsula Along Alaska-Aleutian Subduction Zone”. In: Geophysical Research Letters 45.8, pp. 3453–3460. Li, Shanshan et al. (2016). “Slow slip events and time-dependent variations in locking In: Journal of beneath Lower Cook Inlet of the Alaska-Aleutian subduction zone”. Geophysical Research: Solid Earth 121.2, pp. 1060–1079. Lu, Zhong and Daniel Dzurisin (2014). “InSAR imaging of Aleutian volcanoes”. In: InSAR imaging of Aleutian volcanoes. Springer, Berlin, Heidelberg, pp. 87–345. Lund, Karen (2008). “Geometry of the Neoproterozoic and Paleozoic rift margin of western Laurentia: Implications for mineral deposit settings”. In: Geosphere 4.2, pp. 429–444. Ma, Shutian and Pascal Audet (2017). “Seismic velocity model of the crust in the northern In: Canadian Journal of Canadian Cordillera from Rayleigh wave dispersion data”. Earth Sciences 54.2, pp. 163–172. Mair, John L et al. (2006). “Deformation history of the northwestern Selwyn Basin, Yukon, Canada: Implications for orogen evolution and mid-Cretaceous magmatism”. In: Geo- logical Society of America Bulletin 118.3-4, pp. 304–323. 119 Mao, Ailin et al. (1999). “Noise in GPS coordinate time series”. In: Journal of Geophysical Research: Solid Earth 104.B2, pp. 2797–2816. Marechal, Ana¨ıs et al. (2015). “Indentor-corner tectonics in the Yakutat-St. Elias collision constrained by GPS”. in: Journal of Geophysical Research: Solid Earth 120.5, pp. 3897– 3908. Marone, Chris and CH Scholz (1988). “The depth of seismic faulting and the upper transition from stable to unstable slip regimes”. In: Geophysical Research Letters 15.6, pp. 621– 624. Marsman, Celine P et al. (2021). “The impact of a 3-D Earth structure on glacial isostatic adjustment in Southeast Alaska following the Little Ice Age”. In: Journal of Geophysical Research: Solid Earth 126.12, e2021JB022312. Matmon, Ari et al. (2006). “Denali fault slip rates and Holocene–late Pleistocene kinematics of central Alaska”. In: Geology 34.8, pp. 645–648. Maus, S et al. (2009). “EMAG2: A 2–arc min resolution Earth Magnetic Anomaly Grid com- piled from satellite, airborne, and marine magnetic measurements”. In: Geochemistry, Geophysics, Geosystems 10.8. Mazzotti, St´ephane and Roy D Hyndman (2002). “Yakutat collision and strain transfer across the northern Canadian Cordillera”. In: Geology 30.6, pp. 495–498. McCaffrey, Robert et al. (2002). “Crustal block rotations and plate coupling”. In: Plate Boundary Zones, Geodyn. Ser 30, pp. 101–122. McCann, WR et al. (1979). “Seismic gaps and plate tectonics: Seismic potential for major boundaries”. In: Earthquake prediction and seismicity patterns, pp. 1082–1147. McConeghy, Joseph et al. (2022). “Investigating the Effect of Mantle Flow and Viscos- ity Structure on Surface Velocities in Alaska Using 3-D Geodynamic Models”. In: Journal of Geophysical Research: Solid Earth 127.10. e2022JB024704 2022JB024704, e2022JB024704. Meade, Brendan J and Bradford H Hager (2005). “Block models of crustal motion in south- ern California constrained by GPS measurements”. In: Journal of Geophysical Research: Solid Earth 110.B3. Meade, Brendan J and John P Loveless (2009). “Block modeling with connected fault- network geometries and a linear elastic coupling estimator in spherical coordinates”. In: Bulletin of the Seismological Society of America 99.6, pp. 3124–3139. 120 Mogi, Kiyoo (1958). “Relations between the eruptions of various volcanoes and the defor- mations of the ground surfaces around them”. In: Earthq Res Inst 36, pp. 99–134. Monger, Jim and Ray Price (2002). “The Canadian Cordillera: geology and tectonic evolu- tion”. In: CSEG Recorder 27.2, pp. 17–36. Monger, JWH (1997). “Plate tectonics and northern Cordilleran geology: an unfinished revolution”. In: Geoscience Canada. Morgan, W Jason (1968). “Rises, trenches, great faults, and crustal blocks”. In: Journal of Geophysical Research 73.6, pp. 1959–1982. Mousavi, S Mostafa et al. (2020). “Earthquake transformer—an attentive deep-learning model for simultaneous earthquake detection and phase picking”. In: Nature communi- cations 11.1, pp. 1–12. M¨unchmeyer, Jannes et al. (2022). “Which picker fits my data? A quantitative evaluation of deep learning based seismic pickers”. In: Journal of Geophysical Research: Solid Earth, e2021JB023499. Murphy, DC and JK Mortensen (2003). “Late Paleozoic and Mesozoic features constrain displacement on Tintina Fault and limit large-scale orogen-parallel displacement in the Northern Cordillera”. In: Program with Abstracts. Vol. 28. Nishenko, SP and KH Jacob (1990). “Seismic potential of the Queen Charlotte-Alaska- Aleutian seismic zone”. In: Journal of Geophysical Research: Solid Earth 95.B3, pp. 2511– 2532. Nuttli, Otto W (1973). “Seismic wave attenuation and magnitude relations for eastern North America”. In: Journal of Geophysical Research 78.5, pp. 876–885. Okada, Yoshimitsu (1992). “Internal deformation due to shear and tensile faults in a half- space”. In: Bulletin of the seismological society of America 82.2, pp. 1018–1040. Oldow, John S et al. (1990). “Transpression, orogenic float, and lithospheric balance”. In: Geology 18.10, pp. 991–994. Peltier, W Richard et al. (2015). “Space geodesy constrains ice age terminal deglaciation: The global ICE-6G C (VM5a) model”. In: Journal of Geophysical Research: Solid Earth 120.1, pp. 450–487. Peltier, WR (1985). “The LAGEOS constraint on deep mantle viscosity: results from a new normal mode method for the inversion of viscoelastic relaxation spectra”. In: Journal of Geophysical Research: Solid Earth 90.B11, pp. 9411–9421. 121 Peltier, WR (1998). ““Implicit ice” in the global theory of glacial isostatic adjustment”. In: Geophysical Research Letters 25.21, pp. 3955–3958. Pinet, Nicolas (2021a). “Structural geology of the eastern Richardson Mountains, Yukon and Northwest Territories: Some field observations and a note of caution for palinspastic reconstructions”. In. — (2021b). “Structural geology of the eastern richardson mountains, yukon and northwest territories: Some field observations and a note of caution for palinspastic reconstructions”. In: Yukon Exploration and Geology. Plafker, G et al. (1994). “Neotectonic map of Alaska, scale 1: 2,500,000”. In: Geology of Alaska, Map, GNAG-1, Plate 12a, Geol. Soc. Am., Boulder, Colo. Purcell, Anthony et al. (2016). “An assessment of the ICE6G C (VM5a) glacial isostatic adjustment model”. In: Journal of Geophysical Research: Solid Earth 121.5, pp. 3939– 3950. Ristau, John et al. (2007). “Stress in western Canada from regional moment tensor analysis”. In: Canadian Journal of Earth Sciences 44.2, pp. 127–148. Ross, GM (1991). “Tectonic setting of the Windermere Supergroup revisited”. In: Geology 19.11, pp. 1125–1128. Ross, Zachary E et al. (2019). “PhaseLink: A deep learning approach to seismic phase association”. In: Journal of Geophysical Research: Solid Earth 124.1, pp. 856–869. Roy, Keven and William R Peltier (2018). “Relative sea level in the Western Mediterranean basin: a regional test of the ICE-7G NA (VM7) model and a constraint on late Holocene Antarctic deglaciation”. In: Quaternary Science Reviews 183, pp. 76–87. Roy, Keven and WR Peltier (2017). “Space-geodetic and water level gauge constraints on continental uplift and tilting over North America: regional convergence of the ICE-6G C (VM5a/VM6) models”. In: Geophysical Journal International 210.2, pp. 1115–1142. Ryan, James J et al. (2017). “Landscape antiquity and Cenozoic drainage development of southern Yukon, through restoration modeling of the Tintina Fault”. In: Canadian Journal of Earth Sciences 54.10, pp. 1085–1100. Schaeffer, AJ and S Lebedev (2014). “Imaging the North American continent using waveform inversion of global and USArray data”. In: Earth and Planetary Science Letters 402, pp. 26–41. Schutt, Derek L et al. (2023). “Lithospheric S wave velocity variations beneath the Mackenzie In: Journal of Geophysical Research: Mountains and Northern Canadian Cordillera”. 122 Solid Earth 128.1, e2022JB025517. Shillington, Donna J et al. (2015). “Link between plate fabric, hydration and subduction zone seismicity in Alaska”. In: Nature Geoscience 8.12, pp. 961–964. Sigloch, Karin and Mitchell G Mihalynuk (2013). “Intra-oceanic subduction shaped the assembly of Cordilleran North America”. In: Nature 496.7443, pp. 50–56. Snyder, David B et al. (2005). “Contrasting seismic characteristics of three major faults in northwestern Canada”. In: Canadian Journal of Earth Sciences 42.6, pp. 1223–1237. Spada, Giorgio and P Stocchi (2007). “SELEN: A Fortran 90 program for solving the “sea- level equation””. In: Computers & Geosciences 33.4, pp. 538–562. Steffen, H (2021). Surface Deformations from Glacial Isostatic Adjustment Models with Laterally Homogeneous, Compressible Earth Structure. Suito, Hisashi and Jeffrey T Freymueller (2009). “A viscoelastic and afterslip postseismic deformation model for the 1964 Alaska earthquake”. In: Journal of Geophysical Research: Solid Earth 114.B11. Sykes, Lynn R (1971). “Aftershock zones of great earthquakes, seismicity gaps, and earth- In: Journal of Geophysical Research quake prediction for Alaska and the Aleutians”. 76.32, pp. 8021–8041. Takagi, Ryota et al. (2019). “Along-strike variation and migration of long-term slow slip events in the western Nankai subduction zone, Japan”. In: Journal of Geophysical Re- search: Solid Earth 124.4, pp. 3853–3880. Tape, Carl and Anthony Lomax (2022). “Aftershock regions of Aleutian-Alaska megath- rust earthquakes, 1938–2021”. In: Journal of Geophysical Research: Solid Earth 127.7, e2022JB024336. Thomas, Ian D et al. (2011). “Widespread low rates of Antarctic glacial isostatic adjustment revealed by GPS observations”. In: Geophysical research letters 38.22. Thorkelson, Derek J and Richard P Taylor (1989). “Cordilleran slab windows”. In: Geology 17.9, pp. 833–836. Tushingham, AiMi and WR Peltier (1991). “Ice-3G: A new global model of late Pleistocene deglaciation based upon geophysical predictions of post-glacial relative sea level change”. In: Journal of Geophysical Research: Solid Earth 96.B3, pp. 4497–4523. Van Der Lee, Suzan and Guust Nolet (1997). “Seismic image of the subducted trailing fragments of the Farallon plate”. In: Nature 386.6622, pp. 266–269. 123 Wal, Wouter van der et al. (2015). “Effect of GIA models with 3D composite mantle viscosity on GRACE mass balance estimates for Antarctica”. In: Earth and Planetary Science Letters 414, pp. 134–143. Waldhauser, Felix (2001). “hypoDD—A program to compute double-difference hypocenter locations (hypoDD version 1.0-03/2001)”. In: US Geol. Surv. Open File Rep., 01 113. Wang, Fan et al. (2022). “Subduction fluids control slab slip behaviors and megathrust In: AGU Fall Meeting Abstracts. Vol. 2022, earthquakes at the Alaska Peninsula”. T32C–07. Witter, Robert C et al. (2014). “Little late Holocene strain accumulation and release on the Aleutian megathrust below the Shumagin Islands, Alaska”. In: Geophysical Research Letters 41.7, pp. 2359–2367. Xiao, Zhuohui et al. (2021). “The deep Shumagin gap filled: Kinematic rupture model and slip budget analysis of the 2020 Mw 7.8 Simeonof earthquake constrained by GNSS, global seismic waveforms, and floating InSAR”. in: Earth and Planetary Science Letters 576, p. 117241. Yuan, Huaiyu and Barbara Romanowicz (2010). “Lithospheric layering in the North Amer- ican craton”. In: Nature 466.7310, pp. 1063–1068. Zhang, Haijiang and Clifford Thurber (2003). “User’s manual for tomoDD1. 1 (double- difference tomography) for Determining Event Locations and Velocity Structure from Local Earthquakes and Explosions”. In: Department of Geology and Geophysics, Uni- versity of Wisconsin-Madison, Madison, WI. 124