## 4 Seismological Methods  
We acquired initial event and arrival information from the routine analysis catalog of NKUA-SL, between September 2020 and May 2021. This time period was then narrowed down, isolating the earthquake solutions that--in temporal terms--belonged to the Thiva sequence. After exploring the daily distribution of earthquakes, we concluded that the whole sequence occurred between 3 December 2020 and 3 January 2021. Phase arrivals for these earthquakes were picked anew, to ensure an as detailed as possible catalogue. Recorded waveform data were acquired from the European Integrated Data Archive (EIDA) node at GI-NOA [29], using the integrated Federation of Digital Seismograph Networks (FDSN) data selection service. All available stations across Greece were used, as long as the P and/or S arrivals were identifiable (see Appendix B for more details). We then re-estimated the hypocentral locations of the earthquakes of the sequence. To this end, we used the Hypoinverse location code [30] and the velocity model proposed by [31] for the Eastern Gulf of Corinth and Boeotia. We note the severe limitations in depth estimation imposed by the absence of at least one seismological station near (less than 15.0 km) the epicentral area. Thus, this was constrictive towards any speculation about the vertical distribution of seismicity. Furthermore, stations are ill-distributed azimuthally, as most (in intermediate distances) are located SE of the sequence (Figure 1). Such issues have been present in other areas with either sparsely located or no seismological instruments [32]. It is worth noting that the re-picked catalogue yielded satisfactory results. The final set of 64 events (Figure 5) features a shallow average focal depth of 3.5 km and mean horizontal and vertical errors of 0.7 km and 1.5 km, respectively. Magnitudes (M\({}_{\rm L}\)) were slightly affected by the re-picking process and range between 0.8 and 3.3, the latter for the largest aftershock that occurred approximately two hours after the mainshock. The magnitudes of the following events were weaker than 2.3. To determine the focal mechanism and seismic moment of the main event, we used the latest available version of the ISOLA software [33; 34]. ISOLA is widely established as a reliable program of estimating focal mechanisms of varying magnitudes at local and regional distances [35; 36; 37]. The code searches for the details, through moment tensor inversion, that best fits an initial point source and then generates synthetic seismograms based on the solution, which are compared with the observed ones. This process is repeated for multiple trial sources. The best-fitting solution is determined using a variety of metrics, including similarity between observed and synthetic waveforms [38]. In this process, we used the same velocity model as before [31]. Moreover, waveforms were preprocessed with ObsPy [39] to remove the instrument's response and were converted to velocity. Stations were selected with the best possible azimuthal distribution in mind (see also Appendix B and Figure A6), after visually inspecting recordings for noise or other issues. Multiple band-pass filters were tried to obtain the bandwidth that best removes noise for the frequency content used in the inversion. The optimal filter corners were 0.05 and 0.10 Hz. In all trial filters, the bottom limit imposed by the available instruments' characteristics was considered. We performed a search for 72 sources on a SE-dipping plane, as determined by the initial GI-NOA solution (strike of N111\({}^{\circ}\)E and dip of 35\({}^{\circ}\)). This nodal plane of the initial solution was determined by the results of geodesy, which supports a SW-dipping fault. The source grid started at the mainshock's focal depth, as determined by phase arrivals (0.9 km), and extended down to 7.4 km. The assumed fault was 6 km long (along-strike dimension) and 11 km-wide (along-dip). Due to the great data quality, all components of selected stations were used in the inversion.