ISSN 1794-6190 e-ISSN 2339-3459 https://doi.org/10.15446/esrj.v29n3.119268 EARTH SCIENCES RESEARCH JOURNAL Earth Sci. Res. J. Vol. 29, No. 3 (September, 2025): 333 - 345 G EO PH Y SI C S Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area, North China using ambient seismic noise data Yipei Tan*, Ting Ma, Li Deng, Ke Xu Tianjin Earthquake Agency, Tianjin, China, 300210, China. *Corresponding author: oivertan921@sina.cn Record Manuscript received: 08/03/2025 Accepted for publication: 07/09/2025 How to cite this item: Tan, Y., Ma, T., Deng, L., & Xu, K. (2025). Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area, North China using ambient seismic noise data. Earth Sciences Research Journal, 2(3), 333-345. https://doi.org/10.15446/esrj.v29n3.119268 ABSTRACT It is feasible to continuously monitor the temporal changes of shallow subsurface media by utilizing ambient seismic noise. Here, the authors analyzed three-year continuous records of 117 earthquake early warning stations in the Tianjin area in the North China Plain to explore the seasonal variation of shallow subsurface seismic velocity. They examined the daily horizontal-to-vertical spectral ratio (HVSR) curves and found that 25 stations exhibit clear seasonal variations in peak frequency and peak magnitude within the frequency range of 5-20 Hz. Subsequently, autocorrelation is applied to the seismic noise of these 25 stations in the 5-20 Hz frequency range to calculate the relative velocity changes. Ten stations show the same variation pattern in the peak frequency of HVSR curves and relative velocity changes based on autocorrelation. This pattern exhibits characteristics where the peak frequency and seismic velocity increase rapidly in November, decrease quickly in February of the following year, and change relatively little in other seasons. By compa- ring with meteorological data, the authors infer that the freezing of the soil layer when the temperature drops below 0°C leads to an increase in shear wave velocity and peak frequency in winter. The results suggest that studies on crustal structure detection need to consider temperature effects. Keywords: Seismic noise; Velocity variation; HVSR method; Seismic interferometry; Seasonally frozen soil Características de las variaciones estacionales de la velocidad sísmica del subsuelo superficial y sus factores influyentes en el área de Tianjin, norte de China, utilizando datos de ruido sísmico ambiental RESUMEN Es posible monitorear continuamente los cambios temporales de los medios del subsuelo superficial mediante el uso del ruido sísmico ambiental. En este estudio, los autores analizaron registros continuos de tres años de 117 estaciones de alerta temprana de terremotos en el área de Tianjin, en la llanura del norte de China, para explorar la variación estacional de la velocidad sísmica del subsuelo superficial. Examinaron las curvas diarias de la relación espectral ho- rizontal-vertical (HVSR) y descubrieron que 25 estaciones presentan claras variaciones estacionales en la frecuencia y magnitud pico dentro del rango de frecuencia de 5 a 20 Hz. Posteriormente, se aplicó la autocorrelación al ruido sísmico de estas 25 estaciones en el rango de frecuencia de 5 a 20 Hz para calcular los cambios de velocidad relativa. Diez estaciones muestran el mismo patrón de variación en la frecuencia pico de las curvas HVSR y los cambios de velocidad relativa basados en la autocorrelación. Este patrón presenta características donde la frecuencia pico y la velocidad sísmica aumentan rápidamente en noviembre, disminuyen rápidamente en febrero del año siguiente y varían relativamente poco en otras estaciones. Al comparar datos meteorológicos, los autores infieren que la congelación de la capa de suelo cuando la temperatura desciende por debajo de 0 °C provoca un aumento de la velocidad y la frecuencia pico de las ondas de corte en invierno. Los resultados sugieren que los estudios sobre la detección de la estructura de la corteza deben considerar los efectos de la temperatura. Palabras clave: Ruido sísmico; Variación de la velocidad; Método HVSR; Interferometría sísmica; Suelo congelado estacionalmente; https://doi.org/10.15446/esrj.v29n3.119268 mailto:oivertan921@sina.cn https://doi.org/10.15446/esrj.v29n3.119268 334 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area 1 Introduction Shallow subsurface seismic velocity is crucial in the fields of seismology and earthquake engineering, as it not only serves as a foundation for accurately inverting earthquake source parameters and understanding deeper structures (Dong et al., 2018; Langston, 2011), but also plays a vital role in seismic site response and earthquake loss assessment (Choi & Stewart, 2005; Bendito et al., 2014). Recent studies have revealed that subsurface seismic velocity is unstable, but exhibits rapid changes during strong ground vibrations (Bonilla et al., 2019; Vassallo et al., 2022) and seasonal variations influenced by environmental factors such as temperature (Kula et al., 2018; Luan et al., 2022), rainfall (Miao et al., 2018; Roumelioti et al., 2020), atmospheric pressure (Gradon et al., 2021; Kramer et al., 2023), snow depth (Wang et al., 2017), and groundwater levels (Rigo et al., 2021; Mao et al., 2022; Qin et al., 2022; Delouche & Stehly, 2023). Continuous seismic monitoring of shallow subsurface seismic velocity can elucidate its temporal variation characteristics and the factors influencing these variations, thereby enhancing the accuracy and reliability of seismological and earthquake engineering research outcomes. Monitoring shallow subsurface seismic velocity can rely on vibration signals from natural earthquakes or active seismic experiments (Niu et al., 2008; Bonilla & Ben-Zion, 2021; Luan et al., 2022). However, the sporadic nature of natural earthquakes makes it challenging to obtain continuous temporal variation data. Artificial seismic sources are costly, and their excitation conditions are often constrained by human activities, particularly in urban settings. Consequently, methods that monitor seismic velocity variations using seismic ambient noise have gained widespread adoption. The Horizontal-to- Vertical Spectral Ratio (HVSR) method, introduced by Nakamura (1989), is a straightforward and robust technique widely utilized for mapping soft sediment thickness and studying site effects (Ibs-von Seht & Wohlenberg, 1999; Qin et al., 2022; Chen et al., 2023). The HVSR curve is derived from the ratio of the combined frequency spectra of horizontal components to the vertical component of seismic noise. The peak frequency of the HVSR curve correlates with the depth of geological interfaces and the average shear wave velocity between the surface and the interface (Nakamura, 2019). Recently, the HVSR method has been applied to study shallow subsurface seismic velocity changes across Europe (Lontsi et al., 2022; Seivane et al., 2022), Asia (Chen et al., 2023; Liu et al., 2024), and Africa (Benkaci et al., 2018). Another promising method involves seismic noise correlations, which involves repeatedly extracting the Green’s function through the calculation of autocorrelation or cross- correlation of ambient noise records and measuring minor time-shifts in the Green’s function to detect seismic velocity changes. By filtering noise records in different frequency bands and/or analyzing various time windows, this correlation method can monitor variations from shallow subsurface changes caused by environmental factors (James et al., 2017; Feng et al., 2021; Oakley et al., 2021; Zhang et al., 2023) to deeper earth changes resulting from volcanic eruptions (Liu et al., 2022), earthquakes (Poli et al., 2020; Okubo et al., 2024), and deep tectonic activities (Li et al., 2021; Ashruf & Morelli, 2022). Some researchers (Tchawe et al., 2020; Li & Ben-Zion, 2023) have combined HVSR and autocorrelation to enhance the precision and stability of their measurements. In this study, we aim to investigate the characteristics of shallow subsurface seismic velocity variation and the potential factors influencing these characteristics in Tianjin area, utilizing HVSR and autocorrelation method based on data from earthquake early warning stations established in June 2021. Tianjin, situated in the North China Plain (Figure 1), is predominantly covered by Quaternary strata, with bedrock exposure limited to the northern mountainous region (Tianjin Institute of Geological Survey, 2018). The region has experienced numerous strong earthquakes (M≥7), causing significant damage to major cities in the plain. The findings of this study will provide valuable insights for ground motion simulation, focal parameter inversion, and the detection of temporal variations in deep earth media in the area. 2 Data and methods 2.1 Seismological and meteorological data Within the study area, 117 earthquake early warning stations (refer to Figure 1) were established to facilitate rapid early warning and intensity reporting following significant earthquakes. These stations record the physical quantity of acceleration continuously for 24 hours at a sampling rate of 100 points per second. These stationsstarted operation in June 2021 and continued to operation for three years until May 2024.Among these 117 stations, two are situated in the northern exposed bedrock area, where there is no overlay of soft soil. The remaining 115 stations are founded on sedimentary soil. This paper will utilize the three-component seismic waveforms recorded by these 115 stations over the past three years. Figure 1. Distribution of seismic and meteorological stations in study area. Different colors indicate the age of the surface covering soil layer or rock exposure. NCP signify North china Plain. We have access to the recorded data from 12 meteorological observation stations within the study area, with their site distribution depicted in Figure 1. Each of these stations records temperature, air pressure, and precipitation continuously for 24 hours at a sampling rate of 1 point per minute. The meteorological data will be employed for comparative analysis with the observed variations in seismic velocity. 2.2 HVSR We employ the following formula to calculate the HVSR curve daily at each station: HVSR( ) = ((E( )2+(N( )2)/2 Z( ) (1) Where f denotes the frequency, and Z(f), E(f) and N(f)represent the Fourier spectra of the vertical, east-west and north-south component seismic records, respectively. The continuous waveform data for a 24-hour period is initially processed by removing the mean and trend. We check for any broken or missing data points and mark their locations. Subsequently, we calculate the median absolute deviation (MAD) for the 24-hour waveform and flag data points exceeding 5 times the MAD, along with the 50 data points preceding and following them, which may be affected by human activities near the station, causing significant jumps. We opt not to use the STA/LTA algorithm for earthquake detection for two reasons. First, some studies indicate that the HVSR curve derived from ambient noise is similar to that calculated from seismic records (Nakamura, 2019; Yassminh et al., 2019), rendering it unnecessary to filter out earthquakes when computing the HVSR curve (Rigo et al., 2021). Second, over the three- year period, there were only 31 earthquakes with a magnitude above 1.0 in the study area, averaging less than one per month. The low seismic activity means that seismic records have minimal impact on the HVSR curve. The continuous 24-hour recorded waveform is segmented into 40.96-second lengths (4096 sampling points) with a 50% overlap. With 8,640,000 sampling points per day, this results in 4,218 segments. Any segment marked with issues from the previous steps is discarded. If fewer than 20 335Yipei Tan, Ting Ma, Li Deng, Ke Xu segments remain for the day, the HVSR curve for that day is not computed. We then apply a cosine taper to each data segment and use the Fast Fourier Transform (FFT) to calculate the two horizontal and one vertical spectrums. The average of the horizontal and vertical frequency spectra for all segments over the 24-hour period is calculated separately, and the division of the horizontal by the vertical average frequency spectra yields the 1-day HVSR curve. This curve is smoothed using the Konno-Ohmachi logarithmic window (Konno & Ohmachi, 1998) with a b-value of 30 in the frequency range of 0.05-40Hz, resulting in the 1-day HVSR curve. La Rocca et al. (2020), La Rocca and Chiappetta (2022) suggest that the HVSR curve may vary due to daily temperature and noise source differences, with seismic noise being more stationary during night hours. We tested calculating HVSR curves using data from 0:00-6:00 and found no significant difference in temporal variation characteristics. Since the meteorological data for comparison with seismic velocity changes is the daily average temperature, we calculate the average HVSR curve using ambient noise from the full 24- hour period. Geological survey results (Tianjin Institute of Geological Survey, 2018) indicate that the thickness of Quaternary sediments in the study area ranges from about 40 meters in tectonically uplifted areas to 473 meters in tectonically depressed areas. The Holocene stratum, which is 17-24 meters thick and consists of loose soft soil deposits, is divided into three sections based on lithology. Based on the empirical relationship between the burial depth of Holocene sediments and the peak frequency of the HVSR curve in the North China Plain, h=94.76×fr -1.12 (Qin et al. 2023; Peng et al. 2019), the peak frequency corresponding to the bottom interface of Holocene sediments should be below 4.5Hz. Therefore, we select a frequency range greater than this value. To analyze the variation of seismic velocity in Holocene sediments, we focus on peak frequencies greater than 4.5Hz. Considering the instability of the HVSR curve at frequencies higher than 20Hz (James et al., 2017), we select stations with HVSR curves peaking in the 5-20Hz range. By extracting the peak frequency and amplitude daily within this range, we can track the changes in shallow subsurface seismic velocity. 2.3 Autocorrelation We estimate the Green’s function for each station using the software MSNoise (Lecocq et al., 2014). Our workflow begins by integrating the acceleration waveforms into velocity records, followed by the removal of the mean and trend, and then tamper with the 1-day waveform. Since the autocorrelation calculation is performed at a single station, the instrument response has not been removed. Subsequently, followed by one-bit normalization and spectrum whitening, we compute the daily autocorrelation of the 3-components within the 5-20Hz frequency band every 30 minutes and superimpose the 48 30-minute autocorrelation results to form the daily Green’s function. By linearly superimposing the daily Green’s functions over three years, we obtain the reference Green’s function. To improve the signal-to-noise ratio, we apply a 5-day stacking for each day, covering two days before and two days after, resulting in a 5-day-stacked Green’s function. We employ the stretching method (Sens-Schöfelder & Wegler, 2006) to estimate the changes in seismic velocity over time. For the i-th day, Gi(t) is extracted from the 5-day-stacked Green’s function after 100 times interpolation within the time window from -0.4s to 0.4s. The sampling point interval of Gi(t) is 0.0001s, with a total of 8001 sampling points. We resample the reference Green’s function at an interval of 0.0001ε, and take 8001 sampling points around the zero point to obtain Gstack(t). The cross-correlation cci(ε) of Gi(t) and Gstack(t) is then calculated. During the calculation process, ε ranges from 0.5 to 1.5. When the cross-correlation cci(ε) reaches its maximum, ε is the stretching coefficient of the i-th day, and 1-ε represents the relative change in seismic velocity for the i-th day compared to the reference Green’s function. The uncertainty δ is calculated using the following formula (Weaver et al., 2011): δ = 1−X2 2X 6 π 2T0 ω2(t2 3−t1 3) (2) Where X is the maximum cross-correlation, ω is the angular frequency, T0 is the central period of the reference Green’s function, and t1 and t2 are the start and end times of the measurement window, respectively. 3 Results 3.1 Seasonal variations of HVSR peak frequency We calculated the daily HVSR curves for 115 stations over a period of three years and then averaged these daily curves to obtain the 3-year-average HVSR curves for each station. Among these, two stations showed no significant peak values between 0.1 and 40 Hz in their 3-year-average HVSR curves (e.g., WAK as shown in Figure 2). Of the remaining 113 stations, 106 exhibited peaks in the 0.5-2 Hz range, and 81 displayed peaks within the 5-20 Hz range. Upon examining the resonant frequency of the 3-year-average HVSR curves, we found that 91 stations, accounting for 79.1% of the total, had resonant frequencies in the 0.5-2 Hz range (e.g., AX001 and KB001 in Figure 2). This result is largely consistent with Cao et al. (2018), who reported a resonant frequency of approximately 1 Hz for HVSR curves derived from tremor data in the northern part of our study area. Although 81 stations had peaks in the 5-20 Hz range, only six stations had resonant frequencies in this range (e.g., SJT in Figure 2), with most other stations showing resonant frequencies in the 0.5-2 Hz range and a smaller amplitude peak in the 5-20 Hz range (e.g., KB001 in Figure 2). Figure 2. 3-year-average HVSR cures of four stations represent four types WAK represents no significant peak in 0.1-40 Hz, AX001 represents the resonant frequency is about 1Hz and no significant peak in 5-20 Hz, KB001 represents resonant the frequency is about 1Hz and have significant peak in 5-20 Hz, SJT represents the frequency is in 5-20 Hz and have a peak about 1Hz. 336 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area Taking FJ001 station as a case study, we present the temporal changes in the HVSR curve (Figure 4). Initially, we arranged the daily HVSR curves and used color to represent the amplitude of the curves (Figure 3a). It was observed that from June 2021 to May 2024, the resonant frequency of FJ001 remained around 1 Hz, while its peak amplitude exhibited seasonal variations. The peak amplitude notably decreased at the end of one year and increased at the beginning of the next. Above 4 Hz, the overall amplitude of the HVSR curve decreased in winter. A small peak in the 5-20 Hz frequency band (pink box in Figures 3a and 3c) also showed a significant decrease in amplitude during winter. We extracted the peak frequency and amplitude from the daily HVSR curves within the 5-20 Hz range, and their temporal changes are depicted in Figure 3b, both of which displayed clear seasonal variation characteristics. From spring to autumn each year, the peak frequency was generally stable at around 10 Hz but rose rapidly to 14-16 Hz in November and dropped back to about 10 Hz by February of the following year. The peak amplitude varied inversely to the peak frequency, with small variations from spring to autumn but a rapid decrease in November and increase in February. Compared to summer, the peak frequency in winter increased by about 40% to 60%, while the peak amplitude in summer was reduced to about 40%. Upon reviewing the temporal variation results of HVSR curves at all 115 stations, we categorized them into four groups. (1) 24 stations could not calculate the HVSR curve for more than 200 consecutive days due to data scarcity, significant interference from surrounding human activities, substantial changes in the HVSR curve shape over three years (potentially due to station maintenance, instrument replacement, or environmental changes), or the absence of a clear peak in the 3-year-average HVSR curve. These stations were not suitable for analyzing seasonal variation characteristics. (2) The 3-year- average HVSR curves of 26 stations lacked a clear peak in the 5-20 Hz range, precluding analysis of seasonal variation in this range. (3) The 3-year-average HVSR curves of 40 stations had clear peaks in the 5-20 Hz range but showed no significant seasonal variation in peak frequency or amplitude. (4) The 3-year- average HVSR curves of 25 stations had clear peaks in the 5-20 Hz range, with significant seasonal variations in peak frequency or amplitude. Given the aim of this study is to investigate the characteristics of seasonal variations, we focused our analysis on the fourth category. The temporal variation results of peak frequency and amplitude for the 24 stations in the fourth category, excluding FJ001 (as shown in Figure 3), within the 5-20 Hz range are displayed in Figure 4. Despite minor differences in individual stations, such as the change range of KB001, ND001, and YGZ stations during the winter of 2023-2024 being smaller than in the previous two years, the overall variation characteristics of 21 stations were consistent with those of FJ001 station in Figure 3. In total, 22 out of the 25 stations exhibited similar variation characteristics: a rapid increase in peak frequency in November and a decrease in February of the following year, and a rapid decrease in peak amplitude in November followed by an increase in February of the following year,with the exception of CH001, GS001, and GN001. Figure 3. Temporal variation of H/V ratio at station FJ001 for the 3 years (a) Temporal variation of the HVSR curves for the 3 years, (b) Temporal variation of the peak frequency and peak amplitude in the frequency range of 5-20Hz, (c) 3-year-average HVSR curve. The pink box in (a) and (c) represent 5-20Hz. 337Yipei Tan, Ting Ma, Li Deng, Ke Xu Figure 4. Temporal variation of the peak frequency and peak amplitude in the frequency range of 5-20Hz of 25 stations whose 3-year-average HVSR curves have obvious peaks in 5-20 Hz. 338 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area 3.2 Velocity changes based on autocorrelation For the 25 stations exhibiting clear seasonal variations in their 3-year- average HVSR curves, we employed the autocorrelation method to estimate seismic velocity changes within the 5-20Hz frequency range.Figure 5 illustrates the 5-day-stacked and reference Green’s functions at station MG001. The reference Green’s function in the E component (Figure 5a) closely resembles the N component (Figure 5b), but differs from the vertical component (Figure 5c). We superimposed the 5-day-stacked Green’s functions for July, representing the highest temperatures (summer), and January, representing the lowest temperatures (winter), over the three-year period. In the two horizontal component results (Figures 5d and e), the summer Green’s function (red curve) is similar to the reference Green’s function (black curve), with corresponding peaks and troughs of several sidelobes appearing at similar times, and the amplitude of the summer Green’s function is slightly larger than that of the reference. However, the amplitude of the winter Green’s function (blue curve) is significantly smaller than that of the summer and reference Green’s functions. As the absolute value of the time offset increases, the difference in arrival times between the peaks and troughs of the winter and summer Green’s functions becomes more pronounced. In the vertical component results (Figure 5f), the main lobes of the summer, winter, and reference Green’s functions are almost completely coincident, with the sidelobes’ amplitude decaying rapidly. There is a significant difference in the horizontal component results between summer and winter peaks and troughs, but little difference in the vertical component results, a pattern common to other stations in our analysis. Therefore, we decided to use the horizontal component Green’s functions to calculate seismic velocity changes. Figure 6 presents the seismic velocity variation results at station MG001. The maximum correlation coefficient, obtained by cross-correlating the 5-day- stacked and reference Green’s functions (Figure 6a), decreases in winter but remains above 0.95 in other seasons. The stretch coefficient is determined when the correlation coefficient is at its maximum; if it exceeds 1, the 5-day-stacked Green’s function is compressed on the time axis, and if it is less than 1, it is stretched. The seismic velocity change results, based on the conversion of the stretch coefficient, are shown in Figure 6b, and Figure 6c displays the velocity change error calculated using formula (2). The 3-year velocity variation results for MG001 station indicate that seismic velocity increases rapidly at the end of each year and decreases rapidly at the beginning of the following year. The seismic velocity in winter is 25% to 35% higher than in other seasons. The velocity change error in winter is also larger, approximately 1% in winter and about 0.5% in other seasons, with individual days ranging from 2% to 3%. This error level is significantly smaller than the winter velocity changes, making the seasonal variation pattern of seismic velocity in Figure 6c reliable. Figure 5. The 5-day-stacked and reference Green’s function at the station MG001 (a)-(c) show the 5-day-stacked and 3-year-stacked Green’s function of E, N and Z component respectively, (d)-(f) show the comparison of Green’s functions stacked in July and stacked in January Figure 6. Results of the seismic velocity variation at station MG001 (a) represents changes of autocorrelation coefficient calculated by 5-days-stacked and reference Green’s function, (b) represents velocity changes (dv/v), (c) represents the velocity change error. 339Yipei Tan, Ting Ma, Li Deng, Ke Xu Figure 7 displays the velocity variation results for the 24 stations, excluding MG001 shown in Figure 6. The black dots represent relative velocity variation (dv/v), and the gray vertical lines indicate the error bars. Five stations (GN001, GS001, PM001, QD002, XD003) show no periodic pattern in their velocity variations. The velocity of the other 19 stations exhibits some degree of seasonal variation characteristics. For these 19 stations, we compared the peak frequency variation in the 5-20Hz range of HVSR curves (Figure 4) with the velocity variation based on autocorrelation (Figure 7). The velocity variation at CH001 station consistently decreases from winter to summer and increases from summer to winter, aligning with the variation pattern of the HVSR peak frequency, but it is not the pattern discussed in this paper. The wave velocity at QH001 station shows a pattern of decreasing in winter and being relatively high in other seasons, while PX002 station experiences a sudden rise in a few days during winter followed by a fall to the baseline level, inconsistent with the HVSR peak frequency variation characteristics. Although the seismic velocity at HBT, ND001, and NW001 stations increased in winter, the wave velocity fluctuated greatly, and the error bars indicate large errors in the winter calculation results. The wave velocity at HD001, NH001, and QM001 stations did not significantly increase in winter during 2021-2022, contrasting with the noticeable increase in HVSR peak frequency that year. The wave velocity at DUL station increased in winter but changed greatly from summer to autumn in 2023, with the fluctuation amplitude being essentially the same as that of the winter wave velocity change. By comparison, the seasonal variation pattern of the peak frequency for 10 stations is largely consistent with the velocity variation pattern (CAD, FJ001, JIH, KB001, MG001, PW001, QB001, QC001, WY001, YGZ), all of which show the characteristics of seismic velocity increasing rapidly in November, decreasing rapidly in February of the following year, and changing relatively little in other seasons.This consistent variation pattern indicates that the velocity variation results for these 10 stations are reliable and will be used in the subsequent comparison with meteorological data. Figure 7. Velocity variations at 24 stations except station MG001 showed in Figure 6 340 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area 4 Comparison with meteorological data 4.1 Temperature Correlation We compared the velocity variations of the 10 stations identified in the previous section with the daily average temperature recorded by the nearest meteorological observatory (Figure 8). The comparison results are depicted in Figure 9. The variation patterns of the peak frequency of HVSR and velocity based on autocorrelation are similar in the winter of three consecutive years. It is observable that the timing when the daily average temperature falls below 0°C aligns with the onset of rapid increases in the peak frequency of HVSR and seismic velocity in November. Some stations, such as CAD and PW001, exhibit a velocity increase 1 or 2 days prior to the temperature drop below 0°C, possibly due to the use of 5-day-stacked data in autocorrelation calculations. Conversely, the time when the peak frequency and velocity revert to summer values lags significantly behind the temperature rise above 0°C. This lag cannot be attributed to the 5-day-stacked data in autocorrelation, leading us to infer that the pronounced seasonal increase in seismic velocity during winter may correlate with the freezing of the soil layer as temperatures plummet below zero. The initiation of this increase coincides with the temperature drop below 0°C, while the conclusion is delayed relative to the temperature rise above 0°C, indicating a temporal lag. Figure 8. Distribution of seismic and meteorological stations for velocity variation and temperature comparison Black triangles represent seismic stations, red dots represent meteorological stations. Miao et al. (2019) propose that the velocity changes associated with seasonal soil freezing are governed by cumulative sub-zero temperatures and current temperatures, which can account for the lag in soil wave velocity change relative to winter temperature increases. They introduced a weighted cumulative freezing degree day to encompass these factors and presented an empirical model relating velocity changes to the weighted cumulative freezing degree days. On the nth day following the average daily temperature drop below 0°C, the weighted cumulative freezing degree day SW(n) is expressed as: (3) Where T1(n) is the daily average temperature on the nth day, p is the number of days with daily average temperatures below zero. H is the Heaviside step function (zero for negative arguments, one for nonnegative), andk is an empirical constant related to local site conditions. The relationship between Sw and velocity change dv/v can be formulated as (4) Where γ is the lower threshold for to be considered in the fitting, and Ra is an empirical constant related to the shear wave structure of the sedimentary layer. Ra ₌ α(1−1 λ) τ Vtop , with Vtop representing the seismic wave velocity of the soil layer, α the empirical coefficient linking freezing depth to accumulated sub-zero temperature, λ the ratio of seismic wave velocity change post-soil freezing in winter, and τ the time for seismic waves to propagate through the frozen soil layer. Employing this model, we attempted to establish a quantitative relationship between temperature and velocity changes. Eight stations were included in the calculation; FJ001 was excluded due to poor similarity between peak frequency and velocity change results in the first winter, and the significant distance of YGZ from the nearest meteorological observatory. We averaged the daily average temperatures from June 2021 to May 2024 to derive the average annual temperature curve and the velocity changes calculated in section 3.2 to obtain the annual average velocity variation curve. The fitting calculations based on formulas (3) and (4) required determining the empirical constants k, γ, and Ra for each station, which vary significantly due to site conditions. Lacking data on the seismic velocity structure and winter frozen soil depth under each station, we referred to Miao et al. (2019) and used grid search to identify these parameters, minimizing the least squares difference between the theoretical fitting curve and the velocity change results. The parameter k was varied in 0.2 steps from 0 to 6, γ in 1 steps from 0 to 41, and Ra in 0.01 steps from 0.01 to 5.00 to fit the amplitude of velocity change. This model only considers periods when temperatures fall below zero, so we focused on significant winter velocity increases, excluding other seasons. To align the velocity change results, we adjusted them so that the initial speed change value at the start of the fitting period was set to 0. The fitting results are presented in Figure 10. Except for JIH and QB001, the winter average velocity variations of the other six stations exhibit a “double peak” shape, similar to Steinmann et al. (2021). The velocity variation curves based on the model (red) for CAD, KB001, MG001, PW001, QC001, and WY001 fit this “double peak” pattern well. Although some short-term fluctuations are not perfectly captured, the model curves for JIH and QB001 also align well with the overall annual average velocity variation pattern. 4.2 Comparison with Precipitation In addition to temperature, we have observed that velocity changes at individual stations are also influenced by precipitation. Figure 11 presents a comparison between the velocity changes at the CAD station, the peak frequency of the HVSR curve, and the precipitation data from the nearest meteorological observatory. The study area is characterized by a subtropical monsoon climate, with intense rainfall occurring predominantly in July and August, and minimal rainfall during the winter months. From June 2021 to May 2024, there were a total of eight days when daily precipitation exceeded 50mm. On these heavy precipitation days and in the subsequent days, seismic velocity dropped rapidly (Figure 12a). An exception was observed on July 12, 2023, when precipitation reached 73.5mm, yet the velocity did not decrease significantly. The maximum daily precipitation over the three-year period was 101.6mm on July 30, 2023, the autocorrelation results indicated a seismic velocity decrease of approximately 5.6% in the following five days. On July 29 and 30, 2021, rainfall exceeded 60mm for two consecutive days, leading to a seismic wave velocity decrease of about 4.5% in the subsequent five days. In 2022 and 2023, the peak frequency of the HVSR also decreased on several occasions when daily precipitation exceeded 50mm (Figure 12b). However, the amplitude of the peak frequency decrease was not significantly different from its fluctuation amplitude outside of winter, suggesting that the rapid decline in peak frequency following heavy precipitation was not substantial. 341Yipei Tan, Ting Ma, Li Deng, Ke Xu Figure 9. Comparison of daily average temperature (blue), HVSR peak frequency (black) and relative velocity variation (red). 342 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area Figure 11. Comparison of precipitation with relative velocity variation (a) and HVSR peak frequency (b) at station CAD. Figure 10. Annual average velocity variation (black) and fitting result calculated based on formulas 3 and 4. Figure 12. Comparison the range of relative velocity changes in different time stacking of station KB001 (a) and CAD (b). 343Yipei Tan, Ting Ma, Li Deng, Ke Xu Roumelioti et al. (2020) observed a sharp decrease in the shear wave velocity of shallow sediments during the rainy season. Miao et al. (2018) discovered that significant velocity changes can be observed when rainfall exceeds a certain threshold, with no noticeable changes when it falls below this threshold. Comparing these findings with the results from station CAD in this study, a daily precipitation threshold of 50mm can be identified. In most cases exceeding this threshold, a significant decrease in velocity is observable, whereas it is challenging to detect a rapid decrease in velocity when it is below the threshold. It is also possible that minor velocity changes are obscured by fluctuations in velocity changes caused by temperature and calculation errors, making them difficult to discern. We additionally tested for correlations between velocity variations and air pressure data but did not find any relationship between the patterns of velocity variations and air pressure variations. 5 Discussions 5.1 Velocity Variation Range Numerous studies have reported seasonal relative velocity variations due to soil freezing, utilizing cross-correlation or autocorrelation of seismic noise. For instance, Gassenmeier et al. (2015) observed changes of several tenths of a percent in the 1.0-3.5 Hz frequency range in Germany. James et al. (2017) noted up to a 10 percent variation between summer and winter in the 13-17 Hz frequency range in America. Miao et al. (2019) recorded changes from about 1% to 9% across a frequency range from 1 to 13 Hz in Japan. Lindner et al. (2021) reported about 3% seasonal velocity variations in the 0.1-25 Hz frequency range in Germany. Liu et al. (2024) found about 8% variations in the 4-12 Hz frequency range in China. Cheng et al. (2022) observed 15% shear- wave velocity reductions using fiber-optic seismic networks. In comparison, the seismic velocity variation range obtained in this study is broader, exceeding 5%, with the largest variation range surpassing 20%. The study areas of Miao et al. (2019) and Liu et al. (2024) are further north and experience longer durations of sub-zero temperatures than ours, suggesting that the broader variation range in our results is not solely due to lower temperatures. The difference in stacking time when calculating daily Green’s functions could be one reason for the varying ranges of velocity variation. James et al. (2017) discussed the impact of stacking time on the variation range, finding that longer day-stacks yielded the smallest velocity reduction, while 1 and 2-day stacks yielded the largest reduction. We selected two stations, CAD and KB001, and calculated velocity changes using daily Green’s functions with 5-day, 15- day, and 30-day stacking. The results (Figure 12) indicate that the maximum velocity change amplitude decreases with increased stacking time, with the 5-day stacking results showing a larger maximum change amplitude than the 30-day stacking results over three years for both stations. Particularly during the winter of 2022-2023, the maximum change amplitude of the 5-day stacking results was approximately twice that of the 30-day stacking. The frequency band may also account for the differences observed. Gassenmeier et al. (2015), who used the lowest frequency range in the aforementioned studies, reported the smallest velocity variation range. Other research has noted greater variation ranges in higher frequency bands. In this paper, we used a 5-20 Hz filter when calculating Green’s functions, HVSR curves characteristics, which peak in the 5-10 Hz range in summer and 10-20 Hz in winter. By employing the same frequency band, the variations in peak frequency and seismic velocity reflect almost the same depth interval in the soil. The most dramatic changes in HVSR curves between winter and other seasons occur in the 5-20 Hz frequency range, which may lead to the broader variation range observed in our autocorrelation velocity change results. Laboratory measurements have shown significant velocity changes (>75%) between frozen and thawed soil (Zimmerman & King, 1986). None of the aforementioned studies captured such substantial changes. James et al. (2017) argue that their results only recovered 2-15% velocity decreases because only a portion of the wave energy traveled through the frozen layer. Our use of a broader frequency band may enable our results to encompass velocity changes within the 5-20% range. Chen et al. (2023) demonstrated that the peak frequency of the HVSR curve increased by 3% to 35%, and the amplification factor decreased by 11% to 38%. Kula et al. (2018) also showed that soil freezing significantly impacts the high-frequency part of the HVSR curve, with peak amplitudes notably greater in summer than in winter. Köhler and Weidle (2019) observed an HVSR peak emerging at the beginning of summer and gradually shifting to lower frequencies in Norway. These results are consistent with the seasonal variation characteristics of the HVSR curve found in our study. The reduction in site amplification factor due to soil freezing in winter may imply that the amplification effect of sediment on seismic waves is diminished when temperatures drop below zero. Consequently, the amplitude of ground motion caused by the same source in the higher frequency band (5-20 Hz) is smaller than in other seasons, potentially alleviating earthquake disasters in winter. 5.2 Velocity Variation Characteristic Gassenmeier et al. (2015) and James et al. (2017) observed a sharp increase in seismic velocity in early winter, followed by a drastic decrease in velocity in later spring. Miao et al. (2019) and Liu et al. (2024) noted a smoother velocity rise during the winter months. In this study, we observed a sharp increase in early winter and a decrease in later spring, but with a distinctive “double peak” pattern in the winter velocity variations. There are two potential reasons for this “double peak” shape (Figure 10). First, the daily average temperature in winter does not remain consistently below 0°C but fluctuates around zero. As illustrated in Figure 9 of Steinmann et al. (2021), temperatures rise above 0°C after falling below it for a period, then drop again, leading to the weighted cumulative freezing degree day SW(n) decreasing and then increasing, thus forming a “double peak” in seismic velocity change. Second, the temperature change pattern varies each year. The time for winter temperatures to fall below 0°C in 2022-2023 is significantly later than in the previous two years. The lowest winter temperature in 2022-2023 appears in the second half of winter (end of February), whereas in 2023-2024 it appears early in winter (beginning of January), leading to different peak times for seismic wave velocity changes in these years (this phenomenon is more pronounced at PW001 station). When the velocity changes are averaged over three years, a “double peak” shape emerges. Compared to Miao et al. (2019), who studied velocity changes due to the soil freeze-thaw cycle in Hokkaido, Japan, there is a notable difference in our results. Miao et al. (2019) found that after the temperature dropped below 0°C, velocity slightly increased before falling back, then undergoing a significant change. They interpreted this as the release of heat by water in the soil layer as it froze, keeping the soil layer temperature around 0°C for a short time. In our results, this minor velocity increase is not common but can be observed at station QC001 (Figure 10). Seismic velocity at other stations rises rapidly after the temperature falls below 0°C. Careful comparison of the temperature changes at QC001 with the minor velocity increase reveals that the temperature does not drop monotonically at the beginning of winter but fluctuates around 0°C, rising to 0°C after a few days below, then dropping again. The timing of these temperature fluctuations matches the slight velocity increase. Steinmann et al. (2021) also observed multiple peaks in seismic velocity variation due to temperature fluctuations around 0°C. We speculate that these minor velocity increases may be caused by such temperature variations. Kula et al. (2018)using the HVSR method to study periodic variations in frozen soil characteristics at the Finnish Polar Research Station, compared ground temperature observations at different depths with the HVSR curve’s peak value and found good agreement between the 1m underground temperature variation and the HVSR curve’s peak value. The 1m underground temperature pattern is similar to air temperature but lags in phase, suggesting a lag in the HVSR curve’s peak value change relative to air temperature. Since the HVSR curve’s peak frequency change closely mirrors velocity changes, it can be inferred that velocity changes also lag air temperature. This lag may adversely affect the fitting results of velocity variations and weighted cumulative freezing degree days based on air temperature (Figure 10). In our study area, we have not found underground temperature observation data at several meters depth, preventing a quantitative evaluation of the underground temperature lag relative to air temperature. We plan to conduct low-temperature gradient measurements at different depths to analyze the relationship between ground temperature changes and seismic wave velocity changes in the soil layer. In this paper, the autocorrelation function obtained from the horizontal component shows significant differences between winter and summer, a difference not easily observed in the vertical component results (Figure 5). 344 Seasonal variations characteristics of shallow subsurface seismic velocity and its influencing factors in Tianjin area Steinmann et al. (2021) noted this difference and suggested it was due to the increase in Love wave velocity caused by soil freezing in winter being greater than that of Rayleigh waves. Chen et al. (2023) found that seasonal frozen soil increased the resonant frequency of the horizontal component of the ground motion spectrum but had no significant effect on the vertical component. Therefore, in areas where temperatures drop below zero in winter, horizontal component seismic records should be used to analyze temporal variations in seismic velocity. The seismic velocity of soft soil sediment increases significantly when temperatures fall below 0°C, and the HVSR curve’s shape also changes markedly below 0°C. Therefore, whether calculating the HVSR curve from field ground tremor observations or detecting shallow subsurface structures based on ambient seismic noise from dense arrays, observations should be avoided when temperatures are close to or below 0°C. 6. Conclusions This paper presents an analysis of shallow subsurface seismic velocity variation characteristics in Tianjin, in North China Plain, based on three years of continuous seismic waveform records. Utilizing the Horizontal-to-Vertical Spectral Ratio (HVSR) and autocorrelation methods, we have investigated the possible factors influencing these variations. The primary conclusions are as follows: 1. At 25 stations within the study area, the HVSR curves exhibit seasonal variation in the 5-20 Hz frequency range, with an increase in peak frequency and a decrease in peak amplitude during the winter months. Correspondingly, the seismic velocity calculated using the autocorrelation method at 10 stations also shows an increase during winter. 2. By comparing the peak frequency and seismic velocity variation at 10 stations with temperature records, we find that the seasonal variation patterns of peak frequency and seismic velocity are closely aligned. The rapid increase in November corresponds to the period when the daily average temperature falls below 0°C, and the subsequent decline in February corresponds to when temperatures rise above 0°C. Temperature is deemed the primary factor driving the seasonal variation in velocity. 3. Employing an empirical model that relates velocity variation to weighted cumulative freezing degree days, we find that the model curve effectively captures the overall pattern of annual average velocity variation. We infer that the increase in shear wave velocity and peak frequency in the study area during winter is due to the freezing of the soil layer. 4. The seismic velocity of the shallow subsurface undergoes significant changes in winter compared to other seasons. When detecting deep underground structures, the impact of temperature on the reliability of the results should be taken into consideration. Acknowlegments This work is supported by Tianjin Science and Technology Foundation (No. 22JCYBJC00100 & 23JCYBJC00500). References Ashruf, T. N., & Morelli, A. (2022). The Moho reflectivity of the subduction be- neath the Southwestern Alps from ambient seismic noise autocorrela- tions. Geophysical Journal International, 230(1), 298-316. https://doi. org/10.1093/gji/ggac079 Benkaci, N., Oubaiche, E. H., Chatelain, J. L., Bensalem, R., Benouar, D., & Ab- bes, K. (2018). Non-stability and non-reproducibility of ambient vibra- tion HVSR peaks in Algiers (Algeria). Journal of Earthquake Enginee- ring, 25(5), 853-871. https://doi.org/10.1080/13632469.2018.1537903 Bonilla, L. F., Guéguen, P., & Ben-Zion, Y. (2019). Monitoring coseismic tem- poral changes of shallow material during strong ground motion with interferometry and autocorrelation. Bulletin of the Seismological Society of America, 109(1), 187-198. https://doi.org/10.1785/0120180092 Cao, J., Yan, C., & Zhang, W. (2018). Report of fracture detection and seismic risk assessment of Jiyunhe fault (in Chinese). Tianjin Earthquake Agency, 306-315. Chen, S., Lei, J., & Li, Y. (2023). Microtremor Recording Surveys to Study the Effects of Seasonally Frozen Soil on Site Response. Sensors, 23(12), 5573. https://doi.org/10.3390/s23125573 Cheng, F., Lindsey, N. J., Sobolevskaia, V., Dou, S., Freifeld, B., Wood, T., James, S. R., Wagner, A. M., & Ajo-Franklin, J. B. (2022). Watching the Cr- yosphere Thaw: Seismic Monitoring of Permafrost Degradation Using Distributed Acoustic Sensing During a Controlled Heating Experiment. Geophysical Research Letters, 49(10), e2021GL097195. https://doi.or- g/10.1029/2021GL097195 Delouche, E., & Stehly, L. (2023). Seasonal Seismic Velocity Variations Measured Using Seismic Noise Autocorrelations to Monitor the Dynamic of Aqui- fers in Greece. Journal of Geophysical Research: Solid Earth, 128(12), e2023JB026759. https://doi.org/10.1029/2023JB026759 Dong, Y., Ni, S., Yuen, D. A., & Li, Z. (2018). Crustal rheology from focal depths in the North China Basin. Earth and Planetary Science Letters, 497, 123- 138. https://doi.org/10.1016/j.epsl.2018.06.018 Feng, F., Huang, H., Hsu, J., & Wu, M. (2021). Controls on Seasonal Variations of Crustal Seismic Velocity in Taiwan Using Single-Station Cross-Compo- nent Analysis of Ambient Noise Interferometry. Journal of Geophysical Research: Solid Earth, 126(11), e2021JB022650. https://doi.org/10.1029/ 2021JB022650 Gradon, C., Brenguier, F., Stammeijer, J., Mordret, A., Hindriks, K., Campman, X., Lynch, R., Boué, P., & Chmiel, M. (2021). Seismic velocity response to atmospheric pressure using time-lapse passive seismic interferome- try. Bulletin of the Seismological Society of America, 111(6), 3451–3458. https://doi.org/10.1785/0120210069 Gassenmeier, M., Delatre, M., & Korn, M. (2015). Monitoring of environmental influences on seismic velocity at the geological storage site for CO2 in Ketzin (Germany) with ambient seismic noise. Geophysical Journal In- ternational, 200(1), 524-533. https://doi.org/10.1093/gji/ggu413 James, S. R., Knox, H. A., Abbott, R. E., & Screaton, E. J. (2017). Improved mo- ving window cross-spectral analysis for resolving large temporal seis- mic velocity changes in permafrost. Geophysical Research Letters, 44(9), 4018-4026. https://doi.org/10.1002/2016GL072468 Konno, K., & Ohmachi, T. (1998). Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of mi- crotremor. Bulletin of the Seismological Society of America, 88(1), 228- 241. https://doi.org/10.1785/bssa0880010228 Köhler A., & Weidle, C. (2019). Potentials and pitfalls of permafrost active layer monitoring using the HVSR method: a case study in Svalbard. Earth Surface Dynamics, 7, 1-16. https://doi.org/10.5194/esurf-7-1-2019 Kramer, R., Lu, Y., & Bokelmann, G. (2023). Interaction of Air Pressure and Groundwater as Main Cause of Sub-Daily Relative Seismic Velocity Changes. Geophysical Research Letters, 50(7), e2022GL101298. https:// doi.org/10.1029/2022GL101298 Kula, D., Olszewska, D., Dobiński, W., & Glazer, M. (2018). Horizontal-to-ver- tical spectral ratio variability in the presence of permafrost. Geophysi- cal Journal International, 214(1), 219-231. https://doi.org/10.1093/gji/ ggy118 Langston, C. A. (2011). Wavefield continuation and decomposition for passi- ve seismic imaging under deep unconsolidated sediments. Bulletin of the Seismological Society of America, 101(5), 2176-2190. https://doi. org/10.1785/0120100299 La Rocca, M., Chiappetta, G. D., Gervasi, A., & Festa, R. L. (2020). Non-stability of the noise HVSR at sites near or on topographic heights. Geophysical Journal International, 222(3), 2162-2171. https://doi.org/10.1093/gji/ ggaa297 La Rocca, M., & Chiappetta, G. D. (2022). Day–night cycle of seismic noise HVSR and comparison with body waves and T waves. Geophysical Journal In- ternational, 231(3), 1535-1544. https://doi.org/10.1093/gji/ggac265 Lecocq, T., Caudron, C., & Brenguier, F. (2014). MSNoise, a Python Packa- ge for Monitoring Seismic Velocity Changes Using Ambient Seis- https://doi.org/10.1093/gji/ggac079 https://doi.org/10.1093/gji/ggac079 https://doi.org/10.1080/13632469.2018.1537903 https://doi.org/10.1785/0120180092 https://doi.org/10.3390/s23125573 https://doi.org/10.1029/2021GL097195 https://doi.org/10.1029/2021GL097195 https://doi.org/10.1029/2023JB026759 https://doi.org/10.1016/j.epsl.2018.06.018 https://doi.org/10.1029/2021JB022650 https://doi.org/10.1029/2021JB022650 https://doi.org/10.1785/0120210069 https://doi.org/10.1093/gji/ggu413 https://doi.org/10.1002/2016GL072468 https://doi.org/10.1785/bssa0880010228 https://doi.org/10.5194/esurf-7-1-2019 https://doi.org/10.1029/2022GL101298 https://doi.org/10.1029/2022GL101298 https://doi.org/10.1093/gji/ggy118 https://doi.org/10.1093/gji/ggy118 https://doi.org/10.1785/0120100299 https://doi.org/10.1785/0120100299 https://doi.org/10.1093/gji/ggaa297 https://doi.org/10.1093/gji/ggaa297 https://doi.org/10.1093/gji/ggac265 345Yipei Tan, Ting Ma, Li Deng, Ke Xu mic Noise. Seismological Research Letters, 85(3), 715-726. https://doi. org/10.1785/0220130073 Li, G., & Ben-Zion, Y. (2023). Daily and Seasonal Variations of Shallow Seismic Velocities in Southern California From Joint Analysis of H/V Ratios and Autocorrelations of Seismic Waveforms. Journal of Geophysical Research: Solid Earth, 128(2), e2022JB025682. https://doi.org/10.1029/ 2022JB025682 Li, J., Song, X., Yang, Y., Li, M., Li, J., & Li, Y. (2021). Strong Seasonal Variations of Seismic Velocity in Eastern Margin of Tibetan Plateau and Sichuan Basin from Ambient Noise Interferometry. Journal of Geophysical Re- search: Solid Earth, 126(11), e2021JB022600. https://doi.org/10.1029/ 2021JB022600 Lindner, F., Wassermann, J., & Igel, H. (2021). Seasonal Freeze-Thaw Cycles and Permafrost Degradation on Mt. Zugspitze (German/Austrian Alps) Re- vealed by Single-Station Seismic Monitoring. Geophysical Research Le- tters, 48(18), e2021GL094659. https://doi.org/10.1029/2021GL094659 Liu, H., Li, J., Hu, R., Meng, H., & Lyu, H.  (2024). Quantitatively Monitoring of Seasonal Frozen Ground Freeze–Thaw Cycle Using Ambient Seismic Noise Data. Seismological Research Letters, 96(1), 282-293. https://doi. org/10.1785/0220240201 Liu, Z., Liang, C., Huang, H., Wang, C., & Cao, F. (2022). Seismic Velocity Va- riations at Different Depths Reveal the Dynamic Evolution Associated With the 2018 Kilauea Eruption. Geophysical Research Letters, 49(3), e2021GL093691. https://doi.org/10.1029/2021GL093691 Lontsi, A. M., Hobiger, M., Panzera, F., Sánchez-Sesma, F. J., & Fäh, D. (2022). Seismic characterization of Swiss strong-motion borehole-station sites by inversion of full microtremor horizontal-to-vertical spectral ratios [H/V(z,f)]. Bulletin of the Seismological Society of America, 113, 417- 436. https://doi.org/10.1785/0120210320 Luan, Y., Yang, H., Wang, B., Yang, W., Wang, W., Yang, J., & Li, X. (2022). Ti- me-lapse monitoring of daily velocity changes in Binchuan, southwes- tern China, using large-volume air-gun source array data. Seismological Research Letters, 93, 914-930. DOI: 10.1785/0220210160. Mao, S., Lecointre, A., D., R., & Campillo, M. (2022). Space-time monitoring of groundwater fluctuations with passive seismic interferometry. Na- ture Communications, 13(1), 1-9. https://doi.org/10.1038/s41467-022- 32194-3 Miao, Y., Shi, Y., & Wang, S. (2018). Temporal change of near-surface shear wave velocity associated with rainfall in Northeast Honshu, Japan. Earth, Pla- nets and Space, 70(1), 1-11. https://doi.org/10.1186/s40623-018-0969-3 Miao, Y., Shi, Y., Zhuang, H. Y., Wang, S. Y., Liu, H. B., & Yu, X. B. (2019). In- fluence of Seasonal Frozen Soil on Near-Surface Shear Wave Velocity in Eastern Hokkaido, Japan. Geophysical Research Letters, 46(16), 9497- 9508. https://doi.org/10.1029/2019GL082282 Nakamura, Y. (1989). A method for dynamic characteristics estimation of subsur- face using microtremor on the ground surface. Railway Technical Re- search Institute, 30, 25-33. Nakamura, Y. (2019). What is the Nakamura method? Seismological Research Le- tters, 90(4), 1437–1443. https://doi.org/10.1785/0220180376 S. Oakley, D. O., Forsythe, B., Gu, X., Nyblade, A. A., & Brantley, S. L. (2021). Seismic Ambient Noise Analyses Reveal Changing Temperature and Water Signals to 10s of Meters Depth in the Critical Zone. Journal of Geophysical Research: Earth Surface, 126(2), e2020JF005823. https://doi. org/10.1029/2020JF005823 Okubo, K., Delbridge, B. G., & Denolle, M. A. (2024). Monitoring Velocity Chan- ge Over 20 Years at Parkfield. Journal of Geophysical Research: Solid Ear- th, 129(4), e2023JB028084. https://doi.org/10.1029/2023JB028084 Peng, F., Wang, W., & Kou, H. D. (2020). Microtremer H/V spectral ratio inves- tigation in the Sanhe-Pinggu area: site responses, shallow sedimentary structure, and fault activity revealed. Chinese Journal of Geophysics (in Chinese), 63(10), 3775-3790. DOI: 10.6038/cjg2020O0025 Poli, P., Marguin, V., Wang, Q., & Johnson, P. (2020). Seasonal and Coseismic Velocity Variation in the Region of L’Aquila From Single Station Measu- rements and Implications for Crustal Rheology. Journal of Geophysical Research: Solid Earth, 125(7), e2019JB019316. https://doi.org/10.1029/ 2019JB019316 Qin, L., Steidl, J. H., Qiu, H., Nakata, N., & Ben-Zion, Y. (2022). Monitoring Sea- sonal Shear Wave Velocity Changes in the Top 6 m at Garner Valley in Southern California With Borehole Data. Geophysical Research Letters, 49(23), e2022GL101189. https://doi.org/10.1029/2022GL101189 Qin, T., Lu, L., Ding, Z., Feng, X., & Zhang, Y. (2022). High-Resolution 3D Sha- llow S Wave Velocity Structure of Tongzhou, Subcenter of Beijing, Infe- rred From Multimode Rayleigh Waves by Beamforming Seismic Noise at a Dense Array. Journal of Geophysical Research: Solid Earth, 127(5), e2021JB023689. https://doi.org/10.1029/2021JB023689 Rigo, A., Sokos, E., Lefils, V., & Briole, P. (2021). Seasonal variations in amplitu- des and resonance frequencies of the HVSR amplification peaks linked to groundwater. Geophysical Journal International, 226(1), 1-13. https:// doi.org/10.1093/gji/ggab086 Roumelioti, Z., Hollender, F., & Gueguen, P. (2020). Rainfall-induced variation of seismic waves velocity in soil and implications for soil response: What the ARGONET (Cephalonia, Greece) vertical array data reveal. Bulle- tin of the Seismological Society of America, 110(2), 441-451. https://doi. org/10.1785/0120190183 Seivane, H., García-Jerez, A., Navarro, M., Molina, L., & Navarro-Martínez, F. (2022). On the use of the microtremor HVSR for tracking velocity changes: a case study in Campo de Dalias basin (SE Spain). Geophysical Journal International, 230, 542-564. https://doi.org/10.1093/gji/ggac064 Sens-Schönfelder, C., & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indo- nesia. Geophysical Research Letters, 33(21). https://doi.org/10.1029/ 2006GL027797 Steinmann, R., Hadziioannou, C., & Larose, E. (2020). Effect of centimetric freezing of the near subsurface on Rayleigh and Love wave velocity in ambient seismic noise correlations. Geophysical Journal International, 224(1), 626-636. https://doi.org/10.1093/gji/ggaa406 Tchawe, F. N., Froment, B., Campillo, M., & Margerin, L. (2020). On the use of the coda of seismic noise autocorrelations to compute H/V spectral ratios. Geophysical Journal International, 220(3), 1956-1964. https://doi. org/10.1093/gji/ggz553 Tianjin Institute of Geological Survey. (2018). Regional geology of Tianjin Muni- cipality. Beijing: Geological Publishing House, 279-324. Vassallo, M., Cultrera, G., Giulio, G. D., Cara, F., & Milana, G. (2022). Peak Fre- quency Changes From HV Spectral Ratios in Central Italy: Effects of Strong Motions and Seasonality Over 12 Years of Observations. Journal of Geophysical Research: Solid Earth, 127(5), e2021JB023848. https:// doi.org/10.1029/2021JB023848 Wang, Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., & Aoki, Y. (2017). Seasonal Crustal Seismic Velocity Changes Throughout Japan. Journal of Geophysical Research: Solid Earth, 122(10), 7987-8002. https://doi.or- g/10.1002/2017JB014307 Weaver, R. L., Hadziioannou, C., Larose, E., & Campillo, M. (2011). On the precision of noise correlation interferometry. Geophysical Journal International, 185(3), 1384-1392. https://doi.org/10.1111/j.1365- 246X.2011.05015.x Yassminh, R., Gallegos, A., Sandvol, E., & Ni, J. (2019). Investigation of the re- gional site response in the Central and Eastern United States. Bulletin of the Seismological Society of America, 109, 1005–1024. https://doi. org/10.1785/0120180230 Zhang, S., Luo, B., Ben-Zion, Y., Lumley, D. E., & Zhu, H. (2023). Monitoring Terrestrial Water Storage, Drought and Seasonal Changes in Central Oklahoma With Ambient Seismic Noise. Geophysical Research Letters, 50(17), e2023GL103419. https://doi.org/10.1029/2023GL103419 Zimmerman, R. W., King, M. S. (1986). The effect of the extent of freezing on seismic velocities in unconsolidated permafrost. Geophysics, 51(6), 1285–1290. https://doi.org/10.1190/1.1442181 https://doi.org/10.1785/0220130073 https://doi.org/10.1785/0220130073 https://doi.org/10.1029/2022JB025682 https://doi.org/10.1029/2022JB025682 https://doi.org/10.1029/2021JB022600 https://doi.org/10.1029/2021JB022600 https://doi.org/10.1029/2021GL094659 https://doi.org/10.1785/0220240201 https://doi.org/10.1785/0220240201 https://doi.org/10.1029/2021GL093691 https://doi.org/10.1785/0120210320 https://doi.org/10.1038/s41467-022-32194-3 https://doi.org/10.1038/s41467-022-32194-3 https://doi.org/10.1186/s40623-018-0969-3 https://doi.org/10.1029/2019GL082282 https://doi.org/10.1785/0220180376 https://doi.org/10.1029/2020JF005823 https://doi.org/10.1029/2020JF005823 https://doi.org/10.1029/2023JB028084 https://doi.org/10.6038/cjg2020O0025 https://doi.org/10.1029/2019JB019316 https://doi.org/10.1029/2019JB019316 https://doi.org/10.1029/2022GL101189 https://doi.org/10.1029/2021JB023689 https://doi.org/10.1093/gji/ggab086 https://doi.org/10.1093/gji/ggab086 https://doi.org/10.1785/0120190183 https://doi.org/10.1785/0120190183 https://doi.org/10.1093/gji/ggac064 https://doi.org/10.1029/2006GL027797 https://doi.org/10.1029/2006GL027797 https://doi.org/10.1093/gji/ggaa406 https://doi.org/10.1093/gji/ggz553 https://doi.org/10.1093/gji/ggz553 https://doi.org/10.1029/2021JB023848 https://doi.org/10.1029/2021JB023848 https://doi.org/10.1002/2017JB014307 https://doi.org/10.1002/2017JB014307 https://doi.org/10.1111/j.1365-246X.2011.05015.x https://doi.org/10.1111/j.1365-246X.2011.05015.x https://doi.org/10.1785/0120180230 https://doi.org/10.1785/0120180230 https://doi.org/10.1029/2023GL103419 https://doi.org/10.1190/1.1442181 OLE_LINK530 OLE_LINK531 OLE_LINK532 OLE_LINK533 OLE_LINK9 OLE_LINK10 OLE_LINK11 OLE_LINK12 OLE_LINK1 OLE_LINK2 OLE_LINK526 OLE_LINK525