Angelo Pisconti
email: pisconti@uni-muenster.de

Katrin Plenkers

Joachim Philipp

Christine Thomas

Mapping lithological boundaries in mines with array seismology and in-situ Acoustic Emission monitoring

Keywords: Time series analysis, Body waves, Coda waves, Wave propagation, Wave scattering and diffraction

come

Introduction

In underground exploration the monitoring and analysis of seismic events are common tools for risk and hazard assessment, as underground operations induce low-magnitude seismicity and sometimes trigger rockbursts that could harm miners and damage the mine [START_REF] Hasegawa | Induced seismicity in mines in Canada-An overview[END_REF][START_REF] Potvin | Strategies and tactics to control seismic risks in mines[END_REF]. Since the beginning of last century seismic networks are used to study seismic events in mines [START_REF] Gane | A seismic investigation of the Witwatersrand earth tremors[END_REF]McGarr, 1971a,b;[START_REF] Spottiswoode | Perspectives on seismic and rockburst research in the South African gold mining industry: 1983-1987[END_REF][START_REF] Gibowicz | An Introduction to Mining Seismology[END_REF][START_REF] Collins | High resolution mechanics of a microearthquake sequence[END_REF][START_REF] Plenkers | Observation of seismic events with frequencies f > 25 kHz at Mponeng deep gold mine, South Africa[END_REF]. Microseismic and sometimes in-situ Acoustic Emission (AE) monitoring networks are operated in mines worldwide and provide large amounts of data (see [START_REF] Manthei | Review on In Situ Acoustic Emission Monitoring in the Context of Structural Health Monitoring in Mines[END_REF] for a review).

Until today the analysis is often limited to the estimation of the location, the magnitude and in some cases the faulting mechanism of the seismic events. Analysis of the coda wave field are rare, although the event's waveforms are accessible in modern monitoring networks and coda wave analysis can provide important additional information on structures and scattering properties of the propagation medium (see [START_REF] Sato | Earth heterogeneity and scattering effects on seismic waves[END_REF][START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth[END_REF] for a review of coda waves analysis).

Of great importance in underground excavation is the in-depth knowledge of the internal structures of the rock as this information is an important factor in the hazard assessment (e.g. [START_REF] Milev | Effect of the rock properties on mining-induced seismicity around the Ventersdorp Contact Reef, Witwatersrand Basin, South Africa[END_REF] since it influences our view of the stability of mine structures and directly reflects the production rate. Investigation of mine structure is also crucial in the safety assessment of underground repositories for radioactive nuclear waste. Some methods to study the rock's internal structure already exist, e.g., drilling combined with bore-core analysis, georadar, or seismic profiling. However, data gained by passive seismic monitoring systems are only exploited in rare cases for such purposes (e.g., Olivier et al. 2015a,b).

Gaining information on internal rock structures using already existing waveforms of seismic monitoring arrays is a cost-effective approach that has the potential to provide additional information on the existence and location of geological boundaries, discontinuities, or cavities. Moreover, energy of underground seismic events often travels through a rock volume that may be out of reach for other methods or not suitable for drilling e.g., due to rock integrity requirements.

To exploit seismic arrivals other than simply using first onsets in recorded waveforms, here we apply array seismological techniques. This approach is generally used to gain information on structures in the entire Earth from crust and lithosphere to the Earth's inner core and has provided valuable information on Earth structures (e.g., [START_REF] Gu | Arrays and Array Methods in Global Seismology[END_REF]. With the array techniques, we are able to extract coherent and incoherent signals hidden by noise through stacking, and those techniques also provide directional information of seismic arrivals [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF][START_REF] Gibbons | The Applicability of Incoherent Array Processing to IMS Seismic Arrays[END_REF][START_REF] Husebye | Array seismology-Past, present and future developments[END_REF][START_REF] Rost | Array seismology: methods and applications[END_REF][START_REF] Ringdal | P-Wave Envelope Representation in Event Detection Using Array Data[END_REF][START_REF] Schweitzer | Seismic Arrays[END_REF]. One method is the frequency-wavenumber (fk) analysis [START_REF] Capon | Signal processing and frequency-wavenumber spectrum analysis for a large aperture seismic array[END_REF] that helps to detect waves arriving from out-of-plane directions that can provide information on structures that do not lie on the direct wave path.

To our knowledge, these methods have not yet been applied to mine-scale investigations and picoseismicity, even though, despite the different scales and magnitudes, it is well documented that waveforms of picoseismicity events are well suited for extended waveforms analysis; e.g., source parameter analysis [START_REF] Dahm | Automated moment tensor inversion to estimate source mechanisms of hydraulically induced micro-seismicity in salt rock[END_REF][START_REF] Kwiatek | JAGUARS Research Group. Source Parameters of Picoseismicity Recorded at Mponeng Deep Gold Mine, South Africa: Implications for Scaling Relations[END_REF]. We test the application of these techniques to high-frequency AE picoseismicity that occurs in a frequency range from 1 kHz up to 200 kHz [START_REF] Manthei | Review on In Situ Acoustic Emission Monitoring in the Context of Structural Health Monitoring in Mines[END_REF] on a data set recorded by a sensor network in a salt mine.

We chose the Asse II salt mine in Lower Saxony where seismic events recorded by in-situ AE monitoring are reported [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF]. By analyzing and interpreting these waveforms, we address the following questions: Can seismic array methods be successfully applied to investigate AEs on minescale applications? Is the upper salt dome homogeneous or is there an indication for additional geological structures, presently unknown? Can we exclude the existence of cavities or large open fractures in the upper salt dome? Can we detect and confirm the location of known structures e.g., the southern flank of the salt dome using array seismological techniques?

Dataset, observations, method and data processing

We analyze single component waveforms from 52 events recorded at a network of 16 piezoelectric sensors, deployed above the A3 (western) chamber of the Asse II mine (Figure 1) in Germany. The events define a narrow cluster located in the upper part of the salt dome (Figure 1), with hypocentral distances ranging from 120 to 150 meters. The sensors were installed vertically in upward 1 to 3 m deep boreholes, drilled in the A3 chamber roof [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF], giving a 3D geometry to the network, with sensor locations following the convex roof shape (Figure 1). The network has dimensions of about 37 m x 31 m x 5 m (X x Y x Z, see inset in Fig. 1). The piezoelectric sensors with a sampling rate of 1 MHz allow the recording of very high frequency waves within a broad range from 1 kHz to 100 kHz. To reduce the amount of stored data, the sensors were equipped with a trigger system for the detection and the recordings of events. The recording time for each trigger is 32.768 ms [START_REF] Manthei | Review on In Situ Acoustic Emission Monitoring in the Context of Structural Health Monitoring in Mines[END_REF], allowing the recording of the P-wave and S-wave onsets, necessary for locating the events. A more detailed description of the events and sensors is found in Phillip et al. (2015).

A typical example of the waveforms of the event cluster used in this study is shown in Figure 2. This event occurred at a hypocentral distance of 122 m and shows clear direct P-and S-wave arrivals. The direct P-wave arrival is followed by additional energy packets (pulses) suggesting scattering or reflection of seismic energy somewhere in the propagation medium. Local reverberations close to the sensors can be excluded as possible cause for those observations, since the additional phase appears at all sensors indicating a more coherent structure. The additional signal also shows no clear move-out, when compared to the direct P-wave arrival, suggesting that it likely traveled out of plane and may have been scattered at a structure in the surrounding medium.

Scattering of seismic energy has been extensively used to characterize and probe heterogeneities at different scales, using both stochastic and deterministic approaches (see [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth[END_REF] for a review). Here we assume single scattering for this additional arrival since multiple-scattering processes would scatter the seismic energy continuously in a diffusive manner, therefore generating a smooth decay rate of the coda amplitude. Our observation of individual pulse-like signals, however, suggests the presence of one or more well-defined reflectors or scattering bodies, which can be located in a deterministic way provided that the direction of the incoming scattered wave-front can be measured.

Due to the intrinsic frequency dependency of scattering phenomena [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth[END_REF], we applied different filters before the beamforming process, in order to search for the frequency band which best enhances the amplitude of the scattered wavefield in the P-coda. In most cases, a butterworth bandpass filter of 3 to 8 kHz best reduced the high frequency reverberations while at the same time increased the amplitude of the additional arrival. A typical spectrogram of the analyzed recordings is shown in Figure S1 of the Supporting Information where we find the dominant frequency of the signal in the band between 3 and 8 kHz.

Seismic array techniques allow to measure the direction of the incoming seismic energy at a seismic network. A comprehensive description of array methods is found in [START_REF] Rost | Array seismology: methods and applications[END_REF] and [START_REF] Schweitzer | Seismic Arrays[END_REF]. Array methods are based on shifting and stacking of waveforms with the aim to increase the signal to noise ratio of coherent and incoherent seismic phases, and identify the direction of arrivals [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF][START_REF] Gibbons | The Applicability of Incoherent Array Processing to IMS Seismic Arrays[END_REF][START_REF] Husebye | Array seismology-Past, present and future developments[END_REF][START_REF] Ringdal | P-Wave Envelope Representation in Event Detection Using Array Data[END_REF][START_REF] Rost | Array seismology: methods and applications[END_REF].

The technique most often used for stacking seismic data is beamforming, which is based on the delay and sum method (e.g., [START_REF] Rost | Array seismology: methods and applications[END_REF][START_REF] Schweitzer | Seismic Arrays[END_REF]. This method applies time shifts (τj) to the traces depending on the slowness vector (s) of the incoming planar wave front and on the relative distance (rj) of the sensors with respect to the geometrical center of the array: for the j th -sensor. Expanding the previous equation, we obtain:

𝜏 𝑗 = -𝑥 𝑗 • sin 𝛷 -𝑦 𝑗 • cos 𝛷 𝑣 𝑎𝑝𝑝 + 𝑧 𝑗 • cos 𝑖 𝑣 𝑐 (2)
where the relative distance of the j th -sensor (with respect to the array center) is now expressed in terms of Cartesian coordinates (𝑥 𝑗 , 𝑦 𝑗 , 𝑧 𝑗 ) projected in the direction of wave front propagation along the backazimuth 𝛷 (Figure 3a). 𝑉 𝑎𝑝𝑝 is the apparent velocity of the approaching wave front which is related to the horizontal slowness (Figure 3b). The last term in the equation accounts for additional travel time shifts due to the 3D geometry of the network [START_REF] Schweitzer | Seismic Arrays[END_REF], where 𝑧 𝑗 is the elevation of the j th -sensor (the position in depth) and 𝑣 𝑐 is the velocity of the medium where sensors are located. In our case the velocity of both P-and S-waves is precisely known by ultrasonic transmission measurements performed on-site, using an ultrasonic signal in the same frequency range as the AE events and recorded on the sensors of the network. The ultrasonic transmitter was installed inside a borehole. The measured velocity for P-and S-waves are: 𝑣 𝑃 = 4570 m/s and 𝑣 𝑆 = 2597 m/s [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF]. Using the measured 𝑣 𝑃 and a depth difference of up to 5 m between sensors results in a (vertical) delay time 𝜏 of about 1 ms. This is a significant and relevant additional time shift to take into account in the stacking procedure. In contrast, in global array seismology where the elevation of the station is often much smaller than epicentral distances, these time shifts can often be neglected, however, even there static corrections due to strong topography are important in some cases (e.g., [START_REF] Bokelmann | Azimuth and slowness deviation from the GERESS regional array[END_REF][START_REF] Jacobeit | Influence of station topography and Moho depth on the mislocation vectors for the Kyrgyz Broadband Seismic Network (KNET)[END_REF]).

An approximation used in most array methods is that the wave front approaching the array is planar. This approximation works well at the large distances typical for global seismology but it might not hold for the short hypocentral distances considered in this study. Given the geometrical configuration of the array and its large aperture (~35 m), compared to hypocentral distances (120 -150 m), we only consider the 10 innermost sensor locations closest to the geometrical center. This sub-network has an aperture (~20 m), which is six times smaller than the minimum hypocentral distance (120 m), validating the approximation of planar wave-front [START_REF] Almendros | Array analysis using circular-wavefront geometry: An application to locate the nearby seismo-volcanic source[END_REF]. Relative distances of the sensors with respect to the geometrical center for this sub-network are shown in Figure 3c. This smaller configuration also improves the focusing of the arrivals at the remaining sensors. Due to the size of the source ruptures on the dm-to mm-scale in the considered frequency range [START_REF] Kwiatek | JAGUARS Research Group. Source Parameters of Picoseismicity Recorded at Mponeng Deep Gold Mine, South Africa: Implications for Scaling Relations[END_REF][START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF], near-field effects on the wave-front geometry are largely negligible.

In general, the capability of the arrays in enhancing the seismic signals depends on the interstation distance (i.e., the geometry of the network) compared with the wavelengths of the signals of interest [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF][START_REF] Gibbons | The Applicability of Incoherent Array Processing to IMS Seismic Arrays[END_REF][START_REF] Rost | Array seismology: methods and applications[END_REF], although the variation of signal shape across the array (i.e., its intrinsic incoherency) might still result in a reduction of stacked signal amplitude [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF] especially for high frequencies [START_REF] Ringdal | P-Wave Envelope Representation in Event Detection Using Array Data[END_REF]. In our particular case, considering the dominant high frequencies (3-8kHz) of the signal, the expected wavelength of the plane wave is about 0.6-1.50 m for 𝑣 𝑃 =4570 m/s. In conventional phase-based beamforming, the stacking of coherent signals with this wavelength would then require an interstation distance of the same order, which is significantly smaller than the actual average distance of the order of a few meters between the used sensors (Figure 3). Moreover, the stacked beam might be affected by the poor signal correlation between the sensors [START_REF] Husebye | Array seismology-Past, present and future developments[END_REF]. This problem in classical array processing arises when high frequency signals are recorded at large aperture arrays, and due to waveform dissimilarity between the sensors leads to beam degradation [START_REF] Gibbons | The Applicability of Incoherent Array Processing to IMS Seismic Arrays[END_REF]. However, to exploit high frequency energy, [START_REF] Ringdal | P-Wave Envelope Representation in Event Detection Using Array Data[END_REF] showed that these issues can be overcome by replacing the signals with their envelopes which leads to the so-called incoherent (envelope-based) beamforming [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF][START_REF] Gibbons | The Applicability of Incoherent Array Processing to IMS Seismic Arrays[END_REF][START_REF] Husebye | Array seismology-Past, present and future developments[END_REF][START_REF] Ringdal | P-Wave Envelope Representation in Event Detection Using Array Data[END_REF]. Therefore, using the Hilbert transform we compute the envelope of the filtered signals [START_REF] Kanasewich | Time Sequence Analysis in Geophysics[END_REF], which offers more insight on amplitude information of the seismic trace and enhances the correlation among signals, rather than looking at the complexity of the waveforms which also depends on the phase contained in the high frequency signal itself [START_REF] Farnbach | The complex envelope in seismic signal analysis[END_REF]. In Figure S2 of the Supporting Information we show the comparison between conventional (phase-based) coherent linear beamforming and incoherent (envelope-based) beamforming, where we find higher amplitudes of the stacked signal when using the envelope of the individual waveforms. Phaseindependent envelope-based stacking has also been used at global scale to detected scattering and heterogeneities in the mantle and the core (e.g. [START_REF] Bentham | Fine-scale structure of the mid-mantle characterised by global stacks of PP precursors[END_REF][START_REF] Peng | Inner-core fine-scale structure from scattered waves recorded by LASA[END_REF]Shearer ad Earle, 2004) As shown in Figure 2 and Figure S2, the envelopes of the signals consist of several wave-packet arrivals with a time length of about 1 ms, which is more coherent across our large aperture array with 4-5 meters of interstation distances. However, the use of waveform envelopes might lower the resolution and the resulting ability of the array to locate the origin of signals in the propagation medium.

For the beamforming procedure (following Rost and Thomas, 2001;[START_REF] Schweitzer | Seismic Arrays[END_REF] the stacking is performed on the shifted envelope of the traces 𝑤 𝑗 recorded at the sub-array, as shown in Figure S2 of Supporting Information, based on the following equation:

𝑏(𝑡) = [ 1 𝑀 ∑ √𝑤 𝑗 (𝑡 + 𝜏 𝑗 (𝒔, 𝛷)) 1 𝑛 𝑀 𝑗=1 ] 𝑛 (3)
where M is the number of sensors. The n-root of the single traces reduces the difference in amplitudes, allowing coherent phases with small amplitudes to have a bigger influence on the stacked trace compared with the linear stacking (n = 1), and applying the power of n to the beam restores the amplitude differences, reducing incoherent noise after stacking. The search for out-of-plane signals is implemented using the slowness-backazimuth analysis, where a grid search is carried out by computing the energy of the beams for different slowness-backazimuth values.

Results

We apply the slowness-backazimuth analysis to the array recordings of the 52 events. To explore the directionality of different energy packets arriving at the network, the analysis is carried out using a sliding time-window with a duration of 4 ms. This window is moved along the entire recording with steps of half the size of the window (i.e., 2 ms). Due to the reduced length of ~32 ms for all available recording, we restrict our analysis to the P-coda, since the S-coda is sometimes abruptly interrupted by the end of the recording.

We present the results of the slowness-backazimuth analysis as energy contour plots on polar diagrams, where the direction of the energy arriving at the network can be picked (Figure 4). We clearly identify the direct P-wave for all 52 events. As shown in Figure 4a,c, the direct P-wave arrival reaches the network from the north-east, with the expected theoretical values of slowness and backazimuth. The theoretical values were computed from the known position of the events and the array center, and the ray parameter formula for a homogeneous propagation medium (Udias, 1999) with 𝑣 𝑃 = 4570 m/s. For the 52 events analyzed, the average discrepancy between observed and theoretical slowness values for the direct Pwave is of the order of 0.0073 s/km which corresponds to about 2 degrees of difference in incidence angle between observed and theoretical direct P-wave, while the average residual in backazimuth is of the order of 3 degrees. These small residuals lie within the average standard deviation σ, which is representing the uncertainties in measuring the observed slowness (σ=0.0075 s/km) and backazimuth (σ=4°) values of the direct P-waves. To measure these uncertainties and estimate standard deviation, we consider all points that have a variation of up to 10% in beam power with respect to the maximum (white contour line in the polar plot of Figure 4c) as performed also by [START_REF] Schumacher | Out-of-plane seismic reflections beneath the Pacific and their geophysical implications[END_REF].

The observations of the P-wave arrivals and the match between observed and theoretical slowness values suggests the presence of a homogeneous medium between sources and receivers, at least at the short wavelength signals (~60 cm) we used in our analysis. No other arrivals that deviate from the source to receiver trajectory appear in the time window centered around the P-wave, implying that the early Pcoda (1 -3 ms after the main P-wave arrival) mainly consists of direct energy radiated from the source region.

In many cases the main P-and S-phases are immediately followed after a few milliseconds by a second peak (Figure 2). Considering the distance (1 -3 m) between the sensors and the underlying roof of the A3 chamber and 𝑣 𝑃 = 4570 m/s and 𝑣 𝑆 = 2597m/s, these second peaks are most likely reflections from the top of the chamber, which acts as a total reflector since the chamber is hollow. Indeed, calculation of the two-way-travel time (TWT) yields roughly 0.4 -1.3 ms for the P-wave reflections and 0.8 -2.3 ms for the S-wave reflections for a vertical reflection. These near receiver reflections have roughly the same slowness as the direct P-waves, but they arrive at the network from below. Due to the convex shape of the chamber roof, defocusing of this reflected energy is expected to broaden the reflected P-wave peak in the slowness backazimuth plots, as observed in the relevant time window in Figure 4c. In addition, the use of the envelope in the stacking procedure might cause lower resolution of the slowness and backazimuth values of the detected signals and consequently produce broader peaks. Despite this, the aforementioned small uncertainties in the estimation of both observed slowness and backazimuth suggest a very good capability of the envelope-based stacking method in detecting the travelling direction of the waves.

In the later P-coda, we consistently observe secondary arrivals with propagation directions mainly from the south-east (Figure 4c, middle) and observed backazimuth values ranging from 94° to 200° measured from North. 17 out of 52 events show these late arrivals in the P-coda. These out-of-plane reflections suggest a strong impedance contrast able to reflect/scatter the seismic energy within a wide range of incidence angles. The location of the events showing out-of-plane signals can be found throughout the event cluster, suggesting that these observations do not depend on event location i.e., it is not a local effect.

On average these additional arrivals have a similar lapse time of about 6 (± 0.5) ms after the P-wave and their energy is about half (0.5 ± 0.16) of that of the main P-wave arrivals. Moreover, the additional signals are also characterized by a single energy peak in the energy contour plot (Figure 4c), therefore suggesting that a reflector is responsible for these out of plane observations, rather than multiple-scattered energy, which would appear as a series of peaks in the polar plots and a longer delay of the coda waves [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth[END_REF]. The small variance in the measured directionality (slowness and backazimuth)

and travel time suggests a common source region for the scattered wavefield. No other dominant arrivals are visible in the contour plots of the slowness-backazimuth space at other times (Figure 4), implying that the medium in the vicinity of the monitored rock volume is simple with no other scatterers/reflectors influencing the data in our filtered frequency range.

To locate the position of scatterers/reflector responsible for the observed out-of-plane signals, we use the measured slowness, backazimuth and relative travel time with respect to the direct P-wave. Particularly the slowness and the backazimuth constrain the travelling direction of the out-of-plane arrival while the travel time constrains the position of the scattering body along the ray-path direction (Figure 5). This approach has been used in global seismology for the detection and location of heterogeneities in the Earth mantle that cause scattering and/or reflection of seismic waves (e.g. [START_REF] Kaneshima | Seismic scatterers at the shallowest lower mantle beneath subducted slabs[END_REF][START_REF] Kaneshima | Seismic scatterers in the mid-lower mantle[END_REF][START_REF] Kaneshima | Subparallel dipping heterogeneities in the mid-lower mantle[END_REF][START_REF] Rost | Seismic array detection of subducted oceanic crust in the lower mantle[END_REF][START_REF] Schumacher | Detecting lower-mantle slabs beneath Asia and the Aleutians[END_REF][START_REF] Schumacher | Out-of-plane seismic reflections beneath the Pacific and their geophysical implications[END_REF][START_REF] Weber | Reflections from a distant subduction zone[END_REF][START_REF] Weber | Seismic evidence for a steeply dipping reflector-stagnant slab in the mantle transition zone[END_REF][START_REF] Wicks | Seismic evidence for a fossil subduction zone beneath the Philippine Plate[END_REF][START_REF] Wright | Array studies of seismic waves arriving between P and PP in the distance range 90° to 115°[END_REF][START_REF] Wright | Longitudinal waves from the Novaya Zemlya Nuclear Explosion of October 27, 1966, recorded at the Warramunga seismic array[END_REF].

We run a grid-search procedure aimed at finding the location of the scattering points which best fit the observed travel times, given the observed slowness and backazimuth values (Figure 5). In our backtracing algorithm we assume that both primary and scattered waves travel in a homogeneous medium which is supported by the undisturbed direct P-waves travel paths in the rock volume between the sources and the receivers, for the entire azimuthal coverage available in our data. As stated before, we assume that the observed out of plane signals are due to single P-to-P scattering. P-to-S scattering would have larger travel times, unless the scatterers are located very close to the network but we can exclude this case due to the observed straight ray-paths of the direct P-waves (Figure 4c). Due to the smaller value of the S-wave velocity compared with the P-wave velocity, hence a different move-out, P-to-S scattered waves would also have larger slowness values than those observed for our out-of-plane waves (Figure 4c). We also exclude the possibility of S-to-P scattering because the observed travel times can only be matched by near source scattering, which we can exclude due to the large backazimuth deviations of the out-of-plane waves that are inconsistent with source-side scattering.

Carrying out the search for the reflector locations we find that, as expected by the small variance of the measured values of slowness, backazimuth and travel times, the reflection locations are very close to each other, defining a narrow cluster (Figure 6). As stated above and following [START_REF] Schumacher | Out-of-plane seismic reflections beneath the Pacific and their geophysical implications[END_REF],

we evaluate the uncertainties in estimating the scatter position by considering all points within the area of 90% of the maximum beam power, as shown by the white contour line in the slowness-backazimuth polar plots of Figure 4c and calculate the standard deviation of the maximum slowness and backazimuth values, which yields the angular error on the scatter position. The errors in travel time (i.e., the width of the peak) of the out-of-plane signals were picked as shown for one example in Figure 4c. On average, the travel time uncertainties are of the order of 1ms, which roughly corresponds to a location uncertainty of about 4-5 meters. By using the back-tracing algorithm, the uncertainties in slowness, backazimuth and travel time were then translated into scattering volumes, that yield ellipses of a few meters which we project on the geological map and section of Figure 6.

The cluster of the back-projected reflection points delineates a south-west dipping steep (~70° from horizontal) reflector located at a distance of about 70 m south-east of the network above the A3 chamber (yellow ellipses in Figure 6). It has a depth extension of about 80 m and is about 60 m long, with a strike direction about WNW/ESE. There are also 4 scatter locations which are further away from the cluster (green ellipses in Figure 6). However, their origin is not clear as they do not coincide with any known boundary in the geological section.

Discussion

By using seismic array analysis, originally developed for nuclear test monitoring and to address global seismological problems (e.g. [START_REF] Douglas | Seismometer arrays -Their use in earthquake and test ban seismology[END_REF][START_REF] Husebye | Array seismology-Past, present and future developments[END_REF][START_REF] Rost | Array seismology: methods and applications[END_REF], we were able to locate the origin of the out-of-plane arrivals in the late P-coda of our in-situ acoustic emission observations from 52 events recorded by a network of sensors deployed at the Asse II mine. To our knowledge, the use of array analysis has never been tested on AE in-situ recordings at the small scale of the mines with hypocentral distances on the order of about 100 m. While the principle of the waveform stacking remains the same, there are some significant differences in the processing scheme that we developed, compared with the application of array analysis at global seismological scale: 1) Due to smaller source dimension (mm to dm scale) and because the high frequencies are damped more quickly, in typical AE monitoring the hypocentral distances are usually within a radius of less than 200 m; 2)

Given the typical dimension of the network of a few tens of meters (see [START_REF] Manthei | Review on In Situ Acoustic Emission Monitoring in the Context of Structural Health Monitoring in Mines[END_REF] for a review of AE monitoring), the approximation of a planar wavefront, that is used for most array methods, might not hold. As a rough guideline, the curvature of the wavefront has an impact on the delay time estimation when the source-receiver distance is smaller than about 4-5 times the array aperture [START_REF] Almendros | Array analysis using circular-wavefront geometry: An application to locate the nearby seismo-volcanic source[END_REF]. By excluding the outermost sensors of our network, we restrict the dimension of the used sub-array to be approximately one-fifth of the minimum hypocentral distance in our dataset;

3) The 3D geometry of the network had to be considered in the calculation of the delay time as given by equation ( 2), including the vertical distances between the sensors, while in global seismology the elevation of the stations is often negligible compared to their horizontal relative position of stations. Note also that, in contrast to standard seismological problems, in this experiment the events occur at shallower depths than the receivers with the seismic energy arriving at the array from the top rather than from below; 4) Given the very high frequency content of the AE (of the order of kHz) and due to the large aperture of the array with sparse set of inter-sensor distances, incoherent (envelope-based) beamforming is used to increase the correlation between signals recorded across the used sub-network, in contrast to global seismology where the lower frequencies and the appropriate aperture array (compared also to epicentral distances) allow the use of the classic (phase-based) array analysis.

Our observations of direct P-waves travel times agree with travel times calculated from straight ray-paths from the sources to the receivers, therefore serve also as an independent constraint and validation of the event locations provided by [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF]. Since we use seismic array methods and in particular the slowness-backazimuth analysis, we constrain the directional information of the seismic energy arriving at the network. Figure 4 shows that the observed backazimuth and slowness matches the theoretical values for a homogeneous medium (𝑣 𝑃 = 4570 m/s) computed for each event location from [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF]. If the medium between source and receiver was heterogeneous, the measured slowness and backazimuth would differ from the theoretical values, because the travel path would be bent. Secondly, in a homogeneous medium the backazimuth depends only on the horizontal (that is geographical) position of the source with respect to the sensors and the slowness is related to the moveout of the P-wave at the sensors (i.e., the incidence angle for a given vp), which constrains the direction of the source with respect to the sensors in the three-dimensional space. Then, knowing the differential travel time between the first arrivals of the P-and S-wave from [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF], allows an independent estimation of source location compared with the inversion used by [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF].

Using the observations of straight paths of the direct P-waves we conclude than that the rock volume between events and receivers in the upper salt volume does not contain structures, such as geological heterogeneities or damage zones that are larger than 60 cm. Structures smaller than 60 cm that might be present are not visible in this analysis and cannot be excluded, although we also filtered the data with different frequency bands compared to the one used in the main analysis (3-8 kHz) and find that the direct P-wave still travels relatively undisturbed through the medium. This suggests that there are likely no additional impedance contrasts within the rock body sampled by the direct P-wave, which is in agreement with the geological model of the salt structure (Figure 6). Confirming a homogeneous and undisturbed rock volume in this particular area is an additional valuable information for the mine operator, as it may provide indications on the integrity of rock above the mine's structure, in this case the two chambers.

To verify the robustness of our assumption that the out of plane arrivals are due to P-to-P reflection/scattering, rather than S-to-P or P-to-S conversion phenomena, we carried out the analysis for the latter two cases. We find that it is not possible to fit the travel times using S-to-P reflections, since they show a minimum residual between observed and computed travel times of about 10 ms (Figure S3), much larger than the average travel time difference between the observed out-of-plane wave and the direct P-wave. In contrast, using the P-to-P reflections we have residual travel time of about 0.01 ms.

To test the case of P-to-S reflections, we re-run the slowness-backazimuth analysis using 𝑣 𝑐 = 𝑣 𝑆 in equation ( 2) to account for the different move-out of the S waves. Using 𝑣 𝑆 we find that the S-wave direct energy is focused better while the direct P wave defocusses. By back-tracing the out-of-plane arrivals as potential P-to-S reflections we find that P-to-S reflection points would be located in a tight cluster directly above the sensors with subvertical incidence to the array, as shown in Figure S4 of the Supporting Information. This region is well sampled by the direct waves (Figure S4) and, as previously demonstrated, P-waves travel undisturbed with straight ray-paths from sources to receivers which led us to the conclusion that the rock volume above the sensor appears to be homogeneous (Figure 6). Since the takeoff angles and azimuths of P waves and P-to-S scattered waves near the sensors are very similar (Figure S5), a heterogeneous medium should affect both wave types and we can therefore exclude strong heterogeneity near the sensor array above the A3 chamber. In contrast, the P-to-P reflections leave the source region in a different direction compared with the P-wave (Figure S5), with larger take-off angles and a range of azimuths, allowing the sampling of a larger rock volume outside the source and receiver region. In addition, vertically installed piezoelectric sensors have a minimum sensitivity to vertically travelling S waves, while it is maximum for horizontally travelling S waves. Since P-to-S scattered waves would have a small incidence angle (see above) hence arrive almost vertically, we would not be able to detect them as high-energy arrivals as those that we find in our dataset. Finally, the comparison with the geological section and map (Figure S4) shows that the P-S reflection point locations do not coincide with any known major geological discontinuity or structure within the Rock salt (z3) of the Leine Zechstein sequence. We therefore believe that P-to-P reflections are a robust and valid explanation for the observed out of plane arrivals.

It is unfortunate that the length of the recordings is limited in time by the event triggering which does not allow us to use and analyze the S-coda as an independent proof of the imaged reflector as shown by [START_REF] Weber | Reflections from a distant subduction zone[END_REF] and [START_REF] Weber | Seismic evidence for a steeply dipping reflector-stagnant slab in the mantle transition zone[END_REF] for the detections of heterogeneities at the larger scale of the Earth's mantle. S-to-S reflections would arrive at travel time of at least 35 ms, while our recording last 32.768 ms.

The frequency band (3-8 kHz) used in our analysis places some constraint on the gradient of the impedance contrast across the reflector that is responsible for the observed P-to-P out-of-plane signals.

Considering a 𝑣 𝑃 = 4570 m/s and the highest frequency content of the signal, the impedance contrast must vary within a wavelength of about 60 cm i.e., a sharp reflector. This agrees well with the geological cross section and map of Figure 6, which show that the reflection points are located closely to an abrupt lithological boundary at the southern flank of the salt dome.

The position of our calculated reflection points match the position of the southern flank of the salt dome at depth between approximately -180 and -260 m (below mean sea level, m.s.l.), and define the orientation (WNW-ESE) and dipping (~70° from horizontal) of the contact between the salt dome and the host rock, as shown by the yellow ellipses in Figure 6. This part of the salt dome is characterized by the presence of thin, steep, subparallel layers marking the boundary of the main salt body, where the Permian Rock salt dome from the Leine (z3) and Aller (z4) Zechstein sequence is in contact with the lower Triassic Röt anhydrite of the Buntsandstein formation. The green ellipses in Figure 6 represent scatterers whose location is far from the main cluster of reflection points (yellow ellipses in Figure 6)

and away from any lithological boundary in the geological section. However, the geological map shows that these location points might relate to a minor lithological contrast in the salt dome, where the Anhydrite and salt layers (z3) formation bends and extends over the Rock salt of the Leine sequence (z3)

(see map in Fig. 6).

Based on the geological information reported in Figure 6, the topography of the detected reflector is not expected to be rough. This allows us to calculate reflection coefficients assuming a flat discontinuity in order to compare them with measured P-to-P/P amplitude ratios. In computing the measured amplitude ratios we assume that the direct P-wave amplitude represents an estimate of the amplitude of the incident wave at the reflector. However, differences in take-off angle and azimuth at the source between the observed direct P-wave and the P-to-P reflection (Figure S5) might affect the P-to-P/P amplitude ratio.

Despite this, we measure linear (n=1 in equation ( 3)) beam amplitude ratios of the P-to-P reflections and the direct P-waves and we compare them with the P-to-P reflection coefficient computed for an impedance contrast using the 𝑣 𝑃 , 𝑣 𝑆 , and density values for the lithologies across the salt dome edge (Rock salt and Bundsandstein), as provided by the Asse II mine operators and as found in [START_REF] Karp | Vorplanung der 3D-seismischen Messungen zur Erkundung der Deckgebirgsstruktur, insbesondere der Querstörungen der Schachtanlage Asse II[END_REF]. We show the result of this comparison in Figure S6 of the Supporting Information. For completeness we also computed theoretical P-to-S reflection coefficient values, although, as shown above, we can exclude this possibility based on travel times. Overall the P-to-P reflection coefficient agrees well with the observed amplitude ratios with discrepancies most likely due to source directivity.

It is possible that our out-of-plane waves result from a change in the rheology of the rocks involved, since the location of the reflector agrees with the presence of a narrow branch of seismic activity along the southern flank of the salt dome, at the depth and location where we find our reflections [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF]. In particular, cluster C1 of seismic activity, as reported in [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF], is confined to the edge of the salt rock and does not extend into the host rock. In this study, we do not attempt to discriminate between different causes for the observed reflections based on amplitude, since amplitudes of seismic waves depend on many factors, such as source directivity, radiation pattern, and anisotropy which leads to magnification or reduction of the reflection coefficients in certain travelling directions (e.g., [START_REF] Pisconti | Discriminating between causes of D" anisotropy using reflections and splitting measurements for a single path[END_REF]. Source directivity and radiation pattern might also be a cause for the lack of clear out of plane signals in several events of our dataset, which leaves us with only 17 out of 52 events that show the additional arrival. Furthermore, the true amplitude of the waveforms are difficult to estimate because the instrument response of piezoelectric AE sensors is partially unknown. Accurate waveform modelling would be required to distinguish between these causes, which is beyond the scope of this study but might be of interest in future studies.

Conclusions

In this work we test the application of seismic array methods on in-situ AE recording from the Asse II salt mine in Lower Saxony (Germany). We use recordings of 52 events from a network of 16 piezoelectric sensors deployed above one of the chamber at Asse II. Hypocentral distances range from 120 to 150 m.

Applying a technique developed for monitoring the nuclear test ban treaty and for global seismology to in-mine AE recordings requires some adjustments and complementary processing. To enhance the correlation between the signals in the stacking procedure used in the application of array analysis: (1) we compute the envelope of the bandpass-filtered data (3 -8 kHz), which gives more emphasis to the amplitude of the signal itself rather than its phase;

(2) we reduced the size of the network (~20 m) to be at least one fifth of the smaller hypocentral distances, which allows us to perform the stacking based on a planar wave-front approximation;

(3) we need to take into account that sensors are located at different depth levels.

By applying the slowness-backazimuth analysis we found that all 52 events show direct P-wave arrivals that consistently appear with the expected slowness and backazimuth values for a homogeneous medium, indicating straight ray-paths between source and receivers. This confirms both homogeneity of the sampled rock volume (at wavelength of ~60 cm) between sources and receivers and agrees with the geological maps, and event locations as shown by [START_REF] Philipp | On the potential of In-Situ Acoustic Emission (AE) technology for the monitoring of dynamic processes in salt mines[END_REF].

We also find high amplitude pulse-like additional arrivals in the late P-coda of our AE recordings in 17 out of 52 events. Using seismic array methods we confirm these arrivals to be out-of-plane reflections travelling to the network from the south-east. Assuming P-to-P single scattering and using the measured backazimuth, slowness and travel time we back-project the data and find a narrow cluster of reflection points delineating a SW dipping (~70° from horizontal) steep reflector extending for about 80 meters This study proves that array seismology works at the small scale of mines and can be a valuable additional technique to detect geological heterogeneities and structures underground. To our knowledge this is the first application of array seismology techniques to mines and we find that this additional information could potentially support efforts to reduce hazard risks and might lead to more efficient mining operations and monitoring of underground nuclear waste repository mines. Large amounts of data are often collected in standard routines underground (microseismic monitoring / in-situ AE monitoring) in the framework of structural health monitoring and are available without deploying more instruments. These existing data are able to provide significantly more information on a bigger rock volume, than currently exploited in mine monitoring for hazard and risk purposes using both indirect and direct methods (i.e. drillings).

The analysis routine proposed in this work can be further tested and applied to other AE signals in different mines setting, provided that the size of network is at least one fifth of the hypocentral distances.

The automatization of such technique and its application to large datasets including also the S-coda analysis by increasing the recording's triggered time length, might provide an innovative and detailed passive 3D image of the mine structure in future work.

(https://www.crewes.org/ResearchLinks/FreeSoftware/). Some of the figures were created using Gnuplot (Williams and Kelley, 2011). AP analysed and interpreted the data, extended the methods to acoustic emission data and wrote the manuscript. KP and JP provided the data and information on acoustic emissions and salt mine settings and KP co-supervised the project. CT supervised the project and helped with array processing methods and the writing of the manuscript. Figures S1 to S6

Introduction

The supporting information includes additional explanatory figures that are complementary to the research reported in the main document. In particular: Figure S1 shows a spectrogram displaying the dominant very high frequency nature of a typical AE recordings analyzed in this work; Figure S2 shows delay and sum (beamforming) processing on both the original filtered waveforms and their envelopes, justifying the use of incoherent (envelope based) stacking in our analysis; Figure S3 shows residual travel times (between observational and theoretical out of plane arrivals) obtained when interpreting the observed out of plane arrivals as S-to-P reflections and Figure S4 illustrates the locations of potential P-to-S reflections clustered in the homogeneous medium directly above the sensors, which is also sampled by the undisturbed direct P-waves; Figure S5 shows take-off and azimuth angles at the source for P, P-to-P and P-to-S rays and Figure S6 illustrates P-to-P (and P-to-S) theoretical reflection coefficients computed using the expected impedance contrast across the lithological boundary at the southeastern edge of the main salt body, that we sample with our observed out of plane arrivals, whose amplitude ratio with the direct P wave is also shown.

(

  from about -180 to -260 m below m.s.l.) in depth and 60 m in length. When compared with the available geological information, this reflector matches a known lithological boundary, which has a strike direction WNW/ESE and pertains to the transition from the Permian main salt body (Rock salt from the Leine (z3) and Aller (z4) Zechstein stratigraphic sequence) to the lower Triassic Röt anhydrite of the Buntsandstein formation. The reflector provides an estimation of the position and inclination of the southern flank of the salt dome at a distance of about 70 m from the network.

(

  black triangles) used in this study. The geological cross section shows the main lithological formation with the upper part of the salt dome consisting of the Permian Zechstein evaporitic cycle (z2-z7) and the early Triassic Buntsandstein (map: BGE, Asse II mine, operational state after mining plans from April 1999, last updated on 22 nd of February 2016; section: BGE, Asse II mine, operational state after mining plans from January 2012, last updated on 28 th March 2012). The inset on the bottom left shows the detailed 3D locations of the piezoelectric sensors (black triangles) placed in boreholes (blue lines) on the roof of the A3 chamber of the Asse II mine, whose shape is represented by the gray lines. Depth are in meters below mean sea level (m.s.l.). Record section (distance versus time) of a typical in-situ acoustic emission array recording of a small event that occurred in the upper salt dome of the Asse II mine. As in (a), waveforms were filtered using a bandpass filter of 3 -8 kHz. Positive and negative amplitude lobes of waveforms are filled using yellow and black colors, respectively, while the envelopes are shown as blue lines. Both direct P-and S-waves clearly show their different move-out at the 16 sensors, as marked by the solid red lines. The red box with a question mark indicates additional pulse-like arrivals in the late P-coda.

Figure 3 .

 3 Figure 3. (a) Network configuration, as seen from above looking down into the Earth, and relative position of the sensors with respect to the geometrical center of the array. x15 and y15 mark the distance from the center of the array to sensor 15. A planar wave-front crossing the receivers with backazimuth 𝛷 from North is indicated by dashed lines. (b) Side view of sensor locations projected onto the sagittal plane containing the incident ray and crossing the array with apparent velocity 𝑣 𝑎𝑝𝑝 = 𝑣 𝑐 𝑠𝑖𝑛 𝑖 ⁄ , where 𝑖 is the incidence angle and 𝑣 𝑐 is the propagation velocity of the medium. (c) Euclidean distance (𝑟 𝑗 ) between sensors and array center showing the sensors (in blue) used in the analysis, within a distance range of maximum 15 meters, and the excluded sensors (in red) clustered far away from the network center.

Figure 4 .

 4 Figure 4. (a) source (red circles) to receiver (triangles) configurations and ray trajectory (black line) for the event shown in Figure 2. (b) Record section of the waveform envelopes for this event. The envelopes considered in the analysis are shown in blue, while the discarded envelopes are shown in red (see text and Figure 3). Grey areas denote three selected time windows used for analysis. (c) Polar plots of slowness backazimuth contours for the beam power for the three time windows shown in (b) indicating the incoming direction of the energy as function of backazimuth (baz) and slowness (u) at the A3 network (displayed as reference at the center of each polar plot).The slowness (or the incidence angle at the array, inc) is increasing along the radial direction, while the backazimuth is represented as angle clockwise from North. The beam power is normalized to the power of the direct P-wave energy. Window 1 shows a strong direct P-wave arrival with slowness and backazimuth in the direction of the source. Window 2 shows a dominant out of plane secondary arrival with energy maximum coming from the south-east and traveling with a difference of 80° to the direct P-wave (shown by the solid black circle). This window encompasses the secondary arrival shown in the single recordings (b) and Figure2. Window 3 shows no strong energy arrivals in the late P-coda. The white contour line in the polar plots of time windows 1 and 2 indicates the 0.9 isoline of the maximum for each arrival in the normalized beam power, from which we computed the uncertainties in slowness and backazimuth. Below each polar plot, the relevant beams are shown as computed using the slowness and backazimuth values of the energy maximum. Shown is also the uncertainty on the travel time estimation (dt) for the out of plane arrival of window 2.

Figure 5 .

 5 Figure 5. 3D schematic illustration of the back-tracing algorithm used to find scatter locations. The red star indicates the event, the yellow circle indicates the scatter location and the black triangles represent the network.The ray path direction at the receiver is drawn given the measured backazimuth (baz) from North and incidence angle inc (from zenith).

Figure 6 .

 6 Figure 6. Geological map at a depth of 200 m below m.s.l. (left) and section (right) of the upper part of the salt dome (map: BGE, Asse II mine, operational state after mining plans from April 1999, last updated on 22 nd of February 2016; section: BGE, Asse II mine, operational state after mining plans from January 2012, last updated on 28 th March 2012). Location of the events (red circles), sensors (black triangles) and scattering volumes (yellow and green ellipses) are also indicated. The four green ellipses refer to scatterers located far from the lithological boundary that mark the southern flank of the salt dome. The ellipses represent the scattering volumes as computed from the uncertainties of slowness, backazimuth and travel time (see text for more information). The straight grey lines show the direct P-wave raypaths from the events to the array center, while the white segmented lines indicate the incident and reflected raypaths of the out-of-plane observations.

Acknowledgements

We thank the geophysical working group of Lutz Teichmann of the Asse salt mine (Bundesgesellschaft für Endlagerung (BGE), Schachtanlage Asse II -Federal company for radioactive waste disposal Germany, mine shaft Asse II) for providing the data and fruitful discussions, specifically Harold Kühn and Denise De Nil. We would like to thank the editor Huajian Yao and the reviewers Michael Weber and Keith Koper for their comments and suggestions that greatly improved this research work. We thank James Wookey for useful comments and suggestions that greatly improved the quality of the manuscript.

The work of this study was funded by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 642029 -ITN CREEP. The data were processed using SAC (https://www.iris.edu/hq/), Seismic Handler (Stammler, 1993) and ObsPy, a python toolbox for seismology (Krischer et al., 2015). Theoretical reflection coefficients were computed using the Zoeppritz CREWES Matlab Toolbox (Margrave and Lamoureux, 2019) available at

Supporting Information for

Mapping lithological boundaries in mines with array seismology and in-situ

Acoustic Emission monitoring ). Velocities and densities were provided by the Asse II mine operators and can be found in [START_REF] Karp | Vorplanung der 3D-seismischen Messungen zur Erkundung der Deckgebirgsstruktur, insbesondere der Querstörungen der Schachtanlage Asse II[END_REF]. 𝑣 𝑆 for the Bundsandstein formation was estimated using the Poisson's solid approximation. The observed P-to-P/P amplitude ratios are color coded by their difference in take-off angle at the source, showing that possible discrepancies with theoretical values might be due to source geometry. Only reflection points located closely to the salt dome edge of Figure 6 are shown. Reflection coefficients were computed using the Zoeppritz CREWES Matlab Toolbox [START_REF] Margrave | Numerical Methods of Exploration Seismology: With Algorithms in MATLAB[END_REF].