Seasonal Freeze‐Thaw Cycles and Permafrost Degradation on Mt. Zugspitze (German/Austrian Alps) Revealed by Single‐Station Seismic Monitoring

Thawing of mountain permafrost in response to rising temperatures degrades the stability of rock walls and thereby affects infrastructure integrity in Alpine terrain. In this study, we use 15 yr of passive seismic data from a single station deployed near a known permafrost body on Mt. Zugspitze (Germany), to monitor freeze‐thaw processes. The recordings reveal a persistent cultural seismic noise source, which we utilize to compute single‐station cross‐correlations and extract relative seismic velocity changes. We find that parts of the cross‐correlations show seasonal velocity variations (≈3% peak‐to‐peak amplitude) and a long‐term velocity decrease (≈0.1%/yr). Comparison with meteorological data and a previous electrical resistivity tomography study suggests that these velocity changes are caused by freeze‐thaw cycles and by permafrost degradation, respectively. The results demonstrate the potential of passive seismology for permafrost monitoring and suggest that denser instrumentation will provide detailed spatio‐temporal insights on permafrost dynamics in future studies.

the local solar radiation conditions (Hoelzle et al., 2001), and the exposure to the atmosphere. Steep rock walls are usually free of debris, which promotes a rapid response to changes in the thermal forcing (Gruber et al., 2004b), whereas debris and snow cover have an insulating effect. In addition to the thermal forcing through heat conduction, advective heat transfer through (melt)water percolation can rapidly develop deep thaw corridors in fractured rocks (Kane et al., 2001). The complex interplay between numerous processes results in a heterogeneous three-dimensional distribution of mountain permafrost in lenses rather than layers (Krautblatter & Hauck, 2007). Where permafrost is absent, seasonal freezing may be encountered.
Permafrost bodies can be monitored directly through temperature logging in boreholes (Beniston et al., 2018;Haeberli et al., 1998) or through surface-based geophysical imaging techniques including electrical resistivity tomography (ERT) and active seismics (Hauck, 2013;Kneisel et al., 2008). These methods can be used to differentiate between frozen and unfrozen ground (Harris & Cook, 1986;Hauck, 2002;King, 1977;Kneisel et al., 2008;Timur, 1968), or to even infer the temperature distribution in the case of ERT (Krautblatter et al., 2010;Scandroglio, Draebing, et al., 2021). While boreholes allow continuous permafrost monitoring, their wider applicability is limited by the logistics and costs involved. In contrast, active geophysical imaging techniques offer more flexibility but must be applied repeatedly to obtain temporal resolution. This remains challenging, as an automatic acquisition in harsh Alpine terrain is difficult and manual acquisition, for example, on a monthly basis (Mollaret et al., 2019) is laborious.
To circumvent these limitations, passive seismic methods have been recently explored for permafrost monitoring (Albaric et al., 2021;Guillemot et al., 2020;James et al., 2019;Köhler & Weidle, 2019). Applying ambient-noise based seismic interferometry, that is, extracting repeated seismic impulse responses between pairs of seismic stations, James et al. (2019) find seasonal velocity changes of a few percent over the course of the year at a permafrost site in Alaska, which they attribute to active-layer freeze-thaw cycles. Using an array of continuously recording sensors, James et al. (2019) gain both spatial insights on the thaw depth and high temporal resolution such that they can track permafrost thawing caused by water infiltration from snow melt and rainfall events. Similar results are reported by Guillemot et al. (2020), who find seasonal changes in seismic velocities related to permafrost dynamics in the upper 10 m of a rock glacier. In the present study, we investigate the potential of single-station passive seismology for continuous long-term permafrost monitoring. For this purpose, we apply seismic interferometry to a station deployed on the ridge of Mt. Zugspitze (Germany) close to a known permafrost lens. Data are available for the past 15 yr and we extract seismic velocity change time series, which we compare to meteorological records, borehole logs, and a previous ERT study.

Mt. Zugspitze
Located at the German-Austrian border, Mt. Zugspitze (2,962 masl, WGS84 coordinates of the summit: 47.42119, 10.98634) is the highest peak in Germany. The summit area hosts three cable car stations and various other infrastructure including restaurants and broadcasting facilities (Figure 1). Mt. Zugspitze is composed of Triassic limestone (Wettersteinkalk) being weathered and fractured in the summit area and characterized by a subsurface cave drainage system established by karst dissolution (Gude & Barsch, 2005;Krautblatter et al., 2010). Excavations for cable cars and other constructions revealed that the fractures are filled with up to decameter thick ice lenses and frozen loam (Körner & Ulrich, 1965;Ulrich & King, 1993). In prehistoric times, around 3,700 yr before present, a giant 0.3-0.4 3 km E rock slide occurred from the summit region, which is considered to have been triggered by permafrost degradation at the end of the Holocene climatic optimum (Jerz & von Poschinger, 1995).
In present times, permafrost is found in the north-facing rock, whereas the south-facing slopes are almost all free of permafrost (Gude & Barsch, 2005;Nötzli et al., 2010). In 2007, a borehole was drilled beneath the summit, intersecting the crest entirely from south to north on a length of 44 m with a slope of 20 E  toward north. Temperature logging inside the borehole reveals permafrost temperatures down to about -4 E C and seasonal thaw depths of 4.5 and 1.5 m from the southern and northern side, respectively (Gallemann et al., 2017;Nötzli et al., 2010). Another permafrost body extends several tens of meters along the ridge north of the Schneefernerhaus research station (Figure 1b; Krautblatter et al., 2010), which is evident by perennial ice in a gallery intersecting the ridge. While the mean annual air temperature measured at the summit between 1901 and 2000 was −4.8 E C, it was −3.7 E C between 2001 and 2020, hence more than 1 E C higher compared to the 20th century. The increasing temperatures are reflected by permafrost warming and degradation visible in the borehole temperature logs beneath the summit (Gallemann et al., 2017).

Instrumentation
The Schneefernerhaus accommodates the permanent seismic station BW.ZUGS (Department of Earth and Environmental Sciences, Geophysical Observatory, LMU Munich, 2001) in a vault next to the rock face, which is operational since 2006 ( Figure 1). Initially, the station was equipped with a Mark L4-3D short-period seismometer (natural frequency of 1 Hz) and a Lennartz M24 digitizer. In August 2017, this setup was replaced by a Guralp CMG-3T broadband sensor and a Reftek RT130 digitizer. After a larger data gap starting in May 2018, the Guralp sensor was replaced by a Trillium Compact 120 s seismometer, which is operational since July 2019. The ground velocity output of all three sensors is sampled at 200 Hz. Spectrograms for the three different sensors at different times of the year (Figure 1d) reveal a noise source being persistent over years with strong ground vibrations during the day, bound by the operation hours of the cable cars. Despite lower amplitudes at lower frequencies, the noise source is visible down to 2 Hz, where amplitudes are close to the self-noise level of the Mark L4-3D seismometer. Temporary deployments of two additional stations in February 2019 show that the noise amplitudes are stronger for installations closer to the summit area (not shown), suggesting the latter as excitation area. Significant or systematic long-term variations of the noise amplitudes are not evident.

Single-Station Monitoring and Data Processing
We utilize the single station to compute cross-correlations between the different sensor components (E, N, and Z resulting in EN, EZ, and NZ cross-correlations), which can be considered as an impulse response retrieval for a source and a receiver being colocated (De Plaen et al., 2016;Hobiger et al., 2014;Yates et al., 2019). The single-station cross-correlations contain scattered and reflected waves with the EN component being sensitive to Rayleigh and Love waves, whereas the EZ and NZ components relying on vertical ground motions are not sensitive to Love waves. In addition to surface waves, the cross-correlations may also contain body waves reflected at depth including P-to S-and S-to P-converted phases (Becker & Knapmeyer-Endrun, 2019;Hobiger et al., 2014). We use the MSNoise package (Lecocq et al., 2014) to compute daily cross-correlations, which we form by stacking individual cross-correlations calculated from non-overlapping 30-min windows and the frequency range of 0.1-25 Hz. The preprocessing includes the removal of the instrument response from the raw data, clipping of the seismograms at three times the root mean square amplitude, and spectral whitening to equalize the amplitude of all frequencies (Bensen et al., 2007; all MSNoise processing parameters are provided in Table S1). Spectral whitening increases the robustness of the cross-correlations against noise source variability but cannot be applied for auto-correlations of individual channels as this results in a perfect delta pulse not carrying any information on the medium. The inhibited applicability of spectral whitening is the main disadvantage of auto-correlations (Hobiger et al., 2014). This is confirmed by this study, where the auto-correlation results appear to be similar to the single-station cross-correlations results, but noisier. We therefore focus on the single-station cross-correlations in this work. Furthermore, we note that using only the operation hours during day times for the cross-correlation calculation does not change the results, while using night times only degrades their stability.
To extract seismic velocity changes expressed as travel time changes, we compare the time-lapse cross-correlations against a reference cross-correlation. One common approach for this purpose is the moving-window cross-spectral (MWCS) technique (Clarke et al., 2011), where one employs a sliding window along the coda of bandpass filtered cross-correlations. The sliding window is used to determine the travel time shift E t  as a function of the lag time t, averaged over the width of the sliding window and the frequency range of consideration. In a second step, one determines the slope  t t / through a linear regression, which is equal to the velocity change dv v / relative to the reference, if the latter affects the subsurface uniformly. Being related to the MWCS technique, we here explore the recently introduced wavelet cross-spectrum method (Mao et al., 2020), which also yields travel time shifts relative to the reference cross-correlation, but as a function of lag time and frequency E f, that is, ( , ) E t t f  , hence with increased joint time-frequency resolution. This enables us to investigate specific parts of the cross-correlation. Here, we use the wavelet cross-spectrum implementation of the NoisePy package (Jiang & Denolle, 2020), employing a Morlet wavelet to compare 15-day cross-correlation stacks (calculated in a moving window with a step size of 1 day) against the reference. Regarding the reference, we consider two approaches. (a) We calculate the linear reference stack for the fixed reference period of September 1, 2017 to May 1, 2018 using all available daily cross-correlations. (b) As the seasonal freeze-thaw cycle associated with permafrost can cause such strong velocity changes relative to a fixed reference that the measurement of travel time shifts is affected by cycle skipping (James et al., 2017(James et al., , 2019, we also consider a moving-reference type approach to mitigate this problem. , see Text S1 for details).

Results
Active layer thaw and refreeze are governed by the temperature signal propagating into the rock with a period of 1 yr. If cross-correlation coda waves sample freeze-thaw areas, we expect a periodic velocity change signal with the same period. We thus calculate the amplitude of the 365.25-day periodicity,  time series obtained from the wavelet cross-spectrum method relative to the fixed reference in the lag time of −5 to 5 s and frequency range of 1-20 Hz. Because Fourier analysis is hindered by data gaps in the E t  time series, we employ the Lomb-Scargle periodogram (Lomb, 1976;Scargle, 1982), which enables us to compute power spectra at specific frequencies independent of the sample spacing. Figures 2a-2c shows (365.25 ) LS E A d as a function of lag time and frequency, revealing that parts of the cross-correlations exhibit seasonal changes with a period of 1 yr. This is further emphasized by the three E t  time series (converted to  t t / ) representing specific lag-frequency combinations showing clear seasonal variations (Figure 2d). Complementary to the individual  t t / curves, Figure 2d also depicts the lag and frequency dependent  t t / measurements limits (horizontal dashed lines), beyond which cycle skipping occurs. While the two lower considered frequencies stay within these limits, the higher frequency  t t / curve (EZ, 9.92 Hz) is affected by cycle skipping as the variations are larger than the measurement limits.
To systematically investigate temporal changes, we consider the frequency dependence of  t t / . Because we expect localized changes rather than uniform changes, we refrain from determining  t t / via linear regression, and instead determine the median from individual  t t / curves for each frequency bin (for details see Figure S1). We restrict our analysis to the lag range bound by 1 and 15 times the respective period on the positive and negative branch of the cross-correlations (white dashed lines in Figures 2a-2c), hence using the same number of cycles independent of frequency. In addition, we only consider individual  t t / time series associated with a normalized Lomb-Scargle amplitude of at least 0.15, that is,  (Figures 2a-2c), to focus on coda waves that are subject to seasonal variations. Figures 3a-3c shows that the seasonal pattern is most consistent for the lower frequencies, whereas above about 8 Hz, especially component EZ shows a different behavior. Taking the average over the frequency range of 2-8 Hz (Figure 3d) reveals a similar pattern for all three components with high velocities (high values in  t t / ) in the winter months and low velocities in the summer months. In addition, the time series exhibit a long-term velocity decrease. The results across the components are most consistent after the installation of broadband sensor in 2017, likely due to the sensor itself (lower self noise, Figure 1d), and the temporal proximity to the reference window. To investigate potential artifacts due to cycle skipping when using the fixed reference, Figures 3e-3h shows the results for the moving reference approach (same processing otherwise). In this case, the seasonal pattern emerges more consistently over a broader frequency range, extending beyond 10 Hz. However, the 2-8 Hz averaged  t t / curves are similar as those for the fixed reference and show the same characteristics, that is, high (low) velocities in winter (summer) months and a long-term velocity decrease. While in the case of the fixed reference, component EN is associated with an increased standard deviation (SD; purple shading in Figure 3d), using the moving reference results in increased uncertainty in component EZ. Differences of both reference approaches are present at frequencies above 8 Hz, where only the moving reference shows seasonal variations, which we attribute to the elimination of cycle skipping. However, careful inspection of Figures 3e-3g also reveals some high-velocity notches in summer, best visible on component EZ in 2011EZ in , 2012EZ in , 2013EZ in , and 2015 above 10 Hz. These features are also visible in Figure 2d (fixed reference), where component EZ shows summer drops in  t t / that overshoot the lower cycle skipping limit and subsequently enter the plot again as steep lines from the upper cycle skipping limit. As using the moving reference (in 2016) does not eliminate this type of cycle skipping suggests that some years exhibit changes relative to 2016 that exceed the cycle skipping limit.

Velocity Changes
In most seismic monitoring applications,  t t / is obtained from linear regression assuming a bulk velocity change. Here, we find that only parts of the coda waves at different lag times show clear seasonal changes (Figures 2a-2c), which suggests localized changes and we consequently consider  t t / as a proxy for dv v / . To attach numbers to the seasonal and long-term changes, we fit a velocity change model consisting of the superposition of a sinusoid with a period of 365.25 days and a linear trend to the  t t / time series (Equation S1). Using the 2-8 Hz curves and averaging over the three components yields seasonal peakto-peak velocity changes of 3.3% for the fixed reference and 2.9% for the moving reference and long-term velocity decreases of 0.14%/yr and 0.11%/yr, respectively. In addition to different reference approaches, we further examine different strategies in determining  t t / for component NZ (most consistent component), including linear regression analysis and the classical MWCS technique (see Text S2 and Figure S2 for details). In all cases, we find seasonal velocity changes with high (low) values in late winter (summer) and a velocity decrease over the past 15 yr. Yet, the amplitudes of both characteristics are smaller when using the whole coda wave window independent of the 365.25-day periodicity, which further supports that the velocity changes are localized. The observed 2-8 Hz velocity change characteristics appear to be present also at higher frequencies (Figure 3), however, we refrain from analyzing a frequency dependence, as we are facing a complex setting with steep terrain and heterogeneously distributed medium changes, which is far off from a layered half space typically assumed for example, in surface wave analysis (James et al., 2019). In addition, we encounter cycle skipping at higher frequencies, which partly remains even when using a moving reference. This could hint toward large velocity drops (several percent) between consecutive years due to irreversible thawing of frozen rock. However, this remains speculation and would need to be addressed by an even more rigorous moving reference approach as applied for example, by James et al. (2019), to rule out artifacts of our year-by-year comparison.
The location of the velocity changes can be examined with travel time sensitivity kernels for a colocated source and receiver. In this case, the sensitivity kernels for both of the two end-member scenarios of single scattering (Pacheco & Snieder, 2006) and multiple scattering (Pacheco & Snieder, 2005) peak at the station location and decrease rapidly with distance from the station (Bennington et al., 2021;Sens-Schönfelder & Wegler, 2006). This implies that the velocity changes relate to the direct surroundings of the station (Hobiger et al., 2014) and we hypothesize that they are caused by thaw and refreeze associated with the permafrost lens documented in Krautblatter et al. (2010). The permafrost lens is separated by about 200 m from the station, which is only a fraction of one wavelength at 2 Hz and on the order of one wavelength at 8 Hz (p-wave rock matrix velocities of samples from Mt. Zugspitze measured in the laboratory range between about 2,000 and 6,000 m/s (Draebing & Krautblatter, 2012)). However, we also note that the single station sensitivity kernels may not comprehensively describe the encountered situation with a stationary noise source at the Zugspitze summit in some distance to the station. This situation also admits phases resulting from the cross-correlation of direct waves emitted from the noise source and their singly scattered (laterally or at depth) products, hence source and receiver being not colocated. This is expected to add travel time sensitivity also to the noise source region and the direct path between source and receiver, similar as for two-station cross-correlation sensitivity kernels (Obermann et al., 2019). While parts of these sensitivity areas are affected by permafrost, parts are also subject to seasonal freezing, only. To further pinpoint the velocity changes, we note, that denser instrumentation would be necessary.

Environmental Data and Permafrost Dynamics
We evaluate our hypothesis of freeze-thaw induced velocity changes by considering the recordings from a weather station at the Zugspitze summit run by the Deutscher Wetterdienst (DWD, German weather service). Figure 4 shows the 2-8 Hz velocity changes (fixed reference, same as Figure 3d), as well as the air temperature, snow height, and fluid precipitation measured by the weather station. This reveals that the annual velocity drops starting in April (vertical gray solid lines) occur concurrently with air temperatures rising above the freezing point. This especially holds when adding an offset of around 1 °C to the temperature curve to account for the elevation difference between the summit and the ridge (assuming an atmospheric lapse rate of around −0.6 °C/100 m). Minimum annual velocities are reached in July and August, followed by a velocity increase starting in September, coincidentally with temperatures dropping below the freezing point (vertical gray dashed lines). With temperatures above the freezing point, the period between April and September, where the velocity decreases, is furthermore characterized by snow melt and rain-dominated precipitation (Figure 4c). Considering the long-term trend, the velocity drops on the order of 0.1%/yr, while the temperature rises on average by 0.07 °C/yr in the time period between January 1, 2006 and January 1, 2021. The determined linear trends (using Equation S1) are depicted by the colored dashed lines in Figures 4a and 4b.
The permafrost lens in the ridge to the north of the seismic station ( Figure 1) is monitored by time-lapse temperature-calibrated ERT images (Scandroglio, Rehm, et al., 2021), of which results are documented for 2007 (Krautblatter et al., 2010). Krautblatter et al. (2010) observe pronounced melt from May to August with rock temperature changes being too fast to be solely explained by heat conduction. Coincidentally, they observe water seepage into the gallery and rapid melting along fracture zones suggesting warming and thawing through water percolation. With temperatures dropping below the freezing point in September, the ERT results of Krautblatter et al. (2010) show refreezing from the rock face. Similarly to the electrical resistivity, seismic velocities are different for frozen and unfrozen material and sharply increase at the freezing point (King et al., 1988;Kneisel et al., 2008;Leclaire et al., 1994). Laboratory experiments including samples from Mt. Zugspitze show that this also holds for low-porosity Wetterstein limestones representative for the study site, where p-wave velocity increased by 75% parallel to cleavage upon freezing (Draebing & Krautblatter, 2012). The seasonal velocity changes can thus be explained by the annual heat wave causing progressive thawing to depth starting in spring from the rock face, which will decrease seismic velocities. The observed rapid decline of velocities is presumably enhanced by water percolation from melt and precipitation. Once temperatures drop again below the freezing point in fall, progressive refreezing from the rock face to depth will again increase the velocities (a zoom into a 3 yr period is provided in Figure S3). The immediate response of the velocity to the temperature dropping below and rising above the freezing point (fall and spring, respectively), appear plausible in the light that centimeter-scale ground freezing is sufficient to result in significant surface waves velocity changes (Steinmann et al., 2021). In addition, the seasonal velocity changes are comparable in amplitude with those of freeze-thaw cycles extracted through seismic interferometry in other studies (Guillemot et al., 2020;James et al., 2019;Steinmann et al., 2021). Finally, following the argumentation line from above, the long-term decrease in seismic velocities can be well explained by permafrost degradation, that is, the shrinkage of the perennially frozen rock volume due to rising temperatures. Considering that the seasonal freeze-thaw depths are on the meter scale (Gallemann et al., 2017;Nötzli et al., 2010) and that the seasonal velocity variations are at least one magnitude greater than the long-term velocity decrease, we expect permafrost degradation on a centimeter scale per year. This is in agreement with borehole temperature logging beneath the summit, showing dominant seasonal variations superimposed on a slower permafrost degradation on a centimeter to decameter scale per year (Gallemann et al., 2017). In contrast to air temperature, the annual precipitation (Figure 4c) does not show a long-term increase and is thus not a plausible explanation for the long-term velocity decrease.

Conclusions
Using passive seismic data from Mt. Zugspitze, we find seasonal seismic velocity changes as well as a velocity decrease over the observation period of 15 yr. Comparison of our results with meteorological data and previous ERT and borehole studies, suggest that these velocity changes are caused by seasonal freeze-thaw cycles and permafrost degradation, respectively. Although originally deployed for earthquake monitoring, we were able to exploit the seismic station for long-term permafrost monitoring yielding velocity change values on more than 80% of all days in the 15-yr observation period. This highlights the cost and labor efficient potential of seismology for continuous permafrost monitoring, compared to other methods where long-term monitoring is challenged by manual data acquisition requiring regular field trips. Yet, the single station approach of this study is limited in the spatial resolution. Future studies should therefore extend the instrumentation in order to investigate permafrost dynamics with high spatio-temporal resolution. In this context, the recently introduced distributed acoustic sensing systems, which allow wave propagation sensing on a meter scale along fiber-optic cables are promising for detailed permafrost monitoring.

Data Availability Statement
Data of station BW.ZUGS are available through the EIDA node at http://erde.geophysik.uni-muenchen. de/fdsnws, data of the meteorological station through the DWD Climate Data Center (https://cdc.dwd.de/ portal/, station ID 5792). The scripts to reproduce the figures of this manuscript are available at https://osf. io/jkwat/.