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

S U M M A R Y 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 1-D medium we can verify that surfaces waves of several tens of metres wavelength experience a velocity increase of several percent due to a centimetres 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.


I N T RO D U C T I O N
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 wavefield 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 GF of the medium (Lobkis & Weaver 2001;Campillo & Paul 2003). 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 Poupinet et al. (1984). Sens-Schönfelder & Wegler (2006) 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 wavefields and natural seismic sources to study the coupling between the Earth and the external environment (Larose et al. 2015). Erosion, thermal forcing, precipitation, freezing etc. affect the surface of the Earth and, thus, modify the propagation of a wavefield Effect of centimetric freezing in original form date 627 at the near subsurface. In the last years, many studies used ambient seismic noise to monitor different environmental effects on seismic velocities. Richter et al. (2014) could track the daily and seasonal temperature variations inside the first few metres of the subsurface with PII. James et al. (2017James et al. ( , 2019 used ambient seismic noise to monitor the freezing and thawing cycle of permafrost soils in the Arctic. Prior to a landslide event, Mainsant et al. (2012) observed a decrease of seismic velocity at the base of the sliding layer with PII. Lecocq et al. (2017) used 30 yr 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, for example, temperature, water content or changes of the noise sources (e.g. Tsai 2011;Zhan et al. 2013;Lecocq et al. 2017). Gassenmeier et al. (2014) 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 Fig. 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.

Passive seismic interferometry
The method of PSI is briefly introduced. For a more detailed explanation we refer to Wapenaar et al. (2010). PSI turns passive seismic measurements, for example, 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: 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 GF 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, Lin et al. (2009) retrieved the vertical Rayleigh wave response of a virtual source with uneven distributed noise sources and Hadziioannou et al. (2009) 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: Z-Z, Z-R, Z-T, R-Z, R-R, R-T, T-Z, T-R and T-T. Z denotes the vertical, R the radial and T the transversal components. 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 & 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 ZRTsystem 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 each other (Roux 2009). 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.
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 (Roux 2009). 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 source-receiver 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. Poupinet et al. (1984) 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. Poupinet et al. (1984) interpreted the retrieved velocity variations as a consequence of variations in the subsurface shear wave velocity. Snieder et al. (2002) 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: This implies that one can retrieve dv/v by simply measuring the relative time-shift dt/t. Snieder et al. (2002) called this concept Coda-Wave-Interferometry (CWI). The stretching method, which is used in this study, is one of several methods to estimate the time-shift (Mikesell et al. 2015). 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 ): 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 everyday 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 eq. (3) we have an RCF serving as a reference to every NCF. The grid search over i would then look like: (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. Hadziioannou et al. (2009) 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 crosscorrelation 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 Downloaded from https://academic.oup.com/gji/article/224/1/626/5898674 by UB Frankfurt/Main user on 12 March 2021 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. Bussat (2015) presents an example, where a velocity change of 2 % within 4 hr 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 an NCF and, thus, be misinterpreted as physical changes in the subsurface. Weaver et al. (2011) developed an expression to estimate this kind of error: 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 eq. (5) to verify the apparent velocity changes we retrieve in this study.

DATA A C Q U I S I T I O N A N D P R E P RO C E S S I N G
The data used in this work, were recorded by three three-component Trillium Compact broad-band 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 2018. The data are recorded in UTCtime 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. First, the data were downsampled to 100 Hz to save computation time for all the following work. The aliasing effect was avoided by a low-pass 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 to 100 Hz. Most processing was carried out with the Python package Obspy (Beyreuther et al. 2010).

L O C A L N O I S E F I E L D
In order to characterize the local noise field, a spectrogram of a week in 2018 February 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 during the day in the week and weakening during the night and at the weekends. Frequencies between 3 and 5 Hz seem to be more stable and do not vary as much during the weekend and nights. This is favourable for creating stable NCFs. From the map in Fig. 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 (Bensen et al. 2007).

P R E L I M I N A RY I N V E S T I G AT I O N S
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 and data quality of the NCFs depending on their frequency content: waveform coherence (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. Before performing the cross-correlation, the preprocessed data were processed another time. First the data were rotated from the vertical-north-east-system ZNE into the vertical-radial-transversalsystem 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 and 20 Hz was applied. To exclude earthquakes and other transient events, amplitudes higher than the sixthfold of the standard deviation have been clipped. Then, the slices were spectral whitened between 0.3 and 20 Hz. After the processing was done, a cross-correlation in the time domain with a shift of 6000 samples-which equal 60 s 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 coherence
How is the WFC calculated? First, 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 broad-band NCF contains frequencies between 0.3 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: 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 −25 to 25 s with an interval of 0.5 s, because with the greatest interstation distance of 153 m we do not expect coherent signals after 25 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. A value of 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. Fig. 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 centred asymmetrically between −5 and 5 s. The asymmetry suggests an uneven distribution of noise sources around the stations. Since the station pairs have a similar interstation 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 incoherent waveforms start to dominate, and thus CC decreases.

Signal-to-noise ratio
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: 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 sides. 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. Fig. 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 (Zhan et al. 2013;Bussat 2015). 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 an 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. Fig. 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 d, 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 Fig. 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 d, 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. Fig. 6 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 d with an overlap of 2 d resulting in a time resolution of 2 d. The colour 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. The NCFs are asymmetric due to the uneven distribution of noise sources. The strongest amplitudes can be found for T-T. This agrees with the SNR analysis from Fig. 4. The NCFs time-series shown here is now used to retrieve relative velocity variations with the stretching method.

R E L AT I V E V E L O C I T Y C H A N G E D E R I V E D F RO M P I I A N D E N V I RO N M E N TA L PA R A M E T E R S
For the stretching method, we define a time window on the causal and acausal sides, respectively. Since we want to focus on surface waves, we place the time windows at the first seconds. For the causal side, the time window ranges from 0.1 to 4 s, and for the acausal side, the time window ranges from −4 to −0.1 s. The retrieved velocity variations for all nine channel pairs averaged over all station pairs are displayed in Fig. 7. The grey shaded area represents the error calculated using eq. (5). Two major features, visible on almost all components, are the two velocity increases at mid-February and the beginning of March. The two increases are less prominent on Z-Z, Z-R and T-Z. For all other components the two increases are around 2 %. The asymmetry is quite striking. It is again an indication for an uneven distribution of noise sources with different characteristics around the stations. Fig. 8 shows the velocity variation averaged over all station and channel pairs together with the temperature at the surface, at 5 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. 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 (Zimmerman & King 1986). Gassenmeier et al. (2014) 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, Fig. 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 d. 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 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 eq. (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. velocity increase of roughly 0.5 % happens between mid-March and beginning of April. Here we have five 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.

Setup of the model and first results
The simulation of the dispersion curves is done in Geopsy (Wathelet 2006). For the subsurface, we assume a 1-D 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 1-D model. For the Downloaded from https://academic.oup.com/gji/article/224/1/626/5898674 by UB Frankfurt/Main user on 12 March 2021 Table 1. 1-D 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.
Density (g cm −3 ) 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. Fig. 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 fig. 2 from Miao et al. (2019), an average temperature of −5 • C inside a silt-clay mix results in an increase of 400-600 per cent of the shear wave velocity. 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 & 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 half-space model. The frozen layer is treated as part of the first layer, that is, the thickness of the frozen layer is subtracted from the thickness of the first layer. Fig. 10 depicts the results of the simulation. We calculated the dispersion curve for the fundamental Rayleigh and Love waves 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.
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. For 3 Hz, 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 2-3 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 1-D 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 Foti et al. (2014), 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. Fig. 11 shows the synthetic dv/v curves for Love waves on the left-hand side and Rayleigh waves on the right-hand 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 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.

Variation of the test model
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 fig. 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. Fig. 12 shows 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-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 Fig. 10. The gap between Love and Rayleigh waves in the frequency range of 3 and 6 Hz grew even more.
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.

D I S C U S S I O N A N D C O N C L U S I O N
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 seismometers. The data set ranges from 5th of January to the 14th of April 2018. The data were sliced into 10 min 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. First, a 1-D model from an SH-refraction profile 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 2-3 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. A 1 cm thick layer with an average temperature of −5 • C increases the velocity of Rayleigh and Love waves by more than 1 per cent 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. First, 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 Postma (1955) and Backus (1962). 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 (Ekström & Dziewonski 1998). In Jaxybulatov et al. (2014), horizontally oriented sills underneath a caldera cause 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 (Mao et al. 2020) could deliver more accurate results. Secondly, in our observation we cannot 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 mountainous or polar regions. James et al. (2017James et al. ( , 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, for example, a fault or a storage site, seasonal freezing can mask other signals and be misinterpreted.
In general, the study raises the question whether other smallscale 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.

A C K N O W L E D G E M E N T S
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. Moreover, we want to mention Joachim Buelow, who set up the seismometers and maintains them. 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.