Tracking the Endogenous Dynamics of the Solfatara Volcano (Campi Flegrei, Italy) through the Analysis of Ground Thermal Image Temperatures

Abstract

In the last decades, thermal infrared ground-based cameras have become effective tools to detect significant spatio-temporal anomalies in the hydrothermal/volcanic environment, possibly linked to impending eruptions. In this paper, we analyzed the temperature time-series recorded by the ground-based Thermal Infrared Radiometer permanent network of INGV-OV, installed inside the Solfatara-Pisciarelli area, the most active fluid emission zones of the Campi Flegrei caldera (Italy). We investigated the temperatures’ behavior in the interval 25 June 2016–29 May 2020, with the aim of tracking possible endogenous hydrothermal/volcanic sources. We performed the Independent Component Analysis, the time evolution estimation of the spectral power, the cross-correlation and the Changing Points’ detection. We compared the obtained patterns with the behavior of atmospheric temperature and pressure, of the time-series recorded by the thermal camera of Mt. Vesuvius, of the local seismicity moment rate and of the CO₂ emission flux. We found an overall influence of exogenous, large scale atmospheric effect, which dominated in 2016–2017. Starting from 2018, a clear endogenous forcing overcame the atmospheric factor, and dominated strongly soil temperature variations until the end of the observations. This paper highlights the importance of monitoring and investigating the soil temperature in volcanic environments, as well as the atmospheric parameters.

1. Introduction

In the last decades, thermal remote sensing has been largely used in surveys and in the analysis of thermal behavior of active volcanic areas. Different specific devices, such as portable or ground/airborne sensors and satellites [1,2,3,4], can reveal thermal anomalies linked to volcanic activity, as well as mapping volcanic deposits [5]. Thermal remote explorations are largely used to monitor the state of volcanoes by detecting thermal surface variations likely linked to endogenous dynamics and possible eruption precursors. The thermal infrared ground camera is an effective tool to reveal such variations. Indeed, thermal ground monitoring allowed detection of the intensification of the thermal activity prior to the 18 May 1980 eruption at Mount St. Helen’s [6]. Some days before the Santa Ana volcano lake eruption, Hernández et al. [7] measured a significant increase in the surface lake temperature and an intensification of the fumarolic emissions, by using thermal infrared images acquired with a handheld camera. A handheld thermal camera was also used by Calvari et al. [8] to monitor the Mt. Etna and Stromboli volcanoes. When both the volcanoes erupted in 2002–2003, the acquired thermal images showed a failure on the volcanoes’ flanks before the eruptive fractures’ opening. Hilman et al. [9] used the thermal infrared radiometers (TIR) on board of an Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) to detect thermal surface anomalies related to the gas emissions around Papandayan and Domas craters (Indonesia). Their results suggested that SO₂ emission measurements and land surface temperature corrections are highly correlated. Afterward, those authors [10] applied the satellite imaging to increase the effectiveness of upflow zone detection in those areas.

The surface thermal behavior of Campi Flegrei caldera (CFc; Southern Italy) has been monitored using different thermal infrared (TIR) devices: satellite thermal sensors [11], ground-based TIR permanent network (TIRnet) [12,13,14,15] and handheld TIR cameras [16]. The relative studies are mainly devoted to analyzing the thermal observables at a local scale, inside the caldera. Indeed, the interaction between some atmospheric parameters and TIRs, such as the rainfall episodes and variations in ground temperature in Solfatara, was already explored [17]. At even smaller scales, the heating due to anthropogenic noise is also known [18,19], as well as the effects due to the topography which are expected to cause local wind perturbation that increases the gas flux.

The present work mainly aims at the identification of possible sources at larger spatial scales (tens of kilometers) that might affect the CFc thermal behavior, in order to better characterize the local hydrothermal/volcanic dynamics. Besides the spatial investigation, the analysis was also focused on the variations in TIRs over large time scales (from weeks to years), as at small time scales (of the order of a few days).

Here, the behavior of a Solfatara–Pisciarelli shallow hydrothermal system is investigated in the time span of 25 June 2016–29 May 2020, by processing data acquired by TIR ground cameras belonging to the surveillance TIRnet of INGV in correspondence to selected active fumarole emission areas of CFc (Figure 1). A first comparison was carried out among all the TIRs of CFc, the atmospheric parameters (T_atm and P_atm) and the TIR of Mt. Vesuvius (VES), to identify/discriminate possible exogenous sources at a large spatial scale. Then, the evidenced CFc ground thermal patterns were compared to the spatio-temporal behavior of the seismicity rate [20,21] and CO₂ emission flux [22] in order to better characterize possible endogenous sources acting in the first 3 km of depth [23]. Such endogenous sources are linked to the hydrothermal/volcanic dynamics, such as magma movements, pressured gas injection at depth, regional or local stress field variations and so on.

This work evidenced how TIRs are able to track endogenous hydrothermal/volcanic effects, as well as to identify an overall general atmospheric trend.

2. Materials and Methods

2.1. The Study Area

The Campi Flegrei volcanic area is a resurgent nested caldera about 65 km² wide, located in Southern Italy (Figure 1) and generated by two high magnitude volcanic eruptions: the Campanian Ignimbrite (~40 ka) and the Neapolitan Yellow Tuff eruption (~15 ka) [25,26,27,28]. After the last eruption (Monte Nuovo, 1538 AD [29], CFc was affected by a diffuse subsidence until the second half of the 20th century. After this date, ground uplift started with alternating unrest episodes, characterized by an increase in uplift velocity and seismic crises, termed ‘bradyseism’ [30,31]. The most recent bradyseismic episodes occurred in 1968–1972 and 1982–1985.

In 2004, CFc started a new inflation phase, characterized by an increasing ground uplift [32,33,34,35,36] consisting in a total maximum vertical displacement of ~0.80 m [22] by frequent shallow seismicity [20,37,38,39,40] and by increasing degassing activity [23,41].

According to Chiodini et al. [23], the bradyseismic phenomena are due to the dynamics of the magmatic–hydrothermal system. The authors depict a deep gas accumulation zone located at about 4 km of depth, possibly related to small magma/fluid injection, and a hydrothermal reservoir at about 2 km of depth, where upwelling magmatic fluids mix with meteoric water.

Hydrothermal/volcanic fluid emission is mainly centered in the Solfatara–Pisciarelli area (Figure 1) characterized by the Solfatara ~3.8 ka old tuff-cone. The CO₂ anomalies’ pattern follows the structural alignments, suggesting that the faults and the fractures mainly control the diffuse gas emission [42]. The main fault system is along the NW–SE direction, while the minor one follows the NE–SW trend (Figure 1b). Bocca Grande and Bocca Nuova are the main fumarole vents inside the Solfatara crater and are located at the intersection between the two main fault systems (Figure 1b). The Pisciarelli fumarole field, elongated toward the NW–SE faults system, has shown a significant growth in hydrothermal activity during the ongoing unrest phase, with a recent increase in shallow seismicity [22].

2.2. Thermal Infrared Surface Data and Atmospheric Temperature and Pressure Data

The Osservatorio Vesuviano (OV) is a Section of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) mainly designated to the surveillance of the Neapolitan volcanic area by managing a large number of monitoring networks. Since 2006, the Thermal InfRared Permanent Network (TIRnet) has been monitoring the ground thermal anomalies due to diffuse degassing activity at CFc and the Vesuvio crater [12,13,43,44]. The network is overall composed of five stations at CFc and one at the Mt. Vesuvius crater (Figure 1). In the CFc area, three stations are installed inside the Solfatara cone (SF1, SF2 and SOB), one station in the Pisciarelli area (PIS), one station at the southern slope of Mt. Olibano (OBN). The SF1 station acquires frames of the SE inner slope of the Solfatara crater, including Bocca Grande and Bocca Nuova fumaroles. The SOB station acquires frames of a restricted diffuse degassing area on the top of the SE rim of the Solfatara edifice. The SF2 station targets the cryptodome beside the inner NW slope of Solfatara. The PIS station monitors the area close to the base of the NE Solfatara external slope (Figure 1b). The TIR camera installed at Mt. Vesuvius (VES) is located about 24 km away from the Solfatara–Pisciarelli area and it targets the fumaroles located inside the crater in the western inner slope. The distances between the sensors and the target areas range from 80 m to about 340 m. The TIR stations use FLIR infrared cameras A645SC/A655SC holding a Focal Plane Array (FPA) uncooled Microbolometer detector. The sensor thermal sensitivity is <0.03 °C ± 30 °C, the accuracy ± 2 °C, the spectral range 7.5–14 µm and the resolution 640 × 480 pixels. Every TIR station captures three frames of the target areas, at nighttime (00:00, 02:00, 04:00 a.m.) to avoid the effects of sun heat radiation accumulated by the soil.

UMTS and WiFi transmission systems are used to telemeter the acquired scenes to the INGV-OV monitoring center in order to be processed daily by the Matlab© application ASIRA (Automated System of InfraRed Analysis, Naples, Italy) [44]. Raw TIR images are processed to generate thermal infrared surface temperature time-series with daily maximum TIR temperature values (°C; Figure 2).

In this work, ground temperatures recorded in the time span 25 June 2016–29 May 2020 (Figure 2) were used.

Air temperature (T_atm, in °C) and atmospheric pressure (P_atm, in mbar) time-series are recorded by a VAISALA WXT520 meteorological station, installed inside Capo Miseno lighthouse (Bacoli, Italy). The station is composed of an OTT data logger LogoSens2 which acquires a sample every ten minutes and then transmits the data to INGV-OV center by HIPERLAN (HIgh PErformance Radio LAN).

T_atm shows drastic minima and maxima that interrupt the overall trends. These anomalies correspond to the time range in which the air temperature is unusually cold or hot. For example, in January 2017, southern Italy was affected by an intense cold, with anomalous low altitude snowfalls, while January 2018 was unusually hot.

All the time-series are pre-processed, i.e., de-seasonalized (correction for annual periodicity, i.e., 365 d) and corrected for tidal constituents (Mf lunisolar fortnightly, 1366 days; Msf lunisolar synodic fortnightly, 1476 days; Mm lunar monthly, 2755 days; Msm solar monthly, 3181 days; Ssa solar semiannual, 18,263 days; Sa solar annual, 36,526 days) according to the procedure described in [14,15,44].

The temporal trends of the time-series are best highlighted by a moving average over 31-day-long time windows (about 1 month) shifted by 1 day (red line in the Figure 2).

2.3. ICA Method

In order to reveal which among the waveforms related to the frequencies involved in the TIR time-series can be considered independent, we adopted the Independent Component Analysis (ICA), that is, a decomposition method in the time domain [45].

This technique exploits the statistical independence of the sources, whose linear combination is supposed to form each experimental series (each TIR). At a first approximation, the kurtosis can be used for estimating the independence; a more general approach requires that the negentropy and/or the mutual information be evaluated [45].

The mathematical model can be formulated considering that each TIR is a mixture of some unknown sources, x = As, where x represents the five TIR series, the matrix A provides the coefficients of the linear superposition and basically takes into account the medium effects, and s is a vector of at most five unknown independent sources. The model requires that the number of the recordings, x, must be at least equal to the number of the unknown sources, s. Under these hypotheses, the ICA goal is to obtain a matrix of separation W ≈ A⁻¹, such as the vector y = Wx being an estimate y∼s of the original independent signals. Several approaches are proposed in the literature, such as the Mutual Information, the maximum likelihood, etc., to retrieve the sources and W. To estimate them, we adopted the non-Gaussianity method, which uses as contrast function the negentropy. The latter was approximated by the fixed-point FastICA algorithm [45].

ICA has been successfully applied in the field of volcano seismology [33,46,47,48,49].

3. Results

3.1. ICA Results

The application of ICA was carried out on homogeneous variables, i.e., on the five available CFc TIRs. To get robust results, the series were band-pass filtered in the range (1 year 3 days). The results of the decomposition are reported in Figure 3, where both the separated waveforms and the related power spectral densities (PSD) are shown. As can be seen, some characteristic peaks appear in broader spectra corresponding to specific periodicities, in agreement with those evidenced by spectral analysis reported in Caputo et al. [15].

Consequently, five independent components (ICs) were identified in the following non-overlapped bands: IC1 at (85–115) days; IC2 at (150–160) days; IC3 at (40–50) days; IC4 at (30–40) days. The last extracted component (IC5) has a more complex spectral behavior with two relevant peaks in the following bands: (10–20) days and (55–65) days.

3.2. Time Evolution of the Spectral Power

Based on the information provided by ICA, we estimated the energy contained in the main spectral peaks over the time. The PSDs were evaluated on 180-days-long windows, shifting each 30 days, implying a superposition of 83.3%. The spectral power was, then, calculated by integrating the PSD of the TIR in the following bands: (10–20) days; (30–50) days; (55–65) days; (85–115) days for all the stations. Notice that the periodicity around 150 days was not included in the present analysis due to the shortness of the time-series. The results are reported in Figure 4: the y-axis is almost the same for allowing comparison among the traces.

The overall behavior of TIRs at CFc seems to be independent of the specific band for each station. From a comparison among the different stations, one observes a sort of amplitude attenuation moving away from SOB. Looking at Vesuvius, the behavior at all the time scales is similar to those at CFc in the first part, up to the end of the year 2017 (the day 540 in Figure 4 corresponds to 16 December 2017), implying a common exogenous origin, changing drastically after that.

It is noteworthy that the 10–20 d-band stronger deviates from the overall patterns. In particular, in the time-range late 2019–end of observations, this component shows a spectral power level inside the Solfatara crater (SF1, SF2 and SOB) significantly higher than the other bands/TIRs. We recall that the uncertainty on the date is of about ±90 days, since every spectral power value was estimated over a range of 180 days.

3.3. Changing Points Analysis

To quantify the similarities that emerged in the spectral and ICA analyses, we first compared the TIR time-series of CFc with VES, and all the TIRs with the air temperature by using the cross-correlation function for each of the identified frequency bands. The results are summarized in

Table 1. Cross-correlation values of time-series.

XCorr	10–20 Days	30–50 Days	55–65 Days	85–115 Days
Ves vs. PIS	0.24	0.33	0.70	0.70
PIS vs. Tatm	0.23	0.49	0.49	0.70
Ves vs. Tatm	0.43	0.77	0.72	0.90
SF1 vs. Tatm	0.47	0.53	0.57	0.81
SF2 vs. Tatm	0.39	0.54	0.62	0.89
SOB vs. Tatm	0.28	0.45	0.39	0.73
OBN vs. Tatm	0.33	0.45	0.63	0.84

. In general, it emerges that the longer the period, the higher the correlation coefficient. The high correlation among the CFc TIRs, VES and the air temperature in correspondence in the longer periods could reflect the influence of large scale atmospheric factors on the ground temperature.

To deeper investigate the similarity among the time-series, we performed a Changing Points (CPs) analysis by using Matlab function findchangepts [50]. This function partitions the input time-series into adjacent segments that minimize the sum of the squared residual error of each segment from its local mean, and returns the changing points and, optionally, the trend inside each segment. The algorithm requires one to specify the maximum number of CPs to retrieve, as input parameter. To fix the maximum number of CPs, we let it vary over a wide range of values (from 1 to 20) and then evaluated how many output CPs are estimated. For each series, we found at most 8 CPs and, hence, we set the maximum number of CPs parameter to this value. We checked also the rms option, which evaluates the root-mean-square level of the most significant signal changes, and obtained the same results. The results are shown in Figure 5.

T_atm has 8 CPs, whereas for VES and PIS, we found 7 CPs, as reported in

Table 2. Changing Points of T_atm (CPTs), VES (CPVs) and PIS (CPP) time-series. The minus/plus signs indicate decreasing/increasing trends. The trend is referred to the linear trend between CPs, in °C.

	Tatm CPTs Date	Trend 1 (°C)	Ves CPVs Date	Trend 1 (°C)	PIS CPPs Date	Trend 1 (°C)
1	4 January 2017	+0.81	5 January 2017	+3.30	12 December 2016	+1.52
2	5 February 2017	min, 14.93	22 January 2017	+17.49	21 January 2017	max, 67.38
3	29 December 2017	−2.12	29 September 2017	−7.68	16 January 2018	−7.27
4	28 February 2018	max, 19.03	21 December 2017	−8.88	2 May 2018	+2.70
5	26 April 2019	−0.77	2 April 2018	−5.51	16 April 2019	stable, 0.00
6	5 June 2019	min, 14.52	4 September 2018	−3.50	20 May 2019	min, 59.84
7	22 March 2020	−0.11	2 July 2019	−6.16	10 December 2019	+1.55
8	12 May 2020	min 15.66		+0.68		−0.62
		−6.77

¹ starting from the beginning or the preceding CP.

and represented in Figure 5 (green lines). The linear trends variation in °C between the consecutive CPs is indicated in Table 2 as well. Almost all the time-series show one or two CPs in correspondence with the main air temperature minima: 5 February 2017; 5 June 2019; 12 May 2020. This behavior common to all could be ascribed to the sensibility of soil temperature to strong air temperature anomalies. For all the time-series, further support for that hypothesis comes from the stability of the first 3 CPs’ position as varying the CPs’ number parameter in the Matlab function.

Figure 5 highlights a common pattern of T_atm, VES, PIS, SOB and OBN from the beginning of the observations to the end of year 2017: (1) increasing until the first CP; (2) a minimum between the first and the second CPs; (3) decreasing. The slight differences among homologous CPs throughout the various time-series could be attributed to a different site response to a specific solicitation, due to local rheology and/or topography which modulates the geophysical and geochemical signals [33,34,35,42]. Taking into account the differences and distances among the instruments, the common pattern can be ascribed to an exogenous, large spatial scale effect.

The similarities among T_atm and P_atm and VES and OBN likely indicate a significant influence of external meteorological conditions on these TIRs. VES shows also an overall decreasing trend, in agreement with other observations [22].

Hereinafter, we indicate the CPs of T_atm, VES and PIS as CPT, CPV and CPP following the scheme of Table 2.

Interesting patterns are visible for SOB and, mainly, for PIS (Figure 5). Starting from the end of 2017, T_atm presents a maximum (29 December 2017, CPT-3) that ends in correspondence with CPT-4 (28 February 2018). Otherwise, in the same time range, PIS continued to increase, until CPP-4 (2 May 2018). From this CP, PIS continues to follow a different behavior. In particular, this TIR has a CP on 10 December 2019 (CPP-7) not present in T_atm nor in P_atm. Concerning SOB, even if the CPs of this TIR have the correspondent ones in T_atm, from the beginning of 2018, the time-series amplitude pattern is more like PIS. The difference between PIS and SOB on December 2019 (CPP-7) seems not driven by external factors.

A correlation analysis [51,52] of PIS and SOB TIR data acquired during 2016–2017 and 2018–2020 shows that at zero lag, the time-series are correlated well above the 95% confidence level (blue lines in Figure 6); in addition, the periodicity observed in 2016–2017 (Figure 6, top), and likely related to T_atm, strongly reduces in 2018–2020.

For SF1 and SF2 (Figure 5), we found overall flat trends, both interrupted by 8 and 7 CPs respectively, marking the minima and maxima that exactly correspond to the main minima and maxima of atmospheric temperature (Table 2). SF1 has an overall constant trend in correspondence of about 53 °C. SF2 has a linear slight increasing pattern, starting from about 64 °C to 66.4 °C. These trends are possibly attributable to the presence of several continuous heat sources inside the Solfatara crater, acting as a mechanical system in thermal equilibrium with a heat bath at a fixed temperature.

The described results are compatible with the existence of an overall exogenous large scale effect that significantly influences the ground temperatures. This effect dominates over the trends until the end of 2017. At the beginning of 2018, the CFc TIRs start to significantly deviate by this exogenous pattern (though still present), above all, in correspondence with exceptional long-lasting (several weeks) meteorological phenomena (main, minima and maxima).

CFc TIR Time-Series, Seismicity and CO₂ Flux

We compared the temporal pattern of the TIR PIS with the cumulative seismic moment, which is proportional to the seismic energy [53], and CO₂ flux [22], to check a possible link with endogenous sources.

The seismic moment was obtained by applying the relationship [54] for the CFc area:

LogM₀ = 9.9 + 0.9·Md

(1)

where M₀ is the seismic moment expressed in Nm and Md is the duration magnitude.

As shown in Figure 7, the seismic energy/moment increased starting from 2018. The most significant changes in the rate of energy release occurred in December 2019 and April 2020 in concomitance with Md 3.1 and 3.3 earthquakes. The interval of increasing seismic energy corresponds to the TIRs’ deviation from the atmospheric trend, suggesting that since 2018, TIR patterns could be more affected by endogenous dynamics. Indeed, since 2018, signs of variations were also detected in earthquake number and location. Bellucci Sessa et al. [20] observed an increase of both seismic swarms and background seismicity, particularly marked in 2019. In this last year, the location became shallower, with hypocenters in the first 2 km of depth. Moreover, the authors evidenced a progressive spatial clustering of the seismicity, started in 2016 (and even more evident in 2019), inside the area of Solfatara–Pisciarelli.

In a similar manner, the temporal pattern of PIS temperature is correlated with the CO₂ flux released at Pisciarelli, measured using the accumulation chamber method [42]. The flux is continuously acquired by a multiparametric station, employed 20 m downwind of Pisciarelli main boiling pool. Starting from the beginning of 2018, the CO₂ flux showed a growing oscillating trend, and presented a marked peak in early 2020 [22].

4. Discussion and Conclusions

In this study, we have analyzed time-series over about four years (25 June 2016–29 May 2020) of the surface temperatures extracted from thermal infrared images acquired by the TIR stations operative in the SolfataraPisciarelli volcanic complex. Investigations regarding the thermal behavior of this area, characterized by the main hydrothermal/volcanic emission of CFc, are aimed to contribute to the identification of possible endogenous and/or exogenous sources. The analyses were carried out by comparing the patterns of CFc thermal infrared surface temperatures, with the air temperature and atmospheric pressure measured at Capo Miseno meteo-station. In order to find eventual large-scale influence on TIR trends, a comparison was also carried out with the TIR station at Mt. Vesuvius. CFc TIR time-series are also compared with the cumulative seismic moment release of the CFc seismicity and the CO₂ flux at the Pisciarelli site.

Although seasonal periodicity was removed from all the time series, atmospheric factors clearly still influenced the surface temperatures. Both the spectral-ICA and CP analyses have demonstrated the existence of a non-negligible impact of atmospheric effects. The spectral content of the estimated ICs, which is common to the CFc TIRs, matches the main spectral content of the atmospheric time-series and of VES TIR, in agreement with Caputo et al. [15]. For all the TIRs, we followed the time-evolution of the spectral bands identified by the application of ICA. For all the bands, we found a common pattern for both CFc and Mt. Vesuvius, until late 2017. Then, the CFc TIRs started to differentiate from the VES series’ trend, although some common features were still present. The cross-correlation analysis has shown similarities among the TIRs: the longer the period, the higher the correlation level, in agreement with the influence of a large scale factor.

The common patterns observed in 2016–2017 were confirmed by the CPs analysis, in which we found compatible CPs’ positions and trends among all the TIRs, except for SF1 and SF2. The latter TIRs had, indeed, almost flat trends all over the observation time, suggesting a surface temperatures field inside the Solfatara crater due to a system in thermal equilibrium.

The footprints of a background exogenous large scale influence were mainly testified by the coincidence of the CPs of most series, with the maxima and minima of T_atm, which in turn corresponded to exceptional long-lasting (several weeks) meteorological phenomena.

The interplay between the atmosphere parameters and the ground temperatures is nontrivial. Behind the seasonal periodicity (that we have corrected), the influence of the atmospheric factors on the surface temperatures could be twofold: a direct effect, linked to ground heating and/or the atmospheric loading; and an indirect effect on the flux of hydrothermal fluids, such as CO₂ flux, which in turn affects the surface temperatures [42]. Indeed, the link between atmospheric exogenous effects and gas emissions has been also noticed for CFc and Mt. Vesuvius [42,55,56], as well as in other worldwide cases [57,58].

At the beginning of 2018, while T_atm values decreased after reaching a maximum, PIS temperature continued to increase. Since then, PIS has deviated from the atmospheric trend. The different behavior of PIS could not be simply attributed to a local site effect, since SOB, which points towards the top of the SE slope of the Solfatara crater, shows a pattern similar to PIS. The year 2018 was also peculiar for seismic activity, since an increase of both seismic swarms and background seismicity were observed [20]. According to Bellucci Sessa et al. [20], in this period, seismicity began also to cluster in the study area. Moreover, a consistent fluid transfer was detected starting from 2018 by using polarization analysis of seismic noise [59]. In addition, CO₂ flux measured at the Pisciarelli site showed some anomalies in the flow rate [22]. These evidences strongly support the hypothesis of a common endogenous source to explain the patterns observed in the time span between 2018–2020.

Another strong indication of an endogenous driving of the ground temperatures in the Solfatara–Pisciarelli area appeared in December 2019, when we recorded a significant temperature decrease at PIS (about −4 °C in 20 days). At the same time, the spectral power time-evolution of Solfatara TIRs (SF1, SF2 and SOB) showed an increase of soil temperature in the 10–20 d band. Simultaneously (6 December 2019), a remarkable seismic swarm occurred beneath the Pisciarelli area, causing one of the most significant increases of seismic energy rate of the last years. From this time, the earthquake locations became shallower, with hypocenters in the first 2 km of depth [20]. A significant rising of CO₂ concentration was also recorded after the December 6 swarm, superimposing to an already growing trend, suggesting a consistent increase of the emitted fluid flux at Pisciarelli boiling pool. No evident variation of CO₂ flux was evidenced inside the Solfatara crater [22].

The impact of the hydrothermal/volcanic source over the surface temperatures is also evidenced by the correlation analysis performed for PIS and SOB in 2018–2020. The correlation analysis results are also in agreement with the geophysical/geological models of the shallow hydrothermal system of the Solfatara–Pisciarelli complex reported in recent studies [60,61,62,63]. These studies evidence the presence of a resistive clay-rich cap rock extending underneath the study area, at depths from 0 to 50 m b.s.l. This cap rock is interrupted by conductive gaps, due to fracture zones (where subvertical faults exist). The fractures allow the pressured fluids to up-flow from a deep hydrothermal system to reach the surface, where the highest values of gas flux are recorded, in particular, at Solfatara SE rim fumaroles and Pisciarelli boiling pool. On the other hand, underneath the Pisciarelli area, a shallow moderate resistive zone is observed, interpreted as a zone where the outflow occurs, deepening toward the region where most of the recent seismicity is located.

In conclusion, the results of the analysis performed in the study area allowed one to identify the time periods in which endogenous hydrothermal/volcanic sources strongly influenced TIR surface temperatures from those periods in which exogenous factors prevail. Concerning the endogenous driving, the most interesting patterns occurred at the end of 2019 at the Pisciarelli area. We hypothesize a simple conceptual model to explain the evidenced patterns. The seismic swarm that occurred in December fractured the medium, increasing the permeability and, likely, opening a preferential uprising path for the pressurized fluids towards the boiling pool at the Pisciarelli emission area, where the CO₂ flux is routinely measured. The ejection of the gases could have caused a depressurization of the hydrothermal/volcanic system, which in turn caused the lowering of the ground temperatures.

On the other hand, in the Solfatara crater, the fluid circulation remained almost unchanged after the swarm; thus, a possible moderate fluid injection caused the growth of the ground thermal energy rate.

The present study has confirmed the important contribution of TIR observations in monitoring hydrothermal/volcanic system. TIR datasets are, indeed, useful to detect the effect of the endogenous source when the hydrothermal/volcanic forcing becomes strong enough to overcome the atmospheric effects. Moreover, the use of a network of TIR stations is fundamental to exclude local site effects.

Nevertheless, we stress the importance of monitoring the atmospheric parameters in hydrothermal/volcanic environments, since the ground temperatures, the gas flux and the seismicity patterns have been demonstrated to share evident correspondences with atmospheric behaviors (especially wet–dry seasons’ periodicity) [40,42,56].