We thank the reviewer for the appreciation of our work and the very helpful, careful and detailed review.

R1: GENERAL COMMENTS

This is an interesting paper dealing with the application of array techniques to the problem of earthquake location. It constitutes a very clear demonstration of the fact that when a distributed network around the seismogenic area is not available, seismic arrays are a feasible alternative. The authors deal with earthquakes that are distributed around Fogo and Brava islands. Any seismic network based on land stations will be ill-suited, due to the limited extent of the Cape Verde islands. The authors use time-domain beam-forming at three seismic arrays to estimate the back-azimuths to the source, and design a scheme based on probabilities to pinpoint the epicenter locations.

However, there are some aspects of the work that need a critical revision. I address some of them below.

EPICENTRAL LOCATION

The locations obtained by the multi-array approach are limited to the determination of the epicenters using the estimates of back-azimuth. But sometimes the authors seem to oversee this limitation. At several points in the paper (see detailed comments below) you seem to imply that your method is better that others because you do not need a velocity model, with all the uncertainties that it brings to the calculations. However this is not a real advantage, since you provide only epicentral locations. The determination of depth requires some knowledge or assumptions about the velocity structure.

There is some discussion about source depths, although they are heterogeneous and somewhat confusing. In general, depth estimates are not addressed systematically. Sometimes a velocity model is assumed, either a two-layer model (line 229) or the model of Vales et al. 2014 (line 203) which is more complex. At one instance the depth is fixed (line 230). Depths are obliterated at the beginning but they are brought back during the discussion. I suggest that, if possible, you should exploit fully the information provided by the arrays, not only the back-azimuths but also the apparent slowness results, and include depths as part of the results. Of course the depths would be biased by the choice of velocity model, but at least you would have an estimate.

Answer: It is true, that the method we chose only reveals epicentral locations. This may seem like a disadvantage, but it provides the possibility to determine the epicentral location without any a priori knowledge about the velocity structure. In the case of Fogo and Brava, we decided not to rely on the velocity model, as the available models either provide an average structure for the whole Cape Verde archipelago or of other islands with older structures than Fogo and Brava.

Nevertheless, we agree with the reviewer, that we need to point out more strongly, that we only determine epicentral locations. We will modify this in the revised manuscript.

R1: EARTHQUAKES VS HYBRIDS

The classification of seismic signals in “earthquakes” and “hybrid events” is not clear. Figure 7 compares the seismogram and spectrogram of a volcano-tectonic earthquake and a hybrid event, both recorded at the AF array on Fogo. Why don’t you think that the hybrid event is in reality a small volcano-tectonic earthquake located in Fogo? I don’t see a low-frequency coda, or any hint of a band-limited frequency content. Looking at the spectrum, it is true that there is relatively more energy in the 0-10 Hz band. For the hybrid, the peaks in the 0-10 and 10-20 Hz bands have similar heights. While for the VT earthquake there is definitely more energy in the 10-20 Hz band. Still, this could be an effect of lower signal-to-noise ratio. The noise seems to have bursts of energy in the 0-10 Hz band (before, during and after the event), so the difference could be just an effect of the size of the earthquake. Figure S2 shows a single volcano-tectonic earthquake recorded at the AF array on Fogo and the BR array on Brava. The right plot looks similar to Figure 7b, except for the difference in signal-to-noise ratio.

In any case, I suggest that you justify your classification, showing hybrid events that look clearly different from local earthquakes. This has no impact in the method that you are presenting, which is valid for any earthquake: volcano-tectonic, hybrid, or whatever. But it is important for the interpretation of the activity in the discussion.

A.: In contrast to a volcano-tectonic earthquake, a hybrid event shows a smooth transition from higher to lower frequencies. In the case of a hybrid event there is no S-phase identifiable and in the coda (for about 15 to 20 seconds) there is still more energy available in the 0-10 Hz band than before the event (and in comparison with a VT earthquake). In contrast, even for the small VT earthquakes on Fogo two phases are detectable. To make this point clearer, we will replace the example earthquake displayed in Figure 7a with an earthquake from Fogo. By this, the difference between the event types should become clearer. Additionally, we will add a figure to the supplementary material which shows the spectrogram of all components of a hybrid event. In this figure the low frequency content of the coda becomes evident.

With our discrimination between volcano-tectonic earthquakes and hybrid events we also follow earlier studies on the seismicity of Fogo (e.g. Faria and Fonseca, 2014).

We do not exclude that an earthquake might be involved in the triggering of the hybrid events. McNutt (2000) refers to a mechanism, where an earthquake sets an adjacent fluid-filled cavity into oscillation. This could explain the high-frequency content in the beginning of the signal and the lower frequencies in the coda. Nevertheless, there are also other hypotheses of possible mechanisms for hybrid events, which are also capable to explain their occurrence on Fogo (see discussion in the text).

R1: ARRAY ABERRATIONS

I don’t agree on the use of the term “aberration” referred to the array, e.g. line 64: “seismic arrays can exhibit systematic aberrations of backazimuth and slowness”, line 186: “Effects along the ray path from the source to the array, such as heterogeneities, can result in a systematic aberration of the array”. The heterogeneity of the medium can have an effect on wave propagation, resulting in rays impinging at the array site with back-azimuths that do not point to the source. See for example the interesting results of Garcia-Yeguas et al. GJI 2011. However, this cannot be pinned to the array! The seismic arrays are working ok, providing the directions with which the wavefronts are propagating through the array, these estimates are not biased in any way. The wavefield distortions (or aberrations) are produced in the medium. Therefore, I’d rather talk about ray bending or back-azimuth deviations, instead of array aberrations.

A.: Thank you for pointing this out. In the revised manuscript we will change the terming accordingly.

R1: METHOD

The description of the method, and particularly the uncertainty estimates, is a bit confusing. You take as uncertainty the standard deviation of the maximum energy obtained by varying the start and duration of the stacking window. How do you obtain the standard deviation of back-azimuth and apparent slowness? Are these the values reported in tables S1,S2? Standard deviations of back-azimuth in those tables are large (around 100 degrees), but in the figures the beams seem to be much narrower. The width of the back-azimuth wedges is related to the uncertainty of the solutions, but I don’t fully understand this relationship.

A.: Yes, the standard deviation values are reported in the tables S1 and S2. Regarding for example Figure 6, the beam of array AF (southern beam at Fogo) exhibits a standard deviation of 73°. This corresponds to an error of 20%. The beam width thus corresponds to 80% of the maximum energy of the “original” beam determined during the array analysis (i.e. the beam resulting from the manually chosen stacking window). In this way, the backazimuth wedges are related to the width of the “original” beam. This beam is not necessarily symmetric and the choice of the errors taken from the standard deviation of the backazimuth allows us to include possible asymmetries of the beam. Additionally, in this way also sidelobes are included in the multi-array analysis. This can be seen in Fig. 6, where a small beam at array CG (western array on Fogo) points to the south-southwest.

However, we will clarify this point during the revision and include the details mentioned above.

R1: Additionally, the width of the array response is not considered at all. Even if the energy stack in figure 3 provides a very narrow peak, you should consider the uncertainty produced by the finite array response. In figure S1 I see that the central peak has a diameter of about 0.1 s/km, this should be always added to any uncertainty estimate.

A.: One advantage of the choice of a time-domain array analysis is the implicit inclusion of a broad frequency band in the analysis. The conventional array response is only estimated for a single frequency. Correspondingly, the stacking of (integration over) the array responses for a wide frequency band would also lead to a much narrower peak.

We will add this information to the supplement.

R1: For earthquakes, apparent slowness is smaller at BR. This should imply a larger uncertainty in back-azimuth, which is not clear in the plots.

A.: Yes, this is true, BR shows typically large uncertainties in the backazimuth. Therefore, in many cases it shows the highest values for the standard deviation of the BAZ. However, the higher uncertainties also cause many energy contour plots (used for the BAZ estimations) of the array BR to be rejected, as the results are unstable. For example, if a small variation in the choice of the stacking window or different frequency bands lead to a strong variation of the beam, this result is considered unstable and thus rejected.

R1: Why do you use the time-domain beam-forming method? There are other options that may provide improved estimates of the back-azimuth and apparent slowness. Have you tested any other array method?

A.: We have also tested an f-k analysis (software of NORSAR). However, we think that the implicit inclusion of a broad frequency band, which leads to a strongly narrowed energy peak, and the narrow stacking window around one particular phase at the reference station provide strong advantages compared to the conventional f-k analysis. As described in the text, time-domain analysis has additional benefits with regards to the phase of interest.

The implicit inclusion of a wide frequency band in the time-domain analysis not only reduces the amount of possible sidelobes, but also allows for the analysis of events with different frequency content.

R1: How are the results from the three arrays combined to provide the locations? (as in figure 6). It looks like you simply add the “energy” grids obtained for the arrays. However, if you interpret these maps in terms of probabilities of a grid node being the epicenter location, it make more sense multiplying them. After all this is an “and” operation, you want to determine if the position of a given node complies with the back-azimuths at array AF AND array CG AND array BR.

A.: We do not interpret the maps in terms of probabilities of a grid node. But instead we use the sum of the normalized energy distribution derived from the single arrays. We found this method to give the best results in comparison to the conventional network localization.

R1: NETWORK LOCATIONS

I’d like to have some information about the quality of the network locations. They are used as ground truth in the comparison with the array locations. But the station coverage is poor and I assume that the network locations may have large uncertainties as well, even if you restrict them to the best cases. Perhaps you can show examples with error ellipses (in the supplementary information?).

A.: This is a good point; we will add a map to the supplementary material that shows the locations of the (classical) network-based analysis.

R1: OTHER COMMENTS

- line 54: the method of localization “has the advantage of being independent of velocity models”. This is a bit misleading because you are just providing an epicentral location. If you try to locate the earthquake (including depth) you do need a velocity model.

A.: Thank you for pointing this out. In the revised manuscript we will clarify that our method provides epicentral locations.

R1: - line 55: “The velocity structure is often very complex in a volcanic regime”, I think that the term “regime” is not adequate here. Perhaps “environment”, “setting”, or “medium”. (The same in line 242).

A.: In the revised manuscript we will replace the term “regime” by “environment”.

R1: - lines 57-62: this description of the method is too detailed for the introduction, it could be later in the method section.

A.: We will shorten this appropriately.

R1: - section 2, the arrays are made of different types of seismic stations, short period and broadband. I assume you have taken good care of this and removed the instrument response, in order to make the waveforms comparable?

A.: We remove the instrument response during the magnitude determination. However, during the array analysis we apply appropriate filter functions.

R1: - line 73: “the arrays were designed for events with frequencies between 5 and 10 Hz”. What do you mean? More than frequency of events, what’s important in array design is the wavenumber, that includes the effect of frequency and apparent slowness.

A.: Thank you for pointing this out. We decided to use the array transfer function in terms of frequency and slowness components s_x and s_y (instead of wavenumber). Both formulations are equivalent.

The arrays were designed based on the array transfer function for frequencies between 5 and 10 Hz. We will clarify this in the revised manuscript and change the wording.

R1: - line 83: “absolute slowness”, what do you mean with “absolute”?

A.: This term refers to the estimation of the slowness via s=√s_x²+s_y². However, to make this clear, we will replace the term “absolute slowness” with “magnitude of horizontal slowness”.

R1: - figures 3,4,5: one figure for one array may be enough to show how the method works.

A.: We decided to show all three figures, as their results are combined in Figure 6. However, we will move Figures 4 and 5 to the supplementary material.

R1: - line 103: grid size of 124x124. If the grid search extends from -0.3 to 0.3 s/km, the grid spacing is about 0.005 s/km. This is unnecessarily small, given the shape of the array response.

A.: The array response function is only estimated for one particular frequency. However, the time-domain array analysis implicitly includes a broad frequency band, which narrows the central peak and reduces sidelobes.

Additionally, we performed tests with a grid size of 62x62. This would reduce the computation time. However, it turned out, that the results of the analysis with the finer grid were better.

R1: - line 104: “the resulting energy stack is shown in Fig 3b”. What do you mean by “energy stack”? You delay and sum the signals in the time domain, thus obtaining a seismogram representing the beam for a given apparent slowness and back-azimuth (e.g. green line in figure 3c). But this is still a time series, how do you get a unique value to be plotted in figure 3b? Is this the maximum amplitude of the beam? Is it the rms?

A.: This is the contour plot of the energy, representing the absolute amplitude of the sum trace for each grid node. We will replace the term "energy stack" by "contour plot of the energy" in the revised manuscript, which should clarify this point.

R1: - line 106: the formula BAZ=90-atan(sx/sy) does not follow the usual convention that BAZ=0 for north direction. When sx=0 we obtain BAZ=90. This formula should be corrected, or else explain your convention for azimuths.

A.: Thank you for pointing this out. Accidentally, the given formula is related with the MATLAB coding that we used. We will replace this with the correct formula BAZ=(180°/π)atan(s_x/s_y).

R1: - line 109: “the slowness is related to the angle of incidence by sin(i)=s*vc, with the mean crustal velocity vc”. I have two criticisms. First, this s in the formula is APPARENT slowness (or ray parameter). Slowness (without “apparent”) is just the inverse of velocity, and therefore it is not related to the incidence angle. Second, this formula is valid at points along the seismic ray, with local values of i, s_ap, and velocity. I don’t think that the use of a “mean crustal velocity” is adequate in this context. What are i and s_ap then? If you are thinking of the arrival at the array, i is the incidence angle at the surface, and s_ap the apparent slowness measured by the array. But then v should be the velocity of the shallow layer, not an average of the crust.

A.: Regarding the first point: We will adjust this in the revised manuscript by using the term “horizontal slowness (ray parameter)”. We prefer to use the terms “horizontal slowness” or “ray parameter” as the term apparent is usually used in relation to the velocity.

Regarding the second point: Thank you for pointing this out. In the revised manuscript we will rephrase this by using the term “velocity of the upper layer beneath the array”.

R1: - figure 6: there is a light-blue beam going SSW from the CG array, what is it? (the same happens in figures 12 and 13). The error bars seem underestimated. In the color scale of this plot, 90% of the maximum seems to correspond to the transition from yellow to orange. But in the map, the error bars do not include that transition, I’d say that the error area should include the whole yellow patch, e.g. +-0.03 degrees approx from the epicenter.

A.: The light-blue (small) beam pointing to the south-southwest from array CG in Fig. 6 is an example of a sidelobe, which shows an energy corresponding to the standard deviation of the backazimuth. We will add this information in the captions.

Unfortunately, a mistake occurred during the estimation of the error. In the revised manuscript we will correct for this mistake. The error bars are indeed far larger than displayed in the original figure. Thanks for bringing this to our attention.

R1: - figure 6, caption: “the location of maximum energy estimated” is unclear. You are locating the epicenter of an earthquake. The colors indicate the directions in which the beamforming provides maximum energy, but in what sense the overlapping area has more energy? Also the sentence “The beams correspond to the beams” should be rephrased.

A.: We will rephrase the figure caption in the revised version of the manuscript.

R1: - line 156: “discrepancy of the amount of detected and with the multi-array analysis located earthquakes”, weird sentence, rephrase?

A.: In the revised manuscript this will be rephrased.

R1: - line 164: “on Brava the dominant frequencies of the same event are lower”, why do you mention this, and even show figure S2? Since the earthquakes are closer to Brava than Fogo, is this an indication of a site effect?

A.: In the discussion we go into more detail on this point (lines 240-245). This is the reason why it is also mentioned in the description of the observations.

R1: - line 166: “ray turning point”, that is true for rays propagating downward from the source. For an upward propagating ray there is no turning point.

A.: Thank you for pointing this out. Indeed, this estimate can only be performed for the earthquake signals recorded on Fogo and originating around Brava. This larger distance ensures that the recorded wavefront propagated downwards from the source. However, we will clarify this in the revised manuscript.

R1: - line 180: you superimpose the beams for all hybrids, why would you do that? Are the beams weighted in some way? Why didn’t you do this for earthquakes as well? The white boundary in figure 10b leaves out the largest hybrid events, e.g. the northernmost, magnitude 1.5 event. Do you assume that the sources must be the same? What is the interpretation of this figure 10b?

A.: From our analysis it becomes clear that the hybrid events only occur in a certain area on Fogo. The superposition of the beams of the hybrid events gives an idea of the area with the highest probability of occurrence of these events. Under consideration of the uncertainties of the events (errors and possible BAZ deviations), they could possibly share the same source region. The 80% contour line (white line) leaves out the strongest event. However, its error bars are rather large. Therefore, it is not unlikely, that it could be located closer or within the area of 80% probability of the occurrence of hybrid events.

This approach works well for hybrid events, which are confined to a smaller area than the earthquakes. Additionally, it is known that the earthquakes around Brava frequently change their location. Therefore, such an approach is not useful for the earthquakes we located.

R1: - figure 8: the cumulative number of hybrids seems unrelated to the daily number of hybrids. The curve grows when there are no events (e.g. during April or June 2017), and the peak with 5 events in a day is not coincident with any increase in slope. I think this may be due to building this curve just joining points, when it should have a stair-like aspect (you should have right angles instead of the diagonal). The difference is not significant when there are many earthquakes, as if figure 8a. But it pops out when the number of earthquakes is small, as in figure 8b.

A.: Thank you, this is indeed the case and we will modify the figure accordingly.

R1: - figure 9: very small events (M~0-0.5) near array AF and in areas near Brava, how can you locate these? Do you get good results for these earthquakes even at distances of 30 km from the array? This could be very interesting and emphasize the power of seismic arrays, but of course you need very fine noise conditions to perform the analysis, specially since beam-forming is carried out in the time domain.

A.: Indeed, the noise conditions are the reason, why many of the events with magnitudes below 0.5 are located by using only two arrays. In these cases, the noise condition of the third array typically are too poor. Additionally, they have in common that they occur at nighttime, where man-made noise typically is very low. Regarding the traces, it would indeed be difficult to pick the P-phase at all stations reliably. However, their combination during the beamforming results in a clear P-phase onset of the sum trace. The application of the time-domain analysis is helpful in such a case, as the frequency band of the analysis can be chosen favorably in order to increase the SNR.

We will include this point in the discussion of the revised manuscript.

R1: - figure 10 has error bars in the locations of the hybrid events, but there were no error bars in figure 9.

A.: The reason for leaving out the error bars in figure 9, is that the number of events is larger than in figure 10, which strongly reduces the clarity of the map. However, this is a good point and we will add a map to the supplementary material showing the error bars of the earthquakes.

R1: - section 4.3 is partly related to the method, in fact you refer here to section 3.2. Wouldn’t the part on array parameters fit better in the method section?

A.: We agree that there is a relation to the method section. However, we first decided to include this part in the results section, as we present the outcome of different tests of the method based on recorded events. Nevertheless, we will move this section to the method section.

R1: - around line 195: it is not clear if you are talking about the start of the stacking window or the length of the stacking window.

A.: We refer to the start time of the window and the end time.

R1: - line 206: rms < 0.25 s?

A.: Yes, this is the rms of travel time residuals. We will add the unit.

R1: - line 207: “(theoretical) back-azimuth”. I don’t think you should call “theoretical” to the back-azimuth obtained from the network location. The network location constitutes just another estimate of the epicenter location. And given the sparse network distribution, there is no reason to think that the network location will be better than the array location. So I would pick another word, perhaps “reference” back-azimuth?

A.: Thank you for pointing this out. We will replace the term “theoretical” by “reference”.

R1: - line 207: “slowness values”, you mean apparent slowness? For this you do need a velocity model, I assume you use the same velocity model from Vales et al. (2014) ?

A.: Yes, as for the localization with the classical localization method we also use the velocity model of Vales et al. (2014) for the estimation of the magnitude of horizontal slowness (ray parameter). We will include this information (magnitude of horizontal slowness; velocity model) in the revised manuscript.

R1: - figure 11: this is difficult to understand, since you represent back-azimuth and apparent slowness independently. The red dots are located in the coordinates (baz1,s1) given by the array analysis. But the green dots are not located at the coordinates (baz2,s2) given by the network locations. Instead, they are located at (baz2,s1) in the left panels and (baz1,s2) in the right panels. Why splitting them like this? These green points do not represent real estimates, and I find this figure very confusing, I don’t know what you want to emphasize with this kind of plot. Have you considered x,y plots instead of polar plots?

- line 211: “yields back-azimuths pointing too far to the south by about 7 degrees”, this is difficult to see in figure 11 because all points overlap in a tight cluster.

A.: We have chosen this type of plot because it is actually easier to read than the plot showing the vectors from (baz1,s1) to (baz2, s2).

We do not use x,y plots as especially for the baz components it is easier to see the deviation in a polar plot. Nevertheless, we will modify the figure to make the deviations more obvious.

R1: - figure 12: these examples may serve as a test to understand if your selection of uncertainty levels is adequate. Ideally, the uncertainty region in panels a and c should contain the solution shown in panels b and d. For example, if you choose a level of 0.7 to define the uncertainty region, the solution of panel b would be contained within the uncertainty region of panel a. It is important to understand that even if you use the maximum as best solution (red dots) the source could be anywhere within the uncertainty region (e.g. the orange-yellow areas, sum of energy above 0.7).

A.: This is an interesting point. In fact, the yellow area in e.g. Fig. 12a would extend even further to the west if we would choose a map section with limits further to the west. Therefore, the question arises whether the result can be used at all. If the uncertainty becomes too large, the result is not reliable and is discarded.

These examples (Fig. 12 a,c) are not appropriate to check the level of uncertainty as we would not consider these cases any further for the event localization.

R1: - line 249: “cannot establish a link between them”. I agree, I think that even bringing up the subject is a bit far-reached.

A.: It is true, that it is hard to establish a link between them. However, volcano-tectonic earthquakes are rare on Fogo and therefore it is worth noting that they occur in the same area as the deep earthquakes. Due to the difficulty to establish a link between these events, we deliberately kept this point rather short.

R1: - line 260: “from February to March and from June to September”. Looking at figure S6 I see that there were more than 10 events in June and from September to November. There is a minimum in August. This does not correspond with your description, there is no increment in the hybrids in the period “June to September”.

A.: Yes, indeed, this should read “from September to November”. Thank you for pointing this out. We will correct this.

R1: - line 266: “With the multi-array analysis it is not possible to estimate the depth of the events”. This is not true in general. You can estimate the depth as well, if you exploit the information from apparent slowness and not only back-azimuth, as demonstrated for example by Almendros et al. 2001 and La Rocca et al. 2004. If you want to keep this sentence, you should specify that you are referring to your own implementation of multiarray analysis, focused on epicenter location.

A.: Thank you for pointing this out. We will clarify the statement, that the multi-array method we apply is not suited to estimate the event depth as we do not include a velocity model in our analysis.

R1: - line 281: “Low rupture velocities and strong path effects result in the long low-frequency coda”, regarding this discussion see also the paper by Bean et al. 2014 in Nature Geoscience, who discuss the same issue for long-period events.

A.: Thank you for pointing out this interesting reference. We will include it in the discussion.

R1: - line 283: “we conclude”. It is not clear to me the reason of this whole discussion about the origin of hybrids. It is not based on your results, and I don’t see what is the consequence regarding the array analysis. Perhaps you could explain a little more, is there any evidence in the locations, or the array results, pointing to a fluid-related origin? Is there any evidence pointing to the contrary?

A.: “Conclude” may be a bit strong, but we think that this discussion is important. In the revised manuscript we will rephrase this sentence.

R1: - line 288: “partly strong”, what do you mean?

A.: We will remove the term “partly strong”.

R1: - line 289: “slowness variations”, compared to what?

A.: Slowness variations compared to a simple homogeneous velocity structure as also mentioned in the text.

R1: - line 302: “Being independent of any velocity model and able to locate events”, as commented somewhere above, this is misleading. The multi-array method proposed does not use a velocity model BUT is able to identify just epicenters. You should make this clear at all instances, claiming otherwise is misleading. The location of seismic events (hipocentral location, complete with depth) requires a velocity model.

A.: As stated above, we will specify this at the relevant instances by adding the information that we determine the epicentral location of the events.

R1: - line 313: parallel beams imply over-estimated distance? What do you mean exactly? To determine if the distance is over-estimated, you must know the distance. And this is precisely what you are trying to estimate.

A.: Figure 12 shows what is meant by “over-estimated distance”. In view of these observations we included a simple epicentral distance estimation based on S-P-travel time difference estimations (Figure 13) to give additional information to the analyst. This is explained further above (lines 220-234). This will be clarified in the text.

R1: - line 314: “closer to the expected location”. What is the expected location?

A.: See above. However, we will add a link to Fig. 13 at this point.

R1: - line 326: “This application allows the event localization without assuming a velocity model”. Again, this is misleading. In reality you should talk about epicentral location. You cannot fix the depth without a velocity model.

A.: We will clarify this during revision.

R1: - line 341: hybrids show significantly larger apparent velocities than volcano-tectonic earthquakes? You said that in average, at the Fogo arrays the apparent velocity is 7.1 km/s for volcano-tectonic earthquakes from Brava (line 165), and 7.8-8.4 km/s for hybrids (line 178). I don’t think this is “significantly larger”, specially considering the uncertainties in these estimates.

A.: We will rephrase this.

R1: - line 348: the mention of volcanic tremor in the conclusions of the paper is out of place.

A.: This information is thought as an outlook and to provide motivation for further analyses of the data.

Powered by TCPDF (www.tcpdf.org)

We thank the reviewer for the careful review and comments and for appreciating our efforts.

R2: The manuscript presents an interesting study to locate earthquakes using a number of arrays. Potentially it is a good manuscript for the journal. I However, before this I have several comments and concerns that I think will have to be addressed. Especially I have the feeling there are two different ideas here. I am not sure the multi-array analysis is the important part? Or is it? The discussion seems to focus on the events rather than the method. Maybe that makes the manuscript somewhat difficult to read?

The first main point I have is that the manuscript is fairly difficult to read because of the organisation of the text. For example the introduction talks about the methods and results in far too much detail without the necessary background that then comes in the chapters afterwards (but a bit more introduction and references would be useful in the introduction).

A.: This manuscript presents a relatively new method (use of multi-array and time-domain array analysis) applied to volcanic seismic signals on Fogo and Brava. Therefore, it contains a more detailed description of the method and a discussion of the events analyzed.

We will rephrase the introduction slightly. However, we will focus on previous studies using multi-array analysis at volcanoes and on previous observations in our study area.

R2: In chapter 3.2 the first sentence has nothing to do with the rest of the paragraph, while the remaining paragraph seems to fit better in chapter 3.1. Part of 4.1 should be in methods (it discusses errors of the analysis and should be with methods..). The part on error estimation should follow directly from the methods for better flow before then talking about the locations and the discussion, which focusses on the events rather than the method. Please re-write. Also the discussion introduces new results, so these should be in the results section?

This whole organisation makes the reading hard work. Re-writing would clear up much of this difficulty. Indeed the writing in many places could be much more concise there is often unnecessary detail and repetitions.

A.: The first two sentences of section 3.2 introduce the multi-array analysis. This is essential and all following descriptions are based on that idea. This is what follows after the beamforming, which is performed for each array individually. The section 3.2 contains the information on how the results of the beamforming of each array are combined to a multi-array solution.

However, we agree that the structure can be improved and we will revise the manuscript and move section 4.3 to the method section 3.

R2: Using beamforming, information on the velocity is needed, I am not sure how the authors can say that there is no need for knowledge of a velocity model? Have they tried to run the analysis with a non-fitting velocity model? I assume this would greatly worsen the resolution? This should be shown… (it comes back in line 90, so please comment on the size of these uncertainties and what it means to the location errors.)

A.: In array seismology, during the beamforming, no information about the underlying velocity structure is needed. The method is based on the idea, that a plane wave travels across the array. By shifting the traces and summing them up, the direction and the magnitude of slowness of the wavefront are determined (lines 82-86). For details see Rost and Thomas (2002) and Schweitzer et al. (2012), as also cited in the text. As the beamforming does not rely on a velocity model, the combination of the beams allows us to determine the epicenter of the earthquake from the intersecting beams alone, without knowledge about the velocity structure.

R2: Has stacking been done using a plane wavefront or a circular one? Is the wavefront already plane in this distance? Has this been tested?

A.: This is a good point, we assume a plane wavefront. This approximation is in accordance with conditions given in Schweitzer et al. (2012). We estimated that for distances of more than 6 km the plane wave approximation holds. This is the case for the majority of events we analyzed. However, for the few events being located closer to one of the three arrays, the magnitude of slowness becomes small at this particular array and the backazimuth shows a broad uncertainty due to the broadening of the beam. This uncertainty is accounted for during the localization with the multi-array method.

R2: Has the stacking been done with one backazimuth? How different is the backazimuth for each station? This should be shown.

A.: The stacking is done for a complete range of backazimuths by a grid search over the range of horizontal slowness components. Therefore, it is not determined per station but for the whole array.

R2: How much does out-of-plane travelling influence the results? In the presence of strong velocity changes (as possible in volcanic environments), the waves will not travel on the great circle path. Has this influence been tested? This is different from mislocation vectors as discussed later.

A.: If the velocity changes are strong enough to impact the propagation of the wave, this should also lead to errors in a classical localization method. We constrain the possible errors by comparison with the network localization.

R2: The discussion seems to be disconnected from the rest of the manuscript? Where do all the results come from? They were not shown in results? This needs to be better organized as well. Please re-write

A.: The approximations described in the discussion section serve as additional information on the interpretation of the events analyzed. However, these estimations are only rough approximations (e.g. the estimation of the ray path of hybrid events to estimate the depth) and are often too simple to give a reliable result (for example the depth estimation of the hybrid events results in too large depths, which are in contrast to the observations pointing to a shallow source region). Therefore, such approximations should not be shown in the result section, as this may give the impression that these results are considered reliable. We therefore think that these approximations are best discussed in the discussion section.

R2: What are the systematic aberrations? Is this a mislocation vector? This needs a reference! - ah pages later I find that it is indeed the mislocation vector as given by Krüger and Weber. Please move this part of the manuscript when first discussion aberrations. It is out of place that far back.

A.: Thank you for pointing this out, we will add a reference at the corresponding position in the text.

R2: Has for the estimation of the mislocation vectors in part 4.3 the topography been taken into account? I assume that the stations are located at different heights? Jacobeit et al., 2013 show that topography has a large influence on the mislocation, not only heterogeneity…

A.: It is true, that the station elevation difference can have a strong influence on the result of the array analysis. Therefore, we carefully tested possible influences under the assumption of the different station elevations according to Schweitzer et al. (2012). It turned out that the station elevation differences are small enough to be neglected.

We will add this information to the text.

R2: Is the criterion of 10 times wavelength (line 86) given in this case?

A.: Yes, see above.

R2: How does the chosen filter (which? Please state per event?) influence the resolution?

A.: The influence of the filter has been tested carefully (lines 190-197; Figure S4). It turned out, that the influence of the choice of the cut-off frequencies is very small and entirely covered by the error of the standard deviation of the backazimuth. Therefore, the contribution of the choice of the filter may be neglected.

R2: I am missing a table with station information? Please add to Supp. Mat.

A.: As described in the data availability section, the data including all station locations are available for download at GEOFON (https://geofon.gfz-potsdam.de; https://doi.org/10.14470/4W7562667842).

R2: Where are the locations of all events (2709, line 154) comes from? Have they been located before? By whom? or from catalogue? And how much better is this multiple array beaming? Why can only be 112 of 2709 events be located with the new method? How much better do the locations become when using multi-array methods? This could be part of the discussion

A.: There seems to be a misunderstanding. The 2709 is the total number of events, which we recorded with our network on Fogo and Brava (line 154). Out of these 2709 earthquakes, we located 112 by applying the multi-array analysis. The first reason for this discrepancy is that most earthquakes are of relatively small magnitude and they  occur close to (or beneath) Brava. If they are too small, they are covered by noise at the stations on Fogo. As we operate only one array on Brava, this precludes the usage of multiple arrays. The next reason for this discrepancy are the requirements, which have to be fulfilled during the array analysis. If the result at one array is unstable, this array is not used for the analysis. This further reduces the number of events. (see lines 155-162)

Nevertheless, if the noise conditions are good (especially during nighttime), we are able to determine the epicenters of very small earthquakes (M_L < 0.5), where it can be hard to find the onset of the P-wave at all stations. This is one of the strong advantages of the  array analysis. If we combine multiple arrays, we are able to determine the epicenter of an event without the knowledge of a velocity model.

We will add this point in the revised manuscript.

R2: Line 149-151: how do you know the events are volcano-tectonic or hybrid? Have they been located and classified before? Where does that information come from? Reference?

A.: In contrast to a volcano-tectonic earthquake, a hybrid event shows a smooth transition from higher to lower frequencies. In the case of a hybrid event there is no S-phase identifiable and in the coda (for about 15 to 20 seconds) there is still more energy available in the 0-10 Hz band than before the event (and in comparison with a VT earthquake). In contrast, even for the small VT earthquakes on Fogo two phases are detectable. To make this point clearer, we will replace the example earthquake displayed in Figure 7a with an earthquake from Fogo. By this, the difference between the event types should become clearer. Additionally, we will add a figure to the supplementary material which shows the spectrogram of all components of a hybrid event. In this figure the low frequency content of the coda becomes evident.

A description of such events in previous literature on volcanic seismicity is e.g. found in McNutt (2000) or Wassermann (2012). The references are given in the introduction, lines 50-52.

With our discrimination between volcano-tectonic earthquakes and hybrid events we also follow earlier studies on the seismicity of Fogo (e.g. Faria and Fonseca, 2014, references in lines 255, 298).

R2: Also this part (4. Results) happens abruptly, without much leading into. Has the multi-array method been applied to known events? Or is it used on traces without knowing that there is an event? Please re-write. Much of the discussion seems to have to move to results section for a better flow and a better understanding of what has been done here.

A.: The traces were inspected by eye after application of a trigger mechanism. We will add this information to the results section in the revised manuscript.

As stated above, we will move the part on the error considerations of the multi-array analysis to the method section 3.

R2: Especially the event in Figure 6 is determined by 2 arrays, why do you need the third one there? How much better becomes the location if using three arrays? How much worse is the result when using classical locations?

A.: Some of the events are indeed located using two arrays, if the third array shows e.g. unfavorable noise conditions or if the event cannot be analyzed with the third array due to an unstable result during the beamforming. In such a case, the event will be located at the intersection of the main beams (i.e. by “drawing” a line from the backazimuth for which the amplitude of the sum trace reaches its maximum).

However, the more arrays are incorporated in the analysis, the more information are included, too. For example, Figure 13 shows, that the event is not located at the crossover of two main beams, but rather in between. This is the case, as the highest level of overlapping occurs in the range of the errors of the main beams. This is an advantage of the multi-array method, as the errors of the beamforming are directly taken into account during the localization of the event.

R2: In the introduction, the authors mention that the lack of S-waves makes it a good event for multi array analysis. I think it is more like the lack of S-waves makes it difficult for other location methods and that's why multi-array method can help?

A.: Yes, indeed. However, if the goal is to test the advantages of a multi-array approach it also makes sense to test it on events, which cannot be located reliably with classical localization methods.

We will rephrase this in the revised manuscript.

R2: The description in 3.2 is confusing, please clarify, refer to the example (e.g. Figure 6?) and explain what it means to intersect the broad beam in steps of 1%...

A.: We will clarify this by an additional figure added to the supplements.

R2: If the velocity model of Vales et al., 2014 is used for the classical location technique, why is this not used for the other technique? There a simple generic layer model seems to be used?

A.: The determination of the epicentral distance using the S-P travel time difference is an approximation to serve as reference for the analyst. The velocity model of Vales et al. (2014) is used here to determine Moho depth, mean crustal and mantle velocities. We will include this description in lines 229-230.

R2: Technical Points:

I may be wrong but at least one trace (4) in Figure 3a seems to be reversed in polarity. Have the traces been inspected for polarity changes? How much does the waveform complexity in this event influence the energy stacks?

A.: Trace 4 in Figure 3a is not reversed in polarity. Due to the red line overlaying the traces, it might appear like this. However, the red line is directly above the minimum of this amplitude. In Figure 3c it is better to see, that this trace shows the same polarity as the other traces. The difference between the trace 4 in Figure 3a and 3c is only a time shift, therefore the polarity does not change during this procedure.

In general, different polarities of the traces would make the beamforming impossible. It is based on the assumption, that the same phase at all traces is used for the shifting and stacking. The constructive interference of the phases leads to a reduced SNR of the sum trace. Different polarities, however, would lead to destructive interference.

We perform the array analysis manually, this means, if different polarities would occur, the analyst would recognize this. Therefore, we can assure to only perform the beamforming on the same phases. Additionally, we focus on the first arrival of the P-wave. This ensures, that the same phase is chosen at all arrays. This is important for the multi-array analysis to ensure, that the same phase is analyzed (lines 215-216).

R2: Line 166 needs a reference.

A.: This information is given in most seismology textbooks. Nevertheless, we will include the information, that this only holds true for rays propagating downward from the source.

R2: Line 81: I would add local to earthquakes. For teleseismic events this is not necessarily the case.

A.: Thank you for pointing this out. We will add this term in the revised manuscript.

R2: Figures: Please add the colour bars in slowness-backazimuth plots with relevant information.

A.: We will adjust the figures accordingly.

R2: Lines 109-111 seems to be irrelevant for the manuscript at this point? it does not fit into the flow and it knowledge that should be known.

A.: The angle of incidence is used further below to estimate the ray path of the waves of hybrid events. Therefore, it needs an introduction, which is best placed at the description of the parameters which result from the array analysis.