Permafrost saline water and Early to mid-Holocene permafrost aggradation in Svalbard

Deglaciation in Svalbard was followed by seawater ingression and deposition of marine (deltaic) sediments in fjord valleys, while elastic rebound resulted in fast land uplift and the exposure of these sediments to the atmosphere, whereby the formation of epigenetic permafrost was initiated. This was then followed by the accumulation of aeolian sediments, with syngenetic permafrost formation. Permafrost was studied in the eastern Adventdalen valley, Svalbard, 3–4 km from the maximum up-valley reach of post-deglaciation seawater ingression, and its ground ice was analysed for its chemistry. While ground ice in the syngenetic part is basically fresh, the epigenetic part has a frozen freshwater–saline water interface (FSI), with chloride concentrations increasing from the top of the epigenetic part (at 5.5 m depth) to about 15 % that of seawater at 11 m depth. We applied a one-dimensional freezing model to examine the rate of top-down permafrost formation, which could be accommodated by the observed frozen FSI. The model examined permafrost development under different scenarios of mean average air temperature, water freezing temperature and degree of pore-water freezing. We found that even at the relatively high air temperatures of the Early to mid-Holocene, permafrost could aggrade quite fast down to 20 to 37 m (the whole sediment fill of 25 m at this location) within 200 years. This, in turn, allowed freezing and preservation of the freshwater–saline water interface despite the relatively fast rebound rate, which apparently resulted in an increase in topographic gradients toward the sea. The permafrost aggradation rate could also be enhanced due to non-complete pore-water freezing. We conclude that freezing must have started immediately after the exposure of the marine sediment to atmospheric conditions.

1 Introduction

Cycles of global warming and cooling are well documented in geological history (e.g. Imbrie et al., 1993; Benn and Evans, 2014; Arnscheidt and Rothman, 2020). During the Pleistocene, these cycles followed Northern Hemisphere glaciation and deglaciation, which influenced both marine and land temperatures (Park et al., 2019). This also affected the extent of cryotic conditions in the periglacial environment (e.g. Murton, 2021), i.e. the distribution of permafrost, which currently covers 22 % of Northern Hemisphere land areas (Obu et al., 2019). While temperatures during the Holocene were significantly higher than during the Last Glacial Period, the retreat of glaciers and the follow-up elastic rebound and exposure of new land in the Arctic and the sub-Arctic environment allowed freezing and the aggradation of permafrost (e.g. Landvik et al., 1988). Nevertheless, the relatively high temperatures during the Early and the mid-Holocene warm period raise questions about the timing of initiation and the extent of this process (e.g. Landvik et al., 1988; Humlum, 2005).

Figure 1 (a) The study site, Adventdalen East (ADE), is located at the Adventdalen, Svalbard, on a river terrace. Map provided courtesy of the Norwegian Polar Institute. (b) Drilling at ADE, spring 2017 (top), and core samples before processing at UNIS cold room lab (bottom). Photos: Dotan Rotem and Yishai Weinstein.

In Svalbard (Fig. 1), the fast retreat of glaciers during the end of the Late Pleistocene into the beginning of the Holocene resulted in the ingression of seawater in fjord valleys, which was followed by gradual uplifting and exposure due to elastic rebound. This resulted in epigenetic permafrost aggradation followed by the deposition of fluvial and aeolian sediments and the formation of syngenetic permafrost during the last ca. 4 kyr (Gilbert et al., 2018). In the present study, we use the presence of saline water in the epigenetic permafrost to constrain the timing of freezing.

Permafrost is a soil or rock which has been below 0 ^∘C for at least 2 consecutive years (French, 2017). While winter freezing of the ground is common in land areas to a large extent, the existence of permafrost and its aggradation depend on the annual energy balance between the atmosphere and the land (Black, 1954). Accordingly, permafrost develops when the land heat loss during winter exceeds the gain during the summer for a long enough time. This is controlled by both seasonal solar radiation and the soil and rock thermal properties. Heat exchange between soil and the atmosphere is also strongly affected by land cover, whereby permafrost is usually developed neither under the sea nor beneath warm-based glaciers (Waller et al., 2012). Nevertheless, permafrost can occur beneath lagoons in association with taliks, as well as beneath bottom-fast ice in shallow water (Solomon et al., 2008). The extent and depth of permafrost can be significantly reduced by thick vegetation or snow cover (e.g. Grünberg et al., 2020). During the last glacial cycle, the Barents Sea and the Svalbard area (Fig. 1) were covered by one to three ice caps (Mangerud et al., 2002; Patton et al., 2017). Glacier retreat has been followed since by elastic rebound, which is well documented in Svalbard (Bondevik et al., 1995; Lønne and Nemec, 2004; Sessford et al., 2015), with a land rise of up to 130 m in eastern Svalbard and 65 m in the western part of the archipelago (Forman et al., 2004). In western Svalbard, the focus of this study, research indicates a fast land rise from 19 to 15 mm yr⁻¹ during the Early to the mid-Holocene (11.7–8.2 ka), which decreased from 5 to 4 mm yr⁻¹ toward the end of the mid-Holocene (Salvigsen, 1984; Sessford et al., 2015) and ca. 1 mm yr⁻¹ during the Late Holocene (last 4 kyr; e.g. Forman et al., 2004).

Land uplift and exposure are accompanied by the establishment of a surficial drainage system, as well as the development of a groundwater flow network, which strongly depends on the rate of permafrost deepening (Edmunds et al., 2001). The permeability of frozen soils is greatly reduced (Burt and Williams, 1976; Cochand et al., 2019), such that extensive permafrost prevents penetration of surface water and recharging of groundwater (McEwen and de Marsily, 1991). While in sporadic and discontinuous permafrost, groundwater flow is possible through non-frozen sections or taliks, and flow is practically impossible through continuous-permafrost land areas (Lemieux et al., 2008; Walvoord and Kurylyk, 2016), while it may be active in sub-permafrost zones or along faults in relation to pingos (Hornum et al., 2021). Flow may also provide the necessary conditions for the formation of cryopegs, which hold overcooled liquid brines (Ahonen, 2001). Although cryopegs were identified in Adventdalen (Tavakoli et al., 2021), water paths have not yet been described.

According to the Ghyben–Herzberg approximation (Bear and Dagan, 1964; Verruijt, 1968), the distance from the seawater level to the freshwater–saline water interface (FSI) should be about 1 : 40 of the groundwater head above sea level. This ratio increases with decreasing salinity of saline water. With typical Early to mid-Holocene rebound rates from 15 to 4.5 mm yr⁻¹ (Sessford et al., 2015) and assuming that the groundwater table (saturated conditions) followed the topography, the fresh–saline interface is expected to be pushed downwards from as deep as 120 m to 36 m within 200 years of exposure. A groundwater table of 1 m below the surface would result in a delay of 100 to 200 years, but also in this case a sediment section of tens of metres will be completely flushed within several hundred years unless sediment freezing practically halts flow in the sub-surface.

The objective of this paper is to test the presented hypothesis by studying the ground ice geochemistry of a permafrost core from Adventdalen, Svalbard, and by using a 1-D numerical heat transfer model to simulate permafrost aggradation under various surface temperature conditions and degrees of freezing.

2 Study site

Adventdalen is a U-shaped, glacially eroded valley located in western Spitsbergen, Svalbard, centred at 78.110^∘ N, 16.180^∘ E (Fig. 1a). During the last glacial cycle, the valley was eroded down to the basement, which was then covered by glacial deposits (Elverhøi et al., 1995; Gilbert et al., 2018). This was followed by deglaciation, which was completed ca. 10.5 ka (Mangerud et al., 1992; Svendsen and Mangerud, 1997; Lønne and Lyså, 2005; Farnsworth et al., 2020). Deglaciation was followed by up-valley seawater ingression, up to 13.5 km from the current end of the fjord (Cable et al., 2018; Lønne and Nemec, 2004), and subsequent valley infilling with sediment and a delta front prograding down-valley. Elastic rebound resulted in the exposure of the eastern part of the valley before 9.5 ka, which progressed down-valley, arriving at the current coastline location at about 4 ka (Gilbert et al., 2018). The exposed surface was first covered by fluvial sediments, followed by aeolian deposits between 4–2 ka (Gilbert et al., 2018). Permafrost in Svalbard, both epigenetic and syngenetic, is continuous and is estimated to be > 100 m thick in the valleys (Humlum, 2005). The active-layer thickness is commonly 60–100 cm in the valley bottom sediments (Christiansen and Humlum, 2003; Gilbert et al., 2018; Weinstein et al., 2019; Strand et al., 2021).

Meteorological data have been measured continuously for more than 100 years at the Svalbard Airport (e.g. Nordli et al., 2020), located on the Adventfjorden coast, ca. 12 km north-west of our study site. Mean annual air temperature (hereafter MAAT) was −5.9 ^∘C from 1971 to 2000, although in 2018 it was merely −1.8 ^∘C. Mean annual precipitation between 1971–2000 was 196 mm, while in 2018 it was 177 mm (Strand et al., 2021). MAAT and sea surface temperature during the Early to the mid-Holocene were 2–4 ^∘C higher than today, as suggested by marine molluscs' shells (Mangerud and Svendsen, 2018), lacustrine alkenones (van der Bilt et al., 2018), flora DNA (Alsos et al., 2016) and models incorporating physical and biological considerations (e.g. Park et al., 2019). Since the mid-Holocene, a continuous decline in MAAT has been recorded, which changed into a fast temperature rise during the last several decades (e.g. Christiansen et al., 2013).

The study site, Adventdalen East (ADE), is located on a river terrace (78.1722^∘ N, 16.0613^∘ E), 9.8 km up-valley from the Adventfjorden at 23 m a.s.l. (Fig. 1a and b). We drilled at the same S1 location as Gilbert et al. (2018). The permafrost section (valley-fill sediments) at the ADE site (ca. 20 m) consists of a syngenetic part from 1.0–5.5 m depth, which includes 1.5 m of shallow, fine-grained aeolian deposits and thick ice bodies of segregation ice that intruded above the fluvial gravel channel fill (Gilbert et al., 2018). Below 5.5 m the epigenetic permafrost consists of 3–4 m of fluvial sediments (mud and pebbles, ice-rich). This is underlain by back-delta, deltaic and fore-delta sediments (5.5–17.5 m), which cover glacial sand deposits (17.5–20 m) (Gilbert et al., 2018). The study site was deglaciated by 11.3 ka and emerged above seawater at 9.2 ka (Gilbert et al., 2018), exposing it to atmospheric conditions then, which allowed the development of a groundwater system on one hand and possibly the aggradation of permafrost on the other hand.

3 Methods

Two cores, 0.5 m apart, one 13 m and the other 9 m long, were retrieved at ADE in March 2017, using the UNIS permafrost drill rig (Gilbert et al., 2015), which has core barrels 43 mm in diameter (i.d.). Core length, borehole depth, core condition and gravel content were recorded in the field. Retrieved core sections were placed in plastic bags and marked with a serial number and an arrow pointing towards the core top. Cores were stored at −18 ^∘C in a freezer at UNIS until processing. Cores were sectioned in a cold room (−5 ^∘C) at 0.5 m depth intervals (Fig. 1b). Intervals of the same depth in the two cores were combined to gain enough ground ice per section for Ra isotope measurements. Samples were first scraped and then crushed to small chips, which were placed in 250 mL centrifuge tubes. Ra-free water (up to 40 mL) was added to some of the tubes to facilitate the extraction of pore fluid. Samples were then thawed in a microwave set to 600 W for 2 min, followed by centrifuging for 8 min (11 000 rpm, high g) to separate thawed water from the soil. Extracted water was run through 3 µm filters, followed by 0.45 µm filters. Most of the water was used for Ra isotope analysis (see Weinstein et al., 2019), while 30–60 mL was used for chemistry analyses. Water of the added Ra-free water was analysed to correct for element concentrations. Major elements were analysed at the Geological Survey of Israel (GSI) by inductively coupled plasma atomic emission spectroscopy (ICP-AES; Optima 3000), where Sc was added as an internal standard, whereas Cl⁻ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ was determined by potentiometer titration using a Metrohm 702 SM Titrino titrator connected to a chlorine electrode. The error for all majors is considered less than 5 %.

4 Ground ice chemistry

Major elements of thawed ground ice are presented in Table 1, and concentration profiles of Cl⁻, Na²⁺ and SO${}_{\mathrm{4}}^{-\mathrm{2}}$ are shown in Fig. 2a–c. While the salinity of ground ice in the syngenetic permafrost is that of fresh water (e.g. Cl⁻: 10–74 mg L⁻¹; Na²⁺: 10–33 mg L⁻¹; and SO${}_{\mathrm{4}}^{-\mathrm{2}}$: 9–31 mg L⁻¹), epigenetic permafrost ground ice demonstrates a trend of increasing concentrations down to 9–12 m depth: 440–3600, 80–2700 and 150–740 mg L⁻¹ of Cl⁻, Na²⁺ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, respectively. Between 9–12 m, concentrations are quite scattered, and the increasing pattern is less clear. While the Cl⁻ content (Fig. 2a) of the ground ice is no more than 15 % seawater salinity, and the salinity of a deeper-seated saline water end-member could be significantly higher, the observed increased salinity clearly presents a freshwater–saline water interface.

Table 1 Major elements (mg L⁻¹) of ground ice samples from the ADE core.

^a After Gilbert et al. (2018). ^b After de Baar et al. (2017), salinity of 35 %. ^c Freezing state as recorded in the field.

Figure 2 Major element concentrations in ground ice from Adventdalen East (ADE) deep drill hole: (a) Cl⁻, (b) Na⁺ and (c) SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (in mg L⁻¹). Panels (d)–(f) present selected equivalent ratios along the profile. The dashed line separates the syngenetic and the epigenetic permafrost (Gilbert et al., 2018).

The ionic ratio of Na²⁺ to Cl⁻ in both the epigenetic and the syngenetic permafrost mostly exceeds 1 (Fig. 2d), significantly higher than in seawater (0.86), which is probably the result of sediment dissolution (e.g. of micas), since ion exchange should result in either a conservative behaviour (during freshening, as is the case in the ADE marine section) or in Na⁺ depletion (in the case of salinization; e.g. Russak and Sivan, 2010). On the other hand, SO₄ $/$ Cl in the epigenetic permafrost is close to that of seawater (Fig. 2f), implying a relatively conservative behaviour. Nevertheless, SO₄ $/$ Cl in the syngenetic part is very variable (Fig. 2f), reaching ratios as high as 2, which could be due to shale dissolution (Hindshaw et al., 2016; Cable et al., 2018). A high concentration of SO₄ was also recorded in sub-permafrost (pingo) water, which was attributed to gypsum dissolution (Hodson et al., 2020). Ca $/$ Cl decreases with depth in the syngenetic permafrost and is very low (< 0.01) in the epigenetic permafrost (Fig. 2e), which is in agreement with freshening experiments in freshwater–saline water zones (Russak and Sivan, 2010). The high ratio of Ca $/$ Cl and SO₄ $/$ Cl at −5.45 m depth is enigmatic and may present a paleo-active-layer–permafrost table zone where major and minor elements concentrate (e.g. Cary and Mayland, 1972; Kokelj et al., 2002), but this should be further studied.

5 Model of permafrost formation

5.1 Conceptual model

To study the rate of permafrost formation, we developed a numerical model that solves the temperature distribution in space and time and freezing front progression, i.e. the Stefan solution (e.g. Šarler, 1995). Considering low horizontal temperature variations, the problem was reduced to one-dimensional, depth-dependent heat transfer with a moving internal phase-transition boundary. Various analytical and numerical methods have been developed to obtain a stable solution of the Stefan problem (e.g. Crank, 1984). However, unlike the simple and clean single-component systems, many natural systems, including water-saturated porous rocks, change their phases under a specified temperature range rather than isothermally (Lunardini, 1988; Rühaak et al., 2015). In this case, an evolving “mushy zone” emerges, which separates between the solid and liquid regions, where the thawing or freezing begins and proceeds, accompanied by latent-heat absorbance or release (e.g. Crank, 1984; Yang et al., 2020). Following this approach, instead of the mathematical boundary, we use a narrow mushy zone where phase transition occurs (shaded area in Fig. 3), with pore space consisting of a mixture of ice and water (Rubinstein et al., 1982). The local enthalpy in the mushy zone takes values in the range of those between the pure solid and liquid, and the temperature is approximated by a constant value, T=T_wft, equal to the phase-change temperature (Crank, 1984).

Figure 3 Schematic freezing profile during the top-down (epigenetic) freezing process. Initial and developing geothermal gradients are also shown; the mushy region with constant temperature is the depth zone where phase transition occurs. The dotted line schematically represents the freezing condition, where B=0 stands for fully frozen and B=1 for liquid only. Active-layer seasonality is neglected.

Heat exchange in the sub-surface is controlled by the ground temperature gradient, as well as by the soil thermal properties (i.e. thermal conductivity; Burn, 2011). Aboveground, the main factor is the air temperature, which is measured and reported as MAAT (Luo et al., 2018; Szafraniec and Dobiński, 2020) and was taken to represent the MAGST (mean annual ground surface temperature). Initial surface temperature was defined according to the temperature of the shallow seawater during the mid-Holocene (2 ^∘C; Rasmussen et al., 2012), while the initial temperature profile (solid black line in Fig. 3) was defined using the regional-average geothermal gradient of 0.033 ^∘C m⁻¹, as discussed by Olaussen et al. (2019) and Betlem et al. (2018). The lower boundary of the model was set at 300 m depth, with temperature of 12 ^∘C according to the initial temperature distribution. Throughout the simulation, we searched for the depth- and time-dependent temperature distribution T(z,t), schematically shown as a dashed line in Fig. 3.

Several MAAT values were used for the modelling: (1) the current −5.8 ^∘C (​​​​​​​measured at the Adventdalen “Polygons” site; Christiansen, 2005), 7 km from the fjord; (2) −4 ^∘C, which was taken from climate simulation models for the mid-Holocene (Park et al., 2019; see also Mangerud and Svendsen, 2018; van der Bilt et al., 2019); and (3) −3 and 0 ^∘C assumed by Humlum (2005) for the mid-Holocene. We present simulations for −4 ^∘C and with 0 ^∘C as an extreme high value. While snow may cause differences between MAAT and MAGST due to thermal insulation (Zhang, 2005), it was found that in western Svalbard, specifically in the flat landforms at Adventdalen, the differences between MAAT and MAGST are less than 0.5 ^∘C (Christiansen, 2005; Lüthi, 2010; Etzelmüller et al., 2011; Farnsworth, 2013). Therefore, MAAT was taken to represent MAGST. The amplitude of seasonal temperature oscillation at the surface was set to 12 ^∘C, similar to the current fluctuation (Nordli et al., 2014; Christiansen, 2005; Osuch and Wawrzyniak, 2017).

To distinguish between a frozen and water-saturated cell, we defined a time- and depth-dependent freezing ratio, B(z,t), shown by the dotted line in Fig. 3. B=1 means that the soil is water-saturated, while in the case of B=0 pore space is fully ice-saturated. In the mushy zone (shaded area in Fig. 3), the B value changes between 0 and 1, and the rate of its change defines the amount of energy or latent heat associated with water–ice phase transition. Since it is now well established that permafrost is not necessarily fully frozen (e.g. Keating et al., 2018; Oldenborger and LeBlanc, 2018; see Table 1), we also investigated permafrost aggradation under “partial-freezing” conditions of 25 % and 50 %. Note that our model assumes fully saturated pore-water conditions since freezing starts at sea level soon after exposure; therefore the groundwater level is expected to be at the surface.

Another factor affecting the rate of permafrost formation is the water freezing temperature (hereafter: WFT), which varies with salinity (Morgenstern and Anderson, 1973), as well as due to other environmental factors (Farouki, 1981; Morgenstern and Anderson, 1973). Since salinities at the ADE site down to 12 m do not exceed 15 % that of seawater, freshwater freezing temperature is the most appropriate for our simulations. Nevertheless, we also conducted simulations with a WFT of −2 ^∘C, which is close to that of seawater (T_m = −1.9 ^∘C; Bodnar, 1993; Marion et al., 1999), as well as with −4, −5 and −6 ^∘C, following the reports of Gilbert et al. (2019) and Tavakoli et al. (2021) about high salinities (up to 73 ppt) in cryopegs in Adventdalen.

For simplicity, thermal conductivity was calculated considering three major components, namely water, ice and dry sediments. While the thermal conductivities of water and ice are well defined (0.569 and 2.24 W m⁻¹ K⁻¹, respectively; e.g. Williams and Smith, 1989), the range of conductivity values for dry sediment in the literature is large, from 0.25 for peat to 7.6 W m⁻¹ K⁻¹ for sand, depending mainly on quartz content, water content and other factors (Farouki, 1981; Williams and Smith, 1989; Cosenza et al., 2003; Balland and Arp, 2005; Zhang and Wang, 2017; He et al., 2021). We decided to use the value of 3 W m⁻¹ K⁻¹ for the mineral fraction of the sediment (e.g. Overduin et al., 2019).

Porosity is an important factor in the aggradation or thawing of permafrost (Hornum et al., 2020). It determines both the amount of heat released or required during freezing and thawing and the thermal characteristics of the soil. Porosity of sediments at ADE was mainly set to its average value, 0.3 (Gilbert et al., 2018), although we also tested the effect of other porosity values (Appendix A). Bedrock (≥ 25 m depth) porosity was taken to be 0.1, the value suggested by Hornum et al. (2020), assuming that the bedrock is mainly composed of fractured shales (Benn and Evans, 2014; Grundvåg et al., 2019).

Freezing temperature and degree of freezing were kept uniform in each of the simulations, although as pore water freezes, the remaining fluid becomes saltier, further lowering the freezing temperature of the remaining solute (Herut et al., 1990), such that fully frozen pore space (i.e. 100 % freezing) can only be reached at extremely low temperatures, which are not relevant to our study sites, as well as to most other permafrost areas (Homshaw, 1980; Dobinski, 2011).

The 1-D freezing model provides a good approximation of freezing rate and permafrost aggradation, as shown by many studies (e.g. Harada and Yoshikawa, 1996; Kukkonen and Šafanda, 2001; Farbrot et al., 2007; Etzelmüller et al., 2011; Hornum et al., 2020). Such models actually provide the maximum rates of freezing propagation, as lateral heat transport by groundwater flow is neglected. Neglecting lateral heat transfer is quite justified, considering that (1) upstream shallow groundwater arrives from areas that were exposed earlier and therefore should not be warmer than the ADE groundwater, and (2) experimental data suggest that a temperature drop by a tenth of a degree below 0 ^∘C reduces the hydraulic conductivity by several orders of magnitude (Burt and Williams, 1976; Rühaak et al., 2015). Nevertheless, as stressed above, soon after freezing initiates, the hydraulic conductivity dramatically decreases, and the impact of lateral flow can be neglected.

5.2 Mathematical formulation

The heat conduction equation for time- and depth-dependent temperature profile, T(z,t), is (Crank, 1984)

$$\begin{array}{}\text{(1)}& \mathit{\rho }{C}_p{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\kappa }{\displaystyle \frac{\partial T}{\partial z}}\right)+Q,\end{array}$$

where Q is the energy sink or source, representing the latent heat associated with water–ice phase transition; ρ is the soil density; C_p is the specific heat capacity; and K is heat conductivity. The depth- and time-dependent material properties were calculated assuming linear superposition of the soil, water and ice properties (e.g. Lunardini, 1988). Thus, the depth-dependent density, heat capacity and thermal conductivity were calculated using the porosity, θ, and freezing ratio, B(z,t):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathit{\rho }=(\mathrm{1}-\mathit{\theta }){\mathit{\rho }}_{\mathrm{soil}}+\mathit{\theta }\left(\left(\mathrm{1}-B\right){\mathit{\rho }}_{\mathrm{ice}}+B{\mathit{\rho }}_{\mathrm{water}}\right)\\ & \mathit{\kappa }=(\mathrm{1}-\mathit{\theta }){\mathit{\kappa }}_{\mathrm{soil}}+\mathit{\theta }\left(\left(\mathrm{1}-B\right){\mathit{\kappa }}_{\mathrm{ice}}+B{\mathit{\kappa }}_{\mathrm{water}}\right)\\ & C_p=(\mathrm{1}-\mathit{\theta })C_{p_{\mathrm{soil}}}+\mathit{\theta }\left(\left(\mathrm{1}-B\right)C_{p_{\mathrm{ice}}}+B{C}_{p_{\mathrm{water}}}\right).\end{array}\end{array}$$

The thermal properties and density used for all system components (soil, ice and water) are listed in Table 2.

Table 2 Physical parameters of water, ice and sediment of the 1-D heat transfer model.

^a $D=\mathit{\kappa }/\mathit{\rho }\times C_p$. ^b Value refers to saturated pore space without air.

Out of the mushy zone, for either B = 0 or B = 1 (unfrozen sediments), the heat exchange leads to its cooling below or heating above the freezing temperature. When no latent heat is involved, assuming homogeneous heat conductivity, the heat conduction equation (Eq. 1) is reduced to

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial T}{\partial t}}=D{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}},\end{array}$$

where D (diffusivity) (m² s⁻¹) is defined as

$$D={\displaystyle \frac{\mathit{\kappa }}{\mathit{\rho }{C}_p}}.$$

In the mushy zone, where the water–ice phase transition occurs, both the B value and the thermal properties (C_p and K) are depth-dependent. Accordingly, the complete heat conduction equation (Eq. 1) is solved, including the latent-heat term. The heat source (sink) is equal to the mass of the freezing (thawing) water per unit time multiplied by the latent heat, L. The water mass is equal to the rate of the B-value change times porosity and density. Finally, the source term is

$$\begin{array}{}\text{(4)}& Q=-L\mathit{\theta }\mathit{\rho }{\displaystyle \frac{\partial B}{\partial t}}.\end{array}$$

We neglected the kinetics of the phase transition and assumed that thermodynamic equilibrium is established instantaneously in the mushy zone. This means that the rate of freezing or thawing is defined by the heat flux to and from the mushy zone with 0 < B < 1. Substituting Eq. (4) into Eq. (1) and using $\frac{\partial T}{\partial t}=\mathrm{0}$ leads to

$$\begin{array}{}\text{(5)}& \mathit{\theta }L\mathit{\rho }{\displaystyle \frac{\partial B}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\kappa }{\displaystyle \frac{\partial T}{\partial z}}\right).\end{array}$$

The above equations are solved numerically for two functions, T(z,t) and B(z,t), using the explicit-in-time finite-difference scheme. These functions were approximated using the constant grid steps in depth Δz and in time Δt:

$$\begin{array}{c}{\displaystyle }{T}_{n,m}=T(n\mathrm{\Delta }z,m\mathrm{\Delta }t\\ {\displaystyle }{B}_{n,m}=B(n\mathrm{\Delta }z,m\mathrm{\Delta }t),\end{array}$$

where n is the grid point number (z=nΔz), and m is the time step number (t=mΔt). With this notation, the finite-difference form of the heat conduction equation (Eq. 1) is

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{\rho }}_n{C}_{p_n}{\displaystyle \frac{T_{n,m+\mathrm{1}}-T_{n,m}}{\mathrm{\Delta }t}}=\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{z}^{\mathrm{2}}}}\left[{\displaystyle \frac{{\mathit{\kappa }}_{m+\mathrm{1}}+{\mathit{\kappa }}_m}{\mathrm{2}}}\left(T_{n+\mathrm{1},m}-T_{n,m}\right)\right.\\ & \left.-{\displaystyle \frac{{\mathit{\kappa }}_m+{\mathit{\kappa }}_{m-\mathrm{1}}}{\mathrm{2}}}\left(T_{n,m}-T_{n,m-\mathrm{1}}\right)\right]+L\mathit{\rho }\mathit{\theta }{\displaystyle \frac{B_{n,m+\mathrm{1}}-B_{n,m}}{\mathrm{\Delta }t}},\end{array}\end{array}$$

where density and the thermal properties are calculated using Eq. (2).

Equation (6) is solved using a time step of Δt = 10 800 s (i.e. 3 h) and a depth spacing of Δz=0.25 m. The solution was obtained for the model size down to 300 m depth, summing up to 1200 grid points. These numerical parameters satisfy the von Neumann stability condition for the explicit-in-time numerical scheme (e.g. Ames, 1977) for the material properties of Table 2. The numerical code was written with Python (McKinney, 2012). It allows simulation of the permafrost dynamics and sub-surface sediments freezing under various scenarios of MAAT, water freezing temperatures (WFTs) and freezing extent of pore-space water.

5.3 Model results

We present the results of model runs with various combinations of surface temperature, water freezing temperature, and sediments porosity and thermal conductivity. Complementary modelling results are presented in the Appendix A (Figs. A1–A4). In all cases, simulations started in the spring (May) and maintained an amplitude of 12 ^∘C around the chosen (fixed) MAAT.

Figure 4 One-dimensional freezing model results for (a) 1 year. The model starts (and finishes) in spring (defined as halfway between minimum and maximum surface T; e.g. May), followed by a summertime increase in temperature (solid and dashed lines), cooling and start of freezing during fall (dotted line; e.g. November), and colder winter months (two-point dashed line and wide dashed line). It then concludes in the spring (one-point dashed line). (b) Results for 6 years (starting in spring). (c) Results for 50 to 1000 years. MAAT of −4 ± 12 ^∘C, WFT set at 0 ^∘C, 100 % freezing. Full black circles denote −4 ^∘C. Note the different scale between (a), (b) and (c).

We first present simulations with a MAAT of −4 ± 12 ^∘C, WFT of 0 ^∘C and complete (100 %) freezing. Figure 4a presents results of a 1-year simulation, with temperature profiles shown every second month. The model suggests that freezing under these conditions can reach down to 3 m within the first year. The freezing depth increases to 6 m within 6 years, with a slight but significant deepening of the inflection point (Fig. 4b). After 50 years, freezing arrives at 16–17 m, and within 1000 years the freezing front is already at 73 m (Fig. 4c), within basement rocks (considering that sediment cover at ADE is ca. 25 m). The depth affected by cooling also progresses with time (< 50 m in 50 years, > 150 m in 1000 years), and the T profile approaches linearity.

Lowering the WFT results in a decrease in the freezing rate. For example, with a WFT of 0 ^∘C the freezing front will reach 23 and 37 m after 100 and 200 years, while with a WFT of −2 and −4 ^∘C it will reach 15 and 20 m and 5.5 and 6.5 m during the same periods, respectively (Fig. 5). A thin layer of permafrost may aggrade even under a WFT of −5 ^∘C, which is lower than the MAAT of −4 ^∘C. The freezing front will reach 2.5 m, resulting in 0.5 m of permafrost, considering a hypothetical, conservative, 2 m thick active layer for the high temperatures of the mid-Holocene, which is deeper than present (Humlum et al., 2003; Weinstein et al., 2019). This is because the thermal conductivity of ice is higher than that of water (Farouki, 1981), which results in a deeper advance of the winter freezing front (through ice) than the advance of summer thawing (through water). Last, with a WFT of −6 ^∘C (i.e. significantly lower than the MAAT), there is no apparent permafrost aggradation (not shown).

Figure 5 Results of 1-D freezing simulations with mid-Holocene MAAT of −4 ± 12 ^∘C (Park et al., 2019). Pore-water freezing temperature (WFT) taken to be 0, −2 and −4 ^∘C. Simulations were run for 100 and 200 years.

In Fig. 6, we examine the effect of partial freezing. Partial freezing (25 % and 50 % in our scenarios) results in deepening of the freezing depth, but differences are relatively small for a WFT lower than 0 ^∘C. With a MAAT of −4 ^∘C and WFT of 0 ^∘C, after 200 years freezing depth will reach 37, 42 and 45 m with 100 %, 50 % and 25 % freezing, respectively. Lowering the WFT to −2 ^∘C, freezing depth is indifferent to the freezing degree (Fig. 6). This is because of the trade-off between reducing latent heat and the lower thermal conductivity of the partially frozen pore space. With a WFT of −4 ^∘C, under 100 % freezing 6.5 m of permafrost will aggrade after 200 years, while with partial-freezing scenarios it will aggrade to 3.25 and 4.25 m for 50 % and 25 %, respectively (Fig. 6). With a WFT of −5 ^∘C a thin layer of permafrost may aggrade under the 100 % freezing scenario (0.5 m, not shown), while with partial freezing of 25 %–50 %, it will not aggrade (i.e. no permafrost will develop, assuming an active-layer depth of 1–2 m), which further suggests that when WFT is lower than the average MAAT (−4 ^∘C), aggradation is controlled by the ice thermal conductivity rather than by latent heat.

Figure 6 Mid-Holocene permafrost aggradation with variable freezing degrees. MAAT was set to −4 ^∘C (Park et al., 2019), and freezing proportions were taken to be 100 %, 50 % and 25 %. The shaded greyish zone represents a hypothetical, conservative active layer, which is taken to be 2 m thick. We note that the active layer in Adventdalen is usually ≤ 1 m, and 2 m was chosen due to the higher temperature during the mid-Holocene.

Figure 7 Early and mid-Holocene freezing depth in the first 50 and 200 years with a MAAT of −4 and 0 ^∘C and WFT of 0 and −2 ^∘C: (a) 100 % freezing, (b) 25 % freezing. Also shown are curves of two scenarios of corresponding depths of the freshwater–saline water interface, using a 1 : 40 Ghyben–Herzberg (G–H) approximation for isostasy of 4.5 and 15 mm yr⁻¹. See text for more details. Panels (c) and (d) are close-ups of (a) and (b) for the first 50 years. The fast deepening in (a) and (b) is due to the change in porosity as the freezing front reaches the bedrock (25 m). The legend is the same for all panels.

In general, higher MAAT results in a lower aggradation rate; e.g. with a MAAT of 0 ^∘C, a WFT of 0 ^∘C and 100 % freezing, freezing depth will reach 9.5 m within 200 years. For the same setting, but with a MAAT of −4 ^∘C, the freezing front will result in significant deepening and will reach 37 m (Fig. 7a). Setting the system to 25 % freezing, a MAAT of −4 ^∘C, and a WFT of 0 and −2 ^∘C, freezing depth will reach 47 and 20 m, respectively (Fig. 7b). By using a MAAT of 0 ^∘C, freezing depth will reach 4.25 m for a WFT of 0 ^∘C and 2.75 m for a WFT of −2 ^∘C; in both scenarios freezing depth is constant for 200 years (Fig. 7b). Moreover, we show that even with MAAT > WFT (e.g. 0 and −2 ^∘C, respectively) freezing will arrive at 2.75 m after 10 years and may create a thin layer of permafrost for at least 200 years. As mentioned above, this is due to the asymmetry in the seasonal freezing–thawing process. Higher conductivity during winter (frozen pore space) enhances the loss of heat, while the lower conductivity during summer (warming front goes through thawed pore space) slows the thawing process (Kukkonen and Šafanda, 2001). We note that lower proportions of freezing (e.g. 25 %) will have the effect of reducing this asymmetry due to the higher proportions of liquid water in the cryotic pore space and therefore lower thermal conductivity during freezing. Accordingly, permafrost deepening is hardly observed in the scenario of MAAT = 0 ^∘C (in particular with WFT = −2 ^∘C).

6 Discussion

6.1 Ground ice salinity and the frozen interface

When the freezing front propagates downwards in recently emerged land, and epigenetic permafrost is formed, it might freeze old sub-surface brines (Cascoyne, 2000). As the freezing process proceeds, solute concentrations in the non-frozen residual water commonly increase (e.g. Cocks and Brower, 1974; Herut et al., 1990; El Kadi and Janajreh, 2017). This results in a pore space with high salt concentrations. These brines may then migrate away from the freezing surface, driven by density and capillary forces, and coalesce to form separate saline water lenses (“cryopegs”; Cascoyne, 2000). The level of salinity and water composition will depend on the initial water composition and the extent of freezing.

Complete permafrost freezing can hardly be obtained, since the eutectic point of seawater freezing is at −36 to −54 ^∘C (Gitterman, 1937; Ringer, 1928; Nelson and Thompson, 1954; Marion et al., 1999), while in Adventdalen permafrost temperatures do not usually drop below −6 ^∘C below the depth of zero annual amplitude (ZAA) (Christiansen et al., 2010) and are never lower than −12 ^∘C, even in the shallow permafrost (Christiansen et al., 2020; Isaksen et al., 2007). Although the eutectic point is well below the expected temperature values, the freezing–salt expulsion process still prevails. Under these conditions, the permafrost pore space should hold a small fraction of residual brine solution, which contains most of the solutes originally dissolved in the bulk pore-space water. Partly frozen permafrost has been often observed in the study area during drilling, in particular deeper than a few metres. This was also found in both geophysical and geochemical observations (e.g. Keating et al., 2018; Weinstein et al., 2019). Nevertheless, when ground ice is thawed, the extracted fluid from the relatively large segments used in this study should roughly indicate the salinity of the original in situ pore fluid, assuming no major brine migration had occurred. Pore-water composition may be significantly altered from the original fluid that circulated in the sediments (e.g. seawater) due to ion exchange or dissolution prior to when or even after cryotic conditions occurred, which is reflected in the Na $/$ Cl and SO₄ $/$ Cl ratios in the thawed ground ice (Fig. 2). Nevertheless, Cl⁻ concentration is probably close to and represents the salinity of the original pore fluid. The relatively low salinity (Fig. 2) and the evident mixing profile (Fig. 2a, b and c) suggest that this is not the case at the ADE site and that the observed freshwater–seawater interface is an in situ observation. We refer to the chemistry of the extracted fluid as “ground ice chemistry”, although it could also be that some of it was not actually frozen.

Cable et al. (2018) presented ground ice chemistry of cores from the Adventdalen, albeit closer to the current fjord (< 4 km), west of the ADE site. In these cores, chloride, sodium and sulfate concentrations at depths of 3–11 m were up to 50 % that of seawater. At ADE, farther away from the sea, ground ice in the epigenetic permafrost, 5.5 m from the surface and deeper, shows a gradual increase in salinity (i.e. fresh–saline interface), with Cl⁻ concentrations reaching 15 % that of seawater at 9 m below the surface. Although salinities do not change much between 9–12 m, it is likely that more saline water, close to seawater salinity, either exists today or existed in the past (prior to freezing) at deeper permafrost levels. The existence of a fresh–saline interface in the very shallow permafrost, in the top of the epigenetic permafrost, suggests that freezing at ADE occurred immediately after emergence above seawater. This is further discussed below.

6.2 Rebound, exposure and fresh–saline interface deepening

Assuming that Early Holocene (11–8 ka) precipitation was slightly higher than present (200 mm yr⁻¹; McFarlin et al., 2018; Kjellman et al., 2020) and using a conservative infiltration factor of 0.2 (whether thawing snow or direct rain) and the porosity used in our simulations (0.3), this amounts to an effective annual infiltration of ca. 130 mm yr⁻¹. This could easily keep up with the Early Holocene rebound rates of 15 mm yr⁻¹ (established for the nearby Sassendalen valley; Salvigsen, 1984; Sessford et al., 2015), therefore preserving the groundwater table close to the surface of the emerging land. Using a Ghyben–Herzberg approximation (Bear and Dagan 1964), this would result in an Early Holocene fresh–saline interface deepening of ca. 60 m in 100 years (Fig. 7), assuming the saline waterbody had a common seawater density of 1025 kg m⁻³. FSI deepening could even be faster if the deep waterbody is less saline than seawater. Even if sub-aerial exposure occurred later, during the mid-Holocene, when the rebound rate decreased to 4.5 mm yr⁻¹ (Forman et al., 2004), the fresh–saline interface would still deepen at a rate of 180 mm yr⁻¹, i.e. 18 m in 100 years (Fig. 7).

The existence of a mixing zone at the top of the epigenetic permafrost (from 5.5 m below the current terrain surface), with Cl⁻ content 15 % that of seawater 3.5 m below the Early Holocene surface, suggests that the marine sediment section at ADE was hardly flushed with meteoric water. This further suggests that permafrost aggradation commenced shortly after emergence above the sea (e.g. Kasprzak et al., 2020). Indeed, some of the simulated freezing scenarios can clearly cope with the above fresh–saline interface deepening rates (e.g. MAAT of −4 ^∘C and WFT of 0 ^∘C; Fig. 7). Moreover, assuming that partial freezing (e.g. 25 % or 50 %) can also block flushing, this results in even faster permafrost aggradation (Fig. 7b). However, as permafrost deepens, the freezing rate slows down (e.g. Fig. 7), and none of the scenarios can cope with the assumed deepening of the fresh–saline interface, which should result in flushing of deeper zones.

It seems likely that the key factor in fresh–saline interface fossilization in a continuous-permafrost landscape is the permanent freezing of the very shallow permafrost, which hydraulically disconnects the sub-permafrost zone from the surface and prevents recharge of this zone with meteoric water. As shown (Figs. 6 and 7), freezing of the top 3–5 m can occur within several years even with the relatively high temperatures of the Early to mid-Holocene; e.g. a MAAT of −4 ^∘C or even warmer and a WFT of 0 or −2 ^∘C will result in 5 m of frozen sediments within 5 to 35 years (Fig. 7). Therefore, the fresh–saline interface could effectively be preserved.

We did not include salt diffusion in our model (e.g. Angelopoulos et al., 2019), a process that will reduce the WFT as freezing progresses. This can explain the reason for partially frozen samples extracted from the epigenetic section (Table 1). Including salt diffusion, we assume the freezing front may have advanced somewhat more slowly than suggested by our model.

6.3 Permafrost aggradation during the Holocene

Gilbert et al. (2018) suggested that our drilling site at ADE emerged from the sea at 10 to 9 ka and that the delta front advanced westwards at a rate of 4.4 m yr⁻¹ prior to 9.2 ka, which decreased to 0.9 m yr⁻¹ during the rest of the Holocene. Considering the relatively high rebound rates between 9 and 8 ka (e.g. 15–19 mm yr⁻¹), this suggests that the land surface at ADE reached 3–4 m above sea level and a topographic gradient of 1 %–2 % towards the sea within 200 years. Assuming the groundwater table was close to the surface, this should further result in a good flushing of the sub-surface, unless freezing took control (Fig. 8). The observed mixing zone, which reaches the very top of the pre-Late Holocene surface, suggests that freezing started within just a few years after exposure.

Figure 8 Conceptual presentation of the type of groundwater freezing processes in Adventdalen since the last glacial cycle with a focus on the freshwater–saline water interface (FSI) at the ADE site.

Our simulations suggest that both cryotic conditions (i.e. < 0 ^∘C) and actual ground ice formation started very soon after exposure to the atmosphere (Figs. 5, 6 and 7) and that significant freezing depths of 20–37 m could be achieved in 200 years (Figs. 6 and 7). This is true for both the warmer Early to mid-Holocene period (Kutzbach and Guetter, 1986; McFarlin et al., 2018; Mangerud and Svendsen, 2018; Park et al., 2019; Kjellman et al., 2020) and for any sub-zero MAAT, regardless of the size of the annual fluctuations. We even tested a MAAT of +1 ^∘C and found that some freezing could occur (not shown), which is a seasonal effect derived from the different thermal conductivities of ice and water. We also tested our model with a lower thermal conductivity by an order of magnitude, following recently published research in Adventdalen (Hornum et al., 2020), and found that permafrost aggradation, although shallower, will result in freezing of the FSI under mid-Holocene conditions. Our simulations are in good agreement with Harada and Yoshikawa (1996), who used a 1-D model with a MAGST of −5.7 ^∘C but not completely saturated sediments and found that 533 years is needed to freeze 31.7 m of sediment in Moskuslagoon, slightly to the west of our site, on the Adventfjorden shore.

The ADE site was free of seawater in the Early Holocene prior to 9.2 ka (Gilbert et al., 2018). At that time, an abrupt cooling has been described in Svalbard (Mangerud and Svendsen, 2018; van der Bilt et al., 2018, 2019). The presented model results show that the initiation of permafrost and its gradual aggradation are possible under relatively high temperatures (though a MAAT of ≤ 0 ^∘C) of the mid-Holocene. Christiansen et al. (2013) pointed out that local topographic conditions and winds in Adventdalen can induce lower temperatures in low-altitude depressions, which could enhance the permafrost aggradation during the mid-Holocene.

Our simulations also show that a lower percentage of freezing (e.g. 25 %) may enhance the permafrost aggradation rate (e.g. Fig. 7b); however, this is not true for a lower WFT (e.g. −2 ^∘C) or for a relatively high MAAT (e.g. 0 ^∘C), which as mentioned above is due to the trade-off between latent-heat and thermal-conductivity differences between ice and liquid water. Nevertheless, a MAAT of 0 ^∘C seems unlikely (e.g. van der Bilt et al., 2019).

In summary, the simulations suggest that permafrost aggradation could and did occur immediately following exposure also during the relatively warm period of the Early to mid-Holocene (10–8 ka). This is in disagreement with Hornum et al. (2020), who suggested that while freezing did occur during the Early Holocene, the ground thawed during the mid-Holocene and refroze at about 6.5 ka.

Late Pleistocene to mid-Holocene permafrost dynamics were studied in other permafrost regions. While in some areas there are already records of permafrost degradation and peatland expansion in the Early Holocene (post-deglaciation; e.g. Lenz et al., 2015; Kaufman et al., 2015; Grinter et al., 2018), cumulated evidence indicates that air temperature during this period was highly variable, sometimes higher and sometimes lower than at present (Kaufman et al., 2015). Nevertheless, it is a common observation that during the Holocene thermal maximum (mid-Holocene, 8.2–4.2 ka) permafrost degraded, and thermokarst peaked (e.g. Lenz et al., 2015; Ulrich et al., 2017; Anderson et al., 2019). Permafrost aggradation resumed post-6 ka and mainly during the past 3–4 kyr (e.g. Grinter et al., 2018; Treat and Jones, 2018). As shown above, we believe that this was not the case in Svalbard, and exposed lands continued to freeze throughout the Holocene.

7 Summary and conclusions

Arctic landscapes, including western Svalbard, were rising relatively fast in the Early to mid-Holocene due to glacial isostatic rebound. Accordingly, the preservation of frozen saline water (mixing zone) at a very shallow depth of 5.5 to 12 m is taken to be evidence for fast permafrost aggradation, which could halt the infiltration of fresh meteoric water and the flushing of saline water toward the sea. This is despite the prevailing relatively high temperatures during this period.

Our modelling confirms that freezing could progress relatively fast down the exposed Adventdalen sediments, i.e. to 20–35 m within 200 years, even under the reconstructed relatively high mid-Holocene air temperature of −4 ^∘C used in the 1-D model.

The modelling further suggests that permafrost may aggrade even when water freezing temperature (WFT) is slightly lower than mean annual air temperature (MAAT). This is attributed to the differences in thermal properties between ice and liquid water, where the higher conductivity of ice allows a faster propagation of the freezing front during winter than the propagation of thawing during the summer.

Non-complete freezing of the cryogenic pore space (which is often the case) could result in faster deepening of the freezing front when MAAT is smaller than WFT or even when it is higher, but in the latter case the difference is not large (e.g. MAAT < 0 ^∘C and 0 < WFT > −2 ^∘C). However, when MAAT ≫ WFT (e.g. MAAT = 0 ^∘C and WFT ≤ −2 ^∘C), the presence of liquid water in the pore space and its lower thermal conductivity would result in a halt of permafrost aggradation.

This concept of fast freezing under relatively high air temperatures may suggest that recently exposed areas may still go through permafrost aggradation even under the current global warming. Also, it could imply that a short (years to decades) cooling period could significantly slow down permafrost thawing.

Appendix A: Porosity analyses

Selected simulation results demonstrating the effects of the porosity values on the rate of permafrost formation are presented in Figs. A1–A4. In general, higher porosity (i.e. more pore water to freeze) will result in slower permafrost aggradation due to the higher latent heat involved.

Figure A1 Simulations of freezing front progress with different porosities for a MAAT of −4 ^∘C, a WFT of 0 ^∘C and 100 % freezing. The inset presents results for 1000 years. The fast deepening at a depth of > 25 m is due to the change in porosity as the freezing front reaches the bedrock.

Figure A2 Simulations as in Fig. A1, but with a WFT of −2 ^∘C.

Figure A3 Simulations as in Fig. A1 (MAAT = −4 ^∘C, WFT = 0 ^∘C), but with 25 % freezing.

Figure A4 Simulations as in Fig. A1 (MAAT = −4 ^∘C, WFT = −2 ^∘C), but with 25 % freezing.