## 4 Data and Methods  
### GNSS Data  
In this study, original continuous multi-GNSS measurements of high rate (1 s) were downloaded and processed. The original measurements were downloaded in RINEX (_eng. Receiver Independent Exchange_) 3.03 format from the CORS CROPOS stations for the purpose of post-processing. The CORS stations provide continuous measurements in one-second intervals. However, due to substantial amount of data, authors found it appropriate to use the 30-min interval of the measurements and, therefore, origin data were resampled. For each selected CORS, daily multi-GNSS data files were generated (e.g., 1 January, 00:00:00-23:59:30) with sampling rate of 30 s (one observation per 30 s).  
For each station, data of 10 days before and 10 days after the earthquake were processed. Median value of 10 days prior to earthquake for each of three coordinates (E, N, h) was used as a reference coordinate value for further analysis. Regarding the reference value, the shift in northing, easting, and height direction was computed.  
Displacements (shift) of stations were calculated as the differences between values after and before the earthquake:  
\[Delta\,E =E_{after}-E_{before}, \tag{1a}\] \[Delta\,N =N_{after}-N_{before},\] (1b) \[delta\,h =h_{after}-h_{before}, \tag{1c}\]  
with horizontal shift  
\[displacement_{horizontal}=\sqrt{E_{shift}^{2}+N_{shift}^{2}}, \tag{2}\]  
Each daily GNSS Rinex file consisted of 2880 epochs (24 h files with 30 s interval) and had double-frequency (L1, L2) multi-constellation data collected with Trimble NETR9 and TRM115000.00 Trimble Zephvr 3 Geodetic antennas. Antenna parameters were APC to ARP L1 = 0.064 m, L2 = 0.057 m, whereas ARP to Marker was E = N = h = 0.000 m. The results of processing were X, Y, Z coordinates for each day for each CORS in ITRF2014 (e. 2021.0). The data were processed in Precise Point Positioning (PPP) static mode using both phase and code observations. Only GPS and GLONASS data were used. IGS final orbits and clocks were used. Elevation mask was set to 10\({}^{\circ}\). Ionospheric free linear combination was used to eliminated ionospheric errors. The troposphere correction was modeled using the Zenithal Wet Delay (ZWD) component as a random walk process. Other corrections were applied as well, such as: Ocean Tide Loading (OTL), Earth Body Tide (EBT), Phase Center Variation (PCV) of receiver and satellites, and phase windup. Due to the high quality of measurements, ambiguity was usually fixed after 1 or 1.5 h, with a posteriori accuracy estimation ranging from 2 mm to 5 mm for E and N, and from 5 mm to 1 cm for height component.  
The result of multi-GNSS data processing was three-dimensional position of the CORS, or, more precisely, geocentric Cartesian X, Y, Z coordinates. Resulting coordinates were converted to Easting (E), Northing (N), and height (h) for the purpose of all analyses. For Easting and Northing coordinates, UTM North (zone 33) projection was used. These values were then used as the basis for determination of the space-time sequences of horizontal shift of the Earth's surface on the affected area.  
All observed CROPOS GNSS reference stations were computed with the displacement shift using the adjusted daily GNSS values. To compute the displacement of the station's position for each day before and after the earthquake, computed daily solutions were subtracted from the median value of the time series before and after earthquake values, respectively.