## 2 Materials and Methods
### Methods  
The processing algorithm is shown in the diagram in Figure 2. The subject of analysis was raw daily clock correction data (phase data) obtained from WUM with 30 s sampling intervals and was decimated in the pre-processing phase to a 300 s sampling interval. As an example of data analysis, the processing stages for satellite no. C32 are presented in Figure 3. The analysis was performed on differentiated data (frequency data). Figure 1: Availability of the BeiDou clock product at the WUM analysis centre in 2020. Days with clock data available are marked dark blue. Figure 2: Clock data processing algorithm. Clock metrology uses a form of phase or fractional frequency measurements relative to a reference clock. The phase reading \(\delta t\) is instantaneous (timestamp), while the fractional frequency reading is the difference between the two phases' values divided by their interval, according to [44]:  
\[\frac{d\delta t}{dt}=\frac{\delta t_{i+1}-\delta t_{i}}{t_{i+1}-t_{i}} \tag{1}\]  
The differentiation of the periodic signal changes its amplitude a and shifts phase \(\varphi\) by \(\pi\)/2, but has no impact on signal component frequency \(f\). The amplitudes are scaled by \(2\pi\)_f_, according to the following equations:  
\[\delta t(t)=\sum a_{i}sin(2\pi f_{i}+\varphi_{i}) \tag{2}\]  
\[\frac{d\delta t}{dt}=\sum 2\pi f_{i}a_{i}sin\left(2\pi f_{i}+\varphi_{i}+\frac{ \pi}{2}\right) \tag{3}\]  
The first step of our analysis was to differentiate clock correction data \(\delta t\) (Figure 3a) to produce a time series of clock correction rate \(d\delta t\)/\(dt\) (Figure 3b). The outliers present in the time series were removed two-fold. First, values of \(d\delta t\)/\(dt\) larger than \(10^{-10}\) s/s (arbitrarily chosen) were treated as outliers and were removed from the data set (Figure 3b). The threshold value was chosen based on the analysis of the whole data set and is adopted as \(10^{-10}\) s/s for proper detection of the jumps in the signal (phase data). The jumps that existed in the data series were mainly caused by the day-boundary effect [18]. The maximum number of removed data points at this stage was 0.8%. Such a small amount of removed data has no impact on further data processing. After removal of the primary outliers, jumps in clock correction rate data remained identifiable for several satellites, as is visible in Figure 3C. These discontinuities were removed (Figure 3d) from the data. Next, the third order trend (polynomial) and secondary outlier detection and removal were performed. At this stage, outliers were detected with the aid of the _MAD_ (median absolute deviation) method commonly used in clock stability analysis [44] (Figure 3e). The _MAD_ coefficient is calculated according to:  
\[MAD=med\left(\frac{|x_{i}-med(x)|}{0.6745}\right) \tag{4}\]  
In Equation (4), \(x_{i}\) is the frequency data and \(med(x)\) is a median of \(x\), and then outliers in the time series are identified if:  
\[\frac{|x_{i}-med(x)|}{MAD}>k \tag{5}\]  
In our calculations, the typical value for constant \(k=5\) was adopted [44], and after this stage, the maximum removed data point was 2.2%. Lastly, the Lomb-Scargle [45; 46] power spectrum was calculated (Figure 3f). Herein, we used the Lomb-Scargle method to obtain an approximation of the power spectrum of the analysed signal when gaps in the data occurred due to the removal of outliers or missing observations from the data source. In the resulting L-S spectrum, only period values for which the probability of false identification lower than P (spectrum peak \(\mid\) noise) = 0.00005 were selected as significant. A relatively low value of P was chosen in order to filter out spurious periods present in the analysed clock data. In Figure 3f, meaningful periods are those for which the peaks are above the orange line, i.e., 6.4 h and four peaks indicating periods over 30 h. The relevant periods for each analysed BeiDou satellite are presented as a blue mark in Appendix B.