Impact of the solid Earth mass adjustment by the 2011 Tohoku–Oki earthquake on the regional sea level and hydrological mass change recovery from GRACE

SUMMARY

For more than a decade, GRACE data have provided global mass redistribution measurements due to water cycles, climate change and giant earthquake events. Large earthquakes can yield gravity changes over thousands of kilometres from the epicentre for years to decades, and those solid Earth deformation signals can introduce significant biases in the estimate of regional-scale water and ice mass changes around the epicentres. We suggest a modelling scheme to understand their contribution to the estimates of water and ice mass changes and to remove the earthquake-related solid mass signals from GRACE data. This approach is composed of physics-based earthquake modelling, GRACE data correction and high-resolution surface mass change recovery. In this study, we examined the case of the 2011 Tohoku–Oki earthquake to better estimate the regional sea level and hydrological mass changes in the East Asia. The co- and post-seismic changes from GRACE observations were used to constrain the earthquake model parameters to obtain optimal self-consistent models for the earthquake source and the asthenosphere rheology. The result demonstrated that our earthquake correction model significantly reduced the mass change signals by solid Earth deformation from the time-series of regional surface mass changes on both land and oceans. For example, the apparent climate-related ocean mass increase over the East Sea was 1.59 ± 0.11 mm yr⁻¹ for 2003–2016, significantly lower than the global mean ocean mass trend (2.04 ± 0.10 mm yr⁻¹) due to contamination of the earthquake signals. After accounting for the solid mass changes by the earthquake, the estimate was revised to 1.87 ± 0.11 mm yr⁻¹, that is increased by 20 per cent and insignificantly different from the global estimate.

1 INTRODUCTION

Gravity measurements are sensitive to various mass redistribution processes due to climate change and geophysical events. Water and ice mass loading cycles on the surface of the Earth are the main contributors to currently observed large-scale gravity variations, and the persistent ice mass loss signals found in the polar and high mountain regions are attributed by the global warming (Cazenave & Chen 2010; Velicogna & Wahr 2013). On the other hand, geophysical processes such as Glacial Isostatic Adjustment (GIA) and earthquakes cause mass redistribution within the solid Earth, resulting in significant gravity changes in the affected regions. Since 2002 April, gravity measuring satellites of GRACE (Gravity Recovery and Climate Experiment, Tapley et al. 2004) and GRACE Follow-On (Landerer et al. 2020) missions have monitored comprehensive gravity changes due to a number of geophysical and climate processes. The data are now importantly used in various studies of sea level changes (e.g. Cazenave et al. 2014; Adhikari et al. 2019; Jeon et al. 2021), terrestrial water and ice mass changes (e.g. Kim et al. 2019), planetary motion of the Earth (e.g. Sun et al. 2016; Seo et al. 2021), large earthquakes (e.g. Han et al. 2006; Cambiotti & Sabadini 2012; Han et al.2014), etc.

An important application using GRACE data is to examine water and ice mass variations associated with the contemporary climate change. Such estimates are provided only after removing the contribution by the solid Earth mass redistribution from full mass change signals inferred by GRACE geopotential observations. For example, as one of the important geophysical processes leading to significant gravity changes, the correction for the GIA effect is now routinely included in GRACE data processing. The effect gradually changes both geoid height and the ocean floor, and the contribution has significant impact on both global and regional surface mass change estimate (Peltier et al. 2018). On the other hand, corrections for earthquake-induced signals have not been widely considered, probably because of the small impact on the global-scale estimates as tested in a previous study (Tang et al. 2020). However, the contributions are not negligible in the regional-scale surface mass change estimates. Earthquake-induced gravity changes affect broad region up to a few thousands of kilometres from the epicentre. Further, giant seismic events are often followed by nonlinear temporal changes due to viscoelastic relaxation of the asthenosphere (Han et al. 2013). These signals can last for decades, introducing significant uncertainty to regional surface mass changes around epicentres of large earthquake events.

In particular, the earthquake signals may become a critical error source when applying leakage correction methods for GRACE data. Since GRACE level-2 (L2) monthly solutions contain intense longitudinal stripe-shape errors, post-processing methods like Gaussian smoothing (Wahr et al. 1998) and decorrelation filters (Swenson & Wahr 2006) are commonly used for the reduction. Such processes result in spreading and attenuation of gravity signals (i.e. signal leakage effect), and the low-resolution GRACE data have larger uncertainties in estimates of the mass changes in smaller regions. To overcome this limitation, various methods for the signal recovery have been proposed (Swenson & Wahr 2002; Horwath & Dietrich 2009; Chen et al. 2013; Dobslaw et al. 2020) including GRACE Mascons (Luthcke et al. 2013; Watkins et al. 2015; Save et al. 2016). Such signal recovery methods mainly adjust the distribution of water and ice mass changes inferred from the GRACE data and find an optimal solution of the surface mass distribution with reduced signal leakages. Thus, errors from uncorrected solid Earth signals included in the GRACE data consequently propagate into the refinement version. Several studies analysed residual signals due to important global-scale solid Earth responses including GIA effect and polar motion, and have discussed the impact on a signal recovery process (e.g. Jeon et al. 2018; Kim et al. 2022). Similarly, residual gravity changes due to earthquakes also yield significant regional-scale errors during the process.

For an example, we displayed the impact of the 2011 Tohoku–Oki earthquake signals on an analysis of regional surface mass changes in Fig. 1. Fig. 1(a) shows 1-yr mean gravity difference before and after the event. Specifically we used the Release 06 L2 data from University of Texas Center for Space Research (CSR, Tapley et al. 2019) with the Gaussian smoothing (500 km) and a decorrelation filter (P4M10). Although earthquake signals are attenuated due to the post-processing methods, it is recognizable that the smoothed gravity change in Fig. 1(a) is mainly attributed to the co-seismic change, and a smaller post-seismic effect accumulated for a year. The largest gravity drop shown in the East Sea is closely associated with significant crustal density decrease due to dilatation of hangingwall after the rupture (Han et al. 2014). Without any treatment for this signal, it would contaminate surface mass change estimates over the region as shown in Fig. 1(b). Further, high-resolution surface mass signal recovery through forward modelling (FM, Chen et al. 2013) based on this result gives an unreasonable solution of surface mass changes (Fig. 1c). The mass decrease along the coastline is mostly artefacts from contamination of the earthquake-related negative gravity anomaly.

These earthquake signals included in GRACE data can be effectively reduced via earthquake models for surface mass change recovery. For decades, many studies have proposed numerical models for gravity changes due to earthquake events. Co-seismic gravity change associated with elastic deformation of the solid Earth due to faulting was examined in a homogeneous half-space (Okada 1985; Okubo 1992) and a layered spherical Earth (Pollitz 1996). Post-seismic relaxations had been predicted by many studies that paid attention to the role of asthenosphere (Rundle 1982; Piersanti et al. 1995), and Pollitz (1997) introduced a model for the viscoelastic responses in a layered spherical Earth. Recently, the rheology tests consider 3-D space for a more realistic Earth structure (e.g. Sun et al. 2014; Pollitz 2020).

GRACE observations have been also used to understand giant earthquake events. Previous studies investigated earthquake sources of the 2004 Sumatra–Andaman earthquake (Han et al. 2006; de Linage et al. 2009; Broerse et al. 2011), the 2011 Tohoku–Oki earthquake (Han et al. 2011, 2014; Cambiotti & Sabadini 2012; Panet et al. 2018; Cambiotti 2020; Ji et al. 2022) and other significant events (Heki & Matsuo 2010; Han et al. 2015; Xu et al. 2017; Chao & Liau 2019). Further, there are studies considering earthquake signal correction to obtain improved hydrological mass changes. For example, Deggim et al. (2021) proposed a surface mass change data set considering earthquake corrections for 2004 Sumatra–Andaman and 2011 Tohoku–Oki events by using a Bayesian approach (Einarsson et al. 2010). Han et al. (2013) analysed five major earthquake sources by using GRACE observations, and the estimated co-seismic changes were used in the Mascon solutions of surface mass changes (Luthcke et al. 2013; Watkins et al. 2015). However, the model prediction is supposed to be refined by including the geopotential changes due to viscoelastic relaxations. The post-seismic changes were often analysed by empirical modelling (e.g. logarithmic function fit), and recent studies (Han et al. 2016, 2019a) consider viscoelastic models to explain observed post-seismic gravity variation.

Given these pioneering works using numerical models and GRACE observations for giant earthquakes, we see that a self-consistent earthquake correction for GRACE data and the applications are now available. The methodology includes: (1) analysis of giant earthquake sources constrained by GRACE observations, (2) GRACE data correction by using the GRACE-aided elastic and rheological model prediction, (3) signal recovery based on the reduced data and (4) finally, improved surface mass change estimates with the solid Earth effect removed. Thus in this study, we applied these steps (from modelling to application) to the 2011 Tohoku–Oki earthquake signals. For self-consistent geophysical models of the fault parameters and underlying rheological structure, our method considers a fault slip inversion using the GRACE co-seismic changes and tests various viscosity models for the asthenosphere to find the model in a good agreement with the observed post-seismic gravity relaxation patterns. Here, we used numerical models for elastic and viscoelastic responses due to finite-fault model placed in a layered spherical Earth (Pollitz 1996, 1997), and the detailed method is provided in Section 2. After finding the optimal model parameters from GRACE data, the earthquake-induced monthly gravity changes were computed based on the physical earthquake deformation modelling to correct GRACE data. Finally, we applied a signal recovery method to the post-processed GRACE data and examined the impact of our earthquake correction on the estimate of surface mass changes in the East Asia region. In particular, we show the effect of earthquake correction in terms of regional sea level and hydrological mass changes.

2 DATA AND METHODS

2.1 GRACE data

Due to the stripe errors in monthly GRACE L2 data, earthquake signals inferred from full SH coefficients are subject to large systematic errors. These errors can be suppressed by filtering techniques such as Gaussian smoothing and decorrelation filters, but the application would distort earthquake-induced geophysical signals too. Alternatively, we can use a truncated data set excluding the contaminated coefficients at the higher SH degrees and orders. This method is also able to cause a signal loss (omission error) when discarding too many coefficients. Thus, to find the optimal data processing for earthquake signal analysis, we tested various choices by applying the earthquake analysis method that will be elaborated in Section 2.2 below.

In brief, the test showed that the GRACE data truncated up to SH degree 40 (without applying Gaussian smoothing and/or decorrelation filters) provided the preferable results with minimizing possible complexity and uncertainty associated with the post-processing methods (see Figs S1–S4, Supporting Information). Thus, in Sections 2.2 and 2.3 below, we investigate the co- and post-seismic changes of the 2011 Tohoku–Oki earthquake based on GRACE SH coefficients simply truncated at degree 40. Using the earthquake source and rheological parameters to be examined here, we will compute the global monthly gravity responses from both co- and post-seismic model predictions up to degree 40. As detailed in Section 2.4, the model results will be used for the correction of GRACE data for subsequent water mass load recovery.

2.2 Fault slip model from GRACE co-seismic changes

Earthquake slips are fundamental parameters for co-seismic gravity changes with a given elastic Earth structure, and determine the initial stress field to be applied to the underlying viscoelastic asthenosphere (Han et al. 2019a). The slip models are routinely constructed in terms of a composite of multiple finite faults from the seismic data analysis. However, gravity changes and surface deformations computed from such seismic models often show discrepancy from the geodetic observations (Han et al. 2014). It is in part due to differences in temporal and spatial resolutions between geodetic data and seismic data, in addition to model uncertainty. We also identified significant disagreements between co-seismic gravity changes inferred from various seismic models and GRACE observations (e.g. supporting information of Han et al. 2010), suggesting that GRACE data reduction based on those seismic-only models would be inconsistent. Thus, for a self-consistent analysis, we estimated a fault slip model independently based on GRACE observations, rather than simply adopt the slip vectors from seismic models.

Using the CSR RL06 monthly spherical harmonic coefficients, we analysed gravity changes from 2006 to 2016 March, which covers 5 yr before and after the 2011 Tohoku–Oki earthquake. We excluded the monthly data of March 2011 in which the earthquake occurred. The time-variable monthly gravity changes at a resolution up to SH degree 40 at every grid (1°×1° grids) were analysed by using the least-squares fit. We used the following model to describe the monthly gravity change |${\rm{g}}( t )$|⁠:

for the pre-seismic period, and

for the post-seismic period. Here, t indicates the monthly time interval from the earthquake time. |${\rm{B}}$| is a bias, and |${U}_n$| and |${{\rm{V}}}_{\rm{n}}$| are the amplitudes of periodic signals with the frequency |${f}_n$|⁠; in this study, we considered two frequencies of the annual and semi-annual variations. |${{\rm{T}}}_1$| indicates the linear trends before the earthquake, and it should reflect gravity change effects irrelevant to the earthquake event, such as ocean mass changes, etc. Thus, the parameter |${{\rm{T}}}_1$| can be used as the background gravity change at each position. To model earthquake-related signals, we included |${\rm{C}}$| to obtain the co-seismic jump and |${\rm{P}}$| to fit the curve of time-dependent post-seismic relaxation pattern. Contrary to the previous studies considering multiple logarithmic and exponential terms to model the viscoelastic relaxation curve (e.g. Altamimi et al. 2016), we simply applied a single logarithmic term (with |$\tau $| = 6 months) aided by a post-seismic linear trend |${{\rm{T}}}_2$| considering relatively smooth post-seismic variability of the truncated GRACE data. These parameters represent the transient and long-term viscoelastic change, respectively. Note that this empirical parametrization aids only in accurate quantification of the co-seismic change parameter |${\rm{C}}$|⁠. It is possible to consider other fitting models to obtain reasonable co-seismic jumps (see Fig. S5, Supporting Information). As detailed in the next section, these parameters (e.g. |${{\rm{T}}}_2$| and |${\rm{P}}$|⁠) are not directly adopted as the post-seismic changes, and additional model tests were used for the analysis.

The estimated co-seismic jumps based on the parameter |${\rm{C}}$| are mapped in Fig. 2(a). The spatial distribution includes a typical co-seismic gravity change pattern due to a megathrust event, showing the negative pole and positive arc-shape signals (e.g. Cambiotti & Sabadini 2012; Han et al. 2014). On the other hand, gravity decrease captured in the Northeast China and the increase over the East China Sea (roughly from 125° to 130° of longitude) were from unmodelled inter-annual changes in 2011 March, and the signals are possibly related to the hydrological mass variations over the regions rather than the earthquake effect.

Based on the co-seismic gravity changes shown in Fig. 2(a), we computed the fault slip model using a linear inversion. Gravity changes by an arbitrary slip vector can be linearly decomposed into the contributions due to strike-slip and dip-slip components when the fault geometry is fixed. Here we used the geometry parameters (including strike, dip, depth, along-strike length, along-dip length and location) of a finite-fault model defined by Shao et al. (2011). The original model includes hundreds of fault patches, but we reduced the number of slip vectors to five considering the relatively low spatial resolution of GRACE data (∼500 km, Table 1). Then, co-seismic gravity changes due to unit slip vectors (e.g. 1 m of strike-slip and dip-slip motions) on the five fault patches were computed by using the numerical code STATIC1D (Pollitz 1996), and the unit responses were used as the Green’s functions of co-seismic gravity change. The linear inversion combining the Green’s functions with the GRACE observation in Fig. 2(a) consequently yields the slip vector solution. The computed rakes and slips of the five fault patches are summarized in Table 1, and the vector forms are depicted as the thick black lines in Fig. 2(b). Blue rectangles shown in the figure are horizontal projections of the finite-fault planes that we applied. The directions of the slip vectors indicate relative motion of the hangingwall against the footwall, showing oceanward elastic rebound of the hangingwall in general. From the five slip vectors, the average slip length is estimated to be about 9 m, and the average rake is about 86° (Table 1). The seismic moment (M₀) was estimated by 4.53 × 10²² N·m, corresponding to the moment magnitude (M_w) of 9.07, which is in agreement with the previous studies (Hayes 2011; Simons et al. 2011; Cambiotti & Sabadini 2012; Han et al. 2013). The gravity changes from the model predictions (Fig. 2b) generally shows a similar spatial distribution with the observed changes in Fig. 2(a). Further, the maximum gravity drop of −12 μGal in the centre of the negative pole and gravity increase of +6 μGal in the southern part of the positive arc were well estimated by this fault slip model. Here, we show the result that simply used all data points shown in Fig. 2(a) for the inversion. Other results considering a spherical cap to exclude errors on the northwest and southwest corner of the figure (e.g. by using data within a circular area around the epicentre, etc.) also yield similar estimates. It would indicate that the errors associated with unmodelled periodic signals have a relatively small impact on the fault parameter solution probably due to the distances from the epicentre.

2.3 Post-seismic gravity changes and viscosity of asthenosphere

Using the numerical code called VISCO1D (Pollitz 1997), we computed the time-dependent viscoelastic gravity changes followed by the fault slips shown in Fig. 2(b). The most significant variables governing the post-seismic relaxation are the viscosity of the asthenosphere and the thickness of the overlying elastic layer (i.e. lithosphere). A thin elastic layer and a low-viscosity asthenosphere lead to rapid relaxation pattern (Han et al. 2016). When considering acceptable ranges of both variables, the effect of the asthenosphere viscosity is more significant compared to the elastic layer thickness. Thus, we tested various viscosities for the asthenosphere by taking the average elastic layer thickness of 65 km. Here we tested both Maxwell rheology and Burgers rheology with various Kelvin (transient) and Maxwell (steady state) viscosities ranging from 10¹⁸ to 10¹⁹ Pa·s for the asthenosphere from 65 to 220 km depth. We constrained the upper mantle (220–670 km) and the lower mantle (670–2900 km) to have fixed viscosities of 10²⁰ and 10²¹ Pa·s, respectively, and the other parameters like elastic moduli and density were consistent with Preliminary Reference Earth Model (PREM, Dziewonski & Anderson 1981). These factors have negligible impacts on the post-seismic relaxation pattern.

Time-dependent gravity changes from GRACE and the simulations based on viscosity models are compared in Fig. 3. Grey dots represent monthly gravity changes from GRACE SH degree up to 40, and the thick grey lines indicate the least-squares fit of the data (eqs 1 and 2) from 2006 to 2016 March. Several modelling results from different Maxwell (M) and Kelvin (K) viscosity are shown as curves with different colours and markers. To consider possible background gravity changes, these model predictions additionally included the linear trend before March 2011 (i.e. parameter |${{\rm{T}}}_1$| in eq. 1) at each position. The similarity between the observation (thick grey curves) and the predictions (coloured curves) were examined by using the average root mean square (RMS) of the difference at the marker positions (last months of 2011–2015). This averaging considers all data points within an angular distance of 10° from the epicentre (see Figs S6 and S7, Supporting Information). The subpanels of Fig. 3 show the experiments on six important positions, whose locations are presented in the top panel.

Fig. 3 shows that the most similar post-seismic relaxation patterns with the observed gravity changes were obtained from Maxwell body with the viscosity of 5 × 10¹⁸ Pa·s (black curves with squares) on average. Estimates from Maxwell viscosity of 4 × 10¹⁸  and 6 × 10¹⁸ Pa·s (not shown) yielded favourable results as well (Table S1, Supporting Information). This is different from the result of Han et al. (2014), which inferred the asthenosphere to be Burgers body with Kelvin viscosity of 1 × 10¹⁸ Pa·s and Maxwell viscosity of 10 × 10¹⁸ Pa·s (green curves with triangles). Application of the newest version (RL06) of the GRACE data and analysis for the longer data period (∼5 yr after the event) could explain the difference and constrain better the steady-state change. As shown in Fig. 3, both models from the previous study (green) and this study (black) actually provide nearly identical viscoelastic relaxation patterns at the distant locations showing time-dependent gravity decrease due to the converging motion of mantle toward the epicentre (a, b, e, f). On the other hand, the locations c and d above the fault planes show the most significant gravity increase over time, which is associated with viscoelastic uplift of the lithosphere-asthenosphere boundary (Han et al. 2019a). There, both predictions (green and black curves) are similar to one another in terms of the accumulated gravity change until December 2015 (e.g. ∼8 μGal at d). However, future predictions using the bi-viscous (Burgers) model indicated by the green curves with triangles would make larger long-term errors considering the lower trend in the later period. Further, drastic transient gravity changes associated with Kelvin viscosity are not visible in the truncated monthly CSR RL06 solutions as shown in the figure.

An asthenosphere model with Maxwell viscosity of 5 × 10¹⁸ Pa·s provided the post-seismic gravity changes in general agreement with GRACE observations. However, there are also some important deviations. For example, the viscoelastic gravity changes from our model (black curves with squares) are slightly slower than the observed speed (thick grey curves) at locations a and b, and these points are on the continental plate. On the other hand, the predictions are slightly faster on the oceanic plate (locations e and f). This is probably because our simple 1-D modelling does not consider a laterally heterogeneous Earth structure around the subduction zone. This discrepancy indicates that the actual viscosity would be lower than the average value (5 × 10¹⁸ Pa·s) at locations a and b, and even higher at e and f. Indeed, the result of Fig. 3 shows that the lower value of Maxwell viscosity of 3 × 10¹⁸ Pa·s (blue curves with circles) on the continental plate and the larger viscosity of 10¹⁹ Pa·s (red curves with triangles) on the oceanic plate provided the closer changes to the GRACE observations. For comparison, we analysed the post-seismic surface deformation by using GNSS data of the Korean peninsula (see Fig. S8, Supporting Information). The model test showed that GNSS-observed long-term surface motions were generally in agreement with the predictions using Maxwell viscosity of 3 × 10¹⁸ Pa·s, providing a similar conclusion with GRACE data analysis.

In brief, our experiments show that the asthenosphere of Maxwell rheology with the viscosity ranging from 4 to 6 × 10¹⁸ Pa·s is one of the plausible models explaining the GRACE observations over the study area. Thus, we used the model predictions from Maxwell viscosity of 5 × 10¹⁸ Pa·s and the elastic layer thickness of 65 km for the subsequent GRACE data reduction. Further, considering that the lower viscosity (3 × 10¹⁸ Pa·s) provided the closer post-seismic gravity changes and surface motions on the continental plate, we will additionally discuss the effect of the model difference in Section 3 below.

2.4 Earthquake correction for GRACE data

From both results of the slip vector solution on the fixed fault geometry (Fig. 2) and vertically stratified Earth structure with a Maxwell rheology for the asthenosphere (Fig. 3), we computed the global gravity changes due to the 2011 Tohoku–Oki earthquake. The data are provided monthly at every 1°×1° grid. For example, the data are zero everywhere before 2011 March, and it is equal to the co-seismic change computed from our slip model at 2011 March. Monthly data after 2011 April is equal to the sum of the co-seismic change plus the post-seismic change accumulated until that time. These global monthly gravity changes provide the corresponding time-variable Stokes coefficients associated with solid Earth mass redistribution due to the 2011 event. Using these model earthquake signals, we corrected the CSR RL06 solutions from 2003 January to 2016 August, considering most of GRACE mission period. The data sets from 2011 March are reduced with our monthly earthquake model predictions.

As explained above, we considered the truncated GRACE data at SH degree up to 40 to estimate the earthquake model parameters, and we computed the model predictions via normal mode summations up to SH degree 40. Although our model correction does not include the higher GRACE SH coefficients (from degree 41 to 60), the missing terms are expected to have little impact when considering the subsequent post-processing methods (e.g. Gaussian smoothing). As will be described in the next paragraph, we also consider the filtering techniques to reduce the typical high-frequency noises of GRACE data after the earthquake-related signal correction. Due to the significant signal reduction of the higher SH coefficients in the smoothed GRACE data, both corrections using predictions with SH degrees up to 40 and 60 provide nearly identical results to one another in the smoothed forms (see Fig. S9, Supporting Information). Thus, we simply applied the monthly earthquake model predictions up to SH degree 40 for the GRACE data reduction.

After reducing the earthquake-related effects from GRACE data, we followed ordinary post-processing steps; the GIA effect is corrected by ICE6G-D model (Peltier et al. 2018), and SH degree-1 and -2 coefficients are adjusted by using the alternative estimates (Kim et al. 2019). Finally, to suppress the correlated (stripe-shape) errors in the data, we applied both the 500 km Gaussian smoothing (Wahr et al. 1998) and the P4M10 polynomial fit method (Swenson & Wahr 2006).

To overcome the low spatial resolution of the smoothed GRACE data, we applied a signal recovery process called FM method (Chen et al. 2013), and obtained high-resolution surface mass changes. This method finds surface mass change estimates with reduced signal leakages by iteratively updating the initial estimate. In addition, for more realistic ocean mass changes of the solution, we included sea level fingerprint (SLF) computation dependent on each land mass change distribution in the process (Jeon et al. 2021). As previously shown in Fig. 1(c), this FM method generated strange signals on the adjacent coastal region because of the uncorrected earthquake-related solid Earth signals in the GRACE data. In the next section, we present another FM result aided by our correction for the solid Earth signals due to the 2011 Tohoku–Oki earthquake.

3 RESULT

We demonstrate the effect of our earthquake correction on surface mass change estimates for the study period from 2003 January to 2016 August in Fig. 4. Fig. 4(a) is the improved version of Fig. 1(c), which is the result after removing the solid Earth mass redistribution by the 2011 Tohoku–Oki earthquake. This 1-yr mean difference shows that our correction greatly suppresses peculiar errors along the coast of the East Sea, previously found in Fig. 1(c). The following subpanels show the average time-dependent surface mass changes over the East Sea (b) and the northern and southern coastal lands (c and d, respectively). In Figs 4(b)–(d), the grey solid curves are based on the original data contaminated by the earthquake signals. The black solid curves represent the surface mass changes with our earthquake model correction, which use the co-seismic changes from Fig. 2(b) and the post-seismic relaxation due to Maxwell viscosity of 5 × 10¹⁸ Pa·s for the asthenosphere. For comparison, we additionally included the results considering the smaller viscosity (Maxwell viscosity of 3 × 10¹⁸ Pa·s) as identified by both GRACE and GNSS observations on the continental plate (black dotted curves).

In Fig. 4(b), we presented the average ocean mass changes over the East Sea in terms of SLF computation based on our FM results (Jeon et al. 2021). As shown in the previous figure, the East Sea suffers extreme gravity variations due to co- and post-seismic changes since the 2011 event, and the contribution on average appears as sea level drop. Without correction of the solid Earth mass signals, the average (mass-only) sea level change yields a lower estimate (1.59 ± 0.11 mm yr⁻¹, grey curve) due to the earthquake signals propagated into the FM result as previously shown in Fig. 1(c). In particular, being combined with 2011 La Ninã event (Boening et al. 2012; Cazenave et al. 2014), the earthquake-induced errors create more significant negative anomaly in the East Sea estimate in 2011. With our solid Earth mass corrections, the resultant ocean mass change estimate becomes more consistent with the global average mostly affected by climate change. Indeed, our earthquake correction yields ocean mass change of the East Sea (black curve) with the inter-annual changes in agreement with the global mean (blue curve with an offset, 2.04 ± 0.10 mm yr⁻¹). The climate-related regional mass rate is refined to be 1.87 ± 0.11 mm yr⁻¹ (black solid curve, considering Maxwell viscosity of 5 × 10¹⁸ Pa·s) or 1.86 ± 0.11 mm yr⁻¹ (black dotted curve, 3 × 10¹⁸ Pa·s), which is almost 20 per cent larger than the previous estimate including a bias due to the 2011 earthquake effect. On the other hand, the difference between two different viscosity model corrections (i.e. solid and dotted black curves) had an insignificant impact on the SLF-based ocean mass changes as shown in the figure.

Figs 4(c) and (d) show surface mass changes of our FM results averaged over the northern and southern coastal regions of the East Sea, respectively. The northern coastal region covers most of Korea, and the eastern edge of China and Russia; a dashed line in Fig. 4(a) indicates the area boundary. This area is under the signal leakage of significant gravity decrease by the co- and post-seismic gravity changes. The solid mass adjustment by the earthquake results in the biased surface mass change estimate after 2011 March (grey curve in Fig. 4c), and the estimate is effectively corrected via our earthquake modelling (black solid curves). This estimate is further improved when being aided by a model with the lower-viscosity asthenosphere (3 × 10¹⁸ Pa·s, black dotted curves) inferred from data on the continental plate. Nevertheless, the difference between both corrections is relatively small considering the data variability.

For the southern coastal region (Fig. 4d), we simply investigated average surface mass changes over mainland Japan (i.e. Honshu island). Particularly in the eastern part of the region, there has been the significant co-seismic gravity drop followed by gradual gravity increase due to the post-seismic rebound, and the temporal patterns are also reflected in the uncorrected time-series of FM solution (grey curve). Our earthquake correction with Maxwell viscosity of 5 × 10¹⁸ Pa·s (black solid curve) largely reduces such signals over mainland Japan, resulting in a slight decreasing surface mass trend throughout the study period. In contrast, the estimate corrected by the lower-viscosity model led to a significant negative trend for the post-seismic period (black dotted curves). It means that the model probably overcompensated the post-seismic relaxation effect, and the average viscosity, 5 × 10¹⁸ Pa·s, is more appropriate for the signal correction over this volcanic arc region.

Fig. 4 shows that our co- and post-seismic gravity change model effectively corrects the earthquake-induced solid Earth mass signals from GRACE data. The post-seismic changes due to Maxwell viscosity of 5 × 10¹⁸ Pa·s provide an effective correction in general, and partially the lower viscosity can be considered for the landward locations. Some data points particularly in the South Kyushu (red speckles around 32° N and 130° E in Fig. 4a) seem to have short-term errors. This may indicate the uncertainty of our earthquake modelling and/or GRACE observations, and it will be discussed in the next section.

4 DISCUSSION AND CONCLUSION

This study developed the fault slip model and asthenosphere rheological parameters of the 2011 Tohoku–Oki earthquake by using the monthly GRACE gravity field observations. The constrained physical models of the solid Earth adjustment to the earthquake were then used to make the GRACE data free of the solid Earth mass change signals, and as a result, we demonstrated the improved surface mass change recovery in the East Asia region.

The fault slip models were obtained by using a linear inversion based on the observed GRACE co-seismic gravity changes and the (elastic) Green’s functions computed from a priori fixed fault geometry. Here, we estimated the co-seismic changes with the parameters including the Heaviside step function from GRACE time-series with SH degree up to 40. The uncertainties can arise due to residual stripe errors retained in the lower SH degrees, errors from ocean dynamics and atmospheric pressure model, and any irregular gravity changes that overlap with the earthquake time period. Uncertainties of our co-seismic jumps (denoted by e) were estimated by ∼2 μGal (95 per cent confidence level) over the study area, and we found that the spatial distribution is highly correlated with the locations of the residual stripes and unmodelled interannual changes in the truncated GRACE data (see Fig. S10, Supporting Information). Since it was very difficult to analyse the contributions of these error sources individually, we speculated the comprehensive uncertainties due to these error sources by using the extreme values of co-seismic change parameters (i.e. |${\rm{C}} - e$| and |${\rm{C}} + e$|⁠). These extreme values consequently led to different fault model solutions with M₀ of 5.08 × 10²² and 4.11 × 10²² N·m, or M_w of 9.10 and 9.04, respectively. These would represent the uncertainty range of our fault model solution due to possible errors, but they are all within the range of the previous estimates (Hayes 2011; Simons et al. 2011; Cambiotti & Sabadini 2012; Han et al. 2013).

Rheological modelling suggests that the asthenosphere with Maxwell viscosity of about 5 × 10¹⁸ Pa·s and the vertical distribution from 65 to 220 km in depth provides the best fit to GRACE observations in the study area. However, this configuration is based on 1-D Earth structure, which does not consider 3-D heterogeneity near the subduction zone, such as presence of elastic slab, structural difference beneath oceanic and continental plates and a low-viscosity zone below the volcanic region (Ueda et al. 2003; Sun & Wang 2015). Those effects resulted in the slight positive and negative errors on the continental plate and the oceanic plate, respectively (Fig. 3). This discrepancy can be explained by the heterogeneity of viscosity that may be lower under the continental plate and higher under the oceanic plate than the average value of 5 × 10¹⁸ Pa·s that we used. Indeed, both GRACE and GNSS data on the continental plate showed the better agreement between the observations and predictions when using an asthenosphere viscosity of about 3 × 10¹⁸ Pa·s. However, considering the facts that the lithosphere of continental plate is possibly thicker than our assumption (65 km) and that a thick elastic layer mitigates the viscoelastic relaxation of the asthenosphere, the actual viscosity under the continental margin may be even lower than our estimate. Further, when including the effect of subducting slab, the magnitudes of post-seismic changes based on 1-D model (e.g. top panel of Fig. 3) can be decreased by 10–20 per cent around the East Sea (Tanaka et al. 2015); this may also lead to the lower viscosity estimates. We expect that such estimations associated with structural heterogeneity would be greatly improved via 3-D modelling (e.g. Sun et al. 2014; Tanaka et al. 2015).

Based on the estimated model parameters, we included the earthquake correction (the solid Earth mass redistribution) in ordinary GRACE post-processing steps. Surface mass changes are computed from geopotential changes measured by GRACE assuming that all changes take place at the surface of the Earth, and thus important non-surficial solid Earth signals like earthquakes as well as GIA must be removed a priori. In order to recover the smoothed signals (leakage effect), we examined the surface mass changes after applying a signal recovery method that we refer to as FM (Chen et al. 2013). There, our earthquake correction significantly reduced errors associated with the solid Earth deformation in both land and ocean mass changes (Fig. 4). The residual earthquake signal accounts for about 0.3 mm yr⁻¹ in terms of SLFs over the East Sea, which is equivalent to 20 per cent increase compared to the estimate without the earthquake correction. The terrestrial water and ice mass changes also showed a great improvement, reducing step-like patterns from the earthquake event from the average time-series over the adjacent coastal regions. There are some local-scale errors as shown in the map of Fig. 4(a), and they are likely associated with uncertainties of GRACE data fit, unmodelled periodic and episodic signals, limitations of 1-D earthquake modelling, and/or FM computation (Fig. S9, Supporting Information).

Our result shows both the significant impact of earthquake deformation on regional-scale surface mass change estimates and the possibility of the GRACE data reduction using physical models for earthquake events. Besides the 2011 Tohoku–Oki earthquake examined in this study, there are many earthquakes that have been detected by GRACE (e.g. de Viron et al. 2008; Chao & Liau 2019). For example, the 2004 Sumatra–Andaman (M_w ∼9.3) and 2010 Maule (M_w ∼8.8) events show significant broad-scale impacts as well, and those signals can be corrected by the similar method used in this study. Combined contributions of the individual seismic events will provide the global gravity perturbations due to earthquakes occurred in GRACE period. As discussed in presentations of Han et al. (2019b) and Rietbroek et al. (2020), such a comprehensive product will be an useful tool for users who want to remove earthquake-related signals from GRACE data for improved climate analyses.

SUPPORTING INFORMATION

Figure S1. Left-hand column: observed co-seismic gravity changes from GRACE data with different SH degree truncation. Middle column: fault slip vectors and the co-seismic gravity change predictions, computed based on the GRACE observations on the left by using a linear inversion. Right-hand column: observations minus predictions, divided by maximum gravity drop of the prediction, expressed as a percentile. These results are based on CSR RL06 GRACE data.

Figure S2. Equal to Fig. S1, except that these are based on GFZ RL06 GRACE data.

Figure S3. Equal to Fig. S1, except that these are based on JPL RL06 GRACE data.

Figure S4. Similar to Fig. S1, except that these results consider different Gaussian smoothing range and application of decorrelation filter applied to CSR RL06 GRACE data.

Figure S5. Three different least-squares fits from a GRACE gravity time-series. Grey dots indicate monthly GRACE observations at 38° N and 142° E, with summations of SH degree up to 40. Grey lines in subpanels (a)–(c) represent different fittings as described in the text above. Annual and semi-annual changes are removed from both GRACE data and the fittings.

Figure S6. Accumulated post-seismic gravity change in 2015 December. (a) Post-seismic gravity change inferred from the least-squares fit of GRACE data. (b) Post-seismic gravity change from asthenosphere model with Maxwell viscosity of 5 × 10¹⁸ Pa·s. Circles in both subpanels indicates an angular distance of 10° from the epicentre.

Figure S7. Spatial difference of post-seismic gravity changes from a viscosity model and GRACE observation with respect to time. Difference between Figs S6(a) and (b) is equal to the distribution of (e). RMS of the difference and the mean value is provided in the figure. Here, the model prediction considers Maxwell viscosity of 5 × 10¹⁸ Pa·s, and it is the result showing the lowest mean RMS difference among the cases we tested.

Figure S8. GNSS records of Korea. Top: locations of GNSS stations. (a)–(f) Daily GNSS observations (grey dots) and the least-squares fit (grey curves) are compared with three different estimates from viscosity models. The blue dotted curves with circles use Maxwell body with the viscosity of 3 × 10¹⁸ Pa·s. The black solid curves with squares are from Burgers body with Kelvin viscosity of 1 × 10¹⁸ Pa·s and Maxwell viscosity of 3 × 10¹⁸ Pa·s. The yellow dotted curves with triangles are computed from Maxwell body with the viscosity of 5 × 10¹⁸ Pa·s. The elastic layer thickness is fixed to be 65 km for all cases. The eastward components are shown, and the background trend associated with the plate motion is removed by using the pre-seismic trend.

Figure S9. The co-seismic changes of the 2011 Tohoku–Oki earthquake due to our model parameters with different resolutions and the smoothed form. (a) The co-seismic gravity change prediction with a summation of SH degree up to 40. (b) The co-seismic gravity change prediction with a summation of SH degree up to 60. (c) and (d) The smoothed forms of (a) and (b), respectively, by using 500 km Gaussian smoothing.

Figure S10. (a) Uncertainty distribution of Fig. 2(a): co-seismic changes inferred from parameters of Heaviside step function. The uncertainties are all examined with 95 per cent (2σ) confidence levels. It additionally includes locations of time-series in subpanels (c) and (d). (b) Spatial distribution of gravity disturbance in 2011 March, considering SH degree up to 40. (c) Time-variable gravity changes at the position of 48° N and 131° E. The black curve with dots indicates monthly GRACE data with SH degree up to 40, and the thick grey curve shows the least-squares fit for ± 5 yr from the 2011 event, plotted without the annual and semi-annual variations. (d) Equal to (c), but examined at 30° N and 128° E.

Table S1. Mean RMS differences of observation and prediction from different viscosity models. Ten cases with different Kelvin viscosity (K) and Maxwell viscosity (M) are listed here.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We deeply appreciate the editor Prof Bert Vermeersen and two anonymous reviewers who helped us to improve this manuscript with constructive and important comments. We also thank Jeanne Sauber who reviewed and improved the original manuscript. This work was supported by National Research Foundation of Korea (NRF) grant (2022R1C1C2006586). K-WS was supported by NRF grant (2023R1A2C1004899) and the Korea Institute of Marine Science and Technology Promotion (KIMST) grant funded by the Ministry of Ocean Fisheries (KIMST 20190361). S-CH was supported by Australian Research Council Discovery Program (DP170100224) and NASA GRACE/GRACE Follow-On Science Team Project.

DATA AVAILABILITY

CSR RL06 GRACE data sets are available at GRACE Tellus website (https://grace.jpl.nasa.gov/data/get-data), and GNSS records are from Nevada Geodetic Laboratory (http://geodesy.unr.edu/index.php). The monthly gravity changes due to co- and post-seismic effect of the 2011 Tohoku–Oki event defined in this study are available at Open Science Framework (https://doi.org/10.17605/OSF.IO/8ZYU9). Any other data shown in this paper can be shared on request to the corresponding author.

References