|Home | About | Journals | Submit | Contact Us | Français|
Large earthquakes often trigger viscoelastic adjustment for years to decades depending on the rheological properties and the nature and spatial extent of coseismic stress. The 2006 Mw8.3 thrust and 2007 Mw8.1 normal fault earthquakes of the central Kuril Islands resulted in significant postseismic gravity change in GRACE but without a discernible coseismic gravity change. The gravity increase of ~4 µGal, observed consistently from various GRACE solutions around the epicentral area during 2007–2015, is interpreted as resulting from gradual seafloor uplift by ~6 cm produced by postseismic relaxation. The GRACE data are best fit with a model of 25–35 km for the elastic thickness and ~1018 Pa s for the Maxwell viscosity of the asthenosphere. The large measurable postseismic gravity change (greater than coseismic change) emphasizes the importance of viscoelastic relaxation in understanding tectonic deformation and fault-locking scenarios in the Kuril subduction zone.
On 15 November 2006, the shallow part of the central Kuril Islands megathrust (depths <25 km) failed in a great Mw8.3 earthquake following a swarm of thrust-faulting foreshocks in late September 2006 [Ammon et al., 2008; Steblov et al, 2008; Lay et al., 2009]. Following this event, distributed thrust-faulting aftershocks occurred along the forearc; additionally, a parallel band of primarily normal faulting events occurred in the outer rise region. Based on the predicted evolution of stress redistribution following the 2006 thrust event, Ogata and Toda  hypothesized that continuous triggered slip occurred on the future fault plane of the January 2007 event. On 13 January 2007, the great Mw8.1 normal faulting earthquake ruptured a steeply dipping fault (depths of 4–33 km) beneath the trench slope [Lay et al., 2009]. Various studies of seismic, geodetic, and tsunami data analysis report moment release between 1.6–5.3×1021 N m for the 2006 event and 1.4–2.8×1021 N m for the 2007 event [Steblov et al., 2008; Lay et al., 2009].
Earthquakes with moment magnitudes of Mw8.1–8.3 are nearly at, or below, the detection threshold by Gravity Recovery And Climate Experiment (GRACE). Furthermore, the coseismic gravity changes were dominantly a negative anomaly (due to crustal dilatation) following the 2006 thrust rupture that was then disrupted by a positive anomaly (crustal compression) after the 2007 normal faulting. We calculated a predicted total coseismic gravity change of <1 µGal at GRACE’s spatial resolution. This was certainly below the detection threshold [e.g., Han et al., 2013]. However, following these events, the GRACE Release-05 Level-2 (L2) solutions [Bettadpur et al., 2014] exhibited an increasing gravity trend localized around the epicenters of the 2006 and 2007 earthquakes. This positive trend is consistently observed in various other GRACE solutions processed by different groups (Supporting information). It can be conveniently illustrated by using the CNES/GRGS web tools (http://www.thegraceplotter.com).
Compared to the postseismic gravity change after the 2004 Sumatra-Andaman, 2011 Tohoku-Oki and 2012 Wharton Basin earthquakes [e.g., Han et al., 2008, 2014, 2015], gravity change associated with the 2006–2007 Kuril events is distinctly different; namely, the magnitude of postseismic gravity change becomes by far larger than the coseismic change within a few years. We seek to explain such observations by testing alternate viscoelastic models of the asthenosphere beneath an elastic lithosphere. In this study, we address (i) how the 2006 Mw8.3 thrust and 2007 Mw8.1 normal faulting earthquakes interplay to produce the observed gravity change, (ii) how the rheology of the asthenosphere plays a role in producing postseismic gravity changes, and (iii) how GRACE constrains the rheological structure in the central Kuril Islands. We will present the details of GRACE observations and viscoelastic modeling results with a range of rheological parameters (thickness of lithosphere and asthenosphere viscosity) and of coseismic moment magnitude ratios, and compare with the previous Kuril GPS studies of Kogan et al. [2011, 2013].
We used monthly GRACE gravity data from CSR L2 solutions presenting an effective resolution of 500 km (maximum spherical harmonic degree of 40). Although the GRACE RL05 L2 data are provided up to the maximum degree of 60, the resolution to ~300 km is applicable only at the polar region. The monthly gridded data, computed from the truncated spherical harmonic synthesis, were analyzed by applying the least squares fit of the parameters including annual and semi-annual sinusoids, bias, and a Heaviside step and linear trend after year 2007. Figure 1a presents the linear trend estimated between January 2007 and April 2015. The coseismic and postseismic effects after 2011 Tohoku-Oki earthquake dominate the trend map (Figure 1a). However, a statistically significant positive anomaly of 0.5 µGal/year is found around the epicenters of the Kuril earthquakes.
We examined the finite fault models by Ji [2006, 2007] with the southeast dipping plane for the 2007 normal fault to compute the gravity effect as shown in Figures 1b–d. The synthetic calculation of coseismic (elastic) deformation after the 2006 earthquake results in the negative gravity change northwest of the trench due to crustal dilatation while a positive gravity (due to the hanging wall uplift) is compensated for by the ocean [e.g., Broerse et al., 2011] (Figure 1b). This is a typical pattern in coseismic gravity change after thrust earthquakes [Han et al., 2013]. Broerse et al.  found the horizontal motion of the bathymetry is important to understand the cosesimic gravity change. However, the horizontal postseismic motions during the observation period are ~2 orders of magnitude less than the coseismic horizontal motions, and hence they contribute a much smaller fraction of the net gravity change. The 2007 normal fault yielded instead positive gravity due to compression pushed by the downward motion of the hanging wall (Figure 1c). Again, the gravity change due to density change is greater than the gravity change associated with the elastic vertical motion at this spatial scale. The gravimetric signals produced by the two earthquakes compensate each other and resulted in a gravity anomaly smaller than each of the gravimetric signals in the individual earthquakes (Figure 1d). The elastic gravity change is within ±0.5 µGal and appears away from the epicenters at the spatial scale of 500 km. The solutions by Lay et al.  including the southeast and northwest dipping planes for the 2007 rupture produced slightly different synthetics, but the above conclusion does not change (Figure S1). Elastic deformation alone does not replicate the peaked positive trend near the epicenter in the GRACE data.
The time-series of the gravity change around the peak GRACE signal in Figure 1a (averaged within a spherical cap with a radius of 2.5° centered at 46.5°N, 154.5°E) is shown in Figure 2. The GRACE data (Figure 2, black dotted) presents scatter of 1.8 µGal (root mean square) with respect to the fit (green solid) as described above. The linear trend and annual sinusoid components are statistically significant while the step function and semi-annual components are not. The trend prior to the 2006 earthquake is not statistically different from zero (based on the t-test). The background gravity variation induced by other geophysical processes such as non-tidal atmosphere (6-hourly ECMWF operational data) and ocean mass (OMCT) [Dobslaw et al., 2013] is also shown (magenta solid), which was removed a priori from GRACE gravity solutions. The amplitude of seasonal atmospheric and oceanic variations is estimated to be ~1 µGal. The GRACE data presents the un-/mis-modeled seasonal change with a amplitude of ~0.5 µGal (green solid). Ignoring the two-month time difference between the 2006 and 2007 ruptures, the linear trend after year 2007 is estimated to be 0.45±0.08 µGal/year. The postseismic gravity trend estimate was consistent with the result based on the other GRACE solutions by CNES/GRGS, processed with a different inversion method and different atmosphere and ocean models (Figure S2) and by other groups (Figures S3 and S4).
The gravity increase by ~4 µGal over 8 years from GRACE implies ~6 cm of the seafloor uplift around the area of ~5002 km2 near the epicenters (0.7 µGal of gravity change per 1 cm of seafloor uplift - the Bouguer correction with the density contrast between the crust and ocean). Kogan et al. [2011 and 2013] presented a postseismic uplift rate up to 1.2 cm/year (KURILNET site MATC) from near-field GPS stations along the volcanic arc, which confirmed the postseismic deformation is caused predominantly by viscoelastic relaxation (rather than afterslip). We henceforth assume that the GRACE gravity changes are shaped primarily by viscoelastic relaxation and ignore the possible contribution of afterslip.
For the Kuril earthquake sequence, the postseismic gravity change has emerged as the GRACE time series have lengthened even though the measured, combined coseismic gravity change was insignificant. In this Section, we examined in what circumstances the postseismic gravity changes could be larger than the cosesimic gravity change within a few years after the rupture by quantifying the gravity response functions (Green’s functions) of the viscoelastic Earth. We delineate the overall range of viscoelastic parameters important to control the postseismic gravity change by using a simplified point earthquake source and testing various synthetic structures.
Following Han et al. , the gravity change by an earthquake point source is
where a dislocation source is assumed to be at the north pole and described with the moment tensor components, Mrr, Mrθ, Mrϕ, Mθθ − Mϕϕ, and, Mθϕ. The respective Green’s function is gMrr, gMrθ, gMrϕ, gMθθ−Mϕϕ, and gMθϕ. They are dependent on the Earth’s viscoelastic structure and the depth of dislocation source. The first one gMrr is constructed with the m=0 (monopole source) eigenfunctions of viscoelastic normal modes [Pollitz, 1997; Pollitz et al., 2006]. The pair of and is determined by the m=1 (dipole source) eigenfunctions. The other pair gMθθ−Mϕϕ and gMθϕ is computed by the m=2 (quadrupole source) eigenfunctions. Therefore, there are in principle only three independent functions that govern temporal changes of gravity and displacement fields for each source depth with a given viscoelastic structure.
Using the normal mode approach, we computed the viscoelastic (postseismic) and elastic (coseismic) Green’s functions. The seismic source was assumed to be at 20 km depth within the elastic layer. The viscosity for the upper mantle (220–670 km) and for the lower mantle (670–2891) was fixed to be 1020 Pa s and 1021 Pa s, respectively. The most sensitive parameters to govern temporal evolution of gravity are the asthenosphere viscosity and the lithosphere (elastic) thickness. Three different elastic thicknesses, 25, 40, and 115 km, were tested; below these depths the viscoelastic asthenosphere extends down to 220 km. Two different Maxwell viscosities, 1018 and 1019 Pa s, for the asthenosphere were tested (Table 1 of Supporting Information).
Figure 3 shows the ratio of the postseismic gravity response 5 years after the rupture relative to the coseismic change at each spherical harmonic degree (l) for the m=0, m=1, and m=2 responses, computed from six different cases for the viscoelastic structure. Note that such ratio is not sensitive to the seismic source parameters (except depth), but only to the Earth structure. Generally speaking, the shallower asthenosphere and the lower viscosity yield larger changes in the ratios. The ratio of the m=1 response is significantly smaller than the other response functions for all cases, which implies that the viscoelastic relaxation is governed more effectively by m=0 and m=2 patterns. We found that the asthenosphere viscosity of 1019 Pa s does not produce the postseismic gravity change that is significantly larger than the coseismic change regardless of the thickness of the overlying elastic layer. However, the asthenosphere viscosity of 1018 Pa s and an elastic layer thickness of <40 km can produce an order-of-magnitude larger postseismic gravity change than the coseismic gravity in 5 years. We extend the viscoelastic modeling to constrain such parameters in the following Section.
We exploit the GRACE gravity observation to deduce the optimal parameters that may represent the rheological structure in the central Kuril Islands region. First, it is necessary to fix the seismic sources of the 2006 and 2007 earthquakes. Lay et al.  indicates a large variability in seismic and geodetic source models (point sources as well as finite faults shown in their Tables 1 and 2). The moment release from the 2006 rupture has been estimated to be M0=1.6–5.3×1021 N m including the Global Centroid Moment Tensor (GCMT) M0=3.5×1021 N m [Lay et al., 2009]. For the 2007 normal fault, it is M0=1.4–2.8×1021 N m with GCMT M0=1.8×1021 N m. The moment release ratio of 2006 to 2007 ruptures varies from 1.1 to 3.8. Considering such large uncertainty in the moment release from the two events, we decided to adjust the moment release ratio in addition to the viscoelastic structure parameters. Due to limited detection threshold of the GRACE measurements, it is not possible to constrain the moment release for each earthquake separately.
The synthetic gravity change was computed using the finite fault models of Ji [2006, 2007] and the software VISCO1D for viscoelastic deformation [Pollitz, 1997]. Figure 4 shows the average rate of postseismic gravity change over 2007–2015 for the exemplary case of an elastic thickness of 32 km and the asthenosphere Maxwell viscosity of 1018 Pa s; the deeper (> 220 km) upper mantle and lower mantle are as described in Section 3. The 2006 fault yields the positive gravity change primarily due to vertical motion of various density interfaces (such as seafloor and Moho) because viscoelastic deformation is shown to be nearly incompressible [e.g., Han et al., 2014, 2015]. The top of the asthenosphere is subject to the boundary condition of broad-scale shortening, which triggers postseismic uplift [Melosh, 1983]. In contrast, the 2007 normal fault (for both northwest and southeast dipping planes) yields negative gravity change associated with postseismic subsidence under the top boundary of the asthenosphere subject to extension (See the schematic diagrams in Figure 4). The actual 2007 rupture, along with the pre-slip on this normal faulting rupture plane, in the outer rise caused contrasting vertical deformation to that triggered by the 2006 earthquake at the plate interface.
We write the total postseismic gravity change, that is comparable to GRACE observation, from a given viscoelastic model as follows:
where the synthetic gravity changes g2006 and g2007 are evaluated at time t with the ν viscosity and the elastic thickness da from the 2006 and 2007 ruptures, respectively. The unknown scale factors, α and β, are introduced to accommodate the uncertainty in the moment magnitude in each seismic model, respectively. They are not constrained independently but in terms of a ratio where and is the scalar moment for the 2006 and 2007 earthquake determined from a priori seismic models, e.g., Ji [2006, 2007], respectively. The ratio implies the relative size of the moment release of the 2006 to the 2007 rupture. We tested these three parameters within a reasonable range; i.e., da =26–88 km, ν=1017–1019 Pa s, and γ=1.0–4.0.
We computed numerous cases of synthetic gravity changes and compared with the GRACE data. We quantified the goodness of the fit by computing the R2 statistics (i.e., variance reduction). Figures 5a–c illustrate the R2 statistics with variable viscosities and moment ratios, respectively, for da = 26, 32, and 38 km. The results of da>40 km were ruled out since the fits were significantly degraded within the specified range of ν and γ. The best fits are found consistently when the viscosity is close to 1018 Pa s. The larger elastic layer thickness da requires the larger moment ratio γ (The identical results can be drawn based on the reduced chi-squared χ2 statistics). The long-period seismic estimates of the relative moment ratio between the 2006 and 2007 ruptures are within 2–3 [T. Lay, 2015, personal communication]. With this seismic constraint, the case of da=38 km can also be ruled out. The GRACE data are best fit with a model of 25–35 km for the elastic thickness and ~1018 Pa s for the Maxwell viscosity of the asthenosphere.
The gravity time-series from three best-fit model parameters (color crosses in Figure 5a–c) were compared with the GRACE time series (Figure 5d). All cases reasonably fit the data yielding R2 >0.75 and χ2 ~1.0. We also tested the biviscous (Burgers) rheology that implements the transient (Kelvin) element in series with the steady-state (Maxwell) element [e.g., Han et al., 2014, 2015]. The identical shear modulus was used for the Kelvin and Maxwell elements. We found that the transient viscosity is not as important as the steady-state viscosity to explain the GRACE data in 2007–2015. Any biviscous model with a Maxwell viscosity of 1018 Pa s and a Kelvin viscosity of 1–10×1017 Pa s fits as well as the Maxwell rheology model with 1018 Pa s (Figure 5d). The Burgers or Maxwell model with the steady-state viscosity of 1019 Pa s, however, significantly underpredicts the gravity change (<1 µGal in 8 years) and does not agree with the GRACE measurement (~0.5 µGal/year). The spatial pattern of the average rate of the synthetic gravity changes in 2007–2015 is shown in Figure 5e using the model parameters of da=32 km, ν=1018 Pa s, and γ=3.0. The primary positive gravity around the epicenters is consistent with the GRACE data in terms of the spatial location, the anomaly magnitude, and the temporal pattern, while the secondary smaller anomaly of the negative change is less discernable.
Based on the analysis of 8 years of continuous gravity data, we estimated a Maxwell steady-state viscosity of 1018 Pa s for the asthenosphere, agreeing with Kogan et al. . Based on their analysis of 4 years postseismic GPS displacement data, Kogan et al.  conjectured the Maxwell asthenospheric viscosity is lower than that in Chile and Alaska (1019 Pa s), but these events were much larger in terms of spatial extent and moment magnitude. With the additional seismic constraints on the relative moment between two ruptures, our gravity analysis indicates the elastic layer thickness confined to be <35 km, which may be largely influenced by the mantle wedge beneath the active Kuril volcanic arc. Our preferred elastic thickness is smaller than 62 km fixed (not adjusted) by Kogan et al. [2011, 2013] for their GPS analysis, and rather consistent with the oceanic lithosphere thickness (30–40 km) from the flexural analysis of the bathymetry around the Kuril trench [Hanks, 1971] and with the numerical modeling beneath adjacent northeastern Japan [Muto et al., 2013 and Hu et al., 2014].
Lateral variations in rheology are likely around subduction zones, as differences in elastic plate thickness, elastic moduli and viscosity structure are expected between the landward and oceanic side of the subducting slab (e.g., Billen and Gurnis, ; Abers et al., ). This variation and the presence of the slab can substantially modify predicted postseismic deformation (e.g., Cohen et al., ; Pollitz et al., ; Hu and Wang, ) and presumably the gravity signal as well.
In our modeling, we neglected the lateral heterogeneity of the Earth’s structure that may bias our estimates of the elastic layer thickness. For example, the 2007 outer rise fault occurred in a thicker elastic layer [Lay et al., 2009] and is expected to have produced smaller viscoelastic relaxation than in a thinner elastic layer, which would lead to less disruption of the viscoelastic uplift produced by the 2006 rupture. However, the Maxwell asthenosphere viscosity estimate is consistent regardless of the assigned thickness of the overlying elastic lithosphere. As depicted in Figure 4, the 2006 rupture triggers an initially large viscoelastic uplift while the 2007 rupture adds the extensional boundary condition that pulls down the triggered uplift. What GRACE observed is the remnant of those contrasting postseismic vertical motions.
As the gravity time-series lengthen, the sensitivity increases more drastically for detecting the postseismic trend than the coseismic change. Due to relocking of the shallow segment of the plate interface, a small part (as yet not discernible) of the postseismic gravity change will be constant interseismic deformation as well. While the magnitude of coseismic gravity change is proportional to moment magnitude, postseismic (viscoelastic) gravity change is also dependent on the rheological structures that are variable from region to region. Our analysis of the Mw8.3 2006 and Mw8.1 2007 Kuril earthquakes prompts reexamination of the GRACE data particularly with regard to postseismic gravity changes associated with even smaller earthquakes including the Mw8.1 2009 Samoa-Tonga earthquake doublet. This will become increasingly relevant as the GRACE time-series lengthens and the time series of gravity is extended with GRACE Follow-On starting in 2017 [Flechtner, et al., 2014].
This 2006–2007 Kuril sequence is one of five earthquake sequences that has produced persistent viscoelastic relaxation signals in GRACE gravity; the other cases are those following the 2004, 2010, 2011, and 2012 great earthquakes [e.g., Han et al., 2014, 2015, Broerse et al., 2015]. The prevalence of viscoelastic relaxation was also observed in the seafloor data above the 2011 Tohoku-Oki rupture [Sun et al., 2014]. The continuing analysis of the Earth’s rheological structure is imperative to understanding of the complete earthquake cycle processes. Wang et al. 2013] emphasized the importance of viscoelastic relaxation in drawing an adequate picture of fault-locking distributions at megathrust boundaries, noting that incorrect assessment of seismic and tsunami hazards may result from not considering the Earth’s rheological response. In addition, understanding of the Earth’s rheological response is also crucial for estimating climate signals (water and ice mass change) as the elastic and viscoelastic (un-)loading effects that exist in geodetic data such as GPS and GRACE need to be removed [Ivins et al., 2011].
This work was supported partly by NASA’s GRACE project and Earth Surface and Interior program. We thank Riccardo Riva for the computer codes of the normal mode analysis used in Figure 3. We thank DLR for the GRACE telemetry data and JPL, CSR, GFZ, and GRGS for the high-quality Level-1B and Level-2 products. The GRACE data for this paper are available at podaac.jpl.nasa.gov/GRACE and www.thegraceplotter.com. We thank Jeanne Hardebeck and Ruth Harris for their comments on a preliminary draft. We thank Taco Broerse, an anonymous reviewer, and the editor for their constructive reviews.