Temporal gravity anomalies observed in the Tokai area and a possible relationship with slow slips

Abstract

The water in Earth’s mantle is closely related with plate subduction and volcanism. Recent studies revealed that the mantle wedge corner at approximately 30 km depth holds high-pressure water, where many slow earthquakes occur. To quantify how such water behaves during slow earthquakes helps us understand the mechanisms of these earthquakes and (eventually) a part of the long-term water cycle between the interior and surface of the Earth. However, little evidence has thus far been reported on the transient flows of such deep water. Here, we report anomalous, negative mass anomalies during two recent long-term slow slip events in the Tokai area in Japan, which were detected by absolute gravity measurements over 20 years. We present a poroelastic fluid flow model assuming a localized deformation within the fault fracture zone. The model can reproduce the gravity change with a permeability range between those suggested by laboratory experiments and numerical simulations of slow earthquakes.

Introduction

The observation networks of the global navigation satellite systems (GNSSs) have developed during the recent few decades, which has enabled near real-time monitoring of the surface crustal deformation of the solid Earth (Bock and Melgar 2016). The obtained deformation data have been used for studying various tectonic processes. In convergent plate boundaries, where large earthquakes frequently occur, many observational facts have been reported to be associated with co-, post- and inter-seismic crustal deformations. Such geodetic observations provide data on space–time variations in the inter-plate coupling, which are useful to infer frictional properties that govern diverse patterns of earthquake occurrences, including fast and slow ruptures (Scholz 1998; Obara and Kato 2016).

In addition to crustal deformation, co- and post-seismic gravity changes have been detected by terrestrial and satellite measurements (e.g., Tanaka et al. (2001); Tanaka and Heki (2014)). However, only a small number of papers have reported inter-seismic gravity changes based on μGal-level precisions (Mazzotti et al. 2007; Tanaka et al. 2010; Van Camp et al. 2011) (1 μGal = 10⁻⁸ ms⁻²). Observed temporal gravity changes catch underground density redistributions, which can give us information that cannot be obtained from surface displacements. For instance, in volcano observations, we can monitor height changes in magma head (Okubo et al. 2013). The sensitivity for crustal fluids such as magma and water is a distinct feature of gravity measurement. However, most of the above-mentioned studies for co- and post-seismic deformations have interpreted the observed gravity changes as elastic and viscoelastic deformations caused by fault slips, as for the GNSS data.

It is known that the presence of pore fluids supplied from the subducted slab affects earthquake and volcano activities in subduction zones (e.g., Katayama et al. 2012; Kawamoto et al. 2015; Nakajima and Hasegawa 2016). However, there are still a few observations which suggest relationships between earthquake generation and ‘transient’ flows of deep pore fluids, accompanied by temporal variations in fluid pressure (e.g., Sibson 1992; Ogawa and Heki 2007; Husen and Kissling 2001). Elucidating temporal variations in fluid flow will help us to understand the generation mechanisms of long-term slow slip events (LSSEs) and eventually a part of the water budget in subduction zones (Hirschmann and Kohlstedt 2012).

In this paper, we present temporal variations in gravity in the Tokai area located at the eastern margin of the Nankai Trough over approximately 20 years, detected by absolute gravity measurements (“Methods” section). In this area, the presence of high-pressure pore fluids was confirmed by seismological observations (Kato et al. 2010) and LSSEs occurred twice during the observation period. We will show that unexplainable gravity changes remain after removing the contributions from known natural phenomena (“Results” section). We discuss a temporal coincidence between the residual gravity change and the LSSEs and try to interpret it with a simple poroelastic model (“Discussion” section).

Methods

Observation

Figure 1 shows the Tokai area where the Philippine Sea Plate (PH) is subducting beneath the continental plate (AM) toward the northwest. The LSSEs occurred during 2001–2006 and 2013–2017 at an approximately 30 km depth on the plate interface (Ozawa et al. 2016). The gravity measurement sites are located above the locked zone (OMZ), the transition zone where the LSSEs dominantly occurred (TYH), and between these two sites (KKG); at these three sites, campaign measurements with FG5 absolute gravimeters had been carried out once or twice a year since 1996, 2004 and 2008, respectively. The observation instruments and the data processing are described in Okubo et al. (1997). Figure 1c shows the obtained raw data (g₀). For each data point, the measurement was taken for at least 24 h to cancel local, un-modeled tidal effects with the number of the total drops being at least 4800. Consequently, the measurement precision for each data point becomes better than 1 μGal on average, consistent with the nominal precision of the instrument (Okubo et al. 1997).

Fig. 1

The gravity observation sites and the raw data. a The study area. b The locations of the gravity observation sites. OMZ denotes the crustal deformation monitoring station of the GSI, and KKG and TYH are those of Nagoya University. The gray and red contours represent the long-term slow slips which occurred during 2001–2006 and 2013–2017, respectively [Ozawa et al. 2016; the Geospatial Information Authority of Japan (GSI) (2017)]. The blue circle denotes the asperity of the anticipated great earthquake. Along the thick green line, the presence of high-pressure fluids is suggested (Kato et al. 2010). c The observed gravity changes (g₀) obtained after standard corrections of tides, polar motion and the barometric effect. The vertical axes denote relative gravity change in μGal and offsets are added for visualization. The different colors denote different S/Ns of the instruments. The error bars denote measurement errors (1σ)

Data analysis

The apparent gravity changes due to the vertical displacement at the sites were corrected using daily coordinates at the adjacent GNSS stations (F3 solution), the mean free-air gradient (− 3.086 μGal/cm) and a Bouguer plate with an average crustal density (2.67 g/cm³). The ID numbers of the employed GNSS stations are 93101 (before the year 1997) and 960625 (after 1997) for OMZ, 93093 for KKG, and 93104 for TYH; these stations are distant from the sites by 0.9, 7.0, 6.3 and 2.0 km, respectively. The secular subduction rate at station 93101 nearer OMZ was faster by 0.80 mm/year than that at 960625, which was added to the rate for 960625. Because daily coordinates include atmospheric noise, we adopted 1-year moving averages for the correction of the gravity data (red curves in Fig. 2a).

Fig. 2

The corrections of the gravity data. a The raw data (g ₀ ) at OMZ (blue), KKG (brown) and TYH (green) and the apparent effect caused by vertical displacements (g _h ). The gray and red bars denote the periods of the SSEs. b g ₁  = g ₀  − g _h and the effects of crustal dilatation caused by slip deficit and slow slip on the plate interface (g _d ). c g ₂  = g ₁  − g _d and the effect by the post-seismic relaxation of the 2011 Tohoku earthquake (g _p ). d g ₂ ′ = g ₂  − g _p and the estimated groundwater noises. The blue lines represent the gravity change calculated from the soil parameters determined by the inversion (g _w (LSQ)). The red and gray lines, respectively, denote the average of the groundwater noises estimated for all the allowable parameters and the maximum and the minimum of them. e g ₃  = g ₂ ′ − g _w (ave). f g₃′ = g ₂ ′ − g _w (LSQ)

The effects of the density changes within the crust due to the LSSEs and tectonic loading were estimated based on an elasticity dislocation theory of Okubo (1992) (Additional file 1) with the following fault slip models. Annual slip velocities on the plate interface before the year 2009 were taken from Ochi and Kato (2013). Ozawa et al. (2016) obtained a slip distribution during 2013–2016 by fixing the deformation rate during 2008–2011. We assumed that the slip deficit rate during 2007–2009 obtained by Ochi and Kato (2013) remained at the same rate after 2009 and superimposed the LSSE inferred by Ozawa et al. (2016) on it. For the latter LSSE, a constant slip velocity of 2 cm/year on a rectangular fault was assumed to continue until 2017 since the expected gravity change was only 0.02 μGal/year, which does not affect the conclusion of this study.

The 2011 Tohoku earthquake caused post-seismic deformation in the Tokai area (Ozawa et al. 2016). The apparent effect due to the post-seismic vertical displacement was already removed using the observed data by the method mentioned before. We estimated a gravity change caused by the density redistribution due to the post-seismic deformation based on the theory of Tanaka et al. (2014). The details are described in Additional file 1.

For the gravity data obtained after the above corrections, we estimated a groundwater noise with the GWATER1D (Kazama et al. 2012), incorporating the daily weather data of the AMeDAS during the 20-year observation period (Japan Meteorological Agency, http://www.data.jma.go.jp/obd/stats/etrn/index.php). Weather data at Omaezaki were used for sites OMZ and KKG and those at Hamamatsu for TYH (Fig. 1b). A grid search algorithm was applied to determine the saturated permeability (K _s ) and diffusivity (D _s ) of the soil and the admittance (α) between the gravity change and the water content at each site (see Additional file 1 for more details).

Results

Figure 2a shows the effect of vertical motion. We see that plate subduction increases gravity at OMZ by approximately 15 μGals in 20 years, which is approximately the same as the observed change (g₀) in 20 years. During 2002–2006, when the former LSSE occurred, the apparent effect is temporarily constant. At TYH, plate subduction causes uplift, which decreases gravity. The decreasing rate accelerates during the period of the former LSSE and after the 2011 Tohoku earthquake. The total amount of the correction is – 20 μGals from the year 2005 to 2017, which is twice as large as the change in g₀. At KKG, the vertical motion is relatively small and so is the apparent effect. However, the sign of the correction is opposite. Figure 2b shows the gravity changes obtained by removing those apparent effects. The change at OMZ indicates a small positive trend after year 2000 and transient variations with amplitudes of ~ 5 μGals. At KKG and TYH, positive trends with 0.7 μGal/year are visible, respectively.

The red curves in Fig. 2b shows the contributions from crustal dilatation/compression due to slip on the plate interface (g _d ). At OMZ and KKG close to the Suruga Trough, gravity increases by 3 μGals in the periods after year 2000 and 2008, which can explain 30 and 10% of the secular increase in the gravity change, respectively. The correction at TYH is only 1 μGal from 2005 to 2017. Figure 2c shows the corrected result. The red curves in this figure illustrate the effect due to the co- and post-seismic density changes caused by the 2011 Tohoku earthquake (g _p ). At all sites, calculated gravity increases by approximately 3 μGals from March 2011 to 2017. Figure 2d shows the results with the effect being removed. We see that linear trends over the whole observation periods for g₂′ decrease at all sites than for g₁ in Fig. 2b. However, there still remains a large fluctuation at OMZ and a significant positive trend at KKG and TYH.

The blue lines in Fig. 2d display the modeled groundwater effect. The deviations of the effect are approximately 2–3 μGals. At KKG, a slower fluctuation with periods of ~ 3 years is also superimposed. At OMZ and TYH, we find no such longer-term variations. The result of the inversion implies that the large fluctuation at OMZ and the linear trend at KKG cannot be reproduced by the groundwater effect. Table 1 shows the inferred soil parameters and the admittances. The admittance at OMZ is negative due to the presence of a hill just behind the observation station, where accumulated groundwater attracts the instrument upward more strongly than from the area below the instrument. The admittances are positive at KKG and TYH where the sites are located in flatter areas. Figure 2e, f shows the corrected results obtained from the different groundwater contributions [g _w (ave) and g _w (LSQ), respectively] (Additional file 1). The temporal behaviors of the gravity changes in Fig. 2d–f are basically the same, indicating that there still remain changes which cannot be reproduced by the groundwater model.

Table 1 The soil parameters estimated by the inversion

Discussion

Rate changes in the time-series gravity data and the LSSEs

To identify temporal variations in the gravity changes associated with the LSSEs, we calculated rate changes for the gravity data at all the correction levels (Fig. 3a–c). We see that the rate changes for g₁, g₂, g₂′ and g₃′ show almost similar tendencies. The error bars denote the standard error of the linear regression. We estimated the error of the correction for g _h by assuming that uncertainties for the average crustal density are 10%. Moreover, we calculated the error for g _d by considering that the average uncertainty for the slip velocity distribution of Ochi and Kato (2013) is 10% and that of Ozawa et al. (2016) is 30%. However, these errors are smaller than 0.1 μGal/year, and after all, the total uncertainties are dominated by the deviations of the observed data and the correction levels.

Fig. 3

The rate changes for gravity and model comparisons. a–c Linear trends for the gravity changes obtained by a 4-year moving average. The gray and red colors denote the periods of the LSSEs. The orange circles show the p value for the trend of g ₃ ′. d–f A comparison of the poroelastic models A-C with the gravity change after all the corrections (g₃′)

To confirm a statistical significance of the rates, we plotted the p value for g₃′ where all the corrections were applied (Helsel et al. 2006). At OMZ, gravity significantly decreases before the occurrence of the former LSSE (1997–2001) (p < 0.05). After 2001, significant, positive linear trends appear only during 2007–2012, implying that gravity tends to increase during the periods excluding the LSSEs. A similar feature is seen for KKG, although determination errors of the trends are larger due to the fewer number of the data than for OMZ. At TYH, no significant trends associated with the LSSEs are seen.

The above result indicates a weak temporal correlation between the (relative) decrease in linear trends and the periods with the LSSEs at the shallower part of the subduction zone. A decrease in gravity coinciding with tremor episodes was noted at a site spatially equivalent to KKG in the Cascadia subduction zone although the time scale of the gravity change was much shorter than for long-term SSEs (Fig. 6 of Lambert et al. 2006). The negative trend before year 2000 at OMZ might be caused by unknown, local effects irrelevant to the LSSEs. In addition, we could not identify the reason for the secular increases with 0.5–0.9 μGal/year after year 2000 arising at all the sites.

An interpretation of the transient mass anomaly

The above result implies that the LSSEs cause inelastic deformation, as the elastic deformation model could not quantitatively explain the changes in the trends (g _d in Fig. 2b). In the following, we discuss the possibility that poroelastic deformation occurred associated with high-pressure pore fluids, as a candidate of inelastic deformation.

First, we roughly estimate the amount of the observed mass change in terms of a Bouguer plate. Let us assume that the decrease in gravity during the LSSEs is ~ 4 μGals (Fig. 2f). When a plate has the density of water (1 g/cm³), its thickness to explain the change is 0.1 m. The same change is generated when: (1) a plate with a thickness of 1 km and a density of 1 g/cm³ is compressed by 10⁻⁴; or (2) a plate with a thickness of 30 km and a density of 3 g/cm³ is compressed by 10⁻⁶. Case (1) assumes that deformation occurs only within the 1-km thick fracture zone. For this case, Suzuki and Yamashita (2006) (hereafter abbreviated by SY06) gave a theoretical framework to describe the interactions between slow rupture and fluid flow. Case (2) assumes that the bulk outside the fracture zone also deforms. In this case, the surface vertical displacement occurs, by which the slip distribution on the plate interface and the correction for g_h in Fig. 2a must be reevaluated based on a poroelastic model. In this study, we consider case (1), which can be treated more easily, using the framework of SY06.

The strain change by 10⁻⁴ is approximately two orders of magnitude larger than the elastic strain caused by the slow slip (slip/spatial scale ~ 10 cm/100 km). This means that in order to explain the observed mass change, a stronger driving source to cause poroelastic deformation is necessary than the slow slip (the term including $\varepsilon_{s}$ in SY06). Based on this consideration, we constructed a phenomenological model based on fault-valve action (Sibson 1992). According to this mechanism, permeability temporarily increases with a fault rupture and then recovers by self-sealing. We assume that this holds for LSSEs and identified allowable permeability values from the gravity data.

The one-dimensional governing equation for fluid pressure increment $p_{f}$ is given from Eqs. (11) and (30) in SY06 as

$$ C\frac{{\partial p_{f} }}{\partial t} = \frac{d\kappa }{dx}\frac{{\partial p_{f} }}{\partial x} + \kappa \frac{{\partial^{2} p_{f} }}{{\partial x^{2} }} $$

(1)

$$ C = \eta \left( {\frac{\phi }{{K_{f} }} + \frac{b - \phi }{{K_{s} }}} \right) $$

where the terms associated with temperature and strain changes are neglected, x denotes the up-dip direction (Fig. 3d), and b, K _f , K _s , ϕ, η and $\kappa$ represent the Biot–Willis parameter, bulk moduli for the fluid and solid phases, porosity, and the viscosity of fluid and permeability, respectively. K _f , K _s and η were set to the same values as in SY06 (30 GPa, 3.3 GPa and 2.82 × 10⁻⁴ Pa s, respectively), and ϕ = 0.03. Moreover, b was given by Eq. (38) in SY06, assuming Poisson media. In contrast to SY06, we introduced a spatial gradient for permeability $\kappa$ so that we can obtain a steady-state analytical solution for Eq. (1) by assuming

$$ \kappa \left( x \right) = \kappa_{0} \exp \left( {ax} \right),\;a = \left( {1/L} \right)\ln \left( {\kappa_{L} /\kappa_{0} } \right) $$

(2)

where x = 0 and L corresponds to the depths of 30 and 15 km, respectively, and $\kappa_{0}$ and $\kappa_{L}$ denote the permeability values at x = 0 and L, respectively. Such an exponential decay of permeability with depth is suggested by experiments (David et al. 1994) (note that $\kappa_{L} > \kappa_{0}$). The positions of x = 0 and L were chosen so that an expected gravity change becomes the largest near KKG (Fig. 4a, b). The boundary conditions are $p_{f} = 600$ and 0 MPa at x = 0 and L, respectively, and $p_{f} = 0$ represents hydrostatic pressure. The steady-state solution is then given as

$$ P_{f} \left( x \right) = \frac{{P_{f} \left( {x = 0} \right)}}{{(\kappa_{L} /\kappa_{0} ) - 1}}\left[ {\left( {\frac{{\kappa_{L} }}{{\kappa_{0} }}} \right)\exp \left( { - \,ax} \right) - 1} \right] $$

(3)

(the blue curve in Fig. 4d). The bulk density was obtained from Eqs. (14), (19) and (20) in SY06 as

$$ \frac{{{\text{d}}\rho_{B} }}{{{\text{d}}t}} = \left( {\phi \frac{{\rho_{f} }}{{K_{f} }} + \left( {1 - \phi } \right)\frac{{\rho_{s} }}{{K_{s} }}} \right)\frac{{{\text{d}}p_{f} }}{{{\text{d}}t}} $$

(4)

where $\rho_{f}$ and $\rho_{s}$ denote the density for the fluid and solid phases, respectively, and the gravity change at the surface was calculated with Newton’s law.

We assumed that permeability in the shallower part along the x-axis increases from $\kappa_{L}$ to $\kappa_{L}^{{\prime }}$ instantaneously by an order of magnitude at the onset of the LSSEs and decreases to $\kappa_{L}$ at the end of those events so that upward flow occurs. By this assumption, fluid pressure decreases, and we could explain the observed negative gravity change.

Fig. 4

The poroelastic model. a The rate for gravity change at 3 years after the onset of the LSSE for model B. b The same as in (a) but for the pressure increment. c A model configuration. The dip angle is assumed to be constant with depth. d The pressure increment along the X-axis for model B. The blue and red curves represent the steady-state solution (the left vertical axis) and the pressure change in 3 years from the onset of the LSSE (the right vertical axis), respectively

Figure 3d–f compares the observed gravity change (g₃′) with models for different permeability values: $(\kappa_{0} ,\kappa_{L} ,\kappa_{L}^{{\prime }} )$ = (3 × 10⁻¹⁸, 3 × 10⁻¹⁷, 3 × 10⁻¹⁶), (3 × 10⁻¹⁷, 3 × 10⁻¹⁶, 3 × 10⁻¹⁵) and (10⁻¹⁷, 10⁻¹⁶, 10⁻¹⁵) m² for A, B and C, respectively. We see that model B can best explain the changes in the trends before and after the LSSEs. The difference in reproducibility between models A and B can be confirmed by continuing observations. Permeability lower than 10⁻¹⁸ or higher than 10⁻¹⁵ causes too small or large differences in the average rates. In this comparison, linear trends obtained from the observed data are added on the modeled time series (0.5, 0.5 and 0.9 μGal/year for OMZ, TYH and KKG, respectively), because none of the models could generate long-term positive trends (Figure S2). The spatial distribution of gravity change and the pressure increment for model B are shown in Fig. 4a, b, d, as an example.

From the above comparison, we see that gravity data favors permeability of ~ 10⁻¹⁸–10⁻¹⁵ if the model assumption is correct. This range for permeability is between the values inferred from laboratory experiments (Katayama et al. 2012) and those suggested by numerical simulations, above which slow ruptures could occur (Yamashita 2013).

Conclusions

We reported the results of the campaign of absolute gravity measurements at three sites in the Tokai area for the most recent 20 years. The observed gravity change was corrected for the vertical displacement, crustal dilatation, post-seismic relaxation of the 2011 Tohoku earthquake, and the groundwater noise. Significant mass anomalies remained, which could not be interpreted with these effects. The residual gravity changes at the sites OMZ and KKG closer to the Suruga Trough showed transient decreases during the LSSEs (2001–2006 and 2013–2017), whereas no remarkable change was seen above the transition zone (TYH). In addition, secular trends with approximately − 0.2, 0.9 and 0.5 μGal/year were detected at OMZ, KKG and TYH, respectively, for the whole observation periods. A similar unexplained positive long-term gravity trend of ~ 0.5 μGal/year is seen also in the Cascadia subduction zone (Mazzotti et al. 2007).

To interpret the anomalies, we constructed a poroelastic model based on SY06, assuming that the deformation occurs only within the fault fracture zone. We found that the model could reproduce the observed transient decreases, provided that the initial permeability having a spatial gradient suddenly increases during the LSSEs to trigger upward pore-fluid flow. Our model preferred permeability values of 10⁻¹⁸–10⁻¹⁵ m² for the depths between 30 and 15 km. This range was between the values inferred by laboratory experiments and those suggested by numerical simulations. Our model could not explain the observed positive long-term linear trends, and the reason for this is still unknown. It is important to continue the observations to identify whether these trends are maintained or are a part of longer-term fluctuations.