1 Introduction

Polar motion, which is the motion of the Earth's rotational axis relative to its crust, has been routinely observed using space geodetic observation techniques for more than a century (Anderle, 1973; Anderle & Beuglass, 1970). It is described in an Earth-fixed reference frame, where the origin is the Conventional International Origin and the X and Y coordinate axes point toward the 0° and 90°W longitudes, and it is given in units of milliseconds of arc (mas; 1 mas ≈3 cm). Compared to its motion range (tens to hundreds of mas), its observed uncertainty, which is much less than 1 mas, is negligible. Several recent studies (Adhikari & Ivins, 2016; J. L. Chen et al., 2013; Roy & Peltier, 2011) state that this long-term, high-accuracy record of the Earth's pole provides key information about climate change.

Generally, polar motion is triggered by the mass redistribution and movement of the solid Earth, atmosphere, and hydrosphere on various time scales (Höpfner, 2004). Two dominant periodic oscillations of polar motion, that is, the annual and Chandler wobbles are mostly excited by winds, ocean currents, atmospheric pressure fluctuations, and ocean bottom changes (J. L. Chen et al., 2012; R. S. Gross, 2000; Höpfner, 2004). The low-frequency components of polar motion, especially the linear polar drift over the 20th century, are considered to be mostly related to the solid Earth, especially the glacial isostatic adjustment (GIA) process and mantle convection (Adhikari et al., 2018; R. S. Gross & Vondrák, 1999; McCarthy & Luzum, 1996; Nakada et al., 2015). Beyond these predictable variations, other irregular components of polar motion, especially interannual and interdecadal movements, are likely climate-related and have attracted extensive research interest. R. Gross and Poutanen (2009) first noticed a directional change in secular polar drift in the mid-1990s. Roy and Peltier (2011) further showed that the direction of secular polar drift significantly changed from 68°W (1976–1992) to 58°W (1992–2008). Surprisingly, this change is not caused by either GIA or the mantle convection. Roy and Peltier (2011) suggested that it might be influenced by accelerated ice melting over glacial areas under global warming. However, a lack of key observations, that is, terrestrial water storage (TWS), brings challenges to further studies. TWS represents all forms of water stored on or below the land surface, including snow, ice, groundwater, soil moisture, surface water, and water stored in the canopy and biomass (Hirschi et al., 2007; Kumar et al., 2016). It reflects the net effect of all hydrological flux variables and is, therefore, a key variable for the water cycle and climate change (Long et al., 2017). In 2002, the launch of the twin Gravity Recovery and Climate Experiment (GRACE) satellites (Tapley et al., 2004) represented a breakthrough in TWS mapping. During the GRACE era, J. L. Chen et al. (2013) showed that the accelerated ice melting over Greenland in approximately 2005 drove eastward pole drifting based on GRACE data. Adhikari and Ivins (2016) further pointed out that TWS also plays an important role in decadal-like polar motion: TWS explains ∼ 83% of the peak-to-peak amplitude of polar motion during the period of 2003–2015, and a shift in the global wet-dry pattern dominated a westward drift of polar motion in approximately 2012.

Although climate-driven polar motion became clear during the GRACE era, retrieval of reliable climate-related information on polar motion prior to the GRACE era remains a challenge. Unlike better developed solid Earth, atmosphere, and ocean models, producing actual TWS estimates is a major problem. Many approaches have been proposed to reconstruct long-term TWS (Alkama et al., 2010; J. L. Chen & Wilson, 2005; Hirschi & Seneviratne, 2017; Hirschi, Seneviratne, & Schar, 2006; Hirschi, Viterbo, & Seneviratne, 2006; Mueller et al., 2011; Seneviratne et al., 2004; Tang et al., 2010). Scanlon et al. (2018) reported that TWS trends from hydrological models deviated significantly from the trends obtained from GRACE. Other alternative approaches for estimating TWS are thus valuable in revealing climate-driven polar motion. This study proposes linking the secular polar drift in the 1990s with climate change by testing contributions from simulated TWS under different assumptions. The specific steps include (1) separating the contributions of other excitation sources, that is, the GIA, oceans, and atmosphere, to confirm whether these sources drive secular polar drift in the 1990s; (2) determining the dominant role of TWS in the residual using GRACE and its follow-on mission (GRACE-FO) data; and (3) reconstructing a lengthened TWS based on GRACE/GRACE-FO and reanalysis data and a corrected TWS based on observed glacial change data to assess whether global TWS changes, especially accelerated ice melting over glacial areas, were most likely responsible for polar drift in the 1990s.

2 The Directional Change in Secular Polar Drift in the 1990s

The contribution of multiple excitation sources to polar motion can be quantified by comparing excitations derived from observations and models. Figures 1a and 1b show contributions from the modeled GIA, oceans, and atmosphere to polar motion from January 1981 to June 2020. The geodetic (observed) and geophysical (atmospheric and oceanic) excitations were obtained from the International Earth Rotation and Reference Systems Service (IERS). When manually selecting parameters, the Chandler period of the observed excitation was set as 433 days; the Chandler quality factor was set as 100; the simulation of atmospheric mass redistribution and movement was set based on the European Centre for Medium-Range Weather Forecasts (ECWMF); and the simulation of ocean current and ocean bottom pressure change was set based on the Max Planck Institute Ocean Model (MPIOM). Given that the excitation from the solid Earth is mainly attributed to GIA, this study also takes its effect into account. It is approximately equal to the linear trends of 0.79 mas/yr for χ₁ and −2.95 mas/yr for χ₂ (Adhikari & Ivins, 2016). On a time scale that is longer than the annual and Chandler wobbles, the excitation components χ₁ and χ₂ are approximately equal to the polar motion components X and negative Y (Adhikari & Ivins, 2016). Therefore, the polar drift in this study is computed by applying a 12- and 14- month moving average to excitations. Its direction and rate are thus the moving direction and rate of the linear fitting line of polar drift on the Earth's surface. Polar motion observations show that the Earth's pole drifted generally toward ∼40°W at a rate of ∼2.95 mas/yr during January 1981 and June 2020. Among the three sources of polar motion shown in Figure 1, only the GIA significantly drives the secular polar drift, and it accounts for ∼35% and ∼157% of the variance contributions for χ₁ and χ₂, respectively. Adding atmospheric and oceanic contributions only changes the variance contributions by an additional ∼−14% and ∼−5%. The residuals of the oceanic, atmospheric, GIA-related, and observed excitations (more data processing details are provided in Text S1), which explain ∼81% and ∼−60% of the variance contributions for χ₁ and χ₂, are thus crucial for closing the excitation function budget. Since GIA is assumed to change linearly under the influence of the late Quaternary ice age, the directional change in polar drift in the 1990s is expected to remain in the residuals after subtracting oceanic, atmospheric, and GIA-related excitations. The Breaks for Additive Season and Trend algorithm (Verbesselt et al., 2010) is used to detect breakpoints of the residual because it can effectively detect changes in the long-term trend of the time series in the context of strong periodic changes (Watts & Laffan, 2014). Figures 1c and 1d show that there are breakpoints in χ₁ in approximately October 1998 and in χ₂ in approximately October 1993 for the residual, which is roughly consistent with the conclusions of R. Gross and Poutanen (2009) (in the mid-1990s) and Roy and Peltier (2011) (in approximately 1992). Considering the different results of the two components, the breakpoint of polar drift in the residual is set in October 1995. This study also finds that only for shorter time intervals, for example, 1995–2011, the breakpoints for the two components detected in approximately 2005, similar to the findings of J. L. Chen et al. (2013), which are November 2006 and October 2006 (Figures S2 and S3). The long-term path of observed excitation from 1981 to 2020 further reveals that there is a more significant change in the direction of polar drift that occurred in approximately 1995 (Figure S4). In October 1995, the direction of polar drift in the residual changed from 77°W to 28°E, which is roughly similar to that of the entire observed change, that is, eastward for approximately 71° (from 85°W to 14°W).

3 Comparison of the Excitations of the Residual and GRACE/GRACE-FO-Derived TWS

After the contributions by the GIA, atmosphere, and ocean have been removed, most of the changes in polar motion are likely attributed to changes in the hydrosphere (Adhikari & Ivins, 2016; J. L. Chen et al., 2013; Liu et al., 2018). This study thus further compares the residual and the hydrological excitation at the regional to global scales to help determine the specific drivers of the common directional change in both the secular polar motion and its residual. Limited by insufficient hydrological observations, the comparison is first conducted over the GRACE/GRACE-FO periods. The GRACE/GRACE-FO-derived TWS anomaly data, which are processed (Watkins et al., 2015) at the Jet Propulsion Laboratory (JPL) using the mascon approach (Release 06 Version 02), are employed in this study. These data guarantee that the same atmospheric, oceanic, and GIA simulations used in producing the residual are used to extract TWS from the original GRACE/GRACE-FO data (more data processing details are provided in Text S2). Finally, the GRACE/GRACE-FO-derived TWS data are further converted to the excitation domain as the hydrological excitation (more details are provided in Text S3) to enable the comparison with the residual.

Figure 2 presents the comparison results for the period of April 2002 to June 2020. The monthly hydrological excitation is highly consistent with the residual; the correlation coefficients are as high as 0.92 and 0.82, with the peak-to-peak amplitude accounting for 70% and 66% for χ₁ and χ₂, respectively (Figures 2a and 2b). The TWS-related polar drift, which moves toward 22°E, is also in good agreement with that of the residual (28°E) in direction (Figures 2c and 2d), and the movement rate of the TWS-related polar drift (2.58 mas/yr) accounts for 65% of that of the residual (3.95 mas/yr). According to the study of Adhikari et al. (2018), the discrepancy in direction (6°) and, in particular, rate (1.37 mas/yr) is reasonable when mantle convection is not considered in this study.

Moreover, there is a pronounced difference in the hydrological excitations between polar regions (Greenland and Antarctic) and nonpolar regions. The hydrological excitation in the polar regions is almost linear in time for both components during the study period (Figures 2a and 2b). The mean net trends of the hydrological excitation in the polar regions are 2.78 mas/yr for χ₁ but only −0.43 mas/yr for χ₂. The TWS change in polar regions has driven the pole toward 9°W for approximately 45 mas in 219 months, which mostly explains the moving distance of TWS-related polar drift toward χ₁ (Figures 2c and 2d). In contrast, the hydrological excitation in the nonpolar regions has stronger periodic changes for both components and has driven the pole toward 83°E for approximately 31 mas in 219 months. TWS changes in nonpolar regions mostly excited the residual toward χ₂ in terms of both the drift direction and peak-to-peak amplitude (73%) (Figures 2b and 2d). Previous studies (Adhikari & Ivins, 2016; J. L. Chen et al., 2013) have also found that the TWS in nonpolar regions explains the residual excitation better for χ₂ than for χ₁, and the preliminary explanation is that global lands are mostly distributed along the χ₂ direction. Here, we quantify the contributions from polar regions to χ₁ and from nonpolar regions to χ₂. Considering that the directional change occurred in both components, it is reasonable to infer that TWS probably underwent a significant change in both the polar regions and the nonpolar regions, for example, rapid ice melting or a wet-dry pattern shift across the continents, in the 1990s.

4 Causes of Polar Drift in the 1990s: Accelerated Ice Melting

This section investigates whether TWS underwent a global change and caused a change in polar drift in the 1990s. First, this study assumes that the spatiotemporal characteristics, for example, the rates of ice melting and groundwater depletion, of GRACE/GRACE-FO-derived TWS change are stationary throughout the entire period of 1981–2020. Then, this assumption is tested by explaining the residual.

The assumption is made by lengthening GRACE/GRACE-FO-derived TWS using an improved empirical orthogonal function (EOF) (Yu et al., 2018) analysis. This analysis method decomposes the GRACE/GRACE-FO-derived TWS into several fixed spatial modes and corresponding time series. Given that spatial modes are stationary in time under this assumption, lengthening only requires extending the corresponding time series. An initial TWS estimate throughout the entire study period of 1981–2020, which is the summation of the snow and soil water from ERA5-Land reanalysis, is introduced to help extend the time series. By applying the EOF analysis to this initial TWS estimate, an initial longer-term time series is obtained. The extended GRACE/GRACE-FO-related time series are then computed by correcting this initial longer time series based on the linear relation between time series derived from GRACE/GRACE-FO and ERA5-Land. By multiplying the fixed spatial modes of GRACE/GRACE-FO by the extended GRACE/GRACE-FO-related time series, the lengthened TWS is obtained. More details about lengthening are provided in Text S4. It exhibits the same seasonal fluctuation patterns and linear trends as the GRACE/GRACE-FO-derived pattern at a regional to global scale (Figure S5).

The lengthened TWS is converted into the excitation domain to explain the residual. The hydrological excitation based on the lengthened TWS is similar to the residual between October 1998/1993 and June 2020 (Figure 3). Their correlation coefficients are 0.93 and 0.84 with peak-to-peak amplitudes accounting for 78% and 58% for χ₁ and χ₂, respectively. The consistency with the residual of the lengthened TWS from the 1990s to 2020 is comparable to that of the GRACE/GRACE-FO era, which reveals a high accuracy. Moreover, from the 1990s to present, the lengthened TWS is also significantly correlated with the water level observed by satellite altimetry in four giant lakes (East Aral Sea, West Aral Sea, Lake Ontario and Lake Ladoga) (Figure S6). The surface areas of these lakes are at least larger than one grid cell (∼12,321 km²) of TWS data to guarantee that TWS is mainly affected by lakes in these grids. However, the lengthened TWS cannot explain the residual before the 1990s. The linear trend of the residual before October 1998 (and 1993) is −0.19 mas/yr (and −0.12 mas/yr) for χ₁ (and χ₂), when that of the hydrological excitation based on the lengthened TWS is 2.93 (and 1.26). Beyond the difference in the linear trend, the peak-to-peak amplitude after trend removal is also significant. Before October 1998 (and 1993), hydrological excitation based on the lengthened TWS only accounts for 46% (and 23%) of the peak-to-peak amplitude of χ₁ (and χ₂), which decreased by 22% (and 34%) compared with that after October 1998 (and 1993). The accuracy of the lengthened TWS is challenged between 1981 and the 1990s: that is, the assumption of the TWS lengthening is challenged during this period. The EOF decomposition of GRACE/GRACE-FO is not stationary in time and the spatiotemporal characteristics of the actual TWS probably changed in the 1990s.

To identify specific regions that underwent a significant change and drove polar drift in the 1990s, Figure 4 shows the impacts of the regional TWS change on polar motion during the GRACE/GRACE-FO era. Several areas, for example, Alaska, Greenland, the Southern Andes, Antarctica, the Caucasus, the Middle East, India, and the North China Plain, mainly contributed to the residual during this period. These areas are recent hotspots of TWS change, for example, glacial melting and groundwater depletion. Most areas, such as Alaska, Greenland, the Southern Andes, Antarctica, the Caucasus, and the Middle East, have recently exhibited significant glacier mass changes. According to Zemp et al. (2019), Alaska, Greenland, and the Southern Andes are the world's top three areas of glacier mass loss. Figure 4 also shows that the glacier mass change corresponds to the TWS change after the 1990s. The difference in magnitude can be reasonably explained by the difference in data sources and cover ranges (Text S5 and Table S1). For example, the glacier data cover only ∼1 × 10⁵ km² for Greenland, while the TWS data cover ∼1 × 10⁶ km², which causes the spatial mean glacier mass change to be one order of magnitude greater than that of TWS. In addition, the observational coverage of glacier mass change data only ranges from 1% to 79% of the total glacial area per region. Therefore, the TWS trend of these areas probably mainly attributed to glacier mass change. A previous study (Thomas et al., 2006) suggested that the ice loss in Greenland accelerated sharply after the mid-1990s. According to observed glacier mass change data, the other areas, that is, Alaska, the Southern Andes, Antarctica, the Caucasus, and the Middle East, also experienced sharply accelerated glacier retreats in the 1990s (Figure 4). The lengthened TWS thus cannot capture this rate of change in the 1990s because our hypothesis overestimates the ice melting rates of these areas before the 1990s.

This study attempts to correct the lengthened TWS from January 1981 to December 1995 based on the observed glacier mass change data from Zemp et al. (2019). First, a specific gain factor is computed using the linear regression method for each major glacier area from 1995 to 2018. The input explanatory variable is the yearly observed glacier mass change data, and the input dependent variable is the yearly TWS averaged over all grids that contributed significantly to polar motion (greater than 10⁻³ mas/yr) within this region (Figure 4). To correct the lengthened TWS over these regions, the lengthened TWS time series of the grid contributing more than 10⁻³ mas/yr to polar motion must be multiplied by the region-specific gain factor from January 1981 to December 1995 (more details in Text S5). After correction, the decreasing tendencies of the grids that contribute significantly to polar motion are slowed or altered (Figure S7). Figure 4 shows that the regional TWS change is in better agreement with the observed glacier change after correction. By converting the corrected TWS into the excitation domain, the excitation of TWS to polar motion over these grids also weakens (Figure S8), and the hydrological excitation fits much better with the residual throughout the entire study period after the correction (Figure 3). The linear trends of the hydrological excitation before the 1990s are 0.21 mas/yr and 0.46 mas/yr after correction, which are closer to those of the residual (−0.19 and −0.12) than the uncorrected trend (2.93 and 1.26) (Figure 3). These experiments support the previous suggestion that accelerated ice melting in glacial areas is the main driver of the observed polar drift after the 1990s. In conclusion, the accelerated TWS decline after the 1990s in Alaska, Greenland, the South Andes, Antarctica, the Caucasus and the Middle East, which decreases by an additional approximately −48 mm/yr compared with that before the 1990s, mainly drives a rapid polar drift toward 26°E at a rate of 3.28 mm/yr.

5 Potential of Polar Motion for Revealing Past Global TWS Change

Accelerated ice melting in glacial areas cannot explain the entire polar drift in the 1990s, especially for peak-to-peak amplitudes. The unexplainable part might be attributed to the TWS change in nonpolar regions since it contributes the most to periodic variation. According to previous studies, evidence exists for changes in large-scale climatic conditions in approximately 1995. For example, Kwon et al. (2007) showed that a climatic shift in the summer circulation over East Asia increased the amount of precipitation. Shi et al. (2007) and Y. Chen et al. (2009) reported an enhanced water cycle in Northwest China because of rapid global warming, increased precipitation and runoff. Sutton and Dong (2012) found that a substantial shift in the European climate resulted in a wetter northern Europe and a dryer southern Europe. J. Chen et al. (2002) observed that an increase in the strength of the tropical general circulation caused moist conditions in equatorial convective regions and dry conditions in equatorial and subtropical subsidence regions. These climate changes might alter the spatial pattern of the actual TWS change in nonpolar regions. In addition to possible climate change, Figure 4 shows that the areas of groundwater depletion, for example, northern India and the Middle East, also significantly drive polar drift, at least during the GRACE/GRACE-FO era. The unsustainable consumption of groundwater for irrigation and other anthropogenic activities has been identified as the likely cause of an observed significant decline in the TWS across northern India during the GRACE era (Rodell et al., 2009). This study investigates groundwater withdrawals to determine the existence of a tipping point in Indian TWS changes in the 1990s. The annual groundwater withdrawals of Indian in 1989, 1990, 2000, and 2010 were 194, 229, 235, 341, and 351 billion cubic meters, respectively (Text S6). It is evident that the annual groundwater withdrawals also rose significantly in 2000, that is, after the 1990s. Moreover, the TWS change in the Caucasus and Middle East cannot be totally explained by the glacier mass change. The TWS change is two orders smaller than the glacier mass change when the covered area of the TWS data is three orders larger than that of the glacier mass data. If spreading the glacier mass change evenly over the covered area of the TWS data, it is still one order smaller than the TWS change. This could be related to the incomplete data collection of glacier mass change in the study of Zemp et al. (2019), but previous studies have also mentioned that the TWS decline in the Middle East has resulted from the human exploitation of groundwater. Zavialov et al. (2003) showed that agricultural diversions and unsustainable use of water resources has led to rapid shallowing of the Aral Sea, which began in the 1960s. Voss et al. (2013) further pointed out that 60% of the TWS decline in the Middle East from 2003 to 2009 was due to groundwater depletion. Therefore, the TWS change in the 1990s in the Middle East very likely includes groundwater.

However, the TWS in nonpolar regions always changes over space and time. Compared with polar regions that are dominated by nearly linear changes, the TWS change in nonpolar regions, which could even be a wet-dry pattern shift, is more complex and more difficult to reconstruct and correct. Therefore, although this study suggests a change in TWS across nonpolar regions, the actual change in TWS in these regions before the 1990s cannot be deduced given the limited observations and inadequate hydrological models. Nevertheless, this study offers a clue that the residual could be an alternative factor to assess the accuracy of long-term TWS estimates on a global scale. It should also be noted that, for the same TWS change (1 mm/grid), the regions at low latitudes and high latitudes have a smaller effect on polar motion than those in the mid-latitudes (Figure S9). Thus, polar motion is more sensitive to TWS changes in mid-latitude areas than in other areas.

6 Conclusions

A directional change in polar drift occurred in the 1990s. Given the limited information on global TWS changes before the GRACE/GRACE-FO era, even though there is a hypothesis that accelerated ice melting under global warming is the likely cause, the explanation of the physical process is still a challenge. This study aims to illustrate the primary driver of the polar drift change in the 1990s. A comparison of the geodetic excitation and the geophysical excitations from the atmosphere, oceans, and GIA reveals that only the residual has the same change in the linear trends for both components during the 1990s. The residual is mostly attributed to TWS change. A further investigation into the specific TWS change that mainly drove polar drift after the 1990s is then realized by explaining the residual using two sets of modeled TWS changes in different scenarios, that is, the lengthened GRACE/GRACE-FO-derived TWS and the corrected TWS associated with observed glacier change. The comparison between the residual and the hydrological excitation based on the lengthened TWS shows that the lengthened TWS only adequately explains the residual after the 1990s. A tipping point in the TWS change is thus observed in the 1990s. By correcting the lengthened TWS with glacier change observations, the deviation between the residual and hydrological excitation before the 1990s has been well addressed, with a trend difference decreasing by 87% and 58% for χ₁ and χ₂. The corrected TWS has a slight impact on polar drift before the 1990s but drives a rapid polar drift toward 26°E at a mean rate of 3.26 mas/yr after the 1990s, which is in excellent agreement with that of the residual. The faster ice melting under global warming was the most likely cause of the directional change of the polar drift in the 1990s. The other possible causes are TWS change in non-glacial regions due to climate change and unsustainable consumption of groundwater for irrigation and other anthropogenic activities. The polar drift in 1995 discovered in this study, along with the polar drift in 2005 and 2012 discovered in previous studies, reinforces the suggestion that “TWS is the most plausible causal mechanism for the decadal-like oscillation” (Adhikari & Ivins, 2016).

Erratum

In the Conclusions section, a typographical error regarding a unit of measurement has been fixed from 3.26 mm/yr to 3.26 mas/yr. The present version may be considered the authoritative version of record.