3.2.1 Stress-driven Afterslip Simulation
We carried out stress-driven afterslip simulations using the open-source software RELAX, which solves for the nonlinear time-dependents (x ,t ) in the Fourier domain under the assumption of rate-strengthening friction on faults (equation 1) (Barbot et al., 2009a). The afterslip evolution history on a given patch of the fault is controlled by the rate-strengthening constitutive law (Barbot et al., 2009a).
\(s\left(t\right)=\frac{\Delta\tau_{0}}{G^{*}}[1-\frac{2}{k}\coth^{-1}(e^{\frac{t}{t_{0}}}\coth\coth\ \frac{k}{2}\ )]\)(1)
In equation 1, \(k=\frac{\Delta\tau_{0}}{\text{aσ}}\) is the dimensionless ratio that controls the nonlinearity during the slip, and the time evolution is controlled by k along with the relaxation time \(t_{0}=\frac{1}{2V_{0}}\frac{\text{aσ}}{G^{*}}\text{\ \ }\). Note that the parameter a in the equations of Barbot et al. (2009), and as used here, is more commonly identified as (a -b ) in the context of full rate and state friction. Larger values of k result in models that are more strongly non-linear, with a more rapid decay in slip velocity early in the postseismic period (these models also require shorter time step and thus result in much longer program execution time). \(\Delta\tau_{0}\) refers to the shear stress perturbation due to the earthquake, σ refers to the effective normal stress on fault, and G* is the effective elastic constant per unit area determined by the linear dimension L and the shear modulus. The relaxation time\(t_{0}=\frac{1}{2V_{0}}\frac{\text{aσ}}{G^{*}}\text{\ \ }\)depends on both aσ and the reference sliding velocity on the fault \(V_{0}\). Total slip as goes to infinity is limited to \(\frac{{}_{0}}{{}^{*}}\).
Thus, there are 2 unknown parameters to search for to solve this problem: and \(V_{0}\). Many studies assume a value of and search only for \(V_{0}\), due to the strong parameter tradeoff between the two values when only one time period is considered (e.g., Tian et al., 2020). We first performed a 2-d grid search for and \(V_{0}\) over a relatively large range of parameter values to find the best fit values. We calculated the reduced\(\chi^{2}\) using the three sites on the Alaska Peninsula (AC40, AB13 and AB21) that are most sensitive to the downdip afterslip.
When we consider only one time interval, for example three months, then a very wide range of values, varying by orders of magnitude, yield models that fit the data equally well. Large values of (such as 3MPa suggested by Tian et al. (2020)) produce an afterslip evolution history at GPS sites like the orange curve in Figure 3, showing a low degree of nonlinearity, while small values of (similar to those used by Wang and Bürgmann (2020) or Zhao et al. (2022)) produce models like the blue curve in Figure 3, showing a higher degree of nonlinearity. Because the observations at 3 weeks more closely align with the curve produced by smaller values for (Figure 3 gray star), we limit the range of parameter values to those similar to those of Zhao et al. (2022) and consider displacement predictions for two time windows, 0-3 weeks after the mainshock and 0-3 months after. Based on the total misfit and given the nonlinear nature of the very early afterslip evolution, we fix the value of to be 0.6 MPa. Given that the two time windows we have used are short, using a different value of in our models would produce an equally good fit, with a correspondingly different \(V_{0}\) value. In this study, we vary the \(V_{0}\) value for each different model scenario that we consider in the following sections, and we leave the question of whether it is possible to determine an optimal value of to a future study with a longer time span.