Postseismic gravity changes after the 2011 Tohoku earthquake observed by superconducting gravimeters at Mizusawa, Japan

Abstract

Postseismic gravity changes after the 2011 Tohoku earthquake (Mw9.0) were investigated using the data from superconducting gravimeters (SGs) at Mizusawa, Japan. The data in the period from 2014 to 2021 were used in the analysis. The SG data were first corrected for instrumental drift using the results of absolute gravity measurements. Then, correction for the hydrological effect was applied based on physical modeling of soil moisture. Finally, the effect of vertical displacement of the station (free-air effect) was corrected using GNSS data. After these corrections, residual gravity indicated a long-term increase, with its rate gradually decreasing with time. This fact suggests that viscoelastic relaxation after the earthquake played an important role in producing the long-term gravity changes. Fitting a decaying exponential function of time to the residual series yielded 89.4 ± 4.4 µGal (1 µGal = 10^–8 ms^–2) as the total gravity change and 635 ± 17 days as the characteristic time scale. In addition to the ground-based observations, the data from satellite gravity missions GRACE/GRACE-FO were analyzed to retrieve gravity changes at Mizusawa. Similar analysis of the satellite-based data yielded 15.2 ± 1.6 µGal as the total gravity change and 3444 ± 599 days as the characteristic time scale. The difference in the estimates of the total gravity change, of a factor of about 6, from the ground-based and the satellite-based observations may be attributed to the limited spatial resolution in the latter method. The difference in the estimates of the time scale, of a factor of about 1/5, may originate from the difference in the depth where the two kinds of gravimetry are mainly sensitive. Referring to recent theoretical studies on postseismic deformations after the 2011 Tohoku earthquake, our results can be interpreted consistently by assuming the existence of a layer of viscoelastic materials with viscosity $2\times {10}^{18} \text{Pa s}$ underneath the Tohoku area of Japan.

Graphical Abstract

Introduction

The great 2011 off the Pacific coast of Tohoku Earthquake (Mw9.0; hereafter abbreviated as the 2011 Tohoku earthquake) caused postseismic, as well as coseismic, deformations of the Earth (e.g., Sato et al. 2011; Kido et al. 2011; Ozawa et al. 2012; Munekane 2013; Tobita 2016; Fujiwara et al. 2022). Generally, three kinds of processes are considered to be responsible for inducing such changes; namely, (i) afterslips (e.g., Marone et al. 1991; Heki et al. 1997), (ii) viscoelastic relaxation (e.g., Thatcher et al. 1980; Pollitz et al. 2007), (iii) pore-pressure transients (e.g., Jónsson et al. 2003; Hu et al. 2014). These processes are characterized by different time scales. The pore-pressure transients are considered to have typical time constants around a few days, not exceeding one year (Panuntun et al. 2018), which is the shortest among the three. Effects of the afterslips are commonly considered to decay within a few years (Yamagiwa et al. 2015; Broerse et al. 2015). Viscoelastic relaxation is characterized by time scales of tens to hundreds of years, thus having the longest relaxation time constant among the three processes (Sun et al. 2014; Yamagiwa et al. 2015). Sun et al. (2014) and Suito (2017) studied numerically the process of viscoelastic relaxation of the 2011 Tohoku earthquake using 3-D finite-element methods. Takagi (2018) developed a theory for calculating coseismic and postseismic deformations for general viscoelastic Earth models and applied it to the 2011 Tohoku earthquake. Tobita (2016) and Fujiwara et al. (2022) investigated numerical functions of time to model the postseismic ground motions observed by Japanese GNSS network (GEONET).

Deformations of the Earth caused by large earthquakes also induce changes in the Earth’s gravity field. Coseismic and postseismic changes of gravity for the 2011 Tohoku earthquake have been investigated using the data from the satellite gravity missions GRACE and GRACE-FO (e.g., Han et al. 2011; Matsuo and Heki 2011; Tanaka and Heki 2014; Han, 2014; Cambiotti 2020). In contrast to the satellite gravity studies, there have not been so many publications on the ground-based gravity observations. Geospatial Information Authority of Japan surveyed gravity around Japan to construct Japanese Gravity Standardization Net 2016 (JGSN2016) (Yoshida et al. 2018), but it was not aimed to monitor the postseismic gravity changes after the 2011 Tohoku earthquake. Zhang et al. (2016) studied coseismic gravity changes using ground-based gravity measurements. Takagi (2018) examined his viscoelastic relaxation model using absolute gravity data obtained in the Tohoku area.

In this paper, we investigate the gravity changes after the 2011 Tohoku earthquake recorded by the superconducting gravimeters (SGs) operated at Mizusawa VLBI Observatory, National Astronomical Observatory of Japan (NAOJ; https://www.miz.nao.ac.jp/en.1.html). The SG is an extremely stable and precise relative gravimeter (Prothero and Goodkind 1968; Goodkind 1999; Hinderer et al. 2015). Although the observations with an SG at Mizusawa were started prior to the occurrence of the 2011 Tohoku earthquake, here we focus on the data acquired after the year 2014, because the instrument was seriously disturbed during the period from 2011 to 2014 by the main shock and many aftershocks. The spatial distributions and the temporal evolution of afterslips, viscoelastic relaxation and pore-pressure transients differ significantly depending on the data and the physical models used in various studies (Sun et al. 2014; Yamagiwa et al. 2015; Freed et al. 2017; Agata et al. 2019; Fukuda and Johnson 2021). Here we assume that most of the long-term gravity trend observed at Mizusawa after 2014 is controlled by linear processes of viscoelastic relaxation (Suenaga et al. 2023). To retrieve postseismic gravity signals associated with viscoelastic relaxation, we apply three kinds of corrections to the SG data. First, instrumental drift of the SG is estimated by using absolute gravity measurement data. Second, hydrological effects on gravity are evaluated by employing a physical model based on a vertical diffusion equation of soil moisture. Third, effects of vertical displacement of the gravity station are corrected using GNSS data. After these corrections, we obtain a gradually increasing trend of residual gravity, whose rate is gradually decreasing with time. Fitting a decaying exponential function of time to the residual, we obtain 89.4 ± 4.4 µGal as the total gravity change and 635 ± 17 days as the characteristic time scale. In addition to the SG data, we analyze the data from satellite gravity missions GRACE/GRACE-FO to retrieve gravity changes at Mizusawa in the global gravity field. Analysis of the data gives 15.2 ± 1.6 µGal as the total gravity change and 3444 ± 599 days as the characteristic time scale. We discuss the difference in these results from the viewpoint of different characteristics of the ground-based and the satellite-based methods of gravimetry. We also discuss the time scale in relation with viscoelastic relaxation after the earthquake.

Ground-based gravity observations

Gravity observations with superconducting gravimeters

In this study, we use the data from two SGs. Both SGs are of TT70 type with 200 L Dewar. The first SG with a serial number #007 (hereafter SG007) was operated at the Esashi Earth Tides Station of NAOJ from 1988 to 2008. In December 2008, it was moved to the Earth and Planets Experimental Building (formerly called the Absolute Gravity Measurement Building) in Mizusawa VLBI Observatory (Fig. 1c). The SG was installed on the base which was originally designed for an optical interferometer of a stationary type absolute gravimeter (Sugawa et al. 1979; Fig. 1d). The geographical location of the observation site is 39.1333 N, 141.1334 E and 63 m above the sea level. The compressor of the refrigeration system for the SG, which can be a noise source for gravity observations, was installed in a separate room divided around a buffering space for air conditioning, so that the SG was not seriously affected by ground vibration and waste heat from the compressor.

Fig. 1

a b Location of Mizusawa VLBI Observatory, National Astronomical Observatory of Japan. The red rectangle in a is blown up in the more detailed map (b). c Campus map of Mizusawa VLBI Observatory. The SG is housed inside the building “Obs. Building”. “GNSS1” and “GNSS2” are the GNSS observation sites used in this study. Wn, We, Ww and Ws are the ground water observation wells. d Gravity observation building. “SG” stands for the superconducting gravimeter, and “C” is the point where a compressor and a chiller are installed. Absolute gravity measurements were made at the point “AG”

After the 2011 Tohoku earthquake, the instrumental drift of SG007 became larger, and possibly associated with this, instrumental noise level also became higher than before. In addition, the boil-off rate of liquid helium increased significantly. In October 2017, we replaced SG007 with the second SG having a serial number #016 (hereafter SG016). SG016 was deployed at the Kamioka station of Institute for Cosmic Ray Research, The University of Tokyo, from 2004 to 2016. It was recooled at Tsukuba University in February 2017. After 8 months refurbishment works at Tsukuba, the gravimeter was transported to Mizusawa in October 2017, with the gravity sensor kept at cryogenic temperature. Due to the replacement of the SGs at Mizusawa, there was a gap in the data lasting for about 3 days. SG016 is now in operation at Mizusawa.

The instrumental sensitivity of SG007 was calibrated by comparison with absolute gravity measurements when it was operated at the Esashi station (Tamura et al. 2005). We assume that the estimated scale factor is invariable before and after moving the instrument from Esashi to Mizusawa. Sensitivity calibration for SG016 was also made by absolute gravity measurements in 2005 and 2007 at the Kamioka site (Sato 2007). The estimated scale factor of SG016 was also confirmed by absolute gravity measurements at Mizusawa (Tamura, unpublished results).

The gravity data from the SGs were processed with a standard method as follows. First, the data were checked for glitches, spikes and steps. Next, the original 1-s data were lowpass filtered and decimated to 1-h data. Then, tidal analysis was made using the program BAYTAP-G (Tamura et al. 1991). Estimated parameters of short-period tides are shown in Additional File 1. The residual time series data were further corrected for long-period tides by synthesizing theoretical tidal changes. The effect of polar motion was removed using the earth orientation data provided by IERS Earth Orientation Center (https://hpiers.obspm.fr/eop-pc/index.php).

Figure 2a shows the residual gravity for both SG007 (before October 2017) and SG016 (after October 2017). The offset in October 2017 is due to the replacement of the gravimeters, and its magnitude is arbitrary. Before October 2017, the residual gravity indicates almost linear decrease with time at the rate of approximately − 52 µGal/year. After October 2017, the rate of gravity decrease is significantly reduced to about – 15 µGal/year. Most of these changes are of instrumental, not natural, origin. In the following, we apply three stages of corrections to these data to derive final estimates of long-term gravity changes at Mizusawa, namely (i) instrumental drift, (ii) hydrological effects and (iii) free-air effects.

Fig. 2

a Residual gravity data after tidal analysis for the two SGs (SG007 and SG016) at Mizusawa. The vertical cyan line in October 2017 indicates the time when the gravimeters were replaced. The offset value in October 2017 is arbitrary. b Difference between the nominal gravity values from the SGs and the absolute gravity values from absolute gravimeters (black cross). Fitting a straight line to these data gives the rate of instrumental drift of the SGs. The data before and after the replacement were processed independently. c Residual gravity after correcting the instrumental drift of the SGs. A constant value of 980,146,200 µGal has been subtracted. e Detailed view of the data before and after the replacement of the SGs. d Water level changes observed at ground water observation well Wn

Correction of instrumental drift by absolute gravity measurements

In order to calibrate instrumental drift of the SGs, we performed absolute gravity measurements once a year from 2015 through 2021, using FG5 absolute gravimeters (Niebauer et al. 1995; Amalvict et al. 2001; Rosat et al. 2009) owned by Earthquake Research Institute, The University of Tokyo. The absolute gravimeter was installed in the next room of SG (Fig. 1d). The measurements lasted for at least 38 h for each experiment. The dropping interval was set to every 10 s, thus providing sufficient numbers of drops ~ 10⁴ from 2 days experiments (Table 1). We measured the local vertical gravity gradient at the site, and the obtained result (–3.064 µGal/cm) was input to the controller of the gravimeter. In situ atmospheric pressure was also recorded using either the built-in barometer or an external high-precision barometer (Digiquartz Series 1000, Paroscientific, Inc).

Table 1 Results of absolute gravity measurements with FG5 absolute gravimeters at Mizusawa

The raw gravity values from each drop were corrected for the following contributions: (i) tidal effect, (ii) barometric effect, (iii) polar motion effect. For tidal correction, the tidal parameters estimated from SG observations were used. This means that the tidal correction includes the contribution of oceanic tides. For barometric correction, we used the admittance of − 0.33 µGal/hPa. This is the average of the estimates of the atmospheric admittance obtained from tidal analysis for SG007 and SG016 (Additional File 1).

The results of gravity measurements are listed in Table 1. The values of gravity acceleration refer to the height of 1.30 m above the floor. The standard error of mean gravity arising from statistical fluctuation ranges from 0.1 to 0.3 µGal, all of which are successfully controlled under 2 µGal, the nominal accuracy of FG5 gravimeters (Niebauer et al. 1995; Okubo et al. 1997). Since these measurements employed plural gravimeters with different serial numbers, a special care must be taken for instrumental offsets. Specifically, the instrumental offset of FG5 #109 in 2018 had been confirmed to be + 8.6 µGal from a comparison with other instruments (S. Okubo, personal communications). Because this is substantially larger than the nominal accuracy of FG5, we corrected the data in 2018 by 8.6 µGal. Besides, the instrumental offset between FG5 #109 and #241 was estimated to be $1.6\pm 0.3$ µGal in March 2021, namely in between the last two experiments (Imanishi et al. 2021). Since it may be regarded small enough compared with the nominal accuracy, we did not take into account the instrumental offset except for the year 2018.

Assuming that the SG and the absolute gravimeter sensed the same temporal gravity changes, the difference between the gravity data from the two instruments should correspond to the instrumental drift of the SG. It is known that the instrumental drift of the SG is well approximated by either a linear or a decaying exponential function of time (Van Camp and Francis 2007). Aside from the initial exponential drift that may be seen just after first installation of the instrument (Hinderer et al. 2015), the time constant of the exponential decay is typically several years to ten years. Given the short time span of the present data, an exponential drift is almost indistinguishable from a linear drift. Therefore, here we modeled the drift of the SGs by linear functions of time. The SG data to be compared with the absolute gravity data were prepared by taking moving averages (typically 3 days long) in the period when absolute gravity measurements were performed. We used three data points (2015–2017) and four data points (2018–2021) for estimating the drift rates of SG007 and SG016, respectively.

Figure 2b shows the difference in the gravity data (cross signs) and the results of fitting linear functions of time (magenta lines). The estimated rate of instrumental drift was − 0.1358 ± 0.0042 μGal/day for SG007 and − 0.0330 ± 0.0024 μGal/day for SG016, respectively. The minus sign represents apparent gravity decrease with time. These rates correspond to annual gravity decrease of about 50 µGal (SG007) and about 12 µGal (SG016). Because the magnitude of instrumental drift for the SG of TT70 type ranges typically from 10 to 20 µGal/year (Hinderer et al. 2015), the drift of SG007 is larger than normal, whereas that of SG016 may be regarded as normal.

Once the instrumental drift of the SGs was calibrated by comparison with absolute gravity data, we can treat the continuous gravity recordings from the SGs in terms of absolute values of gravity acceleration. Figure 2c shows the absolute gravity changes (black lines) from the two SGs after correcting for instrumental drift. The data from absolute gravity measurements are also shown in red crosses. The misfit between the SG data and absolute gravity data is typically 2–3 µGal, comparable to the nominal accuracy of FG5. The inset in Fig. 2e shows the blowup of gravimeter replacement in October 2017. The gravity recordings produced by the two SGs are almost continuous before and after the replacement. Note that estimation of instrumental drift for the two SGs was done independently without imposing the condition of continuity in gravity at the point of the replacement. This fact implies that the overall control of the instrumental drift of the SGs with repeated absolute gravity measurements worked almost ideally, and also that SG016 was not seriously affected by initial drift just after the installation at Mizusawa in October 2017.

The gravity time series shown in Fig. 2c indicates long-term decrease in gravity. It is also evident that the data in the period before October 2017 contain more short-term events of gravity changes. Some of such events may be related with the slightly unstable instrumental conditions of SG007. In addition, correlation with the ground water level measured at the well Wn (Fig. 1c) may be noticeable, in particular in the data from SG016, as shown in Fig. 2d. This suggests necessity of correcting hydrological effects on gravity.

Correction of hydrological effects

Our next step of gravity data processing is the correction of hydrological effects (Fig. 3). In the case of the Mizusawa station, Kazama et al. (2012) succeeded in reproducing the hydrological gravity change observed by SG007 from 2009 to 2010 within 1.0 µGal, by solving a non-linear diffusion equation for the spatiotemporal variation in soil moisture numerically. The hydrological gravity change in and after 2011 was also calculated by Kazama et al. (in preparation) based on the physical model of Kazama et al. (2012), so we here utilize their calculation results to correct the hydrological gravity disturbance in the superconducting gravity time series observed in the period from 2014 to 2021.

Fig. 3

Hydrological gravity change at NAOJ Mizusawa from 2018 to 2020. a Cyan bars and blue lines indicate the daily and total annual precipitation ($R(t)$ and ${R}_{tot}(t)$), respectively, observed by a hyetometer at NAOJ Mizusawa. Orange lines indicate the total annual evapotranspiration (${E}_{tot}(t)$), calculated from the Penman (1948)’s model and meteorological data observed at NAOJ Mizusawa. b A black line indicates the water level of the unconfined aquifer ($h(t)$) under the ground of NAOJ Mizusawa. In this panel, $h(t)$ is shown with respect to the ground surface altitude (${z}_{s}$). c Cyan, blue and black lines indicate the calculated soil moisture variations ($\theta (z,t)$) at the depths of 15–20, 45–50 and 95–100 cm from the ground surface, respectively. ${\theta }_{max}$ and ${\theta }_{min}$ are the maximum and minimum soil moisture values for NAOJ Mizusawa (= 0.52 and 0.28 m³/m³), respectively (see Kazama et al. 2012). d Gray and red solid lines indicate the hydrological gravity changes of ${g}_{w}\left(t\right)$ and ${g}_{w}^{out}\left(t\right)$ (Eqs. (1) and (2)), respectively, calculated from the spatial integrals of $\theta (z,t)$. Gray and red dashed lines indicate the averages of ${g}_{w}\left(t\right)$ and ${g}_{w}^{out}\left(t\right)$, and they are oriented at 0 and 10 μGal, respectively, so as not to overlap ${g}_{w}\left(t\right)$ with ${g}_{w}^{out}\left(t\right)$ in this panel. e Black and gray lines indicate the superconducting gravity data, obtained by the SG016 (${g}_{obs}(t)$; Fig. 2c). Yellow triangles indicate the observed hydrological gravity change, which was identified by comparing the time variation in ${g}_{obs}(t)$ with those in $R(t)$, $h(t)$ and $\theta (z,t)$. A red line indicates the superconducting gravimeter data after the hydrological effect of ${g}_{w}^{out}\left(t\right)$ was subtracted from ${g}_{obs}(t)$

Figure 3c illustrates the calculated soil moisture in the period from 2018 to 2020. The underground soil moisture distribution (denoted by $\theta (z,t)$; Fig. 3c) was obtained by numerically solving the non-linear diffusion equation (the so-called Richards equation; e.g., Jury and Horton 2004) under the following two boundary conditions: (i) the vertical Darcy flux on the ground is equal to the observed effective precipitation (= precipitation minus evapotranspiration; see Fig. 3a), and (ii) soil void is saturated with water under the observed height of water table (Fig. 3b). Soil moisture increases during precipitation rapidly and decreases after precipitation exponentially with time. The soil moisture variation is smaller and delayed at the deeper position, because soil moisture infiltrates from the ground surface to the deeper layer through the diffusion process.

The gray line in Fig. 3d shows the calculated hydrological gravity changes ${g}_{w}\left(t\right)$, obtained by spatially integrating the soil moisture distribution assuming horizontal homogeneity:

$$g_{w} \left( t \right) = \,2\pi \rho_{w} G\int {\theta \left( {z,t} \right)dz,}$$

(1)

where ${\rho }_{w}$ and $G$ indicate water density and the gravitational constant, respectively. ${g}_{w}\left(t\right)$ is to be compared with the observed gravity changes ${g}_{obs}(t)$ shown in the black line in Fig. 3e. Although the calculated gravity changes well reproduce several features of the observed gravity changes, ${g}_{w}\left(t\right)$ tends to overestimate the undulations of ${g}_{obs}(t)$. This can be mitigated by considering the umbrella effect of the gravimeter building (e.g., Creutzfeldt et al. 2008; Reich et al. 2019). By limiting the integration area in the right-hand side of Eq. (1) to the outside of the building, we can define the hydrological gravity changes due to the soil moisture outside the gravimetry building as:

$$g_{w}^{out} \left( t \right) = \,2\pi \rho_{w} G\int {a_{out} \left( z \right)\theta \left( {z,t} \right)dz,}$$

(2)

where ${a}_{out}\left(z\right)$ indicates the integration ratio for the outside of the building, determined from the geometry of the building and the position of the gravimeter (Fig. 8 of Kazama et al. 2012). Variation of ${g}_{w}^{out}\left(t\right)$, shown in the red line in Fig. 3d, is smaller than that of ${g}_{w}\left(t\right)$ by approximately 64%. ${g}_{w}^{out}\left(t\right)$ better reproduces the observed gravity changes as indicated by inverted triangles in Fig. 3e.

By subtracting ${g}_{w}^{out}\left(t\right)$ from ${g}_{obs}(t)$, we obtain the time series data shown by the red line in Fig. 3e. Now the short-term undulations in the gravity time series are significantly reduced. The same method of correction is applied to the data from 2014 to 2021.

Correction of free-air effects by GNSS observations

Changes in the vertical position of the gravity station cause changes in the distance from the station to the center of the Earth, thus resulting in changes in Newtonian gravity acceleration (the free-air effect). The gravity trend after hydrological corrections shown in Fig. 3e is dominated by this effect. We make correction of this effect using the data from in situ GNSS observations.

Figure 1c shows the location of the GNSS stations (denoted as GNSS1 and GNSS2) used in this study. The station GNSS1 is closer to the SG site, at the distance of approximately 80 m. Daily coordinates based on ITRF2014 for these stations were routinely determined using GIPSY-OASIS II ver 6.4 software. In November 13, 2020, data acquisition at GNSS1 was stopped because of the failure of the antenna. Therefore, here we use the data from GNSS2 for the period after November 13, 2020. The relative position between the two stations was also determined routinely, so that consistency of the whole dataset was guaranteed.

Figure 4a shows vertical coordinates by GNSS at Mizusawa. A constant value of 105.6 m has been subtracted from the data. Because the daily coordinates have too large scatter, 7-day running means of the daily values were calculated. The coseismic subsidence caused by the 2011 Tohoku earthquake was about 100 mm at Mizusawa (not shown in Fig. 4a). After the earthquake, Mizusawa has been steadily uplifting. As of 2021, the postseismic uplift has already raised the ground of Mizusawa to the level higher than that prior to the earthquake.

Fig. 4

a (black line) Vertical displacement obtained from the daily solutions of GNSS observations at Mizusawa from 2014 to 2021. (red line) Result of fitting a numerical function to the GNSS data based on the method of Tobita (2016). b (gray line) Residual gravity data observed at Mizusawa. (black line) The same but after the correction of free-air effects. The red curve represents the result of fitting Eq. (4) to the data

The effect of this uplift of the ground should be corrected from the residual gravity observed at Mizusawa (gray line in Fig. 4b). In the long time scales, the rate of uplift has been gradually decreasing up to 2021. In the shorter time scales, there is no clear correlation between the vertical coordinates and the residual gravity. Rather than using the GNSS data for directly correcting the gravity data, we modeled the vertical coordinates with a simple function of time to extract the long-term trend of the GNSS observations. Tobita (2016) modeled the postseismic changes in the GEONET (the Japanese nationwide GNSS network) coordinates with a combination of one or two logarithmic function(s), one or two exponential function(s) and a linear function. Tobita (2016) assumed the same time constants of the model functions for all the GNSS stations. After some trials, we found that applying Eq. (5) of Tobita (2016), that is, using a combination of a logarithmic function, an exponential function and a linear function gave reasonable results in the case of our data. Employing two exponential functions or neglecting a logarithmic function changed the results slightly, but the difference in the fitting residual existed mostly in the beginning and the end of the time series. So, our model function of the long-term trend is

$$u_{z} \left( t \right) = a^{\left( Z \right)} + b^{\left( Z \right)} \left( {t - t_{0} } \right) + \zeta_{L}^{\left( Z \right)} \log \left( {1 + \frac{{t - t_{0} }}{{\tau_{L}^{\left( Z \right)} }}} \right) + \zeta_{E}^{\left( Z \right)} \left[ {1 - \exp \left( { - \frac{{t - t_{0} }}{{\tau_{E}^{\left( Z \right)} }}} \right)} \right],$$

(3)

where ${u}_{z}(t)$ is the observed ellipsoidal height at the time $t$, ${t}_{0}$ is the origin time of the 2011 Tohoku earthquake, ${a}^{(Z)}$ and ${b}^{(Z)}$ are the coefficients of constant and linear terms, ${{\zeta }_{L}}^{(Z)}$ and ${{\tau }_{L}}^{(Z)}$ are the magnitude and the time constant of the logarithmic term, and ${{\zeta }_{E}}^{(Z)}$ and ${{\tau }_{E}}^{(Z)}$ are the magnitude and the time constant of the exponential term. For the time constants, we took the same values as used in Tobita (2016), thus ${{\tau }_{L}}^{(Z)}=10.3 \mathrm{days}$ and ${{\tau }_{E}}^{(Z)}=3306.0 \mathrm{days}$. Fitting Eq. (3) to the data with a least-squares method gave ${a}^{(Z)}=374\pm 54 \mathrm{mm}$, ${b}^{(Z)}=- 0.0866\pm 0.0042 \mathrm{mm}/\mathrm{day}$, ${{\zeta }_{L}}^{(Z)}=- 185\pm 17 \mathrm{mm}$ and ${{\zeta }_{E}}^{(Z)}=1410\pm 78 \mathrm{mm}$. The modeled function is shown by the red curve in Fig. 4a.

To correct the free-air effect from the observed gravity changes, we multiplied the modeled vertical coordinates by a coefficient –0.3086 μGal/mm and subtracted them from the residual gravity. This gave the data shown in the black curve in Fig. 4b. This is our final estimate of the ground-based gravity changes obtained by the SGs at Mizusawa.

The overall shape of the residual gravity shown in Fig. 4b indicates a long-term increase with its rate gradually decreasing with time. This fact suggests us to model the temporal changes with a combination of a linear function and an exponentially decaying function of time written as:

$$g^{\left( G \right)} \left( t \right) = a^{\left( G \right)} + b^{\left( G \right)} \left( {t - t_{0} } \right) + \zeta^{\left( G \right)} \left[ {1 - \exp \left( { - \frac{{t - t_{0} }}{{\tau^{\left( G \right)} }}} \right)} \right],$$

(4)

where ${g}^{(G)}(t)$ is the ground-based gravity at the time $t$, ${t}_{0}$ is the origin time of the 2011 Tohoku earthquake, ${a}^{(G)}$ and ${b}^{(G)}$ are the coefficients of constant and linear terms, and ${\zeta }^{(G)}$ and ${\tau }^{(G)}$ are the magnitude and the time constant of exponential decay. ${a}^{(G)}$, ${b}^{(G)}$, ${\zeta }^{(G)}$ and ${\tau }^{(G)}$ are the unknown parameters to be adjusted. Applying a non-linear least-squares method, however, we found that the problem was numerically ill-conditioned, because the correlation between ${b}^{(G)}$ and ${\tau }^{(G)}$ was very high. So, omitting the linear term, we obtained ${\zeta }^{(G)}=89.4\pm 4.4 \mathrm{\mu Gal}$ and ${\tau }^{(G)}=635\pm 17\mathrm{ days}$. This means that postseismic gravity changes as large as about 89 µGal are going on at Mizusawa with a time constant of about 1.7 years. The fitted function is shown by the red curve in Fig. 4b.

Satellite-based gravity observations

To interpret gravity changes from the ground-based observation at Mizusawa, we analyzed the satellite-based gravity changes from Gravity Recovery and Climate Experiment (GRACE) and its follow-on (GRACE-FO). The GRACE/GRACE-FO data used in this study are the Release 6 Level-2 monthly products from three analysis centers: The University of Texas—Center for Space Research, the Jet Propulsion Laboratory, and the GeoForschungZentrum Potsdam. Here we used the ensemble mean of the three GRACE/GRACE-FO datasets to reduce the random errors caused by the different analysis strategy of each analysis center (Sakumura et al. 2014). The data cover the period from April 2002 to December 2022, consisting of 214 monthly datasets with gaps for 35 months. Following Han et al. (2014), we used the spherical harmonic coefficients up to degree and order 40, which represent the Earth’s gravity field at a spatial resolution of approximately 500 km. As commonly done in GRACE/GRACE-FO data processing, we replaced the Earth’s oblateness values (C₂₀) with those determined by Satellite Laser Ranging (Cheng et al. 2013) considering the estimation accuracy. Also, we used the degree-one components (C₁₀, C₁₁ and S₁₁) estimated by the GRACE-OBP approach (Swenson et al. 2008) because GRACE/GRACE-FO alone cannot measure them.

The hydrological effects on the GRACE/GRACE-FO data were corrected using the Global Land Data Assimilation System version 2.1 (GLDAS-2.1) hydrological model (Rodell et al. 2004; Li et al. 2019). We used the ensemble mean of three land surface process models (NOAH, CLSM and VIC) from GLDAS-2.1 to enhance their common hydrological signatures. We corrected for the hydrological effects (soil moisture, snow, and canopy water) by expanding equivalent water depth in GLDAS-2.1 to spherical harmonic coefficients and subtracting them from the GRACE/GRACE-FO data.

We converted the spherical harmonic coefficients to the gravity disturbances using Eq. (2–153) of Heiskanen and Moritz (1967), and then obtained the time-series of monthly gravity changes by subtracting from them the mean gravity disturbance in the period from April 2002 to April 2003. We also derived the errors of the obtained gravity changes from formal errors of spherical harmonic coefficients using Eq. (4) of Wahr et al. (2006). To suppress the correlated noise that arises from aliasing errors, we applied the decorrelation filter proposed by Kusche et al. (2009) to the GRACE/GRACE-FO data. Here we tested the application of eight filters (DDK1–8) and no filter to the GRACE/GRACE-FO data to choose the best filter for capturing earthquake signatures.

Figure 5 shows the time-series of the satellite-based gravity changes at Mizusawa (39.1°N, 141.1°E) from GRACE/GRACE-FO observations with eight different filters (DDK1–8) and no filter. Each time-series shows the clear signatures of annual, coseismic and postseismic gravity changes. The amplitude of the annual gravity changes in each time-series is about 1 μGal. The maximum and minimum peaks of the annual changes were found around February and August, respectively. Such changes are likely to be caused mainly by the snow masses in winter, which were not corrected by using GLDAS-2.1. In the northern part of Japan, especially on the Sea of Japan side, snow accumulations of over 1 m are observed every year (e.g., Fujiwara et al. 2022). The coseismic and postseismic signatures at Mizusawa show the stepwise gravity decrease at the occurrence of the earthquake followed by gradual gravity increase. The coseismic gravity decrease can be explained mainly by dilatation of crust above the down-dip edge of the fault (Han et al. 2006, 2011; Matsuo and Heki 2011), whereas the postseismic gravity increase is interpreted as viscoelastic relaxation of the coseismic stress in the lower crust and upper mantle (e.g., Han et al. 2014).

Fig. 5

Time-series of the satellite-based gravity changes at Mizusawa (39.1°N, 141.1°E) from GRACE/GRACE-FO observation with eight different filters (DDK1–8) and no filter. Open circles denote the gravity changes relative to April 2002–April 2003 time-mean, and error bars denote the observation errors. Blue lines show the results of fitting the functions given by Eq. (5) to the gravity data. Red lines show the same functions but without periodic components

We modeled the data shown in Fig. 5 with a combination of a static, linear, annual, semiannual, 161-day periodic function, and coseismic and postseismic functions of the 2011 Tohoku earthquake (e.g., Cambiotti 2020) as given by

$$g^{\left( S \right)} \left( t \right) = a^{\left( S \right)} + b^{\left( S \right)} \left( {t - t_{0} } \right) + \mathop \sum \limits_{k = 1}^{3} \left( {c_{k} \cos \omega_{k} t + d_{k} \sin \omega_{k} t} \right) + H\left( {t - t_{0} } \right)\left[ {\xi^{\left( S \right)} + \zeta^{\left( S \right)} \left\{ {1 - \exp \left( { - \frac{{t - t_{0} }}{{\tau^{\left( S \right)} }}} \right)} \right\}} \right],$$

(5)

where ${g}^{(S)}(t)$ is the satellite-based gravity at the time $t$, ${t}_{0}$ is the origin time of the 2011 Tohoku earthquake, ${a}^{(S)}$ and ${b}^{(S)}$ are the coefficients of constant and linear terms, ${c}_{k}$ and ${d}_{k}$ are the coefficients of sinusoidal functions of angular frequencies ${\omega }_{k}$, $H(\bullet )$ is the Heaviside step function, ${\xi }^{(S)}$ and ${\zeta }^{(S)}$ are the magnitude of the coseismic and postseismic gravity changes, and ${\tau }^{(S)}$ is the time constant of the postseismic function. For the sinusoidal functions, annual (k = 1), semiannual (k = 2), 161-day periodic (k = 3) functions were used. The 161-day periodic function represents the aliasing effect due to the S₂ tidal wave inherent to the GRACE/GRACE-FO data (Ray and Luthcke 2006). The coefficients and the time constant were determined by a weighted non-linear least-squares method. The errors of the obtained gravity changes were used as weights of the non-linear least-squares fitting.

Table 2 shows the estimates of the coseismic and postseismic gravity changes and the time constant of postseismic gravity changes at Mizusawa from GRACE/GRACE-FO obtained with eight different filters (DDK1–8) and no filter. The coseismic and postseismic gravity changes ${\xi }^{(S)}$ and ${\zeta }^{(S)}$ were well determined by a weighted non-linear least-squares method for any filters. The estimate of ${\zeta }^{(S)}$ varied from + 6.3 to + 19.4 μGal, depending on the choice of the filter. On the other hand, the time constant ${\tau }^{(S)}$ of the postseismic gravity changes was not well constrained for the DDK1 filter. In the case of DDK1, the postseismic signature might be affected by significant omission errors due to the strong filter. For other filters, the estimate of ${\tau }^{(S)}$ was well constrained, ranging from 3094 to 3656 days Considering that the magnitude of the estimation error of ${\tau }^{(S)}$ is smallest in the case of the DDK3 filter, here we adopt the solution for DDK3. For this solution, we obtain ${\zeta }^{(S)}=15.2\pm 1.6\mathrm{ \mu Gal}$ and ${\tau }^{(S)}=3444\pm 599\mathrm{ days}$. In other words, the time constant is about 9.43 years.

Table 2 Results of fitting Eq. (5) to the gravity data at Mizusawa (39.1°N, 141.1°E) from GRACE/GRACE-FO processed with eight different filters (DDK1–8) and no filter

Discussion

Comparison between ground-based and satellite-based gravity observations

We have so far derived postseismic gravity changes at Mizusawa after the 2011 Tohoku earthquake as viewed from ground-based and satellite-based gravimetric methods. Because we have corrected the free-air effects for the ground-based observations, the gravity changes derived by the two approaches can be compared directly. To summarize the results of the previous sections, the ground-based observations yielded ${\zeta }^{(G)}=89.4\pm 4.4 \mathrm{\mu Gal}$ and ${\tau }^{(G)}=635\pm 17\mathrm{ days}$ for the magnitude and the time constant of the exponentially decaying signals, whereas the satellite-based observations yielded ${\zeta }^{(S)}=15.2\pm 1.6 \mathrm{\mu Gal}$ and ${\tau }^{(S)}=3444\pm 599\mathrm{ days}$ for the corresponding parameters. In other words, the magnitude of postseismic gravity change from the ground-based observations is about 6 times larger than that from the satellite-based observations. Conversely, the time constant from the satellite-based observations is about 5 times larger than that from the ground-based observations. Although there are correlations between the estimates of the magnitude and the time constant, it is a robust feature that the ground-based observations give much larger postseismic gravity changes with a much smaller time constant than the satellite-based observations.

We consider that the main cause of the difference in the magnitude of postseismic gravity changes from the ground-based and satellite-based observations is the spatial resolution of each gravimetry. On one hand, the ground-based observations provide information on temporal gravity changes at a particular site on the surface of the Earth, and are, in principle, sensitive to the changes in the gravity field up to infinitely large spatial wavenumbers. Whereas this feature enables a wide variety of investigations using the data from ground-based observations, it is often necessary to correct the observed data for site-specific local effects such as the changes in underground water table, as we did in this study. On the other hand, the satellite gravity observations such as GRACE/GRACE-FO provide information on the global gravity field in terms of spherical harmonic expansion. Truncating the expansion at finite degrees and orders and applying spatial filters to avoid spectral leakage result in limited spatial resolution, as we saw in the previous section. In our case, the spatial resolution is no better than 500 km. Even if there are gravity changes with spatial scales smaller than 500 km, they will be filtered out, resulting in reduced estimates of gravity signals. So, we conclude that the difference in the estimated magnitude of postseismic gravity changes is due to the different spatial resolutions of the two kinds of gravimetry.

To confirm this conclusion, we investigated crustal deformations in the Tohoku area detected by GNSS. Figure 6a shows the vertical displacements in the period from February 2011 to February 2021. The data were taken from the F5 solutions of GEONET (Takamatsu et al. 2023). The east–west sections surrounded by the green rectangle are expanded in Fig. 6c to extract horizontal distributions of vertical motions in the selected periods. As shown by the black line in Fig. 6c, all the stations indicated coseismic subsidence. The amount of subsidence ranged from about 33 mm at 140.4 deg E to about 754 mm at 141.75 deg E. After that, the stations to the west of 140.85 deg E continued to subside, whereas the stations to the east of 140.85 deg E indicated uplift. As a result of these deformations, there are nodal lines around 140.8–141.0 deg E and 141.4–141.6 deg E, both approximately in the north–south direction, in the total vertical displacements from February 2011 to February 2021 as shown in Fig. 6a. The gravity station of Mizusawa (denoted as NAO) is located in the region between the two nodal lines, where the total vertical motion is positive (upward). The typical horizontal scale of the spatial pattern of uplift/subsidence is about 50 km, well smaller than the spatial resolution of satellite gravity missions. Similarly, Fig. 6b, d shows the vertical displacements in the period from February 2014 to February 2021. This period is approximately the same as the period of the SG data used in this study. Although the nodal line around 140.8–141.0 deg E in Fig. 6a still exists, it appears to have shifted slightly westward. The nodal line around 141.4–141.6 deg E in Fig. 6a has disappeared in Fig. 6b, d. Even so, a spatial pattern of uplift/subsidence with horizontal wavelengths of about 100 km can be recognized in Fig. 6b, d.

Fig. 6

Vertical displacements near Mizusawa from GEONET (F5 solutions). a Spatial distribution from February 2011 to February 2021. b Spatial distribution from February 2014 to February 2021. c Longitudinal profile in the area denoted by green rectangle in a. d Longitudinal profile in the area denoted by green rectangle in b

Although crustal deformations observed by GNSS cannot be directly compared with gravity changes, basically surface deformations and gravity changes arise from the same causes, in other words, motions of materials in the deep interior. The data shown in Fig. 6 strongly suggest that postseismic gravity changes at the surface of the Earth also indicate spatial distributions of similar wavelengths (50–100 km). At any place in the Tohoku area, long-term gravity changes similar to those captured at Mizusawa would have been obtained if continuous gravity observations were made, but the magnitude and the sign of such changes would have varied from place to place. Indeed, results of repeated absolute gravity measurements at the Honjo station of Tohoku University (shown by asterisk in Fig. 6) indicate different characteristics from those at Mizusawa (Additional File 2). On the other hand, satellite-based observations sense the spatially smoothed gravity distributions with typical resolution of 500 km. Figure 7 shows the gravity changes near Japan obtained from GRACE/GRACE-FO for the DDK3 filter, in the periods corresponding to those in Fig. 6. Spatial patterns found in Fig. 6 cannot be found in the gravity maps in Fig. 7. This contrast in the spatial resolutions will provide the reason why our analysis of the satellite-based data yielded a much smaller estimate of postseismic gravity changes at Mizusawa compared with the ground-based observations.

Fig. 7

Spatial distribution of gravity changes near Japan obtained from GRACE/GRACE-FO for the DDK3 filter. Green rectangles denote the areas shown in Figure 6 a and b. a From February 2011 to February 2021. b From February 2014 to February 2021

Time constants of viscoelastic relaxation

Next, we discuss the time constants of the exponentially decaying gravity signals determined from the ground-based and the satellite-based observations from the viewpoint of linear viscoelastic materials. In the Maxwell material, the relaxation time ${\tau }_{M}$ is defined as ${\tau }_{M}=\eta /\mu$, where $\eta$ is the coefficient of viscosity and $\mu$ is the rigidity. To interpret our observations in relation with viscoelastic properties of subsurface materials, we refer to two recent studies of numerical calculations by Suito (2017) and Takagi (2018) in the following.

Suito (2017) applied a 3-D finite-element method to model the postseismic deformation after the 2011 Tohoku earthquake observed by GPS-acoustic techniques as well as onshore GNSS. In the model, the mantle wedge and the oceanic mantle were treated as viscoelastic, whereas the crust and the subducting Pacific plate were elastic. By examining several models having different distributions of viscosities, Suito (2017) found that his “Model 3” reproduced the observations fairly well. The viscosities in that model were $2\times {10}^{18} \text{Pa s}$ in the mantle wedge, $1\times {10}^{19} \text{Pa s}$ in the oceanic mantle, and $1\times {10}^{18} \text{Pa s}$ in the thin weak layer below the subducting Pacific plate. To further improve the model, Suito (2017) examined his “Model 4” in which viscosities were variable with depth, but the difference between the reproduced displacement fields from Model 3 and Model 4 is not very dramatic. Suito (2017) did not give calculation results of gravity changes. Based on Model 3 of Suito (2017), there is a contrast in the viscosities of the mantle wedge and the oceanic mantle of a factor of 5. Assuming common rigidities for the whole area, there will be also a contrast in the relaxation time ${\tau }_{M}$ of these regions of a factor of 5. Our ground-based gravity observations (performed in the mantle wedge side) are more sensitive to the physical properties of the materials in the shallower wedge mantle, whereas satellite-based gravity observations are more sensitive to the materials in deeper and wider areas, therefore with higher viscosities (Broerse et al. 2015). This fact will provide a partial explanation of the difference in the time constants determined by the ground-based and the satellite-based observations of a factor of 5.

Takagi (2018) developed a general theory for calculating coseismic and postseismic deformations for the spherically symmetric and viscoelastic (Maxwell) Earth. The theory enables calculation of temporal evolution of internal deformations as well as surface deformations, for arbitrary distribution of viscosities in the Earth. Gravity changes at the surface of the Earth are calculated using internal distribution of volumetric dilatations. Because the theory assumes that the internal structure of the Earth is spherically symmetric, horizontal inhomogeneities at the plate subduction zone as in Suito (2017) are not taken into account. Takagi (2018) applied his theory to the 2011 Tohoku earthquake. In doing so, the outermost layer of 30 km thickness was treated as elastic, and beneath it a viscoelastic layer of 290 km thickness with viscosity $5\times {10}^{18} \text{Pa s}$ was assumed. Takagi (2018) stated that this viscosity was so chosen that the relaxation time of the viscoelastic layer was 2–3 years. As representative results of postseismic gravity rates, we present Fig. 8a, b, which was reproduced from Fig. 5.7d and Fig. 5.8h of Takagi (2018). Effects of vertical displacements on gravity have been corrected in these data. Figure 8a shows the calculated gravity changes in the one-year period from 2014 to 2015. In a wide area of the Tohoku area, the model predicts positive gravity changes. Mizusawa is located in the region between the contours of + 2 and + 5 µGal. On the other hand, our data shown in Fig. 4b indicate a gravity increase of about 7 µGal in the one-year period from 2014 to 2015, larger than the calculation. Also in the two-year period from 2021 to 2023, the calculation shown in Fig. 8b predicts a value between + 2 and + 5 µGal at Mizusawa. If the viscosity used in the model was $2\times {10}^{18} \text{Pa s}$ as in Suito (2017) instead of $5\times {10}^{18} \text{Pa s}$, the relaxation time of the viscoelastic materials was 40% of the original, and the calculation result for the period from 2021 to 2023 can be regarded as that for the 0.8-year period from 2015.2 to 2016.0 (note that the earthquake origin time is 2011.2). Our data shown in Fig. 4b indicate a gravity increase of about 3.5 µGal in the period from 2015.2 to 2016.0, consistent with the theoretical calculation.

Fig. 8

Reproduced from Figs. 5.7(d) and 5.7(h) in Takagi (2018)

Theoretical gravity rates after the 2011 Tohoku earthquake calculated by Takagi (2018). a One-year interval from 3.0 to 4.0 years. b Two-year interval from 10.0 to 12.0 years.

These considerations lead to the conclusion that the magnitude and the characteristic time scale of the postseismic gravity changes from our ground-based gravity observations is consistent with the existence of a layer of viscoelastic materials with viscosity of $2\times {10}^{18} \text{Pa s}$ underneath the Tohoku area. This may be compared with the result of Han et al. (2014), who investigated the postseismic gravity changes using satellite gravity data and found that the signals were best modeled by the biviscous relaxation with a transient and steady state viscosity of 10¹⁸ Pa s and 10¹⁹ Pa s, respectively, for the asthenosphere. Although we do not consider biviscous models in this paper, Fig. 2 of Han et al. (2014) suggests that a Maxwell model with viscosity between $1\times {10}^{18} \text{Pa s}$ and $5\times {10}^{18} \text{Pa s}$ gives reasonable fit to the data. In this sense, we may conclude that our result is consistent with the result of Han et al. (2014).

Concluding remarks

In this paper, we analyzed the data from the superconducting gravimeters at Mizusawa, Japan, as well as the satellite gravity missions GRACE/GRACE-FO for retrieval of postseismic gravity changes after the 2011 Tohoku earthquake. Both methods of gravimetry yielded exponentially decaying signals reflecting viscoelastic relaxation after the quake, but the magnitude and the time constant were much different between the ground-based and satellite-based methods. We interpreted the difference in the magnitude as the consequence of the different spatial resolution of the two methods of gravimetry. The difference in the time constant was ascribed to the difference in the depth where the ground-based and satellite-based gravimetries are mainly sensitive. Our results suggest the existence of a layer of viscoelastic materials with viscosity of $2\times {10}^{18} \text{Pa s}$ underneath the Tohoku area of Japan.

As far as we are aware, this is the first report of continuous gravity changes after a megathrust earthquake recorded by a superconducting gravimeter located not far from the source region. We hope that the present result has proved advantage of ground-based gravity observations, especially with the SG, in that they enable detection of local-scale deformations of the Earth which are invisible by the satellite methods. The gravity signals revealed in this study surely reflect some features of stress relaxation which are taking place underneath Japan, but more detailed theoretical modeling of such processes is beyond the scope of this paper. We plan to extend our investigations to other datasets of gravity monitoring in Japan from superconducting gravimeters as well as relative gravimeters.