Search for a command to run...
Seasonal variations in GNSS coordinate time series are largely driven by environmental mass loadings, which involve not only instantaneous elastic deformation but also time-dependent poroelastic and hydrological responses. Assuming an instantaneous Earth response may bias long-term deformation estimates and affect reference frame stability. We select 4,711 global vertical GNSS coordinate time series spanning at least a decade with > 95% data completeness, together with their corresponding environmental loading deformation time series, both provided by the Nevada Geodetic Laboratory. To address the combined effects of multiple environmental loads, we first identify the dominant loading component at each station — defined as the one contributing more than 75% of the total modeled seasonal deformation amplitude. Then we extract annual signals from both GNSS and the dominant loading deformation time series using Singular Spectrum Analysis and estimate the phase lags between them through cross-correlation algorithm. This approach minimizes interference from secondary loading sources and provides a clean estimate of the time lag associated with the primary driving process. Globally, hydrological loading induces phase delays with significant spatial variability (standard deviation: 41 days). These phase lags exhibit systematic spatial patterns, possibly reflecting diverse hydrological processes across regions: negative delays (GNSS lagging load by ~ 28 days) in the mountainous western United States are possibly associated with unmodeled deep subsurface water retention; widespread positive delays (GNSS leading load by ~ 13 days) in the U.S. Corn Belt and Europe suggest rapid water removal due to human activity; and anti-phase anomalies (difference of > 150 days) in confined aquifers may reflect poroelastic responses. Applying the phase-lag correction by subtracting the time-shifted annual loading signals from GNSS observations improves the consistency between them at ~ 70% of the stations compared to the standard instantaneous correction, with a median reduction in annual signal power increasing from 40.4% to 62.1% (69.2% to 85.2% in regions with strong hydrological loading like the Amazon and the mountainous western U.S.) and a concurrent increase in the median weighted root-mean-square (WRMS) reduction from 2.7% to 4.0% (7.2% to 9.7% in the aforementioned regions). It also mitigates potential bias in site-specific velocity estimates (up to 0.04 mm/yr). Our results demonstrate that accounting for phase lags between GNSS observations and loading models is important for refining loading corrections and thereby enhancing the stability of geodetic reference frames.