Modeling and Forecasting Aftershocks Can Be Improved by Incorporating Rupture Geometry in the ETAS Model

We implemented an extended version of the space‐time Epidemic‐Type Aftershock Sequence (ETAS) model, which simultaneously incorporates earthquake focal depths and rupture geometries of large earthquakes, and applied it to the 2016 Kumamoto earthquake sequence. Results show that the new model corrects the estimation biases of model parameters in the point source ETAS model. The reconstructed patterns of productivity density of aftershocks, along the mainshock rupture plane, form complementary patterns for coseismic slip in space and show significant spatiotemporal migrations. Large aftershocks tend to nucleate at the edges of high productivity density areas. The decay of direct aftershocks near the mainshock rupture is consistent with static stress changes caused by the mainshock. The forecasting capability of the ETAS model can be enhanced by considering the rupture geometries of mainshocks, especially for forecasting aftershocks in the first 1–2 days. In the simulation test, the incorporation of focal depths improves the forecasting resolution.


Introduction
Understanding the patterns of aftershocks following large earthquakes is important because aftershocks can also cause serious damage, and the occurrence times and locations of aftershocks provide us with information about ongoing stress adjustments in the fault zones of mainshocks. With new developments in observational seismology and geodesy over the past few decades, it is possible for us to acquire source process data within hours of a mainshock. Joint inversions of surface deformation, and ground motion data, from multiple sources provide improved constraints and high-resolution information on the distribution and evolution of fault slip that helps resolve the mainshock rupture process. However, past studies on the relation between mainshocks and aftershocks are limited to visual comparisons between aftershock hypocenters and slip distributions of mainshocks (Das & Henry, 2003;Ebel & Chambers, 2016;Mendoza & Hartzell, 1988;Woessner et al., 2006). Therefore, quantitative studies, using sophisticated aftershock models, are needed to elucidate details on the generation of aftershocks and how these are associated with the mainshock ruptures.

10.1029/2019GL084775
The space-time Epidemic-Type Aftershock Sequence (ETAS) model (e.g., Console & Murru, 2001;Console et al., 2003;Ogata, 1998;Zhuang et al., 2002Zhuang et al., , 2005 has proven effective in quantifying the triggering relationship between foreshocks for mainshocks and mainshocks for aftershocks. In this model, earthquake sources are usually treated as point sources and are represented by their epicenters for simplicity. The scaling effect with magnitude is incorporated by introducing an exponential term to the isotropic spatial kernel. Such simplification leads to biases when fitting the model to earthquake sequences originating from large mainshocks. Previous studies show that treating earthquakes as point sources, and assuming an isotropic distribution of aftershocks, results in the underestimation of the productivity of large mainshocks (Bach & Hainzl, 2012;Hainzl et al., 2008). Thus, the space-time ETAS model yields a biased classification of triggering family trees in an aftershock sequence, that is, many direct aftershocks from the mainshock are recognized as high-order aftershocks (aftershocks of aftershocks; Zhuang et al., 2018). One way to improve the modeling of direct aftershocks is by assuming that aftershocks decay with distance to the mainshock fault plane (Bach & Hainzl, 2012;Marsan & Lengliné, 2010). To address this, Guo et al. (2015aGuo et al. ( , 2017 proposed and verified the finite-source ETAS model and demonstrated its feasibility through applications. Besides the finite-source ETAS model, Guo et al. (2015b) proposed a 3-D point source ETAS model, which assumes that the focal depths of triggered aftershocks are decoupled from epicenters and that they follow the Beta distribution. Zhuang et al. (2018) applied three versions of space-time ETAS model (2-D, 3-D, and 2-D finite source) to the modern Italian catalog and found that considering the depth component significantly improved the data fit. They concluded that the focal depth should be incorporated into the finite-source ETAS model.
In this study, we combine the finite-source ETAS model with the 3-D point source ETAS model so that both the focal depth and the rupture geometry are incorporated in model formulation. This model is applied to the Kumamoto aftershock sequence, and the results are compared with those from three other versions of the ETAS model. In addition, an experimental short-term aftershock forecast is carried out to test the forecasting ability of the new ETAS model.

Methodology
In the 2-D and 3-D ETAS models, the conditional intensity function can be written in the forms of with  t representing events before the study period and beyond the study region which may have an influence on events in the target period and region and (x, y, z) denoting the longitude, latitude, and depth coordinates. The background seismicity in both models is assumed to be temporally stationary and spatially inhomogeneous. (M i ) denotes the expected number of direct aftershocks triggered by an event of magnitude M i , given by where M c is the completeness magnitude. The temporal probability density function (p.d.f.) g(t − t i ) is normalized, following the Omori-Utsu law of (Omori, 1894;Utsu, 1969) A specific form for the epicenter p.d.f. provided by Zhuang et al. (2005) can be written as where S i = (x i , y i ) is the point source of event i. Following the rationale discussed in Guo et al. (2015b) and Guo et al. (2018), we use the Beta distribution as the p.d.f. of the depth component h(z; S i ): Note that we add the term h 0 (z) in equation (1) to make it comparable to the 3-D ETAS model, and h 0 (z) is assumed to follow a distribution estimated by the histogram of all depths.
In equations (2), (5), and (7), the S i in the spatial kernel has different meanings for small and large events. Small events are treated as point sources and represented by their hypocenter locations as recorded in the catalog in equation (5). To incorporate the rupture extensions of large events (with magnitude over 6.0 in this study), the spatial p.d.f. in the 2-D and 3-D finite-source ETAS models should be Following Guo et al. (2017), we use an iterative algorithm to simultaneously estimate the background seismicity and productivity distributions of the finite sources. Please refer to the supporting information for detailed steps of this algorithm.

Data and Computation Settings
The Kumamoto earthquake of 15 April 2016 was one of the most notable earthquakes in recent years on the Island of Kyushu, Japan. Two foreshocks, with magnitudes greater than 6.0, occurred in the days preceding the mainshock. The mainshock was a major strike-slip event with a normal slip component (strike/dip/rake: 222 • /77 • /−161 • ) as per the Global Centroid Moment Tensor solution (Ekström et al., 2012). It initially ruptured near the locations of the two large foreshocks, propagated unilaterally to the northeast, and terminated near the volcanic area beneath Mount Aso. We used data from the earthquake catalog of the Japan Meteorological Agency (http://www.data.jma.go.jp/svd/eqev/data/bulletin/hypo.html, last accessed: October 11, 2018) for the period from 2010 to 2018 within the region 130-132 • E, 32-33.5 • N. The epicenter distribution is shown in Figure 1a. The minimum magnitude threshold is taken as magnitude 2.0 so that we can include as many events as possible in our analysis, although some events smaller than magnitude 3.0 seem to be missing in about the first 0.5 hr immediately following the mainshock, as shown by the blank area in Figure 1c.
In setting up the finite sources for the two foreshocks and the mainshock, we take early aftershock clusters as constraints of their boundaries, as shown by the white and red polygons in Figure 1a. Note that the initial source region that contains direct aftershocks from the mainshock is several times larger than the mainshock rupture. Such an extension is needed to model the off-fault aftershocks that are triggered by static or dynamic stress changes. The source volumes of three major earthquakes are divided into 0.01 • × 0.01 • × 1 km grids. Based on the depth distributions of earthquakes before and after the mainshock, the maximum depths of the initial finite sources for the two foreshocks and the mainshock are 15 and 20 km, respectively.

Data Fitting
The fit results for four models are included in Table 1. It shows that the 3-D finite-source model yields the highest likelihood, meaning that its overall fitting is significantly better than the other models. Also, the parameters are quite different between the point source and the finite-source models. Especially, higher q values from the finite-source models mean that the aftershocks are more concentrated around the mainshock source. Another significant difference is related to the parameter pair, and A. The value increases from 1.20 to 1.79 in the 2-D case and from 0.89 to 1.73 in the 3-D case. As a trade-off effect, the A value, which represents the productivity of an event at magnitude M c , decreases. As implied in equation (3), a higher value indicates greater differences between small and large earthquakes in triggering aftershocks. In the point source models, the aftershock productivity of large shocks is underestimated because their triggering ability is limited by an isotropic kernel circling around the focal epicenters. It gets worse in the 3-D point source model because the triggering ability of large shocks is further limited by the depth component. Therefore, many triggered events are classified as higher-order aftershocks in the point source models, as implied by the high A and low values. In this sense, the finite-source models correct the biased classifications of family trees in aftershock sequences (Guo et al., 2017;Zhuang et al., 2018).
It is worth noting that incorporation of finite sources corrects the underestimates in the aftershock productivity of the two foreshocks and the mainshock. However, the triggering abilities of earthquakes with Note. log L = logarithm of likelihood. moderate magnitudes, especially for events between magnitude 5.0 and 6.0, are still underestimated because their rupture geometries are not considered in the calculation, as shown in Figures S1 and S2.

Productivity Density Pattern
The total productivity values of the two foreshocks and the mainshock from the 3-D finite model are 273.0, 196.4, and 1,598.3, respectively. They represent 3.6%, 2.7%, and 23.3%, respectively, of the total number of events that occur after the two foreshocks and the mainshock within the study region. The productivity density pattern of the first foreshock is quite different from that of the second one as shown by the 2-D finite model ( Figure S3). The first foreshock has high productivity density in the middle of the source area, while the second one has high density in the southwestern part of the same source area. From the locations of large events (M5.0+) in Figure 2 (also see Figure S3 for M5.5+ events), we find that such events tend to be located at the edges of areas with high density. This is consistent with the results obtained from other studies (Stallone & Marzocchi, 2019;van der Elst & Shaw, 2015), which find that larger aftershocks tend to nucleate at the edges of the highly clustered aftershock zone. Such a phenomenon can be interpreted as resulting from stress relaxation along either the mainshock rupture or the surrounding subfaults. However, this feature is not included in the framework of current versions of the ETAS model. Although large triggered earthquakes tend to occur at the edges of high-density areas, these events are more likely to be located inside high-density areas based on the spatial kernel of the finite-source ETAS model. However, further study can be carried out on the distribution of aftershock magnitudes, and associated productivity density patterns, to verify this conclusion. This may enable improved capability for the forecasting of large aftershocks.
The productivity density patterns of the mainshock at different depths obtained from the 3-D finite model are shown in Figure 2. The extended source of the mainshock is divided into three subregions according to the locations of aftershock clusters. The high density of aftershock productivity of Region C is mainly distributed between 0 and 10 km, while the high density of productivity in Region B is located at the moderate depth range from 6 to 10 km. In Region A, a high potential for direct aftershocks occurred across almost all depth layers.

Spatial Decay of Direct Aftershocks
To measure the spatial decay of direct aftershocks with distance from the mainshock fault plane, we extract the linear density of direct aftershocks from the productivity density pattern in Figure 2. This is done by calculating the distance of each patch to the nearest subfault (slip > 1 m) of a slip model and summing up all productivity density values of the patches with distances in the same length bin Φ(d) = ∑ i I(d ∈ [d, d + Δd)), where d denotes the distance of patch and I( * ) takes 1 when the statement is true and 0 otherwise. The result is shown in Figures 3a1 and 3a2. Here we consider three slip models to account for the uncertainty in the inversions. Although the modeled slip distributions (Figures 3c1-3c3) are different at fine scales, the decay of the linear aftershock density is similar for each, especially at distances greater than 30 km. In Figure 3a2, the decay of the aftershock density has two stages: (a) with d −0.5 from 2 to 10 km and (b) with d −1.8 approximately from 10 to 30 km. Such results are consistent with the static stress triggering test in Hainzl et al. (2010), in which they find that the existence of a free surface for a crustal earthquake of magnitude 7 can lead to a statically triggered aftershock decay with d −0.5 in the near field (within 60 km). However, we obtain a faster decay in the range of 10-30 km. We explain this discrepancy in the two following ways: First, in Hainzl et al. (2010), the authors present relatively simple simulations of static stress changes without considering the heterogeneity of slip and the focal mechanisms of aftershocks. Second, aftershocks near Mount Aso are in the range of 10-30 km, and the faster decay may be associated with the high heat flow associated with geothermal areas. In Figure 3a1, the two peaks of high productivity values at 50 and 60 km correspond to dynamically triggered earthquakes in Region C (Figure 1a).
The distribution of direct aftershocks interpolated on the rupture faults of the three slip models (Figures 3b1-3b3) forms complementary patterns to coseismic slip. Specifically, in the middle and southwestern parts of the mainshock fault, the high density of aftershock productivity are broadly distributed between 5 and 20 km. Large coseismic slip exists in the up-dip direction of this area. At the northeast end of rupture fault, there is a high productivity asperity where the coseismic rupture terminates and a gap of ∼8 km length, with low productivity values, in between. This gap is consistent with the 10-km-long gap of aftershocks reported in Yue et al. (2017) who attribute this phenomenon to partial melting in this volcanic region.

Forecasting Experiment
To understand the forecasting capability of the improved ETAS models, we carry out a forecast experiment for all four models. Details on the forecast settings, results, and the influence of depth uncertainty can be found in the supporting information. As Figure S6 shows, the 3-D models perform much better than the 2-D models within the source region, indicating that the incorporation of the depth component enhances the resolution of aftershock forecasts. However, the finite-source models yield slightly worse results than the point source models for the forecast period. Such a phenomenon indicates that the finite-source model has a limitation in aftershock forecasts: We need some location information of early aftershocks to reconstruct the productivity pattern first. In this case, the point model can also generate relatively good forecasts because the underestimation in direct aftershocks of the mainshock can be compensated for by the overestimation in high-order aftershocks triggered by aftershocks during the forecasting calibration period. But if , where "tar" and "ref" mean the target and reference models, respectively, and i runs over all the aftershocks following the mainshock.
the mainshock rupture and model parameters are well known, the finite-source model still has the potential for generating improved aftershock forecasts within the first 1-2 days following the mainshock, as the results of theoretical information gain show in Figure 4.

Conclusion
In this study we propose the 3-D finite-source ETAS model, which better fits earthquake data than previous versions of the ETAS model and is able to correct the biased estimates of model parameters in the point source ETAS model. From the productivity density distributions of two large foreshocks and the mainshock, we find that areas of high productivity density form complementary patterns to coseismic slip, and large aftershocks tend to occur near the edges of these high-density areas. Considering that these density patterns exclude higher-order aftershocks and only include the probabilities of events directly triggered by the mainshock, the finite-source model provides a more accurate and reasonable way to compare the aftershock intensity to the coseismic slip of the mainshock. The decay of direct aftershocks in space indicates that the static stress change dominates the occurrence of aftershocks near the mainshock rupture. A comparison in forecast performance among four models indicates that incorporating earthquake focal depths in the ETAS model improves the forecast resolution of aftershocks.