As atmosphere, ocean, ice, and terrestrial water are redistributed, the center of mass (CM) of the Earth's system moves and the accompanying loading yields global surface deformation. In Australia, when GPS surface displacements were corrected for local mass change (hydrology, atmosphere, and ocean) with Gravity Recovery And Climate Experiment (GRACE) data, the residual GPS data reveal a peculiar seasonal mode of continental deformation. During the southern summer, the entire continent coherently shifts northwest by ~1 mm and the southeastern part is uplifted, while the northwestern part is subsided by 2–3 mm and the opposite patterns of deformation are observed during the southern winter. Such characteristic deformation could be understood to be a result of the Earth's elastic response to globally averaged surface mass load, generally heavier in Europe during southern summer and in the South Pacific Ocean during southern winter. It was found that such deformation is even larger than local hydrology‐induced loading effects in horizontal motion over Australia. A simple method of determining locations of the CM was developed by combining GPS and GRACE data; the latter being insensitive to the CM motion but sufficiently accurate to remove the local hydrologic and atmospheric effects in GPS data. The CM signals are pronounced over systematic errors in GPS and GRACE data. The CM coordinates estimated by inversion of the Australian GPS data set and GRACE agree with the geocenter motions determined by satellite tracking analysis. This study suggests an independent way of monitoring the CM motion entirely based on two distinct geodetic measurements of GPS and GRACE.

1 Introduction Mass redistribution within the Earth's system (solid body + fluid envelope including atmosphere) perturbs the Earth's gravity field through direct effect of mass movement and indirect effect of accompanying load deformation [e.g., Farrell, 1972; Chao et al., 1987]. The surface deformations are measured by precision positioning system such as the Global Positioning System (GPS) within a few mm accuracy [e.g., Blewitt et al., 2016] and the global gravity changes are continuously measured by Gravity Recovery And Climate Experiment (GRACE) within a couple of 10−8 m/s2 at a spatial scale of ~400 km [e.g., Han et al., 2013; Watkins et al., 2015]. The GRACE gravity fields are observed with respect to a coordinate system centered at an instantaneous center of mass (CM) of the whole Earth including atmosphere and fluid components. This is implicitly realized by estimating the gravity fields from tracking GRACE satellites orbiting around CM of the Earth. On the contrary, the GPS 3‐D positions are referenced to a coordinate system such as International Terrestrial Reference Frame (ITRF) with the origin that is not necessarily coincident with the “instantaneous” CM. The geocenter motion, a motion of CM with respect to a reference frame realized by establishing the global satellite tracking network, has been best determined by satellite laser ranging (SLR) techniques [Pavlis, 2003]. The origin of such reference frame is (ideally) the center of figure (CF) of the Earth but is practically a center of the network [Wu et al., 2012]. The geocenter motion is simultaneously estimated with the SLR orbit and gravity field solutions [Cheng et al., 2010]. Blewitt et al. [2001] detected the seasonal cycle of deformation induced by CM migrating on the Earth's surface from global GPS data. However, determination of the interannual CM motion and its load deformation was difficult since the interannual GPS signals are largely related to tectonic motion (not surface loading). Wu et al. [2002] found that the GPS analysis of CM could be contaminated by sparse distribution of the GPS stations and aliasing of unknown (uncorrected) “high‐degree” deformation. Davis et al. [2004] found a strong correlation in radial (vertical) displacement between GPS measurements and GRACE‐inferred deformation over the Amazon basin with inclusion of the CM‐induced radial load displacement in GRACE data. Kusche and Schrama [2005] and Wu et al. [2006] determined global surface mass distribution by combining simultaneously GPS, GRACE, and ocean models. More recent analyses reported substantial correlation between GPS and GRACE at different locations due to different kinds of surface mass loads [e.g., Nahmani et al., 2012; Fu et al., 2013, and references therein]. Most of the studies exploit the CM estimates from various techniques including SLR determination [Cheng et al., 2013], GRACE and ocean model estimates [Swenson et al., 2008], geophysical model estimates [Chen et al., 1999], and/or joint inversion of multiple data sets [Wu et al., 2012]. However, there still exists important discrepancies among different estimates from various techniques and geophysical models [Ries, 2013]. In this study, I report the regionally coherent mode of seasonal deformation prevailing over Australia quantified from GPS and GRACE data; the Australian continent undergoes periodic deformation of up to a couple of millimeters in horizontal component and ~3 mm in vertical component seasonally. It is found that the horizontal motion is even larger than local hydrologic and atmospheric load deformation. Such motions are analyzed considering the CM migration on the Earth's surface. This research demonstrates that simple differences between GPS and GRACE can provide meaningful information on load deformation caused by the CM motion. Furthermore, I develop a method to estimate the coordinates of CM motion using a regional set of GPS data combined with GRACE‐inferred displacement. The results are compared with the geocenter motions determined from SLR.

2 Center of Mass of the Earth's System m and its location in each coordinate axis (i.e., linear moment) over the Earth, for example, for the X axis. It can be shown that the CM coordinates are directly related to the “degree 1” gravitational potential field as follows [e.g., Jekeli, 2015 (1a) (1b) (1c) , , and are the degree 1 gravitational potential coefficients. The coefficients are dimensionless and normalized following the convention in Heiskanen and Moritz [ 1967 a, is a radius of the mean Earth's sphere used to expand the gravitational potential field in terms of a series of spherical harmonic functions. The coordinates of CM in a Cartesian coordinates system are defined by a volume integral of product of an infinitesimal mass element dand its location in each coordinate axis (i.e., linear moment) over the Earth, for example,for theaxis. It can be shown that the CM coordinates are directly related to the “degree 1” gravitational potential field as follows [e.g.,]:where, andare the degree 1 gravitational potential coefficients. The coefficients are dimensionless and normalized following the convention in]. The constant,, is a radius of the mean Earth's sphere used to expand the gravitational potential field in terms of a series of spherical harmonic functions. t is written as (2) G is the gravitational constant and M is total mass of the Earth's system. is the fully normalized associated Legendre function of degree 1 and order m at colatitude θ. The degree 1 spherical harmonic functions are formed by multiplying sine and cosine functions of longitude, λ to the associated Legendre function [Heiskanen and Moritz, 1967 , , and ) can be time dependent as mass redistributes within the Earth's system. Most of the time variable mass change occurs on the Earth's surface associated with water, ice, and atmospheric mass redistribution besides solid Earth mass change due to earthquakes, viscous postglacial rebound, and mantle flow. The degree 1 component of the Earth's gravitational potential field at time,is written aswhereis the gravitational constant andis total mass of the Earth's system.is the fully normalized associated Legendre function of degree 1 and orderat colatitude. The degree 1 spherical harmonic functions are formed by multiplying sine and cosine functions of longitude,to the associated Legendre function []. The degree 1 potential coefficients (, and) can be time dependent as mass redistributes within the Earth's system. Most of the time variable mass change occurs on the Earth's surface associated with water, ice, and atmospheric mass redistribution besides solid Earth mass change due to earthquakes, viscous postglacial rebound, and mantle flow. (3) , , and , produce the degree 1 gravitational potential change through [e.g., Kusche and Schrama, 2005 ρ e = 5517 kg/m3) and a density of water (ρ w = 1000 kg/m3). The degree 1 load Love number for gravitational potential is given in the coordinate system that identifies the CM change; e.g., = 0.021 in the CF system [Blewitt, 2003 Assuming mass change occurs entirely as a result of surface mass redistribution and the accompanying loading, the degree 1 component of the surface mass change, which is responsible for gravitational potential change in equation 2 as well as CM change in equation (1), can be expressed in terms of equivalent water height as follows:where the dimensionless degree 1 coefficients of equivalent water height,, and, produce the degree 1 gravitational potential change through[e.g.,] with a mean density of the Earth (= 5517 kg/m) and a density of water (= 1000 kg/m). The degree 1 load Love number for gravitational potentialis given in the coordinate system that identifies the CM change; e.g.,= 0.021 in the CF system []. (4) (5) Expressing the gravitational potential and equivalent water height changes in terms of the coordinates of CM and explicitly representing the associated Legendre functions, equations 2 and 3 are simplified as The only difference between the degree 1 variations of gravitational potential and equivalent water height is a scale factor in front of the brackets in equations 4 and 5. That is, their spatial patterns are identical. It is not necessarily true for higher‐degree (l ≥ 2) components. Farrell [ 1972 u 1 (up), e 1 (east), and n 1 (north) can be expressed as [e.g., Kusche and Schrama, 2005 (6) (7) (8) and are the load Love numbers for vertical and horizontal load displacements; e.g., = –0.269 and = 0.134 in the CF system [Blewitt, 2003 h 1 loading, or alternatively, due to the CM departing from the center of the coordinate system on which the deformation is defined. The surface mass (terrestrial water storage and atmospheric and oceanic mass) applies loads on the lithosphere and yields surface and interior deformation.] developed a theory on the linear response of elastic deformation to surface mass load in a spectral domain and derived also the response (Green's) functions in a spatial domain. The degree 1 load excites only the degree 1 (elastic) displacement and gravitational potential. The degree 1 surface displacements of(up),(east), and(north) can be expressed as [e.g.,whereandare the load Love numbers for vertical and horizontal load displacements; e.g.,= –0.269 and= 0.134 in the CF system []. The degree 1 surface deformation of equations 6-8 are due to the degree 1 surface mass,loading, or alternatively, due to the CM departing from the center of the coordinate system on which the deformation is defined. It is useful to emphasize that spatial patterns of the degree 1 load, h 1 , and its responses in gravitational potential field, V 1 , and in vertical displacement, u 1 , are identical. However, the responses in horizontal displacements, e 1 and n 1 , are distinct each other and from u 1 , and thus, three coordinates of degree 1 displacements can provide full constraints on the CM coordinates through equations 6-8. In principle, 3‐D surface deformation from a single station can provide the complete information on CM motion, if higher‐degree (l ≥ 2) components of the deformation are known and properly removed (such as from GRACE).

3 SLR Observations of the Center of Mass I examined the CM coordinate estimates (or geocenter motions) determined from the analysis of multiple SLR data. The estimates consist of monthly time series of “instantaneous” CM coordinates with respect to the SLR stations' coordinate system that is consistent with ITRF2005 [Cheng et al., 2013]. Here the CM coordinate changes are attributed to mass change on the Earth's surface and accompanying load deformation. Using equation 5, the monthly mean patterns of the degree 1 surface water load were computed from the SLR's CM estimates during 2003 to 2015 (Figure 1). Each panel shows a mean spatial pattern of water load, h 1 , and of horizontal load displacement vector (e 1 , n 1 ) every month. The vertical displacement, u 1 , should be identical to h 1 except for a scale difference of . That is, the degree 1 load corresponding to 40 mm of equivalent water mass yields the degree 1 subsidence of 2 mm. Figure 1 Open in figure viewer PowerPoint Monthly maps of globally averaged (i.e., degree 1) surface mass redistribution quantified in terms of equivalent water height. The SLR‐determined geocenter motions were interpreted as a result of mass redistribution on the surface and the accompanying load deformation. Monthly snapshots were calculated by stacking multiyears of SLR geocenter estimates after removing linear trends in the estimates. The arrows indicate the horizontal displacements induced by the Earth's elastic response to the surface mass load (the longest arrow corresponds to the largest movement of ~1 mm). The vertical displacement is identical to these maps but with a scale factor of −0.05 (see the text). In January–February (Figure 1), for example, there exists a net surplus of surface mass distribution around Europe and north Africa, while a mass deficit is found in the South Pacific Ocean equivalent to 20–30 mm of water. It results in subsidence of 1–1.5 mm over Europe and North Africa and an uplift in the South Pacific Ocean due to the Earth's elastic response to the surface mass load. There are outward lateral motions of <1 mm away from the South Pacific Ocean and inward motions toward Europe. The Australian continent happens to be located where horizontal gradients of the degree 1 load, h 1 , is largest between the positive and negative peaks and thus undergoes the largest degree 1 horizontal displacements in the world. This calculation in Figure 1 based on the SLR's CM estimates predicts that the Australian continent seasonally oscillates by moving toward Europe during (southern) summer, toward North America during autumn, toward the South Pacific Ocean during winter, and toward central Asia during spring, as a result of seasonal migration of global mean (i.e., degree 1) surface mass.

4 GPS and GRACE Observations of the Degree 1 Load Displacement I analyzed GPS displacement data in Australia to verify the seasonally oscillating motions predicted from the SLR's CM estimates. The 3‐D daily position data (north, east, and up) of the Australian GPS network, processed by Nevada Geodetic Laboratory, University of Nevada Reno (UNR), were used in this study. These GPS solutions are referenced to a coordinate system, ITRF08 [Altamimi et al., 2011]. The origin is coincident to the long‐term mean CM realized by averaging SLR solutions, which is different from CF only by a constant offset [Wu et al., 2012]. The GPS positions were corrected for tidal deformations (body, ocean, and pole) and thus imply instantaneous deformation caused by atmosphere, nontidal ocean, and terrestrial water loads as well as tectonics including earthquakes [Tregoning et al., 2013]. The total of 14 Australian stations with longer than 9 years of continuous position solutions to year 2015 was considered. The GPS time series data include biases in all three components introduced by sporadic equipment changes and earthquakes. The UNR solutions provide the information of the timings of such events for all processed stations. I determined and removed those biases empirically by fitting a function of seasonal sinusoids, linear trend, and offset parameters at the reported offset epochs. In addition, monthly GRACE gravity field data of Release‐05 (RL05) Level‐2 (L2) products processed by Center for Space Research (CSR), University of Texas, were used [Bettadpur et al., 2014]. The atmosphere and ocean models that were removed from the CSR L2 products were restored in the GRACE data to make GRACE gravity change consistent with GPS deformation [Dobslaw et al., 2013]. The L2 spherical harmonic coefficients were used up to degree and order 40. (The effects of different truncations are discussed in section 6.) Farrell [ 1972 Davis et al. [ 2004 u 1 is defined as (9) u GPS is the GPS displacement and u GRACE is the computed displacement using the GRACE L2 coefficients (degree l ≥ 2). Because GRACE gravity changes are referenced to a coordinate system centered at the instantaneous CM of the Earth's system, the degree 1 components of GRACE gravity changes in such coordinate system are, by definition, trivial (zeroes). I used the GRACE gravity data to compute 3‐D load deformation using the elastic load Love numbers after] and. [] in the GRACE context. The difference between GPS‐measured and GRACE‐inferred displacements should indicate the origin translation to the instantaneous CM from the origin of the GPS coordinate system. The vectoris defined aswhereis the GPS displacement andis the computed displacement using the GRACE L2 coefficients (degree≥ 2). The vector, u 1 , is the degree 1 displacement presumably related to the CM coordinates through equations 6-8. This is true only if gravity changes and surface displacements are strictly due to an elastic deformation by surface mass load. For example, gravity change and displacement by tectonics and earthquakes violate this assumption. The plate motions are generally the most prominent signals measured in the GPS horizontal displacements. GPS observations have shown that the Australian continent steadily moves northeast by ~7 cm/yr quantified in a time scale of a decade or two [Tregoning, 2003]. The large earthquakes around the world produce measurable coseismic and postseismic deformation in Australia by interrupting GPS time series episodically and gradually [Tregoning et al., 2013]. Surface mass redistribution (water, snow, ice, and atmosphere) is characterized dominantly at the seasonal and interannual time scales. In practice, removing linear trends in GPS and GRACE data can effectively filter out nondesired tectonic/earthquake signals as well as secular changes in surface mass redistribution [Blewitt et al., 2001]. Figure 2 presents the daily observations of east, north, and up components of the degree 1 displacement, u 1 , from GPS‐GRACE, after detrending the time series from all 14 GPS stations considered in this study (GPS site locations are shown in Figure 3). The monthly GRACE data were interpolated into daily intervals of the GPS data. The detrended degree 1 displacements evaluated at those GPS stations using the SLR's CM solutions are compared. The parameters of annual and semiannual sinusoids were used to estimate the seasonal change from each time series. The seasonal displacements of ~1 mm (±0.03 mm) in horizontal components and of ~2 mm (±0.07 mm) in vertical component are evident in the GPS‐GRACE time series and are consistent with the SLR predicted motions, except larger amplitudes estimated from GPS‐GRACE (70, 20, and 40% larger in east, north, and vertical component, respectively). The doubling of the amplitude in vertical CM displacement compared with the horizontal one is consistent with the vertical load Love number, , being double the horizontal number, . Wahr et al. [2013] also found the asymptotic ratio of ~2 between the vertical and horizontal load deformation when the displacement is computed far away from the source mass load. Figure 2 Open in figure viewer PowerPoint l = 1) surface mass were determined using the GPS displacements from 14 stations in Australia after removing the GRACE‐inferred local (l ≥ 2) load displacements (black dots). The detrended monthly degree 1 load displacements evaluated at GPS stations using the SLR geocenter motions are shown in green solid lines. The seasonal components (annual and semiannual sinusoids) were determined from each time series of GPS‐GRACE and SLR, respectively, shown in red and magenta solid lines. The geographic locations of 14 GPS stations are shown in Figure Daily load displacements caused by the degree 1 (= 1) surface mass were determined using the GPS displacements from 14 stations in Australia after removing the GRACE‐inferred local (≥ 2) load displacements (black dots). The detrended monthly degree 1 load displacements evaluated at GPS stations using the SLR geocenter motions are shown in green solid lines. The seasonal components (annual and semiannual sinusoids) were determined from each time series of GPS‐GRACE and SLR, respectively, shown in red and magenta solid lines. The geographic locations of 14 GPS stations are shown in Figure 3 Figure 3 Open in figure viewer PowerPoint Monthly patterns of 3‐D displacements induced by the degree 1 surface mass loading were determined by the difference between GPS and GRACE data. The time series of GPS and GRACE displacement data from 2003 to 2015 were detrended and stacked to estimate the monthly patterns of deformation at each station. The vertical and horizontal displacements are shown with colored dots and the solid black arrows, respectively. I also determined the monthly spatial patterns of the degree 1 displacements, u 1 , from GPS‐GRACE as shown in Figure 3 with arrows indicating horizontal motions and colored dots for vertical motions (positive for uplift). More than 9 years of GPS and GRACE data were stacked and averaged into a monthly interval to estimate monthly snapshots of the 3‐D displacements. The continent‐wide northwest shift in summer (December–February) and southeast shift in winter (June–August) are observed. It is also observed that the continent is subject to seasonal vertical motion with the southeastern Australia elevated and the north of Western Australia lowered by a few millimeters during summer and the opposite during winter, creating the vertical slope (tilt) across the continent by <5 mm over ~3500 km. These horizontal and vertical motions are indicative of the spontaneous (elastic) responses to surface mass loads imposed around Europe in summer and in the South Pacific Ocean in winter following the seasonal migration of the degree 1 load (i.e., CM). The 3‐D CM deformation found from GPS‐GRACE data (Figure 3) is in general consistent with the SLR results (Figure 1). Two independent data sets suggest that the entire Australian continent tilts northwest in summer and southeast in winter and gyrates as a whole in a clockwise direction, making an elliptical motion seasonally by ~1 mm in the northwest‐southeast axis and by a smaller amount in the perpendicular axis, chasing the center of the Earth system's mass. Meanwhile, the displacements by the local (l ≥ 2) hydrology, ocean, and atmosphere loading were computed using the GRACE data and compared with the degree 1 deformation (GPS‐GRACE) in Figure 4. Their monthly patterns are shown in Figure 5. It was found that the CM deformation is larger than the local load deformation in the horizontal component, while it is 2–3 times smaller in the vertical displacement. The CM displacement (not hydrology loading) is the most dominant signal characterizing the seasonal horizontal motion in Australia. While the CM displacements were spatially coherent among different stations, the GRACE local loading data present diverse timings and amplitudes of seasonal variations reflecting spatially variable atmospheric, oceanic and hydrological load cycles. Figure 4 Open in figure viewer PowerPoint l ≥ 2) hydrologic, atmospheric, and oceanic loads, evaluated at 14 GPS stations using GRACE (black). The spherical harmonic degree and order up to 40 were used for GRACE data. They are compared with the estimated seasonal motion of the degree 1 displacements from GPS minus GRACE, averaged from all 14 stations in Australia (red, the same as Figure The time series of load displacements due to local (≥ 2) hydrologic, atmospheric, and oceanic loads, evaluated at 14 GPS stations using GRACE (black). The spherical harmonic degree and order up to 40 were used for GRACE data. They are compared with the estimated seasonal motion of the degree 1 displacements from GPS minus GRACE, averaged from all 14 stations in Australia (red, the same as Figure 2 ). Figure 5 Open in figure viewer PowerPoint Monthly patterns of 3‐D displacements induced by the local hydrologic, and atmospheric, and oceanic loads estimated from GRACE data. The time series of GRACE displacement data from 2003 to 2015 were detrended and stacked to estimate the monthly patterns of deformation at each station. The vertical and horizontal displacements are shown with colored dots and the solid black arrows, respectively.

5 Center of Mass Estimated From GPS and GRACE u 1 = R(θ, λ)X CM where a 3 × 3 matrix, R(θ, λ) is determined by the latitudes and longitudes of GPS stations and the degree 1 load Love numbers, ( , , and ) such as (10) The CM motions were estimated using the degree 1 displacement data of GPS‐GRACE. Equations 6-8 show three independent linear relationships between the CM coordinates and the degree 1 load displacements. They can be rewritten aswhere a 3 × 3 matrix,) is determined by the latitudes and longitudes of GPS stations and the degree 1 load Love numbers, (, and) such as Blewitt and Clarke, 2003 R(θ, λ), as follows: (11) This matrix can be interpreted as a rotation matrix from the geocentric to topocentric coordinate system [also, equation (13)] and is successively scaled by the Love numbers. Inversely, the CM coordinates can be computed by converting the degree 1 displacement vectors epoch by epoch and station by station through inversion of the station‐dependent factor,), as follows: The procedure taken in this paper is as follows: At each GPS station, the 3‐D coordinates of GPS displacements are differenced with the GRACE‐inferred load displacements to estimate the degree 1 load displacements following equation 9. Monthly GRACE data are interpolated to daily interval to match the GPS sampling. Both GPS and GRACE time series are detrended to eliminate tectonic signals particularly in the horizontal components. Then, equation 11 is used to convert the degree 1 load displacements from each station to obtain the CM coordinates. The time series of X CM from multiple GPS stations are combined to estimate the coordinates of the CM motions. Figure 6 shows the daily time series of three components of the CM motion, X CM , converted through equation 11 from GPS minus GRACE data, u 1 . Unlike the degree 1 displacements at each GPS station, the CM motion coordinates are no longer station dependent, allowing simple average of all data from different stations to reduce (mostly) GPS data noise in the CM estimates. The daily estimates present a few centimeters of variations rather uniformly in all three coordinates of CM, in contrast, to those of the surface CM deformation (Figure 2). The monthly averages of the daily estimates were computed and found to be generally within ±10 mm. Figure 6 Open in figure viewer PowerPoint Daily time series of the CM coordinates, X CM , Y CM , and Z CM (i.e., geocenter motion) were estimated by inverting the differences of GPS and GRACE data (black dots). The monthly averages of the daily estimates were also calculated (red dots). I compared total 12 years of the monthly CM motion estimates from GPS‐GRACE with the SLR's geocenter solutions after removing the linear trend in the SLR solutions in Figure 7a. The timings of peak and trough occurrences agree well each other in all three components. The seasonality of the Y CM estimates is better characterized from the GPS/GRACE solutions than the SLR's Y CM solutions (relevant to the geopotential component of in the SLR orbit and gravity analysis) for the period of 2004–2009. However, the GPS/GRACE solutions become worse in 2010 and beyond. The average monthly changes were also compared by stacking all 12 years of the monthly time series (Figure 7b). Two independent solutions are consistent; the peak is observed in January, December, and March, respectively, for X CM , Y CM , and Z CM ; the trough is found in July, May, and August, respectively, for X CM , Y CM , and Z CM ; both data indicate larger variability in X CM and Z CM components than in Y CM . In general, the GPS/GRACE CM solutions present a larger variability than the SLR solutions. Figure 7 Open in figure viewer PowerPoint (a) Comparison of the monthly estimates of the CM coordinates in 2003–2015 from GPS‐GRACE (red) and SLR geocenter solutions (blue). (b) The average monthly variations of the CM coordinates during 2003–2015 from GPS‐GRACE and SLR. The error bar indicates the standard deviation of the variations from the average in each month over 2003–2015. In GPS analysis, nongeophysical signals are often observed at the frequencies of one cycle per the GPS draconitic year (~351.4 days) and of its integer multiples, originating from GPS orbital errors and time modulation of site‐dependent multipath [Ray et al., 2008; Tregoning and Watson, 2009]. The magnitude can reach up to a few millimeters [King and Watson, 2010]. I examined the frequency contents of two CM coordinate time series in a spectral domain. Figure 8 shows the amplitude (square root of power) spectral density of each time series. Both GPS/GRACE and SLR CM solutions present the primary peaks centered at the period of a solar year (~365.25 days) distinguished from the GPS draconitic year. The CM estimates from GPS/GRACE are less contaminated by GPS systematic errors around the annual band. However, GPS/GRACE CM solutions exhibit the secondary peaks at half draconitic year (~175.6 days) as well as at half solar year (for Y CM , only at half draconitic year). The SLR solutions present the peaks at half solar year for X CM and Z CM and none at half draconitic year. Around the semiannual band, the GPS/GRACE CM solutions are significantly corrupted by GPS specific errors. Figure 8 Open in figure viewer PowerPoint The amplitude (square root of power) spectral density of monthly time series of the CM motion estimates from this study and the SLR analysis. The magenta vertical bars indicate the periods of the solar year (365.25 days) and its half (182.63 days). The green vertical bars indicate those of GPS draconitic year (351.5 days) and its half (175.7 days). The CM estimates can be improved by including additional GPS‐GRACE data from other places of the world. The GPS systematic errors may cancel, while the CM motions will become more pronounced by stacking global data through equation 11. Different places exhibit different sensitivities to the CM coordinates. The sensitivity is determined by the matrix R− 1(θ, λ) of equation 11, which is composed of nine elements like ∂X CM /∂u 1 , ∂Y CM /∂u 1 , ∂Z CM /∂u 1 , ∂X CM /∂e 1 , ∂Y CM /∂e 1 , ∂Z CM /∂e 1 , ∂X CM /∂n 1 , ∂Y CM /∂n 1 , and ∂Z CM /∂n 1 (i.e., partial derivatives of three CM coordinates with respect to three degree 1 displacements). The analytic expression of each element can be easily found from equation 10. The results are presented in Figure 9. The CM sensitivity is roughly twice larger with respect to the horizontal u 1 displacement. Figure 9 Open in figure viewer PowerPoint R− 1(θ, λ) of equation The sensitivities (partial derivatives) of the CM coordinates with respect to the degree 1 displacements. All nine components were computed from the matrix) of equation 10 . Note two different scale bars for the sensitivity to the vertical motion and to the horizontal motion. The X CM is most well constrained by the east component of u 1 from stations around longitudes close to 90°E or 90°W (e.g., America and Asia), by the north component from stations in higher latitudes and longitudes close to 0°E or 180°E (e.g., Alaska, the Russian Far East, Iceland, and northern Europe). The vertical component of u 1 from stations in lower latitudes and longitudes close to 0°E or 180°E will also contribute to determine X CM . For Y CM , the east component from stations in Pacific, New Zealand, Africa, UK (longitudes close to 0°E or 180°E), the north component from Canada, Siberia (higher latitudes and longitudes close to 0°E or 180°E), and the vertical component from India and Central America (lower latitudes and longitudes close to 90°E or 90°W) are most sensitive. For Z CM , the north component from the equatorial stations is most sensitive, and the vertical component from the polar regions is sensitive as well; however, the east component is not sensitive at all. For an example of the u 1 displacement data from the east coast Australian stations, the vertical, east, and north components are most sensitive to X CM , Y CM , and Z CM , respectively.

6 Errors in GRACE and GPS Data Processing Systematic errors associated with GRACE and GPS data processing were examined for their impact on the CM estimation. Although GRACE data are provided degree and order up to 90 (e.g., CSR L2 data) [Bettadpur et al., 2014], the coefficients are valid globally up to ~40 and around the polar regions up to 60 or higher. The spherical harmonic truncation will have an effect on computation of local displacement u GRACE used to compute u 1 . Figure 10a shows the square root of power spectrum per degree (i.e., degree amplitude spectrum) and that of cumulative power per degree for vertical load displacement computed from global models of hydrology [Rodell et al., 2004] and atmosphere and ocean mass variation [Dobslaw et al., 2013]. The global average of vertical load displacement is computed to be several millimeters. Figure 10b shows the daily time series of the synthetic displacements at Alice Spring computed from different degree bands. The contribution from degrees higher than 40 is anticipated to be less than 0.1 mm. Figure 10 Open in figure viewer PowerPoint (a) The degree amplitude spectrum and the cumulative spectrum of the vertical load displacement computed from hydrology, atmosphere, and ocean models. (b) Daily time series of the synthetic displacements of atmosphere and ocean mass loads, computed at different degree bands, evaluated at the central Australia (Alice Spring). The other major systematic error in GRACE data is a large uncertainty in the second degree zonal harmonic coefficient (J2) estimate. This component of the geopotential field is highly correlated with the GRACE orbital state parameters and poorly constrained by GRACE intersatellite ranging data [Han et al., 2008]. I used the GRGS RL03 solutions of geopotential fields determined using Laser Geodynamics Satellite (LAGEOS) and GRACE tracking data [Lemoine et al., 2014]. The J2 component was improved by including LAGEOS laser ranging data. This is the set of the solutions most distinct from other project solutions including CSR, JPL, and GFZ in terms of GRACE data processing strategy. The GRGS solutions were constrained within the least squares inversion and available degree and order up to 80. The difference between GRGS solution (out to degree 80) and CSR solution (truncated to 40) may be a good representation of the upper bound of GRACE processing error including J2 uncertainty as well as spherical harmonic series truncation. Figure 11 shows the difference between two GRACE solutions projected to the CM coordinates. The effect of the J2 difference and truncation in GRACE data are insignificant compared to the CM motion signals (Figure 7). Figure 11 Open in figure viewer PowerPoint Similar to Figure 7 but they represent the effects of GPS and GRACE data processing errors (difference between two GPS solutions and between two GRACE solutions) on the CM coordinates. Tregoning and Watson [2009] and Steigenberger et al. [2009] reported that different tropospheric delay modeling may yield systematic difference in GPS position solutions as large as ~1 mm at annual and semiannual periods and found that position solutions can be improved by using Vienna Mapping Function 1 (VMF1) against Global Mapping Function (GMF) [Boehm et al., 2006]. To quantify the effect of GPS tropospheric modeling on the CM estimates, I analyzed the difference between two sets of GPS daily solutions from UNR and JPL; the former uses GMF, while the latter uses VMF1. Figure 11 also presents the difference between two GPS solutions shown as the difference in the CM coordinate estimates. The processing differences of GPS and GRACE data are as large as root‐mean‐square (RMS) of 1.5–2.5 mm in the CM coordinates, while the CM motion signals are significantly larger with RMS of 7–8 mm.

7 Discussion I found that the most dominant signal in horizontal deformation of the Australian continent, beside the (steady) plate tectonic motion, is the elastic load displacement induced by global mean surface mass redistribution (or CM migration). The regular trajectory of seasonal mass migration (a composite of terrestrial water, atmosphere, and ocean loads being heavier over Europe in February and over the South Pacific Ocean in August) excites a peculiar seasonal mode of surface displacements; (1) the Australian continent undergoes seasonal clockwise gyration with an elliptical shape elongated in a northwest‐southeast direction with the major axis radius of ~1 mm and the orthogonal minor axis radius of <0.3 mm. It moves northwest during southern summer and southeast six months later; (2) the continent is tilted northwest in summer and southeast in winter, creating the vertical slope across the continent by <5 mm over ~3500 km. In Australia, the CM deformation is even larger than the local hydrological, oceanic, atmospheric load deformation in horizontal displacement, while it is 2–3 times smaller in vertical displacement. The Australian continent happens to be located between the peaks and troughs of the global mean (degree 1) surface mass or, in other words, where the horizontal gradient of the degree 1 surface mass is largest. Therefore, it is subject to the largest horizontal motion caused by the degree 1 loading among other countries. A simple method was demonstrated to estimate the CM motion based on a combination of GPS and GRACE data by exploiting the fact that GRACE is insensitive to the CM motion but sensitive to other local (high‐degree components) loading, while GPS is sensitive to all components of loading. The same idea has been exploited by Davis et al. [2004] and Swenson et al. [2008]. This study showed that even the regional set of GPS data (3‐D positions) can provide full constraints on the CM motion. The systematic errors caused by different data processing strategies associated with spherical harmonic truncation and the J2 uncertainty in the GRACE data were found relatively insignificant, compared to the size of CM signals, as long as the GRACE data up to degree and order 40 are taken. The GRACE data are accurate enough to remove local (other than degree 1) loading effects from GPS in Australia. The systematic errors related to GPS tropospheric delay modeling are of less concern for the CM determination. However, I found that the GPS/GRACE CM estimates present most likely nongeophysical signal at the frequency of twice per the GPS draconitic year. No signal at such frequency is observed from SLR geocenter solutions. Both GPS/GRACE and SLR solutions reveal the CM solutions largest at the annual cycle (365.25 days) that is distinguished from the GPS draconitic period (351.4 days). The CM motion is the most dominant signal in GPS‐GRACE observations of the degree 1 displacement. Stacking multiple GPS data sets corrected with GRACE data into the CM coordinates effectively reduces data noise mostly in the GPS data. The method can be extended to process GPS data from other regions as well as globally, which will improve the resolution of the CM motion against GPS systematic errors. In this study, GPS displacement and GRACE gravity change are attributed to the result of surface mass loading on the elastic Earth. The tectonic and earthquakes signals existing in GPS data make it difficult to determine particularly the interannual change in the CM variations. The seasonal displacements induced by surface temperature change are computed to be 1–2 mm at a regional spatial scale [Fang et al., 2014] and will bias the CM motion determination. The Earth's variable structures and rheological properties, departing from a simple laterally homogeneous elastic Earth's model used in this study, will influence surface deformations as well as geopotential fields [e.g., Dill et al., 2015]. More complex and comprehensive mechanical models of the realistic Earth's response should be developed to understand geodetic data for global and regional load deformation [e.g., Sauber et al., 2014]. GPS deformation measurements are being increasingly exploited to study terrestrial water and ice mass storage through load deformation calculation [e.g., Fu et al., 2013; Argus et al., 2014; Borsa et al., 2014]. When the surface loading is known such as tides, the GPS measurements can also be used to infer the Earth's structure [e.g., Martens et al., 2016]. These geodetic approaches should be meaningful only with proper understanding of the CM‐induced deformation in GPS data that may amount to be even larger than the local hydrology and ocean loading. Ries [2013] and Melachroinos et al. [2013] emphasized the importance of having models of annual geocenter motions for accurate sea level determination to improve satellite orbits. This study demonstrates a kinematic way of determining seasonal CM motions independently from the dynamic method based on gravity and orbit analysis of SLR satellites.

Acknowledgments This work was funded by University of Newcastle to support NASA's GRACE and GRACE Follow‐On missions. I thank DLR for GRACE telemetry data and JPL, CSR, and GFZ for high‐quality Level‐1B and Level‐2 data. I thank Jürgen Kusche and Xiaoping Wu for clarification on the reference frame and the degree 1 coefficients. I thank Karl Bretreger for proofreading the manuscript. The constructive comments by Paul Tregoning (Editor), Associate Editor, Matt King, and an anonymous reviewer were helpful to improve the original manuscript substantially. GRACE data used in this study are available from http://podaac.jpl.nasa.gov/GRACE and http://grgs.obs‐mip.fr/grace. Two different GPS position solutions are obtained from http://geodesy.unr.edu and https://gipsy‐oasis.jpl.nasa.gov. SLR geocenter data are used from ftp.csr.utexas.edu/pub/slr/geocenter.