Potential Use of Time-Lapse Surface Seismics for Monitoring Thawing of the Terrestrial Arctic

Featured Application

We investigate the feasibility of using seismic data acquired in repeated active experiments for monitoring changes in the degree of freezing in unconsolidated sediments.

Abstract

The terrestrial Arctic is warming rapidly, causing changes in the degree of freezing of the upper sediments, which the mechanical properties of unconsolidated sediments strongly depend upon. This study investigates the potential of using time-lapse surface seismics to monitor thawing of currently (partly) frozen ground utilizing synthetic and real seismic data. First, we construct a simple geological model having an initial temperature of −5 °C, and infer constant surface temperatures of −5 °C, +1 °C, +5 °C, and +10 °C for four years to this model. The geological models inferred by the various thermal regimes are converted to seismic models using rock physics modeling and subsequently seismic modeling based on wavenumber integration. Real seismic data reflecting altered surface temperatures were acquired by repeated experiments in the Norwegian Arctic during early autumn to mid-winter. Comparison of the surface wave characteristics of both synthetic and real seismic data reveals time-lapse effects that are related to thawing caused by varying surface temperatures. In particular, the surface wave dispersion is sensitive to the degree of freezing in unconsolidated sediments. This demonstrates the potential of using surface seismics for Arctic climate monitoring, but inversion of dispersion curves and knowledge of the local near-surface geology is important for such studies to be conclusive.

1. Introduction

Adventdalen is located on Svalbard in the Norwegian Arctic, in an area that traditionally has experienced stable subzero winter temperatures and minimal annual rainfall. However, measurements in recent years show that climate is changing rapidly in this region and in the Arctic in general [1]. The IPCC [2] projects that temperatures and precipitation will continue to increase drastically in the years to come, likely contributing to thawing of currently frozen surfaces, including sea ice, glaciers, snow, and fully or partly frozen ground, often referred to as permafrost. Permafrost is defined as ground that stays at or below 0 °C for at least two years [3], and this purely thermal definition implies that the degree of freezing in permafrost can vary considerably.

Methods for monitoring thawing of permafrost are desirable because thawing of currently frozen ground may have large geomorphic consequences [4]. Frozen ground is stiffer and more rigid than unfrozen ground [5], therefore, subsidence and increased avalanche risk are likely consequences of thawing. Additionally, the large currently frozen land areas in the Northern hemisphere hold huge reserves of carbon that may be released to the atmosphere in the thawing process, inducing further climate change [6]. These methods should preferably be non-intrusive due to efficiency and environmental constraints.

Kneisel, et al. [7] provide an overview of some previously tested geophysical methods for monitoring thawing of permafrost, based on changes in electric, electromagnetic, elastic, or dielectric properties. All methods show limitations, particularly in areas with saline pore water, see [8] for a good overview. Salt causes freezing point depression and gradual freezing of pore water at subzero temperatures because the salinity of the residual pore water increases as salt is expelled when the water freezes [9]. Typical for saline permafrost is therefore unfrozen water within the permafrost, the amount depending on the salinity. If the salinity varies with depth, the amount of ice also varies with depth leading to an irregularly varying stiffness gradient with depth. For such partially frozen soils, the correlation between temperature and resistivity can be complicated [8,10], making electric methods difficult to use. Additionally, refraction seismics cannot be used in such environments, because it will not detect softer layers if these are located below stiffer layers.

Effective elastic properties of a medium depend on the physical properties and quantity of its constituents, as well as their geometrical distribution, which for permafrost typically means sediment grains, water, ice, and air. Thus, since thawing of a frozen material leads to a change in the ice-to-water ratio, referred to as the ice saturation, thawing will change the effective elastic and seismic properties. Laboratory measurements of the seismic properties of permafrost during the freeze-thaw process [8,11] support the idea that elastic properties of unconsolidated sediments decrease with increased thawing, however the relation between the degree of thawing and change in elastic properties is not linear.

Shear modulus of unconsolidated sediments is particularly affected by variation in ice saturation [12], and since surface wave properties strongly depend on the shear modulus, surface wave methods may be a viable option for detecting thawing in permafrost. Multichannel analysis of surface waves (MASW) was first introduced by Park, et al. [13], and has been extensively used for determining the properties of near-surface sediments through inversion of surface wave properties. However, MASW inversion shows limited success for soils with irregular stiffness gradients, and so more sophisticated methods are needed here. Several studies have focused on modified versions of MASW or other surface-wave methods to detect low velocity layers within permafrost with promising results [14,15], since the surface wave properties depend on how the elastic properties of the near-surface sediments vary with depth. For example, Tsuji, et al. [16] detected unfrozen sediments at the base of a glacier on Svalbard using common-midpoint cross-correlation followed by MASW. However, few have used surface wave methods to detect variations in permafrost with time, so-called time-lapse surface seismics. Ajo-Franklin, et al. [17] conducted a controlled warming experiment using distributed acoustic sensing to study the seismic response of such thawing, but we are not aware of any published studies focusing on variations in the seismic signatures of surface waves caused only by natural temperature variations.

The objective of this study is to evaluate the potential of using seismic data in general, and dispersion properties of surface waves in particular, to monitor temporal changes in the permafrost of the Arctic region caused by surface temperature variation. To study this, we make use of both synthetic seismic data generated by seismic modeling, and real time-lapse seismic data obtained by experiments in Adventdalen.

2. Materials and Methods

2.1. Modeling Thawing of Frozen Sediments

Heat transport in the subsurface occurs through conduction or convection, where conduction is commonly believed to be the dominant process in permafrost environments except for during extensive spring thawing [18,19]. Heat is conducted downward into the ground if surface temperature is higher than ground temperature, and opposite. Surface temperature is primarily controlled by air temperature, therefore heat flux into the ground is typically largest during summer, but other factors such as snow cover and vegetation also affect surface temperature [3]. If the ground is initially frozen, which is common in the Arctic, this heat flux may cause thawing. In the Arctic today, winter temperatures are generally cold enough so that summer thawing is completely reversed during winter, but if Arctic temperatures increase as they are expected to do [2,20], we will likely see increased thawing of currently frozen ground in the near future.

By considering conduction only, we can describe the strength of the vertical heat flux by the heat transport equation:

$$C\frac{\partial T}{\partial t}+\nabla ·\left(-\kappa \nabla T\right)=Q.$$

(1)

Here, C (Jkg⁻¹K⁻¹) is the volumetric specific heat capacity (in the following referred to as heat capacity), T ($℃$) is temperature at depth z (m), t is time, κ (Wm⁻¹K⁻¹) is thermal conductivity, and Q (Wm⁻²) is external heat flux. The effective thermal properties of a composite depend on the thermal properties of the individual constituents and their volume fractions. For permafrost that is typically a matrix, water, ice, and air (in the following defined by subscripts m, w, i, and a, respectively). We also need to include latent heat L (kJ/kg) when one or more of the materials involved changes phase (here: water/ice, L = 333.6 kJ/kg), and following [21], we do this by replacing C with an equivalent heat capacity C_eq that includes latent heat:

$$C_{eq}={\theta }_m{\rho }_m{c}_m+{\theta }_a{\rho }_a{c}_a+{\theta }_w{\rho }_w\left(c_w+\frac{\partial Ф}{\partial T}L\right)+{\theta }_i{\rho }_i\left(c_i+\frac{\partial Ф}{\partial T}L\right).$$

(2)

Density, heat capacity, and volume fraction of each constituent are defined by ${\rho }_n$, $c_n$, and ${\theta }_n$, respectively, with subscripts as defined above. Total porosity is $\theta = {\theta }_a+{\theta }_w+{\theta }_i$, and the volume fractions of water and ice are given by ${\theta }_w=(\theta -{\theta }_a)\mathsf{Ф}$ and ${\theta }_i=(\theta -{\theta }_a)-{\theta }_w$, where Ф is a variable changing from 0 to 1 for fully frozen to fully unfrozen material over the phase change temperature interval. Note that saturation here denotes total fluid content in the pore space, i.e., $\frac{{\theta }_w+{\theta }_i}{\theta }$, while ice saturation denotes the amount of the fluid that is frozen, i.e., $\frac{{\theta }_i}{{\theta }_w+{\theta }_i}$. Furthermore, we parameterize thermal conductivity using the approximation of De Vries [22]:

$${\kappa }_{eq}=\frac{f_m{\theta }_m{\kappa }_m + f_a{\theta }_a{\kappa }_a + f_w{\theta }_w{\kappa }_w + f_i{\theta }_i{\kappa }_i}{{\theta }_m{f}_m + {\theta }_a{f}_a + {\theta }_w{f}_w+ {\theta }_i{f}_i}.$$

(3)

The shape factors $f_n$ are given by [19] as:

$$f_n={\left[1+\frac{1}{3}\left(\frac{{\kappa }_n}{{\kappa }_0}-1\right)\right]}^{-1},$$

(4)

where ${\kappa }_0$ denotes the thermal conductivity of the background medium (here: ${\kappa }_0={\kappa }_m$).

By assuming that the only heat source is the surface temperature (Q = 0), we can model how temperature distribution evolves with time and depth by using Equation (1) with thermal properties as given by Equations (2)–(4) in a finite-element scheme as provided by COMSOL [23].

Porosity, saturation, and salinity all affect thawing depth in frozen ground because the effective thermal properties of a composite are sensitive to the volume fractions of each constituent. Heat transport in Arctic terrestrial sediments is particularly sensitive to changes in the amount and chemistry of pore water because it may change phase in the temperature range commonly encountered here. In this paper, we assume full saturation and, therefore, porosity is the only factor determining water content. For a material to change phase from a solid to a liquid (e.g., ice to water), energy is absorbed to change the molecular structure of the material, meaning that heat transport to depth slows down. Non-saline water changes phase instantaneously at 0 °C, while salinity leads to gradual phase change and freezing point depression. Expressions were derived by [9] relating temperature to ice saturation, and using these expressions, we can estimate the ice saturation at a given temperature, assuming known salinity. Combining this knowledge with Equation (1), we can now determine how ice saturation evolves with time and depth.

Figure 1 shows the freezing point isotherm computed using Equation (1) for a simple model consisting of quartz sand with initial temperature T0 = −5 °C throughout the sediment column for three different assumptions of surface temperature: T1 is +1 °C, T2 is +5 °C, and T3 is +10 °C. We use typical values for the physical properties of the constituents from [21,24]. All three scenarios are displayed for a high-porosity (Figure 1a) and a low-porosity sand (Figure 1b) and for high-salinity (solid line) and low-salinity pore water (dashed line). It is evident from Figure 1 that higher surface temperature, increased salinity, and lower porosity all lead to increased thawing depth. This is because higher porosity means that a higher ice content needs to thaw for the thawing depth to progress deeper, and higher salinity means that phase change from ice to water is initiated at a lower temperature.

2.2. 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 Dvorkin 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 µ. We use a rock physics model similar to the two-end member mixing approach used by Minshull, et al. [28] and Dou, 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 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 },$$

(5)

and

$$V_s=\sqrt{\mu /\rho }.$$

(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.

Figure 3 shows how temperature varies with depth at the last time step in Figure 2 for the scenarios T0, T1, and T3, for a 20 m thick model consisting of four layers. Each layer is quartz sand, but the porosity decreases with depth to mimic a compaction effect, as pressure increases with depth due to increased sediment load. A lower porosity means lesser effects of thawing on the elastic properties, because less material changes phase. For T1 and T3, the temperature distribution in the sediment column has significantly altered from the initial state, and the effects of this change are also evident in the corresponding seismic velocity profiles in Figure 3b,c. Since the only factor that varies between the models is the temperature, the large velocity change from T0 to T1 and T3 is solely ascribed to effects caused by the change in the thermal regime.

2.3. 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 \left(3 - 2\eta \right) - 16\left(1 - \eta \right)=0,$$

(7)

where $\mathsf{\zeta }=\frac{c_R^2}{V_s^2}$ and $\eta =\frac{V_s^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 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.

3. Results

To investigate the impact of thawing on seismic data, we use seismic gathers and FV-spectra from synthetic seismic data and real seismic data acquired in the Norwegian Arctic.

3.1. Synthetic Seismic Data

We use an implementation of the wavenumber integration method called OASES [40] to generate synthetic seismic data corresponding to the velocity models shown in Figure 3. We model the vertical particle displacement for a 120 m-long line with geophones spaced 1 m apart and an impulsive source.

Figure 4a–c show synthetic seismic gathers corresponding to the velocity profiles for T0, T1, and T3 of Figure 3. We observe significant differences between the various gathers, and the complexity of the data increases as the temperature variability with depth increases. For T0, only two events are visible: a non-dispersive strong event travelling at approximately 1450 m/s, which is the fundamental mode of the surface wave (A), and one weaker non-dispersive event travelling at approximately 2700 m/s, which is the direct P-wave (B). The seismic gathers for T1 and T3 are more difficult to interpret, so in addition to the seismic gathers showed here, we make use of computed group velocities and modeling of separate modes of surface waves (as done in normal mode summation [41]). For T1, dispersive waves dominate the wavefield. This includes a wide fan-shaped event with group velocities ranging from approximately 200 m/s to 1450 m/s (C), which is associated with the fundamental mode of the surface wave, and a narrower fan-shaped event with group velocities ranging from approximately 500 m/s to 850 m/s (D), which is the first higher mode. Also for T3, dispersive events dominate the wavefield. The lower fan-shaped event (E) has group velocities ranging from approximately 200 m/s to 1450 m/s, and is the fundamental mode of the surface wave. Amplitudes are particularly high at a group velocity of approximately 380 m/s, which we attribute to the fundamental mode being dispersive at low frequencies, but becoming non-dispersive and correspondingly have a constant group velocity of 380 m/s at higher frequencies. The weaker fan-shaped event (F) is the first higher mode of surface wave with group velocities ranging from 380 m/s to 1000 m/s.

Figure 4d–f show the corresponding FV-spectra, with theoretical modal curves computed from the input models overlaid on the spectra for easier mode identification. Just as with the seismic gathers, there are significant differences between the various models. For T0, the peak in the FV-spectrum appears at fairly constant phase velocities of around 1450 m/s for all frequencies, corresponding with the non-dispersive event observed in Figure 4a. T1 and T3 have similar and constant phase velocities of about 380 m/s at frequencies >120 Hz, but at frequencies between approximately 63 and 120 Hz, phase velocities are higher for T1 than T3. By comparing the FV-spectra with the modal curves, the spectrum peak for T0 corresponds well with the fundamental mode at all frequencies. For T1, it corresponds well with the fundamental mode below 31 Hz and above 63 Hz, and with the first higher mode from 31 to 63 Hz. We observe a weak fundamental mode between 31 and 63 Hz as well, but the first higher mode is dominant. For T3, the modal curves correspond well with the strong fundamental mode in most of the frequency range displayed here, while the first higher mode is dominant in a narrow frequency range of approximately 20 to 30 Hz. We also observe two weak higher modes starting at approximately 60 Hz and 110 Hz, respectively. At low frequencies for all models, we find energy at phase velocities higher than the S-wave velocity in the lower half space, which cannot be normal mode Rayleigh waves and, therefore, we interpret these as leaky wave modes.

3.2. Real Seismic Data Example

To further study the effects of thawing and freezing on the seismic wavefield, we conducted three similar seismic experiments on permafrost in Adventdalen on Svalbard in the Norwegian Arctic in October 2013, January 2014, and August 2014. Permafrost in this area typically extends to a depth of around 100 m, and the upper approximately 1.2 m is the active layer that thaws seasonally [1]. There are several boreholes in the area that indicate that the upper approximately 2–3 m of sediments in Adventdalen are eolian with a lower salinity than the deeper sediments of marine origin that extend down to the bedrock at approximately 65 m depth. However, these boreholes were drilled with the objective of obtaining information about the deeper bedrock, and thus lack detailed information about the near-surface soil properties, such as saturation, salinity, and porosity [42]. Gilbert, et al. [43] studied several cores acquired throughout Adventdalen and showed that the near-surface geology is somewhat laterally homogenous in the area, but with some local variations. Temperature measurements from a borehole near the study site show that when the experiments were conducted in October, the surface temperature had been fluctuating around 0$℃$ for a couple of weeks. In January, surface temperature had been subzero since October (average surface temperature −9.81$℃$), while in August, the surface temperature had been stable above zero since mid-May (average surface temperature 0.47$℃$). From this, we presume that ice saturation in the ground was highest in January and lowest in August. We also assume that the ice saturation was more laterally heterogeneous in October than during the other two experiments because the temperatures were close to the freezing point in the weeks prior to this experiment. The soil may also have varied laterally, but this should not affect the time-lapse effect observed in the data, since all three experiments were conducted at the same location.

The objective of the experiments was to acquire reflected, refracted, and surface waves, and we therefore used gimballed geophone strings as seismic receivers in all experiments. Each geophone string consisted of 8 geophones spread out over an interval of 2.5 m, which was also the group interval. Geophone strings reduce the effect of the ground roll (Rayleigh wave), and the actual group length is critical for acquiring surface waves. The array response dampens the ground roll at wavelengths equal to or shorter than the group length, but this effect rapidly decreases for longer wavelengths (lower frequencies). Figure 5 shows that for the seismic setup applied here, we get half of the full response already at around 80 Hz (4 m wavelength) for the air wave traveling at 330 m/s, and the quality improves rapidly with decreasing frequency (increasing wavelength). Therefore, we limit our study to frequencies <100 Hz, where the surface waves should be minorly filtered due to the geophone group length.

The seismic source was firecrackers placed in a 0.2–0.3 m deep hole in the ground. These create an impulsive source signature. The receivers were placed along a 75 m long straight line, with the first receiver 9.5 m from the shot point. Figure 6 shows the location of the study site, and a picture of the seismic spread in October 2013. As the picture and topographic map reveals, the ground below the seismic line was flat, and thus no topography effects were assumed.

Since we lack information about the exact near-surface soil properties and the corresponding elastic properties at the study site, it is not possible to compute theoretical modal curves for these data. However, since we know that salinity is assumed to increase with depth in this area, this implies that the upper sediments freeze at a higher temperature than the deeper sediments due to the freezing point depression in saline pore water. When temperatures decrease gradually with depth, like in the synthetic example, a few possible scenarios are likely: the first is that during summer, we have an irregular velocity gradient in Adventdalen: a thawed upper low velocity layer above a largely frozen low-salinity and high velocity layer, above a less frozen high-salinity and lower velocity layer. A second possibility is that temperatures are high enough during summer so that the low-salinity layer completely thaws, which implies an increasing velocity gradient with depth: a thawed low-salinity and low velocity layer above a partly frozen high-salinity and higher velocity layer. This implies a similar velocity gradient as seen for the synthetic data example. During winter, it is likely that the uppermost low velocity layer disappears, i.e., velocity decreases with depth.

Figure 7a–c show seismic gathers obtained from the three seismic experiments. The data from January are clearly of different nature to the other two gathers, while the differences between the data from October and August are more subtle. A weak event (G) can be seen in all gathers. This is most likely a refraction or reflection from the underlying stiff, consolidated sediments at approximately 65 m depth. The ringing nature of this event seen in the gathers from October and August is probably due to related multiples. Event H in the January data is the Rayleigh wave caused by the frozen, stiff near-surface sediments. The air wave is defined by event I in all gathers. The lack of distinct Rayleigh waves in the data from October and August reveals that the surface sediments are soft. As noted by [44], there could be a strong coupling between the air wave and Rayleigh waves when the upper soil layer supports surface waves with a phase velocity identical to that of the air wave. This requires a soil layer with an S-wave velocity less than 330 m/s at the surface and increasing with depth. The surface waves in the August and October data make a monochromatic wave train (J), which is further indicative of a group velocity lower than phase velocity. This is typical for Rayleigh waves when V_s is increasing with depth (i.e., phase velocity decreasing with increasing frequency).

Figure 7d–f show the corresponding FV-spectra for the three seismic experiments. The energy peak in all spectra is related to the air wave, with a velocity of about 330 m/s, but other events are also seen to vary between the different spectra. The most prominent events are observed in the data from January, where there are several strong amplitude events at increasing phase velocities with increasing frequency. This inversely dispersive trend is probably linked to several higher modes of surface waves, and they are absent in data from October and August. For the August data, some weak events seem to resemble those seen for the synthetic seismic models T1 and T3 at frequencies below 40 Hz. Thus, both the seismic gathers and the FV-spectra show changes in the freezing conditions between the experiments, specifically from January to October and August.

4. Discussion and Conclusions

Both the modeled and the real seismic data reveal that thawing has a significant impact on the seismic wavefield. Seismic properties depend on the ice saturation. Increased surface temperature causes increased heat flux, thawing, and lower seismic velocities of near-surface sediments. The largest change in seismic properties occurs at shallow depths when the surface temperature is high, because this implies the largest heat flux. Since Rayleigh waves are not inherently dispersive, surface wave dispersion means that the elastic properties alter with depth. Correspondingly, thawing or freezing will thus cause temporal changes in the FV-spectra.

In the synthetic seismic data, the fundamental mode of surface wave is dominant for all experiments. Since lower frequencies are more influenced by elastic properties of deeper sediments than higher frequencies, the fundamental mode in the FV-spectra reveals that the seismic velocities in the shallow layers decrease from the T0 to the T1 model and, furthermore, from the T1 to the T3 model. This is consistent with the input data. Since the velocity decrease extends to lower frequencies for T3 than for T1, T3 implies deeper thawing than T1 due to the reduced velocities with increased thawing.

The real seismic data are, as expected, quite different from the modeled seismic data, since the subsurface model used for the synthetic data was not designed to represent local conditions. Since the inversely dispersive trend in the FV-spectra is the strongest in January, weaker in October, and missing in August, this indicates that the presence of higher modes is related to ice saturation. Temperature data indicate that the upper sediments are most frozen in January when the higher modes are strong, and least frozen in August, when the higher modes are absent. Surface waves are generally attenuated more in unfrozen ground than in frozen ground, which also explains why more low-frequency surface waves are seen in the data from January.

We did not consider a varying salinity in the synthetic data example, meaning that seismic velocities here increase steadily with depth. At the study site of the real seismic data, the ground salinity most probably varies with depth, which implies that several velocity depth models are possible. During winter, a stiff upper layer overlies softer sediments, and in this case, strong higher modes are typically generated [14,45,46,47]. This also happens when a low velocity layer is embedded in layers of higher velocities. An analogue to a stiff frozen layer on top of a stack of softer unconsolidated sediments is pavement, which constitutes a thin, very stiff layer overlaying softer sediments. For such velocity depth models, higher order modes of Rayleigh waves have been observed as an inversely dispersive trend in FV-spectra, as well as a large amount of leaky guided waves [48], which are also seen in our synthetic data. During summer, we either have stiffer sediments embedded in softer sediments, or stiffness increasing with depth. For such velocity depth models, higher modes are rarely seen [14,45]. This supports our observation that the fundamental mode is more dominant in October and August, and higher modes are more dominant in January when the uppermost layer is frozen.

Differences in the signatures of the modeled and real seismic data can also be caused by the fact that unlike in the synthetic model, the surface temperature varied frequently between low and high temperatures in the time period between October 2013 and August 2014. Furthermore, variable precipitation throughout the year affects the data. While snow leads to effective insulation of the heat flux, rain alters the moisture content in the uppermost sediments. We also assume a simple horizontally layered geological model, and no lateral soil heterogeneity or anisotropy. Adding complexity to the geological model will certainly affect the seismic data, but will to a lesser degree affect the time-lapse effects in data collected at the same location. The seismic experiments in Adventdalen also included a second seismic data line with the same shot point, but having a 67.5 degree azimuth to the east of the current line. These data display similar FV-spectra and time-lapse effects, indicating gentle variation in the soil properties laterally. The seismic setup of the experiments was not particularly designed for recording surface waves. As full thawing only occurs in the top meter or so in Adventdalen, the receiver interval should in general be minimized to gain more high-frequency data. Although geophone strings were used, the array response retains frequencies up to 100 Hz, which is sufficient for our study.

Our results show significant time-lapse effects in both the seismic gathers and the estimated Rayleigh wave dispersion curves, which can be related to thawing of near-surface Arctic sediments. The synthetic seismic data were made with full control over the physical properties of the models, and show distinct time-lapse effects that are due to changes in the thermal regime. The time-lapse study of the real seismic data strongly indicate that low-frequency time-lapse effect are due to changes in the freezing conditions. Analysis of FV-spectra basically provide trends in the stiffness gradient (increasing/decreasing/irregular). However, to quantify exact changes in the mechanical properties of the subsurface due to thawing requires more sophisticated inversion of FV-spectra to S-wave velocity with depth. Conventional inversion techniques using normal mode methods are not designed for (irregularly) decreasing velocities with depth, or when higher-order modes dominate [45,49,50]. This was seen to occur for the real seismic data example, where higher modes became dominant with increasing frequencies. Full-waveform inversion often works for small modal jumps, but fails for large modal jumps, and where the main energy does not return to the fundamental mode at high frequencies [14,45]. An in-depth discussion or execution of more complex inversion methods is beyond the scope of this study. Our study has been limited to qualitatively assessing the differences between the various seismic gathers and FV-spectra at low frequencies. Our results also highlight the need to know the near-surface soil and pore water properties to estimate the degree of thawing from estimated velocity profiles. The contribution of this study is to demonstrate the potential of using time-lapse surface seismic data in climate monitoring of the terrestrial Arctic.