Correction of GRACE measurements of the Earth’s moment of inertia (MOI)

The widely used 15-year Gravity Recovery and Climate Experiment (GRACE) data sets do not conserve global total mass. They have a spurious decreasing trend of ~ 280 Gt/year. Various regions contribute differently to the global total mass loss error, with the Greenland Ice Sheet (GrIS) generating ~ 10% of the error alone. Atmospheric parameters from reanalysis datasets drive a well-tested ice model to generate mass variation time series over the GrIS for 2002–2015. Because shorter timescale spikes of ~ 10–30 Gt in GRACE measurements are physically based, only the overall trend of ~ 30 Gt/year requires correcting. A more accurate mass loss rate estimate for 2002–2015 is ~ 120 Gt/year, considerably below previous estimates. With the water redistribution to lower latitudes and other effects from a warming climate, the nontidal Earth moment of inertia (MOI) also increases. After rectification, the GRACE measured mass redistribution shows a steady, statistically robust (passed a two-tailed t-test at p = 0.04 for dof = 15) rate of MOI increase reaching ~ 10.1 × 1027 kg m2/year, equivalent to a 10.91 μs/year increase in the length of a day, during 2002–2017.


Introduction
In recent decades, the net top-of-atmosphere energy imbalance reached 1 W/m 2 (Hu and Bates 2018;Trenberth and Fasullo 2013). With the increased energy accumulation in the Earth's climate system, the thermodynamic structure of the atmosphere has been substantially altered. Consequently, the hydrological cycle, the oceanic circulation and the cryosphere dynamics all have been changed. According to the IPCC Special Report on Ocean and Cryosphere in a Changing Climate (2019), the upper (0-200 m depth) ocean heat content had increased at a rate of 8-18 × 10 21 J/ year. Polar ice sheets respond to this climate change by increased creeping and other related total mass shedding mechanisms (Ren and Leslie 2011). As mountain glaciers retreat (Ren and Karoly 2006), significant portions of solid water become liquid water and are redistributed into global oceans or drained into reservoirs, such as soil moisture and groundwater. The kinetic energy of winds is very unevenly distributed spatially and is concentrated in polar and subtropical jet streams. Subtropical jet stream changes have an important influence on the angular momentum of the Earth system (Lee et al. 2019). Polar jet streams influence mid-latitude weather systems, with the storm tracks being essentially a surface expression of the jet stream. Of great importance is that the mid-latitude meridional temperature gradients are being modified by anthropogenic climate change (Vallis et al. 2015; and the jet streams are expected to further adjust in response to these changes Vavrus 2012, 2015;Haarsma et al. 2013). While kinetic energy eventually dissipates in the boundary layer, lasting effects on Earth result from the associated mass redistribution (e.g., precipitation), which will affect the Earth's moment of inertia (MOI).

3
The high latitude Greenland Ice Sheet (GrIS) is a leading contributor to global sea level rise (Box et al. 2006;Shepherd and Wingham 2007;Rahmstorf 2010;Shepherd et al. 2012Shepherd et al. , 2018Hanna et al. 2013;Vaughan et al. 2013;Bamber et al. 2018), accounting for ~ 0.7 mm/year. In a future warming climate, faster ice creeping from reduced ice viscosity, enhances dynamic mass loss (Ren et al. 2011a; Ren and Leslie 2011;Rignot and Kanagaratnam 2006). However, estimates of future changes in net mass balance are highly debatable as studies suggest GrIS snowfall will increase as the climate warms (Ettema et al. 2009;Merz et al. 2013;Fyke et al. 2014;Hanna et al. 2018). There also are studies of precipitation and circulation relationships over the GrIS, under present climate conditions (Box et al. 2006;Cassano et al. 2001;Fettweis et al. 2013) and in a warming climate. Doyle et al. (2015) found that summer storms accelerated ice flow, with consequences for the GrIS in a future warmer climate.
The Gravity Recovery and Climate Experiment (GRACE), launched in 2002, provides large scale, coherent changes in the Earth's mass field. It monitored temporal variations of the Earth's gravity field with monthly temporal resolution (Bentley and Wahr 1998;Famiglietti and Rodell 2013;Ren 2014a;Voss et al. 2013;IMBIE Team 2018). Horizontal mass redistribution in the atmosphere, the oceans, the hydrosphere, the cryosphere or the lithosphere, all are captured by GRACE measurements. Here, we first correct a readily identifiable error in GRACE measurements, for 2002-2017, which provides rectified datasets suitable for studies in a wide range of disciplines. Then, as one application, the total Earth MOI variations and the various contributors are analyzed using the Scalable Extensible Geofluid Modeling framework for ENvironmenTal issues (SEG-MENT, Ren 2014b; Ren et al. 2008Ren et al. , 2012 modeling system. The SEGMENT is the basic modeling system in this study. It has successfully been applied to both the Antarctic Ice Sheet (AIS, Ren et al. 2012) and the GrIS (Ren et al. 2011a, b) for ice sheet total mass balance, to world mountain regions for landslides and debris flow studies (e.g., Ren 2014b), and to various faults for earthquake mechanism studies (Ren and Fu 2019). The primary aim of this study is to quantify the total Earth MOI contributions from its various components. Identifying the persistent contributors during this 15-year period assists in the estimation of the climate-driven MOI evolution, in the transient climate change period.

Method and data
The focus of this study is the MOI fluctuations caused by near surface, climate sensitive mass transport phenomena, such as those from the altered hydrological cycle, atmospheric expansion, and cryospheric mass balances. Multiple data sources and an advanced modeling system are used to provide estimates of the changes in the Earth's total MOI, and in the contributions from the individual MOI components.

GRACE measurements
Because mass changes integrate spatially, GRACE rapidly measures net mass balance over vast regions, including polar ice sheets, thereby providing mass change rates by repeated measurements over the observation period. Measuring the gravity anomaly at monthly temporal frequency and ~ 3° spatial resolution (Vishwakarma et al. 2018), GRACE is ideal for measuring the evolution of the Earth's MOI. Here, data are obtained from the Jet Propulsion Laboratory (JPL) global mascon solution (Watkins et al. 2015), at grace.jpl. nasa.gov, after atmospheric vapor corrections and application of other time-invariant harmonic filters.
Monthly averaged global gravity field data from GRACE then are estimated in the form of spherical harmonics (i.e., Stokes coefficients). For many applications involving mass fluxes, the traditional procedure of direct using of spherical harmonics has its limitations (Luthcke et al. 2013;Save et al. 2016;Tapley et al. 2019). Mass concentration blocks (mascons) basis functions typically are used to better estimate mass fluctuations from GRACE measurements. The JPL mascon used in this study (Watkins et al. 2015) applied explicit partial derivatives with analytical expression for mass concentration to relate the intersatellite range-rate measurements to the individual mascons. Specifically, it is solved on 3 × 3 degrees equal area elements, although the final products are interpolated to 0.5 by 0.5 degrees grids. Water equivalent mass flux data are used here. For this purpose, a separate scale factor dataset is used, derived from the community land surface model predicted terrestrial water storage changes. The scale factors, which also are on 0.5 by 0.5 degrees grids, are multiplied to the mascons to eventually obtain the mass grids in units of equivalent water depth. The post glacial rebound (PGR) effects, also called glacial isostatic adjustment (GIA), are removed from the mascon solution using model predictions. Except for the RL05 release of CSR mascons, an ellipsoidal Earth model usually is adopted for mascon solutions (Table 1).
Mascon solutions are used in this study because they have fewer leakage errors, which is usually achieved by the separation of land and ocean signals. Further processing involves representation mascon solutions into finer grids. Spheric harmonic treatments and filtering operations are spatial operations Swenson and Wahr 2006;Duan et al. 2009;Longuevergne et al. 2010;Frappart et al. 2011;Wouters et al. 2014), and do not affect global total mass conservation. The monthly datasets are the 'original raw data sets' which provide mass grids, in units of equivalent water depth; hence, the mass unit gigatonnes (Gt, or 10 9 tonnes) is equivalent to the volume unit km 3 .
Although the Earth's total mass should be constant, it is not well conserved in the GRACE measurements from the JPL website (Humphrey et al. 2016;Wiese et al. 2018), as shown in Fig. 1a. It also is not conserved in both the CSR and GSFC/NASA mascon solutions. Globally, it is accepted that the negative trend in raw GRACE data is spurious, with no physical meaning. For regional areas, however, physically meaningful information is contaminated by this artificial signal because, globally, the spurious error is clear whereas, at regional scales, it is not clearly separated from the physical signals. In addition to systematic measurement errors, the GRACE data are contaminated by noise. Here, the empirical mode decomposition (EMD) of Wu et al. (2007Wu et al. ( , 2016 determines the spurious trend, and detrends the raw GRACE data. EMD identifies signals of differing inertia (e.g., intrinsic trends and natural variabilities on various time scales). The advantages of EMD over traditional Fourier based decomposition are detailed in Wu et al. (2016). EMD also identifies error sources for cyclical signals and generic trends in raw GRACE data.
In summary, the approach of this study is to: (1) rectify the raw GRACE measurements by removing the artificial (nonphysical) decreasing trend; (2) use the rectified mass variations to train SEGMENT-Ice driven by reanalysis data; and (3) employ the optimized SEGMENT-Ice modeling system for MOI decomposition into its various components. Hence, it is critical to determine the global distribution of GRACE satellite errors: i.e., how much each grid point on Earth contributes to the global total error. From Fig. 2, mass variations are globally uneven. For a 5 mm water equivalence criterion, over 50% of the Earth's surface experiences no mass variations of that magnitude. Strong mass field fluctuations occur in several regions with clear precipitation and/or plate convergence, including the Amazon region, central and southern Africa, monsoonal India, the Bay of Bengal, and south-central Greenland. In these regions, mass fluctuation amplitudes reach 20 cm. It is a dynamic situation in that mass deficit centers can become mass gain centers over time. However, regionally, mass fluctuation magnitudes basically are stable. The absolute error contributing to global total mass fluctuation time series, rather than spread evenly from all global areas, may be proportional to the magnitudes of regional mass fluctuations. A weighting scheme partitions total spurious mass loss into regional areas: Here A i,j is the magnitude of raw GRACE mass variation time series over the grid (i,j), summed over the globe (i = 1, 360; j = 1,180 for 1° degree raw data). GrIS contributes ~ 10% of global total mass variation. This scheme filters GrIS measurements, producing a de-trended mass variation time series, further used for SEGMENT-Ice verification and future GrIS mass loss projections. This linear error ascription is supported by Humphrey et al. (2016), despite differing approaches to detrending and decomposing time series.

Earth's moment of inertia (I)
The MOI (I) is defined as: (1) where is latitude, is longitude, is density, and R is Earth's radius. As the earth has an angular velocity of = 7.29 × 10 −5 rad/s and a moment of inertia I ≈ 8.04 × 10 37 kg m 2 , it therefore has a rotational energy ( 0.5I 2 ) of ~ 2.138 × 10 29 J. The coordinates follow the right-hand rule and all variables are in SI units. The total MOI variations are decomposed into various contributors, for increased understanding long-term trends and short-term variations. Hence, contributions from the polar ice sheets (GrIS and AIS), mountain glaciers (MG), atmospheric expansion (Atm), and precipitation driven hydrological cycles (groundwater and soil moisture, SM) all are included. Estimation formulae for the respective contributors are indicated, below.
a. MOI estimates from GRACE gravity anomaly measurements ( I total ): where subscripts i, j are integer indices of the latitudinal and longitudinal count of the horizontal global grids, W is grid area, Δh is the water equivalent depth (GRACE measurements are in cm, converted into meters), and w is water density. The Earth's MOI is in kg m 2 . b. MOI estimates from the elevated atmospheric mass center of the Earth's atmosphere ( I atm ): (2) Panel a shows the raw GRACE measurements (only area weighting, or a static spatial operation, is applied to the 1-degree gridded data) of mass variations for the northern (a global belt of 0°-60° N) and southern (0°-60° S) hemispheres, and for the entire globe. Data are from the JPL public website https:// grace. jpl. nasa. gov. Panels b-f are decomposed signals, ranging from white noise to slowly varying signals. Global mass conservation is a basic physi-cal requirement. However, this quantity is not conserved in the raw GRACE data, and the associated unrealistic trend should be removed for many applications. Whereas the total global mass loss is spurious (e.g., from errors due, for example, to the ageing of the GRACE satellites), the mass loss trend in both the Northern and Southern Hemisphere, and in regional areas (not shown), all are physical phenomena contaminated by errors where subscribe k is integer count of atmospheric vertical levels, p is atmospheric pressure, z = z(p) is the geopotential height, and g is gravitational acceleration. c. MOI from polar ice sheets and mountain glaciers (I ice ): Essentially the same expression as Eq. (2) is followed except that h is now simulated by an advanced ice sheet modeling system SEGMENT-Ice (Ren et al. 2011b). Some technical details are described in the Appendix. d. MOI from hydrological cycle as a remainder ( I pr ): Essentially the same expression as Eq.
(2) is followed except that h is now simulated by an advanced land surface modeling scheme in SEGMENT-Landslide (Ren et al. 2011c).
Atmospheric variables which are used as inputs to SEGMENT are obtained from a reanalysis dataset. Anomalies (ΔIʹs) of all components are calculated by subtracting from each annual value the value of year 2002. (3)

Reanalysis data and SEGMENT modeling system
The calculation of the MOI components is obtained either directly from atmospheric parameters (e.g., the expansion of Earth atmospheric columns and thus the elevation of their centers of mass) or is obtained using atmospheric parameters as forcing of the SEGMENT modeling system. For example, driven by atmospheric parameters, the SEGMENT system is used to simulate groundwater fluctuations and ice mass loss from the cryosphere. Atmospheric parameters in this study are from reanalysis datasets.  Suarez et al. (2008), and the JRA-55 reanalysis (Kobayashi et al. 2015), at respectively horizontal resolutions of 0.75°, 0.25°, 2.5°, 0.625° × 0.5° and 1.25°. The use of five independently produced reanalysis datasets allows the quantification of the sensitivity of our results to uncertainties in the state of the atmospheric forcing. Six-hourly data from 1980 to 2017 are used, which fully covers the period of GRACE measurements (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). The reason for restricting the temporal coverage to the remote sensing era is due to the large Most areas have no significant mass changes (color shades, in equivalent cm water depth; using 0.5 cm as the threshold). Strong mass field fluctuations are concentrated in several regions with precipitation and/or apparent plate convergences, including the Amazon region, central and southern Africa, the monsoonal region of India and the Bay of Bengal, and south-central Greenland. It is dynamic, as mass deficit centers can become centers of mass gain. However, for a fixed region, magnitude mass fluctuations are relatively stable. The errors in the GRACE data might not be affected by geography. The absolute error introduced in the global total mass fluctuation time series, rather than being contributed evenly, may be proportional to regional mass fluctuation magnitudes. The fluctuation magnitudes, weighted by the representing area, distribute total errors over global regions. This approach filters GrIS mass variation measurements to obtain detrended time series which are further used for ice model verification and for projections of GrIS mass loss uncertainty in the Southern Hemisphere, due to the lack of observational data coverage (Le Marshall et al. 2013). Time series of globally integrated surface pressure are examined to guarantee the total mass conservation is met by each reanalysis dataset. Other details of the data used in this study are in Table 1. The total GrIS mass loss results from complex, interlinked processes, that require a physical modelling system for quantitative estimates of mass loss trends. The SEGMENT-Ice modeling system is a numerical treatment of the physical processes of ice dynamics and has been used in other cryosphere studies (Ren et al. 2011a(Ren et al. , 2013aRen and Leslie 2011). The conservation equations of mass, momentum and energy simulate the ice motion and mass variations. Model settings and spin-up are those in Ren et al. (2011b) and Ren and Leslie (2011). Model tuning uses corrected GRACE measurements for the 2002-2015. Reanalysis data sets provide the atmospheric parameters driving SEGMENT-Ice for the GRACE measurement period. SEGMENT-Ice was used for both GrIS (Ren et al. 2011b) and the AIS (Ren et al. 2013b), and it allows multiple compositions, multiple phases, and multiple rheologies. The SEGMENT-Ice design also includes known physical mechanisms of ocean-land iceatmospheric interactions, and new mechanisms can readily be implemented as they become available. For example, it has a granular bed-sliding enhancement mechanism (Ren et al. 2013a;Ren and Leslie 2011), calving mechanisms (Ren and Leslie 2014), and unique fatigue mechanisms explaining random, unexpected terminal glacier surges.

Results
The Earth's global mass should be time-invariant. However, Fig. 1a shows a decreasing global total mass trend of ~ 280 Gt/year. To diagnose this trend, the EMD technique (Wu et al. 2016) decomposes the total, noisy, signal into six components. Except for the first component, white noise, in Fig. 1 (not shown), all other five components have physical meanings. Components 2, 3 and 4 respectively indicate seasonal, inter-annual, and multiple year cyclical signals. Component 5 shows an overall trend of ~ 200 Gt/year. The cause of this 'apparent' mass loss is unclear but has no physical meaning. Possibly it is due to satellite aging or arises from other technical issues (Watkins et al. 2015). Using the weight function described above, ~ 30 Gt/year of the global mass loss can be credited to the GrIS. After this trend correction, the mass loss rate during 2003-2015 is ~ 120 Gt/year from the GrIS, which still is larger than estimates from InSAR and other measurements, but below that of Velicogna et al. (2005). From Fig. 1a, the global total mass time series also has annular cycles. This cyclical part of the spurious error is largely attributable to tropical regions between 30° S and 30° N, primarily the Amazon region. The in-phase variation of global and Southern Hemispheric belt time series in Fig. 1a is not coincidental. In comparison, the GrIS is insignificant (two orders of magnitude smaller) and is not removed from the GrIS mass loss time series. For other regions of interest, removal of this portion of the spurious signal requires the information in Fig. 1b, c. Removing the overall trend is critical for estimating the GrIS contribution to eustatic sea levels, as shown in Fig. 1f. The two hemispheric belts, benefitting from the net mass loss from the AIS and the GrIS, should experience a net mass gain, as indicated by the most slowly varying components. Unfortunately, due to spurious noise, mass gain trends are insignificant during the observation period (Fig. 1f), for both hemispheric belts. After trend removal, the global signal time series of component 5 are almost zero, whereas for both the Northern and Southern Hemisphere it increased steadily. The Northern Hemisphere (0°-60° N) increase is larger, suggesting that the GrIS net mass loss exceeds the AIS loss (Fig. 1f).
The total mass balance includes dynamic mass balance, due to ice flow convergence/divergence, and surface mass balance, from precipitation and melting. The surface mass balance has strong seasonal variability. The GrIS melting season is May-October. Accumulation occurs year-round but net accumulation is mainly in the frozen season, typically from several snowstorms. Thus, there is steep annual mass accumulation in November-March, followed by a gradual decrease (Fig. 3). The melting season spikes in 2005 and 2007 were snow precipitation from hurricane remnants (TCs Wilma and Noel) reaching high latitudes. Other spikes, primarily in the frozen accumulation season, are largely from polar snowstorms arriving from the south-east, preceded by strong, easterly Fohn winds. The total mass has a significant decreasing trend of ~ 150 Gt/year before correction and ~ 120 Gt/year after rectification (solid and dotted lines in Fig. 3). This is below the estimates of ~ 200 Gt/year of Velicogna and Wahr (2006), and below the ~ 180 Gt/year of Chen et al. (2006). Notably, the seasonal mass variability is small.
SEGMENT-Ice, driven by WRF-ARW model projections, provided encouraging simulations of the GrIS mass loss (Fig. 3). Summer mass loss spikes are captured well, as in year 2005, and are due to the south-east sector of the GrIS being affected by Atlantic TCs approaching Greenland (e.g., TC Wilma 2005; Fettweis et al. 2013). In addition to being a mechanism for the vertical and poleward transport of atmospheric momentum and energy, TCs transport moisture from the tropical oceans to extra-tropical regions (Emanuel and Jagger 2010;Knutson et al. 2010Knutson et al. , 2013. Water vapor transported poleward by TC remnants is precipitated as snow over the ice sheets and glaciers and is critical in maintaining the global eustatic sea level. For example, TCs Wilma (2005) and Noel (2007) both match observational spikes. Because snow precipitation from TC remnants is 'warm' precipitation (i.e., snowfall warmer than surface ice and, for the ablation zone, some liquid precipitation), the accumulated mass is not persistent, enhancing the dynamic ice loss. The net results are spike-followed-by-trough signals in the decomposed time series (order 2 time series in Fig. 1) of monthly GRACE measurements. Comparing simulated and observed total mass loss, the spikes of 10-30 Gt from GRACE measurements are realistic, thereby requiring an adjustment only in the unrealistic trend possibly from GRACE satellite aging. SEGMENT-Ice parameters were optimized using GRACE measurements. Upon convergence, the root-mean-squared error over the 13-year period is below ~ 20 Gt. Hence, the model simulations and observations are acceptably close.
Different mascon solutions have different definitions (spatial resolution, converging schemes, and even the shape of mascon boxes). For example, the CSR mascon solutions are related to the range-rate observations via partial derivatives with respect to a spherical harmonic expansion that is truncated to degree and order 120. Mass anomalies in each mascon are computed from satellite range rate observations through their partial derivatives (e.g., Save et al. 2016). The CSR resolutions also are estimated directly on a geodetic grid of the size equivalent to equatorial 1°, then to the finer resolution of 0.5° as used here. In contrast, the JPL mascon solutions were solved on 3° elements and represented on 0.5° grids. The global total mass balance is well preserved by both mascon solutions. However, the spatial distribution patterns of the mass flux are different, as are the MOIs. The global total MOI trend estimated using CSR mascon solutions is shown in Fig. 4. Importantly, both solutions suggest a steady, statistically robust rate of MOI increase, passing a two-tailed t-test with p = 0.04 for dof = 15, that reaches ~ 10.1 × 10 27 kg m 2 /year, which is equivalent to a 10.91 μs/ year increase in the length of a day, during 2002-2017. In particular, the short-term spikes associated with hydrological patterns linked with major earthquakes [e.g., the 2006 spike associated with the Nepal earthquakes, Ren and Fu (2019)] agree very well. Repeating the above analyses using the CSR mascon solutions suggest, qualitatively, the same conclusions as those from the JPL mascon solutions. Interestingly, the mascon solution from the GSFC, produced a MOI evolution curve closely resembling those from CSR and JPL mascon solutions. The total MOI trend identified by the three mascon solutions is similar, and hence is suggestive of a robust signal. The MOI annual means obtained from the GRACE measurements and the estimated contributing components are estimated (Fig. 5). The linear trends during 2002-2017 for the five contributors: the AIS, GrIS, Atm (atmosphere expansion), MG and SM are, respectively, 5.06 × 10 27 kg m 2 / year, 4.8 × 10 27 kg m 2 /year, 1.1 × 10 27 kg m 2 /year, 1.94 × 10 27 kg m 2 /year, and 0.4 × 10 27 kg m 2 /year. Only the SM has a p-value > 0.1 and fails the Student-t test. Other contributors all have p-values less than 0.1. Specifically, the AIS and GrIS have p-values of ~ 0.035. Hence, except for the SM, other components all contribute steadily to increasing the MOI. In absolute terms, the AIS is the largest contributor. Through the increased creeping of the ice sheet, the AIS loses a similar amount of ice to oceans to the GrIS (Ren et al. 2011b(Ren et al. , 2013aChen et al. 2006;Velicogna and Wahr 2006). However, because it is located at higher latitudes than the GrIS, its mass loss to oceans signifies more salient increase of Earth's moment of inertia. The GrIS is the second largest contributor to the increased MOI, of about an order of magnitude larger than mountain glaciers and that from elevated center of mass of Earth atmosphere. Polar amplification is more salient in the Northern Hemisphere than in the Southern Hemisphere, primarily because the relatively low temperature over Antarctica has not yet activated an albedo feedback (e.g., Charney et al. 1977;Charney and Shukla 1981). When this feedback mechanism is activated in the future due to, for example, increased ground surface exposure from ice cover retreats (or simply from seasonal  Fig. 4 so not plotted here) ice surface melting), the AIS possibly will play an increased role in modulating the Earth's MOI.
Although the contribution from atmospheric warming is the smallest of the four (Fig. 5), the upward trend is steady and started well before 2002 (Fig. 6). For the entire meteorological remote sensing era (1979-present), there is essentially no significant negative contribution. The MOI fluctuations associated with the expansion of the oceans, although present in the GRACE measurements, is not separated out in this study. The magnitude reported by Landerer et al. (2007) is a good reference as it is on the same order of magnitude (− 0.6 μs/year). Similarly, Gross (2000) concluded that the atmospheric fluctuations (surface pressure redistribution) are only about one-half of those caused by ocean bottom pressure fluctuations in affecting the MOI. Although their results revealed the equal importance of atmospheric fluctuations and sea level rise, the upward lift of the center of mass of the Earth's atmosphere (to higher altitudes and away from the axis of rotation; Fig. 6, reaching 1.6 m/year rate at present at polar regions) likely was overlooked. However, this vertical expansion appears to be significant, contributing 0.17 μs/year alone. Horizontally, the mass re-distribution tendency from the polar regions to the mid-latitudes, due to enhanced mid-latitude high level jets (Lee et al. 2019), is the main driver of the MOI increase. During the GRACE observational period, it has a rate of ~ 0.33 μs/year, close to the estimates of Gross (2000). If only the ocean expansion is considered, it opposes the atmospheric effects on the MOI, and they effectively cancel each other. As the climate warms, the regional disparity in mass accumulation may play an increasingly important role in affecting the Earth's MOI. For example, the intensified Antarctic Circumpolar Current, through Coriolis effects, forces ocean water away from the AIS to lower latitudes. The mass accumulation trend over the Southern Ocean hence contributes to an increased MOI. Notably, it opposes the ocean expansion reported by Landerer et al. (2007).
A unique feature of the contribution from the soil moisture and groundwater (Fig. 7) is that its annual variation is in phase with the total MOI as indicated by the GRACE measurements. Consequently, it is the most dynamic contributor to fluctuations in the MOI. For a regional area, enhanced mass can be a result of increased precipitation or enhanced soil moisture/groundwater retention ability. The root mechanism still is not clear, but a plausible explanation is that most of the major earthquakes during this period occurred at lower latitudes (Fig. 7). The energy released then enhanced precipitation distribution in lower latitudes. This might explain the spikes. Comparison of Figs. 4 and 7 indicates that each spike (from earthquakes) elevated the MOI to a new, higher value. For example, the moderate 2004 spike related to the Indian Ocean earthquake (Sumatra Earthquake 2004), would have caused the length of day to shorten by 3 μ s and restored slightly, as the response of upper mantle material. The 2008Wenchuan, 2011Tōhoku, 2012Aceh, and 2015 Nepal earthquakes all can be seen in the MOI time series. Without these spikes, the trend of the soil moisture contribution would be slightly decreased, agreeing with a polar shift of precipitation as the climate warms (Dore 2005;Solomon et al. 2007;Putnam and Broecker 2017). There is post-seismic evidence of enhanced water retention capability of the aquifers, through various hydrologic The shading indicates the spread using different reanalysis atmospheric parameters to drive the SEGMENT modeling system processes induced by earthquake-induced changes in hydrologic and mechanical properties of the crust (e.g., Ingebritsen and Manga 2019). Recently, Zhang et al. (2019) found that earthquakes can increase aquifer and aquitard permeabilities by fracturing/compaction, implying an enhanced groundwater retention. Accordingly, the mass disturbance is ~ 1.1-4.7 km 3 for the 2010 Chilean earthquake. However, the spike in the GRACE curve is > 10 km 3 , suggesting that there are other indirect effects on the fluid sphere from tectonic earthquakes. Over the 15-year GRACE period, the presence of negative MOI contributions that are large enough to cancel the increasing MOI trend from the other components is not supported by the GRACE observations.

Discussion and conclusions
Despite its design lifetime of 8-9 years, the GRACE mission has provided almost 15-years of quality data. Errors in raw GRACE total mass variations have been identified and corrected. The primary error source is satellite aging. Using the error-corrected data, the total mass loss rate from the GrIS during the 2002-2015 is ~ 120 Gt/year, which is ~ 20% less than in previous studies but is much closer to other geodetic measurements.
The rectified GrIS total mass loss time series were wellsimulated by the SEGMENT-Ice model, including small but important spikes produced by snowfalls from tropical cyclone remnants. These seemingly noisy spikes of ~ 10-30 Gt in the GRACE measurements thus are realistic. Hence, only correcting the overall trend is required. Model parameters are optimized with corrected GRACE data and are used to determine the total non-tidal MOI fluctuations for 2002-2017. The modeling system quantifies the contributions from various components. In this study, we found that the mass shed from the polar ice sheets, mountain glaciers retreats, and the atmospheric fluctuations (i.e., the horizontal re-distribution and uplift of the mass center of the atmosphere) all contribute to increasing the Earth's MOI and hence to a small, but identifiable, slowing of the Earth's rotation rate. Components of the GRACE time series were analyzed and the abrupt changes in the MOI correspond well with major earthquakes. Instead of the relocation of crust and upper mantle material, it is the associated groundwater, or precipitation redistribution, that dominate the MOI fluctuations. This may be another facet of the relation between droughts and major earthquakes (Ren and Fu 2019). Notably, the dynamic nature of permeability, as revealed by co-seismic and post-seismic hydrologic phenomena (Ingebritsen et al. 2006), has implications for groundwater systems in the fault vicinity. With advances in earthquake hydrogeology, quantitative detection of the earthquake contribution to the MOI changes is expected. Hitherto, sea level rise and mass redistribution effects on the MOI have been actively studied (Adhikari and Ivins 2016;Munk 2002;Seoane et al. 2012;Mitrovica et al. 2006). However, other components have not yet been studied extensively (e.g., the land hydrological cycle), or remain overlooked (e.g., the expansion of Earth's atmosphere). Although we do not claim a complete decomposition of the total MOI trend, other mechanisms represented in the residual term likely have secondary roles in the Earth's MOI evolution during the GRACE observational period.
It was found that there is a statistically significant (p = 0.039) total increase in the Earth's MOI, during the GRACE 15-year period. This steady increase in MOI has five major contributors. The two leading contributions are from the AIS and GrIS, with trends in both passing twotailed t-tests at p = 0.035 (dof = 15). Considered jointly, the linear trend of the polar ice sheets reaches 8.7 × 10 27 kg m 2 / year (equivalent to a 9.4 μs/year increase in the length-ofa-day). In the foreseeable future, this is the dominant term for MOI changes, as the climate warms. Among the climate warming contributors only the changes in rainfall redistribution has a statistically non-significant trend (p = 0.604). The mass redistribution caused by the hydrological cycle also is the most dynamic, contributing to fluctuation spikes in the global MOI. The ice sheet contribution is steady and is a persistent trend in the recent MOI evolution during the GRACE observational period. Compensating mechanisms, such as the GIA, operate on a much longer time scale. The increasing MOI trend is likely to continue in the transient climate warming period (the upcoming several hundred years), before the GIA and other possible negative factors dominate. Except for the wavier hydrological cycles, other contributing components are responses to a warming climate in synergy, rather than opposing each other. The Earth's MOI variations during the 15-year GRACE period of reliable observational data form a basis for identifying climate change impacts on the Earth's MOI and provide more precise insights into fluid and solid mass transport near the Earth's surface.
(snow precipitation adjusts the surface and melt/sublimation at the ice-bedrock interface) and the dynamic mass balance ice resulting from flow divergence/convergence. Generally, ice creeps from the central region to the peripheral regions, where it either melts or drains to the ocean through the inlets. The oldest ice apparently is located at the bottom of the central region (for GrIS, ~ 50 km north of the Summit) and is of age ~ 110 kyr. The model simulation domain contains the entire ice domain plus 1 km thick bedrock with assumed viscosity of 10 19 Pa s. Thus, the deformation to the bedrock elevations is estimated from the ice weight. The current geometry of the GrIS primarily is a result of snow precipitation climatology (direct precipitation and the turbulent wind-caused re-distribution, especially at the early stage of snow accumulation over the GrIS). Ice flow fields are an integral component of the ice model. With time marching, the ice velocities are updated together with the mass and thermal fields. There is a long-term glacial residual, which is small but negative, in total mass balance. For simulating mass changes of the GrIS over the next several centuries, this glacial residual is assumed to be constant. Spin-up simulations, using paleo-climate records of precipitation and analytical surface air temperature, arguably provide only the general GrIS mass, flow and thermal configuration.
Remote sensing observations, especially for the most recent decade, make available surface ice flow measurements (e.g. NASA's Making Earth System Data Records for Use in Research Environments, MEaSUREs, http:// nsidc. org/ data/ measu res/; Ren et al. 2013a, b) and gravity field fluctuations (e.g. Gravity Recovery and Climate Experiment, GRACE, http:// www. csr. utexas. edu/ grace/ scien ce/ gravi ty_ measu rement. html; Ren et al. 2011a). This information can provide further constraints on 3D ice temperature fields. We used an adjoint-based 4DVar scheme to fully extract the information. Unfortunately, the satellite era also coincides with a salient climate warming period. The topmost tens of meters of ice already has felt the warming and crept faster than under natural variability. To account for the effects from recent warming, it is assumed that the mean net snow precipitation over 1960-1990 perfectly balances the dynamic mass balance term everywhere over the GrIS (after adjusting for glacial residuals). Technical details of SEGMENT-Ice can be found in Ren et al. (2013a). Application to the GrIS are detailed in Ren et al. (2011a, b, c).