GNSS Stations Research Articles

RADARSAT-2 DInSAR and GNSS-Derived Finite Fault Model of the 2012 Mw 7.8 Haida Gwaii Earthquake

The Mw 7.8 Haida Gwaii earthquake occurred on 28 October 2012, generating up to 13 m tsunami waves and 3 m run-up along the British Columbia coastline. Despite the magnitude of the earthquake and tsunami, damages were minor due to the lack of vulnerable infrastructure in the remote area. Previous finite fault models were derived from GNSS, seismic and tsunami data, but the uncertainty remained high due to the limited number of seismic and GNSS stations near the epicenter. In this study, finite fault models were developed using RADARSAT-2 interferograms and previously published GNSS data. These models defined the location of the fault and provided a detailed slip distribution with a high degree of certainty. The results confirmed that the main slip was located on the subduction fault interface between the Pacific and North American plates, west of the Queen Charlotte Fault. The estimated moment magnitude of 7.88–7.93 is slightly larger than the previously reported moment magnitude of around 7.8, due to the capture of postseismic deformation in the interferograms. Overall, the study provides an improved finite fault slip model for the Haida Gwaii earthquake and highlights the importance of utilizing remote sensing data for studying earthquakes in remote areas.

Impact of Tonga volcanic eruption on the ionosphere based on ground-based GNSS and COSMIC occultation data

The eruption of large volcanoes may potentially disturb the upper atmosphere, causing disturbances in ionospheric plasma and triggering irregular changes in the Total Electron Content (TEC). The GNSS ionospheric inversion method has been used to analyze ionospheric anomalies caused by the volcano, while it was limited to the spatial distribution of ground-based GNSS. The joint observation with space-based GNSS can further improve the spatio-temporal distribution of ionospheric monitoring. The Tonga volcano erupted on January 15th, 2022, is the largest volcanic eruption since the beginning of the 21st century and generated intense acoustic, atmospheric gravity waves and ionospheric disturbance. In this paper, we used the GIM products, ground-based GNSS station data around the volcanic eruption point and COSMIC occultation data to analyze the short-term TEC change. The experiments showed that: the average TEC of GIM grid data on eruption day dropped by 7.5 TECu compared to the last day which was substantially larger than other three days. The TEC extracted from nearby tracking stations was smaller than GIM and the difference was inversely proportional to the distance between the station and the volcanic eruption site. The average daily ionospheric TEC obtained from the nearest TONG station was the biggest, as much as 3.44 TECu lower. The peak electron density of COSMIC occultation data gradually diminished and then recovered with TEC about 7.32 TECu lower than GIM products. The experments showed that there was a good consistency in the TEC obtained from space-based and ground-based GNSS, the Tonga volcano produced a huge TEC change over 10 TECu or even bigger. As a global grid product, GIM underestimated ionospheric changes due to its limited spatial resolution, which needs further improvement.

Strain and Deformation Analysis Using 3D Geological Finite Element Modeling with Comparison to Extensometer and Tiltmeter Observations

This study verifies the practicality of using finite element analysis for strain and deformation analysis in regions with sparse GNSS stations. A digital 3D terrain model is constructed using DEM data, and regional rock mass properties are integrated to simulate geological structures, resulting in the development of a 3D geological finite element model (FEM) using the ANSYS Workbench module. Gravity load and thermal constraints are applied to derive directional strain and deformation solutions, and the model results are compared to actual strain and tilt measurements from the Jiujiang Seismic Station (JSS). The results show that temperature variations significantly affect strain and deformation, particularly due to the elevation difference between the mountain base and summit. Higher temperatures increase thermal strain, causing tensile effects, while lower temperatures reduce thermal strain, leading to compressive effects. Strain and deformation patterns are strongly influenced by geological structures, gravity, and topography, with valleys experiencing tensile strain and ridges undergoing compression. The deformation trend indicates a southwestward movement across the study area. A comparison of FEM results with ten years of strain and tiltmeter data from JSS reveals a strong correlation between the model predictions and actual measurements, with correlation coefficients of 0.6 and 0.75 for strain in the NS and EW directions, and 0.8 and 0.9 for deformation in the NS and EW directions, respectively. These findings confirm that the 3D geological FEM is applicable for regional strain and deformation analysis, providing a feasible alternative in areas with limited GNSS monitoring. This method provides valuable insights into crustal deformation in regions with sparse strain and deformation measurement data.

Analysis of the Response Relationship Between PWV and Meteorological Parameters Using Combined GNSS and ERA5 Data: A Case Study of Typhoon Lekima

Precipitable water vapor (PWV) is a crucial parameter of Earth’s atmosphere, with its spatial and temporal variations significantly impacting Earth’s energy balance and weather patterns. Particularly during meteorological disasters such as typhoons, PWV and other meteorological parameters exhibit dramatic changes. Studying the response relationship between PWV and typhoon events, alongside other meteorological parameters, is essential for meteorological and climate analysis and research. To this end, this paper proposes a method for analyzing the response relationship between PWV and meteorological parameters based on Wavelet Coherence (WTC). Specifically, PWV and relevant meteorological parameters were obtained using GNSS and ERA5 data, and the response relationships between PWV and different meteorological parameters before and after typhoon events were studied in time–frequency domain. Considering that many GNSS stations are not equipped with meteorological monitoring equipment, this study interpolated meteorological parameters based on ERA5 data for PWV retrieval. In the experimental section, the accuracy of ERA5 meteorological parameters and the accuracy of PWV retrieval based on ERA5 were first analyzed, verifying the feasibility and effectiveness of this approach. Subsequently, using typhoon Lekima as a case study, data from six GNSS stations affected by the typhoon were selected, and the corresponding PWV was retrieved using ERA5. The WTC method was then employed to analyze the response relationship between PWV and meteorological parameters before and after the typhoon’s arrival. The results show that the correlation characteristics between PWV and pressure can reveal different stages before and after the typhoon passes, while the local characteristics between PWV and temperature better reflect regional precipitation trends.

Modelling geoid height errors for local areas based on data of global models

Abstract The development of global geoid models became feasible following the launch of specialised satellite missions. Today, the root mean square deviation of the heights in global models of high degree & order varies from centimetres to decimetres across different countries. In countries where the accuracy of such models is lower, there is potential to enhance their precision by applying specific corrections. This study presents a novel methodology for locally modelling the height errors of high degree & order global geoid models using levelling sub-benchmarks for GNSS stations. The methodology is based on a combination of optimal interpolation methods, filtering, and the concept of data weighting by gravity anomaly differences. The methodology is aimed at creating a hybrid model that aligns with the local characteristics of the geoid (or quasi-geoid) derived from the traditional levelling network. The advantage of this methodology lies in its ability to reduce the residual height errors of the EGM2008 and EIGEN6C4 models to less than 1 cm when using only four control points. Such results exceed the initial values of the systematic height errors of these models by 90–96 %. For the GECO model, the residual errors are around 2 cm, while for the XGM2019e_2159 model, they reach 3 cm. These results indicate that this methodology can be applied to all global models of high degree & order, although its effectiveness may vary depending on the specifics of a particular model.

An empirical tool for predicting the presence or absence of coseismic displacements at GNSS stations

Unmodeled displacements in GNSS times series, induced by instrumental artifacts or geophysical events, create significant biases in station trajectory parameters that can propagate into the reference frame itself. While non-tectonic ‘jumps’, such as equipment changes, affect only a specific GNSS station, seismically-induced displacements can affect large numbers of sites, severely threatening the frame’s stability. Manually reviewing individual GNSS time series for such effects is highly impractical because there can be thousands of GNSS stations in a frame, and the total number of earthquakes Mw ≥ 6.0 since GPS became fully operational is + 5100. To avoid this time-consuming task, automated methods rely on empirical power-law functions to determine which earthquake-station pairs require coseismic displacement parameters. Still, ‘conservative’ power-law functions tend to add coseismic offsets to stations that do not need them, which can also threaten the stability of the frame. In this work, we present an empirical formulation that was obtained using 809 global seismic events to fit power-law parameters that do not overestimate the region of influence of earthquakes. Our method is based on a two-level selection process: level 1 is isotropic and only considers the epicentral distance between the stations and the earthquake, and level 2 uses the geophysical parameters of the earthquake to predict a ‘tighter’ displacement pattern to select which stations require coseismic trajectory parameters. We applied our level 2 method to a database of ~ 4700 event-station pairs and showed that it removed ~ 55% of the total pairs, all of which had been falsely selected by level 1.

The analysis of the scatters of the filtered residuals of the GNSS stations from the crustal movement observation network of China

Through the analysis of the coordinate time series of the stations in the Crustal Movement Observation Network of China (CMONOC) from 2011 to 2019, it is found that the daily network scatters of the filtered residuals (called as SFR) of the station coordinate time series have obvious seasonal variation characteristics which reaches a maximum during summer and a minimum during winter. This phenomenon has not yet been further explained. Our analysis results demonstrate that the seasonal pattern of the daily SFR solutions is the temperature-dependent, such as small-scale atmospheric perturbation. The pattern of SFR series of the difference between the zenith wet delay (ZWD) calculated with GPS data and that calculated with meteorological data is similar to the SFR solutions computed using GPS data. The results mean that the estimated ZWD using GPS data does not fully capture the total real atmospheric variation. The atmospheric delay model in current GNSS analysis can effectively solve large-scale atmospheric delay estimation, but it may fail to deduct certain small-scale atmospheric disturbances. We think that a considerable portion of these small-scale disturbances are absorbed by the estimates of station positions.

Estimation of Velocities of Ukrainian GNSS Stations in the IGb08 Reference Frame

Estimation of Velocities of Ukrainian GNSS Stations in the IGb08 Reference Frame

Evaluation of the high-degree global gravity field models in the territory of Kazakhstan

The study presents the results of evaluating the accuracy of the high-degree global gravity field models XGM2019e_2159, SGG-UGM-2, SGG-UGM-1, ЕGM2008, GECO and EIGEN-6C4 in the territory of the Republic of Kazakhstan. Based on the results of the study, the optimal model has been determined which best describes the long-wavelength component of the geoid model of the Republic of Kazakhstan. The results of astronomic measurements made at 59 Laplace stations and the results of the GNSS measurements carried out at 154 high-precision levelling in the framework of an investment project of the National Spatial Data Infrastructure of the Republic of Kazakhstan served as basic information for the study. The study results suggest that the high-degree global geopotential models have an approximately common modelling error, the standard deviation of around 11-12 cm, while the range of deviations in the model values of the height anomalies from those obtained by terrain measurements at the GNSS stations is 68-72 cm, which have a negative average value from –0.097 to –0.123 m. This may indicate a regular drift in the system of normal heights and the difference between the real W0 and the normal U0 of potentials at the beginning of height calculation not being equal to zero. The standard deviation of the component of deflection of the vertical in the plane of meridian Δ in the regions with the altitudes less than 500 m varies from 0.93 to 1.13 arcseconds, in the regions with the altitudes greater than 500 m, from 1.26 to 1.54 arcseconds, while the standard deviation of the component of the deflection of the vertical in the plane of the first vertical Δ in the regions with the altitudes less than 500 m varies from 0.70 to 0.84 arcseconds, in the regions with the altitudes greater than 500 m, from 0.62 to 0.81 arcseconds. Based on the statistical characteristics of the evaluatons of the accuracy of modeling height anomalies obtained from 154 stations of high-precision levelling and on the components of the deflection of the vertical from 59 Laplace stations, the global ultra high-degree geopotential model EIGEN-6C4 has proven to be the optimal global gravity field model which best describes the long-wave component of the geoid model for the Republic of Kazakhstan.

Three-dimensional interseismic crustal deformation in the northeastern margin of the Tibetan Plateau using GNSS and InSAR

The nature of the crustal deformation of the northeastern margin of the Tibetan Plateau is vital for elucidating the expansion mechanism of the plateau. We installed 21 continuous GNSS stations and obtained a horizontal velocity field with high spatial resolution. We also acquired a line-of-sight(LOS)velocity field with Sentinel-1 images covering the study area from 2014 to 2022. A multi-scale spherical wavelet method was employed to unify the reference frames of GNSS and InSAR data. After unifying the reference frame, the two data types achieve high consistency. Combining GNSS and InSAR velocities yielded the three-dimensional interseismic velocity field in the northeastern margin of the Tibetan Plateau. Furthermore, we analyzed the crustal deformation characteristics based on the three-dimensional deformation field. The crustal deformation exhibits a pronounced northeastward shift relative to the Ordos block, whereas the interior of the Ordos block remains remarkably stationary, behaving as a rigid unit. The Liupanshan fault is primarily characterized by uplift, devoid of any notable horizontal deformation across the fault. The left-lateral strike-slip mainly characterizes the Haiyuan fault.. The maximum east–west deformation velocity on the southern side of the fault is approximately 3.9 mm/yr, decreasing to about 1.0 mm/yr at the eastern end of the fault. The western segment of the West Qinling fault exhibits a minor east–west motion. Our result provides essential data for further study of the crustal deformation patterns.

Datum Alignment Between GNSS-IR Sea Level Estimations and Tide Gauges in Türkiye: A Vertical Local Tie Approach

Connecting Asia and Europe, Türkiye is vital for maritime transit, which requires continuous and precise sea level monitoring. Traditional tide gauges (TGs), while valuable, often have data gaps. GNSS-IR provides an alternative, though aligning it with TG measurements is challenging. This study introduces a method to estimate Vertical Local Tie (VLT) between the TG zero and the GNSS antenna phase centre, allowing GNSS-IR sea level estimates to be aligned with TG measurements. In disruptions, GNSS-IR can consistently supplement TG data. We validated this approach by evaluating a 1-year data from four GNSS stations namely, ARSZ, GADA, TRBZ, and YLVA involving in the Türkiye National Sea Level Monitoring System (TUDES) network. These stations are located on the coasts of the Mediterranean, Aegean, Black, and Marmara Seas, respectively. Compared with the TG measurements, ARSZ achieved the highest epoch-based correlation as 87.00%. For daily averages, the GADA and TRBZ showed a correlation of over 96% with TG data. Notably, TRBZ recorded the lowest daily average RMSE at 2.2 cm using the VLT-median approach. The results demonstrate that when VLT is utilized, well-configured GNSS stations, ensuring uninterrupted transmissions and appropriate data intervals, provide a reliable alternative or complement to TGs for sea level monitoring.

Improving GNSS-IR Sea Surface Height Accuracy Based on a New Ionospheric Stratified Elevation Angle Correction Model

Approximately 71% of the Earth’s surface is covered by vast oceans. With the exacerbation of global climate change, high-precision monitoring of sea surface height variations is of vital importance for constructing global ocean gravity fields and preventing natural disasters in the marine system. Global Navigation Satellite System Interferometry Reflectometry (GNSS-IR) sea surface altimetry is a method of inferring sea surface height based on the signal-to-noise ratio of satellite signals. It enables the retrieval of sea surface height variations with high precision. However, navigation satellite signals are influenced by the ionosphere during propagation, leading to deviations in the measured values of satellite elevation angles from their true values, which significantly affects the accuracy of GNSS-IR sea surface altimetry. Based on this, the contents of this paper are as follows: Firstly, a new ionospheric stratified elevation angle correction model (ISEACM) was developed by integrating the International Reference Ionosphere Model (IRI) and ray tracing methods. This model aims to improve the accuracy of GNSS-IR sea surface altimetry by correcting the ionospheric refraction effects on satellite elevation angles. Secondly, four GNSS stations (TAR0, PTLD, GOM1, and TPW2) were selected globally, and the corrected sea surface height values obtained using ISEACM were compared with observed values from tide gauge stations. The calculated average Root Mean Square Error (RMSE) and Pearson Correlation Coefficient (PCC) were 0.20 m and 0.83, respectively, indicating the effectiveness of ISEACM in sea surface height retrieval. Thirdly, a comparative analysis was conducted between sea surface height retrieval before and after correction using ISEACM. The optimal RMSE and PCC values with tide gauge station observations were 0.15 m and 0.90, respectively, representing a 20.00% improvement in RMSE and a 4.00% improvement in correlation coefficient compared to traditional GNSS-IR retrieval heights. These experimental results demonstrate that correction with ISEACM can effectively enhance the precision of GNSS-IR sea surface altimetry, which is crucial for accurate sea surface height measurements.

Modeling random isotropic vector fields on the sphere: theory and application to the noise in GNSS station position time series

Modeling random isotropic vector fields on the sphere: theory and application to the noise in GNSS station position time series

Estimating and Analyzing The Stability of Brazilian Gnss Stations for Tectonic Studies

The detailed study of stresses and strains is fundamental for geophysical and geodetic studies, enabling the development of plate motion models, geodynamic processes and several proposals. GNSS positioning has been used to estimate these stresses and strains. However, for the accuracy required, it is necessary a rigorous assessment of the behavior of GNSS stations and their position time series. This paper aims to analyze the influence of the different types of geodetic monuments present in the Brazilian GNSS Continuous Monitoring Network by estimating and characterizing noise in the time series and investigating variables that can affect the stability of the stations to establish criteria for using this data in tectonic and geodynamic studies. The evaluation of 119 GNSS stations throughout Brazil indicates the significant presence of random walk noise in some of these stations, detected by analyzing the magnitude of the noise and the spectral index. It was found that the stations with unstable monuments have been observed for at least eight years and are mostly in coastal areas or under the influence of large watercourses.

Assessment of Satellite Differential Code Biases and Regional Ionospheric Modeling Using Carrier-Smoothed Code of BDS GEO and IGSO Satellites

The geostationary earth orbit (GEO) represents a distinctive geosynchronous orbit situated in the Earth’s equatorial plane, providing an excellent platform for long-term monitoring of ionospheric total electron content (TEC) at a quasi-invariant ionospheric pierce point (IPP). With GEO satellites having limited dual-frequency coverage, the inclined geosynchronous orbit (IGSO) emerges as a valuable resource for ionospheric modeling across a broad range of latitudes. This article evaluates satellite differential code biases (DCB) of BDS high-orbit satellites (GEO and IGSO) and assesses regional ionospheric modeling utilizing data from international GNSS services through a refined polynomial method. Results from a 48-day observation period show a stability of approximately 2.0 ns in BDS satellite DCBs across various frequency signals, correlating with the available GNSS stations and satellites. A comparative analysis between GEO and IGSO satellites in BDS2 and BDS3 reveals no significant systematic bias in satellite DCB estimations. Furthermore, high-orbit BDS satellites exhibit considerable potential for promptly detecting high-resolution fluctuations in vertical TECs compared to conventional geomagnetic activity indicators like Kp or Dst. This research also offers valuable insights into ionospheric responses over mid-latitude regions during the March 2024 geomagnetic storm, utilizing TEC estimates derived from BDS GEO and IGSO satellites.

Constraints on the 2013 Saravan intraslab earthquake using permanent GNSS, InSAR and seismic data

SUMMARY On 16 April 2013, an Mw = 7.7 earthquake struck the border of Iran and Pakistan in the central part of the Makran subduction zone with a reported depth of 80 km by USGS. This rare event in this poorly instrumented region helps to shed light on the kinematics of the subducting slab. We investigate source parameters of the Saravan intraslab normal earthquake using RADARSAT-2 SAR images in three ascending tracks, nine permanent GNSS sites and teleseismic data. The maximum coseismic displacement occurred at the SRVN GNSS station with 54.1 mm southeast horizontal and 42.7 mm upward vertical displacements. The coseismic ascending InSAR displacement maps illustrate a continuous and smooth NE-trending elliptical shape deformation pattern with a maximum of ∼29 cm of displacement away from the satellite. We use 25 broad-band teleseismic P-waveforms to estimate the focal mechanism of the main shock. A joint uniform inversion of InSAR, GNSS and teleseismic data reveals a NW-dipping SW-striking fault and a primarily normal-faulting earthquake with a minor right-lateral strike-slip component. The static slip distribution of the InSAR coseismic maps localizes variable slip at depths between 50 and 81 km with a maximum amplitude >3 m at 60–75.5 km depth, rupturing the oceanic crust of the subducted slab. The kinematic slip distribution exhibits a well-constrained slip pattern with a nucleation depth of 65 km. The source time function indicates that the earthquake reaches its maximum moment tensor release at ∼8 and ∼16 s. The NE-trend of the Saravan earthquake slip pattern, the orientation of the volcanic arc, and the distribution of the intraslab intermediate-depth normal earthquakes provide new insights into slab geometry in the central Makran subduction zone. We suggest that the slab bending at the hinge of subducting Arabian Plate is oblique along a NE–SW direction parallel to the volcanic arc rather than the shoreline or deformation front, and it is likely to be the reason for an oblique volcanic arc in the Makran subduction zone. These new constraints on the Makran slab geometry will help further studies in establishing realistic coupling maps for seismic hazard assessment.

Three-Dimensional Deformation Estimation from Multi-Temporal Real-Scene Models for Landslide Monitoring

This study proposes a three-dimensional (3D) deformation estimation framework based on the integration of shape and texture information for real-scene 3D model matching, effectively addressing the issue of deformation assessment in large-scale geological landslide areas. By extracting and merging the texture and shape features of matched points, correspondences between points in multi-temporal real-scene 3D models are established, resolving the difficulties faced by existing methods in achieving robust and high-precision 3D point matching over landslide areas. To ensure the complete coverage of the geological disaster area while enhancing computational efficiency during deformation estimation, a voxel-based thinning method to generate interest points is proposed. The effectiveness of the proposed method is validated through tests on a dataset from the Lijie north hill geological landslide area in Gansu Province, China. Experimental results demonstrate that the proposed method significantly outperforms existing classic and advanced methods in terms of matching accuracy metrics, and the accuracy of our deformation estimates is close to the actual measurements obtained from GNSS stations, with an average error of only 2.2 cm.

Detecting outliers in local ionospheric model for GNSS precise positioning

Global Navigation Satellite System fast precise positioning can be achieved with accurate ionospheric corrections computed from an adequate number of GNSS stations in a local region. In low-latitude regions, the presence of electron density gradients over short distances can lead to outliers in the map of ionospheric corrections and decrease its accuracy. In this study, we explored outlier detection in ionospheric correction mapping through statistical residuals during a four-month test in 2021. Our findings indicate that the residuals of the local ionospheric model conform to the Laplace distribution. To determine outliers, we use an empirical rule for the Laplace distribution, setting thresholds at μ ± 3b, μ ± 3.5b, and μ ± 5.8b for data retention rates of 95%, 97%, and 99.7%, respectively. Here, μ represents the location parameter, which corresponds to the median of data, and b is the scale parameter, calculated as the medium absolute deviation. We found that while removing outliers can improve model accuracy, it may result in unavailable prediction due to a lack of data in a spare network. For example, applying a μ ± 3.5b threshold for outlier removal led to approximately 2.5% of recording time having no ionospheric corrections map in low-latitude regions, however, the local model has the potential to improve its mean accuracy by up to 50% for both low and mid-latitudes. Therefore, choosing the appropriate percentile threshold depends on the network configuration and the desired accuracy. Removing erroneous satellite data to improve ionospheric accuracy brings positive impacts on precise positioning.

Multi-Instrument Observations of the Ionospheric Response Caused by the 8 April 2024 Total Solar Eclipse

This paper investigates ionospheric response characteristics from multiple perspectives based on globally distributed GNSS data and products, ionosonde data, FORMOSAT-7/COSMIC-2 occultation data, and Swarm satellite observations caused by the total solar eclipse of 8 April 2024 across North and Central America. The results show that both GNSS-derived TEC products have detected the ionospheric TEC degradation triggered by the total solar eclipse, with the maximum degradation exceeding 10 TECU. The TEC data from nine GNSS stations in the path of the maximum eclipse reveal that the intensity of ionospheric TEC degradation is related to the spatial location, with the maximum degradation value of the ionospheric TEC being about 14~23 min behind the moment of the maximum eclipse. Additionally, a negative anomaly of foF2 with a maximum of more than 2.7 MHz is detected by ionosonde. In the eclipse region, NmF2 and hmF2 show trends of decrease and increase, with percentages of variation of 40~70% and 4~16%, respectively. The Ne profile of the Swarm-A satellite is significantly lower than the reference value during the eclipse period, with the maximum negative anomaly value reaching 11.2 × 105 el/cm3, and it failed to show the equatorial ionization anomaly.

Sensitivity of GNSS to vertical land motion over Europe: effects of geophysical loadings and common-mode errors

We perform a statistical sensitivity analysis on a parametric fit to vertical daily displacement time series of 244 European Permanent GNSS stations, with a focus on linear vertical land motion (VLM), i.e., station velocity. We compare two independent corrections to the raw (uncorrected) observed displacements. The first correction is physical and accounts for non-tidal atmospheric, non-tidal oceanic and hydrological loading displacements, while the second approach is an empirical correction for the common-mode errors. For the uncorrected case, we show that combining power-law and white noise stochastic models with autoregressive models yields adequate noise approximations. With this as a realistic baseline, we report improvement rates of about 14% to 24% in station velocity sensitivity, after corrections are applied. We analyze the choice of the stochastic models in detail and outline potential discrepancies between the GNSS-observed displacements and those predicted by the loading models. Furthermore, we apply restricted maximum likelihood estimation (RMLE), to remove low-frequency noise biases, which yields more reliable velocity uncertainty estimates. RMLE reveals that for a number of stations noise is best modeled by a combination of random walk, flicker noise, and white noise. The sensitivity analysis yields minimum detectable VLM parameters (linear velocities, seasonal periodic motions, and offsets), which are of interest for geophysical applications of GNSS, such as tectonic or hydrological studies.