the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Array processing in cryoseismology
Thomas Samuel Hudson
Alex M. Brisbourne
Sofia-Katerina Kufner
J.-Michael Kendall
Andy M. Smith
Abstract. Seismicity at glaciers, ice sheets and ice shelves provides observational constraint of a number of glaciological processes. Detecting and locating this seismicity, specifically icequakes, is a necessary first step in studying processes such as basal slip, crevassing, and imaging ice fabric, for example. Most glacier deployments to date use conventional seismic networks, comprised of seismometers distributed over the entire area of interest. However, smaller aperture seismic arrays can also be used, which are typically sensitive to seismicity distal from the array footprint and require a smaller number of instruments. Here, we investigate the potential of arrays and array-processing methods to detect and locate seismicity in the cryosphere, benchmarking performance against conventional seismic network-based methods. We also provide an array-processing recipe for cryosphere applications. Results from an array and network deployed at Rutford Ice Stream, Antarctica, show that arrays and networks both have strengths and weaknesses. Arrays can detect icequakes from further distances whereas networks outperform arrays for more comprehensive studies of a process within the network extent, due to greater hypocentral constraint and a smaller magnitude of completeness. We also gain new insights into seismic behaviour at Rutford Ice Stream. The array detects basal icequakes in what was previously interpreted to be an aseismic region of the bed, as well as new icequake observations at the ice stream shear-margins, where it would be challenging to deploy instruments. Finally, we make some practical recommendations for future array deployments at glaciers.
Thomas Samuel Hudson et al.
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2023-657', Andreas Köhler, 01 May 2023
In this study the authors apply seismic array processing to data recorded on an ice stream in Antarctica. The goal is to investigate the benefit of array processing (plane wave beamforming) compared to network processing for detecting and locating icequakes (migration and stacking). I do appreciate this study very much since I believe that array processing has a huge potential for analyzing complex and signals-rich cryoseismological data sets. A study comparing array and distributed network deployments on glaciers has not been done previously to my knowledge and is therefore very timely. It can be beneficial for researchers planning future field campaigns on glaciers where the number of sensors is limited, due to logistical reasons, to find the most optimal configuration for their research objectives. This is in contrast to recent large-N deployments on glacier, where multiple processing approaches can be tested after the measurements.
(1) One concern I have is not so much about the methods and results, but how this work is introduced and set in a broader context. The title seen alone suggests either a review article or the very first application of array processing for cryosphere monitoring. This is misleading since several studies have already been using arrays in this context for a while, and these works are not referred to in this manuscript. Array processing has been used on Alpine glaciers, in Greenland, in Svalbard and Antarctica. Please check the list of references below. Not all might be relevant and there might be even more. Some of these studies used permanent arrays (Antarctica; Svalbard). Others used temporary arrays deployments on or close to glaciers for the purpose of studying cryogenic seismic signals (icequakes, calving events, tremors, ...). Some do classic plane wave array processing; others apply matched field processing / general beamforming. The definition of what is array processing and what network processing in these studies may not be so clear always, however, I do believe these studies are related to your topic. Given these previous works and the sole focus on Antarctica in your study, I strongly suggest modifying the title to better reflect the content and contribution of this paper, and also update the text (mainly introduction) and references accordingly.
Examples of papers using array processing, beamforming, and FK analysis in cryoseismology:
https://doi.org/10.1785/0220200280
https://doi.org/10.1029/2021GL095996
https://doi.org/10.1017/aog.2018.25
https://doi.org/10.5194/tc-16-2527-2022
https://doi.org/10.14943/lowtemsci.75.15
https://doi.org/10.1029/2021GL097113
https://doi.org/10.5194/tc-14-1139-2020
https://doi.org/10.1002/2015GL064029
https://doi.org/10.1093/gji/ggac117
https://doi.org/10.5194/tc-13-3117-2019
https://doi.org/10.1002/2016GL070589
https://doi.org/10.3402/polar.v34.26178
(2) In Chapter 3.1.1 the array processing pipeline is introduced. Automatic array processing including phase detection by beamforming, phase classification by F-K analysis, and association using back-azimuth constraints, similar to your method, is routinely done at several national seismic data centers and at the IDC. This should be mentioned. See for example Chapter 9.9 in this one: https://doi.org/10.2312/GFZ.NMSOP_r1_ch9
(3) You use beam power time series P_Z for detecting P waves and P_H for S waves. As far as I understand the slowness is not used to actually confirm that a P or S waves are detected, which is usually done in automatic array processing. Have you considered this? Theoretically, S waves can also be triggered by P_Z, and P waves by P_H, depending on the incident angle of the phase arrivals. Could it for example happen that a horizontally arriving P wave leads to a peak in P_H which you would miss because you are looking only for peaks in P_Z followed by P_H? I acknowledge that your procedure is more sensitive to detect deep icequakes, which is maybe what you intended. Also, there is of course the ray bending due to the firn layer. Some comments on that would be appreciated.
(4) I am not sure about the motivation for Chapter 3.1.4. It is not clear where you do time-domain beamforming at this point. You described that you choose to do array processing in the frequency domain, i.e., F-K analysis in continuous data to measure the slowness vector, but nothing is mentioned about time-domain stacking in the detection and location procedure. Please clarify.
(5) I find it difficult to understand why your 3D location procedure fails so clearly, i.e., all events end up on a vertical line below the array. From the description of the methodology, I had the impression that the synthetic take-of angle and distance PDF should enable reasonable results. Could you add some sort of synthetic resolution test for the 3D method to investigate this? Your approach is by the way similar to how teleseismic events can be located using a single array. Using the 1D velocity models of the Earth, we can predict the horizontal slowness (ray parameter) of a P or S wave at a given distance. With the back-azimuth, we then have the location. Even if your firn model in not perfect, I expect more spread-out locations.
(6) Figure 6: Seeing these waveform examples makes me wonder if there could be many mis-associations in your detection list. Visually it is difficult to identify and relate the indicated P and S arrivals. I acknowledge that the authors have experience in identifying icequakes in that area. However, some proof for that these are not coincidental detections and associations would be appreciated. So how confident are you in the array-detected events?
Minor comments:
Chapter 1: I would already here emphasize that the set of 16 sensors is used for classic network processing later and add the dimension of it in the text, and that the set of 10 sensors used for array processing are located in the centre of the network. This would prepare the reader better for the comparison later. Side note: The 16 sensors together could also be considered as an array for processing longer wavelength signals, for example regional and teleseismic arrivals.
Line 4: Add calving as another source of cryoseismicity.
Line 100: You do you not explicitly write how you measure the take off angle. It is clear that it must come from the observed horizontal slowness and an assumed P/S wave velocity in ice, but good to mention this for the reader.
Line 105: How do you obtain the path-averaged velocities?
Caption figure 3: I think this piece of text has to be deleted: “n of an event alysis (FTAN)”
Line 172: The comparison of network and array processing on the array is very interesting. But you could already mention here that the target signals are of course very different. For plane wave array processing you need events at some distance from the array, whereas the migration stacking is expected to work best for events inside the network or at close distance.
Line 181: Here it would be very helpful to distinguish between the two network processing approaches (all sensor and array only). I would assume that network processing outperforming the array processing for small magnitudes is mainly due to icequakes detected inside the dense array. This is expected since array processing needs plane waves as mentioned above. This discussion is missing and should be added.
Chapter 4.1.1.: You give a good explanation for the peaks in the magnitude distribution. If my guess that the missing small events in the array processing results come mainly from events inside the array is correct, then these peaks could also originate from very close events which arrive with non-planar wavefronts. For evaluating this it would be very helpful to have the sensor location plotted in Fig 4c, and also show a zoom of array-detected events inside the network.
Line 213: double “with”
Chapter 4.2. As you write the comparison of the methods is difficult because not the same qc criteria are applied. One candidate for qc-ing array results would be the coherency or semblance for each detected phase. You could increase the threshold to see of the better constrained events resemble the network locations.
Fig 5: It looks like you show the cumulative and non-cumulative distributions, but the axes labels just says “cumulative”.
Line 265: “… icequakes from” ?
Citation: https://doi.org/10.5194/egusphere-2023-657-RC1 -
RC2: 'Comment on egusphere-2023-657', Anonymous Referee #2, 11 May 2023
In this article, Hudson et al. detect and locate glacial seismicity at Rutford Icestream (Antarctica) using both classical seismic networks and array-processing techniques. After explaining the methodology of both approaches including beamforming, waveform migration, triggering thresholds, etc, the authors compare the resulting event catalogs obtained from their deployment on Rutford Icestream. The main conclusion of this comparison is that the network-based approach detects less events compared to the array analysis. However, this advantage of arrays in terms of event number comes at the price that the individual event locations can be less good constrained, compared to using the network approach.
I do like the idea of applying beamforming to p- and s-wave arrivals separately and to subsequently combine the outcome to locate the events, because it does not require to detect and label the phases initially (the association is done based on the continuous beamforming results). This is a creative approach and an appropriate way to detect coherent arrivals close to or hidden in the background noise. Furthermore, the topic is well suited for the journal and provides an interesting discussion on deployment strategies in the harsh cryospheric environment. However, while the general idea is clear, I think there are some issues regarding the methodology including technical flaws and unclear formulations. For this reason, I am somewhat skeptical about the results and the conclusions drawn. Below, I provide major comments to address these concerns in detail, which are followed by a list of minor comments.
Major comments:
Title: I think the title is a bit misleading, as it suggests a general paper on array processing in cryoseismology or a review paper. As the article focuses on detecting and locating high-frequency basal icequakes, I would suggest something like “Array processing to study icequakes on Rutford Ice Stream, Antarctica” or similar.
Methodology section:
1) The authors state that the length of the time window used for f-k analysis is 0.025s, which is equivalent to 25 samples at 1000 Hz sampling rate. From this short window it will not be possible to get proper frequency spectra needed for the approach, in particular for the low frequencies used, i.e. down to 10 Hz. It is mentioned that this time window length is based on the “slowest time it would take for a seismic wave to travel across the array”, which in my opinion must be 90 m / 1 km/s = 0.09 s, which is still hardly enough to properly resolve 10 Hz. I would appreciate if this is carefully checked and corrected.
2) It is unclear to me, how the f-k analysis is conducted. Lines 74-78, which introduce the basic concept, are confusing in my opinion, e.g. why do you calculate power spectra for each receiver? What does point two mean - only that you calculate the beam power for the 2D slowness space? How do you calculate the results at discrete frequencies? As this is the core of this work, I believe it is necessary to clarify these points and I suggest to also add the formulas on which the f-k analysis is based.
3) The proposed method distinguishes p- and s-waves through the beam power on the vertical and horizontal components, respectively. I would like to note here, that Sv waves may map into the vertical component and p-waves that come from outside the array may map into the horizontal components (though it is noted that due to the firn layer there is a steep incident angle). In addition, Rayleigh waves will show up on both vertical and horizontal components. I think it would thus be helpful to take the p-to-s velocity ratio of an icequake detection as an additional criterion to trigger the location process. There should be a clear difference in p- and s-velocity reflecting the subsurface p-to-s ratio. In any case, I think the manuscript should also contain a time series of the slowness and back azimuth associated with the maximum beam power shown in Figure 3. This would also help the reader to better understand the basic concept.
4) Connected to the previous comment: the f-k method yields slowness/velocity values and large parts of the location analysis requires a velocity model. However, the manuscript does not report any values of seismic velocity, neither in the context of the beamforming results nor in the context of the velocity model. I strongly encourage the authors to include more detail to remove ambiguity.
Results and conclusions
Based on various aspects, I feel that the conclusions of the paper are drawn on somewhat little evidence. The authors admit (which is not a problem per se), that explicit 3D locations of icequakes using the array method are not possible (all events locate on a vertical line beneath the array, Fig. 4b). For this reason, the proposed method prescribes a fixed depth for all events assuming basal seismicity, only. However, repeated surface icequakes producing dominant Rayleigh waves (which can be the dominant event type on glaciers) would be picked up in the array analysis as well (as Rayleigh waves will peak on both vertical and horizontal beam power time series). The authors state, that false triggers are contained in the catalogue, e.g. such surface crevasse fracture events. For these events, the location method based on p-s arrival times and basal seismicity is not justified (or actually violated), which may lead to erroneous event locations and interpretation of the seismicity. Finally, in the presented examples of newly detected events using the array (Figs. 6c and 6d) I cannot discern obvious arrivals, which strengthens my view that the catalogue may contain many false triggers. In total, due to the various issues outlined above, I feel that better evidence needs to be presented to support the conclusions. I think that this may be possible by choosing a more strict trigger criteria/thresholds.
Other comments:
Line 7 (abstract): I think here it is necessary to specify that the method is tailored to high-frequency (>10 Hz) events with distinct p- and s-wave arrivals. For events with dominant Rayleigh waves (which are very numerous in on-ice recordings), the method will not work, as the p-s criterion breaks down and as there will be a peak in beam power on both the vertical and horizontal components.
Line 9 (abstract): Also here, “cryosphere applications” is too general in my opinion (see above).
Line 22-23: I guess with “such icequakes” you mean basal seismicity? Because especially for surface crevassing events, there are several studies that use arrays, e.g. Mikesell et al. (2012), Köhler et al. (2016), Lindner et al. (2019).
Lines 24 and 26: I am having trouble with the term “sensitive”, as this implies that a network would not sense any events outside the network and an array would not sense events within the array, which of course is not true. So maybe “designed” would be more appropriate here?
Line 28: Is it really justified to generally say that arrays enable event detection at greater distances? In the end, a network is a large array of sensors as well.
Line 35: Problem with reference rendering
Line 37: I think there is a problem with the references: only the paper by Klaasen et al. looked into icequakes, but not the other two papers. Please check.
Line 38 and following: Actually, there are a number of other examples of array analysis and beamforming to locate icequeakes, e.g. Mikesell et al. (2012), Köhler et al. (2016), Lindner et al. (2019).
Lines 47-51: Can you provide a few more details on sensor installations and deployment length? Also, it would be nice to get some more context on the study site.
Lines 93-96: I do no understand the third criterion, can you please rephrase it?
Line 106: Maybe add here that the fixed depth method assumes that the events originate at the glacier bed.
Line 107 and following lines: I cannot follow, how the 3D location method works. Also I am afraid that figure 2e does not help to get the concept. I suggest to rewrite this paragraph.
Line 112: Which velocity model are you using here and how does it look like? I think this is important to mention as you later say you “assume that there is negligible uncertainty in the velocity model.”
Line 134: Why “approximate” seismic velocity?
Line 135 and following lines: The provided frequency limits assume slightly different velocity values for the min and max receiver spacings. Please clarify and mention which velocities are used.
Line 140: This sentence implies that you do not use a beamforming technique in the frequency domain? In an earlier part of the manuscript, you refer to this paragraph for the stacking of the beamforming maps in the 2D slowness domain, but here waveforms are mentioned. What are you stacking here, please specify.
Line 148 and following lines: I guess also here, the velocity model is crucial.
Lines 153-155: Why do you need to relocate the events from the QuakeMigrate method, please specify.
Line 178: Why is there a sharp detection limit at ~7.5 km for the network?
Lines 183-185: Does that mean that the results of the array-based method contain false-triggers?
Line 213: Twice “within” at the end of the line.
FIGURES
Figure 1: Most of the subfigures are too small, it is therefore not possible to get the geographical context of the field site. Please make bigger.
Figure 2: Again many fonts are too small. In panel c: numbers are missing on the slowness axis. Maybe also use “beam power” instead of “power” only (throughout the manuscript) to avoid misinterpretation?
Figure 3:
As suggested earlier, it would be great to see also the slowness and back azimuth time series together with the beam power.
Where do the focal mechanisms come from?
Legend: What is “FTAN”? “events.n” > events?
Figure 4: What are the different gray scales of event locations in panel a? Panels c and d: I suggest to also show the locations of the geophones here.
Figure 5: This figure is again too small. I am also having trouble to understand panels a and b: maybe I am misunderstanding something here, but why are there two values of the cumulative event number for a given magnitude? E.g. in panel b, for Mw -0.5, the cumulative event number is both ~3*10**3 and ~6*10**4?
Figure 6: The subpanels showing phase detection are too small.
Legend: “reion” > region
References:
Köhler, A., Nuth, C., Kohler, J., Berthier, E., Weidle, C., & Schweitzer, J. (2016). A 15 year record of frontal glacier ablation rates estimated from seismic data. Geophysical Research Letters, 43(23), 12-155.
Lindner, F., Laske, G., Walter, F., & Doran, A. K. (2019). Crevasse-induced Rayleigh-wave azimuthal anisotropy on Glacier de la Plaine Morte, Switzerland. Annals of Glaciology, 60(79), 96-111.
Mikesell, T. D., van Wijk, K., Haney, M. M., Bradford, J. H., Marshall, H. P., & Harper, J. T. (2012). Monitoring glacier surface seismicity in time and space using Rayleigh waves. Journal of Geophysical Research: Earth Surface, 117(F2).
Citation: https://doi.org/10.5194/egusphere-2023-657-RC2
Thomas Samuel Hudson et al.
Thomas Samuel Hudson et al.
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 162 | 73 | 8 | 243 | 1 | 2 |
- HTML: 162
- PDF: 73
- XML: 8
- Total: 243
- BibTeX: 1
- EndNote: 2
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1