Kristel Chanard
email: kristel.chanard@ign.fr

Luce Fleitout

Eric Calais

Paul Rebischung

Jean-Philippe Avouac

Journal of Geophysical Research: Solid Earth

We model surface displacements induced by variations in continental water, atmospheric pressure, and nontidal oceanic loading, derived from the Gravity Recovery and Climate Experiment (GRACE) for spherical harmonic degrees two and higher. As they are not observable by GRACE, we use at first the degree-1 spherical harmonic coefficients from Swenson et al. (2008, https://doi.org/10.1029/2007JB005338). We compare the predicted displacements with the position time series of 689 globally distributed continuous Global Navigation Satellite System (GNSS) stations. While GNSS vertical displacements are well explained by the model at a global scale, horizontal displacements are systematically underpredicted and out of phase with GNSS station position time series. We then reestimate the degree 1 deformation field from a comparison between our GRACE-derived model, with no a priori degree 1 loads, and the GNSS observations. We show that this approach reconciles GRACE-derived loading displacements and GNSS station position time series at a global scale, particularly in the horizontal components. Assuming that they reflect surface loading deformation only, our degree-1 estimates can be translated into geocenter motion time series. We also address and assess the impact of systematic errors in GNSS station position time series at the Global Positioning System (GPS) draconitic period and its harmonics on the comparison between GNSS and GRACE-derived annual displacements. Our results confirm that surface mass redistributions observed by GRACE, combined with an elastic spherical and layered Earth model, can be used to provide first-order corrections for loading deformation observed in both horizontal and vertical components of GNSS station position time series.

Introduction

Global Navigation Satellite System (GNSS) station position time series show strong seasonal signals, in both horizontal and vertical components [START_REF] Blewitt | A new global mode of Earth deformation: Seasonal cycle detected[END_REF]. These signals primarily result from surface mass redistribution: nontidal oceanic loading, continental hydrology [START_REF] Dong | Anatomy of apparent seasonal variations from GPS-derived site position time series[END_REF]Van Dam et al., 2001), ice and snow [START_REF] Grapenthin | Icelandic rhythmic: Annual modulation of land elevation and plate spreading by snow load[END_REF][START_REF] Jiang | Accelerating uplift in the north Atlantic region as an indicator of ice loss[END_REF][START_REF] Matsuo | Time-variable ice loss in Asian high mountains from satellite gravimetry[END_REF], and atmospheric pressure [START_REF] Kaniuth | Estimating atmospheric pressure loading regression coefficients from GPS observations[END_REF].

Vertical displacements induced by surface loading can be well predicted by the response of an elastic spherical and layered Earth model to mass redistributions derived from the Gravity Recovery and Climate Experiment (GRACE) [START_REF] Bettinelli | Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology[END_REF][START_REF] Chanard | Modeling deformation induced by seasonal variations of continental water in the Himalaya region: Sensitivity to Earth elastic structure[END_REF][START_REF] Davis | Climate-driven deformation of the solid Earth from GRACE and GPS[END_REF]Fu & Freymueller, 2012;Fu et al., 2012[START_REF] Fu | Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE[END_REF][START_REF] Nahmani | Hydrological deformation induced by the West African monsoon: Comparison of GPS, GRACE and loading models[END_REF]. Yet horizontal displacements obtained using the same approach are unsatisfactory. Discrepancies between observed and modeled time series have been attributed to either the poor spatial resolution of GRACE loads estimates, resolving loads of wavelengths larger than 400 km [START_REF] Fu | Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE[END_REF][START_REF] Tregoning | Atmospheric effects and spurious signals in GPS analyses[END_REF] or local variations in the Earth's elastic properties [START_REF] Chanard | Modeling deformation induced by seasonal variations of continental water in the Himalaya region: Sensitivity to Earth elastic structure[END_REF][START_REF] Drouin | Constraints on seasonal load variations and regional rigidity from continuous GPS measurements in Iceland, 1997-2014[END_REF].

Providing an accurate global model for horizontal and vertical surface displacements induced by mass redistribution is of importance. First, there is currently no conventional model recommended to correct nontidal loading deformations in, for example, terrestrial reference frame determination. Second, such a model would help in improving long-term estimates of tectonic velocities [START_REF] Blewitt | GPS and space based geodetic methods[END_REF] and detecting transient Here we investigate how to accurately predict horizontal and vertical seasonal displacements at a global scale. After a description in section 2 of the GRACE and GNSS data used in this study, we model, in section 3, the global elastic surface deformation field induced by surface mass variations derived from GRACE, using the degree 1 spherical harmonics coefficients from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF]. We show that while vertical GNSS displacements are well explained by this model, there is a misfit in phase and amplitude in the horizontal components.

To improve the model, we reestimate, in section 4, the degree 1 deformation field from a global comparison between the displacements derived from GRACE, with no degree 1 loads, and the GNSS observations. Assuming that they reflect surface loading deformation only, the degree 1 coefficients estimated in this approach can be translated into geocenter motion time series, for example, the motion of the center of mass (CM) of the whole Earth system with respect to the center of figure (CF) of the solid Earth surface.

Finally, in the last part of section 4, we discuss the impact of systematic errors in GNSS station position time series on the prediction of seasonal observations. We investigate in particular the systematic errors occurring at the Global Positioning System (GPS) draconitic period and its harmonics, overlaying annual and semiannual geophysical signals for the first and second harmonic.

GRACE and GNSS Data Sets

Surface Load Mass Variations Derived From GRACE Level 2 Solutions

Measurements collected by the GRACE mission allow estimating global mass redistribution at the Earth's surface [START_REF] Tapley | GGM02-An improved Earth gravity field model from GRACE[END_REF]. In this study, we make use of the 10 day Level 2 solutions provided by the CNES/GRGS (http://grgs.obsmip.fr/). Solutions are expressed in terms of Stokes coefficients representing the gravitational effects of nonmodeled phenomena (continental water, sediments displacements, vegetation, oceanic, and atmospheric mass variations, etc.), then converted into geoid and surface mass coefficients, expressed in millimeter of equivalent-water height (EWH) [START_REF] Ramillien | Time variations of land water storage from an inversion of 2 years of GRACE geoids[END_REF]. Grids of EWH are computed from surface mass coefficients and corrected for a time average to express them with respect to the mean surface mass distribution over the time span of analysis. Among the gravitational contributions corrected in the CNES/GRGS solutions, atmospheric and nontidal oceanic loads [START_REF] Carrère | Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing-comparisons with observations[END_REF] are added back to enable comparison with GNSS station position time series (Figures S1 and S2 in the supporting information). The choice of the CNES/GRGS processing center is arbitrary but has no significant impact on the results of our study as differences between processing techniques for seasonal signals (Figure S3) range within the formal error bounds of the GRACE solution itself [START_REF] Sakumura | Ensemble prediction and intercomparison analysis of GRACE time-variable gravity field models[END_REF].

We use detrended time series of EWH from February 2002 to August 2012 as an estimate of the space-and time-varying load, including seasonal signals to first order, with smaller amplitude interannual variability in some regions. During this time interval, EWH time series are affected by large coseismic and postseismic gravity changes due to megaearthquakes [START_REF] Chen | GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earthquake[END_REF][START_REF] Han | Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake[END_REF][START_REF] Han | Implications of postseismic gravity change following the great 2004 Sumatra-Andaman earthquake from the regional harmonic analysis of GRACE intersatellite tracking data[END_REF][START_REF] Ogawa | Slow postseismic recovery of geoid depression formed by the 2004 Sumatra-Andaman earthquake by mantle water diffusion[END_REF]. We empirically model and remove the mass redistribution-induced megaearthquakes and observed in the GRACE data, as follows:

h s (m, t) = h(m, t) - ∑ i 𝛼 i (t -t i ) + ∑ i 𝛽 i (t -t i ) ln ( 1 + t -t i 𝜏 ) (1)
where h s refers to the EWH at point m at the Earth's surface, after earthquake contributions have been removed from the observed EWH h.  is the Heaviside function representing abrupt coseismic changes in the gravity field at time t i , and 𝜏 is a characteristic time of postseismic relaxation proper to each megaearthquake [START_REF] Hoechner | Investigation on afterslip and steady state and transient rheology based on postseismic deformation and geoid change caused by the Sumatra 2004 earthquake[END_REF]. The parameters 𝛼 i and 𝛽 i are constants estimated by a least-squares inversion. The main signal removed from the EWH is the coseismic offset induced by the 2004 Sumatra-Andaman earthquake (Figure S4). Figure 1 shows the peak-to-peak amplitude of the EWH (h s , in mm) over the 2002-2012 time period, once coseismic and postseismic effects have been removed.

As we will be using the Love number formalism to compute the deformation induced by surface mass variations (section 3.1), the EWH grids corrected for coseismic and postseismic effects are finally expanded into spherical harmonic coefficients from degrees 2 to 80. To insure comparison with GNSS station position time series, we first account for the degree 1 spherical harmonic load coefficients, not observable by GRACE, using results from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF].

GNSS Station Position Time Series

Our GNSS station position time series are based on the daily combined terrestrial frame solutions from the International GNSS Service second reprocessing campaign (Rebischung et al., 2016, IGS repro2). Discontinuities related to earthquakes, equipment changes, or other unknown causes were identified as detailed in Altamimi et al. (2016). Postseismic deformations were parametrically modeled and removed as in the ITRF2014 computation (Altamimi et al., 2016). A long-term stacking of the daily combined IGS repro2 solutions was then performed (where no daily scale parameters were estimated) where station positions were modeled as the sum of piecewise linear functions, and annual and semiannual sine waves. The GNSS station position time series used in the following are the residuals from this stacking, to which the estimated annual and semiannual signals were added back. They therefore represent nonlinear station motions, corrected for discontinuities and postseismic deformations with respect to the Center of the IGS08b Network (CN), which can be assumed to approximate the CF.

Among the 1,220 available station position time series, we select a subset of 689 stations based on the formal error of the amplitude of an annual signal estimated simultaneously with draconitic signals up to the sixth harmonic. The selected stations are those with formal errors less than 0.1 mm and 0.3 mm in the horizontal and vertical components, respectively. This selection criterion is directly correlated with the length of the time series (Figure S5). Locations of the selected stations are shown in Figure 1, while discarded stations are also shown in Figure S6.

Initial Modeling of Ground Deformation Induced by Surface Load Variations

Description of the Model

We compute the Earth's surface deformation induced by surface load variations using a numerical model based on a spherical harmonics decomposition of the loads and the Love numbers theory. The surface load 10.1002/2017JB015245

grids 𝜎(t, 𝜙, 𝜆), varying with time (t), longitude (𝜙), and latitude (𝜆) are decomposed into a sum of spherical harmonic coefficients (𝜎 C lm , 𝜎 S lm ) of degree and order (l, m) for each time t as follows:

𝜎(t, 𝜙, 𝜆) = ∞ ∑ l=1 l ∑ m=0 ∑ 𝜓∈{S,C} 𝜎 𝜓 lm (t)Y 𝜓 lm (𝜙, 𝜆) (2) (Y C lm , Y S lm )
are defined in terms of Legendre polynomials P lm as

( Y C lm (𝜙, 𝜆) Y S lm (𝜙, 𝜆) ) = √ (2l + 1)(l -m)! 4𝜋(l + m)! P lm (sin 𝜙) ( cos(m𝜆) sin(m𝜆) ) (3)
Surface displacements at a point P induced by the load at a given time t are obtained by solving a system of equations for the deformation of a self-gravitational spheroidal body, similarly to [START_REF] Farrell | Deformation of the Earth by surface loads[END_REF]. The obtained surface displacements are expressed as

dE(t, 𝜙, 𝜆) = 4𝜋R 3 E M E ∞ ∑ l=0 l ∑ m=0 ∑ 𝜓∈{S,C} l 1 2l + 1 1 cos 𝜙 𝜎 𝜓 lm (t) 𝜕Y 𝜓 lm 𝜕𝜆 (𝜙, 𝜆) dN(t, 𝜙, 𝜆) = 4𝜋R 3 E M E ∞ ∑ l=0 l ∑ m=0 ∑ 𝜓∈{S,C} l 1 2l + 1 𝜎 𝜓 lm (t) 𝜕Y 𝜓 lm 𝜕𝜙 (𝜙, 𝜆) dU(t, 𝜙, 𝜆) = 4𝜋R 3 E M E ∞ ∑ l=0 l ∑ m=0 ∑ 𝜓∈{S,C} h 1 2l + 1 𝜎 𝜓 lm (t)Y 𝜓 lm (𝜙, 𝜆) (4)
where (l 1 , h 1 ) are the tangential and radial load Love numbers, and (l 1 , h 1 ) depend on the reference frame considered. M E and R E are the Earth's mass and radius, respectively. In this study, we use a purely elastic Earth structure based on the Preliminary Reference Earth Model [START_REF] Dziewonski | Preliminary reference Earth model[END_REF], where the top 40 km are replaced by the CRUST 2.0 1-D average continental crust up to 40 km depth [START_REF] Bassin | The current limits of resolution for surface wave tomography in North America[END_REF].

Initial Model Results

Using equation (4), we compute the horizontal and vertical surface displacements at the 689 selected GNSS stations previously described. Figure 2 shows the east, north, and vertical detrended daily station position time series with 1 𝜎 error bars (gray) and their 10 days moving averages (black crosses) for three stations: LHAZ, Tibet (Bundesamt fuer Kartographie) GOLD, California, and BRAZ, Brazil (Jet Propulsion Laboratory). We refer the reader to Figure 1 for their locations. We show the predictions of our initial GRACE-derived model, where degree 1 coefficients from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF] have been used, in blue. While the model explains reasonably well the observed vertical seasonal displacements, the horizontal prediction is unsatisfactory in both phase and amplitude.

To quantify how well the model explains the GNSS observations at a global scale, we compute the WRMS at station i, for component j with N i observations, between the observations d with associated errors 𝜎 and the [START_REF] References Altamimi | ITRF2014: A new release of the international terrestrial reference frame modeling nonlinear station motions[END_REF] and the WRMS reduction (WRMSr), which compares the WRMS results for the model m with the null model:

√ √ √ √ ∑ N i k=1 ( d i,j,k -m i,j,k 𝜎 i,j,k ) 2 ∑ N i k=1 1 𝜎 2 i,j,k
WRMSr i,j = WRMS d-0 i,j -WRMS d-m i,j WRMS d-0 i,j (6)
Figures 3a,3c, and 3e show the WRMSr values obtained by comparing observations at the selected 689 GNSS stations with our initial model for the north, east, and vertical components, respectively. The initial model shows a significant disagreement with horizontal observations at a global scale, while the vertical component is fairly well predicted. We also report in Table 1 mean WRMSr values and the percentages of stations at which the WRMS is reduced by our initial model, confirming the discrepancy between the horizontal and vertical predictions of the observed GNSS signals. A common argument to explain this discrepancy is the spatial resolution of GRACE [START_REF] Fu | Horizontal motion in elastic response to seasonal loading of rain water in the Amazon Basin and monsoon water in Southeast Asia observed by GPS and inferred from GRACE[END_REF], a larger part of horizontal displacements being attributed to high spherical harmonics degree loads, that is, short wavelength loads, that GRACE cannot resolve. However, while the amplitudes of horizontal and vertical displacements in response to unit loads significantly differ for low spherical harmonic degrees, they evolve asymptotically in a similar manner for spherical harmonic degrees higher than 50 (loading wavelengths ≤ 800 km) as shown by Figure 4, indicating that a spatial resolution issue would affect both horizontal and vertical components similarly. [START_REF] Bassin | The current limits of resolution for surface wave tomography in North America[END_REF].

Reconciling GRACE-Derived Loading Deformation With GNSS Station Position Time Series

Understanding differences between the deformation model derived from GRACE and GNSS station position time series is of importance. Indeed, an accurate loading model for both horizontal and vertical components would have implications not only for correcting station position time series but could also provide insights on the nature of seasonal loads and/or help detecting systematic errors in geodetic observations. In this section, we estimate the degree 1 surface load coefficients from the GNSS deformation field itself, instead of using degree 1 loads from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF].

Degree 1 Deformation Approach

In the classical approach used to invert degree 1 deformation from a GNSS deformation field [START_REF] Blewitt | Self-consistency in reference frames, geocenter definition, and surface loading of the solid earth[END_REF][START_REF] Blewitt | A new global mode of Earth deformation: Seasonal cycle detected[END_REF], degree 1 surface load coefficients are estimated together with a net translation (and a net rotation) which aims to accommodate misalignments of the GNSS solution to the theoretical reference frame. The deformation induced by degree 1 surface load coefficients is given by equation (4) as follows:

⎛ ⎜ ⎜ ⎜ ⎜ ⎝ dE 1 (𝜙, 𝜆) dN 1 (𝜙, 𝜆) dU 1 (𝜙, 𝜆) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = 4𝜋R 3 E 3M E ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ l 1 cos 𝜙 𝜕Y C 11 𝜕𝜆 (𝜙, 𝜆) l 1 cos 𝜙 𝜕Y S 11 𝜕𝜆 (𝜙, 𝜆) l 1 cos 𝜙 𝜕Y 10 𝜕𝜆 (𝜙, 𝜆) l 1 𝜕Y C 11 𝜕𝜙 (𝜙, 𝜆) l 1 𝜕Y 10 𝜕𝜙 (𝜙, 𝜆) l 1 𝜕Y 10 𝜕𝜙 (𝜙, 𝜆) h 1 Y C 11 (𝜙, 𝜆) h 1 Y S 11 (𝜙, 𝜆) h 1 Y 10 (𝜙, 𝜆) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜎 C 11 𝜎 S 11 𝜎 10 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (7)
where:

⎛ ⎜ ⎜ ⎜ ⎜ ⎝ Y C 11 (𝜙, 𝜆) Y S 11 (𝜙, 𝜆) Y 10 (𝜙, 𝜆) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ cos𝜙 cos 𝜆 cos 𝜙 sin 𝜆 sin 𝜙 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (8)
and (h 1 , l 1 ) are the spherical harmonic degree 1 load Love numbers appropriate for the reference frame in which the displacements are measured. Combining equations ( 7) and ( 8), the displacement induced by spherical harmonic degree 1 surface loads can be expressed as follows:

⃗ d 1 (𝜙, 𝜆) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ dE 1 (𝜙, 𝜆) dN 1 (𝜙, 𝜆) dU 1 (𝜙, 𝜆) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = 4𝜋R 3 E 3M E ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ -l 1 sin 𝜆 l 1 cos 𝜆 0 -l 1 sin 𝜙 cos 𝜆 -l 1 sin 𝜙 sin 𝜆 l 1 cos 𝜙 h 1 cos 𝜙 cos 𝜆 h 1 cos 𝜙 sin 𝜆 h 1 sin 𝜙 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜎 C 11 𝜎 S 11 𝜎 10 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (9)
10.1002/2017JB015245

Besides, the displacement induced by a net translation ⃗ T(T X , T Y , T Z ) at a point (𝜙, 𝜆) of the Earth's surface can be expressed as:

⃗ d 2 (𝜙, 𝜆) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ dE 2 (𝜙, 𝜆) dN 2 (𝜙, 𝜆) dU 2 (𝜙, 𝜆) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ -sin 𝜆 cos 𝜆 0 -sin 𝜙 cos 𝜆 -sin 𝜙 sin 𝜆 cos 𝜙 cos 𝜙 cos 𝜆 cos 𝜙 sin 𝜆 sin 𝜙 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ T X T Y T Z ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (10)
The sum of both the estimated degree 1 deformation field and net translation

( ⃗ d 1 + ⃗ d 2 )
has a horizontal component and a vertical component equivalent to those which would be induced by a translation ⃗ T H and a translation ⃗ T V , respectively. ⃗ T H and ⃗ T V can be expressed in a geocentric Cartesian frame as follows:

⃗ T H = 4𝜋R 3 E 3M E l 1 ⃗ 𝜎 1 + ⃗ T ( 11
) ⃗ T V = 4𝜋R 3 E 3M E h 1 ⃗ 𝜎 1 + ⃗ T ( 12
)
where

⃗ 𝜎 1 = (𝜎 C 11 , 𝜎 S 11 , 𝜎 10 ). As ( ⃗ T H , ⃗ T V )
is a linear combination of ( ⃗ 𝜎 1 , ⃗ T), the classical degree 1 deformation approach [START_REF] Blewitt | Self-consistency in reference frames, geocenter definition, and surface loading of the solid earth[END_REF][START_REF] Blewitt | A new global mode of Earth deformation: Seasonal cycle detected[END_REF] is exactly equivalent to estimating two translations from the horizontal and vertical station displacements independently. The two approaches being equivalent, we choose to estimate separately a translation (T V X , T V Y , T V Z ) from the vertical displacements, and a translation (

T H X , T H Y , T H Z ) plus a rotation (R H X , R H Y , R H Z )
from the horizontal displacements, rather than using the classical deformation approach, to emphasize the differences in the degree 1 effects on the horizontal and vertical components.

Although geocenter motion is not the goal of our study, we note that, using equations ( 11) and ( 12), degree 1 coefficients can be retrieved from

( ⃗ T H , ⃗ T V ) by ⃗ 𝜎 1 = M E 4𝜋R 3 E ⃗ T V -⃗ T H h 1 -l 1 ( 13
)
and converted into geocenter motion [START_REF] Lavallée | Geocenter motions from GPS: A unified observation model[END_REF], for example, the motion of the CM of the whole Earth system with respect to the CF of the solid Earth surface, assuming that they reflect surface loading only, by

Δ⃗ r CM-CF = ( h CE 1 + 2l CE 1 3 -1 ) ⃗ T V -⃗ T H h 1 -l 1 = 2.53 × ( ⃗ T V -⃗ T H ) ( 14
)
where the difference between degree 1 Love numbers is independent of the reference frame [START_REF] Blewitt | Self-consistency in reference frames, geocenter definition, and surface loading of the solid earth[END_REF].

Adjusted Model Results

We invert the differences between the GNSS station position time series and our initial GRACE-derived model without degree 1 coefficients for translations and rotation parameters, as described in the previous section. Note that we do not invert for a scale parameter that could bias seasonal signals [START_REF] Collilieux | Strategies to mitigate aliasing of loading signals while estimating GPS frame parameters[END_REF]. The inversion is performed for all available stations at each time step (Figure S7). The estimated parameters are shown in Figure 5. Both horizontal and vertical translations include clear annual oscillations, particularly in the Z component, for which both times series are out of phase. Differences between the horizontal and vertical translations highlight the distinct effect of degree 1 loads on both components. The estimated rotations remain small. The geocenter motion time series deduced from the estimated translations are shown in Figure S8. We used a diagonal covariance matrix in our study and not the full network covariance matrix. However, this choice has little influence on the derived geocenter motion that remain small (Figure S9). The annual components of geocenter motion derived in this study are within the range of other geodetic estimates (Wu et al., 2012), with annual amplitudes of 1.4, 2.2, and 6.1 mm in the X, Y, and Z components, respectively. The annual amplitude in Z is however at the high end of geodetic estimates, which could be due to the estimated degree 1 deformation absorbing other effects besides loading deformation, such as the response of the bedrock and GNSS monument to seasonal temperature variations.

We then add the effect of the estimated translations and rotations to our initial GRACE-derived model to account for adjusted degree 1 contributions. Figure 2 shows the north, east, and vertical detrended daily station position time series with 1 𝜎 error bars (gray) and their 10 day moving averages (black crosses) Global Navigation Satellite System stations for horizontal (red) and vertical (blue) components. The number of sites available at each time step is given in Figure S7). Time series of geocenter motion can be derived from the difference between horizontal and vertical translations (see Figure S8).

for three stations (LHAZ, BRAZ, and GOLD), the predictions of our initial model (blue curve), where degree 1 coefficients from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF] have been used, and of our adjusted model (red curve). Compared to the initial model, the adjusted model significantly improves the fit to the GNSS observations, particularly in the horizontal components. shows significantly better WRMS reductions at most stations compared to the initial model. Mean values of WRMS reduction and numbers of stations at which the WRMS is reduced are given in Table 1, confirming a significant global improvement when the degree 1 contribution is adjusted.

We also compare the purely annual signals by fitting a sine function to both the GNSS and GRACE-derived time series and computing the annual WRMS reduction. Results are given in Table 1 and show that our adjusted model reduces the annual WRMS at more than 80% of the stations for all three components, with improved mean values particularly for the horizontal components. Figure 6 shows histograms of the annual WRMS reduction for all 689 GNSS stations, confirming that the adjusted model predicts the annual component of the data considerably better than the initial model, particularly for the east and north components.

While the annual WRMS of the vertical signals are reduced by 55.4% by our adjusted model, the east and north annual signals are only reduced by 16.5% and 35.0%, respectively. The discrepancy between the east and north WRMS reduction is most likely due to spurious north-south stripes arising from the GRACE processing.

Although our adjusted model is a significant improvement over the initial model, as shown in Figure 6, it still misses some of the observed annual signal, particularly in the horizontal components. Unmodeled signals could arise from physical processes not considered in this study, like thermal effects, for example, and also from systematic errors in GNSS station position time series. In particular, the known systematic errors at the first and second draconitic harmonics, of periods 351.6 and 175.8 days (Amiri-Simkooei, 2013), and present in the IGS products [START_REF] Ray | Anomalous harmonics in the spectra of GPS position estimates[END_REF] close to the annual and semiannual periods could explain part of the unmodeled seasonal signals. The source of the draconitic errors has not been clearly identified yet, but it could be due to orbit modeling deficiencies, in particular, due to the Sun-satellite interactions or during eclipse seasons [START_REF] Hugentobler | Identification and mitigation of GNSS errors[END_REF][START_REF] Ray | Anomalous harmonics in the spectra of GPS position estimates[END_REF][START_REF] Rodriguez-Solano | Reducing the draconitic errors in GNSS geodetic products[END_REF], aliasing and/or propagation of site-dependent effects such as multipath or errors in the antenna or radome calibrations [START_REF] Ray | Anomalous harmonics in the spectra of GPS position estimates[END_REF], mismodeling of diurnal and semidiurnal ocean tides, and neglecting of semidiurnal and diurnal atmospheric tides [START_REF] Tregoning | Atmospheric effects and spurious signals in GPS analyses[END_REF], 2011). Whatever their origin, the draconitic errors result in spurious periodic signals that could bias the estimates of seasonal variations of station positions.

10.1002/2017JB015245

To assess the impact of draconitic errors on the comparison of GNSS and GRACE-derived annual displacements, we estimate annual, semiannual, and draconitic signals up to the sixth harmonic from the 689 GNSS station position time series and remove the estimated draconitic signals from the time series. We then reestimate the degree 1 contribution using our two translations approach and compute the annual WRMS reduction of the new adjusted model where the draconitic estimates have been removed from the GNSS station position time series. Figure 7 shows histograms of the annual WRMS reduction obtained with the adjusted model and the new adjusted model where we estimated and removed draconitic errors. We find that estimating and removing draconitic signals does not improve the fit to observations, indicating that the unmodeled signals probably arise from other physical processes. Maps of the observed and predicted annual signals are given in Figure S10, with and without draconitic corrections.

Discussion and Conclusions

Based on an elastic Earth model, we derived the horizontal and vertical displacements induced by the surface mass variations observed by GRACE. These surface mass variations are primarily related to continental water storage, atmospheric pressure, and oceanic circulation. We then compared the predicted displacements with the position time series of 689 globally distributed GNSS stations. We showed that estimating the degree 1 deformation field rather than using the degree 1 coefficients from [START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF] significantly improves the agreement between GNSS and GRACE-derived displacements, particularly in the horizontal components, although misfits still remain. Several effects may contribute to this remaining misfit. Systematic errors in geodetic products may explain part of the misfit. Small-scale loads unresolved by GRACE may also affect GNSS sites at a local scale. Other unmodeled physical effects may finally contribute to explain the remaining discrepancies between GNSS observations and the predictions of our loading model. For instance, the thermal expansion of bedrock likely contributes to the seasonal signals observed in GNSS station position time series [START_REF] Prawirodirdjo | Observation and modeling of thermoelastic strain in Southern California integrated GPS network daily position time series[END_REF]Tsai, 2011), as well as the thermal expansion of the monuments. Regional misfits may also arise from lateral variations of the Earth's elastic parameters, although to a lesser extent than proposed by [START_REF] Chanard | Modeling deformation induced by seasonal variations of continental water in the Himalaya region: Sensitivity to Earth elastic structure[END_REF] once degree 1 contribution has been adjusted. Finally, the hypothesis of a purely elastic Earth model derived from seismic waves (with periods ≤1 h) may not hold at an annual time scale. Indeed, elastic Earth models provide accurate predictions of surface displacements or induced gravity perturbations for periods up to a month, but not for longer periods [START_REF] Benjamin | Constraints on mantle anelasticity from geodetic observations, and implications for the J2 anomaly[END_REF]. At these longer periods, it is now widely accepted that the mantle must be considered as viscoelastic: for instance, the 2010 IERS conventions [START_REF] Petit | IERS conventions[END_REF] introduce viscoelasticity through complex, frequency-dependent, Love numbers. Both [START_REF] Benjamin | Constraints on mantle anelasticity from geodetic observations, and implications for the J2 anomaly[END_REF] and [START_REF] Krásná | Tidal Love and Shida numbers estimated by geodetic VLBI[END_REF], however, suggest that the viscoelastic corrections in the IERS conventions might be insufficient and proposed values for the complex degree 2 tidal Love numbers significantly different. The degree 2 tidal Love numbers are mainly sensitive to the mechanical properties of the lower mantle, while the load Love numbers of higher degrees, which play a major role in the Earth's response to nontidal loads, are mainly sensitive to the rheology of the upper mantle. The validity of using a purely elastic Earth model at an annual time scale, while beyond the scope of this paper, should therefore be addressed in future studies [START_REF] Chanard | Constraints on transient viscoelastic rheology of the asthenosphere from seasonal deformation[END_REF].

Figure 1 .

 1 Figure 1. Peak-to-peak surface load variations, expressed in equivalent-water height (in mm), derived from Gravity Recovery and Climate Experiment for the 2002-2012 period and corrected from detectable coseismic and postseismic contributions. Red dots show the locations of the continuous Global Positioning System stations used in this study. LHAZ, GOLD, and BRAZ are highlighted because time series for these three stations are shown as examples later on.

Figure 2 .Figure 3 .

 23 Figure 2. Daily detrended Global Navigation Satellite System station position time series and associated 1 𝜎 error bars (gray), for stations (a) LHAZ, China (Tibet), (b) BRAZ, Brazil, and (c) GOLD, USA (California) (see locations in Figure1). Black crosses are 10 days moving average of the daily observations. The blue and red lines show the displacements derived from Gravity Recovery and Climate Experiment, where the degree 1 contribution is either taken from[START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF] (blue) or estimated from the Global Navigation Satellite System deformation field (red).

Figure 4 .

 4 Figure 4. Vertical (blue) and horizontal (red) amplitude of surface displacement as a function of spherical harmonic number, in response to a unit 1 year periodic unit harmonic loading function acting on a spherical and layered model for Earth, based on the Preliminary Reference Earth model, proposed by[START_REF] Dziewonski | Preliminary reference Earth model[END_REF] but where the top 40 km have been replaced by a continental crust[START_REF] Bassin | The current limits of resolution for surface wave tomography in North America[END_REF].

Figure 5 .

 5 Figure 5. Time series of inverted translation and rotation parameters derived from observations at a subset of 689 Global Navigation Satellite System stations for horizontal (red) and vertical (blue) components. The number of sites available at each time step is given in FigureS7). Time series of geocenter motion can be derived from the difference between horizontal and vertical translations (see FigureS8).

Figure 6 .

 6 Figure 6. Histograms of annual WRMS reduction in the (a) east, (b) north, and (c) vertical components obtained with our initial Gravity Recovery and Climate Experiment-derived model (degree 1 from Swenson et al. (2008); blue) and with our adjusted model (red).

Figure 7 .

 7 Figure 7. of annual weighted root mean square (WRMS) reduction in the (a) east, (b) north, and (c) vertical components obtained with our adjusted model, when estimating and removing draconitic errors from the Global Navigation Satellite System station position time series (green) or not (red).

Table 1

 1 Weighted Root Mean Square (WRMS) Reduction Statistics From the Comparison of 689 Globally Distributed GNSS StationPosition Time Series and GRACE-Derived Displacements, Where the Degree 1 Contribution Is Either Taken From[START_REF] Swenson | Estimating geocenter variations from a combination of GRACE and ocean model output[END_REF] or Estimated From the GNSS Deformation Field

		East		North		Vertical
		Stations	Mean	Stations	Mean	Stations	Mean
	WRMSr
	Initial	46.8%	-0.09%	47.4%	1.24%	55.6%	9.02%
	Adjusted	82.5%	7.94%	83.4%	14.3%	70.8%	14.7%
	Annual WRMSr
	Initial	56.6%	3.32%	54.1%	5.73%	88.3%	50.3%
	Adjusted	74.4%	16.5%	83.7%	35.0%	90.4%	55.4%
	Annual WRMSr-draconitics removed
	Initial	52.6%	-0.73%	52.7%	-1.17%	85.6%	40.4%
	Adjusted	72.1%	14.2%	80.1%	30.9%	88.3%	48.7%

Note. The percentages of stations with reduced WRMS and the mean WRMS reductions are given for the full time series (WRMSr), and for their annual components in the two cases where draconitic errors were estimated and removed (Annual WRMSr-draconitics removed) or not (Annual WRMSr). GNSS = Global Navigation Satellite System; GRACE = Gravity Recovery and Climate Experiment.

Acknowledgments

This work takes advantage of GRACE data processed by the CNES/GRGS (http://grgs.obsmip.fr/) and the IGS data (Dow et al., 2009). The project was funded by NSF grant EAR 1345136, the Laboratoire de Recherche Commun "Yves Rocard" (ENS-CEA-CNRS), and CNRS/TOSCA grant 2925.