Long term measurements from the M\'atra Gravitational and Geophysical Laboratory

Summary of the long term data taking, related to one of the proposed next generation ground-based gravitational detector's location is presented here. Results of seismic and infrasound noise, electromagnetic attenuation and cosmic muon radiation measurements are reported in the underground Matra Gravitational and Geophysical Laboratory near Gy\"ongy\"osoroszi, Hungary. The collected seismic data of more than two years is evaluated from the point of view of the Einstein Telescope, a proposed third generation underground gravitational wave observatory. Applying our results for the site selection will significantly improve the signal to nose ratio of the multi-messenger astrophysics era, especially at the low frequency regime.


Introduction
Preparation for the next generation ground-based gravitational detectors requires careful pre-analysis at the proposed locations. As for previous gravitational wave detectors, identifying noise sources and constructional risks is a key task. In case of underground installation and improved low frequency operation this must be based on long-term measurements to identify noise types and its origins. Parallel and combined measurements applying the most novel techniques for environmental tests is also beneficiary for the design of the proposed facility and for high-accuracy measurements with it.
The main advantage of underground operation for gravitational observatories is the improved sensitivity at the low frequency regime, 1-10Hz. The observation of gravitational waves with earth-based detectors and the birth of multimessenger astronomy in the last years [1,2,3,4,5,6] revived interest in the scientific and technological challenges due to the underground operation [7,8]. The scientific value of the existing plans of building an all-in-one observation facility increased considerably. Its improved sensitivity, extended frequency range, together with capabilities of measuring polarization and direction offer an unprecedented discovery potential [9].
In order to explore the technological and scientific background of underground operation the Mátra Gravitational and Geophysical Laboratory (MGGL) was established by MTA ), along a horizontal tunnel of the mine, 1280 m from the entrance and 88 m below surface, however, the related research activities are extended to other areas of the mine and the surroundings in the last years. The initial instrumentation and the first results were reported in Refs. [10].
In this paper we report the results of continued research summarizing the geophysical characteristics of the Mátra Mountains as an Einstein Telescope site candidate. In this respect the analysis of the collected seismic noise data of more than two years are probably the most informative. We report them together with geophysical data, rock characteristics, infrasound, electromagnetic and muon flux measurements. The paper is organized as follows. First we introduce the geological and seismological environment of Mátra Mountains. A short report of the mechanical properties of the grey andesite rock of the Mátra is part of the geophysical survey. Then the seismological data collected in the last two years are shown and analysed. Here both the data from broadband seismometers of the Wigner RCP and also the custom made seismometers of the Warsaw University are shown. In the third section underground infrasound measurements are reported. Then the damping of electromagnetic waves in the andesitic rocks of Mátra is estimated using underground and surface electromagnetic measurements. Finally the evaluation of the collected muon flux data is shown, demonstrating the tomographic capabilities of the installed detector technology with rock density maps.
In the following sections we shortly introduce the instruments and elaborate the long term data, whenever it is possible. This includes all seismic information that was collected in MGGL. In the following the complete data collection period will be referred as Run-1 when compared to the previous, shorter, Run-0 data of our previous report [10].

The geophysical environment of MGGL
The lithological composition of Mátra mountain range is moderately blocky andesite, as it is shown in Fig. 1, where the green areas denote various andesite types formed at the same geological era [11]. In the following we characterise the seismic activity of the surrounding area and also the mechanical properties of the typical hard rock from Gyöngyösoroszi mine. In order to enhance comparison with other sites we also show quantitative measures by calculating the related ground displacement and seismic hazard of the Mátra Mountains.
2.1. Seismicity of the Mátra Mountains and the surrounding areas. In general earthquakes are more important for the stability of the underground facility and small ones are not critical regarding the observations. The seismicity of the Carpathian basin as a whole can be considered moderate. The known earthquakes of the area with magnitude larger than 3.0 are shown in Fig. 3. According to the historical collected data on average one earthquake with magnitude M ≥ 5.0 can be observed in the Carpatho-Pannonian region annualy. Also the seismic activity from neighbouring open pit mines can be identified by the ET1H seismological station deployed in the MGGL with the help of the stations of the Hungarian National Seismological Network [14]. The level of seismic activity in the Mátra Mountains is quite low. Fig. 2 shows the epicentres of the known M ≥ 3 earthquakes in the area (19.69-20.18E, 47.8-48.0N) based on the data of the MTA CSFK GGI Earthquake Catalogue [12] and Hungarian National Seismological Bulletins [13]. Only earthquakes with magnitude greater than 3.0 are shown as lower magnitude events can be misclassified quarry explosions. It can be seen that only three small earthquakes (M ≤ 3.5) were ever observed in the area which occurred in 1879, 1895 and 1980.
2.1.1. Ground displacements caused by seismic events. We have selected four events to characterise the ground displacement in the MGGL caused by different seismic events. For the computations of displacements 100 sps streams were used. The data before the instrument correction were detrended and filtered by a second order high pass Butterworth filter with a corner frequency of 1 Hz.
As the seismicity in the Mátra area is mainly determined by the mining activities in the quarries, we computed the ground displacement for characteristic explosions carried out in the three most active mines. Table 1 shows the maximum displacement amplitudes for the three components in the case of the selected events. The magnitude of the displacement is similar for the mines Gyöngyössolymos and Kisnána, while the maximum displacement for the explosion belonging to Kisnána mine is an order of magnitude smaller.
As the recent Tenk earthquake (2013-04-22, M=4.8, epicentral distance=42 km) occurred relatively close to the MGGL it may be of interest. Unfortunately, at the time of the earthquake the ET1H station was not yet installed. However we can use the data of Piszkéstető station, which is located on the surface, 4.88 km from MGGL  and was operational. The fourth row of Table 1. shows the displacements observed there. It can be seen that the Tenk earthquake produced around 20-30 times larger displacements than the strongest explosion. It must be noted that in the MGGLdue to the large subsurface depth -this value probably would have been smaller.  (Table 2.). These correspond to 50%, 39%, 10%, 5% and 2% exceedance in 50 years, respectively.  2.2.1. Measured elastic moduli and rheological parameters. As it is well known, the static elastic moduli of rocks measured in laboratory are different from the dynamic ones determined from wave propagation speeds [17,18,19,20,21,22]. The difference is usually attributed to various heterogeneities, like microcracks, porosity and grain structure. The Kluitenberg-Verhás body of thermodynamic rheology, which is derived from non-equilibrium thermodynamics with internal variables, provides a simple modelling possibility and explanation [23,10]. If the material relaxation times of geometric effects are in the order of the operational domain of the low frequency part of ET, then rheological properties can influence the reliable detection of gravitational waves. Therefore we have performed laboratory measurements in order to determine the rheological properties of the gray andesite of Mátra. The investigated samples were cutted from blocks originated in construction works at the vicinity of MGGL. This is a middle gray, small grained, isotropic piroxen type andesitic rock, which is a differentiated upper miocen (tortona type), about 14.5 million years old and considered typical in this region.
The measurements were performed in the rock mechanics laboratory of Kőmérő Ltd. with the help of a hydraulic instrument (maximal compression is 150kN), HBM C6A (1MN) load sensor and a HBM Spider 8 & CatmanEasy collected the data. The diameter and the lengh of the cylindrical sample was 37.99 mm, and 78.33 mm and its mass was 0.2213 kg. The uploading speed was 0.7 kN/s in every cycles. A hysteresis type measurement with increasing stress amplitudes was chosen to determine rheological parameters. In our calculations only the uploading parts of the last two cycles are considered.
The force and the axial deformation are shown on the left side of Fig. 4 as functions of time. The right side of Fig. 4 shows the axial stress as function of axial and lateral deformations at the positive and negative horizontal axes, respectively. The rheological hysteresis and also the apparent permanent deformation is visible at the end of creep periods. Therefore andesitic rocks from Mátra show clear deviation from ideal elasticity with properly designed laboratory experiments.

2.2.2.
Kluitenberg-Verhás body and static-dynamic elasticity. The isotropic Kluitenberg-Verhás body is given by the following relations between the stress σ and strain : Here the dot denotes time differentiation, the subscripts d and s refer to the deviatoric and spherical stresses and deformations as well as to the related material parameters, like the τ d and τ s deviatoric and spherical relaxation times. In our case a cylindrical laboratory sample was prepared according to ISRM (International Society of Rock Mechanics) standards, therefore the uniaxial stress, σ, and the axial and lateral deformations a and l determine the spherical and deviatoric stresses and strains as σ d = 2σ/3, σ s = σ/3, d = 2( a − l )/3 and s = ( a + 2 l )/3.
The elastic moduli E 0d = 2G and E 0s = 3K are the well known static Lamé coefficients and E 1d and E 1s are the deviatoric and spherical viscoelastic material coefficients. If the E 2d and E 2s dynamic parameters are neglectable, then the Kluitenberg-Verhás body simplifies to the Poynting-Thomson-Zener body [23], which is the standard model of creep and relaxation phenomena in rock mechanics (also with the names generalized Kelvin-Voigt and Hill-Maxwell body [24]). Then it is convenient to transform equations (1) and (2) into a hierarchical form: Here the parameters b d = E 1d /(τ d E 0d ) and b s = E 1s /(τ s E 0s ) measure the deviation from the ideal elastic Hook body. If b d = b s = 1, then the material is apparently completely elastic. The propagation speed of longitudinal and transversal waves determine the material parameters, E d,dyn = b d E 0d and E s,dyn = b s E 0s , respectively. These are called dynamic Lamé coefficients. Both for the static and dynamic cases the Young modulus and the Poisson coefficient are calculated from the Lamé coef-

Material parameters.
One can determine the static and dynamic elastic moduli both from the cyclic loading of the sample given above and the dynamic ones can be determined from direct laboratory measurements of the sound speeds. This direct laboratory measurement of the propagation speed of longitudinal and transversal waves gives the dynamic Young modulus (38.6 ± 1.1) GPa, and the dynamic Poisson coefficient 0.18 ± 0.01. Furthermore, the time dependent data of the chosen rock sample was analysed using the differential equations of the rheological models. In particular we have assumed, that the measured deformation values are given and we have determined the best parameters from the differential equations (1) and (2) of the Kluitenberg-Verhás body to obtain the stress, both for the deviatoric and spherical components. According to these calculations the coefficents of the second derivatives of the deformation, the E 2d and E 2s coefficients, can be neglected and the Poynting-Thomson-Zener model, (3) and (4), was applied. The obtained best fit parameters are shown in Table 3.  Table 3. The static and dynamic elastic coefficients according to rheological data evaluation of the uploading parts of the 3 th and 4 th cycles of the laboratory experiment. The dynamic Young modulus and dynamic Poisson ratio ν are consistent with the values obtained from propagation speeds of the P and S waves.  Table 4. The rheological parameter, b = E 1 E 0 τ , characterising the deviation from the ideal elastic regime both in the deviatoric and spherical cases.   According to our measurements the typical gray andesite of Mátra is not ideally elastic and the deviation from elasticity is of rheological origin. The experimental parameters can be obtained both from the wave propagation speed measurements and from the cyclic loading experiments. These coefficients are consistent as it is shown by our preliminary calculations. The obtained transition regime, which is proportional to the inverse relaxation times, is 0.02-0.1Hz, which is below the low frequency sensitivity of ET. However, our experimental methodology prevents the detection of faster relaxation modes, therefore further investigations are necessary.

Seismological measurements I.
The ET site selection preparation measurements [9,25,26] selected three best candidate sites for the Einstein Telescope in Europe. The selection is based on the low seismic noise level in the critical low frequency range 1-10Hz. The survey investigated several underground locations, each of them about for a week. In order to have a better estimate of the average seismic noise it is reasonable to collect at least two years of data and analyse the annual and seasonal changes in noise level originating from natural sources and human activity. The above mentioned survey was the main motivation in the establishment of MGGL. Due to the ongoing reclamation activity in the Gyöngyösoroszi mine the vicinity of the laboratory was also subject of maintenance works. Therefore in our case both internal and external human activities were contributed to the noise level.
Seismological data collection was performed by two Guralp CMG 3T low noise, broadband seismometers [10], and also by the custom made seismic sensor developed in the Warsaw University. The Guralp instruments are sensitive to ground vibrations with flat velocity response in the frequency range 0.008-50 Hz. The Guralp seismometers were calibrated and cross calibrated, operated according to protocol of the seismometers in the Hungarian National Seismological Network. The self noise of the seismometers was below the low noise model from 0.02 Hz to 10 Hz [27].
The custom made Warsaw seismometer uses one vertical and two horizontal geophone sensors mounted firmly in a single aluminium block placed inside a metal housing along with a data acquisition system. The geophones (LGT-2.5 and LGT-2.5H) used as the sensors have the lower corner frequency of 2.5 Hz. The data acquisition system sampled the analog signal with the frequency of 125 Hz and the 32-bit resolution. The Warsaw seismometer was calibrated by comparing the data with a Trillium seismometer, and also with a Guralp CMG 3T during a data calibration session at the MTA Wigner Research Centre for Physics in January 2015.
One of the Guralp instruments (ET1H) and the Warsaw seismometer (WARS) were permanently installed in the MGGL. The seismometers were deployed on separated concrete piers which were connected to the bedrock. Between the piers and the seismometers a granite plate has been placed. The other Guralp instrument (hereafter GU02) was used in a measurement campaign in the first two weeks of June 2017 in a measurement cabin, constructed next to the main tunnel and prepared for seismometer installation. The mutual performance of the ET1H and WARS seismometers are demonstrated on Fig. 6, where the averages of the day 2017-07-01 and 2017-07-02 are calculated for each instruments and also the ratio of the averages at given frequency range. The agreement is similar to other customary chosen days during the measurement campaign. Therefore the cross-calibration of the instruments is satisfactory.
The data acquisition period of the two instruments slightly differs due to operational problems, however the amount of the collected data is similar. In the following first we analyse the ET1H data, then the WARS data. The methods and the elaboration are different, in order to show more aspects of the noise measurements. -2016-10-16, for 25 days. Therefore 741days of data was collected for the horizontal directions and 716 days in the vertical direction. In our analysis we followed the data processing method of Ref. [26]. The particular problems observed in long term data analysis suggested some methodological improvements, and additional characteristics. These were reported in Ref. [28] where the precise definitions and the justification of the methodology of data processing can be found. In this section some basic definitions are recalled for clarification.
The velocity Power Spectral Density (PSD) is defined as [26,28] where f s is the sampling rate -in our case it is 100 Hz -, N is the length of the analysed data sample -N = 5000 -, and W = 1 . We did not use the advantage of fast Fourier algorithm on the expense of increasing the lowest frequency value. The coefficients , represent the Fourier transform of the deviation of raw velocity data v[n] from its average value v . In our analysis PSD-s were calculated with 50 s data samples with f s = 100 Hz sampling rate. Before further processing, raw data were highpass-filtered with f HP = 0.02 Hz and the overlap is 3/4 due to the window function. To characterize spectral properties the acceleration Amplitude Spectral Density (ASD), will be used and displacement rms will be applied as cumulative property. This is the square root of the integral of displacement PSD between two frequency values where l is the low cut-off index -in Beker's paper this is chosen to be 2 Hz -, the K is the high cut-off index -usually the Nyquist frequency -and T = N fs . In our case l is chosen to be 1 Hz or 2 Hz and K is either the Nyquist-frequency or 10 Hz, because we use rms 2Hz of Beker and also calculate the rms 2−10Hz and rms 1−10Hz values. With these new measures we can drop the irrelevant frequency interval above 10 Hz and consider the technologically already available 1 − 2 Hz region [28]. The commonly used comparative measure for the spectral properties of the sites is a particular value of the amplitude spectral density, the so called Black Forest line: It is worth to give the various rms values corresponding to the the reference spectral density represented by the Black Forest line: rms BF 2Hz ≈ 0.1 nm rms BF 2−10Hz = 0.1 nm and rms BF 1−10Hz = 0.29 nm. An other important aspects of the long term data evaluation are the use of percentiles and the choice of intermediate averaging periods. Previous studies applied the mode of the data for the representation of a characteristic mean value. However, the median, and also the other percentiles, are less sensitive for discretisation and the inevitable averaging. Moreover, naturally select the representative data without a necessity of filtering short large noise bursts. Therefore in the following the analysis of site properties is mostly based on the median of the data. However, we will demonstrate the most important differences in the spectral representation and also calculate the mode related rms 2Hz of Beker for a clear comparison with the previous studies.
It is also worth to mention, that the basic Fourier length used at the Fourier transformation was 50 s for Guralp instruments and 128s for WARS. For the ET1H and GUO2 Guralp instruments we have introduced a short time averaging (STA) period of 300 s and for the two-year data the use of daily averages was convenient. The daily periods were called intrinsic averaging (INA), because it considers the natural periodicity of the data. In other words in order to obtain the long term 10 th , 50 th and 90 th percentiles first the short therm averages (STA) were calculated and then particular chosen periods in each day (INA) enable the comparison of working hours, night time or whole day data. See [28] for more details about the chosen methodology.

3.2.
Long term seismic results. A particular factor in our noise analysis is the ongoing mine reclamation activity in the Gyöngyösoroszi mine. As a result, the investigation and identification of various noise sources originating from external and internal human activity, machine noise, construction works, train noise, etc. proved to be difficult. In November 2016 a three-shift operation period has been started with increased industrial noise, present also during the nights. In order to compare the noise types of these kind of activities in our analysis we have defined three periods for each day for our study: (a) the whole day, (b) night period (20:00 -2:00 UTC) and (c) working period (9:00 -15:00 UTC). In the following we present a comparative analysis of long term data considering seasonal changes, external and internal human noise and also depth dependence for a shorter two weeks long measurement campaign performed by our second Guralp instrument, GUO2, at -404 m depth.
3.2.1. Complete Run-1 results. The acceleration ASD-s for the two-year observation period are shown in Fig. 7. The borderlines of the blue colored area are the 10 th and 90 th percentiles of daily 300 s data. The dotted black lines are the modes of the 1800 s averages, according to the methodology in Ref. [25]. The medians and modes are closer to the Black Forest level in horizontal directions. It is also remarkable that the mode underestimates the median, that is most of the days are noisier than the mode. at a given frequency. This is well represented in the corresponding rms values given in Tab. 5, too. Here the first line rms 2Hz , calculated from the mode of the data. It is worth to recall that it was 0.12 nm in the previous short term measurements of Beker in the same place. As it can be seen, the rms calculated from the median is much larger, about 20% more. It is also worth to mention that the Black forest line reference value of rms 1−10Hz is 0.29 nm, and the rms 1−10Hz normalized to this value is less than the rms 2−10Hz , therefore the site is less noisy at the lower frequencies, as it can be seen in Figs. 7, as well.
The role of the human and industrial noises is shown in Fig. 8, where the spectral densities for the working and night periods is plotted in the North-South and the vertical directions. It is remarkable, that half of the frequency range is below the median of the data (blue line), the asymmetric relative position of the blue area at working period indicates the presence of short noisy periods between 9:00 -15:00 UTC. We have to mention that from the end of 2016 the reclamation works were performed close to the MGGL in a three-shift schedule.  [29].
In order to illustrate the cultural and industrial noises the frequency dependence of the ratio of the working and night periods are plotted in Fig. 8. The cultural noises start at about 0.7 Hz and reach their maxima between 2 − 3 Hz and 10 − 20 Hz. We may suppose that the night shifts cause less noises, the main works are done during the daytime. Then the observed noise level at the night periods can be considered as an upper limit for an operating underground GW detector facility, where the equipments are optimized for a low noise operation.
The corresponding rms values are for the night and working periods are given in Table (6. We can see, that the rms of the noisier working periods can be two times higher than in the calm night ones.
The differences at longer periods in Run-1 were analysed by three different methods. First, the acceleration ASD-s of the night periods were plotted for each year in Fig. 10. The 90 th percentiles at higher frequencies show the increasing industrial noise at night due to the three-shift schedule of the mine works. Interestingly this is not apparent in the mode and in the median of the data, because the particular activity seemingly put up only short noisy periods. In Table 7 the rms 2−10Hz values  Table 6. The rms values in nm for the night and working periods in Run-1, according to the panels of Fig. 8. of the night periods of each year are shown, corresponding to the middle panels of Fig. 10. There is no significant annual variation during Run-1 of MGGL.
The seasonal averages are plotted in Fig. 11 including the rms values in Table 8, which presents that the calmest seasons are the spring and the summer.   Table 8. Seasonal variation of the rms values of the seasons according to the data shown in Fig. 11. The first part is the rms 2Hz of the mode and the second one is the rms 1−10Hz of the median.
Finally, the seasonal changes are represented with a timeline plot of the the daily rms values in Fig. 12. As we can see, there is an annual trend in the curves, the late spring and the early summer are the most silent periods in the 1 − 10 Hz frequency range. Note, this seasonal difference is not apparent in the rms 2Hz timeline, which is at the bottom panel of Fig. 12. This seasonal variation of the lowest frequency noise is in agreement with the observable also in Fig. 11  interval, the percentiles were computed directly from the 300 s averages without INA. The modes were calculated from 1800 s averages as in the previous section. In Fig. 13 the comparison of the two-week and the total Run-1 ET1H data are plotted. One can see, this two-week period is representative, there is no significant differences of the plots, in spite of the INA in case of Run-1.
We have plotted together the ET1H and GU02 data in Fig. 14. There it is apparent that the main attenuation is in the intervall of 1-4 Hz for the horizontal and 1-7 Hz for the vertical direction. This interval is crucial for the low frequency noise budget of the proposed ET, furthermore, this frequency range dominate all the rms values. The corresponding rms values for GU02 are in Tab. 9. In Fig. 15 the night and working periods are shown for the GU02 station with the corresponding rms values in Tab. 10. The ratios of night and working periods are also shown in Fig. 16. 3.3. Summary. We analysed the long term Run-1 data of the ET1H station and compared the data from a representative two-week period of the −88 m deep ET1H station with the −404 m deep GU02 station. For higher than 2 Hz frequencies, there    are no significant annual changes, but for 1 − 2 Hz the spring-summer time shows annual minimum. Comparing the deeper and shallower noise data we have observed that the decrease of seismic noise spectral amplitudes in the 1-8 Hz frequency range is approximately 60%. This range is crucial for the low frequency performance of Einstein Telescope.
We emphasize, that almost in 90% of the observation period detected the noise level below the Black Forest line at night (see Fig. 15). The related average horizontal rms 2−10Hz = 0.0742 nm and the rms 1−10Hz = 0.213 nm. The first value is equal to the mode related rms 2Hz value, which is 0.0745 nm.

Study of seismic noise with the Warsaw seismometer
As it was described in Section 3, the custom made seismometer of the Warsaw University (WARS) is located near to the ET1H seismometer, therefore they should measure the same noise, although at low frequencies there can be differences due to scattering of seismic waves in the mine. The y-axis of the sensor is pointing at the direction of 25 degrees NW. 4.1. The collected and analysed data: The Run-1 data of the WARS seismometer statts from 24 May 2016 till 2 July 2018. In this period 654 days were used for analysis with the rest of the days being unavailable. Days not included in the study will be referred to as "missing-days" henceforth. The data for the study were recorded in each hour. Out of these 654 days, 30 days do not have complete 24 hours and will be referred to as "gap-days" henceforth. After eliminating the missing durationsfrom the gap-days, total of 15 290 hours of complete data were available for analysis.

4.2.
Analysis of the WARS data: The acceleration power spectral density (PSD) [30] was calculated by Welch's method [31]. The data for each hour were divided into 2048 length segments with each segment windowed by Hann window function [30,32,33] using half the length segment as overlap [30]. No zero padding was used. The acceleration amplitude spectral density (ASD) for each hour was calculated by taking square root of the acceleration PSD. The data for each hour were available for each of the three axis for all the 15290 hours so a total of 3×15290 ASDs were generated.

4.3.
Monthly analysis of the WARS data: The ASD generated for each hour were used to study the variation of ASDs with time over the period of each month i.e spectrograms were generated for every month for the period of the study. This was done for all three axis to study the variation in detail. We present an example    acceleration ASD and binned it on a plane spanned by the frequency and the logarithm of acceleration ASD. We then plotted it using the color shading for presenting how often each value occurs in the data. The log ASD axis was divided into twenty bins per decade. We show the results for the two axes in Fig. 18. The variation of the ASD spans about 0.7 dex at low frequencies below 5Hz and goes to almost 1.5 dex at about 20Hz.

4.5.
Complete daily analysis: In order to investigate several types of spectra during different working times of the day in more detail we generated the acceleration ASD for periods corresponding to the working shifts in the mine, 09:00-15:00 UTC and 20:00-02:00 UTC. The acceleration ASDs for these periods were then binned and plotted for each axis for the whole data set. We present these in the panels of Fig. 19 for the working, and in panels of Fig. 20 for the night periods, respectively, like in the case of ET1H station. These plots show less variation than the plots for the entire data set. The lack of characteristic bimodality at higher frequencies indicates lack of noisy periods at night. The working period shows additional lines, that are not visible at the night. The night spectra have very little variability less than 0.3 dex over the entire frequency band. 4.6. Site parameterization. The values of rms 2Hz were calculated for each hour of available data and the cumulative distribution function was plotted for the three axes of the detector in Fig. 21. The median rms 2Hz ≈ 0.18 nm in the vertical direction and ≈ 0.16 nm in the horizontal directions. The two horizontal rms 2Hz distributions are different below the median. It is unclear whether this is an instrumental effect or if this corresponds to the properties of the seismicity in MGGL. 4.7. Conclusion: In conclusion our measurements show that the seismic noise level is very low at the Mátra site. In spite of the spectral differences the median rms 2Hz values are in agreement with the ET1H data of Table 5, measured by a different device.

An electromagnetic study on the signal attenuation in the ELF range
Since the gravitational wave interferometer's mirrors of the proposed ET is going to be stabilized by means of magnetic fields, man-made and natural electromagnetic (EM) signals may also result in contaminated time windows of the gravitational wave observation. Magnetic noise from local sources can be identified based on correlation analysis utilising a network of extremely low frequency (ELF) field observatories. However, global EM components in the frequency range of interest may result in undesirable noise load [34]. The global thunderstorm activity exciting the Earth-ionosphere cavity continuously keeps a background signal at certain resonance frequencies called Schumann resonances [35]. Installation of the gravitational wave detector in a subsurface location may lead to lower EM noise even in the ELF range, 3 − 20Hz.
An EM investigation has been carried out in close vicinity of the MGGL which aims to provide an estimation of the attenuation of EM signal with the depth in the Mátra andesite rock. ELF range geomagnetic observation stations has been installed at the backfilled end of a 140m deep cave and in the close vicinity of the surface projection of the subsurface measurement site. For technical reasons, the subsurface and surface recordings have been run with no overlap in time. Therefore an indirect processing method has been applied utilising the recordings of the very low-noise ELF EM observation site of Hylaty station as reference observation. The amplitude rate between the subsurface and the surface station has been estimated based on the experimental transfer functions between the subsurface-Hylaty and the surface-Hylaty relation.
5.1. The data acquisition system and the observation. The study is based on two mobile observations and an observatory data: time series from a subsurface site, a surface station and from a reference station. The subsurface measurement location could have been accessed from the town Mátraszentimre by an elevator crossing levels down to a tunnel at 140m depth. The observation site has been hundreds of meters away from the vertical mine of the elevator and from all active electronic facilities at that level. Notice that active loads have been in operation during the measurement in other levels of the mine, around 50m and 20m distance vertically. The subsurface data acquisition system consists of a Lemi-423 wide band magnetotelluric station and two Lemi-120 induction coils. The coils have not been fitted to a common plane to increase the separation distance to reduce the inductive coupling effect between the coils. The orientation of the coils have been set to the orientation of the tunnel: NNW-SSE and ENE-WSW, with NNW 22 • to North. The timing of the subsurface station has been initialised at the surface by means of GPS synchronisation and had been transported in the proposed site in operating status with no sensor connected. During the whole measurement session (Period I), the inside clock of the station has not been synchronized due to the lack of GPS or NTP access. The station was powered by a fully charged 84Ah battery. The data logger has been installed 40m away from the sensors along the tunnel. The connector of the coils and the whole instrumentation have been waterproof sealed, since high humidity and leaking water could lead to shortcut of the connector part. The sampling frequency has been set to 2kHz.
The surface measurement has been carried out after the subsurface measurement (Period II) and by means of the same Lemi station which has been set up at the subsurface site before. The coils have been buried 50cm deep in a common horizontal plane in perpendicular orientation. The separation of the sensors was more than 10 meters and the Lemi-423 station has been buried 40m away. The Lemi instruments and the wires have been buried preventing harmful activity of forest animals. The orientation of the coils at the surface stations has been set to NNE-SSW and WNW-ESE adapting to the local environment conditions. The sampling frequency has been set to 2kHz at this mobile station, too.
The Hylaty geophysical station was proved to be the optimal ELF observation facility as reference station based on data quality and geographic aspects too, [36]. The Hylaty station provides very low noise ELF geomagnetic data at around 888Hz sampling rate which is optimal for the proposed study as reference.

Data processing.
5.2.1. Data preparation. The primary goal of the study is to provide an estimation of EM attenuation of the overlying andesite block in the lower ELF range, especially at the first Schumann resonance frequency. In the aspect of the measurement configuration it is identical to the relation of the absolute value of the empiric amplitude transfer functions of the subsurface and the surface stations at certain frequency. Since the surface and subsurface observation time periods were not overlapping the direct comparison of the amplitude spectra was not possible. The concept of the processing is based on utilization of a high-quality, low-noise ELF reference station data covering all the studied time intervals. The indirect method is basically as follows: determination of Tc deep/ref power transfer coefficient in Period I. than the Tc surf ace/ref in Period II and finally the estimation of the attenuation is accounted to the relation of the two coefficient. The amplitude attenuation has been derived as the square root of the power transfer coefficient of subsurface-surface site relation. The observations took 4 and 7 days in Period I. and II., respectively.
After resampling the three datasets to a common frequency of 800Hz both mobile data have been rotated in the coordinate system of the reference site. The magnetic components of the subsurface and the surface observation have been transformed by 22 • and 55 • , respectively. Since the infrastructure of the mine has been operating during the whole observation campaign, the time and frequency domain identification of low magnetic noise intervals has formed the basis of the forthcoming processing. The variation of the x component in time and time-frequency domain is plotted in Fig. 22. An intermittent broadband load with maximum at 50Hz has been demonstrated for the whole observation period and besides of that a peak appears at around 15-17Hz, too for several minutes by no regular timing. These large power load periods are related to the operation phases of the elevator of the mine shaft and the latter one is originated by the horizontal transportation and the lighting on active levels of the mine. Both have been present in the whole observation period. The low-noise segments have been identified based on the dynamic magnetic spectra covering Period I. Finally the remaining subsurface data consists of thirty one 80-100 second long time series of relatively quiet time segments, selected for further processing. The low noise windows are concentrated in the early afternoon hours of 24 March, 2018, 12:00-16:00 UTC. For the reason of having similar relative position between the observation site and the source area the same temporal pattern of time segments have been selected in Period II from the surface magnetic variation time series on 13 April, 2018. From the reference data the same temporal pattern of time sections have been selected for the estimation of the transfer functions in Period I and Period II, too.
Since the inside clock of the data logger has a characteristic delay per day, post synchronization of the subsurface data have been applied first. The method is based on the comparison of the same natural signal recorded on both stations: cross correlation function of the subsurface and the Hylaty data have been estimated on the 31 pair of selected time series sections. Preparation of magnetic data sections by the application of a 5.6-25Hz bandpass Butterworth filter was necessary to attenuate the strong effect of the 50Hz peak of the power grid on the cross correlation function. The time delay according to the cross correlation was 4 samples which  means 5ms for both channels. It has been corrected assuming linear drift of the inside clock time. The total magnetic time series has been correlated in both relation of the subsurface-Hylaty and surface-Hylaty observations. The cross correlation of the bandpass filtered total variations are overlay plotted for each section pairs in Fig. 23 after post synchronization. The cross correlated time sections contain minimum 64.000 samples each at 800Hz sampling rate. After the same processing, the cross correlation between surface registration and reference data has also been demonstrated, see Fig. 23. The sharp and high peak at zero demonstrates that the subsurface station observed the same natural signal as the reference station. In frequency domain the correlated components have also been confirmed by computing the spectral coherence between the magnetic variations in relation of subsurface-Hylaty and surface-Hylaty datasets, see Fig. 24. The spectral coherence plots indicate definite correlation at the first three Schumann resonance frequencies in relation of both subsurface-Hylaty and surface-Hylaty site.
Finally, four dataset has been prepared for the estimation of the transfer functions of subsurface-Hylaty and surface-Hylaty relation and finally to the determination of the surface-subsurface transformation: (1) B I.mine : 31 selected time segments in Period I, subsurface magnetic variation; (2) B I.Hylaty : the same 31 time segments in Period I, reference site magnetic variation; (3) B II.surf ace : 31 selected time segments in Period II with the same daily distribution, surface site magnetic variation; (4) B II.Hylaty : the same 31 time segments in Period II, reference site magnetic variation;

5.3.
Results. The introduced processing steps resulted in four dataset, which are accurately synchronised, transformed to the same coordinate system and resampled to a common sampling rate of 800Hz. The data has been prepared for time domain and spectral analysis. Additionally it has been confirmed that the subsurface observation and the Hylaty reference station has observed the same natural signal which variation is considered to be homogeneous in the scale of their relative distance, 300km. It has been demonstrated in the surface-Hylaty observation, too. The analysis of the distribution of the spectral power values at individual frequencies demonstrated non-normal distribution of the spectral power. The mean, the median and the interquartile range (IQR) of the spectral power have been plotted in Fig. 25.
In the first approach no outlier removal algorithm has been utilised, but the simplest robust statistics, the median of the spectral power has been determined at each frequency values. For the estimation of the signal attenuation with the depth, the median of the power spectral density (PSD) has been computed related to the All four median PSD functions have been determined from the Fourier transform of 1024 sample long consecutive time windows, with 512 sample overlapping between adjoining segments, using Hamming window. For an estimation of the power attenuation with the depth, the power of the fundamental Schumann spectral component has been determined in case of all four observations. The following subsection introduces a method which has been applied to have accurate estimation of the fundamental Schumann resonance spectral power peaks related to the four observations. 5.3.1. Polynomial baseline approximation. The method to provide an estimation of the ELF range attenuation is based on utilisation of an advanced background baseline removal technique. In this method the baseline is estimated by fitting a low-order polynomial, but instead of obtaining parameters by minimising the sum of squares of the residuals a non-quadratic cost function is adapted [37]. Third order polynomials have been considered to fit the median PSDs in a limited frequency range, see left four subplots of Fig. 27. All polynomials have been fitted by means of asymmetric truncated quadratic cost function.  The fitted baseline polynomials have been substracted from the power spectra, plotted on the right four subplots in Fig. 27. The power ratio obtained at 7.81 Hz in the subsurface-Hylaty relation is 1.126, for the surface-Hylaty relation 1.219 which results a power ratio of 0.924. The square root of the power ratio results the estimated amplitude attenuation: 96.1%. The rate of signal amplitude between the surface and at depth h can generally be computed by writing the exponential decay of the magnetic signal amplitude as follows: By substitution of the empiric amplitude rate values and the 140m depth of the subsurface station the characteristic depth of the amplitude attenuation, the so called the δ skin-depth can be estimated. It results δ = 3520 m. The power decrease of the signal with the depth is a consequence of the nonzero electrical conductance of the underlying/overhead andesite rock. Providing an estimation of the bulk resistivity of the andesite block is possible by a simplified approach of homogenious half space approximation. This estimation is considered valid if the homogenious halfspace assumption is realistic. A simple approximation formula between the skindepth and the bulk electrical resistivity is as follows [38]: Extracting the bulk resistivity by substitution of estimated skin-depth value of the Schumann fundamental resonance frequency 387 Ωm is obtained. The estimated resistivity value is in the wide characteristic range of electrical resistivity of the andesite rock in the literature (170-45.000Ωm) [39].

Time domain investigation.
A recent electromagnetic study of ELF range magnetic signal attenuation could have been explained by the assumption of distortion effect of iron artifacts of the mine infrastructure. For this reason the time domain inspection of the bandpass filtered subsurface data has been carried out. Bandpass filter is tuned below common power grid frequency and below upper Schumann frequencies to have comparable pattern of the hodographs. Based on this, 5.6, 10Hz cut off frequencies has been set. About 100 sample long time windows have been picked from the subsurface data and the corresponding reference time series. The magnetic variation of a time segments has been plotted and compared in Fig.  28.
The hodograph demonstrates no strong deviation of the subsurface signal is present confirming that iron parts of the mine infrastructure did not have recognizable effect on the subsurface observation and the estimated attenuation rate at the subsurface measurement site. The noisy environment, the amount of data utilised and the indirect comparation of the signal transformation all increase the deviation of the results and the variance of the power spectra. Considering all these circumstances the obtained attenuation value confirmes that in less noisy intervals, parallel measurement of surface and subsurface relation together can significantly improve the quality of the data and may provide more accurate estimation of the amplitude attenuation in the ELF range. Application of outlier removal may result further improvement on the quality of the estimation, too. The main goal of the study was to provide an estimation of the magnetic signal attenuation with depth. The results confirm that even an indirect estimation method can be promising. The estimation of the attenuation coefficient of the overhead andesite block might be more accurate by means of the data observed in an upcoming EM silent time window and even the more with the possibility of parallel recordings at the surface and subsurface sites. These two circumstances should improve the accuracy of the estimation significantly.
6. Infrasound measurements 6.1. Introduction. Sensitivity of second generation gavitational-wave detectors at low frequencies are limited fundamentally by the thermal noise of the suspension last stage and of the test masses, the seismic noise, and the gravitational gradient noise (so-called Newtonian noise (NN) ) resulted by changes in density of medium near the detector. Seismicity is one source of this type of noise, but changes in density of air caused by pressure waves are also contribute to NN [40].
Inspired by recent theoretical investigations of NN from the infrasound below the Earth surface [40], we have continued our infrasound measurements in MGGL. For these measurements, new infrasound microphone (ISM) was developed by the Institute for Nuclear Research of Hungarian Academy of Sciences (Atomki) [10]. The first instrument was installed at MGGL in 2016 with the purpose on one hand to test its performance on the long-term, and on the other hand to characterize the underground site's ambient infrasound noise.
In this study, Bowman's low-and median noise models were used as a reference for our site characterization. According to Ref. [41], the main sources of infrasound noise on the surface of Earth are microbaroms and the wind. These sources are not present underneath, but other sources related to the miner's work and machines produce infrasound. Our results show that ambient infrasound presumably can be supressed when going under the ground, but the local noise sources and the geometry of the underground site influences the possible extent of the noise supression at different frequencies.
6.2. The instrument. ISM consist of two main parts: an industrially manufactured capacitive sensor and C-V converter, developed by Atomki. The capacitive sensor is a differential pressure sensor with two inlets. One inlet is connected to the environment, while the other is connected to a reference volume. A diaphragm  Table 11. The characteristics of the two infrasound sensors, ISM1 and ISM2, developed in ATOMKI.
is placed in between the two pressure inlets. The diaphragm is deflected when the difference between the environment and the background pressure changes. The deflection is converted into an electrical signal by the sensor, then the signal is amplified to the −10 − 10 V range. There is a leak with a diameter of 250 microns, which connects the background volume to the environment. It makes possible to slowly equalize the background pressure and the environmental pressure, and hence determines the lowest frequency of the measuring range. The main characteristics of two versions of ISM are given in Table 11.
The drawback of the structure of the instrument is that fast changes in temperature dT /dt > 2 • C/h introduce artifacts in spectrum below 1 Hz. When the frequency of the temperature changes is in the measuring range of the detectors, then the temperature change initiated pressure in the background volume cannot be distinguished from the pressure change due to environmental infrasound. Therefore a high-resolution temperature sensor is built into the detector. 6.3. Data acquisition system and data processing. The microphones have analog output. In order to sample, digitalize and store measurement data, a lowcost data acquisition system is added to the sensors, based on a Raspberry Pi 3 model B [10].
The spectral properties of the infrasound background are characterised by the Pressure Amplitude Spectral Density (P ASD) and the analysis of the data is similar to the seismological noise evaluation in Section 3. The collected data was downsampled to 100 Hz, then divided into 163.84 s long segments. For each segment, PASD was computed. For each day, the 5 th , 50 th and 95 th percentile of corresponding PASD data was determined.
6.4. Comparative measurements with MB3 microbarometer. In order to compare performance of ISM1 with the performance of MB3 microbarometers [42], comparative measurements were performed on 2018-04-26 at Infrasound station no. 3 of Piszkéstető observatory in Mátra mountain range, Hungary. This station is part of the Hungarian National Infrasound Network, which is a permanent infrasound network operated by the Kövesligethy Radó Seismological Observatory of the Hungarian Academy of Sciences [43].
At Piszkéstető, the MB3 microbarometers are placed in a reservoir, which is partly covered by soil. MB3 is equipped with a wind reducing system of long pipes attached to the microbarometer. During our measurements, the pipes were detached, and the MB3 and ISM1 were placed next to each other in order to detect similar pressure as much as it is possible. The measurement lasted for three hours.
Our aim was to calibrate our instrument, that is to determine the real PASD values in Pa/ √ Hz from the PASD ISM 1 given in ADCunits/ √ Hz computed from measurement data. During the data processing part, raw signals both of MB3 and ISM1 were cut into 128 seconds length segments and PASD M B3 and PASD ISM 1 of each segment was computed. Then values corresponding to a given frequency of PASD M B3 and PASD ISM 1 pairs from the same time interval were collected. The number of pairs corresponding to a given frequency is equal to the number of segments used for calibration. For a given frequency, the ratios of PASD M B3 values and PASD ISM 1 values were computed (C=PASD M B3 /PASD ISM 1 ). The distribution of these ratios corresponding to a given frequency can be approximated with normal distribution. Standard deviation (σ C ) of the distribution was estimated by using median absolute deviation ( [44]). Two standard deviations were chosen to express the uncertainty of the ratio in Fig. 29.  6.6. Relationship between seismic and infrasound noise in the MGGL. In order to examine the relation between the signal of the seismometer and the signal of ISM1, we have computed the magnitude squared coherence estimate for each 128 s long data segment pairs from the deep measurement campaign of MGGL by using Welch's method [45]. In the case of the ET1H seismometer's signal, we treated the tree directions, East, North and Z (vertical) separately. This way, we got three datasets for the coherence: C v E p , C v N p and C v Z p . We determined the average of the coherence values of each frequency of the three datasets. A peak between 3 and 4 Hz can be observed on each plot on Fig. 31, which indicates a common source has an effect both on the seismometer and ISM1. There were no similar coherences with the farther, GU02 seismometer. The drainage system or the ventilation system may be responsible for the peak. We will prove or refute this hypothesis with the forthcoming experiments in MGGL.
7. Measurements of the inhomogeneity and the rock-density at large-scale with cosmic muon tracking A novel cosmic muon tomography technology were used to explore the precise, large-scale density structure of the Mátra mountain above the MGGL within the Gyöngyösoroszi mine. The idea is to measure the angular distribution of the highly penetrating muon-component flux of the cosmic-rays, which are created in the upper atmosphere and reach the surface of the Earth. This technique, has been applied successfully by our group for various cases such as cave researches, cosmic background measurements for underground laboratories, civil engineering, and volcano muography [47,48,49,50,51,52].
To reach the full, 2π coverage, muon flux measurements were continued, using a portable tracking system, called the Muontomograph. This has been developed by the Wigner Research Centre for Physics and described in details in Ref. [10]. Based on the Run-0 experiences, we made modifications on the Muontomograph in order to get better performance in the MGGL laboratory for the Run-1. The position of the MGGL and the acceptance directions of the detector related to the surface topology are shown in Fig. 32. We adjusted the trigger levels for higher efficiency and a noisy chamber were also removed, thus the number of chambers were reduced to five.  Since we have covered the zenith-and two 90 • -tilted directions in Run-0, the four Run-1 measurements were all tilted to 45 • from the zenith. With this choice we could cover the upper hemisphere fully, without any missing direction. The first measurement was directed to 20.5 • for 77 days, then we continued to 110.5 • , 200.5 • , and 290.5 • directions, respectively for 62 days, 77 days, and 64 days. Further details of the measurements are summarized in Table 12. Data were monitored and downloaded continuously from the detector via Ethernet connection during the data taking period in MGGL for online checking of data quality. Comparing Run-0 and Run-1 date in Table 12, one can see that, although the event number dropped by an order of magnitude, the fraction of the track number increased relatively to the number of events. This means, that the efficiency become better since noise and fake events were suppressed well in Run-1.

7.2.
Results of the rock inhomogeneity measurements. The muon flux distribution can be obtained after applying a track reconstruction algorithm and merging the maps according to the partial overlap and geometry normalization [51]. The Run-0 and Run-1 flux maps are shown in the left and right panels of Fig. 33 respectively. The muon flux distribution is plotted by color-scale contours as a function of azimuth and zenith angles. The cosmic muon (track) rate was 0.005-0.02 Hz for both runs, which is sufficient to provide flux measurement with 5 − 50% statistical errors (depending on zenith angle). The measured flux has a maximum value of 0.7 m -2 sr -1 s -1 at 15 • zenith angle to the West. In Fig. 33 white color 'triangles' at 20 • zenith angle (panel left) and 'dial' at zenith (panel right) represents the directions, which are out of the acceptance of the current run.
The detector-to-surface distance is indicated with dark contour lines in the Fig. 33. For Run-0 in [10], we used Shuttle Radar Topography Mission (SRTM) satellite data for the elevation map with an estimated relative error of 10% (10 m at the zenith). However by comparing the SRTM model with other available elevation data the accuracy of the SRTM model seemed much poorer than it was previously expected, likely because of the different height of the vegetation. Thus the elevation contours obtained from the National Hungarian Grid (EOV) have been used instead. Applying this, new detector-to-surface data, give better correlation for both Run-0 and Run-1. Merging all the data of Run-0 and Run-1 the muon flux distribution can be obtained, which is plotted in Fig. 34 up to zenith angle 60 • . Correlation between muon flux map and the detector-to-distance contour lines is excellent. We found fluctuation is in the order of ±0.05 m -2 sr -1 s -1 as it was obtained in Run-0 as well.
One can observe that the surface distance and the muon flux reasonably match each other in Fig 34. Combining this with the acceptance plotted on the surface map plot on Fig. 32 the maximal flux matches the shortest distance-to-surface directions on the hill's West slopes, while to the East, the muon flux drops at the large zenith angles with the longest traveling lengths. The collected ∼ 340, 000 tracks of the muon flux measurement provided highenough statistics to estimate the rock thickness, based on the muon absorption model [55]. This model has been used successfully in other, similar measurements up to few 100 meters depth. In Fig. 35 the angular distribution of the rock density is presented. The average of the rock density is 2.6 ± 0.1 kg dm -3 , which is the density of the andesite-based host rock. The data do not show any large-scale density inhomogeneities or cavities above the MGGL. However, minor, ∼ 0.2 kg dm -3 , positive density deviation ridge is observed along the 40 • − 50 • zenith angle from the southwest to northwest. This may correspond to several percent metallic ore in the unexplored part of the Gyöngyösoroszi mine. A negative ∼ 0.5 kg dm -3 density anomaly valley was also observed at the lowest 50 • − 60 • zenith angle from the south to north. This is the direction of the known tunnels and caverns of the mine, although here the uncertainty of the measurement is also significant at large zenith angles. Figure 35. The angular distribution of the density based on Run-0 and Run-1 measured from the MGGL is plotted as a function of azimuth and zenith angles from the detector position. Color-scale contours show the rock densities, dark contour lines show the detectorto-surface distance in meters.

7.3.
Conclusions. Based on the muon-flux measurements, the improved detectorto-surface date, and the muon absorption model, we can state, that the rock density map correlates well with the unexplored and mined areas above the MGGL. This unique measurement support to place the Einstein Telescope at Gyöngyösoroszi mine, since the lack of large-scale density inhomogenities would reduce the higherorder tensor-like gravitational corrections during the measurements.

Summary
The Mátra Gravitational and Geophysical Laboratory was established to investigate on a long term basis the conditions and requirements of next generation gravitational wave detectors in case of underground construction and operation. In more general, the aim was to measure various geological and rheological properties of the base rock in addition to test experimentally novel theoretical approaches on the noise penetration and suppression. Also the Mátra mountain range is studied as a possible site of the planned Einstein Telescope.
Our first investigations were published in Ref. [10] including the technical details of the measurements. In this paper we summarize the first two years of measurements of the laboratory by different methods: geophysical environment, electromagnetic attenuation, infrasound noise, cosmic muon tomography of the sorrounding rock mass, and long term seismic noise. The timeline of the data taking periods of the various measurements are shown in Fig. 36. Our topical results are: Geophysical environment: A short survey of the geological and seismological conditions in the Mátra mountains reveals a homogeneous composition of hard andesite rock with low seismic activity. We performed rock mechanic experiments and shown, that a typical grey andesite from the vicinity of the laboratory is not ideal elastic. We characterised and modeled the deviation from ideal elasticity and determined the corresponding rheological parameters (see Tables 3 and 4). The independent measurements of dynamic elastic moduli verify our novel rheological model, which properties are important when calculating the Newtonian noise from rock deformations. Electromagnetic attenuation: With EM noise measurements we estimated the electromagnetic attenuation by the andesite rock mass in the lower ELF range, especially at the frequency of the first Schumann resonance component. Comparing the data of the external surface reference station, we obtained the skin depth 3520 m. For the bulk resistivity at the Schumann resonance frequency 387Ωm was obtained. This value -supporting the validity of this measurement -fits well to the literature value of vulcanic, andesite rocks (170-45000 Ωm). Infrasound noise: Using the custom-made infrasound detector designed by the MTA Atomki, we measured the pressure amplitude spectral distribution in the MGGL (see Fig. 30). This is relatively noisy in the frequency range 1 − 7 Hz. The recorded infrasound noise was found to be in accordance with the seismic noise at 3-4 Hz in the Laboratory. Cosmic muon tomography: During the 404 days long measurement with the Muontomograph merging Run-0 and Run-1 data, led us to map the rock density and its inhomogeneities above MGGL at large scale (see Fig. 35).We verified this novel measurement technique by obtaining the typical andesite rock density (2.6 ± 0.1 kg dm -3 ). According to this investigations, we have found, there were no large scale density inhomogeneities (≤ 0.2 kg dm -3 ) measured in the rock mass. Seismic noise: The long term seismic noise was registered by two seismometers inside the MGGL at depth −88 m almost for the whole two-years period.
In parallel, as a cross check, a third seismometer measured the seismic noise for two-weeks duration. This data record were farther from the MGGL in the mine in −404 m depth. Seismic noise measurement were found to be consistent.
In the MGGL we had the opportunity to combine data taken at the same position, and in addition to this we could test the insights gained from noise analysis. According to our experience the proper evaluation of long term data requires some refinements of the methods applied in the ET Design Report [9]. Our suggestions has been published in Ref. [28] and in the recent study we evaluated the data accordingly. Therefore we calculated and show also the median of the spectra and the integrated noise measures rms 2−10Hz and rms 1−10Hz from the median, here. The analysis of the data from the Guralp ET1H and GU02 stations is based on these characteristics and the evaluation of the data from the WARS station also considered the problems of the more traditional methods.
We found, that the long term data does not show yearly changes in the median, but the intensity of the noisy periods reflects the changes in the activity inside the mine (see Fig. 10 and Tab. 7). The seasonal changes show minimal noise level in the late spring early summer period (see Fig. 11). The comparison of the noise level at the shallower and the deeper locations revealed that the two-weeks measurement is unexpectedly representative for the whole Run-0 and Run-1 period (see Fig. 13). We compared also working and night periods in MGGL in order to estimate an achievable minimum noise level (see Fig. 9). Considering our whole analysis the median spectrum of the night noise at -404m in the mine in Fig. 37 is convincingly representative in this respect. Here the percentiles are calculated from the spectra of the 50s long periods of data. One can see that in 90% of these periods the noise level is below the Black Forest line.
In general the noise levels our long term measurements in the shallower MGGL and also the two-week data from the deeper location are remarkably similar to the noise levels of the previous study [9,25,26], the mode rms 2Hz are almost the same. As we have already mentioned, the median rms values are preferable, because they are more stable, less sensitive to the uncertainties is the data and in the evaluation. The ET1H and WARS median rms 2Hz are close to each other. However, from the point of view of operational requirements, a higher percentile limit is a better estimation of the sensitivity of the detector for continuous operation. Therefore, the 90 th percentile rms values, given in Tab. 13, estimate the noise level for a continuous operation of the Einstein Telescope.
According to the recent study a reliable and comparable characterisation of long term seismic noise data requires the evaluation of various noise measures with a uniform methodology. The separation of reducible and irreducible cultural noise, the seasonal changes and the attenuation by depth are important aspects for site description. Also it is very important to publish the raw seismological data for open   Table 13. The rms values of the 90 th percentile in the E direction with different frequency ranges for the MGGL seismometer ET1H and for the deeper location GU02, too. The first two columns are calculated from the 90 th percentile of the 300 s data of the two-week interval, the last one from the average of daily 90 th percentiles. evaluation and comparison. For the ET1H station of MGGL and for the GU02 station these are available at [56].

Conclusions
During the Run-0 and Run-1 periods of the MGGL in 2016-2018 we performed long-term monitoring the Mátra mountain range as a possible underground site of the planned Einstein Telescope. We used various standard methods in parallel to novel approaches of investigating the geophysical environment, electromagnetic attenuation, infrasound noise, cosmic muon tomography of the surrounding rock mass, and long term seismic noise. The collected data could enable us to cross check and and compare standard measurements and techniques applied in earlier investigations with the new ones. Alongside this, the geological and rheological properties of the base rock were summarized in this paper. In addition to the analysis of the noise background relevant for a next generation, underground-based gravitational wave detector, especially in the low frequency regime, at 1-10 Hz. We strongly believe that applying our results for the site selection will significantly improve the signal to nose ratio of the multi-messenger astrophysics era. Our conclusion was that for the background noise analysis, it is necessary to perform long term data taking and apply the state-of-the-art techniques presented here.