Influence of Farallon Slab Loading on Intraplate Stress and Seismicity in Eastern North America in the Presence of Pre‐Existing Weak Zones

Despite the stability of the continental interior, eastern North America has hosted many significant historical earthquakes. Seismicity concentrates within tectonically inherited structures, which can act as weak zones where stress accumulates. Within these zones, systematic stress rotations may be explained by long‐wavelength sources. We test the hypothesis that mantle‐flow driven by the Farallon slab contributes to intraplate seismicity via the reactivation of pre‐existing faults. We model the stress field using seismically constrained global high‐resolution finite‐element flow models with CitcomS. To isolate the slab's effect, we vary its buoyancy between a case of neutrality and a case with full negative thermal buoyancy derived from tomography. Low‐viscosity lithospheric weak zones located at failed rifts, loaded by a mass anomaly at depth, transmit elevated stresses to the overlying crust. The sinking of the Farallon slab drives localized mantle flow beneath the central‐eastern US, generating a large stress amplification of 100–150 MPa peaking over the New Madrid Seismic Zone (NMSZ). This stress amplification exerts a continent‐wide clockwise rotation on the stress field, which in the presence of weak zones reproduces some observed deviations of the seismically inferred SHmax from the regional borehole SHmax, bringing optimally oriented faults, closer to failure, some of which are associated with major historical earthquakes, including the Reelfoot Fault in the NMSZ and the Timiskaming Fault in Western Quebec. However, stronger lithospheric viscosity gradients, shallower weak zones, or weaker faults are still needed to fully reproduce the observed stress field in some areas.


Introduction
Intraplate seismicity is an enigmatic phenomenon not directly associated with plate boundary tectonics, potentially resulting from a combination of geodynamic processes, and surface loading, with contribution from far-field tectonic forces and pre-existing structures.The continental interior of eastern North America in particular has hosted many significant historical earthquakes and is undergoing both modern day glacial isostatic adjustment (GIA) (Grollimund & Zoback, 2001;Wu & Hasegawa, 1996;Wu & Johnston, 2000) and long-wavelength dynamic subsidence (Forte et al., 2010;Spasojevic et al., 2008).Historic records of major earthquakes in eastern North America date back to at least 1638 (Figure 1a), with most seismicity occurring within clearly defined seismic zones (Figure 1b).Of these zones, the New Madrid, Central Virginia, Charleston South Carolina, Eastern Tennessee, Northern Appalachian, Western Quebec, Charlevoix, and St. Lawrence River Valley Seismic Zones have experienced destructive M 5 to M ≥ 7 earthquakes in recorded history (Figure 1a, Table S1 in Supporting Information S1), most notably the 1811-1812 M 7-8 (Hough et al., 2000;S. Stein, 2007) New Madrid earthquakes (Figure 1a).While these events are infrequent, earthquake hazard remains significant for the centraleastern United States (CEUS) due to limited earthquake preparedness.The M 5.8 Mineral, Virginia earthquake in 2011 and the recent M 4.8 Tewksbury, New Jersey earthquake in April 2024 serve as reminders that intraplate seismicity remains a concern for the continental interior.Moreover, earthquake shaking is potentially just as hazardous as that on plate margins because the stable continental lithosphere can transmit seismic energy more efficiently, and peak ground acceleration may be almost twice that of plate boundary zones for a given earthquake magnitude (S.Stein, 2007).While intraplate seismicity in eastern North America is a well studied phenomena, many unknowns remain about the cause of these earthquakes, in particular the role of long wavelength loading and mantle dynamics; this paper aims to address these unknowns.
Seismicity is even more enigmatic given the intrinsically low strain rates characteristic of intraplate regions and that there is little to no surface expression of active faults (Mazzotti, 2007).For example, earthquake statistics alone suggest a few kilometers of total seismic deformation over the last million years (Mazzotti & Adams, 2005), but seismic reflection surveys of both the Charlevoix and New Madrid seismic zones indicate total accumulated deformation on particular structures to be no more than a few tens of meters (Lamontagne & Ranalli, 1996;Schweig & Ellis, 1994;Van Arsdale, 2000).However, such low strain rates, typically on the order of 10 9 yr 1 or less based on GPS (Calais et al., 2006;Ghosh et al., 2019;Mazzotti & Adams, 2005;Mazzotti et al., 2005), do not preclude the possibility of future large earthquakes.Mechanical models from Kenner and Segall (2000) demonstrate that relaxation of local lithospheric weak zones can transfer stress to the upper crust and trigger slip on faults, resulting in an earthquake sequence that continues until the weak zone fully relaxes, which can be prolonged by cyclic stress transfer from co-and post-seismic slip reloading the lower crust (DiCaprio et al., 2008;Kenner & Segall, 2000;Kenner & Simons, 2005).Such weak zone relaxation and stress cycling can trigger earthquake sequences with large slip events every 500-1,000 years with shear strain rates on the order of 10 7 yr 1 to 5.5 × 10 9 yr 1 , despite surface deformation rates often having been below the GPS detection threshold (Kenner & Segall, 2000).
Likewise, mechanical modeling of the 1811-1812 New Madrid earthquakes has revealed that intraplate seismic zones tend to stay in a stress shadow for hundreds to even thousands of years after large events, owing to the difficulty of full stress restoration due to the strength of the ambient crust and low regional strain rates (Li et al., 2007;S. Stein & Liu, 2009).Such findings imply seismicity in this region may be non-stationary, with mainshocks clustered in time (DiCaprio et al., 2008;Kenner & Simons, 2005; S. Stein & Liu, 2009).This would suggest the New Madrid seismicity is transient in nature and that seismic hazard in the mid-continent may be overestimated (Newman et al., 1999).This idea is supported in part by paleoseismic evidence demonstrating an elevated slip rate on the Reelfoot Fault during the two most recent major earthquake cycles (Van Arsdale, 2000).
In general, intraplate earthquakes tend to be episodic and clustered (S.Stein & Liu, 2009), and seismicity may migrate between similar structures due to long-term deformation (S.Stein, 2007), which may result from regional stress changes due to epeirogenic subsidence (Spasojevic et al., 2008), GIA (Mazzotti et al., 2005;Steffen, Steffen, et al., 2014;Wu & Johnston, 2000), or even denudation and erosion (Calais et al., 2010;Craig et al., 2017).In light of possible stress cycling, however, one must consider whether hazard is highest where seismicity has recently been concentrated or uniform within regions with similar structures.
While local structures are key to controlling intraplate seismicity, one of the more compelling lines of evidence suggestive of a long-wavelength source of stress are rotations in the stress field.Focal mechanism (FM) stress inversion reveals stress orientations mostly in agreement with the regional NE-SW shortening direction, but some aulacogens exhibit rotational deviation, suggesting multiple stress regimes (Hurd & Zoback, 2012;Mazzotti & Townend, 2010;Verdecchia et al., 2022).Within the Central Virginia, Lower St. Lawrence, and Charlevoix seismic zones, there is as much as a 30°-50°statistically significant clockwise rotation of the seismically derived direction of maximum horizontal compressive stress (S Hmax ) relative to the regional borehole-derived S Hmax orientation (Mazzotti & Townend, 2010).Significant depth-dependent rotations of up to 40°-60°are observed within the Charlevoix Seismic Zone as well, which are argued to result from a combination of weak fault zone crust and post-glacial rebound stress (Verdecchia et al., 2022).Similar but slightly smaller rotations are observed in the Northern Appalachian and New Madrid Seismic Zones.These rotations require stress perturbations at midseismogenic depths of at least 160-250 MPa (Mazzotti & Townend, 2010).Mazzotti and Townend (2010) argue that the consistency of these stress rotations over spatial scales of 100-1,000 s of km and the associated seismicity cannot be explained solely by phenomena such as complex fault intersections (Gangopadhyay & Talwani, 2007;Talwani, 1988Talwani, , 1999)), crustal density anomalies (Ghosh et al., 2009;Levandowski et al., 2017), or flexure under local loads (S.Stein et al., 1989).On the other hand, long-wavelength sources like dynamic topography, lithospheric flexure, and GIA, which typically only induce stress perturbations on the order of 10-100 MPa, may result in additional stress amplifications by a factor of 5-10 in the presence of a lithospheric "weak zone" (Grollimund & Zoback, 2001;Wu & Mazzotti, 2007), leading to seismic activity in an otherwise stable continental interior.There is thus cause for considering the role of large-scale continental uplift or subsidence (i.e., epeirogeny) arising from mantle dynamics in augmenting stress within intraplate regions.This idea has been explored in part by some authors, but has yet to be explored thoroughly with the consideration of local-scale intraplate weak zones.Assuming a low-viscosity sub-lithospheric upper mantle, Forte et al. (2010) demonstrated how mantle flow due to the sinking of the Farallon slab can create significant surface depression and localize bending stresses within the NMSZ possibly capable of triggering earthquakes (Forte et al., 2007).Saxena et al. (2021) assess how both gravitational potential energy (GPE) and mantle flow from a putative lithospheric drip influence stresses and seismicity in eastern North America using tomography-constrained regional geodynamic models.Their results suggest upper mantle flow alone, without lithospheric heterogeneity, is not sufficient to reproduce the observed intraplate stress field (Saxena et al., 2021).However, without far-field compressive tectonic stresses or crustal weak zones in their models, their predicted S Hmax directions are somewhat inconsistent with the general NE-SW trending S Hmax from the World Stress Map (WSM) (Heidbach et al., 2018) (Figure 1b), and their predicted regional faulting style is opposite to that obtained from FM stress inversion across the CEUS and parts of Canada (Hurd & Zoback, 2012).Ghosh et al. (2019), on the other hand, addressed the influence of strong lateral viscosity variations between the continental craton and the accreted terranes of the eastern margin and find that such a strength contrast results in a broad stress rotation from NE-SW in the interior to roughly E-W along the coastal margin and in the NMSZ, consistent with the earlier works of M. L. Zoback and Zoback (1980) and the results of FM stress inversion (Ghosh et al., 2019;Levandowski et al., 2018).
At a more local scale, Zhan et al. (2016) do investigate the role of smaller pre-existing weaknesses and different degrees of mantle heterogeneity in their regional models of the NMSZ, where they find that a large low-viscosity mantle region beneath the NMSZ can reload the layer above it, concentrating stresses in the upper mantle through displacement.The strength of the lower crust then enables elevated stresses to transmit to the base of the ancient rift, potentially reactivating faults.Crustal or lithospheric loading resulting from displacements in a mantle lowvelocity zone are likely insignificant under far-field stresses alone, however, but could be catalyzed to a significant degree by unloading processes like GIA or even changes in drainage networks (Zhan et al., 2016).For example, correlation between micro-earthquakes in the NMSZ and annual and multi-annual hydrologic loading cycles in the Mississippi embayment demonstrate that long-wavelength changes in continental water-storage can produce observable crustal deformation and modulate the regional seismicity (Craig et al., 2017).Such observations illustrate the degree to which small changes in loading from long-wavelength processes can promote localization and enhancement of crustal stresses that impact seismicity rates.
These studies highlight the importance of local weak zones in controlling intraplate seismicity, as stresses from tectonics, topography, or density contrasts are typically too small to trigger intraplate failure on their own.The integrated strength of a lithospheric section under a cold geotherm like that of eastern North America is estimated to be 5-10 times larger than plate driving forces (Mazzotti, 2007).Thus, because cratons are too cold and stable to deform easily, even localized weak-zones may require very low viscosities in the lower crust and/or upper mantle (10 19 Pa s or less) to achieve seismic strain rates.Such weakness can be due either to conditions that reduce the strength and viscosity of the entire lithosphere, such as elevated heat flow, or conditions that allow for high stress and strain rates in a weak deforming layer that loads the brittle upper crust.A weak zone may also result from a mechanically weak crust, such as a wet quartz-diorite with low effective viscosity; high pore fluid pressure; or an exceptionally weak fault zone with low-friction fault gouge material (Mazzotti, 2007).
Despite the work reported to date, the possible connections between intraplate seismicity and epeirogeny remain elusive.Geodynamic modeling has demonstrated that sinking slabs, in regions of both active and ancient subduction, create mantle downwelling that generates tractions within the overlying lithosphere, resulting in dynamic subsidence (Forte et al., 2010;Gurnis, 1992;Hager, 1984;Mitrovica et al., 1989;Yang et al., 2016).Initially, sinking slabs can stagnate at the 660 km discontinuity due to the viscosity jump and the phase transition between ringwoodite and bridgmanite (Billen, 2008).The negative buoyancy of the accumulated mass is eventually great enough to break through the boundary and continue sinking, initiating a process called slab avalanche, which causes substantial downward asthenospheric flow (Christensen & Yuen, 1984;Tackley et al., 1993;Yang et al., 2016Yang et al., , 2018)).Dynamic subsidence has been inferred in eastern North America from Cenozoic shorelines (Spasojevic et al., 2008), under which the remnant Farallon slab lies between 410 and 2,200 km depth (Lu et al., 2019;Ren et al., 2007).The presence of the slab predominantly beneath the 660 km discontinuity implies possible slab avalanche, as do models of mantle flow localized beneath the eastern US (Forte et al., 2007), but questions remain as to how this flow-induced continental subsidence affects stress in the lithosphere.Moreover, Geochemistry, Geophysics, Geosystems 10.1029/2024GC011493 the increasingly compressive stress regime moving northeastwards into the Canadian margin (Hurd & Zoback, 2012) is consistent with stress patterns observed above deep slabs far inland of the trench, such as in Sundaland (Yang & Gurnis, 2016;Yang et al., 2018).While this epeirogenic movement may not cause notable folding or faulting (Gurnis, 1992), it may reactivate pre-existing faults.
We explore the hypothesis that mantle flow caused by the sinking of the Farallon slab contributes to intraplate seismicity in eastern North America via stress perturbations and reactivation of pre-existing faults.We develop high-resolution global geodynamic flow models with CitcomS (Moresi et al., 2014;Tan et al., 2006;Zhong et al., 2000), a spherical finite-element thermochemical mantle convection code.The models use an initial thermal structure constrained by seismic tomography and geology so as to accurately reflect lithospheric and mantle structure, adequately capture the Farallon slab at depth, and reproduce observed plate motions.As we are interested in the current state of stress, we calculate the instantaneous flow field and associated stress tensor, with which we compute the S Hmax direction and stress magnitudes in the crust and evaluate our results against the WSM Database (Heidbach et al., 2018) and the data of Mazzotti and Townend (2010).Our work improves upon previous investigations of this topic by (a) using a global, rather than regional scale, convection model based on recent seismic tomography to naturally incorporate far-field tectonic forcing, (b) purposefully testing the impact of the Farallon slab on the stress field by parameterizing its negative buoyancy, (c) explicitly including localscale, low-viscosity lithospheric weak-zones at the locations of the geologically mapped aulacogens, and (d) resolving the modeled CFS onto known faults within the seismic zones of eastern North America.

New Madrid Seismic Zone
The NMSZ is predominantly associated with the Reelfoot Rift-part of a system of intracratonic faults formed from late-stage crustal extension during and following the breakup of super-continent Rodinia and the opening of the Iapetus ocean (Thomas & Powell, 2017).At its northeast end, it merges with the eastward trending Rough Creek Graben, which serves as a strike-slip offset between the Reelfoot Rift and the Rome Trough (Figure 1).The basement of the NMSZ has been deformed by both extensional fracturing during rifting and contraction and reactivation of those fractures within the Reelfoot Rift during the Appalachian-Ouachita orogeny and the assembly of Pangea (Thomas & Powell, 2017).Such repeated deformation and reactivation could result in a weakened lower crustal that helps concentrate seismicity in the NMSZ.However, a high-velocity lower crust, often interpreted as a mafic pillow, is inconsistent with a weakened low-viscosity rheology at these depths.Rather, low seismic velocity zones are identified in the upper mantle beneath the NMSZ, extending from the base of the crust to as deep as 300 km (Chen et al., 2016;Thomas & Powell, 2017).The passage of a mantle plume is also often used to explain the crustal structure and associated concentration of seismicity in the NMSZ (Chu et al., 2013).A combination of thermal variations, compositional heterogeneity, and water content are needed to explain the observed variations in seismic velocity beneath the NMSZ and the Illinois Basin.Rifting and/or the passage of a plume would have introduced iron, water, and basalt into the host peridotite of the region, refertilizing the lithospheric mantle, both reducing its seismic velocities and creating large strength differences between the rift and surrounding region (Chen et al., 2016).The partial melting and mafic material that intruded the lower crust can explain the high velocity layer at the base of the crust below the Moho and above ∼70 km.The weakening of olivine from hydration and greater iron content would result in a low-viscosity weak zone in the upper mantle beneath the rift, which by creeping faster than surrounding areas can transfer stress and load the crust (Kenner & Segall, 2000;Thomas & Powell, 2017).

Eastern Tennessee and Central Virginia Seismic Zones
Unlike the NMSZ, the ETSZ is not associated with an aulacogen or other rift arm but instead with the New York-Alabama (NY-AL) magnetic lineament, which formed as a sinistral continental transform fault between proto-Laurentia and Amazonia during the assembly of supercontinent Rodinia and the Grenville orogeny (Thomas & Powell, 2017).Seismic anisotropy from SKS shear wave splitting indicates a strong shearing component parallel to the lineament (Long et al., 2016;Thomas & Powell, 2017;Wagner et al., 2012).A prominent low velocity zone exists in the basement beneath the ETSZ at depths of 5-25 km, consistent with a sheared crust (Powell et al., 2014) consisting of by soft phyllosilicate minerals like talc, phlogopite, and antigorite (Chen et al., 2016).Such minerals are inherently weak and lower in viscosity due to their layered structure and exhibit significantly reduced shear strengths.In eastern Tennessee, this sheared crust is interpreted to be a weak zone bounded by en echelon conjugate shears and normal faults beneath the NY-AL lineament (Thomas & Powell, 2017).Magnetotelluric imaging of this region (Murphy & Egbert, 2017) also reveals a contrast in lithospheric electrical properties between the Appalachian highlands and the Piedmont and Coastal Plain lowlands to the east; the former, near the ETSZ, being underlain by more conductive material between 25 and 55 km depth.This change in properties between the two domains can be interpreted as a viscosity contrast ranging anywhere from 1 to 6 orders of magnitude (Murphy & Egbert, 2017;Murphy et al., 2019).The NY-AL lineament itself is a narrow feature, at most 25-50 km wide (Thomas & Powell, 2017).The deeper portion of this potential low viscosity structure, however, ranges from 25 to 125 km wide (Murphy & Egbert, 2017).For this reason, we represent the weak zone associated with the ETSZ by a roughly 50-100 km wide low-viscosity zone between 15 and 50 km depth, with narrower widths at shallower depths, which allows us to capture at least two elements width in our models.
MT imaging has also captured the lithospheric structure north of the ETSZ in the vicinity of the Rome Trough and the Central Virginia Seismic Zone (CVSZ).Highly conductive asthenosphere beneath the central Appalachians imaged by the MAGIC array (Evans et al., 2019;Long et al., 2020) reveals thin lithosphere on the order of 80 km thick.Regions of high conductivity extending from 80 km to over 200 km depth across the entire width of the Appalachian Mountains are consistent with earlier shear velocity models showing a pronounced 300 km wide low velocity zone known as the Central Appalachian Anomaly (Schmandt & Lin, 2014).Such an anomaly could be caused by extensive hydration of the mantle due to fluids derived from the relict Farallon slab and/or high concentrations of conductive solid phases, such as sulphides and graphite indicative of metasedimentary shear fabrics developed during continental suturing (Evans et al., 2019), which would lower mantle viscosity.Beneath the Rome Trough, conductivity is sufficiently high to be consistent with some degree of partial melt and hence lower viscosity (Evans et al., 2019).Thus, we include a low viscosity weak zone between 80 and 200 km depth parallel to the Rome Tough and approximately 500 km in width.

Lower Saint Lawrence River, Charlevoix, and Western Quebec Seismic Zones
The Lower Saint Lawrence Rift System (LSLRS) as a whole contains all of the Lower Saint Lawrence Rift, the Ottawa-Bonnechere Graben, and the Saguenay Graben, as well as the Western Quebec, Charlevoix, and Lower Saint Lawrence Seismic Zones (Figure 1).Throughout the LSLRS, earthquakes predominantly occur within the Precambrian Shield between depths of 5-25 km on steeply dipping normal faults that were formed during the opening of the Iapetus Ocean (Lamontagne & Ranalli, 2014).The Lower Saint Lawrence Seismic Zone is the least seismically active of those in the LSLRS.The Charlevoix Seismic Zone (CXSZ), on the other hand, has the highest seismic hazard of continental eastern Canada (Lamontagne & Ranalli, 2014).In the CXSZ, earthquakes cluster along or between the mapped Iapetan rift faults of the St. Lawrence Rift (Lamontagne & Brouillette, 2022), but the CXSZ is also the site of an ancient impact structure.It has been suggested that this comparatively weak damaged volume has lower elastic moduli that influence the stability of the rift faults that intersect it, promoting localization of low level seismicity within the crater and larger events near the perimeter and outside the crater (Baird et al., 2010;Lamontagne & Brouillette, 2022;Thomas & Powell, 2017).The rift faults themselves are believed to have lower Mohr-Coulomb frictional strength (Baird et al., 2010), meaning the predominant source of weakness may lie in the crust, not the mantle.Moreover, the CXSZ exhibits depthdependent crustal stress rotation and strength variations, primarily between 13 and 26 km depth (Verdecchia et al., 2022).The seismicity occurs between 0 and 26 km depth, with the largest stress rotations occurring between 12 and 16 km, yielding a total clockwise rotation of about 30°between shallow and deeper focal mechanisms (Verdecchia et al., 2022).The largest rotation from the regional borehole direction (up to 60°) is at mid-crustal depths between 20 and 26 km.Similar clockwise rotations of 44°-49°are reported for the Lower St. Lawrence Seismic Zone (LSLSZ) north of the CXSZ (Verdecchia et al., 2022).
The WQSZ is likewise associated with Iapetan rift arms and aulacogens, specifically the Ottawa-Bonnechere and Timiskaming Grabens, which are characterized by both NW-SE and NE-SW striking, steeply dipping faults (Rimando & Peace, 2021).The majority of earthquakes in the WQSZ occur at 8-18 km depth and follow a linear NW-SE trend slightly adjacent but parallel to the Ottawa-Bonnechere Graben (Figure 1).In addition to aulacogens, it has been suggested that the passage of the Great Meteor Hotspot during the Mesozoic or an extension of the New England Seamount Chain track, whose path aligns with the orientation of the WQSZ, enabled thermomechanical weakening of ancient faults and emplaced strength contrasts between felsic rocks and mafic intrusions in the mid-crust (Ma & Eaton, 2007;Sykes, 1978).The hotspot theory is supported in part by lithospheric velocity anomalies at around 200 km depth (Ma & Eaton, 2007), but there is no surface geological expression of a hotspot (Lamontagne & Ranalli, 2014).The WQSZ also exhibits thinner crust than the surrounding Grenville Province, which typically ranges from 30 to 48 km thick (Ma & Eaton, 2007).Despite local structural complexities that no doubt play a role in controlling seismicity in the LSLRS, the ancient rift itself likely provides a source of lower crustal or upper mantle weakness.In our models, we position weak zones at the locations of the LSLRS and its rift arms between 25 and 75 km depth.

Geodynamic Modeling in CitcomS
We use the 3D spherical shell finite-element mantle convection code CitcomS to solve the conservation equations of mass and momentum (Moresi et al., 2014;Tan et al., 2006;Zhong et al., 2000Zhong et al., , 2008)).Under the Boussinesq approximation for an incompressible fluid, the conservation equations are: where u is the velocity in the non-rotating mantle reference frame; ε is the strain-rate; p is the dynamic pressure; η, ρ 0 , and g 0 are mantle viscosity, reference density, and gravitational acceleration, respectively; êr is the unit vector in the radial direction.δg is a perturbation to radial gravitational acceleration that if implemented would represent the effect of self-gravitation, and δρ represents perturbations to density, which can arise from both thermal anomalies and composition.
Here, α is the coefficient of thermal expansion, T is the temperature at a given point, T r is the radial mean temperature at a given layer, C is the composition field variable representing compositional and/or rheological heterogeneity and is either 0 or 1, and δρ ch is the density difference between different compositions, which is implemented non-dimensionally using the buoyancy number (Equation 7).It is these density perturbations that give rise to mantle flow.
The stress field in which we are interested arises from these buoyancy forces and is determined by the first two terms in Equation 2. The stress, σ, is inherently dependent on changes in viscosity, η, and strain-rate, ε, through the constitutive relationship and thus-by the definition of strain rate-on the gradients in the velocity of the flow.A dense mass anomaly at depth (e.g., the Farallon slab) will sink and excite downward flow.The presence of the thick high viscosity cratonic lithosphere can also obstruct horizontal asthenospheric flow, driving it downwards (Paul et al., 2023).Such diversion of the flow amplifies the vertical velocity component, meaning the horizontal gradient of the vertical velocity increases approaching the primary source of the downwelling, be it the Farallon slab or the edge of the continent.The increase in this velocity gradient, and thus strain-rate, would induce large tractions on the lithosphere in these regions, resulting in high stress, particularly when scaled by the high viscosity of the lithosphere.Likewise, the presence of intra-lithospheric weak zones introduces velocity gradients between the surrounding rigid high viscosity lithosphere and the faster flowing material within the weak zone, which lead to tractions on the overlying crust.Thus, the presence of lateral gradients in viscosity are key to controlling changes in the strength of the flow that ultimately govern the magnitude and pattern of stress.
Within CitcomS, these equations of motion are solved non-dimensionally and the buoyancy is scaled by the Rayleigh number, Ra (Equation 6), which controls the vigor of convection.In CitcomS, the equations are nondimensionalized using the full planetary radius R instead of the convective layer thickness, so the Rayleigh number is defined using the Earth's radius.η o is the reference viscosity, κ is the thermal diffusivity, and ΔT is the temperature drop from the surface to the core-mantle boundary (CMB) and is the constant used to non-dimensionalize the temperature field.B is the buoyancy ratio (Equation 7), which gives the relative strength between compositional and thermal buoyancy.
CitcomS solves Equations 1 and 2 in a fully global spherical shell geometry with non-dimensional (i.e., scaled by planetary radius) inner and outer radii of r b = 0.55 and r t = 1, respectively.Free-slip boundary conditions and thermal boundary conditions of 0 and 1 (non-dimensional temperature) are applied at the top and bottom boundaries, respectively.Realistic Earth parameters are used for other dimensionalization constants (Table 1).To solve the equations of motion, specification of an initial temperature field is required, which we constrain by seismic tomography.
CitcomS has robust solvers that incorporate temperature, pressure, and composition dependent variable viscosity that can be computed directly from our thermal input (Zhong et al., 2000).Temperature-dependent viscosity is given in non-dimensional form as: where η r (r) is the depth-dependent viscosity prefactor, E is the non-dimensional activation energy that controls the strength of the temperature-dependence of the viscosity, T o is the temperature at which the viscosity equals 1 in the model (set to be that of the upper mantle), and η c is an optional spatially dependent pre-factor that allows for assigning regions of higher viscosity lithosphere or localized weak zones.The mantle is radially stratified in viscosity, with η r defining an average lithospheric layer (0-100 km), the asthenosphere (100-300 km), the upper mantle (300-670 km), and the lower mantle (670-2,870 km).The actual lithosphere, however, is prescribed using the η c prefactor in order to accommodate laterally variable lithospheric thickness, allowing us to better capture regions of thin and thick lithosphere, such as near mid-ocean ridges or under cratonic roots, respectively.
CitcomS makes use of a full multi-grid solver to efficiently solve the discretized matrix equations.The spherical domain is divided into 12 caps, each of which are further divided into N × N × N nearly equal area elements over the surface of the sphere.In order to resolve the narrow rift structures that serve as weak-zones, we need spatial resolutions on the order of at least 50 km.We use meshes with 257 × 257 × 257 elements per cap, which is a ∼25 km horizontal resolution.With radial grid refinement, we can achieve 1 km vertical resolution in the shallowest layers, which is essential for resolving crustal and lithospheric boundaries and incorporating the effect of GPE differences from the crustal density field.

Lithospheric Thermal Structure
In order to incorporate realistic thermal structure and the negative buoyancy of the Farallon slab, we use geologically and seismically constrained thermal input in our models.Global models allow us to fully capture the effect of Farallon slab-driven mantle flow beneath eastern North America and naturally incorporate far-field stresses associated with plate tectonics.This also accounts for the effect of subduction and changes in lithospheric thickness and composition associated with the ocean-continent boundary (Humphreys & Coblentz, 2007;Saxena et al., 2023).We account for the forces exerted by ridges and slabs by using the thermal plate model of the oceanic lithosphere (Turcotte & Schubert, 2014) derived from the seafloor age grid of Seton et al. (2020).We also superimpose a slab thermal structure with geometries constrained by the Slab2.0model (Hayes et al., 2018) and constructed using plate age at the trench and convergence velocity such that temperatures consistent with plate cooling are propagated down the depth of the slab (Hu et al., 2022;Rudi et al., 2022).For the continents, we use the Artemieva (2006) TC1 thermal model of the continental lithosphere, which is constrained primarily by heat flow data, as well as xenolith geotherms and electrical conductivity data, and for which the lithosphere-asthenosphere boundary (LAB) depth is defined as the 1,300°C isotherm.The plate cooling calculation is performed using 1,300°C as the basal temperature of the plate for consistency with the Artemieva ( 2006) model.From the oceanic and continental temperature fields, we obtain both a combined thermal model for the lithosphere and a combined LAB depth, above which we can apply a rigid lithosphere (Figure S1 in Supporting Information S1).

Seismically Constrained Mantle Thermal Structure
The thermal structure of mantle convection models is often constrained by seismic tomography, assuming the seismic anomalies are thermal in origin.It is common to simply scale the velocity to an effective temperature, a method that several studies have used to constrain the temperature and density structure of the mantle beneath eastern North America (Liu et al., 2008;Spasojevic et al., 2009).The more rigorous approach is to invert for temperature using the full anharmonic and anelastic components of the seismic wave-speed equations and the appropriate elastic moduli constrained by mineral physics experimental data or theoretical values for a given composition (Cammarano & Guerri, 2017;Cammarano et al., 2003;Goes, 2002;Goes et al., 2000;Karato, 1993).However, due to the limited resolution of seismic tomography, the increase in compositional influence on seismic velocity with depth, and the uncertainty on the physical properties of deep mantle minerals (Cammarano et al., 2003), inversion tends to become unstable for the lower mantle and drastically over-predict temperatures in the deep earth.Because of this, we adopt the approach of scaling the velocities to temperature using a numerically derived depth dependent scaling factor that, like the inversion technique, incorporates the effect of anelasticity, which is essential because it significantly reduces the sensitivity of shear velocity to temperature.
The scaling factor of velocity to temperature is given as the inverse of the derivative of the natural log of velocity with respect to temperature, where V is either V p or V s .The reader is referred to the Supporting Information S1 for the derivation of the terms in Equation 9and the details of the mineralogical constraints and procedure used to calculate the velocity derivatives and scaling factor.The computed ∂ ln(V)/∂T values (Figure S5c in Supporting Information S1) agree well with those of previous authors (e.g., Steinberger & Calderwood, 2006), and the magnitude of the scaling factors is on the order of those used by Spasojevic et al. (2009) and Liu et al. (2008), who find that a scaling of about 2 × 10 3 °C/ km/s produces flow models consistent with plate motions, stratigraphy, and the history of Farallon subduction beneath North America since the Late Cretaceous.
For the global mantle tomography, we use the TX2019 slab model from Lu et al. (2019).By including a priori 3D slab structure defined by seismicity, this model addresses shortcomings in using seismic tomography to infer density or thermal anomalies that arise from discrepancies between detailed studies of slabs and global tomography models (Lu et al., 2019).In this model, the locations of slabs associated with modern subduction systems are defined by seismicity on a 0.1°× 0.1°grid (Lu et al., 2019).The Farallon slab is identifiable as an elongate high S-wave velocity anomaly in both P and S between 600 and 2,200 km in the mid-lower mantle (Lu et al., 2019;Ren et al., 2007) (Figure 2).P-wave anomalies are determined relative to AK135 (Kennett et al., 1995), and S-wave anomalies are relative to TNA-SNA (Grand & Helmberger, 1984).We first convert all velocity anomalies to absolute velocity using their respective reference profiles, then recompute δV p,s with respect to PREM (Dziewonski & Anderson, 1981) for consistency with the derived scaling factors (see Supporting Information S1).The velocity anomaly at each point in the model is multiplied by the scaling factor corresponding to its depth to yield a temperature anomaly for that point.Absolute temperatures are determined by adding the calculated temperature anomalies to the mantle geotherm (Figure S5d in Supporting Information S1).Temperatures are determined from V P and V S separately and then averaged to get the final temperature anomaly.
The central-eastern North American lithosphere is dominated by high velocities, whereas the eastern Appalachians are characterized by largely average to low velocities, including two strong localized low velocity anomalies in the north-east Appalachians at around 100 km depth (Boyce et al., 2019;Schmandt & Lin, 2014).The high velocities of the North American shield terminate abruptly against the Grenville Front to the east, beyond which low velocities dominate in the Grenville Province.Lithospheric thickness is not observed to vary greatly in this region, but there are some abrupt variations from north to south, which are consistent with previous seismic studies of the LAB in eastern North America (Artemieva, 2006;Hopper & Fischer, 2018).These sharp variations could be explained by metasomatic processes that modified the lithospheric mantle composition during subduction along the Laurentian margin (Boyce et al., 2019) and which could have played a role in weakening the crust and lithosphere.

Combined Temperature and Viscosity Input for CitcomS
The lithospheric thermal model (Figure S1 in Supporting Information S1) is combined with the temperatures obtained from the tomography (Figure 2) by means of a depth-weighted average following a hyperbolic tangent function in order to ensure a smooth temperature field and eliminate any artifacts that may result from local mismatches in temperature estimates.The lithosphere is weighted 100% down to 50 km, which is the starting depth of the tomography model.The weight on the lithospheric model then varies smoothly from about 80% at 50 km depth to near 0% at 240 km depth.Below this, the model is based on the tomography.Temperatures are non-dimensionalized for input to CitcomS assuming a CMB temperature of 4,000 K (Table 1).
Viscosities are computed directly from the temperature field within CitcomS (Equation 8) and then multiplied by an appropriate pre-factor depending on their radius or their location within the lithosphere.Our radial viscosity structure is consistent with both the Haskell constraint (Haskell, 1935), which places the upper mantle at 10 21 Pa s, and the results of previous geodynamic studies that find a 30× viscosity jump at the top of the lower mantle is required to obtain optimal fits to the geoid (Mao & Zhong, 2021;Thoraval & Richards, 1997).We also apply a low viscosity asthenosphere in order to better reproduce plate motions (Mao & Zhong, 2021;Steinberger & Calderwood, 2006).To achieve rigid plates, we also impose 10× viscosity everywhere above the LAB (Figure S1 in Supporting Information S1) in addition to that derived from the low temperatures in the shallow regions of the thermal model.The radial profiles of temperature and viscosity are shown in Figure 3 along with the full range of temperature and viscosity in each layer and the 2σ spread of the values in that layer.The viscosity profile used in CitcomS, with its respective layer-dependent pre-factors, is shown along with other viscosity profiles from the literature in Figure 3b.Lithospheric viscosity structure, including intraplate weak zones, is shown in Figure S8 in Supporting Information S1.

Weak Zones
We explicitly emplace weak zones in the crust and/or lithosphere at the locations of the paleo-rifted structures or seismically inferred low velocity zones in eastern North America.This is achieved by specifying a viscosity prefactor applied at nodal points that fall within the bounds of these structures.Weak zones are placed at an appropriate depth for each seismic zone, based on local seismic and geologic evidence, the details of which are given in Section 2. To summarize, beneath the NMSZ, we emplace a weak zone in the upper mantle between 70 and 150 km along the extent of the Reelfoot Rift as mapped in Figure 1a.In the ETSZ, we emplace a 50-100 km wide weak zone along the NY-AL lineament between 15 and 50 km depth.A bit farther north, along the Rome Trough, we emplace an approximately 500 km wide weak zone between 80 and 200 km depth.Along the entirety of the LSLRS and its associated rift-arms, we emplace a weak zone at the base of the crust between 25 and 75 km depth (Figure S8 in Supporting Information S1).While the viscosities of these weak zones likely differ based on their different conditions of formation, we apply a uniform viscosity reduction to all weak zones for simplicity.
We also test a case where all weak zones are considered shallow and sub-crustal, placed between 25 and 75 km depth.Plate margins are also included as narrow low-viscosity zones, the locations of which are based on the plate margins from Seton et al. (2012).The geometries of weak zones along the subduction interface are also constrained by the slab dip and depth from Slab2.0 (Hayes et al., 2018).

Parameterization of Farallon Slab Buoyancy
The Farallon slab is a significant low temperature anomaly in the mantle beneath eastern North America; it lies mostly between 750 and 2,200 km depth and is believed to be up to 400 km thick in places (Sun et al., 2017).This structure and its negative buoyancy have the potential to induce strong flow in the mantle beneath eastern North America and exert a significant influence on the intraplate stress field.However, buoyancy from the thermal anomaly alone may be overestimated, as compositional buoyancy within the ancient slab may counteract the thermally derived density.The slab's long history of flat-slab subduction under western and central North America during the Laramide Orogeny and period of the Cretaceous Interior Seaway and the passage of the conjugates to the Hess and Shatksy Rises (Sun et al., 2017)-buoyant oceanic plateaus inferred to have been attached to the Farallon plate-indicate the presence of compositional buoyancy.
We explore the degree to which the flow induced by the sinking slab contributes to intraplate stress by parameterizing the buoyancy of the slab into two end-member cases: (a) the slab retains its full negative thermal buoyancy as derived entirely from the tomography (0% of applied buoyancy ratio), and (b) the slab anomaly is entirely neutralized by applying a buoyancy ratio that completely counteracts the thermal buoyancy (100% of the applied buoyancy ratio).Figure S6 in Supporting Information S1 shows a representative longitudinal profile through the effective temperature field (temperature minus the radial mean) and the corresponding profile of buoyancy ratio derived from these effective temperatures.We identify the slab anomaly as that portion of the slab where the temperature anomaly exceeds the 1σ variation of the temperatures at that depth and apply the neutralizing buoyancy ratio only within that portion (i.e., within the filled portion of the temperature profiles in Figure S6a in Supporting Information S1).
Using the same method of calculating the buoyancy ratio from the temperature field, we also apply a neutralizing buoyancy to both the continental lithosphere and the large-low-shear-velocity-provinces (LLSVPs) in the lower mantle.The continental lithosphere, particularly very thick, old cratons, and continental keels, are cold and hence thermally dense; however, cratonic stability over geologic times necessitates a chemically buoyant tectosphere depleted in basaltic components that balances the mass increase due to cold temperatures (Forte & Perry, 2000;Jordan, 1988;Shapiro et al., 1999).Indeed, in models without positive buoyancy applied to the continental lithosphere, the keels start to sink, dramatically altering the global flow field.In a similar fashion, the LLSVPs, while hot and hence thermally buoyant, are believed to be chemically dense and thus neutrally buoyant stable piles at the bottom of the mantle (Ballmer et al., 2016;Ishii & Tromp, 1999;Lau et al., 2017;Vilella et al., 2021).Representative depth slices of the temperature and chemical buoyancy ratio, as applied to the lithosphere, LLSVPs, and Farallon slab, are depicted in Figure S7 in Supporting Information S1.In all models, we apply the full buoyancy ratio required to neutralize the sinking and rising of the lithosphere and LLSVPs, respectively.

Calculation of S Hmax and Coulomb Stress
From the stress tensors computed with CitcomS, we calculate the second invariant of the deviatoric stress, the principal stresses, their directions, and the direction of maximum horizontal compressive stress (S Hmax ) following Lund and Townend (2007), the latter of which is necessary for properly comparing our results to stress data.Accurate calculation of S Hmax , as opposed to using just the principal stress directions, is important because when none of the three principal stress directions are truly vertical, S Hmax does not simply coincide with the largest subhorizontal stress, though it may be an adequate proxy.However, in our predicted stresses, S Hmax , S Hmin , and S V almost always coincide with the principal stresses at crustal depths.We compare estimates of the S Hmax direction with those of the WSM (Heidbach et al., 2018) and other FM stress inversions (Hurd & Zoback, 2012;Mazzotti & Townend, 2010).
Before calculating S Hmax or resolving the stress onto faults, we also add back the lithostatic pressure, which affects the magnitude of the principal stresses and of the normal stress resolved on any given fault plane.Accurately computing such normal stress is essential for computing the CFS on a fault; without the lithostatic component, Coulomb stresses are greatly overestimated.We compute the lithostatic pressure assuming a mean continental crustal density of 2,700 kg/m 3 .For consistency with the theory of the critically stressed crust (Townend & Zoback, 2000; M. D. Zoback & Townend, 2001; M. D. Zoback et al., 2002), the lithostatic pressure is counteracted by the pore-fluid pressure exerted by water and other fluids in the crust such that optimally oriented faults are always on the verge of slipping.Without this reduction in normal stress, lithostatic pressure at seismogenic depth is such that all faults are always far within the stable regime with respect to the Coulomb Failure Criterion (CFC).We use a pore fluid factor of 0.6, which is consistent with intermediate values previously used for conditions of fault slip in intraplate regions like Canada (Rimando & Peace, 2021; M. L. Zoback, 1992).
To better clarify the degree to which Farallon slab-induced mantle flow and/or intraplate weak zones impact the potential for fault reactivation, we must determine whether such stresses are sufficient in both magnitude and orientation to bring intraplate faults closer to failure.The most commonly used criterion for determining whether there will be slip on a fault is the CFC (King et al., 1994), defined by the shear and normal stress on the failure plane, the pore pressure, the coefficient of friction, and the cohesion of the rock (Equation 10).When the shear stress on a fault exceeds a certain threshold, failure may occur, possibly initiating an earthquake (King et al., 1994).The Coulomb failure stress (CFS) (Equation 11) is the difference between the shear stress on a given fault and the failure threshold, where positive Coulomb failure stresses indicate unstable faults (i.e., likely to slip) and negative Coulomb failure stresses indicate stable faults (i.e., unlikely to slip).As is common practice for preexisting faults, we take the cohesion C to be zero, as it is usually only considered when dealing with newly forming faults where the rock actually needs to break.
Because CFS is calculated from the shear and normal stress resolved on an actual fault plane, it is highly dependent on the orientation of that fault (e.g., its strike and dip) with respect to the orientation of the stress tensor.Thus, to assess how intraplate faults might respond to the stress states produced by the models, it is essential to analyze the CFS on realistic faults within the different seismic zones.Using fault orientation data and our computed stress tensors, we compute the CFS on faults in the NMSZ, WQSZ, and CXSZ/LSLRS.The reader is referred to the Supporting Information S1 for more details on the calculation of S Hmax and CFS.

Results
Two main suites of models were computed to explore the influence of low viscosity weak zones on the intraplate stress field of eastern North America: one with the full expression of the Farallon slab in the thermal field as resolved from seismic tomography and the other with the Farallon slab neutralized by means of the buoyancy ratio (Table 2).Within each of these sets, we test a case without weak zones; two cases with weak zones, one with the weak zones at different depths depending on geophysical constraints for their respective locations and one with all weak zones placed at the same depth between 25 and 75 km; and a case without weak zones but in which we have incorporated the effect of GPE due to crustal density and topography variations (see Supporting Information S1).In all cases with weak zones, weak zone viscosity is reduced by a factor of 10 4 .Globally, to first order, the velocity field predicted by the flow models is consistent with global plate motions, and strain rates match observed values from the Global Strain Rate Map (Kreemer et al., 2014) to within 8%-10% error on average.In models with the Farallon slab, strain rates within continental eastern North America are on the order of 1 × 10 9 to 1 × 10 8 yr 1 , and the percent error of the log of the second invariant of deviatoric strain rate between the modeled values and those of Kreemer et al. (2018) is generally less than 10% but can be as high as 20% in the central-eastern US (Figures S9a, S9c, and S9e in Supporting Information S1).In models without the Farallon slab, strain rate is slightly lower, with isolate patches in the central-eastern US reaching 15% error, but the majority of errors remain less than 10% (Figures S9b, S9d, and S9f in Supporting Information S1).The S Hmax directions derived from the CitcomS stress tensors reproduce first order stress provinces globally, including the general NE-SW S Hmax orientation across eastern North America and the radial pattern of stress orientations in China.Across all models, the majority of misfits between the modeled S Hmax orientations and those of the WSM are less than 25°throughout eastern North America (Figure 4a).Misfits are lowest in the central-eastern and north-eastern US in the vicinity of eastern Tennessee and the Rome Trough and up through the Northern Appalachians.Strong fits are also observed in New Madrid and the Canadian seismic zones.The largest misfits are generally isolated to individual points and are not unexpected given the uniformity of the modeled stress field compared to the high frequency spatial variability in the observed S Hmax directions, even across small areas.The presence or lack of either weak zones or the Farallon slab introduces perturbations to the modeled S Hmax direction (Figure 5) and its misfit to the observed data (Figures 4b-4d).Including only intraplate weak zones without the Farallon slab decreases the misfit by only a couple of degrees on average within the major seismic zones (Figure 4, Figure S10 in Supporting Information S1).On the other hand, the inclusion of the Farallon slab alone, without weak zones, changes the S Hmax misfit by up to ±20°(Figures 4c and 4d).Specifically, the fit improves by ∼11°in the NMSZ, ∼17°in CVSZ, ∼13°in the WQSZ, ∼5°in the CXSZ, and ∼13°in the LSLSZ.However, fit worsens by about 8°i n the ETSZ and in South Carolina and by as much as 15°along the Rome Trough (Figure 4 and Figure S10 in Supporting Information S1).Inclusion of both weak zones and the Farallon slab results in further misfit reductions in the NMSZ, CVSZ, LSLSZ, and both on and offshore Nova Scotia, but only by up to 8°.It is important to note that while the Farallon slab has a notable impact on S Hmax orientation, the very small change induced by the inclusion of weak zones, generally less than ±5°, is largely insignificant with respect to the local 90% confidence intervals of the WSM data, which for a single seismic zone often span as much as 25°, with the smallest interval being 8°(Figure S10 in Supporting Information S1).
Overall, the presence of the Farallon slab more accurately reproduces the observed stress field both regionally and in each of the seismic zones, barring the ETSZ.What is more interesting is the direction of rotation induced by the slab effect and the additional small changes induced by weak zones in the presence of the slab's influence.Under only tectonic background forcing, S Hmax directions do generally follow the regional NE-SW trend, but inclusion of the Farallon slab creates a continent wide clockwise rotation of the S Hmax direction by about 10°-15°( Figure 5c).This rotation becomes even greater in the presence of weak zones, ranging from an extra couple of degrees in New Madrid to more than 5°along the Appalachian front, from the ETSZ all the way up through the LSLRS (Figures 5b and 5d).However, weak zones induce counterclockwise rotations in the western part of the WQSZ and in northern New England, as well as on and offshore Nova Scotia (Figure 5b).The largest clockwise stress rotations for models with both the Farallon and weak zones occur in the Central Virginia, Charlevoix, and Lower Saint Lawrence Seismic Zones (Figure 6).
While changes in the S Hmax direction between the different cases are subtle, changes in stress magnitude are substantial.At long-wavelength, the Farallon slab induces a high amplitude stress anomaly across the centraleastern US, the peak of which coincides with the vicinity of the NMSZ and the nearby ETSZ.This peak in stress is largely eliminated in the case without the Farallon slab (dark blue lines in Figures 7a and 7b).Without the slab, mantle flow beneath eastern North America and hence the traction imparted to the base of the lithosphere is greatly reduced (Figures 7e and 7f), resulting in total stress magnitudes on the order of 50-100 MPa and stress perturbations from the inclusion of weak zones on the order of 5-20 MPa (Figure 9).However, when the full expression of the Farallon slab is included (light blue lines in Figures 7a-7d), stress magnitudes jump considerably to 200-250 MPa.Most notably, stress perturbations within the WQSZ, LSLSZ, CXSZ, and ETSZ due to the inclusion of shallow weak zones in the presence of mantle flow (Figures 9a and 9d) are as much as 70-100 MPa-nearly 50 MPa more than in the case of including weak zones alone (Figure 9b).The reason for this is likewise illustrated in Figure 8, which depicts how the Farallon slab both excites and changes the direction of flow (Figures 8a-8d) beneath eastern North America relative to cases without the slab (Figures 8e-8h), localizing  mantle flow beneath the NMSZ.This effect is even more pronounced in the presence of weak zones, particularly those near the base of the lithosphere, where the strong downward flow from the slab propagates into the nearby low viscosity region above it, creating strong localized flow within the weak zone itself (Figures 8c and 8d).
As such, stress magnitudes exhibit a strong dependence on weak zone depth.Models A1 and B1 contain weak zones emplaced at different depths based on where local geophysical data suggests a weak zone may be present (Section 2).In this case, the weak zone in the vicinity of New Madrid and the Rome Trough/West Virginia are deeper than those in eastern Tennessee and southeast Canada.While these weak zones exhibit enhanced flow in response to the forcing of the Farallon slab (Figures 8d and 8c), they are also farther from the surface and thus have a weaker and more diffuse effect on loading the overlying crust.Weak zones that are located near the base of the crust and extend into the upper mantle, between 25 and 75 km depth, but that do not extend to the base of the lithosphere show significantly less flow excitation in response to Farallon loading than those deeper down (Figure 8d), but because they are in contact with the base of the crust, the tractions the weak zone imparts are transmitted directly into the shallow layers, resulting in higher near-surface stresses.In the model with shallow weak zones everywhere (model A1b), stress perturbations in New Madrid become the largest of any of the seismic zones (Figures 7a and 7b).
While differences in stress magnitude between the different model cases demonstrate the pronounced influence of the slab and weak zones, whether those stress perturbations are such that they increase the likelihood of seismicity in eastern North America is dependent on the local fault geometry.To assess the likelihood of fault reactivation for each of the different cases, we resolve CFS on actual faults in the NMSZ, WQSZ, and CXSZ/LSLSZ (Figure 10).Coulomb failure stress is computed for all faults given their strike and dip and the stress tensors along that fault.
In the NMSZ, most faults are steeply dipping, between 70°and 90°, yielding low Coulomb stress values and thus stable faults, but the Reelfoot Fault dips between 30°and 44° (Csontos & Van Arsdale, 2008), placing it in an optimal orientation to be reactivated as a thrust fault in the local stress field.This fault indeed happens to be the site of one of the earthquake epicenters in the 1811-1812 earthquake sequence (Delano et al., 2021).The Reelfoot Fault is also the only fault to exhibit positive Coulomb stress values in each of our models (Figures 10a, 10d, 10g, 10j, and 10m), indicating a higher likelihood of fault reactivation.Similarly, the New Madrid North Fault also shows higher Coulomb stress, though still in the stable regime, and is also the site of one of the historical earthquakes.The NNW striking Joiner Ridge fault, which intersects the Axial Fault, has an unconstrained dip and thus was assigned a dip of 90°.However, this fault is identified as a thrust fault (Delano et al., 2021).If one were to assume a dip between 30°and 45°as in the case of the Reelfoot Fault, this fault would also exhibit high Coulomb stresses, and is in the vicinity of the remaining historical event just to the north on the Axial Fault.While the vertical Axial Fault itself is not optimally oriented and exhibits low Coulomb stress, the intersection of a nearby reverse fault on which high Coulomb stress is localized could explain the occurrence of a major earthquake in this location (Talwani, 1988(Talwani, , 1999)).
The modeled stress state yields results consistent with expectations from the local rupture history, but more interesting are the differences in Coulomb stress we observe between the different models.For non-optimally oriented vertical faults, there is not a pronounced change between different cases, but for the more optimally oriented Reelfoot and New Madrid North Faults, there is an obvious difference.Inclusion of the Farallon slab pushes the Reelfoot Fault farther into the unstable regime by about 75-100 MPa and the New Madrid North Fault closer to failure by about 10-20 MPa.In both cases with and without the Farallon slab, the inclusion of a deep weak zone does not have a notable effect, and actually shifts most faults farther from failure, except for the Reelfoot and New Madrid North Faults.This smaller effect is likely due to the deeper depth of the weak zone.
Without the Farallon slab and without weak zones (Figure 10k), all faults are stable with respect to the Mohr- Coulomb failure criterion assuming a coefficient of friction of 0.4-0.6 and fault dips of 56°, consistent with steeply dipping rift-bounding normal faults of the region (A.L. Bent et al., 2003;Rimando & Peace, 2021).The dips of the Timiskaming fault (northwestern-most purple star in Figures 10b, 10e, 10h, 10k, and 10n) and nearby parallel faults are 45°based on A. L. Bent (1996).Assuming a μ = 0.6, including weak zones without the Farallon slab brings most faults closer to failure (Figure 10n), but not to failure, with the largest changes occurring in and to the northeast of Montreal.If faults are weaker (μ ≈ 0.4), many faults would just pass into the unstable regime, including the four plotted as Mohr circles in Figure S11 in Supporting Information S1.Similar to the NMSZ, influence of the Farallon slab brings all faults closer to failure by at least 15-25 MPa, but only a few faults exhibit high Coulomb stress above the failure criterion-a change again on the order of 75-100 MPa (Figures 10b, 10e,  and 10h).These faults all have similar orientations and dips and include the fault on which the 1935 M 6.1  Timiskaming earthquake occurred, as well as the nearby Cross Lake, Montreal River, and Latchford Faults.The Mohr circle for the Timiskaming location is solidly within the unstable regime for every stress tensor orientation in the map area (Figure S11 in Supporting Information S1).Other faults, though stable for a μ of 0.6, would likewise fall in the unstable regime even for a just slightly weaker fault with μ of 0.55 (Figure S11 in Supporting Information S1).Inclusion of shallow upper mantle/sub-crustal weak zones noticeably shifts faults in the western half of the WQSZ closer to failure by up to 25 MPa, including some that are less optimally oriented.
The LSLR, including the CXSZ, is the only region for which we do not recover any positive Coulomb stress values for faults associated with major historical earthquakes if using μ = 0.6 (Figures 10c,10f,10i,10l,and 10o).
Most faults in the LSLR strike NE-SW, similar to the regional S Hmax orientation and thus are not optimally oriented for thrust reactivation in such a setting.Interestingly, for μ = 0.6, the inclusion of the Farallon slab by itself pushes faults farther from failure.This is because even though the differential stress on the faults increases (i.e., the radius of the Mohr circle expands), the total normal stress on the fault also increases, pushing the Mohr circle to the right along the abscissa, thus increasing its distance from the failure criterion.Nevertheless, similar to the WQSZ, if one were to assume weaker faults with a μ of 0.4 or less, many of the faults, including those associated with major historical earthquakes, would be solidly within the unstable regime.With a μ of even 0.5, most faults would be unstable by 5-10 MPa (Figure S12 in Supporting Information S1).With μ = 0.6, the inclusion of weak zones still brings faults closer to failure by ∼15 MPa in the CXSZ and as much as 20 MPa toward the south near Montreal.For models without the Farallon slab, weak zones increase the Coulomb stress on most faults by ∼10-15 MPa.

Discussion
Using mantle flow models to compute the stress field of eastern North America, we have assessed the influence of Farallon slab buoyancy and intraplate weak zones on the possibility of fault reactivation in intraplate seismic   (2008).Faults with unknown dip are assumed to have a dip of 90°.Purple stars mark earthquake locations; the three northernmost stars mark the locations of the most likely epicenters of the 1811-1812 earthquake sequence (Delano et al., 2021).For the WQSZ, fault locations are from Rimando and Peace (2021) and Lamontagne et al. (2020).Fault strikes approximated using the best fitting great circle through the lineation.Unless otherwise specified, fault dips were assumed to be 56°(A.L. Bent et al., 2003;Rimando & Peace, 2021).Dip of the fault associated with the Timiskaming 1935 earthquake (northwestern-most star) and nearby faults are from A. L. Bent (1996) (strike = 146°, dip = 45°).Other earthquake locations from A. Bent (2022) and A. L. Bent et al. (2003).For the LSLRS, fault strikes approximated for each fault location using the best fitting great circle through the lineation.Fault dips assumed to be 53°, consistent with the values of steeply dipping rift bounding normal faults for the LSLRS as reported by A. L. Bent et al. (2003), A. Bent (1992).Earthquake locations from A. Bent (2022), including the epicenters of the 1663 M 7.0 and 1925 M 6.2 Charlevoix earthquakes and the M 6.3 1732 Montreal earthquake.
zones.Between end member cases of a fully realized negatively buoyant Farallon slab as inferred from seismic tomography (model set A) and that of a neutrally buoyant Farallon slab (model set B), our models reveal a longwavelength stress amplification of up to 100-150 MPa induced by the presence of the slab (Figures 7 and 9), as well as a generally compressive stress state within the continent.To understand the origins of this stress pattern, we examine the change in viscosity and strain rate at the depth of the lithospheric weak zones (Figures 7c-7f).The buoyancy force of the Farallon slab drives localized mantle flow beneath the central eastern US, consistent with previous findings (Forte et al., 2007).This downward flow (Figures 7g and 7h) means an amplified gradient in the vertical velocity in the mantle overlying the downwelling.Because the lithosphere and underlying mantle are viscously coupled, the downward flow in the mantle excites a downward deflection of the overlying lithosphere.This deflection is small due to the high viscosity of the lithosphere, but even these small horizontal gradients in the vertical velocity are enough to produce a broad low amplitude increase in the strain rate across the continent (Figures 7c and 7d) compared to the case with no Farallon slab.While subtle, these higher strain rates, when scaled by the very large viscosity of the lithosphere, result in a strong stress amplification over the central eastern US (Figures 7a and 7b).The peak of this long-wavelength stress perturbation centers over the NMSZ.In Western Quebec and the Lower Saint Lawrence Rift, stress perturbations from the presence of the slab are similarly on the order of 100 MPa.
Superimposed on this broad stress high are peaks in stress due to the inclusion of weak zones.The sharp lateral viscosity gradients imposed by such structures lead to faster flow and higher strain rates within the weak zones (Figure 8), which like the asthenosphere to the lithosphere are viscously coupled to the surrounding rigid lithosphere and crust and impart tractions in these areas.Scaled by the very low viscosity of these weak zones, this leads to lower stresses within the zones themselves, but the tractions in the overlying crust, scaled by the higher viscosity of the crustal lithosphere, yield high stresses localized above the weak zones (Figures 7a and 7b).Stress perturbations arising from weak zones are higher in the models with the Farallon slab by up to 70 MPa relative to cases with no slab because the buoyancy force of this mass anomaly excites more flow in these low viscosity regions (Figure 8).Without the Farallon, weak zone stress perturbations are only on the order of 5-15 MPa, suggesting crustal loading from displacements in a mantle low-viscosity zone are relatively small under far-field tectonic stresses alone, consistent with the findings of Zhan et al. (2016).
However, the depth and extent of the weak zone is also a determining factor in the magnitude of the stress perturbations.Deeper weak zones such as those beneath New Madrid in models A1 and B1 yield stress perturbations at seismogenic depths of only about 10 MPa compared to sub-crustal weak zones such as those in models A1b and B1b, which yield stress perturbations on the order of 75 MPa.This is because the shallower weak zones impart tractions directly to the base of the crust, but the tractions imparted by the deeper low viscosity zones will decay with distance from the weak zone.However, deeper weak zones that are in contact with the sub-lithospheric mantle experience much greater strain rates as the flow excited by the Farallon slab is able to propagate directly into the low viscosity region, concentrating it beneath the seismic zone (Figures 8c and 8d).Thus, with respect to influence on the surface stress field, we identify a trade-off between weak zone depth and coupling to the asthenosphere.Moreover, we find that for weak zones to exert an appreciable influence on intraplate stress magnitude at seismogenic depths, they must be shallow, extending from lower crustal to mid-lithospheric depths and ideally coupled to the crust itself, even when in the presence of full loading from the slab.This is further supported by the pattern of stress orientations predicted from our models.
We also assessed the impact of incorporating gravitational potential energy (GPE) differences (models A0G and B0G) via use of the Crust1.0model (Laske et al., 2013) since lateral variations in density and topography can give rise to variations in deviatoric stress that may help load faults within continental interiors (Ghosh et al., 2009(Ghosh et al., , 2013) ) (see Supporting Information S1).With respect to seismicity, Levandowski et al. (2017) found that for the Great Plains region of North America body forces control the state of stress and style of earthquakes without influence from mantle flow or need for pre-existing crustal weaknesses.Ghosh et al. (2013), on the other hand, find that weak zones at plate boundaries are necessary for accurately reproducing observed stresses and plate motions.Other studies argue for more persistent sources of stress that are concentrated by local structures (Heidbach et al., 2018;Levandowski et al., 2018) since the significant variability and local deviations within some of the seismic zones in the eastern US are generally inconsistent with the continental-scale stress field that would be predicted from far-field tectonic forcing alone (Levandowski et al., 2018;Mazzotti & Townend, 2010).
Incorporating GPE actually reduces the overall stress in the continental interior for both cases with and without the Farallon slab and introduces more high frequency variability in the stress field than we are able to recover from only long-wavelength loading (Figures 9e and 9f and dotted lines in Figures 7a and 7b).Neither model A0G or model B0G have weak zones, but interestingly, they both exhibit a small peak in stress at the location of the NMSZ on the order of 20-30 MPa.This peak is equally pronounced in both models with and without the Farallon slab, and is even greater in amplitude than the stress perturbation caused by the inclusion of a deeper upper mantle weak zone.This suggests that for the NMSZ, crustal density variations could be enough to locally enhance stresses and promote seismicity on pre-existing faults, as has been suggested by other authors (Ghosh et al., 2009;Levandowski et al., 2017).Such a significant density variation could arise from the presence of mafic underplating or high density intrusions originating from the formation of the Reelfoot Rift or the passage of an ancient hotspot (Chu et al., 2013) and is supported in part by the observed seismic structure beneath the NMSZ (Chen et al., 2014(Chen et al., , 2016)).A small peak in stress at around 1,300 km from the center of profile B-B′ in Figure 7a is also more prominent in the models with GPE, reflecting the influence of density and lithospheric viscosity changes at the continent-ocean transition.This is also evident in the viscosity profile at 50 km depth (Figure 7e), where a slight increase in the viscosity gradient at this point leads to slightly higher strains and a small peak in stress magnitude along the passive margin.While small (∼10 MPa), such a stress perturbation in combination with weak zones or structures associated with currently unresolved Mesozoic rift basins along the eastern seaboard could help bring coastal and offshore faults closer to failure.
Similar peaks in stress from the inclusion of GPE are not observed in the WQSZ or CXSZ/LSLR (Figure 7b).However, within the Canadian seismic zones, GPE does induce a positive stress perturbation when not loaded by the Farallon slab, in which case the Northern Appalachian and eastern Canadian region show elevated stress magnitudes.Inclusion of GPE alone without influence from the slab, however, actually increases S Hmax misfits by on average 5°, except in the NMSZ, where it decreases by a couple degrees.On the other hand, including GPE in the presence of Farallon loading reduces misfits by 8°in the Canadian seismic zones, 13°in New Madrid, and as much as 20°in Central Virginia relative to the WSM data (Figure S8 in Supporting Information S1).However, throughout the majority of eastern North America, GPE introduces a counterclockwise rotation of S Hmax in contrast to that observed across the different seismic zones (Mazzotti & Townend, 2010).Overall, while the effect of GPE alone is not insignificant, particularly for the NMSZ, the stress perturbations arising from mantle dynamics and shallow low-viscosity weak zones loaded by the slab are more in line with values necessary to reproduce the observed stress patterns and enhance stress on faults.These results are in agreement with the results of previous studies as well, who find that the mantle flow contribution arising from the history of subduction is the dominant control on the CEUS stress field (Ghosh et al., 2019).
Stress rotations within seismic zones across eastern North America are often considered indicative of long wavelength forcing beyond just tectonic background stress.Mazzotti and Townend (2010) observe 30°-50°c lockwise rotations of the local FM derived S Hmax direction from the regional borehole derived S Hmax direction within several seismic zones in eastern North America.They argue that the consistency of such stress rotations observed in seismic zones separated by large distances, such as the LSLRS and parts of the central-eastern US, suggests a common mechanism.In all cases except the ETSZ, the clockwise rotation in stress induced by the influence of the Farallon slab shifts the S Hmax direction more in line with the observed seismically inferred S Hmax (Figure 6).Most importantly, however, this change is even more pronounced for cases with both the slab and weak zones, and in the case of the CVSZ and the CXSZ, shifts S Hmax from the confidence band of the regional borehole value into that of the seismically derived values.The Charlevoix and Lower Saint Lawrence Seismic Zones have the largest weak zone induced S Hmax rotations of ∼35°and ∼25°, respectively, and the rotation for the CXSZ places S Hmax within a few degrees of the observed value.However, the structural complexities due to the Charlevoix impact crater, as well as variations in S Hmax orientations derived from microseismicity with the CXSZ (Baird et al., 2010), suggest that local stress perturbations may be just as responsible for the anomalous seismicity in this area as is any regional stress perturbation.Nevertheless, the influence of slab loading is shown to induce a clockwise rotation on the intraplate stress field, and the presence of shallow weak zones enhances the magnitude of that stress rotation in line with the observed rotations in some seismic zones.
However, there are places where the predicted S Hmax directions do not match the observed stress field, most notably in New Madrid and Central Virginia, which both exhibit more E-W contraction in the data.One consistent feature that stands out in contrast to earlier works is the lack of spatial heterogeneity of the predicted stresses.Previous studies (e.g., Ghosh et al., 2019; M. L. Zoback & Zoback, 1980) have found stress rotations from NE-SW shortening within the continental interior and cratonic regions to E-W shortening along the coastal margin and in New Madrid, which is amplified if a strong viscosity contrast exists between the craton and the accreted terranes of the eastern margin.Such a rotation yields better fits to some focal mechanisms data (Ghosh et al., 2019;Levandowski et al., 2018) compared to our work, which exhibits a rather uniform and smoothly varying stress field, even in the cases with weak zones.At the regional scale, this smooth stress field is a result of the viscosity structure.Ghosh et al. (2019) find a smoothly varying NE-SW stress field in cases that either lack lateral variations in the continental viscosity or have a smaller degree of variation between the craton and the accreted terranes.The viscosity structure used here, on the other hand, does contain lateral variations as determined from tomography and the Artemieva (2006) thermal model, but they are much more gradual, both laterally and vertically compared to those in the Ghosh et al. (2019) study.For example, in Ghosh et al. (2019), there is a sharp (order of magnitude difference in viscosity) craton boundary traversing the NMSZ and Central Virginia, whereas in our models, the craton boundary manifests as a gradual decrease in viscosity (Figure S8 in Supporting Information S1).As such, the viscosity jump in these regions at any given depth is not an order of magnitude, but closer to a factor of two.The smaller gradients produce more uniform S Hmax orientations across the whole region.However, this smooth viscosity variation may not be heterogeneous enough to adequately recover the spatial heterogeneity of the stress field, suggesting the need for higher lithospheric viscosity gradients, and hence a higher lithospheric activation energy, in future models.
Models with imposed weak zones are the only ones with truly sharp viscosity gradients, and even in these cases, the surface stress field is rather smooth, with differences in angular misfit between models with and without weak zones on the order of ±2°or 3°.As discussed above, weak zone depth is one of the reasons why changes in S Hmax are smaller for some regions than other (e.g., New Madrid), but even for the cases with shallow weak zones, the mean S Hmax orientation in New Madrid is only marginally closer to an E-W direction.However, there are significant rotations of S Hmax , including near E-W orientations in New Madrid, within and on the margins of the weak zones themselves.That is, at the depths of the weak zones rather than at near surface depths (Figure 11).We compute S Hmax misfits at seismogenic depths of around 10 km in the crust above the weak zones because it is those depths that the stress data represents.The S Hmax directions at the depths of the weak zones, on the other hand, are oriented such that they move toward normal to the strike of the boundaries to the viscosity contrasts, meaning the weak zones serve as stress guides, much as the accreted Appalachian terranes do in Ghosh et al. (2019).This would suggest that for the weak zones to rotate S Hmax in such a way as to match the observed stress field, the zone of weakness needs to be shallower in depth, within the upper crust, or to extend throughout the entire lithosphere.If confined to shallower layers, weak zones would also likely need to be of a finer scale requiring higher resolutions than those currently feasible in these models.
There is also an interesting trade-off between the depth of the weak zone and how well connected it is to the underlying asthenosphere, where flow induced by the Farallon slab is much greater.As shown in Figure 8, the New Madrid weak zone in model A1 experiences a much greater excitation of flow because the low viscosity zone is essentially coupled to the stronger asthenospheric downwelling arising from the sinking slab.In contrast, the weak zones under the Canadian seismic zones are a lot shallower and are cut-off from regions of higher mantle flow by stable lithosphere.The Farallon slab does still increase the velocities in these areas, but the effect is more subtle; it nevertheless still has a larger effect on the surface stress field than in New Madrid because those tractions are being applied directly to the base of the crust, much closer to the observed stress field.Thus, near surface stress rotations are lacking in the model results because either (a) the weak zones are too deep, reducing the effect of tractions on the surface stresses, or (b) the weak zones are shallow but cut-off from the deeper flow and excited to a lesser degree.Therefore, the observed S Hmax misfits near the surface provide some constraint on the weak zone depth required for (a) the slab to have an appreciable effect on the weak zones and (b) for the weak zones to have an appreciable effect on the surface stress field, suggesting the need for either shallow weak zones or weak zones that extend through the entire lithosphere.
Stronger influence of weak zones or the slab on the surface stress field, however, will not always yield improved S Hmax fits.In the ETSZ and in the vicinity of the Rome Trough and the Northern Appalachians, the inclusion of weak zones actually reduces the fit to the stress data, rotating it too much relative to the observed local S Hmax .This suggests that despite the prominence of a high conductivity and low seismic velocity zone in this region (Evans et al., 2019;Schmandt & Lin, 2014), corresponding very low viscosities are not essential for reproducing the observed stress field.On the other hand, the ETSZ is one of the few seismic zones for which GPE actually improved the S Hmax fit.The strong shallow density contrast across the NY-AL lineament evident in magnetic data (Thomas & Powell, 2017) suggests lateral buoyancy effects are likely significant in this region.The effects of GPE and of weak zones were tested separately; GPE combined with low-viscosity weak zones may induce even larger or different stress perturbations that more favorably fit the data.For example, in New Madrid, both GPE and the weak zone individually lead to positive rotations in S Hmax , the combined effect of both in the same model would likely enhance that rotation, pushing the S Hmax orientation more in line with observed values.As regards misfits in the CVSZ, there is no discrete zone of weakness underlying that specific area, though the Rome Trough runs nearby.A more localized structure or a broader lower viscosity zone spanning the region as in Ghosh et al. (2019) could indeed help rotate S Hmax into the correct orientation.However, while the NE-SW to E-W stress rotation described above has been claimed to best match observations from FM data along much of the eastern margin (Ghosh et al., 2019;M. L. Zoback & Zoback, 1980), some data challenge this consensus.For example, the recent M 4.8 New Jersey earthquake supports the NE-SW shortening predicted by our models for this location.These results and this earthquake continue to highlight the complexity of the intraplate stress field and the variety of potential mechanisms that can explain it.In general, from these results and the clockwise nature of rotations induced by the slab, we argue that predicted stress orientations that are under-rotated relative to the observed suggest the need for shallower weak zones or greater viscosity contrasts (e.g., NMSZ, CVSZ), and predicted stress orientations that are over-rotated relative to the observed suggest the need for either deeper weak zones or smaller viscosity contrasts relative to surrounding lithosphere (e.g., ETSZ).
The New Madrid, Central Virginia, and the Lower Saint Lawrence River Seismic Zones are those that still require larger rotations to shift S Hmax into the seismically inferred CI band.The observed rotations require stress perturbations on the order of 160-250 MPa, assuming a high friction crust (Mazzotti & Townend, 2010).In our models, the inclusion of zones 10 4 times weaker than the surrounding lithosphere induces stress perturbations on the order of 130 MPa in the presence of slab loading relative to a case with no Farallon loading, yielding total stress magnitudes of 200-250 MPa over the central-eastern US For the NMSZ, only the model with shallow weak zones gets close to the observed S Hmax direction of either the Mazzotti and Townend (2010) or the WSM data, suggesting higher stress perturbations are needed to reach the magnitudes of the observed rotations.This may be achieved by an even shallower or weaker weak zone as discussed or by additional sources of stress, such as GIA or local complexities in the faulting that may help further amplify stresses.For example, Wu and Mazzotti (2007) observe clockwise rotation of GIA stresses in the crust above lithospheric weak zones, and GIA induced S Hmax alone can be nearly E-W in the central-eastern US depending on the radial viscosity structure (Wu, 1997).Strain rates associated with intraplate seismicity are also on the same order of magnitude as, or even slightly smaller than, those of post-glacial rebound (Calais et al., 2006;Mazzotti & Adams, 2005;S. Stein, 2007).In fact, strain rates throughout eastern North America are dominated by GIA (Ghosh et al., 2019;Kreemer et al., 2018), but the orientations of GIA contractions are predominantly N-S, in contrast to the NE-SW shortening direction reflected in the observed stress field (Ghosh et al., 2019), further suggesting that a mechanism other than GIA is responsible for controlling the stress.Nevertheless, regional mechanical models from Mazzotti (2007) suggest that a weak low-viscosity layer in the presence of regional strain consistent with that of GIA can load an upper-crust fault system and amplify the surface velocity gradient, improving the fit to GPS data for the St. Lawrence Valley.Other authors have also found GIA to be a significant driver of seismicity in the St. Lawrence River Valley (James & Bent, 1994;Wu & Hasegawa, 1996) and the NMSZ (Grollimund & Zoback, 2001;Wu & Johnston, 2000).The influence of GIA, in conjunction to that of mantle flow, particularly in the context of lateral viscosity heterogeneity, remains a subject in need of research and is the topic of a companion paper.
Despite the remaining misfits on the predicted stresses in some areas, the rotations induced by the slab and weak zones that we do observe are still important because they reorient the stress tensor into a position more favorable for reactivating faults in some of the seismic zones.In the LSLR and parts of the WQSZ, the majority of faults strike NE-SW in line with the regional S Hmax direction, making them non-optimally oriented for reactivation.When S Hmax rotates into a more E-W orientation, however, these ancient rift bounding faults become more favorably oriented, potentially promoting slip.Even with the rotations we observe (Figure 6), Coulomb stresses on these faults are still relatively low (Figures 10c,10f,10i,10l,and 10o) assuming a μ of 0.6.However, even a small reduction in the coefficient of friction by as little as 0.1 would promote instability on these faults in the modeled stress field (Figure S12 in Supporting Information S1).A μ = 0.5 in the case of both weak zones and the Farallon loading, for example, pushes the majority of faults in the LSLRS into the unstable regime by about 5-10 MPa (Figure S12b in Supporting Information S1).Even without the weak zone, one of the faults nearest the proposed epicenter of the 1732 Montreal earthquake exhibits positive Coulomb stress values.Higher friction values of μ ≥ 0.6 are typical of intraplate crust as a whole (Townend & Zoback, 2000;M. D. Zoback et al., 2002), but locally within pre-existing crustal weak zones, friction may be much lower.In Charlevoix, for example, relationships between seismicity and geological structures suggest the main rift faults of the LSLRS may have abnormally low friction values as small as μ ≈ 0.1 or abnormally high pore fluid pressure (Baird et al., 2010).While such a low value is unlikely, the middle/lower crust in this region is estimated to have a friction coefficient of around 0.5 ± 0.07 based on a joint FM inversion of events within the CXSZ (Verdecchia et al., 2022).Thus, even just a slightly smaller friction coefficient would allow for faults to reach instability under slab loading without the need for such anomalously weak faults.A reduction of only 0.05 in the coefficient of friction would likewise place many faults in the western half of the WQSZ squarely within the unstable regime relative to the Mohr-Coulomb failure criterion (Figure S11 in Supporting Information S1).
Even with μ ≥ 0.6 the likelihood of failure on the Timiskaming fault is clearly resolved in the cases of Farallon loading.Likewise, instability on the Reelfoot Fault associated with the 1811-1812 New Madrid sequence is also apparent in all cases.What these two faults have in common is their optimal orientation for slip within even the background stress field, demonstrating unsurprisingly that the existing fault orientation is a dominant control on its likelihood of reactivation in an intraplate setting.The inclusion of the Farallon slab and/or weak zones, on the other hand, does perturb the magnitudes of stresses in such a way that it brings these faults further into the unstable regime, increasing Coulomb stress on the order of 75-100 MPa (Figure 10).More importantly, assuming weak faults, it can also bring non-optimally oriented faults to failure or closer to failure, as in the case of the WQSZ and the CXSZ/LSLR (Figures S11 and S12 in Supporting Information S1).In this case, the southwest portion of the LSLR is more strongly affected by the Farallon slab than the Charlevoix, reflecting the greater influence of the slab toward the center of eastern North America, beneath which the majority of the Farallon's mass is located (Figure 7).Another peculiarity worth mentioning is that most of the seismicity in the WQSZ is actually concentrated northeast of the old rift, not within the aulacogen itself, yet most of the known faults, including those for which we perform a CFS analysis (Figures 10b,10e,10h,10k,10n), are along the Ottawa-Bonnechere Graben.This would either suggest that the most active faults in this area have simply not been mapped or that the weak zone loading the crust is not bounded by the ancient rift faults but instead located slightly to the northeast, along the suggested position of an ancient transform fault (Sykes, 1978) or the passage of the Great Meteor hotspot (Ma & Eaton, 2007).In either case, if we were to shift the sub-crustal weak zone of Western Quebec slightly to the northeast, we would likely see similar patterns and similar stress magnitude perturbations as those of the current results.
While the substantial influence of the Farallon slab and local weak zones on the intraplate stress field is evident, there are still many local sources of stress perturbations and other factors driving intraplate seismicity to consider.One consideration is the actual viscosity of the weak zones.We only test the presence of weak zones by considering the end member cases of no weak zones and weak zones with a viscosity of 10 4 relative to the surrounding lithosphere.However, weak zone viscosities are likely to vary, possibly by orders of magnitude.There is also likely some trade-off between the impact of the viscosity reduction on the stress perturbations versus the impact of the weak zone depth.Depth is better constrained than viscosity due to the availability of geophysical imaging, but this trade off still needs further investigation.In addition, the more lateral variability that is introduced, the greater the predicted heterogeneity of the stress field should be, potentially improving the fit to the observed S Hmax orientations.This effect can be explored either by varying the activation energy of the lithosphere or by imposing viscosity contrasts at particular locations or structures, as has been explored by Ghosh et al. (2019).Another major consideration is that we do not account for the previous stress history of the fault.In our models, all faults are virtual faults, meaning no physical fault surface is included in the model; rather, the calculated stress tensors at a particular location are applied to a set of faults after the fact to determine their influence on the behavior of a fault of a particular orientation.Were faults present in the model itself, their presence would affect the stress distribution (Steffen, Wu, et al., 2014;Wu et al., 2021).Most importantly, the rupture history of the fault will affect the local change in CFS and the likelihood and location of future earthquakes (R. S. Stein, 1999).This is especially true in light of the long-lived aftershock sequences and temporal clustering that are typical of intraplate settings (DiCaprio et al., 2008;S. Stein & Liu, 2009).The current models cannot explore these processes directly because they are limited in both resolution and the purely viscous rheological formulation of CitcomS.Future models should seek to incorporate true faults, but would require higher lateral resolutions more suited to a regional scale model and elasticity within the crust to properly capture brittle failure.
We also exclusively use free-slip boundary conditions rather than impose plate velocity boundary conditions in order to more naturally recover the motion of plates from mantle flow using a realistic thermal model.Using platelike boundary conditions could lead to a slight continental scale change in the predicted stress direction.However, such a change would likely be uniform across the entire North American plate and would not impact the magnitude and pattern of stress change that results from inclusion of weak zones or the exclusion of the Farallon slab.However, investigating the sensitivity of the stress field to imposed plate boundary velocities versus free slip would be an informative exercise that could also shed light on the degree to which intraplate stresses are controlled by far-field horizontal plate forces versus deeper mantle dynamics.Moreover, future models could also consider the effect of slab strength relative to the surrounding mantle.Slab strength is not likely a significant influence on the present day flow from the Farallon slab because the slab is largely confined to the lower mantle, the bulk of it having penetrated the 660 km discontinuity.However, slab rheology would be relevant for flow models that consider the slab's full subduction history, as stronger or weaker slabs can impact the ability of the structure to pass into the lower mantle, thus impacting the flow and the present day position of slabs (Mao & Zhong, 2021).Slab strength is also globally important for its role in transmitting slab pull stresses and driving plate tectonics in slabs that are still attached to their respective plates (Billen & Hirth, 2007;Saxena et al., 2023;Stadler et al., 2010;Zhong et al., 1998).Changes in their strength could impact global scale mantle flow that affects the global stress field.This impact on intraplate stress has yet to be fully explored.The Farallon slab, however, is no longer attached to any surface plate, so it does not transmit stresses directly, but rather induces stresses through the flow that it excites.It's buoyancy is thus the controlling parameter on its influence on the surface, rather than its rheology.

Conclusions
Using global finite element convection models, we assess the influence of mantle flow induced by the sinking of the Farallon slab on the intraplate stress field of eastern North America.Our work builds upon that of other authors (e.g., Ghosh et al., 2019;Saxena et al., 2021) in several important aspects.While some previous studies (e.g., Ghosh et al., 2019) do use fully global models, many are based on older tomographic models with coarser resolution.We utilize recent global tomography that explicitly constrains the locations of subducting slabs.We also explicitly parameterize the negative buoyancy of the Farallon slab to test the stress field's sensitivity to different buoyancies, rather than solely testing different lithospheric models.Most importantly, we test the effect of local-scale, low viscosity lithospheric weak zones at the locations of geologically mapped inherited tectonic structures, building upon previous work that explores the impact of first order strength contrasts over broad regions of the lithosphere.Finally, we resolve our predicted stress onto known faults in eastern North America to assess the influence of the slab and weak zones on potential fault reactivation.The models help quantify the change in intraplate stress that arises from either the inclusion of weak zones or the exclusion of the Farallon slab.
The Farallon slab introduces a pronounced continent-wide clockwise rotation of the predicted S Hmax direction by as much as ∼20°in places.This rotation is enhanced when lithospheric weak zones are included at the locations of ancient aulacogens and other structures and improves the fit to the observed seismic S Hmax direction within most seismic zones.The presence of weak zones loaded by the Farallon slab at depth can explain in part the pattern of clockwise rotation of the observed FM derived S Hmax direction relative to the regional borehole derived S Hmax as reported by Mazzotti and Townend (2010) in the New Madrid, Central Virginia, Charlevoix, and Lower Saint Lawrence Seismic Zones.However, larger stress perturbations are required to fully reproduce the observed degree of stress rotation in most of the seismic zones, except Charlevoix, suggesting the need for an alternate viscosity structure, weaker weak zones, or another source of stress such as GIA.The magnitudes and S Hmax orientations of the modeled stress perturbations are also dependent on the depth of the weak zones, which have a reduced influence on crustal stresses with greater depth.Misfits on the predicted S Hmax directions are largely due to the smooth viscosity gradient in the lithosphere, but they inform us about the depth of the low viscosity zones needed in order for the slab to have an appreciable effect on the weak zones and for the weak zones to exert an appreciable control on intraplate stress near the surface.Shallow/sub-crustal weak zones that are in contact with the crust are most effective at augmenting stress at seismogenic depths, but deeper weak zones are excited to a greater degree by the negative buoyancy of the Farallon slab.
Even without weak zones, influence of mantle flow from the Farallon slab is enough to explain fault instability at the locations of some significant historical earthquakes, such as the Reelfoot Fault in New Madrid-the site of one of the 1811-1812 M 7-8 earthquakes-and the Timiskaming Fault in Western Quebec, on which the 1935 M 6.2 earthquake occurred.Across all seismic zones, the inclusion of weak zones brings the majority of faults closer to failure.However, fault orientation and fault weakness remain key controls on their reactivation potential.Weaker faults (μ ≤ 0.4) in the NMSZ would not be enough to reactivate them unless optimally oriented, but weak faults (μ ≤ 0.5) could explain fault reactivation in the WQSZ and CXSZ in the modeled stress field.Many previous studies have argued that lithospheric mantle heterogeneity alone is the primary control on intraplate seismicity (Levandowski et al., 2017;Saxena et al., 2021;Zhan et al., 2016).Such heterogeneity in mantle viscosity controls the spatial variability of the velocity gradients by diverting flow to lower viscosity areas.We find that mantle heterogeneity remains a key factor and that stronger viscosity gradients are likely needed to better reproduce the observed state of stress.Nevertheless, the Farallon slab has the effect of augmenting the impact of those gradients, which ultimately moves many intraplate faults closer to failure.While lithospheric mantle heterogeneity may govern the spatial distribution and shortening direction of intraplate stress, epeirogenic processes are strongly influential on the magnitudes of stress required to reactivate some intraplate faults and to explain continent-wide stress rotations.Ultimately, it is challenging to fully recover the complexity of the intraplate stress field with any one model, but the results presented here and the misfits on the predicted stresses inform us about the lithospheric structure and weak zone characteristics necessary to explain the observed stress field in the presence of mantle loading.

Data Availability Statement
The seismic velocity model (TX2019, Lu et al., 2019) used to constrain the temperature field is freely available from IRIS at http://ds.iris.edu/ds/products/emc-tx2019slab/.The BurnMan mineral physics software (Cottaar et al., 2014;Myhill et al., 2021) used in the conversion of seismic velocities to temperature is a Python package freely available from https://github.com/geodynamics/burnman/. The scripts and files used in the velocity to temperature conversion are included within the V2T_Conversion directory in the Supporting Information data available at https://data.caltech.edu/records/wh6q5-80b36and are described in Supporting Information S1 document.The seafloor age grid used to constrain the oceanic lithospheric temperatures is from Seton et al. (2020) and is freely available from https://earthbyte.org/webdav/ftp/earthbyte/agegrid/2020/.The thermal model used to constrain the continental lithospheric temperature field is from Artemieva (2006) and is freely available at http:// www.lithosphere.info/downloads.html.Stress data is available from the World Stress Map Project at https://www.world-stress-map.org/download(Heidbach et al., 2018).Data on rift geometries (polygons) and geology used to construct the weak zone input is available from the Central Eastern United States Seismic Source Characterization for Nuclear Facilities database (http://www.ceus-ssc.com/index.htm)(CEUS-SSC-Project, 2011).CitcomS version 3.3.1 (McNamara & Zhong, 2004;Moresi et al., 2014;Tan et al., 2006;Zhong et al., 2000) is an open-source publicly available geodynamics code maintained by Computational Infrastructure for Geodynamics and available from https://geodynamics.org/resources/citcoms or https://doi.org/10.5281/zenodo.7271920.The code has excellent strong scaling on parallel computers with up to 1,000 s of cores.Our calculations were performed on the NSF ACCESS HPC clusters Stampede2 at the Texas Advanced Computing Center (TACC) and Anvil at Purdue University using 16 nodes and 768 processors per run.Typical compute times per simulation are about 8-12 hr and are highly dependent on the viscosity structure, which can span up to six orders of magnitude.The Supporting Information S1 also includes a copy of the version of CitcomS used, as well as all pre-and postprocessing scripts.
Data sets containing model inputs, outputs, and procedures, as well as scripts for making the figures presented in the paper, are available from the CaltechDATA repository (https://data.caltech.edu/records/wh6q5-80b36).These include the final 3D temperature field input to CitcomS; the 3D viscosity field; the modeled 3D velocity field; the modeled S Hmax orientation, principal stresses, stress magnitudes, and strain rates; and fault segment data taken from sources cited within the text.See Supporting Information S1 for an explanation of each of the included datasets.Figures were constructed with either Python or GMT 6 (or PyGMT) (Wessel et al., 2019).

•
The Farallon slab causes continentwide clockwise rotation of S Hmax ; shallow intraplate weak zones enhance stress rotations and magnitudes • Stress perturbations from both the Farallon slab and intraplate weak zones heighten Coulomb failure stress on some faults associated with historical earthquakes Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Regional setting.(a) Historical seismicity in eastern North America and paleo-tectonic structures.Large orange circles are significant intraplate earthquakes of M ≥ 5.5, and black dots are background seismicity of M ≥ 3 from the USGS-ANSS Comprehensive Earthquake Catalog.Dark green polygons are geologically mapped Proterozoic aulacogens from Whitmeyer and Karlstrom (2007) and the CEUS-SSC-Project (2011).Light green polygons are Mesozoic rift basins from Withjack et al. (1998).Dark blue polygon is the Eastern Tennessee Shear Zone interpreted from Thomas and Powell (2017) and Powell et al. (2014).Individual structures as referred to in the text are labeled in the figure.(b) Stress field (S Hmax orientations) in eastern North America from the World Stress Map (Heidbach et al., 2018), colored by faulting style.Black dots are background seismicity as in (a).Black ellipses outline major named seismic zones.

Figure 2 .
Figure 2. Vertical slices of 3D seismic velocity structure and corresponding mantle thermal structure.(Left) Vertical slices of the TX2019 (Lu et al., 2019) velocity field for both dV s and dV p along two transects.Profile A-A′ trends SW-NE across eastern North America and passes through the New Madrid Seismic Zone (NMSZ) (marked by the orange inverted triangle) and the Western Quebec and Charlevoix Seismic Zones (red inverted triangle).Profile B-B′ trends W-E across the US and passes through the NMSZ.(Right) Vertical slices along the same two transects of the effective temperature field (δT ) calculated from the velocity anomalies.Locations of the seismic zones are marked as in the velocity plots; the Farallon slab stands out as an elongate low temperature thermal anomaly.

Figure 3 .
Figure3.Radial temperature and viscosity profiles for models in CitcomS.(a) Radial temperature profile for the thermal model used in CitcomS; includes both the seismically derived mantle temperatures and lithospheric temperature from the plate cooling model andArtemieva (2006).Light gray shaded region covers the full range of temperature found at a particular depth.The deep extent of very cold temperatures is due to the presence of cold slabs that penetrate deep into the mantle.The teal green band marks the 2σ spread of the temperatures at a given depth, containing 98% of the values.The black line is the geotherm as used in the velocity to temperature conversion, fromSteinberger and Calderwood (2006).(b) Radial viscosity profiles for different models.Dark teal line marks the average radial viscosity structure used in CitcomS, as computed from the temperature model on the left.The teal green shading is likewise the 2σ range of viscosity, and the gray shading covers the full range of viscosity values found at a given depth.Dotted black line is the VM5a viscosity profile ofPeltier et al. (2015); dashed-dotted black line is the seismically derived viscosity profile used in CitcomSVE in the work ofPaulson et al. (2005); the dashed black line is the viscosity profile fromSteinberger and Calderwood (2006) that best reproduces observed plate velocities; and the solid black line is the radial viscosity structure from the best-fitting model inMao and Zhong (2021).The right-hand y-axis shows the depth-dependent pre-factors applied to the temperature-dependent viscosity for the four-layers distinguished by the dark blue lines.

Figure 4 .
Figure 4. Changes in misfit between observed and modeled S Hmax orientations at 10 km depth.(a) S Hmax orientations of model A0 (with Farallon, no weak zones) colored by misfit (|modeled observed|) to the observed stress orientations from the World Stress Map of Heidbach et al. (2018) (light gray lines), where the maximum misfit is 90°.(b) Difference in the S Hmax misfit between model A1 and model A0 (i.e., models with the Farallon slab, but with vs. without weak zones).(c) Difference in S Hmax misfit between model A0 and B0 (i.e., models without weak zones, but with vs. without the Farallon slab).(d) Difference in S Hmax misfit between model A1 and model B1 (that is, between the case with both Farallon and weak zones and the case with weak zones but no Farallon).

Figure 5 .
Figure 5. Rotation in S Hmax orientation between cases with versus without weak zones and/or Farallon loading at 10 km depth.(a) S Hmax orientations of model A0 colored by misfit as in Figure 4. (b) S Hmax rotation between models A1 and A0 (i.e., effect of including weak zones in models with the Farallon slab).Positive values indicate clockwise rotation; negative values indicate counterclockwise rotation.(c) S Hmax rotation between models A0 and B0 (i.e., effect of including the Farallon alone, without weak zones).(d) S Hmax rotation between model A1 and model B0 (i.e., combined effect of including both Farallon and weak zones).(e) S Hmax rotation between model A0 and A0G (i.e., effect of including gravitational potential energy (GPE) in models with Farallon but no weak zones).(f) S Hmax rotation between model B0 and B0G (i.e., effect of including GPE in models with neither Farallon nor weak zones).

Figure 7 .
Figure 7. Stress, strain, and viscosity profiles.(a) Second invariant of deviatoric stress along track B-B′ as in Figure 2, traversing W-E across the U.S. and through the New Madrid Seismic Zone (NMSZ) and Eastern Tennessee Seismic Zone.Light blue lines: models in set A (i.e., with Farallon slab); dark blue lines: models in set B (i.e., neutralized Farallon slab).Solid lines: no weak zones (models A0 and B0).Thick dashed lines: weak zones at variable depths depending on seismic zone (models A1 and B1).Dashed-dotted lines: weak zones all at 25-75 km depth (models A1b and B1b).Dotted lines: with gravitational potential energy but no weak zones (models A0G and B0G).(b) Stress magnitude as in (a) but along track A-A′ as shown in Figure 2, traversing SW-NE across the U.S. and Canada and through the NMSZ, Western Quebec Seismic Zone, and CXSZ/Lower Saint Lawrence Rift System.(c, d) second invariant of deviatoric strain rate at 50 km depth along tracks B-B′ and A-A′, respectively.Lines as in (a, b).(e, f) Viscosity profiles at 50 km depth along B-B′ and A-A′, respectively.Line-styles as in (a, b).(g, h) Vertical slices of the viscosity field along tracks B-B′ and A-A′, respectively, for models with the Farallon slab; black arrows show velocity field.Black boxes denote the region shown in the crosssections in Figure 8 for each profile, respectively.Purple lines: contours of the viscosity field to help highlight the position of the slab.Locations of the seismic zones along the profile are marked as in (a, b).(i, j) Viscosity slices and velocity field as in (g, h) but for models with a neutralized Farallon slab.

Figure 8 .
Figure 8. Zoomed in cross-sections of the viscosity and computed flow field for the section of the arcs outlined by the black boxes in Figures 7g and 7h: zoom ins of profile B-B′ shown in left column; zoom ins of profile A-A′ shown in right column.(a, b) Cross sections for model A0 (Farallon, no weak zones).(c, d) Cross sections for model A1 (Farallon, with weak zones).(e, f) Cross sections for model B0 (no Farallon, no weak zones).(g, h) Cross sections of model B1 (no Farallon, with weak zones).Locations of the seismic zones along the profile are marked by the triangles, as in Figure 7.

Figure 9 .
Figure 9. Deviatoric stress magnitude for eastern North America for different model cases.(a) Second Invariant of deviatoric stress for model A1 with both Farallon and weak zones placed at variable depths.(b) Stress perturbation from including weak zones in models with the Farallon slab (i.e., between models A1 and A0).(c) Stress perturbation from the Farallon slab alone, without weak zones (i.e., between models A0 and B0).(d) Combined stress perturbation from including both the Farallon and weak zones(i.e., between models A1 and B0).(e) Stress perturbation induced by gravitational potential energy (GPE) for models with the Farallon, but no weak zones (i.e., between models A0 and A0G).(f) Stress perturbation induced by GPE for models with neither the Farallon nor weak zones (i.e., between models B0 and B0G).

Figure 10 .
Figure 10.Coulomb failure stress (map inset) resolved onto mapped faults in the New Madrid Seismic Zone (NMSZ) (a, d, g, j, m), Western Quebec Seismic Zone (WQSZ) (b, e, h, k, n), and the Lower Saint Lawrence Rift System (LSLRS) (c, f, i, l, o) calculated using μ = 0.6.(a-c) Results from model A0 (with Farallon, no weak zones).(d-f) Results from model A1 (with Farallon, with weak zones).(g-i) Results from model A1b (with Farallon, with weak zones all at 25-75 km).(j-l) Results from model B0 (no Farallon, no weak zones).(m-o) Results from model B1 (no Farallon, with weak zones).For the NMSZ, fault locations from Thompson et al. (2020).Strike and dip information for the Axial Fault (strike = 46°, dip = 90°), the New Madrid North Fault (strike = 29°, dip = 72°), the Risco Fault (strike = 92°, dip = 82°), and the North Reelfoot Fault (strike = 167°, dip = 30°) are fromCsontos and Van Arsdale (2008).Faults with unknown dip are assumed to have a dip of 90°.Purple stars mark earthquake locations; the three northernmost stars mark the locations of the most likely epicenters of the 1811-1812 earthquake sequence(Delano et al., 2021).For the WQSZ, fault locations are from Rimando and Peace (2021) andLamontagne et al. (2020).Fault strikes approximated using the best fitting great circle through the lineation.Unless otherwise specified, fault dips were assumed to be 56°(A.L.Bent et al., 2003;Rimando & Peace, 2021).Dip of the fault associated with the Timiskaming 1935 earthquake (northwestern-most star) and nearby faults are from A. L.Bent (1996) (strike = 146°, dip = 45°).Other earthquake locations from A.Bent (2022) and A. L.Bent et al. (2003).For the LSLRS, fault strikes approximated for each fault location using the best fitting great circle through the lineation.Fault dips assumed to be 53°, consistent with the values of steeply dipping rift bounding normal faults for the LSLRS as reported by A. L.Bent et al. (2003), A. Bent (1992).Earthquake locations from A.Bent (2022), including the epicenters of the 1663 M 7.0 and 1925 M 6.2 Charlevoix earthquakes and the M 6.3 1732 Montreal earthquake.

Figure 11 .
Figure 11.Viscosity structure and predicted S Hmax direction from Model A1b at (a) 25 km depth, (b) 30 km depth, (c) 50 km depth, and (d) ∼100 km depth, showing the rotations in S Hmax that arise within and at the margins of the weak zones.Stress data was only output for near-surface depths down to 50 km and is not plotted in panel (d).

Table 1
Dimensional Constants and Model Parameters a Values separated by a semi-colon indicate values for upper and lower mantle, respectively.

Table 2
Model Parameters Tested in Different Model Cases Models prepended with A have the full Farallon slab as resolved from seismic tomography.b Models prepended with B have the Farallon slab entirely neutralized by the buoyancy ratio (i.e., Farallon BR 100%).c Gravitational potential energy (GPE) indicates whether the effect of GPE as derived from the Crust1.0model was included.
a Model no.b a d WZ η prefactor is the viscosity reduction applied to the intraplate weak zones.e WZ depth is either variable in accord with the geophysically inferred depths of different structures as discussed in the text (e.g., New Madrid Seismic Zone: 70-150 km, Eastern Tennessee Seismic Zone: 15-50 km, RT: 80-200 km, LSLRS: 25-75 km) or is constant between 25 and 75 km depth.
Mazzotti and Townend (2010)s in different seismic zones.Gray: polar histogram of the local S Hmax orientations from the World Stress Map, binned every 10°.Orange colors: regional borehole (BH) derived mean S Hmax and 90% CI fromMazzotti and Townend (2010).Pink colors: local focal mechanism (FM) derived mean S Hmax and 90% CI from Mazzotti and Townend (2010) (see legend).Values of the mean BH and FM S Hmax are labeled on the theta axis for each seismic zone.blue-green lines: local mean S Hmax from models with the Farallon slab (A0, A1, A1b, A0G).Dark blue lines: local mean S Hmax from models with a neutralized slab (B0, B1, B1b, B0G).Line-styles differ for different weak zone cases (see legend).NMSZ: New Madrid Seismic Zone, ETSZ: Eastern Tennessee Seismic Zone, CVSZ: Central Virginia Seismic Zone, WQSZ: Western Quebec Seismic Zones, CXSZ: Charlevoix Seismic Zone, LSLSZ: Lower Saint Lawrence Seismic Zone.