The Canadian Beaufort Sea continental margin of northwestern Canada is a Cenozoic convergent margin, potentially representing a rare case of incipient subduction. Here, we produce P- and S-wave seismic velocity models of the crust and the uppermost mantle using recordings from regional earthquakes. Our models reveal a northwest-dipping very low-velocity anomaly within the crust (δV up to -15%) beneath the Romanzof Uplift. We interpret this low-velocity feature to correspond to a weaker and thicker crust due to shortening and stacking of igneous and sedimentary rocks. The co-location of the thickened crust and lack of present-day seismicity indicates that north-south compression is accommodated by slow, aseismic deformation in the narrow margin beneath the Romanzof Uplift or more broadly offshore. Neither interpretation requires a subduction initiation process.

The Leech River fault (LRF) zone located on southern Vancouver Island is a major regional seismic source. We investigate potential interactions between earthquake ruptures on the LRF and the neighboring Southern Whidbey Island fault (SWIF), which can be interpreted as a step-over fault system. Using a linear slip-weakening frictional law, we perform 3D finite element simulations to study rupture jumping scenarios from the LRF (source fault) to the SWIF (receiver fault), focusing on the influences of the offset distance, fault initial stress level, and fault burial depth. We find a smaller offset distance, a higher initial stress level on either fault or a shallower fault burial depth will promote rupture jumping. Jumping scenarios can be interpreted as the response of the receiver fault to stress perturbations radiated from the source fault rupture. We demonstrate that the final rupture jumping scenario depends on various parameters, which can be collectively quantified by two keystone variables, the time-averaged Over Stressed Zone (where shear stress exceeds static frictional strength on the receiver fault) size $\overline{R_e}$ and the receiver fault initial stress level. Specifically, a smaller offset distance, a higher initial shear stress level, or a shallower burial depth will lead to a larger $\overline{R_e}$. The seismic moment on the receiver fault increases with increasing $\overline{R_e}$. When $\overline{R_e}$ reaches the threshold dependent on the receiver fault initial stress level, the rupture becomes break-away.

Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent, physics-based framework to connect, interpret, and predict diverse geophysical observations across spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare and verify numerical codes based on boundary-element, finite-element, and finite-difference methods, in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space, where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs for sufficiently large model domains and small grid spacings. However, discrepancies in rupture fronts of the initial event are influenced by the free surface and various computational factors. The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to numerical resolution and domain-size-dependent loading. Despite such variability, key properties of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer a community-based example to improve numerical simulations and reveal sensitivities of model observables, which are important for advancing SEAS models to better understand earthquake system dynamics.

Source processes of injection induced earthquakes involve complex fluid-rock interaction often elusive to regional seismic monitoring. Here we combine observations from a local seismograph array in the Montney basin, northeast British Columbia, and stress modeling to examine the spatial and temporal evolution of the 30 November 2018 M 4.5 hydraulic fracturing induced earthquake sequence. The mainshock occurred at ~ 4.5 km in the crystalline basement two days following injection at ~ 2.5 km, suggesting direct triggering by rapid fluid pressure increase via a high-permeability conduit. Most of the aftershocks are located in the top 2 km sedimentary layers, with focal mechanisms indicating discrete slip along sub-vertical surfaces in a ~ 1 km wide deformation zone. Aftershock distribution is also consistent with static stress triggering from the M 4.5 coseismic slip. Our analysis suggests complex hydraulic and stress transfer between fracture/fault networks needs to be considered in induced seismic hazard assessment.

Crustal velocity variation within impact-related seismic zones is commonly attributed to mechanisms such as pore pressure changes, dense fracture network, and compositional variation. In this study, we combine seismic tomography, rock physics analysis, and potential field modeling to quantitatively investigate the mechanisms that influence crustal velocity variation in the Charlevoix Seismic Zone (CSZ), a meteorite impact-related seismic zone in eastern Canada. Earthquakes in the CSZ align along two broad NE-SW trending clusters related to reactivated paleo-rift faults. Within the impact structure, the earthquakes are diffusely distributed and lower velocity bodies are ubiquitous which can be attributed to crustal damage from tectonic inheritance exacerbated by the meteorite impact. The Bouguer gravity anomaly decreases southeastward across the St. Lawrence River due to density disparity between rocks in the Grenville Province and the Appalachians. We find a higher velocity body northeast of the impact structure that does not exhibit an observable gravity anomaly, which suggests the presence of a rock (e.g. anorthosite) of comparable density but a higher elastic moduli within another rock (e.g. charnockite). Outside the impact structure, compositional variations control velocity changes, whereas inside the impact structure, velocity variations can be explained by porosity enhancement of up to 10% by low (0.1) aspect ratio cracks. Our results suggest that intense fracturing and compositional alteration, rather than pore pressure, control velocity variations, hence earthquake processes in the CSZ.

Intraplate tectonic stress fields are complex due to the imprint of a long geological history. Here we use a new dataset of earthquake focal mechanism solutions and relocated events to investigate the relationship between regional stress, crustal strength, and seismicity in the Charlevoix Seismic Zone (CSZ), the most active seismic zone in eastern Canada. Our stress inversion shows that SHmax gradually rotates clockwise from approximately St. Lawrence River-parallel near the surface to river-perpendicular in the lower crust, as postglacial rebound stress becomes increasingly dominant at greater depth. The stress rotation occurs primarily between ~13 and ~26 km depth, where glacial rebound induced stress perturbation is further amplified by a “weaker” middle crust of an estimated apparent friction coefficient of ~0.5. Finally, depth-dependent b-values confirm the rheological difference between upper and middle crust in the CSZ.

Numerical simulation of rupture dynamics provides critical insights for understanding earthquake physics, while the complex geometry of natural faults makes numerical method development challenging. The discontinuous Galerkin (DG) method is suitable for handling complex fault geometries. In the DG method, the fault boundary conditions can be conveniently imposed through the upwind flux by solving a Riemann problem based on a velocity-strain elastodynamic equation. However, the universal adoption of upwind flux can cause spatial oscillations in cases where elements on adjacent sides of the fault surface are not nearly symmetric. Here we propose a nodal DG method with an upwind/central mixed-flux scheme to solve the spatial oscillation problem, and thus to reduce the dependence on mesh quality. We verify the new method by comparing our results with those from other methods on a series of published benchmark problems with complex fault geometries, heterogeneous materials, off-fault plasticity, roughness, thermal pressurization, and various versions of fault friction laws. Finally, we demonstrate that our method can be applied to simulate the dynamic rupture process of the 2008 Mw 7.9 Wenchuan earthquake along/across multiple fault segments. Our method can achieve high scalability in parallel computing under different orders of accuracy, showing high potential for adaptation to earthquake rupture simulation on natural tectonic faults.