René Steinmann
email: rene.steinmann@univ-grenoble-alpes.fr

Céline Hadziioannou

Eric Larose

Effect of centimetric freezing of the near subsurface on Rayleigh and Love wave velocity in ambient seismic noise correlations

Keywords:

About a decade ago, noise-based monitoring became a key tool in seismology. One of the tools is Passive Image Interferometry (PII), which uses Noise Correlation Functions (NCF) to retrieve seismic velocity variations. Most studies apply PII to vertical components recording oceanic low-frequent ambient noise ( < 1 Hz). In this work, PII is applied to high-frequent urban ambient noise ( > 1 Hz) on three three-component sensors. With environmental sensors inside the subsurface and in the air, we are able to connect observed velocity variations with environmental parameters. Temperatures below 0 • C correlate well with strong shear wave velocity increases. The temperature sensors inside the ground suggest that a frozen layer of less than 5 cm thickness causes apparent velocity increases above 2 %, depending on the channel pair. The observations indicate that the different velocity variation retrieved from the different channel pairs are due to different surface wave responses inherent in the channel pairs. With dispersion curve modelling in a 1D medium we can verify that surfaces waves of several tens of meters wavelength experience a velocity increase of several percent due to a centimeters thick frozen layer. Moreover, the model verifies that Love waves show larger velocity increases than Rayleigh waves. The findings of this study provide new insights for monitoring with PII. A few days with temperature below 0 • C can already mask other potential targets (e.g. faults or storage sites). Here we suggest to use vertical components, which is less sensitive to the frozen layer at the surface. If the target is the seasonal freezing, like in permafrost studies, we suggest to use three-component sensors in order to retrieve the Love wave response. This opens the possibility to study other small scale processes at the shallow subsurface with surface wave responses.

Introduction

For a long time in the history of seismology ambient seismic noise disturbed continuous recordings and seismologists developed methods to eliminate it. In the last two decades this has changed and ambient seismic noise became a valuable source of information when studying the Earth. Nowadays it can happen that signals from earthquakes disturb the recordings of ambient seismic noise and have to be removed from the data. In these cases, either the source of ambient seismic noise is the object of interest or the ambient wave field carries precious information from the subsurface. A method to retrieve the Green's Function (GF) of a medium from noise-based recordings is Passive Seismic Interferometry (PSI).

The cross correlation of two ambient seismic noise recordings from two different stations can converge towards the Green's Function of the medium [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF][START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF]. Due to the permanent presence of ambient seismic noise the retrieval of the GF is arbitrarily repeatable in time. This opened a new way in seismology to monitor the elastodynamic properties of the Earth. One method to monitor the seismic velocity is called Passive Image Interferometry (PII). PII combines PSI with Coda-Wave-Interferometry (CWI), which was first mentioned by [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[END_REF]. [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF] proposed the method of PII and retrieved velocity variations at the Merapi Volcano, which correlated well with changes of the ground water table .   With the new monitoring possibility the field of environmental seismology was born. Environmental seismology uses seismic wave fields and natural seismic sources to study the coupling between the Earth and the external environment [START_REF] Larose | Environmental seismology: What can we learn on earth surface processes with ambient noise[END_REF]. Erosion, thermal forcing, precipitation, freezing etc. affect the surface of the Earth and, thus, modify the propagation of a wave field at the near subsurface. In the last years many studies used ambient seismic noise to monitor different environmental effects on seismic velocities. [START_REF] Richter | Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry[END_REF] could track the daily and seasonal temperature variations inside the first few meters of the subsurface with PII. [START_REF] James | Improved moving window crossspectral analysis for resolving large temporal seismic velocity changes in permafrost[END_REF][START_REF] James | Insights into permafrost and seasonal active-layer dynamics from ambient seismic noise monitoring[END_REF] used ambient seismic noise to monitor the freezing and thawing cycle of permafrost soils in the Arctic. Prior to a landslide event [START_REF] Mainsant | Ambient seismic noise monitoring of a clay landslide: Toward failure prediction[END_REF] observed a decrease of seismic velocity at the base of the sliding layer with PII. [START_REF] Lecocq | Monitoring ground water storage at mesoscale using seismic noise: 30 years of continuous observation and thermo-elastic and hydrological modeling[END_REF] used 30 years of continuous ambient seismic noise recordings of the Gräfenberg array to track hydrological changes in the aquifer.

These studies showed that noise-based monitoring is able to detect minor changes in the near subsurface, making it a good method to study the coupling between the Earth's near subsurface and the external environment. However, often multiple effects overlap and it can be hard to distinguish between them. For example, observed seasonal variations in the seismic velocity can have several causes, e.g. temperature, water content or changes of the noise sources (e.g. [START_REF] Tsai | A model for seasonal changes in GPS positions and seismic wave speeds due to thermoelastic and hydrologic variations[END_REF]; [START_REF] Zhan | Spurious velocity changes caused by temporal variations in ambient noise frequency content[END_REF]; [START_REF] Lecocq | Monitoring ground water storage at mesoscale using seismic noise: 30 years of continuous observation and thermo-elastic and hydrological modeling[END_REF]). [START_REF] Gassenmeier | Monitoring of environmental influences on seismic velocity at the geological storage site for CO2 in Ketzin (Germany) with ambient seismic noise[END_REF] monitored a geological storage site in Germany and seasonal freezing masked the original target. Thus, the detection and understanding of environmental effects are also necessary to remove these signals and focus on deeper targets.

The study, presented here, uses the method of PII in the urban setting of Hamburg to monitor the coupling of the near subsurface with the external environment (see Figure 1). The measuring site is equipped with three seismometers (WM01, WM02 and WM03), temperature sensors in the subsurface and the air, water content sensors in the subsurface and rain gauges. First, the nature of the urban noise field is investigated to assess how well it performs within the framework of PII. Then, relative velocity variations are retrieved and linked to environmental parameters. The last section covers dispersion curve modelling to provide explanations for the observations.

Methods

Passive Seismic Interferometry -PSI

The method of PSI is briefly introduced. For a more detailed explanation we refer to [START_REF] Wapenaar | Tutorial on seismic interferometry: Part 1-basic principles and applications[END_REF]. PSI turns passive seismic measurements, e.g. ambient seismic noise recordings, through cross correlation into deterministic seismic responses, which converge under certain circumstances to the GF of the investigated elastic body. The two stations, which deliver the noise recordings for the cross correlation, need to be evenly surrounded by uncorrelated noise sources. The recordings have to be long enough compared to the dominant period of the noise sources. If these requirements are met, the obtained function from now on called Noise Correlation Function (NCF), can be written as:

N CF = u(X B , t) * u(X A , -t) = [G(X B , X A , t) + G(X A , X B , -t)] * S N (t) (1)
u(X A , t) is the recording of a receiver at X A and u(X B , t) is the recording of a receiver at

X B . G(X B , X A , t)
is the GF with a receiver at X A and a source at X B . G(X A , X B , -t) is the time reversed version of G(X A , X B , t) and S N (t) is the auto correlation of the ambient noise.

Here, the asterisk denotes temporal convolution and . denotes the long term average. Thus, u(X B , t) * u(X A , -t) is the cross correlation of the two recordings u(X B , t) and u(X A , t).

If the noise sources are not evenly distributed, the Green's Function can be retrieved only partly or not at all. Noise sources on Earth, such as oceanic microseism or human-related activity, are far from isotropic and thus violate the principles of PSI. However some studies showed that the isotropic illumination by noise sources may be violated to a certain degree depending on what the purpose of the NCFs will be. For example, [START_REF] Lin | Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array[END_REF] retrieved the vertical Rayleigh wave response of a virtual source with uneven distributed noise sources and [START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments[END_REF] showed that for monitoring temporal changes it is sufficient that the source distribution is stationary.

When three-component sensors are used, nine cross correlations are possible: each other [START_REF] Roux | Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA[END_REF]. All this only holds true in an isotropic medium and when the noise sources are perfectly distributed. If this is not the case, the first arrivals of Rayleigh and Love waves would mix on all channel pairs except for Z-Z, since the vertical component is not sensitive to Love wave motions.

Z-Z, Z-R, Z-T, R-Z, R-R, R-T, T-Z, T-R
Many studies use only the vertical components for PSI and disregard the horizontal components, because of low Signal-to-Noise-Ratio (SNR) of the NCF produced with horizontal components [START_REF] Roux | Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA[END_REF]. In this work we use all three components and show that crucial information about noise source characteristics and velocity variations can be retrieved if all components are considered.

Retrieving relative velocity change with the stretching method

A common concept to measure a change of the seismic velocity is to use the same sourcereceiver-setting before and after the change occurred and compare the two signals with each other. Obviously the characteristics and positions of the source and the receiver need to be exactly the same, otherwise changes in the signal can be wrongly interpreted. [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[END_REF] were the first to conduct this experiment with earthquake doublets. Earthquake doublets occur in the same source region and therefore share the same characteristics and position to some degree. They found out that time shifts in the coda wave were caused by small velocity variations, which were not visible on the direct P-and S-waves. The coda of a seismic record contains scattered waves travelling through the perturbed area several times. Therefore, it is more sensitive to slight changes. [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[END_REF] interpreted the retrieved velocity variations as a consequence of variations in the subsurface shear wave velocity. [START_REF] Snieder | Coda wave interferometry for estimating nonlinear behavior in seismic velocity[END_REF] showed that the time shift dt between individual wiggles in the signals recorded before and after a perturbation of the velocity, increases linearly with the lapse time t, resulting in a constant relative time shift dt/t. They also add that the ratio of time shift to lapse time is equal to the negative of the spatially homogeneous relative velocity change:

dt t = - dv v (2)
This implies that one can retrieve dv/v by simply measuring the relative time shift dt/t. [START_REF] Snieder | Coda wave interferometry for estimating nonlinear behavior in seismic velocity[END_REF] called this concept Coda-Wave-Interferometry. The stretching method, which is used in this study, is one of several methods to estimate the time shift [START_REF] Mikesell | A comparison of methods to estimate seismic phase delays: Numerical examples for coda wave interferometry[END_REF]. Let u 1 (t) represent a seismic response recorded before the velocity change dv took place, and u 2 (t) a seismic response later in time after the velocity change dv took place.

u 2 (t) can be stretched or compressed in time by a constant factor until u 2 (t(1 -)) is brought into the same state as u 1 (t). The stretching factor is in fact the time shift dt/t.

The stretching method performs a grid search over i maximizing the correlation coefficient CC( i ):

CC( i ) = t 2 t 1 u 2 (t[1 -])u 1 (t)dt t 2 t 1 u 2 (t[1 -])dt • t 2 t 1 u 1 (t)dt (3)
Since the method is rather applied on a time window than on the whole signal, the start-time t 1 and end-time t 2 of the time window have to be defined. The stretching factor yielding to the maximum of CC( i ) equals the negative velocity perturbation. Depending on the position of the time window, which is stretched, the retrieved velocity variation can be connected to a certain type of seismic velocity. Most studies focus on the coda and connect their results to the shear wave velocity.

As laid out in the previous remarks the method depends on having a reproducible source.

Earthquakes are far from reproducible and if they are, the time in between an earthquake doublet is unpredictable and random. Explosions and other man-made sources can be costly if the monitoring is supposed to deliver values every day or the monitoring system is set up for several years. PSI fills in the gap here. Depending on which noise sources are available, one can produce NCFs on a daily or even hourly basis. The total amount of NCFs is stacked to a Reference Correlation Function RCF. Instead of the two recordings u 1 (t) and u 2 (t) from equation 3 we have a RCF serving as a reference to every NCF. The grid search over i would then look like:

CC j ( i ) = t 2 t 1 N CF j (t[1 -])RCF (t)dt t 2 t 1 N CF j (t[1 -])dt • t 2 t 1 RCF (t)dt (4)
j denotes the index over all NCFs and therefore a value for dv/v is retrieved for every NCF.

Here, we perform the stretching on the causal and acausal side and then average the obtained stretching factor. [START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments[END_REF] showed with laboratory experiments that the convergence of the NCF towards the GF is not necessary for this method. However, stable noise sources in time and space are important to create stable comparable NCFs. In general more noise in temporal terms improves the stability of the NCF. Though, cross correlating months-long noise recordings is computationally expensive and time consuming. Due to the linearity of the cross correlation the common practice suggests to cross correlate smaller time windows and then stack the NCFs until the stack reaches a stable state. Therefore the temporal resolution of the monitoring system is limited by the amount of data needed to create stable NCFs. Mainly the frequency content and the temporal and spatial stability of the noise sources influence how fast this stability is obtained. If one is not careful about this, the obtained stretching factor can be biased severely. [START_REF] Bussat | Using ambient noise for real-time (hourly) overburden monitoring[END_REF] presents an example, where a velocity change of 2 % within 4 hours was not caused by changes in the subsurface but by a passing storm. This raises the need for tools to estimate the stability of the NCFs before using them for the stretching method. These tools are introduced in section 5.

The combination of PSI and the stretching method is able to build a monitoring system to detect relative velocity changes. The next section introduces a way to estimate the corresponding error to the stretching results.

Error estimation for the stretching method

Changes in the noise sources can introduce a waveform dilation in a NCF and, thus, be misinterpreted as physical changes in the subsurface. [START_REF] Weaver | On the precision of noise correlation interferometry[END_REF] developed an expression to estimate this kind of error:

rms( ) = √ 1 -CC 2 2CC 6 π 2 T ω 2 c [t 3 2 -t 3 1 ] (5)
T is the inverse of the frequency bandwidth, in which the NCFs are filtered, t 1 and t 2 are the start and end time of the defined window for the stretching method, CC is the correlation coefficient obtained by the stretching method and ω c is the central frequency. Furthermore, if the results are averaged over N station and M channel pairs, rms( ) is divided by √ M N .

We use equation 5 to verify the apparent velocity changes we retrieve in this study.

Data acquisition and preprocessing

The data used in this work, were recorded by three three-component Trillium Compact broadband 120 s seismometers and sampled with 200 Hz by Reftek recorders. While those broad band sensors are in general used to investigate deep structure, here we use only the high frequency (short period) content for subsurface applications, such that other cheaper sensors could be used with equal sensitivity for our application.

The seismometers are all placed on the surface. While WM03 is inside a shed, WM01 and WM02 are placed outside underneath a cover. To ensure a better coupling for WM01 and WM02, we dug a 30 cm deep hole, filled it with gravel, put a slab on top and installed the seismometers on top of it. The data set ranges from the 5th of January to the 14th of April of the year 2018. The data is recorded in UTC-time and the time zone of Hamburg is UTC+1. When we mention times in the following remarks, we always refer to the UTC-time.

After retrieving the data from the stations, some preprocessing steps have been carried out.

Firstly the data was downsampled to 100 Hz to save computation time for all the following work. The aliasing effect was avoided by a lowpass filter with a corner frequency of 20 Hz.

Then, the mean value and the linear trend were removed. The instrument response was not removed, since the response function is flat from 0.01 Hz to 100 Hz. Most processing was carried out with the Python package Obspy [START_REF] Beyreuther | Obspy: A python toolbox for seismology[END_REF].

Local Noise Field

In order to characterize the local noise field a spectrogram of a week in February 2018 of a horizontal and vertical component is shown in Fig. 2. Very low frequencies (below 0.1 Hz)

seem to be very prominent on the horizontal component. This is probably due to tilting effects. The mircoseism around 0.1 Hz is also distinctive on both components. Frequencies higher than 0.6 Hz show the typical rhythm of an urban environment: higher amplitudes and nights. This is favourable for creating stable NCFs. From the map in Figure 1 we can infer the possible noise sources: the highway A1, train tracks, roads and a gravel pit.

Since we want to investigate the effect of freezing of the soil, high frequencies are favourable due to their high sensitivity to the near subsurface. However, above 10 Hz the variation in amplitude increases and the vertical component indicates monochromatic noise sources, which appear only during working hours. Temporal instability and monochromatic noise sources are highly unfavourable when performing PII [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF].

Preliminary Investigations

When applying PII and interpreting their results certain investigations have to be done beforehand. First of all, the NCFs have to be filtered in a frequency band, which is suitable for PII. In the following remarks two tools are introduced to check on the temporal stability Before performing the cross correlation the preprocessed data was processed another time.

First the data were rotated from the vertical-north-east-system ZNE into the vertical-radialtransversal-system ZRT. Afterwards, the daily files were sliced into 10 min windows, which then have been tapered to mitigate spectral leaking for the following steps. A bandpass filter between 0.3 Hz and 20 Hz was applied. To exclude earthquakes and other transient events, amplitudes higher than the sixth fold of the standard deviation have been clipped. Then the slices were spectral whitened between 0.3 Hz and 20 Hz. After the processing was done, a cross correlation in the time domain with a shift of 6000 samples -which equal 60 seconds with a sampling rate of 100 Hz -has been performed. In the next section we use the tools SNR and WFC to infer some characteristics of the NCF.

Waveform Coherency (WFC)

How is the WFC calculated? Firstly, we stack 144 NCFs, which equal one day of data and, thus, the stack is named daily NCF. Then, the CC between the daily NCF and the RCF for specific time windows on the time lag are calculated. This is done for every daily NCF and then we average the CC-matrix over all days. The WFC is the CC between the NCF and RCF averaged over the investigation period. As mentioned before, the broadband NCF contains frequencies between 0.3 Hz and 20 Hz. To investigate the stability of different frequencies, the NCF is filtered with octave filters with a minimum frequency of 0.3 Hz and a maximum frequency of 20 Hz. The WFC, depending on the center frequency f c of the octave filter and the center t c of the sliding time window on the time lag, can be written as:

W F C(t c , f c ) = 1 N N k=0 tc+∆t(fc) tc-∆t(fc) N CF fc,k (t)RCF fc (t)dt tc+∆t(fc) tc-∆t(fc) N CF fc,k (t)dt • tc+∆t(fc) tc-∆t(fc) RCF fc (t)dt (6)
N is the number of the daily NCFs and the time window ranges from t c -∆t(f c ) to

t c + ∆t(f c ).
In our case, t c ranges from -20 s to 20 s with an interval of 0.5 s, because with the greatest interstation-distance of 153 m we do not expect coherent signals after 20 s. With ∆t(f c ) we define the time windows around t c according to the frequency content: ∆t(f c ) is set to the largest period in the octave filter. This way, we can also compare the WFC values between different frequencies. The WFC ranges from 0 to 1. 0 means that the waveform is totally unstable and random in time. The stretching method would not give reasonable results. A value of 1 means that the waveform is absolutely coherent in time and no changes at all occur. Optimal would be a value slightly below 1 indicating that there are slight changes in the coherent waveform induced by possible velocity perturbations.

Figure 3 shows the WFC averaged over the channel pairs for each station pair, respectively.

For all frequencies and station pairs the most coherent signals are centered asymmetrically between -5 s and 5 s. The asymmetry suggests an uneven distribution of noise sources around the stations. Since the station pairs have a similar inter-station distance, the time lag positions of coherent signals are similar between the station pairs. Furthermore, lower frequencies produce higher CC values than higher frequencies. With later time lag incoher-ent waveforms start to dominate, and thus CC decreases.

5.2 Signal-to-Noise-Ratio (SNR)

With the help of the WFC-analysis we found out at which time lags the signal is coherent. With this information, we can now calculate the SNR of the NCFs for the different frequencies. There are many definitions for the SNR and in this study we define the SNR as the power ratio between two time windows containing the signal and the trailing noise, respectively. It can be written as:

SN R = t 2 t 1 |N CF (t)| 2 dt t 4 t 3 |N CF (t)| 2 dt (7)
t 1 and t 2 are the boundaries for the window containing the signal and t 3 and t 4 are the boundaries for the noise window. Since we are not interested in the absolute value of the SNR, but in the comparison of the different frequencies, t 1 , t 2 , t 3 and t 4 are fixed for all frequencies. The WFC-analysis showed us that the NCFs are not symmetric. Consequently, the SNR is calculated for the causal and acausal side. The signal window on the causal side is defined from 0.1 to 5 s and on the acausal side from -5 to -0.1 s. The noise window on the causal side is defined from 50 to 55 s and on the acausal side from -55 to -50 s. Instead of using a daily NCF, the RCF is used, since it represents the whole data set. We use the same octave filters as in the WFC-analysis.

Figure 4 depicts the SNR for all nine channel pairs averaged over all station pairs. Almost all nine channel pairs reach their maximum around 3 Hz. Around this frequency, the channel pair T-T even surpasses Z-Z, which mostly shows the highest SNR compared to the other channel pairs. This could be an indication for a high ratio of Love waves in the local noise field at this frequency. The SNR-analysis suggests to filter the NCFs around 3 Hz, since almost all channel pairs peak around this frequency. The WFC-analysis also shows values around 0.8 at this frequency suggesting that the waveforms do not only have a high SNR but are also stable in time.

Therefore, we propose to filter the NCFs between 1 and 6 Hz. In the following remarks, we examine how fast the NCFs reach a stable state in this frequency range.

Convergence of the NCF-stack towards a stable state

Variations in the noise sources can cause spurious velocity perturbations when applying the stretching method [START_REF] Zhan | Spurious velocity changes caused by temporal variations in ambient noise frequency content[END_REF][START_REF] Bussat | Using ambient noise for real-time (hourly) overburden monitoring[END_REF]. Stacking the NCFs mitigates the effect of variation of the noise sources and, thus, is crucial to eliminate spurious velocity perturbations by instability of the noise sources. Here we investigate how the SNR and WFC develop with the amount of stacking. We watch how these two parameters develop in a course of a week. For the SNR, we use the same signal and noise windows as in the previous section. The SNR is calculated every time we add a single NCF to the stack, thus, we see the development of the SNR with the stacking. The WFC between the RCF and the NCF-stack is also calculated every time we add a NCF to the stack. For the WFC, the time window, containing the signal, is defined from 0.1 to 5 s and from -5 to -0.1 s.

Figure 5 illustrates the development of SNR in red and WFC in blue with the stacking time. In the long-term, both parameters increase with stacking time. However, both do not increase constantly. The increase is most intense with the first stacks and the WFC reaches and stays over 0.8 after a stacking time of 1 day, equalling 144 NCFs. The WFC decreases at certain periods, which seem to coincide the working hours. This is a problem for the stretching method, which stretches the signal to increase the WFC. The SNR also shows daily variances in the increase. These features are likely caused by the variation of the urban noise sources. The spectrograms from Figure 2 illustrate variations in the power at the same rhythm.

Monitoring the subsurface on a daily basis does not seem appropriate here. A stack of a week would probably include all temporal variations of the urban noise sources, but the temporal resolution of the PII would be very low. Therefore, we suggest a stacking time of 4 days, which is a compromise between temporal resolution and stability of the NCFs.

The NCF time series

Since the stacking time and the frequency band are set, time series of the NCFs can be calculated. Figure 6 During that time, the temperature at 5 cm depth is not below 0 • C. Since the CWC shows no large variations, the temperature seems to be the cause of the velocity increase. Air temperature below 0 • C causes the subsurface to freeze and increases the rigidity at the surface. This process increases the shear wave velocity [START_REF] Zimmerman | The effect of the extent of freezing on seismic velocities in unconsolidated permafrost[END_REF]. [START_REF] Gassenmeier | Monitoring of environmental influences on seismic velocity at the geological storage site for CO2 in Ketzin (Germany) with ambient seismic noise[END_REF] observed a similar behaviour with noise-based monitoring. However, in our case, the frozen layer does not extend beyond 5 cm and increases the velocity by 2 %. To understand the process better, Figure 9 shows the temperature at the surfaces and at 5 cm depth with a sampling rate of 10 min together with the average over all channel pairs of the dv/v-curve.

We see that the freezing process is not a constant process but interrupted by warmer periods during the day. In the dv/v-curve this diurnal pattern is not visible, since we stack the NCFs over a period of 4 days. However, the trend of the dv/v-curve indicates whether the freezing or thawing process is dominating and whether the extent of the frozen layer is increasing or retreating. Besides the main two velocity increases there are two weaker velocity increases corresponding with consecutive cold nights: the first one is located directly at the beginning of the time series. The second velocity increase of roughly 0.5 % happens between mid March and beginning of April. Here we have 5 consecutive nights with temperatures below 0 • C.

Moreover, the temperature at 5 cm shows that the cold temperatures could not penetrate as deep as during the two colder periods around mid February and beginning of March. This observation suggests that the dv/v measurement is sensitive to the freezing process on a centimeter scale.

Since the stretching window focused on the first arrivals of the surface waves, the surface wave velocity seem to be sensitive to such a thin layer of frozen soil. Z-Z correlograms, containing the Rayleigh wave response, show almost no velocity increase during the cold period. The other channel pairs seem more sensitive to the freezing process than Z-Z. A possible explanation could be the different surface wave responses inherent in the different cross correlations. One would expect that the Rayleigh wave response is represented by Z-Z, Z-R, R-Z and R-R and the Love wave response is represented by T-T. However, with an uneven distribution of noise sources, as we assume in our case, it is possible that Rayleigh and Love wave motions are simultaneously present on the radial and transversal components.

Consequently, Z-Z would be the only channel pair with a clean Rayleigh wave response. To verify that the thin layer of frozen soil is really the cause of the velocity increase of surface waves, we model the dispersion curves of Rayleigh and Love waves in the next section. 7 Simulation of dispersion curves

Setup of the model and first results

The simulation of the dispersion curves is done in Geopsy [START_REF] Wathelet | Geopsy software manual[END_REF]. For the subsurface we assume a 1D layered model. Several SH-refraction profiles were recorded to derive the shear wave velocity and thickness of the layers. Soil samples from the first meter help to estimate the densities at deeper layers. The P-wave velocity is estimated through the Poisson's ratio ν. Unfortunately, the Poisson's ratio itself can only be guessed. During the time of the freezing, the soil is still moist from the winter and the ground water table is roughly at 1.2 m. Thus, we assume a Poisson's ratio of 0.45. For the sake of simplicity the Poisson's ratio is supposed to be depth-independent. This model is far from accurate, but it is acceptable for our purpose. We do not want to derive the absolute true values of the dispersion curves. The aim is to make statements about the relative effect of a thin frozen layer on Love and Rayleigh waves. Table 1 shows the parameters of the 1D model. For the following remarks we call this model the reference model, since it serves as a reference later to the model including the thin layer of ice, which is from now on called test model. For the test model we insert a thin frozen layer into the first layer of the reference model. From the temperature sensors we know that the thickness must be less than 5 cm. Hence, we estimate the shear wave velocity of the frozen layer as the fifth fold of the shear wave velocity in the initial state.

The P-wave velocity of the frozen layer is determined through the Poisson's ratio. After Zimmerman and King (1986) the Poisson's ratio decreases with the ice-to-water-ratio. They took permafrost samples made out of clay and silt and determined a Poisson's ratio of roughly 0.33. Since ice has a lower density than water, we use a density of 1.35 g cm -3 , which is 0.05 g cm -3 less dense than the first layer of our three layer over a halfspace model. The frozen layer is treated as part of the first layer, i.e. the thickness of the frozen layer is subtracted from the thickness of the first layer. The effect of the frozen layer is quite surprising, if we consider the wavelength of the surface waves. For 1 Hz the wavelength equals 350 m and no velocity variation appears. the wavelength of Rayleigh waves is around 300 m and still no velocity variation appears.

For Love waves, however, the wavelength is around 60 m and the velocity increase is more than 2 %. At 7 Hz both surface wave types have a wavelength of ca. 18 m. The velocity of Love waves increases by more than 5 % and the velocity of Rayleigh waves increases by 3 %. The thickness of the frozen soil is two to three magnitudes smaller than the wavelength of the surface waves, and nevertheless it affects spectacularly the velocity. In the following remarks, we examine the influence of some of the estimated parameters more profoundly.

Variation of the Poisson's ratio

In our initial 1D model we assumed a Poisson's ratio of 0.45 to estimate the P-wave velocity.

Now we investigate how the synthetic dv/v curves change with different Poisson's ratio. The shear wave velocities are fixed and, thus, the Poisson's ratio determines the P-wave velocity.

Since Love waves are completely independent of the P-wave velocity, it does not effect Love waves. For Rayleigh waves it is different. After [START_REF] Foti | Surface wave methods for near-surface site characterization[END_REF] the ratio between the Rayleigh wave speed and shear wave velocity of the medium can be expressed as a function of the Poisson's ratio ν. For ν close to 0.5 the ratio is closer to 1. Poisson's ratio. The velocity increase due to the frozen layer is the same for all Poisson's ratio. The velocity variation sensed by the Rayleigh waves do depend on ν. With a lower

Poisson's ratio the frozen layer causes a higher velocity increase for Rayleigh waves, especially for higher frequency. The sensitivity of Rayleigh waves towards the frozen layer seem to increase with lower ν. This implicates, that the sensitivity gap between Rayleigh and Love waves closes with lower ν, which is representative for a "harder" subsurface. Rayleigh waves in the frequency range of 3 and 6 Hz grew even more.

Variation of the test model

The four settings showed that an increase in the thickness or a decrease of the temperature causes higher velocity increases for both surface wave types. For the tested range of parameters, Love waves shows greater velocity increases than Rayleigh waves.

Discussion & Conclusion

Three three-component seismometers were installed in Hamburg (Germany) and velocity variations in the near subsurface were monitored with the PII method. The spectrograms reveal diurnal and weekly variations in the local noise field above 0.5 Hz, indicating an urban origin. Those noise fields are generated by human activity such as traffic, industry, railways etc.

Noise cross correlation was performed on the data from the three components of the three The parameter h defines the thickness of the frozen layer. v f rozen is the shear wave velocity of the frozen layer and v unf rozen is the velocity of the first layer (107 m s -1 ) .

A 1 cm thick layer with an average temperature of -5 • C increases the velocity of Rayleigh and Love waves by more than 1 % for frequencies above 6 Hz. The simulations also confirm that Love waves show larger velocity increases than Rayleigh waves.

At this point, we want to emphasize that we do not try to fit the observed data exactly with the model. The model had two purposes. Firstly, we wanted to know if surface waves are affected by frozen layers with a thickness, which is much smaller than the wavelength of the surface waves. Secondly, we wanted to prove the hypothesis that Love waves are more sensitive to this kind of velocity perturbation than Rayleigh waves are. Both assumptions are confirmed by the simulations.

This phenomena is not completely unknown and was already discussed in theory in [START_REF] Postma | Wave propagation in a stratified medium[END_REF] and [START_REF] Backus | Long-wave elastic anisotropy produced by horizontal layering[END_REF]. A horizontal layering with high velocity contrasts and layer thicknesses, which are much smaller than the observed wavelength, causes radial anisotropy.

This effect causes the well studied shear wave anomalies beneath the pacific plate [START_REF] Ekström | The unique anisotropy of the pacific upper mantle[END_REF]. In [START_REF] Jaxybulatov | A large magmatic sill complex beneath the toba caldera[END_REF] horizontally oriented sills underneath a caldera causes stronger velocity anomalies for Love waves than for Rayleigh waves. Our study shows a similar case on a much smaller scale with only one high velocity layer on top of the surface.

There are still open questions, which would be interesting to address in future work. First of all, the stretching method is a method to retrieve velocity changes from the coda and not from direct arrivals of surface waves. However, the results presented here are still valid to the first order. More novel techniques like the wavelet method [START_REF] Mao | On the measurement of seismic traveltime changes in the time-frequency domain with wavelet cross-spectrum analysis[END_REF] In general, the study raises the question whether other small scale processes can be registered and monitored by the surface wave response retrieved from PII. Especially the Love wave response in the high-frequency domain has the potential to be very sensitive towards environmental processes at the shallow subsurface.

Achknowledgements

We want to thank the Meteorological Institute and the Institute of Soil Science at the University of Hamburg for providing the temperature data set and the soil samples, which we used in this publication. In particular we want to thank Ingo Lange and Sarah Wiesner. Hamburg. The raw seismic data is available under doi:10.5281/zenodo.3992631.

Figure 1 :

 1 Figure 1: Maps illustrating the positions of the measuring site. a) shows the position of Hamburg in Northern Europe. b) shows the greater area of Hamburg. The black shaded rectangle indicates the position of map c). c) shows the seismometers WM01, WM02 and WM03 as red triangles in their surroundings. Next to WM02, a measurement station measures soil temperature and water content. The maps were created with Cartopy (Met Office, 2015) and Stamen. Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

  , T-T. Z denotes the vertical, R the radial and T the transversal component. The set of correlations of the different components is commonly named the Noise Correlation Tensor (NCT). If the noise sources are evenly distributed and produce surface waves, the NCT converges towards the Green's Tensor (GT) for surface waves. Campillo and Paul (2003) showed that the direct arrivals of the surface waves are separated into the Love wave arrivals on T-T and Rayleigh wave arrivals on Z-Z, Z-R, R-Z and R-R, when a rotation into the ZRT-system is performed . Moreover, there is no direct arrival of surface waves on Z-T, R-T, T-Z and T-R because Rayleigh and Love waves do not correlate with

Figure 2 :

 2 Figure 2: Spectrogram over two weeks starting at 8th of January 2018 for two components of WM01 with logarithmic frequency scale

  and data-quality of the NCFs depending on their frequency content: waveform coherency (WFC) and Signal-to-Noise-Ratio (SNR). WFC measures the temporal stability of the NCFs and indicates the arrival time of coherent signals, which are necessary for the stretching method. The WFC measures the coherence of the waveform with respect to the time lag and frequency content of the NCF. The SNR measures the energy of the signal compared to the noise and delivers an estimation how distinguishable the signal is. Both tools are applied on the NCFs.

Figure 3 :

 3 Figure 3: WFC depending on time lag and frequency averaged over all channel pairs for each station pair

  depicts the time series of the NCFs of all nine channel pairs of station pair WM01-WM02. The NCFs are filtered between 1 and 6 Hz. The stacking time is set to 4 days with an overlap of 2 days resulting in a time resolution of 2 days. The color coding illustrates the relative amplitude and the black function in each plot represents the RCF of each channel pair. All components show a coherent signal between -5 and 5 s time lag.

Figure 5 :

 5 Figure 5: The WFC (blue) and SNR (red) averaged over all channel-and station-pairs in the frequency range of 1 to 6 Hz. The investigated time period goes from Monday, the 8th of January 2018, to Sunday, the 14th of January 2018.

Figure 6 :

 6 Figure 6: The NCT for the station pair WM01-WM02 in the frequency range of 1 to 6 Hz. The color coding indicates the amplitude and the black graph represents the RCF.

Figure 8

 8 Figure 8 shows the velocity variation averaged over all station and channel pairs together with the temperature at the surface, at 5 cm and 80 cm depth, and with the cumulative water column (CWC) of the first 120 cm. The environmental data sets are originally recorded with a sampling rate of 10 min, however, here we show the daily averages due to reasons of clarity. Right before the velocity peaks in February and March the temperatures reach very low values. On the 1st of March the temperature at the surface reaches values below -10 • C.

Figure 7 :

 7 Figure7: dv/v time-series of all nine channel pairs averaged over causal and acausal side and all station pairs. The grey shaded area represents the error after equation 5.

Figure 8 :

 8 Figure 8: Averaged dv/v time series together with environmental parameters. Top: averaged dv/v time series with error. Middle: temperature at the surface, 5 and 80 cm depth. Bottom: cumulative water content of the first 120 cm.

Figure 9 :

 9 Figure 9: dv/v time series averaged over all channel pairs together with temperature data. The black line indicates the dv/v time series, the grey area shows the error after equation 5, the red line represents the temperature data at the surface and the blue shows the temperature data at 5 cm depth. The red dashed line marks 0 • C as the freezing point for water.

  Figure 2 from Miao et al. (2019) is used to estimate the shear wave velocity inside a frozen layer of clay, silt or sand. They measured the shear wave velocity inside these three soil types under a low confining pressure of 20 kPa and different moisture and temperature settings. λ(θ v , T )is the ratio between the shear wave velocity at a given temperature T and volumetric water content θ v , and the shear wave velocity at a temperature of 20 • C. The top of the soil at the measuring site in Hamburg is a mix of clay and silt. In our case the temperature is -10.3 • C at the surface and 0.2 • C at 5 cm depth. For the sake of simplification we assume an average temperature of -5 • C. After Figure2from[START_REF] Miao | Influence of seasonal frozen soil on near-surface shear wave velocity in eastern Hokkaido, Japan[END_REF] an average temperature of -5 • C inside a silt-clay mix results in an increase of 400 to 600 % of the shear wave velocity.

Figure 10

 10 Figure10depicts the results of the simulation. We calculated the dispersion curve for the fundamental Rayleigh and Love wave for the reference and test model. The thickness of the frozen layer is 4 cm and shear wave velocity is 535 m s -1 . The dispersion curves are shown at the top two plots. The incorporation of the thin layer of frozen soil alters the dispersion curves for Rayleigh and Love waves slightly. At the bottom we see the ratio between the dispersion curves of the test and reference model for both surface wave types in red. The ratio makes the effect of the frozen layer more visible. Both surface waves experience a frequency-dependent velocity increase. Above 1 Hz Love waves experience a greater velocity increase than Rayleigh waves. For higher frequency Love waves increase their velocity by more than 5 % while Rayleigh waves increase their velocity by less than 4 %. Below 1 Hz both wave types do not experience any velocity perturbation. This simple example already shows that both surface wave types are sensitive to a 4 cm thick layer of frozen soil. It also shows that Love waves with frequencies above 1 Hz are more sensitive to the velocity increase than Rayleigh waves are.

Figure 11

 11 Figure11shows the synthetic dv/v curves for Love waves on the left side and Rayleigh waves on the right side. As Poisson's ratios for the reference model we used 0.35, 0.45 and 0.49. The parameters for the thin layer of frozen soil are the same as in the previous example. As expected, the velocity variation for Love waves are not dependent on the

Figure 11 :

 11 Figure 11: Synthetic dv/v curves for different Poisson's ratio for Love waves on the left side and Rayleigh waves on the right side. Since there is no change for Love waves, the dispersion curves of the Love waves are located at the same place and the different graphs are not distinguishable.

  After considering different Poisson's ratio, we now investigate variations in the test model. More precisely, we look at the dv/v perturbations dependent on the thickness and temperature of the frozen layer. The last test models used a thickness of 4 cm and an average temperature of -5 • C. The next simulations use a thickness of 1 cm and a temperature of -2 • C additionally. According to Figure 2 from Miao et al. (2019), -2 • C inside a silt-clay layer results in a shear wave velocity increase of 300 % compared to a silt-clay layer at room temperature. The two different settings of temperature and thickness result in four different models. The Poisson's ratio of the frozen layer stays at 0.33 while the Poisson's ratio of the unfrozen part stays at 0.45.

Figure 12

 12 Figure12shows the velocity variation due to the frozen layer for Love and Rayleigh waves in the four different models. The top left shows the results for a thickness of 1 cm and a temperature of -2 • C. A velocity increase is noticeable, but it remains below 1 % for both wave types. Moreover Love waves show greater dv/v values than Rayleigh waves. The top right illustrates the results for a thickness of 4 cm and a temperature of -2 • C. Compared to the example from before, the increase of the thickness causes higher dv/v values, especially for frequencies above 1 Hz. Moreover the gap between Love and Rayleigh waves seems to widen, especially in the frequency range of 3 to 6 Hz. The lower left shows the results for a thickness of 1 cm and a temperature of -5 • C. The dv/v curves look similar to the last example. The lower right shows the results for a thickness of 4 cm and a temperature of -5 • C. This example uses the same test model as in Figure10. The gap between Love and

Figure 12 :

 12 Figure 12: The results of four simulations with four different test models in the four subplots. The parameter h defines the thickness of the frozen layer. v f rozen is the shear wave velocity of the frozen layer and v unf rozen is the velocity of the first layer (107 m s -1 ) .

  could deliver more accurate results. Secondly, in our observation we can not claim that we have a clean Love wave response on the T-T component or a clean Rayleigh wave response on R-R. It would be interesting to study this effect with circumstances, which allow to reconstruct cleaner surface wave responses on the horizontal components. The findings of this work are particularly interesting for permafrost monitoring in moun-tainous or polar regions. James et al. (2017, 2019) already used noise-based monitoring to track the freezing and thawing cycle of permafrost soils. Whereas they invert for the spatial distribution of the seismic velocity change, we show the ground truth of the freezing depth. The utilization of three-component seismometers can deliver additional information through the reconstruction of Rayleigh and Love waves. With the velocity information of surface waves, it is possible to track the penetration depth of freezing and thawing at centimeter scale. Through dispersion curve modelling one could invert for the thickness and temperature of the layer of interest. Here Love waves have an advantage, since only the density and shear wave velocity need to be known. Moreover, we show that seasonal freezing can produce strong near-surface velocity perturbations in areas with a moderate climate like in Hamburg at frequencies around 1 -6 Hz. Only a few centimeter thick frozen layer affects the velocity of surface waves by several percent. A few days of temperatures below 0 • C are already enough to freeze the first centimeters of the subsurface. This is important to know when applying PII to monitor other targets. If the target of interest is a geological feature, e.g. a fault or a storage site, seasonal freezing can mask other signals and be misinterpreted.

  Moreover, we want to mention Joachim Buelow, who set up the seismometers and maintains them. Prof. Dr. Christian Hübscher, Charlotte Bruland, David Essing, Lorenz Marten and Anna Hofmann were helping us with the refraction profile. We are also very grateful for the reviewers and editors who helped to improve the manuscript significantly. The work presented here is partly funded through the German Research Foundation (DFG) under Germany's Excellence Strategy -EXC 2037 'CLICCS -Climate, Climatic Change, and Society' -Project Number: 390683824, and the Cluster of Excellence 'CliSAP' (EXC177), contribution to the Center for Earth System Research and Sustainability (CEN) of Universität

Table 1 :

 1 1D model of the subsurface at the measuring site. S-wave velocity v s and thickness are inferred from an SH-refraction profile. P-wave velocity v p and densities are estimated with information from soil samples.

	For 3 Hz

seismometers. The data set ranges from the 5th of January to the 14th of April of the year 2018. The data was sliced into 10 minute windows, which was then cross correlated. The analysis of SNR and WFC suggests that the noise sources are inhomogeneously distributed in space. We also observe high SNR on T-T between 1 and 6 Hz indicating a retrieval of Love waves with higher quality. This may be due to a high ratio of Love waves in the local noise field.

A stack of 576 NCFs, equalling four days of data, produced stable NCFs, which were used for the stretching method. dv/v time-series in the frequency band 1 -6 Hz were created with all components. All time-series show similar results and trends with variation in the amplitude.

The velocity variations correlate well with the temperature data. The freezing and thawing process is clearly visible in the dv/v data. Temperature sensors inside the ground indicate that the penetration depth of the frozen soil is not more than 5 cm. For some channel pairs the velocity increases by more than 2 % due to the freezing. The observations lead to the hypothesis that the type of surface wave inherent in the cross correlations might play a major role.

In order to prove the influence of the centimeter thick layer of frozen soil, dispersion curves of Rayleigh and Love waves have been modelled. Firstly, a 1D model from an SH-refractionprofile was derived (reference model). At the surface, a frozen layer was inserted (test model).

The dispersion curves were calculated for the test and reference model. The ratio of the dispersion curves reveals the relative velocity increase due to the inserted frozen layer. The simulation confirms that a few centimeter thick layer of frozen soil increases the velocity of Rayleigh and Love waves by several percent, even if their wavelength is two to three magnitudes larger than the thickness of the frozen layer. The Poisson's ratio influences the effect of the frozen layer for Rayleigh waves but not for Love waves. A lower Poisson's ratio raises the velocity increase due to the frozen layer for Rayleigh waves. Thickness and temperature of the frozen layer influence the impact on the relative velocity change for both wave types.