Bridging time scales of faulting: From coseismic to postseismic slip of the Mw 6.0 2014 South Napa, California earthquake.

Transient fault slip spans time scales from tens of seconds of earthquake rupture to years of aseismic afterslip. So far, seismic and geodetic recordings of these two phenomena have primarily been studied separately and mostly with a focus on kinematic aspects, which limits our physical understanding of the interplay between seismic and aseismic slip. Here, we use a Bayesian dynamic source inversion method, based on laboratory-derived friction laws, to constrain fault stress and friction properties by joint quantitative modeling of coseismic and postseismic observations. Analysis of the well-recorded 2014 South Napa, California earthquake shows how the stressing and frictional conditions on the fault govern the spatial separation between shallow coseismic and postseismic slip, the progression of afterslip driving deep off-fault aftershocks, and the oblique ribbon-like rupture shape. Such inferences of stress and frictional rheology can advance our understanding of earthquake physics and pave the way for self-consistent cross-scale seismic hazard assessment.

INTRODUCTION

Seismic and aseismic slip are two primary modes of fault behavior whose spatial distribution controls the earthquake potential of a fault and may inform on its mechanical properties. Seismic and aseismic slip tend to occur on separated areas of a fault (), as manifested in the large-scale division of faults into creeping and locked segments and in the modest spatial overlap between coseismic slip and afterslip. Earthquakes occur in the locked portion of faults and originate in the seismogenic zone, surrounded by predominantly aseismic slip at the top and bottom. Several physical mechanisms might determine the seismic or creep behavior of a fault. For example, aseismic behavior close to the surface has been attributed to the presence of fault gouge with low confining stresses (), and the seismogenic depth is bounded by a temperature-controlled transition to plastic sliding (). In addition, changes in lithology () or pore pressure (, ) can influence the preferred type of slip.

A commonly observed form of transient aseismic slip is the afterslip that follows earthquakes in areas adjoining their seismic rupture (, ). There is a large variability in the amount of afterslip following different earthquakes (), which indicates a complicated relationship between physical conditions on the fault and coseismic and postseismic slip. Postseismic slip can occur close to the surface (), in areas that show a coseismic slip deficit (). A well-known example of a fault that generates ample seismic and postseismic slip is the Parkfield segment of the San Andreas fault, which produced a series of Mw (moment magnitude) 6 earthquakes in the 19th and 20th centuries (). The most recent event occurred in 2004, releasing twice as much moment postseismically than coseismically, and was a subject of physics-based studies of transient slip ().

Physics-based modeling, including dynamic source inversion, is one approach to advance our fundamental understanding of the partitioning between seismic and aseismic slip. Rate-and-state friction laws (–), based on laboratory experiments at relatively low slip rates, offer a framework that allows explaining both seismic and aseismic phenomena in dynamic models. A fault can be partitioned into seismic and aseismic portions by its spatially heterogeneous frictional properties (, , ). In particular, steady-state friction can be velocity weakening (potentially unstable, seismic) or velocity strengthening (dominantly aseismic). The framework can be extended to higher slip rates, relevant to coseismic slip, by introducing a fast-velocity-weakening mechanism. High-speed friction experiments (, ) show that a range of fault materials weakens substantially at slip rates above ~0.1 m/s (), which has been attributed to thermally activated processes such as flash heating (). Incorporating fast-velocity-weakening friction into the rate-and-state earthquake model (–) leads to qualitative changes in its behavior as dynamic strength is close to zero, and the difference between prestress and static strength increases (). Moreover, fault areas that are velocity strengthening at low slip rate may switch to velocity-weakening behavior at fast slip rates, as was suggested for the Tohoku earthquake () and observed in laboratory experiments under conditions with increasing normal stress (). A primary goal of the dynamic source inversion proposed here is to infer the friction properties of a fault based on observations of coseismic and postseismic slip.

The 24 August 2014 Mw 6.0 South Napa, California earthquake has a particularly well-documented abundance of coseismic and postseismic slip, making it a good target for dynamic source inversion. It ruptured one of the recently mapped branches of the West Napa fault (). The earthquake’s right-lateral strike-slip mechanism is consistent with the orientation of this fault. The shallow part of the fault (<3 km depth) span two lithological units (): The northern half is positioned on the contact between Cretaceous rocks (sandstone, melange, etc.) from the Franciscan Complex and Cenozoic sediments, while the southern half of the fault goes below the Napa River and is embedded in a layer of Quaternary sediments (–) whose thickness increases in the southward direction from 1.5 km to more than 2 km (see Fig. 1A for the fault position with respect to the regional geology).

Fig. 1.

Maps and fits of coseismic and postseismic data.

(A) Position of the fault with respect to the local geologic conditions (white, Quaternary sediments; green, Cretaceous rocks; purple, Cenozoic volcanic rocks), based on (, ) and seismic (black circles) and GPS (blue squares) stations. (B) Comparison between observed seismograms (black) and our best-fitting model seismograms (red). Kernel density estimates (KDEs) of the posteriors are displayed in blue. Station names and maximum displacements are indicated on the left and right, respectively. (C) Fit between observed coseismic GPS displacements (black arrows) and synthetic data (red arrows); KDEs are displayed in blue. Positions of the stations with their names are shown on the map (black circles) with respect to the fault (white rectangle), with total slip color-coded in white to red. Star denotes the epicenter. (D) Comparison between observed postseismic surface displacement (black) and our best-fitting model GPS (red). KDEs are displayed in blue, while errors of real data are shown as error bars. Station names are indicated on the left, and maximum displacements in centimeters are indicated on the right. (E) Comparison between observed postseismic GPS displacements (black) and our best-fitting model synthetics (red). KDEs are displayed in blue, while errors of real data are shown as error bars. Station names and maximum displacements are indicated on the left and right, respectively.

The 2014 South Napa earthquake is well studied, including measurements of surface slip and afterslip over the whole length of the rupture (, ). Kinematic studies of coseismic and postseismic slip (, –) agree on the main source characteristics. The rupture nucleated at 10 km depth and propagated up-dip and northward for 8 to 10 s along a 13-km distance, generating strong seismic waves amplified toward the north due to the source directivity and a sedimentary basin (). The event produced a 12-km-long surface rupture and rapid shallow afterslip (, ). It was also followed by approximately 1000 aftershocks that occurred mostly below the coseismic rupture (). While most of the shallow coseismic slip was concentrated in the northern half of the rupture, an unusually large shallow afterslip occurred on the southern half, starting 3 hours after the earthquake and continuing over the next several months (). This spatial difference in the release of shallow slip has been attributed to spatial variability of either the local geology (, ) or the coseismic stress changes ().

This paper provides a unifying dynamic model of the 2014 South Napa coseismic rupture and subsequent afterslip based on rate-and-state friction with fast-velocity-weakening effect () included. To this aim, we extend the Bayesian dynamic inversion method () to integrate seismic and geodetic data on diverse time scales from seconds to months. The inferred rupture properties reconcile and refine previous independent studies. In addition, the dynamic inversion results enable previously unexplored physical interpretations of the connection between lithology and friction properties, and insights into the role of enhanced weakening in rupture propagation and on the mechanisms of coexistence of seismic and aseismic slip on a fault. We further examine the effects of heterogeneous dynamic parameters on the rupture propagation and arrest, showing that heightened prestress drove coseismic rupture at depth, while velocity strengthening was the main arresting mechanism near the surface. We also show how spatially heterogeneous frictional rheology affects the development of both coseismic and postseismic slip in the shallow zone, affecting the time scale over which slip is released. In addition, we suggest that spatially limited afterslip had a role in triggering off-fault aftershocks, which were mainly observed below the coseismic rupture.

RESULTS

Our Bayesian dynamic inversion aims at constraining the friction parameters and initial fault stresses that govern the space-time evolution of both seismic and postseismic slip and produce ground motions consistent with seismic and geodetic data. We assume a vertical planar fault of 20 km × 15 km size that reaches the surface and has a strike of 165° (Fig. 1A), which is a simplified representation of the geometry constrained by the position of the surface rupture and relocated aftershocks (). The set of dynamic model parameters determined by the inversion procedure is shear prestress T0, direct effect parameter a, state evolution parameter b, reference friction f0, and characteristic slip distance L as two-dimensional (2D) fields, and weakening velocity and initial velocity as 2D fields on the smaller (velocity-strengthening) portion of the fault. The friction law, simulation techniques, data errors, model parameterization, and sampling of the posterior probability density function (pPDF) are described in Materials and Methods. The result of the inversion is an ensemble of models with spatially varying dynamic rupture parameters, statistically representing samples of the pPDF.

We model data from 10 near-source strong-motion accelerometers, seven continuous GPS stations, and four alignment arrays capturing surface fault offsets (, ) at vineyards crossing the fault (see Fig. 1A for their position with respect to the fault). In addition, forward modeling of a larger dataset of coseismic GPS displacements is used for verification of the inversion results. We use a 1D layered crustal velocity profile based on the GIL-7 model () with added low-velocity surface layers. We consider the frequency range of 0.05 to 0.5 Hz for the seismograms and daily sampled GPS displacements (the original dataset from UNAVCO; URL provided in Acknowledgments). Alignment array measurements were irregular in time, so we use all accessible data points from four sites where substantial afterslip was detected (initial measurements were 2 to 5 days after the earthquake, two more in the first 10 days, and two more between 10 and 30 days). Data from both GPS and alignment arrays are considered in the first 30 days after the earthquake.

Figure 1 (B and C) compares the coseismic data with our best-fitting model, which has a variance reduction of 0.49 for seismograms. The figure also displays the statistical variability of the simulated data due to the model uncertainty using kernel density estimates (KDEs) of the posteriors, representing histograms smoothed by a Gaussian function (). The fit is generally good; we attribute the major portion of the data misfit to unmodeled 3D Earth structure in the velocity model and nonplanar nonvertical geometry of the real fault. Postseismic displacements at GPS stations and alignment arrays are displayed in Fig. 1, D and E, respectively. We note that the displacements recorded by the GPS stations are of the order of 1 cm only due to the rather large distance of the stations from the fault and the moderate size of the earthquake. Nevertheless, the fit is still good despite postseismic displacement amplitudes being much lower than amplitudes of seismograms or alignment array displacements. The fit of the coseismic GPS displacements used for verification is comparable with the fit of those used for inversion (see fig S5). The surface slip measurements provide major constraints on afterslip. They are fitted very well due to their relatively high implicit weights in the inversion and lack of direct trade-offs with the other data. The total variance reduction of all GPS data is 0.63.

Kinematic properties and stress drop

Coseismic ruptures in our model ensemble nucleate at a mean depth of 10.46 ± 0.30 km and propagate upward and to the north. They create two major patches of coseismic slip at 3 and 6 km depths (Fig. 2A), coinciding with the maximum stress drop areas, which locally reach 50 MPa (Fig. 2D). The coseismic rupture propagates for about 8 s at an average speed of ~2.4 km/s, releasing a seismic moment of (1.97 ± 0.10) × 1018 Nm. More than 90% (1.9 × 1018 Nm) of the moment is released within the first 5 s. Rise time fluctuates between 0.5 and 1 s and increases above 1 s in the shallowest 2 km. The rupture reaches the surface, over a length of more than 5 km. The final ruptured area attains a ribbon-like shape of width ~5 km and length ~12 km, and its major axis shows an unusual oblique orientation. Areas of shear stress increase (Fig. 2D) concentrate around the rupture edges.

Fig. 2.

Kinematic rupture parameters and their statistics.

Ensemble averages of (A) coseismic slip, (B) afterslip, and (C) total slip on the fault. Blue lines in the coseismic slip and afterslip map indicate the rupture front and the tip of the shallow afterslip in 1-day increments after the coseismic rupture, respectively. Ensemble averages of (D) coseismic, (E) postseismic, and (F) total stress drop. Contours (threshold of 0.3 m) of slip (red) and afterslip (black) with thinner lines denoting SD are displayed in all six panels. Gray dots represent aftershocks (Northern California Earthquake Data Center) with a fault-perpendicular distance of <5 km.

Postseismic slip evolves continuously after the coseismic slip around most of the rupture area. In particular, shallow afterslip starts within 20 to 24 hours from the southern side of the coseismic rupture (~8 km along strike; see Fig. 2B) and expands rapidly in the first 3 days at ~1.5 km/day toward the south. Expansion continues over the whole modeled period of 30 days, albeit with decreasing rate. This produces a substantial (~14-MPa) postseismic stress drop comparable with coseismic stress drop at the same depths. We also observe ~10 cm of shallow postseismic slip even at the northern (coseismically ruptured) portion of the fault (Fig. 2, A and B), in agreement with the shallow slip measurements ().

Smaller patches of notable afterslip (with a maximum of ~0.4 m) are located at about 7.5 km depth, partially overlapping with coseismic rupture. Some of our models show additional patches of afterslip further away from the earthquake, which we consider unconstrained due to their highly variable occurrence among models and minimal impact on synthetic data. Overall, the postseismic slip increases the total seismic moment of the earthquake by ~40%, with a ~15% increase happening during the first day after the earthquake (see fig. S4). Deep postseismic slip mostly happens in the first week after the earthquake, while shallow slip unfolds over a longer period of time (see fig. S4).

Frictional properties

The rupture properties described above stem from the dynamic rupture models, whose parameters are constrained by the inversion. A parameter of particular interest is (b-a), which quantifies the relative importance between direct and evolution effects of rate-and-state friction, and controls the stability of slip: Positive values are associated with velocity-weakening frictional behavior and unstable slip, while negative values imply velocity strengthening and stable slip. Another important dynamic parameter is initial shear stress T0. We show ensemble averages of spatial distributions of (b-a) and T0 along the fault in Fig. 3 (A and B) and their uncertainties in Fig. 3 (C and D), respectively. Figures S1 and S2 show all the other inverted parameters. We discuss only dynamic parameters in the slip area and closely adjoining regions of the fault, where we can consider them well constrained by data. The along-fault width of this zone of interest expands from 5 to 6 km at depth to 15 km near the surface due to the presence of substantial shallow afterslip. To facilitate discussion about the depth dependence of friction, we also show depth profiles of selected parameters in Fig. 3E, calculated as horizontal averages over the slip region.

Fig. 3.

Selected dynamic parameters and their statistical properties.

(A) Ensemble average of (b-a). Gray dots denote the aftershocks as in Fig. 2. Red and black lines indicate contours of slip and afterslip, respectively. (B) SD of (b-a). (C) Same as (A) but for prestress T0. (D) Same as (B) but for relative SD of T0. (E) Horizontal averages of (b-a), T0, characteristic slip L, a, reference friction f0, and weakening velocity on the ruptured part of the fault. Black dots denote averages of individual ensemble models, while the red line with error bars show ensemble mean and SD, respectively. Vertical black line denotes (b-a) = 0. For the remaining parameters, see the Supplementary Materials.

Stresses in the shallow zone, above 5 km depth, decrease with decreasing depth. This is the case for both the normal stress set a priori (see Materials and Methods for details) and the shear stress constrained by the inversion. The shallow zone hosts a combination of frictional parameters that limit rupture propagation and stabilize the fault, reducing both rupture speed and peak slip rates: velocity-strengthening rheology, increasing characteristic slip distance L up to ~1.5 m, weakening velocity up to 3 m/s, and values of reference friction f0 above 1. We note that the large values of f0 found at shallow depth are unusual for rocks but can be attributed to cohesion (see Materials and Methods). The horizontal transition zone between coseismic and postseismic rupture areas (7 to 12 km along strike) is characterized by low prestress and more velocity-neutral friction (b-a close to zero), and overlaps with both a change in lithology and a geometrical feature where the mapped fault starts to bend more toward the west. The strengthening rheology of the afterslip area is more pronounced in the south with (b-a) ~ −0.01 as opposed to −0.005 in the northern part. The high relative SD of T0 in the shallow postseismic zone is a manifestation of a strong trade-off between T0 and f0 (Fig. 4A)—we note that the two apparent clusters in Fig. 4A are caused by imperfect posterior sampling.

Fig. 4.

Plots documenting various modeling features for discussion.

(A) Scatter plot of local dependence between T0 and f0 at a position located 11.5 km along strike and at 1.5 km depth. (B) A priori estimate of dynamic stress drop calculated as prestress T0 minus steady-state friction with fSS at = 0.1 m/s as friction coefficient. (C) Along-strike distribution of coseismic slip (red), afterslip (black), and total slip (blue) at 200 m depth. Error bars denote the ensemble mean and SD. (D) Ensemble mean and SD of (b-a) at 200 m depth. Circles denote the along-strike position of three points, for which the inset shows the afterslip development. (E) Development of stress rate (error bars showing ensemble mean and SD) and the number of aftershocks per day (black points) in the deep postseismically slipping area denoted by the green rectangle in (B).

At greater depths (Fig. 3E), the main coseismic rupture area is dominated by velocity-weakening friction (b-a > 0), low values of the rate-and-state characteristic slip distance L ~0.25 m, and weakening velocity 0.1 m/s, while reference friction f0 is still relatively high (~0.75). The (b-a) parameter has higher uncertainty here than in other (strengthening) parts of the fault, most likely due to the dominant fast-velocity-weakening effect at high slip rates. On the other hand, the relative SD of prestress (~0.01) is minimal in the coseismic zone, as it is well constrained by seismic waves originating from this area.

Substantial heterogeneity in dynamic parameters exists around the 7.5 km depth overlapping with the patch of notable deep afterslip. Friction becomes velocity strengthening due to the increase in a, while f0 decreases to 0.55. Other dynamic parameters (L, weakening velocity) have values similar to those in the coseismic region. The fracture and radiated energies are (9.2 ± 0.8) and (4.5 ± 0.7) MJ/m2, respectively. The radiation efficiency of the earthquake is thus 0.33 ± 0.11.

DISCUSSION

We have conducted a Bayesian dynamic inversion of the 2014 South Napa earthquake, creating a set of ~7500 models that help explain both coseismic and postseismic data in a unified framework of the rate-and-state fast-velocity-weakening friction law. The model describes frictional behavior over a wide range of time scales, from coseismic seconds to postseismic weeks. The simulations are enabled by a combination of fully dynamic and quasi-dynamic modeling of the coseismic and postseismic phases, respectively. The resulting main source features are consistent with those identified by previous analyses of the coseismic and postseismic data. In particular, the inferred coseismic upward and northward rupture propagation with two main patches of slip and the position of substantial shallow afterslip are consistent with published measurements (, , ) and kinematic models (, , , , ).

The joint modeling of earthquake slip and afterslip allows us to constrain dynamic parameters on larger portions of the fault than only coseismic dynamic inversion would. This is enabled by the fact that inferred coseismic and postseismic slip are spatially complementary, although some afterslip takes place in the coseismic area, especially near its border. The central part of the coseismic zone is dominated by velocity-weakening (b-a > 0) friction. Still, the rupture also propagates through velocity-strengthening (b-a < 0) areas near the free surface and above the hypocenter at about 7.5 km depth. The shallow zone is of particular interest because it hosts a transition from seismic to aseismic slip, which occurs over a short distance of 1 km, in agreement with the surface measurements. In addition, the shallow afterslip rate is spatially heterogeneous, being faster near the coseismic zone than further away. These complexities are encoded in the dynamic parameters, in particular (b-a). The deeper strengthening zone not only ruptured coseismically but also hosted notable afterslip, triggering aftershocks off the fault and below the coseismic rupture. Below we discuss and interpret those important features in detail.

Coseismic rupture arrest

We find evidence for different mechanisms driving rupture arrest at deep and shallow depths. At seismogenic depths, in areas between 5 and 10 km depth that are well within the rupture, slip rates exceed the weakening velocity, and thus, friction drops close to the fully weakened friction coefficient fw. This is not the case close to the rupture edges as we demonstrate in Fig. 4B, which shows an estimate of the dynamic stress drop assuming slip rate lower than the weakening velocity. The large negative stress drop values at the edges suggest that the arrest is primarily driven by low prestress T0 with respect to the residual strength. As the rupture approaches the low-prestress barrier, it slows down, and its peak slip rate diminishes [as expected from theoretical arguments ()], which eventually prevents the fast-velocity-weakening effect. Closer to the surface, the strength excess decreases, and the velocity-strengthening effect gains importance as the rupture arrest mechanism by keeping the peak slip rates below the fast-velocity-weakening limit. This is especially the case in the shallow southern portion of the fault.

The velocity-strengthening zone at 7.5 km depth is an exception to this picture, as the difference between initial stress and reference friction is much lower there (see also the small stress drop estimate in Fig. 4B). This feature not only slows down the coseismic rupture but also produces a patch of large afterslip (Fig. 2E). Low prestress is our preferred rupture arrest mechanism at large depth because the alternative, velocity-strengthening friction, would induce larger deep afterslip that would be inconsistent with the GPS data.

Interplay between coseismic and postseismic ruptures at shallow depths

The unique feature of our modeling is to adopt a single friction law for both the coseismic and postseismic ruptures, in contrast to their independent treatment in previous works [e.g., (–)]. In the case of the South Napa earthquake, the shallow zone above 3 km depth hosts an abrupt horizontal change from seismic to aseismic rupture. The northern portion of the shallow fault ruptured coseismically, switching within ~1 km to the south to primarily postseismic rupture (Fig. 3B). The total shallow slip (coseismic and postseismic) has two local maxima, one in the coseismic zone at around 6 km along strike and one in the postseismic zone at 11 km along strike (Fig. 4C). The local minimum (~9 km along strike) coincides with the border between the coseismic and postseismic slip areas and is associated also with nearly zero total stress drop (Fig. 2F). These characteristics are well constrained by data from the alignment arrays and are in good agreement with previous models of shallow slip [e.g., ()].

The distribution of frictional properties in our results (Figs. 3B and 4B) shows that the whole shallow part of the fault is velocity strengthening, including the coseismic portion. This feature of rate-and-state dynamic models is implied by physical mechanisms (low normal stresses, temperature, and unconsolidated gouge) described in Introduction. Further modeling investigations () suggest that this shallow layer substantially reduces the potential for large coseismic surface rupture and accompanying large seismic wave radiation (unusual for natural earthquakes) in comparison with purely velocity-weakening models.

The along-strike distribution of (b-a) (Fig. 4D) shows a clear difference between the coseismic (~−0.005) and postseismic (~−0.01) areas. This change in (b-a) coincides with the transition between Cretaceous rocks to the north and younger Quaternary sediments in the south (Figs. 1A and 4A). As the unusual properties of the 2014 South Napa earthquake (shallow afterslip, position of the coseismic slip) are at least partially governed by this change in frictional rheology, the rupture propagation was clearly affected by the transition between the two lithological units. This division between Cretaceous rocks and Quaternary sediments happens only in the near surface region, while the rest of coseismic slip occurred at larger depths where the lithology is composed of Cretaceous rocks (). After the coseismic rupture propagates through this deeper area and arrives at the shallow layer, it continues only in the rock (northern) part of the fault, being impeded in the (southern) sedimentary part of the fault where a complementary afterslip develops subsequently. We suggest this mechanism to be responsible for the ribbon-like shape of the coseismic rupture.

Variability in the shallow postseismic slip

The evolution of shallow postseismic slip is spatially heterogeneous. Figure 4D shows the afterslip at three nearby points located from 10 to 15 km along strike. The temporal behavior varies in both amplitude and characteristic decay time. This is well constrained by the surface data and was also identified in kinematic inversions of afterslip (). In our dynamic model, the difference is facilitated by along-strike variations of (b-a) (see Fig. 4D). The value of (b-a) affects the time scales over which afterslip develops, as can be seen from a simple spring slider model (, , ), for which afterslip s(t) develops logarithmically with time t

In addition to (b-a), the temporal evolution of afterslip depends on effective normal stress σn, stiffness k (that scales with shear modulus μ and the inverse of patch size), and initial velocity vi. The normal stress and stiffness can be assumed constant in the horizontal direction (with the potential exception of lateral variations in fluid pressure that are beyond the scope of this paper), while the initial velocity is higher at the northern part, where the slip initiated during the coseismic phase.

We show the development of shallow afterslip in Fig. 4D, as calculated at three points near the surface (at 200 m depth) located from 10 to 15 km along strike. The positions were chosen to show the impact of different values of (b-a) changing from ~0 to −0.01 over 2 km. Afterslip starts much quicker close to the coseismic rupture where (b-a) is close to zero. The characteristic decay time of afterslip then clearly increases further to the south as (b-a) approaches −0.01. The afterslip develops under non–steady-state conditions in 3D models, and therefore does not entirely conform to the simplified logarithmic formula derived for a 1D spring slider, but its basic properties do hold. This short-distance variability in afterslip is a further example of the strong impact of fault lithology on rupture development. Whether it is driven by small-scale changes in mineral composition or pore pressure along the boundary between rocks and sediments remains an open issue.

Interplay between coseismic and postseismic rupture in the deep velocity-strengthening zone

The velocity-strengthening zone at 7.5 km depth (Fig. 3B) is a major finding of our modeling. The zone manages to rupture coseismically due to the lowered friction f0. Coseismic slip (and stress drop) is notably lower here than in other (velocity-weakening) parts, which is consistent with coseismic kinematic inversions (, ). Upward propagating coseismic rupture was followed by substantial deep afterslip (up to 0.4 m; see Fig. 2B) that also expanded out of the coseismic area. It is still concentrated to a relatively small patch, making its signature in the postseismic data relatively weak. Removing this afterslip patch from the model results in only a minimal change of the misfit (1 to 2%). On the basis of this, we suggest that the appearance of this velocity-strengthening zone is constrained by the dynamics of the coseismic rupture, whereas its afterslip is rather a by-product.

The deep afterslip can be indirectly corroborated by the appearance of off-fault aftershocks () that appear below the coseismic rupture with notable concentration around the area (Fig. 2B). Figure 4E shows the time development of the aftershock rate obtained by counting the aftershocks in the area outlined in Fig. 4B. The temporal decay of aftershock rate follows Omori’s law and is very similar to the evolution of stress rate obtained from the middle of the strengthening area, pointing to their possible driving by the deeper afterslip. While we use the aftershock rate to only confirm a stress trend in the strengthening zone, the addition of aftershock rate in the inversion directly as a measure of stress rate can be an additional piece of data to further constrain the postseismic model ().

MATERIALS AND METHODS

Friction law

In our model, coseismic and postseismic slip are governed by rate-and-state friction with fast-velocity-weakening ()

Equation 2 gives a value of friction S for given slip rate and frictional state variable Ψ. It is in the regularized form to avoid divergence at zero slip rate (, ), with only minor difference from the classical formulation for > 0 (). The evolution equation (Eq. 3) for the state variable Ψ is the slip law, in which the time derivative of the state variable is proportional to its distance to a steady-state value ΨSS and ratio of and characteristic slip L. The steady-state value is calculated in Eq. 4 from steady-state friction fSS as an inverse function of Eq. 2. The steady-state friction is defined by Eq. 5, where it decreases from low-velocity friction fLV to fully weakened friction fw with growing slip rate , as for due to the fast-velocity-weakening effect, following the flash-heating model (). The low-velocity steady-state friction coefficient fLV defined by Eq. 6 increases or decreases with slip rate following the sign of the difference between the state evolution (b) and direct effect (a) coefficients. The difference (b-a) in Eq. 6 thus distinguishes the velocity-weakening (b-a > 0) and velocity-strengthening (b-a < 0) modes of friction ().

Forward problem

We simulate the coseismic rupture with the code FD3D_TSN (). It uses a fourth-order finite-difference method to solve the 3D elastodynamic equation. The fault boundary condition (friction) is applied on a vertical fault with the traction-at-split-nodes method (). Free surface conditions are applied using a stress imaging technique (). We use perfectly matched layers () as absorbing boundary conditions. All computationally expensive routines are GPU-accelerated using OpenACC directives, yielding a speedup by a factor of 10 when comparing single GPU and single CPU runs. Accuracy of the code was tested () by using community Southern California Earthquake Center/U.S. Geological Survey (SCEC/USGS) benchmarks for both slip-weakening and fast-velocity-weakening friction laws (). Earthquake nucleation is induced by a second-long gradual increase of prestress in a circular zone. We use a spatial grid size of 100 m, providing a sufficient resolution of the cohesive zone, and a time step of 0.003 s satisfying the Courant–Friedrichs–Lewy stability criterion. The computational domain on one side of the fault is 10 km thick. Synthetic seismograms are obtained by convolving the resulting slip rates with Green’s functions precalculated using the Axitra code ().

Postseismic slip is simulated in a quasi-dynamic approximation, replacing the inertial term of the elastodynamic equation by a radiation damping on the fault (). We use a boundary element approach with a precalculated velocity-stress interaction kernel between fault nodes, assuming a vertical fault in a homogeneous medium (). This reduces the problem to a set of ordinary differential equations for displacements and state variables (–). We solve it by a Runge-Kutta method of fifth order with variable time steps on an undersampled grid with a 400-m spatial step. This quasi-dynamic postseismic modeling is used after the maximum slip rate in the finite-difference coseismic simulation falls below 1 mm/s. We tested the viability of the transition by postponing it by 10 s to 1 min, yielding only a negligible (below 1%) difference in the simulated long-term slip. Both predicted coseismic and postseismic GPS displacements are obtained by convolving the slip with precalculated Green’s functions. We note that the positions of the alignment arrays, NLAR, NWIT, NHNR, and NLOD, that measure the surface slip directly above the fault would not fit with our simplified planar geometry. Therefore, we artificially moved their positions to coincide with the position of the surface rupture on our planar fault, preserving their distance from the epicenter. We model the arrays as if they were GPS stations located at a 50-m distance from the fault with displacement equal to half of the measured slip.

Parameterization

The fast-velocity-weakening rate-and-state friction law involves a challenging number of potentially free parameters in the dynamic inversion, increasing the dimension of the model parameter space and increasing computational requirements. These include parameters of the rate-and-state friction a, b, f0, , and L; additional parameters governing the fast-velocity-weakening effect f and ; stressing conditions at the fault σn and T0; and initial values and Ψini. We assume a purely strike-slip fault so that T0 and are nonzero only in the horizontal direction. All parameters are thus spatially heterogeneous 2D scalar fields across the fault.

We use several relations and assumptions to limit the actual number of model parameters in the inversion and keep the inversion computationally tractable. Normal stress σn is set to be depth dependent, rising from 1 MPa at the surface to 100 MPa at 5 km depth and held constant at greater depth, where further depth increases in pore pressure and hydrostatic pressure are assumed to balance out (). Nonzero normal stress at the surface substitutes the cohesion we did not include directly in the modeling. Models with friction coefficient f are equivalent to models with cohesion c and friction coefficient f′ such that f = f′ + c/σn, provided that cohesion weakens in the same way as friction. At shallow depth (low σn), f0 > 1 can thus be accommodated with reasonable values of c (~1 MPa) and f′ < 1.

The fully weakened friction coefficient fw is set to 0.2, as observed in laboratory experiments (). Any other value of fw can be accommodated a posteriori by a straightforward modification of the results; the adjusted initial stress T0 would, in that case, be calculated by addition of the factor σn[fw(new) − 0.2] to the initial shear stresses constrained by our inversion.

The reference slip velocity is associated with a steady-state friction coefficient equal to f0. Since it is an arbitrary reference, we set it to 10−6 m/s as in other rate-and-state dynamic models [e.g., (, , )]. The initial value of the state variable Ψini is related to T0 and through Eq. 2. We calculate Ψini at the beginning from Eq. 1, following the approach in the SCEC/USGS benchmark TPV104 (). We fixed m/s and m/s in velocity-weakening (b-a > 0) areas of the fault. The former is supported by experiments, and the latter stems from the assumption that the coseismic region is locked before the onset of the earthquake. In contrast, in the velocity-strengthening areas, where the fault is supposed to creep at higher slip rates before the start of the earthquake (at least at ~10−10 m/s) (), we let free. Similarly, we let free in the strengthening zone to allow the rupture to stop.

In the end, the reduced set of dynamic model parameters to be determined by the inversion procedure are T0, a, b, f0, and L as 2D fields, and and as 2D fields on the smaller (velocity-strengthening) portion of the fault. For the purposes of the inversion, we parametrize the spatial distribution on an equidistant grid of 12 × 9 control points, from which the parameters are bilinearly interpolated onto the grids for the dynamic and quasi-dynamic simulations. The 2D fields are supplemented by four more free parameters describing our nucleation procedure realized by a 1-s-long gradual increase of prestress in a circular zone—the position of its center, its radius, and the added stress.

Inverse problem

We formulate the inverse problem in the Bayesian framework (, , ). We assume uniform prior PDFs for the model parameters in wide intervals of permissible values (Table 1). The data are considered to have Gaussian distributions of errors with SDs of 5 cm and 2.5 mm for seismograms and GPS, respectively. We sample the posterior probabilities using the Markov chain Monte Carlo (MCMC) parallel tempering algorithm (), accepting proposed models according to the Metropolis-Hastings rule. We used a modified version of the inversion code fd3d_tsn_pt. This code has been previously validated for slip-weakening friction law and only seismic data using synthetic tests () and applied to the 2016 Amatrice () and 2020 Elazığ earthquakes (). The present application required implementing the modified forward model and parameters.

Table 1.

Minimum and maximum values of prior uniform distributions of inverted parameters.

Note that and have uniform prior distribution in the velocity-strengthening regions only, being constant in the velocity-weakening areas.

Quantity	Label	Minimum value	Maximum value
Shear prestress (horizontal)	T 0	103 Pa	109 Pa
Direct effect parameter	a	0.001	0.1
State evolution parameter	b	0.001	0.1
Reference friction at velocity s·=s·_0= 10−6 m/s	f 0	0.1	2
Characteristic slip distance	L	0.1 m	2 m
Weakening velocity	s·_w	0.1 m/s	3 m/s
Initial velocity (horizontal)	s·_ini	10−13 m/s	10−7 m/s
Along-strike position of the nucleation	hx	14.5 km	16.5 km
Depth of the nucleation	hy	10 km	14 km
Radius of the nucleation patch	r nucl	400 m	1000 m
Stress increase in the nucleation patch	σnucl	1%	20%

We accelerated the inversion progress by starting from a reasonable model that was relatively homogeneous with velocity-weakening friction at the central square-shaped portion of the fault and velocity-strengthening friction on all edges. From there, we allowed the parallel tempering MCMC approach to explore the model space. We manually intervened several times by optimizing the prestress, nucleation, and frictional parameters to find a model with positive variance reduction. After that, we explored the model space by running the MCMC sampling on an IT4I cluster with four Nvidia Tesla V100 GPUs and in-house computers with three GPUs (Nvidia 2080Ti), with each forward simulation taking about 40 s in both cases. The total number of visited models was high (~500,000). The final set consists of ~7500 accepted models with a posterior probability density value larger than 5% of the pPDF maximum.