Combined Geodetic and Seismological Study of the December 2020 M w = 4.6 Thiva (Central Greece) Shallow Earthquake

: On 2 December 2020, a moderate and shallow M w = 4.6 earthquake occurred in Boeotia (Central Greece) near the city of Thiva. Despite its magnitude, the co-seismic ground deformation ﬁeld was detectable and measurable by Sentinel-1, ascending and descending, synthetic aperture interferometry radar (InSAR) acquisitions. The closest available GNSS station to the epicenter, located 11 km west, measured no deformation, as expected. We proceeded to the inversion of the deformation source. Moreover, we reassessed seismological data to identify the activated zone, associated with the mainshock and the aftershock sequence. Additionally, we used the rupture plane information from InSAR to better determine the focal mechanism and the centroid location of the mainshock. We observed that the mainshock occurred at a shallower depth and the rupture then expanded downdip, as revealed by the aftershock distribution. Our geodetic inversion modelling indicated the activation of a normal fault with a small left-lateral component, length of 2.0 km, width of 1.7 km, average slip of 0.2 m, a low dip angle of 33 ◦ , and a SW dip-direction. The inferred fault top was buried at a depth of ~0.5 km, rooted at a depth of ~1.4 km, with its geodetic centroid buried at 1.0 km. It was aligned with the Kallithea fault. In addition, the dip-up projection of the modeled fault to the surface was located very close (~0.4 km SW) to the mapped (by existing geological observations) trace of the Kallithea fault. The ruptured area was settled in a transition zone. We suggest the installation of at least one GNSS and seismological station near Kallithea; as the activated zone (inferred by the aftershock sequence and InSAR results) could yield events with M ≥ 5.0, according to empirical laws relating to rupture zone dimensions and earthquake magnitude.


Introduction
Boeotia is located in Central Greece, an area that has had a mild seismic footprint during the 20th and 21st centuries [1]. Its largest city, Thiva, was founded in antiquity. There are ancient reports of various incidents near the city that could be interpreted as geotectonic [2]. Archaeological evidence suggests earthquake activity, which had an impact on Thiva, between 1350 and 1230 BC [3]. Historical catalogs refer to destructive events since 1321 AD, up to 1914, with estimated magnitudes ranging between 6.0 and 6.5, causing extensive damage and collapses [4][5][6]. It has been suggested by [6] that the causative fault of the 1914 event was an E-W striking structure located between the Kallithea and Asopia villages ( Figure 1). Thus, the city of Thiva is located close to active faults that have not caused a strong earthquake (M ≥ 6.0) in over a century (Figure 1).

Figure 1.
Map of the broader Central Greece area. Epicenters for earthquakes with M ≥ 6.0 (stars) are shown for the whole area, for the historical (white) and instrumental (red) eras. Weaker events, recorded since 1900, are presented (circles) for the study area (brown rectangle). Seismological stations of HUSN (triangles) are also shown. The blue line indicates the boundary between the Central Greece (CGb) and Central Aegean (Cab) blocks proposed by [7]. Faults (black lines) from the NOAFaults v3.0 database [8,9]. Further details about the faults of the study area in Figure 5. Historical earthquakes from [4]. Instrumental seismicity between 1900-2009 according to [1]. Seismicity between 2010 and 2019 was obtained from the reviewed bulletin of the International Seismological Centre. The cities of Athens and Thiva are also shown, as well as the Asopia and Kallithea villages (squares). Inset: location of the map region (black rectangle) in Greece. The trace of the subduction zone is also shown (thin black front). Double arrow indicates the relative movement of the two blocks' barycenter.
Thiva, and its surrounding area, is located in the transition zone between the Central Greece (CGb) and Central Aegean (Cab) blocks, as inferred by [7], with the latter's barycenter moving 3.6 mm·yr −1 eastwards and 9.1 mm·yr −1 southwards, relative to the former one, thus 9.7 mm·yr −1 at an azimuth of N202°Ε. The passage from the structures of the Aegean to the NNE-SSW rifting zone in the Gulf of Corinth [10][11][12] is marked by smaller faults systems, according to the developed stress field [13][14][15].
The area to the east of Thiva is dominated by the WNW-ESE segmented faults [8,9], organized along-strike, as seen in Figure 1. The broader area features N-S slip vectors [13,16,17]. These two characteristics insinuate a clockwise rotation, related to a vertical axis, that affects the fault systems. Tectonic structures construct a series of parallel limestone ridges, bounded by south-and north-dipping faults, on each side [18]. One of these segments is the Kallithea fault, with an inferred length of 6.0 km. Its strike and dip have been observed in escarpments; [18] reported exposed fault surfaces, striking between N80°E and N115°E, with striation vectors in the range of N180°E to N200°E. The authors of [19] carried out field mapping and morphotectonic analysis and estimated a combined length between the faults of Kallithea and Mavrovouni equal to 10.8 km. They proposed that activation of these two segments in tandem (as they comprise a zone) would yield a large earthquake with 6.2 ≤ M ≤ 6.4.
Over the last thirty years, and especially the more recent ones, inversion of synthetic aperture radar interferometry (InSAR) and global navigation satellite system (GNSS) measurements have been extensively applied to study earthquakes in the Euro- Figure 1. Map of the broader Central Greece area. Epicenters for earthquakes with M ≥ 6.0 (stars) are shown for the whole area, for the historical (white) and instrumental (red) eras. Weaker events, recorded since 1900, are presented (circles) for the study area (brown rectangle). Seismological stations of HUSN (triangles) are also shown. The blue line indicates the boundary between the Central Greece (CGb) and Central Aegean (Cab) blocks proposed by [7]. Faults (black lines) from the NOAFaults v3.0 database [8,9]. Further details about the faults of the study area in Figure 5. Historical earthquakes from [4]. Instrumental seismicity between 1900-2009 according to [1]. Seismicity between 2010 and 2019 was obtained from the reviewed bulletin of the International Seismological Centre. The cities of Athens and Thiva are also shown, as well as the Asopia and Kallithea villages (squares). Inset: location of the map region (black rectangle) in Greece. The trace of the subduction zone is also shown (thin black front). Double arrow indicates the relative movement of the two blocks' barycenter.
Thiva, and its surrounding area, is located in the transition zone between the Central Greece (CGb) and Central Aegean (Cab) blocks, as inferred by [7], with the latter's barycenter moving 3.6 mm·yr −1 eastwards and 9.1 mm·yr −1 southwards, relative to the former one, thus 9.7 mm·yr −1 at an azimuth of N202 • E. The passage from the structures of the Aegean to the NNE-SSW rifting zone in the Gulf of Corinth [10][11][12] is marked by smaller faults systems, according to the developed stress field [13][14][15].
The area to the east of Thiva is dominated by the WNW-ESE segmented faults [8,9], organized along-strike, as seen in Figure 1. The broader area features N-S slip vectors [13,16,17]. These two characteristics insinuate a clockwise rotation, related to a vertical axis, that affects the fault systems. Tectonic structures construct a series of parallel limestone ridges, bounded by south-and north-dipping faults, on each side [18]. One of these segments is the Kallithea fault, with an inferred length of 6.0 km. Its strike and dip have been observed in escarpments; [18] reported exposed fault surfaces, striking between N80 • E and N115 • E, with striation vectors in the range of N180 • E to N200 • E. The authors of [19] carried out field mapping and morphotectonic analysis and estimated a combined length between the faults of Kallithea and Mavrovouni equal to 10.8 km. They proposed that activation of these two segments in tandem (as they comprise a zone) would yield a large earthquake with 6.2 ≤ M ≤ 6.4.
Over the last thirty years, and especially the more recent ones, inversion of synthetic aperture radar interferometry (InSAR) and global navigation satellite system (GNSS) mea-surements have been extensively applied to study earthquakes in the Euro-Mediterranean region, and calculate the deformation source parameters [7,[20][21][22][23][24][25]. The GNSS technique can provide accurate information of the ground motion on all three components, but it is limited to point-wise coverage of the study area. Conventional differential SAR interferometry describes the relative spatial deformation of the area (with respect to a reference point) and is limited to the line-of-sight (LOS) direction. GNSS and InSAR can be used jointly on a seismically activated area, to inverse for the parameters of the deformation source and determine the absolute ground deformation field, as well as calculate the E-W and vertical deformation components.
On 2 December 2020 10:54:56 UTC, a M L = 4.5 earthquake occurred east of Thiva. According to the Seismological Laboratory of the National and Kapodistrian University of Athens (NKUA-SL), as initially determined, it was located at 38.3140 • N, 23.4663 • E, characterized by a normal focal mechanism with N137 • E, 56 • and −37 • for strike, dip and rake, respectively. According to the Geodynamics Institute of the National Observatory of Athens (GI-NOA), the revised location was at 38.3226 • N, 23.4489 • E and the focal mechanism was dip-slip normal, with N111 • E, 35 • and −78 • for strike, dip and rake, respectively (see Figure A5 in Appendix B for more details about existing focal mechanism solutions). The aftershocks lasted until 3 January 2021.
Even with a moderate magnitude, the deformation field was measurable by synthetic aperture radar interferometry (InSAR), exploiting Copernicus Seninel-1 data in both orbital directions (i.e., ascending and descending). The closest permanent GNSS station from which we were able to acquire data was THIV, located in the city of Thiva. We invert the geodetic data to estimate the rupture plane with uniform slip. We also reassess seismic waveform data to locate all events in order to identify the activated zone, associated with the mainshock and the aftershock sequence. Moreover, we use the rupture plane information from InSAR to better determine the focal mechanism and the mainshock centroid.

Geodetic Methods
GNSS data from a continuous station (THIV), which is located in the city of Thiva (at 38.3166 • N, 23.3183 • E), about 11 km west of the epicenter, and belongs to the METRICA S.A. commercial network (HexagonSmartNet) were processed (Appendix A). The formed time-series for THIV station ( Figure A1) showed that the moderate mainshock did not caused any observable ground deformations in Thiva. The fluctuation of the two horizontal components (north and east) prior to and after the earthquake remained on the same limits, of about ±3 mm from an average value, whereas the vertical component varied at ±10 mm. There was no evident static co-seismic displacement on the station after the 2 December event, or any alteration on the motion pattern for the period following the event.
We used synthetic aperture radar (SAR) C-band acquisitions in the IW mode, acquired over the study area by the Copernicus European satellite Sentinel-1 from the ascending tracks 102 and the descending track 7. We calculated with ENVI ® SARscape ® software version 5.5 the interferograms between 21 November 2020 and 3 December 2020 in both orbital directions (Figure 2, Figure 6, and Figure A3) and having a LOS unit vector of [−0.61, −0.14, 0.78] and [0.59, −0.14, 0.80] for the ascending and descending geometry, respectively. Precise orbits and the 30-m NASA shuttle radar topography mission (SRTM) 1-s digital elevation model (DEM) were used [26]. The interferograms that we processed had the most coherent deformation patterns, with slightly more than two fringes in the ascending and almost two in the descending geometry. The deformation lobes and the similarity for both orbits declare typical subsidence with a preference to normal mechanism. In addition, the increased density at the NE suggests that the deformation source is located on this side.
As InSAR provides information on only two independent LOS measurements ( Figure 2), in order to retrieve the E-W and vertical components (but not for the N-S one, due to the small sensitivity along that direction, caused by its close alignment with the orbit inclination), the interferograms were unwrapped. Furthermore, their planar trends were removed and their mean value (excluding the deformed area from both tasks) was set to zero. Finally, by using the equation provided in [27], we obtained their components (Figure 3a,b). As InSAR provides information on only two independent LOS measurements ( Figure 2), in order to retrieve the E-W and vertical components (but not for the N-S one, due to the small sensitivity along that direction, caused by its close alignment with the orbit inclination), the interferograms were unwrapped. Furthermore, their planar trends were removed and their mean value (excluding the deformed area from both tasks) was set to zero. Finally, by using the equation provided in [27], we obtained their components (Figure 3a   As InSAR provides information on only two independent LOS measurements ( Figure 2), in order to retrieve the E-W and vertical components (but not for the N-S one, due to the small sensitivity along that direction, caused by its close alignment with the orbit inclination), the interferograms were unwrapped. Furthermore, their planar trends were removed and their mean value (excluding the deformed area from both tasks) was set to zero. Finally, by using the equation provided in [27], we obtained their components (Figure 3a

Inversion of Geodetic Data-Fault Model
As an immediate observation, the distribution, the shape and direction of the fringes, of both orbital geometries, supported a dominant normal mechanism with a fault located in the N-E of the pattern, dipping SW, with a strike parallel to the elongated axis of the fringes and a length comparable with the extension along the low curvature fringes on the N-E side. Accordingly, we inverted the InSAR geodetic data, assuming a half-space elastic model with uniform slip along a rectangular fault surface, using the modified code inverse6 [28]. The same weight was assigned to both the ascending and descending picks. We sampled the ascending track 102 at 34 points (fringe pickings) and the descending track 7 at 29 points. The 63 total data points were compared against the model in Appendix A, Table A1.
The geodetic observations presented in the previous section clearly and unambiguously imply that the seismic fault plane is the one dipping towards the south-west. Indeed, there is no solution with the antithetic plane that can fit correctly the InSAR data. Moreover, the south-dipping plane is consistent with the geological observations in [18], that reported exposed surfaces related to Kallithea fault, very close to the geodetically inferred one, and the tectonic context of the area. Finally, considering the higher normal component and the denser fringes to the NE of the deformed area, we suggest a fault lying beneath, with its top located towards the dense fringes area and, thus, dipping SW.
Since the pattern of the deformation was consistent, well-shaped, with a high coherence, small contribution of tropospheric disturbances and with full spatial coverage, we could invert all the parameters, except dip. We tested a control group with a range of fixed dip angles (35-55 • , at 5 • intervals) and a number of 3553 inversions. The convergence diagrams declare high confidence of the coordinates, the depth (of the upper edge of the fault), the length and the strike (see Appendix A, Figure A2b-f). The latter was calculated as N120 • E and could also be extracted and confirmed by the inclination of the solution field of the diagram of Figure A2a, as well as from [18]. The test diagram on different dip values ( Figure A2g) reveals an overall minimum at~48 • . Existing focal mechanism solutions from seismology provided a dip angle at 33 • . Thus, we preceded with two groups of inversions. For the first group, we locked the azimuth at N120 • E and dip at 48 • for 33 inversions, ending up with a dip angle favored by geodesy, provided in Table 1. For the second group we locked the azimuth at N120 • E and the dip at 33 • for 33 iterations, to acquire the dip angle favored by seismology in Table 1. The uncertainties were estimated from the inversion of the control group ( Figure A2). Table 1. Parameters of the two geodetic models, for the dip angle favoring geodesy and for the dip angle favoring seismology (bold line). The standard deviation of the measured and modeled deformation values (plotting in Figure 4) are mentioned. The uncertainties of the parameters are also included. earthquakes and fault orientations, as was the case of the Samos (Greece) 2020 Earthquake.

Parameter Favored Geodesy Dip Angle Favored Seismology Dip Angle Uncertainties
In a recent publication [20] on an effort to decompose the LOS values, from mean interferograms, using both orbital geometries, the modeled N-S component was considered in the equation. In Figure 4, the modeled and picked deformation values are shown for the ascending and descending orbits, where the misfit of the model and the measurements could be controlled. In terms of earthquake size, the geodetic moment was calculated to 1.6 × 10 16 N·m and 1.9 × 10 16 N·m, for the fault parameters favoring the geodetic and seismological dip angle, respectively, assuming a rigidity of the medium of 3.0 × 10 10 Pa. With the dip angle of 48°, the top-edge of the fault at 0.7 km, and a width of 1.0 km, the root of the fault was located at a depth of 1.5 km and the geodetic centroid depth was at 1.1 km. Regarding the model with the dip angle of 33°, the top-edge of the fault at 0.5 km and a width of 1.7 km, The solution with the geodetically constrained dip angle has a lower standard deviation of 2.6 mm, compared to the 3.9 mm for the one constrained by the seismological dip. The two solutions are similar, mainly differing, apart from the dip angle, at the rake angle and the width. The distance between the two centroids is short,~0.4 km.
As in all elastic models, it was not possible to solve independently for all parameters of the modeled fault. This is the reason why we had to lock the dip and azimuth angles. To test the trade-offs between related parameters during the inversion, the convergence diagrams of the width versus dip angle, rake angle, and slip are considered (see also Figure A2 in Appendix A). After constraining the strike and dip, these trade-offs ceased. Although the azimuth could be sharply estimated, the optimum dip angle is hardly distinguishable ( Figure A2g) compared with the other parameters. Therefore, due to the marginal preference of InSAR for the dip at 48 • , the forward model was calculated accordingly and used for the comparison with the measured values and the calculation of their residuals in this section. We discuss both models in the corresponding section. The forward wrapped interferograms along with the residuals from the measurements are shown in Figure 2. The forward unwrapped decomposed vertical E-W and N-S deformation fields, along with the residuals for the first two components, are shown in Figure 3. Concerning the participation of each deformation component of the model to the observed deformation along the LOS, we use the LOS unit vectors and the maximum values for each component. We obtain a percentage of 26%, 6% and 68% for the E-W, N-S and vertical component, respectively. Thus, as expected, the LOS values are primarily dominated by the vertical component, secondly by the E-W component and least from the north-south one. This is the reason why, in the decomposition process, omitting this component leads to neglectable errors. Of course, this may not be true for other earthquakes and fault orientations, as was the case of the Samos (Greece) 2020 Earthquake. In a recent publication [20] on an effort to decompose the LOS values, from mean interferograms, using both orbital geometries, the modeled N-S component was considered in the equation. In Figure 4, the modeled and picked deformation values are shown for the ascending and descending orbits, where the misfit of the model and the measurements could be controlled.
In terms of earthquake size, the geodetic moment was calculated to 1.6 × 10 16 N·m and 1.9 × 10 16 N·m, for the fault parameters favoring the geodetic and seismological dip angle, respectively, assuming a rigidity of the medium of 3.0 × 10 10 Pa. With the dip angle of 48 • , the top-edge of the fault at 0.7 km, and a width of 1.0 km, the root of the fault was located at a depth of 1.5 km and the geodetic centroid depth was at 1.1 km. Regarding the model with the dip angle of 33 • , the top-edge of the fault at 0.5 km and a width of 1.7 km, the root of the fault was located at a depth of 1.4 km, and the geodetic centroid depth was at 1.0 km.

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 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 • E and dip of 35 • ). 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.
A noticeable difference between the focal coordinates of the mainshock obtained by location and the moment tensor inversion was observed ( Figure 5). The centroid is located approximately 2 km ESE and at the same depth. The focal solution of the candidate fault plane (N106 • E/31 • /−87 • ) shows similarity to the one obtained by geodesy (considering the dip constraints imposed by seismology), with the plane bordering what is considered "low-angle" [40]. The seismic moment M 0 = 8.9 × 10 15 N·m (corresponding to M w = 4.6) is also close. The inversion solution is identical with the one by GI-NOA Appl. Sci. 2021, 11, 5947 8 of 20 ( Figure A6), which was also estimated with ISOLA (but with a different set of stations, different frequency range and assuming a line of sources beneath the epicenter). The focal mechanism offered by the Observatoire de la Côte d'Azur (OCA) outlines a similar plane, in terms of strike and dip, but almost strike-slip faulting ( Figure A6); the same is true for the focal mechanism provided by NKUA-SL. However, all three solutions feature a centroid location incompatible with the sequence's distribution ( Figure A6), whereas the location of the herein determined centroid is compatible with the aftershocks' foci.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 21 plane (N106°E/31°/−87°) shows similarity to the one obtained by geodesy (considering the dip constraints imposed by seismology), with the plane bordering what is considered "low-angle" [40]. The seismic moment M0 = 8.9 × 10 15 N·m (corresponding to Mw = 4.6) is also close. The inversion solution is identical with the one by GI-NOA ( Figure A6), which was also estimated with ISOLA (but with a different set of stations, different frequency range and assuming a line of sources beneath the epicenter). The focal mechanism offered by the Observatoire de la Côte d'Azur (OCA) outlines a similar plane, in terms of strike and dip, but almost strike-slip faulting ( Figure A6); the same is true for the focal mechanism provided by NKUA-SL. However, all three solutions feature a centroid location incompatible with the sequence's distribution ( Figure A6), whereas the location of the herein determined centroid is compatible with the aftershocks' foci. Figure 5. Seismotectonic map, showcasing the obtained location results for the mainshock (star) and the epicenters of the aftershock sequence (circles), colored after the focal depth (see color bar). We also present the centroid's surface projection and the determined focal mechanism (beachball), along with the InSAR-inferred rupture area (brown rectangle) and its projected intersection with the surface (brown line). The city of Thiva (square), where the GNSS station THIV (inverted triangle) is installed, is located to the west of the sequence. The villages of Kallithea and Asopia (squares) are located on nearby faults. Faults (black lines) after [8,9]; Sr: Soros, T: Thiva, Ka: Kallithea, As: Asopia, Mv: Mavrovouni, Ar: Asopos river, MS: Megalos Schinos and Er: Erythres. The boundary between the two blocks proposed in [7] is in the SE corner of the area (blue line). Cross-section (dashed line A-A') is presented in Figure 7. The blue rectangle shows the area of Figures  6 and A4. Figure 5. Seismotectonic map, showcasing the obtained location results for the mainshock (star) and the epicenters of the aftershock sequence (circles), colored after the focal depth (see color bar). We also present the centroid's surface projection and the determined focal mechanism (beachball), along with the InSAR-inferred rupture area (brown rectangle) and its projected intersection with the surface (brown line). The city of Thiva (square), where the GNSS station THIV (inverted triangle) is installed, is located to the west of the sequence. The villages of Kallithea and Asopia (squares) are located on nearby faults. Faults (black lines) after [8,9]; Sr: Soros, T: Thiva, Ka: Kallithea, As: Asopia, Mv: Mavrovouni, Ar: Asopos river, MS: Megalos Schinos and Er: Erythres. The boundary between the two blocks proposed in [7] is in the SE corner of the area (blue line). Cross-section (dashed line A-A') is presented in Figure 7. The blue rectangle shows the area of

Discussion
The projection of the fault to the surface at dip angle is parallel, at a distance of ~0.5 km SW, to the Kallithea fault, as the latter is noted in the NOAFaults v3.0 database (Figures 5 and 6). Horizontally, the seismicity forms a well-defined linear area parallel to the Kallithea fault. Moreover, all epicenters are located SW of the suggested fault trace; this confirms that the earthquakes are related to the SW-dipping fault. Concerning the vertical distribution, foci are generally located in shallow depths, down to approximately 5 km, as shown in the cross-section perpendicular to the rupture plane, (its trace is plotted in Figure 6) in Figure 7a. This plane and, by extension, the aligned Kallithea fault, indicates a very vague SW-dipping structure, with the mainshock located shallow and the largest aftershock (ML = 3.3) being at the deepest end. The mainshock is located close to the InSAR-inferred rupture plane. The seismological centroid is east of the geodetic one, right beneath the top edge of the inferred rupture plane (Figures 6 and 7). The model is plotted with contour lines, every 14 mm (half a fringe), with the same color-coding as the underlying wrapped ascending interferogram ( Figure 6). They coincide in the optimum way.
The released seismic moment is well-defined, as both geodetic (M0 = ~1.7·× 10 16 N·m) and seismological (M0 = 0.9 × 10 16 N·m) estimated values in the same order of magnitude, corresponding to Mw = 4.6. According to the determined magnitude, modern empirical laws indicate that the rupture area has a length of 3.2 km, width of 4.5 km and surface of

Discussion
The projection of the fault to the surface at dip angle is parallel, at a distance of 0.5 km SW, to the Kallithea fault, as the latter is noted in the NOAFaults v3.0 database (Figures 5 and 6). Horizontally, the seismicity forms a well-defined linear area parallel to the Kallithea fault. Moreover, all epicenters are located SW of the suggested fault trace; this confirms that the earthquakes are related to the SW-dipping fault. Concerning the vertical distribution, foci are generally located in shallow depths, down to approximately 5 km, as shown in the cross-section perpendicular to the rupture plane, (its trace is plotted in Figure 6) in Figure 7a. This plane and, by extension, the aligned Kallithea fault, indicates a very vague SW-dipping structure, with the mainshock located shallow and the largest aftershock (M L = 3.3) being at the deepest end. The mainshock is located close to the InSARinferred rupture plane. The seismological centroid is east of the geodetic one, right beneath the top edge of the inferred rupture plane (Figures 6 and 7). The model is plotted with contour lines, every 14 mm (half a fringe), with the same color-coding as the underlying wrapped ascending interferogram ( Figure 6). They coincide in the optimum way.
of [42] suggest smaller rupture dimensions, i.e., a surface length of 2.0 km, down-dip width of 3.0 km and a total ruptured surface of 8.0 km. Our geodetic modelling revealed a smaller area with dimensions of width of 1.3 ± 0.6 km and length of 2.0 ± 0.2 km. The geodynamic regime declares a transition zone with extensional stress at an axis with azimuth of N22°E. The latter, in combination with the strike of the rupture, leads to a normal mechanism with a small left-lateral component, as both methods suggest. The rake from seismology was calculated as −87°, suggesting a dip-slip mechanism.  (Figure 3). (b) The ascending interferogram rotated to the same orientation as the cross section, the projection of the inferred fault (its dip favored by geodesy) and with A-A' the cross-section line. Faults (black lines) and their names as in Figure 5. The focal mechanism is projected to the crosssection plane.
The two fault planes that resulted from geodesy ( Figure 7a) are similar ( Table 1). The weak point of InSAR, i.e., to distinguish a definite dip angle, even if it is possible to be discriminated through a better standard deviation, leads us to adopt the solution favoring the dip angle from seismology (Table 1).  (Figure 3). (b) The ascending interferogram rotated to the same orientation as the cross section, the projection of the inferred fault (its dip favored by geodesy) and with A-A' the cross-section line. Faults (black lines) and their names as in Figure 5. The focal mechanism is projected to the cross-section plane.
The released seismic moment is well-defined, as both geodetic (M 0 =~1.7·× 10 16 N·m) and seismological (M 0 = 0.9 × 10 16 N·m) estimated values in the same order of magnitude, corresponding to M w = 4.6. According to the determined magnitude, modern empirical laws indicate that the rupture area has a length of 3.2 km, width of 4.5 km and surface of 14.7 km 2 , as estimated from the relations proposed by [41]. However, the classic relations of [42] suggest smaller rupture dimensions, i.e., a surface length of 2.0 km, down-dip width of 3.0 km and a total ruptured surface of 8.0 km. Our geodetic modelling revealed a smaller area with dimensions of width of 1.3 ± 0.6 km and length of 2.0 ± 0.2 km. The geodynamic regime declares a transition zone with extensional stress at an axis with azimuth of N22 • E. The latter, in combination with the strike of the rupture, leads to a normal mechanism with a small left-lateral component, as both methods suggest. The rake from seismology was calculated as −87 • , suggesting a dip-slip mechanism.
The two fault planes that resulted from geodesy (Figure 7a) are similar ( Table 1). The weak point of InSAR, i.e., to distinguish a definite dip angle, even if it is possible to be discriminated through a better standard deviation, leads us to adopt the solution favoring the dip angle from seismology (Table 1).
InSAR results showed a centroid source at a depth of~1 km (Table 1), shallow enough to induce 2.0 to 2.5 fringes of deformation in the C-band, even with the moderate magnitude of the earthquake. Until recently, InSAR in the Eastern Mediterranean could detect deformation of events as weak as M = 5.0, as in the case of the Gulpinar earthquake [43], located almost 3 km deeper than the event near Thiva. Depths of the mainshock obtained from location (0.9 km) and moment tensor inversion (1.0 km) are the same. As we have repeatedly mentioned, location depths are not highly accurate, due to the absence of local seismological stations, with an average error of 1.5 km. In any case, both seismological methods and geodesy converge on an uncommonly shallow moderate earthquake. In the city of La Habra (CA, USA), geodetic measurements, along with seismological observations, managed to identify a M w = 5.1 event at a centroid depth of 2 km [44]. Very shallow moderate earthquakes have also been reported in France, M w = 4.9, at a 1-km depth [45] and Ecuador, M w = 5.0 at 1.5 km [46].
It has been shown that aftershock zones of dip-slip events in non-subduction environments are not characterized by a preferential up-or down-dip expansion [47]. In our case, we can assume that the mainshock occurred at a shallower depth and seismicity then expanded downdip. While our depth resolution and accuracy are incapable of forming a comprehensive and reliable view about the actual vertical distribution of seismicity, we infer that the aftershocks occur deeper than the mainshock, as they are located to the SW (i.e., down-dip). It is not possible to support this further, as aftershocks are not always caused by ruptures along the fault plane; instead, activation could occur in the damage zone or even further away [48]. Should the plane grossly defined by the aftershocks (i.e., an along-strike length of~5 km and down-dip width of~6 km) be activated as a whole, it could yield an event with a magnitude ranging between 5.0 [41] and 5.4 [42]. We note such magnitude estimations from possible rupture dimensions should be used cautiously due to their empirical nature [49], even though they are broadly used in seismotectonic [50] and seismic hazard [51,52] research.

Conclusions
Analysis of geodetic and seismological data offered evidence for the occurrence of a M w = 4.6 earthquake at a very shallow depth of approximately 2 km, near the city of Thiva (Boeotia, Central Greece). The mainshock and the aftershock sequence seems to have been caused by ruptures on or adjacent to the SW-dipping Kallithea fault. A shallow source, capable of causing a strong event over M = 5.0, so close to an urban area, i.e., Thiva, which is an economic and agricultural hub for Boeotia, should be better studied. Deployment of at least one backbone seismological and GNSS station nearby could greatly improve the depth resolution and better outline the mechanics of local faulting. There are multiple neotectonic structures in close proximity, including the larger Thiva fault. As past experience has taught us, faults with minor (or, even, no) seismic activity in recent times could still be the cause of strong and dangerous earthquakes [53]. Utilizing smaller seismic sequences to map potential sources is valuable for seismic hazard analysis; to improve the quality of such studies, we strongly suggest the integration of geodetic and seismological observations. Author Contributions: Conceptualization, P.E., G.K. and I.P.; data curation, P.E.; formal analysis, P.E., I.S., A.K. and V.S.; investigation, P.E., I.S. and G.K.; methodology, P.E., I.S. and G.K.; project administration, G.K.; resources, P.E., G.K. and I.P.; software, P.E. and I.S.; supervision, P.E.; validation, P.E., I.S. and G.K.; visualization, P.E., I.S., A.K., T.G. and V.S.; writing-original draft, P.E., I.S., G.K. and V.S.; writing-review and editing, P.E., I.S., G.K. and I.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
The authors are grateful to the personnel of all institutes that collaborate in the operation of the Hellenic Unified Seismological Network. Sentinel-1 satellite imagery was acquired through the European Space Agency and Copernicus. THIV GNSS data were provided by Hexagon SmartNET. Figures 1, 5 and A6 were plotted with General Mapping Tools v.6.2 [54] and its Python wrapper PyGMT v.0.3.1 [55]. Figure A5 was plotted with Matplotlib [56]. Figures A7 and A8 were produced by ISOLA. The authors would like to thank the two reviewers for their constructive feedback on the initial manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Geodetic Methods Details
The acquired data covered a period of approximately two months,~30 days prior to the occurrence of the main event (1 November 2020), extending 30 days after (3 January 2021). The raw GNSS data were processed using Bernese v5.2 GNSS software [57]. Data from several GNSS stations of the European Reference Frame (EUREF) were processed together with the daily 30-s RINEX GNSS data. The precise double-difference method was adopted for the analysis. The absolute antenna phase center corrections were used. Several auxiliary files were introduced on the processing procedure. The precise orbital solutions were obtained from the Center for Orbit Determination in Europe (CODE). For the troposphere modelling, the dry and wet Vienna mapping function (VMF) (http: //ggosatm.hg.tuwien.ac.at, accessed on 1 June 2021) and the Neill mapping function were applied. The FES2004 model (http://holt.oso.chalmers.se/loading, accessed on 1 June 2021) was used for the tide loading corrections. Ambiguities were solved using several strategies, depending on the baseline length between the stations, as is the wide-and narrow-lane, and the quasi-ionosphere-free (QIF) strategy. The final daily coordinates were evaluated for the repeatability error on a weekly basis and values were excluded in cases of large deviations from the weekly solution. The formed time series of the three coordinates of THIV station, on the global ITRF2014 reference frame, for the period 1 November 2020 to 3 January 2021, are presented in Figure A1.
Following, we adduce information concerning the quality of our InSAR processing and solutions, including the inversion-related data (Table A1 and Figures A2 and A3) and the descending model ( Figure A4). Following, we adduce information concerning the quality of our InSAR processing and solutions, including the inversion-related data (Table A1 and Figures A2 and A3) and the descending model ( Figure A4).

Appendix B. Quality of Seismological Solutions
For the newly picked events, we used waveforms from all stations of HUSN (the backbone network in Greece) and the Corinth Rift Laboratory Network (CRLnet), which focuses on the Western Gulf of Corinth. Due to the weak magnitudes of the sequence, the majority of the picked phases were in near-local epicentral distances (i.e., a maximum

Appendix B. Quality of Seismological Solutions
For the newly picked events, we used waveforms from all stations of HUSN (the backbone network in Greece) and the Corinth Rift Laboratory Network (CRLnet), which focuses on the Western Gulf of Corinth. Due to the weak magnitudes of the sequence, the majority of the picked phases were in near-local epicentral distances (i.e., a maximum

Appendix B. Quality of Seismological Solutions
For the newly picked events, we used waveforms from all stations of HUSN (the backbone network in Greece) and the Corinth Rift Laboratory Network (CRLnet), which focuses on the Western Gulf of Corinth. Due to the weak magnitudes of the sequence, the majority of the picked phases were in near-local epicentral distances (i.e., a maximum limit of~100 km). We offer the distribution of station distances with a determined P arrival ( Figure A5a) and the magnitudes for the aftershock sequence ( Figure A5b). Note that we did not pick any S phases without the respective P arrival.
i. 2021, 11, x FOR PEER REVIEW 16 of 21 limit of ~100 km). We offer the distribution of station distances with a determined P arrival ( Figure A5a) and the magnitudes for the aftershock sequence ( Figure A5b). Note that we did not pick any S phases without the respective P arrival. To estimate the focal mechanism of the mainshock, we used the ISOLA software for inverting regional recordings ( [33,34]). To this end, we selected an azimuthally satisfying distribution of broadband stations that belong to HUSN ( Figure A6). Short-period instruments (e.g., STFN, with a corner frequency of 1 Hz) were de facto excluded, as the target frequencies range well below 0.5 Hz. In total, 13 stations were initially considered in inversion, but two were rejected. Inversion was carried out at frequencies ranging between 0.05 and 0.10 Hz; the lower boundary was selected by visually inspecting the spectral signal-to-noise ratio (SNR) at each station. For the upper limit, we tried various frequencies between 0.08 and 0.15 Hz to obtain the best inversion results. To estimate the focal mechanism of the mainshock, we used the ISOLA software for inverting regional recordings [33,34]. To this end, we selected an azimuthally satisfying distribution of broadband stations that belong to HUSN ( Figure A6). Short-period instruments (e.g., STFN, with a corner frequency of 1 Hz) were de facto excluded, as the target frequencies range well below 0.5 Hz. In total, 13 stations were initially considered in inversion, but two were rejected. Inversion was carried out at frequencies ranging between 0.05 and 0.10 Hz; the lower boundary was selected by visually inspecting the spectral signal-to-noise ratio (SNR) at each station. For the upper limit, we tried various frequencies between 0.08 and 0.15 Hz to obtain the best inversion results.
The inversion was carried out for multiple sources on a plane striking N111 • E and dipping 33 • , as indicated by the rupture plane estimated through InSAR and the GI-NOA solution. We used 6 sources along-strike and 12 along-dip at 1-km intervals (total of 72 trials). This grid ( Figure A6) was designed to approach the inferred length and width of the Kallithea fault. We note that the source located by Hypoinverse corresponds to source #4 ( Figure A7). Solutions show spatial independence. Solutions beneath the optimal one (#6) showcase stable~WNW-ESE normal characteristics, similar to the 5th km along-strike, down to the 9th km along-dip. However, deeper sources exhibit a rotated strike and a stronger strike-slip component. The depth (analogous to the along-dip dimension) is strongly related to correlation. Shallower depths are accompanied by higher correlation; the selected solution stands at a high 88% correlation. Figure A6. Sites of the stations used in the inversion (triangles for broadband sensors and squares for short-period) for determining the focal mechanism. The three closest rejected stations (MDRA, STFN and VILL, in faded color) are also shown. Beachballs for the solutions of this study, NKUA-SL, GI-NOA and OCA are plotted, with lines pointing to the centroid coordinates provided by each research team. Blue lines bound the blocks proposed by [7]. The study area of Figure 5 (brown rectangle) is noted. Inset: location of the area shown in the main map (dashed black line) in Greece.
The inversion was carried out for multiple sources on a plane striking N111°E and dipping 33°, as indicated by the rupture plane estimated through InSAR and the GI-NOA solution. We used 6 sources along-strike and 12 along-dip at 1-km intervals (total of 72 trials). This grid ( Figure A6) was designed to approach the inferred length and width of the Kallithea fault. We note that the source located by Hypoinverse corresponds to source #4 ( Figure A7). Solutions show spatial independence. Solutions beneath the optimal one (#6) showcase stable ~WNW-ESE normal characteristics, similar to the 5th km alongstrike, down to the 9th km along-dip. However, deeper sources exhibit a rotated strike and a stronger strike-slip component. The depth (analogous to the along-dip dimension) is strongly related to correlation. Shallower depths are accompanied by higher correlation; the selected solution stands at a high 88% correlation. Figure A6. Sites of the stations used in the inversion (triangles for broadband sensors and squares for short-period) for determining the focal mechanism. The three closest rejected stations (MDRA, STFN and VILL, in faded color) are also shown. Beachballs for the solutions of this study, NKUA-SL, GI-NOA and OCA are plotted, with lines pointing to the centroid coordinates provided by each research team. Blue lines bound the blocks proposed by [7]. The study area of Figure 5 (brown rectangle) is noted. Inset: location of the area shown in the main map (dashed black line) in Greece.
Synthetics were calculated through ISOLA, using the discrete wavenumber method [58] and the regional velocity model of [31]. The similarity with observed signals ( Figure A8) is exceptional. Low amplitude waveforms in MDRA and VILL, probably due to the proximity of the stations to the epicenter, were excluded from the inversion. Figure A7 offers variance reduction (VR) as a metric of similarity. VR = 1 − ∑ r 2 ∑ d 2 , where r the residuals between observed and synthetic amplitudes and d the actual observed data [33]. The average VR from all traces is 0.78. In combination with the high double-couple (DC) percentage contribution to the moment tensor (87%) and the high correlation (88%), the final solution can be considered very reliable. Note that we used a deviatoric-constrained inversion type. Synthetics were calculated through ISOLA, using the discrete wavenumber method [58] and the regional velocity model of [31]. The similarity with observed signals ( Figure  A8) is exceptional. Low amplitude waveforms in MDRA and VILL, probably due to the proximity of the stations to the epicenter, were excluded from the inversion. Figure A7 offers variance reduction ( ) as a metric of similarity. 1 ∑ ∑ , where the residuals between observed and synthetic amplitudes and the actual observed data [33]. The average from all traces is 0.78. In combination with the high double-couple (DC) percentage contribution to the moment tensor (87%) and the high correlation (88%), the final solution can be considered very reliable. Note that we used a deviatoricconstrained inversion type.