Thermodynamic and hydrological drivers of the soil and bedrock thermal regimes in central Spain

An assessment of the soil and bedrock thermal structure of the Sierra de Guadarrama, in central Spain, is provided using subsurface and ground surface temperature data coming from four deep (20 m) monitoring profiles belonging to the Guadarrama Monitoring Network (GuMNet) and two shallow profiles (1 m) from the Spanish Meteorology Service (Agencia Estatal de Meteorología, AEMET) covering the time spans of 2015–2021 and 1989–2018, respectively. An evaluation of air and ground surface temperature coupling showed that soil insulation due to snow cover is the main source of seasonal decoupling, being especially relevant in winter at high-altitude sites. Temperature propagation in the subsurface was characterized by assuming a heat conductive regime by considering apparent thermal diffusivity values derived from the amplitude attenuation and phase shift of the annual cycle with depth. This methodology was further extended to consider the attenuation of all harmonics in the spectral domain, which allowed for analysis of thermal diffusivity from high-frequency changes in the soil near the surface at short timescales. For the deep profiles, the apparent thermal diffusivity ranges from 1 to $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹, which is consistent with values for gneiss and granite, the major bedrock components in the Sierra de Guadarrama. However, thermal diffusivity is lower and more heterogeneous in the soil layers close to the surface (0.4–$\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹). An increase in diffusivity with depth was observed that was generally larger in the soil–bedrock transition at 4–8 m depth. The outcomes are relevant for the understanding of soil thermodynamics in relation to other soil properties. Results with the spectral method suggest that changes in near-surface thermal diffusivity are related to changes in soil moisture content, which makes it a potential tool to gain information about soil drought and water resource availability from soil temperature data.

1 Introduction

Recent decades have witnessed substantial steps forward in knowledge of the climate system and ongoing anthropogenic climate change (Chen et al., 2021). An increase in observational datasets (Hansen et al., 2010; Lawrimore et al., 2011; Jones et al., 2012; Rohde et al., 2013; Lenssen et al., 2019; Osborn et al., 2021) and progressively more realistic climate model ensembles (von Storch, 2010; Taylor et al., 2012; Eyring et al., 2016; Maher et al., 2019) allowed the attribution of a global mean near-surface air temperature rise of more than 1 ^∘C since the end of the 19th century to human activities (Chen et al., 2021; Gulev et al., 2021). This warming has been larger over land than over the ocean, with temperature values in 2011–2020 being 1.59 and 0.88 ^∘C, respectively, higher than 1850–1900 (Chen et al., 2021). It is important to understand how this temperature rise has propagated into the soil and subsoil, which was shown to be important for land surface processes and meteorological extreme events (Miralles et al., 2019; Zhou et al., 2019; Seneviratne et al., 2021).

Ground surface temperature (GST) and surface air temperature (SAT) are coupled at multi-decadal and centennial scales (Pollack and Huang, 2000; Beltrami and Kellman, 2003; Melo-Aguilar et al., 2018), and GST perturbations are subsequently propagated into depth by heat transfer (Carslaw and Jaeger, 1959). A warming of near-surface soil has an impact on its biogeochemical activity (Soong et al., 2021), enhancing the metabolic activity of microorganisms, mineral weathering, and decomposition of soil organic matter (Schlesinger and Emily, 2013). It also affects soil hydrology (Krakauer et al., 2013) and increases soil respiration (Pries et al., 2017), contributing to an acceleration of the terrestrial carbon cycle (Bond-Lamberty and Thomson, 2010; Pries et al., 2017). Moreover, soil warming has relevant impacts on high-latitude regions, where permafrost thawing can produce changes in soil hydrology (Andresen et al., 2020; Burke et al., 2020) and eventually subsurface collapse (Abbott and Jones, 2015; Pegoraro et al., 2021), posing risks to infrastructure (Rotta Loria, 2023) and reshaping landscapes. Moreover, permafrost thawing triggers the release of pollutants (e.g., Schaefer et al., 2020) and massive amounts of carbon (Turetsky et al., 2019) into the atmosphere, potentially affecting the global climate (e.g., de Vrese et al., 2023). The subsurface is sensitive to not only present global warming but also to past climate changes. Deep subsurface temperatures allow recovery of past GST histories by borehole temperature profile (BTP) inversion techniques (Mareschal and Beltrami, 1992; Huang et al., 2000; Pollack and Smerdon, 2004; Jaume-Santero et al., 2016; Cuesta-Valero et al., 2022), which makes them a valuable source for climate reconstruction as far as the two fundamental hypotheses are sustained: SAT and GST changes are coupled at long-term scales and GST changes are propagated downwards by heat conduction, the main heat transfer mechanism in solids.

Support for SAT–GST coupling has been found so far from observational temperature data at some sites in North America (e.g., Putnam and Chapman, 1996; Beltrami and Kellman, 2003; Smerdon et al., 2004; Bartlett et al., 2005, 2006; Smerdon et al., 2006) and Europe (e.g., Cermak et al., 2014; Cermak and Bodri, 2016; Melo-Aguilar et al., 2022; Petersen, 2022; Dorau et al., 2022) as well as from climate simulations (González-Rouco et al., 2006; Smerdon et al., 2009; García-García et al., 2016, 2019; Melo-Aguilar et al., 2018; González-Rouco et al., 2021). These studies showed a strong coupling at long-term scales globally, with some degree of seasonal decoupling due to snow cover (Beltrami and Kellman, 2003) or changes in the surface energy balance (Smerdon et al., 2004; Melo-Aguilar et al., 2018) at regional scales. Examining whether conduction is the main heat transport mechanism in the subsurface was usually tackled by assessing the propagation of temperature perturbations and verifying that they experience an amplitude damping and phase shift with depth to an extent controlled by thermal diffusivity, a parameter that governs heat conduction processes (Carslaw and Jaeger, 1959). Previous research assessed the propagation of temperature trends (Lesperance et al., 2010) and thermal orbits (Smerdon et al., 2009) using a priori defined constant values of thermal diffusivity. Other studies tackled the propagation of single-frequency harmonic signals with depth, i.e., the annual cycle (Hurley and Wiltshire, 1993; Smerdon et al., 2003), in order to derive estimates of the apparent thermal diffusivity. Nevertheless, this approach requires continuous monitoring of subsurface temperature for at least a relatively long lapse of time (at least a few years). This was only achieved in a low number of places in North America (Smerdon et al., 2004; Bartlett et al., 2006) and Europe (Smerdon et al., 2004; Demetrescu et al., 2007; Dědeček et al., 2013; Melo-Aguilar et al., 2022). Moreover, the analysis might be hampered by potential deviations from the conductive regime in the soil near the surface due to biological, chemical, and hydrological processes (Gao et al., 2008; Tong et al., 2017). A deeper insight into the soil thermal structure is required to understand changes in its thermal properties with depth, relevant for land surface models (LSMs) and therefore with implications for land–atmosphere interactions affecting Earth system (Smerdon and Stieglitz, 2006; Ekici et al., 2014; Hermoso de Mendoza et al., 2020; Steinert et al., 2021) and forecast (Miralles et al., 2019) models. A better characterization of subsurface thermal properties can also contribute to more accurate land heat uptake estimates, i.e., energy storage in the subsurface due to industrial warming, which is the second component contributing the most to terrestrial energy partitioning (Beltrami et al., 2002; Cuesta-Valero et al., 2023; von Schuckmann et al., 2023).

This paper assesses the validity of the SAT–GST coupling and conductive thermal propagation hypotheses in the penetration of SAT perturbations into the subsurface in a high-altitude area in central Spain. This is achieved by determining the apparent soil thermal diffusivity from four short (4–6-year) and two longer (ca. 20-year) subsurface temperature records obtained at six sites in the area of the Sierra de Guadarrama and studying temperature variability in time and depth through the soil and down to the first meters of bedrock. The dataset presented in this work is unique since it gathers shallow and deep subsurface information (down to 20 m) from different sites and sensors in a high-mountain area. The complex orography of the Sierra de Guadarrama was shown to dominate the spatial mean and variability of SAT, while temporal variability was found to be very homogeneous over the region (Vegas-Cañas et al., 2020). In a first stage, we explored the SAT–GST coupling at different sites, considering the influence of processes at seasonal timescales. The propagation of the annual wave with depth was subsequently studied. This allowed us to test the conductive hypothesis and to establish a comparison between different sites by using the standard approach of focusing on the propagation of the annual wave with depth. The analysis was further extended by assessing changes in thermal diffusivity in the subsurface temperature profile, which were linked to the soil–bedrock mineral composition transition observed from samples extracted at each of the sites.

The study of the conductive regime also explored changes in apparent thermal diffusivity with time. This was feasible for the longer records for which the analysis can be performed, focusing on different multi-annual temporal intervals. This approach was however hampered by the short time span of observations at most sites considered herein. This problem was overcome by introducing a new method that makes use of the information in the spectral domain. Instead of considering only the annual cycle, the propagation of the complete range of harmonics of different frequencies that contribute to temperature variability was used. This allowed for focusing on intra-annual timescales and using shorter time series to estimate soil thermal diffusivity near the surface. The fact that apparent thermal diffusivity values can be retrieved from shorter time spans also allowed diffusivity changes to be evaluated through time and connected to changes in soil moisture content (Sorour et al., 1990; Fuchs et al., 2021).

Table 1 Name, code, coordinates (longitude, latitude), altitude, and time span of available subsurface temperature data at both the GuMNet and AEMET observational sites.

2 Data

The subsurface data analyzed hereinafter consist of temperature time series from six locations distributed over the Sierra de Guadarrama. Table 1 includes a more detailed description of the name, code, geographical position, and date range of available data of each location. The Sierra de Guadarrama is part of the Sistema Central, a mountain system that splits the Spanish Central Plateau into a southeastern side, with altitudes of around 600 m a.s.l. (above sea level), and a northwestern side, of higher elevation (ca. 750 m a.s.l.). Elevations in the sierra span from 900 at the foothills to ca. 2200–2400 m a.s.l. at the summits. Two sites, Puerto de Navacerrada and Segovia, belong to the Spanish Meteorological Service (Agencia Estatal de Meteorología, AEMET), and the rest are part of the Guadarrama Monitoring Network (GuMNet). The latter was created in 2014 with the aim of gaining further insights into climate variability in mountain environments in central Spain (Vegas-Cañas et al., 2020). GuMNet consists of 10 atmosphere and subsurface monitoring stations, which cover a vertical gradient ranging from 900 to 2200 m a.s.l., spreading over the mountainous terrain of the Sierra de Guadarrama. In this work, four stations were selected that provide information on atmospheric variables and subsurface temperatures (Fig. 1a). In all the cases, vegetation cover at the sites consists of short grass that changes minimally during the year. No site has trees or shrubs within the fenced perimeter of the monitoring station.

Figure 1 (a) Subsurface data availability over the area of the Sierra de Guadarrama used in this work. The map in the inset (upper-left corner) shows the location of the Sierra de Guadarrama in the Iberian Peninsula. The extension covered by the map in panel (a) is represented by the orange-shaded area in the inset. Right: conceptual sketches of subsurface equipment for 1 m trenches (b), 2 m (c), and 20 m boreholes (d). The same distribution of subsurface temperature (referred to as ST in the figure) sensor depths is carried out at every site (c, d). Sensors in trenches are unevenly placed within the first meter of the ground at different sites. Pale yellow (grey) background color refers to soil (bedrock) in all the panels. The transition depth between soil and bedrock (blurred area in panel d) is only intended for illustration and takes place at different depths for each site. Colors and shapes of the symbols are also illustrative. In addition to subsurface temperature data, all the stations provide SAT measurements as well as soil moisture content (referred to as SH in the figure) in the trenches.

The subsurface thermal regime is monitored by two different arrangements. Near-surface soil temperatures are measured at various levels in trenches (TRCH; Fig. 1b) at all the sites. Trenches are excavations forming a slope in which the frontal wall was used to analyze the vertical structure of soil horizons (Brady and Weil, 2017) and to insert soil temperature reflectometer sensors down to a depth of about 1 m. Trenches were refilled with the previously extracted material. Deeper temperatures in the soil and bedrock are monitored within cased boreholes. They consist of a small diameter (ca. 76 mm) cylindrical drill in depth where a casing was inserted and afterward filled up with a silicone gel. Temperatures therein are measured with calibrated platinum resistance thermometers embedded inside the casing (BRH; Fig. 1c, d). For the sake of having a finer vertical resolution in the soil near the surface and some level of redundancy, two boreholes were installed at these sites: a deep one going down to 20 m (BRH 20 m in Fig. 1 and Table 2) and a shallow one of 2 m (BRH 2 m in Fig. 1). Table 2 provides a detailed list of the depths at which subsurface temperature sensors were installed. Levels with malfunctioning sensors are marked with a dash.

Table 2 Depths corresponding to the available individual temperature series in the trenches and boreholes at each of the sites used in this work (see Table 1 for the codes). Dashes indicate missing data due to sensor malfunction.

Figure 2 Subsurface mineral composition of the borehole cores extracted when drilling at CTS, HYS, RPI, and HRR. The top soil layer (red) covers approximately the first soil meter at every site, whilst sediment (brown) depth varies from one site to the other. Bedrock (magenta) underneath is gneiss or granite, normally weathered at the top. In the case of CTS, some photos of the material at different depths are also included for illustration.

During the borehole drilling process at the GuMNet sites, the extracted cores were preserved. They were subsequently sampled and characterized. The result is shown in Fig. 2, where the limits between the soil, sediment, and rock layers were determined. The superficial soil layer is made up of different soil horizons (Brady and Weil, 2017) containing an O and/or A horizon (i.e., topsoil) followed by E and/or B (i.e., subsoil) and several C (i.e., parent material) horizons that reach a depth of approximately 1 m. The underlying sediment layer contains eroded material of loose debris rocks that have accumulated at the corresponding locations. The Cotos (CTS), Hoyas (HYS), and Raso del Pino I (RPI) profiles are located over gneiss rock and Herrería (HRR) on granite rock. In all the profiles, the upper layers have weathered bedrock, and with depth this bedrock is unweathered and compact. This information was used in Sect. 4 to compare changes in apparent thermal diffusivity to changes in mineral composition and textures with depth.

Figure 3 Example of errors in BRH 2 m at CTS in the year 2017. The 10 min GST observations (grey) and those of temperatures for seven levels below the surface are shown. Time series are shifted an offset equivalent to 10 K and ticks are included every 5 K for clearer visualization. The value of 270 K is shown for all levels as a reference. The olive-colored band at the bottom represents the level of 230 K for the 2 m depth temperature time series used to determine nonphysical values. The segment of the level of 1.5 m shown by the arrow indicates a shift in the mean temperature value. Time values after the shift are subsequently eliminated.

Subsurface temperature measurements are taken at a 10 min time resolution at every GuMNet station. As a first step prior to the analysis carried out in this work, subsurface data were subjected to quality-control and resampling procedures. AEMET stations Puerto de Navacerrada and Segovia (NVC and SGV) were independently processed and resampled to daily resolution as described in Melo-Aguilar et al. (2022). The quality control of GuMNet sites focused first on removing data outliers. Since subsurface temperature sensors are connected in series in both boreholes and trenches (see Fig. 1b–d), most often erroneous extreme values are recorded at all levels simultaneously. Thus, to avoid discarding values that might be correct meteorological extremes, this correction was applied for both boreholes and trenches using the deepest time series to detect the erroneous values. The lowermost level was used because it shows the lowest temperature variability due to conductive damping of the surface signal. This reduces the range of high-frequency variability in the series except for the erroneous data, which makes non-reliable outliers easier to detect (Fig. 3). First, outliers flagged as physically implausible were removed, i.e., values falling below 230 K or exceeding 350 K. The detection of outliers below 230 K can be considered redundant since the next step would have detected all the values lying beneath this threshold. Second, the remaining spikes were flagged and discarded. Spikes were detected by calculating the sequential differences of each value minus its preceding one and screening the values in the interval $(Q_{\mathrm{1}}-\mathrm{3}IQR,Q_{\mathrm{3}}+\mathrm{3}IQR)$, where Q₁, Q₃, and IQR are the first and third quartiles and the interquartile range of the distribution of temperature differences (Tukey, 1977). Time steps whose difference values were missing or outside of the screening interval were removed. Figure 3 shows an example of such spike screening for the BRH 2 m at CTS. The outliers depict coinciding failures at all the levels. The procedure is very conservative as the bottom levels have very low temperature variability and all spikes are identified and removed at each level. Note the decrease in high-frequency variability with subsurface depth, which is consistent with a heat conduction process (Carslaw and Jaeger, 1959). Once outliers were removed, time series were subjected to a one-by-one inspection, and periods with clear drifts or strong changes in variability were manually removed (check Table 3 for a list of removed periods at CTS and HYS). Figure 3 also shows an example of one such period for the 1.5 m depth temperatures at CTS (green arrow).

Table 3 Record of periods discarded at CTS and HYS for the analysis subsequently performed in this work. Records are sorted by log depths (descending) and span from the dates indicated up to the present.

After applying the aforementioned quality-control procedures, the GuMNet series were resampled at hourly resolution, which filtered out noisy intra-hourly variability unnecessary for the subsequent analysis of the data. This was attained by calculating pseudo-hourly time series, which was carried out with a nearest-neighbor resampling technique (Brandsma and Können, 2006) by selecting the 10 min data closest to the hour. Finally, the resulting series were averaged to obtain daily mean subsurface temperatures that were subsequently used to analyze the propagation of the annual cycle.

In addition to subsurface temperature data, all the stations include SAT measurements as well as soil moisture content within the trench-monitoring levels. SAT is recorded at 2 m above the ground surface at all the sites using a temperature probe encapsulated within a shield. Snow cover data are also available at CTS, HYS, and RPI (see the codes in Table 1). Snow depth is measured with an ultrasonic snow depth sensor that captures the distance to the surface. Depth values are obtained from the subtraction of the surface level measured from the reference ground surface level value. These additional data were used to explore potential drivers of the SAT–GST seasonal decoupling and temporal changes in soil apparent thermal diffusivity at some sites.

3 Methodology

The SAT–GST seasonal coupling was explored by quantifying the offset and correlation of seasonal (December–January–February, DJF; June–July–August, JJA) and annual temperatures at 2 m above the ground surface with GSTs taken from the top-level temperatures in the BRH 2 m setup. The estimates were obtained from the complete time span of available data at each site (Tables 1 and 2).

Furthermore, the subsurface temperature distribution was analyzed by assuming that GST perturbations are propagated downwards by heat conduction. Considering that every level is at a thermal equilibrium state, horizontal heat transport may be neglected. Thus, the problem is reduced to resolve the one-dimensional time-dependent heat conduction equation, $\partial T/\partial t=\mathit{\alpha }{\partial }^{\mathrm{2}}T/\partial z^{\mathrm{2}}$. If GST is assumed to follow a sinusoidal (e.g., annual) cycle, subsurface temperatures at any depth z can be derived as follows (Hurley and Wiltshire, 1993; Smerdon and Stieglitz, 2006):

$$\begin{array}{}\text{(1)}& T(z,t)=T_{\mathrm{0}}+A_{\mathrm{0}}{e}^{-z\sqrt{\mathit{\pi }{f}_{\text{a}}/\mathit{\alpha }}}\cos (\mathrm{2}\mathit{\pi }{f}_{\text{a}}t-z\sqrt{\mathit{\pi }{f}_{\text{a}}/\mathit{\alpha }}),\end{array}$$

where T₀ and A₀ are the mean temperature and wave amplitude of the annual cycle at the ground surface, respectively, f_a is the frequency of the annual cycle, and α is the thermal diffusivity, which is a direct function of thermal conductivity (k) and the inverse of density (ρ) and heat capacity (c), i.e., $\mathit{\alpha }=k/\mathit{\rho }c$. Equation (1) shows that the annual cycle amplitude is attenuated and phase-shifted with depth. The amplitude attenuation is exponential, and the phase shift is linear with depth and dependent on the thermal diffusivity; i.e., increasing this parameter produces both a greater amplitude attenuation and a phase shift.

Amplitude and phase values of the annual cycle at every level were obtained by least-squares fitting subsurface temperature series at a daily resolution to an annual-period sinusoidal wave. Then, linear regression analysis of the annual-cycle amplitude (phase-shift) values with depth yielded estimates of the apparent thermal diffusivity of the subsurface, since the resulting slope is equal to $\sqrt{\mathit{\pi }{f}_{\text{a}}/\mathit{\alpha }}$ in Eq. (1). Amplitude and phase-shift values coming from every installation (TRCH 1 m, BRH 2 m, and BRH 20 m) were independently adjusted and normalized to prevent any disruption due to blending data from different sources when assessing heat propagation from the ground surface to 20 m depth. Normalization was achieved by dividing the amplitude values of the annual cycle at every level by the amplitude value at the ground surface and by subtracting the phase-shift value at the ground surface from the phase-shift values of the annual cycle at every level. The values at the ground surface were obtained through a linear regression adjustment. Once normalized, all amplitude and phase-shift values were brought together and linearly adjusted to obtain a single estimate of thermal diffusivity for the whole subsurface profile at every station. A similar methodological approach based on the analysis of the annual cycle has been widely used in the literature (e.g., Smerdon et al., 2004; Pollack et al., 2005; González-Rouco et al., 2009), and it will be referred to as the classic analytical approach (CA) hereafter. This CA framework was also used to estimate changes in apparent thermal diffusivity with depth, either by considering changes between pairs of subsurface levels or by identifying depths where significant changes in apparent thermal diffusivity occur. This was achieved using a two-phase regression analysis (Solow, 1987, 1995; Melo-Aguilar et al., 2018), which determines whether there is a change point for which linear fits to the segments before and after the change point significantly improve the results compared to a linear fit to the complete curve.

Beyond the regular nature of the annual cycle, GST variability encloses perturbations at different frequencies. If GST is considered a sum of temperature harmonics, each harmonic is also propagated with depth following Carslaw and Jaeger (1959):

$$\begin{array}{}\text{(2)}& T(z,t)=T_{\mathrm{0}}+\sum _{i=\mathrm{1}}^N{A}_{\mathrm{0}}{e}^{-z\sqrt{\mathit{\pi }{f}_i/\mathit{\alpha }}}\cos (\mathrm{2}\mathit{\pi }{f}_i-z\sqrt{\mathit{\pi }{f}_i/\mathit{\alpha }}),\end{array}$$

where f_i is the frequency of each wave component in a spectral decomposition, as yielded by a Fourier periodogram for instance (Brigham and Morrow, 1967). Equation (2) shows that amplitude attenuation and phase shift are frequency-dependent. More precisely, amplitude attenuation grows exponentially with frequency; i.e., high-frequency harmonics decay more intensively with depth, whilst low-frequency perturbations can propagate more deeply. For instance, the daily cycle hardly reaches a depth of 50 cm, while the annual cycle propagates down to 15 m (Putnam and Chapman, 1996). The temperature spectral attenuation with depth can be expressed as

$$\begin{array}{}\text{(3)}& \mathit{\zeta }=\sqrt{P_z\left(\mathit{\omega }\right)/P_{\text{GST}}\left(\mathit{\omega }\right)}=e^{-z\sqrt{\mathit{\pi }f/\mathit{\alpha }}},\end{array}$$

where P_GST(ω) and P_z(ω) are the spectral power densities of GST and that at a given level z, respectively:

$$\begin{array}{}\text{(4)}& P_z\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{N_{\text{s}}}}{\left|\sum _{t=\mathrm{1}}^{N}T(z,t)e^{-i\mathit{\omega }t}\right|}^{\mathrm{2}},\text{with}\mathit{\omega }=\mathrm{2}\mathit{\pi }f,\end{array}$$

where N_s is the temporal resolution of the temperature time signal at z, T(z,t); f is the frequency; and ω is the angular frequency of the different wave components of T(z,t). In order to minimize the undesirable effect of Gibbs' noise at high frequencies, P was calculated by making use of Welch's modified periodogram (Stoica and Moses, 1997). Because the periodogram calculation requires evenly spaced data, data gaps were previously filled by linear interpolation. As the spectral attenuation must follow an exponential relationship with frequency, the spectral amplitude attenuation curves attained from the analysis of Eq. (3) were subsequently linearized (i.e., taking the natural logarithm) and linearly fitted by least squares to the values obtained at all available depths:

$$\begin{array}{}\text{(5)}& (\ln \mathit{\zeta }{)}^{\mathrm{2}}={\displaystyle \frac{\mathit{\pi }{z}^{\mathrm{2}}}{\mathit{\alpha }}}f.\end{array}$$

Estimates for the apparent thermal diffusivity of the soil layer in between the ground surface and level z can be retrieved by dividing πz² by the regression line slope, $\mathit{\pi }{z}^{\mathrm{2}}/\mathit{\alpha }$, in Eq. (5). To prevent poor signal-to-noise ratios at high frequencies from biasing the analysis, amplitude attenuation values at frequencies higher than the corresponding e²-fold decay cut-off frequency (i.e., spectral harmonics attenuated more than e² times) were filtered out at every level. To determine this frequency at every site and level, near-surface soil thermal diffusivity values coming from the CA are used. However, it is possible to use a physically plausible prior thermal diffusivity value (ca. 10⁻⁶ m² s⁻¹) to calculate this frequency, leading to very similar results and thus making the spectral method (SpM) independent of the CA. This alternative approach to the CA will from now on be called the SpM. While the CA needs at least a few years of data to provide reliable fits of the annual cycle to subsequently retrieve the apparent thermal diffusivity, the SpM can be applied to data at shorter time intervals. This makes it suitable for analyzing short time series, which is the case with some of the time series presented in this work. Consequently, the SpM was used in this work to retrieve apparent thermal diffusivity values in the first meter of the soil at every station. These values were then compared to the estimates coming from the CA to assess the robustness of the results from both methods. Furthermore, the SpM allows some insight into short-term changes in soil apparent thermal diffusivity and its potential relation to changes in snow cover and soil moisture content.

Figure 4 SAT, GST, and subsurface temperatures (K) at HRR. Subsurface temperatures stem from the BRH 20 m (black lines), BRH 2 m (brown lines), and TRCH 1 m (orange lines). Blue dashed lines show the resulting fit to an annual cycle at each depth.

Figure 5 SAT (triangles) and subsurface temperature profiles (continuous line with diamonds for BRH 2 m, squares for BRH 20 m, and circles for TRCH 1 m) at CTS (blue), HYS (green), RPI (purple), and HRR (red) for their respective whole time spans (see Table 1). Annual mean temperatures are shown on the left-hand side of the figure, whilst JJA (DJF) is shown on the right in dark (light) colors.

4 Results

4.1 Temperature variability at the ground surface

In order to illustrate the amplitude attenuation and shift with depth, Fig. 4 compiles all the subsurface temperature series at hourly resolution coming from BRH 20 m, BRH 2 m, and TRCH 1 m at HRR; SAT is also represented to illustrate the coupling at the surface. A first inspection reveals that the only regular temperature variation at all levels, in both timing and magnitude, is the annual cycle. The daily cycle, while regular in occurrence, changes a lot in magnitude due to weather conditions and seasonality. The annual cycle is most prominent at the ground surface (ca. 9 K amplitude) despite the fact that both SAT and GST also show remarkable daily temperature variations (ca. 4 K amplitude). It is noteworthy how the annual cycle decreases in amplitude with depth to be imperceptible at 15 m, where annual temperature variations are on the same order as the accuracy of the sensor, i.e., 0.1 K. Furthermore, the increasing lag in the occurrence of temperature maxima in summer and minima in winter are indicative of the phase shift due to conductive propagation of the temperature variations. Both phenomena can be inferred by observing either the subsurface temperature series or the sinusoidal curves adjusted to an annual cycle (blue dashed lines). Moreover, SAT and GST daily variabilities are similar during the year except in winter, when GST variability is smaller in comparison with SAT. For other soil levels, warm or cold anomalies that last a few weeks are propagated only in the first meter of the soil.

The long-term mean temperatures with depth for CTS, HYS, RPI, and HRR are shown in Fig. 5. The annual vertical mean temperature of the profile, which is observed below the depth at which the annual cycle has been attenuated (ca. 15 m), is remarkably higher at HRR (286.8 K) than at CTS (281.3 K), HYS (279.2 K), and RPI (280.9 K). The JJA and DJF temperature profiles further depict the seasonal temperature variations towards warmer values in summer and colder values in winter as well as the amplitude attenuation of these seasonal variations with depth. JJA and DJF temperature profiles converge at a depth of around 5 m at all the sites to then cross over, which is a sign of the phase shift of the annual cycle. The transition from the BRH 2 m profile to the BRH 20 m profile shows noticeable discontinuities that are on the order of 2.7 and 2.2 K at HYS and about 1.5 and 1 K at HRR for the 1.5 and 2 m temperature levels, respectively. These discontinuities are most likely due to heterogeneous soil and sediment formations. The uppermost soil and sediment layers vary according to the sites and in all cases reach a depth greater than 2 m (Fig. 2). In the case of HYS, for instance, the greatest annual temperature difference is observed between the 2 and 20 m boreholes at 1.5 m depth. The site is on a middle slope with a slope gradient of 10 % to 15 %, and the soil horizons have a sandy loam to sandy clay loam texture. What is most notable is the varying content of rock fragments (15 % to >80 %) within the different soil horizons. This is further observed throughout the sediment layer, where the rock fragment distribution varies with depth and is very weathered. Therefore, temperature variations over short distances and within the uppermost soil and sediment layers are expected due to their different properties, which will also influence the water infiltration and circulation. With regard to SAT–GST coupling, Fig. 5 shows that the annual cycle is less pronounced for SAT than GST due to radiative warming and cooling of the ground surface, resulting in a negative SAT–GST offset in summer and a positive one in winter. The distribution of SAT and GST mean values of the sites agrees with their different altitudes, which yields a vertical gradient of −5.81 K km⁻¹ in SAT given by Vegas-Cañas et al. (2020) for the area of the Sierra de Guadarrama. SATs in Fig. 5 comply well with this value, rendering a gradient of −6.55 K km⁻¹, which is larger in magnitude due to this subset of sites being distributed over the higher altitudes of the region. Subsurface temperatures with depth follow the distribution of SAT at HYS and HRR, with HRR exhibiting higher values for the whole profile. However, this pattern is not observed at RPI and CTS, as RPI shows lower temperatures below 5 m despite SAT and GST being warmer, which might be due to processes within the soil that have not been further investigated in this paper.

Figure 6 SAT, GST, and subsurface temperature anomalies (K) with respect to the annual mean at CTS (a) and HRR (b). (c) SAT (light) and GST (^∘C; dark colors) changes during 10 d in January (blue) and July (red) 2019 and the evolution of the snow cover (m; light pink) at CTS and HRR.

Temperature anomalies with respect to the long-term mean at each depth are shown in Fig. 6 as an example for CTS and HRR. This allows for a clearer visualization of the amplitude attenuation and phase shift with depth. Two cycles are clearly identified again at both sites: the daily cycle, which is visible at GST and hardly noticeable at 0.5 m (Fig. 6a, b), and the annual cycle, with an amplitude reduction and phase shift that are traceable down to 15 m. Nevertheless, there are two differences that are worth mentioning. On the one hand, GST and subsurface temperatures underneath at CTS show some constant signal time periods in winter that are not found in SAT, as shown in Fig. 6c, where SAT and GST are represented for a shared period in January 2019 at CTS and HRR. As soon as a snow layer emerges at CTS on 19 January 2019, GST oscillations are suppressed despite SAT still varying. This is due to the combination of the soil-insulating effect of snow cover (Bartlett et al., 2005; Zhang, 2005) and the zero-curtain effect, which prevents temperatures at the ground surface from dropping far below zero due to latent heat release (Outcalt et al., 1990). This behavior is not shared by HRR, since at this lower altitude a regular seasonal snow layer is less frequent. Hence, HRR is snow-free during the same dates, so SAT and GST variations are strongly coupled. On the other hand, note that SAT variability differs from one site to the other, with the amplitude of the daily cycle greater and the annual cycle smaller at CTS than at HRR, as illustrated by greater differences in SAT variations in July 2019 with respect to those in January 2019. This is consistent with SAT annual variability decreasing and intra-annual increasing with altitude in the same area, as shown by Vegas-Cañas et al. (2020).

Table 4 SAT–GST differences in the mean (K) and Pearson's correlation coefficient using daily data. The annual wave was previously filtered out from both SAT and GST data for this analysis. Non-significant (p<0.05) differences and correlations are shown in bold.

Air–ground temperature seasonal decoupling is further explored in Table 4, where SAT–GST differences and Pearson's correlation coefficient were computed after filtering out the annual cycle. SAT–GST differences show that, on average, GST is warmer than SAT, with annual differences being negative. Correlations using daily data are all significant (p<0.05, Student's t test) and above 0.58 for all the sites. Differences tend to be larger during DJF, likely due to the insulating effect of snow cover, which hampers soil cooling and disrupts SAT–GST coupling. During JJA, significantly warmer GSTs are found at HRR, NVC, and SGV, likely due to surface radiative warming. Additionally, correlation analysis yields the weakest values in DJF at the highest-altitude sites, i.e., RPI, HYS, and CTS, which may be linked to the presence of a steadier seasonal snow layer than at the lowest sites. Low correlation is found in JJA at HRR, which might be due to weaker latent heat fluxes connected to Iberian soils' summer drought (Vicente-Serrano et al., 2013). Nevertheless, the diverse local responses at the various sites could be due to land cover differences, water or latent heat exchanges, or other factors (Cermak and Bodri, 2016; Cermak et al., 2017) that have not been further examined herein.

Figure 7 Phase shift in days (blue) and the logarithm of amplitude ratios (red) vs. depth at CTS, HYS, RPI, HRR, NVC, and SGV (see Table 1). Points in orange, brown, and black stem from a least-squares adjustment for sinusoidal fitting of the annual cycle of TRCH 1 m, BRH 2 m, and BRH 20 m data, respectively. Diffusivity values (10⁻⁶ m² s⁻¹) are retrieved by linearly fitting all the points available at every site (see the text for details). Regression lines, confidence intervals using a significance level of p=0.05, and diffusivity values for the logarithm of the amplitude ratios (phase shift) are included in red (blue).

Table 5 Apparent thermal diffusivity values (10⁻⁶ m² s⁻¹) obtained from the amplitude attenuation, ln($A/A_{\mathrm{0}}$), and phase shift, ϕ, for TRCH 1 m only (first column), the first meter of depth using a blend of TRCH 1 m and BRH 2 m data (second column), and BRH 20 m (third column) at CTS, HYS, RPI, HRR, SGV, and NVC.

4.2 Subsurface thermal regime: the CA

The subsurface thermal structure in the area of the Sierra de Guadarrama is assessed hereafter. The conductive propagation of GST changes is first addressed by applying a CA (Sect. 3). Results for a CA over the whole profile for every site are included in Fig. 7. As can be seen, CTS, HYS, RPI, and HRR integrate data coming from both BRH 2 m, BRH 20 m, and TRCH 1 m down to 20 m, while NVC and SGV provide TRCH 1 m data. Hence, apparent thermal diffusivity values only account for the first meter underground at these two sites. For the rest of the sites, all the information at each location down to 20 m is used to derive apparent thermal diffusivity estimates. Table 5 segregates results for the case in which only TRCH 1 m would be considered as well as for the cases that include information available down to 1 m (TRCH 1 m + BRH 2 m) and down to 20 m using only BRH 20 m. This information complements the estimates provided in Fig. 7.

Figure 8 Two-phase regression of the points shown in Fig. 7. Lines represent the linear trends that are the best fit to the data before and after the estimated point of change for the logarithm of the amplitude ratio (phase shift), depicted in red (blue). Apparent thermal diffusivity values (10⁻⁶ m² s⁻¹) before and after the change point are given as well in every case. The 20 m level at CTS was omitted from the analysis.

Thermal diffusivity estimates, i.e., apparent thermal diffusivities, yielded from the slope of the regression lines for both phase shift and amplitude attenuation are very similar for all the deep sites, generally in the range of 1–$\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹. These values are close to the mineral diffusivity values for quartz (ca. $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹) and feldspar (ca. $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹), two of the main components of gneiss and granite (de Vries, 1963), which are the dominating materials found at the sites of this study (Fig. 2). The observations at SGV render values that fall below this range, with apparent diffusivity values between 0.5 and $\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹. While subsurface temperatures at CTS, HYS, RPI, and HRR extensively cover the soil and bedrock depths, NVC and SGV only include data within the first meter of the soil, which can yield very different thermal diffusivity estimates with respect to the other sites. Near the surface, the different composition, compactness, and soil moisture content can account for the different apparent diffusivity values. In fact, if apparent diffusivity is estimated at all the sites using only TRCH 1 m observations, values in the range of 0.46 (RPI) to $\mathrm{0.88}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ (CTS) are attained (Table 5). Furthermore, if both TRCH 1 m and BRH 2 m temperatures within the first meter are considered, values hardly change and range from 0.42 (RPI) to $\mathrm{0.87}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ (CTS). It is often the case that apparent diffusivity values rendered by amplitude attenuation and phase shift within the first meter of the soil can differ considerably, as shown for NVC in Fig. 7. This is an indication that changes in apparent thermal diffusivity occur towards the surface, as seen at NVC and SGV in Fig. 7. Additionally, changes in the conductive regime may occur within these depths due to soil moisture variability (Melo-Aguilar et al., 2022), freezing and thawing processes (Gao et al., 2008), and near-surface evaporation and latent heat exchanges (Tong et al., 2017) affecting the soil.

Figure 9 Thermal diffusivity (α, 10⁻⁶ m² s⁻¹) and wave velocity (ν, m yr⁻¹) profiles with depth (m). The red curves were obtained from normalized amplitude attenuation values ln($A/A_{\mathrm{0}}$) and the blue ones from phase-shift values with depth. Dashed (solid) lines represent the raw (three-running-layer thickness-weighted average) estimate. Points at the bottom of the profile (shaded yellow) represent the thickness-weighted average of the whole subsurface profile. Shaded areas in red (blue) represent the thermal diffusivity values obtained by linearly fitting the amplitude attenuation (phase shift) in Fig. 7.

The regression spread when considering deep sites in Fig. 7 can be remarkable, scaling up to 25 % at RPI. Even if the logarithm of the amplitude damping and phase shift for the whole soil column can be considered linear in a first approach, values do not seem to distribute randomly around the regression line. For instance, between 2 and 7 m all the sites tend to show larger-amplitude attenuations and phase shifts, hence indicating changes in the regression line slope, i.e., diffusivity, with depth. A further analysis using a two-phase regression model (Fig. 8; Solow, 1987, 1995) confirms there is a tilt in the slope, indicating a significant change in the apparent thermal diffusivity values with depth. Even though the coefficient of determination is very high and significant for every site for both the simple CA (Fig. 7) and two-phase CA (Fig. 8) regression lines, always above 0.9, there is a significant reduction in the square sum of residuals in the second case. For instance, this quantity plunges from 5.67 to 1.02 at HYS and from 8.22 to 0.62 at RPI for amplitude attenuation and from 3.21 to 0.91 and from 11.46 to 1.20 d² for phase shift, respectively, thus indicating that the two-phase regression represents the data more accurately than a single-slope linear regression model. Since amplitude attenuation and phase-shift values at every depth were normalized before applying the CA (see Sect. 3), the change point position is certainly not produced by gathering information from TRCH 1 m and the BRHs. A change in apparent thermal diffusivity might be indicative of the soil–bedrock transition, where porous non-consolidated material is progressively substituted by consolidated and more compacted material, which yields higher diffusivity values with depth (de Vries, 1963; Pollack et al., 2005). Apparent diffusivity values scale up from 0.7–1 to 1.4–$\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ at most of the sites, which is consistent with the aforementioned transition. In Fig. 8, CTS and HYS show similar change point values for the amplitude attenuation and the phase shift analysis around 5 m. However, there is more uncertainty at RPI and HRR, where amplitude attenuation yields shallower change points (7 and 4 m) than the phase shift analysis (11 and 8 m, respectively).

Figure 10 SpM amplitude attenuation analysis at NVC (a) and SGV (b). Monthly time series at both sites were selected for a common interval spanning from 1999 to 2008. Missing data were linearly interpolated. Spectral attenuation curves are calculated as a ratio of amplitudes of Welch's periodogram at the (Stoica and Moses, 1997) 0.1 (orange), 0.2 (red), 0.5 (green), and 1 (blue) m temperature levels with respect to GST. Solid (dashed) lines represent observed (theoretical) curves at every level. Diffusivity values below come from least-squares fitting the observed to theoretical curves. Diffusivity values in grey (top left of each panel) show the result of the CA for the amplitude attenuation of the annual cycle ($\mathrm{1}/\mathrm{12}$ month⁻¹) shown in Fig. 7. Colored points at the annual frequency (vertical grey dashed line) show the corresponding result of the CA for the amplitude attenuation at every level with respect to GST, which can be compared with the corresponding SpM observed (solid) and theoretical (dashed) values.

The previous changes in apparent diffusivity inferred by the CA can be partially traced to the soil–bedrock transitions observed in Fig. 2 (see Sect. 2), where the profiles of subsurface composition have a superficial soil layer which is succeeded by accumulated sediments and crystalline materials of granite and gneiss. Results coming from the CA and the stratigraphy of the samples agree at HYS, placing the soil–bedrock transition at a depth of 5 m. However, if the change points should be due to the bedrock limit, the CA underestimates the observed depth at CTS, where the bedrock starts in the sampled cores at 8 m (Fig. 8). The fact that the amplitude attenuation and phase shift CA results do not coincide at RPI and HRR might also be an indication that other factors influence the occurrence of change points (e.g., soil moisture), the CA failing again to produce consistent amplitude and phase estimates for change points, both above the values of sampled cores that show the soil–bedrock transition at 5 and 2 m, respectively. Considering the combined results from the stratigraphy and the CA for the four sites, it becomes clear that there is some uncertainty in the relationship between subsurface material and apparent thermal diffusivity changes, and other influences play a role in producing the step changes in thermal properties with depth.

The aforementioned sample analysis illustrates that changes in subsurface material do not occur abruptly but rather progressively into depth, contributing to changes in apparent thermal diffusivity with depth as well. For that reason, subsurface diffusivity profiles were derived by applying the CA level by level to CTS, HYS, RPI, and HRR (Fig. 9). This shows that diffusivity increases with depth at every site, which agrees with a previous analysis conducted at Fargo (Smerdon et al., 2003). Whilst RPI and HRR manifest a more linear behavior, apparent thermal diffusivity values at CTS and HYS differ considerably from one layer to the other, with the profiles being particularly variable near the surface. In general, the mean diffusivity values and spread attained in Fig. 7 are well represented by the profiles, being higher at the deepest layers. Since diffusivity estimates in Fig. 9 are weighted by the layer thickness, small values in the thinner soil layers near the surface are down-weighted and the total vertical estimates tend to be greater than those provided in Fig. 7. Overall, Fig. 9 agrees well with results derived from Figs. 2, 7, and 8 and Table 5 and illustrates better the deviations from the linear behavior for the natural logarithm of amplitude attenuation and the phase shift with depth identified in Fig. 7. These deviations are explained by an increase in apparent thermal diffusivity with depth.

4.3 Subsurface thermal regime: the SpM

Heterogeneity near the surface and soil moisture spatial and temporal variability can produce changes in diffusivity. To study this effect, the SpM is introduced in the analysis, taking into account the damping of all temperature perturbation amplitudes in the frequency spectrum but overlooking alterations of the conductive regime due to convection, advection, or latent heat exchange processes (Gao et al., 2008; Tong et al., 2017). Nevertheless, it can provide some hints of the diffusivity changes with depth and mostly with time. An example of its performance is given in Fig. 10, where it was applied to monthly temperature data at SGV and NVC in the 1999–2008 period. Apparent thermal diffusivity was retrieved at the 0.05, 0.1, 0.25, and 0.5 m layers assessing the spectral amplitude attenuation of GST in its propagation to the 0.1, 0.2, 0.5, and 1 m levels, respectively. Results for both sites show that soil thermal diffusivity increases with depth, being higher for NVC than for SGV. The value of the apparent thermal diffusivity for the whole soil column derived from the SpM at SGV (see the points in Fig. 10) agrees with the CA estimate at the annual frequency from Fig. 7 ($\mathrm{1}/\mathrm{12}$ month⁻¹; 0.55 for CA vs. $\mathrm{0.56}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ for the SpM), while it is slightly smaller at NVC (0.82 vs. $\mathrm{0.76}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹). The consistency of the CA and SpM in reproducing the amplitude attenuation of the annual cycle indicates the robustness of the SpM in reproducing CA results. Furthermore, the SpM allows for extension of the amplitude attenuation analysis to the whole harmonic spectrum. This means that shorter periods, i.e., higher frequencies than the annual cycle, can be considered. In fact, the application of the SpM is not limited to data at monthly resolution but is applicable to higher temporal resolutions.

Figure 11 expands the use of the SpM as an example to estimate thermal diffusivity values for weekly (a), 3-daily (b), daily (c), and 3-hourly (d) resolutions at HRR. One can observe how the SpM consistently reproduces the amplitude attenuation of the annual cycle stemming from CA at every resolution. Moreover, Fig. 11 shows the effect of frequency in amplitude decay: whilst in Fig. 11a the observed amplitude attenuation curves show an exponential decay in every case, with values over $\mathrm{1}/e^{\mathrm{2}}$ except for 1 m, all the levels in Fig. 11d are affected by poor signal-to-noise ratios at the high frequencies. Therefore, filtering out unreliable amplitude attenuation values at high frequencies with a suitable criterion (e²-fold decay in our case) becomes crucial for preventing biased thermal diffusivity estimates. If the whole frequency band was considered in the SpM, unbiased estimates would still be attained at shallow soil layers and low frequencies, but poor performance should be expected at the deep and high ones, respectively.

Figure 11 As in Fig. 10, SpM analysis of subsurface temperature attenuation at HRR for weekly (a), 3-daily (b), daily (c), and 3-hourly (d) data. The hollow parts of the observed attenuation curves were excluded in the estimation of apparent thermal diffusivity according to the e²-fold decay criterion explained in Sect. 3. (e) Precipitation, soil moisture content, and estimated near-surface thermal diffusivity from applying the SpM at the 0.1 m (light orange) and 0.2 m (light red) levels with respect to the GST with a 10 d running window at 3-hourly resolution; 11 d low-pass running averages are provided in dark colors. Soil moisture content at 4 cm is shown at daily (10-daily) resolution in light blue (blue). Precipitation at 10-daily resolution is shown in purple. The least-square regression correlation (R) for SpM at every time step is given below in each case. The grey-shaded thick line represents the only time interval (4 d) with snow insulation in Fig. 6c for HRR.

As has already been mentioned, soil properties are especially heterogeneous near the surface due to interactive changes in soil material, texture (Cermak et al., 2017), moisture, and hydrology (Pollack et al., 2005; Gao et al., 2008; Tong et al., 2017). The last factor mostly affects subsurface thermal conductivity and volumetric heat capacity by filling up the empty pore spaces embedded in the soil (Arkhangelskaya and Lukyashchenko, 2018), which in turn modifies the apparent thermal diffusivity. Figure 11e explores the relation between soil moisture content and apparent thermal diffusivity near the surface at HRR.

The SpM at 3-hourly resolution (Fig. 11e) in a 10 d running window is applied to derive an intra-monthly evolution of apparent thermal diffusivity at the shallowest layers (0.1 and 0.2 m), which is subsequently compared to changes in soil moisture and precipitation. It is observed that soil moisture presents a seasonal cycle, with extremely low content values in summer and relatively high but variable values in winter and spring. Transitions between these two states are remarkably abrupt, with soil moisture content increasing substantially by the beginning of fall and responding to enhanced precipitation and with the opposite occurring by the end of spring. In contrast to the other sites in the Sierra de Guadarrama, HRR is rarely snow-covered in winter, which usually occurs after high-impact snow events (e.g., Filomena, which hit the Iberian Peninsula in January 2021; shading in Figs. 6c and 11e; Smart, 2021; Tapiador et al., 2021). Apparent thermal diffusivity also shows a seasonal behavior, which corresponds to soil moisture low-frequency variability. In summer, poor soil moisture content drives apparent thermal diffusivity values down. However, even though diffusivity in winter increases with soil moisture content, the relation between these two variables at the high frequencies is not conclusive. For instance, the evolution of diffusivity in winter–spring 2018 at both soil levels correctly mimics soil moisture content variations (Fig. 11e) yet fails to capture its strong variations in winter 2020. This agrees with a clear increase in apparent thermal diffusivity from dry to wet soils but a rather stable diffusivity in wet soils close to its field capacity observed by Arkhangelskaya and Lukyashchenko (2018). An analogous result for thermal conductivity by Dai et al. (2019) complies with this result as well. Moreover, SpM fails to provide reliable estimates of apparent thermal diffusivity when the ground is insulated by the snow, due to signal loss in the whole frequency band. This is clearly shown in Fig. 11e in January 2021, when the R coefficient, i.e., the square root of the goodness of fit of the amplitude attenuation to an exponential function, becomes extremely poor.

Interannual correlation is hard to establish from our assessment due to the shortness of the time series in Fig. 11. However, previous studies (e.g., Melo-Aguilar et al., 2022) suggest that soil moisture and apparent thermal diffusivity changes are proportionally related to some degree. Nevertheless, the lack of strong evidence of a linear relationship between apparent thermal diffusivity and soil moisture near the surface could also suggest that heat conduction is no longer the most important mechanism of heat transfer. Heat convective transport via soil moisture and latent energy fluxes, water advection, and horizontal heat transport due to uneven temperature changes at the ground surface can also contribute to modifying the vertical heat conduction regime near the surface. The use of a running-window SpM over longer time series or developing an alternate SpM based on a conductive–convective amplitude attenuation (Gao et al., 2008; Tong et al., 2017) solution to the heat diffusion equation might better define this relation in wet soils.

5 Conclusions

The subsurface thermal structure at six sites in the Sierra de Guadarrama is consistent with the hypothesis of conduction being the main heat propagation mechanism of SAT changes into depth. On the one hand, SAT changes are transferred to the near surface at the ground. Although SAT and GST are coupled at long timescales, snow insulation and zero-curtain effects produce a decoupling in winter at the uppermost stations (CTS, HYS, and RPI), which are affected by the presence of a seasonal snow cap (Figs. 5 and 6). This result agrees with previous studies carried out at cold and snowy sites, such as Fargo (North Dakota, US; Smerdon et al., 2003), Emigrant Pass (Utah, US; Bartlett et al., 2006), Bistriţa (Romania; Demetrescu et al., 2007), and Hveravellir (Iceland; Petersen, 2022).

On the other hand, GST changes are propagated to depth by heat conduction to an extent controlled by the apparent thermal diffusivity of the subsurface. A CA method based on the amplitude attenuation and phase shift of the annual cycle yielded values that are between 1 and $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ for the whole soil column (i.e., 20 m) at most of the sites (Fig. 7), which are in the range of typical diffusivity values for granite or gneiss, the main materials found in the Sierra de Guadarrama's subsurface (Fig. 2). These values are higher than those observed at the other sites (e.g., Fargo, Prague, or Cape Henlopen; Smerdon et al., 2004), which include shallower measures. If apparent thermal diffusivity is assessed near the surface, the values are smaller than and closer to the previous literature (0.4 to $\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹; Table 5). A subsequent two-phase adjustment (Fig. 8) and a profile derivation based on CA (Fig. 9) portray an increase in thermal diffusivity with depth at all the sites, being more remarkable near the soil–bedrock transition. Evaluating thermal diffusivity changes with depth using a CA can therefore be a useful approach to determine soil depth at observational sites when no borehole rock samples are available. The deep monitoring profiles at the Sierra de Guadarrama are unique and provide a wider insight into subsurface thermal structure, which may have implications for calibrating soil and bedrock properties in models (Tong et al., 2016; Dai et al., 2019) and monitoring land heat uptake in the terrestrial energy budget (Cuesta-Valero et al., 2022, 2023; von Schuckmann et al., 2023).

Near-surface (i.e., 1 m) soil apparent thermal diffusivity vertical and temporal heterogeneity is further explored by introducing a new SpM, which extends the analysis of temperature amplitude attenuation to the whole frequency spectrum. The SpM reports consistent apparent thermal diffusivity values in the first soil meter at SGV, NVC, and HRR with those previously derived with the CA for each of the respective sites regardless of the time resolution of the input data. Then, it was used in an intra-monthly (10-daily) sliding window to explore apparent thermal diffusivity changes in time at the ground surface at HRR (Fig. 11e). Results suggest that this parameter presents a seasonal cycle that is positively correlated with soil moisture, mostly under soil drought conditions. However, the fact that these changes are poorly correlated under high soil moisture contents might indicate that heat conduction is not the main heat transfer mechanism in wet soils. Heat transport via water advection and phase changes might not be negligible near the surface, so the thermal structure in the shallow soil could be better characterized using a modification of the SpM that complies with a conduction–convection heat transfer equation (Gao et al., 2008; Tong et al., 2017). Even so, the results show the potential of the SpM for inferring soil moisture content changes from soil temperatures, becoming a powerful tool for evaluating soil drought and water resource availability from observational and simulation-based subsurface temperature data.