The Gasperini method (Gasperini et al., 1999; Gasperini et al., 2010) (Step 1 in Figure 2) uses the distribution of the intensity data points to localize the epicentral area of a historical earthquake and estimate the dimensions and orientation of its most likely rupture area and its moment magnitude (M_w). This method generates an oriented rectangle, named “boxer”, which is meant to represent the surface projection of the fault responsible for the earthquake or the portion of the Earth where it would be likely to be located. The boxer’s dimensions derive from the magnitude using empirical relations proposed by several authors (Wells and Coppersmith, 1994; Stirling et al., 2002; Wesnousky, 2008; Hanks and Bakun, 2008).

Comparing the projection of the obtained earthquake area source with the geological information (i.e., active faults) may help to define those faults more suitable to be responsible for the event. This may allow discarding as possible earthquake sources the geological faults that lie too far from the boxer and/or those whose dimensions and strike differ significantly from it. The faults considered as probable seismic sources are then taken as candidates and modeled in the simulation stage (see Seismic Scenario method) (step 2 in Figure 2). The calculations have been done using the Boxer software (Gasperini et al., 1999; Gasperini et al., 2010) calibrated with specific coefficients for the Betics region (Gomez-Capera et al., 2014).

The seismic scenario method (de Pro-Díaz et al., 2022) (Step 2 in Figure 2) builds earthquake scenarios that simulate the distribution of the intensity field related to each of the candidate fault sources identified with the Gasperini method. Then, the obtained intensity fields are compared with the observed intensity field to see which one fits better the observed data.

To model the seismic scenarios this method uses the OpenQuake software (Pagani et al., 2014), which requires as input several source parameters: (1) earthquake magnitude, whose source is explained later in this paragraph; (2) fault geometry and location, which are obtained from the geological maps (e.g., García-Mayordomo et al., 2012) and previous works (e.g., Insua Arévalo, 2008; Martínez-Díaz and Hernández-Enrile, 2004); (3) hypocentral depth and coordinates, whose source is also explained later in this paragraph; and (4) upper and lower seismogenic limits of the candidate fault, which generally range between 0 and 15 km (García-Mayordomo et al., 2012, and references therein). Magnitude is calculated through empirical relations considering the size of the rupture; the relations we used are cited at the end of this section 3.2, as well as the criteria to select them. Hypocentral depth and coordinates are obtained from different bibliographic sources, such as earthquake catalogues (e.g., Martínez Solares and Mezcua Rodríguez, 2002) or other authors’ works searching for the source of the earthquake (e.g., Muñoz and Udías, 1988); whenever we find discrepancies in the bibliography about this parameter, we try several depths in the simulations to test which one produces the best results. With these parameters, OpenQuake calculates the ground motion in the study area, and the ground motion values can then be transformed into intensities using ground-motion-to-intensity conversion equations. To account for amplifications in the ground motion caused by site effects, OpenQuake allows the user to input a Vs.30 model. OpenQuake also allows the user to input a default uniform Vs.30 value for the whole study area. Because each lithology has different Vs.30 values and due to the spatial variability of lithologies in the Betics, we prefer to use the Vs.30 models provided by Allen and Wald (2007) and Rodríguez-Peces et al. (2012).

Finally, to predict ground-motion parameters such as peak ground velocity (PGV) or peak ground acceleration (PGA) it is necessary to use a ground motion model (GMM). Ground motion is usually considered to increase with magnitude and decrease with increasing distance to the seismic source; but it also has an aleatory standard deviation associated with the Gaussian distribution around the mean value predicted by a function of magnitude and distance, the sigma (e.g., Reiter, 1991). Using a GMM, OpenQuake generates a regular grid over the study area, where each cell contains simulated ground motion values for both PGA and PGV. We performed a total of 100 different simulations for each rupture, so each point on the grid contains 100 simulated values for PGA and PGV. We calculated the mean of the logarithms of the 100 values for each point on the grid, since ground motion models show a lognormal distribution. Due to the nonlinearity of the site amplification functions used in GMMs and the aleatory variable sigma, it is better to use 100 different simulations instead of using the mean ground-motion field for each location (de Pro-Díaz et al., 2022). In a single mean simulation we can observe regions with abnormally high or low ground-motion values considering the values in the surrounding area, but we have no way of knowing if this anomaly is caused by the random sigma or by site effects. If the same “anomalous” values appear in that area in 100 simulations, it is much more likely that they are due to actual site effects, so the values sampled from the 100-averaged simulations are considered more representative of the ground-motion generated by the earthquake and its possible site effects amplification. Then, with the help of a ground-motion-to-intensity conversion equation (GMICE) suitable for our area, we transform the grids of simulated ground motion values (PGV and PGA) into grids of intensity values. All the grids resulted from the simulations are crossed with the observed intensity data points, and residuals are obtained as follows: R o b s − r u p = I o b s – I r u p where I o b s is the observed intensity value and I r u p is the spatially coincident simulated intensity value in the grid.

To build accurate seismic scenarios it is essential to use GMM and GMICE developed for the area in which the earthquake occurred, since the crustal composition, the building style and people’s perception of an earthquake can greatly vary from one region to another. However, there are no GMM or GMICE specifically designed for south Spain. Consequently, we needed to carry on a calibration test combining different GMM and GMICE to determine which pair performs best in this area. To do so, we have carried out various seismic scenario simulations for the 2011 Lorca earthquake (Figure 3), which was generated by the Alhama de Murcia Fault (Martínez-Díaz et al., 2012a; Frontera et al., 2012), combining different GMM and GMICE and comparing the simulated intensities to the actual distribution of observed intensities (see section 4.1). On the scenarios of the 1804 Dalías and 1680 Málaga earthquakes, we applied the GMM-GMICE pair that provides better correlations.

The 2011 Lorca earthquake was also used to test several empirical relations for the magnitude calculation of an earthquake from its rupture size. We calculated the magnitude based on the rupture area constrained by Martínez-Díaz et al. (2012a) and Frontera et al. (2012) using equations from Wesnousky (2008), Hanks and Bakun (2008), Stirling et al. (2002) and Wells and Coppersmith (1994), and then compared the results to the actual Mw of the event to see which result was closer to it.

When seismic scenarios from two or more candidate sources show R_obs-rup equally close to 0, de Pro-Díaz et al. (2022) propose an extra step in the methodology to distinguish the better candidate (Step 3 in Figure 2). First, it is necessary to compare the two competing candidate scenarios to identify the areas where the intensity values differ. To do this, we have to subtract one scenario’s raster from the other, resulting a new raster that may show a large area with zero values and smaller areas with positive and negative values (in which one scenario has higher intensity values than the other). From this point onwards, only the areas with positive and negative values, which are automatically identified with GIS software, are considered in the analysis. Then, for each of these areas, a Kolmogorov-Smirnov statistical test (K-S) (Massey, 1951) is done to compare the likeness of two data distributions. The K–S is a nonparametric test that can evaluate whether two datasets belong to the same distribution. This test calculates the maximum distance (D-value) between both empirical cumulative distributions, and compares it with the size-dependent value Dc(α) obtained for the particular significance level α. The null hypothesis (H₀) assumes that the two considered data samples follow the same distribution. In addition, a p-value is calculated and it represents the probability of having the observed D-value from randomness, assuming that H₀ is true. We use the typical value of 5% significance level to evaluate the results. When comparing intensity distributions of the two scenarios, if H₀ is accepted it means that the two scenarios are still too similar in the sampling points and we cannot statistically distinguish the best candidate, so this method cannot be fully applied. However, if H₀ is rejected in this step, we can continue on with the analysis. Finally, the K-S test is carried out again to compare the simulated intensity distribution of each scenario with the observed intensity distribution. Presumably, there will be one scenario were H₀ is accepted and, then, the candidate earthquake fault source after which this scenario was modeled may be the best candidate to explain the earthquake effects and, thus, the closest to the fault that produced that specific event.

Although of moderate magnitude (Mw 5.2), the 2011 Lorca earthquake produced nine casualties and large economic losses in Lorca (Martínez-Díaz et al., 2012a; Benito Oterino et al., 2012; Rodríguez-Pascua et al., 2012), and it was felt in a large area (Figure 3). The shallowness of the earthquake focus (4.6 km), its proximity to the city of Lorca (5.5 km) and the directivity of the rupture propagation that induced maximum PGA values larger than expected, may explain its catastrophic consequences (López-Comino et al., 2012). The analysis of the coseismic deformation using radar interferometry (INSAR), as well as the seismological data, allowed Martínez-Díaz et al. (2012b) and Frontera et al. (2012) to localize and model the causative fault and relate the earthquake to the Alhama de Murcia fault (AMF) (red rectangle in Figure 3B), one of the major tectonic structures of this area, which is known to have caused large earthquakes of M > 6.0 in both historical and pre-historical times (Masana et al., 2004; Martínez-Díaz et al., 2019; Martínez-Díaz et al., 2012a; Ortuño et al., 2012).

Using the Gasperini method, we have modeled the probable rupture area of the 2011 Lorca earthquake from the macroseismic data presented by Martínez-Díaz et al. (2012a) (purple rectangle in Figure 3B and Table 1). The computed epicenter is localized at 37.6815° N and 1.6905° W, centered in the city of Lorca, and the resulting earthquake source has an area of 8.84 km², strikes N75° E and the obtained estimated magnitude is Mw 4.83 ± 0.20. The modeled epicenter is localized 3.6 km towards the SSW of the instrumental epicenter and the rupture area lies on top of the AMF, partially overlapping with the rupture derived from INSAR and seismological data (red rectangle in Figure 3B). Even though the rupture area modeled with the Gasperini method appears to be slightly displaced to the southwest and a bit smaller than the INSAR-derived area, there is a rather good correlation in both size and strike. The discrepancy on the location of both sources may be related to the bias in the distribution of the intensity data points in the epicentral area, most of them restricted to the southeast of the fault where the town of Lorca is located. Considering the good agreement between the source modeled from the intensity data points and the seismological and INSAR data, we consider that the Gasperini method using the coefficients calibrated for the Betics region (Gomez-Capera et al., 2014) may be a suitable first approach to select the candidate faults that could have produced the Malaga 1680 and the Dalías 1804 earthquakes.

In order to select the GMM and the GMICE that better represents the study area, we have built seismic scenarios using the fault parameters corresponding to a simplified AMF (Table 2). According to Quirós Hernández (2017), we have tested the two GMM recommended for this zone, Abrahamson et al. (2014) and Campbell and Bozorgnia (2014). For the GMICE, we have tested five global conversion equations: Caprio et al. (2015), Wald et al. (1999), Tselentis and Danciu (2008), Atkinson and Kaka (2007), and Worden et al. (2012). Residuals for the different combinations of GMM and GMICE for the Lorca earthquake are shown in Figure 4. The closer the residual is to 0, the best the simulated seismic scenario fits the observed intensity field. Both GMMs used show good results (Figure 4), although Campbell and Bozorgnia (2014) shows residuals slightly closer to 0. Both models were built based on similar databases and share almost identical applicability ranges, although Campbell and Bozorgnia (2014) covers a larger range of Vs.30 values for the site amplification: from 150 to 1500 m/s, in contraposition to the 180—900 m/s for Abrahamson et al. (2014). The range of Vs.30 values in the study area is 150–1500 m/s. For the sake of avoiding duplicity of scenarios, for its better performance with the Lorca earthquake (even if the difference with Abrahamson et al. (2014)’s GMM was little) and for its larger range of applicability, we chose to use Campbell and Bozorgnia (2014)’s GMM to generate the ground-motion fields for the studied historical earthquakes. For the GMICE, the best performing equation, the one with R closer to zero, is Worden et al. (2012) (Figure 4). According to these results, we have chosen the GMICE proposed by Worden et al. (2012) to calculate the intensity values in the analysis of the Málaga and Dalías earthquakes.

We also used the Lorca earthquake to test several empirical relations to estimate the magnitude of an earthquake from the size of the rupture. As has been stated before, the Mw for this event was 5.2 (López-Comino et al., 2012). Using the equations in Hanks and Bakun (2008), we obtained Mw 5.0. Wells and Coppersmith (1994) gave Mw 5.7 when using rupture length as input and 5.1 when using rupture area. Stirling et al. (2002)’s equations resulted in Mw 6.3 when using the length of the rupture as the input parameter and Mw 5.8 when using the rupture area. With Wesnousky (2008), the calculated Mw was 6.0. We decided to use Hanks and Bakun (2008), Wells and Coppersmith (1994), and Stirling et al. (2002), using rupture area as the preferred input parameter instead of rupture length. When deciding the input Mw for the Málaga and Dalías simulations, we consider the magnitude proposed in the bibliography as the minimum possible magnitude. If this minimum magnitude does not produce simulations with high enough intensities to match the observed intensity field, we calculate Mw from the size of the rupture modeled in OpenQuake (which does not necessarily have to match the size of the boxer) using the aforementioned size-magnitude scaling equations.

On 25 August 1804, the Dalías earthquake caused important damage as well as geological effects (liquefaction, rock falls and hydrogeological anomalies, among others) all around the Campo de Dalías area (Murphy Corella, 2019) (Figure 1; Figure 5). There are 83 intensity data points available from Murphy Corella (2019) for this earthquake, 47 of them corresponding to geological effects with intensity values assigned in the ESI-07 scale (I_ESI henceforth) and the rest corresponding to EMS-98 intensity values (I_EMS henceforth). Different studies have estimated that its magnitude would be 6.4-6.6 and have located its epicenter either inland (Figure 5) or in an unspecified offshore location, but its seismic fault source remains unclear (Espinar Moreno, 1994; Martínez Solares, 2011; Huerta et al., 2015; Murphy Corella, 2019).

Geologically, Campo de Dalías is one of the Neogene basins controlled by NW-SE faults. The sedimentary sequence, which lies over the Alpujárride metamorphic basement, goes from the Upper Tortonian to the Holocene showing a process of marine regression, and is composed of sedimentary rocks deformed by subtle E-W trending, high-wavelength folds and by neotectonic faulting (Rodríguez-Fernández and Martín-Penela, 1993; Martínez-Díaz, 1999; Sanz de Galdeano and Alfaro, 2004; Marín-Lechado et al., 2007; Pedrera et al., 2012). The Alpujárride basement outcrops outside Campo de Dalías, in Sierra de Gádor, the relief which delimits the basin to the north (Sanz de Galdeano and Alfaro, 2004; Marín-Lechado et al., 2007; Pedrera et al., 2012) (Figure 5C). Tectonically, Campo de Dalías has a complex system of faults and rotating blocks of crust delimited by the CF offshore to the SE and the AFZ inland to the N (Figure 5 in this work; figure 11 in Martínez-Díaz and Hernández-Enrile, 2004). These blocks, delimited by oblique faults, are part of a bigger, wedge-shaped crustal block which may escape westwards (Martínez-Díaz and Hernández-Enrile, 2004). Among these faults is the Loma del Viento Fault (LVF), a NW-SE normal-dextral segmented fault whose inland mapped trace is around 8 km long (Marín-Lechado et al., 2005; García-Mayordomo et al., 2012), but according to Murphy Corella (2019) and Pedrera et al. (2012) it could actually extend both offshore and inland to reach more than 30 km in total length. Mostly parallel to the LVF, but at ∼3 km to the north-east from its trace is the Llano del Águila Fault (LLAF), a 20-km long normal-dextral fault which has been recently mapped in detail by Molins-Vigatà et al. (2022). The Balanegra Fault (BF) is another of these normal-dextral faults delimiting the rotating blocks, which is localized to the south of LVF and with a total inland length of around 9 km (Martínez-Díaz and Hernández-Enrile, 2004; Marín-Lechado et al., 2010; Sanz de Galdeano et al., 2020). The distribution of the intensity data points (Figure 5) suggests that either of these three faults could have generated the 1804 Dalías earthquake.

Since there is a significant amount of I_ESI points in the intensity field of this earthquake, we used the Gasperini method twice, one with the whole intensity field (I_EMS+ESI dataset henceforth) and the other with only the I_EMS points. The aim was to check whether or not the addition of the I_ESI points to the dataset changes the results on this step of the analysis. In both cases, the source location for the Dalías earthquake obtained using the Gasperini method is estimated northwest of the Campo de Dalías, between the LVF and the LLAF (Figure 5). For the I_EMS+ESI dataset, this method has provided an epicenter 36.82° N and 2.85° W, a magnitude of 5.96 ± 0.43, and a rupture area of 94.8 km² with a strike of N111° E, similar to the LVF and the LLAF (boxer A in Figure 5B and Table 1). For the I_EMS dataset, the method provides an epicenter 36.79° N and 2.82° W, a magnitude of 6.06 ± 0.86, and a rupture area of 116.45 km² with a strike of N119° E, also quite close to the strike of the local faults (boxer B in Figure 5B and Table 1). The geometric characteristics, such as length, size and strike, and the location of the modeled source areas suggest that the LVF or the LLAF may be the best fault candidates to have generated the 1804 Dalías earthquake. Even though the BF’s strike is similar to the modeled source areas, its position far to the south may discard it as a possible source for this earthquake. Finally, we have discarded the CF as possible candidate due to the distance from the modeled source and the disagreement in strike. In agreement to all these observations, the preferred rupture candidates are the LVF and/or the LLAF. For LVF we have considered three different possible ruptures scenarios (Table 2; Figure 5) and two of the LLAF (Figure 5):

- Scenario D1: LVF extended towards the sea following descriptions by Murphy Corella (2019) and the model of Martínez-Díaz and Hernández-Enrile (2004) (rupture A).

The seismic scenario simulations for the Dalías earthquake for each candidate fault rupture (Figure 6) show that there is an evident spatial discrepancy between the distribution of I_obs values and the pattern of simulated intensities for Scenario D1 (Figure 6A). In addition, the analysis of the residuals show that the maximum intensity calculated for this scenario D1 is almost 2 degrees lower than the maximum observed intensity on average (Figure 7). These observations suggest that this source may be discarded as a possible source of the earthquake. The scenarios considering inland ruptures (scenarios D2-D5) show better correlation with the observed intensity data points (Figures 6B–E); however, the analysis of the residuals corresponding to the intensities derived from PGA and PGV indicate that overall the modeled intensities are lower than the observed ones (Figure 7). The comparison of the residuals between these four scenarios show that the lower residuals, which are less than 1 intensity degree, are obtained for scenarios D2 and D5 (Figure 7). In contrast, scenarios D3 and, specially, D4 result in intensities a little bit much lower than the previous ones (Figure 7). These differences are consistent with the spatial distribution of intensities in map view (Figures 6B–E), suggesting that scenarios D2 and D5 are the more plausible ones. Accordingly, we selected these two scenarios to carry out the spatial difference analysis to define the best scenario.

The spatial intensity difference map (Figure 8) shows the areas in which the preferred scenarios D2 and D5 differ from each other, as well as the 14 intensity data points located inside these areas (“useful points” henceforth). Any intensity data point outside of these areas is not considered useful for the analysis, as stated in the Methodology section. Scenario D2 produces lower intensities than scenario D5 towards the northeast of the interest area, whereas towards the south the values are rather similar between both scenarios. Histograms on Figure 8 show the distribution of frequencies of the intensity values sampled from the useful points for scenarios D2 (I_rupB) and D5 (I_rupE), as well as the distribution of frequencies of the observed intensities (I_obs). Scenario D5’s distribution seems to be closer to the distribution of I_obs. In contrast, scenario D2 is slightly off towards lower intensity values. However, the distribution of frequencies does not completely match neither for scenario D2, nor D5. Then, visual inspection of the histograms is not conclusive to select between scenarios.

In the areas where the results of scenario D2 and scenario D5 differ we have first used the K-S tests to evaluate if the intensity values sampled in the useful points from both scenarios belong to the same distribution or not. In the case of the scenarios calculated using PGA, the null hypothesis H₀ is rejected at 5% significance level. This means that the two distributions are statistically different enough and thus the last step of the analysis can be carried out. We then compare I_obs with the simulated intensities for each of the two scenarios using the K-S test. In this last step, H₀ is rejected at the 5% significance level for scenario D2 and accepted for scenario D5. These results are presented in Table 3 and may suggest that scenario D5 provides the most consistent distribution of intensities in comparison to the observed data.

The 1680 Málaga earthquake produced damages in the province of Málaga and was felt in a larger area (Figure 9). Its estimated magnitude ranges between 6.6 and 6.8 (Muñoz and Udías, 1988; Martínez Solares and Mezcua Rodríguez, 2002). There is a total of 35 intensity data points available for this earthquake from Silva Barroso et al. (2019), eight based on the evaluation of geological effects in the ESI-07 scale and the others provided by regular macroseismic analysis considering also archaeological effects in the EMS-98 and MSK-64 scales. Even with all this information, there is still uncertainty about the source of the earthquake. Goded et al. (2008) suggested an inland location, although they did not point to any specific fault. Alternatively, Mezcua et al. (2013) pointed at the El Acebuchal and Los Alamillos Faults as possible sources (Figure 9). Several other faults have been mapped in the area, such as the Montes de Málaga Fault (MMF) and the Cártama, Mijas and Villafranco de Guadalhorce (VDG) Faults (Figure 9) (García-Mayordomo et al., 2012). These faults trend approximately E-W and Insua Arévalo (2008) proposed a structural model consisting of blind thrust fault strands whose dip varies along ramps and flats with all strands merging at depth (Figure 9). Although it has not been mapped nor named yet, it has been proposed that there is a third blind thrust to the east of the Mijas system with similar length and strike as the Cártama and VDG Faults (“blind fault” in Figure 9) (Insua Arévalo, 2008). Geologically, the Málaga basin consists of a sequence of regression-transgression-regression in its sedimentary record ranging from Tortonian to Quaternary and lying over the Alpujárride and Maláguide metamorphic basements (Sanz de Galdeano and López-Garrido, 1992; López-Garrido and Sanz de Galdeano, 1999; Guerra-Merchán et al., 2000). The metamorphic basement outcrops to the north and south of this basin, and was formed by a superimposition of the Maláguide and Alpujárride units, both of which have a complex internal structure (Vera et al., 2004; Insua Arévalo, 2008) (Figure 9B).

For the Málaga earthquake we used the Gasperini method with two different datasets. The first one only considered the EMS-98 intensities compiled by Goded et al. (2008) (area A henceforth) and the second one corresponded to the intensities provided in Silva Barroso et al. (2019) (area B henceforth). Both modeled seismic sources overlap and the location of the computed epicenters is rather similar (Figure 9; Table 1). However, there is a significant divergence with the dimensions of the rupture area and its strike. The area obtained for the A dataset measures 216.24 km² and strikes N94° E, whereas for the B dataset measures 303.45 km² and shows a trend of N124° E. In agreement, the change in dimensions of the rupture areas also implies a difference in the resulting estimated magnitude for both ruptures, being 6.36 ± 0.19 and 6.52 ± 0.05 for A and B datasets, respectively. When compared with the faults mapped in the study area, both modeled ruptures partially or totally overlap with the fault system consisting of the VDG, Cártama and Mijas Faults, as well as the southern MMF faults. However, area A’s strike is in better agreement with the general trend of the local faults. Considering the relationship between both modeled ruptures and the mapped faults, we selected seven possible rupture scenarios. In addition, we discarded the El Acebuchal and Los Alamillos Faults as possible fault candidates because they did not match any of the area source solutions neither in size or strike. The modeled fault ruptures are shown in the cross-section in Figure 9D and the used parameters for each rupture are summarized in Table 2. The seven scenarios are (Figure 9D; Figure 10):

- Scenario M1: Conjunct rupture of the Cártama, VFG and the third blind fault to the east considering the flat and frontal ramp geometry, and not involving the deeper thrust ramp corresponding to the Mijas fault, with a hypocentral depth of 2.5 km (rupture A in Figure 9D).

The spatial distribution of simulated intensities (Figure 10) does not allow for a clear rejection or candidate selection from the defined scenarios because of the mismatch of intensities below VII between observed and simulated intensities. The comparison of the distribution between the maximum simulated intensity area and maximum observed intensities (Figure 10) suggests that all the scenarios may seem plausible, with the exception for M4. In addition, the analysis of the residuals for each of the scenarios shows that the simulated average intensities are between 1.6 and 0.3 degrees lower than the observed ones (Figure 11). However, the best residuals are provided by scenarios M6 and M7, which are around 0.5–0.3 degrees lower than the observed intensities. Scenarios M1, M2 and M4 produce average intensities 1.2–1.6 degrees lower than the observed ones, which may suggest that ruptures A, B and D are the less plausible candidates. Finally, scenarios M3 and M5 show average intensities almost 1 degrees lower than the observed intensity, which makes ruptures C and E unfitting to be considered good candidates either. This would point to ruptures G and F as the preferred candidates, in that order. We thus move on to the next step of the analysis with candidate ruptures G and F.

The spatial intensity difference map for this earthquake (Figure 12) shows that in some areas scenario M7 produces higher intensities than M6 (blue polygons in Figure 12) and the contrary in others (red polygons in Figure 12). Unfortunately, there are only four useful points inside these areas from the observed intensities data points, which are not enough to perform a robust statistical test. Thus, the spatial intensity step of the analysis could not be carried out for the Málaga earthquake.

We have used the Lorca earthquake as a calibration event to test the accuracy of both the Gasperini method and several magnitude-rupture size regressions, as well as to establish the best performing combination of GMM-GMICE in the study region. The 2011 Lorca earthquake was the strongest event occurred during the instrumental period at a shallow hypocentral depth. Despite its moderate magnitude (Mw 5.2), the Lorca event generated a notable volume of intensity data that is suitable for intensity distribution analysis. In addition, its seismic source has already been determined in earlier studies allowing the comparison with the results obtained in the present study. Because of this, we consider the Lorca earthquake as the most reliable calibration event.

The magnitude of the Lorca earthquake is significantly lower than the ones estimated for the 1680 Málaga and the 1804 Dalías earthquakes and this could have had some influence on the results of the calibration for choosing the GMICE. The GMICEs of Caprio et al. (2015) and Tselentis and Danciu (2008), which were both designed using Mediterranean data, performed worse than Worden et al. (2012), which was built based mainly on North American data. While the differences between GMICEs’ residuals in the calibration step were notable, the difference between the results of the two GMMs during the calibration was quite subtle (Figure 4), so much that we do not consider an increment on magnitude may aggrandize these differences notably. In addition to this, both GMMs have been recommended for this region and used by Quirós Hernández (2017). Because of this, we consider the difference in magnitude between the Lorca earthquake and the other two studied events has not had a relevant influence in the results of the GMM calibration.

The equations to estimate magnitude from rupture size that performed better with the Lorca earthquake were Hanks and Bakun (2008), Wells and Coppersmith (1994) and Stirling et al. (2002), using rupture area instead of rupture length for the latter. Stirling et al. (2013) did a shortlist of scaling magnitude equations from several authors for each tectonic environment and recommend Hanks and Bakun (2008) for strike-slip faults only (such as the AMF, source of the Lorca earthquake). Stirling et al. (2013) considered Wells and Coppersmith (1994) and Stirling et al. (2002) could be used globally in different geodynamic contexts and applied to any fault regardless of its kinematics. Since our candidate faults were normal for the Dalías earthquake and reverse for the Málaga one, we preferred using global equations instead of equations recommended for particular kinematics. We also chose to rely more on the equations by each author that consider rupture area instead of length, since, especially in the Málaga case, we have several candidate ruptures with the same length but different areas.

In the case of the Dalías earthquake, our results show that the candidates which fit better the observed intensity data are all inland faults, agreeing with Espinar Moreno (1994). The best performing simulations, D2 and D5, were the only ones that generated intensities high enough to match observed ones. The magnitude derived for these two scenarios taking into account the modeled source dimensions ranges between Mw 6.9 and 7.1 for D2 and D5, respectively. Although the obtained magnitude is higher than the Mw 6.4-6.6 proposed previously (Martínez Solares, 2011; Huerta et al., 2015; Murphy Corella, 2019), we prefer the newly estimated magnitudes since the resulting intensity distribution fits better the historical and geological observations.

The source area modeled using the dataset that includes the geological effects (Murphy Corella, 2019) slightly differs in size, strike and position from the one using only macroseismic data (Figure 5; Table 1). Considering that the Gasperini method uses the points with the highest degrees of intensity minus one degree (I_max-1) to calculate position and strike of the source area, these differences may be related to the number and distribution of the maximum intensity data points. Although in both datasets the spatial distributions of the points along Campo de Dalías and the surrounding area is similar, the dataset that combines the macroseismic and the geological effects data includes 24 locations with an assigned intensity of VIII, per 13 locations when using just I_EMS data. Knowing that the Gasperini method results highly depend on the high intensity data point spatial distribution and that earthquake geological effects are usually localized in the areas that have experienced large ground motion, it seems advisable to use the geological effects, when available, to better constrain the modeled earthquake source.

Mw calculated with Boxer for the I_EMS+ESI dataset was 5.96 ± 0.43, while for the I_EMS dataset the result was Mw 6.06 ± 0.86. Both calculated magnitudes are quite similar considering their respective error margins, and they are both also lower than the magnitudes proposed by independent authors in the bibliography. For this reason, we consider Boxer to be clearly underestimating this earthquake’s magnitude. This underestimation could be related to the asymmetry in the observed intensity field: as it is shown in Figure 5, this earthquake occurred in a coastal area, so there is a complete lack of data in the marine area to the southeast. This asymmetry might make it difficult to assess the exact magnitude of this earthquake, and could have also affected the calculations of Huerta et al. (2015) and Martínez Solares (2011), since they used macroseismic intensities to estimate the magnitude as well. Murphy Corella (2019) partially evaded this issue by calculating the magnitude from the length of his proposed source rupture, which is slightly smaller than our preferred source.

Intensities derived from PGV performed similarly to those derived from PGA in the simulations for the Dalías earthquake, although the residuals show a higher standard deviation in the simulations derived from PGV (Figure 7). This tendency is consistent with what we have seen in the control earthquake of Lorca, where scenarios based on PGA showed residuals with overall lower standard deviations than their PGV-based counterparts (Figure 4). In the case of the Lorca earthquake, the PGA-based scenarios also showed average residuals closer to 0 than their PGV-based counterpart (Figure 4), which is also true for the D1, D4 and D5 scenarios in the Dalías earthquake (Figure 7). This could mean that PGA-based scenarios in this area show a more realistic intensity distribution, regardless of their similarity with the observed intensity field.

The spatial analysis of the residuals for rupture B and rupture E scenarios shows that the simulations are quite similar to each other, although there is an area (blue polygons in Figure 8) where intensities produced by rupture E are higher and more similar to the reported values. In addition, the histograms in Figure 8 and the results from the K-S test indicate that the best candidate among the five analyzed candidates may be rupture E. Accordingly, the combined rupture of the Loma del Viento and Llano del Águila Faults may be the closest to the actual seismic source of the 1804 Dalías earthquake. Our results agree with Murphy Corella (2019), who proposed the inland sector of the LVF as the most likely source of this event, and with Molins-Vigatà et al. (2022), who proposed that the LLAF might have also ruptured during the earthquake. As stated in section 3 of this work, this method ranks several candidates and searches for the one most similar to the actual source of the earthquake. This does not mean that the most plausible candidate has to be an exact representation of the earthquake source. Analysis of intensities alone does not provide that level of accuracy; intensity describes the effects of the earthquake, not a direct physical parameter, despite its rather good correlation with ground motion. Unfortunately, because of the intensive agricultural activity and greenhouse constructions, it is currently impossible to carry out paleoseismic studies to accurately identify the source of pre-instrumental earthquakes in Campo de Dalías. This leaves researchers with the analysis of historical earthquakes’ effects (such the one presented in this work) as the only available tool for searching for pre-instrumental seismic sources at the present day in this area. For this reason, and considering how our results agree with Espinar Moreno (1994), Murphy Corella (2019) and Molins-Vigatà et al. (2022), until a better candidate source for the Dalías earthquake is found, we recommend taking into account our results in future characterizations of the LVF and the LLAF, as well as in future seismic hazard assessments of the Campo de Dalías area. In doing so, future researchers should remember that the combined rupture of the LVF and the LLAF seems to be the best rupture candidate, but the single rupture of the LVF is also a plausible source for the Dalías earthquake.

It is important to point out that the 83 available intensity points with which we have worked in the Dalías earthquake are significantly less than the 170 points used during the development of the proposed methodology (de Pro-Díaz et al., 2022). However, despite the difference in observed intensity data points, the same methodology has allowed us to discern the best out of five candidates in the case of the Dalías earthquake. Accordingly, this new study carried out using a lower number of intensity data points highlights the potential of this methodology even with a limited dataset.

The 1804 Dalías earthquake intensity points data set contains 47 sites, out of 83, with the intensity assigned from geological effects (Murphy Corella, 2019). Including these sites in the dataset resulted in having a richer and more complete spatial coverage of intensity data, especially in rural and low-populated areas. According to Guerrieri and Vittori (2007), the EMS-98 scale for macroseismic intensities and the ESI-07 scale for geological effects can be quite well correlated when working from grade IV to grade IX. This is consistent with Ferrario et al. (2022), who verified the good correlation between the two scales and, also, recommended the combined use of EMS and ESI intensities when working with intensity fields. Since the maximum intensity values for the Dalías earthquake are lower than X and the spatial distribution of intensity data is irregular, we consider that the combined use of ESI and EMS intensities contributes to obtain scenarios which are more consistent with all the available data related to the release of seismic energy generated by the seismic source.

In the case of the Málaga earthquake, the difference in size and strike between the area sources obtained from the Gasperini method and using the I_EMS+ESI dataset (A) or only the I_EMS dataset (B) is evident (Figure 9). This shows, once again, how the Gasperini method is biased towards regions with a higher concentration of data points. It is important to remark that none of the candidate ruptures generated intensities high enough to match the observed intensity field, not even the ones with a larger rupture area than the boxer. That could be a sign that the Gasperini method is underestimating the size of the rupture for this earthquake, as it happened with the Dalías case. This underestimation could once again be due to the asymmetry of the intensity field, since the Málaga earthquake also occurred near the coast and we find the same lack of intensity data in the marine area (Figure 9). The lesser amount of available intensity data points in the Málaga case (35 points) might have contributed to underestimating the magnitude as well. Lastly, there is another factor that may have contributed to the underestimation, particularly in the Gasperini method when considering just the macroseismic information based exclusively on the EMS-98 scale (Goded et al., 2008). It has been pointed out that the EMS-98 scale might present some limitations when it comes to studying past earthquakes (e.g., Serva et al., 2015; Silva Barroso et al., 2019). In the case of the area source calculated from the combination of the macroseismic and geological effects data, this issue might have been mitigated by the combination of EMS-98 intensities together with geological and archaeological effects (Silva Barroso et al., 2019). This difference in the source of the input data could also contribute to the difference in size between the two area sources computed using the Gasperini method.

Our methodology identified two plausible candidates to be the source of the Málaga earthquake: a full rupture of the MMF (including both the upper and lower flats, as well as the ramp) and a combined rupture of the full MMF and Cártama, VDG and Mijas systems. The other candidates did not generate intensities high enough to match the observed intensity field. Even scenarios M6 and M7 seem to underestimate slightly the intensities, although by an average of only 0.5 and 0.3 degrees, respectively. Analyzing Figure 12, we can see that the trends we observed in both the Lorca and the Dalías cases, where PGA-based scenarios overall perform better and show lower standard deviations than PGV-based ones, are not reproduced in the Málaga earthquake. The best performing scenarios were built considering Mw 7.3 and 7.4, both higher than the magnitudes previously proposed by other authors for this earthquake, and also higher than the Mw computed by the Gasperini method. We were unable to apply the full extent of the methodology to the Málaga earthquake, so we are aware our results for this event are not as robust as the ones for the Dalías earthquake. Even so, we could rule out five of the seven candidate ruptures. Although we could not discern the best candidate between ruptures F and G, we believe the methodology allowed us to narrow down the possible earthquake source candidates.

We find three issues that might have influenced the results obtained for the Málaga earthquake. The first one is related to the limited amount of available data points. While in the case of the Dalías earthquake we had over 80 observed intensity data points available, in Málaga we only have 35 intensity data points, which are also distributed over a wide area (Figure 9). Even though we consider this an issue, we also have to state that we have been able to point to two possible rupture candidate over the seven proposed ones, which may imply an improvement over the previous knowledge.

The second potential issue is that this is the first time the seismic scenario method is applied to a system of blind thrust faults with complex geometry. Up to this day, the methodology has been used to study earthquakes that could be related to faults with a simpler geometry than that of the MMF and the Cártama, VDG, and Mijas Fault systems. Complex fault geometry in depth might be another limitation, which should be explored in future works. It is also a fact that OpenQuake does not yet compute scenarios with more than one complex rupture surface, so we could not properly model rupture G in this software. Despite this modelling limitation, we use OpenQuake because it is a free, open, user-friendly software that has proven useful before (de Pro-Díaz et al., 2022). The complex geometry issue may also be related to the geology of the Málaga area, where the Maláguide and Alpujárride complexes are predominant. These units have been described to have a complex internal structure (Vera et al., 2004; Insua Arévalo, 2008) and it has been observed before how the local geology and inherited tectonic and lithostratigraphic structures can influence the geometry of an earthquake rupture (e.g., Barchi et al., 2021; Yue et al., 2005; Chiarabba et al., 1997; Zhao and Kanamori, 1995).

The third potential issue is the lack of data about the faults in the study area. We have been working with the models proposed in Insua Arévalo (2008), which were built based on gravimetry, surface geology and geomorphological analysis. As far as we are aware, no public seismic reflection data has been acquired in this area, so the geometry and kinematics of these models have yet to be confirmed. A better constraint is needed in the geometry of the MMF and the Cártama, VDG, and Mijas systems with further geophysical studies, which were not the objective of this work.

Our results point to either a full rupture of the MMF or a conjunct rupture of the MMF and the Cártama, VDG and Mijas Fault system as the best candidate among the considered ones to be the source of the Málaga earthquake. Still, until a better candidate is found and based on our results, we tentatively consider the conjunct rupture as the closest candidate to the source of the Málaga earthquake. Nevertheless, we are aware that more constrains are needed on these results before they can be considered in future seismic hazard assessments of the Málaga area.