Abstract: Objectives Seasonal signals in GNSS coordinate time series exhibit pronounced non-stationary characteristics, as their amplitudes undergo slow and continuous modulation driven by the interannual variability of environmental loading. Conventional methods such as LS-HEM assume fixed seasonal amplitudes throughout the observation period, resulting in systematic residuals when interannual modulation is significant, while SSA reflects amplitude evolution only indirectly through reconstructed components without explicit parametric representation. This study proposes a Gaussian process time-varying amplitude modeling method based on a product kernel, termed GP-MP, to overcome these structural limitations.Methods GP-MP constructs the covariance function as the pointwise product of a Matern 3/2 kernel and an exponential sine-squared periodic kernel. The periodic kernel constrains the seasonal phase structure, while the Matern kernel governs smooth amplitude modulation, enabling the effective seasonal amplitude to evolve continuously over time without any parametric assumption on its trajectory. All model hyperparameters are jointly optimized by maximizing the log marginal likelihood. The method was validated through controlled simulation experiments using 10-year synthetic daily series with known timevarying amplitude profiles, and through real-data analysis of vertical displacement time series from 46 continuous GNSS stations on the northeastern margin of the Tibetan Plateau, processed with GAMIT/GLOBK in the ITRF2014 reference frame.Results In simulation experiments, GP-MP reduced fitting RMSE by 37% relative to LS-HEM and by 18% relative to SSA, and improved missing-data reconstruction accuracy by 26% and 22% respectively, with residuals closely approximating white noise. Under real-data conditions across 46 stations, mean fitting RMSE improved by 32% and 27% over LS-HEM and SSA. At representative stations exhibiting complex signals such as abrupt displacement drops and short-duration amplitude perturbations, GP-MP consistently outperformed both competing methods. For regional velocity field estimation, the mean vertical velocity derived from GP-MP was 30 to 40 percent lower than that from LS-HEM, indicating that unresolved non-stationary seasonal signals are aliased into the long-term trend estimate under constant-amplitude modeling.Conclusions GP-MP provides a principled solution to non-stationary seasonal amplitude modeling in GNSS coordinate time series. The product-kernel structure explicitly encodes amplitude modulation at the covariance function level, achieving consistent accuracy improvements over LS-HEM and SSA in both signal fitting and missing-data reconstruction. Applied to velocity field estimation, GP-MP suppresses the aliasing of seasonal non-tectonic signals into long-term trend estimates, yielding more reliable velocity fields for regional crustal deformation monitoring.