## 2 Materials and Methods
### Shear-Wave Splitting (SWS) Analysis
Figure 3: Example of the automatic determination of the SWS parameters with the Pytheas software, using the Eigenvalue (EV) and Cluster Analysis (CA) methods. Horizontal recordings (**a**) in their initial N-S/E-W co-ordinate system (N-E), after the application of band-pass Butterworth filter between 0.5 and 5.0 Hz, (**b**) rotated into the S\({}_{\text{fast}}\)/S\({}_{\text{slow}}\) co-ordinate system (F-S), having measured \(\varphi\) = N94\({}^{\circ}\)E, then shifted S\({}_{\text{slow}}\) leftwards by \(t_{\text{d}}\) = 60 ms and (**c**) rotated back to the N-S co-ordinate system, after being corrected for anisotropy, where the \(p\) = N58.8\({}^{\circ}\)E polarization direction is measured. The S-wave particle motion in the N-S co-ordinate system is also presented in the hodograms of panels (**d**) for the initial waveforms and (**e**) after correction for anisotropy. (**f**) Contours of the distribution of \(\lambda_{2}\) values during the grid-search of (\(\varphi\), \(t_{\text{d}}\)) pairs performed for the application of the EV method, with the polarization direction displayed in the form of \(\varphi\)-_baz_, i.e., degrees relative to the back-azimuth, wrapped in a range of \(\pm 90^{\circ}\) with respect to the R-axis direction. The selected (\(\varphi\), \(t_{\text{d}}\))  
pair, corresponding to the global minimum of \(\lambda_{2}\), is marked by the dashed cross-hair, pointing to the middle of a dark-shaded region that indicates the 95% confidence interval for \(\lambda_{2}\)[42], from which \(\delta\varphi\) = 0.3\({}^{\circ}\) and \(\delta t_{\rm d}\) = 0.03 ms are measured. (**g-k**) Diagrams of the Cluster Analysis (CA) method for the determination of the best window presenting: (g) the (\(\varphi\), \(t_{\rm d}\)) measurement pairs, (**h**) the performed clustering, (**i**) the selection of the most constrained cluster and measurement with the smallest errors, (**j**) the series of \(\varphi\) and (**k**) the series of \(t_{\rm d}\) measurements per trial window index, with their respective uncertainties. Figure 3 shows an example of the automatic analysis with the combination of the EV and CA methods for an event recorded at station PSAR. Initially, the automatic filtering algorithm decided the optimal band-pass Butterworth filter range between 0.5 and 5.0 Hz. Then, SWS measurements were performed in 102 trial signal windows. For each window, the EV method was applied to determine the \(\varphi\) and \(t_{\rm d}\) values. The horizontal components were first rotated into the Radial-Transverse (R-T). Then, grid-search was performed, rotating the waveforms to the various \(\varphi\) values between R -90\({}^{\circ}\)T and R90\({}^{\circ}\)T, and shifting the transverse component leftwards by a time-delay \(t_{\rm d}\) between 0.0 and 250.0 ms. In this example, the solution pair that yielded the minimum \(\lambda_{2}\) value was \(\varphi\) = N94.0\({}^{\circ}\) E \(\pm\) 0.3\({}^{\circ}\) and \(t_{\rm d}\) = (60.00 \(\pm\) 0.03) ms in a 95% c.i. with grade \"A\", indicating a measurement of the best quality. Figure 3g presents the selection of the optimal solution. It is noteworthy that despite the \(\varphi\) values being spread in a range roughly between 80\({}^{\circ}\) and 120\({}^{\circ}\), the majority are clustered near the selected value of 94\({}^{\circ}\), with the respective \(t_{\rm d}\) at 60 ms (the preferred value). The diagram of Figure 3f shows the global minimum of \(\lambda_{2}\) for the selected signal window, marked by the dashed cross-hair. The validity of the \(\varphi\), \(t_{\rm d}\) measurements is verified in Figure 3e, where the resulting particle motion becomes nearly linear after the removal of anisotropy. We verified the instrument orientation during the study period (2013-2014) using the method of [47], to ensure the validity of the polarization direction measurements.