Waveform arrival picking

Single-station features such as the STA/LTA metric are commonly used to detect impulsive arrivals in continuous seismograms (Allen, Reference Allen1978). However, a shortcoming of the standard STA/LTA approach is that it is sensitive to changing noise levels and generally requires adaptively varying the triggering threshold if the background energy levels are non-stationary, which is common in the glacial environment. In many cases kurtosis-based arrival picking methods appear more robust to statistical fluctuations in the noise level (e.g. Galiana-Merino and others, Reference Galiana-Merino, Rosa-Herranz and Parolai2008; Baillard and others, Reference Baillard, Crawford, Ballu, Hibert and Mangeney2013). Hence, we chose to use a kurtosis phase picker in our work, which was effective at detecting microseismic earthquakes and produced very few false arrival picks, even when energy levels changed dramatically during tremor sequences.

We define our picking algorithm as simply a moving, centered-window estimate of the empirical kurtosis, where at each time step the maximum kurtosis value between the horizontal components is measured. For each ith station we compute

(1a)$$C_i\lpar t\rpar = \max_{\,j \in \lcub 1\comma 2\rcub } {\rm kurtosis}\lpar u_{ij}\lpar t - w\colon t+w\rpar \rpar \comma\; $$

(1b)$${\rm kurtosis}\lpar {\bi x}\rpar = {{1\over n}\sum_{i = 1}^{n}\lpar x_i - \bar{\bi x}\rpar ^4\over \lpar {1\over n}\sum_{i = 1}^{n}\lpar x_i - \bar{\bi x} \rpar ^2\rpar ^2} - 3\comma$$

where u _ij(t) represents the ith station's jth component time-series, the shorthand u _ij(t − w : t + w) represents the centered window around time t of half-width w, x _i represents the ith entry of a vector x and $\bar {{\bi x}}$ is the mean of x. Components are ordered HHE, HHN and HHZ for j = 1, 2 and 3, respectively. We only consider horizontal components in our processing since the vertical channel shows higher noise levels than the horizontal components. The value 3 is subtracted from Eqn (1b) because a normally distributed vector has (unbiased) kurtosis exactly equal to 3, and by subtracting this value the resulting metric measures deviations from a standard normal distribution. Intuitively, when the kurtosis is large the distribution of the data is highly ‘tailed’, implying the presence of many outliers (DeCarlo, Reference DeCarlo1997); a short-duration impulsive arrival surrounded by noise is an example of such a case, however a sustained high energy tremor signal is not. An example detection from using Eqn (1a) is shown in Fig. 2, in which a short-duration impulsive signal (<0.2 s) causes a significant increase in the kurtosis.

Fig. 2.

Example discrete event detection with the kurtosis-based arrival picking method (1a). The window size used in (1a) is w = 1 s, and kurtosis is computed at each 0.2 s increment.

Earthquake location with single stations

The large spatial separation of stations (on average ~8 km), low-magnitude sources and high anelastic attenuation of glacial ice (Peters and others, Reference Peters, Anandakrishnan, Alley and Voigt2012) resulted in a majority of events being detected on only one station at a time, and never more than two stations at a time. A single-station location approach based on waveform polarity and S–P lag times was followed (Magotra and others, Reference Magotra, Ahmed and Chael1989) to estimate locations for events with both discernible compressional (P-) and shear (S-) wave arrivals.

P-wave polarization direction

To implement the location method of Magotra and others (Reference Magotra, Ahmed and Chael1989), the maximum polarization direction of each P wave was used to estimate the propagation direction of the wavefront. This was achieved via eigenvector decomposition of the ground-motion data covariance matrix (following Vidale, Reference Vidale1986; Jurkevics, Reference Jurkevics1988). Specifically, for a three-component slice of data, X, of length n, where X _ji represents the jth component and ith time step of X, the covariance matrix of X is given by

(2)$$\sigma\lpar {\bi X}\rpar = {1\over n} \left(\matrix{\langle {\bi x}_1\comma\; {\bi x}_1\rangle & \langle {\bi x}_1\comma\; {\bi x}_2\rangle & \langle {\bi x}_1\comma\; {\bi x}_3\rangle \cr \langle {\bi x}_2\comma\; {\bi x}_1\rangle & \langle {\bi x}_2\comma\; {\bi x}_2\rangle & \langle {\bi x}_2\comma\; {\bi x}_3\rangle \cr \langle {\bi x}_3\comma\; {\bi x}_1\rangle & \langle {\bi x}_3\comma\; {\bi x}_2\rangle & \langle {\bi x}_3\comma\; {\bi x}_3\rangle }\right)\comma\; $$

where x_j represents the jth row of X, and 〈x_i, x_j〉 is the dot-product between x_i and x_j.

Because Eqn (2) is a symmetric positive semidefinite matrix it has three real, non-negative eigenvalues λ₁ ≥ λ₂ ≥ λ₃ ≥ 0, with corresponding orthonormal eigenvectors v₁, v₂ and v₃ that satisfy σ(X)v_j = λ_jv_j for j=1, 2 and 3. By construction, the eigenvector corresponding to the largest eigenvalue, v₁, approximates the direction in 3-D Cartesian space in which the P wave, which vibrates longitudinally, is propagating. We constrain all v₁ to point downward, by multiplying v₁ by −1 when necessary (which occurs since the sign of the P-wave polarity can vary depending on the source). The incidence angle is estimated with respect to the glacial surface (we approximate the glacier surface as horizontal) for each vector as $\theta _{\rm inc} = \arccos \lpar {\bi v}_{1\comma 3}\rpar$, where v_1,3 is the third component of v₁, and incidence angles of 0° and 90° represent wavefronts propagating perpendicularly and horizontally with respect the glacial surface, respectively.

Hypocentral estimation

Once P-wave propagation directions are estimated with eigenvector decomposition of Eqn (2), event hypocenters are estimated by combining the propagation direction with the inferred source–receiver radial distance. The estimated radial distance is computed from the S–P arrival time difference, t _d. For a homogeneous velocity medium with P- and S-wave velocities V _p and V _s, respectively, the source-to-receiver radial distance Δ and source coordinate h are given by

(3a)$$\Delta = t_{\rm d} {V_{\rm p} V_{\rm s}\over V_{\rm p} - V_{\rm s}}\comma\; $$

(3b)$${\bi h} = {\bi y}_i + \Delta {\bi v}_1\comma\; $$

where y_i is the Cartesian location of the ith station and v₁ is the first eigenvector of the P-wave data covariance matrix (Eqn (2)).

Equations (3a) and (3b) assume that the ice is homogeneous and isotropic. In reality, there are vertical velocity and density gradients due to firn compaction (Vallelonga and others, Reference Vallelonga2014). At NEGIS, there is a ~70 m thick firn layer in the uppermost portions of the glacier (Christianson and others, Reference Christianson2014; Riverman and others, Reference Riverman2019). To account for this we set V _p as varying continuously from 1000 to 3840 m s⁻¹ between 0 and 70 m, and model the continuously bending P-wave with standard finite difference ray tracing (Červenỳ and others, Reference Červenỳ, Molotkov, Molotkov and Pšenčík1977). This allows the trajectory of the P wave to bend away from normal as it goes deeper into the ice and refracts due to the changing velocity structure. The linearly increasing velocity profile we use to account for the firn layer is an approximation to the density structure found in Christianson and others (Reference Christianson2014); however, firn has a broad range of V _p values over many glaciers (Podolskiy and Walter, Reference Podolskiy and Walter2016), and so without further information we select this as a simple first order representation of the increasing velocity profile. We also note that Eqn (3b) uses a simplification, in which free-surface corrections are ignored, however for rays within ~20° of normal incidence this approximation is sufficient (Neuberg and Pointer, Reference Neuberg and Pointer2000).

Lastly, the homogeneous velocity model of this region was somewhat unconstrained at the outset of this study. Christianson and others (Reference Christianson2014) measure an average homogeneous V _p = 3840 m s⁻¹ using data from active seismic surveys, but do not constrain V _s with any additional observations. We follow Christianson and others (Reference Christianson2014) in setting V _p = 3840 m s⁻¹ and compute location estimates over a range of possible V _s values to determine for which range of V _s events locate near the bed. Given our assumptions, the V _s value which localizes events near the bed the most is V _s = 1990 m s⁻¹. This V _s value is very close to the 1930 m s⁻¹ value that is obtained using V _p = 3840 m s⁻¹ and the nominal Poisson ratio of 0.33, and hence does appear to be a reasonable shear wave velocity value of glacial ice.

Clustering arrival waveforms into repeating event groups

Inspection of detected events revealed that many waveforms were nearly identical. In other studies, grouping arrivals into their respective families has been carried out using event templates in parallel with the detection of earthquakes using the template matching (or matched filter) method (Gibbons and Ringdal, Reference Gibbons and Ringdal2006; Allstadt and Malone, Reference Allstadt and Malone2014; Helmstetter and others, Reference Helmstetter, Nicolas, Comon and Gay2015). Without having templates a priori, however, these methods are not directly applicable and a more general approach was required. We used the graph-based Markov clustering algorithm (MCL) (Van Dongen, Reference Van Dongen2000), which is a tool from the wider community detection literature (Lancichinetti and Fortunato, Reference Lancichinetti and Fortunato2009), and which is able to identify non-overlapping group assignments for all events based on pairwise waveform similarity, and hence can be used to extract template waveforms and identify event families.

The MCL method identifies similar populations of events by treating all events as nodes on a mathematical graph and extracting information between the pairwise edge structure of the graph (or adjacency matrix) using an efficient method inspired by random walks on graphs. Between each two nodes in the graph the edge weight is defined as proportional to the maximum cross-correlation coefficient between the waveforms associated with both nodes. That is, the adjacency matrix is given by

(4)$$H_{ij} = c_j\max {1\over 2}\left({{\bi x}_{i\comma 1}\over \Vert {\bi x}_{i\comma 1}\Vert _2} \star {{\bi x}_{\,j\comma 1}\over \Vert {\bi x}_{\,j\comma 1}\Vert _2} + {{\bi x}_{i\comma 2}\over \Vert {\bi x}_{i\comma 2}\Vert _2} \star {{\bi x}_{\,j\comma 2}\over \Vert {\bi x}_{\,j\comma 2}\Vert _2}\right)\comma\; $$

where x_i,1, and x_i,2 denote the ith event's east and north waveform components, respectively, and ⭑ represents the cross-correlation operation. This matrix is calculated separately for each station, and encodes each similar event pair as a high value (maximum of 1), and all dissimilar events as low values (minimum of 0). The matrix (or graph) is additionally sparsified by truncating all entries in H below 0.5 to 0 to ensure that non-zero entries represent actual significant waveform similarities (rather than spurious non-zero correlations). Lastly, in Eqn (4), the column-wise scaling coefficients, c _j, are chosen to ensure the sum over each column is equal to 1 (e.g. $\sum _i H_{ij} = 1$ for $\forall j$). Because the sum over each column is 1, the matrix resembles a stochastic matrix, or in other words a matrix that can be interpreted as encoding discrete transition probabilities between all pairs of nodes, and hence which can be used to simulate random walks on the graph.

The theory behind the MCL method is that by simulating a particular type of oscillatory stochastic flow, which is carried out by iterating two simple matrix operations of ‘inflation’ (given by the pointwise power operation, H^r), and ‘expansion’ (given by matrix–matrix products, H × H), the stochastic flow quickly converges to a sparse, asymmetric matrix with several key properties. In particular, for each node in the resulting matrix there is one outgoing edge, which connects either to itself (i.e. a self-loop), or to a different node in the graph. The nodes with multiple incoming edges from other nodes in the graph are cluster centroids, and generally resemble the most quintessential element of the set of waveforms assigned to a given cluster. Similarly, all of the nodes assigned to a cluster relate to waveforms of high similarity to the other events in the cluster. The subset of nodes which only have self-loop assignments are typically unique or anomalous signals and hence do not appear to be coming from any repeating group.

After the construction of the adjacency matrix there is no other free parameters of the MCL method beyond the single parameter r, which is interpreted to control the granularity or ‘fineness’ of the resulting clusters, where increasing r is generally associated with increasing the number of independent clusters. Since r relates to the pointwise powering operation, H^r, it is usually taken to be small (e.g. within the range r ∈ (1, 10)) to avoid numerical instabilities. To automate the selection of r we scan over a range of r values between (1, 10) and calculate the graph modularity (Newman and Girvan, Reference Newman and Girvan2004), which is a metric that measures how well a clustering assignment naturally conforms to the adjacency matrix of the original graph. Typically, maximizing modularity reveals optimal cluster assignments of a graph (Newman and Girvan, Reference Newman and Girvan2004), however we find that the r value associated with maximum modularity (r _max) has a tendency to create large clusters that combine several similar but distinct families into one. This likely occurs because waveform cross-correlation coefficients of short duration impulsive signals can often be high even for waveforms of different event families. Hence, we opt to set r as r _max + 2, which by our inspection establishes a reasonable balance between revealing a realistic number of families and ensuring clusters do not group many obviously different families into one. The resulting modularity scores as a function of r are shown in Fig. S1, which highlight that the modularity scores obtain a distinct global maxima as r is varied over (1, 10), implying that this procedure for selecting r is stable.

A schematic example of how the clustering algorithm sorts waveforms into common groups from an initially unsorted set is shown in Fig. 3. Here, four earthquake families are present, and MCL reveals the underlying grouping from the initially unsorted pairwise cross-correlation matrix (Eqn (4)). Some strengths of the MCL algorithm are that it is computationally efficient, uses pairwise information between all detections simultaneously, and also automatically infers the number of cluster centroids directly from the data. Nevertheless, many other techniques exist in the field of community detection and may have also been applicable to this problem (for a review, see Lancichinetti and Fortunato, Reference Lancichinetti and Fortunato2009).

Fig. 3.

Schematic example of the waveform clustering method. (a) A small number of initially unsorted waveforms in a dataset. The presence of repeating earthquakes is unknown, potentially present, or not present. (b) Upon clustering with MCL, the resulting edge structure of the new graph assigns each element to a template group (or no group). The new sets of waveforms, once grouped together show high inter-group similarity and low similarity between dissimilar groups. The median stack of the aligned waveforms from each family is taken, to represent the template waveform, as shown in red. For illustration purposes, values in matrices are displayed as binary, however in practice non-zero entries are taken from the continuous interval $[0.5\comma\; 1] \subset {\open R}$.

Continuous tremor detection proxy

A defining characteristic of tremor is that energy levels remain elevated for a long duration (e.g. minutes, hours or days) (Obara, Reference Obara2002; Shelly and others, Reference Shelly, Beroza, Ide and Nakamula2006). As such, an effective proxy for the detection of tremor is to compute statistical measures of the continuous-waveform data distributions over moving time windows (Rouet-Leduc and others, Reference Rouet-Leduc, Hulbert and Johnson2019). If seismic energies are elevated for a sustained duration of time then the interquantile range of the data distribution will increase compared with quiescent intervals. We implement this method and compute the interquantile range of the data, over time t, using a variable window duration w. The metric is define by,

(5)$$C'_i\lpar t\rpar = \max_{\,j\in\lcub 1\comma 2\rcub }[ q_{90}\lpar u_{ij}\lpar t - w\colon t + w\rpar - q_{10}\lpar u_{ij}\lpar t - w\colon t + w\rpar \rpar ] \comma\; $$

which records the maximum value between the interquantile estimates of both horizontal channels. Similar notation is used as in Eqn (1a), and q _p(·) represents the pth quantile of its argument. We remark that the tremor detection proxy (Eqn (5)) is flexible, since short-term and long-term transients can be measured by simply changing w. This metric is not infallible, as impulsive events will also increase C_i′(t), however when the window size of Eqn (5) is much larger than the duration of the transient those contributions are largely suppressed.

Glacial tremor migration

An additional feature of the tremor is that it appears to show slowly migrating moveouts across the seismic network numerous times throughout the study. For example, the tremor proxy curves (Eqn (5)) have pronounced maxima, observed across all of the stations at differing lag times, which are sometimes quite large (e.g. ~300–1200 s). The typical migration speeds of these events are much slower than elastic wave velocities in glacial ice (e.g. ~10's of m s⁻¹ rather than 1000's m s⁻¹), and hence the simplest explanation is that the source of elastic energy is itself migrating across the network.

To parameterize and measure this phenomenon a simple beamforming method was employed to estimate the migration directions and speeds. The tremor source was assumed to migrate across the network as a horizontally oriented plane-wave with a fixed azimuth α and migration speed s. This parameterization is chosen as it is the simplest representation with only two free variables (α and s). Parameter estimates are obtained by minimizing the misfit of such a plane-wave with the observed relative arrival times, τ, of the migrating tremor record across the network, where τ_i is the ith stations relative arrival time of the tremor transient. Using least squares minimization, the solution is given by

(6)$$\min_{\alpha\comma s} \sum_i [ \lpar \tau_i - {\tau}_0\rpar - s^{-1}[ {-}\sin\lpar \alpha\rpar \comma\; \cos\lpar \alpha\rpar ] \cdot\lpar {\bi y}_i - {\bi y}_0\rpar ] ^2\comma\; $$

where · is the dot-product, y_i is the ith stations location and the relative observed and theoretical times are offset by the baseline station, i = 0 (or a different station, if station i = 0 does not have an observed tremor arrival for a given migration event). In a local Cartesian frame where x and y are parallel to the local east and north geographic directions, respectively, the migration direction is given by the vector [−sin(α), cos(α)], where the azimuth α is measured counter-clockwise from geographic north.

In the formulation of Eqn (6), the definition of tremor arrival time (τ) is not exact because the tremor waveforms are typically emergent in nature. Relative arrival times were estimated by visually identifying a characteristic time on the tremor proxy curves (Eqn (5)) that was easily correlated across the network; for example, the maximum amplitude or the onset of the transient was chosen. Observation values were only picked once, and we report source parameters calculated directly from these initial picks (to help minimize bias).

Waveform inversion of tremor

A common interpretation of some types of tectonic tremors is that quasi-continuous signals are emitted from rapidly occurring discrete failures at the source location (Obara, Reference Obara2002; Shelly and others, Reference Shelly, Beroza, Ide and Nakamula2006). When the discrete failures occur at a rate higher than the duration of a typical waveform recording (or support of the Green's function path effect), then the recordings represent an interference of many overlapping waveform arrivals that, if observed individually, would have been impulsive. In accordance, glacial tremor has also been suggested to result from a rapidly occurring set of discrete stick-slip failures at the bed (Winberry and others, Reference Winberry, Anandakrishnan, Wiens and Alley2013; Lipovsky and Dunham, Reference Lipovsky and Dunham2016).

In order to test this model at NEGIS we construct a direct optimization problem where an observed signal d is modeled as the convolution with a fixed Green's function g and an unknown source time function s, by the well-known elastic wave representation theorem (Aki and Richards, Reference Aki and Richards2002, Eqn (3.23)). This formulation is analogous to stacking a fixed template g at a series of (to be determined) time coordinates. The univariate optimization problem is then given by

(7a)$$\min_{{\bi s}} \quad\Vert {\bi d} - {\bi g} \,\ast\, {\bi s} \Vert _{2} + \lambda\Vert {\bi s}\Vert _1$$

(7b)$$\hbox{s.t.} \quad {\bi s} \geq 0\comma\; $$

and the multivariate formulation, with three independent components (i = 1, 2 and 3), by

(8a)$$\min_{{\bi s}} \sum_{i = 1}^{3}\| {\bi d}_i - {\bi g}_i\, \ast\, {\bi s} \|_{2} + \lambda\|{\bi s}\|_1 $$

(8b)$$\hbox{s.t.} \quad {\bi s} \geq 0\comma\;$$

where the * in Eqns (7) and (8) represents the convolution operation. The second terms in Eqns (7a) and (8b) are L1-norm penalties applied to s, weighted by the parameter λ, which is used to promote sparse source activations. Constraints (7b) and (8a) ensures that the source time function is non-negative.

Formulations (7–8) are equivalent to the point source representation of the elastic wave equation (Aki and Richards, Reference Aki and Richards2002, Eqn (3.23)), where the moment tensor is subsumed into the path effect term. This simplification is allowed so long as d is assumed to be produced (or can be approximated) by a single path effect, g, and that the moment tensor is constant for all of the source activations. We discuss this assumption more in the Results section, however note that for repeating earthquake groups this assumption is relatively safe. The source time function is forced to be non-negative in accordance with the assumption that all moment tensors are fixed, and also because empirically we expect more stable and sparse solutions for s when only positive copies of g are forced to stack to reconstruct d.

We minimize Eqn (7) using stochastic gradient descent (Kingma and Ba, Reference Kingma and Ba2014) implemented with automatic differentiation software (Paszke and others, Reference Paszke2017) to iteratively optimize an initial random guess, s₀ ~ Uniform(0, 1), until it converges. We show with synthetic tests that when g is known (or approximately know) this type of inversion is typically stable for a range of λ values (Fig. S3) and accurate reconstructions can be obtained even when event rates are high (e.g. >10 events s⁻¹) and randomly occurring in time. We also give insight into how to evaluate which path effect, g, from a collection of different path effects, {g_i}, is most likely the source of a given tremor recording. In particular, using synthetic tests we show that when all templates are used to reconstruct a synthetic tremor interval the template used to generate the tremor is also recovered as the one with minimum reconstruction error (Fig. S4), indicating that relative reconstruction accuracy can be used to assess which templates are likely the source of a given tremor record.