On the Rainfall Triggering of Phlegraean Fields Volcanic Tremors

: We study whether the shallow volcanic seismic tremors related to the bradyseism observed at the Phlegraean Fields (Campi Flegrei, Pozzuoli, and Naples) from 2008 to 2020 by the Osservatorio Vesuviano could be partially triggered by local rainfall events. We use the daily rainfall record measured at the nearby Meteorological Observatory of San Marcellino in Naples and develop two empirical models to simulate the local seismicity starting from the hypothesized rainfall-water effect under different scenarios. We found statistically signiﬁcant correlations between the volcanic tremors at the Phlegraean Fields and our rainfall model during years of low bradyseism. More speciﬁcally, we observe that large amounts and continuous periods of rainfall could trigger, from a few days to 1 or 2 weeks, seismic swarms with magnitudes up to M = 3. The results indicate that, on long timescales, the seismicity at the Phlegraean Fields is very sensitive to the endogenous pressure from the deep magmatic system causing the bradyseism, but meteoric water inﬁltration could play an important triggering effect on short timescales of days or weeks. Rainfall water likely penetrates deeply into the highly fractured and hot shallow-water-saturated subsurface that characterizes the region, reduces the strength and stiffness of the soil and, ﬁnally, boils when it mixes with the hot hydrothermal magmatic ﬂuids migrating upward. The structural collapse of the saturated fractured soil and the mixing of the meteoric ﬂuid with the hot deep ﬂuids triggers the local seismic activity.


Introduction
The active volcanic area including the Phlegraean Fields (Campi Flegrei) and the Somma-Vesuvius complex in the province of Naples (Italy) is exposed to serious seismic hazards due to its high-density population of about 4 million people [1][2][3]. For this reason, the volcanism in Campania has been extensively studied since Roman times. In this work, we evaluate whether a link exists between meteorological rainfall events and the seismic tremors occurring in the Phlegraean Fields. This area includes the municipality of Pozzuoli, Bacoli, Monte di Procida, the Phlegraean islands (Ischia, Procida, Vivara), and Pisciarelli-Agnano ( Figure 1).
Rainfall events and atmospheric pressure variations have been often claimed to influence seismic and volcanic activity [4][5][6][7][8][9][10][11][12][13][14][15]. For example, the 1924 explosions from Halema'uma'u crater located in Hawaii Volcanoes National Park, and many that occurred between about 1500 and 1790, resulted from interactions between magma or lava with ground-or lake-water [16]. Seismicity and ground deformation of Nevado del Ruiz Volcano (Colombia) are also correlated with weekly rainfall [17]. Other studies correlated rainfall and volcanic and seismic phenomena in Montserrat [10][11][12]18]. Similarly, at Guaga Pichinga (Equador) [19] and Mount Etna (Sicily, Italy) [20], deadly volcanic explosions were associated with rainfalls although the volcanologists did not have enough information to fully document it. Extreme rainfall triggered the 2018 rift eruption at Kīlauea Volcano [15]. Similarly, seismic tremors recorded at Mount Vesuvius present a daily cycle  The seismicity of the Phlegraean Fields has been related to its bradyseism since ancient times [2,23]. Bradyseism is the gradual uplift or descent of a terrain area induced by the filling or emptying of an underground magma chamber and/or by hydrothermal activity. This phenomenon happens particularly in some volcanic calderas such as the Phlegraean Fields and Yellowstone in the USA.
In the recent past, this includes the strong uplifting of some regions of the Phlegraean Fields up to 150-170 cm that occurred in 1970-1972 and of 180 cm that occurred in [1982][1983][1984]. Figure A1 in the appendix shows time-series of changes in altitude from 2000 to 2020 of four different stations placed at Pozzuoli-Rione Terra (RITE), Accademia Aeronautica (ACAE), Solfatara (SOLO), and Pozzuoli-Cemetery (STRZ) (adapted from the Instituto Nazionale di Geofisica e Volcanologia-Osservatorio Vesuviano, INGV-OV). At RITE, the soil was lifted by about 650 mm from 2008 to 2020. From 2008 to 2020, more than 3200 seismic tremors occurred. As discussed later, the number of seismic tremors of the region has increased during years of fast soil uplift.
Modeling suggests that the magmatic-hydrothermal system of the area becomes particularly unstable and produces seismic and bradyseism activity when the impermeable crystalline layer surrounding the magmatic system (6 km depth) fractures [24][25][26]. These fractures facilitate the migration of the deep hot fluids of magmatic origin toward the transition zone between the lithostatic and hydrostatic domain.
The deep fluids push from below the impermeable layer located at about 2.5 km of depth and can cause the lifting of the surface and, therefore, a bradyseism crisis. Subsidence would begin when the upper impermeable layer seriously fractures and the deep hot fluids migrate toward the surface and eventually cool down by mixing with the meteoric and marine waters that feed the near-surface hydrothermal system. This is when most of the near-surface tremors would occur, triggered by the underwater boiling.
The first phase refers to the evolution of the status of the deep magma and has a timeperiod of tens to thousands of years. The second transitory phase is made of cyclical events of uplifting and subsidence of the soil and has a timescale from a year to a few decades.
On shorter timescales (e.g., days to weeks), we suggest that the fracturing of the soil of the area facilitates the infiltration of the meteoric and marine waters toward the deeper hot zones. By mixing with deeper hot fluids and rocks, the input of meteoric water could further destabilize the soil and trigger tremors at a depth of about 0-2.5 km [27,28].
The Phlegraean Fields also experience enhanced rainfalls relative to the surrounding areas, which can also trigger flashfloods [29]. Moreover, the region is characterized by several basins (e.g., Astroni crater) that collect large amounts of meteoric water that can easily penetrate the subsurface. Herein, we develop a statistical model based on the local rainfall record to test whether a link between rainfalls and seismic activity could be physical.
Understanding the triggering phenomenology of the local seismic activity is also socially important because, over the past decades, local inhabitants have been alarmed by the erratic and increasing seismic activity of the area and fear the risk that drilling for geothermal purposes, which is rarely allowed by the authorities, could trigger phreatic explosions, gas emissions, and other harmful phenomena. These events could even culminate with the opening of volcanic fumaroles within the urban areas and seismically destabilize an already highly critical region as extensively reported in the latest book on the topic by De Vivo et al. [2].
For example, in June 2020 a new fumarole suddenly opened in the Pisciarelli-Agnano urban area ( Figure 1) because of geothermal drilling [30,31], and the experiment was immediately halted. Thus, it is important to evaluate the seismic sensitivity of the area to water infiltration.

Data
We analyze and compare the following data. (A): The catalog of earthquake events that occurred at the Phlegraean Fields from 2008 to 2019 available from http://www.ov.ingv.it/ov/doc/ct_FLE_2007_12_dic_2019.txt. The record was updated using the data from January to December 2020, available from the main website of the observatory (http://www.ov.ingv.it/ov/it/campi-flegrei.html, accessed on 1 January 2021). The catalog contains several uncertain events whose magnitude was not recorded (indicated as NaN). In the following analysis, these uncertain events were not considered.
(B): The daily catalog of rainfall measured at the first-class Meteorological Observatory of San Marcellino of the University of Naples Federico II, located in the historical center of Naples. This record is representative also of the daily scale of the rainfall measured at Solfatara-Pisciarelli site because the two sites are about 10 km from each other. In Appendix B, we show that similar rain events occur simultaneously in the two locations, although at the Solfatara, on average, it rains 40% more than in Naples because of the specific orography of the place [29]. The data are available from http://www.meteo.unina. it/davis/reports-mensili-e-annuali, accessed on 1 January 2021.
As discussed below in more detail, the yearly number of days with earthquakes greatly increased from 8 in 2008 to 209 in 2020. The total number of seismic events also greatly increased from 53 events in 2008 to 766 events in 2020. The magnitude of the seismic events at the Phlegraean Fields varied from M = −1.6 to M = 3.3. There were only two events with M > 3: one occurred on 6 December 2019 (M = 3.1), and the other occurred on 26 April 2020 (M = 3.3). This area is characterized by a typical Mediterranean climate that has relatively dry summer months and wet autumns, winters, and springs. The daily rainfall varied up to 70 mm in Naples.

Methods
The hypothesis to be tested is whether rainfalls could trigger seismic activity at the Phlegraean Fields with a delay up to a few days or weeks depending on the rainfall amount, as observed in similar cases in other regions [13].
To do this, we develop an empirical model inspired by Darcy's law [32,33], which describes the volumetric water flow dependency on the pressure difference between the two sides of a soil sample. The proposed model has to simulate a plausible time evolution of the wet-state of the soil as a function of the timing and amount of the rainfall events. We need a time-delayed function that varies with the rainfall amount that can be directly compared against the seismic record. In fact, rainfall water needs some time to percolate inside the soil and reach critical saturated regions where a seismic tremor could occur and, obviously, if the rainfall amount is larger, the effect could last longer. A function with such properties can be obtained in the following way.
If t i (i = 1, 2, 3, . . . N) is the timing of consecutive rainy days, and RA(t i ) is the rainfall amount occurring in the day t i , we propose a time-delayed model, F(t), defined by the following equations: where N is the number of rainfall days, and F i (t, t i ) gives the contribution of the single rain event that occurred at time t i . The parameter τ is a characteristic timescale describing, for example, the water percolation in the soil. The effect is supposed to weaken exponentially following Darcy's law [32,33], which states that the volumetric flow depends on the pressure difference between the two sides of a sample. The variant proposed in Equation (2) describes how water flows vertically through saturated soil by gravity. In our case, the pressure difference between two levels is altered by the addition of meteoric water on the top side, that is, at the surface. If h(t) is the amount of the meteoric water at the surface at a given time t, it induces an additional pressure equal to Stevino's law p(t) = ρgh(t) above the stationary condition, where ρ is the water density and g is the constant of gravity. As this water penetrates the ground, the pressure at the surface diminishes at a rate that is proportional to the top pressure itself, that is, as h(t) = −τ dh(t)/dt, where τ is the characteristic time of the process.
If, at time t 0 , a given amount of meteoric water h(t 0 ) is added to the surface, the solution of the above differential equation is h(t) = h(t 0 ) exp(−(t − t 0 )/τ), which is like Equation (2). For large t, h(t) tends towards 0, which indicates that all meteoric water added to the system at time t 0 has been absorbed by the soil. Then, the system returns to its stationary state unless it is again perturbed by additional meteoric water. This cumulative process is related to all rainfall events and is indicated by Equation (3), which is herein defined in the units of the rainfall amount (mm) in Naples.
In Darcy's law, the characteristic time τ is given by L/K, where L is the specimen length, which is the depth of the soil, and K is the hydraulic conductivity of the soil, which is the speed of the water in the ground. We also have K = krg/m, where m is the water viscosity and k is the permeability of the soil.
A threshold R can be chosen such that, only when F(t) > R, the meteoric perturbation could affect the seismic activity. Thus, we define wet days as those for which F(t) > R and dry days as all the others. With this threshold, Equation (2) means that the larger the rainfall amount RA(t i ) on a given day t i , the longer its ability to trigger tremors lasts.
In the following, we evaluate the statistics of having seismic tremors on wet days or on dry days to determine whether the wet days are more likely to experience seismic activity than the dry days.
The values assigned to R and τ are heuristically chosen because our purpose is to determine whether there exists a rainfall model of the type given by Equations (1)-(3) that correlates rainfall events with the seismic events with a delay from a few days to a few weeks.
By assuming a specimen length of L = 1000 m (which is about the mean depth where the seismic tremors occur at the Phlegraean Fields) and τ = 1 day, the hydraulic conductivity of the soil is K = L/τ = 1.2 cm/s. Such a value of K corresponds to soils made of karst limestones, highly fractured rocks, sand, and/or gravel [34]. Indeed, for the first 2.7 km depth, the Phlegraean Fields are characterized by such a kind of soil [35].
In general, meteoric water can percolate down quite fast through fractured soil [36] and can rapidly affect an already water-saturated and critically unstable underground. With a value of τ equal to about 1 or 2 days, the required time-delayed effect from a few days to 1 or 2 weeks can be modeled with R values around 0.1 or 0.2 mm, which usually correspond to the pluviometry sensitivity. Note that 0.1 mm corresponds to 0.1 kg of rainfall water on a square meter. Thus, these thresholds are heuristically determined.
We run several simulations using different values of R and τ near the values estimated above. In the following, we show our results for four scenarios: (1) R = 0.2 mm and τ = 1 day; (2) R = 0.1 mm and τ = 1 day; (3) R = 0.2 mm and τ = 2 day; (4) R = 0.1 mm and τ = 2 day.
From scenarios #1 to #4, the model allows the rainfall water to have more and more time (from a few days to 1-2 weeks) to affect the seismic activity in the area. For each scenario and each year from 2008 to 2020, we evaluate the following observables: Finally, for each year and scenario, we evaluate the relative statistics made of the following measures: • the probability of having random earthquakes with no relation with wet days (Rand = WD divided by DY), which is our null hypothesis to be tested; • the measured probability of having a seismic day during a wet day (WEQD = WDE divided by DE); • the measured probability of having an earthquake during a wet day (WEQ = EWD divided by E).
The last step of our statistical analysis evaluates whether the measured probabilities differ significantly from the random hypothesis. To do this, we obtained the following exact probability equations derived from an extension of the binomial distribution and considering all possible combinations: where α = 0.5 for i = n 1 , and α = 1 for i > n 1 . The meaning of the two equations is the following. Equation (4) calculates the probability P 1 of obtaining a number of wet days with earthquakes larger than measured n 1 , assuming a random distribution of days with seismic activity. Here, n is the number of days with earthquakes, n 1 is the number of wet days with earthquakes, d is the total number of days (e.g., 365 or 366 for one year), and d 1 is the number of wet days in the same period (e.g., 1 year). The equation assumes that if (d 1 − i) < 0, then 1/(d 1 − i)! = 0 because x! = Γ(x + 1) and, for negative integers, the gamma function diverges to infinity. In particular, in the extreme cases in which, for example, our proposed model (Equation (3)) would cover the entire year with wet days or, alternatively, earthquakes occur all days of the year, the probabilities P 1 and P 2 give values equal to 0.5, which means that no relation could be established between rainfalls and seismic activity.
Equation (5) calculates the probability P 2 of obtaining a number of earthquakes in wet days larger than measured (n 1 ) assuming a random distribution of seismic events. This equation assumes that more than one seismic event could occur on the same day. Here n is the number of earthquakes, n 1 is the number of earthquakes that occur on wet days, d is the total number of days, and d 1 is the number of wet days.
We compare these two statistical tests because the seismic data tend to form swarms and clusters where single events may not be physically independent of each other. The comparison between the two tests is necessary to establish whether and how meteoric water affects the seismic properties of the soil on these short timescales. The statistical test is repeated for each scenario and for each year from 2008 to 2020 to get an ensemble of independent results that are then statistically evaluated.
The above procedure is needed for correctly interpreting the test based on P 2 (Equation (5)). In fact, even one large seismic swarm, which can occur on a wet day or a dry day, could bias the statistics in one way or the other. Thus, finding an extreme value regarding P 2 (e.g., P 2 < 0.05 or P 2 > 0.95) by using only a single record means that in such a period, an excess of events occurred during wet or dry days, respectively. However, because of the clustering nature of the phenomenon, the finding could have been occasional. Thus, running a single test would be inconclusive.
However, by using an ensemble of independent tests using one-year sub-periods from 2008 to 2020, we obtain 13 full and independent 1-year segments that can be statistically tested. For example, if the same result (for example, P 2 < 0.05, which implies that meteoric water could enhance the seismic activity of the area) is found for a significantly large majority of the analyzed years, the working hypothesis cannot be rejected.

Results
The statistics referring to scenarios #1, #2, #3, and #4 are reported in Tables 1-4, which summarize the graphs depicted in Figures 2-5 for each year from 2008 to 2020. These figures show the timing of the earthquake events with their relative magnitude against the four rainfall models given by Equation (3) according to the four scenarios. The final results regarding the distributions of the probabilities P 1 and P 2 reported in the four tables are finally collected in Figure 6.
Tables 1-4 report the statistics associated with each year for the four scenarios listed above. We divided the probability range from 0 to 1 into 5 sectors: p ≤ 0.05, 0.05 < p ≤ 0.30, 0.30 < p < 0.70, 0.70 ≤ p < 0.95, 0.95 ≤ p < 1.00. A value of p ≤ 0.05 means that the null hypothesis has a probability to be true that is less than 5%, which also means that the probability that earthquakes are triggered by rainfall events is higher than 95%. A value of p ≥ 0.95 means the opposite. The alternative values of p mean a moderate bias in one direction or the other if 0.05 < p ≤ 0.30 or 0.70 ≤ p < 0.95, respectively. Values 0.30 < p < 0.70 indicate a neutral condition, which means that no relation in one way or the other can be established between rainfalls and seismic activity. The Tables 1-4 highlight in bold the cases p ≤ 0.05 and 0.05 < p ≤ 0.30, which support a relation between rainfalls and seismic activity.
From scenarios #1 to #4, the model gives the rainfall more and more time (from a few days to 1-2 weeks) to affect the seismic activity in the area. The analysis suggests that the seismic days (that is, days when at least one seismic event has been recorded) are not particularly linked with wet days. In fact, for the analyzed scenarios and for the 13 years from 2008 to 2020, we almost always found 0.05 < P 1 < 0.95. One year (scenarios #2 and #4) and 2 years (scenario #1) had P 1 ≥ 0.95, and one year had P 1 ≤ 0.05 in all scenarios. Therefore, by considering the 13 analyzed years, no bias in one direction or the other can be established according to the first test.
On the contrary, the values P 2 reported in Tables 1-4 suggest that the number of seismic-events could be linked with wet days for most years. That is, seismic swarms are more likely to occur in these conditions. In fact, for each scenario, we found the following statistics: Thus, scenario #4 shows the strongest evidence for demonstrating a time-delayed correlation between the rainfalls and seismic events in the Phlegraean Fields. In fact, the values P 2 , reported in Table 4, suggest that the number of seismic events could be linked with wet days at least for 9 of the 13 analyzed years, that is, in 2008,2009,2010,2011,2012,2013,2014,2017, and 2019.
The statistical results referring to the other three scenarios assume a more rapid seismic response to water infiltration: in these cases, we found that the years 2015 and 2018 are characterized by P 2 ≤ 0.05. Figures 2-5 also show that earthquake events tend to cluster. More specifically, seismic swarms are observed and last up to about one hour. For example, on 23 January 2009, there were 143 tremors (without counting those whose magnitude was unknown) in 36 min from 5:15 a.m. to 5:51 a.m. We could not find tremor swarms extending from a day to the following one. Thus, events occurring on different days belong to different clusters. Several seismic events are also isolated; that is, they do not generate any seismic swarm even if their magnitude is relatively high. This result suggests that the seismic properties of the area are not stationary in time. In fact, the physical properties of the soil are expected to be modified by the addition of meteoric water.
Figures 2-5 also suggest that many seismic events occurred during wet days, that is, up to a few days or 1-2 weeks after large rainfall events. This behavior is more evident for scenario #4 and from 2008 to 2014 than from 2015 to 2000, when the number of seismic events increased significantly.
It is worth noting that, if the seismic events significantly increase, its correlation with the rainfall record would also increase. This makes the association between rainfalls and seismic activity statistically and visually more challenging. However, Equations (4) and (5) already take this problem into account by evaluating the probability of obtaining a number of wet days with earthquakes larger than that measured n 1 , assuming a random distribution of days with seismic activity.  For the four analyzed scenarios #1, #2, #3, and #4, Figure 6 (top panel) shows the number of years when the probabilities P 1 fall within each of the five probability ranges for the years 2008-2020. These histograms do not show a significant bias toward wet or dry days. At most, it is possible to observe a slight skewness toward dry days. Figure 6 (bottom panel) depicts the number of years when the probabilities P 2 fall within each of the five probability ranges for the years 2008-2020. From scenario #1 to #4, the number of years when P 2 ≤ 0.05 increases from 6 to 9; the number of years when P 2 ≥ 0.95 decreases from 3 to 1. Thus, scenario #1 would already support our hypothesis. However, the statistical significance of a positive correlation significantly improves for scenario #4 when the probability that rainfalls could affect seismic events in each given year (P 2 ≤ 0.05) exceeds the alternative possibility by 9 to 1. Two other years (2018 and 2020) had 0.05 < P 2 ≤ 0.30; only one year (2015) had 0.70 ≤ P 2 < 0.95, and one year (2016) had P 2 ≥ 0.95.       Table 3.    Tables 1-4 from 2008-2020. The upper panel does not show any prevailing tendency, which suggests that the seismic activity exists because of the bradyseism. However, the bottom panel shows a strong prevailing tendency toward P2 ≤ 0.05 (blue) that increases from scenario #1 to scenario #4, which suggests that the seismic swarms are mostly linked to rainfalls. The results indicate that wet days enhance the seismic activity in the Phlegraean Fields.
For the four analyzed scenarios #1, #2, #3, and #4, Figure 6 (top panel) shows the number of years when the probabilities P1 fall within each of the five probability ranges for the years 2008-2020. These histograms do not show a significant bias toward wet or dry days. At most, it is possible to observe a slight skewness toward dry days. Figure 6 (bottom panel) depicts the number of years when the probabilities P2 fall within each of the five probability ranges for the years 2008-2020. From scenario #1 to #4, the number of years when P2 ≤ 0.05 increases from 6 to 9; the number of years when P2 ≥ 0.95 decreases from 3 to 1. Thus, scenario #1 would already support our hypothesis. However, the statistical significance of a positive correlation the significantly improves for scenario #4 when the probability that rainfalls could affect seismic events in each given year (P2 ≤ 0.05) exceeds the alternative possibility by 9 to 1. Two other years (2018 and 2020) had 0.05 < P2 ≤ 0.30; only one year (2015) had 0.70 ≤ P2 < 0.95, and one year (2016) had P2 ≥ 0.95.
Thus, for scenario #4, by applying the binomial distribution, we obtain less than 1% probability for obtaining more than 9 "yes" (P2 ≤ 0.05) versus 1 "no" (P2 ≥ 0.95) with 10 independent trials (if the other three years are ignored) or for obtaining more than 9 + 2 = 11 "yes" (P2 ≤ 0.30 or P2 < 0.50) versus 1 + 1 = 2 "no" (P2 ≥ 0.70 or P2 > 0.50) with 13 independent trials (if the leaning years are respectively associated with their correspondent closest hypothesis). We would still get just a probability of 9% for obtaining more than 9 "yes" (P2 ≤ 0.05) versus 4 "no" (P2 > 0.05) with 13 independent trials (if the other four cases are considered contrary to our hypothesis).  Tables 1-4 from 2008-2020. The upper panel does not show any prevailing tendency, which suggests that the seismic activity exists because of the bradyseism. However, the bottom panel shows a strong prevailing tendency toward P 2 ≤ 0.05 (blue) that increases from scenario #1 to scenario #4, which suggests that the seismic swarms are mostly linked to rainfalls. The results indicate that wet days enhance the seismic activity in the Phlegraean Fields.
Thus, for scenario #4, by applying the binomial distribution, we obtain less than 1% probability for obtaining more than 9 "yes" (P 2 ≤ 0.05) versus 1 "no" (P 2 ≥ 0.95) with 10 independent trials (if the other three years are ignored) or for obtaining more than 9 + 2 = 11 "yes" (P 2 ≤ 0.30 or P 2 < 0.50) versus 1 + 1 = 2 "no" (P 2 ≥ 0.70 or P 2 > 0.50) with 13 independent trials (if the leaning years are respectively associated with their correspondent closest hypothesis). We would still get just a probability of 9% for obtaining more than 9 "yes" (P 2 ≤ 0.05) versus 4 "no" (P 2 > 0.05) with 13 independent trials (if the other four cases are considered contrary to our hypothesis).
Scenario #4 also suggests that a link between earthquakes and rainfall events could be stronger when fewer earthquakes occur, which usually happens when the bradyseism of the area is weak. In fact, Figure A1 shows that the years from 2002 to 2014 had a relatively modest vertical land movement, which was momentarily interrupted by the 2012 uplift [37]. On the contrary, from 2015 to 2020 a continuous and strong uplift was observed with just a short pause around 2017. According to Table 1, P 2 ≤ 0.05 from 2008 to 2014 and in 2017 and 2019. Thus, we always found (=100% of the trials) P 2 ≤ 0.05 during the 7 years when the bradyseism was relatively modest (2008)(2009)(2010)(2011)2013, 2014 and 2017).

Discussion and Conclusions
Our analysis suggests that rainfall events in the Phlegraean Fields can trigger, with a delay up to some days or weeks (which depends on the rainfall amount), a majority of the seismic events locally recorded. There should also be an additional input from the seawater infiltration, whose level increases with the reduction in atmospheric pressure occurring when it rains [4,28,38].
The novelty of our result relative to previous studies relating atmospheric events and seismic or volcanic activity [4][5][6][7][8][9][10][11][12][13][14][15] resides in the fact that this area is in a seismic critical state because it is subjected to positive bradyseism, which is the first reason why the region is seismically active. Thus, we found that rainwater infiltration triggers local tremors on short time scales by interacting with the endogenous conditions of the underground driven by the pressure of the active magmatic chamber located as near as 6 km below the surface while the impermeable layer is located at about 2.5 km below the surface.
Our statistical analysis highlighted that the observed seismic events usually occur in large clusters, particularly after large rainfalls lasting for several days or weeks. The link between rainfall and seismic activity appears more evident when the annual number of seismic events was small (up to about 200), as it happened from 2008 to 2014. When the number of seismic events greatly increases, more aftershock tremors could be triggered by the previous tremors, and the process would last longer. In this situation, it is more difficult to track the seismic sequence and to identify their cause. This could have happened in 2016 (288 events).
However, using statistical inference, we could detect an association between rainfall events and seismic activity P 2 ≤ 0.05 also in 2018, which experienced 375 seismic events, and in 2019, with 592 recorded seismic events. Even in 2020, when 766 seismic events were recorded, our statistical approach could still detect a moderate correlation 0.05 < P 2 ≤ 0.30 between rainfall and seismic events. In 2020 it is possible to visually notice a very large increase in seismic activity from September to December, which was also the most rainy period of the year. In fact, from January to August the rainfall amount was 205.7 mm and there were 339 seismic tremors; from September to December the rainfall amount was 571.4 mm and there were 437 tremors.
In conclusion, large seismic swarms often occurred in correspondence, or closely following significant rainy periods. Indeed, as rainfalls reduce the soil suction pressure, a wide area around Naples and Vesuvius is affected by frequent landslide phenomena due to the instability of the pyroclastic cover influenced by wetting-drying cycles resulting in shallow structural collapses, fast landslides, flow liquefaction, and more [39]. An additional effect could derive from seawater variation due to tides and atmospheric pressure variation [28,38].
Several studies have evaluated the hydrogeological cross-sections at the Phlegraean Fields [36,[40][41][42][43]. They found that the terrain below the area centered in the Solfatara crater (see Figure 1) is characterized by a magmatic chamber below about 6 km from the sea level. The magmatic chamber is surmounted by a mixing zone at a temperature above 350 • C until about 2.5 km. This region is followed by a large hydrothermal zone saturated with liquids of meteoric origin and gas up to the sea level, with temperatures ranging from about 100 • C to 220 • C [41,43]. During the last years (2018-2020) the number of seismic events increased significantly, despite not having a significant increase in rainfall as compared to earlier years. In fact, the area is seismically active mostly because of the increased endogenous pressure from the deep magmatic chamber, which accelerated the lifting of the ground.
Thus, rainfall water likely penetrates deeply in the highly fractured and hot shallowwater-saturated subsurface from 0 to 2.5 km depth that characterizes the region, reduces the strength and stiffness of the soil, and, finally, boils when it mixes with the hot hydrothermal magmatic fluids migrating upward. The structural collapse of the saturated soil and the mixing of the meteoric fluid with the hot deep fluids should be what triggers the local seismic activity.
However, we also found a number of days experiencing some seismic events in the so-called dry days, which could have happened independently of rainfall events. Thus, rainfall events do not always trigger earthquakes in the Phlegraean Fields, nor do all seismic events appear to be triggered by water infiltration. In fact, the area is seismically active not because it rains but because of the soil instability induced by the local bradyseism [2]. The rainfall water infiltration simply perturbs an already highly critical system. The overall result suggests that meteoric water infiltrating the soil increases the seismic sensitivity and activity of the local faults by favoring large swarms once a seismic tremor occurs, while the local bradyseism keeps the area in a seismic critical state where even small perturbations can have observable triggering effects.
We speculate that the area is particularly sensitive to rainfalls also because of its numerous and large basins of volcanic origin. For example, large amounts of rainfall water are collected by the Astroni crater (the basin is 2,500,000 m 2 ), where a small lake of about 50,000 m 2 also exists; by the basin of Agnano (about 3,000,000 m 2 ), where a large volcanic lake of about 1,000,000 m 2 formed in the 11th century and existed until 1870 when it was artificially drained; and by the Solfatara crater itself. The first two large basins and the sea surround the Solfatara crater, which is the area most affected by seismic activity (Figure 1). Other craters and basins are on the Western sides of the Phlegraean area, where Avernus lake is located. Moreover, the underground hot aquifers of the area are very shallow, as demonstrated by the numerous thermal springs of the area and by the boiling mud at Solfatara (Figure 1). Thus, the results of this work envisage the possibility of developing specific rainfall-draining systems or other hydrological strategies to possibly mitigate the seismic activity at the Phlegraean Fields.
In conclusion, when the shallow crust is in a critical seismic state, as happens at the Phlegraean Fields because of its bradyseism, even tiny pressure and permeability changes induced in the soil by meteoric water infiltration could trigger earthquakes swarm from near the surface to up to 2.5 km depth [13,44].
It is also interesting to observe that our result confirms the intuitions of Neapolitan scientists of the 18th and 19th centuries. In fact, the phenomenology of lava deposits and fumaroles at the Vesuvius and Solfatara convinced some (e.g., Gaetano De Bottis [45] and Giovanni Maria Della Torre [46]) that meteorological events were correlated with volcanic activity. In addition, the Vesuvian Observatory (Osservatorio Vesuviano, the oldest volcanology observatory in the world: http://www.ov.ingv.it/ov/en.html), was founded in 1841 to study the influence of meteorology on volcanological phenomena. In this regard, we note the work of Luigi Palmieri, who was appointed director of the Osservatorio Vesuviano from 1854 to his death in 1896, and who, simultaneously, was also appointed in 1860 director of the Meteorological Observatory of the University of Naples Federico II. Thus, there is an ancient and consolidated scientific tradition in Naples linking rainfall events to the local volcanic and seismic activity.
Author Contributions: N.S. contributed to the project, the data collection and analysis, and to the writing of the paper; A.M. contributed to the project, the data collection and the writing of the paper. Both authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement:
It is no needed.

Informed Consent Statement:
It is no needed. Data Availability Statement: All data are freely available online at the websites listed in Section 2.

Conflicts of Interest:
The authors declare no conflict of interest.
The figure and its insert demonstrate that similar rain events occurred simultaneously in the two locations, although the amount of rainfall is usually not the same. In 2019, the meteorological station at the Solfatara registered 1385 mm of rain versus 1014 mm measured at San Marcellino, which suggests that on average the former is about 40% larger than the latter, at least in 2019. This result is expected because the rainfall amount at Pozzuoli is enhanced by the local orography and the presence of the nearby island of Ischia, which also induce frequent flash-floods. For additional details on the statistics of the rainfall record in Naples see Ref. [47].