Frost Quakes: Crack Formation by Thermal Stress

Fractures in frozen soils (frost quakes) can cause damage to buildings and other infrastructure, but their formation mechanisms remain poorly understood. A methodology was developed to assess thermal stress on soil due to changes in climate and weather conditions and to investigate the connection between thermal stress and frost quakes in central Finland due to brittle fracturing in uppermost soils. A hydrological model was used to simulate snow accumulation and melt, and a soil temperature model was used to simulate soil temperature at different depths beneath the snow pack. The results of modeling, together with measurements of air temperature, snow cover thickness, and soil temperature, were used to calculate temporal variations in thermal stress in soil. We show that frost quakes occur when thermal stress caused by a rapid decrease in temperature exceeds fracture toughness and strength of the soil‐ice mixture. We compared calculated thermal stress on soil, critical stress intensity factor, and a seismogram recorded in a suburban region in central Finland. Our results suggest that this methodology can be used to predict thermal stresses on soil and identify stress values that may lead to fractures of frozen soils, that is, frost quakes.


Introduction
The freezing of water in saturated soils or rocks, in ice, or on the boundary between ice and rock can cause sudden release of seismic energy in the form of nontectonic seismic events. The events originating from cracking of ice are usually referred to as ice quakes, while events originating in water-saturated soils are called frost quakes (Carmichael et al., 2012;Podolskiy et al., 2018;Zhang et al., 2019). Several studies showed correlation of ice quakes activity with rapid drops in air and surface temperatures. For example, MacAyeal et al. (2019) found hundreds of thousands of small ice quakes on McMurdo Ice Shelf, Antarctica, during the Austral summer melt season. Ice quake activity occurred over 6-12 hr time intervals in the evening when temperatures were falling. Carmichael et al. (2012) also observed seismic response from melting of glacial ice on Taylor Glacier, Antarctica. Other studies have found a correlation between cooling temperatures and the occurrence of ice quakes on frozen lakes in Alberta, Canada (Kavanaugh et al., 2019) and Laohugou Glacier in Tibet (Zhang et al., 2019). In an earlier study of ice quake activity in thin floating ice in Lake Suwa, Japan, Goto et al. (1980) found two clear peaks in seismicity due to ice quake activity: in the morning when the air temperature increased and in the evening when the air temperature decreased. Both periods of ice quake activity coincided with a nonuniform vertical temperature distribution in the ice (Goto et al., 1980). Cracking of ice or frozen soil tends to occur when thermal stress in the ice or frozen soil exceeds its tensile strength. The relationship between thermal shocks, cracking of ice, and tensile strength of ice has been studied to understand the mechanisms of ice quakes. Gold (1963) conducted a laboratory investigation by bringing together two ice blocks of different temperatures, calculating thermal stress and ice strength, and visually observing cracking of the ice blocks. Gold (1963) found that maximum strength of a smooth ice surface was 3 to 4 MPa. In the natural floating ice on Lake Suwa, Japan, Goto et al. (1980) calculated that significantly lower thermal stresses of 0.4 and 0.8 MPa for tensional and compressional stresses, respectively, led to ice quakes.
While there are systematic studies of factors affecting seismic activity of ice quakes, less is known about the factors that modulate frost quakes. A likely explanation for this is that ice quakes occur in locations that attracted particular attention of researchers in relation to climate change and its impact on the Earth's cryosphere (lakes, glaciers, and seas), resulting in a number of systematic instrumented studies of the connection between ice quakes activity and melting and freezing processes in these areas. Frost quake activity, on the other hand, occurs outside areas that are considered climate change indicators such as permafrost areas or areas covered by glaciers. Furthermore, frost quakes seem to occur at random, or less predictable, locations (UTSC Climate Lab, 2020). For this reason, sites of frost quake activity are rarely instrumented; thus, data of the ice conditions before, during, and after the frost quakes are typically not available. The occurrence and evidence of recent frost quakes is mainly based on reports from local citizens (Laine, 2016;Leung et al., 2017).
Frost quakes have been reported in North America (e.g., Battaglia & Changnon, 2016;Draxler, 2014;UTSC Climate Lab, 2020) and eastern Europe (Nikonov, 2010). Battaglia and Changnon (2016) and Nikonov (2010) proposed a mechanism for frost quakes in which liquid water trapped within pores of saturated soil freezes and expands due to a rapid decrease in atmospheric temperature followed by an abrupt decrease in soil temperature. This expansion of the freezing water increases the stress within in the soil, which can be released abruptly in an explosive way. Battaglia and Changnon (2016) found that frost quake episodes are typically preceded by days with temperatures greater than 0°C and occur when the surface and air temperatures drop rapidly by 15°C to 25°C during a period with some snow cover. Nikonov (2010) explained that frost cracking occurred in water-saturated soil when the temperature of the air and uppermost soil decreased by 10°C to 20°C.
The goal of this paper is to reconstruct the conditions that led to a swarm of frost quakes on 6 January 2016, in Talvikangas (Figure 1a), which is a suburban area in the city of Oulu in central Finland located at 65°N and belonging to so-called sub-Arctic or boreal regions. Some of the frost quakes created ruptures in soil, building foundations (Figure 1b), and roads ( Figure 1c). Local citizens reported hearing a loud noise and feeling the ground shake prior to observing the ruptures (Laine, 2016). On the same day, the permanent seismic station (OUL) operated by Sodankylä Geophysical Observatory of the University of Oulu located 14 km from Talvikangas recorded a number of unusual local seismic events, depleted in body wave energy but having large-amplitude Rayleigh waves (Figure 2, Afonin & Kozlovskaya, 2020). Large-amplitude Rayleigh waves suggest that the seismic source was close to the surface (Aki & Richards, 2002). Also on the same day, a rapid decrease in air temperature was recorded (FMI, 2020), suggesting that the frost quake activity was initiated by specific weather conditions that potentially can be repeated in the future.  (Laine, 2016). The damaged house is behind the people.
Because the frost quakes and subsequent damage occurred at an uninstrumented location, as is typical of most frost quakes, we used data from several nearby sources. Thermal stress conditions were inferred from observations of the damage caused by the frost quakes. Calculations of thermal stress are based on temperatures of air and soil. Since soil temperature and soil water content were not available at the Talvikangas site, we developed a hydrological model that can be used to calculate the snow accumulation and melt, the soil ice/water relationship in soils, and the soil temperature. First, we tested the hydrological model and calibrated the soil temperature model using data from the Karhinkangas study site, which is located in central Finland (see Figure 1a). In Karhinkangas, Hendriksson et al. (2018) measured hourly soil water content and soil temperature at different depths during the period of 2011-2015. Hendriksson et al. (2018) measured snow depth, snow water equivalent, and air temperature, and they collected precipitation and air temperature data from the automatic weather station hosted by the Finnish Meteorological Institute (FMI). See more information about the measurements and soil types in Talvikangas and Karhinkangas in section 2.3.
Next, we applied the same hydrological and soil temperature modeling approach to calculate thermal stresses of frozen soil at the Talvikangas site during January 2016, and we compared the results to the stresses estimated from seismic observations of the frost quakes. This approach can be used to predict thermal stress on soils due to changes in weather conditions (e.g., including those accompanying long-term climate variability and change) and subsequently to help communities to prepare for possible frost quake events that may cause damage to buildings and roads.

Seismic Observations
Based on reports and observations by citizens in Talvikangas, we can assume that the most likely fracture mechanism that resulted in the opening of a large fracture on the surface of the road (see Figure 1c) is vertical opening. As the event was recorded at only one seismic station, calculation of source parameters from seismograms was impossible. Moreover, conventional estimates of stress drop obtained for tectonic and mining-induced events in the same magnitude range are often based on Brune source models (Brune, 1970), which presume an idealized simple circular fault and hence are very uncertain. Nevertheless, direct measurement and estimations of crack geometry were conducted.
The width of the crack was evaluated by local citizens to be approximately 0.25 m. The crack length and depth were not measured because the crack occurred in the road, which was repaired quite soon after the event. However, based on observations by local citizens, the depth was between 0.5 and 1.5 m, and the length ranged between 30 and 100 m. Multiplying length, width, and depth of the fracture produces the seismic volume (e.g., the volume occupied by the source), which was between 7.5 and 75 m 3 based on the measured and inferred values stated above.

Theory
In this work, we evaluate the thermal stress in the soil in the vicinity of the Talvikangas frost quakes to determine whether the frost quakes were caused by thermal stress exceeding the strength of the soil. Two simplifying assumptions were made: (a) Thermal stresses and soil cracks occur only in fully frozen soil; and (b) horizontal and vertical displacement of frozen soil is not bounded.
Thermal stress, σ xx (z, t), in the subsurface at depth z from the land surface (z ¼ 0) and at time t can be calculated from the spatiotemporal temperature distribution in the soil as follows (Timoshenko & Goodier, 1951): is the soil temperature as a function of depth z and time t (°C), and d is thickness of the frozen soil layer (m). If the temperature distribution in the medium is known, thermal stress over the thickness of the body can be calculated using Equation 1. However, since frost quakes typically occur at uninstrumented sites, the soil temperature distribution is not known. Thus, we calculate soil temperature with the following approximation, based on the explicit finite difference approximation of the heat equation (Rankinen et al., 2004): where T(z, t + Δt) is soil temperature (°C) at depth z and time t + Δt, K T (Wm −1°C−1 ) is thermal conductivity of the soil, C A (J m −3°C−1 ) is apparent specific heat capacity of soil, which is the sum of specific heat capacity of soil, C s , and ice, C ice , (see details in Rankinen et al., 2004), T a (t) (°C) is the air temperature at time t, D s (m) is the snow depth, f s is an empirical snow parameter that includes the effect of thermal conductivity of snow, and Δt is the time step. This approximation assumes that the change in temperature below the calculation depth z is negligible (Rankinen et al., 2004).
For temperature at the soil surface beneath snow cover (i.e., at z ¼ 0), Equation 2 simplifies to 10.1029/2020JF005616

Journal of Geophysical Research: Earth Surface
Input parameters to Equations 2 and 3 are air temperature and snow depth. Since continuous snow depth measurements are rarely available, we used modeled snow depth, calculated as follows where SWE (m) is snow water equivalent and ρ snow is density of snow (kg m −3 ). Snow water equivalent was calculated as where the liquid water retention capacity of snow pack (WSP) and liquid water of ice in the snow pack (WIP) are obtained as follows: where P r (m) is rainfall, P s (m) is snowfall, M d (m) is daily snow melt, and F d (mm/day) is the freezing in snow pack. The hydrological model used to calculate SWE from 5-7 is explained in detailed in Vehviläinen (1992) and Okkonen and Kløve (2010). Daily snow melt M d (mm) and freezing in the snow pack F d (mm) are calculated as follows: where K m (mm°C −1 ) is the degree day factor. T a (°C) is a daily air temperature, T msnow (°C) is a temperature at which the snow begins melt, K f (mm°C −1 ) is a degree day freezing factor, T fsnow (°C) is a threshold temperature for freezing, and e (-) is a coefficient indicating the nonlinear relationship between refreezing and temperature (Bengtsson, 1982).
Once the soil temperature profile is obtained from Equation 2, it can also be used to calculate the soil ice content, I ce (-), given by Andersland and Ladanyi (2004): where T msoil (°C) and T fsoil (°C) are melting and freezing temperatures of soil water. Rembel (2008) showed that the fitting parameter β (-) can range between 0.19 and 1.15 with a mean of 0.42.

Study Sites and Data Availability
We collected soil samples from the Talvikangas site in the vicinity of frost quakes in 2019. In Talvikangas we excavated a soil pit to a depth of 80 cm, and we collected soil samples from depths of 30 and 50 cm for grain size distribution analysis. We collected only two soil samples because at the first look the soil appeared to be very homogenous, which was confirmed by sieve analysis. We conducted sieve analysis in the technical laboratory in the University of Oulu, Finland. At the depths of 30 and 50 cm, 97% and 91%, respectively, was retained on the sieves of diameters between 1 and 0.125 mm. Air temperature, T a (t) (FMI, 2020), precipitation, P r (t) (FMI, 2020), and snow water equivalent, SWE (SYKE, 2020) were observed at a nearby station (see inset in Figure 1a) and were available from the city of Oulu.

Journal of Geophysical Research: Earth Surface
OKKONEN ET AL.
In Karhinkangas, Hendriksson et al. (2018) excavated a soil pit to a depth of 200 cm and collected soil samples at the depths of 30 and 80 cm. They conducted sieve analysis, and 97% and 89% at the depths of 30 and 80 cm, respectively, was retained on the sieves of diameters between 1 and 0.125 mm. In the soils at both the Talvikangas and Karhinkangas sites, there is organic matter present at 0-10 cm depths, and the soil types are very similar.
Without soil temperature data, it is not possible to calibrate the soil temperature model described in section 2.2. Thus, in this study, the hydrological and soil temperature models were first tested in the Karhinkangas site in central Finland (Figure 1), 200 km south of Oulu, which has similar sandy soil as Talvikangas, and where long-term observations of soil temperature exist.

Karhinkangas Hydrological and Soil Temperature Model
In Karhinkangas soil temperatures at different depths were observed during the period of 2011-2015 (Hendriksson et al., 2018). Thermocouples were installed at depths of 5, 10, 30, 50, and 110 cm, and hourly soil temperatures were measured. Snow water equivalents and snow depth were measured three times in 2012. Daily air temperature and precipitation were obtained from the automatic weather station located 20 km from the Karhinkangas soil station (FMI, 2020). Density of the snow (not measured) used in the model was 225 kg m −3 (Vehviläinen, 1992).
First, we calibrated snow melt model against observed snow water equivalent. Table 1 presents the calibrated parameters for locations at Karhinkangas and Oulu. Simulated snow water equivalents (from Equations 5-9) and snow depths (from Equation 4) match well with observed values (see Figure 3 in Karhinkangas and Figure 5 in Oulu).
Next, the measured air temperature (Figure 4) and calculated snow depths from Figure 3 were used in the soil temperature model, Equation 2. The empirical snow parameter f s was fixed to −2.5, which was obtained from Rankinen et al. (2004). The soil model was calibrated by changing heat capacities of soil, C s , and ice, C ice , and thermal conductivity of the soil K T until simulated and observed soil temperatures matched well ( Figure 4). Table 2 presents the apparent heat capacity of soil, C A , and K T values at different depths. The Nash-Sutcliffe efficiency (NSE) was 0.92, 0.90, 0.91, 0.96, and 0.91 at depths of 5, 10, 30, 50, and 110 cm, respectively. We obtained slightly different parameters for different layers (Table 2). In Karhinkangas, thermal conductivity of soil increases with depth, likely due to organic matter content of the topsoil (0-10 cm)  which reduces the thermal conductivity relative to pure mineral soil. Rankinen et al. (2004) obtained values for C ice from 4 * 10 6 to 15 * 10 6 (J m −3°C−1 ) when C s is from 1 * 10 6 to 1.3 * 10 6 (J m −3°C−1 ). In our simulation, C s is 1 * 10 6 (J m −3°C−1 ) for every layer, and C ice is 5 * 10 6 (J m −3°C−1 ) for depths from 5 to 30 cm and 0 for depths from 50 to 110 cm, indicating no ice in these lower depths.

Talvikangas Soil Temperature Model With Karhinkangas Parameters
The input parameters for the hydrological model in Oulu/Talvikangas were daily temperature and precipitation, which were obtained from the weather station of Finnish Meteorological Institute (FMI, 2020) located about 10 km from the Talvikangas site.
The hydrological model presented in Equations 4-9 was calibrated against observed snow water equivalents (SWE) obtained from the Finnish Environmental Institute (2020) by adjusting degree day factor K m of snowmelt and degree day freezing factor K f . Details of model parameters can be found in Okkonen and Kløve (2010). The final values are shown in Table 1. SWE values were measured during the Period 2008-2013 about 4 km from the Talvikangas. Figure 5 shows the observed SWE values during the Period 2008-  The soil temperatures at different depths were calculated for Talvikangas using Equation 2 with the same parameters as for the Karhinkangas (Tables 2 and 3). In Talvikangas, there is a 10 cm-thick organic layer in the top of the soil column. Thus, to account for this organic layer, we used thermal conductivity of the soil K T ¼ 0.4 Wm −1°C−1 for 0-10 cm of depth and K T ¼ 1.4 Wm −1°C−1 for depths greater than 10 cm (Table 3). Figure 6 shows measured air temperature, simulated temperature beneath the snow cover (top soil), snow depth (from Equation 4), ice content of soil (from Equation 10), and simulated soil temperature (Equation 2) at the depths of 5 cm and 30 cm in January 2016.
Key features of the plots in Figure 6 are as follows. In January 2016, the mean air temperature in Oulu was −13.8°C, the minimum −30.1°C, and the maximum 4.1°C (Figure 6a). Simulated snow depth was largest, 9 cm, at the end of January 2016 (Figure 6b), during which time the top soil temperature reached its maximum of 3°C (Figure 6a). Soil temperature at the depth of 5 cm (blue line in Figure 6c) very closely follows the air/topsoil temperature (Figure 6a), while at greater depths, soil temperature responds to air temperature with a delay and the variation is smoother. At the depth of 5 cm, the water in the soil is frozen (ice content ¼ 1) between the days 2-22 January 2016 (blue line in Figure 6d); at the depth of 15 cm, soil was fully frozen on 4 and 6-8 January 2016 and during the night of 9 January 2016 (red line in Figure 6d), and at the depth of 30 cm, the soil is never frozen (the ice content is between 0.75 and 0.85 (-) and water content 0.25-0.15 (-) during all of January 2016). At the end of January, ice content is nearly uniform with depth, near 0.75, due to mild air temperatures during 23-31 January 2016.
Thermal stresses (Equation 1) of frozen soil were calculated for Talvikangas at depths of 5 cm (Figure 7d Table 3. We sampled the input variables from Equation 1 using the Matlab rand() command. The input parameters E and σ were randomly varied 10,000 times between 20 and 27 GPa and between 0.3 and 0.4 (Table 3), respectively.
At the depth of 5 cm, the maximum calculated thermal stress of 22 MPa (Figure 7d) was reached on 7 and 21 January 2016, while the mean stress was 8.5 MPa. Between 6 and 8 January 2016 and between 17 and 22 January 2016, thermal stress exceeded 12.9 MPa, which is the threshold of the associated thermal stress based on the critical stress intensity factor (K IC ). K IC ¼ 2.69 MPa m 1/2 is calculated with Equation 13 and its associated thermal stress, σ xx ¼ 12.9 MPa, with Equation 11. Equations 11 and 13 are explained in more detail in section 4 below. The frost quake recorded on 6 January 2016 ( Figure 2) occurred during the period of high calculated thermal stress. During the periods of high thermal stress, the air and soil surface temperatures were below −20°C. Thermal stress was also high on days 9, 11, 13, and 15 January 2016. Starting from 23 January 2016 the air and soil surface temperature increased and thermal stress decreased to below 6 MPa. At the depths of 15 and 30 cm, thermal stress (Figures 8d and 9d) was much lower than at 5 cm depth. The maximum calculated stresses were 1.0 and 0.5 MPa at 15 and 30 cm, respectively, which occurred on 7 January 2016.

Discussion
Prior studies have found that frost quakes occur when the air temperature drops abruptly over a short time (Battaglia & Changnon, 2016;Nikonov, 2010). A high negative rate of change of temperature (dT/dt) was found to correlate with high frost quake activity (Goto et al., 1980). On 6 January 2016, when frost quakes occurred in the city of Oulu, the cooling rate of air (and soil surface) temperature was −9°C day −1 (−0.38°C hr −1 ) (Figure 10a), while at the 5 cm depth, it was −0.6°C day −1 (−0.025°C hr −1 ) (Figure 10b), which were the highest (in absolute value) at these depths for the   month. These values are somewhat lower than what has been observed in other studies. For example, Kaminuma and Takahasi (1975) observed ice quakes when air temperature decreased at rates between 1 and 2.5°C hr −1 . Lombardi et al. (2019) found that ice quakes in the Antarctic ice sheet mostly occurred when the daily mean temperature varied more than 10°C. They also found that in colder weather, more ice quakes occurred with the rate of change 0.4-0.6°C hr −1 . Sugawara et al. (1982) observed cracking of asphalt by reducing temperature at rates between 3 and 30°C hr −1 . Our results are not fully comparable with asphalt or ice, but these studies give a range of thermal cycles for which cracks have occurred in different materials exposed to different changes in temperature. As stated above, Battaglia and Changnon (2016) and Nikonov (2010) found that frost quakes occurred under conditions of near freezing air temperatures that rapidly dropped by 10°C or more during a period with some snow cover. We found that the conditions at the Talvikangas site were consistent with these criteria during the occurrence of frost quakes: At the soil surface, a temperature change of 10°C (Figure 7a) was observed during the frost quake event, the soil was covered in snow (calculated snow depth was 2 cm, Figure 6b  linear elastic fracture mechanics (LEFM) can be used to describe the process of fracturing of a frozen soil layer in brittle or quasi-brittle conditions. According to LEFM (Irvin, 1957), fracturing and failure of elastic materials under Mode I loading, which corresponds to fracture opening, starts due to the presence of preexisting cracks or others impurities (Griffith, 1920). The condition for surface crack propagation and material failure is characterized by the value of the stress intensity factor, K I , that is related to the applied tensile stress σ xx (in our case, thermal stress) by the relationship (Tada et al., 2000; Figure 11): where η ¼ a/ω, a is initial vertical crack length that can be of the order of several cm, and ω corresponds to thickness of frozen uppermost layer of the soil estimated in the modeling ( Figure 11). The function f (η) is a geometric shape factor for unconstrained single edge crack plate subject to Mode I remote uniform tension and is given by Tada et al. (2000): f η ð Þ ¼ 1:12 − 0:23η þ 10:56η 2 − 21:74η 3 þ 30:42η 4 : According to Tada et al. (2000), the error of Equation 11 is less than 0.5% for b/ω ≥ 1 and a/ω ≤ 0.6, where b is length of the frozen soil layer. Because we use a 1-D temperature model, b goes to infinity ( Figure 11).
In order to calculate the value of the stress intensity factor K I , we assume that the stress σ xx represents thermal stress from Equation 1. In the region of Oulu, the observed frost depths in forested sandy soil and in  open land sandy soil have been, on average, 0.6 and 0.8 m, respectively, and in January 2016, they were 0.05 and 0.2 m, respectively (SYKE, 2020). In Talvikangas, the exact frost depth has not been measured; thus, for the thickness of frozen uppermost layer ω, we used values between 0.05 and 0.2 m. We sampled the input variables from Equation 11 using the Matlab rand() command. The input parameters, a and ω, were randomly varied 10,000 times between 0.01 and 0.1 m and between 0.05 and 0.2 m, respectively. The value of K I varied from 1.6 to 42.9 MPa m 1/2 from 5 to 7 January 2016, respectively. On 6 January 2016, the day of the frost quake, the minimum value of K I was 2.7 MPa m 1/2 , and the maximum 32.5 MPa m 1/2 Figure 12.
In order to demonstrate that the thermal stress estimated in this study is large enough to initiate brittle fracturing in the uppermost frozen soil (5 cm), we compare the maximum value of K I from 6 January 2016 32.5 MPa m 1/2 to the critical values of K IC (also called fracture roughness). The crack will propagate when the K I ≥ K IC, thus providing a simple failure criterion for brittle material.
As shown by Zhang (2002), the K IC and the tensile strength of rock are empirically related for Mode I fracture; hence, if tensile strength of the material is known from measurements, the K IC can be calculated as where σ t is tensile strength has units of MPa m 1/2 . We obtained values for K IC (Equation 13) from 0.15 to 2.6 MPa m 1/2 with tensile strengths from 1 to 22 MPa (Zhou et al., 2015), respectively, and associated thermal stresses, σ xx , (Equation 11) were 0.17 and 12.9 MPa, respectively. The critical values of K IC calculated from tensile strengths for wide range of rocks vary from 0.2 to 3.0 MPa m 1/2 (Zhang, 2002). Li and Yang (2000) reported that fracture toughness of frozen sand is between 0.25 and 0.7 MPa m 1/2 within the temperature range from −4°C to −10°C. The values of K IC for ice in grounded glaciers with no-slip condition at the base were estimated with modeling by Jimenez and Duddu (2018), who obtained values from 0 to 10 MPa m 1/2 . Our calculations thus suggests that at the depth of 5 cm on 6 January 2016, thermal stress ( Figure 8d) due to intrusion of cool air of insulating snow cover exceeded the K IC ¼ 2.6 MPa m 1/2 which therefore was the probable cause of crack formation and the frost quake swarm.
Calculated thermal stress was also higher than K IC during the Period 17-22 January 2016 (Figure 8d), but people in Talvikangas did not reported soil cracks during those days. On 6 January 2016 at the   (Figures 8d and 9d), and soil temperature gradient (Figures 10c and 10d) were higher and the soil temperature lower (Figures 8c and 9c) on 6 January 2016 than during the Period 17-22 January 2016. The thickness of the frozen layer may have impacted the frost quake event. Lombardi et al. (2019) proposed that the brittle failure and stress release in the case of a thicker frozen layer may produce larger-amplitude ground motion. Hence, it is possible that cracking in the thicker frozen layer on 6 January 2016 produced larger magnitude frost quakes, and consequently more observable damage, than cracking in the thinner frozen layer on 17-22 January 2016. The recent analysis of the seismic signal from the station OUL (Afonin & Kozlovskaya, 2020) suggests that frost quakes with similar seismic signal characteristics occurred during the Period 17-22 January 2016 but were not as large in the magnitude as the event on 6 January 2016.

Conclusions
Frost quake crack formation by a large temperature drop in air and soil is a natural phenomenon that is not well understood. Although formation of a single fracture cannot result in a seismic event comparable in magnitude with tectonic earthquakes, the ground shaking produced by frost quakes can be hazardous in urban areas, as the example of Talvikangas frost quake shows. To understand conditions for frost quake formation, we developed a simple approach to quantify stresses on soil due to changes in air temperature. With the use of hydrological and soil temperature models, we estimated thermal stresses on soil. Our calculations showed that calculated thermal stresses on soil due to rapid decrease of temperature and absence of thick insulating snow cover were greater than the stresses that are needed to initiate cracking of ice and frozen soil. The air temperature dropped by as much as 9°C in the day of the frost quake (6 January 2016), and the calculated thermal stresses on the soil at a depth of 5 cm were up to 22 MPa, which was above the threshold thermal stress calculated with the critical stress intensity factor, K IC . As climate studies show, the present-day sub-Arctic (boreal) weather conditions are generally characterized by temperature and snow cover variability. If such variability were to increase in the future, this would also result in favorable conditions for formation of frost quakes. Our approach for calculating thermal stresses on soils could help communities to anticipate frost quake events that may damage buildings and roads.
Between 17 and 22 January 2016 thermal stress in the upper soil layer was also high due to cool air, but during that time no cracking was observed by the local people in the Talvikangas area. We suggest that initiation of the process of cracking could be related to the growth of the thickness of the frozen layer, in addition to soil temperature and its time derivative. During the frost quake on 6 January 2016, the simulated ice content at depths of 15 and 30 cm was greater, soil temperature lower, and its time derivative higher (in absolute value) than during the period of 17-22 January 2016. For future studies, we suggest that more detailed information of critical stress intensity factor (K IC ) of frozen soil is needed along with the parameters such as Young modulus and Poisson's ratios of frozen soils together with thermal experiments and cracking tests in order to verify the modeling results.