# Seismological Research Letters

- © 2016 by the Seismological Society of America

## ABSTRACT

We present the software package hybridMT which allows performing seismic moment tensor inversion and refinement, optimized for earthquake data recorded by regional‐to‐local seismic networks as well as for acoustic emission activity. The provided software package is designed predominantly for use in MATLAB (see Data and Resources)/shell environments. The algorithm uses first *P*‐wave amplitudes to invert for unconstrained full, deviatoric, and double‐couple constrained moment tensors. Uncertainty assessment is performed by bootstrap resampling. The moment tensor inversion may be performed directly in the shell environment (by a dedicated command‐line tool) or conveniently through the MATLAB interface (m‐functions). In addition to moment tensor inversion, we also provide the MATLAB implementation of the hybrid moment tensor technique. This methodology increases the quality of calculated seismic moment tensors from events forming a spatial cluster by assessing and correcting for poorly known path and site effects. We tested hybridMT on synthetic datasets, acoustic emission data recorded during laboratory rock deformation experiments, and induced seismicity data from a geothermal reservoir. The package is supplemented with extensive documentation, tutorials, and a dedicated website. HybridMT is freely available and distributed under General Public License.

## INTRODUCTION

The seismic moment tensor (MT) is a standard description of earthquake kinematic source processes over a whole range of magnitudes. Seismic MT inversion allows estimating the fault‐plane parameters and the relation between volumetric and nonvolumetric strain in the seismic source (Knopoff and Randall, 1970). The resolved MTs are typically decomposed into volumetric and deviatoric components, using various decomposition schemes allowing for an understanding of the detailed physical kinematic source processes, regardless of the type of seismicity and event magnitude. The MT inversion has been applied to resolve the displacements in the source for large and small natural earthquakes (Vavryčuk *et al.*, 2008; Scognamiglio *et al.*, 2010; Stierle, Bohnhoff, and Vavryčuk, 2014; Stierle, Vavryčuk, *et al.*, 2014), induced microseismicity (Ross *et al.*, 1996; Panza and Saraò, 2000; Šílený and Milev, 2006; Cesca *et al.*, 2013; Guilhem *et al.*, 2014; Johnson, 2014a,b), as well as for acoustic emission activity measured *in situ* (Manthei *et al.*, 2001; Collins *et al.*, 2002) or in laboratory experiments on rocks samples (Sellers *et al.*, 2003; Thompson *et al.*, 2009; Graham *et al.*, 2010; Charalampidou *et al.*, 2011; Kwiatek, Goebel, and Dresen, 2014). The analysis of seismic MTs sheds light on numerous issues of earthquake physics, such as rupture dynamics (McGarr and Fletcher, 2003; McGarr *et al.*, 2010), fault complexity (McLaskey and Glaser, 2011; McGarr, 2012), the role of pore pressure in seismogenic processes (Fischer and Guest, 2011), and damage‐related radiation of seismic energy (Ben‐Zion and Ampuero, 2009; Castro and Ben‐Zion, 2013, among others).

The MT inversion may be performed using a variety of input data from body‐wave amplitudes to full waveforms. Regardless of the method used, the MT inversion is sensitive to input data quality and modeling assumptions. This includes biases in polarity and amplitude readings for low signal‐to‐noise waveform data, synthetic waveform mismodeling due to poor knowledge on the medium (velocity model, rock anisotropy), and finally site effects and sensor characteristics (coupling, limited frequency band, polarity). All of these unfavorable effects are mapped to the resolved MT (see Sellers *et al.*, 2003, for a comprehensive review). To partially suppress these issues, relative MT inversion techniques have been proposed (Oncescu, 1986; Dahm, 1996; Dahm *et al.*, 1998; Andersen, 2001).

The increasing amount of seismic data recorded worldwide, especially in the field of induced seismicity (Ellsworth, 2013) and the huge amount of data recorded during laboratory experiments (Charalampidou *et al.*, 2011; Davi *et al.*, 2013; Kwiatek, Goebel, and Dresen, 2014) impose the development of an efficient, full‐MT‐inversion software that will be capable of automatically performing the inversion of hundreds or more seismic events, enabling statistical analysis of the derived non‐DC properties. In this study, we present a software package for MT inversion and posterior MT refinement designed for MATLAB (see Data and Resources) and shell environments. The shell application *focimt* and its MATLAB wrapper *focimt.m* provide the means to perform different types of MT inversions, capture the inversion results, and generate the graphical output within or outside the MATLAB environment. The MT inversion results may be optionally improved using provided MATLAB implementation of the hybrid moment tensor (HMT) technique. The HMT technique includes the refinement of seismic MTs of events forming a cluster by suppressing the influence of poorly known source–receiver paths and site effects. We test the hybridMT package using a synthetic dataset, acoustic emission data from laboratory experiments, as well as induced seismicity data from a geothermal site.

## PACKAGE CONTENT

The core of the software package is the stand‐alone shell application *focimt* that is precompiled for a Windows environment. The application is supplemented by two MATLAB subroutines, *focimt.m* and *hybridmt.m*, that wrap input or output MT inversion data handling and perform the HMT refinement. The HTML documentation is accessible from the MATLAB help system or any web browser. In addition, the documentation is provided as probability density function (PDF) file. The latest version of the package, with additional documentation and tutorials, is available through the project website (see Data and Resources). The source codes are freely available at GitHub (see Data and Resources).

*FOCIMT*

*focimt* is an application for MT inversion in time domain and designed for local‐to‐regional seismic networks for which Cartesian coordinates may be used. The core of the program is based on the formal description presented in Fitch *et al.* (1980) and Wiejacz (1992). The provided application is essentially a stand‐alone command‐line tool capable of working in batch mode. In addition, the MATLAB function *focimt.m* is provided to handle both input and output parameters in a convenient way, allowing the performance of the seismic MT inversion directly in the MATLAB environment (Fig. 1).

The input data for the MT inversion include information on amplitudes, rise time, and polarity of the *P*‐wave first motions. The main program input for each sensor is the integral of (i.e., the area below) the first *P*‐wave ground‐displacement pulse which is proportional to the seismic moment. The inversion of first *P*‐wave amplitudes is performed to determine the seismic MT. In addition, the user has to provide the station geometry of the seismic network and rock parameters at the source. Two input data formats are available. Either the required data (take‐off angles, angles of incidence, azimuths, and source–receiver distances for each station) is directly provided through an ASCII file, or those values can be calculated by *focimt*, assuming that the network geometry and 1D velocity model are provided to the internal ray‐tracing routine. Table 1 shows an example of the two described forms of input data for a selected earthquake.

The MT inversion is a well‐known procedure that relies on optimizing the following inverse problem:(1)in which **G** is the *n*‐by‐6 matrix containing Green’s function derivatives, **U** is the *n*‐by‐1 matrix of ground displacements observed, **M** is the 1‐by‐6 matrix containing six independent MT components, and *n* is the number of phase observations. The system of equations (1) is overdetermined and solved for in *focimt*, using a least‐squares approach (L2 norm) with the cost function being the sum of squares of residuals. In addition, *focimt* allows using the absolute (L1) norm, which is less sensitive to larger errors at the cost of increased computation time.

Regardless of the norm used to optimize equation (1), the MT inversion is performed systematically in three different ways, assuming that the target MT is (1) unconstrained, (2) constrained deviatoric, and (3) double‐couple (DC). The deviatoric solution is obtained by assuming no volumetric change in the resolved MT during linear inversion. The DC‐constrained MT is resolved by further imposing the determinant of the seismic MT to be zero. The latter constraint makes the MT inversion scheme nonlinear, and the Lagrange multipliers method (Oncescu, 1986) is used to determine the DC MT.

The resulting seismic MT is decomposed into its isotropic (ISO) and deviatoric part, including a compensated linear vector dipole (CLVD) and a DC, following the convention introduced by Knopoff and Randall (1970): in which and and are the eigenvalues of deviatoric moment tensor with the minimum and maximum absolute values, respectively (e.g., Vavryčuk, 2001). The percentages of the individual MT components can be calculated following two different definitions of the scalar seismic moment (compare equation 13 with equation 14 in Vavryčuk, 2014). The ISO part describes changes in the volume in the seismic source region. The deviatoric part (CLVD and DC) is used to estimate the orientation of tensional, compressional and null‐axis directions (P, T, and B axes), fault‐plane orientation, and the slip vector. Furthermore, the covariance matrix of the MT components is also provided. The faulting type is categorized into either strike slip, normal, or thrust faulting, depending on plunges of P, T, and B axes of the resolved MT.

Uncertainties of the estimated MTs can be estimated through the normalized root‐mean‐square (rms) error between theoretical and estimated amplitudes (Stierle, Bohnhoff, and Vavryčuk, 2014; Stierle, Vavryčuk, *et al.*, 2014) following the formula: (2)

Alternatively, the program can perform an uncertainty assessment of the resulting MTs by randomly sampling input data, assuming that (1) a certain percentage of input phase data display false polarity information, (2) a certain percent of sensors are out of operation, or, (3) the amplitude readings are biased. Having the set of solutions, it is possible to calculate additional parameters allowing an assessment of the quality of resolved MTs and MT‐derived parameters.

The *focimt* output data can be exported to ASCII tables formatted according to user preferences. The program is providing a highly customizable graphical representation of the MT solution in a form of the focal‐mechanism plots (Fig. 1). The output files are either raster (PNG) or vector format (PS, PDF, SVG).

In addition, the provided MATLAB function *focimt.m* conveniently handles the input, processing, and output data in a single routine. Numerous parameters can be specified to maintain the input data, define the seismic MT inversion and optimization parameters, as well as handle the properties of the graphical and text output by MATLAB properties (Fig. 1, see package documentation for details). The output data are stored in MAT files in a form of cell or structure arrays. The auxiliary functions for data reading, synthetic data generation, and output data handling are also provided.

## HYBRIDMT

The HMT technique was originally developed by Andersen (2001). The methodology overall aims at decreasing the influence of local path, site, and sensor effects on estimation of seismic MTs. Frequently, the knowledge on the ray path between source and receiver is very poor. However, the frequency and amplitude content of body waves excited from the source is affected by attenuation, scattering, and local site effects, resulting in a generally unknown Green’s function that can hardly be modeled with typically insufficient information on the geological medium. In addition, misinformation on sensor characteristics directly contributes to the actually measured ground displacement values, leading to an additional bias when resolving seismic MTs.

The HMT technique is performed for selected spatial clusters of seismic events with interevent distances being small, in comparison to the source–receiver distances. This allows the assumption that the earthquakes share similar ray paths and thus the invariance of the Green’s function between events forming a cluster and a particular station. This assumption can then be used to investigate and to account for any persistent bias (bad gain, polarity mismatch) to amplitude readings on each sensor.

In the first step, the MT inversion is performed for all seismic events forming the respective cluster, using the original input amplitude‐phase data. In the following, the resulting seismic MTs are used to predict the amplitudes at each sensor for all events. The ratios are formed for the synthetic and observed amplitudes for each station *i* and event *j*. Then, median ratios are calculated for each sensor separately. The median ratio at station *i* is used to update the input ground‐displacement amplitudes: (3)in which *w*_{i} is the weighting factor. The updated ground‐displacement data are used to calculate a new set of seismic MTs. The procedure is then repeated until the ratio correction factor becomes insignificant (Andersen, 2001), for example .

The HMT technique has been implemented in the routine *hybridmt.m*. The function uses *focimt.m* to recalculate and refine the MTs for seismic events forming a spatial cluster. Similarly to *focimt.m*, numerous parameters can be specified to optimize data preparation, processing, and output and to adjust them to personal needs of the user (see package documentation).

The *hybridmt.m* routine takes essentially the same input data as *focimt.m*; that is, an ASCII file containing event and phase data in one of the two available ASCII formats (compare with Table 1). However, in the HMT technique, it is now explicitly assumed that all events in the input file form a cluster with closely collocated events, in comparison to distances between the cluster and sensors, to satisfy the relative invariance of Green’s function at a particular sensor of all cluster‐member events. In the following step, the input data are analyzed and the trial MT inversion is performed. The numerous trial input and output parameters are presented to the user with information and hints regarding how to optimize the input data (i.e., program notes sensors with insufficient amplitude readings, possible wrong polarities etc.). If the user is satisfied with the input data as well as with trial MT run statistics, the HMT refinement is performed. The HMT results typically consist of refined seismic MTs and extended graphical output.

## APPLICATIONS

### Synthetic Data

To demonstrate the performance of the HMT inversion, we generated a synthetic set of 50 focal mechanisms using the shear‐tensile source model (Vavryčuk, 2001). In that model the source kinematics is described by the four parameters strike (Φ), dip (*δ*), and rake (*λ*) of the fault plane and the tensile angle *α*. The tensile angle describes the angle between the actual slip‐vector direction and the slip vector projected on the fault plane (*α*=−90^{°}, *α*=0^{°}, and *α*=90^{°} correspond to the pure tensile closing, pure shear, and pure tensile opening, respectively). The fault‐plane parameters were randomly sampled with Φ_{i}∈*u*(0°,360°), *δ*_{i}∈*u*(0^{°},90^{°}), and *λ*_{i}∈*u*(−180^{°},180^{°}), in which *u*(·) is the random value drawn from the specified interval). The tensile angle was randomly selected between *α*_{i}=−20° and *α*_{i}=20°, *α*_{i}∈*u*(−20^{°},20^{°}), which corresponds to the seismic MTs containing approximately between −20% and +20% of ISO and CLVD component parts and up to 60% of the remaining DC component. It is important to note that for the shear‐tensile source model the following relation holds: (4)in which *λ* and *μ* are Lamé’s constants and *c*_{ISO} and *c*_{CLVD} are the percentages of ISO and CLVD components calculated according to Vavryčuk (2001). It is also worth noting that when assuming an isotropic medium, the ratio on the left side of equation (4) must be positive; that is, the ISO and CLVD percentages must both be either positive or negative.

We assumed the isotropic medium with constant *P*‐wave velocity. For each synthetic fault plane, the expected *P*‐wave amplitudes were calculated using the radiation pattern formula derived in Ou (2008) and Kwiatek and Ben‐Zion (2013) at 24 sensors located at equal distances from hypocenters covering the whole focal sphere (Fig. 2).

In the following, we intentionally modified the recorded *P*‐wave amplitudes at three sensors, S05, S10, and S15, by introducing constant multiplication factors of 5, 0.5, and 10, respectively. The latter procedure can be perceived as an introduction of a constant factor to the Green’s function observed at certain sensors (e.g., specific attenuation, site amplification, or bad sensor‐gain value). Second, for each seismic event, we randomly removed 40% of synthetic input‐phase data to simulate a lack of *P*‐phase arrivals on some sensors (i.e., each event contained always 15 randomly selected phase‐input data). We assumed that all generated events formed a single cluster, and we used the generated data to perform the MT refinement using an HMT algorithm on the biased catalog.

The comparison of MT inversion results using an initial and a biased catalog, as well as the HMT‐refined catalog, is presented in Figure 3 in an equal‐area Hudson plot (Hudson *et al.*, 1989). The Hudson plot provides a convenient way to present the PDF of source types determined from MTs. The source types of the initial catalog are perfectly aligned as predicted from equation (4). The introduction of a persistent bias to amplitude readings from the three sensors, S05, S10, and S15, resulted in a scatter in the estimated ISO, CLVD, and DC components, as seen in Figure 3a. However, the application of the HMT technique on the biased catalog allows the detection of faulty sensors, and suppressing their effect results in an improved quality of the derived MTs. This is clearly seen by decreasing rms values (Fig. 4a) between iteration 1, corresponding to the biased catalog (compare with Fig. 3a), and iteration 40, corresponding to the final refined MTs (compare with Fig. 3b). The initial bias introduced to amplitudes is well suppressed, as can be seen in Figure 4b.

As shown above, the HMT technique is capable of detecting the persistent bias introduced to a recording station, either in the form of sensor (e.g., wrong gain), site (amplification), or path effects (e.g., velocity or attenuation mismodeling). However, it is worth noting that the technique is not capable of correcting for an eventually wrong station polarity. However, during the course of the HMT refinement, the routine provides statistical information on the polarity match for each station by comparing the observed and predicted polarities. This allows the user to identify sensors with wrong polarities and introduce appropriate corrections to amplitude readings from low‐quality sensors and then to repeat the HMT refinement. The examples of the application of the HMT package in various synthetic scenarios are presented in documentation and on the hybridMT website.

### Laboratory Scale: Acoustic Emissions

In the following, we apply and test the HMT method to acoustic emission (AE) data recorded during a rock deformation laboratory experiment performed on Bentheim sandstone. Resolving MTs for AE is important for understanding the detailed microfracturing processes occurring during the sample loading and failure (e.g., see Kwiatek, Goebel, and Dresen, 2014).

During the laboratory experiment, the AE activity is recorded by a set of 16 piezoceramic sensors (1 MHz resonant frequency, recording range 0.3–0.8 MHz) glued directly to the sample surface (Fig. 5a). This allows for a good coverage of the focal sphere. The waveform data are recorded with a 16‐channel recording system at a 10 MHz sampling rate. The AE data are automatically preprocessed using an offline waveform processing software (Stanchits *et al.*, 2006), including automatic *P*‐onset picking based on the Akaike information criterion (Leonard and Kennett, 1999) and hypocenter determination using a time‐dependent anisotropic 1D‐velocity model (repeatedly updated during the experiment). The resolved hypocenter location uncertainty is typically not exceeding ±2 mm.

It was recently found that the quality of AE recordings can be seriously affected by a poor coupling of the sensors that are glued onto the sample surface or due to damage in direct proximity to the sensor introduced during the experiment (Kwiatek, Charalampidou, and Dresen, 2014). The proposed ultrasonic calibration technique (UCT) that uses ultrasonic pulse data allowed the detection of poorly coupled sensors and provided corrections for amplitudes as well as for the incidence angle of incoming waves (see Davi *et al.*, 2013; Kwiatek, Charalampidou, and Dresen, 2014, for details). In the reported case study, it was found that three of the AE sensors (8, 10, and 13; see Fig. 5) display a reduced coupling quality. The corresponding site correction coefficients derived by UCT are 3.24, 3.37, and 2.43, respectively (rhombs in Fig. 4a, see Kwiatek, Charalampidou, and Dresen, 2014). In the following, we tested whether the HMT technique is also capable of detecting poorly coupled sensors. We selected 192 AE events located in the central part of the sample recorded 4500–4600 s after the beginning of the experiment. The limited time window was chosen to suppress relatively small effects related to sample damage introduced during the course of the experiment due to increased loading. We selected events that occurred within a sphere with a radius of 5 mm (marked with a dashed circle in Fig. 5) located at (*x*,*y*,*z*)=(5,0,50) mm to suppress the potential difference in the propagation effects, allowing the HMT algorithm to recover coupling‐related effects.

Figure 6a presents the comparison of amplitude‐correction factors recovered using MT data and the HMT technique with the corresponding correction factors estimated using the UCT method. One can see that the HMT technique (bars in Fig. 6a) was capable of detecting and suppressing the effects of the sensors with bad coupling. However, the estimated correction factors are generally smaller than those estimated by the UTC method (rhombs in Fig. 6c compared to bars for HMT). In summary, the application of the HMT algorithm resulted in improved estimates for the MTs (Fig. 6b). The refined MTs of AE events display strong negative non‐DC components, an expected result for compaction‐band laboratory experiment (Charalampidou *et al.*, 2011).

### Reservoir Scale: Induced Seismicity

The Geysers geothermal field in northern California is the world’s largest producing geothermal field (Fig. 7). It covers an area of 280 km^{2} and is composed of approximately 60 injection and 330 steam production wells. Production‐related activity (water injection, steam production) led to more than 16,800 induced seismic events (*M*_{D}≥1) recorded by a dense seismic network between 2007 and 2012 (e.g., Eberhart‐Phillips and Oppenheimer, 1984; Majer and Peterson, 2007).

Recent studies on one particular seismicity cluster at the northwestern part of the field analyzed nearly 1800 events between 2007 and 2014 with magnitudes *M*_{w} [1.3, 3.2] focused on the local stress field changes in response to injection (Martínez‐Garzón *et al.*, 2013), physical mechanisms inducing the seismicity at different injection rates (Martínez‐Garzón *et al.*, 2014), and aspects related to the long‐term injection record (Kwiatek *et al.*, 2015).

Here, we selected a subset of 171 induced events from this cluster. Because of the large number of regional seismic stations with varying site conditions, we restricted the MT calculations to use only waveform recordings from the local Berkeley‐Geysers (BG) network jointly operated by the Northern California Earthquake Data Center and the Lawrence Berkeley National Laboratory. This way, up to 32 sensors were available for the MT calculation. To ensure the highest quality of the arrival times, polarities, and input amplitude data, all waveforms were manually reviewed and, when necessary, refined.

In the first step, MTs were estimated using all BG sensors available. Take‐off and incidence angles were estimated using the 1D‐velocity model (Eberhart‐Phillips and Oppenheimer, 1984). MTs were then refined by applying the HMT technique.

Figure 8 shows the improvement of the solutions from the HMT technique with respect to the MT inversion. The Hudson plot from Figure 8a shows that the solutions are scattered, with a preference between the explosive and implosive quadrants. The refined MTs from the HMT technique (Fig. 8b) show a generally increased clustering in the center of the plot (DC earthquakes) as well as a reduction in the earthquakes with CLVD and ISO components of opposite sign. Interestingly, no preferential distribution for the events with a large rms is observed for the initial MT inversion (Fig. 8a), while for the HMT, the largest rms errors are situated in the explosive or implosive quadrants and generally display a relatively large CLVD component (>30%). However, we note that the majority of the most explosive or implosive earthquakes do not have large rms errors, which suggests that the observed non‐DC components are probably not an inversion artifact. Also, many of the largest magnitudes occur close to the DC part of the Hudson plot, and therefore they suggest that larger magnitudes tend to have larger shear components than do the lower magnitudes. Finally, some events still are observed to have ISO and CLVD of different sign. These events may be poorly represented by the velocity model used, exhibit different source complexities, or alternatively be representing other physical properties, such as rock anisotropy.

The HMT quantifies the polarity mismatch between the theoretical and the observed amplitudes. For this particular analyzed dataset, we observed that the average polarity matching per station is between 80% and 90% (Fig. 9a). However, one station (DRH) shows a significantly smaller polarity matching, which may be due to an error of the provided sensor orientation, a low signal‐to‐noise ratio, or proximity of the station to one of the two nodal planes for most events, leading to a polarity mismatch. This station was then removed for obtaining refined MT solutions. In addition, the inversion also provides information on the correction factor for the recorded amplitudes that should be applied to each station. In this sense, we observe that the local station FNF needs to be corrected by approximately a factor of 8 (Fig. 9b), likely the result of inappropriate station gain. In the following, we performed a subsequent HMT refinement, introducing the gain correction for FNF station. For each event, an average of 20 local stations with high‐quality first‐motion amplitudes are used. The obtained results do not indicate a significant improvement of the MT quality, which is likely due to the initial good station coverage.

Most of the seismic events for which the MTs were calculated have an approximate DC component from 60% to 90%, (Fig. 10a). Most of the ISO components oscillate between 0% and 30%, which is in good agreement with the expected slight crack opening from fluid injection activities into a very hot geothermal reservoir. Recently, positive ISO components of on average 20% were reported for the seismicity located in the direct vicinity of the cluster investigated in this study (Johnson, 2014b). Finally, the CLVD components oscillate between −20% and 30%, being in general in accordance in sign with the sign of ISO components. Importantly, the HMT inversion decreased slightly the percentage of events with ISO and CLVD components of opposite sign. We speculate small negative CLVD might be reflecting either crack closure (if ISO is also negative), more complex processes such as anisotropy, or larger uncertainties.

When plotting these events in map view, we observe that events with a large non‐DC component are not evenly distributed within the seismic cluster but in fact seem to be mostly localized in a particular area in the southwestern part of the cluster (Fig. 10b). Interestingly, the clustering of these large non‐DC events does not coincide with the nearest area from the open‐hole section of the well injecting fluid nearby.

## SUMMARY

We provide a new software package for moment tensor inversion and refinement. The software package hybridMT is composed of a stand‐alone shell application *focimt*, a MATLAB wrapper *focimt.m* that simplifies the usage of *focimt* application in the MATLAB environment, and the *hybridmt.m* routine that allows for improving the quality of MT solutions for events forming a spatial cluster by accounting for poorly known path, site, and sensor‐related effects. We use synthetic, acoustic emission, and induced‐seismicity data to demonstrate the feasibility of the software package of performing MT inversion and refinement across different spatial scales. The HybridMT software package delivers high‐quality MT solutions, their graphical representation, and provides a flexible output of calculated MT parameters and their uncertainties. The package is delivered with extensive documentation, and it is freely available under GPL license.

## DATA AND RESOURCES

The Geysers waveform data and metadata were accessed through the Northern California Earthquake Data Center (NCEDC), doi: 10.7932/NCEDC. hybridMT package download, additional documentation and tutorials are available through the project website http://www.gfz-potsdam.de/hybridmt, http://induced.pl/hybridmt (last accessed March 2016). The source code is freely available at GitHub (https://github.com/taquart/hybridmt, last accessed March 2016). MATLAB is available at http://www.mathworks.com/products/matlab/ (last accessed March 2016).

## ACKNOWLEDGMENTS

We acknowledge O. Germer for manually reviewing the waveforms used in this study and D. Oppenheimer for helping with interpretation of BG network operations. P. M. G. acknowledge funding from the Helmholtz Association in the frame of the Helmholtz Postdoc Programme. We would like to thank A. Guilhem and A. Schubnel for providing many helpful comments and suggestions regarding the manuscript.