Anomalous temperature regimen in the near-surface soil layer of Tlamacas hill and its relation to activity of Popocatépetl Volcano, Mexico

– This article describes anomalous changes in the diurnal behavior of the temperature measured in the near-surface soil at the Tlamacas monitoring site, Popocatépetl. Results of the statistical analysis show two essential changes for the temperature characteristics observed during the 2007 – 2009 (quiet volcano) and 2013 – 2014 (active volcano) monitoring periods. Under normal conditions, the absolute minimum daily temperature is observed at about 7:40 Local Time (LT) during sunrise for the atmosphere and, with a time delay, at about 8:30 LT, for soil measurements. The absolute temperature maximum is observed about 15:30 LT for the atmosphere and 16:30 LT for in-soil measurements. The dispersion of the residual temperature (24-h running trend of the temperature substituted) is 5.6 times lower for the 2013 – 2014 period in comparison with the 2007 – 2009 period. In other words, in 2013 – 2014, the temperature variability became 5.6 times lower that it was in 2007 – 2009.


Introduction
Monitoring of soil gas temperature, meteorological parameters, radon concentration and certain other associated physical characteristics is an important part of volcanology. Some previous investigations focused on soil temperatures in volcanic areas such as Stromboli, to investigate complex magmatic processes (Brusca et al., 2004;De Gregorio et al., 2007;Cigolini et al., 2009;Revil et al., 2011). Cigolini et al. (2009) showed the results of the combined monitoring of 222 Rn concentration, atmospheric pressure, and soil temperature since April 2007 (after the eruptive crisis of February-April 2007) until May 2008. In general, the variation of temperature seems unremarkable, reaching maximum values in July-August and minimum values in December-February. However, a more detailed analysis of these data is warranted because the maximum temperature corresponds to the minimum pressure and radon concentration. The authors explain this phenomenon as being due to hydrothermal convection (Cigolini et al., 2009). The origin of hydrothermal convective flux variation in the Fossa of Vulcano (Italy) is described in the work of Aubert (1999), Aubert and Alparone (2000) and Aubert et al. (2008). Nick Varley and co-authors describe an explosion activity at Volcán de Colima derived from the rapid release of pressure and compare this process with other volcanoes, such as Santiaguito, Guatemala (Johnson et al., 2008), Sakurajima, Japan and Semeru, Indonesia (Iguchi et al., 2008) and Popocatépetl, Mexico (Varley et al., 2010).
There are different classifications of eruptive activity. One of these is by the volcanic explosivity index (VEI) (Hickson et al., 2013). According to this index moderate eruptions correspond to VEI 2-3 (volume of erupted tephra about 0.01 km 3 ), although major eruptions correspond to VEI 3-4 (volume of erupted tephra about 0.1 km 3 ) and above (Newhall and Self, 1982). Possible scenarios of activity at Popocatépetl could be defined as: small to moderate-size explosions (Alatorre Ibargüengoitia, 2011) some with ash emission; occasionally slight incandescence in the crater observable during the night with a chance of the expulsion of incandescent fragments close to the crater. Since 1994, the volcano has had almost continuous low magnitude activity consisting of small explosions, sometimes with ash. These events occur with a frequency of several to tens of eruptions observed during relatively quiet volcano phases, and up to 50-100 eruptions occurring during a more active phase (González-Pomposo, 2004;Arámbula-Mendoza et al., 2010). Eruptions at Popocatépetl can be preceded by changes in seismic and volcano-magnetic activity, chemical composition of the gases and spring water and sometimes deformation (Martín-Del Pozzo, 2012). A number of reports summarize activity at Popocatépetl during December 2007-December 2009 (Report on Popocatepetl, 2007, 2009, 2012) and June 2013-March 2014(Report on Popocatepetl, 2015. The data as usual came from online daily reports by the Centro Nacional de Prevención de Desastres (CENAPRED). Moderate ash-rich explosions take place from one in several months up to several times per day for quiet and active phases, respectively.
Moderate eruptions can also be accompanied by ballistics, with a frequency of such event varying from several to tens of events per year. Major eruptive activity is not so frequent.
Our earlier studies of the radon emanation in the volcano Popocatépetl area (Fig. 1, latitude 19.07°N, longitude 98.63°W, altitude 5465 m above sea level (a.s.l.)) showed a variety of anomalies in the radon behavior associated with the volcano activity (Kotsarenko et al., 2012). The previous monitoring of soil radon released at three different volcanic sites showed nine cases of decreasing radon concentration when approaching volcanic eruptions for a total of 23 moderate explosions. Surveys of soil radon and gamma ray spectrometry revealed the intensive emanation of radon and high gamma radiation in the area of Tlamacas; the levels of the radon concentration in the Tlamacas area (latitude 19.02°N, longitude 98.63°W, altitude 4000 m a.s.l.) were 10-20 times higher than background volcano values.
Similar radon concentration monitoring and surveys at volcano-tectonic geological structures have been carried out in different active volcanoes all over the world (Martín et al., 2003;Burton et al., 2004;Hernandez et al., 2004;Londoño, 2009). The current study is devoted to the analysis of the temperature variations measured during the radon monitoring (June 2013-May 2014) at the Tlamacas observation site, situated 4 km north of the volcano crater. It coincides with a significant increase in activity of Popocatépetl since April 2012. Results of our statistical analysis reveal the thermal anomaly in this area and we present a qualitative and quantitative description of the observed phenomena.

Geological setting
Popocatépetl is a Pliocene-Quaternary stratovolcano situated 65 km SE of Mexico City and 45 km to the west of the city of Puebla, in the frontal part of the Trans-Mexican Volcanic Belt (TMVB) on a basement of Paleozoic metamorphic and Cretaceous sedimentary rocks ( Fig. 1A) (Siebe et al., 1995a(Siebe et al., , 1996Macías, 2005; Espinasa-Pereña and Martín-Del Pozzo, 2006;Arana-Salinas et al., 2010;Sosa-Ceballos et al., 2012, Padilla y Sánchez, 2017. El Popo was constructed over the remains of a volcanic paleostructure (Sosa-Ceballos et al., 2015) and forms the southern end of the Sierra Nevada volcanic range. Its stratigraphic column consists of pyroclastic deposits (sandy ash, pumice, and ash flow deposits, containing lithic clasts of granodiorite, hornfels, arenite and other xenoliths) erupted during the past 23 000 yr. BP (Siebe et al., 1995b;Siebe and Macías, 2006). Mooser et al. (1996) reported the existence of three regional units that are differentiated by their lithology. The first corresponds to the volcano Popocatepetl, composed of andesitic, dacitic and rhyolitic lavas, alternating with pyroclastic materials. The second unit corresponds to the foothills of the volcano, and is composed of pyroclastic layers as ash and pumice, plus lahar, fluvial and fluvio-glacial deposits. This whole set is called the Tarango Formation. The third unit consists of a different type of monogenetic volcanoes that are located in the SW sector of volcano on its piedmont. These form part of the Chichinautzin Group, where cones of scoria with extensive lava flows of basic and intermediate composition, i.e., basalts and andesites, are widespread.

Data collection and processing
The in-soil temperature measurements were carried out using portable radon detectors SARAD Scout Plus (the primary purpose of which is monitoring of 222 Rn concentration).
The detectors were equipped with simple sensors for temperature, air humidity, atmospheric pressure, and an accelerometer. The devices were installed at the monitoring sites in the soil at a depth of 30-40 cm and operated during a period of 1-2 months before replacement. Deeper installations were not possible because of the presence of solid rock below this depth at the majority of measurement sites. The instruments operated with an exposition time from 1.5 to 2 h, for earlier detector versions of SARAD Scout (the internal memory size was limited to a series of 1023 measurement points) and from 30 min to 1 h for a later version (Scout Plus; 2047 pts).
Two statistical values, "averaged daily temperature" and "residual temperature", were used for subsequent analysis. The averaged diurnal temperature was calculated as follows. Interpolated temperature data were grouped into a matrix N days Â 24 h (the first and last segments representing incomplete 24 h were rejected), arithmetic daily means and their variation (RMS) were calculated. The average daily temperature clearly showed the appearance of a temperature anomaly during the active volcano period (since 2012) compared with a calm one.
The residual temperature is a quantitative temperature value that shows the difference between the maximum and minimum temperature during the day.

Results
The diurnal character of the reference temperature variation measured at a typical medium-altitude (1900 m) meteorological station Juriquilla, Querétaro, averaged for the period of Dec 2007-Dec 2009 is presented in Figure 2a. During the night, in the absence of solar radiation, the temperature continuously decreased until the time of sunrise, so that the temperature minimum is observed about 7:40 LT. After that the temperature monotonically increased reaching its maximum at about 15:30 LT, after which the inclination of the Sun was such that solar radiation could no longer compensate the cooling of both the Earth's surface and its atmosphere.
The character of the average temperature variation at the Tlamacas site for the period of a relatively low activity ( (Levresse et al., 2019). There are slight differences due to the altitude of Tlamacas (4000 m. a.s.l.) and the type of temperature measurements (in-soil method). When the measurements were made in the near-surface soil at a depth of 30-40 cm, we have observed that the soil temperature lags behind the temperature in the open air and, therefore, the minimum and maximum of the soil temperature occurred later, about 8:30 and 16:30, respectively. As well the temperature variability in Tlamacas was much less in comparison with Juriquilla ( Figs. 2a and 2b).
However, since June 2013 when radon monitoring at Popocatepetl was restarted and after the volcano had begun a more active phase (since April 2012), the character of the diurnal temperature variation had changed significantly. The regular diurnal temperature anomaly at Tlamacas was detected, as one can see in Figure 2c (see also Figs. 3 and 4). The maximum of the daily temperature occurred at about 01:00 LT and its minimum at about 14:00 LT.
As described above, the temporal difference between the minimum and maximum daily temperature may vary by ± 1 h for different months of observations or even several hours as was detected during August 2013 (Fig. 3, lower graph). In this way, one can see that a qualitative change in the daily character of the temperature variation at the Tlamacas observation site is neither an instrumental error nor an instrument setup mistake; the phenomena is clear and was regularly observed at all observation periods starting from June 2013. This anomaly was detected using different SARAD Scout detectors.
Another distinctive feature of the observed temperature anomalies for the period-heightened activity is the reduction of the temperature variability. The temperature variability, i.e. the difference between the maximum and minimum temperature, or, in terms of statistics, the amplitude of residual temperature (Fig. 5c) and its dispersion (Fig. 5d), decreased 5.6 times when compared with values for the quiet volcano (Figs. 5a and 5b).

Thermodynamic modeling
The origin of the observed temperature anomaly is still unclear. However, using thermodynamic concepts, the described temperature behaviour can be explained by the presence of a powerful mechanism of heat release, dominating over the natural processes of heating (solar radiation) and cooling of both the Earth's surface and the sub-surface soil layer. Below we present a thermodynamic model of the observed data within the box in which the sensor was placed.
The one-dimensional thermal conductivity equation is used (American Institute of Physics Handbook, 1972): Here T is the temperature of the air within the sensor box, t is time, x is a length coordinate (0 < x < L), and K is the temperature conductivity coefficient of the air. On the right side of equation (1) there is a localized heat source; the function c(x) describes the spatial distribution of the heat source, h(t) is its dependence on time. In addition, there were provided the simulations within a complex three-dimensional model, but the results were qualitatively the same as in this simplest model.
The coefficient K is equal to K = l/(rc v ), where l ≈ 2.2.10 À2 W/(K.m) is the thermal conductivity coefficient of the air, r ≈ 1 kg/m 3 is the density of air, and c v is the heat capacity per 1 kg (American Institute of Physics Handbook, 1972). In turn, c v = (5/2) R/m, where the universal gas constant is R ≈ 8.34 J/(K.mol) and m = 29 g/mol is the molar mass of the air. The estimation for K results in K ≈ 0.3 cm 2 /s; its value changes weakly within the temperature interval 0-30°C.
Below the thermal processes are considered within the sensor box for lengths of several centimeters, l n ∼ 2 cm, whereas the temporal variations are of about several hours, t n ≥ 1 h ≡ 3.6 · 10 3 s. At this scale, the temperature conductivity coefficient is large, and the inertia of the thermal processes can be neglected: K tn /l n 2 ∼ 300 >> 1. Therefore, equation (1) can be simplified: The following boundary conditions apply to: Here ' 1 (t), ' 2 (t) are functions of time. The problem is to reconstruct the dependence of the heat source on time h(t), when the temperature profiles ' 1 (t), ' 2 (t) are known at each end of the system x = 0 and x = L, equation (3), jointly with the measured temperature within the system at x = x*: Also the spatial profile of the heat source c(x) is known. It is assumed to be localized within the system: The problem of the reconstruction can be solved analytically in this simplest treatment. The solution of equation (2) with the given boundary conditions (3) is: Using equation (4), one obtains the following formula for h(t): The reconstruction of the heat source has been performed for the input data presented in Figure 6a. The temperatures at the boundaries are the same: ' 1 (t) = ' 2 (t), as presented in Figure 2b. The temperature c(t) is measured in the centre of the box x * = L/2 = 5 cm.
The result of the reconstruction is given in Figure 6b. The size of the system is L = 10 cm. The heat source is localized in the center of the system x 1 = L/2 = 5 cm and its size is x 0 = 1 cm. The results of the reconstruction are tolerant to the variations of the localization point of the source x 1 , its shape, and its size x 0 . The positive values of h(t) correspond to heat release, whereas the negative values are due to heat absorption. One can see that the experimental temperature profiles can be explained by the heat release within the time interval 21:00-13:00 and the heat absorption during the time interval 13:00-21:00. The total heat release dominates over the total heat absorption.
In the case of neglecting the non-stationary term ∂T/∂t, an influence of other mechanisms of heat transfer, such as thermal convection does not change the results, because it leads only to an increase of the coefficient K. If the thermal processes that release the heat occur at greater spatial scales l n > 20 cm, then the reconstruction should be realized in the general case of equation (1) with the term ∂T/∂t, but with numerical methods only (Samarskii and Vabishchevich, 2007).

Discussion
According to the regularity of the diurnal temperature variation, in the upper part of the sub-surface layer at a depth of 30-40 cm below the surface of the Earth (the temperature sensor immersion depth) the temperature must increase with a certain time delay, after the sunrise to a maximum at approximately the second half of the day (16:00 local time) and then gradually decrease to a minimum near 8:00 local time. Our measurements have shown that this temporal distribution corresponds to the temperature data obtained within a nearby tectonically-passive area (Juriquilla, Queretaro) and within the Popocatepetl area during the inactive volcanic phase (Figs. 2a  and 2b). However, for the active volcanic period, we identify a discrepancy of this variation, whereby there is not just a simple heating/cooling of the material due to the daily solar variation (Fig. 2c). This discrepancy induces three aspects (Figs. 2b and 2c): (1) the displacement of the daily temperature minimum from 8:00 to 14:00-16:00 (LT), i.e., at the time of the usual maximum heating by the Sun; (2) substantially reducing the daily temperature variations from 3°C to 0.5°C; (3) as well as a rise of the temperature as a whole by approximately 0.75°C. As we know that during these periods of time (December 2007-December 2009and June 2013-March 2014 the solar activity (external source of soil heating) showed negligible variation and the daily temperature was regular, we need to seek another explication for the observed phenomenon, e.g. internal sources. Below we discuss this problem with reference to the papers mentioned above and others.
Paleostructure, eruption dynamics and magma mixing, geochronology and geochemistry of El Popo have been discussed in different works (Sosa-Ceballos et al., 2012. These studies show the possible conduit of heat transfer, which could serve as the key to understanding the observed temperature phenomenon. Similar studies have been realized at Volcán de Colima (Varley and Johnson, 2005;Varley et al., 2006Varley et al., , 2010Arámbula-Mendoza et al., 2018) There are different ways of transporting heat: conduction, convection, radiation, and advection. Heat flow and the geothermal gradient are the most important concepts for understanding temperature effects in the Earth's interior. Geochemical evidence for a mantle origin, plus crustal processes in volcanic rocks from Popocatepetl and surrounding monogenetic volcanoes was described by Schaaf et al. (2005). There are many studies concerning the geothermal regime below volcanic areas (e.g. Varley and Johnson, 2005;Pasquale et al., 2014;Salzera et al., 2017). The geothermal gradient in the upper crust in those areas has been the focus of many studies and is characterized by maximum values up to 25-30°C/km (Dobretsov et al., 2001). Although there are different components controlling the geothermal gradient (radiometric and conductive in particular), the thermal evolution of planetary crusts and lithospheres is largely governed by the rate of heat transfer by conduction (Turcotte and Schubert, 2001;Stacey and Davis, 2008). Whittington et al. (2009) showed that the governing physical properties are thermal diffusivity (к) and conductivity (k = кrC P ), where r denotes density and C P denotes specific heat capacity at constant pressure. Vosteen and Schellschmidt (2003) and Mottaghy et al. (2008) revealed that for crustal rocks both к and k decrease above ambient temperature. On the other hand, most thermal models of the Earth's lithosphere assume constant values for k(∼ 1 mm 2 s À1 ) and/or k(∼ 3 to 5 Wm À1 K À1 ) owing to the large experimental uncertainties associated with conventional contact methods at high temperatures (Michaut et al., 2006). Usually, temperature gradients within about 20 m below the surface are strongly dependent on diurnal changes in solar heating (Kearey et al., 2005), but in active volcanic areas, heat flow (the rate of heat energy transfer through a given surface per unit time) has an important role in the sub-surface temperature variations. The rate of heat flow is proportional to the difference in heat between two bodies and may be expressed by: where Q is a heat flow, DT is the temperature difference between the bottom and top of the block, L is a length of the block, A is an area of its cross section, and l is the value of thermal conductivity. We propose that the observed temperature variation in the Tlamacas area is strongly related to temperature variations in Popocatepetl's magma chamber. The magma chamber, which is the source of heat in the Popocatepetl zone, is located at a depth of 5-6 km below the Earth's surface (Sosa-Ceballos et al., 2012. It is filled by dacite magma and has a temperature of about 860°C (Macías, 2005(Macías, , 2007. Hildreth (1981) showed that for magma chambers there are some indications of a temperature gradient of the order of 200°C/ km. This temperature gradient corresponds to a vertical heat flux of 10-15 HFU (HFU = 10 À6 cal cm À2 s À1 = 41.84 mW/m 2 ) (Spera et al., 1982). This calculation does not mean that the  temperature gradient is equal all the way from the magma chamber to the Tlamacas sub-surface, but it may be one of the reasons for the increasing sub-surface temperature. Moreover, a one-dimensional thermal conduction model simulates the repetitive intrusion of basalt sills into the deeper parts of the crust (Annen and Sparks, 2002). So we propose that the abovementioned process, which takes part at the magma chamber boundary, has a periodic change of intensity and direction, i.e. the increasing/decreasing of temperature. Periods of magmatic intrusion create reverse geothermal gradients and thermal anomalies in the crust. Karlstrom and Richards (2011), describing the evolution of large magma chambers for flood basalt eruptions, simulated the evolution of melt within and elastic deviatoric stresses surrounding a deep magma chamber (Karlstrom and Richards, 2011). They assumed an idealized temperature evolution using the equation: where T 0 is the initial temperature of the wall rocks, DT is the temperature change imposed by the magma chamber, R c is the radius and r is distance from the of the sphere, к = 0.5 Â 10 À6 m 2 /s (Whittington et al., 2009) is the thermal diffusivity and t is time. Then they simulated the storage and crystallization of this melt in a magma chamber. Cooling simulations are performed isochorically and isobarically, as end-member fractionation scenarios (Fowler and Spera, 2008). Another way of heat transport is radioactivity. The continental lithosphere contains significant amounts of radioactive elements (U, Th, K) produce heat. In the nearsurface layers, the radiogenic heat flux component is determined by the concentration of radioactive elements. These elements are concentrated in magmas and are carried into the upper crust. The radiogenic heat in the investigated area is closely associated with high values of uranium (Uup to 4.6 ppm) and thorium (Thup to 14-15 ppm) in the mineral composition of the rocks and sediments (Kotsarenko et al., 2012). We carried out gamma ray measurements at the same locations as the temperature determinations using the differential gamma ray spectrometer. The cone of the volcano complex consists of pyroclastic rocks and lava flows of andesite-dacite composition (Schaaf et al., 2003). Their fragments have crystals of plagioclase, hypersthene, augite, olivine, etc. in a glassy microcrystalline matrix. The content of uranium in these volcanic rocks is about 3.5-4 ppm, which approximately corresponds to our data (up to 4.6 ppm). However, our measurements revealed that the thorium content in the Tlamacas near-surface layer reaches up to 14-15 ppm, which is significantly higher than the average level (9-11 ppm). So it is possible that the high level of near-surface radioactivity is associated with the ancient volcanic activity of Popocatepetl. Unfortunately, we did not perform gamma ray monitoring during the study period.
From the analysis of the components of thermal energy in the sub-surface layer of the Tlamacas site, it follows that the observed phenomena of temperature variations during the active and passive phases of Popocatepetl is dominated by the heating/cooling of magma in the magma chamber. Both heating and cooling processes have a strong relation with the radioactive and magmatic generation of heat. During heightened activity the temperature generation is higher due to the heat flux from the magma. However, the influence is apparent when considering the diurnal temperature amplitude, which is relatively low during the active volcanic phase, not more than 0.5°C. Even though melt generation takes place near the mantle-crust boundary, and the magma chamber does not generate more heat during the active volcanic phase, the whole conduit, including fracture zones, becomes more active during this period, and is able to transfer heat faster. Besides, an ascending column of magma in the conduit produces heat.
On the other hand, during the inactive phase, when the solar energy plays a dominant role in heating/cooling of the soil, these values reach 10-12°C, and their variation is regular: maxima at midday and minima during the night hours.
Our model offers hypotheses for the temperature variation in the sub-surface layer of Tlamacas hill: The main source of the heat generation in the near-surface layer during the active volcanic phase is often due to the heat released by the condensation of the gas phase (water vapour and air). This is the most common exothermic process (Diliberto, 2017), and the heat generation depends on the heat transfer manner, which is a complex combination of conductive and mostly convective.
The Solar thermal energy in this period is a secondary factor. This period is characterized by higher temperatures, lower amplitude between diurnal minimum and maximum, and as a result, displacement of the regular diurnal variation.

Conclusions
In this work we described the unusual phenomena of temperature anomaly variation observed in the area of Popocatepetl. The observed anomaly is very likely related to the activity of the volcano, because it was never observed during the period of volcanic quiescence. The described temperature behaviour can be explained by the presence of a powerful mechanism of heat release, i.e. a source of input and output of the energy within the sensor, dominating over the processes of heating (Solar radiation) and cooling of the Earth's surface and the sub-surface soil layer. The heat transfer process is a complex combination of conductive and mostly convective process, because the partial pressure of the gas in porous ground is high. A physical mechanism for the temperature anomaly remains under discussion, but almost certainly, the governing processes are complex in character and for the interpretation more geophysical studies on the presented phenomena (and probably laboratory experiments) are required.