Quick Report of the ML=3.3 on 1 January 2023 Guidonia (Rome, Italy) Earthquake: Evidence of a Seismic Acceleration

1. Introduction

On 1 January 2023 at 13:07:46 Universal Time - UT (14:07 Italian time) an earthquake of magnitude ML=3.3 was localized by INGV North of the Capital of Italy, Rome. In particular the closest town is Guidonia Montecelio, and the epicenter was 3 km South-East from this town at the coordinates of 41.982° N and 12.750° E. The depth of earthquake was shallow with an hypocentral depth estimated equal to 9 km (http://cnt.rm.ingv.it/event/33781261, last access 2 January 2023). The event is relatively slow magnitude, but it got particular interest because it follows another seismic event with a slightly slower magnitude ML=3.1 occurred on 23 December 2022 at 16:35:13 UT (17:35 Italian time) with almost the same epicenter (41.984° N and 12.756° E, distance from the above cited event of about 0.5 km). The hypocentral depth of this event was estimated equal to 10 km (http://terremoti.ingv.it/event/33714971, last access 2 January 2023). Based on the data to INGV website at least 681 people (last update 1 January 2023 at 19:31 Italian time) felt the event on 23 December 2022 filling the form on https://www.hsit.it/ (last access 2 January 2023) and 765 people for the event on 1 January 2023 (last update 2 January 2023 at 3:55 Italian time). Previous earthquakes, around this event, were characterized by different focal mechanisms: toward North-East some thrust earthquakes have been historically recorded, while South-East normal fault and West side strike-slip. In particular, the Guidonia earthquakes that occurred on 23 December 2022 and 1 January 2023 could have been produced by a strike-slip fault proposed by Faccenna et al. [1] starting on Cornicolani Mounts and ending towards Colli Albani volcano district crossing Guidonia city (dashed line in Fig. 1 of their paper).

One basic question is if the events are part of a seismic sequence or a seismic swarm. The seismic sequence is characterized by a higher magnitude event (the mainshock) that is followed by aftershocks and could be preceded by foreshocks both of lower magnitude, while the seismic swarm is a chain of earthquakes of comparable and low magnitude without a clear mainshock [2]. If this is a seismic sequence it’s crucial to know if the event of magnitude ML=3.3 of 1 January 2023 is the mainshock or if a future event of higher magnitude is likely to occur. Furthermore, for seismic sequences the decay of the aftershocks generally follows the Omori-Utsu law [3], contrarywise, a seismic swarm has an almost constant seismic rate. The peculiarities of the seismic earthquakes of 22 December 2022 and 1 January 2023 localised close to Guidonia (Rome, Italy) motivates a further investigation that we report in this short communication, based on our previous experience of the lithosphere, atmosphere and ionosphere analysis prior to seismic events.

2. Materials and Methods

In this communication we investigated the earthquake catalogue, the atmospheric data from climatological archive and the Swarm satellites magnetic data.

2.1. Earthquake catalogue

The earthquake catalogue has been retrieved from INGV ISIDE [4] earthquake web portal (http://cnt.rm.ingv.it/, last access 2 January 2023) selecting the events inside the square box delimited by 41.5°N ≤ latitude ≤ 42.5°N and 12.2°E ≤ longitude ≤ 13.3°E from 1^st January 2018 until 1^st January 2023. No constrains have been selected for magnitude and depth, so all the events have been retrieved and they are 1253 earthquakes and their geographical localization is shown in Figure 1a. In the North-East corner of the map is visible some seismicity associate with the very long Italian Seismic Sequence of Amatrice-Norcia-Capitigliano (Montereale) 2016-2017 [5].

The completeness magnitude “Mc” has been analysed by the ZMap software [6] excluding the events of the above-cited Italian Seismic sequence (i.e., the event which are far more than 50 km from Guidonia 2023 earthquake) to avoid a bias in the analysis. In fact, the Apennines central region of Italy interested by the sequence could be better covered not only with permanent seismic stations but also for temporary stations installed after the start of the same sequence. Completeness magnitude was found equal to Mc = 1.4 as visible in graph in Figure 1b. In the following analyses only the earthquakes with magnitude M ≥ Mc will be taken into account. The cumulative Benioff stress has been computed for the events inside the Dobrovoslky radius. Benioff stress in an indication of the stress that is accumulating on a fault [7]. The Dobrovoslky radius is the distance typically used to estimate the preparation area of an earthquake [8].

In particular, Dobrovolsky radius in kilometres is

$$r_{Dob}={10}^{\left(0.43·M\right)}$$

The Benioff cumulative stress S(t) is

$$S\left(t\right)=\sum _{i=1}^{N\left(t\right)}\sqrt{{10}^{\left(1.5·M_i+4.8\right)}}$$

where “i” is the th-earthquake with magnitude M_i and N(t) is the total number of earthquakes occurred until the time t. According to Mignan et al. [9] and Cianchini et al. [10] the following time-to-failure power law has been fitted to the Benioff cumulative stress:

$$S\left(t\right)=A+B·{\left(t_f-t\right)}^m$$

where A and B parameters are proportional to the initial and final accumulated stress and m is an exponential factor in the range from 0 to 1 and it is typically equal to 0.3 [10]. In this work, m has been fixed to m=3 to guaranty stability of the power law fit. The graphs report t_f which is the time to failure, i.e., the likelihood estimation of the mainshock time.

2.2. Atmospheric data processing

To investigate the atmospheric proprieties the data from the climatological archive MERRA-2 by NASA [11] have been investigated. With the experience of previous investigated earthquake, the MEANS software has been used with the remotion of a linear trend to take into account the possible effect of the global warming as done in De Santis et al. [12] for Ridgecrest, US, 2019 earthquake analysis. The MEANS algorithm has been used the first time by Piscini et al. [13] to investigate volcano eruptions and after applied to several other volcano eruptions and earthquakes [12,14,15,16,17,18]. The algorithm essentially estimates the typical value of the atmospheric parameter for the day and region constructing an historical time series and in particular the mean and standard deviation of the parameter for the area under study. If the parameter in the year with the natural hazard event overpasses two standard deviations of the historical time series is defined as anomalous. In this report we investigated surface air temperature, aerosol and Sulphur dioxide.

2.3. Ionospheric data processing

The ionosphere has been monitored by the European Space Agency (ESA) Swarm constellation. Swarm mission is composed by three identical satellites in low Earth quasi-polar orbits, called Alpha, Bravo and Charlie aiming to measure the geomagnetic field with the best precision available at the state of art [19]. The use of Swarm data to study the preparation of the earthquakes has been explored by several researches in the last years using a single case study approach [18,20,21,22,23,24,25,26,27,28,29,30,31] or even statistically correlating the magnetic and electron density anomalies of Swarm with M5.5+ earthquakes occurred in the first 4.7 years of Swarm mission by De Santis et al. [32], the first 8 years by Marchetti et al. [33] or using Machine Learning by Xiong et al. [34].

In this report we applied the MASS algorithm to Swarm magnetic field data defined with all details in [27,31,32,33]. It analyses the data track by track removing the background with estimating the derivative and subtracting a cubic spline. The residuals are monitored by a moving window and comparing the root mean square of the window with the one of the entire track among −50°S and +50°N of geomagnetic latitude to avoid the influence of polar activity.

3. Results

Here we illustrate the results for seismological atmospheric and ionospheric investigations in respectively sections.

3.1. Seismological investigation

The cumulative Benioff stress has been calculated inside the Dobrovoslky area testing several mainshocks from magnitude 3.3 until 5.5 as listed in Table 1. The lower limit corresponds to the assumption that the event of 1 January 2023 would be the mainshock and the higher magnitude events correspond to the hypothesis that a higher magnitude event is incoming. In case the Dobrovolsky radius overpassed the area represented in Figure 1a, the earthquake with M ≥ Mc have been further retrieved.

For each fit the adjusted R² and the estimated time of failure is reported in Figure 2, Figure 3 and Figure 4 for 5 years, 1 year and 6 months before the ML=3.3, 1 January 2023 Guidonia earthquake, respectively.

In addition the acceleration coefficient C defined by Cianchini et al. [10] has been reported calculating as:

$$C=\frac{{rms}_{powerlaw}}{{rms}_{linearfit}}$$

As much as low is “C” as much the fit represents an acceleration, on the other hand C = 1 indicates that the trend is linear or C > 1 that linear fit is better than power law. Such situation could be due to a trend which is influenced by aftershock of another seismic sequence as the one of Amatrice-Norcia-Capitigliano (Montereale) Central Italy 2016-2017, in particular for long trend and larger radius of analysis.

The investigations in 5 years and 1 year before ML=3.3 Guidonia 1^st January 2023 earthquake don’t provide interesting results. In fact, the trends are dominated by other important seismicity which occurred in the area (for example Mw=3.3 11 April 2022, Ciciliano, Rome event http://cnt.rm.ingv.it/event/30543871, last access 2 January 2023). The tests with larger magnitude are instead dominated by the residual sequence of Apennine cited above. This is the reason why higher magnitude than 4.7 are not reliable to be tested as the crossing with seismicity from Central Apennines region becomes predominant.

On the other hand, the investigation in the 6 months before ML=3.3 Guidonia 1^st January 2023 earthquake reported in Figure 4, provide interesting results. From magnitude 3.3 until 3.9 the estimated time-to-failure is 0.0 days and the maximum acceleration is reported for M=3.3 with C=0.318. Considering that no higher events have been recorded in this area at 1 January 2023, among these 4 analyses the most reliable need to be considered the test in M=3.3, i.e., the hypothesis that the past earthquake was the mainshock. Despite this, the analysis for M=4.1 and M=4.3 provided good fits with a time to failure in few days (M4.1 on 5 January 2023 or M4.3 on 25 February 2023). The investigations of M4.5 and 4.7 provided very goodness of the fits even though they are almost linear (as underlined by their C close to 1). So, selecting the analysis with higher adjusted R² and $\mathrm{C}\ll 1$ the most reliable result is the analysis of M4.1 with a time to failure in 3.5 days after Guidonia 1^st January 2023 earthquake with R²-adj = 0.988 and C = 0.356. Even though this could appear as a prediction, that interpretation is not totally correct, in fact, firstly it’s necessary to estimate the maximum magnitude that the interested fault could generate, and secondary determine if the accumulated stress has been totally or partially released by these events. Furthermore, the time-to-failure power law could not take into account several geological and tectonic factors, so it cannot be used “alone”. Finally, the right interpretation for the authors is that in the past six months in the circle area of radius 58km centered on Guidonia 1 January 2023 earthquake has been recorded a clear seismic acceleration, but further studies are necessary to propose eventual future scenarios.

3.2. Atmospheric investigations

In this section we presented atmospheric investigation in the six months before Guidonia 1^st January 2023 earthquake. The time series ends on 30 November due to December data still not yet available. The results of surface air temperature, aerosol and SO₂ are provided in Supplementary Materials in Figure S1, S2 and S3 respectively. All the time series were calculated in a square box centered on the Guidonia 1^st January 2023 earthquake with 3° side. Surface air temperature have some anomalous values but no-one of them have a persistence which is a criteria that was found effective for Central Italy by Piscini et al. [35]. Nevertheless, aerosol and especially SO2 provided several anomalous days and we provided here their geographical map to investigate better if the increase of these quantities is local or global. The aerosol on 18 August 2022 is particularly anomalous and its map is shown in Figure 5. From the distribution of the aerosol, it seems that all the Tyrrhenian Sea recorded high value of aerosol, probably due to some special weather conditions and in any case, it doesn’t seem related to the earthquake.

Figure 6 reported 3 anomalous days of SO₂ in the investigated area, i.e., the 18 July, 23 October and 1 November 2022. The 18 July and 1 November show a slightly higher value of SO₂ almost on all the coastline, excluding a local phenomenon. Despite this in all the 3 maps the highest values are recorded in the sea in correspondence of an important subduction slab of Tyrrhenian Sea (dashed grey line). The mechanism of emission of this SO₂ could be due to the stress on the subduction place, but also to the important volcanic activity for the several volcanos close to this area (Aeolian Islands with Stromboli, submersed Marsili, or Etna slightly southern in Sicily). These volcanos are not only well known to emit SO₂ but also they contributed significantly to the global budget of SO₂ emissions [36].

3.3. Ionospheric Investigations

The ionosphere has been investigated under the most reliable seismological results of an earthquake of M4.1 and so with the corresponding Dobrovolsky’s area. MASS algorithm has been run for the data from 3 July 2022 until 28 December 2022 (last available data with standard Swarm processor of MAGX_LR data product).

Considering the very low dimension of the research area, Swarm satellites do not cross all days within the area, reducing the chances to detect possible electromagnetic signals. After processing all the tracks of Swarm Alpha, Bravo and Charlie, only 4 tracks present some anomalies inside the research area in geomagnetic quiet conditions (|Dst| ≤ 20 nT and ap ≤ 10 nT), using a window of 1° latitude (found optimal in Central Italy Swarm investigation [37]) and a threshold of rms inside the window to be more than 2.5 times the one of the whole track. Table 2 reports the list of anomalous tracks. No anomalies have been identified in vertical Z component or in absolute scalar intensity of geomagnetic field with the selected parameters.

The most two interesting tracks in Y-East component of geomagnetic field are reported in Figure 7 and 8 and the other two anomalous tracks in Figure S4 and S5 of Supplementary Materials. In particular, Figure 7 shows an anomalous signal of about 3 nT/s peak-peak intensity, but the geomagnetic index Dst is at the upper boundary limit of 20 nT. In fact, looking at the trend of Dst (e.g., on https://wdc.kugi.kyoto-u.ac.jp/dst_realtime/202208/index.html, last access 2 January 2023) it continuous to increase in the following hours and after decreasing reaching −58 nT on 8 August 2022, so this track was acquired in the Sudden Commencement of a very moderate geomagnetic storm. Despite the signal is perfectly center on the incoming earthquake the external source is the more reliable one.

Figure 8 shows instead a track of Swarm Bravo acquired during very quiet geomagnetic conditions (Dst = −1 nT, ap = 2 nT) with an anomalous geomagnetic field Y-East component signal recorded above Central Italy. Despite the track preceded the first shock of ML3.1 of 23 December 2022 of one week and during the seismic acceleration phase the diffuse nature of the anomaly seems difficult to confirm it could be seismo-induced even it cannot be excluded.

4. Discussion and conclusions

In this Communication paper we investigated the ML=3.3 Guidonia (Rome, Italy) earthquake occurred on 1 January 2023 by seismic catalogue, atmospheric and ionospheric data. From a preliminary investigation, a seismological analysis identified an acceleration of the seismicity around such epicenter especially focusing in the last six months. The atmospheric and ionospheric investigations, on the other hand, don’t provide any results which look as seismo-induced nature. The absence of possible coupling of the lithosphere with atmosphere and ionosphere is in agreement with the low magnitude of the seismic event (ML=3.3).

Finally, the seismological investigation suggest that these events are likely part of a seismic sequence and the analysis are compatible with 2 main hypotheses:

• The ML3.3 of 1 January 2023 is the mainshock of the seismic quiescence (R²-adj = 0.977 and acceleration coefficient C = 0.32)
• The mainshock could be an incoming event of magnitude M4.1 (R²-adj = 0.988 and acceleration coefficient C = 0.36).

The second hypothesis seems mathematically slightly more probable, but both the hypotheses describe well the seismicity recorded in the last months around Guidonia (Rome, Italy). In any case, the acceleration nature (C << 1) of all the reliable solutions tend to exclude that the events are part of a seismic swarm.

Even though some speculation on the evolution of the seismic activity can be done, they cannot be considered as a prediction, because no tests have been ever made, so without a validation of the method by an independent organism (e.g., CSEP, https://cseptesting.org/, last access 2 January 2023) the analysis has no prediction value. Furthermore, the correct interpretation of the time-to-failure power law fit is that in the last six months in a circle area of radius of about 58 km has been detected a seismic acceleration but to propose future evolution of the seismic sequence further investigations are absolutely necessary. Due to the peculiarity observed in the past months it could be evaluated the installation of geochemical multi-parametric station to detect eventual gases or fluids released from the ground.

Finally, this example, shows how these techniques could be applied in real time for monitoring (not yet prediction) purpose of lithosphere, atmosphere and ionosphere and the help of remote sensing data to understand and monitor our vulnerable planet.