Analysis of the interannual variability in satellite gravity solutions: detection of climate modes fingerprints in water mass displacements across continents and oceans

This study analyzes the interannual variability of the water mass transport measured by satellite gravity missions in regard to eight major climate modes known to influence the Earth’s climate from regional to global scales. Using sparsity promoting techniques (i.e., LASSO), we automatically select the most relevant predictors of the climate variability among the eight candidates considered. The El Niño–Southern Oscillation, Southern Annular Mode and Arctic Oscillation are shown to account for a large part the interannual variability of the water mass transport observed in extratropical ocean basins (up to 40%) and shallow seas (up to 70%). A combination of three Pacific and one Atlantic modes is needed to account for most (up to 60%) of the interannual variability of the terrestrial water storage observed in the North Amazon, Parana and Zambezi basins. With our technique, the impact of climate modes on water mass changes can be tracked across distinct water reservoirs (oceans, continents and ice-covered regions) and we show that a combination of climate modes is necessary to explain at best the natural variability in water mass transport. The climate modes predictions based on LASSO inversions can be used to reduce the inter-annual variability in satellite gravity measurements and detect processes unrelated with the natural variability of climate but with similar spatio-temporal signatures. However, significant residuals in the satellite gravity measurements remain unexplained at inter-annual time scales and more complex models solving the water mass balance should be employed to better predict the variability of water mass distributions.


Introduction
Climate variability exerts profound influences on the water cycle, and therefore society. Over the past decade, increased evidence of an intensification of the water cycle has been reported (e.g., Eicker et al. 2016), leading to more frequent or more severe climate extremes, such as droughts, floods and wildfires (e.g., AghaKouchak et al. 2014;Turco et al. 2014;Yoon et al. 2015;Mora et al. 2018;Barichivich et al. 2018;Stevens-Rumann et al. 2018). Linking observed changes in climate extremes to anthropogenic greenhouse gases emissions is incredibly difficult due to the limited time-span of historical records and natural climate variability, responsible for large fluctuations in the climate system over a broad range of time and space scales (e.g., Seneviratne et al. 2012;Easterling et al. 2016). In order to adapt to new risks and resources in a changing climate, it is necessary to develop tools able to characterize the natural variability of the climate system and pinpoint the mechanisms triggering changes observed in the climate system.
Climate modes synthesize emergent characteristics of the climate system in various regions of the world. These repeating patterns of the climate variability are usually identified through the statistical analysis (e.g., regional average, empirical orthogonal functions (EOF), principal component analysis) of one or several climate fields, such as sea surface temperature, air temperature, air pressure, sea level pressure, sea surface height, precipitation, wind speed and wind direction. Teleconnections, caused by the transport of moisture and energy through the atmosphere and the ocean, allow propagating the influence of 1 3 regionally defined climate modes over large distances. Notoriously, the El Niño-Southern Oscillation (ENSO), observed by definition in the tropical Pacific, has global impacts on the Earth's climate, ecosystems and society (e.g., McPhaden et al. 2006;Capotondi et al. 2015). Similarly, the Atlantic Multidecadal Oscillation (AMO), has been reported to have nearly global influences on the climate system (e.g., Knight et al. 2006), including Amazonian (e.g., Kayano et al. 2016) or Sahelian (e.g., Mohino et al. 2011) rainfall, Atlantic hurricanes (e.g., Zhang and Delworth 2006), North American (e.g., Hu et al. 2011) and European summer climate (e.g., Zampieri et al. 2017). More and more frequently, it is found that a single climate mode is not sufficient to explain the major characteristics of regional climates, but that a combination of several modes, interacting with one another, is necessary to gain in climate predictability (e.g., Morrow et al. 2010;Li and Wettstein 2012;Kayano and Capistrano 2014;Petrick et al. 2014;Kundzewicz et al. 2019). Besides, it has been suggested that climate modes may not be stationary and that the spatial and temporal characteristics of climate modes, as well as their relationships with one another, may be evolving in a changing climate (e.g., Litzow et al. 2020).
The GRACE (Gravity Recovery And Climate Experiment) mission (e.g., Tapley et al. 2004), and its successor GRACE Follow-On (GRACE-FO) (e.g., Landerer et al. 2020), enabled global monitoring of the mass transport within the Earth's system, leading to unprecedented advances in our understanding of the global water cycle, ice and ocean mass variability in a changing climate (e.g., Tapley et al. 2019). Because they take into account all water bodies from the top of the atmosphere to deep aquifers, GRACE measurements are used to develop reliable indicators of the climate variability, including almost-continuous monitoring of water shortage (e.g., Thomas et al. 2017) and/ or excess (e.g., Chen et al. 2010) since 2002. Numerous studies have shown the influence of climate modes on GRACE and GRACE-FO measurements, but are usually restricted to regional (e.g. Becker et al. 2011;Petrick et al. 2014;Zhang et al. 2015;Ndehedehe et al. 2017;Vissa et al. 2019;Xie et al. 2019) or continental (e.g. de Linage et al. 2013;Anyah et al. 2018;Ndehedehe and Ferreira 2020;Liu et al. 2020) scale studies, and/or to a limited subset of climate indices (e.g. Matsuo and Heki 2012;Bergmann and Dobslaw 2012;Landerer and Volkov 2013). Global analyses of the interannual variability of GRACE data are mostly focused on the impact of ENSO on terrestrial water storage (e.g., Phillips et al. 2012;Boening et al. 2012;Ni et al. 2018;Hamlington et al. 2019). Almost all studies previously cited acknowledge however the importance of the interactions between climate modes and recommend broader scale studies looking at the impact of several climate modes on GRACE and GRACE-FO measurements to better identify the sources of natural variability (e.g., Phillips et al. 2012;Petrick et al. 2014, Ni et al. 2018, Hamlington et al. 2019. The objectives of this study are: (1) identify the climate modes driving significant water mass changes at interannual time-scales and, (2) quantify the water mass changes due to climate modes. Correcting GRACE and GRACE-FO measurements for the contributions of climates modes could in turn unveil weaker signals masked by the natural variability of the climate system, such as deep Earth signals related to mass changes in the fluid core or at the core-mantle boundary (e.g., Mandea et al. 2015). To detect the climate modes responsible for water mass variations at interannual time scales, we carry a LASSO (Least Absolute Shrinkage and Selection Operator: Tibshirani (1996)) regression of eight climate indices, namely MEI (Multivariate ENSO index), PDO (Pacific Decadal Oscillation), NPGO (North Pacific Gyre Oscillation), AMO, NAO (North Atlantic Oscillation), AO (Arctic Oscillation), SAM (Southern Annular Mode) and IOD (Indian Ocean Dipole). Indeed, it was recently suggested that the LASSO regularization can disentangle the contributions of various climate modes in historical reconstructions of steric sea levels to provide simple, easily interpretable models of the natural climate variability (Pfeffer et al. 2018). We test this method on GRACE and GRACE-FO satellite gravity measurements to estimate the contributions of major climate modes to the interannual variability of the water mass distributions. This analysis is the first global climate mode study to include both oceans and continents in the analysis of GRACE and GRACE-FO measurements.

Physical links between climate modes and the global water cycle
The natural variability of the climate system gives rise to a vast array of dynamical climate modes, such as ENSO, AMO, NAO or SAM. Climate modes describe statistically emergent spatio-temporal patterns of regional climate variables, such as the sea/air temperature, precipitation, wind speed and sea/air pressure. Changes in temperature, precipitation and wind speed drive evaporation of water from the oceans and its transport through the atmosphere. During the particularly intense episode of La Niña in 2010/2011, the global mean sea level fell by about 5 mm due to the transport of water from the oceans to the continents, leading to particularly devastating flood events in Australia (e.g., Boening et al. 2012). This event illustrates how climate modes can influence key elements of the global water cycle (precipitation, evapotranspiration and runoff), leading to typical changes in the water storage, summarized with the mass balance equation (Eq. 1): where dS/dt are the water storage variations integrated from the surface to the aquifer, P is the precipitation, ET is the evapotranspiration and R is runoff. Wind forcing can also become extremely significant over oceans, especially at high latitudes and shallow seas, where large mass transport occurs at interannual time-scales (e.g., Piecuch and Ponte 2011;Piecuch et al. 2013). By setting the rhythm of precipitation, evaporation, transpiration, runoff and wind speed, climate modes can exert a profound influence on the variability of water mass distribution across the ocean, atmosphere and continents, which are expected to be observable with GRACE and GRACE-FO satellite gravity missions.

Satellite gravity solutions
Monthly mass variations have been estimated globally using the latest GRACE (Tapley et al. 2004(Tapley et al. , 2019 and GRACE-FO (Landerer et al. 2020) mascon solutions produced by the CSR (Center for Space Research) and JPL (Jet Propulsion Laboratory). In this study, mascon solutions have been preferred over spherical harmonic solutions, because of a better localization of the origin of mass variations and reduced leakage errors. In turn, mascon solutions use an inversion strategy that cannot be easily reproduced or replicated, relying on specific a priori constraints and regularization procedures. While mascon solutions allow less liberty than spherical harmonics in the choice of background models (e.g., Blazquez et al. 2018) and noise reduction techniques (e.g., Horvath et al. 2018;Kusche et al. 2009), they ensure replicable results from one study to another, as long as the same versions of the same products are used.
In this study, we used the CSR and JPL Level 3 Release 6 mascon solutions, consisting of monthly mass grids spanning from April 2002 to July 2020 at the time of the download (September 2020). All mass changes are assumed to occur within a thin layer at the Earth's ellipsoid surface (WGS84) and are expressed in equivalent water heights. The JPL and the CSR provide global mascon solutions, using 3° spherical caps  and 1° equal area elements . Because the CSR and JPL solutions are developed, processed and archived in a shared Science Data System (SDS), they harmonize several necessary processing steps. The degree-1 coefficients are estimated using the methods from Sun et al. (2016), built upon the work of Swenson et al. (2008) suggesting that geocenter motions could be modelled using a combination of GRACE measurements and ocean models. The degree-2, order-zero coefficients (C 20 and after August 2016, the C 30 coefficients (1) dS dt = P − ET − R as well) are replaced with the solutions from Satellite Laser Ranging provided by the GSFC (Goddard Space Flight Center) (Loomis et al. 2019). After August 2016, an accelerometer transplant was carried out to compensate for the loss of one of the two accelerometers (Bandikova et al. 2019). The AOD1B RL06 model (Dobslaw et al. 2017a) is used as a background model for the dealiasing of the ocean and atmosphere. In the CSR and JPL solutions, the effect of atmospheric loading on the ocean surface has been restored using the GAD correction (Dobslaw et al. 2017b) and Glacial Isostatic Adjustment (GIA) has been corrected using the ICE6G-D model (Peltier et al. 2018). Besides, the Coastal Resolution Improvement (CRI) filter was applied on JPL mascon solutions, to reduce leakage in coastal areas ). We did not apply land grid scaling (i.e., gain factors developed for the interpretation of hydrological signals in mascon solutions at higher resolutions; Wiese et al. 2016) to the JPL or CSR mascon solutions, because they are not adapted for fully global analyses including the oceans, continents and ice-covered regions. The CSR and JPL Level 3 Release 6 mascon solutions were downloaded from the PODAAC website (https:// podaac. jpl. nasa. gov), with all corrections listed above applied (Cooley and Landerer 2019). In order to convert ocean bottom pressures to ocean mass variations, we removed the GAD from the CSR and JPL mascon solutions and restored the GAB, accounting for the ocean circulation only. To facilitate the comparison of the two solutions, the data have been converted into freshwater anomalies, using a water density of 1000 kg m −3 and interpolated on a regular 1° × 1° grid using a nearest neighbor interpolation. Both CSR and JPL anomalies refer to the 2004-2009.99 time average. Because this study focuses on interannual signals, a secular trend and seasonal signals have been removed at each grid point, using the outputs of an ordinary least squares regression taking simultaneously into account a constant, linear trend, annual sinusoid and semi-annual sinusoid.

Climate indices
A significant part of the climate variability tends to take place in distinctive patterns repeated in space and time known as climate modes. Climate modes are usually defined as statistically prominent features of the regional climate variability, usually determined by empirical orthogonal functions (EOF) of one or several climate variables in a given region. Here, we consider eight climate indices: the MEI, PDO, NPGO, AMO, NAO, AO, SAM and IOD (Fig. 1).
ENSO is a major actor of the global climate, triggering extreme events such as cyclones, floods, droughts or heat waves in various regions of the world (e.g., Santoso et al. 2017). The MEI is a multi-variable index used to  (Wolter and Timlin 1998). It is defined as the leading component of the combined EOF analysis of sea level pressure (SLP), sea surface temperature (SST), zonal and meridional components of the surface wind, and outgoing longwave radiation (OLR) over the tropical Pacific basin (30° S-30° N and 100° E-70° W).
The PDO consists in a decadal SST oscillation, characterized by a basin-wide dipole dividing the North Pacific along a typical northeast-southwest chevron pattern (Mantua et al. 1997). It is defined as the leading EOF of SST anomalies north of 20° N in the Pacific region (Newman et al. 2016) and estimated with the ERSSTV5 dataset (Huang et al. 2017). The PDO is closely related to the ENSO and is sometimes referred to as a "long-lived ENSO", to characterize the longer periods of the climate index.
The NPGO expresses as a double gyre in the North Pacific (Di Lorenzo et al. 2008). It is defined as the second mode of sea surface height variability in the Northeast Pacific (180° W-110° W; 25° N-62° N). Connected with ENSO and PDO, it contributes to global climate variability, which is evidenced by global sea level trends and sea surface temperature (Di Lorenzo et al. 2010).
The AMO consists in a multidecadal SST oscillation across the whole Atlantic basin, with cool and warm phases that may last for 20-40 years (Enfield et al. 2001). We use the unsmoothed index, provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, from their Web site at https:// psl. noaa. gov and calculated from 1948 to present as the area weighted average of the Kaplan SST V2 dataset (Kaplan et al. 1998) over the North Atlantic (0-70° N).
The NAO index reflects the difference in surface sea level pressure between the Subtropical Azores High and the Subpolar Low (Lamb and Peppler 1987). The NAO impacts the intensity and location of the North Atlantic jet streams and storm tracks (Woollings and Blackburn 2012). It modulates heat and moisture transport, resulting in temperature and precipitation changes that can extend from North America to central Europe (Hurrell et al. 2003). The monthly NAO index is calculated as the first mode of a Rotated Principal Component Analysis (RPCA) applied on monthly mean 500 millibar height anomaly data from 1950 to 2000 over 0-90° N of latitude (Barnston and Livezey 1987).
Strongly connected with NAO (Zhou et al. 2001), the AO is a climate pattern characterized by winds circulating counterclockwise around the Arctic at around 55°N latitude (Baldwin and Dunkerton 1999). During positive AO phases, cold air masses tend to be confined at high latitudes by strong circular winds. During negative AO phases, cold air masses can more easily move southwards, increasing the storminess at mid-latitudes (Higgins et al. 2000). The AO index is defined as the leading EOF mode of monthly mean 1000 millibar height anomaly data from 1979 to 2000 over 20° N-90° N. SAM, also known as the Antarctic oscillation (AAO), is a major driver of climate variability in the Southern Hemisphere (Thompson and Solomon 2002). It has been shown to influence windiness, storminess, rainfall and temperatures from the subtropics to Antarctica (Gillett et al. 2006). Defined as a zonal pressure difference between the latitudes of 40° S and 65° S (Marshall 2003), it reflects the contraction (in a positive phase) or expansion (in a negative phase) of the westerly winds belt circling Antarctica.
Strongly connected with ENSO, the IOD consists in an interannual oscillation of the SST in the equatorial Indian Ocean, forming a typical west-east dipole (Saji et al. 1999). The Dipole Mode Index (DMI) is defined as the difference of SST averaged in the western (50° E-70° E and 10° S-10° N) and eastern (90° E-110° E and 10° S-0° N) equatorial Indian Ocean, and calculated with the HadISST1.1 SST dataset with the 1981-2000 climatology removed.
To be consistent with the processing of satellite gravity solutions, all climate indices have been detrended and deseasoned. Decadal trends and seasonal signals have been removed using the outputs of an ordinary least squares regression of a linear trend, annual sinusoid, semi-annual sinusoid and a constant carried out for the GRACE/GRACE-FO time period (April 2002-September 2020). Besides, to facilitate the comparison of climate modes, all indices have been normalized, so that their time average equals zero and their standard deviation equals one.
The MEI, PDO, AMO and DMI (IOD index) are distributed by the NOAA (National Oceanic and Atmospheric Administration) PSL (Physical Science Laboratory) climate data repository (https:// psl. noaa. gov/ data/ clima teind ices/). NAO and AO are distributed by the NOAA Climate Prediction Center (https:// www. cpc. ncep. noaa. gov/). The NPGO index is available through the NPGO webpage at http:// www. o3d. org/ npgo/ npgo. php. The SAM index is available from Gareth Marshall's webpage: http:// www. nerc-bas. ac. uk/ icd/ gjma/ sam. html. All web pages listed above were accessed for the last time on the 14th June 2021.

LASSO regression
Our inversion procedure consists in a suite of time series analyses, carried out independently from one grid cell to another. To assess the relationships between climate modes and water mass changes, we consider the following regression model (Eq. 2): where h(t) is the change in equivalent water height estimated by GRACE/GRACE-FO at a given location, C i is one of the eight climate indices selected in this study (PDO, MEI, NPGO, NAO, AMO, AO, SAM, DMI), the coefficients β i are the parameters to be estimated, β 0 is a constant and ϵ are the residuals of the regression (i.e., the changes in equivalent water height that are not taken into account by our regression model). A positive β i expresses a positive correlation between the observed signal h(t) and climate indices C i (t), while a negative β i expresses an anti-correlation. Because it is unlikely for all eight indices to contribute simultaneously to the water mass changes measured at a given location, either on land or over the oceans, the model can be considered as sparse. At a given location, the interannual variability of the water mass changes observed by GRACE and GRACE-FO only depends on a subset of climate indices. This subset is expected to change with the geographical location. The LASSO is a regularization technique allowing the selection of relevant features in a multiple variable regression model (Tibshirani 1996). The LASSO minimizes the sum of squared errors, with an upper bound on the sum of the absolute values of the model parameters called the LASSO constraint. The minimization function Ψ is then (Eq. 3), Because of the LASSO constraint, some of the coefficients of the regression model will be set to exactly zero, leading to more interpretable models with less explanatory variables.
The degree of shrinkage of the coefficients β i depends on the value of the penalty λ (Fig. 2b). If λ = 0, the regression becomes an ordinary least squares regression. On the other hand, if λ is too high, all explanatory variables will be eliminated from the regression model, as all β i coefficients will be set to zero (Fig. 2b). The penalty λ is therefore calculated to minimize the standard error of a fivefold cross validation (called CV below, see Fig. 2c.) (e.g., Arlot and Celisse 2010), i.e., each time series is randomly divided in five parts, successively used as testing (one part) and training (four parts) datasets. Detailed information about the search of optimal LASSO penalty through cross-validation procedures can be found in the Sects. 3 and 7 of Hastie et al. (2009). The minimization of the standard error of a K-fold cross-validation has been shown to work well in the presence of correlated variables (Hebiri and Lederer 2013). This procedure is similar to the approach of Pfeffer et al. (2018), who extracted climate mode fingerprints for the PDO, ENSO, NPGO, AMO, IOD and IOBM (Indian Ocean Basin Mode; Yang et al. 2007) in historical reconstructions of steric sea levels.
This procedure allows selecting the most relevant climate indices to reproduce interannual water mass changes at each geographical location sampled by the GRACE/GRACE-FO mascon solutions. However, because of the LASSO constraint, the amplitude of coefficients diminishes (Tibshirani 1996). To recover the optimal values of the regression model coefficients (β i ), an ordinary least squares (OLS) regression is carried out at each geographical location, using a simpler regression model, including only the climate indices selected with the LASSO. This final step avoids the shrinkage of the regression model coefficients (β i ), but keeps the selection capacity of the LASSO (Table 1). We tested our inversion procedure with a synthetic test. Synthetic data were generated using the regression model (Eq. 1) with β i coefficients equal to 10, 7.5 and 5 for MEI, NAO and DMI and 0 for PDO, NPGO, AMO, AO, SAM. We added noise to the synthetic data, using a random sample from a normal distribution centered on zero with a standard deviation of 10. This level of noise is quite strong, as it is the same amplitude of the synthetic ENSO signal, and is significantly larger than the NAO (β NAO = 7.5) and IOD (β IOD = 5) signals. We chose these coefficients to (1) test whether the LASSO is able to set the PDO, NPGO, AMO, AO and SAM coefficients to exactly zero and (2) Table 1. As expected, a classic OLS inversion is not able to select the relevant predictors of the water mass variability, as all coefficients are non-null (Table 1). However, the amplitudes of MEI, NAO and IOD coefficients are close to synthetic values (Table 1). With a LASSO regression, the value of the β i coefficients diminishes with the penalty λ (Fig. 2b). Optimal values of the penalty can be estimated with a CV procedure, estimating the CV error (i.e. the mean squared error evaluated over each test set, in grey in Fig. 2c) as a function of the penalty λ. The penalty minimizing the average CV error (i.e. the CV error averaged over the five folds, in black Fig. 2c) is not sufficient to set the PDO value to zero (Table 1; Fig. 2b). This is likely due to the fact that the PDO is correlated with MEI (online resource Fig. OR9). The standard error, defined as the average CV error plus one standard deviation (see for example Fig. 2.3 in Hastie et al. 2019), takes into account the dispersion due to the random division of the dataset during the CV procedure. The penalty minimizing the standard error is therefore higher, and is able to select the relevant predictors of the regression model (MEI, NAO and IOD) and set all other parameters to exactly zero. The final OLS regression carried out on the selected predictors succeeds in recovering the full amplitude of ENSO, NAO and IOD signals (Table 1).

Climate mode fingerprints
Climate mode fingerprints are defined as the spatial distribution of β i coefficients (Eqs. 2 and 3) resulting from the inversion of climate indices with a suite of LASSO and a final least squares regressions. The climate mode fingerprints of the JPL (Fig. 3) and CSR (Fig. 4) are extremely similar, though the CSR climate mode fingerprints seem to be sparser (more null values), especially across arid areas (Sahara, central Asia, central Australia). For comparison, a classical EOF analysis of the water mass changes observed by GRACE and GRACE-FO and the inferred correlations with climate indices are presented in the online resource (OR).
ENSO has a strong impact on the interannual variability of both continental and oceanic water mass changes. Negative ENSO anomalies (Figs. 3b and 4b) indicate a water shortage during El Niño phases and a water excess during La Niña. For example, negative anomalies in Australia correspond to wet Niña and dry Niño conditions (in red in Figs. 3b and 4b), having the potential to drive extreme rainfall and floods during La Niña (e.g., Evans and Boyer-Soucher 2012;Boening et al. 2012) and disastrous droughts and wildfires during El Niño (e.g., Mariani et al. 2016;Wang and Cai 2020). The largest negative ENSO anomalies are located in West Antarctica (down to − 120 mm), in the North of the Amazon basin (down to − 80 mm) and in the Gulf of Carpentaria (down to − 60 mm). A large (10-15 mm) positive anomaly lies across the Southern Pacific and the Southern Ocean, which was also detected in the EOF analyses of the JPL and CSR solutions (online resources, Figs. OR1, OR2 and OR3). This increased ocean mass (positive anomaly in blue) during positive phases of ENSO is coincident with enhanced ice-mass loss (negative anomaly in red) in West Antarctica, and might be associated with the transport of hot water from the tropical Pacific to the Southern Ocean during El Niño events, which has been shown to accelerate the basal melting of ice-shelves (e.g. Paolo et al. 2018) and may by extension impact grounded ice (e.g. Joughin et al. 2021). This is further discussed in Sect. 4.2. In general, ENSO seems to have a large impact on the interannual variability of water mass distributions in shallow seas (Gulf of Carpentaria and around Indonesia) and continental water storage globally (e.g., North of the Amazon basin, Central America, Okavango and Zambesi basins, Northeast of Australia, Penner basin, Eastern Europe, West Siberia, West coast of Canada, Winnipeg Lake, East coast of USA).
The PDO fingerprint (Figs. 3a and 4a) is correlated with the MEI fingerprint (Figs. 3b and 4b), with some noticeable differences, confirming the tight but complex relationships between the two modes. The PDO mostly impacts land water storage, despite potential influences over the South Pacific and Southern Ocean of large geographical extent but extremely small amplitude (less than 5 mm). The PDO resembles MEI in the North East of the Amazon, over Northwest America, in Africa, North East Australia, large parts of Eastern Europe and the North of Asia. The PDO differs from MEI in the Parana basin, across South East Asia and over the oceans in general. This tends to show that the PDO and MEI have complementary mass transport fingerprints, both representative of ENSO, with correlated but distinct periodicities (see the PDO and MEI times series in Fig. 1 and their spectral content in the online resource: Fig. OR8).
The NPGO fingerprint (Figs. 3c and 4c) adds another element to the complex ensemble of Pacific teleconnections. This regional mode, recently detected in the variability of sea surface heights in the Pacific (Di Lorenzo et al. 2008), has, rather unexpectedly, a lot of common variability with land hydrology over the GRACE/GRACE-FO time period. Some of it might be due to fortuitous correlations of water mass changes with the climate index, especially over arid regions such as the Sahara, the Arabian Peninsula or central Asia, where we do not expect much hydrological signal (Fig. 3c). Remarkably, the CSR and JPL fingerprints differ for the NPGO over the above-mentioned regions. The CSR displays a sparser NPGO fingerprint, with values that are generally null for arid areas (Fig. 4c). The large majority of spatial patterns is however consistent between the two centers. The NPGO fingerprint also seems to be anticorrelated with MEI, over the Zambezi basin, in Northeast Australia, in Northeast America and in central Europe (Fig. 4c). We do expect a physical impact of the NPGO around the North Pacific basin and such signatures are indeed present in Alaska, along the West coast of Canada and in California (Figs. 3c and 4c). AMO (Figs. 3d and 4d) influences the interannual variability of land water storage along the Amazon River, North of California, the extreme South of Alaska and the North of British Columbia, the Middle East and eastern Mediterranean, in between the Ob and Yenisei rivers in Siberia, across Myanmar and Thailand and in the Zambezi River basin. The impact of AMO on interannual ocean mass changes seems to be negligible. It might be noted that the impact of AMO on the interannual hydrological variability remains strong in the Amazon River basin, even for extremely high penalties (not shown here), which makes it extremely robust.
Unlike the PDO, NGPO and AMO fingerprints, AO (Figs. 3e and 4e), NAO (Figs. 3f and 4f) and SAM (Figs. 3 and 4g) have clear and robust influences on the ocean mass transport, but do not impact the interannual variability of the continental water storage. As expected, the AO influences interannual water mass transport over the Arctic, the North  (Figs. 3e and 4e). The AO also influences water mass transport in the Mediterranean Sea (Figs. 3e and 4e). NAO is significantly correlated with AO (R > 0.6 in online resource Fig. OR9). As a consequence, AO seems to have been favored over NAO, and the NAO fingerprint is almost null everywhere, except over the Arctic and North Atlantic (Figs. 3f and 4f), which intensifies the impact of AO (Figs. 3e and 4e). The SAM fingerprint exhibits a clear bipolar anomaly across the Southern Ocean, positive near the coast of Antarctica, and negative away from the Antarctic coastline (Figs. 3 and 4g).
The IOD (Figs. 3 and 4h) seems to have a rather limited impact on the interannual variability of the water mass transport both on land and in the ocean. Positive and negative anomalies are found in the JPL and CSR solutions in southeastern Europe, in California, in the Pacific and Indian Ocean and in the Java Sea.

Selection capacity of the LASSO
We introduced a LASSO regularization to improve the interpretability of our model through the selection of relevant predictors of the climate variability at interannual time scales. The climate mode fingerprints resulting from the LASSO inversion are sparse; for each pair of geographical coordinates most of the β i coefficients (Eqs. 2 and 3) are set to zero. One way to look at this is through the L0-norm, indicating the number of non-zero β i coefficients in the regression model. The L0-norm can range from zero (all climate indices are eliminated from the regression model) to eight (all climate indices are selected in the regression model).
For the JPL and CSR solutions (Fig. 5), the LASSO regularization, coupled with a 5-fold cross validation minimizing the standard error, performed efficient variable selection on a global scale. The L0-norm is set to zero for large parts of the Atlantic Ocean, Indian Ocean and Pacific Ocean. Interannual ocean mass changes are indeed expected to be negligible, except for high-latitude regions and shallow seas (Piecuch and Ponte 2011;Piecuch et al. 2013). SAM, MEI and AO contribute significantly to the interannual variability of ocean mass in the Arctic, Southern Ocean, Mediterranean Sea, Northern Pacific and shallow seas in the North of Australia and around Indonesia (Figs. 3 and 4), leading to non-zero L0-norms in these regions (Fig. 5). We would also expect the L0-norm to be close to zero in arid areas, where not much climate variability is observed. For the CSR, the L0-norm is indeed equal to zero over the Sahara, Arabian Peninsula, central Asia and South West Australia. For the JPL, the outputs of the LASSO regression seem noisier, with non-null coefficients over the same arid areas.
Besides, the L0-norm should not be too high in any region of the world. It seems unlikely for all climate modes to contribute to the interannual water mass transport observed in one place. Some of the climate modes are correlated (such as PDO, MEI and IOD or AO and NAO) and there are some known teleconnections between them. It would therefore seem realistic for 3 or 4 climate modes to contribute together to the interannual variability observed in one place, but more than that seems unlikely. Looking at high values of the L0-norm can therefore give us an indication of the presence of noise in the climate mode fingerprints. Our method is rather successful in this regard on a global scale (Fig. 5), though there is some noise remaining in the outputs of the regression that can be identified in rare locations by high L0-norm values (values above 5 or 6).
High L0-norm values are not found in the same places for the CSR and JPL solutions. For the CSR, high L0-norm values are found near the coast, especially between Greenland and Canada, suggesting noisier solutions in these regions, potentially due to unresolved leakage issues. This issue is not reproduced in the JPL solutions, maybe due to the application of the Coastline Resolution Improvement (CRI) filter ). However, we found high L0-norm values for the JPL solutions in semi-arid to arid areas such as central Asia, West Africa and Southern California (Fig. 5). This may be due to different regularization strategies used to lessen the noise in the mascon solutions, based on geophysical models for the JPL  and on an iterative procedure weighting the signal to noise level for the CSR . If the hydrological model (GLDAS NOAH; Rodell et al. 2004) used to constrain terrestrial water storage is unsatisfactory across semi-arid and arid areas, it might let noise in the JPL mascon solutions that could be retrieved in the results of our inversion.

Performance of climate modes predictions
The model performance can also be assessed using determination coefficients (R 2 ), interpreted as the fraction of variance explained by the regression model compared to the variance measured in locally detrended and deseasoned gravity solutions (Fig. 6). The RMS (root mean square) of locally detrended and deseasoned CSR and JPL solutions ( Fig. 6a and b) ranges from a few millimeters in the tropical Atlantic, to a few centimeters on land, to more than half a meter across fast thinning glaciers located in Greenland and Antarctica.
According to our regression model, the eight climate modes considered in this study contribute less than a centimeter to the interannual variations observed in tropical ocean basins with amplitudes of less than 2 cm. R 2 values reach at most 30% in the tropical Pacific ( Fig. 6e and f), due to the combined effect of the PDO, MEI and SAM (Figs. 3 and  4). In extra-tropical ocean basins, climate modes (RMS up to 25 mm) contribute for about 20-50% to the interannual variability observed in the CSR and JPL solutions (Fig. 6). It is due to the combined effect of SAM and ENSO in the Southern Ocean and NAO and AO in the Arctic (Figs. 3 and  4). In shallow seas, climate modes contribute for several centimeters (RMS up to 7 cm in the Gulf of Carpentaria) to the observed interannual variability (RMS up to 10 cm in the Gulf of Carpentaria), which explains from 20 to 70% of the observed signal (Fig. 6). In the Gulf of Carpentaria, most of the climate variability is explained with ENSO, and for the Mediterranean Sea a significant part (from 30 to 40%) of the interannual variability is explained with AO (Figs. 3 and 4).
The largest contributions of climate modes to interannual terrestrial water storage changes (excluding glaciated regions) are found in the Zambezi, North Amazon and Parana basins (Fig. 6), where the combined effect of the PDO, ENSO, NPGO and AMO (Figs. 3 and 4) generates interannual signals with an RMS reaching 130, 110 and 90 mm respectively. The effect of combined climate modes explains up to 60% of the interannual terrestrial water storage variability in the Zambezi and North Amazon basins and 40% in the Parana basin. Significant interannual variations (with RMS ranging from 30 to 80 mm) in the terrestrial water storage are also generated by climate modes in North Australia, Indonesia, South East Asia, India, along the Tigris and Euphrates rivers, across the Caspian Sea, in the Ob and Yenisei basins, in Eastern Europe, in central America, along the Gulf of Mexico (in Florida in particular), in California, over the Great Lakes and in Canada. The PDO, ENSO, NPGO and AMO are responsible for most of the interannual terrestrial water storage changes (with R 2 values ranging from 20 to 70% depending on the region considered), though AO and NAO were found to contribute around the Mediterranean Sea, and the IOD was found to have a rather surprising influence in Central Europe. However, for several large regions of the world, such as central South America, central North America, West Africa, tropical Africa (including the Congo basin), North India, South East China and North West Australia, significant inter-annual variability of the terrestrial water storage is observed with GRACE and GRACE-FO measurements, that failed to be reproduced with our simple regression model. This is likely due to the occurrence of other processes that may be climate related, but are not adequately represented by climate modes. Geophysical processes (e.g., earthquakes) or anthropogenic influences (e.g., land use change, groundwater pumping) are also known to influence the variability of the mass transport measured by GRACE and GRACE-FO, and are not taken into account by climate modes.
Ice mass changes over glaciated regions (e.g., West Antarctica, Greenland, Alaska, Svalbard and the Andes) also display significant interannual variability, with RMS values ranging between 100 and 500 mm after removal of a linear trend, an annual sinusoid and a semi-annual sinusoid 1 3 Fig. 6 Performance of the regression model: (a, b) Amplitude of the interannual signals the for the CSR and JPL mascons solutions respectively; (c, d) Amplitude of the climate modes predictions for the CSR and JPL solutions respectively; (e, f) Determination coefficients between the climate modes predictions and the locally detrended and deseasoned CSR and JPL solutions respectively ( Fig. 6a and b). Climate modes contribute actively to the inter-annual ice mass variability, by enhancing/reducing precipitation and/or bringing cold/hot air and ocean masses near the glaciers, leading to RMS values ranging between 200 and 500 mm in Svalbard, 100 and 250 mm in Greenland, 50 to 200 mm in West Antarctica and in Alaska. Surprisingly AO and SAM seem to have a negligible effect on the interannual variability of ice mass changes (Figs. 3 and 4), while the PDO, ENSO, NPGO and AMO explain between 30 and 50% of the interannual variability in Greenland and West Antarctica (Fig. 6e and f) and up to 70% in Alaska ( Fig. 6e and f).

Implications for climate applications
Climate modes contribute significantly to the interannual variability in the global water cycle. The GRACE and GRACE-FO measurements offer the unique opportunity to monitor the water mass displacements through the atmosphere, the ocean and the continents. During positive phases of ENSO (El Niño events), we can detect a water mass uptake (negative anomaly in red) in intertropical regions (i.e., Gulf of Carpentaria, shallow seas around Indonesia, North East Australia, North Amazon and Zambezi basin) transported to high latitudes, and resulting in a positive anomaly (in blue) across the Southern Pacific and Southern Ocean. This transport of warm waters during positive phases of ENSO might in turn be responsible for enhanced ice melting (red anomaly) in West Antarctica (in particular along the Pine Island and Thwaites glaciers) and at the Southwest of the Antarctic Peninsula (Figs. 3b and 4b). The opposite happens during negative phases of ENSO (La Niña events), resulting in significant interannual variability (Figs. 6c and d). ENSO signatures have been detected previously in GRACE observations of icemass changes in the same region (Kaitheri et al. 2021). This is in agreement with the study by Paolo et al. 2018, showing that the net ice shelf mass decreases during intense El Niño events in the Antarctic Pacific sector due to basal melting. While the melting of ice-shelves can not be detected by GRACE, it has been shown to significantly enhance glacier flow, leading to significant ice-mass loss on the continent (e.g. Rignot et al. 2019). Such ice dynamics happen over a broad range of time-scales spanning over decades (e.g. Thomas et al. 2004), but also include instantaneous speed-up responses to changing boundary conditions (e.g. ocean warming, grounding line retreat) (e.g. Joughin et al. 2021). Besides, the mass balance of inland glaciers is also impacted by precipitation, that has been shown to contain intermittent ENSO and strong SAM signatures (e.g. Sasgen et al. 2010;Kim et al. 2020). Numerous processes triggering ice-mass changes in West Antarctica are therefore susceptible to contain ENSO signatures. Further studies including multivariate analyses of the ocean temperature, precipitation and state of glaciers (height, mass, discharge) are necessary to determine the mechanisms leading to ENSO signatures in West-Antarctic glaciers ice-mass changes.
The impact of ENSO, depicted with the MEI index and fingerprint, has to be combined with other climate modes such as the PDO, NGPO, AMO and SAM. Indeed, in the Southern Ocean, the combined effect of ENSO (Figs. 3b and 4b) and SAM (Figs. 3 and 4g) contributes moving water masses and accelerating/decelerating ice melting in West Antarctica (e.g., Paolo et al. 2018). More generally, interannual water mass transport both on land and in the ocean depends on complex wind regimes, precipitation and temperature conditions related with several climate modes (e.g., Petrick et al. 2014;Ni et al. 2018;Ndehedehe and Ferreira 2020;Liu et al. 2020). This is evidenced, for example in South America, where a combination of AMO, PDO and ENSO is necessary to explain the interannual variability of the water storage along the Amazon River, in the North of the Amazon basin, along the Atlantic Brazilian coast and within the Parana basin (Figs. 3a,b,d and 4a,b,d). AMO has indeed been shown to interact with ENSO in the Amazon basin, through the exchange of moisture in key tropical regions (e.g., García-García and Ummenhofer 2015).
Besides, it has been widely shown that ENSO is better depicted with a large variety of climate indices, including but not limited to MEI and the PDO, to catch the variety of timescales (from biennial to decadal) and spatial patterns associated with this global climate mode, influencing both the sea level (e.g., Zhang and Church 2012) and terrestrial water storage (Fasullo et al. 2016;Luo et al. 2016). This is clearly shown in the correlation of MEI and PDO fingerprints that are both needed to describe ENSO. More unexpectedly, the NPGO seems to complement the MEI and PDO fingerprints, and contributes significantly to interannual water mass displacements around the globe (Figs. 3a,b,c and 4a,b,c). The three Pacific teleconnections should therefore be interpreted together. Several studies confirm the emerging influence of the NPGO in recent decades (Di Lorenzo et al. 2008;Ye et al. 2016, and its association with the PDO and ENSO (Di Lorenzo et al. 2010;Zhang et al. 2013). The changing nature of the PDO, ENSO and NPGO could also be associated with global warming. The non-stationary nature of climate modes is evidenced, among others, by increased association of the NPGO with the first mode of climate variability in the North Pacific usually associated with the PDO (Litzow et al. 2020). The GRACE and GRACE-FO record being extremely limited in time (here July 2002-August 2020) would therefore only offer a temporary snapshot of the interactions between the water mass displacements and climate modes.
Remarkably, while the PDO, MEI, NGPO and AMO influence the interannual variability of terrestrial water storage; NAO, AO and SAM mostly impact interannual ocean mass changes in extratropical ocean basins (i.e., the Arctic and Southern Ocean) and shallow seas (i.e., the Mediterranean Sea, Gulf of Carpentaria and the ensemble of shallow seas bordering Indonesia). SAM has a clear impact on the Southern Ocean circulation (Figs. 3 and 4g, largely due to the wind regime, which is in line with the findings of Bergmann and Dobslaw (2012), Makowski et al. (2015) and Liau and Chao (2017). As NAO and AO are highly correlated, AO was favored over NAO and explained a large part of the ocean mass variability in the Arctic, North Atlantic and Mediterranean Sea. This is in good agreement with the findings of Peralta-Ferriz et al. (2014), Peralta-Ferriz et al. (2016) and Landerer and Volkov (2013). It also evidences the capacity of the LASSO to choose among several correlated predictors and select the most relevant to predict the observations considered.
Climate modes are therefore important to predict the interannual hydrological and ocean mass variability. However, for large regions of the world, our model is unable to explain the GRACE/GRACE-FO observations. Other mechanisms should be invoked that are not likely to be well represented with a simple empirical model, but rather require to solve the water mass balance with global hydrological models. In future studies, it could also be useful to investigate the possibility of a lagged response of the water mass displacements to climate modes, which has been shown to be significant for various regions of the world (Phillips et al. 2012;Kundzewicz et al. 2019).

Implications for geodetic applications
GRACE and GRACE-FO measurements offer a unique view of the global mass transport within the whole Earth system at interannual scales. The model suggested in this study, based on a suite of robust LASSO inversions of multiple climate indices, provides an efficient way to automatically correct the satellite gravity solutions for relevant modes of the climate variability. Though climate modes often explain a marginal part (less than 50%) of the interannual equivalent water mass transport, the predicted signals reach amplitudes that can be easily confounded with other sources of mass transport and mask weaker processes of geodetic interest. Over extratropical ocean basins and shallow seas, the ocean mass transport due to climate modes can reach several centimeters at inter-annual scales. In large hydrological basins, such as the Amazon, the Zambezi or the Parana, the interannual variability of the terrestrial water storage predicted by climate modes reaches about 10 cm. Over glaciated regions, climate modes generate interannual variations in the ice mass reaching up to 50 cm.
Such signals are significantly larger than the uncertainty on GRACE and GRACE-FO solutions (e.g., Riegger et al. 2012;Scanlon et al. 2016;Ferreira et al. 2016;Blazquez et al. 2018) and need to be corrected in order to unveil small amplitude signals originating, for example, from anthropogenic climate change (e.g., Fasullo et al. 2016) or deep Earth interior processes (e.g., Mandea et al. 2012Mandea et al. , 2015. The increase of global mean ocean mass due to freshwater inputs reaches only 2.3 +/− 0.19 mm/ year over the Jan. 2005-Dec. 2016 time-period , with inter-annual variations of about +/− 5 mm that can be easily masked by climate modes (Hamlington et al. 2019). Trends in total terrestrial water storage range from 0 to +/− 20 mm/year (Rodell et al. 2018), with interannual variations of about +/− 50 mm (Tapley et al. 2019). Changes in the ocean circulation can lead to inter-annual variations in the water mass transport of about +/− 20 mm . Because GRACE and GRACE-FO measurements provide an integrated measurement of the temporal changes in the gravity field, independent models are needed to identify the sources of mass variability. The model designed in this study could therefore be used as a global correction for eight different climate modes leading to significant mass changes in various regions of the world. Special attention has been paid to the robustness of the correction, taking only into account significant predictors of the water mass variability. However, the correction would have a limited scope as the numerous processes contributing to the interannual variability of the climate system cannot all be represented with climate modes.
Besides, climate mode fingerprints provide a framework to compare interannual signals from different solutions, allowing to link subtle choices in the processing strategy (e.g., the regularization, filters applied to reduce leakage) to relatively simple parametric results (i.e., a limited number of maps representative of different climate modes).

Conclusions
The measurements of the GRACE and GRACE-FO satellite gravity missions have now reached enough maturity and precision to allow the detection of typical oscillations of the water mass distributions related to climate modes through (1) the enhancement/reduction of precipitation, evapotranspiration and runoff, (2) transport of moisture across the atmosphere and (3) wind forced ocean circulation. Climate modes contribute significantly to the interannual variations in ocean mass (up to 25 mm in extratropical-basins and up to 70 mm in shallow seas), terrestrial water storage (up to 100 mm in large hydrologic basins such as the Amazon, Zambezi or Parana) and ice thickness (up to 500 mm in equivalent water heights). While climate modes have already been detected in satellite gravity observations by several independent research teams, this study is the first to account for the teleconnections of eight major climate modes in the global water cycle.
Over the oceans, much of the interannual variability is observed in extratropical basins and shallow seas. The principal actors (among the climate modes selected) of the interannual ocean variability are ENSO (nearly all extratropical basins and shallow seas), SAM (Southern Ocean) and AO (Arctic, North Atlantic, Mediterranean Sea).
Over the continents, a large part of the interannual variability of the terrestrial water storage can be attributed to four climate modes, namely AMO, ENSO, PDO and NGPO. When combined together, these four modes explain up to 60% of the interannual storage variability observed in the North Amazon, the Zambezi Basin and up to 40% in the Parana basins, where climate modes predictions generate RMS values up to 10 cm.
In other parts of the world (grey areas in Fig. 6e and f), our climate modes predictions are unable to fully explain the interannual variability of terrestrial water storage observed by GRACE and GRACE-FO measurements. This may be due to the fact that we did not consider a lagged response of the terrestrial water storage to climate modes (e.g., Phillips et al. 2012) or more likely due to the influence of other climatic, anthropogenic and hydrologic processes that require more sophisticated models able to solve the water mass balance depending both on climate conditions and land use changes.
One important result of this study lies in our ability to track interannual water mass displacements across different reservoirs. For example, we can link the transport of water from intertropical regions (both on land and in the tropical Pacific) to the Southern Ocean, where it contributes to the interannual variability of ice mass changes in West Antarctica in connection with ENSO, SAM and PDO. Another significant finding stresses the need to combine several climate indices in order to be able to predict interannual changes in water mass displacements. The variations in terrestrial water storage depend indeed on several factors (wind, temperature, precipitation) influencing the climate regime that are better represented by a combination of climate modes than by a single one. We also found that the NPGO is an important driver of interannual water mass transport globally, which may be due to a recent association of the NPGO with the PDO and ENSO (Di Lorenzo et al. 2010). Such patterns might be transitory, and it is expected that some of the relationships evidenced here will not last because of global warming and global climate change (Litzow et al. 2020).
Our findings show the necessity to maintain our climate observation system, including satellite gravity measurements, in order to extend the length and scope of instrumental records used to diagnose the state of our water resources and predict its future evolution. Finally, the present study will be of great help for better separating, in the GRACE/ GRACE-FO data, the climate-related signals occurring in the surface fluid Earth's envelopes from deep Earth signals due to fluid motions in the outer core and to exchange of matter at the core-mantle boundary (Mandea et al. 2012(Mandea et al. , 2015.