## 2 Materials and Methods
### Modeling Effects of Thawing on Seismic Properties  
To understand the effects of thawing and freezing on the elastic properties of unconsolidated, fully saturated near-surface sediments, we refer to papers of [PERSON] et al. [25; 26; 27] who described the effects of cement (here: ice) forming at various pore-scale configurations on the effective elastic properties of sediments. Ice can form at the grain contacts where it significantly increases the elastic properties, or within the pore space where the effect on effective properties is smaller. Since seismic experiments are dominated by low frequencies (usually \(<\)200 Hz), the relevant wavelengths are orders of magnitude larger than the heterogeneities of the sediments. To convert sediment composition and texture to elastic properties, we can adopt standard rock physics principles. In our case, we will consider isotropic sediments, which then elastically are fully characterized by the bulk modulus \(K\) and the shear modulus \(\mu\). We use a rock physics model similar to the two-end member mixing approach used by [PERSON], et al. [28] and [PERSON], et al. [29] that inherently assumes cementing and pore-filling ice forming simultaneously at all ice saturations. We compute elastic properties of a fully unfrozen composite (\(\theta_{i}=0\)) by first using Hertz-Mindlin contact theory [30] to compute dry rock properties, and subsequently Gassmann fluid substitution [31] to compute fully water-saturated rock properties. Furthermore, we compute elastic properties of a fully frozen composite (\(\theta_{f}=0\)) using  
Figure 1: Freezing point isotherm throughout four years for (**a**) high porosity sand, and (**b**) low porosity sand. Both models are displayed for low salinity (dashed lines, freezing point temperature is \(-\)0.17 \({}^{\circ}\)C) and high salinity (solid lines, freezing point temperature is \(-\)1.43 \({}^{\circ}\)C ) pore water, and for models T1 (top boundary temperature = +1 \({}^{\circ}\)C ), T2 (top boundary temperature = +5 \({}^{\circ}\)C ), and T3 (top boundary temperature = +10 \({}^{\circ}\)C ). the self-consistent approach (SCA) as described by [32; 33] by assuming penny-shaped ice inclusions and spherical sediment grain inclusions in a background medium with initially unknown properties. Lastly, to find elastic properties of a composite with ice saturation varying from 0% to 100%, we use the Hill average of the modified Hashin-Shtrikman bounds to mix the fully unfrozen and the fully frozen end-members [34]. For completeness, the main procedure for computing elastic properties of thawing or freezing sediments is outlined in Appendix A. The resulting elastic properties can together with density \(\rho\) be input to the equations for seismic P- and S-wave velocities \(V_{p}\) and \(V_{s}\), respectively:  
\[V_{p}=\sqrt{\left(K+\frac{4}{3}\;\mu\right)/\rho}, \tag{5}\]  
and  
\[V_{s}=\sqrt{\mu/\rho}. \tag{6}\]  
Combining the resulting velocities with Equation (1), we can determine how seismic velocities vary with time and depth. The complete approach of combing heat flux modeling, rock physics modeling, and the equations for seismic velocities to investigate how seismic velocities vary with time is described in detail in [35]. Figure 2 shows the result of using this approach on high-porosity sand with high-salinity pore water, considering scenarios T1, T2, and T3 (solid lines in Figure 1a). The figure shows how velocities at two stationary points in the subsurface: 2 m depth (Figure 2b,c) and 10 m depth (Figure 2e,f), vary with time. The resulting velocity variations are consistent with the temperature variation with time at these two depths (Figure 2a,d): as temperatures increase, velocities decrease because the water content increases. The largest velocity change is at 2 m depth, consistent with the larger temperature change at 2 m depth than at 10 m depth because less heat is conducted down to 10 m depth.