## 2 Dataset and Methodology  
This study focuses on Attica and its neighboring areas in Central Greece. Several permanent Hellenic Unified Seismological Network (HUSN) stations operate in this area, complemented by 6 temporary stations of the Geodynamics Institute of the National Observatory of Athens (GI-NOA), around the 19 July 2019 aftershock zone. The earthquake catalog consists of 3468 earthquakes (1 January 2019-25 March 2024), with local magnitude (M\({}_{\rm L}\)) ranging between 0.6 and 5.1, recorded by the broad-band and strong motion stations of the Hellenic Unified Seismological Network (HUSN; [[http://eida.gein.noa.gr/](http://eida.gein.noa.gr/)]([http://eida.gein.noa.gr/](http://eida.gein.noa.gr/)) (accessed on 31 March 2024) [30; 31; 32; 33; 34; 35]). The hypocentral parameters were determined using a local 1-D velocity model proposed by [24] (Figure 3).  
The data selection was performed under the following criteria:  
1. The nearest station epicentral distance should be less than 50 km.
2. The number of the P- and S-arrival picks per event should be equal to or larger than 10.  
Figure 3: Location of the initial 3468 events (M\({}_{\rm L}\geq\) 0.6) from January 2019 to March 2024. HUSN stations are depicted in blue triangles. Fault traces as in Figure 1.  
3 Temporal residual should not be larger than 0.4 s.  
The application of these criteria led to a slightly reduced dataset of 3397 seismic events, including 42,122 P and 32,738 S picks (Figure S2).  
In the framework of this study, the analysis was performed using the LOcal TOmographic Software (LOTOS-12) by [36]. The program provides us with two options: either the V\({}_{\rm P}\)-V\({}_{\rm S}\) scheme, using P and S travel-time residuals (dt\({}_{\rm P}\) and dt\({}_{\rm S}\)), or the V\({}_{\rm P}\)-V\({}_{\rm P}\)/V\({}_{\rm S}\) scheme, selecting the available dt\({}_{\rm P}\) and differential residuals, dt\({}_{\rm S}-\)dt\({}_{\rm P}\). During the tomographic inversion, we used both schemes to obtain more information as they matter for the V\({}_{\rm P}\) and V\({}_{\rm S}\) anomalies [36; 37]. Body-wave travel times and ray paths (Figures S2-S6) were calculated based on the bending method as proposed by [38]. The calculations are extensively described in supportive material (Section 2).  
The procedure was performed simultaneously for the body-wave velocities, as well as for the source parameters (dx, dy, dz, and dt) and the station corrections. The LSQR algorithm was used for the inversion [39]. The stability of the procedure was controlled by damping and flattening. To minimize the differences in the solutions from neighboring nodes, we employed a technique to flatten the velocity anomalies. This was achieved by adding equations with two non-zero elements having opposite signs, corresponding to all pairs of neighboring nodes. The best values for the flattening parameters were determined through synthetic modeling for the most accurate model recovery. It is important to note that the number of iterations has a similar effect to damping changes. Therefore, we consistently applied five iterations and only adjusted the flattening parameters. In total, we used five iterations to derive the final models.  
The resolution of the available dataset was evaluated using the checkerboard synthetic test [40] and it is described in the supportive information (Section 3; Figures S7-S12). This method involves applying alternating anomalies of positive and negative velocity perturbations on a 1-D gradient model that is evenly distributed in a checkerboard pattern throughout the model. The optimal values for P- and S-wave amplitude damping and smoothing, as well as the corrections for station and source coordinates and origin times, were determined based on the results of these tests, as described in [36; 41; 42].