Meteorological conditions

The air temperature obtained from the AWS station at 4550 km elevation shows a clear diurnal cycle during the survey period (Fig. 6a). The average air temperature over 9 d was about −5.2°C, with a maximum temperature of 1.4°C and a minimum temperature of −8.1°C. The highest and lowest temperatures appeared at 13 CST and 7 CST, respectively. The temperatures rose rapidly from 9 CST to 12 noon and dropped slowly from 17 CST to 7 CST the following morning. The local sunrise and sunset during our survey period was ~7:30 CST and ~19 CST, respectively. The wind speed (Fig. 6b) and glacier surface humidity also show diurnal cycles, which results in a deviation of the ice surface temperature from the air temperature based on glacier surface energy budget considerations (Sun and others, Reference Sun2014). However, comparison between air and ice surface temperature measurements (Sun and others, Reference Sun2014) in June–September 2011 showed a similar diurnal pattern and same daily hours with maximum and minimum temperatures, although the glacier surface temperature was systematically lower than the air temperature by 2–4°C. In this study, we use the air temperature as a proxy for the air surface temperature, due to the lack of direct ice surface temperature measurements during our survey period.

On 11 October 2015, snowfall lasted from 15 CST to 7 CST the following morning (Fig. 6). Approximately 5 cm of snow covered the LGNT until the end of our observation. We did not collect hydrological data, but note that no sound of water flow was heard in the field and the downstream riverbed was dry during the observation period. We only observed limited surface melting in the daytime close to the glacier fringes. Although the air temperature was mostly below 0°C during our observation period, the local solar radiation, wind, and cloud may lead to slight surface melting at noon. In addition, we found an ice moulin of ~5 m wide at a distance of ~100 m to the northwest of our seismic station S1 (Fig. 1b), close to the glacier terminus. The ice moulin was dry during our survey period.

Long-duration events

We detected a total of 158, 129, 65 and 49 long-duration events for the four stations in the 9 survey days (Fig. 4) based on the parameters in Table 1. Most of the long-duration events occurred between 12 noon and 18 CST in the daytime, when the air temperature ranged from −3 to −1°C. On the last day of our survey (i.e. 11 October), the long-duration events were also identified at night. The amplitudes of long-duration events are relatively low (Figs 2 and 5). Their waveforms lack impulsive onsets, and distinct body or surface wave arrivals. The duration of the long-duration events is tens of seconds to tens of minutes. We selected three representative long-duration events for detailed spectral analysis (Fig. 5), which shows a dominant frequency below 20 Hz. The long-duration events may also be associated with a higher secondary resonance frequency, e.g. 80 Hz in Figure 5a and 40 Hz in Figure 5b. Spectrograms for these two selected events also clearly show a temporal evolution of the resonance frequency (bottom panels of Figs 5a and b).

Fig. 4. Counts of detected long-duration events per hour for stations S1–S4. The correlation coefficients between the hourly event counts and air temperature (Fig. 6a), wind speed (Fig. 6b) and temperature change rate (Fig. 6c) are listed in the top right corners. Date and hour are in CST.

Fig. 5. Waveforms (top), spectrograms (bottom) and power spectra (right) for three representative long-duration events with starting times of (a) 11:35:11 CST, (b) 12:10:05 CST and (c) 15:14:04 CST on 6 October, detected at the seismic stations S2, S2, and S4, respectively.

Fig. 6. (a) Air temperature, (b) wind speed, (c) calculated temperature change rate at a depth of 20 cm and (d–g) hourly short-duration events detected at S1–S4. The correlation coefficients between short-duration events and air temperature, wind speed and temperature change rate are listed in the top right corners. Date and hour are in CST. We did not deploy station S1 on the first day.

Short-duration events

We detected 4514 short-duration events per day on average for an individual station (Fig. 6), which is hundreds of times more than the number of detected long-duration events. At each seismic station, the number of short-duration events was similar in the first 6 d of observation. After a snowfall, the daily number of short-duration events in the last 3 d of observation decreased to ~one-fifth of that during the first 6 d. Most (95%) of the short duration events lasted for 0.04–0.25 s (Fig. 7a). The waveforms of the short-duration events are dominated by Rayleigh waves without visible P- or S-phase arrivals (Fig. 2a). The dominant frequencies of the events range widely from 20 to 100 Hz, with local peaks at 30, 50, 75 and 90 Hz (Fig. 7b). The event counts increase with frequency, and the maximum frequency is close to the Nyquist frequency of 125 Hz.

Fig. 7. Histograms of (a) duration, (b) frequency and (c) local magnitude for the short-duration seismic events. (d) Corresponding cumulative distribution of the seismic magnitude in (c). Red curves correspond to the best-fit to the Gutenberg–Richter distribution.

Fig. 8. Estimates of the short-duration event counts after applying waveform association analysis to exclude the potential repeated counts. The events simultaneously recorded at two stations (light blue), three stations (orange), and four stations (red) are ~ 23, 4 and 0.2% of the events recorded at one station (blue).

The dominance of Rayleigh waves is consistent with the waveforms observed in previous studies (Neave and Savage, Reference Neave and Savage1970; Walter and others, Reference Walter2009; Dalban Canassy and others, Reference Dalban Canassy2013; Röösli and others, Reference Röösli2014; Carmichael and others, Reference Carmichael2015), in which the source mechanism has been attributed to surface tensile cracks within the surface crevasse zone (top ~20 m of the glacier). However, our detected events are shorter in duration, higher in frequency and more frequent than reported from other glaciers. This is mainly due to our chosen shorter time windows for the STA/LTA detector, and less significantly due to repeatedly counting seismic events at multiple stations. To reduce duplicate counts, we performed a waveform association analysis, similar to Carmichael and others (Reference Carmichael2015). This analysis identifies the events detected at different stations within a time interval less than the expected travel time across the network as the same event. However, due to the difficulty of locating the small short-duration events in this study, we can only use the maximum travel time between two stations, which was calculated as the distance between two stations divided by an assumed Rayleigh wave velocity of 1660 m s⁻¹ (Mikesell and others, Reference Mikesell2012). The waveform association analysis thus provides maximum estimates of the event counts that are recorded at two, three and four stations simultaneously, which are ~23, 4 and 0.2% of the total event counts recorded at one station (Fig. 8). The events recorded at multiple stations are expected to be larger in magnitude, and thus the rapid reduction in the event counts is in line with our expectation of an exponential decrease of event counts with magnitude. We plotted the hourly event counts detected simultaneously at two stations in Figure 3c. The diurnal cycle is clearly visible and thus is a robust observation.

The diurnal cycle is the most significant temporal pattern of short-duration events (Figs 6 and 9). Nearly all the short-duration events occurred between 19 CST and 9 CST the following morning, with peak event counts at 20–23 CST. The starting and ending times of the daily burst of seismic events, 19 CST and 9 CST, were confirmed by low p-values in paired t-tests (Fig. 9d). These starting and ending times have an uncertainty of 1 h, as it is based on hourly seismic event counts. After the starting time of 19 CST, the event counts rapidly increased to the maximum occurrence frequency of 363 h⁻¹ at 20 CST (Fig. 9c). The event counts then gradually decreased at night till 112 h⁻¹ at 9 CST the next morning. In the daytime, the event frequency is <40 h⁻¹.

Fig. 9. Nine-day averaged (a) air temperature, (b) temperature change rate at a depth of 20 cm and (c) counts of short-duration events per hour. (d) Probability values (i.e. p-values) of the statistical test (i.e. two-sample t-test) for the hypotheses that the neighboring 2 h have different mean seismic event counts. The notable low p-values at 19 CST and 9 CST indicate the initiation and ending hours for the daily burst of the short-duration events. The gray regions and dashed lines indicate the 95% confidence intervals (i.e. 1.96 times the standard error) for the hourly mean values.

The counts of short-duration events (Figs 6d–g and 9c) are strongly and negatively correlated with the temperature fluctuation (Figs 6a and 9a) with correlation coefficients of −0.43 to −0.32 for the four seismic stations. Most short-duration events occur when the air temperature decreases from −4°C at 19 CST to −7°C at 9 CST the next morning. The occurrences of short-duration events are also negatively correlated with wind speed (Fig. 6b). It has been suggested that wind or water-generated tremor noise in the daytime may reduce the detectability of seismic events (Röösli and others, Reference Röösli2014). We applied a highpass filter at 20 Hz to suppress the low-frequency wind or water-generated noise. Results (Fig. 10b) show that the short-duration events are still concentrated at night due to their high frequency and amplitude. Therefore, the low detection of short-duration events in the daytime may not due to increased tremor noise. Instead, a source mechanism preferentially operating at night is necessary to explain the occurrence of short-duration events.

Fig. 10. (a) One-day plot of the seismic records from 10 to 11 October at S2, and filtered seismograms after applying (b) highpass and (c) lowpass filters at 20 Hz.

Local magnitude of short-duration events

To characterize the energy/amplitude of the short-duration events, we calculated the local magnitudes (M _l) using the method of Richter (Reference Richter1935):

(1)$$M_{\rm l} = {\rm log}\left( {\displaystyle{A \over {A_0}}} \right)$$

where the parameter A is the maximum amplitude of displacement, and A ₀ is a correction factor defined by a reference zero magnitude earthquake that generates a displacement of 1 µm at a 100 km distance. We adopt the empirical A ₀ value provided by Richter (Reference Richter1958), which was originally estimated for southern California but is also used worldwide (Roux and others, Reference Roux, Marsan, Métaxian, O'Brien and Moreau2008).

Our estimated M _l of the short-duration events range from −2 to 0.5 (Fig 7c). Following Podolskiy and others (Reference Podolskiy and Walter2016) (citing Barruol and others, Reference Barruol2013), we found that the magnitude distribution fits a Gutenberg–Richter distribution (Gutenberg and Richter, Reference Gutenberg and Richter1944): log N = a − b × M _l. We estimated the parameters a and b to be 3.75 ± 0.02 and 0.99 ± 0.01, respectively (Fig. 7d). Our estimated b value is close to 1, consistent with previous estimated b-value estimates for surface crack events (Podolskiy and Walter, Reference Podolskiy and Walter2016 and references therein). Our fitting to the Gutenberg–Richter distribution implies a magnitude of completeness of −1.2 (Fig. 7d), below which a large portion of weaker events were not detected. We also grouped the detected short-duration events on an hourly basis to test the evolution of b value. However, no clear temporal evolution trend of the b value was found.

Thermally induced opening of surface cracks

Most of our detected short-duration events are expected to occur near the surface (maximum ~ 20 m depth) due to the dominance of Rayleigh waves (Walter and others, Reference Walter2009; Röösli and others, Reference Röösli2014; Podolskiy and Walter, Reference Podolskiy and Walter2016). In contrast, seismograms of deep or basal events have relatively large P-wave amplitudes without the presence of Rayleigh waves (Deichmann and others, Reference Deichmann2000; Dalban Canassy and others, Reference Dalban Canassy2013). Although it is impossible to determine the fracture mode (shear or tensile) of the seismic events due to the lack of P wave polarity observations, we suggest that the surface events are more likely to be tensile cracks as the tensile strength is much less than the shear strength (Schulson, Reference Schulson2002). Additionally, the tensile strength can be easily overcome at night due to thermal contraction (next paragraph). When local stress near the defects or crack tips exceeds the fracture toughness according to the fracture mechanic considerations (van der Veen, Reference van der Veen1998, Reference van der Veen1999), the formation or propagation of surface cracks can be recorded as seismic events.

The concentration of surface crack opening at night could be induced by thermal stress or refreezing meltwater. Daytime meltwater may accumulate in surface crevasses; refreezing of meltwater at night can trigger crack initiation and propagation (van der Veen, Reference van der Veen2007). An alternative mechanism is that reduced meltwater at night increases effective basal stress and thus induces spatial variations in surface strain (Carmichael and others, Reference Carmichael2015). However, the increased basal stress mechanism requires considerable amounts of meltwater at the base of the glacier, which is less likely to occur at LGNT in October.

Although we cannot exclude the effect of meltwater refreezing at night, we found that thermal contraction alone can induce enough tensile stress in the top tens of centimeters to promote tensile crack opening and seismic events (next section). This thermal contraction effect is supported by the significant negative correlation between the counts of short-duration events and the air temperature, and temperature change rate (Fig. 6). A similar negative correlation has been reported at a Himalayan glacier (Podolskiy and others, Reference Podolskiy, Fujita, Sunako, Tsushima and Kayastha2018) and interpreted as having a thermal origin. In addition, the snowfall on the 7th day of our observations is expected to damp heat conduction. The corresponding reduction in the thermal stress can explain the observed decrease in the icequake event number after the snowfall, which further corroborates the causal relationship between the thermal stress and short-duration events.

Thermal stress evolution

As the temperature decreases, the thermal contraction of the ice exerts tensile stress on the surface ice layer. Once the stress exceeds a threshold value, the opening of surface cracks occurs, which could be recorded as icequakes. To quantify the thermal stress in response to the diurnal temperature fluctuation, we first calculated the propagation of the temperature waves T(z, t) at depth z and time t based on the Fourier equation of heat conduction following Sanderson (Reference Sanderson1978) and Nishio (Reference Nishio1983):

(2)$$\eqalign{T(z,t) & = T_0 + \mathop \sum \limits_{s = 1}^N T_s{\rm exp}\left[ {-z{\left( {\displaystyle{{\omega_s} \over {2k}}} \right)}^{1/2}} \right] \cr & \cos \left[ {\omega_st-\phi_s-z{\left( {\displaystyle{{\omega_s} \over {2k}}} \right)}^{1/2}} \right]\;} $$

where the thermal diffusivity k is assumed to be 1.091 × 10⁻⁶ m² s⁻¹ for a constant ice density of 0.9 g cm⁻³ (Sanderson, Reference Sanderson1978). The temperature amplitudes T ₀ and T _s, frequency ω _s and phase ϕ _s depend on the Fourier transformation of surface temperature (assuming to be equal to the air temperature): $T(0,t) = T_0 + \mathop {\sum}_{s = 1}^N T_s\,{\rm cos}(\omega _st-\phi _s)$. The Fourier transformation was computed using the Matlab Fast Fourier transformation codes fft and ifft, in combination with Fourier coefficient conversion codes complex2real and real2complex (written by Boynton GM for the course of Introduction to Programing for the Behavioral Sciences, Autumn, 2015, University of Washington. Retrieved from http://courses.washington.edu/matlab1/matlab).

The short-duration events mostly occur when the temperature drops at a depth of 20 cm (i.e. when temperature change rate is negative). The correlation of the short-duration event counts with the temperature change rate at 20 cm depth (Figs 6c and 9b) is −0.56 to −0.54 for the four seismic stations, which is more significant than the correlation with the surface temperature fluctuation. The daily-averaged evolution of the temperature and temperature change rate within the top one meter of the glacier ice clearly shows diurnal variations (Figs 11c1 and d1). Figures 11c2 and d2 show depth-dependent temperature and temperature rate profiles at four representative hours. The diurnal temperature cycle only influences the top 1 m of the glacier, with rapidly decreasing amplitude with depth due to attenuation of temperature waves (black curves in Figs 11c2 and d2).

Fig. 11. Hourly short-duration event counts, temperature, temperature change rate and thermal stress averaged for the 9 survey days. (a) Counts of short-duration events per hour (same with Fig. 9c). (b) Calculated thermal stress evolution at a depth of 0 (red), 20 (green) and 40 cm (blue). The solid and dashed curves correspond to (e1–2) the first scenario without background strain rate and (f1–2) the second scenario with background strain rate, respectively. (c1–2) Temperature wave propagation and corresponding depth-dependent temperature profiles at four representative times. Black curve shows the lower enveloping curve of the temperature at various depths. (d1–2) Depth-dependent evolution of temperature change rate and corresponding temperature change rate profiles. Black curve shows the higher enveloping curve. (e1–2) Calculated thermal stress evolution for the first scenario without background strain rate, and depth-dependent profiles. Black curve shows the upper enveloping curve of the tensile stress at various depths. (f1–2) Calculated thermal stress evolution for the second scenario with a background strain rate. The horizontal tension is positive.

Sanderson (Reference Sanderson1978) considered that the surface ice layer cannot move as it is constrained both horizontally and vertically. Internal stress, instead of strain, is thus induced by temperature changes. This thermal stress was calculated by considering the strain rate that would occur if the ice body is free to expand or contract, and then estimating the stress required to reduce the strain rate to zero using Glen's flow law (Glen, Reference Glen1955):

(3)$${\dot{\varepsilon}} ^{\prime}_i = C\tau ^{n-1}\sigma _i^{^\ast}, $$

where τ is the second invariant of the deviatoric stress tensor $\sigma _i^* $ and C is the creep parameter. The best estimate for n is 3 based on laboratory experiments and glacial morphologic studies (Sanderson, Reference Sanderson1978; Cuffey and Paterson, Reference Cuffey and Paterson2010).

The final formula in Sanderson (Reference Sanderson1978) to estimate the horizontal principal stresses due to restrained thermal expansion and contraction are as follows:

(4)$$\sigma _1^{^\ast} = \left( {\displaystyle{{{{\dot{\varepsilon}}^{\prime}}_1} \over {C(1 + a + a^2)}}} \right)^{1/3},$$

(5)$$\sigma _2^{^\ast} = a\sigma _1^{^\ast}, $$

where ${\dot{\varepsilon}} ^{\prime}_1 = \dot{\varepsilon} _1-\alpha ((\partial T)/(\partial t))$ and ${\dot{\varepsilon}} ^{\prime}_2 = \dot{\varepsilon} _2-\alpha ((\partial T)/(\partial t))$ are the sum of the background strain rates, $\dot{\varepsilon} _1$ and $\dot{\varepsilon} _2$, minus the thermal strain rate (extension or contraction) if the ice body is unconstrained. $a = {\dot{\varepsilon}} ^{\prime}_2/{\dot{\varepsilon}} ^{\prime}_1$ is the ratio between the two strain rate components. The expansion coefficient α is assumed to be linearly related to the temperature as α = 5.4 × 10⁻⁵ + 1.8 × 10⁻⁷T (Nishio, Reference Nishio1983) instead of a constant value in Sanderson (Reference Sanderson1978).

The creep factor C depends on the temperature following an Arrhenius relationship: $C = C_*{\rm e}^{-(Q/R)\lpar {(1/T)-(1/T_*)} \rpar }$, where T is the variable temperature that depends on both depth and time, and T _* = 263 K is a transition temperature above which ice softening increases. The gas constant R is equal to 8.314 J mol⁻¹ K⁻¹. The parameter values C _* = 3.5 × 10⁻²⁵ Pa⁻³ s⁻¹ and Q = 1.15 × 10⁵ J mol⁻¹ are based on laboratory experiments and glacial morphologic studies (Cuffey and Paterson, Reference Cuffey and Paterson2010 and reference therein). This C _* value is several orders of magnitude larger than the value used by Sanderson (Reference Sanderson1978), while Q is only slightly larger than Sanderson (Reference Sanderson1978). This updated C _* value incorporates the state-of-the-art lab and field measurements of the glacier creep characteristics, resulting in a lower magnitude of stress in our calculation than Sanderson (Reference Sanderson1978).

We considered two scenarios for the background strain rate: $\dot{\varepsilon} _1 = \dot{\varepsilon} _2 = 0$ and $\dot{\varepsilon} _1$ = 4 × 10⁻¹⁰ s⁻¹ (i.e. 0.013 a⁻¹) and $\dot{\varepsilon} _2$ = −4 × 10⁻¹⁰ s⁻¹. The calculated thermal stress $\sigma _1^* $ has similar diurnal pattern in the top 50 cm of the glacier (Figs 11e–f), but the latter scenario is associated with an additional extensional stress of ~70 kPa at depth due to the background extensional strain rate (Fig. 11f). Our stress calculation for the first scenario shows that the thermal expansion and contraction of the ice leads to diurnal variations in $\sigma _1^* $ from compression to tension with a maximum tensile stress of 48–184 kPa in the top 50 cm (Fig. 11e2). The calculated surface stress becomes tensile after 15 CST as the air temperature starts to decrease (red curves in Fig. 11b). The tensile stress region extends to 20 cm when the hourly short-duration events notably increase at 19 CST. At this time, the surface tensile stress is 156 kPa. After that, the tensile stress region continuously penetrates to greater depth, with reduced amplitude due to depth attenuation of temperature waves. At 8 CST, the air temperature starts to increase and the surface stress becomes compressional. Accordingly, the tensile stress region shrinks. One hour after that, the top 10 cm becomes compressional in stress. This is also the ending time of the daily burst of short-duration events.

For the latter scenario with background strain rate, we assumed a background compressional strain rate $\dot{\varepsilon} _2$ of −4 × 10⁻¹⁰ s⁻¹. This strain rate is based on the surface velocity difference at elevations of ~4400 m and 4500 m measured in autumn and winter by Liu and others (Reference Liu, Qin, Du, Sun and Hou2011). This compressional strain rate is expected to be longitudinal. We further assumed a maximum extensional strain rate $\dot{\varepsilon} _1$ = 4 × 10⁻¹⁰ s⁻¹, presumably in the transverse direction, as the magnitude of the maximum extensional strain is expected to be close to $\dot{\varepsilon} _2$ from previous glacier strain rate measurements (Meier and others, Reference Meier, Kamb, Allen and Sharp1974). The maximum tensional stress for this scenario is the same with the first scenario on the surface, but an increase of 36 kPa at 50 cm depth in comparison with the first scenario is due to the background tensional strain rate. The background strain rate dominantly controls the stress state at one-meter depth or deeper, as the stress asymptotically approaches 73 kPa in this scenario (Fig. 11f2). The enhanced tensile stress at depth can promote downward propagation of surface cracks and continued seismic events.

In either scenario, our thermal stress calculation clearly shows that the daily variations in the glacier thermal state are prominent only in the top tens of centimeters of the glacier. If the thermal effect dominates, it implies that the night-time short-duration events occur or at least initiate within the top tens of centimeters. In addition, the thermal stress is much lower than the stress threshold of 1 MPa that is required to fracture intact ice (Schulson, Reference Schulson1999 and references therein). Instead, the thermal stress is on the same order of magnitude as the critical tensile stress for crack nucleation and propagation if small pre-existing defects or weak planes exist (van der Veen, Reference van der Veen1998).

Comparison with fracture mechanics

Based on the fracture mechanic considerations, the initiation of seismic events implies the moment when the maximum thermal stress reaches the critical tensile stress for fracture formation. The thermal tensile stress during 19 CST to 9 CST the next day with daily bursts of short-duration events is on the order of tens of kilopascals with a maximum tensile stress of 156 kPa on the surface. The calculated thermal tensile stress is consistent with the fracture models of van der Veen (Reference van der Veen1998), who suggests a critical tensile stress range of 90–320 kPa for a densely spaced crack system. The nonexistence of a clear b-value evolution trend implies that the detected short-duration events at different hours of the day likely occur in a similar strain and stress state (Nishio, Reference Nishio1983; Podolskiy and Walter, Reference Podolskiy and Walter2016). As the state of stress and strain differs notably with depth, the short-duration events may not propagate with depth, but mainly correspond to initiation of new surface cracks.

The uncertainties in our stress calculation are primarily due to the applicability of Glen's flow law and a lack of constraints on creep parameter C _* in Glen's flow law. The applicability of Glen's flow law in different temperature and strain rate region is currently under debate (Schulson and Duval, Reference Schulson and Duval2009) and requires further strain and stress measurements. For example, Weiss and others (Reference Weiss, Schulson and Stern2007) found that instead of viscous deformation, the winter and/or perennial sea ice behaves as an elastic-brittle body. In this scenario, the temperature difference between the surface and bottom of the glacier may induce a thermal bending moment at night (Evans and Untersteiner, Reference Evans and Untersteiner1971; Bažant, Reference Bažant1992; Carmichael and others, Reference Carmichael, Pettit, Hoffman, Fountain and Hallet2012), which is later released by microseismic events. For the creep parameter, material properties such as impurities, reduced grain size, increased water content and anisotropic ice fabrics can increase C _* by several fold, and the field measurements of C _* in different glaciers also commonly varies by a factor of 1–3 (Cuffey and Paterson, Reference Cuffey and Paterson2010 and references therein). Material factors are also suggested to influence the fracture toughness (Vaughan, Reference Vaughan1993; Rist and others, Reference Rist1996, Reference Rist1999). Therefore, local field stress and strain measurements are necessary for further calibrating the constitutive relationship and quantitatively understanding icequakes and glacier dynamics.