Seismic Impact of Large Earthquakes on Estimating Global Mean Ocean Mass Change from GRACE

: We analyze the impact of large earthquakes on the estimation of the global mean ocean mass (GMOM) change rate over the 13-year period (January 2003 to December 2015) using the Gravity Recovery and Climate Experiment (GRACE) Release-06 (RL06) monthly gravity solutions released by the Center for Space Research (CSR). We take into account the e ﬀ ects of the December 2004 Mw9.1 and April 2012 Mw8.6 Sumatra earthquakes, the March 2011 Mw9.0 Tohoku-Oki earthquake, and the February 2010 Mw8.8 Chile earthquake. After removing the co- and post-seismic e ﬀ ects of these earthquakes in the oceanic areas by least squares ﬁtting, we estimate the GMOM rate from GRACE monthly observations. Results show that GRACE-observed GMOM rate before the seismic correction is 2.12 ± 0.30 mm / year, while after correction the rate is 2.05 ± 0.30 mm / year. Even though the − 0.07 ± 0.02 mm / year seismic inﬂuence on GRACE GMOM rate is small on a global scale, it is a systematic bias and should be considered for improved quantiﬁcation and understanding of the global sea level change.


Introduction
Ocean mass variation is caused by water exchange between the ocean and other components of the climate system, including polar ice sheets, mountain glaciers, terrestrial water storage (TWS), and atmospheric water vapor [1]. Due to the small contributions of atmospheric water vapor, most ocean mass change studies focused on contributions of polar ice sheets, mountain glaciers, and TWS change. For example, the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE) team [2] shows that the contribution of the Antarctica ice sheet melting to the mean sea level rise is 7.6 ± 3.9 mm during the period 1992-2017. In another recent study the IMBIE team [3] indicates that the Greenland ice sheet melting contributes 10.6 ± 0.9 mm to the mean sea level rise for the period 1992-2018. Wouters et al. [4] demonstrates that the contribution of global glaciers excluding Greenland and Antarctica melting to the mean sea level rise is 8.0 ± 1.3 mm during the period 2002-2016. Reager et al. [5] indicates that climate change has resulted in additional storage of TWS, slowing down the ocean mass change by 0.71 ± 0.20 mm/year during the period 2002-2014.
During the GRACE mission period (2002-2017) several large earthquakes caused observable gravity change signals, which were detected by GRACE [20,21]. The co-and post-seismic effects have been well investigated with GRACE measurements combined with other observations such as Global Navigation Satellite System (GNSS) displacements, for the events including the 2004 Sumatra Mw9.1 [22], 2011 Tohoku-Oki Mw9.0 [23,24], 2010 Chile Mw8.8 [25], 2005 Nias (Sumatra) Mw8.6 [21,22], 2012 Sumatra Mw8.6 [26], 2007 Bengkulu (Sumatra) Mw8.5 [20,21], 2006 Kuril Islands Mw8.3 [27,28], and 2013 Okhotsk Sea Mw8.3 [29] earthquakes. Among these large earthquakes some have induced noticeable mass change signals, which are clearly captured in the global mass change (expressed in equivalent water height, EWH) trends from GRACE ( Figure 1a). Figure 1a shows that the seismic signals within the three regions marked by the red circles are apparent, revealing clear positive-negative spatial patterns corresponding to mass increase and loss induced by the earthquakes. Figure 1b shows the example of the GRACE mass change time series at the point A (2.5°N, 93.5°E) from January 2003 to December 2015, which is located at the oceanic area within the Sumatra earthquake region. The time series exhibit a clear jump in December 2004 (i.e., the occurrence month of the megathrust earthquake) as well as gradual increases thereafter, which are clearly caused by co-and post-seismic deformations, rather than sea water mass changes. Therefore, we need to correct these seismic effects to avoid potential systematic bias in GRACE GMOM rate estimation. The main purpose of this study is to estimate the effect of large offshore earthquakes on the GMOM change rate 2003-2015 observed by GRACE. The earthquakes we select include the December 2004 and April 2012 Sumatra earthquakes, March 2011 Tohoku-Oki earthquake, and February 2010 Chile earthquake, whose effects are visible in Figure 1a. We correct the seismic effect using least squares fitting considering both co-and post-seismic terms. Then we compare the GMOM change rates before and after the correction, in order to show how significantly these large earthquakes affect the determination of GRACE GMOM rate.

Data and Method
We used the Center for Space Research (CSR) Release-06 (RL06) monthly gravity field solutions to estimate the ocean mass change rate during the period January 2003 to December 2015. The The main purpose of this study is to estimate the effect of large offshore earthquakes on the GMOM change rate 2003-2015 observed by GRACE. The earthquakes we select include the December 2004 and April 2012 Sumatra earthquakes, March 2011 Tohoku-Oki earthquake, and February 2010 Chile earthquake, whose effects are visible in Figure 1a. We correct the seismic effect using least squares fitting considering both co-and post-seismic terms. Then we compare the GMOM change rates before and after the correction, in order to show how significantly these large earthquakes affect the determination of GRACE GMOM rate.

Data and Method
We used the Center for Space Research (CSR) Release-06 (RL06) monthly gravity field solutions to estimate the ocean mass change rate during the period January 2003 to December 2015. The solutions contain spherical harmonic (SH) coefficients up to degree and order 60 [30]. We applied the 300 km Gaussian smoothing to suppress the spatial noise of the GRACE observations [31]. Solid earth deformation effect due to the glacial isostatic adjustment (GIA) was removed using the model from A et al. [16], thereafter referred as the A13 GIA model. We replace the GRACE C 20 terms with the Satellite Laser Ranging (SLR) observations [32]. Since the degree-1 terms, representing the geocenter motion were absent in the GRACE monthly solutions, we used the degree-1 coefficients provided by the GRACE project (GRACE Technical Note 13) [33] to include the geocenter motion effect, which are derived from GRACE observations using an improved method [12,34]. Then we calculated the global gridded (1 • × 1 • ) surface mass change fields (with respect to the mean field 2003.01-2015.12) from the 142 monthly solutions. By least squares fitting of the time series at each grid point of these monthly solutions, we obtained the global mass change rate distribution.
During the least squares fitting, the seismic effects of the three regions (as circled by red curves in Figure 1) need to be corrected. The correction is carried out by removing the co-and post-seismic terms with least squares fitting of the GRACE observation time series (e.g., [35,36]). When the time series at the point includes the seismic effect of only one earthquake, the fitting equation can be expressed as: where C 1 and C 2 are the constant terms before and after the earthquake. p is the total post-seismic gravity change reached at the end of relaxation. τ is the e-folding relaxation time of the post-seismic changes. Here we chose τ = 7 months in the Sumatra and Chile earthquake regions fitting. We also chose τ = 4 months in the Tohoku-Oki earthquake region fitting (the selection of the parameter τ on the time series fitting will be discussed later in Section 4). t eq is the earthquake occurrence time. Figure 2a shows an example of the time series fitting at a selected point B (36.5 • N, 145.5 • E) in the Tohoku-Oki earthquake region. When the time series at the point includes the seismic effects of two earthquakes, the fitting equation can be expressed as: where C 1 is the pre-seismic constant of the first earthquake, C 2 is the post-seismic constant of the first earthquake, and C 3 is the post-seismic constant of the second earthquake. t eq1 and t eq2 represent the occurrence times of the first and second earthquakes, respectively. τ 1 and τ 2 are the e-folding relaxation time parameters of the post-seismic change corresponding to the first and second earthquakes, respectively. Note here it is implicitly assumed that the post-seismic deformation of the first event is fully relaxed before the second event occurs. This assumption should make sense for the case of Sumatra earthquake region because the interval between the two events is more than 7 years, during which most post-seismic deformation of the 2004 Sumatra earthquake has been relaxed. Figure 2b shows an example of the time series fitting at a selected point A (2.5 • N, 93.5 • E) in the Sumatra earthquake region. Then we removed the seismic effect from the GRACE observation and obtained the residual time series (GRACE observation minus the co-and post-seismic terms in Equation (1) or Equation (2)) for each grid point over the region. The corrected mass change rate can be figured out from the residual time series with least squares fitting, expressed as: where A c an , A s an , A c se and A s se represent the periodic terms including the annual and semi-annual variations. A tr represents the corrected linear trend and A 0 is the constant term. By removing the seismic terms from the GRACE time series, we corrected the biases from the impact of seismic deformation.
After correcting the seismic effect, we used the global EWH rates to calculate the mean ocean mass change rate. During the calculation, we need to correct the signal leakage effect due to limited spatial resolution of GRACE observations. Since the GRACE monthly solutions are truncated to the SH degree 60 and the spatial smoothing is also applied to suppress noise, significant mass-change signals leak from land to ocean in the coastal regions (e.g., the Gulf of Alaska, and the surrounding regions of Greenland and Antarctica in Figure 1a). In this study we used the forward modeling (FM) to correct the leakage in the continent-ocean boundary areas. The FM method uses iterations to correct the signal leakage between land and ocean [37]. We calculated the GMOM change rates before and after the seismic correction, and evaluated the bias due to the effects of the large earthquakes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 15 After correcting the seismic effect, we used the global EWH rates to calculate the mean ocean mass change rate. During the calculation, we need to correct the signal leakage effect due to limited spatial resolution of GRACE observations. Since the GRACE monthly solutions are truncated to the SH degree 60 and the spatial smoothing is also applied to suppress noise, significant mass-change signals leak from land to ocean in the coastal regions (e.g., the Gulf of Alaska, and the surrounding regions of Greenland and Antarctica in Figure 1a). In this study we used the forward modeling (FM) to correct the leakage in the continent-ocean boundary areas. The FM method uses iterations to correct the signal leakage between land and ocean [37]. We calculated the GMOM change rates before and after the seismic correction, and evaluated the bias due to the effects of the large earthquakes.

Seismic Effect Correction
First, we corrected the seismic effects in the Sumatra earthquake region. Since we considered two earthquakes in this region (the December 2004 Mw9.1 and April 2012 Mw8.6 earthquakes), we adopted Equation (2) (see Section 2) to correct the seismic effect. We noted that the seismic signals of the March 2005 MW 8.6 Nias earthquake were also included in GRACE observations over this region. Nevertheless, its occurrence time was close to the December 2004, so that we implicitly treated the two earthquakes as one event in the time series fitting. The correction region ranged from -13°N to 17°N in latitude and from 81°E to 111°E in longitude. The results are shown in Figure 3. The corrected spatial distribution (Figure 3b) shows that most seismic effects were removed in this region, and the amplitude of the rate from residual variations was basically consistent with that from the surrounding oceanic areas. Figure 3c

Seismic Effect Correction
First, we corrected the seismic effects in the Sumatra earthquake region. Since we considered two earthquakes in this region (the December 2004 Mw9.1 and April 2012 Mw8.6 earthquakes), we adopted Equation (2) (see Section 2) to correct the seismic effect. We noted that the seismic signals of the March 2005 MW 8.6 Nias earthquake were also included in GRACE observations over this region. Nevertheless, its occurrence time was close to the December 2004, so that we implicitly treated the two earthquakes as one event in the time series fitting. The correction region ranged from −13 • N to 17 • N in latitude and from 81 • E to 111 • E in longitude. The results are shown in Figure 3. The corrected spatial distribution (Figure 3b) shows that most seismic effects were removed in this region, and the amplitude of the rate from residual variations was basically consistent with that from the surrounding The linear rates at points D and E were 3.58 cm/year and −1.97 cm/year before the correction, and were 0.28 cm/year and −0.19 cm/year after the correction, respectively. Both examples show significant changes before and after the correction, since the rate values were dramatically reduced. After seismic correction, the linear rates (0.28 cm/year and −0.19 cm/year) were comparable to the amplitudes of the rates at the other grid points in the surrounding oceanic areas, in which all the grid values were within the range −0.3 to +0.3 cm/year. Figure 3b shows that the corrected linear rates at the grid points surrounding the epicenter were consistently small (as those at the grid points in further oceanic areas), suggesting that the seismic impact was well corrected.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 15 the grid values were within the range -0.3 to +0.3 cm/year. Figure 3b shows that the corrected linear rates at the grid points surrounding the epicenter were consistently small (as those at the grid points in further oceanic areas), suggesting that the seismic impact was well corrected. Second, we corrected the seismic effects in the Tohoku-Oki earthquake region. We adopted Equation (1) (see Section 2) to correct the seismic effect in this region. The correction region ranges from 30°N to 42°N in latitude and from 132°E to 152°E in longitude. Figure 4a,b shows the EWH rate spatial distributions before and after the seismic correction. Similar to that in the Sumatra earthquake region, the time series fitting examples at both the points F (36.5°N, 144.5°E) and G (39.5°N, 138.5°E) in Figure 4c,d, also shows significant changes before and after seismic correction (1.82 vs. 0.09 cm/year for the point F, and -1.59 vs. 0.09 cm/year for the point G). The corrected linear rates surrounding the epicenter were small (within the range -0.1 to +0.1 cm/year, Figure 4b), which indicates that the seismic impact was well corrected for the Tohoku-Oki earthquake region. It should be noted that a discontinuity feature was visible at the correction region boundaries, which will be discussed later in Section 4.2.  Second, we corrected the seismic effects in the Tohoku-Oki earthquake region. We adopted Equation (1) (see Section 2) to correct the seismic effect in this region. The correction region ranges from 30 • N to 42 • N in latitude and from 132 • E to 152 • E in longitude. Figure 4a,b shows the EWH rate spatial distributions before and after the seismic correction. Similar to that in the Sumatra earthquake region  Figure 4b), which indicates that the seismic impact was well corrected for the Tohoku-Oki earthquake region. It should be noted that a discontinuity feature was visible at the correction region boundaries, which will be discussed later in Section 4.2.   Finally, we corrected the seismic effects in the Chile earthquake region. Similar to that in the Tohoku-Oki earthquake region, we used Equation (1) (see Section 2) to correct the seismic effect in this region. The correction region ranged from -41°N to -31°N in latitude and from 277°E to 297°E in longitude. Results show that the changes before and after seismic correction were noticeable, not only in the EWH rate distribution (Figure 5a vs. Figure 5b), but also in the time series fitting  Finally, we corrected the seismic effects in the Chile earthquake region. Similar to that in the Tohoku-Oki earthquake region, we used Equation (1)   Finally, we corrected the seismic effects in the Chile earthquake region. Similar to that in the Tohoku-Oki earthquake region, we used Equation (1) (see Section 2) to correct the seismic effect in this region. The correction region ranged from -41°N to -31°N in latitude and from 277°E to 297°E in longitude. Results show that the changes before and after seismic correction were noticeable, not

Impact of Seismic Effect on the GMOM Change Rate
After correcting the seismic effects in the three regions, we calculated the impact of seismic effect on GRACE GMOM rate. We first estimated the GRACE GMOM rate without seismic effect correction. Based on FM with 100 iterations, the recovered GMOM change rates turn out to be convergent (see the convergence curves with respect to the iterations in Figure 6). Then we took the convergent value at 100 iterations as the leakage-corrected GMOM change rate, which was 2.12 mm/year (shown by the red curve in Figure 6).  Figure 5c vs. Figure 5d). The fitted linear rates at points H and I were 0.89 and -1.78 cm/year before the correction, and 0.09 and 0.01 cm/year after the correction, respectively. In the areas surrounding the epicenter, the corrected linear rates (Figure 5b) were within the range -0.1 to +0.1 cm/year, indicating a good seismic impact correction for the Chile earthquake region.

Impact of Seismic Effect on the GMOM Change Rate
After correcting the seismic effects in the three regions, we calculated the impact of seismic effect on GRACE GMOM rate. We first estimated the GRACE GMOM rate without seismic effect correction. Based on FM with 100 iterations, the recovered GMOM change rates turn out to be convergent (see the convergence curves with respect to the iterations in Figure 6). Then we took the convergent value at 100 iterations as the leakage-corrected GMOM change rate, which was 2.12 mm/year (shown by the red curve in Figure 6). We then estimated the GRACE GMOM rate with seismic effect corrections in three cases. When only correcting the seismic effect in the Sumatra earthquake region, the estimated GRACE GMOM rate was 2.08 mm/year (see the blue curve in Figure 6). The spatial distribution of the global EWH rate after correcting the seismic effect in Sumatra region is shown in Figure 7a. When correcting the seismic effects in the Sumatra and Tohoku-Oki earthquake regions, the GMOM rate was 2.07 mm/year (see the magenta curve in Figure 6). The global EWH rate distribution after correcting the seismic effects in these two regions is shown in Figure 7b. When we corrected the seismic effects in all the three regions (considering the Sumatra, Tohoku-Oki and Chile earthquakes), the estimated GMOM rate was 2.05 mm/year (see the black curve in Figure 6). The seismic-corrected global EWH rate distribution is shown in Figure 7c. The corresponding GMOM change rates are summarized in Table 1.  We then estimated the GRACE GMOM rate with seismic effect corrections in three cases. When only correcting the seismic effect in the Sumatra earthquake region, the estimated GRACE GMOM rate was 2.08 mm/year (see the blue curve in Figure 6). The spatial distribution of the global EWH rate after correcting the seismic effect in Sumatra region is shown in Figure 7a. When correcting the seismic effects in the Sumatra and Tohoku-Oki earthquake regions, the GMOM rate was 2.07 mm/year (see the magenta curve in Figure 6). The global EWH rate distribution after correcting the seismic effects in these two regions is shown in Figure 7b. When we corrected the seismic effects in all the three regions (considering the Sumatra, Tohoku-Oki and Chile earthquakes), the estimated GMOM rate was 2.05 mm/year (see the black curve in Figure 6). The seismic-corrected global EWH rate distribution is shown in Figure 7c. The corresponding GMOM change rates are summarized in Table 1.   As a result, the influence of these four large earthquakes on the estimation of GMOM change rate was -0.07 mm/year. Among them, the bias caused by the Sumatra earthquake region

Uncertainty Estimate of GMOM Change Rate
To estimate the uncertainty of the GMOM change rate obtained in this study, we followed the method used in Chen et al. [9,14] based on least squares fitting. We applied the forward modeling (FM) leakage correction (with 100 iterations, as done for the global trend rate field in Section 3.2) to the monthly global EWH change fields (142 months in total, as stated in Section 2), and obtained the GMOM change time series from January 2003 to December 2015. Then we fitted the time series with least squares by taking account of the constant term, linear trend, annual, and semi-annual terms (see Equation (3) in Section 2). We took the 2 sigma of the fitted trend term as the uncertainty of the GMOM rate, which was ± 0.11 mm/year (see Line 2, Column 2 in Table 2). In addition, by correcting  As a result, the influence of these four large earthquakes on the estimation of GMOM change rate was −0.07 mm/year. Among them, the bias caused by the Sumatra earthquake region (including the

Uncertainty Estimate of GMOM Change Rate
To estimate the uncertainty of the GMOM change rate obtained in this study, we followed the method used in Chen et al. [9,14] based on least squares fitting. We applied the forward modeling (FM) leakage correction (with 100 iterations, as done for the global trend rate field in Section 3.2) to the monthly global EWH change fields (142 months in total, as stated in Section 2), and obtained the GMOM change time series from January 2003 to December 2015. Then we fitted the time series with least squares by taking account of the constant term, linear trend, annual, and semi-annual terms (see Equation (3) in Section 2). We took the 2 sigma of the fitted trend term as the uncertainty of the GMOM rate, which was ± 0.11 mm/year (see Line 2, Column 2 in Table 2). In addition, by correcting the seismic effects in the three earthquake regions (with Equations (1) and (2)) from the monthly global EWH change fields and then applying the FM leakage correction, we obtained the seismic-corrected GMOM change time series. The corresponding fitted GMOM rate with seismic correction and its uncertainty are listed in Table 2 (Line 2, Column 3). We also calculated the GMOM rates with the Geoforschungszentrum (GFZ) RL06 [38] and Jet Propulsion Laboratory (JPL) RL06 [39] GRACE monthly solutions from January 2003 to December 2015. All the data processing steps were the same as those done for the CSR RL06 solutions, except that the degree-1 coefficients corresponding to the GFZ and JPL solutions [32] were used, respectively. The GFZ and JPL GMOM rates without and with seismic correction are listed in Table 2 (Lines 3 and 4), together with their uncertainties. The differences between each other among the GMOM rates from CSR, GFZ, and JPL solutions reached up to 0.11 mm/year (i.e., the difference between the GFZ and JPL rates before seismic correction), which was right within the estimated uncertainty range. Our estimated uncertainty ± 0.11 mm/year was close to the one provided in Chen et al. [14], in which they obtained a ± 0.14 mm/year uncertainty level (see their Table 1) by fitting the GMOM time series from January 2005 to December 2015. Our uncertainty estimate was slightly smaller because we used a longer time span (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) of GRACE data in the time series fitting.
It should be noted that the uncertainty estimated from the formal errors by least squares fitting might only represent the contribution from observation noise of the time series. It tends to underestimate the true uncertainty level since other impacts also contribute to the GMOM rate errors. One impact is from the uncertainty of the spatial filtering applied to the GRACE solutions. A recent study in Chen et al. [9] estimates this impact and indicates that different spatial filtering choices on GRACE data may lead to a difference up to 0.2 mm/year in the GMOM rate estimation (see their comparisons among the three cases of applying decorrelation and 300 km Gaussian smoothing, decorrelation only and no filters, in Table 5). Here we took ± 0.20 mm/year as the resulting uncertainty for the GMOM rate from the variety of different filtering choices.
Another impact results from the uncertainty of GIA models used to correct the GIA effect in GRACE observations. Currently there are several GIA models available for the correction in the GMOM rate estimation. Different GIA models may lead to notable differences in the correction. Chen et al. (2018, Table 6) [9] compares the GIA contributions to the GMOM rate from three commonly used models, including the Paulson07 [40], A13 (as used in this study), and ICE6G [41] GIA models. The comparison indicates the GIA contribution difference reaches up to 0.2 mm/year. Thus we took ±0.20 mm/year as the uncertainty corresponding to GIA model errors for the GMOM rate estimation in our study.
Considering the error sources from the observation noise (±0.11 mm/year), the spatial filtering uncertainty (±0.20 mm/year) and the GIA model uncertainty (±0.20 mm/year), we took the combined error estimate ±0.30 mm/year (from the square root of the sum of squares for the three) as a more reasonable uncertainty level for the GMOM rate in our study. As a result, the CSR RL06 GMOM rates should be 2.12 ± 0.30 mm/year and 2.05 ± 0.30 mm/year before and after seismic effect correction, respectively. Here we note that the GRACE's low-degree term errors (e.g., the degree-1 and C 20 terms) might also have notable contribution to the GMOM rate uncertainty. Therefore, the true uncertainty level could be even larger.

Uncertainty Estimate of Seismic Impact on the GMOM Change Rate
For the seismic effect correction, theoretically we could use predictions from earthquake dislocation models to remove the co-and post-seismic terms in the GRACE observations. However, the uncertainty of current fault slip models and Earth models necessary for the model predictions are fairly large, which leads to large errors for the seismic modeling. Different model inputs may lead to distinctly different co-and post-seismic predictions even for the same earthquake. Therefore, in this study we estimated the influence of seismic terms on the GMOM change by least squares fitting of GRACE observations, rather than using the dislocation model predictions.
Nevertheless, the time series fitting method has its own uncertainty on the estimation of seismic impact on the GMOM change rate. The choice of the post-seismic relaxation time parameter τ affects the fitting results, which inevitably leads to errors for the evaluation of the GMOM change rate. The value of τ varies for earthquakes in different regions. For example, de Linage et al. [35] shows that with a non-linear optimization method, the range of the estimated post-seismic parameter τ for the December 2004 Mw9.1 earthquake is 0.5-0.75 year (i.e., 6-9 months) in the region surrounding the epicenter. Li et al. [42] applies a residual minimization standard of the observations and indicates that the value of post-seismic parameter τ for the March 2011 Mw9.0 Tohoku-Oki earthquake is around 0.32 year (i.e., 4 months). Therefore, in this study we used τ = 7 and τ = 4 months to fit the post-seismic terms for the Sumatra and Tohoku-Oki earthquake regions, respectively. Since the post-seismic signal was relatively weaker for the February 2010 Mw8.8 Chile earthquake (compared to the co-seismic jumps and the fluctuations related with other effects and observation noise, see Figure 5c,d), the value of parameter τ should only have a fairly slight impact on the fitting result of post-seismic terms. Here, we empirically chose τ = 7 months to correct the seismic effect for the Chile earthquake region.
We approximately estimated the impact of fitting parameter τ on our GMOM rate by comparing the results with different τ values. Considering that the contribution from the Sumatra earthquake region was dominant among the seismic effect corrections (see Table 1), we took the Sumatra case as an example. We did the time series fitting with the values τ = 6, 7, 8, and 9 months and calculated the GMOM rates for the four cases (following the processing steps in Section 3.2). We found the maximum change from τ = 6 to 9 months was 0.0002 mm/year for the contribution to the GMOM rate. Corresponding results are listed in Table S1 of the Supplementary Material. Thus it should be reasonable to infer that the τ value choice influence on the GMOM rate in this study was at a level of 0.001 mm/year, and should be far less than the range of ± 0.01 mm/year.
The range of the spatial window selected to remove the seismic effects also contributes to the uncertainty of estimating the seismic impact on the GMOM rate. Currently we directly selected a rectangular area surrounding the epicenter as the correction region. However, in fact, an ideal choice would be choosing the correction area based on the spatial characteristics of seismic signals in the near field of the earthquake. The simple rectangular treatment leads to a visible discontinuity feature in the spatial distribution after the seismic correction in the Sumatra, Tohoku-Oki, and Chile earthquake regions (see Figure 3b, Figure 4b, and Figure 5b in Section 3.1). Here we varied the spatial window ranges in the estimation to investigate its influence on the results. Taking also the Sumatra earthquake region as an example, we compared the corrected GMOM rates in three cases with the ranges of 20 • × 20 • , 30 • × 30 • and 40 • × 40 • (see Table S2 in the Supplementary Material). The differences between each other among the three cases were within 0.01 mm/year.
The GRACE monthly solutions from different releasing centers may also lead to discrepancy of the results. According to Table 2, the discrepancies reached up to 0.02 mm/year among the estimates from CSR RL6, GFZ RL06, and JPL RL06 solutions (with −0.07, −0.08, and −0.09 mm/year for CSR, GFZ, and JPL, respectively). The 0.02 mm/year uncertainty related with different data processing centers was well above those contributions from the error sources discussed above. Therefore, we took ± 0.02 mm/year as a relatively reasonable uncertainty estimate for the seismic impact on the GMOM rate 2003-2015 in our study, which came out as −0.07 ± 0.02 mm/year.

Comparison with Estimates from Altimery and Argo
In addition to the GRACE estimation, the GMOM change rate could also be obtained from the Altimetry and Argo measurements (by removing the Argo steric contribution from the Altimetry sea level change observation). Due to the spatial coverage limitation of the Argo floats during its early observation period, the sea level change under a global scale is usually investigated with Argo data after 2005. Recently, Chen et al. [14] obtained a GMOM rate of 2.68 ± 0.15 mm/year for the period 2005.01-2015.12, with 3.79 mm/year (from Altimetry) minus 1.11 mm/year (from Argo). In an earlier study Dieng et al. [10] estimated the GMOM rate (2004.01-2015.12) as 2.35 ± 0.14 mm/year, with 3.49 mm/year (Altimetry) minus 1.14 mm/year (Argo). By comparing the above results one may find that the GMOM rate differs notably for different time spans, which is likely because the data period currently available is not long enough to precisely quantify the long-term trend rate. Moreover, the long-term trend quantification is also contaminated by the strong interannual variations that are difficult to efficiently remove in the current observations of limited data length.
For the period 2003.01-2015.12 considered in this study, the GMOM rate by Altimetry minus Argo is difficult to obtain reliably due to the Argo data availability problem over the earlier years of the period. If one assumes that the GMOM rate 2003.01-2015.12 by Altimetry minus Argo is close to that for the period 2004.01-2015.12, neither of our GRACE estimated GMOM rates before and after seismic correction (2.12 ± 0.30 mm/year and 2.05 ± 0.30 mm/year) agrees with the trend rate from Altimetry and Argo (2.35 ± 0.14 mm/year). The GRACE estimate of GMOM rate is noticeably smaller than that obtained from Altimetry and Argo, which is the same situation as in the previous studies Chen et al. [14] and Dieng et al. [10] (mentioned above). The discrepancy might attribute to the low-degree term errors as well as data filtering uncertainties associated with GRACE observations, model uncertainties in the GIA effect correction, Argo measurement uncertainties (especially for the deep ocean), etc. However, when considering the large uncertainty range of the GMOM rate, our GRACE result after seismic correction (2.05 ± 0.30 mm/year) can overlap with the Altimetry-Argo result (2.35 ± 0.14 mm/year) within the error range.
In an earlier study in 2019, Uebbing et al. [19] 12), their estimated GMOM rate is distinctly smaller than ours. Considering that we used the same A13 model for the GIA effect correction as used in Uebbing et al. [19], the discrepancy might result from the reasons as follows. Besides the data period difference that may contribute to the discrepancy, the GRACE data version difference may contribute notably to the discrepancy. In the calculation we used the RL06 solutions, while Uebbing et al. [19] use the RL05 solutions. Moreover, the applied leakage correction methods are different (i.e., the FM method in our study vs. the buffer zone method in Uebbing et al. [19]). Note that the FM correction method usually leads to a relatively larger GMOM rate (up to 0.2-0.3 mm/yr) than the buffer zone method, as indicated by Chen et al. [9,14].

Seismic Impact for Shorter Time Spans and Regional Scales
For the GMOM change rate estimation, the seismic correction may increase to a large extent when choosing a shorter time period that covers the occurrence time of the earthquake. As an example, we selected the period from January 2003 to December 2007, during which the contribution of the 2004 Sumatra earthquake was apparently dominant. The difference between the 2003-2007 GMOM rates before and after seismic correction reached up to −0.25 mm/year (see Figure S1 in the Supplementary Material), which indicates the seismic effect was significantly larger than the −0.07 mm/year correction for the 2003-2015 period.
The spatial scale of the study area also affects the results. For a global scale, the seismic effect was −0.07 mm/year (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), while for regional scales, the effect might be significantly different. For instance taking the Indian Ocean as the study area, the mass change signals in this region caused by the 2004 and 2012 Sumatra earthquakes will affect the regional ocean the mass rate more dramatically. As a result, the corresponding seismic effect on the Indian Ocean mean mass change rate could be much larger than that on the GMOM rate.

Conclusions
We studied the impact of the seismic effect caused by the December 2004, April 2012 Sumatra earthquakes, the March 2011 Tohoku-Oki earthquake, and the February 2010 Chile earthquake on the GMOM rate estimated by GRACE. We used least squares to fit the co-and post-seismic terms, and applied the FM recovery method to correct the land-ocean leakage effect. Results show that GRACE GMOM rates before and after seismic correction were 2.12 ± 0.30 mm/year and 2.05 ± 0.30 mm/year. The uncertainty level ± 0.30 mm/year was estimated by considering the contributions from the GRACE observation noise, the uncertainty due to various spatial filtering choices and the GIA model errors. The bias caused by the four large earthquakes on GMOM rate was −0.07 ± 0.02 mm/year, with the uncertainty estimated by investigating various error sources. Although the seismic effect on the GMOM rate turned out under the uncertainty level of the current GRACE GMOM rate estimation at a global scale (−0.07 vs. ± 0.30 mm/year), it is a systematic bias. It should be considered for improved quantification and understanding of the global sea level change. For regional scales, this bias could be more significant. Moreover, the seismic effect would affect the GMOM rate more dramatically when focusing on shorter time spans. Consequently, it is necessary to consider the seismic effect correction for more precise estimation of GMOM change by satellite gravimetry observations from GRACE, GRACE Follow-On, and especially for the future satellite gravimetry missions with expected higher resolutions and precisions.