## 5 Discussion  
We designed and developed an algorithm named HYDE, able to estimate the geometry (size, depth, dip and strike angles) of fault planes using the hypocenter locations of main earthquakes and the associated aftershocks, even where the seismogenic sources are fragmented. Our HYDE algorithm is still capable of accurately identifying fault parameters even when the fault geometry is highly non-planar or complex. In cases of fault segmentation, when there are multiple contiguous planar faults with appreciably different strikes and/or dips, the algorithm will first identify the fault segment that contains the largest number of earthquakes. Then, by excluding the earthquakes belonging to this primary fault, the algorithm searches for the next fault cluster, and so on, providing the respective strike, dip and mean location for each identified fault. In cases where fault substructures are orthogonal to the main fault, the algorithm ignores them (along with the associated earthquakes) because they are assigned a weight of zero, as they lie outside the explorer parallelepiped that defines the main fault. When the fault structure takes the form of a cylindrical surface sector, we need to express the earthquake locations in cylindrical coordinates. In this system, the distance from the curved fault surface corresponds to the r-coordinate of the 3D reference system. The analog of a parallelepiped in Cartesian coordinates is a volume contained between two concentric sectors. Many pivots will be attracted to this sector, as in the planar fault case; however, we now need to explore two additional parameters: the cylindrical curvature and the rotation of the symmetry axis. In this scenario, the computing time increases by a factor that depends on the number of discretized curvature values multiplied by the number of discretized rotations for the symmetry axis. For the 2009 L'Aquila seismic sequence, involving approximately 64,000 earthquakes analyzed using Cartesian coordinates, we used 1000 pivots to achieve comprehensive coverage of the parameter space. When run on an Intel i7 PC with 16 cores, the computing time was approximately 10 days. Figure 10: The whole seismicity of the 2009 earthquakes in the L’Aquila area. (**a**) The explorer parallelepiped containing hypocenters. (**b**) Red indicates the hypocenters belonging to the explorer parallelepiped, which individuates the main fault plane also containing the main shock (Mw 6.1). The number of pivots to choose depends on the ratio between the number of earthquakes associated with the fault and the number of localization points outside the defined parallelepiped. However, assuming we have 1000 positive events (within the parallelepiped) and 9000 negative events, the probability of obtaining a null extraction from the fault is approximately \(10^{-5}\), a very low probability. Consequently, the random selection of pivots will likely include around 100 pivots on the fault, statistically distributed according to the density of earthquake hypocenters. In our case (see Figure 8b), about 15,000 to 20,000 earthquakes are associated with the fault, which is a significant fraction of the 64,000 events considered. Using 200 pivots yielded practically the same result. Obviously, determining the optimal number of pivots must be evaluated on a case-by-case basis; however, the fault tends to act as an attractor for the pivots. Different considerations are necessary for substructures with fewer events. In these cases, it is advisable to first isolate the involved volume and then rerun the code--in that subspace, the same considerations still apply. Our algorithm is particularly useful where the earthquake intensity is not so high that the radiation pattern can be easily determined, as often happens for most aftershock sequences or in volcanic areas, where the percentage of earthquakes for which the radiation pattern can be defined is very low. One of the advantages of the proposed method is that it does not require focal mechanisms at all.