## 4 Methods  
### Primary Geodetic Raw Data Analysis  
The extraction of geodetic raw data, collected by permanently installed GNSS stations, is a cutting-edge, instrumental approach for estimating the upper crust strain rate. These primary geodetic raw datasets are partially related to the general, geodetic velocity pattern, calculated throughout the Greek region [93]. Particularly, in the wider Strymon basin, 16 permanent GNSS stations have been installed, recording the aforementioned geodetic datasets; the denser GNSS network within this area is characterized by a uniform distribution. It is mentioned that these stations are distributed within five different permanent networks (HxGN/SmartNet Greece, Hepos, NoaNet, EUREF Permanent Network--EPN and HermesNet; [94; 95]), providing additional value and reliability to the recorded data.  
Regarding the geodetic velocity calculation, a 7-year (2008-2014) interval was analyzed, while the data recording was performed in a 30 s rate during 24-h operation. In addition, the proposed geodetic velocities reference frame is the European Terrestrial Reference Frame 2000 (ETRF2000), according to which the Eurasian plate is considered fixed. The GAMIT/GLOBK software suite (release 10.5) was implemented [96; 97]. Firstly, the raw data were processed for the entire time period in order to estimate daily loosely constrained solutions of sites, while the orbital and Earth Orientation parameters remained fixed to the IGS final and IERS Bulletin A values with weights, respectively. In our strategy, we estimated Zenith Total Delay (ZTD), adjusted every 2 h interval, and the Vienna Mapping Function1 (VMF1) was selected as a mapping function [96]. Concerning the Solid Earth Tides, the IERS2003 conventions were followed, in addition to Ocean Loading, while we applied FES2004 model, as it was considered appropriate to apply Atmospheric Pressure Loading corrections. In the second step, we combined all the individual sessions of the local sub-networks, using a seven-parameter Helmert transformation, with WRMS 2.0, 2.2, 4.8 mm on the east, north, and vertical components, respectively, estimating site coordinates by using the individual daily solution and geodetic velocities produced by the time series analysis with the linear and extended model. The differences between the two models are \(<\)1 mm/yr for the data span, exceeding the 2.5 years, while the extended model is described by the following equation:  
\[y(t_{i})=\alpha+\beta t_{i}+\sum_{k=1}^{n}(c_{k}\sin(2k\pi t_{i})+d_{k}\cos(2k \pi t_{i}))+\sum_{k=1}^{n_{g}}g_{k}H\Big{(}t_{i}-T_{gk}\Big{)}+v_{i} \tag{1}\]  
where \(\alpha\) and \(\beta\) are the position and velocity of monitoring station, \(c_{k}\) and \(d_{k}\) describe periodic motions, and the equation  
\[z(t_{i})=\sum_{k=1}^{n_{g}}g_{k}H\Big{(}t_{i}-T_{gk}\Big{)} \tag{2}\]  
is used to describe the discontinuities caused by the earthquakes and equipment changes, with magnitudes \(g_{k}\) at epochs \(T_{gk}\) and \(H\) is the Heaviside step function.  
On next step, after the outliers' removal, we estimated station coordinates and velocities using the Kalman filtering sequential approach. The realization of the frame was estimated by minimizing the deviations between positions and velocities, relative to the current terrestrial frame with 7-parameter Helmert transformation and Post RMS 2 mm and 0.6 mm/yr for the coordinates and velocities, respectively. Table 1 shows the station velocities of the processed permanent GNSS sites relative to ETRF2000.