3-D fluid channel location from noise tremors using matched field processing

SUMMARY

Presently ongoing geodynamic processes within the intracontinental lithospheric mantle give rise to different natural phenomena in the NW Bohemia/Vogtland region (Czech Republic, Germany), amongst others: earthquake swarms, mineral springs and degassing zones of mantle-derived fluids as well as highly concentrated CO₂ (mofettes). Their interaction mechanisms and relations are not yet fully understood, but fluid pathways within the crust are assumed, that allow efficient fluid transport between the main hypocentral swarm quake region and the degassing areas at the surface. Here, we focus on the location of the presumed fluid channels as well as on the investigation of their near-surface spatio-temporal variability, targeting a depth of a few hundreds of metmetres. We applied a 3-D matched field processing (MFP) approach in the frequency band of 10–20 Hz considering the fluid flow as seismic noise source. Within three campaigns in 2015/2016, we recorded continuous seismic noise data on the Hartoušov Mofette Field within the Cheb Basin (NW Bohemia, CZ), which is a key site to study fluid flow as it is characterized by strong and continuous surface degassing of CO₂. We used temporary arrays varying in extent (70-600 m aperture) and in the amount of stations (25–95 units). Assuming a homogeneous velocity model and applying conventional MFP phase-matching over a 3-D grid search, we located two channel-like structures beneath the test site, which could be traced down to a common source area down to 2000 m depth. We thereby evaluated the influence of amplitude normalization of the measured noise signal on the MFP location considering water-filled or dry mofette channels. Additionally, a spatio-temporal analysis using time windows with a length of 10 min during 5 hr of noise record shows variability of fluid flow activity in space and time and hence, its migration beneath the test site on a short timescale.

Numerical modelling, Time-series analysis, Seismic noise

1 INTRODUCTION

The NW Bohemia/Vogtland region (Czech Republic, Germany, Fig. 1) is an intracontinental and geodynamically active area, which is characterized by the occurence of different natural phenomena such as earthquake swarms and extensive surface degassing of mantle-derived fluids. The area is situated in the western part of the Eger Rift, which was generated due to the interaction of the Post-Alpine crustal extension with alkaline intraplate volcanism and magmatic activity during the Cenozoic. This active rift with an extent of 300 x 50 km² belongs to the European Cenozoic Rift System (Ziegler 1992). There are multiple basins filled with Cenozoic deposits that form its landscape, as for example the Cheb Basin, which is located in the western part of the Eger Rift. This shallow intracontinental basin was formed during the Neogene due to the reactivation of former Hercynian faults and is filled with about 300-m-thick Quaternary and Tertiary deposits. It developed at the intersection of the Eger Rift with the NNW–SSE striking Mariánské Lázne Fault Zone (MLFZ) and the nearly N–S trending Počatky-Plesná Fault Zone (PPFZ, Špičáková 2000; Bankwitz et al. 2003; Ulrych et al. 2011). The MLFZ bounds the Cheb Basin to the east, forming a steep slope with a height of 50–100 m (Bräuer et al. 2011). At the intersection point of the MLFZ with the PPFZ, located in the northern part of the Cheb Basin, the Nový Kostel Focal Zone is situated, which is known for earthquake swarms characterized by the periodical occurrence of microearthquakes with magnitudes smaller than M_L = 5. Considering the entire seismically active NW Bohemia/Vogtland region, about 80 per cent of the seismic energy has been released beneath the Nový Kostel area during the last decades (Fischer & Michálek 2008; Fischer et al. 2014). About 10 km to the south, several degassing centres of highly concentrated CO₂, so-called mofettes, appear (Bräuer et al. 2011), for example the Hartoušov Mofette Field, which is the test site of our study (Fig. 2).

Figure 1.

The NW Bohemia region together with the main tectonical features: the Eger Rift, the Cheb Basin, the Mariánské Lázne Fault Zone (MLFZ) and the Počatky-Plesná Fault Zone (PPFZ). The Nový Kostel Focal Zone is indicated with a white dot, the Hartoušov Mofette Field is represented by a red dot. A digital elevation model (DEM) is underlaid (WGS 84, UTM Zone 33N).

Figure 2.

Temporary array configurations on the Hartoušov Mofette Field. Array1 was installed on 19 May 2015 (25 stations marked as red-bordered white triangles) on the North Hartoušov Mofette Field (NHMF), Array2 on 27 June 2016 (26 stations marked as blue-bordered white triangles) on the South Hartoušov Mofette Field (SHMF). The black-bordered transparent triangles represent Array3, which was installed on 22 November 2016 (95 stations) covering the entire Hartoušov Mofette Field (HMF). The red dots display the locations of mofettes. The white upscaled stations (A, B, C) are part of Array3 and they were selected for an amplitude comparision and a spectral analysis (Fig. 3). Station B was used as the reference station for cross-correlation computations (Fig. 5). The aerial photo is georeferenced (WGS 84, UTM Zone 33N).

The Hartoušov Mofette Field is a key site to study fluid flow and accompanying phenomena as it is characterized by continuous surface degassing of CO₂ forming dry and wet (filled with ground or precipitation water) mofettes. The total CO₂ discharge is estimated to 23–97 tons per day (Nickschick et al. 2015). Fischer et al. (2017) reported a fast increase of the CO₂ flow rate, monitored at a borehole on the Hartoušov Mofette Field after the earthquake sequence during Mai 2014. Simultaneously, an increase in the concentration of CO₂ bubbles in the water column was observed and hence, a clear link between fluid flow, accompanying stress changes and the occurence of earthquake swarms can be drawn. Therefore, it is presumable that fluid pathways (gaseous or water-filled) exist within the upper crust, that feed the mofettes at the surface and allow efficient and permanent fluid transport betweent the Nový Kostel Focal Zone and the degassing area around Hartoušov. According to geochemical studies, the isotope signature reveals an upper-mantle origin of these fluids with CO₂ being the major gas component at the surface (Bräuer et al. 2008). Fischer et al. (2017) claims the fluids to be a two-phased mixture of carbon dioxide and water (dominated by CO₂) in depth. Bussert et al. (2017) reported gas eruptions and blow outs during an ICDP drilling in a depth of about 80 m, indicating a highly overpressured fluidsystem beneath Hartoušov.

CO₂ degassing areas were succesfully located in the near-surface environment (down to 30 m) of the Hartoušov Mofette Field using matched field processing in the course of an ICDP pre-drilling site investigation (Flores-Estrella et al. 2016; Bussert et al. 2017). On the other hand, geophysical investigations on a larger scale reveal an enhanced fluid content within a depth of a few kilometres within the upper crust, such as a seismic tomography by Mousavi et al. (2015) or regional magnetotelluric experiments presented by Muñoz et al. (2018). So far, a connecting image about the fluid migration from depth (Mousavi et al. 2015; Muñoz et al. 2018) up to the visible mofette structures at the surface (red dots in Fig. 2) is missing.

We intend to contribute to the complex process of understanding the fluid flow activity within the upper crust in the Hartoušov area. With passive array seismology and a 3-D matched field processing (MFP) approach, such as Flores-Estrella et al. (2016) and Bussert et al. (2017), we aim at imaging seismic noise sources generated by fluid flow down to several hundreds of metres. We applied MFP to data from three dense temporary arrays installed on the Hartoušov Mofette Field during several hours in 2015/2016. As Nickschick et al. (2015) observed the surface gas flux to be temporally and spatially variable, we intend to study potential variations of the fluid activity in depth and on a short time scale (minutes to hours) as well.

MFP is a novel tool in the field of environmental seismology that is used on the exploration scale (several tenths to hundreds of metres depth) to locate seismic noise sources such as englacial meltwater channels (Walter et al. 2015), hydrothermal vents of geysers (Legaz et al. 2009; Vandemeulebrouck et al. 2010; Cros et al. 2011), a hydrocarbon reservoir (Corciulo et al. 2012) or CO₂ degassing zones (Flores-Estrella et al. 2016; Bussert et al. 2017). Since MFP was originally derived from ocean sciences, it involves the assumption of acoustic waves propagating from multiple point sources from different locations with different source strengths. A set of synthetic Green’s functions (MFP replica vectors) hereby represents potential acoustic sources that will be correlated with the measured array data set on a grid search basis. The highest beampower between the measured and the synthetic wavefield reveals the actual source location. It is essential to evaluate the MFP resolution if we consider seismic noise data from the heterogeneous earth characterized by full elastic wave propagation together with our array configurations. We therefore used finite difference modelling and subsequent MFP inversion which we compare to the MFP outputs derived from the actual field data.

2 DATA

We performed three measurement campaigns using different dense array set-ups that varied in aperture and amount of stations in order to image fluid channels beneath the Hartoušov Mofette Field on different scales and during different times (Table 1). Fig. 2 shows an aerial photo of the test site and the three temporary arrays (Array1–Array3). Array1 and Array2 were small-scale dense configurations (<120 m aperture, 25–26 stations), that covered the main mofette areas to the north and south of the test site respectively. The target depth was z = 100 m, aspiring the location of fluid conduit in the near-surface environment. A larger array configuration was installed later (Array3), which covered an area of 600 m aperture and comprised 95 sensors. The stations were randomly distributed on grassland and meadows and they covered the areas of Array1 and Array2, as we aimed to trace previously detected sources related to fluid flow to greater depth (down to 2000 m). The geometry of the arrays was determined by the exclusion of agricultural areas, streets and swampy parts as well as the forest.

Table 1.

Measurement campaigns on the Hartoušov Mofette Field in 2015/2016 (Fig. 2, Array1–Array3).

	Extent [m2]	Sensors	Period [hrs]	fs [Hz]
Array1 (NHMF)	120 x 100	25	1	100
Array2 (SHMF)	70 x 70	26	1	100
Array3 (HMF)	400 x 600	95	5	200

Array1 was installed on 19 May 2015 on the North Hartoušov Mofette Field (NHMF, Fig. 2) and comprised 25 stations. The array extent measured ≈ 120 x 100 m² and covered wet and dry mofettes at the surface. The second array was installed on 27 June 2016 on the South Hartoušov Mofette Field (SHMF, Array2, Fig. 2). It consisted of 26 stations and measured ≈ 70 x 70 m² covering dry mofettes. Both set-ups were equipped with Reftek Texan recorders and 4.5 Hz vertical geophones. Continuous seismic data were recorded during 1 hr sampled at 100 Hz during the night to reduce anthropogenic noise. On 22 November 2016, the third temporary array was installed, which covered the entire Hartoušov Mofette Field (HMF, Array3, Fig. 2). It comprised 95 stations (50 DSS DATA-CUBE³ with 3-component 4.5 Hz geophones as well as 50 Summit XONE channels with 3-component 10 Hz geophones) and measured ≈ 400 x 600 m². Continuous seismic data were recorded over 5 hrs during the night and sampled at 200 Hz.

In order to determine the seismic velocity of the subsurface, we performed active seismic experiments using several 200-m-long line profiles and a hammer blow source on the test site.

3 METHOD

3.1 Matched field processing

Matched field processing (MFP) is an array-processing method that aims to locate seismic noise sources (Porter 1994; Kuperman & Turek 1997). It is an extension of plane wave beamforming, where the plane wave assumption is violated as dense arrays are used to locate coherent noise sources within, beneath or in the close vicinity of the array. The approach is based on a grid search based correlation of the measured wavefield with a synthetic wavefield, generated by forward modelling (Kuperman & Turek 1997). It aims to estimate phase matches between both wave fields over a certain frequency bandwidth, considering each grid node as a potential source position and assuming the spatial coherence of the wavefield. For each grid point, the beampower between the measured and the synthetic wavefield is estimated. The gridpoint with maximum beampower value represents the most likely source location (Vandemeulebrouck et al. 2010; Walter et al. 2015).

The approach was initially developed in ocean acoustic sciences to locate and study the behaviour of sea mammals or to track submarines (Baggeroer et al. 1993; Thode et al. 2000). Clay (1966) emphasized the close connection between modal propagation of pressure waves and array processing in underwater acoustics to examine the signal-to-noise ratio of hydrophone arrays in a noisy oceanic environment. Based on his findings, Hinrich (1973) introduced the MFP concept using a vertical line array (VLA) of sensors and a linear MFP processor to correlate the measured wavefield with the modelled pressure wavefield. The VLA ideally spans over a significant portion of the water column to sample the wavefield spatially. The measured waveguide contains range and depth information of a source. Waveguide diversity hereby increases the ability to locate a source properly, but also requires a higher model complexity, as the transfer function from a particular source to each of the sensors is unique and different. Bucker (1976) refined the MFP approach by introducing a quadratic processor on the basis of a realistic model of the environment of the signal. He thereby located sound sources in shallow water. An overview of matched field methods in ocean acoustics is presented by Baggeroer et al. (1993).

Geophysical application of matched field processing started during the last decade on the exploration scale using dense seismic monitoring arrays on the heterogeneous Earth. Legaz et al. (2009), Vandemeulebrouck et al. (2010) and Cros et al. (2011) successfully located hydrothermal sources in shallow depth in the Waimangu geothermal valley (New Zealand) and at the Old Faithful Geyser (USA). Gresse et al. (2018) recently presented a multiphysics approach, combining MFP and electrical resistivity tomography, to image fumarolic conduits at Campi Flegrei (Italy). It has already been shown by Corciulo et al. (2012), that it is possible to distinguish between dominant and weaker noise sources with MFP on a hydrocarbon field and to image reservoirs down to 480 m depth. In the field of glacier seismology, MFP also serves as a novel tool. The study of Walter et al. (2015) shows the location of an englacial water tremor at the surface of the ice on a glacier site in Switzerland. Matched field processing also found recent application within our study area, the Cheb Basin (CZ). Flores-Estrella et al. (2016) and Bussert et al. (2017) presented the detection of mofettes (CO₂ degassing zones) close to the Earth’s surface, down to several tens of metres.

Although matched field processing is a relatively new tool in geophysics, in theory it can be applied to the location of any seismic noise source, which is coherent in space and time during a certain time window. The coherency assumption holds by using dense arrays and the simultaneous measurement of continuous seismic noise records at all stations. The recorded time domain signals d_i(t) of the array are transformed into the frequency domain using the discrete fast fourier transform. The cross-correlation of all sensor pairs generates the cross-spectral density matrix (CSDM) K_ij(ω) with i, j ∈ [1, 2, 3..., P] being the array elements using a single temporal snapshot:

$$\begin{eqnarray*} K_{ij}(\omega ) = d_{i}(\omega ) d_{j}^*(\omega ) . \end{eqnarray*}$$

(1)

d_i(ω) is the complex fourier spectrum of the ith sensor and |$d_{j}^*(\omega )$| the complex conjugated spectrum of the jth sensor. The CSDM is the frequency-domain equivalent of the time-domain cross-correlation of the recorded data.

The synthetic wave field is computed on the basis of a 3-D search grid, where each grid node acts as a candidate source position. In the simplest case, the Green’s function at array element i at distance a from a candidate source position is expressed as (Vandemeulebrouck et al. 2010):

$$\begin{eqnarray*} \hat{d}_{i}(\omega ,a) = \frac{1}{4 \pi a} exp^{\frac{-i \omega a}{c}} \end{eqnarray*}$$

(2)

favouring body wave propagation, or as (Corciulo et al. 2012; Walter et al. 2015):

$$\begin{eqnarray*} \hat{d}_{i}(\omega ,a) = \sqrt{\frac{2}{\pi a}} exp^{\frac{-i \pi }{4}} exp^{\frac{-i \omega a}{c}} \end{eqnarray*}$$

(3)

favouring surface wave propagation. c is the phase velocity of the source signal, which can be constant for the entire model space or frequency-dependent in case of surface wave dispersion. If a 2-D or 3-D velocity model is available, eqs (2) and (3) could be replaced by more complex Green’s functions.

The replica vectors in eqs (2) and (3) contain both phase and amplitude information. Phase information is usually regarded to be more stable than amplitudes, as random amplitude fluctuations occur due to, for example unknown propagation or sensor coupling issues. Therefore, we define our replica as:

$$\begin{eqnarray*} \hat{d}_{i}(\omega ,a) = exp^{\frac{-i \omega a}{c}} \end{eqnarray*}$$

(4)

and only include phase information in our synthetic wavefield. Accordingly, we perform amplitude normalization to our data traces.

The CSDM [K_ij(ω)] and the synthetic replica field (represented by |$\hat{d}_{i}(\omega ,a)$|⁠) are correlated at every candidate source position within the model space as well as summed over the frequency band of interest to estimate the beampower using the linear Bartlett processor:

$$\begin{eqnarray*} B_{Bart}(a) = \sum _{\omega } \vert \hat{d}_{i}^*(\omega ,a) K_{ij}(\omega ) \hat{d}_{j}(\omega ,a) \vert . \end{eqnarray*}$$

(5)

Alternatively, Capon (1969) proposed a high-resolution method, the minimum variance distorsionless response (MVDR), for the beampower estimation (using the inverse of the CSDM [K_ij(ω)⁻¹)]. It is an adaptive and non-linear MFP approach which acts as a maximum-likelihood filter between the measured and the synthetic data. MVDR is most efficient for low frequency data, where high SNR and spatial coherency can be assumed, and it is very sensitive to the accuracy of the synthetic wavefield. For higher frequencies, MVDR looses performance because the true heterogeneous Earth can not be represented by a homogeneous replica any longer. In that case, the linear Bartlett approach is recommended (Capon 1969; Tolstoy 2000; Cros et al. 2011; Corciulo et al. 2012).

3.2 Processing

The MFP inversion of the measured noise records was performed considering only the vertical component of each station. Since Array3 consisted of different instruments (Data Cubes and Summit XONE channels), the gains and units of the signal amplitudes were additionally adjusted.

All records reveal high spectral amplitudes between 5 and 25 Hz across the arrays (exemplary records of Array3 and a spectrogram are displayed in Fig. 3). In order to test the sensitivity of MFP within this frequency band, we performed both, calculations for individual frequencies and over a frequency band (Fig. 4). In Fig. 4 a it can be seen, that the beampower, as well as its range across the model space, increases for frequencies above 10 Hz and decreases again above 20 Hz. Figs 4(b) and (c) illustrate the focussing of the derived MFP locations. A cluster of locations can be observed in the southern part of the array, in a depth of z = 500–1000 m beneath active mofettes. That cluster is predominantly composed of locations derived from frequencies between 10 and 20 Hz. All other locations are widely distributed across the model space and mostly refer to sources at the surface or in shallow depth. MFP calculated over the frequency band of 10–20 Hz reveals the highest beampower (Fig. 4a) which is situated within the previously derived cluster composed of locations from individual frequencies (Figs 4b and c). Consequently, all records were bandpass-filtered between 10 and 20 Hz using a 2nd-order butterworth filter and the linear trend was removed from each seismogram.

Figure 3.

(a) Amplitude comparison of the seismic noise signal recorded during 5 min on three exemplary stations of the array (station A–C, positions indicated in Fig. 2). (b) Spectral amplitude of the seismic noise signal of station B.

Figure 4.

(a) MFP output (beampower maxima) derived from 5 min of noise signal from station B (Fig. 3) as a function of 41 frequencies between 5 and 25 Hz (red triangles). The black triangle stands for the beampower maximum derived from MFP performed over the frequency band of 10–20 Hz. The range (lines reaching the triangles) describes the resolution of the MFP beampower across the model space in vertical orientation (minimum - maximum beampower). (b) The dots represent the epicentral locations derived from the 41 individual frequencies, colour-coded according to the source depth. The triangles display Array3, the red dots display the mofettes and the black cross stands for the location derived from MFP performed over the frequency band of 10–20 Hz in a depth of z = 500 m. c The same epicentral locations, but colour-coded as a function of frequency.

The spatial coherence of the seismic noise wavefield across the array and within the chosen frenquency band is essential for the MFP approach. We therefore computed cross-correlations of 10 min of seismic signal, which are shown in Fig. 5. As reference station we used station B, located in the central part of the field. We observe high correlation energy across the array during several tens of seconds, which reveals the coherency of the wavefield. However, no clear phase can be traced. We interpret the rather long window of high correlation energy to be characteristic for the continuous degassing process.

Figure 5.

Cross-correlation of 10 min of seismic signal with station B (sensor 45) as reference station (Fig. 2). The data were bandpass-filtered between 10 and 20 Hz using a 2nd order butterworth filter.

The MFP inversions were computed on the basis of a homogeneous velocity model with the Rayleigh wave velocity c = 443 |$\frac{\textrm {m}}{\textrm {s}}$| for Array1 and Array2 (z = 100 m), which we derived from active seismic experiments, and the P-wave velocity V_p = 3460 |$\frac{\textrm {m}}{\textrm {s}}$| for Array3 (z = 2000 m) derived from the ambient noise tomography fof Mousavi et al. (2015). MFP was performed in the frequency range of f = 10–20 Hz and sampled with Δf = 0.5 Hz. The analysis followed conventional matched field processing using the phase-match replica (eq. (4), see Section 3.3) and the linear Bartlett processor (eq. 5) since it is more robust and recommended for higher frequencies. We computed the beampower for 10-min-long time windows and stacked all outputs for the length of each measurement period. The source grid spacing was set to Δx, Δy, Δz = 1 m, hence, the model space measured [x_min = 0 m, x_max = 700 m, y_min = 0 m, y_max = 700 m, z_min = 0 m, z_max = 2000 m]. The MFP outputs are displayed as relative power plots in 2-D and 3-D (multilayer). A logarithmic decibel scale was used [10*log₁₀(B(a))].

By choosing the phase-match replica, we did not take any amplitude information into account. It is therefore convenient to perform an amplitude normalization with respect to the maximum value of the signal of each station and each time window, as phase information are kept.

3.3 Synthetic test

Matched field processing is a beamforming approach based on the assumption of wave propagation in an acoustic environment (compressional waves only). Applying it to seismic data records from the Earth (full elastic wave equation) involves major approximations. In order to analyse the resolution of matched field processing in an elastic environment and considering the array geometry on the Hartoušov Mofette Field (Array3, Fig. 2), a synthetic array data set was modelled using the Finite Difference Code FD3D (Zehner et al. 2016). Our model mainly addresses geometrical aspects, such as the receiver locations with respect to the potential source distance. Wave propagation was simulated using an explosive point source in a homogeneous medium. The model space comprised [x = 700, y = 700, z = 2000] grid points with 1 m grid spacing for each plane. The receivers were placed at the free surface (analogous to the geometry of Array3, Fig. 2). All other model boundaries were defined as PML (Perfectly Matched Layer) boundaries with 20 grid points boundary width to suppress wave reflections. In accordance with the MFP results on the Hartoušov Mofette Field (see Section 4), the synthetic explosive source was placed at the coordinates [x = 250 m, y = 350 m, z = 1500 m] and emitted a synthetic coloured noise signal (f = 10–20 Hz), which was sampled at f_s = 10 000 Hz. The time of wave propagation was set to t = 2 s. We used a homogeneous velocity model with a P-wave velocity of V_p = 3600 |$\frac{\textrm {m}}{\textrm {s}}$| and a S-wave velocity of V_s = 2080 |$\frac{\textrm {m}}{\textrm {s}}$| (values taken from Mousavi et al. (2015)). The vertical components of the synthetic seismograms were inverted analogous to the field data using phase-matching MFP (eq. 4) with the linear Bartlett processor (eq. 5) since our data reveals a high frequency content (f = 10–20 Hz).

On the basis of the modelled data set, the suitability of different synthetic replicas was evaluated (Fig. 6). Additionally to the conventional phase-match replica (eq. 4), we tested the influence of amplitude terms on the location accuracy. eq. (2) includes an amplitude term, that favours body waves and should be suitable for this test, as we aim to locate a source at depth and assume the synthetic signal to be dominated by body waves. However, inverting for the location of the synthetic depth source using this replica and the P-wave velocity of V_p = 3600 |$\frac{\textrm {m}}{\textrm {s}}$| leads to a shift of the maximum towards the surface and to an offset in the horizontal plane (MFP max coordinates: [x = 292 m, y = 106 m,z = 0 m]). The beampower range is ≈34.42 dB (Fig. 6, blue dashed line). eq. (3) shows a replica with an amplitude term that favours surface waves. We used the Rayleigh wave phase velocity of c ≈ 0.92 * V_s ≈ 1900 |$\frac{\textrm {m}}{\textrm {s}}$| for another inversion test. Although a proper location of a source in depth can not be expected, a maximum is derived, which sticks to the surface of the model space (MFP max coordinates: [x = 310 m, y = 80 m, z = 0 m]). The beampower range is ≈10.64 dB (Fig. 6, blue dotted line) and smaller compared to the usage of eq. (2).

Figure 6.

MFP inversion results of a synthetic explosive source at [x = 250 m, y = 350 m, z = 1500 m] using eq. (4), eq. (2) and eq. (3). The most approximate source location was retrieved using the phase-match replica (eq. (4) with the MFP maximum at [x = 259 m, y = 351 m, z = 1460 m].

Conventional phase-matching (eq. (4)) leads to the most accurate source location (Fig. 6, red line). The MFP output is plotted as a multilayer 3D power distribution map that images the source location (Fig. 7a) representative for five depth slices between the surface and 2000 m depth. The values are given in decibel [10*log₁₀(B(a))]. The highest MFP power value of B = –27.57 dB was located at [x = 259 m, y = 351 m, z = 1460 m]. The power range across the model space shows a variability of only –0.45 dB (2.23 per cent) which is marginal, considering that a proper resolution is described by –3 dB , which is an energy decrease of 50 per cent (⁠|$\frac{\textrm {B}_{\textrm {max}}}{\textrm {2}}$|⁠) and desired to be spatially focused. Considering the medium P-wave velocity of c = V_p = 3600 |$\frac{\textrm {m}}{\textrm {s}}$|⁠, the frequency band of interest (f = 10–20 Hz) and λ = |$\frac{\textrm {c}}{\textrm {f}}$|⁠, the minimum and maximum wavelengths of the measured signal are λ_min = 180 m and λ_max = 360 m. With an array extension of ≈ 400 x 600 m² at most two wavelengths can be sampled properly, which can lead to a MFP resolution loss. Nevertheless, the test shows, that it is generally possible to locate sources in a distance comparable to the array extent with a reasonable accuracy. With an offset of [Δx = 9 m, Δy = 1 m, Δz = 40 m], the MFP maximum was located in the close vicinity of the synthetic explosive source (Fig. 7b). The highest mismatch arises in the vertical direction.

Figure 7.

MFP inversion in the frequency range off = 10–20 Hz of a synthetic explosive source ([x = 250 m, y = 350 m, z = 1500 m]) modelled in an elastic homogeneous environment and using the receiver geometry of Array3, as it is displayed in Fig. 2. (a) 3-D MFP output. (b) 2-D MFP output of the source plane at z = 1460 m.

Due to the reasonable performance of conventional phase-matching on the synthetics, we consistently used the MFP replica as stated in eq. (4) for the data inversion on the Hartoušov Mofette Field.

4 RESULTS

4.1 North Hartoušov Mofette Field (Array1, May 2015)

Fig. 8(a) displays the MFP output derived from 1 hr of continuous seismic record of Array1 (NHMF, Fig. 2) on the North Hartoušov Mofette Field. Five representative depth slices between the surface and 100 m depth are displayed. The values are given as MFP beampower values in decibel (dB). The MFP maxima of each depth level are listed in Table 2. At the surface, the MFP power was B = 10.02 dB. Down to 20 m, the power decreases slightly by 0.17 dB. From 20 m on, it constantly increases again towards the depth and reaches its global MFP maximum of B_max = 11.18 dB at z = 100 m. The geographical position of the respective depth slice maxima migrates from the southern part of the array at the surface towards the northern/central part of the array in 100 m depth.

Figure 8.

(a) MFP output from the North Hartoušov Mofette Field on 19 May 2015. (b) MFP output from the South Hartoušov Mofette Field on June 27, 2016. Both data sets were bandpass-filtered between 10 and 20 Hz and processed following conventional MFP phase-matching using a homogeneous velocity model with the Rayleigh wave phase velocity c = 443 |$\frac{\textrm {m}}{\textrm {s}}$| and a grid spacing of Δx, Δy, Δz = 1 m. The respective source planes are plotted below (NHMF: z = 100 m, SHMF: z = 25 m).

Table 2.

MFP power values at certain depths according to the MFP output of Array1 (NHMF) and Array2 (SHMF) installed on the Hartoušov Mofette Field (Fig. 2). The bold values refer to the beampower maxima.

z [m]	Bmax [dB] NHMF	Bmax [dB] SHMF
0	10.02	5.21
25	9.85	5.41
50	10.92	5.07
75	11.16	4.74
100	11.18	4.52

4.2 South Hartoušov Mofette field (Array2, June 2016)

The MFP output derived from Array2 (1 hr continuous seismic record), that was installed on the South Hartoušov Mofette Field (SHMF, Fig. 2), is displayed in Fig. 8(b). The respective power values are listed in Table 2. The MFP maximum of B_max = 5.41 dB was located beneath the central part of the array, at 25 m depth. Towards the surface, the beampower slightly decreases by 0.2 dB. Towards the depth, it decreases constantly down to 4.52 dB at z = 100 m. Compared to the MFP output of Array1, the beampower range of 0.89 dB is smaller and the source is poorly resolved considering |$\frac{\textrm {B}_{\textrm {max}}}{2}$| = –3 dB, which measures an area of ≈55 x 50 m² in the horizontal plane for z = 25 m and expands over the entire model space in the vertical plane.

4.3 Hartoušov Mofette field (Array3, November 2016)

Array3 (HMF, Fig. 2) covered the entire Hartoušov Mofette Field, including Array1 and Array2. The mean beampower during the measurement period of 5 hr is displayed in Fig. 9(a). The beampower ranges from B = –170.73 dB to B = –170.49 dB between the surface and 2000 m depth. From the surface, the MFP power increases constantly down to 1000 m, focuses between 1000 and 1200 m and then decreases slowly down to 2000 m. Hence, the maximum power is located in the depth of 1000–1200 m (Table 3). Following the maximum of each depth slice of the MFP output, a migration from geographic north towards the south can be observed with increasing depth. With 0.24 dB the beampower range across the model space is marginal.

Figure 9.

(a) MFP output of Array3 installed on the Hartoušov Mofette Field (HMF, Fig. 2) on the 22.11.2016. The data were bandpass-filtered between 10 and 20 Hz and processed according to conventional MFP phase-matching using a homogeneous velocity model with a constant P-wave velocity V_p = 3460 |$\frac{\textrm {m}}{\textrm {s}}$| and a grid spacing of Δx, Δy, Δz = 1 m. (b) MFP output of the same data set, but a previous amplitude normalization was performed for each station and each time window.

Table 3.

MFP power values at certain depths according to the MFP output of Array3 installed on the Hartoušov Mofette Field (HMF, Fig. 2), with and without amplitude normalization. The bold values refer to the beampower maxima.

z [m]	Bmax [dB]	|${B}_{\rm {\small {NORM}}}$| [dB]
0	−170.73	−67.77
200	−170.58	−67.10
400	−170.64	−67.14
600	−170.64	−67.18
800	−170.56	−67.19
1000	−170.49	−67.20
1200	−170.49	−67.21
1400	−170.51	−67.21
1600	−170.52	−67.22
1800	−170.52	−67.22
2000	−170.52	−67.23

We tested an amplitude normalization on the measured signals with respect to their individual maximum value, as we observed, that the amplitude levels differ across the array (Figs 3a and 2). Time records of 5 min lengths from three stations (A–C, Fig. 2) show a decreasing trend in the amplitude levels from north to south across Array3, which we intented to correct by normalization. The resulting MFP output shows a beampower range of 0.67 dB (Fig. 9b), which is enhanced by a factor of 2.8 compared to the MFP inversion without amplitude normalization. The beampower increases constantly from the surface (southern part of the array) down to 200 m depth, where the global MFP maximum is located, and decreases slowly down to 2000 m. A migration of the maxima towards geographic north is observable with increasing depth.

4.4 Spatio-temporal analysis

With matched field processing and continuous seismic noise records, we located sources down to 1200 m. Their epicentral locations coincide with active mofettes at the surface of the test site (Fig. 11). That supports the assumption that fluid pathways/reservoirs within the crust exist. Nickschick et al. (2015) repeatedly measured the gas flux at the surface of the Hartoušov Mofette Field and found that it changes rapidly over space and time (1–100 |$\frac{\textrm {g}}{\textrm {m}^{2}}$| per day). Heinicke et al. (2006) observed similar phenomena on mofettes in Caprese Michelangelo (Italy) and showed that the CO₂ gas flux is temporally variable depending on meteorological changes and the hydrogeological regime. These studies focus on the near surface regime, the first tens of centimetres of the top soil. Since we suppose that the mofettes in Hartoušov are fed by channel-like structures/from reservoirs within the upper crust, corresponding changes of the fluid flow activity in depth can be expected as well.

Instead of stacking all available time windows of 10 min length, we considered each output individually to monitor changes of the source activity in depth and over time. Fig. 10 a shows all single MFP locations (raw and normalized) during the measurement period. Hence, each dot represents the epicentral location of one 10-min-long time window and each time window is represented by two dots (raw and normalized processing (30 sources each, 60 sources in total)). It can be observed, that inverting the raw data predominantly leads to source locations beneath the geographic northern and central part of the array, whereas amplitude normalization results in source locations beneath the southern part of the array, which appear much less distributed in space as the majority clusters within an area of ≈50 x 50 m².

Figure 10.

Array configuration on the Hartoušov Mofette Field, together with the locations of the mofettes at the surface and the epicentral locations of matched field processing maxima derived from 10-min-long time records. (a) The green dots indicate MFP locations using the raw signal (30 sources), the yellow dots display maxima derived using signals with normalized amplitudes corresponding to the maximum value of each seismogram (30 sources). (b) The epicentres are colour-coded with respect to their source depth. Three main depth clusters can be defined (I–III). Exemplary MFP outputs for each cluster are presented respectively to the left.

Fig. 10(b) shows the epicentral locations of the same sources with a colour code according to their source depth. Three main clusters can be observed: one towards geographic north with sources around 1000 m depth (Cluster I, 11 sources), a central cluster with deep sources around 2000 m depth (Cluster II, 14 sources) and another cluster towards geographic south with sources around 300 m depth (Cluster III, 35 sources), which appears to be particularly well focused and is mainly composed of sources derived from normalized signals.

5 DISCUSSION

Noise source location on a mofette field is challenging, as the CO₂ transport and release appears in various forms as, for example bubble collapses at the surface in water-filled sinks, dry gaseous exhalations or fluid flow in subterranean conduits or channel-like structures with unknown floating behaviour at depth. It is most likely that the sources have an extended character and that their mechanisms are complex. However, the true physical properties of a mofette and its adjacent fluid conduits acting as a seismic source are rarely known. Additionally, multiple noise sources can be excited simultaneously. That in mind, our synthetic model can only roughly address the resolution of our MFP results, but provides a reasonable estimate in terms of geometrical parameters (array design versus source location).

We observed low beampower ranges of the MFP outputs for all arrays, meaning that the resolution of our located sources is rather low. Array1 and Array2 show MFP power ranges of 1.16 dB (maximum at 100 m depth) and 0.89 dB (maximum at 25 m depth), probably provoked by the array apertures, considering Rayleigh waves and the target depth of 100 m. The beampower ranges of Array3 (0.24 dB (raw) and 0.67 dB (normalized)) are low in both cases as well, and so is the MFP resolution of the located sources. This is again explainable by the proportion of the array aperture to the distance of the located source.

However, synthetic FD modelling demonstrated, that sources in a large distance compared to the array aperture can be located, but that the MFP resolution decreases accordingly. Therefore, we conclude that our source locations are reasonable. The results are furthermore supported by the congruent locations of active mofettes at the surface and the epicentres of the MFP maxima (Fig. 11) towards geographic north [Array1, Array3 (raw)] and towards geographic south [Array2, Array3 (normalized)] on the test site. Therefore, it can be assumed, that the outputs are linked to the degassing activity and fluid conduits in channel-like structures beneath the Hartoušov Mofette Field. That supports the assumption that the degassing process acts as an ambient seismic noise source and that MFP is a suitable tool to locate fluid flow in depth. Whether the vertical maxima distribution, characterized by rather low beampower ranges across the model spaces, is an artefact emerging from a reservoir source (point source character) or the actual image of a distinct channel/ fluid conduit (extended source), cannot be completely distinguished.

Figure 11.

MFP epicentres (white dots) detected on the NHFM and the SHMF (Fig. 8) as well as on the HMF [Fig. 9 (raw and normalized)]. The red dots indicate the positions of the mofettes. (a) and (b) show exemplary mofettes of the field (wet and dry). The triangles indicate stations (A, B, C), that were selected for an amplitude comparision and a spectral analysis (Fig. 3).

The cross-correlation of the recorded seismic noise data (Fig. 5) clearly shows, that the signal is coherent across the array and within our chosen frequency band. However, there is an ambiguity between the signal character and matched field processing, as it is aimed at matching phases between a data field and a synthetic field, but no clear phase onset can be observed in the data. By choosing a rather long time window of 10 min length for the MFP inversion, we accumulate energy from the most dominant source within that time window. In turn, the search for the most probable match (MFP beamformer) enables to locate the source.

The normalization of amplitudes for each station of Array3 results in the location of different MFP maxima compared to processing the raw data for the same time window. Both locations appear reasonable due to the proximity of their epicentral locations with active mofette structures at the surface (Fig. 11). Fig. 3(a) shows three exemplary stations of Array3 and highlights the differing signal amplitude levels across the test site. It can be observed, that the noise level is higher in the northern part of the array (station A) and relatively low in the southern part of the field (station C), whereas station B, situated in the central part of the field, shows intermediate values. The MFP maximum derived from processing the raw data is situated beneath a wet mofette (Fig. 11) as well as closely beneath station A. This mofette was active during the measurement period, which was observable due to gas release and intense bubble collapses at the surface of the water. It can be assumed that the fluid conduit feeding this mofette was partially filled with water as well, causing high amplitudes due to the enhanced fluid flow (mixture of water and CO₂). That leads to a dominance of stations close by within the MFP output and hence, to the location of the MFP maximum beneath that particular wet mofette. Sources with lower amplitudes therefore appear as relatively small maxima or can not be resolved, depending on the respective time window and the respective fluid flow and effective degassing activity. However, amplitude normalization predominantly leads to the location of MFP maxima beneath the southern part of Array3. At the surface, this area is characterized by dry mofettes. It can be reasoned that the water content within the conduit is lower compared to the norther part of the array, resulting in a lower amplitude level. From our observations we reason that amplitude normalization leads to the removal of station effects, equalizes the amplitude levels across the array and therefore the influence of each particular station on the MFP location. Due to the congruence of both MFP maxima with the locations of active mofettes on the test field, both procedures and results are considered to be reasonable.

The epicentral locations of our derived MFP maxima are displaced by a few metres to decimeters from the actual mofette positions at the surface. This can either be a methodical resolution issue or justified by the natural shape and geometry of a fluid channel/conduit. Since fluid flow within the subsurface is depending on various parameters, such as grain sizes, soil types, bioactivity, water saturation levels, the hydrogeological regimes or even local inhomogeneities as rocks for example, an ideal vertical channel beneath a mofette can not be expected.

Furthermore, we observed that the MFP locations derived from the small-scale configurations (Array1 and Array2) differ from those derived from Array3, even though partly the same mofettes were covered by stations. Mofette a (Fig. 11), for example, was covered by Array1 and Array3, but only located by Array1 and so on. One explanation could be the use of the Raleigh wave phase velocity for the MFP inversion of Array1 and Array2, but the P-wave velocity for Array3. We therefore favour different types of waves for the MFP location, which all contributed to the measured noise signal. Additionally, our campaigns took place between May 2015 and November 2016. It is very likely, that the effective fluid flow changed over time and space, which we already observed on a short time scale during a couple of hours. So, the MFP location derived from each array data set can only represent the fluid flow and/or degassing activity during the respective measurement period.

The spatio-temporal analysis reveals three distinct clusters beneath the test site, in a depth of 1000 m (Cluster I), 2000 m (Cluster II) and 300 m (Cluster III), from north to south, respectively. Due to the larger depth of their MFP maxima, Cluster I and Cluster II show a wider spatial distribution, which might be explained by the resolution loss and location uncertainties provoked by the ratio of the array aperture compared to the target depth. We assume, that these depth clusters might represent fractured areas where fluids from depth potentially accumulate over a certain time and pressure is released where the volume of the fluid grows. Altogether, the analysis reveals that the mofettes at the surface are fed by two individual fluid conduits/ channel-like structures, which seem to have a common origin in depth, situated beneath the central part of the test site.

6 CONCLUSION

The Hartoušov Mofette Field is characterized by multiple degassing areas at the surface (mofettes), appearing as small sinks, which can be dry, or filled with ground or precipitation water. These areas are often identifiable due to vegetation anomalies in the form of extinct grass or the existence of indicator plants for CO₂ degassing (e.g. bog cotton). The aim of this study was to image fluid channels/ reservoirs down to 2000 m depth beneath the Hartoušov Mofette Field and to fill the gap between the visible indicators for CO₂ exhalation at the surface and the indicators for fluid reservoirs/ channels in depth of 12–14 km beneath the test site (Mousavi et al. 2015; Muñoz et al. 2018).

With MFP phase-matching, fluid flow activity within the uppermost crust beneath the Hartoušov Mofette Field could be imaged in more detail. The outcomes of this study support the assumption that fluid channels/reservoirs exist, which could be directly linked to degassing centres (mofettes) at the surface. Altogether, synthetic tests, the stability of the located sources over time and the congruency of their epicentres with mofettes at the surface lead to the assumption, that two individual fluid conduits/ channel-like structures feed the mofettes at the surface and appear to have a common origin in depth, which is situated straight beneath the central part of the test site.

Our study shows that MFP renders the possibility to locate noise sources that are coherent across the array, but do not exhibit clear phase onsets. The ambiguity between the mofette signal character and matched field processing considering the search for phase congruences is challenging and needs to be investigated in further detail. Nevertheless, our results show that MFP is a chance to retrieve reliable source images even though picking of arrival times would be impossible.

ACKNOWLEDGEMENTS

We would like to express our sincere gratitude to the DMT group (Volker Schäpe, Jens Schweitzer) and the IGM GmbH (Steffen Uhlmann) for the supply of equipment and the help during field work. Further equipment was provided by the Geophysical Instrument Pool Potsdam (GIPP). Hortencia Flores-Estrella, Maik Henke and Helen Melaku contributed to this work by supplying results from active seismics and took part during measurement campaigns as well. We thank them very much. We also thank the Institute of Geography (Leipzig University) for providing the DEM of our study area. We would like to express our sincere gratitute to Olaf Hellwig (TU Bergakademie Freiberg) for providing his code FD3D and Katrin Hannemann for proofreading the manuscript and fruitful discussions. This study was funded by the Deutsche Forschungsgemeinschaft (DFG) in the framework of SPP-1006.