Comparison of two time-marching schemes for dynamic rupture simulation with a space-domain BIEM

The boundary integral equation method (BIEM) is one of the important numerical techniques used to simulate geophysical phenomena including dynamic propagation, nucleation, and sequence of earthquake ruptures. We studied the stability and convergence of two time-marching schemes numerically for 2-D problems in Mode I, II, and III conditions. One was a conventional method based on piecewise-constant spatiotemporal distribution of the rate of displacement gap V\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ V $$\end{document} (CM), and the other was a slightly modified scheme from a predictor–corrector method previously applied to a spectral BIEM (NL). In the stability analysis, we simulated behavior of a traction-free fault under uncorrelated random distributions of initial traction. The growth rate of the perturbation is negative in a parameter regime of complex shape with CM, which has two numerical parameters, and the intersection for all the modes is very restricted as reported previously. In contrast, NL has only one parameter and yields simpler and a wide parameter regime of stability, conceivably allowing more flexible meshing on the fault. In the convergence analysis in which a smooth problem was solved, CM resulted in a numerical error scaled as Δx1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta x^{1} $$\end{document} while NL led to the scaling of Δx2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta x^{2} $$\end{document} typically or of Δx1.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Delta x^{1.5} $$\end{document} under certain conditions in Mode II problems. NL requires negligible additional computational costs and modification of the code is quite straightforward relative to CM. Therefore, we conclude that NL is a useful time-marching scheme that has wide applicability in simulations of earthquake ruptures although the reason for the rather complicated convergence behavior and verification of the findings here to more general conditions deserve further study.


Introduction
A boundary integral equation method (BIEM) for dynamic rupture simulation (e.g., Das and Aki 1977) has been widely used in the study of a variety of problems related to earthquake generation processes, from idealized 2-D cases (e.g., Andrews 1985) to 3-D realistic cases dealing with complex geometry composed of many fault planes (e.g., Ando and Kaneko 2018). In addition, the BIEM is a useful tool that can be used to study quasi-static problems (e.g., Okada 1985), the preparation processes of earthquakes (e.g., Dieterich 1992), and their sequences (e.g., Lapusta et al. 2000). It can also be extended to linear viscoelastic materials (Geubelle 1997;Kato 2002;Miyake and Noda 2019).
The BIEM makes use of analytic expressions of Green's function and does not suffer from artificial dispersion during propagation of elastic waves in a medium. Because we do not have to discretize the medium with the BIEM, it is possible to simulate a problem with higher on-fault resolution than with other methods such as finite difference methods and finite element methods, which are based on discretization of volume. Numerical techniques have been developed for or applied to BIEMs that dramatically reduce computational costs, such as fast Fourier transformation for spectral formulation (Lapusta et al. 2000), fast multipole decomposition (Tullis et al. 2000), a Barnes-Hut scheme (Beeler and Tullis 2011), a fast domain partitioning (FDP) (Ando et al. 2007), hierarchical matrices (H-matrices) (Ohtani et al. 2011), and a recently proposed method consisting of FDP, H-matrices, quantization, and averaged reduced time (Sato and Ando 2019). One of the weak points of BIEMs is that it is difficult to deal with a heterogeneous medium using them, but there are some approaches that extend the ability of BIEMs to such problems (Goto et al. 2010;Kame and Kusakabe 2012). Owing to these efforts, BIEMs are currently one of the most efficient and useful methods for studying fault dynamics in a range of problem settings.
There is, however, an issue with numerical stability in BIEMs (e.g., Fukuyama and Madariaga 1998;Tada and Madariaga 2001). For example, Tada and Madariaga (2001) systematically studied a numerical scheme based on spatiotemporally piecewise-constant basis functions for the rate of displacement gap V on the fault (Cochard and Madariaga 1994). They simulated a 2-D self-similar dynamic rupture of Mode I, II, or III, and concluded that the region of stability in the numerical parameter space is rather restricted and depends on the mode of the rupture (Figs. 5, 6 and 7 in Tada and Madariaga 2001). The intersection of the stability regimes for multiple modes is so restricted that the numerical parameters must be chosen carefully for simulation of mixed-mode ruptures. Furthermore, the ability of flexible meshing on the surfaces is very restricted, as the numerical parameters vary locally for a heterogeneous mesh. It is often the case that a small artificial damping factor is added (e.g., Yamashita and Fukuyama 1996;Kame et al. 2003) to enhance the stability of the simulations. The effect of the damping should be small, but could make solutions of insufficient resolution artificially smooth and good-looking if simulations were to be tuned unfairly. Therefore, researchers must be careful to avoid misinterpretation of numerical results. Noda and Lapusta (2010) suggested a different, secondorder time-marching scheme based on a predictor-corrector approach for a spectral, wavenumber-domain BIEM (SBIEM) by slightly modifying a scheme proposed by Lapusta et al. (2000) and Lapusta and Liu (2009). In the SBIEM, temporal convolution is calculated with a midpoint integration scheme, and V at temporal midpoints is estimated by linear interpolation between collocation points. In the present study, we have adapted the midpoint evaluation of V to the space-domain BIEM based on spatiotemporal piecewise-constant basis functions and found that the present scheme was more stable and accurate than the conventional one.
The outline of this paper is as follows. The two timemarching schemes investigated are briefly described in the "Time-marching schemes" section. The procedures of numerical experiments conducted to study the stability and accuracy of the schemes are described in the "Methodology" section, and the results of these experiments are reported in the "Results" section, which is followed by "Discussion" and "Summary" sections.

Time-marching schemes
In the present study, we investigated time-marching schemes in the simulation of dynamic rupture propagation for a fault embedded in an infinite homogeneous linear elastic medium. For simplicity we considered 2-D problems with a planer fault subjected to either Mode I, II, or III loading conditions. In such problems, a relevant component of traction on a fault τ is expressed as spatiotemporal convolution of an integration kernel K and history of the rate of displacement discontinuity V in a boundary integral formulation (Cochard and Madariaga 1994;Geubelle and Rice 1995), where x is the coordinate along the fault, t is time, τ ini is the initial value of the traction on the fault if V = 0 . t − is infinitesimally smaller than t , and the difference of t and t − is crucial in the proposed scheme. K is a hypersingular kernel or can be expressed as a combination of nonhypersingular kernels and derivative operators (Cochard and Madariaga 1994;Tada and Yamashita 1997). η represents the amount of an instantaneous elastodynamic effect due to impedance of the medium extracted explicitly from the convolution. It depends on the mode of the rupture, where ρ is the density, c s is the S-wave speed, and c p is the P-wave speed. Cochard and Madariaga (1994) discretized Eq. (1) based on piecewise-constant spatiotemporal distribution of V . The collocation points where τ and V are defined are indicated in Fig. 1a. The traction after the n th timestep and at the i th element is estimated by the summation of the effects of the sources at the m th timestep and the j th element, n = 0 corresponds to the first snapshot defined at t = e t t . N is the total number of the elements, x is

Conventional time-marching scheme
(2) 2η = ρc p (Mode I) ρc s (Mode II and III) , the length of the elements, h T is the nondimensional timestep h T = c s �t/�x , and e t is the nondimensional timing of the collocation points (Fig. 1a). V n i is defined at t = (n + e t )�t so that the collocation point is in the interior of the support of the corresponding basis function if 0 < e t < 1 . Note that the kernel value depends only on the relative location and time difference between a source and a receiver in our simple problem setting. The instantaneous term in Eq. (1) is included in the summation by assuming that V does not change suddenly at the collocation points. For expression of the kernel, please see Appendix in Tada and Madariaga (2001).
Suppose that the values of V m j are known for all j and m ≤ n − 1 . Then Eq. (3) yields N equations with 2N unknowns ( τ n i and V n i , 1 ≤ i ≤ N ). In problems dealing with rupture propagation, a fault constitutive law is assumed. As it must satisfy the principle of local action, it provides additional N equations that relate τ n i and V n i . Therefore, we can solve for τ n i and V n i , if they exist, and carry on the simulation. This scheme shall be abbreviated as CM in the present paper.
CM has two nondimensional numerical parameters h T and e t , and Tada and Madariaga (2001) demonstrated that the numerical stability depends on the selection of them, and on the mode of the rupture. Their results are compared with our numerical experiments in a later section. Noda and Lapusta (2010) proposed a second-order time-marching scheme based on a temporal midpoint integration scheme of Eq. (1) formulated in a wavenumber domain for use with the SBIEM (e.g., Lapusta et al. 2000). In the SBIEM, Fourier basis is used so that the spatial distribution of V is assumed to be smooth. Such representation is available only for a planer fault, which is recognized as a weakness of the SBIEM. Here, we applied a time-marching scheme similar to that by Noda and Lapusta (2010) to the space-domain BIEM.

Present scheme based on a predictor-corrector method
The spatiotemporal discretization in the present schema is illustrated in Fig. 1b. We approximate Eq. (1) using spatiotemporally piecewise-constant basis functions for V similarly to CM, but estimate the history of V at temporal midpoints indicated by crosses, and allow for sudden change in V at the collocation points indicated by circles.
V n i and V m+1/2 i are defined at time n t and (m + 1/2)�t , respectively. It should be noted that we use the same integration kernel as in CM. 1+ in Eq. (4) is infinitesimally larger than 1 , and explicitly indicate that the collocation point for τ n i and V n i is out of the support of the basis function for V n−1/2 i . Because the timing (4)

a b
Fig. 1 Schematic diagrams of spatiotemporal discretization for time-marching schemes investigated in the present study. a Cochard and Madariaga (1994) (CM) and b the present scheme (NL). Collocation points where V and τ are defined are indicated by circles, and V at the points indicated by crosses are stored and used in numerical integration. Piecewise-constant approximation was adopted in both schemes, and a basis function is indicated by a gray square. Thick dashed line in b indicates that the collocation points (circles) are not included in the support of the basis function of the collocation point relative to the basis function is fixed, this discretization has only one nondimensional parameter h T . After the n − 1 th timestep, we know the V m i values where m ≤ n − 1 from which we can calculate the contribution from the previous history to the functional term, The traction at the end of the current timestep is expressed as where In a predictor step, we assume V n−1/2 j to be identical to V n−1 j , and estimate φ head i , This is a first-order approximation if the solution is smooth. We can then obtain an approximation to τ n i and V n i by solving a friction law and simultaneously. If we adopted τ n * i and V n * i as the final estimate and V n−1/2 j were not updated, the scheme would be CM with e t = 0+.
In a corrector step, we first obtain a better estimate of V n−1/2 j based on linear interpolation, and re-calculate the contribution of the current timestep to the functional term, Then, we update the estimate of τ n i and V n i by solving the friction law and Finally, we adopt these updated values as a numerical solution, and This scheme shall be labeled as NL in this paper. If the corrector step is iterated until the numerical solution converges, the scheme becomes implicit, which solves Eq. (4) and the friction law simultaneously. For a relatively short timestep ( h T c s c p ≤ 0.5 for Modes I and II, and h T ≤ 0.5 for Mode III), φ head i is zero regardless of V n−1/2 j . Therefore, the numerical solution converges after the first corrector step.
The computational cost is dominated by calculation of spatiotemporal convolution. Although this scheme consists of two sub-steps, the main part of the convolution, summing up contributions from the previous steps, is executed only once per timestep. Therefore, the additional numerical resources required are negligible relative to CM.
V is assumed to be constant in the predictor step. We can interpret the corrector step as two-step approximation. Firstly, V is assumed to change linearly with time, which would yield global error scaled as t 2 . Secondly, the mid-point evaluation of V (Eqs. 11 and 16) is applied and the same kernel as for CM is used instead of preparing another integration kernel for the linearly evolving V . This second approximation also yields global error scaled as t 2 where K is smooth enough (Appendix). Just note that in the SBIEM by Noda and Lapusta (2010), both V and the integration kernel are evaluated at the temporal midpoints so that the temporal integration is estimated with the midpoint scheme. In the numerical experiments presented later, NL was shown to be useful owing to its enhanced numerical stability, faster convergence of numerical error, and minimal modification of the simulation code required relative to CM.

Methodology
To investigate the stability and accuracy of the two timemarching schemes, we conducted two types of numerical experiments for Mode I, II, and III problems. The first being the numerical simulation of evolution of a random field mimicking a numerical error, and the second a convergence test for a smooth problem. To avoid additional complexity from the friction law, the fault is assumed to be traction-free ( τ n i = 0 ) in the experiments. Such a simple friction law yields Traction-free surfaces appear frequently in simulations in practice. Hok and Fukuyama (2011) used traction-free elements to include a ground surface in their model. In dynamic rupture simulations with slip-weakening friction laws, a ruptured region of residual, constant frictional strength can be considered traction-free if we focus on the response in the fault motion to perturbation. Therefore, the stability condition obtained in the present numerical experiments can be used as a necessary condition for a wide range of dynamic rupture simulations. The problem settings of the numerical experiments are described in detail in the following subsections.

Numerical stability analysis
For integral-differential equations, an analytic approach to assess numerical stability such as Von Neumann stability analysis may be difficult. In the first numerical experiments, we numerically simulated evolution of initially random perturbation added to a simple analytic solution, Poisson's ratio ν = 0.25 was assumed, which is a good approximation to rocks. The fault was discretized by N = 64 elements, and the fault behavior was simulated until c p t x = 320 . τ ini i was given by Gaussian random distribution around zero with standard deviation of σ , and the same realization of τ ini i was used for all simulations (Fig. 2a). In the simulation, τ and V were nondimensionalized by σ and σ/ρc s , respectively.
During the simulations, the L 2 norm of nondimensionalized V , changes with time. The growth rate was calculated for each simulation by fitting a linear envelope to a plot of log 10 ||V ρc s σ 2 || as a function of c p t x in the latter half of a simulation, 160 ≤ c p t x ≤ 320 , and defined as the slope of the envelope. If the perturbation grows, it is conceivable that the random round-off error increases with time, and the numerical solution cannot be regarded as an approximation to the true solution.
With CM, a 2-D parameter study was conducted for 0.1 ≤ h T ≤ 1.5 and 0 ≤ e t ≤ 1 every 0.01 ( 141 × 101 runs for each mode). With NL, there is only one numerical parameter. Therefore a 1-D parameter study for 0.1 ≤ h T ≤ 1.5 every 0.01 ( 141 runs for each mode) was performed.

Convergence test
Convergence tests were performed for a problem of smooth distribution of τ ini (x) applied to a traction-free fault in Mode I, II, or III conditions. We assumed τ ini to be an infinitely differentiable function, where �τ is the amplitude and is the approximate width of the perturbation (Fig. 2b). Note that the maximum value of the gradient dτ ini /dx is 4�τ/ , and thus it can be anticipated that �x < /4 is required for simulations with a good resolution.
A series of simulations were conducted until tc s exceeded 2 with various x from to /256 with CM and NL. The system size was set long enough that the ends of the fault did not matter. The mesh was refined by adding extra collocation points in the middle of neighboring ones. Therefore, positions of the collocation points in a mesh always form a subset of those in the finer meshes. The numerical parameters e t and h T must be selected from the region of stability obtained in the numerical stability analysis, which is explained in the following section. For CM, we selected a couple of parameter sets, (h T , e t ) = (0.5, 0.4) and (0.9, 0.55) . As h T increases, the amount of memory required for simulation of the same physical time duration decreases. h T = 0.9 is around the upper maximum of the stability condition as explained later. For NL, we tested h T from 0.4 to 1.1 every 0.1.
In addition to CM and NL, we carried out a convergence test with the SBIEM and the second-order time-marching scheme by Noda and Lapusta (2010) for reference with h T = 0.5 . Owing to calculation of the spatial convolution in the wavenumber domain, the SBIEM requires much less numerical resources and thus we can refine the mesh much more than with the space-domain BIEM. Simulations with the SBIEM were conducted with x from to /1024 , and the simulation with the best resolution ( �x = /1024 ) was used as a reference solution in estimating the numerical error.
Numerical error E was defined as the L 2 norm of the difference in V from the reference solution at tc s / = 2 . If the collocation time points did not exist exactly at tc s / = 2 , then V at tc s / = 2 was estimated by temporal linear interpolation between collocation points. E was defined as where V i is the value of V at tc s / = 2 and at the i-th grid point, V R r(i) is the corresponding value of V in the reference solution. E(�x/ ) is a decreasing function for a consistent numerical scheme if x ≪ 1 , and the order of accuracy can be recognized as the slope in a log-log plot. Figure 3 shows representative examples of evolution of V ρc s /σ 2 in the Mode II problem using CM with numerical parameters e t = 0.6, and h T = 0.65, 0.7, 0.75, and 0.8 . The dashed lines in c p t �x > 160 are the envelopes that define the growth rate. Just after initiation of the simulations, the norm of V decreased regardless of the parameters. It began to increase after a while in some cases, showing exponential growth (cases with h T = 0.65 and 0.7 ). In other cases, it sharply decreased at around

Numerical stability analysis
c p t x = 120 , and continued to decrease exponentially. The cases with the smallest growth rate values were almost the same each other and were very close to the plot for h T = 0.8 in Fig. 3.
Spatiotemporal distributions of V for the cases plotted in Fig. 3 are shown in Fig. 4. Note that V at each timestep is normalized by the L 2 norm of V at the same timestep (Fig. 3). The normalization helps us recognize dominant pattern of each numerical solution. With h T = 0.65 (Fig. 4a), alternating oscillation became dominant soon after the initiation of the simulation, and propagation of waves from the initial heterogeneous field was hardly seen. With h T = 0.7 (Fig. 4b), the oscillation grew more slowly, and the waves could be seen until about c p t x = 60 at the center of the fault, and about c p t x = 120 near the ends of the fault. Note that Nc p /c R = 120.57 for ν = 0.25 where c R is the Rayleigh wave speed. After passage of these waves, the alternating oscillation became dominant. With h T = 0.75 (Fig. 4c) and 0.8 (Fig. 4d), it was clear that Rayleigh waves radiated from the initial stress distribution that passed through the fault by c p t x = 120 . Subsequent fault motion apparently consisted of a series of Rayleigh waves evidenced by checkered patterns at large c p t x . The pattern seemed to be richer in high wavenumber and frequency for h T = 0.75 than for h T = 0.8.
The growth rate obtained with CM for each mode and for a range of h T and e t is indicated by color maps in Fig. 5. The color scales are selected from the minimum value to −2 times the minimum value so that variation in largely positive values is not indicated. The shape of the region of stability enclosed by white-black dashed lines is complex, and different for different modes as reported by Tada and Madariaga (2001). Their results are plotted  Fig. 5 for comparison and will be discussed later. Figure 6 shows the stability fields for all combinations of the modes. In problems comprising a mixed-mode rupture or with surfaces where motion in multiple directions is allowed, the numerical scheme must be stable for multiple modes, while the intersection of the stability regimes of all the modes is quite restricted. Therefore, the numerical parameters must be selected carefully for CM. The growth rate obtained with NL is indicated in Fig. 7. It is remarkable that the scheme is stable for a wide range of h T , which is the only parameter in the scheme. For Modes I and III, the scheme is stable almost everywhere as long as studied. For Mode II, it is less stable, but there is a wide region of stability of 0.40 ≤ h T ≤ 1.14 in addition to a smaller one of 0.26 ≤ h T ≤ 0.29 (see Fig. 7b). Intersection of the regions of stability for all modes is identical to that of Mode II. The wide continuous region of stability allows us to vary the size of elements by a factor of 2.85 if the length of the timestep is spatially uniform in a simulation.

Convergence test
The results of the convergence tests for Mode I, II, and III conditions are shown in Fig. 8a, b, c, respectively. For all the modes, CM results in a first-order scheme, while NL shows faster convergence rates. With NL E x depends on h T by about one order of magnitude so that selection of h T is important for good quality simulations if accessible numerical resources allow. It is surprising that NL with h T = 0.5 is slightly better than the SBIEM with the same h T for Modes I and III.
To see detailed behavior with NL, we normalize E x by x 2 indicated by thick gray lines in the upper panels of Fig. 8. Horizontal dependency in the lower panels of Fig. 8 indicates second-order accuracy. It is typically the case that a smaller h T value (i.e., shorter timesteps and larger memory requirement) results in a smaller numerical error. In Mode II, however, the cases with h T = 0.4 and 0.5 show scaling of only x 1.5 , and the numerical error can be worse than in cases with larger h T in simulations with a very good resolution (i.e., small x ).

Discussion
In addition to the results of the present study, numerical stability reported by Tada and Madariaga (2001)  It is apparent that the stability shown by Tada and Madariaga (2001) anticorrelates well with the growth rate calculated in the present study. Therefore, we conclude that the numerical error studied by them is mainly due to growing oscillation in the ruptured, traction-free part of the self-similar ruptures. Note that a smaller h T requires a larger number of timesteps per unit time. Therefore, if two cases have the same growth rate defined in the present study, the case with smaller h T should show less poor stability evaluated by Tada and Madariaga (2001) than the other case. With NL, the observed convergence rate x 1.5 for Mode II with h T ≤ 0.5 was smaller than that of the other cases that showed x 2 . The integration scheme in NL (Eq. (4)) is second-order accurate for smooth enough functions, while the kernel of the space-domain BIEM has singularities at the wave fronts and at the source location. These singularities could be the source of the slower convergence rate, for the singularity of order −1/2 at the wave fronts results in global error scaled as x 1.5 . However, the reason why the other cases show faster rate of x 2 is not clear and deserves further study. This problem may be less severe for 3-D cases in which the singularity is milder. Note that this issue does not exist in the SBIEM, which uses a smooth integration kernel expressed in the wavenumber domain. Figure 8 illuminates additional complex dependency of the numerical error on x and t . Comparison of the cases with the same x and different h T and thus different t shows non-monotonic behavior for Mode II with �x ≤ /16 . The numerical error of the cases along the x 1.5 dependency for Mode II (thick dashed lines in Fig. 8e) decreases by a factor more than 3 with decreasing t by 20% with keeping x , yielding scaling of about t 5 locally. The cases with h T = 0.5 have the same t as those in the cases with half x and h T = 1.0 . Comparison of them almost always indicates that the numerical error increases with decreasing x with keeping t . Those observations cannot be understood based on convergence of numerical integration of a regular function. Further investigations on the accuracy in numerical treatment of the hypersingularity or on the effect of interaction of V and φ on the fault are required to fully understand the present and conventional time-marching schemes to verify, not only experimentally, but also theoretically, the accuracy of dynamic rupture simulations. In addition, the present study deals with the simplest 2-D problem consisting of a planer fault and a trivially simple friction law. Applicability of the results here to more general cases such as 3-D problems with non-planer faults has to be tested.
On the friction law, Lapusta et al. (2000) proposed an adaptive time stepper based on a stability analysis of the quasi-static limit with a rate-and state-dependent friction law (e.g., Ruine 1983), followed by a quasidynamic stability analysis with an exponential time-differencing scheme by Noda and Lapusta (2010). Positive ratedependency, either instantaneous or at steady-state, enhances numerical stability. On the other hand, negative instantaneous rate-dependency may lead to an ill-posed problem (Rice et al. 2001). We would like to point out that the friction law significantly affects numerical stability although investigation of this issue is out of the aim of the present study.
Although the results of numerical experiments reported here have not fully understood, we found NL useful and applicable to many geophysical problems. It requires negligible additional numerical resources compared with CM and allows us more flexible meshing on a fault for problems with a mixed-mode rupture. NL has only one numerical parameter h T , and thus tuning of a simulation is simpler. In addition, NL produces a smaller numerical error than CM, and thus requires less computational costs for simulations of the same quality.

Summary
We tested two time-marching schemes used in numerical simulations of dynamic rupture based on a boundary integral equation method (BIEM). One (CM) is a conventional method based on spatiotemporally piecewiseconstant distribution of V , the stability of which was previously investigated by Tada and Madariaga (2001). The other (NL) is based on a predictor-corrector method previously used in a spectral, wavenumber-domain BIEM (SBIEM) (Lapusta et al. 2000;Noda and Lapusta 2010) and adapted to the space-domain BIEM in this study. CM has two nondimensional numerical parameters e t and h T , which represent the timing of the collocation points and the length of timesteps, respectively, while NL has only one parameter, h T . Two types of numerical experiments have been conducted for 2-D Mode I, II, and III problem settings for a flat traction-free fault.
Numerical stability was numerically evaluated by simulating fault motion due to randomly distributed initial stress τ ini . We found that the growth rate with CM anticorrelates very well with the report by Tada and Madariaga (2001), and that the stable region in the parameter space (h T , e t ) is rather complex. The intersection for all the modes is quite restricted, implying that mixed-mode ruptures are difficult to simulate with CM. However, in contrast, NL shows a wide range of stability in the parameter space. The growth rate is negative for all the modes if 0.26 ≤ h T ≤ 0.29 or 0.40 ≤ h T ≤ 1.14 . Therefore, NL is probably useful in simulating mixed-mode ruptures or in meshing a fault heterogeneously. Convergence tests conducted for initially smooth distribution of τ ini show faster convergence rates for NL than CM. CM shows first-order convergence, and overall NL shows second-order convergence except in a couple of Mode II cases with h T = 0.4 and 0.5 , in which the numerical error is scaled as x 1.5 . Although the reason for such complex dependency deserves further study and applicability of the findings to more general cases has to be confirmed, we conclude that NL is a useful scheme because of its enhanced stability, smaller numerical error, and negligible additional computational cost compared with CM.