Master thesis, 60 hp
Master thesis in Earth Science, 60 credits Vt 2021
Identifying subarctic river thermal and mechanical ice break-up using
seismic sensing
Ștefania Ursică
Supervisor: Lina Polvi Sjöberg
Abstract
River-ice break-up in high-latitude regions, despite its brevity, is a fundamental process, representing the most dynamic and complex period of fluvial processes. Moreover, ice break- up has significant cascading ecological effects, with a different severity for mechanical vs.
thermal break-up, and thus, motivates the importance of monitoring efforts. Classical research methods, such as fieldwork or analysis of photographs and aerial imagery, offer a general perspective on the timing of ice break-up but have safety and logistic issues caused by the dangers of unstable ice cover, the lag times between event occurrence and observation, and the frequent low visibilities. The emerging field of environmental seismology, which studies surface processes through seismic signals, provides an alternative solution to these shortcomings by continuously recording high temporal resolution data. Seismic sensing can potentially record any event within a set distance if the produced signal is powerful enough.
Three geophones had monitored the subarctic Sävar River reach for 185 days to test the efficiency of seismic methods to capture ice-cracking events, and based on their characteristics, to identify thermal vs. mechanical ice break-up. With visual and multivariate analysis, seismic methods provided a conservative set of 2 228 events, detected at milliseconds precision, described, and located. Besides, both trigger lag times and principal component analysis depicted correlations between environmental drivers and ice-cracking events. The automatic picker based on duration and trigger thresholds required manual supervision because of the initial numerous false signals that accounted for 96% of total initial events. Ice-cracking signals as short as 0.2s and frequencies of 8-40 Hz with an average power of -117 dB were statistically defined, classified, and described by case events as two types, associated, based on their spectral and temporal patterns, with the two ice break-up modes. With an estimated Rayleigh wave velocity of 680 m/s, all ice-cracking signals' locations were within the instrumented area.
Trigger lag times analysis improved detection and showed a strong link between ice-cracking events and drivers of lag times less than three hours, including near-immediate responses (<
2s). With multivariate analysis, the lag times showed a mainly climatic control for thermal melting and a primarily fluvial control in mechanical ice break-up. The combination of statistical and seismic analysis provides, despite the considerable manual screening, a valid and potentially site-transferable method to extract and describe ice-cracking signals and thus identify ice break-up modes in northern rivers.
Key words: Ice break-up, River ice, Subarctic, Thermal melting, Mechanical break-up,
Freeze-up, Environmental seismology, Seismic signals
Table of contents
1. Introduction ... 1
1.1 Background ... 1
1.2 Aim ... 3
1.3 Seismic approach and hypotheses ... 3
2. Materials and methods ... 4
2.1 Study site ... 4
2.2 General observations ... 5
2.3 Seismic instrumentation ... 7
2.4 Data processing ... 9
2.4.1 Event detection ... 9
2.4.2 Event location ... 12
2.4.3 Event anticipation and linking to drivers ... 13
2.4.5 Validation ... 14
3. Results ... 15
3.1 Ice cover evolution and break-up modes timing ... 15
3.2 Hydroclimatic context ... 16
3.3 Continuous seismic measurements ... 18
3.4 Ice-cracking events detection ... 19
3.4.1 Manual detection... 19
3.4.2 Automatic detection ... 20
3.5 Seismic characteristics of ice-cracking events ... 21
3.5.1 Waveform characteristics ... 21
3.5.2 Spectral characteristics ... 22
3.5.3 Inner clustering ... 22
3.5.4 Case event A ... 25
3.5.5 Case event B ... 25
3.6 Temporal activity patterns of events ...26
3.7 Location of events ... 27
3.8 Trigger lag times analysis ... 28
3.9 Principal component analysis of drivers ... 30
4. Discussion ... 32
4.1 Seismic signals identify ice break-up modes ... 32
4.1.1 Identifying ice break-up from continuous seismic signals ... 32
4.1.2 Seismic signatures of mechanical vs. thermal ice break-up ... 33
4.1.3 Seismic ice-cracking activity patterns and ice break-up ... 34
4.2 Seismic interpretation of ice break-up trigger mechanism ... 34
4.2.1 Trigger mechanisms for ice-cracking activity ... 34
4.2.2 Environmental drivers of mechanical vs. thermal ice break-up ... 35
4.3 Potential of seismic sensing methods in studying ice break-up ... 36
4.3.1 Seismic results validity and uncertainty ... 36
4.3.2 Seismic methods efficiency in ice break-up study ... 37
5. Conclusions ... 39
Acknowledgments ... 40
References ... 41
Appendix ... 45
1
1. Introduction
1.1 Background
Ice break-up is a fundamental process for subarctic rivers, which form and maintain an ice cover for up to five months (Barrows and Horton, 1907). However, ice cover dynamics vary significantly because of seasonal and interannual hydro-climatic conditions (Prowse et al., 2011, 2015). Ice break-up combined with the spring snowmelt flood, despite their brevity (< 2 weeks), are forceful in shaping the channel form, intensity of fluvial processes, and riparian vegetation (Ettema, 2002; Lind et al., 2014). The rapid ice removal and the high erosional power of ice control mainly the magnitude of the impact. Flow velocities are more variable under ice, with an altered wetted perimeter, and may be pressurized, shifting erosion and deposition sites, with changes in the velocity profile (Demers et al., 2011, 2013; Lotsari et al., 2017, 2019; Polvi et al., 2020). Moreover, the long-term sediment transport patterns are subject to change in their timing, magnitude, and type of mobilized sediment (Kämäri et al., 2015, 2018). Overall, ice break-up with the snowmelt flood represents the most dynamic period for turbulent flow and sediment transport in subarctic rivers. Thus, the sediment flux of a minimum annual of 56.2 ± 0.01 t/km²/a, estimated for the first time for a Fennoscandian catchment during ice break-up by Polvi et al. (2020), is ~5% higher than previous estimations that excluded this period.
Ice-covered, high-latitude rivers will change the most under an altered climate compared to the temperate rivers, as warmer winters will promote a thinner ice cover with low mechanical strength (Prowse et al., 2011, 2015). A later freeze-up, with more frequent and faster ice melting, promotes hydrologic connectivity, erosion, and sediment transport (Prowse et al., 2015). Ice break-up severity and timing are likely to change the most under warmer and wetter climate by shifting to a sooner date (two weeks earlier by the end of the century; Prowse et al., 2002), mainly as thermal decay with potential mid-winter mechanical break-up (Beltaos, 1997). These changes can have domino effects on ecosystems, influencing biological seasonality, mobility, and trophic structures (Huusko et al., 2007; Prowse et al., 2015). For example, the loss of under-ice shelters and access to new, open waters affects aquatic habitats and biotic diversity (Vincent et al., 2009; Stickler et al., 2010).
There are two main types of river ice break-up: thermal (over-mature) and mechanical (premature or dynamic) (Deslauriers, 1968). A thermal break-up is a mild, in situ ice decay, where ice melts and disappears through attrition (Beltaos, 1997). This gradual ice loss results from increased air and water temperatures (Beltaos et al., 1996) and the flow separation above and below ice cover with erosional effects (Eckstaller, 1988). In situ, ice melting can promote erosion as the ice-trapped sediment is released and, thus, increase bedload (Turcotte et al., 2011). On the other hand, since ice decay is slow and starts at the middle of the channel (i.e., where ice is thinner and weaker), ice can protect the banks from erosion (Tananaev, 2013). The more eventful mechanical break-up starts with ice cover movement under increased upstream discharge, followed closely by the fracturing of ice into blocks (Beltaos et al., 1996) and continued with repeated impacts between ice blocks and between ice and banks (Kozitskiy, 1989). As shown by laboratory flume experiments, the thin ice cover edges are the most vulnerable, and if submerged, they quickly break down (Michel and Abdelnour, 1975).
Commonly, thermal decay has already started before the onset of the mechanical break-up.
However, high, premature spring flows can fracture and entrain ice blocks at a maximum
observed thickness of 0.7 m (Barrows and Horton, 1907). The dynamic ice break-up can
promote geomorphic change through abrasion by mobilizing the ice stuck to the streambank
(Turcotte et al., 2011), and jagged ice blocks scar the channel bed and banks. Also, ice blocks
can jam, and gouging events may occur, scouring the river bed (Ettema, 2002). There is often
a blurred boundary between the two main types of ice break-up. Ice-cracking may result from
2
temperature increase and an initial ice decay as in thermal break-up, yet ice blocks and increased erosion follow like in mechanical break-up (Beltaos, 1997).
Although the effects of ice cover are currently relevant in subarctic rivers in the climate change context, field studies are challenging. On thin ice, any investigation would be dangerous, making the ice break-up process understudied. Moreover, since ice-cracking is a spontaneous, quick process, the timing of observations is critical. Despite field shortcomings, there has been an interest in studying river ice-cover in high-latitudes since the early 1900s (Barrows and Horton, 1907). Laboratory flume experiments were the most common methods to study the characteristics and effects of ice cover on fluvial processes (Lau and Khrisnappan, 1985; Urroz and Ettema, 1994; Tsai and Ettema, 1996). Typical approaches to study ice break-up include monitoring with remotely sensed imagery, as well as field cameras (McGinnis and Schneider, 1978; Cooley and Pavelsky, 2016; Beaton et al., 2019). However, these methods have limitations, such as the satellite imagery resolution, which can provide only a general activity pattern. Drone and time-lapse photos require high costs and maintenance, with an often slow and biased visual analysis. None of these established methods have provided continuous, high- resolution data about ice break-up mechanisms and link the two end-member modes (i.e., mechanical vs. thermal) to drivers under natural conditions.
To better understand the effects of ice break-up, improved methods with precise constraints on the information about timing, location, evolution, and triggers of ice cracks are needed. The approaches provided by the emerging field of environmental seismology can potentially solve the knowledge gap. These methods employ the information from seismic noise, traditionally associated with earthquakes, to study surface processes (Larose et al., 2015; Burtin et al., 2016). Portable seismometers placed underground record the data by capturing signals emitted by processes. Seismic approaches significantly narrow the temporal resolution and increase precision, especially if multiple stations form a network (Burtin et al., 2013). For example, subtle processes lasting less than one second are detectable (Dietze, 2016).
Environmental seismology has experienced significant scientific progress in both theory and application, especially over the last decade (Larose et al., 2015). This interest is because of the bridging of seismology with hydrology, geomorphology, meteorology, and other geoscience fields. Also, seismic sensors and data processing tools have recently become more accessible and sensitive (Dietze, 2018). Researchers have successfully used seismometers worldwide in dynamic landscapes to investigate geomorphic processes, such as rockfalls, rock avalanches, landslides, debris flows (e.g., Dammeier et al., 2011; Hibert et al., 2011; Lacroix and Helmstetter, 2011; Burtin et al., 2013; Ekström and Stark, 2013; Dietze et al., 2020) and fluvial processes (Burtin et al., 2008; Gimbert et al., 2014; Roth et al., 2014; Anthony et al., 2018;
Schmandt et al., 2017; Polvi et al., 2020).
Recent studies used signal analysis as a complement to exploring icequakes and crackings in
glaciers (Hammer et al., 2015) and river ice crack activity in the subarctic Sävar River (Polvi et
al., 2020). Previous results show that an ice crack generates high vibration frequencies with a
short duration, and Polvi et al. (2020) has put forward the link between ice cracking activity,
ice break-up type, and environmental drivers. The authors matched a sudden increase in ice-
cracking activity and intensity with a discharge-driven mechanical break-up. In contrast,
evenly spread out, cyclic, and less spiky signals were associated with the thermal break-up
controlled by climatic factors. The commonly used method to pick discrete events from seismic
noise is the ratio between short and long-term average signals (Allen, 1982). Alternative
methods, such as trigger lag times analysis, are now available to better constrain the results
from the classical picker and provide complementary information (Dietze et al., 2017; 2020ab).
3
Despite the rising interest in the study of river ice over the past three decades, there are still significant research gaps to fill in. In this study, I will further explore the potential of seismic sensing to investigate the ice break-up modes, as initiated by Polvi et al. (2020).
1.2 Aim
The primary objective was to provide a method to identify ice break-up using seismic signals in subarctic rivers and link the process to drivers. In this general-purpose context, I aim to answer additional research questions related to the spectral properties and the environmental context of the break-up processes:
1) How distinct are the signal signatures and average patterns of ice-cracking activity at a high temporal resolution for mechanical vs. thermal ice break-up?
2) How do the external triggers of ice break-up differ for thermal vs. mechanical mode?
By applying two distinct methods of seismic data analysis, I present a development of the technique introduced by Polvi et al. (2020) to estimate ice-cracking activity during mechanical vs. thermal break-up from seismic data. This approach gives insight into seismic methods’
performance in partitioning ice-cracking signals from similar discrete events and continuous signals associated with sediment transport and discharge. The first method applies the classical threshold-based algorithm for discrete events, while the second approach tests lag times information between ice-cracking signals and drivers, joined with multivariate analysis.
Thus, the methods contribute to information beyond the properties of individual events, such as the potential links between ice-cracking activity, ice break-up type, and environmental drivers, potentially extendible to a broader research context (i.e., high-latitude rivers).
The study period covered four main transition phases, two of ice cover formation (freeze-up) in late autumn and two major ice break-up phases in autumn and spring. Between these periods, stable, floating ice cover conditions were present. In estimating ice-cracking activity, the assumption is that most events are during transition phases, especially during the ice break-up, but ice-cracking events occur during the entire study period.
1.3 Seismic approach and hypotheses
I test four main hypotheses based on the objectives mentioned above.
H1: There is a correlation between a sudden increase in spectral intensity and the number of ice-cracking signals and mechanical ice break-ups during ice transition periods.
H2: Thermal ice break-up matches an even distribution of ice-cracking signals over the entire ice-covered period with diurnal rhythms and no significant change in the intensity and number of signals during transition phases.
H3: Thermal ice break-up is controlled mainly by climatic factors (i.e., changes in temperature and precipitation) and secondarily by fluvial processes.
H4: Mechanical ice break-up is mainly linked to fluvial processes (i.e., changes in flow discharge magnitude, turbulent flow, and sediment transport) and a lesser degree to climatic factors directly (i.e., storm or wind events and sudden temperature change).
To test these hypotheses, I analyzed seismic signals in time-frequency domains, locating and
correlating them with drivers. When ice cracks, a short, individual, and recognizable seismic
signature appears within a high-frequency band. This template enables the quantification of
the spatiotemporal distribution of the event activity. Since I treat an ice crack here as part of
mechanical ice break-up, events’ activity can distinguish between the two ice break-up types
and tests H1 and H2. Individual ice-cracking signals travel as surface waves with a set velocity
recorded at different times by seismometers (Dammeier et al., 2011; Burtin et al., 2016). As a
result, I can potentially compute the source location on the time offset of signals.
4
Although external forces activate the tension release within ice cover along the weakest path, there is little insight into the triggering mechanism. Adapting the definition of Wieczorek (1996), triggering factors can be hydroclimatic processes that cause an increase in stress and a near-immediate effect as ice-cracking. These drivers contribute to ice break-up differently for thermal vs. mechanical end-member modes (Polvi et al., 2020). However, a clear distinction is challenging because of the suddenness of ice cracks and overlapping environmental factors (Dietze et al., 2017). Thus, I combined the interpretations of lag times between factors and ice- cracking events with multivariate analysis, resulting in patterns of cause and effect relationships (Dietze et al., 2017, 2020a). I expect a dominant near-immediate seismic response to hydro-climatic factors when mechanical ice break-up occurs and longer lag times during thermal decay, supporting H3 and H4. Therefore, seismic data interpretation can give new and explicit information on the ice break-up mechanism.
2. Materials and methods
2.1 Study site
The study reach was preselected and seismically validated by Polvi et al. (2020) in a straight section of the subarctic Sävar River (64.1613° N, 20.3576° E), situated 40 km northwest of Umeå, in northern Sweden (Figure 1a and inset). For the study of river ice break-up, the location is ideal since ice starts to decay on straight reaches first (Beltaos, 1997). Moreover, the river reach is unregulated, with an overall symmetrical cross-section, free from channel bars, islands, large boulders, and significant in-stream vegetation (Polvi et al., 2020, Figures 1b, 2).
Thus, the site is an adequate natural laboratory for ice break-up with sustained and stable ice cover conditions in winter, not hindered or delayed by channel form. Also, the location was suitable for seismic measurements since it was remote, undisturbed by anthropogenic noise during the study period. There were two roads with minimal traffic, one at ~ 200 m from the riverside and a dirt road at ~ 200 m downstream of the study reach (Figures 1b, 3).
Measurements by Polvi et al. (2020) show a bankfull width of ~ 20 m, depth of ~ 2 m, and a bed slope estimated to 0.0045 m/m (Figure 1b). The mild bed slope suggests a uniform and relatively low-turbulent flow regime. The bed sediment is composed mainly of sand and gravel and a smaller proportion of coarser sediment and has a median grain size of 13.5 mm, (Polvi et al., 2020). Sometimes, the transport of gravel or cobbles during high flows can mask ice cracking signals or generate false events (Dietze, 2016).
Coniferous forests (Pinus sylvestris, Picea abies) dominate the land cover around the reach, with a narrow riparian zone (Carex spp., Alnus incana, Betula spp.) (Figures 1-3). The evergreen trees can hold snow on the branches until they sometimes break. Although the temperate species are barren, they are exposed and dry enough to crack under high wind pressure. In both cases, signals that mimic ice-cracking may appear within the seismic profiles.
The climatic context is part of the coastal-northern Sweden, where the average annual
temperature is 2.8° (-6.4°C for winter; average winter monthly temperature ranges from -
13.8°C to -2.7°C) (SMHI, 2020). The annual precipitation, mostly as snow, is 500-600 mm
(SMHI, 2020). Hence, in winter, ice cover is commonly snow veneered. For the study period
(23 October 2019- 24 April 2020), the average air temperature was -0.89°C, with an hourly
average ranging from -18.1°C to 14.7°C and total precipitation was 228 mm. Under these
conditions, ice-cracking occurred over the entire period with stable ice (from 4 November 2019
to 18 April 2021) in episodic and non-uniform distribution.
5
Figure 1. Study site location and context: a) Sävar River (blue) study reach location in Sweden with county borders (inset) and the study reach location within the river valley on hillshaded 2m digital elevation model (DEM from Lantmäteriet ©); managed forest is with green-grey shades; agricultural areas with white patches; roads with dark grey lines; b) zoomed-in study reach location (highlighted); transect and channel profile measured by Polvi et al. (2020) with red line; roads size indicated by line stroke size.
2.2 General observations
A time-lapse camera installed at the study reach was taking pictures every five minutes. The camera, located on the left bank, was in line with two seismic instruments (UME01 and UME03, Figures 2-3). The photos captured ice formation, stable ice cover, and ice break-up events.
However, because of malfunctions, these captured events were only within the period, 23
October 2019 (12.55) to 13 January 2020 (6.05). Although the spring ice removal lacks
photographic evidence, a complete ice break-up event occurred in November 2019 that can
provide seismic method verification. Moreover, Copernicus Sentinel-2 satellite imagery during
6
18-21 April 2020 can pinpoint the spring ice break-up as independent evidence to test the hypothesis that seismic sensing can identify ice break-up.
Figure 2. Study reach captured by the BUSHNELL time-lapse camera from the left bank, aligned with two seismic stations, UME01 and UME03 (time: 06.45.02, 29 November 2019, during the second freeze-up period).
Figure 3. Study site in 3D perspective with the instrumentation (seismometers and the camera) on Sävar River (blue); the highlighted area is the camera-captured reach; the primary and dirt roads with the forest path are with grey lines (stroke size relative to the road size); managed forest area is green-grey shaded, open fields with white patches; Geodata: 2m DEM from Lantmäteriet ©.
7
The weather data for the study period was gathered from the Swedish Meteorological and Hydrological Institute (SMHI). It included precipitation (15-minutes), air temperature, wind speed (hourly on 10-minutes averaged data), and radiation (1-minute and hourly). Two weather stations were data sources: Petriträsk, a station 55 km northwest of the site (64.5669°N, 19.6983°E) and Umeå Airport (63.7947° N, 20.2918° E), 40 km south of the site (SMHI, 2020).
I extracted the start time of meteorological events: wind, precipitation, wide temperature fluctuations (>1°C in 5 minutes), used in validating trigger lag times, from photographic evidence. Hydrologic data were from a calibrated gauge (Ytterträsk Nedre #2236) that measured discharge every 15 minutes. The gauge was 8 km upstream of the study reach, which creates ~ 9 hours time-lag of water flow. Diurnal freezing-melting dynamics specific to the study reach would make time-corrections biased (Polvi et al., 2020). Since the time lag would not significantly bias the results, I used the discharge data from the gauge directly.
2.3 Seismic instrumentation
Three compact seismic sensors sensitive to high-frequency signals captured seismic activity.
The PE6/B three-component 4.5 Hz geophones can record continuous seismic activity data in 3D (vertically, northerly and easterly) (GIPP, 2016), resulting in a large dataset (Figure 4). Two of the geophones were at distances of approximately 10 and 20 m from the right bank in a linear array (UME01 and UME03, respectively). The third deployed station was 40 m upstream from the first two sensors and 30 m from the right bank (UME04), forming a triangular network (Figure 3). Seismic setup and data collection were, in general, the same as in the previous study (see Polvi et al., 2020). The small and nested network (< 50 m sensor spacings) was suitable for ambient monitoring with high data resolution, coherence, and reduced uncertainty levels (Lacroix and Helmstetter, 2011; Burtin et al., 2016; Polvi et al., 2020). A 50 cm deep hole hosted each seismic sensor, fixed by three spikes into the soil. In a plastic box above the ground were data loggers, connected by cables to the sensors, and 12 V AGM batteries. The data logger Omnirecs Digos DataCube3ext recorded ground velocity signals at 400 Hz, with a gain of 32 (GPS flush time 55 minutes; GIPP, 2016; Figure 4).
Seismic recording started on 23 October 2019 (10.00 with a first restart at noon and a second one at 1.00 on 24 October 2019) and ended on 24 April 2020 (6.00). However, not all three sensors could cover the entire monitoring period, with data gaps caused by the battery and logger malfunctions. UME01 has an almost complete record (185 days), with one or two minutes gaps in a total of 7 days (Figure 5). Data gaps are usually around midnight, which might indicate a logger error when a date change occurs. As the total gaps in the UME01 dataset account for ~18 minutes, during the night when ice cracking activity is the lowest, the dataset is reliable for the analysis. UME03 had the second most complete record (116 days), and the UME04 had the least (85 days; Figure 5). In this study, I consider the monitoring period of 185 calendar days and 184 days of recorded time (23 October 2019–24 April 2020).
The three sensors had a coherent response to both fluvial and ice processes. However, there were notable differences in the frequency ranges and degree of overlapping of the recorded signals. These deviations can be attributed to the higher sensor sensitivity closer to the source (Burtin et al., 2013; Anthony et al., 2018). As a result, the UME01 data, with the most continuous record period and closest to the river, might be the most efficient in detecting discrete events with high amplitudes and frequencies (e.g., ice cracks). Even so, an ice-cracking event should be intense enough to be recorded by all stations. Thus, to discriminate between different types of events, where available, the data from all three sensors are concurrently used.
To locate events, I used a 2 m digital elevation model (DEM) obtained from Lantmäteriet.
8
Figure 4. Seismic instruments deployment at Sävar River study site (at the right bank for UME01, UME03, and UME04). Three-component geophones placed underground (~50 cm deep in hand-dug holes) and connected to a data logger and batteries installed on the ground recorded signals from surface processes. Instruments are from the Geophysical Instrument Pool Potsdam (GIPP).
Figure 5. Monitoring period (i.e., study period) for the three seismic stations with the data gaps caused by logger and battery malfunctions (blank sections). UME01 has an almost complete record of the study period (185 days) with minimal data gaps. UME03 had 116 days recorded, and the UME04 had the least of 85 days.
9 2.4 Data processing
The complete record of seismic data was processed using the R statistic software (RCoreTeam, 2020) with the package ‘eseis version 0.6.0’ (Dietze, 2018) (Figure 6).
2.4.1 Event detection
If the sampling rate of one station is 200 samples per second, the amount of recorded data in 185 days (i.e., study period) is over 3.1 billion unique values without data gaps. Thus, an automated method (e.g., script language routine) of identifying potential ice-cracking events is needed.
Signal processing must be efficient yet simplistic enough for a short computation time.
Ice-cracking events are detectable by targeting the seismic signal they generate. The query of seismic signals is a classical method to detect and study earthquakes that dominate seismic profiles below 5 Hz and under 1 Hz for larger-sized events with a distinct signal shape (Dietze et al., 2017; Figure 7a).
Anthropogenic noise shapes are relatively easy to distinguish, with signals gradually rising and falling relative to the source movement. On the other hand, the frequency range can be broad and activity-specific, making automatic detection challenging (Dietze et al., 2017). Rockfalls generate sharpened, often composite signals within a reasonably wide frequency window of 2 to 60 Hz and gradual decay (Lacroix and Helmstetter, 2011;
Hibert et al., 2014; Dietze et al., 2017; Figure 7b), with lower frequencies for rock avalanches (Suriñach et al., 2005). Ice-cracking creates acicular waveforms, shorter than any of the mentioned examples (less than 1s), with sudden onset (Polvi et al., 2020; Figure 7c). However, with many other discrete processes having overlapping signals, accurate estimation of parameter settings can be challenging.
To extract the onset time and duration of ice- cracking events, I ran the classic STA/LTA algorithm (Allen, 1982) on the full seismic record provided by the three geophones. The algorithm continuously computed the ratio between a short- term average (STA) to a long-term average (LTA) moving signal of the seismic envelope (Figure 7 steps 1 to 4). The two window sizes determine the sensitivity of the picker algorithm to sudden rises in ground velocity (Dietze, 2017).
Figure 6. Seismic data analysis workflow; steps 1-7 are part of the event detection, description, and quantification;
steps 8-10 are event location; steps 10-14 are event anticipation; step 15 with manual screening applied to all three processing stages.
10
Figure 7. The seismic portrayal of different types of events with seismic signals (left column) and spectrograms in time-frequency domain (right column); a) Earthquake generated by the Mw7.9 Denali Fault in Alaska and recorded by the broadband station near Parkfield- San Andreas Fault (Peng et al., 2018) (3November 2002, 22.30); b) Rockfall in Lauterbrunnen Valley, Bernese Oberland, Switzerland recorded by Dietze et al. (2017) (12 October 2014, 22.45); c) Ice-cracking recorded by the UME01 station installed by Polvi et al. (2020) on the right bank of Sävar River, Sweden (5 November 2019, 12.49); event start is marked on seismograms (left) by the red line and the end by the orange line; the warmth of color on spectrograms (right) depicts the seismic power (warmer colors show higher power in dB).
When a geophone recorded a discrete event, the STA value is much shorter, while the LTA value remains reasonably constant (Burtin et al., 2016; Figure 8). Hence, the ratio of the two averages increases when an ice crack occurs. After the event, the STA returns to levels close to the LTA value, decreasing the ratio. The algorithm has two additional parameters that need prior setting: a trigger value reflecting the starting point of an event (on-ratio) and a de-trigger value marking the end of the event (off-ratio). Choosing the correct values is critical, as the number and precise timing of picked events depend on it. I manually determined the algorithm parameters from seismic profiles using SWARM 3.2.0 (USGS, 2020), PQL II (PASSCAL, 2019), and SeisGram2K 7.0 (Lomax, 2014) inspection tools for visualizing seismic noise. I used two training seismic time series to manually choose values from the spring ice break-up of 2019 (8 March-28 April) and 2020 (1 March-24 April), with added control events from the dynamic phase of the autumn ice-break up (20-24 November 2019). As a reference point in manual screening, I relied on the experiences of Polvi et al. (2020) for the same study reach. Therefore, in this study, the STA was set to 0.45s and the LTA to 52s (Figure 8). The on-ratio of 6 and an off-ratio of 1.7 were the most efficient settings in detecting all ice-cracking events while minimizing misidentification. Freezing the LTA value was unnecessary, as the ice-cracking events are very short-lasting (< 2s), and the ratio deviates for cases with longer durations (cf.
Burtin et al., 2014).
11
Figure 8. Ice-cracking event detected with the STA/LTA algorithm from the seismic noise detected by UME01 station in the study site (Sävar River) in autumn 2019; end of the event with an orange line; short term moving signal average (STA) window is with yellow shaded in purple (0.45s), long-term moving signal average (LTA) window is black highlighted in blue (52s); on-ratio is of 6 and off-ratio of 1.7.
The optimized algorithm runs on the preprocessed vertical component (.BHZ) of each signal envelope because raw signals contain artifacts and not all ambient seismic signal is helpful.
Before each analysis step, I deconvolved the data (in m/s) to eliminate instrument self-noise from signals (cf. Dietze, 2016). Envelopes are square root values of the squared Hilbert transformed signals, depicting the seismic energy modes (i.e., amplitude) (Hibert et al., 2014;
Dietze, 2016, Figure 9). Given the typical ice-cracking events’ frequency above continuous signals, I filtered the envelopes to the 8-40 Hz band (third-order Butterworth filter) to restrict the analysis to the relevant information as identified from manual screening (Hammer et al., 2015; Dietze, 2016; Polvi et al., 2020).
Figure 9. Calculated hourly envelope on Hilbert transformed (cosine-tapered) signals with seismic energy (amplitude) as modes; seismic data recorded by UME01 station at Sävar River site (28January 2020, 15:00+1hr).
The large number of picked events, comprising 69 515 signals, required post-processing through statistical and manual verification to eliminate spurious results. Station co-detection, when data were available, reduced the number of false events automatically. I set a time window of 0.21s to 1.8s as a clipping mask based on the duration distribution of all identified events. The upper threshold value reflects the estimated maximum S-wave travel times within seismic networks for local rockfalls, seismically similar discrete events (Burtin et al., 2009;
Helmstetter and Garambois, 2010; Dietze et al., 2017). When the algorithm picked two or more overlapping signals with different durations or a time difference less than 2s (i.e., the maximum duration of an ice crack signal), I kept the first event. However, for event stacks (potential crack propagation), manual selection solved coincidence and sequences if co- detection was impossible. Not accounting for these cluster effects would bias results. High fluvial noise or overlapping weather events (e.g., precipitation and wind) are the principal potential erroneous sources.
Statistical variables defined objective exclusion of false events (Figure 6 step 5). The variables
efficiently defined rejection criteria based on the manually selected control events. A summary
of all the statistics used to classify the discrete seismic signals is in Appendix 1. I used
12
multivariate signal analysis to inspect the data fusion and create two major categories: ice- cracking events and false events, with subgroups and hierarchical branches. The principal component analysis (PCA) and hierarchical clustering algorithm (HCPC) estimated the variation, quantified as the standard deviation for each calculated and ranked variable (cf.
‘factoextra’-RCoreTeam, 2020). On the newly defined categorical variable, I performed an analysis of variance (ANOVA) with Tukey’s HSD to test the significance of the groups and subgroups differences to eliminate significantly different clusters (i.e., false results) from the ice crack signals. Finally, a manual inspection of both categories of events allowed plausibility checks. Also, I tracked the spectral evolution of signals identified as ice-cracking events on spectrograms (a.k.a. power spectral densities or PSD) and seismograms (Figure 6 step 7).
Seismograms were useful to translate signal amplitude distribution and variation, describing the shape and position of identified events. I calculated PSD profiles on the 1-90 Hz range, following the Welch (1967) method (Figure 6 step 6). This approach averaged moving windows of 0.7s and 0.4s, each with 90% overlap (Figure 7c). A multitaper (window) method corrected each spectrogram to avoid an erratic rise in spectral power at boundaries (cf. Burtin et al., 2016; Dietze, 2016). To explore event activity, I stacked the filtered hourly time snippets with a buffer of 60s, with 10s-long windows and 6s-long sub-windows with 80% overlap (Figure 10). To obtain profiles of the entire study period, I aggregated the PSDs with 1h-long windows and 6min-long sub-windows, overlapping at 50%. Because ice-cracking signals stand out from the noise on the spectral profiles as bursts of seismic energy above 8 Hz with an exponential decay (Polvi et al., 2020), I used time-frequency to trace picked events and potential cyclicity.
Figure 10. Spectrogram or power spectral density profile (PSD) calculated with Welch’s (1967) method on deconvolved and filtered seismic signals (1-90Hz, on the vertical component .BHZ) recorded at Sävar River (UME01 station on 6 April 2020, 13.00+1hr). The profile is composed of stacked 10s-long windows and 6s-long sub-windows, overlapping (80%). Buffer was 60s-long and with applied multitaper. PSD shows frequency (Hz) and seismic power (dB) by time; the warmth of color correlates positively with seismic power.
2.4.2 Event location
The sudden onset, high frequencies, and amplitude of emitted ice-cracking signals make the seismic location of their sources complicated. I used the signal migration approach to map ice- cracking signals’ origin because estimations are on full, filtered signal envelopes, with cross- correlations, suitable for discrete events (Burtin et al., 2014; Dietze, 2018). The method is probabilistic and calculates locations on the 2 m resolution DEM, based on Rayleigh seismic waves travel time. The setting of wave velocity required prior testing of values between 650 and 780 ms
-1with a 10ms
-1interval (Figure 6 step 8). I selected the range based on the results of Polvi et al. (2020) and on the estimation of apparent wave velocity from spectrograms showing signals amplitude and time delay of a moving animal from UME01 to UME04. I selected the average wave velocity of 680 ms
-1since it has the highest signal cross-correlation between stations in the desired frequency window. As previously done in studies, I assume a spatiotemporally constant velocity (Burtin et al., 2013; Dietze et al., 2017; Polvi et al., 2020).
I computed time offsets, as a function of distance, between event occurrence and signal
recording for station pairs on each grid pixel (Figure 6 step 9). The estimation is in three
dimensions as a straight line, as a sum of segments with sizes based on the DEM pixel
13
resolution and elevation (Dietze et al., 2017). The pixel with the highest pairwise cumulative probability density value is most likely the origin of the ice-cracking event (cf. Dietze et al., 2017; Polvi et al., 2020). I used only the values above 0.95 quantile of the coefficient of determination R
2(i.e., 95% probability) to estimate the spatial distribution of events and measurement uncertainty levels (cf. Burtin et al., 2014; Dietze et al., 2017; Polvi et al., 2020, Figure 6 step 10). Similarly, I removed events outside the instrumented area, and the remaining cases (2 228 events) provided subsequent validation. In the last step, the summed results provided the location map. However, this signal migration approach is suitable only for a general location exploration and validation of ice-cracking events. With data from all three stations, signal migration located 40% of events successfully within the R - ‘eseis’ environment (Dietze, 2018).
2.4.3 Event anticipation and linking to drivers
I placed the cropped event data set into an environmental context by assigning drivers to events. Ice-cracking as a tension release mechanism occurs when the sum of forces that provide cohesion is less than the sum of destabilizing forces (Wieczorek, 1996). The relationship between the two force types can transition towards instability suddenly or gradually (Dietze et al., 2017). Some factors induce a quick ice reaction, with almost equal timing between factor action and ice-cracking event, while drivers with a slow stress buildup on ice have lag times until ice cracks. The significant external stimuli on ice-cracking, selected based on the knowledge of the study site (Figure 6 step 11), are presented in Table 1. Given the low, constant altitude and fairly small study area, I treated the parameters as spatially uniform and neglected any potential local variatbility.
Factor type Factor proxies (variables) Unit ID
Climatic Precipitation event mm/15 minutes P
Wind speed m/s (hourly average) W
Air temperature (positive) °C (hourly average) T Air temperature increase °C (difference in 2 hours) Td
Longwave radiation W/m
2(1 minute) Rl
Direct irradiance Wh/m
2Rd
Hydrologic Discharge m
3/s (15-minutes average) Q
Table 1. Environmental factors and proxies selected for trigger lag times analysis based on the study site characteristics, available hydroclimatic data, and significance level (p-value <0.05 and D-value > 95% CI);
factors are potential ice-cracking event drivers; IDs represent the scaled factors in multivariate analysis (PCA);
units are from SMHI data (2020).
Climatic triggers can have both sudden or delayed ice-cracking seismic responses. The rhythmic impact of snowfall on ice can drive a quick reaction, while the snow accumulation may show a delay until, eventually, the weight is enough to break the ice. I defined the lag times as the time offset between a precipitation event end time and the closest ice-cracking signal start time. A precipitation event end time (trigger phase) was when the amount exceeded 0.1 mm/15 minutes (minimum quarter-hourly value for the study period). A higher threshold included 0.2mm/15 minutes. Wind can fluctuate pressure and thus stress, on exposed ice cover, including sudden block movements that may fracture the ice (Barrows and Horton, 1907). Hence, the reaction should be immediate. I selected wind speed as a proxy trigger and estimated the minimum time difference to the onset of an ice-cracking event.
The ice response-time to thermal stress directly depends on the magnitude of temperature
fluctuation. A quick and large-scaled thermal change (e.g., from negative to positive values)
decreases time lags and prevents ice fissure closing, with potential cracking propagation
14
(Barrows and Horton, 1907). Multiple proxies of the thermal effect were meaningful. Thus, I selected the mean positive-air temperature, air temperature fluctuations, and amount of longwave and direct irradiance (i.e., sunlight exposure duration). Calculated lag times were between the end of a thermal stress episode (i.e., a gradual temperature increase period or a sudden temperature change of minimum 0.1°C within 2 hours) and the onset of seismically detected events. Cumulative daily exposure to sunlight yields a positive response in terms of ice-cracking activity (Beltaos, 1997), and the estimated lag times were between an ice-cracking event and the end of the longwave radiation and direct irradiance increase event. I set this lag time onset from the time the radiation proxies from the STRÅNG model system (SMHI, 2020) have a positive trend (over the minimum hourly direct irradiance of 0.1 Wh/m2 and minutely longwave irradiance value of 181.7 W/m2) to the start time of the closest ice-cracking. The hydrologic driver has a swift mechanical control on ice because of enhanced shear stress and blocks’ movement during high discharge (Ettema, 2002; Turcotte et al., 2011). Therefore, calculations were on the end timing of the closest discharge event based on two threshold values (minimum 15-minute discharge value of 3.258 m3/s and 3.288 m3/s). A further driver can be biotic with animal movement on stable ice cover, increasing the pressure locally on ice, resulting in new ice cracks or expansion of existing ones. According to the joint seismic and visual evidence, there were only a few cases of animal-induced ice-cracking over the study period. Anthropogenic drivers can cause immediate ice-cracking seismic signatures. In this study area, the human influence was minimal, and any time lag discussion would be speculative.
I checked the immediate seismic past (i.e., from one hour to one day prior) of each ice-cracking signal for drivers and their timing, calculating lag times from seismic and independent hydro- meteorological and visual data set (Figure 6 step 12). The set timing was UTC with two hours difference, without daylight saving adjustments for the 30th March 2020. For each ice- cracking event lag times, I studied the distribution of values with kernel density estimates (KDE) (cf. Dietze et al., 2017, Figure 6 step 13). The kernel, expressed as a window applied to decompose data into density curves (Galbraith and Roberts, 2012), must be carefully selected through multiple testing, as it can change the correlations (Dietze et al., 2016). I stacked the obtained kernel density curves on top of each other to find the best-fit values at maximum overlap (cf. Blaauw, 2012; Dietze et al., 2017).
In the last step of getting kinetic details, I applied a principal component analysis on the signals-factors data set (Figure 6 step 14). PCA with the ‘princomp’ and ‘factoextra’ R packages (RCoreTeam, 2020) allows a spectral decomposition to study the correlations between drivers and detected ice-cracking events. The calculated eigenvalues (in a scree plot) were from scaled (i.e., standardized) data to depict the percentage of variance each component (driver) explains.
The graphs of individuals and the biplots (i.e., variables with events) show the potential groupings of data and the driver correlations (cf. Kassambara, 2017). I use the Keiser-Guttman and Broken Stick criteria, based on the eigenvalues, to select the principal components that explain most of the variation in the data.
2.4.5 Validation
Visual analysis of five-minute interval time-lapse photographs and Copernicus Sentinel-2
satellite imagery was the independent approach to confirm the onset and type of ice break-up
(Figure 6 step 15; Appendix 2). Photographs are also evidence for the identity of the isolated
events where the visibility is good. Moreover, the location map shows if the identified signals
are correct, and the time lags analysis constrains the detection criteria. The Monte Carlo
simulation tested the significance of the difference between the empirical lag time distributions
and that of randomly selected events. I defined a set of 2 228 random events within the same
time range as the actual cases and repeated the simulation on lag times 1001 times (cf. Dietze,
2021). The Kolmogorov Smirnov test provided the significance of results through the sample-
15
sized distance values (D at 95% and 99% confidence levels) and probability values (p <0.05). I considered a trigger important for the analysis if both the D and p values were significant.
3. Results
3.1 Ice cover evolution and break-up modes timing
The study reach started to form an ice cover from 30 October 2019 (6.15), and the bankfast ice was 5-15% of the total reach area until 3 November 2019. With the moving frazil ice, the total ice-covered area rose to 60-80% during this early freeze-up (Figure 2). On 4 November 2019, when temperatures decreased to -13°C, the frazil ice was completely consolidated, marking the beginning of the first stable ice-covered period (Figure 11a). The exposed cover suddenly shifted on 5 November 2019 from 13.00 to 14.15, cracking the ice in sequence. The rough ice cover is snow-veneered since 7 November 2019 (Figure 11b). A thermal decay initiated the autumn ice break-up from 14 November 2019, when the air temperature was around 0 °C, and continued for the next eight days, forming spots of thin ice on the cover. The formation of these patches likely also indicates the water flow carving the bottom of the ice cover. Both thermal and mechanical break-ups acted, jointly, on ice on 23 November 2019, reducing the area of ice cover by 85% in one day, with an opened section prograding from the middle of the reach sideways (Figure 11c). Photographs showed thin, fragmented ice edges and many small ice cracks. Thus, from 24 November 2019, the river reach had only bankfast ice that continued to melt in situ and crack.
Figure 11. For Sävar River, as captured by the time-lapse BUSHNELL camera from the left bank, the study period included a) initial freeze-up (4 November 2019), b) first stable ice cover conditions (7 November 2019), c) autumn ice break-up (23 November 2019), d) second freeze-up (3 December 2019), e) long-term stable ice conditions (5 January 2020) and f) open-channel conditions after spring ice break-up (25 April 2020).
From 27 November 2019, moving frazil ice was back on the reach, reducing the open channel
area to ~40% until freeze-up began to stabilize them three days later at air temperatures of
-8 °C (Figure 11d). On 6 December 2019, the nearly closed reach had an area in the middle of
the channel with flow separation (i.e., flow below and above the ice). By the following day, the
finished freeze-up marked the start of the second stable ice cover period. This last freeze-up
for the study period, as photographs showed, experienced frequent ice-cracking events and
occasional melting. The static ice cover conditions, with no visible changes, lasted until the
16
spring of 2020 (Figure 11e). A snow cover was present from 9 December 2019 until the final ice removal and hid ice cracks. Slight ice decay caused by increased air temperatures and water erosion was visible when no fresh snow resided on top of the ice. As identified from Sentinel-2 satellite imagery, on 18 April 2020 (10.12), there was still a complete snow-veneered ice cover on the study reach, as well as upstream and downstream (Figure 12). Although there is no visible imagery for the next two days, the image from 21 April 2020 (10.59) shows an almost ice-free study reach (Figure 12). Thus, the mechanical ice break-up and wash out happened in this time interval. After this spring break-up, there were only open channel conditions.
Figure 12. Copernicus Sentinel-2 L2A true-color satellite imagery showing part of Sävar River valley and the study reach (blue inset) at the onset (18 April) and during (21 April) the spring ice break-up (2020) and ice clearance; used to validate the seismic methods to identify ice break-up; EU Sentinel Hub (Map Tiler) ©.
3.2 Hydroclimatic context
Average hourly temperatures during the autumn freeze-up periods ranged from -14.4 °C to
4.3°C, while during the first ice break-up (14–23 November 2019), the values were consistently
above -4°C (Figure 13a). On the critical day of 23 November 2019, with the fastest ice decay,
the temperature ranged between -0.4 °C and 0.5°C. A broader range was specific to the two
stable ice cover conditions, between -26.8°C and 7.5°C, while during the second ice break-up
(18-21 April 2020), positive temperatures dominated (-0.9°C to 9°C). Direct irradiance was
commonly between 0 and 40 Wh/m2, while longwave irradiance ranged from 181.27 W/m2 to
359.08 W/m2 during the study period. The highest values are mainly during March and April
2020, in the mornings and early afternoons (8.00-14.00), when values exceed 500 Wh/m2
and often 800 Wh/m2 for direct and 300 W/ m2 for longwave irradiance. The lowest values
were during winter and autumn 2019, with 0 Wh/m2 specific to the night and early morning
hours (20.00-5.00). Over the study period, wind speed usually varied between 0-1.2 m/s, with
higher values in April 2020 (> 4 m/s), peaking with 8.2-8.5 m/s on 1 and 13 April 2020
(SMHI,2020; Figure 13b). Values of 3.5-4 m/s were occasionally measured during the autumn
ice break-up (18-22 November 2019). During the stable ice cover periods, the wind speed
overall did not exceed 1.5m/s. Precipitation, mainly as snow, ranged from 0 to 1 mm/15 min.,
with the peak during the autumn ice break-up (20 November 2019) and in general minimum
values (~ 0 mm) in the winter and during the spring ice break-up (Figure 13c). The discharge
values ranged in the study period from 3.258 m³/s to a maximum of 5.242 m³/s, with typical
values around 3.5-4m³/s (Figure 13d). Low discharge levels (< 3.5 m³/s) were during the
17
winter until early April 2020, while the highest levels (> 5 m³/s) were during the two ice break- up periods (21-26 November 2019 and 19-21 April 2020). The discharge peaked at over 5.5 m³/s during 20-21 April 2020 (Figure 13d).
Figure 13. Seismic and hydroclimatic data time series for the study period: a) air temperature (°C); highest values are with orange arrows; b) wind speed (m/s); peak speeds are with green arrows; c)Precipitation (mm); peak events are with light blue arrows; d) Discharge (m³/s; uncorrected for time lag gauge-site); major increases are with dark blue arrows; e) spectrogram calculated with Welch’s (1967) method on deconvolved and filtered (1-90 Hz) seismic data (UME01), with 1h-long windows, 6min.-long sub-windows, each with 50% overlap, multitapered; shows frequency (Hz) and seismic power (dB) by time (color warmth shows seismic power); f) temporal evolution of ice-cracking events with daily rates of detected signals; blue circles are days with only type A signals, grey circles with a majority(>75%) of events of type A; white arrows are days with peak rates during the mechanical phase. Shaded areas are transition periods: white for freeze-up and grey for ice break-up (darker for dynamic phase); the dark blue circle is case event A and the red circle is case event B; conditions present were open-channel (O), freeze-up (F), stable ice cover (S) and ice break-up (B). Discharge and weather data (a-d) was from the Petriträsk meteorological station ~ 55 km from the study site (NW; SMHI, 2020).
18 3.3 Continuous seismic measurements
In the study site, there are an amalgam of sources generating high-frequency (> 1 Hz) seismic signals, depicted by the spectrograms for the study period of the data from the three stations (Appendix 3). The warmth of color correlates to the seismic energy of the signals (i.e., the warmer the color, the more seismic power in the specific frequency band; Figure 13e, Appendices 3-4). The evolution of waveforms in time is comparable for all stations, yet the amplitude and frequency ranges differ, decreasing with the source-receiver distance (e.g., differences can be up to 30 Hz between UME01 and UME03/UME04; Appendix 4). I focus on the data from the seismic station UME01 because it covers the entire study period with minimal data gaps (Figure 13e).
The power spectral density for the 185 days shows three main continuous seismic signals, with only subtle changes over time (Figure 13e). The first steady frequency band is below 2 Hz, excited when the bursts of signals occur. A second consistent frequency band appears between 8-15 Hz and is often in the background of the spikes in seismic energies. The third continuous seismic band is at 28-45 Hz, with slightly lower overall power than the other two bands (up to -20 dB). A temporary, wavy frequency band extends from the third one up to higher frequencies of 50-75 Hz, lasting until the autumn dynamic ice break-up onset. All stations depict the continuous frequency bands, which are easier to identify as the distance river- seismic station increases (Appendices 3-4).
Overall, most of the seismic information is in the 8-45 Hz band (in -120 to -140 dB), which expands to 8-70 Hz with the start of the spring ice break-up. The broadband increases of seismic energy that occasionally overlap the steady seismic bands can rise to 160 Hz (vertical patterns in Figure 13e, Appendices 3-4). The most prominent seismic bursts are on the days when the wind speed peaks (e.g., on 1 and 13 April 2020; Figure 13b) or precipitation increases (e.g., 20 November 2019 and 9 February 2020; Figure 13c). Ice-cracking events also appear as energy bursts mainly in the frequency band of 8-40 Hz, with excursions up to 70 Hz and short durations (< 2s), seen on PSDs on data from all stations (Figure 13e, Appendix 4).
Figure 14. Daily spectrograms calculated with Welch’s (1967) method on deconvolved filtered signals (1-90Hz), multitapered, with 10s-long windows and 6s-long sub-windows, overlapping (80%); show frequency (Hz) and seismic power (dB; showed by the warmth of colors) by time. With black arrows are the continuous seismic bands;
A and B show ice-cracking sequences with the case events A and B; 1 is an individual mimicking signal; 2- broadband bursts of seismic energies, masking ice-cracking events.
19
Figures 14 and 15 show PSDs at a higher temporal resolution of 24 hours and one hour, respectively, with the commonly recorded types of signals. The stationary frequency bands presented above are visible on both panels, together with the superposed energy bursts.
However, as we temporally zoom in, more diverse types of signals appear, such as the tilde- shaped signals (Figures 10, 15) or pennant shaped seismic bursts with frequency ranges of 8- 45 Hz, and increased power (up to +20 dB) in the 8-15 Hz range (Figures 14-15). In the same frequency band, shorter duration vertical signals stand out, identified here as ice-cracking events (Figures 10, 14-15). Thus, ice crack waveforms are potentially extractable from the continuous seismic signals, and among the stations, detection is most facile on UME01 data.
Figure 15. Hourly spectrogram, calculated with Welch’s (1967) method on deconvolved and filtered seismic signals (1-90Hz). The profile has 10s-long windows and 6s-long sub-windows, overlapping (80%). PSD shows frequency (Hz) and seismic power (dB) depicted directly by the warmth of colors) by time. Black arrows indicate the steady seismic bands; a is an ice-cracking event; 1 is a tilde-shaped signal of biotic or anthropogenic source; 2 shows broadband bursts of energies; 3 and 4 are panned-shaped signals, a burst of energy.
3.4 Ice-cracking events detection
3.4.1 Manual detectionThe manual screening of seismic records (hourly; sometimes 10-minutely) from both the previous study of the same site (Polvi et al., 2020) and current data, combined with visual analysis, created a set of 48 control events. The selected waveforms were short, needlelike, and reasonably easy to detect in the seismic profiles (Figure 7c). However, some signals were potential false triggers of the STA/LTA picker if the shape, duration, or frequency range resembles ice-cracking events. For example, an ice crack event is visibly different from the tilde-shaped signal, formed above the continuous frequency bands, with a long duration and gradual rise and fall in seismic amplitude (Figures 10, 15-16a). Thus, the picking algorithm ignores this signal that coincided with an otter's arrival and departure at the river reach (Appendix 5a). However, the evidence is insufficient to conclude that the source is biotic and not anthropogenic (e.g., car transit).
The dune-shaped signal with frequencies at 2-5 Hz, which appeared with an increase in
discharge, is also different enough to be ignored by automatic detection (Figure 16b). Also, the
signals aligned as bar codes over the steady seismic bands are not mimicking ice cracks, and
when its recording time is when photographs show a deer running along the right bank (Figure
16c, Appendix 5b). Although the waveform presented in figure 15 (3-4) is in the same frequency
band as an ice crack signal, the longer duration helps prevent its detection. However, the
broadband seismic bursts are likely sources of false events, as the shape, durations, and on/off-
trigger values mimic ice-cracking signals (Figures 14-15).
20
Figure 16. Spectrograms, calculated with Welch’s (1967) method on deconvolved and filtered signals (1-90 Hz;
vertical component .BHZ from UME01). The profiles have averaged moving 10s-long windows and 6s-long sub- windows, and 80% overlap, with a buffer of 60s and, are multitapered. Power spectral densities (PSDs) show frequency (Hz) and seismic power (dB) on the y axis and time on the x axis; a) comparison between a tilde-shaped signal of human or biotic origin, with a gradual arrival and departure from the instrumented area (1) and an ice- cracking event (A); b) dune-shaped signal within the microseismic band (2) and panned-shaped bursts of seismic energies generated by sediment transport (3); c) barcode-shaped signals which coincided with a deer running along the right bank (4).
3.4.2 Automatic detection
STA/LTA algorithm provided 69 515 initial ice-cracking events for the study period. The highly inflated number of events required the application of automatic rejection criteria. Signal co- detection and time window limitation reduced the number of events to 54 956 potential cases.
However, there were still many spurious signals in the data set, and thus I included more descriptors to eliminate them (Appendix 1). Hence, with the 68 calculated statistical variables, the PCA and HCPC analysis created ten initial large groups (Figure 17, Appendix 6).
Manual screening showed that only three of them had ice-cracking signals (clusters 3, 6 and 7 in Figure 17), which were further clustered in the same way, resulting in two groups with 120 unequal subgroups and 1 046 hierarchical tree branches (Appendix 6c). The clusters, manually checked against the control group and photographs, showed only one of them with ice crack signals, and thus, I used it for further analysis. In this group, five subgroups had mostly ice cracks and, thus, I selected them as additional control clusters after screening the false signals.
The rejected significantly different (p-value <0.05) clusters (n = 690) and subclusters (n = 79)
based on Tukey HSD test, reduced the number of potential events to 13 308. A final manual
check eliminated unnatural sequences of events, missed spurious signals, and events with
coordinates outside the study reach, leaving the data set with 2 228 identified ice-cracking
events (Figure 17). Most of the false signals that were difficult to separate automatically from
ice-cracking events were, as expected from manual detection, picked during days with high
wind speed, precipitation, and occasionally discharge.
21
Figure 17. The clustering of detected ice-cracking events with PCA and HCPC on 68 statistical variables, after the automatic rejection of outliers statistically and with co-detection (stage 2 after 1). Grouping is on the first two dimensions and clusters 3, 6, and 7 have ice-cracking events, with the rest being spurious. The bottom graph shows the workflow of eliminating false events from the initial set of picked signals.
3.5 Seismic characteristics of ice-cracking events
3.5.1 Waveform characteristicsThe remaining signals have specific seismic features which, though varying between events, have an acceptable consistency to characterize an ice crack (Appendix 7). Commonly, ice- cracking events last 0.4s-0.5s and rarely over 1.3s. The duration distribution shows that at 0.9s, the number of events has fallen to 28% of the maximum height (14% at 1.4s; Figure 18a). The shortest duration of an event was 0.21s, and the longest-lasting ice crack lasted 1.76s. The signal-to-noise ratio (SNR), necessary to locate events, is often over 4 for ice crack, but lower values can appear for smaller signals. Ice-cracks are events with sudden onset and rapid decay as the rise time (i.e., time until the signal reaches maximum amplitude) and fall time (i.e., time from maximum amplitude to the end) is commonly below 0.08s. However, some waveforms show values above 1s (a maximum of 1.4s) for both times, which expands the range of rise-to- fall time ratio very much (0-248) and might suggest two types of ice-cracking events.
The detected signals are usually symmetric, with skewness values around 0, and almost all cases are skewed moderately at most, between -1 and 1 (only 22 waveforms are highly skewed).
Of all events, 70% are leptokurtic (kurtosis >3), with the rest of the signals having shapes close
to a normal distribution (kurtosis = 2-3). These results of skewness and kurtosis are similar for
both waveform and calculated envelope. Below 8 Hz, the mean kurtosis is 2.6, while above 20
Hz is over 4. The majority (55%) of events have one principal peak, while the rest have several
small secondary peaks (up to 8). The maximum frequency of all signals is at least 1.5 times
above the median (24-27 Hz) frequencies of the waveform (Figure 18b; Appendix 7). The
seismic energy content had unequal distribution within 98% of the signals, with at least 1.4
times higher energy in the first third of waveforms than the last two-thirds (29 events had equal
22
values between the two sections). There was a group of events with a particularly energetic onset (> 3 times more energy in the first third).
Figure 18. Density curves of the waveform (a-b) and spectral characteristics (c-d) of the 2 228 potential ice- cracking events: a) duration distribution; b) median frequency (Hz); c) median seismic power (dB); d) maximum seismic signal power (dB); kernel bandwidth was 1.
3.5.2 Spectral characteristics
The median spectral power of the signals was -117 dB, higher towards the spectral peaks with -102 dB, and a maximum value of -66.9 dB (Figure 18c-d, Appendix 7). Commonly, the difference between the maximum and median spectral power is 20 dB, but in some cases, as much as 70 dB. The spectral flatness has low values, usually below 0.05 and for all events, less than 0.71, while mean entropy is over 0.67 (of maximum 1), both feature values indicating a significant difference between the signals and the background seismic noise. The spectral energy varies with the frequency, with at least 20% higher energy at 10-30 Hz than below 10 Hz and over 40% higher than the 1-3 Hz band. Within the 10-30 Hz, the higher the frequency, the higher the energy content. The frequency ranges from 6 Hz to 70 Hz, but most commonly (n = 2183) between 8-45 Hz (Figure 18b). Central frequencies gravitate around 10-38 Hz (n = 1774). A summary of the spectral centroid characteristics, power variance, standard deviation, and error, with other descriptors, are presented in Appendix 7.
3.5.3 Inner clustering
The PCA on the 68 statistical variables created 61 potential principal components, from which, based on the eigenvalues, the first five are significant (PCs with above-average eigenvalue and broken stick model bars in Figure 19). The first principal component is the most important, explaining 36% of the variation in the data on its own and 49.4% with the second component.
The first PC explains most of the spectral descriptors positively and negatively, and the second PC is better in explaining waveform features (e.g., kurtosis or envelope variables; Figure 20a).
The HCPC on these dimensions shows two separate groups of ice-cracking events (Figure 20b).
A similar clustering is present when the first PC is against the other significant components
23
(Appendix 8a-c). However, the grouping is no longer visible for the combinations of PCs if the first principal component is excluded (Appendix 8d).
Figure 19. Results of PCA on the detected 2 228 events with 68 statistical descriptors and the 61 principal components: a) Barplots with eigenvalues (%) for all dimensions and the first ten with highest explained variance (inset); the red line shows the average eigenvalue and the threshold cf. Keiser-Guttman criterion; b) barplot with eigenvalues (%) empirical (pale pink) and broken stick modeled (red) with raw values per PC in inset-table. The dimensions above the red line and bars are significant (i.e., the first five).
