## 2 Materials and Methods
### Surface Seismic Analysis  
Seismic data acquired in permafrost environments often show strong surface waves [12] and, therefore, we expect to observe these in both synthetic and real seismic data. Vertical geophones record surface waves of Rayleigh type. These propagate in the uppermost sediments of the subsurface as an interaction phenomenon between P- and S-waves. The phase velocity of a Rayleigh wave depends on the elastic properties of the medium that the wave travels in, which for a homogenous medium corresponds to the solution of the following equation [36]:  
\[\zeta^{3}\,-\,8\zeta^{2}\,+\,8\zeta(3\,-\,2\eta)\,-\,16(1\,-\,\eta)=0, \tag{7}\]  
where \(\zeta=\frac{c_{R}^{2}}{V_{s}^{2}}\) and \(\eta=\frac{V_{g}^{2}}{V_{p}^{2}}\). \(c_{R}\) is phase velocity, and \(V_{p}\) and \(V_{s}\) are the P- and S- wave velocities defined by Equations (5) and (6).  
There is no frequency dependence inherent in Equation (7), meaning that for a wave propagating in a homogenous medium, all frequency components of the surface wave travel with the same phase velocity \(c_{R}\). However, when a wave propagates along the surface, the amplitude of higher frequencies decreases more rapidly with depth than the lower frequencies; thus, low frequency energy is generally dominating at larger depths. If the elastic properties of the subsurface vary with depth, the various frequency components are influenced by different effective elastic properties, and \(c_{R}\) becomes a function of depth and frequency. The equations describing surface wave propagation can then have several  
Figure 3: Variation in temperature and seismic velocities with depth for a specific point in time (after four years) for models T0 (initial model, \(-5\)\({}^{\circ}\)C through the column), T1 (top boundary temperature = +1 \({}^{\circ}\)C), and T3 (top boundary temperature = +10 \({}^{\circ}\)C). (**a**) Temperature profiles, (**b**) P-wave velocity profiles, and (**c**) S-wave velocity profiles. Temperature is equal at all depths for T0, so the observed velocity variation with depth is only due to the porosity decrease: 40% in the upper layer, 35% in the second layer, 30% in the third layer, and 20% in the bottom layer.  
solutions (i.e., possible combinations of frequencies and phase velocities), called modes, where the fundamental mode is the solution with the lowest phase velocity for a given frequency [37]. Since the Rayleigh wave is not inherently dispersive, any observed dispersion must be due to varying elastic properties of the near-surface sediments. Equation (7) infers that \(c_{R}\) depends on the sediment shear modulus through \(V_{s}\), which has been shown to strongly respond to freezing/thawing. Thus, an inversion of frequency-dependent surface wave phase velocity for obtaining S-wave velocity with depth can reveal the degree of freezing with depth if the sediment composition is known [38].  
To analyze surface wave dispersion, we generate frequency-phase velocity spectra (FV-spectra) by using a wavefield transform. This transform builds an FV-spectrum by decomposing the gather into separate frequency components through Fourier transformation, applying amplitude normalization to each component, and scanning for a range of phase velocities to find the magnitude of summed amplitudes for each frequency component [39]. The FV-spectrum is useful for separating body and surface wave energy, and when the surface wave phase velocity is constant for all frequencies, there is no variation in the elastic properties with depth. Correspondingly, if the FV-spectra of the surface wave vary between equal experiments conducted at different times, this implies a temporal variation of the elastic properties. In our case, this time-lapse effect represents a signal of altered freezing conditions.