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 Artic 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 catalogue 10 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 catalogue 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 thermal contraction cracking of ice wedges or other segregated ice bodies. This inference is supported by the correspondence between peaks in observed seismicity with peaks in modelled ground 15 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; boulder producingerosional scarps and a frozen debris/solifluction lobes. Seismic stations providing year-round, high temporal resolution measurements of ground motion may be highly complementary to satellite remote sensing methods, such as InSAR, for studying the dynamics of periglacial 20 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 temperaturess. test contraction transient acceleration spatial and temporal SPITS than a longer and nearly continuous record is available. This allows a more rigorous investigation of the temporal correlation of these events with ground cooling and thermal stress accumulation.


1
Introduction 25 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 pressure stress release (Barosh, 2000;Battaglia et al., 2016). Since this pressure 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 30 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 Cracking typically occurs in response to 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., 35 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 heaveFrost damage (e.g. Rempel, 2010) can be understood as a combination of slow creep (frost heave, e.g., Rempel (2010)) and rapid elastic (frost quakes) 40 deformation of frozen ground and causes damage to roads requiring billions of dollars annually to repair in the United States alone (DiMillio, 1999). 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 days 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 times of yearseasons. As a result, InSAR is more suited to studying thaw subsidence than frost heave in areas such aslike Svalbard that 95 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.

100
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 sporadic 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 hadwith durations of just a few seconds (indicative of a nearby source) and peak amplitudes significantly above the background noise level. For comparison, regional tectonic 105 earthquakes are typically associated with >30 s ground shaking duration since the different velocities of P, S and surface waves causes the wavefield 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 nearregional 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, 110 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 115 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 wavefield sampling of the SPITS array is much coarser than the temporary geophone array they deployed, 120 a much longer and nearly continuous record is available. This allows a more rigorous investigation of the temporal correlation of these events with ground cooling and thermal stress accumulation. Study area and data The small-aperture Spitsbergen seismic array (SPITS) is located on Janssonhaugen, in the Adventdalen valley on the island of Spitsbergen, part of the high-Arctic Svalbard archipelago (Figure 2Figure 1). The SPITS array has been in operation since 125 1992, maintained by the research foundation NORSAR (Schweitzer et al., 2021). At present, it consists of 9 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 -3dB) 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 130 seismometers (Schweitzer et al., 2021). The waveform data following the upgrade is of high quality and well suited to source localization using broadband, coherent 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. Janssonhaugen is a bedrock remnant located in the middle of the Adventdalen valley, with a ~0.2-0.3 m thick weathered 140 sediment crust overlying homogeneous sand-/siltstone bedrock (density 2280kg.m -3 , porosity 20-25 %, >96.5 % SiO2) corresponding to the Ullaberget Member of 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 145 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 formationindicating the presence of sand/ice-wedges. The homogeneous geology at Janssonhaugen, it's flat topography, limited snow cover, relatively large distance to glaciers, rivers, ocean, 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). 150 Figure 32 -Illustration of spatiotemporal borehole temperature history recorded at the PACE 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 seasonal freeze-thaw. The continuous timeline is split across multiple figure panels. Temperature-depth profile was interpolated to regular 10 cm intervals 155 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-Oct to 28-Feb 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 160 surface energy budget (Westermann et al., 2009). As a result of the complex interplay of processes driving temperature variation in this high-Arctic location, Ttemperature 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 hours extending from April 1999 to the present (Isaksen et al., 2001). We 165 focus on the P11 borehole (see Figure 2Figure 1), which that 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 Figure 3Figure 2. 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 inter-annual variability in the magnitude of summer warming and winter cooling. The upper part of the permafrost at the P11 borehole 170 location is sampled by thermistors installed at 2.5, 3, 3.5, 4, 5, 7, 10, 13 and 15 m. A long-term warming trend is observed in the permafrost beneath the active layer (Figure 3Figure 2). Furthermore, the Janssonhaugen Vest weather station (see Figure   2Figure 1), which was installed in September 2019, provides includesan hourly sampled records of air temperature and ground temperature at 0.1 m depth. It therefore and provides a basis to compare depth and temporal sampling effects against the longer duration, more coarsely sampled P11 borehole record. 175

STA/LTA detection of short duration seismic events
Events are detected based on anomalous values of short-time-averaged (STA) amplitude 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, first-pass automatic identification of short duration seismic signals, 180 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 was first de-trended, 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 Figure 4Figure 3a). This passband was selected to eliminate high-frequency random noise which can lead to spurious spikes in STA/LTA. 185 (e.g. Allen, 1982;Trnkoczy, 2009)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 one-second moving-average smoothed trace envelope for each seismogram. The LTA is the STA further smoothed according to a 20 second period moving average. 190 Events are detected based on anomalous values of short-time-averaged amplitude divided by long-time-averaged amplitude, i.e., the classic STA/LTA approach (e.g. Trnkoczy, 2009). In our implementation, the short-term average (STA) from the trace envelope for each seismogram, given by the magnitude of the Hilbert transform and smoothed by taking the one-second moving average. Since we have an array of stations, we represent the array-STA by taking the 80 th percentile station-STA across all 195 stations at a given time sample. By visual inspection of test periods, we found that this emphasizes very local events with large amplitude variation across the array, while still ensuring that there is at least some coherency across the array. If we had chosen the maximum array station STA, we may detect arbitrary noise spikes with large amplitudes registered on a single seismometer.
If the mean or median array 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 200 seismometers.
The STA/LTA threshold was set to 10. Furthermore, after a trigger, no new events are declared within 5 seconds 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, trigger pause) were found to be appropriate for detection of short-duration 205 signals coherent across the array, while avoiding false triggers (noise bursts at single stations), by visually inspecting test periods. We calculate the long-term average (LTA) by simply smoothing the short-time average according to its moving average over a time-span of 20 seconds. This time-span is significantly longer than the events we seek to detect but shorter than the typical duration of the more regional scale seismic events that we want to ignore. We further reject epochs with LTA more than 2.5x 5 times the 2-hourly mean of the LTA in order to filter out large regional eventsearthquakes. Events are then 210 detected by applying a peak finding routine to the STA/LTA ratio. The STA/LTA peaks must have amplitude ≥10 and occur at least 5 seconds apart from one another. An example of event detection based on STA/LTA peaks is shown in Figure 4Figure 3b, which illustrates that short duration events are selected while a longer duration high-amplitude event is not detected as intended.

Source localization by coherent MFP 220
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 wavefield recorded at a receiver array and a series of predicted wavefields 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 225 of wave propagation within the medium. We assume a simple homogeneous medium with amplitude decay according to spherical divergence where the replica column vector, ′ is represented by the theoretical harmonic wave emitted from a test point ( , ) according to Here, = √ −1, is the angular frequency, = [ 1 , 2 , ⋯ , ] is an × 1 vector containing the absolute Euclidean 230 distances between ( , ) and the recording stations where superscript denotes a transpose, ( ) is the medium phase velocity at frequency and the vector is normalised to unit length by dividing by the factor = √∑ (1 2 ⁄ )

=1
. The column data vector is given by where ( ) is a frequency transform of the j-th trace ( ) of N traces recording a specific seismic event. We then form the 235 complex-valued × cross-spectral density matrix (CSDM) by where (•) 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. (wIn this study wee selected the 5-35 Hz band), sampled at an interval of 1 Hz, in order to reject spatially coherent low-frequency 240 background noise while retaining the shared signal/noise high-frequency band with a specific frequency sampling interval (1 Hz) (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, ( , ). Conventionally, the MFP coherence is estimated using the linear Bartlett processor 245 which evaluates the inner product between the recorded and predicted wavefields 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 cross-frequency spatial coherence structures to be exploited to give improved robustness and 250 accuracy. In this scheme, measurement vectors at discrete frequencies are concatenated to form the × 1 measurement super-vector where = [ 1 , 2 , ⋯ , ] is a frequency vector, and the replica vectors are concatenated to form the × 1 replica supervector 255 The super-CSDM ( ) = ( ) ( ) is then composed of × elements and the generalized MFP coherence is given by = ′ ( , ) ( ) ′ ( , ).
The estimated source position, ( , ), and phase velocity, ( ), are those which maximize the coherence measure . Since 260 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.

Ground thermal stress model
Previous studies such as Lachenbruch (1962), Mellon (1997), Maloof et al. (2002), Schulson and Duval (2009)  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 become dominant. Following Mellon (1997), we model the frozen soil as a Maxwellian viscoelastic solid augmented with thermal expansion and contraction. This allows us to decompose the total strain tensor into three components: an elastic ( ), a thermal ( ), and a viscous ( ) component, where subscripts indicate tensor components, , = 1,2,3, where 1 and 2 denote horizontal components, and 3 is the vertical component. The elastic strain tensor is related to the stress tensor by (e.g., Landau and Lifshitz, 1970 where is Poisson's ratio, is Young's modulus, is the Kronecker delta, and Einstein's summation convention is applied throughout this paper. 275 The thermal strain tensor is a measure of the change in volume caused by the thermally driven deformation, and it is expresse d as (e.g., Landau and Lifshitz, 1970) where is the temperature, 0 is a reference temperature for the undeformed state, and is the linear thermal expansion 280 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 where Γ {•} is a nonlinear operator acting on the deviatoric stress = − /3. The chosen parametric form of Γ {•} 285 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 = 11 = 22 and 11 ⁄ = 22 ⁄ = 0. The vertical stress is also neglected, i.e., 33 ≈ 0, corresponding to the assumption that the overburden is negligible in the shallow subsurface (Mellon, 1997). We 290 assume that Poisson's ratio, = ( ), Young's modulus, = ( ), and the coefficient of linear thermal expansion, = ( ) are all temperature dependent, and thus, implicitly time dependent since = ( , ) is the temperature at depth and time . By direct evaluation of ⁄ and collecting terms, we find that the temporal dynamics of horizontal stress ( , ) in a Maxwellian viscoelastic solid driven by thermal expansion and contraction is governed by the following first-order nonlinear and nonhomogeneous differential equation 295 The time-dependent coefficients are found to be and The scientific literature devoted to the rheology of frozen materials favor power-law parametrizations 300 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): In Glen's flow law, R is the universal gas constant, and 0 , , and 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 305 as the 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 ( , ), we specify the initial condition 0 ( ) = ( , = 0) = 0.
If we assume = 0 ⁄ , Eq. (12) reduces to the model proposed by Mellon (1997). By contrast, if we assume 310 , Eq. (12) reduces to the model proposed by Podolskiy et al. (2019). Finally, if we assume ⁄ = ⁄ = ⁄ = 0 = 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.

315
We solved Eq. (12) numerically to obtain a time series of the resulting horizontal stress ( , ), at depth , 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 ( , ) at depth . 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, ( , ), at Janssonhaugen has been 320 logged by a series of thermistors installed in the 15 m deep P11 borehole at 6 hr intervals since April 199 9. We assume the stress at a given depth is decoupled from the stress at adjacent depths (Mellon, 1997) nd solve Eq. (12) for the temperature timeseries at the selected depth of investigation. The set of physical parameters that we assume is described in Table 1. Butkovich (1959) gives 3 rd order polynomial for ice and a 2 nd order polynomial was fitted to Kell (1967) water thermal expansion data Viscous pre-factor 0 , −1 − 1 × 10 −9 Glen's non-Newtonian power law viscous flow (Behn et al., 2021;Glen, 1955;Weertman, 1983 Currier and Schulson (1982), varies according to grain size for randomly oriented polycrystalline ice (finer grained ice stronger), insensitive to temperature (Petrovic, 2003)..

Fracture model
We apply a simplistic model of fracturing by considering the modelled thermal stress as the potential pre-fracture stress. When the potential pre-fracture stress exceeds the tensile strength of the ground (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 330 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. Our aim is to model simple thermal tensional cracking of existing ice wedges or crack filling vein-ice segregated ice bodies and not the initiation or propagation of new cracks into previously undamaged soil/bedrock, where pore scale fluid 335 migration and stress localisation at crack tips become important (e.g. Walder and Hallet, 1985).

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 localisation 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 340 were subsequently used to identify subclasses of events, as detailed in the following section.

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 characterised 345 by significant amplitude variation and arrival time differences across the array seismometers, as illustrated in Figure 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 radius of ~1500 m centred over the array (Figure 5d/e). By contrast, Event class II is characterised by similar amplitudes across the array elements and smaller relative arrival time differences, as illustrated in Figure 5c. Using coherent MFP, we find that these events are associated with more distal inferred source positions (Figure 5f) that also have a 350 consistent azimuth. Figure 5 also illustrates the property that coherent MFP decreases source localisation 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 gives good constraint on source range under the assumed model of amplitude attenuation due to geometrical spreading (Figure 5f). In Figure 5i, we see that incoherent averaging (Eq. 4) has enhanced a sidelobe and produced an incorrect 355 source position that is not consistent with the relative arrival times observed in Figure 5c (earliest arrival at station 7 and latest arrivals with weakest amplitudes at stations 5 and 9). The mean MFP inferred propagation velocities 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 5 Hz bandwidth (see Appendix C), so the distinctive highly dispersive waveforms that one would expect for a wideband dispersive signal were not observed ( Figure 5). By contrast, the mean MFP inferred 370 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 Class I and II are dominated by surface waves and P-waves respectively (see Appendix C).

Results and discussion 375
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 localisation by coherent MFP, only events recorded by at least five seismomet ers 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. Mapping the inferred source positions, as in Figure 6, it is clear that Class II events are spatially coincident with the 385 underground mining areas of Gruve 7. Furthermore, we observe that Class II events occur frequently during all seasons (see Figure 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 localisation performance of the coherent MFP algorithm as applied in this study. The accuracy of the inferred source position s 390 is largely a function of data versus 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 395 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.

Subclasses of short duration seismic events 400
We find that there are two main sub-classes of short duration events recorded by the SPITS array. Event class I is characterised by significant amplitude variation and arrival time differences across the array seismometers, as illustrated in Figure 4a,b.
Using coherent MFP to infer the source positions of these events, shows that they occur in relatively close proximity, within about 1500 m of the centroid of the array (Figure 4d,e). By contrast, Event class II is characterised by similar amplitudes across the array elements and smaller relative arrival time differences, as illustrated in (Figure 4c). Using coherent MFP, we find not 405 only are these events associated with more distant sources (Figure 4f), they also have a consistent azimuth. Figure 4 also illustrates the property that coherent MFP decreases source localisation ambiguity for arrays that coarsely sample the spatial domain when compared to the incoherent scheme (consistent with Michalopoulou, 1998). In Figure 4i, 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 Figure 4c. 410 The mean MFP inferred propagation velocities for Class I events was 1150 m/s with a standard deviation of 1100 m/s, indicating that they are dominated by surface waves. The large standard deviation may indicate the surface waves are dispersive with different frequencies propagating at different phase velocities. 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 415 energy.

420
Mapping the inferred source positions, as in Figure 5, it is clear that Class II events are spatially coincident with the underground mining areas of Gruve 7. Furthermore, we observe that Class II events occur frequently during all seasons (see Figure 5). We infer that Event class II is a result of mining operations and human activity in the underground coal mine, Gruve 7. These events are essentially unwanted noise from the environmental seismology perspective, although they do give a useful indication of the localisation performance of the coherent MFP algorithm as applied in this study. 425

430
We observe a distinct seasonality for Class I events, with the highest detection rates occurring in the winter months (Dec-Feb) as illustrated in Figure 7Figure 6. Detection rates are also high during the cold high-Arctic spring (Mar-May) and are lowest during the thaw seasonsummer (Jun-August). We also observe a cluster of events due south of SPITS, with highest activity in the summer months, decreasing in the autumn (Sep-Nov) 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 435 slides/debris flows or glacier/rock glacier movements. Using InSAR, Rouyet et al. (2019) measured high summer subsidence rates in the river valley, glacier/rock glacier 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 440 area south of the array, we isolated a total of 42,432 Class I events recorded between July 2004 and July 2021By selecting the subset of events with inferred source positions within ~1500 m of the array centroid and excluding the cluster of summerautumn events due 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.

445
Locally, the Class I seismicity is dominated by three source areas (see Figure 8Figure 7) corresponding to areas with erosionalboulder producing scarps and solifluction lobes (Tolgensbakk et al., 2000). These areas may be associated with enhanced ground heat loss, thin or absent snow cover or elevated ground moisture/ice content (e.g. Abolt et al., 2018;Matsuoka, 2008), though we lack the field observations necessary to support this explanation for the anomalous seismicity of these areas. The boulder producingerosional scarps, particularly the one on the northern side of the array (Figure 7Figure 6), 450 in addition to the steep NW flank of Janssonhaugen, were also active during the summer thaw season. These thaw season events may be rockfalls or other mass movements, possibly initiated bydue to 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.  (0.2-15 m). The temperaturedepth profile used to calculate thermal stress was interpolated to regular 10 cm intervals using a spline interpolant.

Modelled thermal stress
Figure 9Figure 8 illustrates the spatiotemporal thermal stress field that was modelled by solving Eq. (12) using the parameters 460 listed in Table 1, and a combination of the 0.1 m ground temperature timeseries 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 frost thermal contraction cracking depth throughout the rest of the study. That peak thermal 465 stress is modelled at 0.2 m is interesting as it corresponds with the 20-30cm 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. suggesting that in-situ frost cracking may have weathered the bedrock 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 470 be very old (Sørbel and Tolgensbakk, 2002) so o there has likely been sufficient time to reach a steady state condition. If the highest thermal stresses occurred frost cracking extended much deeper, we would expect that repeated thermal contraction cracking over thousands of years would have produced a thicker regolith layer than that which is observed.

4.3
Thermal stress associated with seismic events to Jun 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, but source ranges <1500 m and tend to be associated with peaks in thermal stress.

Figure 109 -(a) Ground temperature profiles (solid lines) and bearing (clockwise from North) to coherent MFP localised source positions from centre of SPITS array (blue circles). (b) modelled ground thermal stress from Eq. (12) (solid lines) and ranges to coherent MFP localised sources (blue circles). Source ranges and azimuths are transparent such that denser colours represent clusters of events.
Figure 11Figure 10 shows a detailed comparison of the Class I seismicity and the modelled frequency of tensile fracturing due 485 to thermal stresses exceeding the assumed tensile strength of the mediumpolycrystalline ice (as described in Section 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-icethermal contraction cracking of segregated ice bodies, as modelled by Eq. (12) is a significant process contributing to event Class I seismicity. The clusters of events recorded June-490 August, when thermal stress is low (see Figure 7Figure 6 and Figure 11Figure 10), are most likely mass movements associated with steep terrain (e.g., rockfalls, active layer detachment slides, debris flows etc.), possibly on the steep boulder producing scarps (see Figure 7) i initiated by melting of fracture-filling ice leading to loss of strength or joint lubrication (Matsuoka, 2019;Weber et al., 2017).

500
A similar association between Class I events and peaks in thermal stress persists over the entire study period, from 2004-2021 (see Figure 12Figure 11). Figure 13Figure 12 shows that the modelled and observed frost quake seismicity also matches quite well over the study period (the normalised 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 Jan 2008, 1714-261 Feb 2010, 79-16 Feb, 2012 172 Jan 2016 (see Figure 13Figure 12). In order to explain these anomalies, we note that the Feb 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 510 Matsuoka et al., during a field visit on 28 Feb 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). Dobler et al. (2019); (Serreze and Stroeve, 2015)We consider the use of real ground temperature measurements a strength of 515 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 520 polygons cools significantly faster than the depressed centres 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 stochastic temperature variation via stochastic differential equations may give additional insight into the significance of local scale variability.
While these periods of enhanced seismicity are still associated with peaks in thermal stress, our simple model is unable to 525 explain the intensity of seismicity in these periods relative to, e.g., earlier periods in the same seasons with similar thermal stress but fewer Class I events (see Figure 11). There are likely multiple factors for this deviation, including the inherently stochastic nature of seismicity. However, there are likely also systematic factors that are simply not accounted for in our model like the spatial redistribution of stress between multiple fractures and across multiple depths. In addition, the mild temperatures preceding these episodes of enhanced seismicity (seen as pronounced drops in stress in Figure 11) may have produced unfrozen 530 water and conditions favourable for the formation of new ice lenses (Hallet et al., 1991;Murton et al., 2006;Peppin and Style, 2013;Walder and Hallet, 1985). That is, there are other processes that can drive frost cracking in addition to the thermal contraction cracking of segregated ice bodies and ice wedges that is favoured by our model. We speculate that the inability of the thermal contraction model to explain the anomalously high levels of seismicity in the periods 17-26 Feb 2010, 7-16 Feb, 2012 and8-17 Jan 2016 provides evidence that capillary and frozen-fringe effects are also important (Peppin and Style, 2013). 535 Following this line of reasoning, we were able to connect (see Appendix B) all of the large, anomalous spikes in seismicity 540 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 south-westerly 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 Figure  B1) further strengthens the interpretation that thermal contraction cracking of newly formed surface ice is a plausible 545 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, 550 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 centres 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 stochastic temperature variation via stochastic differenti al equations may 555 give additional insight into the significance of local scale variability.
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 560 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 centres 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 stochastic temperature variation via stochastic differential equations may give additional insight into the significance of local scale variability. 565 Figure 1312 -Histogram binned to 9-day intervals comparing the modelled frost quake frequency (blue) with the recorded frequency of Class I events (black). Anomalous spikes in seismicity following heavy rainfall events (see Figure B1) are annotated in red.

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 575 in Figure 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 580 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 Figure 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 590 ground surface where thermal stresses are largest. For the case of ice wedges, the ice wedge is located below the permafrost table (Figure 1c). There may be vein ice extending through the active layer only if the previous seasons 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 595 cracks could also form in the active layer where ice veins/ice wedges are absent (as in Figure 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
Conclusion 600 We studied the spatial and temporal patterns of the class of seismicity associated withthat lead to short duration ground shaking at the SPITS array in Adventdalen, Svalbard, based on a catalogue of >100 000 events recorded between 2004 and 2021. To the best of our knowledge, this is a uniquely large and long spanning event catalogue 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 605 mostly useful as an indicator of the performance of the coherent-MFP source localisation 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-icethermal contraction cracking of ice wedges or other sub-surface segregated ice bodies. 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; boulder producingerosional scarps and solifluction 610 lobes. Temporally these events are associated with peaks in ground thermal stress as modelled by a simple dynamical thermoviscoelastic 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 615 analyses of the complete, long-term seismic record from SPITS.

Appendix A. 3D perspective models of 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 620 on the north-eastern flank of Janssonhaugen that was also associated with a peak in Class I seismicity.

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 635 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.

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 645 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 Figure C1a). Bandwidth was estimated as the maximum signal bandwidth across all array stations ( Figure C1-b-c). The bandwidth 650 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 don't observe the distinctive highly dispersive waveforms that one would expect for a wideband dispersive signal (see 660 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 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 665 with interpretation as P-wave dominated), whereas the largest amplitudes for the frost quake arrive later (consistent with interpretation as Rayleigh wave dominated). poorer throughout and is progressively more biased towards station number 8 as the passband is narrowed. We can reject this result because Figure D1-a 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 ( Figure D1-i) must be spurious and is likely a result of interfering sidelobes whose position varies with frequency. 680 The study was conceptualised by AK and RR. The model, theory and methodology were developed by RR and AH. RR carried 880 out the data collation and processing and was responsible for analysing and visualising the data. RR drafted the initial manuscript with contributions from all authors, with AH contributing significantly to the development and drafting of the theory sections. AH and AK also provided project supervision.