## 3 Description of the Time Series Statistics Used  
For each station in each time window, four surface tremor statistics are computed for the three ground displacement components. Spectral exponent. Let \(x(t),\ t=1,\ldots,N\) be a sample of the GPS time series within a window of sample lengths \(N\). Denote by \(S_{xx}(\omega_{j})\) the estimate of the signal \(x(t)\) power spectrum \(\omega_{j}\)--a sequence of discrete values of frequencies with a step \(\Delta\omega\), \(\Delta\omega=2\pi/(N_{f}\cdot\Delta s)\), \(j=1,\ldots,N_{f}\), \(\Delta\)s--a time step. Here, \(N_{f}\) is the minimum length, which is greater than or equal to the sample length \(N\): \(N_{f}=2^{n}\geq N\), the index \(j\) in Formula (6) varies up to \(N_{f}/2\) due to the symmetry of the power spectrum of real signals with respect to the Nyquist frequency \(\omega_{N_{f}/2+1}\): \(S_{xx}(\omega_{N_{f}/2+1-k})=S_{xx}(\omega_{N_{f}/2+1+k})\), \(k=1,\ldots,N_{f}/2-1\). We introduce the spectral exponent \(\alpha\) by the following formula:  
\[\log(S_{xx}(T_{j}))=\alpha\cdot\log(T_{j}))+d+\eta_{j},\ j=1,\ldots,N_{f}/2, \tag{1}\]  
where period \(T_{j}=2\pi/\omega_{j}\), the value \(\alpha\) is determined by the least squares method: \(\sum_{j=1}^{N_{f}/2}\eta_{j}^{2}\rightarrow\min\limits_{d,\alpha}\), \(\eta_{j}\) is a sequence of random variables with zero mean. In each time window, the power spectrum was computed after transition to increments using an autoregressive model (maximum entropy estimate) [24] with an autoregressive order of \(N/10\), i.e., 36 for a window of length 365 samples. The transition to increments before calculating the power spectrum aims to free from the dominant effect of low frequencies in the daily time series of GPS. We then address wavelet-based entropy and wavelet-based spectral slope. The minimum normalized wavelet information entropy for the sample \(x(t),\ t=1,..,N\) is determined by the following formula:  
\[En=-\sum_{k=1}^{N}p_{k}\cdot\log(p_{k})/\log(N),\ \ \ p_{k}=c_{k}^{2}/\sum_{j=1}^{N}c_{j}^ {2} \tag{2}\]  
Here, \(c_{k}\) denotes the coefficients of the orthogonal wavelet decomposition. Daubechies bases were used with the number of vanishing moments from \(1\) to \(10\)[25]. We chose such a basis from this set for which the value (2) is minimal for a given time window and for the considered component of the displacement of the Earth's surface. By Formula (2), \(0\leq En\leq 1\). Figure 2: The number of operational GPS stations in a window of the length 365 days with mutual shift of 7 days. For the orthogonal wavelet, which is found from a minimum of entropy, the wavelet power spectrum can be calculated as the average values of the squares of the wavelet coefficients at each level of detail of the decomposition:  
\[S_{k}=\sum_{j=1}^{M_{k}}\bigl{|}c_{j}^{(k)}\bigr{|}^{2}/M_{k}. \tag{3}\]  
In Formula (3), \(c_{j}^{(k)}\) are the wavelet coefficients at the detail level with number \(k\), \(k=1,\ldots k_{\max}\), index \(j=1,\ldots M_{k}\), where \(M_{k}\) is the number of wavelet coefficients at the detail level with number \(k\), which corresponds to the frequency band \([\frac{1}{2^{(k+1)}\Delta s},\ \frac{1}{2^{k}\Delta s}]\) with the central period \(T_{k}=\frac{2\Delta s}{2^{-k}+2^{-(k+1)}}\)[25]. The maximum detail level number \(k_{\max}\) is chosen such that it contains at least a given number \(L_{\min}\) of wavelet coefficients. Recall that the number of wavelet coefficients at the level of detail decreases by a factor of 2 when the level number increases by one. The value \(L_{\min}=4\) was used in the calculations, which gives the value \(k_{\max}=6\) for the window length \(N=365\). Values \(S_{k}=S(T_{k})\), are similar to conventional power spectra. The regression model is used to calculate the wavelet-based spectral slope:  
\[\log(S(T_{k}))=\beta\cdot\log(T_{k}))+c+\varepsilon_{k},\ k=1,\ldots,k_{\max}, \tag{4}\]  
where \(\varepsilon_{k}\) is a sequence of random variables with zero mean.