Geophysical Signal Detection in the Earth’s Oblateness Variation and Its Climate-Driven Source Analysis

: This study analyzes the geophysical and 1.0%, as well as 25.4% and 0.7%, respectively. GRE, ANT, and GLA have ~3 to ~7-year periodic ﬂuctuations and a positive linear trend, excluding GIA.


Introduction
The accurate recovery and analysis of the Earth's oblateness term C 20 (or J 2 ) time series are of great significance for understanding the mass redistribution in the Earth's climate system [1,2]. Since 1976, the C 20 time series derived from the satellite laser ranging (SLR) observations to geodetic satellites by Cheng and Ries [2] from the Center of Space Research (CSR) has been of the highest quality, even though various research institutes have recently taken their efforts to the SLR-derived, GRACE-derived (Gravity Recovery and Climate Experiment satellite mission) or SLR-GRACE-derived C 20 time series, such as the Astronomical Institute, University of Bern (AIUB) [3], the Deutsches Geodätisches Forschungsinstitut, Technische Universität München (DGFI) [4], the NASA Goddard Space Flight Center (GSFC) [5] and the Centre National d'Etudes Spatiales/Groupe de Recherche de Géodésie Spatiale (GRGS) [6]. GRACE is less sensitive to the zonal harmonics, in particular J 2 [7], because of its near-polar orbit inclination (89.3 • ). Though the GRACEderived C 20 has been greatly improved with the improvement of solution strategies and force models in recent years, it is not as reliable as SLR estimates [7] and should be replaced by SLR-derived C 20 for scientific interpretations and geophysical applications of GRACE and GRACE-FO data products. In this paper, the time-variation characteristics of the SLR-derived C 20 time series from CSR [1,2] are intensively analyzed.
During the past four decades, various analyses of SLR data have indicated that the mass redistribution within the Earth's dynamic system has been undergoing significant variations as a response to the tidal and non-tidal changes ranging from hours to decades or longer-periods [2,8], which is characterized by the temporal variations in the Earth's gravity field. The surface mass redistribution concentrated in a thin layer on the Earth's surface is reflected by the variations of the spherical harmonic (SH) coefficients of the Earth's gravity field model. C 20 is the largest SH coefficient, which is often described instead of J 2 with a relationship of J 2 = − √ 5C 20 . J 2 is called the Earth's oblateness and its variation is at the magnitude of 10 −10 .
In the previous studies, the variations of J 2 have been routinely characterized by the following model It is the superposition of a secular trend and accelerations ( . J 2 and .. j 2 ), tidal (∆J T 2 (t)) and non-tidal variations (∆J NT 2 (t)) [1,2]. ∆J T 2 (t) and ∆J NT 2 (t) are typically described by seasonal variations, interannual variations, decadal variations, and longer-period variations [2,9]. Since 1976, Cheng and Ries [2], Cheng et al. [1], Cheng and Tapley [10], and other researchers [9,[11][12][13][14][15] have studied the J 2 time series as its time-span increases year by year, understanding the nature of the secular trend of J 2 , from a linear trend by analyzing a short time series [10,16,17] or a long time series of more than three decades employing a quadratic trend [1,2,7,9,13,15]. Cheng et al. [1] characterized the nature of the secular trend by a quadratic curve imposed on a period tidal signal of 18.6 years. The quadratic term is described by the superposition of a linear decrease induced by glacial isostatic adjustment (GIA) or postglacial rebound (PGR) of the Earth's mantle [11,18,19] and a slowly increasing rate caused by global mass redistribution related to the ice-sheet and glacial melting as well as mass transfer in the atmosphere and ocean [2]. Other effects are also considered as factors influencing the oblateness changes, such as co-seismic effects induced by earthquakes [20,21], artificial water reservoir impoundment [22], and the Earth's spin-down, induced by the tidal forcing [23] which contributes to the deceleration of J 2 at the order of 10 −12 /yr. The global sea-level rise enlarges the above rate by almost the same value. The different secular decrease rates of the Earth's oblateness have been reported by analyzing the SLR data covering different time intervals from 5 years up to four decades, where the first secular decrease is −3 × 10 −11 /yr [12], then −2.6 × 10 −11 /yr derived from 10-year data after 1980 and −2.8 × 10 −11 /yr from 1976 to 1995 [24], −1.95 × 10 −11 /yr from 1980 to 1998 [15], −2.75 × 10 −11 /yr from 1976 to 2004 [10], −5.9 × 10 −11 /yr from 1976 to 2010 by considering it to be a quadratic curve [1].
In addition to the well-known annual and semiannual variations of J 2 , the interannual variations related to the strong EI Nino-Southern Oscillation (ENSO) [10,25] generally with the fluctuations of 4 to 6 years were reposted. "The 1998 anomaly" described by Cox and Chao [15] has drawn significant attention and is believed to be caused by the mass redistribution associated with the atmosphere, ocean, and land water [1]. Besides the interannual period, the decadal variations of 18.6 years and about 10.5 years are the evident signals from the analysis of Cheng and Ries [2]. The period of 18.6 years is a pure tidal harmonic signal, while the possible cause and nature of the variation of about 10.5 years with variable amplitude and phase related to climate change remain unknown. Previous studies [1,2,15] indicated that, apart from the tidal harmonic and linear variations, J 2 variation presents an obvious climate-related nonharmonic behavior, especially the varying amplitude and phase. To process these signals, a variety of methods were used. Cheng and Tapley [1,15] adopted the wavelet analysis to decompose ∆J 2 into high-and low-frequency, where the discrete approximation of the Meyer wavelet is applied as a low pass filter and a sliding average of 360 days is used as a high pass filter, analyzing the long-term variations and removing the seasonal and high-frequency variations in J 2 , respectively. Chao et al. [9] subtracted the quadratic trend and seasonal signals with a constant amplitude from ∆J 2 by using the least-squares fit, ignoring the stochastic behavior of these signals. Cheng and Ries [2] detected the evident variation of 18.6 years in the original ∆J 2 time series by using Meyer wavelet, whereas Ding and Chao [26] and Chao et al. [9] presented a quite small signal of 18.6-year by using the stabilized AR-z spectrum (autoregressive spectrum implemented in the complex z domain) after removing the linear variation and 18.6-year tidal theoretical value. However, Seo et al. [13] failed to find this 18.6-year signal in the residual time series. Cheng and Ries [2] explained that the lack of this signal in the analysis of Seo et al. [13] is attributed to the 11-month sliding average used. Therefore, the existence of the 18.6-year signal from the analysis of different methods has led to different discussions.
Singular spectrum analysis (SSA) is a powerful method for studying nonlinear time series. It was first proposed and used by Colebrook [27] in biological oceanography, then introduced into paleoclimatology by Fraedrich [28] and in the study of social issues by Hassani [29]. It can construct and decompose as well as reconstruct a trajectory matrix based on the time series, realizing the accurate separation of principal components, including trend and periodic signals from the remaining component generally regarded as noise [30,31], even when the time series is incomplete [32,33]. The SSA is able to analyze the structure of the time series and make further predictions. As a powerful non-parametric spectrum estimation technique, SSA allows us to effectively follow the evolution in amplitudes and phases of signals [29,34]. As a special application of the empirical orthogonal function, SSA conducts the construction of the principal components of a time series and its prediction by using the w-correlation criterion [29]. Besides, the Lomb-Scargle (L-S) periodogram is well known for detection and characterization of the periodic signals hidden in an unevenly sampled time series and has been widely used in the analysis of astronomical data. The L-S occupies a unique position in various available methods: it is driven by Fourier analysis [35], but it is equivalent to the least-squares method [36]. It can be generated from the principles of Bayesian probability theory [37] and has a close relationship with the bin-based phase folding technology in some cases [38]. As another method, the L-S periodogram holds a unique correspondence point between many types of methods, so we also attempt to analyze the variations of J 2 by using the L-S periodogram in this study.
In addition to the above spectrum analysis, the time series of ∆J 2 can also be analyzed and reconstructed using geophysical models and other approaches [39]. Many researchers have modeled mass variations of ice sheets, glaciers, terrestrial water, atmosphere, etc., for seeking the geophysical source interpretation in the ∆J 2 time series. Seo et al. [13] modeled the climate contributions to ∆J 2 , accounting for the relation between the glacial mass change over the Arctic and Alaska and the oscillation of 18.6-year tide. Chao et al. [9] added the theoretical terms (i.e., 18.6-year tide) that Cheng et al. [2] has removed back to ∆J 2 and modeled the climate-driven sources of Arctic Oscillation, Pacific Decadal Oscillation (PDO), ENSO, etc. for establishing the connections between these sources and the interannual-todecadal variations of observed ∆J 2 while Cheng and Tapley [10] validated that the 4-to 6year fluctuations are related with the PDO and ENSO by geophysical models. Sun et al. [39] improved a new approach related to fingerprints of the geoid changes patterns to estimate ∆J 2 , explaining the residuals by terrestrial water storage after removal of the atmosphere and ocean variations from ∆J 2 . In this study, we analyze the climate-driven components in ∆J 2 and then conduct the comparison between the reconstructed and observed ∆J 2 time Remote Sens. 2021, 13, 2004 4 of 18 series, which enables us to gain an in-depth understanding of the geophysical interpretation of J 2 variations. Therefore, this study mainly focuses on the analysis and interpretation of Earth's oblateness variations (∆J 2 ) from two different perspectives. One is to analyze the observed ∆J 2 time series from the perspective of spectrum analysis by applying both SSA and L-S periodogram, respectively; the other is to calculate the contributions of climate-driven components to ∆J 2 from the perspective of a geophysical model and then analyze the contribution of each component to the observed ∆J 2 time series and reconstruct the ∆J 2 time series as well as give the corresponding geophysical interpretation. The rest of the paper is organized as follows. Section 2 briefly presents the theory of SSA and L-S periodogram. Section 3 provides the analysis of the secular trend and periodic signals in J 2 variations and their geophysical interpretation. Section 4 gives a brief introduction of the climate-driven components of ∆J 2 and the comparison between the reconstructed and observed ∆J 2 . Conclusions are given in Section 5.

Singular Spectrum Analysis
By using SSA, the original one-dimensional time series can be decomposed into a sum of constituents. Each component can be determined as either a noise component, a quasi-periodic or periodic component, or a trend component [29]. For a J 2 series {J i , i = 1, 2, . . . , M} with length M, the trajectory matrix J can be formed by a window with size W (W < M/2) as follows: Then, the covariance matrix P J , which has a Toeplitz structure, is expressed as where the unbiased autocovariance function a(i) is computed by The covariance matrix P J is decomposed as P J = VΛV T , where Λ is a diagonal matrix with the diagonal values λ q sorted in descending order and V is an orthogonal matrix with row eigenvectors v q . The temporal principal components (PCs) matrix C can be generated as C = VJ after projecting the series J i onto v q . The PCs at kth row and qth column of C is expressed as where 0 ≤ k ≤ M − W, 1 ≤ i ≤ W. The qth Reconstructed Component (RC), which is defined by showing how to optimally extract a series of length M, corresponding to a given set of eigenvalues, can be given by where J k,q is the kth component from the qth PCs. Since λ q is sorted in descending order, the first several numbers of RCs can approximate the original time series, while the remaining RCs can be regarded as noise. Therefore, if one wants to reconstruct seriesĴ q , it can be obtained by the sum of the first q RCs. The window size is a key issue of SSA and is chosen empirically. The optimal window size should maximize the separation of the contained signals [40] and is usually proportional to the periods of main signals. The window length W is selected to best reflect the separability of SSA. The property 'separability' of SSA is also called the weighted correlation or w-correlation [29], which is a natural measure of dependence between two decomposed series T (k) and T (i) in Equation (6), where The smaller absolute w-correlation value means the worse w-orthogonal of two series, and the zero value indicates that the two series are completely separable. On the contrary, the two series will be far from w-orthogonal and strongly dependent. Therefore, we possibly regard the two PCs with large w-correlation as one periodic signal pair.

Lomb-Scargle Periodogram Analysis
The L-S periodogram is a well-known technique to attain periodicity by analyzing the frequency spectrum, especially from unequally spaced data [41,42], which is exactly equivalent to the least-squares fitting of a sine wave. In previous studies, it was mainly developed and used in spectral analysis of the unequally sampled astronomical observation. It has the advantage of maintaining the property of the original data and avoiding the interpolation of uneven data. Here, as a special application, we use it to determine the spectrum of equal data.
For the J 2 series {J k , k = 1, 2, . . . , M} with the length M, the discrete Fourier transform for an arbitrarily sampled data set is defined as Then, the classical periodogram is computed by applying the definition of the Fourier power spectrum, as follows: A generalized form of the periodogram is addressed by Scargle [42] and Zechmeister and Kürster [43] as where τ is defined to ensure time-shift invariance for each t k and is computed by In the least-squares analysis, the shift τ in Equation (10) is considered by Lomb [41] to orthogonalize the normal equations. Because of the strong connection between leastsquares analysis and Fourier analysis, Equation (10) is commonly referred to as the L-S periodogram [43,44].

Geophysical Signal Detection from ∆J 2 Time Series
The J 2 time series from more than 44-year of the satellite laser ranging observations of eight geodetic satellites from 1976 to 2020 derived by Cheng and Ries [2] from the CSR website (http://download.csr.utexas.edu/pub/slr/degree_2/Long_term/C20_Long_ Term.txt, accessed on 20 January 2021) has been employed, where the theoretical value of the 18.6-year tide was deducted. Besides the J 2 time series, the monthly time-variable gravity field up to degree and order 3, plus C 40 , is also co-estimated by using the SLR observations. The J 2 time series and the corresponding uncertainties indicated by the light gray band are shown in Figure 1, exhibiting a superposition of a secular trend and periodic variations and smaller uncertainties after 1992, when observations of two satellites, LAGEOS-2 and Stella, were included.
A generalized form of the periodogram is addressed by Scargle [42] and Zechmeister and Kürster [43] as where τ is defined to ensure time-shift invariance for each k t and is computed by In the least-squares analysis, the shift τ in Equation (10) is considered by Lomb [41] to orthogonalize the normal equations. Because of the strong connection between least-squares analysis and Fourier analysis, Equation (10) is commonly referred to as the L-S periodogram [43,44].

Geophysical Signal Detection from ΔJ2 Time Series
The J2 time series from more than 44-year of the satellite laser ranging observations of eight geodetic satellites from 1976 to 2020 derived by Cheng and Ries [2] from the CSR website (http://download.csr.utexas.edu/pub/slr/de-gree_2/Long_term/C20_Long_Term.txt, access on 20 January 2021) has been employed, where the theoretical value of the 18.6-year tide was deducted. Besides the J2 time series, the monthly time-variable gravity field up to degree and order 3, plus C40, is also co-estimated by using the SLR observations. The J2 time series and the corresponding uncertainties indicated by the light gray band are shown in Figure 1, exhibiting a superposition of a secular trend and periodic variations and smaller uncertainties after 1992, when observations of two satellites, LAGEOS-2 and Stella, were included.

Geophysical Signals Detected by Singular Spectral Analysis
For the 44-year J2 time-series, after the removal of the tidal theoretical value with 542 data points, the selection of the window length W is quite important, which is the only parameter in the SSA technique. On the one hand, if the window length is too short, the periodic signal is not easily separated. On the other hand, since the determination of the appropriate window length depends on preliminary information of the time series and the problem at hand [29], it is generally recommended that the window length should be proportional to the period of the interest but not greater than M/2 [29]. Considering the decadal signals, we mainly focus on signals, such as 10.5 and 18.6 years. Therefore, we select three window lengths with 112 points, 132 points, and 223 points, which correspond to the periodic signals of 9.3, 11.0, and 18.6 years, respectively. To identify periodic signals, we compute the absolute w-correlations of the first 30 reconstructed components (RCs) corresponding to the above three window lengths and show them in Figure 2. The w-

Geophysical Signals Detected by Singular Spectral Analysis
For the 44-year J 2 time-series, after the removal of the tidal theoretical value with 542 data points, the selection of the window length W is quite important, which is the only parameter in the SSA technique. On the one hand, if the window length is too short, the periodic signal is not easily separated. On the other hand, since the determination of the appropriate window length depends on preliminary information of the time series and the problem at hand [29], it is generally recommended that the window length should be proportional to the period of the interest but not greater than M/2 [29]. Considering the decadal signals, we mainly focus on signals, such as 10.5 and 18.6 years. Therefore, we select three window lengths with 112 points, 132 points, and 223 points, which correspond to the periodic signals of 9.3, 11.0, and 18.6 years, respectively. To identify periodic signals, we compute the absolute w-correlations of the first 30 reconstructed components (RCs) corresponding to the above three window lengths and show them in Figure 2. The wcorrelations, shown in Figure 2a,b, are less separable than those of Figure 2c with the window length W = 223. In Figure 2c, the absolute w-correlations of RCs (2,3), RCs (5,6), correlations, shown in Figure 2a,b, are less separable than those of Figure 2c with the window length = 223. In Figure 2c, the absolute w-correlations of RCs (2,3), RCs (5,6), and RCs (7,8) are larger than others, indicating that these pairs are strongly correlated. Moreover, we are also mainly interested in the existence of the 18.6-year signal. Therefore, the window length = 223 is most suitable to separate the various periodic signals in the J2 time-series. The variations in J2 with a window length of 18.6 years are analyzed for the detection of the various periodic signals. The eigenvalues of the first 30 RCs are illustrated in Figure 3 for the window length = 223. The first eigenvalue is much larger than the others, which indicates that RC (1) is much more significant than other components. The first and forth eigenvalues show the general tendency of the J2 time-series better than the first one alone [29]. Therefore, the RCs (1,4) represents the trend with the red solid line, as illustrated in Figure 4. By using the quadratic polynomial fitting to the red solid line (in green dotted line) and directly to the J2 time-series (in black solid line), the derived secular linear decrease rate of about (−5.80 ± 0.08) × 10 −11 /yr and obvious acceleration rate of about (2.38 ± 0.02) × 10 −12 /yr 2 are consistent with the values of −5.90 × 10 −11 /yr and 1.80 × 10 −13 /yr 2 reported by Cheng et al. [1], as listed in Table 1. Moreover, the evolution of the trend of the J2 time-series can be better characterized with the solid red line compared to the black solid line, as shown in Figure 4, indicating that the trend derived by SSA has better performance than that from quadratic polynomial. This quadratic trend is attributed to solid Earth's response to the mass redistribution due to recent significant ice melting [2,5,9,13,45,46].  The eigenvalues of the first 30 RCs are illustrated in Figure 3 for the window length W = 223. The first eigenvalue is much larger than the others, which indicates that RC (1) is much more significant than other components. The first and forth eigenvalues show the general tendency of the J 2 time-series better than the first one alone [29]. Therefore, the RCs (1,4) represents the trend with the red solid line, as illustrated in  Table 1. Moreover, the evolution of the trend of the J 2 time-series can be better characterized with the solid red line compared to the black solid line, as shown in Figure 4, indicating that the trend derived by SSA has better performance than that from quadratic polynomial. This quadratic trend is attributed to solid Earth's response to the mass redistribution due to recent significant ice melting [2,5,9,13,45,46]. correlations, shown in Figure 2a,b, are less separable than those of Figure 2c with the window length = 223. In Figure 2c, the absolute w-correlations of RCs (2,3), RCs (5,6), and RCs (7,8) are larger than others, indicating that these pairs are strongly correlated. Moreover, we are also mainly interested in the existence of the 18.6-year signal. Therefore, the window length = 223 is most suitable to separate the various periodic signals in the J2 time-series. The variations in J2 with a window length of 18.6 years are analyzed for the detection of the various periodic signals. The eigenvalues of the first 30 RCs are illustrated in Figure 3 for the window length = 223. The first eigenvalue is much larger than the others, which indicates that RC (1) is much more significant than other components. The first and forth eigenvalues show the general tendency of the J2 time-series better than the first one alone [29]. Therefore, the RCs (1,4) represents the trend with the red solid line, as illustrated in  Table 1. Moreover, the evolution of the trend of the J2 time-series can be better characterized with the solid red line compared to the black solid line, as shown in Figure 4, indicating that the trend derived by SSA has better performance than that from quadratic polynomial. This quadratic trend is attributed to solid Earth's response to the mass redistribution due to recent significant ice melting [2,5,9,13,45,46].   Besides, the second and third components, the fifth and sixth components, the seventh and eighth components have close eigenvalues at their respective order of magnitude. Thus, RCs (2,3), RCs (5,6), and RCs (7,8) are three pairs of periodic signals. Their characteristics are shown in Figure 4 in detail. Applying the Fast Fourier Transform analysis, we can determine the amplitudes and periods of the three pairs of signals which are the annual, 10.6-year, and semi-annual signals. The annual signal has a much larger amplitude than the 10.6-year and semi-annual signals. As illustrated in Figure 4 and Table 2, the amplitude of the annual signal is almost five times larger than those of 10.6-year and semi-annual periods.

Geophysical Signals Detected by Lomb-Scargle Periodogram
L-S periodogram, as a widely analysis tool used for the astronomical observation, has the simplest statistical behavior, which is important to assess the reliability of the possible  Besides, the second and third components, the fifth and sixth components, the seventh and eighth components have close eigenvalues at their respective order of magnitude. Thus, RCs (2,3), RCs (5,6), and RCs (7,8) are three pairs of periodic signals. Their characteristics are shown in Figure 4 in detail. Applying the Fast Fourier Transform analysis, we can determine the amplitudes and periods of the three pairs of signals which are the annual, 10.6-year, and semi-annual signals. The annual signal has a much larger amplitude than the 10.6-year and semi-annual signals. As illustrated in Figure 4 and Table 2, the amplitude of the annual signal is almost five times larger than those of 10.6-year and semi-annual periods.

Geophysical Signals Detected by Lomb-Scargle Periodogram
L-S periodogram, as a widely analysis tool used for the astronomical observation, has the simplest statistical behavior, which is important to assess the reliability of the possible detection, and is believed to distinguish between the strictly periodic and nonperiodic signals [42,43]. In addition, the study by Reegen [47] indicates that L-S is based on the identification of the highest peak in a chosen frequency interval, while discrete Fourier transform suffers from systematic deviations of peak frequencies [48]. Reegen [47] also elaborates that L-S is more stable against aliasing and shows improved frequency accuracies than Fourier transform. In this study, the L-S periodogram, is applied for the first time to further analyze the Earth's oblateness variation after the removal of the tidal theoretical value. The L-S periodograms of the observed ∆J 2 time series after removing a linear trend and after removing the secular trend with the frequencies from 0 to 5 cycles per year are illustrated in Figures 5 and 6 detection, and is believed to distinguish between the strictly periodic and nonperiodic signals [42,43]. In addition, the study by Reegen [47] indicates that L-S is based on the identification of the highest peak in a chosen frequency interval, while discrete Fourier transform suffers from systematic deviations of peak frequencies [48]. Reegen [47] also elaborates that L-S is more stable against aliasing and shows improved frequency accuracies than Fourier transform. In this study, the L-S periodogram, is applied for the first time to further analyze the Earth's oblateness variation after the removal of the tidal theoretical value. The L-S periodograms of the observed ΔJ2 time series after removing a linear trend and after removing the secular trend with the frequencies from 0 to 5 cycles per year are illustrated in Figures 5 and 6, respectively.  Apart from the significant annual signal, the semi-annual signal and the interannual fluctuations from 4 to 7 years are also detectable in Figure 5. Using this method, a quite small tidal signal with a period of 18.6-year can be detected from the residual time series after removing the 18.6-year theoretical tidal signal and the linear trend from the time series, which was also found by Ding and Chao [26] when the AR-spectrum was analyzed. Unfortunately, the signal of 10.5 years cannot be found, while a broadband signal with a peak of approximately 12 years is visible. The peak of 47.6 years in Figure 5 cannot Apart from the significant annual signal, the semi-annual signal and the interannual fluctuations from 4 to 7 years are also detectable in Figure 5. Using this method, a quite small tidal signal with a period of 18.6-year can be detected from the residual time series after removing the 18.6-year theoretical tidal signal and the linear trend from the time series, which was also found by Ding and Chao [26] when the AR-z spectrum was analyzed. Unfortunately, the signal of 10.5 years cannot be found, while a broadband signal with a peak of approximately 12 years is visible. The peak of 47.6 years in Figure 5 cannot be reasonably explained since it is very challenging to determine whether there is power at wavelengths near the length of the time series (44 years in this paper). For similar reason, no clear explanation can be made by Ding and Chao [26] for the cause of the nominal 56-year peak shown in the time series from 1976 to 2019. The peak of 47.6 years may be from the artificial signal, which is probably attributed to the sensitivity of the used method to the signal hard to be detected and may explain why different methods show varying abilities to detect some particular signals. According to Cheng et al. [1] and as shown in Figure 1, the Lomb-Scargle periodogram is implemented to the ∆J 2 time series after removing the secular trend. The signal of 10.6 years is obvious, and the peak of 47.6-year disappears in Figure 6, indicating that the 47.6-year signal may be part of the secular term. Besides, the disappearance of the 18.6-year signal verified Cheng Minkang's description [1] that the nature of the secular trend is characterized by a quadratic term imposed on a period tidal signal of 18.6 years. In this case, the L-S method does have an advantage over the Fourier transform in detecting the weak 18.6-year signal, which is the same as Chao's research on finding this 18.6-year signal in the de-linear ∆J 2 time series [9]. However, the Fourier transform cannot find this signal.

Discussions
Cheng and Ries [2] detected the dominant period of about 10.4 years from the J 2 time-series, Cheng and Tapley [10] and Seo et al. [13] found a similar period of 10.6 years based on wavelet decomposition, which is determined with the amplitude of 5.01 × 10 −11 in this study using SSA. Seo et al. [13] also showed that solar activity, measured by the number of sunspots, varies for about 11 years, by comparing the modeled J 2 time series and the number of sunspots from 1979 to 2010. Chao et al. [9] suggested a relationship between the solar cycle and climate change by analyzing the wavelet spectrum of the sunspot variations with 10.5 years. Because solar radiation pressure is quite stable in a solar cycle, it is unlikely that solar activity will directly affect the SLR solutions via the variations in solar radiation pressure [13]. The more acceptable explanation may be that solar activity somehow affects the redistribution of atmospheric or water mass by altering climate change or other unknown causes, which are not included in climate models. Consequently, the causes of 10 years variation still need to be further investigated, as well as convincing observations to understand the effects of the solar activities on climate change and the global gravity field variations.
Although the ∆J 2 time series may contain the dominant effects of the 18.6-year tide related to the nutation and Earth-Moon tidal interactions [2] which may be a cause of multiyear variations [10], Cheng and Ries [2] found a significant signal of 18.6 years before removal of the 18.6-year tide from the ∆J 2 time series. Although the residual tidal signal of 18.6 years in the ∆J 2 time series becomes quite small after the theoretical values are removed, it is still expected to be detected in the analysis of Ding and Chao. In this paper, we get the same results of the small residual tidal signal of 18.6 years using the L-S periodogram. However, we cannot find this signal using the SSA approach, even if adopting the length window of 18.6 years, which is the same as that of Seo et al. [13].

∆J 2 Reconstruction with Climate-Driven Sources
Variations in J 2 , due to changes in surface mass density (∆σ 0 (θ, λ, t)), can be computed by [11]: where k 2 is the Earth's load Love number for degree 2; M denotes the Earth's mass and R is the Earth's mean radius; θ, λ, t and P 2 (cos θ) are the colatitude, longitude, time and the Legendre polynomial, respectively. In this study, we derive contributing sources of ∆J 2 including Antarctic (ANT) and Greenland (GRE) ice sheets, atmosphere (ATM), ocean bottom pressure (OBP) as well as continental glaciers (GLA) by using the geophysical models for the first time, similar to the research of geocenter motion by Zhang and Sun [49]. Readers can refer to Zhang and Sun [49] and the appendix of this study for detailed information of the adopted geophysical models. However, models for terrestrial water storage (TWS) and glacial isostatic adjustment (GIA) still have large uncertainties. Therefore, in a way different from that in Zhang and Sun [49], we take solutions from Sun et al. [39,[50][51][52] that estimates contributing sources of ∆J 2 including TWS and GIA using GRACE data combined with a suitable set of fingerprints (fingerprints approach, FPA) covering the period Jan. 2003 to Jun. 2014. As a novel way to interpret temporal variations of the Earth's gravity field, FPA was proposed by Rietbroek et al. [53] and improved by Sun et al. [39] (refer to Appendix A). Figure 7 shows seven climate-driven ∆J 2 components. ATM, OBP, and TWS significantly contribute to seasonal signals in ∆J 2 , while GIA, ANT, GLA, and GRE mainly influence the linear trend in ∆J 2 . Due to ice melting, ANT, GLA, and GRE are all characterized by a positive linear trend. As an ongoing response of the solid Earth to the past glaciations [15,54], GIA is dominated by a negative trend. We then reconstructed ∆J 2 by summing all climate-driven components and comparing them to the observed ones for the period January 2003 to June 2014 ( Figure 8). The reconstructed ∆J 2 is markedly in line with that of the observed value from SLR in terms of the linear trend and the amplitude and phase of the annual term, as shown in Table 3 and Figure 8. We employ the SSA method, as described in Sections 2.1 and 3.1 to extract the annual and semi-annual signals of ATM, OBP and TWS for the analysis of their amplitudes and phases, while the trends for GIA, ANT, GLA and GRE are calculated. The absolute values of the trends of GIA and GRE are equivalent, and ANT and GLA have the same positive trend. The annual amplitude of TWS is slightly larger than those of ATM and OBP with almost the same phase, while the semi-annual amplitude of TWS is almost twice as large as the other two components. In addition, the variance V of the observed ∆J 2 explained by the reconstructed ∆J 2 can be calculated according to

Comparison and Interpretation of Reconstructed and Observed ∆J 2
where <> denotes variance operator; X represents reconstruction; OBS and MOD X represent the observed and reconstructed ∆J 2 , respectively. In order to quantify the contribution of each component to ∆J 2 , the variance of the observed ∆J 2 explained by each contributing source is calculated when X denotes the component and MOD X represents the individual contributing component.    The L-S periodogram of the climate-driven components is illustrated in Figure 9 for further interpreting the corresponding climate or geophysical sources of the observed J2 variabilities. The ATM, OBP, and TWS are the most significant driving sources for annual variation, while they contribute to the semi-annual variations of ΔJ2 time series in a very The results of the variance V (Table 4) show that the observed ∆J 2 is 81.5% explained by the reconstructed solution, suggesting that all the main characteristics of ∆J 2 have been well explained. In addition, the correlation coefficient between the two time series is 92% (Table 3). Both Antarctica ice sheet and continental glaciers account for less than 15% of the variance of the observed ∆J 2 . Greenland ice sheet has a slightly larger contribution, but still explains less than 25% of the variance. The TWS and atmosphere, on the other hand, are the two largest climate sources causing periodic signals in ∆J 2 , which explain up to 31.9% and 40.1% of the variance of the observed ∆J 2 , respectively. Besides, the results of SSA show that the TWS, ATM and OBP can contribute to the annual and semi-annual variations of the observed ∆J 2 up to 30.8% and 1.0%, 30.1% and 1.6%, as well as 25.4% and 0.7%, respectively. The L-S periodogram of the climate-driven components is illustrated in Figure 9 for further interpreting the corresponding climate or geophysical sources of the observed J 2 variabilities. The ATM, OBP, and TWS are the most significant driving sources for annual variation, while they contribute to the semi-annual variations of ∆J 2 time series in a very small order of magnitude. The peaks at the periods of about 3, 4, 6, and 7 years originate from the mass variations of Greenland and Antarctica ice sheets and continental glaciers, while the OBP makes no contributions to these fluctuations. This fact may imply that the strong ENSO [10] with 4-to 7-year fluctuations leads to the mass variations of Greenland and Antarctica ice sheets and continental glaciers to alter J 2 . However, ENSO has no impact on OBP and the 4-7-year signals in GIA could be artificial, resulting from the L-S periodogram. Besides, it is noteworthy that TWS dominates an evident period of 12.2-year, which could either probably explain the existence of this periodic signal in Figure 5, but a further study may be needed. Due to the short ∆J 2 time series during the GRACE era, those long-period variations cannot be verified. The L-S periodogram of the climate-driven components is illustrated in Figure 9 for further interpreting the corresponding climate or geophysical sources of the observed J2 variabilities. The ATM, OBP, and TWS are the most significant driving sources for annual variation, while they contribute to the semi-annual variations of ΔJ2 time series in a very small order of magnitude. The peaks at the periods of about 3, 4, 6, and 7 years originate from the mass variations of Greenland and Antarctica ice sheets and continental glaciers, while the OBP makes no contributions to these fluctuations. This fact may imply that the strong ENSO [10] with 4-to 7-year fluctuations leads to the mass variations of Greenland and Antarctica ice sheets and continental glaciers to alter J2. However, ENSO has no impact on OBP and the 4-7-year signals in GIA could be artificial, resulting from the L-S periodogram. Besides, it is noteworthy that TWS dominates an evident period of 12.2year, which could either probably explain the existence of this periodic signal in Figure 5, but a further study may be needed. Due to the short ΔJ2 time series during the GRACE era, those long-period variations cannot be verified.

Conclusions
The variations in J 2 reflect the mass redistributions within the Earth system, routinely characterized by a superposition of a secular trend, seasonal variations, interannual variations, decadal variations, and longer-period variations. ∆J 2 of more than 44 years from 1976 to 2020 is analyzed by SSA and L-S for the first time. The results show that the secular linearly decreases a rate of approximately (−5.80 ± 0.08) × 10 −11 /yr and an obvious quadratic rate of (2.38 ± 0.02) × 10 −13 /yr 2 derived from SSA is consistent with the results from Cheng et al. [1] by quadratic polynomial fitting. Besides the seasonal variations, the variations at the~10.6-year period, related to the solar activities, can be notably detected, determining its amplitude of 5.01 × 10 -11 for the first time, whereas there is no 18.6-year periodic signal after removal of the theoretical tidal value from J 2 even if the 18.6-year window size is selected by applying SSA. However, in the spectrum analysis of L-S, the 18.6-year periodic variations and~3 to~6-year fluctuations are all detectable in the de-linear time series. This fact suggests that the advantages of different methods used will lead to different sensitivity to the particular signals that are hard to detect. SSA can maintain the stochastic behavior of the signals it separated to the greatest extent, while small long-period signals can be detected in the spectrum analysis of L-S.
In any case, the ∆J 2 time series during the GRACE era is reconstructed based on the geophysical models and FPA by considering the mass variations from GIA, ANT, ATM, GLA, GRE, OBP, and TWS. The reconstructed ∆J 2 , by summing all the climatedriven components, is markedly close to that of the observed ∆J 2 with a high correlation coefficient of 92%. The reconstructed ∆J 2 can also explain up to 81.5% of the variance of the observed value. The most significant climate-driven sources for seasonal signals in the ∆J 2 time series mainly originated from ATM, TWS, and OBP, which explains up to 40.1%, 31.9%, and 26.3% of the variance, respectively. These three components, in total, can account for the annual and semi-annual variations in the observed ∆J 2 , up to 86.3% and 3.3%, respectively. However, GRE, ANT, and GLA dominate the positive linear trend with interannual periodic fluctuations of~3 to~7-years and GIA contributes to a large negative trend.
Author Contributions: H.Y. designed the study, analyzed data, prepared the results and figures, wrote and revised the original manuscript. Q.C. co-designed and led the study. Y.S. computed the climate-driven components of ∆J 2 in Section 4.1. K.S. discussed and given suggestions. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: The C20 time series from 1976 to 2020 can be downloaded from the CSR website (http://download.csr.utexas.edu/pub/slr/degree_2/Long_term/C20_Long_Term.txt, accessed on 20 January 2021). ered and predicted by geophysical models during the GRACE era from August 2002 to June 2014 [49]. The input-output method [55], altimetry measurements [56] and GRACE gravimetry measurements [19] are the popular methodologies to estimate mass variations over continental glaciers and ice sheets. Here, the input-output method is utilized to estimate the mass changes in both Greenland and Antarctica ice sheets [57]. The Regional Atmospheric Climate Model (RACMO) 2.3 surface mass balance (SMB) [58] is used as an input, and the ice discharge acceleration with 6.6 Gt year −1 for Greenland [55] and 9.0 Gt year −1 for Antarctica [59] is the output. The contribution of continental glacier mass balance (GLA) to ∆J 2 is estimated according to Marzeion et al. [60]. In addition, ocean bottom pressure (OBP), inducing the ∆J 2 , is derived from the GAD and the continental atmosphere (ATM) is based on the result of GAC deducting GAD product. GAC and GAD are based on the European Center for Medium-Range Weather Forecasts (ECMWF) [61] and the Ocean Model for Circulation and Tides (OMCT) [62].
TWS and GIA are computed by using the fingerprints approach (FPA) with the GRACE-Level2 RL06 data from the Center for Space Research (CSR). Assuming that a unit mass change occurs at a certain point or a certain area, the global mass variation caused by it can be expressed as f M [63]. When the mass change α f M occurring at a point or local area is not the unit mass, the global mass variation can be expressed as α f M f M . According to the conversion between the mass coefficient and the spherical harmonic coefficient, the global mass change G caused by the sum of the mass variation of each point can be estimated via where F, a matrix of spherical harmonic analysis with m columns, is the set of fingerprints f M . In this study, the C 10 , C 11 and S 11 are co-estimated together with the C 20 , thus G is a matrix with n-4 rows. T is the truncation matrix to eliminate C 10 , C 11 , S 11 , and C 20 in F, which is expressed as follows: where α can be obtained by using the least-squares method. Therefore, the contributing sources to ∆C 20 can be computed by