Reactivation Potential of Intraplate Faults in the Western Quebec Seismic Zone, Eastern Canada

The intraplate western Quebec seismic zone (WQSZ) in eastern Canada experiences moderate seismicity that mainly results from reactivation of inherited structures under the present‐day, NE‐SW‐striking regional maximum horizontal stress (SH) and, possibly to a minor extent, through stress perturbations in response to glacio‐isostatic adjustment. This work comprises the first numerical stress simulation‐based study that predicts the preferred spatial distribution, trends, and sense of slip of contemporary fault reactivation, which may have implications for possible fault segmentation patterns in the WQSZ. We show that mostly NNW‐SSE to NW‐SE‐striking faults exhibit the highest slip tendency values. Spatial patterns of slip tendency and kinematics of reactivation are consistent with the observed seismicity. In an area where Quaternary‐active faults have yet to be systematically identified, we have narrowed down areas to focus on for more detailed, future neotectonic investigations that could provide sound foundation for seismic hazard assessments. This study demonstrates the applicability of slip tendency analysis to identifying potentially active faults in stable continental regions worldwide.

Canada is the fact that most evidence of initiation of fault activity under the current stress regime typically postdates glaciation (Adams, 1989). There has been a debate over whether recent faulting occurs as a result of tectonic stress or glacio-isostatic adjustment, or both (e.g., Adams, 1989;Brooks & Adams, 2020;Wallach et al., 1995). Despite the surficial nature and exposure of some pop-up structures and offset boreholes within quarries in Ontario and Quebec, a tectonic origin is a preferred interpretation for most of these features due to the compatibility of their orientations and kinematics with the current regional stress field in eastern North America (Wallach & Chagnon, 1990;Wallach et al., 1993Wallach et al., , 1995. Besides, detailed studies of the state of stress in eastern Canada (Mazzotti & Townend, 2010) indicate that magnitudes of long-wavelength stress perturbations such as postglacial rebound stresses are an order of magnitude lower, hence minor in comparison to mid-crustal stresses. It is possible, however, for postglacial rebound stress to cause high enough stress perturbations if these are concentrated on faults with an unusually low coefficient of friction (μ ∼ 0.1) RIMANDO AND PEACE 10.1029/2021EA001825 2 of 14 Faults with names are labeled as follows: 1-Timiskaming, 2-Cross Lake, Havelock. The topography is derived from a 30-m-resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer global digital elevation models (https://asterweb.jpl.nasa.gov/gdem. asp) and the fault traces are from the Geological Survey of Canada's WQSZ faults map (Lamontagne et al., 2020). (b) The orientations of the 210 WQSZ faults are summarized in a rose diagram which has 18 bins (10° intervals). (c) S Hmax orientations in the WQSZ (Mazzotti & Townend, 2010); 045° is the average regional stress orientation, while 028° and 073° are the extreme values. (Mazzotti & Townend, 2010). An alternative theory for the enhanced seismicity in the WQSZ, which forms an NW-SE-striking linear cluster of intermediate focal depth (8-18 km) earthquakes across eastern Canada, is the passage of North America over the Great Meteor Hotspot in the Mesozoic. This theory postulates that thermo-mechanical weakening of ancient faults and shear zones was caused by the Great Meteor Hotspot along the interpreted NW-SE-trending hotspot track (e.g., Dineva   (a) Earthquake epicentral plots are color-coded according to depth and scaled to magnitude. (b) Earthquake focal mechanisms with labels of year and magnitude of notable events. Well-localized seismicity is from Adams et al. (1988Adams et al. ( , 1989, Bent (1996), Bent et al. (2002Bent et al. ( , 2003, Du et al. (2003), Horner et al. (1978), Ma and Eaton (2007), Seeber et al. (2002), Wahlström (1987), and the earthquake bulletins of both the Natural Resources Canada (https://earthquakescanada.nrcan.gc.ca/index-en.php) and the United States Geological Survey (https://www.usgs.gov/natural-hazards/earthquake-hazards/earthquakes). references therein), and that present-day seismicity are unusually long-lived aftershock sequences of large past earthquakes (e.g., Stein & Liu, 2009). Shallow seismicity (<8 km), however, appears to be more randomly distributed throughout the region and cannot be accounted for by this theory. In addition, it is not clear why seismicity would be sustained to the present-day through this thermal weakening mechanism.
Knowledge of the nature, distribution, and extent of seismogenic structures in Eastern Canada, including their potential for causing large magnitude earthquakes, is of paramount socio-economic importance (e.g., Morell et al., 2020). A major earthquake in eastern Canada could trigger a chain of events in the insurance industry and have far-reaching economic consequences (Le Pan, 2016). However, there is currently no map of the active structures in the WQSZ, a region that hosts major population centers, including Ottawa and Montréal. As in most intraplate settings, there is a rarity of well-preserved recent fault/fold scarps in the WQSZ (McCalpin, 2009), probably as a result of long recurrence intervals on slow-moving faults (Stein, 2007), and due to glacial peneplanation (Dyke et al., 2002). Similarly, the short temporal coverage of instrumental seismicity records is not enough to create an inventory of detailed earthquake source models and to determine earthquake recurrence intervals in intraplate regions. However, we know from global examples that large earthquakes are far more common than previously thought in "stable continental regions" (Calais et al., 2016;Rimando et al., 2021). In eastern North America, instrumental, historical, and paleoseismic records show that the region has been shaken by devastating earthquakes, such as the 1929 M S 7.2 Grand Banks and the 1811-1812 M W 7.2-8.2 New Madrid earthquakes (Bent, 1995;Hasegawa & Kanamori, 1987;Tuttle et al., 2002). There is a wealth of studies on the occurrence of paleoearthquakes in the WQSZ, inferred from subaqueous mass transport deposits and liquefaction features. However, these have yet to be linked with specific earthquake-generating faults (Aylsworth & Hunter, 2003;Brooks, 2013Brooks, , 2014Brooks, , 2015Brooks & Adams, 2020).
Ottawa is the seat of Canada's national government and Montréal has the second-largest economy (based on GDP) in Canada (Statistics Canada, 2020a). Both are major national and global commercial, financial, technological, and cultural hubs. As of 2020, the total population in the Ottawa-Gatineau and Montréal metropolitan areas is over 5.5 million (Statistics Canada, 2020b). Despite the economic and political importance of these areas, infrastructure in these cities is mostly old and is likely to be ill-prepared for possible future strong ground shaking during an earthquake (Goda, 2019;Kovacs, 2010;Mirza & Ali, 2017). Mapping of active faults is therefore critical to estimating the seismic hazard in an area and, consequently, to formulating a reliable building code.
In the absence of a map of active faults, a good first step is to determine the potential of pre-existing faults to be reactivated under the current stress field. A good correlation between relatively high slip tendency and evidence of Quaternary activity has been shown in different tectonic settings with wide-ranging levels of seismic activity (e.g., Worum et al., 2004;Yukutake et al., 2015). For instance, in Japan, most faults that exhibited relatively high slip tendency values (slip tendencies of ≥ 0.6) were consistent with faults that were classified independently as active (Yukutake et al., 2015). While previous studies have described the likely orientation and kinematics of fault reactivation conceptually under the current stress field (Daneshfar & Benn, 2002;Rimando, 1994;Rimando & Benn, 2005), no work has characterized in detail the slip tendency and expected slip kinematics of specific faults in the WQSZ. This work presents the first study which uses 3D numerical stress simulations to explore the preferred spatial distribution and trends, and predicted sense of slip of reactivated pre-existing structures in the WQSZ under the current tectonic stress field. We show how fault reactivation potential studies can be useful for identifying key areas or fault populations to focus on for more detailed seismic hazard assessment.

Data and Methods
Slip tendency (Morris et al., 1996) is the ratio of the shear stress to the normal stress on a fault surface, which is expressed as the following equation: where T S is the slip tendency, τ is the shear stress, and σ n is the normal stress.
RIMANDO AND PEACE

10.1029/2021EA001825
Slip tendency analysis has been used in different tectonic settings worldwide to characterize the relative likelihood of populations of faults to slip under current or past stress fields (e.g., Morris et al., 1996;Peace et al., 2018;Worum et al., 2004). Faults that are likely to slip are those with a higher ratio of shear stress to normal stress. Faults subject to a certain stress field are characterized as "optimally oriented" if the set of strikes and dips yield high slip tendency values.
It should be noted, however, that we make implicit assumptions in using this technique, which introduce some limitations on how closely our model represents reality. Following the Wallace (1951) and Bott (1959) hypothesis wherein slip on faults is expected to occur along the direction of the maximum resolved shear stress, this method assumes simple planar faults and a relatively uniform stress field, and it neglects fault interaction, fault block rotation, and internal deformation. Despite these caveats, it has been shown through numerical studies to be a good first-order approximation (Dupin et al., 1993;Pollard et al., 1993), as the deviation between actual and theoretical slip kinematics are on average less than 10°. A strong match between modeled slip tendencies and directions and natural reactivated fault plane orientations and slip kinematics from global geological and seismological data sets lends support to the use of slip tendency analysis as a useful prediction tool (Collettini & Trippetta, 2007;Lisle & Srivastava, 2004).
We modeled the reactivation potential of faults in the WQSZ using the "Slip Tendency" function in the "Stress Analysis" module of the software Move™ by Petroleum Experts Limited (https://www.petex.com/), and determined the predicted sense of fault slip using the "Slicken 1.0" software (Xu et al., 2017). In both software packages, the stress tensor and fault plane orientations were the required inputs.
We primarily used the average regional stress orientation of 045° but also ran simulations using the extreme values (028° and 073°) determined for the Montréal and Gatineau zones by Mazzotti and Townend (2010). These were based on Bayesian inversion of earthquake focal mechanisms and calculation of the weighted averages (weights based on quality) of borehole breakout stress measurements within 250 km of the area of interest. The orientations of S H measured from both sources were consistent and roughly parallel. Focal mechanism inversions revealed a nearly vertical σ 3 and nearly horizontal σ 1 and σ 2 , defining a reverse faulting stress regime, which is consistent with findings of previous crustal stress orientation studies (Heidbach et al., 2010;Reiter et al., 2014;Snee & Zoback, 2020). Knowing that σ 1 is horizontal enabled us to use the S H from borehole measurements as an approximation for σ 1. We estimated regional stress magnitudes by assuming a critically stressed crust (Townend & Zoback, 2000;M. D. Zoback & Zoback, 2002), in which differential stress values are at a level such that optimally oriented Andersonian faults are on the verge of slipping. In a reverse-faulting stress regime, the principal stresses can be computed using the following equations: where σ 1 is the maximum horizontal stress, σ 2 is the minimum horizontal stress, σ 3 is the vertical stress, λ is the pore-fluid factor (pore fluid pressure, P f divided by the σ 3 ), ρ is the average crustal density (2,700 g/m 3 ), g is the gravitational acceleration (9.8 m/s 2 ), z is the depth (in m), R is the principal stress difference ratio (typically 0.6 in this intraplate continental region), F is the frictional parameter, and μ is the coefficient of friction.
We calculated the stress magnitudes at a mid-seismogenic zone depth of 10 km. Our choice of coefficient of friction (μ) and pore fluid factor (λ) values considers the fact that the maximum possible differential stress RIMANDO AND PEACE 10.1029/2021EA001825 (σ 1 −σ 3 ) values at mid-crustal depths in this region are unlikely to exceed 200 megapascals (MPa), which is based on previous modeling studies and extrapolation of in situ field measurements at shallower depths in Canada and in intraplate settings in general (e.g., Hasegawa et al., 1985;Lamontagne & Ranalli, 1996, and references therein). This condition of relatively low differential stress necessitates the μ and λ values to be lower and higher, respectively, than the values that are commonly assumed (e.g., Byerlee, 1978;Townend & Zoback, 2000). For the purposes of this study, we used a σ 1 −σ 3 value close to the 200 MPa upper limit. We adapted a μ of 0.5 and a λ of 0.6, which are intermediate to the values previously determined from modeling of the conditions for slip of earthquakes in southeastern Canada (M. L. Zoback, 1992). A low coefficient of friction at shallow depths likely results from the presence of thick phyllosilicate-rich fault gouges (e.g., den Hartog et al., 2020) and in the mid-crust through the presence of a dense network of fractures and faults (Ito & Zoback, 2000). The upward migration of mantle-derived mixed H 2 O-CO 2 fluids is a proposed source of high pore fluid pressure that enables the reactivation of high-angle faults in eastern Canada (Sibson, 1989). Hence, we applied the following stress tensors in our analyses: σ 1 = 028°-073° = 277.09 MPa, σ 2 = 118°-163° = 174.34 MPa, and σ 3 = vertical = 105.84 MPa (Figure 1c; Mazzotti & Townend, 2010).
With the exception of a few faults (e.g., Doughty et al., 2012), information on the deeper structure of most faults is lacking in the WQSZ. Where available, information on the geometry of faults from subsurface imaging techniques is confined to depths of less than 100 m at best (e.g., Doughty et al., 2012). However, most papers, reports, and geological cross-sections on the area indicate steep fault dips with an average of around 60° (e.g., Lovell & Caine, 1970;Rimando, 1994;Rocher & Tremblay, 2001). We therefore assumed a dip of 60° to build 3D fault surface models. The 3D models were built by projecting surfaces from shapefiles of the Geological Survey of Canada's latest WQSZ faults map (Lamontagne et al., 2020). We also ran simulations over a range of fault dips (at 15° increments) to test the effect of underestimating or overestimating the actual fault dip on the calculated slip potential and sense of slip of faults.

Results
Our slip tendency analysis reveals that under the present-day regional stress field (Figure 1c), and assuming a dip of 60° on faults in the WQSZ, mostly NNW-SSE to NW-SE striking faults tend to have relatively higher slip tendencies (compared to the more NE-SW-striking faults in the eastern sector), and are therefore, considered optimally oriented to be reactivated (Figure 3a). Faults are predicted to slip mostly either as pure reverse faults or as oblique reverse faults (Figure 3b). The magnitude and spatial distribution of the slip tendency varies, albeit insignificantly, over the range of the assumed possible orientations of σ 1 (i.e., 028°, 045°, and 073°). For instance, the average slip tendency values of NW-SE-striking and NE-SW-striking faults (for 60° dip) vary on average by ∼5% and ∼10% of each other, respectively ( Figure 3a and Table S1). On the other hand, the preferred kinematics of slip has a more noticeable spatial variation. Faults that are oriented perpendicularly to σ 1 are expected to be reactivated with a pure reverse slip. Consequently, the distribution of pure reverse faults appears to rotate clockwise as the σ 1 azimuth is rotated clockwise, and coincidentally, decreases in abundance due to the decreasing number of faults oriented at a high angle to the σ 1 in the WQSZ (Figure 3b). While the magnitudes of slip tendency appear to vary as a function of the assumed fault dip, the NW-SE striking faults of the WQSZ mainly exhibits generally relatively higher slip tendencies, with average values being 20% higher than NE-SW-striking faults for σ 1 = 045° (Figure 4; Table S1). The predicted slip in the 15°-60° dip scenarios all exhibit either pure reverse or oblique reverse faulting. It is notable that in the 75° dip scenario, there tends to be significantly more faults with a dominant strike-slip component, and in the 90° dip scenario, faults are expected to be reactivated almost entirely as pure strike-slip faults ( Figure 5). Additionally, for faults dipping 75° and 90° with σ 1 = 028° and 045°, NE-SW-striking faults tend to have slightly higher average slip tendency values than NW-SE-striking faults, possibly due to these faults being optimally oriented to be reactivated with a dominant strike-slip component (Figures S1-S4; Table S1).
It is worth noting that, with the exception of the 30° dip scenario (which seems to be the near-optimal dip for fault reactivation in the WQSZ), slip tendency values decrease as the assumed fault dip increases (Figure 4). Most of these observations on the variations of slip tendency and expected slip kinematics as a function of dip and as a function of σ 1 orientation, apply as well to other σ 1 orientation scenarios (σ 1 = 028° and σ 1 = 073°) and to the different dip scenarios (15°, 30°, 45°, 75°, and 90°), respectively ( Figures S1-S4).

Discussion and Conclusions
Our results show that the NNW-SSE-striking to NW-SE-striking faults of the WQSZ consistently exhibit relatively higher slip tendency values. However, although we considered the effect of adapting different dip values, we modeled dips uniformly across the WQSZ during each simulation. In reality, however, it is unlikely the case that the dips of all the faults in the WQSZ are uniform. Nonetheless, unless the faults in the west sector are all dipping vertically or nearly vertically, and assuming that the Wallace (1951) and Bott (1959) hypothesis applies and pore fluid pressure conditions and fault frictional properties vary insignificantly throughout the WQSZ, then the slip tendency values of the NW-SE-striking faults, mostly to the west of the WQSZ, are expected to be higher for most conceivable variable-dip scenarios.
There is also a good correspondence between the modeled kinematics of fault reactivation (i.e., predominantly reverse and oblique reverse) (Figures 3b and 5) under the assumed range of stress tensors and the actual kinematics of recent natural seismicity as indicated by earthquake focal mechanisms (Figure 2b). Additionally, in upstate New York and southeastern Ontario, previous work has demonstrated that recent natural seismicity is spatially associated dominantly with the NW-SE-striking population of faults (Daneshfar & Benn, 2002;Rimando, 1994). Therefore, there seems to be multiple lines of evidence that give credence to the results of our modeling.
In reality, however, the assumptions of Wallace (1951) and Bott (1959) are rarely entirely met (Lisle, 2013), which could feasibly cause deviations between our modeling results and the actual fault activity and kinematics. The results of this study, nonetheless, should provide a first-order approximation of the distribution of fault activity in the WQSZ. It can be argued, however, that the orientation of stress is fairly homogenous in the region, and the effects of fault block rotations and fault interactions are unlikely to significantly affect our results given the scale of our study area and the resolution of our modeling. After all, it is possible that the objections to the assumptions of the Wallace-Bott hypothesis may apply only at a relatively local scale (Lisle, 2013).
Future numerical stress simulations would benefit from more detailed 3D fault geometry models. For instance, some faults in the region have been shown to exhibit flat-ramp-flat and listric fault geometries (Busch et al., 1996;Rimando & Benn, 2005), which could possibly result in more complex along-strike and down-dip patterns of slip tendency and slip kinematics. Future studies should consider improving 3D fault models by integrating data from geophysical subsurface imaging, natural and trenching exposures, and seismicity. Having such data will not only reduce the uncertainty in our modeling, but will also help address concerns regarding the validity of the approach used in this study.
While fault slip in the WQSZ could very well be mainly tectonic in nature, glacial unloading could play an important role as a trigger (Wu & Hasegawa, 1996). In reverse/thrust faulting stress regimes such as in the WQSZ (M. L. Zoback, 1992), glacial unloading could possibly decrease the magnitude of the vertical σ 3 , resulting in a higher differential stress (e.g., Muir-Wood, 1989). The extent to which glacio-isostatic adjustment can influence the regional stress field and, consequently, the modeled distribution of slip tendency of faults in the WQSZ should also be further explored in future work. Ongoing isostatic rebound in eastern Canada (George et al., 2012) can be modeled by progressively increasing the differential stress.
Our findings on the range of conditions in which faults in the WQSZ may experience reactivation have implications for the assessment of seismic hazards in the area. For instance, possible fault segmentation can be inferred from our slip tendency and slip kinematics maps (e.g., Figure 3). In the absence of documented actual earthquake rupture segments, which is the case in the WQSZ, structural, geologic, and geometric characteristics such as bends, stepovers, gaps, branching, intersections, and changes in the orientation and sense of slip of faults are typically used as preliminary criteria for assessing fault segmentation (e.g., DePolo, 1989;DePolo et al., 1991;Knuepfer, 1989;Rimando & Knuepfer, 2006). Assuming that the surface geometries of faults resemble the geometries at depth, the length of inferred segments can be utilized to compute possible earthquake magnitudes from empirical relationships between moment magnitude (M W ) and fault length by Wells and Coppersmith (1994). Such criteria for rupture endpoints, however, are not as clear in normal faults as they are in strike-slip and reverse faults (Knuepfer, 1989). In the WQSZ, where structures are mostly rift structures reactivated as reverse faults, interpreting segments may not be as straightforward using these criteria since some faults may still exhibit segmentation associated with earlier periods of normal faulting.
Alternatively, examination of patterns on slip tendency and predicted slip maps, which were modeled using the present-day stress field, could allow us to identify fault segments. Endpoints of fault segments could be inferred where there are abrupt changes in slip tendency values or fault kinematics. These fault lengths can then be used to provide initial estimates of potential earthquake magnitudes from scaling relationships by Wells and Coppersmith (1994). The modeled slip type (i.e., dominantly strike-slip, reverse, or normal) could also provide further constraints for potential earthquake magnitude estimates, as calculations differ according to the type of slip (Wells & Coppersmith, 1994). Different ranges of possible earthquake magnitudes are associated with similar-length faults that have different slip-types. A 70-km-long strike-slip fault, for instance, is capable of generating M W 7-7.4 earthquakes, while a normal fault of the same length is capable of generating a wider possible range of earthquake magnitudes of M W 6.6-7.8 (Wells & Coppersmith, 1994). Additionally, the areal extent over which long-period ground motions occur has also been demonstrated to vary as a function of the style of faulting (Aagaard et al., 2004), and our results place constraints on the style of faulting.
While a detailed analysis of fault segmentation and estimation of the associated possible earthquake magnitudes associated with individual interpreted segments is beyond the scope of this study, we look at the example of the NNW-SSE-striking major faults that comprise the Timiskaming graben to illustrate the potential applicability of slip tendency analysis for this purpose. The Timiskaming, Cross Lake, Montréal River, and Latchford faults (faults 1-4 in Figure 1), including an unnamed fault which lies very close to the epicentral location of the 1935 M W 6.1 Timiskaming earthquake (Figure 2), exhibit relatively high slip tendency values and slip type distributions (for σ 1 = 045°, fault dip 60°) that are fairly similar such that each fault can be considered as single segments. Assuming that the shortest mapped major fault in the OBG, the Montréal River fault (fault 3 in Figure 1), ruptures along its entire ∼35-km-length, the possible associated earthquake magnitude would be at least M W 6.6, which is much higher than the largest recorded earthquake (1935 M W 6.1) in the vicinity. Similarly, if the other faults rupture even one-third our one-fourth of their entire length, then higher magnitudes would still be expected. More detailed mapping of late Quaternary fault traces and documentation of evidence of coseismic late Quaternary fault offsets along these optimally oriented faults, however, would be necessary to confirm the presence of such long fault segments in fact exist in the OBG.
The NW-SE-striking to WNW-ESE-striking faults starting at westernmost reach of the Ottawa-Bonnechere graben (OBG) near Crystal Falls Fault (fault 5 in Figure 1), down to the southeast at the intersection of an unnamed NW-SE-striking OBG fault with the Rideau Lakes Fault (fault 26 in Figure 1), also exhibit relatively high slip tendency values. This area is particularly promising for future neotectonic investigations as the area exhibits fault-bounded topographic relief. Detailed characterization of the tectonic geomorphology of the area could possibly reveal observable recent morphotectonic features.
The WNW-ESE-striking faults in the central OBG, from the Hazeldean Fault in Ottawa (fault 24 in Figure 1) in the west to the Ile Bizard Fault in Montréal (fault 42 in Figure 1) in the east, are also optimally oriented and, consequently, of interest for follow-up field studies. Exploring these may be more challenging, though, due to the location of faults in the center of the Ottawa River basin. Fault geometry may have to be characterized from inversions of potential field geophysical data sets (e.g., Yan et al., 2019), which have provided constraints at a variety of complex geological settings, while offset may have to be established through a combination of drilling and analysis of deep (several kilometers) seismic reflection profiles.
Finally, while this work brings us one step forward in our efforts to assess the activity of faults in the WQSZ, caution should be taken with mistaking slip tendency as the risk for fault rupture (Yukutake et al., 2015). The probability of an earthquake occurring on certain faults is a topic that is best tackled by more appropriate neotectonic analyses. Indeed, knowledge of which faults are more likely to be reactivated in the WQSZ, and in other stable continental regions worldwide where slip tendency analysis can be applied, can help earthquake geologists narrow down which faults/areas to prioritize and focus on for more detailed active fault mapping and paleoseismic studies that can reveal earthquake magnitude and recurrence interval estimates.