Long-term analysis of cryoseismic events and associated ground thermal stress in Adventdalen, Svalbard

The small-aperture Spitsbergen seismic array (SPITS) has been in continuous operation at Janssonhaugen for decades. The high-Arctic location in the Svalbard archipelago makes SPITS an ideal laboratory for the study of cryoseisms, a nontectonic class of seismic events caused by freeze processes in ice, ice–soil, and ice–rock materials. We extracted a catalog of > 100 000 events from the nearly continuous observation period between 2004 and 2021, characterized by short-duration ground shaking of just a few seconds. This catalog contains two main subclasses where one subclass is related to underground coal-mining activities and the other is inferred to be dominated by frost quakes resulting from thermal-contraction cracking of ice wedges and crack-filling vein ice. This inference is supported by the correspondence between peaks in observed seismicity with peaks in modeled ground thermal stress, based on a Maxwellian thermo-viscoelastic model constrained by borehole observations of ground temperature. The inferred frost quakes appear to be dominated by surface wave energy and SPITS proximal source positions, with three main areas that are associated with dynamic geomorphological features, i.e., erosional scarps and a frozen-debris/solifluction lobe. Seismic stations providing year-round, high-temporal-resolution measurements of ground motion may be highly complementary to satellite remote sensing methods, such as InSAR (interferometric synthetic aperture radar), for studying the dynamics of periglacial environments. The long-term observational record presented in this study, containing tens of thousands of cryoseismic events, in combination with a detailed record of borehole ground temperature observations, provides a unique insight into the spatiotemporal patterns of cryoseisms. The observed patterns may guide the development of models that can be used to understand future changes to cryoseismicity based on projected temperatures.

1 Introduction

Cryoseisms are a nontectonic class of seismic events caused by freeze processes in ice, ice–soil, and ice–rock materials (Lacroix, 1980). For example, the buildup of thermal stress in frozen soils during intense periods of cooling can lead to tensional fracturing and explosive stress release (Barosh, 2000; Battaglia et al., 2016). Since this stress release triggers the propagation of seismic waves, these events are sometimes referred to as frost quakes (e.g., Okkonen et al., 2020). Observations of frost quakes are limited because seismic amplitudes decay rapidly with distance from the point of rupture, but they have been felt at distances of several hundred meters to several kilometers and are usually accompanied by cracking or booming noises, resembling falling trees, gunshots, or underground thunder (Leung et al., 2017; Nikonov, 2010). Frost quakes have typically been observed in association with rapid air and ground cooling, in the absence of an insulating snow layer and where sufficient moisture is present for ice to form (Barosh, 2000; Battaglia et al., 2016; Matsuoka et al., 2018; Nikonov, 2010). Ice can initiate ground cracking by different mechanisms including thermal contraction and tensile failure (e.g., Lachenbruch, 1962) or by the growth of segregation ice driven by the capillary migration of water to a freezing front (Peppin and Style, 2013; Walder and Hallet, 1985). The resulting fractures may potentially cause damage to buildings and other infrastructure in cold regions (Okkonen et al., 2020). Frost damage can be understood as a combination of slow creep (frost heave, e.g., Rempel, 2010) and rapid elastic (frost quake) deformation of frozen ground and causes damage to roads requiring billions of dollars annually to repair in the United States alone (DiMillio, 1999).

Figure 1 Modified after the Lachenbruch (1962) model of ice wedge formation. (a) Thermal contraction of the frozen active layer initiates a tensional crack that penetrates into permafrost when the stress, σ, exceeds the tensile strength of frozen ground, σ_T(ground). (b) Meltwater infiltrates and refreezes during the thaw season. (c) Ice with lower tensile strength, σ_T(ice), than surrounding ground forms a plane of weakness and cracks repeatedly over many years. (d) The crack–infill–freeze cycle causes ice wedge growth in the permafrost and ground surface deformation that organizes into ice wedge polygons in a two-dimensional plan view (e.g., Plug and Werner, 2002).

In some cases, frost quakes are involved in forming various geomorphic features on the surface. For example, ice wedge or sand wedge polygons are a widely observed geomorphic feature in the periglacial environment (e.g., Black, 1976; Matsuoka et al., 2004). These wedges form when water that infiltrates and freezes to ice or wind-transported sand grains hold open an initial thermal-contraction crack that subsequently becomes an enduring plane of structural weakness (Lachenbruch, 1962; Mackay, 1984; Matsuoka et al., 2004; Plug and Werner, 2002; Sørbel and Tolgensbakk, 2002). As shown in Fig. 1, repeated cracking, infilling, and refreezing along these planes of weakness causes the ice wedges to grow laterally in permafrost. This forces the displaced ground upwards and results in a series of troughs and ridges in a polygonal arrangement that are one of the most recognizable landforms in permafrost environments (Christiansen et al., 2016; Lachenbruch, 1962; Matsuoka et al., 2018; Plug and Werner, 2002). For this reason, sand wedge polygons, which were formed at sea level on the paleo-equator during late Neoproterozoic glacial episodes, have been used as evidence supporting the snowball Earth hypothesis (Maloof et al., 2002). Small-scale polygonal features observed on the surface of Mars have also been inferred to result from thermal-contraction cracking (Mellon, 1997).

Frost cracking driven by segregation ice growth is also an important agent of bedrock erosion in cold mountainous areas, where rockfall, active screes, and high headwall erosion rates are observed in areas where frost action is most intense (Hales and Roering, 2009, 2007; Scherler, 2014). An important mode of crack growth in water permeable bedrock is the migration of water to form segregation ice bodies (Hallet et al., 1991; Murton et al., 2006; Walder and Hallet, 1985), which is similar to the mechanism by which ice lenses develop in freezing soil (Peppin and Style, 2013). On slopes, segregation ice growth, frost heaving, and creep lead to the development of solifluction lobes and sheets (Cable et al., 2018; Matsuoka, 2001). Solifluction is broadly defined as the slow mass wasting resulting from freeze–thaw action in fine-textured soils (French, 2017; Matsuoka, 2001) and occurs due to the asymmetry between frost heaving perpendicular with the sloped ground surface and vertical subsidence upon thawing under the force of gravity.

Spatiotemporal records with high resolution are essential to better understand the processes mentioned above and predict their evolution, for example, in a warming climate. Interferometric synthetic aperture radar (InSAR) is a satellite remote sensing technique capable of resolving down to millimeter scale ground surface displacements at a spatial scale of tens of meters (Hanssen, 2001; Rosen et al., 2000). Consequently, InSAR has been used to resolve seasonal patterns of frost heave and thaw subsidence (Chen et al., 2020; Rouyet et al., 2021; Wahr et al., 2008). The broad spatial coverage and high spatial resolution of InSAR have also been used to show that the spatial patterns of ground displacement are related to specific geomorphological units, consisting of specific frost-prone sediment types and periglacial landforms (Liu et al., 2018; Rouyet et al., 2019). On the other hand, the temporal resolution is determined by the repeat cycle of the satellite (∼ 6 d for Sentinel-1), which means that it is not possible to distinguish between slow frost creep and rapid elastic-deformation (frost quake) processes. In addition, snow cover causes a loss of radar coherence between repeat satellite passes that limits InSAR to snow-free areas and seasons. As a result, InSAR is more suited to studying thaw subsidence than frost heave in areas like Svalbard that are snow covered for the majority of the freezing season (e.g., Rouyet et al., 2019). Seismic stations providing year-round, high-temporal-resolution, and highly sensitive measurements of ground motion, at the expense of more limited spatial source resolution and coverage, may be highly complementary to InSAR methods for studying the dynamics of periglacial environments.

The Arctic archipelago of Svalbard is an ideal place to study these processes, given its rapidly warming climate and strong multi-decadal focus on environmental research that has yielded some unique observational datasets. This study was motivated by the intermittent observation of clusters of events recorded by the small-aperture Spitsbergen seismic array (SPITS), located on the main island of the archipelago. These seismic events had durations of just a few seconds (indicative of a nearby source) and peak amplitudes significantly above the background noise level. For comparison, regional tectonic earthquakes are typically associated with a ground-shaking duration of > 30 s, since the different velocities of P, S, and surface waves cause the wave field to spread out (dispersion) as the source distance increases. SPITS has previously been used to study a class of cryoseismic signals associated with iceberg calving at the termini of grounded tidewater glaciers at local to near-regional distances of up to 200 km (Köhler et al., 2015). These calving-related seismic signals occurred most frequently during the melt season, and the ground motion lasts ∼ 15–20 s. Signals with intermediate durations have also been observed at SPITS, originating from nearby mountain glaciers (Albaric et al., 2021) and coal-mining operations (Gibbons and Ringdal, 2006).

Based on previous work further down-valley in Adventdalen, e.g., Matsuoka et al. (2018) and Romeyn et al. (2021), we hypothesized that the short-duration events at SPITS might be dominated by frost quakes initiated by thermal-contraction cracking in the local vicinity of the array. The aim of this study was to test this hypothesis by analyzing the spatial and temporal occurrence of these events. We focus on the thermal-contraction cracking mechanism because it is consistent with the association of transient ground acceleration events with rapid cooling episodes and cold winter temperatures reported by Matsuoka et al. (2018) and consistent with previous descriptions of frost quakes (Barosh, 2000; Battaglia et al., 2016; Nikonov, 2010; Okkonen et al., 2020). This study is highly complementary to the previous work of Romeyn et al. (2021). While the spatial and temporal wave field sampling of the SPITS array is much coarser than the temporary geophone array they deployed, a much longer and nearly continuous record is available. This allows for a more rigorous investigation of the temporal correlation of these events with ground cooling and thermal-stress accumulation.

2 Study area and data

The small-aperture Spitsbergen seismic array (SPITS) is located on Janssonhaugen, in Adventdalen on the island of Spitsbergen, part of the high-Arctic Svalbard archipelago (Fig. 2). The SPITS array has been in operation since 1992, maintained by the research foundation NORSAR (Norwegian Seismic Array; Schweitzer et al., 2021). At present, it consists of nine CMG-3T seismometers with an aperture of 1 km and interstation distances > 250 m (e.g., Gibbons et al., 2011; Köhler et al., 2015) installed in shallow boreholes below the permafrost active layer. The standard frequency response of the CMG-3T seismometer is flat (within −3 dB) for the range from 0.0083 to 50 Hz. We chose to limit the present study to the period following August 2004 when the SPITS array was upgraded to a full broadband array with an increase in sampling rate from 40 to 80 Hz for all seismometers (Schweitzer et al., 2021). The waveform data following the upgrade are of high quality and well suited to source localization using matched-field processing. Waveform data for the SPITS array were retrieved from the European Integrated Data Archive (EIDA), maintained by the University of Bergen and NORSAR.

Figure 2 Overview of the location of the SPITS array located on Janssonhaugen in Adventdalen, Svalbard, up-valley from the settlement of Longyearbyen and the temporary seismic array of Romeyn et al. (2021). The P11 temperature logging borehole and the underground mining areas and tunnels of the operational coal mine, Gruve 7, are also shown. The http://met.no (14 January 2022) weather station “Janssonhaugen Vest” is marked by its station number, SN99874. Orthophoto © Norwegian Polar Institute (http://npolar.no, 17 January 2022).

Janssonhaugen is a bedrock remnant located in the middle of Adventdalen, with a ∼ 0.2–0.3 m thick weathered sediment crust overlying homogeneous sand-/siltstone bedrock (density of 2280 kg m⁻³, porosity of 20 %–25 %, > 96.5 % SiO₂) corresponding to the Ullaberget Member of the Lower Cretaceous Rurikfjellet Formation (Dypvik et al., 1991; Isaksen et al., 2001). The bedrock at Janssonhaugen has been drilled to a depth of 102 m, and the permafrost zone is estimated to extend down to 220 m, based on downward extrapolation of the measured temperature gradient (Isaksen et al., 2001). Despite annual precipitation of 300–500 mm yr⁻¹, the snow cover on Janssonhaugen is typically thin or completely absent due to the scouring effect of the prevailing winds (Isaksen et al., 2001). The surface topography is generally flat, but loose surface material is sorted into polygons (Isaksen et al., 2001), due to active-layer cryoturbation and/or ice wedge formation. The homogeneous geology at Janssonhaugen; its flat topography; limited snow cover; relatively large distance to glaciers, rivers, ocean, and human activity such as coal mining; and position well above the Holocene marine limit make it a good location to study permafrost processes (Isaksen et al., 2000).

Figure 3 Illustration of spatiotemporal borehole temperature recorded at the PACE (Permafrost and Climate in Europe) P11 borehole on Janssonhaugen. A long-term warming trend is observed in the permafrost below the ∼ 2 m thick active layer that is subject to a seasonal freeze–thaw cycle. The continuous timeline is split across multiple figure panels. The temperature–depth profile was interpolated to regular 10 cm intervals using a spline interpolant.

At this location, the polar night, where no shortwave solar radiation is received at the ground surface, lasts from around 1 October to 28 February each year but is nonetheless the season associated with the largest diurnal temperature range (Przybylak et al., 2014). This is because intense winter storms promote the advection of warmth to the region, driven 95 % by atmospheric circulation and 5 % by oceanic circulation (Bednorz, 2011). Snow cover and latent heat fluxes also contribute to a complex surface energy budget (Westermann et al., 2009). As a result of the complex interplay of processes driving temperature variation in this high-Arctic location, temperature-monitoring boreholes installed at Janssonhaugen are critical to our ability to accurately model the subsurface buildup of thermal stress. Thermistors installed in the GTN-P (Global Terrestrial Network for Permafrost) boreholes P10 (102 m deep) and P11 (15 m deep) provide a continuous record of the subsurface temperature field at Janssonhaugen, with a sampling interval of 6 h extending from April 1999 to the present (Isaksen et al., 2001). We focus on the P11 borehole (see Fig. 2), which was drilled 13 m horizontally offset from P10, gives the most detailed record of the near-surface temperature field and is least disturbed by installations at the ground surface (Isaksen et al., 2000). The temperature field measured by the P11 borehole is illustrated in Fig. 3. The ∼ 2 m thick active layer (Christiansen et al., 2020) is sampled by thermistors at 0.2, 0.4, 0.8, 1.2, 1.6 and 2 m, and there is significant interannual variability in the magnitude of summer warming and winter cooling. The upper part of the permafrost at the P11 borehole location is sampled by thermistors installed at 2.5, 3, 3.5, 4, 5, 7, 10, 13 and 15 m. These measurements record a long-term warming trend in the permafrost beneath the active layer (see Fig. 3). Furthermore, the Janssonhaugen Vest weather station (see Fig. 2), which was installed in September 2019, includes hourly sampled records of air temperature and ground temperature at 0.1 m depth. It therefore provides a basis to compare depth and temporal sampling effects against the longer-duration, more coarsely sampled P11 borehole record.

3 Methods

3.1 STA / LTA detection of short-duration seismic events

Events are detected based on anomalous values of short-time-averaged amplitude (STA) divided by long-time-averaged amplitude (LTA), i.e., the classic STA / LTA approach widely used in seismology (e.g., Allen, 1982; Trnkoczy, 2009). The purpose of the event detector is to make an initial, coarse, automatic identification of short-duration seismic signals, which should be distinguished from both background noise and longer-duration local and regional seismic events that may be high amplitude. The raw data from each seismometer were first detrended, corrected according to the calibrated instrument sensitivity, and bandpass-filtered to the range 2.5–20 Hz using a delay-compensated minimum-phase filter with a stopband attenuation of 60 dB (see Fig. 4a). This passband was selected to eliminate high-frequency random noise which can lead to spurious spikes in STA / LTA.

The STA window length should be comparable to the target signal duration, while the LTA represents the background noise level. When the STA / LTA ratio exceeds a given threshold, an event is triggered. In our implementation, the STA is given by the 1 s moving-average smoothed trace envelope for each seismogram. The LTA is the STA further smoothed according to a 20 s period moving average. Since we have an array of stations, we represent the array STA by taking the 80th percentile of the station STA across all stations at a given time sample. If we had chosen the maximum station STA, we may detect arbitrary noise spikes with large amplitudes registered on a single seismometer. If the mean or median station STA were chosen, we would preferentially detect larger regional events with more consistent amplitudes across the array and suppress smaller local events with high amplitudes limited to a small subset of seismometers.

The STA / LTA threshold was set to 10. Furthermore, after a trigger, no new events are declared within 5 s after the STA / LTA threshold was exceeded to avoid detecting the same events multiple times. All processing parameters (STA and LTA lengths, STA / LTA threshold, STA percentile, and trigger pause) were found to be appropriate for detection of short-duration signals coherent across the array, while avoiding false triggers (noise bursts at single stations), by visually inspecting test periods. We further reject epochs with LTA more than 2.5 times the 2 h mean of the LTA in order to filter out large regional earthquakes. An example of event detection based on STA / LTA peaks is shown in Fig. 4b, which illustrates that short-duration events are selected, while a longer-duration high-amplitude event is not detected as intended.

Figure 4 Example of (a) a time series of vertical-component ground motion and (b) event detection using the STA / LTA detector. Short-duration events with sufficient amplitude and array coherence are selected, while longer events such as the high-amplitude example at 04:50 UTC are ignored. Inset boxes show detailed views for a specific detection. Please note that the date format used in this figure is month day, year.

3.2 Source localization by coherent MFP

Matched-field processing (MFP) is an established technique for localizing the position of seismic sources recorded by passive seismic arrays (Chmiel et al., 2016; Cros et al., 2011; Harley and Moura, 2014; Sergeant et al., 2020; Walter et al., 2015). MFP proceeds via an evaluation of the similarity between the wave field recorded at a receiver array and a series of predicted wave fields calculated for a grid of test source locations and a theoretical model of the source–receiver wave propagation. The MFP coherence is estimated by comparing the recorded data vector with a “replica” vector, which is a model representation of wave propagation within the medium. We assume a simple homogeneous medium with amplitude decay according to spherical divergence where the replica column vector R^′ is represented by the theoretical harmonic wave emitted from a test point p(x,y) according to

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathbf{R}}^{\prime }\left(\mathit{\omega },\mathit{d}\right)=\\ & {\left[{\displaystyle \frac{\mathrm{1}}{d_{\mathrm{1}}}}{e}^{i\mathit{\omega }{d}_{\mathrm{1}}/c\left(\mathit{\omega }\right)},{\displaystyle \frac{\mathrm{1}}{d_{\mathrm{2}}}}{e}^{i\mathit{\omega }{d}_{\mathrm{2}}/c\left(\mathit{\omega }\right)},\mathrm{\cdots },{\displaystyle \frac{\mathrm{1}}{d_N}}{e}^{i\mathit{\omega }{d}_N/c\left(\mathit{\omega }\right)}\right]}^T/\mathit{\eta },\end{array}\end{array}$$

where $i=\sqrt{-\mathrm{1}},\mathit{\omega }$ is the angular frequency; $\mathit{d}={\left[d_{\mathrm{1}},d_{\mathrm{2}},\mathrm{\cdots },d_N\right]}^{\mathrm{T}}$ is an N×1 vector containing the absolute Euclidean distances between p(x,y) and the N recording stations, where superscript T denotes a transpose; c(ω) is the medium-phase velocity at frequency ω; and the vector is normalized to unit length by dividing by the factor $\mathit{\eta }=\sqrt{{\sum }_{j=\mathrm{1}}^N\left(\mathrm{1}/d_j^{\mathrm{2}}\right)}$. The column data vector is given by

$$\begin{array}{}\text{(2)}& \mathbf{R}\left(\mathit{\omega }\right)={\left[R_{\mathrm{1}}\left(\mathit{\omega }\right),R_{\mathrm{2}}\left(\mathit{\omega }\right),\mathrm{\cdots },R_N\left(\mathit{\omega }\right)\right]}^{\mathrm{T}},\end{array}$$

where R_j(ω) is a frequency transform of the jth trace r_j(t) of N traces recording a specific seismic event. We then form the complex-valued N×N cross-spectral density matrix (CSDM) by

$$\begin{array}{}\text{(3)}& \mathbf{K}\left(\mathit{\omega }\right)=\mathbf{R}\left(\mathit{\omega }\right){\mathbf{R}}^{\mathrm{H}}\left(\mathit{\omega }\right),\end{array}$$

where (⋅)^H denotes the Hermitian (conjugate transpose) operator. The frequency domain transform was implemented via the chirp z transform (Rabiner et al., 1969), which provides a convenient means to evaluate the band-limited transform. In this study we selected the 5–35 Hz band, sampled at an interval of 1 Hz, in order to reject spatially coherent low-frequency background noise while retaining the shared signal–noise high-frequency band (Appendix D shows that the results are relatively insensitive to the selected frequency band). This allows us to efficiently and compactly represent the CSDM, whose size is a significant factor in the speed of computation, in addition to the number of test source points of p(x,y). Conventionally, the MFP coherence is estimated using the linear Bartlett processor as

$$\begin{array}{}\text{(4)}& G={\sum }_{\mathit{\omega }}{\mathbf{R}}^{{}^{\prime }\mathrm{H}}\left(\mathit{\omega },\mathit{d}\right)\mathbf{K}\left(\mathit{\omega }\right){\mathbf{R}}^{\prime }(\mathit{\omega },\mathit{d}),\end{array}$$

which evaluates the inner product between the recorded and predicted wave fields before summing incoherently across frequency. The matrix operations in Eq. (4) are quadratic forms that formally guarantee a non-negative real-valued output.

In this study we implement the coherent-MFP scheme developed by Michalopoulou (1998), which is an elegant addition to the conventional approach, allowing for cross-frequency spatial coherence structures to be exploited to give improved robustness and accuracy. In this scheme, measurement vectors at L discrete frequencies are concatenated to form the NL×1 measurement super-vector of

$$\begin{array}{}\text{(5)}& \mathbf{R}\left(\mathit{\omega }\right)={\left[\mathbf{R}\left({\mathit{\omega }}_{\mathrm{1}}\right),\mathbf{R}\left({\mathit{\omega }}_{\mathrm{2}}\right)\mathrm{\cdots },R\left({\mathit{\omega }}_L\right)\right]}^{\mathrm{T}},\end{array}$$

where $\mathit{\omega }={\left[{\mathit{\omega }}_{\mathrm{1}},{\mathit{\omega }}_{\mathrm{2}},\mathrm{\cdots },{\mathit{\omega }}_L\right]}^{\mathrm{T}}$ is a frequency vector and the replica vectors are concatenated to form the NL×1 replica super-vector of

$$\begin{array}{}\text{(6)}& {\mathbf{R}}^{\prime }\left(\mathit{\omega },\mathit{d}\right)={\left[{\mathbf{R}}^{\prime }\left({\mathit{\omega }}_{\mathrm{1}},\mathit{d}\right),{\mathbf{R}}^{\prime }\left({\mathit{\omega }}_{\mathrm{2}},\mathit{d}\right),\mathrm{\cdots },{\mathbf{R}}^{\prime }\left({\mathit{\omega }}_L,\mathit{d}\right)\right]}^{\mathrm{T}}.\end{array}$$

The super-CSDM K(ω)=R(ω)R^H(ω) is then composed of NL×NL elements, and the generalized MFP coherence is given by

$$\begin{array}{}\text{(7)}& \mathcal{G}={\mathbf{R}}^{{}^{\prime }\mathrm{H}}(\mathit{\omega },\mathit{d})\mathbf{K}\left(\mathit{\omega }\right){\mathbf{R}}^{\prime }(\mathit{\omega },\mathit{d}).\end{array}$$

The estimated source position of p(x,y) and phase velocity of c(ω) are those which maximize the coherence measure 𝒢. Since some of the seismic events we wish to locate are dominated by surface waves and others are dominated by body waves, we simply assume that phase velocity is a constant and scan over the range 250–6000 m s⁻¹.

3.3 Ground thermal-stress model

Previous studies such as Lachenbruch (1962), Mellon (1997), Maloof et al. (2002), Schulson and Duval (2009), and Podolskiy et al. (2019) have demonstrated that ice and frozen ground deform elastically on short timescales and viscously on long timescales. Thermal loading due to temperature changes acts as an external driving agent, and the resulting dynamical balance between the elastic and viscous response governs whether creep or fracture becomes dominant. Following Mellon (1997), we neglect the layered structure of the ground and model the frozen ground as a simple homogeneous Maxwellian viscoelastic solid augmented with thermal expansion and contraction. This allows us to decompose the total strain tensor ε_ij into three components, an elastic (${\mathit{\epsilon }}_{ij}^{\mathrm{e}}$), a thermal (${\mathit{\epsilon }}_{ij}^{\mathrm{T}}$), and a viscous $\left({\mathit{\epsilon }}_{ij}^{\mathrm{V}}\right)$ component, as

$$\begin{array}{}\text{(8)}& {\mathit{\epsilon }}_{ij}={\mathit{\epsilon }}_{ij}^{\mathrm{e}}+{\mathit{\epsilon }}_{ij}^{\mathrm{T}}+{\mathit{\epsilon }}_{ij}^{\mathrm{V}},\end{array}$$

where subscripts ij indicate tensor components of $i,j=\mathrm{1},\mathrm{2},\mathrm{3}$, where 1 and 2 denote horizontal components and 3 is the vertical component. The elastic strain tensor is related to the stress tensor σ_ij by (e.g., Landau and Lifshitz, 1970)

$$\begin{array}{}\text{(9)}& {\mathit{\epsilon }}_{ij}^{\mathrm{e}}={\displaystyle \frac{\mathrm{1}+\mathit{\nu }}{E}}{\mathit{\sigma }}_{ij}-{\displaystyle \frac{\mathit{\nu }}{E}}{\mathit{\sigma }}_{kk}{\mathit{\delta }}_{ij},\end{array}$$

where ν is Poisson's ratio, E is Young's modulus, and δ_ij is the Kronecker delta; Einstein's summation convention is applied throughout this paper.

The thermal strain tensor is a measure of the change in volume caused by the thermally driven deformation, and it is expressed as (e.g., Landau and Lifshitz, 1970)

$$\begin{array}{}\text{(10)}& {\mathit{\epsilon }}_{ij}^{\mathrm{T}}=\mathit{\alpha }\left(T-T_{\mathrm{0}}\right){\mathit{\delta }}_{ij},\end{array}$$

where T is the temperature, T₀ is a reference temperature for the undeformed state, and α is the linear thermal-expansion coefficient. The viscous-strain tensor is a complicated topic in itself, and a wide range of phenomenological and heuristic parametric models for the viscous-strain rate exist (e.g., Bingham, 1922; Carreau, 1972; Glen, 1955; Herschel, 1926; Saramito, 2007). To encompass this generality, we formulate the viscous-strain rate as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial {\mathit{\epsilon }}_{ij}^{\mathrm{V}}}{\partial t}}={\mathrm{\Gamma }}_N\left\{s_{ij}\right\},\end{array}$$

where Γ_N{⋅} is a nonlinear operator acting on the deviatoric stress $s_{ij}={\mathit{\sigma }}_{ij}-{\mathit{\sigma }}_{kk}/\mathrm{3}$. The chosen parametric form of Γ_N{⋅} determines how induced stress relaxes in the medium, and both Newtonian and non-Newtonian behavior can be incorporated in this formulation.

Following Mellon (1997), we assume $\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{11}}={\mathit{\sigma }}_{\mathrm{22}}$ and $\partial {\mathit{\epsilon }}_{\mathrm{11}}/\partial t=\partial {\mathit{\epsilon }}_{\mathrm{22}}/\partial t=\mathrm{0}$. The vertical stress is also neglected, i.e., σ₃₃≈0, corresponding to the assumption that the overburden is negligible in the shallow subsurface (Mellon, 1997). We assume that Poisson's ratio ν=ν(T), Young's modulus E=E(T), and the coefficient of linear thermal expansion α=α(T) are all temperature dependent and, thus, implicitly time dependent, since $T=T(z,t)$ is the temperature at depth z and time t. By direct evaluation of $\partial {\mathit{\epsilon }}_{ij}/\partial t$ and collecting terms, we find that the temporal dynamics of horizontal stress σ(z,t) in a Maxwellian viscoelastic solid driven by thermal expansion and contraction is governed by the following first-order nonlinear and nonhomogeneous differential equation:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial \mathit{\sigma }}{\partial t}}+\mathit{\beta }\left(t\right)\mathit{\sigma }+\mathrm{\Gamma }\left\{\mathit{\sigma }\right\}=\mathit{\kappa }\left(t\right).\end{array}$$

The time-dependent coefficients are found to be

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\beta }\left(t\right)=-\left({\displaystyle \frac{\mathrm{1}}{E}}{\displaystyle \frac{\partial E}{\partial T}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\nu }}}{\displaystyle \frac{\partial \mathit{\nu }}{\partial T}}\right){\displaystyle \frac{\partial T}{\partial t}},\text{(14)}& {\displaystyle }\mathit{\kappa }\left(t\right)=-{\displaystyle \frac{E}{\mathrm{1}-\mathit{\nu }}}\left[\mathit{\alpha }+{\displaystyle \frac{\partial \mathit{\alpha }}{\partial T}}(T-T_{\mathrm{0}})\right]{\displaystyle \frac{\partial T}{\partial t}},\end{array}$$

and $\mathrm{\Gamma }\left\{\cdot \right\}=\frac{E}{\mathrm{1}-\mathit{\nu }}{\mathrm{\Gamma }}_N\left\{\cdot \right\}$. The scientific literature devoted to the rheology of frozen materials favor power law parametrizations for the viscous term (e.g., Schulson and Duval, 2009). In this paper, we apply the heuristic temperature-dependent power law proposed by Glen (1955) in the form used by Mellon (1997) and Maloof et al. (2002):

$$\begin{array}{}\text{(15)}& \mathrm{\Gamma }\left\{\mathit{\sigma }\right\}={\displaystyle \frac{E}{\mathrm{1}-\mathit{\nu }}}{A}_{\mathrm{0}}{\left|{\displaystyle \frac{\mathit{\sigma }}{\mathrm{2}}}\right|}^n\mathrm{sign}\left(\mathit{\sigma }\right)\exp \left(-Q/RT\right).\end{array}$$

In Glen's flow law R is the universal gas constant, and A₀, Q, and n are empirical parameters that need to be chosen. The temperature-dependent Arrhenius exponential term in this particular choice of Γ{⋅} is included to model the increasing ductility as temperature increases (e.g., Glen, 1955). To sum up, the first two terms on the left-hand side in Eq. (12) are connected to the elastic response of the solid; the third term models viscous relaxation; and the right-hand side is the thermal driving term. In order to solve Eq. (12) for σ(z,t), we specify the initial condition ${\mathit{\sigma }}_{\mathrm{0}}\left(z\right)=\mathit{\sigma }\left(z,t=\mathrm{0}\right)=\mathrm{0}$.

Table 1 Physical parameters used in the thermal-stress model.

If we assume $\partial \mathit{\nu }/\partial T=\mathrm{0}$, Eq. (12) reduces to the model proposed by Mellon (1997). By contrast, if we assume $\partial \mathit{\nu }/\partial T=\partial E/\partial T=\partial \mathit{\alpha }/\partial T=\mathrm{0}$, Eq. (12) reduces to the model proposed by Podolskiy et al. (2019). Finally, if we assume $\partial \mathit{\nu }/\partial T=\partial E/\partial T=\partial \mathit{\alpha }/\partial T=A_{\mathrm{0}}=\mathrm{0}$, Eq. (12) reduces to the Timoshenko and Goodier (1951) model for thermal stress as applied by Okkonen et al. (2020), excluding the boundary correction terms that represent compressive stresses at the free surfaces of the finite thickness plate assumed by the latter authors.

We solved Eq. (12) numerically to obtain a time series of the resulting horizontal stress σ(z,t), at depth z, using the standard MATLAB solver ode45, based on the well-known fifth-order Runge–Kutta method (Dormand and Prince, 1980). The crucial input to the forward model is the measured temperature time series T(z,t) at depth z. Notably, most previous studies of subsurface thermal stress have relied on thermal-conduction models in order to infer subsurface temperature variation from measurements of air temperature. A significant novelty of this study is that the ground temperature profile T(z,t) at Janssonhaugen has been logged by a series of thermistors installed in the 15 m deep P11 borehole at 6 h intervals since April 1999. The set of physical parameters that we assume is described in Table 1 and gives a homogenized representation of frozen regolith and soft sandstone interspersed with veins of ice that form planes of structural weakness.

Figure 5 (a, b) Examples of class I events with significant amplitude variation across the array and relatively close source position inferred by panels (d) and (e) as coherent MFP (Eq. 7). (c) Example of a class II event with little amplitude variation across the array and a more distal (f) coherent-MFP-inferred source position. Panels (g), (h), and (i) show corresponding incoherent-MFP results (Eq. 4), demonstrating the improvement gained by coherent MFP. Station numbers in panels (a), (b), and (c) correspond to the labels annotated in the MFP panels, where the seismometer locations are marked with magenta crosses.

3.3.1 Fracture model

Our aim is to model simple thermal tensional cracking of existing ice wedges or crack-filling vein ice and not the initiation or propagation of new cracks into previously undamaged ground, where pore scale fluid migration and stress localization at crack tips become important (e.g., Walder and Hallet, 1985). To this end, we apply a simplistic model of fracturing by considering the modeled thermal stress as the pre-fracture stress. When the pre-fracture stress exceeds the tensile strength of ice (see Table 1), a frost quake is assumed to occur and dissipate stress corresponding to 100 % of the tensile strength. We keep a tally of the number of frost quakes over time and assume that multiple frost quakes occur if the tensile strength is exceeded by an integer factor greater than one. In reality, a fraction of the stress would be redistributed elastically rather than completely lost to friction as in this simple model. Since ground temperature is a measured quantity, we do not account for warming by frictional heat as a result of frost quake movement.

Figure 6 Spatially binned (150 × 150 m) distribution of coherent-MFP-inferred seismic source positions plotted by season for events recorded between August 2004 and July 2021. Positions of SPITS seismometers are indicated by white crosses. The distal-event cluster corresponds with the location of underground mining operations at Gruve 7 (yellow polygon). Orthophoto © Norwegian Polar Institute (https://www.npolar.no/, 17 January 2022).

4 Results and discussion

A total of 137 532 short-duration seismic events were detected by our STA / LTA detector between July 2004 and June 2021. In order to improve the precision of the source localization by coherent MFP, only events recorded by at least five seismometers (max four inoperative stations) were located, for a total of 137 456 located short-duration events. The estimated source locations were subsequently used to identify subclasses of events, as detailed in the following section.

4.1 Subclasses of short-duration seismic events

We find that there are two main sub-classes of short-duration events recorded by the SPITS array. Event class I is characterized by significant amplitude variation and arrival time differences across the array seismometers, as illustrated in Fig. 5a, b. Using coherent MFP to infer the source positions of these events shows that they occur in relatively close proximity and are clustered inside a circle with a radius of ∼ 1500 m centered over the array (Fig. 5d and e). By contrast, event class II is characterized by similar amplitudes across the array elements and smaller relative arrival time differences, as illustrated in Fig. 5c. Using coherent MFP, we find that these events are associated with more distal inferred source positions (Fig. 5f) that also have a consistent azimuth. Figure 5 also illustrates the property that coherent MFP decreases source localization ambiguity for arrays that sample the spatial domain coarsely for a given range of observed wavelengths when compared to the incoherent scheme (consistent with Michalopoulou, 1998). We also observed that the well-calibrated amplitude responses of the SPITS seismometers give good constraint on source range under the assumed model of amplitude attenuation due to geometrical spreading (Fig. 5f). In Fig. 5i, we see that incoherent averaging (Eq. 4) has enhanced a sidelobe and produced an incorrect source position that is not consistent with the relative arrival times observed in Fig. 5c (earliest arrival at station 7 and latest arrivals with weakest amplitudes at stations 5 and 9).

Figure 7 Detailed view of Janssonhaugen overlaid with spatially binned (50 × 50 m) distribution of coherent-MFP-inferred seismic source positions plotted by season for events recorded between August 2004 and July 2021. Positions of SPITS seismometers are indicated by white crosses. Orthophoto © Norwegian Polar Institute (https://www.npolar.no/, 17 January 2022).

Figure 8 (a) December–February events as plotted in Fig. 7 with orthophotograph details illustrating the geomorphologic features associated with the most seismically active areas. Erosional scarps (b) and (d) are annotated in red with associated boulder fields in yellow, and frozen-debris/solifluction lobes (c) are annotated in blue. A faintly visible area of polygonally patterned ground is annotated in green in panel (c). Orthophoto © Norwegian Polar Institute (https://www.npolar.no/, 17 January 2022).

The mean MFP inferred propagation velocity for class I events was 1150 m s⁻¹ with a standard deviation of 1100 m s⁻¹, implying a relatively shallow propagation path. The large standard deviation may indicate the surface waves are dispersive with different frequencies propagating at different phase velocities. However, the signal for a given event was dominated by a relatively narrow band of frequencies rarely exceeding the 5 Hz bandwidth (see Appendix C), so the distinctive highly dispersive waveforms that one would expect for a wideband dispersive signal were not observed (Fig. 5). By contrast, the mean MFP inferred propagation velocity for class II events was 5750 m s⁻¹ with a standard deviation of 400 m s⁻¹, indicating that this event class is dominated by P-wave energy. Three-component seismograms also support the interpretation that event classes I and II are dominated by surface waves and P waves, respectively (see Appendix C).

Mapping the inferred source positions, it is clear that class II events are spatially coincident with the underground mining areas of Gruve 7, as illustrated in Fig. 6. Furthermore, we observe that class II events occur frequently during all seasons (see Fig. 6). We infer that event class II is a result of mining operations and human activity in the underground coal mine Gruve 7. This is also supported by the dominant P wave of the class II signals which indicates an explosive source. These events are essentially unwanted noise from the environmental seismology perspective, although they do give a useful indication of the localization performance of the coherent-MFP algorithm as applied in this study. The accuracy of the inferred source positions is largely a function of a data–model phase velocity mismatch and amplitude attenuation that deviates from the model. The error in the estimated source position will then be proportional to the amplitude-spectrum-weighted model–data phase velocity mismatch integrated over the signal bandwidth, convolved with the deviation from a purely geometrical spreading model of amplitude attenuation. It is very difficult to quantify this deviation a priori, so the best way to assess MFP accuracy is by locating known test sources. A set of test sources covering the entire study domain would have been ideal, but the fact that the seismic events associated with mining activities at Gruve 7 show reasonable spatial correspondence with the location of the underground mine workings is encouraging and indicates that the model assumptions give a reasonable approximation of reality.

We observe a distinct seasonality for class I events, with the highest detection rates occurring in the winter months (December–February) as illustrated in Fig. 7. Detection rates are also high during the cold high-Arctic spring (March–May) and are lowest during the thaw season. We also observe a cluster of events south of SPITS, with activity highest in the summer months, decreasing in the fall (September–November), and absent during the winter and spring. These events may be rockfalls related to fluvial undercutting of steep river cliffs, rockfalls from steep mountain flanks, active-layer detachment slides/debris flows, or glacier/rock glacier movements. Using InSAR, Rouyet et al. (2019) measured high summer subsidence rates in the river valley, glaciers/rock glaciers, and mountain flank areas south of Janssonhaugen corresponding to the inferred source locations of these seismic events. However, since the dynamics of these processes are not represented by our model, we were careful to exclude these events when spatially isolating the cluster interpreted as class I events. By selecting the subset of events with inferred source positions within ∼ 1500 m of the array centroid, excluding the river valley/rock glacier area south of the array, we isolated a total of 42 432 class I events recorded between July 2004 and July 2021. Class I events occur most frequently during the cold winter and spring seasons, suggesting a relation to freezing processes.

Locally, the class I seismicity is dominated by three source areas (see Fig. 8) corresponding to areas with erosional scarps and solifluction lobes (Tolgensbakk et al., 2000). The erosional scarps, particularly the one on the northern side of the array (Fig. 7), in addition to the steep northwestern flank of Janssonhaugen, were also active during the summer thaw season. These thaw season events may be rockfalls or other mass movements, possibly initiated by the melting of ice causing loss of strength or joint lubrication (Matsuoka, 2019; Weber et al., 2017). Additional three-dimensional perspective views illustrating the three dominant source areas are provided in Appendix A.

Figure 9 Illustration of the spatiotemporal thermal-stress field modeled according to Eq. (12) and constrained by temperature measurements recorded by sensors installed at Janssonhaugen Vest (0.1 m) and in the P11 borehole (0.2–15 m). The temperature–depth profile used to calculate thermal stress was interpolated to regular 10 cm intervals using a spline interpolant. Sign convention is positive for tensional stress.

4.2 Modeled thermal stress

Figure 9 illustrates the spatiotemporal thermal-stress field that was modeled by solving Eq. (12) using the parameters listed in Table 1, the 0.1 m ground temperature time series recorded at the Janssonhaugen Vest meteorological station, and the 0.2–15 m temperatures recorded in the P11 borehole. The largest thermal stresses occur in the active layer, which is subject to large-amplitude winter cooling cycles. Since the peak annual stresses occur at 0.2 m depth and we have a much longer record from the P11 borehole compared to the Janssonhaugen Vest meteorological station, we focus on this horizon as the dominant thermal-contraction cracking depth throughout the rest of the study. That peak thermal stress modeled at 0.2 m is interesting, as it corresponds with the 20–30 cm thick regolith layer at Janssonhaugen (Isaksen et al., 2001), suggesting that cryoturbation within the active layer may have weathered the bedrock over time to produce the surficial layer. Aggradation of weathering material would be another explanation, but this is less likely, since the top of Janssonhaugen is a mountainous plateau where erosion is expected to dominate over deposition. Frost polygons in this region are interpreted to be very old (Sørbel and Tolgensbakk, 2002), so there has likely been sufficient time to reach a steady-state condition. If the highest thermal stresses occurred much deeper, we would expect that repeated thermal-contraction cracking over thousands of years would have produced a thicker regolith layer than is observed.

4.3 Thermal stress associated with seismic events

Figure 10 illustrates the shallowest measured temperature profiles and corresponding thermal stress for the September 2019 to June 2021 period. The class II events, corresponding to mining activities at Gruve 7, have source ranges of 6000–7000 m with a consistent azimuth of ∼ 210^∘ and no correspondence with thermal stress. On the other hand, class I events have variable azimuths and source ranges < 1500 m and tend to be associated with peaks in thermal stress.

Figure 10 (a) Ground temperature profiles (solid lines) and bearing (clockwise from north) to coherent-MFP-localized source positions from the center of the SPITS array (blue circles). (b) Modeled ground thermal stress from Eq. (12) (solid lines) and ranges to coherent-MFP-localized sources (blue circles). Source ranges and azimuths are transparent such that denser colors represent clusters of events.

Figure 11 shows a detailed comparison of the class I seismicity and the modeled frequency of tensile fracturing due to thermal stresses exceeding the assumed tensile strength of polycrystalline ice (as described in Sect. 3.3.1). Given that seismicity is a complex, stochastic process in time and space, our simple thermal-stress-based fracture model does a reasonably good job of capturing the time periods and approximate frequency of the class I events. This leads us to infer that thermal-contraction cracking of ice wedges and crack-filling vein ice, as modeled by Eq. (12), is a significant process contributing to event class I seismicity. The clusters of events recorded June–August, when thermal stress is low (see Figs. 7 and 11), are most likely mass movements associated with steep terrain (rockfalls, active-layer detachment slides, debris flows, etc.), possibly initiated by melting of fracture-filling ice leading to loss of strength or joint lubrication (Matsuoka, 2019; Weber et al., 2017).

Figure 11 (a) Pre-fracture thermal stress (black) given by the solution to Eq. (12) and post-fracture thermal stress accounting for stress release by cracking (red) with a blue color gradient indicating the 6-hourly binned detection rate of class I events. (b) Modeled number of frost quakes (red) according to the fracture model (see Sect. 3.3.1) and a 6-hourly binned frequency histogram of class I events (black).

A similar association between class I events and peaks in thermal stress persists over the entire study period, from 2004–2021 (see Fig. 12). Figure 13 shows that the modeled and observed frost quake seismicity also matches quite well over the study period (the normalized cross correlation is 0.61), though the observed seismicity has a tendency to occur in more defined periods than predicted by the model, resulting in a relatively spiky frequency histogram. Examples of anomalous seismicity not explained by the model are the periods 8–9 January 2008, 14–21 February 2010, 9–16 February 2012, and 5–12 January 2016 (see Fig. 13). In order to explain these anomalies, we note that the February 2010 episode corresponds temporally to the “C10” major cracking episode identified by Matsuoka et al. (2018) at a field site down-valley in Adventdalen. Matsuoka et al. (2018) observed that this cracking episode was preceded by a highly unusual period of mild weather, accompanied by rain, positive air temperatures, significant snowmelt, and surface water pooling. This surficial water subsequently froze when air temperatures dropped and an extensive series of fresh cracks were observed in the surficial ice by Matsuoka et al. (2018), during a field visit on 28 February 2010. This period of anomalous seismicity can therefore be explained as thermal-contraction cracking of the surficial ice formed after a period of heavy rain, and this explains why they were not accurately predicted by our subsurface thermal-stress model (since they occur in response to air temperature rather than ground temperature).

Figure 12 Modeled thermal stress from Eq. (12) at 0.2 m depth (black lines) and ranges to coherent-MFP-localized sources (blue circles). Source ranges and azimuths are plotted as transparent points such that denser colors represent clusters of events. Event classes I and II are associated with small (< 1500 m) and large (> 5000 m) source ranges, respectively.

Following this line of reasoning, we were able to connect (see Appendix B) all of the large, anomalous spikes in seismicity that were not predicted by our model with rare, heavy-rainfall events reported on by Dobler et al. (2019). These unusual winter rainfall events are driven by strong south to southwesterly atmospheric flows with advection of water vapor from warmer areas and are often linked to “atmospheric river” features in the precipitable water anomaly field (Serreze and Stroeve, 2015). That the anomalous spikes in seismicity consistently occur in the wake of significant mild weather and rain events (see Fig. B1) further strengthens the interpretation that thermal-contraction cracking of newly formed surface ice is a plausible explanation for the anomalous spikes in seismicity we observed at the SPITS array.

We consider the use of real ground temperature measurements a strength of this study, since the subsurface temperature field is a complex product of sensible, latent, and convective heat fluxes, depending on matrix and pore-filling material properties as well as surface properties like the thermal conductivity of snow (e.g., Badache et al., 2016; Rankinen et al., 2004). However, it is important to recognize that the P11 borehole and Janssonhaugen Vest ground temperature measurements are point samples of a temperature field that may vary spatially according to local geomorphology and variation in snow cover. For example, Abolt et al. (2018) have demonstrated that ground under the elevated rims of frost polygons cools significantly faster than the depressed centers due to decreased snow cover. Accounting for variation on the local scale fell outside the scope of the present study, but extending the model framework to allow for stochastic temperature variation via stochastic differential equations may give additional insights into the significance of local-scale variability.

Figure 13 Histogram binned to 9 d intervals comparing the modeled frost quake frequency (blue) with the recorded frequency of class I events (black). Anomalous spikes in seismicity following heavy-rainfall events (see Fig. B1) are annotated in red.

Figure 14 Mass movements such as (a) erosional scarps and (b) frozen-debris/solifluction flows initiate transverse cracks/fissures which become infiltrated by water that freezes to ice. (c) Thermal contraction of the surrounding ground and downslope gravitational stress causes the ice to crack under tension when the stress, σ, exceeds the tensile strength of ice, σ_T(ice).

4.4 Interpretation of event class I cryoseismic source

Polygonally patterned ground indicative of cryoturbation of the active layer and/or ice wedges in the underlying permafrost are observed extensively across Janssonhaugen, except where downslope mass movements destroy or interrupt the formation of polygonal networks (Sørbel and Tolgensbakk, 2002). In this study, we observed anomalously high cryoseismicity associated with erosional scarps and a frozen-debris/solifluction lobe, which are all associated with downslope mass wasting. As shown in Fig. 14, we suggest that these downslope mass movements initiate or precondition transverse cracks or fissures to open (Darrow et al., 2016; Price, 1974), particularly where there is a transition in slope, i.e., convex terrain. Water from rain or snowmelt infiltrates these fissures and freezes when temperatures drop below freezing (e.g., Darrow et al., 2016). Rapidly falling ground temperatures during winter then cause the surrounding ground to thermally contract. This causes the vein ice filling the frozen fissures to crack under tension, since the tensile strength of ice is lower than the surrounding ground and therefore constitutes a plane of weakness (e.g., Plug and Werner, 2002). Cracking relieves the accumulated thermal stress.

This mechanism is analogous to the Lachenbruch (1962) model of thermal-contraction cracking of ice wedges in permafrost (see Fig. 1). However, the cracking may occur more frequently or under milder surface cooling because (1) the frozen fissures may be pre-stressed by downslope gravitational forces making them more prone to failure and (2) the fissures extend to the ground surface where thermal stresses are largest. For the case of ice wedges, the ice wedge is located below the permafrost table (Fig. 1c). There may be vein ice extending through the active layer only if the previous season's thermal-contraction crack has remained open through the summer thaw season. If the crack has closed so that vein ice is not formed during the early freezing season, initiation of thermal-contraction cracking of the ice wedge would require accumulation of thermal-contraction stress at the level of the permafrost table, requiring a longer or more extreme surface cooling episode. Tensile cracks could also form in the active layer where ice veins/ice wedges are absent (as in Fig. 1a), but this would require thermal-contraction stresses exceeding the tensile strength of frozen ground, which is greater than that of polycrystalline ice (e.g., Plug and Werner, 2002) and may therefore occur less frequently.

5 Conclusions

We studied the spatial and temporal patterns of the class of seismicity associated with short-duration ground shaking at the SPITS array in Adventdalen, Svalbard, based on a catalog of > 100 000 events recorded between 2004 and 2021. To the best of our knowledge, this is a uniquely large and long-spanning event catalog amongst studies with a focus on cryoseisms. We find that these short-duration seismic events, with ground motion lasting just a few seconds, can be grouped into two main subclasses. One class is clearly associated with mining activities at the underground coal mine Gruve 7 and is mostly useful as an indicator of the performance of the coherent-MFP source localization algorithm. The other subclass of short-duration seismic events is interpreted to be dominated by frost quakes produced by thermal-contraction cracking of ice wedges and crack-filling vein ice. These events appear to be dominated by surface wave energy and source positions that are proximal to the SPITS array, particularly three areas that are associated with dynamic geomorphological features, erosional scarps and solifluction lobes. Temporally these events are associated with peaks in ground thermal stress as modeled by a simple dynamical thermo-viscoelastic model constrained by borehole measurements of ground temperature. The long-term continuous observational record, containing tens of thousands of inferred cryoseismic events, in close proximity to high-quality borehole temperature observations, provides a unique insight into the spatiotemporal patterns of cryoseismicity. Further experiments using controlled test sources in known locations could be used to calibrate MFP model assumptions that could subsequently benefit future MFP analyses of the complete, long-term seismic record from SPITS.

Appendix A: Three-dimensional perspective models of the study area

Three-dimensional perspective views of the study area were rendered by draping a composite orthophoto over a digital elevation model. Figure A1 shows the erosional scarp and frozen-debris/solifluction lobe on the southern flank of Janssonhaugen that were associated with spatial peaks in class I seismicity. Figure A2 illustrates the second erosional scarp on the northeastern flank of Janssonhaugen that was also associated with a peak in class I seismicity.

Figure A1 Three-dimensional perspective view of the southern flank of Janssonhaugen highlighting two of the three areas associated with anomalous seismicity. The three-dimensional model was constructed by draping an orthophoto on top of a 1.5× vertically exaggerated digital elevation model (DEM). Red crosses mark locations of SPITS seismometer stations. Orthophoto and DEM © Norwegian Polar Institute (https://npolar.no, 17 January 2022).

Figure A2 Three-dimensional perspective view of the northeastern flank of Janssonhaugen highlighting the third area associated with anomalous seismicity. The three-dimensional model was constructed by draping an orthophoto on top of a 1.5× vertically exaggerated DEM. Red crosses mark locations of SPITS seismometer stations. Orthophoto and DEM © Norwegian Polar Institute (https://npolar.no, 17 January 2022).

Appendix B: Meteorological conditions associated with anomalous seismicity

Figure B1 illustrates the meteorological conditions preceding anomalous spikes in seismicity that were not captured by the subsurface thermal stress and simple fracture model presented in this study. Meteorological conditions are represented by temperature and precipitation observations from Svalbard Airport (∼ 21 km from Janssonhaugen), since these variables were not recorded at Janssonhaugen over the relevant time period. The meteorological observations show that the anomalous periods of seismicity were preceded by > 0 ^∘C temperatures and rain in all instances. We expect that these rainfall events led to surface water pooling that subsequently froze to ice when temperatures again fell below 0 ^∘C. The periods of rapid cooling occurring synchronously with the anomalous seismicity would have then promoted intense thermal-contraction cracking of the exposed surficial ice.

Figure B1 Svalbard Airport air temperature (line) and daily precipitation (bar) records (http://met.no, 14 January 2022) illustrate that anomalous spikes in cracking-related seismicity were preceded by unseasonably mild temperatures > 0 ^∘C (red line) and rain, presumably leading to snowmelt, water pooling, and surface ice formation as observed in the field at Adventdalen by Matsuoka et al. (2018) in February 2010.

Appendix C: Additional details on signal characteristics of detected seismic events

Two classes of short-duration seismic events were identified. Event class I was interpreted to be dominated by cryoseisms, while event class II was associated with mining activities at the operational coal mine Gruve 7. Figure C1 shows that the dominant signal frequency is (on average) lower for class II than for the SPITS proximal class I events. We also see clustering at distinct dominant frequencies for event class II, probably related to differences in the specific mining activities that trigger these events. Event class I spans a wide range of dominant frequencies, though the most common is ∼ 17 Hz (see Fig. C1a). Bandwidth was estimated as the maximum signal bandwidth across all array stations (Fig. C1b–c). The bandwidth containing 99 % of the signal energy is quite high, since it is typically around 30 Hz, and the waveforms are pre-filtered with a 5–35 Hz passband. On the other hand, the half-power (3 dB) bandwidth is quite low, rarely exceeding 5 Hz. This matches our observation that even though class I events appear to be dominated by surface wave energy (which is typically dispersive), we do not observe the distinctive highly dispersive waveforms that one would expect for a wideband dispersive signal (see Fig. 5).

Figure C1 Illustration of the relationship between the dominant frequency estimated by Thomson's multitaper spectral amplitude peak (Hanssen, 1997) and the (a) MFP estimated range to source, (b) bandwidth containing 99 % of signal power estimated from the periodogram power spectral density, and (c) half-power bandwidth determined as the frequency range where power spectral density is within 3 dB of the maximum. The two event classes are most readily identified by the range to source, as annotated by the red dashed line.

Figure C2 shows illustrative examples of three component seismograms for the two identified event classes. There is some indication of phase rotation between the vertical and radial components of ground motion for the highest-amplitude arrival of event class I (frost quake), which would indicate a Rayleigh wave. The class II event (mining activity at the Gruve 7 coal mine) that is interpreted to be dominated by P-wave energy has vertical and radial components that are more similar to one another, as expected for body waves. In addition, the largest amplitudes for the Gruve 7 event correspond to the first arrival (consistent with interpretation as dominated by P waves), whereas the largest amplitudes for the frost quake arrive later (consistent with interpretation as Rayleigh wave dominated).

Figure C2 Three-component seismograms from SPITS station SPA0 for (a) an example of an event class I (frost quake) and (b) an example of an event class II (Gruve 7 mining activity). The signals are bandpass-filtered to 2.5–35 Hz, and the horizontal components of ground motion (measured in north–south and east–west orientations) were rotated to radial and transverse to the source bearings estimated by MFP (see Fig. 5e and f).

Appendix D: Frequency bandwidth and MFP

The MFP results in this study were found to be insensitive to the frequency band analyzed, particularly for the coherent-MFP formulation, as in Eq. (7). Figure D1 shows that as the passband is narrowed and resolution decreases somewhat for coherent MFP, but the estimated source position remains more or less constant. The incoherent MFP (classical formulation, Eq. 4) performs poorly throughout and is progressively more biased towards station number 8 as the passband is narrowed. We can reject this result because Fig. D1a clearly shows that the signal arrives at station 9 around the same time or slightly before station 8, so the incoherent-MFP maximum in the vicinity of station 8 (Fig. D1i) must be spurious and is likely a result of interfering sidelobes whose position varies with frequency.

Figure D1 Example of MFP source location for an event with a dominant frequency of 10.2 Hz and half-power bandwidth of 6.7 Hz (same event as shown in Fig. 5b) filtered to passbands of (a) 5–35 Hz, (b) 5–15 Hz, and (c) 7.5–12.5 Hz prior to MFP. Panels (d), (e), and (f) show corresponding coherent-MFP ambiguity surfaces, while panels (g), (h), and (i) show incoherent-MFP ambiguity surfaces.