The Preparation Phase of the 2022 M L 5.7 Offshore Fano (Italy) Earthquake: A Multiparametric–Multilayer Approach

: This paper presents an analysis of anomalies detected during the preparatory phase of the 9 November 2022 M L = 5.7 earthquake, occurring approximately 30 km off the coast of the Marche region in the Adriatic Sea (Italy). It was the largest earthquake in Italy in the last 5 years. According to lithosphere–atmosphere–ionosphere coupling (LAIC) models, such earthquake could induce anomalies in various observable variables, from the Earth’s surface to the ionosphere. Therefore, a multiparametric and multilayer approach based on ground and satellite data collected in each geolayer was adopted. This included the revised accelerated moment release method, the identification of anomalies in atmospheric parameters, such as Skin Temperature and Outgoing Longwave Radiation, and ionospheric signals, such as Es and F2 layer parameters from ionosonde measurements, magnetic field from Swarm satellites, and energetic electron precipitations from NOAA satellites. Several anomalies were detected in the days preceding the earthquake, revealing that their cumulative occurrence follows an exponential trend from the ground, progressing towards the upper atmosphere and the ionosphere. This progression of anomalies through different geolayers cannot simply be attributed to chance and is likely associated with the preparation phase of this earthquake, supporting the LAIC approach.


Introduction
The study of a possible coupling mechanism between the lithosphere, atmosphere, and ionosphere, also known as lithosphere-atmosphere-ionosphere coupling (LAIC), before an earthquake (EQ) or a volcanic eruption, is becoming increasingly relevant within the scientific community.Prior to the occurrence of such geophysical events, the Earth emits transient signals, sometimes strong but more often subtle and fleeting, which can manifest as local variations in the magnetic field, electromagnetic emissions across a wide range of frequencies (mostly 0.001 Hz-100 KHz, i.e., ULF-ELF-VLF), and a variety of atmospheric and ionospheric phenomena.This is a rather complex mechanism and there is considerable uncertainty about the nature of the processes that could produce these signals, both within the Earth's crust and on its surface.Over the years, various models have been developed to suggest a connection, i.e., coupling, among the geolayers following different channels.
Hayakawa [1] proposed three channels to establish connections between different observations in each layer: the chemical channel (also known as the electric field channel), the acoustic gravity wave (AGW) channel, and the electromagnetic (EM) channel.Freund [2] proposed an electrostatic channel in which stressed rocks release positive charge carriers known as "positive holes".Dahlgren et al. [3] contested this theory, in particular studies [4][5][6].Scoville et al. [7] pointed out some pitfalls in the experiment.
Pulinets and Ouzounov [8] proposed a radon release from the crust in the earthquake preparation zone; radon is the radioactive gas produced in the decay chain of uranium or thorium emitted from the ground, affecting the electric field in the troposphere-ionosphere electric circuit.Surkov et al. [9] further investigated the potential impact of radon emissions on the atmosphere and ionosphere, finding that localized changes in atmospheric currents due to radon have minimal effects on the electron distribution in the ionosphere.Similarly, the study proposed by Schekotov et al. [10] does not find clear connections between electromagnetic variations and changes in pre-seismic temperatures, casting doubt on the hypothesis that radon emissions influence the ionosphere.
To explain seismo-ionospheric effects, other authors [11][12][13] suggested atmospheric processes producing acoustic and/or gravity waves in the seismic preparation region.Investigating LAIC processes, the multiparametric analysis around earthquakes has become a highly debated topic, involving various parameters across layers and utilizing different observation technologies within the same layer.A comprehensive study [14] on investigating the LAIC mechanism of the Mw = 7.2 Haiti earthquake on 14 August 2021, considered 52 precursors, including GPS Total Electron Content (TEC), 4 from CSES-01 satellite data, 7 lithospheric and atmospheric precursors from AIRS and OMI sensors, and 40 from the Swarm satellite constellation.This author observed a significant amount of anomalous values, suggesting a sequence possibly due to ion radiation from the Earth, i.e., a thin layer of particles transferring the electric field to the upper atmosphere and then to the ionosphere.Another study [15] examined two major channels of LAIC mechanisms, acoustic and electromagnetic, to analyze pre-seismic irregularities of the 2020 Samos (Greece) earthquake of M = 6.9.This study considered TEC, AGW, bursts of energetic particles in the radiation belt, magnetic field, electron density, and temperature, obtaining significant anomalies from 10 to 1 day before the seismic event.Pre-seismic low values of TEC were observed in regions with lower b-values in a study focused on the analysis [16] of the Mw = 7.7 Colima (Mexico) earthquake on 19 September 2022, indicating a higher probability of larger earthquakes.
Several studies focused on the correlation between the ionospheric TEC anomalies and the earthquakes [17][18][19][20].Specifically, to explain the coupling of TEC anomalies and seismic events [18], they confirmed the existence of an anomalous electric field in seismogenic zones triggered by stressed rocks in earthquake regions associated with fault lineaments, while Tachema [20] proposed that the space between the lithosphere and ionosphere is occupied by a coherent structure of electrons and protons, transmitting electromagnetic waves generated during the seismic nucleation of rocks at depth.
De Santis et al. [21] analyzed several parameters from the lithosphere, atmosphere, and ionosphere on the occasion of the 2019 M = 7.1 Ridgecrest EQ, finding an accelerated progression of the cumulative number of all anomalies as the mainshock was approaching.Wang et al. [22] detected anomalous changes in the lithosphere, atmosphere, and ionosphere near the epicenter before the 2021 M = 7.4 Madoi EQ; meanwhile, Akhoondzadeh and Marchetti [23] analyzed the behavior of more than 50 different lithosphere-atmosphereionosphere anomalies during the preparation phase of the 2023 Turkey EQ, identifying a progressive increase in the number of anomalies starting about 10 days before, with the major peak the day before the mainshock.
On 9 November 2022, at 06:07:25 UTC, an M L = 5.7 (Mw = 5.5) EQ occurred approximately 30 km offshore the Marche coast in the Adriatic Sea, at latitude 43.984 • N, longitude 13.324 • E, and 5 km of depth.Just a minute after the mainshock, a strong aftershock of M L = 5.2 occurred about 8 km south of the main event, at latitude 43.913 • N, longitude 13.345 • E, and 8 km of depth.The two major shocks activated a seismic sequence of about 400 aftershocks lasting a week, thirteen of them with M L = 3.5.The seismic sequence oc-curred in correspondence with the frontal fault systems of the Northern Apennines, where ongoing convergence is accommodated on a series of buried faults, still poorly understood.The moment tensor solution of the mainshock indicates a reverse mechanism on a NW-SE trending fault plane, as well as the moment tensor solutions of the M L ≥ 3.5 events of the sequence.However, no moment tensor solution has been computed for the M L = 5.2 event due to the overlap and interference of phases from the two events [24].
According to Pezzo et al. [24], the slip occurred along a thrust fault dipping approximately 24 • SSW over a length of about 15 km, consistent with seismic reflection data propagating downward from the mainshock hypocenter, confirming the ongoing seismotectonic activity of this sector of the Apennines that is still propagating towards the foreland, approximately in a piggy-back thrust sequence.The area of greatest impact is along the coastal stretch between Fano and Ancona, where a maximum intensity of 5 European Macroseismic Scale (EMS-98) [25] has been estimated.Within this zone, very slight sporadic damages have been observed.Particularly, the district outskirts of Ancona have reported damages to recent reinforced concrete buildings, likely due to local amplification effects.The estimated maximum intensity reached 5 EMS-98 in some locations; however, macroseismic effects rapidly decreased to 4 EMS-98 inland, at a short distance from the coast [26].
This study aims to provide a comprehensive view of the effects observed on various geological and atmospheric layers within the framework of lithosphere-atmosphereionosphere coupling (LAIC) models, during the preparatory phase of the seismic event recorded offshore the coast of Marche on 9 November 2022.Previous studies related to this seismic event [24,[27][28][29][30][31] have predominantly focused on its tectonic and seismological aspects, while we adopted the LAIC approach as a novelty for this EQ.Additionally, the characteristics of this seismic event, such as the magnitude slightly below 6 and its occurrence at sea, piqued our interest in investigating whether this approach could yield fruitful results.In particular, the presence of the sea conductive layer could limit the occurrence of one kind of LAIC instead of another.
Figure 1 describes the parameters studied across different geological layers from bottom to top, looking for anomalous perturbations potentially associated with the preparation phase of the M L = 5.7 EQ.Specifically, seismicity acceleration preceding the mainshock is investigated in the lithosphere using the revised accelerated moment release (R-AMR) method.Moving upwards, the atmosphere is studied through parameters such as Skin Temperature (SKT) and Outgoing Longwave Radiation (OLR).Lastly, ionospheric analysis includes data from the Rome AIS-INGV ionosonde, the European Space Agency (ESA) Swarm satellite mission, and the National Oceanic and Atmospheric Administration (NOAA) satellites.

Seismotectonic Settings
The Adriatic Sea appears as a result of a good variety of structural and stratigraphic processes (Figure 2) guided by fault-related anticlines formed in the Plio-Miocene connected to the main Apennine thrust chain, deeper carbonate structures developed in In the next section, the data used will be introduced, followed by a presentation of the applied methods alongside the main results.Finally, the work will conclude with a discussion and a conclusion.
Given the proximity of Fano town to the epicenter, i.e., approximately 29 km, this event will be simply referred to as the "Fano EQ".

Seismotectonic Settings
The Adriatic Sea appears as a result of a good variety of structural and stratigraphic processes (Figure 2) guided by fault-related anticlines formed in the Plio-Miocene connected to the main Apennine thrust chain, deeper carbonate structures developed in the south, and a very shallow structure in the Late Pliocene to Quaternary in the central area [32].During the Mesozoic, this area was affected by an extensional tectonic phase in the Middle Liassic and a compressional paleoinversion in the Lowermost Cretaceous [33,34].The development of the Alps and the Apennines started from the Middle Eocene onwards in the African continental margin [35][36][37][38][39]. Subsequently, a flexure of the lithosphere belonging to the Adria margin concerned the most internal areas and migrated eastward through time, forming foredeep basins oriented sub-parallel to the belts.The Adriatic domain corresponds to the youngest part of the belt, strictly connected to the evolution of the Apennine fold and thrust belt and to the interaction with the Dinarides, sub-parallel orogenic belts with opposing vergences [32].In detail, this thrust front is buried beneath Early Pliocene-Quaternary synorogenic deposits.

Seismotectonic Settings
The Adriatic Sea appears as a result of a good variety of structural and stratigraphic processes (Figure 2) guided by fault-related anticlines formed in the Plio-Miocene connected to the main Apennine thrust chain, deeper carbonate structures developed in the south, and a very shallow structure in the Late Pliocene to Quaternary in the central area [32].During the Mesozoic, this area was affected by an extensional tectonic phase in the Middle Liassic and a compressional paleoinversion in the Lowermost Cretaceous [33,34].The development of the Alps and the Apennines started from the Middle Eocene onwards in the African continental margin [35][36][37][38][39]. Subsequently, a flexure of the lithosphere belonging to the Adria margin concerned the most internal areas and migrated eastward through time, forming foredeep basins oriented sub-parallel to the belts.The Adriatic domain corresponds to the youngest part of the belt, strictly connected to the evolution of the Apennine fold and thrust belt and to the interaction with the Dinarides, sub-parallel orogenic belts with opposing vergences [32].In detail, this thrust front is buried beneath Early Pliocene-Quaternary synorogenic deposits.The entire seismic sequence of November 2022 unfolds along the outermost structure of the Apennine orogeny, characterized by a series of NW-SE trending, NE verging folds forming the easternmost edge of the Apennine thrust front [31,[41][42][43].Over the Tertiary-Quaternary period, this front has gradually migrated towards the east-northeast (e.g., [44]).Geological evidence suggests the ongoing growth of these folds, indicating the continued activity of blind thrust fronts (e.g., [41]).
The presence of historical and instrumental earthquakes with Mw ≥ 5.5, as shown in Figure 3, suggests that thrust faults are also seismogenic [30].
area, which is experiencing a similar but opposite (southwestward) movement.The two hypocenters belong to the blind thrust fault system running parallel to the Marche Coast.This compressive front is located approximately 25-35 km offshore, with a length of about 70 km, and it is referred to as ITCS10 in the Database of Individual Seismogenic Sources (DISS) [41].This thrust system has therefore been identified as responsible for this seismic sequence (see Figure 3).; the first event is marked with a yellow star and the second event is shown with a green star.Historical and instrumental earthquakes from CPTI15 [45] are indicated with colored squares, with earthquakes of Mw ≥ 5.5 highlighted in red.The surface projections of seismogenic zones are depicted with orange ribbons [41].The focal mechanisms of the 9 November 2022 earthquake and the event of 30 October 1930, represented by the grey and white balls, come from TDMT (Time Domain Moment Tensor) and Vannoli et al. [46], respectively (modified from [30]).

Data and Methods
To study the effects of LAIC, several datasets are required since each geolayer being investigated demands specific data from various sources.As the analysis is conducted separately in each layer, there are different resolutions both in time and in space, which will be described below.[46], respectively (modified from [30]).
The events recorded on 9 November 2022 represent a manifestation of the ongoing contraction between the Apennine chain, moving towards the northeast, and the Balkan area, which is experiencing a similar but opposite (southwestward) movement.The two hypocenters belong to the blind thrust fault system running parallel to the Marche Coast.This compressive front is located approximately 25-35 km offshore, with a length of about 70 km, and it is referred to as ITCS10 in the Database of Individual Seismogenic Sources (DISS) [41].This thrust system has therefore been identified as responsible for this seismic sequence (see Figure 3).

Data and Methods
To study the effects of LAIC, several datasets are required since each geolayer being investigated demands specific data from various sources.As the analysis is conducted separately in each layer, there are different resolutions both in time and in space, which will be described below.

Lithospheric Data
From a lithospheric point of view, a seismological analysis was carried out to characterize the seismicity of the area affected by the M L = 5.7 EQ under inspection, paying particular attention to the area mainly affected by the imminent seismic event.The seismic data provided by the Italian INGV Catalog [47] were used by selecting a circular area with a radius of 150 km around the epicenter and imposing a depth limit of 100 km, covering the period from 1 January 2012 to 8 November 2022 (Figure 4).The temporal and spatial reso-lutions of the used catalog in terms of earthquake detection (with associated information) are of the order of seconds (or even fractions of a second) and a few km, respectively.
From a lithospheric point of view, a seismological analysis was carried out to characterize the seismicity of the area affected by the ML = 5.7 EQ under inspection, paying particular attention to the area mainly affected by the imminent seismic event.The seismic data provided by the Italian INGV Catalog [47] were used by selecting a circular area with a radius of 150 km around the epicenter and imposing a depth limit of 100 km, covering the period from 1 January 2012 to 8 November 2022 (Figure 4).The temporal and spatial resolutions of the used catalog in terms of earthquake detection (with associated information) are of the order of seconds (or even fractions of a second) and a few km, respectively.To evaluate R-AMR [48,49], the magnitude completeness (Mc) of the catalog, which represents the minimum value of magnitude for detection, was estimated using the method of maximum curvature [50].The estimation of Mc was performed as a function of time by sliding time windows, each containing 150 EQs and stepping by five events.The uneven distribution of events within the selected circular area, due to the lack of seismometers on the Adriatic seabed, was taken into consideration.However, the presence of a seismic station located in Banja Luka (Bosnia and Herzegovina), named BLY, allowed for precise event localization.Based on these considerations, an Mc = 2.2 was To evaluate R-AMR [48,49], the magnitude completeness (Mc) of the catalog, which represents the minimum value of magnitude for detection, was estimated using the method of maximum curvature [50].The estimation of Mc was performed as a function of time by sliding time windows, each containing 150 EQs and stepping by five events.The uneven distribution of events within the selected circular area, due to the lack of seismometers on the Adriatic seabed, was taken into consideration.However, the presence of a seismic station located in Banja Luka (Bosnia and Herzegovina), named BLY, allowed for precise event localization.Based on these considerations, an Mc = 2.2 was chosen, resulting in a catalog of 9.099 events.Figure 5 depicts the events (blue and red points) that contributed to the acceleration identified by the R-AMR analysis.
Several studies (e.g., [51,52]) suggested that before significant EQs, there is an accelerated seismic activity under specific conditions.This phenomenon, explained by the Critical Point Theory, likens the main EQ to a phase transition occurring at a "time-to-failure", t f .chosen, resulting in a catalog of 9.099 events.Figure 5 depicts the events (blue and points) that contributed to the acceleration identified by the R-AMR analysis.
Several studies (e.g., [51,52]) suggested that before significant EQs, there is accelerated seismic activity under specific conditions.This phenomenon, explained by Critical Point Theory, likens the main EQ to a phase transition occurring at a "time failure",  .The seismicity preceding the mainshock, often hidden in catalogs, can be revea through methods like Accelerated Moment Release (AMR), particularly its revi version, known as R-AMR [48,49].The R-AMR algorithm was applied to EQ data in Marche region before the Fano EQ.Acceleration in seismicity is measured by examin the accumulation of seismic Benioff strain  =  , where each event releases str proportional to the square root of its energy  = 10 ( . . ) J. The cumulative str () = ∑  , is known as Cumulative Benioff Strain.The regional increase in cumulative Benioff deformation before a large shock is expressed by a power-law ti to-failure functional relation: where  is the time-to-failure (i.e., the occurrence of the mainshock) and m is an inve measure of how quickly the acceleration grows around  .To evaluate the quality seismic acceleration compared to a linear trend representing the background seismic Bowman et al. (1998) [51] introduced the C factor, which is the ratio of the sum of squares of the residuals of the fit of s = s(t) and the same quantity of a linear fit; a C va less than 1 indicates acceleration, with lower values indicating more promin acceleration.De Santis et al. [48] improved the technique by focusing on strain depos The seismicity preceding the mainshock, often hidden in catalogs, can be revealed through methods like Accelerated Moment Release (AMR), particularly its revised version, known as R-AMR [48,49].The R-AMR algorithm was applied to EQ data in the Marche region before the Fano EQ.Acceleration in seismicity is measured by examining the accumulation of seismic Benioff strain s i = √ E i , where each event releases strain proportional to the square root of its energy E i = 10 (1.5M i +4.8)J.The cumulative strain, s(t) = ∑ s i , is known as Cumulative Benioff Strain.The regional increase in the cumulative Benioff deformation before a large shock is expressed by a power-law time-to-failure functional relation: where t f is the time-to-failure (i.e., the occurrence of the mainshock) and m is an inverse measure of how quickly the acceleration grows around t f .To evaluate the quality of seismic acceleration compared to a linear trend representing the background seismicity, Bowman et al. (1998) [51] introduced the C factor, which is the ratio of the sum of the squares of the residuals of the fit of s = s(t) and the same quantity of a linear fit; a C value less than 1 indicates acceleration, with lower values indicating more prominent acceleration.De Santis et al. [48] improved the technique by focusing on strain deposited on the mainshock fault by surrounding seismicity, corrected by a damping function with distance.Cianchini et al. [49] further enhanced the algorithm based on 14 case studies, confirming its ability to reveal hidden acceleration in seismic sequences and provide estimates of t f and expected magnitude based on parameters A and B of the functional relation.In this analysis, the algorithm is automatic: considering the s = s(t) time series backward and excluding the mainshock, it detects the time when the acceleration starts and the minimum (no attenuation) and maximum circles (with some attenuation).Therefore, for the characterization of seismicity, the R-AMR algorithm was applied to the seismic catalog without taking into account any potential heterogeneity.

Atmospheric Data
For the atmospheric analyses, data from the ECMWF (European Center for Mediumrange Weather Forecasts) ERA5 climatological reanalysis dataset were utilized.This dataset provides comprehensive reanalysis from 1940 up to 5 days prior to the current date, assimilating a wide range of observations in the upper atmosphere and near-surface regions.In our study, however, we focused on parameters dating back to 1980 [53].This dataset is known for its consistent coverage in both space and time and is minimally affected by observational conditions such as cloud cover in satellite observations.Nighttime values were specifically considered due to their reduced susceptibility to local meteorological changes.Specifically, we analyzed the parameters of SKT and OLR, which are typically reported to be influenced by impending earthquakes.Several studies [54][55][56] have demonstrated how these parameters are directly affected by the "thermodynamic channel" in LAIC models.The ECMWF time series for each atmospheric parameter were collected and pre-processed to apply the Climatological Analysis for Seismic Precursor Identification (CAPRI) algorithm [54,55].This algorithm compares daily time series of the current year with a historical dataset spanning forty-two years (1980-2021), within a temporal window preceding the seismic event, in our case of 90 days.An anomaly is identified if the observed value persistently exceeds the mean of the historical series by two standard deviations.However, for the current study, anomalies were defined as values exceeding 1.5 standard deviations, considering that the earthquake magnitude was below 6.Additionally, the geographic area investigated was determined based on the circular earthquake preparation region (or Dobrovolsky area) centered on the epicenter [57], resulting in the selection of a geographical area of 2 • latitude and 2 • longitude.

Ionospheric Data
To investigate the ionosphere for potential disturbances associated with the Fano EQ, data from ionosonde and satellite sources were analyzed.From the ground, ionosonde measurements can detect the critical frequency of the F2 layer (foF2), the height (hmF2) of the main electron density peak, and also the information of the sporadic E (Es) layer (such as its height, h'Es; its critical frequency, foEs; and the blanketing frequency, fbEs), which can represent the variations of the corresponding layers.Low Earth Orbit (LEO) satellites can detect the in situ plasma parameters and electromagnetic field.The ionospheric station of Rome (43.98 • N; 13.32 • E) is located 251.44 km from the epicenter, so within the EQ preparation zone according to the formula by Dobrovolsky et al. [57].Hourly data manually scaled by an experienced operator from the ionograms recorded with the Advanced Ionospheric Sounder (AIS-INGV) [58] were used in this study [59].Ionospheric anomalies are defined by significant deviations in the parameters h'Es, fbEs, and foF2 compared to a specified background level determined by 27-day hourly running medians centered on the observation day.Specifically, these deviations are calculated according to the following expressions and must satisfy the following criteria, provided they occur within a few hours under geomagnetically quiet conditions specified by daily geomagnetic index values Ap < 9 nT: where ∆ indicates absolute deviations and δ indicates relative deviations.The satellite data, on the other hand, are collected by the ESA Swarm three-satellite constellation (Alpha, Bravo, and Charlie).Alpha and Charlie orbit side by side at an altitude of around 460 km, while Bravo orbits at 510 km.Satellite magnetic data have been analyzed using the MASS (MAgnetic Swarm anomaly detection by Spline analysis) methodology for the magnetic data method (see, e.g., [60]).This technique is used to detect electromagnetic anomalies from Swarm magnetic field data (Level 1B, low resolution of 1 Hz) from the analysis of the three components of the geomagnetic field (X, Y, and Z) and intensity (F) for every track of each satellite (Swarm A, Swarm B, and Swarm C) recorded over the EQ preparation area [57].The analysis consists of the determination of the first differences of the time series (dX/dt, dY/dt, dZ/dt, and dF/dt) and the removal of the long trend using a cubic spline.Moreover, it is important to evaluate the quality of the data using the quality flags provided by ESA, considering only magnetic quiet times (|Dst| ≤ 20 nT and ap ≤ 10 nT) and excluding polar regions (±50 • geomagnetic latitudes) because they are very disturbed.After these steps, the root mean square (rms) of sliding windows of 7 • shifting every 1.4 • was compared with the whole root mean square of the track (RMS).When rms is greater than kt times the value of RMS, the corresponding window is classified as anomalous.For satellite magnetic data, kt is established as 2.5 [60].This analysis covered a time period ranging from 90 days before the earthquake occurrence to 10 days after and was confined in an area comparable with the EQ preparation region [57].
Since the 1980s, NOAA satellites have continuously circled Earth's orbit, employing a shared instrument for detecting charged particles since 1998.This instrument, the Medium-Energy Proton and Electron Detector (MEPED), is integrated into NOAA satellites and features eight solid-state detectors meticulously crafted to gauge proton and electron counting rates (CRs) within the 30 keV-200 MeV range with a sampling of 2s [61].These measurements provide invaluable insights into various phenomena, including radiation belt populations, energetic solar proton events, and the low-energy segment of the galactic cosmic-ray population.Specifically tailored for electron detection, the two telescopes within MEPED delineate three energy bands ranging from 30 keV to 2.5 MeV.The first telescope, angled at 9 • towards the local zenith, captures one perspective, while the second telescope, positioned orthogonally at 90 • along the satellite's motion, provides a complementary view.Each compact solid-state detector boasts nominal geometric acceptances of 0.1 cm 2 sr and opening angle apertures of ±15 • .
Statistical correlations between strong EQs and electron bursts (EBs) detected by NOAA were mainly observed considering seismic events around the equator in both the West and East Pacific [62,63].It was in agreement with the bouncing points of the inner Van Allen Belt (VAB), points where electrons go down close to the earth's surface, which are also around the equator and cross the ring of fire twice.However, the inner VAB does not generally overhang latitudes above 35 • , while mid-latitude EQs are located below the slot region between the internal and external VABs.Seismic events occurring at mid latitudes have therefore apparently minor probability to interact with trapped particles.The representation of a possible interaction scenario between mid-latitude seismic events and trapped electrons is depicted in Figure 6, where the possible area to observe EBs connected to the Fano EQ is shown in pink color.These areas have a longitude range of a few tenths of degrees around those of the seismic event.In fact, the longitude of the seismic event is close to the South Atlantic Anomaly (SAA) border longitude, and the electrons' mirror points that are far from the SAA longitude will hardly descend to the satellite altitude.Thus, the probability of bouncing electrons to cross the satellite decreases getting away from the SAA.For what concerns the lower latitude EBs, electrons are thought to be coming (direction indicated with pink arrow) from the inner VAB (in yellow), with the interaction (in blue) running along the magnetic field lines.Instead, for what concerns mid-high latitude EBs, electrons are thought to be escaping the trapped conditions (direction indicated with another pink arrow) from the external VAB (still in yellow), with the interaction (still in blue) connecting the two phenomena along the minimum path.The geomagnetic lines (in green) of the internal VAB can cross the lithosphere up to 30 • -35 • in latitude, so being unlikely to be affected by tectonic activity, whereas the radial propagation of LAIC (in blue) should be able to intercept with greater probability the external belts.
conditions (direction indicated with another pink arrow) from the external VAB (still in yellow), with the interaction (still in blue) connecting the two phenomena along the minimum path.The geomagnetic lines (in green) of the internal VAB can cross the lithosphere up to 30°-35° in latitude, so being unlikely to be affected by tectonic activity, whereas the radial propagation of LAIC (in blue) should be able to intercept with greater probability the external belts.Given the pivotal role of geomagnetic activity in perturbing the ionosphere, instances where the ap and Dst indices exceeded some predefined thresholds were excluded from the ionospheric data analysis.We also verified that no X-class flares occurred during the detected ionospheric anomalies.

Lithospheric Data Analysis
The acceleration of seismicity prior to the mainshock was analyzed by applying the R-AMR method [48,49] to the INGV Catalog, in order to highlight a diverging power-law function over time for the cumulative value of Benioff strain.The focal search area was centered on the epicenter of the event, whose responsible fault system has a length of approximately 70 km [41].Figure 7 shows the result of the R-AMR analysis, excluding the mainshock and its aftershocks.An evident acceleration is observed, characterized by the C value indicating the onset of "critical" behavior relative to the background, which is 0.598; however, the estimated critical time tf is 200 days after the mainshock.This accelerated behavior was observed in the fault area within a radius ranging from 0 to 50 km from the epicenter.Additionally, the application of this method provides two expected magnitudes, M(A) = 4.8 and M(B) = 4.7, which, although underestimating the real mainshock magnitude, predict an impending EQ that significantly exceeds the background seismicity.The result of this automatic analysis allows us to place the initial progression of seismic acceleration in the lithosphere around September 2019, i.e., 1155 days before the mainshock, resuming in August 2021, i.e., 455 days before the event.These accelerations are visible from the change in slope in the cumulative curve.Around these Given the pivotal role of geomagnetic activity in perturbing the ionosphere, instances where the ap and Dst indices exceeded some predefined thresholds were excluded from the ionospheric data analysis.We also verified that no X-class flares occurred during the detected ionospheric anomalies.

Lithospheric Data Analysis
The acceleration of seismicity prior to the mainshock was analyzed by applying the R-AMR method [48,49] to the INGV Catalog, in order to highlight a diverging power-law function over time for the cumulative value of Benioff strain.The focal search area was centered on the epicenter of the event, whose responsible fault system has a length of approximately 70 km [41].Figure 7 shows the result of the R-AMR analysis, excluding the mainshock and its aftershocks.An evident acceleration is observed, characterized by the C value indicating the onset of "critical" behavior relative to the background, which is 0.598; however, the estimated critical time t f is 200 days after the mainshock.This accelerated behavior was observed in the fault area within a radius ranging from 0 to 50 km from the epicenter.Additionally, the application of this method provides two expected magnitudes, M(A) = 4.8 and M(B) = 4.7, which, although underestimating the real mainshock magnitude, predict an impending EQ that significantly exceeds the background seismicity.The result of this automatic analysis allows us to place the initial progression of seismic acceleration in the lithosphere around September 2019, i.e., 1155 days before the mainshock, resuming in August 2021, i.e., 455 days before the event.These accelerations are visible from the change in slope in the cumulative curve.Around these dates, it is possible to attempt to identify the establishment of different energy transmission channels towards the ionosphere, as hypothesized in LAIC models.
dates, it is possible to attempt to identify the establishment of different energy transmission channels towards the ionosphere, as hypothesized in LAIC models.At the bottom of the main figure, the magnitudes of the involved events are represented: red is used for EQs falling within 37 km from the fault and green those outside that limit.

Atmospheric Data Analysis
SKT and OLR were analyzed to conduct the climatological study.We examined 90 days of ECMWF data preceding the Fano EQ, comparing it with a historical time series spanning the previous 42 years, i.e., from 1980 to 2021.The analysis revealed two highly anomalous days for each atmospheric parameter where the 2022 time series reached the limit of the two standard deviation bands of the historical series.The first anomalous day for SKT occurred on 18 August, i.e., 82 days before the EQ, while the second occurred on 15 September, i.e., 54 days before the EQ (as shown in Figure 8).Furthermore, Figure 9 shows the spatial distribution of these anomalous values, confirming their proximity to the epicenter, especially for the main anomaly.Please note that SKT is defined only on land.
Figure 7. Outcome of the R-AMR algorithm applied to the extracted seismic dataset.The red points represent EQs that are closer to the fault (within 37 km) than those represented by the blue points.At the bottom of the main figure, the magnitudes of the involved events are represented: red is used for EQs falling within 37 km from the fault and green those outside that limit.

Atmospheric Data Analysis
SKT and OLR were analyzed to conduct the climatological study.We examined 90 days of ECMWF data preceding the Fano EQ, comparing it with a historical time series spanning the previous 42 years, i.e., from 1980 to 2021.The analysis revealed two highly anomalous days for each atmospheric parameter where the 2022 time series reached the limit of the two standard deviation bands of the historical series.The first anomalous day for SKT occurred on 18 August, i.e., 82 days before the EQ, while the second occurred on 15 September, i.e., 54 days before the EQ (as shown in Figure 8).Furthermore, Figure 9 shows the spatial distribution of these anomalous values, confirming their proximity to the epicenter, especially for the main anomaly.Please note that SKT is defined only on land.
Geosciences 2024, 14, x FOR PEER REVIEW 11 of 21 dates, it is possible to attempt to identify the establishment of different energy transmission channels towards the ionosphere, as hypothesized in LAIC models.At the bottom of the main figure, the magnitudes of the involved events are represented: red is used for EQs falling within 37 km from the fault and green those outside that limit.

Atmospheric Data Analysis
SKT and OLR were analyzed to conduct the climatological study.We examined 90 days of ECMWF data preceding the Fano EQ, comparing it with a historical time series spanning the previous 42 years, i.e., from 1980 to 2021.The analysis revealed two highly anomalous days for each atmospheric parameter where the 2022 time series reached the limit of the two standard deviation bands of the historical series.The first anomalous day for SKT occurred on 18 August, i.e., 82 days before the EQ, while the second occurred on 15 September, i.e., 54 days before the EQ (as shown in Figure 8).Furthermore, Figure 9 shows the spatial distribution of these anomalous values, confirming their proximity to the epicenter, especially for the main anomaly.Please note that SKT is defined only on land.Regarding the OLR parameter, Figure 10 shows two anomalies exceeding the historical mean by two standard deviations on 5 and 12 September 2022, i.e., respectively, 65 and 58 days before the Fano EQ.Interestingly, both anomalies occurred after the first SKT anomaly but before the second one.Figure 11 depicts the spatial distribution of these OLR anomalies, with the first map (a) revealing an extended structure located to the north of the EQ epicenter.Regarding the OLR parameter, Figure 10 shows two anomalies exceeding the historical mean by two standard deviations on 5 and 12 September 2022, i.e., respectively, 65 and 58 days before the Fano EQ.Interestingly, both anomalies occurred after the first SKT anomaly but before the second one.Figure 11 depicts the spatial distribution of these OLR anomalies, with the first map (a) revealing an extended structure located to the north of the EQ epicenter.Regarding the OLR parameter, Figure 10 shows two anomalies exceeding the historical mean by two standard deviations on 5 and 12 September 2022, i.e., respectively, 65 and 58 days before the Fano EQ.Interestingly, both anomalies occurred after the first SKT anomaly but before the second one.Figure 11 depicts the spatial distribution of these OLR anomalies, with the first map (a) revealing an extended structure located to the north of the EQ epicenter.The application of the ionosonde multiparametric approach [64,65], that takes into account the variations of three ionosonde characteristics, h'Es, fbEs, and foF2, manually  The application of the ionosonde multiparametric approach [64,65], that takes into account the variations of three ionosonde characteristics, h'Es, fbEs, and foF2, manually scaled from hourly ionograms of the AIS-INGV ionosonde of Rome [66], revealed a single anomaly occurring 9 days prior to the Fano EQ (i.e., on 31 October) at 06:00 UT (Table 1; Figure 12).Since the anomaly occurred on a day with the geomagnetic index Ap = 11 nT, the criterion was not strictly satisfied; however, we preferred to take this anomaly into consideration, given the moderate magnitude of the Fano EQ.As shown in Figure 13, it is worth noting that the anomaly is consistent with the relationship between ∆T•R and M previously found by the analysis of the most powerful Central Italian EQs since 1984 [66], with ∆T representing the anticipation time (in days), R the distance (in km) between the epicenter and the ionosonde, and M the EQ magnitude.

Swarm Satellite Data Analysis
The analysis of satellite magnetic field data from the Swarm constellation has reported two electromagnetic anomalies, which could be correlated with the Fano EQ.Specifically, applying the MASS algorithm to the Swarm A satellite, considering the 90 days before the EQ and 10 days after, an anomaly was detected 4 days after the EQ (Figure 14a), potentially associated with the aftershock period.Another anomaly was also detected 75 days before the EQ, very close to the edge of the analyzed track (Figure 14b).The analysis of satellite magnetic field data from the Swarm constellation has reported two electromagnetic anomalies, which could be correlated with the Fano EQ.Specifically, applying the MASS algorithm to the Swarm A satellite, considering the 90 days before the EQ and 10 days after, an anomaly was detected 4 days after the EQ (Figure 14a), potentially associated with the aftershock period.Another anomaly was also detected 75 days before the EQ, very close to the edge of the analyzed track (Figure 14b).Looking at the geomagnetic indices, the geomagnetic activity of these days was quite calm.Then, some candidate EBs were observed over the expected regions.For example, a flux of about 600 electrons cm −2 s −1 str −1 was observed in the NOAA-15 track on November 8 at a latitude of 50-55 • around 20:40 UT, as shown by a red circle in Figure 15.Slight perturbations of around 400 electrons cm −2 s −1 str −1 were also observed at the previous eastward satellite trajectory, all occurring at around 19:00 UT.Thus, these two perturbations anticipated the Central Italy EQ by 9.5 and 11.2 h, respectively.

Discussion
In this paper, precursor anomalies possibly associated with the 2022 M L = 5.7 Fano EQ (Marche, Italy) were studied using a multiparametric and multilayer approach, including seismic, atmospheric, and ionospheric parameters.The purpose of this multidisciplinary approach is to gather various contributions and connect them to identify the best LAIC model.
Although the resolutions of data are different, we are confident that, when we assimilate them at daily intervals (see Table 2), we can reconstruct a reliable cumulative number of anomalies.The tracking of the cumulative number of anomalies in chronological order reveals a distinctive behavior, as illustrated in Figure 16.
Table 2. List of anomalies detected for the case study of the Fano EQ.From right to left, it shows the analyzed parameter, the day when the anomaly was identified, and the days of occurrence relative to the mainshock.All anomalies appear from bottom (lithosphere) to top (atmosphere and ionosphere).There is only one exception (indicated in bold): a satellite anomaly appears among the atmospheric anomalies.multidisciplinary approach is to gather various contributions and connect them to identify the best LAIC model.Although the resolutions of data are different, we are confident that, when we assimilate them at daily intervals (see Table 2), we can reconstruct a reliable cumulative number of anomalies.The tracking of the cumulative number of anomalies in chronological order reveals a distinctive behavior, as illustrated in Figure 16.In a comprehensive approach of the anomalies, the cumulative number of anomalies for Fano EQ is shown here.It is possible to notice that the anomalies appear in time mostly from below (seismic data in the lithosphere) to above (atmosphere and ionosphere).The red curve is an exponential fit of the data.

Figure 16.
In a comprehensive approach of the anomalies, the cumulative number of anomalies for Fano EQ is shown here.It is possible to notice that the anomalies appear in time mostly from below (seismic data in lithosphere) to above (atmosphere and ionosphere).The red curve is an exponential fit of the data.
An exponential fit (indicated by the red curve) accurately represents the overall acceleration of anomalies.This collective progression of anomalies from different geolayers cannot simply be attributed to chance and is probably associated with the preparation phase of the Fano EQ.Furthermore, as highlighted in Table 2, most anomalies appear chronologically from the lithosphere to the atmosphere and ionosphere.This pattern suggests these can be defined as "thermodynamic anomalies", related to a diffusivedelayed coupling model likely driven by thermodynamic processes.Notably, there is an ionospheric satellite anomaly (indicated in bold in Table 2) amidst the atmospheric anomalies, which could be due to direct electromagnetic coupling between the lithosphere and the ionosphere.
Based on the multiparametric and multilayer analysis, a long-term precursor in the lithosphere was identified.For example, the R-AMR value showed its anomalous acceleration starting about three years before the EQ.From the analysis of atmospheric and ionospheric parameters, anomalies were detected starting 82 days before the Fano EQ.Specifically, the energetic particle signal from NOAA showed an anomaly one day before the seismic event.The Swarm satellites detected an anomaly 75 days before the event and another 4 days after, likely associated with aftershocks.An anomaly was recorded in the ionosonde 9 days before the EQ, and atmospheric anomalies were mainly detected at −82 days and −54 days.The number of anomalies in the atmosphere and ionosphere for this EQ is comparable.This similarity suggests a LAIC behavior of the "thermodynamic" or "diffusive-delayed" coupling, with progression from the lithosphere through the atmosphere to the ionosphere.This progression might correspond to the chemical channel or the acoustic gravity channel, as described in the Hayakawa model [1].

Conclusions
In conclusion, this study presents results from the retrospective analysis of the major earthquake that occurred in Italy in November 2022, applying a multiparametric and multilayer approach.This approach involved analyzing data from different geophysical layers (lithosphere, atmosphere, and ionosphere) engaged in the coupling process during the earthquake preparation phase.In the ca.1200 days preceding the Fano EQ, anomalies appeared primarily from the lowest level (seismic data in the lithosphere) to the higher levels (atmosphere and ionosphere), following an overall acceleration pattern.This confirms that the observed anomalies, which originated during the EQ preparation phase and progressed thermodynamically from the lithosphere to the atmosphere and ionosphere, seem to be consistent with the delayed coupling model.However, the presence of a satellite anomaly between other atmospheric anomalies also seems to confirm the possibility of another direct coupling.Therefore, the overall results would confirm a two-way LAIC model.

Geosciences 2024 , 21 Figure 1 .
Figure 1.Summary map of the study conducted: the different layers analyzed are observed, from bottom to top.For each layer, the types of parameters considered are indicated.

Figure 1 .
Figure 1.Summary map of the study conducted: the different layers analyzed are observed, from bottom to top.For each layer, the types of parameters considered are indicated.

Figure 1 .
Figure 1.Summary map of the study conducted: the different layers analyzed are observed, from bottom to top.For each layer, the types of parameters considered are indicated.

Figure 2 .
Figure 2. Simplified structural map of Italy, with the epicenter of the 9 November 2022 Fano EQ indicated by a small yellow star (modified from [40]).

Figure 2 .
Figure 2. Simplified structural map of Italy, with the epicenter of the 9 November 2022 Fano EQ indicated by a small yellow star (modified from [40]).

Figure 3 .
Figure 3. Seismotectonic framework of the coastal area of the Marche region.The light blue squares represent the seismic sequence from 9 November 2022 to 14 February 2023; the first event is marked with a yellow star and the second event is shown with a green star.Historical and instrumental earthquakes from CPTI15 [45] are indicated with colored squares, with earthquakes of Mw ≥ 5.5 highlighted in red.The surface projections of seismogenic zones are depicted with orange ribbons [41].The focal mechanisms of the 9 November 2022 earthquake and the event of 30 October 1930, represented by the grey and white balls, come from TDMT (Time Domain Moment Tensor) and Vannoli et al.[46], respectively (modified from[30]).

Figure 3 .
Figure 3. Seismotectonic framework of the coastal area of the Marche region.The light blue squares represent the seismic sequence from 9 November 2022 to 14 February 2023; the first event is marked with a yellow star and the second event is shown with a green star.Historical and instrumental earthquakes from CPTI15 [45] are indicated with colored squares, with earthquakes of Mw ≥ 5.5 highlighted in red.The surface projections of seismogenic zones are depicted with orange ribbons [41].The focal mechanisms of the 9 November 2022 earthquake and the event of 30 October 1930, represented by the grey and white balls, come from TDMT (Time Domain Moment Tensor) and Vannoli et al.[46], respectively (modified from[30]).

Figure 4 .
Figure 4. Spatial distribution of the 174.723 events extracted from the INGV Catalog during the period 2012-2022 within a circular radius of 150 km from the epicenter of the main EQ, highlighted by the yellow star.The grey and white sphere represents the focal mechanism of the earthquake on 9 November 2022.The chosen radius includes the Central Italy sequence (2016), identifiable by the cluster of events to the south near the edge of the area.

Figure 4 .
Figure 4. Spatial distribution of the 174.723 events extracted from the INGV Catalog during the period 2012-2022 within a circular radius of 150 km from the epicenter of the main EQ, highlighted by the yellow star.The grey and white sphere represents the focal mechanism of the earthquake on 9 November 2022.The chosen radius includes the Central Italy sequence (2016), identifiable by the cluster of events to the south near the edge of the area.

Figure 5 .
Figure 5. Spatial distribution of ten years of seismicity around the mainshock, in a radius of 150 from the epicenter (largest green circle).Blue and red dots (confined within blue and red cir respectively) represent the events contributing to the acceleration found by the R-AMR analy The red dots are the events closer to the seismogenic fault.

Figure 5 .
Figure 5. Spatial distribution of ten years of seismicity around the mainshock, in a radius of 150 km from the epicenter (largest green circle).Blue and red dots (confined within blue and red circles, respectively) represent the events contributing to the acceleration found by the R-AMR analysis.The red dots are the events closer to the seismogenic fault.

Figure 6 .
Figure 6.Scenario of hypothesized pre-EQ coupling processes between the lithosphere of Central Italy and areas where possible EBs could be detected by LEO satellites.

Figure 6 .
Figure 6.Scenario of hypothesized pre-EQ coupling processes between the lithosphere of Central Italy and areas where possible EBs could be detected by LEO satellites.

Figure 7 .
Figure 7. Outcome of the R-AMR algorithm applied to the extracted seismic dataset.The red points represent EQs that are closer to the fault (within 37 km) than those represented by the blue points.At the bottom of the main figure, the magnitudes of the involved events are represented: red is used for EQs falling within 37 km from the fault and green those outside that limit.

Figure 7 .
Figure 7. Outcome of the R-AMR algorithm applied to the extracted seismic dataset.The red points represent EQs that are closer to the fault (within 37 km) than those represented by the blue points.At the bottom of the main figure, the magnitudes of the involved events are represented: red is used for EQs falling within 37 km from the fault and green those outside that limit.

Figure 8 .
Figure 8. Analysis of the SKT parameter for the Fano EQ with comparison between the 2022 time series (dashed red line) and the historical time series (1980-2021, blue line).Evidenced by red circles there are two quite anomalous values near the second standard deviations from the mean: the first one refers to 18 August, and the second one to 15 September.

Figure 8 .
Figure 8. Analysis of the SKT parameter for the Fano EQ with comparison between the 2022 time series (dashed red line) and the historical time series (1980-2021, blue line).Evidenced by red circles there are two quite anomalous values near the second standard deviations from the mean: the first one refers to 18 August, and the second one to 15 September.

Figure 9 .
Figure 9. Maps of the SKT anomalous days in terms of difference with respect to the historical mean: (a) 18 August; (b) 15 September.The EQ epicenter is indicated by the central star.SKT is defined only on land.

Figure 10 .
Figure 10.Analysis of the OLR for the Fano EQ with the identification of two anomalous days that exceed the historical average calculated from 1980 to 2021 by two standard deviations.

Figure 9 .
Figure 9. Maps of the SKT anomalous days in terms of difference with respect to the historical mean: (a) 18 August; (b) 15 September.The EQ epicenter is indicated by the central star.SKT is defined only on land.

Figure 8 .
Figure 8. Analysis of the SKT parameter for the Fano EQ with comparison between the 2022 time series (dashed red line) and the historical time series (1980-2021, blue line).Evidenced by red circles there are two quite anomalous values near the second standard deviations from the mean: the first one refers to 18 August, and the second one to 15 September.

Figure 9 .
Figure 9. Maps of the SKT anomalous days in terms of difference with respect to the historical mean: (a) 18 August; (b) 15 September.The EQ epicenter is indicated by the central star.SKT is defined only on land.

Figure 10 .
Figure 10.Analysis of the OLR for the Fano EQ with the identification of two anomalous days that exceed the historical average calculated from 1980 to 2021 by two standard deviations.

Figure 10 . 21 Figure 11 .
Figure 10.Analysis of the OLR for the Fano EQ with the identification of two anomalous days that exceed the historical average calculated from 1980 to 2021 by two standard deviations.Geosciences 2024, 14, x FOR PEER REVIEW 13 of 21

Figure 11 .
Figure 11.Maps of OLR anomalous days maps in terms of difference with respect to the historical mean: (a) 5 September; (b) 12 September.The epicenter is indicated by the central star.

Figure 12 .
Figure 12.The anomaly observed 9 days before the 9 November 2022 Fano EQ using ∆h'Es, δfbEs, and δfoF2 variations, along with 3 h Kp index values given as a reference of geomagnetic activity.

Figure 13 .
Figure 13.Ionosonde anomaly for the 9 November 2022 M5.7 Fano EQ (red square), compared to the relationship between ∆T•R and M previously found by the analysis of the most powerful Central Italian EQs since 1984 (red line and black squares).

Figure 14 .
Figure 14.Anomalies found 4 days after (a) and 75 days before the Fano EQ (b) by means of an automatic search for magnetic anomalies 90 days before and 10 days after the EQ; MASS algorithm (kt = 2.5) applied to Swarm A satellite.The anomalies are evidenced by coloured rectangles.The vertical red line on the geographical map represents the satellite track.

4. 3 Figure 14 .
Figure 14.Anomalies found 4 days after (a) and 75 days before the Fano EQ (b) by means of an automatic search for magnetic anomalies 90 days before and 10 days after the EQ; MASS algorithm (kt = 2.5) applied to Swarm A satellite.The anomalies are evidenced by coloured rectangles.The vertical red line on the geographical map represents the satellite track.

4. 3 . 3 .
NOAA Satellite Data Analysis: Electron Burst Data Analysis NOAA electron fluxes were analyzed two days before the EQ and on the event day.

Figure 15 .
Figure 15.Three-dimensional representation of the NOAA-15 semi-orbits on 8 November 2022; EB evidenced by a red circle while the star identified the EQ epicenter.

Figure 16 .
Figure16.In a comprehensive approach of the anomalies, the cumulative number of anomalies for Fano EQ is shown here.It is possible to notice that the anomalies appear in time mostly from below (seismic data in the lithosphere) to above (atmosphere and ionosphere).The red curve is an exponential fit of the data.

Table 1 .
Anomaly detected at the ionospheric station of Rome from ionosonde measurements and possibly related to the Fano EQ.