Abstract
The surveillance of geothermal seismicity is typically conducted using seismic networks, deployed around the power plants and subject to noise conditions in often highly urbanized areas. In contrast, seismic arrays can be situated at greater distances and allow monitoring of different power plants from one central location, less affected by noise interference. However, the effectiveness of arrays to monitor geothermal reservoirs is not well investigated and the increased distance to the source coincides with a decreased accuracy of the earthquake localizations. It is therefore essential to establish robust data processing and to obtain precise estimates of the location uncertainties. Here, we use time-domain array data processing and solve for the full 3-D slowness vector using robust linear regression. The approach implements a Biweight M-estimator, which yields stable parameter estimates and is well suited for real-time applications. We compare its performance to conventional least squares regression and frequency wavenumber analysis. Additionally, we implement a statistical approach based on changepoint analysis to automatically identify P- and S-wave arrivals within the recorded waveforms. The method can be seen as a simplification of autoregressive prediction. The estimated onsets facilitate reliable calculations of epicentral distances. We assess the performance of our methodology by comparison to network localizations for 77 induced earthquakes from the Landau and Insheim deep-geothermal reservoirs, situated in Rhineland-Palatinate, Germany. Our results demonstrate that we can differentiate earthquakes originating from both reservoirs and successfully localize the majority of events within the magnitude range of ML -0.2 to ML 1.3. The discrepancy between the two localization methods is mostly less than 1 km, which falls within the statistical errors. However, a few localizations deviate significantly, which can be attributed to poor observations during the winter of 2021/2022.
1 Introduction
Geothermal energy plays an important role in the transition of the energy sector towards sustainable resources. Unfortunately, high-pressure injection of geothermal fluids is often associated with weak to moderate seismicity (e.g., ; or ). To minimize the seismic hazard, it is crucial to continuously monitor the injection and production processes and localize associated induced seismicity. A reliable and transparent monitoring also helps to increase the public acceptance of existing and future geothermal projects.
Induced seismicity usually relates to man-made stress perturbations of the subsurface, frequently interfering with the local tectonic stress field, and resulting in earthquake activity (e.g., ). It can occur in the context of mining, hydrocarbon or shale gas extraction, wastewater disposal, and geothermal energy production (see, e.g., ; ; ; ). Mechanisms that drive induced seismicity in geothermal environments include pore-pressure and temperature increase, volume change due to fluid withdrawal or injection, and chemical alteration of fracture surfaces (; Zang et al., 2014). Commercial geothermal energy production requires a high geothermal gradient and is therefore often located in active tectonic regions (). The size and rate of seismicity is then defined by the injection volume (and rate), the orientation of the tectonic stress field relative to the pore pressure increase and the extent of the deviatoric stress field within the local fault system (; ; ). analyzes the annual frequency-magnitude distribution of induced geothermal seismicity in central Europe in the period from 2000 to 2011 and compares it to the natural earthquake activity in the region. The results show that induced geothermal events with local magnitudes above ML 2.0 are rare if compared to tectonic earthquakes. However, the intensity of micro-seismicity (ML < 2.0) is significant, with a b-value of 1.94 (±0.21).
To understand the physical processes during a geothermal stimulation and to establish a reliable seismic hazard assessment, the detection and localization of induced geothermal micro-seismicity has gained more and more relevance (). Unfortunately, the preferred locations of geothermal reservoirs lie beneath sedimentary basins, which are often densely populated. This significantly aggravates the detection process due to seismic wave attenuation and an increased level of seismic background noise (; ).
Conventional short-/long-term triggers (STA/LTA) will likely fail under complex noise conditions and advanced detection methods are required (). apply a template correlation trigger to detect micro-seismicity related to a stimulation test in the Landau deep-geothermal reservoir. A similar approach is introduced by , who use the envelope of recorded seismograms to establish a real-time detector based on template correlation. and use (unsupervised) sonogram pattern recognition for event detection. More recent applications focus on deep-learning approaches, such as convolutional neural networks, to monitor induced (micro-) seismicity and to establish automatized seismic phase picking (e.g., Zhu & Beroza, 2019; ; ; ; ). Further methods to obtain automatized seismic phase arrivals include autoregressive (AR) prediction (e.g., ; ), sometimes combined with the Akaike-Information-Criterion (AR-AIC; ; ; ), higher-order statistics (e.g., ) or relative travel-time determination via multi-channel cross correlation (e.g., ).
In contrast to seismic networks, seismic arrays are located outside the source region and can be used to measure the back azimuth and horizontal apparent velocity of an incoming seismic signal, even without clear phase onsets (). Seismic arrays have been frequently used for earthquake detection on a global, regional, and local scale. This includes studies on the Earth’s (fine-scale) structure, detection of human-induced seismicity, volcano monitoring (cf. or and references therein), and ocean-bottom arrays (). Following its initial purpose of detecting nuclear explosions (e.g., ), different studies utilize seismic arrays for seismic risk assessment. , for example, use autoregressive prediction and narrowband f-k analysis for a case study to monitor mining blasts. use a distributed acoustic sensing array and template matching to detect induced geothermal seismicity in the Brady field. Further examples include real-time infrasound monitoring at the Alaska Volcano Observatory () and the real-time array data processing software RETREAT (), developed with a focus on volcano monitoring and volcanic tremor.
Most standard array processing methods apply beamforming in the time- (beam power analysis; see ; ) or frequency-wavenumber domain (f-k analysis; see ). Both approaches perform calculations of the beam power over a predefined slowness grid in the horizontal (-) plane and search for its maximum. Further established methods are progressive multi-channel correlation (), which evaluates a travel time closure condition over narrow frequency bands and for varying combinations of array-station triplets, and incoherent beamforming (; ). and use least squares regression to fit vectors of observed inter-station delay times to obtain estimates for the back azimuth and horizontal apparent velocity of seismic and infrasound signals. The method was adapted by , and . It is computationally efficient (), which makes it suitable for real-time applications. However, least squares regression is sensitive to outliers in the response and predictor variables. Therefore, adopt different robust estimators, including L1-norm regression, weighted M-estimation and Least Trimmed Squares (LTS), and apply them to infrasound data. Their results show significant improvement for different examples with limited data quality, especially for the LTS estimator.
Seismic arrays are less frequently used for distance estimation. For instance, and use manually picked P- and S-wave onsets and a 2-D velocity model to estimate the epicentral distance of events at the Central Indian Ridge and volcanic events near Fogo and Brava, Cabo Verde. They further implement a multi-array analysis, which allows for epicentral localizations without assumptions about the velocity model ().
Our study focuses on developing a computationally efficient and robust solver to determine the slowness vector of seismic phases, with application to induced geothermal seismicity in the Landau and Insheim deep-geothermal reservoirs. We use linear regression to fit observed delay times and, for the first time, implement robust regression estimators for seismic array processing. consider a hypothetically multidimensional array configuration, but practically the method was never applied outside the horizontal plane. We demonstrate that the inclusion of inter-station elevation differences into the regression model yields estimates for the full slowness vector. The regression approaches are subsequently compared to the widely used frequency-wavenumber analysis. We further introduce statistical changepoint analysis as a tool to obtain automatized P- and S-phase arrivals. The approach minimizes the deviation of individual data points from two empirical statistical parameters. This corresponds to a maximization of the likelihood function and can be seen as a simplification of the autoregressive Akaike-Information-Criterion. We evaluate our methodology by a comparison to 77 network localizations from the Landau and Insheim geothermal reservoirs.
2 Study area and array design
The Upper Rhine Graben (URG) is part of the European Cenozoic rift system and is one of Central Europe’s most active tectonic regions with a small to moderate seismic risk (; ). The Landau and Insheim geothermal reservoirs are located near the western rim of the URG in southwestern Germany (cf.Figure 1A). The geological setting includes a crystalline basement, covered by up to 3 km of Paleozoic, Mesozoic and Tertiary sediments and unconsolidated Quaternary sequences (; ). The region has a geothermal gradient of 150 () and features water-bearing sediment layers between two and 3 km depth (). The URG is a densely populated area, and the seismic noise level reaches average values of 10 ground velocity at frequencies of 1–40 ().
FIGURE 1
The Landau and Insheim geothermal power plants are located 4 apart (cf.Figure 1) and are equipped with boreholes of 3,340 and 3,800 depth, respectively (). Both power plants are enhanced geothermal systems with fluid injections in different horizons in the transition zone between Mesozoic sediments and crystalline basement (; ). provide an analysis of the associated seismicity for the years 2006–2013. detect and locate more than 600 events in the Insheim reservoir between 2013 and 2016. In total, more than 2,200 induced micro-events were detected for the period between 2006 and 2016 (). perform re-localizations for events after 2012, using Nonlinloc () and a 1-D velocity model (). We show their results in Figure 1B.
Since 2013, the seismic activity in both reservoirs is monitored by the Federal Institute for Geosciences and Natural Resources (BGR) and the Geological Survey and Mining Authority of Rhineland-Palatinate (LGB-RLP). In addition to the permanent network, a temporary network of seismic stations is operational since 2020 (Südpfalz network). Real-time event detection is implemented using template correlation (); localizations are performed manually by the LGB-RLP ().
The seismic array was installed in June 2021 in the Palatinate Forest, a small mountain range at the western border of the URG, which is characterized by Buntsandstein formations (e.g., ). The distances to the power plants in Landau and Insheim are 12.5 and 14 , respectively (cf.Figure 1A). The array includes 10 seismic stations; it has an aperture of 1.1 (Figure 2) and the maximum distance in elevation between two sites is 71 . The instrumentation includes nine 120 s broadband and one 10 s short-period seismometer. All stations have continuous real-time data streaming of all 3 seismic components (vertical, North, East) at a sampling rate of 200 . The average seismic noise level at frequencies of 5–25 is usually below 0.1 ground velocity.
FIGURE 2
3 Methods
The slowness vector of a seismic phase, traversing an array as a plane wave, relates the inter-station delay time to the position vector (see, e.g., ):
Here, the position vector is defined by the spatial distance between two sites and j () and the inverse of the absolute value of the slowness vector equals the average P- or S-wave velocity of the medium beneath the array (). In Eq. 1.2, the slowness vector is written in terms of the back azimuth angle (BAZ, ) and the horizontal and vertical components of , referred to as the horizontal and vertical apparent velocity ( and ), respectively.
Most array processing techniques estimate the slowness vector in the horizontal plane exclusively. This requires an array setup with marginal differences in elevation. state that elevation correction factors should be applied if deviations in time delay become larger than ¼ of the dominant signal period. These correction terms involve assumptions about the usually unknown subsurface velocity beneath the array () and the vertical incidence angle .
A common procedure to estimate the wavefront parameters and relates to a maximization of the beam power in the horizontal slowness plane . Applying Parseval’s theorem, the beam energy is defined as (e.g., Kelly, 1967):where is the Fourier transform of the signal at site , the position vector of site , the number of sites in the array and the angular frequency. The vectors and are wavenumber vectors, defined by the unknown horizontal slowness vector of the plane wave and the horizontal beam steering vector . The second squared term of the integrand in (2.1) defines the array response function:
The array response function characterizes the array pattern in the wavenumber space at a given frequency. In Figure 2 it is shown for the array in the Palatinate Forest, for a slowness range from −0.4 to 0.4 and a frequency of 10 .
Eq. 2.1 can be evaluated over a grid in the horizontal slowness plane, where the location of the energy maximum provides an estimate for the back azimuth and horizontal apparent velocity of the incident plane wave. The method is referred to as frequency-wavenumber (f-k) analysis.
In our work, we use observed inter-station delay times to estimate the slowness vector from Eq. 1.1 using linear regression (; or ). This requires an accurate estimation of the delay times and a reliable and robust regression approach, which is addressed in detail in the following sections.
3.1 Delay time estimation
We use a normalized cross-correlation function to obtain estimates for the delay times related to an incoming seismic wavefront (e.g., ; ). The normalized cross correlation-function () for two signals and with length and centered at the time , is defined as dependent on the time shift between the signals (e.g., ):with being the expected value, the asterisk (*) denoting the complex conjugate, and and defining the mean and standard deviation of and , respectively. In practice, Eq. 3 can be solved either in the time- or in the frequency-domain.
The travel time difference (delay time, ) of a seismic phase, recorded at two points of observation, is then given by the argument of the maximum () of within a predefined time shift interval of the two signals. Here, the choice of the maximum time shift should account for the aperture of the array and the minimum of the apparent velocity range of interest. In our case we choose , which is appropriate for an aperture of and a minimum of .
The function can be evaluated for all pairs of array stations ( and , with being the number of sites), resulting in a cross-correlation matrix of the seismic signal at the time :
The associated delay time vector forms the basis for our further analysis.
Eqs. 3, 4 are evaluated in fixed time steps, which results in continuous functions of the median cross correlation matrix and the delay time vector (and subsequently back azimuth, and horizontal apparent velocity), with time . In our case, we use a temporal resolution of 10 samples (0.05 ). The values of the time dependent cross-correlation matrix further provide useful parameters for event detection (). We suggest using the median as a robust estimator of the mean of , which must exceed a defined threshold ():
The choice of depends on the local noise conditions, the subsurface characteristics beneath the array sites, and on the size of the correlation window. In our case, proves to be a good trade-off to maintain the sensitivity to earthquakes with small magnitudes while minimizing the number of false detections. It is important to note that at this stage, the term event does not imply a defined source of the signal (e.g., induced/natural earthquake or correlated seismic noise). This is different to, e.g., , who use a continuous correlation with master events to detect seismicity in the reservoirs.
Figure 3 demonstrates the principles of the method for an exemplary event from the Insheim reservoir, recorded on 19 November 2021. It has a local magnitude (ML) of 0.5 and the network localization involves a theoretical back azimuth angle of 97.5° at the array (see Figure 3A). The time series (Figure 3B) are band pass filtered (between 5 Hz and 25 Hz) and a 1.5 time window is used for the correlation analysis. Results are shown for six different site combinations in terms of their normalized cross-correlation functions ( and the corresponding delay times (see Figure 3C), for a time window centered at the point of maximum correlation ). The plane wavefront traverses the array from east to west (BAZ 97.5°), resulting in large delay times for station pairs with significant location differences in east-west direction.
FIGURE 3
The normalized cross-correlation function in Equation 3 performs well for adequate signal to noise conditions, but outliers in single observations can significantly bias subsequent calculations. We therefore recommend using a continuous evaluation of the median of the cross-correlation matrix () as a robust trigger function (bottom Figure 3B). This function (usually) takes a maximum when the correlation window includes the P-wave onsets but excludes the P-wave coda (cf.Figure 3B, dashed red and black lines). It remains unaffected by correlated noise between single station pairs or degraded signal to noise conditions at individual sites. If a set of observed delay times includes outliers, the use of robust array processing techniques will be essential (see section 3.2.1 and supplements).
3.2 Estimating the full 3-D slowness vector using linear regression
and use a linear regression model to fit a vector of inter-station delay times () and obtain estimates for the horizontal slowness vector . The concept can be extended to account for differences in elevation:where is the number of independent observations in the regression model, , and are inter-station distance and elevation difference vectors, includes the unknown error terms and the matrix is the design matrix, defined by the predictor variables. The full 3-D slowness vector corresponds to the coefficient (or parameter) vector of the regression model.
The solution of the regression model is crucial to get accurate parameter estimates. M-estimators (maximum likelihood-type; ) provide a broad class of extremum estimators and allow for the inclusion of robust statistics (see section 3.2.1). They are a generalization of the objective function in L1- (least absolute deviation) and L2-norm (least squares) regression and estimate the maximum of the likelihood function for a parameter vector and a sample distribution with probability density function (, following ). The solution is usually implemented as a minimization of a cost function ():
Here, is a symmetric function of the regression residuals , which has a unique zero (; following ). In Eq. 6, the sample distribution is given by the design matrix and the minimization problem can be written as:
With and being the th row of and the th component of the delay time vector , respectively. The estimated parameter vector minimizes the cost function and predicts the slowness vector of the incident plane wave.
If the errors in are normally distributed, the least squares estimator () recovers the optimum parameter vector that minimizes the squared residuals. The ordinary least squares (OLS) formulation assumes a linear relation between response and predictor variables. It is given as (see, e.g., ):where denotes the transpose of a matrix. Eq. 9 has a unique solution if the matrix is full rank.
The deviation between observed ( and predicted () delay times defines the root mean squared error () of the regression model:with being the number of independent predictor variables (here three). Further, the coefficient variances and covariances are calculated from the mean squared error () of the dependent variable and the inverse covariance matrix of the predictor variables:
The square-root of the diagonal variances are the standard errors of the regression coefficients (, and ). They can be used to derive confidence intervals (CI; e.g., ):where is the percentile of the t-distribution with degrees of freedom. For the 95% confidence interval (significance level ) and 87 degrees of freedom, takes a value of 1.9913 (see any statistical table for the t-distribution).
The back azimuth angle and the horizontal and vertical apparent velocity ( and ) of the plane wave can be calculated from the components of the estimated slowness vector:
We calculate the associated errors using error propagation (see ; ), neglecting the coefficient co-variances, which are in average ten times smaller than the variances. This follows the assumption of uncorrelated errors in the predictor variables.
Figure 4 shows regression results for the exemplary event from the Insheim reservoir, with delay times obtained at the point of maximum correlation (, cf.Figure 3). The regression model reveals a strong dependence in east-west () direction, whereas the dependence in north-south () direction seems much smaller (see Figure 4A). Figure 4B examines the effect of the predictor variables , and on the regression result through partial regression leverage plots (see, e.g., ). It shows that most of the observed delay time (adjusted ) is covered by variations in ; the influences from and are comparable. At the same time, the three regression coefficients (, and ) clearly reject the null-hypothesis at a significance level of 0.05 (), which indicates a substantial contribution to the regression results for all spatial coordinates. However, the standard errors of the regression coefficients reveal that uncertainties within the vertical component of the slowness vector () are by a factor of ten larger if compared to the horizontal components (0.065 for , 0.005 for and ). This is a consequence of the much smaller variance in elevation in comparison to the horizontal inter-station distances (cf. Eq. 11). The estimated 3-D slowness vector and the corresponding 95% confidence ellipsoid are visualized in Figure 4C.
FIGURE 4
Our results show that the inclusion of elevation differences in the regression model allows for an estimation of the vertical slowness component (). However, in our case, the maximum difference in elevation is 71 , which is less than one-fifteenth of the maximum horizontal extension. Therefore, the statistical uncertainties of the results must be evaluated with caution.
3.2.1 Robustness
Outliers in the response and predictor variables can significantly bias ordinary least squares regression and require careful consideration. The statistical definition of an outlier refers to observations that deviate significantly from other members of the underlying data distribution (e.g., or ). In case of seismic arrays, outliers are delay time observations that are inconsistent with the plane wave model (). They can relate to a low signal to noise ratio at one or multiple sites, timing errors (clock drift or failure) or strong subsurface distortion at individual sites.
The effects from outliers on linear regression are well studied and robust estimators are designed to weaken or eliminate their influence (e.g., ). We tested and compared the performance of a Biweight M-estimator, implemented via iteratively reweighted least squares (IRLS), and least trimmed squares (LTS) regression (see Supplementary Text S1, Supplementary Figure S1 and Supplementary Figure S2 for details). The results show that the Biweight M-estimator yields stable and consistent results for BAZ and apparent velocity, even for large quantities of outliers (>25%). IRLS proves to be particularly well suited, as it diminishes effects from outliers by dragging them towards a normal distribution, whilst the parameter estimates are still defined by the mean of the data. Here, the algorithm minimizes the exclusion of data. In this regard it is similar to the limited sensor pair correlation approach (), which improves the robustness of an f-k analysis by excluding weakly correlated sensor pairs.
Our implementation of the Biweight M-estimator follows , and . It adds a weighting term () to the error function () in Equation 7 (e.g.; ):
There is no analytical solution to this equation and we use a reweighted least squares algorithm () to iteratively adjust the weighting function. The algorithm starts from an ordinary least squares regression () and successively reweights observations causing untypically large residuals. For every iteration and each observation , weights () are calculated using the Biweight function (; ):with being the standardized adjusted residuals, calculated from the residuals and the leverage values of the previous iteration:
Here, the leverage values are defined as the diagonal elements of the projection matrix , where defines the diagonal weight matrix. The term is a robust estimate of scale, derived as the median absolute deviation () of the residuals from their median, and with the value related to the inverse of the cumulative distribution function. It makes the estimate unbiased for a normal distribution (see, e.g., ). The factor is a tuning constant, where smaller values result in stronger down-weighting of outliers (we chose ).
The solution of the weighted least squares problem is defined as:
After each iteration, the loss function from Eq. 15 is evaluated:
The algorithm stops if the solution converges or if the iteration limit is reached. We chose the loss function to be a type L2-normalization (), however, it can also contain a least absolute deviation model (L1-norm, ).
3.3 Comparison of methods
In this section, we compare results from iteratively reweighted least squares (IRLS) and ordinary least squares (OLS) regression to the widely used frequency-wavenumber (f-k) analysis (). The performance and stability of an f-k analysis highly depends on the applied frequency band (see, e.g., or ). Generally, the application of suitable, fixed frequency bands is supposed to yield superior results in comparison to a wide-frequency band approach (). At the same time, the width of the adapted frequency bands is crucial and should not be too small (; ).
We calculate the energy of the array beam according to Eq. 2.1 for a slowness-grid in the range from −0.4 to 0.4 and with a resolution of 0.01 s/km. The fixed time window length is 1.5 (identical to the regression approaches). The data are band-pass filtered and tapered prior to the transformation in the frequency domain.
Figure 5 evaluates the performance and stability of the OLS and IRLS regression methods and compares them to frequency-wavenumber analysis. Figure 5A shows the time dependent cross-correlation function for the exemplary event from the Insheim reservoir (ML 0.5, BAZ network localization: 97.5°). Figure 5B compares results for back azimuth and horizontal apparent velocity depending on the shift relative to the point of maximum correlation. It shows that the results are generally unreliable for negative time shifts. In this case, the analysis window does not include enough signal component (cf. Position of the time window at in Figure 3B). However, for positive time shifts, IRLS regression yields reliable and consistent values, not much influenced by the position of the time window. OLS regression is stable up to a + 0.4 s. The f-k results, on the other hand, are unstable and only partially reliable (e.g., at + 0.15 s). Figure 5C examines the frequency band parameters for the f-k analysis with a correlation window shifted by 0.15 from . It shows that the results heavily rely on the frequency band. In this case, the band between 5 and 11 Hz (8 Hz ± 3 Hz) is preferable (BAZ closest to the network localization); however, this conclusion cannot be generalized. The optimal settings can vary significantly depending on the position of the time window and on the source field and ray paths characteristics of the signal (). Figure 5D shows the energy grid of the f-k analysis at , derived for optimized frequency settings (8 Hz ± 3 Hz). The energy maximum is close to the results from IRLS and OLS regression, but there is a nearby secondary maximum, which indicates an instability of the solution. In comparison to the network localization, the back azimuth of the f-k analysis is slightly too small. The deviation most likely relates to the neglection of elevation differences between the array sites. If the regression approaches are restricted to the horizontal plane, the resulting back azimuth is similar to the f-k result.
FIGURE 5
The computation time is of great relevance regarding real-time applications. It is compared for all three approaches using a modern desktop computer (Figure 5E). All methods take less than 0.1 s for one calculation. However, OLS regression is almost 50 times and IRLS approximately 4 times faster if compared to an f-k analysis. Considering that IRLS results are significantly more stable than OLS, the IRLS algorithm appears to be a good trade-off between computational efficiency and accuracy, which makes it an excellent choice for real-time application.
3.4 Distance estimation
The determination of P- and S-wave arrival times is a key task in localizing earthquakes. Many current approaches apply deep learning algorithms, which perform very efficient in real-time applications, but require extensive training data sets (usually 10th of thousands of events, e.g., or ). Another well-established method is the autoregressive Akaike-Information-Criterion (AR-AIC; e.g., ), which uses autoregressive filtering to estimate the wave onset as a maximization of the likelihood function in dependence of the division point between two locally stationary signal segments.
Here, we apply a statistical changepoint approach, which, similar to the AC-AIC, divides a time series signal into two locally stationary segments and evaluates a global statistical parameter for each part (see, e.g., ; or ). The changepoint is then derived through a minimization of a loss function, defined by the residuals from the individual samples with reference to the global parameters.
For decent signal to noise ratios, the onset of a seismic signal usually involves an increase in the signal’s standard deviation. Assuming a time series at site and with length , a function () can be formulated in dependence on the division point :
The signal onset is then defined for the point that minimizes (19)and simultaneously fulfillsi.e., the introduction of a changepoint must improve the cost function. We implement the cost function as a least absolute deviation model (a L1-normalization), which is usually less sensitive to outliers if compared to least squares. Equation 20.1 can be solved by a penalized likelihood approach (see Yao, 1988 or ) and the application of a suitable information criterion (here, Bayesian Information Criterion, BIC). The method can be generalized to search for multiple changepoints using, e.g., binary segmentation (), the segmented neighborhood approach (), or the OP and PELT methods (; ). For a scenario where the approximate travel time difference between P- and S-wave is predictable (e.g., when monitoring reservoirs), it is more efficient to apply two single changepoint searches within appropriate time windows. In general, the search for a statistical changepoint can be extended to the spectral domain (e.g., ), which might yield improved results for limited signal to noise conditions.
Similar to the AR-AIC, the determination of wave onsets using changepoint analysis depends on an initial estimate of the P-wave arrival. Here, the maximum of the median of the cross-correlation function () provides a suitable reference point. For the P-wave arrival we choose a 3.5 time window, starting at . The window for the S-wave arrival starts 0.5 after the estimated P-wave arrival and has a length of 5 . Equations 20.1 and 20.2 are evaluated for all array sites, using the vertical and horizontal components for the P- and S-wave arrivals, respectively. The individual arrival times are subsequently corrected for the inter-station delay times, derived from the array analysis, and a robust estimate for is calculated as:
The median absolute deviation of the individual arrival times from their median is used for error estimation (Hampel, 1986):where the value relates to the inverse of the 0.75th quantile of the cumulative distribution.
Figures 6A, B show continuous evaluations of (Eq. 19) within fixed time windows, for the vertical and East component, respectively. The minima of the functions yield consistent estimates for the P- and S-wave arrivals, related to the exemplary event from the Insheim reservoir.
FIGURE 6
Assuming a homogeneous velocity distribution, the epicentral distance is calculated using a P-wave velocity and a fixed ratio (:
Associated errors are derived using error propagation:
The final localization of an event is defined by the distance and the back azimuth angle , resulting from the array analysis (Eq. 13.1). It can be transformed to Cartesian coordinates, with and . The corresponding errors are:
Latitude and Longitude are calculated with reference to the geometrical mean of the array coordinates.
4 Results for the Insheim and Landau deep-geothermal reservoirs
We compare our results to a data catalogue of 77 induced seismic events from the Insheim and Landau deep-geothermal reservoirs. The catalogue was provided by the Geological Survey and Mining Authority of Rhineland-Palatinate (). It covers a period from July 2021 to May 2022 and includes events in the magnitude range from ML -0.2 to ML 1.3. The events were detected using the Südpfalz network () and a template correlation detector (). Localizations were performed using Seismic Handler () and an optimized minimum 1-D velocity model ().
We re-localize all events using the seismic array in the Palatinate Forest and the methods introduced in section 3. The data from all 10 sites are taken in 20 s windows, starting 5 s ahead of the source times defined by the data catalogue. The waveforms are bandpass filtered between 5 and 25 , and the instrument response is removed. Afterwards, the data are processed in 1.5 windows using the IRLS algorithm. They are continuously shifted by 10 samples (0.05 ), which results in a continuous function of the median cross correlation , the back azimuth θ, and the horizontal apparent velocity with time. Results for each event are obtained for the time window that minimizes the root mean squared error of the linear regression (. Distance estimates are calculated using a constant ratio of 1.76 ± 0.03 and constant P-wave velocities () of 5.15 ± 0.2 and 5.25 ± 0.2 for events from the Landau and Insheim reservoirs, respectively. The values for were determined empirically, by minimizing the deviation from the network localizations. They are in very good agreement with a Granite layer (depth: >3 , : 1.76, : 5.2 ) in the 1-D VSP velocity model for the Insheim reservoir (see ). The homogeneous velocity models are optimized to localize events from the two reservoirs and are not adequate for differing source regions. The results from the array analysis and the corresponding localizations are summarized in Supplementary Tables S1, S2 in the supplements.
Figure 7A shows processing results for the exemplary event from the Insheim reservoir (19 November 2021; ML 0.5). The solutions for back azimuth and horizontal apparent velocity are very consistent during the period of increased correlation, which is associated with the seismic phases traversing the array. The upper plot additionally shows in dependence of the position of the analysis window. Figure 7B compares array and network localization. The distance between both methods is 0.8 , which is well within the statistical errors.
FIGURE 7
Figure 8 compares array and network localizations for the entire data catalogue (Figure 8A). Most array localizations form distinct clusters and can be clearly attributed to either the Insheim or the Landau reservoir (Figure 8B). However, there are some outliers that do not yield reliable results and are far from the network localizations (small map in Figure 8B). Figure 8C gives an overview of the statistics for the array and network localizations. Here, the BAZ and distance of the network localizations are given with reference to the position of the array. The distance estimates resulting from the array analysis are very consistent and usually within a few hundred meters from the network localization. The BAZ values, however, reveal a small and systematic misdirection of +4.1° for the Insheim and – 4.7° for the Landau events.
FIGURE 8
To investigate the quality of the array analysis, Figure 9A visualizes the array localizations, color-coded by the weighted root mean squared error of the regression analysis (cf. equations 10 and 18). The outlying localizations clearly correspond to low-quality linear regression models and involve large standard errors ( and ). This validates the error estimation. For a final comparison between network and array localizations, we correct the systematic misdirection of the back azimuth angle derived from the array analysis (Figure 9B). The resulting median distance between the two methods is 0.9 , with a median absolute deviation () of 0.45 .
FIGURE 9
5 Discussion
Most of the array localizations show a remarkable agreement with the network localizations (median deviation <1 km) and especially the distance estimates are highly consistent. The BAZ calculations feature a systematic misdirection, which can be attributed to either an inadequate velocity model (i.e., 2-D/3-D effects) or local subsurface heterogeneities at the array. The assumption of a uniform velocity distribution is a simplification, and the velocity models are designed to localize seismicity from the two reservoirs. Therefore, the implementation of a more accurate velocity model could further improve the localization accuracy within the reservoirs.
Examining the outliers reveals a direct link between low-quality regression results (i.e., large ) and low signal to noise conditions (Figure 10A). Most significant outliers occur during the period between December 2021 and February 2022, which might relate to an increased level of seismic noise in the northern hemisphere during winter times (see, e.g., ). On the other hand, the corresponding network localizations are also located away from the main clusters of the reservoirs (Figure 10B), either indicating a reduced quality of the network localizations or a slightly different source region within the reservoir. The latter might involve differing source field characteristics, probably hampering an accurate signal recognition at the location of the array. This scenario is supported by a re-evaluation of a low-quality localization from 15 December 2021 (ML 0.7), using an adapted bandpass filter between 25 and 35 (Supplementary Figure S3). The adjusted frequency band yields a more reliable regression result (though the errors remain large), which suggests wavefield characteristics that differ from a typical event from the reservoirs. In such a case, advanced methods in the spectral domain could yield improved results. , for example, analyze the spatial coherence of the seismic wavefield by an eigenvalue decomposition of the covariance matrix. Incoherent noise is then minimized through a reduction of the signal to components related to the dominant eigenvalue. The PMCC algorithm () calculates cross-correlation functions within narrow frequency bands, thus offering an improved separation between frequency bands with and without noise. On the downside, those methods are either computationally more expensive or require manual adjustments to the signal and noise characteristics.
FIGURE 10
The statistical errors for the array and network localizations consistently increase with decreasing localization quality (Figure 10C). This is a good validation of the error calculation, which is crucial for the quality assessment during real-time processing. The standard errors associated with low-quality array localizations, derived in the winter of 2021/2022, are exceptionally large and coincide with large errors for the corresponding network localizations. It is important to note that the errors from the array analysis are generally smaller when compared to errors resulting for the network localizations. This might partly be due to the different calculation approaches, but more likely it reflects the superior location characteristics at the array. Here, the seismic noise level is in average 0.1 (at 1–25 ), which is about 100 times smaller if compared to the Upper Rhine Graben (10 at 1–40 ; ). Further, the network analysis and associated statistical errors are likely distorted by local 3-D velocity anomalies. For the array approach, the similarity of the ray paths between source and receivers involves smaller statistical errors (provided the data quality at the array sites is good). In this case, the expected errors have a systematic origin and relate to the uncertainty of the adapted uniform velocity model.
Figure 10D examines the distribution of the seismicity within the two reservoirs, separately for array and network localizations. It shows that the array results cluster more distinctly, especially for the Landau reservoir. Here, the events focus between injection and production well. The corresponding network localizations also mainly occur between injection and production side, but they scatter more widely. In the case of the Insheim reservoir, seismicity is concentrated southwest of the injection wells for both localization methods. Again, the horizontal variation is smaller for the array localizations. At this point it is difficult to conclude which results are more accurate. Looking at source-receiver distances exclusively, the 1-D velocity model used in the network analysis is not superior to an optimized uniform model. However, the missing depth dependence in the array approach also affects the epicentral localization.
Our regression approach includes inter-station elevation differences, which allows for estimates of the full 3-D slowness vector. In Figure 11 we investigate results for the horizontal and vertical apparent velocities ( and ), and the vertical angle of incidence , for the entire data catalogue. The horizontal apparent velocity is mostly between 6 and 8 (median: 6.6 , MAD: 0.4 ) and the vertical apparent velocity between 3 and 7 (median: 4.1 , MAD: 1 ). The angles of incidence mostly range between 20° and 40°. For the Insheim reservoir, some events feature increased vertical apparent velocity values. The associated angles of incidence are larger than 50°. Low-quality regression results for the period between December 2021 and February 2022 involve very small values for the vertical apparent velocity and the vertical angle of incidence. Consistent with the results for the exemplary event (Section 3.2), errors associated with are distinctly larger if compared to the horizontal components. To estimate the vertical slowness component with sufficient accuracy, differences in elevation should be of the same order of magnitude as the horizontal distances between the array sites. In this case, the vertical angle of incidence could be used for depth estimation, e.g., using a vertical (borehole) array in combination with a conventional horizontal layout.
FIGURE 11
Our analysis shows that conventional array processing techniques, such as f-k analysis or OLS regression, are highly sensitive to outlying data points and can heavily rely on manual adjustments of the evaluation parameters (window size and position, frequency band). Here, the IRLS algorithm, in combination with the Biweight function, must be preferred. It yields stable and consistent results, even in the presence of corrupted data. In comparison to a regular f-k analysis, it is less sensitive to the frequency band (but requires SNR >1 for the P wave onset) and is computationally more efficient. The application of the algorithm to real-time data involves a continuous evaluation of the cross-correlation function and continuous robust estimations of the slowness vector and associated uncertainties. If the correlation function exceeds a certain threshold, distance estimates are calculated. In case of insufficient processing results, indicated by the statistical errors, the event can be re-evaluated using an adapted frequency band. Here, an automatized choice of the filter settings can be based on the cross-spectral matrix of the signals. The computational efficiency of the algorithm would also allow for a second evaluation stream (i.e., secondary continuous calculations of the cross-correlation, the slowness vector and associated uncertainties) within a different frequency band.
6 Conclusion
We investigate the suitability of seismic arrays for monitoring multiple geothermal reservoirs from one central remote location. Here, the increased distance to the source requires accurate processing techniques to receive reliable earthquake localizations. We therefore employ robust linear regression to estimate the slowness vector of seismic phases and use statistical changepoint analysis to obtain automatized P- and S-wave arrival times, which can be used for distance calculations. The comparison to standard array processing tools, such as ordinary least-squares regression and f-k analysis, demonstrates that a robust approach is crucial to achieve localization accuracy suitable for geothermal monitoring. We further validate our results using a data catalogue of 77 network localizations for the Landau and Insheim deep-geothermal reservoirs, located in the Upper Rhine Graben. It shows that we can clearly separate earthquakes originating from the two reservoirs and the quality of the array localizations is at least comparable to those from the seismic network. Moreover, the remote location of the array involves a significantly lower level of seismic noise compared to the seismic network. This enhances the sensitivity to small magnitude events and ensures surveillance during noisy episodes.
Estimating the slowness vector of a seismic phase using linear regression relies on observed delay times, derived from inter-station cross-correlation functions. Here, we recommend using the median of the cross-correlation matrix as a robust trigger function as it remains unaffected by correlated noise between single station pairs or degraded signal to noise conditions at specific array sites.
We further demonstrate that incorporating elevation differences into the regression model allows for an estimation of the vertical slowness component (). This separates the method from those limited to the horizontal plane, but the statistical significance and accuracy of the results must be evaluated with caution. If an array is to be used to accurately estimate the vertical apparent velocity (and subsequently the vertical angle of incidence), the elevation differences between the array sites should be of the same order of magnitude as the horizontal inter-station distances. However, elevation differences should be included in the array analysis if they are expected to contribute to the observed delay times. This eliminates the need for elevation correction terms, and the impact on the calculation time is only marginal.
When a set of observed delay times includes outliers, the use of robust array processing techniques is crucial. We therefore implement and test robust regression estimators for seismic array data. Here, iteratively reweighted least squares in combination with a Biweight function yields reliable parameter estimates, that are significantly more stable compared to conventional least squares regression and f-k analysis. The algorithm is computationally efficient, making it well suited for real-time applications.
To obtain P- and S-wave arrivals by an automated approach, we introduce statistical changepoint analysis as an alternative to autoregressive prediction. The determination of a statistical changepoint only relies on the calculation of a basic statistical parameter and not on autoregressive filtering. This makes it computationally more efficient, at least when the search problem is restricted to a single changepoint. The quality of the estimated arrival times is remarkable, resulting in highly accurate distance estimates for the array localizations.
The final comparison between network and array localizations shows that the results are very consistent. Most array localizations form distinct clusters that can be clearly attributed to either the Insheim or the Landau reservoir. A few outliers for the array localizations in the period between December 2021 and February 2022, coincide with low-quality network localizations in the outer reservoir domains. Upon closer examination of the seismicity within the two reservoirs, it becomes evident that the array results cluster more distinctly than the network results. Furthermore, the statistical errors from the array analysis are generally smaller compared to those from the network localizations. This reflects the superior location characteristics at the array, where the average seismic noise level is about 100 times lower than in the Upper Rhine Graben. As a result, the quality of the epicentral array localizations is at least comparable with those derived from the network.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
PH was responsible for the planning and setup of the seismic array. He developed the methodology and the software used to obtain the results and wrote the initial draft of the manuscript. ML was involved in the planning and realisation of the project. He participated in the determination of the array layout and setup and contributed to the final draft of the manuscript. GR conceived and initiated the project, taking responsibility for its execution. All authors contributed to the article and approved the submitted version.
Funding
This work was funded through the Federal Ministry for Economic Affairs and Climate Action of the Federal Republic of Germany as part of the SEIGER (Seismic monitoring of deep geothermal power plants and possible seismic impact) research project (FKZ 03EE4003F).
Acknowledgments
We thank the Geological Survey and Mining Authority of Rhineland-Palatinate (LGB-RLP), and particularly Bernd Schmidt and Helmuth Winter, for their support during the setup of the array and for providing the earthquake catalogue. We thank KH, FN, VT, LDS for their careful and constructive comments, which improved the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1217587/full#supplementary-material
References
1
AkaikeH. (1973). Maximum likelihood identification of Gaussian autoregressive moving average models. Biometrika60 (2), 255–265. 10.1093/biomet/60.2.255
2
AugerI. E.LawrenceC. E. (1989). Algorithms for the optimal identification of segment neighborhoods. Bull. Math. Biol.51 (1), 39–54. 10.1016/S0092-8240(89)80047-3
3
BartzJ. (1974). “Die Mächtigkeit des Quartärs im Oberrheingraben,” in Approaches to taphrogenesis. Editors IlliesJ. H.FuchsK. (Stuttgart, Germany: Schweitzerbarth), 78–87.
4
BeatonA. E.TukeyJ. W. (1974). The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. Technometrics16 (2), 147–185. 10.1080/00401706.1974.10489171
5
BGR (2023). Bundesanstalt für Geowissenschaften und Rohstoffe. https://www.bgr.bund.de/DE/Themen/Erdbeben-Gefaehrdungsanalysen/Seismologie/Seismologie/GERSEIS/gerseis_node.html (Accessed June 02, 2023).
6
BishopJ. W.FeeD.SzuberlaC. A. (2020). Improved infrasound array processing with robust estimators. Geophys. J. Int.221 (3), 2058–2074. 10.1093/gji/ggaa110
7
BruneJ.ThatcherW. (2002). “Strength and energetics of active fault zones,” International handbook of earthquake and engineering seismology, Part A. Editors LeeW. H. K.KanamoriH.PaulC. J.KisslingerH. (Int. Assoc. Seismol. Phys. Earth's Interior, Comm. Educ), 81, 569–588.
8
CansiY. (1995). An automatic seismic event processing for detection and location: The PMCC method. Geophys. Res. Lett.22 (9), 1021–1024. 10.1029/95GL00468
9
CaponJ. (1969). High-resolution frequency-wavenumber spectrum analysis. Proc. IEEE57 (8), 1408–1418. 10.1109/PROC.1969.7278
10
ChenJ.GuptaA. K. (1997). Testing and locating variance changepoints with application to stock prices. J. Am. Stat. Assoc.92 (438), 739–747. 10.1080/01621459.1997.10474026
11
ChenJ.GuptaA. K. (2001). On change point detection and estimation. Commun. statistics-simulation Comput.30 (3), 665–697. 10.1081/SAC-100105085
12
ClaerboutJ. F. (1986). Fundamentals of geophysical data processing, 2nd edn. J. R. Astron. Soc.86, 217–219. 10.1111/j.1365-246X.1986.tb01085.x86
13
CoombsM. L.WechA. G.HaneyM. M.LyonsJ. J.SchneiderD. J.SchwaigerH. F.et al (2018). Short-term forecasting and detection of explosions during the 2016–2017 eruption of Bogoslof volcano, Alaska. Front. Earth Sci.6, 122. 10.3389/feart.2018.00122
14
CornetF. H.HelmJ.PoitrenaudH.EtchecoparA. (1997). Seismic and aseismic slips induced by large-scale fluid injections. Pure Appl. Geophys.150, 563–583. 10.1007/978-3-0348-8814-1_12
15
CornetF. H.JianminY. (1995). Analysis of induced seismicity for stress field determination and pore pressure mapping. PAGEOPH145, 677–700. 10.1007/BF00879595
16
CornetF. H.JulienPh. (1989). Stress determination from hydraulic test data and focal mechanisms of induced seismicity. Int. J. Rock Mech. Min. Sci. Geomechanics Abstr.26 (3-4), 235–248. 10.1016/0148-9062(89)91973-6
17
CuenotN.DorbathC.DorbathL. (2008). Analysis of the microseismicity induced by fluid injections at the EGS site of soul-sous-forêts (alae, France): Implications for the characterization of the geothermal reservoir properties. Pure Appl. Geophys.165, 797–828. 10.1007/s00024-008-0335-7
18
De AngelisS.HaneyM. M.LyonsJ. J.WechA.FeeD.Diaz-MorenoA.et al (2020). Uncertainty in detection of volcanic activity using infrasound arrays: Examples from Mt. Etna, Italy. Front. Earth Sci.8, 169. 10.3389/feart.2020.00169
19
Del PezzoE.GiudicepietroF. (2002). Plane wave fitting method for a plane, small aperture, short period seismic array: A MATHCAD program. Comput. geosciences28 (1), 59–64. 10.1016/S0098-3004(01)00076-0
20
DoeblF.OlbrechtW. (1974). “An isobath map of the Tertiary base in the Rhinegraben,” in Approaches to taphrogenesis. Editors IlliesJ. H.FuchsK. (Stuttgart, Germany: Schweitzerbarth), 71–72.
21
DornstadterJ.KappelmeyerO.WelterM. (1999). “The geothermal potential in the Upper Rhine Graben valley,” in Proceedings of the European Geothermal Conference Basel '99 (Basel, Switzerland) 2, 77–85.
22
DouglasA.BowersD.MarshallP. D.YoungJ. B.PorterD.WallisN. J. (1999). Putting nuclear-test monitoring to the test. Nature398 (6727), 474–475. 10.1038/19000
23
DumouchelW.O’brienF. (1989). “Integrating a robust option into a multiple regression computing environment,” in Computer science and statistics: Proceedings of the 21st symposium on the interface (Alexandria, Egypt: American Statistical Association), 297–302.
24
EvansK. F.ZapponeA.KraftT.DeichmannN.MoiaF. (2012). A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe. Geothermics41, 30–54. 10.1016/j.geothermics.2011.08.002
25
FarahbodA. M.KaoH.CassidyJ. F.WalkerD. (2015). How did hydraulic-fracturing operations in the Horn River Basin change seismicity patterns in northeastern British Columbia, Canada?Lead. Edge34 (6), 658–663. 10.1190/tle34060658.1,
26
GibbonsS. J.KværnaT.RingdalF. (2010). Considerations in phase estimation and event location using small-aperture regional seismic arrays. Pure Appl. Geophys.167, 381–399. 10.1007/s00024-009-0024-1
27
GibbonsS. J.KværnaT.RingdalF. (2005). Monitoring of seismic events from a specific source region using a single regional array: A case study. J. Seismol.9, 277–294. 10.1007/s10950-005-5746-7
28
GibbonsS. J.NäsholmS. P.RuigrokE.KværnaT. (2018). Improving slowness estimate stability and visualization using limited sensor pair correlation on seismic arrays. Geophys. J. Int.213 (1), 447–460. 10.1093/gji/ggx550
29
GibbonsS. J.RingdalF.KværnaT. (2008). Detection and characterization of seismic phases using continuous spectral estimation on incoherent and partially coherent arrays. Geophys. J. Int.172 (1), 405–421. 10.1111/j.1365-246X.2007.03650.x
30
GroosJ. C.RitterJ. R. (2014). Verbundprojekt MAGS - konzepte zur Begrenzung der mikroseismischen Aktivität bei der energetischen Nutzung geothermischer Systeme im tiefen Untergrund, Einzelprojekt 1: Quantifizierung und Charakterisierung des induzierten seismischen Volumens im Bereich Landau/Südpfalz. Germany: Karlsruher Institut für Technologie, Final Report of the Project MAGS.
31
GrubbsF. E. (1969). Procedures for detecting outlying observations in samples. Technometrics11 (1), 1–21. 10.1080/00401706.1969.10490657
32
GrünthalG. (2014). Induced seismicity related to geothermal projects versus natural tectonic earthquakes and other types of induced seismic events in Central Europe. Geothermics52, 22–35. 10.1016/j.geothermics.2013.09.009
33
GrünthalG.WahlströmR. (2012). The European-Mediterranean earthquake catalogue (EMEC) for the last millennium. J. Seismol.16, 535–570. 10.1007/s10950-012-9302-y
34
HanekeJ.WeidenfellerM. (2010). “Die geologischen baueinheiten der Pfalz,” in Geographie der Pfalz (Landau/Pfalz, Germany: Pfälzische Landeskunde).
35
HaneyM. M.Van EatonA. R.LyonsJ. J.KramerR. L.FeeD.IezziA. M. (2018). Volcanic thunder from explosive eruptions at Bogoslof volcano, Alaska. Geophys. Res. Lett.45 (8), 3429–3435. 10.1002/2017GL076911
36
HollandP. W.WelschR. E. (1977). Robust regression using iteratively reweighted least-squares. Commun. Statistics-theory Methods6 (9), 813–827. 10.1080/03610927708827533
37
HuberP. J. (1973). Robust regression: Asymptotics conjectures and Monte Carlo. Ann. Statistics1 (5), 799–821. 10.1214/aos/1176342503
38
HuberP. J. (1981). Robust statistics. New York, NY, USA: John Wiley and Sons, Inc.
39
HurtigE.CermakV.HaenelR.ZuiV. (1992). Geothermal atlas of Europe. Vienna, Austria: International Atomic Energy Agency.
40
IlliesJ. H. (1972). The Rhine graben rift system-plate tectonics and transform faulting. Geophys. Surv.1 (1), 27–60. 10.1007/BF01449550
41
JacksonB.ScargleJ. D.BarnesD.ArabhiS.AltA.GioumousisP.et al (2005). An algorithm for optimal partitioning of data on an interval. IEEE Signal Process. Lett.12 (2), 105–108. 10.1109/LSP.2001.838216
42
JohnsonS. W.ChambersD. J.BoltzM. S.KoperK. D. (2021). Application of a convolutional neural network for seismic phase picking of mining-induced seismicity. Geophys. J. Int.224 (1), 230–240. 10.1093/gji/ggaa449
43
JoswigM. (1990). Pattern recognition for earthquake detection. Bull. Seismol. Soc. Am.80 (1), 170–186. 10.1785/BSSA0800010170
44
KillickR.FearnheadP.EckleyI. A. (2012). Optimal detection of changepoints with a linear computational cost. J. Am. Stat. Assoc.107 (500), 1590–1598. 10.1080/01621459.2012.737745
45
KingD. W.HaddonR. A. W.HusebyeE. S. (1975). Precursors to PP. Phys. Earth Planet. Interiors10 (2), 103–127. 10.1016/0031-9201(75)90029-1
46
KingD. W.HusebyeE. S.HaddonR. A. W. (1976). Processing of seismic precursor data. Phys. Earth Planet. Interiors12 (2-3), 128–134. 10.1016/0031-9201(76)90042-X
47
KrügerF.DahmT.HannemannK. (2020). Mapping of Eastern North Atlantic Ocean seismicity from Po/So observations at a mid-aperture seismological broad-band deep sea array. Geophys. J. Int.221 (2), 1055–1080. 10.1093/gji/ggaa054
48
KüperkochL.MeierT.BrüstleA.LeeJ.FriederichW. (2012). Automated determination of S-phase arrival times using autoregressive prediction: Application to local and regional distances. Geophys. J. Int.188 (2), 687–702. 10.1111/j.1365-246X.2011.05292.x
49
KüperkochL.MeierT.LeeJ.FriederichW.Egelados Working Group (2010). Automated determination of P-phase arrival times at regional and local distances using higher order statistics. Geophys. J. Int.181 (2), 1159–1170. 10.1111/j.1365-246X.2010.04570.x
50
KüperkochL.OlbertK.MeierT. (2018). Long‐term monitoring of induced seismicity at the Insheim geothermal site, Germany. Bull. Seismol. Soc. Am.108 (6), 3668–3683. 10.1785/0120170365
51
KværnaT.DoornbosD. J. (1991). Scattering of regionalP n by moho topography. Geophys. Res. Lett.18 (7), 1273–1276. 10.1029/91GL01292
52
KværnaT.RingdalF. (1986). NORSAR scientific rep 1-86/87. Stab. Var. fk Estim. Tech., 29–40.
53
LaiT. L.RobbinsH.WeiC. Z. (1978). Strong consistency of least squares estimates in multiple regression. PNAS75 (7), 3034–3036. 10.1073/pnas.75.7.3034
54
LeonardM.KennettB. L. N. (1999). Multi-component autoregressive techniques for the analysis of seismograms. Phys. Earth Planet. Interiors113 (1-4), 247–263. 10.1016/S0031-9201(99)00054-0
55
LevaC.RümpkerG.WölbernI. (2022). Multi-array analysis of volcano-seismic signals at Fogo and Brava, Cape Verde. Solid earth.13, 1243–1258. 10.5194/se-13-1243-2022
56
LevaC.RümpkerG.WölbernI. (2020). Remote monitoring of seismic swarms and the August 2016 seismic crisis of Brava, Cabo Verde, using array methods. Nat. Hazards Earth Syst. Sci.20 (12), 3627–3638. 10.5194/nhess-20-3627-2020
57
LGB-RLP (2022). Landesamt für Geologie und Bergbau Rheinland-Pfalz. https://www.lgb-rlp.de/de/fachthemen/landeserdbebendienstrlp/forschungsprojekt-seiger.html (Accessed February 08, 2023).
58
LiW.ChakrabortyM.FennerD.FaberJ.ZhouK.RümpkerG.et al (2022). EPick: Attention-based multi-scale UNet for earthquake detection and seismic phase picking. Front. Earth Sci.10. 10.3389/feart.2022.953007
59
LiZ.ZhanZ. (2018). Pushing the limit of earthquake detection with distributed acoustic sensing and template matching: A case study at the Brady geothermal field. Geophys. J. Int.215 (3), 1583–1593. 10.1093/gji/ggy359
60
LomaxA.VirieuxJ.VolantP.Berge-ThierryC. (2000). “Probabilistic earthquake location in 3D and layered models,” in Advances in seismic event location. Modern approaches in geophysics. Editors ThurberC.RabinowitzN. (Dordrecht, Netherlands: Springer), 101–134. 10.1007/978-94-015-9536-0_5
61
MajerE. L.BariaR.StarkM.OatesS.BommerJ.SmithB.et al (2007). Induced seismicity associated with enhanced geothermal systems. Geothermics36 (3), 185–222. 10.1016/j.geothermics.2007.03.003
62
MousaviS. M.EllsworthW. L.ZhuW.ChuangL. Y.BerozaG. C. (2020). Earthquake transformer—An attentive deep-learning model for simultaneous earthquake detection and phase picking. Nat. Commun.11, 3952. 10.1038/s41467-020-17591-w
63
OlsonJ. V.SzuberlaC. A. (2005). Distribution of wave packet sizes in microbarom wave trains observed in Alaska. J. Acoust. Soc. Am.117 (3), 1032–1037. 10.1121/1.1854651
64
PicardD. (1985). Testing and estimating change-points in time series. Adv. Appl. Probab.17 (4), 841–867. 10.2307/1427090
65
PlenkersK.RitterJ. R. R.SchindlerM. (2013). Low signal-to-noise event detection based on waveform stacking and cross-correlation: Application to a stimulation experiment. J. Seismol.17, 27–49. 10.1007/s10950-012-9284-9
66
RitterJ. R.SudhausH. (2007). Characterization of small local noise sources with array seismology. Near Surf. Geophys.5 (4), 253–261. 10.3997/1873-0604.2007007
67
RostS.ThomasC. (2002). Array seismology: Methods and applications. Rev. Geophys.40 (3), 2–27. 10.1029/2000RG000100
68
RousseeuwP. J. (1984). Least median of squares regression. J. Am. Stat. Assoc.79 (388), 871–880. 10.1080/01621459.1984.10477105
69
RousseeuwP. J.LeroyA. M. (1987). Robust regression and outlier detection. Hoboken, New Jersey, United States: John Wiley and Sons, Inc.
70
SchweitzerJ.FyenJ.MykkeltveitS.GibbonsS. J.PirliM.KühnD.et al (2012). Seismic arrays. New Man. Seismol. observatory Pract. 2 (NMSOP-2), 1–80. 10.2312/GFZ.NMSOP-2_ch9
71
SenA.SrivastavaM. S. (1975). On tests for detecting change in mean. Ann. Statistics3, 98–108. 10.1214/aos/1176343001
72
SeydouxL.ShapiroN. M.De RosnyJ.LandèsM. (2016). Spatial coherence of the seismic wavefield continuously recorded by the USArray. Geophys. Res. Lett.43 (18), 9644–9652. 10.1002/2016GL070320
73
ShiX.GallagherC.LundR.KillickR. (2022). A comparison of single and multiple changepoint techniques for time series data. Comput. Statistics Data Analysis170, 107433. 10.1016/j.csda.2022.107433
74
SickB.GuggenmosM.JoswigM. (2015). Chances and limits of single-station seismic event clustering by unsupervised pattern recognition. Geophys. J. Int.201 (3), 1801–1813. 10.1093/gji/ggv126
75
SinghM.RümpkerG. (2020). Seismic gaps and intraplate seismicity around Rodrigues Ridge (Indian Ocean) from time domain array analysis. Solid earth.11 (6), 2557–2568. 10.5194/se-11-2557-2020
76
SleemanR.Van EckT. (1999). Robust automatic P-phase picking: An on-line implementation in the analysis of broadband seismogram recordings. Phys. earth Planet. interiors113 (1-4), 265–275. 10.1016/S0031-9201(99)00007-2
77
SmithP. J.BeanC. J. (2020). Retreat: A REal-time TREmor analysis tool for seismic arrays, with applications for volcano monitoring. Front. Earth Sci.8, 586955. 10.3389/feart.2020.586955
78
StammlerK. (1993). SeismicHandler—Programmable multichannel data handler for interactive and automatic processing of seismological analyses. Comput. geosciences19 (2), 135–140. 10.1016/0098-3004(93)90110-Q
79
SteinbergA.GaeblerP. (2023). Automatisch erstellter Erdbeben Katalog des Seiger Projektes. Zenodo. 10.5281/zenodo.7973593
80
StutzmannE.SchimmelM.PatauG.MaggiA. (2009). Global climate imprint on seismic noise. Geochem. Geophys. Geosystems10 (11). 10.1029/2009GC002619
81
SuckaleJ. (2010). Moderate-to-large seismicity induced by hydrocarbon production. Lead. Edge29 (3), 310–319. 10.1190/1.3353728
82
SzuberlaC. A.OlsonJ. V. (2004). Uncertainties associated with parameter estimation in atmospheric infrasound arrays. J. Acoust. Soc. Am.115 (1), 253–258. 10.1121/1.1635407
83
TakanamiT.KitagawaG. (1988). A new efficient procedure for the estimation of onset times of seismic waves. J. Phys. Earth36 (6), 267–290. 10.4294/jpe1952.36.267
84
VanDecarJ. C.CrossonR. S. (1990). Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares. Bull. Seismol. Soc. Am.80 (1), 150–169. 10.1785/BSSA0800010150
85
VasterlingM.WeglerU.BeckerJ.BrüstleA.BischoffM. (2017). Real-time envelope cross-correlation detector: Application to induced seismicity in the Insheim and Landau deep geothermal reservoirs. J. Seismol.21, 193–208. 10.1007/s10950-016-9597-1
86
VellemanP.WelschR. (1981). Efficient computing of regression diagnostics. Effic. Comput. Regres. Diagnostics. Am. Statistician35 (4), 234–242. 10.2307/2683296
87
WaldA. (1943). Tests of statistical hypotheses concerning several parameters when the number of observations is large. Trans. Am. Math. Soc.54 (3), 426–482. 10.1090/s0002-9947-1943-0012401-3
88
WangR.SchmandtB.ZhangM.GlasgowM.KiserE.RysanekS.et al (2020). Injection‐induced earthquakes on complex fault zones of the Raton Basin illuminated by machine‐learning phase picker and dense nodal array. Geophys. Res. Lett.47 (14). 10.1029/2020GL088168
89
WeingartenM.GeS.GodtJ. W.BekinsB. A.RubinsteinJ. L. (2015). High-rate injection is associated with the increase in US mid-continent seismicity. Science348, 1336–1340. 10.1126/science.aab1345
90
WilsonD.LeonJ.AsterR.NiJ.SchlueJ.GrandS.et al (2002). Broadband seismic background noise at temporary seismic stations observed on a regional scale in the southwestern United States. Bull. Seismol. Soc. Am.92 (8), 3335–3342. 10.1785/0120010234
91
WithersM.AsterR.YoungC.BeirigerJ.HarrisM.MooreS.et al (1998). A comparison of select trigger algorithms for automated global seismic phase and event detection. Bull. Seismol. Soc. Am.88 (1), 95–106. 10.1785/BSSA0880010095
92
YaoY. C. (1988). Estimating the number of change-points via Schwarz'criterion. Statistics Probab. Lett.6 (3), 181–189. 10.1016/0167-7152(88)90118-6
93
ZangA.OyeV.JoussetP.DeichmannN.GrittoR.McGarrA.et al (2014). Analysis of induced seismicity in geothermal reservoirs – an overview. Geothermics52, 6–21. 10.1016/j.geothermics.2014.06.005
94
ZhuW.BerozaG. C. (2019). PhaseNet: A deep-neural-network-based seismic arrival-time picking method. Geophys. J. Int.216 (1), 261–273. 10.1093/gji/ggy423
Summary
Keywords
seismic arrays, induced seismicity, geothermal energy, epicentral localization, robust processing
Citation
Hering P, Lindenfeld M and Rümpker G (2023) Automatized localization of induced geothermal seismicity using robust time-domain array processing. Front. Earth Sci. 11:1217587. doi: 10.3389/feart.2023.1217587
Received
05 May 2023
Accepted
22 June 2023
Published
06 July 2023
Volume
11 - 2023
Edited by
Luca De Siena, Johannes Gutenberg University Mainz, Germany
Reviewed by
Vladimir Tcheverda, Institute of Petroleum Geology and Geophysics (RAS), Russia
Katrin Hannemann, University of Münster, Germany
Ferdinando Napolitano, University of Salerno, Italy
Updates
Copyright
© 2023 Hering, Lindenfeld and Rümpker.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Philip Hering, phhering@geophysik.uni-frankfurt.de
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.