Integrated Investigation on Heterogeneous Lower Crust Rheology in Kyushu and Afterslip Behavior Following the 2016 Mw7.1 Kumamoto Earthquake

The viscoelastic lower crust beneath Kyushu Island, influenced by the volcanic arc, interplays with active crustal faults in this region and helps to shape local tectonics. In this study, we employed a three‐dimensional viscoelastic finite element model to gain insights into the lithospheric rheology and crustal faulting kinematics, through modeling the postseismic deformation processes of the 2016 Mw 7.1 Kumamoto earthquake. Our model reveals a viscosity of 2 × 1020 Pa s for the lower crust and 2 × 1019 Pa s for the upper mantle. A reduced lower crust viscosity of 2 × 1019 Pa s in the volcanic arc area is required for better reproducing the Global Positioning System data. The stress‐driven afterslip decays rapidly over time and is up to 0.3 m within 5 years after the earthquake. We propose additional normal‐component afterslip to better explain the complex postseismic deformation in the near field, which may be due to the interaction between the fault and volcano Aso.


Introduction
The Mw 7.1 Kumamoto earthquake on 16 April 2016 occurred along the Futagawa-Hinagu fault zone (Figure 1a), which is part of the Median Tectonic Line (MTL) in the central Kyushu Island (Toda et al., 2016).The earthquake was preceded by two significant foreshocks (Mw ≥ 6.0) on 14 April 2016 and caused widespread structural damage.Previous coseismic rupture models and focal mechanism of the earthquake showed a dominant rightlateral strike-slip motion with significant normal-slip component (Ekström et al., 2012;Z. He et al., 2020;Himematsu & Furuya, 2016;Milliner et al., 2020).The low-velocity and high-attenuation magma chamber below the Aso caldera play an important role in decelerating and terminating the rupture of this earthquake (Hao et al., 2017;Yagi et al., 2016;Yue et al., 2017).Additionally, the opening of the Okinawa Trough (OT) and the upwelling flow beneath the active volcanoes in the Beppu-Shimabara Graben (BSG) also profoundly influence the stress field and tectonic environment in the region (Takayama & Yoshida, 2007;Yu et al., 2019).No earthquake with a magnitude greater than 7.0 had been recorded on the MTL in Kyushu before the 2016 Kumamoto earthquake (S.Matsumoto et al., 2015), making it a unique opportunity to better understand the lithospheric rheology beneath the island.
Understanding the lithospheric rheology may help us explore the tectonics in situ, for example, the interaction between fault and volcanos (Yue et al., 2017), the generation of the BSG (D. Zhao et al., 2018), and so on.Previous studies on the Kumamoto earthquake mainly focused on its early-period postseismic deformation.Using 8.5-month Global Positioning System (GPS) observations, Pollitz et al. (2017) proposed that the steady-state viscosity of the lower crust in Central Kyushu (∼2 × 10 19 Pa s) is lower than the surrounding area (∼2 × 10 20 Pa s), and suggested a weak upper mantle in fore arc that may relates to the shallow fluid release in North Kyushu.Himematsu and Furuya (2020) analyzed 2-year postseismic deformation with interferometric synthetic aperture radar data, and concluded that the complex near-field deformation requires consideration of multiple postseismic processes, along with potential influences from geological heterogeneous structure.Moore et al. (2017) investigated the transient rheology of the lower crust over 3-month GPS observations, and reported an abnormal lowviscosity area beneath the subduction-related volcanic arc, characterized by a transient viscosity of ∼10 17 Pa s.
In this study, we utilize a three-dimensional (3D) viscoelastic finite element model (FEM) to study the 5-year postseismic deformation (Figure 1 and Figure S1 in Supporting Information S1), considering both the afterslip and viscoelastic relaxation.We first construct test models with uniform lower crust and upper mantle to constrain the first-order viscosities using far-field GPS observations.Then, we consider the rheological heterogeneities in the BSG and volcanic arc region.To better fit the complex near-field deformation, we propose a potential additional afterslip and emphasize the contribution of the Mt.Aso in this process.

GPS Data
We extract initial 5-year postseismic deformation (Text S1 and Figure S2 in Supporting Information S1) at more than 200 GPS stations managed by the Global Navigation Satellite System Earth Observation Network System (GEONET).These stations are located within 250 km to the epicenter of the 2016 Kumamoto mainshock.Note that the postseismic deformation at farther stations is less than 1 mm/year, below the resolution of the GPS data.The obtained horizontal displacements exhibit a four-quadrant distribution, consistent with the deformation characteristics of a typical strike-slip event; Widespread uplift dominates the vertical deformation, except in the near field with much more complex pattern (Figure 1c).Considering the significantly complex vertical deformation pattern and potential effects of the Aso volcano on vertical component (Himematsu & Furuya, 2020;Mochizuki & Mitsui, 2016), we use only the horizontal displacements to constrain the model parameters, while the vertical deformation is presented as well to further verify the model performance.

Finite Element Model
Based on previously published studies on intraplate earthquakes (Chen et al., 2022;Pollitz et al., 2017) and tomography studies of this region (Lallemand, 2016;M. Nakamura & Umedu, 2009;Salah & Zhao, 2004;Yoshioka & Ito, 2001), we have constructed a 3D viscoelastic FEM (Figure 1b, Figure S1 and Text S1 in Supporting Information S1) to study the processes responsible for the postseismic deformation of the 2016 Kumamoto earthquake sequence, that is, the viscoelastic relaxation and afterslip.The coseismic slip model utilized in this study is from Yue et al. (2017), which considered multiple data sets and revealed a detailed slip pattern (Table S1 in Supporting Information S1).
Viscoelastic units of the model include the continental (Eurasian) lower crust, the continental upper mantle, the oceanic (Philippine Sea) asthenosphere and the oceanic upper mantle.To capture both the long-term and rapid short-term response after an earthquake, here we use the bi-viscous Burgers rheology (a Maxwell element in series with a Kelvin element) to represent viscoelastic relaxation (Muller, 1986).Following previous studies (Hu et al., 2016;Zhang et al., 2024), we assume that the transient viscosity (η K ) is one order of magnitude lower than the steady-state viscosity (η M ).Hereafter, we refer to the η M of each unit as their "viscosity." In the FEM, We use a viscoelastic shear zone attached to the main fault to simulate the stress-driven afterslip (Hu et al., 2016).The afterslip amount is calculated by differentiating displacements of the upper and bottom planes of the shear zone.Viscous creep in a shear zone has been demonstrated to generate a surface deformation pattern similar to that due to velocity-strengthening frictional afterslip (Hearn et al., 2002).To balance the computational efficiency and physical representation of fault slip, the thickness of the shear zone is fixed to be 2 km.

Uniform Continental Lower Crust
In this model, we assume that the continental lower crust is uniform in rheology.Because surface deformation on land where GPS stations locate is not sensitive to the stress relaxation in the oceanic units that are far from the rupture area, we assume the same rheological property in the oceanic side as that in Hu et al. (2016), that is, an 80km asthenosphere with a viscosity of 2 × 10 18 Pa s overlying the oceanic upper mantle with a viscosity of 10 20 Pa s (Figure S3d in Supporting Information S1).We thus grid search only the viscosities of the shear zone (η S ), lower crust (η L ) and (continental) upper mantle (η U ).
We have constructed more than 600 test models, each evaluated by comparing their predictions with GPS observations at three time windows: 0-1 year, 2-3 years and 4-5 years after the Kumamoto earthquake.The complex near-field surface deformation is likely dominated by afterslip, concurrently modulated by multiple secondary effects, including volcano activities of Mt.Aso, poroelastic rebound, possible effects from the viscoelastic relaxation of a weaker lower crust in the BSG and volcanic arc region and so on (Hashimoto, 2020;Moore et al., 2017;Pollitz et al., 2017), which may result in substantial fitting residuals of GPS data in near field (<50 km from the rupture fault) for the uniform-lower-crust model.To minimize contamination from these secondary-order mechanisms, we use only the far-field GPS data (>50 km from the rupture fault) to constrain the three model parameters (In contrast, both far-and near-field GPS data are used to constrain the heterogeneouslower-crust model in Section 4.2, which give a better constraint on η S ).The average misfit (Text S1 in Supporting Information S1) is used to evaluate the fit of test models.
Compared with η U that well constrained to be ∼1-4 × 10 19 Pa s with an optimal value 2 × 10 19 Pa s, the η L is less constrained to be 0.5-5 × 10 20 Pa s, with an optimal value 2 × 10 20 Pa s (Figure S4 in Supporting Information S1).The broader range for η L may indicate the existence of heterogeneous structures in the lower crust.These preferred viscosities, featuring same magnitude as that in previous studies (Pollitz et al., 2017), enable our uniform-lower-crust model to reproduce the first-order pattern and overall temporal change of the far-field GPS observations (Figure S5a in Supporting Information S1).However, the uniform model cannot fit the GPS observations in northern near field well and underestimates the deformation around the volcanic arc region in southern Kyushu (Figure S6a in Supporting Information S1).To explore whether the misfit can be alleviated by varied coseismic slip models, we test two additional models with coseismic models modified from Himematsu and Furuya (2016) and H. Kobayashi et al. (2017) (Figure S7a in Supporting Information S1).All models display similar misfits in these regions (Figure S7b in Supporting Information S1), indicating the presence of additional factors influencing near-field deformation.

Heterogeneous Continental Lower Crust
Tomography results show that widely distributed heterogeneous bodies related to active volcanoes and magmatic fluids exist in and around the rupture area of the earthquake (Komatsu et al., 2017;Shito et al., 2017), indicating the necessity to further consider a heterogeneous lower crust.Here we construct a 50-km-wide volcanic arc region based on the quaternary volcanic front associated with the subduction and dehydration of the Philippine Sea plate (Nakajima & Hasegawa, 2007), and a ∼30-km-wide western BSG (WBSG) (Toda et al., 2016) since its eastern portion overlap with the volcanic arc and is likely dominated by the volcanos (D.Zhao et al., 2018) (Figure 1a).This setup allows us to test whether this WBSG exhibits lower viscosity (Pollitz et al., 2017), or not (Moore et al., 2017).The viscosity of the remaining lower crust and upper mantle are fixed at the values obtained in Section 4.1.We then take the same grid search approach as that in the previous section and construct more than 500 test models to refine the viscosities of the WBSG (η B ), volcanic arc (η V ) and shear zone.
In the heterogeneous-lower-crust model, the η V is well constrained to be ∼1-3 × 10 19 Pa s, about one order lower than the η L , consistent with regional heat flow and attenuation anomaly studies (T.Matsumoto, 2007;D. Zhao et al., 2018).In contrast, the η B is loosely constrained with only a lower boundary ∼10 20 Pa s (Figure 2).For convenience, we set η B equal to η L at 2 × 10 20 Pa s.Despite different GPS data set for constraint, the preferred η S in uniform-and heterogeneous-lower-crust models show similar results (10 18 Pa s for the uniform one and 2 × 10 18 Pa s for the heterogeneous one) (Figure 2 and Figure S4 in Supporting Information S1).Their slight discrepancy may be due to the trade-off between the shear zone and volcanic arc.Compared to the model in Section 4.1, the optimal heterogeneous-lower-crust model has a better performance in South Kyushu within the volcanic arc region, abating the underestimation problems in last section (Figures S6a and S6b in Supporting Information S1).

Preferred Model
The preferred model, as identified in Section 4.2, features a weaker lower crust in the volcanic arc region and a "stiff" WBSG, reproducing the first-order pattern of GPS observations in the horizontal components (Figure 3, Figures S5b and S8 in Supporting Information S1).To further quantify the individual contribution of these two rheological units, we vary η V and η B in the preferred model and compare the model performance along an example profile in each region (Figures 3b and 3c).In the volcanic arc area, the model with η V = 2 × 10 19 Pa s better fits the GPS data, especially in the South-North direction (Figure 3b and Figure S6d in Supporting Information S1).For the WBSG, the model is insensitive to η B ≥ 10 20 Pa s, while a lower η B will damage the fit to GPS observations (Figure 3c and Figure S6d in Supporting Information S1).Consequently, the results do not support a weaker WBSG, consistent with tomographic studies that reported no low velocity or high attenuation anomalies in lower crust within non-volcanic areas (Shito et al., 2017;Wang et al., 2017;D. Zhao et al., 2018).
Compared to the horizontal deformation, the preferred model successfully reproduces the uplift in the far field and subsidence in the southern near field but with less precision in magnitude (Figure S8 in Supporting Information S1).The predicted opposite direction at a few stations in the northern and eastern near field may be due to the influence of the complex volcanic activities (Ozawa et al., 2016), such as the processes relate to the eruption of the Mt.Aso 6 months after the mainshock (Sato et al., 2018).

Afterslip Distribution and Evolution
The spatial pattern of postseismic deformation exhibit a wide variety in the near field and are believed to be caused by a combination of multiple processes (Hashimoto, 2020;Himematsu & Furuya, 2020).Here we first assess to what extent the stress-driven afterslip can be accounted for the observed deformation.The preferred model shows a dominant right-lateral afterslip of up to ∼300 mm, most of which concentrates downdip of the rupture zone over the Futagawa fault.This afterslip decays rapidly from up to ∼200 mm in the first year to ∼10 mm in the fifth year (Figure S9a in Supporting Information S1).
The model with stress-driven afterslip roughly reproduces the first-order near-field GPS observations (Figure 3).However, a systematic misfit persists between the predicted and observed displacements at many northern nearfield stations.This misfit persists within 5 years after the earthquake (Figures S5b and S8 in Supporting Information S1) and cannot be explained by a short-lived fluid influence (poroelastic rebound etc.) or volcano activity, for example, the magma chamber inflation beneath Aso volcano (the resulting deformation pattern also deviates from residuals; Text S2 and Figure S10a in Supporting Information S1).Adjusting the viscosity of the shear zone (Figure S3a in Supporting Information S1) or two heterogeneous units in lower crust (Figure S11 in Supporting Information S1) fails to reduce the residuals either.Furthermore, introducing a more heterogeneous lower crust, that is, incorporating a 20-km-wide weakened block with a viscosity of 2 × 10 16 Pa s in region displaying significant near-field deformation residuals and featuring geothermal gradient anomaly (Figure 2 in Pollitz et al. (2017)) also fails to alleviate the residuals (Text S2 and Figure S10b in Supporting Information S1).These findings may suggest a complex, higher-order heterogeneity in the afterslip.

Rheological Heterogeneity of the Continental Lower Crust
In the volcanic arc region, the viscosity of lower crust is ∼2 × 10 19 Pa s, much lower than that of other areas.To explore if this weakened volcanic arc is necessary, we isolate the effects of viscoelastic relaxation in the upper mantle and uniform lower crust from the GPS data, then invert residuals to get a kinematic afterslip (B.Zhao et al., 2022).However, this approach cannot resolve the underestimation in southern Kyushu (Figure S6c in Supporting Information S1), indicating the necessity of including the volcanic arc, whose contribution thus cannot be replaced by the afterslip.This weakened lower crust, coinciding with the low-velocity and high-attenuation anomalies identified in seismic surveys (Sudo & Kong, 2001;Tsutsui & Sudo, 2004;Unglert et al., 2011), is consistent with the region of low transient viscosity reported in Moore et al. (2017).However, the viscosity obtained in this study is about one order of magnitude higher.This discrepancy could be due to two factors: (a) We apply the 5-year GPS observations while Moore et al. (2017) used 3-month GPS data and focused on transient viscosity.The initial rapid decay of the effective viscosity with time may not be reflected by the proportional relation between Maxwell and Kelvin bodies in our model.(b) We calculate η V as an average viscosity over the 50-km-wide volcanic arc region, while Moore et al. (2017) imaged a finer-scale distribution of the viscosity in the lower crust.In the WBSG region, the viscosity of lower crust exceeds 10 20 Pa s, also consistent with Moore et al. ( 2017) that reported no low-viscosity anomaly over here, even beneath the non-subduction-related Unzen volcano (Figure 1a).This is also supported by the tomography studies that show no low-velocity anomaly in the WBSG (Wang et al., 2017;D. Zhao et al., 2018).However, Pollitz et al. (2017) proposed a weak BSG through modeling 8.5-month GPS data.This discrepancy is mostly due to the different model setup.Pollitz et al. (2017) used a simplified model with the whole central Kyushu as the BSG region and revealed the first-order rheology contrast.Their central Kyushu unit may include the collective effects of the volcanic arc and WBSG, and thus lost the resolution on the heterogeneity in the BSG region.The "stiff" WBSG in our model declines the assumption that hot mantle upwelling dominates the generation of BSG (Logatchev et al., 1983).Other factors, such as, northward extension of the opening OT (Nakahigashi et al., 2004) and westward extension of the MTL (D. Zhao et al., 2018), may also play important roles in the formation of the BSG.Pollitz et al. (2017) inverted the residual between GPS observations and the predictions of viscoelastic relaxation models for afterslip.They reported that most afterslip distributed on the eastern Hinagu fault, overlapping the coseismic rupture area and contradicting the stress-driven afterslip pattern.Moreover, they derived the afterslip to be predominantly right-lateral, while significant coseismic normal-component slip was reported over the Futagawa Fault (Himematsu & Furuya, 2020;Yagi et al., 2016;Yue et al., 2017).In addition, a lot of historical earthquakes and aftershocks in the middle Futagawa fault also feature major normal-slip component (Hao et al., 2017;T. Kobayashi, 2017;T. Nakamura & Aoi, 2017).Our preferred model underestimates the deformation in the northeastern near field (Figure 3a and Figure S6b in Supporting Information S1).To address this, we invert the residuals with the constraint of allowing only right-lateral strike slip over the Hinagu fault but both right-lateral and normal slip over the Futagawa Fault (B.Zhao et al., 2022).The inverted results reveal ∼0.44 m normal-component slip on the Futagawa fault (Figure 4a and Figure S9 in Supporting Information S1), which cannot be captured by our stress-driven shear zone afterslip.Incorporating this inverted afterslip substantially decreases the misfit in the near field (Figure 4b,Figures S12 and S13a in Supporting Information S1).This is consistent with Hashimoto (2020) who reported substantial normal slip on the same fault.Even if we did not constrain the location of the afterslip, the inverted additional slip spontaneously fringes the coseismic rupture zone (Figure S9b in Supporting Information S1).The combined afterslip distributes complementary with aftershocks (Figure 4a), which may indicate mixed strain release process over the fault.

Secondary-Order Deformation Processes
Additionally, We explore the alternative approach, eliminating the influence of viscoelastic relaxation (lower crust, upper mantle, and volcanic arc) from the GPS observations and then straightly inverting the residuals for kinematic afterslip.This alternative method yields comparable and only slightly superior model performance to that of the hybrid afterslip model (shear zone afterslip plus inverted additional afterslip) (Figure S13 in Supporting Information S1).Notably, this approach also gets a substantial normal-component afterslip (∼0.35 m) on the Futagawa Fault.However, straightly inverted kinematic afterslip lacks constraints from physical processes such as coseismic stress changes and/or friction conditions on the fault (Freed et al., 2006;K. He et al., 2022).Furthermore, the GPS observations are collective effects of all the tectonic processes.That means, the afterslip inverted from the kinematic approach may actually include contributions from other processes and/or structure heterogeneities, which may bias the interpretation of model results.Therefore, we prefer the approach that combines the shear zone and additional kinematic slip, which not only provides comparable data fitting, but also allows for the identification of the stress-driven afterslip, facilitating subsequent discussions on the additional afterslip mechanisms.One possible explanation for this additional afterslip is that the significant topography contrast near the Aso volcano, approximately 1,000 m, could generate gravitationally induced shear forces leading to prominent normal-faulting stress accumulation (Yue et al., 2017).This shear stress exists persistently over earthquake cycles, thereby promoting additional normal-component afterslip over the Futagawa fault.
Alternatively, ductile material beneath the Aso volcano, characterized by velocity-strengthening behavior, could contribute to the termination of rupture and cause stress concentration at the rupture boundary.Subsequent relaxation of the concentrated stress leads to large localized afterslip (Yue et al., 2017).The observed complementary between the additional slip and coseismic rupture zone (Figure S9b in Supporting Information S1) agrees with the expected behavior caused by such ductile material.

Conclusion
In this study, we have derived the first 5-year postseismic deformation at more than 200 GPS stations following the 2016 Kumamoto earthquake to constrain the rheological properties beneath the Kyushu Island and the afterslip distribution.Model results have determined the viscosities of the lower crust and upper mantle to be 2 × 10 20 Pa s and 2 × 10 19 Pa s, respectively.A weakened lower crust in the volcanic arc is required to better reproduce the GPS observations in the local area and has a viscosity of 2 × 10 19 Pa s.The viscosity of the lower crust in the WBSG is no less than 10 20 Pa s, supporting the muti-generation mechanisms of BSG.An additional normal-component afterslip over the Futagawa fault helps better explain the complex surface deformation pattern in the northern near field, which may indicate mixed stress accumulation/release process over this volcano-fault system.

Figure 1 .
Figure 1.Map of the Kyushu Island and the conceptual finite element model (FEM).(a) Beach balls represent the focal mechanism of the mainshock (red) and two foreshocks (blue).The Coral belt represents the volcanic arc.The blue shaded area represents the western Beppu-Shimabara Graban (WBSG), while the blue dashed lines mark the overlap of Beppu-Shimabara graban and volcanic arc.Black triangles in inset denote the volcanos used to constrain the geometry of the volcanic arc.(b) Conceptual illustration of the FEM.η Μ , η Κ and μ represent the steady-state viscosity, transient viscosity and shear modulus, respectively.WBSG and volcanic arc are same as that in (a).(c) Five-year postseismic Global Positioning System (GPS) observations.Green arrows and colored dots represent horizontal and vertical GPS observations in first 5 years after the earthquake, respectively.The rectangles with black edge correspond to the stations shown in Figure S2 in Supporting Information S1.MTL, Median Tectonic Line; OT, Okinawa Trough; EP, Eurasian Plate; PHS, Philippine Sea Plate; KPR, the Kyushu-Palau Ridge; WPB, Western Philippine Sea Basin; SB, Shikoku Basin.

Figure 2 .
Figure 2. Misfit of test models with heterogeneous lower crust and the trade-off between rheological parameters.(a) Misfit of test models.Each cube represents one test model.(b) Trade-off between η S and η V while η B is fixed at 2 × 10 20 Pa s.The yellow star denotes the optimal model parameters.(c) Trade-off between η B and η S while η V is fixed at 2 × 10 19 Pa s.(d) Trade-off between η B and η V while η S is fixed at 2 × 10 18 Pa s.

Figure 3 .
Figure 3.Comparison of Global Positioning System (GPS) observations with model predictions.(a) Comparison of the 5-year postseismic observations with the predictions of the preferred model.Green and black arrows represent GPS observations and model predictions, respectively.The red and blue lines denote the location of the example profiles in (b) and (c).Coral and blue shadows represent the volcanic arc and western Beppu-Shimabara Graban (WBSG), respectively.(b) Comparison of GPS observations with model predictions along the profile A-B-C whose location is shown as a thick solid red line in (a).Colored solid lines represent the predictions of the test models with different viscosities of the volcanic arc.Solid black dots with error bars represent the east component (upper panel) and north component (bottom panel) of the GPS observations, from stations located within a 30-km width of the profile.(c) Similar to (b) but for the profile D-E whose location is shown as a thick solid blue line in (a).Colored lines represent the predictions of test models with varied viscosity of the WBSG.

Figure 4 .
Figure 4. Performance of the preferred model combined the inverted afterslip.(a) The combination of the stress-driven shear zone afterslip and the inverted additional afterslip within the first 5 years.Blue and green curves represent coseismic slip of the mainshock and two foreshocks, respectively.Gray dots represent relocated aftershocks within 15 days after the earthquake from Yue et al. (2017).(b) Comparison of the 5-year postseismic observations with the predictions combined the preferred model and the inverted afterslip.Green and black arrows represent observations and predictions, respectively.